chem/module_lightning_nox_decaria.F

   1 !WRF:MODEL_LAYER:CHEMISTRY
   2 !
   3 ! Contains subroutine for converting flash rates into NO emissions
   4 ! based on Decaria 2000 vertical distirbutions.
   5 !
   6 ! Input: flashes (#/s)
   7 ! Output: tendency (ppmv/s)
   8 !
   9 ! The output will be muliplied by timestep and used to incremeent NO
  10 ! concentration and the respective passive tracer in lightning_nox_driver.
  11 !
  12 ! See module_lightning_nox_driver for more info.
  13 !
  14 ! Contact: M. Barth <barthm@ucar.edu>, J. Wong <johnwong@ucar.edu>
  15 !
  16 !**********************************************************************
  17  MODULE module_lightning_nox_decaria
  18
  19  IMPLICIT NONE
  20
  21  CONTAINS
  22
  23 !**********************************************************************
  24 !
  25 ! DeCaria et al, 2000
  26 !
  27 ! DeCaria, A. J., K. E. Pickering, G. L. Stenchikov, and L. E. Ott (2005),
  28 ! Lightning-generated NOX and its impact on tropospheric ozone production:
  29 ! A three-dimensional modeling study of a Stratosphere-Troposphere Experiment:
  30 ! Radiation, Aerosols and Ozone (STERAO-A) thunderstorm, J. Geophys. Res.,
  31 ! 110, D14303, doi:10.1029/2004JD005556.
  32 !
  33 !**********************************************************************
  34  SUBROUTINE lightning_nox_decaria ( &
  35                           ! Frequently used prognostics
  36                             dx, dy, xland, ht, t, rho, z, p,      &
  37                             ic_flashrate, cg_flashrate,           & ! flashes (#/s)
  38                           ! Scheme specific prognostics
  39                             refl,                                 &
  40                           ! Namelist inputs
  41                             N_IC, N_CG,                           &
  42                             ltng_temp_upper,ltng_temp_lower,      &
  43                             cellcount_method,                     &
  44                           ! Order dependent args for domain, mem, and tile dims
  45                             ids, ide, jds, jde, kds, kde,         &
  46                             ims, ime, jms, jme, kms, kme,         &
  47                             ips, ipe, jps, jpe, kps, kpe,         &
  48                           ! outputs
  49                             lnox_ic_tend, lnox_cg_tend            & ! tendency (ppmv/s)
  50                           )
  51 !-----------------------------------------------------------------
  52 ! Framework
  53  USE module_state_description
  54
  55 ! Model layer
  56  USE module_model_constants
  57  USE module_wrf_error
  58
  59  USE module_dm, only: wrf_dm_max_real, wrf_dm_min_real, wrf_dm_sum_real
  60
  61 ! Lightning method
  62  USE module_lightning_driver, only: countCells
  63
  64  IMPLICIT NONE
  65 !-----------------------------------------------------------------
  66
  67 ! Frequently used prognostics
  68  REAL,    INTENT(IN   )    ::       dx, dy
  69
  70  REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) :: xland, ht
  71  REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) :: t, rho, z, p
  72  REAL,    DIMENSION( ims:ime,          jms:jme ), INTENT(IN   ) :: ic_flashrate  , cg_flashrate ! #/sec
  73
  74
  75 ! Scheme specific prognostics
  76  REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   ) :: refl
  77
  78 ! Scheme specific namelist inputs
  79  REAL,    INTENT(IN   )    ::       N_IC, N_CG
  80  REAL,    INTENT(IN   )    ::       ltng_temp_lower, ltng_temp_upper
  81  INTEGER, INTENT(IN   )    ::       cellcount_method
  82
  83 ! Order dependent args for domain, mem, and tile (patch) dims
  84  INTEGER, INTENT(IN   )    ::       ids,ide, jds,jde, kds,kde
  85  INTEGER, INTENT(IN   )    ::       ims,ime, jms,jme, kms,kme
  86  INTEGER, INTENT(IN   )    ::       ips,ipe, jps,jpe, kps,kpe
  87
  88 ! Mandatory outputs for all quantitative schemes
  89  REAL,    DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(  OUT) :: lnox_ic_tend,lnox_cg_tend
  90
  91 ! Local variables
  92  INTEGER :: i,k,j
  93  INTEGER :: ktop,kbtm,kupper,klower
  94  REAL :: ic_fr, cg_fr, delta ! reconsolidated flashrates
  95  REAL :: reflmax, cellmax
  96  REAL :: term2, B, B_denom
  97  CHARACTER (LEN=250) :: message
  98  REAL, DIMENSION( kps:kpe ) :: cellcount
  99  REAL, DIMENSION( kps:kpe ) :: z_average, t_average, p_average, rho_average, conv
 100  REAL, DIMENSION( kps:kpe ) :: fd, fd2, dz ! fd = distribution
 101
 102  REAL, PARAMETER :: refl_threshold = 20.
 103
 104 !-----------------------------------------------------------------
 105
 106  lnox_ic_tend (ips:ipe,kps:kpe,jps:jpe ) = 0.
 107  lnox_cg_tend (ips:ipe,kps:kpe,jps:jpe ) = 0.
 108
 109 ! Determine cloud extents. Also calculated in physics but cellcount is not persistent
 110  CALL countCells( &
 111       ! Inputs
 112         refl, refl_threshold, cellcount_method,     &
 113       ! Order dependent args for domain, mem, and tile dims
 114         ids, ide, jds, jde, kds, kde,              &
 115         ims, ime, jms, jme, kms, kme,              &
 116         ips, ipe, jps, jpe, kps, kpe,              &
 117       ! Output
 118         cellcount )
 119
 120 ! Reconsolidate flash counts
 121  ic_fr = sum(ic_flashrate(ips:ipe,jps:jpe))
 122  cg_fr = sum(cg_flashrate(ips:ipe,jps:jpe))
 123  if ( cellcount_method .eq. 2 ) then
 124     ic_fr = wrf_dm_sum_real(ic_fr)
 125     cg_fr = wrf_dm_sum_real(cg_fr)
 126  ENDIF
 127  reflmax = maxval(refl(ips:ipe,kps:kpe,jps:jpe))
 128  cellmax = maxval(cellcount(kps:kpe))
 129  WRITE(message, * ) ' LNOx tracer: max_refl, max_cellcount, ic_fr = ',  reflmax, cellmax, ic_fr
 130  CALL wrf_debug ( 100, message )
 131
 132 !-----------------------------------------------------------------
 133
 134 ! Average z, t, p, rho
 135  CALL horizontalAverage( z( ips:ipe,kps:kpe,jps:jpe ), ips, ipe, kps, kpe, jps, jpe, z_average )
 136  CALL horizontalAverage( t( ips:ipe,kps:kpe,jps:jpe ), ips, ipe, kps, kpe, jps, jpe, t_average )
 137  CALL horizontalAverage( p( ips:ipe,kps:kpe,jps:jpe ), ips, ipe, kps, kpe, jps, jpe, p_average )
 138  CALL horizontalAverage( rho( ips:ipe,kps:kpe,jps:jpe ), ips, ipe, kps, kpe, jps, jpe, rho_average )
 139
 140 ! molesofair(kps:kpe) = rho_average(kps:kpe) * 1E3 * dx * dy / .02897 ! # moles per km in z
 141 ! term2 = 30 * 8.3145E6/dx/dy/28.96/100./100.
 142
 143  conv(kps:kpe) = 8.314 *t_average(kps:kpe) / (dx * dy)                !  conversion term with units J/(mol-m2)
 144
 145
 146  CALL  kfind ( cellcount, t_average,            &
 147                ltng_temp_upper,ltng_temp_lower, cellcount_method, &
 148                 ips, ipe, jps, jpe, kps, kpe,              &
 149               ! Outputs
 150                 ktop,kbtm,kupper,klower )
 151
 152 ! Calculates IC distribution
 153 !IF (( ic_fr .gt. 0 ) .and. (( ktop .gt. klower ) .and. (kbtm .lt. ktop) ) )THEN
 154  IF (( ic_fr > 0 ) .and. (( ktop > klower ) .and. (kbtm < klower) ) )THEN
 155    call bellcurve(kbtm,ktop,klower,z_average, kps,kpe, fd, dz)
 156    if (ktop .gt. kupper) then
 157      call bellcurve(kbtm,ktop,kupper,z_average, kps,kpe, fd2, dz)
 158      fd(kbtm:ktop) = 0.5*( fd(kbtm:ktop) + fd2(kbtm:ktop) )         ! unitless
 159    endif
 160 !   B = N_IC/sum(f(kbtm:ktop)*p_average(kbtm:ktop))             ! *** used in calculating NO
 161    B_denom = DOT_PRODUCT( fd(kbtm:ktop),p_average(kbtm:ktop) )      ! N/m2
 162
 163
 164    DO k=kbtm,ktop
 165      if ( cellcount(k) .gt. 0. ) THEN
 166 !       delta = term2*B*fd(k)*t_average(k)*ic_fr/cellcount(k)/dz(k)/100.
 167         !*  implementation note: 1) ic_fr * N_IC/cellcount gives moles of NO in the column
 168         !*                       2) Multiplying by fd gives the # moles per level
 169         !*                       3) Convert to mol NO/mol air per minute
 170         !*                       4) Multiply by 1E6 gives ppmv
 171
 172         delta = (ic_fr * N_IC / cellcount(k)) * fd(k) / B_denom * conv(k)/dz(k) * 1E6
 173 !units:      flash/sec * mol/flash /()       * m2/ N      * J/(mol-m2) / m * ppmv/(mol NO/mol air)
 174 !units:      flash/sec * mol/flash           * m2/ N      * N-m/(mol-m3) * ppmv/(mol NO/mol air)
 175 !units:      ppmv/sec
 176
 177         where(refl(ips:ipe,k,jps:jpe) .gt. refl_threshold )
 178           lnox_ic_tend(ips:ipe,k,jps:jpe) = delta
 179         endwhere
 180      ENDIF
 181    ENDDO
 182
 183  ENDIF ! IC lightning
 184
 185 !-----------------------------------------------------------------
 186 ! Calculates CG distribution
 187 !IF ((cg_fr .gt. 0 ) .and. (( ktop .gt. klower ) .and. (kbtm .lt. ktop) ) ) THEN
 188  IF ((cg_fr > 0 ) .and. (( ktop > klower ) .and. (kbtm < klower) ) ) THEN
 189    call bellcurve(kps,ktop,klower,z_average, kps,kpe, fd, dz)
 190 !   B = N_CG/(sum(fd(kps:ktop)*p_average(kps:ktop)))
 191
 192    B_denom = DOT_PRODUCT( fd(kbtm:ktop),p_average(kbtm:ktop) )      ! N/m2
 193
 194    k = ktop
 195
 196    DO WHILE (k .ge. kps)
 197      IF (cellcount(k) .gt. 0) THEN
 198
 199 !      delta = (cg_fr * N_CG / cellcount(k)) * fd(k) / (molesofair(k)*dz(k)) * 1E6
 200       delta = (cg_fr * N_CG / cellcount(k)) * fd(k) / B_denom * conv(k)/dz(k) * 1E6
 201 !units:      flash/sec * mol/flash /()       * m2/ N      * J/(mol-m2) / m * ppmv/(mol NO/mol air)
 202 !units:      flash/sec * mol/flash           * m2/ N      * N-m/(mol-m3) * ppmv/(mol NO/mol air)
 203 !units:      ppmv/sec
 204
 205        where( refl(ips:ipe,k,jps:jpe) .gt. refl_threshold )
 206           lnox_cg_tend(ips:ipe,k,jps:jpe) = delta
 207        endwhere
 208
 209      ENDIF
 210
 211      k = k - 1      !07/23/14 KAC added
 212
 213    ENDDO
 214
 215  ENDIF
 216
 217  END SUBROUTINE lightning_nox_decaria
 218
 219
 220
 221 !************************************************************************
 222 ! This subroutine prepares a normal distribution between k_min and
 223 ! k_max centered at k_mu. Distribution for each level is
 224 ! normalized to \int^{z_at_w(k_max}_{z_at_w(k_min-1)}f(z)dz
 225 !
 226 ! Unit of f is fraction of total column
 227 !
 228 ! Modified from v3.4.1 module_ltng_crm.F, kept the math but changed
 229 ! the implementation for better clarity. Removed patch-wide averaging
 230 ! of z.
 231 !************************************************************************
 232
 233  SUBROUTINE bellcurve ( k_min, k_max, k_mu, z, kps,kpe, f, dz )
 234 !-----------------------------------------------------------------
 235
 236  IMPLICIT NONE
 237  INTEGER,                      INTENT(IN   ) :: k_min, k_max, k_mu
 238  REAL,   DIMENSION( kps:kpe ), INTENT(IN   ) :: z       ! at phy
 239  INTEGER,                      INTENT(IN   ) :: kps,kpe
 240
 241  REAL,   DIMENSION( kps:kpe ), INTENT(  OUT) :: f, dz
 242
 243  INTEGER :: i,j,k
 244  REAL, DIMENSION( kps:kpe ) :: ex
 245  REAL :: sigma, z_mu, cuml_f_dist
 246  REAL, PARAMETER :: two_pi = 6.2831854
 247
 248 !-----------------------------------------------------------------
 249
 250  f(kps:kpe) = 0.
 251  z_mu = z(k_mu)
 252  sigma = AMIN1(z(k_max)-z_mu,z_mu-z(k_min))/3.0
 253
 254  ! distance from mean
 255  ex(k_min:k_max) = (z(k_min:k_max)-z_mu)/sigma
 256
 257  ! Truncated Gaussian at 3 sigma
 258  f(k_min:k_max) = (1.0/(sqrt(two_pi)*sigma))*exp(-ex(k_min:k_max)*ex(k_min:k_max)/2.0)
 259
 260 !++mcb   We do need dz at bottom and top of domain
 261 !   dz(kps) = 0. ! safe as long as k_min != kps
 262 !   dz(kpe) = 0. ! safe as long as k_max != kpe
 263  dz(kps) = (z(kps+1) - z(kps))*.5
 264  dz(kpe) = (z(kpe) - z(kpe-1))*.5
 265  DO k=kps+1,kpe-1
 266 !  dz(k) = (z(k+1)+z(k))/2. - (z(k)+z(k-1))/2.
 267    dz(k) = (z(k+1) - z(k-1))*.5
 268  ENDDO
 269
 270  ! Normalize
 271  cuml_f_dist = DOT_PRODUCT(dz(k_min:k_max),f(k_min:k_max))
 272  f(k_min:k_max) = f(k_min:k_max)*dz(k_min:k_max)/cuml_f_dist
 273
 274  END SUBROUTINE bellcurve
 275
 276 !************************************************************************
 277 ! This subroutine finds the vertical indices (phy grid) of the follow
 278 ! within a column:
 279 ! 1) ktop - reflectivity cloud top
 280 ! 2) kbtm - reflectivity cloud bottom
 281 ! 3) kupper - upper mode isotherm
 282 ! 3) klower - lower mode isotherm
 283 !************************************************************************
 284  SUBROUTINE kfind ( &
 285               ! Prognostics
 286                 cellcount, t,                         &
 287               ! Namelist settings
 288                 ltng_temp_upper,ltng_temp_lower,      &
 289                 cellcount_method,                     &
 290               ! Order dependent args for domain, mem, and tile dims
 291                 ips, ipe, jps, jpe, kps, kpe,          &
 292               ! Outputs
 293                 ktop,kbtm,kupper,klower               &
 294             )
 295 !-----------------------------------------------------------------
 296 ! Framework
 297  USE module_state_description
 298
 299 ! Model layer
 300  USE module_model_constants
 301
 302  USE module_dm, only: wrf_dm_max_real, wrf_dm_min_real, wrf_dm_sum_real
 303
 304  IMPLICIT NONE
 305 !-----------------------------------------------------------------
 306
 307 ! Prognostics
 308  REAL, DIMENSION( kps:kpe ), INTENT(IN   ) :: cellcount
 309  REAL, DIMENSION( kps:kpe ), INTENT(IN   ) :: t
 310
 311 ! Namelist settings
 312  REAL,    INTENT(IN   )    ::       ltng_temp_lower, ltng_temp_upper
 313  INTEGER, INTENT(IN   )    ::       cellcount_method
 314
 315 ! Order dependent args for domain, mem, and tile (patch) dims
 316  INTEGER, INTENT(IN   )    ::       ips,ipe, jps,jpe, kps,kpe
 317
 318 ! Outputs
 319  INTEGER, INTENT(  OUT)    ::       ktop,kbtm,kupper,klower
 320
 321 ! Local vars
 322  CHARACTER (LEN=250) :: message
 323  REAL    :: ktop_r, kbtm_r, kupper_r, klower_r
 324  INTEGER :: k
 325
 326 !-----------------------------------------------------------------
 327  ktop = kps
 328  kbtm = kps
 329  kupper = kps
 330  klower = kps
 331
 332  ! look for ktop
 333  k = kpe
 334  DO WHILE ( cellcount(k) .eq. 0 .and. k .gt. kps)
 335   k = k-1
 336  ENDDO
 337  ktop = k
 338
 339  ! Look for kbtm
 340  k = kps
 341  DO WHILE( cellcount(k) .eq. 0 .and. k .le. ktop )
 342   k = k+1
 343  ENDDO
 344  kbtm = k
 345 ! if (kbtm .eq. kps) kbtm = kpe
 346
 347  ! Look for kupper
 348  k = kps
 349  DO WHILE ( t(k) .gt. ltng_temp_lower + 273.15 .and. k .lt. kpe )
 350    k = k + 1
 351  ENDDO
 352  klower = k
 353
 354  DO WHILE ( t(k) .gt. ltng_temp_upper + 273.15 .and. k .lt. kpe )
 355    k = k + 1
 356  ENDDO
 357  kupper = k
 358
 359  WRITE(message, * ) ' LNOx_driver: kbtm, ktop, klower, kupper = ', kbtm, ktop, klower, kupper
 360  CALL wrf_debug ( 100, message )
 361
 362  IF ( cellcount_method .eq. 2 ) THEN
 363    kbtm_r = real(kbtm)
 364    ktop_r = real(ktop)
 365    klower_r = real(klower)
 366    kupper_r = real(kupper)
 367    kbtm = nint(wrf_dm_min_real(kbtm_r))
 368    ktop = nint(wrf_dm_max_real(ktop_r))
 369    klower = nint(wrf_dm_max_real(klower_r))
 370    kupper = nint(wrf_dm_max_real(kupper_r))
 371
 372    WRITE(message, * ) ' lightning_driver: kbtm, ktop, klower, kupper = ', kbtm, ktop, klower, kupper
 373    CALL wrf_debug ( 100, message )
 374  endif
 375
 376  END SUBROUTINE kfind
 377
 378 !************************************************************************
 379 ! This function averages out a 3D array into 1D in the 2nd dimension
 380 !************************************************************************
 381  SUBROUTINE horizontalAverage( array3D, ips, ipe, kps, kpe, jps, jpe, array1D )
 382 !-----------------------------------------------------------------
 383  IMPLICIT NONE
 384  REAL, DIMENSION(ips:ipe,kps:kpe,jps:jpe), INTENT(IN) :: array3D
 385  INTEGER, INTENT(IN) :: ips,ipe,kps,kpe,jps,jpe
 386
 387  INTEGER :: k
 388  REAL, DIMENSION(kps:kpe), INTENT(OUT) :: array1D
 389 !-----------------------------------------------------------------
 390  DO k=kps,kpe
 391    array1D(k) = sum(array3D(ips:ipe,k,jps:jpe))/((ipe-ips+1)*(jpe-jps+1))
 392  ENDDO
 393
 394  END SUBROUTINE horizontalAverage
 395
 396  END MODULE module_lightning_nox_decaria