ecearth3/sources/tm5mp/proj/output/user_output_column.F90

!#################################################################
!
! MODULE: user_output_column *************** DUMMY TM5-MP VERSION *******************
!
! PUBLIC SUBROUTINES: 
!   user_output_column_init
!   user_output_column_accum
!   user_output_column_evaluate
!   user_output_column_write
!   user_output_column_reset
!   user_output_column_done
!
! DESCRIPTION: 
!   Write column tracer fields. Called from user_output.F90. 
!
! REVISION HISTORY: 
!   Original module written by:
!     Peter Bergamaschi, Maarten Krol, Wouter Peters, others?
!
!   Mike Trudeau, Feb 2012
!     Modified to write NetCDF4 output in NOAA CarbonTracker
!     release format. 
!
!     Required rc file keys
!       output.column : [T/F]
!       output.column.dtsec : [integer] 
!       output.column.infile   : [/path/myColumnList.txt]
!     Obsolete rc file keys
!       output.column.filename  : [column.nc]
!       inputdir.column : [/path]
!       output.column.infilename : [myColumnList.txt]
!
!   Mike Trudeau, Feb 2013
!     1. Modified pressure units (hectopascals to pascals)
!     2. Removed mean pressure, wind speed, wind direction from
!        the meteo output fields.
!     3. CO2 tracers written with background subtracted.
!     4. Total CO2 tracer added to output.
!
!  Andy Jacobson, April 2013
!     Change to ndyn weighting
!
!### macro's #####################################################
!
#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 user_output_column

  use GO,         only : gol, goErr, goPr
!   use dims,       only : nregions, itaur

!   use dims,       only : im, jm, lm, dx, dy, xref, yref, xbeg, ybeg, xend, yend, zbeg, zend
!   use chem_param, only : ntracet, ntrace, names

  implicit none

  character(len=*), parameter   ::  mname = 'user_output_column'

  private

  ! --- public routines ------------------------

  public :: user_output_column_init
  public :: user_output_column_done
  public :: user_output_column_write
  public :: user_output_column_accum
  public :: user_output_column_reset
  public :: user_output_column_evaluate

!   ! --- constants ------------------------------

!   integer, parameter           :: nmethods = 3 ! number of column interpolation methods
!   integer, parameter           :: nmeteoLvl = 4 
!   integer, parameter           :: nmeteoBnd = 2
!   integer, parameter           :: writeMethod = 2 ! interpolation method output to file
!                                                   ! 1: no interpolation ('grd')
!                                                   ! 2: interpolation based on slopes ('slp')
!                                                   ! 3: 3D linear interpolations ('int')

!   character(len=70), dimension(5), parameter :: comment_species = & 
!     (/'background - component due to initial condition (at 2000-01-01)', &
!       'component due to fossil fuel burning (since 2000-01-01)', &
!       'component due to air-sea gas exchange (since 2000-01-01)', &
!       'component due to terrestrial biosphere exchange (since 2000-01-01)', &
!       'component due to direct fire emissions (since 2000-01-01)'/)

!   character(len=14), dimension(nmeteoLvl), parameter :: name_meteoLvl = &
!     (/'temperature', &
!       'u', &
!       'v', &
!       'blh'/)
!   character(len=13), dimension(nmeteoLvl), parameter :: units_meteoLvl = &
!     (/'Kelvin', &
!       'meters/second', &
!       'meters/second', &
!       'meters'/)
!   character(len=50), dimension(nmeteoLvl), parameter :: standard_name_meteoLvl = &
!     (/'temperature at level center', &
!       'zonal component of the wind velocity', &
!       'meridional component of the wind velocity', &
!       'boundary layer height'/)
!   character(len=14), dimension(nmeteoLvl), parameter :: long_name_meteoLvl = &
!     (/'temperature', &
!       'u', &
!       'v', &
!       'blh'/)

!   character(len=14), dimension(nmeteoBnd), parameter :: name_meteoBnd = &
!     (/'pressure', 'gph'/)
!   character(len=13), dimension(nmeteoBnd), parameter :: units_meteoBnd = &
!     (/'pascals', 'meters'/)
!   character(len=50), dimension(nmeteoBnd), parameter :: standard_name_meteoBnd = &
!     (/'pressure at layer interface', 'geopotential height of layer interface'/)
!   character(len=14), dimension(nmeteoBnd), parameter :: long_name_meteoBnd = &
!     (/'pressure', 'gph'/)

!   ! --- variables ------------------------------
  
!   integer, dimension(ntrace)    :: species
!   integer                       :: nspecies
!   integer, save                 :: nstations 
!   logical                       :: evaluated

!   type fieldsType1 
!     character(len=13)                 :: id        ! column identifier
!     character(len=50)                 :: desc      ! column name
!     real                              :: lat       ! latitude of column
!     real                              :: lon       ! longitude of column
!     real                              :: height    ! height of column
!     real                              :: oro       ! ground height at station 
!     integer                           :: index     ! index of column
!     integer                           :: region 
!     integer                           :: is
!     integer                           :: js
!     integer                           :: ls
!     real,dimension(nmeteoLvl,lm(1))   :: meteoLvl  ! meteo (interpolated) at columns
!     real,dimension(:,:,:),pointer     :: c_avg     ! species, interpolation
!     real,dimension(:,:,:),pointer     :: c_std
!     integer                           :: ndata    ! for averages
!     real	                      :: accum_weight    ! for averages
!     real,dimension(nmeteoBnd,lm(1)+1) :: meteoBnd  ! meteo (interpolated) at columns
!   end type fieldsType1 

!   type(fieldsType1), dimension(:), allocatable :: stations
!   type(fieldsType1), dimension(:), allocatable :: stationsRegion

!   type fieldsType2 
!     integer, dimension(6)                     :: idate_start
!     integer, dimension(6)                     :: idate_end
!   end type fieldsType2 

!   type(fieldsType2), dimension(nregions), save     :: timestamp 

!   ! --- begin -------------------------------

  contains

  subroutine user_output_column_init( status )
  
!     use ParTools ,   only : myid, root
!     use dims,        only : nregions, idatei

!     implicit none

!     ! --- in/out ------------------------
 
!     integer                    :: region 
    integer, intent(out)       :: status
 
    ! --- constants ---------------------
    
!     character(len=*), parameter  ::  rname = mname//'/user_output_column_init'
    
!     ! --- begin -------------------------------

!     if (myid /= root) return

!     evaluated = .FALSE.

!     ! grab the date
!     do region = 1, nregions
!       timestamp(region)%idate_start = idatei
!     end do 

!     ! read user supplied station list 
!     call read_columnlist( status )
!     IF_NOTOK_RETURN( status=1 )

!     ! initialize the column output
!     call user_output_column_reset( status )

    ! ok
    status = 0

  end subroutine user_output_column_init

!   subroutine read_columnlist( status )
!     use toolbox,      only : escape_tm
! #ifdef MPI  
!     use mpi_const, only: myid, root, mpi_integer, mpi_character, &
!          my_real, mpi_comm_world, ierr
! #else
!     use ParTools, only: myid, root
! #endif
!     use GO,          only : TrcFile, Init, Done, ReadRc
!     use global_data, only : rcfile
    

!     implicit none

!     ! --- variables ---------------------
    
!     integer, intent(out) :: status
!     character(len=1024)               :: infile
!     integer                           :: nchar !string length
!     character(len=200)                :: dummy
!     integer                           :: q
!     integer                           :: region
!     integer,dimension(6)              :: idate_sim
!     integer                           :: itau_sim, isim, ispec
!     real                              :: dxr, dyr 
!     character(len=8)                  :: name_spec    !CMK changed from 6-->8
!     integer                           :: sunit0=110
!     type(TrcFile)                     :: rcF

