WebSVN – lagranto.ecmwf – Diff – /trunk/caltra/caltra.f


      PROGRAM caltra
 
C     ********************************************************************
C     *                                                                  *
C     * Calculates trajectories                                          *
C     *                                                                  *
C     *	Heini Wernli	   first version:	April 1993               *
C     * Michael Sprenger   major upgrade:       2008-2009                *
C     *                                                                  *
C     ********************************************************************
 
      implicit none      
 
c     --------------------------------------------------------------------
c     Declaration of parameters
c     --------------------------------------------------------------------
 
c     Maximum number of levels for input files
      integer   nlevmax
      parameter	(nlevmax=200)
 
c     Maximum number of input files (dates, length of trajectories)
      integer   ndatmax
      parameter	(ndatmax=500)
 
c     Numerical epsilon (for float comparison)
      real      eps
      parameter (eps=0.001)
 
c     Distance in m between 2 lat circles 
      real	deltay
      parameter	(deltay=1.112E5)
      
c     Numerical method for the integration (0=iterative Euler, 1=Runge-Kutta)
      integer   imethod
      parameter (imethod=1)
 
c     Number of iterations for iterative Euler scheme
      integer   numit
      parameter (numit=3)
 
c     Input and output format for trajectories (see iotra.f)
      integer   inpmode
      integer   outmode
 
c     Filename prefix for primary and secondary files (typically 'P' and 'S')
      character*1 prefix
      parameter   (prefix='P')
      character*1 srefix
      parameter   (srefix='S')
 
c     Physical constants - needed to compute potential temperature
      real      rdcp,tzero
      data      rdcp,tzero /0.286,273.15/
 
c     --------------------------------------------------------------------
c     Declaration of variables
c     --------------------------------------------------------------------
 
c     Input parameters
      integer                                fbflag          ! Flag for forward/backward mode
      integer                                numdat          ! Number of input files
      character*11                           dat(ndatmax)    ! Dates of input files
      real                                   timeinc         ! Time increment between input files
      real                                   per             ! Periodicity (=0 if none)
      integer                                ntra            ! Number of trajectories
      character*80                           cdfname         ! Name of output files
      real                                   ts              ! Time step
      real                                   tst,ten         ! Shift of start and end time relative to first data file
      integer                                deltout         ! Output time interval (in minutes)
      integer                                jflag           ! Jump flag (if =1 ground-touching trajectories reenter atmosphere)
      real                                   wfactor         ! Factor for vertical velocity field
      character*80                           strname         ! File with start positions
      character*80                           timecheck       ! Either 'yes' or 'no'
      character*80                           isen            ! Isentropic trajectories ('yes' or 'no')
      integer                                thons           ! Isentropic mode: is TH availanle on S
      character*80                           modlev          ! 2D (model level) trajectories ('yes' or 'no')
 
c     Trajectories
      integer                                ncol            ! Number of columns for insput trajectories
      real,allocatable, dimension (:,:,:) :: trainp          ! Input start coordinates (ntra,1,ncol)
      real,allocatable, dimension (:,:,:) :: traout          ! Output trajectories (ntra,ntim,4)
      integer                                reftime(6)      ! Reference date
      character*80                           vars(200)       ! Field names
      real,allocatable, dimension (:)     :: xx0,yy0,pp0     ! Position of air parcels
      integer,allocatable, dimension (:)  :: leftflag        ! Flag for domain-leaving
      real                                   xx1,yy1,pp1     ! Updated position of air parcel
      integer                                leftcount       ! Number of domain leaving trajectories
      integer                                ntim            ! Number of output time steps
      real,allocatable, dimension (:)     :: theta           ! Potential temperature for isentropic trajectories
      real,allocatable, dimension (:)     :: zindex          ! Vertical index for model-level (2D) trajectories
 
c     Meteorological fields
      real,allocatable, dimension (:)     :: spt0,spt1       ! Surface pressure
      real,allocatable, dimension (:)     :: uut0,uut1       ! Zonal wind
      real,allocatable, dimension (:)     :: vvt0,vvt1       ! Meridional wind
      real,allocatable, dimension (:)     :: wwt0,wwt1       ! Vertical wind
      real,allocatable, dimension (:)     :: p3t0,p3t1       ! 3d-pressure 
      real,allocatable, dimension (:)     :: tht0,tht1       ! 3d potential temperature
      real,allocatable, dimension (:)     :: sth0,sth1       ! Surface potential temperature
 
c     Grid description
      real                                   pollon,pollat   ! Longitude/latitude of pole
      real                                   ak(nlevmax)     ! Vertical layers and levels
      real                                   bk(nlevmax) 
      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
      real                                   mdv             ! Missing data value
 
c     Auxiliary variables                 
      real                                   delta,rd
      integer	                             itm,iloop,i,j,k,filo,lalo
      integer                                ierr,stat
      integer                                cdfid,fid
      real	                                 tstart,time0,time1,time
      real                                   reltpos0,reltpos1
      real                                   xind,yind,pind,pp,sp,stagz
      character*80                           filename,varname
      integer                                reftmp(6)
      character                              ch
      real                                   frac,tload
      integer                                itim
      integer                                wstep
      real                                   x1,y1,p1
      real                                   thetamin,thetamax
      real                                   zindexmin,zindexmax
 
c     Externals
      real                                   int_index4
      external                               int_index4
 
c     --------------------------------------------------------------------
c     Start of program, Read parameters
c     --------------------------------------------------------------------
 
c     Write start message
      print*,'========================================================='
      print*,'              *** START OF PROGRAM CALTRA ***'
      print*
 
c     Open the parameter file
      open(9,file='caltra.param')
 
c     Read flag for forward/backward mode (fbflag)
      read(9,*) fbflag
 
c     Read number of input files (numdat)
      read(9,*) numdat
      if (numdat.gt.ndatmax) then
        print*,' ERROR: too many input files ',numdat,ndatmax
        goto 993
      endif
 
c     Read list of input dates (dat, sort depending on forward/backward mode)
      if (fbflag.eq.1) then
        do itm=1,numdat
          read(9,'(a11)') dat(itm)
        enddo
      else
        do itm=numdat,1,-1
          read(9,'(a11)') dat(itm)
        enddo
      endif
 
c     Read time increment between input files (timeinc)
      read(9,*) timeinc
 
C     Read if data domain is periodic and its periodicity
      read(9,*) per
 
c     Read the number of trajectories and name of position file
      read(9,*) strname
      read(9,*) ntra
      read(9,*) ncol 
      if (ntra.eq.0) goto 991
 
C     Read the name of the output trajectory file and set the constants file
      read(9,*) cdfname
 
C     Read the timestep for trajectory calculation (convert from minutes to hours)
      read(9,*) ts
      ts=ts/60.        
 
C     Read shift of start and end time relative to first data file
      read(9,*) tst
      read(9,*) ten
 
C     Read output time interval (in minutes)
      read(9,*) deltout
 
C     Read jumpflag (if =1 ground-touching trajectories reenter the atmosphere)
      read(9,*) jflag
 
C     Read factor for vertical velocity field
      read(9,*) wfactor
 
c     Read the reference time and the time range
      read(9,*) reftime(1)                  ! year
      read(9,*) reftime(2)                  ! month
      read(9,*) reftime(3)                  ! day
      read(9,*) reftime(4)                  ! hour
      read(9,*) reftime(5)                  ! min
      read(9,*) reftime(6)                  ! time range (in min)
 
c     Read flag for 'no time check'
      read(9,*) timecheck
 
c     Read flag for isentropic trajectories
      read(9,*) isen, thons
 
c     Read flag for model-level trajectories (2D mode)
      read(9,*) modlev
 
c     Close the input file
      close(9)
 
c     Calculate the number of output time steps
      ntim = abs(reftime(6)/deltout) + 1
 
c     Set the formats of the input and output files
      call mode_tra(inpmode,strname)
      call mode_tra(outmode,cdfname)
      if (outmode.eq.-1) outmode=1
 
c     Write some status information
      print*,'---- INPUT PARAMETERS -----------------------------------'
      print* 
      print*,'  Forward/Backward       : ',fbflag
      print*,'  #input files           : ',numdat
      print*,'  First/last input file  : ',trim(dat(1)),' ... ',
     >                                     trim(dat(numdat))
      print*,'  time increment         : ',timeinc
      print*,'  Output file            : ',trim(cdfname)
      print*,'  Time step (min)        : ',60.*ts
      write(*,'(a27,f7.2,f7.2)') '   Time shift (start,end) : ',tst,ten
      print*,'  Output time interval   : ',deltout
      print*,'  Jump flag              : ',jflag
      print*,'  Vertical wind (scale)  : ',wfactor
      print*,'  Trajectory pos file    : ',trim(strname)
      print*,'  # of trajectories      : ',ntra
      print*,'  # of output timesteps  : ',ntim
      if ( inpmode.eq.-1) then
         print*,'  Input format           : (lon,lat,p)-list'
      else
         print*,'  Input format           : ',inpmode
      endif
      print*,'  Output format          : ',outmode
      print*,'  Periodicity            : ',per
      print*,'  Time check             : ',trim(timecheck)
      print*,'  Isentropic trajectories: ',trim(isen),thons
      print*,'  Model-level trajs (2D) : ',trim(modlev)
      print*
 
      if ( (isen.eq.'yes').and.(modlev.eq.'yes') ) then
         print*,
     >    '  WARNING: isentropic and 2D mode chosen -> 2D accepted'
         print*
         isen = 'no'
      endif
 
c     Init missing data value
      mdv = -999.
 
      print*,'---- FIXED NUMERICAL PARAMETERS -------------------------'
      print*
      print*,'  Numerical scheme       : ',imethod
      print*,'  Number of iterations   : ',numit
      print*,'  Filename prefix        : ',prefix
      print*,'  Missing data value     : ',mdv
      print*
 
c     --------------------------------------------------------------------
c     Read grid parameters, checks and allocate memory
c     --------------------------------------------------------------------
 
c     Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)
c     The negative <-fid> of the file identifier is used as a flag for parameter retrieval  
      filename = prefix//dat(1)
      varname  = 'U'
      nx       = 1
      ny       = 1
      nz       = 1
      tload    = -tst
      call input_open (fid,filename)
      call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >                 tload,pollon,pollat,rd,rd,nz,rd,rd,rd,timecheck)
      call input_close(fid)
 
C     Check if the number of levels is too large
      if (nz.gt.nlevmax) goto 993
 
C     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
            goto 992
        else if (abs(delta).lt.eps) then              ! Periodic and hemispheric
           hem=1
           per=360.
        endif
      else                                            ! Periodic and hemispheric
         hem=1
      endif
      
