WebSVN – lagranto.wrf – Diff – /branches/met.nps.monterey/trace/trace.f90


 
PROGRAM trace_label_circle
 
  ! ********************************************************************
  ! *                                                                  *
  ! * Trace fields along trajectories                                  *
  ! *                                                                  *
  ! * Heini Wernli       first version:       April 1993               *
  ! * Michael Sprenger   major upgrade:       2008-2009                *
  ! * Sebastian Schemm WRF                                             *
  ! ********************************************************************
 
  implicit none
 
  ! --------------------------------------------------------------------
  ! Declaration of parameters
  ! --------------------------------------------------------------------
 
  ! Maximum number of levels for input files
  integer :: nlevmax
  parameter (nlevmax=100)
 
  ! Maximum number of input files (dates, length of trajectories)
  integer :: ndatmax
  parameter (ndatmax=500)
 
  ! Numerical epsilon (for float comparison)
  real :: eps
  parameter (eps=0.001)
 
  ! Prefix for primary and secondary fields
  character :: charp
  character :: chars
  parameter (charp='P')
  parameter (chars='S')
 
  ! --------------------------------------------------------------------
  ! Declaration of variables
  ! --------------------------------------------------------------------
 
  ! Input and output format for trajectories (see iotra.f)
  integer :: inpmode
  integer :: outmode
 
  ! Input parameters
  character(len=80) :: inpfile         ! Input trajectory file
  character(len=80) :: outfile         ! Output trajectory file
  integer :: ntra                      ! Number of trajectories
  integer :: ncol                      ! Number of columns (including time, lon, lat, p)
  integer :: ntim                      ! Number of times per trajectory
  integer :: ntrace0                   ! Number of trace variables
  character(len=80) :: tvar0(200)      ! Tracing variable (with mode specification)
  character(len=80) :: tvar(200)       ! Tracing variable name (only the variable)
  character(len=1)  :: tfil(200)       ! Filename prefix
  real :: fac(200)                     ! Scaling factor
  integer :: compfl(200)               ! Computation flag (1=compute)
  integer :: numdat                    ! Number of input files
  character(len=11) :: dat(ndatmax)    ! Dates of input files
  real :: timeinc                      ! Time increment between input files
  real :: tst                          ! Time shift of start relative to first data file
  real :: ten                          ! Time shift of end relatiev to first data file
  character(len=20) :: startdate       ! First time/date on trajectory
  character(len=20) :: enddate         ! Last time/date on trajectory
  character(len=80) :: timecheck       ! Either 'yes' or 'no'
  character(len=80) :: intmode         ! Interpolation mode ('normal', 'nearest')
 
  ! Trajectories
  real,allocatable, dimension (:,:,:) :: trainp          ! Input trajectories (ntra,ntim,ncol)
  real,allocatable, dimension (:,:,:) :: traint          ! Internal trajectories (ntra,ntim,ncol+ntrace0)
  real,allocatable, dimension (:,:,:) :: traout          ! Output trajectories (ntra,ntim,ncol+ntrace0)
  integer :: reftime(6)                                  ! Reference date
  character(len=80) :: varsinp(100)                      ! Field names for input trajectory
  character(len=80) :: varsint(100)                      ! Field names for internal trajectory
  character(len=80) :: varsout(100)                      ! Field names for output trajectory
  integer :: fid,fod                                     ! File identifier for inp and out trajectories
  real :: x0,y0,p0                                       ! Position of air parcel (physical space)
  real :: reltpos0                                       ! Relative time of air parcel
  real :: xind,yind,pind                                 ! Position of air parcel (grid space)
  integer :: fbflag                                      ! Flag for forward (1) or backward (-1) trajectories
  integer :: fok(100)                                    ! Flag whether field is ready
 
  ! Meteorological fields
  real,allocatable, dimension (:)     :: spt0,spt1       ! Surface pressure
  real,allocatable, dimension (:)     :: p3t0,p3t1       ! 3d-pressure
  real,allocatable, dimension (:)     :: f3t0,f3t1       ! 3d field for tracing
  character(len=80) :: svars(100)                        ! List of variables on S file
  character(len=80) :: pvars(100)                        ! List of variables on P file
  integer :: n_svars                                     ! Number of variables on S file
  integer :: n_pvars                                     ! Number of variables on P file
 
  ! Grid description
  real :: pollon,pollat   ! Longitude/latitude of pole
  real :: xmin,xmax       ! Zonal grid extension
  real :: ymin,ymax       ! Meridional grid extension
  integer :: nx,ny,nz     ! Grid dimensions
  real :: dx,dy           ! Horizontal grid resolution
  integer :: hem          ! Flag for hemispheric domain
  integer :: per          ! Flag for periodic domain
  real :: stagz           ! Vertical staggering
  real :: mdv             ! Missing data value
 
  ! Auxiliary variables
  integer :: i,j,k,l,m,n
  real :: rd
  character(len=80) :: filename,varname
  real :: time0,time1,reltpos
  integer :: itime0,itime1
  integer :: stat
  real :: tstart
  integer :: iloaded0,iloaded1
  real :: f0
  real :: frac
  real :: tload,tfrac
  integer :: isok
  character :: ch
  integer :: ind
  integer :: ind1,ind2,ind3,ind4,ind5
  integer :: ind6,ind7,ind8,ind9,ind0
  integer :: noutside
  real :: delta
  integer :: itrace0
  character(len=80) :: string
 
  ! Externals
  real :: int_index4
  external int_index4
  real :: sdis   
  external sdis 
 
  ! --------------------------------------------------------------------
  ! Start of program, Read parameters, get grid parameters
  ! --------------------------------------------------------------------
 
  ! Write start message
  print*,'========================================================='
  print*,'              *** START OF PROGRAM TRACE ***'
  print*
 
  ! Read parameters
  open(10,file='trace.param')
   read(10,*) inpfile
   read(10,*) outfile
   read(10,*) startdate
   read(10,*) enddate
   read(10,*) fbflag
   read(10,*) numdat
   if ( fbflag.eq.1) then
      do i=1,numdat
         read(10,'(a11)') dat(i)
      enddo
   else
      do i=numdat,1,-1
         read(10,'(a11)') dat(i)
      enddo
   endif
   read(10,*) timeinc
   read(10,*) tst
   read(10,*) ten
   read(10,*) ntra
   read(10,*) ntim
   read(10,*) ncol
   read(10,*) ntrace0
   do i=1,ntrace0
      read(10,*) tvar(i), fac(i), compfl(i), tfil(i)
   enddo
   read(10,*) n_pvars
   do i=1,n_pvars
      read(10,*) pvars(i)
   enddo
   read(10,*) n_svars
   do i=1,n_svars
      read(10,*) svars(i)
   enddo
   read(10,*) timecheck
   read(10,*) intmode
  close(10)
 
  ! Remove commented tracing fields
  itrace0 = 1
  do while ( itrace0.le.ntrace0)
     string = tvar(itrace0)
     if ( string(1:1).eq.'#' ) then
        do i=itrace0,ntrace0-1
           tvar(i)   = tvar(i+1)
           fac(i)    = fac(i+1)
           compfl(i) = compfl(i+1)
           tfil(i)   = tfil(i+1)
        enddo
        ntrace0 = ntrace0 - 1
     else
        itrace0 = itrace0 + 1
     endif
  enddo
 
  ! Save the tracing variable (including all mode specifications)
  do i=1,ntrace0
     tvar0(i) = tvar(i)
  enddo
 
  ! Set the formats of the input and output files
  call mode_tra(inpmode,inpfile)
  if (inpmode.eq.-1) inpmode=1
  call mode_tra(outmode,outfile)
  if (outmode.eq.-1) outmode=1
 
  ! Convert time shifts <tst,ten> from <hh.mm> into fractional time
  call hhmm2frac(tst,frac)
  tst = frac
  call hhmm2frac(ten,frac)
  ten = frac
 
  ! Set the time for the first data file (depending on forward/backward mode)
  if (fbflag.eq.1) then
    tstart = -tst
  else
    tstart = tst
  endif
 
  ! Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)
  ! The negative <-fid> of the file identifier is used as a flag for parameter retrieval
  filename = charp//dat(1)
  varname  = 'U'
  call input_open (fid,filename)
  call input_grid_wrf(fid,xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz)
  call input_close(fid)
 
  ! Allocate memory for some meteorological arrays
  allocate(spt0(nx*ny),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array spt0 ***'   ! Surface pressure
  allocate(spt1(nx*ny),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array spt1 ***'
  allocate(p3t0(nx*ny*nz),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array p3t0 ***'   ! Pressure
  allocate(p3t1(nx*ny*nz),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array p3t1 ***'
  allocate(f3t0(nx*ny*nz),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array p3t0 ***'   ! Tracing field
  allocate(f3t1(nx*ny*nz),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array p3t1 ***'
 
  ! Get memory for trajectory arrays
  allocate(trainp(ntra,ntim,ncol),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array tra      ***'
 
  ! Write some status information
  print*,'---- INPUT PARAMETERS -----------------------------------'
  print*
  print*,'  Input trajectory file  : ',trim(inpfile)
  print*,'  Output trajectory file : ',trim(outfile)
  print*,'  Format of input file   : ',inpmode
  print*,'  Format of output file  : ',outmode
  print*,'  Forward/backward       : ',fbflag
  print*,'  #tra                   : ',ntra
  print*,'  #col                   : ',ncol
  print*,'  #tim                   : ',ntim
  print*,'  No time check          : ',trim(timecheck)
  print*,'  Interpolation mode     : ',trim(intmode)
  do i=1,ntrace0
     if (compfl(i).eq.0) then
        print*,'  Tracing field          : ', trim(tvar(i)), fac(i), ' 0 ', tfil(i)
     else
        print*,'  Tracing field          : ', trim(tvar(i)), fac(i), '  1 ', tfil(i)
     endif
  enddo
  print*
  print*,'---- INPUT DATA FILES -----------------------------------'
  print*
  call frac2hhmm(tstart,tload)
  print*,'  Time of 1st data file  : ',tload
  print*,'  #input files           : ',numdat
  print*,'  time increment         : ',timeinc
  call frac2hhmm(tst,tload)
  print*,'  Shift of start         : ',tload
  call frac2hhmm(ten,tload)
  print*,'  Shift of end           : ',tload
  print*,'  First/last input file  : ',trim(dat(1)), ' ... ', trim(dat(numdat))
  print*,'  Primary variables      : ',trim(pvars(1))
  do i=2,n_pvars
     print*,'                         : ',trim(pvars(i))
  enddo
  if ( n_svars.ge.1 ) then
     print*,'  Secondary variables    : ',trim(svars(1))
     do i=2,n_svars
        print*,'                         : ',trim(svars(i))
     enddo
  endif
  print*
  print*,'---- CONSTANT GRID PARAMETERS ---------------------------'
  print*
  print*,'  xmin,xmax     : ',xmin,xmax
  print*,'  ymin,ymax     : ',ymin,ymax
  print*,'  dx,dy         : ',dx,dy
  print*,'  pollon,pollat : ',pollon,pollat
  print*,'  nx,ny,nz      : ',nx,ny,nz
  print*,'  per, hem      : ',per,hem
  print*
 
  ! --------------------------------------------------------------------
  ! Load the input trajectories
  ! --------------------------------------------------------------------
 
  ! Read the input trajectory file
  call ropen_tra(fid,inpfile,ntra,ntim,ncol,reftime,varsinp,inpmode)
  call read_tra (fid,trainp,ntra,ntim,ncol,inpmode)
  call close_tra(fid,inpmode)
 
  ! Check that first four columns correspond to time,x,y,z
  if ( (varsinp(1).ne.'time' ).or. &
       (varsinp(2).ne.'x'    ).or. &
       (varsinp(3).ne.'y'    ).or. &
       (varsinp(4).ne.'z'   ) )    &
  then
     print*,' ERROR: problem with input trajectories ...'
     stop
  endif
  varsinp(1) = 'time'
  varsinp(2) = 'x'
  varsinp(3) = 'y'
  varsinp(4) = 'z'
 
  ! Transform start coordinates from index space to WRF grid
  do i=1,ntra
     do j=1,ntim
         trainp(i,j,2) = xmin + ( trainp(i,j,2) - 1. ) * dx
         trainp(i,j,3) = ymin + ( trainp(i,j,3) - 1. ) * dy
     enddo
  enddo
 
  ! Write some status information of the input trajectories
  print*,'---- INPUT TRAJECTORIES ---------------------------------'
  print*
  print*,' Start date             : ',trim(startdate)
  print*,' End date               : ',trim(enddate)
  print*,' Reference time (year)  : ',reftime(1)
  print*,'                (month) : ',reftime(2)
  print*,'                (day)   : ',reftime(3)
  print*,'                (hour)  : ',reftime(4)
  print*,'                (min)   : ',reftime(5)
  print*,' Time range (min)       : ',reftime(6)
  do i=1,ncol
     print*,' Var                    :',i,trim(varsinp(i))
  enddo
  print*
 
  ! Check that first time is 0 - otherwise the tracing will produce
  ! wrong results because later in the code only absolute times are
  ! considered: <itime0   = int(abs(tfrac-tstart)/timeinc) + 1>. This
  ! will be changed in a future version.
  if ( abs( trainp(1,1,1) ).gt.eps ) then
     print*,' ERROR: First time of trajectory must be 0, i.e. '
     print*,'     correspond to the reference date. Otherwise'
     print*,'     the tracing will give wrong results... STOP'
     stop
  endif
 
  ! Save the full variable name (including shift specification)
  do i=1,ncol
     varsint(i)      = varsinp(i)
  enddo
  do i=1,ntrace0
     varsint(i+ncol) = tvar(i)
  enddo
 
  ! Split the variable name and set flags
  do i=1,ntrace0
 
     ! Set the prefix of the file name
     if ( compfl(i).eq.0  ) then
        tfil(i)=' '
        do j=1,n_pvars
           if ( tvar(i).eq.pvars(j) ) tfil(i)=charp
        enddo
        do j=1,n_svars
           if ( tvar(i).eq.svars(j) ) tfil(i)=chars
        enddo
     elseif ( tvar(i).eq.'Z' ) then
        tfil(i) = 'P'
     elseif ( tvar(i).eq.'ZB' ) then
        tfil(i) = 'P'
     else
        tfil(i) = '*' 
     endif
 
  enddo
 
  ! Write status information
  print*
  print*,'---- COMPLETE TABLE FOR TRACING -------------------------'
  print*
  do i=1,ntrace0
      write(*,'(i4,a4,a8,10x,8x,3x,a1,i5)') &
             i,' : ',trim(tvar(i)),tfil(i),compfl(i)
  enddo
 
  ! Allocate memory for internal trajectories
  allocate(traint(ntra,ntim,ncol+ntrace0),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array traint   ***'
 
  ! Copy input to output trajectory
  do i=1,ntra
     do j=1,ntim
        do k=1,ncol
           traint(i,j,k)=trainp(i,j,k)
        enddo
     enddo
  enddo
 
  ! --------------------------------------------------------------------
  ! Trace the fields (fields available on input files)
  ! --------------------------------------------------------------------
 
  print*
  print*,'---- TRACING FROM PRIMARY AND SECONDARY DATA FILES ------'
 
  ! Loop over all tracing fields
  do i=1,ntrace0
 
      ! Write some status information
      print*
      print*,' Now tracing             : ', trim(tvar(i)),' ',trim(tfil(i))
 
      ! Reset flags for load manager
      iloaded0 = -1
      iloaded1 = -1
 
      ! Reset the counter for fields outside domain
      noutside = 0
 
      ! Loop over all times
      do j=1,ntim
 
         ! Convert trajectory time from hh.mm to fractional time
         call hhmm2frac(trainp(1,j,1),tfrac)
 
         ! Get the times which are needed
         itime0   = int(abs(tfrac-tstart)/timeinc) + 1
         time0    = tstart + fbflag * real(itime0-1) * timeinc
         itime1   = itime0 + 1
         time1    = time0 + fbflag * timeinc
         if ( itime1.gt.numdat ) then
            itime1 = itime0
            time1  = time0
         endif
 
         ! Load manager: Check whether itime0 can be copied from itime1
         if ( itime0.eq.iloaded1 ) then
            f3t0     = f3t1
            p3t0     = p3t1
            spt0     = spt1
            iloaded0 = itime0
         endif
 
         ! Load manager: Check whether itime1 can be copied from itime0
         if ( itime1.eq.iloaded0 ) then
            f3t1     = f3t0
            p3t1     = p3t0
            spt1     = spt0
            iloaded1 = itime1
         endif
 
 
         ! Load manager:  Load first time (tracing variable and grid)
         if ( itime0.ne.iloaded0 ) then
 
            filename = tfil(i)//dat(itime0)
            call frac2hhmm(time0,tload)
            varname  = tvar(i)
            write(*,'(a23,a20,a3,a5,f7.2)') '    ->  loading          : ', trim(filename),'  ',trim(varname),tload
            call input_open (fid,filename)
 
	        call input_var_wrf(fid,varname,f3t0,nx,ny,nz)	! field to trace
 
	        varname='geopot'					            ! geopot.height
	        call input_var_wrf(fid,varname,p3t0,nx,ny,nz)
 
	        varname='geopots'					            ! surface geopot.height
	        call input_var_wrf(fid,varname,spt0,nx,ny,1 )
 
            call input_close(fid)
 
            iloaded0 = itime0
 
         endif
 
         ! Load manager: Load second time (tracing variable and grid)
         if ( itime1.ne.iloaded1 ) then
 
            filename = tfil(i)//dat(itime1)
            call frac2hhmm(time1,tload)
            varname  = tvar(i)
            write(*,'(a23,a20,a3,a5,f7.2)') '    ->  loading          : ', trim(filename),'  ',trim(varname),tload
            call input_open (fid,filename)
 
            call input_var_wrf(fid,varname,f3t1,nx,ny,nz)	        ! field to trace
 
	        varname='geopot'					                    ! geopot.height
	        call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
 
	        varname='geopots'					                    ! geopot.height
	        call input_var_wrf(fid,varname,spt1,nx,ny,1)
 
            call input_close(fid)
 
            iloaded1 = itime1
 
         endif
 
         ! Loop over all trajectories
         do k=1,ntra
 
            ! Set the horizontal position where to interpolate to
            x0       = traint(k,j,2)                          ! Longitude
            y0       = traint(k,j,3)                          ! Latitude
 
            ! Set the vertical position where to interpolate to: 
            if ( nz.gt.1 ) then
               p0 = traint(k,j,4)                            ! Pressure (3D tracing)
            else
               p0 = 1                                        ! Lowest level (2D tracing)
            endif
 
            ! Get the index where to interpolate (x0,y0,p0)
            if ( (abs(x0-mdv).gt.eps).and. (abs(y0-mdv).gt.eps) ) then
               call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0, &
                    p3t0,p3t1,spt0,spt1,3,nx,ny,nz,xmin,ymin,dx,dy,mdv)
            else
               xind = mdv
               yind = mdv
               pind = mdv
            endif
 
 
            ! Check if point is within grid (keep indices if ok)
            if ( (xind.ge.1.).and.(xind.le.real(nx)).and. &
                 (yind.ge.1.).and.(yind.le.real(ny)).and. &
                 (pind.ge.1.).and.(pind.le.real(nz)) )    &
            then
		           xind = xind
		           yind = yind
		           pind = pind
 
            endif
 
            ! ------------------------ NEAREST mode ------------------------------- 
 
            if ( intmode.eq.'nearest') then
 
               xind = real( nint(xind) )
               yind = real( nint(yind) )
               pind = real( nint(pind) )
 
               if ( xind.lt.1.  ) xind = 1.
               if ( xind.gt.nx  ) xind = real(nx)
               if ( yind.lt.1.  ) yind = 1.
               if ( yind.gt.ny  ) yind = real(ny)
 
               if ( pind.lt.1.  ) pind = 1.
               if ( pind.gt.nz  ) pind = real(nz)
 
               if ( abs(reltpos0).ge.eps ) then
                  print*,'ERROR (nearest): reltpos != 0',reltpos0
                  stop
               endif
 
               ! interpolate
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
 
            ! ------------------------ NORMAL mode -------------------------------
            else 
 
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
 
            endif
 
            ! Save the new field
            if ( abs(f0-mdv).gt.eps) then
               traint(k,j,ncol+i) = f0
            else
               traint(k,j,ncol+i) = mdv
            endif
 
         enddo ! end loop over all trajectories
 
      enddo ! end loop over all times
 
      ! Exit point for loop over all tracing variables
 110   continue
 
  enddo ! end loop over all tracing varoables
 
  ! --------------------------------------------------------------------
  ! Write output to output trajectory file
  ! --------------------------------------------------------------------
 
  ! Write status information
  print*
  print*,'---- WRITE OUTPUT TRAJECTORIES --------------------------'
  print*
 
  ! Allocate memory for internal trajectories
  allocate(traout(ntra,ntim,ncol+ntrace0),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array traout   ***'
 
  ! Copy input to output trajectory (apply scaling of output)
  do i=1,ntra
     do j=1,ntim
        do k=1,ncol+ntrace0
           if ( k.le.ncol ) then
              traout(i,j,k) = traint(i,j,k)
           elseif ( abs(traint(i,j,k)-mdv).gt.eps ) then
              traout(i,j,k) = fac(k-ncol) * traint(i,j,k)
           else
              traout(i,j,k) = mdv
           endif
        enddo
     enddo
  enddo
 
  ! Set the variable names for output trajectory
  do i=1,ncol+ntrace0
     varsout(i)      = varsint(i)
  enddo
 
 ! Transform coordinates into index space
   do i=1,ntra
     do j=1,ntim
         traout(i,j,2) = ( traout(i,j,2) - xmin ) / dx + 1.
         traout(i,j,3) = ( traout(i,j,3) - ymin ) / dy + 1.
     enddo
  enddo
 
  ! Write trajectories
  call wopen_tra(fod,outfile,ntra,ntim,ncol+ntrace0, &
       reftime,varsout,outmode)
  call write_tra(fod,traout ,ntra,ntim,ncol+ntrace0,outmode)
  call close_tra(fod,outmode)
 
  ! Write some status information, and end of program message
  print*
  print*,'---- STATUS INFORMATION --------------------------------'
  print*
  print*,' ok'
  print*
  print*,'              *** END OF PROGRAM TRACE ***'
  print*,'========================================================='
 
 
end program trace_label_circle
 

Rev 2	Rev 10
1		1
2	`PROGRAM trace_label_circle`	2	`PROGRAM trace_label_circle`
3		3
4	`! ********************************************************************`	4	`! ********************************************************************`
5	`! * *`	5	`! * *`
6	`! * Trace fields along trajectories *`	6	`! * Trace fields along trajectories *`
7	`! * *`	7	`! * *`
8	`! * Heini Wernli first version: April 1993 *`	8	`! * Heini Wernli first version: April 1993 *`
9	`! * Michael Sprenger major upgrade: 2008-2009 *`	9	`! * Michael Sprenger major upgrade: 2008-2009 *`
10	`! * Sebastian Schemm WRF *`	10	`! * Sebastian Schemm WRF *`
11	`! ********************************************************************`	11	`! ********************************************************************`
12		12
13	`implicit none`	13	`implicit none`
14		14
15	`! --------------------------------------------------------------------`	15	`! --------------------------------------------------------------------`
16	`! Declaration of parameters`	16	`! Declaration of parameters`
17	`! --------------------------------------------------------------------`	17	`! --------------------------------------------------------------------`
18		18
19	`! Maximum number of levels for input files`	19	`! Maximum number of levels for input files`
20	`integer :: nlevmax`	20	`integer :: nlevmax`
21	`parameter (nlevmax=100)`	21	`parameter (nlevmax=100)`
22		22
23	`! Maximum number of input files (dates, length of trajectories)`	23	`! Maximum number of input files (dates, length of trajectories)`
24	`integer :: ndatmax`	24	`integer :: ndatmax`
25	`parameter (ndatmax=500)`	25	`parameter (ndatmax=500)`
26		26
27	`! Numerical epsilon (for float comparison)`	27	`! Numerical epsilon (for float comparison)`
28	`real :: eps`	28	`real :: eps`
29	`parameter (eps=0.001)`	29	`parameter (eps=0.001)`
30		30
31	`! Prefix for primary and secondary fields`	31	`! Prefix for primary and secondary fields`
32	`character :: charp`	32	`character :: charp`
33	`character :: chars`	33	`character :: chars`
34	`parameter (charp='P')`	34	`parameter (charp='P')`
35	`parameter (chars='S')`	35	`parameter (chars='S')`
36		36
37	`! --------------------------------------------------------------------`	37	`! --------------------------------------------------------------------`
38	`! Declaration of variables`	38	`! Declaration of variables`
39	`! --------------------------------------------------------------------`	39	`! --------------------------------------------------------------------`
40		40
41	`! Input and output format for trajectories (see iotra.f)`	41	`! Input and output format for trajectories (see iotra.f)`
42	`integer :: inpmode`	42	`integer :: inpmode`
43	`integer :: outmode`	43	`integer :: outmode`
44		44
45	`! Input parameters`	45	`! Input parameters`
46	`character(len=80) :: inpfile ! Input trajectory file`	46	`character(len=80) :: inpfile ! Input trajectory file`
47	`character(len=80) :: outfile ! Output trajectory file`	47	`character(len=80) :: outfile ! Output trajectory file`
48	`integer :: ntra ! Number of trajectories`	48	`integer :: ntra ! Number of trajectories`
49	`integer :: ncol ! Number of columns (including time, lon, lat, p)`	49	`integer :: ncol ! Number of columns (including time, lon, lat, p)`
50	`integer :: ntim ! Number of times per trajectory`	50	`integer :: ntim ! Number of times per trajectory`
51	`integer :: ntrace0 ! Number of trace variables`	51	`integer :: ntrace0 ! Number of trace variables`
52	`character(len=80) :: tvar0(200) ! Tracing variable (with mode specification)`	52	`character(len=80) :: tvar0(200) ! Tracing variable (with mode specification)`
53	`character(len=80) :: tvar(200) ! Tracing variable name (only the variable)`	53	`character(len=80) :: tvar(200) ! Tracing variable name (only the variable)`
54	`character(len=1) :: tfil(200) ! Filename prefix`	54	`character(len=1) :: tfil(200) ! Filename prefix`
55	`real :: fac(200) ! Scaling factor`	55	`real :: fac(200) ! Scaling factor`
56	`integer :: compfl(200) ! Computation flag (1=compute)`	56	`integer :: compfl(200) ! Computation flag (1=compute)`
57	`integer :: numdat ! Number of input files`	57	`integer :: numdat ! Number of input files`
58	`character(len=11) :: dat(ndatmax) ! Dates of input files`	58	`character(len=11) :: dat(ndatmax) ! Dates of input files`
59	`real :: timeinc ! Time increment between input files`	59	`real :: timeinc ! Time increment between input files`
60	`real :: tst ! Time shift of start relative to first data file`	60	`real :: tst ! Time shift of start relative to first data file`
61	`real :: ten ! Time shift of end relatiev to first data file`	61	`real :: ten ! Time shift of end relatiev to first data file`
62	`character(len=20) :: startdate ! First time/date on trajectory`	62	`character(len=20) :: startdate ! First time/date on trajectory`
63	`character(len=20) :: enddate ! Last time/date on trajectory`	63	`character(len=20) :: enddate ! Last time/date on trajectory`
64	`character(len=80) :: timecheck ! Either 'yes' or 'no'`	64	`character(len=80) :: timecheck ! Either 'yes' or 'no'`
65	`character(len=80) :: intmode ! Interpolation mode ('normal', 'nearest')`	65	`character(len=80) :: intmode ! Interpolation mode ('normal', 'nearest')`
66		66
67	`! Trajectories`	67	`! Trajectories`
68	`real,allocatable, dimension (:,:,:) :: trainp ! Input trajectories (ntra,ntim,ncol)`	68	`real,allocatable, dimension (:,:,:) :: trainp ! Input trajectories (ntra,ntim,ncol)`
69	`real,allocatable, dimension (:,:,:) :: traint ! Internal trajectories (ntra,ntim,ncol+ntrace0)`	69	`real,allocatable, dimension (:,:,:) :: traint ! Internal trajectories (ntra,ntim,ncol+ntrace0)`
70	`real,allocatable, dimension (:,:,:) :: traout ! Output trajectories (ntra,ntim,ncol+ntrace0)`	70	`real,allocatable, dimension (:,:,:) :: traout ! Output trajectories (ntra,ntim,ncol+ntrace0)`
71	`integer :: reftime(6) ! Reference date`	71	`integer :: reftime(6) ! Reference date`
72	`character(len=80) :: varsinp(100) ! Field names for input trajectory`	72	`character(len=80) :: varsinp(100) ! Field names for input trajectory`
73	`character(len=80) :: varsint(100) ! Field names for internal trajectory`	73	`character(len=80) :: varsint(100) ! Field names for internal trajectory`
74	`character(len=80) :: varsout(100) ! Field names for output trajectory`	74	`character(len=80) :: varsout(100) ! Field names for output trajectory`
75	`integer :: fid,fod ! File identifier for inp and out trajectories`	75	`integer :: fid,fod ! File identifier for inp and out trajectories`
76	`real :: x0,y0,p0 ! Position of air parcel (physical space)`	76	`real :: x0,y0,p0 ! Position of air parcel (physical space)`
77	`real :: reltpos0 ! Relative time of air parcel`	77	`real :: reltpos0 ! Relative time of air parcel`
78	`real :: xind,yind,pind ! Position of air parcel (grid space)`	78	`real :: xind,yind,pind ! Position of air parcel (grid space)`
79	`integer :: fbflag ! Flag for forward (1) or backward (-1) trajectories`	79	`integer :: fbflag ! Flag for forward (1) or backward (-1) trajectories`
80	`integer :: fok(100) ! Flag whether field is ready`	80	`integer :: fok(100) ! Flag whether field is ready`
81		81
82	`! Meteorological fields`	82	`! Meteorological fields`
83	`real,allocatable, dimension (:) :: spt0,spt1 ! Surface pressure`	83	`real,allocatable, dimension (:) :: spt0,spt1 ! Surface pressure`
84	`real,allocatable, dimension (:) :: p3t0,p3t1 ! 3d-pressure`	84	`real,allocatable, dimension (:) :: p3t0,p3t1 ! 3d-pressure`
85	`real,allocatable, dimension (:) :: f3t0,f3t1 ! 3d field for tracing`	85	`real,allocatable, dimension (:) :: f3t0,f3t1 ! 3d field for tracing`
86	`character(len=80) :: svars(100) ! List of variables on S file`	86	`character(len=80) :: svars(100) ! List of variables on S file`
87	`character(len=80) :: pvars(100) ! List of variables on P file`	87	`character(len=80) :: pvars(100) ! List of variables on P file`
88	`integer :: n_svars ! Number of variables on S file`	88	`integer :: n_svars ! Number of variables on S file`
89	`integer :: n_pvars ! Number of variables on P file`	89	`integer :: n_pvars ! Number of variables on P file`
90		90
91	`! Grid description`	91	`! Grid description`
92	`real :: pollon,pollat ! Longitude/latitude of pole`	92	`real :: pollon,pollat ! Longitude/latitude of pole`
93	`real :: xmin,xmax ! Zonal grid extension`	93	`real :: xmin,xmax ! Zonal grid extension`
94	`real :: ymin,ymax ! Meridional grid extension`	94	`real :: ymin,ymax ! Meridional grid extension`
95	`integer :: nx,ny,nz ! Grid dimensions`	95	`integer :: nx,ny,nz ! Grid dimensions`
96	`real :: dx,dy ! Horizontal grid resolution`	96	`real :: dx,dy ! Horizontal grid resolution`
97	`integer :: hem ! Flag for hemispheric domain`	97	`integer :: hem ! Flag for hemispheric domain`
98	`integer :: per ! Flag for periodic domain`	98	`integer :: per ! Flag for periodic domain`
99	`real :: stagz ! Vertical staggering`	99	`real :: stagz ! Vertical staggering`
100	`real :: mdv ! Missing data value`	100	`real :: mdv ! Missing data value`
101		101
102	`! Auxiliary variables`	102	`! Auxiliary variables`
103	`integer :: i,j,k,l,m,n`	103	`integer :: i,j,k,l,m,n`
104	`real :: rd`	104	`real :: rd`
105	`character(len=80) :: filename,varname`	105	`character(len=80) :: filename,varname`
106	`real :: time0,time1,reltpos`	106	`real :: time0,time1,reltpos`
107	`integer :: itime0,itime1`	107	`integer :: itime0,itime1`
108	`integer :: stat`	108	`integer :: stat`
109	`real :: tstart`	109	`real :: tstart`
110	`integer :: iloaded0,iloaded1`	110	`integer :: iloaded0,iloaded1`
111	`real :: f0`	111	`real :: f0`
112	`real :: frac`	112	`real :: frac`
113	`real :: tload,tfrac`	113	`real :: tload,tfrac`
114	`integer :: isok`	114	`integer :: isok`
115	`character :: ch`	115	`character :: ch`
116	`integer :: ind`	116	`integer :: ind`
117	`integer :: ind1,ind2,ind3,ind4,ind5`	117	`integer :: ind1,ind2,ind3,ind4,ind5`
118	`integer :: ind6,ind7,ind8,ind9,ind0`	118	`integer :: ind6,ind7,ind8,ind9,ind0`
119	`integer :: noutside`	119	`integer :: noutside`
120	`real :: delta`	120	`real :: delta`
121	`integer :: itrace0`	121	`integer :: itrace0`
122	`character(len=80) :: string`	122	`character(len=80) :: string`
123		123
124	`! Externals`	124	`! Externals`
125	`real :: int_index4`	125	`real :: int_index4`
126	`external int_index4`	126	`external int_index4`
127	`real :: sdis`	127	`real :: sdis`
128	`external sdis`	128	`external sdis`
129		129
130	`! --------------------------------------------------------------------`	130	`! --------------------------------------------------------------------`
131	`! Start of program, Read parameters, get grid parameters`	131	`! Start of program, Read parameters, get grid parameters`
132	`! --------------------------------------------------------------------`	132	`! --------------------------------------------------------------------`
133		133
134	`! Write start message`	134	`! Write start message`
135	`print*,'========================================================='`	135	`print*,'========================================================='`
136	`print,' START OF PROGRAM TRACE *'`	136	`print,' START OF PROGRAM TRACE *'`
137	`print*`	137	`print*`
138		138
139	`! Read parameters`	139	`! Read parameters`
140	`open(10,file='trace.param')`	140	`open(10,file='trace.param')`
141	`read(10,*) inpfile`	141	`read(10,*) inpfile`
142	`read(10,*) outfile`	142	`read(10,*) outfile`
143	`read(10,*) startdate`	143	`read(10,*) startdate`
144	`read(10,*) enddate`	144	`read(10,*) enddate`
145	`read(10,*) fbflag`	145	`read(10,*) fbflag`
146	`read(10,*) numdat`	146	`read(10,*) numdat`
147	`if ( fbflag.eq.1) then`	147	`if ( fbflag.eq.1) then`
148	`do i=1,numdat`	148	`do i=1,numdat`
149	`read(10,'(a11)') dat(i)`	149	`read(10,'(a11)') dat(i)`
150	`enddo`	150	`enddo`
151	`else`	151	`else`
152	`do i=numdat,1,-1`	152	`do i=numdat,1,-1`
153	`read(10,'(a11)') dat(i)`	153	`read(10,'(a11)') dat(i)`
154	`enddo`	154	`enddo`
155	`endif`	155	`endif`
156	`read(10,*) timeinc`	156	`read(10,*) timeinc`
157	`read(10,*) tst`	157	`read(10,*) tst`
158	`read(10,*) ten`	158	`read(10,*) ten`
159	`read(10,*) ntra`	159	`read(10,*) ntra`
160	`read(10,*) ntim`	160	`read(10,*) ntim`
161	`read(10,*) ncol`	161	`read(10,*) ncol`
162	`read(10,*) ntrace0`	162	`read(10,*) ntrace0`
163	`do i=1,ntrace0`	163	`do i=1,ntrace0`
164	`read(10,*) tvar(i), fac(i), compfl(i), tfil(i)`	164	`read(10,*) tvar(i), fac(i), compfl(i), tfil(i)`
165	`enddo`	165	`enddo`
166	`read(10,*) n_pvars`	166	`read(10,*) n_pvars`
167	`do i=1,n_pvars`	167	`do i=1,n_pvars`
168	`read(10,*) pvars(i)`	168	`read(10,*) pvars(i)`
169	`enddo`	169	`enddo`
170	`read(10,*) n_svars`	170	`read(10,*) n_svars`
171	`do i=1,n_svars`	171	`do i=1,n_svars`
172	`read(10,*) svars(i)`	172	`read(10,*) svars(i)`
173	`enddo`	173	`enddo`
174	`read(10,*) timecheck`	174	`read(10,*) timecheck`
175	`read(10,*) intmode`	175	`read(10,*) intmode`
176	`close(10)`	176	`close(10)`
177		177
178	`! Remove commented tracing fields`	178	`! Remove commented tracing fields`
179	`itrace0 = 1`	179	`itrace0 = 1`
180	`do while ( itrace0.le.ntrace0)`	180	`do while ( itrace0.le.ntrace0)`
181	`string = tvar(itrace0)`	181	`string = tvar(itrace0)`
182	`if ( string(1:1).eq.'#' ) then`	182	`if ( string(1:1).eq.'#' ) then`
183	`do i=itrace0,ntrace0-1`	183	`do i=itrace0,ntrace0-1`
184	`tvar(i) = tvar(i+1)`	184	`tvar(i) = tvar(i+1)`
185	`fac(i) = fac(i+1)`	185	`fac(i) = fac(i+1)`
186	`compfl(i) = compfl(i+1)`	186	`compfl(i) = compfl(i+1)`
187	`tfil(i) = tfil(i+1)`	187	`tfil(i) = tfil(i+1)`
188	`enddo`	188	`enddo`
189	`ntrace0 = ntrace0 - 1`	189	`ntrace0 = ntrace0 - 1`
190	`else`	190	`else`
191	`itrace0 = itrace0 + 1`	191	`itrace0 = itrace0 + 1`
192	`endif`	192	`endif`
193	`enddo`	193	`enddo`
194		194
195	`! Save the tracing variable (including all mode specifications)`	195	`! Save the tracing variable (including all mode specifications)`
196	`do i=1,ntrace0`	196	`do i=1,ntrace0`
197	`tvar0(i) = tvar(i)`	197	`tvar0(i) = tvar(i)`
198	`enddo`	198	`enddo`
199		199
200	`! Set the formats of the input and output files`	200	`! Set the formats of the input and output files`
201	`call mode_tra(inpmode,inpfile)`	201	`call mode_tra(inpmode,inpfile)`
202	`if (inpmode.eq.-1) inpmode=1`	202	`if (inpmode.eq.-1) inpmode=1`
203	`call mode_tra(outmode,outfile)`	203	`call mode_tra(outmode,outfile)`
204	`if (outmode.eq.-1) outmode=1`	204	`if (outmode.eq.-1) outmode=1`
205		205
206	`! Convert time shifts <tst,ten> from <hh.mm> into fractional time`	206	`! Convert time shifts <tst,ten> from <hh.mm> into fractional time`
207	`call hhmm2frac(tst,frac)`	207	`call hhmm2frac(tst,frac)`
208	`tst = frac`	208	`tst = frac`
209	`call hhmm2frac(ten,frac)`	209	`call hhmm2frac(ten,frac)`
210	`ten = frac`	210	`ten = frac`
211		211
212	`! Set the time for the first data file (depending on forward/backward mode)`	212	`! Set the time for the first data file (depending on forward/backward mode)`
213	`if (fbflag.eq.1) then`	213	`if (fbflag.eq.1) then`
214	`tstart = -tst`	214	`tstart = -tst`
215	`else`	215	`else`
216	`tstart = tst`	216	`tstart = tst`
217	`endif`	217	`endif`
218		218
219	`! Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)`	219	`! Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)`
220	`! The negative <-fid> of the file identifier is used as a flag for parameter retrieval`	220	`! The negative <-fid> of the file identifier is used as a flag for parameter retrieval`
221	`filename = charp//dat(1)`	221	`filename = charp//dat(1)`
222	`varname = 'U'`	222	`varname = 'U'`
223	`call input_open (fid,filename)`	223	`call input_open (fid,filename)`
224	`call input_grid_wrf(fid,xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz)`	224	`call input_grid_wrf(fid,xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz)`
225	`call input_close(fid)`	225	`call input_close(fid)`
226		226
227	`! Allocate memory for some meteorological arrays`	227	`! Allocate memory for some meteorological arrays`
228	`allocate(spt0(nx*ny),stat=stat)`	228	`allocate(spt0(nx*ny),stat=stat)`
229	`if (stat.ne.0) print,' error allocating array spt0 *' ! Surface pressure`	229	`if (stat.ne.0) print,' error allocating array spt0 *' ! Surface pressure`
230	`allocate(spt1(nx*ny),stat=stat)`	230	`allocate(spt1(nx*ny),stat=stat)`
231	`if (stat.ne.0) print,' error allocating array spt1 *'`	231	`if (stat.ne.0) print,' error allocating array spt1 *'`
232	`allocate(p3t0(nxnynz),stat=stat)`	232	`allocate(p3t0(nxnynz),stat=stat)`
233	`if (stat.ne.0) print,' error allocating array p3t0 *' ! Pressure`	233	`if (stat.ne.0) print,' error allocating array p3t0 *' ! Pressure`
234	`allocate(p3t1(nxnynz),stat=stat)`	234	`allocate(p3t1(nxnynz),stat=stat)`
235	`if (stat.ne.0) print,' error allocating array p3t1 *'`	235	`if (stat.ne.0) print,' error allocating array p3t1 *'`
236	`allocate(f3t0(nxnynz),stat=stat)`	236	`allocate(f3t0(nxnynz),stat=stat)`
237	`if (stat.ne.0) print,' error allocating array p3t0 *' ! Tracing field`	237	`if (stat.ne.0) print,' error allocating array p3t0 *' ! Tracing field`
238	`allocate(f3t1(nxnynz),stat=stat)`	238	`allocate(f3t1(nxnynz),stat=stat)`
239	`if (stat.ne.0) print,' error allocating array p3t1 *'`	239	`if (stat.ne.0) print,' error allocating array p3t1 *'`
240		240
241	`! Get memory for trajectory arrays`	241	`! Get memory for trajectory arrays`
242	`allocate(trainp(ntra,ntim,ncol),stat=stat)`	242	`allocate(trainp(ntra,ntim,ncol),stat=stat)`
243	`if (stat.ne.0) print,' error allocating array tra *'`	243	`if (stat.ne.0) print,' error allocating array tra *'`
244		244
245	`! Write some status information`	245	`! Write some status information`
246	`print*,'---- INPUT PARAMETERS -----------------------------------'`	246	`print*,'---- INPUT PARAMETERS -----------------------------------'`
247	`print*`	247	`print*`
248	`print*,' Input trajectory file : ',trim(inpfile)`	248	`print*,' Input trajectory file : ',trim(inpfile)`
249	`print*,' Output trajectory file : ',trim(outfile)`	249	`print*,' Output trajectory file : ',trim(outfile)`
250	`print*,' Format of input file : ',inpmode`	250	`print*,' Format of input file : ',inpmode`
251	`print*,' Format of output file : ',outmode`	251	`print*,' Format of output file : ',outmode`
252	`print*,' Forward/backward : ',fbflag`	252	`print*,' Forward/backward : ',fbflag`
253	`print*,' #tra : ',ntra`	253	`print*,' #tra : ',ntra`
254	`print*,' #col : ',ncol`	254	`print*,' #col : ',ncol`
255	`print*,' #tim : ',ntim`	255	`print*,' #tim : ',ntim`
256	`print*,' No time check : ',trim(timecheck)`	256	`print*,' No time check : ',trim(timecheck)`
257	`print*,' Interpolation mode : ',trim(intmode)`	257	`print*,' Interpolation mode : ',trim(intmode)`
258	`do i=1,ntrace0`	258	`do i=1,ntrace0`
259	`if (compfl(i).eq.0) then`	259	`if (compfl(i).eq.0) then`
260	`print*,' Tracing field : ', trim(tvar(i)), fac(i), ' 0 ', tfil(i)`	260	`print*,' Tracing field : ', trim(tvar(i)), fac(i), ' 0 ', tfil(i)`
261	`else`	261	`else`
262	`print*,' Tracing field : ', trim(tvar(i)), fac(i), ' 1 ', tfil(i)`	262	`print*,' Tracing field : ', trim(tvar(i)), fac(i), ' 1 ', tfil(i)`
263	`endif`	263	`endif`
264	`enddo`	264	`enddo`
265	`print*`	265	`print*`
266	`print*,'---- INPUT DATA FILES -----------------------------------'`	266	`print*,'---- INPUT DATA FILES -----------------------------------'`
267	`print*`	267	`print*`
268	`call frac2hhmm(tstart,tload)`	268	`call frac2hhmm(tstart,tload)`
269	`print*,' Time of 1st data file : ',tload`	269	`print*,' Time of 1st data file : ',tload`
270	`print*,' #input files : ',numdat`	270	`print*,' #input files : ',numdat`
271	`print*,' time increment : ',timeinc`	271	`print*,' time increment : ',timeinc`
272	`call frac2hhmm(tst,tload)`	272	`call frac2hhmm(tst,tload)`
273	`print*,' Shift of start : ',tload`	273	`print*,' Shift of start : ',tload`
274	`call frac2hhmm(ten,tload)`	274	`call frac2hhmm(ten,tload)`
275	`print*,' Shift of end : ',tload`	275	`print*,' Shift of end : ',tload`
276	`print*,' First/last input file : ',trim(dat(1)), ' ... ', trim(dat(numdat))`	276	`print*,' First/last input file : ',trim(dat(1)), ' ... ', trim(dat(numdat))`
277	`print*,' Primary variables : ',trim(pvars(1))`	277	`print*,' Primary variables : ',trim(pvars(1))`
278	`do i=2,n_pvars`	278	`do i=2,n_pvars`
279	`print*,' : ',trim(pvars(i))`	279	`print*,' : ',trim(pvars(i))`
280	`enddo`	280	`enddo`
281	`if ( n_svars.ge.1 ) then`	281	`if ( n_svars.ge.1 ) then`
282	`print*,' Secondary variables : ',trim(svars(1))`	282	`print*,' Secondary variables : ',trim(svars(1))`
283	`do i=2,n_svars`	283	`do i=2,n_svars`
284	`print*,' : ',trim(svars(i))`	284	`print*,' : ',trim(svars(i))`
285	`enddo`	285	`enddo`
286	`endif`	286	`endif`
287	`print*`	287	`print*`
288	`print*,'---- CONSTANT GRID PARAMETERS ---------------------------'`	288	`print*,'---- CONSTANT GRID PARAMETERS ---------------------------'`
289	`print*`	289	`print*`
290	`print*,' xmin,xmax : ',xmin,xmax`	290	`print*,' xmin,xmax : ',xmin,xmax`
291	`print*,' ymin,ymax : ',ymin,ymax`	291	`print*,' ymin,ymax : ',ymin,ymax`
292	`print*,' dx,dy : ',dx,dy`	292	`print*,' dx,dy : ',dx,dy`
293	`print*,' pollon,pollat : ',pollon,pollat`	293	`print*,' pollon,pollat : ',pollon,pollat`
294	`print*,' nx,ny,nz : ',nx,ny,nz`	294	`print*,' nx,ny,nz : ',nx,ny,nz`
295	`print*,' per, hem : ',per,hem`	295	`print*,' per, hem : ',per,hem`
296	`print*`	296	`print*`
297		297
298	`! --------------------------------------------------------------------`	298	`! --------------------------------------------------------------------`
299	`! Load the input trajectories`	299	`! Load the input trajectories`
300	`! --------------------------------------------------------------------`	300	`! --------------------------------------------------------------------`
301		301
302	`! Read the input trajectory file`	302	`! Read the input trajectory file`
303	`call ropen_tra(fid,inpfile,ntra,ntim,ncol,reftime,varsinp,inpmode)`	303	`call ropen_tra(fid,inpfile,ntra,ntim,ncol,reftime,varsinp,inpmode)`
304	`call read_tra (fid,trainp,ntra,ntim,ncol,inpmode)`	304	`call read_tra (fid,trainp,ntra,ntim,ncol,inpmode)`
305	`call close_tra(fid,inpmode)`	305	`call close_tra(fid,inpmode)`
306		306
307	`! Check that first four columns correspond to time,x,y,z`	307	`! Check that first four columns correspond to time,x,y,z`
308	`if ( (varsinp(1).ne.'time' ).or. &`	308	`if ( (varsinp(1).ne.'time' ).or. &`
309	`(varsinp(2).ne.'x' ).or. &`	309	`(varsinp(2).ne.'x' ).or. &`
310	`(varsinp(3).ne.'y' ).or. &`	310	`(varsinp(3).ne.'y' ).or. &`
311	`(varsinp(4).ne.'z' ) ) &`	311	`(varsinp(4).ne.'z' ) ) &`
312	`then`	312	`then`
313	`print*,' ERROR: problem with input trajectories ...'`	313	`print*,' ERROR: problem with input trajectories ...'`
314	`stop`	314	`stop`
315	`endif`	315	`endif`
316	`varsinp(1) = 'time'`	316	`varsinp(1) = 'time'`
317	`varsinp(2) = 'x'`	317	`varsinp(2) = 'x'`
318	`varsinp(3) = 'y'`	318	`varsinp(3) = 'y'`
319	`varsinp(4) = 'z'`	319	`varsinp(4) = 'z'`
320		320
321	`! Transform start coordinates from index space to WRF grid`	321	`! Transform start coordinates from index space to WRF grid`
322	`do i=1,ntra`	322	`do i=1,ntra`
323	`do j=1,ntim`	323	`do j=1,ntim`
324	`trainp(i,j,2) = xmin + ( trainp(i,j,2) - 1. ) * dx`	324	`trainp(i,j,2) = xmin + ( trainp(i,j,2) - 1. ) * dx`
325	`trainp(i,j,3) = ymin + ( trainp(i,j,3) - 1. ) * dy`	325	`trainp(i,j,3) = ymin + ( trainp(i,j,3) - 1. ) * dy`
326	`enddo`	326	`enddo`
327	`enddo`	327	`enddo`
328		328
329	`! Write some status information of the input trajectories`	329	`! Write some status information of the input trajectories`
330	`print*,'---- INPUT TRAJECTORIES ---------------------------------'`	330	`print*,'---- INPUT TRAJECTORIES ---------------------------------'`
331	`print*`	331	`print*`
332	`print*,' Start date : ',trim(startdate)`	332	`print*,' Start date : ',trim(startdate)`
333	`print*,' End date : ',trim(enddate)`	333	`print*,' End date : ',trim(enddate)`
334	`print*,' Reference time (year) : ',reftime(1)`	334	`print*,' Reference time (year) : ',reftime(1)`
335	`print*,' (month) : ',reftime(2)`	335	`print*,' (month) : ',reftime(2)`
336	`print*,' (day) : ',reftime(3)`	336	`print*,' (day) : ',reftime(3)`
337	`print*,' (hour) : ',reftime(4)`	337	`print*,' (hour) : ',reftime(4)`
338	`print*,' (min) : ',reftime(5)`	338	`print*,' (min) : ',reftime(5)`
339	`print*,' Time range (min) : ',reftime(6)`	339	`print*,' Time range (min) : ',reftime(6)`
340	`do i=1,ncol`	340	`do i=1,ncol`
341	`print*,' Var :',i,trim(varsinp(i))`	341	`print*,' Var :',i,trim(varsinp(i))`
342	`enddo`	342	`enddo`
343	`print*`	343	`print*`
344		344
345	`! Check that first time is 0 - otherwise the tracing will produce`	345	`! Check that first time is 0 - otherwise the tracing will produce`
346	`! wrong results because later in the code only absolute times are`	346	`! wrong results because later in the code only absolute times are`
347	`! considered: <itime0 = int(abs(tfrac-tstart)/timeinc) + 1>. This`	347	`! considered: <itime0 = int(abs(tfrac-tstart)/timeinc) + 1>. This`
348	`! will be changed in a future version.`	348	`! will be changed in a future version.`
349	`if ( abs( trainp(1,1,1) ).gt.eps ) then`	349	`if ( abs( trainp(1,1,1) ).gt.eps ) then`
350	`print*,' ERROR: First time of trajectory must be 0, i.e. '`	350	`print*,' ERROR: First time of trajectory must be 0, i.e. '`
351	`print*,' correspond to the reference date. Otherwise'`	351	`print*,' correspond to the reference date. Otherwise'`
352	`print*,' the tracing will give wrong results... STOP'`	352	`print*,' the tracing will give wrong results... STOP'`
353	`stop`	353	`stop`
354	`endif`	354	`endif`
355		355
356	`! Save the full variable name (including shift specification)`	356	`! Save the full variable name (including shift specification)`
357	`do i=1,ncol`	357	`do i=1,ncol`
358	`varsint(i) = varsinp(i)`	358	`varsint(i) = varsinp(i)`
359	`enddo`	359	`enddo`
360	`do i=1,ntrace0`	360	`do i=1,ntrace0`
361	`varsint(i+ncol) = tvar(i)`	361	`varsint(i+ncol) = tvar(i)`
362	`enddo`	362	`enddo`
363		363
364	`! Split the variable name and set flags`	364	`! Split the variable name and set flags`
365	`do i=1,ntrace0`	365	`do i=1,ntrace0`
366		366
367	`! Set the prefix of the file name`	367	`! Set the prefix of the file name`
368	`if ( compfl(i).eq.0 ) then`	368	`if ( compfl(i).eq.0 ) then`
369	`tfil(i)=' '`	369	`tfil(i)=' '`
370	`do j=1,n_pvars`	370	`do j=1,n_pvars`
371	`if ( tvar(i).eq.pvars(j) ) tfil(i)=charp`	371	`if ( tvar(i).eq.pvars(j) ) tfil(i)=charp`
372	`enddo`	372	`enddo`
373	`do j=1,n_svars`	373	`do j=1,n_svars`
374	`if ( tvar(i).eq.svars(j) ) tfil(i)=chars`	374	`if ( tvar(i).eq.svars(j) ) tfil(i)=chars`
375	`enddo`	375	`enddo`
376	`elseif ( tvar(i).eq.'Z' ) then`	376	`elseif ( tvar(i).eq.'Z' ) then`
377	`tfil(i) = 'P'`	377	`tfil(i) = 'P'`
378	`elseif ( tvar(i).eq.'ZB' ) then`	378	`elseif ( tvar(i).eq.'ZB' ) then`
379	`tfil(i) = 'P'`	379	`tfil(i) = 'P'`
380	`else`	380	`else`
381	`tfil(i) = '*'`	381	`tfil(i) = '*'`
382	`endif`	382	`endif`
383		383
384	`enddo`	384	`enddo`
385		385
386	`! Write status information`	386	`! Write status information`
387	`print*`	387	`print*`
388	`print*,'---- COMPLETE TABLE FOR TRACING -------------------------'`	388	`print*,'---- COMPLETE TABLE FOR TRACING -------------------------'`
389	`print*`	389	`print*`
390	`do i=1,ntrace0`	390	`do i=1,ntrace0`
391	`write(*,'(i4,a4,a8,10x,8x,3x,a1,i5)') &`	391	`write(*,'(i4,a4,a8,10x,8x,3x,a1,i5)') &`
392	`i,' : ',trim(tvar(i)),tfil(i),compfl(i)`	392	`i,' : ',trim(tvar(i)),tfil(i),compfl(i)`
393	`enddo`	393	`enddo`
394		394
395	`! Allocate memory for internal trajectories`	395	`! Allocate memory for internal trajectories`
396	`allocate(traint(ntra,ntim,ncol+ntrace0),stat=stat)`	396	`allocate(traint(ntra,ntim,ncol+ntrace0),stat=stat)`
397	`if (stat.ne.0) print,' error allocating array traint *'`	397	`if (stat.ne.0) print,' error allocating array traint *'`
398		398
399	`! Copy input to output trajectory`	399	`! Copy input to output trajectory`
400	`do i=1,ntra`	400	`do i=1,ntra`
401	`do j=1,ntim`	401	`do j=1,ntim`
402	`do k=1,ncol`	402	`do k=1,ncol`
403	`traint(i,j,k)=trainp(i,j,k)`	403	`traint(i,j,k)=trainp(i,j,k)`
404	`enddo`	404	`enddo`
405	`enddo`	405	`enddo`
406	`enddo`	406	`enddo`
407		407
408	`! --------------------------------------------------------------------`	408	`! --------------------------------------------------------------------`
409	`! Trace the fields (fields available on input files)`	409	`! Trace the fields (fields available on input files)`
410	`! --------------------------------------------------------------------`	410	`! --------------------------------------------------------------------`
411		411
412	`print*`	412	`print*`
413	`print*,'---- TRACING FROM PRIMARY AND SECONDARY DATA FILES ------'`	413	`print*,'---- TRACING FROM PRIMARY AND SECONDARY DATA FILES ------'`
414		414
415	`! Loop over all tracing fields`	415	`! Loop over all tracing fields`
416	`do i=1,ntrace0`	416	`do i=1,ntrace0`
417		417
418	`! Write some status information`	418	`! Write some status information`
419	`print*`	419	`print*`
420	`print*,' Now tracing : ', trim(tvar(i)),' ',trim(tfil(i))`	420	`print*,' Now tracing : ', trim(tvar(i)),' ',trim(tfil(i))`
421		421
422	`! Reset flags for load manager`	422	`! Reset flags for load manager`
423	`iloaded0 = -1`	423	`iloaded0 = -1`
424	`iloaded1 = -1`	424	`iloaded1 = -1`
425		425
426	`! Reset the counter for fields outside domain`	426	`! Reset the counter for fields outside domain`
427	`noutside = 0`	427	`noutside = 0`
428		428
429	`! Loop over all times`	429	`! Loop over all times`
430	`do j=1,ntim`	430	`do j=1,ntim`
431		431
432	`! Convert trajectory time from hh.mm to fractional time`	432	`! Convert trajectory time from hh.mm to fractional time`
433	`call hhmm2frac(trainp(1,j,1),tfrac)`	433	`call hhmm2frac(trainp(1,j,1),tfrac)`
434		434
435	`! Get the times which are needed`	435	`! Get the times which are needed`
436	`itime0 = int(abs(tfrac-tstart)/timeinc) + 1`	436	`itime0 = int(abs(tfrac-tstart)/timeinc) + 1`
437	`time0 = tstart + fbflag * real(itime0-1) * timeinc`	437	`time0 = tstart + fbflag * real(itime0-1) * timeinc`
438	`itime1 = itime0 + 1`	438	`itime1 = itime0 + 1`
439	`time1 = time0 + fbflag * timeinc`	439	`time1 = time0 + fbflag * timeinc`
440	`if ( itime1.gt.numdat ) then`	440	`if ( itime1.gt.numdat ) then`
441	`itime1 = itime0`	441	`itime1 = itime0`
442	`time1 = time0`	442	`time1 = time0`
443	`endif`	443	`endif`
444		444
445	`! Load manager: Check whether itime0 can be copied from itime1`	445	`! Load manager: Check whether itime0 can be copied from itime1`
446	`if ( itime0.eq.iloaded1 ) then`	446	`if ( itime0.eq.iloaded1 ) then`
447	`f3t0 = f3t1`	447	`f3t0 = f3t1`
448	`p3t0 = p3t1`	448	`p3t0 = p3t1`
449	`spt0 = spt1`	449	`spt0 = spt1`
450	`iloaded0 = itime0`	450	`iloaded0 = itime0`
451	`endif`	451	`endif`
452		452
453	`! Load manager: Check whether itime1 can be copied from itime0`	453	`! Load manager: Check whether itime1 can be copied from itime0`
454	`if ( itime1.eq.iloaded0 ) then`	454	`if ( itime1.eq.iloaded0 ) then`
455	`f3t1 = f3t0`	455	`f3t1 = f3t0`
456	`p3t1 = p3t0`	456	`p3t1 = p3t0`
457	`spt1 = spt0`	457	`spt1 = spt0`
458	`iloaded1 = itime1`	458	`iloaded1 = itime1`
459	`endif`	459	`endif`
460		460
461		461
462	`! Load manager: Load first time (tracing variable and grid)`	462	`! Load manager: Load first time (tracing variable and grid)`
463	`if ( itime0.ne.iloaded0 ) then`	463	`if ( itime0.ne.iloaded0 ) then`
464		464
465	`filename = tfil(i)//dat(itime0)`	465	`filename = tfil(i)//dat(itime0)`
466	`call frac2hhmm(time0,tload)`	466	`call frac2hhmm(time0,tload)`
467	`varname = tvar(i)`	467	`varname = tvar(i)`
468	`write(*,'(a23,a20,a3,a5,f7.2)') ' -> loading : ', trim(filename),' ',trim(varname),tload`	468	`write(*,'(a23,a20,a3,a5,f7.2)') ' -> loading : ', trim(filename),' ',trim(varname),tload`
469	`call input_open (fid,filename)`	469	`call input_open (fid,filename)`
470		470
471	`call input_var_wrf(fid,varname,f3t0,nx,ny,nz) ! field to trace`	471	`call input_var_wrf(fid,varname,f3t0,nx,ny,nz) ! field to trace`
472		472
473	`varname='geopot' ! geopot.height`	473	`varname='geopot' ! geopot.height`
474	`call input_var_wrf(fid,varname,p3t0,nx,ny,nz)`	474	`call input_var_wrf(fid,varname,p3t0,nx,ny,nz)`
475		475
476	`varname='geopots' ! surface geopot.height`	476	`varname='geopots' ! surface geopot.height`
477	`call input_var_wrf(fid,varname,spt0,nx,ny,1 )`	477	`call input_var_wrf(fid,varname,spt0,nx,ny,1 )`
478		478
479	`call input_close(fid)`	479	`call input_close(fid)`
480		480
481	`iloaded0 = itime0`	481	`iloaded0 = itime0`
482		482
483	`endif`	483	`endif`
484		484
485	`! Load manager: Load second time (tracing variable and grid)`	485	`! Load manager: Load second time (tracing variable and grid)`
486	`if ( itime1.ne.iloaded1 ) then`	486	`if ( itime1.ne.iloaded1 ) then`
487		487
488	`filename = tfil(i)//dat(itime1)`	488	`filename = tfil(i)//dat(itime1)`
489	`call frac2hhmm(time1,tload)`	489	`call frac2hhmm(time1,tload)`
490	`varname = tvar(i)`	490	`varname = tvar(i)`
491	`write(*,'(a23,a20,a3,a5,f7.2)') ' -> loading : ', trim(filename),' ',trim(varname),tload`	491	`write(*,'(a23,a20,a3,a5,f7.2)') ' -> loading : ', trim(filename),' ',trim(varname),tload`
492	`call input_open (fid,filename)`	492	`call input_open (fid,filename)`
493		493
494	`call input_var_wrf(fid,varname,f3t1,nx,ny,nz) ! field to trace`	494	`call input_var_wrf(fid,varname,f3t1,nx,ny,nz) ! field to trace`
495		495
496	`varname='geopot' ! geopot.height`	496	`varname='geopot' ! geopot.height`
497	`call input_var_wrf(fid,varname,p3t1,nx,ny,nz)`	497	`call input_var_wrf(fid,varname,p3t1,nx,ny,nz)`
498		498
499	`varname='geopots' ! geopot.height`	499	`varname='geopots' ! geopot.height`
500	`call input_var_wrf(fid,varname,spt1,nx,ny,1)`	500	`call input_var_wrf(fid,varname,spt1,nx,ny,1)`

Subversion Repositories lagranto.wrf

(root)/branches/met.nps.monterey/trace/trace.f90 – Rev 2 → 10