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

#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  ! *************** DUMMY TM5-MP VERSION *******************
!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:
!
!     "obspack_id" : record dimension, unique 100-character string across obspacks
!     "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.
!  
!   Andy Jacobson, April 2013
!  
!   Change to ndyn weighting.
!
!   Mike Trudeau, Apr 2013
!
!   Removed 'obspack_num' references and added netcdf writes of 'obspack_id'.
!
!EOP
!-------------------------------------------------------------------------


! #ifdef MPI
!   use mpi_const,			only : my_real,mpi_comm_world,ierr,mpi_integer,mpi_character
!   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_NETCDF, MDF_REPLACE, MDF_GLOBAL, MDF_CHAR, MDF_INT, MDF_FLOAT, MDF_UNLIMITED, MDF_READ
!   use MDF      , 			only : MDF_Put_Att
!   use MDF      , 			only : MDF_Def_Dim
!   use MDF      , 			only : MDF_Inq_Dimid, MDF_Inquire_Dimension
!   use MDF      , 			only : MDF_Inq_VarID
!   use MDF      , 			only : MDF_Def_Var, MDF_Put_Var, MDF_Get_Var

  
  implicit none

  private

  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'

!   ! "obspack_id" is meant to be a flask identifier guaranteed to be unique, for the 
!   ! purpose of matching input file records with output file records.

!   type flask_sample
!      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=100)			:: obspack_id ! unique sample identifier
!      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, goTranslate

