LCOV - code coverage report
Current view: | top level - physics/cam - nucleate_ice.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 213 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 6 | 0.0 %

              Line data    Source code
           1             : module nucleate_ice
           2             :
           3             : !-------------------------------------------------------------------------------
           4             : ! Purpose:
           5             : !  A parameterization of ice nucleation.
           6             : !
           7             : !  *** This module is intended to be a "portable" code layer.  Ideally it should
           8             : !  *** not contain any use association of modules that belong to the model framework.
           9             : !
          10             : !
          11             : ! Method:
          12             : !  The current method is based on Liu & Penner (2005) & Liu et al. (2007)
          13             : !  It related the ice nucleation with the aerosol number, temperature and the
          14             : !  updraft velocity. It includes homogeneous freezing of sulfate & immersion
          15             : !  freezing on mineral dust (soot disabled) in cirrus clouds, and
          16             : !  Meyers et al. (1992) deposition nucleation in mixed-phase clouds
          17             : !
          18             : !  The effect of preexisting ice crystals on ice nucleation in cirrus clouds is included,
          19             : !  and also consider the sub-grid variability of temperature in cirrus clouds,
          20             : !  following X. Shi et al. ACP (2014).
          21             : !
          22             : !  Ice nucleation in mixed-phase clouds now uses classical nucleation theory (CNT),
          23             : !  follows Y. Wang et al. ACP (2014), Hoose et al. (2010).
          24             : !
          25             : ! Authors:
          26             : !  Xiaohong Liu, 01/2005, modifications by A. Gettelman 2009-2010
          27             : !  Xiangjun Shi & Xiaohong Liu, 01/2014.
          28             : !
          29             : !  With help from C. C. Chen and B. Eaton (2014)
          30             : !-------------------------------------------------------------------------------
          31             :
          32             : use wv_saturation,  only: svp_water, svp_ice
          33             :
          34             : implicit none
          35             : private
          36             : save
          37             :
          38             : integer, parameter :: r8 = selected_real_kind(12)
          39             :
          40             : public :: nucleati_init, nucleati
          41             :
          42             : logical  :: use_preexisting_ice
          43             : logical  :: use_hetfrz_classnuc
          44             : logical  :: use_incloud_nuc
          45             : integer  :: iulog
          46             : real(r8) :: pi
          47             : real(r8) :: mincld
          48             :
          49             : real(r8), parameter :: Shet   = 1.3_r8     ! het freezing threshold
          50             : real(r8), parameter :: rhoice = 0.5e3_r8   ! kg/m3, Wpice is not sensitive to rhoice
          51             : real(r8), parameter :: minweff= 0.001_r8   ! m/s
          52             : real(r8), parameter :: gamma1=1.0_r8
          53             : real(r8), parameter :: gamma2=1.0_r8
          54             : real(r8), parameter :: gamma3=2.0_r8
          55             : real(r8), parameter :: gamma4=6.0_r8
          56             :
          57             : real(r8) :: ci
          58             :
          59             : !===============================================================================
          60             : contains
          61             : !===============================================================================
          62             :
          63           0 : subroutine nucleati_init( &
          64             :    use_preexisting_ice_in, use_hetfrz_classnuc_in, use_incloud_nuc_in, iulog_in, pi_in, &
          65             :    mincld_in)
          66             :
          67             :    logical,  intent(in) :: use_preexisting_ice_in
          68             :    logical,  intent(in) :: use_hetfrz_classnuc_in
          69             :    logical,  intent(in) :: use_incloud_nuc_in
          70             :    integer,  intent(in) :: iulog_in
          71             :    real(r8), intent(in) :: pi_in
          72             :    real(r8), intent(in) :: mincld_in
          73             :
          74           0 :    use_preexisting_ice = use_preexisting_ice_in
          75           0 :    use_hetfrz_classnuc = use_hetfrz_classnuc_in
          76           0 :    use_incloud_nuc     = use_incloud_nuc_in
          77           0 :    iulog               = iulog_in
          78           0 :    pi                  = pi_in
          79           0 :    mincld              = mincld_in
          80             :
          81           0 :    ci = rhoice*pi/6._r8
          82             :
          83           0 : end subroutine nucleati_init
          84             :
          85             : !===============================================================================
          86             :
          87           0 : subroutine nucleati(  &
          88             :    wbar, tair, pmid, relhum, cldn,      &
          89             :    qc, qi, ni_in, rhoair,               &
          90             :    so4_num, dst_num, soot_num, subgrid, &
          91             :    nuci, onihf, oniimm, onidep, onimey, &
          92             :    wpice, weff, fhom, regm, &
          93             :    oso4_num, odst_num, osoot_num, &
          94             :    call_frm_zm_in, add_preexisting_ice_in)
          95             :
          96             :    ! Input Arguments
          97             :    real(r8), intent(in) :: wbar        ! grid cell mean vertical velocity (m/s)
          98             :    real(r8), intent(in) :: tair        ! temperature (K)
          99             :    real(r8), intent(in) :: pmid        ! pressure at layer midpoints (pa)
         100             :    real(r8), intent(in) :: relhum      ! relative humidity with respective to liquid
         101             :    real(r8), intent(in) :: cldn        ! new value of cloud fraction    (fraction)
         102             :    real(r8), intent(in) :: qc          ! liquid water mixing ratio (kg/kg)
         103             :    real(r8), intent(in) :: qi          ! grid-mean preexisting cloud ice mass mixing ratio (kg/kg)
         104             :    real(r8), intent(in) :: ni_in       ! grid-mean preexisting cloud ice number conc (#/kg)
         105             :    real(r8), intent(in) :: rhoair      ! air density (kg/m3)
         106             :    real(r8), intent(in) :: so4_num     ! so4 aerosol number (#/cm^3)
         107             :    real(r8), intent(in) :: dst_num     ! total dust aerosol number (#/cm^3)
         108             :    real(r8), intent(in) :: soot_num    ! soot (hydrophilic) aerosol number (#/cm^3)
         109             :    real(r8), intent(in) :: subgrid     ! subgrid saturation scaling factor
         110             :
         111             :    ! Output Arguments
         112             :    real(r8), intent(out) :: nuci       ! ice number nucleated (#/kg)
         113             :    real(r8), intent(out) :: onihf      ! nucleated number from homogeneous freezing of so4
         114             :    real(r8), intent(out) :: oniimm     ! nucleated number from immersion freezing
         115             :    real(r8), intent(out) :: onidep     ! nucleated number from deposition nucleation
         116             :    real(r8), intent(out) :: onimey     ! nucleated number from deposition nucleation  (meyers: mixed phase)
         117             :    real(r8), intent(out) :: wpice      ! diagnosed Vertical velocity Reduction caused by preexisting ice (m/s), at Shom
         118             :    real(r8), intent(out) :: weff       ! effective Vertical velocity for ice nucleation (m/s); weff=wbar-wpice
         119             :    real(r8), intent(out) :: fhom       ! how much fraction of cloud can reach Shom
         120             :    real(r8), intent(out) :: regm       ! nucleation regime indiator
         121             :    real(r8), intent(out) :: oso4_num   ! so4 aerosol number (#/cm^3)
         122             :    real(r8), intent(out) :: odst_num   ! total dust aerosol number (#/cm^3)
         123             :    real(r8), intent(out) :: osoot_num  ! soot (hydrophilic) aerosol number (#/cm^3)
         124             :
         125             :    ! Optional Arguments
         126             :    logical,  intent(in), optional :: call_frm_zm_in ! true if called from ZM convection scheme
         127             :    logical,  intent(in), optional :: add_preexisting_ice_in ! only false if called with pumas_v1.21+
         128             :
         129             :    ! Local workspace
         130             :    real(r8) :: nihf                      ! nucleated number from homogeneous freezing of so4
         131             :    real(r8) :: niimm                     ! nucleated number from immersion freezing
         132             :    real(r8) :: nidep                     ! nucleated number from deposition nucleation
         133             :    real(r8) :: nimey                     ! nucleated number from deposition nucleation (meyers)
         134             :    real(r8) :: n1, ni                    ! nucleated number
         135             :    real(r8) :: tc, A, B                  ! work variable
         136             :    real(r8) :: esl, esi, deles           ! work variable
         137             :    real(r8) :: wbar1, wbar2
         138             :
         139             :    ! used in SUBROUTINE Vpreice
         140             :    real(r8) :: Ni_preice        ! cloud ice number conc (1/m3)
         141             :    real(r8) :: lami,Ri_preice   ! mean cloud ice radius (m)
         142             :    real(r8) :: Shom             ! initial ice saturation ratio; if <1, use hom threshold Si
         143             :    real(r8) :: detaT,RHimean    ! temperature standard deviation, mean cloudy RHi
         144             :    real(r8) :: wpicehet   ! diagnosed Vertical velocity Reduction caused by preexisting ice (m/s), at shet
         145             :
         146             :    real(r8) :: weffhet    ! effective Vertical velocity for ice nucleation (m/s)  weff=wbar-wpicehet
         147             :
         148             :    logical  :: call_frm_zm, add_preexisting_ice
         149             :    !-------------------------------------------------------------------------------
         150             :
         151           0 :    nuci = 0._r8
         152             :
         153           0 :    RHimean = relhum*svp_water(tair)/svp_ice(tair)*subgrid
         154             :
         155             :    ! temp variables that depend on use_preexisting_ice
         156           0 :    wbar1 = wbar
         157           0 :    wbar2 = wbar
         158             :
         159             :    ! If not using prexisting ice, the homogeneous freezing happens in the
         160             :    ! entire gridbox.
         161           0 :    fhom = 1._r8
         162             :
         163           0 :    if (present(call_frm_zm_in)) then
         164           0 :      call_frm_zm = call_frm_zm_in
         165             :    else
         166             :      call_frm_zm = .false.
         167             :    end if
         168             :
         169           0 :    if (present(add_preexisting_ice_in)) then
         170           0 :      add_preexisting_ice = add_preexisting_ice_in
         171             :    else
         172             :      add_preexisting_ice = .true.
         173             :    end if
         174             :
         175           0 :    if (use_preexisting_ice .and. (.not. call_frm_zm)) then
         176             :
         177           0 :       Ni_preice = ni_in*rhoair                    ! (convert from #/kg -> #/m3)
         178           0 :       Ni_preice = Ni_preice / max(mincld,cldn)   ! in-cloud ice number density
         179             :
         180           0 :       if (Ni_preice > 10.0_r8 .and. qi > 1.e-10_r8) then    ! > 0.01/L = 10/m3
         181           0 :          Shom = -1.5_r8   ! if Shom<1 , Shom will be recalculated in SUBROUTINE Vpreice, according to Ren & McKenzie, 2005
         182           0 :          lami = (gamma4*ci*ni_in/qi)**(1._r8/3._r8)
         183           0 :          Ri_preice = 0.5_r8/lami                  ! radius
         184           0 :          Ri_preice = max(Ri_preice, 1e-8_r8)       ! >0.01micron
         185           0 :          call Vpreice(pmid, tair, Ri_preice, Ni_preice, Shom, wpice)
         186           0 :          call Vpreice(pmid, tair, Ri_preice, Ni_preice, Shet, wpicehet)
         187             :       else
         188           0 :          wpice    = 0.0_r8
         189           0 :          wpicehet = 0.0_r8
         190             :       endif
         191             :
         192           0 :       weff     = max(wbar-wpice, minweff)
         193           0 :       wpice    = min(wpice, wbar)
         194           0 :       weffhet  = max(wbar-wpicehet,minweff)
         195           0 :       wpicehet = min(wpicehet, wbar)
         196             :
         197           0 :       wbar1 = weff
         198           0 :       wbar2 = weffhet
         199             :
         200           0 :       detaT   = wbar/0.23_r8
         201           0 :       if (use_incloud_nuc) then
         202           0 :         call frachom(tair, 1._r8, detaT, fhom)
         203             :       else
         204           0 :         call frachom(tair, RHimean, detaT, fhom)
         205             :       end if
         206             :    end if
         207             :
         208           0 :    ni = 0._r8
         209           0 :    tc = tair - 273.15_r8
         210             :
         211             :    ! initialize
         212           0 :    niimm = 0._r8
         213           0 :    nidep = 0._r8
         214           0 :    nihf  = 0._r8
         215           0 :    deles = 0._r8
         216           0 :    esi   = 0._r8
         217           0 :    regm  = 0._r8
         218             :
         219           0 :    oso4_num  = 0._r8
         220           0 :    odst_num  = 0._r8
         221           0 :    osoot_num = 0._r8
         222             :
         223           0 :    if ((so4_num >= 1.0e-10_r8 .or. (soot_num+dst_num) >= 1.0e-10_r8) .and. cldn > 0._r8) then
         224             :
         225           0 :       if (RHimean.ge.1.2_r8) then
         226             :
         227           0 :          if ( ((tc.le.0.0_r8).and.(tc.ge.-37.0_r8).and.(qc.lt.1.e-12_r8)).or.(tc.le.-37.0_r8)) then
         228             :
         229           0 :             A = -1.4938_r8 * log(soot_num+dst_num) + 12.884_r8
         230           0 :             B = -10.41_r8  * log(soot_num+dst_num) - 67.69_r8
         231           0 :             regm = A * log(wbar1) + B
         232             :
         233             :             ! heterogeneous nucleation only
         234           0 :             if (tc .gt. regm .or. so4_num < 1.0e-10_r8) then
         235             :
         236           0 :                if(tc.lt.-40._r8 .and. wbar1.gt.1._r8 .and. so4_num >= 1.0e-10_r8) then ! exclude T<-40 & W>1m/s from hetero. nucleation
         237             :
         238           0 :                   call hf(tc,wbar1,relhum*subgrid,so4_num,nihf)
         239           0 :                   niimm=0._r8
         240           0 :                   nidep=0._r8
         241             :
         242             :                   ! If some homogeneous nucleation happened, assume all of the that heterogeneous
         243             :                   ! and coarse mode sulfate particles nucleated.
         244           0 :                   if (nihf.gt.1e-3_r8) then ! hom occur,  add preexisting ice
         245           0 :                      niimm     = dst_num + soot_num       ! assuming dst_num freeze firstly
         246           0 :                      odst_num  = dst_num
         247           0 :                      osoot_num = soot_num
         248             :
         249           0 :                      oso4_num  = nihf
         250             :                   endif
         251             :
         252           0 :                   nihf      = nihf * fhom
         253           0 :                   oso4_num  = oso4_num * fhom
         254             :
         255           0 :                   n1        = nihf + niimm
         256             :                else
         257             :
         258           0 :                   call hetero(tc,wbar2,soot_num+dst_num,niimm,nidep)
         259             :
         260           0 :                   nihf = 0._r8
         261           0 :                   n1   = niimm + nidep
         262             :
         263           0 :                   osoot_num = soot_num * (niimm + nidep) / (soot_num + dst_num)
         264           0 :                   odst_num  = dst_num  * (niimm + nidep) / (soot_num + dst_num)
         265             :                endif
         266             :
         267             :             ! homogeneous nucleation only
         268           0 :             else if (tc.lt.regm-5._r8 .or. (soot_num+dst_num) < 1.0e-10_r8) then
         269             :
         270           0 :                call hf(tc,wbar1,relhum*subgrid,so4_num,nihf)
         271           0 :                niimm=0._r8
         272           0 :                nidep=0._r8
         273             :
         274             :                ! If some homogeneous nucleation happened, assume all of the that
         275             :                ! heterogeneous and coarse mode sulfate particles nucleated.
         276           0 :                if (nihf.gt.1e-3_r8) then !  hom occur,  add preexisting ice
         277           0 :                   niimm     = dst_num + soot_num       ! assuming dst_num freeze firstly
         278           0 :                   odst_num  = dst_num
         279           0 :                   osoot_num = soot_num
         280             :
         281           0 :                   oso4_num  = nihf
         282             :                endif
         283             :
         284           0 :                nihf      = nihf * fhom
         285           0 :                oso4_num  = oso4_num * fhom
         286             :
         287           0 :                n1        = nihf + niimm
         288             :
         289             :             ! transition between homogeneous and heterogeneous: interpolate in-between
         290             :             else
         291             :
         292           0 :                if (tc.lt.-40._r8 .and. wbar1.gt.1._r8) then ! exclude T<-40 & W>1m/s from hetero. nucleation
         293             :
         294           0 :                   call hf(tc, wbar1, relhum*subgrid, so4_num, nihf)
         295           0 :                   niimm = 0._r8
         296           0 :                   nidep = 0._r8
         297             :
         298             :                   ! If some homogeneous nucleation happened, assume all of the
         299             :                   ! that heterogeneous and coarse mode sulfate particles nucleated.
         300           0 :                   if (nihf.gt.1e-3_r8) then ! hom occur,  add preexisting ice
         301           0 :                      niimm     = dst_num + soot_num       ! assuming dst_num freeze firstly
         302           0 :                      odst_num  = dst_num
         303           0 :                      osoot_num = soot_num
         304             :
         305           0 :                      oso4_num  = nihf
         306             :                   endif
         307             :
         308           0 :                   nihf      = nihf * fhom
         309           0 :                   oso4_num  = oso4_num * fhom
         310             :
         311           0 :                   n1        = nihf + niimm
         312             :
         313             :                else
         314             :
         315           0 :                   call hf(regm-5._r8,wbar1,relhum*subgrid,so4_num,nihf)
         316           0 :                   call hetero(regm,wbar2,soot_num+dst_num,niimm,nidep)
         317             :
         318             :                   ! If some homogeneous nucleation happened, assume all of the
         319             :                   ! heterogeneous particles nucleated and add in a fraction of
         320             :                   ! the homogeneous freezing.
         321           0 :                   if (nihf.gt.1e-3_r8) then ! hom occur,  add preexisting ice
         322           0 :                      oso4_num  = nihf
         323             :                   endif
         324             :
         325           0 :                   osoot_num = soot_num * (niimm + nidep) / (soot_num + dst_num)
         326           0 :                   odst_num  = dst_num  * (niimm + nidep) / (soot_num + dst_num)
         327             :
         328           0 :                   nihf      = nihf      * fhom * ((regm - tc) / 5._r8)**2
         329           0 :                   oso4_num  = oso4_num  * fhom * ((regm - tc) / 5._r8)**2
         330             :
         331           0 :                   n1 = niimm + nidep + nihf
         332             :
         333             :                end if
         334             :             end if
         335             :
         336             :             ! Scale the rates for in-cloud number, since this is what
         337             :             ! MG is expecting to find.
         338           0 :             ni = n1
         339             :
         340             :             ! If using prexsiting ice, and allowed to add, then add it to the total.
         341           0 :             if (use_preexisting_ice) then
         342           0 :                if (add_preexisting_ice .and. (.not. call_frm_zm)) then
         343           0 :                   ni = ni + Ni_preice * 1e-6_r8
         344             :                end if
         345             :             end if
         346             :          end if
         347             :       end if
         348             :    end if
         349             :
         350             :    ! deposition/condensation nucleation in mixed clouds (-37<T<0C) (Meyers, 1992)
         351           0 :    if(tc.lt.0._r8 .and. tc.gt.-37._r8 .and. qc.gt.1.e-12_r8) then
         352           0 :       esl = svp_water(tair)     ! over water in mixed clouds
         353           0 :       esi = svp_ice(tair)     ! over ice
         354           0 :       deles = (esl - esi)
         355           0 :       nimey=1.e-3_r8*exp(12.96_r8*deles/esi - 0.639_r8)
         356             :    else
         357             :       nimey=0._r8
         358             :    endif
         359             :
         360           0 :    if (use_hetfrz_classnuc) nimey = 0._r8
         361             :
         362           0 :    nuci=ni + nimey
         363             :
         364           0 :    if(nuci.gt.9999._r8.or.nuci.lt.0._r8) then
         365           0 :       write(iulog, *) 'Warning: incorrect ice nucleation number (nuci reset =0)'
         366           0 :       write(iulog, *) ni, tair, relhum, wbar, nihf, niimm, nidep,deles,esi,dst_num,so4_num
         367           0 :       nuci=0._r8
         368             :    endif
         369             :
         370           0 :    nuci   = nuci*1.e+6_r8/rhoair    ! change unit from #/cm3 to #/kg
         371           0 :    onimey = nimey*1.e+6_r8/rhoair
         372           0 :    onidep = nidep*1.e+6_r8/rhoair
         373           0 :    oniimm = niimm*1.e+6_r8/rhoair
         374           0 :    onihf  = nihf*1.e+6_r8/rhoair
         375           0 : end subroutine nucleati
         376             :
         377             : !===============================================================================
         378             :
         379           0 : subroutine hetero(T,ww,Ns,Nis,Nid)
         380             :
         381             :     real(r8), intent(in)  :: T, ww, Ns
         382             :     real(r8), intent(out) :: Nis, Nid
         383             :
         384             :     real(r8) A11,A12,A21,A22,B11,B12,B21,B22
         385             :     real(r8) B,C
         386             :
         387             : !---------------------------------------------------------------------
         388             : ! parameters
         389             :
         390           0 :       A11 = 0.0263_r8
         391           0 :       A12 = -0.0185_r8
         392           0 :       A21 = 2.758_r8
         393           0 :       A22 = 1.3221_r8
         394           0 :       B11 = -0.008_r8
         395           0 :       B12 = -0.0468_r8
         396           0 :       B21 = -0.2667_r8
         397           0 :       B22 = -1.4588_r8
         398             :
         399             : !     ice from immersion nucleation (cm^-3)
         400             :
         401           0 :       B = (A11+B11*log(Ns)) * log(ww) + (A12+B12*log(Ns))
         402           0 :       C =  A21+B21*log(Ns)
         403             :
         404           0 :       Nis = exp(A22) * Ns**B22 * exp(B*T) * ww**C
         405           0 :       Nis = min(Nis,Ns)
         406             :
         407           0 :       Nid = 0.0_r8    ! don't include deposition nucleation for cirrus clouds when T<-37C
         408             :
         409           0 : end subroutine hetero
         410             :
         411             : !===============================================================================
         412             :
         413           0 : subroutine hf(T,ww,RH,Na,Ni)
         414             :
         415             :       real(r8), intent(in)  :: T, ww, RH, Na
         416             :       real(r8), intent(out) :: Ni
         417             :
         418             :       real(r8)    A1_fast,A21_fast,A22_fast,B1_fast,B21_fast,B22_fast
         419             :       real(r8)    A2_fast,B2_fast
         420             :       real(r8)    C1_fast,C2_fast,k1_fast,k2_fast
         421             :       real(r8)    A1_slow,A2_slow,B1_slow,B2_slow,B3_slow
         422             :       real(r8)    C1_slow,C2_slow,k1_slow,k2_slow
         423             :       real(r8)    regm
         424             :       real(r8)    A,B,C
         425             :       real(r8)    RHw
         426             :
         427             : !---------------------------------------------------------------------
         428             : ! parameters
         429             :
         430           0 :       A1_fast  =0.0231_r8
         431           0 :       A21_fast =-1.6387_r8  !(T>-64 deg)
         432           0 :       A22_fast =-6.045_r8   !(T<=-64 deg)
         433           0 :       B1_fast  =-0.008_r8
         434           0 :       B21_fast =-0.042_r8   !(T>-64 deg)
         435           0 :       B22_fast =-0.112_r8   !(T<=-64 deg)
         436           0 :       C1_fast  =0.0739_r8
         437           0 :       C2_fast  =1.2372_r8
         438             :
         439           0 :       A1_slow  =-0.3949_r8
         440           0 :       A2_slow  =1.282_r8
         441           0 :       B1_slow  =-0.0156_r8
         442           0 :       B2_slow  =0.0111_r8
         443           0 :       B3_slow  =0.0217_r8
         444           0 :       C1_slow  =0.120_r8
         445           0 :       C2_slow  =2.312_r8
         446             :
         447           0 :       Ni = 0.0_r8
         448             :
         449             : !----------------------------
         450             : !RHw parameters
         451           0 :       A = 6.0e-4_r8*log(ww)+6.6e-3_r8
         452           0 :       B = 6.0e-2_r8*log(ww)+1.052_r8
         453           0 :       C = 1.68_r8  *log(ww)+129.35_r8
         454           0 :       RHw=(A*T*T+B*T+C)*0.01_r8
         455             :
         456           0 :       if((T.le.-37.0_r8) .and. ((RH).ge.RHw)) then
         457             :
         458           0 :         regm = 6.07_r8*log(ww)-55.0_r8
         459             :
         460           0 :         if(T.ge.regm) then    ! fast-growth regime
         461             :
         462           0 :           if(T.gt.-64.0_r8) then
         463             :             A2_fast=A21_fast
         464             :             B2_fast=B21_fast
         465             :           else
         466           0 :             A2_fast=A22_fast
         467           0 :             B2_fast=B22_fast
         468             :           endif
         469             :
         470           0 :           k1_fast = exp(A2_fast + B2_fast*T + C2_fast*log(ww))
         471           0 :           k2_fast = A1_fast+B1_fast*T+C1_fast*log(ww)
         472             :
         473           0 :           Ni = k1_fast*Na**(k2_fast)
         474           0 :           Ni = min(Ni,Na)
         475             :
         476             :         else       ! slow-growth regime
         477             :
         478           0 :           k1_slow = exp(A2_slow + (B2_slow+B3_slow*log(ww))*T + C2_slow*log(ww))
         479           0 :           k2_slow = A1_slow+B1_slow*T+C1_slow*log(ww)
         480             :
         481           0 :           Ni = k1_slow*Na**(k2_slow)
         482           0 :           Ni = min(Ni,Na)
         483             :
         484             :         endif
         485             :
         486             :       end if
         487             :
         488           0 : end subroutine hf
         489             :
         490             : !===============================================================================
         491             :
         492           0 : SUBROUTINE Vpreice(P_in, T_in, R_in, C_in, S_in, V_out)
         493             :
         494             :    !  based on  Karcher et al. (2006)
         495             :    !  VERTICAL VELOCITY CALCULATED FROM DEPOSITIONAL LOSS TERM
         496             :
         497             :    ! SUBROUTINE arguments
         498             :    REAL(r8), INTENT(in)  :: P_in       ! [Pa],INITIAL AIR pressure
         499             :    REAL(r8), INTENT(in)  :: T_in       ! [K] ,INITIAL AIR temperature
         500             :    REAL(r8), INTENT(in)  :: R_in       ! [m],INITIAL MEAN  ICE CRYSTAL NUMBER RADIUS
         501             :    REAL(r8), INTENT(in)  :: C_in       ! [m-3],INITIAL TOTAL ICE CRYSTAL NUMBER DENSITY, [1/cm3]
         502             :    REAL(r8), INTENT(in)  :: S_in       ! [-],INITIAL ICE SATURATION RATIO;; if <1, use hom threshold Si
         503             :    REAL(r8), INTENT(out) :: V_out      ! [m/s], VERTICAL VELOCITY REDUCTION (caused by preexisting ice)
         504             :
         505             :    ! SUBROUTINE parameters
         506             :    REAL(r8), PARAMETER :: ALPHAc  = 0.5_r8 ! density of ice (g/cm3), !!!V is not related to ALPHAc
         507             :    REAL(r8), PARAMETER :: FA1c    = 0.601272523_r8
         508             :    REAL(r8), PARAMETER :: FA2c    = 0.000342181855_r8
         509             :    REAL(r8), PARAMETER :: FA3c    = 1.49236645E-12_r8
         510             :    REAL(r8), PARAMETER :: WVP1c   = 3.6E+10_r8
         511             :    REAL(r8), PARAMETER :: WVP2c   = 6145.0_r8
         512             :    REAL(r8), PARAMETER :: FVTHc   = 11713803.0_r8
         513             :    REAL(r8), PARAMETER :: THOUBKc = 7.24637701E+18_r8
         514             :    REAL(r8), PARAMETER :: SVOLc   = 3.23E-23_r8    ! SVOL=XMW/RHOICE
         515             :    REAL(r8), PARAMETER :: FDc     = 249.239822_r8
         516             :    REAL(r8), PARAMETER :: FPIVOLc = 3.89051704E+23_r8
         517             :    REAL(r8) :: T,P,S,R,C
         518             :    REAL(r8) :: A1,A2,A3,B1,B2
         519             :    REAL(r8) :: T_1,PICE,FLUX,ALP4,CISAT,DLOSS,VICE
         520             :
         521           0 :    T = T_in          ! K  , K
         522           0 :    P = P_in*1e-2_r8  ! Pa , hpa
         523             :
         524           0 :    IF (S_in.LT.1.0_r8) THEN
         525           0 :       S = 2.349_r8 - (T/259.0_r8) ! homogeneous freezing threshold, according to Ren & McKenzie, 2005
         526             :    ELSE
         527             :       S = S_in                    ! INPUT ICE SATURATION RATIO, -,  >1
         528             :    ENDIF
         529             :
         530           0 :    R     = R_in*1e2_r8   ! m  => cm
         531           0 :    C     = C_in*1e-6_r8  ! m-3 => cm-3
         532           0 :    T_1   = 1.0_r8/ T
         533           0 :    PICE  = WVP1c * EXP(-(WVP2c*T_1))
         534           0 :    ALP4  = 0.25_r8 * ALPHAc
         535           0 :    FLUX  = ALP4 * SQRT(FVTHc*T)
         536           0 :    CISAT = THOUBKc * PICE * T_1
         537           0 :    A1    = ( FA1c * T_1 - FA2c ) * T_1
         538           0 :    A2    = 1.0_r8/ CISAT
         539           0 :    A3    = FA3c * T_1 / P
         540           0 :    B1    = FLUX * SVOLc * CISAT * ( S-1.0_r8 )
         541           0 :    B2    = FLUX * FDc * P * T_1**1.94_r8
         542           0 :    DLOSS = FPIVOLc * C * B1 * R**2 / ( 1.0_r8+ B2 * R )
         543           0 :    VICE  = ( A2 + A3 * S ) * DLOSS / ( A1 * S )  ! 2006,(19)
         544           0 :    V_out = VICE*1e-2_r8  ! cm/s => m/s
         545             :
         546           0 : END SUBROUTINE Vpreice
         547             :
         548           0 : subroutine frachom(Tmean,RHimean,detaT,fhom)
         549             :    ! How much fraction of cirrus might reach Shom
         550             :    ! base on "A cirrus cloud scheme for general circulation models",
         551             :    ! B. Karcher and U. Burkhardt 2008
         552             :
         553             :    real(r8), intent(in)  :: Tmean, RHimean, detaT
         554             :    real(r8), intent(out) :: fhom
         555             :
         556             :    real(r8), parameter :: seta = 6132.9_r8  ! K
         557             :    integer,  parameter :: Nbin=200          ! (Tmean - 3*detaT, Tmean + 3*detaT)
         558             :
         559             :    real(r8) :: PDF_T(Nbin)    ! temperature PDF;  ! PDF_T=0  outside (Tmean-3*detaT, Tmean+3*detaT)
         560             :    real(r8) :: Sbin(Nbin)     ! the fluctuations of Si that are driven by the T variations
         561             :    real(r8) :: Sihom, deta
         562             :    integer  :: i
         563             :
         564           0 :    Sihom = 2.349_r8-Tmean/259.0_r8   ! homogeneous freezing threshold, according to Ren & McKenzie, 2005
         565           0 :    fhom  = 0.0_r8
         566             :
         567           0 :    do i = Nbin, 1, -1
         568             :
         569           0 :       deta     = (i - 0.5_r8 - Nbin/2)*6.0_r8/Nbin   ! PDF_T=0  outside (Tmean-3*detaT, Tmean+3*detaT)
         570           0 :       Sbin(i)  = RHimean*exp(deta*detaT*seta/Tmean**2.0_r8)
         571           0 :       PDF_T(i) = exp(-deta**2.0_r8/2.0_r8)*6.0_r8/(sqrt(2.0_r8*Pi)*Nbin)
         572             :
         573             :
         574           0 :       if (Sbin(i).ge.Sihom) then
         575           0 :          fhom = fhom + PDF_T(i)
         576             :       else
         577             :          exit
         578             :       end if
         579             :    end do
         580             :
         581           0 :    fhom = min(1.0_r8, fhom/0.997_r8)   ! accounting for the finite limits (-3 , 3)
         582           0 : end subroutine frachom
         583             :
         584             : end module nucleate_ice
         585             :

Generated by: LCOV version 1.14