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