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

#include "tm5.inc"
!
!#################################################################
module user_output_forecast

  use GO        , only : gol, goErr, goPr, goLabel, goUpCase
  use GO,         only : TrcFile, Init, Done, ReadRc
  use global_data,only : rcfile
  use chem_param, only : ntrace
  use ParTools,   only : myid, root, ntracetloc, par_barrier

  implicit none


  character(len=200)           :: dasinputdir = ''
  character(len=200)           :: ObslistFFilename = ''
  character(len=200)           :: SamplelistFFilename = ''
  integer                      :: n_fcast

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

  ! ____ for rc interface ____

  type(TrcFile)                      :: rcF

  ! ____ local data structures _____

  type point_observation

     !WP! The following information is associated with the observation itself, i.e.,
     !WP! data values and auxiliary information about the sample.

     integer,dimension(6)            :: idate = (/-9,-9,-9,-9,-9,-9/)           ! date of observation in integer*6 format (YY,MM,DD,hh,mm,ss)
     real                            :: decdate = -999.9     ! date of observation in decimal format
     integer                         :: itau = -999          ! date of observation in "seconds since YYYY/MM/DD format"
     real                            :: t_window = -999.9    ! number of seconds to sample around central time itau
     character(len=24)               :: ident = 'none'       ! station identifier
     real                            :: lat = -999.9         ! latitude of station
     real                            :: lon = -999.9         ! longitude of station
     real                            :: height = -999.9      ! height of station
     character(len=24)               :: eventnumber = 'none' ! eventnumber comes from laboratory database
     integer                         :: identifier = -999    ! identifier is a unique key assigned during this run
     character(len=10)               :: species = 'none'     ! tracer to sample for this obs
     real                            :: obs = -999.9         ! observed value (from file)
     integer                         :: strategy = 999       ! strategy of sampling
     character(len=3)                :: ccg_flags='...'      ! flag string from the NOAA ccg database

     !WP! The following information is not associated with the observation, but gives
     !WP! tracer values and auxiliary info from the TM5 model sample

     real,dimension(ntrace)          :: model = -999.9       ! model predicted background value
     real                            :: gph = -999.9         ! the geopotential height of the lowest model level, averaged like obs
     integer                         :: region = -999        ! TM5 zoom region sampled
     integer                         :: is = -999            ! the grid box index i (zonal index)
     integer                         :: js = -999            ! the grid box index j (meridional index)
     integer                         :: ls = -999            ! the grid box index l (vertical index)   
     real                            :: ris = -999.9         ! the relative x position inside the grid box (-0.5 to +0.5)
     real                            :: rjs = -999.9         ! the relative y position inside the grid box (-0.5 to +0.5)
     integer                         :: n_accumulate = 0     ! nr of times sample was added in get_forecast       

  end type point_observation


  type(point_observation), dimension(:), allocatable  :: sample

  private

  public :: user_output_forecast_init, get_sample, user_output_forecast_done, write_sample_info
  public :: read_sample_info

