WebSVN – lagranto.um – Diff – /tags/2.2/trace/trace.f


      PROGRAM trace
 
C     ********************************************************************
C     *                                                                  *
C     * Trace fields along trajectories                                  *
C     *                                                                  *
C     * Heini Wernli       first version:       April 1993               *
C     * Michael Sprenger   major upgrade:       2008-2009                *
C     *                                                                  *
C     ********************************************************************
 
      implicit none
      
c     --------------------------------------------------------------------
c     Declaration of parameters
c     --------------------------------------------------------------------
 
c     Maximum number of levels for input files
      integer   nlevmax
      parameter (nlevmax=100)
 
c     Maximum number of input files (dates, length of trajectories)
      integer   ndatmax
      parameter (ndatmax=500)
 
c     Numerical epsilon (for float comparison)
      real      eps
      parameter (eps=0.001)
 
c     Conversion factors
      real      pi180                                   ! deg -> rad
      parameter (pi180=3.14159/180.)
      real      deg2km                                  ! deg -> km (at equator)
      parameter (deg2km=111.)
 
c     Prefix for primary and secondary fields
      character charp
      character chars
      parameter (charp='P')
      parameter (chars='S')
 
c     --------------------------------------------------------------------
c     Declaration of variables
c     --------------------------------------------------------------------
 
c     Input and output format for trajectories (see iotra.f)
      integer   inpmode
      integer   outmode
 
c     Input parameters
      character*80                           inpfile         ! Input trajectory file
      character*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*80                           tvar0(200)      ! Tracing variable (with mode specification)
      character*80                           tvar(200)       ! Tracing variable name (only the variable)
      character*1                            tfil(200)       ! Filename prefix 
      real                                   fac(200)        ! Scaling factor 
      real                                   shift_val(200)  ! Shift in space and time relative to trajectory position
      character*80                           shift_dir(200)  ! Direction of shift
      integer                                compfl(200)     ! Computation flag (1=compute)
      integer                                numdat          ! Number of input files
      character*13                           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*20                           startdate       ! First time/date on trajectory
      character*20                           enddate         ! Last time/date on trajectory
      integer                                ntrace1         ! Count trace and additional variables
      character*80                           timecheck       ! Either 'yes' or 'no'
 
c     Trajectories
      real,allocatable, dimension (:,:,:) :: trainp          ! Input trajectories (ntra,ntim,ncol)
      real,allocatable, dimension (:,:,:) :: traint          ! Internal trajectories (ntra,ntim,ncol+ntrace1)
      real,allocatable, dimension (:,:,:) :: traout          ! Output trajectories (ntra,ntim,ncol+ntrace0)
      integer                                reftime(6)      ! Reference date
      character*80                           varsinp(100)    ! Field names for input trajectory
      character*80                           varsint(100)    ! Field names for internal trajectory
      character*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
 
c     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*80                           svars(100)      ! List of variables on S file
      character*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
      
c     Grid description
      real                                   pollon,pollat   ! Longitude/latitude of pole
      real                                   ak(100)         ! Vertical layers and levels
      real                                   bk(100) 
      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
 
c     Auxiliary variables
      integer                                i,j,k,l,n
      real                                   rd
      character*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*80                           string
      real                                   lon,lat,rlon,rlat
 
c     Externals 
      real                                   int_index4
      external                               int_index4
 
      real                                   lmtolms 
      real                                   phtophs    
      real                                   lmstolm
      real                                   phstoph        
      external                               lmtolms,phtophs
      external                               lmstolm,phstoph
 
c     --------------------------------------------------------------------
c     Start of program, Read parameters, get grid parameters
c     --------------------------------------------------------------------
 
c     Write start message
      print*,'========================================================='
      print*,'              *** START OF PROGRAM TRACE ***'
      print*
 
c     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,'(a13)') dat(i)
          enddo
       else
          do i=numdat,1,-1
             read(10,'(a13)') 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
      close(10)
 
c     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
      
c     Save the tracing variabel (including all mode specifications)
      do i=1,ntrace0
         tvar0(i) = tvar(i)
      enddo
 
c     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
 
C     Convert time shifts <tst,ten> from <hh.mm> into fractional time
      call hhmm2frac(tst,frac)
      tst = frac
      call hhmm2frac(ten,frac)
      ten = frac
 
c     Set the time for the first data file (depending on forward/backward mode)
      if (fbflag.eq.1) then
        tstart = -tst
      else
        tstart = tst
      endif
 
c     Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)
c     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 (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >                 tstart,pollon,pollat,rd,rd,nz,rd,rd,rd,timecheck)
      call input_close(fid)
 
C     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 ***'
 
C     Get memory for trajectory arrays
      allocate(trainp(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra      ***' 
 
c     Set the flags for periodic domains
      if ( abs(xmax-xmin-360.).lt.eps ) then
         per = 1
      elseif ( abs(xmax-xmin-360.+dx).lt.eps ) then
         per = 2
      else
         per = 0
      endif
 
C     Set logical flag for periodic data set (hemispheric or not)
      if (per.eq.0.) then
         delta=xmax-xmin-360.
         if (abs(delta+dx).lt.eps) then               ! Program aborts: arrays must be closed
            print*,' ERROR: arrays must be closed... Stop'
        else if (abs(delta).lt.eps) then              ! Periodic and hemispheric
           hem=1
           per=360.
        endif
      else                                            ! Periodic and hemispheric
         hem=1
      endif
 
c     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)
      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*
 
c     --------------------------------------------------------------------
c     Load the input trajectories
c     --------------------------------------------------------------------
 
c     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)
 
