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

! ********************* DUMMY for TM5-MP **************************


!#################################################################
!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       
!      real                       :: weight ! accumulated time weighting, replaces n_accumulate in evaluation

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

     STATUS=0


     !     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

     STATUS=0

     
!     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)%weight = 0.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
!     use dims,		only : ndyn,ndyn_max
!     implicit none

!     ! input/output
     integer,intent(in)              :: region
!     ! idate_f: date for which output required...
     integer(kind=8),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
!     real      :: this_weight

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

!     this_weight=real(ndyn)/real(ndyn_max)

!     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
!              concentration(n)%weight=concentration(n)%weight + this_weight

! !!             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%weight
!     enddo


!   end subroutine evaluate_noaa

end module user_output_noaa