!     ! --- constants -----------------------

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

!     ! --- begin -----------------------

!     ! start
!     if ( myid == root ) then
!        ! read rc file
!        call Init( rcF, rcfile, status)

!        call ReadRc( rcF, 'output.column.infile', infile, status )
!        IF_NOTOK_RETURN( status=1 )

!        write (gol, *) trim(rname)//'/column file: ', trim(infile); call goPr

!        call Done( rcF, status)

!        ! read columnlist file
!        nchar = len_trim(infile) 
!        if ( nchar > 0 ) then
!           !=========================================
!           ! open columnFile (with list of columns)
!           !=========================================
!           open(UNIT=sunit0, FORM='FORMATTED', &
!                FILE=trim(infile), &
!                status='OLD', ERR=1000)      
!           read(sunit0, '(a)') dummy

!           ! count number of columns
!           nstations=0
!           stations_loop: DO
!              read(sunit0,  '(a)', END=100) dummy
!              if (dummy(1:7) /=  'species') then   
!                 nstations = nstations+1 
!              else
!                 nspecies = 0
!                 do
!                    read(sunit0, '(a8)' , end=100) name_spec    ! CMK changed a6-->a8
!                    !CMK not allowed yet to write short-lived: TODO
!                    specloop: do ispec = 1, ntracet
!                       if (name_spec == names(ispec) ) then
!                          nspecies = nspecies + 1
!                          species(nspecies) = ispec
!                          exit specloop
!                       endif
!                    enddo specloop
!                 enddo
!                 exit stations_loop
!              end if
!           end do stations_loop
! 100       continue

!           print *, '_________________________________________________' &
!                //'_______________________________'
!           print *, nstations, ' columns read from file: ', infile
!           print *, nspecies, ' concentrations will be gathered:'
!           if (nspecies == 0) then 
!              print *, 'ERROR user_output_column'
!              print *, 'ERROR no species for output'
!              print *, 'ERROR Set column_output to .F. or give valid species name'
!              call escape_tm('Error reading column file')
!           endif

!           do ispec = 1, nspecies
!              print '(2i4,2x,a8)', ispec, species(ispec), &   ! CMK change a6--->a8
!                   names(species(ispec))
!           enddo
!           print *, '_________________________________________________' &
!                //'_______________________________'


!           allocate(stations(nstations))

!           do q = 1, nstations
!              ! wow, this lm(1) assumes column in region 1, or that 
!              ! child regions all have lm(n) <= lm(1)
!              allocate(stations(q)%c_avg(nspecies,nmethods,lm(1)))
!              allocate(stations(q)%c_std(nspecies,nmethods,lm(1)))
!           enddo
!           dyr = dy/yref(1) ! again, as above comment: region 1 assumed.
!           dxr = dx/xref(1) !  "  "  "

!           !===============
!           ! read columns
!           !===============

!           rewind(sunit0)
!           read(sunit0, '(a)') dummy
!           do q=1, nstations 
!              read(sunit0,  '(a7,1x,f8.2,f8.2,f8.1,1x,a50)') &
!                   stations(q)%id, &
!                   stations(q)%lat, &
!                   stations(q)%lon, &
!                   stations(q)%height, &
!                   stations(q)%desc

!              stations(q)%oro = -99999.0

!              stations(q)%index = q

!              ! assume global region as default for average mixing ratios
!              stations(q)%region = 1
!              stations(q)%is = &
!                   int((stations(q)%lon-float(xbeg(1)))/dxr + 0.99999) 
!              stations(q)%js = &
!                   int((stations(q)%lat-float(ybeg(1)))/dyr + 0.99999)
!           end do

!           ! close columnfile
!           close(sunit0)
!           do q=1, nstations
!             do region=1, nregions
      
!               dyr = dy/yref(region)
!               dxr = dx/xref(region)

!               if ( (stations(q)%lon .ge. xbeg(region) .and. &
!                     stations(q)%lon .le. xend(region)) .and.      &
!                    (stations(q)%lat .ge. ybeg(region) .and. &
!                     stations(q)%lat .le. yend(region) ) ) then
!                    !=====================
!                    ! column is in region
!                    !=====================

!                 ! determine region with hightes refinement factor
!                 if (xref(region) > xref(stations(q)%region)) then
!                   stations(q)%region = region
!                   stations(q)%is = &
!                     int((stations(q)%lon-float(xbeg(region)))/dxr + 0.99999) 
!                   stations(q)%js = &
!                     int((stations(q)%lat-float(ybeg(region)))/dyr + 0.99999)
!                 end if
!               end if
!             end do
!           end do
!           do q=1, nstations 
!              print '(a7,1x,f8.2,f8.2,f8.1,1x,i4,1x,a50)', &
!                   stations(q)%id, &
!                   stations(q)%lat, &
!                   stations(q)%lon, &
!                   stations(q)%height, &
!                   stations(q)%region, &
!                   stations(q)%desc
!           end do
!        end if
!     end if  ! myid = root  

!     ! ok
!     status = 0

!     return