C     Allocate memory for some meteorological arrays
      allocate(spt0(nx*ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array spt0 ***'   ! Surface pressure
      allocate(spt1(nx*ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array spt1 ***'
      allocate(uut0(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array uut0 ***'   ! Zonal wind
      allocate(uut1(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array uut1 ***'
      allocate(vvt0(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array vvt0 ***'   ! Meridional wind
      allocate(vvt1(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array vvt1 ***'
      allocate(wwt0(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array wwt0 ***'   ! Vertical wind
      allocate(wwt1(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array wwt1 ***'
      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(tht0(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tht0 ***'   ! Potential temperature
      allocate(tht1(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tht1 ***'
      allocate(sth0(nx*ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array spt0 ***'   ! Surface potential temperature
      allocate(sth1(nx*ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array sth1 ***'
 
C     Get memory for trajectory arrays
      allocate(trainp(ntra,1,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array trainp   ***' ! Input start coordinates
      allocate(traout(ntra,ntim,4),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array traout   ***' ! Output trajectories
      allocate(xx0(ntra),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array xx0      ***' ! X position (longitude)
      allocate(yy0(ntra),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array yy0      ***' ! Y position (latitude)
      allocate(pp0(ntra),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array pp0      ***' ! Pressure
      allocate(leftflag(ntra),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array leftflag ***' ! Leaving-domain flag
      allocate(theta(ntra),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array theta ***'    ! Potential temperature for isentropic trajectories
      allocate(zindex(ntra),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array kind ***'     ! Vertical index for model-level trajectories
 
c     Write some status information
      print*,'---- CONSTANT GRID PARAMETERS ---------------------------'
      print*
      print*,'  xmin,xmax     : ',xmin,xmax
      print*,'  ymin,ymax     : ',ymin,ymax
      print*,'  dx,dy         : ',dx,dy
      print*,'  pollon,pollat : ',pollon,pollat
      print*,'  nx,ny,nz      : ',nx,ny,nz
      print*,'  per, hem      : ',per,hem
      print*
 
c     --------------------------------------------------------------------
c     Initialize the trajectory calculation
c     --------------------------------------------------------------------
 
c     Read start coordinates from file - Format (lon,lat,lev)
      if (inpmode.eq.-1) then
         open(fid,file=strname)
          do i=1,ntra
             read(fid,*) xx0(i),yy0(i),pp0(i)
          enddo
         close(fid)
 
c     Read start coordinates from trajectory file - check consistency of ref time
      else
         call ropen_tra(cdfid,strname,ntra,1,ncol,reftmp,vars,inpmode)
         call read_tra (cdfid,trainp,ntra,1,ncol,inpmode)
         do i=1,ntra
            time   = trainp(i,1,1)
            xx0(i) = trainp(i,1,2) 
            yy0(i) = trainp(i,1,3) 
            pp0(i) = trainp(i,1,4) 
         enddo
         call close_tra(cdfid,inpmode)
 
         if ( ( reftime(1).ne.reftmp(1) ).or.
     >        ( reftime(2).ne.reftmp(2) ).or.
     >        ( reftime(3).ne.reftmp(3) ).or.
     >        ( reftime(4).ne.reftmp(4) ).or.
     >        ( reftime(5).ne.reftmp(5) ) )
     >   then
            print*,' WARNING: Inconsistent reference times'
            write(*,'(5i8)') (reftime(i),i=1,5)
            write(*,'(5i8)') (reftmp (i),i=1,5)
            print*,'Enter a key to proceed...'
            stop
         endif
      endif
 
c     Set sign of time range
      reftime(6) = fbflag * reftime(6)
         
c     Write some status information
      print*,'---- REFERENCE DATE---------- ---------------------------'
      print*
      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             :',reftime(6),' min'
      print*
 
C     Save starting positions 
      itim = 1
      do i=1,ntra
         traout(i,itim,1) = 0.
         traout(i,itim,2) = xx0(i)
         traout(i,itim,3) = yy0(i)
         traout(i,itim,4) = pp0(i)
      enddo
 
c     Init the flag and the counter for trajectories leaving the domain
      leftcount=0
      do i=1,ntra
         leftflag(i)=0
      enddo
 
C     Convert time shifts <tst,ten> from <hh.mm> into fractional time
      call hhmm2frac(tst,frac)
      tst = frac
      call hhmm2frac(ten,frac)
      ten = frac
 
c     Get 3D and surface pressure from first data file 
      filename = prefix//dat(1)
      call input_open (fid,filename)
      if ( modlev.eq.'no' ) then 
         varname = 'P'
         call input_grid
     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
      else  
         varname = 'P.ML'
         call input_grid
     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
      endif
      call input_close(fid)
      
c     Check that all starting positions are above topography    
      do i=1,ntra
 
C       Interpolate surface pressure to actual position (from first input file)
        x1 = xx0(i)
        y1 = yy0(i)
        call get_index4 (xind,yind,pind,x1,y1,1050.,0.,
     >                   p3t1,p3t1,spt1,spt1,3,
     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
        sp = int_index4 (spt1,spt1,nx,ny,1,xind,yind,1.,0.,mdv)
 
c       Decide whether to keep the trajectory
        if ( pp0(i).gt.sp ) then
            write(*,'(a30,3f10.2)')
     >               'WARNING: starting point below topography ',
     >               xx0(i),yy0(i),pp0(i)
            leftflag(i) = 1
        endif
 
      enddo
 
c     Special handling for isentropic trajectories - read potential
c     temperature from S file or calculate it based on temperature and
c     pressure from P file; then, calculate for each trajectory its
c     potential temperature - which will stay fixed over time
      if ( isen.eq.'yes' ) then
 
c         Get potential temperature from S file
          if ( thons.eq.1 ) then
              filename = srefix//dat(1)
              print*,' TH <- ',trim(filename)
              call input_open (fid,filename)
              varname='TH'
              call input_wind
     >            (fid,varname,tht1,tload,stagz,mdv,
     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
              call input_close(fid)
 
c         Calculate potential temperature from P file
          else
              filename = prefix//dat(1)
              print*,' TH = T * (1000/P)^RDCP <- ',trim(filename)
              call input_open (fid,filename)
              varname='T'
              call input_wind
     >            (fid,varname,tht1,tload,stagz,mdv,
     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
              call input_close(fid)
              do i=1,nx*ny*nz
                if (tht1(i).lt.100.) then
                   tht1(i)=(tht1(i)+tzero)*( (1000./p3t1(i))**rdcp )
                else
                   tht1(i)=tht1(i)*( (1000./p3t1(i))**rdcp )
                endif
              enddo
          endif
 
c         Take surface potential temperature from lowest level
          do i=1,nx*ny
            sth1(i) = tht1(i)
          enddo
 
c         Calculate now the potential temperature of all trajectories
          do i=1,ntra
 
             x1 = xx0(i)
             y1 = yy0(i)
             p1 = pp0(i)
             call get_index4 (xind,yind,pind,x1,y1,p1,0.,
     >                 p3t1,p3t1,spt1,spt1,3,
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
             theta(i) =
     >        int_index4(tht1,tht1,nx,ny,nz,xind,yind,pind,0.,mdv)
 
          enddo
 
c         Write info about theta range of starting positions
          thetamin = theta(1)
          thetamax = theta(1)
          do i=2,ntra
            if ( theta(i).gt.thetamax ) thetamax = theta(i)
            if ( theta(i).lt.thetamin ) thetamin = theta(i)
          enddo
 
c         Write some status information
          print*
          print*,
     >     '---- THETA RANGE OF ISENTROPIC TRAJECTORIES -------------'
          print*
          print*,' Theta(min)             :',thetamin
          print*,' Theta(max)             :',thetamax
 
      endif
 
c     Special handling for model-level (2D) trajectories - get the
c     vertical index for each trajectory - which will remain fixed
      if ( modlev.eq.'yes' ) then
          do i=1,ntra
             x1 = xx0(i)
             y1 = yy0(i)
             p1 = pp0(i)
             call get_index4 (xind,yind,pind,x1,y1,p1,0.,
     >                 p3t1,p3t1,spt1,spt1,3,
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
             zindex(i) = pind
          enddo
 
          do i=1,nz
            print*,i,p3t1(189+(141-1)*nx+(i-1)*nx*ny)
          enddo
             print*,x1,y1,p1
             print*,xind,yind,pind
 
c         Write info about zindex range of starting positions
          zindexmin = zindex(1)
          zindexmax = zindex(1)
          do i=2,ntra
            if ( zindex(i).gt.zindexmax ) zindexmax = zindex(i)
            if ( zindex(i).lt.zindexmin ) zindexmin = zindex(i)
          enddo
 
c         Write some status information
          print*
          print*,
     >     '---- INDEX RANGE OF MODEL-LEVEL TRAJECTORIES ------------'
          print*
          print*,' Zindex(min)            :',zindexmin
          print*,' Zindex(max)            :',zindexmax
 
 
      endif
 
c     -----------------------------------------------------------------------
c     Loop to calculate trajectories
c     -----------------------------------------------------------------------
 
c     Write some status information
      print*
      print*,'---- TRAJECTORIES ----------- ---------------------------'
      print*    
 
C     Set the time for the first data file (depending on forward/backward mode)
      if (fbflag.eq.1) then
        tstart = -tst
      else
        tstart = tst
      endif
 
c     Set the minute counter for output
      wstep = 0
 
c     Read wind fields and vertical grid from first file
      filename = prefix//dat(1)
 
      call frac2hhmm(tstart,tload)
 
      write(*,'(a16,a20,f9.2)') '  (file,time) : ',
     >                       trim(filename),tload
 
      call input_open (fid,filename)
 
      varname='U'                                      ! U
      call input_wind 
     >    (fid,varname,uut1,tload,stagz,mdv,
     >     xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
 
      varname='V'                                      ! V
      call input_wind 
     >    (fid,varname,vvt1,tload,stagz,mdv,
     >     xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
 
      if ( (modlev.eq.'no').and.(isen.eq.'no') ) then
         varname='OMEGA'                                  ! OMEGA
         call input_wind
     >       (fid,varname,wwt1,tload,stagz,mdv,
     >        xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
      endif
 
      if ( modlev.eq.'no' ) then
         call input_grid                                  ! GRID - AK,BK -> P
     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
      else
         varname='P.ML'                                   ! GRID - P,PS
         call input_grid                                  !
     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
      endif
 
      call input_close(fid)
 
c     Special handling for isentropic trajectories - read potential
c     temperature from S file or calculate it based on temperature and
c     pressure from P file
      if ( isen.eq.'yes' ) then
 
c         Get potential temperature from S file
          if ( thons.eq.1 ) then
              filename = srefix//dat(1)
              print*,' TH <- ',trim(filename)
              call input_open (fid,filename)
              varname='TH'
              call input_wind
     >            (fid,varname,tht1,tload,stagz,mdv,
     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
              call input_close(fid)
 
c         Calculate potential temperature from P file
          else
              filename = prefix//dat(1)
              print*,' TH = T * (1000/P)^RDCP <- ',trim(filename)
              call input_open (fid,filename)
              varname='T'
              call input_wind
     >            (fid,varname,tht1,tload,stagz,mdv,
     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
              call input_close(fid)
              do i=1,nx*ny*nz
                if (tht1(i).lt.100.) then
                   tht1(i)=(tht1(i)+tzero)*( (1000./p3t1(i))**rdcp )
                else
                   tht1(i)=tht1(i)*( (1000./p3t1(i))**rdcp )
                endif
              enddo
          endif
 
c         Take surface potential temperature from lowest level
          do i=1,nx*ny
            sth1(i) = tht1(i)
          enddo
      endif
 
c     Loop over all input files (time step is <timeinc>)
      do itm=1,numdat-1
 
c       Calculate actual and next time
        time0 = tstart+real(itm-1)*timeinc*fbflag
        time1 = time0+timeinc*fbflag
 
c       Copy old velocities and pressure fields to new ones       
        do i=1,nx*ny*nz
           uut0(i)=uut1(i)
           vvt0(i)=vvt1(i)
           wwt0(i)=wwt1(i)
           p3t0(i)=p3t1(i)
           tht0(i)=tht1(i)
        enddo
        do i=1,nx*ny
           spt0(i)=spt1(i)
           sth0(i)=sth1(i)
        enddo
 
c       Read wind fields and surface pressure at next time
        filename = prefix//dat(itm+1)
 
        call frac2hhmm(time1,tload)
        write(*,'(a16,a20,f9.2)') '  (file,time) : ',
     >                          trim(filename),tload
 
        call input_open (fid,filename)
 
        varname='U'                                     ! U
        call input_wind 
     >       (fid,varname,uut1,tload,stagz,mdv,
     >        xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
 
        varname='V'                                     ! V
        call input_wind 
     >       (fid,varname,vvt1,tload,stagz,mdv,
     >        xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
 
        if ( (modlev.eq.'no').and.(isen.eq.'no') ) then
           varname='OMEGA'                              ! OMEGA
           call input_wind
     >          (fid,varname,wwt1,tload,stagz,mdv,
     >           xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
        endif
 
        if ( modlev.eq.'no' ) then
           call input_grid                                  ! GRID - AK,NK -> P
     >          (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >           tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
        else
           varname='P.ML'                                   ! GRID - P,PS
           call input_grid                                  !
     >       (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >        tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)
        endif
 
        call input_close(fid)
 
c     Special handling for isentropic trajectories - read potential
c     temperature from S file or calculate it based on temperature and
c     pressure from P file
      if ( isen.eq.'yes' ) then
 
c         Get TH from S file
          if ( thons.eq.1 ) then
              filename = srefix//dat(itm+1)
              print*,' TH <- ',trim(filename)
              call input_open (fid,filename)
              varname='TH'
              call input_wind
     >            (fid,varname,tht1,tload,stagz,mdv,
     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
              call input_close(fid)
 
c         Calculate potential temperature from P file
          else
              filename = prefix//dat(itm+1)
              print*,' TH = T * (1000/P)^RDCP <- ',trim(filename)
              call input_open (fid,filename)
              varname='T'
              call input_wind
     >            (fid,varname,tht1,tload,stagz,mdv,
     >              xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
              call input_close(fid)
              do i=1,nx*ny*nz
                if (tht1(i).lt.100.) then
                   tht1(i)=(tht1(i)+tzero)*( (1000./p3t1(i))**rdcp )
                else
                   tht1(i)=tht1(i)*( (1000./p3t1(i))**rdcp )
                endif
              enddo
          endif
 
c         Take surface potential temperature from lowest level
          do i=1,nx*ny
            sth1(i) = tht1(i)
          enddo
      endif
        
C       Determine the first and last loop indices
        if (numdat.eq.2) then
          filo = nint(tst/ts)+1
          lalo = nint((timeinc-ten)/ts) 
        elseif ( itm.eq.1 ) then
          filo = nint(tst/ts)+1
          lalo = nint(timeinc/ts)
        else if (itm.eq.numdat-1) then
          filo = 1
          lalo = nint((timeinc-ten)/ts)
        else
          filo = 1
          lalo = nint(timeinc/ts)
        endif
 
c       Split the interval <timeinc> into computational time steps <ts>
        do iloop=filo,lalo
 
C         Calculate relative time position in the interval timeinc (0=beginning, 1=end)
          reltpos0 = ((real(iloop)-1.)*ts)/timeinc
          reltpos1 = real(iloop)*ts/timeinc
 
c         Timestep for all trajectories
          do i=1,ntra
 
C           Check if trajectory has already left the data domain
            if (leftflag(i).ne.1) then	
 
c             3D: Iterative Euler timestep (x0,y0,p0 -> x1,y1,p1)
              if ( (imethod.eq.1   ).and.
     >             (isen   .eq.'no').and.
     >             (modlev .eq.'no') )
     >        then
                 call euler_3d(
     >               xx1,yy1,pp1,leftflag(i),
     >               xx0(i),yy0(i),pp0(i),reltpos0,reltpos1,
     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz)
 
c             3D: Runge-Kutta timestep (x0,y0,p0 -> x1,y1,p1)
              else if ( (imethod.eq.2   ).and.
     >                  (isen   .eq.'no').and.
     >                  (modlev .eq.'no') )
     >        then
                 call runge(
     >               xx1,yy1,pp1,leftflag(i),
     >               xx0(i),yy0(i),pp0(i),reltpos0,reltpos1,
     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz)
 
c             ISENTROPIC: Iterative Euler timestep (x0,y0,p0 -> x1,y1,p1)
              else if ( (imethod.eq.1    ).and.
     >                  (isen   .eq.'yes').and.
     >                  (modlev .eq.'no' ) )
     >        then
                 call euler_isen(
     >               xx1,yy1,pp1,leftflag(i),
     >               xx0(i),yy0(i),pp0(i),theta(i),reltpos0,reltpos1,
     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,
     >               sth0,sth1,tht0,tht1,
     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz)
 
c             MODEL-LEVEL (2D): Iterative Euler timestep (x0,y0,p0 -> x1,y1,p1)
              else if ( (imethod.eq.1    ).and.
     >                  (isen   .eq.'no' ).and.
     >                  (modlev .eq.'yes') )
     >        then
                 call euler_2d(
     >               xx1,yy1,pp1,leftflag(i),
     >               xx0(i),yy0(i),pp0(i),zindex(i),reltpos0,reltpos1,
     >               ts*3600,numit,jflag,mdv,wfactor,fbflag,
     >               spt0,spt1,p3t0,p3t1,uut0,uut1,vvt0,vvt1,
     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz)
 
              endif
 
c             Update trajectory position, or increase number of trajectories leaving domain
              if (leftflag(i).eq.1) then
                leftcount=leftcount+1
                if ( leftcount.lt.10 ) then
                   print*,'     -> Trajectory ',i,' leaves domain'
                elseif ( leftcount.eq.10 ) then
                   print*,'     -> N>=10 trajectories leave domain'
                endif
              else
                xx0(i)=xx1      
                yy0(i)=yy1
                pp0(i)=pp1
              endif
 
c          Trajectory has already left data domain (mark as <mdv>)
           else
              xx0(i)=mdv
              yy0(i)=mdv
              pp0(i)=mdv
              
           endif
 
          enddo
 
C         Save positions only every deltout minutes
          delta = aint(iloop*60*ts/deltout)-iloop*60*ts/deltout
          if (abs(delta).lt.eps) then
c          wstep = wstep + abs(ts)
c          if ( mod(wstep,deltout).eq.0 ) then
            time = time0+reltpos1*timeinc*fbflag
            itim = itim + 1
            if ( itim.le.ntim ) then
              do i=1,ntra
                 call frac2hhmm(time,tload)
                 traout(i,itim,1) = tload
                 traout(i,itim,2) = xx0(i)
                 traout(i,itim,3) = yy0(i)
                 traout(i,itim,4) = pp0(i)
              enddo
            endif
          endif
 
        enddo
 
      enddo
 
c     Write trajectory file
      vars(1)  ='time'
      vars(2)  ='lon'
      vars(3)  ='lat'
      vars(4)  ='p'
      call wopen_tra(cdfid,cdfname,ntra,ntim,4,reftime,vars,outmode)
      call write_tra(cdfid,traout,ntra,ntim,4,outmode)
      call close_tra(cdfid,outmode)   
 
c     Write some status information, and end of program message
      print*  
      print*,'---- STATUS INFORMATION --------------------------------'
      print*
      print*,'  #leaving domain    ', leftcount
      print*,'  #staying in domain ', ntra-leftcount
      print*
      print*,'              *** END OF PROGRAM CALTRA ***'
      print*,'========================================================='
 
      stop
 
c     ------------------------------------------------------------------
c     Exception handling
c     ------------------------------------------------------------------
 
 991  write(*,*) '*** ERROR: all start points outside the data domain'
      call exit(1)
      
 992  write(*,*) '*** ERROR: close arrays on files (prog. closear)'
      call exit(1)
 
 993  write(*,*) '*** ERROR: problems with array size'
      call exit(1)
 
      end 
 
 
c     *******************************************************************
c     * Time step : either Euler or Runge-Kutta                         *
c     *******************************************************************
 
C     Time-step from (x0,y0,p0) to (x1,y1,p1)
C
C     (x0,y0,p0) input	coordinates (long,lat,p) for starting point
C     (x1,y1,p1) output	coordinates (long,lat,p) for end point
C     deltat	 input	timestep in seconds
C     numit	 input	number of iterations
C     jump	 input  flag (=1 trajectories don't enter the ground)
C     left	 output	flag (=1 if trajectory leaves data domain)
 
c     -------------------------------------------------------------------
c     Iterative Euler time step (KINEMATIC 3D TRAJECTORIES)
c     -------------------------------------------------------------------
 
      subroutine euler_3d(x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
     >                deltat,numit,jump,mdv,wfactor,fbflag,
     >		          spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
     >                xmin,ymin,dx,dy,per,hem,nx,ny,nz)
 
      implicit none
 
c     Declaration of subroutine parameters
      integer      nx,ny,nz
      real         x1,y1,p1
      integer      left
      real	   x0,y0,p0
      real         reltpos0,reltpos1
      real   	   deltat
      integer      numit
      integer      jump
      real         wfactor
      integer      fbflag
      real     	   spt0(nx*ny)   ,spt1(nx*ny)
      real         uut0(nx*ny*nz),uut1(nx*ny*nz)
      real 	   vvt0(nx*ny*nz),vvt1(nx*ny*nz)
      real         wwt0(nx*ny*nz),wwt1(nx*ny*nz)
      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
      real         xmin,ymin,dx,dy
      real         per
      integer      hem
      real         mdv
 
c     Numerical and physical constants
      real         deltay
      parameter    (deltay=1.112E5)  ! Distance in m between 2 lat circles
      real         pi                       
      parameter    (pi=3.1415927)    ! Pi
 
c     Auxiliary variables
      real         xmax,ymax
      real	   xind,yind,pind
      real	   u0,v0,w0,u1,v1,w1,u,v,w,sp
      integer	   icount
      character    ch
 
c     Externals    
      real         int_index4
      external     int_index4
 
c     Reset the flag for domain-leaving
      left=0
 
c     Set the esat-north bounray of the domain
      xmax = xmin+real(nx-1)*dx
      ymax = ymin+real(ny-1)*dy
 
C     Interpolate wind fields to starting position (x0,y0,p0)
      call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0,
     >                 p3d0,p3d1,spt0,spt1,3,
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
      u0 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
      v0 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
      w0 = int_index4(wwt0,wwt1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
 
c     Force the near-surface wind to zero
      if (pind.lt.1.) w0=w0*pind
 
C     For first iteration take ending position equal to starting position
      x1=x0
      y1=y0
      p1=p0
 
C     Iterative calculation of new position
      do icount=1,numit
 
C        Calculate new winds for advection
         call get_index4 (xind,yind,pind,x1,y1,p1,reltpos1,
     >                    p3d0,p3d1,spt0,spt1,3,
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
         u1 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
         v1 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
         w1 = int_index4(wwt0,wwt1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
 
c        Force the near-surface wind to zero
         if (pind.lt.1.) w1=w1*pind
 
c        Get the new velocity in between
         u=(u0+u1)/2.
         v=(v0+v1)/2.
         w=(w0+w1)/2.
         
C        Calculate new positions
         x1 = x0 + fbflag*u*deltat/(deltay*cos(y0*pi/180.))
         y1 = y0 + fbflag*v*deltat/deltay
         p1 = p0 + fbflag*wfactor*w*deltat/100.
 
c       Handle pole problems (crossing and near pole trajectory)
        if ((hem.eq.1).and.(y1.gt.90.)) then
          y1=180.-y1
          x1=x1+per/2.
        endif
        if ((hem.eq.1).and.(y1.lt.-90.)) then
          y1=-180.-y1
          x1=x1+per/2.
        endif
        if (y1.gt.89.99) then
           y1=89.99
        endif
 
c       Handle crossings of the dateline
        if ((hem.eq.1).and.(x1.gt.xmin+per-dx)) then
           x1=xmin+amod(x1-xmin,per)
        endif
        if ((hem.eq.1).and.(x1.lt.xmin)) then
           x1=xmin+per+amod(x1-xmin,per)
        endif
 
C       Interpolate surface pressure to actual position
        call get_index4 (xind,yind,pind,x1,y1,1050.,reltpos1,
     >                   p3d0,p3d1,spt0,spt1,3,
     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
        sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos1,mdv)
 
c       Handle trajectories which cross the lower boundary (jump flag)
        if ((jump.eq.1).and.(p1.gt.sp)) p1=sp-10.
 
C       Check if trajectory leaves data domain
        if ( ( (hem.eq.0).and.(x1.lt.xmin)    ).or.
     >       ( (hem.eq.0).and.(x1.gt.xmax-dx) ).or.
     >         (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.gt.sp) )
     >  then
          left=1
          goto 100
        endif
 
      enddo
 
c     Exit point for subroutine
 100  continue
 
      return
 
      end
 
c     -------------------------------------------------------------------
c     Iterative Euler time step (ISENTROPIC)
c     -------------------------------------------------------------------
 
      subroutine euler_isen(x1,y1,p1,left,x0,y0,p0,theta,
     >                reltpos0,reltpos1,
     >                deltat,numit,jump,mdv,wfactor,fbflag,
     >                spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,
     >                sth0,sth1,tht0,tht1,
     >                xmin,ymin,dx,dy,per,hem,nx,ny,nz)
 
      implicit none
 
c     Declaration of subroutine parameters
      integer      nx,ny,nz
      real         x1,y1,p1
      integer      left
      real         x0,y0,p0
      real         reltpos0,reltpos1
      real         deltat
      integer      numit
      integer      jump
      real         wfactor
      integer      fbflag
      real         spt0(nx*ny)   ,spt1(nx*ny)
      real         sth0(nx*ny)   ,sth1(nx*ny)
      real         uut0(nx*ny*nz),uut1(nx*ny*nz)
      real         vvt0(nx*ny*nz),vvt1(nx*ny*nz)
      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
      real         tht0(nx*ny*nz),tht1(nx*ny*nz)
      real         xmin,ymin,dx,dy
      real         per
      integer      hem
      real         mdv
      real         theta
 
c     Numerical and physical constants
      real         deltay
      parameter    (deltay=1.112E5)  ! Distance in m between 2 lat circles
      real         pi
      parameter    (pi=3.1415927)    ! Pi
 
c     Auxiliary variables
      real         xmax,ymax
      real         xind,yind,pind
      real         u0,v0,w0,u1,v1,w1,u,v,w,sp
      integer      icount
      character    ch
 
c     Externals
      real         int_index4
      external     int_index4
 
c     Reset the flag for domain-leaving
      left=0
 
c     Set the esat-north bounray of the domain
      xmax = xmin+real(nx-1)*dx
      ymax = ymin+real(ny-1)*dy
 
C     Interpolate wind fields to starting position (x0,y0,p0)
      call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0,
     >                 p3d0,p3d1,spt0,spt1,3,
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
      u0 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
      v0 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
 
C     For first iteration take ending position equal to starting position
      x1=x0
      y1=y0
      p1=p0
 
C     Iterative calculation of new position
      do icount=1,numit
 
C        Calculate new winds for advection
         call get_index4 (xind,yind,pind,x1,y1,p1,reltpos1,
     >                    p3d0,p3d1,spt0,spt1,3,
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
         u1 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
         v1 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
 
c        Get the new velocity in between
         u=(u0+u1)/2.
         v=(v0+v1)/2.
 
C        Calculate new positions
         x1 = x0 + fbflag*u*deltat/(deltay*cos(y0*pi/180.))
         y1 = y0 + fbflag*v*deltat/deltay
 
c        Get the pressure on the isentropic surface at the new position
         call get_index4 (xind,yind,pind,x1,y1,theta,reltpos1,
     >                    tht0,tht1,sth0,sth1,1,
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
         p1 = int_index4(p3d0,p3d1,nx,ny,nz,xind,yind,pind,reltpos1,mdv)
 
c       Handle pole problems (crossing and near pole trajectory)
        if ((hem.eq.1).and.(y1.gt.90.)) then
          y1=180.-y1
          x1=x1+per/2.
        endif
        if ((hem.eq.1).and.(y1.lt.-90.)) then
          y1=-180.-y1
          x1=x1+per/2.
        endif
        if (y1.gt.89.99) then
           y1=89.99
        endif
 
c       Handle crossings of the dateline
        if ((hem.eq.1).and.(x1.gt.xmin+per-dx)) then
           x1=xmin+amod(x1-xmin,per)
        endif
        if ((hem.eq.1).and.(x1.lt.xmin)) then
           x1=xmin+per+amod(x1-xmin,per)
        endif
 
C       Interpolate surface pressure to actual position
        call get_index4 (xind,yind,pind,x1,y1,1050.,reltpos1,
     >                   p3d0,p3d1,spt0,spt1,3,
     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
        sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos1,mdv)
 
c       Handle trajectories which cross the lower boundary (jump flag)
        if ((jump.eq.1).and.(p1.gt.sp)) p1=sp-10.
 
C       Check if trajectory leaves data domain
        if ( ( (hem.eq.0).and.(x1.lt.xmin)    ).or.
     >       ( (hem.eq.0).and.(x1.gt.xmax-dx) ).or.
     >         (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.gt.sp) )
     >  then
          left=1
          goto 100
        endif
 
      enddo
 
c     Exit point for subroutine
 100  continue
 
      return
 
      end
 
c     -------------------------------------------------------------------
c     Iterative Euler time step (MODEL-LEVEL, 2D)
c     -------------------------------------------------------------------
 
      subroutine euler_2d(x1,y1,p1,left,x0,y0,p0,zindex,
     >                reltpos0,reltpos1,
     >                deltat,numit,jump,mdv,wfactor,fbflag,
     >                spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,
     >                xmin,ymin,dx,dy,per,hem,nx,ny,nz)
 
      implicit none
 
c     Declaration of subroutine parameters
      integer      nx,ny,nz
      real         x1,y1,p1
      integer      left
      real         x0,y0,p0
      real         reltpos0,reltpos1
      real         deltat
      integer      numit
      integer      jump
      real         wfactor
      integer      fbflag
      real         spt0(nx*ny)   ,spt1(nx*ny)
      real         uut0(nx*ny*nz),uut1(nx*ny*nz)
      real         vvt0(nx*ny*nz),vvt1(nx*ny*nz)
      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
      real         xmin,ymin,dx,dy
      real         per
      integer      hem
      real         mdv
      real         zindex
 
c     Numerical and physical constants
      real         deltay
      parameter    (deltay=1.112E5)  ! Distance in m between 2 lat circles
      real         pi
      parameter    (pi=3.1415927)    ! Pi
      real         eps
      parameter    (eps=0.001)
 
c     Auxiliary variables
      real         xmax,ymax
      real         xind,yind,pind
      real         u0,v0,w0,u1,v1,w1,u,v,w,sp
      integer      icount
      character    ch
 
c     Externals
      real         int_index4
      external     int_index4
 
c     Reset the flag for domain-leaving
      left=0
 
c     Set the esat-north bounray of the domain
      xmax = xmin+real(nx-1)*dx
      ymax = ymin+real(ny-1)*dy
 
C     Interpolate wind fields to starting position (x0,y0,p0)
      call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0,
     >                 p3d0,p3d1,spt0,spt1,3,
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
      u0 = int_index4(uut0,uut1,nx,ny,nz,xind,yind,zindex,reltpos0,mdv)
      v0 = int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,zindex,reltpos0,mdv)
 
C     For first iteration take ending position equal to starting position
      x1=x0
      y1=y0
      p1=p0
 
C     Iterative calculation of new position
      do icount=1,numit
 
C        Calculate new winds for advection
         call get_index4 (xind,yind,pind,x1,y1,p1,reltpos1,
     >                    p3d0,p3d1,spt0,spt1,3,
     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
         u1=int_index4(uut0,uut1,nx,ny,nz,xind,yind,zindex,reltpos1,mdv)
         v1=int_index4(vvt0,vvt1,nx,ny,nz,xind,yind,zindex,reltpos1,mdv)
         if ( abs(u1-mdv).lt.eps ) then
             left = 1
             goto 100
         endif
         if ( abs(v1-mdv).lt.eps ) then
             left = 1
             goto 100
         endif
 
c        Get the new velocity in between
         u=(u0+u1)/2.
         v=(v0+v1)/2.
 
C        Calculate new positions
         x1 = x0 + fbflag*u*deltat/(deltay*cos(y0*pi/180.))
         y1 = y0 + fbflag*v*deltat/deltay
 
c        Get the pressure on the model surface at the new position
         xind = (x1 - xmin ) / dx + 1.
         yind = (y1 - ymin ) / dy + 1.
         p1 =
     >     int_index4(p3d0,p3d1,nx,ny,nz,xind,yind,zindex,reltpos1,mdv)
         if ( abs(p1-mdv).lt.eps ) then
             left = 1
             goto 100
         endif
 
c       Handle pole problems (crossing and near pole trajectory)
        if ((hem.eq.1).and.(y1.gt.90.)) then
          y1=180.-y1
          x1=x1+per/2.
        endif
        if ((hem.eq.1).and.(y1.lt.-90.)) then
          y1=-180.-y1
          x1=x1+per/2.
        endif
        if (y1.gt.89.99) then
           y1=89.99
        endif
 
c       Handle crossings of the dateline
        if ((hem.eq.1).and.(x1.gt.xmin+per-dx)) then
           x1=xmin+amod(x1-xmin,per)
        endif
        if ((hem.eq.1).and.(x1.lt.xmin)) then
           x1=xmin+per+amod(x1-xmin,per)
        endif
 
C       Interpolate surface pressure to actual position
        call get_index4 (xind,yind,pind,x1,y1,1050.,reltpos1,
     >                   p3d0,p3d1,spt0,spt1,3,
     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
        sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos1,mdv)
 
c       Handle trajectories which cross the lower boundary (jump flag)
        if ((jump.eq.1).and.(p1.gt.sp)) p1=sp-10.
 
C       Check if trajectory leaves data domain
        if ( ( (hem.eq.0).and.(x1.lt.xmin)    ).or.
     >       ( (hem.eq.0).and.(x1.gt.xmax-dx) ).or.
     >         (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.gt.sp) )
     >  then
          left=1
          goto 100
        endif
 
      enddo
 
c     Exit point for subroutine
 100  continue
 
      return
 
      end
 
c     -------------------------------------------------------------------
c     Runge-Kutta (4th order) time-step
c     -------------------------------------------------------------------
 
      subroutine runge(x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
     >                 deltat,numit,jump,mdv,wfactor,fbflag,
     >		       spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
     >                 xmin,ymin,dx,dy,per,hem,nx,ny,nz)
 
      implicit none
 
c     Declaration of subroutine parameters
      integer      nx,ny,nz
      real         x1,y1,p1
      integer      left
      real	   x0,y0,p0
      real         reltpos0,reltpos1
      real   	   deltat
      integer      numit
      integer      jump
      real         wfactor
      integer      fbflag
      real     	   spt0(nx*ny)   ,spt1(nx*ny)
      real         uut0(nx*ny*nz),uut1(nx*ny*nz)
      real 	   vvt0(nx*ny*nz),vvt1(nx*ny*nz)
      real         wwt0(nx*ny*nz),wwt1(nx*ny*nz)
      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
      real         xmin,ymin,dx,dy
      real         per
      integer      hem
      real         mdv
 
c     Numerical and physical constants
      real         deltay
      parameter    (deltay=1.112E5)  ! Distance in m between 2 lat circles
      real         pi                       
      parameter    (pi=3.1415927)    ! Pi
 
c     Auxiliary variables
      real         xmax,ymax
      real	   xind,yind,pind
      real	   u0,v0,w0,u1,v1,w1,u,v,w,sp
      integer	   icount,n
      real         xs,ys,ps,xk(4),yk(4),pk(4)
      real         reltpos
 
c     Externals    
      real         int_index4
      external     int_index4
 
c     Reset the flag for domain-leaving
      left=0
 
c     Set the esat-north bounray of the domain
      xmax = xmin+real(nx-1)*dx
      ymax = ymin+real(ny-1)*dy
 
c     Apply the Runge Kutta scheme
      do n=1,4
 
c       Get intermediate position and relative time
        if (n.eq.1) then
          xs=0.
          ys=0.
          ps=0.
          reltpos=reltpos0
        else if (n.eq.4) then
          xs=xk(3)
          ys=yk(3)
          ps=pk(3)
          reltpos=reltpos1
        else
          xs=xk(n-1)/2.
          ys=yk(n-1)/2.
          ps=pk(n-1)/2.
          reltpos=(reltpos0+reltpos1)/2.
        endif
        
C       Calculate new winds for advection
        call get_index4 (xind,yind,pind,x0+xs,y0+ys,p0+ps,reltpos,
     >                   p3d0,p3d1,spt0,spt1,3,
     >                   nx,ny,nz,xmin,ymin,dx,dy,mdv)
        u = int_index4 (uut0,uut1,nx,ny,nz,xind,yind,pind,reltpos,mdv)
        v = int_index4 (vvt0,vvt1,nx,ny,nz,xind,yind,pind,reltpos,mdv)
        w = int_index4 (wwt0,wwt1,nx,ny,nz,xind,yind,pind,reltpos,mdv)
         
c       Force the near-surface wind to zero
        if (pind.lt.1.) w1=w1*pind
 
c       Update position and keep them
        xk(n)=fbflag*u*deltat/(deltay*cos(y0*pi/180.))
        yk(n)=fbflag*v*deltat/deltay
        pk(n)=fbflag*w*deltat*wfactor/100.
 
      enddo
 
C     Calculate new positions
      x1=x0+(1./6.)*(xk(1)+2.*xk(2)+2.*xk(3)+xk(4))
      y1=y0+(1./6.)*(yk(1)+2.*yk(2)+2.*yk(3)+yk(4))
      p1=p0+(1./6.)*(pk(1)+2.*pk(2)+2.*pk(3)+pk(4))
 
c     Handle pole problems (crossing and near pole trajectory)
      if ((hem.eq.1).and.(y1.gt.90.)) then
         y1=180.-y1
         x1=x1+per/2.
      endif
      if ((hem.eq.1).and.(y1.lt.-90.)) then
         y1=-180.-y1
         x1=x1+per/2.
      endif
      if (y1.gt.89.99) then
         y1=89.99
      endif
      
c     Handle crossings of the dateline
      if ((hem.eq.1).and.(x1.gt.xmin+per-dx)) then
         x1=xmin+amod(x1-xmin,per)
      endif
      if ((hem.eq.1).and.(x1.lt.xmin)) then
         x1=xmin+per+amod(x1-xmin,per)
      endif
      
C     Interpolate surface pressure to actual position
      call get_index4 (xind,yind,pind,x1,y1,1050.,reltpos1,
     >                 p3d0,p3d1,spt0,spt1,3,
     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
      sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1,reltpos,mdv)
 
c     Handle trajectories which cross the lower boundary (jump flag)
      if ((jump.eq.1).and.(p1.gt.sp)) p1=sp-10.
      
C     Check if trajectory leaves data domain
      if ( ( (hem.eq.0).and.(x1.lt.xmin)    ).or.
     >     ( (hem.eq.0).and.(x1.gt.xmax-dx) ).or.
     >       (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.gt.sp) )
     >then
         left=1
         goto 100
      endif
      
c     Exit point fdor subroutine
 100  continue
 
      return
      end

Rev 19	Rev 27
1	`PROGRAM caltra`	1	`PROGRAM caltra`
2		2
3	`C ********************************************************************`	3	`C ********************************************************************`
4	`C * *`	4	`C * *`
5	`C * Calculates trajectories *`	5	`C * Calculates trajectories *`
6	`C * *`	6	`C * *`
7	`C * Heini Wernli first version: April 1993 *`	7	`C * Heini Wernli first version: April 1993 *`
8	`C * Michael Sprenger major upgrade: 2008-2009 *`	8	`C * Michael Sprenger major upgrade: 2008-2009 *`
9	`C * *`	9	`C * *`
10	`C ********************************************************************`	10	`C ********************************************************************`
11		11
12	`implicit none`	12	`implicit none`
13		13
14	`c --------------------------------------------------------------------`	14	`c --------------------------------------------------------------------`
15	`c Declaration of parameters`	15	`c Declaration of parameters`
16	`c --------------------------------------------------------------------`	16	`c --------------------------------------------------------------------`
17		17
18	`c Maximum number of levels for input files`	18	`c Maximum number of levels for input files`
19	`integer nlevmax`	19	`integer nlevmax`
20	`parameter (nlevmax=100)`	20	`parameter (nlevmax=200)`
21		21
22	`c Maximum number of input files (dates, length of trajectories)`	22	`c Maximum number of input files (dates, length of trajectories)`
23	`integer ndatmax`	23	`integer ndatmax`
24	`parameter (ndatmax=500)`	24	`parameter (ndatmax=500)`
25		25
26	`c Numerical epsilon (for float comparison)`	26	`c Numerical epsilon (for float comparison)`
27	`real eps`	27	`real eps`
28	`parameter (eps=0.001)`	28	`parameter (eps=0.001)`
29		29
30	`c Distance in m between 2 lat circles`	30	`c Distance in m between 2 lat circles`
31	`real deltay`	31	`real deltay`
32	`parameter (deltay=1.112E5)`	32	`parameter (deltay=1.112E5)`
33		33
34	`c Numerical method for the integration (0=iterative Euler, 1=Runge-Kutta)`	34	`c Numerical method for the integration (0=iterative Euler, 1=Runge-Kutta)`
35	`integer imethod`	35	`integer imethod`
36	`parameter (imethod=1)`	36	`parameter (imethod=1)`
37		37
38	`c Number of iterations for iterative Euler scheme`	38	`c Number of iterations for iterative Euler scheme`
39	`integer numit`	39	`integer numit`
40	`parameter (numit=3)`	40	`parameter (numit=3)`
41		41
42	`c Input and output format for trajectories (see iotra.f)`	42	`c Input and output format for trajectories (see iotra.f)`
43	`integer inpmode`	43	`integer inpmode`
44	`integer outmode`	44	`integer outmode`
45		45
46	`c Filename prefix for primary and secondary files (typically 'P' and 'S')`	46	`c Filename prefix for primary and secondary files (typically 'P' and 'S')`
47	`character*1 prefix`	47	`character*1 prefix`
48	`parameter (prefix='P')`	48	`parameter (prefix='P')`
49	`character*1 srefix`	49	`character*1 srefix`
50	`parameter (srefix='S')`	50	`parameter (srefix='S')`
51		51
52	`c Physical constants - needed to compute potential temperature`	52	`c Physical constants - needed to compute potential temperature`
53	`real rdcp,tzero`	53	`real rdcp,tzero`
54	`data rdcp,tzero /0.286,273.15/`	54	`data rdcp,tzero /0.286,273.15/`
55		55
56	`c --------------------------------------------------------------------`	56	`c --------------------------------------------------------------------`
57	`c Declaration of variables`	57	`c Declaration of variables`
58	`c --------------------------------------------------------------------`	58	`c --------------------------------------------------------------------`
59		59
60	`c Input parameters`	60	`c Input parameters`
61	`integer fbflag ! Flag for forward/backward mode`	61	`integer fbflag ! Flag for forward/backward mode`
62	`integer numdat ! Number of input files`	62	`integer numdat ! Number of input files`
63	`character*11 dat(ndatmax) ! Dates of input files`	63	`character*11 dat(ndatmax) ! Dates of input files`
64	`real timeinc ! Time increment between input files`	64	`real timeinc ! Time increment between input files`
65	`real per ! Periodicity (=0 if none)`	65	`real per ! Periodicity (=0 if none)`
66	`integer ntra ! Number of trajectories`	66	`integer ntra ! Number of trajectories`
67	`character*80 cdfname ! Name of output files`	67	`character*80 cdfname ! Name of output files`
68	`real ts ! Time step`	68	`real ts ! Time step`
69	`real tst,ten ! Shift of start and end time relative to first data file`	69	`real tst,ten ! Shift of start and end time relative to first data file`
70	`integer deltout ! Output time interval (in minutes)`	70	`integer deltout ! Output time interval (in minutes)`
71	`integer jflag ! Jump flag (if =1 ground-touching trajectories reenter atmosphere)`	71	`integer jflag ! Jump flag (if =1 ground-touching trajectories reenter atmosphere)`
72	`real wfactor ! Factor for vertical velocity field`	72	`real wfactor ! Factor for vertical velocity field`
73	`character*80 strname ! File with start positions`	73	`character*80 strname ! File with start positions`
74	`character*80 timecheck ! Either 'yes' or 'no'`	74	`character*80 timecheck ! Either 'yes' or 'no'`
75	`character*80 isen ! Isentropic trajectories ('yes' or 'no')`	75	`character*80 isen ! Isentropic trajectories ('yes' or 'no')`
76	`integer thons ! Isentropic mode: is TH availanle on S`	76	`integer thons ! Isentropic mode: is TH availanle on S`
77	`character*80 modlev ! 2D (model level) trajectories ('yes' or 'no')`	77	`character*80 modlev ! 2D (model level) trajectories ('yes' or 'no')`
78		78
79	`c Trajectories`	79	`c Trajectories`
80	`integer ncol ! Number of columns for insput trajectories`	80	`integer ncol ! Number of columns for insput trajectories`
81	`real,allocatable, dimension (:,:,:) :: trainp ! Input start coordinates (ntra,1,ncol)`	81	`real,allocatable, dimension (:,:,:) :: trainp ! Input start coordinates (ntra,1,ncol)`
82	`real,allocatable, dimension (:,:,:) :: traout ! Output trajectories (ntra,ntim,4)`	82	`real,allocatable, dimension (:,:,:) :: traout ! Output trajectories (ntra,ntim,4)`
83	`integer reftime(6) ! Reference date`	83	`integer reftime(6) ! Reference date`
84	`character*80 vars(200) ! Field names`	84	`character*80 vars(200) ! Field names`
85	`real,allocatable, dimension (:) :: xx0,yy0,pp0 ! Position of air parcels`	85	`real,allocatable, dimension (:) :: xx0,yy0,pp0 ! Position of air parcels`
86	`integer,allocatable, dimension (:) :: leftflag ! Flag for domain-leaving`	86	`integer,allocatable, dimension (:) :: leftflag ! Flag for domain-leaving`
87	`real xx1,yy1,pp1 ! Updated position of air parcel`	87	`real xx1,yy1,pp1 ! Updated position of air parcel`
88	`integer leftcount ! Number of domain leaving trajectories`	88	`integer leftcount ! Number of domain leaving trajectories`
89	`integer ntim ! Number of output time steps`	89	`integer ntim ! Number of output time steps`
90	`real,allocatable, dimension (:) :: theta ! Potential temperature for isentropic trajectories`	90	`real,allocatable, dimension (:) :: theta ! Potential temperature for isentropic trajectories`
91	`real,allocatable, dimension (:) :: zindex ! Vertical index for model-level (2D) trajectories`	91	`real,allocatable, dimension (:) :: zindex ! Vertical index for model-level (2D) trajectories`
92		92
93	`c Meteorological fields`	93	`c Meteorological fields`
94	`real,allocatable, dimension (:) :: spt0,spt1 ! Surface pressure`	94	`real,allocatable, dimension (:) :: spt0,spt1 ! Surface pressure`
95	`real,allocatable, dimension (:) :: uut0,uut1 ! Zonal wind`	95	`real,allocatable, dimension (:) :: uut0,uut1 ! Zonal wind`
96	`real,allocatable, dimension (:) :: vvt0,vvt1 ! Meridional wind`	96	`real,allocatable, dimension (:) :: vvt0,vvt1 ! Meridional wind`
97	`real,allocatable, dimension (:) :: wwt0,wwt1 ! Vertical wind`	97	`real,allocatable, dimension (:) :: wwt0,wwt1 ! Vertical wind`
98	`real,allocatable, dimension (:) :: p3t0,p3t1 ! 3d-pressure`	98	`real,allocatable, dimension (:) :: p3t0,p3t1 ! 3d-pressure`
99	`real,allocatable, dimension (:) :: tht0,tht1 ! 3d potential temperature`	99	`real,allocatable, dimension (:) :: tht0,tht1 ! 3d potential temperature`
100	`real,allocatable, dimension (:) :: sth0,sth1 ! Surface potential temperature`	100	`real,allocatable, dimension (:) :: sth0,sth1 ! Surface potential temperature`
101		101
102	`c Grid description`	102	`c Grid description`
103	`real pollon,pollat ! Longitude/latitude of pole`	103	`real pollon,pollat ! Longitude/latitude of pole`
104	`real ak(nlevmax) ! Vertical layers and levels`	104	`real ak(nlevmax) ! Vertical layers and levels`
105	`real bk(nlevmax)`	105	`real bk(nlevmax)`
106	`real xmin,xmax ! Zonal grid extension`	106	`real xmin,xmax ! Zonal grid extension`
107	`real ymin,ymax ! Meridional grid extension`	107	`real ymin,ymax ! Meridional grid extension`
108	`integer nx,ny,nz ! Grid dimensions`	108	`integer nx,ny,nz ! Grid dimensions`
109	`real dx,dy ! Horizontal grid resolution`	109	`real dx,dy ! Horizontal grid resolution`
110	`integer hem ! Flag for hemispheric domain`	110	`integer hem ! Flag for hemispheric domain`
111	`real mdv ! Missing data value`	111	`real mdv ! Missing data value`
112		112
113	`c Auxiliary variables`	113	`c Auxiliary variables`
114	`real delta,rd`	114	`real delta,rd`
115	`integer itm,iloop,i,j,k,filo,lalo`	115	`integer itm,iloop,i,j,k,filo,lalo`
116	`integer ierr,stat`	116	`integer ierr,stat`
117	`integer cdfid,fid`	117	`integer cdfid,fid`
118	`real tstart,time0,time1,time`	118	`real tstart,time0,time1,time`
119	`real reltpos0,reltpos1`	119	`real reltpos0,reltpos1`
120	`real xind,yind,pind,pp,sp,stagz`	120	`real xind,yind,pind,pp,sp,stagz`
121	`character*80 filename,varname`	121	`character*80 filename,varname`
122	`integer reftmp(6)`	122	`integer reftmp(6)`
123	`character ch`	123	`character ch`
124	`real frac,tload`	124	`real frac,tload`
125	`integer itim`	125	`integer itim`
126	`integer wstep`	126	`integer wstep`
127	`real x1,y1,p1`	127	`real x1,y1,p1`
128	`real thetamin,thetamax`	128	`real thetamin,thetamax`
129	`real zindexmin,zindexmax`	129	`real zindexmin,zindexmax`
130		130
131	`c Externals`	131	`c Externals`
132	`real int_index4`	132	`real int_index4`
133	`external int_index4`	133	`external int_index4`
134		134
135	`c --------------------------------------------------------------------`	135	`c --------------------------------------------------------------------`
136	`c Start of program, Read parameters`	136	`c Start of program, Read parameters`
137	`c --------------------------------------------------------------------`	137	`c --------------------------------------------------------------------`
138		138
139	`c Write start message`	139	`c Write start message`
140	`print*,'========================================================='`	140	`print*,'========================================================='`
141	`print,' START OF PROGRAM CALTRA *'`	141	`print,' START OF PROGRAM CALTRA *'`
142	`print*`	142	`print*`
143		143
144	`c Open the parameter file`	144	`c Open the parameter file`
145	`open(9,file='caltra.param')`	145	`open(9,file='caltra.param')`
146		146
147	`c Read flag for forward/backward mode (fbflag)`	147	`c Read flag for forward/backward mode (fbflag)`
148	`read(9,*) fbflag`	148	`read(9,*) fbflag`
149		149
150	`c Read number of input files (numdat)`	150	`c Read number of input files (numdat)`
151	`read(9,*) numdat`	151	`read(9,*) numdat`
152	`if (numdat.gt.ndatmax) then`	152	`if (numdat.gt.ndatmax) then`
153	`print*,' ERROR: too many input files ',numdat,ndatmax`	153	`print*,' ERROR: too many input files ',numdat,ndatmax`
154	`goto 993`	154	`goto 993`
155	`endif`	155	`endif`
156		156
157	`c Read list of input dates (dat, sort depending on forward/backward mode)`	157	`c Read list of input dates (dat, sort depending on forward/backward mode)`
158	`if (fbflag.eq.1) then`	158	`if (fbflag.eq.1) then`
159	`do itm=1,numdat`	159	`do itm=1,numdat`
160	`read(9,'(a11)') dat(itm)`	160	`read(9,'(a11)') dat(itm)`
161	`enddo`	161	`enddo`
162	`else`	162	`else`
163	`do itm=numdat,1,-1`	163	`do itm=numdat,1,-1`
164	`read(9,'(a11)') dat(itm)`	164	`read(9,'(a11)') dat(itm)`
165	`enddo`	165	`enddo`
166	`endif`	166	`endif`
167		167
168	`c Read time increment between input files (timeinc)`	168	`c Read time increment between input files (timeinc)`
169	`read(9,*) timeinc`	169	`read(9,*) timeinc`
170		170
171	`C Read if data domain is periodic and its periodicity`	171	`C Read if data domain is periodic and its periodicity`
172	`read(9,*) per`	172	`read(9,*) per`
173		173
174	`c Read the number of trajectories and name of position file`	174	`c Read the number of trajectories and name of position file`
175	`read(9,*) strname`	175	`read(9,*) strname`
176	`read(9,*) ntra`	176	`read(9,*) ntra`
177	`read(9,*) ncol`	177	`read(9,*) ncol`
178	`if (ntra.eq.0) goto 991`	178	`if (ntra.eq.0) goto 991`
179		179
180	`C Read the name of the output trajectory file and set the constants file`	180	`C Read the name of the output trajectory file and set the constants file`
181	`read(9,*) cdfname`	181	`read(9,*) cdfname`
182		182
183	`C Read the timestep for trajectory calculation (convert from minutes to hours)`	183	`C Read the timestep for trajectory calculation (convert from minutes to hours)`
184	`read(9,*) ts`	184	`read(9,*) ts`
185	`ts=ts/60.`	185	`ts=ts/60.`
186		186
187	`C Read shift of start and end time relative to first data file`	187	`C Read shift of start and end time relative to first data file`
188	`read(9,*) tst`	188	`read(9,*) tst`
189	`read(9,*) ten`	189	`read(9,*) ten`
190		190
191	`C Read output time interval (in minutes)`	191	`C Read output time interval (in minutes)`
192	`read(9,*) deltout`	192	`read(9,*) deltout`
193		193
194	`C Read jumpflag (if =1 ground-touching trajectories reenter the atmosphere)`	194	`C Read jumpflag (if =1 ground-touching trajectories reenter the atmosphere)`
195	`read(9,*) jflag`	195	`read(9,*) jflag`
196		196
197	`C Read factor for vertical velocity field`	197	`C Read factor for vertical velocity field`
198	`read(9,*) wfactor`	198	`read(9,*) wfactor`
199		199
200	`c Read the reference time and the time range`	200	`c Read the reference time and the time range`
201	`read(9,*) reftime(1) ! year`	201	`read(9,*) reftime(1) ! year`
202	`read(9,*) reftime(2) ! month`	202	`read(9,*) reftime(2) ! month`
203	`read(9,*) reftime(3) ! day`	203	`read(9,*) reftime(3) ! day`
204	`read(9,*) reftime(4) ! hour`	204	`read(9,*) reftime(4) ! hour`
205	`read(9,*) reftime(5) ! min`	205	`read(9,*) reftime(5) ! min`
206	`read(9,*) reftime(6) ! time range (in min)`	206	`read(9,*) reftime(6) ! time range (in min)`
207		207
208	`c Read flag for 'no time check'`	208	`c Read flag for 'no time check'`
209	`read(9,*) timecheck`	209	`read(9,*) timecheck`
210		210
211	`c Read flag for isentropic trajectories`	211	`c Read flag for isentropic trajectories`
212	`read(9,*) isen, thons`	212	`read(9,*) isen, thons`
213		213
214	`c Read flag for model-level trajectories (2D mode)`	214	`c Read flag for model-level trajectories (2D mode)`
215	`read(9,*) modlev`	215	`read(9,*) modlev`
216		216
217	`c Close the input file`	217	`c Close the input file`
218	`close(9)`	218	`close(9)`
219		219
220	`c Calculate the number of output time steps`	220	`c Calculate the number of output time steps`
221	`ntim = abs(reftime(6)/deltout) + 1`	221	`ntim = abs(reftime(6)/deltout) + 1`
222		222
223	`c Set the formats of the input and output files`	223	`c Set the formats of the input and output files`
224	`call mode_tra(inpmode,strname)`	224	`call mode_tra(inpmode,strname)`
225	`call mode_tra(outmode,cdfname)`	225	`call mode_tra(outmode,cdfname)`
226	`if (outmode.eq.-1) outmode=1`	226	`if (outmode.eq.-1) outmode=1`
227		227
228	`c Write some status information`	228	`c Write some status information`
229	`print*,'---- INPUT PARAMETERS -----------------------------------'`	229	`print*,'---- INPUT PARAMETERS -----------------------------------'`
230	`print*`	230	`print*`
231	`print*,' Forward/Backward : ',fbflag`	231	`print*,' Forward/Backward : ',fbflag`
232	`print*,' #input files : ',numdat`	232	`print*,' #input files : ',numdat`
233	`print*,' First/last input file : ',trim(dat(1)),' ... ',`	233	`print*,' First/last input file : ',trim(dat(1)),' ... ',`
234	`> trim(dat(numdat))`	234	`> trim(dat(numdat))`
235	`print*,' time increment : ',timeinc`	235	`print*,' time increment : ',timeinc`
236	`print*,' Output file : ',trim(cdfname)`	236	`print*,' Output file : ',trim(cdfname)`
237	`print,' Time step (min) : ',60.ts`	237	`print,' Time step (min) : ',60.ts`
238	`write(*,'(a27,f7.2,f7.2)') ' Time shift (start,end) : ',tst,ten`	238	`write(*,'(a27,f7.2,f7.2)') ' Time shift (start,end) : ',tst,ten`
239	`print*,' Output time interval : ',deltout`	239	`print*,' Output time interval : ',deltout`
240	`print*,' Jump flag : ',jflag`	240	`print*,' Jump flag : ',jflag`
241	`print*,' Vertical wind (scale) : ',wfactor`	241	`print*,' Vertical wind (scale) : ',wfactor`
242	`print*,' Trajectory pos file : ',trim(strname)`	242	`print*,' Trajectory pos file : ',trim(strname)`
243	`print*,' # of trajectories : ',ntra`	243	`print*,' # of trajectories : ',ntra`
244	`print*,' # of output timesteps : ',ntim`	244	`print*,' # of output timesteps : ',ntim`
245	`if ( inpmode.eq.-1) then`	245	`if ( inpmode.eq.-1) then`
246	`print*,' Input format : (lon,lat,p)-list'`	246	`print*,' Input format : (lon,lat,p)-list'`
247	`else`	247	`else`
248	`print*,' Input format : ',inpmode`	248	`print*,' Input format : ',inpmode`
249	`endif`	249	`endif`
250	`print*,' Output format : ',outmode`	250	`print*,' Output format : ',outmode`
251	`print*,' Periodicity : ',per`	251	`print*,' Periodicity : ',per`
252	`print*,' Time check : ',trim(timecheck)`	252	`print*,' Time check : ',trim(timecheck)`
253	`print*,' Isentropic trajectories: ',trim(isen),thons`	253	`print*,' Isentropic trajectories: ',trim(isen),thons`
254	`print*,' Model-level trajs (2D) : ',trim(modlev)`	254	`print*,' Model-level trajs (2D) : ',trim(modlev)`
255	`print*`	255	`print*`
256		256
257	`if ( (isen.eq.'yes').and.(modlev.eq.'yes') ) then`	257	`if ( (isen.eq.'yes').and.(modlev.eq.'yes') ) then`
258	`print*,`	258	`print*,`
259	`> ' WARNING: isentropic and 2D mode chosen -> 2D accepted'`	259	`> ' WARNING: isentropic and 2D mode chosen -> 2D accepted'`
260	`print*`	260	`print*`
261	`isen = 'no'`	261	`isen = 'no'`
262	`endif`	262	`endif`
263		263
264	`c Init missing data value`	264	`c Init missing data value`
265	`mdv = -999.`	265	`mdv = -999.`
266		266
267	`print*,'---- FIXED NUMERICAL PARAMETERS -------------------------'`	267	`print*,'---- FIXED NUMERICAL PARAMETERS -------------------------'`
268	`print*`	268	`print*`
269	`print*,' Numerical scheme : ',imethod`	269	`print*,' Numerical scheme : ',imethod`
270	`print*,' Number of iterations : ',numit`	270	`print*,' Number of iterations : ',numit`
271	`print*,' Filename prefix : ',prefix`	271	`print*,' Filename prefix : ',prefix`
272	`print*,' Missing data value : ',mdv`	272	`print*,' Missing data value : ',mdv`
273	`print*`	273	`print*`
274		274
275	`c --------------------------------------------------------------------`	275	`c --------------------------------------------------------------------`
276	`c Read grid parameters, checks and allocate memory`	276	`c Read grid parameters, checks and allocate memory`
277	`c --------------------------------------------------------------------`	277	`c --------------------------------------------------------------------`
278		278
279	`c Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)`	279	`c Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)`
280	`c The negative <-fid> of the file identifier is used as a flag for parameter retrieval`	280	`c The negative <-fid> of the file identifier is used as a flag for parameter retrieval`
281	`filename = prefix//dat(1)`	281	`filename = prefix//dat(1)`
282	`varname = 'U'`	282	`varname = 'U'`
283	`nx = 1`	283	`nx = 1`
284	`ny = 1`	284	`ny = 1`
285	`nz = 1`	285	`nz = 1`
286	`tload = -tst`	286	`tload = -tst`
287	`call input_open (fid,filename)`	287	`call input_open (fid,filename)`
288	`call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`	288	`call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`
289	`> tload,pollon,pollat,rd,rd,nz,rd,rd,rd,timecheck)`	289	`> tload,pollon,pollat,rd,rd,nz,rd,rd,rd,timecheck)`
290	`call input_close(fid)`	290	`call input_close(fid)`
291		291
292	`C Check if the number of levels is too large`	292	`C Check if the number of levels is too large`
293	`if (nz.gt.nlevmax) goto 993`	293	`if (nz.gt.nlevmax) goto 993`
294		294
295	`C Set logical flag for periodic data set (hemispheric or not)`	295	`C Set logical flag for periodic data set (hemispheric or not)`
296	`hem = 0`	296	`hem = 0`
297	`if (per.eq.0.) then`	297	`if (per.eq.0.) then`
298	`delta=xmax-xmin-360.`	298	`delta=xmax-xmin-360.`
299	`if (abs(delta+dx).lt.eps) then ! Program aborts: arrays must be closed`	299	`if (abs(delta+dx).lt.eps) then ! Program aborts: arrays must be closed`
300	`goto 992`	300	`goto 992`
301	`else if (abs(delta).lt.eps) then ! Periodic and hemispheric`	301	`else if (abs(delta).lt.eps) then ! Periodic and hemispheric`
302	`hem=1`	302	`hem=1`
303	`per=360.`	303	`per=360.`
304	`endif`	304	`endif`
305	`else ! Periodic and hemispheric`	305	`else ! Periodic and hemispheric`
306	`hem=1`	306	`hem=1`
307	`endif`	307	`endif`
308		308
309	`C Allocate memory for some meteorological arrays`	309	`C Allocate memory for some meteorological arrays`
310	`allocate(spt0(nx*ny),stat=stat)`	310	`allocate(spt0(nx*ny),stat=stat)`
311	`if (stat.ne.0) print,' error allocating array spt0 *' ! Surface pressure`	311	`if (stat.ne.0) print,' error allocating array spt0 *' ! Surface pressure`
312	`allocate(spt1(nx*ny),stat=stat)`	312	`allocate(spt1(nx*ny),stat=stat)`
313	`if (stat.ne.0) print,' error allocating array spt1 *'`	313	`if (stat.ne.0) print,' error allocating array spt1 *'`
314	`allocate(uut0(nxnynz),stat=stat)`	314	`allocate(uut0(nxnynz),stat=stat)`
315	`if (stat.ne.0) print,' error allocating array uut0 *' ! Zonal wind`	315	`if (stat.ne.0) print,' error allocating array uut0 *' ! Zonal wind`
316	`allocate(uut1(nxnynz),stat=stat)`	316	`allocate(uut1(nxnynz),stat=stat)`
317	`if (stat.ne.0) print,' error allocating array uut1 *'`	317	`if (stat.ne.0) print,' error allocating array uut1 *'`
318	`allocate(vvt0(nxnynz),stat=stat)`	318	`allocate(vvt0(nxnynz),stat=stat)`
319	`if (stat.ne.0) print,' error allocating array vvt0 *' ! Meridional wind`	319	`if (stat.ne.0) print,' error allocating array vvt0 *' ! Meridional wind`
320	`allocate(vvt1(nxnynz),stat=stat)`	320	`allocate(vvt1(nxnynz),stat=stat)`
321	`if (stat.ne.0) print,' error allocating array vvt1 *'`	321	`if (stat.ne.0) print,' error allocating array vvt1 *'`
322	`allocate(wwt0(nxnynz),stat=stat)`	322	`allocate(wwt0(nxnynz),stat=stat)`
323	`if (stat.ne.0) print,' error allocating array wwt0 *' ! Vertical wind`	323	`if (stat.ne.0) print,' error allocating array wwt0 *' ! Vertical wind`
324	`allocate(wwt1(nxnynz),stat=stat)`	324	`allocate(wwt1(nxnynz),stat=stat)`
325	`if (stat.ne.0) print,' error allocating array wwt1 *'`	325	`if (stat.ne.0) print,' error allocating array wwt1 *'`
326	`allocate(p3t0(nxnynz),stat=stat)`	326	`allocate(p3t0(nxnynz),stat=stat)`
327	`if (stat.ne.0) print,' error allocating array p3t0 *' ! Pressure`	327	`if (stat.ne.0) print,' error allocating array p3t0 *' ! Pressure`
328	`allocate(p3t1(nxnynz),stat=stat)`	328	`allocate(p3t1(nxnynz),stat=stat)`
329	`if (stat.ne.0) print,' error allocating array p3t1 *'`	329	`if (stat.ne.0) print,' error allocating array p3t1 *'`
330	`allocate(tht0(nxnynz),stat=stat)`	330	`allocate(tht0(nxnynz),stat=stat)`
331	`if (stat.ne.0) print,' error allocating array tht0 *' ! Potential temperature`	331	`if (stat.ne.0) print,' error allocating array tht0 *' ! Potential temperature`
332	`allocate(tht1(nxnynz),stat=stat)`	332	`allocate(tht1(nxnynz),stat=stat)`
333	`if (stat.ne.0) print,' error allocating array tht1 *'`	333	`if (stat.ne.0) print,' error allocating array tht1 *'`
334	`allocate(sth0(nx*ny),stat=stat)`	334	`allocate(sth0(nx*ny),stat=stat)`
335	`if (stat.ne.0) print,' error allocating array spt0 *' ! Surface potential temperature`	335	`if (stat.ne.0) print,' error allocating array spt0 *' ! Surface potential temperature`
336	`allocate(sth1(nx*ny),stat=stat)`	336	`allocate(sth1(nx*ny),stat=stat)`
337	`if (stat.ne.0) print,' error allocating array sth1 *'`	337	`if (stat.ne.0) print,' error allocating array sth1 *'`
338		338
339	`C Get memory for trajectory arrays`	339	`C Get memory for trajectory arrays`
340	`allocate(trainp(ntra,1,ncol),stat=stat)`	340	`allocate(trainp(ntra,1,ncol),stat=stat)`
341	`if (stat.ne.0) print,' error allocating array trainp *' ! Input start coordinates`	341	`if (stat.ne.0) print,' error allocating array trainp *' ! Input start coordinates`
342	`allocate(traout(ntra,ntim,4),stat=stat)`	342	`allocate(traout(ntra,ntim,4),stat=stat)`
343	`if (stat.ne.0) print,' error allocating array traout *' ! Output trajectories`	343	`if (stat.ne.0) print,' error allocating array traout *' ! Output trajectories`
344	`allocate(xx0(ntra),stat=stat)`	344	`allocate(xx0(ntra),stat=stat)`
345	`if (stat.ne.0) print,' error allocating array xx0 *' ! X position (longitude)`	345	`if (stat.ne.0) print,' error allocating array xx0 *' ! X position (longitude)`
346	`allocate(yy0(ntra),stat=stat)`	346	`allocate(yy0(ntra),stat=stat)`
347	`if (stat.ne.0) print,' error allocating array yy0 *' ! Y position (latitude)`	347	`if (stat.ne.0) print,' error allocating array yy0 *' ! Y position (latitude)`
348	`allocate(pp0(ntra),stat=stat)`	348	`allocate(pp0(ntra),stat=stat)`
349	`if (stat.ne.0) print,' error allocating array pp0 *' ! Pressure`	349	`if (stat.ne.0) print,' error allocating array pp0 *' ! Pressure`
350	`allocate(leftflag(ntra),stat=stat)`	350	`allocate(leftflag(ntra),stat=stat)`
351	`if (stat.ne.0) print,' error allocating array leftflag *' ! Leaving-domain flag`	351	`if (stat.ne.0) print,' error allocating array leftflag *' ! Leaving-domain flag`
352	`allocate(theta(ntra),stat=stat)`	352	`allocate(theta(ntra),stat=stat)`
353	`if (stat.ne.0) print,' error allocating array theta *' ! Potential temperature for isentropic trajectories`	353	`if (stat.ne.0) print,' error allocating array theta *' ! Potential temperature for isentropic trajectories`
354	`allocate(zindex(ntra),stat=stat)`	354	`allocate(zindex(ntra),stat=stat)`
355	`if (stat.ne.0) print,' error allocating array kind *' ! Vertical index for model-level trajectories`	355	`if (stat.ne.0) print,' error allocating array kind *' ! Vertical index for model-level trajectories`
356		356
357	`c Write some status information`	357	`c Write some status information`
358	`print*,'---- CONSTANT GRID PARAMETERS ---------------------------'`	358	`print*,'---- CONSTANT GRID PARAMETERS ---------------------------'`
359	`print*`	359	`print*`
360	`print*,' xmin,xmax : ',xmin,xmax`	360	`print*,' xmin,xmax : ',xmin,xmax`
361	`print*,' ymin,ymax : ',ymin,ymax`	361	`print*,' ymin,ymax : ',ymin,ymax`
362	`print*,' dx,dy : ',dx,dy`	362	`print*,' dx,dy : ',dx,dy`
363	`print*,' pollon,pollat : ',pollon,pollat`	363	`print*,' pollon,pollat : ',pollon,pollat`
364	`print*,' nx,ny,nz : ',nx,ny,nz`	364	`print*,' nx,ny,nz : ',nx,ny,nz`
365	`print*,' per, hem : ',per,hem`	365	`print*,' per, hem : ',per,hem`
366	`print*`	366	`print*`
367		367
368	`c --------------------------------------------------------------------`	368	`c --------------------------------------------------------------------`
369	`c Initialize the trajectory calculation`	369	`c Initialize the trajectory calculation`
370	`c --------------------------------------------------------------------`	370	`c --------------------------------------------------------------------`
371		371
372	`c Read start coordinates from file - Format (lon,lat,lev)`	372	`c Read start coordinates from file - Format (lon,lat,lev)`
373	`if (inpmode.eq.-1) then`	373	`if (inpmode.eq.-1) then`
374	`open(fid,file=strname)`	374	`open(fid,file=strname)`
375	`do i=1,ntra`	375	`do i=1,ntra`
376	`read(fid,*) xx0(i),yy0(i),pp0(i)`	376	`read(fid,*) xx0(i),yy0(i),pp0(i)`
377	`enddo`	377	`enddo`
378	`close(fid)`	378	`close(fid)`
379		379
380	`c Read start coordinates from trajectory file - check consistency of ref time`	380	`c Read start coordinates from trajectory file - check consistency of ref time`
381	`else`	381	`else`
382	`call ropen_tra(cdfid,strname,ntra,1,ncol,reftmp,vars,inpmode)`	382	`call ropen_tra(cdfid,strname,ntra,1,ncol,reftmp,vars,inpmode)`
383	`call read_tra (cdfid,trainp,ntra,1,ncol,inpmode)`	383	`call read_tra (cdfid,trainp,ntra,1,ncol,inpmode)`
384	`do i=1,ntra`	384	`do i=1,ntra`
385	`time = trainp(i,1,1)`	385	`time = trainp(i,1,1)`
386	`xx0(i) = trainp(i,1,2)`	386	`xx0(i) = trainp(i,1,2)`
387	`yy0(i) = trainp(i,1,3)`	387	`yy0(i) = trainp(i,1,3)`
388	`pp0(i) = trainp(i,1,4)`	388	`pp0(i) = trainp(i,1,4)`
389	`enddo`	389	`enddo`
390	`call close_tra(cdfid,inpmode)`	390	`call close_tra(cdfid,inpmode)`
391		391
392	`if ( ( reftime(1).ne.reftmp(1) ).or.`	392	`if ( ( reftime(1).ne.reftmp(1) ).or.`
393	`> ( reftime(2).ne.reftmp(2) ).or.`	393	`> ( reftime(2).ne.reftmp(2) ).or.`
394	`> ( reftime(3).ne.reftmp(3) ).or.`	394	`> ( reftime(3).ne.reftmp(3) ).or.`
395	`> ( reftime(4).ne.reftmp(4) ).or.`	395	`> ( reftime(4).ne.reftmp(4) ).or.`
396	`> ( reftime(5).ne.reftmp(5) ) )`	396	`> ( reftime(5).ne.reftmp(5) ) )`
397	`> then`	397	`> then`
398	`print*,' WARNING: Inconsistent reference times'`	398	`print*,' WARNING: Inconsistent reference times'`
399	`write(*,'(5i8)') (reftime(i),i=1,5)`	399	`write(*,'(5i8)') (reftime(i),i=1,5)`
400	`write(*,'(5i8)') (reftmp (i),i=1,5)`	400	`write(*,'(5i8)') (reftmp (i),i=1,5)`
401	`print*,'Enter a key to proceed...'`	401	`print*,'Enter a key to proceed...'`
402	`stop`	402	`stop`
403	`endif`	403	`endif`
404	`endif`	404	`endif`
405		405
406	`c Set sign of time range`	406	`c Set sign of time range`
407	`reftime(6) = fbflag * reftime(6)`	407	`reftime(6) = fbflag * reftime(6)`
408		408
409	`c Write some status information`	409	`c Write some status information`
410	`print*,'---- REFERENCE DATE---------- ---------------------------'`	410	`print*,'---- REFERENCE DATE---------- ---------------------------'`
411	`print*`	411	`print*`
412	`print*,' Reference time (year) :',reftime(1)`	412	`print*,' Reference time (year) :',reftime(1)`
413	`print*,' (month) :',reftime(2)`	413	`print*,' (month) :',reftime(2)`
414	`print*,' (day) :',reftime(3)`	414	`print*,' (day) :',reftime(3)`
415	`print*,' (hour) :',reftime(4)`	415	`print*,' (hour) :',reftime(4)`
416	`print*,' (min) :',reftime(5)`	416	`print*,' (min) :',reftime(5)`
417	`print*,' Time range :',reftime(6),' min'`	417	`print*,' Time range :',reftime(6),' min'`
418	`print*`	418	`print*`
419		419
420	`C Save starting positions`	420	`C Save starting positions`
421	`itim = 1`	421	`itim = 1`
422	`do i=1,ntra`	422	`do i=1,ntra`
423	`traout(i,itim,1) = 0.`	423	`traout(i,itim,1) = 0.`
424	`traout(i,itim,2) = xx0(i)`	424	`traout(i,itim,2) = xx0(i)`
425	`traout(i,itim,3) = yy0(i)`	425	`traout(i,itim,3) = yy0(i)`
426	`traout(i,itim,4) = pp0(i)`	426	`traout(i,itim,4) = pp0(i)`
427	`enddo`	427	`enddo`
428		428
429	`c Init the flag and the counter for trajectories leaving the domain`	429	`c Init the flag and the counter for trajectories leaving the domain`
430	`leftcount=0`	430	`leftcount=0`
431	`do i=1,ntra`	431	`do i=1,ntra`
432	`leftflag(i)=0`	432	`leftflag(i)=0`
433	`enddo`	433	`enddo`
434		434
435	`C Convert time shifts <tst,ten> from <hh.mm> into fractional time`	435	`C Convert time shifts <tst,ten> from <hh.mm> into fractional time`
436	`call hhmm2frac(tst,frac)`	436	`call hhmm2frac(tst,frac)`
437	`tst = frac`	437	`tst = frac`
438	`call hhmm2frac(ten,frac)`	438	`call hhmm2frac(ten,frac)`
439	`ten = frac`	439	`ten = frac`
440		440
441	`c Get 3D and surface pressure from first data file`	441	`c Get 3D and surface pressure from first data file`
442	`filename = prefix//dat(1)`	442	`filename = prefix//dat(1)`
443	`call input_open (fid,filename)`	443	`call input_open (fid,filename)`
444	`if ( modlev.eq.'no' ) then`	444	`if ( modlev.eq.'no' ) then`
445	`varname = 'P'`	445	`varname = 'P'`
446	`call input_grid`	446	`call input_grid`
447	`> (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`	447	`> (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`
448	`> tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)`	448	`> tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)`
449	`else`	449	`else`
450	`varname = 'P.ML'`	450	`varname = 'P.ML'`
451	`call input_grid`	451	`call input_grid`
452	`> (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`	452	`> (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`
453	`> tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)`	453	`> tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,timecheck)`
454	`endif`	454	`endif`
455	`call input_close(fid)`	455	`call input_close(fid)`
456		456
457	`c Check that all starting positions are above topography`	457	`c Check that all starting positions are above topography`
458	`do i=1,ntra`	458	`do i=1,ntra`
459		459
460	`C Interpolate surface pressure to actual position (from first input file)`	460	`C Interpolate surface pressure to actual position (from first input file)`
461	`x1 = xx0(i)`	461	`x1 = xx0(i)`
462	`y1 = yy0(i)`	462	`y1 = yy0(i)`
463	`call get_index4 (xind,yind,pind,x1,y1,1050.,0.,`	463	`call get_index4 (xind,yind,pind,x1,y1,1050.,0.,`
464	`> p3t1,p3t1,spt1,spt1,3,`	464	`> p3t1,p3t1,spt1,spt1,3,`
465	`> nx,ny,nz,xmin,ymin,dx,dy,mdv)`	465	`> nx,ny,nz,xmin,ymin,dx,dy,mdv)`
466	`sp = int_index4 (spt1,spt1,nx,ny,1,xind,yind,1.,0.,mdv)`	466	`sp = int_index4 (spt1,spt1,nx,ny,1,xind,yind,1.,0.,mdv)`
467		467
468	`c Decide whether to keep the trajectory`	468	`c Decide whether to keep the trajectory`
469	`if ( pp0(i).gt.sp ) then`	469	`if ( pp0(i).gt.sp ) then`
470	`write(*,'(a30,3f10.2)')`	470	`write(*,'(a30,3f10.2)')`
471	`> 'WARNING: starting point below topography ',`	471	`> 'WARNING: starting point below topography ',`
472	`> xx0(i),yy0(i),pp0(i)`	472	`> xx0(i),yy0(i),pp0(i)`
473	`leftflag(i) = 1`	473	`leftflag(i) = 1`
474	`endif`	474	`endif`
475		475
476	`enddo`	476	`enddo`
477		477
478	`c Special handling for isentropic trajectories - read potential`	478	`c Special handling for isentropic trajectories - read potential`
479	`c temperature from S file or calculate it based on temperature and`	479	`c temperature from S file or calculate it based on temperature and`
480	`c pressure from P file; then, calculate for each trajectory its`	480	`c pressure from P file; then, calculate for each trajectory its`
481	`c potential temperature - which will stay fixed over time`	481	`c potential temperature - which will stay fixed over time`
482	`if ( isen.eq.'yes' ) then`	482	`if ( isen.eq.'yes' ) then`
483		483
484	`c Get potential temperature from S file`	484	`c Get potential temperature from S file`
485	`if ( thons.eq.1 ) then`	485	`if ( thons.eq.1 ) then`
486	`filename = srefix//dat(1)`	486	`filename = srefix//dat(1)`
487	`print*,' TH <- ',trim(filename)`	487	`print*,' TH <- ',trim(filename)`
488	`call input_open (fid,filename)`	488	`call input_open (fid,filename)`
489	`varname='TH'`	489	`varname='TH'`
490	`call input_wind`	490	`call input_wind`
491	`> (fid,varname,tht1,tload,stagz,mdv,`	491	`> (fid,varname,tht1,tload,stagz,mdv,`
492	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)`	492	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)`
493	`call input_close(fid)`	493	`call input_close(fid)`
494		494
495	`c Calculate potential temperature from P file`	495	`c Calculate potential temperature from P file`
496	`else`	496	`else`
497	`filename = prefix//dat(1)`	497	`filename = prefix//dat(1)`
498	`print,' TH = T (1000/P)^RDCP <- ',trim(filename)`	498	`print,' TH = T (1000/P)^RDCP <- ',trim(filename)`
499	`call input_open (fid,filename)`	499	`call input_open (fid,filename)`
500	`varname='T'`	500	`varname='T'`

Subversion Repositories lagranto.ecmwf

(root)/trunk/caltra/caltra.f – Rev 19 → 27