phys/module_ltng_crmpr92.F

   1 ! WRF:MODEL_LAYER:PHYSICS
   2 !
   3 ! Lightning flash rate prediction based on max vert. verlocity. Implemented
   4 ! for resolutions permitting resolved deep convection.
   5 !
   6 ! Price, C., and D. Rind (1992), A Simple Lightning Parameterization for Calculating
   7 !   Global Lightning Distributions, J. Geophys. Res., 97(D9), 9919-9933, doi:10.1029/92JD00719.
   8 !
   9 ! Wong, J., M. Barth, and D. Noone, 2013: Evaluating a lightning parameterization
  10 !   based on cloud-top height for mesoscale numerical model simulations,
  11 !   Geosci. Model Dev., 6, 429–443, https://doi.org/10.5194/gmd-6-429-2013.
  12 !
  13 ! Unlike previous implementation, this version will produce slightly inconsistent
  14 ! IC and CG grid-flash rates against NO emission after production via calling
  15 ! lightning_nox_decaria.
  16 !
  17 !**********************************************************************
  18
  19  MODULE module_ltng_crmpr92
  20  CONTAINS
  21
  22  SUBROUTINE ltng_crmpr92w ( &
  23                           ! Frequently used prognostics
  24                             dx, dy, xland, ht, z, t,              &
  25                           ! Scheme specific prognostics
  26                             w, refl, reflthreshold, cellcount,    &
  27                           ! Scheme specific namelist inputs
  28                             cellcount_method,                     &
  29                           ! Order dependent args for domain, mem, and tile dims
  30                             ids, ide, jds, jde, kds, kde,         &
  31                             ims, ime, jms, jme, kms, kme,         &
  32                             ips, ipe, jps, jpe, kps, kpe,         &
  33                           ! Mandatory output for all quantitative schemes
  34                             total_flashrate                       &
  35                           )
  36 !-----------------------------------------------------------------
  37 ! Framework
  38  USE module_state_description
  39
  40 ! Model layer
  41  USE module_model_constants
  42  USE module_wrf_error
  43
  44  USE module_dm, only: wrf_dm_max_real
  45
  46  IMPLICIT NONE
  47 !-----------------------------------------------------------------
  48
  49 ! Frequently used prognostics
  50  REAL,    INTENT(IN   )    ::       dx, dy
  51
  52  REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) :: xland, ht
  53  REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) :: z, t
  54
  55 ! Scheme specific prognostics
  56  REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) :: w
  57  REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) :: refl
  58  REAL,                                            INTENT(IN   ) :: reflthreshold
  59  REAL,    DIMENSION(          kms:kme          ), INTENT(IN   ) :: cellcount
  60
  61 ! Scheme specific namelist inputs
  62  INTEGER, INTENT(IN   )    ::       cellcount_method
  63
  64 ! Order dependent args for domain, mem, and tile (patch) dims
  65  INTEGER, INTENT(IN   )    ::       ids,ide, jds,jde, kds,kde
  66  INTEGER, INTENT(IN   )    ::       ims,ime, jms,jme, kms,kme
  67  INTEGER, INTENT(IN   )    ::       ips,ipe, jps,jpe, kps,kpe
  68
  69 ! Mandatory outputs for all quantitative schemes
  70  REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(  OUT) :: total_flashrate
  71
  72 ! Local variables
  73  REAL :: wmax            ! max w in patch or domain
  74  REAL :: total_fr,ave_fr ! cloud flash rate
  75  INTEGER :: i,k,j
  76  INTEGER :: k_maxcount
  77  REAL :: maxcount
  78  CHARACTER (LEN=250) :: message
  79
  80 !-----------------------------------------------------------------
  81
  82  total_flashrate( ips:ipe,jps:jpe ) = 0.
  83
  84  IF ( maxval(cellcount(kps:kpe)) .eq. 0 ) RETURN
  85
  86 ! Compute flash rate across cell
  87  wmax = maxval(w(ips:ipe,kps:kpe,jps:jpe))
  88  IF ( cellcount_method .eq. 2 ) THEN
  89    wmax = wrf_dm_max_real(wmax)
  90  ENDIF
  91
  92  total_fr = 5.7e-6 * wmax**4.5
  93
  94 ! Locating widest part of convective core
  95  k_maxcount = kps
  96  maxcount = cellcount(kps)
  97  DO k=kps+1,kpe
  98    IF ( cellcount(k) .gt. maxcount ) THEN
  99      k_maxcount = k
 100      maxcount = cellcount(k)
 101    ENDIF
 102  ENDDO
 103
 104 ! Distributing across convective core
 105  ave_fr = total_fr/maxcount/60.
 106  WHERE( refl(ips:ipe,k_maxcount,jps:jpe) .gt. reflthreshold )
 107    total_flashrate(ips:ipe,jps:jpe) = ave_fr
 108  ENDWHERE
 109
 110  END SUBROUTINE ltng_crmpr92w
 111
 112  SUBROUTINE ltng_crmpr92z ( &
 113                           ! Frequently used prognostics
 114                             dx, dy, xland, ht, z, t,              &
 115                           ! Scheme specific prognostics
 116                             refl, reflthreshold, cellcount,       &
 117                           ! Scheme specific namelist inputs
 118                             cellcount_method,                     &
 119                           ! Order dependent args for domain, mem, and tile dims
 120                             ids, ide, jds, jde, kds, kde,         &
 121                             ims, ime, jms, jme, kms, kme,         &
 122                             ips, ipe, jps, jpe, kps, kpe,         &
 123                           ! Mandatory output for all quantitative schemes
 124                             total_flashrate                       &
 125                           )
 126 !-----------------------------------------------------------------
 127 ! Framework
 128  USE module_state_description
 129
 130 ! Model layer
 131  USE module_model_constants
 132  USE module_wrf_error
 133
 134  USE module_dm, only: wrf_dm_max_real
 135
 136  IMPLICIT NONE
 137 !-----------------------------------------------------------------
 138
 139 ! Frequently used prognostics
 140  REAL,    INTENT(IN   )    ::       dx, dy
 141
 142  REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) :: xland, ht
 143  REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) :: z, t
 144
 145 ! Scheme specific prognostics
 146  REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) :: refl
 147  REAL,                                            INTENT(IN   ) :: reflthreshold
 148  REAL,    DIMENSION(          kms:kme          ), INTENT(IN   ) :: cellcount
 149
 150 ! Scheme specific namelist inputs
 151  INTEGER, INTENT(IN   )    ::       cellcount_method
 152
 153 ! Order dependent args for domain, mem, and tile (patch) dims
 154  INTEGER, INTENT(IN   )    ::       ids,ide, jds,jde, kds,kde
 155  INTEGER, INTENT(IN   )    ::       ims,ime, jms,jme, kms,kme
 156  INTEGER, INTENT(IN   )    ::       ips,ipe, jps,jpe, kps,kpe
 157
 158 ! Mandatory outputs for all quantitative schemes
 159  REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(  OUT) :: total_flashrate
 160
 161 ! Local variables
 162  REAL :: zmax            ! max w in patch or domain
 163  REAL :: total_fr,ave_fr ! cloud flash rate
 164  INTEGER :: i,k,j
 165  INTEGER :: k_maxcount, count
 166  REAL :: maxcount, mostlyLand
 167  CHARACTER (LEN=250) :: message
 168
 169 !-----------------------------------------------------------------
 170
 171  total_flashrate( ips:ipe,jps:jpe ) = 0.
 172
 173  IF ( maxval(cellcount(kps:kpe)) .eq. 0 ) RETURN
 174
 175 ! Compute flash rate across cell
 176  k = kpe
 177  do while ( cellcount(k) .eq. 0 .and. k .gt. kps)
 178    k = k-1
 179  ENDDO
 180  zmax = 0.
 181  mostlyland = 0.
 182  count = 0
 183  DO i=ips,ipe
 184    DO j=jps,jpe
 185      IF ( (refl(i,k,j) .gt. reflthreshold) .and. (t(i,k,j) .lt. 273.15) ) THEN
 186        IF (z(i,k,j)-ht(i,j) .gt. zmax) THEN
 187          zmax = z(i,k,j)-ht(i,j)
 188        ENDIF
 189        count = count + 1
 190        mostlyland = mostlyland + xland(i,j)
 191      ENDIF
 192    ENDDO
 193  ENDDO
 194  mostlyland = mostlyland/count
 195
 196  zmax = zmax * 1.e-3
 197  WRITE(message, * ) ' ltng_crmpr92z: reflectivity cloud top height: ', zmax
 198  CALL wrf_debug ( 15, message )
 199
 200  if ( cellcount_method .eq. 2 ) THEN
 201    zmax = wrf_dm_max_real(zmax)
 202  endif
 203
 204  if ( mostlyLand .lt. 1.5 ) then
 205     total_fr = 3.44E-5 * (zmax**4.9)  ! PR 92 continental eq
 206  else
 207     total_fr = 6.57E-6 * (zmax**4.9)  ! Michalon 99 marine eq
 208  ENDIF
 209
 210 ! Locating widest part of convective core
 211  k_maxcount = kps
 212  maxcount = cellcount(kps)
 213  DO k=kps+1,kpe
 214    IF ( cellcount(k) .gt. maxcount ) THEN
 215      k_maxcount = k
 216      maxcount = cellcount(k)
 217    ENDIF
 218  ENDDO
 219
 220 ! Distributing across convective core
 221  ave_fr = total_fr/maxcount/60.
 222  WHERE( refl(ips:ipe,k_maxcount,jps:jpe) .gt. reflthreshold  )
 223    total_flashrate(ips:ipe,jps:jpe) = ave_fr
 224  ENDWHERE
 225
 226  END SUBROUTINE ltng_crmpr92z
 227
 228 !**********************************************************************
 229 !
 230 ! Price and Rind 1993 base on cold cloud depth (CCD)
 231 !
 232 ! Price, C. and D. Rind (1993), What determines the cloud-to-ground lightning
 233 ! fraction in thunderstorms?, Geophys. Res. Lett., 20(6), 463-466, doi:10.1029/93GL00226.
 234 !
 235 ! Valid range of CCD is set to 5.5-14 km. Beyond this range CCD is assumed
 236 ! to be 5.5 or 14 for continuity.
 237 !
 238 !**********************************************************************
 239  SUBROUTINE iccg_crm_pr93( &
 240                             refl, reflthreshold, t, z,                 &
 241                           ! Order dependent args for domain, mem, and tile dims
 242                             ids, ide, jds, jde, kds, kde,              &
 243                             ims, ime, jms, jme, kms, kme,              &
 244                             ips, ipe, jps, jpe, kps, kpe,              &
 245                           ! Input
 246                             total_flashrate,                           &
 247                           ! Output
 248                             ic_flashrate, cg_flashrate                 &
 249                         )
 250 !-----------------------------------------------------------------
 251  IMPLICIT NONE
 252 !-----------------------------------------------------------------
 253 ! Inputs
 254  REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) :: refl, t, z
 255  REAL,                                            INTENT(IN   ) :: reflthreshold
 256
 257 ! Order dependent args for domain, mem, and tile dims
 258  INTEGER, INTENT(IN   )    ::       ids,ide, jds,jde, kds,kde
 259  INTEGER, INTENT(IN   )    ::       ims,ime, jms,jme, kms,kme
 260  INTEGER, INTENT(IN   )    ::       ips,ipe, jps,jpe, kps,kpe
 261
 262 ! Primary inputs and outpus
 263  REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) :: total_flashrate
 264  REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(  OUT) :: ic_flashrate, cg_flashrate
 265
 266 ! Local variables
 267  INTEGER :: kfreeze, ktop
 268
 269  INTEGER :: i,j,k
 270  REAL    :: ratio, cgfrac, depth
 271
 272  REAL, PARAMETER :: dH_min = 5.5
 273  REAL, PARAMETER :: dH_max = 14.
 274
 275  REAL, PARAMETER :: coef_A = 0.021
 276  REAL, PARAMETER :: coef_B = -0.648
 277  REAL, PARAMETER :: coef_C = 7.493
 278  REAL, PARAMETER :: coef_D = -36.54
 279  REAL, PARAMETER :: coef_E = 63.09
 280 !-----------------------------------------------------------------
 281
 282  ic_flashrate(ips:ipe,jps:jpe) = 0.
 283  cg_flashrate(ips:ipe,jps:jpe) = 0.
 284
 285  jloop: DO j=jps,jpe
 286     iloop: DO i=ips,ipe
 287     IF ( total_flashrate(i,j) .gt. 0.) THEN
 288         ktop = kpe
 289         do while ( refl(i,ktop,j) .lt. reflthreshold .and. ktop .gt. kps)
 290           ktop = ktop-1
 291         enddo
 292
 293         kfreeze = ktop
 294         DO WHILE ( t(i,kfreeze,j) .lt. 273.15 .and. kfreeze .gt. kps )
 295             kfreeze = kfreeze - 1
 296         ENDDO
 297
 298         depth = ( z(i,ktop,j) - z(i,kfreeze,j) ) * 1E-3
 299         IF (depth .le. 0.) CONTINUE
 300         depth = max( dH_min, min( dH_max, depth ))
 301
 302         ratio = (((coef_A*depth+coef_B )*depth+coef_C)*depth+coef_D)*depth+coef_E
 303         cgfrac = 1./(ratio+1.)
 304
 305         cg_flashrate(i,j) = total_flashrate(i,j) * cgfrac
 306         ic_flashrate(i,j) = total_flashrate(i,j) - cg_flashrate(i,j)
 307     ENDIF
 308     ENDDO iloop
 309  ENDDO jloop
 310
 311  END SUBROUTINE iccg_crm_pr93
 312
 313  END MODULE module_ltng_crmpr92