#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
#include "output.inc"  
!
!------------------------------------------------------------------------------
!                    TM5                                                      !
!------------------------------------------------------------------------------
!BOP
!
! !MODULE:      EMISSION
!
! !DESCRIPTION: wrappers around various emissions (init/declare/apply/done)
!                routines, needed for TM5 CBM4 version.
!               Also hold emissions budget variables.
! 
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION
  !
  ! !USES:
  !
  USE GO,             ONLY : gol, goErr, goPr
  USE GO,             ONLY : GO_Timer_Def, GO_Timer_End, GO_Timer_Start
  USE TM5_DISTGRID,   ONLY : dgrid, Get_DistGrid, scatter, gather
  USE dims,           ONLY : nregions, okdebug
  
  USE emission_data,  ONLY : plandr, emis2D

#ifdef with_budgets
  USE emission_data,  ONLY : budemi_dat, budemi_data, sum_emission
  USE budget_global,  ONLY : nbud_vg,nbudg
  USE chem_param,     ONLY : ntracet
#endif

  use emission_co2   , only : Emission_CO2_Init   , Emission_CO2_Done    , emission_co2_declare   , emission_co2_apply

  IMPLICIT NONE  
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  PUBLIC  :: Emission_Init     ! allocate/init budget var; call other emiss-related init
  PUBLIC  :: Emission_Done     ! gather/write final budget
  PUBLIC  :: Declare_Emission  ! allocate emiss var (new run), read emiss data (new month)
  PUBLIC  :: Emission_Apply    ! 
  !
  ! !PRIVATE DATA MEMBERS:
  !
  ! These budgets should be communicated to chemistry
#ifdef with_budgets
  REAL, DIMENSION(nbudg, nbud_vg, ntracet)         :: budemig
  REAL, DIMENSION(nbudg, nbud_vg, ntracet), PUBLIC :: budemig_all  ! for buggy MPI (see budget_global.F90)
#endif

  integer ::  itim_appl, itim_co2

  CHARACTER(len=*), PARAMETER  ::  mname = 'emission'  
  !
  ! !REVISION HISTORY: 
  !   20 Aug 2010 - A. Strunk   - Adapted to AR5 emissions + various other changes. 
  !   27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !   14 May 2014 - T. van Noije- created from chemistry version
  !                               production of CO2 from other C emissions is
  !                               currently not included
  !
  !   - NO MORE REVISION HISTORY. This is gathered from repository from now on.
  !EOP
  !----------------------------------------------------------------------

CONTAINS

  !----------------------------------------------------------------------
  !                    TM5                                              !
  !----------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_INIT
  !
  ! !DESCRIPTION: Initialise emission fields and parameters by reading 
  !               the rc-file. Allocate and initialize budget variables.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_INIT( status )
    !
    ! !USES:
    !
    use GO,              only : TrcFile, Init, Done, ReadRc
    use meteodata,       only : Set,     oro_dat
    use dims,            only : iglbsfc, nregions, lm
    use dims,            only : idate,   ndyn_max, tref
    use global_data,     only : inputdir, rcfile
    use emission_data,   only : emis_input_dir
    use emission_data,   only : emis_input_year_co2
    use emission_data,   only : LCMIP6
    use emission_data,   only : LAR5
    use emission_data,   only : emis_input_dir_cmip6
    use emission_data,   only : emis_input_dir_ar5
    use emission_read,   only : emission_read_init
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)           ::  status
    !
    ! !REVISION HISTORY: 
    !   27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER    ::  rname = mname//'/Emission_Init'
    INTEGER                        ::  region, i1, i2, j1, j2, ntim, lmr, imode
    TYPE(TrcFile)                  ::  rcF
    REAL                           ::  dtime
    integer                        ::  emis_input_year

    ! -----------------------------------
    ! read settings from rcfile
    ! -----------------------------------
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    if (okdebug) then
       write(gol,*) "EMISS-INFO - running year : ", idate(1) ; call goPr
    end if

    ! emission base year (assumption: no run overlapping more than one year)
    call ReadRc( rcF, 'input.emis.year', emis_input_year, status, default=idate(1) )
    IF_ERROR_RETURN(status=1)

    write(gol,*) 'EMISS-INFO - Anthropogenic emissions'; call goPr
    write(gol,*) 'EMISS-INFO - requested for the following years:'; call goPr

    call ReadRc( rcF, 'input.emis.year.co2', emis_input_year_co2, status, default=emis_input_year )
    IF_ERROR_RETURN(status=1)

    write(gol,*) 'EMISS-INFO - CO2:   ', emis_input_year_co2; call goPr

    ! default directory for emissions data is "standard input files" dir
    emis_input_dir=trim(inputdir)

    ! directory of each data provider
    call ReadRc( rcF, 'input.emis.dir.CMIP6', emis_input_dir_cmip6, status, default=emis_input_dir )
    IF_ERROR_RETURN(status=1)

    call ReadRc( rcF, 'input.emis.dir.AR5',   emis_input_dir_ar5,   status, default=emis_input_dir )
    IF_ERROR_RETURN(status=1)

    ! Flags
    call ReadRc( rcF, 'use_cmip6',       LCMIP6,  status, default=.false. )
    IF_ERROR_RETURN(status=1)

    call ReadRc( rcF, 'use_ar5',         LAR5,    status, default=.false. )
    IF_ERROR_RETURN(status=1) 

    ! very basic checks
    if (count((/ LCMIP6, LAR5 /)) > 1) then
       write(gol,*) 'ERROR: You use more than one ANTHROPOGENIC inventory'; call goErr
       status=1; TRACEBACK; return
    end if

    ! init providers info 
    call emission_read_init( rcF, status )

    ! used by vertical distribution:
    CALL Set( oro_dat(iglbsfc), status, used=.TRUE. )


    ! Allocate data
    ! -------------
    DO region=1,nregions

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

       ALLOCATE( plandr(region)%surf(i1:i2, j1:j2))
       ALLOCATE( emis2D(region)%surf(i1:i2, j1:j2))    