contains

  subroutine user_output_forecast_init(status)
  use Meteo,               only : Set
  use Meteo,               only : gph_dat
  use dims,                only: nregions

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

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

  ! local

  integer :: region


  ! Set the flags to use Meteo data that we will need in this routine

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

  !WP! Read site information rc-file. I propose to move this to the higher level
  !WP! shell considering the quite complex manipulations that will be required
  !WP! based on the site information. The results of this (parameters such as
  !WP! duplicate, site_move, may_localize, etc) can then be added to the
  !WP! observation.nc file by the scripts, and saved as an intermediate step between
  !WP! the raw observation.nc file coming from Ken, and the observation++.nc file
  !WP! read by TM5.

  !call read_sites(status)
  !IF_NOTOK_RETURN(status=1)

  ! Read observation value NetCDF file

  call read_sample_info(status)
  IF_NOTOK_RETURN(status=1)

  call add_tm5_metadata(status)
  IF_NOTOK_RETURN(status=1)

  ! all done, all ok

  status = 0

  end subroutine user_output_forecast_init

  subroutine read_sample_info(status)
    use dims,       only : idate
    use toolbox,    only : escape_tm
    use go_string,  only : goUpCase, goNum2str   
    use file_hdf

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

    !_____ for HDF interface ____

    type(THdfFile)               :: io_hdf      
    type(TSds)                   :: sds

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

    !___  local ___

    integer, dimension(1)                       :: data_dims
    integer                                     :: n
    integer, dimension(:,:),allocatable         :: dummy_idate
    character(len=10), dimension(:),allocatable :: dummy_species
    character(len=24), dimension(:),allocatable :: dummy_sites


    ! start

    ! Get name of obs list file

    call Init( rcF, trim(rcfile), status)
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'carbontracker.obs.inputfile', ObslistFFilename, status,default='obs_forecast.list')
    write(gol,'(a)') 'Observation input data will be read from: '//trim(ObslistFFilename) ; call goPr
    IF_NOTOK_RETURN(status=1)
    ObslistFFilename = trim(ObslistFFilename)// &
              !'.'//         &
              !trim(goNum2str(idate(1),'(i4.4)'))//  &
              !trim(goNum2str(idate(2),'(i2.2)'))//  &
              !trim(goNum2str(idate(3),'(i2.2)'))//  &
              '.nc'
    write(gol,'(a)') 'Observation input data will be read from: '//trim(ObslistFFilename) ; call goPr

    call Done( rcF, status )

    !=======================================================
    ! open StationFile (with list of stations)
    !   "station" means "observation" in this context (arj)
    !=======================================================

    ! start access to NetCDF observation file
    call Init( io_hdf, trim(ObslistFFilename), 'read', status )
    IF_ERROR_RETURN(status=1)
    write(gol,'(a)') '...File opened' ; call goPr

    ! access data set with identifiers and get dimension of observation array
    call Init(sds,io_hdf,'id',status)
    IF_ERROR_RETURN(status=1)
    call GetInfo(sds,status,data_dims=data_dims)
    IF_ERROR_RETURN(status=1)
    write(gol,'(a,i3)') '......Dimension of observation array read: ', data_dims(1) ; call goPr

    ! assign dimension to n_fcast and create array of type "sample" for
    ! further reading of NetCDF data

    n_fcast=data_dims(1)
    if(n_fcast.eq.0) call escape_tm( 'read_sample_info: no obs in file ' // trim(ObslistFFilename)) 
    allocate(sample(n_fcast))  

    write(gol,'(a)') '......Data structure for observations allocated' ; call goPr

    ! read the observation identifier array
    call ReadData( sds, sample%identifier, status )
    IF_ERROR_RETURN(status=1)
    write(gol,'(a)') '......Observation identifiers read' ; call goPr
    call Done(sds,status)
    IF_ERROR_RETURN(status=1)

    ! read date array
    call Init(sds,io_hdf,'decimal_date',status)
    IF_ERROR_RETURN(status=1)
    call ReadData( sds, sample%decdate, status )
    IF_ERROR_RETURN(status=1)
    write(gol,'(a)') '......Observation decimal date' ; call goPr
    call Done(sds,status)
    IF_ERROR_RETURN(status=1)

    ! read latitude array
    call Init(sds,io_hdf,'lat',status)
    IF_ERROR_RETURN(status=1)
    call ReadData( sds, sample%lat, status )
    IF_ERROR_RETURN(status=1)
    write(gol,'(a)') '......Observation latitudes read' ; call goPr
    call Done(sds,status)
    IF_ERROR_RETURN(status=1)

    ! read longitude array
    call Init(sds,io_hdf,'lon',status)
    IF_ERROR_RETURN(status=1)
    call ReadData( sds, sample%lon, status )
    IF_ERROR_RETURN(status=1)
    write(gol,'(a)') '......Observation longitudes read' ; call goPr
    call Done(sds,status)
    IF_ERROR_RETURN(status=1)

    ! read height array
    call Init(sds,io_hdf,'alt',status)
    IF_ERROR_RETURN(status=1)
    call ReadData( sds, sample%height, status )
    IF_ERROR_RETURN(status=1)
    write(gol,'(a)') '......Observation heights read' ; call goPr
    call Done(sds,status)
    IF_ERROR_RETURN(status=1)

    ! read sampling stragegy array
    call Init(sds,io_hdf,'sampling_strategy',status)
    IF_ERROR_RETURN(status=1)
    call ReadData( sds, sample%strategy, status )
    IF_ERROR_RETURN(status=1)
    write(gol,'(a)') '......Observation strategy read' ; call goPr
    call Done(sds,status)
    IF_ERROR_RETURN(status=1)

    ! read idate array
    call Init(sds,io_hdf,'date_components',status)
    IF_ERROR_RETURN(status=1)
    allocate(dummy_idate(6,n_fcast))
    call ReadData( sds, dummy_idate, status)
    IF_ERROR_RETURN(status=1)
    do n=1,n_fcast
        sample(n)%idate(:) = dummy_idate(:,n)
    end do
    deallocate(dummy_idate)
    write(gol,'(a)') '......Observation date_components read' ; call goPr
    call Done(sds,status)
    IF_ERROR_RETURN(status=1)

    ! read species array
    !call Init(sds,io_hdf,'species',status)
    !IF_ERROR_RETURN(status=1)
    !allocate(dummy_species(n_fcast))
    !call ReadData( sds, dummy_species, status)
    !IF_ERROR_RETURN(status=1)
    do n=1,n_fcast
        sample(n)%species = 'dummy'
    end do
    !deallocate(dummy_species)
    write(gol,'(a)') '......Observation species array read' ; call goPr
    !call Done(sds,status)
    !IF_ERROR_RETURN(status=1)

    ! read site name array
    !call Init(sds,io_hdf,'site',status)
    !IF_ERROR_RETURN(status=1)
    !allocate(dummy_sites(n_fcast))
    !call ReadData( sds, dummy_sites, status)
    !IF_ERROR_RETURN(status=1)
    do n=1,n_fcast
        sample(n)%ident = 'dummy_site'
    end do
    !deallocate(dummy_sites)
    write(gol,'(a)') '......Observation sites array read' ; call goPr
    !call Done(sds,status)
    !IF_ERROR_RETURN(status=1)

    call Done(io_hdf,status)
    IF_ERROR_RETURN(status=1)

    write(gol,'(a)') '...File closed' ; call goPr

    !WP! Print summary of all samples read

    !do n=1,n_fcast
    !    
    !      write(gol,'(i10,f8.2,f8.2,f8.1,1x,6i4,1x,i3,1x,a24,1x,a10)'), &
    !      !write(gol,*), &
    !           sample(n)%identifier, &
    !           sample(n)%lat, &
    !           sample(n)%lon, &
    !           sample(n)%height, &
    !           sample(n)%idate, &
    !           sample(n)%strategy, &
    !           trim(sample(n)%ident), &
    !           trim(sample(n)%species)
    !       call goPr

    !end do

  end subroutine read_sample_info

  subroutine add_tm5_metadata(status)
    use dims,       only : nregions
    use dims,       only : im, jm, lm, dx, dy, xref, yref, xbeg, ybeg, xend, yend
    use datetime,   only : date2tau

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

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

    !___  local ___

    real                                        :: dxr, dyr 
    integer                                     :: i_stat, region
    integer                                     :: is,js
    real                                        :: flon, flat, ris, rjs 

    !WP! Determine to which zoom region this sample belongs, and initialize some
    !WP! values needed later

    do i_stat=1, n_fcast

       ! assume global region as default for average mixing ratios
       sample(i_stat)%region = 1
       sample(i_stat)%n_accumulate = 0 ! no model values addded yet
       sample(i_stat)%t_window = 4*3600 ! to be read from file later on?

       call date2tau(sample(i_stat)%idate,sample(i_stat)%itau)  ! fill itau

       do region=1, nregions

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

             sample(i_stat)%region = region

            !WP! Figure out the (x,y) gridbox in which the sample is located

             flon=sample(i_stat)%lon
             flat=sample(i_stat)%lat

             dyr = dy/yref(region)
             dxr = dx/xref(region)
             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
             ris = ris-is-0.5 ! fraction from the center of the is-box  (-0.5---+0.5)
             rjs = rjs-js-0.5 ! idem js

             sample(i_stat)%is = is
             sample(i_stat)%js = js
             sample(i_stat)%ris = ris
             sample(i_stat)%rjs = rjs

            !WP! Sample height is calculated during the sampling and not here as gph changes
            !WP! in time so we need the hour-specific information

          end if
       end do  ! regions
    end do   ! samples

    status = 0

  end subroutine add_tm5_metadata

  subroutine get_sample(region,status)
    !
    ! This subroutine samples the model at given locations 
    !
    ! Note: the values accumulate so care should be taken to clear the forecast
    ! and n_accumulate tags between function calls. 
    !
    ! Extension by Andy, spring 2007: see evaluate_forecast routine.  It does
    ! the unit conversion and n_accumulate averaging.  Is supposed to be robust
    ! in face of accidentally calling the routine multiple times.

    use global_data, only: mass_dat, region_dat 
    use meteodata,   only: m_dat
    use Meteo,       only: gph_dat
    use ParTools,    only: myid,ntracet_ar
    use dims,        only: im,jm,lm, itaur, xcyc
    use chem_param,  only: ntracet, fscale

    implicit none
    ! --- in/out -----------------------------
    
    integer, intent(inout)          ::  status
    integer,intent(in)              :: region

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

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

    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, rmf, window


    ! 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
    gph => gph_dat(region)%data

    lmr=lm(region)

    fcast_loop: do n=1,n_fcast 
       !  0. Is the site in the area?
       !  1. Is te sampling window straddling the current time?
       !  2. Determine gridbox
       !  3. Use slopes to determine sample at the site.

       if (sample(n)%region /= region) cycle   !WP! Only sites in this zoom region

       window = sample(n)%t_window

       if ( sample(n)%itau.lt.itaur(region)-window/2.or.&    
            sample(n)%itau.ge.itaur(region)+window/2 ) cycle !WP! Only sites in this time interval

       if (sample(n)%n_accumulate == 0) then 
            write(gol,'(a,a)') 'New sample window was entered, recording mixing ratios for ',trim(sample(n)%ident)
            call goPr
	    sample(n)%model = 0.0 !WP! Always start from zero
       endif


       falt = sample(n)%height
       is   = sample(n)%is
       js   = sample(n)%js
       ris  = sample(n)%ris
       rjs  = sample(n)%rjs

       !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 vertical model structure to sample 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 ls=0,lm(region)   ! selects layer where the model height exceends the sample height
                           ! note that we start from second layer from surface
                           ! to avoid sampling biases near the surface
          if(height(ls).gt.falt) exit
       enddo

       if (ls==0) then
     !       write(gol,'("get_sample: Warning..., forecast altitude ",&
     !            "is below the surface height",f8.2," > ", f8.2, a)'),height(0),falt,trim(sample(n)%ident)

    !        write(gol,'("get_sample: Warning..., forecast altitude ","is below the surface height",f," > ", f, a)'),height(0),falt,trim(sample(n)%ident)


            call goPr
            ls = 1
       endif

       sample(n)%ls = ls

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

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

          sample(n)%model(j)=sample(n)%model(j)+rmf*fscale(j) 

     enddo

     sample(n)%n_accumulate=sample(n)%n_accumulate+1

     enddo fcast_loop

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

    status = 0

  end subroutine get_sample

  subroutine evaluate_sample(status)
    use chem_param, only : ntracet
    implicit none
    ! --- in/out -----------------------------
    
    integer, intent(inout)          ::  status

    ! convert concentrations to ppm and average over n_accumulate

    integer :: n

    do n=1,ntracet
       where(sample%n_accumulate.ge.1) 
          sample%model(n)=sample%model(n)/sample%n_accumulate
       end where
    enddo

    sample%n_accumulate = 0 

    status = 0

  end subroutine evaluate_sample


  subroutine write_sample_info(status)
    use chem_param, only : my_member_no
    use dims,       only : idate
    use dims,       only : nregions
    use toolbox,    only : escape_tm
    use go_string,  only : goUpCase, goNum2str   
    use global_data,only : outdir
    use file_hdf
    use chem_param, only : my_member_no

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

    !_____ for HDF interface ____

    type(TSdsDim)                :: dim
    type(TTimeSeriesHDF)         :: F
    type(TSds)                   :: sds

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

    !___  local ___

    integer :: n, data_rank

    !WP! average samples before writing

    call evaluate_sample(status) 

    !WP! testing output to NetCDF through Records interface

!    SamplelistFFilename = trim(outdir)//'/samples.'//         &
!                          trim(goNum2str(idate(1),'(i4.4)'))//  &
!                          trim(goNum2str(idate(2),'(i2.2)'))//  &
!                          trim(goNum2str(idate(3),'(i2.2)'))//  &
!                          '.hdf'

    SamplelistFFilename = trim(outdir)//'/samples.'//trim(my_member_no)//'.nc'

    write(gol,'(a)') 'Full sample output list filename is '//trim(SamplelistFFilename) ; call goPr

    call Init(F,trim(SamplelistFFilename), 400, status)
    IF_NOTOK_RETURN(status=1)

    !WP! Write one record of each of the variables to file

    do n=1,n_fcast
          call AddRecord(F,'id','sample ID','none','real(4)',sample(n)%identifier, status)
          IF_NOTOK_RETURN(status=1)
          call AddRecord(F,'decimal_date','decimal_date','seconds since yesterday','real(8)',sample(n)%decdate, status)
          IF_NOTOK_RETURN(status=1)
          call AddRecord(F,'lon','longitude','degrees_east','real(4)',sample(n)%lon, status)
          IF_NOTOK_RETURN(status=1)
          call AddRecord(F,'lat','latitude','degrees_north','real(4)',sample(n)%lat, status)
          IF_NOTOK_RETURN(status=1)
          call AddRecord(F,'date_components','integer*6 idate','yyyy/mm/dd/HH/MM/SS','integer',(/6/),sample(n)%idate(:), status)
          IF_NOTOK_RETURN(status=1)
          call AddRecord(F,'sampled_mole_frac','sampled_model_mole_fractions','mol/mol','real(8)',(/ntrace/),sample(n)%model(:), status)
          IF_NOTOK_RETURN(status=1)
          call AddRecord(F,'n_accumulate','n_samples_in_average','none','integer',sample(n)%n_accumulate, status)
          IF_NOTOK_RETURN(status=1)
    enddo

    !WP! Clean up the output file by changing names of dimensions and adding
    !WP! useful attributes

    !WP! Change name of the unlimited dimension to 'time' for each dataset

    do n=0,F%nfield-1

          call Init(sds,status)
          IF_NOTOK_RETURN(status=1)
          call Select(sds, F%hdf , n , status)
          IF_NOTOK_RETURN(status=1)
          call GetInfo( sds, status, data_rank=data_rank)
          IF_NOTOK_RETURN(status=1)
          call Init(dim,status)
          IF_NOTOK_RETURN(status=1)
          call Select(dim, sds , data_rank-1, status)
          IF_NOTOK_RETURN(status=1)
          call SetName( dim, 'time',status)
          IF_NOTOK_RETURN(status=1)
          call Done(dim,status)
          IF_NOTOK_RETURN(status=1)
          call Done(sds,status)
          IF_NOTOK_RETURN(status=1)

    enddo

    !WP! Close file

    call Done(F,status)
    IF_NOTOK_RETURN(status=1)

    status = 0

  end subroutine write_sample_info

  subroutine user_output_forecast_done(status)
  use dims,                only: nregions

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

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

  ! local

  integer :: region


  ! Set the flags to use Meteo data that we will need in this routine

  call write_sample_info(status)
  IF_NOTOK_RETURN(status=1)

  deallocate(sample)

  ! all done, all ok

  status = 0

  end subroutine user_output_forecast_done

end module user_output_forecast