c     Coordinate rotation
      do i=1,ntra
         do j=1,ntim
 
            lon = trainp(i,j,2)
            lat = trainp(i,j,3)
            
            if ( abs(pollat-90.).gt.eps) then
               rlon = lmtolms(lat,lon,pollat,pollon)
               rlat = phtophs(lat,lon,pollat,pollon)
            else
               rlon = lon
               rlat = lat
            endif
            
            trainp(i,j,2) = rlon
            trainp(i,j,3) = rlat
            
         enddo
      enddo
 
 
c     Check that first four columns correspond to time,lon,lat,p
      if ( (varsinp(1).ne.'time' ).or.
     >     (varsinp(2).ne.'xpos' ).and.(varsinp(2).ne.'lon' ).or.
     >     (varsinp(3).ne.'ypos' ).and.(varsinp(3).ne.'lat' ).or.
     >     (varsinp(4).ne.'ppos' ).and.(varsinp(4).ne.'p'   ) )
     >then
         print*,' ERROR: problem with input trajectories ...'
         stop
      endif
      varsinp(1) = 'time'
      varsinp(2) = 'lon'
      varsinp(3) = 'lat'
      varsinp(4) = 'p'
 
c     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*
 
c     --------------------------------------------------------------------
c     Check dependencies for trace fields which must be calculated
c     --------------------------------------------------------------------
 
c     Set the counter for extra fields
      ntrace1 = ntrace0
 
c     Loop over all tracing variables
      i = 1
      do while (i.le.ntrace1)
 
c        Skip fields which must be available on the input files
         if (i.le.ntrace0) then
            if (compfl(i).eq.0) goto 100
         endif
 
c        Get the dependencies for potential temperature (TH)
         if ( tvar(i).eq.'TH' ) then
            varname='P'                                   ! P
            call add2list(varname,tvar,ntrace1)
            varname='T'                                   ! T
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for potential temperature (TH)
         elseif ( tvar(i).eq.'RHO' ) then
            varname='P'                                   ! P
            call add2list(varname,tvar,ntrace1)
            varname='T'                                   ! T
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for relative humidity (RH)
         elseif ( tvar(i).eq.'RH' ) then
            varname='P'                                   ! P
            call add2list(varname,tvar,ntrace1)
            varname='T'                                   ! T
            call add2list(varname,tvar,ntrace1)
            varname='Q'                                   ! Q
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for equivalent potential temperature (THE)
         elseif ( tvar(i).eq.'THE' ) then
            varname='P'                                   ! P
            call add2list(varname,tvar,ntrace1)
            varname='T'                                   ! T
            call add2list(varname,tvar,ntrace1)
            varname='Q'                                   ! Q
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for latent heating rate (LHR)
         elseif ( tvar(i).eq.'LHR' ) then
            varname='P'                                   ! P
            call add2list(varname,tvar,ntrace1)
            varname='T'                                   ! T
            call add2list(varname,tvar,ntrace1)
            varname='Q'                                   ! Q
            call add2list(varname,tvar,ntrace1)
            varname='OMEGA'                               ! OMEGA
            call add2list(varname,tvar,ntrace1)
            varname='RH'                                  ! RH
            call add2list(varname,tvar,ntrace1)
         
