C-----------------------------------------------------------------------
C----------------------------------------------------------------------- C
C       DROG3DDT  -  ***  DROG3DDT.f  
C                    ***  TIME-DEPENDENT VELOCITY FIELDS
C                    ***  PASSIVE BEHAVIORAL CODING DEFAULT;
C                    ***     SAMPLE BEHAVIOR SUPPLIED
C
C                    ***  BRIAN O. BLANTON 
C                    ***  UNIVERSITY OF NORTH CAROLINA - CHAPEL HILL
C                    ***  DEPARTMENT OF MARINE SCIENCES
C                    ***  12-7 VENABLE HALL
C                    ***  CB #3300
C                    ***  CHAPEL HILL, NORTH CAROLINA
C                    ***  27599-3300
C                    ***  
C                    ***  FALL 1998
C
C-----------------------------------------------------------------------
C
C   AUTHORS OF 2-DIMENSIONAL CODE DEVELOPMENT :
C          ANTONIO M. BAPTISTA, MIT, OGC
C          DOUGLAS COSLER, MIT
C          PAUL J. TURNER, OGC
C          E. ERIC ADAMS, MIT
C          RICHARD KOSSICK, MIT
C
C   EXPANSION TO 3D OF HARMONIC CODE :
C          BRIAN O. BLANTON, FALL 92, SKIO
C
C   CONVERSION FROM HARMONIC STAND-ALONE TO TIME-DEPENDENT SUBROUTINE
C   OF THE QUODDY SERIES OF MODELS:
C          BRIAN O. BLANTON, SUMMER 94, UNC-CH
C
C   INCLUSION OF ADDITIONAL QUODDY VARIABLES FOR TURBULENT
C   KICKS, AND ADDITIONAL TIMING INFORMATION FOR STAGGERED 
C   DROGUE RELEASES.  THIS IS PART OF THE REAL-TIME DATA ANALYSIS
C   Observational System Simulation Experiments (OSSE)
C          BRIAN O. BLANTON, 13 December, 1999, UNC-CH
C
C      CODE ALTERATIONS/CHANGES:
C        * TIME INPUTS IN SECONDS!!; AS A CONSEQUENCE, TIME VARIABLES
C          NOW LACK THE HOUR/SEC DESIGNATION IN PREVIOUS VERSIONS.
C        * ITER0 ARRAY TO CONTROL STAGGERED DROGUE RELEASE
C          DIMENSIONS 1:MAXDROG; THIS IS COMPUTED BY
C          DRIVER AND PASSED TO DROG3DDT
C        * INCLUSION OF ADDITIONAL QUODDY VARIABLES FOR TURB
C          KICK:
C                Q2,Q2L,RHO,T,S,ENZM
C
C       10        20        30        40        50        60        7012
C23456789012345678901234567890123456789012345678901234567890123456789012 
C-----------------------------------------------------------------------
C
      SUBROUTINE DROG3DDT(NMND,NMEL,NNV,T2Q,T1Q,ITER,ITER0,
     +                    Z2,U2,V2,W2,Q22,Q2L2,RHO2,
     +                    TMP2,SAL2,ENZM2,ENZH2,ENZQ2,
     +                    DOMAINNAME,INDNAME)
C
C-----------------------------------------------------------------------
C
C     SET MAXIMUM DIMENSIONING OF ARRAYS THROUGH EXTERNAL INCLUDE
C
      INCLUDE 'DCMS.DIM'
C      IMPLICIT NONE
C
C *** THE ONLY INTERNAL (TO DROG3DDT) PARAMETER, MAXDROG - MAX NUMBER 
C *** OF DROGUES.   THIS IS THE ONLY PLACE IT IS DEFINED.  IT IS 
C *** PASSED IF IT IS NEEDED.
C
      INTEGER MAXDROG
      PARAMETER(MAXDROG=5000)
C
C-----------------------------------------------------------------------
C
C *** DECLARE INCOMING VARIABLES;  A VARIABLE WITH A "2"
C *** POST-FIXED MEANS A VALUE FOR THE END OF THE TIMESTEP, CORRESPONDING
C *** TO T2Q; I.E, Q22 IUS THE NAME OF THE ARRAY CONTAINING THE QUANTITY
C *** q2 AT THE TIME T2Q
C
C     NMND     - NUMBER OF NODES IN DOMAIN
C     NMEL     - NUMBER OF ELEMENTS IN DOMAIN
C     NNV      - NUMBER OF VERTICAL LEVELS IN DOMAIN
C     T2Q      - TIME AT END OF ITER TIMESTEP, IN SECONDS!!
C     T1Q      - TIME AT START OF ITER TIMESTEP, IN SECONDS!!
C     ITER     - TIME-STEPPER (QOUODDY4) ITERATION NUMBER
C     ITER0    - ITERATION TO START INDIVIDUAL DROGUE a-trackin'
C     Z2       - ARRAY CONTAINING VERTICAL GRID AT TIME T2Q      
C     U2       - U-VELOCITY AT TIME T2Q;    U IN QUODDY4
C     V2       - V-VELOCITY AT TIME T2Q;    V IN QUODDY4
C     W2       - W-VELOCITY AT TIME T2Q;    W IN QUODDY4
C     Q22      - MODEL Q2 AT TIME T2Q;      Q2 IN QUODDY4
C     Q2L2     - MODEL Q2L AT TIME T2Q;     Q2L IN QUODDY4
C     RHO2     - MODEL DENSITY AT TIME T2Q; RHO IN QUODDY4
C     TMP2     - MODEL TEMP AT TIME T2Q;    T IN QUODDY4
C     SAL2     - MODEL SALT AT TIME T2Q;    S IN QUODDY4 
C     ENZM2    - MODEL ENZM AT TIME T2Q;    ENZM IN QUODDY4 
C     ENZH2    - MODEL ENZH AT TIME T2Q;    ENZH IN QUODDY4 
C     ENZQ2    - MODEL ENZQ AT TIME T2Q;    ENZQ IN QUODDY4 
C     DOMAINNAME - NAME OF DOMAIN (GRID); CASNAM IN QUODDY4
C     INDNAME    - NAME OF DROG INPUT FILE; FILDR3 IN QUODDY4
C
      INTEGER NMND,NMEL,NNV
      REAL T2Q,T1Q
      INTEGER ITER,ITER0(1:MAXDROG) ! ITERATION-BASED STAGGERED DROGUE 
C                                   ! RELEASE ARRAY
      REAL Z2(1:NNDIM,*)
      REAL U2(1:NNDIM,*)
      REAL V2(1:NNDIM,*)
      REAL W2(1:NNDIM,*)
      REAL Q22(1:NNDIM,*)
      REAL Q2L2(1:NNDIM,*)
      REAL RHO2(1:NNDIM,*)
      REAL TMP2(1:NNDIM,*)
      REAL SAL2(1:NNDIM,*)
      REAL ENZM2(1:NNDIM,*)
      REAL ENZH2(1:NNDIM,*)
      REAL ENZQ2(1:NNDIM,*)
      CHARACTER*(*) DOMAINNAME,INDNAME
C 
C *** Z1   - ARRAY CONTAINING VERTICAL GRID  AT TIME T1Q    
C *** U1   - U-VELOCITY AT TIME T1Q;        U IN QUODDY4
C *** V1   - V-VELOCITY AT TIME T1Q;        V IN QUODDY4
C *** W1   - W-VELOCITY AT TIME T1Q;        W IN QUODDY4
C *** Q21      - MODEL Q2 AT TIME T1Q;      Q2 IN QUODDY4
C *** Q2L1     - MODEL Q2L AT TIME T1Q;     Q2L IN QUODDY4
C *** RHO1     - MODEL DENSITY AT TIME T1Q; RHO IN QUODDY4
C *** TMP1     - MODEL TEMP AT TIME T1Q;    T IN QUODDY4
C *** SAL1     - MODEL SALT AT TIME T1Q;    S IN QUODDY4
C *** ENZM1    - MODEL ENZM AT TIME T1Q;    ENZM IN QUODDY4
C *** ENZH1    - MODEL ENZH AT TIME T1Q;    ENZH IN QUODDY4
C
      REAL Z1(1:NNDIM,1:NNVDIM)
      REAL U1(1:NNDIM,1:NNVDIM)
      REAL V1(1:NNDIM,1:NNVDIM)
      REAL W1(1:NNDIM,1:NNVDIM)
      REAL Q21(1:NNDIM,1:NNVDIM)
      REAL Q2L1(1:NNDIM,1:NNVDIM)
      REAL RHO1(1:NNDIM,1:NNVDIM)
      REAL TMP1(1:NNDIM,1:NNVDIM)
      REAL SAL1(1:NNDIM,1:NNVDIM)
      REAL ENZM1(1:NNDIM,1:NNVDIM)
      REAL ENZH1(1:NNDIM,1:NNVDIM)
      REAL ENZQ1(1:NNDIM,1:NNVDIM)
C
C *** DROGUE POSITION ARRAYS
C
      REAL XDR(1:MAXDROG),YDR(1:MAXDROG),ZDR(1:MAXDROG)
      INTEGER IDR(1:MAXDROG),JJDR(1:MAXDROG),LLDR(1:MAXDROG)
C
C *** LOCAL DROGUE INFO PASSED TO TRACK
C *** BEFORE CALLING TRACK, DROG3DDT MOVES INFORMATION ABOUT THE DROGUE
C *** CURRENTLY BEING TRACKED INTO SCALAR VARIABLES
C
C *** NACTIVE IS THE NUMBER OF DROGUES STILL ACTIVE!!, not INactive
C
      REAL U0,V0,W0
      REAL X0,Y0,Z0
      INTEGER J0,L0,NDR,JNEW,IBOUN,JTEMP,LTEMP
      REAL XNEW,YNEW,ZNEW
      INTEGER LNEW,NACTIVE
      DATA IBOUN/0/
C
C *** MODEL TIME-STEP, DROGUE TRACKING TIME-STEP
C
      REAL T1DR,T2DR
      REAL DT_TS,DT_DR
C
C *** DROGUE OUTPUT INTERVAL AND DROGUE TRACKING INTERVAL; THE OUTPUT 
C *** ROUTINE WILL BE CALLED ONLY EVERY IPRINT ITERATIONS. DROGUE 
C *** POSITIONS WILL BE UPDATED EVERY ITRACK MODEL DTS.
C
      INTEGER IPRINT,ITRACK
C
C *** TRACKFLAG - INTERNAL TRACKING FLAG; SET TO NO IF ANY INPUT PROBLEMS
C *** ARE ENCOUNTERED SO THAT NO TRACKING WILL BE PERFORMED AND 
C *** QUODDY2 WILL NOT BE TERMINATED DUE TO A PROBLEM IN DROG3DDT. 
C *** TRACKFLAG IS INITIALLY SET TO 'YES' IN THE FIRST IF BLOCK BELOW.
C
      CHARACTER*3 TRACKFLAG
      COMMON TRACKFLAG
C 
C *** LOOP COUNTERS
C
      INTEGER I,II,L,IND
C
C *** SAVE CERTAIN VARIABLES TO AVOID STATIC ARRAY COMPILATION
C *** ON CERTAIN ARCHITECTURES
C
      SAVE U1,V1,W1,Z1,T1DR,T2DR
      SAVE XDR,YDR,ZDR
      SAVE IDR,JJDR,LLDR,NDR
      SAVE ITRACK,IPRINT,NACTIVE
      SAVE DT_TS,DT_DR
C
C *** PARAMETERIZE HORIZOTAL KICKS: 
C *** HTURBFLAG - 0=NO HORZ KICK, 1=INCLUDE HORZ KICK
C *** K_X,K_Y   - East/West and North/South eddy-diffusivity
C *** DEVIATE   - 'GAUSS'=GAUSSIAN, 'UNIFO'=UNIFORM
C *** ALPHA     - ALPHA IS A FACTOR DEPENDING IN DEVIATE CHOICE
C ***             SEE DDOCIMENTATION;  THIS IS SET BY THE SOFTWARE
C *** BETA_X,BETA_Y - East/West and North/South Random deviate 
C ***                 from normal distribution (0,1), "gasdev"
C ***                 from NRFv2.0
C *** KICK_X,KICK_Y - THE ACTUAL HORZ KICKS, COMPUTED BELOW
C 
      INTEGER HTURBFLAG
      PARAMETER (HTURBFLAG=0)
      REAL K_X,K_Y
      PARAMETER (K_X=50.)
      PARAMETER (K_Y=50.) 
CBOB      CHARACTER*5 DEVIATE
CBOB      PARAMETER (DEVIATE='GAUSS')
      REAL ALPHA
      PARAMETER (ALPHA=2.)
      REAL KICK_X,KICK_Y
      REAL BETA_X,BETA_Y
C
C *** VERTICAL KICK VARIABLES: 
C *** VTURBFLAG - 0=NO VERT KICK, 1=INCLUDE VERT KICK
      INTEGER VTURBFLAG
      PARAMETER (VTURBFLAG=0)
C
C *** NRFv2.0 ROUTINES
      REAL RAN1,RAN2,GASDEV
C
C *** MISC DECLARATIONS
C
      LOGICAL BOTFLAG,INDLL
      REAL ZLOW
      REAL ZNOW(1:3,1:NNVDIM)
      REAL XHOR,YHOR
      INTEGER ISEED
      DATA ISEED/-80002/
C
C *** 
C *** INPUT TIMES ARE NOW IN SECONDS  !!!!!!!!!!!!!!
C *** 
C
C
C *** CALL INPUT AND INITIALIZATION ROUTINE BEFORE INITIAL TIMESTEP ONLY
C *** I.E., ITER=0;  THIS MEANS DROG3DDT MUST BE CALLED BEFORE TIME-STEPPING 
C *** BEGINS IN THE TIME-STEPPER.  DROG3DDT MUST BE PASSED A 0 (INTEGER ZERO)
C *** IN THE 'ITER' PLACE.  THIS IS THE VARIABLE ITER IN OUTPUTQ2.
C
      IF(ITER.EQ.0)THEN  
         TRACKFLAG='YES'
         WRITE(2,*)'DROGUE TRACKING STARTED AT T = ',T1Q
         DO I=1,NMND
            DO L=1,NNV
               U1(I,L)=U2(I,L)
               V1(I,L)=V2(I,L)
               W1(I,L)=W2(I,L)
               Z1(I,L)=Z2(I,L)
               Q21(I,L)=Q22(I,L)
               Q2L1(I,L)=Q2L2(I,L)
               RHO1(I,L)=RHO2(I,L)
               TMP1(I,L)=TMP2(I,L)
               ENZM1(I,L)=ENZM2(I,L)
               ENZH1(I,L)=ENZH2(I,L)
               ENZQ1(I,L)=ENZQ2(I,L)
            END DO
         END DO
