WebSVN – lagranto.icon – Diff – /trunk/select/special_stt.f


 
      SUBROUTINE special_stt (flag,cmd,tra,ntim,ncol,
     >                    vars,times,param,nparam,outfile)
 
c     ***************************************************************************
c     *                                                                         *
c     * OUTPUT:  flag           -> 1 if trajectory is selected, 0 if not        *
c     *                                                                         *
c     * INPUT:   cmd            <- command string                               *
c     *          tra(ntim,ncol) <- single trajectory: indices time,column       *
c     *          ntim           <- number of times                              *
c     *          ncol           <- number of columns (including time,lon,lat,p) *
c     *          vars(ncol)     <- names of columns                             *
c     *          times(ntim)    <- List of times
c     *          param(nparam)  <- parameter values                             *
c     *          nparam         <- number of parameters                         *
c     *                                                                         *
c     ***************************************************************************
 
      implicit none
      
c     ---------------------------------------------------------------------------
c     Declaration of subroutine parameters
c     ---------------------------------------------------------------------------
 
      integer       flag           ! Boolean flag whether trajectory is selected
      character*80  cmd            ! Command string
      integer       ntim,ncol      ! Dimension of single trajectory
      real          tra(ntim,ncol) ! Single trajectory
      character*80  vars(ncol)     ! Name of columns
      real          times(ntim)    ! List of times
      integer       nparam         ! # parameters
      real          param(1000)    ! List of parameters
      character*80  outfile        ! NAme of output file
 
      integer       i,n		   ! Loop indices
 
c     ---------------------------------------------------------------------------
c     Local variables for SPECIAL:WCB
c     ---------------------------------------------------------------------------
 
      integer       ip,i0,i1
 
c     ---------------------------------------------------------------------------
c     Local variables for SPECIAL:STT
c     ---------------------------------------------------------------------------
 
      integer       debugflag
      integer       tb,ta					!tb and ta are the minimum residence times before and after the exchange
      integer	    stat, ex, outputflag			!Variables for I/O files
      integer	    labelflag					!Flag to handle label field
      integer	    pvflag, thflag, zflag, pflag		!Flags to signal presence of various fields
      integer	    pvcol,thcol,labelcol,zcol,pcol		!Positions of PV, TH, LABEL, z and p columns
      real	    mean1, mean2, sd1, sd2			!Mean and standard deviations of PV and TH
      real	    deltat					!Trajectory timestep
      real	    extime, exlon, exlat, exth, exheight, expres!Position at exchange time	
      character*80  namefile            			!Name of output file
      real	    a,b						!Normalized tropopause distance for PV and TH
      real	    thtrop, pvtrop				!Tropopause level in K and pvu
      integer	    count1, count2
      integer	    sttcount, p				!Ancillary variables in STT do while loop
      real,allocatable, dimension (:)   :: tropdist	! Vector for normalized distance from tropopause
      real,allocatable, dimension (:)   :: exflag	! Vector for normalized distance from tropopause
      integer, dimension (ntim)		:: pos		! Vector to store exchange indices
 
 
 
c     --------------------------------------------------------------------------  %)
c     SPECIAL:STT:pvtrop,thtrop,tmin_b,tmin_f,outputflag 		                  	  %)
                                     
c     --------------------------------------------------------------------------- %)
 
      if ( cmd.eq.'STT' ) then
 
c        Reset the flag for selection
         flag = 0
c	 Reset the flag for debugging
	 debugflag = 0
c	 Set other flags
         labelflag = 0
	 pvflag = 0
	 thflag = 0
C	 Set data at the exchange
         extime = 0
	 exlon = 0 
	 exlat = 0
	 exheight = 0
	 expres = 0
 
c	 PV and TH threshold read from command line
	 pvtrop = param(1)
	 thtrop = param(2)
 
 
c	Determine number of time steps in the trajectory 
c	Handles backward OR forward trajectories
	
	if ( (ntim.eq.0).or.(ntim.eq.1) ) then
          print*,' ERROR: zero length trajectory in SPECIAL:STT... Stop'	!If there are problems, stop computation.
            stop
         endif
 
	deltat = (times(ntim)-times(1))/(ntim-1)				!By default array indexing starts in Fortran at 1
 
 
C 	Now some checks for essential fields : label, potential vorticity, potential temperature
 
c	Check for the presence of the LABEL field
	labelcol = 1
	do i = 1,ncol
	    if((vars(i) .eq. 'LABEL').or.(vars(i) .eq. 'label')) then
		labelflag = 1
		labelcol = i
	    if(debugflag .eq. 1) then
		print*, 'Found label'
	    endif
	    endif
	enddo
 
 
c	Check for the presence of the variables PV and TH
	
	pvcol = 1	
	thcol = 1
	zcol = 1
	pcol = 1			!In case th or pv are not available, tra(timestep,thcol) is always defined
 
	
	do i = 1,ncol
	    if((vars(i) .eq. 'PV') .or. (vars(i) .eq. 'POT_VORTI')) then
		pvflag = 1
		pvcol = i
	    endif
	enddo
 
 
	do i = 1,ncol
	    if((vars(i) .eq. 'TH') .or. (vars(i) .eq. 'POT_T')) then
		thflag = 1
		thcol = i
	    endif
	enddo
 
	do i = 1,ncol
	    if((vars(i) .eq. 'z') .or. (vars(i) .eq. 'Z')) then
		zcol = i
		zflag = 1
	    endif
	enddo
 
	do i = 1,ncol
	    if((vars(i) .eq. 'p') .or. (vars(i) .eq. 'P')) then
		pcol = i
		pflag = 1
	    endif
	enddo
 
c	No PV and THETA available: not possible determining tropopause. Program stop. 
 
	if ( (thflag .eq. 0).and.(pvflag .eq. 0) ) then
            print*,' ERROR: no PV and no TH in SPECIAL:STT..'
	    print*,' ERROR: not possible determining tropopause'	
            stop
         endif
 
c	Check for missing data in the trajectories
 
  	do n=1,ntim
     		if ((tra(n,pvcol).lt.-500.).and.(pvflag.eq.1.)) then
	          if(debugflag .eq. 1) then
            		print*,'WARNING: PV missing data',tra(n,pvcol)
			print*,'Trajectory skipped'
		  endif
			goto 600
    		 endif 
		if ((tra(n,thcol) .lt.-500.).and.(thflag .eq. 1.)) then
		  if(debugflag .eq. 1) then
            		print*,'WARNING: TH missing data',tra(n,thcol)
			print*,'Trajectory skipped'
		  endif
			goto 600
    		 endif 
  	enddo
 
c	Computing mean of PV and TH along a trajectory
 
	mean1=0.
  	count1=0
  	mean2=0.
  	count2=0
  	do n=1,ntim						
     	    if ((tra(n,pvcol).gt.-100.).and.(pvflag .eq. 1)) then
        	count1=count1+1
        	mean1=mean1+tra(n,pvcol)
	    else
	    if(debugflag .eq. 1) then
		print*,'WARNING: PV is missing or non-physical, mean1 = 0'
	    endif
		mean1 = 0
	    endif
     	    if ((tra(n,thcol) .gt. 1E-6) .and. (thflag .eq. 1)) then
        	count2=count2+1
        	mean2=mean2+tra(n,thcol)
	    else
	    if(debugflag .eq. 1) then
		print*,'WARNING: TH is missing or non-physical, mean1 = 0'
	    endif
		mean2 = 0
	    endif
  	enddo
 
