Subversion Repositories lagranto.wrf

PROGRAM trace_label_circle

  ! ********************************************************************
  ! *                                                                  *
  ! * Trace fields along trajectories                                  *
  ! *                                                                  *
  ! * Heini Wernli       first version:       April 1993               *
  ! * Michael Sprenger   major upgrade:       2008-2009                *
  ! * Sebastian Schemm WRF                                             *
  ! ********************************************************************

  implicit none

  ! --------------------------------------------------------------------
  ! Declaration of parameters
  ! --------------------------------------------------------------------

  ! Maximum number of levels for input files
  integer :: nlevmax
  parameter (nlevmax=100)

  ! Maximum number of input files (dates, length of trajectories)
  integer :: ndatmax
  parameter (ndatmax=500)

  ! Numerical epsilon (for float comparison)
  real :: eps
  parameter (eps=0.001)
  real         mdv                              ! missing data value
  parameter    (mdv=-999.)

  ! Prefix for primary and secondary fields
  character :: charp
  character :: chars
  parameter (charp='P')
  parameter (chars='S')

  ! --------------------------------------------------------------------
  ! Declaration of variables
  ! --------------------------------------------------------------------

  ! Input and output format for trajectories (see iotra.f)
  integer :: inpmode
  integer :: outmode

  ! Input parameters
  character(len=80) :: inpfile         ! Input trajectory file
  character(len=80) :: outfile         ! Output trajectory file
  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(len=80) :: tvar0(200)      ! Tracing variable (with mode specification)
  character(len=80) :: tvar(200)       ! Tracing variable name (only the variable)
  character(len=1)  :: tfil(200)       ! Filename prefix
  real :: fac(200)                     ! Scaling factor
  integer :: compfl(200)               ! Computation flag (1=compute)
  integer :: numdat                    ! Number of input files
  character(len=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(len=20) :: startdate       ! First time/date on trajectory
  character(len=20) :: enddate         ! Last time/date on trajectory
  character(len=80) :: timecheck       ! Either 'yes' or 'no'
  character(len=80) :: intmode         ! Interpolation mode ('normal', 'nearest')

  ! Trajectories
  real,allocatable, dimension (:,:,:) :: trainp          ! Input trajectories (ntra,ntim,ncol)
  real,allocatable, dimension (:,:,:) :: traint          ! Internal trajectories (ntra,ntim,ncol+ntrace0)
  real,allocatable, dimension (:,:,:) :: traout          ! Output trajectories (ntra,ntim,ncol+ntrace0)
  integer :: reftime(6)                                  ! Reference date
  character(len=80) :: varsinp(100)                      ! Field names for input trajectory
  character(len=80) :: varsint(100)                      ! Field names for internal trajectory
  character(len=80) :: varsout(100)                      ! Field names for output trajectory
  integer :: fid,fod                                     ! File identifier for inp and out trajectories
  real :: x0,y0,p0                                       ! Position of air parcel (physical space)
  real :: reltpos0                                       ! Relative time of air parcel
  real :: xind,yind,pind                                 ! Position of air parcel (grid space)
  integer :: fbflag                                      ! Flag for forward (1) or backward (-1) trajectories
  integer :: fok(100)                                    ! Flag whether field is ready

  ! Meteorological fields
  real,allocatable, dimension (:)     :: spt0,spt1       ! Surface pressure
  real,allocatable, dimension (:)     :: p3t0,p3t1       ! 3d-pressure
  real,allocatable, dimension (:)     :: f3t0,f3t1       ! 3d field for tracing
  character(len=80) :: svars(100)                        ! List of variables on S file
  character(len=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

  ! Grid description
  real :: pollon,pollat   ! Longitude/latitude of pole
  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

  ! Auxiliary variables
  integer :: i,j,k,l,m,n
  real :: rd
  character(len=80) :: filename,varname
  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(len=80) :: string

  ! Externals
  real :: int_index4
  external int_index4
  real :: sdis   
  external sdis 

  ! --------------------------------------------------------------------
  ! Start of program, Read parameters, get grid parameters
  ! --------------------------------------------------------------------

  ! Write start message
  print*,'========================================================='
  print*,'              *** START OF PROGRAM TRACE ***'
  print*

  ! Read parameters
  open(10,file='trace.param')
   read(10,*) inpfile
   read(10,*) outfile
   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
  close(10)

  ! 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

  ! Save the tracing variable (including all mode specifications)
  do i=1,ntrace0
     tvar0(i) = tvar(i)
  enddo

  ! Set the formats of the input and output files
  call mode_tra(inpmode,inpfile)
  if (inpmode.eq.-1) inpmode=1
  call mode_tra(outmode,outfile)
  if (outmode.eq.-1) outmode=1

  ! Convert time shifts <tst,ten> from <hh.mm> into fractional time
  call hhmm2frac(tst,frac)
  tst = frac
  call hhmm2frac(ten,frac)
  ten = frac

  ! Set the time for the first data file (depending on forward/backward mode)
  if (fbflag.eq.1) then
    tstart = -tst
  else
    tstart = tst
  endif
 
  ! Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)
  ! 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_wrf(fid,xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz)
  call input_close(fid)

  ! 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 ***'   ! Tracing field
  allocate(f3t1(nx*ny*nz),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array p3t1 ***'

  ! Get memory for trajectory arrays
  allocate(trainp(ntra,ntim,ncol),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array tra      ***'

  ! Write some status information
  print*,'---- INPUT PARAMETERS -----------------------------------'
  print*
  print*,'  Input trajectory file  : ',trim(inpfile)
  print*,'  Output trajectory file : ',trim(outfile)
  print*,'  Format of input file   : ',inpmode
  print*,'  Format of output file  : ',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)), fac(i), '  1 ', tfil(i)
     endif
  enddo
  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*

  ! --------------------------------------------------------------------
  ! Load the input trajectories
  ! --------------------------------------------------------------------

  ! 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)
 
  ! Check that first four columns correspond to time,x,y,z
  if ( (varsinp(1).ne.'time' ).or. &
       (varsinp(2).ne.'x'    ).or. &
       (varsinp(3).ne.'y'    ).or. &
       (varsinp(4).ne.'z'   ) )    &
  then
     print*,' ERROR: problem with input trajectories ...'
     stop
  endif
  varsinp(1) = 'time'
  varsinp(2) = 'x'
  varsinp(3) = 'y'
  varsinp(4) = 'z'

  ! Transform start coordinates from index space to WRF grid
  do i=1,ntra
     do j=1,ntim
         if ( ( abs(trainp(i,j,2)-mdv).gt.eps).and.( abs(trainp(i,j,3)-mdv).gt.eps) ) then
            trainp(i,j,2) = xmin + ( trainp(i,j,2) - 1. ) * dx
            trainp(i,j,3) = ymin + ( trainp(i,j,3) - 1. ) * dy
         endif
     enddo
  enddo

  ! 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*

  ! Check that first time is 0 - otherwise the tracing will produce
  ! wrong results because later in the code only absolute times are
  ! considered: <itime0   = int(abs(tfrac-tstart)/timeinc) + 1>. This
  ! 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

  ! Save the full variable name (including shift specification)
  do i=1,ncol
     varsint(i)      = varsinp(i)
  enddo
  do i=1,ntrace0
     varsint(i+ncol) = tvar(i)
  enddo

  ! Split the variable name and set flags
  do i=1,ntrace0

     ! Set the prefix of the file name
     if ( compfl(i).eq.0  ) then
        tfil(i)=' '
        do j=1,n_pvars
           if ( tvar(i).eq.pvars(j) ) tfil(i)=charp
        enddo
        do j=1,n_svars
           if ( tvar(i).eq.svars(j) ) tfil(i)=chars
        enddo
     elseif ( tvar(i).eq.'Z' ) then
        tfil(i) = 'P'
     elseif ( tvar(i).eq.'ZB' ) then
        tfil(i) = 'P'
     else
        tfil(i) = '*' 
     endif

  enddo

  ! Write status information
  print*
  print*,'---- COMPLETE TABLE FOR TRACING -------------------------'
  print*
  do i=1,ntrace0
      write(*,'(i4,a4,a8,10x,8x,3x,a1,i5)') &
             i,' : ',trim(tvar(i)),tfil(i),compfl(i)
  enddo

  ! Allocate memory for internal trajectories
  allocate(traint(ntra,ntim,ncol+ntrace0),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array traint   ***'

  ! Copy input to output trajectory
  do i=1,ntra
     do j=1,ntim
        do k=1,ncol
           traint(i,j,k)=trainp(i,j,k)
        enddo
     enddo
  enddo

  ! --------------------------------------------------------------------
  ! Trace the fields (fields available on input files)
  ! --------------------------------------------------------------------

  print*
  print*,'---- TRACING FROM PRIMARY AND SECONDARY DATA FILES ------'

  ! Loop over all tracing fields
  do i=1,ntrace0

      ! Write some status information
      print*
      print*,' Now tracing             : ', trim(tvar(i)),' ',trim(tfil(i))

      ! Reset flags for load manager
      iloaded0 = -1
      iloaded1 = -1

      ! Reset the counter for fields outside domain
      noutside = 0

      ! Loop over all times
      do j=1,ntim

         ! Convert trajectory time from hh.mm to fractional time
         call hhmm2frac(trainp(1,j,1),tfrac)

         ! 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

         ! Load manager: Check whether itime0 can be copied from itime1
         if ( itime0.eq.iloaded1 ) then
            f3t0     = f3t1
            p3t0     = p3t1
            spt0     = spt1
            iloaded0 = itime0
         endif

         ! Load manager: Check whether itime1 can be copied from itime0
         if ( itime1.eq.iloaded0 ) then
            f3t1     = f3t0
            p3t1     = p3t0
            spt1     = spt0
            iloaded1 = itime1
         endif


         ! 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_var_wrf(fid,varname,f3t0,nx,ny,nz)   ! field to trace

                varname='geopot'                                                    ! geopot.height
                call input_var_wrf(fid,varname,p3t0,nx,ny,nz)

                varname='geopots'                                                   ! surface geopot.height
                call input_var_wrf(fid,varname,spt0,nx,ny,1 )

            call input_close(fid)

            iloaded0 = itime0

         endif

         ! 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_var_wrf(fid,varname,f3t1,nx,ny,nz)               ! field to trace

                varname='geopot'                                                            ! geopot.height
                call input_var_wrf(fid,varname,p3t1,nx,ny,nz)

                varname='geopots'                                                           ! surface geopot.height 
                call input_var_wrf(fid,varname,spt1,nx,ny,1)

            call input_close(fid)

            iloaded1 = itime1

         endif

         ! Loop over all trajectories
         do k=1,ntra

            ! Set the horizontal position where to interpolate to
            x0       = traint(k,j,2)                          ! Longitude
            y0       = traint(k,j,3)
                                   ! Latitude
            ! Set the relative time
            call hhmm2frac(traint(k,j,1),tfrac)
            reltpos0 = fbflag * (tfrac-time0)/timeinc

            ! Set the vertical position where to interpolate to: 
            if ( nz.gt.1 ) then
               p0 = traint(k,j,4)                            ! Pressure (3D tracing)
            else
               p0 = 1                                        ! Lowest level (2D tracing)
            endif

            ! 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


            ! Check if point is within grid (keep indices if 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
                           xind = xind
                           yind = yind
                           pind = pind

            endif

            ! ------------------------ NEAREST mode ------------------------------- 

            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)

               if ( abs(reltpos0).ge.eps ) then
                  print*,'ERROR (nearest): reltpos != 0',reltpos0
                  stop
               endif

               ! interpolate
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)

            ! ------------------------ NORMAL mode -------------------------------
            else 

               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
               
            endif

            ! Save the new field
            if ( abs(f0-mdv).gt.eps) then
               traint(k,j,ncol+i) = f0
            else
               traint(k,j,ncol+i) = mdv
            endif

         enddo ! end loop over all trajectories

      enddo ! end loop over all times

      ! Exit point for loop over all tracing variables
 110   continue

  enddo ! end loop over all tracing variables

  ! --------------------------------------------------------------------
  ! Write output to output trajectory file
  ! --------------------------------------------------------------------

  ! Write status information
  print*
  print*,'---- WRITE OUTPUT TRAJECTORIES --------------------------'
  print*

  ! Allocate memory for internal trajectories
  allocate(traout(ntra,ntim,ncol+ntrace0),stat=stat)
  if (stat.ne.0) print*,'*** error allocating array traout   ***'

  ! Copy input to output trajectory (apply scaling of output)
  do i=1,ntra
     do j=1,ntim
        do k=1,ncol+ntrace0
           if ( k.le.ncol ) then
              traout(i,j,k) = traint(i,j,k)
           elseif ( abs(traint(i,j,k)-mdv).gt.eps ) then
              traout(i,j,k) = fac(k-ncol) * traint(i,j,k)
           else
              traout(i,j,k) = mdv
           endif
        enddo
     enddo
  enddo

  ! Set the variable names for output trajectory
  do i=1,ncol+ntrace0
     varsout(i)      = varsint(i)
  enddo

 ! Transform coordinates into index space
   do i=1,ntra
     do j=1,ntim

         if ( ( abs(traout(i,j,2)-mdv).gt.eps ).and.( abs(traout(i,j,3)-mdv).gt.eps ) ) then
            traout(i,j,2) = ( traout(i,j,2) - xmin ) / dx + 1.
            traout(i,j,3) = ( traout(i,j,3) - ymin ) / dy + 1.
         endif

     enddo
  enddo

  ! Write trajectories
  call wopen_tra(fod,outfile,ntra,ntim,ncol+ntrace0, &
       reftime,varsout,outmode)
  call write_tra(fod,traout ,ntra,ntim,ncol+ntrace0,outmode)
  call close_tra(fod,outmode)

  ! Write some status information, and end of program message
  print*
  print*,'---- STATUS INFORMATION --------------------------------'
  print*
  print*,' ok'
  print*
  print*,'              *** END OF PROGRAM TRACE ***'
  print*,'========================================================='


end program trace_label_circle