ecearth3/sources/tm5mp/proj/cbm4/emission_ch4.F90

!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') 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_CH4
!
! !DESCRIPTION: Perform CH4 emissions needed for TM5 CBM4 version.
!               For AR5 runs, fix the CH4 concentration and keep track
!               of the 3D field where CH4 is added
!               (with_ch4_emissions). 
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_CH4
  !
  ! !USES:
  !
  use GO,            only : gol, goPr, goErr
  use tm5_distgrid,  only : dgrid,     get_distgrid, scatter, scatter_i_band
  use partools,      only : isRoot,    par_broadcast
  use dims,          only : nregions,  idate,   okdebug
  use global_types,  only : emis_data, d3_data, d2_data, d23_data
  use emission_data, only : emis_ch4_single,    emis_ch4_fix3d
  use emission_data, only : emis_ch4_fixed_ppb
  use emission_read, only : used_providers_ch4, has_ch4_emis

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public   :: Emission_CH4_Init
  public   :: Emission_CH4_Done
  public   :: emission_ch4_declare
  public   :: emission_ch4_apply
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'emission_ch4'

#ifdef with_ch4_emis

  type(emis_data), dimension(:,:), allocatable :: ch4_emis_2d
  type(d3_data),   dimension(:,:), allocatable :: ch4_emis_3d

  integer :: ch4_2dsec, ch4_3dsec

#endif

  logical, allocatable :: has_data_3d(:), has_data_2d(:)

  type(d2_data),  target :: zch4(nregions), tmp2d(nregions)
  type(d23_data), target :: wrk_c(nregions)
  type(d2_data),  target :: zch4_p(nregions) ! previous month
  type(d23_data), target :: wrk_p(nregions)
  type(d2_data),  target :: zch4_n(nregions) ! next month
  type(d23_data), target :: wrk_n(nregions)
  !
  ! !REVISION HISTORY: 
  !    1 Oct 2010 - Achim Strunk - overhaul for AR5
  !   28 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
  !                              - made zch4 of type d2_data
  !   20 Nov 2012 - Ph. Le Sager - fix runs longer than a month
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_CH4_INIT
  !
  ! !DESCRIPTION: Allocate memory.
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Emission_CH4_Init( status )
    !
    ! !USES:
    !
    use dims,          only : jm, lm
    use emission_read, only : providers_def, numb_providers, ed42_nsect_ch4
    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
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - adapted for AR5
    !   28 Jun 2012 - Ph. Le Sager - always allocate zch4
    !                              - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    ! --- begin --------------------------------------

