LCOV - code coverage report
Current view: | top level - chemistry/mozart - mo_lightning.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 90 | 167 | 53.9 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 4 | 4 | 100.0 %

              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    40207752 :        do k = 1,pver
         307   646510176 :           zmid(:ncol,k)   = state(lchnk)%zm(:ncol,k) + zsurf(:ncol)
         308   647999352 :           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    49139217 :           do
         321    72515817 :              if( cldbind <= cldtind .or. state(lchnk)%t(icol,cldbind) < t0 ) then
         322             :                 exit
         323             :              end if
         324    49139217 :              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    11125221 :              cldhgt(icol,lchnk)   = m2km*max( 0._r8,zint(icol,cldtind,lchnk) )
         331    11125221 :              dchgz = cldhgt(icol,lchnk) - m2km*zmid(icol,cldbind)
         332    11125221 :              dchgzone(icol,lchnk) = dchgz
         333             :              !--------------------------------------------------------------------------------
         334             :              !       ... compute flash frequency for given cloud top height
         335             :              !           (flashes storm^-1 min^-1)
         336             :              !--------------------------------------------------------------------------------
         337    11125221 :              flash_freq_land = 3.44e-5_r8 * cldhgt(icol,lchnk)**4.9_r8
         338    11125221 :              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    11125221 :                                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    11125221 :              cgic(icol,lchnk) = 1._r8/((((ca*dchgz + cb)*dchgz + cc) *dchgz + cd)*dchgz + ce)
         350    11125221 :              if( dchgz < 5.5_r8 ) then
         351     9180248 :                 cgic(icol,lchnk) = 0._r8
         352     1944973 :              else if( dchgz > 14._r8 ) then
         353          76 :                 cgic(icol,lchnk) = .02_r8
         354             :              end if
         355             :
         356    11125221 :              cld2grnd_flash_freq(icol) = cam_in(lchnk)%landfrac(icol)*flash_freq_land*cgic(icol,lchnk) ! cld-to-grnd flash frq (per min)
         357             :
         358    11125221 :              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

Generated by: LCOV version 1.14