Changeset 757 for trunk/LMDZ.MARS/libf/phymars/nlte_aux.F

Timestamp:

Aug 7, 2012, 3:14:07 PM (12 years ago)

Author:

emillour

Message:

Mars GCM:

• Improvement of the NLTE 15um scheme (for running with nltemodel = 2); now MUCH faster than previously (by a factor 5 or so):
• Improvements included to the parameterization:
  • Cool-to-space calculation included above P(atm)=1e-10, with a soft merging to the full result (without the CTS approximation) below that level
  • exhaustive cleaning of the code, including FTNCHK and reordering of loops, subroutines and internal calls
  • simplification of the precomputed tables of CO2 bands' atmospheric transmittances
  • the two internal grids (the one used in the CTS region and the one below) have been , in order to reduce the CPU time consumption
  • reading of the spectroscopic histograms is made only once, at the beginning of the GCM, to avoid repetitive readings of ASCII files
  • F90 matrix operations (matmul,...) included.
• Changes in routines:
  • removed nlte_leedat.F
  • updated nlte_calc.F, nlte_tcool.F, nlte_aux.F
  • updated nlte_commons.h, nlte_paramdef.h
  • added nlte_setup.F
• Important: The input files (in the NLTEDAT directory) read as input by these routines have changed. the NLTEDAT directory should now on contain only the following files:

deltanu26.dat enelow27.dat hid26-3.dat parametp_Tstar_IAA1204.dat
deltanu27.dat enelow28.dat hid26-4.dat parametp_VC_IAA1204.dat
deltanu28.dat enelow36.dat hid27-1.dat
deltanu36.dat hid26-1.dat hid28-1.dat
enelow26.dat hid26-2.dat hid36-1.dat

FGG+MALV

File:

• trunk/LMDZ.MARS/libf/phymars/nlte_aux.F (modified) (16 diffs)

trunk/LMDZ.MARS/libf/phymars/nlte_aux.F

@@ +1 @@
+c**********************************************************************
+
+c     Includes the following old 1-D model files/subroutines
+
+c     -MZTCRSUB_dlvr11.f
+c     *dinterconnection
+c     *planckd
+c     *leetvt
+c     -MZTFSUB_dlvr11_02.f
+c     *initial
+c     *intershphunt
+c     *interstrhunt
+c     *intzhunt
+c     *intzhunt_cts
+c     *rhist
+c     *we_clean
+c     *mztf_correccion
+c     *mzescape_normaliz
+c     *mzescape_normaliz_02
+c     -interdpESCTVCISO_dlvr11.f
+c     -hunt_cts.f
+c     -huntdp.f
+c     -hunt.f
+c     -interdp_limits.f
+c     -interhunt2veces.f
+c     -interhunt5veces.f
+c     -interhuntdp3veces.f
+c     -interhuntdp4veces.f
+c     -interhuntdp.f
+c     -interhunt.f
+c     -interhuntlimits2veces.f
+c     -interhuntlimits5veces.f
+c     -interhuntlimits.f
+c     -lubksb_dp.f
+c     -ludcmp_dp.f
+c     -LUdec.f
+c     -mat_oper.f
+c     *unit
+c     *diago
+c     *invdiag
+c     *samem
+c     *mulmv
+c     *trucodiag
+c     *trucommvv
+c     *sypvmv
+c     *mulmm
+c     *resmm
+c     *sumvv
+c     *sypvvv
+c     *zerom
+c     *zero4m
+c     *zero3m
+c     *zero2m
+c     *zerov
+c     *zero4v
+c     *zero3v
+c     *zero2v
+c     -suaviza.f
+
+c**********************************************************************
+
+
+c     *** Old MZTCRSUB_dlvr11.f ***
+
+!************************************************************************
+
+!      subroutine dinterconnection ( v, vt )
+
+
+************************************************************************
+
+!      implicit none
+!      include 'nlte_paramdef.h'
+
+c     argumentos
+!      real*8 vt(nl), v(nl)
+
+c     local variables
+!      integer  i
+
+c     *************
+!
+!      do i=1,nl
+!         v(i) = vt(i)
+!      end do
+
+!      return
+!      end
+
 c***********************************************************************
-c     File with all subroutines required by mztf      
-c     Subroutines previously included in mztfsub_overlap.F
-c                                               
-c     jan 98    malv            basado en mztfsub_solar        
-c     jul 2011 malv+fgg   adapted to LMD-MGCM
-c                                               
-c contiene:                                     
-c     initial                                 
-c     intershape                              
-c     interstrength                           
-c     intz                                    
-c     rhist                                   
-c     we                                      
-c     simrul                                 
-c     fi                                      
-c     f                                       
-c     findw                                   
-c     voigtf                                  
+      function planckdp(tp,xnu)
 c***********************************************************************
-                                                
-c     ****************************************************************
-      subroutine initial                             
-                                                
-c     ma & crs  !evita troubles 16-july-96           
-c     ****************************************************************
-                                                
-      implicit none                                  
-                                                
+
+      implicit none
+
+      include 'nlte_paramdef.h'
+
+      real*8 planckdp
+      real*8 xnu
+      real tp
+
+      planckdp = gamma*xnu**3.0d0 / exp( ee*xnu/dble(tp) )
+                                !erg cm-2.sr-1/cm-1.
+
+c     end
+      return
+      end
+
+c***********************************************************************
+      subroutine leetvt
+
+c***********************************************************************
+
+      implicit none
+
       include 'nlte_paramdef.h'
       include 'nlte_commons.h'
-                                                
-c local variables                               
-      integer   i                                     
-                                                
-c     ***************                               
-                                                
-      eqw = 0.0d00                                   
-      aa = 0.0d00                                    
-      bb = 0.0d00                                    
-      cc = 0.0d00                                    
-      dd = 0.0d00                                    
-                                                
-      do i=1,nbox                                    
-         ua(i) = 0.0d0                                
-         ccbox(i) = 0.0d0                             
-         ddbox(i) = 0.0d0                             
-      end do                                         
-                                                
-      return                                         
-      end                                            
-                                                
-c     **********************************************************************
-      subroutine intershape(alsx,alnx,adx,xtemp)     
-c     interpolates the line shape parameters at a temperature xtemp from    
-c     input histogram data.                         
-c     **********************************************************************
-                                                
-      implicit none                                  
-      
+
+c     local variables
+      integer i
+      real*8    zld(nl), zyd(nzy)
+      real*8  xvt11(nzy), xvt21(nzy), xvt31(nzy), xvt41(nzy) 
+
+c***********************************************************************
+
+      do i=1,nzy
+         zyd(i) = dble(zy(i))
+         xvt11(i)= dble( ty(i) )
+         xvt21(i)= dble( ty(i) )
+         xvt31(i)= dble( ty(i) )
+         xvt41(i)= dble( ty(i) )
+      end do
+
+      do i=1,nl
+         zld(i) = dble( zl(i) )
+      enddo
+      call interhuntdp4veces ( v626t1,v628t1,v636t1,v627t1, zld,nl,
+     $     xvt11, xvt21, xvt31, xvt41, zyd,nzy, 1 )
+
+
+c     end
+      return
+      end
+
+
+c     *** MZTFSUB_dlvr11_02.f ***
+
+
+c     ****************************************************************
+      subroutine initial
+
+c     ****************************************************************
+
+      implicit none
+
       include 'nlte_paramdef.h'
       include 'nlte_commons.h'
-                                                
-c arguments                                     
-      real*8 alsx(nbox_max),alnx(nbox_max),adx(nbox_max),
-     &     xtemp(nbox_max)     
-                                                
-c local variables                               
-      integer   i, k                                  
-                                                
-c     ***********                                   
-                                                
-!     write (*,*)  'intershape  xtemp =', xtemp                      
-                                                
-      do 1, k=1,nbox     
-         if (xtemp(k).gt.tmax) then 
-            write (*,*) ' WARNING !  Tpath,tmax= ',xtemp(k),tmax
-            xtemp(k) = tmax        
-         endif
-         if (xtemp(k).lt.tmin) then
-            write (*,*) ' WARNING !  Tpath,tmin= ',xtemp(k),tmin
-            xtemp(k) = tmin        
-         endif   
-               
-         i = 1                                        
-         do while (i.le.mm)                           
-            i = i + 1                                  
-            
-            if (abs(xtemp(k)-thist(i)) .lt. 1.0d-4) then !evita troubles     
-               alsx(k)=xls1(i,k) !16-july-1996      
-               alnx(k)=xln1(i,k)                        
-               adx(k)=xld1(i,k)                         
-               goto 1                                   
-            elseif ( thist(i) .le. xtemp(k) ) then     
-               alsx(k) = (( xls1(i,k)*(thist(i-1)-xtemp(k)) +       
-     @              xls1(i-1,k)*(xtemp(k)-thist(i)) )) /
-     $              (thist(i-1)-thist(i))
-               alnx(k) = (( xln1(i,k)*(thist(i-1)-xtemp(k)) +       
-     @              xln1(i-1,k)*(xtemp(k)-thist(i)) )) /
-     $              (thist(i-1)-thist(i))
-               adx(k)  = (( xld1(i,k)*(thist(i-1)-xtemp(k)) +       
-     @              xld1(i-1,k)*(xtemp(k)-thist(i)) )) /
-     $              (thist(i-1)-thist(i))
-               goto 1                                   
-            end if                                     
-         end do                                       
-         write (*,*)  
-     @        ' error in xtemp(k). it should be between tmin and tmax'
- 1    continue                                      
-                                                
-      return                                        
-      end                                            
+
+c     local variables
+      integer   i
+
+c     ***************
+
+      eqw = 0.0d00
+      aa = 0.0d00
+      cc = 0.0d00
+      dd = 0.0d00
+
+      do i=1,nbox
+         ccbox(i) = 0.0d0
+         ddbox(i) = 0.0d0
+      end do
+
+      return
+      end
+
 c     **********************************************************************
-      subroutine interstrength (stx, ts, sx, xtemp)  
-c     interpolates the line strength at a temperature xtemp from            
-c     input histogram data.                         
+
+      subroutine intershphunt (i, alsx,adx,xtemp)
+
 c     **********************************************************************
-                                                
-      implicit none                                  
-                                                
+
+      implicit none
+
       include 'nlte_paramdef.h'
       include 'nlte_commons.h'
-                                                
-c arguments                                     
-      real*8            stx     ! output, total band strength    
-      real*8            ts      ! input, temp for stx              
-      real*8            sx(nbox_max) ! output, strength for each box  
-      real*8            xtemp(nbox_max) ! input, temp for sx        
-                                                
-c local variables                               
-      integer   i, k                                  
-                                                
-c       ***********                                   
-                                                
-      do 1, k=1,nbox                                 
-!          if(xtemp(k).lt.ts)then
-!             write(*,*)'***********************'
-!             write(*,*)'mztfsub_overlap/EEEEEEH!',xtemp(k),ts,k
-!             write(*,*)'***********************'
-!          endif
-         if (xtemp(k).gt.tmax) xtemp(k) = tmax        
-         if (xtemp(k).lt.tmin) xtemp(k) = tmin        
-         i = 1                                        
-         do while (i.le.mm-1)                         
-            i = i + 1         
-!            write(*,*)'mztfsub_overlap/136',i,xtemp(k),thist(i)
-            if ( abs(xtemp(k)-thist(i)) .lt. 1.0d-4 ) then         
-               sx(k) = sk1(i,k)                         
-!              write(*,*)'mztfsub_overlap/139',sx(k),k,i
-               goto 1                                   
-            elseif ( thist(i) .le. xtemp(k) ) then     
-               sx(k) = ( sk1(i,k)*(thist(i-1)-xtemp(k)) + sk1(i-1,k)*           
-     @              (xtemp(k)-thist(i)) ) / (thist(i-1)-thist(i))  
-!              write(*,*)'mztfsub_overlap/144',sx(k),k,i
-               goto 1                                   
-            end if                                     
-         end do                                       
-         write (*,*)  ' error in xtemp(kr) =', xtemp(k),               
-     @        '. it should be between '                    
-         write (*,*)  ' tmin =',tmin, '   and   tmax =',tmax           
-         stop                                         
- 1    continue                                     
-                                                
-      stx = 0.d0                                     
-      if (ts.gt.tmax) ts = dble( tmax )              
-      if (ts.lt.tmin) ts = dble( tmin )              
-      i = 1                                          
-      do while (i.le.mm-1)                           
-         i = i + 1                                    
-!     write(*,*)'mztfsub_overlap/160',i,ts,thist(i)
-         if ( abs(ts-thist(i)) .lt. 1.0d-4 ) then     
-            do k=1,nbox                                
-               stx = stx + no(k) * sk1(i,k)    
-!     write(*,*)'mztfsub_overlap/164',stx
-            end do                                     
-            return                                     
-         elseif ( thist(i) .le. ts ) then             
-            do k=1,nbox                                
-               stx = stx + no(k) * (( sk1(i,k)*(thist(i-1)-ts) +    
-     @              sk1(i-1,k)*(ts-thist(i)) )) / (thist(i-1)-thist(i))      
-!              write(*,*)'mztfsub_overlap/171',stx
-            end do                                     
-!     stop
-            return                                     
-         end if                                       
-      end do  
-                                                
-      return                                        
-      end
-
-                                            
+
+c     arguments
+      real*8 alsx(nbox_max),adx(nbox_max) ! Output
+      real*8 xtemp(nbox_max)    ! Input
+      integer    i              ! I , O
+
+c     local variables
+      integer   k
+      real*8          factor
+      real*8    temperatura     ! para evitar valores ligeramnt out of limits
+
+c     ***********
+
+      do 1, k=1,nbox_max
+         temperatura = xtemp(k)
+         if (abs(xtemp(k)-thist(1)).le.0.01d0) then
+            temperatura=thist(1)
+         elseif (abs(xtemp(k)-thist(nhist)).le.0.01d0) then
+            temperatura=thist(nhist)
+         endif
+         call huntdp ( thist,nhist, temperatura, i )
+         if ( i.eq.0 .or. i.eq.nhist ) then
+            write (*,*) ' HUNT/ Limits input grid:',
+     @           thist(1),thist(nhist)
+            write (*,*) ' HUNT/ location in new grid:', xtemp(k)
+            stop ' INTERSHP/ Interpolation error. T out of Histogram.'
+         endif
+         factor = 1.d0 /  (thist(i+1)-thist(i))
+         alsx(k) = (( xls1(i,k)*(thist(i+1)-xtemp(k)) +
+     @        xls1(i+1,k)*(xtemp(k)-thist(i)) )) * factor
+         adx(k)  = (( xld1(i,k)*(thist(i+1)-xtemp(k)) +
+     @        xld1(i+1,k)*(xtemp(k)-thist(i)) )) * factor
+ 1    continue
+
+      return
+      end
+
 c     **********************************************************************
-      subroutine intz(h,aco2,ap,amr,at, con)         
-c     return interp. concentration, pressure,mixing ratio and temperature   
-c     for a input height h                          
+
+      subroutine interstrhunt (i, stx, ts, sx, xtemp )
+
 c     **********************************************************************
-                                                
-      implicit none                                  
+
+      implicit none
+
       include 'nlte_paramdef.h'
       include 'nlte_commons.h'
-                                                
-c arguments                                     
-      real              h       ! i 
-      real*8            con(nzy) ! i                          
-      real*8            aco2, ap, at, amr ! o                  
-                                                
-c local variables                               
-      integer           k                                     
-                                                
-c     ************                                  
-                                                
-      if ( ( h.lt.zy(1) ).and.( h.le.-1.e-5 ) ) then 
-         write (*,*) ' zp= ',h,' zy(1)= ',zy(1)                         
-         stop'from intz: error in interpolation, z < minimum height'
-      elseif (h.gt.zy(nzy)) then                      
-         write (*,*) ' zp= ',h,' zy(nzy)= ',zy(nzy)                       
-         stop'from intz: error in interpolation, z > maximum height'
-      end if                                         
-                                                
-      if (h.eq.zy(nzy)) then                          
-         ap  = dble( py(nzy)  )                       
-         aco2= con(nzy)                               
-         at  = dble( ty(nzy)  )                         
-         amr = dble( mr(nzy) )                          
-         return                                        
-      end if                                         
-                                                
-      do k=1,nzy-1                                    
-         if( abs( h-zy(k) ).le.( 1.e-5 ) ) then        
-            ap  = dble( py(k)  )                       
-            aco2= con(k)                               
-            at  = dble( ty(k)  )                         
-            amr = dble( mr(k) )                          
-            return                                       
-         elseif(h.gt.zy(k).and.h.lt.zy(k+1))then       
-            ap = dble( exp( log(py(k)) + log(py(k+1)/py(k)) *         
-     @           (h-zy(k)) / (zy(k+1)-zy(k)) ) )             
-            aco2 = exp( log(con(k)) + log( con(k+1)/con(k) ) *        
-     @           (h-zy(k)) / (zy(k+1)-zy(k)) )               
-            at = dble( ty(k)+(ty(k+1)-ty(k))*(h-zy(k))/
-     @           (zy(k+1)-zy(k)) )
-            amr = dble( mr(k)+(mr(k+1)-mr(k))*(h-zy(k))/
-     @           (zy(k+1)-zy(k)) ) 
-            return                                       
-         end if                                        
-      end do                                         
-                                                
-      return                                         
-      end                                            
-                                                
-                                            
-                                                
+
+c     arguments
+      real*8            stx     ! output, total band strength
+      real*8            ts      ! input, temp for stx
+      real*8            sx(nbox_max) ! output, strength for each box
+      real*8            xtemp(nbox_max) ! input, temp for sx
+      integer   i
+
+c     local variables
+      integer   k
+      real*8          factor
+      real*8    temperatura
+
+c     ***********
+
+      do 1, k=1,nbox
+         temperatura = xtemp(k)
+         if (abs(xtemp(k)-thist(1)).le.0.01d0) then
+            temperatura=thist(1)
+         elseif (abs(xtemp(k)-thist(nhist)).le.0.01d0) then
+            temperatura=thist(nhist)
+         endif
+         call huntdp ( thist,nhist, temperatura, i )
+         if ( i.eq.0 .or. i.eq.nhist ) then
+            write(*,*)'HUNT/ Limits input grid:',thist(1),thist(nhist)
+            write(*,*)'HUNT/ location in new grid:',xtemp(k)
+            stop'INTERSTR/1/ Interpolation error. T out of Histogram.'
+         endif
+         factor = 1.d0 /  (thist(i+1)-thist(i))
+         sx(k) = ( sk1(i,k)   * (thist(i+1)-xtemp(k))
+     @        + sk1(i+1,k) * (xtemp(k)-thist(i))  ) * factor
+ 1    continue
+
+
+      temperatura = ts
+      if (abs(ts-thist(1)).le.0.01d0) then
+         temperatura=thist(1)
+      elseif (abs(ts-thist(nhist)).le.0.01d0) then
+         temperatura=thist(nhist)
+      endif
+      call huntdp ( thist,nhist, temperatura, i )
+      if ( i.eq.0 .or. i.eq.nhist ) then
+         write (*,*) ' HUNT/ Limits input grid:',
+     @        thist(1),thist(nhist)
+         write (*,*) ' HUNT/ location in new grid:', ts
+         stop ' INTERSTR/2/ Interpolat error. T out of Histogram.'
+      endif
+      factor = 1.d0 /  (thist(i+1)-thist(i))
+      stx = 0.d0
+      do k=1,nbox
+         stx = stx + no(k) * ( sk1(i,k)*(thist(i+1)-ts) +
+     @        sk1(i+1,k)*(ts-thist(i)) ) * factor
+      end do
+
+
+      return
+      end
+
 c     **********************************************************************
