!### macros #####################################################
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
                                                                                                                                                     
#include "tm5.inc"
!#################################################################

module user_output_flask
!BOP
!
! !MODULE:  user_output_flask
! 
! !PUBLIC TYPES: none
!
! !PUBLIC MEMBER FUNCTIONS:
!    user_output_flask_init
!    user_output_flask_sample
!    user_output_flask_done
!
! !PUBLIC DATA MEMBERS:
!
! !DESCRIPTION: 
!
!   Routine for simulating flask samples in TM5. 
!
!   "Flask" is shorthand for any x,y,z,t sampling of the model.  This
!   routine is intended for producing model output of low-frequency,
!   intermittent, or sporadic sampling (like NOAA Cooperative Sampling
!   Network flasks), but is also useful for other observations like
!   aircraft profiles or daily averages from continuous station data.
!
!   I/O is via netCDF files, using the MDF interface.  Observation
!   list from input file is truncated to begin/end times of the current run.
!
!   Input file format:
!     Routine expects to read a netCDF (v3) file with:
!
!     "obs_num"   : record dimension, contains unique integers
!     "obs_id"    : 30-character strings also identifying obs (redundant with obs_num)
!     "latitude"  : 1-d real array, length id, of latitudes (degrees N)
!     "longitude" : 1-d real array, length id, of longitudes (degrees E)
!     "altitude"  : 1-d real array, length id, of altitudes (m ASL)
!
!     "calendar_components": 2-D variable with dims (id,6),
!        with time of observations in 6-integer 'idate' format
!        (year, month, day, hour, minute, second in UTC time.
!
!     "sampling_strategy" 1-D integer array designed to represent
!        distinct schemes for sampling TM5 to represent the
!        observations in this file.  Add to the following list if
!        you define a new scheme:
!
!        VALUE                  MEANING
!
!          1     Four-hour averages, centered on observation time.
!          2     One-hour averages, centered on observation time.
!          3     90-minute averages, centered on observation time.
!
!   Controlling rc-file keys:
!
!     output.flask  (logical, default F)
!         Whether or not this routine is active.  Note that
!         unlike other rc-file keys, this one is read in 
!         user_output.F90.  Sets flask_data (a public logical
!         in that module).
!
!     output.flask.infile  (string, required, no default value)
!         Path to input netCDF file, e.g. "/path/to/input/flask_obs.nc"
!
!     output.flask.verbose (logical, optional, default F)
!         If true, some extra information is printed to stdout.
!
!     output.flask.replicate.forecast (logical, optional, default F)
!         If true, mimic the sampling scheme of user_output_forecast.F90:
!
!         (1) obs are in-window for this simulation if itaui <= obs < itaue,
!             where itaui and itaue are the initial and ending simulation
!             time values, and obs is the observation center time.  The 
!             default behavior of this routine is instead to sample all
!             obs meeting the criterion itaui-window/2 <= obs <= itaue+window/2.  
!             This will result in obs being sampled in multiple simulations
!             when their sampling windows cross itaui or itaue.  This
!             must be dealt with in post-processing.
!
!
! !REVISION HISTORY: 
!
!   Andy Jacobson, Aug 2010
!
!   Adapted from code in user_output_forecast, user_output_station,
!   and user_output_noaa, mostly if not all originally written by
!   Wouter Peters.  Also used Arjo Segers' file_MDF example code
!   for netCDF I/O.
!
!   Andy Jacobson, Aug 2011
!
!   Changed variable names to be consistent with Ken Masarie's new
!   ObsPack data distribution system.
!  
!
!EOP
!-------------------------------------------------------------------------


#ifdef MPI
  use mpi_const,			only : my_real,mpi_comm_world,ierr,mpi_integer,mpi_character, com_trac
  use mpi_const,			only : MPI_SUCCESS, MPI_INFO_NULL
  use mpi_comm,				only : barrier
#endif
  
  use ParTools,				only : ntracetloc, myid, root, npes, ntracet_ar, par_barrier
  use chem_param,			only : ntracet
  use MDF      , 			only : MDF_Open
  use MDF      , 			only : MDF_Create, MDF_Close, MDF_EndDef
  use MDF      , 			only : MDF_NETCDF4, MDF_REPLACE, MDF_GLOBAL, MDF_CHAR, MDF_INT, MDF_FLOAT, MDF_UNLIMITED, MDF_READ
  use MDF      , 			only : MDF_Put_Att
  use MDF      , 			only : MDF_Def_Dim
  use MDF      , 			only : MDF_Inq_Dimid, MDF_Inquire_Dimension
  use MDF      , 			only : MDF_Inq_VarID
  use MDF      , 			only : MDF_Def_Var, MDF_Put_Var, MDF_Get_Var

  
  implicit none

  private

  integer :: communicator

  public 				:: user_output_flask_init
  public 				:: user_output_flask_done
  public 				:: user_output_flask_sample

  character(len=1024)			:: inFile = ''
  character(len=1024)			:: outFile = ''

  logical				:: flask_replicate_forecast=.false.
  logical				:: flask_sample_meteo=.false.
  logical				:: flask_verbose=.false.

  integer				:: nflasks
  real*4, parameter			:: flask_missing_value=-1.0e34

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

  ! "obs_num" is meant to be a flask identifier guaranteed to be unique, for the 
  ! purpose of matching input file records with output file records.  In 
  ! CarbonTracker, it is generated when the observations netCDF file is generated,
  ! and is the record dimension.
  !
  ! "obs_id" is a string descriptor meant to give humans an indication of what
  ! each observation represents.  For CarbonTracker, this can mean NOAA database
  ! eventnumbers or strings made up to represent 4-hour averages from continuous
  ! instruments.  It is only present to help identify observations.

  type flask_sample
     integer 				:: obs_num  ! unique indentifier placed in input file
     real				:: lat ! sample latitude
     real				:: lon ! sample longitude
     real				:: alt ! sample altitude (meters above sea level)
     integer				:: itau_start ! sampling start time
     integer				:: itau_center ! sampling center time
     integer				:: itau_end ! sampling end time
     character(len=30)			:: obs_id ! another sample identifier; could be NOAA's CCG database event number
     real, dimension(:), allocatable	:: mix ! mixing ratio of interest sampled using x-y-z slopes
     real, dimension(:), allocatable	:: mix_grd ! mixing ratio of interest sampled using grid-box values
     integer				:: nsamples ! number of samples accumulated
     real				:: accum_weight ! accumulated weight
     integer				:: region ! zoom region index from which this flask is sampled
     integer				:: ifr, jfr ! i,j region indices for flask's grid cell
     integer				:: ifn, jfn ! i,j region indices for flask's "next" grid cell
     real				:: rif, rjf ! fractions from center of ifr,jfr box
     real				:: surface_height ! surface height in meters
     integer				:: lfr ! vertical level number of the sample
     real				:: wcx, wcy ! x and y weighting factors for slopes interpolation
     logical				:: evaluated ! whether or not the average has been computed from mix/nsamples
     real				:: u,v,blh,q,pressure,temperature ! meteorological variables
     logical            :: below_surface_warning = .False.
  end type flask_sample

  type(flask_sample), dimension(:), allocatable  :: flasks




contains

!BOP
!
! !IROUTINE:  user_output_flask_init
!
! !INPUT PARAMETERS: 
!   rcF - handle to open rc file
!   status - integer status indicator
!
! !OUTPUT PARAMETERS: none
!
! !DESCRIPTION: 
!   Initialization tasks for user_output_flask:
!    - reads rc file for flask keys
!    - allocates flasks array
!    - reads input netCDF file
!    - computes x,y grid indices and weights for slopes interpolation
!    - transmits flasks array to all PEs
!
!
!EOP

  subroutine user_output_flask_init(rcF, status)

    use GO,			only : TrcFile, ReadRc
    use GO,			only : pathsep
    use GO,			only : gol, goErr, goPr, goBug

    use global_data, 		only : outdir

    use Meteo,			only : Set
    use Meteo,			only : gph_dat
    use binas,			only : ae, twopi, grav
    use dims,			only : im, jm, lm, dx, dy, xref, yref, xbeg, ybeg, xend, yend
    use dims,			only : nregions, region_name, itaui, itaue, xcyc, idate
    use chem_param,		only : ntrace, ntracet, names, fscale
    use toolbox,	 	only : escape_tm
    use datetime,		only : date2tau, tau2date

    use ParTools

    implicit none

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

    type(TrcFile), intent(in)			:: rcF
    integer, intent(out)			:: status
    
    ! local
    integer					:: hid,dimid
    character(len=256)				:: name
    logical,dimension(:), allocatable		:: mask
    integer,dimension(:), allocatable		:: obs_num
    integer,dimension(:,:), allocatable		:: idate_f
    real, dimension(:), allocatable		:: lat
    real, dimension(:), allocatable		:: lon
    real, dimension(:), allocatable		:: alt
    integer, dimension(:), allocatable		:: sampling_strategy
    integer, dimension(:), allocatable		:: itau_center
    integer, dimension(:), allocatable		:: itau_start
    integer, dimension(:), allocatable		:: itau_end
    character*30, dimension(:), allocatable	:: obs_id
    integer, dimension(6)			:: idatei,idatee,idatef
    logical					:: new_region
    integer					:: region
    integer					:: iflask
    integer					:: varid
    character(len=15)				:: idstr
    logical					:: input_file_exists

    ! x,y resolution (in degrees) for current region 
    real					:: dxr, dyr 

    ! --- const ------------------------------

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

    do region=1,nregions

       call Set(	 gph_dat(region), status, used=.true. )
       IF_NOTOK_RETURN(status=1)


    enddo

    if(myid==root) then

       call ReadRc(rcF, 'output.flask.infile', inFile, status)
       IF_NOTOK_RETURN(status=1)

       inFile=trim(inFile)

       call ReadRc(rcF, 'output.flask.verbose', flask_verbose, status, default = .false.)
       IF_NOTOK_RETURN(status=1)

       call ReadRc(rcF, 'output.flask.meteo', flask_sample_meteo, status, default = .false.)
       IF_NOTOK_RETURN(status=1)

       call ReadRc(rcF, 'output.flask.replicate.forecast', flask_replicate_forecast, status, default = .false.)
       IF_NOTOK_RETURN(status=1)

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

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

       write (gol,'("[user_output_flask_init]  input from ",a,".")') trim(inFile); call goPr
       write (gol,'("[user_output_flask_init]  output to ",a,".")') trim(outFile); call goPr

       if(flask_verbose) then
          write (gol,'("[user_output_flask_init]  verbose output requested.")'); call goPr
       endif

       if(flask_sample_meteo) then
          write (gol,'("[user_output_flask_init]  meteo variables will also be sampled.")'); call goPr
       endif

       if(flask_replicate_forecast) then
          write (gol,'("[user_output_flask_init]  sampling will replicate user_output_forecast scheme.")'); call goPr
       endif

    endif

    ! read input file

    ! get from nc input file:
    !	lat, lon, alt, date_components, sampling_strategy, obs_id
    !   restrict to flasks between itaui and itaue
    !   for each flask, determine region, ifr, jfr, ...things that don't change.
    if(myid==root) then

       inquire(file=inFile,exist=input_file_exists)

       if (input_file_exists) then

          call MDF_Open(filename=inFile, ftype=MDF_NETCDF4, mode=MDF_READ, hid=hid, status=status)
          IF_NOTOK_RETURN(status=1)

          call MDF_Inq_DimID(hid, "obs", dimid, status)
          if (status /= 0) write (gol,'("[user_output_flask_init] WARNING: dimension in the observation input file not recognized")') ; call goPr
          if (status /= 0) write (gol,'("[user_output_flask_init] WARNING: are py-CTDAS and tm5-CTDAS  compatible?")') ; call goPr
          if (status /= 0) write (gol,'("[user_output_flask_init] WARNING: check method WriteSampleInfo in obs.py or obspack.py in da/ct/")') ; call goPr
          IF_NOTOK_RETURN(status=1)

          call MDF_Inquire_Dimension(hid, dimid, status, length=nflasks)
          IF_NOTOK_RETURN(status=1)

          write (gol,'("[user_output_flask_init] ",i," obs in input file.")') nflasks; call goPr

          allocate(idate_f(6,nflasks))
          allocate(mask(nflasks))
          allocate(obs_num(nflasks))
          allocate(itau_start(nflasks))
          allocate(itau_center(nflasks))
          allocate(itau_end(nflasks))
          allocate(lat(nflasks))
          allocate(lon(nflasks))
          allocate(alt(nflasks))
          allocate(sampling_strategy(nflasks))
          allocate(obs_id(nflasks))

          call MDF_Inq_VarID(hid, "obs_num", varid, status)
          IF_NOTOK_RETURN(status=1)

          call MDF_Get_Var(hid,varid,obs_num,status)
          IF_NOTOK_RETURN(status=1)

          call MDF_Inq_VarID(hid, "date_components", varid, status)
          IF_NOTOK_RETURN(status=1)

          call MDF_Get_Var(hid,varid,idate_f,status)
          IF_NOTOK_RETURN(status=1)

          call MDF_Inq_VarID(hid, "latitude", varid, status)
          IF_NOTOK_RETURN(status=1)

          call MDF_Get_Var(hid,varid,lat,status)
          IF_NOTOK_RETURN(status=1)

          call MDF_Inq_VarID(hid, "longitude", varid, status)
          IF_NOTOK_RETURN(status=1)

          call MDF_Get_Var(hid,varid,lon,status)
          IF_NOTOK_RETURN(status=1)
          
          call MDF_Inq_VarID(hid, "altitude", varid, status)
          IF_NOTOK_RETURN(status=1)
          
          call MDF_Get_Var(hid,varid,alt,status)
          IF_NOTOK_RETURN(status=1)
          
          call MDF_Inq_VarID(hid, "sampling_strategy", varid, status)
          IF_NOTOK_RETURN(status=1)
          
          call MDF_Get_Var(hid,varid,sampling_strategy,status)
          IF_NOTOK_RETURN(status=1)
          
          call MDF_Inq_VarID(hid, "obs_id", varid, status)
          IF_NOTOK_RETURN(status=1)
          
          call MDF_Get_Var(hid,varid,obs_id,status)
          IF_NOTOK_RETURN(status=1)
          
          call MDF_Close(hid,status)
          IF_NOTOK_RETURN(status=1)
          
          ! determine valid obs and subset arrays; fill flasks structure
          do iflask = 1,nflasks

             call date2tau(idate_f(:,iflask),itau_center(iflask))

             select case (sampling_strategy(iflask))

             case (1) ! 4-hour window
                
                itau_start(iflask) = itau_center(iflask)-2*3600
                itau_end(iflask) = itau_center(iflask)+2*3600

             case (2) ! 1-hour window
                
                itau_start(iflask) = itau_center(iflask)-1800
                itau_end(iflask) = itau_center(iflask)+1800

             case (3) ! 90-minute window
                
                itau_start(iflask) = itau_center(iflask)-2700
                itau_end(iflask) = itau_center(iflask)+2700

             case default
                
                write (gol,'("[user_output_flask_init] Flask with obs_num ",i," and obs_id string ",a,":")') &
                     obs_num(iflask), trim(adjustl(obs_id(iflask)))
                call goPr
                write (gol, '("  Unknown sampling strategy = ",i10,".")') sampling_strategy(iflask)
                call goErr
                
                status=1
                return
                
             end select


          enddo

          if(flask_replicate_forecast) then
             mask=((itau_center .ge. itaui) .and. (itau_center .lt. itaue))  ! sample only itaui <= obs < itaue
          else
             mask=((itau_end .gt. itaui) .and. (itau_start .lt. itaue))  ! sample all obs whose sampling windows fall inside (itaui,itaue)
          endif

          nflasks=count(mask)  ! note that this changes the value of nflasks
          write (gol,'("[user_output_flask_init] ",i," obs in time range ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," to ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2,".")') &
               nflasks,idatei(1:5),idatee(1:5); call goPr

       else

          nflasks = 0
          write (gol,'("[user_output_flask_init]  input file ",a," does not exist: nflasks set to zero.")') trim(inFile) ; call goPr

       endif
    endif ! myid is root
    
    call par_barrier

    ! jump out to broadcast the final nflasks and get all PEs to allocate array
    
#ifdef MPI	 
    communicator = com_trac
    call MPI_BCAST(nflasks, 1, MPI_INTEGER, root, communicator, ierr)  
#endif

    if(nflasks .eq. 0) then
       if(allocated(idate_f)) deallocate(idate_f)
       if(allocated(mask)) deallocate(mask)
       if(allocated(obs_num)) deallocate(obs_num)
       if(allocated(itau_center)) deallocate(itau_center)
       if(allocated(itau_start)) deallocate(itau_start)
       if(allocated(itau_end)) deallocate(itau_end)
       if(allocated(lat)) deallocate(lat)
       if(allocated(lon)) deallocate(lon)
       if(allocated(alt)) deallocate(alt)
       if(allocated(sampling_strategy)) deallocate(sampling_strategy)
       if(allocated(obs_id)) deallocate(obs_id)
       if(allocated(flasks)) deallocate(flasks)
       status=0
       return
    endif

    if(allocated(flasks)) deallocate(flasks)  ! avoid double allocation

    allocate(flasks(nflasks))
        
    !AB? Make flag .false. at startup
    flasks%evaluated=.false.
    
    do iflask=1,nflasks
       allocate(flasks(iflask)%mix(ntracetloc))
       allocate(flasks(iflask)%mix_grd(ntracetloc))
       flasks(iflask)%mix = 0.0
       flasks(iflask)%mix_grd = 0.0
       flasks(iflask)%nsamples = 0
       flasks(iflask)%accum_weight = 0.0
       if(flask_sample_meteo) then
          flasks(iflask)%u = 0.0
          flasks(iflask)%v = 0.0
          flasks(iflask)%blh = 0.0
          flasks(iflask)%q = 0.0
          flasks(iflask)%pressure = 0.0
          flasks(iflask)%temperature = 0.0
       endif
    enddo

    ! root only does these next computations, will be broadcast later

    if(myid == root) then

       flasks%obs_num=pack(obs_num,mask)
       flasks%lat=pack(lat,mask)
       flasks%lon=pack(lon,mask)
       flasks%alt=pack(alt,mask)
       flasks%obs_id=pack(obs_id,mask)
       flasks%itau_start=pack(itau_start,mask)
       flasks%itau_center=pack(itau_center,mask)
       flasks%itau_end=pack(itau_end,mask)

       ! initialize structure with default and undefined values
       flasks%region=-1
       flasks%ifr = -1
       flasks%jfr = -1
       flasks%rif = -1e12
       flasks%rjf = -1e12
       flasks%ifn = -1
       flasks%jfn = -1
       flasks%wcx = -1e12
       flasks%wcy = -1e12

       do iflask=1,nflasks

          write(idstr,'(i8)') flasks(iflask)%obs_num
          call tau2date(flasks(iflask)%itau_center,idatef)

          if(flask_replicate_forecast) then
             ! per user_output_flask, move obs within one second of itaui and itaue towards the middle of the
             ! simulation by one hour.
             if(abs(flasks(iflask)%itau_center-itaui).lt.1) then
                write (gol,'("[user_output_flask_init]  attention: flask ",a," (obs_id ",a,"): ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," at start; moved forward one hour.")') &
                     trim(adjustl(idstr)), trim(adjustl(flasks(iflask)%obs_id)), idatef(1:5)
                call goPr
                flasks(iflask)%itau_center = flasks(iflask)%itau_center+3600
                flasks(iflask)%itau_start = flasks(iflask)%itau_start+3600
                flasks(iflask)%itau_end = flasks(iflask)%itau_end+3600
             endif
             if(abs(flasks(iflask)%itau_center-itaue).lt.1) then
                write (gol,'("[user_output_flask_init]  attention: flask ",a," (obs_id ",a,"): ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," at end; moved back one hour.")') &
                     trim(adjustl(idstr)), trim(adjustl(flasks(iflask)%obs_id)), idatef(1:5)
                call goPr
                flasks(iflask)%itau_center=flasks(iflask)%itau_center-3600
                flasks(iflask)%itau_start = flasks(iflask)%itau_start-3600
                flasks(iflask)%itau_end = flasks(iflask)%itau_end-3600
             endif
          endif

          if(flasks(iflask)%itau_start .lt. itaui) then

             if(flask_replicate_forecast) then
                write (gol,'("[user_output_flask_init]  attention: flask ",a," (obs_id ",a,"): ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," close to start; sampling truncated.")') &
                     trim(adjustl(idstr)), trim(adjustl(flasks(iflask)%obs_id)), idatef(1:5)
             else
                write (gol,'("[user_output_flask_init]  attention: flask ",a," (obs_id ",a,"): ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," close to start; sampling continued from previous run.")') &
                     trim(adjustl(idstr)), trim(adjustl(flasks(iflask)%obs_id)), idatef(1:5)
             endif
             call goPr

             ! change itau_start so that averaging time in output file reflects actual averaging period
             flasks(iflask)%itau_start = itaui

          endif

          if(flasks(iflask)%itau_end .gt. itaue) then

             call tau2date(flasks(iflask)%itau_center,idatef)
             
             if(flask_replicate_forecast) then
                write (gol,'("[user_output_flask_init]  attention: flask ",a," (obs_id ",a,"): ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," close to end; sampling truncated.")') &
                     trim(adjustl(idstr)), trim(adjustl(flasks(iflask)%obs_id)), idatef(1:5)
             else
                write (gol,'("[user_output_flask_init]  attention: flask ",a," (obs_id ",a,"): ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," close to end; sampling will continue in subsequent run.")') &
                     trim(adjustl(idstr)), trim(adjustl(flasks(iflask)%obs_id)), idatef(1:5)
             endif
             call goPr

             ! change itau_end so that averaging time in output file reflects actual averaging period
             flasks(iflask)%itau_end = itaue

          endif

          ! assign region and compute indices and slopes weighting factors

          do region=1, nregions
             
             dyr = dy/yref(region)
             dxr = dx/xref(region)

             if ((flasks(iflask)%lon .gt. xbeg(region)) .and. &
                  (flasks(iflask)%lon .le. xend(region)) .and. &
                  (flasks(iflask)%lat .gt. ybeg(region)) .and. &
                  (flasks(iflask)%lat .le. yend(region))) then


                ! Flask is in region

                ! Check whether an existing region assignment exists.
                ! If so, assign this region only if it is has a finer
                ! grid than existing region.

                new_region = .FALSE.

                if (flasks(iflask)%region .eq. -1) then

                   flasks(iflask)%region = region
                   new_region = .TRUE.

                else if(xref(region) > xref(flasks(iflask)%region)) then

                   flasks(iflask)%region = region
                   new_region = .TRUE.

                endif
                
                if(new_region) then

                   ! compute indices and weighting factors

                   flasks(iflask)%rif = (flasks(iflask)%lon-float(xbeg(region)))/dxr + 0.99999
                   flasks(iflask)%rjf = (flasks(iflask)%lat-float(ybeg(region)))/dyr + 0.99999

                   flasks(iflask)%ifr  = int(flasks(iflask)%rif)   ! i-index of grid cell in which observation is located
                   flasks(iflask)%jfr  = int(flasks(iflask)%rjf)   ! j-index of grid cell in which observation is located

!                   write(gol,'("   ",i8," ",f8.2," ",f8.2," (",i3,")",f8.2," ",f8.2," (",i3,")")') \
!                   flasks(iflask)%obs_num, \
!                   flasks(iflask)%lon,flasks(iflask)%rif,flasks(iflask)%ifr, \
!                   flasks(iflask)%lat,flasks(iflask)%rjf,flasks(iflask)%jfr
!                   call goPr
                   
                   !fraction from the center of the is-box  (-0.5---+0.5)
                   flasks(iflask)%rif = flasks(iflask)%rif-flasks(iflask)%ifr-0.5
                   !idem js
                   flasks(iflask)%rjf = flasks(iflask)%rjf-flasks(iflask)%jfr-0.5

                   !the neighbour for x interpolation
                   if(flasks(iflask)%rif .gt. 0) then 
                      flasks(iflask)%ifn = flasks(iflask)%ifr+1      
                   else
                      flasks(iflask)%ifn = flasks(iflask)%ifr-1
                   endif

                   !the neighbour for y interpolation
                   if(flasks(iflask)%rjf .gt. 0) then 
                      flasks(iflask)%jfn = flasks(iflask)%jfr+1      
                   else
                      flasks(iflask)%jfn = flasks(iflask)%jfr-1
                   endif
                   
                   ! x- / y-weighting of grid cell in which observation is located
                   flasks(iflask)%wcx = (1.0-abs(flasks(iflask)%rif))    ! 1.0 ... 0.5
                   flasks(iflask)%wcy = (1.0-abs(flasks(iflask)%rjf))    ! 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 (flasks(iflask)%jfn < 1) flasks(iflask)%jfn=1
                   if (flasks(iflask)%jfn > jm(region) ) flasks(iflask)%jfn=jm(region)
                   if (xcyc(region) == 0 ) then
                      ! non-cyclic boundaries
                      if (flasks(iflask)%ifn < 1) flasks(iflask)%ifn=1
                      if (flasks(iflask)%ifn > im(region) ) flasks(iflask)%ifn=im(region)
                   else
                      ! cyclic x-boundaries
                      if (flasks(iflask)%ifn < 1 ) flasks(iflask)%ifn=im(region)
                      if (flasks(iflask)%ifn > im(region) ) flasks(iflask)%ifn=1
                   endif
                   
                   
                endif ! if new_region
                
             endif ! if in region

          enddo ! regions

       enddo ! flasks

       if(flask_verbose) then
          write(gol,'("[user_output_flask_init]  list of observations to be sampled during this simulation:")');call goPr
          write(gol,'("    flask obs_num region  longitude (i)  latitude (j)  altitude         date               obs_id")');call goPr
         do iflask = 1,nflasks
            if(flasks(iflask)%region > -1) then
               call tau2date(flasks(iflask)%itau_center,idatef)
               write(gol,'("   ",i5," ",i8," ",a," ",f8.2," (",i3,")",f8.2," (",i3,")",f9.1,"   ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2,":",i2.2, " UTC   ",a)') & ! use '&' character on huygens  
               iflask,flasks(iflask)%obs_num,trim(adjustl(region_name(flasks(iflask)%region))),flasks(iflask)%lon,flasks(iflask)%ifr,flasks(iflask)%lat, &
               flasks(iflask)%jfr,flasks(iflask)%alt,idatef,trim(adjustl(flasks(iflask)%obs_id))
               call goPr
            endif
          enddo
       endif

    endif ! myid is root

    call par_barrier

#ifdef MPI	 
    communicator = com_trac
    do iflask=1, nflasks
       call MPI_BCAST(flasks(iflask)%obs_num,1,	MPI_INTEGER , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%lat,1,	MY_REAL , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%lon,1,	MY_REAL , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%alt,1,	MY_REAL , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%obs_id ,30,	MPI_CHARACTER, &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%itau_start,1,	MPI_INTEGER , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%itau_center,1,	MPI_INTEGER , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%itau_end,1,	MPI_INTEGER , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%region,1,	MPI_INTEGER , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%nsamples,1,	MPI_INTEGER , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%accum_weight,1,	MY_REAL , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%ifr,1,	MPI_INTEGER , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%jfr,1,	MPI_INTEGER , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%ifn,1,	MPI_INTEGER , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%jfn,1,	MPI_INTEGER , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%rif,1,	MY_REAL , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%rjf,1,	MY_REAL , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%wcx,1,	MY_REAL , &
            root ,communicator,ierr)  
       call MPI_BCAST(flasks(iflask)%wcy,1,	MY_REAL , &
            root ,communicator,ierr)  
    enddo
#endif

    call Par_Barrier

    if(allocated(idate_f)) deallocate(idate_f)
    if(allocated(mask)) deallocate(mask)
    if(allocated(obs_num)) deallocate(obs_num)
    if(allocated(itau_center)) deallocate(itau_center)
    if(allocated(itau_start)) deallocate(itau_start)
    if(allocated(itau_end)) deallocate(itau_end)
    if(allocated(lat)) deallocate(lat)
    if(allocated(lon)) deallocate(lon)
    if(allocated(alt)) deallocate(alt)
    if(allocated(sampling_strategy)) deallocate(sampling_strategy)
    if(allocated(obs_id)) deallocate(obs_id)

    status=0
    return

  end subroutine user_output_flask_init

!BOP
!
! !IROUTINE:  user_output_flask_sample
!
! !INPUT PARAMETERS: 
!   region - region number being sampled
!   itau_mod - current model time
!
! !OUTPUT PARAMETERS: none
!
! !DESCRIPTION: 
!    - sample model at flask locations valid for this region and time
!
!
!EOP


  subroutine user_output_flask_sample(region,itau_mod,status)

    ! This subroutine samples the model at given locations in flasks array

    use global_data, only				: mass_dat, region_dat, conv_dat 
    use Meteo, only					: gph_dat, m_dat, humid_dat, temper_dat, pu_dat, pv_dat, phlb_dat
    use dims, only 					: lm,im,jm
    use GO, only 					: gol, goErr, goPr, goBug
    use chem_param, only 				: fscale
    use datetime, only					: tau2date
    use dims, only 					: dy, dx, xref, yref, gtor
    use dims, only 					: ybeg, ndyn, ndyn_max
    use binas, only 					: ae

    implicit none

    ! input/output
    integer,intent(in)					:: region
    integer,intent(in)					:: itau_mod
    integer,intent(inout)				:: status

    real,dimension(:,:,:), pointer			:: m,gph
    real,dimension(:,:,:,:), pointer			:: rm, rxm, rym, rzm
    real,dimension(:,:), pointer			:: blh ! boundary layer height [m]
    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			:: q ! specific humidity [kg/kg]

    real,dimension(0:lm(region))			:: height
    integer						:: iflask, offsetj, lfrt, lfrtn, lmr, lfrn
    integer						:: l,itr
    integer						:: n,i,j
    integer						:: i0,i1,j0,j1,l0,l1
    integer						:: ifr,jfr,lfr,ifn,jfn
    real						:: alt,rlf, wcz
    real						:: rmf
    real						:: wcx,wcy,rif,rjf
    real						:: this_weight

    logical 						:: dbug = .false.
    logical 						:: in_window

    integer, dimension(6)				:: idatem

    real 						:: yres,xres,dyy
    real, dimension(2) 					:: lat
    real,dimension(2,2,2) 				:: u,v
    real, dimension(2) 					:: dxx
    
    ! --- const ------------------------------

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

    this_weight=real(ndyn)/real(ndyn_max)

    ! pointers to global arrays

    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
    gph => gph_dat(region)%data
    blh => conv_dat(region)%blh
    t => temper_dat(region)%data
    pu => pu_dat(region)%data
    pv => pv_dat(region)%data
    phlb => phlb_dat(region)%data
    q => humid_dat(region)%data

    lmr=lm(region)

    do iflask=1,nflasks

       !  1. Is the site in this zoom region?
       !  2. Is model time in the sampling window?
       !  3. Use slopes to determine tracer concentrations at the site.

       if (.not. flasks(iflask)%region .eq. region) cycle

       in_window = .false.

       if ((flask_replicate_forecast) .and. &
            (itau_mod .gt. flasks(iflask)%itau_start) .and. &
            (itau_mod .le. flasks(iflask)%itau_end)) then
          in_window = .true.
       endif

       if ((.not. flask_replicate_forecast) .and. &
            (itau_mod .ge. flasks(iflask)%itau_start) .and. &
            (itau_mod .le. flasks(iflask)%itau_end)) then
          in_window = .true.
       endif


       if(in_window)  then

          dbug=.false.

          alt = flasks(iflask)%alt

          ifr = flasks(iflask)%ifr
          ifn = flasks(iflask)%ifn
          jfr = flasks(iflask)%jfr
          jfn = flasks(iflask)%jfn
          rif = flasks(iflask)%rif
          rjf = flasks(iflask)%rjf
          wcx = flasks(iflask)%wcx
          wcy = flasks(iflask)%wcy

          ! interpolate the altitude to site position...
          lfr = 1 !layer
          do l=0,lm(region)
             height(l) = &
                        wcx *      wcy  * gph(ifr,jfr,l+1)  + &
                  (1.0-wcx) *      wcy  * gph(ifn,jfr,l+1) + &
                        wcx * (1.0-wcy) * gph(ifr,jfn,l+1) + &
                  (1.0-wcx) * (1.0-wcy) * gph(ifn,jfn,l+1)

             if(l==0) flasks(iflask)%surface_height = height(0)

!             write(gol,'("height(",i2,")=",f10.2," masl.")') l,height(l);call goPr

          enddo
          do l=2,lm(region)	! selects layer , note that we start from second layer from surface
             if(height(l).gt.alt) exit
          enddo

          flasks(iflask)%lfr = l

          select case(l)
          case(0)
             if(myid==root) then
                if (.not.flasks(iflask)%below_surface_warning) then
                    write (gol,'("WARNING:  For flask with obs_id ",a,":")') trim(flasks(iflask)%obs_id)
                    call goPr
                    write (gol,'("  Sample altitude of ",f8.2," masl is below surface height of ",f8.2," masl.")') &
                         alt,height(0)
                    call goPr
                    write (gol,'("  Will sample at surface.")') 
                    call goPr
                    flasks(iflask)%below_surface_warning = .True.
                endif
                lfr = 1
                rlf = -0.5	!surface...
             endif
          case default
             lfr = l  !the site layer
             ! the offset from the center of the layer (-0.5--->+0.5)
             ! (interpolation is in (m))
             rlf = (alt-height(l-1))/(height(l)-height(l-1)) - 0.5
          end select

!          write(gol,'("sample with alt=",f10.2," we choose lfr=",i2," with height ",f10.2," masl.")') alt,lfr,height(lfr);call goPr
          !=================================
          !the neighbour for z interpolation
          !=================================
          if (rlf .gt. 0 ) then 
             lfrn = lfr+1
          else
             lfrn = lfr-1
          endif
          ! z-weighting of grid cell in which observation is located 
          wcz = (1.0-abs(rlf))  !.0 ... 0.5

          !=========================================================
          ! if vertical neighbor is 0 (which does not exist)
          ! take vertical layer with l=2 for EXTRApolation to ground
          !=========================================================

          IF(lfrn == 0) THEN
             lfrn=2
             wcz=1.0-rlf  ! 1.0 ... 1.5
          ENDIF
          IF(lfrn == lmr+1) THEN
             !=========================================================
             ! if vertical neighbor is lmr+1 (which does not exist)
             ! -> no interpolation
             !=========================================================
             lfrn=lmr ! no interpolation 
             wcz=1.0
          ENDIF

          ! sample tracers

          do itr=1,ntracetloc

             ! rm-value is obtained from rm + slopes. 
             ! slope = rxm = (rm*dX/dx *deltaX/2)

             ! offsetj + itr gives the *absolute* tracer number;
             ! i.e., across all PEs.

             offsetj = sum(ntracet_ar(0:myid-1))

             ! full x,y,z slopes sampling
             rmf = (         rm(ifr,jfr,lfr,itr) + &
                  2.0*(rif*rxm(ifr,jfr,lfr,itr) + &
                  rjf*rym(ifr,jfr,lfr,itr) + &
                  rlf*rzm(ifr,jfr,lfr,itr)))/ &
                  m(ifr,jfr,lfr)*fscale(offsetj+itr) 

             flasks(iflask)%mix(itr)=flasks(iflask)%mix(itr)+this_weight*rmf !*fscale(offsetj+itr)

             ! grid-box sampling 
             rmf = rm(ifr,jfr,lfr,itr)/m(ifr,jfr,lfr)*fscale(offsetj+itr) 

             flasks(iflask)%mix_grd(itr)=flasks(iflask)%mix_grd(itr)+this_weight*rmf !*fscale(offsetj+itr)

          enddo

          ! sample meteo
          if((myid == root) .and. flask_sample_meteo) then

             rmf =       wcx  *      wcy  * blh(ifr,jfr)  + &
                    (1.0-wcx) *      wcy  * blh(ifn,jfr) + &
                         wcx  * (1.0-wcy) * blh(ifr,jfn) + &
                    (1.0-wcx) * (1.0-wcy) * blh(ifn,jfn)

             flasks(iflask)%blh=flasks(iflask)%blh + this_weight*rmf


             rmf = &
                       wcx  *      wcy  *      wcz  * T(ifr,jfr,lfr)  + & 
                  (1.0-wcx) *      wcy  *      wcz  * T(ifn,jfr,lfr)  + &
                       wcx  * (1.0-wcy) *      wcz  * T(ifr,jfn,lfr)  + &
                  (1.0-wcx) * (1.0-wcy) *      wcz  * T(ifn,jfn,lfr)  + &
                       wcx  *      wcy  * (1.0-wcz) * T(ifr,jfr,lfrn) + & 
                  (1.0-wcx) *      wcy  * (1.0-wcz) * T(ifn,jfr,lfrn) + &
                       wcx  * (1.0-wcy) * (1.0-wcz) * T(ifr,jfn,lfrn) + &
                  (1.0-wcx) * (1.0-wcy) * (1.0-wcz) * T(ifn,jfn,lfrn) 

             flasks(iflask)%temperature=flasks(iflask)%temperature + this_weight*rmf

             rmf = &
                       wcx  *      wcy  *      wcz  * q(ifr,jfr,lfr)  + & 
                  (1.0-wcx) *      wcy  *      wcz  * q(ifn,jfr,lfr)  + &
                       wcx  * (1.0-wcy) *      wcz  * q(ifr,jfn,lfr)  + &
                  (1.0-wcx) * (1.0-wcy) *      wcz  * q(ifn,jfn,lfr)  + &
                       wcx  *      wcy  * (1.0-wcz) * q(ifr,jfr,lfrn) + & 
                  (1.0-wcx) *      wcy  * (1.0-wcz) * q(ifn,jfr,lfrn) + &
                       wcx  * (1.0-wcy) * (1.0-wcz) * q(ifr,jfn,lfrn) + &
                  (1.0-wcx) * (1.0-wcy) * (1.0-wcz) * q(ifn,jfn,lfrn) 

             flasks(iflask)%q=flasks(iflask)%q + this_weight*rmf

             yres = dy/yref(region)
             xres = dx/xref(region)

             do j = jfr,jfn
                lat(j-jfr+1) = ybeg(region) + 0.5 * yres + yres * (j-1)
             enddo

             ! Convert mass fluxes pu and pv to winds, following code
             ! in diffusion.F90.  -Andy Jacobson 4 Oct 12

             dxx=0.
             u=0.
             v=0.
             lat=0.

             dxx(:) = ae * xres * gtor * cos(lat(:)*gtor)
             dyy = ae * yres * gtor

             l0=min(lfr,lfrn)
             l1=max(lfr,lfrn)
             j0=min(jfr,jfn)
             j1=max(jfr,jfn)
             i0=min(ifr,ifn)
             i1=max(ifr,ifn)
             do l=l0,l1
                do j=j0,j1
                   do i=i0,i1
                      u(i+1-i0,j+1-j0,l+1-l0) = dxx(j+1-j0)*(pu(i,j,l) + pu(i-1,j,l))*0.5 / m(i,j,l)
                      v(i+1-i0,j+1-j0,l+1-l0) = dyy*  (pv(i,j,l) + pv(i,j+1,l))*0.5 / m(i,j,l)
                   enddo
                enddo
             enddo
             
             rmf = & 
                  (((0.5-rlf) * u(1,1,1) + (0.5+rlf) * u(1,1,2)) *      wcx  *      wcy  + &
                   ((0.5-rlf) * u(2,1,1) + (0.5+rlf) * u(2,1,2)) * (1.0-wcx) *      wcy  + &
                   ((0.5-rlf) * u(1,2,1) + (0.5+rlf) * u(1,2,2)) *      wcx  * (1.0-wcy) + &
                   ((0.5-rlf) * u(2,2,1) + (0.5+rlf) * u(2,2,2)) * (1.0-wcx) * (1.0-wcy))

             flasks(iflask)%u=flasks(iflask)%u + rmf

             rmf = & 
                  (((0.5-rlf) * v(1,1,1) + (0.5+rlf) * v(1,1,2)) *      wcx  *      wcy  + &
                   ((0.5-rlf) * v(2,1,1) + (0.5+rlf) * v(2,1,2)) * (1.0-wcx) *      wcy  + &
                   ((0.5-rlf) * v(1,2,1) + (0.5+rlf) * v(1,2,2)) *      wcx  * (1.0-wcy) + &
                   ((0.5-rlf) * v(2,2,1) + (0.5+rlf) * v(2,2,2)) * (1.0-wcx) * (1.0-wcy))

             flasks(iflask)%v=flasks(iflask)%v + rmf

             rmf = & 
                  (((0.5-rlf) * phlb(ifr,jfr,lfr) + (0.5+rlf) * phlb(ifr,jfr,lfrn)) *      wcx  *      wcy  + &
                   ((0.5-rlf) * phlb(ifn,jfr,lfr) + (0.5+rlf) * phlb(ifn,jfr,lfrn)) * (1.0-wcx) *      wcy  + &
                   ((0.5-rlf) * phlb(ifr,jfn,lfr) + (0.5+rlf) * phlb(ifr,jfn,lfrn)) *      wcx  * (1.0-wcy) + &
                   ((0.5-rlf) * phlb(ifn,jfn,lfr) + (0.5+rlf) * phlb(ifn,jfn,lfrn)) * (1.0-wcx) * (1.0-wcy))

             flasks(iflask)%pressure=flasks(iflask)%pressure + this_weight*rmf


          endif

          flasks(iflask)%nsamples=flasks(iflask)%nsamples+1

          flasks(iflask)%accum_weight=flasks(iflask)%accum_weight+this_weight

       end if
    end do

    nullify(m)
    nullify(rm)
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
    nullify(gph)
    nullify(T)
    nullify(blh)
    nullify(pu)
    nullify(pv)
    nullify(phlb)

  end subroutine user_output_flask_sample

!BOP
!
! !IROUTINE:  user_output_flask_evaluate
!
! !INPUT PARAMETERS: none
!
! !OUTPUT PARAMETERS: none
!
! !DESCRIPTION: 
!    - convert accumulated model mixing ratios to averages
!
!EOP

   subroutine user_output_flask_evaluate

    use GO,			only : gol, goErr, goPr, goBug
    use datetime,		only : date2tau, tau2date
    use dims,			only : jm

    ! average flask samples over nsamples

    integer 					:: itrace
    integer 					:: iflask
    character(len=15)				:: idstr
    integer,dimension(6)			:: idatef

    ! --- const ------------------------------

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

    if (nflasks .eq. 0) then
       return
    end if

    do iflask=1,nflasks
       if((.not. flasks(iflask)%evaluated) .and. (flasks(iflask)%accum_weight .gt. 0.)) then
          do itrace=1,ntracetloc
             flasks(iflask)%mix(itrace)=flasks(iflask)%mix(itrace)/flasks(iflask)%accum_weight
             flasks(iflask)%mix_grd(itrace)=flasks(iflask)%mix_grd(itrace)/flasks(iflask)%accum_weight
          enddo
          if(flasks(iflask)%region .eq. 1 .and. ((flasks(iflask)%jfr .eq. 1) .or. (flasks(iflask)%jfr .eq. jm(1)))) then
             write(idstr,'(i)') flasks(iflask)%obs_num
             call tau2date(flasks(iflask)%itau_center,idatef)
             write (gol,'("[user_output_flask_evaluate]  attention: flask ",a," (obs_id ",a,") at ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," at latitude ",f6.2," degrees; rejecting slopes sampling:")') &
                  trim(adjustl(idstr)), trim(adjustl(flasks(iflask)%obs_id)), idatef(1:5),flasks(iflask)%lat
             call goPr
             do itrace=1,ntracetloc
                write (gol,'("[user_output_flask_evaluate]     tracer",i4,"    slopes (not used): ",f6.2," ppm;    gridbox (used): ",f6.2," ppm.")') &
                     itrace,1.e6*flasks(iflask)%mix(itrace),1.e6*flasks(iflask)%mix_grd(itrace)  ! 1.e6 assumes CO2 here.
                call goPr

                ! This is where grid sampling is substituted for slopes sampling
                flasks(iflask)%mix(itrace) = flasks(iflask)%mix_grd(itrace)

             end do
          end if

          if((myid==root) .and. flask_sample_meteo) then
             flasks(iflask)%u=flasks(iflask)%u/flasks(iflask)%accum_weight
             flasks(iflask)%v=flasks(iflask)%v/flasks(iflask)%accum_weight
             flasks(iflask)%blh=flasks(iflask)%blh/flasks(iflask)%accum_weight
             flasks(iflask)%q=flasks(iflask)%q/flasks(iflask)%accum_weight
             flasks(iflask)%pressure=flasks(iflask)%pressure/flasks(iflask)%accum_weight
             flasks(iflask)%temperature=flasks(iflask)%temperature/flasks(iflask)%accum_weight
          endif
       else
          write(idstr,'(i8)') flasks(iflask)%obs_num
          call tau2date(flasks(iflask)%itau_center,idatef)
          write (gol,'("[user_output_flask_evaluate]  attention: flask ",a," (obs_id ",a,") at ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," not evaluated; nsamples is ",4i,".")') &
               trim(adjustl(idstr)), trim(adjustl(flasks(iflask)%obs_id)), idatef(1:5),flasks(iflask)%nsamples
          call goPr
          
          flasks(iflask)%mix=flask_missing_value
          flasks(iflask)%mix_grd=flask_missing_value
          if((myid==root) .and. flask_sample_meteo) then
             flasks(iflask)%u=flask_missing_value
             flasks(iflask)%v=flask_missing_value
             flasks(iflask)%blh=flask_missing_value
             flasks(iflask)%q=flask_missing_value
             flasks(iflask)%pressure=flask_missing_value
             flasks(iflask)%temperature=flask_missing_value
          endif
       endif
    enddo

    flasks%evaluated = .true.  ! prevents double evaluation

    
  end subroutine user_output_flask_evaluate

!BOP
!
! !IROUTINE:  user_output_flask_done
!
! !OUTPUT PARAMETERS: status integer
!
! !DESCRIPTION: 
!    - call evaluate routine
!    - gather results from other PEs
!    - write final mixing ratios to output file
!
!EOP

  subroutine user_output_flask_done(status)

    use GO,			only : gol, goErr, goPr, goBug
    use chem_param,		only : maxtrace, names, tracer_name_len, fscale
    use dims,			only : itaui, itaue, region_name
    use datetime,		only : tau2date

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

    integer, intent(out)			::  status

    ! --- const ------------------------------

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

    ! local
    integer					:: hid
    integer 					:: dim_obs, dim_ntracet, dim_tracer_name_len
    integer					:: dim_region_name_len,dim_indices
    integer 					:: var_obs_num, var_ntracet !, var_tracer_name_len
    integer 					:: var_tracer_names, var_flask
    integer 					:: var_nsamples, var_avetime, var_surface_height
    integer 					:: var_u,var_v,var_blh,var_q,var_pressure,var_temperature
    integer 					:: var_region_name,var_region_indices
    real, dimension(:,:), allocatable		:: mix_all
    real, dimension(:,:), allocatable 		:: mix_local
    integer					:: iflask, ipe, itrace
    integer, dimension(:), allocatable		:: recvcounts
    integer, dimension(:), allocatable		:: offsets
    integer, dimension(:), allocatable		:: nsamples
    character(len=1024)				:: attstring
    integer, dimension(tracer_name_len)		:: dimvar_name_len
    integer, dimension(ntracet)			:: dimvar_ntracet
    integer, dimension(6)			:: idatei,idatee,idatef
    real, dimension(:), allocatable		:: avetime

    call user_output_flask_evaluate

    call par_barrier

    if (nflasks .eq. 0 ) then
       return
    end if
    ! write results to output file

    do iflask = 1,tracer_name_len
       dimvar_name_len(iflask)=iflask
    enddo

    do iflask = 1,ntracet
       dimvar_ntracet(iflask)=iflask
    enddo

    ! recall that myid and root are published from module ParTools, and do 
    ! not depend on the MPI macro being defined.
    if (myid == root ) then

      ! new file:                                                                                                                                                                                                     
      call MDF_Create(trim(outFile), MDF_NETCDF4, MDF_REPLACE, hid, status )
      IF_NOTOK_RETURN(status=1)

      ! define dimensions:                                                                                                                                                                                            
      call MDF_Def_Dim(hid, 'obs_num', MDF_UNLIMITED, dim_obs, status )
      IF_NOTOK_RETURN(status=1)

      call MDF_Def_Dim(hid, 'tracer', ntracet, dim_ntracet, status )
      IF_NOTOK_RETURN(status=1)

      call MDF_Def_Dim(hid, 'tracer_name_len', tracer_name_len, dim_tracer_name_len, status )
      IF_NOTOK_RETURN(status=1)

      ! dimension variables:                                                                                                                                                                                                    
      call MDF_Def_Var(hid, 'obs_num', MDF_INT, (/dim_obs/), var_obs_num, status )
      IF_NOTOK_RETURN(status=1)

      call MDF_Def_Var(hid, 'tracer', MDF_INT, (/dim_ntracet/), var_ntracet, status )
      IF_NOTOK_RETURN(status=1)
      ! in dims_grid.F90, 10 chars is the region name length
      call MDF_Def_Dim(hid, 'region_name_len', 10, dim_region_name_len, status )
      IF_NOTOK_RETURN(status=1)

      call MDF_Def_Dim(hid, 'region_indices', 3, dim_indices, status )
      IF_NOTOK_RETURN(status=1)

      ! variables:                                                                                                                                                                                                    
      call MDF_Def_Var(hid, 'flask', MDF_FLOAT, (/dim_ntracet,dim_obs/), var_flask, status )
      IF_NOTOK_RETURN(status=1)

      call MDF_Def_Var(hid, 'nsamples', MDF_INT, (/dim_obs/), var_nsamples, status )
      IF_NOTOK_RETURN(status=1)

      call MDF_Def_Var(hid, 'tracer_names', MDF_CHAR, (/dim_tracer_name_len,dim_ntracet/), var_tracer_names, status )
      IF_NOTOK_RETURN(status=1)

      call MDF_Def_Var(hid, 'averaging_time', MDF_INT, (/dim_obs/), var_avetime, status )
      IF_NOTOK_RETURN(status=1)

      call MDF_Def_Var(hid, 'surface_height', MDF_FLOAT, (/dim_obs/), var_surface_height, status )
      IF_NOTOK_RETURN(status=1)

      call MDF_Def_Var(hid, 'region_name', MDF_CHAR, (/dim_region_name_len,dim_obs/), var_region_name, status )
      IF_NOTOK_RETURN(status=1)

      call MDF_Def_Var(hid, 'region_indices', MDF_INT, (/dim_indices,dim_obs/), var_region_indices, status )
      IF_NOTOK_RETURN(status=1)

      if(flask_sample_meteo) then
         call MDF_Def_Var( hid, 'u', MDF_FLOAT, (/dim_obs/), var_u, status )
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_u,"long_name",values="x-wind",status=status)
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_u,"units",values="m s^-1",status=status)
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_u,"_FillValue",values=flask_missing_value,status=status)
         IF_NOTOK_RETURN(status=1)

         call MDF_Def_Var( hid, 'v', MDF_FLOAT, (/dim_obs/), var_v, status )
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_v,"long_name",values="y-wind",status=status)
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_v,"units",values="m s^-1",status=status)
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_v,"_FillValue",values=flask_missing_value,status=status)
         IF_NOTOK_RETURN(status=1)

         call MDF_Def_Var( hid, 'blh', MDF_FLOAT, (/dim_obs/), var_blh, status )
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_blh,"long_name",values="atmosphere_boundary_layer_thickness",status=status)
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_blh,"units",values="m",status=status)
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_blh,"_FillValue",values=flask_missing_value,status=status)
         IF_NOTOK_RETURN(status=1)

         call MDF_Def_Var( hid, 'q', MDF_FLOAT, (/dim_obs/), var_q, status )
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_q,"long_name",values="mass_fraction_of_water_in_air",status=status)
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_q,"units",values="kg water (kg air)^-1",status=status)
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_q,"_FillValue",values=flask_missing_value,status=status)
         IF_NOTOK_RETURN(status=1)

         call MDF_Def_Var( hid, 'pressure', MDF_FLOAT, (/dim_obs/), var_pressure, status )
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_pressure,"long_name",values="air_pressure",status=status)
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_pressure,"units",values="Pa",status=status)
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_pressure,"_FillValue",values=flask_missing_value,status=status)
         IF_NOTOK_RETURN(status=1)

         call MDF_Def_Var( hid, 'temperature', MDF_FLOAT, (/dim_obs/), var_temperature, status )
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_temperature,"long_name",values="air_temperature",status=status)
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_temperature,"units",values="K",status=status)
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Att(hid, var_temperature,"_FillValue",values=flask_missing_value,status=status)
         IF_NOTOK_RETURN(status=1)

      endif

      ! add attributes
      call MDF_Put_Att(hid, var_flask,"long_name",values="mole_fraction_of_trace_gas_in_air",status=status)
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Att(hid, var_flask,"units",values="mol tracer (mol air)^-1",status=status)
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Att(hid, var_flask,"_FillValue",values=flask_missing_value,status=status)
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Att(hid, var_obs_num,"comment",values="Use this identifier to match observations with data from input file.",status=status)
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Att(hid, var_avetime,"units",values="seconds",status=status)
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Att(hid, var_avetime,"comment",values="Amount of model time over which this sample is averaged.",status=status)
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Att(hid, var_nsamples,"comment",values="Number of discrete samples in flask average.",status=status)
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Att(hid, var_surface_height,"comment",values="Height of the TM5 surface in the flask gridbox",status=status)
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Att(hid, var_surface_height,"units",values="meter",status=status)
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Att(hid, var_region_name,"comment",values="Name of TM5 zoom region containing sample.",status=status)
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Att(hid, var_region_name,"units",values="unitless",status=status)
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Att(hid, var_region_indices,"comment",values="Zonal, meridional, and level indices within the model region for grid cell containing sample.",status=status)
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Att(hid, var_region_indices,"units",values="unitless",status=status)
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Att(hid, MDF_GLOBAL,"input_file",values=inFile,status=status)
      IF_NOTOK_RETURN(status=1)
      
      call tau2date(itaui,idatei)
      call tau2date(itaue,idatee)

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

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

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

      ! mix_all is root's array to receive all data.  Note that
      ! ntracet is ntracetloc * npes, so this is implicitly
      ! dimensioned by no. PEs.  To be completely general, the output
      ! tracer array could be 3-dimensional:  no. types of tracers *
      ! no. flasks * no. ensemble members.  In practice, however,
      ! we're only interested in two of those dimensions.  For inverse
      ! runs, we have ensembles of total CO2; for forward runs we
      ! instead have one version each of CO2 components (e.g. CO2 due
      ! to air-sea exchange).  Both of these cases are 2-D, ntracet
      ! by nflasks.  If you want ensembles of the components, this
      ! code will still work, but the two dimensions (components and
      ! ensemble) will need to be unwrapped by hand by code reading
      ! the flask output file.

      allocate(mix_all(ntracet,nflasks))


      ! N.B.  every PE has a flasks structure, and for nsamples,
      ! itaui, and itaue, root's values should be the same as
      ! all other PEs, so there's no need for MPI communication.

      allocate(avetime(nflasks))
      avetime=flasks(:)%itau_end-flasks(:)%itau_start

      allocate(nsamples(nflasks))
      nsamples=flasks(:)%nsamples



   endif ! root only

   allocate(mix_local(ntracetloc,nflasks)) ! all PEs including root allocate a send buffer...

   ! ...and fill it
   do iflask=1,nflasks
      mix_local(:,iflask)=flasks(iflask)%mix
   enddo