#ifdef with_budgets
       ALLOCATE( budemi_dat(region)%emi(i1:i2, j1:j2, nbud_vg, ntracet) )
       budemi_dat(region)%emi = 0.0
       sum_emission(region)   = 0.0
#endif       
    ENDDO
        
    ! Done
    ! -------------
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

    ! define timers:
    call GO_Timer_Def( itim_appl, 'emission appl', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_co2, 'emission co2', status )
    IF_NOTOK_RETURN(status=1)


    status = 0

  END SUBROUTINE EMISSION_INIT
  !EOC
  
  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_DONE
  !
  ! !DESCRIPTION: calculate and write final budgets
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_DONE( status )
    !
    ! !USES:
    !
    USE dims,          ONLY : nregions, im, jm

#ifdef with_budgets
    USE chem_param,    ONLY : ntracet, names
    USE budget_global, ONLY : budget_file_global, nbud_vg, budg_dat, nbudg, NHAB
#ifdef with_hdf4          
    USE file_hdf,      ONLY : THdfFile, TSds
    USE file_hdf,      ONLY : Init, Done, WriteAttribute, WriteData, SetDim
#endif
    USE Dims,          ONLY : region_name
    USE partools,      ONLY : isRoot, par_reduce, par_reduce_element
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)           ::  status
    !
    ! !REVISION HISTORY: 
    !   16 Jul 2010 - A. Strunk - 
    !   27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER :: rname = mname//'/Emission_Done'
    INTEGER                     :: region, i1, i2, j1, j2

#ifdef with_budgets
#ifdef with_hdf4          
    TYPE(THdfFile)                        :: io
    TYPE(TSds)                            :: sds
#endif          
    REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: collect_emissions
    INTEGER                               :: nsend,j,i,n,nzone,nzone_v
    real, dimension(nregions) :: sum_emission_all
#endif

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

#ifdef with_budgets

    ! add up contribution from all proc
    DO region = 1, nregions
       
       CALL PAR_REDUCE(sum_emission(region), 'SUM', sum_emission_all(region), status)
       IF_NOTOK_RETURN(status=1)
       
    END DO
    
    ! Write global budget of tracer #1 
    IF ( isRoot ) THEN 

       write (gol,'("EMISS-INFO - ----------------------------------------------")'); call goPr
       write (gol,'("EMISS-INFO -  Budget of tracer ",a," (kg) ")') trim(names(1))  ; call goPr 
       write (gol,'("EMISS-INFO - ----------------------------------------------")'); call goPr

       do region = 1, nregions
          write (gol,'(A,E13.6)') 'EMISS-INFO - mass emitted   :            ',sum_emission_all(region); call goPr
       enddo

#ifdef with_hdf4          
       CALL Init(io, budget_file_global, 'write', status)
       IF_NOTOK_RETURN(status=1)
       
       CALL WriteAttribute(io, 'sum_emission', sum_emission_all, status)
       IF_NOTOK_RETURN(status=1)
#endif
       budemig = 0.0
    END IF

    ! Gather budgets
    REG: DO region = 1, nregions

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

       if (isRoot) then 
          ALLOCATE(collect_emissions(im(region), jm(region), nbud_vg, ntracet))
       else
          ALLOCATE(collect_emissions(1,1,1,1) )
       end if

       CALL GATHER( dgrid(region), budemi_dat(region)%emi, collect_emissions, 0, status)
       IF_NOTOK_RETURN(status=1)
       