c        Get the dependencies for wind speed (VEL)
         elseif ( tvar(i).eq.'VEL' ) then
            varname='U'                                   ! U
            call add2list(varname,tvar,ntrace1)
            varname='V'                                   ! V
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for wind direction (DIR)
         elseif ( tvar(i).eq.'DIR' ) then
            varname='U'                                   ! U
            call add2list(varname,tvar,ntrace1)
            varname='V'                                   ! V
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for du/dx (DUDX)
         elseif ( tvar(i).eq.'DUDX' ) then
            varname='U:+1DLON'                            ! U:+1DLON
            call add2list(varname,tvar,ntrace1)
            varname='U:-1DLON'                            ! U:-1DLON
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for dv(dx (DVDX)
         elseif ( tvar(i).eq.'DVDX' ) then
            varname='V:+1DLON'                            ! V:+1DLON
            call add2list(varname,tvar,ntrace1)
            varname='V:-1DLON'                            ! V:-1DLON
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for du/dy (DUDY)
         elseif ( tvar(i).eq.'DUDY' ) then
            varname='U:+1DLAT'                            ! U:+1DLAT
            call add2list(varname,tvar,ntrace1)
            varname='U:-1DLAT'                            ! U:-1DLAT
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for dv/dy (DVDY)
         elseif ( tvar(i).eq.'DVDY' ) then
            varname='V:+1DLAT'                            ! V:+1DLAT
            call add2list(varname,tvar,ntrace1)
            varname='V:-1DLAT'                            ! V:-1DLAT
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for du/dp (DUDP)
         elseif ( tvar(i).eq.'DUDP' ) then
            varname='U:+1DP'                              ! U:+1DP
            call add2list(varname,tvar,ntrace1)
            varname='U:-1DP'                              ! U:-1DP
            call add2list(varname,tvar,ntrace1)
            varname='P:+1DP'                              ! P:+1DP
            call add2list(varname,tvar,ntrace1)
            varname='P:-1DP'                              ! P:-1DP
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for dv/dp (DVDP)
         elseif ( tvar(i).eq.'DVDP' ) then
            varname='V:+1DP'                              ! V:+1DP
            call add2list(varname,tvar,ntrace1)
            varname='V:-1DP'                              ! V:-1DP
            call add2list(varname,tvar,ntrace1)
            varname='P:+1DP'                              ! P:+1DP
            call add2list(varname,tvar,ntrace1)
            varname='P:-1DP'                              ! P:-1DP
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for dt/dx (DTDX)
         elseif ( tvar(i).eq.'DTDX' ) then
            varname='T:+1DLON'                            ! T:+1DLON
            call add2list(varname,tvar,ntrace1)
            varname='T:-1DLON'                            ! T:-1DLON
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for dth/dy (DTHDY)
         elseif ( tvar(i).eq.'DTHDY' ) then
            varname='T:+1DLAT'                            ! T:+1DLON
            call add2list(varname,tvar,ntrace1)
            varname='T:-1DLAT'                            ! T:-1DLON
            call add2list(varname,tvar,ntrace1)
            varname='P'                                   ! P
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for dth/dx (DTHDX)
         elseif ( tvar(i).eq.'DTHDX' ) then
            varname='T:+1DLON'                            ! T:+1DLON
            call add2list(varname,tvar,ntrace1)
            varname='T:-1DLON'                            ! T:-1DLON
            call add2list(varname,tvar,ntrace1)
            varname='P'                                   ! P
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for dt/dy (DTDY)
         elseif ( tvar(i).eq.'DTDY' ) then
            varname='T:+1DLAT'                            ! T:+1DLON
            call add2list(varname,tvar,ntrace1)
            varname='T:-1DLAT'                            ! T:-1DLON
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for dt/dp (DTDP)
         elseif ( tvar(i).eq.'DTDP' ) then
            varname='T:+1DP'                              ! T:+1DP
            call add2list(varname,tvar,ntrace1)
            varname='T:-1DP'                              ! T:-1DP
            call add2list(varname,tvar,ntrace1)
            varname='P:+1DP'                              ! P:+1DP
            call add2list(varname,tvar,ntrace1)
            varname='P:-1DP'                              ! P:-1DP
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for dth/dp (DTHDP)
         elseif ( tvar(i).eq.'DTHDP' ) then
            varname='T:+1DP'                              ! T:+1DP
            call add2list(varname,tvar,ntrace1)
            varname='T:-1DP'                              ! T:-1DP
            call add2list(varname,tvar,ntrace1)
            varname='T'                                   ! T
            call add2list(varname,tvar,ntrace1)
            varname='P:+1DP'                              ! P:+1DP
            call add2list(varname,tvar,ntrace1)
            varname='P:-1DP'                              ! P:-1DP
            call add2list(varname,tvar,ntrace1)
            varname='P'                                   ! P
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for squared Brunt Vaiäläa frequency (NSQ)
         elseif ( tvar(i).eq.'NSQ' ) then
            varname='DTHDP'                                ! DTHDP
            call add2list(varname,tvar,ntrace1)
            varname='TH'                                   ! TH
            call add2list(varname,tvar,ntrace1)
            varname='RHO'                                  ! RHO
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for relative vorticity (RELVORT)
         elseif ( tvar(i).eq.'RELVORT' ) then
            varname='U'                                    ! U
            call add2list(varname,tvar,ntrace1)
            varname='DUDY'                                 ! DUDY
            call add2list(varname,tvar,ntrace1)
            varname='DVDX'                                 ! DVDX
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for relative vorticity (ABSVORT)
         elseif ( tvar(i).eq.'ABSVORT' ) then
            varname='U'                                    ! U
            call add2list(varname,tvar,ntrace1)
            varname='DUDY'                                 ! DUDY
            call add2list(varname,tvar,ntrace1)
            varname='DVDX'                                 ! DVDX
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for divergence (DIV)
         elseif ( tvar(i).eq.'DIV' ) then
            varname='V'                                    ! U
            call add2list(varname,tvar,ntrace1)
            varname='DUDX'                                 ! DUDX
            call add2list(varname,tvar,ntrace1)
            varname='DVDY'                                 ! DVDY
            call add2list(varname,tvar,ntrace1)
 
c        Get the dependencies for deformation (DEF)
         elseif ( tvar(i).eq.'DEF' ) then
            call add2list(varname,tvar,ntrace1)
            varname='DUDX'                                 ! DUDX
            call add2list(varname,tvar,ntrace1)
            varname='DVDY'                                 ! DVDY
            call add2list(varname,tvar,ntrace1)
            varname='DUDY'                                 ! DUDY
            call add2list(varname,tvar,ntrace1)
            varname='DVDX'                                 ! DVDX
            call add2list(varname,tvar,ntrace1)
       
c        Get the dependencies for potential vorticity (PV)
         elseif ( tvar(i).eq.'PV' ) then
            varname='ABSVORT'                              ! ABSVORT
            call add2list(varname,tvar,ntrace1)
            varname='DTHDP'                                ! DTHDP
            call add2list(varname,tvar,ntrace1)
            varname='DUDP'                                 ! DUDP
            call add2list(varname,tvar,ntrace1)
            varname='DVDP'                                 ! DVDP
            call add2list(varname,tvar,ntrace1)
            varname='DTHDX'                                ! DTHDX
            call add2list(varname,tvar,ntrace1)
            varname='DTHDY'                                ! DTHDY
            call add2list(varname,tvar,ntrace1)
       
c        Get the dependencies for Richardson number (RI)
         elseif ( tvar(i).eq.'RI' ) then
            varname='DUDP'                                 ! DUDP
            call add2list(varname,tvar,ntrace1)
            varname='DVDP'                                 ! DVDP
            call add2list(varname,tvar,ntrace1)
            varname='NSQ'                                  ! NSQ
            call add2list(varname,tvar,ntrace1)
            varname='RHO'                                  ! RHO
            call add2list(varname,tvar,ntrace1)
            
c        Get the dependencies for Ellrod&Knapp's turbulence index (TI)
         elseif ( tvar(i).eq.'TI' ) then
            varname='DEF'                                  ! DEF
            call add2list(varname,tvar,ntrace1)       
            varname='DUDP'                                 ! DUDP
            call add2list(varname,tvar,ntrace1)       
            varname='DVDP'                                 ! DVDP
            call add2list(varname,tvar,ntrace1)       
            varname='RHO'                                  ! RHO
            call add2list(varname,tvar,ntrace1)       
 
         endif
 
c        Exit point for handling additional fields
 100     continue
         i = i + 1
 
      enddo
 
c     Save the full variable name (including shift specification)
      do i=1,ncol
         varsint(i)      = varsinp(i)
      enddo
      do i=1,ntrace1
         varsint(i+ncol) = tvar(i)
      enddo
 
c     Split the variable name and set flags
      do i=1,ntrace1
 
c        Set the scaling factor
         if ( i.gt.ntrace0 ) fac(i) = 1.
 
c        Set the base variable name, the shift and the direction
         call splitvar(tvar(i),shift_val(i),shift_dir(i) )
 
c        Set the prefix of the file name
         if ( ( compfl(i).eq.0 ).or.(i.gt.ntrace0) ) 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
         else
            tfil(i) = '*'
         endif
         
c        Set the computational flag
         if ( (tvar(i).eq.'P'   ).or.
     >        (tvar(i).eq.'PLAY').or.
     >        (tvar(i).eq.'PLEV') )
     >   then    
            compfl(i) = 0
            tfil(i)   = charp
         elseif ( ( tfil(i).eq.charp ).or.( tfil(i).eq.chars ) ) then
            compfl(i) = 0
         else
            compfl(i) = 1
         endif
 
      enddo
 
c     Check whether the shift modes are supported
      do i=1,ntrace1
         if ( ( shift_dir(i).ne.'nil'     ).and.
     >        ( shift_dir(i).ne.'DLON'    ).and.
     >        ( shift_dir(i).ne.'DLAT'    ).and.
     >        ( shift_dir(i).ne.'DP'      ).and.
     >        ( shift_dir(i).ne.'HPA'      ).and.
     >        ( shift_dir(i).ne.'KM(LON)' ).and.
     >        ( shift_dir(i).ne.'KM(LAT)' ).and.
     >        ( shift_dir(i).ne.'H'       ).and.
     >        ( shift_dir(i).ne.'MIN'     ).and.
     >        ( shift_dir(i).ne.'INDP'    ) )
     >   then
            print*,' ERROR: shift mode ',trim(shift_dir(i)),
     >             ' not supported'
            stop
         endif
      enddo
 
c     Write status information
      print*
      print*,'---- COMPLETE TABLE FOR TRACING -------------------------'   
      print*
      do i=1,ntrace1
         if ( ( shift_dir(i).ne.'nil' ) ) then
            write(*,'(i4,a4,a8,f10.2,a8,3x,a1,i5)') 
     >           i,' : ',trim(tvar(i)),shift_val(i),trim(shift_dir(i)),
     >           tfil(i),compfl(i)
         else
            write(*,'(i4,a4,a8,10x,8x,3x,a1,i5)') 
     >           i,' : ',trim(tvar(i)),tfil(i),compfl(i)
         endif
      enddo
      if ( i.eq.ntrace0 ) print*
      
c     --------------------------------------------------------------------
c     Prepare the internal and output trajectories
c     --------------------------------------------------------------------
 
c     Allocate memory for internal trajectories
      allocate(traint(ntra,ntim,ncol+ntrace1),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array traint   ***'
 
c     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
 
c     Set the flags for ready fields/colums - at begin only read-in fields are ready
      do i=1,ncol
         fok(i) = 1
      enddo
      do i=ncol+1,ntrace1
         fok(i) = 0
      enddo
 
c     --------------------------------------------------------------------
c     Trace the fields (fields available on input files)
c     --------------------------------------------------------------------
 
      print*
      print*,'---- TRACING FROM PRIMARY AND SECONDARY DATA FILES ------'     
 
c     Loop over all tracing fields
      do i=1,ntrace1
         
C         Skip fields which must be computed (compfl=1), will be handled later
          if (compfl(i).ne.0)  goto 110
 
c         Write some status information
          print*
          print*,' Now tracing             : ',
     >         trim(tvar(i)),shift_val(i),trim(shift_dir(i)),
     >         compfl(i),' ',trim(tfil(i))
 
c         Set the flag for ready field/column
          fok(ncol+i) = 1
 
c         Reset flags for load manager
          iloaded0 = -1
          iloaded1 = -1
 
c         Reset the counter for fields outside domain
          noutside = 0
 
c         Loop over all times
          do j=1,ntim
 
c            Convert trajectory time from hh.mm to fractional time
             call hhmm2frac(trainp(1,j,1),tfrac)
 
c            Shift time if requested
             if ( shift_dir(i).eq.'H' ) then
                tfrac = tfrac + shift_val(i)
             elseif ( shift_dir(i).eq.'MIN' ) then
                tfrac = tfrac + shift_val(i)/60.
             endif
 
c            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
 
c            Load manager: Check whether itime0 can be copied from itime1
             if ( itime0.eq.iloaded1 ) then
                f3t0     = f3t1
                p3t0     = p3t1
                spt0     = spt1
                iloaded0 = itime0
             endif
 
c            Load manager: Check whether itime1 can be copied from itime0
             if ( itime1.eq.iloaded0 ) then
                f3t1     = f3t0
                p3t1     = p3t0
                spt1     = spt0
                iloaded1 = itime1
             endif
 
c            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_wind 
     >               (fid,varname,f3t0,tload,stagz,mdv,
     >               xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)     
                call input_grid      
     >               (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >                tload,pollon,pollat,p3t0,spt0,nz,ak,bk,stagz,
     >                timecheck)
                call input_close(fid)
                
                iloaded0 = itime0
                
             endif
 
c            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_wind 
     >               (fid,varname,f3t1,tload,stagz,mdv,
     >               xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)      
                call input_grid      
     >               (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >               tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,
     >               timecheck)
                call input_close(fid)
                
                iloaded1 = itime1
                
             endif
 
c            Loop over all trajectories
             do k=1,ntra
                
c               Set the horizontal position where to interpolate to
                x0       = traint(k,j,2)                          ! Longitude
                y0       = traint(k,j,3)                          ! Latitude
 
c               Set the vertical position where to interpolate to
                if ( nz.gt.1 ) then
                   p0 = traint(k,j,4)                             ! Pressure (3D tracing)
                else
                   p0 = 1050.                                     ! Lowest level (2D tracing)
                endif
 
c               Set negative pressures to mdv
                if (p0.lt.0.) traint(k,j,4) = mdv
 
c               Set the relative time
                call hhmm2frac(traint(k,j,1),tfrac)
                reltpos0 = fbflag * (tfrac-time0)/timeinc    
 
c               Make adjustments depending on the shift flag
                if ( shift_dir(i).eq.'DLON' ) then                         ! DLON
                   x0 = x0 + shift_val(i)
 
                elseif  ( shift_dir(i).eq.'DLAT' ) then                    ! DLAT
                   y0 = y0 + shift_val(i)
 
                elseif ( shift_dir(i).eq.'KM(LON)' ) then                  ! KM(LON)
                   x0 = x0 + shift_val(i)/deg2km * 1./cos(y0*pi180)
 
                elseif ( shift_dir(i).eq.'KM(LAT)' ) then                  ! KM(LAT)
                   y0 = y0 + shift_val(i)/deg2km 
 
                elseif ( shift_dir(i).eq.'HPA' ) then                      ! HPA
                   p0 = p0 + shift_val(i)
 
                elseif ( shift_dir(i).eq.'DP' ) then                       ! DP
                   call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0,
     >                              p3t0,p3t1,spt0,spt1,3,
     >                              nx,ny,nz,xmin,ymin,dx,dy,mdv)
                   pind = pind - shift_val(i)
                   p0   = int_index4(p3t0,p3t1,nx,ny,nz,
     >                              xind,yind,pind,reltpos0,mdv)
 
                elseif ( shift_dir(i).eq.'INDP' ) then         
                   p0   = int_index4(p3t0,p3t1,nx,ny,nz,
     >                            xind,yind,shift_val(i),reltpos0,mdv)
 
                endif
                   
c               Handle periodic boundaries in zonal direction
                if ( (x0.gt.xmax).and.(per.ne.0) ) x0 = x0 - 360.
                if ( (x0.lt.xmin).and.(per.ne.0) ) x0 = x0 + 360.
 
c               Handle pole problems for hemispheric data (taken from caltra.f)
                if ((hem.eq.1).and.(y0.gt.90.)) then
                   y0=180.-y0
                   x0=x0+per/2.
                endif
                if ((hem.eq.1).and.(y0.lt.-90.)) then
                   y0=-180.-y0
                   x0=x0+per/2.
                endif
                if (y0.gt.89.99) then
                   y0=89.99
                endif
 
C               Get the index where to interpolate (x0,y0,p0)
                if ( (abs(x0-mdv).gt.eps).and.
     >               (abs(x0-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
 
c               Do the interpolation
                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           
                   f0 = int_index4(f3t0,f3t1,nx,ny,nz,
     >                             xind,yind,pind,reltpos0,mdv)
                elseif (noutside.lt.10) then
                   print*,' ',trim(tvar(i)),' @ ',x0,y0,p0,'outside'
                   f0       = mdv
                   noutside = noutside + 1
                elseif (noutside.eq.10) then
                   print*,' ...more than 10 outside...'
                   f0       = mdv
                   noutside = noutside + 1
                else
                   f0       = mdv
                endif
 
c               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
 
          enddo
 
c         Exit point for loop over all tracing variables
 110      continue
 
       enddo
 
c     --------------------------------------------------------------------
c     Calculate additional fields along the trajectories
c     --------------------------------------------------------------------
 
      print*
      print*,'---- CALCULATE ADDITIONAL FIELDS FROM TRAJECTORY TABLE --'
 
c     Loop over all tracing fields
      do i=ntrace1,1,-1
 
C         Skip fields which must not be computed (compfl=0)
          if (compfl(i).eq.0) goto 120
 
c         Write some status information
          print*
          write(*,'(a10,f10.2,a5,i3,3x,a2)')
     >         trim(tvar(i)),shift_val(i),trim(shift_dir(i)),
     >         compfl(i),trim(tfil(i))
 
c         Loop over trajectories and times
          do j=1,ntra
          do k=1,ntim
                
c            Potential temperature (TH)
             if  ( varsint(i+ncol).eq.'TH' ) then
                
                varname='T'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='p'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
 
                call calc_TH  (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                             traint(j,k,ind2) )
 
c            Density (RHO)
             elseif  ( varsint(i+ncol).eq.'RHO' ) then
                
                varname='T'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='p'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
 
                call calc_RHO (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                             traint(j,k,ind2) )
                
c            Relative humidity (RH)
             elseif  ( varsint(i+ncol).eq.'RH' ) then
                
                varname='T'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='p'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='Q'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)               
 
                call calc_RH (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                            traint(j,k,ind2),
     >                                            traint(j,k,ind3) )
 
c            Equivalent potential temperature (THE)
             elseif  ( varsint(i+ncol).eq.'THE' ) then
                
                varname='T'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='p'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='Q'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)                
 
                call calc_THE (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                             traint(j,k,ind2),
     >                                             traint(j,k,ind3) )
 
c            Latent heating rate (LHR)
             elseif  ( varsint(i+ncol).eq.'LHR' ) then
                
                varname='T'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='p'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='Q'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)                
                varname='OMEGA'
                call list2ind (ind4,varname,varsint,fok,ncol+ntrace1)                
                varname='RH'
                call list2ind (ind5,varname,varsint,fok,ncol+ntrace1)                
 
                call calc_LHR (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                             traint(j,k,ind2),
     >                                             traint(j,k,ind3),
     >                                             traint(j,k,ind4),
     >                                             traint(j,k,ind5) )
 
c            Wind speed (VEL)
             elseif  ( varsint(i+ncol).eq.'VEL' ) then
                
                varname='U'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='V'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
 
                call calc_VEL (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                             traint(j,k,ind2) )
 
c            Wind direction (DIR)
             elseif  ( varsint(i+ncol).eq.'DIR' ) then
                
                varname='U'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='V'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DIR (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                             traint(j,k,ind2) )
 
c            Zonal gradient of U (DUDX)
             elseif  ( varsint(i+ncol).eq.'DUDX' ) then
                
                varname='U:+1DLON'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='U:-1DLON'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='lat'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DUDX (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                              traint(j,k,ind2),
     >                                              traint(j,k,ind3) )
                
c            Zonal gradient of V (DVDX)
             elseif  ( varsint(i+ncol).eq.'DVDX' ) then
                
                varname='V:+1DLON'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='V:-1DLON'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='lat'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DVDX (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                              traint(j,k,ind2),
     >                                              traint(j,k,ind3) )
                
c            Zonal gradient of T (DTDX)
             elseif  ( varsint(i+ncol).eq.'DVDX' ) then
                
                varname='T:+1DLON'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='T:-1DLON'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='lat'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DTDX (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                              traint(j,k,ind2),
     >                                              traint(j,k,ind3) )
 
c            Zonal gradient of TH (DTHDX)
             elseif  ( varsint(i+ncol).eq.'DTHDX' ) then
                
                varname='T:+1DLON'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='T:-1DLON'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='P'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)
                varname='lat'
                call list2ind (ind4,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DTHDX (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                               traint(j,k,ind2),
     >                                               traint(j,k,ind3),
     >                                               traint(j,k,ind4) )
 
c            Meridional gradient of U (DUDY)
             elseif  ( varsint(i+ncol).eq.'DUDY' ) then
                
                varname='U:+1DLAT'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='U:-1DLAT'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DUDY (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                              traint(j,k,ind2) )
 
c            Meridional gradient of V (DVDY)
             elseif  ( varsint(i+ncol).eq.'DVDY' ) then
                
                varname='V:+1DLAT'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='V:-1DLAT'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DVDY (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                              traint(j,k,ind2) )
 
c            Meridional gradient of T (DTDY)
             elseif  ( varsint(i+ncol).eq.'DTDY' ) then
                
                varname='T:+1DLAT'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='T:-1DLAT'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DTDY (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                              traint(j,k,ind2) )
 
c            Meridional gradient of TH (DTHDY)
             elseif  ( varsint(i+ncol).eq.'DTHDY' ) then
                
                varname='T:+1DLAT'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='T:-1DLAT'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='P'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DTDY (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                              traint(j,k,ind2),
     >                                              traint(j,k,ind3) )
 
 
c            Vertical wind shear DU/DP (DUDP)
             elseif  ( varsint(i+ncol).eq.'DUDP' ) then
                
                varname='U:+1DP'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='U:-1DP'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='P:+1DP'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1) 
                varname='P:-1DP'
                call list2ind (ind4,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DUDP (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                              traint(j,k,ind2),
     >                                              traint(j,k,ind3),
     >                                              traint(j,k,ind4) )
 
c            Vertical wind shear DV/DP (DVDP)
             elseif  ( varsint(i+ncol).eq.'DVDP' ) then
                
                varname='V:+1DP'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='V:-1DP'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='P:+1DP'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1) 
                varname='P:-1DP'
                call list2ind (ind4,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DVDP (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                              traint(j,k,ind2),
     >                                              traint(j,k,ind3),
     >                                              traint(j,k,ind4) )
 
c            Vertical derivative of T (DTDP)
             elseif  ( varsint(i+ncol).eq.'DTDP' ) then
                
                varname='T:+1DP'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='T:-1DP'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='P:+1DP'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1) 
                varname='P:-1DP'
                call list2ind (ind4,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DTDP (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                              traint(j,k,ind2),
     >                                              traint(j,k,ind3),
     >                                              traint(j,k,ind4) )
 
c            Vertical derivative of TH (DTHDP)
             elseif  ( varsint(i+ncol).eq.'DTHDP' ) then
                
                varname='T:+1DP'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='T:-1DP'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='P:+1DP'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1) 
                varname='P:-1DP'
                call list2ind (ind4,varname,varsint,fok,ncol+ntrace1)
                varname='P'
                call list2ind (ind5,varname,varsint,fok,ncol+ntrace1)
                varname='T'
                call list2ind (ind6,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DTHDP (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                               traint(j,k,ind2),
     >                                               traint(j,k,ind3),
     >                                               traint(j,k,ind4),
     >                                               traint(j,k,ind5),
     >                                               traint(j,k,ind6) )
 
c            Squared Brunt-Vaisäla frequency (NSQ)
             elseif  ( varsint(i+ncol).eq.'NSQ' ) then
                
                varname='DTHDP'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='TH'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='RHO'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)
 
                call calc_NSQ (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                             traint(j,k,ind2),
     >                                             traint(j,k,ind3))
 
