!#################################################################
!WP! 
!WP! Routine samples specific events in space and time, as read from a list of
!WP! events. These events are created from the NOAA ESRL Flask database, and can
!WP! thus be matched to observations using the eventnumbers. Co-sampling is based on
!WP! a 4-hour window around the observations time and date, and slopes sampling is
!WP! currently used for the interpolation. Contact Wouter Peters
!WP! (Wouter.Peters@noaa.gov) to get the lists of NOAA ESRL event numbers in the
!WP! proper format. Or copy and paste the following lines to a file (remove comment, no header):
!
!ASK_01D0                    23.18    5.42  2728.0 2000   1   1  13  30   0                           9703
!CHR_01D1                     1.70 -157.17     3.0 2000   1   1   3  40   0                          42746
!SPO_01D0                   -89.98  -24.80  2810.0 2000   1   1   5  55   0                         129319
!STM_01D0                    66.00    2.00     7.0 2000   1   1  11  15   0                         134065
!KZD_01D0                    44.45   75.57   412.0 2000   1   3   4  20   0                          69908
!RPB_01D0                    13.17  -59.43    45.0 2000   1   3  19  30   0                         108606
!KEY_01D0                    25.67  -80.20     3.0 2000   1   3  18   0   0                          63211
!MLO_01D0                    19.53 -155.58  3397.0 2000   1   3  20  37   0                          84745
!
!WP! To turn this routine on, add the following RC items to your tm5.rc file:
!
!       output.noaa         : T
!       input.noaa.filename : noaa.list            # get from Wouter
!       inputdir.noaa       : /p73/co2/input/1x1/  # your noaa.list file location here
!       output.noaa.filename: noaa.hdf 
!WP!
!WP! Routine is off as default and will thus not bother other users
!WP! 
!WP! Future improvements:
!WP! - window size can be read from rc-file
!WP! - tracer list can be specific
!WP! - interpolation mode could be an rc-file item
!WP! - or all 3 interpolations could be used
!WP! - more extensive output on model sample (gph, grid, ...)
!WP!
!WP! Created June 9th, 2007 by Wouter Peters
!WP! 
!### 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_noaa

  use GO        , only : gol, goErr, goPr, goLabel, goUpCase
  use dims,       only : im, jm, lm, dx, xref, dy, yref, xbeg, xend, ybeg, yend
  use dims,       only : nregions, region_name, itaui,itaue, xcyc 
  use chem_param, only : fscale, ntracet
  use file_hdf,   only : THdfFile, Init, Done, TSds, SetDim, WriteData, WriteAttribute
  use ParTools
#ifdef MPI
    use mpi_const,  only : my_real,mpi_comm_world,ierr,mpi_integer,mpi_character
    use mpi_comm,   only : barrier
#endif

  implicit none

  public :: init_noaa_events, get_noaa, write_noaa_events, noaa_data

  logical :: noaa_data=.false.

  private

  character(len=100)           :: EventlistFilename = ''
  character(len=100)           :: EventlistOutFilename = ''
  character(len=100)           :: InputDir = './'
  integer                      :: n_event, n_temp
  integer,dimension(6)         :: idate_event
  integer                      :: itau_event
  integer                      :: window = 4*3600 ! 2 hrs around obs hour
  type(THdfFile)               :: io_hdf          ! contains HDF output file

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

  type observation
     integer,dimension(6)       :: idate     ! date of observation
     character(len=24)           :: ident     ! station identifier
     real                       :: lat       ! latitude of station
     real                       :: lon       ! longitude of station
     real                       :: height    ! height of station
     integer                    :: eventnumber       ! unique reference to NOAA ESRL database
     real,dimension(ntracet)    :: mmix=0.0

     integer                    :: itau=0
     integer                    :: n_accumulate=0  ! nr of times concentration was added in get_noaa       

     integer                    :: region 
     integer                    :: is
     integer                    :: js
     integer                    :: lsmin
     integer                    :: lsmax
  end type observation

  type(observation), dimension(:), allocatable  :: concentration

  ! private 

  integer,parameter           :: nf_trace =  1
  real                        :: rmf
  ! number of locations to be calculated for 1 model time
  integer,dimension(nf_trace) :: if_trace = (/ 1 /)


