Line data Source code
1 : module mo_lightning
2 : !----------------------------------------------------------------------
3 : ! ... the lightning module
4 : !----------------------------------------------------------------------
5 :
6 : use shr_kind_mod, only : r8 => shr_kind_r8
7 : use ppgrid, only : begchunk, endchunk, pcols, pver
8 : use phys_grid, only : ngcols_p => num_global_phys_cols
9 : use cam_abortutils, only : endrun
10 : use cam_logfile, only : iulog
11 : use spmd_utils, only : masterproc, mpicom
12 :
13 : use physics_buffer, only : pbuf_get_index, physics_buffer_desc, pbuf_get_field, pbuf_get_chunk
14 : use physics_buffer, only : pbuf_add_field, pbuf_set_field, dtype_r8
15 :
16 : implicit none
17 :
18 : private
19 :
20 : public :: lightning_readnl
21 : public :: lightning_register
22 : public :: lightning_init
23 : public :: lightning_no_prod
24 : public :: prod_no
25 :
26 : real(r8),protected, allocatable :: prod_no(:,:,:)
27 :
28 : real(r8) :: factor = 0.1_r8 ! user-controlled scaling factor to achieve arbitrary no prod.
29 : real(r8) :: geo_factor = -huge(1._r8) ! grid cell area factor
30 : real(r8), allocatable :: vdist(:,:) ! vertical distribution of lightning
31 :
32 : logical :: calc_nox_prod = .false.
33 : logical :: calc_lightning = .false.
34 :
35 : integer :: flsh_frq_ndx = -1
36 : integer :: cldtop_ndx = -1, cldbot_ndx = -1
37 :
38 : ! namelist parameter
39 : real(r8) :: lght_no_prd_factor = -huge(1._r8)
40 :
41 : contains
42 :
43 : !-------------------------------------------------------------------------
44 : ! Read namelist options
45 : !-------------------------------------------------------------------------
46 370944 : subroutine lightning_readnl(nlfile)
47 : use namelist_utils, only : find_group_name
48 : use spmd_utils, only : mpicom, masterprocid, mpi_real8, mpi_success
49 :
50 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
51 :
52 : integer :: unitn, ierr
53 : character(len=*), parameter :: subname = 'lightning_readnl'
54 :
55 : ! ===================
56 : ! Namelist definition
57 : ! ===================
58 : namelist /lightning_nl/ lght_no_prd_factor
59 :
60 : ! =============
61 : ! Read namelist
62 : ! =============
63 1536 : if (masterproc) then
64 2 : open( newunit=unitn, file=trim(nlfile), status='old' )
65 2 : call find_group_name(unitn, 'lightning_nl', status=ierr)
66 2 : if (ierr == 0) then
67 0 : read(unitn, lightning_nl, iostat=ierr)
68 0 : if (ierr /= 0) then
69 0 : call endrun(subname // ':: ERROR reading namelist')
70 : end if
71 : end if
72 2 : close(unitn)
73 : end if
74 :
75 : ! ============================
76 : ! Broadcast namelist variables
77 : ! ============================
78 1536 : call mpi_bcast(lght_no_prd_factor, 1, mpi_real8, masterprocid, mpicom, ierr)
79 1536 : if (ierr/=mpi_success) then
80 0 : call endrun(subname//': MPI_BCAST ERROR: lght_no_prd_factor')
81 : end if
82 :
83 1536 : if (masterproc) then
84 2 : write(iulog,*) subname,' lght_no_prd_factor: ',lght_no_prd_factor
85 : end if
86 :
87 1536 : if( lght_no_prd_factor /= 1._r8 ) then
88 1536 : factor = factor*lght_no_prd_factor
89 : end if
90 :
91 1536 : end subroutine lightning_readnl
92 :
93 : !-------------------------------------------------------------------------
94 : ! register phys buffer field for cloud to ground lightning flash frequency
95 : ! to pass to the mediator for land model
96 : !-------------------------------------------------------------------------
97 1536 : subroutine lightning_register()
98 :
99 1536 : call pbuf_add_field('LGHT_FLASH_FREQ','global',dtype_r8,(/pcols/),flsh_frq_ndx) ! per minute
100 :
101 1536 : end subroutine lightning_register
102 :
103 : !-------------------------------------------------------------------------
104 : !-------------------------------------------------------------------------
105 1536 : subroutine lightning_init( pbuf2d )
106 : !----------------------------------------------------------------------
107 : ! ... initialize the lightning module
108 : !----------------------------------------------------------------------
109 : use mo_constants, only : pi
110 :
111 : use cam_history, only : addfld, add_default, horiz_only
112 : use phys_control, only : phys_getopts
113 : use time_manager, only : is_first_step
114 :
115 : !----------------------------------------------------------------------
116 : ! ... dummy args
117 : !----------------------------------------------------------------------
118 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
119 :
120 : !----------------------------------------------------------------------
121 : ! ... local variables
122 : !----------------------------------------------------------------------
123 : integer :: astat, err
124 : logical :: history_cesm_forcing
125 : character(len=*),parameter :: prefix = 'lightning_init: '
126 :
127 1536 : cldtop_ndx = pbuf_get_index('CLDTOP',errcode=err)
128 1536 : cldbot_ndx = pbuf_get_index('CLDBOT',errcode=err)
129 1536 : calc_lightning = cldtop_ndx>0 .and. cldbot_ndx>0
130 :
131 1536 : if (.not.calc_lightning) return
132 :
133 1536 : calc_nox_prod = lght_no_prd_factor>0._r8
134 :
135 1536 : if (calc_nox_prod) then
136 :
137 0 : if (masterproc) write(iulog,*) prefix,'lightning no production scaling factor = ',factor
138 :
139 : !----------------------------------------------------------------------
140 : ! ... vdist(kk,itype) = % of lightning nox between (kk-1) and (kk)
141 : ! km for profile itype
142 : !----------------------------------------------------------------------
143 0 : allocate(vdist(16,3),stat=astat)
144 0 : if( astat /= 0 ) then
145 0 : write(iulog,*) prefix,'failed to allocate vdist; error = ',astat
146 0 : call endrun(prefix//'failed to allocate vdist')
147 : end if
148 0 : vdist(:,1) = (/ 3.0_r8, 3.0_r8, 3.0_r8, 3.0_r8, 3.4_r8, 3.5_r8, 3.6_r8, 4.0_r8, & ! midlat cont
149 0 : 5.0_r8, 7.0_r8, 9.0_r8, 14.0_r8, 16.0_r8, 14.0_r8, 8.0_r8, 0.5_r8 /)
150 0 : vdist(:,2) = (/ 2.5_r8, 2.5_r8, 2.5_r8, 2.5_r8, 2.5_r8, 2.5_r8, 2.5_r8, 6.1_r8, & ! trop marine
151 0 : 17.0_r8, 15.4_r8, 14.5_r8, 13.0_r8, 12.5_r8, 2.8_r8, 0.9_r8, 0.3_r8 /)
152 0 : vdist(:,3) = (/ 2.0_r8, 2.0_r8, 2.0_r8, 1.5_r8, 1.5_r8, 1.5_r8, 3.0_r8, 5.8_r8, & ! trop cont
153 0 : 7.6_r8, 9.6_r8, 11.0_r8, 14.0_r8, 14.0_r8, 14.0_r8, 8.2_r8, 2.3_r8 /)
154 :
155 0 : allocate( prod_no(pcols,pver,begchunk:endchunk),stat=astat )
156 0 : if( astat /= 0 ) then
157 0 : write(iulog,*) prefix, 'failed to allocate prod_no; error = ',astat
158 0 : call endrun(prefix//'failed to allocate prod_no')
159 : end if
160 0 : geo_factor = ngcols_p/(4._r8*pi)
161 :
162 0 : call addfld( 'LNO_COL_PROD', horiz_only, 'I', 'Tg N yr-1', 'lightning column NO source' )
163 0 : call addfld( 'LNO_PROD', (/ 'lev' /), 'I', 'molecules/cm3/s', 'lightning insitu NO source' )
164 0 : call addfld( 'FLASHENGY', horiz_only, 'I', 'J', 'lightning flash energy' ) ! flash energy
165 :
166 0 : call phys_getopts( history_cesm_forcing_out = history_cesm_forcing )
167 0 : if ( history_cesm_forcing ) then
168 0 : call add_default('LNO_COL_PROD',1,' ')
169 : endif
170 :
171 0 : if (is_first_step()) then
172 0 : call pbuf_set_field(pbuf2d, flsh_frq_ndx, 0.0_r8)
173 : endif
174 :
175 : endif
176 :
177 1536 : call addfld( 'FLASHFRQ', horiz_only, 'I', 'min-1', 'lightning flash rate' ) ! flash frequency in grid box per minute (PPP)
178 1536 : call addfld( 'CLDHGT', horiz_only, 'I', 'km', 'cloud top height' ) ! cloud top height
179 1536 : call addfld( 'DCHGZONE', horiz_only, 'I', 'km', 'depth of discharge zone' ) ! depth of discharge zone
180 1536 : call addfld( 'CGIC', horiz_only, 'I', '1', 'ratio of cloud-ground/intracloud discharges' ) ! ratio of cloud-ground/intracloud discharges
181 1536 : call addfld( 'LGHTNG_CLD2GRND', horiz_only, 'I', 'min-1', 'clound-to-ground lightning flash rate') ! clound to ground flash frequency
182 :
183 1536 : end subroutine lightning_init
184 :
185 : !-------------------------------------------------------------------------
186 : !-------------------------------------------------------------------------
187 369408 : subroutine lightning_no_prod( state, pbuf2d, cam_in )
188 : !----------------------------------------------------------------------
189 : ! ... set no production from lightning
190 : !----------------------------------------------------------------------
191 1536 : use physics_types, only : physics_state
192 : use physconst, only : rga
193 : use phys_grid, only : get_rlat_all_p, get_wght_all_p
194 : use cam_history, only : outfld
195 : use camsrfexch, only : cam_in_t
196 : use shr_reprosum_mod, only : shr_reprosum_calc
197 : use mo_constants, only : rearth, d2r
198 :
199 : !----------------------------------------------------------------------
200 : ! ... dummy args
201 : !----------------------------------------------------------------------
202 : type(physics_state), intent(in) :: state(begchunk:endchunk) ! physics state
203 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
204 : type(cam_in_t), intent(in) :: cam_in(begchunk:endchunk) ! physics state
205 :
206 : !----------------------------------------------------------------------
207 : ! ... local variables
208 : !----------------------------------------------------------------------
209 : real(r8), parameter :: land = 1._r8
210 : real(r8), parameter :: secpyr = 365._r8 * 8.64e4_r8
211 :
212 : integer :: cldtind ! level index for cloud top
213 : integer :: cldbind ! level index for cloud base > 273k
214 : integer :: k, kk, zlow_ind, zhigh_ind, itype
215 : real(r8) :: glob_flashfreq ! global flash frequency [s-1]
216 : real(r8) :: glob_noprod ! global rate of no production [as tgn/yr]
217 : real(r8) :: frac_sum ! work variable
218 : real(r8) :: zlow
219 : real(r8) :: zhigh
220 : real(r8) :: zlow_scal
221 : real(r8) :: zhigh_scal
222 : real(r8) :: fraction
223 : real(r8) :: dchgz
224 738816 : real(r8) :: dchgzone(pcols,begchunk:endchunk) ! depth of discharge zone [km]
225 738816 : real(r8) :: cldhgt(pcols,begchunk:endchunk) ! cloud top height [km]
226 738816 : real(r8) :: cgic(pcols,begchunk:endchunk) ! cloud-ground/intracloud discharge ratio
227 738816 : real(r8) :: flash_energy(pcols,begchunk:endchunk) ! energy of flashes per second
228 738816 : real(r8) :: prod_no_col(pcols,begchunk:endchunk) ! global no production rate for diagnostics
229 : real(r8) :: wrk, wrk1, wrk2(1)
230 : integer :: icol ! column index
231 : integer :: ncol ! columns per chunk
232 : integer :: lchnk ! chunk index
233 369408 : real(r8),pointer :: cldtop(:) ! cloud top level index
234 369408 : real(r8),pointer :: cldbot(:) ! cloud bottom level index
235 : real(r8) :: zmid(pcols,pver) ! geopot height above surface at midpoints (m)
236 738816 : real(r8) :: zint(pcols,pver+1,begchunk:endchunk) ! geopot height above surface at interfaces (m)
237 : real(r8) :: zsurf(pcols) ! geopot height above surface at interfaces (m)
238 : real(r8) :: rlats(pcols) ! column latitudes in chunks
239 : real(r8) :: wght(pcols)
240 :
241 738816 : real(r8) :: glob_prod_no_col(pcols,begchunk:endchunk)
242 738816 : real(r8) :: flash_freq(pcols,begchunk:endchunk)
243 :
244 : !----------------------------------------------------------------------
245 : ! ... parameters to determine cg/ic ratio [price and rind, 1993]
246 : !----------------------------------------------------------------------
247 : real(r8), parameter :: ca = .021_r8
248 : real(r8), parameter :: cb = -.648_r8
249 : real(r8), parameter :: cc = 7.49_r8
250 : real(r8), parameter :: cd = -36.54_r8
251 : real(r8), parameter :: ce = 64.09_r8
252 : real(r8), parameter :: t0 = 273._r8
253 : real(r8), parameter :: m2km = 1.e-3_r8
254 : real(r8), parameter :: km2cm = 1.e5_r8
255 : real(r8), parameter :: lat25 = 25._r8*d2r ! 25 degrees latitude in radians
256 :
257 : real(r8) :: flash_freq_land, flash_freq_ocn
258 369408 : real(r8), pointer :: cld2grnd_flash_freq(:)
259 :
260 369408 : if (.not.calc_lightning) return
261 :
262 369408 : nullify(cld2grnd_flash_freq)
263 :
264 : !----------------------------------------------------------------------
265 : ! ... initialization
266 : !----------------------------------------------------------------------
267 :
268 25685400 : flash_freq(:,:) = 0._r8
269 25685400 : cldhgt(:,:) = 0._r8
270 25685400 : dchgzone(:,:) = 0._r8
271 25685400 : cgic(:,:) = 0._r8
272 25685400 : flash_energy(:,:) = 0._r8
273 :
274 369408 : if (calc_nox_prod) then
275 0 : prod_no(:,:,:) = 0._r8
276 0 : prod_no_col(:,:) = 0._r8
277 0 : glob_prod_no_col(:,:) = 0._r8
278 : end if
279 :
280 : !--------------------------------------------------------------------------------
281 : ! ... estimate flash frequency and resulting no emissions
282 : ! [price, penner, prather, 1997 (jgr)]
283 : ! lightning only occurs in convective clouds with a discharge zone, i.e.
284 : ! an altitude range where liquid water, ice crystals, and graupel coexist.
285 : ! we test this by examining the temperature at the cloud base.
286 : ! it is assumed that only one thunderstorm occurs per grid box, and its
287 : ! flash frequency is determined by the maximum cloud top height (not the
288 : ! depth of the discharge zone). this is somewhat speculative but yields
289 : ! reasonable results.
290 : !
291 : ! the cg/ic ratio is determined by an empirical formula from price and
292 : ! rind [1993]. the average energy of a cg flash is estimated as 6.7e9 j,
293 : ! and the average energy of a ic flash is assumed to be 1/10 of that value.
294 : ! the no production rate is assumed proportional to the discharge energy
295 : ! with 1e17 n atoms per j. the total number of n atoms is then distributed
296 : ! over the complete column of grid boxes.
297 : !--------------------------------------------------------------------------------
298 1858584 : Chunk_loop : do lchnk = begchunk,endchunk
299 1489176 : ncol = state(lchnk)%ncol
300 1489176 : call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk), flsh_frq_ndx, cld2grnd_flash_freq )
301 1489176 : call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk), cldtop_ndx, cldtop )
302 1489176 : call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk), cldbot_ndx, cldbot )
303 24865776 : zsurf(:ncol) = state(lchnk)%phis(:ncol)*rga
304 1489176 : call get_wght_all_p(lchnk, pcols, wght)
305 :
306 139982544 : do k = 1,pver
307 2312517168 : zmid(:ncol,k) = state(lchnk)%zm(:ncol,k) + zsurf(:ncol)
308 2314006344 : zint(:ncol,k,lchnk) = state(lchnk)%zi(:ncol,k) + zsurf(:ncol)
309 : end do
310 24865776 : zint(:ncol,pver+1,lchnk) = state(lchnk)%zi(:ncol,pver+1) + zsurf(:ncol)
311 :
312 25315992 : cld2grnd_flash_freq(:) = 0.0_r8
313 :
314 24865776 : col_loop : do icol = 1,ncol
315 : !--------------------------------------------------------------------------------
316 : ! ... find cloud top and bottom level above 273k
317 : !--------------------------------------------------------------------------------
318 23376600 : cldtind = nint( cldtop(icol) )
319 23376600 : cldbind = nint( cldbot(icol) )
320 54513377 : do
321 77889977 : if( cldbind <= cldtind .or. state(lchnk)%t(icol,cldbind) < t0 ) then
322 : exit
323 : end if
324 54513377 : cldbind = cldbind - 1
325 : end do
326 24865776 : cloud_layer : if( cldtind < pver .and. cldtind > 0 .and. cldtind < cldbind ) then
327 : !--------------------------------------------------------------------------------
328 : ! ... compute cloud top height and depth of charging zone
329 : !--------------------------------------------------------------------------------
330 1653102 : cldhgt(icol,lchnk) = m2km*max( 0._r8,zint(icol,cldtind,lchnk) )
331 1653102 : dchgz = cldhgt(icol,lchnk) - m2km*zmid(icol,cldbind)
332 1653102 : dchgzone(icol,lchnk) = dchgz
333 : !--------------------------------------------------------------------------------
334 : ! ... compute flash frequency for given cloud top height
335 : ! (flashes storm^-1 min^-1)
336 : !--------------------------------------------------------------------------------
337 1653102 : flash_freq_land = 3.44e-5_r8 * cldhgt(icol,lchnk)**4.9_r8
338 1653102 : flash_freq_ocn = 6.40e-4_r8 * cldhgt(icol,lchnk)**1.7_r8
339 : flash_freq(icol,lchnk) = cam_in(lchnk)%landfrac(icol)*flash_freq_land + &
340 1653102 : cam_in(lchnk)%ocnfrac(icol) *flash_freq_ocn
341 :
342 : !--------------------------------------------------------------------------------
343 : ! cgic = proportion of cloud-to-ground flashes
344 : ! NOx from lightning 1. Global distribution based on lightning physics, C Price et al
345 : ! JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 102, NO. D5, PAGES 5929-5941, MARCH 20, 1997
346 : ! (https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/96JD03504)
347 : ! eq 14
348 : !--------------------------------------------------------------------------------
349 1653102 : cgic(icol,lchnk) = 1._r8/((((ca*dchgz + cb)*dchgz + cc) *dchgz + cd)*dchgz + ce)
350 1653102 : if( dchgz < 5.5_r8 ) then
351 1461991 : cgic(icol,lchnk) = 0._r8
352 191111 : else if( dchgz > 14._r8 ) then
353 0 : cgic(icol,lchnk) = .02_r8
354 : end if
355 :
356 1653102 : cld2grnd_flash_freq(icol) = cam_in(lchnk)%landfrac(icol)*flash_freq_land*cgic(icol,lchnk) ! cld-to-grnd flash frq (per min)
357 :
358 1653102 : if (calc_nox_prod) then
359 : !--------------------------------------------------------------------------------
360 : ! ... compute flash energy (cg*6.7e9 + ic*6.7e8)
361 : ! and convert to total energy per second
362 : ! set ic = cg
363 : !--------------------------------------------------------------------------------
364 0 : flash_energy(icol,lchnk) = 6.7e9_r8 * flash_freq(icol,lchnk)/60._r8
365 : !--------------------------------------------------------------------------------
366 : ! ... LKE Aug 23, 2005: scale production to account for different grid
367 : ! box sizes. This requires a reduction in the overall fudge factor
368 : ! (e.g., from 1.2 to 0.5)
369 : !--------------------------------------------------------------------------------
370 0 : flash_energy(icol,lchnk) = flash_energy(icol,lchnk) * wght(icol) * geo_factor
371 : !--------------------------------------------------------------------------------
372 : ! ... compute number of n atoms produced per second
373 : ! and convert to n atoms per second per cm2 and apply fudge factor
374 : !--------------------------------------------------------------------------------
375 0 : prod_no_col(icol,lchnk) = 1.e17_r8*flash_energy(icol,lchnk)/(1.e4_r8*rearth*rearth*wght(icol)) * factor
376 :
377 : !--------------------------------------------------------------------------------
378 : ! ... compute global no production rate in tgn/yr:
379 : ! tgn per second: * 14.00674 * 1.65979e-24 * 1.e-12
380 : ! nb: 1.65979e-24 = 1/avo
381 : ! tgn per year: * secpyr
382 : !--------------------------------------------------------------------------------
383 : glob_prod_no_col(icol,lchnk) = 1.e17_r8*flash_energy(icol,lchnk) &
384 0 : * 14.00674_r8 * 1.65979e-24_r8 * 1.e-12_r8 * secpyr * factor
385 : end if
386 : end if cloud_layer
387 : end do Col_loop
388 :
389 1858584 : call outfld( 'LGHTNG_CLD2GRND', cld2grnd_flash_freq, pcols, lchnk )
390 : end do Chunk_loop
391 :
392 1858584 : do lchnk = begchunk,endchunk
393 1489176 : call outfld( 'FLASHFRQ', flash_freq(:,lchnk), pcols, lchnk )
394 1489176 : call outfld( 'CGIC', cgic(:,lchnk), pcols, lchnk )
395 1489176 : call outfld( 'CLDHGT', cldhgt(:,lchnk), pcols, lchnk )
396 1858584 : call outfld( 'DCHGZONE', dchgzone(:,lchnk), pcols, lchnk )
397 : enddo
398 :
399 369408 : if (.not.calc_nox_prod) return
400 :
401 : !--------------------------------------------------------------------------------
402 : ! ... Accumulate global total, convert to flashes per second
403 : ! ... Accumulate global NO production rate
404 : !--------------------------------------------------------------------------------
405 0 : kk = pcols*(endchunk-begchunk+1)
406 0 : call shr_reprosum_calc( flash_freq, wrk2,kk,kk,1, commid=mpicom)
407 0 : glob_flashfreq=wrk2(1)/60._r8
408 0 : call shr_reprosum_calc( glob_prod_no_col, wrk2,kk,kk,1, commid=mpicom)
409 0 : glob_noprod = wrk2(1)
410 0 : if( masterproc ) then
411 0 : write(iulog,*) ' '
412 : write(iulog,'(''Global flash freq (/s), lightning NOx (TgN/y) = '',2f10.4)') &
413 0 : glob_flashfreq, glob_noprod
414 : end if
415 :
416 0 : if( glob_noprod > 0._r8 ) then
417 : !--------------------------------------------------------------------------------
418 : ! ... Distribute production up to cloud top [Pickering et al., 1998 (JGR)]
419 : !--------------------------------------------------------------------------------
420 0 : do lchnk = begchunk,endchunk
421 0 : call get_rlat_all_p(lchnk, pcols, rlats)
422 0 : ncol = state(lchnk)%ncol
423 0 : call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk), cldtop_ndx, cldtop )
424 0 : do icol = 1,ncol
425 0 : cldtind = nint( cldtop(icol) )
426 0 : if( prod_no_col(icol,lchnk) > 0._r8 ) then
427 0 : if( cldhgt(icol,lchnk) > 0._r8 ) then
428 0 : if( abs( rlats(icol) ) > lat25 ) then
429 : itype = 1 ! midlatitude continental
430 0 : else if( nint( cam_in(lchnk)%landfrac(icol) ) == land ) then
431 : itype = 3 ! tropical continental
432 : else
433 0 : itype = 2 ! topical marine
434 : end if
435 0 : frac_sum = 0._r8
436 0 : do k = cldtind,pver
437 0 : zlow = zint(icol,k+1,lchnk) * m2km ! lower interface height (km)
438 0 : zlow_scal = zlow * 16._r8/cldhgt(icol,lchnk) ! scale to 16 km convection height
439 0 : zlow_ind = max( 1,INT(zlow_scal)+1 ) ! lowest vdist index to include in layer
440 0 : zhigh = zint(icol,k,lchnk) * m2km ! upper interface height (km)
441 0 : zhigh_scal = zhigh * 16._r8/cldhgt(icol,lchnk) ! height (km) scaled to 16km convection height
442 0 : zhigh_ind = max( 1,MIN( 16,INT(zhigh_scal)+1 ) ) ! highest vdist index to include in layer
443 0 : do kk = zlow_ind,zhigh_ind
444 0 : wrk = kk
445 0 : wrk1 = kk - 1
446 : fraction = min( zhigh_scal,wrk ) & ! fraction of vdist in this model layer
447 0 : - max( zlow_scal,wrk1 )
448 0 : fraction = max( 0._r8, min( 1._r8,fraction ) )
449 0 : frac_sum = frac_sum + fraction*vdist(kk,itype)
450 0 : prod_no(icol,k,lchnk) = prod_no(icol,k,lchnk) & ! sum the fraction of column NOx in layer k
451 0 : + fraction*vdist(kk,itype)*.01_r8
452 : end do
453 0 : prod_no(icol,k,lchnk) = prod_no_col(icol,lchnk) * prod_no(icol,k,lchnk) & ! multiply fraction by column amount
454 0 : / (km2cm*(zhigh - zlow)) ! and convert to atom N cm^-3 s^-1
455 : end do
456 : end if
457 : end if
458 : end do
459 : end do
460 : end if
461 :
462 : !--------------------------------------------------------------------------------
463 : ! ... output lightning no production to history file
464 : !--------------------------------------------------------------------------------
465 0 : do lchnk = begchunk,endchunk
466 0 : call outfld( 'LNO_PROD', prod_no(:,:,lchnk), pcols, lchnk )
467 0 : call outfld( 'LNO_COL_PROD', glob_prod_no_col(:,lchnk), pcols, lchnk )
468 0 : call outfld( 'FLASHENGY', flash_energy(:,lchnk), pcols, lchnk )
469 : enddo
470 :
471 738816 : end subroutine lightning_no_prod
472 :
473 : end module mo_lightning
|
