#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
                                                                                                                                                     
#include "tm5.inc"
!#################################################################

module emission_c13_biodis
  !-----------------------------------------------------------------------
  !       purpose
  !       -------
  !       perform c13 biosphere disequilibrium emissions 
  !
  !       EMISSIONS
  !          interface:     call emission_apply_bio
  !          method:        subroutine is called from its general parent routine apply_emission
  !
  !
  ! Table of contents:
  !   module emission_c13_biodis
  !     subroutine declare_emission_c13_biodis
  !     end subroutine declare_emission_c13_biodis
  !     subroutine calc_emission_c13_biodis
  !     end subroutine calc_emission_c13_biodis
  !     subroutine emission_apply_c13_biodis(region)
  !     end subroutine emission_apply_c13_biodis
  !     subroutine free_emission_c13_biodis
  !     end subroutine free_emission_c13_biodis
  !   end module emission_c13_biodis
  !
  !-----------------------------------------------------------------------
  
  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,gridbox_per_month_to_flux
  use emission_common,     only: enkf_is_forward,enkf_is_inverse, EmisInputDir,C13InputDir
  use emission_common,     only: c13_biodis, c13_biodis_m, c13_biodis_n, flux_means_are_priors, fluxmultiplier

  use chem_param,          only: nmembersloc

  implicit none

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

  private

  ! public routines:
  public   :: declare_emission_c13_biodis, free_emission_c13_biodis, calc_emission_c13_biodis, emission_apply_c13_biodis

  logical                                :: use_fires
  character(len=500)              :: terdis_prefix

  ! c13_biodis_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_biodis_regxmem

contains

  subroutine declare_emission_c13_biodis
    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_biodis'

    call Init( rcF, rcfile, status )
    call ReadRc( rcF, 'emis.use.fires', use_fires ,status ,default=.false.)
    call Done( rcF, status )

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

    call Init( rcF, rcfile, status )
    call ReadRc( rcF, 'terdis.prefix', terdis_prefix ,status)
    call Done( rcF ,status)


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

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

  end subroutine declare_emission_c13_biodis

  subroutine calc_emission_c13_biodis
    use dims,    only        : nlon360, nlat180, im,jm, sec_month, newsrun, idate, itau, staggered, nread, newmonth, iglbsfc
    use emission_common, only        : ReadFromFile3D,ReadFromFile
    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,xm13co2
    use go_string, only      : goNum2str
    implicit none
    integer                                      :: status
    integer                                      :: region, imr, jmr, mem
    integer, parameter                           :: annual_field = 0, rank2 = 2, level1 = 1
    real,dimension(nlon360,nlat180)              :: p0,r0
    real,dimension(nlon360,nlat180,12),save      :: field3d_dis 
    real,dimension(nlon360,nlat180),save      :: field2d_dis 
    real,dimension(nlon360,nlat180)              :: c13_biodis_mem ! flux for a given ensemble member
    integer, parameter                           :: add_field=0
    real,parameter                               :: delta_atm=-8.

    real,pointer,dimension(:,:)                  :: t2m,ssr
    character(len=10)                            :: cdate

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

    !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
       cdate=trim(goNum2str(idate(1),'(i4.4)'))//trim(goNum2str(idate(2),'(i2.2)'))
      
       !Read biosphere disequilibrium
       
       write(gol,'(a)') 'BIODIS data to be read '//trim(cdate) ; call goPr
       call ReadFromFile(trim(terdis_prefix)//trim(cdate)//'.nc','ter_dis',field2d_dis,status,netcdf4=.True.)

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

    endif

    write(gol,'(a)') 'CarbonTracker BIODIS data to be recalculated ' ; call goPr
    c13_biodis=field2d_dis(:,:)     !kg 13CO2/cell/month

    call gridbox_per_month_to_flux(c13_biodis)   ! to kg 13CO2 /m2/s
    c13_biodis=c13_biodis*1.e3/xm13co2  ! to mol/m2/s

    ! accumulate weekly average
    if(allocated(c13_biodis_m)) then
       if(flux_means_are_priors .eqv. .TRUE.) then
          c13_biodis_m=c13_biodis_m+c13_biodis ! mol/m2/s
       else
          c13_biodis_m=c13_biodis_m+c13_biodis
       endif
       c13_biodis_n=c13_biodis_n+1
    endif

    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_biodis_mem=c13_biodis  !*fluxmultiplier(:,:,mem) ! parameter applied, mol/m2/s          
       c13_biodis_mem=c13_biodis_mem*xmc13/1.e3 ! from mol/m2/s to kgC13/m2/s
       call flux_to_gridbox_per_month(c13_biodis_mem) ! to kgC13/area/month
       call coarsen_emission_2d('c13_biodis',360,180,c13_biodis_mem,c13_biodis_regxmem(:,mem),add_field)
    enddo

  end subroutine calc_emission_c13_biodis

  subroutine emission_apply_c13_biodis(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: myid,tracer_active
    use chem_param, only         : ic13_biodis, xmc13, xm13co2


    implicit none

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

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

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

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

  end subroutine emission_apply_c13_biodis

  subroutine free_emission_c13_biodis
    use dims,        only: nregions

    implicit none
    integer               :: region,n

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

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

  end subroutine free_emission_c13_biodis

end module emission_c13_biodis