C
C ***  COMPUTE MODEL TIME-STEP
C         
         DT_TS=T2Q-T1Q
         T1DR=T2Q
         CALL INPUT(NMND,NMEL,NNV,MAXDROG,T1Q,T2Q,DT_TS,DT_DR,
     +              IPRINT,ITRACK,NDR,Z1,Z2,
     +              IDR,JJDR,LLDR,XDR,YDR,ZDR,
     +              DOMAINNAME,INDNAME)
         IF(NDR.EQ.0)RETURN
         NACTIVE=NDR
         WRITE(*,7000)
         WRITE(13,7000)
         RETURN
      END IF 
C 
C *** DETERMINE IF ENOUGH MODEL TIMESTEPS (DT_TS'S) HAVE ELAPSED TO 
C *** COMPLETE A TRACKING FROM T1 TO T2 AT 'ITRACK*DT_TS' TIMESTEP
C
C *** THE PARAMETER ITRACK DETERMINES WHEN TO UPDATE DROGUE POSITIONS,
C *** RELATIVE TO THE MODEL TIME-STEP.  IF DT_TS (THE MODEL TIME-STEP) 
C *** IS 15 MINUTES AND ITRACK=4, THEN THE DROG3DDT UPDATES PARTICLE
C *** POSITIONS EVERY 4 DT_TS'S, OR EVERY HOUR.  DT_DR (THE TRACKING 
C *** TIME-STEP) IS 1 HOUR.  
C *** IF TIMESTEPS ARE NOT EQUAL, WAIT UNTIL 'ITER' IS A 
C *** MULTIPLE OF ITRACK TO TRACK
C
      T2DR=T2Q
      IF(TRACKFLAG.EQ.'NO')RETURN
      IF(MOD(ITER,ITRACK).EQ.0) THEN
         GOTO 15
      ELSE
         RETURN
      END IF        
C
C *** BEGIN TRACKING OF PARTICLES FOR THIS TRACKING TIME-STEP
C
C *** LOOP OVER EACH DROGUE
C
 15   CONTINUE



c
c kws 11/26/2003
c

      NN=NMND
      NE=NMEL
      if (rwtype(1:4).eq.'SMAG')then
       do n=1,nn
        U(n)=UT1(n,1)
        V(n)=VT1(n,1)
       enddo 
       call smagor1(U,V,AR,X,Y,IN,NE,NN,AHIN,SMAGR,AHMIN)
c
c Output smagorinsky diffusion field every 
c
       IF (MOD(ITQ,10).EQ.0) THEN
        CALL INTTOCHR(ITQ,TLABEL,SSLOC)
            FILSMAG = 'SMAGD'//TLABEL(1:SSLOC)//'ahi.s2r'
            JLOC = INDEX(FILSMAG,' ') - 1
            OPEN (UNIT=66, FILE=FILSMAG(1:JLOC), STATUS='UNKNOWN')
           WRITE (66,*) 'MESHNAME'
            WRITE (66,*) 'Smagorinsky Diffusivity # '
     &                    //TLABEL(1:SSLOC)       
            WRITE (66,'(I7,E12.4)') (I,AHIN(I),I=1,NN)
            CLOSE (66)
      ENDIF

      endif
      DO 160 II = 1,NDR
C
C  ***IF THIS ELEMENT WAS ELIMINATED ON A PREVIOUS STEP THEN SKIP
C
         IF (IDR(II).EQ.0) GO TO 160
C
C *** IF IT'S NOT TIME TO RELEASE THIS DROGUE, SKIP TRACKING;
C
         IF (ITER.LT.ITER0(II))GOTO 160
C
C *** X0, Y0,Z0,J0,L0 ARE THE STARTING POSITIONS FOR DROGUE AT THIS TIME-STEP
C
         J0 = JJDR(II)
         L0 = LLDR(II)
         X0 = XDR(II)
         Y0 = YDR(II)
         Z0 = ZDR(II)
C
C *** TRACK DROGUE FROM T1DR TO T2DR
C
         CALL RK4(NNV,T1DR,T2DR,DT_DR,
     +            Z1,Z2,
     +            U1,U2,
     +            V1,V2,
     +            W1,W2,
     +            J0,JNEW,L0,LNEW,
     +            X0,Y0,Z0,U0,V0,W0,
     +            XNEW,YNEW,ZNEW,IBOUN)
C
C *** IF PARTICLE HITS A BOUNDARY, RESTORE PREVIOUS POSITION
C *** AND ELIMINATE IT
C
         IF (IBOUN.EQ.0) THEN
C
C ***       VERTICAL KICK SEGMENT
C         
            IF(VTURBFLAG.EQ.0)GOTO 665
C
C ***       HORIZONTAL KICK SEGMENT
C
 665        IF(HTURBFLAG.EQ.0)GOTO 666
	    BETA_X=GASDEV(ISEED)
	    BETA_Y=GASDEV(ISEED)
	    KICK_X=SQRT(ALPHA*K_X*DT_DR)
	    KICK_Y=SQRT(ALPHA*K_Y*DT_DR)
	    XHOR=XNEW+BETA_X*KICK_X
	    YHOR=YNEW+BETA_Y*KICK_Y
            WRITE(*,*)XNEW,XHOR,YNEW,YHOR
	    CALL BELEL(JNEW,XHOR,YHOR,IND)
	    IF (IND.EQ.0) THEN
	       CALL FIND_ELEMENT(JNEW,JTEMP,XHOR,YHOR,IND)
	       IF (IND.EQ.0) THEN
		  IBOUN=1
	      WRITE(13,*)
     +  	  'HOR-KICK DROGUE EXITED HORZ AT T (DAYS) =',T1Q/86400.
                  GOTO 667
	       END IF
	       JNEW = JTEMP
	    END IF
	    CALL FIND_LEVEL(NNV,T1Q,T2Q,Z1,Z2,
     +  		    JNEW,XHOR,YHOR,ZNEW,T2Q,
     +                      LNEW,LTEMP,INDLL,BOTFLAG)
	    
	    IF(LTEMP.GT.NNV)THEN
	       PRINT*,'DROGUE EXITED FREE SURFACE.'
	       PRINT*,'IDR = ',IDR(II)
	       PRINT*,'ZO = ',ZNEW
	       PRINT*,'LNEW = ',LNEW
	       STOP
	    ELSEIF(LTEMP.LT.1)THEN
	       PRINT*,'2-DROGUE EXITED BOTTOM.'
	       PRINT*,'IDR = ',IDR(II)
	       PRINT*,'ZO = ',ZNEW
	       PRINT*,'LNEW = ',LNEW
	       CALL DEPTH_INTERP(NNV,T1DR,T2DR,Z1,Z2,JNEW,T2DR,ZNOW)
	       CALL GET_DEPTH(ZNOW,JNEW,XHOR,YHOR,LTEMP,ZLOW)	    
	       ZNEW=ZLOW+.1
	    ENDIF
            XNEW=XHOR
            YNEW=YHOR

 666        CONTINUE
C
C ***       UPDATE THE PARTICLE'S NEW POSITION;  
	    JJDR(II) = JNEW
	    LLDR(II) = LNEW
	    XDR(II)  = XNEW
	    YDR(II)  = YNEW
	    ZDR(II)  = ZNEW
C
C *** ELSE,ELIMINATE PARTICLE
C 
         ELSE
C
C *** SET IBOUN TO 0 FOR NEXT!! DROGUE; I.E., DROGUE WITHIN DOMAIN
C *** AND RESTORE STARTING LOCATIONS INTO ARRAYS
C
 667        IBOUN = 0
            JJDR(II) = J0
            LLDR(II) = L0
            XDR(II)  = X0
            YDR(II)  = Y0
            ZDR(II)  = Z0
            IDR(II)  = 0
            NACTIVE=NACTIVE-1  ! THIS IS THE NUMBER STILL ACTIVE!!
            WRITE(13,8010)II,ITER,T1DR/86400.E0
            WRITE(*,8010)II,ITER,T1DR/86400.E0
         END IF

  160 CONTINUE
C
C *** CALL SUBROUTINE OUTPUT
C
 165  IF(MOD(ITER,IPRINT*ITRACK).EQ.0) CALL OUTPUT(NNV,T1DR,T2DR,
     +                                     Z1,Z2,IDR,JJDR,LLDR,
     +                                     NDR,XDR,YDR,ZDR)
C
C *** MOVE T2 ARRAYS TO T1 ARRAYS FOR NEXT TIMESTEP
C
      DO I=1,NMND
         DO L=1,NNV
            U1(I,L)   =U2(I,L)
            V1(I,L)   =V2(I,L)
            W1(I,L)   =W2(I,L)
            Z1(I,L)   =Z2(I,L)
            Q21(I,L)  =Q22(I,L)
            Q2L1(I,L) =Q2L2(I,L)
            RHO1(I,L) =RHO2(I,L)
            TMP1(I,L) =TMP2(I,L)
            ENZM1(I,L)=ENZM2(I,L)
            ENZH1(I,L)=ENZH2(I,L)
            ENZQ1(I,L)=ENZQ2(I,L)
         END DO
      END DO  
      T1DR=T2DR
      IF(NACTIVE.EQ.0)THEN
         WRITE(13,*)'ALL DROGUES ELIMINATED BEFORE TIME ELAPSED.'
         WRITE(*,*)'ALL DROGUES ELIMINATED BEFORE TIME ELAPSED.'
         CLOSE(12)
         CLOSE(13)
CBOB         TRACKFLAG='NO'
CBOB         STOP 'ALL DROGUES ELIMINATED!! QUODDY TERRMINATED!!'
      END IF
C
C *** RETURN TO TIMESTEPPER
C
      RETURN
C
C *** FORMATS
C
 7000 FORMAT('START OF TIME-STEPPING AND DROGUE TRACKING')
 8000 FORMAT(I4,' DROGUE(S) WERE ELIMINATED DURING TRACKING')
 8010 FORMAT('DROGUE # ',I4,' ELIMINATED AT ITER = ',I7,' TIME(DAYS) ',
     +F10.4)
      END
C
C************************************************************************
C
C *** INPUT.F  READS ALL NECESSARY FILES FROM THE CURRENT WORKING 
C *** DIRECTORY 
C
C TIME-DEPENDENT VERSION SPRING 1994
C
C************************************************************************
      SUBROUTINE INPUT(NMND,NMEL,NNV,MAXDROG,T1SEC,T2SEC,DT_TS,DT_DR,
     +                 IPRINT,ITRACK,NDR,DEPTH1,DEPTH2,
     +                 IDR,JJDR,LLDR,XDR,YDR,ZDR,
     +                 DOMAINNAME,INDNAME)
C************************************************************************
C
C *** INCLUDE QUODDY PARAMETER LIST 
C
      INCLUDE 'DCMS.DIM'
C
C *** DECLARE VARIABLES COMING THROUGH PARAMETER LIST
C
C     NMND    - NUMBER OF NODES IN DOMAIN
C     NMEL    - NUMBER OF ELEMENTS IN DOMAIN
C     NNV     - NUMBER OF VERTICAL LEVELS IN DOMAIN
C     MAXDROG - NUMBER OF VERTICAL LEVELS IN DOMAIN
C     T1SEC   - TIME AT START OF TIMESTEP 
C     T2SEC   - TIME AT END OF TIMESTEP 
C     DT_TS   - DELTA T (T2-T1) USED BY MODEL 
C     DT_DR   - DELTA T (T2-T1) USED BY DROG3DDT
C     IPRINT  - OUTPUT TIMESTEP INTERVAL
C     ITRACK  - NUMBER OF MODEL DTS TO SKIP BEFORE TRACKING 
C     NDR     - NUMBER OF DROGUES AFTER INITIAL LOCATION
C     IDR     - ARRAY CONTAINING ACTIVE/INACTIVE LIST
C     DEPTH1  - ARRAY CONTAINING VERTICAL GRID POSITIONS AT TIME T1
C     DEPTH2  - ARRAY CONTAINING VERTICAL GRID POSITIONS AT TIME T2      
C     JJDR    - ARRAY CONTAINING CURRENT DROG ELEMENT 
C     LLDR    - ARRAY CONTAINING CURRENT DROGUE LEVEL
C     XDR     - ARRAY CONTAINING X-POSITIONS OF DROGUES
C     YDR     - ARRAY CONTAINING Y-POSITIONS OF DROGUES
C     ZDR     - ARRAY CONTAINING Z-POSITIONS OF DROGUES
C     DOMAINNAME - NAME OF DOMAIN (GRID) 
C     INDNAME    - NAME OF DROG INPUT FILE; FILDR3 IN QUODDY2
C     TRACKFLAG  - INTERNAL TRACKING FLAG
C
      INTEGER NMND,NMEL,NNV,IPRINT,NDR,ITRACK,MAXDROG
      INTEGER IDR(*),JJDR(*),LLDR(*)
      REAL T1SEC,T2SEC,DT_TS,DT_DR
      REAL DEPTH1(1:NNDIM,1:NNVDIM),DEPTH2(1:NNDIM,1:NNVDIM)
      REAL XDR(*),YDR(*),ZDR(*)
      CHARACTER*(*) DOMAINNAME,INDNAME
C
C *** COORDINATES OF THE NODES FROM .GR2 FILE
C
      COMMON /COORDS/XND,YND
      REAL XND(1:NNDIM),YND(1:NNDIM)
      INTEGER NMNDGR2,NMELGR2
C
C *** TABLE OF ELEMENTS FROM .GR2 FILE
C
      COMMON /ELEM/ELEMS
      INTEGER ELEMS(1:NEDIM,1:3)
C
C *** CONNECTIVITY MATRICES - ELEMENT-ELEMENT & NODE-ELEMENT FROM .GR2 FILE
C
      COMMON /CONEC/ICEE,ICNE
      INTEGER ICEE(1:NEDIM,1:3),ICNE(1:NNDIM,1:12)
C
C *** SCALING FACTORS FOR INIT DROGUE POSITIONS AND GRID DIMENSIONS
C 
      REAL SCDRX,SCDRY,SCDRZ,SCNDX,SCNDY,SCNDZ
C
C *** DROGUE POSITION INFORMATION
C
      REAL XD,YD,ZD,ZBOT
      REAL DEPTHNOW(1:3,1:NNVDIM)
C
C *** VARIABLES FOR FILENAMES AND FILE HANDLING
C
      COMMON TRACKFLAG
      CHARACTER*3 TRACKFLAG
      CHARACTER*72 COMMENT,CASENAME,LINE
      CHARACTER*72 FILEIND,FILEGRID,FILEPATH,FILEDIAG
      CHARACTER*4 GRID,VEL,IND,PATH
      DATA GRID/'.gr2'/,VEL/'.vel'/,IND/'.ind'/,PATH/'.pth'/
      CHARACTER*5 DIAG
      DATA DIAG/'.diag'/
      CHARACTER*1 BLANK
      DATA BLANK/' '/
      INTEGER IOS
C
C *** LOOP COUNTERS
C
      INTEGER I,J,N,II
