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


      PROGRAM trace
 
C     ********************************************************************
C     *                                                                  *
C     * Pseudo-lidar plots along trajectories                             *
C     *                                                                  *
C     * Heini Wernli       first version:       April 1993               *
C     * Michael Sprenger   major upgrade:       2008-2009                *
C     *                                                                  *
C     ********************************************************************
 
      implicit none
      
c     --------------------------------------------------------------------
c     Declaration of parameters
c     --------------------------------------------------------------------
 
c     Maximum number of levels for input files
      integer   nlevmax
      parameter (nlevmax=100)
 
c     Maximum number of input files (dates, length of trajectories)
      integer   ndatmax
      parameter (ndatmax=5000)
 
c     Numerical epsilon (for float comparison)
      real      eps
      parameter (eps=0.001)
 
c     Conversion factors
      real      pi180                                   ! deg -> rad
      parameter (pi180=3.14159/180.)
      real      deg2km                                  ! deg -> km (at equator)
      parameter (deg2km=111.)
 
c     Prefix for primary and secondary fields
      character charp
      character chars
      parameter (charp='P')
      parameter (chars='S')
 
c     --------------------------------------------------------------------
c     Declaration of variables
c     --------------------------------------------------------------------
 
c     Input and output format for trajectories (see iotra.f)
      integer   inpmode
 
c     Input parameters
      character*80                           inpfile         ! Input trajectory file
      character*80                           outfile         ! Output netCDF file
      character*80                           outmode         ! Output mode (sum,mean)
      integer                                ntra            ! Number of trajectories
      integer                                ncol            ! Number of columns (including time, lon, lat, p)
      integer                                ntim            ! Number of times per trajectory
      integer                                ntrace0         ! Number of trace variables
      character*80                           tvar(200)       ! Tracing variable name (only the variable)
      character*1                            tfil(200)       ! Filename prefix 
      real                                   fac(200)        ! Scaling factor 
      integer                                compfl(200)     ! Computation flag (1=compute)
      integer                                numdat          ! Number of input files
      character*11                           dat(ndatmax)    ! Dates of input files
      real                                   timeinc         ! Time increment between input files
      real                                   tst             ! Time shift of start relative to first data file
      real                                   ten             ! Time shift of end relatiev to first data file  
      character*20                           startdate       ! First time/date on trajectory
      character*20                           enddate         ! Last time/date on trajectory
      character*80                           timecheck       ! Either 'yes' or 'no'
      character*80                           intmode         ! Interpolation mode ('normal', 'nearest')
	  real                                   pmin,pmax       ! Pressure range for output grid
	  integer                                npre            ! Number of pressure levels in output grid
	  character*80                           centering       ! Centering around trajectory position ('yes','no')
	  character*80                       direction       ! Direction of lidar (vertical,lat,lon,normal)
      character*80                           dumpcoord       ! Dumping coordinates ('yes','no')
 
c     Trajectories
      real,allocatable, dimension (:,:,:) :: trainp                ! Input trajectories (ntra,ntim,ncol)
      integer                                reftime(6)            ! Reference date
      character*80                           varsinp(100)          ! Field names for input trajectory
      integer                                fid,fod               ! File identifier for inp and out trajectories
      real                                   x0_tra,y0_tra,p0_tra  ! Position of air parcel (physical space)
      real                                   reltpos0              ! Relative time of air parcel
      real                                   xind,yind,pind        ! Position of air parcel (grid space)
      integer                                fbflag                ! Flag for forward (1) or backward (-1) trajectories
 
c     Meteorological fields from input file
      real,allocatable, dimension (:)     :: spt0,spt1       ! Surface pressure
      real,allocatable, dimension (:)     :: p3t0,p3t1       ! 3d-pressure 
      real,allocatable, dimension (:)     :: f3t0,f3t1       ! 3d field for tracing 
      character*80                           svars(100)      ! List of variables on S file
      character*80                           pvars(100)      ! List of variables on P file
      integer                                n_svars         ! Number of variables on S file
      integer                                n_pvars         ! Number of variables on P file
      
c     Input grid description
      real                                   pollon,pollat   ! Longitude/latitude of pole
      real                                   ak(100)         ! Vertical layers and levels
      real                                   bk(100) 
      real                                   xmin,xmax       ! Zonal grid extension
      real                                   ymin,ymax       ! Meridional grid extension
      integer                                nx,ny,nz        ! Grid dimensions
      real                                   dx,dy           ! Horizontal grid resolution
      integer                                hem             ! Flag for hemispheric domain
      integer                                per             ! Flag for periodic domain
      real                                   stagz           ! Vertical staggering
      real                                   mdv             ! Missing data value
 
c	  Output grid  and fields
      real                                   levels(1000)    ! Ouput levels
      real                                   times (5000)    ! Output times
	  real,allocatable, dimension (:,:)   :: out_pos         ! Position of trajectories
	  real,allocatable, dimension (:,:)   :: out_val         ! Output lidar field
	  real,allocatable, dimension (:,:)   :: out_cnt         ! # output lidar sum ups
 
 
