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 *** CURRICULUM IN MARINE SCIENCE C *** 12-7 VENABLE HALL C *** CB #3300 C *** CHAPEL HILL, NORTH CAROLINA C *** 27599-3300 C *** C *** MARCH 1995 C C #$%#$%#$% THIS IS A GAMMA-TESTING VERSION FOR 3-DIMENSIONAL TIME- C #$%#$%#$% DEPENDENT DROGUE TRACKING. 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 EXTRAPOLATION TO 3D OF HARMONIC CODE : C BRIAN O. BLANTON, FALL 92, SKIO C C CONVERSION FROM HARMONIC TO TIME-DEPENDENT VELOCITY INPUT: C BRIAN O. BLANTON, SUMMER 94, UNC-CH C C------------------------------------------------------------------------ C SUBROUTINE DROG3DDT(NMND,NMEL,NNV,T2HOURSQ2,T1HOURSQ2,ITQ, + DEPTH2,UT2,VT2,WT2,DOMAINNAME,INDNAME) C C----------------------------------------------------------------------- C C SET MAXIMUM DIMENSIONING OF ARRAYS THROUGH EXTERNAL INCLUDE C INCLUDE 'QUODDY.DIM' 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=1000) C C----------------------------------------------------------------------- C C *** DECLARE VARIABLES INCOMING 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 T2HOURSQ2- TIME AT END OF ITQ TIMESTEP, IN HOURS C T1HOURSQ2- TIME AT START OF ITQ TIMESTEP, IN HOURS C ITQ - TIME-STEPPER ITERATION NUMBER C DEPTH2 - ARRAY CONTAINING VERTICAL GRID POSITIONS AT TIME T2HOURSQ2 C UT2 - U-VELOCITY COMPONENT AT TIME T2HOURSQ2; UZMID IN QUODDY2 C VT2 - V-VELOCITY COMPONENT AT TIME T2HOURSQ2; VZMID IN QUODDY2 C WT2 - W-VELOCITY COMPONENT AT TIME T2HOURSQ2; WZMID IN QUODDY2 C DOMAINNAME - NAME OF DOMAIN (GRID); CASNAM IN QUODDY2 C INDNAME - NAME OF DROG INPUT FILE; FILDR3 IN QUODDY2 C INTEGER NMND,NMEL,NNV REAL T2HOURSQ2,T1HOURSQ2 INTEGER ITQ REAL DEPTH2(1:NNDIM,*) REAL UT2(1:NNDIM,*) REAL VT2(1:NNDIM,*) REAL WT2(1:NNDIM,*) CHARACTER*(*) DOMAINNAME,INDNAME C C *** T1SECQ2 - MODEL TIME AT BEGINNING OF ITQ TIMESTEP, IN SECONDS C *** T2SECQ2 - MODEL TIME AT END OF ITQ TIMESTEP, IN SECONDS C *** DEPTH1 - ARRAY CONTAINING VERTICAL GRID POSITIONS AT TIME T1HOURSQ2 C *** UT1 - U-VELOCITY COMPONENT AT TIME T1HOURSQ2 C *** VT1 - V-VELOCITY COMPONENT AT TIME T1HOURSQ2 C *** WT1 - W-VELOCITY COMPONENT AT TIME T1HOURSQ2 C REAL T1SECQ2,T2SECQ2 REAL DEPTH1(1:NNDIM,1:NNVDIM) REAL UT1(1:NNDIM,1:NNVDIM) REAL VT1(1:NNDIM,1:NNVDIM) REAL WT1(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 *** SCALAR DROGUE INFO PASSED TO TRACK C *** BEFORE CALLING TRACK, DROG3DDT MOVES INFORMATION ABOUT THE DROGUE C *** CURRENTLY BEING TRACKED INTO SCALAR VARIABLES C REAL U0,V0,W0 REAL X0,Y0,Z0 INTEGER J0,L0,NDR,JNEW,IBOUN REAL XNEW,YNEW,ZNEW INTEGER LNEW DATA IBOUN/0/ C C *** MODEL TIME-STEP, DROGUE TRACKING TIME-STEP C REAL T1SECDR,T2SECDR 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 C C *** SAVE CERTAIN VARIABLES TO AVOID -K COMPILE-LINE C *** OPTION ON HP9000 C SAVE UT1,VT1,WT1,DEPTH1,T1SECDR,T2SECDR,T1SECQ2,T2SECQ2 SAVE XDR,YDR,ZDR SAVE IDR,JJDR,LLDR,NDR SAVE ITRACK,IPRINT,NACTIVE SAVE DT_TS,DT_DR C C *** CONVERT INPUT TIMES TO SECONDS; T1SECQ2 TO T2SECQ2 IS THE MODEL C *** INTERVAL; THE TRACKING INTERVAL SHOULD BE T1SECQ2 TO T1SECQ2+DT_DR C *** OR T1SECQ2+DT_TS*ITRACK C T1SECQ2=T1HOURSQ2*3600. T2SECQ2=T2HOURSQ2*3600. C C CALL INPUT AND INITIALIZATION ROUTINE BEFORE INITIAL TIMESTEP ONLY C I.E., IT=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 'ITQ' PLACE. THIS IS THE VARIABLE ITER IN OUTPUTQ2. C IF(ITQ.EQ.0)THEN TRACKFLAG='YES' WRITE(2,*)'DROGUE TRACKING STARTED AT T = ',HOURS DO I=1,NMND DO L=1,NNV UT1(I,L)=UT2(I,L) VT1(I,L)=VT2(I,L) WT1(I,L)=WT2(I,L) DEPTH1(I,L)=DEPTH2(I,L) END DO END DO C C *** COMPUTE MODEL TIME-STEP C DT_TS=T2SECQ2-T1SECQ2 T1SECDR=T2SECQ2 CALL INPUT(NMND,NMEL,NNV,MAXDROG,T1SECQ2,T2SECQ2,DT_TS,DT_DR, + IPRINT,ITRACK,NDR,DEPTH1,DEPTH2, + IDR,JJDR,LLDR,XDR,YDR,ZDR, + DOMAINNAME,INDNAME) IF(NDR.EQ.0)RETURN NACTIVE=NDR 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 'ITQ' IS A C *** MULTIPLE OF ITRACK TO TRACK C T2SECDR=T2SECQ2 IF(TRACKFLAG.EQ.'NO')RETURN IF(MOD(ITQ,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 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 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 T1SECDR TO T2SECDR C CALL RK4(NNV,T1SECDR,T2SECDR,DEPTH1,DEPTH2, + UT1,UT2,VT1,VT2,WT1,WT2, + J0,JNEW,L0,LNEW,X0,Y0,Z0,U0,V0,W0,DT_DR, + XNEW,YNEW,ZNEW,IBOUN) C C *** UPDATE THE PARTICLE'S NEW POSITION C JJDR(II) = JNEW LLDR(II) = LNEW XDR(II) = XNEW YDR(II) = YNEW ZDR(II) = ZNEW C C *** IF PARTICLE HITS A BOUNDARY, RESTORE PREVIOUS POSITION C *** AND ELIMINATE IT C IF (IBOUN.EQ.1) THEN C C *** SET IBOUN TO 0 FOR NEXT DROGUE; I.E., DROGUE WITHIN DOMAIN C IBOUN = 0 JJDR(II) = J0 LLDR(II) = L0 XDR(II) = X0 YDR(II) = Y0 ZDR(II) = Z0 IDR(II) = 0 NACTIVE=NACTIVE-1 WRITE(13,8010)II,ITQ,T1SECDR/86400.E0 WRITE(*,8010)II,ITQ,T1SECDR/86400.E0 END IF 160 CONTINUE C C *** CALL SUBROUTINE OUTPUT C 165 IF(MOD(ITQ,IPRINT*ITRACK).EQ.0) CALL OUTPUT(NNV,T1SECDR,T2SECDR, + DEPTH1,DEPTH2,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 UT1(I,L)=UT2(I,L) VT1(I,L)=VT2(I,L) WT1(I,L)=WT2(I,L) DEPTH1(I,L)=DEPTH2(I,L) END DO END DO T1SECDR=T2SECDR 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 'QUODDY.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) 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 QUODDY.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)'QUODDY.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,J),J=1,3) 20 CONTINUE C C *** ELEMENT-ELEMENT ADJACENCY INFORMATION C DO 30 I = 1,NMEL READ(9,*) N, (ICEE(N,J),J=1,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 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 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) CALL GET_DEPTH(DEPTHNOW,JJDR(II),XDR(II),YDR(II),1,ZBOT) WRITE(12,9005) XDR(II),YDR(II),ZDR(II),ZBOT 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 'QUODDY.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 'QUODDY.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,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 'QUODDY.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 'QUODDY.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 'QUODDY.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 'QUODDY.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 'QUODDY.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 'QUODDY.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 'QUODDY.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,DEPTH1,DEPTH2, + UT1,UT2,VT1,VT2,WT1,WT2, + J,JOUT,L,LOUT,XX,YY,ZZ,U,V,W,DT, + XOUT,YOUT,ZOUT,IBOUN) C*********************************************************************** C C----------------------------------------------------------------------- INCLUDE 'QUODDY.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 'QUODDY.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 'QUODDY.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*********************************************************************** 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 'QUODDY.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 'QUODDY.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 'QUODDY.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 'QUODDY.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,DEPTH1,DEPTH2, + IDR,JJDR,LLDR,NDR,XDR,YDR,ZDR) C************************************************************************** C INCLUDE 'QUODDY.DIM' C C *** NNV - NUMBER OF VERTICAL NODES C *** T1,T2 - BEGINNING AND ENDING TIME OF THIS TIME-STEP IN SEC C *** DEPTH1,DEPTH2 - 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 DEPTH1(1:NNDIM,*),DEPTH2(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 DEPTHNOW(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,DEPTH1,DEPTH2,JJDR(I),T2,DEPTHNOW) CALL GET_DEPTH(DEPTHNOW,JJDR(I),XDR(I),YDR(I),1,ZBOT) WRITE(12,9000) XDR(I),YDR(I),ZDR(I),ZBOT WRITE(17,9000) XDR(I),YDR(I),ZDR(I) END DO C C *** CLOSE THE DROGUE HOT START FILE C CLOSE(UNIT=17) RETURN 9000 FORMAT (5(2X,F12.4),I4) END C******************* END OF EVERYTHING *************************C C******************* END OF EVERYTHING *************************C C******************* END OF EVERYTHING *************************C