WebSVN – lagranto.um – Diff – /trunk/startf/startf.f/


      PROGRAM startf
 
c     **************************************************************
c     * Create a <startfile> for <lagrangto>. It can be chosen     *
c     * whether to start from an isentropic or an isobaric         *
c     * surface. The starting points are equidistantly distributed *
c     * Michael Sprenger / Autumn 2004                             *
c     **************************************************************
 
      implicit none
 
 
c     --------------------------------------------------------------
c     Set parameters
c     --------------------------------------------------------------
 
c     Maximum number of starting positions
      integer           nmax
      parameter         (nmax=1000000)
 
c     Maximum number of model levels
      integer           nlevmax
      parameter         (nlevmax=200)
      
c     Grid constant (distance in km corresponding to 1 deg at the equator)
      real              deltat
      parameter         (deltat=111.)
 
c     Mathematical constant (conversion degree -> radian)
      real              pi180
      parameter         (pi180=3.14159/180.)
 
c     Numerical epsilon
      real              eps
      parameter         (eps=0.00001)
 
c     --------------------------------------------------------------
c     Set variables
c     --------------------------------------------------------------
 
c     Filenames and output format
      character*80      pfile0,pfile1                  ! P filenames
      character*80      sfile0,sfile1                  ! S filenames
      character*80      ofile                          ! Output filename
      integer           oformat                        ! Output format
      real              timeshift                      ! Time shift relative to data files <*0>
      real              timeinc                        ! Time increment between input files
 
c     Horizontal grid
      character*80      hmode                          ! Horizontale mode
      real              lat1,lat2,lon1,lon2            ! Lat/lon boundaries
      real              ds,dlon,dlat                   ! Distance and lat/lon shifts
      character*80      hfile                          ! Filename
      integer           hn                             ! Number of entries in lat/lon list 
      real              latlist(nmax)                  ! List of latitudes
      real              lonlist(nmax)                  ! List of longitudes
      integer           pn                             ! Number of entries in lat/lon poly
      real              latpoly(500)                   ! List of polygon latitudes
      real              lonpoly(500)                   ! List of polygon longitudes
      real              loninpoly,latinpoly            ! Lon/lat inside polygon
      character*80      regionf                        ! Region file
      integer           iregion                        ! Region number
      real              xcorner(4),ycorner(4)          ! Vertices of region
 
c     Vertical grid
      character*80      vmode                          ! Vertical mode
      real              lev1,lev2,levlist(nmax)        ! Single levels, and list of levels
      character*80      vfile                          ! Filename
      integer           vn                             ! Number of entries
 
c     Unit of vertical axis
      character*80      umode                          ! Unit of vertical axis
 
c     Flag for 'no time check'
      character*80      timecheck                      ! Either 'no' or 'yes'
 
c     List of all starting positions
      integer           start_n                        ! Number of coordinates
      real              start_lat(nmax)                ! Latitudes
      real              start_lon(nmax)                ! Longitudes
      real              start_lev(nmax)                ! Levels (depending on vertical unit)
      real              start_pre(nmax)                ! Level in hPa
      integer           reftime(6)                     ! Reference time
      character*80      vars(10)                       ! Name of output fields (time,lon,lat,p)
      real,allocatable, dimension (:,:,:) :: tra       ! Trajectories (ntra,ntim,ncol)
      real              latmin,latmax
      real              lonmin,lonmax
      real              premin,premax
 
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
      real,allocatable, dimension (:,:,:) :: pr        ! 3d pressure
      real,allocatable, dimension (:,:)   :: prs       ! surface pressure
      real,allocatable, dimension (:,:,:) :: th        ! 3d potential temperature
      real,allocatable, dimension (:,:)   :: ths       ! surface poential temperature
      real,allocatable, dimension (:,:,:) :: pv        ! 3d potential vorticity
      real,allocatable, dimension (:,:)   :: pvs       ! surface potential vorticiy
      character*80      varname                        ! Name of input variable
      integer           fid                            ! File identifier
      real              stagz                          ! Vertical staggering
      real              mdv                            ! Missing data values
      real              tstart,tend                    ! Time on P and S file
      real              rid,rjd,rkd                    ! Real grid position
      integer           rotategrid                     ! Flag whether to rotate grid
 
c     Auxiliary variable
      integer           i,j,k
      real              lon,lat,rlon,rlat
      real              rd
      integer           stat,flag
      real              tmp1,tmp2
      real              tfrac,frac
      real              radius,dist
      character*80      string
      character*80      selectstr  
      integer           noutside 
      real,allocatable, dimension (:,:,:) :: fld0
      real,allocatable, dimension (:,:,:) :: fld1
      real,allocatable, dimension (:,:  ) :: sfc0
      real,allocatable, dimension (:,:)   :: sfc1
 
c     Externals 
      real              int_index3     ! 3d interpolation
      external          int_index3
      real              sdis           ! Speherical distance
      external          sdis
      integer           inregion       ! In/out of region
      external          inrehion
      real              lmtolms        ! Grid rotation
      real              phtophs    
      real              lmstolm
      real              phstoph        
      external          lmtolms,phtophs,lmstolm,phstoph
 
c     ------------------------------------------------------------------
c     Start of program, Read parameters
c     ------------------------------------------------------------------
 
c     Write start message
      print*,'========================================================='
      print*,'         *** START OF PROGRAM CREATE_STARTF ***'
      print*
 
c     Read parameter file
      open(10,file='startf.param')
 
c      Input P and S file
       read(10,*) pfile0,pfile1
       read(10,*) sfile0,sfile1
       read(10,*) ofile
 
c      Read name of region file
       read(10,*) regionf
 
c      Reference time
       do i=1,6
          read(10,*) reftime(i)
       enddo
 
c      Time shift relative to data files <pfile0,sfile0> - format (hh.mm)
       read(10,*) timeshift
 
c      Read timeincrement between input files
       read(10,*) timeinc
 
c      Parameters for horizontal grid
       read(10,*) hmode
       if ( hmode.eq.'file' ) then                  ! from file
          read(10,*) hfile
       elseif ( hmode.eq.'line' ) then              ! along a line
          read(10,*) lon1,lon2,lat1,lat2,hn
       elseif ( hmode.eq.'box.eqd' ) then           ! box: 2d equidistant
          read(10,*) lon1,lon2,lat1,lat2,ds
       elseif ( hmode.eq.'box.grid' ) then          ! box: 2d grid
          read(10,*) lon1,lon2,lat1,lat2
       elseif ( hmode.eq.'point' ) then             ! single point
          read(10,*) lon1,lat1
       elseif ( hmode.eq.'shift' ) then             ! centre + shifted
          read(10,*) lon1,lat1,dlon,dlat
       elseif ( hmode.eq.'polygon.eqd' ) then       ! polygon: 2d equidistant
          read(10,*) hfile,ds
       elseif ( hmode.eq.'polygon.grid' ) then      ! polygon: 2d grid
          read(10,*) hfile
       elseif ( hmode.eq.'circle.eqd' ) then        ! circle: 2d equidistant
          read(10,*) lon1,lat1,radius,ds 
       elseif ( hmode.eq.'circle.grid' ) then       ! circle: 2d grid
          read(10,*) lon1,lat1,radius
       elseif ( hmode.eq.'region.eqd' ) then        ! region: 2d equidistant
          read(10,*) iregion,ds 
       elseif ( hmode.eq.'region.grid' ) then       ! iregion: 2d grid
          read(10,*) iregion
       else
          print*,' ERROR: horizontal mode not supported ',trim(hmode)
          stop
       endif
 
c      Parameters for vertical grid
       read(10,*) vmode
       if ( vmode.eq.'file') then                   ! from file
          read(10,*) vfile
       elseif ( vmode.eq.'level' ) then             ! single level (explicit command)
          read(10,*) lev1                         
       elseif ( vmode.eq.'list') then               ! a list
          read(10,*) vn
          read(10,*) (levlist(i),i=1,vn)
       elseif ( vmode.eq.'profile') then            ! a profile
          read(10,*) lev1,lev2,vn
       elseif ( vmode.eq.'nil') then                ! assume that vertzical positions are on file
          continue
       else
          print*,' ERROR: vertical mode not supported ',trim(vmode)
          stop
       endif
 
c      Read units of vertical axis
       read(10,*) umode
       if ( ( umode.ne.'hPa'     ).and.
     >      ( umode.ne.'hPa,agl' ).and.   
     >      ( umode.ne.'K'       ).and.
     >      ( umode.ne.'PVU'     ) ) 
     > then  
          print*,' ERROR: unit not supported ',trim(umode)
          stop
       endif
 
c     Read select string (dummy read)
      read(10,*) selectstr
 
c     Read flag for 'no time check'
      read(10,*) timecheck 
 
c     Close parameter file
      close(10)
 
c     Decide which output format is used (1..4: trajectory format, -1: triple list)
      call mode_tra(oformat,ofile)
      
c     Decide whether all lat/lon/lev coordaintes are read from one file
      if ( (hmode.eq.'file').and.(vmode.eq.'nil') ) then
         hmode='file3'
      elseif ( (hmode.eq.'file').and.(vmode.ne.'nil') ) then
         hmode='file2'
      endif
 
c     Convert timeshift (hh.mm) into a fractional time shift
      call hhmm2frac(timeshift,tfrac)
      if (tfrac.gt.0.) then
         tfrac=tfrac/timeinc
      else
         tfrac=0.
      endif
 
c     Read the region coordinates if needed
      if ( (hmode.eq.'region.eqd' ).or.
     >     (hmode.eq.'region.grid') ) then
         
         open(10,file=regionf)
          
 50       read(10,*,end=51) string
          
          if ( string(1:1).ne.'#' ) then
             call regionsplit(string,i,xcorner,ycorner)
             if ( i.eq.iregion ) goto 52
          endif
 
          goto 50
 
 51      close(10)
         
         print*,' ERROR: region ',iregion,' not found on ',trim(regionf)
         stop
 
 52      continue
          
      endif
 
c     Write some status information
      print*,'---- INPUT PARAMETERS -----------------------------------'
      print*
      if ( timeshift.gt.0. ) then
         print*,'  P file                     : ',trim(pfile0),
     >                                            ' ',
     >                                            trim(pfile1)
         print*,'  S file                     : ',trim(sfile0),
     >                                            ' ',
     >                                            trim(sfile1)
      else
         print*,'  P file                     : ',trim(pfile0)
         print*,'  S file                     : ',trim(sfile0)
      endif
      print*,'  Output file                : ',trim(ofile) 
      print*
      if (oformat.eq.-1) then
         print*,'  Output format              : (lon,lat,lev)-list'
      else
         print*,'  Output format              : ',oformat
      endif
      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)
      print*
      print*,'  Time shift                 : ',timeshift,' + ',
     >                                         trim(pfile0)
      print*,'  Region file                : ',trim(regionf)
      print*
      print*,'  hmode                      : ',trim(hmode)
      if ( hmode.eq.'file2' ) then        
        print*,'      filename [lat/lon]      : ',trim(hfile)
      elseif ( hmode.eq.'file3' ) then        
        print*,'      filename [lat/lon/lev]  : ',trim(hfile)
      elseif ( hmode.eq.'line' ) then   
        write(*,'(a30,4f10.2,i4)') 
     >         '      lon1,lon2,lat1,lat2,n   : ',lon1,lon2,lat1,lat2,hn
      elseif ( hmode.eq.'box.eqd' ) then     
        write(*,'(a30,5f10.2)') 
     >         '      lon1,lon2,lat1,lat2,ds  : ',lon1,lon2,lat1,lat2,ds
      elseif ( hmode.eq.'box.grid' ) then    
        write(*,'(a30,4f10.2)') 
     >         '      lon1,lon2,lat1,lat2     : ',lon1,lon2,lat1,lat2 
      elseif ( hmode.eq.'point' ) then   
        print*,'      lon,lat                 : ',lon1,lat1
      elseif ( hmode.eq.'shift' ) then   
        write(*,'(a30,4f10.2)') 
     >         '      lon,lat,dlon,dlat       : ',lon1,lat1,dlon,dlat
      elseif ( hmode.eq.'polygon.eqd' ) then   
        write(*,'(a30,a10,f10.2)') 
     >         '      hfile, ds               : ',trim(hfile),ds
      elseif ( hmode.eq.'polygon.grid' ) then   
        write(*,'(a30,a10)') 
     >         '      hfile                   : ',trim(hfile)
      elseif ( hmode.eq.'circle.eqd' ) then   
        write(*,'(a30,4f10.2)') 
     >         '      lonc,latc,radius, ds    : ',lon1,lat1,radius,ds
      elseif ( hmode.eq.'circle.grid' ) then   
        write(*,'(a30,3f10.2)') 
     >         '      lonc,latc,radius        : ',lon1,lat1,radius
      elseif ( hmode.eq.'region.eqd' ) then   
        write(*,'(a30,i4,1f10.2)') 
     >         '      iregion, ds             : ',iregion,ds
        write(*,'(a30,4f10.2)')
     >         '      xcorner                 : ',(xcorner(i),i=1,4)
        write(*,'(a30,4f10.2)')
     >         '      ycorner                 : ',(ycorner(i),i=1,4)
      elseif ( hmode.eq.'region.grid' ) then   
        write(*,'(a30,i4)') 
     >         '      iregion                 : ',iregion
        write(*,'(a30,4f10.2)')
     >         '      xcorner                 : ',(xcorner(i),i=1,4)
        write(*,'(a30,4f10.2)')
     >         '      ycorner                 : ',(ycorner(i),i=1,4)
      endif
      print*
      print*,'  vmode                      : ',trim(vmode)
      if ( vmode.eq.'file') then 
         print*,'      filename               : ',trim(vfile)
      elseif ( vmode.eq.'level' ) then 
         print*,'      level                  : ',lev1                         
      elseif ( vmode.eq.'list') then 
         print*,'      n                      : ',vn
         print*,'      level(i)               : ',(levlist(i),i=1,vn)
      elseif ( vmode.eq.'profile') then 
         print*,'      lev1,lev2,n            : ',lev1,lev2,vn
      endif
      print* 
      print*,'  umode                      : ',trim(umode)
      print*
      print*,'  time check                 : ',trim(timecheck)
      print*
 