C
C *** DEFINE PI.
C
      REAL PI
      PI=ACOS(-1.E0)
      WRITE(*,*)NMND,NMEL,NNV,MAXDROG,T1SEC,T2SEC,DT_TS,DT_DR
C 
C *** EXTRACT CASENAME OF DROGUE RUN FROM INDNAME
C
      CASENAME=INDNAME(1:INDEX(INDNAME,BLANK)-5)
C
C *** OPEN .DIAG FILE FIRST; ALL DIAGNOSTICS WILL BE
C *** WRITTEN TO '.DIAG' FILE
C
      FILEDIAG=CASENAME(1:INDEX(CASENAME,BLANK)-1)//DIAG
      OPEN (UNIT=13,FILE=FILEDIAG,IOSTAT=IOS)
      IF(IOS.NE.0)THEN
         WRITE(*,103)FILEDIAG(:INDEX(FILEDIAG,BLANK)-1)
         WRITE(*,8050)
         TRACKFLAG='NO'
         RETURN
      END IF
C 
C *** ECHO DROGUE CASENAME TO .DIAG FILE
C
      WRITE(13,6004)CASENAME
      WRITE(13,6005)DOMAINNAME(1:INDEX(DOMAINNAME,BLANK)-1)
      WRITE(13,*)
C
C *** OPEN DROGUE DATA FILE, SUFFIXED '.IND'
C
      FILEIND=CASENAME(:INDEX(CASENAME,BLANK)-1)//IND
      OPEN (UNIT=11,FILE=FILEIND,IOSTAT=IOS,STATUS='OLD')
      IF(IOS.NE.0)THEN
         WRITE(13,103)FILEIND(:INDEX(FILEIND,BLANK)-1)
         WRITE(13,8050)
         WRITE(*,103)FILEIND(:INDEX(FILEIND,BLANK)-1)
         WRITE(*,8050)
         TRACKFLAG='NO'
         CLOSE(13)
         RETURN
      END IF   
C
C *** OPEN FILE FOR OUTPUT, WITH SAME NAME AS DROGUE CASE NAME ;
C *** THE '.IND' FILE WILL BE ECHOED INTO THE FIRST LINES OF THE 
C *** OUTPUT FILE
C
      FILEPATH=CASENAME(1:INDEX(CASENAME,BLANK)-1)//PATH
      OPEN (UNIT=12,FILE=FILEPATH,IOSTAT=IOS)
      IF(IOS.NE.0)THEN
         WRITE(13,103)FILEPATH(:INDEX(FILEPATH,BLANK)-1)
         WRITE(13,8050)
         WRITE(*,103)FILEPATH(:INDEX(FILEPATH,BLANK)-1)
         WRITE(*,8050)
         TRACKFLAG='NO'
         CLOSE(13)
         RETURN
      END IF
C
C *** ECHO COMMON_BLOCK INCLUDE FILE TO .DIAG FILE
C *** THE FOLLOWING PARAMETERS ARE SET IN THE DCMS.DIM FILE 
C *** UPON COMPILATION OF THE CODE.  THE CURRENT VALUES ARE 
C *** WRITTEN TO THE .DIAG FILE BELOW.
C
C *** NNDIM  - MAXIMUM NUMBER OF HORIZONTAL NODES
C *** NEDIM  -  "        "    "     "      ELEMENTS
C *** NNVDIM -  "        "    " VERTICAL LEVELS
C *** MAXDROG -  "        "    " DROGUES
C *** DTOL  - DIMENSIONLESS TOLERANCE VALUE USED TO FIND ELEMENT OF 
C ***          CURRENT DROGUE 
C
      WRITE(13,*)' ' 
      WRITE(13,9000)'DCMS.DIM PARAMETER VALUES:'
      WRITE(13,6007)'MAX NUMBER OF NODES             = ',NNDIM
      WRITE(13,6007)'MAX NUMBER OF ELEMENTS          = ',NEDIM
      WRITE(13,6007)'MAX NUMBER OF VERT LEVELS       = ',NNVDIM
      WRITE(13,9000)'DROG3DDT PARAMETER VALUES:'
      WRITE(13,6007)'MAX NUMBER OF STARTING DROUGES  = ',MAXDROG
C
C *** WRITE GRIDNAME AT TOP OF '.PTH' FILE
C
      WRITE(12,9000) DOMAINNAME
      WRITE(12,*)'RECORD ABOVE THIS LINE IS THE DOMAIN NAME ON WHICH'
      WRITE(12,*)'VELOCITIES WERE COMPUTED.'      
C
C READ RUN PARAMETERS FROM '.IND' FILE
C
      READ(11,9020) COMMENT
      READ(11,9020) COMMENT
      READ(11,9020) COMMENT
      READ(11,*) ITRACK
      READ(11,9020) COMMENT
      READ(11,*) IPRINT
C *** MAKE SURE ITRACK AND IPRINT ARE INTEGER 1 OR >.
      IF(IPRINT.LT.1)THEN
         WRITE(13,5010)
         WRITE(13,8050)
         WRITE(*,8050)
         TRACKFLAG='NO'
         CLOSE(13)
         RETURN
      ELSEIF(ITRACK.LT.1)THEN 
         WRITE(13,5015)
         WRITE(13,8050)
         WRITE(*,8050)
         TRACKFLAG='NO'
         RETURN
      END IF
      READ(11,9020) COMMENT
      READ(11,*) SCNDX,SCNDY,SCNDZ
      FILEGRID=DOMAINNAME(:INDEX(DOMAINNAME,BLANK)-1)//GRID
      OPEN (UNIT=9,FILE=FILEGRID,IOSTAT=IOS,STATUS='OLD')
      IF(IOS.NE.0)THEN
         WRITE(13,103)FILEGRID(:INDEX(FILEGRID,BLANK)-1)
         WRITE(13,8050)
         WRITE(*,103)FILEGRID(:INDEX(FILEGRID,BLANK)-1)
         WRITE(*,8050)
         TRACKFLAG='NO'
         CLOSE(13)
         RETURN
      END IF
C
C *** READ GRID DATA FROM FILE '.GR2'
C
C *** CHECK FOR CONSISTANCY BETWEEN NMND AND NMEL AS SUPPLIED FROM 
C *** TIMESTEPPER VERSUS NMNDGR2 AND NMELGR2 FROM THE TOP OF THE .GR2 FILE
C
      READ(9,*) NMELGR2,NMNDGR2
      IF(NMEL.NE.NMELGR2.AND.NMND.NE.NMNDGR2)THEN
         WRITE(13,5000)
         WRITE(13,8050)
         WRITE(*,5000)
         WRITE(*,8050)
         TRACKFLAG='NO'
         CLOSE(13)
         RETURN
      ELSEIF (NMEL.GT.NEDIM) THEN
         WRITE(13,6000)'ELEMENTS',NEDIM
         WRITE(13,8050)
         WRITE(*,6000)'ELEMENTS',NEDIM
         WRITE(*,8050)
         TRACKFLAG='NO'
         CLOSE(13)
         RETURN
      ELSEIF (NMND.GT.NNDIM) THEN
         WRITE(13,6000)'NODES',NNDIM
         WRITE(13,8050)
         WRITE(*,6000)'NODES',NNDIM
         WRITE(*,8050)
         TRACKFLAG='NO'
         CLOSE(13)
         RETURN
      END IF
C
C *** READ COORDINATES OF NODES
C
      DO 10 I = 1,NMND
          READ(9,*) N,XND(N),YND(N)
          XND(N)=XND(N)*SCNDX
          YND(N)=YND(N)*SCNDY
   10 CONTINUE
C
C *** READ THE TABLE OF ELEMENTS
C
      DO 20 I = 1,NMEL
         READ(9,*) N,ELEMS(N,1),ELEMS(N,2),ELEMS(N,3)
   20 CONTINUE
C
C *** ELEMENT-ELEMENT ADJACENCY INFORMATION
C
      DO 30 I = 1,NMEL
         READ(9,*) N,ICEE(N,1),ICEE(N,2),ICEE(N,3)
   30 CONTINUE
C
C *** NODE-ELEMENT ADJACENCY INFORMATION
C
      DO 40 I = 1,NMND
         READ(9,*) N,ICNE(N,1),(ICNE(N,J+1),J=1,ICNE(N,1))
   40 CONTINUE

      WRITE(13,8070) 
C 
C *** COMPUTE DROGUE TRACKING TIME-STEP
C
      DT_DR=ITRACK*DT_TS
C
C *** REPORT RUN CHARACTERISTICS TO .DIAG FILE
C
      WRITE(13,*) ' '
      WRITE(13,9000)'TOTAL RUN PARAMETERS'
      WRITE(13,7015)T2SEC/3600.E0
      WRITE(13,7025)DT_TS
      WRITE(13,7026)DT_DR
      WRITE(13,7027)DT_DR*IPRINT
C
C *** SCALING FACTORS FOR DROG COORDINATES
C
      READ(11,9020) COMMENT
      READ(11,*) SCDRX,SCDRY,SCDRZ
C
C ***     TOTAL NUMBER OF DROGUES AT START
C
      READ(11,9020) COMMENT
      READ(11,*) NDR
      IF (NDR.GT.MAXDROG) THEN
          WRITE(13,6000)'DROGS.',MAXDROG
          WRITE(13,8050)
          WRITE(*,6000)'DROGS.',MAXDROG
          WRITE(*,8050)
          TRACKFLAG='NO'
          CLOSE(13)
          RETURN
      END IF
C
C *** READ IN DROGUE COORDINATES
C
      READ(11,9020) COMMENT
      DO 90 I = 1,NDR
          READ(11,*) XD,YD,ZD
          XDR(I)=XD*SCDRX
          YDR(I)=YD*SCDRY
          ZDR(I)=ZD*SCDRZ
  90  CONTINUE
C
      WRITE(13,7070)NNV
      WRITE(*,7070)NNV
C
C *** COMPUTE AREA COORDINATES FOR ELEMENT INTERPOLATION FUNCTIONS
C
      CALL COMP_AREAS(NMEL)
C
C *** INITIALIZE DATA FOR BELEL
C
      CALL BELELINIT(NMEL)
C
C *** CALL LOCDROGINIT 
C
      WRITE(13,7075)NDR
      WRITE(*,*)NMND,NMEL,NNV,T1SEC,T2SEC
      CALL LOCDROGINIT(NMND,NMEL,NNV,T1SEC,T2SEC,DEPTH1,DEPTH2,
     +                 SCDRX,SCDRY,SCDRZ,NDR,XDR,YDR,ZDR,
     +                 IDR,JJDR,LLDR)
C
C *** CLOSE '.GR2' FILE, AND REWIND, ECHO, AND CLOSE '.IND' FILE 
C
      CLOSE(UNIT=9)
      WRITE(12,*)'************** BEGIN .IND FILE ECHO *************'
      REWIND(UNIT=11)
 55   READ(11,6006,END=56)LINE
         WRITE(12,6006)LINE
         GOTO 55
 56   CONTINUE 
      CLOSE(UNIT=11)
      WRITE(12,*)'************** END .IND FILE ECHO ***************' 
C
C *** OUTPUT OTHER INFO TO .PTH FILE
C
      WRITE(12,7025)DT_TS
      WRITE(12,7026)DT_DR      
      WRITE(12,7027)DT_DR *IPRINT     
C
C *** WRITE INITIAL PATH INFO TO .PTH FILE, THE DROGUE POSITIONS AT
C *** ITERATION 0.  
C
      WRITE(12,9000) 'XXXX'
      WRITE(12,9002) '0000',DT_DR*IPRINT,NDR
      IF(NDR.EQ.0)THEN
         WRITE(12,8040)
         WRITE(12,8050)
         WRITE(*,8040)
         TRACKFLAG='NO'
      END IF
      DO II=1,NDR
         CALL DEPTH_INTERP(NNV,T1SEC,T2SEC,DEPTH1,DEPTH2,
     +                  JJDR(II),T1SEC,DEPTHNOW)
         WRITE(12,9005) XDR(II),YDR(II),ZDR(II),T2SEC
      END DO
      WRITE(13,8090)
      WRITE(13,*) 
       
      RETURN
C       10        20        30        40        50        60        70     
C23456789.123456789.123456789.123456789.123456789.123456789.123456789.12
 103  FORMAT('ERROR FINDING/OPENING ',A,' FILE; NO TRACKING')
 
 5000 FORMAT('GRID DIMENSIONS SUPPLIED BY TIMESTEPPER DO NOT MATCH ',
     +'THOSE FOUND IN THE .GR2 FILE')
 5010 FORMAT('IPRINT IN .IND FILE MUST BE >= 1')
 5015 FORMAT('ITRACK IN .IND FILE MUST BE >= 1')
 6000 FORMAT('TOO MANY',A13,'!!  MAXIMUM =',I6,'   CHECK PARAMETER ',
     +'STATEMENTS')
 6004 FORMAT('CASENAME   = ',A30)
 6005 FORMAT('DOMAINNAME = ',A30)
 6006 FORMAT(A72)
 6007 FORMAT(A34,I10)
 6070 FORMAT('NUMBER OF VERTICAL NODES IS           : ',I5)
 7010 FORMAT('LENGTH OF TRACKING (HOURS)            : ',F10.3)
 7015 FORMAT('TIME AT START OF TRACKING (HOURS)     : ',F10.3)
 7019 FORMAT('NUMBER OF MODEL INTERVALS             : ',I10)
 7020 FORMAT('NUMBER OF TRACKING INTERVALS          : ',I10)
 7025 FORMAT('MODEL TIME-STEP (SECONDS)             : ',F14.4)
 7026 FORMAT('DROG3DDT  TIME-STEP (SECONDS)         : ',F14.4)
 7027 FORMAT('DROG3DDT OUTPUT INTERVAL (SECONDS)    : ',F14.4)
 7070 FORMAT('NUMBER OF VERTICAL NODES IS           : ',I5)
 7075 FORMAT('NUMBER OF DROGUES READ FROM .IND FILE : ',I5)
 7080 FORMAT('WARNING:  DTMIN IS 0.0; MAXSTP SET TO 1000')
 7095 FORMAT('COMPONENT ',I4,' IN .IND LIST NOT USED')
 8040 FORMAT('NO DROGUES INITIALLY LOCATED WITHIN DOMAIN.')
 8050 FORMAT('TRACKING FLAG SET TO NO IN DROG3DDT INPUT ROUTINE.')
 8070 FORMAT('GRID-INPUT COMPLETE')
 8090 FORMAT('INPUT ROUTINE COMPLETE')
 9000 FORMAT(A)
 9002 FORMAT(A4,1X,F20.6,1X,I10)
 9005 FORMAT (5(2X,F12.4),I4)
 9020 FORMAT(A72)
      END
C
C***********************************************************************
      SUBROUTINE LOCDROGINIT(NMND,NMEL,NNV,T1,T2,DEPTH1,DEPTH2,
     +                       SCDRX,SCDRY,SCDRZ,
     +                       NDR,XDR,YDR,ZDR,IDR,JJDR,LLDR)