-      real*8 function iaa_we(ig,me,pe,plaux,idummy,nt_local,p_local,
-     $     Desp,wsL)  
-c     icls=5 -->para mztf                           
-c     icls=1,2,3-->para fot, line shape (v=1,l=2,d=3) (only use if wr=2)    
-c     calculates an approximate equivalent width for an error estimate.     
-c                                               
-c     ioverlap = 0  ....... no correction for overlaping        
-c     1  ....... "lisat" first correction (see overlap_box.
-c     2  .......    "      "    "  plus "supersaturation"  
-                                           
-c     idummy=0   do nothing
-c     1   write out some values for diagnostics
-c     2   correct the Strong Lorentz behaviour for SZA>90
-c     3   casos 1 & 2
-     
-c     malv   nov-98    add overlaping's corrections       
+
+      subroutine intzhunt (k, h, aco2,ap,amr,at, con)
+
+c     k lleva la posicion de la ultima llamada a intz , necesario para
+c     que esto represente una aceleracion real.
 c     **********************************************************************
-                                                
-      implicit none                                  
-      
+
+      implicit none
       include 'nlte_paramdef.h'
       include 'nlte_commons.h'
-                                                
-c arguments                 
-      integer         ig        ! ADDED FOR TRACEBACK
-      real*8          me      ! I. path's absorber amount  
-      real*8          pe        ! I. path's presion total
-      real*8          plaux     ! I. path's partial pressure of CO2
-      real*8          nt_local  ! I. needed for strong limit of Lorentz profil
-      real*8          p_local   ! I.    "          "              "
-      integer         idummy    ! I. indica varias opciones
-      real*8          wsL       ! O. need this for strong Lorentz correction
-      real*8          Desp      ! I. need this for strong Lorentz correction
-      
-c local variables                               
-      integer         i                                     
-      real*8          y,x,wl,wd                    
-      real*8          cn(0:7),dn(0:7)                        
-      real*8          pi, xx                                
-      real*8          f_sat_box                      
-      real*8          dv_sat_box, dv_corte_box       
-      real*8          area_core_box, area_wing_box   
-      real*8          wlgood , parentesis , xlor
-      real*8          wsl_grad
-  
-                                                       
-c data blocks                                   
-      data cn/9.99998291698d-1,-3.53508187098d-1,9.60267807976d-2,           
-     @     -2.04969011013d-2,3.43927368627d-3,-4.27593051557d-4,    
-     @     3.42209457833d-5,-1.28380804108d-6/          
-      data dn/1.99999898289,5.774919878d-1,-5.05367549898d-1,    
-     @     8.21896973657d-1,-2.5222672453,6.1007027481, 
-     @     -8.51001627836,4.6535116765/                 
-                                                
-c     ***********                                   
-                                                
-c     equivalent width of atmospheric line.         
-                                                
-      pi = acos(-1.d0)                               
-
-      if ( idummy.gt.9 )
-     @     write (*,*) ' S, m, alsa, pp =', ka(kr), me, alsa(kr), plaux
-      
-      y=ka(kr)*me                                
-!     x=y/(2.0*pi*(alsa(kr)*pl+alna(kr)*(pe-pl)))    
-      x=y/(2.0d0*pi* alsa(kr)*plaux) !+alna(kr)*(pe-pl)))           
-
-! Strong limit of Lorentz profile:  WL = 2 SQRT( S * m * alsa*pl )
-! Para anular esto, comentar las siguientes 5 lineas
-!        if ( x .gt. 1.e6 ) then 
-!           wl = 2.0*sqrt( y * alsa(kr)*pl )
-!        else
-!          wl=y/sqrt(1.0d0+pi*x/2.0d0)                        
-!        endif
-
-      wl=y/sqrt(1.0d0+pi*x/2.0d0)                        
-
-      if (wl .le. 0.d0) then 
-         write(*,*)'mztfsub_overlap/496',ig,y,ka(kr),me,kr
-         stop'WE/Lorentz EQW zero or negative!/498' !,ig
-      endif
-
-      if ( idummy.gt.9 ) 
-     @     write (*,*) ' y, x =', y, x
-
-      xlor = x
-      if ( (idummy.eq.2 .or. idummy.eq.12) .and. xlor.gt.1e5 ) then
-                                         ! en caso que estemos en el regimen
-                                         ! Strong Lorentz y la presion local 
-                                         ! vaya disminuyendo, corregimos la EQW
-                                         ! con un gradiente analitico (notebook)
-         wsL = 2.0*sqrt( y * alsa(kr)*plaux ) 
-         wsl_grad = - 2.d0 * ka(kr)*alsa(kr) * nt_local*p_local / wsL
-         wlgood = w_strongLor_prev(kr) + wsl_grad * Desp
-         if (idummy.eq.12) 
-     @        write (*,*) ' W(wrong), W_SL, W_SL prev, W_SL corrected=', 
-     @        wl, wsL, w_strongLor_prev(kr), wlgood
-         wl = wlgood
-      endif
-        ! wsL = wl  pero esto no lo hacemos todavia, porque necesitamos 
-        !           el valor que ahora mismo tiene wsL para corregir la
-        !           expresion R&W below
-
-!        write (*,*) 'WE arguments me,pe,pl =', me,pe,pl
-!        write (*,*) 'WE/ wl,ka(kr),alsa(kr) =', 
-!     @       wl, ka(kr),alsa(kr)
-
-
-!>>>>>>>
- 500  format (a,i3,3(2x,1pe15.8))
- 600  format (a,2(2x,1pe16.9))
- 700  format (a,3(1x,1pe16.9))
-!        if (kr.eq.8 .or. kr.eq.13) then  
-!           write (*,500) 'WE/kr,m,pt,pl=', kr, me, pe, pl
-!           write (*,700) '  /aln,als,d_x=', alna(kr),alsa(kr),
-!     @                2.0*pi*( alsa(kr)*pl + alna(kr)*(pe-pl) )
-!           write (*,600) '  /alsa*p_CO2, alna*p_n2 :', 
-!     @             alsa(kr)*pl, alna(kr)*(pe-pl)
-!           write (*,600) '  a*p, S =', 
-!     @                 alsa(kr)*pl + alna(kr)*(pe-pl)  , ka(kr)
-!           write (*,600) '  /S*m, x =', y, x
-!           write (*,600) '  /aprox, WL =',  
-!     @         2.*sqrt( y*( alsa(kr)*pl+alna(kr)*(pe-pl) ) ), WL
-!        endif
-        !                                              
-        ! corrections to lorentz eqw due to overlaping and super-saturation    
-        !                                              
-                                                
-      i_supersat = 0                                 
-                                                
-      if ( icls.eq.5 .and. ioverlap.gt.0 ) then      
-           ! for the moment, only consider overlaping for mztf.f, not fot.f    
-                                                
-           ! definition of saturation in the lisat model           
-           !                                           
-         asat_box = 0.99d0                           
-         f_sat_box = 2.d0 * x                        
-         xx = f_sat_box / log( 1./(1-asat_box) )     
-         if ( xx .lt. 1.0d0 ) then                   
-            dv_sat_box = 0.0d0                       
-            asat_box = 1.0d0 - exp( - f_sat_box )    
-         else                                        
-            dv_sat_box = alsa(kr) * sqrt( xx - 1.0d0 )           
-              ! approximation: only use of alsa in mars and venus  
-         endif                                       
-                                                
-           ! area of saturated line                    
-           !                                           
-         area_core_box = 2.d0 * dv_sat_box * asat_box            
-         area_wing_box = 0.5d0 * ( wl - area_core_box )          
-         dv_corte_box = dv_sat_box + 2.d0*area_wing_box/asat_box 
-                                                       
-           ! super-saturation or simple overlaping?    
-           !                                           
-!          i_supersat = 0                             
-         xx = dv_sat_box - ( 0.5d0 * dist(kr) )      
-         if ( xx .ge. 0.0       ! definition of supersaturation  
-     @        .and. dv_sat_box .gt. 0.0 ! definition of saturation 
-     @        .and. (dist(kr).gt.0.0) ) ! box contains more than 1 line 
-     @                              ! and not too far apart       
-     @        then                                         
-                                                
-            i_supersat = 1                           
-                                                
-         else                                        
-           ! no super-saturation, then use "lisat + first correction", i.e.,   
-           ! correct for line products                 
-           !                                           
-                                                
-            wl = wl                                  
-                                                
-         endif                                       
-                                                
-      end if                    ! end of overlaping loop           
-
-      if (icls.eq.2) then 
-         iaa_we = wl              
-         return
-      endif
-
-cc  doppler limit:    
-      if ( idummy.gt.9 ) 
-     @     write (*,*) ' S*m, alf_dop =', y, alda(kr)*sqrt(pi)
-
-      x = y / (alda(kr)*sqrt(pi))  
-      if ( x.lt.1.e-10 ) then   ! to avoid underflow
-         wd = y
-      else
-         wd=alda(kr)*sqrt(4.0*pi*x**2*(1.0+log(1.0+x))/(4.0+pi*x**2))
-      endif
-      if ( idummy.gt.9 )
-     @     write (*,*) ' wd =', wd
-                                     
-cc  doppler weak limit                          
-c       wd = ka(kr) * me                              
-                                                
-cc  good doppler                                
-      if(icls.eq.5) then        !para mztf                  
-         !write (*,*) 'para mztf, icls=',icls                           
-         if (x.lt.5.) then                             
-            wd = 0.d0                                    
-            do i=0,7                                     
-               wd = wd + cn(i) * x**i                     
-            end do                                       
-            wd = alda(kr) * x * sqrt(pi) * wd            
-         elseif (x.gt.5.) then                         
-            wd = 0.d0                                    
-            do i=0,7                                     
-               wd = wd + dn(i) / (log(x))**i              
-            end do                                       
-            wd = alda(kr) * sqrt(log(x)) * wd            
-         else                                          
-            stop ' x should not be less than zero'        
-         end if                                        
-      end if                                         
-                                                
-
-      if ( i_supersat .eq. 0 ) then                  
-
-         parentesis = wl**2+wd**2-(wd*wl/y)**2
-                                ! changed +(wd*wl/y)**2 to -...14-3-84      
-
-         if ( parentesis .lt. 0.0 ) then 
-            if ((idummy.eq.2 .or. idummy.eq.12) .and. xlor.gt.1e5) then
-               parentesis = wl**2+wd**2-(wd*wsL/y)**2
-                                ! este cambio puede ser necesario cuando se hace
-                                ! correccion Strong Lor, para evitar valores
-                                ! negativos del parentesis en sqrt( )
-            else
-               stop ' WE/ Error en las EQW  wl,wl,y '
-            endif
-         endif
-
-         iaa_we = sqrt( parentesis )
-!          write (*,*)  ' from iaa_we: xdop,alda,wd', sngl(x),alda(kr),sngl(wd)
-!          write (*,*)  ' from iaa_we: we', iaa_we                             
-
-      else                                           
-
-         iaa_we = wl                                     
-          ! if there is supersaturation we can ignore wd completely;           
-          ! mztf.f will compute the eqw of the whole box afterwards            
-
-      endif                                          
-                                                
-      if (icls.eq.3) iaa_we = wd                         
-      
-      if ( idummy.gt.9 )
-     @     write (*,*) ' wl,wd,w =', wl,wd,iaa_we
-      
-      wsL = wl 
-
-      return                                         
-      end                                            
-                                                
-      
+
+c     arguments
+      real              h       ! i
+      real*8            con(nzy) ! i
+      real*8            aco2, ap, at, amr ! o
+      integer           k       ! i
+
+c     local variables
+      real          factor
+
+c     ************
+
+      call hunt ( zy,nzy, h, k )
+      factor =  (h-zy(k)) /  (zy(k+1)-zy(k))
+      ap = dble( exp( log(py(k)) + log(py(k+1)/py(k)) * factor ) )
+      aco2 = dlog(con(k)) + dlog( con(k+1)/con(k) ) * dble(factor)
+      aco2 = exp( aco2 )
+      at = dble( ty(k) + (ty(k+1)-ty(k)) * factor )
+      amr = dble( mr(k) + (mr(k+1)-mr(k)) * factor )
+
+
+      return
+      end
+
 c     **********************************************************************
-      real*8 function simrul(a,b,fsim,c,acc)         
-c     adaptively integrates fsim from a to b, within the criterion acc.     
+
+      subroutine intzhunt_cts (k, h, nzy_cts_real, 
+     @     aco2,ap,amr,at, con)
+
+c     k lleva la posicion de la ultima llamada a intz , necesario para
+c     que esto represente una aceleracion real.
 c     **********************************************************************
-        
-      implicit none
-                          
-      real*8 res,a,b,g0,g1,g2,g3,g4,d,a0,a1,a2,h,x,acc,c,fsim
-      real*8 s1(70),s2(70),s3(70)
-      real*8 c1, c2
-      integer*4 m,n,j                                
-                                                
-      res=0.                                         
-      c=0.                                           
-      m=0                                            
-      n=0                                            
-      j=30                                           
-      g0=fsim(a)                                     
-      g2=fsim((a+b)/2.)                              
-      g4=fsim(b)                                     
-      a0=(b-a)*(g0+4.0*g2+g4)/2.0                    
- 1    d=2.0**n                                  
-      h=(b-a)/(4.0*d)                                
-      x=a+(4.0*m+1.0)*h                              
-      g1=fsim(x)                                     
-      g3=fsim(x+2.0*h)                               
-      a1=h*(g0+4.0*g1+g2)                            
-      a2=h*(g2+4.0*g3+g4)                            
-      if ( abs(a1+a2-a0).gt.(acc/d)) goto 2          
-      res=res+(16.0*(a1+a2)-a0)/45.0                 
-      m=m+1                                          
-      c=a+m*(b-a)/d                                  
- 6    if (m.eq.(2*(m/2))) goto 4                
-      if ((m.ne.1).or.(n.ne.0)) goto 5               
- 8    simrul=res                                
-      return                                         
- 2    m=2*m                                     
-      n=n+1                                          
-      if (n.gt.j) goto 3                             
-      a0=a1                                          
-      s1(n)=a2                                       
-      s2(n)=g3                                       
-      s3(n)=g4                                       
-      g4=g2                                          
-      g2=g1                                          
-      goto 1                                         
- 3    c1=c-(b-a)/d                              
-      c2=c+(b-a)/d                                   
-      write(2,7) c1,c,c2,fsim(c1),fsim(c),fsim(c2)    
-      write(*,7) c1,c,c2,fsim(c1),fsim(c),fsim(c2)      
- 7    format(2x,'17hsimrule fails at ',/,3e15.6,/,3e15.6)     
-      goto 8                                         
- 5    a0=s1(n)                                  
-      g0=g4                                          
-      g2=s2(n)                                       
-      g4=s3(n)                                       
-      goto 1                                         
- 4    m=m/2                                     
-      n=n-1                                          
-      goto 6                                         
-      end                                            
-                                                
-c     **********************************************************************
-      subroutine findw(ig,iirw,idummy,c1,p1, Desp, wsL)                         
-c     this routine sets up accuracy criteria and calls simrule between limit
-c     that depend on the number of atmospheric and cell paths. it gives eqw.
-
-c     Add correction for EQW in Strong Lorentz regime and SZA>90 
-c     **********************************************************************
-                                                
-      implicit none                                  
+
+      implicit none
       include 'nlte_paramdef.h'
       include 'nlte_commons.h'
