!### 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_bio
  !-----------------------------------------------------------------------
  !       purpose
  !       -------
  !       perform c13 biosphere emissions 
  !
  !       EMISSIONS
  !          interface:     call emission_apply_bio
  !          method:        subroutine is called from its general parent routine apply_emission
  !
  !
  ! Table of contents:
  !   module emission_c13_bio
  !     subroutine declare_emission_c13_bio
  !     end subroutine declare_emission_c13_bio
  !     subroutine calc_emission_c13_bio
  !     end subroutine calc_emission_c13_bio
  !     subroutine emission_apply_c13_bio(region)
  !     end subroutine emission_apply_c13_bio
  !     subroutine free_emission_c13_bio
  !     end subroutine free_emission_c13_bio
  !   end module emission_c13_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: field2d_eps,field2d_eps_m,field2d_eps_n,co2_bio,c13_bio, c13_bio_m, c13_bio_n, flux_means_are_priors, fluxmultiplier,fluxmultiplier_epsbio,c13_sig
  use emission_common,     only: c13_gpp, c13_gpp_m, c13_gpp_n
  use emission_common,     only: c13_res, c13_res_m, c13_res_n
  use emission_common,     only: field2d_eps,field2d_eps_m,field2d_eps_n,co2_gpp, co2_gpp_m, co2_gpp_n

  use chem_param,          only: nmembersloc

  implicit none

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

  character(len=500)              :: eps_prefix,bio_prefix, bio_prefix2, bio_postfix, bio_g_varname , bio_r_varname, bio_p_varname, bio_n_varname,bio_c13_n_varname,bio_w_varname

  private

  ! public routines:
  public   :: declare_emission_c13_bio, free_emission_c13_bio, calc_emission_c13_bio, emission_apply_c13_bio

  logical                                :: use_fires

  ! c13_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       :: c13_bio_regxmem

