WebSVN – lagranto.ecmwf – Diff – /trunk/trace/trace.f90

-PROGRAM trace_label_circle
+PROGRAM trace
   ! ********************************************************************
   ! *                                                                  *
   ! * Trace fields along trajectories                                  *
   ! *                                                                  *
-  ! * Heini Wernli       first version:       April 1993               *
+  ! * April 1993: First version (Heini Wernli)                         *
-  ! * Michael Sprenger   major upgrade:       2008-2009                *
+  ! * 2008-2009 : Major upgrades (Michael Sprenger)                    *
-  ! * Bojan Skerlak      clustering option:   March 2012               *
+  ! * Mar 2012: Clustering option (Bojan Skerlak)                      *
-  ! * Bojan Skerlak      circle option:       November 2012            *
+  ! * Nov 2012: Circle options (")                                     *
+  ! * Jul 2013: user-defined PV,TH @ clustering mode (")               *
   ! *                                                                  *
   ! ********************************************************************
   implicit none
   real :: pi180                                   ! deg -> rad
   parameter (pi180=3.14159/180.)
   real :: deg2km                                  ! deg -> km (at equator)
   parameter (deg2km=111.)
   real :: pir
-  parameter (pir=255032235.95489)                 ! 2*Pi*R^2 (Bojan)
+  parameter (pir=255032235.95489)                 ! 2*Pi*R^2 Bojan
   ! Prefix for primary and secondary fields
   character :: charp
   character :: chars
   parameter (charp='P')
   integer :: ind6,ind7,ind8,ind9,ind0
   integer :: noutside
   real :: delta
   integer :: itrace0
   character(len=80) :: string
+  integer err_c1,err_c2,err_c3
   ! Bojan
   real    :: dist,circlesum,circlemax,circlemin,circleavg,radius       ! distance (great circle), sum/max/min/avg in circle, radius of circle
   integer :: ist,jst,kst,sp,ml,mr,nd,nu                                ! ijk in stack, sp=stack counter, ml (left), mr (right), nd (down), nu (up)
   integer :: lci,lcj,xindb,xindf                                       ! label count i and j, xind back, xind forward
   integer :: yindb,yindf,pindb,pindf                                   ! yind back, yind forward, pind back, pind forward
   integer :: pvpos,thpos                                               ! position of variables PV and TH in trajectory
+  real    :: tropo_pv,tropo_th                                         ! values of PV and TH at dynamical tropopause
   integer, allocatable, dimension (:)   :: stackx,stacky               ! lon/lat of stack
   integer, allocatable, dimension (:)   :: lblcount                    ! counter for label
   integer, allocatable, dimension (:,:) :: connect                     ! array that keeps track of the visited grid points
   real, allocatable, dimension (:)      :: lbl                         ! label
   real, allocatable, dimension (:)      :: circlelon,circlelat         ! value of f, lon and lat in circle
   ! Bojan
   ! Externals
   real :: int_index4
-  external int_index4
+  external                               int_index4
-  real :: sdis
+  real :: sdis ! Bojan: need function sdis (calculates great circle distance)
-  external sdis
+  external                               sdis ! Bojan: need function sdis
   ! --------------------------------------------------------------------
   ! Start of program, Read parameters, get grid parameters
   ! --------------------------------------------------------------------
       read(10,*) svars(i)
    enddo
    read(10,*) timecheck
    read(10,*) intmode
    read(10,*) radius
+   read(10,*) tropo_pv
+   read(10,*) tropo_th
   close(10)
   ! Bojan: error if radius < 0
-  if (radius .lt. 0) then
+  if (((intmode .eq. "circle_avg") .or. (intmode .eq. "circle_min") .or. (intmode .eq. "circle_max")) .and. (radius .lt. 0)) then
      print*,'ERROR (circle): radius < 0!'
      stop
   endif
   ! Remove commented tracing fields
      allocate(lbl(8),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lbl      ***'
      allocate(lblcount(5),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lblcount ***'
   endif
   ! allocate memory for circle mode
   if ( (intmode.eq.'circle_avg') .or. (intmode.eq.'circle_min') .or. (intmode.eq.'circle_max') ) then
      allocate(connect(nx,ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating connect ***'
      allocate(stackx(nx*ny),stat=stat)
   print*,'  #tra                   : ',ntra
   print*,'  #col                   : ',ncol
   print*,'  #tim                   : ',ntim
   print*,'  No time check          : ',trim(timecheck)
   print*,'  Interpolation mode     : ',trim(intmode)
+  ! Bojan
+  if (trim(intmode) .eq. "clustering") then
+     print*,'  Tropopause PV [pvu]    : ',tropo_pv
+     print*,'  Tropopause TH [K]      : ',tropo_th
+  endif
   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)
         stop
      endif
   endif
   ! Bojan
+  ! Split the tracing variables
+  do i=1,ntrace0
+    call splitvar(tvar(i),shift_val(i),shift_dir(i) )
+  enddo
   ! Split the variable name and set flags
-  do i=1,ntrace1
+  do i=ntrace0+1,ntrace1
      ! Set the scaling factor
-     if ( i.gt.ntrace0 ) fac(i) = 1.
+     fac(i) = 1.
      ! Set the base variable name, the shift and the direction
      call splitvar(tvar(i),shift_val(i),shift_dir(i) )
-     ! Set the prefix of the file name
+     ! Set the prefix of the file name for additional fields
-     if ( ( compfl(i).eq.0 ).or.(i.gt.ntrace0) ) then
-        tfil(i)=' '
+     tfil(i)='*'
-        do j=1,n_pvars
+     do j=1,n_pvars
-           if ( tvar(i).eq.pvars(j) ) tfil(i)=charp
+        if ( tvar(i).eq.pvars(j) ) tfil(i)=charp
-        enddo
+     enddo
-        do j=1,n_svars
+     do j=1,n_svars
-           if ( tvar(i).eq.svars(j) ) tfil(i)=chars
+        if ( tvar(i).eq.svars(j) ) tfil(i)=chars
-        enddo
+     enddo
-     else
-        tfil(i) = '*'
-     endif
      ! Set the computational flag
      if ( (tvar(i).eq.'P'   ).or. &
           (tvar(i).eq.'PLAY').or. &
           (tvar(i).eq.'PLEV') ) then
       iloaded0 = -1
       iloaded1 = -1
       ! Reset the counter for fields outside domain
       noutside = 0
+      err_c1   = 0
+      err_c2   = 0
+      err_c3   = 0
       ! Loop over all times
       do j=1,ntim
          ! Convert trajectory time from hh.mm to fractional time
             else
                p0 = 1050.                                     ! Lowest level (2D tracing)
             endif
             ! Set negative pressures to mdv
-            if (p0.lt.0.) traint(k,j,4) = mdv
+            if (p0.lt.0.) then
+	        f0 = mdv
+                goto 109
+            endif
             ! Set the relative time
             call hhmm2frac(traint(k,j,1),tfrac)
             reltpos0 = fbflag * (tfrac-time0)/timeinc
                print*,'WARNING: y0>90 ',y0,' => setting to 180-y0 ',180.-y0
                y0=180.-y0
                x0=x0+per/2.
             endif
             if ((hem.eq.1).and.(y0.lt.-90.)) then
-               print*,'WARNING: y0<90 ',y0,' => setting to -180-y0 ',-180.-y0
+               print*,'WARNING: y0<-90 ',y0,' => setting to -180-y0 ',-180.-y0
                y0=-180.-y0
                x0=x0+per/2.
             endif
             ! Get the index where to interpolate (x0,y0,p0)
 		 yind = yind
 		 pind = pind
             ! Check if pressure is outside, but rest okay => adjust to lowest or highest level
             elseif ( (xind.ge.1.).and.(xind.le.real(nx)).and. (yind.ge.1.).and.(yind.le.real(ny)) ) then ! only vertical problem
                if ( pind.gt.nz ) then ! pressure too low, index too high
+                 err_c1 = err_c1 + 1
+                 if ( err_c1.lt.10 ) then
-                 print('Af5.3A'),' WARNING: pressure too low (pind = ',pind,') => adjusted to highest level (pind=nz.)'
+                    write(*,'(Af5.3A)') ' WARNING: pressure too low (pind = ',pind,') => adjusted to highest level (pind=nz.)'
-                 print*,'(x0,y0,p0)=',x0,y0,p0
+                    print*,'(x0,y0,p0)=',x0,y0,p0
-                 pind = real(nz)
+                    pind = real(nz)
+                 elseif ( err_c1.eq.10 ) then
+                    print*,' WARNING: more pressures too low -> adjusted to highest level '
+                    pind = real(nz)
+                 else
+                    pind = real(nz)
+                 endif
                elseif (pind.lt.1.) then ! pressure too high, index too low
+                 err_c2 = err_c2 + 1
+                 if ( err_c2.lt.10 ) then
-                 print('Af5.3A'),' WARNING: pressure too high (pind = ',pind,') => adjusted to lowest level, (pind=1.)'
+                    write(*,'(Af5.3A)') ' WARNING: pressure too high (pind = ',pind,') => adjusted to lowest level, (pind=1.)'
-                 print*,'(x0,y0,p0)=',x0,y0,p0
+                    print*,'(x0,y0,p0)=',x0,y0,p0
-                 pind = 1.
+                    pind = 1.
+                 elseif ( err_c2.eq.10 ) then
+                    print*,' WARNING: more pressures too high -> adjusted to lowest level '
+                    pind = 1.
+                 else
+                    pind = 1.
+                 endif
               endif
+            ! Grid point is outside!
             else
+               err_c3 = err_c3 + 1
+               if ( err_c3.lt.10 ) then
-               print*,'ERROR: point is outside grid (horizontally)'
+                  print*,'ERROR: point is outside grid (horizontally)'
+                  print*,'   Trajectory # ',k
+                  print*,'   Position     ',x0,y0,p0
-               print*,'Grid (nx,ny): ',nx,ny,' and (xind,yind): ',xind,yind
+                  print*,'  (xind,yind):  ',xind,yind
+                  xind          = mdv
+                  yind          = mdv
+                  pind          = mdv
+                  traint(k,j,2) = mdv
+                  traint(k,j,3) = mdv
+                  traint(k,j,4) = mdv
+               elseif ( err_c3.eq.10 ) then
+                  print*,'ERROR: more points outside grid (horizontally)'
+                  xind          = mdv
+                  yind          = mdv
+                  pind          = mdv
+                  traint(k,j,2) = mdv
+                  traint(k,j,3) = mdv
+                  traint(k,j,4) = mdv
-               stop
+               else
+                  xind          = mdv
+                  yind          = mdv
+                  pind          = mdv
+                  traint(k,j,2) = mdv
+                  traint(k,j,3) = mdv
+                  traint(k,j,4) = mdv
-            endif
+               endif
-            ! ------------------------ NEAREST mode -------------------------------
+            endif
+            ! ------------------------ NEAREST mode -------------------------------
+            ! Interpolate to nearest grid point
             if ( intmode.eq.'nearest') then
                xind = real( nint(xind) )
                yind = real( nint(yind) )
                pind = real( nint(pind) )
                endif
                ! interpolate
                f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
-            ! ------------------------ CLUSTERING mode ------------------------------- Bojan
+            ! ------------------------ end NEAREST mode -------------------------------
+            ! ------------------------ CLUSTERING mode ------------------------------- Bojan
             elseif (intmode.eq.'clustering') then
-            ! Clustering mode only works with LABEL
-            if (varname.ne.'LABEL' ) then
+               if (varname.ne.'LABEL' ) then
-               print*,'ERROR (clustering): varname is not LABEL'
+                  print*,'ERROR (clustering): varname is not LABEL'
-               stop
+                  stop
-            endif
+               endif
                   pindb=nz-1
                else
                   pindf=pindb+1
                endif
             endif
             ! Give warnings and stop if problems occur
             if ( xindb.eq.xindf ) then
                print*,'ERROR (clustering): xindb=xindf'
                print*,xind,xindb,xindf
                stop
             ! Set to -9 to detect if no label was assigned in the end
             f0=-9
             ! Stratosphere (PV)
-            if ( abs(traint(k,j,pvpos)) .ge. 2.0 ) then
+            if ( abs(traint(k,j,pvpos)) .ge. tropo_pv ) then
                if ( (lblcount(2).ge.lblcount(3)).and. (lblcount(2).ge.lblcount(5)) ) then
                   f0=2
                elseif ( lblcount(3).ge.lblcount(5) ) then
                   f0=3
                elseif ( lblcount(5).gt.lblcount(3) ) then
                   f0=5
                endif
             endif
             ! Troposphere (PV)
-            if ( abs(traint(k,j,pvpos)) .lt. 2.0 ) then
+            if ( abs(traint(k,j,pvpos)) .lt. tropo_pv ) then
                if ( lblcount(1).ge.lblcount(4) ) then
                f0=1
                elseif ( lblcount(4).gt.lblcount(1) ) then
                f0=4
                endif
             endif
             ! Stratosphere (TH)
-            if ( traint(k,j,thpos) .ge. 380.0 ) then
+            if ( traint(k,j,thpos) .ge. tropo_th ) then
                f0=2
             endif
             if (f0.eq.-9) then
                print*,'ERROR (Clustering): No label assigned!'
                stop
             endif
+            ! ------------------------ end CLUSTERING mode -------------------------------
             ! ------------------------ CIRCLE modes ------------------------------- Bojan
+            ! elseif (not clustering but one of the possible circle modes)
             elseif ( (intmode.eq.'circle_avg') .or. (intmode.eq.'circle_min') .or. (intmode.eq.'circle_max') ) then
             ! reset arrays for this point
-            connect    = 0
+            connect=0
-            stackx     = 0
+            stackx=0
             ! Initialize stack with nint(xind),nint(yind)
-            kst       = 0            ! counts the number of points in circle
+            kst=0 ! counts the number of points in circle
-            stackx(1) = nint(xind)
+            stackx(1)=nint(xind)
-            stacky(1) = nint(yind)
+            stacky(1)=nint(yind)
-            sp        = 1            ! stack counter
+            sp=1 ! stack counter
-            ! Repeat as long as the stack is not empty
             do while (sp.ne.0)
             ! Get an element from stack
              ist=stackx(sp)
              jst=stacky(sp)
              endif ! endif radius is smaller => end of updating stack
             end do ! end working on stack
-	    ! Choose output depending on intmode
             if (kst.ge.1) then
-               ! calculate area-weighted average of f in circle
+               ! Choose output depending on intmode
                if ( intmode .eq. 'circle_avg' ) then
+                  ! calculate area-weighted average of f in circle
                   circlesum=0.
                   do l=1,kst
                      circlesum=circlesum+circlef(l)*circlearea(l)
                   enddo
                   circleavg=circlesum/sum(circlearea(1:kst))
+                  !print*,'area-weighted average of f in circle=',circleavg
                   f0=circleavg
-               ! calculate minimum in circle
                elseif ( intmode .eq. 'circle_min' ) then
+                  ! calculate minimum in circle
                   circlemin=circlef(1)
                   do l=1,kst
                      if (circlef(l) .lt. circlemin) then
                         circlemin=circlef(l)
                      endif
                   enddo
+                  !print*,'minimum of f in circle=',circlemin
                   f0=circlemin
-               ! calculate maximum in circle
                elseif ( intmode .eq. 'circle_max' ) then
+                  ! calculate maximum in circle
                   circlemax=circlef(1)
                   do l=1,kst
                      if (circlef(l) .gt. circlemax) then
                         circlemax=circlef(l)
                      endif
                   enddo
+                  !print*,'maximum of f in circle=',circlemax
                   f0=circlemax
                else
                   print*,'ERROR (circle): intmode not valid!'
                   stop
                endif
             else
                print*,'ERROR (circle): Search radius is too small... (2). r =',radius
                stop
             endif
+            ! ------------------------ end CIRCLE modes -------------------------------
             ! ------------------------ NORMAL mode -------------------------------
-            else
+            else ! not clustering nor circle: NORMAL mode
+            ! Check if point is within grid
+!            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
+!
+            ! Do the interpolation: everthing is ok
                f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
+!            ! Check if pressure is outside, but rest okay: adjust to lowest or highest level
+!            elseif ( (xind.ge.1.).and.(xind.le.real(nx)).and. (yind.ge.1.).and.(yind.le.real(ny)) ) then ! only vertical problem
+!               if ( pind.gt.nz ) then ! pressure too low, index too high
+!                 pind = real(nz)
+!                 print*,' Warning: pressure too low -> adjusted to highest level, pind=nz.'
+!                 print*,'(x0,y0,p0)=',x0,y0,p0
+!               elseif (pind.lt.1.) then ! pressure too high, index too low
+!                 pind = 1.
+!                 print*,' Warning: pressure too high -> adjusted to lowest level, pind=1.'
+!                 print*,'(x0,y0,p0)=',x0,y0,p0
+!              endif
+!              f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
+!            ! Less than 10 outside
+!            elseif (noutside.lt.10) then
+!               print*,' ',trim(tvar(i)),' @ ',x0,y0,p0,'outside'
+!               f0       = mdv
+!               noutside = noutside + 1
+!
+!            ! More than 10 outside
+!            elseif (noutside.eq.10) then
+!               print*,' ...more than 10 outside...'
+!               f0       = mdv
+!               noutside = noutside + 1
+!            ! Else (not everything okay and also not 'tolerated cases') set to missing data
+!            else
+!               f0       = mdv
-            endif
+!            endif
+            ! ------------------------ end NORMAL mode -------------------------------
+            endif ! end if nearest case
+           ! Exit for loop over all trajectories and times -save interpolated value
+continue
             ! Save the new field
             if ( abs(f0-mdv).gt.eps) then
                traint(k,j,ncol+i) = f0
             else
   print*
   print*,'              *** END OF PROGRAM TRACE ***'
   print*,'========================================================='
-end program trace_label_circle
+end program trace
 ! ******************************************************************
 ! * SUBROUTINE SECTION                                             *

Subversion Repositories lagranto.ecmwf

(root)/trunk/trace/trace.f90 – Rev 3 → 5