ecearth3/sources/tm5mp/proj/co2/user_output_c4mip.F90

#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') 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
#define IF_NOTOK_MDF(action) if (status/=0) then; TRACEBACK; action; call MDF_CLose(fid,status); status=1; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      USER_OUTPUT_C4MIP
!
! !DESCRIPTION: 
!   
!    This module provides the code needed to produce the CMIP6 C4MIP
!    output from TM5. Code is based on the user_output_aerocom.
!
!    output_c4mip_init:
!     - initialise the list of variables (allvars)
!     - initialise the data holder within each variable (2Ddata/3Ddata,...)
!     - initialise the output netcdf files, one for eacht variable
!    output_c4mip_accumulate:
!     - do accumulation for all variables and save data to either
!       2Ddata or 3Ddata holder of the variable (excluding optical vars)
!    output_c4mip_write
!     - write the monthly variable data to netcdf-file
!     - initialise the dataholders for accumulation of new month
!    output_c4mip_write_hourly
!     - write the hourly variable data to netcdf-file
!     - initialise the dataholders for accumulation of new hour
!    output_c4mip_write_daily
!     - write the daily variable data to netcdf-file
!     - initialise the dataholders for accumulation of new day
!    write_var
!    
! Tommi Bergman 1.11.2019
!\\
!\\
! !INTERFACE: 
!
MODULE USER_OUTPUT_C4MIP
  !
  ! !USES:
  !
  use go,           only : gol, goErr, goPr, goLabel
  use GO,             ONLY : GO_Timer_Def, GO_Timer_End, GO_Timer_Start
  use dims,         only : nregions !=1, support for zooming with larger values not available for C4MIP
  use meteodata,    only : global_lli, levi
  use MDF
  use TM5_DISTGRID, only : dgrid, Get_DistGrid, update_halo, update_halo_iband
  use chem_param
!#ifdef with_m7
!  use m7_data,     only : h2o_mode
!#endif
  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: output_c4mip_init
  !public :: output_c4mip_step
  public :: output_c4mip_write_monthly
  public :: output_c4mip_write_hourly
  public :: output_c4mip_write_6hourly
  public :: output_c4mip_write_daily
  public :: output_c4mip_done
  public :: accumulate_c4mip_data
!  public :: wdep_out
  character(len=*), parameter ::  mname = 'user_output_c4mip'
  ! max indices, at maximum 7, one for each mode
  integer,parameter :: n_indices=11

  type varfile
     integer                         :: itm5             !
     character(len=16)               :: vname            !
     character(len=64)               :: lname            !
     character(len=11)               :: unit             !
     character(len=10)               :: positive         !
     character(len=130)              :: standard_name    !
     real,dimension (:,:),pointer    :: data2D           !
     real,dimension (:,:,:),pointer  :: data3D           !
