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     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"
 
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 
 
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
      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     --------------------------------------------------------------------
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     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           : (x,y,z)-list'
      else
         print*,'  Input format           : ',inpmode
      endif
      print*,'  Output format          : ',outmode
      print*,'  Periodicity            : ',per
      print*,'  Time check             : ',trim(timecheck)
      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     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...'
            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 
      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     -----------------------------------------------------------------------
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)
 
      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       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)
 
        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)
              call euler(
     >               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)
 
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
           traout(i,j,2) = ( traout(i,j,2) - xmin ) / dx + 1.
           traout(i,j,3) = ( traout(i,j,3) - ymin ) / dy + 1.
        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     ------------------------------------------------------------------
 
 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
c     -------------------------------------------------------------------
 
      subroutine euler(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,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
 100  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

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

Subversion Repositories lagranto.wrf

(root)/trunk/caltra/caltra.f – Rev 11 → 15