ecearth3/sources/tm5mp/proj/output/user_output_station.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"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      USER_OUTPUT_STATION
!
! !DESCRIPTION: Routines to deal with output for specific stations.
!
!               This version adds code to the original version in order to
!               write AOD and alpha to the station files (cf #ifdef with_optics)
!\\
!\\
! !INTERFACE: 
!
MODULE USER_OUTPUT_STATION
  !
  ! !USES:
  !
  use GO,         only : gol, goErr, goPr
  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, itaur, xcyc
  use chem_param, only : ntrace, ntracet, names, fscale
  use file_hdf,   only : THdfFile, Init, Done, TSds, SetDim
  use ParTools,   only : isRoot, myid

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: read_stationlist               ! read files with station definitions
  public :: init_station_output            ! init output files

  public :: reset_stationconc_accumulator  ! set to 0 acummulation arrays
  public :: output_stationconc             ! accumulate arrays
  public :: evaluate_stationconc           ! compute mean and std dev
  public :: write_stationconc              ! write mean conc to file
  
  public :: free_stationfields             ! close output files, free arrays

#if defined(with_pm) || defined(with_optics)
  public :: n_stat2
  public :: io_hdf
  ! Optics wants to write wavelength information here.
  ! Optics can do it more easily than user_output.
  ! Optics uses user_output_station. This means that user_output_station 
  ! cannot use optics.
  ! PM also uses this.
#endif
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=100)           :: StationlistFFilename2  ! list of 'validation' stations
  character(len=100)           :: InputDir = './'
  logical                      :: output_stations_meteo = .true.
  character(len=300)           :: StationSpeclist 
  integer                      :: StationlistFileFormat  ! 1: old format, 2: new format, 3: newer format 

  integer                      :: n_stat2           ! number of 'validation' stations (not used for 4DVAR)

  integer, dimension(ntrace)   :: station_speclist  ! list of species for output   (may differ per PE)
  integer                      :: station_nspec     ! number of species for output (may differ per PE)

  integer, parameter           :: n_metdat = 7      ! max number of meteo data (for output in station output file)
  integer, parameter           :: n_interpol = 3    ! number of station interpolation methods
  integer, parameter           :: n_neighbors = 6   ! number of neighboring grid cells used for estimate of represent. error

  character(len=4), dimension(n_interpol), parameter :: c_interpol = (/'_grd','_slp','_int'/)
  character(len=12),dimension(n_metdat  ), parameter :: C_metdat   =  &
       (/'Temperature', &
       'Pressure   ',&
       'uwind      ',&
       'vwind      ',&
       'windspeed  ',&
       'winddir    ',&
       'blh        '/)

  character(len=12),dimension(n_metdat  ), parameter :: U_metdat   =  &
       (/'K           ', &
       'hPa         ', &
       'm/s         ', &
       'm/s         ', &
       'm/s         ', &
       'degrees wrtN', &
       'm           '/)

  type(THdfFile)               :: io_hdf          ! contains HDF output file
  integer                      :: io_record = 0
  ! 'open' datasets of the io_hdf output file:
  type(TSds), dimension(:,:), allocatable :: Sds
  type(TSds), dimension(n_metdat)         :: Sds_meteo
  type(TSds)                              :: Sds_idate

#ifdef with_optics
  Type(TSds), public :: sds_station_aod
  Type(TSds), public :: sds_station_alpha
#endif
#ifdef with_pm
  Type(TSds), public :: sds_station_pm
