WebSVN – lagranto.icon – 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                                vector(200)     ! Flag for vectorial transformation
      integer                                nvector         ! Counter of vectorial variables
      integer                                numdat          ! Number of input files
      character*13                           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                                   zmin,zmax       ! Pressure range for output grid
	  integer                                nlev            ! Number of pressure levels in output grid
	  character*80                           centering       ! Centering around trajectory position ('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,y0,z0        ! Position of air parcel (physical space)
      real                                   reltpos0        ! Relative time of air parcel
      real                                   xind,yind,zind  ! 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 (:)     :: zbt0,zbt1       ! Topography
      real,allocatable, dimension (:)     :: z3t0,z3t1       ! 3d height
      real,allocatable, dimension (:)     :: f3t0,f3t1       ! 3d field for tracing 
      real,allocatable, dimension (:)     :: v3t0,v3t1       ! second component for vector field
      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                                   polgam          ! Grid rotation
      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,m
      real                                   rd
      character*80                           filename
      real                                   time0,time1,reltpos
      integer                                itime0,itime1
      integer                                stat
      real                                   tstart
      integer                                iloaded0,iloaded1
      real                                   f0,v0
      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                                   lon,lat,rlon,rlat
      character*80                           name
      integer                                ipoint
      real                                   urot,vrot,u,v
 
c     Externals 
      real                                   int_index4
      external                               int_index4
 
      real                                   lmtolms
      real                                   phtophs
      real                                   lmstolm
      real                                   phstoph
      external                               lmtolms,phtophs
      external                               lmstolm,phstoph
 
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,'(a)') dat(i)
          enddo
       else
          do i=numdat,1,-1
             read(10,'(a)') 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,*) zmin,zmax,nlev
       read(10,*) centering
      close(10)
    
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     Check whether it is a vectorial variable - if yes, split it into
c     its two components and mark it as vectorial
      nvector = 0
      i       = 1
      do while ( i.le.ntrace0 )
          name   = tvar(i)
          ipoint = 0
          do j=1,80
             if ( name(j:j).eq.'.' ) ipoint = j
     	  enddo
     	  if ( ipoint.ne.0 ) then
             nvector   = nvector + 1
     	  	 tvar(i)   = trim ( name(1:ipoint-1) )
     	  	 vector(i) = nvector
     	  	 do j=i+1,ntrace0
     	  	 	tvar  (j+1) = tvar  (j)
     	  	    fac   (j+1) = fac   (j)
     	  	    compfl(j+1) = compfl(j)
     	  	    tfil  (j+1) = tfil  (j)
     	  	 enddo
     	  	 tvar  (i+1) = trim( name(ipoint+1:80) )
     	  	 fac   (i+1) = fac   (i)
     	  	 compfl(i+1) = compfl(i)
     	  	 tfil  (i+1) = tfil  (i)
     	  	 vector(i+1) = nvector
     	  	 ntrace0     = ntrace0 + 1
     	  	 i = i + 1
     	  else
     	     vector(i) = 0
     	  endif
     	  i = i + 1
      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'
      nx       = 1
      ny       = 1
      nz       = 1
      call input_open (fid,filename)
      call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
     >                 tload,pollon,pollat,polgam,rd,rd,nz,rd,timecheck)
      call input_close(fid)
 
