!### macro's #####################################################
!
#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"
!
!#################################################################

module emission
  !BOP
  !-----------------------------------------------------------------------
  !       purpose
  !       -------
  !       perform emissions needed for TM5 CBM4 version
  !
  !       interface
  !       ---------
  !       call declare_emission
  !       call emission_apply
  !       call free_emission
  !
  !       method
  !       ------
  !       subroutine declare_emission is called from trace0
  !       subroutine emission_apply   is called from source1
  !       subroutine free_emission    is called from trace_end
  !-----------------------------------------------------------------------

  ! Table of Contents:
  !
  !  module emission
  !    subroutine Emis_Init(status)
  !    end subroutine Emis_Init
  !    subroutine emission_apply(region)
  !    end subroutine emission_apply
  !    subroutine Emis_Done(budget_file,status)
  !    end subroutine Emis_Done
  !  end module emission
  !EOP
  use GO,                 only: goLoCase,gol, goErr, goPr, goLabel
#ifdef with_budgets
  use emission_data,      only: budemi_dat
#endif
  use emission_data,      only: plandr, emis2D
  use emission_data,      only: sum_emission
  use emission_common,    only: declare_emission_averagers, free_emission_averagers, calc_emission_averagers, save_emission_averagers
  use emission_common,    only: EmisInputDir, ParamInputFile, C13InputDir,GetParameterMaps
  use emission_co2_ff,    only: declare_emission_co2_ff, free_emission_co2_ff, calc_emission_co2_ff, emission_apply_co2_ff
  use emission_co2_ocean, only: declare_emission_co2_ocean, free_emission_co2_ocean, calc_emission_co2_ocean, emission_apply_co2_ocean
  use emission_co2_bio,   only: declare_emission_co2_bio, free_emission_co2_bio, calc_emission_co2_bio, emission_apply_co2_bio
  use emission_co2_fires, only: declare_emission_co2_fires, free_emission_co2_fires, calc_emission_co2_fires, emission_apply_co2_fires
  use emission_co2_all,   only: declare_emission_co2_all, free_emission_co2_all, emission_apply_co2_all, calc_emission_co2_all
  use emission_c13_bio,   only: declare_emission_c13_bio, free_emission_c13_bio, calc_emission_c13_bio, emission_apply_c13_bio
  use emission_c13_ocean,   only: declare_emission_c13_ocean, free_emission_c13_ocean, calc_emission_c13_ocean, emission_apply_c13_ocean
  use emission_c13_ff,   only: declare_emission_c13_ff, free_emission_c13_ff, calc_emission_c13_ff, emission_apply_c13_ff
  use emission_c13_fires,   only: declare_emission_c13_fires, free_emission_c13_fires, calc_emission_c13_fires, emission_apply_c13_fires
  use emission_c13_biodis,   only: declare_emission_c13_biodis, free_emission_c13_biodis, calc_emission_c13_biodis, emission_apply_c13_biodis
  use emission_c13_oceandis,   only: declare_emission_c13_oceandis, free_emission_c13_oceandis, calc_emission_c13_oceandis, emission_apply_c13_oceandis
  use emission_c13_all,   only: declare_emission_c13_all, free_emission_c13_all, emission_apply_c13_all, calc_emission_c13_all

  private

  ! The default inverse behavior is to use emission_averagers to
  ! accumulate weekly average fluxes, to be written to the 
  ! savestate.hdf file.  If the rc file specifies flux1x1 output,
  ! then the averaging window may be different, and the output
  ! is to flux1x1*hdf files.  However, both facilities use
  ! the emission_averagers and they cannot both be turned on.  If
  ! flux1x1 output is requested, this logical will be .true. and
  ! the weekly average output to savestate will be suppressed.
  ! Note that we cannot use value of flux1x1_data from user_output
  ! to determine this override, because it is manipulated in
  ! main__inv.F90.
  logical :: flux1x1_override = .False.

  ! public routines:
  public   :: Emis_Init, emission_apply, Emis_Done

  ! --- const --------------------------------------

  character(len=*), parameter  ::  mname = 'emission'



contains

  subroutine Emis_Init(status)
    !BOP
    ! !ROUTINE: Emis_Init
    ! !DESCRIPTION:
    !
    !       purpose
    !       -------
    !       allocate space needed to handle the emissions
    !
    !       interface
    !       ---------
    !       call declare_emission is called from trace0 (start and newmonth)
    !
    !       method
    !       ------
    !       dynamic allocation, depending on dimensions problem
    !       opening of daily input files
    !       reading monthly/yearly emissions
    !       
    !---------------------------------------------------------------------   
    !EOP
    use dims,                 only: im,jm,lm,newsrun,nregions, nlon360, nlat180
    use chem_param,           only: names,ntracet, ico2_ff, ico2_fires, ico2_oce, ico2_bio, nmembersloc
#ifdef with_budgets
    use budget_global,      only: nbud_vg