c     ------------------------------------------------------------------
c     Read grid parameters from inital files
c     ------------------------------------------------------------------
 
c     Init missing data value
      mdv = -999.0
 
c     Get the time of the first and second data file
      tstart = -timeshift                              ! Format hh.mm
      call hhmm2frac(tstart,frac)
      frac   = frac + timeinc          
      call frac2hhmm(frac,tend)                        ! Format hh.mm                  
 
c     Convert timeshift (hh.mm) into a fractional time shift
      tfrac=real(int(timeshift))+
     >          100.*(timeshift-real(int(timeshift)))/60.
      if (tfrac.gt.0.) then
         tfrac=tfrac/timeinc
      else
         tfrac=0.
      endif
 
c     Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,
c     pollon,pollat) The negative <-fid> of the file identifier is used 
c     as a flag for parameter retrieval
      varname  = 'U'
      call input_open (fid,pfile0)
      call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >                 tstart,pollon,pollat,rd,rd,nz,rd,rd,rd,timecheck)
      call input_close(fid)
 
c     Decide whether the grid has to be rotated
      if ( abs(pollat-90.).gt.eps) then
         rotategrid = 1
      else
         rotategrid = 0
      endif
 
c     Check whether region coordinates are within the domain
      if ( (hmode.eq.'region.eqd' ).or.
     >     (hmode.eq.'region.grid') ) then
 
c        Rotate coordinates
         if ( rotategrid.eq.1) then
            do i=1,4
               lon = xcorner(i)
               lat = ycorner(i)
               if ( rotategrid.eq.1) then
                  rlon = lmtolms(lat,lon,pollat,pollon)
                  rlat = phtophs(lat,lon,pollat,pollon)
               else
                  rlon = lon
                  rlat = lat
               endif 
               xcorner(i) = rlon
               ycorner(i) = rlat
            enddo
         endif
         
c        Check 
         do i=1,4
            if ( (xcorner(i).lt.xmin).or.
     >           (ycorner(i).lt.ymin).or.
     >           (xcorner(i).gt.xmax).or.
     >           (ycorner(i).gt.ymax) )
     >      then
               print*,' ERROR: region not included in data domain...'
               print*,'       ',trim(string)
               print*,'  x:   ',(xcorner(j),j=1,4)
               print*,'  y:   ',(ycorner(j),j=1,4)
               stop
            endif
               
         enddo
 
c        Write rotated coordinates to logfile
         if ( rotategrid.eq.1) then
            write(*,'(a30,i4)') 
     >         '  Rotated region              : ',iregion
            write(*,'(a30,4f10.2)')
     >         '      xcorner (rotated)       : ',(xcorner(i),i=1,4)
            write(*,'(a30,4f10.2)')
     >         '      ycorner (rotated)       : ',(ycorner(i),i=1,4)
            print*
         endif
 
c        Rotate coordinates back
         if ( rotategrid.eq.1 ) then
            do i=1,4
               rlon = xcorner(i)
               rlat = ycorner(i)
               if ( rotategrid.eq.1) then
                  lon = lmstolm(rlat,rlon,pollat,pollon)
                  lat = phstoph(rlat,rlon,pollat,pollon)
               else
                  lon = rlon
                  lat = rlat
               endif 
               xcorner(i) = lon
               ycorner(i) = lat
            enddo
         endif
 
      endif
 
C     Check if the number of levels is too large
      if (nz.gt.nlevmax) goto 993
 