C     Allocate memory for some meteorological arrays
      allocate(zbt0(nx*ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array zbt0 ***'   ! Topography
      allocate(zbt1(nx*ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array zbt1 ***'
      allocate(z3t0(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array z3t0 ***'   ! 3d model height
      allocate(z3t1(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array z3t1 ***'
      allocate(f3t0(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array f3t0 ***'   ! Lidar field
      allocate(f3t1(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array f3t1 ***'
      allocate(v3t0(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array v3t0 ***'   ! Second component for vector field
      allocate(v3t1(nx*ny*nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array v3t1 ***'
 
c	  Allocate memory for output field
	  allocate(out_pos(ntim,nlev),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array out_pos ***'
	  allocate(out_val(ntim,nlev),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array out_val ***'
	  allocate(out_cnt(ntim,nlev),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
            write(*,'(a20,a10,f10.2,a3,1x,a1,a20,i2)')
     >                 '   Tracing field          : ',
     >                 trim(tvar(i)), fac(i), ' 0 ', tfil(i),
     >                 ' -> vector/scalar : ',vector(i)
         else
            print*,'  Tracing field          : ',
     >                trim(tvar(i)),' : online calc not supported'
         endif
      enddo
      print*,'  Output (zmin,zmax,n)   : ',zmin,zmax,nlev
      print*,'  Centering              : ',trim(centering)
      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*,'  polgam        : ',polgam
      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     Coordinate rotation
      do i=1,ntra
         do j=1,ntim
 
            lon = trainp(i,j,2)
            lat = trainp(i,j,3)
 
            if ( abs(pollat-90.).gt.eps) then
               rlon = lmtolms(lat,lon,pollat,pollon)
               rlat = phtophs(lat,lon,pollat,pollon)
            else
               rlon = lon
               rlat = lat
            endif
 
            trainp(i,j,2) = rlon
            trainp(i,j,3) = rlat
 
         enddo
      enddo
 
c     Check that first four columns correspond to time,lon,lat,z
      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.'zpos' ).and.(varsinp(4).ne.'z'   ) )
     >then
         print*,' ERROR: problem with input trajectories ...'
         stop
      endif
      varsinp(1) = 'time'
      varsinp(2) = 'lon'
      varsinp(3) = 'lat'
      varsinp(4) = 'z'
 
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     Run a simple consistency check
      call do_timecheck(dat,numdat,timeinc,fbflag,reftime)
 
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,nlev
	      		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         Special handling for vectorial fields: we need the second component
c	      of the vector field and must transform them. Remember the  two vector
c	      components in <ind0> and <ind1>.
	      if ( vector(i).ne.0 ) then
	         ind0 = i
	         do m=1,ntrace0
	           if ( (vector(i).eq.vector(m)).and.(i.ne.m) ) then
	              ind1 = m
	           endif
	         enddo
	      endif
 
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
                v3t0     = v3t1
                z3t0     = z3t1
                zbt0     = zbt1
                iloaded0 = itime0
             endif
 
c            Load manager: Check whether itime1 can be copied from itime0
             if ( itime1.eq.iloaded0 ) then
                f3t1     = f3t0
                v3t1     = v3t0
                z3t1     = z3t0
                zbt1     = zbt0
                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,polgam,z3t0,zbt0,nz,stagz,timecheck)
                call input_close(fid)
                
                if ( vector(i).ne.0 ) then
                   filename = tfil(i)//dat(itime0)
                   call frac2hhmm(time0,tload)
                   varname  = tvar(ind1)
                   write(*,'(a23,a20,a3,a5,f7.2)')
     >                     '    -> loading (vector) : ',
     >                      trim(filename),'  ',trim(varname),tload
                   call input_open (fid,filename)
                   call input_wind
     >                  (fid,varname,v3t0,tload,stagz,mdv,
     >                  xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
                   call input_close(fid)
                endif
 
                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,polgam,z3t1,zbt1,nz,stagz,timecheck)
                call input_close(fid)
                
                if ( vector(i).ne.0 ) then
                    filename = tfil(i)//dat(itime1)
                    call frac2hhmm(time1,tload)
                    varname  = tvar(ind1)
                    write(*,'(a23,a20,a3,a5,f7.2)')
     >                      '    -> loading (vector) : ',
     >                      trim(filename),'  ',trim(varname),tload
                    call input_open (fid,filename)
                    call input_wind
     >                   (fid,varname,v3t1,tload,stagz,mdv,
     >                    xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
                    call input_close(fid)
                endif
 
                iloaded1 = itime1
                
             endif
 
c            Loop over all trajectories
             do k=1,ntra
                
c               Set the horizontal position where to interpolate to
                x0       = trainp(k,j,2)                          ! Longitude
                y0       = trainp(k,j,3)                          ! Latitude
 
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 height profile
	            do l=1,nlev
 
c                 Set the pressure where to interpolate to
	              z0 = zmin + real(l-1)/real(nlev-1) * (zmax-zmin)
	              if ( centering.eq.'yes' )then
	              	  z0 = z0 + trainp(k,j,4)
	              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,zind,x0,y0,z0,reltpos0,
     >                                z3t0,z3t1,zbt0,zbt1,3,
     >                                nx,ny,nz,xmin,ymin,dx,dy,mdv)
                  else
                     xind = mdv
                     yind = mdv
                     zind = mdv
                  endif
 
c                 If requested, apply nearest-neighbor interpolation
                  if ( intmode.eq.'nearest') then
                   
                     xind = real( nint(xind) )
                     yind = real( nint(yind) )
                     zind = real( nint(zind) )
                   
                     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 ( zind.lt.1.  ) zind = 1.
                     if ( zind.gt.nz  ) zind = 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.
     >                 (zind.ge.1.).and.(zind.le.real(nz)) )
     >            then
                     f0 = int_index4(f3t0,f3t1,nx,ny,nz,
     >                               xind,yind,zind,reltpos0,mdv)
	              else
	              	 f0 = mdv
	              endif
 
c	              For a vector field, we need the second component
	              if((vector(i).ne.0).and.(abs(f0-mdv).gt.eps))then
	              	   v0 = int_index4(v3t0,v3t1,nx,ny,nz,
     >                                  xind,yind,zind,reltpos0,mdv)
	              else
	              	   v0 = mdv
	              endif
 
c	              Let's do a rotation for vector fields - we need the
c	              second vector component, but will not save it!
	              if ( (abs(f0-mdv).gt.eps).and.
     >                 (abs(v0-mdv).gt.eps) ) then
     >
                      if ( ind0.lt.ind1 ) then
                      	  urot = f0
                      	  vrot = v0
                      	  call uvrot2uv (urot,vrot,y0,x0,
     >                               	 pollat,pollon,u,v)
                      	  f0 = u
                      	  v0 = v
                      else
                      	  urot = v0
                      	  vrot = f0
                      	  call uvrot2uv (urot,vrot,y0,x0,
     >                               	 pollat,pollon,u,v)
                      	  f0 = v
                      	  v0 = u
                      endif
 
                  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
	           ind_time = j
	           if ( centering.eq.'no' ) then
	              ind_pre  = nint( real(nlev) *
     >          	       ( (trainp(k,j,4) - zmin)/(zmax-zmin) ) + 1.)
	           else
	           	  ind_pre  = nint( real(nlev) *
     >          	       ( (0.            - zmin)/(zmax-zmin) ) + 1.)
	           endif
 
	           if ( (ind_time.ge.1).and.(ind_time.le.ntim).and.
     >              (ind_pre .ge.1).and.(ind_pre .le.nlev) )
     >         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,nlev
              	levels(k) = zmin + real(k-1)/real(nlev-1) * (zmax-zmin)
              enddo
              do k=1,ntim
                 times(k) = trainp(1,k,1)
              enddo
              call writecdf2D_cf
     >            (cdfname,varname,longname,unit,out_pos,time,levels,
     >             times,nlev,ntim,1,1)
	      endif
 
c	      If no valid lidar count: set the field to missing data
          do k=1,ntim
          	do l=1,nlev
          		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,nlev
	      			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,nlev,ntim,0,1)
 
	  	  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,nlev,ntim,0,1)
 
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     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     ********************************************************************https://questionpro.com/t/CNFr6ZHxARs
 
c     --------------------------------------------------------------------
c     Subroutines to write the CF netcdf output file
c     --------------------------------------------------------------------
 
      subroutine writecdf2D_cf
     >          (cdfname,varname,longname,unit,arr,time,levels,times,
     >           npre,ntim,crefile,crevar)
 
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
 
c     Local variables
      integer      ierr
      integer      ncID
      integer      LevDimId,    varLevID
      integer      TimeDimID,   varTimeID
      real         timeindex
      integer      i
      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
      ierr = nf90_def_var(ncID,'level',NF90_FLOAT,
     >     (/ LevDimID /),varLevID)
      ierr = nf90_put_att(ncID, varLevID, "standard_name","level")
      ierr = nf90_put_att(ncID, varLevID, "units"        ,"m")
 
      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
 
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
 
 
         
         
      
 
 
      

Rev 3	Rev 5
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=500)`
25		25
26	`c Numerical epsilon (for float comparison)`	26	`c Numerical epsilon (for float comparison)`
27	`real eps`	27	`real eps`
28	`parameter (eps=0.001)`	28	`parameter (eps=0.001)`
29		29
30	`c 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 vector(200) ! Flag for vectorial transformation`	61	`integer vector(200) ! Flag for vectorial transformation`
62	`integer nvector ! Counter of vectorial variables`	62	`integer nvector ! Counter of vectorial variables`
63	`integer numdat ! Number of input files`	63	`integer numdat ! Number of input files`
64	`character*13 dat(ndatmax) ! Dates of input files`	64	`character*13 dat(ndatmax) ! Dates of input files`
65	`real timeinc ! Time increment between input files`	65	`real timeinc ! Time increment between input files`
66	`real tst ! Time shift of start relative to first data file`	66	`real tst ! Time shift of start relative to first data file`
67	`real ten ! Time shift of end relatiev to first data file`	67	`real ten ! Time shift of end relatiev to first data file`
68	`character*20 startdate ! First time/date on trajectory`	68	`character*20 startdate ! First time/date on trajectory`
69	`character*20 enddate ! Last time/date on trajectory`	69	`character*20 enddate ! Last time/date on trajectory`
70	`character*80 timecheck ! Either 'yes' or 'no'`	70	`character*80 timecheck ! Either 'yes' or 'no'`
71	`character*80 intmode ! Interpolation mode ('normal', 'nearest')`	71	`character*80 intmode ! Interpolation mode ('normal', 'nearest')`
72	`real zmin,zmax ! Pressure range for output grid`	72	`real zmin,zmax ! Pressure range for output grid`
73	`integer nlev ! Number of pressure levels in output grid`	73	`integer nlev ! Number of pressure levels in output grid`
74	`character*80 centering ! Centering around trajectory position ('yes','no')`	74	`character*80 centering ! Centering around trajectory position ('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,y0,z0 ! Position of air parcel (physical space)`	81	`real x0,y0,z0 ! 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,zind ! Position of air parcel (grid space)`	83	`real xind,yind,zind ! 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 (:) :: zbt0,zbt1 ! Topography`	87	`real,allocatable, dimension (:) :: zbt0,zbt1 ! Topography`
88	`real,allocatable, dimension (:) :: z3t0,z3t1 ! 3d height`	88	`real,allocatable, dimension (:) :: z3t0,z3t1 ! 3d height`
89	`real,allocatable, dimension (:) :: f3t0,f3t1 ! 3d field for tracing`	89	`real,allocatable, dimension (:) :: f3t0,f3t1 ! 3d field for tracing`
90	`real,allocatable, dimension (:) :: v3t0,v3t1 ! second component for vector field`	90	`real,allocatable, dimension (:) :: v3t0,v3t1 ! second component for vector field`
91	`character*80 svars(100) ! List of variables on S file`	91	`character*80 svars(100) ! List of variables on S file`
92	`character*80 pvars(100) ! List of variables on P file`	92	`character*80 pvars(100) ! List of variables on P file`
93	`integer n_svars ! Number of variables on S file`	93	`integer n_svars ! Number of variables on S file`
94	`integer n_pvars ! Number of variables on P file`	94	`integer n_pvars ! Number of variables on P file`
95		95
96	`c Input grid description`	96	`c Input grid description`
97	`real pollon,pollat ! Longitude/latitude of pole`	97	`real pollon,pollat ! Longitude/latitude of pole`
98	`real polgam ! Grid rotation`	98	`real polgam ! Grid rotation`
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 (1000) ! 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,v0`	125	`real f0,v0`
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 lon,lat,rlon,rlat`	144	`real lon,lat,rlon,rlat`
145	`character*80 name`	145	`character*80 name`
146	`integer ipoint`	146	`integer ipoint`
147	`real urot,vrot,u,v`	147	`real urot,vrot,u,v`
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	`real lmtolms`	153	`real lmtolms`
154	`real phtophs`	154	`real phtophs`
155	`real lmstolm`	155	`real lmstolm`
156	`real phstoph`	156	`real phstoph`
157	`external lmtolms,phtophs`	157	`external lmtolms,phtophs`
158	`external lmstolm,phstoph`	158	`external lmstolm,phstoph`
159		159
160	`c --------------------------------------------------------------------`	160	`c --------------------------------------------------------------------`
161	`c Start of program, Read parameters, get grid parameters`	161	`c Start of program, Read parameters, get grid parameters`
162	`c --------------------------------------------------------------------`	162	`c --------------------------------------------------------------------`
163		163
164	`c Write start message`	164	`c Write start message`
165	`print*,'========================================================='`	165	`print*,'========================================================='`
166	`print,' START OF PROGRAM LIDAR *'`	166	`print,' START OF PROGRAM LIDAR *'`
167	`print*`	167	`print*`
168		168
169	`c Read parameters`	169	`c Read parameters`
170	`open(10,file='trace.param')`	170	`open(10,file='trace.param')`
171	`read(10,*) inpfile`	171	`read(10,*) inpfile`
172	`read(10,*) outfile`	172	`read(10,*) outfile`
173	`read(10,*) outmode`	173	`read(10,*) outmode`
174	`read(10,*) startdate`	174	`read(10,*) startdate`
175	`read(10,*) enddate`	175	`read(10,*) enddate`
176	`read(10,*) fbflag`	176	`read(10,*) fbflag`
177	`read(10,*) numdat`	177	`read(10,*) numdat`
178	`if ( fbflag.eq.1) then`	178	`if ( fbflag.eq.1) then`
179	`do i=1,numdat`	179	`do i=1,numdat`
180	`read(10,'(a)') dat(i)`	180	`read(10,'(a)') dat(i)`
181	`enddo`	181	`enddo`
182	`else`	182	`else`
183	`do i=numdat,1,-1`	183	`do i=numdat,1,-1`
184	`read(10,'(a)') dat(i)`	184	`read(10,'(a)') dat(i)`
185	`enddo`	185	`enddo`
186	`endif`	186	`endif`
187	`read(10,*) timeinc`	187	`read(10,*) timeinc`
188	`read(10,*) tst`	188	`read(10,*) tst`
189	`read(10,*) ten`	189	`read(10,*) ten`
190	`read(10,*) ntra`	190	`read(10,*) ntra`
191	`read(10,*) ntim`	191	`read(10,*) ntim`
192	`read(10,*) ncol`	192	`read(10,*) ncol`
193	`read(10,*) ntrace0`	193	`read(10,*) ntrace0`
194	`do i=1,ntrace0`	194	`do i=1,ntrace0`
195	`read(10,*) tvar(i), fac(i), compfl(i), tfil(i)`	195	`read(10,*) tvar(i), fac(i), compfl(i), tfil(i)`
196	`enddo`	196	`enddo`
197	`read(10,*) n_pvars`	197	`read(10,*) n_pvars`
198	`do i=1,n_pvars`	198	`do i=1,n_pvars`
199	`read(10,*) pvars(i)`	199	`read(10,*) pvars(i)`
200	`enddo`	200	`enddo`
201	`read(10,*) n_svars`	201	`read(10,*) n_svars`
202	`do i=1,n_svars`	202	`do i=1,n_svars`
203	`read(10,*) svars(i)`	203	`read(10,*) svars(i)`
204	`enddo`	204	`enddo`
205	`read(10,*) timecheck`	205	`read(10,*) timecheck`
206	`read(10,*) intmode`	206	`read(10,*) intmode`
207	`read(10,*) zmin,zmax,nlev`	207	`read(10,*) zmin,zmax,nlev`
208	`read(10,*) centering`	208	`read(10,*) centering`
209	`close(10)`	209	`close(10)`
210		210
211	`c Remove commented tracing fields`	211	`c Remove commented tracing fields`
212	`itrace0 = 1`	212	`itrace0 = 1`
213	`do while ( itrace0.le.ntrace0)`	213	`do while ( itrace0.le.ntrace0)`
214	`string = tvar(itrace0)`	214	`string = tvar(itrace0)`
215	`if ( string(1:1).eq.'#' ) then`	215	`if ( string(1:1).eq.'#' ) then`
216	`do i=itrace0,ntrace0-1`	216	`do i=itrace0,ntrace0-1`
217	`tvar(i) = tvar(i+1)`	217	`tvar(i) = tvar(i+1)`
218	`fac(i) = fac(i+1)`	218	`fac(i) = fac(i+1)`
219	`compfl(i) = compfl(i+1)`	219	`compfl(i) = compfl(i+1)`
220	`tfil(i) = tfil(i+1)`	220	`tfil(i) = tfil(i+1)`
221	`enddo`	221	`enddo`
222	`ntrace0 = ntrace0 - 1`	222	`ntrace0 = ntrace0 - 1`
223	`else`	223	`else`
224	`itrace0 = itrace0 + 1`	224	`itrace0 = itrace0 + 1`
225	`endif`	225	`endif`
226	`enddo`	226	`enddo`
227		227
228	`c Check whether it is a vectorial variable - if yes, split it into`	228	`c Check whether it is a vectorial variable - if yes, split it into`
229	`c its two components and mark it as vectorial`	229	`c its two components and mark it as vectorial`
230	`nvector = 0`	230	`nvector = 0`
231	`i = 1`	231	`i = 1`
232	`do while ( i.le.ntrace0 )`	232	`do while ( i.le.ntrace0 )`
233	`name = tvar(i)`	233	`name = tvar(i)`
234	`ipoint = 0`	234	`ipoint = 0`
235	`do j=1,80`	235	`do j=1,80`
236	`if ( name(j:j).eq.'.' ) ipoint = j`	236	`if ( name(j:j).eq.'.' ) ipoint = j`
237	`enddo`	237	`enddo`
238	`if ( ipoint.ne.0 ) then`	238	`if ( ipoint.ne.0 ) then`
239	`nvector = nvector + 1`	239	`nvector = nvector + 1`
240	`tvar(i) = trim ( name(1:ipoint-1) )`	240	`tvar(i) = trim ( name(1:ipoint-1) )`
241	`vector(i) = nvector`	241	`vector(i) = nvector`
242	`do j=i+1,ntrace0`	242	`do j=i+1,ntrace0`
243	`tvar (j+1) = tvar (j)`	243	`tvar (j+1) = tvar (j)`
244	`fac (j+1) = fac (j)`	244	`fac (j+1) = fac (j)`
245	`compfl(j+1) = compfl(j)`	245	`compfl(j+1) = compfl(j)`
246	`tfil (j+1) = tfil (j)`	246	`tfil (j+1) = tfil (j)`
247	`enddo`	247	`enddo`
248	`tvar (i+1) = trim( name(ipoint+1:80) )`	248	`tvar (i+1) = trim( name(ipoint+1:80) )`
249	`fac (i+1) = fac (i)`	249	`fac (i+1) = fac (i)`
250	`compfl(i+1) = compfl(i)`	250	`compfl(i+1) = compfl(i)`
251	`tfil (i+1) = tfil (i)`	251	`tfil (i+1) = tfil (i)`
252	`vector(i+1) = nvector`	252	`vector(i+1) = nvector`
253	`ntrace0 = ntrace0 + 1`	253	`ntrace0 = ntrace0 + 1`
254	`i = i + 1`	254	`i = i + 1`
255	`else`	255	`else`
256	`vector(i) = 0`	256	`vector(i) = 0`
257	`endif`	257	`endif`
258	`i = i + 1`	258	`i = i + 1`
259	`enddo`	259	`enddo`
260		260
261	`c Set the formats of the input files`	261	`c Set the formats of the input files`
262	`call mode_tra(inpmode,inpfile)`	262	`call mode_tra(inpmode,inpfile)`
263	`if (inpmode.eq.-1) inpmode=1`	263	`if (inpmode.eq.-1) inpmode=1`
264		264
265	`C Convert time shifts <tst,ten> from <hh.mm> into fractional time`	265	`C Convert time shifts <tst,ten> from <hh.mm> into fractional time`
266	`call hhmm2frac(tst,frac)`	266	`call hhmm2frac(tst,frac)`
267	`tst = frac`	267	`tst = frac`
268	`call hhmm2frac(ten,frac)`	268	`call hhmm2frac(ten,frac)`
269	`ten = frac`	269	`ten = frac`
270		270
271	`c Set the time for the first data file (depending on forward/backward mode)`	271	`c Set the time for the first data file (depending on forward/backward mode)`
272	`if (fbflag.eq.1) then`	272	`if (fbflag.eq.1) then`
273	`tstart = -tst`	273	`tstart = -tst`
274	`else`	274	`else`
275	`tstart = tst`	275	`tstart = tst`
276	`endif`	276	`endif`
277		277
278	`c Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)`	278	`c Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)`
279	`c The negative <-fid> of the file identifier is used as a flag for parameter retrieval`	279	`c The negative <-fid> of the file identifier is used as a flag for parameter retrieval`
280	`filename = charp//dat(1)`	280	`filename = charp//dat(1)`
281	`varname = 'U'`	281	`varname = 'U'`
282	`nx = 1`	282	`nx = 1`
283	`ny = 1`	283	`ny = 1`
284	`nz = 1`	284	`nz = 1`
285	`call input_open (fid,filename)`	285	`call input_open (fid,filename)`
286	`call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`	286	`call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,`
287	`> tload,pollon,pollat,polgam,rd,rd,nz,rd,timecheck)`	287	`> tload,pollon,pollat,polgam,rd,rd,nz,rd,timecheck)`
288	`call input_close(fid)`	288	`call input_close(fid)`
289		289
290	`C Allocate memory for some meteorological arrays`	290	`C Allocate memory for some meteorological arrays`
291	`allocate(zbt0(nx*ny),stat=stat)`	291	`allocate(zbt0(nx*ny),stat=stat)`
292	`if (stat.ne.0) print,' error allocating array zbt0 *' ! Topography`	292	`if (stat.ne.0) print,' error allocating array zbt0 *' ! Topography`
293	`allocate(zbt1(nx*ny),stat=stat)`	293	`allocate(zbt1(nx*ny),stat=stat)`
294	`if (stat.ne.0) print,' error allocating array zbt1 *'`	294	`if (stat.ne.0) print,' error allocating array zbt1 *'`
295	`allocate(z3t0(nxnynz),stat=stat)`	295	`allocate(z3t0(nxnynz),stat=stat)`
296	`if (stat.ne.0) print,' error allocating array z3t0 *' ! 3d model height`	296	`if (stat.ne.0) print,' error allocating array z3t0 *' ! 3d model height`
297	`allocate(z3t1(nxnynz),stat=stat)`	297	`allocate(z3t1(nxnynz),stat=stat)`
298	`if (stat.ne.0) print,' error allocating array z3t1 *'`	298	`if (stat.ne.0) print,' error allocating array z3t1 *'`
299	`allocate(f3t0(nxnynz),stat=stat)`	299	`allocate(f3t0(nxnynz),stat=stat)`
300	`if (stat.ne.0) print,' error allocating array f3t0 *' ! Lidar field`	300	`if (stat.ne.0) print,' error allocating array f3t0 *' ! Lidar field`
301	`allocate(f3t1(nxnynz),stat=stat)`	301	`allocate(f3t1(nxnynz),stat=stat)`
302	`if (stat.ne.0) print,' error allocating array f3t1 *'`	302	`if (stat.ne.0) print,' error allocating array f3t1 *'`
303	`allocate(v3t0(nxnynz),stat=stat)`	303	`allocate(v3t0(nxnynz),stat=stat)`
304	`if (stat.ne.0) print,' error allocating array v3t0 *' ! Second component for vector field`	304	`if (stat.ne.0) print,' error allocating array v3t0 *' ! Second component for vector field`
305	`allocate(v3t1(nxnynz),stat=stat)`	305	`allocate(v3t1(nxnynz),stat=stat)`
306	`if (stat.ne.0) print,' error allocating array v3t1 *'`	306	`if (stat.ne.0) print,' error allocating array v3t1 *'`
307		307
308	`c Allocate memory for output field`	308	`c Allocate memory for output field`
309	`allocate(out_pos(ntim,nlev),stat=stat)`	309	`allocate(out_pos(ntim,nlev),stat=stat)`
310	`if (stat.ne.0) print,' error allocating array out_pos *'`	310	`if (stat.ne.0) print,' error allocating array out_pos *'`
311	`allocate(out_val(ntim,nlev),stat=stat)`	311	`allocate(out_val(ntim,nlev),stat=stat)`
312	`if (stat.ne.0) print,' error allocating array out_val *'`	312	`if (stat.ne.0) print,' error allocating array out_val *'`
313	`allocate(out_cnt(ntim,nlev),stat=stat)`	313	`allocate(out_cnt(ntim,nlev),stat=stat)`
314	`if (stat.ne.0) print,' error allocating array out_cnt *'`	314	`if (stat.ne.0) print,' error allocating array out_cnt *'`
315		315
316	`C Get memory for trajectory arrays`	316	`C Get memory for trajectory arrays`
317	`allocate(trainp(ntra,ntim,ncol),stat=stat)`	317	`allocate(trainp(ntra,ntim,ncol),stat=stat)`
318	`if (stat.ne.0) print,' error allocating array tra *'`	318	`if (stat.ne.0) print,' error allocating array tra *'`
319		319
320	`c Set the flags for periodic domains`	320	`c Set the flags for periodic domains`
321	`if ( abs(xmax-xmin-360.).lt.eps ) then`	321	`if ( abs(xmax-xmin-360.).lt.eps ) then`
322	`per = 1`	322	`per = 1`
323	`elseif ( abs(xmax-xmin-360.+dx).lt.eps ) then`	323	`elseif ( abs(xmax-xmin-360.+dx).lt.eps ) then`
324	`per = 2`	324	`per = 2`
325	`else`	325	`else`
326	`per = 0`	326	`per = 0`
327	`endif`	327	`endif`
328		328
329	`C Set logical flag for periodic data set (hemispheric or not)`	329	`C Set logical flag for periodic data set (hemispheric or not)`
330	`hem = 0`	330	`hem = 0`
331	`if (per.eq.0.) then`	331	`if (per.eq.0.) then`
332	`delta=xmax-xmin-360.`	332	`delta=xmax-xmin-360.`
333	`if (abs(delta+dx).lt.eps) then ! Program aborts: arrays must be closed`	333	`if (abs(delta+dx).lt.eps) then ! Program aborts: arrays must be closed`
334	`print*,' ERROR: arrays must be closed... Stop'`	334	`print*,' ERROR: arrays must be closed... Stop'`
335	`else if (abs(delta).lt.eps) then ! Periodic and hemispheric`	335	`else if (abs(delta).lt.eps) then ! Periodic and hemispheric`
336	`hem=1`	336	`hem=1`
337	`per=360.`	337	`per=360.`
338	`endif`	338	`endif`
339	`else ! Periodic and hemispheric`	339	`else ! Periodic and hemispheric`
340	`hem=1`	340	`hem=1`
341	`endif`	341	`endif`
342		342
343	`c Write some status information`	343	`c Write some status information`
344	`print*,'---- INPUT PARAMETERS -----------------------------------'`	344	`print*,'---- INPUT PARAMETERS -----------------------------------'`
345	`print*`	345	`print*`
346	`print*,' Input trajectory file : ',trim(inpfile)`	346	`print*,' Input trajectory file : ',trim(inpfile)`
347	`print*,' Format of input file : ',inpmode`	347	`print*,' Format of input file : ',inpmode`
348	`print*,' Output netCDF file : ',trim(outfile)`	348	`print*,' Output netCDF file : ',trim(outfile)`
349	`print*,' Format of output file : ',trim(outmode)`	349	`print*,' Format of output file : ',trim(outmode)`
350	`print*,' Forward/backward : ',fbflag`	350	`print*,' Forward/backward : ',fbflag`
351	`print*,' #tra : ',ntra`	351	`print*,' #tra : ',ntra`
352	`print*,' #col : ',ncol`	352	`print*,' #col : ',ncol`
353	`print*,' #tim : ',ntim`	353	`print*,' #tim : ',ntim`
354	`print*,' No time check : ',trim(timecheck)`	354	`print*,' No time check : ',trim(timecheck)`
355	`print*,' Interpolation mode : ',trim(intmode)`	355	`print*,' Interpolation mode : ',trim(intmode)`
356	`do i=1,ntrace0`	356	`do i=1,ntrace0`
357	`if (compfl(i).eq.0) then`	357	`if (compfl(i).eq.0) then`
358	`write(*,'(a20,a10,f10.2,a3,1x,a1,a20,i2)')`	358	`write(*,'(a20,a10,f10.2,a3,1x,a1,a20,i2)')`
359	`> ' Tracing field : ',`	359	`> ' Tracing field : ',`
360	`> trim(tvar(i)), fac(i), ' 0 ', tfil(i),`	360	`> trim(tvar(i)), fac(i), ' 0 ', tfil(i),`
361	`> ' -> vector/scalar : ',vector(i)`	361	`> ' -> vector/scalar : ',vector(i)`
362	`else`	362	`else`
363	`print*,' Tracing field : ',`	363	`print*,' Tracing field : ',`
364	`> trim(tvar(i)),' : online calc not supported'`	364	`> trim(tvar(i)),' : online calc not supported'`
365	`endif`	365	`endif`
366	`enddo`	366	`enddo`
367	`print*,' Output (zmin,zmax,n) : ',zmin,zmax,nlev`	367	`print*,' Output (zmin,zmax,n) : ',zmin,zmax,nlev`
368	`print*,' Centering : ',trim(centering)`	368	`print*,' Centering : ',trim(centering)`
369	`print*`	369	`print*`
370	`print*,'---- INPUT DATA FILES -----------------------------------'`	370	`print*,'---- INPUT DATA FILES -----------------------------------'`
371	`print*`	371	`print*`
372	`call frac2hhmm(tstart,tload)`	372	`call frac2hhmm(tstart,tload)`
373	`print*,' Time of 1st data file : ',tload`	373	`print*,' Time of 1st data file : ',tload`
374	`print*,' #input files : ',numdat`	374	`print*,' #input files : ',numdat`
375	`print*,' time increment : ',timeinc`	375	`print*,' time increment : ',timeinc`
376	`call frac2hhmm(tst,tload)`	376	`call frac2hhmm(tst,tload)`
377	`print*,' Shift of start : ',tload`	377	`print*,' Shift of start : ',tload`
378	`call frac2hhmm(ten,tload)`	378	`call frac2hhmm(ten,tload)`
379	`print*,' Shift of end : ',tload`	379	`print*,' Shift of end : ',tload`
380	`print*,' First/last input file : ',trim(dat(1)),`	380	`print*,' First/last input file : ',trim(dat(1)),`
381	`> ' ... ',`	381	`> ' ... ',`
382	`> trim(dat(numdat))`	382	`> trim(dat(numdat))`
383	`print*,' Primary variables : ',trim(pvars(1))`	383	`print*,' Primary variables : ',trim(pvars(1))`
384	`do i=2,n_pvars`	384	`do i=2,n_pvars`
385	`print*,' : ',trim(pvars(i))`	385	`print*,' : ',trim(pvars(i))`
386	`enddo`	386	`enddo`
387	`if ( n_svars.ge.1 ) then`	387	`if ( n_svars.ge.1 ) then`
388	`print*,' Secondary variables : ',trim(svars(1))`	388	`print*,' Secondary variables : ',trim(svars(1))`
389	`do i=2,n_svars`	389	`do i=2,n_svars`
390	`print*,' : ',trim(svars(i))`	390	`print*,' : ',trim(svars(i))`
391	`enddo`	391	`enddo`
392	`endif`	392	`endif`
393	`print*`	393	`print*`
394	`print*,'---- CONSTANT GRID PARAMETERS ---------------------------'`	394	`print*,'---- CONSTANT GRID PARAMETERS ---------------------------'`
395	`print*`	395	`print*`
396	`print*,' xmin,xmax : ',xmin,xmax`	396	`print*,' xmin,xmax : ',xmin,xmax`
397	`print*,' ymin,ymax : ',ymin,ymax`	397	`print*,' ymin,ymax : ',ymin,ymax`
398	`print*,' dx,dy : ',dx,dy`	398	`print*,' dx,dy : ',dx,dy`
399	`print*,' pollon,pollat : ',pollon,pollat`	399	`print*,' pollon,pollat : ',pollon,pollat`
400	`print*,' polgam : ',polgam`	400	`print*,' polgam : ',polgam`
401	`print*,' nx,ny,nz : ',nx,ny,nz`	401	`print*,' nx,ny,nz : ',nx,ny,nz`
402	`print*,' per, hem : ',per,hem`	402	`print*,' per, hem : ',per,hem`
403	`print*`	403	`print*`
404		404
405	`c --------------------------------------------------------------------`	405	`c --------------------------------------------------------------------`
406	`c Load the input trajectories`	406	`c Load the input trajectories`
407	`c --------------------------------------------------------------------`	407	`c --------------------------------------------------------------------`
408		408
409	`c Read the input trajectory file`	409	`c Read the input trajectory file`
410	`call ropen_tra(fid,inpfile,ntra,ntim,ncol,reftime,varsinp,inpmode)`	410	`call ropen_tra(fid,inpfile,ntra,ntim,ncol,reftime,varsinp,inpmode)`
411	`call read_tra (fid,trainp,ntra,ntim,ncol,inpmode)`	411	`call read_tra (fid,trainp,ntra,ntim,ncol,inpmode)`
412	`call close_tra(fid,inpmode)`	412	`call close_tra(fid,inpmode)`
413		413
414	`c Coordinate rotation`	414	`c Coordinate rotation`
415	`do i=1,ntra`	415	`do i=1,ntra`
416	`do j=1,ntim`	416	`do j=1,ntim`
417		417
418	`lon = trainp(i,j,2)`	418	`lon = trainp(i,j,2)`
419	`lat = trainp(i,j,3)`	419	`lat = trainp(i,j,3)`
420		420
421	`if ( abs(pollat-90.).gt.eps) then`	421	`if ( abs(pollat-90.).gt.eps) then`
422	`rlon = lmtolms(lat,lon,pollat,pollon)`	422	`rlon = lmtolms(lat,lon,pollat,pollon)`
423	`rlat = phtophs(lat,lon,pollat,pollon)`	423	`rlat = phtophs(lat,lon,pollat,pollon)`
424	`else`	424	`else`
425	`rlon = lon`	425	`rlon = lon`
426	`rlat = lat`	426	`rlat = lat`
427	`endif`	427	`endif`
428		428
429	`trainp(i,j,2) = rlon`	429	`trainp(i,j,2) = rlon`
430	`trainp(i,j,3) = rlat`	430	`trainp(i,j,3) = rlat`
431		431
432	`enddo`	432	`enddo`
433	`enddo`	433	`enddo`
434		434
435	`c Check that first four columns correspond to time,lon,lat,z`	435	`c Check that first four columns correspond to time,lon,lat,z`
436	`if ( (varsinp(1).ne.'time' ).or.`	436	`if ( (varsinp(1).ne.'time' ).or.`
437	`> (varsinp(2).ne.'xpos' ).and.(varsinp(2).ne.'lon' ).or.`	437	`> (varsinp(2).ne.'xpos' ).and.(varsinp(2).ne.'lon' ).or.`
438	`> (varsinp(3).ne.'ypos' ).and.(varsinp(3).ne.'lat' ).or.`	438	`> (varsinp(3).ne.'ypos' ).and.(varsinp(3).ne.'lat' ).or.`
439	`> (varsinp(4).ne.'zpos' ).and.(varsinp(4).ne.'z' ) )`	439	`> (varsinp(4).ne.'zpos' ).and.(varsinp(4).ne.'z' ) )`
440	`>then`	440	`>then`
441	`print*,' ERROR: problem with input trajectories ...'`	441	`print*,' ERROR: problem with input trajectories ...'`
442	`stop`	442	`stop`
443	`endif`	443	`endif`
444	`varsinp(1) = 'time'`	444	`varsinp(1) = 'time'`
445	`varsinp(2) = 'lon'`	445	`varsinp(2) = 'lon'`
446	`varsinp(3) = 'lat'`	446	`varsinp(3) = 'lat'`
447	`varsinp(4) = 'z'`	447	`varsinp(4) = 'z'`
448		448
449	`c Write some status information of the input trajectories`	449	`c Write some status information of the input trajectories`
450	`print*,'---- INPUT TRAJECTORIES ---------------------------------'`	450	`print*,'---- INPUT TRAJECTORIES ---------------------------------'`
451	`print*`	451	`print*`
452	`print*,' Start date : ',trim(startdate)`	452	`print*,' Start date : ',trim(startdate)`
453	`print*,' End date : ',trim(enddate)`	453	`print*,' End date : ',trim(enddate)`
454	`print*,' Reference time (year) : ',reftime(1)`	454	`print*,' Reference time (year) : ',reftime(1)`
455	`print*,' (month) : ',reftime(2)`	455	`print*,' (month) : ',reftime(2)`
456	`print*,' (day) : ',reftime(3)`	456	`print*,' (day) : ',reftime(3)`
457	`print*,' (hour) : ',reftime(4)`	457	`print*,' (hour) : ',reftime(4)`
458	`print*,' (min) : ',reftime(5)`	458	`print*,' (min) : ',reftime(5)`
459	`print*,' Time range (min) : ',reftime(6)`	459	`print*,' Time range (min) : ',reftime(6)`
460	`do i=1,ncol`	460	`do i=1,ncol`
461	`print*,' Var :',i,trim(varsinp(i))`	461	`print*,' Var :',i,trim(varsinp(i))`
462	`enddo`	462	`enddo`
463	`print*`	463	`print*`
464		464
465	`c Check that first time is 0 - otherwise the tracing will produce`	465	`c Check that first time is 0 - otherwise the tracing will produce`
466	`c wrong results because later in the code only absolute times are`	466	`c wrong results because later in the code only absolute times are`
467	`c considered: <itime0 = int(abs(tfrac-tstart)/timeinc) + 1>. This`	467	`c considered: <itime0 = int(abs(tfrac-tstart)/timeinc) + 1>. This`
468	`c will be changed in a future version.`	468	`c will be changed in a future version.`
469	`if ( abs( trainp(1,1,1) ).gt.eps ) then`	469	`if ( abs( trainp(1,1,1) ).gt.eps ) then`
470	`print*,' ERROR: First time of trajectory must be 0, i.e. '`	470	`print*,' ERROR: First time of trajectory must be 0, i.e. '`
471	`print*,' correspond to the reference date. Otherwise'`	471	`print*,' correspond to the reference date. Otherwise'`
472	`print*,' the tracing will give wrong results... STOP'`	472	`print*,' the tracing will give wrong results... STOP'`
473	`stop`	473	`stop`
474	`endif`	474	`endif`
475		475
476	`c Run a simple consistency check`	476	`c Run a simple consistency check`
477	`call do_timecheck(dat,numdat,timeinc,fbflag,reftime)`	477	`call do_timecheck(dat,numdat,timeinc,fbflag,reftime)`
478		478
479	`c --------------------------------------------------------------------`	479	`c --------------------------------------------------------------------`
480	`c Trace the fields (fields available on input files)`	480	`c Trace the fields (fields available on input files)`
481	`c --------------------------------------------------------------------`	481	`c --------------------------------------------------------------------`
482		482
483	`print*`	483	`print*`
484	`print*,'---- LIDAR FROM PRIMARY AND SECONDARY DATA FILES ------'`	484	`print*,'---- LIDAR FROM PRIMARY AND SECONDARY DATA FILES ------'`
485		485
486	`c Loop over all tracing fields`	486	`c Loop over all tracing fields`
487	`do i=1,ntrace0`	487	`do i=1,ntrace0`
488		488
489	`c Skip all fields marked for online calculation`	489	`c Skip all fields marked for online calculation`
490	`if ( compfl(i).eq.1 ) goto 110`	490	`if ( compfl(i).eq.1 ) goto 110`
491		491
492	`c Init the output fields: position and lidar field`	492	`c Init the output fields: position and lidar field`
493	`do k=1,ntim`	493	`do k=1,ntim`
494	`do l=1,nlev`	494	`do l=1,nlev`
495	`out_pos(k,l) = 0.`	495	`out_pos(k,l) = 0.`
496	`out_val(k,l) = 0.`	496	`out_val(k,l) = 0.`
497	`out_cnt(k,l) = 0.`	497	`out_cnt(k,l) = 0.`
498	`enddo`	498	`enddo`
499	`enddo`	499	`enddo`
500		500

Subversion Repositories lagranto.icon

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