phys/module_ra_effective_radius.F

   1 !
   2 !-------------------------------------------------------------------------------
   3    module module_ra_effective_radius
   4 !-------------------------------------------------------------------------------
   5 !
   6 !  abstract :
   7 !    calculate effective radius of cloud water, cloud ice, and snow
   8 !             based on WSM/DM microphysics scheme
   9 !
  10 !  history log:
  11 !   2015-07-01  soo ya bae   initial setup
  12 !   2015-10-01  soo ya bae   move effective radius to radiation
  13 !                            add subroutine for consistency among rad-cps-mps
  14 !   2016-01-01  soo ya bae   bug fix
  15 !
  16 !  reference :
  17 !    bae et al. (2016, adv meteorol)
  18 !
  19 !-------------------------------------------------------------------------------
  20    use module_model_constants,  only: denr=>rhowater, dens=>rhosnow
  21 !-------------------------------------------------------------------------------
  22    real, parameter, private :: n0s = 2.e6      ! temperature dependent
  23                                                ! intercept parameter of snow
  24    real, parameter, private :: n0g = 4.e6
  25    real, parameter, private :: n0smax =  1.e11 ! maximum n0s (t=-90C unlimited)
  26    real, parameter, private :: nc0 = 3.e8      ! number conc. of cloud droplet
  27    real, parameter, private :: deng = 500.0    ! density of graupel
  28    real, parameter, private :: alpha = .12     ! .122 exponen factor for n0s
  29    real, parameter, private :: pi = 4.*atan(1.)
  30 !
  31    contains
  32 !-------------------------------------------------------------------------------
  33 !
  34 !
  35 !-------------------------------------------------------------------------------
  36    real function rgmma(x)
  37 !-------------------------------------------------------------------------------
  38 !
  39 !  abstract : rgmma function
  40 !             use infinite product form
  41 !
  42 !-------------------------------------------------------------------------------
  43 !
  44    implicit none
  45 !
  46    real, parameter :: euler = 0.577215664901532
  47    real            :: x, y
  48    integer         :: i
  49 !-------------------------------------------------------------------------------
  50    if (x.eq.1.) then
  51      rgmma = 0.
  52    else
  53      rgmma = x*exp(euler*x)
  54      do i = 1,10000
  55        y = real(i)
  56        rgmma = rgmma*(1.000+x/y)*exp(-x/y)
  57      enddo
  58      rgmma = 1./rgmma
  59    endif
  60 !
  61    end function rgmma
  62 !-------------------------------------------------------------------------------
  63 !
  64 !
  65 !-------------------------------------------------------------------------------
  66    subroutine effectRad(t, qc, nc, qi, qs, qg, rho, qmin, t0c,                 &
  67                         qccps, f_qnc,                                          &
  68                         re_qc, re_qi, re_qs, kts, kte)
  69 !-------------------------------------------------------------------------------
  70 !
  71 !  abstract :
  72 !    compute radiation effective radii of cloud water, ice, and snow
  73 !             for double-moment microphysics.
  74 !
  75 !  input  :
  76 !    kts, kte   - dimension
  77 !    ii, jj
  78 !    qmin       - minimum value of hydrometeor in kg/kg
  79 !    t0c        - 273.15 K
  80 !    t          - temperature
  81 !    qc         - mixing ratio of cloud water in kg/kg
  82 !    qi         - mixing ratio of cloud ice in kg/kg
  83 !    qs         - mixing ratio of snow in kg/kg
  84 !    qg         - mixing ratio of graupel in kg/kg
  85 !    nc         - number concentratino of cloud water in 1/m3
  86 !    rho        - air density of kg/m3
  87 !
  88 !  inout :
  89 !    re_qc     - effective radius of cloud water in m
  90 !    re_qi     - effective radius of cloud ice in m
  91 !    re_qs     - effective radius of snow in m
  92 !
  93 !-------------------------------------------------------------------------------
  94 !
  95    implicit none
  96 !
  97    integer                 , intent(in   ) :: kts, kte
  98    real                    , intent(in   ) :: qmin, t0c
  99    real, dimension(kts:kte), intent(in   ) :: t
 100    real, dimension(kts:kte), intent(in   ) :: qc
 101    real, dimension(kts:kte), intent(in   ) :: qi
 102    real, dimension(kts:kte), intent(in   ) :: qs
 103    real, dimension(kts:kte), intent(in   ) :: qg
 104    real, dimension(kts:kte), intent(in   ) :: nc
 105    real, dimension(kts:kte), intent(in   ) :: qccps
 106    real, dimension(kts:kte), intent(in   ) :: rho
 107    real, dimension(kts:kte), intent(inout) :: re_qc
 108    real, dimension(kts:kte), intent(inout) :: re_qi
 109    real, dimension(kts:kte), intent(inout) :: re_qs
 110    logical                 , intent(in   ) :: f_qnc
 111 !
 112 ! local variables
 113 !
 114    integer                  :: i,k
 115    integer                  :: kte_in
 116    integer                  :: index
 117    real                     :: cdm2
 118    real                     :: temp
 119    real                     :: supcol, n0sfac, lamdas
 120    real                     :: diai      ! diameter of ice in m
 121    real                     :: corr
 122    real                     :: pidnc, pidn0s
 123    real                     :: pidn0g
 124    real                     :: lamdag
 125    double precision         :: lamc
 126    double precision         :: lammps
 127    logical                  :: has_qc, has_qi, has_qs
 128    real, dimension(kts:kte) :: ni
 129    real, dimension(kts:kte) :: rqc
 130    real, dimension(kts:kte) :: rnc
 131    real, dimension(kts:kte) :: rqi
 132    real, dimension(kts:kte) :: rni
 133    real, dimension(kts:kte) :: rqs
 134    real, dimension(kts:kte) :: rqg
 135    real, dimension(kts:kte) :: qsqg
 136    real, dimension(kts:kte) :: re_qg
 137    real, dimension(kts:kte) :: qcmps
 138    real, dimension(kts:kte) :: rqcmps
 139    real, dimension(kts:kte) :: nccps
 140 !
 141 ! minimum microphys values
 142 !
 143    real, parameter          :: r0 = 0.
 144    real, parameter          :: r1 = 1.e-12
 145    real, parameter          :: r2 = 1.e-6
 146 !
 147 ! mass power law relations:  mass = am*d**bm
 148 !
 149    real, parameter          :: bm_r = 3.0
 150    real, parameter          :: obmr = 1.0/bm_r
 151    real, parameter          :: cdm  = 5./3.
 152 !-------------------------------------------------------------------------------
 153    has_qc = .false.
 154    has_qi = .false.
 155    has_qs = .false.
 156
 157    ni=1.e3
 158    rnc=0.
 159    rni=0.
 160    rqc=0.
 161    rqi=0.
 162    rqs=0.
 163    rqg=0.
 164    qsqg=0.
 165    re_qg=25.e-6
 166    qcmps=0.
 167    rqcmps=0.
 168    nccps=0.
 169
 170 !
 171    kte_in = kte
 172    pidnc  = pi*denr/6.
 173    pidn0s = pi*dens*n0s
 174    pidn0g = pi*deng*n0g
 175    cdm2 = rgmma(cdm)
 176 !
 177    do k = kts,kte_in
 178 !
 179 ! for cloud
 180 !
 181      rqc(k) = max(r1,qc(k)*rho(k))
 182      qcmps(k)  = qc(k)-qccps(k)
 183      rqcmps(k) = max(r1,qcmps(k)*rho(k))
 184      if(f_qnc) then
 185        lammps    = 2.*cdm2*(pidnc*MAX(r0,nc(k))/rqcmps(k))**obmr
 186      else
 187        lammps   = (pidnc*nc0/rqcmps(k))**obmr
 188      endif
 189      nccps(k)  = qccps(k)*6./pi*denr/rho(k)*lammps**3.
 190      rnc(k)    = max(r2,(nc(k)+nccps(k))*rho(k))
 191      if (rqc(k).gt.r1.and.rnc(k).gt.r2) has_qc = .true.
 192 !
 193 ! for ice
 194 !
 195      rqi(k) = max(r1,qi(k)*rho(k))
 196      temp   = (rho(k)*max(qi(k),qmin))
 197      temp   = sqrt(sqrt(temp*temp*temp))
 198      ni(k)  = min(max(5.38e7*temp,1.e3),1.e6)
 199      rni(k) = max(r2,ni(k)*rho(k))
 200      if (rqi(k).gt.r1.and.rni(k).gt.r2) has_qi = .true.
 201 !
 202 ! for snow
 203 !
 204      rqs(k) = max(r1,qs(k)*rho(k))
 205      rqg(k) = max(r1,qg(k)*rho(k))
 206      if (rqs(k).gt.r1.or.rqg(k).gt.r1) has_qs = .true.
 207    enddo
 208 !
 209    if (has_qc) then
 210      do k = kts,kte_in
 211        if (rqc(k).le.r1.or.rnc(k).le.r2) cycle
 212        lamc = 2.*cdm2*(pidnc*(MAX(r0,nc(k))+nccps(k))/rqc(k))**obmr
 213        re_qc(k) =  max(2.51e-6,min(sngl(3.0/lamc),50.e-6))
 214      enddo
 215    endif
 216 !
 217    if (has_qi) then
 218      do k = kts,kte_in
 219        if (rqi(k).le.r1.or.rni(k).le.r2) cycle
 220        diai = 11.9*sqrt(rqi(k)/ni(k))
 221        re_qi(k) = max(10.01e-6,min(0.75*0.163*diai,125.e-6))
 222      enddo
 223    endif
 224 !
 225    if (has_qs) then
 226      do k = kts,kte_in
 227        if (rqs(k).le.r1) cycle
 228        supcol = t0c-t(k)
 229        n0sfac = max(min(exp(alpha*supcol),n0smax/n0s),1.)
 230        lamdas = sqrt(sqrt(pidn0s*n0sfac/rqs(k)))
 231        re_qs(k) = max(25.e-6,min(0.5*(1./lamdas),999.e-6))
 232      enddo
 233 !
 234      do k = kts,kte_in
 235        if (rqg(k).le.r1) cycle
 236        lamdag = sqrt(sqrt(pidn0g*n0g/rqg(k)))
 237        re_qg(k) = max(25.e-6,min(1.5*(1./lamdag),999.e-6))
 238      enddo
 239 !
 240      do k = kts,kte_in
 241        qsqg(k) = max(r1,qs(k)+qg(k))
 242        re_qs(k) = (re_qs(k)*max(r0,qs(k))+re_qg(k)*max(r0,qg(k)))/qsqg(k)
 243      enddo
 244    endif
 245 !
 246    end subroutine effectRad
 247 !-------------------------------------------------------------------------------
 248 !
 249 !
 250 !-------------------------------------------------------------------------------
 251    end module module_ra_effective_radius
 252 !-------------------------------------------------------------------------------