#endif

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

  character(len=*), parameter  ::  mname = 'user_output_station'
  !
  ! !PUBLIC TYPES:
  !
  type, public :: stations
     character(len=13)          :: ident     ! station identifier
     character(len=50)          :: name      ! station name
     real                       :: lat       ! latitude of station
     real                       :: lon       ! longitude of station
     real                       :: height    ! height of station
     real                       :: height_sample   ! sampling height
     real                       :: height_surface  ! height of model surface
     real                       :: height_sample_above_surf  ! height_sample-height_surface
     integer                    :: index     ! index of station
     integer                    :: region 
     integer                    :: is
     integer                    :: js
     integer                    :: ls
     character(len=2)           :: type      ! sampling type (FM = flask; CM = continuous)
     integer                    :: land      ! land or ocean site
     integer                    :: nonsurf   ! 1: nonsurface sites (e.g. towers); 0: 'regular' sites
     real,dimension(n_metdat)   :: metdat    ! meteodata (interpolated) at stations

     real,dimension(:,:),pointer :: c_mean   ! c_mean(species, interpolation)
     real,dimension(:,:),pointer :: c_std    ! c_std(species, interpolation)
     integer                     :: ncount   ! counter
     logical                     :: ontile   ! is located on current tile of domain decomposition

     ! In the chem version, the following comment was found:
     ! "No one is estimating representativeness errors. You can compare the
     ! different interpolations to acquire something similar." And the code
     ! with *_nb was commented.
     ! That code is kept live here in the user_output version until further decision.
          
     ! grid neighbors (for estimate of representativeness error)
     real,dimension(:,:),pointer :: c_mean_nb           ! mean concentrations of neighbors (species, neighb. index)
     real,dimension(:,:),pointer :: c_std_nb            ! standard deviation (species, neighb. index)
     integer,dimension(:,:),pointer :: ncount_nb        ! counter (species, neighb. index)

      real,dimension(:),pointer   :: dc_mod              ! estimate of representativeness error (species)
                                                         !(based on 3D gradient to grid cell neighbors)

   end type stations
  !
  ! !PUBLIC DATA MEMBERS:
  !
  type(stations), dimension(:), public, allocatable  :: stat  ! collection of stations
  character(len=100), public   :: StationlistFFilename        ! list of 'regular' stations
  logical,            public   :: write_Stationfiles = .true. ! do station 
  integer, save,      public   :: n_stat         ! number of 'regular' stations (used for 4DVAR) 
  integer, save,      public   :: io_nrecords    ! # record = # timesteps
  !
  ! !REVISION HISTORY:
  ! 
  ! Modified by Peter Beramaschi.
  ! Parallel version Maarten Krol Jul 2003
  ! jun 2005: switch to HDF output file
  ! jul 2007: paralel writing to increase performance..
  !   10 Jul 2012 - Ph. Le Sager - merged with chem base version to account
  !                                for with_optics, with_pm 
  !                              - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    FREE_STATIONFIELDS
  !
  ! !DESCRIPTION: close output files and deallocate arrays
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE FREE_STATIONFIELDS( status )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(inout)   :: status
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer :: i_stat

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

    !__START_SUBROUTINE______________________________________________________

    if(write_Stationfiles) then 
       call close_stationfiles( status )
       IF_NOTOK_RETURN(status=1)
    endif

    write(gol,*) '************************************', station_nspec ; call goPr

    !cmk : the following lines do not work on huygens????

    if (station_nspec > 0 ) then 
       do i_stat = 1, n_stat+n_stat2

          deallocate(stat(i_stat)%c_mean)
          deallocate(stat(i_stat)%c_std )

       enddo

       do i_stat = 1, n_stat
          deallocate(stat(i_stat)%c_mean_nb)
          deallocate(stat(i_stat)%c_std_nb)
          deallocate(stat(i_stat)%ncount_nb)
          deallocate(stat(i_stat)%dc_mod)
       enddo

       deallocate( stat )
    end if
    
    status = 0

  END SUBROUTINE FREE_STATIONFIELDS
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    READ_STATIONLIST
  !
  ! !DESCRIPTION: Reads input files with stations definition. Allocate and
  !               fill station objects.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE READ_STATIONLIST( status )
    !
    ! !USES:
    !
    use GO,           only : TrcFile, Init, Done, ReadRc
    use global_data,  only : rcfile
    use chem_param,   only : ntracet, names
    use tm5_distgrid, only : dgrid, get_distgrid
    use ParTools,     only : isRoot
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(inout)            :: status
    !
    ! !REVISION HISTORY: 
    !   16 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer                 :: slen1, slen2 !string length
    character(len=200)      :: dummy
    integer                 :: fstatus 
    integer                 :: i_stat
    integer                 :: region
    integer, dimension(6)   :: idate_sim
    integer                 :: itau_sim, isim, ispec
    integer                 :: ipos, i, n
    integer                 :: i1, j1, i2, j2, is, js ! tile bounds and station indices
    real                    :: dxr, dyr   ! x,y resolution (in degrees) for current region 
    character(len=8)        :: name_spec  
    integer                 :: sunit1=110
    integer                 :: sunit2=111
                           
    type(TrcFile)           :: rcF

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

    !__START_SUBROUTINE______________________________________________________

    call Init( rcF, rcfile, status)
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'inputdir.station', InputDir,status)
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'input.station.filename', StationlistFFilename,status)
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'input.station.filename2', StationlistFFilename2,status)
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'input.station.fileformat', StationlistFileFormat,status)   
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'station.speclist', StationSpeclist,status)
    IF_NOTOK_RETURN(status=1)
    call Done( rcF ,status)
    IF_NOTOK_RETURN(status=1)

    ! -------- count # of regular and validation stations
    n_stat=0
    n_stat2=0

    slen1=len_trim(StationlistFFilename) 
    if ( slen1 > 0 ) then
       !=========================================
       ! open StationFile (with list of stations)
       !=========================================
       open(UNIT=sunit1, FORM='FORMATTED', &
            FILE=trim(InputDir)//trim(StationlistFFilename), &
            STATUS='OLD', ERR=1000)      
       read(sunit1, '(a)') dummy

       ! count number of stations
       do
          read(sunit1,  '(a)', end=100) dummy
          !if((dummy(1:1) /= ' ') .and. (dummy(1:7) /=  'species')) then   
          n_stat = n_stat+1 
       enddo
100    continue
    end if


    slen2=len_trim(StationlistFFilename2) 
    if ( slen2 > 0 ) then
       !=========================================
       ! open StationFile (with list of stations)
       !=========================================
       open(UNIT=sunit2, FORM='FORMATTED', &
            FILE=trim(InputDir)//trim(StationlistFFilename2), &
            STATUS='OLD', ERR=1001)      
       read(sunit2, '(a)') dummy

       ! count number of stations
       do
          read(sunit2,  '(a)', END=200) dummy
          if((dummy(1:1) /= ' ') .and. (dummy(1:7) /=  'species')) then   

             n_stat2 = n_stat2+1 
          else
             goto 200
          endif
       enddo

