ecearth3/sources/tm5mp/proj/cb05/emission_nox.F90

!
#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      EMISSION_NOX
!
! !DESCRIPTION: hold data and methods for NOx emissions. 
!               -----------------
!               AR5 emissions
!               -----------------
!               For each month, arrays emis_nox_2d/3d have to be filled.
!               It follows the following settins:
!               - take emiss_ene/dom/ind/wst/agr/awb/slv/tra/shp/
!                 /air/forestfire/grassfire from AR5 data sets (NO GFED!!) 
!               - use natural emissions from MACC data sets 
!                 (emiss_nat/soil/bio/oc)
!               - vertical distribution is done via emission_vdist_by_sector 
!                 (emission_data.F90)
!               - lightning is done online (eminox_lightning)
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_NOX
  !
  ! !USES:
  !
  use GO,              only : gol, goErr, goPr, goBug
  use tm5_distgrid,    only : dgrid, get_distgrid, scatter, gather
  use dims,            only : nregions, idate, dy, okdebug
  use global_types,    only : emis_data, d3_data,d23_data
  use emission_read,   only : used_providers, has_emis

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public   ::  Emission_NOx_Init
  public   ::  Emission_NOx_Done
  public   ::  Emission_NOx_Declare
  public   ::  Emission_NOx_bb_daily_cycle
#ifndef without_convection
  public   ::  lightningNOX
#endif
  public   ::  nox_emis_3d, nox_emis_3d_bb_app
  public   ::  eminox_lightning
  public   ::  flashrate_lightning
  !
  ! !DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'emission_nox'

  type(d3_data), dimension(nregions), target  :: nox_emis_3d, nox_emis_3d_bb, nox_emis_3d_bb_app
  type(d3_data), dimension(nregions), target  :: eminox_lightning
  type(d23_data), dimension(nregions), target  :: flashrate_lightning

  integer :: nox_2dsec, nox_3dsec
  real    :: fscalelig   ! scaling used in lightning NOX production to get 5.98 Tg for 2006
  !
  ! !REVISION HISTORY: 
  !    1 Oct 2010 - Achim Strunk  - overhaul for AR5
  !    1 Dec 2011 - Narcisa Banda - added EDGAR 4
  !   27 Jun 2012 - P. Le Sager   - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !    NOx emissions are added directly in chemistry, instead of apart from it. 
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_NOX_INIT
  !
  ! !DESCRIPTION: allocate memory
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Emission_NOx_Init( status )
    !
    ! !USES:
    !
    use dims,          only : lm
    use emission_read, only : providers_def, numb_providers, ed42_nsect_nox
    use emission_data, only : LAR5BMB
    use emission_read, only : n_ar5_ant_sec, n_ar5_shp_sec, n_ar5_air_sec, n_ar5_bmb_sec
    use emission_read, only : ar5_cat_ant,   ar5_cat_shp,   ar5_cat_air,   ar5_cat_bmb
#ifndef without_convection
    use meteodata,     only : set, gph_dat, temper_dat, cp_dat
    use emission_data, only : use_tiedkte
#endif 
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - adapted for AR5
    !   27 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !   10 Jul 2013 - Ph. Le Sager - init lightning when no inventory is selected
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/Emission_NOx_Init'
    integer                     ::  region, i1, j1, i2, j2
    integer                     ::  imr, jmr, lmr, lsec, lprov

    ! --- begin --------------------------------------
    
    status = 0
    
#ifndef without_convection
    ! Meteo used for LightningNOx
    do region=1,nregions
       call Set( temper_dat(region), status, used=.true. )
       call Set(    gph_dat(region), status, used=.true. )
       call Set(     cp_dat(region), status, used=.true. )
    enddo
    !
    ! Scaling parameter for LiNOx
    !
    ! Set to get 5.98 Tg N with 2006 EI met fields. This is resolution
    ! and met fields dependent. The factor has been estimated for:
    !  - @3x2 and 34 levels, Tiedkte : fscalelig=13.715
    !  - @1x1 and 34 levels, Tiedkte : fscalelig=17.051
    !  - @3x2 and 34 levels, EI conv : fscalelig=13.715*0.786
    !  - @1x1 and 34 levels, EI conv : fscalelig=17.051*0.649
    !
    if (use_tiedkte) then ! convective fluxes computed from T/rh/wind (Tiedkte)
       fscalelig=13.715              ! 3x2-34L, Tiedkte scheme
       if (dy == 1) fscalelig=17.051 ! 1x1-34L, Tiedkte scheme
    else
       fscalelig=10.78               ! 3x2-34L, EI convective fluxes
       if (dy == 1) fscalelig=11.066 ! 1x1-34L, EI convective fluxes
    endif

#endif    

    ! allocate information arrays (2d and 3d)
    do region=1,nregions

       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )          
       lmr = lm(region)

       allocate(eminox_lightning(region)%d3(i1:i2,j1:j2,lmr) )
       eminox_lightning(region)%d3=0.
       allocate(flashrate_lightning(region)%d23(i1:i2,j1:j2) )
       flashrate_lightning(region)%d23=0.

    enddo


    ! Check if any inventory is used
    if(.not. has_emis) return
    
    ! nb of sectors
    nox_2dsec = 0
    nox_3dsec = 0
    do lprov = 1, numb_providers
       if (count(used_providers.eq.providers_def(lprov)%name)/=0) then 

          if (trim(providers_def(lprov)%name) .eq. 'AR5') then
             
             ! nb of available sectors in AR5 depends on category
             nox_2dsec = nox_2dsec + n_ar5_ant_sec*count('NO'.eq.ar5_cat_ant) + &
                                     n_ar5_shp_sec*count('NO'.eq.ar5_cat_shp)

             if (LAR5BMB) nox_2dsec = nox_2dsec + n_ar5_bmb_sec*count('NO'.eq.ar5_cat_bmb)

             nox_3dsec = nox_3dsec + n_ar5_air_sec*count('NO'.eq.ar5_cat_air)