c     Auxiliary variables
      integer                                i,j,k,l,n,m
      real                                   rd
      character*80                           filename
      real                                   time0,time1,reltpos
      integer                                itime0,itime1
      integer                                stat
      real                                   tstart
      integer                                iloaded0,iloaded1
      real                                   f0
      real                                   frac
      real                                   tload,tfrac
      integer                                isok
      character                              ch
      integer                                ind
      integer                                ind1,ind2,ind3,ind4,ind5
      integer                                ind6,ind7,ind8,ind9,ind0
      integer                                noutside
      real                                   delta
      integer                                itrace0
      character*80                           string
      character*80                           cdfname
      character*80                           varname
      real                                   time
      character*80                           longname
      character*80                           unit
      integer                                ind_time
      integer                                ind_pre
      real                                   rlat,rlon
      real                                   x0,y0,p0
      real                                   vx0,vy0,vx1,vy1
      real                                   rotation,lon,lat
 
c     Externals 
      real                                   int_index4
      external                               int_index4
 
c     --------------------------------------------------------------------
c     Start of program, Read parameters, get grid parameters
c     --------------------------------------------------------------------
 
c     Write start message
      print*,'========================================================='
      print*,'              *** START OF PROGRAM LIDAR ***'
      print*
 
c     Read parameters
      open(10,file='trace.param')
       read(10,*) inpfile
       read(10,*) outfile
       read(10,*) outmode
       read(10,*) startdate
       read(10,*) enddate 
       read(10,*) fbflag
       read(10,*) numdat
       if ( fbflag.eq.1) then
          do i=1,numdat
             read(10,'(a11)') dat(i)
          enddo
       else
          do i=numdat,1,-1
             read(10,'(a11)') dat(i)
          enddo
       endif
       read(10,*) timeinc
       read(10,*) tst
       read(10,*) ten
       read(10,*) ntra
       read(10,*) ntim
       read(10,*) ncol
       read(10,*) ntrace0
       do i=1,ntrace0
          read(10,*) tvar(i), fac(i), compfl(i), tfil(i)
       enddo
       read(10,*) n_pvars
       do i=1,n_pvars
          read(10,*) pvars(i)
       enddo
       read(10,*) n_svars
       do i=1,n_svars
          read(10,*) svars(i)
       enddo
       read(10,*) timecheck
       read(10,*) intmode
       read(10,*) pmin,pmax,npre
       read(10,*) centering
       read(10,*) direction
       read(10,*) dumpcoord
      close(10)
 
c     Check that the direction is ok
      if ( ( direction.ne.'vertical' ).and.
     >     ( direction.ne.'lat'      ).and.
     >     ( direction.ne.'lon'      ).and.
     >     ( direction.ne.'normal'   ) ) 
     >then 
         print*,' ERROR: invalid direction ',trim(direction)
         stop
      endif
 
c     Remove commented tracing fields
      itrace0 = 1
      do while ( itrace0.le.ntrace0) 
         string = tvar(itrace0)
         if ( string(1:1).eq.'#' ) then
            do i=itrace0,ntrace0-1
               tvar(i)   = tvar(i+1)
               fac(i)    = fac(i+1)
               compfl(i) = compfl(i+1)
               tfil(i)   = tfil(i+1)
            enddo
            ntrace0 = ntrace0 - 1
         else
            itrace0 = itrace0 + 1
         endif
      enddo
 
c     Set the formats of the input  files
      call mode_tra(inpmode,inpfile)
      if (inpmode.eq.-1) inpmode=1
 
C     Convert time shifts <tst,ten> from <hh.mm> into fractional time
      call hhmm2frac(tst,frac)
      tst = frac
      call hhmm2frac(ten,frac)
      ten = frac
 
c     Set the time for the first data file (depending on forward/backward mode)
      if (fbflag.eq.1) then
        tstart = -tst
      else
        tstart = tst
      endif
 
c     Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)
c     The negative <-fid> of the file identifier is used as a flag for parameter retrieval  
      filename = charp//dat(1)
      varname  = tvar(1)
      call input_open (fid,filename)
      call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >                 tstart,pollon,pollat,rd,rd,nz,rd,rd,rd,timecheck)
      call input_close(fid)
 