c            Relative vorticity (RELVORT)
             elseif  ( varsint(i+ncol).eq.'RELVORT' ) then
 
                varname='DUDY'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='DVDX'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='U'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)
                varname='lat'
                call list2ind (ind4,varname,varsint,fok,ncol+ntrace1)
 
                call calc_RELVORT (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                                 traint(j,k,ind2),
     >                                                 traint(j,k,ind3),
     >                                                 traint(j,k,ind4))
 
c            Absolute vorticity (ABSVORT)
             elseif  ( varsint(i+ncol).eq.'ABSVORT' ) then
 
                varname='DUDY'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='DVDX'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='U'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)
                varname='lat'
                call list2ind (ind4,varname,varsint,fok,ncol+ntrace1)
                varname='lon'
                call list2ind (ind5,varname,varsint,fok,ncol+ntrace1)
 
                call calc_ABSVORT (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                                 traint(j,k,ind2),
     >                                                 traint(j,k,ind3),
     >                                                 traint(j,k,ind4),
     >                                                 traint(j,k,ind5),
     >                                                 pollon,pollat  )
 
c            Divergence (DIV)
             elseif  ( varsint(i+ncol).eq.'DIV' ) then
 
                varname='DUDX'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='DVDY'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='V'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)
                varname='lat'
                call list2ind (ind4,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DIV (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                             traint(j,k,ind2),
     >                                             traint(j,k,ind3),
     >                                             traint(j,k,ind4))
 
c            Deformation (DEF)
             elseif  ( varsint(i+ncol).eq.'DEF' ) then
 
                varname='DUDX'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='DVDX'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='DUDY'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)
                varname='DVDY'
                call list2ind (ind4,varname,varsint,fok,ncol+ntrace1)
 
                call calc_DEF (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                             traint(j,k,ind2),
     >                                             traint(j,k,ind3),
     >                                             traint(j,k,ind4))
 