200    continue

    end if

    ! ---- determine index of species in good order
    station_nspec = 0
    ipos = 1
    do
       i = index(StationSpeclist(ipos:),' ')
       if (i == 0) then
          name_spec = trim(StationSpeclist(ipos:))
       else 
          name_spec = trim(StationSpeclist(ipos:ipos+i-1))
          ipos = ipos + i 
       endif
       do n=1,ntracet
          if ((trim(name_spec).eq.trim(names(n))).or.(trim(name_spec).eq.'all')) then
             station_nspec = station_nspec + 1
             station_speclist(station_nspec) = n
          endif
       enddo
       if(trim(name_spec).eq.'') exit
    enddo

    write(gol,*) trim(mname),':station_nspec; station_speclist', &
         station_nspec, station_speclist(1:station_nspec)
    call goPr


    ! ---- allocate station objects
    if (station_nspec > 0) then

       allocate(stat(n_stat+n_stat2)) ! all PE have all station, but only some will fill data
       
       do i_stat = 1, n_stat+n_stat2
          allocate(stat(i_stat)%c_mean(station_nspec,n_interpol))
          allocate(stat(i_stat)%c_std(station_nspec,n_interpol))
       enddo

       do i_stat = 1, n_stat
          allocate(stat(i_stat)%c_mean_nb(station_nspec,n_neighbors))
          allocate(stat(i_stat)%c_std_nb(station_nspec,n_neighbors))
          allocate(stat(i_stat)%ncount_nb(station_nspec,n_neighbors))
          allocate(stat(i_stat)%dc_mod(station_nspec))
       enddo

       
       ! ---- Read/fill stations definition, assuming first they are in region #1
       
       dyr = dy/yref(1)
       dxr = dx/xref(1)
       call get_distgrid( dgrid(1), i_strt=i1, i_stop=i2, j_strt=j1, j_stop=j2)

       if ( slen1 > 0 ) then

          !=========================
          ! read 'regular' stations
          !=========================
          rewind(sunit1)
          read(sunit1, '(a)') dummy

          do i_stat=1, n_stat 

             select case(StationlistFileFormat)
             case(1)
                read(sunit1,  '(a7,1x,f8.2,f8.2,f8.1,1x,a50)') &
                     stat(i_stat)%ident, &
                     stat(i_stat)%lat, &
                     stat(i_stat)%lon, &
                     stat(i_stat)%height, &
                     stat(i_stat)%name
                stat(i_stat)%nonsurf = 0
             case(2)
                read(sunit1,  '(a11,f8.2,4x,f8.2,4x,f8.1,1x,a2,1x,i2,1x,a50)') &
                     stat(i_stat)%ident, &
                     stat(i_stat)%lat, &
                     stat(i_stat)%lon, &
                     stat(i_stat)%height, &
                     stat(i_stat)%type, &
                     stat(i_stat)%nonsurf, &
                     stat(i_stat)%name
             case(3)
                read(sunit1,  '(a13,f6.2,4x,f8.2,4x,f8.0,4x,i1,5x,i1,a50)') &
                     stat(i_stat)%ident, &
                     stat(i_stat)%lat, &
                     stat(i_stat)%lon, &
                     stat(i_stat)%height, &
                     stat(i_stat)%land, &
                     stat(i_stat)%nonsurf, &
                     stat(i_stat)%name
                stat(i_stat)%ident = adjustl(stat(i_stat)%ident)
             endselect

             stat(i_stat)%index = i_stat

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

             ! is it on current PE?
             if (is>=i1.and.is<=i2.and.js>=j1.and.js<=j2) then
                stat(i_stat)%ontile=.true.
             else
                stat(i_stat)%ontile=.false.
             end if

          end do

          ! close stationfile
          close(sunit1)

       endif

       if ( slen2 > 0 ) then
          
          !=========================
          ! read 'validation' stations
          !=========================
          rewind(sunit2)
          read(sunit2, '(a)') dummy

          do i_stat=n_stat+1, n_stat+n_stat2 

             select case(StationlistFileFormat)
             case(1)
                read(sunit2,  '(a7,1x,f8.2,f8.2,f8.1,1x,a50)') &
                     stat(i_stat)%ident, &
                     stat(i_stat)%lat, &
                     stat(i_stat)%lon, &
                     stat(i_stat)%height, &
                     stat(i_stat)%name

                stat(i_stat)%nonsurf = 0
             case(2)
                read(sunit2,  '(a11,f8.2,f8.2,f8.1,1x,a2,1x,i2,1x,a50)') &
                     stat(i_stat)%ident, &
                     stat(i_stat)%lat, &
                     stat(i_stat)%lon, &
                     stat(i_stat)%height, &
                     stat(i_stat)%type, &
                     stat(i_stat)%nonsurf, &
                     stat(i_stat)%name
             case(3)
                read(sunit2,  '(a13,f6.2,4x,f8.2,4x,f8.0,4x,i1,5x,i1,a50)') &
                     stat(i_stat)%ident, &
                     stat(i_stat)%lat, &
                     stat(i_stat)%lon, &
                     stat(i_stat)%height, &
                     stat(i_stat)%land, &
                     stat(i_stat)%nonsurf, &
                     stat(i_stat)%name
                stat(i_stat)%ident = adjustl(stat(i_stat)%ident)
             endselect

             stat(i_stat)%index = i_stat

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

             if (is>=i1.and.is<=i2.and.js>=j1.and.js<=j2) then
                stat(i_stat)%ontile=.true.
             else
                stat(i_stat)%ontile=.false.
             end if


          end do

          ! close stationfile
          close(sunit2)

       end if

       
       !=================================================
       ! determine region with highest refinement factor
       !=================================================
       
       do i_stat=1, n_stat+n_stat2
          do region=1, nregions

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

             if ( (stat(i_stat)%lon .ge. xbeg(region) .and. &
                  stat(i_stat)%lon .le. xend(region)) .and.      &
                  (stat(i_stat)%lat .ge. ybeg(region) .and. &
                  stat(i_stat)%lat .le. yend(region) ) ) then
                !=====================
                ! station is in region
                !=====================

                ! determine region with hightes refinement factor
                if (xref(region) > xref(stat(i_stat)%region)) then
                   stat(i_stat)%region = region
                   is = int((stat(i_stat)%lon-float(xbeg(region)))/dxr + 0.99999)
                   js = int((stat(i_stat)%lat-float(ybeg(region)))/dyr + 0.99999)
                   stat(i_stat)%is = is
                   stat(i_stat)%js = js

                   call get_distgrid( dgrid(region), i_strt=i1, i_stop=i2, j_strt=j1, j_stop=j2)
                   if (is>=i1.and.is<=i2.and.js>=j1.and.js<=j2) then
                      stat(i_stat)%ontile=.true.
                   else
                      stat(i_stat)%ontile=.false.
                   end if

                end if
             end if
          end do
       end do

       if ( isRoot ) then
          write(gol,*)  '==========================' ; call goPr
          write(gol,*)  ' regular stations         ' ; call goPr
          write(gol,*)  '==========================' ; call goPr
          do i_stat=1, n_stat 
             write(gol, '(a13,1x,f8.2,f8.2,f8.1,1x,i4,1x,i2,1x,a50)') &
                  stat(i_stat)%ident, &
                  stat(i_stat)%lat, &
                  stat(i_stat)%lon, &
                  stat(i_stat)%height, &
                  stat(i_stat)%region, &
                  stat(i_stat)%nonsurf, &
                  stat(i_stat)%name ; call goPr
          end do

          write(gol,*)  '==========================' ; call goPr
          write(gol,*)  ' validation stations      ' ; call goPr
          write(gol,*)  '==========================' ; call goPr
          
          do i_stat=n_stat+1, n_stat+n_stat2
             write(gol, '(a11,1x,f8.2,f8.2,f8.1,1x,i4,1x,i2,1x,a50)') &
                  stat(i_stat)%ident, &
                  stat(i_stat)%lat, &
                  stat(i_stat)%lon, &
                  stat(i_stat)%height, &
                  stat(i_stat)%region, &
                  stat(i_stat)%nonsurf, &
                  stat(i_stat)%name ; call goPr
          end do
       endif

    endif
    call reset_stationconc_accumulator

    ! ok
    status = 0
    return

    ! error handling
