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