
      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