c	Compute average, otherwise skip problematic trajectories
 
  	if (count1.ne.0) then
     		mean1=mean1/count1
  	endif
 
  	if (count2.ne.0) then
     		mean2=mean2/count2
  	endif
 
        if ( (count1.eq.0) .and. (count2.eq.0) )then
	    if(debugflag .eq. 1) then
		print*, 'Trajectory with only missing data skipped' 
	    endif       	
		goto 600
        endif
 
c	Computing standard deviation of PV and TH along a trajectory
 
	sd1=0.
  	count1=0.
  	sd2=0.
  	count2=0.
  	do n=1,ntim						
     	    if (tra(n,pvcol).gt.-100.) then
        	count1=count1+1.
        	sd1=sd1+(tra(n,pvcol)-mean1)**2
	    endif
     	    if (tra(n,thcol).gt. 0.) then
        	count2=count2+1.
        	sd2=sd2+(tra(n,thcol)-mean2)**2
	    endif
  	enddo
	sd1=sqrt(sd1/(count1-1.))
  	sd2=sqrt(sd2/(count2-1.))
 
c     Allocate distance array
 
	allocate(tropdist(1:ntim),stat=stat)
	if (stat.ne.0) then 
	    print*,'ERROR: allocating array tropdist 
     &not possible in SPECIAL:STT.. '
	    stop
	endif
 
c	If LABEL is not available, go to section without usage of label
 
	if(labelflag .ne. 1) goto 300
 
 
C-------------------------------------------------------------------------------------------------------------------------
c	BEGINNING OF THE SECTION WITH LABEL
C-------------------------------------------------------------------------------------------------------------------------
 
	
c	Check for label field
	do i = 1,ntim
	    if((tra(i,labelcol) .ne.1) .and. (tra(i,labelcol) .ne.2)) then
	      if(debugflag .eq. 1) then
	    	print*,'WARNING: unexpected label field in SPECIAL:STT.. '
	      endif
		print*,'WARNING: unused label field in SPECIAL:STT'
	    	goto 300
	    endif
	end do
 
 
c     Distance from the '2 PVU, 380 K' tropopause (-:trop/+:strat) : pvtrop, thtrop if PV and TH and LABEL are available
	
	if ( (pvflag .ne. 0) .and. (thflag .ne. 0) .and. 
     &(labelflag .eq. 1)) then 
	    if(debugflag .eq. 1) then
	      print*,'Tropopause detected with PV, TH and LABEL'
            endif
	  do n=1,ntim
     	    a=abs(pvtrop-abs(tra(n,pvcol)))/sd1 	! a = tropdist(PV)
     	    b=abs(thtrop-tra(n,thcol))/sd2 		! b = tropdist(TH)
     	    if (tra(n,labelcol) .eq. 1) then ! Troposphere (LABEL)  
        	tropdist(n)=-min(a,b)
     	    else
         	tropdist(n)=+min(a,b)
     	    endif
  	  enddo
	endif
 
c	If PV is not available, use only TH
	
	if  (pvflag .eq. 0) then 
	  if(debugflag .eq. 1) then
	    print*,'Tropopause detected with TH only'
	  endif
	  do n=1,ntim
     	    b=abs(thtrop-tra(n,thcol))/sd2 		! b = tropdist(TH)
     	    if(tra(n,labelcol) .eq. 1) then 		! Troposphere (TH)
        	tropdist(n)= -b
	    else
	        tropdist(n)= b				! Stratosphere (TH)
     	    endif
  	  enddo
	endif
 
c	If TH is not available, use only PV
	if  (thflag .eq. 0) then
	  if(debugflag .eq. 1) then
	    print*,'Tropopause detected with PV only'	
	  endif
	  do n=1,ntim
     	    a=abs(pvtrop-tra(n,pvcol))/sd2 		! a = tropdist(PV)
     	    if( (tra(n,pvcol).lt.pvtrop) .and. 
     &(tra(n,labelcol) .eq. 1)) then 		! Troposphere (PV)
        	tropdist(n)= -a
     	    else
	        tropdist(n)= a				! Stratosphere (PV)
     	    endif
  	  enddo
	endif
 
C	Detecting STT with LABEL field
c	If PV is lower than pvtrop despite the label is 2, then in next step this potential STT is not considered.
	
	sttcount = 0
	
	do i=2,ntim
	    if((tropdist(i)*tropdist(i-1) .lt. 0) .and. 
     &(tropdist(i-1) .gt. 0)) then
C		Adding label check for STT
		if((tra(i,labelcol).ne.1).and.(tra(i-1,labelcol).ne.2))then
		  if(debugflag .eq. 1) then
	            print*, 'STT label check failed:' 
		    print*, 'moving to next STT candidate'
	          endif
		    cycle
		else
		    sttcount = sttcount+1
		    pos(sttcount) = i-1
		endif
	    endif
	enddo
	
c	If no STT is detected, move to next trajectory
 
	if(sttcount .lt. 1 ) then
	  if(debugflag .eq. 1) then
	    print*, 'No STT detected: moving to next trajectory'
	  endif
	    goto 600
	endif
	
c	Re-connection with normal steps: label has been used
 
	goto 350
 
C-------------------------------------------------------------------------------------------------------------------------
C	END OF LABEL SECTION
c-------------------------------------------------------------------------------------------------------------------------
 
 
 
c     Distance from the '2 PVU, 380 K' tropopause (-:trop/+:strat) : pvtrop, thtrop if PV and TH are available
300	continue							!Here the code continues if no LABEL has been found
	
	if ( (pvflag .ne. 0) .and. (thflag .ne. 0)) then
	    if(debugflag .eq. 1) then
	      print*,'Tropopause detected with PV and TH'
            endif
	  do n=1,ntim
     	    a=abs(pvtrop-abs(tra(n,pvcol)))/sd1 	! a = tropdist(PV)
     	    b=abs(thtrop-tra(n,thcol))/sd2 		! b = tropdist(TH)
     	    if ((tra(n,thcol).lt.thtrop).and.
     &(abs(tra(n,pvcol)).lt.pvtrop)) then ! Troposphere (TH and PV)	 ---------WHY ABS??---------
        	tropdist(n)=-min(a,b)
     	    else
       	    	if ( (abs(tra(n,pvcol)).ge.pvtrop).and.
     &(tra(n,thcol).ge.thtrop) ) then ! Stratosphere (TH and PV)
           	    tropdist(n)=+min(a,b)
            	endif
 
            	if ( (abs(tra(n,pvcol)).ge.pvtrop).and.
     &(tra(n,thcol).lt.thtrop) ) then ! Stratosphere (PV)
           	    tropdist(n)=a
            	endif
 
            	if ( (abs(tra(n,pvcol)).lt.pvtrop).and.
     &(tra(n,thcol).ge.thtrop) ) then ! Stratosphere (TH)
           	    tropdist(n)=b
            	endif
     	    endif
  	  enddo
	endif
 
c	If PV is not available, use only TH
	
	if  (pvflag .eq. 0) then
	    if(debugflag .eq. 1) then
	      print*,'Tropopause detected with TH only'
            endif
	  do n=1,ntim
     	    b=abs(thtrop-tra(n,thcol))/sd2 		! b = tropdist(TH)
     	    if (tra(n,thcol).lt.thtrop) then 		! Troposphere (TH)
        	tropdist(n)= -b
     	    else
	        tropdist(n)= b				! Stratosphere (TH)
     	    endif
  	  enddo
	endif
 
