!### 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_common

  ! While the averagers could be put into each flux module, managing
  ! these routines is easier when all the repetitive code is gathered
  ! here in one place.

  use GO,                  only: gol, goErr, goPr, goLabel
  use Dims,                only: nlon360,nlat180

  public   :: co2_ocean, co2_fires, co2_ff, co2_bio, co2_all, co2_gpp, field2d_eps
  public   :: co2_ocean_m, co2_fires_m, co2_ff_m, co2_bio_m, co2_all_m, co2_gpp_m, field2d_eps_m
  public   :: co2_ocean_n, co2_fires_n, co2_ff_n, co2_bio_n, co2_all_n, co2_gpp_n, field2d_eps_n
  public   :: c13_bio, c13_sig,c13_gpp, c13_res, c13_ocean, c13_ff, c13_fires, c13_biodis, c13_oceandis, c13_all
  public   :: c13_bio_m, c13_gpp_m, c13_res_m, c13_ocean_m, c13_ff_m, c13_fires_m, c13_biodis_m, c13_oceandis_m, c13_all_m
  public   :: c13_bio_n, c13_gpp_n, c13_res_n, c13_ocean_n, c13_ff_n, c13_fires_n, c13_biodis_n, c13_oceandis_n, c13_all_n  
  public   :: declare_emission_averagers,free_emission_averagers,reset_emission_averagers, calc_emission_averagers, save_emission_averagers
  public   :: enkf_is_inverse, enkf_is_forward, enkf_set_inverse, enkf_set_forward
  public   :: flux_means_are_priors, EmisInputDir,C13InputDir
  public   :: ReadFromFile,ReadFromFile3D
  public   :: fluxmultiplier, fluxmultiplier_biodis,fluxmultiplier_epsbio,GetParameterMaps, ParamInputFile

  character(len=*), parameter   ::  mname = 'emission_common'
  private:: mname

  real,dimension(nlon360,nlat180),target :: co2_ocean
  real,dimension(nlon360,nlat180),target :: co2_fires
  real,dimension(nlon360,nlat180),target :: co2_bio
  real,dimension(nlon360,nlat180),target :: co2_ff
  real,dimension(nlon360,nlat180),target :: co2_all
  real,dimension(nlon360,nlat180),target :: co2_gpp
  real,dimension(nlon360,nlat180),target :: c13_all
  real,dimension(nlon360,nlat180),target :: c13_bio
  real,dimension(nlon360,nlat180),target :: c13_gpp
  real,dimension(nlon360,nlat180),target :: c13_res
  real,dimension(nlon360,nlat180),target :: c13_ocean    
  real,dimension(nlon360,nlat180),target :: c13_ff
  real,dimension(nlon360,nlat180),target :: c13_fires
  real,dimension(nlon360,nlat180),target :: c13_biodis
  real,dimension(nlon360,nlat180),target :: c13_oceandis
  real,dimension(nlon360,nlat180),target :: field2d_eps
  real,dimension(nlon360,nlat180),target :: c13_sig
  ! weekly average flux accumulators.  "_m" is the flux field itself,
  ! "_n" is the number of accumulations.

  real,allocatable,dimension(:,:) :: co2_ocean_m
  real,allocatable,dimension(:,:) :: co2_fires_m
  real,allocatable,dimension(:,:) :: co2_bio_m
  real,allocatable,dimension(:,:) :: co2_ff_m
  real,allocatable,dimension(:,:) :: co2_all_m
  real,allocatable,dimension(:,:) :: co2_gpp_m
  real,allocatable,dimension(:,:) :: c13_all_m
  real,allocatable,dimension(:,:) :: c13_bio_m
  real,allocatable,dimension(:,:) :: c13_gpp_m
  real,allocatable,dimension(:,:) :: c13_res_m
  real,allocatable,dimension(:,:) :: c13_ocean_m
  real,allocatable,dimension(:,:) :: c13_ff_m
  real,allocatable,dimension(:,:) :: c13_fires_m
  real,allocatable,dimension(:,:) :: c13_biodis_m
  real,allocatable,dimension(:,:) :: c13_oceandis_m
  real,allocatable,dimension(:,:) :: field2d_eps_m

  !WP! This will hold multiplier values for each member of the inverse run, or
  !WP! for the mean of the forward run

  real,allocatable,dimension(:,:,:) :: fluxmultiplier
  real,allocatable,dimension(:,:,:) :: fluxmultiplier_biodis
  real,allocatable,dimension(:,:,:) :: fluxmultiplier_epsbio

  integer                                :: co2_ocean_n
  integer                                :: co2_fires_n
  integer                                :: co2_bio_n
  integer                                :: co2_ff_n
  integer                                :: co2_all_n
  integer                                :: co2_gpp_n
  integer                                :: c13_all_n
  integer                                :: c13_bio_n
  integer                                :: c13_gpp_n
  integer                                :: c13_res_n
  integer                                :: c13_ocean_n  
  integer                                :: c13_ff_n
  integer                                :: c13_fires_n
  integer                                :: c13_biodis_n
  integer                                :: c13_oceandis_n
  integer                                :: field2d_eps_n

  ! Are we running in inverse or forward mode?  Set enkf mode in Emis_Init().
  integer,parameter             :: enkf_mode_unknown = 0
  integer,parameter             :: enkf_mode_inverse = 1
  integer,parameter             :: enkf_mode_forward = 2
  integer                       :: enkf_mode = enkf_mode_unknown

  character(len=256)            :: EmisInputDir
  character(len=256)            :: ParamInputFile
  character(len=256)            :: C13InputDir

  character(len=1024)			:: FluxOutFile = ''

  ! Weekly mean fluxes can either be means of the raw priors, or means
  ! of the optimized fluxes.  Default is priors.  This was the
  ! behavior in previous versions, since one needs to verify what the
  ! flux modules are doing, and it is possible to retrieve the scaled
  ! fluxes in postprocessing.  This switch removes the need for that
  ! postprocessing, and may be useful in a future with a more complex
  ! parameter model.
  !
  ! flux_means_are_priors = .TRUE.  -> report weekly means of prior fluxes
  !
  ! flux_means_are_priors = .FALSE.  -> report weekly means of scaled/parameterized/optimized fluxes
  !
  logical                        :: flux_means_are_priors = .TRUE.