C***********************************************************************
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
C *** DROGUE POSITION INFORMATION
C
      REAL XDR(*),YDR(*),ZDR(*)
      REAL XSTART,YSTART,ZSTART
      INTEGER IDR(*),JJDR(*),LLDR(*),NDR,LL
      LOGICAL LLIND,BOTFLAG
      INTEGER IND,LSTART,JJ,NOTFND
      REAL SCDRX,SCDRY,SCDRZ,T1,T2
      REAL DEPTH1(1:NNDIM,1:NNVDIM),DEPTH2(1:NNDIM,1:NNVDIM)
C
C *** NUMBER OF NODES, NUMBER OF ELEMENTS, NUMBER OF VERTICAL LEVELS
C
      INTEGER NMND,NMEL,NNV
      INTEGER III,I
C
C *** FIND_LEVEL NOW NEEDS A STARTING SEED, LIKE FIND_ELEMENT;
C *** SO ON THE INITIAL LOCATION CALLS, USE 
C ***                 LSTART =NINT(FLOAT(NNV)/2.D0).
C 
      LSTART = NINT(FLOAT(NNV)/2.D0)
C
C *** FIND STARTING ELEMENT FOR EACH DROGUE
C
C
C *** LOCATE DROGUE COORDINATES IN HORIZONTAL (ELEMENT JJ) AND
C *** VERTICAL (LEVEL IMMEDIATELY BELOW DROGUE POSITION, LLDR)
C
      NOTFND=0
      DO 130 III = 1,NDR
          XSTART = XDR(III)
          YSTART = YDR(III)
          ZSTART = ZDR(III)
          IF (NOTFND.GT.0) THEN
              XDR(III-NOTFND)=XDR(III)
              YDR(III-NOTFND)=YDR(III)
              ZDR(III-NOTFND)=ZDR(III)
          ENDIF
          DO 120 I = 1,NMEL
              CALL BELEL(I,XSTART,YSTART,IND)
              IF (IND.EQ.1) THEN
                  JJDR(III-NOTFND) = I
                  GO TO 125
              END IF
  120     CONTINUE
          WRITE(13,8010) III
          NOTFND=NOTFND+1
	  IF ((NDR-NOTFND).EQ.0)GOTO 140
          GOTO 130

  125     JJ=JJDR(III-NOTFND)
C
          CALL FIND_LEVEL(NNV,T1,T2,DEPTH1,DEPTH2,
     +                    JJ,XSTART,YSTART,ZSTART,T1,
     +                    LSTART,LL,LLIND,BOTFLAG)
	  IF (.NOT.LLIND) THEN 
	     WRITE(13,8020) III
             NOTFND=NOTFND+1
	     IF ((NDR-NOTFND).EQ.0)GOTO 140
	  ELSEIF (BOTFLAG) THEN 
	     WRITE(13,8020) III
             NOTFND=NOTFND+1
	     IF ((NDR-NOTFND).EQ.0)GOTO 140
         ELSE
            LLDR(III-NOTFND)=LL
         END IF 
  130 CONTINUE
C
C *** INITIALIZE DROGUE ACTIVE/INACTIVE ARRAY
C
      DO  III = 1,NDR
         IDR(III) = 1    
      END DO
       
      NDR=NDR-NOTFND
      WRITE(13,8030) NDR
      WRITE(13,8035)
      RETURN
 140  CONTINUE      
      NDR=0
      RETURN
 8010 FORMAT('COULD NOT FIND STARTING ELEMENT FOR DROGUE ',I6)
 8020 FORMAT('DROGUE # ',I6,' VERTICALLY OUT OF BOUNDS.')
 8030 FORMAT('NUMBER OF DROGUES INITIALLY WITHIN DOMAIN = ',I6)
 8035 FORMAT('LOCDROGINIT COMPLETED')
 9009 FORMAT(6(F14.4,1X))
      END
C
C***********************************************************************
      SUBROUTINE BELELINIT(NMEL)
C***********************************************************************
C
C *** INITIALIZE DATA FOR SUBROUTINE BELEL
C *** PTURNER 8-9-88
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
      COMMON /ARJ/AR
      COMMON /ABA0/A,B,A0
      COMMON /ARBEL/T
      INTEGER NMEL
      REAL*8 T(1:NEDIM,1:3)
      REAL*8 AR(1:NEDIM)
      REAL*8 A(1:NEDIM,1:3)
      REAL*8 B(1:NEDIM,1:3)
      REAL*8 A0(1:NEDIM,1:2)
      INTEGER J
C
C *** THE SHAPE FUNCTION, S, FOR ELEMENT "J" IS USED TO DETERMINE WHETHE
C *** POINT (XP,YP) IS WITHIN OR ON THE ELEMENT'S BOUNDARY.  IF THE POIN
C *** IS WITHIN OR ON THE ELEMENT, ALL 3 SHAPE FUNCTIONS WILL RANGE IN
C *** VALUE BETWEEN [0,1].  IF (XP,YP) LIES OUTSIDE OF THE ELEMENT, ONE
C *** OR MORE OF THE ELEMENT'S SHAPE FUNCIONS WILL ATTAIN A VALUE EITHER
C *** LESS THAN 0 OR GREATER THAN 1 ( S > 1, OR S < 0 ).
C
C *** COMPUTE EACH OF THE 3 SHAPE FUNCTIONS [S(1), S(2), S(3)] BASED ON
C *** THE COORDINATES (XP,YP).
C
C
      DO 90 J=1,NMEL
         T(J,1) = A0(J,1)*2.D0
         T(J,2) = A0(J,2)*2.D0
         T(J,3) = (2.D0*AR(J)-T(J,1)-T(J,2))
 90   CONTINUE   
      WRITE(13,6000) 
      RETURN
 6000 FORMAT('BELELINIT COMPLETED')
      END

C
C***********************************************************************
      SUBROUTINE BELEL(J,XP,YP,NFLAG)
C***********************************************************************
C
C *** DETERMINE WHETHER THE POINT (XP,YP) IS WITHIN ELEMENT "J"
C *** (OR ON ITS BOUNDARIES).  IF SO, IND=1; IF NOT, IND=0
C
C *** THIS IS THE CROSS-PRODUCT METHOD OF CHRISTOPHER NAIMIE,
C *** NUMERICAL METHODS LABORATORY, DARTMOUTH COLLEGE, HANOVER NH
C
C PURPOSE: THIS SUBROUTINE DETERMINES IF A POINT (XVAL,YVAL) IS ON A 
C            LOCAL 2-D TRIANGULAR LINEAR ELEMENT
C    
C RESTRICTIONS: APPLICABLE ONLY FOR  2-D TRIANGULAR LINEAR ELEMENT
C
C INPUTS:  J         - CURRENT ELEMENT BEING TESTED
C          XVAL,YVAL - LOCATION OF THE POINT 
C
C OUTPUTS: NFLAG IS SET EQUAL TO 1 IF THE POINT (XVAL,YVAL) IS IN/ON 
C            THE LOCAL ELEMENT J.
C          NFLAG IS SET EQUAL TO 0 IF THE POINT (XVAL,YVAL) IS NOT ON 
C            THE LOCAL ELEMENT J.
C
C HISTORY:  WRITTEN BY CHRISTOPHER E. NAIMIE
C           DARTMOUTH COLLEGE
C           22 APRIL 1992
C 
C           IMPLEMENTED IN DROG3DDT BY BRIAN BLANTON
C           UNC-CH
C           10 Aug 1995
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
      COMMON /COORDS/XND,YND
      REAL XND(1:NNDIM),YND(1:NNDIM)
      COMMON /ELEM/ELEMS
      INTEGER ELEMS(1:NEDIM,1:3)
C
C *** INCOMING ARGUMENTS
C
      REAL XP,YP
      INTEGER J,NFLAG
C
C *** LOCAL DECLARATIONS
C
      INTEGER I,K
      REAL XLOCAL(3),YLOCAL(3),CROSSPROD
      REAL VX(3),VY(3),DELX,DELY,D,THETA
C
C *** EXTRACT LOCAL NODE COORDINATES
C
      XLOCAL(1)=XND(ELEMS(J,1))
      XLOCAL(2)=XND(ELEMS(J,2))
      XLOCAL(3)=XND(ELEMS(J,3))
      YLOCAL(1)=YND(ELEMS(J,1))
      YLOCAL(2)=YND(ELEMS(J,2))
      YLOCAL(3)=YND(ELEMS(J,3))
C
C  CALCULATE THE X AND Y COMPONENTS OF VECTORS POINTING FROM (XP,YP)
C  TO EACH NODE ON THE ELEMENT
C
      DO 10 I=1,3
         DELX=XLOCAL(I)-XP
         DELY=YLOCAL(I)-YP
         D=SQRT(DELX**2.0+DELY**2.0)
         THETA=ATAN2(DELY,DELX)
         VX(I)=D*COS(THETA)
         VY(I)=D*SIN(THETA)
 10   CONTINUE

C
C  DETERMINE IF THE POINT IS ON THE ELEMENT BY CALCULATING THE
C    CROSSPRODUCTS OF NEIGHBORING VECTORS IN A DIRECTION
C    WHICH WILL YIELD ALL NON-NEGATIVE NUMBERS IF X,Y IS ON THE ELEMENT.
C    (IE: V1XV2>OR=0, V2XV3>OR=0, AND V3XV1>OR=0 => X,Y IS ON ELEMENT)
      DO 20 I=1,3
         K=I+1
         IF(I.EQ.3)K=1
         CROSSPROD=VX(I)*VY(K)-VY(I)*VX(K)
         IF(CROSSPROD.LT.0.0)THEN
            NFLAG=0
            GO TO 21
         ENDIF
 20   CONTINUE
      NFLAG=1
 21   CONTINUE
C
C END OF ROUTINE
      RETURN
      END
C
C***********************************************************************
      SUBROUTINE BELELOLD(J,XP,YP,IND)
C***********************************************************************
C
C *** DETERMINE WHETHER THE POINT (XP,YP) IS WITHIN ELEMENT "J"
C *** (OR ON ITS BOUNDARIES).  IF SO, IND=1; IF NOT, IND=0
C
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
      REAL*8 TU,TL,DTOL
      PARAMETER(DTOL=0.001D0)
      PARAMETER (TU=1.0D0+DTOL,TL=-DTOL)
      COMMON /ARJ/AR
      COMMON /ABA0/A,B,A0
      COMMON /ARBEL/T
      REAL*8 T(1:NEDIM,1:3),AR(1:NEDIM),S1,S2,S3
      REAL XP,YP
      REAL*8 A(1:NEDIM,1:3)
      REAL*8 B(1:NEDIM,1:3)
      REAL*8 A0(1:NEDIM,1:2)
      INTEGER J,IND
C
C *** THE SHAPE FUNCTION, S, FOR ELEMENT "J" IS USED TO DETERMINE WHETHE
C *** POINT (XP,YP) IS WITHIN OR ON THE ELEMENT'S BOUNDARY.  IF THE POINT
C *** IS WITHIN OR ON THE ELEMENT, ALL 3 SHAPE FUNCTIONS WILL RANGE IN
C *** VALUE BETWEEN [0,1].  IF (XP,YP) LIES OUTSIDE OF THE ELEMENT, ONE
C *** OR MORE OF THE ELEMENT'S SHAPE FUNCIONS WILL ATTAIN A VALUE EITHER
C *** LESS THAN 0 OR GREATER THAN 1 ( S > 1, OR S < 0 ).
C
C *** COMPUTE EACH OF THE 3 SHAPE FUNCTIONS [S1, S2, S3] BASED ON
C *** THE COORDINATES (XP,YP).
C
      S1 = (T(J,1)+B(J,1)*DBLE(XP)+A(J,1)*DBLE(YP))*0.5D0/AR(J)
      IF (S1.GT.TU.OR.S1.LT.TL) GOTO 10 
      S2 = (T(J,2)+B(J,2)*DBLE(XP)+A(J,2)*DBLE(YP))*0.5D0/AR(J)
      IF (S2.GT.TU.OR.S2.LT.TL) GOTO 10
      S3 = (T(J,3)+B(J,3)*DBLE(XP)+A(J,3)*DBLE(YP))*0.5D0/AR(J)
      IF (S3.GT.TU.OR.S3.LT.TL) GOTO 10
C
      IND = 1
      RETURN
 
 10   IND = 0
      RETURN

      END
C
C***********************************************************************
      SUBROUTINE COMP_AREAS(NMEL)
C***********************************************************************
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
      COMMON /COORDS/XND,YND
      REAL XND(1:NNDIM),YND(1:NNDIM)
      COMMON /ARJ/AR
      REAL*8 AR(1:NEDIM)
      COMMON /ABA0/A,B,A0
      REAL*8 A(1:NEDIM,1:3)
      REAL*8 B(1:NEDIM,1:3)
      REAL*8 A0(1:NEDIM,1:2)
      COMMON /ELEM/ELEMS
      INTEGER ELEMS(1:NEDIM,1:3)
      INTEGER NMEL
      COMMON TRACKFLAG
      CHARACTER*3 TRACKFLAG
C
C *** DEFINE LOCAL VARIABLES
C
      INTEGER J,N1,N2,N3
C
C *** COMPUTE AREA COORDINATES FOR ELEMENT INTERPOLATION FUNCTIONS
C
        DO 100 J = 1,NMEL
          N1 = ELEMS(J,1)
          N2 = ELEMS(J,2)
          N3 = ELEMS(J,3)
          A(J,1) = DBLE(XND(N3)) - DBLE(XND(N2))
          A(J,2) = DBLE(XND(N1)) - DBLE(XND(N3))
          A(J,3) = DBLE(XND(N2)) - DBLE(XND(N1))
          B(J,1) = DBLE(YND(N2)) - DBLE(YND(N3))
          B(J,2) = DBLE(YND(N3)) - DBLE(YND(N1))
          B(J,3) = DBLE(YND(N1)) - DBLE(YND(N2))
          A0(J,1) = 0.5D0 * (DBLE(XND(N2))*DBLE(YND(N3))
     +                      -DBLE(XND(N3))*DBLE(YND(N2)))
          A0(J,2) = 0.5D0 * (DBLE(XND(N3))*DBLE(YND(N1))
     +                      -DBLE(XND(N1))*DBLE(YND(N3)))
          AR(J) = 0.5D0 * (A(J,2)*B(J,1)-A(J,1)*B(J,2))
C
C *** CHECK FOR NON-POSITIVE ELEMENT AREAS
C
          IF (AR(J).LE.0.E0) THEN
            WRITE(13,8000)J
            WRITE(13,8050)
            WRITE(*,8050)
            TRACKFLAG='NO'
            RETURN
          END IF
  100   CONTINUE
      WRITE(13,8010)
 8000 FORMAT('NON-POSITIVE AREA AT ELEMENT ',I6,' *** DROG3DDT HALTED')
 8010 FORMAT('ELEMENT AREAS COMPUTED')
 8050 FORMAT('TRACKING FLAG SET TO NO IN DROG3DDT INPUT ROUTINE.')
      END