-                                                
-c arguments                
-      integer         ig        ! ADDED FOR TRACEBACK
-      integer           iirw        
-      integer         idummy    ! I. indica varias opciones
-      real*8          c1        ! I. needed for strong limit of Lorentz profil
-      real*8          p1        ! I.    "          "              "
-      real*8          wsL       ! O. need this for strong Lorentz correction
-      real*8          Desp      ! I. need this for strong Lorentz correction
-      
-c local variables                               
-      real*8            ept,eps,xa                            
-      real*8            acc,  c                                
-      real*8            iaa_we                                    
-      real*8            iaa_f, iaa_fi, simrul                          
-                                                
-      external iaa_f,iaa_fi                                  
-                                                
-c       ********** *********** *********                                     
-
-      if(icls.eq.5) then        !para mztf                  
-!           if(ig.eq.1682)write(*,*)'mztfsub_overlap/768',ua(kr),iirw
-         if (iirw.eq.2) then    !iirw=icf=2 ==> we use the w&r formula      
-            w = iaa_we(ig,ua(kr),pt,pp, idummy, c1,p1, Desp, wsL )  
-            return                                        
-         end if                                      
-         ept=iaa_we(ig,ua(kr),pt,pp, idummy,c1,p1, Desp, wsL)
-      else                      !para fot               
-         if (iirw.eq.2) then    ! icf=2 ==> we use the w&r formula 
-            w = iaa_we(ig,sl_ua,pt,pp, idummy,c1,p1, Desp, wsL)
-            return                                        
-         end if                                      
-         ept=iaa_we(ig,sl_ua,pt,pp, idummy,c1,p1, Desp, wsL)
-      end if                                         
-                                                
-c the next block is a modification to avoid nul we.         
-c this situation appears for weak lines and low path temperature, but   
-c there is not any loss of accuracy. first july 1986        
-      if (ept.eq.0.) then       ! for weak lines sometimes we=0       
-         ept=1.0e-18                                   
-         write (*,*)  'ept =',ept                                       
-         write (*,*) 'from we: we=0.0'                                  
-         return                                        
-      end if                                         
-                                                
-      acc = 4.d0                                     
-      acc = 10.d0**(-acc)                            
-                                                
-      eps = acc * ept           !accuracy 10-4 atmospheric eqw.  
-      xa=0.5*ept/iaa_f(0.d0)        !width of doppler shifted atmospheric line.    
-      w = 2.0*( simrul(0.0d0,xa,iaa_f,c,eps)
-     .        + simrul(0.1d0,1.0/xa,iaa_fi,c,eps) )
-!no shift.                                      
-                                                
-      return                                     
-      end                                            
-                                                
-                                                
+
+c     arguments
+      real              h       ! i
+      real*8            con(nzy_cts) ! i
+      real*8            aco2, ap, at, amr ! o
+      integer           k       ! i
+      integer         nzy_cts_real ! i
+
+c     local variables
+      real          factor
+
+c     ************
+
+      call hunt_cts ( zy_cts,nzy_cts, nzy_cts_real, h, k )
+      factor =  (h-zy_cts(k)) /  (zy_cts(k+1)-zy_cts(k))
+      ap = dble( exp( log(py_cts(k)) +
+     @     log(py_cts(k+1)/py_cts(k)) * factor ) )
+      aco2 = dlog(con(k)) + dlog( con(k+1)/con(k) ) * dble(factor)
+      aco2 = exp( aco2 )
+      at = dble( ty_cts(k) + (ty_cts(k+1)-ty_cts(k)) * factor )
+      amr = dble( mr_cts(k) + (mr_cts(k+1)-mr_cts(k)) * factor )
+
+
+      return
+      end
+
+
 c     **********************************************************************
-      double precision function iaa_fi(y)                
-c     returns the value of f(1/y)                   
+
+      real*8 function we_clean  ( y,pl, xalsa, xalda )
+
 c     **********************************************************************
-                                                
-      implicit none                                  
-      real*8 iaa_f, y                                    
-                                                
-      iaa_fi=iaa_f(1.0/y)/y**2                               
-      return                                         
-      end                                            
-                                                
-                                                
-c     **********************************************************************
-      double precision function iaa_f(nuaux)                
-c     calculates 1-exp(-k(nu)u) for all series paths or combinations thereof
-c     **********************************************************************
-                                                
-      implicit none                                  
+
+      implicit none
+
+      include 'nlte_paramdef.h'
+
+c     arguments
+      real*8            y       ! I. path's absorber amount * strength
+      real*8          pl        ! I. path's partial pressure of CO2
+      real*8          xalsa     ! I.  Self lorentz linewidth for 1 isot & 1 box
+      real*8          xalda     ! I.  Doppler linewidth        "           "
+
+c     local variables
+      integer   i
+      real*8            x,wl,wd,wvoigt
+      real*8            cn(0:7),dn(0:7)
+      real*8          factor, denom
+      real*8          pi, pi2, sqrtpi
+
+c     data blocks
+      data cn/9.99998291698d-1,-3.53508187098d-1,9.60267807976d-2,
+     @     -2.04969011013d-2,3.43927368627d-3,-4.27593051557d-4,
+     @     3.42209457833d-5,-1.28380804108d-6/
+      data dn/1.99999898289,5.774919878d-1,-5.05367549898d-1,
+     @     8.21896973657d-1,-2.5222672453,6.1007027481,
+     @     -8.51001627836,4.6535116765/
+
+c     ***********
+
+      pi = 3.141592
+      pi2= 6.28318531
+      sqrtpi = 1.77245385
+
+      x=y / ( pi2 * xalsa*pl )
+
+
+c     Lorentz
+      wl=y/sqrt(1.0d0+pi*x/2.0d0)
+
+c     Doppler
+      x = y / (xalda*sqrtpi)
+      if (x .lt. 5.0d0) then
+         wd = cn(0)
+         factor = 1.d0
+         do i=1,7
+            factor = factor * x
+            wd = wd + cn(i) * factor
+         end do
+         wd = xalda * x * sqrtpi * wd
+      else
+         wd = dn(0)
+         factor = 1.d0 / log(x)
+         denom = 1.d0
+         do i=1,7
+            denom = denom * factor
+            wd = wd + dn(i) * denom
+         end do
+         wd = xalda * sqrt(log(x)) * wd
+      end if
+
+c     Voigt
+      wvoigt = wl*wl + wd*wd - (wd*wl/y)*(wd*wl/y)
+
+      if ( wvoigt .lt. 0.0d0 ) then
+         write (*,*) ' Subroutine WE/ Error in Voift EQS calculation '
+         write (*,*) '  WL, WD, X, Y = ', wl, wd, x, y
+         stop '  ERROR : Imaginary EQW. Revise spectral data. '
+      endif
+
+      we_clean = sqrt( wvoigt )
+
+
+      return
+      end
+
+
+c     ***********************************************************************
+
+      subroutine mztf_correccion (coninf, con, ib )
+
+c     ***********************************************************************
+
+      implicit none
+
       include 'nlte_paramdef.h'
       include 'nlte_commons.h'
-      
-      double precision tra,xa,ya,za,yy,nuaux
-      double precision voigtf                        
-      tra=1.0d0                                      
-                                                
-      yy=1.0d0/alda(kr)                          
-      xa=nuaux*yy                                       
-      ya= ( alsa(kr)*pp + alna(kr)*(pt-pp) ) * yy                        
-      za=ka(kr)*yy                                   
-                                                
-      if(icls.eq.5) then        !para mztf                  
-          ! write (*,*) 'icls=',icls                                    
-         tra=za*ua(kr)*voigtf(sngl(xa),sngl(ya))     
+
+c     arguments
+      integer           ib
+      real*8            con(nzy), coninf
+
+!     local variables
+      integer   i, isot
+      real*8    tvt0(nzy), tvtbs(nzy), zld(nl),zyd(nzy)
+      real*8  xqv, xes, xlower, xfactor
+
+c     *********
+
+      isot = 1
+      nu11 = dble( nu(1,1) )
+
+      do i=1,nzy
+         zyd(i) = dble(zy(i))
+      enddo
+      do i=1,nl
+         zld(i) = dble( zl(i) )
+      end do
+
+!     tvtbs
+      call interhuntdp (tvtbs,zyd,nzy, v626t1,zld,nl, 1 )
+
+!     tvt0
+      if (ib.eq.2 .or. ib.eq.3 .or. ib.eq.4) then
+         call interhuntdp (tvt0,zyd,nzy, v626t1,zld,nl, 1 )
       else
-         tra=za*sl_ua*voigtf(sngl(xa),sngl(ya))         
-      end if                                         
-                                                
-      if (tra.gt.50.0) then                          
-         tra=1.0                !2.0e-22 overflow cut-off.          
-      else if (tra.gt.1.0e-4) then                   
-         tra=1.0-exp(-tra)                              
-      end if                                         
-                                                
-      iaa_f=tra                                          
-      return                                         
-      end                                            
-                                                
-c     **********************************************************************
-      double precision function voigtf(x1,y)         
-c     computes voigt function for any value of x1 and any +ve value of y.   
-c     where possible uses modified lorentz and modified doppler approximatio
-c     otherwise uses a rearranged rybicki routine.  
-c     c(n) = exp(-(n/h)**2)/(pi*sqrt(pi)), with h = 2.5 .       
-c     accurate to better than 1 in 10000.           
-c     **********************************************************************
-
-      implicit none
-      
-      real x1, y
-      real x, xx, xxyy, xh,xhxh, yh,yhyh, f1,f2, p, q, xn,xnxn, voig
-      
-      real*8 b,g0,g1,g2,g3,g4,d1,d2,d3,d4,c          
-      integer*4 n                             
-                                                
-      dimension c(10)                                
-      complex xp,xpp,z                               
-                                                
-      data c(1)/0.15303405/                          
-      data c(2)/0.94694928e-1/                       
-      data c(3)/0.42549174e-1/                       
-      data c(4)/0.13882935e-1/                       
-      data c(5)/0.32892528e-2/                       
-      data c(6)/0.56589906e-3/                       
-      data c(7)/0.70697890e-4/                       
-      data c(8)/0.64135678e-5/                       
-      data c(9)/0.42249221e-6/                       
-      data c(10)/0.20209868e-7/                      
-                                                
-      x=abs(x1)                                      
-      if (x.gt.7.2) goto 1                           
-      if ((y+x*0.3).gt.5.4) goto 1                   
-      if (y.gt.0.01) goto 3                          
-      if (x.lt.2.1) goto 2                           
-      goto 3                                         
-c     here uses modified lorentz approx.            
- 1    xx=x*x                                    
-      xxyy=xx+y*y                                    
-      b=xx/xxyy                                      
-      voigtf=y*(1.+(2.*b-0.5+(0.75-(9.-12.*b)*b)/xxyy)/          
-     *     xxyy)/(xxyy*3.141592654)                 
-      return                                         
-c     here uses modified doppler approx.            
- 2    xx=x*x                                    
-      voigtf=0.56418958*exp(-xx)*(1.-y*(1.-0.5*y)*(1.1289-xx*(1.1623+        
-     *     xx*(0.080812+xx*(0.13854-xx*(0.033605-0.0073972*xx))))))         
-      return                                         
-c     here uses a rearranged rybicki routine.       
- 3    xh=2.5*x                                  
-      xhxh=xh*xh                                     
-      yh=2.5*y                                       
-      yhyh=yh*yh                                     
-      f1=xhxh+yhyh                                   
-      f2=f1-0.5*yhyh                                 
-      if (y.lt.0.1) goto 20                          
-      p=-y*7.8539816            !7.8539816=2.5*pi            
-      q=x*7.8539816                                  
-      xpp=cmplx(p,q)                                 
-      z=cexp(xpp)                                    
-      d1=xh*aimag(z)                                 
-      d2=-d1                                         
-      d3=yh*(1.-real(z))                             
-      d4=-d3+2.*yh                                   
-      voig=0.17958712*(d1+d3)/f1                     
-      goto 30                                        
- 20   p=x*7.8539816                             
-      q=y*7.8539816                                  
-      xp=cmplx(p,q)                                  
-      z=ccos(xp)                                     
-      d1=xh*aimag(z)                                 
-      d2=-d1                                         
-      d3=yh*(1.-real(z))                             
-      d4=-d3+2.*yh                                   
-      voig=0.56418958*exp(y*y-x*x)*cos(2.*x*y)+0.17958712*(d1+d3)/f1         
- 30   xn=0.                                     
-      do 55 n=1,10,2                             
-         xn=xn+1.                                   
-         xnxn=xn*xn                                 
-         g1=xh-xn                                   
-         g2=g1*(xh+xn)                              
-         g3=0.5*g2*g2                               
-         voig=voig+c(n)*(d2*(g2+yhyh)+d4*(f1+xnxn))/
-     &        (g3+yhyh*(f2+xnxn))     
-         xn=xn+1.                                   
-         xnxn=xn*xn                                 
-         g1=xh-xn                                   
-         g2=g1*(xh+xn)                              
-         g3=0.5*g2*g2                               
-         voig=voig+c(n+1)*(d1*(g2+yhyh)+d3*(f1+xnxn))/
-     @        (g3+yhyh*(f2+xnxn))   
- 55   continue                              
-      voigtf=voig                                    
-      return                                         
-      end  
-
-
-
-c     **********************************************************************
-c     elimin_mz1d.F (includes smooth_cf)
-c     ************************************************************************
-      subroutine elimin_mz1d (c,vc, ilayer,nanaux,itblout, nwaux)
-
-c     Eliminate anomalous negative numbers in c(nl,nl) according to "nanaux":
-
-c     nanaux = 0 -> no eliminate 
-c              @       -> eliminate all numbers with absol.value<abs(max(c(n,r)))/300.
-c              2 -> eliminate all anomalous negative numbers in c(n,r).
-c              3 -> eliminate all anomalous negative numbers far from the main
-c                   diagonal. 
-c              8 -> eliminate all non-zero numbers outside the main diagonal,
-c                   and the contibution of lower boundary. 
-c              9 -> eliminate all non-zero numbers outside the main diagonal. 
-c              4 -> hace un smoothing cuando la distancia de separacion entre
-c                   el valor maximo y el minimo de cf > 50 capas.
-c              5 -> elimina valores menores que 1.0d-19
-c              6 -> incluye los dos casos 4 y 5 
-c              7 -> llama a lisa: smooth con width=nw & elimina mejorado
-c              78-> incluye los dos casos 7 y 8
-c              79-> incluye los dos casos 7 y 9
-
-c     itblout (itableout in calling program) is the option for writing 
-c     out or not the purged c(n,r) matrix: 
-c     itblout = 0 -> no write  
-c               = 1 -> write out in curtis***.out according to ilayer
-
-c     ilayer is the index for the layer selected to write out the matrix:
-c     ilayer = 0  => matrix elements written out cover all the altitudes
-c                                                     with 5 layers steps
-c              > 0  =>   "        "      "      "  are  c(ilayer,*)
-c     NOTA: 
-c       EXISTE LA POSIBILIDAD DE SACAR TODAS LAS CAPAS (TODA LA MATRIZ)
-c       UTILIZANDO itableout=30 EN MZTUD
-
-c     jul 2011        malv+fgg       adapted to LMD-MGCM
-c     Sep-04          FGG+MALV        Correct include and call parameters
-c     cristina  25-sept-1996   y  27-ene-1997
-c     JAN 98            MALV            Version for mz1d
-c     ************************************************************************
+         do i=1,nzy
+            tvt0(i) = dble( ty(i) )
+         end do
+      end if
+
+c     factor
+      do i=1,nzy
+
+         xlower = exp( ee*dble(elow(isot,ib)) *
+     @        ( 1.d0/dble(ty(i))-1.d0/tvt0(i) ) )
+         xes = 1.0d0
+         xqv = ( 1.d0-exp( -ee*nu11/tvtbs(i) ) ) /
+     @        (1.d0-exp( -ee*nu11/dble(ty(i)) ))
+         xfactor = xlower * xqv**2.d0 * xes
+
+         con(i) = con(i) * xfactor
+         if (i.eq.nzy) coninf = coninf * xfactor
+
+      end do
+
+
+      return
+      end
+
+
+c     ***********************************************************************
+
+      subroutine mzescape_normaliz ( taustar, istyle )
+
+c     ***********************************************************************
+
+      implicit none
+      include 'nlte_paramdef.h'
+
+c     arguments
+      real*8            taustar(nl) ! o
+      integer         istyle    ! i
+
+c     local variables and constants
+      integer   i, imaximum
+      real*8          maximum
+
+c     ***************
+
+      taustar(nl) = taustar(nl-1)
+
+      if ( istyle .eq. 1 ) then
+         imaximum = nl
+         maximum = taustar(nl)
+         do i=1,nl-1
+            if (taustar(i).gt.maximum) taustar(i) = taustar(nl)
+         enddo
+      elseif ( istyle .eq. 2 ) then
+         imaximum = nl
+         maximum = taustar(nl)
+         do i=nl-1,1,-1
+            if (taustar(i).gt.maximum) then
+               maximum = taustar(i)
+               imaximum = i
+            endif
+         enddo
+         do i=imaximum,nl
+            if (taustar(i).lt.maximum) taustar(i) = maximum
+         enddo
+      endif
+
+      do i=1,nl
+         taustar(i) = taustar(i) / maximum
+      enddo
+
+
+c     end
+      return
+      end
+
+c     ***********************************************************************
+
+      subroutine mzescape_normaliz_02 ( taustar, nn, istyle )
+
+c     ***********************************************************************
+
+      implicit none
+
+c     arguments
+      real*8            taustar(nn) ! i,o
+      integer         istyle    ! i
+      integer         nn        ! i
+
+c     local variables and constants
+      integer   i, imaximum
+      real*8          maximum
+
+c     ***************
+
+      taustar(nn) = taustar(nn-1)
+
+      if ( istyle .eq. 1 ) then
+         imaximum = nn
+         maximum = taustar(nn)
+         do i=1,nn-1
+            if (taustar(i).gt.maximum) taustar(i) = taustar(nn)
+         enddo
+      elseif ( istyle .eq. 2 ) then
+         imaximum = nn
+         maximum = taustar(nn)
+         do i=nn-1,1,-1
+            if (taustar(i).gt.maximum) then
+               maximum = taustar(i)
+               imaximum = i
+            endif
+         enddo
+         do i=imaximum,nn
+            if (taustar(i).lt.maximum) taustar(i) = maximum
+         enddo
+      endif
+
+      do i=1,nn
+         taustar(i) = taustar(i) / maximum
+      enddo
+
+
+c     end
+      return
+      end
+
+
+c     *** interdp_ESCTVCISO_dlvr11.f ***
+
+c***********************************************************************
+
+      subroutine interdp_ESCTVCISO
+
+c***********************************************************************
 
       implicit none
