!### 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_bio
  !-----------------------------------------------------------------------
  !       purpose
  !       -------
  !       perform co2 biosphere emissions 
  !
  !       EMISSIONS
  !          interface:     call emission_apply_bio
  !          method:        subroutine is called from its general parent routine apply_emission
  !
  !
  ! Table of contents:
  !   module emission_co2_bio
  !     subroutine declare_emission_co2_bio
  !     end subroutine declare_emission_co2_bio
  !     subroutine calc_emission_co2_bio
  !     end subroutine calc_emission_co2_bio
  !     subroutine emission_apply_co2_bio(region)
  !     end subroutine emission_apply_co2_bio
  !     subroutine free_emission_co2_bio
  !     end subroutine free_emission_co2_bio
  !   end module emission_co2_bio
  !
  !-----------------------------------------------------------------------
  
  use GO,                  only: gol, goErr, goPr, goLabel
  use Dims,                only: nregions, nlon360,nlat180
  use Meteo,               only: t2m_dat,ssr_dat
  use Meteo,               only: Set
  use global_types,        only: emis_data
  use emission_data,       only: do_add_2d, flux_to_gridbox_per_month
  use emission_common,     only: enkf_is_forward,enkf_is_inverse, EmisInputDir
  use emission_common,     only: co2_bio, co2_bio_m, co2_bio_n, flux_means_are_priors, fluxmultiplier
  use chem_param,          only: nmembersloc

  implicit none

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

  character(len=500)              :: bio_prefix, bio_prefix2, bio_postfix, bio_r_varname, bio_p_varname, bio_n_varname

  private

  ! public routines:
  public   :: declare_emission_co2_bio, free_emission_co2_bio, calc_emission_co2_bio, emission_apply_co2_bio

  logical                                :: use_fires

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

contains

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

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

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

    call Init( rcF, rcfile, status )
    call ReadRc( rcF, 'emis.use.fires', use_fires ,status ,default=.false.)
    bio_prefix = trim(EmisInputDir)//'/biofireparams_'
    call ReadRc( rcF, 'bio.file.prefix', bio_prefix ,status ,default=bio_prefix)
    call ReadRc( rcF, 'bio.file.prefix2', bio_prefix2 ,status ,default=bio_prefix2)
    call ReadRc( rcF, 'bio.file.postfix', bio_postfix ,status ,default='.hdf')
    call ReadRc( rcF, 'bio.r.varname', bio_r_varname ,status ,default='rparam1')
    call ReadRc( rcF, 'bio.p.varname', bio_p_varname ,status ,default='pparam1')
    call ReadRc( rcF, 'bio.n.varname', bio_n_varname ,status ,default='nep')
    call Done( rcF, status )

    call Set(   t2m_dat(iglbsfc), status, used=.true. )
    call Set(   ssr_dat(iglbsfc), status, used=.true. )

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

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

  end subroutine declare_emission_co2_bio

  subroutine calc_emission_co2_bio
    use dims,    only        : nlon360, nlat180, im,jm, sec_month, newsrun, idate, itau, staggered, nread, newmonth, iglbsfc
    use emission_common, only        : ReadFromFile3D, ReadFromFile, get_ntimesteps
    use toolbox, only        : escape_tm
    use toolbox, only        : coarsen_emission_2d
    use emission_data, only  : chardate
    use Meteo, only          : Set
    use chem_param, only     : xmc
    use go_string, only      : goNum2str
    implicit none
    integer                                      :: status
    integer                                      :: region, imr, jmr, mem, ntimesteps
    integer, parameter                           :: annual_field = 0, rank2 = 2, level1 = 1
    real,dimension(nlon360,nlat180),save         :: p0,r0
    real,dimension(nlon360,nlat180,248),save     :: n0
    real,dimension(nlon360,nlat180)              :: co2_bio_mem ! flux for a given ensemble member
    integer, parameter                           :: add_field=0

    real,pointer,dimension(:,:)                  :: t2m,ssr
    character(len=10)                            :: cdate
    integer,save                                 :: counter  
    character(len=*), parameter        :: rname = mname//'/calc_emission_co2_bio'

    !WP! CAUTION, this statement will prevent accumulation of flux1x1 means on timescale that are not integer multiples of 
    !WP! nread. Perhaps in the future change this to hourly, or every time step for instance when interoplation of meteorology
    !WP! is used???!

    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)'))
      if(use_fires) then
       write (gol,'("[3-hourly biofireparams] filename : ",a,".")') (trim(bio_prefix)//trim(cdate)//trim(bio_postfix)) ; call goPr
          call ReadFromFile(trim(bio_prefix)//trim(cdate)//trim(bio_postfix),trim(bio_r_varname),r0,status,netcdf4=.True.)
          call ReadFromFile(trim(bio_prefix)//trim(cdate)//trim(bio_postfix),trim(bio_p_varname),p0,status,netcdf4=.True.)
          call ReadFromFile3D(trim(bio_prefix)//trim(cdate)//trim(bio_postfix),trim(bio_n_varname),n0(:,:,:ntimesteps),status,netcdf4=.True.)
       else 
          call ReadFromFile(trim(EmisInputDir)//'/bioparams_'//trim(cdate)//'.hdf','rparam1',r0,status,sds_date=(/idate(1),idate(2),1,0,0,0/))
          call ReadFromFile(trim(EmisInputDir)//'/bioparams_'//trim(cdate)//'.hdf','pparam1',p0,status,sds_date=(/idate(1),idate(2),1,0,0,0/))
       endif
      
       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_bio=n0(:,:,counter)

    ! accumulate weekly average
    if(allocated(co2_bio_m)) then
       if(flux_means_are_priors .eqv. .TRUE.) then
          co2_bio_m=co2_bio_m+co2_bio ! mol/m2/s
       else
          co2_bio_m=co2_bio_m+co2_bio*fluxmultiplier(:,:,1) ! mol/m2/s, tracer 0 on PE0 is mean
       endif
       co2_bio_n=co2_bio_n+1
    endif

    ! 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.  When parameters are used beyond simple
    ! scaling of the net flux, some piece of logic from above will
    ! have to be subsumed by this loop to get the parameter
    ! deviations.

    do mem=1,nmembersloc
       co2_bio_mem=co2_bio*fluxmultiplier(:,:,mem) ! parameter applied, mol/m2/s
       co2_bio_mem=co2_bio_mem*xmc/1.e3 ! from mol/m2/s to kgC/m2/s
       call flux_to_gridbox_per_month(co2_bio_mem) ! to kgC/area/month
       call coarsen_emission_2d('co2_bio',360,180,co2_bio_mem,co2_bio_regxmem(:,mem),add_field)
    enddo

  end subroutine calc_emission_co2_bio

  subroutine emission_apply_co2_bio(region)

    ! This routine should only be called by the forward enkf code.  It
    ! applies the bio co2 flux by member and region to the ico2_bio
    ! 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_bio, xmc, xmco2


    implicit none

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

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

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

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

  end subroutine emission_apply_co2_bio

  subroutine free_emission_co2_bio
    use dims,        only: nregions

    implicit none
    integer               :: region,n

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

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

  end subroutine free_emission_co2_bio

end module emission_co2_bio