C
C***********************************************************************
      SUBROUTINE FIND_ELEMENT(J,JJ,XX,YY,INDD)
C***********************************************************************
C
C *** FINDS NEW ELEMENT "JJ" WHICH CONTAINS THE POINT (XX,YY).  "J" IS T
C *** OLD ELEMENT.  INDD=1 IF THE NEW ELEMENT IS FOUND; INDD=0 IF IT
C *** CANNOT BE FOUND.
C *** HORIZONTAL EVALUATION ONLY !!
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
      INTEGER ICEE(1:NEDIM,1:3),ICNE(1:NNDIM,1:12),ELEMS(1:NEDIM,1:3)
      INTEGER I,J,IDONE,N,IND,INDD,JJ,L1,L2,K,KK,NL,ICOUNT,ISUM
      REAL XX,YY
      COMMON /CONEC/ICEE,ICNE
      COMMON /ELEM/ELEMS
      INTEGER NEC(1:NEDIM),NECN(1:NEDIM),NDONE(1:NEDIM)

      IDONE = 1
      NDONE(IDONE) = J
C
C *** USE THE LOGIC DESCRIBED IN SUBROUTINE "LOCATE"
C *** NEXT, CHECK NEIGHBORING ELEMENTS
C
      DO 10 I = 1,3
          N = ICEE(J,I)
          IF (N.EQ.0) GO TO 10
          CALL BELEL(N,XX,YY,IND)
          IF (IND.EQ.1) THEN
              JJ = N
              INDD = 1
              GO TO 130

          END IF

          IDONE = IDONE + 1
          NDONE(IDONE) = N
   10 CONTINUE
      ISUM = 0
      DO 40 L1 = 1,3
          L2 = ELEMS(J,L1)
          DO 30 KK = 1,ICNE(L2,1)
              N = ICNE(L2,KK+1)
              DO 20 I = 1,IDONE
                  IF (N.EQ.NDONE(I)) GO TO 30
   20         CONTINUE
              ISUM = ISUM + 1
              IDONE = IDONE + 1
              NDONE(IDONE) = N
              NEC(ISUM) = N
   30     CONTINUE
   40 CONTINUE
C
C *** PROGRESSIVELY LOOP OVER MORE AND MORE ELEMENTS IN SEARCH OF THE EL
C *** MENT CONTAINING (XX,YY)
C
      ICOUNT = ISUM
      ISUM = 0
      IF (ICOUNT.EQ.0) THEN
          DO 50 I = 2,IDONE
              ICOUNT = ICOUNT + 1
              NEC(ICOUNT) = NDONE(I)
   50     CONTINUE
      END IF

   60 CONTINUE
      DO 70 K = 1,ICOUNT
          N = NEC(K)
          CALL BELEL(N,XX,YY,IND)
          IF (IND.EQ.1) THEN
              JJ = N
              INDD = 1
              GO TO 130

          END IF

   70 CONTINUE
      IF (IDONE.GT.70) THEN
          INDD = 0
          GO TO 130

      END IF
C
C *** ACCUMULATE LIST OF NEW ELEMENTS TO BE CHECKED ON NEXT LOOP
C
      DO 110 K = 1,ICOUNT
          N = NEC(K)
          DO 100 L1 = 1,3
              L2 = ELEMS(N,L1)
              DO 90 KK = 1,ICNE(L2,1)
                  NL = ICNE(L2,KK+1)
                  DO 80 I = 1,IDONE
                      IF (NL.EQ.NDONE(I)) GO TO 90
   80             CONTINUE
                  ISUM = ISUM + 1
                  IDONE = IDONE + 1
                  NDONE(IDONE) = NL
                  NECN(ISUM) = NL
   90         CONTINUE
  100     CONTINUE
  110 CONTINUE
      ICOUNT = ISUM
      IF (ICOUNT.EQ.0) THEN
          INDD = 0
          GO TO 130

      END IF

      DO 120 I = 1,ICOUNT
          NEC(I) = NECN(I)
  120 CONTINUE
      ISUM = 0
      GO TO 60

  130 CONTINUE
      RETURN

      END
C
C*************************************************************************
      SUBROUTINE FIND_LEVEL(NNV,T1,T2,DEPTH1,DEPTH2,
     +                      J,XDR,YDR,ZDR,TNOW,
     +                      LIN,LOUT,INDLL,BOTFLAG)
C*************************************************************************
C
C *** FINDS THE LEVEL IMMEDIATELY BELOW THE DROG WITH HORIZONTAL POSITION
C *** X,Y AND CURRENT DEPTH Z, AND IN ELEMENT J. IND INDICATES IF A
C *** LOWER LEVEL BOUNDING THE DROG AT (X,Y,Z) WAS FOUND. IF SO, LOUT 
C *** CONTAINS THE INTEGER VAULE OF THE LOWER LEVEL.
C
C------------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
      REAL T1,T2,TNOW
      INTEGER NNV
      COMMON /COORDS/XND,YND
      REAL XND(1:NNDIM),YND(1:NNDIM)
      REAL DEPTH1(1:NNDIM,1:NNVDIM),DEPTH2(1:NNDIM,1:NNVDIM)
      REAL DEPTHNOW(1:3,1:NNVDIM)
C
C *** LOCAL VARIABLE DECLARATION
C
      REAL XDR,YDR,ZDR
      REAL ZIL,ZPUP,ZPUPP1,ZPDN,ZOLD
      INTEGER J,LIN,LOUT,LPUP,LPDN,IL
      LOGICAL INDLL,BOTFLAG
      INDLL=.TRUE.
C
      CALL DEPTH_INTERP(NNV,T1,T2,DEPTH1,DEPTH2,J,TNOW,DEPTHNOW)
C
C *** CALLS TO SUBROUTINE GET_DEPTH FIND THE DEPTH H BELOW POINT 
C *** XDR,YDR AT LEVEL L.
C
C *** ASSUME THAT MORE OFTEN THAN NOT, A DROGUE'S LEVEL WILL CHANGE 
C *** BY NO MORE THAN ONE.  THEREFORE, BEFORE LOOPING OVER ALL LEVELS, 
C *** SEARCH LIN+1 AND LIN-1 FIRST.
C
C *** DETERMINATION OF LEVEL UNDER DROG @ XDR,YDR
C
      ZIL=0.
      LPUP=LIN+1
      LPDN=LIN-1
      CALL GET_DEPTH(DEPTHNOW,J,XDR,YDR,LIN,ZOLD)
      CALL GET_DEPTH(DEPTHNOW,J,XDR,YDR,LPUP,ZPUP)
      IF(ZDR.GE.ZOLD.AND.ZDR.LT.ZPUP) THEN 
         LOUT=LIN
         GOTO 11       ! LEVEL REMAINS THE SAME
      END IF
      IF(LPDN.GE.1)THEN
         CALL GET_DEPTH(DEPTHNOW,J,XDR,YDR,LPDN,ZPDN)
         IF(ZDR.GE.ZPDN.AND.ZDR.LT.ZOLD) THEN
            LOUT=LPDN 
            GOTO 11       ! DROGUE DROPS ONE LEVEL, IF POSSIBLE
         END IF
      ENDIF
      IF(LPUP.LT.NNV)THEN
         CALL GET_DEPTH(DEPTHNOW,J,XDR,YDR,LPUP+1,ZPUPP1)
         IF(ZDR.GE.ZPUP.AND.ZDR.LT.ZPUPP1) THEN
            LOUT=LPUP 
            GOTO 11       ! DROGUE RISES ONE LEVEL, IF POSSIBLE
         END IF
      ENDIF
        
C
C *** ELSE, DETERMINE DROGUES LOWER LEVEL BY RAMBO APPROACH
C *** DETERMINATION OF LEVEL UNDER DROG @ XDR,YDR
C
  	DO 10, IL = 1 , NNV
      	CALL GET_DEPTH(DEPTHNOW,J,XDR,YDR,IL,ZIL)
           IF(ZIL.GT.ZDR) THEN 
              LOUT = IL-1
              GOTO 11
           END IF
 10 	CONTINUE
C
C *** IF NO LEVEL BETWEEN 1 AND NNV IS GREATER THAN ZDR, THEN THE DROG
C *** HAS EXITED THRU THE FREE-SURFACE. LLEV IS SET TO NNV TO FLAG
C *** APPROPRIATE DIAGNOSTIC
C
	LOUT = NNV
C
C *** ENSURE LLEV IS WITHIN DOMAIN
C
 11   IF(LOUT.LT.NNV .AND. LOUT.GT.0) THEN 
   	 BOTFLAG = .FALSE.
         INDLL = .TRUE.
      ELSE
	 IF(LOUT.LE.0) THEN
	    BOTFLAG = .TRUE.
	    INDLL = .TRUE.
	 ELSE
	    BOTFLAG = .FALSE.
            INDLL = .FALSE.
         END IF
      END IF
      RETURN
      END
C
C*************************************************************************
      SUBROUTINE GET_DEPTH(DEPTHNOW,J,XDR,YDR,LL,H0)
C*************************************************************************
C
C *** RETURNS THE DEPTH H OF LEVEL LL DIRECTLY BELOW DROGUE POSITION
C *** XDR,YDR, IN ELEMENT J.
C
C------------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
      COMMON /ELEM/ELEMS
      INTEGER ELEMS(1:NEDIM,1:3)
      COMMON /COORDS/XND,YND
      REAL XND(1:NNDIM),YND(1:NNDIM)
      REAL XJ(1:3),YJ(1:3),ZJ(1:3),DEPTHNOW(1:3,1:NNVDIM)
C
C *** LOCAL VARIABLE DECLARATION
C
      REAL XDR,YDR,H0
      REAL A,B,C,D,TRASH
      INTEGER J,N1,N2,N3,LL
C
C *** CALL SUBROUTINE NORMAL TO CALCULATE THE COEFFICIENTS OF THE PLANE
C *** CONTAINING ELEMENT J.   A,B,C,D ARE THE COEFFICIENTS OF THE 
C *** PLANE AX+BY+CZ+D=0 AS WELL AS THE COEFFICIENTS OF A VECTOR
C *** AI+BJ+CK NORMAL TO THE PLANE AX+BY+CZ+D=0.  THE SCRATCH VARIABLE
C *** 'TRASH' CONTAINS THE MAGNITUDE OF THE NORMAL AI+BJ+CK, WHICH IS
C *** NOT USED HERE.
C
      N1=ELEMS(J,1)  
      N2=ELEMS(J,2)  
      N3=ELEMS(J,3)  
      XJ(1)=XND(N1)
      XJ(2)=XND(N2)
      XJ(3)=XND(N3)
      YJ(1)=YND(N1)
      YJ(2)=YND(N2)
      YJ(3)=YND(N3)
      ZJ(1)=DEPTHNOW(1,LL)
      ZJ(2)=DEPTHNOW(2,LL)
      ZJ(3)=DEPTHNOW(3,LL)
      
      CALL NORMAL(XJ,YJ,ZJ,A,B,C,D,TRASH)
C
C *** COMPUTATION OF DEPTH H BY EQUATION OF PLANE CONTAINING ELEMENT J
C *** N1 IS A NODE IN ELEMENT J, XJ(1),YJ(1),ZJ(1) ARE THE COORDINATES 
C *** OF N1 IN J AT LEVEL LL.
C
      H0=-(A*(XDR-XJ(1))+B*(YDR-YJ(1)))/C+ZJ(1)
      RETURN
      END
C
C*************************************************************************
C THIS SUBROUTINE IS THE CODE
C ADDITION FOR DYNAMIC VERTICAL-GRID LEVEL FINDING.  
C 14 MAR 1994, BRIAN BLANTON
      SUBROUTINE DEPTH_INTERP(NNV,T1,T2,DEPTH1,DEPTH2,J,TNOW,DEPTHNOW)
C
C*************************************************************************
      INCLUDE 'DCMS.DIM'
      REAL DEPTH1(1:NNDIM,1:NNVDIM),DEPTH2(1:NNDIM,1:NNVDIM)
      REAL T1,T2,TNOW,TSLIDE
      INTEGER NNV
      COMMON /ELEM/ELEMS
      INTEGER ELEMS(1:NEDIM,1:3)
      INTEGER J,L,N1,N2,N3
      REAL DEPTHNOW(1:3,1:NNVDIM)
C
C *** GIVEN ELEMENT J'S VERTICAL STRUCTURE AT T1 AND T2, INTERPOLATE 
C *** THE VERTICAL STRUCTURE TO TNOW.  
C
      N1=ELEMS(J,1)
      N2=ELEMS(J,2)
      N3=ELEMS(J,3)
      TSLIDE=(TNOW-T1)/(T2-T1) 
C
C *** LOOP OVER NNV FOR EACH OF THE THREE NODES IN ELEMENT J AND 
C *** INTERPOLATE THE VERTICAL NODE POSITIONS TO TNOW.  
C
      DO L=1,NNV
         DEPTHNOW(1,L)=(1.E0-TSLIDE)*DEPTH1(N1,L)+
     +                 TSLIDE*DEPTH2(N1,L) 
         DEPTHNOW(2,L)=(1.E0-TSLIDE)*DEPTH1(N2,L)+
     +                 TSLIDE*DEPTH2(N2,L) 
         DEPTHNOW(3,L)=(1.E0-TSLIDE)*DEPTH1(N3,L)+
     +                 TSLIDE*DEPTH2(N3,L) 
      END DO
      RETURN
      END
C
C***************************************************************************
      SUBROUTINE NORMAL(XJ,YJ,ZJ,AAA,BBB,CCC,DDD,UNITNORM)
C***************************************************************************
C
C THIS SUBROUTINE COMPUTES THE COEFFICIENTS OF A VECTOR NORMAL
C TO THE PLANE CONTAINING THE THREE VERTICES OF ELEMENT J.  TO DO SO,
C THE COEFFICIENTS OF THE PLANE CONTAINING ELEMENT J ARE CALCULATED TO GIVE
C                   AAAX + BBBY + CCCZ + DDD = 0.
C THE NUMBERS AAA,BBB,CCC,DDD ARE ALSO THE COEFFICIENTS OF THE VECTOR 
C                   AAAI + BBBJ + CCCK,
C WHICH IS NORMAL TO AAAX + BBBY + CCCZ + DDD = 0.
C NORMAL ALSO RETURNS THE MAGNITUDE OF AAAI + BBBJ + CCCK, CALLED UNITNORM,  
C SO THAT A UNIT VECTOR MAY BE CALCULATED.
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
      REAL XJ(1:3),YJ(1:3),ZJ(1:3)
      REAL NEAR_ZERO
      PARAMETER (NEAR_ZERO=1.0E-4)
      REAL AAA,BBB,CCC,DDD,UNITNORM
      COMMON TRACKFLAG
      CHARACTER*3 TRACKFLAG