@@ -833 +622 @@
       include 'nlte_commons.h'
 
-      integer   nanaux,j,i,itblout,kk,k,ir,in
-      integer   ilayer,jmin, jmax,np,nwaux,ntimes,ntimes2
-!*      real*8    c(nl,nl), vc(nl), amax, cmax, cmin, cs(nl,nl), mini
-      real*8    c(nl,nl), vc(nl), amax, cmax, cmin, mini
-      real*8 aux(nl), auxs(nl)
-      character layercode*3
-
-      ntimes=0
-      ntimes2=0
-!       type *,'from elimin_mz4: nan, nw',nan,nw
-
-      if (nanaux .eq. 0) goto 200 
-
-      if(nanaux.eq.1)then
-         do i=1,nl
-            amax=1.0d-36
-            do j=1,nl
-               if(abs(c(i,j)).gt.amax)amax=abs(c(i,j))
-            end do
-            do j=1,nl
-               if(abs(c(i,j)).lt.amax/300.0d0)c(i,j)=0.0d0
-            end do
-         enddo
-      elseif(nanaux.eq.2)then
-         do i=1,nl
-            do j=1,nl
-               if( ( j.le.(i-2) .or. j.gt.(i+2) ).and.
-     @              ( c(i,j).lt.0.0d0 ) ) c(i,j)=0.0d0
-            end do
-         enddo
-      elseif(nanaux.eq.3)then
-         do i=1,nl
-            do j=1,nl
-               if (abs(i-j).ge.10) c(i,j)=0.0d0
-            end do
-         enddo
-      elseif(nanaux.eq.8)then
-         do i=1,nl
-            do j=1,i-1 
-               c(i,j)=0.0d0
-            enddo
-            do j=i+1,nl 
-               c(i,j)=0.0d0
-            enddo
-            vc(i)= 0.d0
-         enddo
-      elseif(nanaux.eq.9)then
-         do i=1,nl
-            do j=1,i-1 
-               c(i,j)=0.0d0
-            enddo
-            do j=i+1,nl 
-               c(i,j)=0.0d0
-            enddo
-         enddo
-!       elseif(nan.eq.7.or.nan.eq.78.or.nan.eq.79)then
-!               call lisa(c, vc, nl, nw)
-      end if
-      if(nanaux.eq.78)then
-         do i=1,nl
-            do j=1,i-1 
-               c(i,j)=0.0d0
-            enddo
-            do j=i+1,nl 
-               c(i,j)=0.0d0
-            enddo
-            vc(i)= 0.d0
-         enddo
-      endif
-      if(nanaux.eq.79)then
-         do i=1,nl
-            do j=1,i-1 
-               c(i,j)=0.0d0
-            enddo
-            do j=i+1,nl 
-               c(i,j)=0.0d0
-            enddo
-         enddo
-      endif
-
-      if(nanaux.eq.5.or.nanaux.eq.6)then
-         do i=1,nl
-            mini = 1.0d-19
-            do j=1,nl
-               if(abs(c(i,j)).le.mini.and.c(i,j).ne.0.d0) then
-                  ntimes2=ntimes2+1
-               end if
-               if ( abs(c(i,j)).le.mini) c(i,j)=0.d0
-            end do
-         enddo
-      end if
-
-      if(nanaux.eq.4.or.nanaux.eq.6)then
-         do i=1,nl
-            do j=1,nl
-               aux(j)=c(i,j)
-               auxs(j)=c(i,j)
-            end do
-                        !call maxdp_2(aux,nl,cmax,jmax)
-            cmax=maxval(aux)
-            jmax=maxloc(aux,dim=1)
-            if(abs(jmax-i).ge.50) then
-               call smooth_cf(aux,auxs,i,nl,3)
-                                !!!call smooth_cf(aux,auxs,i,nl,5)
-               ntimes=ntimes+1
-            end if
-            do j=1,nl
-               c(i,j)=auxs(j)
-            end do
-         end do
-      end if
-
-!          type *, 'elimin_mz4: c(n,r) procesed for elimination. '
-!          type *, ' '
-!          if(nan.eq.4.or.nan.eq.6) type *, '    call smoothing:',ntimes
-!          if(nan.eq.5.or.nan.eq.6) type *, '    call elimina:  ',ntimes2
-!          if(nan.eq.7)   type *, '    from elimin: lisa w=',nw
-!          type *, ' '
-
-
- 200  continue
-
-c       writting out of c(n,r) in ascii file 
-
-!       if(itblout.eq.1) then
-
-!         if (ilayer.eq.0) then
-
-!          open (unit=2, status='new', 
-!     @    file=dircurtis//'curtis_gnu.out', recl=1024)
-!           write(2,'(a)') 
-!     @    ' curtis matrix:     table with   1.e+7 * acf(n,r) '
-!           write(2,114) 'n,r', ( i, i=nl,1,-5 )
-!           do in=nl,1,-5
-!             write(2,*)
-!             write(2,115) in, ( c(in,ir)*1.d7, ir=nl,1,-5 )
-!           end do 
-!          close(2)
-
-
-!          write (*,*)  ' '
-!          write (*,*)  '  curtis.out has been created. '
-!          write (*,*)  ' '
-
-!         else
-
-!            write (layercode,132) ilayer
-!           open (2, status='new', 
-!     @    file=dircurtis//'curtis'//layercode//'.out')
-!           write(2,'(a)') 
-!     @    ' curtis matrix:     table with   1.e+7 * acf(n,r) '
-!           write(2,116) ' layer x       c(',layercode,
-!     @    ',x)           c(x,', layercode,')' 
-!           do in=nl,1,-1
-!            if (c(ilayer,ilayer).ne.0.d0) then 
-!             write(2,117) in, c(ilayer,in), c(in,ilayer), 
-!     @        c(ilayer,in)/c(ilayer,ilayer),
-!     @        c(in,ilayer)/c(ilayer,ilayer)
-!            else
-!             write(2,118) in, c(ilayer,in), c(in,ilayer)
-!            end if
-!           end do 
-!           close(2)
-!           write (*,*)  ' '
-!           write (*,*)  dircurtis//'curtis'//layercode//'.out', 
-!     @ ' has been created.'
-!           write (*,*) ' '
-
-!         end if
-
-!       elseif(itblout.eq.0)then
-
-!         continue
-
-!       else
-
-!         write (*,*) ' error from elimin: ', 
-!     @      ' itblout should be 1 or 0;   itblout= ',itblout
-!         stop
-
-!       end if
-        
-      return
-
- 112  format(10x,10(i3,9x))
- 113  format(1x,i3,2x,9(1pe9.2,2x))
+c     local variables
+      integer    i
+      real*8     lnpnb(nl)
+
+
+c***********************************************************************
+
+c     Use pressure in the NLTE grid but in log and in nb
+      do i=1,nl
+         lnpnb(i) = log( dble( pl(i) * 1013.25 * 1.e6) )
+      enddo
+
+c     Interpolations
+
+      call interhuntdp3veces
+     @     ( taustar21,taustar31,taustar41,    lnpnb, nl,
+     @     tstar21tab,tstar31tab,tstar41tab, lnpnbtab, nztabul,
+     @     1 )
+
+      call interhuntdp3veces ( vc210,vc310,vc410, lnpnb, nl,
+     @     vc210tab,vc310tab,vc410tab, lnpnbtab, nztabul, 2 )
+
+c     end
+      return
+      end
+
+
+c     *** hunt_cts.f ***
+
+cccc  
+      SUBROUTINE hunt_cts(xx,n,n_cts,x,jlo)  
+c     
+c     La dif con hunt es el uso de un indice superior (n_cts) mas bajito que (n)
+c     
+c     Arguments
+      INTEGER jlo               ! O
+      INTEGER n                 ! I
+      INTEGER n_cts             ! I
+      REAL  xx(n)               ! I
+      REAL  x                   ! I
+
+c     Local variables
+      INTEGER inc,jhi,jm  
+      LOGICAL ascnd  
+c     
+cccc  
+c     
+      ascnd=xx(n_cts).ge.xx(1)  
+      if(jlo.le.0.or.jlo.gt.n_cts)then  
+         jlo=0  
+         jhi=n_cts+1  
+         goto 3  
+      endif  
+      inc=1  
+      if(x.ge.xx(jlo).eqv.ascnd)then  
+ 1       jhi=jlo+inc  
+!     write (*,*) jlo
+         if(jhi.gt.n_cts)then  
+            jhi=n_cts+1  
+!     write (*,*) jhi-1
+         else if(x.ge.xx(jhi).eqv.ascnd)then  
+            jlo=jhi  
+            inc=inc+inc  
+!     write (*,*) jlo
+            goto 1  
+         endif  
+      else  
+         jhi=jlo  
+ 2       jlo=jhi-inc  
+!     write (*,*) jlo
+         if(jlo.lt.1)then  
+            jlo=0  
+         else if(x.lt.xx(jlo).eqv.ascnd)then  
+            jhi=jlo  
+            inc=inc+inc  
+            goto 2  
+         endif  
+      endif  
+ 3    if(jhi-jlo.eq.1)then  
+         if(x.eq.xx(n_cts))jlo=n_cts-1  
+         if(x.eq.xx(1))jlo=1  
+!     write (*,*) jlo
+         return  
+      endif  
+      jm=(jhi+jlo)/2  
+      if(x.ge.xx(jm).eqv.ascnd)then  
+         jlo=jm  
+      else  
+         jhi=jm  
+      endif  
+!     write (*,*) jhi-1
+      goto 3  
+c     
+      END  
+
       
- 114  format(1x,a3, 11(8x,i3))
- 115  format( 1x,i3, 2x, 11(1pe10.3))
- 116  format( 1x,a17,a2,a18,a2,a1 ) 
- 117  format( 3x,i3, 4(8x,1pe10.3) )
- 118  format( 3x,i3, 2(8x,1pe10.3) )
- 120  format( 1x,i3, 1x,i3, 2x, 11(1pe10.3))
-
- 132  format(i3)
-
-!  cambio: los formatos 114, 115 , 117 y 118
-!  cambio: al cambia nl de 51 a 140 hay que cambiar el formato i2-->i3
-!          y ahora en vez de 11 capas de 5 en 5, hay 28
-!
-      end
-c**************************************************************************
-      subroutine smooth_cf( c, cs, i, nl, w )
-c     hace un smoothing de c(i,*), de la contribucion de todas las capas
-c     menos de la capa en cuestion, la i.
-c     opcion w (width): el tamanho de la ventana del smoothing.
-c     output values: cs
-c**************************************************************************
-
-      implicit none
+c     *** huntdp.f ***
+
+cccc  
+      SUBROUTINE huntdp(xx,n,x,jlo)  
+c     
+c     Arguments
+      INTEGER jlo               ! O
+      INTEGER n                 ! I
+      REAL*8  xx(n)             ! I
+      REAL*8  x                 ! I
+
+c     Local variables
+      INTEGER inc,jhi,jm  
+      LOGICAL ascnd  
+c     
+cccc  
+c     
+      ascnd=xx(n).ge.xx(1)  
+      if(jlo.le.0.or.jlo.gt.n)then  
+         jlo=0  
+         jhi=n+1  
+         goto 3  
+      endif  
+      inc=1  
+      if(x.ge.xx(jlo).eqv.ascnd)then  
+ 1       jhi=jlo+inc  
+         if(jhi.gt.n)then  
+            jhi=n+1  
+         else if(x.ge.xx(jhi).eqv.ascnd)then  
+            jlo=jhi  
+            inc=inc+inc  
+            goto 1  
+         endif  
+      else  
+         jhi=jlo  
+ 2       jlo=jhi-inc  
+         if(jlo.lt.1)then  
+            jlo=0  
+         else if(x.lt.xx(jlo).eqv.ascnd)then  
+            jhi=jlo  
+            inc=inc+inc  
+            goto 2  
+         endif  
+      endif  
+ 3    if(jhi-jlo.eq.1)then  
+         if(x.eq.xx(n))jlo=n-1  
+         if(x.eq.xx(1))jlo=1  
+         return  
+      endif  
+      jm=(jhi+jlo)/2  
+      if(x.ge.xx(jm).eqv.ascnd)then  
+         jlo=jm  
+      else  
+         jhi=jm  
+      endif  
+      goto 3  
+c     
+      END  
+
       
-      integer  j,np,i,nl,w
-      real*8   c(nl), cs(nl)
-
-      if(w.eq.0) then
-         do j=1,nl
-            cs(j)=c(j)
-         end do
-         
-      elseif(w.eq.3) then
-
-!       write (*,*) 'smoothing w=3'
-         do j=1,i-4
-            if(j.eq.1) then
-               cs(j)=c(j)
-            else
-               cs(j)=1/3.d0*(c(j-1)+c(j)+c(j+1))
-            end if
-         end do
-         do j=i+4,nl-1
-            if(j.eq.nl) then
-               cs(j)=c(j)
-            else
-               cs(j)=1/3.d0*(c(j-1)+c(j)+c(j+1))
-            end if
-         end do
-      elseif(w.eq.5) then
-
-!       type *,'smoothing w=5'
-         do j=3,i-4
-            if(j.eq.1) then
-               cs(j)=c(j)
-            else
-               cs(j)=1/5.d0*(c(j-2)+c(j-1)+c(j)+c(j+1)+c(j+2))
-            end if
-         end do
-         do j=i+4,nl-2
-            if(j.eq.nl) then
-               cs(j)=c(j)
-            else
-               cs(j)=1/5.d0*(c(j-2)+c(j-1)+c(j)+c(j+1)+c(j+2))
-            end if
-         end do
-      end if
-      return
-      end
-
-
-
-c*****************************************************************************
-c     suaviza
-c*****************************************************************************
-c
-      subroutine suaviza ( x, n, ismooth, y )
-c
-c     x - input and return values 
-c     y - auxiliary vector
-c     ismooth = 0  --> no smoothing is performed
-c     ismooth = 1  --> weak smoothing (5 points, centred weighted)
-c     ismooth = 2  --> normal smoothing (3 points, evenly weighted)
-c     ismooth = 3  --> strong smoothing (5 points, evenly weighted)
-
-
-c     malv  august 1991
-c*****************************************************************************
-
-      implicit none
-
-      integer   n, imax, imin, i, ismooth
-      real*8    x(n), y(n)
-c*****************************************************************************
-
-      imin=1
-      imax=n
-
-      if (ismooth.eq.0) then
-
-         return 
-
-      elseif (ismooth.eq.1) then ! 5 points, with central weighting 
-
-         do i=imin,imax 
-            if(i.eq.imin)then
-               y(i)=x(imin)
-            elseif(i.eq.imax)then
-               y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0
-            elseif(i.gt.(imin+1) .and. i.lt.(imax-1) )then
-               y(i) = ( x(i+2)/4.d0 + x(i+1)/2.d0 + 2.d0*x(i)/3.d0 + 
-     &              x(i-1)/2.d0 + x(i-2)/4.d0 )* 6.d0/13.d0
-            else
-               y(i)=(x(i+1)/2.d0+x(i)+x(i-1)/2.d0)/2.d0
-            end if
-         end do 
-         
-      elseif (ismooth.eq.2) then ! 3 points, evenly spaced
-
-         do i=imin,imax 
-            if(i.eq.imin)then
-               y(i)=x(imin)
-            elseif(i.eq.imax)then
-               y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0
-            else
-               y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0
-            end if
-         end do 
-          
-      elseif (ismooth.eq.3) then ! 5 points, evenly spaced
-
-         do i=imin,imax 
-            if(i.eq.imin)then
-               y(i) = x(imin)
-            elseif(i.eq.(imin+1) .or. i.eq.(imax-1))then
-               y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0 
-            elseif(i.eq.imax)then
-               y(i) = ( x(imax-1) + x(imax-1) + x(imax-2) ) / 3.d0
-            else
-               y(i) = ( x(i+2)+x(i+1)+x(i)+x(i-1)+x(i-2) )/5.d0
-            end if
-         end do 
-
-      else
-
-         write (*,*) ' Error in suaviza.f   Wrong ismooth value.'
-         stop
-
-      endif
-
-c rehago el cambio, para devolver x(i)
-      do i=imin,imax 
-         x(i)=y(i)
-      end do
-
-      return
-      end
-
-
-
-
-c*****************************************************************************
-c     LUdec.F (includes lubksb_dp and ludcmp_dp subroutines)
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-c Solution of linear equation without inverting matrix
-c using LU decomposition: 
-c        AA * xx = bb         AA, bb: known
-c                                 xx: to be found
-c AA and bb are not modified in this subroutine
-c                               
-c MALV , Sep 2007
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      subroutine LUdec(xx,aa,bb,m,n)
-
-      implicit none
-
-! Arguments 
-      integer,intent(in) ::     m, n
-      real*8,intent(in) ::      aa(m,m), bb(m)
-      real*8,intent(out) ::     xx(m)
-
-
-! Local variables
-      real*8      a(n,n), b(n), x(n), d
-      integer    i, j, indx(n)      
-
-
-! Subrutinas utilizadas
-!     ludcmp_dp, lubksb_dp
-
-!!!!!!!!!!!!!!! Comienza el programa !!!!!!!!!!!!!!
+c     *** hunt.f ***
+
+cccc  
+      SUBROUTINE hunt(xx,n,x,jlo)  
+c     
+c     Arguments
+      INTEGER jlo               ! O
+      INTEGER n                 ! I
+      REAL  xx(n)               ! I
+      REAL  x                   ! I
+
+c     Local variables
+      INTEGER inc,jhi,jm  
+      LOGICAL ascnd  
+c     
+cccc  
+c     
+      ascnd=xx(n).ge.xx(1)  
+      if(jlo.le.0.or.jlo.gt.n)then  
+         jlo=0  
+         jhi=n+1  
+         goto 3  
+      endif  
+      inc=1  
+      if(x.ge.xx(jlo).eqv.ascnd)then  
+ 1       jhi=jlo+inc  
+!     write (*,*) jlo
+         if(jhi.gt.n)then  
+            jhi=n+1  
+!     write (*,*) jhi-1
+         else if(x.ge.xx(jhi).eqv.ascnd)then  
+            jlo=jhi  
+            inc=inc+inc  
+!     write (*,*) jlo
+            goto 1  
+         endif  
+      else  
+         jhi=jlo  
+ 2       jlo=jhi-inc  
+!     write (*,*) jlo
+         if(jlo.lt.1)then  
+            jlo=0  
+         else if(x.lt.xx(jlo).eqv.ascnd)then  
+            jhi=jlo  
+            inc=inc+inc  
+            goto 2  
+         endif  
+      endif  
+ 3    if(jhi-jlo.eq.1)then  
+         if(x.eq.xx(n))jlo=n-1  
+         if(x.eq.xx(1))jlo=1  
+!     write (*,*) jlo
+         return  
+      endif  
+      jm=(jhi+jlo)/2  
+      if(x.ge.xx(jm).eqv.ascnd)then  
+         jlo=jm  
+      else  
+         jhi=jm  
+      endif  
+!     write (*,*) jhi-1
+      goto 3  
+c     
+      END  
+
       
-      do i=1,n
-        b(i) = bb(i+1)
-        do j=1,n
-           a(i,j) = aa(i+1,j+1)
-        enddo
-      enddo
-
-      ! Descomposicion de auxm1
-      call ludcmp_dp ( a, n, n, indx, d)
-
-      ! Sustituciones foward y backwards para hallar la solucion
-      do i=1,n
-           x(i) = b(i)
-      enddo
-      call lubksb_dp( a, n, n, indx, x )
-
-      do i=1,n
-        xx(i+1) = x(i)
-      enddo
-
-      return
-      end
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      subroutine ludcmp_dp(a,n,np,indx,d)
-
-c       jul 2011 malv+fgg
-
-      implicit none
-
-      integer,intent(in) :: n, np
-      real*8,intent(inout) :: a(np,np)
-      real*8,intent(out) :: d
-      integer,intent(out) :: indx(n)
-      
-      integer i, j, k, imax
-      real*8,parameter :: tiny=1.0d-20                                       
-      real*8 vv(n), aamax, sum, dum
-
-
-      d=1.0d0
-      do 12 i=1,n                                                             
-        aamax=0.0d0
-        do 11 j=1,n                                                           
-          if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))                         
-11      continue                                                              
-        if (aamax.eq.0.0) then
-          write(*,*) 'ludcmp_dp: singular matrix!'
-          stop
-        endif                         
-        vv(i)=1.0d0/aamax                          
-12    continue                                                                
-      do 19 j=1,n                                                             
-        if (j.gt.1) then                                                      
-          do 14 i=1,j-1                                                       
-            sum=a(i,j)                                                        
-            if (i.gt.1)then                                                   
-              do 13 k=1,i-1                                                   
-                sum=sum-a(i,k)*a(k,j)                                         
-13            continue                                                        
-              a(i,j)=sum                                                      
-            endif                                                             
-14        continue                                                            
-        endif                                                                 
-        aamax=0.0d0                                                           
-        do 16 i=j,n                                                           
-          sum=a(i,j)                                                          
-          if (j.gt.1)then                                                     
-            do 15 k=1,j-1                                                     
-              sum=sum-a(i,k)*a(k,j)                                           
-15          continue                                                          
-            a(i,j)=sum                                                        
-          endif                                                               
-          dum=vv(i)*abs(sum)                                                  
-          if (dum.ge.aamax) then                                              
-            imax=i                                                            
-            aamax=dum                                                         
-          endif                                                               
-16      continue                                                              
-        if (j.ne.imax)then                                                    
-          do 17 k=1,n                                                         
-            dum=a(imax,k)                                                     
-            a(imax,k)=a(j,k)                                                  
-            a(j,k)=dum                                                        
-17        continue                                                            
-          d=-d                                                                
-          vv(imax)=vv(j)                                                      
-        endif                                                                 
-        indx(j)=imax                                                          
-        if(j.ne.n)then                                                        
-          if(a(j,j).eq.0.0)a(j,j)=tiny                                     
-          dum=1.0d0/a(j,j)                                                  
-          do 18 i=j+1,n                                                       
-            a(i,j)=a(i,j)*dum                                                 
-18        continue                                                            
-        endif                                                                 
-19    continue                                                                
-      if(a(n,n).eq.0.0)a(n,n)=tiny                                     
-      return                                                                  
-      end                                                                     
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-      subroutine lubksb_dp(a,n,np,indx,b)                              
-
-c     jul 2011 malv+fgg
-
-      implicit none
-
-      integer,intent(in) :: n,np
-      real*8,intent(in) ::  a(np,np)
-      integer,intent(in) :: indx(n)
-      real*8,intent(out) :: b(n) 
-
-      real*8 sum
-      integer ii, ll, i, j
-
-      ii=0                                                              
-      do 12 i=1,n                                                             
-        ll=indx(i)                                                            
-        sum=b(ll)                                                             
-        b(ll)=b(i)                                                            
-        if (ii.ne.0)then                                                      
-          do 11 j=ii,i-1                                                      
-            sum=sum-a(i,j)*b(j)                                               
-11        continue                                                            
-        else if (sum.ne.0.0) then                       
-          ii=i                                                                
-        endif                                                                 
-        b(i)=sum                                                              
-12    continue                                                                
-      do 14 i=n,1,-1                                                          
-        sum=b(i)                                                              
-        if(i.lt.n)then                                                        
-          do 13 j=i+1,n                                                       
-            sum=sum-a(i,j)*b(j)                                               
-13        continue                                                            
-        endif                                                                 
-        b(i)=sum/a(i,i)                                                       
-14    continue                                                                
-      return                                                                  
-      end 
-
-
-
-
-c*****************************************************************************
-c     intersp
-c     ***********************************************************************
-      subroutine intersp(yy,zz,m,y,z,n,opt)
-c     interpolation soubroutine. input values: y(n) at z(n). 
-c     output values: yy(m) at zz(m). options: 1 -> lineal; 2 -> logarithmic
-
-c     jul 2011 malv+fgg
-c     ***********************************************************************
-
-      implicit none
-
-      integer   n,m,i,j,opt
-      real      zz(m),yy(m),z(n),y(n)
-      real      zmin,zzmin,zmax,zzmax
-
-!       write(*,*) ' interpolating'
-!       call minsp(z,n,zmin)
-      zmin=minval(z)
-!       call minsp(zz,m,zzmin)
-      zzmin=minval(zz)
-!       call maxsp(z,n,zmax)
-      zmax=maxval(z)
-!       call maxsp(zz,m,zzmax)
-      zzmax=maxval(zz)
-
-      if(zzmin.lt.zmin)then
-         write(*,*) 'from interp: new variable out of limits'
-         write(*,*) zzmin,'must be .ge. ',zmin
-         stop
-!       elseif(zzmax.gt.zmax)then
-!         write(*,*)'from interp: new variable out of limits'
-!         write(*,*)zzmax, 'must be .le. ',zmax
-!         stop
-      end if
-
-      do 1,i=1,m
-
-         do 2,j=1,n-1
-            if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3
- 2       continue
-c       in this case (zz(m).ge.z(n)) and j leaves the loop with j=n-1+1=n
-         if(opt.eq.1)then
-            yy(i)=y(n-1)+(y(n)-y(n-1))*(zz(i)-z(n-1))/(z(n)-z(n-1))
-         elseif(opt.eq.2)then
-            if(y(n).eq.0.0.or.y(n-1).eq.0.0)then
-               yy(i)=0.0
-            else
-               yy(i)=exp(log(y(n-1))+log(y(n)/y(n-1))*
-     @              (zz(i)-z(n-1))/(z(n)-z(n-1)))
-            end if
-         else
-            write(*,*)'from interp error: opt must be 1 or 2, opt= ',opt
-         end if
-         goto 1
- 3       continue
-         if(opt.eq.1)then
-            yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j))
-         elseif(opt.eq.2)then
-            if(y(j+1).eq.0.0.or.y(j).eq.0.0)then
-               yy(i)=0.0
-            else
-               yy(i)=exp(log(y(j))+log(y(j+1)/y(j))*
-     @              (zz(i)-z(j))/(z(j+1)-z(j)))
-            end if
-         else
-            write(*,*)'from interp error: opt must be 1 or 2, opt= ',opt
-         end if
- 1    continue
-
-      return
-      end
-
-
-
-c*****************************************************************************
-c     interdp
-c     ***********************************************************************
-      subroutine interdp(yy,zz,m,y,z,n,opt)
-c     interpolation soubroutine. input values: y(n) at z(n). 
-c     output values: yy(m) at zz(m). options: 1 -> lineal; 2 -> logarithmic
-c     jul 2011:  malv+fgg   Adapted to LMD-MGCM
-c     ***********************************************************************
-      implicit none
-      integer n,m,i,j,opt
-      real*8 zz(m),yy(m),z(n),y(n), zmin,zzmin,zmax,zzmax
-
-!       write (*,*) ' d interpolating '
-!       call mindp (z,n,zmin)
-      zmin=minval(z)
-!       call mindp (zz,m,zzmin)
-      zzmin=minval(zz)
-!       call maxdp (z,n,zmax)
-      zmax=maxval(z)
-!       call maxdp (zz,m,zzmax)
-      zzmax=maxval(zz)
+c     *** interdp_limits.f ***
+
+c     ***********************************************************************
+
+      subroutine interdp_limits ( yy, zz, m,   i1,i2, 
+     @     y,  z, n,   j1,j2,  opt)
+
+c     Interpolation soubroutine. 
+c     Returns values between indexes i1 & i2, donde  1 =< i1 =< i2 =< m
+c     Solo usan los indices de los inputs entre j1,j2, 1 =< j1 =< j2 =< n    
+c     Input values: y(n) , z(n)  (solo se usarann los valores entre j1,j2)
+c     zz(m) (solo se necesita entre i1,i2)
+c     Output values: yy(m) (solo se calculan entre i1,i2)
+c     Options:    opt=1 -> lineal ,,  opt=2 -> logarithmic
+c     Difference with interdp:  
+c     here interpolation proceeds between indexes i1,i2 only 
+c     if i1=1 & i2=m, both subroutines are exactly the same
+c     thus previous calls to interdp or interdp2 could be easily replaced
+
+c     JAN 98    MALV            Version for mz1d
+c     ***********************************************************************
+
+      implicit none
+
+!     Arguments 
+      integer   n,m             ! I. Dimensions
+      integer   i1, i2, j1, j2, opt ! I
+      real*8            zz(m)   ! I
+      real*8            yy(m)   ! O
+      real*8            z(n),y(n) ! I
+
+!     Local variables
+      integer   i,j
+      real*8            zmin,zzmin,zmax,zzmax
+
+c     *******************************
+
+!     write (*,*) ' d interpolating '
+!     call mindp_limits (z,n,zmin, j1,j2)
+!     call mindp_limits (zz,m,zzmin, i1,i2)
+!     call maxdp_limits (z,n,zmax, j1,j2)
+!     call maxdp_limits (zz,m,zzmax, i1,i2)
+      zmin=minval(z(j1:j2))
+      zzmin=minval(zz(i1:i2))
+      zmax=maxval(z(j1:j2))
+      zzmax=maxval(zz(i1:i2))
 
       if(zzmin.lt.zmin)then