c            Potential Vorticity (PV)
             elseif  ( varsint(i+ncol).eq.'PV' ) then
 
                varname='ABSVORT'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='DTHDP'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='DUDP'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)
                varname='DVDP'
                call list2ind (ind4,varname,varsint,fok,ncol+ntrace1)
                varname='DTHDX'
                call list2ind (ind5,varname,varsint,fok,ncol+ntrace1)
                varname='DTHDY'
                call list2ind (ind6,varname,varsint,fok,ncol+ntrace1)
 
                call calc_PV (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                            traint(j,k,ind2),
     >                                            traint(j,k,ind3),
     >                                            traint(j,k,ind4),
     >                                            traint(j,k,ind5),
     >                                            traint(j,k,ind6) )
 
c            Richardson number (RI)
             elseif  ( varsint(i+ncol).eq.'RI' ) then
 
                varname='DUDP'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='DVDP'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='NSQ'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)
                varname='RHO'
                call list2ind (ind4,varname,varsint,fok,ncol+ntrace1)
 
                call calc_RI (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                            traint(j,k,ind2),
     >                                            traint(j,k,ind3),
     >                                            traint(j,k,ind4) )
 
c            Ellrod and Knapp's turbulence idicator (TI)
             elseif  ( varsint(i+ncol).eq.'TI' ) then
 
                varname='DEF'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='DUDP'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                varname='DVDP'
                call list2ind (ind3,varname,varsint,fok,ncol+ntrace1)
                varname='RHO'
                call list2ind (ind4,varname,varsint,fok,ncol+ntrace1)
 
                call calc_TI (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                            traint(j,k,ind2),
     >                                            traint(j,k,ind3),
     >                                            traint(j,k,ind4) )
 
