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

!### 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"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:     USER_OUTPUT_FLIGHT
!
! !DESCRIPTION: 
!\\
!\\
! !INTERFACE: 
!
MODULE USER_OUTPUT_FLIGHT
  !
  ! !USES:
  !
  use dims,       only : lm, dx, xref, dy, yref, xbeg, xend, ybeg, yend
  use dims,       only : nregions, region_name!, meteodir 
  use chem_param, only : fscale
  use PARTOOLS,   only : isRoot, PAR_REDUCE_ELEMENT, PAR_BROADCAST
  use GO,         only : GOL, goErr

  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: get_flightdata
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter, private :: mname='ParTools'
  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, rmf_final
  ! 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
  !
  ! !REVISION HISTORY: 
  !   10 Jul 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:   GET_FLIGHTDATA
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE GET_FLIGHTDATA( region, idate_f, status)
    !
    ! !USES:
    !
    use global_data,  only : region_dat
    use tracer_data,  only : mass_dat
    use meteodata  ,  only : m_dat, phlb_dat
    use global_data,  only : outdir
    use tm5_distgrid, only : dgrid, get_distgrid, update_halo
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(in)               :: region
    integer, intent(in), dimension(6) :: idate_f ! date for which output required
    integer, intent(out)              :: status
    !
    ! !REVISION HISTORY: 
    !   10 Jul 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    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, i1, i2, j1, j2
    real      :: flon,flat,fpres,ris,rjs,dxr,dyr,wcx,wcy,rls

    ! --- START
    
    if ( end_file ) return
    
    if (isroot) then 
      if ( .not. filer_open ) then   ! open input file...
        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.
      endif
    endif
    call PAR_BROADCAST(idate_flight, status)

    if (isroot) then
      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
    end if

    ! bounds on this PE for this region
    call get_distgrid( dgrid(region), i_strt=i1, i_stop=i2, j_strt=j1, j_stop=j2)

    ! Assume that the halo of phlb_dat has not updated yet
    call update_halo( dgrid(region), phlb_dat(region)%data, phlb_dat(region)%halo, status)
    
    
    ! ---------- Method ----------------
    !  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.

    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
    rxm  => mass_dat(region)%rxm
    rym  => mass_dat(region)%rym
    rzm  => mass_dat(region)%rzm

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

    do n = 1,nsamples

      if (isroot) then 
        read(funit0+region,*) flon,flat,fpres
      endif
      call PAR_BROADCAST(flon, status)
      IF_NOTOK_RETURN(status=1)
      call PAR_BROADCAST(flat, status)
      IF_NOTOK_RETURN(status=1)
      call PAR_BROADCAST(fpres, status)
      IF_NOTOK_RETURN(status=1)
      
      rmf=-1.0
       
       ! if in region
       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)

          ! if on current PE
          if (is>=i1.and.is<=i2.and.js>=j1.and.js<=j2) then
          
             !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
           endif
       end if

       ! Barrier/gather on root
       call PAR_REDUCE_ELEMENT(rmf, 'MAX', rmf_final, status)
       IF_NOTOK_RETURN(status=1)
       
       if (isroot) then
         write(funit0+region+nregions,*) flon,flat,fpres,rmf_final
       endif
       
    enddo   ! nsamples

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

    ! try to read another date in the input file
    if (isroot) then
      read(funit0+region,*,end=999) idate_flight

      write(funit0+region+nregions,'(6i6)') idate_flight
      end_file = .false.      
      goto 1000
      
999   close(funit0+region+nregions)
      end_file = .true.
      
1000  continue
    endif

    call PAR_BROADCAST(end_file,status)
    IF_NOTOK_RETURN(status=1)
    status=0
    
  end subroutine get_flightdata
  !EOC


END MODULE USER_OUTPUT_FLIGHT