!     ! error handling
! 1000 call escape_tm( 'read_columnlist: could not open '//trim(infile))

!   end subroutine read_columnlist

  subroutine user_output_column_accum( region, status )

! #ifdef MPI  
!     use mpi_const !, only: myid, root
!     use mpi_comm, only: gather_tracer_k, gather_tracer_t, gather_tracer_k_3d
! #endif
!     use ParTools, only: myid, root, root_t, root_k
!     use global_data, only : mass_dat,  conv_dat
!     use datetime,    only : tau2date
!     use toolbox,     only : escape_tm
!     use Meteo,       only : gph_dat, oro_dat, temper_dat, pu_dat, pv_dat, m_dat, phlb_dat
!     use dims,        only : xcyc,ndyn,ndyn_max
!     use chem_param,   only : fscale, uscale
!     use binas,      only : ae, twopi

!     implicit none

    ! --- in/out ------------------------

    integer, intent(out)             :: status
  
    ! --- variables -----------------------

    integer, intent(in)              :: region

!     real,dimension(:,:,:,:),pointer :: rm                   ! tracer mass
!     real,dimension(:,:,:,:),pointer :: rxm,rym,rzm          ! slopes of tracer mass
!     real,dimension(:,:,:)  ,pointer :: m                    ! air mass
!     real,dimension(:,:,:)  ,pointer :: gph                  ! geopotential height
!     real,dimension(:,:)    ,pointer :: oro                  ! orography

!     real, dimension(:,:,:,:), allocatable, target      :: rmg                      ! MPI arrays to gather fields.
!     real, dimension(:,:,:,:), allocatable, target      :: rxmg, rymg, rzmg 
!     real, dimension(:,:,:  ), allocatable, target      :: mg

!     integer                         :: i_interpol           ! interpolation procedures:   
!     ! 1: no interpolation
!     ! 2: interpolation based on slopes
!     ! 3: 3D linear interpolations


!     integer                         :: i,j,l                ! grid cell indices
!     integer                         :: is,js,ls             ! i,j,l index of grid cell in which column lis located
!     integer                         :: isn,jsn,lsn          ! i,j,l index of grid cell which is taken as neighbour for linear interpolation
!     integer                         :: lst, lstn            ! MPI implementation
!     integer                         :: js_north, js_south   ! j index of grid cell for neighbours for windspeed_v interpolation
!     integer                         :: is_west              ! i index of grid cell for neighbour for windspeed_u interpolation
!     integer                         :: n                    ! index of tracer
!     integer                         :: q               ! index of column
!     real                            :: ris,rjs,rls          ! deviation of column from center of the grid cell (-0.5 ... +0.5)
!     real                            :: wcx,wcy,wcz          ! x,y,z-weight of grid cell (1.0 ... 0.5) 
!     real                            :: dxr, dyr             ! x,y resolution (in degrees) for current region 

!     real,dimension(0:lm(region))    :: height               ! height of vertical grid boundaries 
!     integer, dimension(6)           :: idate_f              ! current idate for region 
!     integer                         :: sunit                ! unit number for column output file
!     real,dimension(lm(1))           :: rm_interpol          ! tracer mass, interpolated to column location 
!     integer                         :: ispec                ! index over activated tracers

!     !=============
!     ! meteoLvl output
!     !=============
!     real,dimension(nstations,nmeteoLvl,lm(1)) :: meteoLvl   ! meteo (interpolated) at columns
!     real,dimension(:,:,:), pointer  :: T                    ! temperature
!     real,dimension(:,:,:), pointer  :: phlb                 ! pressure grid boundaries
!     real,dimension(:,:,:), pointer  :: pu                   ! mass flux x-direction [kg/s]
!     real,dimension(:,:,:), pointer  :: pv                   ! mass flux y-direction [kg/s]
!     real,dimension(:,:), pointer    :: blh                  ! boundary layer height [m]
!     real,dimension(lm(1))           :: T_interpol           ! temperature, interpolated to column location 
!     real,dimension(lm(1))           :: p_interpol           ! pressure, interpolated to column location 
!     real,dimension(lm(1)+1)         :: p_interpol_bnd       ! boundary pressure, interpolated to column location 
!     real,dimension(lm(1))           :: windspeed_u_interpol ! wind speed [m/s] x-direction , interpolated to column location
!     real,dimension(lm(1))           :: windspeed_v_interpol ! wind speed [m/s] y-direction , interpolated to column location
!     real,dimension(lm(1))           :: windspeed_interpol   ! wind speed [m/s]             , interpolated to column location
!     real,dimension(lm(1))           :: winddir_interpol     ! wind direction               , interpolated to column location 
!     real                            :: blh_interpol         ! boundary layer height        , interpolated to column location 
!     real                            :: lat_j                ! latitude of grid cell center js  [degrees]
!     real                            :: lat_jn               ! latitude of grid cell center jsn [degrees]
!     real                            :: ddx_j                ! x-extension of grid cell js    [m]
!     real                            :: ddx_jn               ! x-extension of grid cell jsn   [m]
!     real                            :: ddy                  ! y-extension of grid cell js    [m]
!     real,dimension(lm(1))           :: vmr                  ! volume mi
!     real			    :: this_weight
!     integer                         :: communicator,root_id,lmr,imr,jmr

!     ! start

!     ! assign pointers depending on the paralel status
!     ! (over levels after chemistry!)
!     imr = im(region) ; jmr = jm(region) ; lmr = lm(region)

!     this_weight=real(ndyn)/real(ndyn_max)

! #ifdef MPI

!     if ( root_k /= root_t ) print *, 'Warning from MPI column output: ', &
!          'root_k and root_t are not equal'

!     allocate( rmg (-1:imr+2,-1:jmr+2,lmr, ntracet) ) ; rmg  = 0.0
!     allocate( rxmg(-1:imr+2,-1:jmr+2,lmr, ntracet) ) ; rxmg = 0.0
!     allocate( rymg(-1:imr+2,-1:jmr+2,lmr, ntracet) ) ; rymg = 0.0
!     allocate( rzmg(-1:imr+2,-1:jmr+2,lmr, ntracet) ) ; rzmg = 0.0
!     allocate( mg  (-1:imr+2,-1:jmr+2,lmr         ) ) ; mg   = 0.0
!     if ( which_par == 'tracer' ) then 
!        m => m_dat(region)%data 
!        rm => mass_dat(region)%rm_t
!        rxm => mass_dat(region)%rxm_t
!        rym => mass_dat(region)%rym_t
!        rzm => mass_dat(region)%rzm_t
!        call gather_tracer_t(rmg,imr,jmr,lmr,2,2,0,ntracet,rm,.false.)
!        call gather_tracer_t(rxmg,imr,jmr,lmr,2,2,0,ntracet,rxm,.false.)
!        call gather_tracer_t(rymg,imr,jmr,lmr,2,2,0,ntracet,rym,.false.)
!        call gather_tracer_t(rzmg,imr,jmr,lmr,2,2,0,ntracet,rzm,.false.)
!        nullify(rm)   ! reassign pointers
!        nullify(rxm)
!        nullify(rym)
!        nullify(rzm)
!        rm => rmg
!        rxm => rxmg
!        rym => rymg
!        rzm => rzmg
!        communicator=com_trac  !WP! assign com_trac as communicator
!        root_id=root_t         !root_t is the PE that knows the result.
!     else
!        m => m_dat(region)%data 
!        rm => mass_dat(region)%rm_k
!        rxm => mass_dat(region)%rxm_k
!        rym => mass_dat(region)%rym_k
!        rzm => mass_dat(region)%rzm_k
!        call gather_tracer_k(rmg,imr,jmr,lmr,2,2,0,ntracet,rm,.false.)
!        call gather_tracer_k(rxmg,imr,jmr,lmr,2,2,0,ntracet,rxm,.false.)
!        call gather_tracer_k(rymg,imr,jmr,lmr,2,2,0,ntracet,rym,.false.)
!        call gather_tracer_k(rzmg,imr,jmr,lmr,2,2,0,ntracet,rzm,.false.)
!        call gather_tracer_k_3d( mg, imr,jmr,lmr, 2,2,0, m, .false.)
       
!        nullify(m)   ! reassign pointers
!        nullify(rm) 
!        nullify(rxm)
!        nullify(rym)
!        nullify(rzm)
!        m => mg
!        rm => rmg
!        rxm => rxmg
!        rym => rymg
!        rzm => rzmg
!        communicator=com_lev  !WP! assign com_lev as communicator
!        root_id=root_k
!     end if
! #else
!     m => m_dat(region)%data 
!     rm => mass_dat(region)%rm_t
!     rxm => mass_dat(region)%rxm_t
!     rym => mass_dat(region)%rym_t
!     rzm => mass_dat(region)%rzm_t
! #endif
!     ! assign pointers (METEO)
!     gph => gph_dat(region)%data
!     oro => oro_dat(region)%data(:,:,1)
!     t => temper_dat(region)%data
!     pu => pu_dat(region)%data
!     pv => pv_dat(region)%data
!     blh => conv_dat(region)%blh
!     ! point to pressure paralel over tracers (contains all layers!)
!     phlb => phlb_dat(region)%data

!     if ( myid == root_t ) then 
!        ! perform only on root: tracers are already gathered here
!        ! x and y resolution [degrees] for current region

!        dyr = dy/yref(region)
!        dxr = dx/xref(region)

!        ! idate for current region
!        call tau2date(itaur(region),idate_f)


!        !====================
!        ! loop over columns
!        !====================

!        do q=1, nstations

!           IF (stations(q)%region == region)THEN

!              !avoid edges!
!              ris = (stations(q)%lon-float(xbeg(region)))/dxr + 0.99999
!              rjs = (stations(q)%lat-float(ybeg(region)))/dyr + 0.99999
!              is  = int(ris)   ! i-index of grid cell in which column is located
!              js  = int(rjs)   ! j-index of grid cell in which column is located
!              ! x-deviation from the center of the grid cell  (-0.5...+0.5)
!              ris = ris-is-0.5
!              ! y-deviation from the center of the grid cell  (-0.5...+0.5)
!              rjs = rjs-js-0.5

!              !=================================
!              !the neighbour for x interpolation
!              !=================================
!              if ( ris .gt. 0 ) then 
!                 isn = is+1     
!              else
!                 isn = is-1
!              endif

!              !=================================
!              !the neighbour for y interpolation
!              !=================================
!              if ( rjs .gt. 0 ) then 
!                 jsn = js+1     
!              else
!                 jsn = js-1
!              endif

!              ! x- / y-weighting of grid cell in which column is located
!              wcx = (1.0-abs(ris))    ! 1.0 ... 0.5
!              wcy = (1.0-abs(rjs))    ! 1.0 ... 0.5


!              !=================================================================
!              ! if index of neighbour is exceeding range of region set 
!              ! neighbour = current cell (i.e. no interpolation)
!              ! in case of cyclic x-boundaries take corresponding cyclic i index
!              !=================================================================
             
!              if ( jsn == 0) jsn=1
!              if ( jsn == jm(region)+1 ) jsn=jm(region)  ! isn-->jsn (wouter Peters)
!              if ( xcyc(region) == 0 ) then
!                 ! non-cyclic boundaries
!                 if ( isn == 0) isn=1
!                 if ( isn == im(region)+1 ) isn=im(region)
!              else
!                 ! cyclic x-boundaries
!                 if ( isn == 0 ) isn=im(region)
!                 if ( isn == im(region)+1 ) isn=1
!              end if

!              !============================================
!              ! interpolate the gph to xy position of column
!              !============================================
!              !ls = 1   !layer
!              do l=0,lm(region)
!                 !CMK note: gph now from 1-->lm+1 (dec 2002)
!                 height(l) =     wcx *      wcy*  gph(is,js,l+1)  + &
!                            (1.0-wcx)*      wcy*  gph(isn,js,l+1) + &
!                                 wcx *(1.0-wcy)*  gph(is,jsn,l+1) + &
!                            (1.0-wcx)*(1.0-wcy)*  gph(isn,jsn,l+1)
!              end do

!              !===========================
!              ! determine layer of column
!              !===========================
!              !do l=0,lm(region)
!              !   if(height(l) .gt. stations(q)%height) exit
!              !end do

!              !select case(l)
!              !case(0)
!              !   ! force column to model surface
!              !   ls = 1
!              !   rls = -0.5  
!              !case default
!              !   ls = l       ! column's layer
!              !   ! the off-set from the center of the layer (-0.5--->+0.5) 
!              !   ! (interpolation is in (m))
!              !   rls = (stations(q)%height - height(l-1)) / &
!              !        (height(l)-height(l-1)) - 0.5
!              !end select
!              !stations(q)%ls = ls  ! store for output to file

!              rls = 0.0 ! no interpolation
!              wcz = 1.0 ! no interpolation

!              !============================================
!              ! if not yet done, interpolate the orography
!              ! to xy position of column
!              !============================================
!              if (stations(q)%oro < -10000) then
!                  stations(q)%oro =     wcx *      wcy*  oro(is,js)  + &
!                            (1.0-wcx)*      wcy*  oro(isn,js) + &
!                                 wcx *(1.0-wcy)*  oro(is,jsn) + &
!                            (1.0-wcx)*(1.0-wcy)*  oro(isn,jsn)
!              endif

             
!              do ispec=1,nspecies
!                 n = species(ispec)   
!                 lst = ls 
!                 lstn = lsn 
!                 !========================================
!                 ! value of grid box without interpolation
!                 !========================================
!                 vmr = rm(is,js,:,n) / m(is,js,:) *  fscale(n) * uscale(n)
!                 stations(q)%c_avg(ispec,1,:) = &
!                    stations(q)%c_avg(ispec,1,:) + this_weight*vmr
!                 stations(q)%c_std(ispec,1,:) = &
!                    stations(q)%c_std(ispec,1,:) + this_weight*vmr*vmr
!                 !========================================
!                 ! interpolation using slopes  
!                 !========================================
!                 !rm-value is obtained from rm + slopes. 
!                 !slope = rxm = (rm*dX/dx *deltaX/2)
!                 rm_interpol = rm(is,js,:,n) + 2.0*(ris*rxm(is,js,:,n) + &
!                      rjs*rym(is,js,:,n) + &
!                      rls*rzm(is,js,:,n) )
!                 rm_interpol = max(0.0,rm_interpol) 
!                 vmr = rm_interpol/m(is,js,:) * fscale(n) * uscale(n)
!                 stations(q)%c_avg(ispec,2,:) = &
!                    stations(q)%c_avg(ispec,2,:) + this_weight*vmr
!                 stations(q)%c_std(ispec,2,:) = &
!                    stations(q)%c_std(ispec,2,:) + this_weight*vmr*vmr

!                 !========================================
!                 ! 3D linear interpolation 
!                 !========================================
!                 ! this PE handles also the neighbour layer
! !PB             if ( lstn >= 1 .and. lstn <= lmr) then
! !PBi
!                 rm_interpol = wcx      *wcy          *rm(is,js,:,n)   /m(is,js,:)   + & 
!                                  (1.0-wcx)*wcy          *rm(isn,js,:,n)  /m(isn,js,:)  + &
!                                  wcx      *(1.0-wcy)    *rm(is,jsn,:,n)  /m(is,jsn,:)  + &
!                                  (1.0-wcx)*(1.0-wcy)    *rm(isn,jsn,:,n) /m(isn,jsn,:) 

!                 vmr = rm_interpol * fscale(n) * uscale(n)
!                 stations(q)%c_avg(ispec,3,:) = &
!                    stations(q)%c_avg(ispec,3,:) + this_weight*vmr
!                 stations(q)%c_std(ispec,3,:) = &
!                    stations(q)%c_std(ispec,3,:) + this_weight*vmr*vmr


!              end do  !ispec



!                 !======================
!                 !meteoLvl data at columns
!                 !======================
!                 !if (output_columns_meteoLvl) then              

!                    !============================================
!                    ! gph already interpolated to xy position of column
!                    ! Note that height is defined with index 0:lm(region)
!                    ! whereas meteoLvl(8,:) is 1:lm(region)
!                    !============================================
!                    stations(q)%meteoBnd(2,:) = &
!                         stations(q)%meteoBnd(2,:) + this_weight*height(0:(lm(region)))

!                    !============================================
!                    ! interpolate the blh to xy position of column
!                    !============================================
!                    blh_interpol =     wcx *      wcy*  blh(is,js)  + &
!                                      (1.0-wcx)*      wcy*  blh(isn,js) + &
!                                      wcx *(1.0-wcy)*  blh(is,jsn) + &
!                                      (1.0-wcx)*(1.0-wcy)*  blh(isn,jsn)
!                    stations(q)%meteoLvl(4,:) = &
!                       stations(q)%meteoLvl(4,:) + this_weight*blh_interpol
!                    !====================================
!                    ! 3D linear interpolation Temperature
!                    !====================================
!                    T_interpol = &
!                         wcx      *wcy      *wcz       *T(is,js,:)    + & 
!                         (1.0-wcx)*wcy      *wcz       *T(isn,js,:)   + &
!                         wcx      *(1.0-wcy)*wcz       *T(is,jsn,:)   + &
!                         (1.0-wcx)*(1.0-wcy)*wcz       *T(isn,jsn,:) 

!                    stations(q)%meteoLvl(1,:) = &
!                       stations(q)%meteoLvl(1,:) + this_weight*T_interpol
!                    !====================================
!                    ! 3D linear interpolation pressure
!                    !====================================
!                    p_interpol =(((0.5-rls) * phlb(is,js,1:lm(1))  + (0.5+rls) * phlb(is,js,2:lm(1)+1))  * wcx         * wcy       + &
!                         ((0.5-rls) * phlb(isn,js,1:lm(1)) + (0.5+rls) * phlb(isn,js,2:lm(1)+1)) * (1.0-wcx)   * wcy       + &
!                         ((0.5-rls) * phlb(is,jsn,1:lm(1)) + (0.5+rls) * phlb(is,jsn,2:lm(1)+1)) * wcx         * (1.0-wcy) + &
!                         ((0.5-rls) * phlb(isn,jsn,1:lm(1))+ (0.5+rls) * phlb(isn,jsn,2:lm(1)+1))* (1.0-wcx)   * (1.0-wcy) ) ![Pa]

!                    !stations(q)%meteoLvl(2,:) = &
!                    !   stations(q)%meteoLvl(2,:) + p_interpol

!                    p_interpol_bnd =(phlb(is,js,1:lm(1)+1) * wcx * wcy + &
!                                     phlb(isn,js,1:lm(1)+1) * (1.0-wcx) * wcy + &
!                                     phlb(is,jsn,1:lm(1)+1) * wcx * (1.0-wcy) + &
!                                     phlb(isn,jsn,1:lm(1)+1) * (1.0-wcx) * (1.0-wcy)) ![Pa]

!                    stations(q)%meteoBnd(1,:) = &
!                       stations(q)%meteoBnd(1,:) + this_weight*p_interpol_bnd 

!                    !====================================
!                    ! 3D windspeed_u_interpol (x-direction)
!                    !====================================

!                    ! pu is eastward flux through east grid box boundary
!                    ! windspeed_u_interpol wind component from east     

!                    ! latitude of grid cell center js  [degrees]
!                    lat_j =ybeg(region)+(js -0.5)*dyr
!                    ! latitude of grid cell center jsn [degrees]
!                    lat_jn=ybeg(region)+(jsn-0.5)*dyr
!                    ! x-extension of grid cell js   [m]
!                    ddx_j =ae * twopi * dxr / 360.0 * cos(lat_j*twopi/360.0)
!                    ! x-extension of grid cell jsn  [m]
!                    ddx_jn=ae * twopi * dxr / 360.0 * cos(lat_jn*twopi/360.0)

!                    is_west=is-1
!                    if ( xcyc(region) == 0 ) then
!                       ! non-cyclic boundaries
!                       if ( is_west == 0 ) is_west=1
!                    else
!                       ! cyclic x-boundaries
!                       if ( is_west == 0 ) is_west=im(region)
!                    end if

!                    !west border                         !east border  
!                    windspeed_u_interpol=-(((0.5-ris) * pu(is_west,js,:)/m(is_west,js,:)                                   + &
!                         (0.5+ris) * pu(is,js,:)/ m(is,js,:) )      * &
!                         ddx_j  *      wcy *      wcz          + &  
!                         ((0.5-ris) * pu(is_west,jsn,:)/m(is_west,jsn,:)                                 + & 
!                         (0.5+ris) * pu(is,jsn,:)/m(is,jsn,:) )     * &
!                         ddx_jn * (1.0-wcy)*      wcz          + &  
!                         ((0.5-ris) * pu(is_west,js,:)/m(is_west,js,:)                                 + &
!                         (0.5+ris) * pu(is,js,:)/ m(is,js,:))    * &
!                         ddx_j  *      wcy * (1.0-wcz)         + &  
!                         ((0.5-ris) * pu(is_west,jsn,:)/m(is_west,js,:)                                + &
!                         (0.5+ris) * pu(is,jsn,:)/ m(is,jsn,:) ) * &
!                         ddx_jn * (1.0-wcy)* (1.0-wcz)         )
!                    !====================================
!                    ! 3D windspeed_v_interpol (y-direction)
!                    !====================================

!                    ! pv northward flux through south grid box boundary 
!                    ! windspeed_v_interpol wind component from north    

!                    ddy =ae * twopi * dyr / 360.0 ! y-extension of grid cell [m]

!                    js_south = js-1
!                    if (js_south==0) js_south=1

!                    js_north = js+1
!                    if (js_north==(jm(region)+1)) js_north=jm(region)

!                    !south border                  !north border  
!                    windspeed_v_interpol=-(( &
!                         (0.5-rjs)*pv(is,js,:)/m(is,js_south,:)       + &
!                         (0.5+rjs)*pv(is,js_north,:)/ m(is,js,:) )     * &
!                         wcx *      wcz           + &  
!                         ((0.5-rjs)*pv(isn,js,:)/m(isn,js_south,:)    + &
!                         (0.5+rjs)*pv(isn,js_north,:)/ m(isn,js,:) )   * &
!                         (1.0 -  wcx)*      wcz           + &     
!                         ((0.5-rjs)*pv(is,js,:)/m(is,js_south,:)    + &
!                         (0.5+rjs)*pv(is,js_north,:)/ m(is,js,:) )   * &
!                         wcx * (1.0-wcz)          + &  
!                         ((0.5-rjs)*pv(isn,js,:)/m(isn,js_south,:)  + &
!                         (0.5+rjs)*pv(isn,js_north,:)/ m(isn,js,:) ) * &
!                         (1.0 -  wcx)* (1.0-wcz)        ) * &
!                         ddy 


!                    !windspeed_interpol = &
!                    !     sqrt(windspeed_u_interpol**2 + windspeed_v_interpol**2)

!                    stations(q)%meteoLvl(2,:) = &
!                       stations(q)%meteoLvl(2,:) + this_weight*windspeed_u_interpol
!                    stations(q)%meteoLvl(3,:) = &
!                       stations(q)%meteoLvl(3,:) + this_weight*windspeed_v_interpol

!                 !endif  ! meteoLvl out...
!                 stations(q)%ndata = stations(q)%ndata + 1
!                 stations(q)%accum_weight = stations(q)%accum_weight + this_weight
!            end if
!        end do ! end q loop

!     endif  ! myid = root_t

!     nullify(m)
!     nullify(rm)
!     nullify(rxm)
!     nullify(rym)
!     nullify(rzm)
!     nullify(gph)
!     nullify(oro)
!     nullify(t)
!     nullify(pu)
!     nullify(pv)
!     nullify(phlb)
!     nullify(blh)
! #ifdef MPI    
!     deallocate(rmg)
!     deallocate(rxmg)
!     deallocate(rymg)
!     deallocate(rzmg)
!     deallocate(mg)
! #endif

    ! ok
    status = 0

  end subroutine user_output_column_accum

  subroutine user_output_column_evaluate(status)

    !==================================================================
    ! evaluate mean concentration and std
    !==================================================================

!     use dims, only  : idate
!     use ParTools, only: myid, root
!     use toolbox,     only : escape_tm
!     use binas,      only : twopi

!     implicit none
  
    integer, intent(out)       :: status

!     integer q, i, ii,iii
!     integer ndata, isds, nsds, itp
!     real				:: accum_weight
!     real temp
!     real, dimension(:,:), allocatable    :: cout                 ! to get output
!     integer, dimension(:,:), allocatable :: iidate
!     real                                 :: windspeed_u_interpol ! wind speed [m/s] x-direction , interpolated to column location
!     real                                 :: windspeed_v_interpol ! wind speed [m/s] y-direction , interpolated to column location
!     real                                 :: windspeed_interpol   ! wind speed [m/s]             , interpolated to column location
!     real                                 :: winddir_interpol     ! wind direction               , interpolated to column location 
!     if(myid /= root) return

!     if(evaluated) return

!     do q=1, nstations
!        accum_weight = stations(q)%accum_weight
!        if (accum_weight > 0.0) then
!           stations(q)%c_avg(:,:,:) = stations(q)%c_avg(:,:,:) / accum_weight
!           stations(q)%meteoLvl(:,:)   = stations(q)%meteoLvl(:,:)   / accum_weight
!           stations(q)%meteoBnd(:,:)   = stations(q)%meteoBnd(:,:)   / accum_weight
!           if (ndata > 1) then
!              do i=1, nspecies
!                 do ii= 1, nmethods
!                 do iii= 1, lm(1)
!                    temp = (((stations(q)%c_std(i,ii,iii)/accum_weight) -  real(ndata) * ((stations(q)%c_avg(i,ii,iii))**2) ) / (real(ndata)-1.0) )
!                    if(temp > 0) then 
!                       stations(q)%c_std(i,ii,iii) = sqrt(temp)
!                    else
!                       stations(q)%c_std(i,ii,iii) = 0.0
!                    endif
!                 enddo
!                 enddo
!              enddo
!           else
!              stations(q)%c_std(:,:,:)=-1.0
!           end if
!        else
!           stations(q)%c_avg(:,:,:)= -1.0
!           stations(q)%c_std(:,:,:) = -1.0
!           stations(q)%meteoLvl(:,:)  = -999.
!           stations(q)%meteoBnd(:,:)  = -999.
!        end if

!       !do i=1,lm(1)
!       !   windspeed_u_interpol = stations(q)%meteoLvl(3,i)
!       !   windspeed_v_interpol = stations(q)%meteoLvl(4,i)
!       !   windspeed_interpol   = stations(q)%meteoLvl(5,i)

!       !   !====================================
!       !   ! wind direction
!       !   !====================================

!       !   ! definition of winddirection:          
!       !   ! east     :  90
!       !   ! south    : 180
!       !   ! west     : 270
!       !   ! north    : 360
!       !   if ((windspeed_u_interpol > 0) .and. (windspeed_v_interpol > 0))then              ! wind from north east (0...90)
!       !      winddir_interpol = atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi
!       !   else if ((windspeed_u_interpol > 0) .and. (windspeed_v_interpol < 0)) then     ! wind from south east (90...180)
!       !      winddir_interpol = 180.0+atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi
!       !   else if ((windspeed_u_interpol > 0) .and. (windspeed_v_interpol == 0)) then    ! wind from east (90)
!       !      winddir_interpol = 90.0
!       !   else if ((windspeed_u_interpol < 0) .and. (windspeed_v_interpol > 0)) then     ! wind from north west (270... 360)
!       !      winddir_interpol = 360.0+atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi
!       !   else if ((windspeed_u_interpol < 0) .and. (windspeed_v_interpol < 0)) then     ! wind from south west (180...270)
!       !      winddir_interpol = 180.0+atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi
!       !   else if ((windspeed_u_interpol < 0) .and. (windspeed_v_interpol == 0)) then    ! wind from west (270)
!       !      winddir_interpol = 270.0
!       !   else if ((windspeed_u_interpol == 0) .and. (windspeed_v_interpol > 0)) then    ! wind from north (360)
!       !      winddir_interpol = 360.0
!       !   else if ((windspeed_u_interpol == 0) .and. (windspeed_v_interpol < 0)) then    ! wind from south (180)
!       !      winddir_interpol = 180.0
!       !   else if ((windspeed_u_interpol == 0) .and. (windspeed_v_interpol == 0)) then   ! no wind
!       !      winddir_interpol = -1.0
!       !   else
!       !      call escape_tm('output_columnconc: error wind direction')
!       !   end if
!       !   stations(q)%meteoLvl(6,i) = winddir_interpol
!       ! end do
!     end do

!     evaluated = .TRUE.

    ! ok
    status = 0

  end subroutine user_output_column_evaluate

  subroutine user_output_column_write( region, status )

!     use ParTools,    only : myid, root
!     use global_data, only : outdir
!     use dims,        only : region_name
!     use datetime,    only : tau2date, date2tau, idate2ddate
!     use MDF,         only : MDF_Create, MDF_Close, MDF_EndDef
!     use MDF,         only : MDF_NETCDF, MDF_REPLACE, MDF_GLOBAL, MDF_INT, MDF_FLOAT, MDF_UNLIMITED, MDF_DOUBLE, MDF_CHAR
!     use MDF,         only : MDF_Put_Att, MDF_Def_Dim, MDF_Def_Var, MDF_Put_Var

!     implicit none
  
    ! --- in/out ------------------------

    integer, intent(out)             :: status
    integer, intent(in)              :: region
  
    ! --- variables -----------------------

!     character(len=1024)                    :: outfile 
!     integer                                :: fid
!     integer                                :: i, ii, q
!     integer                                :: dimid_stations, dimid_date, dimid_cal, dimid_nchar
!     integer                                :: dimid_lvl, dimid_bnd
!     integer                                :: varid_desc, varid_id, varid_date_int, varid_column_int
!     integer                                :: varid_elapsed, varid_dec, varid_cal_int, varid_oro
!     integer                                :: varid_lvl, varid_bnd, varid_lon, varid_lat, varid_hgt, varid_total
!     integer, dimension(nspecies, nmethods) :: varid_species
!     integer, dimension(nmeteoLvl)          :: varid_meteoLvl
!     integer, dimension(nmeteoBnd)          :: varid_meteoBnd
!     integer                                :: itau_start
!     integer                                :: itau_end
!     integer                                :: itau_avg
!     integer                                :: itau_ref
!     integer                                :: itau_elapsed
!     real*8                                 :: ddate_avg
!     real*8                                 :: days_elapsed 
!     integer, dimension(6)                  :: idate_ref
!     integer, dimension(6)                  :: idate_avg
!     character(len=256)                     :: notes 
!     character(len=1024)                    :: disclaimer
!     character(len=256)                     :: email
!     character(len=256)                     :: url
!     character(len=256)                     :: institution
!     character(len=256)                     :: conventions
!     character(len=256)                     :: history
!     character(len=256)                     :: source 
!     character(len=256)                     :: sysdate 
!     character(len=256)                     :: progstring
!     integer, dimension(8)                  :: isysdate

!     integer                                :: nstationsRegion
!     integer, dimension(:), allocatable     :: idxRegion
!     real, dimension(:,:,:), allocatable    :: ncArr ! collect output for ncdf writes
!     real, dimension(:,:,:), allocatable    :: ncArr_co2_total ! collect total co2 for ncdf writes

!     ! --- constants -----------------------

!     character(len=*), parameter            :: rname = mname//'/user_output_column_write'
!     real*4                                 :: fillval_r4 = -1e34

!     ! --- begin -------------------------------

!     if(myid /= root) return

!     ! subset column list by region
!     nstationsRegion = count( stations%region == region )

!     allocate(idxRegion(nstationsRegion))
!     allocate(stationsRegion(nstationsRegion))

!     ii = 0 
!     do i = 1, nstations ! there must be a better way to do this... 
!       if (stations(i)%region == region) then
!         ii = ii + 1
!         idxRegion(ii) = i
!       end if
!     end do

!     stationsRegion = stations(idxRegion)

!     call tau2date(itaur(region), timestamp(region)%idate_end)

!     call date_and_time(values = isysdate)
!     write (sysdate, '(i4.4,"-",i2.2,"-",i2.2," ",i2.2,":",i2.2,":",i2.2," UTC")') &
!            isysdate(1), isysdate(2), isysdate(3), isysdate(5), isysdate(6), isysdate(7)

!     ! Standard file information for CarbonTracker output.
!     notes = "This file contains CarbonTracker time averaged mole fractions for vertical profiles at select sampling locations."
!     disclaimer = "CarbonTracker is an open product of the NOAA Earth System Research \n"//&
!                  "Laboratory using data from the Global Monitoring Division greenhouse \n"//&
!                  "gas observational network and collaborating institutions.  Model results \n"//&
!                  "including figures and tabular material found on the CarbonTracker \n"//&
!                  "website may be used for non-commercial purposes without restriction,\n"//&
!                  "but we request that the following acknowledgement text be included \n"//&
!                  "in documents or publications made using CarbonTracker results: \n\n"//&
!                  "     CarbonTracker results provided by NOAA/ESRL,\n"//&
!                  "     Boulder, Colorado, USA, http://carbontracker.noaa.gov\n\n"//&
!                  "Since we expect to continuously update the CarbonTracker product, it\n"//&
!                  "is important to identify which version you are using.  To provide\n"//&
!                  "accurate citation, please include the version of the CarbonTracker\n"//&
!                  "release in any use of these results.\n\n"//&
!                  "The CarbonTracker team welcomes special requests for data products not\n"//&
!                  "offered by default on this website, and encourages proposals for\n"//&
!                  "collaborative activities.  Contact us at carbontracker.team@noaa.gov.\n"
!     email = "carbontracker.team@noaa.gov"
!     url = "http://carbontracker.noaa.gov"
!     institution = "NOAA Earth System Research Laboratory"
!     conventions = "CF-1.5"
!     call getarg (0, progstring, status)
!     write(history,'("Created ",a," by ",a,".")') trim(sysdate),trim(progstring)
!     source = "CarbonTracker release CT2012"

!     ! date/time conversions
!     call date2tau(timestamp(region)%idate_start, itau_start)
!     call date2tau(timestamp(region)%idate_end, itau_end)

!     itau_avg = nint(5.0D-1 * dble(itau_start + itau_end))

!     call tau2date(itau_avg, idate_avg)
!     ddate_avg = idate2ddate(idate_avg)

!     idate_ref = (/2000, 1, 1, 0, 0, 0/)

!     call date2tau(idate_ref, itau_ref)
!     itau_elapsed = itau_avg - itau_ref
!     days_elapsed = dble(itau_elapsed) / 86400.0D+0
!     write (outfile, '(a,"/column_",a,"_",i4.4,4i2.2,"_",i4.4,4i2.2,".nc")') &
!            trim(outdir), trim(region_name(region)), timestamp(region)%idate_start(1:5), timestamp(region)%idate_end(1:5)

!     ! create new file
!     call MDF_Create( trim(outfile), MDF_NETCDF, MDF_REPLACE, fid, status )
!     IF_ERROR_RETURN(status=1)

!     ! global attributes 
!     call MDF_Put_Att( fid, MDF_GLOBAL, "notes", values=trim(notes), status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, MDF_GLOBAL, "disclaimer", values=trim(disclaimer), status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, MDF_GLOBAL, "email", values=trim(email), status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, MDF_GLOBAL, "url", values=trim(url), status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, MDF_GLOBAL, "institution", values=trim(institution), status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, MDF_GLOBAL, "conventions", values=trim(conventions), status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, MDF_GLOBAL, "history", values=trim(history), status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, MDF_GLOBAL, "source", values=trim(source), status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, MDF_GLOBAL, "region_name", values=trim(region_name(region)), status=status )
!     IF_ERROR_RETURN(status=1)

!     ! define dimensions
!     call MDF_Def_Dim( fid, 'date', MDF_UNLIMITED, dimid_date, status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Def_Dim( fid, 'calendar_components', 6, dimid_cal, status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Def_Dim( fid, "level", lm(region), dimid_lvl, status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Def_Dim( fid, "boundary", lm(region)+1, dimid_bnd, status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Def_Dim( fid, "station", nstationsRegion, dimid_stations, status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Def_Dim( fid, "character_length", 50, dimid_nchar, status ) ! must be GE character lengths defined in column list derived type 
!     IF_ERROR_RETURN(status=1)

!     ! dimension variables
!     call MDF_Def_Var( fid, 'date', MDF_DOUBLE, (/dimid_date/), varid_elapsed, status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_elapsed, "units", values="days since 2000-01-01 00:00:00 UTC", status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_elapsed, "long_name", values="date", status=status )
!     IF_ERROR_RETURN(status=1)

!     call MDF_Def_Var( fid, 'decimal_date', MDF_DOUBLE, (/dimid_date/), varid_dec, status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_dec, "units", values="years", status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_dec, "_FillValue", values=-1e34, status=status )
!     IF_ERROR_RETURN(status=1)

!     call MDF_Def_Var( fid, 'calendar_components', MDF_INT, (/dimid_cal/), varid_cal_int, status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_cal_int, "units", values="none", status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_cal_int, "long_name", values="calendar components", status=status )
!     IF_ERROR_RETURN(status=1)

!     call MDF_Def_Var( fid, 'date_components', MDF_INT, (/dimid_cal, dimid_date/), varid_date_int, status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_date_int, "units", values="none", status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_date_int, "long_name", values="integer value calendar components of UTC date", status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_date_int, "comment", values="year, month, day, hour, minute, second", status=status )
!     IF_ERROR_RETURN(status=1)

!     call MDF_Def_Var( fid, 'station_index', MDF_INT, (/dimid_stations/), varid_column_int, status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_column_int, "units", values="none", status=status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_column_int, "long_name", values="station index", status=status )
!     IF_ERROR_RETURN(status=1)

!     call MDF_Def_Var( fid, 'lon', MDF_DOUBLE, (/dimid_stations/), varid_lon, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lon, "units", values="degrees east", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lon, "long_name", values="longitude", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lon, "actual_range", values=(/xbeg(region), xend(region)/), status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'lat', MDF_DOUBLE, (/dimid_stations/), varid_lat, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lat, "units", values="degrees north", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lat, "long_name", values="latitude", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lat, "actual_range", values=(/ybeg(region), yend(region)/), status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'alt', MDF_DOUBLE, (/dimid_stations/), varid_hgt, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_hgt, "units", values="meters", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_hgt, "long_name", values="altitude", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_hgt, "actual_range", values=(/zbeg(region), zend(region)/), status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'level', MDF_INT, (/dimid_lvl/), varid_lvl, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lvl, "units", values="none", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lvl, "long_name", values="level", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_lvl, "positive", values="up", status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'boundary', MDF_INT, (/dimid_bnd/), varid_bnd, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_bnd, "units", values="none", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_bnd, "long_name", values="boundary", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_bnd, "positive", values="up", status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'station_names', MDF_CHAR, (/dimid_nchar, dimid_stations/), varid_desc, status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_desc, "long_name", values="station names", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_desc, "standard_name", values="station description", status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'station_id', MDF_CHAR, (/dimid_nchar, dimid_stations/), varid_id, status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_id, "long_name", values="station id", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_id, "standard_name", values="station code", status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'oro', MDF_FLOAT, (/dimid_stations/), varid_oro, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_oro, "units", values="m^s/s^2", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_oro, "long_name", values="orography", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_oro, "standard_name", values="surface geopotential", status=status )
!     IF_NOTOK_RETURN(status=1)

!     call MDF_Def_Var( fid, 'co2', MDF_FLOAT, (/dimid_stations, dimid_lvl, dimid_date/), varid_total, status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_total, "units", values="micromol mol-1", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_total, "long_name", values="mole fraction of carbon dioxide in air", status=status )
!     IF_NOTOK_RETURN(status=1)
!     call MDF_Put_Att( fid, varid_total, "comment", values="total carbon dioxide estimate", status=status )
!     IF_NOTOK_RETURN(status=1)

!     do i = 1, nspecies 
!       call MDF_Def_Var( fid, trim(names(species(i))), MDF_FLOAT, (/dimid_stations, dimid_lvl, dimid_date/), varid_species(i, writeMethod), status )
!       IF_NOTOK_RETURN(status=1)
!       call MDF_Put_Att( fid, varid_species(i, writeMethod), "units", values='micromol mol-1', status=status )
!       IF_NOTOK_RETURN(status=1)
!       call MDF_Put_Att( fid, varid_species(i, writeMethod), "_FillValue", values=fillval_r4, status=status )
!       IF_NOTOK_RETURN(status=1)
!       call MDF_Put_Att( fid, varid_species(i, writeMethod), "long_name", values='mole fraction of carbon dioxide in air', status=status )
!       IF_NOTOK_RETURN(status=1)
!       call MDF_Put_Att( fid, varid_species(i, writeMethod), "comment", values=comment_species(i), status=status )
!       IF_NOTOK_RETURN(status=1)
!     enddo

!     do i = 1, nmeteoLvl
!       if (i .eq. 4) then
!         call MDF_Def_Var( fid, name_meteoLvl(i), MDF_FLOAT, (/dimid_stations, dimid_date/), varid_meteoLvl(i), status )
!       else
!         call MDF_Def_Var( fid, name_meteoLvl(i), MDF_FLOAT, (/dimid_stations, dimid_lvl, dimid_date/), varid_meteoLvl(i), status )
!       endif
!       IF_NOTOK_RETURN(status=1)
!       call MDF_Put_Att( fid, varid_meteoLvl(i), "units", values=units_meteoLvl(i), status=status )
!       IF_NOTOK_RETURN(status=1)
!       call MDF_Put_Att( fid, varid_meteoLvl(i), "long_name", values=long_name_meteoLvl(i), status=status )
!       IF_NOTOK_RETURN(status=1)
!       call MDF_Put_Att( fid, varid_meteoLvl(i), "standard_name", values=standard_name_meteoLvl(i), status=status )
!       IF_NOTOK_RETURN(status=1)
!     enddo

!     do i = 1, nmeteoBnd
!       call MDF_Def_Var( fid, name_meteoBnd(i), MDF_FLOAT, (/dimid_stations, dimid_bnd, dimid_date/), varid_meteoBnd(i), status )
!       IF_NOTOK_RETURN(status=1)
!       call MDF_Put_Att( fid, varid_meteoBnd(i), "units", values=units_meteoBnd(i), status=status )
!       IF_NOTOK_RETURN(status=1)
!       call MDF_Put_Att( fid, varid_meteoBnd(i), "long_name", values=long_name_meteoBnd(i), status=status )
!       IF_NOTOK_RETURN(status=1)
!       call MDF_Put_Att( fid, varid_meteoBnd(i), "standard_name", values=standard_name_meteoBnd(i), status=status )
!       IF_NOTOK_RETURN(status=1)
!     enddo

!     ! finished definition
!     call MDF_EndDef( fid, status )
!     IF_ERROR_RETURN(status=1)

!     ! write variables
!     call MDF_Put_Var( fid, varid_elapsed, (/(days_elapsed)/), status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_dec, (/(ddate_avg)/), status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_cal_int, (/(i,i=1,size(idate_avg))/), status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_date_int, idate_avg, status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_desc, (/(trim(stationsRegion(i)%desc),i=1,nstationsRegion)/), status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_id, (/(trim(stationsRegion(i)%id),i=1,nstationsRegion)/), status )
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_column_int, (/(i,i=1,size(stationsRegion))/), status)
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_lat, stationsRegion%lat, status)
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_lon, stationsRegion%lon, status)
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_lvl, (/(i,i=1,lm(region))/), status)
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_bnd, (/(i,i=1,lm(region)+1)/), status)
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_hgt, stationsRegion%height, status)
!     IF_ERROR_RETURN(status=1)
!     call MDF_Put_Var( fid, varid_oro, stationsRegion%oro, status)
!     IF_ERROR_RETURN(status=1)

