Subversion Repositories lagranto.ecmwf

Rev

Rev 3 | Rev 15 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 3 Rev 5
Line 1... Line 1...
1 
PROGRAM trace_label_circle
 1 
PROGRAM trace
2 
 
 2 
 
3 
  ! ********************************************************************
 3 
  ! ********************************************************************
4 
  ! *                                                                  *
 4 
  ! *                                                                  *
5 
  ! * Trace fields along trajectories                                  *
 5 
  ! * Trace fields along trajectories                                  *
6 
  ! *                                                                  *
 6 
  ! *                                                                  *
7 
  ! * Heini Wernli       first version:       April 1993               *
 7 
  ! * April 1993: First version (Heini Wernli)                         *
8 
  ! * Michael Sprenger   major upgrade:       2008-2009                *
 8 
  ! * 2008-2009 : Major upgrades (Michael Sprenger)                    *
9 
  ! * Bojan Skerlak      clustering option:   March 2012               *
 9 
  ! * Mar 2012: Clustering option (Bojan Skerlak)                      *
10 
  ! * Bojan Skerlak      circle option:       November 2012            *
 10 
  ! * Nov 2012: Circle options (")                                     *
- 
 
 11 
  ! * Jul 2013: user-defined PV,TH @ clustering mode (")               *
11 
  ! *                                                                  *
 12 
  ! *                                                                  *
12 
  ! ********************************************************************
 13 
  ! ********************************************************************
13 
 
 14 
 
14 
  implicit none
 15 
  implicit none
15 
 
 16 
 
Line 33... Line 34...
33 
  real :: pi180                                   ! deg -> rad
 34 
  real :: pi180                                   ! deg -> rad
34 
  parameter (pi180=3.14159/180.)
 35 
  parameter (pi180=3.14159/180.)
35 
  real :: deg2km                                  ! deg -> km (at equator)
 36 
  real :: deg2km                                  ! deg -> km (at equator)
36 
  parameter (deg2km=111.)
 37 
  parameter (deg2km=111.)
37 
  real :: pir
 38 
  real :: pir
38 
  parameter (pir=255032235.95489)                 ! 2*Pi*R^2 (Bojan)
 39 
  parameter (pir=255032235.95489)                 ! 2*Pi*R^2 Bojan
39 
 
 40 
 
40 
  ! Prefix for primary and secondary fields
 41 
  ! Prefix for primary and secondary fields
41 
  character :: charp
 42 
  character :: charp
42 
  character :: chars
 43 
  character :: chars
43 
  parameter (charp='P')
 44 
  parameter (charp='P')
Line 132... Line 133...
132 
  integer :: ind6,ind7,ind8,ind9,ind0
 133 
  integer :: ind6,ind7,ind8,ind9,ind0
133 
  integer :: noutside
 134 
  integer :: noutside
134 
  real :: delta
 135 
  real :: delta
135 
  integer :: itrace0
 136 
  integer :: itrace0
136 
  character(len=80) :: string
 137 
  character(len=80) :: string
- 
 
 138 
  integer err_c1,err_c2,err_c3
137 
 
 139 
 
138 
  ! Bojan
 140 
  ! Bojan
139 
  real    :: dist,circlesum,circlemax,circlemin,circleavg,radius       ! distance (great circle), sum/max/min/avg in circle, radius of circle
 141 
  real    :: dist,circlesum,circlemax,circlemin,circleavg,radius       ! distance (great circle), sum/max/min/avg in circle, radius of circle
140 
  integer :: ist,jst,kst,sp,ml,mr,nd,nu                                ! ijk in stack, sp=stack counter, ml (left), mr (right), nd (down), nu (up)
 142 
  integer :: ist,jst,kst,sp,ml,mr,nd,nu                                ! ijk in stack, sp=stack counter, ml (left), mr (right), nd (down), nu (up)
141 
  integer :: lci,lcj,xindb,xindf                                       ! label count i and j, xind back, xind forward
 143 
  integer :: lci,lcj,xindb,xindf                                       ! label count i and j, xind back, xind forward
142 
  integer :: yindb,yindf,pindb,pindf                                   ! yind back, yind forward, pind back, pind forward
 144 
  integer :: yindb,yindf,pindb,pindf                                   ! yind back, yind forward, pind back, pind forward
143 
  integer :: pvpos,thpos                                               ! position of variables PV and TH in trajectory
 145 
  integer :: pvpos,thpos                                               ! position of variables PV and TH in trajectory
- 
 
 146 
  real    :: tropo_pv,tropo_th                                         ! values of PV and TH at dynamical tropopause
144 
  integer, allocatable, dimension (:)   :: stackx,stacky               ! lon/lat of stack
 147 
  integer, allocatable, dimension (:)   :: stackx,stacky               ! lon/lat of stack
145 
  integer, allocatable, dimension (:)   :: lblcount                    ! counter for label
 148 
  integer, allocatable, dimension (:)   :: lblcount                    ! counter for label
146 
  integer, allocatable, dimension (:,:) :: connect                     ! array that keeps track of the visited grid points
 149 
  integer, allocatable, dimension (:,:) :: connect                     ! array that keeps track of the visited grid points
147 
  real, allocatable, dimension (:)      :: lbl                         ! label
 150 
  real, allocatable, dimension (:)      :: lbl                         ! label
148 
  real, allocatable, dimension (:)      :: circlelon,circlelat         ! value of f, lon and lat in circle
 151 
  real, allocatable, dimension (:)      :: circlelon,circlelat         ! value of f, lon and lat in circle
Line 150... Line 153...
150 
  real, allocatable, dimension (:)      :: longrid,latgrid             ! arrays of longitude and latitude of grid points
 153 
  real, allocatable, dimension (:)      :: longrid,latgrid             ! arrays of longitude and latitude of grid points
151 
  ! Bojan
 154 
  ! Bojan
152 
 
 155 
 
153 
  ! Externals
 156 
  ! Externals
154 
  real :: int_index4
 157 
  real :: int_index4
155 
  external int_index4
 158 
  external                               int_index4
156 
  real :: sdis
159 
  real :: sdis ! Bojan: need function sdis (calculates great circle distance)
157 
  external sdis
160 
  external                               sdis ! Bojan: need function sdis
158 
 
 161 
 
159 
  ! --------------------------------------------------------------------
 162 
  ! --------------------------------------------------------------------
160 
  ! Start of program, Read parameters, get grid parameters
 163 
  ! Start of program, Read parameters, get grid parameters
161 
  ! --------------------------------------------------------------------
 164 
  ! --------------------------------------------------------------------
162 
 
 165 
 
Line 201... Line 204...
201 
      read(10,*) svars(i)
 204 
      read(10,*) svars(i)
202 
   enddo
 205 
   enddo
203 
   read(10,*) timecheck
 206 
   read(10,*) timecheck
204 
   read(10,*) intmode
 207 
   read(10,*) intmode
205 
   read(10,*) radius
 208 
   read(10,*) radius
- 
 
 209 
   read(10,*) tropo_pv
- 
 
 210 
   read(10,*) tropo_th
206 
  close(10)
 211 
  close(10)
207 
 
 212 
 
208 
  ! Bojan: error if radius < 0
 213 
  ! Bojan: error if radius < 0
209 
  if (radius .lt. 0) then
 214 
  if (((intmode .eq. "circle_avg") .or. (intmode .eq. "circle_min") .or. (intmode .eq. "circle_max")) .and. (radius .lt. 0)) then
210 
     print*,'ERROR (circle): radius < 0!'
 215 
     print*,'ERROR (circle): radius < 0!'
211 
     stop
 216 
     stop
212 
  endif
 217 
  endif
213 
 
 218 
 
214 
  ! Remove commented tracing fields
 219 
  ! Remove commented tracing fields
Line 284... Line 289...
284 
     allocate(lbl(8),stat=stat)
 289 
     allocate(lbl(8),stat=stat)
285 
     if (stat.ne.0) print*,'*** error allocating array lbl      ***'
 290 
     if (stat.ne.0) print*,'*** error allocating array lbl      ***'
286 
     allocate(lblcount(5),stat=stat)
 291 
     allocate(lblcount(5),stat=stat)
287 
     if (stat.ne.0) print*,'*** error allocating array lblcount ***'
 292 
     if (stat.ne.0) print*,'*** error allocating array lblcount ***'
288 
  endif
 293 
  endif
289 
 
 - 
 
290 
  ! allocate memory for circle mode
 294 
  ! allocate memory for circle mode
291 
  if ( (intmode.eq.'circle_avg') .or. (intmode.eq.'circle_min') .or. (intmode.eq.'circle_max') ) then
 295 
  if ( (intmode.eq.'circle_avg') .or. (intmode.eq.'circle_min') .or. (intmode.eq.'circle_max') ) then
292 
     allocate(connect(nx,ny),stat=stat)
 296 
     allocate(connect(nx,ny),stat=stat)
293 
     if (stat.ne.0) print*,'*** error allocating connect ***'
 297 
     if (stat.ne.0) print*,'*** error allocating connect ***'
294 
     allocate(stackx(nx*ny),stat=stat)
 298 
     allocate(stackx(nx*ny),stat=stat)
Line 350... Line 354...
350 
  print*,'  #tra                   : ',ntra
 354 
  print*,'  #tra                   : ',ntra
351 
  print*,'  #col                   : ',ncol
 355 
  print*,'  #col                   : ',ncol
352 
  print*,'  #tim                   : ',ntim
 356 
  print*,'  #tim                   : ',ntim
353 
  print*,'  No time check          : ',trim(timecheck)
 357 
  print*,'  No time check          : ',trim(timecheck)
354 
  print*,'  Interpolation mode     : ',trim(intmode)
 358 
  print*,'  Interpolation mode     : ',trim(intmode)
- 
 
 359 
  ! Bojan
- 
 
 360 
  if (trim(intmode) .eq. "clustering") then
- 
 
 361 
     print*,'  Tropopause PV [pvu]    : ',tropo_pv
- 
 
 362 
     print*,'  Tropopause TH [K]      : ',tropo_th
- 
 
 363 
  endif
355 
  do i=1,ntrace0
 364 
  do i=1,ntrace0
356 
     if (compfl(i).eq.0) then
 365 
     if (compfl(i).eq.0) then
357 
        print*,'  Tracing field          : ', trim(tvar(i)), fac(i), ' 0 ', tfil(i)
 366 
        print*,'  Tracing field          : ', trim(tvar(i)), fac(i), ' 0 ', tfil(i)
358 
     else
 367 
     else
359 
        print*,'  Tracing field          : ', trim(tvar(i)), fac(i), '  1 ', tfil(i)
 368 
        print*,'  Tracing field          : ', trim(tvar(i)), fac(i), '  1 ', tfil(i)
Line 748... Line 757...
748 
        stop
 757 
        stop
749 
     endif
 758 
     endif
750 
  endif
 759 
  endif
751 
  ! Bojan
 760 
  ! Bojan
752 
 
 761 
 
- 
 
 762 
  ! Split the tracing variables
- 
 
 763 
  do i=1,ntrace0
- 
 
 764 
    call splitvar(tvar(i),shift_val(i),shift_dir(i) )
- 
 
 765 
  enddo
- 
 
 766 
 
- 
 
 767 
 
753 
  ! Split the variable name and set flags
 768 
  ! Split the variable name and set flags
754 
  do i=1,ntrace1
 769 
  do i=ntrace0+1,ntrace1
755 
 
 770 
 
756 
     ! Set the scaling factor
 771 
     ! Set the scaling factor
757 
     if ( i.gt.ntrace0 ) fac(i) = 1.
 772 
     fac(i) = 1.
758 
 
 773 
 
759 
     ! Set the base variable name, the shift and the direction
 774 
     ! Set the base variable name, the shift and the direction
760 
     call splitvar(tvar(i),shift_val(i),shift_dir(i) )
 775 
     call splitvar(tvar(i),shift_val(i),shift_dir(i) )
761 
 
 776 
 
762 
     ! Set the prefix of the file name
 777 
     ! Set the prefix of the file name for additional fields
763 
     if ( ( compfl(i).eq.0 ).or.(i.gt.ntrace0) ) then
 - 
 
764 
        tfil(i)=' '
 778 
     tfil(i)='*'
765 
        do j=1,n_pvars
 779 
     do j=1,n_pvars
766 
           if ( tvar(i).eq.pvars(j) ) tfil(i)=charp
 780 
        if ( tvar(i).eq.pvars(j) ) tfil(i)=charp
767 
        enddo
 781 
     enddo
768 
        do j=1,n_svars
 782 
     do j=1,n_svars
769 
           if ( tvar(i).eq.svars(j) ) tfil(i)=chars
 783 
        if ( tvar(i).eq.svars(j) ) tfil(i)=chars
770 
        enddo
 - 
 
771 
     else
 - 
 
772 
        tfil(i) = '*'
 - 
 
773 
     endif
 784 
     enddo
774 
 
 785 
 
775 
     ! Set the computational flag
 786 
     ! Set the computational flag
776 
     if ( (tvar(i).eq.'P'   ).or. &
 787 
     if ( (tvar(i).eq.'P'   ).or. &
777 
          (tvar(i).eq.'PLAY').or. &
 788 
          (tvar(i).eq.'PLAY').or. &
778 
          (tvar(i).eq.'PLEV') ) then
 789 
          (tvar(i).eq.'PLEV') ) then
Line 866... Line 877...
866 
      iloaded0 = -1
 877 
      iloaded0 = -1
867 
      iloaded1 = -1
 878 
      iloaded1 = -1
868 
 
 879 
 
869 
      ! Reset the counter for fields outside domain
 880 
      ! Reset the counter for fields outside domain
870 
      noutside = 0
 881 
      noutside = 0
- 
 
 882 
      err_c1   = 0
- 
 
 883 
      err_c2   = 0
- 
 
 884 
      err_c3   = 0
871 
 
 885 
 
872 
      ! Loop over all times
 886 
      ! Loop over all times
873 
      do j=1,ntim
 887 
      do j=1,ntim
874 
 
 888 
 
875 
         ! Convert trajectory time from hh.mm to fractional time
 889 
         ! Convert trajectory time from hh.mm to fractional time
Line 963... Line 977...
963 
            else
 977 
            else
964 
               p0 = 1050.                                     ! Lowest level (2D tracing)
 978 
               p0 = 1050.                                     ! Lowest level (2D tracing)
965 
            endif
 979 
            endif
966 
 
 980 
 
967 
            ! Set negative pressures to mdv
 981 
            ! Set negative pressures to mdv
968 
            if (p0.lt.0.) traint(k,j,4) = mdv
 982 
            if (p0.lt.0.) then
- 
 
 983 
	        f0 = mdv
- 
 
 984 
                goto 109
- 
 
 985 
            endif
969 
 
 986 
 
970 
            ! Set the relative time
 987 
            ! Set the relative time
971 
            call hhmm2frac(traint(k,j,1),tfrac)
 988 
            call hhmm2frac(traint(k,j,1),tfrac)
972 
            reltpos0 = fbflag * (tfrac-time0)/timeinc
 989 
            reltpos0 = fbflag * (tfrac-time0)/timeinc
973 
 
 990 
 
Line 1007... Line 1024...
1007 
               print*,'WARNING: y0>90 ',y0,' => setting to 180-y0 ',180.-y0
 1024 
               print*,'WARNING: y0>90 ',y0,' => setting to 180-y0 ',180.-y0
1008 
               y0=180.-y0
 1025 
               y0=180.-y0
1009 
               x0=x0+per/2.
 1026 
               x0=x0+per/2.
1010 
            endif
 1027 
            endif
1011 
            if ((hem.eq.1).and.(y0.lt.-90.)) then
 1028 
            if ((hem.eq.1).and.(y0.lt.-90.)) then
1012 
               print*,'WARNING: y0<90 ',y0,' => setting to -180-y0 ',-180.-y0
 1029 
               print*,'WARNING: y0<-90 ',y0,' => setting to -180-y0 ',-180.-y0
1013 
               y0=-180.-y0
 1030 
               y0=-180.-y0
1014 
               x0=x0+per/2.
 1031 
               x0=x0+per/2.
1015 
            endif
 1032 
            endif
1016 
 
 1033 
 
1017 
            ! Get the index where to interpolate (x0,y0,p0)
 1034 
            ! Get the index where to interpolate (x0,y0,p0)
Line 1022... Line 1039...
1022 
               xind = mdv
 1039 
               xind = mdv
1023 
               yind = mdv
 1040 
               yind = mdv
1024 
               pind = mdv
 1041 
               pind = mdv
1025 
            endif
 1042 
            endif
1026 
 
 1043 
 
1027 
            ! Check if point is within grid (keep indices if ok)
 1044 
           ! Check if point is within grid (keep indices if ok)
1028 
            if ( (xind.ge.1.).and.(xind.le.real(nx)).and. &
 1045 
            if ( (xind.ge.1.).and.(xind.le.real(nx)).and. &
1029 
                 (yind.ge.1.).and.(yind.le.real(ny)).and. &
 1046 
                 (yind.ge.1.).and.(yind.le.real(ny)).and. &
1030 
                 (pind.ge.1.).and.(pind.le.real(nz)) ) then
 1047 
                 (pind.ge.1.).and.(pind.le.real(nz)) ) then
1031 
		 xind = xind
 1048 
		 xind = xind
1032 
		 yind = yind
 1049 
		 yind = yind
1033 
		 pind = pind
 1050 
		 pind = pind
1034 
 
 1051 
 
1035 
            ! Check if pressure is outside, but rest okay => adjust to lowest or highest level
 1052 
            ! Check if pressure is outside, but rest okay => adjust to lowest or highest level
1036 
            elseif ( (xind.ge.1.).and.(xind.le.real(nx)).and. (yind.ge.1.).and.(yind.le.real(ny)) ) then ! only vertical problem
 1053 
            elseif ( (xind.ge.1.).and.(xind.le.real(nx)).and. (yind.ge.1.).and.(yind.le.real(ny)) ) then ! only vertical problem
- 
 
 1054 
 
1037 
               if ( pind.gt.nz ) then ! pressure too low, index too high
 1055 
               if ( pind.gt.nz ) then ! pressure too low, index too high
- 
 
 1056 
                 err_c1 = err_c1 + 1
- 
 
 1057 
                 if ( err_c1.lt.10 ) then
1038 
                 print('Af5.3A'),' WARNING: pressure too low (pind = ',pind,') => adjusted to highest level (pind=nz.)'
 1058 
                    write(*,'(Af5.3A)') ' WARNING: pressure too low (pind = ',pind,') => adjusted to highest level (pind=nz.)'
1039 
                 print*,'(x0,y0,p0)=',x0,y0,p0
 1059 
                    print*,'(x0,y0,p0)=',x0,y0,p0
- 
 
 1060 
                    pind = real(nz)
- 
 
 1061 
                 elseif ( err_c1.eq.10 ) then
- 
 
 1062 
                    print*,' WARNING: more pressures too low -> adjusted to highest level '
- 
 
 1063 
                    pind = real(nz)
- 
 
 1064 
                 else
1040 
                 pind = real(nz)
 1065 
                    pind = real(nz)
- 
 
 1066 
                 endif
- 
 
 1067
1041 
               elseif (pind.lt.1.) then ! pressure too high, index too low
 1068 
               elseif (pind.lt.1.) then ! pressure too high, index too low
- 
 
 1069 
                 err_c2 = err_c2 + 1
- 
 
 1070 
                 if ( err_c2.lt.10 ) then
1042 
                 print('Af5.3A'),' WARNING: pressure too high (pind = ',pind,') => adjusted to lowest level, (pind=1.)'
 1071 
                    write(*,'(Af5.3A)') ' WARNING: pressure too high (pind = ',pind,') => adjusted to lowest level, (pind=1.)'
1043 
                 print*,'(x0,y0,p0)=',x0,y0,p0
 1072 
                    print*,'(x0,y0,p0)=',x0,y0,p0
- 
 
 1073 
                    pind = 1.
- 
 
 1074 
                 elseif ( err_c2.eq.10 ) then
- 
 
 1075 
                    print*,' WARNING: more pressures too high -> adjusted to lowest level '
- 
 
 1076 
                    pind = 1.
- 
 
 1077 
                 else
1044 
                 pind = 1.
 1078 
                    pind = 1.
- 
 
 1079 
                 endif
- 
 
 1080 
 
1045 
              endif
 1081 
              endif
1046 
 
 1082 
 
- 
 
 1083 
            ! Grid point is outside!
1047 
            else
 1084 
            else
- 
 
 1085
- 
 
 1086 
               err_c3 = err_c3 + 1
- 
 
 1087 
               if ( err_c3.lt.10 ) then
1048 
               print*,'ERROR: point is outside grid (horizontally)'
 1088 
                  print*,'ERROR: point is outside grid (horizontally)'
- 
 
 1089 
                  print*,'   Trajectory # ',k
- 
 
 1090 
                  print*,'   Position     ',x0,y0,p0
1049 
               print*,'Grid (nx,ny): ',nx,ny,' and (xind,yind): ',xind,yind
 1091 
                  print*,'  (xind,yind):  ',xind,yind
- 
 
 1092 
                  xind          = mdv
- 
 
 1093 
                  yind          = mdv
- 
 
 1094 
                  pind          = mdv
- 
 
 1095 
                  traint(k,j,2) = mdv
- 
 
 1096 
                  traint(k,j,3) = mdv
- 
 
 1097 
                  traint(k,j,4) = mdv
- 
 
 1098 
               elseif ( err_c3.eq.10 ) then
- 
 
 1099 
                  print*,'ERROR: more points outside grid (horizontally)'
- 
 
 1100 
                  xind          = mdv
- 
 
 1101 
                  yind          = mdv
- 
 
 1102 
                  pind          = mdv
- 
 
 1103 
                  traint(k,j,2) = mdv
- 
 
 1104 
                  traint(k,j,3) = mdv
- 
 
 1105 
                  traint(k,j,4) = mdv
1050 
               stop
 1106 
               else
- 
 
 1107 
                  xind          = mdv
- 
 
 1108 
                  yind          = mdv
- 
 
 1109 
                  pind          = mdv
- 
 
 1110 
                  traint(k,j,2) = mdv
- 
 
 1111 
                  traint(k,j,3) = mdv
- 
 
 1112 
                  traint(k,j,4) = mdv
- 
 
 1113 
               endif
- 
 
 1114 
 
1051 
            endif
 1115 
            endif
1052 
 
 1116 
 
1053 
            ! ------------------------ NEAREST mode -------------------------------
1117 
            ! ------------------------ NEAREST mode -------------------------------
1054 
 
 - 
 
- 
 
 1118 
            ! Interpolate to nearest grid point
1055 
            if ( intmode.eq.'nearest') then
 1119 
            if ( intmode.eq.'nearest') then
1056 
 
 1120 
 
1057 
               xind = real( nint(xind) )
 1121 
               xind = real( nint(xind) )
1058 
               yind = real( nint(yind) )
 1122 
               yind = real( nint(yind) )
1059 
               pind = real( nint(pind) )
 1123 
               pind = real( nint(pind) )
Line 1072... Line 1136...
1072 
               endif
 1136 
               endif
1073 
 
 1137 
 
1074 
               ! interpolate
 1138 
               ! interpolate
1075 
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
 1139 
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
1076 
 
 1140 
 
1077 
            ! ------------------------ CLUSTERING mode ------------------------------- Bojan
 1141 
            ! ------------------------ end NEAREST mode -------------------------------
1078 
 
 1142 
 
- 
 
 1143 
            ! ------------------------ CLUSTERING mode ------------------------------- Bojan
1079 
            elseif (intmode.eq.'clustering') then
 1144 
            elseif (intmode.eq.'clustering') then
1080 
 
 - 
 
1081 
            ! Clustering mode only works with LABEL
 - 
 
1082 
            if (varname.ne.'LABEL' ) then
1145 
               if (varname.ne.'LABEL' ) then
1083 
               print*,'ERROR (clustering): varname is not LABEL'
 1146 
                  print*,'ERROR (clustering): varname is not LABEL'
1084 
               stop
 1147 
                  stop
1085 
            endif
 1148 
               endif
1086 
 
 1149 
 
1087 
            ! Get indices of box around the point
 1150 
            ! Get indices of box around the point
1088 
            xindb = floor(xind)
 1151 
            xindb=floor(xind)
1089 
            xindf = ceiling(xind)
 1152 
            xindf=ceiling(xind)
1090 
            yindb = floor(yind)
 1153 
            yindb=floor(yind)
1091 
            yindf = ceiling(yind)
 1154 
            yindf=ceiling(yind)
1092 
            pindb = floor(pind)
 1155 
            pindb=floor(pind)
1093 
            pindf = ceiling(pind)
 1156 
            pindf=ceiling(pind)
1094 
 
 1157 
 
1095 
            ! Make sure all points are within grid
 1158 
            ! Make sure all points are within grid
1096 
            if ( xindb.lt.1  ) xindb = 1
 1159 
            if ( xindb.lt.1 ) xindb = 1
1097 
            if ( xindf.lt.1  ) xindf = 1
 1160 
            if ( xindf.lt.1 ) xindf = 1
1098 
            if ( xindb.gt.nx ) xindb = nx
 1161 
            if ( xindb.gt.nx ) xindb = nx
1099 
            if ( xindf.gt.nx ) xindf = nx
 1162 
            if ( xindf.gt.nx ) xindf = nx
1100 
            if ( yindb.lt.1  ) yindb = 1
 1163 
            if ( yindb.lt.1 ) yindb = 1
1101 
            if ( yindf.lt.1  ) yindf = 1
 1164 
            if ( yindf.lt.1 ) yindf = 1
1102 
            if ( yindb.gt.ny ) yindb = ny
 1165 
            if ( yindb.gt.ny ) yindb = ny
1103 
            if ( yindf.gt.ny ) yindf = ny
 1166 
            if ( yindf.gt.ny ) yindf = ny
1104 
            if ( pindb.lt.1  ) pindb = 1
 1167 
            if ( pindb.lt.1 ) pindb = 1
1105 
            if ( pindf.lt.1  ) pindf = 1
 1168 
            if ( pindf.lt.1 ) pindf = 1
1106 
            if ( pindb.gt.nz ) pindb = nz
 1169 
            if ( pindb.gt.nz ) pindb = nz
1107 
            if ( pindf.gt.nz ) pindf = nz
 1170 
            if ( pindf.gt.nz ) pindf = nz
1108 
 
 1171 
 
1109 
            ! Shift one point if they are equal
 1172 
            ! Shift one point if they are equal
1110 
            if ( xindb.eq.xindf ) then
 1173 
            if ( xindb.eq.xindf ) then
Line 1126... Line 1189...
1126 
                  pindb=nz-1
 1189 
                  pindb=nz-1
1127 
               else
 1190 
               else
1128 
                  pindf=pindb+1
 1191 
                  pindf=pindb+1
1129 
               endif
 1192 
               endif
1130 
            endif
 1193 
            endif
1131 
 
 - 
 
1132 
            ! Give warnings and stop if problems occur
 1194 
            ! Give warnings and stop if problems occur
1133 
            if ( xindb.eq.xindf ) then
 1195 
            if ( xindb.eq.xindf ) then
1134 
               print*,'ERROR (clustering): xindb=xindf'
 1196 
               print*,'ERROR (clustering): xindb=xindf'
1135 
               print*,xind,xindb,xindf
 1197 
               print*,xind,xindb,xindf
1136 
               stop
 1198 
               stop
Line 1188... Line 1250...
1188 
 
 1250 
 
1189 
            ! Set to -9 to detect if no label was assigned in the end
 1251 
            ! Set to -9 to detect if no label was assigned in the end
1190 
            f0=-9
 1252 
            f0=-9
1191 
 
 1253 
 
1192 
            ! Stratosphere (PV)
 1254 
            ! Stratosphere (PV)
1193 
            if ( abs(traint(k,j,pvpos)) .ge. 2.0 ) then
 1255 
            if ( abs(traint(k,j,pvpos)) .ge. tropo_pv ) then
1194 
               if ( (lblcount(2).ge.lblcount(3)).and. (lblcount(2).ge.lblcount(5)) ) then
 1256 
               if ( (lblcount(2).ge.lblcount(3)).and. (lblcount(2).ge.lblcount(5)) ) then
1195 
                  f0=2
 1257 
                  f0=2
1196 
               elseif ( lblcount(3).ge.lblcount(5) ) then
 1258 
               elseif ( lblcount(3).ge.lblcount(5) ) then
1197 
                  f0=3
 1259 
                  f0=3
1198 
               elseif ( lblcount(5).gt.lblcount(3) ) then
 1260 
               elseif ( lblcount(5).gt.lblcount(3) ) then
1199 
                  f0=5
 1261 
                  f0=5
1200 
               endif
 1262 
               endif
1201 
            endif
 1263 
            endif
1202 
 
 1264 
 
1203 
            ! Troposphere (PV)
 1265 
            ! Troposphere (PV)
1204 
            if ( abs(traint(k,j,pvpos)) .lt. 2.0 ) then
 1266 
            if ( abs(traint(k,j,pvpos)) .lt. tropo_pv ) then
1205 
               if ( lblcount(1).ge.lblcount(4) ) then
 1267 
               if ( lblcount(1).ge.lblcount(4) ) then
1206 
               f0=1
 1268 
               f0=1
1207 
               elseif ( lblcount(4).gt.lblcount(1) ) then
 1269 
               elseif ( lblcount(4).gt.lblcount(1) ) then
1208 
               f0=4
 1270 
               f0=4
1209 
               endif
 1271 
               endif
1210 
            endif
 1272 
            endif
1211 
 
 1273 
 
1212 
            ! Stratosphere (TH)
 1274 
            ! Stratosphere (TH)
1213 
            if ( traint(k,j,thpos) .ge. 380.0 ) then
 1275 
            if ( traint(k,j,thpos) .ge. tropo_th ) then
1214 
               f0=2
 1276 
               f0=2
1215 
            endif
 1277 
            endif
1216 
 
 1278 
 
1217 
            if (f0.eq.-9) then
 1279 
            if (f0.eq.-9) then
1218 
               print*,'ERROR (Clustering): No label assigned!'
 1280 
               print*,'ERROR (Clustering): No label assigned!'
1219 
               stop
 1281 
               stop
1220 
            endif
 1282 
            endif
- 
 
 1283 
            ! ------------------------ end CLUSTERING mode -------------------------------
1221 
 
 1284 
 
1222 
            ! ------------------------ CIRCLE modes ------------------------------- Bojan
 1285 
            ! ------------------------ CIRCLE modes ------------------------------- Bojan
1223 
 
 - 
 
- 
 
 1286 
            ! elseif (not clustering but one of the possible circle modes)
1224 
            elseif ( (intmode.eq.'circle_avg') .or. (intmode.eq.'circle_min') .or. (intmode.eq.'circle_max') ) then
 1287 
            elseif ( (intmode.eq.'circle_avg') .or. (intmode.eq.'circle_min') .or. (intmode.eq.'circle_max') ) then
1225 
 
 1288 
 
1226 
            ! reset arrays for this point
 1289 
            ! reset arrays for this point
1227 
            connect    = 0
 1290 
            connect=0
1228 
            stackx     = 0
 1291 
            stackx=0
1229 
            stacky     = 0
 1292 
            stacky=0
1230 
            circlelon  = 0
 1293 
            circlelon=0
1231 
            circlelat  = 0
 1294 
            circlelat=0
1232 
            circlef    = 0
 1295 
            circlef=0
1233 
            circlearea = 0
 1296 
            circlearea=0
1234 
 
 1297 
 
1235 
            ! Get indices of one coarse grid point within search radius (nint=round to next integer)
 1298 
            ! Get indices of one coarse grid point within search radius (nint=round to next integer)
1236 
            if ( sdis(x0,y0,longrid(nint(xind)),latgrid(nint(yind))) .gt. radius) then
 1299 
            if ( sdis(x0,y0,longrid(nint(xind)),latgrid(nint(yind))) .gt. radius) then
1237 
               print*,'ERROR (circle): Search radius is too small... (1). r =',radius
 1300 
               print*,'ERROR (circle): Search radius is too small... (1). r =',radius
1238 
               print*,'Distance to nint grid point (minimum search radius)=',sdis(x0,y0,longrid(nint(xind)),latgrid(nint(yind)))
 1301 
               print*,'Distance to nint grid point (minimum search radius)=',sdis(x0,y0,longrid(nint(xind)),latgrid(nint(yind)))
1239 
               stop
 1302 
               stop
1240 
            endif
 1303 
            endif
1241
1304
1242 
            ! Initialize stack with nint(xind),nint(yind)
 1305 
            ! Initialize stack with nint(xind),nint(yind)
1243 
            kst       = 0            ! counts the number of points in circle
 1306 
            kst=0 ! counts the number of points in circle
1244 
            stackx(1) = nint(xind)
 1307 
            stackx(1)=nint(xind)
1245 
            stacky(1) = nint(yind)
 1308 
            stacky(1)=nint(yind)
1246 
            sp        = 1            ! stack counter
 1309 
            sp=1 ! stack counter
1247 
 
 - 
 
1248 
            ! Repeat as long as the stack is not empty
 - 
 
1249 
            do while (sp.ne.0)
 1310 
            do while (sp.ne.0)
1250 
 
 1311 
 
1251 
            ! Get an element from stack
 1312 
            ! Get an element from stack
1252 
             ist=stackx(sp)
 1313 
             ist=stackx(sp)
1253 
             jst=stacky(sp)
 1314 
             jst=stacky(sp)
Line 1313... Line 1374...
1313 
 
 1374 
 
1314 
             endif ! endif radius is smaller => end of updating stack
 1375 
             endif ! endif radius is smaller => end of updating stack
1315 
 
 1376 
 
1316 
            end do ! end working on stack
1377 
            end do ! end working on stack
1317 
 
 1378 
 
1318 
	    ! Choose output depending on intmode
- 
 
1319 
            if (kst.ge.1) then
 1379 
            if (kst.ge.1) then
1320 
 
 - 
 
1321 
               ! calculate area-weighted average of f in circle
1380 
               ! Choose output depending on intmode
1322 
               if ( intmode .eq. 'circle_avg' ) then
 1381 
               if ( intmode .eq. 'circle_avg' ) then
- 
 
 1382 
                  ! calculate area-weighted average of f in circle
1323 
                  circlesum=0.
 1383 
                  circlesum=0.
1324 
                  do l=1,kst
 1384 
                  do l=1,kst
1325 
                     circlesum=circlesum+circlef(l)*circlearea(l)
 1385 
                     circlesum=circlesum+circlef(l)*circlearea(l)
1326 
                  enddo
 1386 
                  enddo
1327 
                  circleavg=circlesum/sum(circlearea(1:kst))
 1387 
                  circleavg=circlesum/sum(circlearea(1:kst))
- 
 
 1388 
                  !print*,'area-weighted average of f in circle=',circleavg
1328 
                  f0=circleavg
 1389 
                  f0=circleavg
1329 
 
 - 
 
1330 
               ! calculate minimum in circle
 - 
 
1331 
               elseif ( intmode .eq. 'circle_min' ) then
 1390 
               elseif ( intmode .eq. 'circle_min' ) then
- 
 
 1391 
                  ! calculate minimum in circle
1332 
                  circlemin=circlef(1)
 1392 
                  circlemin=circlef(1)
1333 
                  do l=1,kst
 1393 
                  do l=1,kst
1334 
                     if (circlef(l) .lt. circlemin) then
 1394 
                     if (circlef(l) .lt. circlemin) then
1335 
                        circlemin=circlef(l)
 1395 
                        circlemin=circlef(l)
1336 
                     endif
 1396 
                     endif
1337 
                  enddo
 1397 
                  enddo
- 
 
 1398 
                  !print*,'minimum of f in circle=',circlemin
1338 
                  f0=circlemin
 1399 
                  f0=circlemin
1339 
 
 - 
 
1340 
               ! calculate maximum in circle
 - 
 
1341 
               elseif ( intmode .eq. 'circle_max' ) then
1400 
               elseif ( intmode .eq. 'circle_max' ) then
- 
 
 1401 
                  ! calculate maximum in circle
1342 
                  circlemax=circlef(1)
 1402 
                  circlemax=circlef(1)
1343 
                  do l=1,kst
 1403 
                  do l=1,kst
1344 
                     if (circlef(l) .gt. circlemax) then
 1404 
                     if (circlef(l) .gt. circlemax) then
1345 
                        circlemax=circlef(l)
 1405 
                        circlemax=circlef(l)
1346 
                     endif
 1406 
                     endif
1347 
                  enddo
 1407 
                  enddo
- 
 
 1408 
                  !print*,'maximum of f in circle=',circlemax
1348 
                  f0=circlemax
 1409 
                  f0=circlemax
1349 
 
 - 
 
1350 
               else
1410 
               else
1351 
                  print*,'ERROR (circle): intmode not valid!'
 1411 
                  print*,'ERROR (circle): intmode not valid!'
1352 
                  stop
 1412 
                  stop
1353 
               endif
 1413 
               endif
1354 
 
 - 
 
1355 
            else
 1414 
            else
1356 
               print*,'ERROR (circle): Search radius is too small... (2). r =',radius
 1415 
               print*,'ERROR (circle): Search radius is too small... (2). r =',radius
1357 
               stop
 1416 
               stop
1358 
            endif
 1417 
            endif
1359 
 
 1418 
 
- 
 
 1419 
            ! ------------------------ end CIRCLE modes -------------------------------
1360 
 
 1420 
 
1361 
            ! ------------------------ NORMAL mode -------------------------------
 1421 
            ! ------------------------ NORMAL mode -------------------------------
1362 
            else
1422 
            else ! not clustering nor circle: NORMAL mode
1363 
 
 1423 
 
- 
 
 1424 
            ! Check if point is within grid
- 
 
 1425 
!            if ( (xind.ge.1.).and.(xind.le.real(nx)).and. &
- 
 
 1426 
!                 (yind.ge.1.).and.(yind.le.real(ny)).and. &
- 
 
 1427 
!                 (pind.ge.1.).and.(pind.le.real(nz)) ) then
- 
 
 1428 
!
- 
 
 1429 
            ! Do the interpolation: everthing is ok
1364 
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
 1430 
               f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
1365 
 
 1431 
 
- 
 
 1432 
!            ! Check if pressure is outside, but rest okay: adjust to lowest or highest level
- 
 
 1433 
!            elseif ( (xind.ge.1.).and.(xind.le.real(nx)).and. (yind.ge.1.).and.(yind.le.real(ny)) ) then ! only vertical problem
- 
 
 1434 
!               if ( pind.gt.nz ) then ! pressure too low, index too high
- 
 
 1435 
!                 pind = real(nz)
- 
 
 1436 
!                 print*,' Warning: pressure too low -> adjusted to highest level, pind=nz.'
- 
 
 1437 
!                 print*,'(x0,y0,p0)=',x0,y0,p0
- 
 
 1438 
!               elseif (pind.lt.1.) then ! pressure too high, index too low
- 
 
 1439 
!                 pind = 1.
- 
 
 1440 
!                 print*,' Warning: pressure too high -> adjusted to lowest level, pind=1.'
- 
 
 1441 
!                 print*,'(x0,y0,p0)=',x0,y0,p0
- 
 
 1442 
!              endif
- 
 
 1443 
!              f0 = int_index4(f3t0,f3t1,nx,ny,nz,xind,yind,pind,reltpos0,mdv)
- 
 
 1444 
 
- 
 
 1445 
!            ! Less than 10 outside
- 
 
 1446 
!            elseif (noutside.lt.10) then
- 
 
 1447 
!               print*,' ',trim(tvar(i)),' @ ',x0,y0,p0,'outside'
- 
 
 1448 
!               f0       = mdv
- 
 
 1449 
!               noutside = noutside + 1
- 
 
 1450 
!
- 
 
 1451 
!            ! More than 10 outside
- 
 
 1452 
!            elseif (noutside.eq.10) then
- 
 
 1453 
!               print*,' ...more than 10 outside...'
- 
 
 1454 
!               f0       = mdv
- 
 
 1455 
!               noutside = noutside + 1
- 
 
 1456 
 
- 
 
 1457 
!            ! Else (not everything okay and also not 'tolerated cases') set to missing data
- 
 
 1458 
!            else
- 
 
 1459 
!               f0       = mdv
1366 
            endif
 1460 
!            endif
- 
 
 1461 
 
- 
 
 1462 
            ! ------------------------ end NORMAL mode -------------------------------
- 
 
 1463 
            endif ! end if nearest case
- 
 
 1464 
 
- 
 
 1465 
           ! Exit for loop over all trajectories and times -save interpolated value
- 
 
 1466 
 109        continue
1367 
 
 1467 
 
1368 
            ! Save the new field
 1468 
            ! Save the new field
1369 
            if ( abs(f0-mdv).gt.eps) then
 1469 
            if ( abs(f0-mdv).gt.eps) then
1370 
               traint(k,j,ncol+i) = f0
 1470 
               traint(k,j,ncol+i) = f0
1371 
            else
 1471 
            else
Line 1889... Line 1989...
1889 
  print*
 1989 
  print*
1890 
  print*,'              *** END OF PROGRAM TRACE ***'
 1990 
  print*,'              *** END OF PROGRAM TRACE ***'
1891 
  print*,'========================================================='
 1991 
  print*,'========================================================='
1892 
 
 1992 
 
1893 
 
 1993 
 
1894 
end program trace_label_circle
 1994 
end program trace
1895 
 
 1995 
 
1896 
 
 1996 
 
1897 
 
 1997 
 
1898 
! ******************************************************************
 1998 
! ******************************************************************
1899 
! * SUBROUTINE SECTION                                             *
 1999 
! * SUBROUTINE SECTION                                             *