WebSVN – lagranto.wrf – 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=100)
 c     Maximum number of input files (dates, length of trajectories)
       integer   ndatmax
       parameter	(ndatmax=500)
 c     Numerical epsilon (for float comparison)
       real      eps
       parameter (eps=0.001)
 c     Distance in m between 2 lat circles
       real	deltay
       parameter	(deltay=1.112E5)
 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     Missing data value
       real      mdv
       parameter (mdv = -999. )
 c     Filename prefix (typically 'P')
       character*1 prefix
       parameter   (prefix='P')
 c     Externals
       real       sdis
       external   sdis
 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
       integer                                per             ! Periodicity flag (=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*10		              	     gridflag        ! "ideal" or "real"
+      real                                   aglconst        ! constant level above ground
+      real                                   isen            ! isentropic level
+      character*80                           vertmode        ! Vertical mode (kinematic,aglconst)
 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
 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       ! 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,allocatable, dimension (:,:)   :: mpx,mpy         ! Map scale factor in x,y direction
       real,allocatable, dimension (:,:)   :: lat,lon         ! Map scale factor in x,y direction
 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
       character*80                           radunit
       real                                   lon1,lat1,lon2,lat2
       real                                   scale_max,scale_min
+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 trajectories at constant level above ground
+      read(9,*) aglconst
+c     Read flag for trajectories at constant isentropic level
+      read(9,*) isen
 c     Close the input file
       close(9)
+c     Set the vertical mode of the trajectory calculation
+      vertmode = 'kinematic'
+      if ( aglconst.gt.0. ) then
+         vertmode = 'aglconst'
+      endif
+      if ( isen.gt.0. ) then
+         vertmode = 'isen'
+      endif
 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           : (x,y,z)-list'
       else
          print*,'  Input format           : ',inpmode
       endif
       print*,'  Output format          : ',outmode
       print*,'  Periodicity            : ',per
       print*,'  Time check             : ',trim(timecheck)
+      if ( vertmode.eq.'aglconst' ) then
+         print*,'  Z = const AGL          : ',aglconst
+      endif
+      if ( vertmode.eq.'isen' ) then
+         print*,'  isentropic             : ',isen
+      endif
       print*
       print*,'---- FIXED NUMERICAL PARAMETERS -------------------------'
       print*
       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_wrf(fid,xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz)
       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 hemispheric mode (not yet supported)
       hem = 0
 C     Allocate memory for some meteorological arrays
       allocate(spt0(nx*ny),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array spt0 ***'   ! Surface geopotential height
       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 ***'   ! geopotential height
       allocate(p3t1(nx*ny*nz),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array p3t1 ***'
       allocate(mpx(nx,ny),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array mpx * **'   ! Map scale factor in x
       allocate(mpy(nx,ny),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array mpy ***'    ! Map scale factor in y
       allocate(lon(nx,ny),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array lon ***'    ! Grid -> lon
       allocate(lat(nx,ny),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array lat ***'    ! Gridy -> lat
 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
+c     Get memory for potential temperature (isentropic trajectories)
+      if ( vertmode.eq.'isen' ) then
+          allocate(tht0(nx*ny*nz),stat=stat)
+          if (stat.ne.0) print*,'*** error allocating array tht0 ***'
+          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 sth0 ***'
+          allocate(sth1(nx*ny),stat=stat)
+          if (stat.ne.0) print*,'*** error allocating array sth1 ***'
+      endif
 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*,'  nx,ny,nz      : ',nx,ny,nz
       print*,'  per, hem      : ',per,hem
       print*
 c     --------------------------------------------------------------------
 c     Calculate the map scale factors
 c     --------------------------------------------------------------------
 c     Read lon/lat on the model grid from first data file
       filename = prefix//dat(1)
       call input_open (fid,filename)
       varname='XLONG'
       call input_var_wrf(fid,varname,lon,nx,ny,1)
       varname='XLAT'
       call input_var_wrf(fid,varname,lat,nx,ny,1)
       call input_close(fid)
 c     Get map scale for interior points
       do i=2,nx-1
          do j=2,ny-1
 c           Map scale in x direction
             lon1     = lon(i-1,j)
             lat1     = lat(i-1,j)
             lon2     = lon(i+1,j)
             lat2     = lat(i+1,j)
             radunit  = 'km.haversine'
             mpx(i,j) = 0.5 * 1000. * sdis(lon1,lat1,lon2,lat2,radunit)
 c           Map scale in y direction
             lon1     = lon(i,j-1)
             lat1     = lat(i,j-1)
             lon2     = lon(i,j+1)
             lat2     = lat(i,j+1)
             radunit  = 'km.haversine'
             mpy(i,j) = 0.5 * 1000. * sdis(lon1,lat1,lon2,lat2,radunit)
          enddo
       enddo
 c     Copy map scale for boundary points
       do i=2,nx-1
         mpx(i, 1) = mpx(i,   2)
         mpx(i,ny) = mpx(i,ny-1)
         mpy(i, 1) = mpy(i,   2)
         mpy(i,ny) = mpy(i,ny-1)
       enddo
       do j=2,ny-1
         mpx(1, j) = mpx(2,   j)
         mpx(nx,j) = mpx(nx-1,j)
         mpy(1, j) = mpy(2,   j)
         mpy(nx,j) = mpy(nx-1,j)
       enddo
       mpx(1,1)   = mpx(2,2)
       mpy(1,1)   = mpy(2,2)
       mpx(nx,1)  = mpx(nx-1,2)
       mpy(nx,1)  = mpy(nx-1,2)
       mpx(nx,ny) = mpx(nx-1,ny-1)
       mpy(nx,ny) = mpy(nx-1,ny-1)
       mpx(1,ny)  = mpx(2,ny-1)
       mpy(1,ny)  = mpy(2,ny-1)
 c	  Write some status information
       print*,'---- MAPPING SCALE FACTORS -----------------------------'
       print*
       scale_max = mpx(1,1)
       scale_min = mpx(1,1)
       do i=1,nx
       	do j=1,ny
       	   if ( mpx(i,j).gt.scale_max ) scale_max = mpx(i,j)
       	   if ( mpx(i,j).lt.scale_min ) scale_min = mpx(i,j)
       	 enddo
       enddo
       print*,'  map scale x (min,max) : ',scale_min,scale_max
       scale_max = mpy(1,1)
       scale_min = mpy(1,1)
       do i=1,nx
       	do j=1,ny
       	   if ( mpy(i,j).gt.scale_max ) scale_max = mpy(i,j)
       	   if ( mpy(i,j).lt.scale_min ) scale_min = mpy(i,j)
       	 enddo
       enddo
       print*,'  map scale y (min,max) : ',scale_min,scale_max
       print*
 c     --------------------------------------------------------------------
 c     Initialize the trajectory calculation
 c     --------------------------------------------------------------------
 c     Read start coordinates from file - Format (x,y,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 ( ( vars(2).ne.'x').and.(vars(3).ne.'y') ) then
             print*,' ERROR: starting positions must be in x/y/z format'
             stop
          endif
          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...'
+            print*,' Enter a key to proceed...'
             stop
          endif
       endif
 c     Transform start coordinates from index space to WRF grid
       do i=1,ntra
          xx0(i) = xmin + ( xx0(i) - 1. ) * dx
          yy0(i) = ymin + ( yy0(i) - 1. ) * dy
       enddo
 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
+C     Save starting positions - vertical position will be corrected below for <aglconst> and <isen>
       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     Adjust the starting heights in mode <aglconst>
+      if ( vertmode.eq.'aglconst' ) then
+        filename = prefix//dat(1)
+        call input_open (fid,filename)
+	      varname='geopot'				! geopot.height
+	      call input_var_wrf(fid,varname,p3t0,nx,ny,nz)
+	      varname='geopots'				! geopot.height at surface
+	      call input_var_wrf(fid,varname,spt0,nx,ny,1)
+        call input_close(fid)
+        filename = prefix//dat(2)
+        call input_open (fid,filename)
+	      varname='geopot'				! geopot.height
+	      call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
+	      varname='geopots'				! geopot.height at surface
+	      call input_var_wrf(fid,varname,spt1,nx,ny,1)
+        call input_close(fid)
+        do i=1,ntra
+           call get_index4 (xind,yind,pind,xx0(i),yy0(i),0.,reltpos0,
+     >                      p3t0,p3t1,spt0,spt1,3,
+     >                      nx,ny,nz,xmin,ymin,dx,dy,mdv)
+           sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos0,mdv)
+           traout(i,1,4) = sp + aglconst
+        enddo
+      endif
+c     Adjust the starting heights in mode <isen>
+      if ( vertmode.eq.'isen' ) then
+        filename = prefix//dat(1)
+        call input_open (fid,filename)
+	      varname='geopot'				! geopot.height
+	      call input_var_wrf(fid,varname,p3t0,nx,ny,nz)
+	      varname='geopots'				! geopot.height at surface
+	      call input_var_wrf(fid,varname,spt0,nx,ny,1)
+          varname='TH'					! TH
+          call input_var_wrf(fid,varname,tht0,nx,ny,nz)
+          do i=1,nx*ny
+              sth0(i) = tht0(i)
+          enddo
+        call input_close(fid)
+        filename = prefix//dat(2)
+        call input_open (fid,filename)
+	      varname='geopot'				! geopot.height
+	      call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
+	      varname='geopots'				! geopot.height at surface
+	      call input_var_wrf(fid,varname,spt1,nx,ny,1)
+          varname='TH'					! TH
+          call input_var_wrf(fid,varname,tht1,nx,ny,nz)
+          do i=1,nx*ny
+              sth1(i) = tht1(i)
+          enddo
+        call input_close(fid)
+        do i=1,ntra
+          call get_index4 (xind,yind,pind,xx0(i),yy0(i),isen,reltpos0,
+     >                    tht0,tht1,sth0,sth1,1,
+     >                    nx,ny,nz,xmin,ymin,dx,dy,mdv)
+          traout(i,1,4) = int_index4
+     >          (p3t0,p3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
+        enddo
+      endif
 c     -----------------------------------------------------------------------
 c     Loop to calculate trajectories
 c     -----------------------------------------------------------------------
 c     Write some status information
       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,f7.2)') '  (file,time) : ',
      >                       trim(filename),tload
       call input_open (fid,filename)
       varname='U'					! U
  	  call input_var_wrf(fid,varname,uut1,nx,ny,nz)
 	  varname='V'					! V
 	  call input_var_wrf(fid,varname,vvt1,nx,ny,nz)
 	  varname='W'					! W
 	  call input_var_wrf(fid,varname,wwt1,nx,ny,nz)
 	  varname='geopot'				! geopot.height
 	  call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
 	  varname='geopots'				! geopot.height at surface
 	  call input_var_wrf(fid,varname,spt1,nx,ny,1)
+      if (vertmode.eq.'isen') then
+           varname='TH'					! TH
+ 	       call input_var_wrf(fid,varname,tht1,nx,ny,nz)
+           do i=1,nx*ny
+              sth1(i) = tht1(i)
+           enddo
+      endif
       call input_close(fid)
 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)
         enddo
         do i=1,nx*ny
            spt0(i)=spt1(i)
         enddo
+c       Copy potential temperature
+        if ( vertmode.eq.'isen' ) then
+           do i=1,nx*ny*nz
+              tht0(i)=tht1(i)
+           enddo
+           do i=1,nx*ny
+             sth0(i)=sth1(i)
+           enddo
+        endif
 c       Read wind fields and surface pressure at next time
         filename = prefix//dat(itm+1)
         call frac2hhmm(time1,tload)
         write(*,'(a16,a20,f7.2)') '  (file,time) : ',
      >                          trim(filename),tload
         call input_open (fid,filename)
   	    varname='U'					! U
 	    call input_var_wrf(fid,varname,uut1,nx,ny,nz)
 	    varname='V'					! V
 	    call input_var_wrf(fid,varname,vvt1,nx,ny,nz)
 	    varname='W'					! W
 	    call input_var_wrf(fid,varname,wwt1,nx,ny,nz)
 	    varname='geopot'				! geopot.height
 	    call input_var_wrf(fid,varname,p3t1,nx,ny,nz)
 	    varname='geopots'				! geopot.height
 	    call input_var_wrf(fid,varname,spt1,nx,ny,1)
+        if (vertmode.eq.'isen') then
+           varname='TH'					! TH
+ 	       call input_var_wrf(fid,varname,tht1,nx,ny,nz)
+           do i=1,nx*ny
+              sth1(i) = tht1(i)
+           enddo
+        endif
         call input_close(fid)
 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         Initialize counter for domain leaving trajectories
           leftcount=0
 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             Iterative Euler timestep (x0,y0,p0 -> x1,y1,p1)
+              if ( vertmode.eq.'kinematic' ) then
-              call euler(
+                 call euler_kinematic(
      >               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,mpx,mpy)
+              elseif ( vertmode.eq.'aglconst' ) then
+                 call euler_aglconst(
+     >               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,
+     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy,aglconst)
+              elseif ( vertmode.eq.'isen' ) then
+                 call euler_isen(
+     >               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,
+     >               tht0,tht1,sth0,sth1,
+     >               xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy,isen)
+              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
             time = time0+reltpos1*timeinc*fbflag
             itim = itim + 1
             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
         enddo
       enddo
 c     Transform trajectory position from WRF grid do grid index
       do i=1,ntra
         do j=1,ntim
            if ( abs(traout(i,j,2)-mdv).gt.eps ) then
               traout(i,j,2) = ( traout(i,j,2) - xmin ) / dx + 1.
            else
               traout(i,j,2) = mdv
            endif
            if ( abs(traout(i,j,3)-mdv).gt.eps ) then
               traout(i,j,3) = ( traout(i,j,3) - ymin ) / dy + 1.
            else
               traout(i,j,3) = mdv
            endif
         enddo
       enddo
 c     Write trajectory file
       vars(1)  ='time'
       vars(2)  ='x'
       vars(3)  ='y'
       vars(4)  ='z'
       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     ------------------------------------------------------------------
 write(*,*) '*** ERROR: all start points outside the data domain'
       call exit(1)
 write(*,*) '*** ERROR: close arrays on files (prog. closear)'
       call exit(1)
 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
+c     Iterative Euler time step - kinematic trajectory
 c     -------------------------------------------------------------------
+      subroutine euler_kinematic
-      subroutine euler(x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
+     >          (x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
-     >                 deltat,numit,jump,mdv,wfactor,fbflag,
+     >           deltat,numit,jump,mdv,wfactor,fbflag,
-     >		           spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
+     >		     spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,wwt0,wwt1,
-     >                 xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy)
+     >           xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy)
       implicit none
 c     Flag for test mode
       integer      test
       parameter    (test=0)
 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
       integer      per
       integer      hem
       real         mdv
       real         mpx(nx*ny),mpy(nx*ny)
 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
       real         mpsc_x,mpsc_y
 c     Externals
       real         int_index4
       external     int_index4
 c     Reset the flag for domain-leaving
       left=0
 c     Set the esat-north boundary 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        Get the mapping scale factors for this position
          mpsc_x = int_index4 (mpx,mpx,nx,ny,1,xind,yind,1.,reltpos1,mdv)
          mpsc_y = int_index4 (mpy,mpy,nx,ny,1,xind,yind,1.,reltpos1,mdv)
 C        Calculate new positions (adapted for cartesian grid)
          x1 = x0 + fbflag * deltat * u * dx/mpsc_x
          y1 = y0 + fbflag * deltat * v * dy/mpsc_y
          p1 = p0 + fbflag * deltat * w * wfactor
 c        Write the update positions for tests
          if ( (icount.eq.3).and.(test.eq.1) ) then
             write(*,'(10f10.2)') x0,u,x1-x0,p0,w,p1-p0
          endif
 C       Interpolate surface geop.height to actual position
         call get_index4 (xind,yind,pind,x1,y1,0.,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.lt.sp)) p1=sp+10.
 c       Handle periodioc boundaries for 'ideal' mode
         if ( per.eq.1 ) then
 	        if ( x1.gt.xmax ) x1=x1-xmax
 	        if ( x1.lt.xmin ) x1=x1+xmax
         endif
 C       Check if trajectory leaves data domain
         if ( (x1.lt.xmin).or.(x1.gt.xmax).or.
      >       (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.lt.sp) )   ! geopot : .lt. ; pressure: .gt.
      >  then
           left=1
           print*,x1,y1,p1
           print*,xmin,ymin
 	      print*,xmax,ymax
 	      print*,sp
           goto 100
         endif
       enddo
 c     Exit point for subroutine
 continue
       return
       end
+c     -------------------------------------------------------------------
+c     Iterative Euler time step - constant level above ground
+c     -------------------------------------------------------------------
+      subroutine euler_aglconst
+     >          (x1,y1,p1,left,x0,y0,p0,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,mpx,mpy,aglconst)
+      implicit none
+c     Flag for test mode
+      integer      test
+      parameter    (test=0)
+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
+      integer      per
+      integer      hem
+      real         mdv
+      real         mpx(nx*ny),mpy(nx*ny)
+      real         aglconst
+c     Auxiliary variables
+      real         xmax,ymax
+      real	       xind,yind,pind
+      real	       u0,v0,u1,v1,u,v,w,sp
+      integer	   icount
+      character    ch
+      real         mpsc_x,mpsc_y
+c     Externals
+      real         int_index4
+      external     int_index4
+c     Reset the flag for domain-leaving
+      left=0
+c     Set the esat-north boundary of the domain
+      xmax = xmin+real(nx-1)*dx
+      ymax = ymin+real(ny-1)*dy
+C     Set parcel height to aglconst above topography
+      call get_index4 (xind,yind,pind,x1,y1,0.,reltpos0,
+     >                 p3d0,p3d1,spt0,spt1,3,
+     >                 nx,ny,nz,xmin,ymin,dx,dy,mdv)
+      sp = int_index4 (spt0,spt1,nx,ny,1,xind,yind,1.,reltpos0,mdv)
+      p0 = sp + aglconst
+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        Get the mapping scale factors for this position
+         mpsc_x = int_index4 (mpx,mpx,nx,ny,1,xind,yind,1.,reltpos1,mdv)
+         mpsc_y = int_index4 (mpy,mpy,nx,ny,1,xind,yind,1.,reltpos1,mdv)
+C        Calculate new positions (adapted for cartesian grid)
+         x1 = x0 + fbflag * deltat * u * dx/mpsc_x
+         y1 = y0 + fbflag * deltat * v * dy/mpsc_y
+c        Write the update positions for tests
+         if ( (icount.eq.3).and.(test.eq.1) ) then
+            write(*,'(10f10.2)') x0,u,x1-x0,p0,w,p1-p0
+         endif
+C       Set trajectory height to aglconst above topography
+        call get_index4 (xind,yind,pind,x1,y1,0.,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)
+        p1 = sp + aglconst
+c       Handle trajectories which cross the lower boundary (jump flag)
+        if ((jump.eq.1).and.(p1.lt.sp)) then
+            print*,'jump'
+            p1=sp+10.
+        endif
+c       Handle periodioc boundaries for 'ideal' mode
+        if ( per.eq.1 ) then
+	        if ( x1.gt.xmax ) x1=x1-xmax
+	        if ( x1.lt.xmin ) x1=x1+xmax
+        endif
+C       Check if trajectory leaves data domain
+        if ( (x1.lt.xmin).or.(x1.gt.xmax).or.
+     >       (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.lt.sp) )   ! geopot : .lt. ; pressure: .gt.
+     >  then
+          left=1
+          goto 100
+        endif
+      enddo
+c     Exit point for subroutine
+continue
+      return
+      end
+c     -------------------------------------------------------------------
+c     Iterative Euler time step - isentropic
+c     -------------------------------------------------------------------
+      subroutine euler_isen
+     >          (x1,y1,p1,left,x0,y0,p0,reltpos0,reltpos1,
+     >           deltat,numit,jump,mdv,wfactor,fbflag,
+     >		     spt0,spt1,p3d0,p3d1,uut0,uut1,vvt0,vvt1,
+     >           tht0,tht1,sth0,sth1,
+     >           xmin,ymin,dx,dy,per,hem,nx,ny,nz,mpx,mpy,isen)
+      implicit none
+c     Flag for test mode
+      integer      test
+      parameter    (test=0)
+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         tht0(nx*ny*nz),tht1(nx*ny*nz),sth0(nx*ny),sth1(nx*ny)
+      real         p3d0(nx*ny*nz),p3d1(nx*ny*nz)
+      real         xmin,ymin,dx,dy
+      integer      per
+      integer      hem
+      real         mdv
+      real         mpx(nx*ny),mpy(nx*ny)
+      real         isen
+c     Auxiliary variables
+      real         xmax,ymax
+      real	       xind,yind,pind
+      real	       u0,v0,u1,v1,u,v,w,sp
+      integer	   icount
+      character    ch
+      real         mpsc_x,mpsc_y
+c     Externals
+      real         int_index4
+      external     int_index4
+c     Reset the flag for domain-leaving
+      left=0
+c     Set the esat-north boundary 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        Get the mapping scale factors for this position
+         mpsc_x = int_index4 (mpx,mpx,nx,ny,1,xind,yind,1.,reltpos1,mdv)
+         mpsc_y = int_index4 (mpy,mpy,nx,ny,1,xind,yind,1.,reltpos1,mdv)
+C        Calculate new positions (adapted for cartesian grid)
+         x1 = x0 + fbflag * deltat * u * dx/mpsc_x
+         y1 = y0 + fbflag * deltat * v * dy/mpsc_y
+c        Write the update positions for tests
+         if ( (icount.eq.3).and.(test.eq.1) ) then
+            write(*,'(10f10.2)') x0,u,x1-x0,p0,w,p1-p0
+         endif
+c        Set trajectory height to isentropic level
+         call get_index4 (xind,yind,pind,x1,y1,isen,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       Interpolate surface geop.height to actual position
+        call get_index4 (xind,yind,pind,x1,y1,0.,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.lt.sp)) then
+            print*,'jump'
+            p1=sp+10.
+        endif
+c       Handle periodioc boundaries for 'ideal' mode
+        if ( per.eq.1 ) then
+	        if ( x1.gt.xmax ) x1=x1-xmax
+	        if ( x1.lt.xmin ) x1=x1+xmax
+        endif
+C       Check if trajectory leaves data domain
+        if ( (x1.lt.xmin).or.(x1.gt.xmax).or.
+     >       (y1.lt.ymin).or.(y1.gt.ymax).or.(p1.lt.sp) )   ! geopot : .lt. ; pressure: .gt.
+     >  then
+          left=1
+          goto 100
+        endif
+      enddo
+c     Exit point for subroutine
+continue
+      return
+      end
 c     ----------------------------------------------------------------------
 c     Get spherical distance between lat/lon points
 c     ----------------------------------------------------------------------
       real function sdis(xp,yp,xq,yq,unit)
 c     Calculates spherical distance (in km) between two points given
 c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
       real,parameter ::       re=6371.229
       real,parameter ::       rad2deg=180./3.14159265
       real,parameter ::       deg2rad=3.141592654/180.
       real         xp,yp,xq,yq,arg
       character*80 unit
       real         dlon,dlat,alpha,rlat1,ralt2
       if ( unit.eq.'km' ) then
          arg=sin(yp*deg2rad)*sin(yq*deg2rad)+
      >              cos(yp*deg2rad)*cos(yq*deg2rad)*cos((xp-xq)*deg2rad)
          if (arg.lt.-1.) arg=-1.
          if (arg.gt.1.) arg=1.
          sdis=re*acos(arg)
       elseif ( unit.eq.'deg' ) then
          dlon = xp-xq
          if ( dlon.gt. 180. ) dlon = dlon - 360.
          if ( dlon.lt.-180. ) dlon = dlon + 360.
          sdis = sqrt( dlon**2 + (yp-yq)**2 )
       elseif ( unit.eq.'km.haversine' ) then
         dlat   =  (yp - yq)*deg2rad
         dlon   =  (xp - xq)*deg2rad
         rlat1   =  yp*deg2rad
         rlat2   =  yq*deg2rad
         alpha  = ( sin(0.5*dlat)**2 ) +
      >           ( sin(0.5*dlon)**2 )*cos(rlat1)*cos(rlat2)
         sdis = 4. * re * atan2( sqrt(alpha),1. + sqrt( 1. - alpha ) )
       endif
       end

Subversion Repositories lagranto.wrf

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