c            Spherical distance from starting position (DIST0)
             elseif  ( varsint(i+ncol).eq.'DIST0' ) then
 
                varname='lon'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='lat'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                
                call calc_DIST0 (traint(j,k,ncol+i), traint(j,k,ind1), 
     >                                               traint(j,k,ind2),
     >                                               traint(j,1,ind1), 
     >                                               traint(j,1,ind2) )
 
c            Spherical distance length of trajectory (DIST)
             elseif  ( varsint(i+ncol).eq.'DIST' ) then
 
                varname='lon'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='lat'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                
                if ( k.eq.1 ) then
                   traint(j,k,ncol+i) = 0.
                else
                   call calc_DIST0 (delta, traint(j,k  ,ind1), 
     >                                     traint(j,k  ,ind2),
     >                                     traint(j,k-1,ind1), 
     >                                     traint(j,k-1,ind2) )
                   traint(j,k,ncol+i) = traint(j,k-1,ncol+i) + delta
                endif
 
c            Heading of the trajectory (HEAD)
             elseif  ( varsint(i+ncol).eq.'HEAD' ) then
 
                varname='lon'
                call list2ind (ind1,varname,varsint,fok,ncol+ntrace1) 
                varname='lat'
                call list2ind (ind2,varname,varsint,fok,ncol+ntrace1)
                 
                if (k.eq.ntim) then
                   traint(j,k,ncol+i) = mdv
                else
                   call calc_HEAD (traint(j,k,ncol+i), 
     >                                         traint(j,k  ,ind1), 
     >                                         traint(j,k  ,ind2),
     >                                         traint(j,k+1,ind1), 
     >                                         traint(j,k+1,ind2) )
                endif
 
                   
