Ignore:
Timestamp:
Apr 14, 2026, 10:31:30 AM (4 weeks ago)
Author:
jbclement
Message:

Mars PCM:
Cosmetic cleaning of "improvedclouds".
JBC

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r4161 r4186  
    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
    … …  
    5050!------------------------------------------------------------------
    5151!     Inputs/outputs:
    52 INTEGER, INTENT(IN) :: ngrid,nlay
    53 INTEGER, INTENT(IN) :: nq ! nombre de traceurs
    54 REAL, INTENT(IN) :: ptimestep ! pas de temps physique (s)
    55 REAL, dimension(ngrid,nlay), INTENT(IN) :: pplay ! pression au milieu des couches (Pa)           
    56 REAL, dimension(ngrid,nlay), INTENT(IN) :: pt ! temperature at the middle of the layers (K)
    57 REAL, dimension(ngrid,nlay), INTENT(IN) :: pdt ! temperature tendency (K/s)
    58 REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pq ! tracer (kg/kg)
    59 REAL, dimension(ngrid,nlay,nq), INTENT(IN) :: pdq ! tracer tendency (kg/kg/s)
    60 REAL, dimension(ngrid), INTENT(IN) :: tauscaling ! Convertion factor for qdust and Ndust
    61 INTEGER, INTENT(IN) :: imicro ! nb. microphy calls(retrocompatibility)
    62 
    63 REAL, dimension(ngrid,nlay,nq), INTENT(OUT) :: zq ! tracers post microphy (kg/kg) 
    64 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)
    6565
    6666!------------------------------------------------------------------
    6767!     Local variables:
    68 LOGICAL, SAVE ::  firstcall = .true.
     68LOGICAL, SAVE               ::  firstcall = .true.
    6969!$OMP THREADPRIVATE(firstcall)
    70 REAL*8 :: derf ! Error function
     70REAL*8                      :: derf ! Error function
    7171!external derf
    72 INTEGER :: ig,l,i
    73 REAL, dimension(ngrid,nlay) :: zqsat ! saturation
    74 REAL :: lw ! Latent heat of sublimation (J.kg-1)
    75 REAL :: cste
    76 REAL :: dMice ! mass of condensed ice
    77 REAL :: dMicetot ! JN
    78 REAL :: dMice_hdo ! mass of condensed HDO ice
    79 REAL, dimension(ngrid,nlay) :: alpha ! HDO equilibrium fractionation coefficient (Saturation=1)
    80 REAL, dimension(ngrid,nlay) :: alpha_c ! HDO real fractionation coefficient
    81 REAL*8 :: ph2o ! Water vapor partial pressure (Pa)
    82 REAL*8 :: satu ! Water vapor saturation ratio over ice
    83 REAL*8 :: Mo,No, Rn, Rm, dev2, n_derf, m_derf
    84 REAL*8, dimension(nbin_cld) :: n_aer ! number conc. of particle/each size bin
    85 REAL*8, dimension(nbin_cld) :: m_aer ! mass mixing ratio of particle/each size bin
     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
     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
     81REAL*8                      :: ph2o      ! Water vapor partial pressure (Pa)
     82REAL*8                      :: satu      ! Water vapor saturation ratio over ice
     83REAL*8                      :: Mo, No, Rn, Rm, dev2, n_derf, m_derf
     84REAL*8, dimension(nbin_cld) :: n_aer      ! number conc. of particle/each size bin
     85REAL*8, dimension(nbin_cld) :: m_aer      ! mass mixing ratio of particle/each size bin
    8686REAL :: dN, dM, seq
    87 REAL, dimension(nbin_cld) :: rate ! nucleation rate
    88 REAL, dimension(ngrid,nlay) :: rice ! Ice mass mean radius (m) (r_c in montmessin_2004)
    89 REAL, dimension(ngrid,nlay) :: rhocloud ! Cloud density (kg.m-3)
    90 REAL, dimension(ngrid,nlay) :: rdust ! Dust geometric mean radius (m)
    91 REAL :: res ! Resistance growth
    92 REAL :: Dv,Dv_hdo ! Water/HDO vapor diffusion coefficient
     87REAL, dimension(nbin_cld)   :: rate       ! Nucleation rate
     88REAL, dimension(ngrid,nlay) :: rice       ! Ice mass mean radius (m) (r_c in montmessin_2004)
     89REAL, dimension(ngrid,nlay) :: rhocloud   ! Cloud density (kg.m-3)
     90REAL, dimension(ngrid,nlay) :: rdust      ! Dust geometric mean radius (m)
     91REAL                        :: res        ! Resistance growth
     92REAL                        :: Dv, Dv_hdo ! Water/HDO vapor diffusion coefficient
    9393
    9494! Parameters of the size discretization used by the microphysical scheme
    95 DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
    96 DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
    97 DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6 ! Minimum boundary radius (m)
    98 DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
    99 DOUBLE PRECISION :: vrat_cld ! Volume ratio
    100 DOUBLE PRECISION, dimension(nbin_cld+1), SAVE :: rb_cld ! boundary values of each rad_cld bin (m)
     95DOUBLE PRECISION, PARAMETER                     :: rmin_cld = 0.1e-6    ! Minimum radius (m)
     96DOUBLE PRECISION, PARAMETER                     :: rmax_cld = 10.e-6    ! Maximum radius (m)
     97DOUBLE PRECISION, PARAMETER                     :: rbmin_cld = 0.0001e-6 ! Minimum boundary radius (m)
     98DOUBLE PRECISION, PARAMETER                     :: rbmax_cld = 1.e-2    ! Maximum boundary radius (m)
     99DOUBLE PRECISION                                :: vrat_cld              ! Volume ratio
     100DOUBLE PRECISION, dimension(nbin_cld + 1), SAVE :: rb_cld                ! Boundary values of each rad_cld bin (m)
    101101!$OMP THREADPRIVATE(rb_cld)
    102 DOUBLE PRECISION, dimension(nbin_cld) :: dr_cld ! width of each rad_cld bin (m)
    103 DOUBLE PRECISION, dimension(nbin_cld) :: vol_cld ! particle volume for each bin (m3)
    104 REAL, SAVE :: sigma_ice ! Variance of the ice and CCN distributions
     102DOUBLE PRECISION, dimension(nbin_cld) :: dr_cld    ! Width of each rad_cld bin (m)
     103DOUBLE PRECISION, dimension(nbin_cld) :: vol_cld   ! Particle volume for each bin (m3)
     104REAL, SAVE                            :: sigma_ice ! Variance of the ice and CCN distributions
    105105!$OMP THREADPRIVATE(sigma_ice)
    106106
    107 !----------------------------------     
    108 ! JN : used in subtimestep
    109 REAL :: microtimestep! subdivision of ptimestep (s)
    110 REAL, dimension(ngrid,nlay) :: subpdtcloud ! Temperature variation due to latent heat
    111 REAL, dimension(ngrid,nlay,nq) :: zq0 ! local initial value of tracers
    112 INTEGER, dimension(ngrid,nlay) :: zimicro ! Subdivision of ptimestep
    113 INTEGER, dimension(ngrid,nlay) :: count_micro ! Number of microphys calls
    114 REAL, dimension(ngrid,nlay) :: zpotcond ! maximal condensable water (previous two)
    115 REAL, dimension(1) :: zqsat_tmp ! maximal condensable water (previous two)
    116 REAL :: spenttime ! time spent in while loop
    117 REAL :: zdq ! used to compute adaptative timestep
    118 LOGICAL :: ending_ts ! Condition to end while loop
    119 
    120 !----------------------------------     
     107!----------------------------------
     108! JN: used in subtimestep
     109REAL                              :: microtimestep ! Subdivision of ptimestep (s)
     110REAL,    dimension(ngrid,nlay)    :: subpdtcloud  ! Temperature variation due to latent heat
     111REAL,    dimension(ngrid,nlay,nq) :: zq0           ! Local initial value of tracers
     112INTEGER, dimension(ngrid,nlay)    :: zimicro      ! Subdivision of ptimestep
     113INTEGER, dimension(ngrid,nlay)    :: count_micro  ! Number of microphys calls
     114REAL,    dimension(ngrid,nlay)    :: zpotcond      ! Maximal condensable water (previous two)
     115REAL,    dimension(1)             :: zqsat_tmp     ! Maximal condensable water (previous two)
     116REAL                              :: spenttime     ! Time spent in while loop
     117REAL                              :: zdq           ! Used to compute adaptative timestep
     118LOGICAL                           :: ending_ts    ! Condition to end while loop
     119
     120!----------------------------------
    121121! TESTS
    122 INTEGER :: countcells
    123 LOGICAL, SAVE :: test_flag    ! flag for test/debuging outputs
     122INTEGER       :: countcells
     123LOGICAL, SAVE :: test_flag ! Flag for test/debuging outputs
    124124!$OMP THREADPRIVATE(test_flag)
    125 REAL, dimension(ngrid,nlay) :: satubf,satuaf, res_out
     125REAL, dimension(ngrid,nlay) :: satubf, satuaf, res_out
    126126
    127127!------------------------------------------------------------------
    … …  
    133133!=============================================================
    134134! rad_cld is the primary radius grid used for microphysics computation.
    135 ! The grid spacing is computed assuming a constant volume ratio
    136 ! between two consecutive bins; i.e. vrat_cld.
    137 ! vrat_cld is determined from the boundary values of the size grid:
    138 ! rmin_cld and rmax_cld.
     135! The grid spacing is computed assuming a constant volume ratio between two consecutive bins; i.e. vrat_cld.
     136! vrat_cld is determined from the boundary values of the size grid: rmin_cld and rmax_cld.
    139137! The rb_cld array contains the boundary values of each rad_cld bin.
    140138! dr_cld is the width of each rad_cld bin.
    141139
    142   ! Volume ratio between two adjacent bins
    143   ! vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
    144   ! vrat_cld = exp(vrat_cld)
    145   vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
    146   vrat_cld = exp(vrat_cld)
    147   write(*,*) "vrat_cld", vrat_cld
    148 
    149   rb_cld(1)  = rbmin_cld
    150   rad_cld(1) = rmin_cld
    151   vol_cld(1) = 4./3. * dble(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.
     140    ! Volume ratio between two adjacent bins
     141    ! vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
     142    ! vrat_cld = exp(vrat_cld)
     143    vrat_cld = log(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
     144    vrat_cld = exp(vrat_cld)
     145    write(*,*) "vrat_cld", vrat_cld
     146
     147    rb_cld(1)  = rbmin_cld
     148    rad_cld(1) = rmin_cld
     149    vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
     150
     151    do i=1,nbin_cld-1
     152        rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
     153        vol_cld(i+1) = vol_cld(i) * vrat_cld
     154    enddo
     155
     156    do i=1,nbin_cld
     157        rb_cld(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) * rad_cld(i)
     158        dr_cld(i)  = rb_cld(i+1) - rb_cld(i)
     159    enddo
     160    rb_cld(nbin_cld+1) = rbmax_cld
     161    dr_cld(nbin_cld)   = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
     162
     163    write(*,*) ' '
     164    write(*,*)'Microphysics: size bin information:'
     165    write(*,*)'i,rb_cld(i), rad_cld(i),dr_cld(i)'
     166    write(*,*)'-----------------------------------'
     167    do i=1,nbin_cld
     168        write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i), dr_cld(i)
     169    enddo
     170    write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
     171    write(*,*)'-----------------------------------'
     172
     173    ! we save that so that it is not computed at each timestep and gridpoint
     174    rb_cld(:) = log(rb_cld(:))
     175
     176    ! Contact parameter of water ice on dust ( m=cos(theta) )
     177    ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     178    ! mteta is initialized in conf_phys
     179    write(*,*) 'water_param contact parameter:', mteta
     180
     181    ! Volume of a water molecule (m3)
     182    vo1 = mh2o / dble(rho_ice)
     183    ! Variance of the ice and CCN distributions
     184    sigma_ice = sqrt(log(1.+nuice_sed))
     185
     186    write(*,*) 'Variance of ice & CCN distribs:', sigma_ice
     187    write(*,*) 'nuice for sedimentation       :', nuice_sed
     188    write(*,*) 'Volume of a water molecule    :', vo1
     189
     190    test_flag = .false.
     191    firstcall = .false.
    194192ENDIF
    195193
    … …  
    197195! 1. Initialisation
    198196!=============================================================
    199 
    200197res_out(:,:) = 0
    201198rice(:,:) = 1.e-8
    … …  
    235232! 2. Compute saturation
    236233!=============================================================
    237 
    238234dev2 = 1. / ( sqrt(2.) * sigma_ice )
    239235
    … …  
    248244! Main loop over the GCM's grid
    249245DO l=1,nlay
    250   DO ig=1,ngrid
    251     ! Subtimestep : here we go
    252     IF (cloud_adapt_ts) call adapt_imicro(ptimestep,zpotcond(ig,l), zimicro(ig,l))
    253     spenttime = 0.
    254     dMicetot= 0.
    255     ending_ts=.false.
    256     DO while (.not.ending_ts)
    257       call watersat(1,(/zt(ig,l)/),(/pplay(ig,l)/),zqsat_tmp)
    258       zqsat(ig,l)=zqsat_tmp(1)
    259       ! Get the partial pressure of water vapor and its saturation ratio
    260       ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
    261       satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
    262       microtimestep=ptimestep/real(zimicro(ig,l))
    263      
    264       ! Initialize tracers for scavenging + hdo computations (JN)
    265       zq0(ig,l,:)=zq(ig,l,:)
    266 
    267       ! Check if we are integrating over ptimestep
    268       if (spenttime+microtimestep.ge.ptimestep) then
    269         microtimestep=ptimestep-spenttime
    270         ! If so : last call !
    271         ending_ts=.true.
    272       endif! (spenttime+microtimestep.ge.ptimestep) then
     246    DO ig=1,ngrid
     247        ! Subtimestep: here we go
     248        IF (cloud_adapt_ts) call adapt_imicro(ptimestep,zpotcond(ig,l),zimicro(ig,l))
     249        spenttime = 0.
     250        dMicetot= 0.
     251        ending_ts=.false.
     252        DO while (.not.ending_ts)
     253            call watersat(1,(/zt(ig,l)/),(/pplay(ig,l)/),zqsat_tmp)
     254            zqsat(ig,l)=zqsat_tmp(1)
     255            ! Get the partial pressure of water vapor and its saturation ratio
     256            ph2o = zq(ig,l,igcm_h2o_vap) * (mmean(ig,l)/18.) * pplay(ig,l)
     257            satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
     258            microtimestep=ptimestep/real(zimicro(ig,l))
     259
     260            ! Initialize tracers for scavenging + hdo computations (JN)
     261            zq0(ig,l,:) = zq(ig,l,:)
     262
     263            ! Check if we are integrating over ptimestep
     264            if (spenttime+microtimestep >= ptimestep) then
     265                microtimestep = ptimestep-spenttime
     266                ! If so: last call !
     267                ending_ts = .true.
     268            endif! (spenttime+microtimestep.ge.ptimestep) then
    273269
    274270!=============================================================
    275271! 3. Nucleation
    276272!=============================================================
    277 
    278       IF ( satu .ge. 1. ) THEN ! if there is condensation
    279         call updaterccn(zq(ig,l,igcm_dust_mass),zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
    280 
    281         ! Expand the dust moments into a binned distribution
    282         Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)   + 1.e-30
    283         No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
    284         Rn = rdust(ig,l)
    285         Rn = -log(Rn)
    286         Rm = Rn - 3. * sigma_ice*sigma_ice 
    287         n_derf = derf( (rb_cld(1)+Rn) *dev2)
    288         m_derf = derf( (rb_cld(1)+Rm) *dev2)
    289         do i = 1, nbin_cld
    290           n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
    291           m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
    292           n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
    293           m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
    294           n_aer(i) = n_aer(i) + 0.5 * No * n_derf
    295           m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
    296         enddo
    297      
    298         ! Get the rates of nucleation
    299         call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
    300      
    301         dN = 0.
    302         dM = 0.
    303         do i = 1, nbin_cld
    304           dN = dN + n_aer(i)*(exp(-rate(i)*microtimestep)-1.)
    305           dM = dM + m_aer(i)*(exp(-rate(i)*microtimestep)-1.)
    306         enddo
    307 
    308         ! Update Dust particles
    309         zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    310         zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    311         ! Update CCNs
    312         zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    313         zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
    314      
    315       ENDIF ! of is satu >1
     273            IF ( satu >= 1. ) THEN ! if there is condensation
     274                call updaterccn(zq(ig,l,igcm_dust_mass),zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
     275
     276                ! Expand the dust moments into a binned distribution
     277                Mo = zq(ig,l,igcm_dust_mass)  *tauscaling(ig) + 1.e-30
     278                No = zq(ig,l,igcm_dust_number)*tauscaling(ig) + 1.e-30
     279                Rn = rdust(ig,l)
     280                Rn = -log(Rn)
     281                Rm = Rn - 3. * sigma_ice*sigma_ice
     282                n_derf = derf( (rb_cld(1)+Rn) *dev2)
     283                m_derf = derf( (rb_cld(1)+Rm) *dev2)
     284                do i = 1, nbin_cld
     285                    n_aer(i) = -0.5 * No * n_derf ! this ith previously computed
     286                    m_aer(i) = -0.5 * Mo * m_derf ! this ith previously computed
     287                    n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
     288                    m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
     289                    n_aer(i) = n_aer(i) + 0.5 * No * n_derf
     290                    m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
     291                enddo
     292
     293                ! Get the rates of nucleation
     294                call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
     295
     296                dN = 0.
     297                dM = 0.
     298                do i = 1, nbin_cld
     299                    dN = dN + n_aer(i)*(exp(-rate(i)*microtimestep)-1.)
     300                    dM = dM + m_aer(i)*(exp(-rate(i)*microtimestep)-1.)
     301                enddo
     302
     303                ! Update Dust particles
     304                zq(ig,l,igcm_dust_mass)   = zq(ig,l,igcm_dust_mass)   + dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
     305                zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
     306                ! Update CCNs
     307                zq(ig,l,igcm_ccn_mass)   = zq(ig,l,igcm_ccn_mass)   - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
     308                zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
     309            ENDIF ! of is satu >1
    316310
    317311!=============================================================
    … …  
    321315! Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
    322316! to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
    323       IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.) THEN ! we trigger crystal growth
    324         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))
    325 
    326         No  = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
    327 
    328         ! saturation at equilibrium
    329         ! rice should not be too small, otherwise seq value is not valid
    330         seq  = exp(2.*sig_h2o(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)*max(rice(ig,l),1.e-7)))
    331      
    332         ! get resistance growth
    333         call growthrate(zt(ig,l),pplay(ig,l),real(ph2o/satu),rice(ig,l),res,Dv)
    334 
    335         res_out(ig,l) = res
    336 
    337         ! implicit scheme of mass growth
    338         ! cste here must be computed at each step
    339         cste = 4*pi*rho_ice*microtimestep
    340 
    341         dMice = (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))/(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.)
    342 
    343         ! With the above scheme, dMice cannot be bigger than vapor, but can be bigger than all available ice.
    344         dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
    345         dMice = min(dMice,zq(ig,l,igcm_h2o_vap))
    346 
    347         zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)+dMice
    348         zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)-dMice
    349 
    350         countcells = countcells + 1
    351      
    352         ! latent heat release
    353         lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3
    354         subpdtcloud(ig,l)= dMice*lw/cpp/microtimestep
    355      
    356         ! DIFF: trend is enforce in a range, stabilize the scheme ?
    357         if (subpdtcloud(ig,l)*microtimestep.gt.5.0) then
    358           subpdtcloud(ig,l)=5./microtimestep
    359         endif! (subpdtcloud(ig,l)*microtimestep.gt.5.0) then
    360         if (subpdtcloud(ig,l)*microtimestep.lt.-5.0) then
    361             subpdtcloud(ig,l)=-5./microtimestep
    362         endif! (subpdtcloud(ig,l)*microtimestep.gt.5.0) then
    363 
    364         ! Special case of the isotope of water HDO   
    365         if (hdo) then
    366           ! condensation
    367           if (dMice.gt.0.0) then
    368             ! do we use fractionation?
    369             if (hdofrac) then
    370               ! Calculation of the HDO vapor coefficient
    371               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) )
    372               ! Calculation of the fractionnation coefficient at equilibrium
    373               alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
    374               ! Calculation of the 'real' fractionnation coefficient
    375               alpha_c(ig,l) = (alpha(ig,l)*satu)/( (alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.)
    376             else
    377               alpha_c(ig,l) = 1.d0
    378             endif
    379             if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then
    380               dMice_hdo = dMice*alpha_c(ig,l)*( zq0(ig,l,igcm_hdo_vap)/zq0(ig,l,igcm_h2o_vap) )
    381             else
    382               dMice_hdo=0.
    383             endif
    384           !! sublimation
    385           else
    386             if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then
    387               dMice_hdo=dMice*( zq0(ig,l,igcm_hdo_ice)/zq0(ig,l,igcm_h2o_ice) )
    388             else
    389               dMice_hdo=0.
    390             endif
    391           endif !if (dMice.gt.0.0)
    392 
    393           dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice))
    394           dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap))
    395 
    396           zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice)+dMice_hdo
    397           zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap)-dMice_hdo
    398 
    399         endif ! if (hdo)       
     317            IF ( zq(ig,l,igcm_ccn_number)*tauscaling(ig) >= 1.) THEN ! we trigger crystal growth
     318                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))
     319
     320                No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
     321
     322                ! saturation at equilibrium
     323                ! rice should not be too small, otherwise seq value is not valid
     324                seq  = exp(2.*sig_h2o(zt(ig,l))*mh2o / (rho_ice*rgp*zt(ig,l)*max(rice(ig,l),1.e-7)))
     325
     326                ! get resistance growth
     327                call growthrate(zt(ig,l),pplay(ig,l),real(ph2o/satu),rice(ig,l),res,Dv)
     328
     329                res_out(ig,l) = res
     330
     331                ! implicit scheme of mass growth
     332                ! cste here must be computed at each step
     333                cste = 4*pi*rho_ice*microtimestep
     334
     335                dMice = (zq(ig,l,igcm_h2o_vap)-seq*zqsat(ig,l))/(res*zqsat(ig,l)/(cste*No*rice(ig,l)) + 1.)
     336
     337                ! With the above scheme, dMice cannot be bigger than vapor, but can be bigger than all available ice.
     338                dMice = max(dMice,-zq(ig,l,igcm_h2o_ice))
     339                dMice = min(dMice,zq(ig,l,igcm_h2o_vap))
     340
     341                zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) + dMice
     342                zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) - dMice
     343
     344                countcells = countcells + 1
     345
     346                ! latent heat release
     347                lw = (2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)*(zt(ig,l)-To))*1.e+3
     348                subpdtcloud(ig,l) = dMice*lw/cpp/microtimestep
     349
     350                ! DIFF: trend is enforce in a range, stabilize the scheme?
     351                if (subpdtcloud(ig,l)*microtimestep > 5.0) then
     352                    subpdtcloud(ig,l) = 5./microtimestep
     353                endif
     354                if (subpdtcloud(ig,l)*microtimestep < -5.0) then
     355                    subpdtcloud(ig,l) = -5./microtimestep
     356                endif
     357
     358                ! Special case of the isotope of water HDO
     359                if (hdo) then
     360                    ! condensation
     361                    if (dMice > 0.) then
     362                        ! do we use fractionation?
     363                        if (hdofrac) then
     364                            ! Calculation of the HDO vapor coefficient
     365                            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) )
     366                            ! Calculation of the fractionnation coefficient at equilibrium
     367                            alpha(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
     368                            ! Calculation of the 'real' fractionnation coefficient
     369                            alpha_c(ig,l) = (alpha(ig,l)*satu) / ((alpha(ig,l)*(Dv/Dv_hdo)*(satu-1.)) + 1.)
     370                        else
     371                            alpha_c(ig,l) = 1.d0
     372                        endif
     373                        if (zq0(ig,l,igcm_h2o_vap) > qperemin) then
     374                            dMice_hdo = dMice*alpha_c(ig,l)*( zq0(ig,l,igcm_hdo_vap)/zq0(ig,l,igcm_h2o_vap) )
     375                        else
     376                            dMice_hdo = 0.
     377                        endif
     378                    ! sublimation
     379                    else
     380                        if (zq0(ig,l,igcm_h2o_ice) > qperemin) then
     381                            dMice_hdo = dMice*( zq0(ig,l,igcm_hdo_ice)/zq0(ig,l,igcm_h2o_ice) )
     382                        else
     383                            dMice_hdo = 0.
     384                        endif
     385                    endif ! if (dMice > 0.)
     386
     387                    dMice_hdo = max(dMice_hdo,-zq(ig,l,igcm_hdo_ice))
     388                    dMice_hdo = min(dMice_hdo,zq(ig,l,igcm_hdo_vap))
     389
     390                    zq(ig,l,igcm_hdo_ice) = zq(ig,l,igcm_hdo_ice) + dMice_hdo
     391                    zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) - dMice_hdo
     392
     393                endif ! if (hdo)
    400394
    401395!=============================================================
    402396! 5. Dust cores released, tendancies, latent heat, etc ...
    403397!=============================================================
    404 
    405         ! If all the ice particles sublimate, all the condensation
    406         ! nuclei are released:
    407         if (zq(ig,l,igcm_h2o_ice).le.1.e-28) then
    408           ! Water
    409           zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) + zq(ig,l,igcm_h2o_ice)
    410           zq(ig,l,igcm_h2o_ice) = 0.
    411           if (hdo) then
    412             zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) + zq(ig,l,igcm_hdo_ice)
    413             zq(ig,l,igcm_hdo_ice) = 0.
    414           endif
    415           ! Dust particles
    416           zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccn_mass)
    417           zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccn_number)
    418           ! CCNs
    419           zq(ig,l,igcm_ccn_mass) = 0.
    420           zq(ig,l,igcm_ccn_number) = 0.
    421         endif
    422      
    423       ELSE
    424         ! Initialization of dMice when it's not computed
    425         dMice=0 ! no condensation/sublimation to account for
    426         subpdtcloud(ig,l)=0 ! no condensation/sublimation to account for
    427       ENDIF !of if Nccn>1
    428      
    429       ! No more getting tendency : we increment tracers & temp directly
    430 
    431       ! If not activice, the tendency from latent heat release is set to zero
    432       IF (.not.activice) subpdtcloud(ig,l)=0.
    433 
    434       ! Temperature change as a feedback from latent heat release
    435       zt(ig,l) = zt(ig,l)+subpdtcloud(ig,l)*microtimestep
    436 
    437       ! Prevent negative tracers ! JN
    438       WHERE (zq(ig,l,:) < 1.e-30) zq(ig,l,:) = 1.e-30
    439 
    440       IF (cloud_adapt_ts) then
    441         ! Estimation of how much is actually condensing/subliming
    442          dMicetot=dMicetot+abs(dMice) ! total accumulated ice formation since the beginning
    443         IF (spenttime.ne.0) then
    444           zdq=(dMicetot/spenttime)!*(ptimestep-spenttime)
    445         ELSE
    446           ! Initialization for spenttime=0
    447           zdq=zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep)
    448         ENDIF
    449         zdq=abs(zdq)
    450         call adapt_imicro(ptimestep,zdq,zimicro(ig,l))
    451       ENDIF! (cloud_adapt_ts) then
    452       ! Increment time spent in here
    453       spenttime=spenttime+microtimestep
    454       count_micro(ig,l)=count_micro(ig,l)+1
    455     ENDDO ! while (.not. ending_ts)
    456   ENDDO ! of ig loop
     398                ! If all the ice particles sublimate, all the condensation
     399                ! nuclei are released:
     400                if (zq(ig,l,igcm_h2o_ice) <= 1.e-28) then
     401                    ! Water
     402                    zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap) + zq(ig,l,igcm_h2o_ice)
     403                    zq(ig,l,igcm_h2o_ice) = 0.
     404                    if (hdo) then
     405                        zq(ig,l,igcm_hdo_vap) = zq(ig,l,igcm_hdo_vap) + zq(ig,l,igcm_hdo_ice)
     406                        zq(ig,l,igcm_hdo_ice) = 0.
     407                    endif
     408                    ! Dust particles
     409                    zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccn_mass)
     410                    zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccn_number)
     411                    ! CCNs
     412                    zq(ig,l,igcm_ccn_mass) = 0.
     413                    zq(ig,l,igcm_ccn_number) = 0.
     414                endif
     415
     416            ELSE
     417                ! Initialization of dMice when it's not computed
     418                dMice = 0 ! no condensation/sublimation to account for
     419                subpdtcloud(ig,l) = 0 ! no condensation/sublimation to account for
     420            ENDIF ! of if Nccn>1
     421
     422            ! No more getting tendency: we increment tracers & temp directly
     423
     424            ! If not activice, the tendency from latent heat release is set to zero
     425            IF (.not.activice) subpdtcloud(ig,l)=0.
     426
     427            ! Temperature change as a feedback from latent heat release
     428            zt(ig,l) = zt(ig,l) + subpdtcloud(ig,l)*microtimestep
     429
     430            ! Prevent negative tracers ! JN
     431            WHERE (zq(ig,l,:) < 1.e-30) zq(ig,l,:) = 1.e-30
     432
     433            IF (cloud_adapt_ts) then
     434                ! Estimation of how much is actually condensing/subliming
     435                dMicetot=dMicetot+abs(dMice) ! total accumulated ice formation since the beginning
     436                IF (spenttime /= 0) then
     437                    zdq = (dMicetot/spenttime)!*(ptimestep-spenttime)
     438                ELSE
     439                    ! Initialization for spenttime=0
     440                    zdq = zpotcond(ig,l)*((ptimestep-spenttime)/ptimestep)
     441                ENDIF
     442                    zdq = abs(zdq)
     443                    call adapt_imicro(ptimestep,zdq,zimicro(ig,l))
     444            ENDIF ! (cloud_adapt_ts) then
     445            ! Increment time spent in here
     446            spenttime = spenttime + microtimestep
     447            count_micro(ig,l) = count_micro(ig,l) + 1
     448        ENDDO ! while (.not. ending_ts)
     449    ENDDO ! of ig loop
    457450ENDDO ! of nlayer loop
    458451
    … …  
    461454call write_output("count_micro","count_micro after microphysics","integer",count_micro(:,:))
    462455
    463 !!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS 
     456!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
    464457IF (test_flag) then
    465   DO l=1,nlay
    466     DO ig=1,ngrid
    467       satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
    468       satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
     458    DO l=1,nlay
     459        DO ig=1,ngrid
     460            satubf(ig,l) = zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l)
     461            satuaf(ig,l) = zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
     462        ENDDO
    469463    ENDDO
    470   ENDDO
    471   print*, 'count is ',countcells, ' i.e. ', countcells*100/(nlay*ngrid), '% for microphys computation'
     464    write(*,*) 'count is ',countcells, ' i.e. ', countcells*100/(nlay*ngrid), '% for microphys computation'
    472465ENDIF ! endif test_flag
    473466
    … …  
    483476! efficiency.
    484477
    485 real,intent(in) :: ptimestep ! total duration of physics (sec)
    486 real,intent(in) :: potcond ! condensible vapor / sublimable ice(kg/kg)
    487 real :: alpha, beta ! Coefficients for power law
    488 real :: defstep,coef ! Default ptimestep of 7.5 mins (iphysiq=5)
    489 integer,intent(out) :: zimicro ! number of ptimestep division
    490 
    491 ! Default ptimestep : defstep (7.5 mins)
    492 defstep=88775.*5./960.
    493 coef=ptimestep/defstep
     478real,    intent(in)  :: ptimestep ! Total duration of physics (sec)
     479real,    intent(in)  :: potcond   ! Condensible vapor / sublimable ice(kg/kg)
     480integer, intent(out) :: zimicro   ! Number of ptimestep division
     481real :: alpha, beta   ! Coefficients for power law
     482real :: defstep, coef ! Default ptimestep of 7.5 mins (iphysiq=5)
     483
     484! Default ptimestep: defstep (7.5 mins)
     485defstep = 88775.*5./960.
     486coef = ptimestep/defstep
    494487! Conservative coefficients :
    495 ! alpha=1.81846504e+11
    496 ! beta=1.54550140e+00
    497 alpha=1.88282793e+05 ! Latest values for high obliquity
    498 beta=4.57764370e-01  ! Latest values for high obliquity
    499 !alpha=1.72198978e+10 ! Present day Mars
    500 !beta=1.88473210e+00 
    501 zimicro=ceiling(coef*min(max(alpha*abs(potcond)**beta,5.),7000.))
    502 !zimicro=2*zimicro ! Prediction times two, allow to complete year at high obliquity
     488!alpha = 1.81846504e+11
     489!beta = 1.54550140e+00
     490alpha = 1.88282793e+05 ! Latest values for high obliquity
     491beta = 4.57764370e-01  ! Latest values for high obliquity
     492!alpha = 1.72198978e+10 ! Present day Mars
     493!beta = 1.88473210e+00 
     494zimicro = ceiling(coef*min(max(alpha*abs(potcond)**beta,5.),7000.))
     495!zimicro = 2*zimicro ! Prediction times two, allow to complete year at high obliquity
    503496
    504497END SUBROUTINE adapt_imicro
Note: See TracChangeset for help on using the changeset viewer.