1000 write(gol,*) 'read_stationlist: could not open ' // trim(InputDir)//trim(StationlistFFilename) ; call goErr
1001 write(gol,*) 'read_stationlist: could not open ' // trim(InputDir)//trim(StationlistFFilename2) ; call goErr
    status = 1

  END SUBROUTINE READ_STATIONLIST
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    INIT_STATION_OUTPUT
  !
  ! !DESCRIPTION: initialize hdf files for station output.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE INIT_STATION_OUTPUT( status )
    !
    ! !USES:
    !
    use Dims,        only : itaui, itaue, ndyn_max, idatei, idatee
    use file_hdf,    only : Init, WriteAttribute, WriteData, TSds, Done
    use GO,          only : TrcFile, Init, Done, ReadRc
    use global_data, only : rcfile, outdir
    use ParTools,    only : isRoot
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)      ::  status
    !
    ! !REVISION HISTORY: 
    !   16 Mar 2012 - P. Le Sager -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer                  :: i_stat
    integer                  :: io_status
    type(TrcFile)            :: rcF
    character(len=12)        :: xname
    integer                  :: i, ii, j
    type(TSds)               :: sds_temp
    character(len=200)       :: station_filename
    character(len=3)         :: smyid
    character(len=10)        :: sidatei, sidatee

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

    !__START_SUBROUTINE______________________________________________________

    if (write_Stationfiles .and. station_nspec > 0)then

       call Init( rcF, rcfile ,status)
       IF_NOTOK_RETURN(status=1)
       call ReadRc( rcF, 'output.station.filename', station_filename,status)
       IF_NOTOK_RETURN(status=1)
       call Done( rcF ,status)
       IF_NOTOK_RETURN(status=1)
#ifdef MPI 
       i = index(station_filename,'.hdf')
       write(smyid,'(a1,i2.2)') '_',myid
       write(sidatei,'(i4.4,3i2.2)') idatei(1:4)
       write(sidatee,'(i4.4,3i2.2)') idatee(1:4)
       if(i /= 0) then
          station_filename = station_filename(1:i-1)//smyid//'_'//sidatei//'_'//sidatee//'.hdf'
       else
          station_filename = trim(station_filename)//smyid//'_'//sidatei//'_'//sidatee//'.hdf'
       endif
#else
       i = index(station_filename,'.hdf')
       write(sidatei,'(i4.4,3i2.2)') idatei(1:4)
       write(sidatee,'(i4.4,3i2.2)') idatee(1:4)
       if(i /= 0) then
          station_filename = station_filename(1:i-1)//'_'//sidatei//'_'//sidatee//'.hdf'
       else
          station_filename = trim(station_filename)//'_'//sidatei//'_'//sidatee//'.hdf'
       endif
#endif
       call Init( io_hdf, trim(outdir)//'/'//trim(station_filename),'create', status )       
       IF_NOTOK_RETURN(status=1)
       io_record = 0
       Write(*,*) "Using itaui in stations"
       io_nrecords =abs(itaue-itaui)/ndyn_max
       call WriteAttribute( io_hdf, 'io_nrecords', io_nrecords,status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'n_stat', n_stat,status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'n_stat2', n_stat2,status)         ! validation stations
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'n_spec', station_nspec,status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'station_speclist', station_speclist(1:station_nspec), status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'station_specnames', names(station_speclist(1:station_nspec)), status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'station_lat', stat(:)%lat, status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'station_lon', stat(:)%lon, status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'station_height', stat(:)%height, status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'station_index', stat(:)%index, status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'station_nonsurf', stat(:)%nonsurf, status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'station_region', stat(:)%region, status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'station_is', stat(:)%is, status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'station_js', stat(:)%js, status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'station_names', stat(:)%name, status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'station_ident', stat(:)%ident, status)
       IF_NOTOK_RETURN(status=1)

       allocate(Sds(station_nspec, n_interpol))
       do i=1, station_nspec
          do ii= 1,n_interpol
             xname = trim(names(station_speclist(i)))//c_interpol(ii)
             call init(Sds(i,ii), io_hdf, xname, (/n_stat+n_stat2, io_nrecords/), 'real(4)', status)
             IF_NOTOK_RETURN(status=1)
             call SetDim( sds(i,ii), 0, 'n_stations', 'station number', (/ (j, j=1,n_stat+n_stat2) /) , status)
             IF_NOTOK_RETURN(status=1)
             call SetDim( sds(i,ii), 1, 'n_records ', 'time slot #', (/ (j, j=1,io_nrecords) /) , status)
             IF_NOTOK_RETURN(status=1)
             call WriteAttribute(Sds(i,ii), 'Unit', 'Volume mixing ration (ppb)', status)
             IF_NOTOK_RETURN(status=1)
          enddo
       enddo

       if (.not. isRoot) output_stations_meteo = .false.
       if (output_stations_meteo) then 
          do i=1,n_metdat
             call init(sds_meteo(i), io_hdf, c_metdat(i), (/n_stat+n_stat2, io_nrecords/), 'real(4)', status)
             IF_NOTOK_RETURN(status=1)
             call SetDim( sds_meteo(i), 0, 'n_stations', 'station number', (/ (j, j=1,n_stat+n_stat2) /) , status)
             IF_NOTOK_RETURN(status=1)
             call SetDim( sds_meteo(i), 1, 'n_records ', 'time slot #', (/ (j, j=1,io_nrecords) /) , status)
             IF_NOTOK_RETURN(status=1)
             call WriteAttribute(Sds_meteo(i), 'Unit', u_metdat(i), status)
             IF_NOTOK_RETURN(status=1)
          enddo
       endif
       call init(sds_idate, io_hdf, 'idate', (/6, io_nrecords/), 'integer(2)', status)
       IF_NOTOK_RETURN(status=1)
       call SetDim( sds_idate, 1, 'n_records ', 'time slot #', (/ (j, j=1,io_nrecords) /), status )
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(Sds_idate, 'Unit', 'year month day hour minutes seconds', status)
       IF_NOTOK_RETURN(status=1)
    end if   ! write stationfiles = T
    ! ok
    status = 0
    
  END SUBROUTINE INIT_STATION_OUTPUT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CLOSE_STATIONFILES
  !
  ! !DESCRIPTION: close ouput files.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CLOSE_STATIONFILES( status )
    !
    ! !USES:
    !
    use file_hdf,  only : Done
    use ParTools,  only : isRoot
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(inout)            :: status
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer                           :: i, ii
    character(len=*), parameter       :: rname = mname//'/close_station_files'

    !__START_SUBROUTINE______________________________________________________

    if (write_Stationfiles .and. station_nspec > 0)then

       do i=1, station_nspec
          do ii= 1,n_interpol
             call Done(sds(i,ii), status)
             IF_NOTOK_RETURN(status=1)
          enddo
       enddo
       
       if (.not.isRoot) output_stations_meteo = .false.

       if (output_stations_meteo) then 
          do i=1,n_metdat
             call Done(sds_meteo(i),status)
             IF_NOTOK_RETURN(status=1)
          enddo
       endif
       call Done(sds_idate, status)
       IF_NOTOK_RETURN(status=1)

       deallocate(Sds)
       
       call Done(io_hdf, status)
       IF_NOTOK_RETURN(status=1)
       write(gol,*) 'Station file closed ' ; call goPr
    endif

    ! ok
    status = 0

  END SUBROUTINE CLOSE_STATIONFILES
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OUTPUT_STATIONCONC
  !
  ! !DESCRIPTION: accumulate mean and std dev arrays
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE OUTPUT_STATIONCONC(region, status)
    !
    ! !USES:
    !
    use global_data, only : mass_dat
