WebSVN – lagranto.ecmwf – Diff – /trunk/trace/trace.f90


PROGRAM trace
 
  ! ********************************************************************
  ! *                                                                  *
  ! * Trace fields along trajectories                                  *
  ! *                                                                  *
  ! * April 1993: First version (Heini Wernli)                         *
  ! * 2008-2009 : Major upgrades (Michael Sprenger)                    *
  ! * Mar 2012: Clustering option (Bojan Skerlak)                      *  
  ! * Nov 2012: Circle options (")                                     *
  ! * Jul 2013: user-defined PV,TH @ clustering mode (")               *
  ! *                                                                  *
  ! ********************************************************************
 
  implicit none
 
  ! --------------------------------------------------------------------
  ! Declaration of parameters
  ! --------------------------------------------------------------------
 
  ! Maximum number of levels for input files
  integer :: nlevmax
  parameter (nlevmax=100)
 
  ! Maximum number of input files (dates, length of trajectories)
  integer :: ndatmax
  parameter (ndatmax=500)
 
  ! Numerical epsilon (for float comparison)
  real :: eps
  parameter (eps=0.001)
 
  ! Conversion factors
  real :: pi180                                   ! deg -> rad
  parameter (pi180=3.14159/180.)
  real :: deg2km                                  ! deg -> km (at equator)
  parameter (deg2km=111.)
  real :: pir
  parameter (pir=255032235.95489)                 ! 2*Pi*R^2 Bojan
 
  ! Prefix for primary and secondary fields
  character :: charp
  character :: chars
  parameter (charp='P')
  parameter (chars='S')
 
  ! --------------------------------------------------------------------
  ! Declaration of variables
  ! --------------------------------------------------------------------
 
  ! Input and output format for trajectories (see iotra.f)
  integer :: inpmode
  integer :: outmode
 
  ! Input parameters
  character(len=80) :: inpfile         ! Input trajectory file
  character(len=80) :: outfile         ! Output trajectory file
  integer :: ntra                      ! Number of trajectories
  integer :: ncol                      ! Number of columns (including time, lon, lat, p)
  integer :: ntim                      ! Number of times per trajectory
  integer :: ntrace0                   ! Number of trace variables
  character(len=80) :: tvar0(200)      ! Tracing variable (with mode specification)
  character(len=80) :: tvar(200)       ! Tracing variable name (only the variable)
  character(len=80) :: tfil(200)       ! Filename prefix
  real :: fac(200)                     ! Scaling factor
  real :: shift_val(200)               ! Shift in space and time relative to trajectory position
  character(len=80) :: shift_dir(200)  ! Direction of shift
  character(len=80) :: shift_rel(200)  ! Operator/relator for variable
  integer :: compfl(200)               ! Computation flag (1=compute)
  integer :: numdat                    ! Number of input files
  character(len=11) :: dat(ndatmax)    ! Dates of input files
  real :: timeinc                      ! Time increment between input files
  real :: tst                          ! Time shift of start relative to first data file
  real :: ten                          ! Time shift of end relatiev to first data file
  character(len=20) :: startdate       ! First time/date on trajectory
  character(len=20) :: enddate         ! Last time/date on trajectory
  integer :: ntrace1                   ! Count trace and additional variables
  character(len=80) :: timecheck       ! Either 'yes' or 'no'
  character(len=80) :: intmode         ! Interpolation mode ('normal', 'nearest') Bojan ('clustering','circle_avg','circle_max','circle_min')
 
  ! 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(len=80) :: varsinp(200)                      ! Field names for input trajectory
  character(len=80) :: varsint(200)                      ! Field names for internal trajectory
  character(len=80) :: varsout(200)                      ! 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(200)                                    ! Flag whether field is ready
 
  ! Meteorological fields
  real,allocatable, dimension (:)     :: spt0,spt1       ! Surface pressure
  real,allocatable, dimension (:)     :: p3t0,p3t1       ! 3d-pressure
  real,allocatable, dimension (:)     :: f3t0,f3t1       ! 3d field for tracing
  character(len=80) :: svars(100)                        ! List of variables on S file
  character(len=80) :: pvars(100)                        ! List of variables on P file
  integer :: n_svars                                     ! Number of variables on S file
  integer :: n_pvars                                     ! Number of variables on P file
 
  ! Grid description
  real :: pollon,pollat   ! Longitude/latitude of pole
  real :: 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
 
  ! Auxiliary variables
  integer :: i,j,k,l,m,n
  real :: rd
  character(len=80) :: filename,varname
  real :: time0,time1,reltpos
  integer :: itime0,itime1
  integer :: stat
  real :: tstart
  integer :: iloaded0,iloaded1
  real :: f0
  real :: frac
  real :: tload,tfrac
  integer :: isok
  character :: ch
  integer :: ind
  integer :: ind1,ind2,ind3,ind4,ind5
  integer :: ind6,ind7,ind8,ind9,ind0
  integer :: noutside
  real :: delta
  integer :: itrace0
  character(len=80) :: string
  integer err_c1,err_c2,err_c3
 
  ! Bojan
  real    :: dist,circlesum,circlemax,circlemin,circleavg,radius       ! distance (great circle), sum/max/min/avg in circle, radius of circle
  integer :: ist,jst,kst,sp,ml,mr,nd,nu                                ! ijk in stack, sp=stack counter, ml (left), mr (right), nd (down), nu (up)
  integer :: lci,lcj,xindb,xindf                                       ! label count i and j, xind back, xind forward
  integer :: yindb,yindf,pindb,pindf                                   ! yind back, yind forward, pind back, pind forward
  integer :: pvpos,thpos                                               ! position of variables PV and TH in trajectory
  real    :: tropo_pv,tropo_th                                         ! values of PV and TH at dynamical tropopause
  integer, allocatable, dimension (:)   :: stackx,stacky               ! lon/lat of stack
  integer, allocatable, dimension (:)   :: lblcount                    ! counter for label
  integer, allocatable, dimension (:,:) :: connect                     ! array that keeps track of the visited grid points
  real, allocatable, dimension (:)      :: lbl                         ! label
  real, allocatable, dimension (:)      :: circlelon,circlelat         ! value of f, lon and lat in circle
  real, allocatable, dimension (:)      :: circlef,circlearea          ! value of f, lon and lat in circle
  real, allocatable, dimension (:)      :: longrid,latgrid             ! arrays of longitude and latitude of grid points
  ! Bojan
 
  ! Externals
  real :: int_index4
  external                               int_index4
  real :: sdis ! Bojan: need function sdis (calculates great circle distance)
  external                               sdis ! Bojan: need function sdis
 
  ! --------------------------------------------------------------------
  ! Start of program, Read parameters, get grid parameters
  ! --------------------------------------------------------------------
 
  ! Write start message
  print*,'========================================================='
  print*,'              *** START OF PROGRAM TRACE ***'
  print*
 
  ! Read parameters
  open(10,file='trace.param')
   read(10,*) inpfile
   read(10,*) outfile
   read(10,*) startdate
   read(10,*) enddate
   read(10,*) fbflag
   read(10,*) numdat
   if ( fbflag.eq.1) then
      do i=1,numdat
         read(10,'(a11)') dat(i)
      enddo
   else
      do i=numdat,1,-1
         read(10,'(a11)') dat(i)
      enddo
   endif
   read(10,*) timeinc
   read(10,*) tst
   read(10,*) ten
   read(10,*) ntra
   read(10,*) ntim
   read(10,*) ncol
   read(10,*) ntrace0
   do i=1,ntrace0
      read(10,*) tvar(i), fac(i), compfl(i), tfil(i)
   enddo
   read(10,*) n_pvars
   do i=1,n_pvars
      read(10,*) pvars(i)
   enddo
   read(10,*) n_svars
   do i=1,n_svars
      read(10,*) svars(i)
   enddo
   read(10,*) timecheck
   read(10,*) intmode
   read(10,*) radius
   read(10,*) tropo_pv
   read(10,*) tropo_th
  close(10)
 
  ! Bojan: error if radius < 0
  if (((intmode .eq. "circle_avg") .or. (intmode .eq. "circle_min") .or. (intmode .eq. "circle_max")) .and. (radius .lt. 0)) then
     print*,'ERROR (circle): radius < 0!'
     stop
  endif
 
  ! Remove commented tracing fields
  itrace0 = 1
  do while ( itrace0.le.ntrace0)
     string = tvar(itrace0)
     if ( string(1:1).eq.'#' ) then
        do i=itrace0,ntrace0-1
           tvar(i)   = tvar(i+1)
           fac(i)    = fac(i+1)
           compfl(i) = compfl(i+1)
           tfil(i)   = tfil(i+1)
        enddo
        ntrace0 = ntrace0 - 1
     else
        itrace0 = itrace0 + 1
     endif
  enddo
 
  ! Save the tracing variable (including all mode specifications)
  do i=1,ntrace0
     tvar0(i) = tvar(i)
  enddo
 
  ! Set the formats of the input and output files
  call mode_tra(inpmode,inpfile)
  if (inpmode.eq.-1) inpmode=1
  call mode_tra(outmode,outfile)
  if (outmode.eq.-1) outmode=1
 
  ! Convert time shifts <tst,ten> from <hh.mm> into fractional time
  call hhmm2frac(tst,frac)
  tst = frac
  call hhmm2frac(ten,frac)
  ten = frac
 
  ! Set the time for the first data file (depending on forward/backward mode)
  if (fbflag.eq.1) then
    tstart = -tst
  else
    tstart = tst
  endif
 
  ! Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)
  ! The negative <-fid> of the file identifier is used as a flag for parameter retrieval
  filename = charp//dat(1)
  varname  = 'U'
  call input_open (fid,filename)
  call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,tstart,pollon,pollat,rd,rd,nz,rd,rd,rd,timecheck)
  call input_close(fid)
 
  ! Allocate memory for some meteorological arrays
  allocate(spt0(nx*ny),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array spt0 ***'   ! Surface pressure
  allocate(spt1(nx*ny),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array spt1 ***'
  allocate(p3t0(nx*ny*nz),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array p3t0 ***'   ! Pressure
  allocate(p3t1(nx*ny*nz),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array p3t1 ***'
  allocate(f3t0(nx*ny*nz),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array p3t0 ***'   ! Tracing field
  allocate(f3t1(nx*ny*nz),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array p3t1 ***'
 
  ! Get memory for trajectory arrays
  allocate(trainp(ntra,ntim,ncol),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array tra      ***'
 
  ! Bojan
  ! allocate memory for clustering mode 
  if (intmode .eq. 'clustering') then
     allocate(lbl(8),stat=stat)
     if (stat.ne.0) print*,'*** error allocating array lbl      ***'
     allocate(lblcount(5),stat=stat)
     if (stat.ne.0) print*,'*** error allocating array lblcount ***'
  endif
  ! allocate memory for circle mode
  if ( (intmode.eq.'circle_avg') .or. (intmode.eq.'circle_min') .or. (intmode.eq.'circle_max') ) then
     allocate(connect(nx,ny),stat=stat)
     if (stat.ne.0) print*,'*** error allocating connect ***'
     allocate(stackx(nx*ny),stat=stat)
     if (stat.ne.0) print*,'*** error allocating stackx ***'
     allocate(stacky(nx*ny),stat=stat)
     if (stat.ne.0) print*,'*** error allocating stacky ***'
     allocate(circlelon(nx*ny),stat=stat)
     if (stat.ne.0) print*,'*** error allocating circlelon ***'
     allocate(circlelat(nx*ny),stat=stat)
     if (stat.ne.0) print*,'*** error allocating circlelat ***'
     allocate(circlef(nx*ny),stat=stat)
     if (stat.ne.0) print*,'*** error allocating circlef ***'
     allocate(circlearea(nx*ny),stat=stat)
     if (stat.ne.0) print*,'*** error allocating circlearea ***'
     allocate(longrid(nx),stat=stat)
     if (stat.ne.0) print*,'*** error allocating longrid ***'
     allocate(latgrid(ny),stat=stat)
     if (stat.ne.0) print*,'*** error allocating latgrid ***'
     do m=1,nx
        longrid(m)=xmin+dx*(m-1)
     enddo
     do n=1,ny
        latgrid(n)=ymin+dy*(n-1)
     enddo
  endif
  ! Bojan
 
  ! 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
 
  ! Set logical flag for periodic data set (hemispheric or not)
  hem = 0
  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
 
  ! Write some status information
  print*,'---- INPUT PARAMETERS -----------------------------------'
  print*
  print*,'  Input trajectory file  : ',trim(inpfile)
  print*,'  Output trajectory file : ',trim(outfile)
  print*,'  Format of input file   : ',inpmode
  print*,'  Format of output file  : ',outmode
  print*,'  Forward/backward       : ',fbflag
  print*,'  #tra                   : ',ntra
  print*,'  #col                   : ',ncol
  print*,'  #tim                   : ',ntim
  print*,'  No time check          : ',trim(timecheck)
  print*,'  Interpolation mode     : ',trim(intmode)
  ! Bojan
  if (trim(intmode) .eq. "clustering") then
     print*,'  Tropopause PV [pvu]    : ',tropo_pv
     print*,'  Tropopause TH [K]      : ',tropo_th
  endif
  do i=1,ntrace0
     if (compfl(i).eq.0) then
        print*,'  Tracing field          : ', trim(tvar(i)), fac(i), ' 0 ', tfil(i)
     else
        print*,'  Tracing field          : ', trim(tvar(i)), fac(i), '  1 ', tfil(i)
     endif
  enddo
  print*
  print*,'---- INPUT DATA FILES -----------------------------------'
  print*
  call frac2hhmm(tstart,tload)
  print*,'  Time of 1st data file  : ',tload
  print*,'  #input files           : ',numdat
  print*,'  time increment         : ',timeinc
  call frac2hhmm(tst,tload)
  print*,'  Shift of start         : ',tload
  call frac2hhmm(ten,tload)
  print*,'  Shift of end           : ',tload
  print*,'  First/last input file  : ',trim(dat(1)), ' ... ', trim(dat(numdat))
  print*,'  Primary variables      : ',trim(pvars(1))
  do i=2,n_pvars
     print*,'                         : ',trim(pvars(i))
  enddo
  if ( n_svars.ge.1 ) then
     print*,'  Secondary variables    : ',trim(svars(1))
     do i=2,n_svars
        print*,'                         : ',trim(svars(i))
     enddo
  endif
  print*
  print*,'---- CONSTANT GRID PARAMETERS ---------------------------'
  print*
  print*,'  xmin,xmax     : ',xmin,xmax
  print*,'  ymin,ymax     : ',ymin,ymax
  print*,'  dx,dy         : ',dx,dy
  print*,'  pollon,pollat : ',pollon,pollat
  print*,'  nx,ny,nz      : ',nx,ny,nz
  print*,'  per, hem      : ',per,hem
  print*
 
  ! --------------------------------------------------------------------
  ! Load the input trajectories
  ! --------------------------------------------------------------------
 
  ! Read the input trajectory file
  call ropen_tra(fid,inpfile,ntra,ntim,ncol,reftime,varsinp,inpmode)
  call read_tra (fid,trainp,ntra,ntim,ncol,inpmode)
  call close_tra(fid,inpmode)
 
  ! Check that first four columns correspond to time,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'
 
  ! Write some status information of the input trajectories
  print*,'---- INPUT TRAJECTORIES ---------------------------------'
  print*
  print*,' Start date             : ',trim(startdate)
  print*,' End date               : ',trim(enddate)
  print*,' Reference time (year)  : ',reftime(1)
  print*,'                (month) : ',reftime(2)
  print*,'                (day)   : ',reftime(3)
  print*,'                (hour)  : ',reftime(4)
  print*,'                (min)   : ',reftime(5)
  print*,' Time range (min)       : ',reftime(6)
  do i=1,ncol
     print*,' Var                    :',i,trim(varsinp(i))
  enddo
  print*
 
  ! Check that first time is 0 - otherwise the tracing will produce
  ! wrong results because later in the code only absolute times are
  ! considered: <itime0   = int(abs(tfrac-tstart)/timeinc) + 1>. This
  ! will be changed in a future version.
  if ( abs( trainp(1,1,1) ).gt.eps ) then
     print*,' ERROR: First time of trajectory must be 0, i.e. '
     print*,'     correspond to the reference date. Otherwise'
     print*,'     the tracing will give wrong results... STOP'
     stop
  endif
 
  ! --------------------------------------------------------------------
  ! Check dependencies for trace fields which must be calculated
  ! --------------------------------------------------------------------
 
  ! Set the counter for extra fields
  ntrace1 = ntrace0
 
  ! Loop over all tracing variables
  i = 1
  do while (i.le.ntrace1)
 
     ! Skip fields which must be available on the input files
     if (i.le.ntrace0) then
        if (compfl(i).eq.0) goto 100
     endif
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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)
 
     ! 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
 
     ! Exit point for handling additional fields
 100   continue
     i = i + 1
 
  enddo
 
  ! 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
 
  ! Bojan: check that PV and TH are on trajectory
  if (intmode .eq. 'clustering') then
     pvpos=-1
     thpos=-1
     do i=1,ncol+ntrace1
        if (varsint(i) .eq. 'TH') then
        thpos=i
        print*,'Clustering: Found TH at position:',thpos
        endif
        if (varsint(i) .eq. 'PV') then
        pvpos=i
        print*,'Clustering: Found PV at position:',pvpos
        endif
     enddo
     if (thpos .eq. -1) then
        print*,'WARNING (clustering): Did not find TH'
        stop
     endif
     if (pvpos .eq. -1) then
        print*,'WARNING (clustering): Did not find PV'
        stop
     endif
  endif
  ! Bojan
 
  ! Split the tracing variables
  do i=1,ntrace0
    call splitvar(tvar(i),shift_val(i),shift_dir(i),shift_rel(i) )
  enddo
 
 
  ! Split the variable name and set flags
  do i=ntrace0+1,ntrace1
 
     ! Set the scaling factor
     fac(i) = 1.
 
     ! Set the base variable name, the shift and the direction
     call splitvar(tvar(i),shift_val(i),shift_dir(i),shift_rel(i) )
 
     ! Set the prefix of the file name for additional fields
     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
 
     ! 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
 
  ! 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.'HPA(ABS)').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'    ).and. & 
          ( shift_dir(i).ne.'PMIN'    ).and. &
          ( shift_dir(i).ne.'PMAX'    ) ) then
        print*,' ERROR: shift mode ',trim(shift_dir(i)), ' not supported'
        stop
     endif
  enddo
 
  ! 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,a8,3x,a4,i5)') i,' : ',trim(tvar(i)), &
             shift_val(i),trim(shift_dir(i)),trim(shift_rel(i)),tfil(i),compfl(i)
     else
        write(*,'(i4,a4,a8,10x,8x,3x,a4,i5)') &
             i,' : ',trim(tvar(i)),tfil(i),compfl(i)
     endif
  enddo
 
  ! --------------------------------------------------------------------
  ! Prepare the internal and output trajectories
  ! --------------------------------------------------------------------
 
  ! Allocate memory for internal trajectories
  allocate(traint(ntra,ntim,ncol+ntrace1),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array traint   ***'
 
  ! Copy input to output trajectory
  do i=1,ntra
     do j=1,ntim
        do k=1,ncol
           traint(i,j,k)=trainp(i,j,k)
        enddo
     enddo
  enddo
 
  ! 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
 
  ! --------------------------------------------------------------------
  ! Trace the fields (fields available on input files)
  ! --------------------------------------------------------------------
 
  print*
  print*,'---- TRACING FROM PRIMARY AND SECONDARY DATA FILES ------'
 
  ! Loop over all tracing fields
  do i=1,ntrace1
 
      ! Skip fields which must be computed (compfl=1), will be handled later
      if (compfl(i).ne.0)  goto 110
 
      ! Write some status information
      print*
      print*,' Now tracing             : ', trim(tvar(i)),shift_val(i),trim(shift_dir(i)),compfl(i),' ',trim(tfil(i))
 
      ! Set the flag for ready field/column
      fok(ncol+i) = 1
 
      ! Reset flags for load manager
      iloaded0 = -1
      iloaded1 = -1
 
      ! Reset the counter for fields outside domain
      noutside = 0
      err_c1   = 0
      err_c2   = 0
      err_c3   = 0
 
      ! Loop over all times
      do j=1,ntim
 
         ! Convert trajectory time from hh.mm to fractional time
         call hhmm2frac(trainp(1,j,1),tfrac)
 
         ! 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
 
         ! Get the times which are needed
         itime0   = int(abs(tfrac-tstart)/timeinc) + 1
         time0    = tstart + fbflag * real(itime0-1) * timeinc
         itime1   = itime0 + 1
         time1    = time0 + fbflag * timeinc
         if ( itime1.gt.numdat ) then
            itime1 = itime0
            time1  = time0
         endif
 
         ! Load manager: Check whether itime0 can be copied from itime1
         if ( itime0.eq.iloaded1 ) then
            f3t0     = f3t1
            p3t0     = p3t1
            spt0     = spt1
            iloaded0 = itime0
         endif
 
         ! Load manager: Check whether itime1 can be copied from itime0
         if ( itime1.eq.iloaded0 ) then
            f3t1     = f3t0
            p3t1     = p3t0
            spt1     = spt0
            iloaded1 = itime1
         endif
 
         ! Load manager:  Load first time (tracing variable and grid)
         if ( itime0.ne.iloaded0 ) then
 
            filename = trim(tfil(i))//trim(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)
            mdv = -999.
            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
 
         ! Load manager: Load second time (tracing variable and grid)
         if ( itime1.ne.iloaded1 ) then
 
            filename = trim(tfil(i))//trim(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
 
         ! Loop over all trajectories
         do k=1,ntra
 
            ! Set the horizontal position where to interpolate to
            x0       = traint(k,j,2)                          ! Longitude
            y0       = traint(k,j,3)                          ! Latitude
 
            ! Set the vertical position where to interpolate to
            if ( nz.gt.1 ) then
               p0 = traint(k,j,4)                             ! Pressure (3D tracing)
            else
               p0 = 1050.                                     ! Lowest level (2D tracing)
            endif
 
            ! Set negative pressures to mdv
            if (p0.lt.0.) then
	        f0 = mdv
                goto 109
            endif
 
            ! Set the relative time
            call hhmm2frac(traint(k,j,1),tfrac)
            reltpos0 = fbflag * (tfrac-time0)/timeinc
 
            ! 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.'HPA(ABS)' ) then                 ! HPA(ABS)
               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
 
            ! 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.
 
            ! Handle pole problems for hemispheric data (taken from caltra.f)
            if ((hem.eq.1).and.(y0.gt.90.)) then
               print*,'WARNING: y0>90 ',y0,' => setting to 180-y0 ',180.-y0
               y0=180.-y0
               x0=x0+per/2.
            endif
            if ((hem.eq.1).and.(y0.lt.-90.)) then
               print*,'WARNING: y0<-90 ',y0,' => setting to -180-y0 ',-180.-y0
               y0=-180.-y0
               x0=x0+per/2.
            endif
 
            ! Get the index where to interpolate (x0,y0,p0)
            if ( (abs(x0-mdv).gt.eps).and. (abs(y0-mdv).gt.eps) ) then
               call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0, &
                    p3t0,p3t1,spt0,spt1,3,nx,ny,nz,xmin,ymin,dx,dy,mdv)
            else
               xind = mdv
               yind = mdv
               pind = mdv
            endif
 
           ! Check if point is within grid (keep indices if ok)
            if ( (xind.ge.1.).and.(xind.le.real(nx)).and. &
                 (yind.ge.1.).and.(yind.le.real(ny)).and. &
                 (pind.ge.1.).and.(pind.le.real(nz)) ) then
		 xind = xind
		 yind = yind
		 pind = pind
 
            ! Check if pressure is outside, but rest okay => adjust to lowest or highest level
            elseif ( (xind.ge.1.).and.(xind.le.real(nx)).and. (yind.ge.1.).and.(yind.le.real(ny)) ) then ! only vertical problem
 
               if ( pind.gt.nz ) then ! pressure too low, index too high
                 err_c1 = err_c1 + 1
                 if ( err_c1.lt.10 ) then
                    write(*,'(a,f7.3,a)') ' WARNING: pressure too low (pind = ',pind,') => adjusted to highest level (pind=nz.)'
                    print*,'(x0,y0,p0)=',x0,y0,p0
                    pind = real(nz)
                 elseif ( err_c1.eq.10 ) then
                    print*,' WARNING: more pressures too low -> adjusted to highest level '
                    pind = real(nz)
                 else
                    pind = real(nz)
                 endif
                 
               elseif (pind.lt.1.) then ! pressure too high, index too low
                 err_c2 = err_c2 + 1
                 if ( err_c2.lt.10 ) then
                    write(*,'(a,f7.3,a)') ' WARNING: pressure too high (pind = ',pind,') => adjusted to lowest level, (pind=1.)'
                    print*,'(x0,y0,p0)=',x0,y0,p0
                    pind = 1.
                 elseif ( err_c2.eq.10 ) then
                    print*,' WARNING: more pressures too high -> adjusted to lowest level '
                    pind = 1.
                 else
                    pind = 1.
                 endif
 
              endif
 
            ! Grid point is outside!
            else
               
               err_c3 = err_c3 + 1
               if ( err_c3.lt.10 ) then
                  print*,'ERROR: point is outside grid (horizontally)'
                  print*,'   Trajectory # ',k
                  print*,'   Position     ',x0,y0,p0
                  print*,'  (xind,yind):  ',xind,yind
                  xind          = mdv
                  yind          = mdv
                  pind          = mdv
                  traint(k,j,2) = mdv  
                  traint(k,j,3) = mdv  
                  traint(k,j,4) = mdv  
               elseif ( err_c3.eq.10 ) then
                  print*,'ERROR: more points outside grid (horizontally)'
                  xind          = mdv
                  yind          = mdv
                  pind          = mdv
                  traint(k,j,2) = mdv  
                  traint(k,j,3) = mdv  
                  traint(k,j,4) = mdv  
               else
                  xind          = mdv
                  yind          = mdv
                  pind          = mdv
                  traint(k,j,2) = mdv  
                  traint(k,j,3) = mdv  
                  traint(k,j,4) = mdv  
               endif
 
            endif
 
            ! ------------------------ NEAREST mode ------------------------------- 
            ! Interpolate to nearest grid point
            if ( intmode.eq.'nearest') then
 
               xind = real( nint(xind) )
               yind = real( nint(yind) )
               pind = real( nint(pind) )
 
               if ( xind.lt.1.  ) xind = 1.
               if ( xind.gt.nx  ) xind = real(nx)
               if ( yind.lt.1.  ) yind = 1.
               if ( yind.gt.ny  ) yind = real(ny)
 
               if ( pind.lt.1.  ) pind = 1.
               if ( pind.gt.nz  ) pind = real(nz)
 
               if ( abs(reltpos0).ge.eps ) then
                  print*,'ERROR (nearest): reltpos != 0',reltpos0
                  stop
               endif
 
               ! interpolate
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
 
            ! ------------------------ end NEAREST mode ------------------------------- 
 
            ! ------------------------ CLUSTERING mode ------------------------------- Bojan
            elseif (intmode.eq.'clustering') then
               if (varname.ne.'LABEL' ) then
                  print*,'ERROR (clustering): varname is not LABEL'
                  stop
               endif
 
            ! Get indices of box around the point
            xindb=floor(xind)
            xindf=ceiling(xind)
            yindb=floor(yind)
            yindf=ceiling(yind)
            pindb=floor(pind)
            pindf=ceiling(pind)
 
            ! Make sure all points are within grid
            if ( xindb.lt.1 ) xindb = 1
            if ( xindf.lt.1 ) xindf = 1
            if ( xindb.gt.nx ) xindb = nx
            if ( xindf.gt.nx ) xindf = nx
            if ( yindb.lt.1 ) yindb = 1
            if ( yindf.lt.1 ) yindf = 1
            if ( yindb.gt.ny ) yindb = ny
            if ( yindf.gt.ny ) yindf = ny
            if ( pindb.lt.1 ) pindb = 1
            if ( pindf.lt.1 ) pindf = 1
            if ( pindb.gt.nz ) pindb = nz
            if ( pindf.gt.nz ) pindf = nz
 
            ! Shift one point if they are equal
            if ( xindb.eq.xindf ) then
               if ( xindf.eq.nx ) then
                  xindb=nx-1
               else
                  xindf=xindb+1
               endif
            endif
            if ( yindb.eq.yindf ) then
               if ( yindf.eq.ny ) then
                  yindb=ny-1
               else
                  yindf=yindb+1
               endif
            endif
            if ( pindb.eq.pindf ) then
               if ( pindf.eq.nz ) then
                  pindb=nz-1
               else
                  pindf=pindb+1
               endif
            endif
            ! Give warnings and stop if problems occur
            if ( xindb.eq.xindf ) then
               print*,'ERROR (clustering): xindb=xindf'
               print*,xind,xindb,xindf
               stop
            endif
            if ( yindb.eq.yindf ) then
               print*,'ERROR (clustering): yindb=yindf'
               print*,yind,yindb,yindf
               stop
            endif
            if ( pindb.eq.pindf ) then
               print*,'ERROR (clustering): pindb=pindf'
               print*,pind,pindb,pindf
               stop
            endif
            if ( ( xindb.lt.1 ).or.( xindf.gt.nx ) ) then
               print*,'ERROR (clustering): xindb/f outside'
               print*,xind,xindb,xindf
               stop
            endif
            if ( ( yindb.lt.1 ).or.( yindf.gt.ny ) ) then
               print*,'ERROR (clustering): yindb/f outside'
               print*,yind,yindb,yindf
               stop
            endif
            if ( ( pindb.lt.1 ).or.( pindf.gt.nz ) ) then
               print*,'ERROR (clustering): pindb/f outside'
               print*,pind,pindb,pindf
               stop
            endif
            if ( abs(reltpos0).ge.eps ) then
               print*,'ERROR (clustering): reltpos != 0',reltpos0
               stop
            endif
 
            ! Get Value in Box
            lblcount=(/0,0,0,0,0/)
 
            lbl(1) = f3t0( xindb + nx*(yindb-1) + nx*ny*(pindb-1) )
            lbl(2) = f3t0( xindf + nx*(yindb-1) + nx*ny*(pindb-1) )
            lbl(3) = f3t0( xindb + nx*(yindf-1) + nx*ny*(pindb-1) )
            lbl(4) = f3t0( xindf + nx*(yindf-1) + nx*ny*(pindb-1) )
            lbl(5) = f3t0( xindb + nx*(yindb-1) + nx*ny*(pindf-1) )
            lbl(6) = f3t0( xindf + nx*(yindb-1) + nx*ny*(pindf-1) )
            lbl(7) = f3t0( xindb + nx*(yindf-1) + nx*ny*(pindf-1) )
            lbl(8) = f3t0( xindf + nx*(yindf-1) + nx*ny*(pindf-1) )
 
            ! Count the number of times every label appears
            do lci=1,5
               do lcj=1,8
                  if ( abs(lbl(lcj)-lci).lt.eps ) then
                     lblcount(lci)=lblcount(lci)+1
                  endif
               enddo
            enddo
 
            ! Set to -9 to detect if no label was assigned in the end
            f0=-9
 
            ! Stratosphere (PV)
            if ( abs(traint(k,j,pvpos)) .ge. tropo_pv ) then
               if ( (lblcount(2).ge.lblcount(3)).and. (lblcount(2).ge.lblcount(5)) ) then
                  f0=2
               elseif ( lblcount(3).ge.lblcount(5) ) then
                  f0=3
               elseif ( lblcount(5).gt.lblcount(3) ) then
                  f0=5
               endif
            endif
 
            ! Troposphere (PV)
            if ( abs(traint(k,j,pvpos)) .lt. tropo_pv ) then
               if ( lblcount(1).ge.lblcount(4) ) then
               f0=1
               elseif ( lblcount(4).gt.lblcount(1) ) then
               f0=4
               endif
            endif
 
            ! Stratosphere (TH)
            if ( traint(k,j,thpos) .ge. tropo_th ) then
               f0=2
            endif
 
            if (f0.eq.-9) then
               print*,'ERROR (Clustering): No label assigned!'
               stop
            endif
            ! ------------------------ end CLUSTERING mode -------------------------------
 
            ! ------------------------ CIRCLE modes ------------------------------- Bojan
            ! elseif (not clustering but one of the possible circle modes)
            elseif ( (intmode.eq.'circle_avg') .or. (intmode.eq.'circle_min') .or. (intmode.eq.'circle_max') ) then
 
            ! reset arrays for this point
            connect=0
            stackx=0
            stacky=0
            circlelon=0
            circlelat=0
            circlef=0
            circlearea=0
 
            ! Get indices of one coarse grid point within search radius (nint=round to next integer)
            if ( sdis(x0,y0,longrid(nint(xind)),latgrid(nint(yind))) .gt. radius) then
               print*,'ERROR (circle): Search radius is too small... (1). r =',radius
               print*,'Distance to nint grid point (minimum search radius)=',sdis(x0,y0,longrid(nint(xind)),latgrid(nint(yind)))
               stop
            endif
 
            ! Initialize stack with nint(xind),nint(yind)
            kst=0 ! counts the number of points in circle
            stackx(1)=nint(xind)
            stacky(1)=nint(yind)
            sp=1 ! stack counter
            do while (sp.ne.0)
 
            ! Get an element from stack
             ist=stackx(sp)
             jst=stacky(sp)
             sp=sp-1
 
            ! Get distance from reference point
             dist=sdis(x0,y0,longrid(ist),latgrid(jst))
 
            ! Check whether distance is smaller than search radius: connected
             if (dist.lt.radius) then
 
            ! Increase total stack index
              kst=kst+1
              circlelon(kst)=longrid(ist)
              circlelat(kst)=latgrid(jst)
  
            ! Interpolate field to position of point (interpolation in time!) 
              circlef(kst) = int_index4(f3t0,f3t1,nx,ny,nz,real(ist),real(jst),pind,reltpos0,mdv)
 
            ! Calculate area of point (for circle_avg mode only)
              if ( intmode .eq. 'circle_avg' ) then
                 circlearea(kst) = pir/(nx-1)*(sin(pi180*abs(circlelat(kst)))-sin(pi180*(abs(circlelat(kst))-dy)))
              endif
 
            ! Mark this point as visited
              connect(ist,jst)=1
 
            ! Get coordinates of neighbouring points and implement periodicity
              mr=ist+1
              if (mr.gt.nx) mr=1
              ml=ist-1
              if (ml.lt.1) ml=nx
              nu=jst+1
              if (nu.gt.ny) nu=ny
              nd=jst-1
              if (nd.lt.1) nd=1
 
            ! Update stack with neighbouring points
              if (connect(mr,jst).ne. 1) then
                 connect(mr,jst)=1
                 sp=sp+1
                 stackx(sp)=mr
                 stacky(sp)=jst
              endif
              if (connect(ml,jst).ne. 1) then
                 connect(ml,jst)=1
                 sp=sp+1
                 stackx(sp)=ml
                 stacky(sp)=jst
              endif
              if (connect(ist,nd).ne. 1) then
                 connect(ist,nd)=1
                 sp=sp+1
                 stackx(sp)=ist
                 stacky(sp)=nd
              endif
              if (connect(ist,nu).ne. 1) then
                 connect(ist,nu)=1
                 sp=sp+1
                 stackx(sp)=ist
                 stacky(sp)=nu
              endif
 
             endif ! endif radius is smaller => end of updating stack
 
            end do ! end working on stack 
 
            if (kst.ge.1) then
               ! Choose output depending on intmode
               if ( intmode .eq. 'circle_avg' ) then
                  ! calculate area-weighted average of f in circle
                  circlesum=0.
                  do l=1,kst
                     circlesum=circlesum+circlef(l)*circlearea(l)
                  enddo
                  circleavg=circlesum/sum(circlearea(1:kst))
                  !print*,'area-weighted average of f in circle=',circleavg
                  f0=circleavg
               elseif ( intmode .eq. 'circle_min' ) then
                  ! calculate minimum in circle
                  circlemin=circlef(1)
                  do l=1,kst
                     if (circlef(l) .lt. circlemin) then
                        circlemin=circlef(l)
                     endif
                  enddo
                  !print*,'minimum of f in circle=',circlemin       
                  f0=circlemin
               elseif ( intmode .eq. 'circle_max' ) then             
                  ! calculate maximum in circle
                  circlemax=circlef(1)
                  do l=1,kst
                     if (circlef(l) .gt. circlemax) then
                        circlemax=circlef(l)
                     endif
                  enddo
                  !print*,'maximum of f in circle=',circlemax
                  f0=circlemax
               else 
                  print*,'ERROR (circle): intmode not valid!'
                  stop
               endif
            else
               print*,'ERROR (circle): Search radius is too small... (2). r =',radius
               stop
            endif
 
            ! ------------------------ end CIRCLE modes -------------------------------
 
            ! ------------------------ NORMAL mode -------------------------------
            else ! not clustering nor circle: NORMAL mode
 
            ! Check if point is within grid
!            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
!
            ! Do the interpolation: everthing is ok
            if ( ( shift_dir(i).ne.'PMIN' ).and.( shift_dir(i).ne.'PMAX' ) ) then
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
               
            elseif  ( shift_dir(i).eq.'PMIN' ) then
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
               
               ! Handle > and < operator
               if ( (shift_rel(i).eq.'>').and.(f0.gt.shift_val(i)) ) then
                 do while ( (f0.gt.shift_val(i)).and.(pind.lt.nz) )
                     pind = pind + 0.1
                     f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
                 enddo
                 f0 = int_index4(p3t0,p3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
               
               elseif ( (shift_rel(i).eq.'<').and.(f0.lt.shift_val(i)) ) then
                 do while ( (f0.lt.shift_val(i)).and.(pind.lt.nz) )
                     pind = pind + 0.1
                     f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
                 enddo
                 f0 = int_index4(p3t0,p3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
                 
              else
                 f0 = mdv
              endif
              
              
           elseif  ( shift_dir(i).eq.'PMAX' ) then
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
               
               ! Handle > and < operator
               if ( (shift_rel(i).eq.'>').and.(f0.gt.shift_val(i)) ) then
                 do while ( (f0.gt.shift_val(i)).and.(pind.gt.1) )
                     pind = pind - 0.1
                     f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
                 enddo
                 f0 = int_index4(p3t0,p3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
              
               elseif ( (shift_rel(i).eq.'<').and.(f0.lt.shift_val(i)) ) then
                 do while ( (f0.lt.shift_val(i)).and.(pind.gt.1) )
                     pind = pind - 0.1
                     f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
                 enddo
                 f0 = int_index4(p3t0,p3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
              
              else
                 f0 = mdv
              endif
               
          endif  
 
!            ! Check if pressure is outside, but rest okay: adjust to lowest or highest level
!            elseif ( (xind.ge.1.).and.(xind.le.real(nx)).and. (yind.ge.1.).and.(yind.le.real(ny)) ) then ! only vertical problem
!               if ( pind.gt.nz ) then ! pressure too low, index too high
!                 pind = real(nz)
!                 print*,' Warning: pressure too low -> adjusted to highest level, pind=nz.'
!                 print*,'(x0,y0,p0)=',x0,y0,p0
!               elseif (pind.lt.1.) then ! pressure too high, index too low
!                 pind = 1.
!                 print*,' Warning: pressure too high -> adjusted to lowest level, pind=1.'
!                 print*,'(x0,y0,p0)=',x0,y0,p0
!              endif
!              f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
 
!            ! Less than 10 outside
!            elseif (noutside.lt.10) then
!               print*,' ',trim(tvar(i)),' @ ',x0,y0,p0,'outside'
!               f0       = mdv
!               noutside = noutside + 1
!
!            ! More than 10 outside
!            elseif (noutside.eq.10) then
!               print*,' ...more than 10 outside...'
!               f0       = mdv
!               noutside = noutside + 1
 
!            ! Else (not everything okay and also not 'tolerated cases') set to missing data
!            else
!               f0       = mdv
!            endif
 
            ! ------------------------ end NORMAL mode -------------------------------
            endif ! end if nearest case
 
           ! Exit for loop over all times -save interpolated value
 109        continue
 
            ! Save the new field
            if ( abs(f0-mdv).gt.eps) then
               traint(k,j,ncol+i) = f0
            else
               traint(k,j,ncol+i) = mdv
            endif
 
         enddo ! end loop over all trajectories
 
      enddo ! end loop over all times
 
      ! Exit point for loop over all tracing variables
 110   continue
 
   enddo ! end loop over all variables
 
  ! --------------------------------------------------------------------
  ! Calculate additional fields along the trajectories
  ! --------------------------------------------------------------------
 
  print*
  print*,'---- CALCULATE ADDITIONAL FIELDS FROM TRAJECTORY TABLE --'
 
  ! Loop over all tracing fields
  do i=ntrace1,1,-1
 
      ! Skip fields which must not be computed (compfl=0)
      if (compfl(i).eq.0) goto 120
 
      ! 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))
 
      ! Loop over trajectories and times
      do j=1,ntra
      do k=1,ntim
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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))
 
         ! 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))
 
         ! 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)
 
            call calc_ABSVORT (traint(j,k,ncol+i), traint(j,k,ind1), &
                 traint(j,k,ind2),traint(j,k,ind3),traint(j,k,ind4))
 
         ! 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))
 
         ! 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))
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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) )
 
         ! 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
 
         ! 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
 
 
         ! Invalid tracing variable
         else
 
            print*,' ERROR: invalid tracing variable ', trim(varsint(i+ncol))
            stop
 
 
         endif
 
      ! End loop over all trajectories and times
      enddo
      enddo
 
      ! Set the flag for a ready field/column
      fok(ncol+i) = 1
 
 
      ! Exit point for loop over all tracing fields
 120   continue
 
   enddo
 
  ! --------------------------------------------------------------------
  ! Write output to output trajectory file
  ! --------------------------------------------------------------------
 
  ! Write status information
  print*
  print*,'---- WRITE OUTPUT TRAJECTORIES --------------------------'
  print*
 
  ! Allocate memory for internal trajectories
  allocate(traout(ntra,ntim,ncol+ntrace0),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array traout   ***'
 
  ! Copy input to output trajectory (apply scaling of output)
  do i=1,ntra
     do j=1,ntim
        do k=1,ncol+ntrace0
           if ( k.le.ncol ) then
              traout(i,j,k) = traint(i,j,k)
           elseif ( abs(traint(i,j,k)-mdv).gt.eps ) then
              traout(i,j,k) = fac(k-ncol) * traint(i,j,k)
           else
              traout(i,j,k) = mdv
           endif
        enddo
     enddo
  enddo
 
  ! Set the variable names for output trajectory
  do i=1,ncol+ntrace0
        varsout(i)      = varsint(i)
  enddo
  do i=1,ntrace0
     if ( (shift_dir(i).eq.'PMIN').or.(shift_dir(i).eq.'PMAX') ) then
        varsout(ncol+i) = trim(tvar(i))//':'//trim(shift_dir(i))
     endif
  enddo
  
  ! Write trajectories
  call wopen_tra(fod,outfile,ntra,ntim,ncol+ntrace0, &
       reftime,varsout,outmode)
  call write_tra(fod,traout ,ntra,ntim,ncol+ntrace0,outmode)
  call close_tra(fod,outmode)
 
  ! Write some status information, and end of program message
  print*
  print*,'---- STATUS INFORMATION --------------------------------'
  print*
  print*,' ok'
  print*
  print*,'              *** END OF PROGRAM TRACE ***'
  print*,'========================================================='
 
 
end program trace
 
 
 
! ******************************************************************
! * SUBROUTINE SECTION                                             *
! ******************************************************************
 
! ------------------------------------------------------------------
! Add a variable to the list if not yet included in this list
! ------------------------------------------------------------------
 
subroutine add2list (varname,list,nlist)
 
  implicit none
 
  ! Declaration of subroutine parameters
  character(len=80) :: varname
  character(len=80) :: list(200)
  integer :: nlist
 
  ! Auxiliray variables
  integer :: i,j
  integer :: isok
 
  ! 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
 
  ! Check for too large number of fields
  if ( nlist.ge.200) then
     print*,' ERROR: too many additional fields for tracing ...'
     stop
  endif
 
end subroutine add2list
 
 
! ------------------------------------------------------------------
! Get the index of a variable in the list
! ------------------------------------------------------------------
 
subroutine list2ind (ind,varname,list,fok,nlist)
 
  implicit none
 
  ! Declaration of subroutine parameters
  integer :: ind
  character(len=80) :: varname
  character(len=80) :: list(200)
  integer :: fok(200)
  integer :: nlist
 
  ! Auxiliray variables
  integer :: i,j
  integer :: isok
 
  ! 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
 
  ! Exit point
 100   continue
 
  ! 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 subroutine list2ind
 
 
! ------------------------------------------------------------------
! Split the variable name into name, shift and direction
! ------------------------------------------------------------------
 
subroutine splitvar (tvar,shift_val,shift_dir,shift_rel)
 
  implicit none
 
  ! Declaration of subroutine parameters
  character(len=80) :: tvar
  real :: shift_val
  character(len=80) :: shift_dir
  character(len=80) :: shift_rel
 
  ! Auxiliary variables
  integer :: i,j
  integer :: icolon,inumber,irelator
  character(len=80) :: name
  character :: ch
  integer      isabsval
 
  ! Save variable name
  name = tvar
  shift_rel = 'nil'
  shift_dir = 'nil'
 
  ! 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
 
  ! If there is a colon, split the variable name
  if ( icolon.ne.0 ) then
 
     tvar = name(1:(icolon-1))
 
     ! Get the index for number
     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
 
     ! If the variable name is e.g. PMIN:UMF>0, the variable to be read
     ! is UMF, the value is 0, and the direction is 'PMIN'
     shift_rel = 'nil'
     if ( (tvar.eq.'PMIN').or.(tvar.eq.'PMAX') ) then
         shift_dir = tvar
         irelator = 0
         do i=icolon+1,80
            ch = name(i:i)
            if ( (ch.eq.'>').or.(ch.eq.'<') ) then
               irelator = i
            endif
         enddo
         if ( irelator.eq.0 ) then
             print*,' ERROR: dont know how to interpret ',trim(name)
             stop
         endif
         tvar = name(icolon+1:irelator-1)
         read(name( (irelator+1):80 ),*) shift_val
         shift_rel = name(irelator:irelator)
           
         goto 90
  
     endif
    
     ! Get the number
     read(name( (icolon+1):(inumber-1) ),*) shift_val
 
     ! Decide whether it is a shift relatiev to trajectory or absolute value
     ! If the number starts with + or -, it is relative to the trajectory
     isabsval = 1
     do i=icolon+1,inumber-1
       ch = name(i:i)
       if ( (ch.eq.'+').or.(ch.eq.'-') ) isabsval = 0
     enddo
 
     ! Get the unit/shift axis
     shift_dir = name(inumber:80)
     if ( isabsval.eq.1 ) then
       shift_dir=trim(shift_dir)//'(ABS)'
     endif
 
  else
 
     shift_dir = 'nil'
     shift_val = 0.
 
  endif
 
 90 continue
 
    return
 
  ! Error handling
 100   continue
 
  print*,' ERROR: cannot split variable name ',trim(tvar)
  stop
 
end subroutine splitvar

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

Subversion Repositories lagranto.ecmwf

(root)/trunk/trace/trace.f90 – Rev 15 → 39