c            Invalid tracing variable
             else
 
                print*,' ERROR: invalid tracing variable ',
     >                   trim(varsint(i+ncol))
                stop
 
 
             endif
 
c         End loop over all trajectories and times   
          enddo
          enddo
 
c         Set the flag for a ready field/column
          fok(ncol+i) = 1
             
 
c         Exit point for loop over all tracing fields
 120      continue
          
       enddo
 
c     --------------------------------------------------------------------
c     Write output to output trajectory file
c     --------------------------------------------------------------------
 
c     Write status information
      print*
      print*,'---- WRITE OUTPUT TRAJECTORIES --------------------------'   
      print*
 
c     Allocate memory for internal trajectories
      allocate(traout(ntra,ntim,ncol+ntrace0),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array traout   ***'
 
c     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
 
c     Set the variable names for output trajectory
      do i=1,ncol+ntrace0
         varsout(i)      = varsint(i)
      enddo  
 
c     Coordinate rotation
      do i=1,ntra
         do j=1,ntim
 
            rlon = traout(i,j,2)
            rlat = traout(i,j,3)
            
            if ( abs(pollat-90.).gt.eps) then
               lon = lmstolm(rlat,rlon,pollat,pollon)
               lat = phstoph(rlat,rlon,pollat,pollon)
            else
               lon = rlon
               lat = rlat
            endif
            
            traout(i,j,2) = lon
            traout(i,j,3) = lat
            
         enddo
      enddo
 
 
c     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)
 