!             nox_2dsec = nox_2dsec + providers_def(lprov)%nsect2d
!             nox_3dsec = nox_3dsec + count('NO'.eq.ar5_cat_air)

          elseif (trim(providers_def(lprov)%name) .eq. 'ED42') then

             nox_2dsec = nox_2dsec + ed42_nsect_nox
             ! no 3d sectors in EDGAR 4.2

          else
             nox_2dsec = nox_2dsec + providers_def(lprov)%nsect2d
             nox_3dsec = nox_3dsec + providers_def(lprov)%nsect3d
          endif
       endif
    enddo

    ! allocate information arrays (2d and 3d)
    do region=1,nregions

       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )          
       lmr = lm(region)

       allocate( nox_emis_3d       (region)%d3(i1:i2,j1:j2,lmr) )
       allocate( nox_emis_3d_bb    (region)%d3(i1:i2,j1:j2,lmr) )
       allocate( nox_emis_3d_bb_app(region)%d3(i1:i2,j1:j2,lmr) )
       
    enddo

    status = 0

  END SUBROUTINE EMISSION_NOX_INIT
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_NOX_DONE
  !
  ! !DESCRIPTION: Free memory
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_NOX_DONE( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)           ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - adapted for AR5
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter    :: rname = mname//'/Emission_NOx_Done'
    integer                        :: region, lsec

    ! --- begin ---------------------------------
    
    status = 0
    if(.not. has_emis) return

    do region = 1, nregions
       deallocate( nox_emis_3d       (region)%d3 )
       deallocate( nox_emis_3d_bb    (region)%d3 )
       deallocate( nox_emis_3d_bb_app(region)%d3 )
       deallocate( eminox_lightning  (region)%d3 )
       deallocate( flashrate_lightning  (region)%d23 )
    end do

    status = 0

  END SUBROUTINE EMISSION_NOX_DONE
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_NOX_DECLARE
  !
  ! !DESCRIPTION: Opens, reads and evaluates input files (per month). 
  !               Provides emissions on 2d/3d-arrays which are then added
  !               in the chemistry routine (no *apply !).
  !               Vertically distribute the 2D dataset according to sector.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_NOX_DECLARE( status )
    !
    ! !USES:
    !
    use toolbox,       only : coarsen_emission
    use partools,      only : isRoot, par_broadcast
    use dims,          only : im, jm, lm, nlon360, nlat180, iglbsfc
    use dims,          only : newsrun, idate, sec_month
    use chem_param,    only : xmn, xmno2, xmno
    use emission_data, only : msg_emis, emission_vdist_by_sector, LAR5BMB
    ! ---------------- AR5 - EDGAR 4 - ETC. --------------------
    use emission_data, only : emis_input_year_nox, emis_input_year_nat
    use emission_data, only : emis_input_dir_mac,   emis_input_dir_ed4
    use emission_data, only : emis_input_dir_retro, emis_input_dir_gfed

    use emission_read, only : emission_ar5_regrid_aircraft

    use emission_read, only : emission_cmip6_ReadSector
    use emission_read, only : emission_cmip6bmb_ReadSector
    use emission_read, only : emission_ar5_ReadSector,  emission_macc_ReadSector
    use emission_read, only : emission_ed4_ReadSector,  emission_gfed_ReadSector
    use emission_read, only : emission_retro_ReadSector

    use emission_read, only : sectors_def, numb_sectors
    use emission_read, only : ar5_dim_3ddata
    use emission_read, only : ed42_nox_sectors
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)    ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk  - adapted for AR5
    !    1 Dec 2011 - Narcisa Banda - added EDGAR 4
    !   27 Jun 2012 - Ph. Le Sager  - adapted for lon-lat MPI domain decomposition
    !   25 Feb 2014 - Jason Williams - separate array for BMB so that burning daily cycle can be applied
    !
    ! !REMARKS:
    !   (1) Because we do not use an apply method, the vertical distribution
    !        is done here. However this is a bug, since this is time dependent.
    !        Possible solution: do vert dist in chemistry like the BMB cycle,
    !        or in the more general BMB cycle routine.
    !        
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/emission_nox_declare'

    integer                     :: region, hasData
    integer,parameter           :: add_field=0
    integer,parameter           :: amonth=2
    integer                     :: imr, jmr, lmr, lsec, i1, i2, j1, j2

    ! AR5 temporary arrays 
    real, dimension(:,:,:), allocatable          :: field3d !, field3d2
    type(d3_data),   dimension(nregions), target :: emis3d, work, work3d
    type(emis_data), dimension(nregions), target :: emis2d, wrk2D
    ! defensive
    integer :: seccount2d, seccount3d

    ! --- begin -----------------------------------------
    
    status = 0
    if(.not. has_emis) return

    write(gol,'(" EMISS-INFO ------------- read NOx emissions -------------")'); call goPr

    ! reset emissions, allocate work array
    do region = 1, nregions
       nox_emis_3d(region)%d3   = 0.0 ; nox_emis_3d_bb(region)%d3   = 0.0

       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )          
       lmr = lm(region)

       allocate( work3d(region)%d3  (i1:i2,j1:j2, ar5_dim_3ddata) ) ; work3d(region)%d3   = 0.0
       allocate( emis3d(region)%d3  (i1:i2,j1:j2, lmr           ) ) ; emis3d(region)%d3   = 0.0
       allocate( emis2d(region)%surf(i1:i2,j1:j2)                 ) ; emis2d(region)%surf = 0.0
    end do

    ! global arrays for coarsening    
    do region = 1, nregions
       if (isRoot)then
          allocate(work(region)%d3(im(region),jm(region),ar5_dim_3ddata))
       else
          allocate(work(region)%d3(1,1,1))
       end if
    enddo

    do region = 1, nregions
       wrk2D(region)%surf => work(region)%d3(:,:,1)
    end do

    ! --------------------------------
    ! do a loop over available sectors
    ! --------------------------------
    
    ! count 2d and 3d sectors
    seccount2d = 0
    seccount3d = 0

    ! always allocate here 3d data set (for 2d sectors it will be filled in first layer only)
    if (isRoot) then
       allocate( field3d( nlon360, nlat180, ar5_dim_3ddata ) ) ; field3d = 0.0
    else
       allocate( field3d( 1, 1, 1 ) )
    end if

    sec: do lsec = 1, numb_sectors

       if (count(used_providers.eq.sectors_def(lsec)%prov).eq.0) cycle 

       if ((trim(sectors_def(lsec)%prov).eq.'ED42') .and. (count(ed42_nox_sectors.eq.sectors_def(lsec)%name) .eq. 0)) cycle 

       if (associated(sectors_def(lsec)%species)) then  ! AR5 check
          if (count('NO'.eq.sectors_def(lsec)%species).eq.0) cycle
          if ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) cycle
       endif

       if( sectors_def(lsec)%f3d ) then 
          seccount3d = seccount3d + 1
       else 
          seccount2d = seccount2d + 1
       end if
       
       field3d = 0.0

