Date: Fri, 24 May 91 12:36:19 +0200 From: jarleb@sentral.imr.no (Jarle Berntsen) Message-Id: <9105241036.AA15320@imr.no> To: miller@rcf.rsmas.miami.edu Status: RO Jerry, Thank you for your answer. About tracing. I may have differnt collections of tracers, and for each internal step I compute new positions by an ordinary Euler forward step. The difficult part (if there is one) is to take care that all tracers remain in the sea and to keep track of which tracers leave our region. I include my routines at the end of this letters. TRACE is really doing the job and is best documented. TRACSU just calls TRACE. I also send you my RIVER code. The river runoffs are read from file and may be a function of time. I do not know what ECOM3D comes from, but this is the start of the main program from Blumberg. PROGRAM OPEC1 C*********************************************************************** C : C Y C GENERAL CIRCULATION MODEL Y C ORTHOGONAL CURVILINEAR COORDINATE VERSION Y C ECOM3D Y C Y C made by Alan F. Blumberg. ( HydroQual Inc.) : C Y C : C PREOPERATIONAL VERSION Jarle SUBROUTINE TRACE(TRACER,NUMTRA,U,V,W,IMAX,JMAX,KMAX, + DT,DX,DY,H,E,DZ, + FILENR,IOUT,IPER,IOPT,IFAIL) C***BEGIN PROLOGUE TRACE C***DATE WRITTEN 910107 (YYMMDD) C***REVISION DATE 910213 (YYMMDD) C***AUTHOR C Jarle Berntsen, Institute of Marine Research, C Postboks 1870, C N-5024 Bergen-Nordnes, Norway C Email.. jarleb@imr.no C C***PURPOSE TRACE writes the coordinates of TRACER to file, C or propagates the tracers and writes the coordinates C to file if MOD(IOUT,IPER) = 0. C C This version of TRACE is based on the C datastructures used in the ECOM3D model. C TRACE assumes constant grid spacing. C That is DX and DY must both be constant. C C***DESCRIPTION C C The coordinates are defined in a 3-D rectangular region C defined C from 0 to IMAX in X-direction, C from 0 to JMAX in Y-direction and C from -H to the water level in Z-direction. C C The coordinates in Z-direction are transformed to C the sigma coordinate system through the transformation C SIGMA = (Z - ETA)/(H + ETA) C before the tracers are propagated. C After the propagation we transform back to the C physical coordinate system. C C If IOPT = 1, TRACE will write the coordinates of the C NUMTRA tracers in TRACER to a file, opened in the C calling program, with unit number FILENR. C C If IOPT = 2, TRACE will propagate the NUMTRA tracers C by computing the new positions as functions of U, V and W. C After the tracers have been propagated, IOPT is incremented C by 1 and the coordinates of each tracer is written to file C if MOD(IOPT,IPER) = 0. C C If a tracer is outside the region of integration on C entry and IOPT = 2, we exit with IFAIL = 2. C C If a tracer is on land on entry and IOPT=2, C we exit with IFAIL = 3. C C If a tracer is in sea on entry and on land C on exit (IOPT=2), we exit with IFAIL = 4. C This should in principle never happen, and the C author should be contacted if we exit from TRACE C with IFAIL = 4. C C If a tracer moves out of the region of integration, C we mark this tracer as dead by giving it the C coordinates RUDEF. C This tracer will not be further propagated. C C If the velocities U, V and W indicate that the tracer C should move from a sea cell to a land cell, C the tracer is moved in the indicated direction C as far as possible INSIDE the sea cell. C C ON ENTRY C C TRACER Real array of dimension (3,NUMTRA). C The coordinates of NUMTRA tracers. C IF IOPT = 1, TRACER will be unchanged on exit. C IF IOPT = 2, the tracers will be moved and C the coordinates may be changed. C The leading dimension of TRACER MUST be declared C 3 in the calling program. C NUMTRA Integer. C The number of tracers. C U Real array of dimension (IMAX,JMAX,KMAX). C Velocity in X-direction defined in U points (m/s). C The leading two dimensions of U MUST be declared C IMAX+1,JMAX in the calling program. C V Real array of dimension (IMAX,JMAX,KMAX). C Velocity in Y-direction defined in V points (m/s). C The leading two dimensions of V MUST be declared C IMAX,JMAX+1 in the calling program. C W Real array of dimension (IMAX,JMAX,KMAX). C Velocity in Z-direction defined in W points (m/s). C The leading two dimensions of W MUST be declared C IMAX,JMAX in the calling program. C IMAX Integer. C Number of cells in X-direction. C JMAX Integer. C Number of cells in Y-direction. C KMAX Integer. C Number of cells in Z-direction. C DT Real. C The timestep (s). C DX Real. C The width of a cell in X-direction (m). C DY Real. C The width of a cell in Y-direction (m). C H Real array of dimension (IMAX,JMAX). C The depth matrix of the region. C E Real array of dimension (IMAX,JMAX). C The water elevation. C DZ Real array of dimension KMAX. C The thickness of each sigma layer. C FILENR Integer. C Unit number of file. C IOUT Integer. C IOUT is increased by 1 each time the tracers C are propagated, and the new coordinates are C written to file if MOD(IOPT,IPER) = 0. C IPER Integer. C IPER defines how frequently the positions C of the tracers shall be written to file. C IOPT Integer. C If IOPT = 1, TRACE will write the coordinates of the C NUMTRA tracers in TRACER to the file with unit number C FILENR. C If IOPT = 2, TRACE will propagate the NUMTRA tracers C by computing the new positions as functions of U, V and W. C C ON RETURN C C TRACER Real array of dimension (3,NUMTRA). C The coordinates of NUMTRA tracers. C IF IOPT = 2, TRACER will contain C the coordinates of the tracers after C they have been moved a timestep DT. C IOUT Integer. C If IOPT = 2 and the tracers are successfully C propagated, IOUT is incremented by 1 on return. C IFAIL Integer. C IFAIL is 0 on exit if no errors C has been detected. C IFAIL is 1 on exit if C NUMTRA <= 0 or C IMAX <= 0 or C JMAX <= 0 or C KMAX <= 0 or C FILENR <= 6 or C IPER = 0 and IOPT = 2 or C IOPT < 1 or C IOPT > 2. C IFAIL is 2 on exit if C a tracer is outside the region on entry. C IFAIL is 3 on exit if C a tracer is on land on entry. C IFAIL is 4 on exit if C a tracer is in sea on entry and on land C on exit. C C C***ROUTINES CALLED-NONE C***END PROLOGUE TRACE C C Global variables. C INTEGER NUMTRA,IOUT,IPER,FILENR,IMAX,JMAX,KMAX,IOPT,IFAIL REAL TRACER(3,NUMTRA),U(IMAX,JMAX,KMAX), + V(IMAX,JMAX,KMAX),W(IMAX,JMAX,KMAX), + DT,DX,DY,E(IMAX,JMAX),DZ(KMAX),H(IMAX,JMAX) C C Local variables. C C DEPTHP Depth of bottom of actual cell. C PARTX Relative position in X-direction in present cell. C PARTY Relative position in Y-direction in present cell. C PARTZ Relative position in Z-direction in present cell. C UPOS Estimated value of U at this position. C VPOS Estimated value of V at this position. C WPOS Estimated value of W at this position. C XPOS Esitmated new X-coordinate of tracer. C YPOS Esitmated new Y-coordinate of tracer. C SPOS Esitmated new Z-coordinate of tracer. C ICELL Index in X-direction of the west side of present cell. C JCELL Index in Y-direction of the south side of present cell. C KCELL Index in Z-direction of the top of present cell. C XDIST Distance in X-direction between new and old position. C YDIST Distance in Y-direction between new and old position. C DISTX Distance in X-direction between old position and nearest C X-axis. C DISTY Distance in Y-direction between old position and nearest C Y-axis. C RUDEF Undefined REAL value. C PI The mathematical constant PI. C ALPHA Computed direction of the tracer. (0 to 2*PI). C REAL DEPTHP,PARTX,PARTY,PARTZ,UPOS,VPOS,WPOS,XPOS,YPOS,SPOS INTEGER I,ITRACE,ICELL,JCELL,KCELL,IXNEW,JYNEW REAL RUDEF,DISTX,DISTY,XDIST,YDIST,PI,ALPHA,EPSIL,SIGMA PARAMETER (RUDEF = 1.E12, EPSIL = 0.001) C C***FIRST EXECUTABLE STATEMENT TRACE C C C We first check on legal input. C IF (NUMTRA.LE.0 .OR. + IMAX.LE.0 .OR. + JMAX.LE.0 .OR. + KMAX.LE.0 .OR. + FILENR.LE.6 .OR. + (IPER.EQ.0 .AND. IOPT.EQ.2) .OR. + IOPT.LT.1 .OR. + IOPT.GT.2 ) THEN IFAIL = 1 RETURN ENDIF C C Define PI. C PI = 4.*ATAN(1.0) C C If IOPT = 1, we write the current positions to file C and exit. C IF (IOPT.EQ.1) THEN DO 10 ITRACE = 1,NUMTRA DO 10 I = 1,3 WRITE(FILENR,*)TRACER(I,ITRACE) 10 CONTINUE IFAIL = 0 RETURN ENDIF C C If IOPT = 2, we propagate the tracers and write the C new positions to file if MOD(IOUT,IPER) = 0. C IF (IOPT.EQ.2) THEN C C Main loop for computing new positions for C particles. C DO 500 ITRACE = 1,NUMTRA C C If the tracer is dead, we propagate the next tracer. C IF (TRACER(1,ITRACE).EQ.RUDEF) THEN GO TO 500 ENDIF C C Compute actual box. C ICELL = INT(TRACER(1,ITRACE)) + 1 JCELL = INT(TRACER(2,ITRACE)) + 1 IF (TRACER(1,ITRACE).EQ.FLOAT(IMAX)) THEN ICELL = IMAX ENDIF IF (TRACER(2,ITRACE).EQ.FLOAT(JMAX)) THEN JCELL = JMAX ENDIF C C If the tracer is outside the actual region, C we return with IFAIL = 2. C IF (ICELL.LT.1 .OR. ICELL.GT.IMAX .OR. + JCELL.LT.1 .OR. JCELL.GT.JMAX) THEN IFAIL = 2 RETURN ENDIF SIGMA = (TRACER(3,ITRACE) - E(ICELL,JCELL))/ + (H(ICELL,JCELL) + E(ICELL,JCELL)) C C The tracer may have moved into the air (SIGMA > 0.), C if the water level has fallen since the last call to C TRACE. In this case the tracer is moved into the top C level of the sea again. C IF (SIGMA.GT.0.) THEN SIGMA = -EPSIL TRACER(3,ITRACE) = SIGMA* + (H(ICELL,JCELL) + E(ICELL,JCELL)) + E(ICELL,JCELL) ENDIF IF ((H(ICELL,JCELL).LE.0.).OR. + (SIGMA.LT.-1.)) THEN IFAIL = 3 WRITE(*,*)ICELL,JCELL WRITE(*,*)H(ICELL,JCELL),SIGMA WRITE(*,*)TRACER(3,ITRACE),E(ICELL,JCELL) RETURN ENDIF KCELL = 1 DEPTHP = -DZ(KCELL) 50 IF (KCELL.LT.KMAX-1) THEN IF (SIGMA.LT.DEPTHP) THEN KCELL = KCELL + 1 DEPTHP = DEPTHP - DZ(KCELL) GO TO 50 ENDIF ENDIF C C Compute relative position in actual box. C PARTX = TRACER(1,ITRACE) + 1 - ICELL PARTY = TRACER(2,ITRACE) + 1 - JCELL PARTZ = (SIGMA - DEPTHP + DZ(KCELL))/DZ(KCELL) C C Compute new position by linear interpolation and Euler forward. C IF (ICELL.LT.IMAX) THEN UPOS = (1. -PARTX)*U(ICELL,JCELL,KCELL) + + PARTX*U(ICELL+1,JCELL,KCELL) ELSE UPOS = U(ICELL,JCELL,KCELL) ENDIF XPOS = TRACER(1,ITRACE) + DT*UPOS/DX IF (JCELL.LT.JMAX) THEN VPOS = (1. -PARTY)*V(ICELL,JCELL,KCELL) + + PARTY*V(ICELL,JCELL+1,KCELL) ELSE VPOS = V(ICELL,JCELL,KCELL) ENDIF YPOS = TRACER(2,ITRACE) + DT*VPOS/DY WPOS = (1. -PARTZ)*W(ICELL,JCELL,KCELL) + + PARTZ*W(ICELL,JCELL,KCELL+1) SPOS = SIGMA + DT*WPOS SPOS = MAX(SPOS,-1.0+EPSIL) SPOS = MIN(SPOS,-EPSIL) C C Define new cell. C IXNEW = INT(XPOS) + 1 JYNEW = INT(YPOS) + 1 IF (XPOS.EQ.FLOAT(IMAX)) THEN IXNEW = IMAX ENDIF IF (YPOS.EQ.FLOAT(JMAX)) THEN JYNEW = JMAX ENDIF C C If the tracer moves out of our region, we mark the C tracer as dead by giving the coordinates the values C RUDEF. C IF (IXNEW.LT.1 .OR. IXNEW.GT.IMAX .OR. + JYNEW.LT.1 .OR. JYNEW.GT.JMAX) THEN TRACER(1,ITRACE) = RUDEF TRACER(2,ITRACE) = RUDEF TRACER(3,ITRACE) = RUDEF GO TO 500 ENDIF C C If we have entered a land cell, C we move the tracer in the old cell as far as C possible in the computed direction. C IF (H(IXNEW,JYNEW).LE.0.) THEN XDIST = XPOS - TRACER(1,ITRACE) YDIST = YPOS - TRACER(2,ITRACE) C C Compute the direction of the tracer. C IF (XDIST.EQ.0.) THEN IF (YDIST.EQ.0.) THEN ALPHA = 0. ELSEIF (YDIST.LT.0.) THEN ALPHA = 3.*PI/2. ELSEIF (YDIST.GT.0) THEN ALPHA = PI/2. ENDIF ELSEIF (YDIST.EQ.0.) THEN IF (XDIST.LT.0.) THEN ALPHA = PI ELSEIF (XDIST.GT.0.) THEN ALPHA = 0. ENDIF ELSE ALPHA = ATAN2(YDIST,XDIST) ENDIF C C Move ALPHA to (0,2*PI) C IF (ALPHA.LT.0.) THEN ALPHA = ALPHA +2*PI ENDIF C C Then follow this direction almost to the edge C of the cell. C IF (ALPHA.EQ.0.) THEN XPOS = FLOAT(ICELL) - EPSIL YPOS = TRACER(2,ITRACE) ELSEIF (ALPHA.GT.0. .AND. + ALPHA.LT.PI/2) THEN DISTY = (JCELL - TRACER(2,ITRACE))/SIN(ALPHA) DISTX = (ICELL - TRACER(1,ITRACE))/COS(ALPHA) IF (DISTY.LE.DISTX) THEN YPOS = FLOAT(JCELL) - EPSIL XPOS = TRACER(1,ITRACE) + DISTY*COS(ALPHA) ELSE XPOS = FLOAT(ICELL) - EPSIL YPOS = TRACER(2,ITRACE) + DISTX*SIN(ALPHA) ENDIF ELSEIF (ALPHA.EQ.PI/2) THEN XPOS = TRACER(1,ITRACE) YPOS = FLOAT(JCELL) - EPSIL ELSEIF (ALPHA.GT.PI/2 .AND. + ALPHA.LT.PI) THEN DISTY = (JCELL - TRACER(2,ITRACE))/SIN(PI-ALPHA) DISTX = (TRACER(1,ITRACE) - (ICELL-1))/COS(PI-ALPHA) IF (DISTY.LE.DISTX) THEN YPOS = FLOAT(JCELL) - EPSIL XPOS = TRACER(1,ITRACE) - DISTY*COS(PI-ALPHA) ELSE XPOS = FLOAT(ICELL) - (1.0 - EPSIL) YPOS = TRACER(2,ITRACE) + DISTX*SIN(PI-ALPHA) ENDIF ELSEIF (ALPHA.EQ.PI) THEN XPOS = FLOAT(ICELL) - (1.0 - EPSIL) YPOS = TRACER(2,ITRACE) ELSEIF (ALPHA.GT.PI .AND. + ALPHA.LT.3*PI/2) THEN DISTY = (TRACER(2,ITRACE) - (JCELL-1))/SIN(ALPHA-PI) DISTX = (TRACER(1,ITRACE) - (ICELL-1))/COS(ALPHA-PI) IF (DISTY.LE.DISTX) THEN YPOS = FLOAT(JCELL) - (1.0 - EPSIL) XPOS = TRACER(1,ITRACE) - DISTY*COS(ALPHA-PI) ELSE XPOS = FLOAT(ICELL) - (1.0 - EPSIL) YPOS = TRACER(2,ITRACE) - DISTX*SIN(ALPHA-PI) ENDIF ELSEIF (ALPHA.EQ.3*PI/2) THEN YPOS = FLOAT(JCELL) - (1.0 - EPSIL) XPOS = TRACER(1,ITRACE) ELSEIF (ALPHA.GT.3*PI/2 .AND. + ALPHA.LT.2*PI) THEN DISTY = (TRACER(2,ITRACE) - (JCELL-1))/SIN(2*PI-ALPHA) DISTX = (ICELL - TRACER(1,ITRACE))/COS(2*PI-ALPHA) IF (DISTY.LE.DISTX) THEN YPOS = FLOAT(JCELL) - (1.0 - EPSIL) XPOS = TRACER(1,ITRACE) + DISTY*COS(2*PI-ALPHA) ELSE XPOS = FLOAT(ICELL) - EPSIL YPOS = TRACER(2,ITRACE) - DISTX*SIN(2*PI-ALPHA) ENDIF ENDIF C C Redefine IXNEW and JYNEW. C IXNEW = INT(XPOS) + 1 JYNEW = INT(YPOS) + 1 IF (XPOS.EQ.FLOAT(IMAX)) THEN IXNEW = IMAX ENDIF IF (YPOS.EQ.FLOAT(JMAX)) THEN JYNEW = JMAX ENDIF C C End if we have entered a land cell. C ENDIF C C We test on whether we still are in the sea. C IF (H(IXNEW,JYNEW).LE.0.) THEN IFAIL = 4 RETURN ENDIF C C The horisontal coordinates will now be in a sea cell. C C C The tracer is now placed inside a sea cell, C and we update the position. C TRACER(1,ITRACE) = XPOS TRACER(2,ITRACE) = YPOS TRACER(3,ITRACE) = SPOS*(H(IXNEW,JYNEW)+E(IXNEW,JYNEW)) + + E(IXNEW,JYNEW) C C End main loop. C 500 CONTINUE C C We increment IOUT, and test on printout. C IOUT = IOUT + 1 IF (MOD(IOUT,IPER).EQ.0) THEN DO 600 ITRACE = 1,NUMTRA DO 600 I = 1,3 WRITE(FILENR,*)TRACER(I,ITRACE) 600 CONTINUE ENDIF C C END IF IOPT=2. C ENDIF IFAIL = 0 C RETURN END subroutine tracsu(trace1,numtr1, + maxtra,u,v,w,imax,jmax,kmax, + dt,dx,dy,H,E,DZ,iopt) c integer imax,jmax,kmax real u(imax,jmax,kmax),v(imax,jmax,kmax) real w(imax,jmax,kmax) real dt,dx,dy real H(IMAX,JMAX),E(IMAX,JMAX),DZ(KMAX) integer maxtra real trace1(3,maxtra) integer numtr1 c integer iopt,iout1,iper,ifail logical first save first,iout1,iper iopt = 2 data first/.TRUE./ if (first) then iper = 12 open (61,file='trace1.dat') iout1 = 0 first = .FALSE. numtr1 = 0 do 100 i=1,imax,2 do 100 j = 1,jmax,2 if (h(i,j).gt.10.0) then numtr1 = numtr1 + 1 trace1(1,numtr1) = float(i) - 0.5 trace1(2,numtr1) = float(j) - 0.5 trace1(3,numtr1) = -5.0 endif 100 continue iopt = 1 write(*,*)'tracsu-numtri1=',numtr1 endif c call trace(trace1,numtr1,u,v,w,imax,jmax,kmax, + dt,dx,dy,H,E,DZ, + 61,iout1,iper,iopt,ifail) if (ifail.gt.0) then write(*,*)'after trace1' write(*,*)'ifail =',ifail write(*,*)'iout1=',iout1,iopt stop endif c return end