c     Write some status information, and end of program message
      print*  
      print*,'---- STATUS INFORMATION --------------------------------'
      print*
      print*,' ok'
      print*
      print*,'              *** END OF PROGRAM TRACE ***'
      print*,'========================================================='
 
 
      end 
 
 
 
c     ******************************************************************
c     * SUBROUTINE SECTION                                             *
c     ******************************************************************
 
c     ------------------------------------------------------------------
c     Add a variable to the list if not yet included in this list
c     ------------------------------------------------------------------
     
      subroutine add2list (varname,list,nlist)
 
      implicit none
      
c     Declaration of subroutine parameters
      character*80    varname
      character*80    list(200)
      integer         nlist
 
c     Auxiliray variables
      integer i,j
      integer isok
 
c     Expand the list, if necessary
      isok = 0
      do i=1,nlist      
         if ( list(i).eq.varname ) isok = 1
      enddo
      if ( isok.eq.0 ) then
         nlist       = nlist + 1
         list(nlist) = varname
      endif
 
c     Check for too large number of fields
      if ( nlist.ge.200) then
         print*,' ERROR: too many additional fields for tracing ...'
         stop
      endif
 
      end
 
 
c     ------------------------------------------------------------------
c     Get the index of a variable in the list
c     ------------------------------------------------------------------
 
      subroutine list2ind (ind,varname,list,fok,nlist)
 
      implicit none
      
c     Declaration of subroutine parameters
      integer         ind
      character*80    varname
      character*80    list(200)
      integer         fok(200)
      integer         nlist
 
c     Auxiliray variables
      integer i,j
      integer isok   
 
c     Get the index - error message if not found
      ind = 0
      do i=1,nlist
         if ( varname.eq.list(i) ) then
            ind = i
            goto 100
         endif
      enddo
      
      if ( ind.eq.0) then
         print*
         print*,' ERROR: cannot find ',trim(varname),' in list ...'
         do i=1,nlist
            print*,i,trim(list(i))
         enddo
         print*
         stop
      endif
 
c     Exit point
 100  continue
 
c     Check whether the field/column is ready
      if ( fok(ind).eq.0 ) then
         print*
         print*,' ERROR: unresolved dependence : ',trim(list(ind))
         print*
         stop
      endif
 
      end
 
 
c     ------------------------------------------------------------------
c     Split the variable name into name, shift and direction 
c     ------------------------------------------------------------------
 
      subroutine splitvar (tvar,shift_val,shift_dir)
         
      implicit none
      
c     Declaration of subroutine parameters
      character*80 tvar
      real         shift_val
      character*80 shift_dir
 
c     Auxiliary variables
      integer      i,j
      integer      icolon,inumber
      character*80 name
      character    ch
 
c     Save variable name
      name = tvar
 
c     Search for colon
      icolon=0
      do i=1,80
         if ( (name(i:i).eq.':').and.(icolon.ne.0) ) goto 100
         if ( (name(i:i).eq.':').and.(icolon.eq.0) ) icolon=i
      enddo
 
c     If there is a colon, split the variable name
      if ( icolon.ne.0 ) then
 
         tvar = name(1:(icolon-1))
         
         do i=icolon+1,80
            ch = name(i:i)
            if ( ( ch.ne.'0' ).and. ( ch.ne.'1' ).and.( ch.ne.'2' ).and.
     >           ( ch.ne.'3' ).and. ( ch.ne.'4' ).and.( ch.ne.'5' ).and.
     >           ( ch.ne.'6' ).and. ( ch.ne.'7' ).and.( ch.ne.'8' ).and.
     >           ( ch.ne.'9' ).and. ( ch.ne.'+' ).and.( ch.ne.'-' ).and.
     >           ( ch.ne.'.' ).and. ( ch.ne.' ' )  )
     >      then
               inumber = i
               exit
            endif
         enddo
 
         read(name( (icolon+1):(inumber-1) ),*) shift_val
 
         shift_dir = name(inumber:80)
 
      else
         
         shift_dir = 'nil'
         shift_val = 0.
 
      endif
            
      return
 
c     Error handling
 100  continue
      
      print*,' ERROR: cannot split variable name ',trim(tvar)
      stop
      
      end
      
      
         
         
         
      
 
 
      

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

Subversion Repositories lagranto.um

(root)/tags/2.2/trace/trace.f – Rev 5 → 14