C     Allocate memory for some meteorological arrays
      allocate(spt0(nx*ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array spt0 ***'   ! Surface pressure
      allocate(spt1(nx*ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array spt1 ***'
      allocate(p3t0(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array p3t0 ***'   ! Pressure
      allocate(p3t1(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array p3t1 ***'
      allocate(f3t0(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array p3t0 ***'   ! Lidar field
      allocate(f3t1(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array p3t1 ***'
 
c	  Allocate memory for output field
	  allocate(out_pos(ntim,npre),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array out_pos ***'
	  allocate(out_val(ntim,npre),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array out_val ***'
	  allocate(out_cnt(ntim,npre),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array out_cnt ***'
 
C     Get memory for trajectory arrays
      allocate(trainp(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra      ***' 
 
c     Set the flags for periodic domains
      if ( abs(xmax-xmin-360.).lt.eps ) then
         per = 1
      elseif ( abs(xmax-xmin-360.+dx).lt.eps ) then
         per = 2
      else
         per = 0
      endif
 
C     Set logical flag for periodic data set (hemispheric or not)
      hem = 0
      if (per.eq.0.) then
         delta=xmax-xmin-360.
         if (abs(delta+dx).lt.eps) then               ! Program aborts: arrays must be closed
            print*,' ERROR: arrays must be closed... Stop'
         else if (abs(delta).lt.eps) then ! Periodic and hemispheric
           hem=1
           per=360.
        endif
      else                                            ! Periodic and hemispheric
         hem=1
      endif
 
c     Write some status information
      print*,'---- INPUT PARAMETERS -----------------------------------'
      print*
      print*,'  Input trajectory file  : ',trim(inpfile)
      print*,'  Format of input file   : ',inpmode
      print*,'  Output netCDF    file  : ',trim(outfile)
      print*,'  Format of output file  : ',trim(outmode)
      print*,'  Forward/backward       : ',fbflag
      print*,'  #tra                   : ',ntra
      print*,'  #col                   : ',ncol
      print*,'  #tim                   : ',ntim
      print*,'  No time check          : ',trim(timecheck)
      print*,'  Interpolation mode     : ',trim(intmode)
      do i=1,ntrace0
         if (compfl(i).eq.0) then
            print*,'  Tracing field          : ',
     >                 trim(tvar(i)), fac(i), ' 0 ', tfil(i)
         else
            print*,'  Tracing field          : ',
     >                trim(tvar(i)),' : online calc not supported'
         endif
      enddo
      print*,'  Output (pmin,pmax,n)   : ',pmin,pmax,npre
      print*,'  Centering              : ',trim(centering)
      print*,'  Orientation            : ',trim(direction)
      print*,'  Coordinate Dump        : ',trim(dumpcoord)
      print*
      print*,'---- INPUT DATA FILES -----------------------------------'
      print*
      call frac2hhmm(tstart,tload)
      print*,'  Time of 1st data file  : ',tload
      print*,'  #input files           : ',numdat
      print*,'  time increment         : ',timeinc
      call frac2hhmm(tst,tload)
      print*,'  Shift of start         : ',tload
      call frac2hhmm(ten,tload)
      print*,'  Shift of end           : ',tload
      print*,'  First/last input file  : ',trim(dat(1)),
     >                                     ' ... ',
     >                                     trim(dat(numdat)) 
      print*,'  Primary variables      : ',trim(pvars(1))
      do i=2,n_pvars
         print*,'                         : ',trim(pvars(i))
      enddo
      if ( n_svars.ge.1 ) then
         print*,'  Secondary variables    : ',trim(svars(1))
         do i=2,n_svars
            print*,'                         : ',trim(svars(i))
         enddo
      endif
      print*
      print*,'---- CONSTANT GRID PARAMETERS ---------------------------'
      print*
      print*,'  xmin,xmax     : ',xmin,xmax
      print*,'  ymin,ymax     : ',ymin,ymax
      print*,'  dx,dy         : ',dx,dy
      print*,'  pollon,pollat : ',pollon,pollat
      print*,'  nx,ny,nz      : ',nx,ny,nz
      print*,'  per, hem      : ',per,hem
      print*
 
c     --------------------------------------------------------------------
c     Load the input trajectories
c     --------------------------------------------------------------------
 
c     Read the input trajectory file
      call ropen_tra(fid,inpfile,ntra,ntim,ncol,reftime,varsinp,inpmode)
      call read_tra (fid,trainp,ntra,ntim,ncol,inpmode)
      call close_tra(fid,inpmode)
 
c     Check that first four columns correspond to time,lon,lat,p
      if ( (varsinp(1).ne.'time' ).or.
     >     (varsinp(2).ne.'xpos' ).and.(varsinp(2).ne.'lon' ).or.
     >     (varsinp(3).ne.'ypos' ).and.(varsinp(3).ne.'lat' ).or.
     >     (varsinp(4).ne.'ppos' ).and.(varsinp(4).ne.'p'   ) )
     >then
         print*,' ERROR: problem with input trajectories ...'
         stop
      endif
      varsinp(1) = 'time'
      varsinp(2) = 'lon'
      varsinp(3) = 'lat'
      varsinp(4) = 'p'
 
c     Write some status information of the input trajectories
      print*,'---- INPUT TRAJECTORIES ---------------------------------'
      print*
      print*,' Start date             : ',trim(startdate)
      print*,' End date               : ',trim(enddate)
      print*,' Reference time (year)  : ',reftime(1)
      print*,'                (month) : ',reftime(2)
      print*,'                (day)   : ',reftime(3)
      print*,'                (hour)  : ',reftime(4)
      print*,'                (min)   : ',reftime(5)
      print*,' Time range (min)       : ',reftime(6)
      do i=1,ncol
         print*,' Var                    :',i,trim(varsinp(i))
      enddo
      print*
 
c     Check that first time is 0 - otherwise the tracing will produce
c     wrong results because later in the code only absolute times are
c     considered: <itime0   = int(abs(tfrac-tstart)/timeinc) + 1>. This 
c     will be changed in a future version.
      if ( abs( trainp(1,1,1) ).gt.eps ) then
         print*,' ERROR: First time of trajectory must be 0, i.e. '
         print*,'     correspond to the reference date. Otherwise'
         print*,'     the tracing will give wrong results... STOP'
         stop
      endif
 
c     If requested, open the coordinate dump file
      if ( dumpcoord.eq.'yes' ) then
         open(10,file=trim(outfile)//'.coord')
      endif
 
c     --------------------------------------------------------------------
c     Trace the fields (fields available on input files)
c     --------------------------------------------------------------------
 
      print*
      print*,'---- LIDAR FROM PRIMARY AND SECONDARY DATA FILES ------'
 
c     Loop over all tracing fields
      do i=1,ntrace0
 
c	      Skip all fields marked for online calculation
          if ( compfl(i).eq.1 ) goto 110
         
c	      Init the output fields: position and lidar field
	      do k=1,ntim
	      	do l=1,npre
	      		out_pos(k,l) = 0.
	      	    out_val(k,l) = 0.
	      	    out_cnt(k,l) = 0.
	      	 enddo
	      enddo
 
c         Write some status information
          print*
          print*,' Now lidaring           : ',
     >         trim(tvar(i)),compfl(i),' ',trim(tfil(i))
 
c         Reset flags for load manager
          iloaded0 = -1
          iloaded1 = -1
 
c         Reset the counter for fields outside domain
          noutside = 0
 
c         Loop over all times
          do j=1,ntim
 
c            Convert trajectory time from hh.mm to fractional time
             call hhmm2frac(trainp(1,j,1),tfrac)
 
c            Get the times which are needed
             itime0   = int(abs(tfrac-tstart)/timeinc) + 1
             time0    = tstart + fbflag * real(itime0-1) * timeinc
             itime1   = itime0 + 1
             time1    = time0 + fbflag * timeinc
             if ( itime1.gt.numdat ) then
                itime1 = itime0
                time1  = time0
             endif
 
c            Load manager: Check whether itime0 can be copied from itime1
             if ( itime0.eq.iloaded1 ) then
                f3t0     = f3t1
                p3t0     = p3t1
                spt0     = spt1
                iloaded0 = itime0
             endif
 
c            Load manager: Check whether itime1 can be copied from itime0
             if ( itime1.eq.iloaded0 ) then
                f3t1     = f3t0
                p3t1     = p3t0
                spt1     = spt0
                iloaded1 = itime1
             endif
 
c            Load manager:  Load first time (tracing variable and grid)
             if ( itime0.ne.iloaded0 ) then
 
                filename = tfil(i)//dat(itime0)
                call frac2hhmm(time0,tload)
                varname  = tvar(i) 
                write(*,'(a23,a20,a3,a5,f7.2)') 
     >               '    ->  loading          : ',
     >               trim(filename),'  ',trim(varname),tload
                call input_open (fid,filename)                
                call input_wind 
     >               (fid,varname,f3t0,tload,stagz,mdv,
     >               xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck) 
    
                call input_grid      
     >               (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >                tload,pollon,pollat,p3t0,spt0,nz,ak,bk,stagz,
     >                timecheck)
                call input_close(fid)
                
                iloaded0 = itime0
 
             endif
 
c            Load manager: Load second time (tracing variable and grid)
             if ( itime1.ne.iloaded1 ) then
                
                filename = tfil(i)//dat(itime1)
                call frac2hhmm(time1,tload)
                varname  = tvar(i) 
                write(*,'(a23,a20,a3,a5,f7.2)') 
     >               '    ->  loading          : ',
     >               trim(filename),'  ',trim(varname),tload
                call input_open (fid,filename)
                call input_wind 
     >               (fid,varname,f3t1,tload,stagz,mdv,
     >               xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)      
                call input_grid      
     >               (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >               tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,
     >               timecheck)
                call input_close(fid)
                
                iloaded1 = itime1
                
             endif
 
c            Loop over all trajectories
             do k=1,ntra
                
c               Set the trajectory position 
                x0_tra = trainp(k,j,2)                          ! Longitude
                y0_tra = trainp(k,j,3)                          ! Latitude
                p0_tra = trainp(k,j,4)                          ! Pressure
 
c               Get rotation angle - orient normal to trajectory
                if ( direction.eq.'normal' ) then
 
                   vx0 = 1.                     
                   vy0 = 0.
 
                   if ( j.lt.ntim ) then
                      lat =  0.5 * ( trainp(k,j,3) + trainp(k,j+1,3) )
                      vx1 = ( trainp(k,j+1,2) - trainp(k,j,2) ) * 
     >                      cos( lat * pi180 )
                      vy1 = ( trainp(k,j+1,3) - trainp(k,j,3) )
                   else
                      lat =  0.5 * ( trainp(k,j,3) + trainp(k,j-1,3) )
                      vx1 = ( trainp(k,j,2) - trainp(k,j-1,2) ) * 
     >                      cos( lat * pi180 )
                      vy1 = ( trainp(k,j,3) - trainp(k,j-1,3) )
                   endif
 
                   if ( vx1.gt.180  ) vx1 = vx1 - 360
                   if ( vx1.lt.-180 ) vx1 = vx1 + 360.
 
                   call getangle (vx0,vy0,vx1,vy1,rotation)
                   rotation = -rotation
                   
                else
                   rotation = 0.
                endif
 
c               Set the relative time
                call hhmm2frac(trainp(k,j,1),tfrac)
                reltpos0 = fbflag * (tfrac-time0)/timeinc    
 
c               Loop over pressure profile (or other positions for horizontal mode)
	        do l=1,npre
 
c                     Vertical
                      if ( direction.eq.'vertical' ) then
                        x0 = x0_tra
                        y0 = y0_tra
                        p0 = pmin + real(l-1)/real(npre-1) * (pmax-pmin)
                        if ( centering.eq.'yes' )then
                           p0 = p0 + trainp(k,j,4)
                        endif
 
c                     Longitude
                      elseif ( direction.eq.'lon' ) then
                        x0 = pmin + real(l-1)/real(npre-1) * (pmax-pmin)
                        y0 = y0_tra
                        p0 = p0_tra
                        if ( centering.eq.'yes' )then
                           x0 = x0 + x0_tra
                        endif
                         
c                     Latitude
                      elseif ( direction.eq.'lat' ) then
                        x0 = x0_tra
                        y0 = pmin + real(l-1)/real(npre-1) * (pmax-pmin)
                        p0 = p0_tra
                        if ( centering.eq.'yes' )then
                           y0 = y0 + y0_tra
                        endif
 
c                     Normal to trajerctory
                      elseif ( direction.eq.'normal' ) then
 
c                        Set the coordinate in the rotated system
                         rlat = pmin + 
     >                             real(l-1)/real(npre-1) * (pmax-pmin)
                         rlon = 0.
 
c                        Transform it back to geographical lon/lat
                         call getenvir_b (x0_tra,y0_tra,rotation,
     >                                              x0,y0,rlon,rlat,1)
 
c                        Pressure unchanged
                         p0 = p0_tra
 
                      endif
 
c                     Handle periodic boundaries in zonal direction
                      if ( (x0.gt.xmax).and.(per.ne.0) ) x0 = x0 - 360.
                      if ( (x0.lt.xmin).and.(per.ne.0) ) x0 = x0 + 360.
 
c                     Handle pole problems for hemispheric data (taken from caltra.f)
                      if ((hem.eq.1).and.(y0.gt.90.)) then
                         y0=180.-y0
                         x0=x0+per/2.
                      endif
                      if ((hem.eq.1).and.(y0.lt.-90.)) then
                         y0=-180.-y0
                         x0=x0+per/2.
                      endif
                      if (y0.gt.89.99) then
                         y0=89.99
                      endif     
 
c                 If requested, dump the lidar coordinates
                  if ( (dumpcoord.eq.'yes').and.(i.eq.1) ) then
                     write(10,'(3f10.2)') x0,y0,trainp(k,j,1)
                     write(10,'(3f10.2)') x0_tra,y0_tra,5.
                  endif
 
C                 Get the index where to interpolate (x0,y0,p0)
                  if ( (abs(x0-mdv).gt.eps).and.
     >                 (abs(y0-mdv).gt.eps) )
     >            then
                     call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0,
     >                                p3t0,p3t1,spt0,spt1,3,
     >                                nx,ny,nz,xmin,ymin,dx,dy,mdv)
                  else
                     xind = mdv
                     yind = mdv
                     pind = mdv
                  endif
 
c                 If requested, apply nearest-neighbor interpolation
                  if ( intmode.eq.'nearest') then
                   
                     xind = real( nint(xind) )
                     yind = real( nint(yind) )
                     pind = real( nint(pind) )
                   
                     if ( xind.lt.1.  ) xind = 1.
                     if ( xind.gt.nx  ) xind = real(nx)
                     if ( yind.lt.1.  ) yind = 1.
                     if ( yind.gt.ny  ) yind = real(ny)
 
                     if ( pind.lt.1.  ) pind = 1.
                     if ( pind.gt.nz  ) pind = real(nz)
 
                  endif
 
c                 Do the interpolation: everthing is ok
                  if ( (xind.ge.1.).and.(xind.le.real(nx)).and.
     >                 (yind.ge.1.).and.(yind.le.real(ny)).and.
     >                 (pind.ge.1.).and.(pind.le.real(nz)) )
     >            then
                     f0 = int_index4(f3t0,f3t1,nx,ny,nz,
     >                               xind,yind,pind,reltpos0,mdv)
 
c                 Set to missing data
                  else
                     f0       = mdv
                  endif
 
c	              Save result to output array
                  if (abs(f0-mdv).gt.eps) then
                     out_val(j,l) = out_val(j,l) + f0 * fac(i)
                     out_cnt(j,l) = out_cnt(j,l) + 1.
 
	              endif
 
c              End loop over all pressure levels
	           enddo
 
c	           Save output - time index
	           ind_time = j
 
c                  Save output - space index for 'no centering'
	           if ( centering.eq.'no' ) then
                      if ( direction.eq.'vertical') then 
                         ind_pre  = nint( real(npre) *
     >          	       ( (p0_tra - pmin)/(pmax-pmin) ) + 1.)
                      elseif ( direction.eq.'lon') then 
                         ind_pre  = nint( real(npre) *
     >          	       ( (x0_tra - pmin)/(pmax-pmin) ) + 1.)
                      elseif ( direction.eq.'lat') then 
                         ind_pre  = nint( real(npre) *
     >          	       ( (y0_tra - pmin)/(pmax-pmin) ) + 1.)
                      endif
 
c                  Save output - space index for 'centering'
	           else
	           	  ind_pre  = nint( real(npre) *
     >          	       ( (0.            - pmin)/(pmax-pmin) ) + 1.)
	           endif
 
c                  Update the output array
	           if ( (ind_time.ge.1).and.(ind_time.le.ntim).and.
     >              (ind_pre .ge.1).and.(ind_pre .le.npre) )
     >         then
                    out_pos(ind_time,ind_pre) =
     >                          	out_pos(ind_time,ind_pre) + 1.
	           endif
 
c	         End loop over all trajectories
             enddo
 
c	      End loop over all times
          enddo
 
c	      Write the trajectory position to netCDF file - only once
	      if ( i.eq.1 ) then
	      	  cdfname  = outfile
	      	  varname  = 'POSITION'
	      	  longname = 'position of trajectory points'
	      	  unit     = 'none'
	      	  time     = 0.
              do k=1,npre
              	levels(k) = pmin + real(k-1)/real(npre-1) * (pmax-pmin)
              enddo
              do k=1,ntim
                 times(k) = trainp(1,k,1)
              enddo
              call writecdf2D_cf
     >            (cdfname,varname,longname,unit,out_pos,time,levels,
     >             times,npre,ntim,1,1,direction)
	      endif
 
c	      If no valid lidar count: set the field to missing data
          do k=1,ntim
          	do l=1,npre
          		if (abs(out_cnt(k,l)).lt.eps) then
          			out_val(k,l) = mdv
          	    endif
          	 enddo
          enddo
 
c	      If requested, calculate the mean of the lidar field
	      if ( outmode.eq.'mean' ) then
	      	do k=1,ntim
	      		do l=1,npre
	      			if ( (abs(out_val(k,l)-mdv).gt.eps).and.
     >                   (abs(out_cnt(k,l)    ).gt.0. ) )
     >              then
	      				out_val(k,l) = out_val(k,l) / out_cnt(k,l)
	      		    endif
	      		 enddo
	          enddo
	      endif
 
c	      Write the lidar field and count
	      cdfname  = outfile
	      if (outmode.eq.'sum' ) then
	         varname  = trim(tvar(i))//'_SUM'
	      elseif (outmode.eq.'mean' ) then
	         varname  = trim(tvar(i))//'_MEAN'
	      endif
	      longname = 'sum over all '//trim(tvar(i))//' profiles'
	      unit     = 'not given'
	      time     = 0.
          call writecdf2D_cf
     >            (cdfname,varname,longname,unit,out_val,time,levels,
     >             times,npre,ntim,0,1,direction)
 
	  	  cdfname  = outfile
	      varname  = trim(tvar(i))//'_CNT'
	      longname = 'counts of all '//trim(tvar(i))//' profiles'
	      unit     = 'not given'
	      time     = 0.
          call writecdf2D_cf
     >            (cdfname,varname,longname,unit,out_cnt,time,levels,
     >             times,npre,ntim,0,1,direction)
 
c         Exit point for loop over all tracing variables
 110      continue
 
c	   End loop over all lidar variables
       enddo
 
 
c     --------------------------------------------------------------------
c     Write output to netCDF file
c     --------------------------------------------------------------------
 
c     Write status information
      print*
      print*,'---- WRITE OUTPUT LIDAR FIELDS --------------------------'
      print*
 
c     Close coord dump file
      print*,' LIDAR written to      : ',trim(outfile)
      if ( dumpcoord.eq.'yes' ) then
        print*,' Coordinates dumped to : ',trim(outfile)//'.coord'
      endif
 
c     Write some status information, and end of program message
      print*  
      print*,'---- STATUS INFORMATION --------------------------------'
      print*
      print*,' ok'
      print*
      print*,'              *** END OF PROGRAM LIDAR ***'
      print*,'========================================================='
 
 
      end 
 
 
c     ********************************************************************
c     * INPUT / OUTPUT SUBROUTINES                                       *
c     ********************************************************************
 
c     --------------------------------------------------------------------
c     Subroutines to write 2D CF netcdf output file
c     --------------------------------------------------------------------
 
      subroutine writecdf2D_cf
     >          (cdfname,varname,longname,unit,arr,time,levels,times,
     >           npre,ntim,crefile,crevar,direction)
 
c     Create and write to the CF netcdf file <cdfname>. The variable
c     with name <varname> and with time <time> is written. The data
c     are in the two-dimensional array <arr>.  The flags <crefile> and
c     <crevar> determine whether the file and/or the variable should
c     be created.
 
      USE netcdf
 
      IMPLICIT NONE
 
c     Declaration of input parameters
      character*80 cdfname
      character*80 varname,longname,unit
      integer      npre,ntim
      real         arr(ntim,npre)
      real         levels(npre)
      real         times (ntim)
      real         time
      integer      crefile,crevar
      character*80 direction
 
c     Numerical epsilon
      real         eps
      parameter    (eps=1.e-5)
 
c     Local variables
      integer      ierr
      integer      ncID
      integer      LevDimId,    varLevID
      integer      TimeDimID,   varTimeID
      real         timeindex
      integer      i,j
      integer      nvars,varids(100)
      integer      ndims,dimids(100)
      real         timelist(1000)
      integer      ntimes
      integer      ind
      integer      varID
 
c     Quick an dirty solution for fieldname conflict
      if ( varname.eq.'time' ) varname = 'TIME'
 
c     Initially set error to indicate no errors.
      ierr = 0
 
c     ---- Create the netCDF - skip if <crefile=0> ----------------------
      if ( crefile.ne.1 ) goto 100
 
c     Create the file
      ierr = nf90_create(trim(cdfname), NF90_CLOBBER, ncID)
 
c     Define dimensions
      ierr=nf90_def_dim(ncID,'level',npre, LevDimID )
      ierr=nf90_def_dim(ncID,'time' ,ntim, TimeDimID)
 
c     Define space coordinate 
      ierr = nf90_def_var(ncID,'level',NF90_FLOAT,
     >     (/ LevDimID /),varLevID)
      if ( direction.eq.'vertical' ) then
        ierr = nf90_put_att(ncID, varLevID, "standard_name","level")
        ierr = nf90_put_att(ncID, varLevID, "units"        ,"hPa")
      elseif ( direction.eq.'lat' ) then
        ierr = nf90_put_att(ncID, varLevID, "standard_name","latitude")
        ierr = nf90_put_att(ncID, varLevID, "units"        ,"deg")
      elseif ( direction.eq.'lon' ) then
        ierr = nf90_put_att(ncID, varLevID, "standard_name","longitude")
        ierr = nf90_put_att(ncID, varLevID, "units"        ,"deg")
      elseif ( direction.eq.'normal' ) then
        ierr = nf90_put_att(ncID, varLevID, "standard_name","normal")
        ierr = nf90_put_att(ncID, varLevID, "units"        ,"deg")
      endif
 
c     Define time coordinate
      ierr = nf90_def_var(ncID,'time',NF90_FLOAT,
     >     (/ TimeDimID /), varTimeID)
      ierr = nf90_put_att(ncID, varTimeID, "long_name",    "time")
      ierr = nf90_put_att(ncID, varTimeID, "units",       "hours")
 
c     Write global attributes
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'Conventions', 'CF-1.0')
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'title',
     >     'pseudo-lidar from trajectory file')
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'source',
     >     'Lagranto Trajectories')
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'institution',
     >     'ETH Zurich, IACETH')
 
c     Check whether the definition was successful
      ierr = nf90_enddef(ncID)
      if (ierr.gt.0) then
         print*, 'An error occurred while attempting to ',
     >        'finish definition mode.'
         stop
      endif
 
c     Write coordinate data
      ierr = nf90_put_var(ncID,varLevID  ,levels)
      ierr = nf90_put_var(ncID,varTimeID ,times )
 
c     Close netCDF file
      ierr = nf90_close(ncID)
 
 100  continue
 
c     ---- Define a new variable - skip if <crevar=0> -----------------------
 
      if ( crevar.ne.1 ) goto 110
 
      print*,'Now defining ',trim(varname)
 
c     Open the file for read/write access
      ierr = nf90_open  (trim(cdfname), NF90_WRITE  , ncID)
 
c     Get the IDs for dimensions
      ierr = nf90_inq_dimid(ncID,'level', LevDimID )
      ierr = nf90_inq_dimid(ncID,'time' , TimeDimID)
 
c     Enter define mode
      ierr = nf90_redef(ncID)
 
c     Write definition and add attributes
      ierr = nf90_def_var(ncID,varname,NF90_FLOAT,
     >                   (/ TimeDimID, LevDimID /),varID)
      ierr = nf90_put_att(ncID, varID, "long_name" , longname )
      ierr = nf90_put_att(ncID, varID, "units"     , unit     )
      ierr = nf90_put_att(ncID, varID, '_FillValue', -999.99  )
 
c     Check whether definition was successful
      ierr = nf90_enddef(ncID)
      if (ierr.gt.0) then
         print*, 'An error occurred while attempting to ',
     >           'finish definition mode.'
         stop
      endif
 
c     Close netCDF file
      ierr = nf90_close(ncID)
 
 110  continue
 
c     ---- Write data --------------------------------------------------
 
c     Open the file for read/write access
      ierr = nf90_open  (trim(cdfname), NF90_WRITE , ncID)
 
c     Get the varID
      ierr = nf90_inq_varid(ncID,varname, varID )
      if (ierr.ne.0) then
         print*,'Variable ',trim(varname),' is not defined on ',
     >          trim(cdfname)
         stop
      endif
 
c     Write data block
      ierr = nf90_put_var(ncID,varID,arr,
     >                    start = (/ 1, 1 /),
     >                    count = (/ ntim, npre/) )
 
c     Check whether writing was successful
      ierr = nf90_close(ncID)
      if (ierr.ne.0) then
         write(*,*) trim(nf90_strerror(ierr))
         write(*,*) 'An error occurred while attempting to ',
     >              'close the netcdf file.'
         write(*,*) 'in clscdf_CF'
      endif
 
      end
 
c     ********************************************************************************
c     * Coordinate rotation - lidar normal to trajectory                             *
c     ********************************************************************************
 
c     --------------------------------------------------------------------------------
c     Backward coordinate transformation (Rotated lon/lat -> True lon/lat)
c     --------------------------------------------------------------------------------
 
      SUBROUTINE getenvir_b (clon,clat,rotation,
     >                       lon,lat,rlon,rlat,n)
 
      implicit none
 
c     Declaration of input and output parameters
      integer     n
      real        clon,clat,rotation
      real        lon(n), lat(n)
      real        rlon(n),rlat(n)
 
c     Auxiliary variables 
      real         pollon,pollat
      integer      i
      real         rlon1(n),rlat1(n)
 
c     Externals
      real         lmstolm,phstoph
      external     lmstolm,phstoph
 
c     First coordinate transformation (make the local coordinate system parallel to equator)
      pollon=-180.
      pollat=90.+rotation
      do i=1,n
         rlon1(i)=90.+lmstolm(rlat(i),rlon(i)-90.,pollat,pollon)
         rlat1(i)=phstoph(rlat(i),rlon(i)-90.,pollat,pollon)            
      enddo
 
c     Second coordinate transformation (make the local coordinate system parallel to equator)
      pollon=clon-180.
      if (pollon.lt.-180.) pollon=pollon+360.
      pollat=90.-clat
      do i=1,n
         lon(i)=lmstolm(rlat1(i),rlon1(i),pollat,pollon)
         lat(i)=phstoph(rlat1(i),rlon1(i),pollat,pollon)            
      enddo
 
      END
 
c     ---------------------------------------------------------------------
c     Determine the angle between two vectors
c     ---------------------------------------------------------------------
 
      SUBROUTINE getangle (vx1,vy1,vx2,vy2,angle)
 
c     Given two vectors <vx1,vy1> and <vx2,vy2>, determine the angle (in deg) 
c     between the two vectors.
 
      implicit none
 
c     Declaration of subroutine parameters
      real vx1,vy1
      real vx2,vy2
      real angle
 
c     Auxiliary variables and parameters
      real len1,len2,len3
      real val1,val2,val3
      real pi
      parameter (pi=3.14159265359)
 
      len1=sqrt(vx1*vx1+vy1*vy1)
      len2=sqrt(vx2*vx2+vy2*vy2)
 
      if ((len1.gt.0.).and.(len2.gt.0.)) then
         vx1=vx1/len1
         vy1=vy1/len1
         vx2=vx2/len2
         vy2=vy2/len2
         
         val1=vx1*vx2+vy1*vy2
         val2=-vy1*vx2+vx1*vy2
         
         len3=sqrt(val1*val1+val2*val2)
         
         if ( (val1.ge.0.).and.(val2.ge.0.) ) then
            val3=acos(val1/len3)
         else if ( (val1.lt.0.).and.(val2.ge.0.) ) then
            val3=pi-acos(abs(val1)/len3)
         else if ( (val1.ge.0.).and.(val2.le.0.) ) then
            val3=-acos(val1/len3)
         else if ( (val1.lt.0.).and.(val2.le.0.) ) then
            val3=-pi+acos(abs(val1)/len3)
         endif
      else
         val3=0.
      endif
      
      angle=180./pi*val3                                                                                     
 
      END
 
c     --------------------------------------------------------------------------------
c     Transformation routine: LMSTOLM and PHSTOPH from library gm2em            
c     --------------------------------------------------------------------------------
 
      REAL FUNCTION LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
C
C**** LMSTOLM  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   LAM = LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C**                IM ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
C**   AUSGABE-
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   D.MAJEWSKI
 
      REAL        LAMS,PHIS,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL = ZPIR18*POLLAM
      ZPHIS   = ZPIR18*PHIS
      ZLAMS   = LAMS
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
      ZLAMS   = ZPIR18*ZLAMS
 
      ZARG1   = SIN(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
     1                          ZCOSPOL*           SIN(ZPHIS)) -
     2          COS(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
      ZARG2   = COS(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
     1                          ZCOSPOL*           SIN(ZPHIS)) +
     2          SIN(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
      IF (ABS(ZARG2).LT.1.E-30) THEN
        IF (ABS(ZARG1).LT.1.E-30) THEN
          LMSTOLM =   0.0
        ELSEIF (ZARG1.GT.0.) THEN
              LMSTOLAM =  90.0
            ELSE
              LMSTOLAM = -90.0
            ENDIF
      ELSE
        LMSTOLM = ZRPI18*ATAN2(ZARG1,ZARG2)
      ENDIF
 
      RETURN
      END
 
 
      REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
C
C**** PHSTOPH  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C****                 ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   PHI = PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
C**   AUSGABE-
C**   PARAMETER:   WAHRE GEOGRAPHISCHE BREITE ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   D.MAJEWSKI
 
      REAL        LAMS,PHIS,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      SINPOL = SIN(ZPIR18*POLPHI)
      COSPOL = COS(ZPIR18*POLPHI)
      ZPHIS  = ZPIR18*PHIS
      ZLAMS  = LAMS
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
      ZLAMS  = ZPIR18*ZLAMS
      ARG     = COSPOL*COS(ZPHIS)*COS(ZLAMS) + SINPOL*SIN(ZPHIS)
 
      PHSTOPH = ZRPI18*ASIN(ARG)
 
      RETURN
      END
 
 
      REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
C
C**** LMTOLMS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
C**   AUSGABE-
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   G. DE MORSIER
 
      REAL        LAM,PHI,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL =     ZPIR18*POLLAM
      ZPHI    =     ZPIR18*PHI
      ZLAM    = LAM
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
      ZLAM    = ZPIR18*ZLAM
 
      ZARG1   = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI)
      ZARG2   = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI)
      IF (ABS(ZARG2).LT.1.E-30) THEN
        IF (ABS(ZARG1).LT.1.E-30) THEN
          LMTOLMS =   0.0
        ELSEIF (ZARG1.GT.0.) THEN
              LMTOLMS =  90.0
            ELSE
              LMTOLMS = -90.0
            ENDIF
      ELSE
        LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
      ENDIF
 
      RETURN
      END
 
 
      REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
C
C**** PHTOPHS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
C**   AUSGABE-
C**   PARAMETER:   ROTIERTE BREITE PHIS ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   G. DE MORSIER
 
      REAL        LAM,PHI,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL = ZPIR18*POLLAM
      ZPHI    = ZPIR18*PHI
      ZLAM    = LAM
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
      ZLAM    = ZPIR18*ZLAM
      ZARG    = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI)
 
      PHTOPHS = ZRPI18*ASIN(ZARG)
 
      RETURN
      END
      

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

Subversion Repositories lagranto.ecmwf

(root)/trunk/lidar/lidar.f – Rev 13 → 29