#ifndef without_convection
    use global_data, only : conv_dat
#endif
    use MeteoData,   only : m_dat, gph_dat, temper_dat, pu_dat, pv_dat, phlb_dat
    use datetime,    only : tau2date
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in)   :: region
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer,intent(inout):: status
    !
    ! !REVISION HISTORY: 
    !   16 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

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


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

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

    !=============
    ! meteo output
    !=============
    real,dimension(n_stat,n_metdat) :: metdat               ! meteodata (interpolated) at stations
    real,dimension(:,:,:), pointer  :: T                    ! temperature
    real,dimension(:,:,:), pointer  :: phlb                 ! pressure grid boundaries
    real,dimension(:,:,:), pointer  :: pu                   ! mass flux x-direction [kg/s]
    real,dimension(:,:,:), pointer  :: pv                   ! mass flux y-direction [kg/s]
    real,dimension(:,:), pointer    :: blh                  ! boundary layer height [m]
    real                            :: T_interpol           ! temperature, interpolated to station location 
    real                            :: p_interpol           ! pressure, interpolated to station location 
    real                            :: windspeed_u_interpol ! wind speed [m/s] x-direction , interpolated to station location
    real                            :: windspeed_v_interpol ! wind speed [m/s] y-direction , interpolated to station location
    real                            :: windspeed_interpol   ! wind speed [m/s]             , interpolated to station location
    real                            :: winddir_interpol     ! wind direction               , interpolated to station location 
    real                            :: blh_interpol         ! boundary layer height        , interpolated to station location 
    real                            :: lat_j                ! latitude of grid cell center js  [degrees]
    real                            :: lat_jn               ! latitude of grid cell center jsn [degrees]
    real                            :: ddx_j                ! x-extension of grid cell js    [m]
    real                            :: ddx_jn               ! x-extension of grid cell jsn   [m]
    real                            :: ddy                  ! y-extension of grid cell js    [m]
    real                            :: vmr                  ! volume mi

    integer   :: communicator,root_id,lmr,imr,jmr

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

    !__START_SUBROUTINE______________________________________________________

    if (station_nspec == 0) then
       status = 0
       return
    endif

    imr = im(region) ; jmr = jm(region) ; lmr = lm(region)

    m => m_dat(region)%data
    rm => mass_dat(region)%rm
    rxm => mass_dat(region)%rxm
    rym => mass_dat(region)%rym
    rzm => mass_dat(region)%rzm

    ! assign pointers (METEO)

    gph => gph_dat(region)%data
    t => temper_dat(region)%data
    pu => pu_dat(region)%data
    pv => pv_dat(region)%data
#ifndef without_convection
    blh => conv_dat(region)%blh
#endif
    phlb => phlb_dat(region)%data

    ! perform now on all PEs
    ! x and y resolution [degrees] for current region

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

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


    !====================
    ! loop over stations
    !====================

    do i_stat=1, n_stat+n_stat2

       if (stat(i_stat)%region == region) then

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

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

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

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


          !=================================================================
          ! if index of neighbour is exceeding range of region set 
          ! neighbour = current cell (i.e. no interpolation)
          ! in case of cyclic x-boundaries take corresponding cyclic i index
          !=================================================================

          if ( jsn == 0) jsn=1
          if ( jsn == jm(region)+1 ) jsn=jm(region)  ! isn-->jsn (wouter Peters)