#ifdef MPI

   if (myid .eq. root) then
      allocate(recvcounts(npes))
      allocate(offsets(npes))
      offsets(1)=0
      do ipe = 1, npes
         recvcounts(ipe)=ntracet_ar(ipe-1)
         if (ipe .gt. 1) then
            offsets(ipe)=sum(recvcounts(1:(ipe-1)))
         endif
      enddo
   else

      ! ifort's "-check all" complains severely if the mix_all, offsets, and recvcounts buffers
      ! are unallocated, despite the MPI spec's insistence that these buffers are only used
      ! by the root process.
 
     ! allocate(mix_all(1,1))
     ! allocate(recvcounts(1))
     ! allocate(offsets(1))

      allocate(mix_all(ntracet,nflasks))
      allocate(recvcounts(npes))
      allocate(offsets(npes))




   endif

   do iflask=1,nflasks
      ! MPI_GATHERV is how root gets all the mix_local values from other
      ! PEs.  Note that they are in rank-order, i.e. sorted by the myid
      ! value, and they include the values from the root process also.
      ! 
      ! Note that we use MPI_GATHERV (v for variable data amount for each
      ! node) because PEs can have differing numbers of tracers.
      communicator = com_trac
      call MPI_GATHERV(&
           mix_local(:,iflask), ntracetloc, MY_REAL, &
           mix_all(:,iflask), recvcounts, offsets, MY_REAL, & 
           root, communicator, ierr)
      if (.not. ierr .eq. MPI_SUCCESS) then
      endif
    !print*,iflask,mix_local(:,iflask)  
   enddo