contains

  subroutine declare_emission_c13_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_c13_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='resp13')
    call ReadRc( rcF, 'bio.p.varname', bio_p_varname ,status ,default='gpp13')
    call ReadRc( rcF, 'bio.g.varname', bio_g_varname ,status ,default='gpp')
    call ReadRc( rcF, 'bio.c13.n.varname', bio_c13_n_varname ,status ,default='nep13')
    call ReadRc( rcF, 'bio.n.varname', bio_n_varname ,status ,default='nep')
    call ReadRc( rcF, 'bio.w.varname', bio_w_varname ,status ,default='wk')
    call ReadRc( rcF, 'eps.prefix', eps_prefix ,status)

    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(c13_bio_regxmem(nregions,nmembersloc))

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

  end subroutine declare_emission_c13_bio

  subroutine calc_emission_c13_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,xmc13
    use go_string, only      : goNum2str
    use rs_read, only        : da
    use partools,   only : myid
    implicit none
    integer                                      :: status,i,j,k
    integer                                      :: region, imr, jmr, mem, ntimesteps
    integer, parameter                           :: annual_field = 0, rank2 = 2, level1 = 1
    real,dimension(nlon360,nlat180,248),save         :: p0,r0,g0,c13_signature
    real,dimension(nlon360,nlat180,248),save     :: n0,n13c0,w0
    real,dimension(nlon360,nlat180),save         :: Delta,Delta2,alpha !,field2d_epsbio
    real,dimension(nlon360,nlat180)              :: c13_bio_mem, c13_gpp_mem, c13_res_mem,co2_gpp_mem ! flux for a given ensemble member
    real,dimension(4,nlon360,nlat180),save       :: field2d_epsbio
    integer, parameter                           :: add_field=0

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



    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 ReadFromFile3D(trim(bio_prefix)//trim(cdate)//trim(bio_postfix),trim(bio_r_varname),r0(:,:,:ntimesteps),status,netcdf4=.True.)
          call ReadFromFile3D(trim(bio_prefix)//trim(cdate)//trim(bio_postfix),trim(bio_p_varname),p0(:,:,:ntimesteps),status,netcdf4=.True.)
          call ReadFromFile3D(trim(bio_prefix)//trim(cdate)//trim(bio_postfix),trim(bio_n_varname),n0(:,:,:ntimesteps),status,netcdf4=.True.)
          call ReadFromFile3D(trim(bio_prefix)//trim(cdate)//trim(bio_postfix),trim(bio_g_varname),g0(:,:,:ntimesteps),status,netcdf4=.True.)
          call ReadFromFile3D(trim(bio_prefix)//trim(cdate)//trim(bio_postfix),trim(bio_w_varname),w0(:,:,:ntimesteps),status,netcdf4=.True.)
          call ReadFromFile3D(trim(bio_prefix)//trim(cdate)//trim(bio_postfix),trim(bio_c13_n_varname),n13c0(:,:,:ntimesteps),status,netcdf4=.True.)
          call ReadFromFile3D(trim(eps_prefix)//trim(cdate)//'.nc','epsilon',field2d_epsbio,status,netcdf4=.True.)



       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

    counter2 = (int((idate(3)-1)/7) + 1)
    if (idate(3) .gt. 28) then
      counter2 = 4
    endif



 
    co2_gpp=g0(:,:,counter)
    co2_bio=n0(:,:,counter)
    field2d_eps = field2d_epsbio(counter2,:,:)
    alpha = 1+field2d_eps/1000
    c13_bio = co2_bio*(0.011112*(1.+da/1000.))*alpha



    ! accumulate weekly average
    if(allocated(c13_bio_m)) then
       if(flux_means_are_priors .eqv. .TRUE.) then
          c13_bio_m=c13_bio_m+c13_bio ! mol/m2/s
          c13_gpp_m=c13_gpp_m+co2_gpp*field2d_eps ! mol/m2/s
          co2_gpp_m=co2_gpp_m+co2_gpp ! mol/m2/s
          field2d_eps_m=field2d_eps_m+field2d_eps

       else
          c13_bio_m=c13_bio_m+c13_bio*fluxmultiplier(:,:,1) ! mol/m2/s, tracer 0 on PE0 is mean
          c13_gpp_m=c13_gpp_m+co2_gpp*fluxmultiplier_epsbio(:,:,1) ! mol/m2/s, tracer 0 on PE0 is mean
          co2_gpp_m=co2_gpp_m+co2_gpp !fluxmultiplier(:,:,1) ! mol/m2/s, tracer 0 on PE0 is mean
          field2d_eps_m=field2d_eps_m+field2d_eps*fluxmultiplier_epsbio(:,:,1)

       endif
       c13_bio_n=c13_bio_n+1
       c13_gpp_n=c13_gpp_n+1
       co2_gpp_n=co2_gpp_n+1
       field2d_eps_n=field2d_eps_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
       c13_bio_mem=c13_bio*fluxmultiplier(:,:,mem) ! parameter applied, mol/m2/s
       c13_bio_mem=c13_bio_mem*xmc13/1.e3 ! from mol/m2/s to kgC/m2/s
       call flux_to_gridbox_per_month(c13_bio_mem) ! to kgC13/area/month
       call coarsen_emission_2d('c13_bio',360,180,c13_bio_mem,c13_bio_regxmem(:,mem),add_field)


    enddo

  end subroutine calc_emission_c13_bio

  subroutine emission_apply_c13_bio(region)

    ! This routine should only be called by the forward enkf code.  It
    ! applies the bio c13 flux by member and region to the ic13_bio
    ! 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_bio, ic13_gpp, ic13_res, ico2_gpp, xmc13, xm13co2,xmc,xmco2


    implicit none

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

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

    if(enkf_is_forward()) then
       if(tracer_active(ic13_bio)) then
             call do_add_2d(region,ic13_bio,1,c13_bio_regxmem(region,1)%surf,xm13co2,xmc13)

       endif
    endif

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

  end subroutine emission_apply_c13_bio

  subroutine free_emission_c13_bio
    use dims,        only: nregions

    implicit none
    integer               :: region,n

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

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

          end do
       enddo
       
       deallocate(c13_bio_regxmem)

    endif

  end subroutine free_emission_c13_bio

end module emission_c13_bio
