Subversion Repositories lagranto.icon

Rev

Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed 


      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