Changeset 3939


Ignore:
Timestamp:
Oct 29, 2025, 10:52:34 AM (29 hours ago)
Author:
jmauxion
Message:

MARS PCM:
A stable version of the watercloud microphysic scheme with adaptative
timestep. Now it should run 1 year at high obliquity with no crash.
JM

Location:
trunk/LMDZ.MARS
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/changelog.txt

    r3924 r3939  
    49814981Formating "improveclouds" into F90 and simplifying some instructions.
    49824982
    4983 == 06/09/2025 == JM
     4983== 06/10/2025 == JM
    49844984Adding a slow_diagfi flag to the run.def file for 1D model only. When False, the netcdf
    49854985file is opened/closed once, thus saving significant computing time. When true,
    49864986the opening frequency is at output frequency (required in debug mode).
    49874987
    4988 == 06/09/2025 == JM
     4988== 06/10/2025 == JM
    49894989Fixing revision 3923: adding OMP_THREADPRIVATE to new save variables.
     4990
     4991== 29/10/2025 == JM
     4992A stable version of the watercloud microphysic scheme with adaptative
     4993timestep. Now it should run 1 year at high obliquity with no crash.

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

    r3919 r3939  
    11MODULE improvedclouds_mod
    22
    3 implicit none
     3IMPLICIT NONE
    44
    55CONTAINS
    … …  
    99SUBROUTINE improvedclouds(ngrid,nlay,ptimestep,pplay,pt,pdt,pq,pdq,nq,tauscaling,imicro,zt,zq)
    1010
    11 use updaterad,        only: updaterice_micro, updaterccn
    12 use watersat_mod,     only: watersat
    13 use tracer_mod,       only: rho_ice, nuice_sed, igcm_h2o_vap, igcm_h2o_ice, igcm_dust_mass, &
    14                             igcm_dust_number, igcm_ccn_mass, igcm_ccn_number, igcm_hdo_vap, &
    15                             igcm_hdo_ice, qperemin
    16 use conc_mod,         only: mmean
    17 use comcstfi_h,       only: pi, cpp
    18 use microphys_h,      only: nbin_cld, rad_cld, mteta, kbz, nav, rgp
    19 use microphys_h,      only: mco2, vo1, mh2o, mhdo, molco2, molhdo, To
    20 use nuclea_mod,       only: nuclea
    21 use sig_h2o_mod,      only: sig_h2o
    22 use growthrate_mod,   only: growthrate
     11use updaterad, only: updaterice_micro, updaterccn
     12use watersat_mod, only: watersat
     13use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap, igcm_h2o_ice, igcm_dust_mass, &
     14                      igcm_dust_number, igcm_ccn_mass, igcm_ccn_number, igcm_hdo_vap, &
     15                      igcm_hdo_ice,qperemin
     16use conc_mod, only: mmean
     17use comcstfi_h, only: pi, cpp
     18use microphys_h, only: nbin_cld, rad_cld, mteta, kbz, nav, rgp
     19use microphys_h, only: mco2, vo1, mh2o, mhdo, molco2, molhdo, To
     20use nuclea_mod, only: nuclea
     21use sig_h2o_mod, only: sig_h2o
     22use growthrate_mod, only: growthrate
    2323use write_output_mod, only: write_output
    24 use callkeys_mod,     only: activice, scavenging, cloud_adapt_ts, hdo, hdofrac
     24use callkeys_mod, only: activice, scavenging, cloud_adapt_ts, hdo, hdofrac
    2525
    2626implicit none
    … …  
    4747!           J. Naar, adaptative subtimestep now done here (June 2023)
    4848!------------------------------------------------------------------
     49
     50!------------------------------------------------------------------
    4951!     Inputs/outputs:
    50 INTEGER,                        INTENT(IN) :: ngrid,nlay
    51 INTEGER,                        INTENT(IN) :: nq        ! nombre de traceurs
    52 REAL,                           INTENT(IN) :: ptimestep ! pas de temps physique (s)
    53 REAL, dimension(ngrid,nlay),    INTENT(IN) :: pplay      ! pression au milieu des couches (Pa)
    54 REAL, dimension(ngrid,nlay),    INTENT(IN) :: pt         ! temperature at the middle of the layers (K)
    55 REAL, dimension(ngrid,nlay),    INTENT(IN) :: pdt        ! temperature tendency (K/s)
    56 REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pq         ! tracer (kg/kg)
    57 REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pdq        ! tracer tendency (kg/kg/s)
    58 REAL, dimension(ngrid),         INTENT(IN) :: tauscaling ! Convertion factor for qdust and Ndust
    59 INTEGER,                        INTENT(IN) :: imicro    ! nb. microphy calls(retrocompatibility)
    60 
    61 REAL, dimension(ngrid,nlay,nq), INTENT(OUT) :: zq ! tracers post microphy (kg/kg)
    62 REAL, dimension(ngrid,nlay),    INTENT(OUT) :: zt ! temperature post microphy (K)
     52INTEGER, INTENT(IN) :: ngrid,nlay
     53INTEGER, INTENT(IN) :: nq ! nombre de traceurs
     54REAL, INTENT(IN) :: ptimestep ! pas de temps physique (s)
     55REAL, dimension(ngrid,nlay), INTENT(IN) :: pplay ! pression au milieu des couches (Pa)           
     56REAL, dimension(ngrid,nlay), INTENT(IN) :: pt ! temperature at the middle of the layers (K)
     57REAL, dimension(ngrid,nlay), INTENT(IN) :: pdt ! temperature tendency (K/s)
     58REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pq ! tracer (kg/kg)
     59REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pdq ! tracer tendency (kg/kg/s)
     60REAL, dimension(ngrid), INTENT(IN) :: tauscaling ! Convertion factor for qdust and Ndust
     61INTEGER, INTENT(IN) :: imicro ! nb. microphy calls(retrocompatibility)
     62
     63REAL, dimension(ngrid,nlay,nq), INTENT(OUT) :: zq ! tracers post microphy (kg/kg) 
     64REAL, dimension(ngrid,nlay), INTENT(OUT) :: zt ! temperature post microphy (K)
    6365
    6466!------------------------------------------------------------------
    6567!     Local variables:
    66 LOGICAL, SAVE               ::  firstcall = .true.
     68LOGICAL, SAVE ::  firstcall = .true.
    6769!$OMP THREADPRIVATE(firstcall)
    68 REAL*8                      :: derf ! Error function
     70REAL*8 :: derf ! Error function
    6971!external derf
    70 INTEGER                     :: ig, l, i
    71 REAL, dimension(ngrid,nlay) :: zqsat ! Saturation
    72 REAL                        :: lw    ! Latent heat of sublimation (J.kg-1)
    73 REAL                        :: cste
    74 REAL                        :: dMice     ! mass of condensed ice
    75 REAL                        :: dMice_hdo ! mass of condensed HDO ice
    76 REAL, dimension(ngrid,nlay) :: alpha     ! HDO equilibrium fractionation coefficient (Saturation=1)
    77 REAL, dimension(ngrid,nlay) :: alpha_c   ! HDO real fractionation coefficient
     72INTEGER :: ig,l,i
     73REAL, dimension(ngrid,nlay) :: zqsat ! saturation
     74REAL :: lw ! Latent heat of sublimation (J.kg-1)
     75REAL :: cste
     76REAL :: dMice ! mass of condensed ice
     77REAL :: dMicetot ! JN
     78REAL :: dMice_hdo ! mass of condensed HDO ice
     79REAL, dimension(ngrid,nlay) :: alpha ! HDO equilibrium fractionation coefficient (Saturation=1)
     80REAL, dimension(ngrid,nlay) :: alpha_c ! HDO real fractionation coefficient
    7881! REAL :: sumcheck
    79 REAL*8                      :: ph2o ! Water vapor partial pressure (Pa)
    80 REAL*8                      :: satu ! Water vapor saturation ratio over ice
    81 REAL*8                      :: Mo, No, Rn, Rm, dev2, n_derf, m_derf
     82REAL*8 :: ph2o ! Water vapor partial pressure (Pa)
     83REAL*8 :: satu ! Water vapor saturation ratio over ice
     84REAL*8 :: Mo,No, Rn, Rm, dev2, n_derf, m_derf
    8285REAL*8, dimension(nbin_cld) :: n_aer ! number conc. of particle/each size bin
    8386REAL*8, dimension(nbin_cld) :: m_aer ! mass mixing ratio of particle/each size bin
    8487REAL :: dN, dM, seq
    85 REAL, dimension(nbin_cld)   :: rate ! Nucleation rate
    86 REAL, dimension(ngrid,nlay) :: rice ! Ice mass mean radius (m)
    87                                     ! (r_c in montmessin_2004)
    88 REAL, dimension(ngrid,nlay) :: rhocloud   ! Cloud density (kg.m-3)
    89 REAL, dimension(ngrid,nlay) :: rdust      ! Dust geometric mean radius (m)
    90 REAL                        :: res        ! Resistance growth
    91 REAL                        :: Dv, Dv_hdo ! Water/HDO vapor diffusion coefficient
    92 
    93 ! Parameters of the size discretization used by the microphysical scheme
    94 DOUBLE PRECISION, PARAMETER                     :: rmin_cld = 0.1e-6     ! Minimum radius (m)
    95 DOUBLE PRECISION, PARAMETER                     :: rmax_cld = 10.e-6     ! Maximum radius (m)
    96 DOUBLE PRECISION, PARAMETER                     :: rbmin_cld = 0.0001e-6 ! Minimum boundary radius (m)
    97 DOUBLE PRECISION, PARAMETER                     :: rbmax_cld = 1.e-2     ! Maximum boundary radius (m)
    98 DOUBLE PRECISION                                :: vrat_cld ! Volume ratio
    99 DOUBLE PRECISION, dimension(nbin_cld + 1), SAVE :: rb_cld   ! Boundary values of each rad_cld bin (m)
     88REAL, dimension(nbin_cld) :: rate ! nucleation rate
     89REAL, dimension(ngrid,nlay) :: rice ! Ice mass mean radius (m) (r_c in montmessin_2004)
     90REAL, dimension(ngrid,nlay) :: rhocloud ! Cloud density (kg.m-3)
     91REAL, dimension(ngrid,nlay) :: rdust ! Dust geometric mean radius (m)
     92REAL :: res ! Resistance growth
     93REAL :: Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient
     94
     95!     Parameters of the size discretization used by the microphysical scheme
     96DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
     97DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
     98DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6 ! Minimum boundary radius (m)
     99DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
     100DOUBLE PRECISION :: vrat_cld ! Volume ratio
     101DOUBLE PRECISION, dimension(nbin_cld+1), SAVE :: rb_cld ! boundary values of each rad_cld bin (m)
    100102!$OMP THREADPRIVATE(rb_cld)
    101 DOUBLE PRECISION, dimension(nbin_cld) :: dr_cld    ! Width of each rad_cld bin (m)
    102 DOUBLE PRECISION, dimension(nbin_cld) :: vol_cld   ! Particle volume for each bin (m3)
    103 REAL, SAVE                            :: sigma_ice ! Variance of the ice and CCN distributions
     103DOUBLE PRECISION, dimension(nbin_cld) :: dr_cld ! width of each rad_cld bin (m)
     104DOUBLE PRECISION, dimension(nbin_cld) :: vol_cld ! particle volume for each bin (m3)
     105REAL, SAVE :: sigma_ice ! Variance of the ice and CCN distributions
    104106!$OMP THREADPRIVATE(sigma_ice)
    105107
    106 ! JN: used in subtimestep
    107 REAL                              :: microtimestep ! Subdivision of ptimestep (s)
    108 REAL,    dimension(ngrid,nlay)    :: subpdtcloud   ! Temperature variation due to latent heat
    109 REAL,    dimension(ngrid,nlay,nq) :: zq0           ! Local initial value of tracers
    110 !REAL,    dimension(ngrid,nlay,nq) :: zt0           ! Local initial value of temperature
    111 INTEGER, dimension(ngrid,nlay)    :: zimicro       ! Subdivision of ptimestep
    112 INTEGER, dimension(ngrid,nlay)    :: count_micro   ! Number of microphys calls
    113 REAL,    dimension(ngrid,nlay)    :: zpotcond_inst ! Condensable water at the beginning of physics (kg/kg)
    114 REAL,    dimension(ngrid,nlay)    :: zpotcond_full ! Condensable water with integrated tendancies (kg/kg)
    115 REAL,    dimension(ngrid,nlay)    :: zpotcond      ! Maximal condensable water (previous two)
    116 REAL,    dimension(1)             :: zqsat_tmp     ! Maximal condensable water (previous two)
    117 REAL                              :: spenttime     ! Time spent in while loop
    118 REAL                              :: zdq           ! Used to compute adaptative timestep
    119 LOGICAL                           :: ending_ts     ! Condition to end while loop
    120 
    121 !----------------------------------
     108!----------------------------------     
     109! JN : used in subtimestep
     110REAL :: microtimestep! subdivision of ptimestep (s)
     111REAL, dimension(ngrid,nlay) :: subpdtcloud ! Temperature variation due to latent heat
     112REAL, dimension(ngrid,nlay,nq) :: zq0 ! local initial value of tracers
     113! REAL, dimension(ngrid,nlay,nq) :: zt0 ! local initial value of temperature
     114INTEGER, dimension(ngrid,nlay) :: zimicro ! Subdivision of ptimestep
     115INTEGER, dimension(ngrid,nlay) :: count_micro ! Number of microphys calls
     116REAL, dimension(ngrid,nlay) :: zpotcond_inst ! condensable water at the beginning of physics (kg/kg)
     117REAL, dimension(ngrid,nlay) :: zpotcond_full ! condensable water with integrated tendancies (kg/kg)
     118REAL, dimension(ngrid,nlay) :: zpotcond ! maximal condensable water (previous two)
     119REAL, dimension(1) :: zqsat_tmp ! maximal condensable water (previous two)
     120REAL :: spenttime ! time spent in while loop
     121REAL :: zdq ! used to compute adaptative timestep
     122LOGICAL :: ending_ts ! Condition to end while loop
     123
     124!----------------------------------     
    122125! TESTS
    123 INTEGER       :: countcells
    124 LOGICAL, SAVE :: test_flag ! Flag for test/debuging outputs
     126INTEGER :: countcells
     127LOGICAL, SAVE :: test_flag    ! flag for test/debuging outputs
    125128!$OMP THREADPRIVATE(test_flag)
    126 REAL, dimension(ngrid,nlay) :: satubf, satuaf, res_out
     129REAL, dimension(ngrid,nlay) :: satubf,satuaf, res_out
    127130
    128131!------------------------------------------------------------------
    … …  
    134137!=============================================================
    135138! rad_cld is the primary radius grid used for microphysics computation.
    136 ! The grid spacing is computed assuming a constant volume ratio between two consecutive bins; i.e. vrat_cld.
    137 ! vrat_cld is determined from the boundary values of the size grid: rmin_cld and rmax_cld.
     139! The grid spacing is computed assuming a constant volume ratio
     140! between two consecutive bins; i.e. vrat_cld.
     141! vrat_cld is determined from the boundary values of the size grid:
     142! rmin_cld and rmax_cld.
    138143! The rb_cld array contains the boundary values of each rad_cld bin.
    139144! dr_cld is the width of each rad_cld bin.
    140145
    141 ! Volume ratio between two adjacent bins
    142 !    vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
    143 !    vrat_cld = exp(vrat_cld)
    144     vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
    145     vrat_cld = exp(vrat_cld)
    146     write(*,*) "vrat_cld", vrat_cld
    147 
    148     rb_cld(1)  = rbmin_cld
    149     rad_cld(1) = rmin_cld
    150     vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
    151     ! vol_cld(1) = 4./3. * pi * rmin_cld*rmin_cld*rmin_cld
    152 
    153     do i=1,nbin_cld-1
    154         rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
    155         vol_cld(i+1) = vol_cld(i) * vrat_cld
    156     enddo
    157 
    158     do i=1,nbin_cld
    159         rb_cld(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * rad_cld(i)
    160         dr_cld(i)  = rb_cld(i+1) - rb_cld(i)
    161     enddo
    162     rb_cld(nbin_cld+1) = rbmax_cld
    163     dr_cld(nbin_cld)   = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
    164 
    165     print*, ' '
    166     print*,'Microphysics: size bin information:'
    167     print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
    168     print*,'-----------------------------------'
    169     do i=1,nbin_cld
    170         write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i), dr_cld(i)
    171     enddo
    172     write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
    173     print*,'-----------------------------------'
    174 
    175     ! we save that so that it is not computed at each timestep and gridpoint
    176     rb_cld(:) = log(rb_cld(:))
    177 
    178     ! Contact parameter of water ice on dust ( m=cos(theta) )
    179     ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    180     ! mteta is initialized in conf_phys
    181     write(*,*) 'water_param contact parameter:', mteta
    182 
    183     ! Volume of a water molecule (m3)
    184     vo1 = mh2o / dble(rho_ice)
    185     ! Variance of the ice and CCN distributions
    186     sigma_ice = sqrt(log(1.+nuice_sed))
    187 
    188     write(*,*) 'Variance of ice & CCN distribs:', sigma_ice
    189     write(*,*) 'nuice for sedimentation       :', nuice_sed
    190     write(*,*) 'Volume of a water molecule    :', vo1
    191 
    192     test_flag = .false.
    193     firstcall = .false.
     146  ! Volume ratio between two adjacent bins
     147  ! vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
     148  ! vrat_cld = exp(vrat_cld)
     149  vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
     150  vrat_cld = exp(vrat_cld)
     151  write(*,*) "vrat_cld", vrat_cld
     152
     153  rb_cld(1)  = rbmin_cld
     154  rad_cld(1) = rmin_cld
     155  vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
     156  ! vol_cld(1) = 4./3. * pi * rmin_cld*rmin_cld*rmin_cld
     157
     158  do i=1,nbin_cld-1
     159    rad_cld(i+1)  = rad_cld(i) * vrat_cld**(1./3.)
     160    vol_cld(i+1)  = vol_cld(i) * vrat_cld
     161  enddo
     162     
     163  do i=1,nbin_cld
     164    rb_cld(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * rad_cld(i)
     165    dr_cld(i)  = rb_cld(i+1) - rb_cld(i)
     166  enddo
     167  rb_cld(nbin_cld+1) = rbmax_cld
     168  dr_cld(nbin_cld)   = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
     169
     170  print*, ' '
     171  print*,'Microphysics: size bin information:'
     172  print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
     173  print*,'-----------------------------------'
     174  do i=1,nbin_cld
     175    write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i), dr_cld(i)
     176  enddo
     177  write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
     178  print*,'-----------------------------------'
     179
     180  ! we save that so that it is not computed at each timestep and gridpoint
     181  ! rb_cld(i) = log(rb_cld(i)) 
     182  rb_cld(:) = log(rb_cld(:))
     183
     184  ! Contact parameter of water ice on dust ( m=cos(theta) )
     185  !  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     186  !  mteta is initialized in conf_phys
     187  write(*,*) 'water_param contact parameter:', mteta
     188
     189  ! Volume of a water molecule (m3)
     190  vo1 = mh2o / dble(rho_ice)
     191  ! Variance of the ice and CCN distributions
     192  sigma_ice = sqrt(log(1.+nuice_sed))
     193     
     194  write(*,*) 'Variance of ice & CCN distribs :', sigma_ice
     195  write(*,*) 'nuice for sedimentation:', nuice_sed
     196  write(*,*) 'Volume of a water molecule:', vo1
     197
     198  test_flag = .false.
     199  firstcall=.false.
    194200ENDIF
    195201
    … …  
    197203! 1. Initialisation
    198204!=============================================================
     205
    199206res_out(:,:) = 0
    200207rice(:,:) = 1.e-8
    201208
    202 ! Initialize the temperature, tracers
    203 zt(:,:) = pt(:,:)
    204 zq(:,:,:) = pq(:,:,:)
    205 subpdtcloud(:,:) = 0
    206 
    207 WHERE ( zq(:,:,:) < 1.e-30 ) zq(:,:,:) = 1.e-30
     209!     Initialize the temperature, tracers
     210zt(:,:)=pt(:,:)
     211zq(:,:,:)=pq(:,:,:)
     212subpdtcloud(:,:)=0
     213
     214WHERE( zq(:,:,:) < 1.e-30 ) zq(:,:,:) = 1.e-30
    208215
    209216!=============================================================
    210217! 2. Compute saturation
    211218!=============================================================
     219
    212220dev2 = 1. / ( sqrt(2.) * sigma_ice )
    213221
    214222! Compute the condensable amount of water vapor or the sublimable water ice (if negative)
    215 call watersat(ngrid*nlay,zt+pdt*ptimestep,pplay,zqsat)
    216 zpotcond_full = (zq(:,:,igcm_h2o_vap)+pdq(:,:,igcm_h2o_vap)*ptimestep) - zqsat
    217 zpotcond_full = max(zpotcond_full,-zq(:,:,igcm_h2o_ice))
     223! Attention ici pdt*ptimestep peut etre severe
     224call watersat(ngrid*nlay,max(1.,zt+pdt*ptimestep),pplay,zqsat) ! DIFF: "temp+trend" at least 1
     225zpotcond_full=(zq(:,:,igcm_h2o_vap)+pdq(:,:,igcm_h2o_vap)*ptimestep) - zqsat
     226zpotcond_full=max(zpotcond_full,-zq(:,:,igcm_h2o_ice))
    218227call watersat(ngrid*nlay,zt,pplay,zqsat)
    219 zpotcond_inst = zq(:,:,igcm_h2o_vap) - zqsat
    220 zpotcond_inst = max(zpotcond_inst,-zq(:,:,igcm_h2o_ice))
    221 ! zpotcond is the most strict criterion between the two
     228zpotcond_inst=zq(:,:,igcm_h2o_vap) - zqsat
     229zpotcond_inst=max(zpotcond_inst,-zq(:,:,igcm_h2o_ice))
     230! zpotcond is the most strict criterion between the two 
    222231DO l=1,nlay
    223     DO ig=1,ngrid
    224         if (zpotcond_full(ig,l) > 0.) then
    225             zpotcond(ig,l) = max(zpotcond_inst(ig,l),zpotcond_full(ig,l))
    226         else if (zpotcond_full(ig,l) <= 0.) then
    227             zpotcond(ig,l) = min(zpotcond_inst(ig,l),zpotcond_full(ig,l))
    228         endif ! (zpotcond_full.gt.0.) then
    229     ENDDO !ig=1,ngrid
     232  DO ig=1,ngrid
     233    if (zpotcond_full(ig,l).gt.0.) then
     234      zpotcond(ig,l)=max(zpotcond_inst(ig,l),zpotcond_full(ig,l))
     235    else if (zpotcond_full(ig,l).le.0.) then
     236      zpotcond(ig,l)=min(zpotcond_inst(ig,l),zpotcond_full(ig,l))
     237    endif ! (zpotcond_full.gt.0.) then
     238    ! zpotcond(ig,l)=1.e-2
     239  ENDDO !ig=1,ngrid
    230240ENDDO !l=1,nlay
    231 
     241     
    232242countcells = 0
    233 zimicro(1:ngrid,1:nlay)=imicro
     243zimicro(:,:)=imicro
    234244count_micro(:,:)=0
    235245
    236246! Main loop over the GCM's grid
    237247DO l=1,nlay
    238     DO ig=1,ngrid
    239         ! Subtimestep: here we go
    240         IF (cloud_adapt_ts) call adapt_imicro(ptimestep,zpotcond(ig,l),zimicro(ig,l))
    241         spenttime = 0.
    242         ending_ts=.false.
    243         DO while (.not.ending_ts)
    244             microtimestep=ptimestep/real(zimicro(ig,l))
    245             ! Initialize tracers for scavenging + hdo computations (JN)
    246             zq0(ig,l,:) = zq(ig,l,:)
    247 
    248             ! Check if we are integrating over ptimestep
    249             if (spenttime+microtimestep >= ptimestep) then
    250                 microtimestep = ptimestep-spenttime
    251                 ! If so : last call !
    252                 ending_ts = .true.
    253             endif ! (spenttime+microtimestep.ge.ptimestep) then
    254             call watersat(1,(/zt(ig,l)/),(/pplay(ig,l)/),zqsat_tmp)
    255             zqsat(ig,l) = zqsat_tmp(1)
    256             ! Get the partial pressure of water vapor and its saturation ratio
    257             ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
    258             satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
     248  DO ig=1,ngrid
     249    ! Subtimestep : here we go
     250    IF (cloud_adapt_ts) call adapt_imicro(ptimestep,zpotcond(ig,l), zimicro(ig,l))
     251    spenttime = 0.
     252    dMicetot= 0. ! DIFF: dMicetot=new var
     253    ending_ts=.false.
     254    DO while (.not.ending_ts)
     255      call watersat(1,(/zt(ig,l)/),(/pplay(ig,l)/),zqsat_tmp)
     256      zqsat(ig,l)=zqsat_tmp(1)
     257      ! Get the partial pressure of water vapor and its saturation ratio
     258      ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
     259      satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
     260      microtimestep=ptimestep/real(zimicro(ig,l))
     261      ! if (satu.lt.1.0) then
     262      !   microtimestep=ptimestep/30.
     263      !   write(*,*) "saturation correction" !JN
     264      ! endif
     265     
     266      ! Initialize tracers for scavenging + hdo computations (JN)
     267      zq0(ig,l,:)=zq(ig,l,:)
     268
     269      ! Check if we are integrating over ptimestep
     270      if (spenttime+microtimestep.ge.ptimestep) then
     271        microtimestep=ptimestep-spenttime
     272        ! If so : last call !
     273        ending_ts=.true.
     274      endif! (spenttime+microtimestep.ge.ptimestep) then
    259275
    260276!=============================================================
    261277! 3. Nucleation
    262278!=============================================================
    263             IF ( satu >= 1. ) THEN ! if there is condensation
    264                 call updaterccn(zq(ig,l,igcm_dust_mass),zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
    265 
    266                 ! Expand the dust moments into a binned distribution
    267                 Mo = zq(ig,l,igcm_dust_mass)  *tauscaling(ig) + 1.e-30
    268                 No = zq(ig,l,igcm_dust_number)*tauscaling(ig) + 1.e-30
    269                 Rn = rdust(ig,l)
    270                 Rn = -log(Rn)
    271                 Rm = Rn - 3. * sigma_ice*sigma_ice
    272                 n_derf = derf( (rb_cld(1)+Rn) *dev2)
    273                 m_derf = derf( (rb_cld(1)+Rm) *dev2)
    274                 do i = 1, nbin_cld
    275                     n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
    276                     m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
    277                     n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
    278                     m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
    279                     n_aer(i) = n_aer(i) + 0.5 * No * n_derf
    280                     m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
    281                 enddo
    282 
    283 !                sumcheck = 0
    284 !                do i = 1, nbin_cld
    285 !                    sumcheck = sumcheck + n_aer(i)
    286 !                enddo
    287 !                sumcheck = abs(sumcheck/No - 1)
    288 !                if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
    289 !                    print*, "WARNING, No sumcheck PROBLEM"
    290 !                    print*, "sumcheck, No",sumcheck, No
    291 !                    print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
    292 !                    print*, "Dust binned distribution", n_aer
    293 !                endif
    294 !
    295 !                sumcheck = 0
    296 !                do i = 1, nbin_cld
    297 !                    sumcheck = sumcheck + m_aer(i)
    298 !                enddo
    299 !                sumcheck = abs(sumcheck/Mo - 1)
    300 !                if ((sumcheck .gt. 1e-5) .and.  (1./Rn .gt. rmin_cld)) then
    301 !                    print*, "WARNING, Mo sumcheck PROBLEM"
    302 !                    print*, "sumcheck, Mo",sumcheck, Mo
    303 !                    print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
    304 !                    print*, "Dust binned distribution", m_aer
    305 !                endif
    306 
    307                 ! Get the rates of nucleation
    308                 call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
    309 
    310                 dN = 0.
    311                 dM = 0.
    312                 do i = 1, nbin_cld
    313                     dN = dN + n_aer(i)*(exp(-rate(i)*microtimestep)-1.)
    314                     dM = dM + m_aer(i)*(exp(-rate(i)*microtimestep)-1.)
    315                 enddo
    316 
    317                 ! Update Dust particles
    318                 zq(ig,l,igcm_dust_mass)   = zq(ig,l,igcm_dust_mass)   + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    319                 zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    320                 ! Update CCNs
    321                 zq(ig,l,igcm_ccn_mass)   = zq(ig,l,igcm_ccn_mass)   - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    322                 zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    323 
    324             ENDIF ! of is satu >1
     279
     280      IF ( satu .ge. 1. ) THEN ! if there is condensation
     281        call updaterccn(zq(ig,l,igcm_dust_mass),zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
     282
     283        ! Expand the dust moments into a binned distribution
     284        Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)   + 1.e-30
     285        No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
     286        Rn = rdust(ig,l)
     287        Rn = -log(Rn)
     288        Rm = Rn - 3. * sigma_ice*sigma_ice 
     289        n_derf = derf( (rb_cld(1)+Rn) *dev2)
     290        m_derf = derf( (rb_cld(1)+Rm) *dev2)
     291        do i = 1, nbin_cld
     292          n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
     293          m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
     294          n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
     295          m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
     296          n_aer(i) = n_aer(i) + 0.5 * No * n_derf
     297          m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
     298        enddo
     299     
     300        ! sumcheck = 0
     301        ! do i = 1, nbin_cld
     302        !   sumcheck = sumcheck + n_aer(i)
     303        ! enddo
     304        ! sumcheck = abs(sumcheck/No - 1)
     305        ! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
     306        !   print*, "WARNING, No sumcheck PROBLEM"
     307        !   print*, "sumcheck, No",sumcheck, No
     308        !   print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
     309        !   print*, "Dust binned distribution", n_aer
     310        ! endif
     311        !       
     312        ! sumcheck = 0
     313        ! do i = 1, nbin_cld
     314        !   sumcheck = sumcheck + m_aer(i)
     315        ! enddo
     316        ! sumcheck = abs(sumcheck/Mo - 1)
     317        ! if ((sumcheck .gt. 1e-5) .and.  (1./Rn .gt. rmin_cld)) then
     318        !   print*, "WARNING, Mo sumcheck PROBLEM"
     319        !   print*, "sumcheck, Mo",sumcheck, Mo
     320        !   print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
     321        !   print*, "Dust binned distribution", m_aer
     322        ! endif
     323
     324        ! Get the rates of nucleation
     325        call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
     326     
     327        dN = 0.
     328        dM = 0.
     329        do i = 1, nbin_cld
     330          dN = dN + n_aer(i)*(exp(-rate(i)*microtimestep)-1.)
     331          dM = dM + m_aer(i)*(exp(-rate(i)*microtimestep)-1.)
     332        enddo
     333
     334        ! Update Dust particles
     335        zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
     336        zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
     337        ! Update CCNs
     338        zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
     339        zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
     340     
     341      ENDIF ! of is satu >1
    325342
    326343!=============================================================
    … …  
    330347! Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
    331348! to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
    332             IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig) >= 1.) THEN ! we trigger crystal growth
    333 
    334                 call updaterice_micro(zq(ig,l,igcm_h2o_ice),zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),tauscaling(ig),rice(ig,l),rhocloud(ig,l))
    335 
    336                 No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
    337 
    338                 ! saturation at equilibrium
    339                 ! rice should not be too small, otherwise seq value is not valid
    340                 seq  = exp(2.*sig_h2o(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)*max(rice(ig,l),1.e-7)))
    341 
    342                 ! get resistance growth
    343                 call growthrate(zt(ig,l),pplay(ig,l),real(ph2o/satu),rice(ig,l),res,Dv)
    344 
    345                 res_out(ig,l) = res
    346 
    347                 ! implicit scheme of mass growth
    348                 ! cste here must be computed at each step
    349                 cste = 4*pi*rho_ice*microtimestep
    350 
    351                 dMice = (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))/(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.)
    352 
    353 
    354                 ! With the above scheme, dMice cannot be bigger than vapor, but can be bigger than all available ice.
    355                 dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
    356                 dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
    357 
    358                 zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) + dMice
    359                 zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) - dMice
    360 
    361                 countcells = countcells + 1
    362 
    363                 ! latent heat release
    364                 lw = (2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3
    365                 subpdtcloud(ig,l) = dMice*lw/cpp/microtimestep
    366 
    367                 ! Special case of the isotope of water HDO
    368                 if (hdo) then
    369                     ! condensation
    370                     if (dMice > 0.) then
    371                         ! do we use fractionation?
    372                         if (hdofrac) then
    373                             ! Calculation of the HDO vapor coefficient
    374                             Dv_hdo = 1./3. * sqrt( 8*kbz*zt(ig,l)/(pi*mhdo/nav) )* kbz * zt(ig,l) / (pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo)* sqrt(1.+mhdo/mco2) )
    375                             ! Calculation of the fractionnation coefficient at equilibrium
    376                             alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
    377                             ! alpha = exp(13525./zt(ig,l)**2.-5.59d-2)  !Lamb
    378                             ! Calculation of the 'real' fractionnation coefficient
    379                             alpha_c(ig,l) = (alpha(ig,l)*satu) / ((alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.)
    380                             ! alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics
    381                         else
    382                             alpha_c(ig,l) = 1.d0
    383                         endif
    384                         if (zq0(ig,l,igcm_h2o_vap) > qperemin) then
    385                             dMice_hdo = dMice*alpha_c(ig,l)*( zq0(ig,l,igcm_hdo_vap)/zq0(ig,l,igcm_h2o_vap) )
    386                         else
    387                             dMice_hdo = 0.
    388                         endif
    389                     ! sublimation
    390                     else
    391                         if (zq0(ig,l,igcm_h2o_ice) > qperemin) then
    392                             dMice_hdo = dMice*( zq0(ig,l,igcm_hdo_ice)/zq0(ig,l,igcm_h2o_ice) )
    393                         else
    394                             dMice_hdo = 0.
    395                         endif
    396                     endif ! if (dMice > 0.)
    397 
    398                     dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice))
    399                     dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap))
    400 
    401                     zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo
    402                     zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo
    403 
    404                 endif ! if (hdo)
     349      IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
     350        call updaterice_micro(zq(ig,l,igcm_h2o_ice),zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),tauscaling(ig),rice(ig,l),rhocloud(ig,l))
     351
     352        No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
     353
     354        ! saturation at equilibrium
     355        ! rice should not be too small, otherwise seq value is not valid
     356        seq  = exp(2.*sig_h2o(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)*max(rice(ig,l),1.e-7)))
     357     
     358        ! get resistance growth
     359        call growthrate(zt(ig,l),pplay(ig,l),real(ph2o/satu),rice(ig,l),res,Dv)
     360
     361        res_out(ig,l) = res
     362
     363        ! implicit scheme of mass growth
     364        ! cste here must be computed at each step
     365        cste = 4*pi*rho_ice*microtimestep
     366
     367        dMice = (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))/(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.)
     368
     369        ! With the above scheme, dMice cannot be bigger than vapor, but can be bigger than all available ice.
     370        ! dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
     371        dMice = max(dMice,-zq(ig,l,igcm_h2o_ice) - pdq(ig,l,igcm_h2o_ice)*microtimestep) ! DIFF: take trend into account
     372        ! dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
     373        dMice = min(dMice,zq(ig,l,igcm_h2o_vap) + pdq(ig,l,igcm_h2o_vap)*microtimestep)
     374
     375        zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
     376        zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
     377
     378        countcells = countcells + 1
     379     
     380        ! latent heat release
     381        lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3
     382        subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
     383     
     384        ! DIFF: trend is enforce in a range, stabilize the scheme ?
     385        if (subpdtcloud(ig,l)*microtimestep.gt.5.0) then
     386          subpdtcloud(ig,l)=5./microtimestep
     387        endif! (subpdtcloud(ig,l)*microtimestep.gt.5.0) then
     388        if (subpdtcloud(ig,l)*microtimestep.lt.-5.0) then
     389            subpdtcloud(ig,l)=-5./microtimestep
     390        endif! (subpdtcloud(ig,l)*microtimestep.gt.5.0) then
     391
     392        ! Special case of the isotope of water HDO   
     393        if (hdo) then
     394          ! condensation
     395          if (dMice.gt.0.0) then
     396            ! do we use fractionation?
     397            if (hdofrac) then
     398              ! Calculation of the HDO vapor coefficient
     399              Dv_hdo = 1./3. * sqrt( 8*kbz*zt(ig,l)/(pi*mhdo/nav) ) * kbz * zt(ig,l) / ( pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo) * sqrt(1.+mhdo/mco2) )
     400              ! Calculation of the fractionnation coefficient at equilibrium
     401              alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
     402              ! alpha = exp(13525./zt(ig,l)**2.-5.59d-2)  !Lamb
     403              ! Calculation of the 'real' fractionnation coefficient
     404              alpha_c(ig,l) = (alpha(ig,l)*satu)/( (alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.)
     405              ! alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics
     406            else
     407              alpha_c(ig,l) = 1.d0
     408            endif
     409            if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then
     410              dMice_hdo = dMice*alpha_c(ig,l)*( zq0(ig,l,igcm_hdo_vap)/zq0(ig,l,igcm_h2o_vap) )
     411            else
     412              dMice_hdo=0.
     413            endif
     414          !! sublimation
     415          else
     416            if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then
     417              dMice_hdo=dMice*( zq0(ig,l,igcm_hdo_ice)/zq0(ig,l,igcm_h2o_ice) )
     418            else
     419              dMice_hdo=0.
     420            endif
     421          endif !if (dMice.gt.0.0)
     422
     423          dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice))
     424          dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap))
     425
     426          zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo
     427          zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo
     428
     429        endif ! if (hdo)       
    405430
    406431!=============================================================
    407432! 5. Dust cores released, tendancies, latent heat, etc ...
    408433!=============================================================
    409                 ! If all the ice particles sublimate, all the condensation
    410                 ! nuclei are released:
    411                 if (zq(ig,l,igcm_h2o_ice) <= 1.e-28) then
    412                     ! Water
    413                     zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) + zq(ig,l,igcm_h2o_ice)
    414                     zq(ig,l,igcm_h2o_ice) = 0.
    415                     if (hdo) then
    416                         zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) + zq(ig,l,igcm_hdo_ice)
    417                         zq(ig,l,igcm_hdo_ice) = 0.
    418                     endif
    419                     ! Dust particles
    420                     zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccn_mass)
    421                     zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccn_number)
    422                     ! CCNs
    423                     zq(ig,l,igcm_ccn_mass) = 0.
    424                     zq(ig,l,igcm_ccn_number) = 0.
    425                 endif
    426 
    427             ELSE
    428                 ! Initialization of dMice when it's not computed
    429                 dMice = 0 ! no condensation/sublimation to account for
    430             ENDIF ! of if Nccn>1
    431 
    432             ! No more getting tendency: we increment tracers & temp directly
    433 
    434             ! Increment tracers & temp for following microtimestep
    435             ! Dust:
    436             ! Special treatment of dust_mass / number for scavenging?
    437             IF (scavenging) THEN
    438                 zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + pdq(ig,l,igcm_dust_mass)*microtimestep
    439                 zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + pdq(ig,l,igcm_dust_number)*microtimestep
    440             ELSE
    441                 zq(ig,l,igcm_dust_mass) = zq0(ig,l,igcm_dust_mass)
    442                 zq(ig,l,igcm_dust_number) = zq0(ig,l,igcm_dust_number)
    443             ENDIF !(scavenging) THEN
    444             zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + pdq(ig,l,igcm_ccn_mass)*microtimestep
    445             zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + pdq(ig,l,igcm_ccn_number)*microtimestep
    446 
    447             ! Water:
    448             zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) + pdq(ig,l,igcm_h2o_ice)*microtimestep
    449             zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) + pdq(ig,l,igcm_h2o_vap)*microtimestep
    450 
    451             ! HDO (if computed):
    452             IF (hdo) THEN
    453                 zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice) + pdq(ig,l,igcm_hdo_ice)*microtimestep
    454                 zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) + pdq(ig,l,igcm_hdo_vap)*microtimestep
    455             ENDIF ! hdo
    456 
    457             ! Could also set subpdtcloud to 0 if not activice to make it simpler or change name of the flag
    458             IF (.not.activice) subpdtcloud(ig,l) = 0.
    459 
    460             ! Temperature:
    461             zt(ig,l) = zt(ig,l) + (pdt(ig,l)+subpdtcloud(ig,l))*microtimestep
    462 
    463             ! Prevent negative tracers ! JN
    464             WHERE (zq(ig,l,:) < 1.e-30) zq(ig,l,:) = 1.e-30
    465 
    466             IF (cloud_adapt_ts) then
    467                 ! Estimation of how much is actually condensing/subliming
    468                 IF (spenttime /= 0) then
    469                     zdq = (dMice/spenttime)*(ptimestep-spenttime)
    470                 ELSE
    471                     ! Initialization for spenttime=0
    472                     zdq = zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep)
    473                 ENDIF
    474                     ! Threshold with powerlaw (sanity check)
    475                     zdq = min(abs(zdq),abs(zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep)))
    476                     call adapt_imicro(ptimestep,zdq,zimicro(ig,l))
    477             ENDIF ! (cloud_adapt_ts) then
    478             ! Increment time spent in here
    479             spenttime = spenttime + microtimestep
    480             count_micro(ig,l) = count_micro(ig,l) + 1
    481         ENDDO ! while (.not. ending_ts)
    482     ENDDO ! of ig loop
     434
     435        ! If all the ice particles sublimate, all the condensation
     436        ! nuclei are released:
     437        if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
     438          ! Water
     439          zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) + zq(ig,l,igcm_h2o_ice)
     440          zq(ig,l,igcm_h2o_ice) = 0.
     441          if (hdo) then
     442            zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) + zq(ig,l,igcm_hdo_ice)
     443            zq(ig,l,igcm_hdo_ice) = 0.
     444          endif
     445          ! Dust particles
     446          zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccn_mass)
     447          zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccn_number)
     448          ! CCNs
     449          zq(ig,l,igcm_ccn_mass) = 0.
     450          zq(ig,l,igcm_ccn_number) = 0.
     451        endif
     452     
     453      ELSE
     454        ! Initialization of dMice when it's not computed
     455        dMice=0 ! no condensation/sublimation to account for
     456        subpdtcloud(ig,l)=0 ! no condensation/sublimation to account for
     457      ENDIF !of if Nccn>1
     458     
     459      ! No more getting tendency : we increment tracers & temp directly
     460
     461      ! Increment tracers & temp for following microtimestep
     462      ! Dust :
     463      ! Special treatment of dust_mass / number for scavenging ?
     464      IF (scavenging) THEN
     465        zq(ig,l,igcm_dust_mass) =zq(ig,l,igcm_dust_mass)+pdq(ig,l,igcm_dust_mass)*microtimestep
     466        zq(ig,l,igcm_dust_number) =zq(ig,l,igcm_dust_number)+pdq(ig,l,igcm_dust_number)*microtimestep
     467      ELSE
     468        zq(ig,l,igcm_dust_mass) =zq0(ig,l,igcm_dust_mass)
     469        zq(ig,l,igcm_dust_number) =zq0(ig,l,igcm_dust_number)
     470      ENDIF !(scavenging) THEN
     471      zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + pdq(ig,l,igcm_ccn_mass)*microtimestep
     472      zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + pdq(ig,l,igcm_ccn_number)*microtimestep
     473
     474      ! Water :
     475      zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*microtimestep
     476      zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*microtimestep
     477
     478      ! HDO (if computed) :
     479      IF (hdo) THEN
     480        zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+pdq(ig,l,igcm_hdo_ice)*microtimestep
     481        zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)+pdq(ig,l,igcm_hdo_vap)*microtimestep
     482      ENDIF ! hdo
     483
     484      ! Could also set subpdtcloud to 0 if not activice to make it simpler or change name of the flag
     485      IF (.not.activice) subpdtcloud(ig,l)=0.
     486
     487      ! Temperature :
     488      zt(ig,l) = zt(ig,l)+(pdt(ig,l)+subpdtcloud(ig,l))*microtimestep
     489
     490      ! Prevent negative tracers ! JN
     491      WHERE (zq(ig,l,:) < 1.e-30) zq(ig,l,:) = 1.e-30
     492
     493      IF (cloud_adapt_ts) then
     494        ! Estimation of how much is actually condensing/subliming
     495         dMicetot=dMicetot+abs(dMice) ! DIFF: accumulation of absolute dMice
     496        ! dMicetot=dMicetot+dMice
     497        ! write(*,*) "dMicetot= ", dMicetot
     498        ! write(*,*) "we are in (l,ig):", l,ig !JN
     499        IF (spenttime.ne.0) then
     500          zdq=(dMicetot/spenttime)!*(ptimestep-spenttime)
     501        ELSE
     502          !Initialization for spenttime=0
     503          zdq=zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep)
     504        ENDIF
     505        ! Threshold with powerlaw (sanity check)
     506        ! zdq=min(abs(zdq),
     507        ! & abs(zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep)))
     508        zdq=abs(zdq)
     509        call adapt_imicro(ptimestep,zdq,zimicro(ig,l))
     510      ENDIF! (cloud_adapt_ts) then
     511      ! Increment time spent in here
     512      spenttime=spenttime+microtimestep
     513      count_micro(ig,l)=count_micro(ig,l)+1
     514    ENDDO ! while (.not. ending_ts)
     515  ENDDO ! of ig loop
    483516ENDDO ! of nlayer loop
    484517
    485518!------ Useful outputs to check how it went
    486 call write_output("zpotcond_inst","zpotcond_inst "//"microphysics","(kg/kg)",zpotcond_inst(:,:))
    487 call write_output("zpotcond_full","zpotcond_full "//"microphysics","(kg/kg)",zpotcond_full(:,:))
    488 call write_output("zpotcond","zpotcond "//"microphysics","(kg/kg)",zpotcond(:,:))
    489 call write_output("count_micro","count_micro "//"after microphysics","integer",count_micro(:,:))
    490 
    491 !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
     519call write_output("zpotcond_inst","zpotcond_inst microphysics","(kg/kg)",zpotcond_inst(:,:))
     520call write_output("zpotcond_full","zpotcond_full microphysics","(kg/kg)",zpotcond_full(:,:))
     521call write_output("zpotcond","zpotcond microphysics","(kg/kg)",zpotcond(:,:))
     522call write_output("count_micro","count_micro after microphysics","integer",count_micro(:,:))
     523
     524!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
    492525IF (test_flag) then
    493 !    error2d(:) = 0.
    494     DO l=1,nlay
    495         DO ig=1,ngrid
    496 !            error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
    497             satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
    498             satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
    499         ENDDO
     526  ! error2d(:) = 0.
     527  DO l=1,nlay
     528    DO ig=1,ngrid
     529      ! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
     530      satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
     531      satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
    500532    ENDDO
    501     write(*,*) 'count is ',countcells, ' i.e. ',countcells*100/(nlay*ngrid), '% for microphys computation'
     533  ENDDO
     534
     535  print*, 'count is ',countcells, ' i.e. ', countcells*100/(nlay*ngrid), '% for microphys computation'
    502536ENDIF ! endif test_flag
    503537
    504538END SUBROUTINE improvedclouds
    505539
    506 !=======================================================================
    507 
    508 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    509 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     540!=============================================================
    510541! The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
    511 ! It is an analytical solution to the ice radius growth equation,
    512 ! with the approximation of a constant 'reduced' cunningham correction factor
    513 ! (lambda in growthrate.F) taken at radius req instead of rice
    514 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    515 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     542! It is an analytical solution to the ice radius growth equation,
     543! with the approximation of a constant 'reduced' cunningham correction factor
     544! (lambda in growthrate.F) taken at radius req instead of rice   
     545!=============================================================
     546
     547! subroutine phi(rice,req,coeff1,coeff2,time)
     548!     
     549! implicit none
     550!     
     551! ! inputs
     552! real rice ! ice radius
     553! real req  ! ice radius at equilibirum
     554! real coeff1  ! coeff for the log
     555! real coeff2  ! coeff for the arctan
    516556!
    517 ! subroutine phi(rice,req,coeff1,coeff2,time)
    518 !
    519 ! implicit none
    520 !
    521 ! ! inputs
    522 ! real, intent(in) :: rice ! ice radius
    523 ! real, intent(in) :: req  ! ice radius at equilibirum
    524 ! real, intent(in) :: coeff1  ! coeff for the log
    525 ! real, intent(in) :: coeff2  ! coeff for the arctan
    526 !
    527 ! ! output
    528 ! real, intent(out) :: time
    529 !
     557! ! output     
     558! real time
     559!     
    530560! !local
    531 ! real :: var
    532 !
    533 ! ! 1.73205 is sqrt(3)
    534 !
    535 ! var = max(abs(rice-req) / sqrt(rice*rice + rice*req  + req*req),1e-30)
    536 !
    537 ! time = coeff1 * log( var ) + coeff2 * 1.73205 * atan( (2*rice+req) / (1.73205*req) )
    538 !
    539 ! end subroutine phi
     561! real var
     562!     
     563!  1.73205 is sqrt(3)
     564!     
     565! var = max(
     566! &  abs(rice-req) / sqrt(rice*rice + rice*req  + req*req),1e-30)
     567!           
     568! time =
     569! &   coeff1 *
     570! &   log( var )
     571! & + coeff2 * 1.73205 *
     572! &   atan( (2*rice+req) / (1.73205*req) )
     573!     
     574! return
     575! end
    540576
    541577!=======================================================================
    … …  
    550586! efficiency.
    551587
    552 real, intent(in) :: ptimestep ! total duration of physics (sec)
    553 real, intent(in) :: potcond ! condensible vapor / sublimable ice(kg/kg)
    554 integer, intent(out) :: zimicro ! number of ptimestep division
     588real,intent(in) :: ptimestep ! total duration of physics (sec)
     589real,intent(in) :: potcond ! condensible vapor / sublimable ice(kg/kg)
    555590real :: alpha, beta ! Coefficients for power law
    556591real :: defstep,coef ! Default ptimestep of 7.5 mins (iphysiq=5)
    557 
    558 ! Default ptimestep: defstep (7.5 mins)
    559 defstep = 88775.*5./960.
    560 coef = ptimestep/defstep
     592integer,intent(out) :: zimicro ! number of ptimestep division
     593
     594! Default ptimestep : defstep (7.5 mins)
     595defstep=88775.*5./960.
     596coef=ptimestep/defstep
    561597! Conservative coefficients :
    562 alpha = 1.81846504e+11
    563 beta = 1.54550140e+00
    564 zimicro = ceiling(coef*min(max(alpha*abs(potcond)**beta,5.),7000.))
     598! alpha=1.81846504e+11
     599! beta=1.54550140e+00
     600alpha=1.88282793e+05 ! DIFF: new power law coefficients
     601beta=4.57764370e-01
     602zimicro=ceiling(coef*min(max(alpha*abs(potcond)**beta,50.),7000.))
     603zimicro=2*zimicro ! DIFF: prefiction times two, allow to complete year
    565604
    566605END SUBROUTINE adapt_imicro
Note: See TracChangeset for help on using the changeset viewer.