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

module user_output_station

  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 : myid, root, tracer_active, tracer_loc, ntracet_ar

  implicit none

  private

  public :: read_stationlist, stat, stations
  public :: init_station_output
  public :: evaluate_stationconc
  public :: reset_stationconc_accumulator
  public :: write_Stationfiles, StationlistFFilename , output_stationconc
  public :: free_stationfields, write_stationconc
  public :: n_stat, io_nrecords

#ifdef with_optics
  public :: sds_station_aod,sds_station_alpha
#endif
#ifdef with_pm
  public :: sds_station_pm
#endif
! and now a very advanced construction thanks to internet
#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

  ! 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 

  character(len=100)           :: StationlistFFilename   ! list of 'regular' stations
  character(len=100)           :: StationlistFFilename2  ! list of 'validation' stations

!  character(len=100)           :: StationOutFFilename 
  character(len=100)           :: InputDir = './'
  logical                      :: output_stations_meteo = .true.
  logical                      :: write_Stationfiles = .true.
  character(len=300)           :: StationSpeclist 
  integer                      :: StationlistFileFormat  ! 1 : old format
                                                         ! 2:  new format 

  integer,save                 :: n_stat            ! number of 'regular' stations (used for 4DVAR) 
  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
  integer,save                 :: io_nrecords 
  ! 'open' datasets of the io_hdf output file:
  type(TSds), dimension(:,:), allocatable :: Sds
  type(TSds), dimension(n_metdat)         :: Sds_meteo
  type(TSds)                              :: Sds_idate

  type 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

     ! 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

  type(stations), dimension(:), allocatable  :: stat

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

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

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

contains

!===========================================================================================================
!===========================================================================================================

  subroutine free_stationfields(status)

    implicit none

    !__IO____________________________________________________________________

    integer, intent(inout)   :: status

    !__LOCAL_VARIABLES_______________________________________________________

    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 )
    endif
    !ok
    status = 0
  end subroutine free_stationfields

!===========================================================================================================
!===========================================================================================================

  subroutine read_stationlist( status )

    use GO,          only : TrcFile, Init, Done, ReadRc
    use global_data, only : rcfile
    use chem_param,  only : ntracet, names

    implicit none

    !__IO____________________________________________________________________

    integer, intent(inout)            :: status

    !__LOCAL_VARIABLES_______________________________________________________

    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
    ! x,y resolution (in degrees) for current region 
    real                              :: dxr, dyr 
    character(len=8)                  :: name_spec    !CMK changed from 6-->8
    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)

    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
    print*,n_stat

    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 Which Pe handles the output of the required species:

    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.'all'.and.tracer_active(n) ) then
             station_nspec = station_nspec + 1
             station_speclist(station_nspec) = tracer_loc(n)
          else if ((trim(name_spec).eq.trim(names(n)) .and. tracer_active(n)) ) then
             station_nspec = station_nspec + 1
             station_speclist(station_nspec) = tracer_loc(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
    
    if (station_nspec > 0) then
       allocate(stat(n_stat+n_stat2))

       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

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


       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
             stat(i_stat)%is = &
                  int((stat(i_stat)%lon-float(xbeg(1)))/dxr + 0.99999) 
             stat(i_stat)%js = &
                  int((stat(i_stat)%lat-float(ybeg(1)))/dyr + 0.99999)
          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(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
             stat(i_stat)%is = &
                  int((stat(i_stat)%lon-float(xbeg(1)))/dxr + 0.99999) 
             stat(i_stat)%js = &
                  int((stat(i_stat)%lat-float(ybeg(1)))/dyr + 0.99999)


          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
                   stat(i_stat)%is = &
                        int((stat(i_stat)%lon-float(xbeg(region)))/dxr + 0.99999) 
                   stat(i_stat)%js = &
                        int((stat(i_stat)%lat-float(ybeg(region)))/dyr + 0.99999)
                end if
             end if
          end do
       end do

       if (myid == 0) 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 ! station nspec > 0
    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

!===========================================================================================================
!===========================================================================================================

  subroutine init_station_output(status)
    !=====================================
    ! HDF output for stations 
    ! (instantaneous model results)
    !=====================================
    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

    implicit none

    ! --- in/out -----------------------------
    
    integer, intent(out)      ::  status

    !__LOCAL_VARIABLES_______________________________________________________

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

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

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

    !__START_SUBROUTINE______________________________________________________

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

    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
       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(offsetn+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(offsetn + 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 ratio (ppb)', status)
             IF_NOTOK_RETURN(status=1)
          enddo
       enddo
       if (myid /= 0) 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

!===========================================================================================================
!===========================================================================================================

  subroutine close_stationfiles( status )

    use file_hdf,  only: Done
    implicit none

    !__IO____________________________________________________________________

    integer, intent(inout)            :: status

    !__LOCAL_VARIABLES_______________________________________________________

    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 (myid /= 0) 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

!===========================================================================================================
!===========================================================================================================

  subroutine output_stationconc(region, status)

    use global_data, only : mass_dat
#ifndef without_convection
    use global_data, only : conv_dat
#endif
    use Meteo,       only : m_dat, gph_dat, temper_dat, pu_dat, pv_dat, phlb_dat
    use datetime,    only : tau2date

    implicit none


    !__IO____________________________________________________________________

    integer,intent(in)   :: region
    integer,intent(inout):: status

    !__LOCAL_VARIABLES_______________________________________________________

    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
    integer                         :: offsetn              ! for paralel processing

    !=============
    ! 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_t
    rxm => mass_dat(region)%rxm_t
    rym => mass_dat(region)%rym_t
    rzm => mass_dat(region)%rzm_t

    ! 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
    !====================

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

    do i_stat=1, n_stat+n_stat2

       if (stat(i_stat)%region == region)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(offsetn+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(offsetn+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(offsetn+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(offsetn+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(offsetn+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(offsetn+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(offsetn+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(offsetn+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(offsetn+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
    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

!===========================================================================================================
!===========================================================================================================

  subroutine reset_stationconc_accumulator

    !==================================================================
    ! initialize accumulation arrays for average station concentration
    !==================================================================

    implicit none

    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
    endif

  end subroutine reset_stationconc_accumulator

!===========================================================================================================
!===========================================================================================================

  subroutine evaluate_stationconc(status)

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

    use dims,        only : idate
    use file_hdf,    only : Init, WriteData, TSds, Done, select
    use file_hdf,    only : CheckInfo


    implicit none

    !__IO____________________________________________________________________

    integer, intent(inout) :: status

    !__LOCAL_VARIABLES_______________________________________________________


    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(myid /= root) return
    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
          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
    endif

    !ok
    status = 0
  end subroutine evaluate_stationconc

!===========================================================================================================
!===========================================================================================================

  subroutine write_stationconc(status)

    !==================================================================
    ! write mean concentration
    !==================================================================

    use dims,        only : idate
    use file_hdf,    only : Init, WriteData, TSds, Done
    use file_hdf,    only : WriteAttribute

    implicit none


    !__IO____________________________________________________________________

    integer, intent(inout)     :: status

    !__LOCAL_VARIABLES_______________________________________________________


    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))

       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
       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
    endif
    ! ok
    status = 0
  end subroutine write_stationconc

!===========================================================================================================
!===========================================================================================================

end module user_output_station