#ifdef with_online_nox
       ! skip natural nox in case it is calculated online
       if( trim(sectors_def(lsec)%namecat == 'natural') ) then 
          write (gol,'(80("-"))'); call goPr
          write (gol,'("INFO: skipping sector `",a,"` due to `with_online_nox` ")') trim(sectors_def(lsec)%name); call goPr
          cycle
       end if
#endif 

       root: if (isRoot) then ! READ

          select case( trim(sectors_def(lsec)%prov) )

          case( 'CMIP6' )
            call emission_cmip6_ReadSector( 'NOx', emis_input_year_nox, idate(2), lsec, field3d, status )
            IF_NOTOK_RETURN(status=1;deallocate(field3d))

            ! convert from (kg NO2)/s to (kg N)/s
            field3d = field3d * xmn / xmno2                 

          case( 'CMIP6BMB' )
            call emission_cmip6bmb_ReadSector( 'NOx', emis_input_year_nox, idate(2), lsec, field3d, status )
            IF_NOTOK_RETURN(status=1;deallocate(field3d))

            ! convert from (kg NO)/s to (kg N)/s
            ! both for the historical period and the future scenarios
            ! see http://www.falw.vu/~gwerf/GFED/GFED4/ancill/GFED4_Emission_Factors.txt
            ! and the following attribute in the SSP NetCDF files: 
            ! reporting_unit = "Mass flux of NOx, reported as NO"
            field3d = field3d * xmn / xmno

          case( 'AR5' )

             ! AR5 emissions included NO and NO2 aircraft emissions, but they are duplicates. So 
             ! only take into account one of the sets. (in TM5: NO, skip NO2). 

             ! Screen out solvent sector for NO, 
             ! because it is zero in the RCPs
             ! and not present in the historical files.
             if (trim(sectors_def(lsec)%name) .ne. 'emiss_slv') then
               call emission_ar5_ReadSector( 'NO', emis_input_year_nox, idate(2), lsec, field3d, status )
               IF_NOTOK_RETURN(status=1)

               ! convert from (kg NO)/s to (kg N)/s
               field3d = field3d * xmn / xmno
             endif 

          case( 'MACC' )

             ! screen out sectors w/o NOx (bio, oc, nat)
             if ( (trim(sectors_def(lsec)%name) .eq. 'emiss_soil'  ) .or. &
                  (trim(sectors_def(lsec)%name) .eq. 'emiss_anthro') .or. &
                  (trim(sectors_def(lsec)%name) .eq. 'emiss_air'   ) ) then

                if (trim(sectors_def(lsec)%catname) .eq. 'natural') then
                   call emission_macc_ReadSector( emis_input_dir_mac, 'NO', emis_input_year_nat, idate(2), &
                                                  '0.5x0.5_kg.nc', sectors_def(lsec)%name, 'kg NO / s', field3d, status )
                   IF_NOTOK_RETURN(status=1)
                else
                   call emission_macc_ReadSector( emis_input_dir_mac, 'NO', emis_input_year_nox, idate(2), &
                                                  '0.5x0.5_kg.nc', sectors_def(lsec)%name, 'kg NO / s', field3d, status )
                   IF_NOTOK_RETURN(status=1)
                endif
                
                ! convert from (kg NO)/s to (kg N)/s
                field3d = field3d * xmn / xmno
                
             endif

          case( 'ED41' )
             select case(trim(sectors_def(lsec)%name))
             case ('1A3a','1A3b_c_e','1A3d_SHIP','1A3d1')

                ! anthropogenic sources
                call emission_ed4_ReadSector( emis_input_dir_ed4, 'NOx','nox', emis_input_year_nox, idate(2), &
                                              lsec, trim(sectors_def(lsec)%prov), 'kg / s', field3d, status )
                IF_NOTOK_RETURN(status=1;deallocate(field3d))

             end select
             
          case( 'ED42' )

             ! biomass burning (GFED/RETRO/AR5BMB) and transport (ED41) are excluded through ED42_NOX_SECTORS definition
             call emission_ed4_ReadSector( emis_input_dir_ed4, 'NOx', 'nox', emis_input_year_nox, idate(2), &
                                           lsec, trim(sectors_def(lsec)%prov), 'kg / s', field3d, status )
             IF_NOTOK_RETURN(status=1)
                
          case('GFEDv3')
          
             call emission_gfed_ReadSector( emis_input_dir_gfed, 'nox', emis_input_year_nox, idate(2), &
                                             sectors_def(lsec)%name, 'kg NO2 / s', field3d(:,:,1), status )
             IF_NOTOK_RETURN(status=1)

             ! convert from (kg NO2)/s to (kg N)/s
             field3d = field3d * xmn / xmno2
     
          case('RETRO')
             
             call emission_retro_ReadSector( emis_input_dir_retro, 'NOX', emis_input_year_nox, idate(2), &
                                             sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
             IF_NOTOK_RETURN(status=1)

             ! in the file kg(species)/m2/s - what does this mean?? by the numbers I assume kg NO2
             ! convert from (kg NO2)/s to (kg N)/s
             field3d = field3d * xmn / xmno2               
            
          case('MEGAN')      
             !
             !     use soil emissions from MACC
             !
          case('DUMMY')
             
          case default
             write(gol,*) "Error in building list of providers USED_PROVIDERS"; call goErr
             status=1; TRACEBACK; return

          end select

          ! verbose
          if(sum(field3d) < 100.*TINY(1.0) ) then
             if (okdebug) then
                write(gol,'("EMISS-INFO - no NOx emissions found for ",a," ",a," for month ",i2 )') &
                     trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, idate(2) ; call goPr
             endif
             hasData=0
          else
             if (okdebug) then 
                write(gol,'("EMISS-INFO - found NOx emissions for ",a," ",a," for month ",i2 )') &
                     trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, idate(2) ; call goPr
             endif
             ! scale from kg/s to kg/month
             field3d = field3d * sec_month ! kg / month
             hasData=1
          endif

       end if root

       call Par_broadcast(hasData, status)
       IF_NOTOK_RETURN(status=1)

       if (hasData == 0) cycle sec
       
       ! reset temporary arrays
       do region = 1, nregions
          emis3d(region)%d3 = 0.0   
          work3d(region)%d3 = 0.0
          emis2d(region)%surf = 0.0
       end do

       ! distinguish 2d/3d sectors
       if( sectors_def(lsec)%f3d ) then 
          
          ! ---------------------------
          ! 3d data (AIRCRAFT)
          ! ---------------------------
          if (isRoot) then

             ! write some numbers
             call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, 'NOx', xmn, sum(field3d) )

             ! distribute to work arrays in regions
             call Coarsen_Emission( 'NOX '//trim(sectors_def(lsec)%name), nlon360,  nlat180, ar5_dim_3ddata,   &
                                                                          field3d,  work,    add_field, status )
             IF_NOTOK_RETURN(status=1)
 
          end if

          ! scatter, sum up on target array
          do region = 1, nregions

             call scatter(dgrid(region), work3d(region)%d3, work(region)%d3, 0, status)
             IF_NOTOK_RETURN( status=1 )
             
             CALL GET_DISTGRID( dgrid(region), I_STRT=i1, J_STRT=j1)

             ! aircraft data: regrid vertically to model layers
             call emission_ar5_regrid_aircraft( region, i1, j1, work3d(region)%d3, emis3d(region)%d3, status )
             IF_NOTOK_RETURN( status=1 )

             nox_emis_3d(region)%d3 = nox_emis_3d(region)%d3 + emis3d(region)%d3
             
          end do

       else ! ar5_sector is 2d

          ! ---------------------------
          ! 2d data (Anthropogenic, Ships, Biomassburning, ...)
          ! ---------------------------
          
          if (isRoot) then ! print total & regrid

             call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, 'NOx', xmn, sum(field3d(:,:,1)) )

             call coarsen_emission( 'NOx '//sectors_def(lsec)%name, nlon360, nlat180,   field3d(:,:,1), &
                                                                    wrk2D,   add_field, status          )
             IF_NOTOK_RETURN(status=1)

          end if
          
          ! scatter, distribute vertically according to sector, and sum up on target array
          do region = 1, nregions 

             call scatter(dgrid(region), emis2d(region)%surf, work(region)%d3(:,:,1), 0, status)
             IF_NOTOK_RETURN(status=1)
             
             call emission_vdist_by_sector( sectors_def(lsec)%vdisttype, 'NOx', region, emis2d(region), emis3d(region), status )
             IF_NOTOK_RETURN(status=1)

             if ( trim(sectors_def(lsec)%catname) .eq. 'biomassburning') then
                nox_emis_3d_bb(region)%d3 = nox_emis_3d_bb(region)%d3 + emis3d(region)%d3
             else
                nox_emis_3d(region)%d3 = nox_emis_3d(region)%d3 + emis3d(region)%d3
             endif

          end do

       end if ! sectors_def
       
    end do sec ! sectors


    deallocate( field3d )

    do region = 1, nregions
       if (associated(wrk2D(region)%surf)) nullify(wrk2D(region)%surf)
       deallocate( emis3d(region)%d3, emis2d(region)%surf )
       deallocate( work(region)%d3 )
       deallocate( work3d(region)%d3 )
    end do
    
    ! check sectors found
    if( seccount2d /= nox_2dsec ) then 
       write(gol,'(80("-"))') ; call goPr
       write(gol,'("ERROR: 2d sectors do not equal total number:",i4," /= ",i4," !")') seccount2d, nox_2dsec ; call goErr
       write(gol,'(80("-"))') ; call goPr
       status=1; return
    end if
    if( seccount3d /= nox_3dsec ) then 
       write(gol,'(80("-"))') ; call goPr
       write(gol,'("ERROR: 3d sectors do not equal total number:",i4," /= ",i4," !")') seccount3d, nox_3dsec ; call goErr
       write(gol,'(80("-"))') ; call goPr
       status=1; return
    end if

    ! ok
    status = 0

  END SUBROUTINE EMISSION_NOX_DECLARE
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_NOX_BB_DAILY_CYCLE
  !
  ! !DESCRIPTION: Impose daily burning cycle to BMB NOx emissions for current
  !               time step. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_NOX_BB_DAILY_CYCLE( status )
    !
    ! !USES:
    !
    use dims,          only : itaur, nsrce, tref, lm
    use dims,          only : dx, xref, xbeg, yref, ybeg, ndyn_max  
    use partools,      only : myid
    use emission_data, only : emis_bb_trop_cycle, bb_cycle    
    use datetime,      only : tau2date
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)    ::  status
    !
    ! !REVISION HISTORY: 
    !   23 Jan 2014 - Jason Williams - V0
    !                                  
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
    character(len=*), parameter ::  rname = mname//'/emission_nox_bb_daily_cycle'
    
    integer                     :: i,j,l,region, lmr, itim, ntim
    integer                     :: i1, i2, j1, j2, ipos, sec_in_day    
    integer, dimension(6)       :: idater
    real                        :: dtime, dtime2, xlon, xlat
    !
    REG: do region = 1, nregions
       !
       ! Re-initialize the bb NOx array
       !
       nox_emis_3d_bb_app(region)%d3 = 0.0

       call tau2date(itaur(region),idater)
       
       dtime  = float(nsrce)/(2*tref(region))    ! emissions are added in two steps...XYZECCEZYX.
       dtime2 = float(ndyn_max)/(2*tref(region))
       ntim   = 86400/nint(dtime2)               ! number of timesteps in 24 hours.

       sec_in_day = idater(4)*3600 + idater(5)*60 + idater(6) ! elapsed seconds this day
       itim       = sec_in_day/nint(dtime2)+1 ! time interval 
       
       if(okdebug) then
          write(gol,*)'emission_nox_bb_daily_cycle in region ',region,' at date: ',idater, ' with time step:', dtime,' on ',myid 
          call goPr
       end if

       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )          
       lmr = lm(region)
       
       if( emis_bb_trop_cycle) then       
          do l=1,lmr
          do j=j1,j2
          do i=i1,i2
             xlat = ybeg(region) + (j-0.5)*dy/yref(region)
             if (xlat .gt. -20 .and. xlat .lt. 20) then
                ! apply emission cycle in tropics only
                ! itim = 1 and lon = -180 --->position 1
                ! itim = ntim ant lon = -180 --->position ntim, etc.
                ! itim = 1 and lon = 0 ---->position ntim/2
                xlon = xbeg(region) + (i-0.5)*dx/xref(region)
                ipos = 1 + mod(int((xlon+180.)*ntim/360.) + (itim-1),ntim)  !position in array depending on time and lon.
                nox_emis_3d_bb_app(region)%d3(i,j,l) = nox_emis_3d_bb(region)%d3(i,j,l)*bb_cycle(region)%scalef(ipos)
             else
                nox_emis_3d_bb_app(region)%d3(i,j,l) = nox_emis_3d_bb(region)%d3(i,j,l)
             endif
          enddo
          enddo
          enddo
       else
          nox_emis_3d_bb_app(region)%d3 = nox_emis_3d_bb(region)%d3
       endif
           
    end do REG
  
    if(okdebug) then
       write(gol,*) 'end of emission_nox_bb_daily_cycle'; call goPr
    end if
    
    status=0

  END SUBROUTINE EMISSION_NOX_BB_DAILY_CYCLE

  
