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=500)
 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)
+      real,allocatable, dimension (:,:,:) :: trainp                ! Input trajectories (ntra,ntim,ncol)
-      integer                                reftime(6)      ! Reference date
+      integer                                reftime(6)            ! Reference date
-      character*80                           varsinp(100)    ! Field names for input trajectory
+      character*80                           varsinp(100)          ! Field names for input trajectory
-      integer                                fid,fod         ! File identifier for inp and out trajectories
+      integer                                fid,fod               ! File identifier for inp and out trajectories
-      real                                   x0,y0,p0        ! Position of air parcel (physical space)
+      real                                   x0_tra,y0_tra,p0_tra  ! Position of air parcel (physical space)
-      real                                   reltpos0        ! Relative time of air parcel
+      real                                   reltpos0              ! Relative time of air parcel
-      real                                   xind,yind,pind  ! Position of air parcel (grid space)
+      real                                   xind,yind,pind        ! Position of air parcel (grid space)
-      integer                                fbflag          ! Flag for forward (1) or backward (-1) trajectories
+      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 (1000)    ! 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
       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  = 'U'
       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 horizontal position where to interpolate to
+c               Set the trajectory position
-                x0       = trainp(k,j,2)                          ! Longitude
+                x0_tra = trainp(k,j,2)                          ! Longitude
-                y0       = trainp(k,j,3)                          ! Latitude
+                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               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               Loop over pressure profile
+c               Loop over pressure profile (or other positions for horizontal mode)
-	            do l=1,npre
+	        do l=1,npre
+c                     Vertical
-c                 Set the pressure where to interpolate to
+                      if ( direction.eq.'vertical' ) then
+                        x0 = x0_tra
+                        y0 = y0_tra
-	              p0 = pmin + real(l-1)/real(npre-1) * (pmax-pmin)
+                        p0 = pmin + real(l-1)/real(npre-1) * (pmax-pmin)
-	              if ( centering.eq.'yes' )then
+                        if ( centering.eq.'yes' )then
-	              	  p0 = p0 + trainp(k,j,4)
+                           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
+                  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 the output trajectory position
+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) *
+                         ind_pre  = nint( real(npre) *
-     >          	       ( (trainp(k,j,4) - pmin)/(pmax-pmin) ) + 1.)
+     >          	       ( (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)
+     >             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)
+     >             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)
+     >             times,npre,ntim,0,1,direction)
 c         Exit point for loop over all tracing variables
 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)
+     >           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 coordinate Variables
+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, "standard_name","level")
-      ierr = nf90_put_att(ncID, varLevID, "units"        ,"hPa")
+        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)
 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,
-     >                   (/ varTimeID, varLevID /),varID)
+     >                   (/ 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)
 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)  +
+ZCOSPOL*           SIN(ZPHIS)) -
+COS(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
+      ZARG2   = COS(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
+ZCOSPOL*           SIN(ZPHIS)) +
+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

Subversion Repositories lagranto.ecmwf

(root)/trunk/lidar/lidar.f – Rev 3 → 5