!### 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_co2_ff
  !-----------------------------------------------------------------------
  !       purpose
  !       -------
  !       perform co2 fossil fuel emissions 
  !
  !       EMISSIONS
  !          interface:     call emission_apply_ff
  !          method:        subroutine is called from its general parent routine apply_emission
  !
  !       COVARIANCE
  !          interface:     call covariance_co2_compute_ff
  !          method:        subroutine is called from MakePriorUncertainty in enkf_tools.F90
  !
  ! Table of contents:
  !   module emission_co2_ff
  !     subroutine declare_emission_co2_ff
  !     end subroutine declare_emission_co2_ff
  !     subroutine calc_emission_co2_ff
  !     end subroutine calc_emission_co2_ff
  !     subroutine emission_apply_co2_ff(region)
  !     end subroutine emission_apply_co2_ff
  !     subroutine free_emission_co2_ff
  !     end subroutine free_emission_co2_ff
  !   end module emission_co2_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, co2_ff_m, co2_ff_n
  use emission_data,       only: do_add_2d, flux_to_gridbox_per_month

  implicit none

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

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

  ! public routines:
  public   :: declare_emission_co2_ff, free_emission_co2_ff, calc_emission_co2_ff, emission_apply_co2_ff


  type(emis_data),allocatable,dimension(:,:),target       :: co2_ff_regxmem

contains

  subroutine declare_emission_co2_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_co2_ff'

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

       do region=1,nregions
          imr=im(region)
          jmr=jm(region)
          do n=1,nmembersloc
             allocate(co2_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_co2_ff

  subroutine calc_emission_co2_ff

    use dims,          only: nlon360, nlat180, im,jm, sec_month, newsrun, idate, newmonth, sec_year, itau, nread
    use emission_common, only: ReadFromFile3D
    use emission_data, only: chardate
    use chem_param,    only: xmc
    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)    :: co2_ff_mem ! flux for a given ensemble member
    real,dimension(:,:,:),allocatable  :: emis3d 
    integer, parameter                 :: add_field=0
    integer                            :: mem
    integer :: hid
    type(TrcFile)        ::  rcF

    character(len=*), parameter        :: rname = mname//'/calc_emission_co2_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

       ! accumulate weekly average
       !
       ! Since FF fluxes are not optimized, the following 
       ! if statement is extraneous; we want to use co2_ff
       ! regardless of the setting of flux_means_are_priors.
       !   -arj, 7 Feb 07.
       !
       !    if(flux_means_are_priors .eqv. .TRUE.) then

    !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
    if(allocated(co2_ff_m)) then
       co2_ff_m=co2_ff_m+co2_ff ! mol/m2/s
       co2_ff_n=co2_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
          co2_ff_mem=co2_ff ! mol/m2/s (this is where parameters would be applied for optimized fluxes)
          co2_ff_mem=co2_ff_mem*xmc/1.e3 ! from mol/m2/s to kgC/m2/s
          call flux_to_gridbox_per_month(co2_ff_mem) ! to kgC/area/month
          call coarsen_emission_2d('co2_ff',360,180,co2_ff_mem,co2_ff_regxmem(:,mem),add_field)
       enddo

    endif



  end subroutine calc_emission_co2_ff

  subroutine emission_apply_co2_ff(region)

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

    use dims,          only: okdebug
    use ParTools,      only: tracer_active
    use chem_param, only         : ico2_ff, xmc, xmco2


    implicit none

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

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

    if(enkf_is_forward()) then
       if(tracer_active(ico2_ff)) then
             call do_add_2d(region,ico2_ff,1,co2_ff_regxmem(region,1)%surf,xmco2,xmc)
       endif
    endif

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

  end subroutine emission_apply_co2_ff

  subroutine free_emission_co2_ff
    use dims,        only: nregions

    implicit none
    integer               :: region,n

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

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

       deallocate(co2_ff_regxmem)
    end if

  end subroutine free_emission_co2_ff

end module emission_co2_ff