#ifdef with_ch4_emis

    ! --- count sectors
    if(has_ch4_emis) then 

       ch4_2dsec = 0
       ch4_3dsec = 0

       do lprov = 1, numb_providers
          if (count(used_providers_ch4.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
                ch4_2dsec = ch4_2dsec + n_ar5_ant_sec*count('CH4'.eq.ar5_cat_ant) + &
                                        n_ar5_shp_sec*count('CH4'.eq.ar5_cat_shp)

                if (LAR5BMB) ch4_2dsec = ch4_2dsec + n_ar5_bmb_sec*count('CH4'.eq.ar5_cat_bmb)

                ch4_3dsec = ch4_3dsec + n_ar5_air_sec*count('CH4'.eq.ar5_cat_air)

             elseif (trim(providers_def(lprov)%name) .eq. 'ED42') then
                
                ch4_2dsec = ch4_2dsec + ed42_nsect_ch4
                ! no 3d sectors in EDGAR 4.2
             else
                ch4_2dsec = ch4_2dsec + providers_def(lprov)%nsect2d
                ch4_3dsec = ch4_3dsec + providers_def(lprov)%nsect3d
             endif
          endif
       enddo

       allocate( ch4_emis_2d( nregions, ch4_2dsec ) )
       allocate( ch4_emis_3d( nregions, ch4_3dsec ) )

       allocate( has_data_2d(ch4_2dsec)) ; has_data_2d=.false.
       allocate( has_data_3d(ch4_3dsec)) ; has_data_3d=.false.

    end if
#endif


    REG: do region=1,nregions

       lmr = lm(region) ; jmr = jm(region)

       CALL GET_DISTGRID( dgrid(region), i_strt=i1, j_strt=j1, i_stop=i2, j_stop=j2)

#ifdef with_ch4_emis
       if(has_ch4_emis) then
          
          ! --- allocate information arrays (2d and 3d)
       
          do lsec=1,ch4_2dsec
             allocate( ch4_emis_2d(region,lsec)%surf(i1:i2, j1:j2) )
          end do

          do lsec=1,ch4_3dsec
             allocate( ch4_emis_3d(region,lsec)%d3(i1:i2, j1:j2, lmr) )
          end do
       end if
#endif
       
       ! 1d-constraining surface array
       allocate( zch4(region)%d2(j1:j2) )
       allocate( zch4_p(region)%d2(j1:j2) )
       allocate( zch4_n(region)%d2(j1:j2) )
       allocate( tmp2d(region)%d2(j2-j1+1))
       zch4(region)%d2  = 0.0
       zch4_p(region)%d2  = 0.0
       zch4_n(region)%d2  = 0.0

       ! work arrays
       if(isRoot) then
          allocate( wrk_c(region)%d23(1,jmr))
          allocate( wrk_p(region)%d23(1,jmr))
          allocate( wrk_n(region)%d23(1,jmr))
       else
          allocate( wrk_c(region)%d23(1,1))
          allocate( wrk_p(region)%d23(1,1))
          allocate( wrk_n(region)%d23(1,1))
       end if

    enddo REG
    
    status = 0

  END SUBROUTINE EMISSION_CH4_INIT
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_CH4_DONE
  !
  ! !DESCRIPTION: Free memory
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Emission_CH4_Done( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - adapted to new structures
    !   28 Jun 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    ! --- begin --------------------------------------
    
    reg: do region = 1, nregions

       deallocate(  zch4(region)%d2  )
       deallocate(  zch4_n(region)%d2  )
       deallocate(  zch4_p(region)%d2  )
       deallocate( tmp2d(region)%d2  )
       deallocate(   wrk_c(region)%d23 )
       deallocate(   wrk_p(region)%d23 )
       deallocate(   wrk_n(region)%d23 )

#ifdef with_ch4_emis
       if (has_ch4_emis) then
          do lsec=1,ch4_2dsec
             deallocate( ch4_emis_2d(region,lsec)%surf )
          end do
          do lsec=1,ch4_3dsec
             deallocate( ch4_emis_3d(region,lsec)%d3 )
          end do
          deallocate( has_data_2d, has_data_3d)
       end if
#endif

    end do reg

#ifdef with_ch4_emis
    if (has_ch4_emis) then 
       deallocate( ch4_emis_2d )
       deallocate( ch4_emis_3d )
    end if
#endif

    status = 0

  end subroutine Emission_CH4_Done
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_CH4_DECLARE
  !
  ! !DESCRIPTION: Opens, reads and evaluates input files (per month). 
  !               Provides emissions on 2d/3d-arrays which are then added
  !               to mixing ratios in routine *apply. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_CH4_DECLARE( status )
    !
    ! !USES:
    !
    use toolbox,       only : coarsen_emission
    use dims,          only : im, jm, lm, idate, sec_month, nlon360, nlat180, iglbsfc
    use chem_param,    only : xmch4,    ich4
    use emission_data, only : msg_emis, LAR5BMB
    ! ---------------- AR5 - EDGAR4 - ETC
    use emission_data, only : emis_input_year
    use emission_data, only : emis_input_dir_retro
    use emission_data, only : emis_input_dir_gfed
    use emission_data, only : emis_input_dir_ed4
#ifdef with_ch4_emis
    use emission_data, only : emis_input_dir_natch4
#endif

    use emission_read, only : emission_ar5_regrid_aircraft

    use emission_read, only : emission_ar5_ReadSector
    use emission_read, only : emission_ed4_ReadSector
    use emission_read, only : emission_gfed_ReadSector
    use emission_read, only : emission_retro_ReadSector
    use emission_read, only : emission_lpj_ReadSector
    use emission_read, only : emission_hymn_ReadSector

    use emission_read, only : sectors_def, numb_sectors
    use emission_read, only : ar5_dim_3ddata
    use emission_read, only : ed42_ch4_sectors

    use Grid,          only : FillGrid
    use meteodata,     only : lli_z
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)     ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk  - adapted for AR5
    !    1 Dec 2011 - Narcisa Banda - added EDGAR 4.1
    !   28 Jun 2012 - Ph. Le Sager  - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter :: rname = mname//'/emission_ch4_declare'
    integer                     :: region, hasData
    integer, parameter          :: add_field=0
    integer, parameter          :: amonth=2
    integer                     :: lsec, nyear, nmon

    ! Emissions arrays
    real, dimension(:,:,:), allocatable :: field3d, field3d_p, field3d_n
    type(d3_data)                       :: emis3d,     work(nregions)
    type(emis_data)                     :: wrk2D(nregions)
    integer                             :: seccount2d, seccount3d

    ! --- begin ----------------------------------

#ifdef with_ch4_emis
    
    write(gol,'(" EMISS-INFO ------------- read CH4 emissions -------------")'); call goPr
    
    if(has_ch4_emis) then 

       do region = 1, nregions
          do lsec=1,ch4_2dsec
             ch4_emis_2d(region,lsec)%surf = 0.0
          end do
          do lsec=1,ch4_3dsec
             ch4_emis_3d(region,lsec)%d3   = 0.0
          end do
       end do

       ! global arrays for coarsening    
       do region = 1, nregions
          if (isRoot)then
             allocate(work(region)%d3(im(region),jm(region),lm(region)))
          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_ch4 == sectors_def(lsec)%prov).eq.0) cycle 

          if ((trim(sectors_def(lsec)%prov).eq.'ED42') .and. (count(ed42_ch4_sectors.eq.sectors_def(lsec)%name) .eq. 0)) cycle
          
          if (associated(sectors_def(lsec)%species)) then  ! AR5 checks
             if (count('CH4'.eq.sectors_def(lsec)%species).eq.0) cycle
             if ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) cycle
          endif

          field3d = 0.0

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

          if (isRoot) then ! READ

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

             case( 'AR5' )

                ! Screen out solvent sector for CH4, 
                ! 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( 'CH4', emis_input_year, idate(2), lsec, field3d, status )
                  IF_NOTOK_RETURN(status=1)
                endif

             case( 'ED41' )
                
                ! only transport sector (others provided by ED42)
                select case(trim(sectors_def(lsec)%name))
                case ('1A3b_c_e','1A3d_SHIP','1A3d1')

                   ! anthropogenic sources (only 2d)
                   call emission_ed4_ReadSector( emis_input_dir_ed4, 'CH4', 'ch4', emis_input_year, 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_CH4_SECTORS definition
                call emission_ed4_ReadSector( emis_input_dir_ed4, 'CH4', 'ch4', emis_input_year, 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, 'ch4', emis_input_year, idate(2), &
                                               sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
                IF_NOTOK_RETURN(status=1;deallocate(field3d))

             case('RETRO')

                call emission_retro_ReadSector( emis_input_dir_retro, 'CH4', emis_input_year, idate(2), &
                                                sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
                IF_NOTOK_RETURN(status=1;deallocate(field3d))

             case( 'LPJ' )

                ! this here is for natural sources (only 2d)
                call emission_lpj_ReadSector( trim(emis_input_dir_natch4)//'/LPJdata-jan2011', emis_input_year, idate(2), &
                                              sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
                IF_NOTOK_RETURN(status=1;deallocate(field3d))

             case( 'HYMN' )

                ! this here is for natural sources (only 2d)
                call emission_hymn_ReadSector( trim(emis_input_dir_natch4), &
                                               sectors_def(lsec)%name, 'kg / s', field3d(:,:,1), status )
                IF_NOTOK_RETURN(status=1;deallocate(field3d))

             case default
                write(gol,*) "Error in buidling list of providers USED_PROVIDERS_CH4"; call goErr
                status=1; TRACEBACK; return

             end select

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

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

          if (hasData == 0) then
             cycle sec
          else
             if ( sectors_def(lsec)%f3d ) then
                has_data_3d(seccount3d)=.true.
             else
                has_data_2d(seccount2d)=.true.
             end if
          end if
       
          ! distinguish 2d/3d sectors
          if( sectors_def(lsec)%f3d ) then 

             write(gol,'("EMISS-ERROR - Unexpected 3D data ")'); call goErr
             status=1; TRACEBACK; return

          else          
             ! ---------------------------
             ! 2d data (Anthropogenic, Biomassburning, Natural)
             ! ---------------------------

             if (isRoot) then ! regrid

                ! write some numbers
                call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, 'CH4', xmch4, &
                               sum(field3d(:,:,1)) )

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

             end if

             do region = 1, nregions
                call scatter(dgrid(region), ch4_emis_2d(region,seccount2d)%surf, work(region)%d3(:,:,1), 0, status)
                IF_NOTOK_RETURN(status=1)
             end do


          end if

       end do sec ! sectors


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

    !
    ! Read in either the zonal mean surface field for this year/month, or set to fixed value
    !
    
    if (isRoot) then

       allocate( field3d(1,nlat180,1) )
       allocate( field3d_p(1,nlat180,1) )
       allocate( field3d_n(1,nlat180,1) )
       field3d = 0.0
       field3d_p = 0.0
       field3d_n = 0.0

       if (emis_ch4_single) then ! fixed value

          field3d(:,:,:)=emis_ch4_fixed_ppb*1.0e-9

       else ! read ERA Interim monthly surface latitudinal distribution

          call read_bkglat_ch4(field3d(1,:,1),idate(1),idate(2), status)

          call prev_mon(idate(1),idate(2),nyear,nmon)
          call read_bkglat_ch4(field3d_p(1,:,1),nyear,nmon, status)

          call next_mon(idate(1),idate(2),nyear,nmon)
          call read_bkglat_ch4(field3d_n(1,:,1),nyear,nmon, status)

       endif

       ! coarsen or distribute to zoom regions:
       do region = 1, nregions

          call FillGrid( lli_z(region),  'n', wrk_c(region)%d23(:,:), &
               lli_z(iglbsfc), 'n', field3d(:,:,1), 'area-aver', status )
          IF_NOTOK_RETURN(status=1)

          call FillGrid( lli_z(region),  'n', wrk_p(region)%d23(:,:), &
               lli_z(iglbsfc), 'n', field3d_p(:,:,1), 'area-aver', status )
          IF_NOTOK_RETURN(status=1)

          call FillGrid( lli_z(region),  'n', wrk_n(region)%d23(:,:), &
               lli_z(iglbsfc), 'n', field3d_n(:,:,1), 'area-aver', status )
          IF_NOTOK_RETURN(status=1)

       end do

       deallocate( field3d, field3d_p, field3d_n )
    end if

    ! scatter along latitude direction, then broadcast to all longitudes
    DO region = 1, nregions

       call SCATTER_I_BAND( dgrid(region), zch4(region)%d2, wrk_c(region)%d23(1,:), status )
       IF_NOTOK_RETURN(status=1)
       call SCATTER_I_BAND( dgrid(region), zch4_p(region)%d2, wrk_p(region)%d23(1,:), status )
       IF_NOTOK_RETURN(status=1)
       call SCATTER_I_BAND( dgrid(region), zch4_n(region)%d2, wrk_n(region)%d23(1,:), status )
       IF_NOTOK_RETURN(status=1)

       tmp2d(region)%d2 = zch4(region)%d2
       call PAR_BROADCAST( tmp2d(region)%d2, status, row=.true. )
       IF_NOTOK_RETURN(status=1)
       zch4(region)%d2 = tmp2d(region)%d2

       tmp2d(region)%d2 = zch4_p(region)%d2
       call PAR_BROADCAST( tmp2d(region)%d2, status, row=.true. )
       IF_NOTOK_RETURN(status=1)
       zch4_p(region)%d2 = tmp2d(region)%d2

       tmp2d(region)%d2 = zch4_n(region)%d2
       call PAR_BROADCAST( tmp2d(region)%d2, status, row=.true. )
       IF_NOTOK_RETURN(status=1)
       zch4_n(region)%d2 = tmp2d(region)%d2

    END DO
    
    ! ok
    status = 0

  END SUBROUTINE EMISSION_CH4_DECLARE
  !EOC
  
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_CH4_APPLY
  !
  ! !DESCRIPTION: Takes monthly emissions, and
  !               - split them vertically
  !               - apply time splitting factors 
  !               - add them up (add_3d) 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_CH4_APPLY( region, status)
    !
    ! !USES:
    !
    use binas,         only : xmair
    use dims,          only : itaur, nsrce, tref
    use dims,          only : im, jm, lm, at, bt
    use datetime,      only : tau2date
    use emission_data, only : emission_vdist_by_sector,   LAR5BMB
    use emission_data, only : do_add_3d, do_add_3d_cycle, bb_cycle
    use emission_data, only : emis_bb_trop_cycle
    use emission_read, only : sectors_def, numb_sectors
    use emission_read, only : ed42_ch4_sectors
    use chem_param,    only : ich4, xmch4
    ! NB use chem_param,    only : ch4_ps

    use meteodata,     only : m_dat
    use global_data,   only : region_dat, mass_dat
    use toolbox,       only : lvlpress
#ifdef with_budgets
    use budget_global, only : budg_dat, nzon_vg 
    use emission_data, only : sum_emission, budemi_dat
#endif
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)         :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        :: status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - adapted to AR5 emissions structures
    !                              - added nudging if emissions
    !   28 Jun 2012 - Ph. Le Sager - nudging in case of with_ch4_emis no more
    !                                limited to 3x2 runs
    !                              - optimized loop order for case of NO emissions
    !                              - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !---------------------------------------------------------------------------
    !BOC

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

    integer, dimension(6)       :: idater
    real                        :: dtime, fraction
    integer                     :: imr, jmr, lmr, lsec
    integer                     :: seccount2d, seccount3d
    type(d3_data)               :: emis3d
    type(emis_data)             :: emis2d

    integer :: i1, j1, i2, j2

    real, dimension(:,:,:), allocatable :: field3d
    real, dimension(:,:,:,:),   pointer :: rm
    real, dimension(:,:,:),     pointer :: m

    real                :: ch4ppb, toadd
    integer             :: nzone_v
    real                :: rmobs,  fmin,  mrch4
    integer             :: i, j, l, lm_ch4, idateline
    real,   allocatable :: rmmod(:)
    
    real, parameter :: gnud_ch4  = 4.0e-6 ! ~3day time scale, see below


    ! --- begin ---------------------------

    call get_distgrid( dgrid(region), i_strt=i1, j_strt=j1, i_stop=i2, j_stop=j2)

    nullify(rm)
    nullify(m)

    rm => mass_dat(region)%rm
    m  => m_dat(region)%data

    allocate(rmmod(j1:j2))

    !---------------------------------------------------------------
    ! CH4 emissions 
    !--------------------------------------------------------------    
#ifdef with_ch4_emis

    if (has_ch4_emis) then
       
       call tau2date(itaur(region),idater)
       dtime=float(nsrce)/(2*tref(region))    !emissions are added in two steps...XYZECCEZYX.

       if(okdebug) then
          write(gol,*)'emission_ch4_apply in region ',region,' at date: ',idater, ' with time step:', dtime
          call goPr
       end if

       ! get a working structure
       lmr = lm(region)
       allocate( emis3d%d3  (i1:i2, j1:j2, lmr) ) ; emis3d%d3 = 0.0
       allocate( emis2d%surf(i1:i2, j1:j2     ) ) ; emis2d%surf = 0.0

       ! count 2d and 3d sectors
       seccount2d = 0
       seccount3d = 0

       ! cycle over sectors
       do lsec = 1, numb_sectors

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

          if ((trim(sectors_def(lsec)%prov).eq.'ED42') .and. (count(ed42_ch4_sectors.eq.sectors_def(lsec)%name) .eq. 0)) cycle 
          
          if (associated(sectors_def(lsec)%species)) then  ! AR5 checks
             if (count('CH4'.eq.sectors_def(lsec)%species).eq.0) cycle
             if ((trim(sectors_def(lsec)%catname) .eq. 'biomassburning').and.(.not.LAR5BMB)) cycle
          endif

          ! default: no additional splitting
          fraction = 1.0

          ! ----------------------------------------------------------------------------------------
          ! distinguish here between sectors and whether they should have additional splitting 
          ! if( sectors_def(lsec)%catname == 'biomassburning' ) fraction = fraction * bb_frac etc...
          ! ----------------------------------------------------------------------------------------
          
          ! distinguish between 2d/3d sectors
          if( sectors_def(lsec)%f3d ) then 

             seccount3d = seccount3d + 1
             if (.not.has_data_3d(seccount3d)) cycle
             emis3d%d3 = ch4_emis_3d(region,seccount3d)%d3

          else
             
             seccount2d = seccount2d + 1
             if (.not.has_data_2d(seccount2d)) cycle
             emis2d%surf = ch4_emis_2d(region,seccount2d)%surf

             ! account for soil sink scale with respect to surface [CH4] over i,j
             if (trim(sectors_def(lsec)%name).eq.'soilconsumption') then

                do j=j1,j2
                   do i=i1,i2

                      if(emis2d%surf(i,j)>0.) &
                           emis2d%surf(i,j) = ( - emis2d%surf(i,j)) * (rm(i,j,1,ich4)/m(i,j,1)/1.0e-6) * (xmair/xmch4)

                   enddo
                enddo

             end if

             emis3d%d3 = 0.0

             ! vertically distribute according to sector
             call emission_vdist_by_sector( sectors_def(lsec)%vdisttype, 'CH4', region, emis2d, emis3d, status )
             IF_NOTOK_RETURN(status=1;deallocate(emis3d%d3))

          end if

          ! add dataset according to sector and category
          if( emis_bb_trop_cycle .and. trim(sectors_def(lsec)%catname) == "biomassburning" ) then 
             call do_add_3d_cycle( region, ich4, i1, j1, emis3d%d3, bb_cycle(region)%scalef, xmch4, xmch4, status,fraction )
             IF_NOTOK_RETURN(status=1)
          else
             call do_add_3d( region, ich4, i1, j1, emis3d%d3, xmch4, xmch4, status, fraction )
             IF_NOTOK_RETURN(status=1)
          end if

       end do

       deallocate( emis3d%d3 )
       deallocate( emis2d%surf )

    end if
    
#endif /* with_ch4_emis */


    !-----------------------------------------------------
    ! Fix methane mixing ratio
    !-----------------------------------------------------
    
    ! Set highest level for nudging
    !------------------------------
#ifdef with_ch4_emis
    !apply nuddging up to 550 hpa
    lm_ch4 = lvlpress(region,55000.,98400.)
#else
    ! default: nudging in the lowest layer
    lm_ch4 = 1
    if( emis_ch4_fix3d ) then
       lm_ch4 = lm(region)
    else 
       ! TvN: in case of a shallow first layer (incl. L60, L91), nudge the lowest 2 layers
       if (at(2)+bt(2)*101325. > 101000.) lm_ch4 = 2
    endif
#endif

    
    ! nudging if NO emissions ( optimized loop order)
    !------------------------------
#ifndef with_ch4_emis

    do l = 1,  lm_ch4
    do j = j1, j2
    do i = i1, i2

       ! budget for fixing:
       ch4ppb = 1e9 * (xmair/xmch4) * rm(i,j,l,ich4) / m(i,j,l) 
       !mvw for ERA Interim we do not apply fac_ch4 because we have zonal monthly mean input files per year
       !mvw unit of zch4 in netcdf files is in ppb
       !             toadd  = (zch4(j)*fac_ch4 *1e9 - ch4ppb) * 1e-9 * m(i,j,l) * 1e3 / xmair  ! moles ch4
       toadd  = (zch4(region)%d2(j) - ch4ppb) * 1e-9 * m(i,j,l) * 1e3 / xmair  ! moles ch4 
       !PLS                ch4_ps(region)%d3(i,j,l) = ch4_ps(region)%d3(i,j,l) + toadd
       ! fix value:  (mol air) * (mol ch4)/(mol air) (kg ch4)/(mol ch4)
       !             rm(i,j,l,ich4) = m(i,j,l)/xmair * ch4_fixed_ppb*1e-9 * xmch4
       !mvw for ERA Interim we do not apply fac_ch4 because we have zonal monthly mean input files per year
       !mvw unit of zch4 in netcdf files is in ppb
       !             rm(i,j,l,ich4) = m(i,j,l)/xmair * zch4(j)*fac_ch4 *1e-9 * xmch4
       rm(i,j,l,ich4) = m(i,j,l)/xmair * zch4(region)%d2(j) *1e-9 * xmch4

#ifdef with_budgets
       nzone_v=nzon_vg(l) 
       budemi_dat(region)%emi(i,j,nzone_v,ich4) =  budemi_dat(region)%emi(i,j,nzone_v,ich4) + toadd
       if(ich4 ==1) sum_emission(region) = sum_emission(region) + toadd*xmch4/1e3  !in kg 
#endif
             
    end do
    end do
    end do

#endif


    ! nudging if emissions
    !------------------------------
#ifdef with_ch4_emis


    ! use arbitary large value for initial rmmod
    rmmod=3.0e-6
    !
    ! Use value at the dateline away from major CH4 sources
    !
    idateline = 1

    ! if current processor include dateline
    if (i1==1) then
       do j=j1,j2
          mrch4 = (rm(idateline,j,1,ich4)/m(idateline,j,1))*(xmair/xmch4)
          rmmod(j) = min(rmmod(j),mrch4) ! mimimum CH4 at latitude J
       end do
    endif

    ! broadcast (from proc with i1==1, ie from root of row_communicator, ie 0, ie default broadcaster)
    call par_broadcast(rmmod, status, row=.true.)
    IF_NOTOK_RETURN(status=1)

        
    do l =  1, lm_ch4
       do j = j1, j2
          do i = i1, i2
             ! nudging full 3D field using surface observations in background atmosphere
             ! the observations represent the minimum of rm in the zonal band
             ! the unit of zch4 in the netcdf files is in ppb
             ! nudging like in the stratosphere for ozone would be rm = (rmold+dtime*gnud*rmobs)/(1.+dtime*gnud)
             ! nudging used here is done without the approximation: 1 - exp(-dtime*gnud) = dtime*gnud
             ! we determine per zonal band the ratio between the minimum surface model field and the observation
             ! rmobs = rmmod * fmin 
             ! we apply the surface ratio to the total 3d field with a slow nudging time scale to distribute the
             ! corrections at this latitude band over all heights 
             ! we use gnud_ch4 = 4.e-7 (1/sec; for ~1 month nudging time scale )
             ! The nudging equation is:
             ! rm = rmold * exp(-dtime*gnud) + rmobs * (1 - exp(-dtime*gnud)) or:
             ! rm = rmold - (1 - exp(-dtime*gnud))(rmold - rmobs)  or:
             ! rm = rmold + rmold*(fmin-1)*(1 - exp(-dtime*gnud))  with: fmin = rmobs/rmmod
             ! toadd = rmold*(fmin-1)*(1 - exp(-dtime*gnud))  with: fmin = rmobs/rmmod
             dtime=nsrce/(2*tref(region))
             call interp_zch4(region,j,idater,itaur(region),rmobs) 

             ! trap bad cases (0. as initial conditions,...)
             if(rmmod(j)==0.) then
                fmin=1.
             else
                fmin = rmobs/rmmod(j)  ! fmin should be always very close to 1 (forcing fmin to 1 removes the nudging)
             end if
             
             ! expected behavior is nudging within 10%
             if(fmin>1.1 .or. fmin<0.9) then
                !PLS: commented since it appears too often and fills log file
                !write(gol,*)"WARNING: CH4 nudging larger than 10% for region/ijl",region,i,j,l; call goPr
             end if

             toadd  = rm(i,j,l,ich4) * (fmin - 1.) * (1. - exp(- dtime*gnud_ch4)) *1e3 / xmair  ! moles ch4

             
             !NB ch4_ps(region)%d3(i,j,l) = ch4_ps(region)%d3(i,j,l) + toadd
             rm(i,j,l,ich4) = rm(i,j,l,ich4) + toadd * xmch4/1e3
#ifdef with_budgets
             nzone_v=nzon_vg(l) 
             ! add the nudging adjustments to the emission budget (now full 3D field filled)
             ! the emissions are not added here
             budemi_dat(region)%emi(i,j,nzone_v,ich4) =  budemi_dat(region)%emi(i,j,nzone_v,ich4) + toadd
#endif            
          end do
       end do
    end do
#endif
    
    !-----------------------------------------------------
    ! Done
    !-----------------------------------------------------
    nullify(rm)
    nullify(m)

    deallocate(rmmod)

    status = 0

  END SUBROUTINE EMISSION_CH4_APPLY
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  READ_BKGLAT_CH4
  !
  ! !DESCRIPTION: Read ERA Interim background surface zonal CH4 for specified
  !                year/month 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE READ_BKGLAT_CH4(field, year, month, status)
    !
    ! !USES:
    !
    USE MDF, ONLY : MDF_OPEN, MDF_NETCDF, MDF_READ, MDF_Get_Var, MDF_Close, MDF_Inq_VarID
    use emission_data, only : emis_zch4_fname
    use dims,          only : nlat180
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,                  intent(out) :: status
    real, dimension(nlat180), intent(out) ::  field
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: year, month
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter :: rname = mname//'/Read_BkgLat_CH4'
    character(len=256)          :: filemon, emis_zch4_fullname
    integer                     :: fid, varid, year_valid
    integer, dimension(2), parameter :: era_interim_avail = (/1989, 2014/)

    ! valid year
    year_valid = MIN(era_interim_avail(2),MAX(year, era_interim_avail(1)))
    
    ! target file name with year e.g. surface_ch4_zm_1989.nc
    write (emis_zch4_fullname,'(a,"/surface_ch4_zm_",i4.4,".nc")') trim(emis_zch4_fname), year_valid

    CALL MDF_Open( TRIM(emis_zch4_fullname), MDF_NETCDF, MDF_READ, fid, status )
    IF_NOTOK_RETURN(status=1)
    
    ! search for the record for requested month:
    if (month < 10) then
       write(filemon,'(a,i1)')'CH4_zm', month
    else
       write(filemon,'(a,i2)')'CH4_zm', month
    endif

    CALL MDF_Inq_VarID( fid, TRIM(filemon), varid, status )
    IF_ERROR_RETURN(status=1)

    if ( varid < 0 ) then
       write (gol,'("WARNING - no surface CH4 concentrations `",a,"`")') trim(emis_zch4_fullname) ; call goPr
       status=1; TRACEBACK; return
    else         
       CALL MDF_Get_Var( fid, varid, field, status )
       IF_NOTOK_RETURN(status=1)
    endif
    
    CALL MDF_Close( fid, status )
    IF_NOTOK_RETURN(status=1)
    
    status = 0

  END SUBROUTINE READ_BKGLAT_CH4
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  INTERP_ZCH4
  !
  ! !DESCRIPTION:
  !    Interpolates in time the latitudinal monthly background CH4
  !    concentration. The data are considered representative of the
  !    16th of each month, except February when it is 15th.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE INTERP_ZCH4(region,j,idater,itau, rmobs)
    !
    ! !USES:
    !
    use datetime, only : date2tau
    !
    ! !INPUT PARAMETERS:
    !
    integer,               intent(in) :: region,j
    integer, dimension(6), intent(in) :: idater
    integer,               intent(in) :: itau
    !
    ! !OUTPUT PARAMETERS:
    !
    real, intent(out) :: rmobs
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer, dimension(6) :: idate_c, idate_p, idate_n
    integer               :: nmon, nyear
    integer               :: itau_c, itau_p, itau_n   

    idate_c(1) = idater(1)
    idate_c(2) = idater(2)
    if (idater(2) == 2) then 
       idate_c(3) = 15
    else
       idate_c(3) = 16
    endif
    idate_c(4:6) = 0
    call date2tau(idate_c,itau_c)

    call prev_mon(idater(1),idater(2),nyear,nmon)
    idate_p(1) = nyear
    idate_p(2) = nmon
    if (idate_p(2) == 2) then 
       idate_p(3) = 15
    else
       idate_p(3) = 16
    endif
    idate_p(4:6) = 0
    call date2tau(idate_p,itau_p)

    call next_mon(idater(1),idater(2),nyear,nmon)
    idate_n(1) = nyear
    idate_n(2) = nmon
    if (idate_n(2) == 2) then 
       idate_n(3) = 15
    else
       idate_n(3) = 16
    endif
    idate_n(4:6) = 0
    call date2tau(idate_n,itau_n)

    if (itau.lt.itau_c) then
       rmobs = zch4_p(region)%d2(j) * float(itau_c-itau)/float(itau_c-itau_p) + zch4(region)%d2(j) * float(itau-itau_p)/float(itau_c-itau_p)
    else
       rmobs = zch4(region)%d2(j) * float(itau_n-itau)/float(itau_n-itau_c) + zch4_n(region)%d2(j) * float(itau-itau_c)/float(itau_n-itau_c)
    endif

    rmobs = rmobs *1.e-9

  END SUBROUTINE INTERP_ZCH4
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  NEXT_MON
  !
  ! !DESCRIPTION: Return month number and year of next month
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE NEXT_MON(year,mon,nyear,nmon)
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)    :: year, mon
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)   :: nyear, nmon
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    if (mon.lt.12) then
       nyear = year
       nmon = mon+1  
    else
       nyear = year+1
       nmon = 1
    endif

  END SUBROUTINE NEXT_MON
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  PREV_MON
  !
  ! !DESCRIPTION: Return month number and year of previous month
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PREV_MON(year,mon,nyear,nmon)
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)    :: year, mon
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)   :: nyear, nmon
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    if (mon.gt.1) then
       nyear = year
       nmon = mon-1  
    else
       nyear = year-1
       nmon = 12
    endif

  END SUBROUTINE PREV_MON
  !EOC

END MODULE EMISSION_CH4