@@ -1462 +895 @@
          write (*,*) zzmin,'must be .ge. ',zmin
          stop
-!       elseif(zzmax.gt.zmax)then
-!               write (*,*) 'from interp: new variable out of limits'
-!               write (*,*) zzmax, 'must be .le. ',zmax
-!               stop
-      end if
-
-      do 1,i=1,m
-
-         do 2,j=1,n-1
-            if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3
- 2       continue
-c       in this case (zz(m).eq.z(n)) and j leaves the loop with j=n-1+1=n
-         if(opt.eq.1)then
-            yy(i)=y(n-1)+(y(n)-y(n-1))*(zz(i)-z(n-1))/(z(n)-z(n-1))
-         elseif(opt.eq.2)then
-            if(y(n).eq.0.0d0.or.y(n-1).eq.0.0d0)then
-               yy(i)=0.0d0
-            else
-               yy(i)=dexp(dlog(y(n-1))+dlog(y(n)/y(n-1))*
-     @              (zz(i)-z(n-1))/(z(n)-z(n-1)))
-            end if
-         else
-            write (*,*) 
-     @           ' from d interp error: opt must be 1 or 2, opt= ',opt
-         end if
-         goto 1
- 3       continue
-         if(opt.eq.1)then
-            yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j))
-!       write (*,*) ' '
-!       write (*,*) ' z(j),z(j+1) =', z(j),z(j+1)
-!       write (*,*) ' t(j),t(j+1) =', y(j),y(j+1)
-!       write (*,*) ' zz, tt =  ', zz(i), yy(i)
-         elseif(opt.eq.2)then
-            if(y(j+1).eq.0.0d0.or.y(j).eq.0.0d0)then
-               yy(i)=0.0d0
-            else
-               yy(i)=dexp(dlog(y(j))+dlog(y(j+1)/y(j))*
-     @              (zz(i)-z(j))/(z(j+1)-z(j)))
-            end if
-         else
-            write (*,*) ' from interp error: opt must be 1 or 2, opt= ',
-     @           opt
-         end if
- 1    continue
-      return
-      end
-
-
-c*****************************************************************************
-c     interdp_limits.F
-c     ***********************************************************************
-
-      subroutine interdp_limits ( yy,zz,m, i1,i2, y,z,n, j1,j2, opt)
-
-c     Interpolation soubroutine. 
-c     Returns values between indexes i1 & i2, donde  1 =< i1 =< i2 =< m
-c     Solo usan los indices de los inputs entre j1,j2, 1 =< j1 =< j2 =< n    
-c     Input values: y(n) , z(n)  (solo se usan los valores entre j1,j2)
-c                     zz(m) (solo se necesita entre i1,i2)
-c     Output values: yy(m) (solo se calculan entre i1,i2)
-c     Options:    opt=1 -> lineal ,,  opt=2 -> logarithmic
-c     Difference with interdp:  
-c          here interpolation proceeds between indexes i1,i2 only 
-c          if i1=1 & i2=m, both subroutines are exactly the same
-c          thus previous calls to interdp or interdp2 could be easily replaced
-
-c     JAN 98    MALV            Version for mz1d
-c     jul 2011 malv+fgg       Adapted to LMD-MGCM
-c     ***********************************************************************
-
-      implicit none
-
-! Arguments 
-      integer   n,m             ! I. Dimensions
-      integer   i1, i2, j1, j2, opt ! I
-      real*8            zz(m),yy(m) ! O
-      real*8            z(n),y(n) ! I
-
-! Local variables
-      integer   i,j
-      real*8            zmin,zzmin,zmax,zzmax
-
-c     *******************************
-
-!       type *, ' d interpolating '
-!       call mindp_limits (z,n,zmin, j1,j2)
-      zmin=minval(z(j1:j2))
-!       call mindp_limits (zz,m,zzmin, i1,i2)
-      zzmin=minval(zz(i1:i2))
-!       call maxdp_limits (z,n,zmax, j1,j2)
-      zmax=maxval(z(j1:j2))
-!       call maxdp_limits (zz,m,zzmax, i1,i2)
-      zzmax=maxval(zz(i1:i2))
-
-      if(zzmin.lt.zmin)then
-         write (*,*) 'from d interp: new variable out of limits'
-         write (*,*) zzmin,'must be .ge. ',zmin
-         stop
-!       elseif(zzmax.gt.zmax)then
-!               type *,'from interp: new variable out of limits'
-!               type *,zzmax, 'must be .le. ',zmax
-!               stop
+!     elseif(zzmax.gt.zmax)then
+!     type *,'from interp: new variable out of limits'
+!     type *,zzmax, 'must be .le. ',zmax
+!     stop
       end if
 
