Changeset 633


Ignore:
Timestamp:
Apr 25, 2012, 5:38:49 PM (13 years ago)
Author:
tnavarro
Message:

last scheme in commit r626 led to a wrong physical behaviour. This version uses a new subtimestep for microphysics that should be faster than the last one.

Location:
trunk/LMDZ.MARS
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r627 r633  
    16221622>> Some cleanup on messages and comments in code about the reference pressure
    16231623   for dust opacity which is now 610 Pa.
     1624
     1625== 25/04/12 == TN
     1626>> The new scheme does not work for now. Instead, use of a subtimestep for nucleation and growth.

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

    r626 r633  
    1       subroutine growthrate(temp,pmid,psat,seq,rcrystal,coeff1,coeff2)
     1      subroutine growthrate(temp,pmid,psat,rcrystal,res)
    22
    33      IMPLICIT NONE
    … …  
    3131      REAL temp     ! temperature in the middle of the layer (K)
    3232      REAL pmid     ! pressure in the middle of the layer (K)
    33       REAL*8 psat   ! water vapor saturation pressure (Pa)
    34       REAL*8 seq    ! Equilibrium saturation ratio
     33      REAL psat   ! water vapor saturation pressure (Pa)
    3534      REAL rcrystal ! crystal radius before condensation (m)
    3635
    3736c     Output
    38       REAL coeff1,coeff2     ! coefficients for the analytical relation between time and radius
     37      REAL res      ! growth resistance (res=Rk+Rd)
    3938
    4039
    … …  
    4241c   ------
    4342
    44       REAL*8 k,Lv                 
    45       REAL*8 knudsen           ! Knudsen number (gas mean free path/particle radius)
    46       REAL*8 afactor,Dv,lambda       ! Intermediate computations for growth rate
    47       REAL*8 Rk,Rd
     43      REAL k,Lv                 
     44      REAL knudsen           ! Knudsen number (gas mean free path/particle radius)
     45      REAL afactor,Dv,lambda       ! Intermediate computations for growth rate
     46      REAL Rk,Rd
    4847     
    4948     
    … …  
    8685
    8786c     - Compute Rk
    88       Rk = Lv**2 * rho_ice * mh2o / (k*rgp*temp*temp)
     87      Rk = Lv*Lv* rho_ice * mh2o / (k*rgp*temp*temp)
    8988c     - Compute Rd
    9089      Rd = rgp * temp *rho_ice / (Dv*psat*mh2o)
    9190     
    9291     
    93       coeff1 = real(Rk + Rd*(1. + lambda * knudsen))
    94       coeff2 = real(Rk + Rd*(1. - lambda * knudsen))
     92      res = Rk + Rd*(1. + lambda * knudsen)
     93     
     94      !coeff1 = real(Rk + Rd*(1. + lambda * knudsen))
     95      !coeff2 = real(Rk + Rd*(1. - lambda * knudsen))
    9596     
    9697c Below are growth rate used for other schemes

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

    r626 r633  
    11      subroutine improvedclouds(ngrid,nlay,ptimestep,
    2      &             pplev,pplay,pt,pdt,
     2     &             pplay,pt,pdt,
    33     &             pq,pdq,pdqcloud,pdtcloud,
    4      &             nq,tauscaling,rdust,rice,nuice,
    5      &             rsedcloud,rhocloud)
     4     &             nq,tauscaling)
    65! to use  'getin'
    76      USE ioipsl_getincom
    87      implicit none
     8     
     9     
    910c------------------------------------------------------------------
    1011c  This routine is used to form clouds when a parcel of the GCM is
    … …  
    2122c    values which are then used by the sedimentation and advection
    2223c    schemes.
    23 c  A word about the radius growth ...
    24 c  A word about nucleation and ice growth strategies ...
    2524
    2625c  Authors: J.-B. Madeleine, based on the work by Franck Montmessin
    … …  
    4342      integer nq                 ! nombre de traceurs
    4443      REAL ptimestep             ! pas de temps physique (s)
    45       REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
    4644      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
    4745           
    … …  
    5452                                 !   (kg/kg.s-1)
    5553      REAL tauscaling(ngridmx)     ! Convertion factor for qdust and Ndust
    56       REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
    5754
    5855c     Outputs:
    59       REAL rice(ngrid,nlay)      ! Ice mass mean radius (m)
    60                                  ! (r_c in montmessin_2004)
    61       REAL nuice(ngrid,nlay)     ! Estimated effective variance
    62                                  !   of the size distribution
    63       REAL rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
    64       REAL rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
    6556      REAL pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
    6657                                   !   H2O(kg/kg.s-1)
    … …  
    8980      REAL zqsat(ngridmx,nlayermx)    ! saturation
    9081      REAL lw                         !Latent heat of sublimation (J.kg-1)
    91       REAL Cste
    92       REAL dMice           ! mass of condensated ice
    93       REAL sumcheck
     82      REAL cste
     83      REAL dMice           ! mass of condensed ice
     84!      REAL sumcheck
    9485      REAL*8 ph2o          ! Water vapor partial pressure (Pa)
    9586      REAL*8 satu          ! Water vapor saturation ratio over ice
    9687      REAL*8 Mo,No
    97       REAL*8 dN,dM,newvap
    98       REAL*8 Rn, Rm, dev2, yeah, n_derf, m_derf
     88      REAL*8 Rn, Rm, dev2, n_derf, m_derf
    9989      REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
    10090      REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
    101       REAL*8 rate(nbin_cld)  ! nucleation rate
    102       !REAL*8 up,dwn,Ctot,gr
    103       REAl*8 seq
     91
    10492      REAL*8 sig      ! Water-ice/air surface tension  (N.m)
    10593      EXTERNAL sig
     94
     95      REAL dN,dM
     96      REAL rate(nbin_cld)  ! nucleation rate
     97      REAL seq
     98
     99      REAL rice(ngrid,nlay)      ! Ice mass mean radius (m)
     100                                 ! (r_c in montmessin_2004)
     101      REAL rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
     102      REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
     103
     104      REAL res      ! Resistance growth
     105
    106106
    107107c     Parameters of the size discretization
    … …  
    118118      DOUBLE PRECISION vol_cld(nbin_cld)  ! particle volume for each bin (m3)
    119119
     120
    120121      REAL sigma_ice ! Variance of the ice and CCN distributions
    121122      SAVE sigma_ice
    122      
    123       REAL tdicho, tmax, rmin, rmax, req, rdicho
    124       REAL coeff0, coeff1, coeff2
    125       REAL error_out(ngridmx,nlayermx)
    126       REAL error2d(ngridmx)
    127      
    128       REAL var1,var2,var3
    129      
    130       REAL rccn, epsilon
     123
     124
    131125     
    132126c----------------------------------     
    133 c some outputs for 1D -- TESTS
    134       REAL satu_out(ngridmx,nlayermx) ! satu ratio for output
    135       REAL dN_out(ngridmx,nlayermx)  ! mass variation for output
    136       REAL dM_out(ngridmx,nlayermx)  ! number variation for output
    137       REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)
    138       REAL gr_out(ngridmx,nlayermx)   ! for 1d output     
    139       REAL rice_out(ngridmx,nlayermx) ! ice radius change
    140       REAL rate_out(ngridmx,nlayermx) ! nucleation rate
    141       REAL satubf(ngridmx,nlayermx),satuaf(ngridmx,nlayermx)
    142       REAL ccnbf(ngridmx,nlayermx),ccnaf(ngridmx,nlayermx)
    143       INTEGER count
     127c TESTS
     128
     129      INTEGER countcells
    144130     
    145131      LOGICAL test_flag    ! flag for test/debuging outputs
    146      
    147       test_flag = .false.
    148 c----------------------------------     
    149 c----------------------------------     
     132      SAVE    test_flag     
     133 
    150134
    151135c------------------------------------------------------------------
    … …  
    164148
    165149c       Volume ratio between two adjacent bins
     150   !     vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
     151   !     vrat_cld = exp(vrat_cld)
    166152        vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
    167153        vrat_cld = dexp(vrat_cld)
    … …  
    171157        rad_cld(1) = rmin_cld
    172158        vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
     159   !     vol_cld(1) = 4./3. * pi * rmin_cld*rmin_cld*rmin_cld
    173160
    174161        do i=1,nbin_cld-1
    … …  
    197184
    198185        do i=1,nbin_cld+1
     186    !        rb_cld(i) = log(rb_cld(i)) 
    199187            rb_cld(i) = dlog(rb_cld(i))  !! we save that so that it is not computed
    200188                                         !! at each timestep and gridpoint
    … …  
    210198c       Variance of the ice and CCN distributions
    211199        sigma_ice = sqrt(log(1.+nuice_sed))
    212        
     200       
     201 
    213202        write(*,*) 'Variance of ice & CCN distribs :', sigma_ice
    214203        write(*,*) 'nuice for sedimentation:', nuice_sed
    215204        write(*,*) 'Volume of a water molecule:', vo1
    216205
     206
     207        test_flag = .false.
     208 
    217209        firstcall=.false.
    218210      END IF
    219211
     212
    220213!=============================================================
    221214! 1. Initialisation
    222215!=============================================================
     216      cste = 4*pi*rho_ice*ptimestep
    223217
    224218c     Initialize the tendencies
    225219      pdqcloud(1:ngrid,1:nlay,1:nq)=0
    226220      pdtcloud(1:ngrid,1:nlay)=0
    227      
    228 c     Update the needed variables
    229       do l=1,nlay
    230         do ig=1,ngrid
    231 c         Temperature
    232           zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
    233 c         Dust mass mixing ratio
    234           zq(ig,l,igcm_dust_mass) =
    235      &      pq(ig,l,igcm_dust_mass) +
    236      &      pdq(ig,l,igcm_dust_mass) * ptimestep
    237           zq0(ig,l,igcm_dust_mass)=zq(ig,l,igcm_dust_mass)
    238 c         Dust particle number
    239           zq(ig,l,igcm_dust_number) =
    240      &      pq(ig,l,igcm_dust_number) +
    241      &      pdq(ig,l,igcm_dust_number) * ptimestep
    242           zq0(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)
    243 c         Update rdust from last tendencies
    244           rdust(ig,l)=
    245      &      CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
    246      &      max(zq(ig,l,igcm_dust_number),0.01))
    247           rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)         
    248 c         CCN mass mixing ratio
    249           zq(ig,l,igcm_ccn_mass)=
    250      &      pq(ig,l,igcm_ccn_mass)+pdq(ig,l,igcm_ccn_mass)*ptimestep
    251           zq(ig,l,igcm_ccn_mass)=max(zq(ig,l,igcm_ccn_mass),1.E-30)
    252           zq0(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
    253 c         CCN particle number
    254           zq(ig,l,igcm_ccn_number)=
    255      &      pq(ig,l,igcm_ccn_number)+pdq(ig,l,igcm_ccn_number)*ptimestep
    256           zq(ig,l,igcm_ccn_number)=max(zq(ig,l,igcm_ccn_number),1.E-30)
    257           zq0(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
    258 c         Water vapor
    259           zq(ig,l,igcm_h2o_vap)=
    260      &      pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep
    261           zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004
    262           zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)
    263 c         Water ice
    264           zq(ig,l,igcm_h2o_ice)=
    265      &      pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
    266           zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),1e-30) ! FF 12/2004
    267           zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
    268         enddo
    269       enddo
    270 
     221   
     222c     Initialize the tendencies
     223      pdqcloud(1:ngrid,1:nlay,1:nq)=0
     224      pdtcloud(1:ngrid,1:nlay)=0
     225     
     226c     Initialize the tendencies
     227      pdqcloud(1:ngrid,1:nlay,1:nq)=0
     228      pdtcloud(1:ngrid,1:nlay)=0
     229     
     230     
     231      zt(1:ngrid,1:nlay) =
     232     &      pt(1:ngrid,1:nlay) +
     233     &      pdt(1:ngrid,1:nlay) * ptimestep
     234
     235      zq(1:ngrid,1:nlay,1:nq) =
     236     &      pq(1:ngrid,1:nlay,1:nq) +
     237     &      pdq(1:ngrid,1:nlay,1:nq) * ptimestep
     238     
     239      zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
     240     
     241      WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
     242     &       zq(1:ngrid,1:nlay,1:nq) = 1.e-30
     243     
    271244!=============================================================
    272245! 2. Compute saturation
    273246!=============================================================
    274247
     248
    275249      dev2 = 1. / ( sqrt(2.) * sigma_ice )
    276       error_out(:,:) = 0.
    277 
    278250
    279251      call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
    280252           
    281       count = 0
     253      countcells = 0
    282254
    283255c     Main loop over the GCM's grid
    … …  
    288260        ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l)
    289261        satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
    290                
    291         IF (( satu .ge. 1. )       .or.                             ! if there is condensation
    292      &      ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.)) THEN ! or sublimation
    293 
    294262
    295263!=============================================================
    … …  
    297265!=============================================================
    298266
     267       IF ( satu .ge. 1. ) THEN         ! if there is condensation
     268
     269c         Update rdust from last tendencies
     270        rdust(ig,l)=
     271     &     CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
     272     &      (zq(ig,l,igcm_dust_number)+1.e-30))
     273        rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
     274
     275
    299276c       Expand the dust moments into a binned distribution
    300         Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)   !+ 1.e-30
     277        Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)   + 1.e-30
    301278        No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
    302279        Rn = rdust(ig,l)
    … …  
    313290          m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
    314291        enddo
    315 !!! MORE EFFICIENT COMPUTATIONALLY THAN
    316 !        Rn = rdust(ig,l)
    317 !        Rm = Rn * exp( 3. * sigma_ice**2. )
    318 !        Rn = 1. / Rn
    319 !        Rm = 1. / Rm
    320 !        do i = 1, nbin_cld
    321 !          n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)*Rn) * dev2 )
    322 !     &                           -derf( dlog(rb_cld(i) * Rn) * dev2 ) )
    323 !          m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)*Rm) * dev2 )
    324 !     &                           -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
    325 !        enddo
    326 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    327292       
    328293!        sumcheck = 0
    … …  
    349314!          print*, "Dust binned distribution", m_aer
    350315!        endif
    351                
     316
     317 
    352318c       Get the rates of nucleation
    353319        call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
    354320       
    355 
    356321        dN = 0.
    357322        dM = 0.
    358         rate_out(ig,l) = 0
    359323        do i = 1, nbin_cld
    360           n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
    361           m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
     324          n_aer(i) = n_aer(i)/( 1. + rate(i)*ptimestep)
     325          m_aer(i) = m_aer(i)/( 1. + rate(i)*ptimestep)
    362326          dN       = dN + n_aer(i) * rate(i) * ptimestep
    363327          dM       = dM + m_aer(i) * rate(i) * ptimestep
    364           !rate_out(ig,l)=rate_out(ig,l)+rate(i)
    365328        enddo
    366329
    367 c       Update CCNs, can also be done after the radius growth
    368 c       Dust particles
     330
     331c       Update Dust particles
    369332        zq(ig,l,igcm_dust_mass)   =
    370333     &         zq(ig,l,igcm_dust_mass)   - dM/ tauscaling(ig)
    371334        zq(ig,l,igcm_dust_number) =
    372335     &         zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig)
    373 c       CCNs
     336c       Update CCNs
    374337        zq(ig,l,igcm_ccn_mass)   =
    375338     &         zq(ig,l,igcm_ccn_mass)   + dM/ tauscaling(ig)
    … …  
    377340     &         zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig)
    378341       
     342        ENDIF ! of is satu >1
    379343
    380344!=============================================================
    … …  
    382346!=============================================================
    383347
     348c We trigger crystal growth if and only if there is at least one nuclei (N>1).
     349c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
     350c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
     351
     352       IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
     353
     354
    384355        Mo   = zq(ig,l,igcm_h2o_ice) +
    385      &           zq(ig,l,igcm_ccn_mass)* tauscaling(ig)  + 1.e-30
    386         No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
     356     &           zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
     357        No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
     358
    387359
    388360        rhocloud(ig,l) =  zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
    … …  
    391363        rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
    392364         
    393 c       nuclei radius 
    394         rccn  = CBRT(zq(ig,l,igcm_ccn_mass)/
    395      &                (pi*rho_dust*zq(ig,l,igcm_ccn_number)*4./3.))
    396365
    397366c       ice crystal radius   
    398367        rice (ig,l) =
    399      &      CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) )
    400 
    401 c       enforce physical value : crystal cannot be smaller than its nuclei !
    402         rice(ig,l) = max(rice(ig,l), rccn)   
     368     &      CBRT( real(Mo/No) * 0.75 / pi / rhocloud(ig,l) )
     369
    403370
    404371c       saturation at equilibrium
    … …  
    406373     &           (rho_ice*rgp*zt(ig,l)*rice(ig,l)) )
    407374
    408 
    409 c       If there is more than on nuclei, we peform ice growth 
    410         var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN
    411         IF (var1 .ge. -1) THEN
    412 
    413 
    414       if (test_flag) then
    415        print*, ' '
    416        print*, ptimestep
    417        print*, 'satu,seq', real(satu), real(seq), ig,l
    418        print*, 'dN,dM', real(dN),real(dM)           
    419        print*,'rccn', rccn
    420        print*, 'Nccn, Mccn', zq(ig,l,igcm_ccn_number)*tauscaling(ig),
    421      &  zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
    422       endif
    423      
    424 c         crystal radius to reach saturation at equilibrium (i.e. satu=seq)         
    425           req = ( zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap)
    426      &          + zq(ig,l,igcm_ccn_mass)*tauscaling(ig)*
    427      &           rho_ice/rho_dust - seq * zqsat(ig,l))
    428      &         / ( zq(ig,l,igcm_ccn_number)*tauscaling(ig)*
    429      &             pi*rho_ice*4./3. )
    430           req = CBRT(req)
    431          
    432 c         compute ice radius growth resistances (diffusive and latent heat resistancea)       
    433           call growthrate(zt(ig,l),pplay(ig,l),
    434      &                    ph2o/satu,seq,req,coeff1,coeff2)
    435                
    436           coeff0 = -zqsat(ig,l) / (4.*pi*req*rho_ice
    437      &               * zq(ig,l,igcm_ccn_number)*tauscaling(ig))
    438                
    439 c         compute tmax, the time needed to reach req           
    440           call phi(rice(ig,l),req,coeff1,coeff2,tmax)
    441          
    442       if (test_flag) then
    443           print*, 'coeffs', coeff0,coeff1,coeff2
    444           print*, 'req,tmax', req,tmax*coeff0
    445           print*, 'i,rmin,rdicho,rmax,tdicho'
    446       endif
    447 
    448 c         rmin is rice if r increases (satu >1) or req if it decreases (satu<1)
    449 c         if req is lower than rccn (ie not enough ice to reach saturation), rmin is forced to rccn
    450           if (satu .ge. seq) then
    451             ! crystal size is increasing
    452             rmin = max(min(rice(ig,l),req),rccn)
    453             rmax = max(rice(ig,l),req)
    454           else
    455             !  crystal size is decreasing
    456             rmax = max(min(rice(ig,l),req),rccn)
    457             rmin = max(rice(ig,l),req)
    458           endif
    459                          !rmax = min(rmax,1.e-3) ! au max on a des rayons de 1 mm  pour la dicho ...
    460           rdicho = 0.5*(rmin+rmax)
    461          
    462           ! for output
    463           var1 = rice(ig,l)
    464           var2 = rmin
    465           var3 = rmax
    466          
    467 c         Given the phi function is monotonous, we perform a simple dichotomy to find the raidus at t+1       
    468           do i = 1,10 ! dichotomy loop
    469              
    470 c            compute tdicho, the time needed to reach rdicho
    471              call phi(rdicho,req,coeff1,coeff2,tdicho)
    472              !print*, tdicho,tmax
    473              tdicho = coeff0*(tdicho - tmax)
    474              
    475              if (test_flag) print*, i,rmin,rdicho,rmax,tdicho
    476                        
    477              if (tdicho .ge. ptimestep) then
    478                 rmax = rdicho
    479              else
    480                 rmin = rdicho
    481              endif
    482              
    483              rdicho = 0.5*(rmin+rmax)
    484          
    485           enddo  ! of dichotomy loop
    486          
    487        if (test_flag) then
    488         if (abs(rdicho - rccn) .ge. 1e-15) then ! to avoid infinite values
    489           epsilon = (rmax - rmin)/(2**10)
    490           error_out(ig,l) =
    491      &    100*(epsilon**3 +3*epsilon**2*rdicho +3*epsilon*rdicho**2)
    492      &    / (rdicho**3-rccn**3)     
    493         endif
    494         print*, 'error masse glace %', error_out(ig,l)
    495         print*, 'rice,ice,vap bf',
    496      &   rice(ig,l),zq0(ig,l,igcm_h2o_ice),zq0(ig,l,igcm_h2o_vap)
    497        endif
    498          
    499 c         If the initial condition is subsaturated and there is not enough ice available for sublimation
    500 c         to reach equilibrium, req is neagtive. Therefore, enforce physical value.
    501           rice(ig,l) = max(rdicho,rccn)
    502          
    503 !!!!! Water ice mass is computed with radius at t+1 and their current number
    504 !!!!! Nccn is at t or t+1, depending on what has been done before
    505 !          zq(ig,l,igcm_h2o_ice) =
    506 !     &       (pi*rho_ice*zq(ig,l,igcm_ccn_number)*4./3.
    507 !     &         * rdicho*rdicho*rdicho -
    508 !     &         zq(ig,l,igcm_ccn_mass)*rho_ice/rho_dust)
    509 !     &       * tauscaling(ig)
    510 
    511 !!!!! Water ice mass is computed with radius at t+1 and their number at t+1
    512 !!!!! that is dirty, but this enforces the use of Nccn at t+1, whatever is done before
    513 !!!!! TO BE CLEANED ONE DAY
    514           var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN
    515           var2 = zq0(ig,l,igcm_ccn_mass)*tauscaling(ig)   + dM
    516           zq(ig,l,igcm_h2o_ice) =
    517      &       (pi*rho_ice*var1*4./3.
    518      &         * rdicho*rdicho*rdicho -
    519      &         var2*rho_ice/rho_dust)
    520      
    521      
    522       !!!! enforce realistic values (if the case of growth on Nccn(t) and condensation on Nccn(t+1))
    523           zq(ig,l,igcm_h2o_ice) =
    524      &     min(zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap),
    525      &     zq(ig,l,igcm_h2o_ice))
    526           zq(ig,l,igcm_h2o_ice) =
    527      &     max(1e-30,zq(ig,l,igcm_h2o_ice))
    528      
    529           zq(ig,l,igcm_h2o_vap) =
    530      &         zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap)
    531      &       - zq(ig,l,igcm_h2o_ice)
    532      
    533      
    534        if (test_flag) then
    535         print*, 'rice,ice,vap af',
    536      &     rice(ig,l),zq(ig,l,igcm_h2o_ice),zq(ig,l,igcm_h2o_vap)
    537         print*, 'satu bf, af',
    538      &     zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l),
    539      &     zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
    540        endif
    541      
    542      
    543 !!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST           
    544           if ((zq(ig,l,igcm_h2o_ice).le. -1e-8)
    545      &    .or. (zq(ig,l,igcm_h2o_vap).le. -1e-8)) then
    546            print*, 'NEGATIVE WATER'
    547            print*, 'ig,l', ig,l
    548            print*, 'satu', satu
    549            print*, 'vap, ice bf',
    550      &          zq0(ig,l,igcm_h2o_vap), zq0(ig,l,igcm_h2o_ice)
    551            print*, 'vap, ice af',
    552      &          zq(ig,l,igcm_h2o_vap), zq(ig,l,igcm_h2o_ice)
    553            print*, 'ccn N,M bf',
    554      &          zq0(ig,l,igcm_ccn_number), zq0(ig,l,igcm_ccn_mass)
    555            print*, 'ccn N,M af',
    556      &          zq(ig,l,igcm_ccn_number), zq(ig,l,igcm_ccn_mass)
    557            print*, 'tauscaling',
    558      &          tauscaling(ig)
    559            print*, 'req,rccn,rice bf,rice af',
    560      &         req,rccn,var1,rice(ig,l)
    561            print*, 'rmin, rmax', var2,var3
    562            print*, 'error_out,tdicho,ptimestep',
    563      &          error_out(ig,l),tdicho,ptimestep
    564            print*, ' '
    565           endif
    566 !!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST           
     375c       get resistance growth
     376        call growthrate(zt(ig,l),pplay(ig,l),
     377     &             real(ph2o/satu),rice(ig,l),res)
     378
     379
     380ccccccc  implicit scheme of mass growth
     381
     382        dMice =
     383     & (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
     384     &   /(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.)
     385
     386
     387! With the above scheme, dMice cannot be bigger than vapor,
     388! but can be bigger than all available ice.
     389       dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
     390       dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
     391
     392       zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
     393       zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
     394
     395
     396       countcells = countcells + 1
     397       
     398       ! latent heat release
     399       lw=(2834.3-0.28*(zt(ig,l)-To)-
     400     &     0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3
     401       pdtcloud(ig,l)= dMice*lw/cpp/ptimestep
    567402         
    568403         
    569         ENDIF !of if Nccn >1
    570          
    571404     
    572405!=============================================================
    … …  
    576409c         If all the ice particles sublimate, all the condensation
    577410c         nuclei are released:
    578           if (zq(ig,l,igcm_h2o_ice).le.1e-10) then
    579 !           for coherence
    580 !            dM = 0
    581 !            dN = 0
    582 !            dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig)
    583 !            dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
     411          if (zq(ig,l,igcm_h2o_ice).le.1.e-8) then
     412         
    584413c           Water
    585414            zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
    … …  
    594423            zq(ig,l,igcm_ccn_mass) = 0.
    595424            zq(ig,l,igcm_ccn_number) = 0.
     425
    596426          endif
    597          
    598 
    599 !        dN = dN/ tauscaling(ig)
    600 !        dM = dM/ tauscaling(ig)
    601 !c       Dust particles
    602 !        zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM
    603 !        zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN
    604 !c       CCNs
    605 !        zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
    606 !        zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
    607 
    608                
    609         pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass)
    610      &                         -zq0(ig,l,igcm_dust_mass))/ptimestep
    611         pdqcloud(ig,l,igcm_dust_number)=(zq(ig,l,igcm_dust_number)
    612      &                         -zq0(ig,l,igcm_dust_number))/ptimestep
    613         pdqcloud(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
    614      &                         -zq0(ig,l,igcm_ccn_mass))/ptimestep
    615         pdqcloud(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
    616      &                         -zq0(ig,l,igcm_ccn_number))/ptimestep
    617         pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap)
    618      &                         -zq0(ig,l,igcm_h2o_vap))/ptimestep
    619         pdqcloud(ig,l,igcm_h2o_ice)=(zq(ig,l,igcm_h2o_ice)
    620      &                         -zq0(ig,l,igcm_h2o_ice))/ptimestep
    621      
    622      
    623         count = count +1
    624        
    625        
    626         ELSE ! for coherence (rdust, rice computations etc ..)
    627           zq(ig,l,igcm_dust_mass)   = zq0(ig,l,igcm_dust_mass)
    628           zq(ig,l,igcm_dust_number) = zq0(ig,l,igcm_dust_number)
    629           zq(ig,l,igcm_ccn_mass)    = zq0(ig,l,igcm_ccn_mass)
    630           zq(ig,l,igcm_ccn_number)  = zq0(ig,l,igcm_ccn_number)
    631           zq(ig,l,igcm_h2o_ice)     = zq0(ig,l,igcm_h2o_ice)
    632           zq(ig,l,igcm_h2o_vap)     = zq0(ig,l,igcm_h2o_vap)
    633          
    634 ! pour les sorties de test
    635 !          satu_out(ig,l) = satu
    636 !          gr_out(ig,l) = 0
    637 !          dN_out(ig,l) = 0
    638 !          dM_out(ig,l) = 0
    639427         
    640         ENDIF ! end if (saturation ratio > 1) or (there is h2o_ice)
    641        
    642 !        ccnbf(ig,l) = zq0(ig,l,igcm_ccn_number) * tauscaling(ig)
    643 !        ccnaf(ig,l) = zq(ig,l,igcm_ccn_number) * tauscaling(ig)
    644 !
    645 !        satubf(ig,l) = zq0(ig,l,igcm_h2o_vap) / zqsat(ig,l)
    646 !        satuaf(ig,l) = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
    647        
    648        
    649 c-----update temperature (latent heat relase)     
    650           lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
    651           pdtcloud(ig,l)= -pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
    652          
    653 c----- update rice & rhocloud AFTER microphysic
    654           Mo   = zq(ig,l,igcm_h2o_ice) +
    655      &           zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
    656           No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
    657           rhocloud(ig,l) = zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
    658      &                     +zq(ig,l,igcm_ccn_mass)* tauscaling(ig)
    659      &                      / Mo * rho_dust
    660           rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
    661          
    662           if ((Mo.lt.1.e-20) .or. (No.le.1)) then
    663               rice(ig,l) = 1.e-8
    664           else
    665               rice(ig,l) =
    666      &          CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
    667           endif
    668          
    669           nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
    670 
    671 c----- update rdust and sedimentation radius
    672           rdust(ig,l)=
    673      &      CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
    674      &      max(zq(ig,l,igcm_dust_number),0.01))
    675           rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
    676          
    677           rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)*
    678      &                         (1.+nuice_sed)*(1.+nuice_sed),
    679      &                         rdust(ig,l) )
    680           rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
    681 
     428          ENDIF !of if Nccn>1
     429         
    682430        ENDDO ! of ig loop
    683431      ENDDO ! of nlayer loop
     432     
     433     
     434      ! Get cloud tendencies
     435      pdqcloud(1:ngrid,1:nlay,1:nq) =
     436     & (zq(1:ngrid,1:nlay,1:nq)-zq0(1:ngrid,1:nlay,1:nq))/ptimestep
     437     
    684438     
    685439     
    … …  
    689443      IF (test_flag) then
    690444     
    691        error2d(:) = 0.
    692        DO l=1,nlay
    693        DO ig=1,ngrid
    694          error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
    695        ENDDO
    696        ENDDO
    697 
    698       print*, 'count is ',count, ' i.e. ',
    699      &     count*100/(nlay*ngrid), '% for microphys computation'     
    700 
    701       IF (ngrid.ne.1) THEN ! 3D
     445!       error2d(:) = 0.
     446!       DO l=1,nlay
     447!       DO ig=1,ngrid
     448!         error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
     449!       ENDDO
     450!       ENDDO
     451
     452      print*, 'count is ',countcells, ' i.e. ',
     453     &     countcells*100/(nlay*ngrid), '% for microphys computation'
     454
     455!      IF (ngrid.ne.1) THEN ! 3D
    702456!         call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
    703457!     &                    satu_out)
    … …  
    706460!         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
    707461!     &                    dN_out)
    708          call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
    709      &                    error2d)
     462!         call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
     463!     &                    error2d)
    710464!         call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
    711465!     &                    zqsat)
    712       ENDIF
    713 
    714       IF (ngrid.eq.1) THEN ! 1D
    715          call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
    716      &                    error_out)
    717          call WRITEDIAGFI(ngrid,"satu_bf","satu before","kg/kg",1,
    718      &                    satubf)
    719          call WRITEDIAGFI(ngrid,"satu_af","satu after","kg/kg",1,
    720      &                    satuaf)
    721          call WRITEDIAGFI(ngrid,"vapbf","h2ovap before","kg/kg",1,
    722      &                    zq0(1,:,igcm_h2o_vap))
    723          call WRITEDIAGFI(ngrid,"vapaf","h2ovap after","kg/kg",1,
    724      &                    zq(1,:,igcm_h2o_vap))
    725          call WRITEDIAGFI(ngrid,"icebf","h2oice before","kg/kg",1,
    726      &                    zq0(1,:,igcm_h2o_ice))
    727          call WRITEDIAGFI(ngrid,"iceaf","h2oice after","kg/kg",1,
    728      &                    zq(1,:,igcm_h2o_ice))
    729          call WRITEDIAGFI(ngrid,"ccnbf","ccn before","/kg",1,
    730      &                    ccnbf)
    731          call WRITEDIAGFI(ngrid,"ccnaf","ccn after","/kg",1,
    732      &                    ccnaf)
     466!      ENDIF
     467
     468!      IF (ngrid.eq.1) THEN ! 1D
     469!         call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
     470!     &                    error_out)
     471!         call WRITEDIAGFI(ngrid,"satu_bf","satu before","kg/kg",1,
     472!     &                    satubf)
     473!         call WRITEDIAGFI(ngrid,"satu_af","satu after","kg/kg",1,
     474!     &                    satuaf)
     475!         call WRITEDIAGFI(ngrid,"vapbf","h2ovap before","kg/kg",1,
     476!     &                    zq0(1,:,igcm_h2o_vap))
     477!         call WRITEDIAGFI(ngrid,"vapaf","h2ovap after","kg/kg",1,
     478!     &                    zq(1,:,igcm_h2o_vap))
     479!         call WRITEDIAGFI(ngrid,"icebf","h2oice before","kg/kg",1,
     480!     &                    zq0(1,:,igcm_h2o_ice))
     481!         call WRITEDIAGFI(ngrid,"iceaf","h2oice after","kg/kg",1,
     482!     &                    zq(1,:,igcm_h2o_ice))
     483!         call WRITEDIAGFI(ngrid,"ccnbf","ccn before","/kg",1,
     484!     &                    ccnbf)
     485!         call WRITEDIAGFI(ngrid,"ccnaf","ccn after","/kg",1,
     486!     &                    ccnaf)
    733487c         call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
    734488c     &                    gr_out)
    … …  
    739493c         call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
    740494c     &                    dN_out)
    741          call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
    742      &                    zqsat)
    743          call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
    744      &                    satu_out)
    745          call WRITEDIAGFI(ngrid,"rice_sca","ice radius","m",1,
    746      &                    rice)
    747          call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
    748      &                    rdust)
    749          call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
    750      &                    rsedcloud)
    751          call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
    752      &                    rhocloud)
    753       ENDIF
     495!         call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
     496!     &                    zqsat)
     497!         call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
     498!     &                    satu_out)
     499!         call WRITEDIAGFI(ngrid,"rice_sca","ice radius","m",1,
     500!     &                    rice)
     501!         call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
     502!     &                    rdust)
     503!         call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
     504!     &                    rsedcloud)
     505!         call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
     506!     &                    rhocloud)
     507!      ENDIF
    754508     
    755509      ENDIF ! endif test_flag
    … …  
    772526cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    773527
    774       subroutine phi(rice,req,coeff1,coeff2,time)
    775      
    776       implicit none
    777      
    778       ! inputs
    779       real rice ! ice radius
    780       real req  ! ice radius at equilibirum
    781       real coeff1  ! coeff for the log
    782       real coeff2  ! coeff for the arctan
    783 
    784       ! output     
    785       real time
    786      
    787       !local
    788       real var
    789      
    790       ! 1.73205 is sqrt(3)
    791      
    792       var = max(
    793      &  abs(rice-req) / sqrt(rice*rice + rice*req  + req*req),1e-30)
    794            
    795        time =
    796      &   coeff1 *
    797      &   log( var )
    798      & + coeff2 * 1.73205 *
    799      &   atan( (2*rice+req) / (1.73205*req) )
    800      
    801       return
    802       end
    803      
    804      
    805      
     528c      subroutine phi(rice,req,coeff1,coeff2,time)
     529c     
     530c      implicit none
     531c     
     532c      ! inputs
     533c      real rice ! ice radius
     534c      real req  ! ice radius at equilibirum
     535c      real coeff1  ! coeff for the log
     536c      real coeff2  ! coeff for the arctan
     537c
     538c      ! output     
     539c      real time
     540c     
     541c      !local
     542c      real var
     543c     
     544c      ! 1.73205 is sqrt(3)
     545c     
     546c      var = max(
     547c     &  abs(rice-req) / sqrt(rice*rice + rice*req  + req*req),1e-30)
     548c           
     549c       time =
     550c     &   coeff1 *
     551c     &   log( var )
     552c     & + coeff2 * 1.73205 *
     553c     &   atan( (2*rice+req) / (1.73205*req) )
     554c     
     555c      return
     556c      end
     557     
     558     
     559     

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

    r627 r633  
    252252         nircorr=0      !default value. this is OK below 60 km.
    253253#else
    254          nircorr=1      !default value
     254         nircorr=0      !default value
    255255#endif
    256256         call getin("nircorr",nircorr)
    … …  
    531531         call getin("tituscap",tituscap)
    532532         write(*,*) "tituscap",tituscap
    533 
    534 !!!!!!!!!!!!!!!! TEMP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    535 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!     
    536          write(*,*) "temp tag for water caps ?"
    537          temptag=.false.
    538          call getin("temptag",temptag)
    539          write(*,*) " temptag = ",temptag
    540 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    541 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!             
     533                     
    542534
    543535         write(*,*) "photochemistry: include chemical species"

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

    r530 r633  
    2626
    2727c     Output
    28       DOUBLE PRECISION nucrate(nbin_cld)
     28   !   DOUBLE PRECISION nucrate(nbin_cld)
     29      REAL nucrate(nbin_cld)
    2930
    3031c     Local variables
    … …  
    119120            nucrate(i) = 0.
    120121          else
    121             nucrate(i)= sqrt ( fistar /
     122            nucrate(i)= real(sqrt ( fistar /
    122123     &               (3.*pi*kbz*temp*(gstar*gstar)) )
    123124     &                  * kbz * temp * rstar
    … …  
    126127     &                  * ( nh2o*rad_cld(i) )
    127128     &                  / ( zefshape * nus * m0 )
    128      &                  * dexp (deltaf)
     129     &                  * dexp (deltaf))
    129130          endif
    130131

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

    r627 r633  
    321321                               !   (direct comparison with TES)
    322322
    323 
    324323c Test 1d/3d scavenging
    325324      real h2otot(ngridmx)
    … …  
    355354      REAL zu2(ngridmx)
    356355c=======================================================================
     356
     357      print*,'physiq0', nq,nqmx
     358
    357359
    358360c 1. Initialisation:
    … …  
    851853     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
    852854     $ dtke_th,hfmax_th,wstar)
    853  
     855    
    854856         DO l=1,nlayer
    855857           DO ig=1,ngrid
    … …  
    985987
    986988           call watercloud(ngrid,nlayer,ptimestep,
    987      &                pplev,pplay,pdpsrf,zzlev,zzlay, pt,pdt,
     989     &                pplev,pplay,pdpsrf,zzlay, pt,pdt,
    988990     &                pq,pdq,zdqcloud,zdtcloud,
    989991     &                nq,tau,tauscaling,rdust,rice,nuice,
    990992     &                rsedcloud,rhocloud)
     993     
     994c Temperature variation due to latent heat release
    991995           if (activice) then
    992 c Temperature variation due to latent heat release
    993            DO l=1,nlayer
    994              DO ig=1,ngrid
    995                pdt(ig,l)=pdt(ig,l)+zdtcloud(ig,l)
    996              ENDDO
    997            ENDDO
     996               pdt(1:ngrid,1:nlayer) =
     997     &          pdt(1:ngrid,1:nlayer) +
     998     &          zdtcloud(1:ngrid,1:nlayer)
    998999           endif
    999 
     1000           
    10001001! increment water vapour and ice atmospheric tracers tendencies
    1001            IF (water) THEN
    1002           DO l=1,nlayer
    1003                DO ig=1,ngrid
    1004                  pdq(ig,l,igcm_h2o_vap)=
    1005      &                     pdq(ig,l,igcm_h2o_vap)+
    1006      &                    zdqcloud(ig,l,igcm_h2o_vap)
    1007                  pdq(ig,l,igcm_h2o_ice)=
    1008      &                    pdq(ig,l,igcm_h2o_ice)+
    1009      &                    zdqcloud(ig,l,igcm_h2o_ice)
    1010              ! There are negative values issues with some tracers (17/04)
    1011                  if (((pq(ig,l,igcm_h2o_ice) +
    1012      &                 pdq(ig,l,igcm_h2o_ice)*ptimestep)) .le. 0)
    1013      &             pdq(ig,l,igcm_h2o_ice) =               
    1014      &                 - pq(ig,l,igcm_h2o_ice)/ptimestep + 1e-30 
    1015                  if (((pq(ig,l,igcm_h2o_vap) +
    1016      &                 pdq(ig,l,igcm_h2o_vap)*ptimestep)) .le. 0)
    1017      &             pdq(ig,l,igcm_h2o_vap) =               
    1018      &                 - pq(ig,l,igcm_h2o_vap)/ptimestep + 1e-30
    1019      
    1020                  IF (scavenging) THEN
    1021                    pdq(ig,l,igcm_ccn_mass)=
    1022      &                      pdq(ig,l,igcm_ccn_mass)+
    1023      &                      zdqcloud(ig,l,igcm_ccn_mass)
    1024                    pdq(ig,l,igcm_ccn_number)=
    1025      &                      pdq(ig,l,igcm_ccn_number)+
    1026      &                      zdqcloud(ig,l,igcm_ccn_number)
    1027                    pdq(ig,l,igcm_dust_mass)=
    1028      &                      pdq(ig,l,igcm_dust_mass)+
    1029      &                      zdqcloud(ig,l,igcm_dust_mass)
    1030                    pdq(ig,l,igcm_dust_number)=
    1031      &                      pdq(ig,l,igcm_dust_number)+
    1032      &                      zdqcloud(ig,l,igcm_dust_number)
    1033 !!!!!!!!!!!!!!!!!!!!! We need to check that we have Nccn & Ndust > 0
    1034 !!!!!!!!!!!!!!!!!!!!! This is due to single precision rounding problems 
    1035                    if (((pq(ig,l,igcm_ccn_number) +
    1036      &                 pdq(ig,l,igcm_ccn_number)*ptimestep)) .le. 0)
    1037      &             then
    1038                        pdq(ig,l,igcm_ccn_number) =
    1039      &                 - pq(ig,l,igcm_ccn_number)/ptimestep + 1e-30
    1040                        pdq(ig,l,igcm_ccn_mass) =
    1041      &                 - pq(ig,l,igcm_ccn_mass)/ptimestep + 1e-30
    1042                    endif
    1043                    if (((pq(ig,l,igcm_dust_number) +
    1044      &                 pdq(ig,l,igcm_dust_number)*ptimestep)) .le. 0)
    1045      &             then
    1046                        pdq(ig,l,igcm_dust_number) =
    1047      &                 - pq(ig,l,igcm_dust_number)/ptimestep + 1e-30
    1048                        pdq(ig,l,igcm_dust_mass) =
    1049      &                 - pq(ig,l,igcm_dust_mass)/ptimestep + 1e-30
    1050                    endif
    1051 !!!!!!!!!!!!!!!!!!!!!
    1052 !!!!!!!!!!!!!!!!!!!!!
    1053                  ENDIF
    1054                ENDDO
    1055              ENDDO
    1056            ENDIF ! of IF (water) THEN
     1002           pdq(1:ngrid,1:nlayer,1:nq) =
     1003     &       pdq(1:ngrid,1:nlayer,1:nq) +
     1004     &       zdqcloud(1:ngrid,1:nlayer,1:nq)
     1005           
     1006! We need to check that we have Nccn & Ndust > 0
     1007! This is due to single precision rounding problems
     1008         if (scavenging) then
     1009
     1010           where (pq(:,:,igcm_ccn_mass) +
     1011     &      ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
     1012              pdq(:,:,igcm_ccn_mass) =
     1013     &        - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
     1014              pdq(:,:,igcm_ccn_number) =
     1015     &        - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
     1016           end where
     1017           where (pq(:,:,igcm_ccn_number) +
     1018     &      ptimestep*pdq(:,:,igcm_ccn_number) < 0.)
     1019              pdq(:,:,igcm_ccn_mass) =
     1020     &        - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
     1021              pdq(:,:,igcm_ccn_number) =
     1022     &        - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
     1023           end where
     1024           where (pq(:,:,igcm_dust_mass) +
     1025     &      ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
     1026              pdq(:,:,igcm_dust_mass) =
     1027     &        - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
     1028              pdq(:,:,igcm_dust_number) =
     1029     &        - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
     1030           end where
     1031           where (pq(:,:,igcm_dust_number) +
     1032     &      ptimestep*pdq(:,:,igcm_dust_number) < 0.)
     1033              pdq(:,:,igcm_dust_mass) =
     1034     &        - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
     1035              pdq(:,:,igcm_dust_number) =
     1036     &        - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
     1037           end where
     1038           
     1039         endif ! of if scavenging
     1040                     
    10571041
    10581042         END IF  ! of IF (water)
    … …  
    11011085     
    11021086     
    1103       !print*, 'h2o_ice zdqsed ds physiq', zdqsed(1,:,igcm_h2o_ice)
    1104 
    11051087           DO iq=1, nq
    11061088             DO l=1,nlayer
    … …  
    21212103     &                'Nbr.m-2.s-1',1,zdqsed(1,:,igcm_dust_number))
    21222104     
    2123         call WRITEDIAGFI(ngridmx,'zdqcloud*dt ice','cloud ice',
    2124      &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)*ptimestep)
    2125         call WRITEDIAGFI(ngridmx,'zdqcloud*dt vap','cloud vap',
    2126      &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap)*ptimestep)
    2127         call WRITEDIAGFI(ngridmx,'zdqdif*dt ice','dif ice',
    2128      &            'kg.m-2.s-1',1,zdqdif(1,:,igcm_h2o_ice)*ptimestep)
    2129         call WRITEDIAGFI(ngridmx,'zdqdif*dt vap','dif vap',
    2130      &            'kg.m-2.s-1',1,zdqdif(1,:,igcm_h2o_vap)*ptimestep)
    2131         call WRITEDIAGFI(ngridmx,'zdqdif*dt vap 0','dif vap',
    2132      &            'kg.m-2.s-1',0,zdqdif(1,1,igcm_h2o_vap)*ptimestep)
    2133 
    2134 !       if (scavenging) then
    2135 !         call WRITEDIAGFI(ngridmx,'zdqsed_ccnq','sedimentation q',
    2136 !     &                   'kg.m-2.s-1',1,zdqsed(1,:,igcm_ccn_mass))
    2137 !         call WRITEDIAGFI(ngridmx,'zdqsed_ccnN','sedimentation N',
    2138 !    &                 'Nbr.m-2.s-1',1,zdqsed(1,:,igcm_ccn_number))
    2139 !      endif
    2140 !       call WRITEDIAGFI(ngridmx,'zdqsed_ice','sedimentation q',
    2141 !     &                   'kg.m-2.s-1',1,zdqsed(1,:,igcm_h2o_ice))
    2142 !     
     2105        call WRITEDIAGFI(ngridmx,'zdqcloud_ice','cloud ice',
     2106     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
     2107        call WRITEDIAGFI(ngridmx,'zdqcloud_vap','cloud vap',
     2108     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_vap))
     2109        call WRITEDIAGFI(ngridmx,'zdqcloud','cloud ice',
     2110     &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice)
     2111     &                          +zdqcloud(1,:,igcm_h2o_vap))
     2112     
    21432113     
    21442114          call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,

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

    r520 r633  
    11      subroutine simpleclouds(ngrid,nlay,ptimestep,
    2      &             pplev,pplay,pzlev,pzlay,pt,pdt,
     2     &             pplev,pplay,pzlay,pt,pdt,
    33     &             pq,pdq,pdqcloud,pdtcloud,
    4      &             nq,tau,rice,nuice,rsedcloud)
     4     &             nq,tau,rice)
    55      implicit none
    66c------------------------------------------------------------------
    … …  
    4343      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
    4444      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
    45       REAL pzlev(ngrid,nlay+1)   ! altitude at layer boundaries
    4645      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
    4746      REAL pt(ngrid,nlay)        ! temperature at the middle of the
    … …  
    5756      REAL rice(ngrid,nlay)      ! Ice mass mean radius (m)
    5857                                 ! (r_c in montmessin_2004)
    59       REAL nuice(ngrid,nlay)     ! Estimated effective variance
    60                                  !   of the size distribution
    61       real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
    6258      real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
    6359                                   !   H2O(kg/kg.s-1)
    … …  
    181177     &      +ccntyp(ig,l)*(4./3.)*pi*rdusttyp(ig,l)**3.)
    182178     &      /(ccntyp(ig,l)*4./3.*pi) ), rdusttyp(ig,l))
    183 c         Effective variance of the size distribution
    184           nuice(ig,l)=nuice_ref
    185 
    186 c         Sedimentation radius:
    187 c         ~~~~~~~~~~~~~~~~~~~~
    188 
    189           rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)**3.,
    190      &                         rdusttyp(ig,l) )
    191           rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
     179
    192180
    193181        enddo ! of do ig=1,ngrid

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

    r528 r633  
    2323      REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
    2424      REAL icedryness ! ice dryness
     25     
     26      ! longwatercaptag is watercaptag. Trick for some compilers
     27      LOGICAL, DIMENSION(100000) :: longwatercaptag
    2528
    2629      EXTERNAL ISMIN,ISMAX
    … …  
    3033! i.e based only on lat & lon values of each physical point.
    3134! Choose true to get a carefully choosen distribution for GCM resolutions 32x24 or 64x48
    32 ! For vizualisation : soon ...
     35! For vizualisation : /u/tnalmd/bin/watercaps gcm_txt_output
    3336       LOGICAL,SAVE :: improvedicelocation = .true.
    3437c
    … …  
    5558
    5659         alternate = 0
    57          temptag = .true.
     60         
     61         if (ngridmx .ne. 1) then
     62           watercaptag(:) = .false.
     63           longwatercaptag(:) = .false.
     64         endif
    5865         
    5966         write(*,*) "Ice dryness ?"
    … …  
    6168         call getin("icedryness",icedryness)
    6269         write(*,*) " icedryness = ",icedryness
    63 
     70         
    6471       
    6572#ifdef MESOSCALE
     73
    6674      do ig=1,ngridmx
     75
    6776         !write(*,*) "all qsurf to zero. dirty."
    6877         do iq=1,nqmx
    … …  
    8190                 watercaptag(ig)  = .false.
    8291                 dryness(ig)      = 1.
    83          endif 
    84       enddo
     92         endif
     93         
     94      endif
    8595#else
    86        ! print*,'ngridmx',ngridmx
    8796
    8897      if (improvedicelocation) then
    8998     
    9099        if (ngridmx .eq. 738) then ! hopefully 32x24
     100           
    91101          print*,'water ice caps distribution for 32x24 resolution:'
    92           watercaptag(1:9)    = .true. ! central cap - core
    93           watercaptag(26:33)  = .true. ! central cap
     102          longwatercaptag(1:9)    = .true. ! central cap - core
     103          longwatercaptag(26:33)  = .true. ! central cap
    94104           
    95105        else if (ngridmx .eq. 3010) then ! hopefully 64x48
     106
     107! For latitudes:
     108! [ 67.5   71.25  75.    78.75  82.5   86.25]
     109! The water ice caps should be (according to TES temperatures):
     110! [  8.63e-03   1.02e+00   5.99e+00   2.66e+01   5.65e+01]
     111
    96112          print*,'water ice caps distribution for 64x48 resolution:'
    97           watercaptag(1:65)   = .true. ! central cap - core
    98 !          watercaptag(85)     = .true. ! central cap
    99           watercaptag(83:85)  = .true. ! central cap - BIGGER
    100           watercaptag(93:114) = .true. ! central cap
    101 !          watercaptag(254)    = .true. ! Korolev crater
    102 !          watercaptag(255)    = .true. ! Korolev crater --DECALE
    103           watercaptag(242:257)= .true. ! Korolev crater -- VERY BIG
    104 !--------------------------------------------------------------------         
    105 !          watercaptag(136)    = .true. ! outlier, lat = 78.75
    106 !          watercaptag(138)    = .true. ! outlier, lat = 78.75
    107 !          watercaptag(140)    = .true. ! outlier, lat = 78.75
    108 !          watercaptag(142)    = .true. ! outlier, lat = 78.75
    109 !          watercaptag(183)    = .true. ! outlier, lat = 78.75
    110 !          watercaptag(185)    = .true. ! outlier, lat = 78.75
    111 !          watercaptag(187)    = .true. ! outlier, lat = 78.75
    112 !          watercaptag(189)    = .true. ! outlier, lat = 78.75
    113 !          watercaptag(191)    = .true. ! outlier, lat = 78.75
    114 !          watercaptag(193)    = .true. ! outlier, lat = 78.75
    115 !--------------------------------------------------------------------         
    116           watercaptag(135:142)= .true. ! BIG outlier, lat = 78.75
    117           watercaptag(181:193)= .true. ! BIG outlier, lat = 78.75
     113          longwatercaptag(1:65)   = .true. ! central cap - core
     114          longwatercaptag(75:85)  = .true. ! central cap
     115          longwatercaptag(93:114) = .true. ! central cap
     116!---------------------   OUTLIERS  ----------------------------         
     117          longwatercaptag(136)    = .true. ! outlier, lat = 78.75
     118          longwatercaptag(138)    = .true. ! outlier, lat = 78.75
     119          longwatercaptag(140)    = .true. ! outlier, lat = 78.75
     120          longwatercaptag(142)    = .true. ! outlier, lat = 78.75
     121          longwatercaptag(161)    = .true. ! outlier, lat = 78.75
     122          longwatercaptag(163)    = .true. ! outlier, lat = 78.75
     123          longwatercaptag(165)    = .true. ! outlier, lat = 78.75
     124          longwatercaptag(183)    = .true. ! outlier, lat = 78.75
     125          longwatercaptag(185)    = .true. ! outlier, lat = 78.75
     126          longwatercaptag(187)    = .true. ! outlier, lat = 78.75
     127          longwatercaptag(189)    = .true. ! outlier, lat = 78.75
     128          longwatercaptag(191)    = .true. ! outlier, lat = 78.75
     129          longwatercaptag(193)    = .true. ! outlier, lat = 78.75
     130          longwatercaptag(194)    = .true. ! outlier, lat = 75
     131          longwatercaptag(203)    = .true. ! outlier, lat = 75
     132          longwatercaptag(207)    = .true. ! outlier, lat = 75
     133          longwatercaptag(242)    = .true. ! outlier, lat = 75
     134          longwatercaptag(244)    = .true. ! outlier, lat = 75
     135          longwatercaptag(246)    = .true. ! outlier, lat = 75
     136          longwatercaptag(248)    = .true. ! outlier, lat = 75
     137          longwatercaptag(250)    = .true. ! outlier, lat = 75
     138          longwatercaptag(252)    = .true. ! outlier, lat = 75
     139          longwatercaptag(254)    = .true. ! outlier, lat = 75
     140          longwatercaptag(256:257)= .true. ! outlier, lat = 75
     141!!-------------------------   OLD   ----------------------------         
     142!!          longwatercaptag(83:85)  = .true.
     143!!          longwatercaptag(135:142)  = .true.
     144!!          longwatercaptag(181:193)  = .true.
     145!!          longwatercaptag(242:257)  = .true.
     146
    118147           
    119148        else if (ngridmx .ne. 1) then
    … …  
    129158        do ig=1,ngridmx
    130159          dryness (ig) = icedryness
    131           !!! OU alors tu mets dryness = icedrag sur outliers et 1 ailleurs
    132           !!! OU remettre comme c'était avant Aymeric et pis c'est tout
    133           if (watercaptag(ig)) then
     160
     161          if (longwatercaptag(ig)) then
     162             watercaptag(ig) = .True.
    134163             print*,'ice water cap', lati(ig)*180./pi,
    135164     .              long(ig)*180./pi, ig
    … …  
    137166        enddo
    138167
    139 !!------------------------ test -----------------------------
    140 !!-----------------------------------------------------------   
    141 !!      else if (1) then ! DUMMY !!, CAREFUL !!, TOUT CA ....
    142 !!       
    143 !!         do ig=1,ngridmx
    144 !!          dryness (ig) = 1
    145 !!          if (albedodat(ig) .ge. 0.35) then
    146 !!             watercaptag(ig)    = .true.
    147 !!             print*,' albedo criteria ice water cap', lati(ig)*180./pi,
    148 !!     .              long(ig)*180./pi, ig
    149 !!          endif
    150 !!        enddo
    151 !!-----------------------------------------------------------   
    152 !!-----------------------------------------------------------   
    153      
     168
    154169      else
    155170
    … …  
    213228     .         then
    214229     
    215              if (temptag) then
    216230             
    217231               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
    218232                if (ngridmx.ne.1) watercaptag(ig)=.true.
    219 !                  print*,'ice water cap', lati(ig)*180./pi,
    220 !     .              long(ig)*180./pi, ig
    221                   !dryness(ig) = 1.
     233                  print*,'ice water cap', lati(ig)*180./pi,
     234     .              long(ig)*180./pi, ig
     235                  dryness(ig) = 1.
    222236                  alternate = 1
    223237                else
    224238                  alternate = 0
    225239                endif !end if alternate = 0
    226                
    227              else
    228              
    229 !               if (ngridmx.ne.1) watercaptag(ig)=.true.
    230 !                  write(*,*) "outliers ", lati(ig)*180./pi,
    231 !     .              long(ig)*180./pi
    232      
    233              endif ! end if temptag
     240
    234241             
    235242           endif
    … …  
    259266     .            long(ig)*180./pi, ig
    260267              endif
    261               !dryness(ig) = 1.
     268             dryness(ig) = 1.
    262269           endif
    263270
    … …  
    266273c--------------------------------------------------------
    267274
    268            if (lati(ig)*180./pi.gt.84) then
     275           if (abs(lati(ig)*180./pi).gt.80) then
    269276             print*,'ice water cap', lati(ig)*180./pi,
    270277     .         long(ig)*180./pi, ig

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

    r631 r633  
    122122     &            ccn0*(4./3.)*pi*rdust(ig,l)**3.) /
    123123     &            (ccn0*4./3.*pi)),rdust(ig,l) )
     124                rice(ig,l)=min(max(rice(ig,l),1.e-10),500.e-6)
    124125                nuice(ig,l) = nuice_ref
    125126              ENDDO
    … …  
    142143                rice(ig,l) =
    143144     &           CBRT( Mo/No * 0.75 / pi / rhocloud(ig,l))
     145                rice(ig,l)=min(max(rice(ig,l),1.e-10),500.e-6)
    144146                nuice(ig,l) = nuice_ref
    145147              ENDDO

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

    r626 r633  
    1        SUBROUTINE watercloud(ngrid,nlay, ptimestep,
    2      &                pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
     1       SUBROUTINE watercloud(ngrid,nlay,ptimestep,
     2     &                pplev,pplay,pdpsrf,pzlay,pt,pdt,
    33     &                pq,pdq,pdqcloud,pdtcloud,
    44     &                nq,tau,tauscaling,rdust,rice,nuice,
    55     &                rsedcloud,rhocloud)
     6! to use  'getin'
     7      USE ioipsl_getincom
    68      IMPLICIT NONE
     9
    710
    811c=======================================================================
    … …  
    1316c    - An improved microphysical scheme (see improvedclouds.F)
    1417c
     18c  There is a time loop specific to cloud formation
     19c  due to timescales smaller than the GCM integration timestep.
     20c
    1521c  Authors: Franck Montmessin, Francois Forget, Ehouarn Millour,
    1622c           J.-B. Madeleine, Thomas Navarro
    1723c
    18 c  2004 - Oct. 2011
     24c  2004 - 2012
    1925c=======================================================================
    2026
    … …  
    3541
    3642      INTEGER ngrid,nlay
    37       integer nq                 ! nombre de traceurs
     43      INTEGER nq                 ! nombre de traceurs
    3844      REAL ptimestep             ! pas de temps physique (s)
    3945      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
    4046      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
    41       REAL pdpsrf(ngrid)         ! tendance surf pressure
    42       REAL pzlev(ngrid,nlay+1)   ! altitude at layer boundaries
     47      REAL pdpsrf(ngrid)         ! tendence surf pressure
    4348      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
    4449      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
    45       REAL pdt(ngrid,nlay)       ! tendance temperature des autres param.
     50      REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
    4651
    4752      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
    48       real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
    49 
    50       REAL tau(ngridmx,naerkind)   ! Column dust optical depth at each point
    51       REAL tauscaling(ngridmx)     ! Convertion factor for dust amount
    52       real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
     53      real pdq(ngrid,nlay,nq)    ! tendence avant condensation  (kg/kg.s-1)
     54
     55      REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
     56      REAL tauscaling(ngridmx)   ! Convertion factor for dust amount
     57      real rdust(ngridmx,nlay)  ! Dust geometric mean radius (m)
    5358
    5459c   Outputs:
    5560c   -------
    5661
    57       real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation H2O(kg/kg.s-1)
    58       REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
    59                                    !   a la chaleur latente
     62      real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
     63      REAL pdtcloud(ngrid,nlay)    ! tendence temperature due
     64                                   ! a la chaleur latente
    6065
    6166      REAL rice(ngrid,nlay)    ! Ice mass mean radius (m)
    … …  
    6368      REAL nuice(ngrid,nlay)   ! Estimated effective variance
    6469                               !   of the size distribution
    65       real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
    66       real rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)
     70      real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius
     71      real rhocloud(ngridmx,nlay)  ! Cloud density (kg.m-3)
    6772
    6873c   local:
    6974c   ------
    70 
    71       INTEGER ig,l
     75     
     76      ! for ice radius computation
     77      REAL Mo,No
     78      REAl ccntyp
     79     
     80      ! for time loop
     81      INTEGER microstep  ! time subsampling step variable
     82      INTEGER imicro     ! time subsampling for coupled water microphysics & sedimentation
     83      SAVE imicro
     84      REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
     85      SAVE microtimestep
     86     
     87      ! tendency given by clouds (inside the micro loop)
     88      REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud
     89      REAL subpdtcloud(ngrid,nlay)    ! cf. pdtcloud
     90
     91      ! global tendency (clouds+physics)
     92      REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
     93      REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
     94     
     95      REAL CBRT
     96      EXTERNAL CBRT
     97
     98
     99      INTEGER iq,ig,l
    72100      LOGICAL,SAVE :: firstcall=.true.
    73101
    … …  
    92120        write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
    93121        write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice
     122               
     123        write(*,*) "time subsampling for microphysic ?"
     124#ifdef MESOSCALE
     125        imicro = 2
     126#else
     127        imicro = 15
     128#endif
     129        call getin("imicro",imicro)
     130        write(*,*)"imicro = ",imicro
     131       
     132        microtimestep = ptimestep/real(imicro)
     133        write(*,*)"Physical timestep is",ptimestep
     134        write(*,*)"Microphysics timestep is",microtimestep
    94135
    95136        firstcall=.false.
    96137      ENDIF ! of IF (firstcall)
    97138     
    98 
    99 c     Main call to the different cloud schemes:
    100       IF (microphys) THEN
    101         CALL improvedclouds(ngrid,nlay,ptimestep,
    102      &             pplev,pplay,pt,pdt,
    103      &             pq,pdq,pdqcloud,pdtcloud,
    104      &             nq,tauscaling,rdust,rice,nuice,
    105      &             rsedcloud,rhocloud)
    106       ELSE
    107         CALL simpleclouds(ngrid,nlay,ptimestep,
    108      &             pplev,pplay,pzlev,pzlay,pt,pdt,
    109      &             pq,pdq,pdqcloud,pdtcloud,
    110      &             nq,tau,rice,nuice,rsedcloud)
    111       ENDIF
     139c-----Initialization
     140      subpdq(1:ngrid,1:nlay,1:nq) = 0
     141      subpdt(1:ngrid,1:nlay)      = 0
     142     
     143      ! default value if no ice
     144      rhocloud(1:ngrid,1:nlay) = rho_dust
     145
     146
     147
     148c------------------------------------------------------------------
     149c Time subsampling for microphysics
     150c------------------------------------------------------------------
     151      DO microstep=1,imicro
     152     
     153c-------------------------------------------------------------------
     154c   1.  Tendencies:
     155c------------------
     156
     157
     158c------ Temperature tendency subpdt
     159        ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
     160        ! If imicro=1 subpdt is the same as pdt
     161        DO l=1,nlay
     162          DO ig=1,ngrid
     163             subpdt(ig,l) = subpdt(ig,l)
     164     &        + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
     165          ENDDO
     166        ENDDO
     167c------ Tracers tendencies subpdq
     168c------ At each micro timestep we add pdq in order to have a stepped entry
     169        IF (microphys) THEN
     170          DO l=1,nlay
     171            DO ig=1,ngrid
     172              subpdq(ig,l,igcm_dust_mass) =
     173     &            subpdq(ig,l,igcm_dust_mass)
     174     &          + pdq(ig,l,igcm_dust_mass)
     175              subpdq(ig,l,igcm_dust_number) =
     176     &            subpdq(ig,l,igcm_dust_number)
     177     &          + pdq(ig,l,igcm_dust_number)
     178              subpdq(ig,l,igcm_ccn_mass) =
     179     &            subpdq(ig,l,igcm_ccn_mass)
     180     &          + pdq(ig,l,igcm_ccn_mass)
     181              subpdq(ig,l,igcm_ccn_number) =
     182     &            subpdq(ig,l,igcm_ccn_number)
     183     &          + pdq(ig,l,igcm_ccn_number)
     184            ENDDO
     185          ENDDO
     186        ENDIF
     187        DO l=1,nlay
     188          DO ig=1,ngrid
     189            subpdq(ig,l,igcm_h2o_ice) =
     190     &          subpdq(ig,l,igcm_h2o_ice)
     191     &        + pdq(ig,l,igcm_h2o_ice)
     192            subpdq(ig,l,igcm_h2o_vap) =
     193     &          subpdq(ig,l,igcm_h2o_vap)
     194     &        + pdq(ig,l,igcm_h2o_vap)
     195          ENDDO
     196        ENDDO
     197       
     198       
     199c-------------------------------------------------------------------
     200c   2.  Main call to the different cloud schemes:
     201c------------------------------------------------
     202        IF (microphys) THEN
     203           CALL improvedclouds(ngrid,nlay,microtimestep,
     204     &             pplay,pt,subpdt,
     205     &             pq,subpdq,subpdqcloud,subpdtcloud,
     206     &             nq,tauscaling)
     207
     208        ELSE
     209           CALL simpleclouds(ngrid,nlay,microtimestep,
     210     &             pplay,pzlay,pt,subpdt,
     211     &             pq,subpdq,subpdqcloud,subpdtcloud,
     212     &             nq,tau)
     213        ENDIF
     214
     215
     216c-------------------------------------------------------------------
     217c   3.  Updating tendencies after cloud scheme:
     218c-----------------------------------------------
     219
     220        IF (microphys) THEN
     221          DO l=1,nlay
     222            DO ig=1,ngrid
     223              subpdq(ig,l,igcm_dust_mass) =
     224     &            subpdq(ig,l,igcm_dust_mass)
     225     &          + subpdqcloud(ig,l,igcm_dust_mass)
     226              subpdq(ig,l,igcm_dust_number) =
     227     &            subpdq(ig,l,igcm_dust_number)
     228     &          + subpdqcloud(ig,l,igcm_dust_number)
     229              subpdq(ig,l,igcm_ccn_mass) =
     230     &            subpdq(ig,l,igcm_ccn_mass)
     231     &          + subpdqcloud(ig,l,igcm_ccn_mass)
     232              subpdq(ig,l,igcm_ccn_number) =
     233     &            subpdq(ig,l,igcm_ccn_number)
     234     &          + subpdqcloud(ig,l,igcm_ccn_number)
     235            ENDDO
     236          ENDDO
     237        ENDIF
     238        DO l=1,nlay
     239          DO ig=1,ngrid
     240            subpdq(ig,l,igcm_h2o_ice) =
     241     &          subpdq(ig,l,igcm_h2o_ice)
     242     &        + subpdqcloud(ig,l,igcm_h2o_ice)
     243            subpdq(ig,l,igcm_h2o_vap) =
     244     &          subpdq(ig,l,igcm_h2o_vap)
     245     &        + subpdqcloud(ig,l,igcm_h2o_vap)
     246          ENDDO
     247        ENDDO
     248     
     249 
     250      ENDDO ! of DO microstep=1,imicro
     251     
     252c-------------------------------------------------------------------
     253c   6.  Compute final tendencies after time loop:
     254c------------------------------------------------
     255
     256c------ Temperature tendency pdtcloud
     257       DO l=1,nlay
     258         DO ig=1,ngrid
     259             pdtcloud(ig,l) =
     260     &         subpdt(ig,l)/real(imicro)-pdt(ig,l)
     261          ENDDO
     262       ENDDO
     263c------ Tracers tendencies pdqcloud
     264       DO iq=1,nq
     265        DO l=1,nlay
     266         DO ig=1,ngrid
     267            pdqcloud(ig,l,iq) = subpdq(ig,l,iq)/real(imicro)
     268     &       - pdq(ig,l,iq)
     269         ENDDO
     270        ENDDO
     271       ENDDO
     272
     273c------- Due to stepped entry, other processes tendencies can add up to negative values
     274c------- Therefore, enforce positive values and conserve mass
     275
     276       IF(microphys) THEN
     277        DO l=1,nlay
     278         DO ig=1,ngrid
     279          IF (pq(ig,l,igcm_ccn_number) +
     280     &      ptimestep* (pdq(ig,l,igcm_ccn_number) +
     281     &        pdqcloud(ig,l,igcm_ccn_number)) .le. 0.) THEN
     282         pdqcloud(ig,l,igcm_ccn_number) =
     283     &     - pq(ig,l,igcm_ccn_number)/ptimestep
     284     &     - pdq(ig,l,igcm_ccn_number)
     285         pdqcloud(ig,l,igcm_dust_number) = 
     286     &     -pdqcloud(ig,l,igcm_ccn_number)
     287         pdqcloud(ig,l,igcm_ccn_mass) =
     288     &     - pq(ig,l,igcm_ccn_mass)/ptimestep
     289     &     - pdq(ig,l,igcm_ccn_mass)
     290         pdqcloud(ig,l,igcm_dust_mass) =
     291     &     -pdqcloud(ig,l,igcm_ccn_mass)
     292          ENDIF
     293         ENDDO
     294        ENDDO
     295       ENDIF
     296
     297        DO l=1,nlay
     298         DO ig=1,ngrid
     299          IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
     300     &       (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice))
     301     &       .le. 1.e-8) THEN
     302           pdqcloud(ig,l,igcm_h2o_ice) =
     303     &     - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
     304           pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
     305          ENDIF
     306          IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
     307     &       (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
     308     &       .le. 1.e-8) THEN
     309           pdqcloud(ig,l,igcm_h2o_vap) =
     310     &     - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
     311           pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
     312          ENDIF
     313         ENDDO
     314        ENDDO
     315
     316
     317c------Update the ice and dust particle size "rice" for output or photochemistry
     318c------Only rsedcloud is used for the water cycle
     319
     320      IF(scavenging) THEN
     321        DO l=1, nlay
     322         DO ig=1,ngrid
     323
     324          rdust(ig,l)=
     325     &      CBRT(r3n_q*(pq(ig,l,igcm_dust_mass)+
     326     &           (pdq(ig,l,igcm_dust_mass)
     327     &           +pdqcloud(ig,l,igcm_dust_mass))*ptimestep)/
     328     &      (pq(ig,l,igcm_dust_number)+
     329     &           (pdq(ig,l,igcm_dust_number)+
     330     &            pdqcloud(ig,l,igcm_dust_number)*ptimestep)))
     331          rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
     332
     333          Mo = pq(ig,l,igcm_h2o_ice)
     334     &         + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep
     335     &         + (pq(ig,l,igcm_ccn_mass)
     336     &          + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
     337     &         *tauscaling(ig) + 1.e-30
     338          No = (pq(ig,l,igcm_ccn_number)
     339     &          + pdqcloud(ig,l,igcm_ccn_number)*ptimestep)
     340     &         *tauscaling(ig) + 1.e-30
     341          rhocloud(ig,l) = (pq(ig,l,igcm_h2o_ice) + 
     342     &                   pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)
     343     &                   / Mo * rho_ice
     344     &                  + (pq(ig,l,igcm_ccn_mass)
     345     &                  + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
     346     &                   *tauscaling(ig)/ Mo * rho_dust
     347          rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
     348          if ((Mo.lt.1.e-15) .or. (No.le.1.)) then
     349              rice(ig,l) = 1.e-8
     350          else
     351              !! AS: only perform computations if conditions not met
     352              rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.)
     353              rice(ig,l)=min(max(rice(ig,l),1.e-10),500.e-6)
     354          endif
     355          rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)*
     356     &                         (1.+nuice_sed)*(1.+nuice_sed),
     357     &                         rdust(ig,l) )
     358          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
     359         ENDDO
     360        ENDDO
     361       ELSE
     362        DO l=1,nlay
     363          DO ig=1,ngrid
     364
     365          rdust(ig,l)=
     366     &      CBRT(r3n_q*(pq(ig,l,igcm_dust_mass)+
     367     &           (pdq(ig,l,igcm_dust_mass)
     368     &           +pdqcloud(ig,l,igcm_dust_mass))*ptimestep)/
     369     &      (pq(ig,l,igcm_dust_number)+
     370     &           (pdq(ig,l,igcm_dust_number)+
     371     &            pdqcloud(ig,l,igcm_dust_number)*ptimestep)))
     372          rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
     373
     374            ccntyp =
     375     &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-pzlay(ig,l)/10000.)
     376            ccntyp = ccntyp /ccn_factor
     377            rice(ig,l)=max( CBRT ( ((pq(ig,l,igcm_h2o_ice)
     378     &       + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)/rho_ice
     379     &      +ccntyp*(4./3.)*pi*rdust(ig,l)*rdust(ig,l)*rdust(ig,l))
     380     &      /(ccntyp*4./3.*pi) ), rdust(ig,l))
     381            rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)*
     382     &                         (1.+nuice_sed)*(1.+nuice_sed),
     383     &                         rdust(ig,l) )
     384            rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
     385          ENDDO
     386        ENDDO
     387       ENDIF ! of IF(scavenging)
     388       
     389       
     390! used for rad. transfer calculations
     391! nuice is constant because a lognormal distribution is prescribed
     392      nuice(1:ngrid,1:nlay)=nuice_ref
     393
     394
     395!--------------------------------------------------------------
     396!--------------------------------------------------------------
    112397     
    113398
    … …  
    123408      end do
    124409
    125       RETURN
     410c=======================================================================
     411
    126412      END
    127413 
Note: See TracChangeset for help on using the changeset viewer.