Changeset 645


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

Location:
trunk/LMDZ.MARS
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r635 r645  
    16451645- Minor correction in "callsedim": compute "rdust" and/or "rice" only when
    16461646  it makes sense.
     1647
     1648== 03/05/12 == JYC+EM
     1649- Update of molecular diffusion routines (diffusion.h, moldiff_red.F90
     1650  and moldiffcoeff_red.F): optimized to run faster and diffuse only
     1651  specific tracers (if they are present).
     1652- bugs fixed in simpleclouds.F (extra arguments) and watercloud.F (wrong
     1653  arguments in call to simpleclouds and rdust() should only be recomputed
     1654  if dust_mass & dust_number tracers are available).

  • trunk/LMDZ.MARS/libf/aeronomars/diffusion.h

    r606 r645  
    55!**********************************************************************
    66
    7       integer ncompdiff
    8       integer nmajdiff,nmindiff
    9       integer idifflow
     7      real*8 Pdiff     
    108      real*8 tdiffmin
    119      real*8 dzres
    1210
    13 !      parameter (ncomptot=ncomp+1)     
    14       parameter (ncompdiff=16)
    15 !      parameter (nmajdiff=4)
    16 !      parameter (nmindiff=ncompdiff-nmajdiff)
    17       parameter (idifflow=20)
    18       parameter (tdiffmin=5d0)
    19       parameter (dzres=2d0)     ! grid resolution (km), for diffusion
     11      parameter (Pdiff=15.)      ! pressure below which diffusion is computed
     12      parameter (tdiffmin=10d0)
     13      parameter (dzres=2d0)      ! grid resolution (km) for diffusion
     14     

  • 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       

  • trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff_red.F

    r517 r645  
    1       subroutine moldiffcoeff_red(dij)
     1      subroutine moldiffcoeff_red(dij,indic,gcmind,ncompdiff2)
    22
    33       IMPLICIT NONE
    … …  
    2424c    ------------
    2525c       integer,parameter :: ncompmoldiff = 12
    26       real dij(ncompdiff,ncompdiff)
     26        integer ncompdiff2
     27        integer gcmind(ncompdiff2)
     28      real dij(ncompdiff2,ncompdiff2)
     29      integer indic(nqmx)
    2730
    2831c    Local variables:
    … …  
    3336cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    3437
    35       integer,parameter :: i_o   = 1
    36       integer,parameter :: i_n2   = 2
    37       integer,parameter :: i_co   = 3
    38       integer,parameter :: i_ar  = 4
    39       integer,parameter :: i_h2   = 5
    40       integer,parameter :: i_h    = 6
    41       integer,parameter :: i_o2   = 7
    42       integer,parameter :: i_oh  = 8
    43       integer,parameter :: i_ho2  = 9
    44       integer,parameter :: i_h2o = 10
    45       integer,parameter :: i_h2o2  = 11
    46       integer,parameter :: i_o1d   = 12
     38      real :: dijh2co,dijh2n2,dijh2co2,dijh2o2,dijho,dijref
     39!       integer :: i_h2,i_h,i_o
     40        integer :: g_h2,g_h,g_o
     41!      integer,parameter :: i_o   = 1
     42!      integer,parameter :: i_n2   = 2
     43!      integer,parameter :: i_co   = 3
     44!      integer,parameter :: i_ar  = 4
     45!      integer,parameter :: i_h2   = 5
     46!      integer,parameter :: i_h    = 6
     47!      integer,parameter :: i_o2   = 7
     48!      integer,parameter :: i_oh  = 8
     49!      integer,parameter :: i_ho2  = 9
     50!      integer,parameter :: i_h2o = 10
     51!      integer,parameter :: i_h2o2  = 11
     52!      integer,parameter :: i_o1d   = 12
    4753!      integer,parameter :: i_o3   = 13
    48       integer,parameter :: i_n    = 13
    49       integer,parameter :: i_no   = 14
    50       integer,parameter :: i_no2  = 15
     54!      integer,parameter :: i_n    = 13
     55!      integer,parameter :: i_no   = 14
     56!      integer,parameter :: i_no2  = 15
    5157!      integer,parameter :: i_n2d  = 17
    5258!      integer,parameter :: i_oplus = 18
    53       integer,parameter :: i_co2    = 16
     59!      integer,parameter :: i_co2    = 16
    5460!      integer,parameter :: i_oplus = 17
    5561!      integer,parameter :: i_hplus = 18
    5662
    5763! Tracer indexes in the GCM:
    58       integer,save :: g_co2=0
    59       integer,save :: g_co=0
    60       integer,save :: g_o=0
    61       integer,save :: g_o1d=0
    62       integer,save :: g_o2=0
    63       integer,save :: g_o3=0
    64       integer,save :: g_h=0
    65       integer,save :: g_h2=0
    66       integer,save :: g_oh=0
    67       integer,save :: g_ho2=0
    68       integer,save :: g_h2o2=0
    69       integer,save :: g_n2=0
    70       integer,save :: g_ar=0
    71       integer,save :: g_h2o=0
    72       integer,save :: g_n=0
    73       integer,save :: g_no=0
    74       integer,save :: g_no2=0
    75       integer,save :: g_n2d=0
     64!      integer,save :: g_co2=0
     65!      integer,save :: g_co=0
     66!      integer,save :: g_o=0
     67!      integer,save :: g_o1d=0
     68!      integer,save :: g_o2=0
     69!      integer,save :: g_o3=0
     70!      integer,save :: g_h=0
     71!      integer,save :: g_h2=0
     72!      integer,save :: g_oh=0
     73!      integer,save :: g_ho2=0
     74!      integer,save :: g_h2o2=0
     75!      integer,save :: g_n2=0
     76!      integer,save :: g_ar=0
     77!      integer,save :: g_h2o=0
     78!      integer,save :: g_n=0
     79!      integer,save :: g_no=0
     80!      integer,save :: g_no2=0
     81!      integer,save :: g_n2d=0
    7682!      integer,save :: g_oplus=0
    7783!      integer,save :: g_hplus=0
    7884
    79       integer,save :: gcmind(ncompdiff)
     85!      integer,save :: gcmind(ncompdiff)
    8086
    8187      real dnh
    … …  
    8591      if (firstcall) then
    8692        ! identify the indexes of the tracers we'll need
    87         g_co2=igcm_co2
    88         if (g_co2.eq.0) then
    89           write(*,*) "moldiffcoeff: Error; no CO2 tracer !!!"
    90           stop
    91         endif
    92         g_n2=igcm_n2
    93         if (g_n2.eq.0) then
    94           write(*,*) "moldiffcoeff: Error; no N2 tracer !!!"
    95           stop
    96         endif
    97         g_ar=igcm_ar
    98         if (g_ar.eq.0) then
    99           write(*,*) "moldiffcoeff: Error; no Ar tracer !!!"
    100           stop
    101         endif       
    102         g_h2=igcm_h2
    103         if (g_h2.eq.0) then
    104           write(*,*) "moldiffcoeff: Error; no H2 tracer !!!"
    105           stop
    106         endif
    107         g_h=igcm_h
    108         if (g_h.eq.0) then
    109           write(*,*) "moldiffcoeff: Error; no H tracer !!!"
    110           stop
    111         endif
    112         g_co=igcm_co
    113         if (g_co.eq.0) then
    114           write(*,*) "moldiffcoeff: Error; no CO tracer !!!"
    115           stop
    116         endif
    117         g_o2=igcm_o2
    118         if (g_o2.eq.0) then
    119           write(*,*) "moldiffcoeff: Error; no O2 tracer !!!"
    120           stop
    121         endif
    122         g_oh=igcm_oh
    123         if (g_oh.eq.0) then
    124           write(*,*) "moldiffcoeff: Error; no OH tracer !!!"
    125           stop
    126         endif
    127         g_ho2=igcm_ho2
    128         if (g_ho2.eq.0) then
    129           write(*,*) "moldiffcoeff: Error; no HO2 tracer !!!"
    130           stop
    131         endif
    132         g_h2o=igcm_h2o_vap
    133         if (g_h2o.eq.0) then
    134           write(*,*) "moldiffcoeff: Error; no H2O tracer !!!"
    135           stop
    136         endif
    137         g_h2o2=igcm_h2o2
    138         if (g_h2o2.eq.0) then
    139           write(*,*) "moldiffcoeff: Error; no H2O2 tracer !!!"
    140           stop
    141         endif
    142         g_o1d=igcm_o1d
    143         if (g_h.eq.0) then
    144           write(*,*) "moldiffcoeff: Error; no O1d tracer !!!"
    145           stop
    146         endif
    147         g_o3=igcm_o3
    148         if (g_o3.eq.0) then
    149           write(*,*) "moldiffcoeff: Error; no O3 tracer !!!"
    150           stop
    151         endif
    152         g_n=igcm_n
    153         if (g_n.eq.0) then
    154           write(*,*) "moldiffcoeff: Error; no N tracer !!!"
    155           stop
    156         endif
    157         g_no=igcm_no
    158         if (g_no.eq.0) then
    159           write(*,*) "moldiffcoeff: Error; no NO tracer !!!"
    160           stop
    161         endif
    162         g_no2=igcm_no2
    163         if (g_no2.eq.0) then
    164           write(*,*) "moldiffcoeff: Error; no NO2 tracer !!!"
    165           stop
    166         endif
    167         g_n2d=igcm_n2d
    168         if (g_n2d.eq.0) then
    169           write(*,*) "moldiffcoeff: Error; no N2(D) tracer !!!"
    170           stop
    171         endif
     93!        g_co2=igcm_co2
     94!        if (g_co2.eq.0) then
     95!          write(*,*) "moldiffcoeff: Error; no CO2 tracer !!!"
     96!          stop
     97!        endif
     98!        g_n2=igcm_n2
     99!        if (g_n2.eq.0) then
     100!          write(*,*) "moldiffcoeff: Error; no N2 tracer !!!"
     101!          stop
     102!        endif
     103!        g_ar=igcm_ar
     104!        if (g_ar.eq.0) then
     105!          write(*,*) "moldiffcoeff: Error; no Ar tracer !!!"
     106!          stop
     107!        endif       
     108!        g_h2=igcm_h2
     109!        if (g_h2.eq.0) then
     110!          write(*,*) "moldiffcoeff: Error; no H2 tracer !!!"
     111!          stop
     112!        endif
     113!        g_h=igcm_h
     114!        if (g_h.eq.0) then
     115!          write(*,*) "moldiffcoeff: Error; no H tracer !!!"
     116!          stop
     117!        endif
     118!        g_co=igcm_co
     119!        if (g_co.eq.0) then
     120!          write(*,*) "moldiffcoeff: Error; no CO tracer !!!"
     121!          stop
     122!        endif
     123!        g_o2=igcm_o2
     124!        if (g_o2.eq.0) then
     125!          write(*,*) "moldiffcoeff: Error; no O2 tracer !!!"
     126!          stop
     127!        endif
     128!        g_oh=igcm_oh
     129!        if (g_oh.eq.0) then
     130!          write(*,*) "moldiffcoeff: Error; no OH tracer !!!"
     131!          stop
     132!        endif
     133!        g_ho2=igcm_ho2
     134!        if (g_ho2.eq.0) then
     135!          write(*,*) "moldiffcoeff: Error; no HO2 tracer !!!"
     136!          stop
     137!        endif
     138!        g_h2o=igcm_h2o_vap
     139!        if (g_h2o.eq.0) then
     140!          write(*,*) "moldiffcoeff: Error; no H2O tracer !!!"
     141!          stop
     142!        endif
     143!        g_h2o2=igcm_h2o2
     144!        if (g_h2o2.eq.0) then
     145!          write(*,*) "moldiffcoeff: Error; no H2O2 tracer !!!"
     146!          stop
     147!        endif
     148!        g_o1d=igcm_o1d
     149!        if (g_h.eq.0) then
     150!          write(*,*) "moldiffcoeff: Error; no O1d tracer !!!"
     151!          stop
     152!        endif
     153!        g_o3=igcm_o3
     154!        if (g_o3.eq.0) then
     155!          write(*,*) "moldiffcoeff: Error; no O3 tracer !!!"
     156!          stop
     157!        endif
     158!        g_n=igcm_n
     159!        if (g_n.eq.0) then
     160!          write(*,*) "moldiffcoeff: Error; no N tracer !!!"
     161!          stop
     162!        endif
     163!        g_no=igcm_no
     164!        if (g_no.eq.0) then
     165!          write(*,*) "moldiffcoeff: Error; no NO tracer !!!"
     166!          stop
     167!        endif
     168!        g_no2=igcm_no2
     169!        if (g_no2.eq.0) then
     170!          write(*,*) "moldiffcoeff: Error; no NO2 tracer !!!"
     171!          stop
     172!        endif
     173!        g_n2d=igcm_n2d
     174!        if (g_n2d.eq.0) then
     175!          write(*,*) "moldiffcoeff: Error; no N2(D) tracer !!!"
     176!          stop
     177!        endif
    172178!        g_oplus=igcm_oplus
    173179!        if (g_oplus .eq. 0) then
    … …  
    180186!        stop
    181187!        endif
    182         g_o=igcm_o
    183         if (g_o.eq.0) then
    184           write(*,*) "moldiffcoeff: Error; no O tracer !!!"
    185           stop
    186         endif
     188!        g_o=igcm_o
     189!        if (g_o.eq.0) then
     190!          write(*,*) "moldiffcoeff: Error; no O tracer !!!"
     191!          stop
     192!        endif
    187193       
    188194c
    … …  
    191197cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    192198
    193         gcmind(i_co2)  =   g_co2
    194         gcmind(i_n2)  =   g_n2
    195         gcmind(i_ar)  =   g_ar
    196         gcmind(i_h2) =   g_h2
    197         gcmind(i_h)  =   g_h
    198         gcmind(i_co)   =   g_co
    199         gcmind(i_o2) =   g_o2
    200         gcmind(i_oh)=   g_oh
    201         gcmind(i_ho2)  =   g_ho2
    202         gcmind(i_h2o) = g_h2o
    203         gcmind(i_h2o2)= g_h2o2
    204         gcmind(i_o1d) = g_o1d
     199!        gcmind(i_co2)  =   g_co2
     200!        gcmind(i_n2)  =   g_n2
     201!        gcmind(i_ar)  =   g_ar
     202!        gcmind(i_h2) =   g_h2
     203!        gcmind(i_h)  =   g_h
     204!        gcmind(i_co)   =   g_co
     205!        gcmind(i_o2) =   g_o2
     206!        gcmind(i_oh)=   g_oh
     207!        gcmind(i_ho2)  =   g_ho2
     208!        gcmind(i_h2o) = g_h2o
     209!        gcmind(i_h2o2)= g_h2o2
     210!        gcmind(i_o1d) = g_o1d
    205211!        gcmind(i_o3) = g_o3
    206         gcmind(i_n)= g_n
    207         gcmind(i_no) = g_no
    208         gcmind(i_no2) = g_no2
     212!        gcmind(i_n)= g_n
     213!        gcmind(i_no) = g_no
     214!        gcmind(i_no2) = g_no2
    209215!        gcmind(i_n2d) = g_n2d
    210216!        gcmind(i_oplus) =  g_oplus
    211217!        gcmind(i_hplus) = g_hplus
    212         gcmind(i_o)   =   g_o
     218!        gcmind(i_o)   =   g_o
    213219
    214220c
    … …  
    217223      endif ! of if (firstcall)
    218224
    219 
    220       dij(i_h2,i_co)   = 0.0000651
    221       dij(i_h2,i_n2)   = 0.0000674
    222       dij(i_h2,i_o2)   = 0.0000697
    223       dij(i_h2,i_co2)  = 0.0000550
    224       dij(i_h2,i_h2)   = 0.0
    225       dij(i_h2,i_h)    = 0.0
    226       dij(i_h2,i_h2o)  = 0.0    !0003
    227       dij(i_h2,i_h2o2) = 0.0    !0003
     225        dijh2co = 0.0000651
     226        dijh2n2 = 0.0000674
     227        dijh2o2 = 0.0000697
     228        dijh2co2 = 0.0000550
     229        dijho = 0.000114
     230
     231!      dij(i_h2,i_co)   = 0.0000651
     232!      dij(i_h2,i_n2)   = 0.0000674
     233!      dij(i_h2,i_o2)   = 0.0000697
     234!      dij(i_h2,i_co2)  = 0.0000550
     235!      dij(i_h2,i_h2)   = 0.0
     236!      dij(i_h2,i_h)    = 0.0
     237!      dij(i_h2,i_h2o)  = 0.0   !0003
     238!      dij(i_h2,i_h2o2) = 0.0   !0003
    228239!      dij(i_h2,i_o3)   = 0.0   !0003
    229       dij(i_h2,i_o)    = 0.0
    230       dij(i_h2,i_ar)   = 0.0
    231       dij(i_h2,i_n)    = 0.0
    232 
    233 c      dij(i_h,i_o)     = 0.0000144
    234       dij(i_h,i_o)     = 0.000114
     240!      dij(i_h2,i_o)    = 0.0
     241!      dij(i_h2,i_ar)   = 0.0
     242!      dij(i_h2,i_n)    = 0.0
     243
     244!c      dij(i_h,i_o)     = 0.0000144
     245!      dij(i_h,i_o)     = 0.000114
     246
     247! find h2, h and o index in gcm
     248! these species are used to define the diffusion coefficients
     249
     250        do n=1,nqmx
     251        if (noms(n) .eq. 'h2') g_h2=n
     252        if (noms(n) .eq. 'h') g_h=n
     253        if (noms(n) .eq. 'o') g_o=n
     254        enddo
     255        print*,'gh2',g_h2,g_h,g_o
    235256
    236257       print*,'moldiffcoef: COEFF CALC'
    237258       open(56,file='coeffs.dat',status='unknown')
    238       do n=1,ncompdiff
    239         if (dij(i_h2,n).gt.0.0) then
    240           do nn=n,ncompdiff
    241             dij(nn,n)=dij(i_h2,n)
     259      do n=1,ncompdiff2
     260        dijref=0.
     261        if (noms(gcmind(n)) .eq. 'co') dijref=dijh2co
     262        if (noms(gcmind(n)) .eq. 'n2') dijref=dijh2n2
     263        if (noms(gcmind(n)) .eq. 'o2') dijref=dijh2o2
     264        if (noms(gcmind(n)) .eq. 'co2') dijref=dijh2co2
     265!       print*,'test',n,dijref
     266        if (dijref .gt. 0.0) then
     267          do nn=n,ncompdiff2
     268            dij(nn,n)=dijref
    242269     &                  *sqrt(mmol(g_h2)/mmol(gcmind(nn)))
    243270            if(n.eq.nn) dij(nn,n)=1.0
    … …  
    245272          enddo
    246273        endif
    247         if (dij(i_h2,n).eq.0.0) then
    248           dnh=dij(i_h,i_o)*sqrt(mmol(g_o)/mmol(gcmind(n)))
    249           do nn=n,ncompdiff
     274        if (dijref .eq. 0.0) then
     275        dijref=dijho
     276          dnh=dijref*sqrt(mmol(g_o)/mmol(gcmind(n)))
     277          do nn=n,ncompdiff2
    250278            dij(nn,n)=dnh*sqrt(mmol(g_h)/mmol(gcmind(nn)))
    251279            if(n.eq.nn) dij(nn,n)=1.0
    … …  
    255283      enddo
    256284
    257       do n=1,ncompdiff
    258         do nn=n,ncompdiff
     285      do n=1,ncompdiff2
     286        do nn=n,ncompdiff2
    259287          write(56,*) n,nn,dij(n,nn)    !*1.e5/1.381e-23/(273**1.75)
    260288        enddo
    … …  
    265293      return   
    266294      end
     295     

  • trunk/LMDZ.MARS/libf/phymars/simpleclouds.F

    r633 r645  
    11      subroutine simpleclouds(ngrid,nlay,ptimestep,
    2      &             pplev,pplay,pzlay,pt,pdt,
     2     &             pplay,pzlay,pt,pdt,
    33     &             pq,pdq,pdqcloud,pdtcloud,
    44     &             nq,tau,rice)
    … …  
    4141      integer nq                 ! nombre de traceurs
    4242      REAL ptimestep             ! pas de temps physique (s)
    43       REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
     43!      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
    4444      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
    4545      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers

  • trunk/LMDZ.MARS/libf/phymars/watercloud.F

    r633 r645  
    210210     &             pplay,pzlay,pt,subpdt,
    211211     &             pq,subpdq,subpdqcloud,subpdtcloud,
    212      &             nq,tau)
     212     &             nq,tau,rice)
    213213        ENDIF
    214214
    … …  
    360360        ENDDO
    361361       ELSE
    362         DO l=1,nlay
     362       
     363        if ((igcm_dust_mass.ne.0).and.(igcm_dust_number.ne.0)) then
     364         ! recompute rdust(), if we have dust_mass & dust_number tracers
     365         DO l=1,nlay
    363366          DO ig=1,ngrid
    364 
    365           rdust(ig,l)=
     367           rdust(ig,l)=
    366368     &      CBRT(r3n_q*(pq(ig,l,igcm_dust_mass)+
    367369     &           (pdq(ig,l,igcm_dust_mass)
    … …  
    370372     &           (pdq(ig,l,igcm_dust_number)+
    371373     &            pdqcloud(ig,l,igcm_dust_number)*ptimestep)))
    372           rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
    373 
     374           rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
     375          ENDDO
     376         ENDDO
     377        endif ! of if ((igcm_dust_mass.ne.0).and.(igcm_dust_number.ne.0))
     378       
     379        DO l=1,nlay
     380          DO ig=1,ngrid
    374381            ccntyp =
    375382     &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-pzlay(ig,l)/10000.)
Note: See TracChangeset for help on using the changeset viewer.