Subversion Repositories lagranto.um

      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