!     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		:: 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*100, dimension(:), allocatable	:: obspack_id
!     integer, dimension(6)			:: idatei,idatee,idatef
!     logical					:: new_region
!     integer					:: region
!     integer					:: iflask
!     integer					:: varid
!     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.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 (inFile,'(a,"/flask_input.",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 (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]  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, obspack_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_NETCDF, mode=MDF_READ, hid=hid, status=status)
!           IF_NOTOK_RETURN(status=1)

!           call MDF_Inq_DimID(hid, "obs", dimid, status)
!           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(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(obspack_id(nflasks))

!           call MDF_Inq_VarID(hid, "time_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, "obspack_id", varid, status)
!           IF_NOTOK_RETURN(status=1)
          
!           call MDF_Get_Var(hid, varid, obspack_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 goTranslate(obspack_id(iflask),char(0),' ',status)
!              IF_NOTOK_RETURN(status=1)

!              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 obspack_id string ",a,":")') trim(adjustl(obspack_id(iflask)))
!                 call goPr
!                 write (gol, '("  Unknown sampling strategy = ",i,".")') 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	 
!     call MPI_BCAST(nflasks, 1, MPI_INTEGER, root, MPI_COMM_WORLD, ierr)  
! #endif

!     if(nflasks .eq. 0) then
!        if(allocated(idate_f)) deallocate(idate_f)
!        if(allocated(mask)) deallocate(mask)
!        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(obspack_id)) deallocate(obspack_id)
!        if(allocated(flasks)) deallocate(flasks)
!        status=0
!        return
!     endif

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

!     allocate(flasks(nflasks))

!     do iflask=1,nflasks
!        allocate(flasks(iflask)%mix(ntracetloc))
!        allocate(flasks(iflask)%mix_grd(ntracetloc))
!        flasks(iflask)%obspack_id = " "
!        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%lat=pack(lat,mask)
!        flasks%lon=pack(lon,mask)
!        flasks%alt=pack(alt,mask)
!        flasks%obspack_id=pack(obspack_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

!           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: obspack_id ",a,": ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," at start; moved forward one hour.")') &
!                      trim(adjustl(flasks(iflask)%obspack_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: obspack_id ",a,": ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," at end; moved back one hour.")') &
!                      trim(adjustl(flasks(iflask)%obspack_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: obspack_id ",a,": ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," close to start; sampling truncated.")') &
!                      trim(adjustl(flasks(iflask)%obspack_id)), idatef(1:5)
!              else
!                 write (gol,'("[user_output_flask_init]  attention: obspack_id ",a,": ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," close to start; sampling continued from previous run.")') &
!                      trim(adjustl(flasks(iflask)%obspack_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: obspack_id ",a,": ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," close to end; sampling truncated.")') &
!                      trim(adjustl(flasks(iflask)%obspack_id)), idatef(1:5)
!              else
!                 write (gol,'("[user_output_flask_init]  attention: obspack_id ",a,": ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," close to end; sampling will continue in subsequent run.")') &
!                      trim(adjustl(flasks(iflask)%obspack_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)%obspack_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 region  longitude (i)  latitude (j)  altitude         date               obspack_id")');call goPr
!           do iflask = 1,nflasks
!              if(flasks(iflask)%region > -1) then
!                 call tau2date(flasks(iflask)%itau_center,idatef)
!                 write(gol,'("   ",i5," ",a," ",f8.2," (",i3,")",f8.2," (",i3,")",f9.1,"   ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2,":",i2.2, " UTC   ",a)') &
!                      iflask,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)%obspack_id))
!                 call goPr
!              endif
!           enddo
!        endif
       
!     endif ! myid is root

!     call par_barrier

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

!     call Par_Barrier

!     if(allocated(idate_f)) deallocate(idate_f)
!     if(allocated(mask)) deallocate(mask)
!     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(obspack_id)) deallocate(obspack_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,jm,ndyn,ndyn_max
!     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
!     use binas, only 					: ae

!     implicit none

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

     STATUS=0
     
!     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 obspack_id ",a,":")') trim(flasks(iflask)%obspack_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 accumulated weight

!     integer 					:: itrace
!     integer 					:: iflask
!     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
!              call tau2date(flasks(iflask)%itau_center,idatef)
!              write (gol,'("[user_output_flask_evaluate]  attention: obspack_id ",a," at ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," at latitude ",f6.2," degrees; rejecting slopes sampling:")') &
!                   trim(adjustl(flasks(iflask)%obspack_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
!           call tau2date(flasks(iflask)%itau_center,idatef)
!           write (gol,'("[user_output_flask_evaluate]  attention: obspack_id ",a," at ",i4.4,"/",i2.2,"/",i2.2," ",i2.2,":",i2.2," not evaluated; nsamples is ",4i,".")') &
!                trim(adjustl(flasks(iflask)%obspack_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, goTranslate
!     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


     status=0

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

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

!     ! local
!     integer					:: hid
!     integer 					:: dim_obs, dim_ntracet, dim_tracer_name_len, dim_char100
!     integer					:: dim_region_name_len,dim_indices
!     integer 					:: var_obspack_id, 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_NETCDF, MDF_REPLACE, hid, status )
!       IF_NOTOK_RETURN(status=1)

!       ! define dimensions:                                                                                                                                                                                            
!       call MDF_Def_Dim(hid, 'obs', 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)

!       call MDF_Def_Dim( hid, "char100", 100, dim_char100, 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, 'grid_indices', 3, dim_indices, status )
!       IF_NOTOK_RETURN(status=1)

!       ! dimension variables: 
!       call MDF_Def_Var(hid, 'obspack_id', MDF_CHAR, (/dim_char100, dim_obs/), var_obspack_id, status )
!       IF_NOTOK_RETURN(status=1)

! !      call MDF_Def_Var(hid, 'tracer', MDF_INT, (/dim_ntracet/), var_ntracet, status )
! !      IF_NOTOK_RETURN(status=1)

! !      call MDF_Def_Var(hid, 'tracer_name_len', MDF_INT, (/dim_tracer_name_len/), var_tracer_name_len, 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_obspack_id,"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))

!    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.
!       call MPI_GATHERV(&
!            mix_local(:,iflask), ntracetloc, MY_REAL, &
!            mix_all(:,iflask), recvcounts, offsets, MY_REAL, & 
!            root, MPI_COMM_WORLD, ierr)
!       if (.not. ierr .eq. MPI_SUCCESS) then
!       endif
!    enddo

! #else

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

! #endif

!    if (myid == root ) then

!       do iflask=1,nflasks
!          call goTranslate( flasks(iflask)%obspack_id,' ',char(0),status)
!          IF_NOTOK_RETURN(status=1)
!       enddo

!       call MDF_Put_Var(hid, var_obspack_id, (/(flasks(iflask)%obspack_id, iflask=1, nflasks)/), 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