ecearth3/sources/tm5mp/proj/cb05/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 : LCMIP6_CH4 
  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. 'CMIP6') then
                ! Historical emissions of CH4 from aircraft are zero everywhere
                ! and it is assumed that this will also be the case for the future scenarios
                ch4_2dsec = ch4_2dsec + providers_def(lprov)%nsect2d
             elseif (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
    ! ---------------- CMIP6 - AR5 - EDGAR4 - ETC
    use emission_data, only : ch4_fixyear
    use emission_data, only : emis_input_year_ch4, emis_input_year_nat
    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 : read_cmip6_zch4
    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
    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
    integer                     :: target_year, target_month
    integer                     :: lmr, i1, i2, j1, j2

    ! Emissions arrays
    real, dimension(:,:,:), allocatable :: field3d, field3d_p, field3d_n
    type(d3_data), dimension(nregions)  :: emis3d, work, work3D
    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

          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
       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_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

          ! Historical emissions of CH4 from aircraft are zero everywhere
          ! and it is assumed that this will also be the case for the future scenarios
          if ((trim(sectors_def(lsec)%prov).eq.'CMIP6') .and. (trim(sectors_def(lsec)%catname) .eq. 'aircraft')) cycle

          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( 'CMIP6' )

               call emission_cmip6_ReadSector( 'CH4', emis_input_year_ch4, idate(2), lsec, field3d, status )
               IF_NOTOK_RETURN(status=1;deallocate(field3d))

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

             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_ch4, 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_ch4, 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_ch4, 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_ch4, 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_ch4, 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_nat, 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 

             ! ---------------------------------------
             ! 3d data (AIRCRAFT), available for CMIP6
             ! ---------------------------------------

             if (isRoot) then

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

                ! distribute to work arrays in regions
                call Coarsen_Emission( 'CH4 '//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 )

                ch4_emis_3d(region,seccount3d)%d3 = ch4_emis_3d(region,seccount3d)%d3 + emis3d(region)%d3

             end do
 
          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 )
          deallocate( work3d(region)%d3 )
          deallocate( emis3d(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 (LCMIP6_CH4) then
          write(gol,'("Reading CMIP6 monthly zonal mean CH4 mixing ratios")'); call goPr 

          ! previous month
          target_month = idate(2)-1
          target_year = MIN(2500,MAX(emis_input_year_ch4,1850))
          if (idate(2) .eq. 1) then
             target_month = 12
             if (.not.ch4_fixyear) then
                target_year = MIN(2500,MAX(idate(1)-1,1850))
             endif
          endif
          write (gol,*) '   Use year ', target_year,' and month ', target_month,' for previous month'; call goPr
          call read_cmip6_zch4(field3d_p(1,:,1),target_year,target_month, status)

          ! current month 
          target_month = idate(2)
          target_year = MIN(2500,MAX(emis_input_year_ch4,1850))
          write (gol,*) '   Use year ', target_year,' and month ', target_month,' for current month'; call goPr 
          call read_cmip6_zch4(field3d(1,:,1),target_year,target_month, status) 

          ! next month
          target_month = idate(2)+1
          target_year = MIN(2500,MAX(emis_input_year_ch4,1850))
          if (idate(2) .eq. 12) then
             target_month = 1
             if (.not.ch4_fixyear) then
                target_year = MIN(2500,MAX(emis_input_year_ch4+1,1850))
             endif
          endif
          write (gol,*) '   Use year ', target_year,' and month ', target_month,' for next month'; call goPr
          call read_cmip6_zch4(field3d_n(1,:,1),target_year,target_month, status)

          !write (gol,'("STOP AFTER TESTING READING OF CMIP6 CH4")') ; call goErr
          !status=1; TRACEBACK; return

       else if (emis_ch4_single) then ! fixed value
 
          ! TvN: bug fix
          ! unit should be in ppbv as in CMIP6 and NOAA/GMD files
          !field3d(:,:,:)=emis_ch4_fixed_ppb*1.0e-9
          field3d(:,:,:)=emis_ch4_fixed_ppb
          

       else ! read NOAA/GMD 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
    use partools,      only : par_reduce_element
#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                :: zvmr_obs, f_ratio, vmr_min 
    integer             :: i, j, l, lm_ch4, idateline
    real, allocatable   :: vmr_mod(:,:), zvmrsum_local(:)
    real, allocatable   :: zvmr_mod(:)
    
    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(zvmr_mod(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

          ! Historical emissions of CH4 from aircraft are zero everywhere
          ! and it is assumed that this will also be the case for the future scenarios
          if ((trim(sectors_def(lsec)%prov).eq.'CMIP6') .and. (trim(sectors_def(lsec)%catname) .eq. 'aircraft')) cycle

          ! 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_single .and. 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

    ! In simulations without CH4 emissions,
    ! the CH4 mixing ratios are prescribed (not nudged).
    !---------------------------------------------------
#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) 
       ! unit of zch4 is in ppb:
       toadd  = (zch4(region)%d2(j) - ch4ppb) * 1e-9 * m(i,j,l) * 1e3 / xmair  ! moles ch4 
       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

    ! If CH4 emissions are used 
    ! the CH4 mixing ratios are nudged,
    ! to either zonal means from CMIP6 or
    ! zonal background values from NOAA/GMD.
    !---------------------------------------
#ifdef with_ch4_emis

    if (LCMIP6_CH4) then
       ! Nudging for CMIP6 based on zonal means
       ! adapted from code for nudging to HALOE

       allocate(vmr_mod(i1:i2,j1:j2))
       allocate(zvmrsum_local(j1:j2))

       ! convert CH4 mass to volume mixing ratios (mol/mol)
       vmr_mod = (rm(i1:i2,j1:j2,1,ich4)/m(i1:i2,j1:j2,1))*(xmair/xmch4)
       zvmrsum_local=sum(vmr_mod,dim=1) ! zonal total for this subdomain

       CALL PAR_REDUCE_ELEMENT(zvmrsum_local, 'SUM', zvmr_mod, status, all=.true., row=.true.)
       IF_NOTOK_RETURN(status=1)

       ! zonal average
       zvmr_mod(:)=zvmr_mod(:)/im(1)

       deallocate(vmr_mod)
       deallocate(zvmrsum_local)

    else
       ! Nudging to zonal background values from NOAA/GMD
       ! 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

       ! zonal_mod is the zonal minimum volume mixing ratio (mol/mol) in the model
       ! use arbitary large initial value
       zvmr_mod(:)=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
             vmr_min = (rm(idateline,j,1,ich4)/m(idateline,j,1))*(xmair/xmch4)
             zvmr_mod(j) = min(zvmr_mod(j),vmr_min) ! 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(zvmr_mod, status, row=.true.)
       IF_NOTOK_RETURN(status=1)

    endif ! LCMIP6_CH4
   
    ! nudging like in the stratosphere for ozone would be:
    ! rm = (rmold+dtime*gnud*zvmr_obs)/(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 CMIP6 mean/observations
    ! and the zonal mean/minimum surface model field: 
    ! f_ratio = zvmr_obs / zvmr_mod
    ! 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:
    ! vmr = vmrold * exp(-dtime*gnud) + zvmr_obs * (1 - exp(-dtime*gnud)) or:
    ! vmr = vmrold - (1 - exp(-dtime*gnud))(vmrold - vmrobs)  or:
    ! vmr = vmrold + vmrold*(f_ratio-1)*(1 - exp(-dtime*gnud))
    ! toadd = vmrold*(f_ratio-1)*(1 - exp(-dtime*gnud)) 

    dtime=nsrce/(2*tref(region))    

    do j = j1, j2
 
       ! interpolate in time and convert to volume mixing ratio (mol/mol)
       call interp_zch4(region,j,idater,itaur(region),zvmr_obs)

       ! trap bad cases (0. as initial conditions,...)
       if (zvmr_mod(j)==0.) then
          f_ratio=1.
       else
          f_ratio = zvmr_obs/zvmr_mod(j)  ! f_ratio should be always very close to 1 
                                          ! (forcing f_ratio to 1 removes the nudging)
       end if

       ! expected behavior is nudging within 10%
       !PLS: commented since it appears too often and fills log file 
       !if (f_ratio>1.1 .or. f_ratio<0.9) then
           !write(gol,*)"WARNING: CH4 nudging larger than 10% for region/ijl",region,i,j,l; call goPr
       !end if
       
       do l = 1, lm_ch4
          do i = i1, i2
             
             toadd  = rm(i,j,l,ich4) * (f_ratio - 1.) * (1. - exp(- dtime*gnud_ch4)) *1e3 / xmair  ! moles ch4

             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(zvmr_mod)

    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, zch4_obs)
    !
    ! !USES:
    !
    use datetime, only : date2tau
    !
    ! !INPUT PARAMETERS:
    !
    integer,               intent(in) :: region,j
    integer, dimension(6), intent(in) :: idater
    integer(kind=8),       intent(in) :: itau
    !
    ! !OUTPUT PARAMETERS:
    !
    real, intent(out) :: zch4_obs
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer, dimension(6) :: idate_c, idate_p, idate_n
    integer               :: nmon, nyear
    integer(kind=8)       :: 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
       zch4_obs = 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
       zch4_obs = 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

    zch4_obs = zch4_obs *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