@@ -1572 +906 @@
             if(zz(i).ge.z(j).and.zz(i).lt.z(j+1)) goto 3
  2       continue
-c       in this case (zz(i2).eq.z(j2)) and j leaves the loop with j=j2-1+1=j2
+c     in this case (zz(i2).eq.z(j2)) and j leaves the loop with j=j2-1+1=j2
          if(opt.eq.1)then
-            yy(i)=y(j2-1)+(y(j2)-y(j2-1))*
-     $           (zz(i)-z(j2-1))/(z(j2)-z(j2-1))
+            yy(i)=y(j2-1)+(y(j2)-y(j2-1))*(zz(i)-z(j2-1))/
+     $           (z(j2)-z(j2-1))
          elseif(opt.eq.2)then
             if(y(j2).eq.0.0d0.or.y(j2-1).eq.0.0d0)then
@@ -1590 +924 @@
          if(opt.eq.1)then
             yy(i)=y(j)+(y(j+1)-y(j))*(zz(i)-z(j))/(z(j+1)-z(j))
-!       type *, ' '
-!       type *, ' z(j),z(j+1) =', z(j),z(j+1)
-!       type *, ' t(j),t(j+1) =', y(j),y(j+1)
-!       type *, ' zz, tt =  ', zz(i), yy(i)
+!     type *, ' '
+!     type *, ' z(j),z(j+1) =', z(j),z(j+1)
+!     type *, ' t(j),t(j+1) =', y(j),y(j+1)
+!     type *, ' zz, tt =  ', zz(i), yy(i)
          elseif(opt.eq.2)then
             if(y(j+1).eq.0.0d0.or.y(j).eq.0.0d0)then
@@ -1610 +944 @@
 
 
-
-c*****************************************************************************
-c     Subroutines previously included in tcrco2_subrut.F
+c     *** interhunt2veces.f ***
+
+c     ***********************************************************************
+
+      subroutine interhunt2veces ( y1,y2,  zz,m, 
+     @     x1,x2,  z,n,  opt)
+
+c     interpolation soubroutine basada en Numerical Recipes HUNT.FOR
+c     input values: y(n) at z(n) 
+c     output values: yy(m) at zz(m)
+c     options: 1 -> lineal
+c     2 -> logarithmic
+c     ***********************************************************************
+
+      implicit none
+
+!     Arguments
+      integer   n,m,opt         ! I
+      real      zz(m),z(n)      ! I
+      real    y1(m),y2(m)       ! O
+      real    x1(n),x2(n)       ! I
+
+
+!     Local variables 
+      integer i, j 
+      real    factor 
+      real    zaux
+
+!!!!  
+
+      j = 1                     ! initial first guess (=n/2 is anothr pssblty)
+
+      do 1,i=1,m                ! 
+
+                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
+         zaux = zz(i)
+         if (abs(zaux-z(1)).le.0.01) then 
+            zaux=z(1)
+         elseif (abs(zaux-z(n)).le.0.01) then
+            zaux=z(n)
+         endif
+         call hunt ( z,n, zaux, j )
+         if ( j.eq.0 .or. j.eq.n ) then 
+            write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
+            write (*,*) ' HUNT/ location in new grid:', zz(i)
+            stop ' interhunt2/ Interpolat error. zz out of limits.'
+         endif
+
+                                ! Perform interpolation
+         factor = (zz(i)-z(j))/(z(j+1)-z(j))
+         if (opt.eq.1) then
+            y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
+            y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
+         else
+            y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
+            y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
+         end if
+
+ 1    continue
+
+      return
+      end
+
+
+c     *** interhunt5veces.f ***
+
+c     ***********************************************************************
+
+      subroutine interhunt5veces ( y1,y2,y3,y4,y5,  zz,m, 
+     @     x1,x2,x3,x4,x5,  z,n,  opt)
+
+c     interpolation soubroutine basada en Numerical Recipes HUNT.FOR
+c     input values: y(n) at z(n) 
+c     output values: yy(m) at zz(m)
+c     options: 1 -> lineal
+c     2 -> logarithmic
+c     ***********************************************************************
+
+      implicit none
+!     Arguments
+      integer   n,m,opt         ! I
+      real      zz(m),z(n)      ! I
+      real    y1(m),y2(m),y3(m),y4(m),y5(m) ! O
+      real    x1(n),x2(n),x3(n),x4(n),x5(n) ! I
+
+
+!     Local variables 
+      integer i, j 
+      real    factor 
+      real    zaux
+
+!!!!  
+
+      j = 1                     ! initial first guess (=n/2 is anothr pssblty)
+
+      do 1,i=1,m                ! 
+
+                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
+         zaux = zz(i)
+         if (abs(zaux-z(1)).le.0.01) then 
+            zaux=z(1)
+         elseif (abs(zaux-z(n)).le.0.01) then
+            zaux=z(n)
+         endif
+         call hunt ( z,n, zaux, j )
+         if ( j.eq.0 .or. j.eq.n ) then 
+            write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
+            write (*,*) ' HUNT/ location in new grid:', zz(i)
+            stop ' interhunt5/ Interpolat error. zz out of limits.'
+         endif
+
+                                ! Perform interpolation
+         factor = (zz(i)-z(j))/(z(j+1)-z(j))
+         if (opt.eq.1) then
+            y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
+            y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
+            y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor
+            y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor
+            y5(i) = x5(j) + (x5(j+1)-x5(j)) * factor
+         else
+            y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
+            y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
+            y3(i) = exp( log(x3(j)) + log(x3(j+1)/x3(j)) * factor )
+            y4(i) = exp( log(x4(j)) + log(x4(j+1)/x4(j)) * factor )
+            y5(i) = exp( log(x5(j)) + log(x5(j+1)/x5(j)) * factor )
+         end if
+
+ 1    continue
+
+      return
+      end
+
+
+
+c     *** interhuntdp3veces.f ***
+
+c     ***********************************************************************
+
+      subroutine interhuntdp3veces ( y1,y2,y3, zz,m, 
+     @     x1,x2,x3,  z,n,  opt)
+
+c     interpolation soubroutine basada en Numerical Recipes HUNT.FOR
+c     input values: x(n) at z(n) 
+c     output values: y(m) at zz(m)
+c     options: opt = 1 -> lineal
+c     opt=/=1 -> logarithmic
+c     ***********************************************************************
+!     Arguments
+      integer   n,m,opt         ! I
+      real*8    zz(m),z(n)      ! I
+      real*8    y1(m),y2(m),y3(m) ! O
+      real*8    x1(n),x2(n),x3(n) ! I
+
+
+!     Local variables 
+      integer i, j 
+      real*8    factor 
+      real*8    zaux
+
+!!!!  
+
+      j = 1                     ! initial first guess (=n/2 is anothr pssblty)
+
+      do 1,i=1,m                ! 
+
+                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
+         zaux = zz(i)
+         if (abs(zaux-z(1)).le.0.01d0) then 
+            zaux=z(1)
+         elseif (abs(zaux-z(n)).le.0.01d0) then
+            zaux=z(n)
+         endif
+         call huntdp ( z,n, zaux, j )
+         if ( j.eq.0 .or. j.eq.n ) then 
+            write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
+            write (*,*) ' HUNT/ location in new grid:', zz(i)
+            stop ' INTERHUNTDP3/ Interpolat error. zz out of limits.'
+         endif
+
+                                ! Perform interpolation
+         factor = (zz(i)-z(j))/(z(j+1)-z(j))
+         if (opt.eq.1) then
+            y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
+            y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
+            y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor
+         else
+            y1(i) = dexp( dlog(x1(j)) + dlog(x1(j+1)/x1(j)) * factor )
+            y2(i) = dexp( dlog(x2(j)) + dlog(x2(j+1)/x2(j)) * factor )
+            y3(i) = dexp( dlog(x3(j)) + dlog(x3(j+1)/x3(j)) * factor )
+         end if
+
+ 1    continue
+
+      return
+      end
+
+
+c     *** interhuntdp4veces.f ***
+
+c     ***********************************************************************
+
+      subroutine interhuntdp4veces ( y1,y2,y3,y4, zz,m, 
+     @     x1,x2,x3,x4,  z,n,  opt)
+
+c     interpolation soubroutine basada en Numerical Recipes HUNT.FOR
+c     input values: x1(n),x2(n),x3(n),x4(n) at z(n) 
+c     output values: y1(m),y2(m),y3(m),y4(m) at zz(m)
+c     options: 1 -> lineal
+c     2 -> logarithmic
+c     ***********************************************************************
+
+      implicit none
+
+!     Arguments
+      integer   n,m,opt         ! I
+      real*8    zz(m),z(n)      ! I
+      real*8    y1(m),y2(m),y3(m),y4(m) ! O
+      real*8    x1(n),x2(n),x3(n),x4(n) ! I
+
+
+!     Local variables 
+      integer i, j 
+      real*8    factor 
+      real*8    zaux
+
+!!!!  
+
+      j = 1                     ! initial first guess (=n/2 is anothr pssblty)
+
+      do 1,i=1,m                ! 
+
+                                ! Caza del indice j donde ocurre que zz(i) esta entre [z(j),z(j+1)]
+         zaux = zz(i)
+         if (abs(zaux-z(1)).le.0.01d0) then 
+            zaux=z(1)
+         elseif (abs(zaux-z(n)).le.0.01d0) then
+            zaux=z(n)
+         endif
+         call huntdp ( z,n, zaux, j )
+         if ( j.eq.0 .or. j.eq.n ) then 
+            write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
+            write (*,*) ' HUNT/ location in new grid:', zz(i)
+            stop ' INTERHUNTDP4/ Interpolat error. zz out of limits.'
+         endif
+
+                                ! Perform interpolation
+         factor = (zz(i)-z(j))/(z(j+1)-z(j))
+         if (opt.eq.1) then
+            y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
+            y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor
+            y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor
+            y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor
+         else
+            y1(i) = dexp( dlog(x1(j)) + dlog(x1(j+1)/x1(j)) * factor )
+            y2(i) = dexp( dlog(x2(j)) + dlog(x2(j+1)/x2(j)) * factor )
+            y3(i) = dexp( dlog(x3(j)) + dlog(x3(j+1)/x3(j)) * factor )
+            y4(i) = dexp( dlog(x4(j)) + dlog(x4(j+1)/x4(j)) * factor )
+         end if
+
+ 1    continue
+
+      return
+      end
+
+
+c     *** interhuntdp.f ***
+
+c     ***********************************************************************
+
+      subroutine interhuntdp ( y1, zz,m, 
+     @     x1,  z,n,  opt)
+
+c     interpolation soubroutine basada en Numerical Recipes HUNT.FOR
+c     input values: x1(n) at z(n) 
+c     output values: y1(m) at zz(m)
+c     options: 1 -> lineal
+c     2 -> logarithmic
+c     ***********************************************************************
+
+      implicit none 
+
+!     Arguments
+      integer   n,m,opt         ! I
+      real*8    zz(m),z(n)      ! I
+      real*8    y1(m)           ! O
+      real*8    x1(n)           ! I
+
+
+!     Local variables 
+      integer i, j 
+      real*8    factor 
+      real*8    zaux
+
+!!!!  
+
+      j = 1                     ! initial first guess (=n/2 is anothr pssblty)
+
+      do 1,i=1,m                ! 
+
+                                ! Caza del indice j donde ocurre que zz(i) esta entre [z(j),z(j+1)]
+         zaux = zz(i)
+         if (abs(zaux-z(1)).le.0.01d0) then 
+            zaux=z(1)
+         elseif (abs(zaux-z(n)).le.0.01d0) then
+            zaux=z(n)
+         endif
+         call huntdp ( z,n, zaux, j )
+         if ( j.eq.0 .or. j.eq.n ) then 
+            write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
+            write (*,*) ' HUNT/ location in new grid:', zz(i)
+            stop ' INTERHUNT/ Interpolat error. zz out of limits.'
+         endif
+
+                                ! Perform interpolation
+         factor = (zz(i)-z(j))/(z(j+1)-z(j))
+         if (opt.eq.1) then
+            y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
+         else
+            y1(i) = dexp( dlog(x1(j)) + dlog(x1(j+1)/x1(j)) * factor )
+         end if
+
+ 1    continue
+
+      return
+      end
+
+
+c     *** interhunt.f ***
+
 c***********************************************************************
-c     tcrco2_subrut.f                              
-c                                               
-c     jan 98    malv    version for mz1d. copied from solar10/mz4sub.f         
-c     jul 2011 malv+fgg   adapted to LMD-MGCM
+
+      subroutine interhunt ( y1, zz,m, 
+     @     x1,  z,n,  opt)
+
+c     interpolation soubroutine basada en Numerical Recipes HUNT.FOR
+c     input values: x1(n) at z(n) 
+c     output values: y1(m) at zz(m)
+c     options: 1 -> lineal
+c     2 -> logarithmic
 c***********************************************************************
-                                                
-************************************************************************
-                                                
-      subroutine dinterconnection ( v, vt )          
-                                                
-*  input: vib. temp. from che*.for programs, vt(nl)         
-*  output: test vibrational temp. for other che*.for, v(nl) 
-! iconex_smooth=1  ==>  with smoothing          
-! iconex_smooth=0  ==>  without smoothing       
-! iconex_tk=40  ==>  with forced lte up to 40 km            
-! iconex_tk=20  ==>  with forced lte up to 20 km            
-************************************************************************
-                                                
-      implicit none                           
-      include 'nlte_paramdef.h'
-      include 'nlte_commons.h'
-                                                
-c argumentos                                    
-      real*8 vt(nl), v(nl)                           
-                                                
-c local variables                               
-      integer   i                                     
-                                                
-c   *************                               
-                                                
-      do i=1,nl                                      
-         v(i) = vt(i)                                 
-      end do                                         
-                                                
-! lo siguiente se utilizaba en solar10, pero es mejor introducirlo en   
-! la driver. por ahora no lo uso todavia.       
-!       call fluctua(v,iconex_fluctua)                
-!       call smooth_nl(v,iconex_smooth,nl)               
-!       call forzar_tk(v,iconex_tk)                   
-                                                
-      return                                         
-      end                  
-                                                
+
+      implicit none
+
+!     Arguments
+      integer   n,m,opt         ! I
+      real      zz(m),z(n)      ! I
+      real    y1(m)             ! O
+      real    x1(n)             ! I
+
+
+!     Local variables 
+      integer i, j 
+      real    factor 
+      real    zaux
+
+!!!!  
+
+      j = 1                     ! initial first guess (=n/2 is anothr pssblty)
+
+      do 1,i=1,m                ! 
+
+                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
+         zaux = zz(i)
+         if (abs(zaux-z(1)).le.0.01) then 
+            zaux=z(1)
+         elseif (abs(zaux-z(n)).le.0.01) then
+            zaux=z(n)
+         endif
+         call hunt ( z,n, zaux, j )
+         if ( j.eq.0 .or. j.eq.n ) then 
+            write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
+            write (*,*) ' HUNT/ location in new grid:', zz(i)
+            stop ' interhunt/ Interpolat error. z out of limits.'
+         endif
+
+                                ! Perform interpolation
+         factor = (zz(i)-z(j))/(z(j+1)-z(j))
+         if (opt.eq.1) then
+            y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
+         else
+            y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
+         end if
+
+
+ 1    continue
+
+      return
+      end
+
+
+c     *** interhuntlimits2veces.f ***
+
 c***********************************************************************
-      function planckdp(tp,xnu)                      
-c     returns the black body function at wavenumber xnu and temperature t.  
+
+      subroutine interhuntlimits2veces 
+     @     ( y1,y2, zz,m, limite1,limite2, 
+     @     x1,x2,  z,n, opt)
+
+c     Interpolation soubroutine basada en Numerical Recipes HUNT.FOR
+c     Input values:  x1,x2(n) at z(n) 
+c     Output values: 
+c     y1,y2(m) at zz(m)   pero solo entre los indices de zz
+c     siguientes: [limite1,limite2]
+c     Options: 1 -> linear in z and linear in x
+c     2 -> linear in z and logarithmic in x
+c     3 -> logarithmic in z and linear in x
+c     4 -> logarithmic in z and logaritmic in x
+c     
+c     NOTAS: Esta subrutina extiende y generaliza la usual  
+c     "interhunt5veces" en 2 direcciones: 
+c     - la condicion en los limites es que zz(limite1:limite2)
+c     esté dentro de los limites de z (pero quizas no todo zz)
+c     - se han añadido 3 opciones mas al caso de interpolacion
+c     logaritmica, ahora se hace en log de z, de x o de ambos.
+c     Notese que esta subrutina engloba a la interhunt5veces 
+c     ( esta es reproducible haciendo  limite1=1  y  limite2=m 
+c     y usando una de las 2 primeras opciones  opt=1,2 )
 c***********************************************************************