!     real,dimension (:,:,:),pointer  :: budgetdata          !
     
     integer                         :: varid            !
     integer                         :: fileunit         ! file unit number 
     character(len=200)              :: filename         !
     integer                         :: dimensions       ! 
     integer                         :: lon_varid       ! 
     integer                         :: lat_varid       ! 
     integer                         :: lev_varid       ! 
     integer                         :: time_varid
     integer                         :: hyam_varid
     integer                         :: hybm_varid
     integer                         :: hyai_varid
     integer                         :: hybi_varid
     integer                         :: bounds_varid
     integer                         :: dims
     character(len=10)               :: freq
     character(len=9)                :: class            ! which class of variable :emi, ddep, wdep,conc,aod,met,crescendo
     integer,dimension(n_indices)         :: indices
     real                            :: xmgas
     character(len=20)               :: table_id
  end type varfile

  type dimdata
     integer               :: nlon             ! size of x dimension of requested field
     integer               :: nlat             ! size of y dimension of requested field
     integer               :: nlev             ! size of z dimension of requested field
     real,dimension(:),pointer                :: lon             ! x dimension of requested field
     real,dimension(:),pointer                :: lat             ! y dimension of requested field
     real,dimension(:),pointer               :: lev             ! z dimension of requested field
     integer               :: lonid             ! x dimension id in nc
     integer               :: latid             ! y dimension id in nc
     integer               :: levid             ! z dimension id in nc
     integer               :: timeid             ! time dimension id in nc
     integer               :: time_varid

  end type dimdata
  type(dimdata)::dimension_data


  !!!!
  integer::test_fileunit
  !!!!

  integer :: n_vars=0
  type(varfile), dimension(:), allocatable :: allvars
  type(varfile), dimension(:), allocatable :: fixedvars
  integer :: n_var_max=300
  integer :: n_fixed=3
  integer, public :: n_days_in_month
  character(len=20), public :: c4mip_exper ! AeroCom experiment name
  
  integer, save             :: od550aer, &
                               areacella,&
                               sftlf,&
                               orog

  integer                                   ::  fid ! file id for IF_NOTOK_MDF macro
  integer                                   :: access_mode
  integer                                   :: accumulation_mon,accumulation_day,accumulation_hr,accumulation_6hr
  integer                                   :: timeidx_mon,timeidx_hr,timeidx_day,timeidx_6hr

  integer,parameter::iemiunit=1
  integer,parameter::iddepunit=1 !same dimensions as emi
  integer,parameter::iwdepunit=1 !same dimensions as emi
  integer,parameter::iprod3dunit=2
  integer,parameter::immrunit=3
  integer,parameter::idimensionlessunit=4
  integer,parameter::iheightunit=5
  integer,parameter::itempunit=6
  integer,parameter::io3unit=7
  integer,parameter::ipresunit=8
  integer,parameter::ivmrunit=9
  integer,parameter::irateunit=10
  integer,parameter::iloadunit=11
  integer,parameter::iextunit=12
  integer,parameter::iccunit=13
  integer,parameter::im3unit=14
  integer,parameter::imassunit=15


  character(len=11),dimension(15),parameter::units=(/'kg m-2 s-1','kg m-3 s-1','kg kg-1','1','m','K','DU','Pa','mole mole-1',&
  's-1','kg m-2','m-1','cm-3','m-3','kg'/)  
  character (len=11), parameter::unit='m-3'


  Character(len=8),dimension(3),parameter   :: monthly_var=(/'ps','co2','co2mass'/)
  character(len=11),dimension(3),parameter:: monthly_varunit=(/units(ipresunit), units(ivmrunit), units(iloadunit)/)
  real,dimension(3),parameter   :: xmmonthly_var=(/1.0,xmco2,xmco2/)
  integer,dimension(3),parameter::monthly_dim=(/2,3,2/)


  !SPECIAL
  !6HOURLY
  !character(len=8),dimension(1),parameter:: ps6hr=(/'ps'/)
  !character(len=11),dimension(1),parameter:: ps6hrunit=(/units(ipresunit)/)

  !HOURLY
  character(len=8),dimension(3),parameter:: hourly_var=(/'ps','co2','co2mass'/) !,'co2mass1'/)
  character(len=11),dimension(3),parameter:: hourly_varunit=(/units(ipresunit), units(ivmrunit), units(iloadunit)/) ! , 'kg(co2)'/)
  real,dimension(3),parameter ::xmhourly=(/1.0,xmco2,xmco2/) !,xmco2/)
  integer,dimension(3),parameter::hourly_dim=(/2,3,2/) !,0/)

  !DAILY
  character(len=8),dimension(3),parameter:: daily_var=(/'ps','co2','co2mass'/)
  character(len=11),dimension(3),parameter:: daily_varunit=(/ units(ipresunit),units(ivmrunit), units(iloadunit)/)
  real,dimension(3),parameter ::xmdaily=(/-1.0,xmco2,xmco2/)
  integer,dimension(3),parameter::daily_dim=(/2,3,2/)


  ! global attributes that might change with run or something else
  character(len=3),parameter::grid='3x2'!'250 km'
  character(len=3),parameter::grid_label='gn'!'gnz' for zonal means
  character(len=300),parameter::c4mip_source='EC-Earth3-CC (2017): atmosphere: IFS cy36r4 (TL255, linearly &
       &reduced Gaussian grid equivalent to 512 x 256, 91 levels, top level: 0.01 hPa);atmospheric_chemistry: &
       &TM5 (3 deg. (long.) x 2 deg. (lat.), 34 levels, top level: 0.1 hPa; aerosol: TM5'
  character(len=17),parameter::c4mip_source_id='EC-Earth3-CC'
  character(len=20),public ::c4mip_source_type!='AOGCM CHEM AER' !or 'AGCM CHEM AER' for amip-runs
  character(len=60),public ::c4mip_realm
  character(len=60),public::c4mip_experiment_id
  character(len=60),public::c4mip_experiment
  character(len=1),public::realization_i='1'
  character(len=1),public::physics_i='1'
  character(len=1),public::forcing_i='1'
  character(len=1),public::initialization_i='1'
  integer, public:: c4mip_dhour
 ! Timers
  integer :: itim_init, itim_addvar, itim_write, itim_accu, itim_attr, itim_accu_opt, itim_write_hour, itim_write_day, &
       itim_write_mon, itim_write_gather
contains

  subroutine output_c4mip_init(status)
    ! Open files
    ! allocate dataholders
    use dims,            only : newsrun,itau,mlen
    use global_data,     only : outdir
    use datetime,        only : tau2date, date2tau
    use partools,        only : MPI_INFO_NULL, localComm
#ifdef with_m7
    !use optics,          only : Optics_Init
    !use sedimentation,   only : ised,nsed
#endif
    use partools    ,    only : isRoot,myid
    use global_data,     only : region_dat
    use tm5_distgrid,    only : gather
    use meteodata      , only : lsmask_dat,oro_dat
    use Binas      ,     only : grav
    
    implicit none

    !OUTPUT parameters
    integer, intent(out)     :: status

    !LOCAL parameters
    integer  :: region !iterator for regions
    integer  :: nlon_region
    integer  :: nlat_region
    integer  :: nlev_region ! also global

    integer  :: j_var
    !integer :: nlev_region ! also global
    !integer :: nlev_region ! also global
    integer  :: itrac
    integer  :: i_sed
    integer  :: i,i1,i2,j1,j2,k,j,imr,jmr
    character(len=*), parameter ::  rname = mname//'/output_c4mip_init'
    !FIXED  VARS
    real, dimension(:),pointer :: dxyp 
    real, allocatable :: arr2d(:,:)
    real ::xmcomp

    call goLabel(rname)
    
    ! define timers:
    call GO_Timer_Def( itim_init, 'output c4mip init', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_write, 'output c4mip write', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_write_gather, 'output c4mip write gather', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_write_day, 'output c4mip write day', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_write_hour, 'output c4mip write hour', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_write_mon, 'output c4mip write mon', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_accu, 'output c4mip accu', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_attr, 'output c4mip attr', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_addvar, 'output c4mip addvar', status )
    IF_NOTOK_RETURN(status=1)

    call Get_DistGrid( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2)
    
    if (newsrun) then
       
       accumulation_mon=0.0
       accumulation_hr=0.0
       accumulation_6hr=0.0
       accumulation_day=0.0
       region=1
       ! intermediate structures for budgets

       imr = global_lli(1)%nlon
       jmr = global_lli(1)%nlat   
       ! for areacella,orog,sftlf
       if (isRoot) then
          allocate( arr2d(imr,jmr)       )
       else
          allocate( arr2d(1,1)   )
       endif
       arr2d(:,:)=0.0
       ! for monthly output
       ! initialise with 31 for january
       n_days_in_month=31

    end if

    call GO_Timer_Start( itim_init, status )
    IF_NOTOK_RETURN(status=1)
    ! C4MIP only available for global-> region=1
    region=1
    
    !Initialise grid definitions
    nlon_region = global_lli(region)%nlon
    nlat_region = global_lli(region)%nlat
    nlev_region = levi%nlev
    dimension_data%nlon= nlon_region
    dimension_data%nlat= nlat_region
    dimension_data%nlev= nlev_region
    allocate(dimension_data%lon(nlon_region))       
    allocate(dimension_data%lat(nlat_region))       
    allocate(dimension_data%lev(nlev_region))       
    dimension_data%lon=global_lli(region)%lon_deg
    dimension_data%lat=global_lli(region)%lat_deg

    ! initialise output timeidx used for keeping track which output steps is written
    timeidx_mon=1
    timeidx_day=1
    timeidx_hr=1
    timeidx_6hr=1
    ! allocate room for variables
    allocate(allvars(n_var_max))       
    allocate(fixedvars(n_fixed))       


!!$    do i=1,size(ps6hr)
!!$        call add_variable(-1,trim(ps6hr(i)),trim(ps6hr(i)),ps6hrunit(i),2,status,'ps6h','AER6hr',-1.0)
!!$    end do

   
    ! Gas-phase species volume mixingratios
    do i=1,size(monthly_var)
       write(gol,*) 'monthly_var add,',trim(monthly_var(i))
       
       if (xmmonthly_var(i)>0.0) then
          call add_variable(-1,trim(monthly_var(i)),'volume mixing ratio of '//trim(monthly_var(i)), hourly_varunit(i),monthly_dim(i),status,'monthly','AERmon',xmmonthly_var(i))
       else
          write(gol,*) 'monthly_var with negative molar mass'
       end if
    end do

    ! add hourly output
    do i=1,size(hourly_var)
       call add_variable(-1,trim(hourly_var(i)),trim(hourly_var(i)),hourly_varunit(i),hourly_dim(i),status,'hourly','AERhr',xmhourly(i))

    end do
    ! add daily fields
    do i=1,size(daily_var)
       call add_variable(-1,trim(daily_var(i)),trim(daily_var(i)),daily_varunit(i),daily_dim(i),status,'daily','AERday',xmdaily(i))
    end do


    call add_variable(-1,'areacella','cell area','m2',2,status,'fixed','AERfx',-1.0)
    call add_variable(-1,'orog','surface_altitude','m',2,status,'fixed','AERfx',-1.0)
    call add_variable(-1,'sftlf','land_area_fraction','1',2,status,'fixed','AERfx',-1.0)
    ! 
    
   
    do j_var = 1, n_vars
       ! initialise a single file for each variables as per C4MIP specification
       
       ! overwrite existing files (clobber)
       if (isroot)call MDF_Create( allvars(j_var)%filename, MDF_NETCDF4, MDF_REPLACE, allvars(j_var)%fileunit, status )
       IF_NOTOK_RETURN(status=1)
       
       
       !For each file 
       ! write grid dimension attributes
       if (isroot) call write_dimensions(status, j_var)
       
       IF_NOTOK_RETURN(status=1)
       ! write global attributes
       if (isroot) call write_attributes(status, j_var)
       IF_NOTOK_RETURN(status=1)
       !write dimension variables
       if (isroot) call write_var(status,j_var)
       IF_NOTOK_RETURN(status=1)
       if (allvars(j_var)%table_id=='AERfx')then
          if (trim(allvars(j_var)%vname)=='areacella')then
             ! Gridbox area
             dxyp   => region_dat(1)%dxyp
             
             do j=j1,j2
                allvars(j_var)%data2D(i1:i2,j)=dxyp(j)
             end do
             call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0,status)
             if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1/), count=(/imr,jmr/) )
             IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
        
          else  if (trim(allvars(j_var)%vname)=='orog')then
             ! Gridbox area
             allvars(j_var)%data2D(i1:i2,j1:j2) =oro_dat(region)%data(i1:i2,j1:j2,1)/grav
             call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0,status)
             if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1/), count=(/imr,jmr/) )
             IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)


          else  if (trim(allvars(j_var)%vname)=='sftlf')then
             ! Gridbox area

             allvars(j_var)%data2D(i1:i2,j1:j2)=lsmask_dat(1)%data(i1:i2,j1:j2,1)/100.
             call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0,status)
             if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1/), count=(/imr,jmr/) )
             IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

          end if
       end if
    end do
    
    deallocate(arr2d)
    call GO_Timer_End( itim_init, status )
    IF_NOTOK_RETURN(status=1)       
    call goLabel()
    status = 0

  end subroutine output_c4mip_init
  
  subroutine output_c4mip_write_monthly(region,newhour,status)
    use MeteoData,      only : temper_dat,sst_dat,albedo_dat,ci_dat,lsp_dat,cp_dat,ssr_dat,str_dat,&
         blh_dat,gph_dat,lwc_dat,iwc_dat,cco_dat,cc_dat,humid_dat, m_dat,phlb_dat,sp_dat   !
    use global_data,    only : conv_dat
    use GO,             only : TDate, NewDate
    use go_date,only: days_in_month!
    use datetime,       only : tau2date,date2tau,julday   !
    use dims,           only : itau,iyear0          !current time
    !use ebischeme,      only : AC_diag_prod,AC_O3_lp,AC_loss
    use tm5_distgrid, only : dgrid,    Get_DistGrid ,gather
    use partools    , only : isRoot,myid
  
    !    use domain_decomp, only: gather
    implicit none
    logical,intent(in) ::newhour
    integer,intent(out)::status
    integer::region
    integer:: j_var
    integer:: imr,jmr,i,i1,i2,j1,j2,lmr
    character(len=*), parameter ::  rname = mname//'/output_c4mip_write_monthly'
    real, allocatable ::  arr3d(:,:,:),arr3dh(:,:,:),arr2d(:,:)
    integer, dimension(6) :: curdate    
    ! reference time:
    integer, parameter            ::  time_reftime6(6) = (/1750,01,01,00,00,00/)
    integer(kind=8)                       ::  itauref  ! reftime in seconds
    real                          ::  reftime  ! seconds from reference time
    real                          :: rangemon
    type(Tdate)::curdate_tdate
    call goLabel(rname)

    call GO_Timer_Start( itim_write_mon, status )
    IF_NOTOK_RETURN(status=1)

    if (region > 1) then
       write(gol,*) 'output_c4mip_write_monthly: region >1, only available in global mode!'
       call goErr
       status=1
       return
    end if

    call Get_DistGrid( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2)
    
    ! entire region grid size
    imr = global_lli(1)%nlon
    jmr = global_lli(1)%nlat
    lmr = levi%nlev

    ! define the reference time in seconds (itauref)
    ! (for now in days since 1850-01-01 00:00, local variable) 
    ! returns the difference to current time, beginning of new step
    ! 
    call date2tau( time_reftime6, itauref )
    
    ! calculate reftime as fractional days from itauref, hence division by 86400
    ! call date2tau( idater, itaucur )
    ! itau is the 1st day of new month at 00:00 so we need to fix the reftime back half a month (15th day)
    ! ((cursecond - reftimesecond) / seconds_in_day) - days_in_last_month + 15days
    
    !reftime = (itau - itauref -n_days_in_month*24*3600 + 15*24*3600) / 86400.  
    reftime = (itau - itauref ) / 86400.  
    !get current date in integers
    call tau2date(itau, curdate)
    ! create a TDATE date variable of the previous month (curdate(3)-1)
    
    curdate_tdate=newdate(curdate(1),curdate(2),curdate(3)-1,curdate(4),curdate(5),curdate(6),calender='greg')
    ! get days in month and save for next step
    n_days_in_month=days_in_month(curdate_tdate)
    ! change reftime to beginning of last month (the month data is from)
    reftime=reftime-n_days_in_month
    !length of the month-1s(in days) for the time_bounds
    rangemon=n_days_in_month !-(1.0/86400.0)

    ! allocate 3D and 4D global arrays for gathering data
    ! only root needs the full array, but it must be allocated in all tasks
    if (isRoot) then
       allocate( arr3d(imr,jmr,lmr)       )
       allocate( arr3dh(imr,jmr,lmr+1)       )
       allocate( arr2d(imr,jmr)       )
    else
       allocate( arr3d(1,1,1)   )
       allocate( arr3dh(1,1,1)   )
       allocate( arr2d(1,1)   )
    endif
    arr2d(:,:)=0.0
    arr3d(:,:,:)=0.0
    arr3dh(:,:,:)=0.0
    do j_var =1,n_vars
       ! hourly and daily variables are handled separately
       if (allvars(j_var)%table_id=='AERhr'.or.allvars(j_var)%table_id=='AER6hr'.or.&
            allvars(j_var)%table_id=='AERday'.or.allvars(j_var)%table_id=='AERfx')then

          cycle
       end if
       if (allvars(j_var)%dims==2)then !2D
          

          if (trim(allvars(j_var)%vname)/='minpblz'.and.trim(allvars(j_var)%vname)/='tasmin'.and. &
               trim(allvars(j_var)%vname)/='maxpblz'.and.trim(allvars(j_var)%vname)/='tasmax')then

             allvars(j_var)%data2D(i1:i2,j1:j2)=allvars(j_var)%data2D(i1:i2,j1:j2)/real(accumulation_mon)
          end if
          call GO_Timer_Start( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)
          
          call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0,status)
          call GO_Timer_End( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)
          
          IF_NOTOK_RETURN(status=1)
          
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1,timeidx_mon/), &
               count=(/imr,jmr,1/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          ! Zero out the accumulated and written data for the next interval
          if (trim(allvars(j_var)%vname)=='minpblz' .or. trim(allvars(j_var)%vname)=='tasmin') then
             ! put high number so simple comparison is needed for finding minimum
             allvars(j_var)%data2D(i1:i2,j1:j2)=1e10
          else
             allvars(j_var)%data2D(i1:i2,j1:j2)=0.0
          end if
       else !3D
          if (trim( allvars(j_var)%vname)=='phalf') then !lmr+1
             allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr+1)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr+1)/real(accumulation_mon) 
             call GO_Timer_Start( itim_write_gather, status )
             IF_NOTOK_RETURN(status=1)
             call gather( dgrid(1), allvars(j_var)%data3D , arr3dh(:,:,:), 0, status)
             IF_NOTOK_RETURN(status=1)
             call GO_Timer_End( itim_write_gather, status )
             IF_NOTOK_RETURN(status=1)

             if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr3dh , status, start=(/1,1,1,timeidx_mon/), &
                  count=(/imr,jmr,lmr+1,1/) )
             IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          
          else
             allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)/real(accumulation_mon)
             call GO_Timer_Start( itim_write_gather, status )
             IF_NOTOK_RETURN(status=1)
             call gather( dgrid(1), allvars(j_var)%data3D , arr3d(:,:,:), 0, status)
             IF_NOTOK_RETURN(status=1)
             call GO_Timer_End( itim_write_gather, status )
             IF_NOTOK_RETURN(status=1)
             
             if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr3d , status, start=(/1,1,1,timeidx_mon/), &
                  count=(/imr,jmr,lmr,1/) )
             IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
             
          end if

          ! Zero out the accumulated and written data for the next interval
          allvars(j_var)%data3D(i1:i2,j1:j2,:)=0.0
          
       end if
       !end if
       ! write the date for timestep
       if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime+real(rangemon/2)/) , status, start=(/timeidx_mon/), count=(/1/) )
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%bounds_varid,(/ reftime,reftime+rangemon/) , status, &
            start=(/1,timeidx_mon/), count=(/2,1/) )
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

       

    end do
    deallocate( arr3d,arr3dh,arr2d)
    ! change time index
    timeidx_mon= timeidx_mon + 1
    
    ! accululated time to zero
    accumulation_mon=0

    call GO_Timer_End( itim_write_mon, status )
    IF_NOTOK_RETURN(status=1)


    call goLabel()
    status = 0
  end subroutine output_c4mip_write_monthly
 subroutine output_c4mip_write_daily(region,newday,status)
    use MeteoData,      only : temper_dat, sst_dat, albedo_dat, ci_dat, lsp_dat, cp_dat, ssr_dat, str_dat, &
         blh_dat, gph_dat, lwc_dat, iwc_dat, cco_dat, cc_dat, humid_dat, m_dat, phlb_dat, sp_dat   !
    use meteodata   , only : global_lli, levi
    use partools    , only : isRoot,myid
    use GO,             only : TDate, NewDate!
    use datetime,       only : tau2date,date2tau,julday   !
    use dims,           only : itau,iyear0          !current time
    use tm5_distgrid, only : dgrid,    Get_DistGrid ,gather
    implicit none
    logical,intent(in) ::newday
    integer,intent(out)::status
    integer::region
    integer:: imr,jmr,i,i1,i2,j1,j2,lmr
    character(len=*), parameter ::  rname = mname//'/output_c4mip_write_daily'
    integer:: j_var

    ! reference time:
    integer, parameter            ::  time_reftime6(6) = (/1750,01,01,00,00,00/)
    integer(kind=8)                       ::  itauref  ! reftime in seconds
    real                          ::  reftime  ! seconds from reference time
    real                          :: rangeday ! for bounds
    ! root writes netcdf arrays
    real, allocatable :: arr3d(:,:,:), arr2d(:,:)    