contains

    logical function enkf_is_inverse()
      enkf_is_inverse = .FALSE.
      if(enkf_mode .EQ. enkf_mode_inverse) then
         enkf_is_inverse = .TRUE.
      endif
    end function enkf_is_inverse

    logical function enkf_is_forward()
      enkf_is_forward = .FALSE.
      if(enkf_mode .EQ. enkf_mode_forward) then
         enkf_is_forward = .TRUE.
      endif
    end function enkf_is_forward

    ! Comment on flux_means_are_priors  (Andy Jacobson, 29 Sep 2007)
    !
    ! Wouter and I decided that the inverse run will output weekly-average
    ! prior fluxes by default, and that a follow-on forward run (which is 
    ! already pretty much a requirement for various output) will write out 
    ! the posterior fluxes.  Thus the default values for flux_means_as_priors
    ! below.  Note that in the inv and fwd versions of emission.F90, there
    ! is a commented-out ReadRC call that could be used to change this
    ! hard-coded default behavior.  Use of old RC files with flux_means_are_priors
    ! lines may cause some confusion, but the variables written to the 
    ! output file actually are given explicit names so it will be tough
    ! to make a real mistake.
    !
    ! It is convenient to look at posterior fluxes while an inverse run is
    ! underway.  For this, we have two suggestions.  First, the current model
    ! for fluxes is that the posterior is a scaled version of the prior.  The
    ! scale factors are the EnKF parameters, and are written to the savestate
    ! file along with the prior fluxes in the inverse run.  Thus, the weekly
    ! mean posterior fluxes can be constructed from those outputs.  This will
    ! not work if the parameter model changes such that we are optimizing
    ! parameters which affect the flux in a nonlinear fashion.  Second,
    ! the forward run is generally much, much cheaper to run and should be
    ! able to generate "real" posterior fluxes quite quickly.


    subroutine enkf_set_inverse
      enkf_mode=enkf_mode_inverse
      flux_means_are_priors=.TRUE.
      ! print *,'EnKF mode set to inverse'
    end subroutine enkf_set_inverse

    subroutine enkf_set_forward
      enkf_mode=enkf_mode_forward
      ! print *,'EnKF mode set to forward'
      flux_means_are_priors=.TRUE.
    end subroutine enkf_set_forward

  subroutine declare_emission_averagers