#else

   ! not MPI: root needs to copy its mix_local values to mix_all
   mix_all = mix_local

#endif

   if (myid == root ) then

      call MDF_Put_Var(hid, var_obs_num, flasks%obs_num, status)
      IF_NOTOK_RETURN(status=1)

      call MDF_Put_Var(hid, var_ntracet, dimvar_ntracet, status)
      IF_NOTOK_RETURN(status=1)

!      call MDF_Put_Var(hid, var_tracer_name_len, dimvar_name_len, status)
!      IF_NOTOK_RETURN(status=1)

      call MDF_Put_Var(hid, var_flask, mix_all, status)
      IF_NOTOK_RETURN(status=1)

      call MDF_Put_Var(hid, var_nsamples, nsamples, status)
      IF_NOTOK_RETURN(status=1)

      call MDF_Put_Var(hid, var_avetime, avetime, status)
      IF_NOTOK_RETURN(status=1)

      call MDF_Put_Var(hid, var_surface_height, flasks%surface_height, status)
      IF_NOTOK_RETURN(status=1)
 
     
      call MDF_Put_Var(hid, var_region_indices, flasks%ifr, status,start=(/1,1/),count=(/1,nflasks/))
      IF_NOTOK_RETURN(status=1)

      call MDF_Put_Var(hid, var_region_indices, flasks%jfr, status,start=(/2,1/),count=(/1,nflasks/))
      IF_NOTOK_RETURN(status=1)

      call MDF_Put_Var(hid, var_region_indices, flasks%lfr, status,start=(/3,1/),count=(/1,nflasks/))
      IF_NOTOK_RETURN(status=1)

      do iflask=1,nflasks
         call MDF_Put_Var(hid, var_region_name, region_name(flasks(iflask)%region),status,start=(/1,iflask/),count=(/10,1/))
         IF_NOTOK_RETURN(status=1)
      enddo

      call MDF_Put_Var(hid, var_tracer_names, names(1:ntracet), status)
      IF_NOTOK_RETURN(status=1)

      if(flask_sample_meteo) then
         call MDF_Put_Var(hid, var_u, flasks%u, status )
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Var(hid, var_v, flasks%v, status )
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Var(hid, var_blh, flasks%blh, status )
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Var(hid, var_q, flasks%q, status )
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Var(hid, var_pressure, flasks%pressure, status )
         IF_NOTOK_RETURN(status=1)
         call MDF_Put_Var(hid, var_temperature, flasks%temperature, status )
         IF_NOTOK_RETURN(status=1)
      endif

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

   endif

   write (gol,'("[user_output_flask_done]  Deallocating arrays and closing output file.")'); call goPr

   do iflask = 1,nflasks
      if(allocated(flasks(iflask)%mix)) then
         deallocate(flasks(iflask)%mix)
      endif
      if(allocated(flasks(iflask)%mix_grd)) then
         deallocate(flasks(iflask)%mix_grd)
      endif
   enddo
   
   if(allocated(flasks)) then
      deallocate(flasks)
   endif

   if(allocated(mix_local)) then
      deallocate(mix_local)
   endif

   if(allocated(mix_all)) then
      deallocate(mix_all)
   endif

   if(allocated(avetime)) then
      deallocate(avetime)
   endif

   if(allocated(nsamples)) then
      deallocate(nsamples)
   endif

   if(allocated(recvcounts)) then
      deallocate(recvcounts)
   endif

   if(allocated(offsets)) then
      deallocate(offsets)
   endif

 end subroutine user_output_flask_done

end module user_output_flask

