!
#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_CO2
!
! !DESCRIPTION: Hold data and methods for CO2 emissions.
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION_CO2
  !
  ! !USES:
  !
  use GO,            only : gol, goErr, goPr
  use dims,          only : nregions, idate, okdebug
  use global_types,  only : emis_data, d3_data
  use emission_read, only : used_providers, has_emis

  use tm5_distgrid,  only : dgrid, get_distgrid, scatter
  use partools,      only : isRoot, par_broadcast

  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public   :: Emission_CO2_Init     ! allocate dataset
  public   :: Emission_CO2_Done     ! deallocate dataset
  public   :: Emission_CO2_Declare  ! read monthly input
  public   :: Emission_CO2_Apply    ! distribute & add emissions to tracer array
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'emission_co2'

  type( emis_data ), dimension(:,:), allocatable :: co2_emis_2d
  type( d3_data   ), dimension(:,:), allocatable :: co2_emis_3d

  logical, allocatable :: has_data_3d(:), has_data_2d(:)
  integer              :: co2_2dsec, co2_3dsec
  !
  ! !REVISION HISTORY: 
  !    1 Oct 2010 - Achim Strunk  - revamped for AR5
  !   26 Jun 2012 - P. Le Sager   - adapted for lon-lat MPI domain decomposition
  !   14 May 2014 - T. van Noije  - created CO2 version from CO version
  !                                 to be done: convert CMIP5 CO2 input files 
  !                                 for emissions from fossil fuel use and land use
  !                                 into the same format as for the other trace species
  !   23 March 2018 - T. van Noije - included CMIP6 anthropogenic emissions from CEDS
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_CO2_INIT
  !
  ! !DESCRIPTION: Allocate memory to handle the emissions.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_CO2_INIT( status )
    !
    ! !USES:
    !
    use dims,           only : lm
    use emission_read,  only : providers_def, numb_providers
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - v0
    !   26 Jun 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    status = 0
    if(.not. has_emis) return

    ! nb of sectors
    co2_2dsec = 0
    co2_3dsec = 0
    do lprov = 1, numb_providers
         if (count(used_providers.eq.providers_def(lprov)%name)/=0) then 
              co2_2dsec = co2_2dsec + providers_def(lprov)%nsect2d
              co2_3dsec = co2_3dsec + providers_def(lprov)%nsect3d
         endif
    enddo

    allocate( co2_emis_2d( nregions, co2_2dsec ) )
    allocate( co2_emis_3d( nregions, co2_3dsec ) )

    allocate( has_data_2d(co2_2dsec)) ; has_data_2d=.false.
    allocate( has_data_3d(co2_3dsec)) ; has_data_3d=.false.

    ! 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)

       do lsec=1,co2_2dsec
          allocate( co2_emis_2d(region,lsec)%surf(i1:i2,j1:j2) )
       end do
       
       do lsec=1,co2_3dsec
          allocate( co2_emis_3d(region,lsec)%d3(i1:i2,j1:j2,lmr) )
       end do
       
    enddo

    ! ok
    status = 0

  END SUBROUTINE EMISSION_CO2_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_CO2_DONE
  !
  ! !DESCRIPTION: Free memory after handling of the emissions.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_CO2_DONE( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - v0
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    status = 0
    if(.not. has_emis) return

    do region = 1, nregions
       do lsec=1,co2_2dsec
          deallocate( co2_emis_2d(region,lsec)%surf )
       end do
       do lsec=1,co2_3dsec
          deallocate( co2_emis_3d(region,lsec)%d3 )
       end do
    end do

    deallocate( co2_emis_2d )
    deallocate( co2_emis_3d )
    deallocate( has_data_2d, has_data_3d)

    ! ok
    status = 0

  END SUBROUTINE EMISSION_CO2_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_CO2_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_CO2_DECLARE( status )
    !
    ! !USES:
    !
    use toolbox,       only : coarsen_emission
    use dims,          only : im, jm, lm, idate, sec_month, nlon360, nlat180, iglbsfc
    use chem_param,    only : xmco2
    use emission_data, only : msg_emis
    ! ---------------- AR5 - ETC. --------------------
    use emission_data, only : emis_input_year_co2

    use emission_read, only : emission_ar5_regrid_aircraft

    use emission_read, only : emission_cmip6_ReadSector
    use emission_read, only : emission_ar5_ReadCO2

    use emission_read, only : sectors_def, numb_sectors
    use emission_read, only : ar5_dim_3ddata
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk  - adapted for AR5
    !   26 Jun 2012 - P. Le Sager   - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter :: rname = mname//'/emission_co2_declare'
    integer                     :: region, hasData
    integer, parameter          :: add_field=0
    integer, parameter          :: amonth=2
    integer                     :: imr, jmr, lmr, lsec, i1, i2, j1, j2
    ! AR5
    real,dimension(:,:,:),  allocatable  :: field3d,    field3d2
    type(d3_data), dimension(nregions)   :: emis3d, work, work3d
    type(emis_data)                      :: wrk2D(nregions)
    integer                              :: seccount2d, seccount3d
    real                        :: emiss_total=0.

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

    status = 0
    if(.not. has_emis) return

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

    do region = 1, nregions
       do lsec=1,co2_2dsec
          co2_emis_2d(region,lsec)%surf = 0.0
       end do
       do lsec=1,co2_3dsec
          co2_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.eq.sectors_def(lsec)%prov).eq.0) 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( 'CO2', emis_input_year_co2, idate(2), lsec, field3d, status )
            IF_NOTOK_RETURN(status=1;deallocate(field3d))

          case( 'AR5' )

             call emission_ar5_ReadCO2( 'CO2', emis_input_year_co2, idate(2), lsec, field3d, status )
             IF_NOTOK_RETURN(status=1;deallocate(field3d))

          case('DUMMY')

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

          end select

          ! nothing found??? 
          if( abs(sum(field3d)) < 100.*TINY(1.0) ) then
             if (okdebug) then
                write(gol,'("EMISS-INFO - no CO2 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 CO2 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, 'CO2', xmco2, sum(field3d) )
             emiss_total =  emiss_total + sum(field3d)

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

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

          end do

        else ! ar5_sector is 2d

           ! ---------------------------
           ! 2d data
           ! ---------------------------

           if (isRoot) then ! print total & regrid

              call msg_emis( amonth, trim(sectors_def(lsec)%prov), sectors_def(lsec)%name, 'CO2', xmco2, sum(field3d(:,:,1)) )
              emiss_total =  emiss_total + sum(field3d(:,:,1))

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

           end if

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

        endif

    end do sec ! sectors

    ! write some numbers
    call msg_emis( amonth, 'TOTAL', '', 'CO2', xmco2, emiss_total )

    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 /= co2_2dsec ) then
       write(gol,'(80("-"))') ; call goPr
       write(gol,'("ERROR: 2d sectors do not equal total number:",i4," /= ",i4," !")') seccount2d, co2_2dsec ; call goErr
       write(gol,'(80("-"))') ; call goPr
       status=1; return
    end if
    if( seccount3d /= co2_3dsec ) then
       write(gol,'(80("-"))') ; call goPr
       write(gol,'("ERROR: 3d sectors do not equal total number:",i4," /= ",i4," !")') seccount3d, co2_3dsec ; call goErr
       write(gol,'(80("-"))') ; call goPr
       status=1; return
    end if

    ! ok
    status = 0

  END SUBROUTINE EMISSION_CO2_DECLARE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_CO2_APPLY
  !
  ! !DESCRIPTION: Take monthly emissions, and
  !               - split them vertically
  !               - apply time splitting factors 
  !               - add them to tracers (add_3d) 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_CO2_APPLY( region, status )
    !
    ! !USES:
    !
    use dims,          only : itaur, nsrce, tref
    use dims,          only : im, jm, lm
    use datetime,      only : tau2date
    use emission_data, only : emission_vdist_by_sector
    use emission_data, only : do_add_3d
    use chem_param,    only : ico2, xmco2
    use emission_read, only : sectors_def, numb_sectors
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)         ::  region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - AR5
    !   26 Jun 2012 - P. Le Sager  - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter :: rname = mname//'/emission_co2_apply'
    integer, dimension(6)       :: idater
    real                        :: dtime, fraction
    integer                     :: imr, jmr, lmr, lsec, i1, i2, j1, j2
    integer                     :: seccount2d, seccount3d
    type(d3_data)               :: emis3d

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

    if( okdebug ) then
       write(gol,*) 'start of emission_co2_apply'; call goPr
    end if

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

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

    ! get a working structure for 3d emissions
    call get_distgrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    allocate( emis3d%d3(i1:i2,j1:j2,lm(region)) ) ; emis3d%d3 = 0.0

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

    ! cycle over sectors
    do lsec = 1, numb_sectors

       if (count(used_providers.eq.sectors_def(lsec)%prov).eq.0) cycle 
          
       ! default: no additional splitting
       fraction = 1.0

       ! distinguish between 2d/3d sectors
       if( sectors_def(lsec)%f3d ) then

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

       else

          seccount2d = seccount2d + 1
          if (.not.has_data_2d(seccount2d)) cycle
          emis3d%d3 = 0.0
          ! vertically distribute according to sector
          call emission_vdist_by_sector( sectors_def(lsec)%vdisttype, 'CO2', region, co2_emis_2d(region,seccount2d), emis3d, status )
          IF_NOTOK_RETURN(status=1;deallocate(emis3d%d3))
       end if

       ! add dataset according to sector and category
       call do_add_3d( region, ico2, i1, j1, emis3d%d3, xmco2, xmco2, status, fraction )
       IF_NOTOK_RETURN(status=1)

    end do

    deallocate( emis3d%d3 )

    if(okdebug) then
       write(gol,*) 'end of emission_co2_apply'; call goPr
    end if

    ! OK 
    status = 0

  END SUBROUTINE EMISSION_CO2_APPLY
  !EOC
  
END MODULE EMISSION_CO2