!### macros #####################################################
#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_c13_ff
  !-----------------------------------------------------------------------
  !       purpose
  !       -------
  !       perform c13 fossil fuel emissions 
  !
  !       EMISSIONS
  !          interface:     call emission_apply_ff
  !          method:        subroutine is called from its general parent routine apply_emission
  !
  !       COVARIANCE
  !          interface:     call covariance_c13_compute_ff
  !          method:        subroutine is called from MakePriorUncertainty in enkf_tools.F90
  !
  ! Table of contents:
  !   module emission_c13_ff
  !     subroutine declare_emission_c13_ff
  !     end subroutine declare_emission_c13_ff
  !     subroutine calc_emission_c13_ff
  !     end subroutine calc_emission_c13_ff
  !     subroutine emission_apply_c13_ff(region)
  !     end subroutine emission_apply_c13_ff
  !     subroutine free_emission_c13_ff
  !     end subroutine free_emission_c13_ff
  !   end module emission_c13_ff
  !
  !-----------------------------------------------------------------------

  use GO,                  only: gol, goErr, goPr, goLabel
  use Dims,                only: nregions, nlon360,nlat180
  use emission_data,       only: gridbox_per_month_to_flux
  use global_types,        only: emis_data
  use chem_param,          only: nmembersloc
  use emission_common,     only: enkf_is_forward, enkf_is_inverse, EmisInputDir
  use emission_common,     only: co2_ff, c13_ff, c13_ff_m, c13_ff_n
  use emission_data,       only: do_add_2d, flux_to_gridbox_per_month

  implicit none

  ! --- const -------------------------
  character(len=*), parameter  ::  mname = 'emission_c13_ff'

  private
  character(len=2000)              :: ff_prefix, ff_varname, ff_postfix


  ! public routines:
  public   :: declare_emission_c13_ff, free_emission_c13_ff, calc_emission_c13_ff, emission_apply_c13_ff

  ! c13_ff_regxmem is the result of the calculate routine.
  ! It holds the coarsened flux per ensemble member and region.
  type(emis_data),allocatable,dimension(:,:),target       :: c13_ff_regxmem

contains

  subroutine declare_emission_c13_ff

    use dims,       only: im,jm,nregions, idate
    use global_data,   only: rcfile
    use GO,            only: TrcFile, Init, Done, ReadRc
    implicit none
    integer                            :: status
    integer                            :: region, imr, jmr, n
    type(TrcFile)        :: rcF

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

    if(enkf_is_forward()) then
       allocate(c13_ff_regxmem(nregions, nmembersloc))

       do region=1,nregions
          imr=im(region)
          jmr=jm(region)
          do n=1,nmembersloc
             allocate(c13_ff_regxmem(region,n)%surf(imr,jmr))
          end do
       enddo
    end if

    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'ff.file.prefix', ff_prefix ,status)
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'ff.file.postfix', ff_postfix ,status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'ff.varname', ff_varname ,status)
    IF_NOTOK_RETURN(status=1)
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

  end subroutine declare_emission_c13_ff

  subroutine calc_emission_c13_ff

    use dims,          only: nlon360, nlat180, im,jm, sec_month, newsrun, idate, newmonth, sec_year, itau, nread
    use emission_common,       only: ReadFromFile,ReadFromFile3D
    use emission_data, only: chardate
    use chem_param,    only: xmc,xmc13
    use toolbox,       only: coarsen_emission_2d, escape_tm
    use go_string,     only: goNum2str
    use global_data,   only: rcfile
    use GO,            only: TrcFile, Init, Done, ReadRc


    implicit none
    integer                            :: status
    integer                            :: region, imr, jmr
    integer, parameter                 :: annual_field = 0, rank2 = 2, level1 = 1
    character(len=10)                  :: cdate
    real                               :: year2month
    real,dimension(nlon360,nlat180)    :: c13_ff_mem ! flux for a given ensemble member
    real,dimension(:,:,:),allocatable  :: emis3d 
    integer, parameter                 :: add_field=0
    integer                            :: mem
    type(TrcFile)        ::  rcF
    real,parameter     :: deltaff=-30.

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

    if(mod(itau,nread) /= 0) return  ! only every nread hours

    if(newmonth) then

        allocate(emis3d(nlon360,nlat180,12))
        cdate=trim(goNum2str(idate(1),'(i4.4)'))
        call ReadFromFile3D(trim(ff_prefix)//trim(cdate)//trim(ff_postfix),trim(ff_varname),emis3d,status,netcdf4=.True.)
        co2_ff = emis3d(:,:,idate(2))*1e3/xmc ! from kgC/m2/s to mol/m2/s

        deallocate(emis3d)
        if (status /= 0) call escape_tm('stop in '//rname)
     endif



    !WP! this block has been moved to be outside the (newmonth) flag to ensure that the averaging is done at each call.
    !WP! This ensures that the n-hourly average fossil fuel emissions written by a forward code are non-zero even though we only read
    !WP! and coarsen these emissions each month

       c13_ff=co2_ff*(0.011112*(1.+deltaff/1000.))

    if(allocated(c13_ff_m)) then
       c13_ff_m=c13_ff_m+c13_ff ! mol/m2/s
       c13_ff_n=c13_ff_n+1
    endif

  
    if(newmonth) then

       ! The rest only useful to the forward code
       if(enkf_is_inverse()) return

       ! Apply net flux scaling by ensemble member, change units, and
       ! coarsen to the regions.

       do mem=1,nmembersloc
          c13_ff_mem=c13_ff ! mol/m2/s (this is where parameters would be applied for optimized fluxes)
          c13_ff_mem=c13_ff_mem*xmc13/1.e3 ! from mol/m2/s to kgC/m2/s
          call flux_to_gridbox_per_month(c13_ff_mem) ! to kgC/area/month
          call coarsen_emission_2d('c13_ff',360,180,c13_ff_mem,c13_ff_regxmem(:,mem),add_field)
       enddo

    endif



  end subroutine calc_emission_c13_ff

  subroutine emission_apply_c13_ff(region)

    ! This routine should only be called by the forward enkf code.  It
    ! applies the ff c13 flux by member and region to the ic13ff
    ! tracer.  In the inverse code, the calculated flux from each
    ! module is applied to the (single) c13 tracer.

    use dims,          only: okdebug
    use ParTools,      only: tracer_active
    use chem_param, only         : ic13_ff, xmc13, xm13co2


    implicit none

    integer, intent(in)        :: region
    integer                    :: mem

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

    if(enkf_is_forward()) then
       if(tracer_active(ic13_ff)) then
             call do_add_2d(region,ic13_ff,1,c13_ff_regxmem(region,1)%surf,xm13co2,xmc13)
       endif
    endif

    if(okdebug) write(*,*) 'end of emission_apply_c13_ff'

  end subroutine emission_apply_c13_ff

  subroutine free_emission_c13_ff
    use dims,        only: nregions

    implicit none
    integer               :: region,n

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

    if(enkf_is_forward()) then
       do region=1,nregions
          do n=1,nmembersloc
             deallocate(c13_ff_regxmem(region,n)%surf)
          end do
       enddo

       deallocate(c13_ff_regxmem)
    end if

  end subroutine free_emission_c13_ff

end module emission_c13_ff