#endif
    use emission_common,      only: flux_means_are_priors, enkf_set_inverse, fluxmultiplier,fluxmultiplier_biodis,fluxmultiplier_epsbio
    use GO,                   only: goLoCase,TrcFile, Init, Done, ReadRc
    use global_data,          only: rcfile
    use ParTools,             only: ntracetloc,ntracet_ar,myid

    implicit none
    ! --- in/out --------------------------------

    integer, intent(out)           ::  status

    ! --- const ------------------------------

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

    ! local
    integer                   :: region, imr, jmr, lmr,offsetn

    type(TrcFile)        :: rcF

    !WP! Set the run mode to inverse, use local tracer numbers as local ensemble
    !WP! size

    call enkf_set_inverse
    !print*,'enkf_set_inverse',enkf_set_inverse(),'emission_inv'
    nmembersloc = ntracetloc !/2
  !  nmembersloc = ntracetloc/2

    write(gol,*) 'CarbonTracker run mode set to INVERSE' ; call goPr

    allocate(fluxmultiplier(nlon360,nlat180,nmembersloc))
    allocate(fluxmultiplier_biodis(nlon360,nlat180,nmembersloc))
    allocate(fluxmultiplier_epsbio(nlon360,nlat180,nmembersloc))

    fluxmultiplier = 1.0
    fluxmultiplier_biodis = 1.0
    fluxmultiplier_epsbio = 1.0

    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    ! This ReadRC call commented out but not removed in case someone cares to reinstate it later.
    ! Today, changed enkf_set_inverse routine (in emission_common.F90) to set the default
    ! value of flux_means_are_priors.  For the inverse code, the default value in .TRUE.;
    ! posterior fluxes are available from a follow-on forward run or by scaling the prior
    ! fluxes by the parameter values (presuming that post = prior * scalefac).
    !                                                        -Andy Jacobson, 29 Sep 2007.
    !
    ! Reinstated!  15 Nov 2008, ARJ.  Defaults work just fine, no need to change the rc file.
    ! If you want the inverse run to have posterior fluxes in its savestate output, then
    ! set flux_means_are_priors : F in the rc file.  Otherwise, if that key is missing in the
    ! rc file, or if it is set to T, then the "standard" scheme will prevail.
    call ReadRc( rcF, 'flux_means_are_priors', flux_means_are_priors ,status ,default=.True.)
    IF_NOTOK_RETURN(status=1)

    call ReadRc( rcF, 'emis.input.dir', EmisInputDir ,status ,default='dir.not.specified')
    IF_NOTOK_RETURN(status=1)

    call ReadRc( rcF, 'c13.input.dir', C13InputDir ,status ,default='dir.not.specified')
    IF_NOTOK_RETURN(status=1)

    call ReadRc( rcF, 'ct.params.input.file', ParamInputFile, status ,default='./')
    IF_NOTOK_RETURN(status=1)

    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

    !WP! Get the values for the fluxmultipliers from file

    call GetParameterMaps(status)
    IF_NOTOK_RETURN(status=1)

    ! (maybe) initialize the average flux accumulators for weekly output
    if(.not. flux1x1_override) then
       call declare_emission_averagers
    endif

#ifdef with_budgets
    ! on all PEs to be made available:
    do region=1,nregions
       imr = im(region) ;jmr = jm(region) ; lmr = lm(region)
       allocate(budemi_dat(region)%emi(im(region),jm(region),nbud_vg,ntracet))
       budemi_dat(region)%emi     = 0.0
       sum_emission(region)       = 0.0
    enddo