-                                                
-      implicit none                                  
-
-      include 'nlte_paramdef.h'
-      include 'nlte_commons.h'
-
-!        common/datis/ pi, vlight, ee, hplanck, gamma, ab,
-!     @       n_avog, GG, R0, cte_sb, kboltzman,  raddeg
-!        real*8  pi, vlight, ee, hplanck, gamma, ab,
-!     @       n_avog, GG, R0, cte_sb, kboltzman,  raddeg
-
-      real*8 planckdp                                
-      real*8 xnu                                     
-      real tp                                        
-                                                
-      planckdp = gamma*xnu**3.0 / exp( ee*xnu/dble(tp) )             
-      !erg cm-2.sr-1/cm-1.                           
-                                                
-      return                                         
-      end                                            
-
-c     ****************************************************************
-      function bandid (ib)                           
-c     returns the 2 character code of the band            
-c     ****************************************************************
-      implicit none                           
-        
-      integer ib                              
-      character*2 bandid                      
-                                                
- 132  format(i2)                              
-!     encode (2,132,bandid) ib                
-      write ( bandid, 132) ib                
-                                                
-      if ( ib .eq. 1 ) bandid = '01'          
-      if ( ib .eq. 2 ) bandid = '02'          
-      if ( ib .eq. 3 ) bandid = '03'          
-      if ( ib .eq. 4 ) bandid = '04'          
-      if ( ib .eq. 5 ) bandid = '05'          
-      if ( ib .eq. 6 ) bandid = '06'          
-      if ( ib .eq. 7 ) bandid = '07'          
-      if ( ib .eq. 8 ) bandid = '08'          
-      if ( ib .eq. 9 ) bandid = '09'          
-      if ( ib .eq. 0 ) bandid = '00'          
-                                                
-c end                                           
-      return                                  
+
+      implicit none
+
+!     Arguments
+      integer   n,m,opt, limite1,limite2 ! I
+      real      zz(m),z(n)      ! I
+      real    y1(m),y2(m)       ! O
+      real    x1(n),x2(n)       ! I
+
+
+!     Local variables 
+      integer i, j 
+      real    factor 
+      real    zaux
+
+!!!!  
+
+      j = 1                     ! initial first guess (=n/2 is anothr pssblty)
+
+      do 1,i=limite1,limite2              
+
+                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
+         zaux = zz(i)
+         if (abs(zaux-z(1)).le.0.01) then 
+            zaux=z(1)
+         elseif (abs(zaux-z(n)).le.0.01) then
+            zaux=z(n)
+         endif
+         call hunt ( z,n, zaux, j )
+         if ( j.eq.0 .or. j.eq.n ) then 
+            write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
+            write (*,*) ' HUNT/ location in new grid:', zz(i)
+            stop ' interhuntlimits/ Interpolat error. z out of limits.'
+         endif
+
+                                ! Perform interpolation
+         if (opt.eq.1) then
+            factor = (zz(i)-z(j))/(z(j+1)-z(j))
+            y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
+            y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor 
+
+         elseif (opt.eq.2) then 
+            factor = (zz(i)-z(j))/(z(j+1)-z(j))
+            y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
+            y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
+         elseif (opt.eq.3) then 
+            factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
+            y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
+            y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor 
+         elseif (opt.eq.4) then 
+            factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
+            y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
+            y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
+         end if
+
+
+ 1    continue
+
+      return
+      end
+
+
+c     *** interhuntlimits5veces.f ***
+
+c***********************************************************************
+
+      subroutine interhuntlimits5veces 
+     @     ( y1,y2,y3,y4,y5, zz,m, limite1,limite2, 
+     @     x1,x2,x3,x4,x5,  z,n, opt)
+
+c     Interpolation soubroutine basada en Numerical Recipes HUNT.FOR
+c     Input values:  x1,x2,..,x5(n) at z(n) 
+c     Output values: 
+c     y1,y2,...,y5(m) at zz(m)   pero solo entre los indices de zz
+c     siguientes: [limite1,limite2]
+c     Options: 1 -> linear in z and linear in x
+c     2 -> linear in z and logarithmic in x
+c     3 -> logarithmic in z and linear in x
+c     4 -> logarithmic in z and logaritmic in x
+c     
+c     NOTAS: Esta subrutina extiende y generaliza la usual  
+c     "interhunt5veces" en 2 direcciones: 
+c     - la condicion en los limites es que zz(limite1:limite2)
+c     esté dentro de los limites de z (pero quizas no todo zz)
+c     - se han añadido 3 opciones mas al caso de interpolacion
+c     logaritmica, ahora se hace en log de z, de x o de ambos.
+c     Notese que esta subrutina engloba a la interhunt5veces 
+c     ( esta es reproducible haciendo  limite1=1  y  limite2=m 
+c     y usando una de las 2 primeras opciones  opt=1,2 )
+c***********************************************************************
+
+      implicit none
+
+!     Arguments
+      integer   n,m,opt, limite1,limite2 ! I
+      real      zz(m),z(n)      ! I
+      real    y1(m),y2(m),y3(m),y4(m),y5(m) ! O
+      real    x1(n),x2(n),x3(n),x4(n),x5(n) ! I
+
+
+!     Local variables 
+      integer i, j 
+      real    factor 
+      real    zaux
+
+!!!!  
+
+      j = 1                     ! initial first guess (=n/2 is anothr pssblty)
+
+      do 1,i=limite1,limite2              
+
+                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
+         zaux = zz(i)
+         if (abs(zaux-z(1)).le.0.01) then 
+            zaux=z(1)
+         elseif (abs(zaux-z(n)).le.0.01) then
+            zaux=z(n)
+         endif
+         call hunt ( z,n, zaux, j )
+         if ( j.eq.0 .or. j.eq.n ) then 
+            write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
+            write (*,*) ' HUNT/ location in new grid:', zz(i)
+            stop ' interhuntlimits/ Interpolat error. z out of limits.'
+         endif
+
+                                ! Perform interpolation
+         if (opt.eq.1) then
+            factor = (zz(i)-z(j))/(z(j+1)-z(j))
+            y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
+            y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor 
+            y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor 
+            y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor 
+            y5(i) = x5(j) + (x5(j+1)-x5(j)) * factor 
+         elseif (opt.eq.2) then 
+            factor = (zz(i)-z(j))/(z(j+1)-z(j))
+            y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
+            y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
+            y3(i) = exp( log(x3(j)) + log(x3(j+1)/x3(j)) * factor )
+            y4(i) = exp( log(x4(j)) + log(x4(j+1)/x4(j)) * factor )
+            y5(i) = exp( log(x5(j)) + log(x5(j+1)/x5(j)) * factor )
+         elseif (opt.eq.3) then 
+            factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
+            y1(i) = x1(j) + (x1(j+1)-x1(j)) * factor 
+            y2(i) = x2(j) + (x2(j+1)-x2(j)) * factor 
+            y3(i) = x3(j) + (x3(j+1)-x3(j)) * factor 
+            y4(i) = x4(j) + (x4(j+1)-x4(j)) * factor 
+            y5(i) = x5(j) + (x5(j+1)-x5(j)) * factor 
+         elseif (opt.eq.4) then 
+            factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
+            y1(i) = exp( log(x1(j)) + log(x1(j+1)/x1(j)) * factor )
+            y2(i) = exp( log(x2(j)) + log(x2(j+1)/x2(j)) * factor )
+            y3(i) = exp( log(x3(j)) + log(x3(j+1)/x3(j)) * factor )
+            y4(i) = exp( log(x4(j)) + log(x4(j+1)/x4(j)) * factor )
+            y5(i) = exp( log(x5(j)) + log(x5(j+1)/x5(j)) * factor )
+         end if
+
+
+ 1    continue
+
+      return
+      end
+
+
+
+c     *** interhuntlimits.f ***
+
+c***********************************************************************
+
+      subroutine interhuntlimits ( y, zz,m, limite1,limite2, 
+     @     x,  z,n, opt)
+
+c     Interpolation soubroutine basada en Numerical Recipes HUNT.FOR
+c     Input values:  x(n) at z(n) 
+c     Output values: y(m) at zz(m)   pero solo entre los indices de zz
+c     siguientes: [limite1,limite2]
+c     Options: 1 -> linear in z and linear in x
+c     2 -> linear in z and logarithmic in x
+c     3 -> logarithmic in z and linear in x
+c     4 -> logarithmic in z and logaritmic in x
+c     
+c     NOTAS: Esta subrutina extiende y generaliza la usual  "interhunt"
+c     en 2 direcciones: 
+c     - la condicion en los limites es que zz(limite1:limite2)
+c     esté dentro de los limites de z (pero quizas no todo zz)
+c     - se han añadido 3 opciones mas al caso de interpolacion
+c     logaritmica, ahora se hace en log de z, de x o de ambos.
+c     Notese que esta subrutina engloba a la usual interhunt 
+c     ( esta es reproducible haciendo  limite1=1  y  limite2=m 
+c     y usando una de las 2 primeras opciones  opt=1,2 )
+c***********************************************************************
+
+      implicit none
+
+!     Arguments
+      integer   n,m,opt, limite1,limite2 ! I
+      real      zz(m),z(n)      ! I
+      real    y(m)              ! O
+      real    x(n)              ! I
+
+
+!     Local variables 
+      integer i, j 
+      real    factor 
+      real    zaux
+
+!!!!  
+
+      j = 1                     ! initial first guess (=n/2 is anothr pssblty)
+
+      do 1,i=limite1,limite2              
+
+                                ! Busca indice j donde ocurre q zz(i) esta entre [z(j),z(j+1)]
+         zaux = zz(i)
+         if (abs(zaux-z(1)).le.0.01) then 
+            zaux=z(1)
+         elseif (abs(zaux-z(n)).le.0.01) then
+            zaux=z(n)
+         endif
+         call hunt ( z,n, zaux, j )
+         if ( j.eq.0 .or. j.eq.n ) then 
+            write (*,*) ' HUNT/ Limits input grid:', z(1),z(n)
+            write (*,*) ' HUNT/ location in new grid:', zz(i)
+            stop ' interhuntlimits/ Interpolat error. z out of limits.'
+         endif
+
+                                ! Perform interpolation
+         if (opt.eq.1) then
+            factor = (zz(i)-z(j))/(z(j+1)-z(j))
+            y(i) = x(j) + (x(j+1)-x(j)) * factor 
+         elseif (opt.eq.2) then 
+            factor = (zz(i)-z(j))/(z(j+1)-z(j))
+            y(i) = exp( log(x(j)) + log(x(j+1)/x(j)) * factor )
+         elseif (opt.eq.3) then 
+            factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
+            y(i) = x(j) + (x(j+1)-x(j)) * factor 
+         elseif (opt.eq.4) then 
+            factor = (log(zz(i))-log(z(j)))/(log(z(j+1))-log(z(j)))
+            y(i) = exp( log(x(j)) + log(x(j+1)/x(j)) * factor )
+         end if
+
+
+ 1    continue
+
+      return
+      end
+
+
+c     *** lubksb_dp.f ***
+
+      subroutine lubksb_dp(a,n,np,indx,b)     
+
+      implicit none
+
+      integer,intent(in) :: n,np
+      real*8,intent(in) :: a(np,np) 
+      integer,intent(in) :: indx(n)
+      real*8,intent(out) :: b(n) 
+
+      real*8 sum
+      integer ii, ll, i, j
+
+      ii=0            
+      do 12 i=1,n           
+         ll=indx(i)          
+         sum=b(ll)           
+         b(ll)=b(i)          
+         if (ii.ne.0)then    
+            do 11 j=ii,i-1        
+               sum=sum-a(i,j)*b(j)     
+ 11         continue                  
+         else if (sum.ne.0.0) then  
+            ii=i                      
+         endif                       
+         b(i)=sum                    
+ 12   continue                      
+      do 14 i=n,1,-1                
+         sum=b(i)                    
+         if(i.lt.n)then              
+            do 13 j=i+1,n             
+               sum=sum-a(i,j)*b(j)     
+ 13         continue                  
+         endif                       
+         b(i)=sum/a(i,i)             
+ 14   continue                      
+      return   
+      end      
+
+
+c     *** ludcmp_dp.f ***
+
+      subroutine ludcmp_dp(a,n,np,indx,d)
+
+      implicit none
+
+      integer,intent(in) :: n, np
+      real*8,intent(inout) :: a(np,np)
+      real*8,intent(out) :: d
+      integer,intent(out) :: indx(n)
+      
+      integer nmax, i, j, k, imax
+      real*8 tiny
+      parameter (nmax=100,tiny=1.0d-20)   
+      real*8 vv(nmax), aamax, sum, dum
+
+
+      d=1.0d0
+      do 12 i=1,n             
+         aamax=0.0d0
+         do 11 j=1,n           
+            if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))             
+ 11      continue              
+         if (aamax.eq.0.0) then
+            write(*,*) 'ludcmp_dp: singular matrix!'
+            stop
+         endif            
+         vv(i)=1.0d0/aamax  
+ 12   continue                
+      do 19 j=1,n             
+         if (j.gt.1) then      
+            do 14 i=1,j-1       
+               sum=a(i,j)        
+               if (i.gt.1)then               
+                  do 13 k=1,i-1               
+                     sum=sum-a(i,k)*a(k,j)     
+ 13               continue    
+                  a(i,j)=sum      
+               endif             
+ 14         continue            
+         endif                 
+         aamax=0.0d0           
+         do 16 i=j,n           
+            sum=a(i,j)          
+            if (j.gt.1)then     
+               do 15 k=1,j-1                     
+                  sum=sum-a(i,k)*a(k,j)                     
+ 15            continue              
+               a(i,j)=sum            
+            endif                   
+            dum=vv(i)*abs(sum)      
+            if (dum.ge.aamax) then  
+               imax=i                
+               aamax=dum             
+            endif                   
+ 16      continue                  
+         if (j.ne.imax)then        
+            do 17 k=1,n             
+               dum=a(imax,k)         
+               a(imax,k)=a(j,k)      
+               a(j,k)=dum            
+ 17         continue                
+            d=-d                    
+            vv(imax)=vv(j)          
+         endif                     
+         indx(j)=imax              
+         if(j.ne.n)then            
+            if(a(j,j).eq.0.0)a(j,j)=tiny               
+            dum=1.0d0/a(j,j)      
+            do 18 i=j+1,n           
+               a(i,j)=a(i,j)*dum     
+ 18         continue                
+         endif                     
+ 19   continue                    
+      if(a(n,n).eq.0.0)a(n,n)=tiny               
+      return                      
       end   
 
 
-
-c*****************************************************************************
-c     Subroutines previously included in mat_oper.F
-c*****************************************************************************
-c set of subroutines for the cz*.for programs:
-!     subroutine unit(a,n)
-!     subroutine diago(a,v,n)             diagonal matrix with v
-!     subroutine invdiag(a,b,n)           inverse of diagonal matrix 
-!     subroutine sypvvv(a,b,c,d,n)        suma y prod de 3 vectores, muy comun
-!     subroutine sypvmv(v,w,b,u,n)        suma y prod de 3 vectores, muy comun
-!     subroutine mulmvv(w,b,u,v,n)        prod matriz vector vector
-!     subroutine muymvv(w,b,u,v,n)        prod matriz (inv.vector) vector
-!     subroutine samem (a,m,n)
-!     subroutine mulmv(a,b,c,n)
-!     subroutine mulmm(a,b,c,n)
-!     subroutine resmm(a,b,c,n)
-!     subroutine mulvv(a,b,c,n)
-!     subroutine sumvv(a,b,c,n)
-!     subroutine zerom(a,n)
-!     subroutine zero4m(a,b,c,d,n)
-!     subroutine zero3m(a,b,c,n)
-!     subroutine zero2m(a,b,n)
-!     subroutine zerov(a,n)
-!     subroutine zero4v(a,b,c,d,n)
-!     subroutine zero3v(a,b,c,n)
-!     subroutine zero2v(a,b,n)
-
-!
-!
-!   May-05 Sustituimos todos los zerojt de cristina por las subrutinas
-!          genericas zerov***
-!
+c     *** LUdec.f ***
+
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c     
+c     Solution of linear equation without inverting matrix
+c     using LU decomposition: 
+c     AA * xx = bb         AA, bb: known
+c     xx: to be found
+c     AA and bb are not modified in this subroutine
+c     
+c     MALV , Sep 2007
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      subroutine LUdec(xx,aa,bb,m,n)
+
+      implicit none
+
+!     Arguments 
+      integer,intent(in) ::     m, n
+      real*8,intent(in) ::      aa(m,m), bb(m)
+      real*8,intent(out) ::     xx(m)
+
+
+!     Local variables
+      real*8      a(n,n), b(n), x(n), d
+      integer    i, j, indx(n)      
+
+
+!     Subrutinas utilizadas
+!     ludcmp_dp, lubksb_dp
+
+!!!!!!!!!!!!!!!Comienza el programa !!!!!!!!!!!!!!
+      
+      do i=1,n
+         b(i) = bb(i+1)
+         do j=1,n
+            a(i,j) = aa(i+1,j+1)
+         enddo
+      enddo
+
+                                ! Descomposicion de auxm1
+      call ludcmp_dp ( a, n, n, indx, d)
+
+                                ! Sustituciones foward y backwards para hallar la solucion
+      do i=1,n
+         x(i) = b(i)
+      enddo
+      call lubksb_dp( a, n, n, indx, x )
+
+      do i=1,n
+         xx(i+1) = x(i)
+      enddo
+
+      return
+      end
+
+
+c     *** mat_oper.f ***
+
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
 c     ***********************************************************************
       subroutine unit(a,n)
 c     store the unit value in the diagonal of a 
 c     ***********************************************************************
+      implicit none
       real*8 a(n,n)
       integer n,i,j,k
@@ -1795 +1843 @@
 
 c     ***********************************************************************