!    print *,"DEBUG: in declare_emission_averagers"

    allocate(co2_bio_m(nlon360,nlat180))
    allocate(co2_fires_m(nlon360,nlat180))
    allocate(co2_ocean_m(nlon360,nlat180))
    allocate(co2_ff_m(nlon360,nlat180))
    allocate(co2_all_m(nlon360,nlat180))
    allocate(co2_gpp_m(nlon360,nlat180))
    allocate(c13_all_m(nlon360,nlat180))
    allocate(c13_bio_m(nlon360,nlat180))
    allocate(c13_gpp_m(nlon360,nlat180))
    allocate(c13_res_m(nlon360,nlat180))
    allocate(c13_ocean_m(nlon360,nlat180))
    allocate(c13_ff_m(nlon360,nlat180))
    allocate(c13_fires_m(nlon360,nlat180))
    allocate(c13_biodis_m(nlon360,nlat180))
    allocate(c13_oceandis_m(nlon360,nlat180))
    allocate(field2d_eps_m(nlon360,nlat180))

    call reset_emission_averagers

  end subroutine declare_emission_averagers

  subroutine free_emission_averagers

    if(allocated(co2_all_m)) then
       deallocate(co2_all_m)
    endif
    

    if(allocated(co2_ocean_m)) then
       deallocate(co2_ocean_m)
    endif

    if(allocated(co2_fires_m)) then
       deallocate(co2_fires_m)
    endif

    if(allocated(co2_bio_m)) then
       deallocate(co2_bio_m)
    endif

    if(allocated(field2d_eps_m)) then
       deallocate(field2d_eps_m)
    endif

    if(allocated(co2_ff_m)) then
       deallocate(co2_ff_m)
    endif

    if(allocated(co2_gpp_m)) then
       deallocate(co2_gpp_m)
    endif

    if(allocated(c13_all_m)) then
       deallocate(c13_all_m)
    endif

    if(allocated(c13_bio_m)) then
       deallocate(c13_bio_m)
    endif

    if(allocated(c13_gpp_m)) then
       deallocate(c13_gpp_m)
    endif

    if(allocated(c13_res_m)) then
       deallocate(c13_res_m)
    endif

    if(allocated(c13_ocean_m)) then
       deallocate(c13_ocean_m)
    endif

    if(allocated(c13_ff_m)) then
       deallocate(c13_ff_m)
    endif

    if(allocated(c13_fires_m)) then
       deallocate(c13_fires_m)
    endif
   
    if(allocated(c13_biodis_m)) then
       deallocate(c13_biodis_m)
    endif

    if(allocated(c13_oceandis_m)) then
       deallocate(c13_oceandis_m)
    endif



 
  end subroutine free_emission_averagers

  subroutine calc_emission_averagers

    if(allocated(co2_all_m)) then       
       if(co2_all_n .GT. 0) then
         co2_all_m=co2_all_m/co2_all_n
      endif
    endif

    if(allocated(co2_ff_m)) then       
       if(co2_ff_n .GT. 0) then
         co2_ff_m=co2_ff_m/co2_ff_n
      endif
    endif

    if(allocated(co2_bio_m)) then       
       if(co2_bio_n .GT. 0) then
         co2_bio_m=co2_bio_m/co2_bio_n
      endif
    endif

    if(allocated(field2d_eps_m)) then
       if(field2d_eps_n .GT. 0) then
         field2d_eps_m=field2d_eps_m/field2d_eps_n
      endif
    endif



    if(allocated(co2_fires_m)) then       
       if(co2_fires_n .GT. 0) then
         co2_fires_m=co2_fires_m/co2_fires_n
      endif
    endif

    if(allocated(co2_ocean_m)) then       
       if(co2_ocean_n .GT. 0) then
         co2_ocean_m=co2_ocean_m/co2_ocean_n
      endif
    endif
   
    if(allocated(co2_gpp_m)) then
       if(co2_gpp_n .GT. 0) then
         co2_gpp_m=co2_gpp_m/co2_gpp_n
      endif
    endif


 
    if(allocated(c13_all_m)) then       
       if(c13_all_n .GT. 0) then
         c13_all_m=c13_all_m/c13_all_n
      endif
    endif
    
    if(allocated(c13_bio_m)) then       
       if(c13_bio_n .GT. 0) then
         c13_bio_m=c13_bio_m/c13_bio_n
      endif
    endif        

    if(allocated(c13_gpp_m)) then
       if(c13_gpp_n .GT. 0) then
         c13_gpp_m=c13_gpp_m/c13_gpp_n
      endif
    endif

    if(allocated(c13_res_m)) then
       if(c13_res_n .GT. 0) then
         c13_res_m=c13_res_m/c13_res_n
      endif
    endif


    if(allocated(c13_ocean_m)) then
       if(c13_ocean_n .GT. 0) then
         c13_ocean_m=c13_ocean_m/c13_ocean_n
      endif
    endif

    if(allocated(c13_ff_m)) then
       if(c13_ff_n .GT. 0) then
         c13_ff_m=c13_ff_m/c13_ff_n
      endif
    endif


    if(allocated(c13_fires_m)) then
       if(c13_fires_n .GT. 0) then
         c13_fires_m=c13_fires_m/c13_fires_n
      endif
    endif

    if(allocated(c13_biodis_m)) then
       if(c13_biodis_n .GT. 0) then
         c13_biodis_m=c13_biodis_m/c13_biodis_n
      endif
    endif

    if(allocated(c13_oceandis_m)) then
       if(c13_oceandis_n .GT. 0) then
         c13_oceandis_m=c13_oceandis_m/c13_oceandis_n
      endif
    endif


  end subroutine calc_emission_averagers

  subroutine save_emission_averagers(status)

    use ParTools,  only: myid, root
    use MDF      , only : MDF_Open
    use MDF      , only : MDF_Create, MDF_Close, MDF_EndDef
    use MDF      , only : MDF_NETCDF, MDF_REPLACE, MDF_GLOBAL, MDF_CHAR, MDF_INT, MDF_FLOAT, MDF_UNLIMITED, MDF_READ
    use MDF      , only : MDF_Put_Att
    use MDF      , only : MDF_Def_Dim
    use MDF      , only : MDF_Inq_Dimid, MDF_Inquire_Dimension
    use MDF      , only : MDF_Inq_VarID
    use MDF      , only : MDF_Def_Var, MDF_Put_Var, MDF_Get_Var
    use global_data, only : outdir
    use dims,        only : itaui, itaue
    use datetime,    only : tau2date

    integer,intent(inout)               :: status

    character(len=200)    :: dsname,attname,comp,pp

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

    ! local
    integer :: i
    integer :: hid
    integer :: dim_lon_id, dim_lat_id
    integer :: var_flux_c13_all, var_flux_c13_bio, var_flux_c13_gpp, var_flux_co2_gpp, var_flux_c13_res, var_flux_c13_ocean, var_flux_c13_ff, var_flux_c13_fires, var_flux_c13_biodis, var_flux_c13_oceandis, var_flux_all,var_flux_ff, var_flux_bio, var_eps,var_flux_ocean, var_flux_fires, var_lon, var_lat, var_multiplier,var_multiplier_epsbio
    integer, dimension(6) :: idatei,idatee,idatef
    character(len=1024)   :: attstring
    
    if(myid .ne. root) return
    pp='post'
    if(flux_means_are_priors .eqv. .True.) then
       pp='prior'
    endif

    call tau2date(itaui,idatei)
    call tau2date(itaue,idatee)

    write (FluxOutFile,'(a,"/flux1x1_",i4.4,3i2.2,"_",i4.4,3i2.2,".nc")') &
                trim(outdir), idatei(1:4), idatee(1:4)

    call MDF_Create( trim(FluxOutFile), MDF_NETCDF, MDF_REPLACE, hid, status )
    IF_NOTOK_RETURN(status=1)

    ! Global Attributes 

    write(attstring,'(i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2,":",i2.2, " UTC")') idatei
    call MDF_Put_Att(hid, MDF_GLOBAL,"model_start_date",values=trim(attstring),status=status)
    IF_NOTOK_RETURN(status=1)

    write(attstring,'(i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2,":",i2.2, " UTC")') idatee
    call MDF_Put_Att(hid, MDF_GLOBAL,"model_end_date",values=trim(attstring),status=status)
    IF_NOTOK_RETURN(status=1)

    ! define dimensions:                                                                                                                                                                                            
    call MDF_Def_Dim( hid, 'lon', 360 , dim_lon_id, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Def_Dim( hid, 'lat', 180 , dim_lat_id, status )
    IF_NOTOK_RETURN(status=1)

    ! dimension variables:                                                                                                                                                                                                    
    call MDF_Def_Var( hid, 'lon', MDF_FLOAT, (/dim_lon_id/), var_lon, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid, var_lon,"units",values="degrees East",status=status)
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid, var_lon,"comment",values="center gridbox",status=status)
    IF_NOTOK_RETURN(status=1)

    call MDF_Def_Var( hid, 'lat', MDF_FLOAT, (/dim_lat_id/), var_lat, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid, var_lat,"units",values="degrees North",status=status)
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid, var_lat,"comment",values="center gridbox",status=status)
    IF_NOTOK_RETURN(status=1)

    ! define variables:                                                                                                                                                                                                    
    dsname='flux_multiplier_m'

    call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_multiplier, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid, var_multiplier,"units",values="dimensionless",status=status)
    IF_NOTOK_RETURN(status=1)

    dsname='flux_multiplier_epsbio_m'
    call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_multiplier_epsbio, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att(hid, var_multiplier_epsbio,"units",values="dimensionless",status=status)
    IF_NOTOK_RETURN(status=1)



    if(allocated(co2_all_m)) then

       comp='all'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_all, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_all,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_all,attname,values=co2_all_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif
    if(allocated(co2_ff_m)) then

       comp='ff'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_ff, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_ff,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_ff,attname,values=co2_ff_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif
    if(allocated(co2_bio_m)) then

       comp='bio'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_bio, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_bio,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_bio,attname,values=co2_bio_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif
    if(allocated(field2d_eps_m)) then

       comp='eps'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_eps, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_eps,"units",values="permil",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_eps,attname,values=field2d_eps_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif




    if(allocated(co2_fires_m)) then

       comp='fires'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_fires, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_fires,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_fires,attname,values=co2_fires_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif
    if(allocated(co2_ocean_m)) then

       comp='ocean'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_ocean, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_ocean,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_ocean,attname,values=co2_ocean_n,status=status)
       IF_NOTOK_RETURN(status=1)
 
    endif
    if(allocated(c13_all_m)) then

       comp='c13_all'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_all, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_all,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_all,attname,values=c13_all_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif
    
           

     if(allocated(c13_bio_m)) then

       comp='c13_bio'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_bio, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_bio,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_bio,attname,values=c13_bio_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif

     if(allocated(c13_gpp_m)) then

       comp='c13_gpp'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_gpp, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_gpp,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_gpp,attname,values=c13_gpp_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif

     if(allocated(co2_gpp_m)) then

       comp='gpp'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_co2_gpp, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_co2_gpp,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_co2_gpp,attname,values=co2_gpp_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif


     if(allocated(c13_res_m)) then

       comp='c13_res'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_res, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_res,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_res,attname,values=c13_res_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif



     if(allocated(c13_ocean_m)) then

       comp='c13_ocean'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_ocean, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_ocean,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_ocean,attname,values=c13_ocean_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif

     if(allocated(c13_ff_m)) then

       comp='c13_ff'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_ff, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_ff,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_ff,attname,values=c13_ff_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif

     if(allocated(c13_fires_m)) then

       comp='c13_fires'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_fires, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_fires,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_fires,attname,values=c13_fires_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif

     if(allocated(c13_biodis_m)) then

       comp='c13_biodis'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_biodis, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_biodis,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_biodis,attname,values=c13_biodis_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif

     if(allocated(c13_oceandis_m)) then

       comp='c13_oceandis'
       dsname='flux_'//trim(comp)//'_'//trim(pp)//'_mean'
       attname='flux_'//trim(comp)//'_'//trim(pp)//'_n'

       call MDF_Def_Var( hid, trim(dsname), MDF_FLOAT, (/dim_lon_id,dim_lat_id/), var_flux_c13_oceandis, status )
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_oceandis,"units",values="mol m-2 s-1",status=status)
       IF_NOTOK_RETURN(status=1)
       call MDF_Put_Att(hid, var_flux_c13_oceandis,attname,values=c13_oceandis_n,status=status)
       IF_NOTOK_RETURN(status=1)

    endif





    ! finished definition:                                                                                                                                                                                          
    call MDF_EndDef( hid, status )
    IF_NOTOK_RETURN(status=1)

    call MDF_Put_Var( hid, var_lon, (/(i+0.5,i=-180,179)/), status)
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Var( hid, var_lat, (/(i+0.5,i=-90,89)/), status)
    IF_NOTOK_RETURN(status=1)

    call MDF_Put_Var( hid, var_multiplier, fluxmultiplier(:,:,1), status)
    IF_NOTOK_RETURN(status=1)

    call MDF_Put_Var( hid, var_multiplier_epsbio, fluxmultiplier_epsbio(:,:,1), status)
    IF_NOTOK_RETURN(status=1)


    if(allocated(co2_all_m)) then
       call MDF_Put_Var( hid, var_flux_all, co2_all_m, status)
       IF_NOTOK_RETURN(status=1)
    endif
    if(allocated(co2_ff_m)) then
       call MDF_Put_Var( hid, var_flux_ff, co2_ff_m, status)
       IF_NOTOK_RETURN(status=1)
    endif
    if(allocated(co2_bio_m)) then
       call MDF_Put_Var( hid, var_flux_bio, co2_bio_m, status)
       IF_NOTOK_RETURN(status=1)
    endif
    if(allocated(field2d_eps_m)) then
       call MDF_Put_Var( hid, var_eps, field2d_eps_m, status)
       IF_NOTOK_RETURN(status=1)
    endif
    if(allocated(co2_ocean_m)) then
       call MDF_Put_Var( hid, var_flux_ocean, co2_ocean_m, status)
       IF_NOTOK_RETURN(status=1)
    endif
    if(allocated(co2_fires_m)) then
       call MDF_Put_Var( hid, var_flux_fires, co2_fires_m, status)
       IF_NOTOK_RETURN(status=1)
    endif
    
    if(allocated(c13_all_m)) then
       call MDF_Put_Var( hid, var_flux_c13_all, c13_all_m, status)
       IF_NOTOK_RETURN(status=1)
    endif
    
    if(allocated(c13_bio_m)) then
       call MDF_Put_Var( hid, var_flux_c13_bio, c13_bio_m, status)
       IF_NOTOK_RETURN(status=1)
    endif        

    if(allocated(c13_gpp_m)) then
       call MDF_Put_Var( hid, var_flux_c13_gpp, c13_gpp_m, status)
       IF_NOTOK_RETURN(status=1)
    endif

    if(allocated(co2_gpp_m)) then
       call MDF_Put_Var( hid, var_flux_co2_gpp, co2_gpp_m, status)
       IF_NOTOK_RETURN(status=1)
    endif


    if(allocated(c13_res_m)) then
       call MDF_Put_Var( hid, var_flux_c13_res, c13_res_m, status)
       IF_NOTOK_RETURN(status=1)
    endif


    if(allocated(c13_ocean_m)) then
       call MDF_Put_Var( hid, var_flux_c13_ocean, c13_ocean_m, status)
       IF_NOTOK_RETURN(status=1)
       print*,'PUT VAR',sum(c13_ocean_m)
    endif

    if(allocated(c13_ff_m)) then
       call MDF_Put_Var( hid, var_flux_c13_ff, c13_ff_m, status)
       IF_NOTOK_RETURN(status=1)
       print*,'PUT VAR',sum(c13_ff_m)
    endif

    if(allocated(c13_fires_m)) then
       call MDF_Put_Var( hid, var_flux_c13_fires, c13_fires_m, status)
       IF_NOTOK_RETURN(status=1)
       print*,'PUT VAR',sum(c13_fires_m)
    endif

    if(allocated(c13_biodis_m)) then
       call MDF_Put_Var( hid, var_flux_c13_biodis, c13_biodis_m, status)
       IF_NOTOK_RETURN(status=1)
       print*,'PUT VAR',sum(c13_biodis_m)
    endif

    if(allocated(c13_oceandis_m)) then
       call MDF_Put_Var( hid, var_flux_c13_oceandis, c13_oceandis_m, status)
       IF_NOTOK_RETURN(status=1)
       print*,'PUT VAR',sum(c13_oceandis_m)
    endif




    ! close file:                                                                                                                                                                                                   
    call MDF_Close( hid, status )
    IF_NOTOK_RETURN(status=1)

    write (gol,'("[save_emission_averagers]  Done writing arrays and closing output file.")'); call goPr

  end subroutine save_emission_averagers

  subroutine reset_emission_averagers

    if(allocated(co2_all_m)) then
       co2_all_m=0.
       co2_all_n=0
    endif

    if(allocated(co2_ocean_m)) then
       co2_ocean_m=0.
       co2_ocean_n=0
    endif

    if(allocated(co2_fires_m)) then
       co2_fires_m=0.
       co2_fires_n=0
    endif

    if(allocated(co2_bio_m)) then
       co2_bio_m=0.
       co2_bio_n=0
    endif

    if(allocated(field2d_eps_m)) then
       field2d_eps_m=0.
       field2d_eps_n=0
    endif


    if(allocated(co2_ff_m)) then
       co2_ff_m=0.
       co2_ff_n=0
    endif
    
    
    if(allocated(c13_all_m)) then
       c13_all_m=0.
       c13_all_n=0
    endif
    
    if(allocated(c13_bio_m)) then
       c13_bio_m=0.
       c13_bio_n=0
    endif        

    if(allocated(c13_gpp_m)) then
       c13_gpp_m=0.
       c13_gpp_n=0
    endif

    if(allocated(co2_gpp_m)) then
       co2_gpp_m=0.
       co2_gpp_n=0
    endif

    if(allocated(c13_res_m)) then
       c13_res_m=0.
       c13_res_n=0
    endif


    if(allocated(c13_ocean_m)) then
       c13_ocean_m=0.
       c13_ocean_n=0
    endif

    if(allocated(c13_ff_m)) then
       c13_ff_m=0.
       c13_ff_n=0
    endif

    if(allocated(c13_fires_m)) then
       c13_fires_m=0.
       c13_fires_n=0
    endif

    if(allocated(c13_biodis_m)) then
       c13_biodis_m=0.
       c13_biodis_n=0
    endif

    if(allocated(c13_oceandis_m)) then
       c13_oceandis_m=0.
       c13_oceandis_n=0
    endif


  end subroutine reset_emission_averagers

    subroutine ReadFromFile3D( filename, sds_name , data_array, status, netcdf4)
    use GO        , only : gol, goErr, goPr, goLabel, goUpCase
    !use file_hdf
    use MDF

    !I/O
    character(len=*),intent(in)  :: filename  
    character(len=*),intent(in)  :: sds_name
    real,dimension(:,:,:),intent(out):: data_array
    logical,optional:: netcdf4
    integer, intent (inout)      :: status
   
    character(len=64) :: varname
    integer :: hid
    integer             ::  n, k, l
    character(len=*), parameter ::  rname = mname//'/ReadFromFile'

    if (present(netcdf4)) then 
        call MDF_Open(filename=trim(filename), ftype=MDF_NETCDF4, mode=MDF_READ, hid=hid, status=status)
        IF_NOTOK_RETURN(status=1)
    else
        call MDF_Open(filename=trim(filename), ftype=MDF_HDF4, mode=MDF_READ, hid=hid, status=status)
        IF_NOTOK_RETURN(status=1)
    endif


    call MDF_Inquire(hid, status, nVariables = n )
    IF_NOTOK_RETURN(status=1)

    k = 0
    do
      k = k + 1
      if ( k > n ) then 
          write(gol,*) 'SDS not found' ; call goPr
          write(gol,*) 'filename   : ', trim(filename) ; call goPr
          IF_NOTOK_RETURN(status=1)
      endif
      call MDF_Inquire_Variable( hid, k, status, name=varname )
      IF_NOTOK_RETURN(status=1)

      if (trim(sds_name) .ne. trim(varname) ) cycle  ! name not correct

      ! all lights green, go read
      call MDF_Get_Var( hid, k , data_array, status )
      IF_ERROR_RETURN(status=1)
      exit ! break from loop

    end do
    call MDF_Close( hid, status )
    IF_ERROR_RETURN(status=1)
    status=0

    end subroutine ReadFromFile3D
    subroutine ReadFromFile( filename, sds_name , data_array, status, sds_date, netcdf4)
    use GO        , only : gol, goErr, goPr, goLabel, goUpCase
    !use file_hdf
    use MDF

    !I/O
    character(len=*),intent(in)  :: filename  
    character(len=*),intent(in)  :: sds_name
    real,dimension(:,:),intent(out):: data_array
    integer,dimension(6),optional:: sds_date
    logical,optional:: netcdf4
    integer, intent (inout)      :: status
   
    character(len=64) :: varname
    integer :: hid
    integer             ::  n, k, l
    character(len=*), parameter ::  rname = mname//'/ReadFromFile'
    integer,dimension(6) :: idate_f

    if (present(netcdf4)) then 
        call MDF_Open(filename=trim(filename), ftype=MDF_NETCDF4, mode=MDF_READ, hid=hid, status=status)
        IF_NOTOK_RETURN(status=1)
    else
        call MDF_Open(filename=trim(filename), ftype=MDF_HDF4, mode=MDF_READ, hid=hid, status=status)
        IF_NOTOK_RETURN(status=1)
    endif


    call MDF_Inquire(hid, status, nVariables = n )
    IF_NOTOK_RETURN(status=1)

    k = 0
    do
      k = k + 1
      if ( k > n ) then 
          write(gol,*) 'SDS not found' ; call goPr
          write(gol,*) 'filename   : ', trim(filename) ; call goPr
          write(gol,*) 'sds_name   : ', trim(sds_name) ; call goPr
          if (present(sds_date))   write(gol,*) 'sds_date   :',sds_date ; call goPr
          IF_NOTOK_RETURN(status=1)
      endif
      call MDF_Inquire_Variable( hid, k, status, name=varname )
      IF_NOTOK_RETURN(status=1)

      if (trim(sds_name) .ne. trim(varname) ) cycle  ! name not correct

      if (present(sds_date)) then
           call MDF_Get_Att( hid, k, 'idate', idate_f, status)
           IF_NOTOK_RETURN(status=1)
           do l = 1,6
               if (sds_date(l) .ne. idate_f(l) ) status = 1
           enddo
           if (status == 1) cycle
      endif


      ! all lights green, go read
      call MDF_Get_Var( hid, k , data_array, status )
      IF_ERROR_RETURN(status=1)
      exit ! break from loop

    end do
    call MDF_Close( hid, status )
    IF_ERROR_RETURN(status=1)
    status=0

    end subroutine ReadFromFile

    subroutine GetParameterMaps(status)

    use GO     , only : goLoCase,gol, goErr, goPr, goBug, goNum2str
    use MDF    , only : MDF_Open
    use MDF    , only : MDF_Create, MDF_Close, MDF_EndDef
    use MDF    , only : MDF_NETCDF4, MDF_REPLACE, MDF_GLOBAL, MDF_CHAR, MDF_INT, MDF_FLOAT, MDF_UNLIMITED, MDF_READ
    use MDF    , only : MDF_Put_Att
    use MDF    , only : MDF_Def_Dim
    use MDF    , only : MDF_Inq_Dimid, MDF_Inquire_Dimension
    use MDF    , only : MDF_Inq_VarID
    use MDF    , only : MDF_Def_Var, MDF_Put_Var, MDF_Get_Var

    use Partools, only: ntracetloc, ntracet_ar, myid
    use chem_param, only : nmembersloc,names

    implicit none

    ! ------------- I/O ------------------
    integer, intent(out) :: status
    ! ------------- local ------------------
    character(len=*), parameter  ::  rname = mname//'/GetParameterMaps'

    character(len=1024)          :: inFile = ''
    logical                      :: file_exists

    integer :: n, offsetn, nloop
    integer :: hid,dimid, varid


    !WP! This routine opens a file with 1x1 degree parameters to be used in
    !WP! CarbonTracker. For now, these are assumed to be scaling factors that are
    !WP! applied to the ocean and biosphere at each time step. We read one file 
    !WP! for every ensemble member (inv) or one just for the mean multipliers
    !WP! (fwd).
    !WP!
    !WP! The file names are constructed from a root file name that comes from
    !WP! the rc-file and is read into the ParamInputFile variable during Emis_Init. 
    !WP! We add a qualifier with the tracer number to the filename to read one
    !WP! file per member for the inv runs.

    if (enkf_is_inverse()) then
        write (gol,'("[GetParameterMaps] Parameter files will be read for each member ")') ; call goPr
        nloop=nmembersloc !ntracetloc
        !nloop=ntracetloc 
    else
        write (gol,'("[GetParameterMaps] Parameter file will be read for mean ")') ; call goPr
        nloop=1
    endif


    do n=1,nloop  ! loop over tracers (members) on this processor
        !IV! This version of carbontracker uses two kinds of tracers: CO2 and 13CO2. 
        !IV! It is designed to use an even number of processors. The first half of it are used 
        !IV! for CO2 ensembles, 10 ensembles per processors, and the other half are used for 13CO2 ensembles. 
        !IV! That requires that same parameters files must be used for both kinds of tracers. As such, for first 13CO2
        !IV! ensemble the offsetn is set back to zero. 

        offsetn = sum(ntracet_ar(0:myid-1))
        if (goLoCase(trim(names(offsetn+1)))=='co2c13'.or.goLoCase(trim(names(offsetn+1)))=='co2c13_bg') then
            offsetn=offsetn-sum(ntracet_ar)/2
        endif
        if (trim(ParamInputFile)=='None') then
            write (gol,'("[GetParameterMaps] ..... File specified is: ",a,".")') trim(ParamInputFile); call goPr
            write (gol,'("[GetParameterMaps] ..... Using parametervalues of 1.0 as default")'); call goPr
            fluxmultiplier = 1.0
            fluxmultiplier_biodis = 1.0
            fluxmultiplier_epsbio = 1.0
            status=0
            return
        endif
        
        inFile = trim(ParamInputFile)//'.'//trim(goNum2str(offsetn+n-1,fmt='(i3.3)'))//'.nc'
        INQUIRE(FILE=trim(inFile), EXIST=file_exists) 
        if (.not.file_exists) then
            write (gol,'("[GetParameterMaps] ..... File to open does not exist ",a,".")') trim(inFile); call goErr
            status=1
            IF_ERROR_RETURN(status=1)
        endif  

        ! open existing file:
        call MDF_Open( trim(inFile), MDF_NETCDF4, MDF_READ, hid=hid, status=status )
        IF_NOTOK_RETURN(status=1)

        call MDF_Inq_VarID(hid, "parametermap", varid, status)
        if (status /= 0) then
            write (gol,'("[GetParameterMaps] ..... required variable parametermap does not exist")') ; call goErr
            status=1
            IF_NOTOK_RETURN(status=1)
        endif

        call MDF_Get_Var(hid,varid,fluxmultiplier(:,:,n),status)
        if (status /= 0) then
            write (gol,'("[GetParameterMaps] ..... failed to read variable parametermap")') ; call goErr
            status=1
            IF_NOTOK_RETURN(status=1)
        endif