#ifndef without_convection  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  lightningNOx
  !
  ! !DESCRIPTION: Calculates NOx emissions from lightning as input for
  !     photochemistry module. NOx lightning is calculated using a linear
  !     relationship between lightning flashes and convective precipitation
  !     
  !     * total annual production is approximately 5 Tg(N)/yr
  !     * marine lightning is ten times less active 
  !     * fraction of cloud-to-ground over total flashes is determined by
  !       4th order polynomial fit of the cold cloud thickness (Price and
  !       Rind, GRL 1993).
  !     * NOx production per IC and CG flash is according to Price et al, 
  !       JGR, 1997.
  !     * vertical NOx profile is an approximation of the 'outflow' profile
  !       adopted from Pickering et al., JGR 1998.
  !     
  !     Calculate distribution of lightning using cloudtop heights
  !     of deep convection, cloud cover and convective precipitation
  !     
  !     Reference: E. Meijer, KNMI.
  !                Physics and Chemistry of the earth, Manuscript ST6.03-4
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE LIGHTNINGNOX(region, I1, J1, emilig, status)
    !
    ! !USES:
    !    
    USE dims,          only : im,jm,lm,ybeg,yref
    use Dims,          only : CheckShape
    USE Binas,         only : Avog
    use chem_param,    only : xmn
    use partools,      only : isRoot,     par_reduce
    USE toolbox,       only : ltropo_ifs, lvlpress
    USE meteodata,     only : m_dat,      phlb_dat
    use meteodata,     only : temper_dat, gph_dat, cp_dat
    USE global_data,   only : region_dat, conv_dat
    USE emission_data, only : plandr
    !
    ! !INPUT PARAMETERS:
    !
    integer,     intent(in) :: region, i1, j1
    !
    ! !OUTPUT PARAMETERS:
    !    
    real,       intent(out) :: emilig(i1:,j1:,:)  ! lighting emissions (kg N/s)
    integer,    intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !    ? ??? 2001 - Ernst Meijer - Set up
    !    ? ??? 2002 - Olivie, van Weele - Revisions
    !    ? Jul 2002 - Frank Dentener - adapted for TM5
    !    ? Jan 2003 - Maarten krol - adapted for NEW TM5
    !    1 Oct 2010 - Achim Strunk - protex
    !   27 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !   21 Aug 2013 - Ph. Le Sager - update fscalelig for 1x1
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/lightningNOx'

    real, dimension(:,:,:), pointer    :: phlb  ! pressure (Pa) (1:lm+1)
    real, dimension(:,:,:), pointer    :: m     ! air mass (kg)
    real, dimension(:,:,:), pointer    :: gph   ! height (incl. oro)
    integer,dimension(:,:), pointer    :: lutop       ! cloud_top
    integer,dimension(:,:), pointer    :: lubottom    ! cloud base
    real,dimension(:),      pointer    :: dxyp        ! area (m)
    real,dimension(:,:,:),  pointer    :: t           ! temperature (K)
    real,dimension(:,:,:),  pointer    :: cp          ! convective precipitation (mm/day)
    real,dimension(:,:),    pointer    :: plandregion ! landfraction (0-1)
    !
    real,dimension(:),allocatable      :: gphx,tx
    real                  :: top                ! cloudtop (km)
    real                  :: cldd               ! cold cloud extension (km) (0 deg - tep)
    integer               :: ibase              ! base of cloud layer
    integer               :: itropo             ! tropopuase layer
    integer               :: itop               ! top of cloud layer
    integer               :: itmin15            ! layer with the t=-15 isotherm
    integer               :: itmin24            ! layer with the t=-24 isotherm
    integer               :: it0                ! layer with the t=  0 isotherm
    ! 
    real                  :: cpc         ! convective precipitation => (m s^-1)  
    !
    real                  :: fl                 ! flash rate (1/s)
    real                  :: cg                 ! cloud-to-ground flashe rate (1/s)
    real                  :: ic                 ! intra-cloud flash rate (1/s)
    real                  :: pnocg,pnoic        ! molecules NO produced per CG, IC flash
    ! DIAGNOSTIC variables
    real                  :: flashglob          ! global flash frequency (1/s)
    real                  :: flashtrop          ! tropical flash frequency (1/s)
    real                  :: ic_flashglob       ! global flash frequency (1/s)
    real                  :: ic_flashtrop       ! tropical flash frequency (1/s)
    real                  :: noxltrop           ! tropical NOxL (kg/s)

    !
    logical               :: cld         ! flag for initialisation and cloud
    !
    integer               :: i,j,l,ll,nlay,i2,j2
    real                  :: dsum
    real                  :: airmass
    integer               :: level_pblh
    real                  :: xlat, dlat, tot_emilig
    logical               :: lightning_output = .false.  ! switch for diagnostic output
    logical               :: lightning_output_2 = .false.

    ! for buggy MPI (see comment in budget_global.F90 for details)
    real :: flashglob_all, ic_flashglob_all, flashtrop_all, ic_flashtrop_all, noxltrop_all, tot_emilig_all

    
    
    ! --- begin --------------------------------
    
    CALL GET_DISTGRID( dgrid(region), I_STOP=i2, J_STOP=j2 )
    
    call CheckShape( (/i2-i1+1,j2-j1+1,lm(region)/), shape(emilig), status )
    IF_NOTOK_RETURN(status=1)
    
    m           =>      m_dat (region)%data
    phlb        =>   phlb_dat (region)%data
    cp          =>     cp_dat (region)%data
    plandregion =>     plandr (region)%surf
    dxyp        => region_dat (region)%dxyp
    lubottom    =>   conv_dat (region)%cloud_base
    lutop       =>   conv_dat (region)%cloud_top
    t           => temper_dat (region)%data
    gph         =>    gph_dat (region)%data
    !
    !
    allocate(gphx(0:lm(region)))   ! note now from 0-->lm
    allocate(tx(lm(region)))
    !
    ! FD region coordinates are needed for determining statistics in tropics
    ! (-30 N,30N) and for excluding polar lightning (75N, 75 S) the emissions
    ! in parent regions containing a zoom are only set to zero after budget
    ! calculations using zoomed
    !
    ! Initialising statistics
    !  
    flashglob    = 0.
    flashtrop    = 0.
    ic_flashglob = 0.
    ic_flashtrop = 0.
    noxltrop     = 0.
    !  
    ! initialising lightning emission  
    !
    emilig(:,:,:) = 0. ! (im,jm,lm)
    dlat = dy/yref(region) 

    !
    do j=j1,j2

       xlat = float(ybeg(region)) + j*dlat   ! southern edge of gridbox
       if(xlat < -75.0 .or. (xlat+dlat) > 75.0) cycle   ! exclude poles....

       do i=i1,i2
          !   
          fl      = 0.
          cldd    = 0.
          cg      = 0.
          ic      = 0.
          !
          ibase   = 0
          itop    = 0
          itmin24 = 0
          itmin15 = 0
          it0     = 0
          !
          cld     = .false.
          !
          cpc = cp(i,j,1)    
          !old data    cpc = cp(i,j,1)     / 1000./86400.   ! mm/day => m/s
          !
          if (cpc.gt.0.) then       
             tx(:)=t(i,j,:)
             gphx(0:lm(region))=gph(i,j,1:lm(region)+1)   !note the bounds
             ibase  = lubottom(i,j)
             itop   = lutop(i,j)
             itropo = ltropo_ifs(region,i,j,tx,lm(region))

             if (ibase.ne.0.and.itop.ne.0) then
                do l = itop, ibase, -1
                   if (tx(l).le.249.15) itmin24=l
                   if (tx(l).le.258.15) itmin15=l
                   if (tx(l).le.273.15) it0 = l
                enddo !l
                top = (gphx(itop)+gphx(itop-1))/2000. ! cloud top (km)
                if (itmin24.ne.0.and.top.gt.5) then
                   cld = .true.     ! IF CLOUD REGIONS, IT IS A DEEP CLOUD
                   if (itop.gt.itropo) itop    = itropo
                   if (it0 .gt.itropo) it0     = itropo
                endif ! itmin24
             endif !ibase ne 0
             !
             if (cld) then
                !fd         top = (gphx(itop)+gphx(itop-1))/2000.     ! cloud top (km)
                cldd= top - (gphx(it0)+gphx(it0-1))/2000. ! cold top (km)
                fl = fscalelig *4.e6 * cpc * dxyp(j)*1.e-12
                fl = (0.9*plandregion(i,j) + 0.1) * fl 
                if (cldd.ge.5.5) then
                   cg = fl / (.021*cldd**4-.648*cldd**3+7.493*cldd**2-36.54*cldd+64.09)
                   ic = fl - cg
                else
                   ! changed from [0.;fl] to [fl; fl-cg] (TvN, PLS, 22-04-2013)
                   ! this increases the LiNOx by ~8.2 TgN/yr
                   cg = fl
                   ic = fl-cg
                endif !cldd

                ! Price et al. (JGR, 1997) assumed that cloud-ground flashed
                ! are 10 times more efficient in producing NOx than intraground flashes.
                ! They proposed a production efficiency of 6.7e26 and 6.7e25 molecules NO
                ! per CG and IC flash, respectively.
                ! These values were also adopted by Pickering et al. (JGR, 1998).
                ! Ridley et al. (Atm. Env., 2005) argued that the production efficiency
                ! is similar for both types, and Ott et al. (JGR, 2010)
                ! later come to a similar conclusion.
                !
                ! changed pnocg factor from 6.7e9 to 6.7e8 (TvN, PLS, 22-04-2013),
                ! the value used for pnoic
                !
                ! The new value of 6.7e25 molecules NO per flash corresponds
                ! to about 112 mol NO per flash, which is lower than current estimates.
                ! Ott et al. find a mean of 500 mol per flash, 
                ! while Finney et al. (ACP, 2016) use 250 mol per flash in their model.
                ! Assuming these estimates to be realistic,
                ! this implies that the resolution dependent scale factor fscalelig
                ! can be written as the product of two factors:
                ! - a resolution independent factor that increases the production efficiency
                !   to more realistic value. 
                !   Assuming a production efficiency between 250 and 500 mol per flash,
                !   this factor is in the range 2.2 to 4.5.
                ! - a resolution dependent factor that scales the flash rates.

                ! Currently we use an target total NOx production of 6.0 Tg N/yr,
                ! using meteorology for the reference year 2006.
                ! When ERA-Interim is used (also for the convective fluxes)
                ! this results in a scale factor of about 11 (see above). 
                ! Then the scale factor for the flash rate is in the range 2.5 to 5.
                !
                ! Reducing the target (e.g to 3.0 TgN/yr),
                ! the scale factor for the flash rate would be reduced by proportionally.
                !
                ! The estimated range for the global NOx production from lightning
                ! based on the review by Schumann and Huntrieser (ACP, 2007)
                ! is 5 +- 3 Tg N/yr (or 2 to 8 Tg N/yr).
                !
                pnocg = 1e17*6.7e8*cg *xmn*1e-3/Avog
                pnoic = 1e17*6.7e8*ic *xmn*1e-3/Avog

                ! DISTRIBUTION of LNOx over the COLUMN
                !
                ! assume all IC-LNOx and 70% of CG-LNOx betweem t=-15 and cloudtop;
                !     assume                 10% of CG-LNOx between EARTH SURFACE and t=-15
                !     assume                 20% of CG-LNOx in BOUNDARY LAYER
                !         
                ! To avoid dependency on the vertical resolution:            
                ! - surface emission in lowest layers : boundary layer height;
                ! - LNOx within one of these three regions is distributed proportional to the mass
                !   of each layer.
                ! - the algorithm assumes that itropo is 3 or more. If it is less (typically 1 if
                !   not found and surface pressure is so low that the default at 100hPa is in the
                !   first layer), then we set it to 3, distributing 70% of GC in layer 2 and 10% in
                !   layer 1.
                dsum = 0.7*pnocg+pnoic
                ! determining nlay
                itropo=max(3,itropo) ! assumption of the algorithm
                itop = min(itop,itropo-1)
                nlay = itop - itmin15  
                if (nlay.le.0) then
                   itmin15=itop-1
                   nlay = 1
                   if (lightning_output_2) write(6,*) 'WARNING noxlight_cvp: itmin15>=itropo: ',i,j,itropo,itmin15
                endif !nlay le 0
                ! distributing according to airmass
                airmass = 0.
                do l=itmin15+1,itop
                   airmass = airmass + m(i,j,l)
                enddo
                do l=itmin15+1,itop
                   emilig(i,j,l) = dsum*m(i,j,l)/airmass
                enddo

                ! distributing 10% of CG LNOx between EARTH SURFACE and t=-15
                dsum = 0.1 * pnocg
                ! distributing according to air mass
                airmass = 0.
                do l=1,itmin15
                   airmass = airmass + m(i,j,l)
                enddo
                do l=1,itmin15
                   emilig(i,j,l) = emilig(i,j,l) + dsum*m(i,j,l)/airmass
                enddo

                ! distributing 20% of CG LNOx between ground pressure and 0.8*ground pressure
                dsum = 0.2 * pnocg
                level_pblh = lvlpress(region,0.8*phlb(i,j,1),phlb(i,j,1))
                ! distributing according to airmass
                airmass = 0.
                do l=1,level_pblh
                   airmass = airmass + m(i,j,l)
                enddo
                do l=1,level_pblh
                   emilig(i,j,l) = emilig(i,j,l) + dsum*m(i,j,l)/airmass 
                enddo
                flashrate_lightning(region)%d23(i,j)=fl
                if (lightning_output) then     ! CALCULATE GLOBAL/TROPICAL flash, NOxL rates  
                   flashglob = flashglob + fl
                   ic_flashglob = ic_flashglob + ic
                   select case(nint(xlat))   ! xlat is the southern edge of box j....
                   case(-30:29)
                      flashtrop = flashtrop + fl
                      ic_flashtrop = ic_flashtrop + ic
                      do l = 1, lm(region)
                         noxltrop = noxltrop + emilig(i,j,l)
                      enddo
                   case default 
                   end select
                endif   ! lightning_output

                
                
             endif  !cld = .true.
          endif !cpc.gt.0.
       enddo !i
    enddo !j

    if (lightning_output) then

       call par_reduce( flashglob    , 'sum', flashglob_all    , status )
       call par_reduce( ic_flashglob , 'sum', ic_flashglob_all , status )
       call par_reduce( flashtrop    , 'sum', flashtrop_all    , status )
       call par_reduce( ic_flashtrop , 'sum', ic_flashtrop_all , status )
       call par_reduce( noxltrop ,     'sum', noxltrop_all     , status )
       tot_emilig = sum(emilig)
       call par_reduce( tot_emilig,    'sum', tot_emilig_all,    status )
       
       if (isRoot) then
          write(gol,*) 'EMISS-INFO - global lightning frequency     = ',flashglob_all,' s-1'                        ; call goPr
          write(gol,*) 'EMISS-INFO - ic global lightning frequency  = ',ic_flashglob_all,' s-1'                     ; call goPr
          write(gol,*) 'EMISS-INFO - tropical lightning frequency   = ',flashtrop_all,' s-1'                        ; call goPr
          write(gol,*) 'EMISS-INFO - ic tropical lightning frequency= ',ic_flashtrop_all,' s-1'                     ; call goPr
          write(gol,*) 'EMISS-INFO - global lightning emission      : ',tot_emilig_all*1.e-9*365.*86400.,' Tg[N]/a' ; call goPr
          write(gol,*) 'EMISS-INFO - tropical lightning emission    : ',noxltrop_all*1.e-9*365.*86400.,' Tg[N]/a'   ; call goPr
       end if
    endif

    nullify(m)
    nullify(phlb)
    nullify(gph)
    nullify(cp)
    nullify(plandregion)
    nullify(lubottom)
    nullify(lutop)
    nullify(dxyp)
    nullify(t)
    deallocate(gphx)
    deallocate(tx)

    status=0

  END SUBROUTINE LIGHTNINGNOX
#endif
  
END MODULE EMISSION_NOX