!     allocate(ncArr(nstationsRegion, lm(region), 1))
!     allocate(ncArr_co2_total(nstationsRegion, lm(region), 1))

!     ncArr_co2_total = 0.0

!     do i = 1, nspecies 
!       do q = 1, nstationsRegion
!         ncArr(q,:,1) = stationsRegion(q)%c_avg(i,writeMethod,:)
!         ! subtract the background tracer from the others
!         if (i .gt. 1) then
!           ncArr(q,:,1) = ncArr(q,:,1) - stationsRegion(q)%c_avg(1,writeMethod,:)
!         endif
!         ncArr_co2_total(q,:,1) = ncArr_co2_total(q,:,1) + ncArr(q,:,1)
!       enddo
!       call MDF_Put_Var( fid, varid_species(i,writeMethod), ncArr, status )
!       IF_ERROR_RETURN(status=1)
!     enddo

!     call MDF_Put_Var( fid, varid_total, ncArr_co2_total, status )
!     IF_ERROR_RETURN(status=1)

!     deallocate(ncArr_co2_total)

!     do i = 1, nmeteoLvl
!       do q = 1, nstationsRegion
!         ncArr(q,:,1) = stationsRegion(q)%meteoLvl(i,:)
!       enddo
!       if (i .eq. 4) then
!         call MDF_Put_Var( fid, varid_meteoLvl(i), ncArr(:,1,:), status )
!       else
!         call MDF_Put_Var( fid, varid_meteoLvl(i), ncArr, status )
!       endif
!       IF_ERROR_RETURN(status=1)
!     enddo

