Ignore:
Timestamp:
May 3, 2012, 8:44:26 AM (13 years ago)
Author:
emillour
Message:

Mars GCM:

  • Update of molecular diffusion routines (diffusion.h, moldiff_red.F90 and moldiffcoeff_red.F): optimized to run faster and diffuse only specific tracers (if they are present).
  • bugs fixed in simpleclouds.F (extra arguments) and watercloud.F (wrong arguments in call to simpleclouds and rdust() should only be recomputed if dust_mass & dust_number tracers are available).

JYC+EM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90

    r552 r645  
    3333     
    3434
    35       real hco2(ncompdiff),ho
    36 
    37       integer ig,nz,l,k,n,nn,iq,kn,p,il0,kl,igen,igen2,kl2,ngen,iter,il1,ij0
     35!      real hco2(ncompdiff),ho
     36
     37      integer,dimension(nqmx) :: indic_diff
     38      integer ig,nz,l,k,n,nn,p,ij0
    3839      integer istep,il,gcn,ntime,nlraf
    3940      real*8 masse
    4041      real*8 masseU,kBolt,RgazP,Rmars,Grav,Mmars
    41       real*8 rho0,D0,T0,H0,time0,dZ,time,dZraf,tdiff,Hup,Zmin,Zmax
    42       real*8 hh,dcoef,dcoef1,ptfac, ntot, dens,mgcn,FacEsc
     42      real*8 rho0,D0,T0,H0,time0,dZ,time,dZraf,tdiff,Zmin,Zmax
     43      real*8 FacEsc,invsgmu
    4344      real*8 hp(nlayermx)
    4445      real*8 pp(nlayermx)
    … …  
    4647      real*8 tt(nlayermx),tnew(nlayermx),tint(nlayermx)
    4748      real*8 zz(nlayermx)
    48       real*8 qq(nlayermx,ncompdiff)
    49       real*8 qnew(nlayermx,ncompdiff)
    50       real*8 qint(nlayermx,ncompdiff),FacMass(nlayermx,ncompdiff)
    51       real*8 rhoK(nlayermx,ncompdiff),rhoT(nlayermx)
    52       real*8 rhoKinit(nlayermx,ncompdiff)
     49      real*8,dimension(:,:),allocatable :: qq,qnew,qint,FacMass
     50      real*8,dimension(:,:),allocatable :: rhoK,rhokinit
     51      real*8 rhoT(nlayermx)
    5352      real*8 dmmeandz(nlayermx)
    5453      real*8 massemoy(nlayermx)
    … …  
    6160      real*8,dimension(:),allocatable :: Atri,Btri,Ctri,Dtri,Xtri,Tad,Dad,Zad,rhoad
    6261      real*8,dimension(:),allocatable :: alpha,beta,gama,delta,eps
    63       real*8 wi(ncompdiff),Wad(ncompdiff),Uthermal(ncompdiff)
    64       real*8 Lambdaexo(ncompdiff)
    65       real*8 MToT1(ncompdiff),Mtot2(ncompdiff)
    66       real*8 MRaf1(ncompdiff),Mraf2(ncompdiff)
     62      real*8,dimension(:),allocatable ::  wi,Wad,Uthermal,Lambdaexo,Hspecie
     63      real*8,dimension(:),allocatable :: Mtot1,Mtot2,Mraf1,Mraf2
     64      character(len=20),dimension(14) :: ListeDiff
    6765!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    6866!     tracer numbering in the molecular diffusion
    6967!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    70 !  Atomic oxygen must always be the LAST species of the list as
    71 ! it is the dominant species at high altitudes. 
    72       integer,parameter :: i_o   = 1
    73       integer,parameter :: i_n2   = 2
    74       integer,parameter :: i_co   = 3
    75       integer,parameter :: i_ar  = 4
    76       integer,parameter :: i_h2   = 5
    77       integer,parameter :: i_h    = 6
    78       integer,parameter :: i_o2   = 7
    79       integer,parameter :: i_oh  = 8
    80       integer,parameter :: i_ho2  = 9
    81       integer,parameter :: i_h2o = 10
    82       integer,parameter :: i_h2o2  = 11
    83       integer,parameter :: i_o1d   = 12
    84 !      integer,parameter :: i_o3   = 13
    85       integer,parameter :: i_n    = 13
    86       integer,parameter :: i_no   = 14
    87       integer,parameter :: i_no2  = 15
    88 !      integer,parameter :: i_n2d  = 17
    89 !      integer,parameter :: i_oplus = 18
    90       integer,parameter :: i_co2    = 16
    91 !      integer,parameter :: i_oplus = 17
    92 !      integer,parameter :: i_hplus = 18
     68!  We need the index of escaping species
     69
     70      integer,save :: i_h2 
     71      integer,save :: i_h
     72! vertical index limit for the molecular diffusion
     73      integer,save :: il0 
    9374
    9475! Tracer indexes in the GCM:
    95       integer,save :: g_co2=0
    96       integer,save :: g_co=0
    97       integer,save :: g_o=0
    98       integer,save :: g_o1d=0
    99       integer,save :: g_o2=0
    100       integer,save :: g_o3=0
    101       integer,save :: g_h=0
    102       integer,save :: g_h2=0
    103       integer,save :: g_oh=0
    104       integer,save :: g_ho2=0
    105       integer,save :: g_h2o2=0
    106       integer,save :: g_n2=0
    107       integer,save :: g_ar=0
    108       integer,save :: g_h2o=0
    109       integer,save :: g_n=0
    110       integer,save :: g_no=0
    111       integer,save :: g_no2=0
    112       integer,save :: g_n2d=0
     76!      integer,save :: g_co2=0
     77!      integer,save :: g_co=0
     78!      integer,save :: g_o=0
     79!      integer,save :: g_o1d=0
     80!      integer,save :: g_o2=0
     81!      integer,save :: g_o3=0
     82!      integer,save :: g_h=0
     83!      integer,save :: g_h2=0
     84!      integer,save :: g_oh=0
     85!      integer,save :: g_ho2=0
     86!      integer,save :: g_h2o2=0
     87!      integer,save :: g_n2=0
     88!      integer,save :: g_ar=0
     89!      integer,save :: g_h2o=0
     90!      integer,save :: g_n=0
     91!      integer,save :: g_no=0
     92!      integer,save :: g_no2=0
     93!      integer,save :: g_n2d=0
    11394!      integer,save :: g_oplus=0
    11495!      integer,save :: g_hplus=0
    11596
    116       integer,save :: gcmind(ncompdiff) ! array of GCM indexes
     97      integer,save :: ncompdiff
     98      integer,dimension(:),allocatable,save :: gcmind ! array of GCM indexes
    11799      integer ierr
    118100
    119101      logical,save :: firstcall=.true.
    120       real abfac(ncompdiff)
    121       real,save :: dij(ncompdiff,ncompdiff)
     102!      real abfac(ncompdiff)
     103      real,dimension(:,:),allocatable,save :: dij
    122104      real,save :: step
    123105
    124106! Initializations at first call
    125107      if (firstcall) then
    126         call moldiffcoeff_red(dij)
     108
     109! Liste des especes qui sont diffuses si elles sont presentes
     110!       ListeDiff=['co2','o','n2','ar','co','h2','h','d2','d','hd','o2','oh','ho2','h2o_vap','h2o2','o1d','n','no','no2']
     111        ListeDiff(1)='co2'
     112        ListeDiff(2)='o'
     113        ListeDiff(3)='n2'
     114        ListeDiff(4)='ar'
     115        ListeDiff(5)='co'
     116        ListeDiff(6)='h2'
     117        ListeDiff(7)='h'
     118        ListeDiff(8)='d2'
     119        ListeDiff(9)='hd'
     120        ListeDiff(10)='d'
     121        ListeDiff(11)='o2'
     122        ListeDiff(12)='h2o_vap'
     123        ListeDiff(13)='o3'
     124        ListeDiff(14)='n'
     125        i_h=1000
     126        i_h2=1000
     127! On regarde quelles especes diffusables sont presentes
     128
     129        ncompdiff=0
     130        indic_diff(1:nqmx)=0
     131       
     132        do nn=1,nqmx
     133        do n=1,14
     134        if (ListeDiff(n) .eq. noms(nn)) then
     135        indic_diff(nn)=1
     136        if (noms(nn) .eq. 'h') i_h=n
     137        if (noms(nn) .eq. 'h2') i_h2=n
     138        endif
     139
     140        enddo
     141        if (indic_diff(nn) .eq. 1) then
     142        print*,'specie', noms(nn), 'diffused in moldiff_red'
     143        ncompdiff=ncompdiff+1
     144        endif
     145        enddo
     146        print*,'number of diffused species:',ncompdiff
     147        allocate(gcmind(ncompdiff))
     148
     149! Store gcm indexes in gcmind
     150        n=0
     151        do nn=1,nqmx
     152        if (indic_diff(nn) .eq. 1) then
     153        n=n+1
     154        gcmind(n)=nn
     155        endif
     156        enddo
     157
     158        print*,'gcmind',gcmind
     159
     160! find vertical index above which diffusion is computed
     161
     162        do l=1,nlayermx
     163        if (pplay(1,l) .gt. Pdiff) then
     164        il0=l
     165        endif
     166        enddo
     167        il0=il0+1
     168        print*,'vertical index for diffusion',il0,pplay(1,il0)
     169
     170        allocate(dij(ncompdiff,ncompdiff))
     171
     172        call moldiffcoeff_red(dij,indic_diff,gcmind,ncompdiff)
    127173        print*,'MOLDIFF  EXO'
    128 
    129         ! identify the indexes of the tracers we'll need
    130         g_co2=igcm_co2
    131         if (g_co2.eq.0) then
    132           write(*,*) "moldiff: Error; no CO2 tracer !!!"
    133           stop
    134         endif
    135         g_n2=igcm_n2
    136         if (g_n2.eq.0) then
    137           write(*,*) "moldiff: Error; no N2 tracer !!!"
    138           stop
    139         endif
    140         g_ar=igcm_ar
    141         if (g_ar.eq.0) then
    142           write(*,*) "moldiff: Error; no Ar tracer !!!"
    143           stop
    144         endif
    145         g_h2=igcm_h2
    146         if (g_h2.eq.0) then
    147           write(*,*) "moldiff: Error; no H2 tracer !!!"
    148           stop
    149         endif
    150         g_h=igcm_h
    151         if (g_h.eq.0) then
    152           write(*,*) "moldiff: Error; no H tracer !!!"
    153           stop
    154         endif
    155         g_co=igcm_co
    156         if (g_co.eq.0) then
    157           write(*,*) "moldiff: Error; no CO tracer !!!"
    158           stop
    159         endif
    160         g_o2=igcm_o2
    161         if (g_o2.eq.0) then
    162           write(*,*) "moldiff: Error; no O2 tracer !!!"
    163           stop
    164         endif
    165         g_oh=igcm_oh
    166         if (g_oh.eq.0) then
    167           write(*,*) "moldiff: Error; no OH tracer !!!"
    168           stop
    169         endif
    170         g_ho2=igcm_ho2
    171         if (g_ho2.eq.0) then
    172           write(*,*) "moldiff: Error; no HO2 tracer !!!"
    173           stop
    174         endif
    175         g_h2o=igcm_h2o_vap
    176         if (g_h2o.eq.0) then
    177           write(*,*) "moldiff: Error; no H2O tracer !!!"
    178           stop
    179         endif
    180         g_h2o2=igcm_h2o2
    181         if (g_h2o2.eq.0) then
    182           write(*,*) "moldiff: Error; no H2O2 tracer !!!"
    183           stop
    184         endif
    185         g_o1d=igcm_o1d
    186         if (g_o1d.eq.0) then
    187           write(*,*) "moldiff: Error; no O(1D) tracer !!!"
    188           stop
    189         endif
    190         g_o3=igcm_o3
    191         if (g_o3.eq.0) then
    192           write(*,*) "moldiff: Error; no O3 tracer !!!"
    193           stop
    194         endif
    195         g_n=igcm_n
    196         if (g_n.eq.0) then
    197           write(*,*) "moldiff: Error; no N tracer !!!"
    198           stop
    199         endif
    200         g_no=igcm_no
    201         if (g_no.eq.0) then
    202           write(*,*) "moldiff: Error; no NO tracer !!!"
    203           stop
    204         endif
    205         g_no2=igcm_no2
    206         if (g_no2.eq.0) then
    207           write(*,*) "moldiff: Error; no NO2 tracer !!!"
    208           stop
    209         endif
    210         g_n2d=igcm_n2d
    211         if (g_n2d.eq.0) then
    212           write(*,*) "moldiff: Error; no N2(D) tracer !!!"
    213           stop
    214         endif
    215         g_o=igcm_o
    216         if (g_o.eq.0) then
    217           write(*,*) "moldiff: Error; no O tracer !!!"
    218           stop
    219         endif
    220 
    221 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    222 !   fill array to relate local indexes to gcm indexes
    223 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    224 
    225         gcmind(i_co2)= g_co2
    226         gcmind(i_n2)= g_n2
    227         gcmind(i_ar)= g_ar
    228         gcmind(i_h2)= g_h2
    229         gcmind(i_h)= g_h
    230         gcmind(i_co)= g_co
    231         gcmind(i_o2)= g_o2
    232         gcmind(i_oh)= g_oh
    233         gcmind(i_ho2)= g_ho2
    234         gcmind(i_h2o)= g_h2o
    235         gcmind(i_h2o2)= g_h2o2
    236         gcmind(i_o1d)= g_o1d
    237 !        gcmind(i_o3)=g_o3
    238         gcmind(i_n)=g_n
    239         gcmind(i_no)=g_no
    240         gcmind(i_no2)=g_no2
    241 !        gcmind(i_n2d)=g_n2d
    242         gcmind(i_o)=g_o
    243174
    244175        firstcall= .false.
    … …  
    257188        Grav=6.67d-11
    258189        Mmars=6.4d23
    259       il0=idifflow ! Diffusion limited to regions above il0 JYC 
    260         ij0=302 ! For test
     190        ij0=6000 ! For test
     191        invsgmu=1d0/g/masseU
    261192       
     193! allocatation des tableaux dependants du nombre d especes diffusees
     194        allocate(qq(nlayermx,ncompdiff))
     195        allocate(qnew(nlayermx,ncompdiff))
     196        allocate(qint(nlayermx,ncompdiff))
     197        allocate(FacMass(nlayermx,ncompdiff))
     198        allocate(rhok(nlayermx,ncompdiff))
     199        allocate(rhokinit(nlayermx,ncompdiff))
     200
     201        allocate(wi(ncompdiff))
     202        allocate(wad(ncompdiff))
     203        allocate(uthermal(ncompdiff))
     204        allocate(lambdaexo(ncompdiff))
     205        allocate(Hspecie(ncompdiff))
     206        allocate(Mtot1(ncompdiff))
     207        allocate(Mtot2(ncompdiff))
     208        allocate(Mraf1(ncompdiff))
     209        allocate(Mraf2(ncompdiff))
     210
     211!       print*,'moldiff',i_h2,i_h,ncompdiff
    262212      do ig=1,ngridmx
    263 
    264213        pp=dble(pplay(ig,:))
    265214
    … …  
    310259        enddo
    311260
    312         Zmin=zz(idifflow)
     261        Zmin=zz(il0)
    313262        Zmax=zz(nlayermx)
    314263
    … …  
    348297        allocate(alpha(nlraf),beta(nlraf),gama(nlraf),delta(nlraf),eps(nlraf))
    349298
    350 ! before beginning, I use a better vertical resolution above il1,
     299! before beginning, I use a better vertical resolution above il0,
    351300! altitude grid reinterpolation
    352301! The diffusion is solved on an altitude grid with constant step dz
    … …  
    404353!       enddo
    405354
    406 ! The diffusion is computed above il0 (defined in diffusion.h)
     355! The diffusion is computed above il0 computed from Pdiff in diffusion.h
    407356! No change below il0
    408357       
    … …  
    435384! Loop on species
    436385
     386        T0=Traf(nlraf)
     387        rho0=1d0
     388
    437389        do nn=1,ncompdiff
    438390        masse=dble(mmol(gcmind(nn)))
    … …  
    443395        CALL HSCALEREAL(nn,Nrafk,Dzraf,Hraf,nlraf,ncompdiff)
    444396        Hrafmol(:,nn)=Hraf(:)
     397!       Hspecie(nn)=kbolt*T0/masse*invsgmu
    445398! Computation of the diffusion vertical velocity of the specie
    446399        CALL VELVERT(nn,Traf,Hraf,Draf,Dzraf,masse,Wraf,nlraf)
    … …  
    451404        enddo
    452405! We use a lower time step
    453         Tdiff=minval(Tdiffrafmol)/10d0
    454 
     406        Tdiff=minval(Tdiffrafmol)
     407        Tdiff=minval(Tdiffrafmol(nlraf,:))*Mraf(nlraf)
    455408! Some problems when H is dominant
    456409! The time step is chosen function of atmospheric mass at higher level
    457410! It is not very nice
    458411
    459         if (tdiff .lt. tdiffmin) tdiff=tdiffmin*Mraf(nlraf)
     412!       if (ig .eq. ij0) then
     413!       print*,'test',ig,tdiff,tdiffmin,minloc(Tdiffrafmol),minloc(Tdiffrafmol(nlraf,:))
     414!       endif
     415        if (tdiff .lt. tdiffmin*Mraf(nlraf)) tdiff=tdiffmin*Mraf(nlraf)
    460416
    461417! Number of time step
    462418        ntime=anint(dble(ptimestep)/tdiff)
     419!       print*,'ptime',ig,ptimestep,tdiff,ntime,tdiffmin,Mraf(nlraf)
     420! Adimensionned temperature
     421
     422        do l=1,nlraf
     423        Tad(l)=Traf(l)/T0
     424        enddo
    463425
    464426        do istep=1,ntime
    … …  
    469431! Parameters to adimension the problem
    470432
    471         rho0=1d0
    472         T0=Traf(nlraf)
    473         H0=kbolt*T0/dble(mmol(gcmind(nn)))/g/masseU
     433        H0=kbolt*T0/dble(mmol(gcmind(nn)))*invsgmu
    474434        D0=Draf(nlraf)
    475435        Time0=H0*H0/D0
    … …  
    479439
    480440        do l=1,nlraf
    481         Tad(l)=Traf(l)/T0
    482441        Dad(l)=Draf(l)/D0
    483442        Zad(l)=Zraf(l)/H0
    484         RhoAd(l)=Rrafk(l,nn)/rho0
    485443        enddo
    486444        Wad(nn)=wi(nn)*Time0/H0
    487445        DZ=Zad(2)-Zad(1)
     446        FacEsc=exp(-wad(nn)*DZ/Dad(nlraf-1))
     447
     448        do l=1,nlraf
     449        RhoAd(l)=Rrafk(l,nn)/rho0
     450        enddo
    488451
    489452! Compute intermediary coefficients
    … …  
    494457! Compute escape factor (exp(-ueff*dz/D(nz)))
    495458
    496         FacEsc=exp(-wad(nn)*DZ/Dad(nlraf-1))
    497 
    498459! Compute matrix coefficients
    499460
    … …  
    509470!       Xtri=rhoAd
    510471
    511 !       if (nn .eq. 6) then
     472!       if (ig .eq. ij0 .and. (nn .eq. 1 .or. nn .eq. 5 .or. nn .eq. 6 .or. nn .eq. 16)) then
    512473!       do l=1,nlraf
    513474!       if (Xtri(l) .lt. 0.) then
    514 !       print*,'l',l,rhoAd(l)*rho0,Xtri(l)*rho0,nn
     475!       print*,'l',l,rhoAd(l)*rho0,Xtri(l)*rho0,nn,Tad(l),Zad(l),Dad(l)
    515476!       stop
    516477!       endif
    … …  
    533494        print*,'pb moldiff',istep,ig,l,nn,Rrafk(l,nn),tdiff,&
    534495     &  rho0*Rhoad(l),Zmin,Zmax,ntime
    535 !       print*,'Atri',Atri
    536 !       print*,'Btri',Btri
    537 !       print*,'Ctri',Ctri
    538 !       print*,'Dtri',Dtri
    539 !       print*,'Xtri',Xtri
    540 !       print*,'alpha',alpha
    541 !       print*,'beta',beta
    542 !       print*,'gama',gama
    543 !       print*,'delta',delta
    544 !       print*,'eps',eps
    545 !       print*,'Dad',Dad
     496        print*,'Atri',Atri
     497        print*,'Btri',Btri
     498        print*,'Ctri',Ctri
     499        print*,'Dtri',Dtri
     500        print*,'Xtri',Xtri
     501        print*,'alpha',alpha
     502        print*,'beta',beta
     503        print*,'gama',gama
     504        print*,'delta',delta
     505        print*,'eps',eps
     506        print*,'Dad',Dad
     507        print*,'Temp',Traf
     508        print*,'alt',Zraf
     509        print*,'Mraf',Mraf
     510        stop
    546511!       pdqdiff(1:ngridmx,1:nlayermx,1:nqmx)=0.
    547512!       return
    … …  
    593558        CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayermx,ncompdiff)
    594559
    595 !       if (ig .eq. ij0) then
    596 !       do l=il0,nlayermx
    597 !        write(*,'(i2,1x,19(e12.4,1x))') l,zz(l),tt(l),RHOK(l,16)/sum(RHOK(l,:)),RHOKINIT(l,16)/sum(RHOKINIT(l,:))
    598 !        enddo
    599 !       endif
     560        if (ig .eq. ij0) then
     561        do l=il0,nlayermx
     562        write(*,'(i2,1x,19(e12.4,1x))') l,zz(l),tt(l),RHOK(l,1)/sum(RHOK(l,:)),RHOKINIT(l,1)/sum(RHOKINIT(l,:)),&
     563     &  RHOK(l,5)/sum(RHOK(l,:)),RHOKINIT(l,5)/sum(RHOKINIT(l,:)),&
     564     &  RHOK(l,3)/sum(RHOK(l,:)),RHOKINIT(l,3)/sum(RHOKINIT(l,:)),&
     565     &  RHOK(l,6)/sum(RHOK(l,:)),RHOKINIT(l,6)/sum(RHOKINIT(l,:)),&
     566     &  RHOK(l,7)/sum(RHOK(l,:)),RHOKINIT(l,7)/sum(RHOKINIT(l,:))
     567        enddo
     568        endif
    600569
    601570! Compute total mass of each specie on the GCM vertical grid
    … …  
    640609
    641610       enddo  ! ig loop         
     611
     612        deallocate(qq,qnew,qint)
     613        deallocate(FacMass)
     614        deallocate(rhok,rhokinit)
     615        deallocate(wi,wad,uthermal,lambdaexo)
     616        deallocate(Hspecie)
     617        deallocate(Mtot1,Mtot2,Mraf1,Mraf2)
     618
    642619      return
    643620      end
    … …  
    12421219        INTEGER :: nl,l
    12431220        REAL*8,DIMENSION(nl) :: T,D
    1244         REAL*8 :: DZ
     1221        REAL*8 :: DZ,DZinv
    12451222        REAL*8,DIMENSION(nl) :: alpha,beta,gama,delta,eps
    12461223
    12471224! Compute alpha,beta and delta
    12481225! lower altitude values
    1249 
    1250         alpha(1)=1D0/T(1)*(-3D0*T(1)+4D0*T(2)-T(3))/(2D0*DZ)
    1251         delta(1)=(-3D0*D(1)+4D0*D(2)-D(3))/(2D0*DZ)
    1252         beta(1)=1D0/T(1)
    1253 
    1254         do l=2,nl-1
    1255         alpha(l)=1D0/T(l)*(T(l+1)-T(l-1))/(2D0*DZ)
    1256         delta(l)=(D(l+1)-D(l-1))/(2D0*DZ)
     1226        DZinv=1d0/(2D0*DZ)
     1227
     1228        beta(1)=1d0/T(1)
     1229        alpha(1)=beta(1)*(-3D0*T(1)+4D0*T(2)-T(3))*Dzinv
     1230        delta(1)=(-3D0*D(1)+4D0*D(2)-D(3))*Dzinv
     1231
     1232        beta(2)=1d0/T(2)
     1233        alpha(2)=beta(2)*(T(3)-T(1))*Dzinv
     1234        delta(2)=(D(3)-D(1))*Dzinv
     1235
     1236!       do l=2,nl-1
     1237!       beta(l)=1D0/T(l)
     1238!       alpha(l)=beta(l)*(T(l+1)-T(l-1))*Dzinv
     1239!       delta(l)=(D(l+1)-D(l-1))*Dzinv
     1240!       end do   
     1241
     1242        do l=3,nl-1
    12571243        beta(l)=1D0/T(l)
    1258         end do   
     1244        alpha(l)=beta(l)*(T(l+1)-T(l-1))*Dzinv
     1245        delta(l)=(D(l+1)-D(l-1))*Dzinv
     1246        gama(l-1)=(beta(l)-beta(l-2))*Dzinv
     1247        eps(l-1)=(alpha(l)-alpha(l-2))*Dzinv
     1248        enddo
    12591249
    12601250! Upper altitude values
    12611251
    1262         alpha(nl)=1D0/T(nl)*(3D0*T(nl)-4D0*T(nl-1)+T(nl-2))/(2D0*DZ)     
    1263         delta(nl)=(3D0*D(nl)-4D0*D(nl-1)+D(nl-2))/(2D0*DZ)
    12641252        beta(nl)=1D0/T(nl)
     1253        alpha(nl)=beta(nl)*(3D0*T(nl)-4D0*T(nl-1)+T(nl-2))*Dzinv     
     1254        delta(nl)=(3D0*D(nl)-4D0*D(nl-1)+D(nl-2))*Dzinv
    12651255
    12661256! Compute the gama and eps coefficients
    12671257! Lower altitude values
    12681258
    1269         gama(1)=(-3D0*beta(1)+4D0*beta(2)-beta(3))/(2d0*dz)
    1270         eps(1)=(-3D0*alpha(1)+4D0*alpha(2)-alpha(3))/(2d0*dz)
    1271 
    1272         do l=2,nl-1
    1273         gama(l)=(beta(l+1)-beta(l-1))/(2d0*Dz)
    1274         eps(l)=(alpha(l+1)-alpha(l-1))/(2d0*Dz)
    1275         end do
    1276 
    1277         gama(nl)=(3D0*beta(nl)-4D0*beta(nl-1)+beta(nl-2))/(2D0*DZ)
    1278         eps(nl)=(3D0*alpha(nl)-4D0*alpha(nl-1)+alpha(nl-2))/(2D0*DZ)
    1279 
     1259        gama(1)=(-3D0*beta(1)+4D0*beta(2)-beta(3))*Dzinv
     1260        eps(1)=(-3D0*alpha(1)+4D0*alpha(2)-alpha(3))*Dzinv
     1261
     1262        gama(nl-1)=(beta(nl)-beta(nl-2))*Dzinv
     1263        eps(nl-1)=(alpha(nl)-alpha(nl-2))*Dzinv
     1264
     1265!       do l=2,nl-1
     1266!       gama(l)=(beta(l+1)-beta(l-1))*Dzinv
     1267!       eps(l)=(alpha(l+1)-alpha(l-1))*Dzinv
     1268!       end do
     1269
     1270        gama(nl)=(3D0*beta(nl)-4D0*beta(nl-1)+beta(nl-2))*Dzinv
     1271        eps(nl)=(3D0*alpha(nl)-4D0*alpha(nl-1)+alpha(nl-2))*Dzinv
     1272
     1273!       do l=1,nl
     1274!       print*,'test diffparam',alpha(l),beta(l),delta(l),gama(l),eps(l)
     1275!       enddo
     1276!       stop
    12801277
    12811278        END
    … …  
    15411538
    15421539        END
     1540       
Note: See TracChangeset for help on using the changeset viewer.