#ifdef with_hdf4          
       ! write Not-Horizontally-Aggregated-Budgets
       IF (isRoot.and.NHAB) THEN
          CALL Init(Sds,io, 'budemi_dat_'//region_name(region),(/im(region),jm(region),nbud_vg,ntracet/), 'real(4)', status)
          CALL SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status)
          CALL SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status)
          CALL SetDim(Sds, 2, 'nbud_vg','vertical layer', (/(j, j=1,nbud_vg)/), status)
          CALL SetDim(Sds, 3, 'ntracet','tracer number', (/(j, j=1,ntracet)/), status)
          IF_NOTOK_RETURN(status=1)
          CALL WriteData(Sds,collect_emissions,status)
          IF_NOTOK_RETURN(status=1)
          CALL Done(Sds, status)
          IF_NOTOK_RETURN(status=1)
       ENDIF
#endif
       ! horizontally aggregates budgets
       DO n=1,ntracet
       DO nzone_v=1,nbud_vg
       DO j=j1,j2
       DO i=i1,i2
          nzone = budg_dat(region)%nzong(i,j)
          budemig(nzone,nzone_v,n) = budemig(nzone,nzone_v,n) + budemi_dat(region)%emi(i,j,nzone_v,n)
       END DO
       END DO !j
       END DO !nzone_v
       END DO !nt
       
       DEALLOCATE( collect_emissions      )
       DEALLOCATE( budemi_dat(region)%emi )

    ENDDO REG


    CALL PAR_REDUCE_ELEMENT( budemig, 'SUM', budemig_all, status)
    IF_NOTOK_RETURN(status=1)
    
#ifdef with_hdf4          
    ! Write horizontally aggregated budget
    IF ( isRoot ) THEN 
       CALL Init(Sds,io, 'budemi',(/nbudg,nbud_vg,ntracet/), 'real(8)', status)
       IF_NOTOK_RETURN(status=1)
       CALL SetDim(Sds, 0, 'nbudg','horizontal region', (/(j, j=1,nbudg)/), status)
       CALL SetDim(Sds, 1, 'nbud_vg','vertical layer', (/(j, j=1,nbud_vg)/), status)
       CALL SetDim(Sds, 2, 'ntracet','tracer number', (/(j, j=1,ntracet)/), status)
       IF_NOTOK_RETURN(status=1)
       CALL WriteData(Sds,budemig_all,status)
       IF_NOTOK_RETURN(status=1)
       CALL Done(Sds, status)
       IF_NOTOK_RETURN(status=1)
       CALL Done(io, status)
       IF_NOTOK_RETURN(status=1)
    ENDIF
#endif
    