contains

  subroutine write_noaa_events(status)
    !=====================================
    ! HDF output for noaa events 
    ! (averaged model results)
    !=====================================
    use ParTools 
    use Dims,        only : itaui, itaue, ndyn_max, idate, 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 chem_param,  only : ntracet,names
    implicit none
    ! i/o
    integer   :: status

    integer                           :: i_stat, tag,n
    type(TrcFile)                     :: rcF
    character(len=12)                 :: xname
    integer                           :: i, ii, j
    type(TSds), dimension(:), allocatable :: Sds
    type(TSds) :: ds
    real,dimension(:),allocatable   :: data_array
    integer, dimension(:,:), allocatable :: iidate
    character(len=10)                 :: sidatei, sidatee
    ! --- const ------------------------------

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

    call evaluate_noaa()

    if(myid==root) then  !WP! Create output file and write header

       call Init( rcF, rcfile,status)
       call ReadRc( rcF, 'output.noaa.filename', EventlistOutFilename,status)
       call Done( rcF ,status)
       IF_NOTOK_RETURN(status=1)

!WP! Add the date to the filename

       i = index(EventlistOutFilename,'.hdf')
       write(sidatei,'(i4.4,3i2.2)') idatei(1:4)
       write(sidatee,'(i4.4,3i2.2)') idatee(1:4)
       if(i /= 0) then
          EventlistOutFilename = EventlistOutFilename(1:i-1)//'_'//sidatei//'_'//sidatee//'.hdf'
       else
          EventlistOutFilename = trim(EventlistOutFilename)//'_'//sidatei//'_'//sidatee//'.hdf'
       endif

       call Init( io_hdf, trim(OutDir)//'/'// trim(EventlistOutFilename),'create' , status)
       IF_NOTOK_RETURN(status=1)
       call Init(ds, io_hdf, 'event_lat'    ,(/n_event/),'real(4)',status)
       call SetDim( ds, 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status)
       call WriteAttribute(ds, 'Unit', 'deg N',status)
       call WriteData( ds, concentration(1:n_event)%lat,status)
       call Done(ds,status)
       IF_NOTOK_RETURN(status=1)
       call Init(ds, io_hdf, 'event_lon'    ,(/n_event/),'real(4)',status)
       call SetDim( ds, 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status)
       call WriteAttribute(ds, 'Unit', 'deg E',status)
       call WriteData( ds, concentration(1:n_event)%lon,status)
       call Done(ds,status)
       IF_NOTOK_RETURN(status=1)
       call Init(ds, io_hdf, 'event_height' ,(/n_event/),'real(4)',status)
       call SetDim( ds, 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status)
       call WriteData( ds, concentration(1:n_event)%height,status)
       call WriteAttribute(ds, 'Unit', 'meters a.s.l.',status)
       call Done(ds,status)
       IF_NOTOK_RETURN(status=1)
       call Init(ds, io_hdf, 'event_region' ,(/n_event/),'integer(2)',status)
       call SetDim( ds, 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status)
       call WriteData( ds, concentration(1:n_event)%region,status)
       call Done(ds,status)
       IF_NOTOK_RETURN(status=1)
       call Init(ds, io_hdf, 'event_is'     ,(/n_event/),'integer(2)',status)
       call SetDim( ds, 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status)
       call WriteData( ds, concentration(1:n_event)%is,status)
       call Done(ds,status)
       IF_NOTOK_RETURN(status=1)
       call Init(ds, io_hdf, 'event_js'     ,(/n_event/),'integer(2)',status)
       call SetDim( ds, 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status)
       call WriteData( ds, concentration(1:n_event)%js,status)
       call Done(ds,status)
       IF_NOTOK_RETURN(status=1)
       call Init(ds, io_hdf, 'idate', (/6, n_event/), 'integer(2)',status)
       call SetDim( ds, 0, 'n_idate ', 'yymmddhhssmm', (/ (j, j=1,6) /) ,status)
       call SetDim( ds, 1, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status)
       call WriteAttribute(ds, 'Unit', 'year month day hour minutes seconds',status)
       allocate(iidate(6,n_event))
       do j=1,n_event
           iidate(:,j) = concentration(j)%idate
       enddo
       call WriteData( ds, iidate,status)
       IF_NOTOK_RETURN(status=1)
       deallocate(iidate)
       call Done(ds,status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute( io_hdf, 'ntracet', ntracet,status)

       call WriteAttribute( io_hdf, 'n_event', n_event,status)
       call WriteAttribute( io_hdf, 'event_ident'  , concentration(1:n_event)%ident,status)
       call WriteAttribute( io_hdf, 'tracer_names'  , names(1:ntracet),status)
       call WriteAttribute( io_hdf, 'event_nsample', concentration(1:n_event)%n_accumulate,status)
       IF_NOTOK_RETURN(status=1)

       allocate(Sds(ntracet))
       do i=1, ntracet
             xname = trim(names(i))
             call init(Sds(i), io_hdf, xname, (/n_event/), 'real(4)',status)
             call SetDim( Sds(i), 0, 'eventnumber', 'event number', (/ (concentration(j)%eventnumber, j=1,n_event) /) ,status)
             call WriteAttribute(Sds(i), 'Unit', 'Volume mixing ratio',status)
       enddo
       IF_NOTOK_RETURN(status=1)
    endif ! root

    !WP! Now write datasets from each PE to file

    allocate(data_array(n_event))
    do i=1,ntracet
         data_array=concentration(1:n_event)%mmix(i)  ! copy tracer to data_array, proper PE will send to root for I/O

#ifdef MPI
    call MPI_BCAST(data_array(1),n_event,MY_REAL,tracer_id(i),mpi_comm_world,ierr)
#endif
         if(myid==root) then
             call WriteData(Sds(i),data_array,status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds(i),status)
             IF_NOTOK_RETURN(status=1)
         endif
     enddo

     deallocate(data_array)
     if(myid==root)deallocate(Sds)
     deallocate(concentration)

     if(myid==root) then
         call Done(io_hdf,status)
         IF_NOTOK_RETURN(status=1)
     endif

  end subroutine write_noaa_events

  subroutine init_noaa_events(status)
    use Meteo,      only : Set
    use Meteo,      only : temper_dat, humid_dat, gph_dat
    use binas,      only : ae, twopi, grav
    use dims,       only : im, jm, lm, dx, dy, xref, yref, xbeg, ybeg, xend, yend
    use dims,       only : nregions, region_name, itaur, xcyc, idate
    use chem_param, only : ntrace, ntracet, names, fscale
    use toolbox,    only : escape_tm
    use datetime,   only : date2tau
    use go_string,  only : goUpCase   
    use go_string,   only : goNum2Str
    use GO,          only : TrcFile, Init, Done, ReadRc
    use global_data, only : rcfile
    implicit none
    ! i/o
    integer           ::  status
    ! local
    character(len=200)                :: dummy
    integer                           :: fstatus 
    integer                           :: i_stat, j_stat,n
    integer                           :: sunit0=123
    integer                           :: region
    integer,dimension(6)              :: idate_sim
    integer                           :: itau_sim, isim, ispec
    ! x,y resolution (in degrees) for current region 
    real                              :: dxr, dyr 
    character(len=8)                  :: name_spec    !CMK changed from 6-->8
    logical                           :: site_found
    type(TrcFile)                     :: rcF
    integer, dimension(6)     :: temp_idate
    integer                   :: temp_itau,temp_eventnumber
    real                      :: temp_lon,temp_lat,temp_height
    character(len=30)         :: temp_ident
    ! --- const ------------------------------

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

    ! start

    do region=1,nregions

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

    enddo

    if(myid==root) then

       call Init( rcF, rcfile,status)
       call ReadRc( rcF, 'input.noaa.filename', EventlistFilename,status)
       call ReadRc( rcF, 'inputdir.noaa', InputDir,status)
       call Done( rcF,status)
       !=======================================================
       ! open StationFile (with list of stations)
       !   "station" means "observation" in this context (arj)
       !=======================================================

       open(UNIT=sunit0, FORM='FORMATTED', &
             FILE=trim(InputDir)//'/'//trim(EventlistFilename)//trim(goNum2Str(idate(1),'(i4.4)'))//trim(goNum2Str(idate(2),'(i2.2)'))//trim(goNum2Str(1,'(i2.2)')) , &
            ERR=1000)      

       ! count number of stations
       n_temp=0
       stat_loop: DO
          read(sunit0,  '(a)', END=100) dummy
          n_temp = n_temp+1 
       end do stat_loop
100    continue
       print *, '_________________________________________________' &
            //'_____________________________'
       print *, ' preparing to read ', n_temp, ' NOAA events from file: ' //trim(InputDir)//'/'//trim(EventlistFilename)//trim(goNum2Str(idate(1),'(i4.4)'))//trim(goNum2Str(idate(2),'(i2.2)'))//trim(goNum2Str(idate(3),'(i2.2)')) 
       print *, '_________________________________________________' &
            //'_____________________________'

    endif

#ifdef MPI    
    call MPI_BCAST(n_temp           ,1,      MPI_INTEGER, root ,MPI_COMM_WORLD,ierr)  
#endif

    if(n_temp.eq.0) call escape_tm( 'read_noaa_events: no obs in file ' // EventlistFilename) 
    if(allocated(concentration)) deallocate(concentration)  ! avoid double allocation

    allocate(concentration(n_temp))  ! very long array, only partly used


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

    !===============
    ! read stations
    !===============
    n_event=0
    if(myid==root) then
       rewind(sunit0)
       do i_stat=1, n_temp 

          read(sunit0,  '(a24,1x,f8.2,f8.2,f8.1,1x,6i4,1x,i30)') &
               temp_ident, &
               temp_lat, &
               temp_lon, &
               temp_height, &
               temp_idate, &
               temp_eventnumber


          call date2tau(temp_idate,temp_itau)  ! fill itau

          if(abs(temp_itau-itaui).lt.1) temp_itau=temp_itau+3600
          if(abs(temp_itau-itaue).lt.1) temp_itau=temp_itau-3600

          if(temp_itau<itaui.or.temp_itau>itaue) then
            cycle ! exclude events not within this run
          endif
          if(temp_lat<-90.or.temp_lat>90) then
            cycle ! exclude events not within this run
          endif
          if(temp_lon<-180.or.temp_lon>180) then
            cycle ! exclude events not within this run
          endif

          !WP! Proceed to fill array concentration with obs within start and end of this run

          n_event=n_event+1

          ! goUpCase is *probably* unnecessary
          concentration(n_event)%ident=goUpCase(trim(adjustl(temp_ident)))
          concentration(n_event)%lat=temp_lat
          concentration(n_event)%lon=temp_lon
          concentration(n_event)%itau=temp_itau
          concentration(n_event)%height=temp_height
          concentration(n_event)%idate=temp_idate
          concentration(n_event)%eventnumber=temp_eventnumber

          ! Note that the is and js components for the observation locations are
          ! set *before* the site moves are applied.  Is this a bug?
          !  -Andy, 7 Nov 2006

          concentration(n_event)%is = &
               int((concentration(n_event)%lon-float(xbeg(1)))/dxr + 0.99999) 
          concentration(n_event)%js = &
               int((concentration(n_event)%lat-float(ybeg(1)))/dyr + 0.99999)
          concentration(n_event)%lsmin = 999
          concentration(n_event)%lsmax = 0

          print('(a24,1x,f8.2,f8.2,f8.1,1x,6i4,x,i30)'), &
               concentration(n_event)%ident, &
               concentration(n_event)%lat, &
               concentration(n_event)%lon, &
               concentration(n_event)%height, &
               concentration(n_event)%idate, &
               concentration(n_event)%eventnumber

       end do

       ! close stationfile
       close(sunit0)


    endif

! share the number of real events: n_event
!
#ifdef MPI    
    call MPI_BCAST(n_event           ,1,      MPI_INTEGER, root ,MPI_COMM_WORLD,ierr)  
#endif

    if(n_event==0) then ! no data in range
        write(gol,*) 'WARNING: No NOAA ESRL events within this run, skipping output ' ; call goPr
        deallocate(concentration)
        noaa_data=.false.
        status=0
        return
    else
        write(gol,*) 'Found ',n_event,' NOAA ESRL events within this run' ; call goPr
    endif

    ! determine assignment of station to region 
    ! (with highest refinement factor)

#ifdef MPI    
    do i_stat=1, n_event 
       call MPI_BCAST(concentration(i_stat)%ident ,24,        MPI_CHARACTER, &
            root ,MPI_COMM_WORLD,ierr)  
       call MPI_BCAST(concentration(i_stat)%lat   ,1,        MY_REAL      , &
            root ,MPI_COMM_WORLD,ierr)  
       call MPI_BCAST(concentration(i_stat)%lon   ,1,        MY_REAL      , &
            root ,MPI_COMM_WORLD,ierr)  
       call MPI_BCAST(concentration(i_stat)%height,1,        MY_REAL      , &
            root ,MPI_COMM_WORLD,ierr)  
       call MPI_BCAST(concentration(i_stat)%idate,6,        MPI_INTEGER      , &
            root ,MPI_COMM_WORLD,ierr)  
       call MPI_BCAST(concentration(i_stat)%eventnumber,1,        MPI_INTEGER      , &
            root ,MPI_COMM_WORLD,ierr)  
       call MPI_BCAST(concentration(i_stat)%itau,1,        MPI_INTEGER      , &
            root ,MPI_COMM_WORLD,ierr)  
       call MPI_BCAST(concentration(i_stat)%is    ,1,        MPI_INTEGER,   &
            root ,MPI_COMM_WORLD,ierr)  
       call MPI_BCAST(concentration(i_stat)%js    ,1,        MPI_INTEGER,   &
            root ,MPI_COMM_WORLD,ierr)  
       call MPI_BCAST(concentration(i_stat)%lsmin ,1,        MPI_INTEGER,   &
            root ,MPI_COMM_WORLD,ierr)  
       call MPI_BCAST(concentration(i_stat)%lsmax ,1,        MPI_INTEGER,   &
            root ,MPI_COMM_WORLD,ierr)  
    end do
#endif
    do i_stat=1, n_event

       ! assume global region as default for average mixing ratios
       concentration(i_stat)%region = 1
       concentration(i_stat)%is=1
       concentration(i_stat)%js=1
       concentration(i_stat)%n_accumulate = 0 ! no model values addded yet
       concentration(i_stat)%mmix(:) = 0.0 ! no model values addded yet

       do region=1, nregions

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

          if ( (concentration(i_stat)%lon .gt. xbeg(region) .and. &
               concentration(i_stat)%lon .lt. xend(region)) .and.      &
               (concentration(i_stat)%lat .gt. ybeg(region) .and. &
               concentration(i_stat)%lat.lt.yend(region) ) ) then
             !=====================
             ! station is in region
             !=====================

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

    call Par_Barrier

    return

    ! error handling
1000 call escape_tm( 'read_noaa_events: could not open ' // trim(InputDir)//'/'//EventlistFilename) 

  end subroutine init_noaa_events

  subroutine get_noaa(region,itau_f,force)
    !
    ! This subroutine samples the model at given locations in array
    ! concentrations and puts the values in the tag: forecast
    ! This is done only in a certain time window, or when the keyword /force is
    ! used. Thus, a model field can be sampled at sites instantaneously (force),
    ! or at appropriate time ( omit force)
    !
    ! Note: the values accumulate so care should be taken to clear the forecast
    ! and n_accumulate tags between functions
    !
    !
    use global_data, only: mass_dat, region_dat 
    use Meteo,       only: gph_dat, m_dat
    use ParTools,    only: myid,ntracet_ar

    implicit none

    ! input/output
    integer,intent(in)              :: region
    ! idate_f: date for which output required...
    integer,intent(in)              :: itau_f
    logical,optional                :: force

    ! local
    real,dimension(:,:,:), pointer      :: m,gph
    real,dimension(:,:,:,:), pointer    :: rm, rxm, rym, rzm

    real,dimension(0:lm(region)) :: height
    integer   :: i,is,js,l,n,isn,jsn,ls,j,offsetj, lst, lstn, lmr, lsn
    real      :: flon,flat,falt,ris,rjs,dxr,dyr,wcx,wcy,rls, wcz


    ! start

    m => m_dat(region)%data    !pointers to global arrays...
    rm => mass_dat(region)%rm_t
    rxm => mass_dat(region)%rxm_t
    rym => mass_dat(region)%rym_t
    rzm => mass_dat(region)%rzm_t
    gph => gph_dat(region)%data

    lmr=lm(region)

    fcast_loop: do n=1,n_event 

       !  0. Is idate equal to idate_site
       !  1. Is the site in the area?--->no, then put -1 in c
       !  2. Determine gridbox
       !  3. Use slopes to determine concentration at the site.

       if ( (concentration(n)%itau.ge.itau_f-window/2.and.&
            concentration(n)%itau.lt.itau_f+window/2).or.force ) then


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

          flon=concentration(n)%lon
          flat=concentration(n)%lat
          falt=concentration(n)%height

          if (concentration(n)%region == region) then
             ris = (flon-float(xbeg(region)))/dxr + 0.99999
             rjs = (flat-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
             !fraction from the center of the is-box  (-0.5---+0.5)
             ris = ris-is-0.5
             !idem js
             rjs = rjs-js-0.5

             !the neighbour for pressure 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 < 1) jsn=1
             if ( jsn > jm(region) ) jsn=jm(region)  ! isn-->jsn (wouter Peters)
             if ( xcyc(region) == 0 ) then
                ! non-cyclic boundaries
                if ( isn < 1) isn=1
                if ( isn > im(region) ) isn=im(region)
             else
                ! cyclic x-boundaries
                if ( isn < 1 ) isn=im(region)
                if ( isn > im(region) ) isn=1
             end if
             
             ! interpolate the altitude to site position...
             ls = 1   !layer
             do l=0,lm(region)
                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)
             enddo
             do l=2,lm(region)   ! selects layer , note that we start from second layer from surface
                if(height(l).gt.falt) exit
             enddo

             select case(l)
             case(0)
                if(myid==root) print*,'get_noaa: Warning..., forecast altitude ',&
                     'is below the surface height',height(0),' > ',falt,concentration(n)%ident
                ls = 1
                rls = -0.5   !surface...
             case default
                ls = l       !the site layer
                ! the off-set from the center of the layer (-0.5--->+0.5)
                ! (interpolation is in (m))
                rls = (falt-height(l-1))/(height(l)-height(l-1)) - 0.5
             end select
             !=================================
             !the neighbour for z interpolation
             !=================================
             if ( rls .gt. 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


             !from is,js,ls, ris,rjs,rls, determine the mixing ratio of each tracer  ...
             do j=1,ntracetloc
                ! rm-value is obtained from rm + slopes. 
                ! slope = rxm = (rm*dX/dx *deltaX/2)

                ! For EnKF forward, not sure that sum over a zero
                ! length subarray works properly.  Let's be specific
                ! anyway.
                offsetj = sum(ntracet_ar(0:myid-1))

                rmf = (      rm(is,js,ls,j) + &
                     2.0*(ris*rxm(is,js,ls,j) + &
                          rjs*rym(is,js,ls,j) + &
                          rls*rzm(is,js,ls,j) ) )/m(is,js,ls)*fscale(offsetj+j) 

                concentration(n)%mmix(offsetj+j)=concentration(n)%mmix(offsetj+j)+rmf !*fscale(offsetj+j)

             enddo

             concentration(n)%n_accumulate=concentration(n)%n_accumulate+1
!!             if(myid==root) then
!                print*,'DEBUG Check ',myid, concentration(n)%ident,concentration(n)%n_accumulate,concentration(n)%model(:)/concentration(n)%n_accumulate
!             endif

          endif
       endif ! force/time

    enddo fcast_loop

    nullify(m)
    nullify(rm)
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
    nullify(gph)


  end subroutine get_noaa

  subroutine evaluate_noaa

    ! convert concentrations to ppm and average over n_accumulate

    integer :: n

    do n=1,ntracet
       where(concentration%n_accumulate.ge.1) &
            concentration%mmix(n)=1.e6*concentration%mmix(n)/concentration%n_accumulate
    enddo


  end subroutine evaluate_noaa

end module user_output_noaa