4    integer:: imr_f,jmr_f,lmr_f

    call goLabel(rname)

    call GO_Timer_Start( itim_write_day, status )
    IF_NOTOK_RETURN(status=1)

    if (region > 1) then
       write(gol,*) 'output_c4mip_write_daily: region >1, only available in global mode!'
       call goErr
       status=1
       return
    end if


    call Get_DistGrid( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2)
    

    ! entire region grid size
    imr = global_lli(1)%nlon
    jmr = global_lli(1)%nlat
    lmr = levi%nlev
    
    ! allocate 3D and 4D global arrays for gathering data
    if (isRoot) then
       allocate( arr3d(imr,jmr,lmr)       )
       allocate( arr2d(imr,jmr)       )
    else
       allocate( arr3d(1,1,1)   )
       allocate( arr2d(1,1)   )
    endif
    arr2d(:,:)=0.0
    arr3d(:,:,:)=0.0

    ! define the reference time in seconds (itauref)
    ! (for now in days since 1850-01-01 00:00, local variable) 
    
    call date2tau( time_reftime6, itauref )
    ! calculate reftime as fractional days from itauref, hence division by 86400
    ! call date2tau( idater, itaucur )
    reftime = (itau - itauref) / 86400. - 1.0
    !23h59 as days
    rangeday=1.0!(23.0*3600.0+59.0*60.0+59.0)/86400.0
    do j_var =1,n_vars
       if (allvars(j_var)%table_id/='AERday')then 
          cycle
       end if
       if (allvars(j_var)%dims==2)then
          
          if ( trim(allvars(j_var)%vname)/='minpblz' .and. trim(allvars(j_var)%vname)/='tasmin'.and.trim(allvars(j_var)%vname)/='maxpblz'.and. trim(allvars(j_var)%vname)/='tasmax')then
             allvars(j_var)%data2D(i1:i2,j1:j2)=allvars(j_var)%data2D(i1:i2,j1:j2)/real(accumulation_day)
          end if
          call GO_Timer_Start( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)
          call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0, status)
          IF_NOTOK_RETURN(status=1)
          call GO_Timer_End( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)

          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1,timeidx_day/), count=(/imr,jmr,1/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          
          if (trim(allvars(j_var)%vname)=='minpblz' .or. trim(allvars(j_var)%vname)=='tasmin') then
             ! put high number so simple comparison is needed for finding minimum
             allvars(j_var)%data2D(i1:i2,j1:j2)=1e10
          else
             ! Zero out the accumulated and written data for the next interval
             allvars(j_var)%data2D(i1:i2,j1:j2)=0.0
          end if
       else
          
          allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)/real(accumulation_day)
          !end if
          call GO_Timer_Start( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)
          call gather( dgrid(1), allvars(j_var)%data3D , arr3d(:,:,:), 0, status)
          call GO_Timer_End( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)          !if (trim(allvars(j_var)%vname)=='od5503ddust')then
          IF_NOTOK_RETURN(status=1)
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr3d, status, start=(/1,1,1,timeidx_day/), count=(/imr,jmr,lmr,1/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          
          
          allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=0.0
          
       end if

       
       ! write the date for timestep
        if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime+0.5/) , status, start=(/timeidx_day/), count=(/1/) )
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%bounds_varid,(/ reftime,reftime+ rangeday/) , status, start=(/1,timeidx_day/), count=(/2,1/) )
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

    end do
    deallocate(arr3d, arr2d)
    ! Timeindex to next day
    timeidx_day= timeidx_day + 1
    ! daily accumulated time to zero
    accumulation_day=0.0
    !status=1
    !return
    call GO_Timer_End( itim_write_day, status )
    IF_NOTOK_RETURN(status=1)


    call goLabel()
    status = 0
  end subroutine output_c4mip_write_daily


 subroutine output_c4mip_write_hourly(region,newhour,status)
    use MeteoData,      only : temper_dat,sst_dat,albedo_dat,ci_dat,lsp_dat,cp_dat,ssr_dat,str_dat,blh_dat,gph_dat,lwc_dat,iwc_dat,cco_dat,cc_dat,humid_dat, m_dat,phlb_dat,sp_dat   !
    use GO,             only : TDate, NewDate!
    use datetime,       only : tau2date,date2tau,julday   !
    use dims,           only : itau,iyear0          !current time
    use tm5_distgrid, only : dgrid,    Get_DistGrid ,gather
    use partools    , only : isRoot,myid
    implicit none
    logical,intent(in) ::newhour
    integer,intent(out)::status
    integer:: j_var
    integer::region
    integer:: imr,jmr,i,i1,i2,j1,j2,lmr
    character(len=*), parameter ::  rname = mname//'/output_c4mip_write_hourly'
    real                          :: rangehr,range6hr ! hour in days for bounds in NC-file

    ! reference time:
    integer, parameter            ::  time_reftime6(6) = (/1750,01,01,00,00,00/)
    integer(kind=8)                       ::  itauref  ! reftime in seconds
    real                          ::  reftime  ! seconds from reference time
    ! root writes netcdf arrays
    real, allocatable :: arr3d(:,:,:) , arr2d(:,:)
    call goLabel(rname)

    call GO_Timer_Start( itim_write_hour, status )
    IF_NOTOK_RETURN(status=1)

    if (region > 1) then
       write(gol,*) 'output_c4mip_write_hourly: region >1, only available in global mode!'
       call goErr
       status=1
       return
    end if

    call Get_DistGrid( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2)
    
    ! entire region grid size
    imr = global_lli(1)%nlon
    jmr = global_lli(1)%nlat
    lmr = levi%nlev
    
    ! allocate 3D and 4D global arrays for gathering data
    if (isRoot) then
       allocate( arr3d(imr,jmr,lmr)       )
       allocate( arr2d(imr,jmr)       )
    else
       ! other than root need the variable, but no space
       allocate( arr3d(1,1,1)   )
       allocate( arr2d(1,1)   )
    endif
    arr2d(:,:)=0.0
    arr3d(:,:,:)=0.0

    ! define the reference time in seconds (itauref)
    ! (for now in days since 1850-01-01 00:00, local variable) 
    
    call date2tau( time_reftime6, itauref )
    ! call date2tau( idater, itaucur )
    
    rangehr=1.0/24.0!(3600)/86400.0
    do j_var =1,n_vars
       
       if (allvars(j_var)%table_id/='AERhr' .and. allvars(j_var)%table_id/='AER6hr' )then 
          cycle
       end if
       if (allvars(j_var)%dims==0 .and.trim(allvars(j_var)%vname)=='co2mass1' )then
          reftime = (itau - itauref) / 86400. - (1./24)
          allvars(j_var)%data2D(i1:i2,j1:j2)=allvars(j_var)%data2D(i1:i2,j1:j2)/real(accumulation_hr)
          call gather( dgrid(1), allvars(j_var)%data2D , arr3d(:,:,1), 0, status)
          IF_NOTOK_RETURN(status=1)
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid,(/sum(arr3d(:,:,1))/), status, start=(/timeidx_hr/), count=(/1/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       else if (allvars(j_var)%dims==2)then
          if ( trim(allvars(j_var)%table_id)=='AERhr') then
             reftime = (itau - itauref) / 86400. - (1./24)
             
               allvars(j_var)%data2D(i1:i2,j1:j2)=allvars(j_var)%data2D(i1:i2,j1:j2)/real(accumulation_hr)
               call GO_Timer_Start( itim_write_gather, status )
               IF_NOTOK_RETURN(status=1)
               
               call gather( dgrid(1), allvars(j_var)%data2D , arr3d(:,:,1), 0, status)
               IF_NOTOK_RETURN(status=1)
               call GO_Timer_End( itim_write_gather, status )
               IF_NOTOK_RETURN(status=1)               

               if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid,arr3d(:,:,1), status, start=(/1,1,timeidx_hr/), count=(/imr,jmr,1/) )
               IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
               
               
               ! write the date for timestep
                if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime+(0.5/24.0)/) , status, start=(/timeidx_hr/), count=(/1/) )
                IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
                if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%bounds_varid,(/ reftime,reftime+rangehr/) , status, start=(/1,timeidx_hr/), count=(/2,1/) )
                IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
               ! Zero out the accumulated and written data for the next interval
               allvars(j_var)%data2D(i1:i2,j1:j2)=0.0
            end if
       else

          if ( trim(allvars(j_var)%table_id)=='AERhr') then
             reftime = (itau - itauref) / 86400. - (1./24)
             allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)/real(accumulation_hr)
             call GO_Timer_Start( itim_write_gather, status )
             IF_NOTOK_RETURN(status=1)
             call gather( dgrid(1), allvars(j_var)%data3D , arr3d(:,:,:), 0, status)
             call GO_Timer_End( itim_write_gather, status )
             IF_NOTOK_RETURN(status=1)
             if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid,arr3d, status, start=(/1,1,1,timeidx_hr/), count=(/imr,jmr,lmr,1/) )
             
             ! write the date for timestep
              if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime+(0.5/24.0)/) , status, start=(/timeidx_hr/), count=(/1/) )
             IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
             if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%bounds_varid,(/ reftime,reftime+(1./24.)/) , status, start=(/1,timeidx_hr/), count=(/2,1/) )
             IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

             ! Zero out the accumulated and written data for the next interval
             allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=0.0
          end if
       end if
       
       !end if

    end do
    deallocate(arr3d, arr2d)
    ! changed index to next hour
    timeidx_hr= timeidx_hr + 1
    ! accumulated hours to zero, actually this will always be 1h
    accumulation_hr=0.0
    
    
    !status=1
    !return
    call GO_Timer_End( itim_write_hour, status )
    IF_NOTOK_RETURN(status=1)
    

    call goLabel()
    status = 0
  end subroutine output_c4mip_write_hourly
 subroutine output_c4mip_write_6hourly(region,newhour,status)
    use MeteoData,      only : temper_dat,sst_dat,albedo_dat,ci_dat,lsp_dat,cp_dat,ssr_dat,str_dat,blh_dat,gph_dat,lwc_dat,iwc_dat,cco_dat,cc_dat,humid_dat, m_dat,phlb_dat,sp_dat   !
    use GO,             only : TDate, NewDate!
    use datetime,       only : tau2date,date2tau,julday   !
    use dims,           only : itau,iyear0          !current time
    use tm5_distgrid, only : dgrid,    Get_DistGrid ,gather
    use partools    , only : isRoot,myid
    !use ebischeme,      only : AC_diag_prod,iprod_soa2dhour
    implicit none
    logical,intent(in) ::newhour
    integer,intent(out)::status
    integer::region
    integer:: j_var
    integer:: imr,jmr,i,i1,i2,j1,j2,lmr
    character(len=*), parameter ::  rname = mname//'/output_c4mip_write_6hourly'
    
    ! reference time:
    integer, parameter            ::  time_reftime6(6) = (/1750,01,01,00,00,00/)
    integer(kind=8)                       ::  itauref  ! reftime in seconds
    real                          ::  reftime  ! seconds from reference time
    ! root writes netcdf arrays
    real, allocatable :: arr3d(:,:,:) , arr2d(:,:)
    call goLabel(rname)

    call GO_Timer_Start( itim_write_hour, status )
    IF_NOTOK_RETURN(status=1)

    if (region > 1) then
       write(gol,*) 'output_c4mip_write_6hourly: region >1, only available in global mode!'
       call goErr
       status=1
       return
    end if

    call Get_DistGrid( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2)
    
    ! entire region grid size
    imr = global_lli(1)%nlon
    jmr = global_lli(1)%nlat
    lmr = levi%nlev
    
    ! allocate 3D and 4D global arrays for gathering data
    if (isRoot) then
       allocate( arr3d(imr,jmr,lmr)       )
       allocate( arr2d(imr,jmr)       )
    else
       ! other than root need the variable, but no space
       allocate( arr3d(1,1,1)   )
       allocate( arr2d(1,1)   )
    endif
    arr2d(:,:)=0.0
    arr3d(:,:,:)=0.0

    ! define the reference time in seconds (itauref)
    ! (for now in days since 1850-01-01 00:00, local variable) 
    
    call date2tau( time_reftime6, itauref )
    ! call date2tau( idater, itaucur )
    reftime = (itau - itauref) / 86400. 
    
    do j_var =1,n_vars
       
       if ( allvars(j_var)%table_id/='AER6hr' )then 
          cycle
       end if
       if (allvars(j_var)%dims==2)then
          !allvars(j_var)%data2D(i1:i2,j1:j2)=allvars(j_var)%data2D(i1:i2,j1:j2)/real(accumulation_6hr)
             call GO_Timer_Start( itim_write_gather, status )
             IF_NOTOK_RETURN(status=1)
             call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0, status)
             IF_NOTOK_RETURN(status=1)
             call GO_Timer_End( itim_write_gather, status )
             IF_NOTOK_RETURN(status=1)
             if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid,arr2d, status, start=(/1,1,timeidx_6hr/), count=(/imr,jmr,1/) )
             IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

             ! write the date for timestep
             if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime/) , status, start=(/timeidx_6hr/), count=(/1/) )
             IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
             
             ! Zero out the accumulated and written data for the next interval
             allvars(j_var)%data2D(i1:i2,j1:j2)=0.0
          
       else
        
             !allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)/real(accumulation_6hr)
             call GO_Timer_Start( itim_write_gather, status )
             IF_NOTOK_RETURN(status=1)
             call gather( dgrid(1), allvars(j_var)%data3D , arr3d(:,:,:), 0, status)
             call GO_Timer_End( itim_write_gather, status )
             IF_NOTOK_RETURN(status=1)
             if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid,arr3d, status, start=(/1,1,1,timeidx_6hr/), count=(/imr,jmr,lmr,1/) )
             
             ! write the date for timestep
             if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime/) , status, start=(/timeidx_6hr/), count=(/1/) )
             IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

             !start=(/i1,j1,1,timeidx_mon/), count=(/imr,jmr,lmr,1/) )
             IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
             ! Zero out the accumulated and written data for the next interval
             allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=0.0

       end if
    end do
    deallocate(arr3d,arr2d)
    ! changed index to next 6hour
    timeidx_6hr= timeidx_6hr + 1
    ! exception for one 6hr field, no need to do another subroutine for it
    accumulation_6hr=0.0

    
    call GO_Timer_End( itim_write_hour, status )
    IF_NOTOK_RETURN(status=1)

    call goLabel()
    status = 0
  end subroutine output_c4mip_write_6hourly

  subroutine accumulate_c4mip_data(dhour,status)
    use GO         , only : TDate, NewDate, rTotal, operator(-)
    use Grid       , only : FPressure,HPressure
    use binas,       only : rgas, rol,xmair,Dobs,Avog
    USE toolbox,       only : ltropo_ifs,     lvlpress
    !use datetime,    only : tau2date
    use MeteoData,      only : temper_dat, sst_dat, albedo_dat, ci_dat, lsp_dat, cp_dat, ssr_dat, str_dat, blh_dat, &
         gph_dat, lwc_dat, iwc_dat, cco_dat, cc_dat, humid_dat, m_dat, phlb_dat, sp_dat, pu_dat, pv_dat,pw_dat
    !use photolysis_data,only:phot_dat   !
    use global_data, only : mass_dat, region_dat,conv_dat
    use dims,        only : lm,sec_month
    use partools    , only : isRoot,myid
    use dims,    only: gtor, dx, dy, ybeg, xref, yref,ndyn
    use binas,   only: ae

    implicit none

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

    !    integer                 :: indices(7)
    integer                 :: itrac,gasind
    integer                 :: i_sed
    integer                 :: i_emi
    integer                 :: i, j, k, imr, jmr, lmr, lwl, dtime,index,imode,m
    integer                 :: i1, i2, j1, j2
    integer                 :: ie, je ! indices for subdomain extended with halo cells
    integer,intent(in)      :: dhour
    integer                 :: status
    integer                 :: j_var,region,i_var,i_wdep,sedindex,icomp

    real                    :: dens
    real                    :: vol

    real                    :: tempbud,xmcomp,temp
    real, dimension(:,:,:,:), pointer :: tracers     ! transported tracers
    real, dimension(:),       pointer :: dxyp      
    integer, dimension(n_indices) :: indices
    real::xmgas

    real, dimension(:,:,:), pointer :: t     ! temperature (K)
    !

    call goLabel(rname)

    call GO_Timer_Start( itim_accu, status )
    IF_NOTOK_RETURN(status=1)
    region=1
    ! grid size
    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2)
    imr = i2-i1+1
    jmr = j2-j1+1
    lmr = levi%nlev

    t           => temper_dat (region)%data
    !accumulation_6hr=0.0!accumulation_6hr+dtime
    ! Gridbox area
    dxyp   => region_dat(region)%dxyp

    ! mass_dat keep data in kg/gridbox
    ! 
    tracers => mass_dat(region)%rm 

    dtime = dhour*3600
    accumulation_mon=accumulation_mon+dtime
    accumulation_hr=accumulation_hr+dtime
    accumulation_day=accumulation_day+dtime
    do j_var = 1, n_vars
       indices(:)=allvars(j_var)%indices(:)
       !Here we use the tm5-indices to collect data for output
  
       if (trim(allvars(j_var)%class)=='fixed') then
          cycle
  
       else if (trim(allvars(j_var)%class)=='monthly')then
          index=indices(1)
          do k=1,lmr
             do j=j1,j2
                do i=i1,i2
                   if (trim(allvars(j_var)%vname)=='ps')then
                      if (k .eq. 1) allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*sp_dat(1)%data(i,j,1)!Pa
                   else if(trim(allvars(j_var)%vname)=='co2')then
                      xmcomp=ra(index)
                      allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)/m_dat(region)%data(i,j,k)*xmair/xmcomp
                   else if(trim(allvars(j_var)%vname)=='co2mass')then
                      index= allvars(j_var)%indices(1)
                      allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,k,index)/dxyp(j)! kg/m2
                   else if (index<=0) then
                      ! you should not end up here!!!
                      cycle
                   end if
                end do
             end do
          end do
       else if (trim(allvars(j_var)%class)=='ps6h')then
          do i=i1,i2
             do j=j1,j2
                if (trim(allvars(j_var)%vname)=='ps')then
                   allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*sp_dat(1)%data(i,j,1)!Pa
                end if
             end do
          end do
       ! 1 hourly surface variables
       else if (trim(allvars(j_var)%class)=='hourly')then
          do i=i1,i2
             do j=j1,j2
                if (trim(allvars(j_var)%vname)=='ps')then
                   allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*sp_dat(1)%data(i,j,1)!Pa

                else if (trim(allvars(j_var)%vname)=='co2') then

                   index= indices(1)
                   xmcomp=ra(index)
                   do k=1,lmr
                      allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)/m_dat(region)%data(i,j,k)*xmair/xmcomp
                   end do
                else if(trim(allvars(j_var)%vname)=='co2mass')then
                   do k=1,lmr
                      index= allvars(j_var)%indices(1)

                      allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,k,index)/dxyp(j)!kg/m2
                   end do
                else if(trim(allvars(j_var)%vname)=='co2mass1')then
                   do k=1,lmr
                      index= allvars(j_var)%indices(1)
                      allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,k,index)/dxyp(j)!kg/m2
                   end do
                else if (trim(allvars(j_var)%vname)=='tas')then
                   allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*temper_dat(1)%data(i,j,1)!K

                end if
             end do
          end do
       ! surface daily variables
       else if (trim(allvars(j_var)%class)=='daily')then
          index= indices(1)
          do i=i1,i2
             do j=j1,j2
                if (trim(allvars(j_var)%vname)=='ps')then
                   allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*sp_dat(1)%data(i,j,1)!Pa
                else if (trim(allvars(j_var)%vname)=='co2') then

                   xmcomp=ra(index)
                   do k=1,lmr
                      allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)/m_dat(region)%data(i,j,k)*xmair/xmcomp

                   end do
                else if(trim(allvars(j_var)%vname)=='co2mass')then
                   do k=1,lmr
                      allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,k,index)/dxyp(j)! kg/m2
                   end do
                end if
             end do
          end do
       else
          write(gol,*) 'output_c4mip_accumulate: output class not found!!!',trim(allvars(j_var)%vname),trim(allvars(j_var)%class)
          !call goPr
          call goErr
          status=1
          return
          
       end if
    end do

    ! zero accumulated budget variables for the amount between output steps


    call GO_Timer_End( itim_accu, status )
    IF_NOTOK_RETURN(status=1)
    call goLabel()
    !status = 1
  end subroutine accumulate_c4mip_data

  subroutine output_c4mip_done(status)
    use partools, only: isRoot,myid
    implicit none
    integer :: status
    character(len=*), parameter ::  rname = mname//'/output_c4mip_done'
    integer :: j_var, region

    call goLabel(rname)
    region = 1
    if (isroot) then
       do j_var=1,n_vars
          call MDF_Close( allvars(j_var)%fileunit, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       end do
    end if


    deallocate(dimension_data%lon) 
    deallocate(dimension_data%lat) 
    deallocate(dimension_data%lev) 
 
    do j_var=1,n_vars
      deallocate(allvars(j_var)%data2D)
      deallocate(allvars(j_var)%data3D)
    end do
    deallocate(allvars)
    deallocate(fixedvars)
    
    call goLabel()
    status = 0
  end subroutine output_c4mip_done

  subroutine add_variable(itm5,shortname,longname,unit,data_dims,status,class,table,pxmgas)
#ifdef with_m7
    use chem_param,  only: mode_end_so4,mode_end_pom,mode_end_bc,mode_end_ss,mode_end_dust
#else
    use chem_param, only: ico2
#endif
    use global_data, only: outdir
    use datetime,    only: tau2date, date2tau
    use dims,        only: itau,itaue,itaut
  
    implicit none
    integer:: itm5
    character(len=*),intent(in)::shortname
    character(len=*),intent(in)::longname
    character(len=*)::unit
    character(len=30)::standardname

    character(len=*)::table
    character(len=*),optional::class
    real,optional::pxmgas
    integer:: data_dims
    integer,intent(out)::status
    !LOCAL
    character(len=4)::positive=''
    integer:: fileunit=-1 !defined only when creating a file
    integer:: varid=-1! defined when opening ncfile
    !character(len=120)::filename
    character(len=30)::table_id
    !character(len=30)::c4mip_source_id
    !character(len=30)::c4mip_experiment_id
    character(len=30)::member_id
    !character(len=30)::grid_label
    character(len=30)::time_range
    character(len=200)::filename1
    character(len=10)::freq
    real,dimension(:,:),pointer::data2D
!    real,dimension(:,:),pointer::dataZonal
    real,dimension(:,:,:),pointer::data3D
!    real,dimension(:,:,:),pointer::budget

    character(len=*), parameter ::  rname = mname//'/output_c4mip_add_variable'
    integer ::i1,i2,j1,j2,jmr,imr,lmr
    integer, dimension(6) :: idater, idateend, idatett
    integer               :: endmonth, endday
    character(len=30)     :: idates

    call tau2date(itau,idater)

    ! define frequency from table
    if (trim(table)=='AERhr')then
    !table id depends on variable
       table_id=table
       freq='1hr'
    else if (trim(table)=='AER6hr')then
    !table id depends on variable
       table_id=table
       freq='6hr'
    else if( trim(table)=='AERmon'.or.trim(table)=='AERmonZ'.or.trim(table)=='Emon')then
       !table id depends on variable
       table_id=table
       freq='mon'
    else if(trim(table)=='AERday')then
       !table id depends on variable
       table_id=table
       freq='day'
    else if(trim(table)=='AERfx')then
       !table id depends on variable
       table_id=table
       freq='na'
    else
       freq='freq-nondefined'
       table_id='table-nondefined'
    end if
    
    ! CREATE date string for output
    ! 
    ! ATM assumed that the output is initilised at the beginninh of year
    endmonth=12
    endday=31
    !
    if (freq=='mon')then
       ! By default CO2 runs are done by 1-year chunks -> idater(2) -  idater(2)+11
       write(idates,  '(i4,i2.2,a,i4,i2.2)') idater(1), idater(2),'-', idater(1),endmonth
    else if(freq=='day')then
       !time range form Jan-1 ->Dec-31x
       write(idates,  '(i4,2i2.2,a,i4,2i2.2)') idater(1), idater(2), idater(3),'-', idater(1), endmonth, endday

    else if(freq=='1hr')then
       write(idates,  '(i4,2i2.2,2a2,a,i4,2i2.2,2a2)') idater(1), idater(2), idater(3),'00','00','-', idater(1), endmonth, endday, '23', '59'    
    else if (freq=='6hr')then

       write(idates,  '(i4,2i2.2,2a2,a,i4,2i2.2,2a2)') idater(1), idater(2), idater(3),'00','00','-', idater(1), endmonth, endday,'18','00'    
    end if
       call goLabel(rname)
    
    call GO_Timer_Start( itim_addvar, status )
    IF_NOTOK_RETURN(status=1)
    !outdir='output'


    ! temporary
    standardname=longname
    ! c4mip_source_id constant
    !c4mip_source_id='EC-EARTH-CC'
    ! experiment depends on run
    !experiment_id='exp_dynamic'
    member_id='r'//trim(realization_i)//'i'//trim(initialization_i)//'p'//trim(physics_i)//'f'//trim(forcing_i)
    !grid_label='3x2_degrees'
    
    ! time range has divider in place since it can be omitted alltogether with non-time dependendent variables
    if (trim(table)=='AERfx')then
       time_range=''
    else
       time_range='_'//trim(idates)
    end if

    ! for fixed variables time range should not be written
    n_vars=n_vars+1
    call Get_DistGrid( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

    ! define sizes for arrays
    imr=i2-i1+1
    jmr=j2-j1+1
    lmr = levi%nlev
!    if (trim(shortname)=='phalf') then
!       allocate(budget(i1:i2,j1:j2,1:lmr+1))
!    else
!       allocate(budget(i1:i2,j1:j2,1:lmr))
!    end if
!    budget(:,:,:)=0.0
    !write(2004,*)shortname
    ! 2D variables
    if (data_dims==2) then
       !Allocate holders for data
       allocate(allvars(n_vars)%data2D(i1:i2,j1:j2))
       allocate(allvars(n_vars)%data3D(1,1,1))

       ! allocate local variables
       allocate(data2D(i1:i2,j1:j2))
       allocate(data3D(1,1,1))
       ! zero local data holders
       data2D(:,:)=0.0
       !       dataZonal(:,:)=0.0
       data3D(:,:,:)=0.0
       !init variable

       ! set default for minimum variables to high value
       if (shortname=='minpblz' .or. shortname=='tasmin')then
          data2D(:,:)=1000000.0
       end if
       !create filename
     
       if (trim(class)=='crescendo')then
 
!          filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(table_id)//'_'//trim(c4mip_source_id)//'_'//trim(experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//'_'//trim(time_range)//trim('.nc')
          filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(class)//'_'//trim(table_id)//'_'//trim(c4mip_source_id)//'_'//trim(c4mip_experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//trim(time_range)//trim('.nc')
       else
          filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(table_id)//'_'//trim(c4mip_source_id)//'_'//trim(c4mip_experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//trim(time_range)//trim('.nc')
       end if
       allvars(n_vars)=varfile(itm5,shortname,longname,unit,positive,standardname,data2D,data3D,-1,-1,&
            filename1,2,-1,-1,-1,-1,-1,-1,-1,-1,-1,data_dims,freq,class,(/0,0,0,0,0,0,0,0,0,0,0/),pxmgas,table_id)


!! LEFT HERE on purpose to see what variables go where in above statement
!!$       allvars(n_vars)%itm5=itm5
!!$       allvars(n_vars)%vname=shortname
!!$       allvars(n_vars)%lname=longname
!!$       allvars(n_vars)%unit=unit
!!$       allvars(n_vars)%positive=positive
!!$       allvars(n_vars)%standard_name=standardname
!!$       allvars(n_vars)%data2D=data2D
!!$       allvars(n_vars)%data3D=data3D
!!$       allvars(n_vars)%budgetdata=data3D
!!$       allvars(n_vars)varid=-1
!!$       allvars(n_vars)%filenunit=-1
!!$       allvars(n_vars)%filename=filename1
!!$       allvars(n_vars)%dimensions=2
!!$       allvars(n_vars)%lon_varid=-1
!!$       allvars(n_vars)%lat_varid=-1
!!$       allvars(n_vars)%lev_varid=-1
!!$       allvars(n_vars)%time_varid=-1
!!$       allvars(n_vars)%bounds_varid=-1
!!$       allvars(n_vars)%dims=dims
!!$       allvars(n_vars)%class=class
!!$       allvars(n_vars)%indices=(/0,0,0,0,0,0,0/)) 
!!$       allvars(n_vars)%xmgas=molarweight
!!$       allvars(n_vars)%table_id=


    ! 3D variables
    else if (data_dims==3) then
       ! allocate holders for data
       allocate(allvars(n_vars)%data2D(1,1))
       
       if (trim(shortname)=='phalf') then
          allocate(allvars(n_vars)%data3D(i1:i2,j1:j2,1:lmr+1))
          allocate(data3D(i1:i2,j1:j2,1:lmr+1))
  
       else
          allocate(allvars(n_vars)%data3D(i1:i2,j1:j2,1:lmr))
          allocate(data3D(i1:i2,j1:j2,1:lmr))
          
   
       end if
  

       ! allocate local variables
       ! maybe remove these
       allocate(data2D(1,1))
       !allocate(data3D(i1:i2,j1:j2,1:lmr))
       ! zero local data holders
       data2D(:,:)=0.0
       data3D(:,:,:)=0.0

       !init variable

         filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(table_id)//'_'//trim(c4mip_source_id)//'_'//trim(c4mip_experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//trim(time_range)//trim('.nc')


         allvars(n_vars)=varfile(itm5,shortname,longname,unit,positive,standardname,data2D,data3D,-1,-1,&
            filename1,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,data_dims,freq,class,(/0,0,0,0,0,0,0,0,0,0,0/),pxmgas,table)
    end if
    ! add chemical info also: (vars beginning with emi,wet,dry)
    select case (trim(shortname(4:)))

    case ('so2')
       allvars(n_vars)%indices(1)=iso2

    end select
    select case (trim(shortname))

    case('areacella')
        allvars(n_vars)%indices(:)=0
        areacella=n_vars


     case('co2','co2mass','co2mass1')
        write(2000,*) ico2,n_vars,trim(shortname),table
        
        allvars(n_vars)%indices(1)=ico2


    end select


      call goLabel()
      status = 0
      call GO_Timer_End( itim_addvar, status )
      IF_NOTOK_RETURN(status=1)
    end subroutine add_variable

  subroutine write_attributes(status,j_var)
    implicit none
    integer     :: j_var
    integer,intent(out)::status
    character(len=*), parameter ::  rname = mname//'/output_c4mip_writeattr'
    call goLabel(rname)
    call GO_Timer_Start( itim_attr, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att( allvars(j_var)%fileunit, MDF_GLOBAL, 'title',               'Model output for C4mip', status )
    IF_NOTOK_MDF(fid= allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lon_varid , 'units',          'degrees_east', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lon_varid , 'axis',          'X', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lon_varid , 'long_name',          'longitude', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lon_varid , 'standard_name',          'longitude', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lat_varid , 'units',          'degrees_north', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lat_varid , 'axis',          'Y', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lat_varid , 'long_name',          'latitude', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lat_varid , 'standard_name',          'latitude', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
    ! allvars(j_var)%lev_varid
     if (allvars(j_var)%dims==3) then
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'units',          'level', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'axis',          'Z', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'long_name',          'hybrid model level at layer midpoints', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'standard_name',          'hybrid_model_level', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    

       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'formula',          'a+b*ps', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    

       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'positive',          'up', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
       

       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hyam_varid , 'long_name',          'hybrid A coefficient at layer midpoints', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hyam_varid , 'units',          'Pa', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hybm_varid , 'long_name',          'hybrid B coefficient at layer midpoints', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hybm_varid , 'units',          '1', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    

       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hyai_varid , 'long_name',          'hybrid A coefficient at layer interfaces', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hyai_varid , 'units',          'Pa', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hybi_varid , 'long_name',          'hybrid B coefficient at layer interfaces', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hybi_varid , 'units',          '1', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
    END if
    if (allvars(j_var)%table_id/='AERfx')then
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'units',          'days since 1750-01-01 00:00:00', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       !    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'axis',          'X', status)
       !    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'calendar',          'gregorian', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'long_name',          'time', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'axis',          'T', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    

       !time bounds
       call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'bounds',          'time_bounds', status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
    end if
   
    !experiment=
    !CMIP6 global attributes:
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,    'Conventions', 'CF-1.7 CMIP-6.0 UGRID-0.9', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'activity_id',  'C4mip', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'branch_method','', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'creation_date','', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'data_specs_version','1.0.0', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'experiment',trim(c4mip_experiment), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'experiment_id',trim(c4mip_experiment_id), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'forcing_index',trim(forcing_i), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'frequency',trim(allvars(j_var)%freq), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'further_info_url','MISSING', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'grid','native '//trim(grid)//' degree grid', status)!module variables
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'grid_label',trim(grid_label), status)!module variables
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'initialization_index',trim(initialization_i), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'institution','KNMI, The Netherlands; SMHI, Sweden; DMI, Denmark; &
         &AEMET, Spain; Met Eireann, Ireland; CNR-ISAC, Italy; Instituto de Meteorologia, Portugal; FMI, Finland; BSC, Spain; &
         &Centro de Geofisica, University of Lisbon, Portugal; ENEA, Italy; Geomar, Germany; Geophysical Institute, University of Bergen, Norway; &
         &ICHEC, Ireland; ICTP, Italy; IMAU, The Netherlands; IRV, Sweden;  Lund University, Sweden; &
         &Meteorologiska Institutionen, Stockholms University, Sweden; Niels Bohr Institute, University of Copenhagen, Denmark; &
         &NTNU, Norway; SARA, The Netherlands; Unite ASTR, Belgium; Universiteit Utrecht, The Netherlands; &
         &Universiteit Wageningen, The Netherlands; University College Dublin, Ireland; Vrije Universiteit Amsterdam, the Netherlands; &
         &University of Helsinki, Finland; KIT, Karlsruhe, Germany; USC, University of Santiago de Compostela, Spain; &
         &Uppsala Universitet, Sweden; NLeSC, Netherlands eScience Center, The Netherlands', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'institution_id','EC-Earth-Consortium', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'license','NEEDS DEFINING', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'mip_era','CMIP6', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'nominal_resolution','250 km', status)!dmax
    !dmax=r*phi/2*(1+((phi**2+lamb**2)/(phi*lamb))*np.arctan(lamb/phi))=348 r=6371, phi=2(lat), lamb=3(long)
    !CMIP6 global attributes: 100 < dmax < 350 -> nominal resolution 250 km
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'physics_index',trim(physics_i), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'product','output', status)!only choice
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'realization_index',trim(realization_i), status)!1 for primary or single realization
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'realm',trim(c4mip_realm), status)!  depends on run, some are AGCM
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'source',trim(c4mip_source), status)!
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'source_id',trim(c4mip_source_id), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'source_type',trim(c4mip_source_type), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'sub_experiment','', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'sub_experiment_id','', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'table_id',trim(allvars(j_var)%table_id), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'tracking_id','', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'variable_id',trim(allvars(j_var)%vname), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'variant_label','', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'variant_label','', status)
    
    call GO_Timer_End( itim_attr, status )
    IF_NOTOK_RETURN(status=1)

    call goLabel() 
    status = 0
    
  end subroutine write_attributes
  
  subroutine write_dimensions(status,j_var)
    use dims,           only : iyear0          !current year
    use go_date,        only : days_in_year,newDate
    use partools    , only : isRoot,myid
    implicit none
    integer     :: j_var
    integer,intent(out)::status
    integer     :: i1,i2,j1,j2,imr,jmr,lmr
    integer     :: lon_varid,lonid,lon_dimid
    integer     :: lat_varid,latid,lat_dimid
    integer     :: lev_varid,levid,lev_dimid
    integer     :: hym_varid,hym_dimid
    integer     :: hyi_varid,hyi_dimid
    integer     :: time_varid,timeid,time_dimid,bounds_dimid,bounds_varid,boudid
    ! most of data  is monthly mean, but change to dynamic number of output steps needed
    integer     :: nout_steps!=12 
    integer     :: nhym
    integer     :: nhyi 

    character(len=*), parameter ::  rname = mname//'/output_c4mip_write_dim'
    call goLabel(rname)

    ! define dimensions
    !call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    imr=dimension_data%nlon
    jmr=dimension_data%nlat
    lmr=dimension_data%nlev
    nhym=lmr
    nhyi=lmr+1

    ! With parallel netcdf whole netcdf must be reserved at the time of initialisation
    ! therefore we need to know the number of output steps per file.

    ! Define number of output steps in a file depending on the output frequency
    ! use newdate to create TDate structure, and use that in days_in_year
    if (allvars(j_var)%table_id=='AERhr')then
       nout_steps=24*days_in_year(newdate(iyear0))
    else if (allvars(j_var)%table_id=='AER6hr')then
       nout_steps=4*days_in_year(newdate(iyear0))
    else if (allvars(j_var)%table_id=='AERday')then
       nout_steps=days_in_year(newdate(iyear0))
    else if (allvars(j_var)%table_id=='AERmon'.or. (allvars(j_var)%table_id=='AERmonZ'))then
       nout_steps=12
    end if

    if (isroot) then
       !DEFINE DIMENSIONS
       if (allvars(j_var)%dims>0) then
          call MDF_Def_Dim( allvars(j_var)%fileunit, 'lon',     imr,lon_dimid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          call MDF_Def_Dim( allvars(j_var)%fileunit, 'lat',      jmr, lat_dimid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       end if
       if (allvars(j_var)%dims==3) then
          if (trim(allvars(j_var)%vname)=='phalf') then
             call MDF_Def_Dim( allvars(j_var)%fileunit, 'lev',        lmr+1, lev_dimid, status )
             IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          else 
             call MDF_Def_Dim( allvars(j_var)%fileunit, 'lev',        lmr, lev_dimid, status )
             IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          end if
       end if

       if (allvars(j_var)%table_id/='AERfx')then
          !call MDF_Def_Dim( allvars(j_var)%fileunit, 'time',       nout_steps, time_dimid, status )
          call MDF_Def_Dim( allvars(j_var)%fileunit, 'time',       MDF_UNLIMITED, time_dimid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          call MDF_Def_Dim( allvars(j_var)%fileunit, 'bounds',       2, bounds_dimid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       end if
       if (allvars(j_var)%dims==3) then
          call MDF_Def_Dim( allvars(j_var)%fileunit, 'nhym',       nhym, hym_dimid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          call MDF_Def_Dim( allvars(j_var)%fileunit, 'nhyi',       nhyi, hyi_dimid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       end if

       ! define dimension variables 
       ! dim 0= global sum
       if (allvars(j_var)%dims>0) then
          ! longitude
          call MDF_Def_Var( allvars(j_var)%fileunit, 'lon', MDF_DOUBLE, &
               (/ lon_dimid/), allvars(j_var)%lon_varid, status )
          
          
          ! define latitude variable
          call MDF_Def_Var( allvars(j_var)%fileunit, 'lat', MDF_DOUBLE, &
               (/ lat_dimid/), allvars(j_var)%lat_varid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
                 
       end if
       ! level
       if (allvars(j_var)%dims==3) then
          call MDF_Def_Var( allvars(j_var)%fileunit, 'lev', MDF_DOUBLE, &
               (/ lev_dimid/), allvars(j_var)%lev_varid, status )
       end if

       if (allvars(j_var)%table_id/='AERfx')then
    
          call MDF_Def_Var( allvars(j_var)%fileunit, 'time', MDF_DOUBLE, &
            (/ time_dimid/), allvars(j_var)%time_varid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
          call MDF_Def_Var( allvars(j_var)%fileunit, 'time_bounds', MDF_DOUBLE, &
            (/ bounds_dimid,time_dimid/), allvars(j_var)%bounds_varid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
       end if
       if (allvars(j_var)%dims==3) then
          ! define hybm variable
          call MDF_Def_Var( allvars(j_var)%fileunit, 'hybm', MDF_DOUBLE, &
               (/ hym_dimid/), allvars(j_var)%hybm_varid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
          ! define hyam variable
          call MDF_Def_Var( allvars(j_var)%fileunit, 'hyam', MDF_DOUBLE, &
               (/ hym_dimid/), allvars(j_var)%hyam_varid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
          ! define hybi variable
          call MDF_Def_Var( allvars(j_var)%fileunit, 'hybi', MDF_DOUBLE, &
               (/ hyi_dimid/), allvars(j_var)%hybi_varid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
          ! define hyai variable
          call MDF_Def_Var( allvars(j_var)%fileunit, 'hyai', MDF_DOUBLE, &
               (/ hyi_dimid/), allvars(j_var)%hyai_varid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
        end if
        
       ! * Write out dimension variable values
 
       ! Write out hybm 
       if (allvars(j_var)%dims==3) then
          ! midpoints
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%hybm_varid,levi%fb, status)
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          ! Write out hyam
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%hyam_varid,levi%fa, status)
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          ! interfaces 
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%hybi_varid,levi%b, status)
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          ! Write out hyai
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%hyai_varid,levi%a, status)
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       end if

       ! Write out the longitudes
       if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%lon_varid, dimension_data%lon, status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

       !write out the latitude to variable
       if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,  allvars(j_var)%lat_varid, dimension_data%lat, status)
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       
       if (allvars(j_var)%dims==3) then
          if (trim(allvars(j_var)%vname)=='phalf') then
             !N= levi%fb?
             if (lmr+1==35)then
                ! Write out the levels
                if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%lev_varid, (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35/), status)
                IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
             else if (lmr==11)then
                if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%lev_varid, (/1,2,3,4,5,6,7,8,9,10,11/), status)
                IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
             else
                write(gol,*) 'Number of levels not supported for c4mip output: ,',lmr
                IF_ERROR_RETURN(status=1)
             end if

!             end if
          else

             !N= levi%b
             if (lmr==34) then
                if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%lev_varid, (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34/), status)
                IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
             else if(lmr==10) then
                if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%lev_varid, (/1,2,3,4,5,6,7,8,9,10/), status)
                IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
             else

                write(gol,*) 'Number of levels not supported for c4mip output: ,',lmr

                IF_ERROR_RETURN(status=1)
             end if
          end if

       end if
       ! Time will be written during output write -steps

    end if
    call goLabel()
    status = 0
    
  end subroutine write_dimensions
  
  subroutine write_var(status,j_var)
    implicit none
    integer     :: j_var
    integer,intent(out)::status
    integer     :: i1,i2,j1,j2,imr,jmr,lmr
    integer     :: lon_varid,lonid
    integer     :: lat_varid,latid
    integer     :: lev_varid,levid
    integer     :: time_varid,timeid
    

    character(len=*), parameter ::  rname = mname//'/output_c4mip_write_dim'
    call goLabel(rname)

    ! define dimensions
    !call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    imr=dimension_data%nlon
    jmr=dimension_data%nlat
    lmr=dimension_data%nlev
    ! define dimension variables 
    ! longitude

    if (allvars(j_var)%dims==2.and.allvars(j_var)%table_id/='AERfx') then 
       call MDF_Def_Var( allvars(j_var)%fileunit, allvars(j_var)%vname, MDF_DOUBLE, &
            (/allvars(j_var)%lon_varid, allvars(j_var)%lat_varid, allvars(j_var)%time_varid/), allvars(j_var)%varid, status )
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

    else  if (allvars(j_var)%dims==2.and.allvars(j_var)%table_id=='AERfx') then 
       call MDF_Def_Var( allvars(j_var)%fileunit, allvars(j_var)%vname, MDF_DOUBLE, &
            (/allvars(j_var)%lon_varid, allvars(j_var)%lat_varid/), allvars(j_var)%varid, status )
       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
       else
!(allvars(j_var)%dims==3)


       call MDF_Def_Var( allvars(j_var)%fileunit, allvars(j_var)%vname, MDF_DOUBLE, &
            (/allvars(j_var)%lon_varid, allvars(j_var)%lat_varid, allvars(j_var)%lev_varid,allvars(j_var)%time_varid/), allvars(j_var)%varid, status )

       IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    end if
    ! Write out the longitudes

   
   call MDF_Put_Att( allvars(j_var)%fileunit, allvars(j_var)%varid, 'long_name',     trim(allvars(j_var)%lname), status )
   IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
   call MDF_Put_Att(allvars(j_var)%fileunit,allvars(j_var)%varid, 'standard_name', trim(allvars(j_var)%standard_name), status )
   IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
   call MDF_Put_Att(allvars(j_var)%fileunit , allvars(j_var)%varid, 'units',         trim(allvars(j_var)%unit), status )
   IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

   call MDF_EndDef( allvars(j_var)%fileunit, status )
   IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

   call goLabel()
    status = 0

  end subroutine write_var


end MODULE USER_OUTPUT_C4MIP