c	If TH is not available, use only PV
	    if(debugflag .eq. 1) then
	      print*,'Tropopause detected with PV only'
            endif	
	if  (thflag .eq. 0) then
	  do n=1,ntim
     	    a=abs(pvtrop-tra(n,pvcol))/sd2 		! a = tropdist(PV)
     	    if (tra(n,pvcol).lt.pvtrop) then 		! Troposphere (PV)
        	tropdist(n)= -a
     	    else
	        tropdist(n)= a				! Stratosphere (PV)
     	    endif
  	  enddo
	endif
 
 
c	Detect the various STT occurrences and save positions in traj file in array pos
	
	sttcount = 0
	
	do i=2,ntim
	    if((tropdist(i)*tropdist(i-1) .lt. 0) .and.
     & (tropdist(i-1) .gt. 0)) then
		sttcount = sttcount+1
		pos(sttcount) = i-1
	    endif
	enddo
	
c	If no STT are detected, go to next trajectory
 
	if(sttcount .lt. 1 ) then
	  if(debugflag .eq. 1) then
	    print*, 'No STT detected: moving to next trajectory'
	  endif
	    goto 600
	endif
 
 
 
350	continue								!The same if LABEL is present or not.
	
c	In timesteps, minimum residence time before and after the exchange
 
	tb = nint( param(3)/deltat )
	ta = nint( param(4)/deltat )
 
 
c	Now a big while loop analyzing every candidate STT
 
	do while(sttcount>0)
 
	    p = pos(sttcount)	
	
c	    Check whether the trajectory is still in the troposphere after ta timesteps
 
	   if(ta .ne. 0) then
	    	do n=1,ta
	            if((p+n) .le. ntim) then
		    	if (tropdist(p)*tropdist(p+n).gt.0.) goto 500		!Try with new STT candidate
		    else
		    	goto 500						!Out of range. Try with new STT candidate							
	    	    endif
  	    	enddo
	   endif
 
c	    Check whether the trajectory had been in the stratosphere for at least tb timesteps
 
	    if(tb .ne. 0) then
	        do n=1,tb
	            if((p-n) .gt. 0) then
		    	if (tropdist(p)*tropdist(p-n).lt.0.) goto 500		!Try with new STT candidate
		    else
		    	goto 500						!Out of range. Try with new STT candidate
	    	    endif
  	    	enddo
	    endif
 
c	   If the loops have gone fine, the trajectory is considered, flag = 1, values at the exchange are written and the loop stops.
 
	    flag = 1
	    extime = tra(p,1)
	    exlon = tra(p,2)
	    exlat = tra(p,3)
	    if(zflag .eq. 1) then
	    	exheight = tra(p,zcol)
	    endif
	    if(pflag .eq. 1) then
	    	expres = tra(p,pcol)
	    endif
	    if(thflag .eq. 1) then
	    	exth = tra(p,thcol)
	    endif
	    goto 600
	
500	    continue
	    sttcount = sttcount-1 
 
	enddo		! End of do while loop
 
	    goto 600
 
 
 
600   continue
 
      outputflag = param(5)
 
c 	Read the name of the file
 
	    namefile = trim(outfile)//'.stt'
 
 	if( (flag .eq. 1) .and. (outputflag .eq. 1) ) then
 
      if ( nint( param(nparam+1) ).eq.0) then
		open(15,file=namefile,status='replace',	
     &position='append', iostat=stat)
 
c        ---------------------------------------------------------------
       write(15,*) 'time      lon     lat     height  p       th'
       write(15,*) '--------------------------------------------'
       write(15,*) 
            	write(15,'(f7.2,1x,2f8.2,f8.0,9f10.3,5x,i1,4x)') 
     &extime,exlon,exlat,exheight,expres,exth
c        ---------------------------------------------------------------
 
        	close(15)
            param(nparam+1) = 1.
            
      else
 
		open(15,file=namefile,status='old',	
     &position='append', iostat=stat)
 
c        ---------------------------------------------------------------
            	write(15,'(f7.2,1x,2f8.2,f8.0,9f10.3,5x,i1,4x)') 
     &extime,exlon,exlat,exheight,expres,exth
c        ---------------------------------------------------------------
 
        	close(15)
      endif
 
 
 
 
 	endif 
 
 
c	Write output if outputflag = 1
	
c	outputflag = param(5)
 
c	Check for correct inserted parameter
 
c	if(debugflag .eq. 1) then
c	if( (outputflag .ne. 0) .and. (outputflag .ne. 1) ) then
c	  print*, 'WARNING: last parameter must be either 0 or 1'
c	  print*, 'NO data file produced'
c	endif
c	endif
	
c        if( (flag .eq. 1) .and. (outputflag .eq. 1) ) then
 
c	Read the name of the file
c	    namefile = "data_exch_stt"							!Maybe write in the command 'name'
c	Open file (already opened before)
c	
c            inquire(file=namefile, exist=ex)
c	    if (ex) then
c		open(15,file=namefile, status = 'old', 
c     &position='append',iostat=stat)
c	    if(stat .eq. -1) then							!If the file is empty behave as it was non-existent
c	      goto 125
c	    endif
c        Write output for trajectory making STT
c        ---------------------------------------------------------------
c            	write(15,'(f7.2,1x,2f8.2,f6.0,9f10.3,5x,i1,4x)') 
c     &extime,exlon,exlat,exheight,expres,exth
c        ---------------------------------------------------------------
c            	write(15,*) ! empty line
c
c
c	    else
c
c	First time opening
c
c		open(15,file=namefile,status='replace',
c     &position='append', iostat=stat)
c
c        Write title
c        ---------------------------------------------------------------
c 125           	write(15,"(A)") "time	lon	lat	z	p	th"
c          write(15,"(A)") "-------------------------------------------"
    
c        ---------------------------------------------------------------
c            	write(15,*) ! empty line
 
c        Write output for trajectory making STT
c        ---------------------------------------------------------------
c            	write(15,'(f7.2,1x,2f8.2,f6.0,9f10.3,5x,i1,4x)')
c     &extime,exlon,exlat,exheight,expres,exth
c        ---------------------------------------------------------------
c            	write(15,*) ! empty line
c	    endif
 
 
c	    if (stat .ne. 0) then
c		print *,"ERROR: Problem with output file",namefile
c	        stop
c	    endif
c	close(15)
c	endif
 
c	"Closing" the STT criterion
	endif 
 
       end

