Changeset 2966 for trunk


Ignore:
Timestamp:
May 19, 2023, 12:49:25 PM (2 years ago)
Author:
jnaar
Message:

Architecture changes in watercloud_mod, improvedclouds_mod in preparation for adaptative timestep of clouds (not yet implemented). Changes prevent bit by bit correspondance with previous revision. "simpleclouds" routine broken.

Location:
trunk/LMDZ.MARS
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/changelog.txt

    r2965 r2966  
    40434043+ fixed unit attribute of surface/soil water densities in field_def_physics_mars.xml
    40444044
     4045== 19/05/2023 == JN
     4046+ Architecture change in watercloud_mod.F, improvedclouds_mod.F :
     4047Instead of computing all subtimesteps simultaneously, we now loop on
     4048(ngrid,nlay) first. This is to allow for a future adaptative timestep.
     4049+ Second architecture change in the computations of tendencies within
     4050watercloud_mod.F : we now increment directly tracers and temperature instead
     4051of tendencies. This may allow future developments for optimizations.
     4052+ Because of second change, no bit by bit correspondance with previous
     4053revision (truncature stuff)
     4054+ "simpleclouds" routine should be broken in this build (microphys=.false.)
     4055+ Introducing a flag "cloud_adapt_ts" (false by default) to set the adaptative timestep (not yet ready). Otherwise fully retrocompatible.

  • trunk/LMDZ.MARS/libf/phymars/callkeys.h

    r2823 r2966  
    1515     &   ,calltherm,callrichsl,callslope,tituscap,callyamada4,co2clouds &
    1616     &   ,co2useh2o,meteo_flux,activeco2ice,CLFvaryingCO2,spantCO2      &
    17      &   ,CLFvarying,satindexco2,rdstorm,topflows,calllott_nonoro        &
     17     &   ,CLFvarying,satindexco2,rdstorm,topflows,calllott_nonoro       &
    1818     &   ,latentheat_surfwater,gwd_convective_source,startphy_file      &
    19      &   ,hdo,hdofrac,cst_cap_albedo,temp_dependant_m,refill_watercap
     19     &   ,hdo,hdofrac,cst_cap_albedo,temp_dependant_m,refill_watercap   &
     20     &   ,cloud_adapt_ts
    2021!$OMP THREADPRIVATE(/callkeys_l/)
    2122
    … …  
    7677      logical temp_dependant_m ! temperature-dependant water contact parameter
    7778      logical refill_watercap ! h2o_ice_s is converted to watercap when above threshold
     79      logical cloud_adapt_ts ! adaptative timestep for cloud microphysics
    7880      logical sedimentation
    7981      logical activice,tifeedback,supersat,caps

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

    r2916 r2966  
    741741           write(*,*) "mteta = ", mteta
    742742        endif
     743
     744! Adaptative timestep for cloud microphysics (JN 2023)
     745         write(*,*)"Adaptative timestep for cloud",
     746     &              " microphysics ? (default is false)"
     747         cloud_adapt_ts=.false. ! default value
     748         call getin_p("cloud_adapt_ts",cloud_adapt_ts)
     749         write(*,*)"cloud_adapt_ts= ",cloud_adapt_ts
    743750       
    744751! scavenging

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

    r2932 r2966  
    55      CONTAINS
    66
    7       subroutine improvedclouds(ngrid,nlay,microtimestep,
    8      &             pplay,pteff,sum_subpdt,
    9      &             pqeff,sum_subpdq,subpdqcloud,subpdtcloud,
    10      &             nq,tauscaling)
     7      subroutine improvedclouds(microtimestep,
     8     &             pplay,pt,pq,
     9     &             subpdqcloud,subpdtcloud,
     10     &             nq,tauscaling,mmean)
    1111      USE updaterad, ONLY: updaterice_micro, updaterccn
    1212      USE watersat_mod, ONLY: watersat
    … …  
    1717     &                      igcm_hdo_vap,igcm_hdo_ice,
    1818     &                      qperemin
    19       use conc_mod, only: mmean
    2019      use comcstfi_h, only: pi, cpp
    2120      use write_output_mod, only: write_output
    … …  
    4241c           T. Navarro, debug,correction, new scheme (October-April 2011)
    4342c           A. Spiga, optimization (February 2012)
     43c           J. Naar, from global to local (no more ngrid, nlay) to allow
     44c           different microphysical timesteps (May 2023)
    4445c------------------------------------------------------------------
    4546#include "callkeys.h"
    … …  
    4849c     Inputs/outputs:
    4950
    50       INTEGER, INTENT(IN) :: ngrid,nlay
    5151      INTEGER, INTENT(IN) :: nq                 ! nombre de traceurs
    5252      REAL, INTENT(IN) :: microtimestep         ! pas de temps physique (s)
    53       REAL, INTENT(IN) :: pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)           
    54       REAL, INTENT(IN) :: pteff(ngrid,nlay)    ! temperature at the middle of the
     53      REAL, INTENT(IN) :: pplay     ! pression au milieu des couches (Pa)           
     54      REAL, INTENT(IN) :: pt ! temperature at the middle of the
    5555                                                ! layers (K)
    56       REAL, INTENT(IN) :: sum_subpdt(ngrid,nlay)! tendance temperature des autres
    5756                                                !   param.
    58       REAL, INTENT(IN) :: pqeff(ngrid,nlay,nq)  ! traceur (kg/kg)
    59       REAL, INTENT(IN) :: sum_subpdq(ngrid,nlay,nq)    ! tendance avant condensation
     57      REAL, INTENT(IN) :: pq(nq)  ! traceur (kg/kg)
    6058                                                       !   (kg/kg.s-1)
    61       REAL, INTENT(IN) :: tauscaling(ngrid)     ! Convertion factor for qdust and Ndust
    62      
    63       REAL, INTENT(OUT) :: subpdqcloud(ngrid,nlay,nq)  ! tendance de la condensation
     59      REAL, INTENT(IN) :: tauscaling     ! Convertion factor for qdust and Ndust
     60      REAL, INTENT(IN) :: mmean ! Mean atmospheric mass
     61     
     62      REAL, INTENT(OUT) :: subpdqcloud(nq)  ! tendance de la condensation
    6463                                                       !   H2O(kg/kg.s-1)
    65       REAL, INTENT(OUT) :: subpdtcloud(ngrid,nlay)     ! tendance temperature due
     64      REAL, INTENT(OUT) :: subpdtcloud     ! tendance temperature due
    6665                                                       !   a la chaleur latente
    6766
    … …  
    8079      INTEGER ig,l,i
    8180     
    82       REAL zq(ngrid,nlay,nq)  ! local value of tracers
    83       REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
    84       REAL zt(ngrid,nlay)       ! local value of temperature
    85       REAL zqsat(ngrid,nlay)    ! saturation
     81      REAL zq(nq)  ! local value of tracers
     82      REAL zq0(nq) ! local initial value of tracers
     83      REAL zt       ! local value of temperature
     84      REAL zqsat_tmp(1)    ! saturation
     85      REAL zqsat    ! saturation
    8686      REAL lw                         !Latent heat of sublimation (J.kg-1)
    8787      REAL cste
    8888      REAL dMice           ! mass of condensed ice
    8989      REAL dMice_hdo       ! mass of condensed HDO ice
    90       REAL alpha(ngrid,nlay) ! HDO equilibrium fractionation coefficient (Saturation=1)
    91       REAL alpha_c(ngrid,nlay) ! HDO real fractionation coefficient
     90      REAL alpha ! HDO equilibrium fractionation coefficient (Saturation=1)
     91      REAL alpha_c ! HDO real fractionation coefficient
    9292!      REAL sumcheck
    9393      REAL*8 ph2o          ! Water vapor partial pressure (Pa)
    … …  
    105105      REAL seq
    106106
    107       REAL rice(ngrid,nlay)      ! Ice mass mean radius (m)
     107      REAL rice      ! Ice mass mean radius (m)
    108108                                 ! (r_c in montmessin_2004)
    109       REAL rhocloud(ngrid,nlay)  ! Cloud density (kg.m-3)
    110       REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
     109      REAL rhocloud  ! Cloud density (kg.m-3)
     110      REAL rdust ! Dust geometric mean radius (m)
    111111
    112112      REAL res      ! Resistance growth
    … …  
    148148
    149149
    150       REAL satubf(ngrid,nlay),satuaf(ngrid,nlay)
    151       REAL res_out(ngrid,nlay)
     150      REAL satubf,satuaf
     151      REAL res_out
    152152 
    153153
    … …  
    236236      cste = 4*pi*rho_ice*microtimestep
    237237
    238       res_out(:,:) = 0
    239       rice(:,:) = 1.e-8
     238      res_out = 0
     239      rice = 1.e-8
    240240
    241241c     Initialize the tendencies
    242       subpdqcloud(1:ngrid,1:nlay,1:nq)=0
    243       subpdtcloud(1:ngrid,1:nlay)=0
     242      subpdqcloud(1:nq)=0
     243      subpdtcloud=0
    244244   
    245      
    246       zt(1:ngrid,1:nlay) =
    247      &      pteff(1:ngrid,1:nlay) +
    248      &      sum_subpdt(1:ngrid,1:nlay) * microtimestep
    249 
    250       zq(1:ngrid,1:nlay,1:nq) =
    251      &      pqeff(1:ngrid,1:nlay,1:nq) +
    252      &      sum_subpdq(1:ngrid,1:nlay,1:nq) * microtimestep
    253      
    254      
    255       WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
    256      &       zq(1:ngrid,1:nlay,1:nq) = 1.e-30
    257 
    258       zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
     245c     temperature and tracers previously incremented here     
     246c     now done outside, in watercloud_mod.F     
     247c      zt =  pteff + sum_subpdt * microtimestep
     248c      zq(1:nq) = pqeff(1:nq) + sum_subpdq(1:nq) * microtimestep
     249     
     250      zq(1:nq)=pq(1:nq)
     251      zt=pt
     252     
     253      WHERE( zq(1:nq) < 1.e-30 )
     254     &       zq(1:nq) = 1.e-30
     255
     256      zq0(1:nq) = zq(1:nq)
    259257     
    260258!=============================================================
    … …  
    265263      dev2 = 1. / ( sqrt(2.) * sigma_ice )
    266264
    267       call watersat(ngrid*nlay,zt,pplay,zqsat)
    268            
     265      call watersat(1,(/zt/),(/pplay/),zqsat_tmp)
     266      zqsat=zqsat_tmp(1)     
    269267      countcells = 0
    270268
    271 c     Main loop over the GCM's grid
    272       DO l=1,nlay
    273         DO ig=1,ngrid
    274 
    275269c       Get the partial pressure of water vapor and its saturation ratio
    276         ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
    277         satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
     270      ph2o = zq(igcm_h2o_vap) * (mmean/18.) * pplay
     271      satu = zq(igcm_h2o_vap) / zqsat
    278272
    279273!=============================================================
    … …  
    281275!=============================================================
    282276
    283        IF ( satu .ge. 1. ) THEN         ! if there is condensation
    284 
    285         call updaterccn(zq(ig,l,igcm_dust_mass),
    286      &          zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
     277      IF ( satu .ge. 1. ) THEN         ! if there is condensation
     278
     279       call updaterccn(zq(igcm_dust_mass),
     280     &         zq(igcm_dust_number),rdust,tauscaling)
    287281
    288282
    289283c       Expand the dust moments into a binned distribution
    290         Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)   + 1.e-30
    291         No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
    292         Rn = rdust(ig,l)
     284        Mo = zq(igcm_dust_mass)* tauscaling   + 1.e-30
     285        No = zq(igcm_dust_number)* tauscaling + 1.e-30
     286        Rn = rdust
    293287        Rn = -log(Rn)
    294288        Rm = Rn - 3. * sigma_ice*sigma_ice 
    … …  
    330324 
    331325c       Get the rates of nucleation
    332         call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
     326        call nuclea(ph2o,zt,satu,n_aer,rate)
    333327       
    334328        dN = 0.
    … …  
    341335
    342336c       Update Dust particles
    343         zq(ig,l,igcm_dust_mass)   =
    344      &  zq(ig,l,igcm_dust_mass)   + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    345         zq(ig,l,igcm_dust_number) =
    346      &  zq(ig,l,igcm_dust_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
     337        zq(igcm_dust_mass)   =
     338     &  zq(igcm_dust_mass)   + dM/ tauscaling !max(tauscaling,1.e-10)
     339        zq(igcm_dust_number) =
     340     &  zq(igcm_dust_number) + dN/ tauscaling !max(tauscaling,1.e-10)
    347341c       Update CCNs
    348         zq(ig,l,igcm_ccn_mass)   =
    349      &  zq(ig,l,igcm_ccn_mass)   - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    350         zq(ig,l,igcm_ccn_number) =
    351      &  zq(ig,l,igcm_ccn_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
     342        zq(igcm_ccn_mass)   =
     343     &  zq(igcm_ccn_mass)   - dM/ tauscaling !max(tauscaling,1.e-10)
     344        zq(igcm_ccn_number) =
     345     &  zq(igcm_ccn_number) - dN/ tauscaling !max(tauscaling,1.e-10)
    352346       
    353         ENDIF ! of is satu >1
     347       ENDIF ! of is satu >1
    354348
    355349!=============================================================
    … …  
    361355c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
    362356
    363        IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
     357       IF ( zq(igcm_ccn_number)*tauscaling.ge. 1.) THEN ! we trigger crystal growth
    364358
    365359     
    366         call updaterice_micro(zq(ig,l,igcm_h2o_ice),
    367      &    zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
    368      &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
    369 
    370         No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
     360        call updaterice_micro(zq(igcm_h2o_ice),
     361     &    zq(igcm_ccn_mass),zq(igcm_ccn_number),
     362     &    tauscaling,rice,rhocloud)
     363
     364        No   = zq(igcm_ccn_number)* tauscaling + 1.e-30
    371365
    372366c       saturation at equilibrium
    373367c       rice should not be too small, otherwise seq value is not valid
    374         seq  = exp(2.*sig(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)*
    375      &             max(rice(ig,l),1.e-7)))
     368        seq  = exp(2.*sig(zt)*mh2o / (rho_ice*rgp*zt*
     369     &             max(rice,1.e-7)))
    376370       
    377371c       get resistance growth
    378         call growthrate(zt(ig,l),pplay(ig,l),
    379      &             real(ph2o/satu),rice(ig,l),res,Dv)
    380 
    381         res_out(ig,l) = res
     372        call growthrate(zt,pplay,
     373     &             real(ph2o/satu),rice,res,Dv)
     374
     375        res_out = res
    382376
    383377ccccccc  implicit scheme of mass growth
    384378
    385379        dMice =
    386      & (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
    387      &   /(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.)
     380     & (zq(igcm_h2o_vap)-seq*zqsat)
     381     &   /(res*zqsat/(cste*No*rice) + 1.)
    388382
    389383
    390384! With the above scheme, dMice cannot be bigger than vapor,
    391385! but can be bigger than all available ice.
    392        dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
    393        dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
    394 
    395        zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
    396        zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
     386       dMice = max(dMice,-zq(igcm_h2o_ice))
     387       dMice = min(dMice,zq(igcm_h2o_vap)) ! this should be useless...
     388
     389       zq(igcm_h2o_ice) = zq(igcm_h2o_ice)+dMice
     390       zq(igcm_h2o_vap) = zq(igcm_h2o_vap)-dMice
    397391
    398392
    … …  
    400394       
    401395       ! latent heat release
    402        lw=(2834.3-0.28*(zt(ig,l)-To)-
    403      &     0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3
    404        subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
     396       lw=(2834.3-0.28*(zt-To)-
     397     &     0.004*(zt-To)*(zt-To))*1.e+3
     398       subpdtcloud= dMice*lw/cpp/microtimestep
    405399         
    406400c Special case of the isotope of water HDO   
    … …  
    411405           if (hdofrac) then
    412406             !! Calculation of the HDO vapor coefficient
    413              Dv_hdo = 1./3. * sqrt( 8*kbz*zt(ig,l)/(pi*mhdo/nav) )
    414      &          * kbz * zt(ig,l) /
    415      &          ( pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo)
     407             Dv_hdo = 1./3. * sqrt( 8*kbz*zt/(pi*mhdo/nav) )
     408     &          * kbz * zt /
     409     &          ( pi * pplay * (molco2+molhdo)*(molco2+molhdo)
    416410     &          * sqrt(1.+mhdo/mco2) )
    417411             !! Calculation of the fractionnation coefficient at equilibrium
    418 c             alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)  ! Merlivat and Nief et al. 1967
    419              alpha = exp(13525./zt(ig,l)**2.-5.59d-2)  ! Lamb et al. 2017
     412c             alpha = exp(16288./zt**2.-9.34d-2)  ! Merlivat and Nief et al. 1967
     413             alpha = exp(13525./zt**2.-5.59d-2)  ! Lamb et al. 2017
    420414             !! Calculation of the 'real' fractionnation coefficient
    421              alpha_c(ig,l) = (alpha(ig,l)*satu)/
    422      &          ( (alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.)
    423 c             alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics
     415             alpha_c = (alpha*satu)/
     416     &          ( (alpha*(Dv/Dv_hdo)*(satu-1.)) + 1.)
     417c             alpha_c = alpha ! to test without the effect of cinetics
    424418           else
    425              alpha_c(ig,l) = 1.d0
     419             alpha_c = 1.d0
    426420           endif
    427            if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then
     421           if (zq0(igcm_h2o_vap).gt.qperemin) then
    428422              dMice_hdo=
    429      &          dMice*alpha_c(ig,l)*
    430      &     ( zq0(ig,l,igcm_hdo_vap)
    431      &             /zq0(ig,l,igcm_h2o_vap) )
     423     &          dMice*alpha_c*
     424     &     ( zq0(igcm_hdo_vap)
     425     &             /zq0(igcm_h2o_vap) )
    432426           else
    433427             dMice_hdo=0.
    … …  
    435429         !! sublimation
    436430         else
    437            if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then
     431           if (zq0(igcm_h2o_ice).gt.qperemin) then
    438432             dMice_hdo=
    439433     &            dMice*
    440      &     ( zq0(ig,l,igcm_hdo_ice)
    441      &             /zq0(ig,l,igcm_h2o_ice) )
     434     &     ( zq0(igcm_hdo_ice)
     435     &             /zq0(igcm_h2o_ice) )
    442436           else
    443437             dMice_hdo=0.
    … …  
    445439         endif !if (dMice.gt.0.0)
    446440
    447        dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice))
    448        dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap))
    449 
    450        zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo
    451        zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo
     441       dMice_hdo = max(dMice_hdo,-zq(igcm_hdo_ice))
     442       dMice_hdo = min(dMice_hdo,zq(igcm_hdo_vap))
     443
     444       zq(igcm_hdo_ice) = zq(igcm_hdo_ice)+dMice_hdo
     445       zq(igcm_hdo_vap) = zq(igcm_hdo_vap)-dMice_hdo
    452446
    453447       endif ! if (hdo)       
    … …  
    459453c         If all the ice particles sublimate, all the condensation
    460454c         nuclei are released:
    461           if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
     455          if (zq(igcm_h2o_ice).le.1.e-28) then
    462456         
    463457c           Water
    464             zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
    465      &                            + zq(ig,l,igcm_h2o_ice)
    466             zq(ig,l,igcm_h2o_ice) = 0.
     458            zq(igcm_h2o_vap) = zq(igcm_h2o_vap)
     459     &                            + zq(igcm_h2o_ice)
     460            zq(igcm_h2o_ice) = 0.
    467461            if (hdo) then
    468               zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)
    469      &                            + zq(ig,l,igcm_hdo_ice)
    470               zq(ig,l,igcm_hdo_ice) = 0.
     462              zq(igcm_hdo_vap) = zq(igcm_hdo_vap)
     463     &                            + zq(igcm_hdo_ice)
     464              zq(igcm_hdo_ice) = 0.
    471465            endif
    472466c           Dust particles
    473             zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
    474      &                              + zq(ig,l,igcm_ccn_mass)
    475             zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
    476      &                                + zq(ig,l,igcm_ccn_number)
     467            zq(igcm_dust_mass) = zq(igcm_dust_mass)
     468     &                              + zq(igcm_ccn_mass)
     469            zq(igcm_dust_number) = zq(igcm_dust_number)
     470     &                                + zq(igcm_ccn_number)
    477471c           CCNs
    478             zq(ig,l,igcm_ccn_mass) = 0.
    479             zq(ig,l,igcm_ccn_number) = 0.
     472            zq(igcm_ccn_mass) = 0.
     473            zq(igcm_ccn_number) = 0.
    480474
    481475          endif
    … …  
    483477          ENDIF !of if Nccn>1
    484478         
    485         ENDDO ! of ig loop
    486       ENDDO ! of nlayer loop
    487      
    488479     
    489480      ! Get cloud tendencies
    490         subpdqcloud(1:ngrid,1:nlay,igcm_h2o_vap) =
    491      &   (zq(1:ngrid,1:nlay,igcm_h2o_vap) -
    492      &    zq0(1:ngrid,1:nlay,igcm_h2o_vap))/microtimestep
    493         subpdqcloud(1:ngrid,1:nlay,igcm_h2o_ice) =
    494      &   (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
    495      &    zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
     481        subpdqcloud(igcm_h2o_vap) =
     482     &   (zq(igcm_h2o_vap) -
     483     &    zq0(igcm_h2o_vap))/microtimestep
     484        subpdqcloud(igcm_h2o_ice) =
     485     &   (zq(igcm_h2o_ice) -
     486     &    zq0(igcm_h2o_ice))/microtimestep
    496487        if (hdo) then
    497           subpdqcloud(1:ngrid,1:nlay,igcm_hdo_vap) =
    498      &     (zq(1:ngrid,1:nlay,igcm_hdo_vap) -
    499      &      zq0(1:ngrid,1:nlay,igcm_hdo_vap))/microtimestep
    500           subpdqcloud(1:ngrid,1:nlay,igcm_hdo_ice) =
    501      &     (zq(1:ngrid,1:nlay,igcm_hdo_ice) -
    502      &      zq0(1:ngrid,1:nlay,igcm_hdo_ice))/microtimestep
     488          subpdqcloud(igcm_hdo_vap) =
     489     &     (zq(igcm_hdo_vap) -
     490     &      zq0(igcm_hdo_vap))/microtimestep
     491          subpdqcloud(igcm_hdo_ice) =
     492     &     (zq(igcm_hdo_ice) -
     493     &      zq0(igcm_hdo_ice))/microtimestep
    503494        endif
    504         subpdqcloud(1:ngrid,1:nlay,igcm_ccn_mass) =
    505      &   (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
    506      &    zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
    507         subpdqcloud(1:ngrid,1:nlay,igcm_ccn_number) =
    508      &   (zq(1:ngrid,1:nlay,igcm_ccn_number) -
    509      &    zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
     495        subpdqcloud(igcm_ccn_mass) =
     496     &   (zq(igcm_ccn_mass) -
     497     &    zq0(igcm_ccn_mass))/microtimestep
     498        subpdqcloud(igcm_ccn_number) =
     499     &   (zq(igcm_ccn_number) -
     500     &    zq0(igcm_ccn_number))/microtimestep
    510501     
    511502      if (scavenging) then
    512503     
    513         subpdqcloud(1:ngrid,1:nlay,igcm_dust_mass) =
    514      &   (zq(1:ngrid,1:nlay,igcm_dust_mass) -
    515      &    zq0(1:ngrid,1:nlay,igcm_dust_mass))/microtimestep
    516         subpdqcloud(1:ngrid,1:nlay,igcm_dust_number) =
    517      &   (zq(1:ngrid,1:nlay,igcm_dust_number) -
    518      &    zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
     504        subpdqcloud(igcm_dust_mass) =
     505     &   (zq(igcm_dust_mass) -
     506     &    zq0(igcm_dust_mass))/microtimestep
     507        subpdqcloud(igcm_dust_number) =
     508     &   (zq(igcm_dust_number) -
     509     &    zq0(igcm_dust_number))/microtimestep
    519510         
    520511      endif
    … …  
    528519     
    529520!       error2d(:) = 0.
    530        DO l=1,nlay
    531        DO ig=1,ngrid
    532 !         error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
    533           satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
    534           satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
    535        ENDDO
    536        ENDDO
    537 
    538       print*, 'count is ',countcells, ' i.e. ',
    539      &     countcells*100/(nlay*ngrid), '% for microphys computation'
     521!         error2d(ig) = max(abs(error_out),error2d(ig))
     522        satubf = zq0(igcm_h2o_vap)/zqsat
     523        satuaf = zq(igcm_h2o_vap)/zqsat
     524
     525c      print*, 'count is ',countcells, ' i.e. ',
     526c     &     countcells*100/(nlay*ngrid), '% for microphys computation'
    540527
    541528#ifndef MESOSCALE
    … …  
    557544!     &                    error_out)
    558545         call write_output("resist","resistance","s/m2",
    559      &                    res_out(:,:))
     546     &                    res_out)
    560547         call write_output("satu_bf","satu before","kg/kg",
    561      &                    satubf(:,:))
     548     &                    satubf)
    562549         call write_output("satu_af","satu after","kg/kg",
    563      &                    satuaf(:,:))
     550     &                    satuaf)
    564551         call write_output("vapbf","h2ovap before","kg/kg",
    565      &                    zq0(:,:,igcm_h2o_vap))
     552     &                    zq0(igcm_h2o_vap))
    566553         call write_output("vapaf","h2ovap after","kg/kg",
    567      &                    zq(:,:,igcm_h2o_vap))
     554     &                    zq(igcm_h2o_vap))
    568555         call write_output("icebf","h2oice before","kg/kg",
    569      &                    zq0(:,:,igcm_h2o_ice))
     556     &                    zq0(igcm_h2o_ice))
    570557         call write_output("iceaf","h2oice after","kg/kg",
    571      &                    zq(:,:,igcm_h2o_ice))
     558     &                    zq(igcm_h2o_ice))
    572559         call write_output("ccnbf","ccn before","/kg",
    573      &                    zq0(:,:,igcm_ccn_number))
     560     &                    zq0(igcm_ccn_number))
    574561         call write_output("ccnaf","ccn after","/kg",
    575      &                    zq(:,:,igcm_ccn_number))
     562     &                    zq(igcm_ccn_number))
    576563c         call write_output("growthrate","growth rate","m^2/s",
    577564c     &                    gr_out)
    … …  
    583570c     &                    dN_out)
    584571         call write_output("zqsat","p vap sat","kg/kg",
    585      &                    zqsat(:,:))
     572     &                    zqsat)
    586573!         call write_output("satu","ratio saturation","",
    587574!     &                    satu_out(:,:))
    588575         call write_output("rice","ice radius","m",
    589      &                    rice(:,:))
     576     &                    rice)
    590577!         call write_output("rdust_sca","rdust","m",
    591578!     &                    rdust)

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

    r2934 r2966  
    2727     &                      qperemin
    2828      use dimradmars_mod, only: naerkind
     29      use conc_mod, only: mmean
    2930      use write_output_mod, only: write_output
    3031      IMPLICIT NONE
    … …  
    4546c
    4647c  2004 - 2012
     48c
     49c     2023: J. Naar, adding different subtimestep on each grid cell
     50c          plus not doing microphysics if no water present
     51c          plus simpleclouds no longer in the loop for microphysics
    4752c=======================================================================
    4853
    … …  
    105110
    106111      ! global tendency (clouds+physics)
     112      ! JN : keeping this for simpleclouds scheme
    107113      REAL sum_subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
    108114      REAL sum_subpdt(ngrid,nlay)         ! cf. pdtcloud
    … …  
    137143     
    138144      !$OMP THREADPRIVATE(flagcloud)
     145 
     146! Scheme for adaptative timestep J. Naar 2023
     147c      LOGICAL :: computed_micro(ngrid,nlay) ! Check if microphy was done in this cell
     148      REAL :: computed_micro(ngrid,nlay) ! Check if microphy was done inthis cell (logical)
     149      REAL :: zt_micro(ngrid,nlay) ! Temperature during microphysics (K)
     150      REAL :: zq_micro(ngrid,nlay,nq) ! Tracers during microphysics (kg/kg)
     151      REAL :: zqsat_micro(ngrid,nlay) ! Theoritical q(h2o_vap) at saturation (kg/kg)
     152      INTEGER :: zimicro(ngrid,nlay)    ! Number of ptimestep divisions
     153      REAL :: zpotcond_inst(ngrid,nlay) ! condensable water at the beginning of physics (kg/kg)
     154      REAL :: zpotcond_full(ngrid,nlay) ! condensable water with integrated tendancies (kg/kg)
     155      REAL :: zpotcond(ngrid,nlay) ! maximal condensable water, used to
     156compute adaptative subdivision of ptimestep
     157
    139158
    140159c    ** un petit test de coherence
    … …  
    282301c------------------------------------------------------------------
    283302      rhocloud(1:ngrid,1:nlay) = rho_dust
    284       DO microstep=1,imicro
    285      
     303
     304
     305c     Initialisation of all the stuff JN2023
     306c      computed_micro(1:ngrid,1:nlay)=.false.
     307      computed_micro(1:ngrid,1:nlay)=0.
     308      zt_micro(:,:)=pt(:,:)
     309      zq_micro(:,:,:)=pq(:,:,:)
     310      zq_micro(:,:,:)=pq(:,:,:)
     311      call watersat(ngrid*nlay,zt_micro,pplay,zqsat_micro)
     312      zpotcond_inst=zq_micro(:,:,igcm_h2o_vap) - zqsat_micro
     313      call watersat(ngrid*nlay,zt_micro+pdt*ptimestep,pplay,zqsat_micro)
     314      zpotcond_full=(zq_micro(:,:,igcm_h2o_vap)+
     315     &             pdq(:,:,igcm_h2o_vap)*ptimestep) - zqsat_micro
     316      zimicro(1:ngrid,1:nlay)=imicro
     317      if (cloud_adapt_ts) then
     318              call adapt_imicro(ptimestep,zpotcond(ig,l),
     319     $                   zimicro(ig,l))
     320      endif! (cloud_adapt_ts) then
     321      DO l=1,nlay
     322        DO ig=1,ngrid
     323c         Start by computing the condensable water vapor amount
     324          if (zpotcond_full(ig,l).gt.0.) then
     325            zpotcond(ig,l)=max(zpotcond_inst(ig,l),zpotcond_full(ig,l))
     326          else if (zpotcond_full(ig,l).le.0.) then
     327            zpotcond(ig,l)=min(zpotcond_inst(ig,l),zpotcond_full(ig,l))
     328          endif! (zpotcond_full.gt.0.) then
     329          microtimestep=ptimestep/real(zimicro(ig,l))
     330c         Check if microphysics is even needed, that is if enough action
     331c         is happening water-wise
     332          if ((pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
     333     &      .gt.1e-22) .or. (abs(zpotcond(ig,l)).gt.1e-22)) then
     334c         Eventuellement sortir simpleclouds de la boucle egalement
     335          computed_micro(ig,l)=1.
     336          DO microstep=1,zimicro(ig,l)
     337     
     338
     339c JN : incrementing after main microphysics scheme
     340c Previously we were incrementing tendencies, we now
     341c increment tracers and temperature directly
     342c We are thus starting at the end of the first iteration
     343c
    286344c-------------------------------------------------------------------
    287 c   1.  Tendencies:
    288 c------------------
    289 
    290 
    291 c------ Temperature tendency subpdt
    292         ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
    293         ! If imicro=1 subpdt is the same as pdt
    294         DO l=1,nlay
    295           DO ig=1,ngrid
    296              sum_subpdt(ig,l) = sum_subpdt(ig,l)
    297      &        + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
    298           ENDDO
    299         ENDDO
    300 c------ Tracers tendencies subpdq are additionned
    301 c------ At each micro timestep we add pdq in order to have a stepped entry
    302         IF (microphys) THEN
    303           DO l=1,nlay
    304             DO ig=1,ngrid
    305               sum_subpdq(ig,l,igcm_dust_mass) =
    306      &            sum_subpdq(ig,l,igcm_dust_mass)
    307      &          + pdq(ig,l,igcm_dust_mass)
    308               sum_subpdq(ig,l,igcm_dust_number) =
    309      &            sum_subpdq(ig,l,igcm_dust_number)
    310      &          + pdq(ig,l,igcm_dust_number)
    311               sum_subpdq(ig,l,igcm_ccn_mass) =
    312      &            sum_subpdq(ig,l,igcm_ccn_mass)
    313      &          + pdq(ig,l,igcm_ccn_mass)
    314               sum_subpdq(ig,l,igcm_ccn_number) =
    315      &            sum_subpdq(ig,l,igcm_ccn_number)
    316      &          + pdq(ig,l,igcm_ccn_number)
    317             ENDDO
    318           ENDDO
    319         ENDIF
    320         DO l=1,nlay
    321           DO ig=1,ngrid
    322             sum_subpdq(ig,l,igcm_h2o_ice) =
    323      &          sum_subpdq(ig,l,igcm_h2o_ice)
    324      &        + pdq(ig,l,igcm_h2o_ice)
    325             sum_subpdq(ig,l,igcm_h2o_vap) =
    326      &          sum_subpdq(ig,l,igcm_h2o_vap)
    327      &        + pdq(ig,l,igcm_h2o_vap)
    328             IF (hdo) THEN
    329             sum_subpdq(ig,l,igcm_hdo_ice) =
    330      &          sum_subpdq(ig,l,igcm_hdo_ice)
    331      &        + pdq(ig,l,igcm_hdo_ice)
    332             sum_subpdq(ig,l,igcm_hdo_vap) =
    333      &          sum_subpdq(ig,l,igcm_hdo_vap)
    334      &        + pdq(ig,l,igcm_hdo_vap)
    335             ENDIF !hdo
    336           ENDDO
    337         ENDDO     
    338        
    339 c-------------------------------------------------------------------
    340 c   2.  Main call to the different cloud schemes:
     345c   1.  Main call to the different cloud schemes:
    341346c------------------------------------------------
    342347        IF (microphys) THEN
    343            CALL improvedclouds(ngrid,nlay,microtimestep,
    344      &             pplay,pteff,sum_subpdt,
    345      &             pqeff,sum_subpdq,subpdqcloud,subpdtcloud,
    346      &             nq,tauscaling)
     348           CALL improvedclouds(microtimestep,
     349     &          pplay(ig,l),zt_micro(ig,l),
     350     &          zq_micro(ig,l,:),subpdqcloud(ig,l,:),
     351     &          subpdtcloud(ig,l),nq,tauscaling(ig),mmean(ig,l))
    347352
    348353        ELSE
     354c Simpleclouds should maybe be taken out and put in a specific loop ?
    349355           CALL simpleclouds(ngrid,nlay,microtimestep,
    350356     &             pplay,pzlay,pteff,sum_subpdt,
    … …  
    354360
    355361c-------------------------------------------------------------------
    356 c   3.  Updating tendencies after cloud scheme:
     362c   2.  Updating tracers and temperature after cloud scheme:
    357363c-----------------------------------------------
    358364
    359365        IF (microphys) THEN
    360           DO l=1,nlay
    361             DO ig=1,ngrid
    362               sum_subpdq(ig,l,igcm_dust_mass) =
    363      &            sum_subpdq(ig,l,igcm_dust_mass)
    364      &          + subpdqcloud(ig,l,igcm_dust_mass)
    365               sum_subpdq(ig,l,igcm_dust_number) =
    366      &            sum_subpdq(ig,l,igcm_dust_number)
    367      &          + subpdqcloud(ig,l,igcm_dust_number)
    368               sum_subpdq(ig,l,igcm_ccn_mass) =
    369      &            sum_subpdq(ig,l,igcm_ccn_mass)
    370      &          + subpdqcloud(ig,l,igcm_ccn_mass)
    371               sum_subpdq(ig,l,igcm_ccn_number) =
    372      &            sum_subpdq(ig,l,igcm_ccn_number)
    373      &          + subpdqcloud(ig,l,igcm_ccn_number)
    374             ENDDO
    375           ENDDO
     366              zq_micro(ig,l,igcm_dust_mass) =
     367     &         zq_micro(ig,l,igcm_dust_mass)+(pdq(ig,l,igcm_dust_mass)
     368     &         +subpdqcloud(ig,l,igcm_dust_mass))*microtimestep
     369              zq_micro(ig,l,igcm_dust_number) =
     370     &         zq_micro(ig,l,igcm_dust_number)
     371     &         +(pdq(ig,l,igcm_dust_number)
     372     &         + subpdqcloud(ig,l,igcm_dust_number))*microtimestep
     373              zq_micro(ig,l,igcm_ccn_mass) =
     374     &         zq_micro(ig,l,igcm_ccn_mass) +
     375     &         (pdq(ig,l,igcm_ccn_mass)
     376     &         +subpdqcloud(ig,l,igcm_ccn_mass))*microtimestep
     377              zq_micro(ig,l,igcm_ccn_number) =
     378     &          zq_micro(ig,l,igcm_ccn_number) +
     379     &         (pdq(ig,l,igcm_ccn_number)
     380     &          + subpdqcloud(ig,l,igcm_ccn_number))*microtimestep
    376381        ENDIF
    377         DO l=1,nlay
    378           DO ig=1,ngrid
    379             sum_subpdq(ig,l,igcm_h2o_ice) =
    380      &          sum_subpdq(ig,l,igcm_h2o_ice)
    381      &        + subpdqcloud(ig,l,igcm_h2o_ice)
    382             sum_subpdq(ig,l,igcm_h2o_vap) =
    383      &          sum_subpdq(ig,l,igcm_h2o_vap)
    384      &        + subpdqcloud(ig,l,igcm_h2o_vap)
     382            zq_micro(ig,l,igcm_h2o_ice) =
     383     &       zq_micro(ig,l,igcm_h2o_ice)+
     384     &         (pdq(ig,l,igcm_h2o_ice)
     385     &        + subpdqcloud(ig,l,igcm_h2o_ice))*microtimestep
     386            zq_micro(ig,l,igcm_h2o_vap) =
     387     &       zq_micro(ig,l,igcm_h2o_vap)+
     388     &         (pdq(ig,l,igcm_h2o_vap)
     389     &        + subpdqcloud(ig,l,igcm_h2o_vap))*microtimestep
    385390
    386391            IF (hdo) THEN
    387             sum_subpdq(ig,l,igcm_hdo_ice) =
    388      &          sum_subpdq(ig,l,igcm_hdo_ice)
    389      &        + subpdqcloud(ig,l,igcm_hdo_ice)
    390             sum_subpdq(ig,l,igcm_hdo_vap) =
    391      &          sum_subpdq(ig,l,igcm_hdo_vap)
    392      &        + subpdqcloud(ig,l,igcm_hdo_vap)
     392            zq_micro(ig,l,igcm_hdo_ice) =
     393     &       zq_micro(ig,l,igcm_hdo_ice)+
     394     &         (pdq(ig,l,igcm_hdo_ice)
     395     &        + subpdqcloud(ig,l,igcm_hdo_ice))*microtimestep
     396            zq_micro(ig,l,igcm_hdo_vap) =
     397     &       zq_micro(ig,l,igcm_hdo_vap)+
     398     &         (pdq(ig,l,igcm_hdo_vap)
     399     &        + subpdqcloud(ig,l,igcm_hdo_vap))*microtimestep
    393400            ENDIF ! hdo
    394401
    395           ENDDO
    396         ENDDO
    397         
     402c  Could also set subpdtcloud to 0 if not activice to make it simpler
     403            zt_micro(ig,l) = zt_micro(ig,l)+
     404     &           pdt(ig,l)*microtimestep
    398405        IF (activice) THEN
    399           DO l=1,nlay
    400             DO ig=1,ngrid
    401               sum_subpdt(ig,l) =
    402      &            sum_subpdt(ig,l) + subpdtcloud(ig,l)
    403             ENDDO
    404           ENDDO
     406              zt_micro(ig,l) = zt_micro(ig,l)+
     407     &           subpdtcloud(ig,l)*microtimestep
    405408        ENDIF
    406 
    407409!      !! Example of how to use writediagmicrofi useful to
    408410!      !! get outputs at each microphysical sub-timestep (better to be used in 1D)
    … …  
    411413!     &      'subpdtcloud','K/s',1,subpdtcloud(:,:))     
    412414 
    413       ENDDO ! of DO microstep=1,imicro
    414      
     415          ENDDO ! of DO microstep=1,imicro
     416          endif! (zq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
     417!     &      .gt.1e-22).or.(abs(zpotcond).gt.1e-22) then
     418        ENDDO ! ig=1,ngrid
     419      ENDDO ! l=1,nlay
     420
     421     
     422c------ Useful outputs to check how it went
     423      call write_output("computed_micro","computed_micro "//
     424     &   "after microphysics","logical",computed_micro(:,:))
     425      call write_output("zimicro","Used number of subtimestep "//
     426     &   "in cloud microphysics","integer",real(zimicro(:,:)))
    415427c-------------------------------------------------------------------
    416 c   6.  Compute final tendencies after time loop:
     428c   3.  Compute final tendencies after time loop:
    417429c------------------------------------------------
    418430
    … …  
    420432       DO l=1,nlay
    421433         DO ig=1,ngrid
    422              pdtcloud(ig,l) =
    423      &         sum_subpdt(ig,l)/real(imicro)-pdt(ig,l)
    424           ENDDO
     434         pdtcloud(ig,l) = -pdt(ig,l)+
     435     &        (zt_micro(ig,l)-pt(ig,l)) / ptimestep
     436         ENDDO
    425437       ENDDO
    426        
     438
    427439c------ Tracers tendencies pdqcloud
    428440       DO l=1,nlay
    429441         DO ig=1,ngrid
    430             pdqcloud(ig,l,igcm_h2o_ice) =
    431      &        sum_subpdq(ig,l,igcm_h2o_ice)/real(imicro)
    432      &       - pdq(ig,l,igcm_h2o_ice)
    433             pdqcloud(ig,l,igcm_h2o_vap) =
    434      &        sum_subpdq(ig,l,igcm_h2o_vap)/real(imicro)
    435      &       - pdq(ig,l,igcm_h2o_vap)
     442            pdqcloud(ig,l,igcm_h2o_ice) =
     443     &        -pdq(ig,l,igcm_h2o_ice)+
     444     &       (zq_micro(ig,l,igcm_h2o_ice) -
     445     &        pq(ig,l,igcm_h2o_ice)) / ptimestep
     446            pdqcloud(ig,l,igcm_h2o_vap) =
     447     &        -pdq(ig,l,igcm_h2o_vap)+
     448     &       (zq_micro(ig,l,igcm_h2o_vap) -
     449     &        pq(ig,l,igcm_h2o_vap)) / ptimestep
    436450            IF (hdo) THEN
    437             pdqcloud(ig,l,igcm_hdo_ice) =
    438      &        sum_subpdq(ig,l,igcm_hdo_ice)/real(imicro)
    439      &       - pdq(ig,l,igcm_hdo_ice)
    440             pdqcloud(ig,l,igcm_hdo_vap) =
    441      &        sum_subpdq(ig,l,igcm_hdo_vap)/real(imicro)
    442      &       - pdq(ig,l,igcm_hdo_vap)
     451            pdqcloud(ig,l,igcm_hdo_ice) =
     452     &        -pdq(ig,l,igcm_hdo_ice)+
     453     &       (zq_micro(ig,l,igcm_hdo_ice) -
     454     &        pq(ig,l,igcm_hdo_ice)) / ptimestep
     455            pdqcloud(ig,l,igcm_hdo_vap) =
     456     &        -pdq(ig,l,igcm_hdo_vap)+
     457     &       (zq_micro(ig,l,igcm_hdo_vap) -
     458     &        pq(ig,l,igcm_hdo_vap)) / ptimestep
    443459            ENDIF !hdo
    444460         ENDDO
    445        ENDDO
    446        
     461       ENDDO       
     462
    447463       IF(microphys) THEN
    448464        DO l=1,nlay
    449465         DO ig=1,ngrid
    450             pdqcloud(ig,l,igcm_ccn_mass) =
    451      &        sum_subpdq(ig,l,igcm_ccn_mass)/real(imicro)
    452      &       - pdq(ig,l,igcm_ccn_mass)
    453             pdqcloud(ig,l,igcm_ccn_number) =
    454      &        sum_subpdq(ig,l,igcm_ccn_number)/real(imicro)
    455      &       - pdq(ig,l,igcm_ccn_number)
     466            pdqcloud(ig,l,igcm_ccn_mass) =
     467     &        -pdq(ig,l,igcm_ccn_mass)+
     468     &       (zq_micro(ig,l,igcm_ccn_mass) -
     469     &        pq(ig,l,igcm_ccn_mass)) / ptimestep
     470            pdqcloud(ig,l,igcm_ccn_number) =
     471     &        -pdq(ig,l,igcm_ccn_number)+
     472     &       (zq_micro(ig,l,igcm_ccn_number) -
     473     &        pq(ig,l,igcm_ccn_number)) / ptimestep
    456474         ENDDO
    457475        ENDDO
    458476       ENDIF
    459        
     477
    460478       IF(scavenging) THEN
    461479        DO l=1,nlay
    462480         DO ig=1,ngrid
    463             pdqcloud(ig,l,igcm_dust_mass) =
    464      &        sum_subpdq(ig,l,igcm_dust_mass)/real(imicro)
    465      &       - pdq(ig,l,igcm_dust_mass)
    466             pdqcloud(ig,l,igcm_dust_number) =
    467      &        sum_subpdq(ig,l,igcm_dust_number)/real(imicro)
    468      &       - pdq(ig,l,igcm_dust_number)
     481            pdqcloud(ig,l,igcm_dust_mass) =
     482     &        -pdq(ig,l,igcm_dust_mass)+
     483     &       (zq_micro(ig,l,igcm_dust_mass) -
     484     &        pq(ig,l,igcm_dust_mass)) / ptimestep
     485            pdqcloud(ig,l,igcm_dust_number) =
     486     &        -pdq(ig,l,igcm_dust_number)+
     487     &       (zq_micro(ig,l,igcm_dust_number) -
     488     &        pq(ig,l,igcm_dust_number)) / ptimestep
    469489         ENDDO
    470490        ENDDO
    … …  
    764784       end subroutine end_watercloud_mod
    765785
     786
     787
     788       SUBROUTINE adapt_imicro(ptimestep,potcond,
     789     $                     zimicro)
     790
     791c Pas de temps adaptatif pour les nuages
     792
     793      real,intent(in) :: ptimestep ! total duration of physics (sec)
     794      real,intent(in) :: potcond ! total duration of physics (sec)
     795      real :: alpha, beta ! total duration of physics (sec)
     796      integer,intent(out) :: zimicro ! number of ptimestep division
     797
     798c       zimicro = ptimestep*alpha*potcond**beta
     799       zimicro = 30
     800
     801       END SUBROUTINE adapt_imicro
     802
     803
    766804      END MODULE watercloud_mod
Note: See TracChangeset for help on using the changeset viewer.