!        !IV 2nd set parameter files

        call MDF_Inq_VarID(hid, "parametermap_epsbio", varid, status)
        if (status /= 0) then
            write (gol,'("[GetParameterMaps] ..... required variable parametermap does not exist")') ; call goErr
            status=1
            IF_NOTOK_RETURN(status=1)
        endif

        call MDF_Get_Var(hid,varid,fluxmultiplier_epsbio(:,:,n),status)
         if (status /= 0) then
            write (gol,'("[GetParameterMaps] ..... failed to read variable parametermap")') ; call goErr
            status=1
            IF_NOTOK_RETURN(status=1)
        endif

        call MDF_Close(hid,status)
        IF_NOTOK_RETURN(status=1)



    end do

    write (gol,'("[GetParameterMaps] ..... All parametermaps read successfully")') ; call goPr

    !ok
    status = 0

    end subroutine GetParameterMaps

    subroutine get_ntimesteps(ntimesteps)
      use dims,        only: idate
      implicit none
      integer, intent(out) :: ntimesteps
      if(idate(2) .eq. 1 .or. idate(2) .eq. 3 .or. idate(2) .eq. 5 .or. idate(2) .eq. 7 .or. idate(2) .eq. 8 .or. idate(2) .eq. 10 .or. idate(2) .eq. 12) then
        ntimesteps = 248
      else if (idate(2) .eq. 4 .or. idate(2) .eq. 6 .or. idate(2) .eq. 9 .or. idate(2) .eq. 11) then
        ntimesteps = 240
      else if (idate(2) .eq. 2) then
        ntimesteps = 224
      endif
    end subroutine get_ntimesteps


end module emission_common