#endif /* BUDGETS */

    ! call other emission_*_done routines
    CALL FREE_EMISSION(status)
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0

  END SUBROUTINE EMISSION_DONE
  !EOC

  !---------------------------------------------------------------------------
  !                    TM5                                                   !
  !---------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DECLARE_EMISSION
  !
  ! !DESCRIPTION: Called at run start (init/allocate emiss data) and then at
  !                beginning of every month to just read in data.
  !               Called from SS_MONTHLY_UPDATE.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DECLARE_EMISSION( status )
    !
    ! !USES:
    !
    USE Grid,            ONLY : FillGrid
    USE MDF,             ONLY : MDF_Open, MDF_NETCDF, MDF_READ, MDF_Inq_VarID, MDF_Get_Var, MDF_Close
    USE dims,            ONLY : im, jm, lm, newsrun
    USE dims,            ONLY : nregions, iglbsfc, nlat180, nlon360
    USE chem_param        
    USE partools,        ONLY : isRoot
    USE global_data,     ONLY : emis_data
    USE meteodata,       ONLY : global_lli
    ! AR5/EDGAR4
    use emission_data,   only : emis_input_dir
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)      ::  status
    !
    ! !REVISION HISTORY: 
    !   16 Jul 2010 - A. Strunk - Adapted to revised emission_*.F90 routines
    !   27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !   - anything that is done only if it's a newrun, and that does not require meteo data, should go in INIT
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER          :: rname = mname//'/declare_emission'
    INTEGER                              :: region, imode, hid, varid
    REAL, DIMENSION(:,:), ALLOCATABLE    :: pland
    type(emis_data), dimension(nregions) :: wrk
    
    ! ---------------------------------------------------------------
    ! ** land fraction
    ! ---------------------------------------------------------------
    IF( newsrun ) THEN

       if(isRoot)then
          ALLOCATE( pland(nlon360,nlat180) )
          DO region=1,nregions
             allocate( wrk(region)%surf(im(region),jm(region)) )
          end DO
       else
          ALLOCATE( pland(1,1) )
          DO region=1,nregions
             allocate( wrk(region)%surf(1,1))
          end DO
       end if

       if (isRoot)then
          CALL MDF_Open( TRIM(emis_input_dir)//'/land/landfraction.nc4', MDF_NETCDF, MDF_READ, hid, status )
          IF_NOTOK_RETURN(status=1)
          CALL MDF_Inq_VarID( hid, 'LANDFRACTION', varid, status ) 
          IF_NOTOK_RETURN(status=1)
          CALL MDF_Get_Var( hid, varid, pland, status )
          IF_NOTOK_RETURN(status=1)
          CALL MDF_Close( hid, status )
          IF_NOTOK_RETURN(status=1)

          ! coarsen or distribute to zoom regions:
          DO region = 1, nregions
             ! convert grid:
             CALL FillGrid( global_lli(region),  'n', wrk(region)%surf,                     &
                            global_lli(iglbsfc), 'n', pland,            'area-aver', status )
             IF_NOTOK_RETURN(status=1)
          END DO
       end if

       DO region = 1, nregions
          call scatter( dgrid(region), plandr(region)%surf, wrk(region)%surf, 0, status)
          IF_NOTOK_RETURN(status=1)
          DEALLOCATE( wrk(region)%surf )
       END DO
       
       DEALLOCATE( pland )       

    ENDIF
    

    ! ---------------------------------------------------------------
    ! ** init each constituent
    ! ---------------------------------------------------------------

    ! ** 1st time, initialise emissions

    IF ( newsrun ) THEN
       CALL Emission_CO2_Init( status )
       IF_NOTOK_RETURN(status=1)
    END IF
    

    ! ** every month, read and re-grid
    
    CALL emission_co2_declare( status )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0

  END SUBROUTINE DECLARE_EMISSION
  !EOC


  !-------------------------------------------------------------------
  !                    TM5                                           !
  !-------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_APPLY
  !
  ! !DESCRIPTION: Call emission_apply methods of constituent modules. 
  !               --> add current emissions to tracers array. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_APPLY( region, status )
    !
    ! !USES:
    !
    USE chem_param
    !
    ! !INPUT PARAMETERS:
    !
    INTEGER, INTENT(in)        ::  region
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)       ::  status
    !
    ! !REVISION HISTORY: 
    !   16 Jul 2010 - A. Strunk   - Adapted to revised emission_*.F90 routines
    !   27 Mar 2012 - P. Le Sager - cleanup for lat-lon mpi decomposition
    !
    !EOP
    !-----------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER ::  rname = mname//'/emission_apply'

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

    ! start timing:
    call GO_Timer_Start( itim_appl, status )
    IF_NOTOK_RETURN(status=1)

    IF (okdebug) then
       WRITE(gol,*) 'start of emission_apply for region:',region ; call goPr
    END IF

    ! CO2 emissions 
    call GO_Timer_Start( itim_co2, status )
    IF_NOTOK_RETURN(status=1)
    CALL emission_co2_apply( region, status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_End( itim_co2, status )
    IF_NOTOK_RETURN(status=1)

    IF(okdebug) then
       WRITE(gol,*) 'End of adding emission '; call goPr
    END IF

    ! end timing:
    call GO_Timer_End( itim_appl, status )
    IF_NOTOK_RETURN(status=1)
    ! ok
    status = 0

  END SUBROUTINE EMISSION_APPLY
  !EOC

  !----------------------------------------------------------------------------
  !                    TM5                                                    !
  !----------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    FREE_EMISSION
  !
  ! !DESCRIPTION: Deallocate space needed to handle the emissions  by calling 
  !               *done methods of constituents' modules.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE FREE_EMISSION( status )
    !
    ! !USES:
    !
    USE dims,          ONLY : nregions
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !   16 Jul 2010 - A. Strunk   - Adapted to revised emission_*.F90 routines
    !   27 Mar 2012 - P. Le Sager - Adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !---------------------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER ::  rname = mname//'/free_emission'

    INTEGER                     ::  region, imode

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

    DO region = 1, nregions

       DEALLOCATE(plandr(region)%surf)
       DEALLOCATE(emis2D(region)%surf)
    ENDDO

    CALL Emission_CO2_Done( status )
    IF_NOTOK_RETURN(status=1)

    ! done
    status = 0

  END SUBROUTINE FREE_EMISSION
  !EOC


END MODULE EMISSION