C
C *** CALCULATION OF COEFFICIENTS AAA,BBB,CCC,DDD FOR PLANE 
C *** AAAX+BBBY+CCCZ+DDD = 0 CONTAINING ELEMENT J.
C
      AAA = (YJ(2)-YJ(1))*(ZJ(3)-ZJ(1))-(YJ(3)-YJ(1))*(ZJ(2)-ZJ(1))
      BBB = (ZJ(2)-ZJ(1))*(XJ(3)-XJ(1))-(ZJ(3)-ZJ(1))*(XJ(2)-XJ(1))
      CCC = (XJ(2)-XJ(1))*(YJ(3)-YJ(1))-(XJ(3)-XJ(1))*(YJ(2)-YJ(1))
      IF(CCC.LT.NEAR_ZERO) THEN
	  WRITE(13,2000) 
          WRITE(13,8050)
          WRITE(*,8050)
          TRACKFLAG='NO'
          RETURN
      END IF	
      DDD = -AAA*XJ(1)-BBB*YJ(1)-CCC*ZJ(1)
C
C *** VECTOR NORMAL TO PLANE AX+BY+CZ+D = 0 IS N=AAAI+BBBJ+CCCK 
C
C *** UNIT NORMAL IS N/|N|;  |N| IS CALLED 'UNITNORM'
C
      UNITNORM = SQRT(AAA*AAA+BBB*BBB+CCC*CCC)
      RETURN
 2000 FORMAT('INTERPOLATION MATRIX SINGULAR IN ROUTINE GET_DEPTH.',
     +' DIVISION BY ~ ZERO ON NEXT CALCULATION. PROGRAM HALTED')
 8050 FORMAT('TRACKING FLAG SET TO NO IN DROG3DDT INPUT ROUTINE.')

	END
C
C***********************************************************************
      SUBROUTINE RK4(NNV,T1,T2,DT,DEPTH1,DEPTH2,
     +               UT1,UT2,VT1,VT2,WT1,WT2,      
     +               J,JOUT,L,LOUT,XX,YY,ZZ,U,V,W,
     +               XOUT,YOUT,ZOUT,IBOUN)
C***********************************************************************
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
C *** USES 4TH ORDER RUNGE-KUTTA SCHEME TO ADVANCE SOLUTION OVER A TIME
C *** INTERVAL "DT" AND RETURNS THE RESULTING POINT (XOUT,YOUT,ZOUT) 
C *** AS THE COORDINATE OF THE ENDPOINT OF THE INTEGRATION STEP.
C
C *** LOCAL VARIABLE TYPE DECLARATION 
C
      REAL UT1(1:NNDIM,1:NNVDIM),UT2(1:NNDIM,1:NNVDIM)
      REAL VT1(1:NNDIM,1:NNVDIM),VT2(1:NNDIM,1:NNVDIM)
      REAL WT1(1:NNDIM,1:NNVDIM),WT2(1:NNDIM,1:NNVDIM)
      REAL XX,YY,ZZ,U,V,W,DT,XOUT,YOUT,ZOUT,DTH,DT6
      REAL XT,YT,ZT,UT,VT,WT,TH,UM,VM,WM,U4,V4,W4,T1,T2
      REAL DEPTH1(1:NNDIM,1:NNVDIM),DEPTH2(1:NNDIM,1:NNVDIM)
      INTEGER IND,INDD,J,JOUT,L,LOUT,JT,LT,JJ,LL,NNV,IBOUN
      LOGICAL INDLL,BOTFLAG  
C
      INDLL = .TRUE.
      DTH = DT/2.E0
      DT6 = DT/6.E0
      TH = T1 + DTH
      JT = J
      LT = L
C
C *** FIRST STEP
C
C
C GET COMPONENTS OF FLOW U,V,W AT XX,YY,ZZ.  THIS IS THE INITIAL
C CONDITION FOR THE RUNGE-KUTTA TIMESTEP
C
      CALL VEL_INTERP(NNV,T1,T2,DEPTH1,DEPTH2,
     +                UT1,UT2,VT1,VT2,WT1,WT2,      
     +                JT,LT,XX,YY,ZZ,U,V,W,T1)
      
      XT = XX + DTH*U
      YT = YY + DTH*V
      ZT = ZZ + DTH*W
C
C *** CHECK WHETHER POINT STILL LIES IN ELEMENT "J";  
C *** IF NOT, FIND NEW ELEMENT
C
      CALL BELEL(JT,XT,YT,IND)
      IF (IND.EQ.0) THEN
          CALL FIND_ELEMENT(JT,JJ,XT,YT,INDD)
          IF (INDD.EQ.0) THEN
             IBOUN=1
	     WRITE(13,*)
     +       'RK-1 DROGUE EXITED HORZ AT T (DAYS) =',T1/86400.
             RETURN
          END IF
          JT = JJ
      END IF
      CALL FIND_LEVEL(NNV,T1,T2,DEPTH1,DEPTH2,
     +                JT,XT,YT,ZT,TH,LT,LL,INDLL,BOTFLAG)
      IF(.NOT.INDLL) THEN
	 XOUT=XX
	 YOUT=YY
	 ZOUT=ZZ
	 JOUT=JT
	 LOUT=LT
	 IBOUN=1
	 WRITE(13,*)
     +   'RK-1 DROGUE EXITED FREE SURFACE AT T (DAYS) =',T1/86400.
         RETURN
      ENDIF
      IF(BOTFLAG)THEN
         GOTO 1000
      END IF
      LT=LL
C
C *** SECOND STEP
C
      CALL VEL_INTERP(NNV,T1,T2,DEPTH1,DEPTH2,
     +                UT1,UT2,VT1,VT2,WT1,WT2,      
     +                JT,LT,XT,YT,ZT,UT,VT,WT,TH)
      XT = XX + DTH*UT
      YT = YY + DTH*VT
      ZT = ZZ + DTH*WT
      CALL BELEL(JT,XT,YT,IND)
      IF (IND.EQ.0) THEN
          CALL FIND_ELEMENT(JT,JJ,XT,YT,INDD)
          IF (INDD.EQ.0) THEN
             IBOUN=1
	     WRITE(13,*)
     +       'RK-2 DROGUE EXITED HORZ AT T (DAYS) =',T1/86400.
             RETURN
          END IF
          JT = JJ
      END IF
      CALL FIND_LEVEL(NNV,T1,T2,DEPTH1,DEPTH2,
     +                JT,XT,YT,ZT,TH+DTH,LT,LL,INDLL,BOTFLAG)
      IF(.NOT.INDLL) THEN
         XOUT=XT
         YOUT=YT
         ZOUT=ZT
         JOUT=JT
         LOUT=LT
	 IBOUN=1
	 WRITE(13,*)
     +'RK-2 DROGUE EXITED FREE SURFACE AT T (DAYS) =',T1/86400.
         RETURN
      ENDIF
      IF(BOTFLAG)THEN
         GOTO 1000
      END IF
      LT=LL
C
C *** THIRD STEP
C
      CALL VEL_INTERP(NNV,T1,T2,DEPTH1,DEPTH2,
     +                UT1,UT2,VT1,VT2,WT1,WT2,      
     +                JT,LT,XT,YT,ZT,UM,VM,WM,TH)
      XT = XX + DT*UM
      YT = YY + DT*VM
      ZT = ZZ + DT*WM
C
      CALL BELEL(JT,XT,YT,IND)
      IF (IND.EQ.0) THEN
          CALL FIND_ELEMENT(JT,JJ,XT,YT,INDD)
          IF (INDD.EQ.0) THEN
             IBOUN=1
	     WRITE(13,*)
     +       'RK-3 DROGUE EXITED HORZ AT T (DAYS) =',T1/86400.
             RETURN
          END IF
          JT = JJ
      END IF
      CALL FIND_LEVEL(NNV,T1,T2,DEPTH1,DEPTH2,
     +                JT,XT,YT,ZT,TH+DT,LT,LL,INDLL,BOTFLAG)
      IF(.NOT.INDLL) THEN
	 XOUT=XT
         YOUT=YT
         ZOUT=ZT
	 JOUT=JT
	 LOUT=LT
	 IBOUN=1
	 WRITE(13,*)
     +   'RK-3 DROGUE EXITED FREE SURFACE AT T (DAYS) =',T1/86400.
         RETURN
      END IF
      IF(BOTFLAG)THEN
         GOTO 1000
      END IF
      LT=LL
C
C *** FOURTH STEP
C
      CALL VEL_INTERP(NNV,T1,T2,DEPTH1,DEPTH2,
     +                UT1,UT2,VT1,VT2,WT1,WT2,      
     +                JT,LT,XT,YT,ZT,U4,V4,W4,T1+DT)
C
C *** ACCUMULATE INCREMENTS WITH PROPER WEIGHTS
C
      XOUT = XX + (U+2.E0* (UT+UM)+U4)*DT6
      YOUT = YY + (V+2.E0* (VT+VM)+V4)*DT6
      ZOUT = ZZ + (W+2.E0* (WT+WM)+W4)*DT6

      CALL BELEL(JT,XOUT,YOUT,IND)
      IF (IND.EQ.0) THEN
          CALL FIND_ELEMENT(JT,JJ,XOUT,YOUT,INDD)
          IF (INDD.EQ.0) THEN
             IBOUN=1
	     WRITE(13,*)
     +       'RK-4 DROGUE EXITED HORZ AT T (DAYS) =',T1/86400.
             RETURN
          END IF
          JT = JJ
      END IF
      CALL FIND_LEVEL(NNV,T1,T2,DEPTH1,DEPTH2,
     +                JT,XOUT,YOUT,ZOUT,T1+DT,LT,LL,INDLL,BOTFLAG)
      IF(.NOT.INDLL) THEN
         XOUT=XT
         YOUT=YT
         ZOUT=ZT
         JOUT=JT
         LOUT=LT
	 IBOUN=1
	 WRITE(13,*)
     +   'RK-4 DROGUE EXITED FREE SURFACE AT T (DAYS) ',T1/86400.
         RETURN
      END IF
      IF(BOTFLAG)THEN
         GOTO 1000
      END IF

      JOUT = JT
      LOUT = LT

      RETURN
 1000 CONTINUE
      WRITE(13,*)'BOTTOM ENCOUNTERED AT T=',T1/86400.
      CALL HANDLE_BOTTOM(NNV,T1,T2,DEPTH1,DEPTH2,
     +            UT1,UT2,VT1,VT2,WT1,WT2,      
     +            J,JOUT,L,LOUT,XX,YY,ZZ,
     +            XOUT,YOUT,ZOUT)
      RETURN         
       
      END
C
C************************************************************************
      SUBROUTINE VEL_INTERP(NNV,T1,T2,DEPTH1,DEPTH2,
     +                      UT1,UT2,VT1,VT2,WT1,WT2,
     +                      J,L,XDR,YDR,ZDR,UDR,VDR,WDR,T)
C************************************************************************
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
C *** LOCAL VARIABLE TYPE DECLARATION
C
      REAL DEPTH1(1:NNDIM,*),DEPTH2(1:NNDIM,*)
      REAL XDR,YDR,ZDR,UDR,VDR,WDR,T1,T2,DZ,T
      REAL ULOW,UUP,VLOW,VUP,WUP,WLOW
      REAL ZLOW,ZUP,WEIGHTUP,WEIGHTLOW
      REAL UFUN,VFUN,WFUN
      REAL DEPTHNOW(1:3,NNVDIM)
      INTEGER J,L,NNV
      REAL UT1(1:NNDIM,*),UT2(1:NNDIM,*)
      REAL VT1(1:NNDIM,*),VT2(1:NNDIM,*)
      REAL WT1(1:NNDIM,*),WT2(1:NNDIM,*)
C
C *** GET U,V,W COMPONENTS OF FLOW AT LEVELS ABOVE AND BELOW (X,Y,Z)
C
      ULOW =UFUN(T1,T2,UT1,UT2,J,L,XDR,YDR,ZDR,T)
      UUP  =UFUN(T1,T2,UT1,UT2,J,L+1,XDR,YDR,ZDR,T)
      VLOW =VFUN(T1,T2,VT1,VT2,J,L,XDR,YDR,ZDR,T)
      VUP  =VFUN(T1,T2,VT1,VT2,J,L+1,XDR,YDR,ZDR,T)
      WLOW =WFUN(T1,T2,WT1,WT2,J,L,XDR,YDR,ZDR,T)
      WUP  =WFUN(T1,T2,WT1,WT2,J,L+1,XDR,YDR,ZDR,T)
C
C *** GET DEPTHS AT LEVELS DIRECTLY ABOVE (L+1) 
C *** AND BELOW(L) DROGUE POSITION
C
      CALL DEPTH_INTERP(NNV,T1,T2,DEPTH1,DEPTH2,J,T,DEPTHNOW)
      CALL GET_DEPTH(DEPTHNOW,J,XDR,YDR,L,ZLOW)	
      CALL GET_DEPTH(DEPTHNOW,J,XDR,YDR,L+1,ZUP)	
C
C *** COMPUTE VERTICAL SPACING BETWEEN LEVELS 
C *** AT THIS (XDR,YDR) LOCATION
C
      DZ=ABS(ZUP-ZLOW)
C
C *** COMPUTE "WEIGHTS" FOR LINEAR INTERPOLATION OF VELOCITY COMPONENTS 
C
      WEIGHTUP=ABS(ZDR-ZLOW)/DZ
      WEIGHTLOW=ABS(ZDR-ZUP)/DZ
C
C *** LINEARLY INTERPOLATE FIELD VELOCITIES TO (XDR,YDR,ZDR)
C 
      UDR = WEIGHTUP*UUP+WEIGHTLOW*ULOW
      VDR = WEIGHTUP*VUP+WEIGHTLOW*VLOW
      WDR = WEIGHTUP*WUP+WEIGHTLOW*WLOW
 1000 FORMAT(A3,1X,6(E14.6,1X),1X,F10.3)
C
      RETURN
      END
C
C************************************************************************
      SUBROUTINE BOT_VEL(T1,T2,UT1,UT2,VT1,VT2,WT1,WT2,
     +                   J,XDR,YDR,ZDR,UDR,VDR,WDR,T)
C************************************************************************
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
C *** LOCAL VARIABLE TYPE DECLARATION
C
      REAL XDR,YDR,ZDR,UDR,VDR,WDR,T
      REAL UFUN,VFUN,WFUN,T1,T2
      INTEGER J
      REAL UT1(1:NNDIM,1:NNVDIM),UT2(1:NNDIM,1:NNVDIM)
      REAL VT1(1:NNDIM,1:NNVDIM),VT2(1:NNDIM,1:NNVDIM)
      REAL WT1(1:NNDIM,1:NNVDIM),WT2(1:NNDIM,1:NNVDIM)
      
      UDR=UFUN(T1,T2,UT1,UT2,J,1,XDR,YDR,ZDR,T)
      VDR=VFUN(T1,T2,VT1,VT2,J,1,XDR,YDR,ZDR,T)
      WDR=WFUN(T1,T2,WT1,WT2,J,1,XDR,YDR,ZDR,T)
      RETURN
      END
