!### 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_fires
  !-----------------------------------------------------------------------
  !       purpose
  !       -------
  !       perform co2 fire emissions 
  !
  !       interface
  !       ---------
  !       call emission_apply_fires
  !
  !       method
  !       ------
  !       subroutine is called from its general parent routine apply_emission
  !-----------------------------------------------------------------------
  use GO,                  only: gol, goErr, goPr, goLabel
  use Dims,                only: nregions, nlon360,nlat180
  use emission_data,       only: gridbox_per_month_to_flux
  use chem_param ,         only: xmco2, xmc
  use global_types,        only: emis_data
  use emission_common,     only: enkf_is_forward, enkf_is_inverse, EmisInputDir
  use emission_common,     only: co2_fires, co2_fires_m, co2_fires_n, flux_means_are_priors
  use emission_data,       only: do_add_2d, flux_to_gridbox_per_month
  use chem_param,          only: nmembersloc

  implicit none

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

  private

  ! public routines:
  public   :: declare_emission_co2_fires, free_emission_co2_fires, calc_emission_co2_fires, emission_apply_co2_fires

  character(len=500)              :: fire_prefix, fire_postfix, fire_varname
  logical                         :: use_fires


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

contains
  subroutine declare_emission_co2_fires
    use GO,                  only : TrcFile, Init, Done, ReadRc
    use global_data,         only : rcfile
    use dims,                only: im,jm,nregions

    implicit none
    type(TrcFile)        :: rcF
    integer              :: status
    integer              :: region, imr, jmr, n

    call Init( rcF, rcfile, status )
    call ReadRc( rcF, 'emis.use.fires', use_fires ,status ,default=.false.)
    call ReadRc( rcF, 'fire.file.prefix', fire_prefix ,status,default=trim(EmisInputDir)//'/BBguido1x1' )
    call ReadRc( rcF, 'fire.file.postfix', fire_postfix ,status ,default='.hdf')
    call ReadRc( rcF, 'fire.varname', fire_varname ,status ,default='BB_co2')
    call Done( rcF, status )

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

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

  end subroutine declare_emission_co2_fires

  subroutine calc_emission_co2_fires

    use dims,          only: nlon360, nlat180, im,jm, sec_month, newsrun, idate, newmonth, sec_year, itau, nread
    use emission_common,       only: ReadFromFile3D, get_ntimesteps
    use emission_data, only: chardate
    use toolbox,       only: coarsen_emission_2d, escape_tm
    use go_string,     only: goNum2str

    implicit none
    integer                                      :: status
    integer                                      :: region, imr, jmr, mem, ntimesteps
    integer, parameter                           :: annual_field = 0, rank2 = 2, level1 = 1
    character(len=10)                            :: cdate
    real                                         :: year2month
    real,dimension(nlon360,nlat180,248),save     :: f0
    real,dimension(nlon360,nlat180)              :: co2_fires_mem  ! flux for a given ensemble member
    integer, parameter                           :: add_field=0
    integer,save                                 :: counter 
    character(len=*), parameter                  :: rname = mname//'/calc_emission_co2_fires'

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

    if(newmonth) then
        call get_ntimesteps(ntimesteps)
        cdate=trim(goNum2str(idate(1),'(i4.4)'))//trim(goNum2str(idate(2),'(i2.2)'))
        call ReadFromFile3D(trim(fire_prefix)//trim(cdate)//trim(fire_postfix),trim(fire_varname),f0(:,:,:ntimesteps),status,netcdf4=.True.)
        if (status /= 0) call escape_tm('stop in '//rname)
    endif

    counter = ((idate(3)-1)*8 + int(idate(4)/3) + 1)
    if (idate(2) .eq. 2 .and. idate(3) .eq. 29) then
      counter = counter - 8
    endif
    co2_fires=f0(:,:,counter)

    !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 fire emissions written by a forward code are non-zero even though we only read
    !WP! and coarsen these emissions each month
    
    ! accumulate weekly average
    if(allocated(co2_fires_m)) then
       co2_fires_m=co2_fires_m+co2_fires ! mol/m2/s
       co2_fires_n=co2_fires_n+1
    endif

    if(newmonth.or.newsrun) then

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

       ! Change units, and coarsen to the regions.

       do mem=1,nmembersloc
          co2_fires_mem=co2_fires*xmc/1.e3 ! from mol/m2/s to kgC/m2/s
          call flux_to_gridbox_per_month(co2_fires_mem) ! to kgC/area/month
          call coarsen_emission_2d('co2_fires',360,180,co2_fires_mem,co2_fires_regxmem(:,mem),add_field)
       enddo
       
    endif ! if newmonth .or. newrun

  end subroutine calc_emission_co2_fires
     
  subroutine emission_apply_co2_fires(region)

    ! This routine should only be called by the forward enkf code.  It
    ! applies the fires co2 flux by member and region to the ico2_fires
    ! 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_fires


    implicit none

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

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

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

  end subroutine emission_apply_co2_fires

  subroutine free_emission_co2_fires
    use dims,        only: nregions

    implicit none
    integer               :: region,n

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

end module emission_co2_fires