#endif


    !offsetn = sum(ntracet_ar(0:myid-1))
    !if (goLoCase(trim(names(offsetn+1)))=='co2'.or.goLoCase(trim(names(offsetn+1)))=='co2_bg') then


    call declare_emission_co2_ff  
    call declare_emission_co2_fires  
    call declare_emission_co2_ocean
    call declare_emission_co2_bio
    call declare_emission_co2_all
    !endif

    !if (goLoCase(trim(names(offsetn+1)))=='co2c13'.or.goLoCase(trim(names(offsetn+1)))=='co2c13_bg') then
    call declare_emission_c13_bio
    call declare_emission_c13_ocean
    call declare_emission_c13_ff
    call declare_emission_c13_fires
    call declare_emission_c13_biodis
    call declare_emission_c13_oceandis
    call declare_emission_c13_all
    !endif

    status=0

  end subroutine Emis_Init

  ! ---------------------------------

  subroutine emission_apply(region, status)
    use Dims,          only: okdebug
    use chem_param,         only: ntracet, ico2_ff, ico2_fires, ico2_oce, ico2_bio,names
    use ParTools,      only: ntracetloc,ntracet_ar, tracer_active, myid, root_t

    implicit none
    ! --- const ---------------------------------------

    character(len=*), parameter ::  rname = mname//'/emission_apply'
    integer                    :: offsetn 
    integer, intent(in)        :: region
    integer, intent(out)       ::  status

    if(okdebug) write(*,*) 'start of emission_apply'
    if(ntracetloc == 0) return      !if no tracer on this processor do nothing


    if(region == 1) then  !WP! only needed once, on 1x1 grid, coarsened array used after that
       call calc_emission_co2_ff
       call calc_emission_co2_bio
       call calc_emission_co2_fires
       call calc_emission_co2_ocean
       call calc_emission_c13_bio
       call calc_emission_c13_ocean
       call calc_emission_c13_ff
       call calc_emission_c13_fires
       call calc_emission_c13_biodis
       call calc_emission_c13_oceandis
    endif
    
    call calc_emission_co2_all
    call emission_apply_co2_all(region)
    call calc_emission_c13_all
    call emission_apply_c13_all(region)

    
    if(okdebug) write(*,*) 'End of adding emission ' 

  !ok
  status = 0

  end subroutine emission_apply

  subroutine Emis_Done(budget_file,status)
    !
    !       purpose
    !       -------
    !       deallocate space needed to handle the emissions
    !
    !       interface
    !       ---------
    !       call free_emission is called from trace_end
    !
    !       method
    !       ------
    !       collecting budgets on coarser budget regions
    !       writing budget file
    !       deallocation of memory
    !       closing of input files
    !---------------------------------------------------------------------   
    use dims,          only : im,jm, nregions
    use chem_param,      only: names,ntracet, ico2_ff, ico2_fires, ico2_oce, ico2_bio
    use emission_common, only: fluxmultiplier,fluxmultiplier_biodis,    fluxmultiplier_epsbio
#ifdef with_budgets
    use budget_global, only : nbudg, nbud_vg, budg_dat, apply_budget_global, budget_hdf_names_unique
#endif
    use io_hdf,        only : DFACC_WRITE, io_write
#ifdef MPI
    use mpi_const, only: ntracet_ar, lmloc, tracer_active, ntracetloc, myid, root, my_real, mpi_sum, com_trac, ierr, root_t
#else
    use ParTools, only: lmloc, ntracet_ar,tracer_active, myid, root, ntracetloc, root_t
#endif    

    implicit none
    ! --- in/out --------------------------------

    integer, intent(out)           ::  status
    character(len=*)               :: budget_file

    ! --- const ------------------------------

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

    ! local
    integer                                  :: region, imr, jmr
    integer                                  :: i,j,n, nzone,offsetn
#ifdef with_budgets
    real, dimension(nbudg, nbud_vg, ntracet) :: budemig
    real, dimension(:,:,:,:),allocatable     :: budget4d
#endif

    integer                                  :: ifail,sfend, nsend, io, sfstart

#ifdef with_budgets
    if(apply_budget_global) then

       !handle budget emission:

       if(myid == root) then
          io = sfstart(budget_file, DFACC_WRITE)
       endif

       budemig(:,:,:) = 0.0
       do region = 1, nregions
          allocate(budget4d(im(region),jm(region),nbud_vg,ntracet))
#ifdef MPI
          nsend = im(region)*jm(region)*nbud_vg*ntracet
          call mpi_allreduce(budemi_dat(region)%emi,budget4d, &
               nsend,my_real,mpi_sum,com_trac,ierr)
#else
          budget4d = budemi_dat(region)%emi
#endif
          call io_write(io,im(region),'im',jm(region),'jm', nbud_vg, 'nbud_vg', ntracet,'ntracet', &
               budget4d ,'budemi_dat')
          do n=1,ntracet 
             do j=1,jm(region)
                do i=1,im(region)
                   nzone = budg_dat(region)%nzong(i,j)
                   budemig(nzone,:,n) = budemig(nzone,:,n) + &
                        budget4d(i,j,:,n) 
                end do
             end do!j
          end do !
          deallocate(budget4d)
       enddo  ! regions

       if(myid == root) then
          call io_write(io,nbudg,'nbud',nbud_vg,'nbud_v', &
               ntracet,'ntracet',budemig,'budemi')
          print *, 'Close budget file:', sfend(io)
       endif

    endif ! apply_budget_global

    do region = 1, nregions
       deallocate(budemi_dat(region)%emi)
    enddo
#endif

    if(.not. flux1x1_override) then
       call calc_emission_averagers
       call save_emission_averagers(status)
       IF_NOTOK_RETURN(status=1)
    endif

    call free_emission_co2_ff  
    call free_emission_co2_fires  
    call free_emission_co2_ocean
    call free_emission_co2_bio
    call free_emission_co2_all
    !endif

    call free_emission_c13_bio
    call free_emission_c13_ocean
    call free_emission_c13_ff
    call free_emission_c13_fires
    call free_emission_c13_biodis
    call free_emission_c13_oceandis
    call free_emission_c13_all
   ! endif

    if(.not. flux1x1_override) then
       call free_emission_averagers
    endif

    deallocate(fluxmultiplier)   
    deallocate(fluxmultiplier_biodis)
    deallocate(fluxmultiplier_epsbio)
    status=0

  end subroutine Emis_Done

end module emission
