source: trunk/LMDZ.MARS/libf/phymars/ludcmp_dp.F @ 481

Last change on this file since 481 was 414, checked in by aslmd, 14 years ago

LMDZ.MARS : NEW NLTE MODEL FROM GRANADA AMIGOS

23/11/11 == FGG + MALV

New parameterization of the NLTE 15 micron cooling. The old parameterization is kept as an option, including or not variable atomic oxygen concentration. A new flag is introduced in callphys.def, nltemodel, to select which parameterization wants to be used (new one, old one with variable [O], or old one with fixed [O], see below). Includes many new subroutines and commons in phymars. Some existing routines are also modified:

-physiq.F. Call to the new subroutine NLTE_leedat in first call. Call to nltecool modified to allow for variable atomic oxygen. Depending on the value of nltemodel, the new subroutine NLTEdlvr09_TCOOL is called instead of nltecool.

-inifis.F. Reading of nltemodel is added.

-callkeys.h Declaration of nltemodel is added.

The following lines should be added to callphys.def (ideally after setting callnlte):

# NLTE 15um scheme to use.
# 0-> Old scheme, static oxygen
# 1-> Old scheme, dynamic oxygen
# 2-> New scheme
nltemodel = 2

A new directory, NLTEDAT, has to be included in datagcm.

Improvements into NLTE NIR heating parameterization to take into account variability with O/CO2 ratio and SZA. A new subroutine, NIR_leedat.F, and a new common, NIRdata.h, are included. A new flag, nircorr, is added in callphys.def, to include or not these corrections. The following files are modified:

-nirco2abs.F: nq and pq are added in the arguments. The corrections factors are interpolated to the GCM grid and included in the clculation. A new subroutine, interpnir, is added at the end of the file.

-physiq.F: Call to NIR_leedat added at first call. Modified call to nirco2abs

-inifis: Reading new flag nircorr.

-callkeys.h: Declaration of nircorr.

The following lines have to be added to callphys.def (ideally after callnirco2):

# NIR NLTE correction for variable SZA and O/CO2?
# matters only if callnirco2=T
# 0-> no correction
# 1-> include correction
nircorr=1

A new data file, NIRcorrection_feb2011.dat, has to be included in datagcm.

Small changes to the molecular diffusion scheme to fix the number of species considered, to avoid problems when compiling with more than 15 tracers (for example, when CH4 is included). Modified subroutines: aeronomars/moldiff.F and aeronomars/moldiffcoeff.F
File size: 4.6 KB
Line 
1      subroutine ludcmp_dp(a,n,np,indx,d)
2
3c       jul 2011 malv+fgg
4
5      implicit none
6
7      integer n, np
8      real*8 a(np,np), d
9      integer indx(n)
10     
11      integer nmax, i, j, k, imax
12      real*8 tiny
13      parameter (nmax=100,tiny=1.0d-20)                                       
14      real*8 vv(nmax), aamax, sum, dum
15
16
17      d=1.0d0
18      do 12 i=1,n                                                             
19        aamax=0.0d0
20        do 11 j=1,n                                                           
21          if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))                         
2211      continue                                                             
23        if (aamax.eq.0.0) pause 'singular matrix.'                         
24        vv(i)=1.0d0/aamax                         
2512    continue                                                               
26      do 19 j=1,n                                                             
27        if (j.gt.1) then                                                     
28          do 14 i=1,j-1                                                       
29            sum=a(i,j)                                                       
30            if (i.gt.1)then                                                   
31              do 13 k=1,i-1                                                   
32                sum=sum-a(i,k)*a(k,j)                                         
3313            continue                                                       
34              a(i,j)=sum                                                     
35            endif                                                             
3614        continue                                                           
37        endif                                                                 
38        aamax=0.0d0                                                           
39        do 16 i=j,n                                                           
40          sum=a(i,j)                                                         
41          if (j.gt.1)then                                                     
42            do 15 k=1,j-1                                                     
43              sum=sum-a(i,k)*a(k,j)                                           
4415          continue                                                         
45            a(i,j)=sum                                                       
46          endif                                                               
47          dum=vv(i)*abs(sum)                                                 
48          if (dum.ge.aamax) then                                             
49            imax=i                                                           
50            aamax=dum                                                         
51          endif                                                               
5216      continue                                                             
53        if (j.ne.imax)then                                                   
54          do 17 k=1,n                                                         
55            dum=a(imax,k)                                                     
56            a(imax,k)=a(j,k)                                                 
57            a(j,k)=dum                                                       
5817        continue                                                           
59          d=-d                                                               
60          vv(imax)=vv(j)                                                     
61        endif                                                                 
62        indx(j)=imax                                                         
63        if(j.ne.n)then                                                       
64          if(a(j,j).eq.0.0)a(j,j)=tiny                                     
65          dum=1.0d0/a(j,j)                                                 
66          do 18 i=j+1,n                                                       
67            a(i,j)=a(i,j)*dum                                                 
6818        continue                                                           
69        endif                                                                 
7019    continue                                                               
71      if(a(n,n).eq.0.0)a(n,n)=tiny                                     
72      return                                                                 
73      end                                                                     
Note: See TracBrowser for help on using the repository browser.