Changeset 3919 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Sep 26, 2025, 12:30:38 PM (7 months ago)
Author:
jbclement
Message:

Mars PCM:
Formating "improveclouds" into F90 and simplifying some instructions.
JBC

File:
1 moved

Legend:

Unmodified
Added
Removed

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

    r3918 r3919  
    1       MODULE improvedclouds_mod
    2 
    3       IMPLICIT NONE
    4 
    5       CONTAINS
    6 
    7       subroutine improvedclouds(ngrid,nlay,ptimestep,
    8      &             pplay,pt,pdt,pq,pdq,nq,tauscaling,
    9      &             imicro,zt,zq)
    10       USE updaterad, ONLY: updaterice_micro, updaterccn
    11       USE watersat_mod, ONLY: watersat
    12       use tracer_mod, only: rho_ice, nuice_sed, igcm_h2o_vap,
    13      &                      igcm_h2o_ice, igcm_dust_mass,
    14      &                      igcm_dust_number, igcm_ccn_mass,
    15      &                      igcm_ccn_number,
    16      &                      igcm_hdo_vap,igcm_hdo_ice,
    17      &                      qperemin
    18       use conc_mod, only: mmean
    19       use comcstfi_h, only: pi, cpp
    20       use microphys_h, only: nbin_cld, rad_cld, mteta, kbz, nav, rgp
    21       use microphys_h, only: mco2, vo1, mh2o, mhdo, molco2, molhdo, To
    22       use nuclea_mod, only: nuclea
    23       use sig_h2o_mod, only: sig_h2o
    24       use growthrate_mod, only: growthrate
    25       use write_output_mod, only: write_output
    26       use callkeys_mod, only: activice, scavenging, cloud_adapt_ts,
    27      &                        hdo, hdofrac
    28       implicit none
    29      
    30      
    31 c------------------------------------------------------------------
    32 c  This routine is used to form clouds when a parcel of the GCM is
    33 c    saturated. It includes the ability to have supersaturation, a
    34 c    computation of the nucleation rates, growthrates and the
    35 c    scavenging of dust particles by clouds.
    36 c  It is worth noting that the amount of dust is computed using the
    37 c    dust optical depth computed in aeropacity.F. That's why
    38 c    the variable called "tauscaling" is used to convert
    39 c    pq(dust_mass) and pq(dust_number), which are relative
    40 c    quantities, to absolute and realistic quantities stored in zq.
    41 c    This has to be done to convert the inputs into absolute
    42 c    values, but also to convert the outputs back into relative
    43 c    values which are then used by the sedimentation and advection
    44 c    schemes.
    45 
    46 c  Authors: J.-B. Madeleine, based on the work by Franck Montmessin
    47 c           (October 2011)
    48 c           T. Navarro, debug,correction, new scheme (October-April 2011)
    49 c           A. Spiga, optimization (February 2012)
    50 c           J. Naar, adaptative subtimestep now done here (June 2023)
    51 c------------------------------------------------------------------
    52 c     Inputs/outputs:
    53 
    54       INTEGER, INTENT(IN) :: ngrid,nlay
    55       INTEGER, INTENT(IN) :: nq                 ! nombre de traceurs
    56       REAL, INTENT(IN) :: ptimestep             ! pas de temps physique (s)
    57       REAL, INTENT(IN) :: pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)           
    58       REAL, INTENT(IN) :: pt(ngrid,nlay)        ! temperature at the middle of the
    59                                                 ! layers (K)
    60       REAL, INTENT(IN) :: pdt(ngrid,nlay)       ! temperature tendency (K/s)
    61       REAL, INTENT(IN) :: pq(ngrid,nlay,nq)     ! tracer (kg/kg)
    62       REAL, INTENT(IN) :: pdq(ngrid,nlay,nq)    ! tracer tendency (kg/kg/s)
    63       REAL, INTENT(IN) :: tauscaling(ngrid)     ! Convertion factor for qdust and Ndust
    64       INTEGER, INTENT(IN) :: imicro             ! nb. microphy calls(retrocompatibility)
    65      
    66       REAL, INTENT(OUT) :: zq(ngrid,nlay,nq)  ! tracers post microphy (kg/kg)
    67       REAL, INTENT(OUT) :: zt(ngrid,nlay)     ! temperature post microphy (K)
    68 
    69 c------------------------------------------------------------------
    70 c     Local variables:
    71 
    72       LOGICAL, SAVE ::  firstcall = .true.
     1MODULE improvedclouds_mod
     2
     3implicit none
     4
     5CONTAINS
     6
     7!=======================================================================
     8
     9SUBROUTINE improvedclouds(ngrid,nlay,ptimestep,pplay,pt,pdt,pq,pdq,nq,tauscaling,imicro,zt,zq)
     10
     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
     23use write_output_mod, only: write_output
     24use callkeys_mod,     only: activice, scavenging, cloud_adapt_ts, hdo, hdofrac
     25
     26implicit none
     27
     28!------------------------------------------------------------------
     29!  This routine is used to form clouds when a parcel of the GCM is
     30!    saturated. It includes the ability to have supersaturation, a
     31!    computation of the nucleation rates, growthrates and the
     32!    scavenging of dust particles by clouds.
     33!  It is worth noting that the amount of dust is computed using the
     34!    dust optical depth computed in aeropacity.F. That's why
     35!    the variable called "tauscaling" is used to convert
     36!    pq(dust_mass) and pq(dust_number), which are relative
     37!    quantities, to absolute and realistic quantities stored in zq.
     38!    This has to be done to convert the inputs into absolute
     39!    values, but also to convert the outputs back into relative
     40!    values which are then used by the sedimentation and advection
     41!    schemes.
     42
     43!  Authors: J.-B. Madeleine, based on the work by Franck Montmessin
     44!           (October 2011)
     45!           T. Navarro, debug,correction, new scheme (October-April 2011)
     46!           A. Spiga, optimization (February 2012)
     47!           J. Naar, adaptative subtimestep now done here (June 2023)
     48!------------------------------------------------------------------
     49!     Inputs/outputs:
     50INTEGER,                        INTENT(IN) :: ngrid,nlay
     51INTEGER,                        INTENT(IN) :: nq         ! nombre de traceurs
     52REAL,                           INTENT(IN) :: ptimestep  ! pas de temps physique (s)
     53REAL, dimension(ngrid,nlay),    INTENT(IN) :: pplay      ! pression au milieu des couches (Pa)
     54REAL, dimension(ngrid,nlay),    INTENT(IN) :: pt         ! temperature at the middle of the  layers (K)
     55REAL, dimension(ngrid,nlay),    INTENT(IN) :: pdt        ! temperature tendency (K/s)
     56REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pq         ! tracer (kg/kg)
     57REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pdq        ! tracer tendency (kg/kg/s)
     58REAL, dimension(ngrid),         INTENT(IN) :: tauscaling ! Convertion factor for qdust and Ndust
     59INTEGER,                        INTENT(IN) :: imicro     ! nb. microphy calls(retrocompatibility)
     60
     61REAL, dimension(ngrid,nlay,nq), INTENT(OUT) :: zq ! tracers post microphy (kg/kg)
     62REAL, dimension(ngrid,nlay),    INTENT(OUT) :: zt ! temperature post microphy (K)
     63
     64!------------------------------------------------------------------
     65!     Local variables:
     66LOGICAL, SAVE               ::  firstcall = .true.
    7367!$OMP THREADPRIVATE(firstcall)
    74 
    75       REAL*8   derf ! Error function
    76       !external derf
    77      
    78       INTEGER ig,l,i
    79      
    80       REAL zqsat(ngrid,nlay)    ! saturation
    81       REAL lw                   !Latent heat of sublimation (J.kg-1)
    82       REAL cste
    83       REAL dMice           ! mass of condensed ice
    84       REAL dMice_hdo       ! mass of condensed HDO ice
    85       REAL alpha(ngrid,nlay) ! HDO equilibrium fractionation coefficient (Saturation=1)
    86       REAL alpha_c(ngrid,nlay) ! HDO real fractionation coefficient
    87 !      REAL sumcheck
    88       REAL*8 ph2o          ! Water vapor partial pressure (Pa)
    89       REAL*8 satu          ! Water vapor saturation ratio over ice
    90       REAL*8 Mo,No
    91       REAL*8 Rn, Rm, dev2, n_derf, m_derf
    92       REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
    93       REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
    94 
    95       REAL dN,dM
    96       REAL rate(nbin_cld)  ! nucleation rate
    97       REAL seq
    98 
    99       REAL rice(ngrid,nlay)      ! Ice mass mean radius (m)
    100                                  ! (r_c in montmessin_2004)
    101       REAL rhocloud(ngrid,nlay)  ! Cloud density (kg.m-3)
    102       REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
    103 
    104       REAL res      ! Resistance growth
    105       REAL Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient
    106      
    107 
    108 c     Parameters of the size discretization
    109 c       used by the microphysical scheme
    110       DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
    111       DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
    112       DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
    113                                            ! Minimum boundary radius (m)
    114       DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
    115       DOUBLE PRECISION vrat_cld ! Volume ratio
    116       DOUBLE PRECISION, SAVE :: rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
     68REAL*8                      :: derf ! Error function
     69!external derf
     70INTEGER                     :: ig, l, i
     71REAL, dimension(ngrid,nlay) :: zqsat ! Saturation
     72REAL                        :: lw    ! Latent heat of sublimation (J.kg-1)
     73REAL                        :: cste
     74REAL                        :: dMice     ! mass of condensed ice
     75REAL                        :: dMice_hdo ! mass of condensed HDO ice
     76REAL, dimension(ngrid,nlay) :: alpha     ! HDO equilibrium fractionation coefficient (Saturation=1)
     77REAL, dimension(ngrid,nlay) :: alpha_c   ! HDO real fractionation coefficient
     78! REAL :: sumcheck
     79REAL*8                      :: ph2o ! Water vapor partial pressure (Pa)
     80REAL*8                      :: satu ! Water vapor saturation ratio over ice
     81REAL*8                      :: Mo, No, Rn, Rm, dev2, n_derf, m_derf
     82REAL*8, dimension(nbin_cld) :: n_aer ! number conc. of particle/each size bin
     83REAL*8, dimension(nbin_cld) :: m_aer ! mass mixing ratio of particle/each size bin
     84REAL :: dN, dM, seq
     85REAL, dimension(nbin_cld)   :: rate ! Nucleation rate
     86REAL, dimension(ngrid,nlay) :: rice ! Ice mass mean radius (m)
     87                                    ! (r_c in montmessin_2004)
     88REAL, dimension(ngrid,nlay) :: rhocloud   ! Cloud density (kg.m-3)
     89REAL, dimension(ngrid,nlay) :: rdust      ! Dust geometric mean radius (m)
     90REAL                        :: res        ! Resistance growth
     91REAL                        :: Dv, Dv_hdo ! Water/HDO vapor diffusion coefficient
     92
     93! Parameters of the size discretization used by the microphysical scheme
     94DOUBLE PRECISION, PARAMETER                     :: rmin_cld = 0.1e-6     ! Minimum radius (m)
     95DOUBLE PRECISION, PARAMETER                     :: rmax_cld = 10.e-6     ! Maximum radius (m)
     96DOUBLE PRECISION, PARAMETER                     :: rbmin_cld = 0.0001e-6 ! Minimum boundary radius (m)
     97DOUBLE PRECISION, PARAMETER                     :: rbmax_cld = 1.e-2     ! Maximum boundary radius (m)
     98DOUBLE PRECISION                                :: vrat_cld ! Volume ratio
     99DOUBLE PRECISION, dimension(nbin_cld + 1), SAVE :: rb_cld   ! Boundary values of each rad_cld bin (m)
    117100!$OMP THREADPRIVATE(rb_cld)
    118 
    119       DOUBLE PRECISION dr_cld(nbin_cld)   ! width of each rad_cld bin (m)
    120       DOUBLE PRECISION vol_cld(nbin_cld)  ! particle volume for each bin (m3)
    121 
    122       REAL, SAVE :: sigma_ice ! Variance of the ice and CCN distributions
     101DOUBLE PRECISION, dimension(nbin_cld) :: dr_cld    ! Width of each rad_cld bin (m)
     102DOUBLE PRECISION, dimension(nbin_cld) :: vol_cld   ! Particle volume for each bin (m3)
     103REAL, SAVE                            :: sigma_ice ! Variance of the ice and CCN distributions
    123104!$OMP THREADPRIVATE(sigma_ice)
    124105
    125 
    126 c----------------------------------     
    127 c JN : used in subtimestep
    128       REAL :: microtimestep! subdivision of ptimestep (s)
    129       REAL :: subpdtcloud(ngrid,nlay)  ! Temperature variation due to latent heat
    130       REAL :: zq0(ngrid,nlay,nq) ! local initial value of tracers
    131 c      REAL zt0(ngrid,nlay,nq) ! local initial value of temperature
    132       INTEGER :: zimicro(ngrid,nlay) ! Subdivision of ptimestep
    133       INTEGER :: count_micro(ngrid,nlay) ! Number of microphys calls
    134       REAL :: zpotcond_inst(ngrid,nlay) ! condensable water at the beginning of physics (kg/kg)
    135       REAL :: zpotcond_full(ngrid,nlay) ! condensable water with integrated tendancies (kg/kg)
    136       REAL :: zpotcond(ngrid,nlay) ! maximal condensable water (previous two)
    137       REAL :: zqsat_tmp(1) ! maximal condensable water (previous two)
    138       REAL :: spenttime ! time spent in while loop
    139       REAL :: zdq ! used to compute adaptative timestep
    140       LOGICAL :: ending_ts ! Condition to end while loop
    141 
    142      
    143 c----------------------------------     
    144 c TESTS
    145 
    146       INTEGER countcells
    147      
    148       LOGICAL, SAVE :: test_flag    ! flag for test/debuging outputs
     106! JN: used in subtimestep
     107REAL                              :: microtimestep ! Subdivision of ptimestep (s)
     108REAL,    dimension(ngrid,nlay)    :: subpdtcloud   ! Temperature variation due to latent heat
     109REAL,    dimension(ngrid,nlay,nq) :: zq0           ! Local initial value of tracers
     110!REAL,    dimension(ngrid,nlay,nq) :: zt0           ! Local initial value of temperature
     111INTEGER, dimension(ngrid,nlay)    :: zimicro       ! Subdivision of ptimestep
     112INTEGER, dimension(ngrid,nlay)    :: count_micro   ! Number of microphys calls
     113REAL,    dimension(ngrid,nlay)    :: zpotcond_inst ! Condensable water at the beginning of physics (kg/kg)
     114REAL,    dimension(ngrid,nlay)    :: zpotcond_full ! Condensable water with integrated tendancies (kg/kg)
     115REAL,    dimension(ngrid,nlay)    :: zpotcond      ! Maximal condensable water (previous two)
     116REAL,    dimension(1)             :: zqsat_tmp     ! Maximal condensable water (previous two)
     117REAL                              :: spenttime     ! Time spent in while loop
     118REAL                              :: zdq           ! Used to compute adaptative timestep
     119LOGICAL                           :: ending_ts     ! Condition to end while loop
     120
     121!----------------------------------
     122! TESTS
     123INTEGER       :: countcells
     124LOGICAL, SAVE :: test_flag ! Flag for test/debuging outputs
    149125!$OMP THREADPRIVATE(test_flag)
    150 
    151 
    152       REAL satubf(ngrid,nlay),satuaf(ngrid,nlay)
    153       REAL res_out(ngrid,nlay)
    154  
    155 
    156 c------------------------------------------------------------------
    157 
    158       ! AS: firstcall OK absolute
    159       IF (firstcall) THEN
     126REAL, dimension(ngrid,nlay) :: satubf, satuaf, res_out
     127
     128!------------------------------------------------------------------
     129
     130! AS: firstcall OK absolute
     131IF (firstcall) THEN
    160132!=============================================================
    161133! 0. Definition of the size grid
    162134!=============================================================
    163 c       rad_cld is the primary radius grid used for microphysics computation.
    164 c       The grid spacing is computed assuming a constant volume ratio
    165 c       between two consecutive bins; i.e. vrat_cld.
    166 c       vrat_cld is determined from the boundary values of the size grid:
    167 c       rmin_cld and rmax_cld.
    168 c       The rb_cld array contains the boundary values of each rad_cld bin.
    169 c       dr_cld is the width of each rad_cld bin.
    170 
    171 c       Volume ratio between two adjacent bins
    172 !        vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
    173 !        vrat_cld = exp(vrat_cld)
    174         vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
    175         vrat_cld = exp(vrat_cld)
    176         write(*,*) "vrat_cld", vrat_cld
    177 
    178         rb_cld(1)  = rbmin_cld
    179         rad_cld(1) = rmin_cld
    180         vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
    181 !        vol_cld(1) = 4./3. * pi * rmin_cld*rmin_cld*rmin_cld
    182 
    183         do i=1,nbin_cld-1
    184           rad_cld(i+1)  = rad_cld(i) * vrat_cld**(1./3.)
    185           vol_cld(i+1)  = vol_cld(i) * vrat_cld
    186         enddo
    187        
    188         do i=1,nbin_cld
    189           rb_cld(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) *
    190      &      rad_cld(i)
    191           dr_cld(i)  = rb_cld(i+1) - rb_cld(i)
    192         enddo
    193         rb_cld(nbin_cld+1) = rbmax_cld
    194         dr_cld(nbin_cld)   = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
    195 
    196         print*, ' '
    197         print*,'Microphysics: size bin information:'
    198         print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
    199         print*,'-----------------------------------'
    200         do i=1,nbin_cld
    201           write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
    202      &      dr_cld(i)
    203         enddo
    204         write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
    205         print*,'-----------------------------------'
    206 
    207         do i=1,nbin_cld+1
    208 !           rb_cld(i) = log(rb_cld(i)) 
    209             rb_cld(i) = log(rb_cld(i))  !! we save that so that it is not computed
    210                                          !! at each timestep and gridpoint
    211         enddo
    212 
    213 c       Contact parameter of water ice on dust ( m=cos(theta) )
    214 c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    215 !       mteta is initialized in conf_phys
    216         write(*,*) 'water_param contact parameter:', mteta
    217 
    218 c       Volume of a water molecule (m3)
    219         vo1 = mh2o / dble(rho_ice)
    220 c       Variance of the ice and CCN distributions
    221         sigma_ice = sqrt(log(1.+nuice_sed))
    222        
    223  
    224         write(*,*) 'Variance of ice & CCN distribs :', sigma_ice
    225         write(*,*) 'nuice for sedimentation:', nuice_sed
    226         write(*,*) 'Volume of a water molecule:', vo1
    227 
    228 
    229         test_flag = .false.
    230  
    231         firstcall=.false.
    232       END IF
    233 
     135! 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.
     138! The rb_cld array contains the boundary values of each rad_cld bin.
     139! dr_cld is the width of each rad_cld bin.
     140
     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.
     194ENDIF
    234195
    235196!=============================================================
    236197! 1. Initialisation
    237198!=============================================================
    238 
    239       res_out(:,:) = 0
    240       rice(:,:) = 1.e-8
    241 
    242 c     Initialize the temperature, tracers
    243       zt(1:ngrid,1:nlay)=pt(1:ngrid,1:nlay)
    244       zq(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)
    245       subpdtcloud(1:ngrid,1:nlay)=0
    246      
    247       WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
    248      &       zq(1:ngrid,1:nlay,1:nq) = 1.e-30
    249 
    250      
     199res_out(:,:) = 0
     200rice(:,:) = 1.e-8
     201
     202! Initialize the temperature, tracers
     203zt(:,:) = pt(:,:)
     204zq(:,:,:) = pq(:,:,:)
     205subpdtcloud(:,:) = 0
     206
     207WHERE ( zq(:,:,:) < 1.e-30 ) zq(:,:,:) = 1.e-30
     208
    251209!=============================================================
    252210! 2. Compute saturation
    253211!=============================================================
    254 
    255 
    256       dev2 = 1. / ( sqrt(2.) * sigma_ice )
    257 
    258 c     Compute the condensable amount of water vapor
    259 c     or the sublimable water ice (if negative)
    260       call watersat(ngrid*nlay,zt+pdt*ptimestep,pplay,zqsat)
    261       zpotcond_full=(zq(:,:,igcm_h2o_vap)+
    262      &             pdq(:,:,igcm_h2o_vap)*ptimestep) - zqsat
    263       zpotcond_full=max(zpotcond_full,-zq(:,:,igcm_h2o_ice))
    264       call watersat(ngrid*nlay,zt,pplay,zqsat)
    265       zpotcond_inst=zq(:,:,igcm_h2o_vap) - zqsat
    266       zpotcond_inst=max(zpotcond_inst,-zq(:,:,igcm_h2o_ice))
    267 c     zpotcond is the most strict criterion between the two
    268       DO l=1,nlay
    269         DO ig=1,ngrid
    270           if (zpotcond_full(ig,l).gt.0.) then
    271             zpotcond(ig,l)=max(zpotcond_inst(ig,l),zpotcond_full(ig,l))
    272           else if (zpotcond_full(ig,l).le.0.) then
    273             zpotcond(ig,l)=min(zpotcond_inst(ig,l),zpotcond_full(ig,l))
    274           endif! (zpotcond_full.gt.0.) then
    275         ENDDO !ig=1,ngrid
    276       ENDDO !l=1,nlay
    277            
    278       countcells = 0
    279       zimicro(1:ngrid,1:nlay)=imicro
    280       count_micro(1:ngrid,1:nlay)=0
    281 
    282 c     Main loop over the GCM's grid
    283       DO l=1,nlay
    284         DO ig=1,ngrid
    285 c       Subtimestep : here we go
    286         IF (cloud_adapt_ts) then
    287                 call adapt_imicro(ptimestep,zpotcond(ig,l),
    288      &             zimicro(ig,l))
    289         ENDIF! (cloud_adapt_ts) then
     212dev2 = 1. / ( sqrt(2.) * sigma_ice )
     213
     214! Compute the condensable amount of water vapor or the sublimable water ice (if negative)
     215call watersat(ngrid*nlay,zt+pdt*ptimestep,pplay,zqsat)
     216zpotcond_full = (zq(:,:,igcm_h2o_vap)+pdq(:,:,igcm_h2o_vap)*ptimestep) - zqsat
     217zpotcond_full = max(zpotcond_full,-zq(:,:,igcm_h2o_ice))
     218call watersat(ngrid*nlay,zt,pplay,zqsat)
     219zpotcond_inst = zq(:,:,igcm_h2o_vap) - zqsat
     220zpotcond_inst = max(zpotcond_inst,-zq(:,:,igcm_h2o_ice))
     221! zpotcond is the most strict criterion between the two
     222DO 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
     230ENDDO !l=1,nlay
     231
     232countcells = 0
     233zimicro(1:ngrid,1:nlay)=imicro
     234count_micro(:,:)=0
     235
     236! Main loop over the GCM's grid
     237DO 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))
    290241        spenttime = 0.
    291242        ending_ts=.false.
    292243        DO while (.not.ending_ts)
    293           microtimestep=ptimestep/real(zimicro(ig,l))
    294 c         Initialize tracers for scavenging + hdo computations (JN)
    295           DO i=1,nq
    296              zq0(ig,l,i)=zq(ig,l,i)
    297           ENDDO !i=1,nq
    298 
    299           ! Check if we are integrating over ptimestep
    300           if (spenttime+microtimestep.ge.ptimestep) then
    301                   microtimestep=ptimestep-spenttime
    302           !       If so : last call !
    303                   ending_ts=.true.
    304           endif! (spenttime+microtimestep.ge.ptimestep) then
    305           call watersat(1,(/zt(ig,l)/),(/pplay(ig,l)/),zqsat_tmp)
    306           zqsat(ig,l)=zqsat_tmp(1)
    307 c       Get the partial pressure of water vapor and its saturation ratio
    308         ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
    309         satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
     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)
    310259
    311260!=============================================================
    312261! 3. Nucleation
    313262!=============================================================
    314 
    315        IF ( satu .ge. 1. ) THEN         ! if there is condensation
    316 
    317         call updaterccn(zq(ig,l,igcm_dust_mass),
    318      &          zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
    319 
    320 
    321 c       Expand the dust moments into a binned distribution
    322         Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)   + 1.e-30
    323         No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
    324         Rn = rdust(ig,l)
    325         Rn = -log(Rn)
    326         Rm = Rn - 3. * sigma_ice*sigma_ice 
    327         n_derf = derf( (rb_cld(1)+Rn) *dev2)
    328         m_derf = derf( (rb_cld(1)+Rm) *dev2)
    329         do i = 1, nbin_cld
    330           n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
    331           m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
    332           n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
    333           m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
    334           n_aer(i) = n_aer(i) + 0.5 * No * n_derf
    335           m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
    336         enddo
    337        
    338 !        sumcheck = 0
    339 !        do i = 1, nbin_cld
    340 !          sumcheck = sumcheck + n_aer(i)
    341 !        enddo
    342 !        sumcheck = abs(sumcheck/No - 1)
    343 !        if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
    344 !          print*, "WARNING, No sumcheck PROBLEM"
    345 !          print*, "sumcheck, No",sumcheck, No
    346 !          print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
    347 !          print*, "Dust binned distribution", n_aer
    348 !        endif
    349 !       
    350 !        sumcheck = 0
    351 !        do i = 1, nbin_cld
    352 !          sumcheck = sumcheck + m_aer(i)
    353 !        enddo
    354 !        sumcheck = abs(sumcheck/Mo - 1)
    355 !        if ((sumcheck .gt. 1e-5) .and.  (1./Rn .gt. rmin_cld)) then
    356 !          print*, "WARNING, Mo sumcheck PROBLEM"
    357 !          print*, "sumcheck, Mo",sumcheck, Mo
    358 !          print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
    359 !          print*, "Dust binned distribution", m_aer
    360 !        endif
    361 
    362  
    363 c       Get the rates of nucleation
    364         call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
    365        
    366         dN = 0.
    367         dM = 0.
    368         do i = 1, nbin_cld
    369           dN       = dN + n_aer(i)*(exp(-rate(i)*microtimestep)-1.)
    370           dM       = dM + m_aer(i)*(exp(-rate(i)*microtimestep)-1.)
    371         enddo
    372 
    373 
    374 c       Update Dust particles
    375         zq(ig,l,igcm_dust_mass)   =
    376      &  zq(ig,l,igcm_dust_mass)   + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    377         zq(ig,l,igcm_dust_number) =
    378      &  zq(ig,l,igcm_dust_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    379 c       Update CCNs
    380         zq(ig,l,igcm_ccn_mass)   =
    381      &  zq(ig,l,igcm_ccn_mass)   - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    382         zq(ig,l,igcm_ccn_number) =
    383      &  zq(ig,l,igcm_ccn_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    384        
    385         ENDIF ! of is satu >1
     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
    386325
    387326!=============================================================
    388327! 4. Ice growth: scheme for radius evolution
    389328!=============================================================
    390 
    391 c We trigger crystal growth if and only if there is at least one nuclei (N>1).
    392 c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
    393 c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
    394 
    395        IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
    396 
    397         call updaterice_micro(zq(ig,l,igcm_h2o_ice),
    398      &    zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
    399      &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
    400 
    401         No   = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
    402 
    403 c       saturation at equilibrium
    404 c       rice should not be too small, otherwise seq value is not valid
    405         seq  = exp(2.*sig_h2o(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)*
    406      &             max(rice(ig,l),1.e-7)))
    407        
    408 c       get resistance growth
    409         call growthrate(zt(ig,l),pplay(ig,l),
    410      &             real(ph2o/satu),rice(ig,l),res,Dv)
    411 
    412         res_out(ig,l) = res
    413 
    414 ccccccc  implicit scheme of mass growth
    415 c         cste here must be computed at each step
    416         cste = 4*pi*rho_ice*microtimestep
    417 
    418         dMice =
    419      & (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))
    420      &   /(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.)
    421 
    422 
    423 ! With the above scheme, dMice cannot be bigger than vapor,
    424 ! but can be bigger than all available ice.
    425        dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
    426        dMice = min(dMice,zq(ig,l,igcm_h2o_vap)) ! this should be useless...
    427 
    428        zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
    429        zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
    430 
    431 
    432        countcells = countcells + 1
    433        
    434        ! latent heat release
    435        lw=(2834.3-0.28*(zt(ig,l)-To)-
    436      &     0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3
    437        subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
    438          
    439 c Special case of the isotope of water HDO   
    440        if (hdo) then
    441          !! condensation
    442          if (dMice.gt.0.0) then
    443          !! do we use fractionation?
    444            if (hdofrac) then
    445              !! Calculation of the HDO vapor coefficient
    446              Dv_hdo = 1./3. * sqrt( 8*kbz*zt(ig,l)/(pi*mhdo/nav) )
    447      &          * kbz * zt(ig,l) /
    448      &          ( pi * pplay(ig,l) * (molco2+molhdo)*(molco2+molhdo)
    449      &          * sqrt(1.+mhdo/mco2) )
    450              !! Calculation of the fractionnation coefficient at equilibrium
    451              alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
    452 c             alpha = exp(13525./zt(ig,l)**2.-5.59d-2)  !Lamb
    453              !! Calculation of the 'real' fractionnation coefficient
    454              alpha_c(ig,l) = (alpha(ig,l)*satu)/
    455      &          ( (alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.)
    456 c             alpha_c(ig,l) = alpha(ig,l) ! to test without the effect of cinetics
    457            else
    458              alpha_c(ig,l) = 1.d0
    459            endif
    460            if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then
    461               dMice_hdo=
    462      &          dMice*alpha_c(ig,l)*
    463      &     ( zq0(ig,l,igcm_hdo_vap)
    464      &             /zq0(ig,l,igcm_h2o_vap) )
    465            else
    466              dMice_hdo=0.
    467            endif
    468          !! sublimation
    469          else
    470            if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then
    471              dMice_hdo=
    472      &            dMice*
    473      &     ( zq0(ig,l,igcm_hdo_ice)
    474      &             /zq0(ig,l,igcm_h2o_ice) )
    475            else
    476              dMice_hdo=0.
    477            endif
    478          endif !if (dMice.gt.0.0)
    479 
    480        dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice))
    481        dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap))
    482 
    483        zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo
    484        zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo
    485 
    486        endif ! if (hdo)       
    487      
     329! We trigger crystal growth if and only if there is at least one nuclei (N>1).
     330! Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
     331! 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)
     405
    488406!=============================================================
    489407! 5. Dust cores released, tendancies, latent heat, etc ...
    490408!=============================================================
    491 
    492 c         If all the ice particles sublimate, all the condensation
    493 c         nuclei are released:
    494           if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
    495          
    496 c           Water
    497             zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
    498      &                            + zq(ig,l,igcm_h2o_ice)
    499             zq(ig,l,igcm_h2o_ice) = 0.
    500             if (hdo) then
    501               zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)
    502      &                            + zq(ig,l,igcm_hdo_ice)
    503               zq(ig,l,igcm_hdo_ice) = 0.
    504             endif
    505 c           Dust particles
    506             zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
    507      &                              + zq(ig,l,igcm_ccn_mass)
    508             zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
    509      &                                + zq(ig,l,igcm_ccn_number)
    510 c           CCNs
    511             zq(ig,l,igcm_ccn_mass) = 0.
    512             zq(ig,l,igcm_ccn_number) = 0.
    513 
    514           endif
    515          
    516           ELSE
    517           ! Initialization of dMice when it's not computed
    518             dMice=0 ! no condensation/sublimation to account for
    519           ENDIF !of if Nccn>1
    520 
    521          
    522       ! No more getting tendency : we increment tracers & temp directly
    523 
    524       ! Increment tracers & temp for following microtimestep
    525       ! Dust :
    526       ! Special treatment of dust_mass / number for scavenging ?
     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?
    527437            IF (scavenging) THEN
    528               zq(ig,l,igcm_dust_mass) =zq(ig,l,igcm_dust_mass)+
    529      &         pdq(ig,l,igcm_dust_mass)*microtimestep
    530               zq(ig,l,igcm_dust_number) =zq(ig,l,igcm_dust_number)+
    531      &         pdq(ig,l,igcm_dust_number)*microtimestep
     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
    532440            ELSE
    533               zq(ig,l,igcm_dust_mass) =zq0(ig,l,igcm_dust_mass)
    534               zq(ig,l,igcm_dust_number) =zq0(ig,l,igcm_dust_number)
     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)
    535443            ENDIF !(scavenging) THEN
    536               zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) +
    537      &         pdq(ig,l,igcm_ccn_mass)*microtimestep
    538               zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) +
    539      &          pdq(ig,l,igcm_ccn_number)*microtimestep
    540 
    541       ! Water :
    542             zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+
    543      &         pdq(ig,l,igcm_h2o_ice)*microtimestep
    544             zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)+
    545      &         pdq(ig,l,igcm_h2o_vap)*microtimestep
    546 
    547             ! HDO (if computed) :
     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):
    548452            IF (hdo) THEN
    549             zq(ig,l,igcm_hdo_ice) =
    550      &       zq(ig,l,igcm_hdo_ice)+
    551      &         pdq(ig,l,igcm_hdo_ice)*microtimestep
    552             zq(ig,l,igcm_hdo_vap) =
    553      &       zq(ig,l,igcm_hdo_vap)+
    554      &         pdq(ig,l,igcm_hdo_vap)*microtimestep
     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
    555455            ENDIF ! hdo
    556456
    557 c  Could also set subpdtcloud to 0 if not activice to make it simpler
    558 c  or change name of the flag
    559             IF (.not.activice) THEN
    560                     subpdtcloud(ig,l)=0.
    561             ENDIF
    562       ! Temperature :
    563             zt(ig,l) = zt(ig,l)+(pdt(ig,l)+
    564      &          subpdtcloud(ig,l))*microtimestep
    565 
    566 c         Prevent negative tracers ! JN
    567           DO i=1,nq
    568             IF(zq(ig,l,i).lt.1.e-30) zq(ig,l,i)=1.e-30
    569           ENDDO !i=1,nq
    570 
    571           IF (cloud_adapt_ts) then
    572 c           Estimation of how much is actually condensing/subliming
    573                 IF (spenttime.ne.0) then
    574                    zdq=(dMice/spenttime)*(ptimestep-spenttime)
    575                 ELSE
    576                 !Initialization for spenttime=0
    577                    zdq=zpotcond(ig,l)*((ptimestep-spenttime)/
    578      &          ptimestep)
    579                 ENDIF
    580 c           Threshold with powerlaw (sanity check)
    581                 zdq=min(abs(zdq),
    582      &            abs(zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep)))
    583                 call adapt_imicro(ptimestep,zdq,
    584      &             zimicro(ig,l))
    585           ENDIF! (cloud_adapt_ts) then
    586 c         Increment time spent in here
    587           spenttime=spenttime+microtimestep
    588           count_micro(ig,l)=count_micro(ig,l)+1
    589           ENDDO ! while (.not. ending_ts)
    590         ENDDO ! of ig loop
    591       ENDDO ! of nlayer loop
    592      
    593      
    594 c------ Useful outputs to check how it went
    595       call write_output("zpotcond_inst","zpotcond_inst "//
    596      &   "microphysics","(kg/kg)",zpotcond_inst(:,:))
    597       call write_output("zpotcond_full","zpotcond_full "//
    598      &   "microphysics","(kg/kg)",zpotcond_full(:,:))
    599       call write_output("zpotcond","zpotcond "//
    600      &   "microphysics","(kg/kg)",zpotcond(:,:))
    601       call write_output("count_micro","count_micro "//
    602      &   "after microphysics","integer",count_micro(:,:))
    603      
    604      
    605      
    606 !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
    607 !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
     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
     483ENDDO ! of nlayer loop
     484
     485!------ Useful outputs to check how it went
     486call write_output("zpotcond_inst","zpotcond_inst "//"microphysics","(kg/kg)",zpotcond_inst(:,:))
     487call write_output("zpotcond_full","zpotcond_full "//"microphysics","(kg/kg)",zpotcond_full(:,:))
     488call write_output("zpotcond","zpotcond "//"microphysics","(kg/kg)",zpotcond(:,:))
     489call write_output("count_micro","count_micro "//"after microphysics","integer",count_micro(:,:))
     490
    608491!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
    609       IF (test_flag) then
    610      
    611 !       error2d(:) = 0.
    612        DO l=1,nlay
    613        DO ig=1,ngrid
    614 !         error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
    615           satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
    616           satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
    617        ENDDO
    618        ENDDO
    619 
    620       print*, 'count is ',countcells, ' i.e. ',
    621      &     countcells*100/(nlay*ngrid), '% for microphys computation'
    622      
    623       ENDIF ! endif test_flag
    624 
    625      
    626 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc     
    627 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc     
    628 c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
    629 c It is an analytical solution to the ice radius growth equation,
    630 c with the approximation of a constant 'reduced' cunningham correction factor
    631 c (lambda in growthrate.F) taken at radius req instead of rice   
    632 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc     
    633 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    634 
    635 c      subroutine phi(rice,req,coeff1,coeff2,time)
    636 c     
    637 c      implicit none
    638 c     
    639 c      ! inputs
    640 c      real rice ! ice radius
    641 c      real req  ! ice radius at equilibirum
    642 c      real coeff1  ! coeff for the log
    643 c      real coeff2  ! coeff for the arctan
    644 c
    645 c      ! output     
    646 c      real time
    647 c     
    648 c      !local
    649 c      real var
    650 c     
    651 c      ! 1.73205 is sqrt(3)
    652 c     
    653 c      var = max(
    654 c     &  abs(rice-req) / sqrt(rice*rice + rice*req  + req*req),1e-30)
    655 c           
    656 c       time =
    657 c     &   coeff1 *
    658 c     &   log( var )
    659 c     & + coeff2 * 1.73205 *
    660 c     &   atan( (2*rice+req) / (1.73205*req) )
    661 c     
    662 c      return
    663 c      end
    664      
    665      
    666      
    667       END SUBROUTINE improvedclouds
    668 
    669       SUBROUTINE adapt_imicro(ptimestep,potcond,
    670      $                     zimicro)
    671 
    672 c Adaptative timestep for water ice clouds.
    673 c Works using a powerlaw to compute the minimal duration of subtimestep
    674 c (in s) should all the avalaible vapor (resp. ice) condenses (resp.sublimates)
    675 c Then, we use the instantaneous vap (ice) gradient extrapolated to the
    676 c rest of duration to increase subtimestep duration, for computing
    677 c efficiency.
    678 
    679       real,intent(in) :: ptimestep ! total duration of physics (sec)
    680       real,intent(in) :: potcond ! condensible vapor / sublimable ice(kg/kg)
    681       real :: alpha, beta ! Coefficients for power law
    682       real :: defstep,coef ! Default ptimestep of 7.5 mins (iphysiq=5)
    683       integer,intent(out) :: zimicro ! number of ptimestep division
    684 
    685 c      Default ptimestep : defstep (7.5 mins)
    686        defstep=88775.*5./960.
    687        coef=ptimestep/defstep
    688 c      Conservative coefficients :
    689        alpha=1.81846504e+11
    690        beta=1.54550140e+00
    691        zimicro=ceiling(coef*min(max(alpha*abs(potcond)**beta,5.),7000.))
    692 
    693       END SUBROUTINE adapt_imicro
    694      
    695       END MODULE improvedclouds_mod
     492IF (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
     500    ENDDO
     501    write(*,*) 'count is ',countcells, ' i.e. ',countcells*100/(nlay*ngrid), '% for microphys computation'
     502ENDIF ! endif test_flag
     503
     504END SUBROUTINE improvedclouds
     505
     506!=======================================================================
     507
     508!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     509!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     510! 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
     516!
     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!
     530! !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
     540
     541!=======================================================================
     542
     543SUBROUTINE adapt_imicro(ptimestep,potcond,zimicro)
     544
     545! Adaptative timestep for water ice clouds.
     546! Works using a powerlaw to compute the minimal duration of subtimestep
     547! (in s) should all the avalaible vapor (resp. ice) condenses (resp.sublimates)
     548! Then, we use the instantaneous vap (ice) gradient extrapolated to the
     549! rest of duration to increase subtimestep duration, for computing
     550! efficiency.
     551
     552real, intent(in) :: ptimestep ! total duration of physics (sec)
     553real, intent(in) :: potcond ! condensible vapor / sublimable ice(kg/kg)
     554integer, intent(out) :: zimicro ! number of ptimestep division
     555real :: alpha, beta ! Coefficients for power law
     556real :: defstep,coef ! Default ptimestep of 7.5 mins (iphysiq=5)
     557
     558! Default ptimestep: defstep (7.5 mins)
     559defstep = 88775.*5./960.
     560coef = ptimestep/defstep
     561! Conservative coefficients :
     562alpha = 1.81846504e+11
     563beta = 1.54550140e+00
     564zimicro = ceiling(coef*min(max(alpha*abs(potcond)**beta,5.),7000.))
     565
     566END SUBROUTINE adapt_imicro
     567
     568END MODULE improvedclouds_mod
Note: See TracChangeset for help on using the changeset viewer.