C
C************************************************************************
      SUBROUTINE SCA_INTERP(NNV,T1,T2,Z1,Z2,S1,S2,
     +                      J,L,XDR,YDR,ZDR,SDR,T)
C************************************************************************
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
C *** LOCAL VARIABLE TYPE DECLARATION
C
      REAL Z1(1:NNDIM,*),Z2(1:NNDIM,*)
      REAL S1(1:NNDIM,*),S2(1:NNDIM,*)
      REAL XDR,YDR,ZDR,SDR,T1,T2,DZ,T
      REAL SLOW,SUP
      REAL ZLOW,ZUP,WEIGHTUP,WEIGHTLOW
      REAL SFUN
      REAL ZNOW(1:3,NNVDIM)
      INTEGER J,L,NNV
C
C *** GET S COMPONENTS OF SCALAR AT LEVELS ABOVE AND BELOW (X,Y,Z)
C
      SLOW =SFUN(T1,T2,S1,S2,J,L,XDR,YDR,ZDR,T)
      SUP  =SFUN(T1,T2,S1,S2,J,L+1,XDR,YDR,ZDR,T)
C
C *** GET DEPTHS AT LEVELS DIRECTLY ABOVE (L+1) 
C *** AND BELOW(L) DROGUE POSITION
C
      CALL DEPTH_INTERP(NNV,T1,T2,Z1,Z2,J,T,ZNOW)
      CALL GET_DEPTH(ZNOW,J,XDR,YDR,L,ZLOW)	
      CALL GET_DEPTH(ZNOW,J,XDR,YDR,L+1,ZUP)	
C
C *** COMPUTE VERTICAL SPACING BETWEEN LEVELS 
C *** AT THIS (XDR,YDR) LOCATION
C
      DZ=ABS(ZUP-ZLOW)
C
C *** COMPUTE "WEIGHTS" FOR LINEAR INTERPOLATION OF VELOCITY COMPONENTS 
C
      WEIGHTUP=ABS(ZDR-ZLOW)/DZ
      WEIGHTLOW=ABS(ZDR-ZUP)/DZ
C
C *** LINEARLY INTERPOLATE FIELD SCALAR TO (XDR,YDR,ZDR)
C 
      SDR = WEIGHTUP*SUP+WEIGHTLOW*SLOW
C
      RETURN
      END
C
C***********************************************************************
      REAL FUNCTION SFUN(T1,T2,S1,S2,J,L,X,Y,Z,T)
C***********************************************************************
C
C       DETERMINES THE LOCAL VALUE OF S AT (X,Y,L)
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
      REAL T1,T2
      COMMON /ARJ/AR
      COMMON /ABA0/A,B,A0
      REAL*8 A(1:NEDIM,3),B(1:NEDIM,3),A0(1:NEDIM,2)
      REAL*8 AR(1:NEDIM)
      COMMON /ELEM/ELEMS
      INTEGER ELEMS(1:NEDIM,3)
      REAL S1(1:NNDIM,1:NNVDIM),S2(1:NNDIM,1:NNVDIM)
      INTEGER N1,N2,N3,J,L
      REAL X,Y,Z,SDR1,SDR2,SDR3,TSLIDE,A03,ARI,C1,C2,C3,T
      
      N1=ELEMS(J,1)
      N2=ELEMS(J,2)
      N3=ELEMS(J,3)
      A03 = AR(J) - A0(J,1) - A0(J,2)
      ARI = 0.5/AR(J)
C
C *** LINEARLY WEIGHT THE VELOCITY AT TIME T BETWEEN VELOCITIES
C *** AT TIME T1 AND TIME T2.
C
      TSLIDE=(T-T1)/(T2-T1)
      SDR1 = (1.E0-TSLIDE)*S1(N1,L)+TSLIDE*S2(N1,L)
      SDR2 = (1.E0-TSLIDE)*S1(N2,L)+TSLIDE*S2(N2,L)
      SDR3 = (1.E0-TSLIDE)*S1(N3,L)+TSLIDE*S2(N3,L)
        
      C1 = ARI* (B(J,1)*SDR1+B(J,2)*SDR2+B(J,3)*U3)
      C2 = ARI* (A(J,1)*SDR1+A(J,2)*SDR2+A(J,3)*U3)
      C3 = 2*ARI* (A0(J,1)*SDR1+A0(J,2)*SDR2+A03*U3)
        
      SFUN = C1*X + C2*Y + C3
      RETURN

      END
C
C***********************************************************************
      REAL FUNCTION UFUN(T1,T2,UT1,UT2,J,L,X,Y,Z,T)
C***********************************************************************
C
C       DETERMINES U COMPONENT OF VELOCITY AT (X,Y,L)
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
      REAL T1,T2
      COMMON /ARJ/AR
      COMMON /ABA0/A,B,A0
      REAL*8 A(1:NEDIM,3),B(1:NEDIM,3),A0(1:NEDIM,2)
      REAL*8 AR(1:NEDIM)
      COMMON /ELEM/ELEMS
      INTEGER ELEMS(1:NEDIM,3)
      REAL UT1(1:NNDIM,1:NNVDIM),UT2(1:NNDIM,1:NNVDIM)
      INTEGER N1,N2,N3,J,L
      REAL X,Y,Z,U1,U2,U3,TSLIDE,A03,ARI,C1,C2,C3,T
      
      N1=ELEMS(J,1)
      N2=ELEMS(J,2)
      N3=ELEMS(J,3)
      A03 = AR(J) - A0(J,1) - A0(J,2)
      ARI = 0.5/AR(J)
C
C *** LINEARLY WEIGHT THE VELOCITY AT TIME T BETWEEN VELOCITIES
C *** AT TIME T1 AND TIME T2.
C
      TSLIDE=(T-T1)/(T2-T1)
      U1 = (1.E0-TSLIDE)*UT1(N1,L)+TSLIDE*UT2(N1,L)
      U2 = (1.E0-TSLIDE)*UT1(N2,L)+TSLIDE*UT2(N2,L)
      U3 = (1.E0-TSLIDE)*UT1(N3,L)+TSLIDE*UT2(N3,L)
        
      C1 = ARI* (B(J,1)*U1+B(J,2)*U2+B(J,3)*U3)
      C2 = ARI* (A(J,1)*U1+A(J,2)*U2+A(J,3)*U3)
      C3 = 2*ARI* (A0(J,1)*U1+A0(J,2)*U2+A03*U3)
        
      UFUN = C1*X + C2*Y + C3
      RETURN

      END

C
C***********************************************************************
      REAL FUNCTION VFUN(T1,T2,VT1,VT2,J,L,X,Y,Z,T)
C***********************************************************************
C
C       COMPUTES V COMPONENT OF VELOCITY AT (X,Y,L)
C
C
C ***   SET MAXIMUM ARRAY DIMENSIONS
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
      REAL T1,T2
      COMMON /ARJ/AR
      COMMON /ABA0/A,B,A0
      REAL*8  A(1:NEDIM,3),B(1:NEDIM,3),A0(1:NEDIM,2)
      REAL*8 AR(1:NEDIM)
      COMMON /ELEM/ELEMS
      INTEGER ELEMS(1:NEDIM,3)
      REAL VT1(1:NNDIM,1:NNVDIM),VT2(1:NNDIM,1:NNVDIM)
      INTEGER N1,N2,N3,J,L
      REAL X,Y,Z,V1,V2,V3,TSLIDE,A03,ARI,D1,D2,D3,T
      
      N1=ELEMS(J,1)
      N2=ELEMS(J,2)
      N3=ELEMS(J,3)
      A03 = AR(J) - A0(J,1) - A0(J,2)
      ARI = 0.5/AR(J)
C
C *** LINEARLY WEIGHT THE VELOCITY AT TIME T BETWEEN VELOCITIES
C     AT TIME T1 AND TIME T2.
C
      TSLIDE=(T-T1)/(T2-T1)
      V1 = (1.E0-TSLIDE)*VT1(N1,L)+TSLIDE*VT2(N1,L)
      V2 = (1.E0-TSLIDE)*VT1(N2,L)+TSLIDE*VT2(N2,L)
      V3 = (1.E0-TSLIDE)*VT1(N3,L)+TSLIDE*VT2(N3,L)

      D1 = ARI* (B(J,1)*V1+B(J,2)*V2+B(J,3)*V3)
      D2 = ARI* (A(J,1)*V1+A(J,2)*V2+A(J,3)*V3)
      D3 = 2*ARI* (A0(J,1)*V1+A0(J,2)*V2+A03*V3)
       
      VFUN = D1*X + D2*Y + D3
      RETURN

      END
C
C***********************************************************************
      REAL FUNCTION WFUN(T1,T2,WT1,WT2,J,L,X,Y,Z,T)
C***********************************************************************
C
C       COMPUTES W COMPONENT OF VELOCITY AT (X,Y,L)
C
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
C-----------------------------------------------------------------------
C
      REAL T1,T2
      COMMON /ARJ/AR
      COMMON /ABA0/A,B,A0
      REAL*8 A(1:NEDIM,3),B(1:NEDIM,3),A0(1:NEDIM,2)
      REAL*8 AR(1:NEDIM)
      COMMON /ELEM/ELEMS
      INTEGER ELEMS(1:NEDIM,3)
      REAL WT1(1:NNDIM,1:NNVDIM),WT2(1:NNDIM,1:NNVDIM)
      INTEGER N1,N2,N3,J,L
      REAL X,Y,Z,W1,W2,W3,TSLIDE,A03,ARI,E1,E2,E3,T
      
      N1=ELEMS(J,1)
      N2=ELEMS(J,2)
      N3=ELEMS(J,3)
      A03 = AR(J) - A0(J,1) - A0(J,2)
      ARI = 0.5/AR(J)
C
C *** LINEARLY WEIGHT THE VELOCITY AT TIME T BETWEEN VELOCITIES
C     AT TIME T1 AND TIME T2.
C
      TSLIDE=(T-T1)/(T2-T1)
      W1 = (1.E0-TSLIDE)*WT1(N1,L)+TSLIDE*WT2(N1,L)
      W2 = (1.E0-TSLIDE)*WT1(N2,L)+TSLIDE*WT2(N2,L)
      W3 = (1.E0-TSLIDE)*WT1(N3,L)+TSLIDE*WT2(N3,L)

      E1 = ARI* (B(J,1)*W1+B(J,2)*W2+B(J,3)*W3)
      E2 = ARI* (A(J,1)*W1+A(J,2)*W2+A(J,3)*W3)
      E3 = 2*ARI* (A0(J,1)*W1+A0(J,2)*W2+A03*W3)
       
      WFUN = E1*X + E2*Y + E3
      RETURN

      END
C
C**************************************************************************
      SUBROUTINE HANDLE_BOTTOM(NNV,T1,T2,DEPTH1,DEPTH2,
     +                         UT1,UT2,VT1,VT2,WT1,WT2,      
     +                         J,JOUT,L,LOUT,
     +                         XIN,YIN,ZIN,XOUT,YOUT,ZOUT) 
C**************************************************************************
C
      INCLUDE 'DCMS.DIM'     
      REAL XIN,YIN,ZIN,XOUT,YOUT,ZOUT,T1,T2
      INTEGER NNV,J,JOUT,L,LOUT
      REAL UT1(1:NNDIM,1:NNVDIM),UT2(1:NNDIM,1:NNVDIM)
      REAL VT1(1:NNDIM,1:NNVDIM),VT2(1:NNDIM,1:NNVDIM)
      REAL WT1(1:NNDIM,1:NNVDIM),WT2(1:NNDIM,1:NNVDIM)
      REAL DEPTH1(1:NNDIM,1:NNVDIM),DEPTH2(1:NNDIM,1:NNVDIM)
      
      LOGICAL BOTFLAG,LLIND 
      REAL UB,VB,WB
      INTEGER IND,INDD,IBOUN
      REAL DEPTHNOW(1:3,1:NNVDIM),ZBOT
              
      CALL BOT_VEL(T1,T2,UT1,UT2,VT1,VT2,WT1,WT2,
     +             J,XIN,YIN,ZIN,UB,VB,WB,T1)
      XOUT = XIN + (T2-T1) * UB
      YOUT = YIN + (T2-T1) * VB
      
      CALL BELEL(J,XOUT,YOUT,IND)
      IF (IND.EQ.0) THEN
          CALL FIND_ELEMENT(J,JOUT,XOUT,YOUT,INDD)
          IF (INDD.EQ.0) THEN
             IBOUN=1
	     WRITE(13,*)
     +       'HB- DROGUE EXITED HORZ AT T (DAYS) =',T1/86400.
             RETURN
          END IF
      END IF
         
      CALL DEPTH_INTERP(NNV,T1,T2,DEPTH1,DEPTH2,JOUT,T2,DEPTHNOW)
      CALL GET_DEPTH(DEPTHNOW,JOUT,XOUT,YOUT,1,ZBOT)
      ZOUT=ZBOT+0.01
      CALL FIND_LEVEL(NNV,T1,T2,DEPTH1,DEPTH2,
     +                JOUT,XOUT,YOUT,ZOUT,T2,
     +                1,LOUT,LLIND,BOTFLAG)
      IF(BOTFLAG)THEN
         PRINT*,'CODE STOPPED IN HANDLE_BOTTOM'
         STOP
      END IF
      RETURN
      END        
C
C**************************************************************************
      SUBROUTINE OUTPUT(NNV,T1,T2,Z1,Z2,
     +                  IDR,JJDR,LLDR,NDR,XDR,YDR,ZDR)
C**************************************************************************
C
      INCLUDE 'DCMS.DIM'
C
C *** NNV - NUMBER OF VERTICAL NODES
C *** T1,T2 - BEGINNING AND ENDING TIME OF THIS TIME-STEP IN SEC
C *** Z1,Z2 - VERTICAL NODE POSITIONS AT T1 AND T2
C *** IDR - ARRAY CONTAINING THE I-TH DROGUE'S CURRENT STATUS; 
C           I.E.,  1 IF BEING TRACKED, 0 IF NO LONGER IN THE DOMAIN
C *** JJDR - ARRAY CONTAINING CURRENT HORIZONTAL ELEMENT NUMBER OF DROGUE I
C *** LLDR - ARRAY CONTAINING LEVEL NUMBER DIRECTLY BELOW DROGUE I;
C           I.E., IF DROGUE I IS IN LOWEST LEVEL, LLDR(I) = 1, AND IF
C           DROGUE I IS IN THE SURFACE LEVEL, LLDR(I) = NNV-1
C *** NDR - NUMBER OF DROGUES AT START OF COMPUTATION
C *** XDR,YDR,ZDR - ARRAYS CONTAINING CURRENT (X,Y,Z) LOCATION OF DROGUE I
C
      INTEGER NNV,NDR
      REAL T1,T2
      REAL Z1(1:NNDIM,*),Z2(1:NNDIM,*)
      INTEGER IDR(*),JJDR(*),LLDR(*)
      REAL XDR(*),YDR(*),ZDR(*)