!           if ( xcyc(region) == 0 ) then
!              ! non-cyclic boundaries
!              if ( isn == 0) isn=1
!              if ( isn == im(region)+1 ) isn=im(region)
!           else
!              ! cyclic x-boundaries
!              if ( isn == 0 ) isn=im(region)
!              if ( isn == im(region)+1 ) isn=1
!           end if

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

          !===========================
          ! determine layer of station
          !===========================
          if(stat(i_stat)%nonsurf == 1) then
             sampleheight=height(0)+stat(i_stat)%height 
          else
             sampleheight=stat(i_stat)%height 
          endif
          ! do not allow sampleheight below model surface
          sampleheight = max(sampleheight, height(0))

          do l=0,lm(region)
             if(height(l) > sampleheight) exit
          end do

          select case(l)
          case(0)
             !PB this should no longer occur !!
             ! force station to model surface
             ls = 1
             rls = -0.5  
          case default
             ls = l       ! station's layer
             ! the off-set from the center of the layer (-0.5--->+0.5) 
             ! (interpolation is in (m))
             rls = (sampleheight - height(l-1)) / &
                  (height(l)-height(l-1)) - 0.5
          end select
          stat(i_stat)%ls = ls  ! store for output to file
          stat(i_stat)%height_sample = sampleheight  ! store for output to file
          stat(i_stat)%height_surface = height(0)  ! store for output to file
          stat(i_stat)%height_sample_above_surf = sampleheight-height(0)  ! store for output to file


          !=================================
          !the neighbour for z interpolation
          !=================================
          if ( rls > 0 ) then 
             lsn = ls+1     
          else
             lsn = ls-1
          end if
          ! z-weighting of grid cell in which station is located 
          wcz = (1.0-abs(rls))  !.0 ... 0.5

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

          if(lsn == 0) then
             lsn=2
             wcz=1.0-rls  ! 1.0 ... 1.5
          endif
          if(lsn == lmr+1) then
             !=========================================================
             ! if vertical neighbor is lmr+1 (which does not exist)
             ! -> no interpolation
             !=========================================================
             lsn=lmr      ! no interpolation 
             wcz=1.0
          endif

          do ispec=1,station_nspec
                n = station_speclist(ispec)   
                lst = ls 
                lstn = lsn 
                !========================================
                ! value of grid box without interpolation
                !========================================
                vmr = rm(is,js,lst,n) / m(is,js,lst) *  fscale(n)
                stat(i_stat)%c_mean(ispec,1) = &
                     stat(i_stat)%c_mean(ispec,1) + vmr

                stat(i_stat)%c_std(ispec,1) = &
                     stat(i_stat)%c_std(ispec,1) + vmr*vmr


                if(i_stat <= n_stat) then

                   ! for 'regular' stations only

                   !=================================================
                   ! values of neigboring grid boxes 
                   ! (used as estimate for representativeness error
                   !=================================================

                   ! lower neighbor
                   if(ls > 1)then
                      vmr = rm(is,js,ls-1,n) / m(is,js,ls-1) * fscale(n)
                      stat(i_stat)%c_mean_nb(ispec,1) = stat(i_stat)%c_mean_nb(ispec,1) + vmr
                      stat(i_stat)%c_std_nb(ispec,1) = stat(i_stat)%c_std_nb(ispec,1) + vmr*vmr
                      stat(i_stat)%ncount_nb(ispec,1) = stat(i_stat)%ncount_nb(ispec,1) + 1
                   endif

                   ! upper neighbor
                   if(ls < lm(region)) then
                      vmr = rm(is,js,ls+1,n) / m(is,js,ls+1) * fscale(n)
                      stat(i_stat)%c_mean_nb(ispec,2) = stat(i_stat)%c_mean_nb(ispec,2) + vmr
                      stat(i_stat)%c_std_nb(ispec,2) = stat(i_stat)%c_std_nb(ispec,2) + vmr*vmr
                      stat(i_stat)%ncount_nb(ispec,2) = stat(i_stat)%ncount_nb(ispec,2) + 1
                   endif

                   ! western neighbor
                   if(is > 1) then
                      vmr = rm(is-1,js,ls,n) / m(is-1,js,ls) * fscale(n)
                      stat(i_stat)%c_mean_nb(ispec,3) = stat(i_stat)%c_mean_nb(ispec,3) + vmr
                      stat(i_stat)%c_std_nb(ispec,3) = stat(i_stat)%c_std_nb(ispec,3) + vmr*vmr
                      stat(i_stat)%ncount_nb(ispec,3) = stat(i_stat)%ncount_nb(ispec,3) + 1
                   endif

                   ! eastern neighbor
                   if(is < im(region)) then
                      vmr = rm(is+1,js,ls,n) / m(is+1,js,ls) * fscale(n)
                      stat(i_stat)%c_mean_nb(ispec,4) = stat(i_stat)%c_mean_nb(ispec,4) + vmr
                      stat(i_stat)%c_std_nb(ispec,4) = stat(i_stat)%c_std_nb(ispec,4) + vmr*vmr
                      stat(i_stat)%ncount_nb(ispec,4) = stat(i_stat)%ncount_nb(ispec,4) + 1
                   endif

                   ! southern neighbor
                   if(js > 1) then
                      vmr = rm(is,js-1,ls,n) / m(is,js-1,ls) * fscale(n)
                      stat(i_stat)%c_mean_nb(ispec,5) = stat(i_stat)%c_mean_nb(ispec,5) + vmr
                      stat(i_stat)%c_std_nb(ispec,5) = stat(i_stat)%c_std_nb(ispec,5) + vmr*vmr
                      stat(i_stat)%ncount_nb(ispec,5) = stat(i_stat)%ncount_nb(ispec,5) + 1
                   endif

                   ! northern neighbor
                   if(js < jm(region)) then
                      vmr = rm(is,js+1,ls,n) / m(is,js+1,ls) * fscale(n)
                      stat(i_stat)%c_mean_nb(ispec,6) = stat(i_stat)%c_mean_nb(ispec,6) + vmr
                      stat(i_stat)%c_std_nb(ispec,6) = stat(i_stat)%c_std_nb(ispec,6) + vmr*vmr
                      stat(i_stat)%ncount_nb(ispec,6) = stat(i_stat)%ncount_nb(ispec,6) + 1
                   endif

                end if
                


                
                !========================================
                ! interpolation using slopes  
                !========================================
                !rm-value is obtained from rm + slopes. 
                !slope = rxm = (rm*dX/dx *deltaX/2)
                rm_interpol = rm(is,js,lst,n) + 2.0*(ris*rxm(is,js,lst,n) + &
                     rjs*rym(is,js,lst,n) + &
                     rls*rzm(is,js,lst,n) )
                rm_interpol = max(0.0,rm_interpol) 
                vmr = rm_interpol/m(is,js,lst) *fscale(n)

                stat(i_stat)%c_mean(ispec,2) = &
                     stat(i_stat)%c_mean(ispec,2) + vmr
                stat(i_stat)%c_std(ispec,2) = &
                     stat(i_stat)%c_std(ispec,2) + vmr*vmr

                !========================================
                ! 3D linear interpolation 
                !========================================
                ! this PE handles also the neighbour layer
                if (lstn <= lmr) then    ! lstn = 0 is already excluded (apply vertical EXTRApolation with upper neighbor; see above)
                   rm_interpol = wcx      *wcy      *wcz         *rm(is,js,lst,n)   /m(is,js,lst)   + & 
                                 (1.0-wcx)*wcy      *wcz         *rm(isn,js,lst,n)  /m(isn,js,lst)  + &
                                 wcx      *(1.0-wcy)*wcz         *rm(is,jsn,lst,n)  /m(is,jsn,lst)  + &
                                 (1.0-wcx)*(1.0-wcy)*wcz         *rm(isn,jsn,lst,n) /m(isn,jsn,lst) + &
                                 wcx      *wcy      *(1.0-wcz)   *rm(is,js,lstn,n)  /m(is,js,lstn)  + & 
                                 (1.0-wcx)*wcy      *(1.0-wcz)   *rm(isn,js,lstn,n) /m(isn,js,lstn) + &
                                 wcx      *(1.0-wcy)*(1.0-wcz)   *rm(is,jsn,lstn,n) /m(is,jsn,lstn) + &
                                 (1.0-wcx)*(1.0-wcy)*(1.0-wcz)   *rm(isn,jsn,lstn,n)/m(isn,jsn,lstn) 
                else   !other PE calculates neighbouring contribution
                   !PB                   do NOT apply factor wcz here (no interpolation in z direction in this case) !!!!!

                   rm_interpol = wcx      *wcy          *rm(is,js,lst,n)   /m(is,js,lst)   + & 
                                 (1.0-wcx)*wcy          *rm(isn,js,lst,n)  /m(isn,js,lst)  + &
                                 wcx      *(1.0-wcy)    *rm(is,jsn,lst,n)  /m(is,jsn,lst)  + &
                                 (1.0-wcx)*(1.0-wcy)    *rm(isn,jsn,lst,n) /m(isn,jsn,lst) 

                endif
                vmr = rm_interpol*fscale(n)
                stat(i_stat)%c_mean(ispec,3) = stat(i_stat)%c_mean(ispec,3) + vmr
                stat(i_stat)%c_std(ispec,3)  = stat(i_stat)%c_std(ispec,3) + vmr*vmr


             end do  !ispec


          !======================
          !meteo data at stations
          !======================
          if (output_stations_meteo) then              


             !====================================
             ! 3D linear interpolation Temperature
             !====================================
             T_interpol = &
                  wcx      *wcy      *wcz       *T(is,js,ls)    + & 
                  (1.0-wcx)*wcy      *wcz       *T(isn,js,ls)   + &
                  wcx      *(1.0-wcy)*wcz       *T(is,jsn,ls)   + &
                  (1.0-wcx)*(1.0-wcy)*wcz       *T(isn,jsn,ls)  + &
                  wcx      *wcy      *(1.0-wcz) *T(is,js,lsn)   + & 
                  (1.0-wcx)*wcy      *(1.0-wcz) *T(isn,js,lsn)  + &
                  wcx      *(1.0-wcy)*(1.0-wcz) *T(is,jsn,lsn)  + &
                  (1.0-wcx)*(1.0-wcy)*(1.0-wcz) *T(isn,jsn,lsn) 

             stat(i_stat)%metdat(1) = stat(i_stat)%metdat(1) + T_interpol
             
             !====================================
             ! 3D linear interpolation pressure
             !====================================
             p_interpol =(((0.5-rls) * phlb(is,js,ls)  + (0.5+rls) * phlb(is,js,ls+1))  * wcx         * wcy       + &
                          ((0.5-rls) * phlb(isn,js,ls) + (0.5+rls) * phlb(isn,js,ls+1)) * (1.0-wcx)   * wcy       + &
                          ((0.5-rls) * phlb(is,jsn,ls) + (0.5+rls) * phlb(is,jsn,ls+1)) * wcx         * (1.0-wcy) + &
                          ((0.5-rls) * phlb(isn,jsn,ls)+ (0.5+rls) * phlb(isn,jsn,ls+1))* (1.0-wcx)   * (1.0-wcy) ) *0.01 ![Pa] -> [hPa]

             stat(i_stat)%metdat(2) = stat(i_stat)%metdat(2) + p_interpol

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

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

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

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

             !west border                         !east border  
             windspeed_u_interpol=((0.5-ris) * pu(is_west,js,ls)/m(is_west,js,ls)                                   + &
                                   (0.5+ris) * pu(is,js,ls)/ m(is,js,ls) )      * &
                                   ddx_j  *      wcy *      wcz          + &  
                                   ((0.5-ris) * pu(is_west,jsn,ls)/m(is_west,jsn,ls)                                 + & 
                                   (0.5+ris) * pu(is,jsn,ls)/m(is,jsn,ls) )     * &
                                   ddx_jn * (1.0-wcy)*      wcz          + &  
                                   ((0.5-ris) * pu(is_west,js,lsn)/m(is_west,js,lsn)                                 + &
                                   (0.5+ris) * pu(is,js,lsn)/ m(is,js,lsn))    * &
                                   ddx_j  *      wcy * (1.0-wcz)         + &  
                                   ((0.5-ris) * pu(is_west,jsn,lsn)/m(is_west,js,lsn)                                + &
                                   (0.5+ris) * pu(is,jsn,lsn)/ m(is,jsn,lsn) ) * &
                                   ddx_jn * (1.0-wcy)* (1.0-wcz)

             !====================================
             ! 3D windspeed_v_interpol (y-direction)
             !====================================

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

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

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

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

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


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

             stat(i_stat)%metdat(3) = &
                  stat(i_stat)%metdat(3) + windspeed_u_interpol
             stat(i_stat)%metdat(4) = &
                  stat(i_stat)%metdat(4) + windspeed_v_interpol
             stat(i_stat)%metdat(5) = &
                  stat(i_stat)%metdat(5) + windspeed_interpol

#ifndef without_convection
             !============================================
             ! interpolate the blh to xy position of station
             !============================================
             blh_interpol = &
                  wcx *      wcy*  blh(is,js)  + &
                  (1.0-wcx)*      wcy*  blh(isn,js) + &
                  wcx *(1.0-wcy)*  blh(is,jsn) + &
                  (1.0-wcx)*(1.0-wcy)*  blh(isn,jsn)

             stat(i_stat)%metdat(7) = &
                  stat(i_stat)%metdat(7) + blh_interpol
#endif

          endif  ! meteo out...
          
          stat(i_stat)%ncount = stat(i_stat)%ncount + 1

       end if ! on current tile
    end if ! on current region
    end do ! end i_stat loop

    nullify(m)
    nullify(rm)
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
    nullify(gph)
    nullify(t)
    nullify(pu)
    nullify(pv)
    nullify(phlb)
    nullify(blh)
    ! ok
    status = 0
    
  END SUBROUTINE OUTPUT_STATIONCONC
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RESET_STATIONCONC_ACCUMULATOR
  !
  ! !DESCRIPTION: initialize to 0 accumulation arrays for average station
  !               concentration 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RESET_STATIONCONC_ACCUMULATOR
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer i_stat

    if (station_nspec > 0) then 
       do i_stat=1, n_stat+n_stat2
          stat(i_stat)%c_mean(:,:) = 0.0
          stat(i_stat)%c_std(:,:)  = 0.0
          if (output_stations_meteo) stat(i_stat)%metdat(:)   = 0.0
          stat(i_stat)%ncount  = 0
       end do
    end if
    ! You should also set the _nb-variables to zero. Otherwise, there is a random chance of crashing. Now I am not using them.
    ! do i_stat=1, n_stat
    !    stat(i_stat)%c_mean_nb = 0.0
    !    stat(i_stat)%c_std_nb = 0.0
    !    stat(i_stat)%ncount_nb = 0
    ! end do

  END SUBROUTINE RESET_STATIONCONC_ACCUMULATOR
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EVALUATE_STATIONCONC
  !
  ! !DESCRIPTION: evaluate mean concentration and std
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EVALUATE_STATIONCONC( status )
    !
    ! !USES:
    !
    use dims,        only : idate
    use file_hdf,    only : Init, WriteData, TSds, Done, select
    use file_hdf,    only : CheckInfo
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(inout) :: status
    !
    ! !REVISION HISTORY: 
    !   16 Mar 2012 - P. Le Sager -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer i_stat, i, ii, io_status
    integer n_data, isds, nsds, itp
    real temp
    character(len=12)                    :: xname
    real, dimension(:,:), allocatable    :: cout                 ! to get output
    integer, dimension(:,:), allocatable :: iidate
    real                                 :: windspeed_u_interpol ! wind speed [m/s] x-direction , interpolated to station location
    real                                 :: windspeed_v_interpol ! wind speed [m/s] y-direction , interpolated to station location
    real                                 :: windspeed_interpol   ! wind speed [m/s]             , interpolated to station location
    real                                 :: winddir_interpol     ! wind direction               , interpolated to station location 

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

    !__START_SUBROUTINE______________________________________________________

    if (station_nspec > 0) then
       do i_stat=1, n_stat+n_stat2
          n_data = stat(i_stat)%ncount
          if (n_data > 0) then
             stat(i_stat)%c_mean(:,:) = stat(i_stat)%c_mean(:,:) / n_data
             if (output_stations_meteo) stat(i_stat)%metdat(:)   = stat(i_stat)%metdat(:)   / n_data
             if (n_data > 1) then
                do i=1, station_nspec
                   do ii= 1, n_interpol
                      temp = ((stat(i_stat)%c_std(i,ii) -  n_data * ((stat(i_stat)%c_mean(i,ii))**2) ) / (n_data-1) )
                      if(temp > 0) then 
                         stat(i_stat)%c_std(i,ii) = sqrt(temp)
                      else
                         stat(i_stat)%c_std(i,ii) = 0.0
                      endif
                   enddo
                enddo
             else
                stat(i_stat)%c_std(:,:)=-1.0
             end if
          else
             stat(i_stat)%c_mean(:,:)= -1.0
             stat(i_stat)%c_std(:,:) = -1.0
             if(output_stations_meteo) stat(i_stat)%metdat(:)  = -999.
          end if

          ! Write down wind directions
          if(output_stations_meteo) then 
             windspeed_u_interpol = stat(i_stat)%metdat(3)
             windspeed_v_interpol = stat(i_stat)%metdat(4)
             windspeed_interpol   = stat(i_stat)%metdat(5)

             !====================================
             ! wind direction
             !====================================
             
             ! definition of winddirection:          
             ! east     :  90
             ! south    : 180
             ! west     : 270
             ! north    : 360
             if ((windspeed_u_interpol < 0) .and. (windspeed_v_interpol < 0))then              ! wind from north east (0...90)
                winddir_interpol = atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi
             else if ((windspeed_u_interpol < 0) .and. (windspeed_v_interpol > 0)) then     ! wind from south east (90...180)
                winddir_interpol = 180.0+atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi
             else if ((windspeed_u_interpol < 0) .and. (windspeed_v_interpol == 0)) then    ! wind from east (90)
                winddir_interpol = 90.0
             else if ((windspeed_u_interpol > 0) .and. (windspeed_v_interpol < 0)) then     ! wind from north west (270... 360)
                winddir_interpol = 360.0+atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi
             else if ((windspeed_u_interpol > 0) .and. (windspeed_v_interpol > 0)) then     ! wind from south west (180...270)
                winddir_interpol = 180.0+atan(windspeed_u_interpol / windspeed_v_interpol) * 360.0 / twopi
             else if ((windspeed_u_interpol > 0) .and. (windspeed_v_interpol == 0)) then    ! wind from west (270)
                winddir_interpol = 270.0
             else if ((windspeed_u_interpol == 0) .and. (windspeed_v_interpol < 0)) then    ! wind from north (360)
                winddir_interpol = 360.0
             else if ((windspeed_u_interpol == 0) .and. (windspeed_v_interpol > 0)) then    ! wind from south (180)
                winddir_interpol = 180.0
             else if ((windspeed_u_interpol == 0) .and. (windspeed_v_interpol == 0)) then   ! no wind
                winddir_interpol = -1.0
             else
                write(gol,*)'output_stationconc: error wind direction'; call goErr
                TRACEBACK
                status=1; return
             end if
             stat(i_stat)%metdat(6) = winddir_interpol
          endif
       end do
    end if
    
    !ok
    status = 0
    
  END SUBROUTINE EVALUATE_STATIONCONC
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:   WRITE_STATIONCONC
  !
  ! !DESCRIPTION: write mean concentration
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE WRITE_STATIONCONC( status )
    !
    ! !USES:
    !
    use dims,        only : idate
    use file_hdf,    only : Init, WriteData, TSds, Done
    use file_hdf,    only : WriteAttribute
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(inout)     :: status
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer i, ii, io_status, i_stat
    integer n_data, isds, nsds, itp
    real,    dimension(:,:), allocatable :: cout                 ! to get output
    integer, dimension(:,:), allocatable :: iidate

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

    !__START_SUBROUTINE______________________________________________________

    ! output
    if (station_nspec > 0 ) then 

       if (io_record == 0) then
          call WriteAttribute( io_hdf, 'station_ls', stat(:)%ls, status)
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( io_hdf, 'height_sample', stat(:)%height_sample, status)
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( io_hdf, 'height_surface', stat(:)%height_surface, status)
          IF_NOTOK_RETURN(status=1)
          call WriteAttribute( io_hdf, 'height_sample_above_surf', stat(:)%height_sample_above_surf, status)
          IF_NOTOK_RETURN(status=1)
       endif

       allocate(cout(n_stat+n_stat2,1))

       If (station_nspec > 0) Then
          do i=1, station_nspec
             do ii= 1,n_interpol
                do i_stat = 1, n_stat+n_stat2
                   cout(i_stat,1) = stat(i_stat)%c_mean(i,ii)
                enddo
                call writedata(sds(i,ii), cout, status, start=(/0,io_record/))
                IF_NOTOK_RETURN(status=1)
             enddo
          enddo
       End If
       if (output_stations_meteo) then
          do i=1, n_metdat
             do i_stat = 1, n_stat+n_stat2
                cout(i_stat,1) = stat(i_stat)%metdat(i)
             enddo
             call writedata(sds_meteo(i), cout, status, start=(/0,io_record/))
             IF_NOTOK_RETURN(status=1)
          enddo
       endif

       deallocate(cout)

       allocate(iidate(6,1))
       iidate(:,1) = idate
       call writedata(sds_idate, iidate, status, start=(/0,io_record/) )
       IF_NOTOK_RETURN(status=1)
       deallocate(iidate)

       io_record = io_record + 1

    end if

    ! ok
    status = 0

  END SUBROUTINE WRITE_STATIONCONC
  !EOC

END MODULE USER_OUTPUT_STATION