!### 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_flight

  use dims,       only : lm, dx, xref, dy, yref, xbeg, xend, ybeg, yend
  use dims,       only : nregions, region_name!, meteodir 
  use chem_param, only : fscale

  implicit none

  private

  public :: get_flightdata

  ! private 

  logical                     :: end_file   = .false. 
  ! filer_open: signal file open for flight input 
  logical                     :: filer_open = .false.
  ! file_open: signal file open for output
  logical,dimension(nregions) :: file_open  = .false.
  ! funit0: base unit for writing formatted output
  integer, parameter          :: funit0 = 210
  integer,parameter           :: nf_trace =  1
  real,dimension(nf_trace)    :: rmf
  ! number of locations to be calculated for 1 model time
  integer,parameter           :: nsamples =  2
  integer,dimension(nf_trace) :: if_trace = (/ 1 /)
  integer,dimension(6)        :: idate_flight


contains



  subroutine get_flightdata(region,idate_f)
    !
    !
    !
    use global_data, only : region_dat
    use tracer_data, only : mass_dat
    use meteodata  , only : m_dat, phlb_dat
    use global_data, only : outdir

    ! input/output
    integer,intent(in)              :: region
    ! idate_f: date for which output required...
    integer,intent(in),dimension(6) :: idate_f

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

    real,dimension(0:lm(region)) :: presh
    integer   :: i,is,js,l,n,isn,jsn,ls,j
    real      :: flon,flat,fpres,ris,rjs,dxr,dyr,wcx,wcy,rls

    ! start

    if ( .not. filer_open ) then   ! open input file...
       !i = len_trim(meteodir)
       !open( unit=funit0+region, form = 'formatted', &
       !     file = meteodir(1:i)//'flight.data', status = 'OLD')
       open( unit=funit0+region, form = 'formatted', &
            file = trim(outdir)//'/flight.data', status = 'OLD')
       read(funit0+region,*) idate_flight
       print *,'get_flightdata: Initial idate_flight read as:',idate_flight
       filer_open = .true.
    end if

    if ( .not. file_open(region) ) then   ! open output file
       open(unit = funit0+region+nregions,form = 'formatted', &
            file = trim(outdir)//'/flight_'//region_name(region)//'.out', &
            status = 'unknown')
       file_open(region)  = .true.
       write(funit0+region+nregions,'(6i6)') idate_flight
    end if

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

    if ( end_file ) return

    do i = 1,6 
       if (idate_flight(i).ne.idate_f(i)) return
    enddo

    !pointers to global arrays...
    m    =>    m_dat(region)%data
    phlb => phlb_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

    dyr = dy/yref(region)
    dxr = dx/xref(region)
    do n = 1,nsamples
       read(funit0+region,*) flon,flat,fpres
       if ( (flon .ge. xbeg(region) .and. (flon .le. xend(region)) ) .and. &
            (flat .ge. ybeg(region) .and. (flat .le. yend(region)) ) ) then
          ris = (flon-float(xbeg(region)))/dxr + 1.0
          rjs = (flat-float(ybeg(region)))/dyr + 1.0
          !is,js is the box where we want the mixing ratio
          is  = int(ris)
          js  = int(rjs)
          !fraction from the center of the is-box  (-0.5---+0.5)
          ris = ris-is-0.5
          !idem js
          rjs = rjs-js-0.5
          if(ris.gt.0) then 
             isn = is+1      !the neighbour for pressure interpolation
          else
             isn = is-1
          endif
          if(rjs.gt.0) then 
             jsn = js+1      !the neighbour for pressure interpolation
          else
             jsn = js-1
          endif
          wcx = (1.-abs(ris))
          wcy = (1.-abs(rjs))
          ! interpolate the pressure to flight position...
          ls = 1   !layer
          do l=0,lm(region)
             presh(l) = wcx*wcy*          phlb(is,js,l+1)+ &
                  (1.-wcx)*wcy*      phlb(isn,js,l+1) + &
                  wcx*(1.-wcy)*      phlb(is,jsn,l+1) + &
                  (1.-wcx)*(1.-wcy)* phlb(isn,jsn,l+1)
          enddo
          do l=0,lm(region)   ! selects layer 
             if(presh(l).lt.fpres) exit
          enddo

          select case(l)
          case(0)
             print*,'get_flightdata: Warning..., flight pressure ',&
                  'is below the surface pressure'
             ls = 1
             rls = -0.5   !surface...
          case default
             ls = l       !the flight layer
             ! the off-set from the center of the layer (-0.5--->+0.5)
             ! (interpolation is in (m))
             rls = (presh(l-1)-fpres)/(presh(l-1)-presh(l)) - 0.5
          end select
          !from is,js,ls, ris,rjs,rls, determine the mixing ratio   ...
          do j=1,nf_trace
             i = if_trace(j)
             ! rm-value is obtained from rm + slopes. 
             ! slope = rxm = (rm*dX/dx *deltaX/2)
             rmf(j) = rm(is,js,ls,i) + 2.0*(ris*rxm(is,js,ls,i) + &
                  rjs*rym(is,js,ls,i) + &
                  rls*rzm(is,js,ls,i) )
             rmf(j) = rmf(j)/m(is,js,ls) *fscale(i) 
          enddo
       else  !station not in region
          rmf(:) = -1.0 
       endif
       write(funit0+region+nregions,*) flon,flat,fpres,rmf

    enddo   ! nsamples

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

    read(funit0+region,*,end=999) idate_flight   !read new time....
    write(funit0+region+nregions,'(6i6)') idate_flight
    close(funit0+region+nregions)
    return

999 continue
    end_file = .true.

  end subroutine get_flightdata


end module user_output_flight