C
C *** LOCAL VARIABLES
C
C *** I - LOOP COUNTER FOR DROGUES
C *** ZBOT - DEPTH OF GRID AT LOWER MOST LEVEL (NNV=1) OVER (XDR(I),YDR(I))
C *** DEPTHNOW - INTERPOLATED DEPTHS OF NODES AT ELEMENT JJDR AT TIME T2
C
      INTEGER I
      REAL ZBOT
      REAL ZNOW(1:3,1:NNVDIM)
C
C *** OPEN '.LTS' (LAST TIME-STEP) FILE FOR HOT START, UNIT 17 
C
      OPEN(UNIT=17,FILE='DROG.LTS')
C
C *** LOOP OVER DROGUES TO OUTPUT
C
      DO I=1,NDR
        CALL DEPTH_INTERP(NNV,T1,T2,Z1,Z2,JJDR(I),T2,ZNOW)
        WRITE(12,9000) XDR(I),YDR(I),ZDR(I),T2
        WRITE(17,9000) XDR(I),YDR(I),ZDR(I),T2
      END DO
C
C *** CLOSE THE DROGUE HOT START FILE
C
      CLOSE(UNIT=17)
      RETURN
 9000 FORMAT (5(2X,F12.4),I4)
      END
C
C-----------------------------------------------------------------------
C
      FUNCTION GASDEV(IDUM)
C
C GASDEV - GASDEV RETURNS A NORMALLY DISTRIBUTED DEVIATE
C          WITH ZERO MEAN AND UNIT VARIANCE, USING RAN1(IDUM) AS
C          THE SOURCE FOR THE UNIFORM DEVIATES.
C
      INTEGER IDUM
      REAL GASDEV
CU    USES RAN1 or RAN2
      INTEGER ISET
      REAL FAC,GSET,RSQ,V1,V2,RAN1,RAN2
      SAVE ISET,GSET
      DATA ISET/0/
      IF (ISET.EQ.0) THEN
C1       V1=2.*RAN1(IDUM)-1.
C        V2=2.*RAN1(IDUM)-1.
1       V1=2.*RAN2(IDUM)-1.
        V2=2.*RAN2(IDUM)-1.
        RSQ=V1**2+V2**2
        IF(RSQ.GE.1..OR.RSQ.EQ.0.)GOTO 1
        FAC=SQRT(-2.*LOG(RSQ)/RSQ)
        GSET=V1*FAC
        GASDEV=V2*FAC
        ISET=1
      ELSE
        GASDEV=GSET
        ISET=0
      ENDIF
      RETURN
      END
C  (C) COPR. 1986-92 NUMERICAL RECIPES SOFTWARE &OL`.
C
C-----------------------------------------------------------------------
C
      FUNCTION RAN1(IDUM)
C
C RAN1   - RAN1 RETURNS A UNIFORM RANDOM DEVIATE BETWEEN
C          0.0 AND 1.0, ENDPOINTS NOT INCLUDED.  IDUM IS SET BY A
C          CALL TO THE SYSTEM TIME FUNCTION (MULTIPLIED BY -1) TO
C          INITIALIZE THE RANDOM NUMBER SEQUENCE.  IDUM SHOULD NOT
C          BE CHANGED BETWEEN SUCCESSIVE DEVIATES IN A SEQUENCE.
C          THIS IS A MINIMAL RANDOM NUMBER GENERATOR.
C
      INTEGER IDUM,IA,IM,IQ,IR,NTAB,NDIV
      REAL RAN1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2E-7,RNMX=1.-EPS)
      INTEGER J,K,IV(NTAB),IY
      SAVE IV,IY
      DATA IV /NTAB*0/, IY /0/
      IF (IDUM.LE.0.OR.IY.EQ.0) THEN
        IDUM=MAX(-IDUM,1)
        DO 11 J=NTAB+8,1,-1
          K=IDUM/IQ
          IDUM=IA*(IDUM-K*IQ)-IR*K
          IF (IDUM.LT.0) IDUM=IDUM+IM
          IF (J.LE.NTAB) IV(J)=IDUM
11      CONTINUE
        IY=IV(1)
      ENDIF
      K=IDUM/IQ
      IDUM=IA*(IDUM-K*IQ)-IR*K
      IF (IDUM.LT.0) IDUM=IDUM+IM
      J=1+IY/NDIV
      IY=IV(J)
      IV(J)=IDUM
      RAN1=MIN(AM*IY,RNMX)
      RETURN
      END
C  (C) COPR. 1986-92 NUMERICAL RECIPES SOFTWARE &OL`.
c-----------------------------------------------------------------------
c
      FUNCTION ran2(idum)
c
c RAN2   - For situations requiring a longer random
c          sequence than that provided by ran1, ran2 returns a
c          uniform random deviate, between 0.0 and 1.0 (endpoints
c          not included), from a sequence with period ~2X10e18.
c
      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      REAL ran2,AM,EPS,RNMX
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
     *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
     *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
      if (idum.le.0) then
        idum=max(-idum,1)
        idum2=idum
        do 11 j=NTAB+8,1,-1
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if (idum.lt.0) idum=idum+IM1
      k=idum2/IQ2
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if (idum2.lt.0) idum2=idum2+IM2
      j=1+iy/NDIV
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+IMM1
      ran2=min(AM*iy,RNMX)
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software &OL`.
C
C******************* END OF EVERYTHING *************************C
C******************* END OF EVERYTHING *************************C
C******************* END OF EVERYTHING *************************C




CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C SUBROUTINE RWALKNODEAH(DT,XIN,YIN,JIN,AHI,XOUT,YOUT,JOUT,IBOUN,
C                        TIME,X,Y,EDGENODE,EDGENUM,IBC1,IBC2,IN,AR)
C
C     DT             - time step (seconds).
C     XIN,YIN        - input advected and undiffused location.
cfin     JIN            - input element # of drogue.
C     AHI            - uniform diffusivity from .inp file. 
C     XOUT,YOUT,JOUT - output advected and diffused location. 
C     JOUT           - output element # of drogue.
C     IBOUN          - Drogue within domain switch (0(in)/1(out)) 
C     TIME           - Simulation time.
C     X,Y            - Grid node coordinates.
C     EDGENODE       - Numbers of nodes forming a boundary edge.
C     EDGENUM        - Number of boundary edge.
C     IBC1           - Domain code left of segment (ccw).  
C     IBC2           - Domain code right of segment (ccw).
C   
C Author: J.A. Proehl, May 2003
C
C Modification by Ata Bilgili, May 16th, 2003.
C              - Convert to Double Precision to work with BELLAMY.
C              - Add element testing routine. 
c Keston Smith 11/26/2003
c              changed for non-uniform diffusion field
C              this version corrects random walk for 
C              gradient of diffusivity.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcC
C
      SUBROUTINE RWALKNODEAH(DT,XIN,YIN,JIN,AHI,XOUT,YOUT,JOUT,
     &                IBOUN,TIME,X,Y,EDGENODE,EDGENUM,IBC1,IBC2,
     &                IN,ar)

      IMPLICIT NONE
      INCLUDE 'DCMS.DIM'
C
      INTEGER JIN,JOUT,INDDRW,INDRW,NIXRW,EDGENUM,
     &          EDGENODE(2,*),IBC1(*),IBC2(*),IBOUN,IN(3,*)
      REAL DT,XIN,YIN,XOUT,YOUT,AHI(*),XSTEP,YSTEP,GASDEV,RX,RY,
     &       STDDEVRWX,STDDEVRWY,TIME,X(NNDIM),Y(NNDIM),ar(*)
      integer i1,i2,i3
      real AH
      real phi(3),xl(3),yl(3),dx(3),dy(3)
      real xstar,ystar,DAhDx,DAhDy
C
C... Initialize JOUT with JIN.
C
      JOUT=JIN
C
C... Random number generation. Uses Statpak2.f.
C 
      RX=gasdev()  
      RY=gasdev()  
      STDDEVRWX=1.d0
      STDDEVRWY=1.d0
      
      
      i1=in(1,Jin)
      i2=in(2,Jin)
      i3=in(3,Jin)

      xl(1)=x(i1)
      xl(2)=x(i2)
      xl(3)=x(i3)

      yl(1)=y(i1)
      yl(2)=y(i2)
      yl(3)=y(i3)
c
c compute gradient of AH with in the element
c
      DX(1)=X(I2)-X(I3)          
      DX(2)=X(I3)-X(I1)          
      DX(3)=X(I1)-X(I2)          

      DY(1)=Y(I2)-Y(I3)          
      DY(2)=Y(I3)-Y(I1)          
      DY(3)=Y(I1)-Y(I2)

      DAhDx= (AHI(i1)*Dy(1)+AHI(i2)*Dy(2)+AHI(i3)*Dy(3))/2./AR(Jin)          
      DAhDy=-(AHI(i1)*Dx(1)+AHI(i2)*Dx(2)+AHI(i3)*Dx(3))/2./AR(Jin)          
c
c done gradient calculation 
c      

C
C First order gradient corrected random Walk
C      

      Xstar=XIN+DAhDx*DT/2.D0
      Ystar=YIN+DAhDy*DT/2.D0
c
c kws compute local basis function
c
      call PHI2D(xl,yl,Xstar,Ystar,phi) 
c
c interpolate to local AH
c
      AH=AHI(i1)*phi(1)+AHI(i2)*phi(2)+AHI(I3)*phi(3)
      
      XSTEP=RX*DSQRT(2.d0*AH*DT/STDDEVRWX)+DAhDx*DT
      YSTEP=RY*DSQRT(2.d0*AH*DT/STDDEVRWY)+DAhDy*DT

      XOUT=XIN+XSTEP
      YOUT=YIN+YSTEP

C
C... Find the new element the drogue is in and update its number.
C... Also make sure that the element leaves the domain if it 
C... crosses an open ocean boundary and not a land boundary.
C... Ata Bilgili, May 16, 2003.
C
      CALL BELEL(JIN,XOUT,YOUT,INDRW)
       IF (INDRW.EQ.0) THEN
          CALL FIND_ELEMENT(JIN,JOUT,XOUT,YOUT,INDDRW)
          IF (INDDRW.EQ.0) THEN
      CALL OPEN_BND_CHECK(X,Y,XIN,YIN,XOUT,YOUT,NIXRW,EDGENODE,IBC1,
     &                    IBC2,EDGENUM)
		IF(NIXRW.EQ.3)THEN 
             IBOUN=1
	     WRITE(13,*)
     + 'RANDOM WALK DROGUE EXITED HORZ AT T (DAYS) =',TIME/86400.D0
             RETURN
                ELSE
             IBOUN=2 
                END IF
          JIN = JOUT
          END IF
      END IF
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SMAGOR1(U,V,AR,X,Y,IN,NE,NN,AHI,SMAGR,AHMIN)
C-----------------------------------------------------------------------
C  Purpose:  This subroutine determines the node-based nonlinear 
C            horizontal diffusion matrix (AGPGP) and node-based values
C            of horizontal viscosity/diffusivity (AHI).  The formulation
C            is similar in spirit to Smagorinsky (1963) (i.e. horizontal
C            viscosity/diffusivity is parameterized in terms of the 
C            local shear and grid scale).  The exact Laplacian form 
C            implemented herein is as presented in Lynch et al (1996) 
C            [CSR, pp 875-906].
C
C  Notes: This version differs from the 19 JUN 98 version of SUBROUTINE
C         SMAGOR in that there is a depth adjustment added to the 
C         integration rule as in Mellor et al 198X.
C
C         Original Version of SMAGOR: 
C           Daniel R. Lynch and Justin T.C. Ip,                 3 AUG 93
C         Modified Version of SMAGOR: Christopher E. Naimie,   19 JUN 98
C         Original Version of SMAGOR1:
C           Daniel R. Lynch and Justin T.C. Ip,                 1 OCT 98
C          KWS- Returned coeficient - DP for bellamy drouge tracking
C-----------------------------------------------------------------------
C
C Global arrays
C
      implicit none
      include 'DCMS.DIM'
      REAL U(*),V(*),AR(*),X(*),Y(*),AHI(*),SMAGR,AHMIN
      INTEGER IN(3,*)
      REAL SV(NNDIM)
C
C Local dimensioning and assignments
C
      REAL DELY(3),DELX(3),HAVG,D1,D2,AH,ARR
      INTEGER I,L,i1,i2,i3,ii,NN,NE
C
C Initialize AGPGP and AHI
C
      DO I=1,NN
         AHI(I)=0.D0
      ENDDO
C
C  Laplacian matrix:
C     AGPGP = -<AH*grad(phi_i) . grad(phi_j)> 
C     AGPGP*V = + <DEL.(AH GRAD V), phi_i>
C
      DO L=1,NE       
         ARR=AR(L) 
         I1=IN(1,L)
         I2=IN(2,L)
         I3=IN(3,L)
         DELX(1)=X(I2)-X(I3)          
         DELX(2)=X(I3)-X(I1)          
         DELX(3)=X(I1)-X(I2)          
         DELY(1)=Y(I2)-Y(I3)          
         DELY(2)=Y(I3)-Y(I1)          
         DELY(3)=Y(I1)-Y(I2)          
         D1=(U(I1)*DELY(1)+U(I2)*DELY(2)+U(I3)*DELY(3)  
     &      +V(I1)*DELX(1)+V(I2)*DELX(2)+V(I3)*DELX(3))/2./ARR    
         D2=(V(I1)*DELY(1)+V(I2)*DELY(2)
     &      +V(I3)*DELY(3)  
     &      -U(I1)*DELX(1)-U(I2)*DELX(2)-U(I3)*DELX(3))/2./ARR
C
C  Notes: 2.*ARR is a measure of the local grid scale squared
C         SQRT(D1**2+D2**2) is a measure of the local shear
C
         AH=SMAGR*2.*ARR*SQRT(D1**2+D2**2)    
         DO I=1,3
            II=IN(I,L)
            AHI(II)=AHI(II)+AH*ARR/3.
         ENDDO
      ENDDO

      DO I=1,NN
       SV(I)=0.
      ENDDO
      DO L=1,NE       
         ARR=AR(L) 
         DO I=1,3
            II=IN(I,L)
            SV(II)=SV(II)+ARR/3.
         ENDDO
      ENDDO
C
C Divide by mass matrix and nodal total depth to compute AHI and compare
C to minimum value
C
      DO I=1,NN
         AHI(I)=AHI(I)/SV(I)
         AHI(I)=MAX(AHI(I),AHMIN)
      ENDDO
C
      RETURN
      END
C***********************************************************************
C
C******************* END OF EVERYTHING *************************C
C******************* END OF EVERYTHING *************************C
C******************* END OF EVERYTHING *************************C