c     Allocate memory for 3d arrays: pressure, theta, pv
      allocate(pr(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array pr ***'
      allocate(prs(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array prs **'
      allocate(th(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array th ***'
      allocate(ths(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array ths **'
      allocate(pv(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array pv ***'
      allocate(pvs(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array pvs **'
 
c     Allocate memory for temporary arrays for time interpolation
      allocate(fld0(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tmp0 ***'
      allocate(fld1(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tmp1 ***'
      allocate(sfc0(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array sfc0 ***'
      allocate(sfc1(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array sfc1 ***'
 
c     ------ Pressure --------------------------------------------------
 
c     Read pressure from first data file (pfile0) on U-grid 
      call input_open (fid,pfile0)
      varname='U'
      call input_grid                                      
     >     (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >      tstart,pollon,pollat,fld0,sfc0,nz,ak,bk,stagz,timecheck)
      call input_close(fid)
 
c     Read or set pressure for second data file (pfile1)
      if ( timeshift.ne.0.) then
         call input_open (fid,pfile1)
         varname='U'
         call input_grid                                      
     >        (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >        tend,pollon,pollat,fld1,sfc1,nz,ak,bk,stagz,timecheck)
         call input_close(fid)
      else
         do i=1,nx
            do j=1,ny
               do k=1,nz
                  fld1(i,j,k) = fld0(i,j,k)
               enddo
               sfc1(i,j) = sfc0(i,j)
            enddo
         enddo
      endif
 
c     Time interpolation to get the final pressure field
      do i=1,nx
         do j=1,ny
            do k=1,nz
               pr(i,j,k) = (1.-tfrac) * fld0(i,j,k) +
     >                     tfrac      * fld1(i,j,k)
            enddo
            prs(i,j) = (1.-tfrac) * sfc0(i,j) +
     >                 tfrac      * sfc1(i,j)
         enddo
      enddo
 
c     ------ Potential temperature -------------------------------------
 
       if ( (umode.eq.'K').or.(umode.eq.'PVU') ) then
 
c         Read potential temperature from first data file <sfile0>
          call input_open (fid,sfile0)
          varname='TH'                                      ! Theta                                      
          call input_wind
     >         (fid,varname,fld0,tstart,stagz,mdv,
     >         xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
          call input_close(fid)
 
c         Read or set potential temperature for second data file (sfile1)
          if ( timeshift.ne.0.) then
             call input_open (fid,sfile1)
             varname='TH' 
             call input_wind
     >            (fid,varname,fld1,tend,stagz,mdv,
     >            xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
             call input_close(fid)
          else
             do i=1,nx
                do j=1,ny
                   do k=1,nz
                      fld1(i,j,k) = fld0(i,j,k)
                   enddo
                enddo
             enddo
          endif
 
c         Time interpolation to get the final potential temperature field
          do i=1,nx
             do j=1,ny
                do k=1,nz
                   th(i,j,k) = (1.-tfrac) * fld0(i,j,k) +
     >                         tfrac      * fld1(i,j,k)
                enddo
             enddo
          enddo
          
c         Set the surface potential temperature
          do i=1,nx                            
             do j=1,ny
                ths(i,j)=th(i,j,1)
             enddo
          enddo
          
       endif
 
 
c     ------ Potential vorticity -----------------------------------------
 
       if ( (umode.eq.'PVU') ) then
 
c         Read potential vorticity from first data file <sfile0>
          call input_open (fid,sfile0)
          varname='PV'            
          call input_wind
     >         (fid,varname,fld0,tstart,stagz,mdv,
     >         xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
          call input_close(fid)
          
c         Read or set potential vorticity for second data file (sfile1)
          if ( timeshift.ne.0.) then
             call input_open (fid,sfile1)
             varname='PV' 
             call input_wind
     >            (fid,varname,fld1,tend,stagz,mdv,
     >            xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
             call input_close(fid)
          else
             do i=1,nx
                do j=1,ny
                   do k=1,nz
                      fld1(i,j,k) = fld0(i,j,k)
                   enddo
                enddo
             enddo
          endif
 
c         Time interpolation to get the final potential vorticity field
          do i=1,nx
             do j=1,ny
                do k=1,nz
                   pv(i,j,k) = (1.-tfrac) * fld0(i,j,k) +
     >                         tfrac      * fld1(i,j,k)
                enddo
             enddo
          enddo
 
c         Set the surface potential vorticity
          do i=1,nx                          
             do j=1,ny
                pvs(i,j)=pv(i,j,1)
             enddo
          enddo
       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
      if ( rotategrid.eq.0)  then
         print*,'  pollon,pollat    : ',pollon,pollat
      else
         print*,'  pollon,pollat    : ',pollon,pollat, 'GRID ROTATION'
      endif
      print*,'  nx,ny,nz         : ',nx,ny,nz
      print*
      print*,'  Pressure loaded  : ',trim(pfile0),' ',trim(pfile1)
      if ( (umode.eq.'K').or.(umode.eq.'PVU') ) then
         print*,'  Theta loaded     : ',trim(sfile0),' ',trim(sfile1)
      endif
      if ( (umode.eq.'PVU') ) then
         print*,'  PV loaded        : ',trim(sfile0),' ',trim(sfile1)
      endif
      print*
 
c     ------------------------------------------------------------------
c     Determine the expanded list of starting coordinates
c     ------------------------------------------------------------------
 
c     Write some status information
      print*,'---- EXPAND LIST OF STARTING POSITIONS  -----------------'
      print*
 
c     ------ Read lat/lon/lev from <hfile> -----------------------------
      if ( hmode.eq.'file3' ) then
         start_n = 0
         open(10,file=hfile)
 100       continue
              start_n = start_n + 1
              read(10,*,end=101) start_lon(start_n),
     >                           start_lat(start_n),
     >                           start_lev(start_n)
              goto 100
 101       continue
           start_n = start_n - 1
         close(10)
         goto 400
      endif
 
c     ------ Get lat/lon (horizontal) coordinates ---------------------
 
c     Read lat/lon from <hfile> 
      if ( hmode.eq.'file2' ) then
         hn = 0
         open(10,file=hfile)
 200       continue
              hn = hn + 1
              read(10,*,end=201) lonlist(hn),
     >                           latlist(hn)
              goto 200
 201       continue
           hn = hn - 1
         close(10)
      endif
 
c     Get lat/lon along a line (linear in lat/lon space)
      if ( hmode.eq.'line' ) then
         do i=1,hn
            lonlist(i) = lon1 + real(i-1)/real(hn-1)*(lon2-lon1)
            latlist(i) = lat1 + real(i-1)/real(hn-1)*(lat2-lat1)
         enddo
      endif
 
c     Lat/lon box: equidistant
      if ( hmode.eq.'box.eqd' ) then
         hn  = 0
         lat = lat1
         do while ( lat.le.lat2 )      
           lon = lon1
           do while ( lon.le.lon2 ) 
             hn          = hn+1
             lonlist(hn) = lon
             latlist(hn) = lat
             lon         = lon+ds/(deltat*cos(pi180*lat))
           enddo
           lat = lat+ds/deltat
        enddo
      endif
 
c     Lat/lon box: grid
      if ( hmode.eq.'box.grid' ) then
         hn = 0
         do j=1,ny
            do i=1,nx
               rlon = xmin + real(i-1) * dx
               rlat = ymin + real(j-1) * dy
               if ( rotategrid.eq.1) then
                  lon = lmstolm(rlat,rlon,pollat,pollon)
                  lat = phstoph(rlat,rlon,pollat,pollon)
               else
                  lon = rlon
                  lat = rlat
               endif
               if ( (lon.ge.lon1).and.(lon.le.lon2).and.
     >              (lat.ge.lat1).and.(lat.le.lat2) )
     >         then
                  hn          = hn+1
                  lonlist(hn) = lon
                  latlist(hn) = lat
               endif
            enddo
         enddo
      endif
 
c     Get single starting point
      if ( hmode.eq.'point' ) then
         hn          = 1
         lonlist(hn) = lon1
         latlist(hn) = lat1
      endif
 
c     Get shifted and central starting point
      if ( hmode.eq.'shift' ) then
         hn         = 5
         lonlist(1) = lon1
         latlist(1) = lat1
         lonlist(2) = lon1+dlon
         latlist(2) = lat1
         lonlist(3) = lon1-dlon
         latlist(3) = lat1
         lonlist(4) = lon1
         latlist(4) = lat1+dlat
         lonlist(5) = lon1
         latlist(5) = lat1-dlat
      endif
 
c     Lat/lon polygon: grid
      if ( hmode.eq.'polygon.grid' ) then
 
c        Read list of polygon coordinates
         pn = 0
         open(10,file=hfile)
           read(10,*) loninpoly,latinpoly
 210       continue
              pn = pn + 1
              read(10,*,end=211) lonpoly(pn),
     >                           latpoly(pn)
 
              print*,pn,lonpoly(pn),latpoly(pn)
              
              goto 210
 211       continue
           pn = pn - 1
         close(10)
 
c        Define the polygon boundaries
         call DefSPolyBndry(latpoly,lonpoly,pn,latinpoly,loninpoly)
 
c        Get the grid points inside the polygon
         hn = 0
         do j=1,ny
            do i=1,nx
 
               rlon = xmin + real(i-1) * dx
               rlat = ymin + real(j-1) * dy
               if ( rotategrid.eq.1) then
                  lon = lmstolm(rlat,rlon,pollat,pollon)
                  lat = phstoph(rlat,rlon,pollat,pollon)
               else
                  lon = rlon
                  lat = rlat
               endif
 
               call LctPtRelBndry(lat,lon,flag)
 
               if ( (flag.eq.1).or.(flag.eq.2) ) then
                  hn          = hn+1
                  lonlist(hn) = lon
                  latlist(hn) = lat
               endif
 
            enddo
         enddo
         
      endif
 
c     Lat/lon polygon: equidistant
      if ( hmode.eq.'polygon.eqd' ) then
 
c        Read list of polygon coordinates
         pn = 0
 
         open(10,file=hfile)
           read(10,*) loninpoly,latinpoly
 220       continue
              pn = pn + 1
              read(10,*,end=221) lonpoly(pn),
     >                           latpoly(pn)
              goto 220
 221       continue
           pn = pn - 1
         close(10)
 
 
c        Define the polygon boundaries
         call DefSPolyBndry(latpoly,lonpoly,pn,latinpoly,loninpoly)
 
c        Get the grid points inside the polygon
         hn  = 0
         lat = -90.
         do while ( lat.le.90. )      
           lon = -180.
           do while ( lon.lt.180. )
 
              call LctPtRelBndry(lat,lon,flag)
               
               if ( (flag.eq.1).or.(flag.eq.2) ) then
                  hn          = hn+1
                  lonlist(hn) = lon
                  latlist(hn) = lat
 
               endif
               
               lon = lon+ds/(deltat*cos(pi180*lat))
           enddo
           lat = lat+ds/deltat
 
        enddo
 
      endif
 
c     Circle: equidistant
      if ( hmode.eq.'circle.eqd' ) then
         hn  = 0
         lat = -90.
         do while ( lat.le.90. )      
           lon = -180.
           do while ( lon.le.180. ) 
              dist = sdis(lon1,lat1,lon,lat)
              if ( dist.le.radius ) then
                 hn          = hn+1
                 lonlist(hn) = lon
                 latlist(hn) = lat
              endif
              lon = lon+ds/(deltat*cos(pi180*lat))
           enddo
           lat = lat+ds/deltat
        enddo
      endif
 
c     Circle: grid
      if ( hmode.eq.'circle.grid' ) then
         hn = 0
         do j=1,ny
            do i=1,nx
 
               rlon = xmin + real(i-1) * dx
               rlat = ymin + real(j-1) * dy
               if ( rotategrid.eq.1) then
                  lon = lmstolm(rlat,rlon,pollat,pollon)
                  lat = phstoph(rlat,rlon,pollat,pollon)
               else
                  lon = rlon
                  lat = rlat
               endif
 
               dist = sdis(lon1,lat1,lon,lat)
               if ( dist.le.radius ) then
                  hn          = hn+1
                  lonlist(hn) = lon
                  latlist(hn) = lat
               endif
            enddo
         enddo
         
      endif
 
c     Region: equidistant
      if ( hmode.eq.'region.eqd' ) then
         hn  = 0
         lat = -90.
         do while ( lat.le.90. )      
           lon = -180.
           do while ( lon.le.180. ) 
              flag = inregion(lon,lat,xcorner,ycorner)
              if ( flag.eq.1 ) then
                 hn          = hn+1
                 lonlist(hn) = lon
                 latlist(hn) = lat
              endif
              lon = lon+ds/(deltat*cos(pi180*lat))
           enddo
           lat = lat+ds/deltat
        enddo
      endif
 
c     Region: grid
      if ( hmode.eq.'region.grid' ) then
         hn = 0
         do j=1,ny
            do i=1,nx
 
               rlon  = xmin + real(i-1) * dx
               rlat  = ymin + real(j-1) * dy
               if ( rotategrid.eq.1) then
                  lon = lmstolm(rlat,rlon,pollat,pollon)
                  lat = phstoph(rlat,rlon,pollat,pollon)
               else
                  lon = rlon
                  lat = rlat
               endif
               flag = inregion(lon,lat,xcorner,ycorner)
               if ( flag.eq.1 ) then
                  hn          = hn+1
                  lonlist(hn) = lon
                  latlist(hn) = lat
               endif
            enddo
         enddo
         
      endif
 
c     ------ Get lev (vertical) coordinates -------------------------
 
c     Read level list from file
      if ( vmode.eq.'file' ) then
         vn = 0
         open(10,file=vfile)
 300       continue
              vn = vn + 1
              read(10,*,end=301) levlist(vn)
              goto 300
 301       continue
           vn = vn - 1
         close(10)
      endif
      
c     Get single starting level
      if ( vmode.eq.'level' ) then
         vn          = 1
         levlist(vn) = lev1
      endif
      
c     Get level profile
      if ( vmode.eq.'profile' ) then
         do i=1,vn
            levlist(i) = lev1 + real(i-1)/real(vn-1)*(lev2-lev1)
         enddo
      endif
 
c     ------ Compile the complete list of starting positions  ------
 
c     Get all starting points in specified vertical coordinate system
      start_n = 0
      do i=1,vn
         do j=1,hn
            start_n = start_n + 1
            start_lon(start_n) = lonlist(j)
            start_lat(start_n) = latlist(j)
            start_lev(start_n) = levlist(i)
         enddo
      enddo
 
c     ------- Rotate coordinates -----------------------------------
      do i=1,start_n
 
         lon = start_lon(i)
         lat = start_lat(i)
         
         if ( rotategrid.eq.1) then
            rlon = lmtolms(lat,lon,pollat,pollon)
            rlat = phtophs(lat,lon,pollat,pollon)
         else
            rlon = lon
            rlat = lat
         endif
         
         start_lon(i) = rlon
         start_lat(i) = rlat
 
      enddo
 
c     ------ Exit point of this section -----------------------------
 400  continue
 
c     Wriet status information
      print*,'  # expanded points : ', start_n
      print*
 
c     ------------------------------------------------------------------
c     Transform starting levels into pressure
c     ------------------------------------------------------------------
 
c     Write some status information
      print*,'---- STARTING POSITIONS ---------------------------------'
      print*
 
c     Vertical mode <hPa,asl> or simply <hPa>
      if ( (umode.eq.'hPa,asl').or.(umode.eq.'hPa') ) then
 
         do i=1,start_n
            start_pre(i) = start_lev(i)
         enddo
 
c     Vertical mode <hPa,agl>
      elseif ( umode.eq.'hPa,agl' ) then
 
         do i=1,start_n
            call get_index3(rid,rjd,rkd,start_lon(i),start_lat(i),1050.,
     >                      3,pr,prs,nx,ny,nz,xmin,ymin,dx,dy)
            tmp1 = int_index3 (prs,nx,ny,1,rid,rjd,1,mdv)
            start_pre(i) = tmp1 - start_lev(i) 
         enddo
         
c     Vertical mode <K>
      elseif ( umode.eq.'K' ) then
 
         do i=1,start_n
            call get_index3(rid,rjd,rkd,start_lon(i),start_lat(i),
     >                start_lev(i),1,th,ths,nx,ny,nz,xmin,ymin,dx,dy)
            tmp1 = int_index3 (pr,nx,ny,nz,rid,rjd,rkd,mdv)
            start_pre(i) = tmp1
         enddo
         
c     Vertical mode <PVU>
      elseif ( umode.eq.'PVU' ) then
 
         do i=1,start_n
            call get_index3(rid,rjd,rkd,start_lon(i),start_lat(i),
     >               start_lev(i),2,pv,pvs,nx,ny,nz,xmin,ymin,dx,dy)
            tmp1 = int_index3 (pr,nx,ny,nz,rid,rjd,rkd,mdv)
            start_pre(i) = tmp1
         enddo
 
      endif
 
c     Check whether the starting levels are valid (in data domain)
      do i=1,start_n
 
         call get_index3(rid,rjd,rkd,start_lon(i),start_lat(i),1050.,
     >                   3,pr,prs,nx,ny,nz,xmin,ymin,dx,dy)
         tmp1 = int_index3 (prs,nx,ny, 1,rid,rjd,real( 1),mdv)           ! Surface
         tmp2 = int_index3 (pr ,nx,ny,nz,rid,rjd,real(nz),mdv)           ! Top of domain
 
         if ( (start_pre(i).gt.tmp1).or.
     >        (start_pre(i).lt.tmp2).or.
     >        (start_lon(i).lt.xmin).or. 
     >        (start_lon(i).gt.xmax).or. 
     >        (start_lat(i).lt.ymin).or.
     >        (start_lat(i).gt.ymax) )
     >   then
            start_lon(i) = mdv
            start_lat(i) = mdv
            start_pre(i) = mdv
         endif
 
      enddo
      
c     Remove all starting points outside the domain
      i        = 1
      noutside = 0
      do while ( i.le.start_n )
         
         if ( abs(start_pre(i)-mdv).lt.eps ) then
            
            noutside = noutside +1
            if ( noutside.lt.10 ) then
                print*,'  Outside ', 
     >           start_lon(i),start_lat(i),start_lev(i)
            else if ( noutside.eq.10 ) then
                print*,'  Outside > 10 starting points'
            endif
 
            do j=i,start_n
               start_lon(j) = start_lon(j+1)
               start_lat(j) = start_lat(j+1)
               start_pre(j) = start_pre(j+1)
               start_lev(j) = start_lev(j+1)
            enddo
            start_n = start_n - 1
 
         else
         
            i = i + 1
 
         endif
         
      enddo
 
c     Write some status information
      latmin = start_lat(1)
      latmax = start_lat(1)
      lonmin = start_lon(1)
      lonmax = start_lon(1)
      premin = start_pre(1)
      premax = start_pre(1)
      do i=1,start_n
         if (start_lat(i).lt.latmin) latmin = start_lat(i)
         if (start_lat(i).gt.latmax) latmax = start_lat(i)
         if (start_lon(i).lt.lonmin) lonmin = start_lon(i)
         if (start_lon(i).gt.lonmax) lonmax = start_lon(i)
         if (start_pre(i).lt.premin) premin = start_pre(i)
         if (start_pre(i).gt.premax) premax = start_pre(i)
      enddo
      print*,'  min(lat),max(lat) : ', latmin,latmax
      print*,'  min(lon),max(lon) : ', lonmin,lonmax
      print*,'  min(pre),max(pre) : ', premin,premax
      print*
      print*,'  # starting points : ', start_n
      print*
 
      
c     ------------------------------------------------------------------
c     Write starting positions to output file
c     ------------------------------------------------------------------
 
c     Rotate coordinates
      do i=1,start_n
 
         rlon = start_lon(i)
         rlat = start_lat(i)
         
         if ( rotategrid.eq.1) then
            lon = lmstolm(rlat,rlon,pollat,pollon)
            lat = phstoph(rlat,rlon,pollat,pollon)
         else
            lon = rlon
            lat = rlat
         endif
         
         start_lon(i) = lon
         start_lat(i) = lat
 
      enddo
 
c     Output as a trajectory file (with only one time == 0)
      if (oformat.ne.-1) then
 
         allocate(tra(start_n,1,5),stat=stat)
 
         vars(1)  ='time'
         vars(2)  ='lon'
         vars(3)  ='lat'
         vars(4)  ='p'
         vars(5)  ='level'
         call wopen_tra(fid,ofile,start_n,1,5,reftime,vars,oformat)
 
         do i=1,start_n
            tra(i,1,1) = 0.
            tra(i,1,2) = start_lon(i)
            tra(i,1,3) = start_lat(i)
            tra(i,1,4) = start_pre(i)
            tra(i,1,5) = start_lev(i)
         enddo
         call write_tra(fid,tra,start_n,1,5,oformat)
         
         call close_tra(fid,oformat)
 
c     Output as a triple list (corresponding to <startf> file)
      else
         
         fid = 10
         open(fid,file=ofile)
          do i=1,start_n
             write(fid,'(3f10.3)') start_lon(i),start_lat(i),
     >                             start_pre(i) 
          enddo
         close(fid)
         
      endif
 
c     Write some status information, and end of program message
      print*  
      print*,'---- STATUS INFORMATION --------------------------------'
      print*
      print*,'ok'
      print*
      print*,'       *** END OF PROGRAM CREATE_STARTF ***'
      print*,'========================================================='
        
c     ------------------------------------------------------------------
c     Exception handling
c     ------------------------------------------------------------------
 
      stop
 
 993  write(*,*) '*** ERROR: problems with array size'
      call exit(1)
 
      end
 
c     --------------------------------------------------------------------------
c     Split a region string and get corners of the domain
c     --------------------------------------------------------------------------
 
      subroutine regionsplit(string,iregion,xcorner,ycorner)
 
c     The region string comes either as <lonw,lone,lats,latn> or as <lon1,lat1,
c     lon2,lat2,lon3,lat3,lon4,lat4>: split it into ints components and get the
c     four coordinates for the region
      
      implicit none
 
c     Declaration of subroutine parameters
      character*80    string
      real            xcorner(4),ycorner(4)
      integer         iregion
 
c     Local variables
      integer         i,n
      integer         il,ir
      real            subfloat (80)
      integer         stat
      integer         len
 
c     ------- Split the string
      i    = 1
      n    = 0
      stat = 0
      il   = 1
      len  = len_trim(string)
 
 100  continue
 
c     Find start of a substring
      do while ( stat.eq.0 )
         if ( string(i:i).ne.' ' ) then
            stat = 1
            il   = i
         else
            i = i + 1
         endif
      enddo
 
c     Find end of substring
      do while ( stat.eq.1 )         
         if ( ( string(i:i).eq.' ' ) .or. ( i.eq.len ) ) then
            stat = 2
            ir   = i
         else
            i    = i + 1
         endif
      enddo
 
c     Convert the substring into a number
      if ( stat.eq.2 ) then
         n = n + 1
         read(string(il:ir),*) subfloat(n)
         stat = 0
      endif
 
      if ( i.lt.len ) goto 100
 
 
c     -------- Get the region number
      
      iregion = nint(subfloat(1))
 
c     -------- Get the corners of the region
      
      if ( n.eq.5 ) then     ! lonw(2),lone(3),lats(4),latn(5)
 
         xcorner(1) = subfloat(2)
         ycorner(1) = subfloat(4)
 
         xcorner(2) = subfloat(3)
         ycorner(2) = subfloat(4)
 
         xcorner(3) = subfloat(3)
         ycorner(3) = subfloat(5)
         
         xcorner(4) = subfloat(2)
         ycorner(4) = subfloat(5)
        
      elseif ( n.eq.9 ) then     ! lon1,lat1,lon2,lat2,lon3,lon4,lat4
 
         xcorner(1) = subfloat(2)
         ycorner(1) = subfloat(3)
 
         xcorner(2) = subfloat(4)
         ycorner(2) = subfloat(5)
 
         xcorner(3) = subfloat(6)
         ycorner(3) = subfloat(7)
         
         xcorner(4) = subfloat(8)
         ycorner(4) = subfloat(9)
 
      else
         
         print*,' ERROR: invalid region specification '
         print*,'     ',trim(string)
         stop
         
      endif
         
 
      end
 
c     --------------------------------------------------------------------------
c     Decide whether lat/lon point is in or out of region
c     --------------------------------------------------------------------------
      
      integer function inregion (lon,lat,xcorner,ycorner)
      
c     Decide whether point (lon/lat) is in the region specified by <xcorner(1..4),
c     ycorner(1..4).
      
      implicit none
      
c     Declaration of subroutine parameters
      real    lon,lat
      real    xcorner(4),ycorner(4)
 
c     Local variables
      integer flag
      real    xmin,xmax,ymin,ymax
      integer i
 
c     Reset the flag
      flag = 0
 
c     Set some boundaries
      xmax = xcorner(1)
      xmin = xcorner(1)
      ymax = ycorner(1)
      ymin = ycorner(1)
      do i=2,4
        if (xcorner(i).lt.xmin) xmin = xcorner(i)
        if (xcorner(i).gt.xmax) xmax = xcorner(i)
        if (ycorner(i).lt.ymin) ymin = ycorner(i)
        if (ycorner(i).gt.ymax) ymax = ycorner(i)
      enddo
 
c     Do the tests - set flag=1 if all tests pased
      if (lon.lt.xmin) goto 970
      if (lon.gt.xmax) goto 970
      if (lat.lt.ymin) goto 970
      if (lat.gt.ymax) goto 970
      
      if ((lon-xcorner(1))*(ycorner(2)-ycorner(1))-
     >    (lat-ycorner(1))*(xcorner(2)-xcorner(1)).gt.0.) goto 970
      if ((lon-xcorner(2))*(ycorner(3)-ycorner(2))-
     >    (lat-ycorner(2))*(xcorner(3)-xcorner(2)).gt.0.) goto 970
      if ((lon-xcorner(3))*(ycorner(4)-ycorner(3))-
     >    (lat-ycorner(3))*(xcorner(4)-xcorner(3)).gt.0.) goto 970
      if ((lon-xcorner(4))*(ycorner(1)-ycorner(4))-
     >    (lat-ycorner(4))*(xcorner(1)-xcorner(4)).gt.0.) goto 970
 
      flag = 1
 
c     Return the value
 970  continue
      
      inregion = flag
      
      return
      
      end
 
c     --------------------------------------------------------------------------
c     Spherical distance between lat/lon points                                                       
c     --------------------------------------------------------------------------
 
      real function sdis(xp,yp,xq,yq)
c
c     calculates spherical distance (in km) between two points given
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
c
      real      re
      parameter (re=6370.)
      real      pi180
      parameter (pi180=3.14159/180.)
      real      xp,yp,xq,yq,arg
 
      arg=sin(pi180*yp)*sin(pi180*yq)+
     >    cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(xp-xq))
      if (arg.lt.-1.) arg=-1.
      if (arg.gt.1.) arg=1.
 
      sdis=re*acos(arg)
 
      end
 
 
c     ****************************************************************
c     * Given some spherical polygon S and some point X known to be  *
c     * located inside S, these routines will determine if an arbit- *
c     * -rary point P lies inside S, outside S, or on its boundary.  *
c     * The calling program must first call DefSPolyBndry to define  *
c     * the boundary of S and the point X. Any subsequent call to    *
c     * subroutine LctPtRelBndry will determine if some point P lies *
c     * inside or outside S, or on its boundary. (Usually            *
c     * DefSPolyBndry is called once, then LctPrRelBndry is called   *
c     * many times).                                                 *
c     *                                                              * 
c     * REFERENCE:            Bevis, M. and Chatelain, J.-L. (1989)  * 
c     *                       Maflaematical Geology, vol 21.         *
c     * VERSION 1.0                                                  *
c     ****************************************************************
 
      Subroutine DefSPolyBndry(vlat,vlon,nv,xlat, xlon)
 
c     ****************************************************************
c     * This mmn entry point is used m define ~e spheric~ polygon S  *
c     * and the point X.                                             *
c     * ARGUMENTS:                                                   *
c     * vlat,vlon (sent) ... vectors containing the latitude and     * 
c     *                      longitude of each vertex of the         *
c     *                      spherical polygon S. The ith.vertex is  *
c     *                      located at [vlat(i),vlon(i)].           *
c     * nv        (sent) ... the number of vertices and sides in the *
c     *                      spherical polygon S                     *
c     * xlat,xlon (sent) ... latitude and longitude of some point X  *
c     *                      located inside S. X must not be located *
c     *                      on any great circle that includes two   *
c     *                      vertices of S.                          *
c     *                                                              *
c     * UNITS AND SIGN CONVENTION:                                   *
c     *  Latitudes and longitudes are specified in degrees.          *
c     *  Latitudes are positive to the north and negative to the     *
c     *  south.                                                      *
c     *  Longitudes are positive to the east and negative to the     *
c     *  west.                                                       *
c     *                                                              * 
c     * VERTEX ENUMERATION:                                          * 
c     * The vertices of S should be numbered sequentially around the *
c     * border of the spherical polygon. Vertex 1 lies between vertex*
c     * nv and vertex 2. Neighbouring vertices must be seperated by  *
c     * less than 180 degrees. (In order to generate a polygon side  *
c     * whose arc length equals or exceeds 180 degrees simply        *
c     * introduce an additional (pseudo)vertex). Having chosen       *
c     * vertex 1, the user may number the remaining vertices in      *
c     * either direction. However if the user wishes to use the      *
c     * subroutine SPA to determine the area of the polygon S (Bevis *
c     * & Cambareri, 1987, Math. Geol., v.19, p. 335-346) then he or *
c     * she must follow the convention whereby in moving around the  *
c     * polygon border in the direction of increasing vertex number  *
c     * clockwise bends occur at salient vertices. A vertex is       *
c     * salient if the interior angle is less than 180 degrees.      *
c     * (In the case of a convex polygon this convention implies     *
c     * that vertices are numbered in clockwise sequence).           *
c     ****************************************************************
 
      implicit none
      
      integer mxnv,nv
 
c     ----------------------------------------------------------------
c     Edit next statement to increase maximum number of vertices that 
c     may be used to define the spherical polygon S               
c     The value of parameter mxnv in subroutine LctPtRelBndry must match
c     that of parameter mxnv in this subroutine, as assigned above.
c     ----------------------------------------------------------------
      parameter (mxnv=500)
 
      real  vlat(nv),vlon(nv),xlat,xlon,dellon
      real  tlonv(mxnv),vlat_c(mxnv),vlon_c(mxnv),xlat_c,xlon_c
      integer i,ibndry,nv_c,ip
 
      data ibndry/0/
      
      common/spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
 
      if (nv.gt.mxnv) then
         print *,'nv exceeds maximum allowed value'
         print *,'adjust parameter mxnv in subroutine DefSPolyBndry'
         stop
      endif
 
      ibndry=1                  ! boundary defined at least once (flag)
      nv_c=nv                   ! copy for named common
      xlat_c=xlat               ! . . . .
      xlon_c=xlon               !
 
      do i=1,nv
         vlat_c(i)=vlat(i)      ! "
         vlon_c(i)=vlon(i)      !
 
         call TrnsfmLon(xlat,xlon,vlat(i),vlon(i),tlonv(i))
 
         if (i.gt.1) then
            ip=i-1
         else
            ip=nv
         endif
         
         if ((vlat(i).eq.vlat(ip)).and.(vlon(i).eq.vlon(ip))) then
            print *,'DefSPolyBndry detects user error:'
            print *,'vertices ',i,' and ',ip,' are not distinct'
            print*,'lat ',i,ip,vlat(i),vlat(ip)
            print*,'lon ',i,ip,vlon(i),vlon(ip)            
            stop
         endif
 
         if (tlonv(i).eq.tlonv(ip)) then
            print *,'DefSPolyBndry detects user error:'
            print *,'vertices ',i,' & ',ip,' on same gt. circle as X'
            stop
         endif
 
         if (vlat(i).eq.(-vlat(ip))) then
            dellon=vlon(i)-vlon(ip)
            if (dellon.gt.+180.) dellon=dellon-360.
            if (dellon.lt.-180.) dellon=dellon-360.
            if ((dellon.eq.+180.0).or.(dellon.eq.-180.0)) then
               print *,'DefSPolyBndry detects user error:'
               print *,'vertices ',i,' and ',ip,' are antipodal'
               stop
            endif
         endif
      enddo
 
      return
      
      end
 
 
c     ****************************************************************
 
      Subroutine LctPtRelBndry(plat,plon,location)
 
c     ****************************************************************
 
c     ****************************************************************
c     * This routine is used to see if some point P is located       *
c     * inside, outside or on the boundary of the spherical polygon  *
c     * S previously defined by a call to subroutine DefSPolyBndry.  *
c     * There is a single restriction on point P: it must not be     *
c     * antipodal to the point X defined in the call to DefSPolyBndry*
c     * (ie.P and X cannot be seperated by exactly 180 degrees).     *
c     * ARGUMENTS:                                                   *  
c     * plat,plon (sent)... the latitude and longitude of point P    *
c     * location (returned)... specifies the location of P:          *
c     *                        location=0 implies P is outside of S  *
c     *                        location=1 implies P is inside of S   *
c     *                        location=2 implies P on boundary of S *
c     *                        location=3 implies user error (P is   *
c     *                                     antipodal to X)          *
c     * UNFfS AND SIGN CONVENTION:                                   * 
c     *  Latitudes and longitudes are specified in degrees.          *
c     *  Latitudes are positive to the north and negative to the     *
c     *  south.                                                      *    
c     *  Longitudes are positive to the east and negative to the     *
c     *  west.                                                       *
c     ****************************************************************
      
      implicit none
      
      integer mxnv
 
c     ----------------------------------------------------------------
c     The statement below must match that in subroutine DefSPolyBndry
c     ----------------------------------------------------------------
 
      parameter (mxnv=500)
 
      real tlonv(mxnv),vlat_c(mxnv),vlon_c(mxnv),xlat_c,xlon_c
      real plat,plon,vAlat,vAlon,vBlat,vBlon,tlonA,tlonB,tlonP
      real tlon_X,tlon_P,tlon_B,dellon
      integer i,ibndry,nv_c,location,icross,ibrngAB,ibrngAP,ibrngPB
      integer ibrng_BX,ibrng_BP,istrike
 
      common/spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
 
      if (ibndry.eq.0) then     ! user has never defined the bndry
         print*,'Subroutine LctPtRelBndry detects user error:'
         print*,'Subroutine DefSPolyBndry must be called before'
         print*,'subroutine LctPtRelBndry can be called'
         stop
      endif
 
      if (plat.eq.(-xlat_c)) then
         dellon=plon-xlon_c
         if (dellon.lt.(-180.)) dellon=dellon+360.
         if (dellon.gt.+180.) dellon=dellon-360.
         if ((dellon.eq.+180.0).or.(dellon.eq.-180.)) then
            print*,'Warning: LctPtRelBndry detects case P antipodal
     >           to X'
            print*,'location of P relative to S is undetermined'
            location=3
            return
         endif
      endif 
 
      location=0                ! default ( P is outside S)
      icross=0                  ! initialize counter
 
      if ((plat.eq.xlat_c).and.(plon.eq.xlon_c)) then
         location=1
         return
      endif
 
      
      call TrnsfmLon (xlat_c,xlon_c,plat,plon,tlonP)
 
      do i=1,nv_c              ! start of loop over sides of S 
 
         vAlat=vlat_c(i)
         vAlon=vlon_c(i)
         tlonA=tlonv(i)
 
         if (i.lt.nv_c) then
            vBlat=vlat_c(i+1)
            vBlon=vlon_c(i+1)
            tlonB=tlonv(i+1)
         else
            vBlat=vlat_c(1)
            vBlon=vlon_c(1)
            tlonB=tlonv(1)
         endif
         
         istrike=0
         
         if (tlonP.eq.tlonA) then
            istrike=1
         else
            call EastOrWest(tlonA,tlonB,ibrngAB)
            call EastOrWest(tlonA,tlonP,ibrngAP)
            call EastOrWest(tlonP,tlonB,ibrngPB)
            
 
            if((ibrngAP.eq.ibrngAB).and.(ibrngPB.eq.ibrngAB)) istrike=1
         endif
 
         
         if (istrike.eq.1) then
 
            if ((plat.eq.vAlat).and.(plon.eq.vAlon)) then
               location=2       ! P lies on a vertex of S
               return
            endif
            call TrnsfmLon(vAlat,vAlon,xlat_c,xlon_c,tlon_X)
            call TrnsfmLon(vAlat,vAlon,vBlat,vBlon,tlon_B)
            call TrnsfmLon(vAlat,vAlon,plat,plon,tlon_P)
            
            if (tlon_P.eq.tlon_B) then
               location=2       ! P lies on side of S
               return 
            else
               call EastOrWest(tlon_B,tlon_X,ibrng_BX)
               call EastOrWest(tlon_B,tlon_P,ibrng_BP)
               if(ibrng_BX.eq.(-ibrng_BP)) icross=icross+1
            endif
            
         endif
      enddo                     ! end of loop over the sides of S
 
 
c     if the arc XP crosses the boundary S an even number of times then P
c     is in S
 
      if (mod(icross,2).eq.0) location=1
 
      return
 
      end
 
 
c     ****************************************************************
      
      subroutine TrnsfmLon(plat,plon,qlat,qlon,tranlon)
 
c     ****************************************************************
c     * This subroutine is required by subroutines DefSPolyBndry &   *
c     * LctPtRelBndry. It finds the 'longitude' of point Q in a      *
c     * geographic coordinate system for which point P acts as a     *
c     * 'north pole'. SENT: plat,plon,qlat,qlon, in degrees.         *
c     * RETURNED: tranlon, in degrees.                               *
c     ****************************************************************
 
      implicit none
 
      real pi,dtr,plat,plon,qlat,qlon,tranlon,t,b
      parameter (pi=3.141592654,dtr=pi/180.0)
 
      if (plat.eq.90.) then
         tranlon=qlon
      else
         t=sin((qlon-plon)*dtr)*cos(qlat*dtr)
         b=sin(dtr*qlat)*cos(plat*dtr)-cos(qlat*dtr)*sin(plat*dtr)
     >    *cos((qlon-plon)*dtr)
         tranlon=atan2(t,b)/dtr
      endif
 
      return
      end
 
c     ****************************************************************
 
      subroutine EastOrWest(clon,dlon,ibrng)
 
c     ****************************************************************
c     * This subroutine is required by subroutine LctPtRelBndry.     *
c     * This routine determines if in travelling the shortest path   *
c     * from point C (at longitude clon) to point D (at longitude    *
c     * dlon) one is heading east, west or neither.                  *
c     * SENT: clon,dlon; in degrees. RETURNED: ibrng                 *
c     * (1=east,-1=west, 0=neither).                                 *
c     ****************************************************************
 
      implicit none
      real clon,dlon,del
      integer ibrng
      del=dlon-clon
      if (del.gt.180.) del=del-360.
      if (del.lt.-180.) del=del+360.
      if ((del.gt.0.0).and.(del.ne.180.)) then
         ibrng=-1               ! (D is west of C)
      elseif ((del.lt.0.0).and.(del.ne.-180.)) then
         ibrng=+1               ! (D is east of C)
      else
         ibrng=0                ! (D north or south of C)
      endif
      return
      end
 
 

Rev 7	Rev 16
1	`PROGRAM startf`	1	`PROGRAM startf`
2		2
3	`c **************************************************************`	3	`c **************************************************************`
4	`c * Create a <startfile> for <lagrangto>. It can be chosen *`	4	`c * Create a <startfile> for <lagrangto>. It can be chosen *`
5	`c * whether to start from an isentropic or an isobaric *`	5	`c * whether to start from an isentropic or an isobaric *`
6	`c * surface. The starting points are equidistantly distributed *`	6	`c * surface. The starting points are equidistantly distributed *`
7	`c * Michael Sprenger / Autumn 2004 *`	7	`c * Michael Sprenger / Autumn 2004 *`
8	`c **************************************************************`	8	`c **************************************************************`
9		9
10	`implicit none`	10	`implicit none`
11		11
12		12
13	`c --------------------------------------------------------------`	13	`c --------------------------------------------------------------`
14	`c Set parameters`	14	`c Set parameters`
15	`c --------------------------------------------------------------`	15	`c --------------------------------------------------------------`
16		16
17	`c Maximum number of starting positions`	17	`c Maximum number of starting positions`
18	`integer nmax`	18	`integer nmax`
19	`parameter (nmax=1000000)`	19	`parameter (nmax=1000000)`
20		20
21	`c Maximum number of model levels`	21	`c Maximum number of model levels`
22	`integer nlevmax`	22	`integer nlevmax`
23	`parameter (nlevmax=200)`	23	`parameter (nlevmax=200)`
24		24
25	`c Grid constant (distance in km corresponding to 1 deg at the equator)`	25	`c Grid constant (distance in km corresponding to 1 deg at the equator)`
26	`real deltat`	26	`real deltat`
27	`parameter (deltat=111.)`	27	`parameter (deltat=111.)`
28		28
29	`c Mathematical constant (conversion degree -> radian)`	29	`c Mathematical constant (conversion degree -> radian)`
30	`real pi180`	30	`real pi180`
31	`parameter (pi180=3.14159/180.)`	31	`parameter (pi180=3.14159/180.)`
32		32
33	`c Numerical epsilon`	33	`c Numerical epsilon`
34	`real eps`	34	`real eps`
35	`parameter (eps=0.00001)`	35	`parameter (eps=0.00001)`
36		36
37	`c --------------------------------------------------------------`	37	`c --------------------------------------------------------------`
38	`c Set variables`	38	`c Set variables`
39	`c --------------------------------------------------------------`	39	`c --------------------------------------------------------------`
40		40
41	`c Filenames and output format`	41	`c Filenames and output format`
42	`character*80 pfile0,pfile1 ! P filenames`	42	`character*80 pfile0,pfile1 ! P filenames`
43	`character*80 sfile0,sfile1 ! S filenames`	43	`character*80 sfile0,sfile1 ! S filenames`
44	`character*80 ofile ! Output filename`	44	`character*80 ofile ! Output filename`
45	`integer oformat ! Output format`	45	`integer oformat ! Output format`
46	`real timeshift ! Time shift relative to data files <*0>`	46	`real timeshift ! Time shift relative to data files <*0>`
47	`real timeinc ! Time increment between input files`	47	`real timeinc ! Time increment between input files`
48		48
49	`c Horizontal grid`	49	`c Horizontal grid`
50	`character*80 hmode ! Horizontale mode`	50	`character*80 hmode ! Horizontale mode`
51	`real lat1,lat2,lon1,lon2 ! Lat/lon boundaries`	51	`real lat1,lat2,lon1,lon2 ! Lat/lon boundaries`
52	`real ds,dlon,dlat ! Distance and lat/lon shifts`	52	`real ds,dlon,dlat ! Distance and lat/lon shifts`
53	`character*80 hfile ! Filename`	53	`character*80 hfile ! Filename`
54	`integer hn ! Number of entries in lat/lon list`	54	`integer hn ! Number of entries in lat/lon list`
55	`real latlist(nmax) ! List of latitudes`	55	`real latlist(nmax) ! List of latitudes`
56	`real lonlist(nmax) ! List of longitudes`	56	`real lonlist(nmax) ! List of longitudes`
57	`integer pn ! Number of entries in lat/lon poly`	57	`integer pn ! Number of entries in lat/lon poly`
58	`real latpoly(500) ! List of polygon latitudes`	58	`real latpoly(500) ! List of polygon latitudes`
59	`real lonpoly(500) ! List of polygon longitudes`	59	`real lonpoly(500) ! List of polygon longitudes`
60	`real loninpoly,latinpoly ! Lon/lat inside polygon`	60	`real loninpoly,latinpoly ! Lon/lat inside polygon`
61	`character*80 regionf ! Region file`	61	`character*80 regionf ! Region file`
62	`integer iregion ! Region number`	62	`integer iregion ! Region number`
63	`real xcorner(4),ycorner(4) ! Vertices of region`	63	`real xcorner(4),ycorner(4) ! Vertices of region`
64		64
65	`c Vertical grid`	65	`c Vertical grid`
66	`character*80 vmode ! Vertical mode`	66	`character*80 vmode ! Vertical mode`
67	`real lev1,lev2,levlist(nmax) ! Single levels, and list of levels`	67	`real lev1,lev2,levlist(nmax) ! Single levels, and list of levels`
68	`character*80 vfile ! Filename`	68	`character*80 vfile ! Filename`
69	`integer vn ! Number of entries`	69	`integer vn ! Number of entries`
70		70
71	`c Unit of vertical axis`	71	`c Unit of vertical axis`
72	`character*80 umode ! Unit of vertical axis`	72	`character*80 umode ! Unit of vertical axis`
73		73
74	`c Flag for 'no time check'`	74	`c Flag for 'no time check'`
75	`character*80 timecheck ! Either 'no' or 'yes'`	75	`character*80 timecheck ! Either 'no' or 'yes'`
76		76
77	`c List of all starting positions`	77	`c List of all starting positions`
78	`integer start_n ! Number of coordinates`	78	`integer start_n ! Number of coordinates`
79	`real start_lat(nmax) ! Latitudes`	79	`real start_lat(nmax) ! Latitudes`
80	`real start_lon(nmax) ! Longitudes`	80	`real start_lon(nmax) ! Longitudes`
81	`real start_lev(nmax) ! Levels (depending on vertical unit)`	81	`real start_lev(nmax) ! Levels (depending on vertical unit)`
82	`real start_pre(nmax) ! Level in hPa`	82	`real start_pre(nmax) ! Level in hPa`
83	`integer reftime(6) ! Reference time`	83	`integer reftime(6) ! Reference time`
84	`character*80 vars(10) ! Name of output fields (time,lon,lat,p)`	84	`character*80 vars(10) ! Name of output fields (time,lon,lat,p)`
85	`real,allocatable, dimension (:,:,:) :: tra ! Trajectories (ntra,ntim,ncol)`	85	`real,allocatable, dimension (:,:,:) :: tra ! Trajectories (ntra,ntim,ncol)`
86	`real latmin,latmax`	86	`real latmin,latmax`
87	`real lonmin,lonmax`	87	`real lonmin,lonmax`
88	`real premin,premax`	88	`real premin,premax`
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	`real,allocatable, dimension (:,:,:) :: pr ! 3d pressure`	98	`real,allocatable, dimension (:,:,:) :: pr ! 3d pressure`
99	`real,allocatable, dimension (:,:) :: prs ! surface pressure`	99	`real,allocatable, dimension (:,:) :: prs ! surface pressure`
100	`real,allocatable, dimension (:,:,:) :: th ! 3d potential temperature`	100	`real,allocatable, dimension (:,:,:) :: th ! 3d potential temperature`
101	`real,allocatable, dimension (:,:) :: ths ! surface poential temperature`	101	`real,allocatable, dimension (:,:) :: ths ! surface poential temperature`
102	`real,allocatable, dimension (:,:,:) :: pv ! 3d potential vorticity`	102	`real,allocatable, dimension (:,:,:) :: pv ! 3d potential vorticity`
103	`real,allocatable, dimension (:,:) :: pvs ! surface potential vorticiy`	103	`real,allocatable, dimension (:,:) :: pvs ! surface potential vorticiy`
104	`character*80 varname ! Name of input variable`	104	`character*80 varname ! Name of input variable`
105	`integer fid ! File identifier`	105	`integer fid ! File identifier`
106	`real stagz ! Vertical staggering`	106	`real stagz ! Vertical staggering`
107	`real mdv ! Missing data values`	107	`real mdv ! Missing data values`
108	`real tstart,tend ! Time on P and S file`	108	`real tstart,tend ! Time on P and S file`
109	`real rid,rjd,rkd ! Real grid position`	109	`real rid,rjd,rkd ! Real grid position`
110	`integer rotategrid ! Flag whether to rotate grid`	110	`integer rotategrid ! Flag whether to rotate grid`
111		111
112	`c Auxiliary variable`	112	`c Auxiliary variable`
113	`integer i,j,k`	113	`integer i,j,k`
114	`real lon,lat,rlon,rlat`	114	`real lon,lat,rlon,rlat`
115	`real rd`	115	`real rd`
116	`integer stat,flag`	116	`integer stat,flag`
117	`real tmp1,tmp2`	117	`real tmp1,tmp2`
118	`real tfrac,frac`	118	`real tfrac,frac`
119	`real radius,dist`	119	`real radius,dist`
120	`character*80 string`	120	`character*80 string`
121	`character*80 selectstr`	121	`character*80 selectstr`
-		122	`integer noutside`
122	`real,allocatable, dimension (:,:,:) :: fld0`	123	`real,allocatable, dimension (:,:,:) :: fld0`
123	`real,allocatable, dimension (:,:,:) :: fld1`	124	`real,allocatable, dimension (:,:,:) :: fld1`
124	`real,allocatable, dimension (:,: ) :: sfc0`	125	`real,allocatable, dimension (:,: ) :: sfc0`
125	`real,allocatable, dimension (:,:) :: sfc1`	126	`real,allocatable, dimension (:,:) :: sfc1`
126		127
127	`c Externals`	128	`c Externals`
128	`real int_index3 ! 3d interpolation`	129	`real int_index3 ! 3d interpolation`
129	`external int_index3`	130	`external int_index3`
130	`real sdis ! Speherical distance`	131	`real sdis ! Speherical distance`
131	`external sdis`	132	`external sdis`
132	`integer inregion ! In/out of region`	133	`integer inregion ! In/out of region`
133	`external inrehion`	134	`external inrehion`
134	`real lmtolms ! Grid rotation`	135	`real lmtolms ! Grid rotation`
135	`real phtophs`	136	`real phtophs`
136	`real lmstolm`	137	`real lmstolm`
137	`real phstoph`	138	`real phstoph`
138	`external lmtolms,phtophs,lmstolm,phstoph`	139	`external lmtolms,phtophs,lmstolm,phstoph`
139		140
140	`c ------------------------------------------------------------------`	141	`c ------------------------------------------------------------------`
141	`c Start of program, Read parameters`	142	`c Start of program, Read parameters`
142	`c ------------------------------------------------------------------`	143	`c ------------------------------------------------------------------`
143		144
144	`c Write start message`	145	`c Write start message`
145	`print*,'========================================================='`	146	`print*,'========================================================='`
146	`print,' START OF PROGRAM CREATE_STARTF *'`	147	`print,' START OF PROGRAM CREATE_STARTF *'`
147	`print*`	148	`print*`
148		149
149	`c Read parameter file`	150	`c Read parameter file`
150	`open(10,file='startf.param')`	151	`open(10,file='startf.param')`
151		152
152	`c Input P and S file`	153	`c Input P and S file`
153	`read(10,*) pfile0,pfile1`	154	`read(10,*) pfile0,pfile1`
154	`read(10,*) sfile0,sfile1`	155	`read(10,*) sfile0,sfile1`
155	`read(10,*) ofile`	156	`read(10,*) ofile`
156		157
157	`c Read name of region file`	158	`c Read name of region file`
158	`read(10,*) regionf`	159	`read(10,*) regionf`
159		160
160	`c Reference time`	161	`c Reference time`
161	`do i=1,6`	162	`do i=1,6`
162	`read(10,*) reftime(i)`	163	`read(10,*) reftime(i)`
163	`enddo`	164	`enddo`
164		165
165	`c Time shift relative to data files <pfile0,sfile0> - format (hh.mm)`	166	`c Time shift relative to data files <pfile0,sfile0> - format (hh.mm)`
166	`read(10,*) timeshift`	167	`read(10,*) timeshift`
167		168
168	`c Read timeincrement between input files`	169	`c Read timeincrement between input files`
169	`read(10,*) timeinc`	170	`read(10,*) timeinc`
170		171
171	`c Parameters for horizontal grid`	172	`c Parameters for horizontal grid`
172	`read(10,*) hmode`	173	`read(10,*) hmode`
173	`if ( hmode.eq.'file' ) then ! from file`	174	`if ( hmode.eq.'file' ) then ! from file`
174	`read(10,*) hfile`	175	`read(10,*) hfile`
175	`elseif ( hmode.eq.'line' ) then ! along a line`	176	`elseif ( hmode.eq.'line' ) then ! along a line`
176	`read(10,*) lon1,lon2,lat1,lat2,hn`	177	`read(10,*) lon1,lon2,lat1,lat2,hn`
177	`elseif ( hmode.eq.'box.eqd' ) then ! box: 2d equidistant`	178	`elseif ( hmode.eq.'box.eqd' ) then ! box: 2d equidistant`
178	`read(10,*) lon1,lon2,lat1,lat2,ds`	179	`read(10,*) lon1,lon2,lat1,lat2,ds`
179	`elseif ( hmode.eq.'box.grid' ) then ! box: 2d grid`	180	`elseif ( hmode.eq.'box.grid' ) then ! box: 2d grid`
180	`read(10,*) lon1,lon2,lat1,lat2`	181	`read(10,*) lon1,lon2,lat1,lat2`
181	`elseif ( hmode.eq.'point' ) then ! single point`	182	`elseif ( hmode.eq.'point' ) then ! single point`
182	`read(10,*) lon1,lat1`	183	`read(10,*) lon1,lat1`
183	`elseif ( hmode.eq.'shift' ) then ! centre + shifted`	184	`elseif ( hmode.eq.'shift' ) then ! centre + shifted`
184	`read(10,*) lon1,lat1,dlon,dlat`	185	`read(10,*) lon1,lat1,dlon,dlat`
185	`elseif ( hmode.eq.'polygon.eqd' ) then ! polygon: 2d equidistant`	186	`elseif ( hmode.eq.'polygon.eqd' ) then ! polygon: 2d equidistant`
186	`read(10,*) hfile,ds`	187	`read(10,*) hfile,ds`
187	`elseif ( hmode.eq.'polygon.grid' ) then ! polygon: 2d grid`	188	`elseif ( hmode.eq.'polygon.grid' ) then ! polygon: 2d grid`
188	`read(10,*) hfile`	189	`read(10,*) hfile`
189	`elseif ( hmode.eq.'circle.eqd' ) then ! circle: 2d equidistant`	190	`elseif ( hmode.eq.'circle.eqd' ) then ! circle: 2d equidistant`
190	`read(10,*) lon1,lat1,radius,ds`	191	`read(10,*) lon1,lat1,radius,ds`
191	`elseif ( hmode.eq.'circle.grid' ) then ! circle: 2d grid`	192	`elseif ( hmode.eq.'circle.grid' ) then ! circle: 2d grid`
192	`read(10,*) lon1,lat1,radius`	193	`read(10,*) lon1,lat1,radius`
193	`elseif ( hmode.eq.'region.eqd' ) then ! region: 2d equidistant`	194	`elseif ( hmode.eq.'region.eqd' ) then ! region: 2d equidistant`
194	`read(10,*) iregion,ds`	195	`read(10,*) iregion,ds`
195	`elseif ( hmode.eq.'region.grid' ) then ! iregion: 2d grid`	196	`elseif ( hmode.eq.'region.grid' ) then ! iregion: 2d grid`
196	`read(10,*) iregion`	197	`read(10,*) iregion`
197	`else`	198	`else`
198	`print*,' ERROR: horizontal mode not supported ',trim(hmode)`	199	`print*,' ERROR: horizontal mode not supported ',trim(hmode)`
199	`stop`	200	`stop`
200	`endif`	201	`endif`
201		202
202	`c Parameters for vertical grid`	203	`c Parameters for vertical grid`
203	`read(10,*) vmode`	204	`read(10,*) vmode`
204	`if ( vmode.eq.'file') then ! from file`	205	`if ( vmode.eq.'file') then ! from file`
205	`read(10,*) vfile`	206	`read(10,*) vfile`
206	`elseif ( vmode.eq.'level' ) then ! single level (explicit command)`	207	`elseif ( vmode.eq.'level' ) then ! single level (explicit command)`
207	`read(10,*) lev1`	208	`read(10,*) lev1`
208	`elseif ( vmode.eq.'list') then ! a list`	209	`elseif ( vmode.eq.'list') then ! a list`
209	`read(10,*) vn`	210	`read(10,*) vn`
210	`read(10,*) (levlist(i),i=1,vn)`	211	`read(10,*) (levlist(i),i=1,vn)`
211	`elseif ( vmode.eq.'profile') then ! a profile`	212	`elseif ( vmode.eq.'profile') then ! a profile`
212	`read(10,*) lev1,lev2,vn`	213	`read(10,*) lev1,lev2,vn`
213	`elseif ( vmode.eq.'nil') then ! assume that vertzical positions are on file`	214	`elseif ( vmode.eq.'nil') then ! assume that vertzical positions are on file`
214	`continue`	215	`continue`
215	`else`	216	`else`
216	`print*,' ERROR: vertical mode not supported ',trim(vmode)`	217	`print*,' ERROR: vertical mode not supported ',trim(vmode)`
217	`stop`	218	`stop`
218	`endif`	219	`endif`
219		220
220	`c Read units of vertical axis`	221	`c Read units of vertical axis`
221	`read(10,*) umode`	222	`read(10,*) umode`
222	`if ( ( umode.ne.'hPa' ).and.`	223	`if ( ( umode.ne.'hPa' ).and.`
223	`> ( umode.ne.'hPa,agl' ).and.`	224	`> ( umode.ne.'hPa,agl' ).and.`
224	`> ( umode.ne.'K' ).and.`	225	`> ( umode.ne.'K' ).and.`
225	`> ( umode.ne.'PVU' ) )`	226	`> ( umode.ne.'PVU' ) )`
226	`> then`	227	`> then`
227	`print*,' ERROR: unit not supported ',trim(umode)`	228	`print*,' ERROR: unit not supported ',trim(umode)`
228	`stop`	229	`stop`
229	`endif`	230	`endif`
230		231
231	`c Read select string (dummy read)`	232	`c Read select string (dummy read)`
232	`read(10,*) selectstr`	233	`read(10,*) selectstr`
233		234
234	`c Read flag for 'no time check'`	235	`c Read flag for 'no time check'`
235	`read(10,*) timecheck`	236	`read(10,*) timecheck`
236		237
237	`c Close parameter file`	238	`c Close parameter file`
238	`close(10)`	239	`close(10)`
239		240
240	`c Decide which output format is used (1..4: trajectory format, -1: triple list)`	241	`c Decide which output format is used (1..4: trajectory format, -1: triple list)`
241	`call mode_tra(oformat,ofile)`	242	`call mode_tra(oformat,ofile)`
242		243
243	`c Decide whether all lat/lon/lev coordaintes are read from one file`	244	`c Decide whether all lat/lon/lev coordaintes are read from one file`
244	`if ( (hmode.eq.'file').and.(vmode.eq.'nil') ) then`	245	`if ( (hmode.eq.'file').and.(vmode.eq.'nil') ) then`
245	`hmode='file3'`	246	`hmode='file3'`
246	`elseif ( (hmode.eq.'file').and.(vmode.ne.'nil') ) then`	247	`elseif ( (hmode.eq.'file').and.(vmode.ne.'nil') ) then`
247	`hmode='file2'`	248	`hmode='file2'`
248	`endif`	249	`endif`
249		250
250	`c Convert timeshift (hh.mm) into a fractional time shift`	251	`c Convert timeshift (hh.mm) into a fractional time shift`
251	`call hhmm2frac(timeshift,tfrac)`	252	`call hhmm2frac(timeshift,tfrac)`
252	`if (tfrac.gt.0.) then`	253	`if (tfrac.gt.0.) then`
253	`tfrac=tfrac/timeinc`	254	`tfrac=tfrac/timeinc`
254	`else`	255	`else`
255	`tfrac=0.`	256	`tfrac=0.`
256	`endif`	257	`endif`
257		258
258	`c Read the region coordinates if needed`	259	`c Read the region coordinates if needed`
259	`if ( (hmode.eq.'region.eqd' ).or.`	260	`if ( (hmode.eq.'region.eqd' ).or.`
260	`> (hmode.eq.'region.grid') ) then`	261	`> (hmode.eq.'region.grid') ) then`
261		262
262	`open(10,file=regionf)`	263	`open(10,file=regionf)`
263		264
264	`50 read(10,*,end=51) string`	265	`50 read(10,*,end=51) string`
265		266
266	`if ( string(1:1).ne.'#' ) then`	267	`if ( string(1:1).ne.'#' ) then`
267	`call regionsplit(string,i,xcorner,ycorner)`	268	`call regionsplit(string,i,xcorner,ycorner)`
268	`if ( i.eq.iregion ) goto 52`	269	`if ( i.eq.iregion ) goto 52`
269	`endif`	270	`endif`
270		271
271	`goto 50`	272	`goto 50`
272		273
273	`51 close(10)`	274	`51 close(10)`
274		275
275	`print*,' ERROR: region ',iregion,' not found on ',trim(regionf)`	276	`print*,' ERROR: region ',iregion,' not found on ',trim(regionf)`
276	`stop`	277	`stop`
277		278
278	`52 continue`	279	`52 continue`
279		280
280	`endif`	281	`endif`
281		282
282	`c Write some status information`	283	`c Write some status information`
283	`print*,'---- INPUT PARAMETERS -----------------------------------'`	284	`print*,'---- INPUT PARAMETERS -----------------------------------'`
284	`print*`	285	`print*`
285	`if ( timeshift.gt.0. ) then`	286	`if ( timeshift.gt.0. ) then`
286	`print*,' P file : ',trim(pfile0),`	287	`print*,' P file : ',trim(pfile0),`
287	`> ' ',`	288	`> ' ',`
288	`> trim(pfile1)`	289	`> trim(pfile1)`
289	`print*,' S file : ',trim(sfile0),`	290	`print*,' S file : ',trim(sfile0),`
290	`> ' ',`	291	`> ' ',`
291	`> trim(sfile1)`	292	`> trim(sfile1)`
292	`else`	293	`else`
293	`print*,' P file : ',trim(pfile0)`	294	`print*,' P file : ',trim(pfile0)`
294	`print*,' S file : ',trim(sfile0)`	295	`print*,' S file : ',trim(sfile0)`
295	`endif`	296	`endif`
296	`print*,' Output file : ',trim(ofile)`	297	`print*,' Output file : ',trim(ofile)`
297	`print*`	298	`print*`
298	`if (oformat.eq.-1) then`	299	`if (oformat.eq.-1) then`
299	`print*,' Output format : (lon,lat,lev)-list'`	300	`print*,' Output format : (lon,lat,lev)-list'`
300	`else`	301	`else`
301	`print*,' Output format : ',oformat`	302	`print*,' Output format : ',oformat`
302	`endif`	303	`endif`
303	`print*`	304	`print*`
304	`print*,' Reference time (year) : ',reftime(1)`	305	`print*,' Reference time (year) : ',reftime(1)`
305	`print*,' (month) : ',reftime(2)`	306	`print*,' (month) : ',reftime(2)`
306	`print*,' (day) : ',reftime(3)`	307	`print*,' (day) : ',reftime(3)`
307	`print*,' (hour) : ',reftime(4)`	308	`print*,' (hour) : ',reftime(4)`
308	`print*,' (min) : ',reftime(5)`	309	`print*,' (min) : ',reftime(5)`
309	`print*,' Time range : ',reftime(6)`	310	`print*,' Time range : ',reftime(6)`
310	`print*`	311	`print*`
311	`print*,' Time shift : ',timeshift,' + ',`	312	`print*,' Time shift : ',timeshift,' + ',`
312	`> trim(pfile0)`	313	`> trim(pfile0)`
313	`print*,' Region file : ',trim(regionf)`	314	`print*,' Region file : ',trim(regionf)`
314	`print*`	315	`print*`
315	`print*,' hmode : ',trim(hmode)`	316	`print*,' hmode : ',trim(hmode)`
316	`if ( hmode.eq.'file2' ) then`	317	`if ( hmode.eq.'file2' ) then`
317	`print*,' filename [lat/lon] : ',trim(hfile)`	318	`print*,' filename [lat/lon] : ',trim(hfile)`
318	`elseif ( hmode.eq.'file3' ) then`	319	`elseif ( hmode.eq.'file3' ) then`
319	`print*,' filename [lat/lon/lev] : ',trim(hfile)`	320	`print*,' filename [lat/lon/lev] : ',trim(hfile)`
320	`elseif ( hmode.eq.'line' ) then`	321	`elseif ( hmode.eq.'line' ) then`
321	`write(*,'(a30,4f10.2,i4)')`	322	`write(*,'(a30,4f10.2,i4)')`
322	`> ' lon1,lon2,lat1,lat2,n : ',lon1,lon2,lat1,lat2,hn`	323	`> ' lon1,lon2,lat1,lat2,n : ',lon1,lon2,lat1,lat2,hn`
323	`elseif ( hmode.eq.'box.eqd' ) then`	324	`elseif ( hmode.eq.'box.eqd' ) then`
324	`write(*,'(a30,5f10.2)')`	325	`write(*,'(a30,5f10.2)')`
325	`> ' lon1,lon2,lat1,lat2,ds : ',lon1,lon2,lat1,lat2,ds`	326	`> ' lon1,lon2,lat1,lat2,ds : ',lon1,lon2,lat1,lat2,ds`
326	`elseif ( hmode.eq.'box.grid' ) then`	327	`elseif ( hmode.eq.'box.grid' ) then`
327	`write(*,'(a30,4f10.2)')`	328	`write(*,'(a30,4f10.2)')`
328	`> ' lon1,lon2,lat1,lat2 : ',lon1,lon2,lat1,lat2`	329	`> ' lon1,lon2,lat1,lat2 : ',lon1,lon2,lat1,lat2`
329	`elseif ( hmode.eq.'point' ) then`	330	`elseif ( hmode.eq.'point' ) then`
330	`print*,' lon,lat : ',lon1,lat1`	331	`print*,' lon,lat : ',lon1,lat1`
331	`elseif ( hmode.eq.'shift' ) then`	332	`elseif ( hmode.eq.'shift' ) then`
332	`write(*,'(a30,4f10.2)')`	333	`write(*,'(a30,4f10.2)')`
333	`> ' lon,lat,dlon,dlat : ',lon1,lat1,dlon,dlat`	334	`> ' lon,lat,dlon,dlat : ',lon1,lat1,dlon,dlat`
334	`elseif ( hmode.eq.'polygon.eqd' ) then`	335	`elseif ( hmode.eq.'polygon.eqd' ) then`
335	`write(*,'(a30,a10,f10.2)')`	336	`write(*,'(a30,a10,f10.2)')`
336	`> ' hfile, ds : ',trim(hfile),ds`	337	`> ' hfile, ds : ',trim(hfile),ds`
337	`elseif ( hmode.eq.'polygon.grid' ) then`	338	`elseif ( hmode.eq.'polygon.grid' ) then`
338	`write(*,'(a30,a10)')`	339	`write(*,'(a30,a10)')`
339	`> ' hfile : ',trim(hfile)`	340	`> ' hfile : ',trim(hfile)`
340	`elseif ( hmode.eq.'circle.eqd' ) then`	341	`elseif ( hmode.eq.'circle.eqd' ) then`
341	`write(*,'(a30,4f10.2)')`	342	`write(*,'(a30,4f10.2)')`
342	`> ' lonc,latc,radius, ds : ',lon1,lat1,radius,ds`	343	`> ' lonc,latc,radius, ds : ',lon1,lat1,radius,ds`
343	`elseif ( hmode.eq.'circle.grid' ) then`	344	`elseif ( hmode.eq.'circle.grid' ) then`
344	`write(*,'(a30,3f10.2)')`	345	`write(*,'(a30,3f10.2)')`
345	`> ' lonc,latc,radius : ',lon1,lat1,radius`	346	`> ' lonc,latc,radius : ',lon1,lat1,radius`
346	`elseif ( hmode.eq.'region.eqd' ) then`	347	`elseif ( hmode.eq.'region.eqd' ) then`
347	`write(*,'(a30,i4,1f10.2)')`	348	`write(*,'(a30,i4,1f10.2)')`
348	`> ' iregion, ds : ',iregion,ds`	349	`> ' iregion, ds : ',iregion,ds`
349	`write(*,'(a30,4f10.2)')`	350	`write(*,'(a30,4f10.2)')`
350	`> ' xcorner : ',(xcorner(i),i=1,4)`	351	`> ' xcorner : ',(xcorner(i),i=1,4)`
351	`write(*,'(a30,4f10.2)')`	352	`write(*,'(a30,4f10.2)')`
352	`> ' ycorner : ',(ycorner(i),i=1,4)`	353	`> ' ycorner : ',(ycorner(i),i=1,4)`
353	`elseif ( hmode.eq.'region.grid' ) then`	354	`elseif ( hmode.eq.'region.grid' ) then`
354	`write(*,'(a30,i4)')`	355	`write(*,'(a30,i4)')`
355	`> ' iregion : ',iregion`	356	`> ' iregion : ',iregion`
356	`write(*,'(a30,4f10.2)')`	357	`write(*,'(a30,4f10.2)')`
357	`> ' xcorner : ',(xcorner(i),i=1,4)`	358	`> ' xcorner : ',(xcorner(i),i=1,4)`
358	`write(*,'(a30,4f10.2)')`	359	`write(*,'(a30,4f10.2)')`
359	`> ' ycorner : ',(ycorner(i),i=1,4)`	360	`> ' ycorner : ',(ycorner(i),i=1,4)`
360	`endif`	361	`endif`
361	`print*`	362	`print*`
362	`print*,' vmode : ',trim(vmode)`	363	`print*,' vmode : ',trim(vmode)`
363	`if ( vmode.eq.'file') then`	364	`if ( vmode.eq.'file') then`
364	`print*,' filename : ',trim(vfile)`	365	`print*,' filename : ',trim(vfile)`
365	`elseif ( vmode.eq.'level' ) then`	366	`elseif ( vmode.eq.'level' ) then`
366	`print*,' level : ',lev1`	367	`print*,' level : ',lev1`
367	`elseif ( vmode.eq.'list') then`	368	`elseif ( vmode.eq.'list') then`
368	`print*,' n : ',vn`	369	`print*,' n : ',vn`
369	`print*,' level(i) : ',(levlist(i),i=1,vn)`	370	`print*,' level(i) : ',(levlist(i),i=1,vn)`
370	`elseif ( vmode.eq.'profile') then`	371	`elseif ( vmode.eq.'profile') then`
371	`print*,' lev1,lev2,n : ',lev1,lev2,vn`	372	`print*,' lev1,lev2,n : ',lev1,lev2,vn`
372	`endif`	373	`endif`
373	`print*`	374	`print*`
374	`print*,' umode : ',trim(umode)`	375	`print*,' umode : ',trim(umode)`
375	`print*`	376	`print*`
376	`print*,' time check : ',trim(timecheck)`	377	`print*,' time check : ',trim(timecheck)`
377	`print*`	378	`print*`
378		379
379	`c ------------------------------------------------------------------`	380	`c ------------------------------------------------------------------`
380	`c Read grid parameters from inital files`	381	`c Read grid parameters from inital files`
381	`c ------------------------------------------------------------------`	382	`c ------------------------------------------------------------------`
382		383
-		384	`c Init missing data value`
-		385	`mdv = -999.0`
-		386
383	`c Get the time of the first and second data file`	387	`c Get the time of the first and second data file`
384	`tstart = -timeshift ! Format hh.mm`	388	`tstart = -timeshift ! Format hh.mm`
385	`call hhmm2frac(tstart,frac)`	389	`call hhmm2frac(tstart,frac)`
386	`frac = frac + timeinc`	390	`frac = frac + timeinc`
387	`call frac2hhmm(frac,tend) ! Format hh.mm`	391	`call frac2hhmm(frac,tend) ! Format hh.mm`
388		392
389	`c Convert timeshift (hh.mm) into a fractional time shift`	393	`c Convert timeshift (hh.mm) into a fractional time shift`
390	`tfrac=real(int(timeshift))+`	394	`tfrac=real(int(timeshift))+`
391	`> 100.*(timeshift-real(int(timeshift)))/60.`	395	`> 100.*(timeshift-real(int(timeshift)))/60.`
392	`if (tfrac.gt.0.) then`	396	`if (tfrac.gt.0.) then`
393	`tfrac=tfrac/timeinc`	397	`tfrac=tfrac/timeinc`
394	`else`	398	`else`
395	`tfrac=0.`	399	`tfrac=0.`
396	`endif`	400	`endif`
397		401
398	`c Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,`	402	`c Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,`
399	`c pollon,pollat) The negative <-fid> of the file identifier is used`	403	`c pollon,pollat) The negative <-fid> of the file identifier is used`
400	`c as a flag for parameter retrieval`	404	`c as a flag for parameter retrieval`
401	`varname = 'U'`	405	`varname = 'U'`
402	`call input_open (fid,pfile0)`	406	`call input_open (fid,pfile0)`
403	`call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`	407	`call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`
404	`> tstart,pollon,pollat,rd,rd,nz,rd,rd,rd,timecheck)`	408	`> tstart,pollon,pollat,rd,rd,nz,rd,rd,rd,timecheck)`
405	`call input_close(fid)`	409	`call input_close(fid)`
406		410
407	`c Decide whether the grid has to be rotated`	411	`c Decide whether the grid has to be rotated`
408	`if ( abs(pollat-90.).gt.eps) then`	412	`if ( abs(pollat-90.).gt.eps) then`
409	`rotategrid = 1`	413	`rotategrid = 1`
410	`else`	414	`else`
411	`rotategrid = 0`	415	`rotategrid = 0`
412	`endif`	416	`endif`
413		417
414	`c Check whether region coordinates are within the domain`	418	`c Check whether region coordinates are within the domain`
415	`if ( (hmode.eq.'region.eqd' ).or.`	419	`if ( (hmode.eq.'region.eqd' ).or.`
416	`> (hmode.eq.'region.grid') ) then`	420	`> (hmode.eq.'region.grid') ) then`
417		421
418	`c Rotate coordinates`	422	`c Rotate coordinates`
419	`if ( rotategrid.eq.1) then`	423	`if ( rotategrid.eq.1) then`
420	`do i=1,4`	424	`do i=1,4`
421	`lon = xcorner(i)`	425	`lon = xcorner(i)`
422	`lat = ycorner(i)`	426	`lat = ycorner(i)`
423	`if ( rotategrid.eq.1) then`	427	`if ( rotategrid.eq.1) then`
424	`rlon = lmtolms(lat,lon,pollat,pollon)`	428	`rlon = lmtolms(lat,lon,pollat,pollon)`
425	`rlat = phtophs(lat,lon,pollat,pollon)`	429	`rlat = phtophs(lat,lon,pollat,pollon)`
426	`else`	430	`else`
427	`rlon = lon`	431	`rlon = lon`
428	`rlat = lat`	432	`rlat = lat`
429	`endif`	433	`endif`
430	`xcorner(i) = rlon`	434	`xcorner(i) = rlon`
431	`ycorner(i) = rlat`	435	`ycorner(i) = rlat`
432	`enddo`	436	`enddo`
433	`endif`	437	`endif`
434		438
435	`c Check`	439	`c Check`
436	`do i=1,4`	440	`do i=1,4`
437	`if ( (xcorner(i).lt.xmin).or.`	441	`if ( (xcorner(i).lt.xmin).or.`
438	`> (ycorner(i).lt.ymin).or.`	442	`> (ycorner(i).lt.ymin).or.`
439	`> (xcorner(i).gt.xmax).or.`	443	`> (xcorner(i).gt.xmax).or.`
440	`> (ycorner(i).gt.ymax) )`	444	`> (ycorner(i).gt.ymax) )`
441	`> then`	445	`> then`
442	`print*,' ERROR: region not included in data domain...'`	446	`print*,' ERROR: region not included in data domain...'`
443	`print*,' ',trim(string)`	447	`print*,' ',trim(string)`
444	`print*,' x: ',(xcorner(j),j=1,4)`	448	`print*,' x: ',(xcorner(j),j=1,4)`
445	`print*,' y: ',(ycorner(j),j=1,4)`	449	`print*,' y: ',(ycorner(j),j=1,4)`
446	`stop`	450	`stop`
447	`endif`	451	`endif`
448		452
449	`enddo`	453	`enddo`
450		454
451	`c Write rotated coordinates to logfile`	455	`c Write rotated coordinates to logfile`
452	`if ( rotategrid.eq.1) then`	456	`if ( rotategrid.eq.1) then`
453	`write(*,'(a30,i4)')`	457	`write(*,'(a30,i4)')`
454	`> ' Rotated region : ',iregion`	458	`> ' Rotated region : ',iregion`
455	`write(*,'(a30,4f10.2)')`	459	`write(*,'(a30,4f10.2)')`
456	`> ' xcorner (rotated) : ',(xcorner(i),i=1,4)`	460	`> ' xcorner (rotated) : ',(xcorner(i),i=1,4)`
457	`write(*,'(a30,4f10.2)')`	461	`write(*,'(a30,4f10.2)')`
458	`> ' ycorner (rotated) : ',(ycorner(i),i=1,4)`	462	`> ' ycorner (rotated) : ',(ycorner(i),i=1,4)`
459	`print*`	463	`print*`
460	`endif`	464	`endif`
461		465
462	`c Rotate coordinates back`	466	`c Rotate coordinates back`
463	`if ( rotategrid.eq.1 ) then`	467	`if ( rotategrid.eq.1 ) then`
464	`do i=1,4`	468	`do i=1,4`
465	`rlon = xcorner(i)`	469	`rlon = xcorner(i)`
466	`rlat = ycorner(i)`	470	`rlat = ycorner(i)`
467	`if ( rotategrid.eq.1) then`	471	`if ( rotategrid.eq.1) then`
468	`lon = lmstolm(rlat,rlon,pollat,pollon)`	472	`lon = lmstolm(rlat,rlon,pollat,pollon)`
469	`lat = phstoph(rlat,rlon,pollat,pollon)`	473	`lat = phstoph(rlat,rlon,pollat,pollon)`
470	`else`	474	`else`
471	`lon = rlon`	475	`lon = rlon`
472	`lat = rlat`	476	`lat = rlat`
473	`endif`	477	`endif`
474	`xcorner(i) = lon`	478	`xcorner(i) = lon`
475	`ycorner(i) = lat`	479	`ycorner(i) = lat`
476	`enddo`	480	`enddo`
477	`endif`	481	`endif`
478		482
479	`endif`	483	`endif`
480		484
481	`C Check if the number of levels is too large`	485	`C Check if the number of levels is too large`
482	`if (nz.gt.nlevmax) goto 993`	486	`if (nz.gt.nlevmax) goto 993`
483		487
484	`c Allocate memory for 3d arrays: pressure, theta, pv`	488	`c Allocate memory for 3d arrays: pressure, theta, pv`
485	`allocate(pr(nx,ny,nz),stat=stat)`	489	`allocate(pr(nx,ny,nz),stat=stat)`
486	`if (stat.ne.0) print,' error allocating array pr *'`	490	`if (stat.ne.0) print,' error allocating array pr *'`
487	`allocate(prs(nx,ny),stat=stat)`	491	`allocate(prs(nx,ny),stat=stat)`
488	`if (stat.ne.0) print,' error allocating array prs '`	492	`if (stat.ne.0) print,' error allocating array prs '`
489	`allocate(th(nx,ny,nz),stat=stat)`	493	`allocate(th(nx,ny,nz),stat=stat)`
490	`if (stat.ne.0) print,' error allocating array th *'`	494	`if (stat.ne.0) print,' error allocating array th *'`
491	`allocate(ths(nx,ny),stat=stat)`	495	`allocate(ths(nx,ny),stat=stat)`
492	`if (stat.ne.0) print,' error allocating array ths '`	496	`if (stat.ne.0) print,' error allocating array ths '`
493	`allocate(pv(nx,ny,nz),stat=stat)`	497	`allocate(pv(nx,ny,nz),stat=stat)`
494	`if (stat.ne.0) print,' error allocating array pv *'`	498	`if (stat.ne.0) print,' error allocating array pv *'`
495	`allocate(pvs(nx,ny),stat=stat)`	499	`allocate(pvs(nx,ny),stat=stat)`
496	`if (stat.ne.0) print,' error allocating array pvs '`	500	`if (stat.ne.0) print,' error allocating array pvs '`

Subversion Repositories lagranto.um

(root)/trunk/startf/startf.f/ – Rev 7 → 16