!     deallocate(ncArr)

!     allocate(ncArr(nstationsRegion, lm(region)+1, 1))

!     do i = 1, nmeteoBnd
!       do q = 1, nstationsRegion
!         ncArr(q,:,1) = stationsRegion(q)%meteoBnd(i,:)
!       enddo
!       call MDF_Put_Var( fid, varid_meteoBnd(i), ncArr, status )
!       IF_ERROR_RETURN(status=1)
!     enddo

!     deallocate(ncArr)

!     ! close file
!     call MDF_Close( fid, status )
!     IF_ERROR_RETURN(status=1)

!     write (gol, '("'//trim(rname)//' Done writing arrays and closing output file.")'); call goPr

!     timestamp(region)%idate_start = timestamp(region)%idate_end ! end of interval becomes start of next interval  

!     deallocate(idxRegion)
!     deallocate(stationsRegion)

    ! ok
    status = 0

  end subroutine user_output_column_write

  subroutine user_output_column_reset( region )

!     use ParTools, only: root_t, myid

!     implicit none
  
    integer, intent(in)        :: region

!     integer                    :: q

!     if (myid /= root_t) return

!     do q=1, nstations
!        stations(q)%c_avg(:,:,:) = 0.0
!        stations(q)%c_std(:,:,:) = 0.0
!        stations(q)%meteoLvl(:,:) = 0.0
!        stations(q)%meteoBnd(:,:) = 0.0
!        stations(q)%ndata = 0
!        stations(q)%accum_weight = 0.0
!     end do

!     evaluated = .FALSE.

  end subroutine user_output_column_reset 

  subroutine user_output_column_done( status )

!     use Partools,         only : myid, root
!     use dims,             only : nregions
 
!     implicit none

    ! --- in/out ------------------------

    integer, intent(out)  :: status

    ! --- variables ---------------------------
    
!     integer               :: region
!     integer               :: q 

!     ! --- constants ---------------------------
    
!     character(len=*), parameter  ::  rname = mname//'/user_output_column_done'

!     ! --- begin -------------------------------

!     if (myid /= root) return

!     do region = 1, nregions
!        call user_output_column_evaluate( status )
!        IF_NOTOK_RETURN(status=1)
!        call user_output_column_write( region, status )
!        IF_NOTOK_RETURN(status=1)
!     end do

!     do q = 1, nstations
!        deallocate(stations(q)%c_avg)
!        deallocate(stations(q)%c_std)
!     enddo

!     deallocate(stations)

    ! ok
    status = 0

  end subroutine user_output_column_done

end module user_output_column