Rev 3	Rev 5
1		1
2	`SUBROUTINE special_stt (flag,cmd,tra,ntim,ncol,`	2	`SUBROUTINE special_stt (flag,cmd,tra,ntim,ncol,`
3	`> vars,times,param,nparam,outfile)`	3	`> vars,times,param,nparam,outfile)`
4		4
5	`c ***************************************************************************`	5	`c ***************************************************************************`
6	`c * *`	6	`c * *`
7	`c * OUTPUT: flag -> 1 if trajectory is selected, 0 if not *`	7	`c * OUTPUT: flag -> 1 if trajectory is selected, 0 if not *`
8	`c * *`	8	`c * *`
9	`c * INPUT: cmd <- command string *`	9	`c * INPUT: cmd <- command string *`
10	`c * tra(ntim,ncol) <- single trajectory: indices time,column *`	10	`c * tra(ntim,ncol) <- single trajectory: indices time,column *`
11	`c * ntim <- number of times *`	11	`c * ntim <- number of times *`
12	`c * ncol <- number of columns (including time,lon,lat,p) *`	12	`c * ncol <- number of columns (including time,lon,lat,p) *`
13	`c * vars(ncol) <- names of columns *`	13	`c * vars(ncol) <- names of columns *`
14	`c * times(ntim) <- List of times`	14	`c * times(ntim) <- List of times`
15	`c * param(nparam) <- parameter values *`	15	`c * param(nparam) <- parameter values *`
16	`c * nparam <- number of parameters *`	16	`c * nparam <- number of parameters *`
17	`c * *`	17	`c * *`
18	`c ***************************************************************************`	18	`c ***************************************************************************`
19		19
20	`implicit none`	20	`implicit none`
21		21
22	`c ---------------------------------------------------------------------------`	22	`c ---------------------------------------------------------------------------`
23	`c Declaration of subroutine parameters`	23	`c Declaration of subroutine parameters`
24	`c ---------------------------------------------------------------------------`	24	`c ---------------------------------------------------------------------------`
25		25
26	`integer flag ! Boolean flag whether trajectory is selected`	26	`integer flag ! Boolean flag whether trajectory is selected`
27	`character*80 cmd ! Command string`	27	`character*80 cmd ! Command string`
28	`integer ntim,ncol ! Dimension of single trajectory`	28	`integer ntim,ncol ! Dimension of single trajectory`
29	`real tra(ntim,ncol) ! Single trajectory`	29	`real tra(ntim,ncol) ! Single trajectory`
30	`character*80 vars(ncol) ! Name of columns`	30	`character*80 vars(ncol) ! Name of columns`
31	`real times(ntim) ! List of times`	31	`real times(ntim) ! List of times`
32	`integer nparam ! # parameters`	32	`integer nparam ! # parameters`
33	`real param(1000) ! List of parameters`	33	`real param(1000) ! List of parameters`
34	`character*80 outfile ! NAme of output file`	34	`character*80 outfile ! NAme of output file`
35		35
36	`integer i,n ! Loop indices`	36	`integer i,n ! Loop indices`
37		37
38	`c ---------------------------------------------------------------------------`	38	`c ---------------------------------------------------------------------------`
39	`c Local variables for SPECIAL:WCB`	39	`c Local variables for SPECIAL:WCB`
40	`c ---------------------------------------------------------------------------`	40	`c ---------------------------------------------------------------------------`
41		41
42	`integer ip,i0,i1`	42	`integer ip,i0,i1`
43		43
44	`c ---------------------------------------------------------------------------`	44	`c ---------------------------------------------------------------------------`
45	`c Local variables for SPECIAL:STT`	45	`c Local variables for SPECIAL:STT`
46	`c ---------------------------------------------------------------------------`	46	`c ---------------------------------------------------------------------------`
47		47
48	`integer debugflag`	48	`integer debugflag`
49	`integer tb,ta !tb and ta are the minimum residence times before and after the exchange`	49	`integer tb,ta !tb and ta are the minimum residence times before and after the exchange`
50	`integer stat, ex, outputflag !Variables for I/O files`	50	`integer stat, ex, outputflag !Variables for I/O files`
51	`integer labelflag !Flag to handle label field`	51	`integer labelflag !Flag to handle label field`
52	`integer pvflag, thflag, zflag, pflag !Flags to signal presence of various fields`	52	`integer pvflag, thflag, zflag, pflag !Flags to signal presence of various fields`
53	`integer pvcol,thcol,labelcol,zcol,pcol !Positions of PV, TH, LABEL, z and p columns`	53	`integer pvcol,thcol,labelcol,zcol,pcol !Positions of PV, TH, LABEL, z and p columns`
54	`real mean1, mean2, sd1, sd2 !Mean and standard deviations of PV and TH`	54	`real mean1, mean2, sd1, sd2 !Mean and standard deviations of PV and TH`
55	`real deltat !Trajectory timestep`	55	`real deltat !Trajectory timestep`
56	`real extime, exlon, exlat, exth, exheight, expres!Position at exchange time`	56	`real extime, exlon, exlat, exth, exheight, expres!Position at exchange time`
57	`character*80 namefile !Name of output file`	57	`character*80 namefile !Name of output file`
58	`real a,b !Normalized tropopause distance for PV and TH`	58	`real a,b !Normalized tropopause distance for PV and TH`
59	`real thtrop, pvtrop !Tropopause level in K and pvu`	59	`real thtrop, pvtrop !Tropopause level in K and pvu`
60	`integer count1, count2`	60	`integer count1, count2`
61	`integer sttcount, p !Ancillary variables in STT do while loop`	61	`integer sttcount, p !Ancillary variables in STT do while loop`
62	`real,allocatable, dimension (:) :: tropdist ! Vector for normalized distance from tropopause`	62	`real,allocatable, dimension (:) :: tropdist ! Vector for normalized distance from tropopause`
63	`real,allocatable, dimension (:) :: exflag ! Vector for normalized distance from tropopause`	63	`real,allocatable, dimension (:) :: exflag ! Vector for normalized distance from tropopause`
64	`integer, dimension (ntim) :: pos ! Vector to store exchange indices`	64	`integer, dimension (ntim) :: pos ! Vector to store exchange indices`
65		65
66		66
67		67
68	`c -------------------------------------------------------------------------- %)`	68	`c -------------------------------------------------------------------------- %)`
69	`c SPECIAL:STT:pvtrop,thtrop,tmin_b,tmin_f,outputflag %)`	69	`c SPECIAL:STT:pvtrop,thtrop,tmin_b,tmin_f,outputflag %)`
70		70
71	`c --------------------------------------------------------------------------- %)`	71	`c --------------------------------------------------------------------------- %)`
72		72
73	`if ( cmd.eq.'STT' ) then`	73	`if ( cmd.eq.'STT' ) then`
74		74
75	`c Reset the flag for selection`	75	`c Reset the flag for selection`
76	`flag = 0`	76	`flag = 0`
77	`c Reset the flag for debugging`	77	`c Reset the flag for debugging`
78	`debugflag = 0`	78	`debugflag = 0`
79	`c Set other flags`	79	`c Set other flags`
80	`labelflag = 0`	80	`labelflag = 0`
81	`pvflag = 0`	81	`pvflag = 0`
82	`thflag = 0`	82	`thflag = 0`
83	`C Set data at the exchange`	83	`C Set data at the exchange`
84	`extime = 0`	84	`extime = 0`
85	`exlon = 0`	85	`exlon = 0`
86	`exlat = 0`	86	`exlat = 0`
87	`exheight = 0`	87	`exheight = 0`
88	`expres = 0`	88	`expres = 0`
89		89
90	`c PV and TH threshold read from command line`	90	`c PV and TH threshold read from command line`
91	`pvtrop = param(1)`	91	`pvtrop = param(1)`
92	`thtrop = param(2)`	92	`thtrop = param(2)`
93		93
94		94
95	`c Determine number of time steps in the trajectory`	95	`c Determine number of time steps in the trajectory`
96	`c Handles backward OR forward trajectories`	96	`c Handles backward OR forward trajectories`
97		97
98	`if ( (ntim.eq.0).or.(ntim.eq.1) ) then`	98	`if ( (ntim.eq.0).or.(ntim.eq.1) ) then`
99	`print*,' ERROR: zero length trajectory in SPECIAL:STT... Stop' !If there are problems, stop computation.`	99	`print*,' ERROR: zero length trajectory in SPECIAL:STT... Stop' !If there are problems, stop computation.`
100	`stop`	100	`stop`
101	`endif`	101	`endif`
102		102
103	`deltat = (times(ntim)-times(1))/(ntim-1) !By default array indexing starts in Fortran at 1`	103	`deltat = (times(ntim)-times(1))/(ntim-1) !By default array indexing starts in Fortran at 1`
104		104
105		105
106	`C Now some checks for essential fields : label, potential vorticity, potential temperature`	106	`C Now some checks for essential fields : label, potential vorticity, potential temperature`
107		107
108	`c Check for the presence of the LABEL field`	108	`c Check for the presence of the LABEL field`
109	`labelcol = 1`	109	`labelcol = 1`
110	`do i = 1,ncol`	110	`do i = 1,ncol`
111	`if((vars(i) .eq. 'LABEL').or.(vars(i) .eq. 'label')) then`	111	`if((vars(i) .eq. 'LABEL').or.(vars(i) .eq. 'label')) then`
112	`labelflag = 1`	112	`labelflag = 1`
113	`labelcol = i`	113	`labelcol = i`
114	`if(debugflag .eq. 1) then`	114	`if(debugflag .eq. 1) then`
115	`print*, 'Found label'`	115	`print*, 'Found label'`
116	`endif`	116	`endif`
117	`endif`	117	`endif`
118	`enddo`	118	`enddo`
119		119
120		120
121	`c Check for the presence of the variables PV and TH`	121	`c Check for the presence of the variables PV and TH`
122		122
123	`pvcol = 1`	123	`pvcol = 1`
124	`thcol = 1`	124	`thcol = 1`
125	`zcol = 1`	125	`zcol = 1`
126	`pcol = 1 !In case th or pv are not available, tra(timestep,thcol) is always defined`	126	`pcol = 1 !In case th or pv are not available, tra(timestep,thcol) is always defined`
127		127
128		128
129	`do i = 1,ncol`	129	`do i = 1,ncol`
130	`if((vars(i) .eq. 'PV') .or. (vars(i) .eq. 'POT_VORTI')) then`	130	`if((vars(i) .eq. 'PV') .or. (vars(i) .eq. 'POT_VORTI')) then`
131	`pvflag = 1`	131	`pvflag = 1`
132	`pvcol = i`	132	`pvcol = i`
133	`endif`	133	`endif`
134	`enddo`	134	`enddo`
135		135
136		136
137	`do i = 1,ncol`	137	`do i = 1,ncol`
138	`if((vars(i) .eq. 'TH') .or. (vars(i) .eq. 'POT_T')) then`	138	`if((vars(i) .eq. 'TH') .or. (vars(i) .eq. 'POT_T')) then`
139	`thflag = 1`	139	`thflag = 1`
140	`thcol = i`	140	`thcol = i`
141	`endif`	141	`endif`
142	`enddo`	142	`enddo`
143		143
144	`do i = 1,ncol`	144	`do i = 1,ncol`
145	`if((vars(i) .eq. 'z') .or. (vars(i) .eq. 'Z')) then`	145	`if((vars(i) .eq. 'z') .or. (vars(i) .eq. 'Z')) then`
146	`zcol = i`	146	`zcol = i`
147	`zflag = 1`	147	`zflag = 1`
148	`endif`	148	`endif`
149	`enddo`	149	`enddo`
150		150
151	`do i = 1,ncol`	151	`do i = 1,ncol`
152	`if((vars(i) .eq. 'p') .or. (vars(i) .eq. 'P')) then`	152	`if((vars(i) .eq. 'p') .or. (vars(i) .eq. 'P')) then`
153	`pcol = i`	153	`pcol = i`
154	`pflag = 1`	154	`pflag = 1`
155	`endif`	155	`endif`
156	`enddo`	156	`enddo`
157		157
158	`c No PV and THETA available: not possible determining tropopause. Program stop.`	158	`c No PV and THETA available: not possible determining tropopause. Program stop.`
159		159
160	`if ( (thflag .eq. 0).and.(pvflag .eq. 0) ) then`	160	`if ( (thflag .eq. 0).and.(pvflag .eq. 0) ) then`
161	`print*,' ERROR: no PV and no TH in SPECIAL:STT..'`	161	`print*,' ERROR: no PV and no TH in SPECIAL:STT..'`
162	`print*,' ERROR: not possible determining tropopause'`	162	`print*,' ERROR: not possible determining tropopause'`
163	`stop`	163	`stop`
164	`endif`	164	`endif`
165		165
166	`c Check for missing data in the trajectories`	166	`c Check for missing data in the trajectories`
167		167
168	`do n=1,ntim`	168	`do n=1,ntim`
169	`if ((tra(n,pvcol).lt.-500.).and.(pvflag.eq.1.)) then`	169	`if ((tra(n,pvcol).lt.-500.).and.(pvflag.eq.1.)) then`
170	`if(debugflag .eq. 1) then`	170	`if(debugflag .eq. 1) then`
171	`print*,'WARNING: PV missing data',tra(n,pvcol)`	171	`print*,'WARNING: PV missing data',tra(n,pvcol)`
172	`print*,'Trajectory skipped'`	172	`print*,'Trajectory skipped'`
173	`endif`	173	`endif`
174	`goto 600`	174	`goto 600`
175	`endif`	175	`endif`
176	`if ((tra(n,thcol) .lt.-500.).and.(thflag .eq. 1.)) then`	176	`if ((tra(n,thcol) .lt.-500.).and.(thflag .eq. 1.)) then`
177	`if(debugflag .eq. 1) then`	177	`if(debugflag .eq. 1) then`
178	`print*,'WARNING: TH missing data',tra(n,thcol)`	178	`print*,'WARNING: TH missing data',tra(n,thcol)`
179	`print*,'Trajectory skipped'`	179	`print*,'Trajectory skipped'`
180	`endif`	180	`endif`
181	`goto 600`	181	`goto 600`
182	`endif`	182	`endif`
183	`enddo`	183	`enddo`
184		184
185	`c Computing mean of PV and TH along a trajectory`	185	`c Computing mean of PV and TH along a trajectory`
186		186
187	`mean1=0.`	187	`mean1=0.`
188	`count1=0`	188	`count1=0`
189	`mean2=0.`	189	`mean2=0.`
190	`count2=0`	190	`count2=0`
191	`do n=1,ntim`	191	`do n=1,ntim`
192	`if ((tra(n,pvcol).gt.-100.).and.(pvflag .eq. 1)) then`	192	`if ((tra(n,pvcol).gt.-100.).and.(pvflag .eq. 1)) then`
193	`count1=count1+1`	193	`count1=count1+1`
194	`mean1=mean1+tra(n,pvcol)`	194	`mean1=mean1+tra(n,pvcol)`
195	`else`	195	`else`
196	`if(debugflag .eq. 1) then`	196	`if(debugflag .eq. 1) then`
197	`print*,'WARNING: PV is missing or non-physical, mean1 = 0'`	197	`print*,'WARNING: PV is missing or non-physical, mean1 = 0'`
198	`endif`	198	`endif`
199	`mean1 = 0`	199	`mean1 = 0`
200	`endif`	200	`endif`
201	`if ((tra(n,thcol) .gt. 1E-6) .and. (thflag .eq. 1)) then`	201	`if ((tra(n,thcol) .gt. 1E-6) .and. (thflag .eq. 1)) then`
202	`count2=count2+1`	202	`count2=count2+1`
203	`mean2=mean2+tra(n,thcol)`	203	`mean2=mean2+tra(n,thcol)`
204	`else`	204	`else`
205	`if(debugflag .eq. 1) then`	205	`if(debugflag .eq. 1) then`
206	`print*,'WARNING: TH is missing or non-physical, mean1 = 0'`	206	`print*,'WARNING: TH is missing or non-physical, mean1 = 0'`
207	`endif`	207	`endif`
208	`mean2 = 0`	208	`mean2 = 0`
209	`endif`	209	`endif`
210	`enddo`	210	`enddo`
211		211
212	`c Compute average, otherwise skip problematic trajectories`	212	`c Compute average, otherwise skip problematic trajectories`
213		213
214	`if (count1.ne.0) then`	214	`if (count1.ne.0) then`
215	`mean1=mean1/count1`	215	`mean1=mean1/count1`
216	`endif`	216	`endif`
217		217
218	`if (count2.ne.0) then`	218	`if (count2.ne.0) then`
219	`mean2=mean2/count2`	219	`mean2=mean2/count2`
220	`endif`	220	`endif`
221		221
222	`if ( (count1.eq.0) .and. (count2.eq.0) )then`	222	`if ( (count1.eq.0) .and. (count2.eq.0) )then`
223	`if(debugflag .eq. 1) then`	223	`if(debugflag .eq. 1) then`
224	`print*, 'Trajectory with only missing data skipped'`	224	`print*, 'Trajectory with only missing data skipped'`
225	`endif`	225	`endif`
226	`goto 600`	226	`goto 600`
227	`endif`	227	`endif`
228		228
229	`c Computing standard deviation of PV and TH along a trajectory`	229	`c Computing standard deviation of PV and TH along a trajectory`
230		230
231	`sd1=0.`	231	`sd1=0.`
232	`count1=0.`	232	`count1=0.`
233	`sd2=0.`	233	`sd2=0.`
234	`count2=0.`	234	`count2=0.`
235	`do n=1,ntim`	235	`do n=1,ntim`
236	`if (tra(n,pvcol).gt.-100.) then`	236	`if (tra(n,pvcol).gt.-100.) then`
237	`count1=count1+1.`	237	`count1=count1+1.`
238	`sd1=sd1+(tra(n,pvcol)-mean1)**2`	238	`sd1=sd1+(tra(n,pvcol)-mean1)**2`
239	`endif`	239	`endif`
240	`if (tra(n,thcol).gt. 0.) then`	240	`if (tra(n,thcol).gt. 0.) then`
241	`count2=count2+1.`	241	`count2=count2+1.`
242	`sd2=sd2+(tra(n,thcol)-mean2)**2`	242	`sd2=sd2+(tra(n,thcol)-mean2)**2`
243	`endif`	243	`endif`
244	`enddo`	244	`enddo`
245	`sd1=sqrt(sd1/(count1-1.))`	245	`sd1=sqrt(sd1/(count1-1.))`
246	`sd2=sqrt(sd2/(count2-1.))`	246	`sd2=sqrt(sd2/(count2-1.))`
247		247
248	`c Allocate distance array`	248	`c Allocate distance array`
249		249
250	`allocate(tropdist(1:ntim),stat=stat)`	250	`allocate(tropdist(1:ntim),stat=stat)`
251	`if (stat.ne.0) then`	251	`if (stat.ne.0) then`
252	`print*,'ERROR: allocating array tropdist`	252	`print*,'ERROR: allocating array tropdist`
253	`&not possible in SPECIAL:STT.. '`	253	`&not possible in SPECIAL:STT.. '`
254	`stop`	254	`stop`
255	`endif`	255	`endif`
256		256
257	`c If LABEL is not available, go to section without usage of label`	257	`c If LABEL is not available, go to section without usage of label`
258		258
259	`if(labelflag .ne. 1) goto 300`	259	`if(labelflag .ne. 1) goto 300`
260		260
261		261
262	`C-------------------------------------------------------------------------------------------------------------------------`	262	`C-------------------------------------------------------------------------------------------------------------------------`
263	`c BEGINNING OF THE SECTION WITH LABEL`	263	`c BEGINNING OF THE SECTION WITH LABEL`
264	`C-------------------------------------------------------------------------------------------------------------------------`	264	`C-------------------------------------------------------------------------------------------------------------------------`
265		265
266		266
267	`c Check for label field`	267	`c Check for label field`
268	`do i = 1,ntim`	268	`do i = 1,ntim`
269	`if((tra(i,labelcol) .ne.1) .and. (tra(i,labelcol) .ne.2)) then`	269	`if((tra(i,labelcol) .ne.1) .and. (tra(i,labelcol) .ne.2)) then`
270	`if(debugflag .eq. 1) then`	270	`if(debugflag .eq. 1) then`
271	`print*,'WARNING: unexpected label field in SPECIAL:STT.. '`	271	`print*,'WARNING: unexpected label field in SPECIAL:STT.. '`
272	`endif`	272	`endif`
273	`print*,'WARNING: unused label field in SPECIAL:STT'`	273	`print*,'WARNING: unused label field in SPECIAL:STT'`
274	`goto 300`	274	`goto 300`
275	`endif`	275	`endif`
276	`end do`	276	`end do`
277		277
278		278
279	`c Distance from the '2 PVU, 380 K' tropopause (-:trop/+:strat) : pvtrop, thtrop if PV and TH and LABEL are available`	279	`c Distance from the '2 PVU, 380 K' tropopause (-:trop/+:strat) : pvtrop, thtrop if PV and TH and LABEL are available`
280		280
281	`if ( (pvflag .ne. 0) .and. (thflag .ne. 0) .and.`	281	`if ( (pvflag .ne. 0) .and. (thflag .ne. 0) .and.`
282	`&(labelflag .eq. 1)) then`	282	`&(labelflag .eq. 1)) then`
283	`if(debugflag .eq. 1) then`	283	`if(debugflag .eq. 1) then`
284	`print*,'Tropopause detected with PV, TH and LABEL'`	284	`print*,'Tropopause detected with PV, TH and LABEL'`
285	`endif`	285	`endif`
286	`do n=1,ntim`	286	`do n=1,ntim`
287	`a=abs(pvtrop-abs(tra(n,pvcol)))/sd1 ! a = tropdist(PV)`	287	`a=abs(pvtrop-abs(tra(n,pvcol)))/sd1 ! a = tropdist(PV)`
288	`b=abs(thtrop-tra(n,thcol))/sd2 ! b = tropdist(TH)`	288	`b=abs(thtrop-tra(n,thcol))/sd2 ! b = tropdist(TH)`
289	`if (tra(n,labelcol) .eq. 1) then ! Troposphere (LABEL)`	289	`if (tra(n,labelcol) .eq. 1) then ! Troposphere (LABEL)`
290	`tropdist(n)=-min(a,b)`	290	`tropdist(n)=-min(a,b)`
291	`else`	291	`else`
292	`tropdist(n)=+min(a,b)`	292	`tropdist(n)=+min(a,b)`
293	`endif`	293	`endif`
294	`enddo`	294	`enddo`
295	`endif`	295	`endif`
296		296
297	`c If PV is not available, use only TH`	297	`c If PV is not available, use only TH`
298		298
299	`if (pvflag .eq. 0) then`	299	`if (pvflag .eq. 0) then`
300	`if(debugflag .eq. 1) then`	300	`if(debugflag .eq. 1) then`
301	`print*,'Tropopause detected with TH only'`	301	`print*,'Tropopause detected with TH only'`
302	`endif`	302	`endif`
303	`do n=1,ntim`	303	`do n=1,ntim`
304	`b=abs(thtrop-tra(n,thcol))/sd2 ! b = tropdist(TH)`	304	`b=abs(thtrop-tra(n,thcol))/sd2 ! b = tropdist(TH)`
305	`if(tra(n,labelcol) .eq. 1) then ! Troposphere (TH)`	305	`if(tra(n,labelcol) .eq. 1) then ! Troposphere (TH)`
306	`tropdist(n)= -b`	306	`tropdist(n)= -b`
307	`else`	307	`else`
308	`tropdist(n)= b ! Stratosphere (TH)`	308	`tropdist(n)= b ! Stratosphere (TH)`
309	`endif`	309	`endif`
310	`enddo`	310	`enddo`
311	`endif`	311	`endif`
312		312
313	`c If TH is not available, use only PV`	313	`c If TH is not available, use only PV`
314	`if (thflag .eq. 0) then`	314	`if (thflag .eq. 0) then`
315	`if(debugflag .eq. 1) then`	315	`if(debugflag .eq. 1) then`
316	`print*,'Tropopause detected with PV only'`	316	`print*,'Tropopause detected with PV only'`
317	`endif`	317	`endif`
318	`do n=1,ntim`	318	`do n=1,ntim`
319	`a=abs(pvtrop-tra(n,pvcol))/sd2 ! a = tropdist(PV)`	319	`a=abs(pvtrop-tra(n,pvcol))/sd2 ! a = tropdist(PV)`
320	`if( (tra(n,pvcol).lt.pvtrop) .and.`	320	`if( (tra(n,pvcol).lt.pvtrop) .and.`
321	`&(tra(n,labelcol) .eq. 1)) then ! Troposphere (PV)`	321	`&(tra(n,labelcol) .eq. 1)) then ! Troposphere (PV)`
322	`tropdist(n)= -a`	322	`tropdist(n)= -a`
323	`else`	323	`else`
324	`tropdist(n)= a ! Stratosphere (PV)`	324	`tropdist(n)= a ! Stratosphere (PV)`
325	`endif`	325	`endif`
326	`enddo`	326	`enddo`
327	`endif`	327	`endif`
328		328
329	`C Detecting STT with LABEL field`	329	`C Detecting STT with LABEL field`
330	`c If PV is lower than pvtrop despite the label is 2, then in next step this potential STT is not considered.`	330	`c If PV is lower than pvtrop despite the label is 2, then in next step this potential STT is not considered.`
331		331
332	`sttcount = 0`	332	`sttcount = 0`
333		333
334	`do i=2,ntim`	334	`do i=2,ntim`
335	`if((tropdist(i)*tropdist(i-1) .lt. 0) .and.`	335	`if((tropdist(i)*tropdist(i-1) .lt. 0) .and.`
336	`&(tropdist(i-1) .gt. 0)) then`	336	`&(tropdist(i-1) .gt. 0)) then`
337	`C Adding label check for STT`	337	`C Adding label check for STT`
338	`if( (tra(i,labelcol) .ne. 1) .and. (tra(i-1,labelcol) .ne. 2)) then`	338	`if((tra(i,labelcol).ne.1).and.(tra(i-1,labelcol).ne.2))then`
339	`if(debugflag .eq. 1) then`	339	`if(debugflag .eq. 1) then`
340	`print*, 'STT label check failed:'`	340	`print*, 'STT label check failed:'`
341	`print*, 'moving to next STT candidate'`	341	`print*, 'moving to next STT candidate'`
342	`endif`	342	`endif`
343	`cycle`	343	`cycle`
344	`else`	344	`else`
345	`sttcount = sttcount+1`	345	`sttcount = sttcount+1`
346	`pos(sttcount) = i-1`	346	`pos(sttcount) = i-1`
347	`endif`	347	`endif`
348	`endif`	348	`endif`
349	`enddo`	349	`enddo`
350		350
351	`c If no STT is detected, move to next trajectory`	351	`c If no STT is detected, move to next trajectory`
352		352
353	`if(sttcount .lt. 1 ) then`	353	`if(sttcount .lt. 1 ) then`
354	`if(debugflag .eq. 1) then`	354	`if(debugflag .eq. 1) then`
355	`print*, 'No STT detected: moving to next trajectory'`	355	`print*, 'No STT detected: moving to next trajectory'`
356	`endif`	356	`endif`
357	`goto 600`	357	`goto 600`
358	`endif`	358	`endif`
359		359
360	`c Re-connection with normal steps: label has been used`	360	`c Re-connection with normal steps: label has been used`
361		361
362	`goto 350`	362	`goto 350`
363		363
364	`C-------------------------------------------------------------------------------------------------------------------------`	364	`C-------------------------------------------------------------------------------------------------------------------------`
365	`C END OF LABEL SECTION`	365	`C END OF LABEL SECTION`
366	`c-------------------------------------------------------------------------------------------------------------------------`	366	`c-------------------------------------------------------------------------------------------------------------------------`
367		367
368		368
369		369
370	`c Distance from the '2 PVU, 380 K' tropopause (-:trop/+:strat) : pvtrop, thtrop if PV and TH are available`	370	`c Distance from the '2 PVU, 380 K' tropopause (-:trop/+:strat) : pvtrop, thtrop if PV and TH are available`
371	`300 continue !Here the code continues if no LABEL has been found`	371	`300 continue !Here the code continues if no LABEL has been found`
372		372
373	`if ( (pvflag .ne. 0) .and. (thflag .ne. 0)) then`	373	`if ( (pvflag .ne. 0) .and. (thflag .ne. 0)) then`
374	`if(debugflag .eq. 1) then`	374	`if(debugflag .eq. 1) then`
375	`print*,'Tropopause detected with PV and TH'`	375	`print*,'Tropopause detected with PV and TH'`
376	`endif`	376	`endif`
377	`do n=1,ntim`	377	`do n=1,ntim`
378	`a=abs(pvtrop-abs(tra(n,pvcol)))/sd1 ! a = tropdist(PV)`	378	`a=abs(pvtrop-abs(tra(n,pvcol)))/sd1 ! a = tropdist(PV)`
379	`b=abs(thtrop-tra(n,thcol))/sd2 ! b = tropdist(TH)`	379	`b=abs(thtrop-tra(n,thcol))/sd2 ! b = tropdist(TH)`
380	`if ((tra(n,thcol).lt.thtrop).and.`	380	`if ((tra(n,thcol).lt.thtrop).and.`
381	`&(abs(tra(n,pvcol)).lt.pvtrop)) then ! Troposphere (TH and PV) ---------WHY ABS??---------`	381	`&(abs(tra(n,pvcol)).lt.pvtrop)) then ! Troposphere (TH and PV) ---------WHY ABS??---------`
382	`tropdist(n)=-min(a,b)`	382	`tropdist(n)=-min(a,b)`
383	`else`	383	`else`
384	`if ( (abs(tra(n,pvcol)).ge.pvtrop).and.`	384	`if ( (abs(tra(n,pvcol)).ge.pvtrop).and.`
385	`&(tra(n,thcol).ge.thtrop) ) then ! Stratosphere (TH and PV)`	385	`&(tra(n,thcol).ge.thtrop) ) then ! Stratosphere (TH and PV)`
386	`tropdist(n)=+min(a,b)`	386	`tropdist(n)=+min(a,b)`
387	`endif`	387	`endif`
388		388
389	`if ( (abs(tra(n,pvcol)).ge.pvtrop).and.`	389	`if ( (abs(tra(n,pvcol)).ge.pvtrop).and.`
390	`&(tra(n,thcol).lt.thtrop) ) then ! Stratosphere (PV)`	390	`&(tra(n,thcol).lt.thtrop) ) then ! Stratosphere (PV)`
391	`tropdist(n)=a`	391	`tropdist(n)=a`
392	`endif`	392	`endif`
393		393
394	`if ( (abs(tra(n,pvcol)).lt.pvtrop).and.`	394	`if ( (abs(tra(n,pvcol)).lt.pvtrop).and.`
395	`&(tra(n,thcol).ge.thtrop) ) then ! Stratosphere (TH)`	395	`&(tra(n,thcol).ge.thtrop) ) then ! Stratosphere (TH)`
396	`tropdist(n)=b`	396	`tropdist(n)=b`
397	`endif`	397	`endif`
398	`endif`	398	`endif`
399	`enddo`	399	`enddo`
400	`endif`	400	`endif`
401		401
402	`c If PV is not available, use only TH`	402	`c If PV is not available, use only TH`
403		403
404	`if (pvflag .eq. 0) then`	404	`if (pvflag .eq. 0) then`
405	`if(debugflag .eq. 1) then`	405	`if(debugflag .eq. 1) then`
406	`print*,'Tropopause detected with TH only'`	406	`print*,'Tropopause detected with TH only'`
407	`endif`	407	`endif`
408	`do n=1,ntim`	408	`do n=1,ntim`
409	`b=abs(thtrop-tra(n,thcol))/sd2 ! b = tropdist(TH)`	409	`b=abs(thtrop-tra(n,thcol))/sd2 ! b = tropdist(TH)`
410	`if (tra(n,thcol).lt.thtrop) then ! Troposphere (TH)`	410	`if (tra(n,thcol).lt.thtrop) then ! Troposphere (TH)`
411	`tropdist(n)= -b`	411	`tropdist(n)= -b`
412	`else`	412	`else`
413	`tropdist(n)= b ! Stratosphere (TH)`	413	`tropdist(n)= b ! Stratosphere (TH)`
414	`endif`	414	`endif`
415	`enddo`	415	`enddo`
416	`endif`	416	`endif`
417		417
418	`c If TH is not available, use only PV`	418	`c If TH is not available, use only PV`
419	`if(debugflag .eq. 1) then`	419	`if(debugflag .eq. 1) then`
420	`print*,'Tropopause detected with PV only'`	420	`print*,'Tropopause detected with PV only'`
421	`endif`	421	`endif`
422	`if (thflag .eq. 0) then`	422	`if (thflag .eq. 0) then`
423	`do n=1,ntim`	423	`do n=1,ntim`
424	`a=abs(pvtrop-tra(n,pvcol))/sd2 ! a = tropdist(PV)`	424	`a=abs(pvtrop-tra(n,pvcol))/sd2 ! a = tropdist(PV)`
425	`if (tra(n,pvcol).lt.pvtrop) then ! Troposphere (PV)`	425	`if (tra(n,pvcol).lt.pvtrop) then ! Troposphere (PV)`
426	`tropdist(n)= -a`	426	`tropdist(n)= -a`
427	`else`	427	`else`
428	`tropdist(n)= a ! Stratosphere (PV)`	428	`tropdist(n)= a ! Stratosphere (PV)`
429	`endif`	429	`endif`
430	`enddo`	430	`enddo`
431	`endif`	431	`endif`
432		432
433		433
434	`c Detect the various STT occurrences and save positions in traj file in array pos`	434	`c Detect the various STT occurrences and save positions in traj file in array pos`
435		435
436	`sttcount = 0`	436	`sttcount = 0`
437		437
438	`do i=2,ntim`	438	`do i=2,ntim`
439	`if((tropdist(i)*tropdist(i-1) .lt. 0) .and.`	439	`if((tropdist(i)*tropdist(i-1) .lt. 0) .and.`
440	`& (tropdist(i-1) .gt. 0)) then`	440	`& (tropdist(i-1) .gt. 0)) then`
441	`sttcount = sttcount+1`	441	`sttcount = sttcount+1`
442	`pos(sttcount) = i-1`	442	`pos(sttcount) = i-1`
443	`endif`	443	`endif`
444	`enddo`	444	`enddo`
445		445
446	`c If no STT are detected, go to next trajectory`	446	`c If no STT are detected, go to next trajectory`
447		447
448	`if(sttcount .lt. 1 ) then`	448	`if(sttcount .lt. 1 ) then`
449	`if(debugflag .eq. 1) then`	449	`if(debugflag .eq. 1) then`
450	`print*, 'No STT detected: moving to next trajectory'`	450	`print*, 'No STT detected: moving to next trajectory'`
451	`endif`	451	`endif`
452	`goto 600`	452	`goto 600`
453	`endif`	453	`endif`
454		454
455		455
456		456
457	`350 continue !The same if LABEL is present or not.`	457	`350 continue !The same if LABEL is present or not.`
458		458
459	`c In timesteps, minimum residence time before and after the exchange`	459	`c In timesteps, minimum residence time before and after the exchange`
460		460
461	`tb = nint( param(3)/deltat )`	461	`tb = nint( param(3)/deltat )`
462	`ta = nint( param(4)/deltat )`	462	`ta = nint( param(4)/deltat )`
463		463
464		464
465	`c Now a big while loop analyzing every candidate STT`	465	`c Now a big while loop analyzing every candidate STT`
466		466
467	`do while(sttcount>0)`	467	`do while(sttcount>0)`
468		468
469	`p = pos(sttcount)`	469	`p = pos(sttcount)`
470		470
471	`c Check whether the trajectory is still in the troposphere after ta timesteps`	471	`c Check whether the trajectory is still in the troposphere after ta timesteps`
472		472
473	`if(ta .ne. 0) then`	473	`if(ta .ne. 0) then`
474	`do n=1,ta`	474	`do n=1,ta`
475	`if((p+n) .le. ntim) then`	475	`if((p+n) .le. ntim) then`
476	`if (tropdist(p)*tropdist(p+n).gt.0.) goto 500 !Try with new STT candidate`	476	`if (tropdist(p)*tropdist(p+n).gt.0.) goto 500 !Try with new STT candidate`
477	`else`	477	`else`
478	`goto 500 !Out of range. Try with new STT candidate`	478	`goto 500 !Out of range. Try with new STT candidate`
479	`endif`	479	`endif`
480	`enddo`	480	`enddo`
481	`endif`	481	`endif`
482		482
483	`c Check whether the trajectory had been in the stratosphere for at least tb timesteps`	483	`c Check whether the trajectory had been in the stratosphere for at least tb timesteps`
484		484
485	`if(tb .ne. 0) then`	485	`if(tb .ne. 0) then`
486	`do n=1,tb`	486	`do n=1,tb`
487	`if((p-n) .gt. 0) then`	487	`if((p-n) .gt. 0) then`
488	`if (tropdist(p)*tropdist(p-n).lt.0.) goto 500 !Try with new STT candidate`	488	`if (tropdist(p)*tropdist(p-n).lt.0.) goto 500 !Try with new STT candidate`
489	`else`	489	`else`
490	`goto 500 !Out of range. Try with new STT candidate`	490	`goto 500 !Out of range. Try with new STT candidate`
491	`endif`	491	`endif`
492	`enddo`	492	`enddo`
493	`endif`	493	`endif`
494		494
495	`c If the loops have gone fine, the trajectory is considered, flag = 1, values at the exchange are written and the loop stops.`	495	`c If the loops have gone fine, the trajectory is considered, flag = 1, values at the exchange are written and the loop stops.`
496		496
497	`flag = 1`	497	`flag = 1`
498	`extime = tra(p,1)`	498	`extime = tra(p,1)`
499	`exlon = tra(p,2)`	499	`exlon = tra(p,2)`
500	`exlat = tra(p,3)`	500	`exlat = tra(p,3)`

Subversion Repositories lagranto.icon

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