+      subroutine invdiag(a,b,n)
+c     inverse of a diagonal matrix 
+c     ***********************************************************************
+      implicit none
+
+      integer n,i,j,k
+      real*8 a(n,n),b(n,n)
+
+      do 1,i=2,n-1
+         do 2,j=2,n-1
+            if (i.eq.j) then
+               a(i,j) = 1.d0/b(i,i)
+            else
+               a(i,j)=0.0d0
+            end if
+ 2       continue
+ 1    continue
+      do k=1,n
+         a(n,k) = 0.0d0
+         a(1,k) = 0.0d0
+         a(k,1) = 0.0d0
+         a(k,n) = 0.0d0
+      end do
+      return
+      end
+
+
+c     ***********************************************************************
       subroutine samem (a,m,n)
 c     store the matrix m in the matrix a 
 c     ***********************************************************************
+      implicit none
       real*8 a(n,n),m(n,n)
       integer n,i,j,k
@@ -1813 +1890 @@
       return
       end 
+
+
 c     ***********************************************************************
       subroutine mulmv(a,b,c,n)
@@ -1825 +1904 @@
          sum=0.0d0
          do 2,j=2,n-1
-            sum=sum+ (b(i,j)) * (c(j))
+            sum = sum + b(i,j) * c(j)
  2       continue
          a(i)=sum
@@ -1834 +1913 @@
       end
 
-c     ***********************************************************************
-      subroutine mulmm(a,b,c,n)
-c     ***********************************************************************
-      real*8 a(n,n), b(n,n), c(n,n)
+
+c     ***********************************************************************
+      subroutine trucodiag(a,b,c,d,e,n)
+c     inputs: matrices b,c,d,e
+c     output: matriz diagonal a
+c     Operacion a realizar:  a = b * c^(-1) * d + e
+c     La matriz c va a ser invertida
+c     Todas las matrices de entrada son diagonales excepto b
+c     Aprovechamos esa condicion para invertir c, acelerar el calculo, y 
+c     ademas, para forzar que a sea diagonal
+c     ***********************************************************************
+      implicit none
+      real*8 a(n,n),b(n,n),c(n,n),d(n,n),e(n,n), sum
       integer n,i,j,k
-
-!       do i=2,n-1
-!         do j=2,n-1
-!           a(i,j)= 0.d00
-!           do k=2,n-1
-!               a(i,j) = a(i,j) + b(i,k) * c(k,j)
-!           end do
-!         end do
-!       end do
-      do j=2,n-1
-         do i=2,n-1
-            a(i,j)=0.d0
-         enddo
-         do k=2,n-1
-            do i=2,n-1
-               a(i,j)=a(i,j)+b(i,k)*c(k,j)
-            enddo
-         enddo
-      enddo
+      do 1,i=2,n-1
+         sum=0.0d0
+         do 2,j=2,n-1
+            sum=sum+ b(i,j) * d(j,j)/c(j,j)
+ 2       continue
+         a(i,i) = sum + e(i,i)
+ 1    continue
       do k=1,n
          a(n,k) = 0.0d0
@@ -1864 +1940 @@
          a(k,n) = 0.0d0
       end do
-
-      return
-      end
-
-c     ***********************************************************************
-      subroutine resmm(a,b,c,n)
-c     ***********************************************************************
-      real*8 a(n,n), b(n,n), c(n,n)
-      integer n,i,j,k
-
-      do i=2,n-1
-         do j=2,n-1
-            a(i,j)= b(i,j) - c(i,j)
-         end do
-      end do
-      do k=1,n
-         a(n,k) = 0.0d0
-         a(1,k) = 0.0d0
-         a(k,1) = 0.0d0
-         a(k,n) = 0.0d0
-      end do
-
-      return
-      end
+      return
+      end
+
+
+c     ***********************************************************************
+      subroutine trucommvv(v,b,c,u,w,n)
+c     inputs: matrices b,c , vectores u,w
+c     output: vector v 
+c     Operacion a realizar:  v = b * c^(-1) * u + w
+c     La matriz c va a ser invertida 
+c     c es diagonal, b no
+c     Aprovechamos esa condicion para invertir c, y acelerar el calculo
+c     ***********************************************************************
+      implicit none
+      real*8 v(n),b(n,n),c(n,n),u(n),w(n), sum
+      integer n,i,j
+      do 1,i=2,n-1
+         sum=0.0d0
+         do 2,j=2,n-1
+            sum=sum+ b(i,j) * u(j)/c(j,j)
+ 2       continue
+         v(i) = sum + w(i)
+ 1    continue
+      v(1) = 0.d0
+      v(n) = 0.d0
+      return
+      end
+
+
+c     ***********************************************************************
+      subroutine sypvmv(v,u,c,w,n)
+c     inputs: matriz diagonal c , vectores u,w
+c     output: vector v 
+c     Operacion a realizar:  v = u + c * w 
+c     ***********************************************************************
+      implicit none
+      real*8 v(n),u(n),c(n,n),w(n)
+      integer n,i
+      do 1,i=2,n-1
+         v(i)= u(i) + c(i,i) * w(i)
+ 1    continue
+      v(1) = 0.0d0
+      v(n) = 0.0d0
+      return
+      end
+
 
 c     ***********************************************************************
@@ -1899 +1997 @@
 
       do 1,i=2,n-1
-         a(i)= (b(i)) + (c(i))
+         a(i)= b(i) + c(i)
  1    continue
       a(1) = 0.0d0
@@ -1906 +2004 @@
       end
 
-c     ***********************************************************************
-      subroutine zerom(a,n)
+
+c     ***********************************************************************
+      subroutine sypvvv(a,b,c,d,n)
+c     a(i)=b(i)+c(i)*d(i)
+c     ***********************************************************************
+      implicit none
+      real*8 a(n),b(n),c(n),d(n)
+      integer n,i
+      do 1,i=2,n-1
+         a(i)= b(i) + c(i) * d(i)
+ 1    continue
+      a(1) = 0.0d0
+      a(n) = 0.0d0
+      return
+      end
+
+
+c     ***********************************************************************
+!      subroutine zerom(a,n)
 c     a(i,j)= 0.0
 c     ***********************************************************************
-
-      implicit none
-
-      integer n,i,j
-      real*8 a(n,n)
-
-      do 1,i=1,n
-         do 2,j=1,n
-            a(i,j) = 0.0d0
- 2       continue
- 1    continue
-      return
-      end
+!      implicit none
+!      integer n,i,j
+!      real*8 a(n,n)
+
+!      do 1,i=1,n
+!         do 2,j=1,n
+!           a(i,j) = 0.0d0
+! 2       continue
+! 1    continue
+!      return
+!      end
+
 
 c     ***********************************************************************
@@ -1928 +2042 @@
 c     a(i,j) = b(i,j) = c(i,j) = d(i,j) = 0.0 
 c     ***********************************************************************
+      implicit none
       real*8 a(n,n), b(n,n), c(n,n), d(n,n)
-      integer n,i,j
-      do 1,i=1,n
-         do 2,j=1,n
-            a(i,j) = 0.0d0
-            b(i,j) = 0.0d0
-            c(i,j) = 0.0d0
-            d(i,j) = 0.0d0
- 2       continue
- 1    continue
-      return
-      end
+      integer n
+      a(1:n,1:n)=0.d0
+      b(1:n,1:n)=0.d0
+      c(1:n,1:n)=0.d0
+      d(1:n,1:n)=0.d0
+!      do 1,i=1,n
+!         do 2,j=1,n
+!           a(i,j) = 0.0d0
+!           b(i,j) = 0.0d0
+!           c(i,j) = 0.0d0
+!           d(i,j) = 0.0d0
+! 2       continue
+! 1    continue
+      return
+      end
+
 
 c     ***********************************************************************
@@ -1945 +2065 @@
 c     a(i,j) = b(i,j) = c(i,j) = 0.0 
 c     **********************************************************************
+      implicit none
       real*8 a(n,n), b(n,n), c(n,n)
-      integer n,i,j
-      do 1,i=1,n
-         do 2,j=1,n
-            a(i,j) = 0.0d0
-            b(i,j) = 0.0d0
-            c(i,j) = 0.0d0
- 2       continue
- 1    continue
-      return
-      end
+      integer n
+      a(1:n,1:n)=0.d0
+      b(1:n,1:n)=0.d0
+      c(1:n,1:n)=0.d0
+!      do 1,i=1,n
+!         do 2,j=1,n
+!           a(i,j) = 0.0d0
+!           b(i,j) = 0.0d0
+!           c(i,j) = 0.0d0
+! 2       continue
+! 1    continue
+      return
+      end
+
 
 c     ***********************************************************************
@@ -1961 +2086 @@
 c     a(i,j) = b(i,j) = 0.0 
 c     ***********************************************************************
+      implicit none
       real*8 a(n,n), b(n,n)
-      integer n,i,j
-      do 1,i=1,n
-         do 2,j=1,n
-            a(i,j) = 0.0d0
-            b(i,j) = 0.0d0
- 2       continue
- 1    continue
-      return
-      end
-c     ***********************************************************************
-      subroutine zerov(a,n)
+      integer n
+      a(1:n,1:n)=0.d0
+      b(1:n,1:n)=0.d0
+!      do 1,i=1,n
+!         do 2,j=1,n
+!           a(i,j) = 0.0d0
+!           b(i,j) = 0.0d0
+! 2       continue
+! 1    continue
+      return
+      end
+
+
+c     ***********************************************************************
+!      subroutine zerov(a,n)
 c     a(i)= 0.0
 c     ***********************************************************************
-      real*8 a(n)
-      integer n,i
-      do 1,i=1,n
-         a(i) = 0.0d0
- 1    continue
-      return
-      end
+!      implicit none
+!      real*8 a(n)
+!      integer n,i
+!      do 1,i=1,n
+!         a(i) = 0.0d0
+! 1    continue
+!      return
+!      end
+
+
 c     ***********************************************************************
       subroutine zero4v(a,b,c,d,n)
 c     a(i) = b(i) = c(i) = d(i,j) = 0.0
 c     ***********************************************************************
+      implicit none
       real*8 a(n), b(n), c(n), d(n)
-      integer n,i
-      do 1,i=1,n
-         a(i) = 0.0d0
-         b(i) = 0.0d0
-         c(i) = 0.0d0
-         d(i) = 0.0d0
- 1    continue
-      return
-      end
+      integer n
+      a(1:n)=0.d0
+      b(1:n)=0.d0
+      c(1:n)=0.d0
+      d(1:n)=0.d0
+!      do 1,i=1,n
+!         a(i) = 0.0d0
+!         b(i) = 0.0d0
+!         c(i) = 0.0d0
+!         d(i) = 0.0d0
+! 1    continue
+      return
+      end
+
+
 c     ***********************************************************************
       subroutine zero3v(a,b,c,n)
 c     a(i) = b(i) = c(i) = 0.0
 c     ***********************************************************************
+      implicit none
       real*8 a(n), b(n), c(n) 
-      integer n,i
-      do 1,i=1,n
-         a(i) = 0.0d0
-         b(i) = 0.0d0
-         c(i) = 0.0d0
- 1    continue
-      return
-      end
+      integer n
+      a(1:n)=0.d0
+      b(1:n)=0.d0
+      c(1:n)=0.d0
+!      do 1,i=1,n
+!         a(i) = 0.0d0
+!         b(i) = 0.0d0
+!         c(i) = 0.0d0
+! 1    continue
+      return
+      end
+
+
 c     ***********************************************************************
       subroutine zero2v(a,b,n)
 c     a(i) = b(i) = 0.0
 c     ***********************************************************************
+      implicit none
       real*8 a(n), b(n) 
-      integer n,i
-      do 1,i=1,n
-         a(i) = 0.0d0
-         b(i) = 0.0d0
- 1    continue
-      return
-      end
-
-
+      integer n
+      a(1:n)=0.d0
+      b(1:n)=0.d0
+!      do 1,i=1,n
+!         a(i) = 0.0d0
+!         b(i) = 0.0d0
+! 1    continue
+      return
+      end
+
+c     ***********************************************************************
+
+
+c****************************************************************************
+
+c     *** suaviza.f ***
+
+c*****************************************************************************
+c     
+      subroutine suaviza ( x, n, ismooth, y )
+c     
+c     x - input and return values 
+c     y - auxiliary vector
+c     ismooth = 0  --> no smoothing is performed
+c     ismooth = 1  --> weak smoothing (5 points, centred weighted)
+c     ismooth = 2  --> normal smoothing (3 points, evenly weighted)
+c     ismooth = 3  --> strong smoothing (5 points, evenly weighted)
+
+
+c     august 1991
+c*****************************************************************************
+
+      implicit none
+
+      integer   n, imax, imin, i, ismooth
+      real*8    x(n), y(n)
+c*****************************************************************************
+
+      imin=1
+      imax=n
+
+      if (ismooth.eq.0) then
+
+         return 
+
+      elseif (ismooth.eq.1) then ! 5 points, with central weighting 
+
+         do i=imin,imax 
+            if(i.eq.imin)then
+               y(i)=x(imin)
+            elseif(i.eq.imax)then
+               y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0
+            elseif(i.gt.(imin+1) .and. i.lt.(imax-1) )then
+               y(i) = ( x(i+2)/4.d0 + x(i+1)/2.d0 + 2.d0*x(i)/3.d0 + 
+     @              x(i-1)/2.d0 + x(i-2)/4.d0 )* 6.d0/13.d0
+            else
+               y(i)=(x(i+1)/2.d0+x(i)+x(i-1)/2.d0)/2.d0
+            end if
+         end do 
+
+      elseif (ismooth.eq.2) then ! 3 points, evenly spaced
+
+         do i=imin,imax 
+            if(i.eq.imin)then
+               y(i)=x(imin)
+            elseif(i.eq.imax)then
+               y(i)=x(imax-1)+(x(imax-1)-x(imax-3))/2.d0
+            else
+               y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0
+            end if
+         end do 
+         
+      elseif (ismooth.eq.3) then ! 5 points, evenly spaced
+
+         do i=imin,imax 
+            if(i.eq.imin)then
+               y(i) = x(imin)
+            elseif(i.eq.(imin+1) .or. i.eq.(imax-1))then
+               y(i) = ( x(i+1)+x(i)+x(i-1) )/3.d0 
+            elseif(i.eq.imax)then
+               y(i) = ( x(imax-1) + x(imax-1) + x(imax-2) ) / 3.d0
+            else
+               y(i) = ( x(i+2)+x(i+1)+x(i)+x(i-1)+x(i-2) )/5.d0
+            end if
+         end do 
+
+      else
+
+         write (*,*) ' Error in suaviza.f   Wrong ismooth value.'
+         stop
+
+      endif
+
+c     rehago el cambio, para devolver x(i)
+      do i=imin,imax 
+         x(i)=y(i)
+      end do
+
+      return
+      end
+
+
+c     ***********************************************************************
+      subroutine mulmmf90(a,b,c,n)
+c     ***********************************************************************
+      implicit none
+      real*8 a(n,n), b(n,n), c(n,n)
+      integer n
+
+      a=matmul(b,c)
+      a(1,:)=0.d0
+      a(:,1)=0.d0
+      a(n,:)=0.d0
+      a(:,n)=0.d0
+
+      return
+      end
+
+
+c     ***********************************************************************
+      subroutine resmmf90(a,b,c,n)
+c     ***********************************************************************
+      implicit none
+      real*8 a(n,n), b(n,n), c(n,n)
+      integer n
+
+      a=b-c
+      a(1,:)=0.d0
+      a(:,1)=0.d0
+      a(n,:)=0.d0
+      a(:,n)=0.d0
+
+      return
+      end
+
+
+c*******************************************************************
+
+      subroutine gethist_03 (ihist) 
+
+c*******************************************************************
+
+      implicit none
+
+      include   'nlte_paramdef.h'
+      include   'nlte_commons.h'
+
+
+c     arguments 
+      integer         ihist
+
+c     local variables
+      integer   j, r, mm
+      real*8          xx
+
+c     ***************
+
+      nbox = nbox_stored(ihist)
+      do j=1,mm_stored(ihist)
+         thist(j) = thist_stored(ihist,j) 
+         do r=1,nbox_stored(ihist)
+            no(r) = no_stored(ihist,r) 
+            sk1(j,r) = sk1_stored(ihist,j,r) 
+            xls1(j,r) = xls1_stored(ihist,j,r)
+            xld1(j,r) = xld1_stored(ihist,j,r)
+         enddo
+      enddo
+
+
+      return
+      end
+
+
+c     *******************************************************************
+
+      subroutine rhist_03 (ihist) 
+
+c     *******************************************************************
+
+      implicit none
+
+      include   'nlte_paramdef.h'
+      include   'nlte_commons.h'
+
+
+c     arguments 
+      integer         ihist
+
+c     local variables
+      integer   j, r, mm
+      real*8          xx
+
+c     ***************
+
+      open(unit=3,file=hisfile,status='old')
+
+      read(3,*)
+      read(3,*)
+      read(3,*) mm_stored(ihist)
+      read(3,*)
+      read(3,*) nbox_stored(ihist)
+      read(3,*)
+
+      if ( nbox_stored(ihist) .gt. nbox_max ) then
+         write (*,*) ' nbox too large in input file ', hisfile
+         stop ' Check maximum number nbox_max in mz1d.par '
+      endif
+
+      do j=1,mm_stored(ihist)
+         read(3,*) thist_stored(ihist,j) 
+         do r=1,nbox_stored(ihist)
+            read(3,*) no_stored(ihist,r), 
+     &           sk1_stored(ihist,j,r), 
+     &           xls1_stored(ihist,j,r),
+     &           xx,
+     &           xld1_stored(ihist,j,r)
+         enddo
+
+      enddo
+
+      close(unit=3)
+
+
+      return
+      end