C*********************************************************************** C pom.f for the North Eastern Atlantic Ocean, especially C forcasting and data assimilation for the east coast of U.S.A.. C This code contains two versions of S.R. BCOND. Check C the subroutines for details. C --------------------------------------------------------- C Keep record of subroutines that are "commented out". C none C --------------------------------------------------------- C calculation of the monthly mean of T,S,E,PSI,U and V. C 08/04/97 Lee C consider only the surface heat flux without SSP C---------------------------------------------------------- PROGRAM MAIN C************************************************************************ C TEST WITH SIMPLE FLAT BOTTOMED BASIN * C * C THIS IS A VERSION OF THE THREE DIMENSIONAL, TIME DEPENDENT, * C PRIMITIVE EQUATION, OCEAN MODEL DEVELOPED BY ALLAN BLUMBERG AND ME * C WITH SUBSEQUENT CONTRIBUTIONS BY LEO OEY AND OTHERS. TWO REFERENCES * C ARE: * C * C BLUMBERG, A.F. AND G.L. MELLOR, DIAGNOSTIC AND PROGNOSTIC * C NUMERICAL CIRCULATION STUDIES OF THE SOUTH ATLANTIC BIGHT * C J. GEOPHYS. RES. 88, 4579-4592, 1983. * C * C BLUMBERG, A.F. AND G.L. MELLOR, A DESCRIPTION OF A THREE * C COASTAL OCEAN CIRCULATION MODEL, THREE DIMENSIONAL SHELF * C MODELS, COASTAL AND ESTUARINE SCIENCES, 5, N.HEAPS, ED., * C AMERICAN GEOPHYSICAL UNION, 1987. * C * C IN SUBROUTINE PROFQ THE MODEL MAKES USE OF THE TURBULENCE CLOSURE * C SUB-MODEL MOST RECENTLY DESCRIBED IN: * C * C MELLOR, G.L. AND T. YAMADA, DEVELOPMENT OF A TURBULENCE CLOSURE * C MODEL FOR GEOPHYSICAL FLUID PROBLEMS, REV. GEOPHYS. SPACE PHYS., * C 851-879, 1982. * C * C GEORGE MELLOR, 03/87 * C * C IN THIS VERSION THE HORIZONTAL COORDINATES ARE GENERALIZED CURVI- * C LINEAR ORTHOGANAL COORDINATES AND THE MODEL HAS BEEN CONFIGURED * C TO RUN IN HALF PRECISION (32 BITS ) ON THE CYBER 205. * C TO REDUCE ROUNDOFF ERROR, SALINTIES AND TEMPERATURES ARE REDUCED * C BY 35. AND 10. RESPECTIVELY. 64 BIT PRECISION IS USED LOCALLY * C IN SUBROUTINES DENS (WHEREIN THE T AND S OFFSET HAS BEEN TAKEN * C INTO ACCOUNT) AND PROFT. * C GM, 09/87 * C * C THIS VERSION HAS BEEN CONVERTED TO STANDARD FORTRAN 77 BY STEVE * C BRENNER WITH SOME ADDITIONAL CHANGES BY ME IN PROFQ WHICH IS NOW * C SOMEWHAT SIMPLER. CODE RESTORED TO SINGLE PRECISION. * C A USER'S GUIDE IS AVAILABLE: * C * C MELLOR, G.L., A USER'S GUIDE FOR A THREE-DIMENSIONAL, C ERICAL OCEAN MODEL. PRINCETON UNIVERSITY REPORT, 1990. * C GM, 04/90 * C * C Previous to this date the code was called pmod.f and all changes * C were recorded in pmod.change. Henceforth, the code will be called * C pom.f and changes will be recorded in pom.change. * C GM, 11/92 * C * C THE LAST CODE CHANGE, as recorded in pom.change, WAS ON March 5, 1997 * C * C * C Copyright 1996 The Trustees of Princeton University * C * C Princeton makes no warranties, expressed or implied, as to any matter* C whatsoever, including, without limitation, the condition of the POM * C software, including, but not limited to, ownership, merchantability, * C or fitness for a particular purpose. Princeton shall not be liable * C for any direct, consequential, or other damages suffered by any user * C or any others resulting from the use of this research software. * C * C This code was converted to the barotropic experiment. C H.C.Lee 06/97 C************************************************************************ C INCLUDE 'cof_a.h' LOGICAL SEAMT DIMENSION UTB(IM,JM),VTB(IM,JM),UTF(IM,JM),VTF(IM,JM), 1 ADVUA(IM,JM),ADVVA(IM,JM), 2 TSURF(IM,JM),SSURF(IM,JM), 2 DRX2D(IM,JM),DRY2D(IM,JM), 3 ADX2D(IM,JM),ADY2D(IM,JM), 4 SWRAD(IM,JM),X(IM),Y(JM), 7 ENDM(12),sem(jm) C-------------------------------------------------------------------- C++++++++++++++LHC REAL hh,dd,mm,yy, h_rs,d_rs,m_rs,y_rs,h_b,d_b,m_b,y_b & ,h_s,d_s,m_s,y_s C++++++++++++++LHC REAL ISPI,ISP2I DATA PI/3.141592654/,SMALL/1.E-10/ DATA (ENDM(M),m=1,12)/ & 31.,28.,31.,30.,31.,30.,31.,31.,30.,31.,30.,31./ GRAV=9.806 TBIAS=0. SBIAS=0. IINT=0 TIME=0. C--------overwrited by params--------------------------------------- C DAYS=20. C PRTD1=10. C MODE=3 C TPRNI=1.0 C DTE=15. C ISPLIT=45 C SEAMT=.F. C NREAD=0 C HORCON=0.05 C------------------------------------------------------------------- c print *, 'POINT 10' C ISPADV=5 SMOTH=0.10 UMOL=2.E-5 C ISWTCH=100000 IPRTD2=1 C For an island, DELH>1.0; for a seamount DELH<1.0. DELX=24.E3 DELH=0.9 RA=25.E3 H0=4500. C Flow magnitude VEL=0.2 C One can create a 'params' file in a runscript to overwrite parameters. INCLUDE 'params' C C----------------------------------------------------------------------- C Herewith is a list of parameters that are user discretionary. C----------------------------------------------------------------------- C In unst.h C IM, JM, KB = grid dimension C C In MAIN C SEAMT=.T. The ready-to-run seamount probelem will be executed. C SEAMT=.F. Code expect user prepared input for READ(40) which C will overwrite seamount problem. C MODE = 2; 2-D CALCULATION (BOTTOM STRESS CALCULATED IN ADVAVE) C 3; 3-D CALCULATION (BOTTOM STRESS CALCULATED IN PROFU,V) C 4; 3-D CALCULATION WITH T AND S HELD FIXED C DTE = EXTERNAL TIME STEP C ISPLIT = DTI/DTE C DAYS = LENGTH OF RUN IN DAYS C PRTD1 = PRINT INTERVAL IN DAYS; AFTER ISWTCH DAYS, PRINT C INTERVAL = IPRTD2 C NREAD=0, NO RESTART INPUT FILE; NREAD=1, RESTART C ISPADV = STEP INTERVAL WHEREIN EXTERNAL MODE ADVECTIVE TERMS ARE C NOT UPDATED C HORCON = CONSTANT IN SMAGORINSKY HORIZONTAL VISCOSITY C TPRNI = HORIZONTAL DIFFUSIVITY INVERSE PRANDTL NUMBER C (see S.R. ADVT) C UMOL = Background viscosity used in S.R. PROFT, PROFU and PROFV. C SMOTH = CONSTANT IN TIME SMOOTHER TO PREVENT SOLUTION SPLITTING C TBIAS, SBIAS = Value set in data statement. For 32 bit arithmetic, C setting S=35. and T=10. will reduce round-off error. C RA, HO, DELH = Island/seamount topographical parameters C RAMP = Device to slowly ramp up wind stress and baroclinic forcing. C It can be used to eliminate initial inertial oscillations. C However, RAMP is set to a constant 1.0 below and is C therefor disabled. C C In BCOND C HMAX = maximum depth C C In DEPTH C KL1, KL2 = # of logarithmic grid points in surface and bottom C respectively C C In PROFT C NBC = choice of surface boundary conditions C NTP = choice of radiative penetration parameter. C C In SLPMIN C SLMIN = see subroutine C-------------------------------------------------------------------------- DTI=DTE*FLOAT(ISPLIT) DTE2=DTE*2 DTI2=DTI*2 Clh01/14 c IEND=DAYS*24*3600/DTI+1 IEND=DAYS*24*3600/IFIX(DTI) IPRINT=PRTD1*24*3600/DTI ISWTCH=ISWTCH*24*3600/INT(DTI) c print *, 'POINT 11' write(10,'(''DAYS ='',F10.2)') DAYS WRITE(10,7030) MODE,DTE,DTI,IEND,ISPLIT,ISPADV,IPRINT,SMOTH, 1 HORCON,TPRNI 7030 FORMAT(//,' MODE =',I3, 1 ' DTE =',F7.2,' DTI =',F9.1,' IEND =',I6, 2 ' ISPLIT =',I6,' ISPADV =',I6,' IPRINT =',I6,/, 3 ' SMOTH =',F7.2,' HORCON =',F7.3,' TPRNI ='F7.3,/) WRITE(10,'('' NREAD ='',I5,//)') NREAD C ISPI=1.E0/FLOAT(ISPLIT) ISP2I=1.E0/(2.E0*FLOAT(ISPLIT)) C---------------------------------------------------------------------- C ESTABLISH PROBLEM CHARACTERISTICS C ****** ALL UNITS IN M.K.S. SYSTEM ****** C F,BLANK AND B REFERS TO FORWARD,CENTRAL AND BACKWARD TIME LEVELS. C---------------------------------------------------------------------- C C IP=IM ISK=5 JSK=5 C--------------------------------------------------------------------- jmd2=jm/2 C******************************************************************** C READ IN USER SUPPLIED GRID DATA AND INITIAL CONDITIONS C----------------------------------------------------------------------- REWIND(40) READ(40) KBB,Z,ZZ,DZ,DZZ,IMM,JMM,ALON,ALAT,DX,DY,H,FSM, 1 DUM,DVM,ART,ARU,ARV,COR C----------------------------------------------------------------------- C C In this application, the initial condition for density is C a function of z (cartesian). When this is not so, make sure C that RMEAN has been area averaged before transfer to sigma c print *, 'POINT 12' C coordinates. DO 51 J=1,JM DO 51 I=1,IM UAB(I,J)=0.0 VAB(I,J)=0.0 ELB(I,J)=0.0 ETB(I,J)=0.0 AAM2D(I,J)=500.0 51 CONTINUE C******************************************************************** C IF(SEAMT) GOTO 7005 TIME0=0. C C In this user application, the north,south and west boundaries are C closed whereas the east boundaries are open. C DO I=1,IM H(I,JM)=1. H(I,1)=1. ENDDO DO J=1,JM H(1,J)=1. H(IM,J)=H(IMM1,J) ENDDO C++++++++++++++++++++++ FSM(:,jm)=0.0 C++++++++++++++++++++++ C------------------------------------------------ C Set lateral boundary conditions for use in S.R.BCOND. As seen C in that S.R., if, for example, on the eastern boundary, RFE C =0. and UABE=0., the vertically mean velocity will be nil, but C the internal (baroclinic) velocities might not be. C------------------------------------------------ RFE=1. DO J=1,JM UABE(J)=0. ELE(J)=0. ENDDO 7005 CONTINUE C******************************************************************** do j=86,87 fsm(138,j)=0.0 enddo do i=280,282 fsm(i,238)=0.0 enddo do i=323,325 fsm(i,240)=0.0 enddo C************************ DO 35 J=2,JM DO 35 I=2,IM DUM(I,J)=FSM(I,J)*FSM(I-1,J) DVM(I,J)=FSM(I,J)*FSM(I,J-1) 35 CONTINUE C************************ for linearization do j=1,jm do i=1,im cbc(i,j)=0.0025e0 enddo enddo C************************ c print *, 'POINT 13' C C A very empirical specification of the bottom roughness C parameter follows c DO 45 J=1,JM c DO 45 I=1,IM c DT(I,J)=H(I,J) c Z0B=.01 c CBCMIN=.0025E0 c 45 CBC(I,J)=MAX(CBCMIN,.16E0/LOG((ZZ(KBM1)-Z(KB))*H(I,J) c 1 /Z0B)**2)*FSM(I,J) C c CALL PRXY(' TOPOGRAPHY 2 ',TIME, H,IP,1,JM,1,10.) c CALL PRXY(' FSM ',TIME,FSM,IP,1,JM,1,1.) c CALL PRXY(' DUM ',TIME,DUM,IP,1,JM,1,1.) c CALL PRXY(' DVM ',TIME,DVM,IP,1,JM,1,1.) c CALL PRXY(' CBC ',TIME,CBC,IP,1,JM,1,0.) C c DO 47 J=1,JM c DO 47 I=1,IM c 47 TPS(I,J)=0.5E0/SQRT(1.E0/DX(I,J)**2+1.E0/DY(I,J)**2) c 1 /SQRT(GRAV*H(I,J))*FSM(I,J) c CALL PRXY(' EXT. CFL TIME STEP ',TIME,TPS,IP,ISK,JM,JSK,1.E0) c CALL PRXY(' COR ',TIME,COR,IP,ISK,JM,1,0.E0) DO 62 I=1,IM DO 62 J=1,JM UA(I,J)=UAB(I,J) VA(I,J)=VAB(I,J) EL(I,J)=ELB(I,J) ET(I,J)=ETB(I,J) ETF(I,J)=ET(I,J) D(I,J)=H(I,J)+EL(I,J) DT(I,J)=H(I,J)+ET(I,J) DRX2D(I,J)=0. DRY2D(I,J)=0. 62 CONTINUE C C RAMP=1. C C++++++++++++++LHC DO 81 J=1,JM DO 81 I=1,IM D(I,J)=H(I,J)+EL(I,J) 81 DT(I,J)=H(I,J)+ET(I,J) c print *, 'POINT 14' C TIME=TIME0 c CALL PRXY(' FREE SURFACE ',TIME, ELB,IP,ISK,JM,JSK,1.E-3) c CALL PRXY(' UAB ',TIME, UAB,IP,ISK,JM,JSK,0.) c CALL PRXY(' VAB ',TIME, VAB,IP,ISK,JM,JSK,0.) PERIOD=(2.E0*PI)/COR(IM/2,JM/2)/86400. C*********************************************************************** C * C BEGIN NUMERICAL INTEGRATION * C * C*********************************************************************** C C++++++++++++++LHC c print *, 'POINT 15' C++++++++++++++LHC C*********************************************************************** NO_TIDE=1 IF(NO_TIDE .NE. 0) CALL TIDE_SETUP C*********************************************************************** DO 9000 IINT=1,IEND C IF(TIME.GT. 5.) MODE=3 C time=float(iint)*dte*float(isplit)/86400. period=10.0 RAMP=TIME/PERIOD IF (RAMP.GT.1.0) RAMP=1.0 C++++++++++++++LHC write(*,*) time C----------------------------------------------------------------------- C SET TIME DEPENDENT, SURFACE AND LATERAL BOUNDARY CONDITIONS. C THE LATTER WILL BE USED IN SUBROUTINE BCOND. Users may C wish to create a subroutine to supply WUSURF, WVSURF, WTSURF, C WSSURF and SWRAD. C----------------------------------------------------------------------- C++++++++++++++LHC c write(*,*) 'POINT 1' C----------------------------------------------------------------------- c write(*,*) 'POINT 2' C----------------------------------------------------------------------- C WIND FORCING C----------------------------------------------------------------------- C INTRODUCE SIMPLE WIND STRESS, VALUE IS NEGATIVE FOR WESTERLY WIND C The following wind stress has been tapered along the boundary to suppress C numerically induced oscilations near the boundary (Jamart&Ozer, C J.G.R.,91, 10621-10631) cd cd---------------------------------------------------------------------- WTSURF=0.0 WUSURF=0.0 WVSURF=0.0 SWRAD =0.0 EAT =0.0 cd---------------------------------------------------------------------- C SURFACE ELIVATION ALONG THE EASTERN BOUNDARY C----------------------------------------------------------------------- ELE(1)=0.0 DO J=2,JMM1 ELE(J)=ELE(J-1)-COR(IMM1,J)*UABE(J)/GRAV*DY(IMM1,J) ENDDO DO J=2,JMM1 c ELE(J)=(ELE(J)-ELE(JMM1/2))*FSM(IM,J)*RAMP ELE(J)=(ELE(J)-ELE(JMM1/2))*FSM(IM,J) ENDDO c write(*,*) 'POINT 6' C----------------------------------------------------------------------- C TIDAL FORCING C----------------------------------------------------------------------- C c write(*,*) 'POINT 7' IF(MODE.EQ.2) GOTO 8001 8001 CONTINUE C------------------------------------------------------------------------ DO 399 J=1,JM DO 399 I=1,IM 399 EGF(I,J)=EL(I,J)*ISPI DO 400 J=2,JM DO 400 I=2,IM UTF(I,J)=UA(I,J)*(D(I,J)+D(I-1,J))*ISP2I 400 VTF(I,J)=VA(I,J)*(D(I,J)+D(I,J-1))*ISP2I write(*,*) 'IINT=', iint C********** BEGIN EXTERNAL MODE *************************************** DO 8000 IEXT=1,ISPLIT c write(*,*) 'POINT 100' c write(10,'('' IEXT,TIME ='',I5,F9.2)') IEXT,TIME C********************************************************************* hh=0.0 dd=1.0 mm=1.0 yy=95.0 IF(NO_TIDE .NE. 0) > CALL TIDE_BODY(hh,dd,mm,yy,DTE,IEXT,EEQ,ELF,ALAT,ALON) C********************************************************************* cd------------------- for tidal response c EEQ=0.0 c********************************************************************* DO 405 J=2,JM DO 405 I=2,IM FLUXUA(I,J)=.25E0*(D(I,J)+D(I-1,J))*(DY(I,J)+DY(I-1,J))*UA(I,J) 405 FLUXVA(I,J)=.25E0*(D(I,J)+D(I,J-1))*(DX(I,J)+DX(I,J-1))*VA(I,J) C c write(*,*) 'POINT 101' DO 410 J=2,JMM1 DO 410 I=2,IMM1 410 ELF(I,J)=ELB(I,J) 1 -DTE2*(FLUXUA(I+1,J)-FLUXUA(I,J)+FLUXVA(I,J+1)-FLUXVA(I,J)) 2 /ART(I,J) C CALL BCOND(1,IEXT) C cd IF(MOD(IEXT,ISPADV).EQ.0) CALL ADVAVE(ADVUA,ADVVA,MODE) cd C*********************** Linearization ADVUA=0.0 ADVVA=0.0 C*********************** C Note that ALPHA = 0. is perfectly acceptable. The value, ALPHA = .225 C permits a longer time step. ALPHA=0.225 C*********************** if (mode.eq.2) then ADX2D=0.0 ADY2D=0.0 endif DO 420 J=2,JMM1 DO 420 I=2,IM cd c advua(i,j)=0.0 c wusurf(i,j)=0.0 c wubot(i,j)=0.002*uab(i,j) cd 420 UAF(I,J)=ADX2D(I,J)+ADVUA(I,J) 1 -ARU(I,J)*.25*( COR(I,J)*D(I,J)*(VA(I,J+1)+VA(I,J)) 2 +COR(I-1,J)*D(I-1,J)*(VA(I-1,J+1)+VA(I-1,J)) ) 3 +.25E0*GRAV*(DY(I,J)+DY(I-1,J))*(D(I,J)+D(I-1,J)) 4 *( (1.E0-2.E0*ALPHA)*(EL(I,J)-EL(I-1,J)) c11/97before 4 +ALPHA*(ELB(I,J)-ELB(I-1,J)+ELF(I,J)-ELF(I-1,J)) ) 4 +ALPHA*(ELB(I,J)-ELB(I-1,J)+ELF(I,J)-ELF(I-1,J)) & -(EAT(I,J)-EAT(I-1,J))*dum(i,j) C equilibrium tides 4 -(EEQ(I,J)-EEQ(I-1,J))*RAMP ) ct & ) 5 +DRX2D(I,J) c11/97before 6 +ARU(I,J)*( WUSURF(I,J)-WUBOT(I,J) ) 6 +ARU(I,J)*( 0.5*(WUSURF(I,J)+WUSURF(I-1,J))-WUBOT(I,J) ) c write(*,*) 'POINT 105' DO 425 J=2,JMM1 DO 425 I=2,IM cd c advva(i,j)=0.0 c wvsurf(i,j)=0.0 c wvbot(i,j)=0.002*vab(i,j) cd 425 UAF(I,J)= 1 ((H(I,J)+ELB(I,J)+H(I-1,J)+ELB(I-1,J))*ARU(I,J)*UAB(I,J) 2 -4.E0*DTE*UAF(I,J)) 3 /((H(I,J)+ELF(I,J)+H(I-1,J)+ELF(I-1,J))*ARU(I,J)) c write(*,*) 'POINT 106' DO 430 J=2,JM DO 430 I=2,IMM1 VAF(I,J)=ADY2D(I,J)+ADVVA(I,J) 1 +ARV(I,J)*.25*( COR(I,J)*D(I,J)*(UA(I+1,J)+UA(I,J)) 2 +COR(I,J-1)*D(I,J-1)*(UA(I+1,J-1)+UA(I,J-1)) ) 3 +.25E0*GRAV*(DX(I,J)+DX(I,J-1))*(D(I,J)+D(I,J-1)) 4 *( (1.E0-2.E0*ALPHA)*(EL(I,J)-EL(I,J-1)) c11/97before 4 +ALPHA*(ELB(I,J)-ELB(I,J-1)+ELF(I,J)-ELF(I,J-1)) ) 4 +ALPHA*(ELB(I,J)-ELB(I,J-1)+ELF(I,J)-ELF(I,J-1)) & -(EAT(I,J)-EAT(I,J-1))*dvm(i,j) C equilibrium tides 4 -(EEQ(I,J)-EEQ(I,J-1))*RAMP ) ct & ) 5 +DRY2D(I,J) c11/97before 6 + ARV(I,J)*( WVSURF(I,J)-WVBOT(I,J) ) 6 + ARV(I,J)*( 0.5*(WVSURF(I,J)+WVSURF(I,J-1))-WVBOT(I,J) ) 430 CONTINUE DO 435 J=2,JM DO 435 I=2,IMM1 VAF(I,J)= 1 ((H(I,J)+ELB(I,J)+H(I,J-1)+ELB(I,J-1))*VAB(I,J)*ARV(I,J) 2 -4.E0*DTE*VAF(I,J)) 3 /((H(I,J)+ELF(I,J)+H(I,J-1)+ELF(I,J-1))*ARV(I,J)) 435 CONTINUE CALL BCOND(2,IEXT) C IF(IEXT.LT.(ISPLIT-2)) GO TO 440 IF(IEXT.EQ.(ISPLIT-2))THEN DO 431 J=1,JM DO 431 I=1,IM 431 ETF(I,J)=.25*SMOTH*ELF(I,J) ENDIF IF(IEXT.EQ.(ISPLIT-1)) THEN DO 432 J=1,JM DO 432 I=1,IM 432 ETF(I,J)=ETF(I,J)+.5*(1.-.5*SMOTH)*ELF(I,J) ENDIF IF(IEXT.EQ.(ISPLIT-0)) THEN DO 433 J=1,JM DO 433 I=1,IM 433 ETF(I,J)=(ETF(I,J)+.5*ELF(I,J))*FSM(I,J) ENDIF 440 CONTINUE C C TEST FOR CFL VIOLATION. IF SO, PRINT AND STOP C VMAXL=100. VAMAX=0. DO 442 J=1,JM DO 442 I=1,IM IF(ABS(VAF(I,J)).GE.VAMAX) THEN VAMAX=ABS(VAF(I,J)) IMAX=I JMAX=J ENDIF 442 CONTINUE IF(VAMAX.GT.VMAXL) GO TO 9001 C C APPLY FILTER TO REMOVE TIME SPLIT. RESET TIME SEQUENCE. DO 445 J=1,JM DO 445 I=1,IM UA(I,J)=UA(I,J)+.5E0*SMOTH*(UAB(I,J)-2.E0*UA(I,J)+UAF(I,J)) VA(I,J)=VA(I,J)+.5E0*SMOTH*(VAB(I,J)-2.E0*VA(I,J)+VAF(I,J)) EL(I,J)=EL(I,J)+.5E0*SMOTH*(ELB(I,J)-2.E0*EL(I,J)+ELF(I,J)) ELB(I,J)=EL(I,J) EL(I,J)=ELF(I,J) D(I,J)=H(I,J)+EL(I,J) UAB(I,J)=UA(I,J) UA(I,J)=UAF(I,J) VAB(I,J)=VA(I,J) VA(I,J)=VAF(I,J) 445 CONTINUE C IF(IEXT.EQ.ISPLIT) GO TO 8000 DO 450 J=2,JM DO 450 I=2,IM EGF(I,J)=EGF(I,J)+EL(I,J)*ISPI UTF(I,J)=UTF(I,J)+UA(I,J)*(D(I,J)+D(I-1,J))*ISP2I 450 VTF(I,J)=VTF(I,J)+VA(I,J)*(D(I,J)+D(I,J-1))*ISP2I C C 8000 CONTINUE IF(VAMAX.GT.VMAXL) GO TO 9001 C--------------------------------------------------------------------- C END EXTERNAL (2-D) MODE CALCULATION C AND CONTINUE WITH INTERNAL (3-D) MODE CALCULATION C--------------------------------------------------------------------- IF(IINT.EQ.1.AND.TIME0.EQ.0.) GO TO 8200 IF(MODE.EQ.2) GO TO 8200 C--------------------------------------------------------------------- C C 8200 CONTINUE C DO 8210 J=1,JM DO 8210 I=1,IM EGB(I,J)=EGF(I,J) ETB(I,J)=ET(I,J) ET(I,J)=ETF(I,J) DT(I,J)=H(I,J)+ET(I,J) UTB(I,J)=UTF(I,J) 8210 VTB(I,J)=VTF(I,J) C C C-------------------------------------------------------------------- C BEGIN PRINT SECTION C-------------------------------------------------------------------- C----------------- SEA LEVEL OUTPUT if(mod(iint,5).eq.0) then nunit=47 CALL SAVEL(NUNIT,hh,dd,mm,yy,EL,IM,JM) endif C----------------- SEA LEVEL OUTPUT for check ETIDE sem(:)=el(im-1,:) write(48) hh,dd,mm,yy,sem C----------------- SEA LEVEL OUTPUT for UTIDE sem(:)=UAB(im-1,:) write(49) hh,dd,mm,yy,sem C-------------------------------------------------------------------- cd if(mod(iint,60).eq.0) then c CALL FINDPSI c dpsi=(psi(141,92)-psi(152,88))*1.e-6 c write(46,'(f7.4,3f5.0,1x,f20.10)') hh,dd,mm,yy,dpsi cd endif Chc------------------------------------------------------------------ ct IF(MOD(IINT,IPRINT).NE.0) GO TO 7000 IMPR=1 IF(IMPR .EQ. 0) GOTO 7000 Chc------------------------------------------------------------------ 9001 CONTINUE IF(IINT.GE.ISWTCH) IPRINT=IPRTD2*24*3600/INT(DTI) WRITE(10,'(/,'' TIME ='',F10.2,'' IINT ='',I8, 1 '' IEXT ='',I8,'' IPRINT ='',I5,//)') TIME,IINT,IEXT,IPRINT C IK=304 JK=200 IF(VAMAX.GT.VMAXL) THEN CALL PRXYM(' EL ',TIME, EL ,IM,IK,JM,JK,0.E0) CALL PRXYM(' ELF',TIME, ELF,IM,IK,JM,JK,0.E0) CALL PRXYM(' UA ',TIME,UA ,IM,IK,JM,JK,0.E0) CALL PRXYM(' UAF ',TIME,UAF,IM,IK,JM,JK,0.E0) CALL PRXYM(' VA ',TIME,VA ,IM,IK,JM,JK,0.E0) CALL PRXYM(' VAF ',TIME,VAF,IM,IK,JM,JK,0.E0) c CALL PRXY(' EL ',TIME, EL ,IP,ISK,JM,1,0.E0) c CALL PRXY(' ELF',TIME, ELF,IP,ISK,JM,1,0.E0) c CALL PRXY(' UA ',TIME,UA ,IP,ISK,JM,JSK,0.E0) c CALL PRXY(' UAF ',TIME,UAF,IP,ISK,JM,JSK,0.E0) c CALL PRXY(' VA ',TIME,VA ,IP,ISK,JM,JSK,0.E0) c CALL PRXY(' VAF ',TIME,VAF,IP,ISK,JM,JSK,0.E0) WRITE(10,'(//////////////////)') WRITE(10,'(''************************************************'')') WRITE(10,'(''************ ABNORMAL JOB END ******************'')') WRITE(10,'(''************* USER TERMINATED ******************'')') WRITE(10,'(''************************************************'')') WRITE(10,'('' VAMAX ='',E12.3,'' IMAX,JMAX ='',2I5)') 1 VAMAX,IMAX,JMAX WRITE(10,'('' IINT,IEXT ='',2I6)') IINT,IEXT WRITE(41) TIME,UA,VA,EL STOP ENDIF 7000 CONTINUE c IF(MOD(IINT,IPRINT/2).EQ.0) c IF(MOD(IINT,IPRINT).EQ.0) Clhc 1/20; in every 15th day save output c if(HH.eq.0.0 .and. DD.eq.15.0) C-------------- OUTPUT in the last day C*********************** c pday=days-1.0 c prt=(pday*86400./(dte*float(isplit))) c if(float(iint).ge.prt) then c IF(MOD(IINT,5).EQ.0) c 1 WRITE(41) TIME,UA,VA,EL c endif C*********************** pindex=0.0 ptime=time+1.0 if (pindx.eq.1.0 .or. mod(int(ptime),30).eq. 0) then if (mod(int(time),30) .eq. 0) then pindex=0.0 IF(MOD(IINT,5).EQ.0) WRITE(41) TIME,UA,VA,EL else pindex=1.0 IF(MOD(IINT,5).EQ.0) WRITE(41) TIME,UA,VA,EL endif endif C*********************** IF(MOD(IINT,5).EQ.0) then dumke=0.0 dumpe=0.0 dumv=0.0 duma=0.0 do j=1,jm do i=1,im dumke=dumke+0.5*((0.5*(UA(i,j)+UA(i+1,j)))**2 & +(0.5*(VA(i,j)+VA(i,j+1)))**2) & *dx(i,j)*dy(i,j)*(h(i,j)+el(i,j))*fsm(i,j) dumpe=dumpe+0.5*9.8*el(i,j)*el(i,j)*dx(i,j)*dy(i,j)*fsm(i,j) dumv=dumv+dx(i,j)*dy(i,j)*(h(i,j)+el(i,j))*fsm(i,j) duma=duma+dx(i,j)*dy(i,j)*fsm(i,j) enddo enddo tke=dumke/duma tpe=dumpe/duma write(42,'(f10.3,2e20.7)') time, tke,tpe ENDIF C-------------------------------------------------------------------- C END PRINT SECTION C-------------------------------------------------------------------- 9000 CONTINUE C*********************************************************************** C * C END NUMERICAL INTEGRATION * C * C*********************************************************************** C STOP END C SUBROUTINE ADVAVE(ADVUA,ADVVA,MODE) INCLUDE 'cof_a.h' DIMENSION ADVUA(IM,JM),ADVVA(IM,JM),CURV2D(IM,JM) EQUIVALENCE (TPS,CURV2D) C--------------------------------------------------------------------- C CALCULATE U-ADVECTION & DIFFUSION C--------------------------------------------------------------------- C C-------- ADVECTIVE FLUXES ------------------------------------------- C DO 200 J=1,JM DO 200 I=1,IM 200 ADVUA(I,J)=0. DO 300 J=2,JM DO 300 I=2,IMM1 300 FLUXUA(I,J)=.125E0*((D(I+1,J)+D(I,J))*UA(I+1,J) 1 +(D(I,J)+D(I-1,J))*UA(I,J)) 2 *(UA(I+1,J)+UA(I,J)) DO 400 J=2,JM DO 400 I=2,IM 400 FLUXVA(I,J)=.125E0*((D(I,J)+D(I,J-1))*VA(I,J) 1 +(D(I-1,J)+D(I-1,J-1))*VA(I-1,J)) 2 *(UA(I,J)+UA(I,J-1)) C----------- ADD VISCOUS FLUXES --------------------------------- DO 460 J=2,JM DO 460 I=2,IMM1 460 FLUXUA(I,J)=FLUXUA(I,J) 1 -D(I,J)*2.E0*AAM2D(I,J)*(UAB(I+1,J)-UAB(I,J))/DX(I,J) DO 470 J=2,JM DO 470 I=2,IM TPS(I,J)=.25E0*(D(I,J)+D(I-1,J)+D(I,J-1)+D(I-1,J-1)) 1 *(AAM2D(I,J)+AAM2D(I,J-1)+AAM2D(I-1,J)+AAM2D(I-1,J-1)) 2 *((UAB(I,J)-UAB(I,J-1)) 3 /(DY(I,J)+DY(I-1,J)+DY(I,J-1)+DY(I-1,J-1)) 4 +(VAB(I,J)-VAB(I-1,J)) 5 /(DX(I,J)+DX(I-1,J)+DX(I,J-1)+DX(I-1,J-1)) ) FLUXUA(I,J)=FLUXUA(I,J)*DY(I,J) FLUXVA(I,J)=(FLUXVA(I,J)-TPS(I,J)) 1 *.25E0*(DX(I,J)+DX(I-1,J)+DX(I,J-1)+DX(I-1,J-1)) 470 CONTINUE C---------------------------------------------------------------- DO 480 J=2,JMM1 DO 480 I=2,IMM1 480 ADVUA(I,J)=FLUXUA(I,J)-FLUXUA(I-1,J) 1 +FLUXVA(I,J+1)-FLUXVA(I,J) C---------------------------------------------------------------- C CALCULATE V-ADVECTION & DIFFUSION C---------------------------------------------------------------- C DO 481 J=1,JM DO 481 I=1,IM 481 ADVVA(I,J)=0. C C---------ADVECTIVE FLUXES ---------------------------- DO 700 J=2,JM DO 700 I=2,IM 700 FLUXUA(I,J)=.125E0*((D(I,J)+D(I-1,J))*UA(I,J) 1 +(D(I,J-1)+D(I-1,J-1))*UA(I,J-1))* 2 (VA(I-1,J)+VA(I,J)) DO 800 J=2,JMM1 DO 800 I=2,IM 800 FLUXVA(I,J)=.125E0*((D(I,J+1)+D(I,J)) 1 *VA(I,J+1)+(D(I,J)+D(I,J-1))*VA(I,J)) 2 *(VA(I,J+1)+VA(I,J)) C------- ADD VISCOUS FLUXES ----------------------------------- DO 860 J=2,JMM1 DO 860 I=2,IM 860 FLUXVA(I,J)=FLUXVA(I,J) 1 -D(I,J)*2.E0*AAM2D(I,J)*(VAB(I,J+1)-VAB(I,J))/DY(I,J) DO 870 J=2,JM DO 870 I=2,IM FLUXVA(I,J)=FLUXVA(I,J)*DX(I,J) 870 FLUXUA(I,J)=(FLUXUA(I,J)-TPS(I,J)) 3 *.25E0*(DY(I,J)+DY(I-1,J)+DY(I,J-1)+DY(I-1,J-1)) C--------------------------------------------------------------- DO 880 J=2,JMM1 DO 880 I=2,IMM1 880 ADVVA(I,J)=FLUXUA(I+1,J)-FLUXUA(I,J) 1 +FLUXVA(I,J)-FLUXVA(I,J-1) C C--------------------------------------------------------------- IF(MODE.NE.2) GO TO 5000 DO 100 J=2,JMM1 DO 100 I=2,IMM1 WUBOT(I,J)=-0.5E0*(CBC(I,J)+CBC(I-1,J)) 1 *SQRT(UAB(I,J)**2+(.25E0*(VAB(I,J) 2 +VAB(I,J+1)+VAB(I-1,J)+VAB(I-1,J+1)))**2)*UAB(I,J) 100 CONTINUE DO 102 J=2,JMM1 DO 102 I=2,IMM1 WVBOT(I,J)=-0.5E0*(CBC(I,J)+CBC(I,J-1)) 1 *SQRT((.25E0*(UAB(I,J)+UAB(I+1,J) 2 +UAB(I,J-1)+UAB(I+1,J-1)))**2+VAB(I,J)**2)*VAB(I,J) 102 CONTINUE DO 120 J=2,JMM1 DO 120 I=2,IMM1 CURV2D(I,J)=.25*((VA(I,J+1)+VA(I,J))*(DY(I+1,J)-DY(I-1,J)) 1 -(UA(I+1,J)+UA(I,J))*(DX(I,J+1)-DX(I,J-1)) ) 2 /(DX(I,J)*DY(I,J)) 120 CONTINUE DO 130 J=3,JMM2 DO 130 I=3,IMM2 ADVUA(I,J)=ADVUA(I,J) 1 -ARU(I,J)*.25*( CURV2D(I,J)*D(I,J)*(VA(I,J+1)+VA(I,J)) 2 +CURV2D(I-1,J)*D(I-1,J)*(VA(I-1,J+1)+VA(I-1,J)) ) ADVVA(I,J)=ADVVA(I,J) 1 +ARV(I,J)*.25*( CURV2D(I,J)*D(I,J)*(UA(I+1,J)+UA(I,J)) 2 +CURV2D(I,J-1)*D(I,J-1)*(UA(I+1,J-1)+UA(I,J-1)) ) 130 CONTINUE 5000 CONTINUE C RETURN END C C SUBROUTINE BCOND(IDX,IEXT) INCLUDE 'cof_a.h' REAL hh,dd,mm,yy PARAMETER (NC=6) REAL MASK_T COMMON/TIDE1/AMP(NC),OME(NC),CHI(NC),MASK_T(NC) COMMON/TIDE2/EAMPE(JM,NC),EPHAE(JM,NC), 1 UAMPE(JM,NC),UPHAE(JM,NC), 2 VAMPE(JM,NC),VPHAE(JM,NC) C C Closed boundary conditions are automatically enabled through C specification of the masks, DUM, DVM and FSM in which case open C boundary condition, included below, will be overwritten. C C C SPECIFICATION OF OPEN BOUNDARY CONDITIONS C BRACKETS OF * INDICATE INTERIOR (NON-BOUNDARY) GRID POINTS. C HORIZONTAL LOCATIONS OF EL, T AND S ARE COINCIDENT. C C I-1 I I+1 C C C U(JM) EL(JM) =. U(JM) BC C . C V(JM) . BC C . C *U(JMM1)* *EL(JMM1)* *U(JMM1)* C C *V(JMM1)* C . *V(IMM1)* V(IM) C . C . *U(IMM1)* *EL(IMM1)* U(IM) EL(IM) C . ". . . . . . . ." C . C . *V(IMM1)* V(IM) C . C . BC BC C C *V(3)* C C *U(2)* *EL(2)* *U(2)* C . C V(2) . BC C . C U(1) EL(1)=. U(1) BC C -------------------------------------------------------------------- C V(1) IS NOT USED C C DATA PI/3.14167/,GEE/9.807/,RAD/0.017453292/ C C DXX=(XXTIME(hh,dd,mm,yy)-XXTIME(0.,1.,1.,yy))/24.0 dxx=float(iint*isplit)*dte C TIS=DXX+IEXT*DTE/86400. c TIS=(DXX+IEXT*DTE)/86400. c TIB=(DXX+IEXT*DTE-DTE)/86400. TIS=(float(IINT*isplit+IEXT)*dte)/86400.0 TIB=(float(IINT*isplit+IEXT-1)*dte)/86400.0 if(TIB .le.0.0) TIB=0.0 C period=10.0 RAMP=TIS/PERIOD IF (RAMP.GT.1.0) RAMP=1.0 C GO TO (10,20,30,40,50,60), IDX C C----------------------------------------------------------------------- C EXTERNAL B.C.'S C----------------------------------------------------------------------- C C 10 CONTINUE C---------ELEVATION -------------------- c write(*,*) TIS C----------------------------- C In this application, elevation is not a primary B.C. DO 120 J=1,JM ELF(1,J)=ELF(2,J) c ELF(IM,J)=ELF(IMM1,J) ETIDE= EAMPE(J,1)*COS(OME(1)*TIS+CHI(1)-RAD*EPHAE(J,1))*MASK_T(1) 1 +EAMPE(J,2)*COS(OME(2)*TIS+CHI(2)-RAD*EPHAE(J,2))*MASK_T(2) 1 +EAMPE(J,3)*COS(OME(3)*TIS+CHI(3)-RAD*EPHAE(J,3))*MASK_T(3) 1 +EAMPE(J,4)*COS(OME(4)*TIS+CHI(4)-RAD*EPHAE(J,4))*MASK_T(4) 1 +EAMPE(J,5)*COS(OME(5)*TIS+CHI(5)-RAD*EPHAE(J,5))*MASK_T(5) 1 +EAMPE(J,6)*COS(OME(6)*TIS+CHI(6)-RAD*EPHAE(J,6))*MASK_T(6) ELF(IM,J)=ETIDE*RAMP ELF(IM-1,J)=ELF(IM,J) 120 CONTINUE DO 130 I=1,IM ELF(I,1)=ELF(I,2) ELF(I,JM)=ELF(I,JMM1) 130 CONTINUE C DO 140 J=1,JM DO 140 I=1,IM 140 ELF(I,J)=ELF(I,J)*FSM(I,J) RETURN C C 20 CONTINUE C---------- VELOCITY -------------- C--- Governing external B.C.s ------------- ramp=1.0 C----------------------------------------------------------------------- C ***** EAST *********** DO 220 J=1,JM c UTIDE= UAMPE(J,1)*COS(OME(1)*TIS+CHI(1)-RAD*UPHAE(J,1))*MASK_T(1) c 1 +UAMPE(J,2)*COS(OME(2)*TIS+CHI(2)-RAD*UPHAE(J,2))*MASK_T(2) c 1 +UAMPE(J,3)*COS(OME(3)*TIS+CHI(3)-RAD*UPHAE(J,3))*MASK_T(3) c 1 +UAMPE(J,4)*COS(OME(4)*TIS+CHI(4)-RAD*UPHAE(J,4))*MASK_T(4) c 1 +UAMPE(J,5)*COS(OME(5)*TIS+CHI(5)-RAD*UPHAE(J,5))*MASK_T(5) c 1 +UAMPE(J,6)*COS(OME(6)*TIS+CHI(6)-RAD*UPHAE(J,6))*MASK_T(6) ETIDE= EAMPE(J,1)*COS(OME(1)*TIS+CHI(1)-RAD*EPHAE(J,1))*MASK_T(1) 1 +EAMPE(J,2)*COS(OME(2)*TIS+CHI(2)-RAD*EPHAE(J,2))*MASK_T(2) 1 +EAMPE(J,3)*COS(OME(3)*TIS+CHI(3)-RAD*EPHAE(J,3))*MASK_T(3) 1 +EAMPE(J,4)*COS(OME(4)*TIS+CHI(4)-RAD*EPHAE(J,4))*MASK_T(4) 1 +EAMPE(J,5)*COS(OME(5)*TIS+CHI(5)-RAD*EPHAE(J,5))*MASK_T(5) 1 +EAMPE(J,6)*COS(OME(6)*TIS+CHI(6)-RAD*EPHAE(J,6))*MASK_T(6) ETIDB= EAMPE(J,1)*COS(OME(1)*TIb+CHI(1)-RAD*EPHAE(J,1))*MASK_T(1) 1 +EAMPE(J,2)*COS(OME(2)*TIb+CHI(2)-RAD*EPHAE(J,2))*MASK_T(2) 1 +EAMPE(J,3)*COS(OME(3)*TIb+CHI(3)-RAD*EPHAE(J,3))*MASK_T(3) 1 +EAMPE(J,4)*COS(OME(4)*TIb+CHI(4)-RAD*EPHAE(J,4))*MASK_T(4) 1 +EAMPE(J,5)*COS(OME(5)*TIb+CHI(5)-RAD*EPHAE(J,5))*MASK_T(5) 1 +EAMPE(J,6)*COS(OME(6)*TIb+CHI(6)-RAD*EPHAE(J,6))*MASK_T(6) cd change the sign 6/23 c utide=ua(im-1,j)-(dx(im,j)/dte)*(etide-etidb) c & /(el(imm1,j)+h(imm1,j)) utide=0.0 cd UAF(IM,J)=RAMP*(UABE(J)+UTIDE) 1 +SQRT(GRAV/H(IMM1,J))*(EL(IMM1,J)-ELE(J)-ETIDE) UMID=.50*(UA(IM,J)+UA(IM,J-1)) VAF(IM,J)=VA(IM,J)-DTE/(DX(IM,J)+DX(IMM1,J)) 1 *( (UMID+ABS(UMID))*(VA(IM,J)-VA(IMM1,J)) 2 +(UMID-ABS(UMID))*(0.0 -VA(IM,J)) ) VAF(IM,J)=0.0 220 CONTINUE C----------------------------------------------------------------------- DO J=2,JM UAF(IM,J)=0.0 VAF(IM,J)=0.0 ENDDO C--------------- BEFORE 7/14/97 c UAF(IM,JM)=0.0 c VAF(IM,JM)=0.0 c UAF(IM,1)=0.0 c VAF(IM,1)=0.0 C--------------------------- C DO 240 J=1,JM DO 240 I=1,IM UAF(I,J)=UAF(I,J)*DUM(I,J) 240 VAF(I,J)=VAF(I,J)*DVM(I,J) RETURN C C 30 CONTINUE C----------------------------------------------------------------------- RETURN C C 40 CONTINUE C----------------------------------------------------------------------- RETURN C C 50 CONTINUE C---------------VERTICAL VEL. B. C.'S -------------------------------- RETURN C C 60 CONTINUE C---------------- Q2 AND Q2L B.C.'S ----------------------------------- RETURN END C C C C SUBROUTINE DEPTH(Z,ZZ,DZ,DZZ,KB,KBM1) C >>> DIMENSION Z(KB),ZZ(KB),DZ(KB),DZZ(KB) KL1=6 KL2=KB-2 C *** C THIS SUBROUTINE ESTABLISHES THE VERTICAL RESOLUTION WITH LOG C DISTRIBUTIONS AT THE TOP AND BOTTOM AND A LINEAR DISTRIBUTION C BETWEEN. CHOOSE KL1 AND KL2. THE DEFAULT KL1 = .3*KB AND KL2 = KB-2 C YIELDS A LOG DISTRIBUTION AT THE TOP AND NONE AT THE BOTTOM. C *** BB=FLOAT(KL2-KL1)+4. CC=FLOAT(KL1)-2. DEL1=2./BB/EXP(.693147*FLOAT(KL1-2)) DEL2=2./BB/EXP(.693147*FLOAT(KB-KL2-1)) Z(1)=0. ZZ(1)=-DEL1/2. DO 30 K=2,KL1-2 Z(K)=-DEL1*EXP(.693147*FLOAT(K-2)) 30 ZZ(K)=-DEL1*EXP(.693147*(FLOAT(K)-1.5)) DO 40 K=KL1-1,KL2+1 Z(K)=-(FLOAT(K)-CC)/BB 40 ZZ(K)=-(FLOAT(K)-CC+0.5)/BB DO 50 K=KL2+1,KBM1 Z(K)=(1.0-DEL2*EXP(.693147*FLOAT(KB-K-1)))*(-1.) 50 ZZ(K)=(1.0-DEL2*EXP(.693147*(FLOAT(KB-K)-1.5)))*(-1.) Z(KB)=-1.0 ZZ(KBM1)=-1.*(1.-DEL2/2.) ZZ(KB)=-1.*(1.+DEL2/2.) C *** C DO 55 K=1,KB C DENOM=FLOAT(KB-1) C Z(K)=-FLOAT(K-1)/DENOM C ZZ(K)=Z(K)-0.5/DENOM C 55 CONTINUE C *** DO 60 K=1,KBM1 DZ(K)=Z(K)-Z(K+1) 60 DZZ(K)=ZZ(K)-ZZ(K+1) DZ(KB)=0. DZZ(KB)=0. WRITE(10,70) 70 FORMAT(/2X,'K',7X,'Z',9X,'ZZ',9X,'DZ',9X,'DZZ',/) DO 90 K=1,KB WRITE(10,80) K,Z(K),ZZ(K),DZ(K),DZZ(K) 80 FORMAT(' ',I5,4F10.3) 90 CONTINUE RETURN END C C SUBROUTINE FINDPSI C This subroutine calculates the stream function first assuming C zero on the southern boundary and then, using the values on the C the western boundary, the stream function is calculated again. C If the elevation field is near steady state, the two calculations C should agree; otherwise not. INCLUDE 'cof_a.h' DO 9004 J=1,JM DO 9004 I=1,IM 9004 PSI(I,J)=0.E0 C Sweep northward. DO 9006 J=2,JMM1 DO 9006 I=2,IM 9006 PSI(I,J+1)=PSI(I,J)+.25E0*UAB(I,J)*(D(I,J)+D(I-1,J)) 1 *(DY(I,J)+DY(I-1,J)) c CALL PRXY(' PSI FROM U ',TIME,PSI,IM,2,JM,2,1.E6) C Sweep eastward. DO 9005 J=2,JM DO 9005 I=2,IMM1 9005 PSI(I+1,J)=PSI(I,J)-VAB(I,J)*.25E0*(D(I,J)+D(I,J-1)) 1 *(DX(I,J)+DX(I,J-1)) c CALL PRXY(' PSI FROM V ',TIME,PSI,IM,2,JM,2,1.E6) RETURN END C C C SUBROUTINE SLPMIN(H,IM,JM,FSM,SL) DIMENSION H(IM,JM),FSM(IM,JM) DIMENSION SL(IM,JM) C This subroutine limits the maximum difference in H divided C by twice the average of H of adjacent cells. C The maximum possible value is unity. C SLMIN=0.2 C DO 100 LOOP=1,10 C sweep right DO 3 J=2,JM-1 DO 1 I=2,IM-1 IF(FSM(I,J).EQ.0..OR.FSM(I+1,J).EQ.0.) GOTO 1 SL(I,J)=ABS(H(I+1,J)-H(I,J))/(H(I,J)+H(I+1,J)) IF(SL(I,J).LT.SLMIN) GOTO 1 DH=0.5*(SL(I,J)-SLMIN)*(H(I,J)+H(I+1,J)) SN=-1. IF(H(I+1,J).GT.H(I,J)) SN=1. H(I+1,J)=H(I+1,J)-SN*DH H(I,J)=H(I,J)+SN*DH 1 CONTINUE C sweep left DO 2 I=IM-1,2,-1 IF(FSM(I,J).EQ.0..OR.FSM(I+1,J).EQ.0.) GOTO 2 SL(I,J)=ABS(H(I+1,J)-H(I,J))/(H(I,J)+H(I+1,J)) IF(SL(I,J).LT.SLMIN) GOTO 2 DH=0.5*(SL(I,J)-SLMIN)*(H(I,J)+H(I+1,J)) SN=-1. IF(H(I+1,J).GT.H(I,J)) SN=1. H(I+1,J)=H(I+1,J)-SN*DH H(I,J)=H(I,J)+SN*DH 2 CONTINUE 3 CONTINUE C sweep up DO 13 I=2,IM-1 DO 11 J=2,JM-1 IF(FSM(I,J).EQ.0..OR.FSM(I,J+1).EQ.0.) GOTO 11 SL(I,J)=ABS(H(I,J+1)-H(I,J))/(H(I,J)+H(I,J+1)) IF(SL(I,J).LT.SLMIN) GOTO 11 DH=0.5*(SL(I,J)-SLMIN)*(H(I,J)+H(I,J+1)) SN=-1. IF(H(I,J+1).GT.H(I,J)) SN=1. H(I,J+1)=H(I,J+1)-SN*DH H(I,J)=H(I,J)+SN*DH 11 CONTINUE C sweep down DO 12 J=JM-1,2,-1 IF(FSM(I,J).EQ.0..OR.FSM(I,J+1).EQ.0.) GOTO 12 SL(I,J)=ABS(H(I,J+1)-H(I,J))/(H(I,J)+H(I,J+1)) IF(SL(I,J).LT.SLMIN) GOTO 12 DH=0.5*(SL(I,J)-SLMIN)*(H(I,J)+H(I,J+1)) SN=-1. IF(H(I,J+1).GT.H(I,J)) SN=1. H(I,J+1)=H(I,J+1)-SN*DH H(I,J)=H(I,J)+SN*DH 12 CONTINUE 13 CONTINUE C 100 CONTINUE RETURN END C SUBROUTINE PRXZ(LABEL,TIME,A,IM,ISKP,JM,J1,J2,KB,SCALA,DT,ZZ) DIMENSION A(IM,JM,KB),NUM(100),LINE(100),IDT(100), 1 ZZ(KB),DT(IM,JM) CHARACTER LABEL*(*) DATA ZERO/1.E-10/ C C THIS SUBROUTINE WRITES VERTICAL SECTIONS OF A 3-D FIELD C C TIME=TIME IN DAYS C A= ARRAY(IM,JM,KB) TO BE PRINTED C ISPL=PRINT SKIP FOR I C JSPL=PRINT SKIP FOR J C SCALE=DIVISOR FOR VALUES OF A C SCALE=SCALA IF (SCALE.GT.ZERO) GO TO 150 AMX=ZERO DO 140 K=1,KB DO 140 J=1,JM DO 140 I=1,IM,ISKP AMX=MAX(ABS(A(I,J,K)),AMX) 140 CONTINUE IF(AMX.EQ.0.) THEN SCALEI=0. GOTO 165 ENDIF SCALE=10.E0**(INT(LOG10(AMX)+1.E2)-103) 150 CONTINUE SCALEI=1./SCALE 165 CONTINUE WRITE(10,9980) LABEL 9980 FORMAT(' ',A30) WRITE(10,9981) TIME,SCALE 9981 FORMAT(' TIME = ',F9.3,' DAYS MULTIPLY ALL VALUES BY ',1PE8.2) C23456 | | | | | | | xxxxxxx| DO 10 JPP=1,2 J=J1 IF(JPP.EQ.2) J=J2 IF(J.EQ.0) RETURN WRITE(10,30) J 30 FORMAT(3X,/' SECTION J =',I3) IB=1 IE=IB+23*ISKP 50 CONTINUE IF(IE.GT.IM) IE=IM DO 100 I=IB,IE,ISKP IDT(I)=DT(I,J) 100 NUM(I)=I WRITE(10,9990) (NUM(I),I=IB,IE,ISKP) 9990 FORMAT(/,' I = ',24I5,/) WRITE(10,9991) (IDT(I),I=IB,IE,ISKP) 9991 FORMAT(8X,'D =',24I5.0,/,' ZZ ') DO 200 K=1,KB DO 190 I=IB,IE,ISKP 190 LINE(I)=SCALEI*A(I,J,K) WRITE(10,9966) K,ZZ(K),(LINE(I),I=IB,IE,ISKP) 9966 FORMAT(1X,I2,2X,F6.3,24I5) 200 CONTINUE WRITE(10,9984) IF(IE.EQ.IM) GO TO 10 IB=IB+24*ISKP IE=IB+23*ISKP GO TO 50 9984 FORMAT(//) 10 CONTINUE RETURN END C SUBROUTINE PRYZ(LABEL,TIME,A,IM,ISKP,JM,JSKP,KB,SCALA,DT,ZZ) DIMENSION A(IM,JM,KB),NUM(200),PLINE(200),DT(IM,JM),ZZ(KB) DIMENSION IP(2) CHARACTER LABEL*(*) DATA ZERO/1.E-10/ IP(1)=IM/2 IP(2)=IM-3 C C THIS SUBROUTINE WRITES VERTICAL SECTIONS OF A 3-D FIELD C C TIME=TIME IN DAYS C A= ARRAY(IM,JM,KB) TO BE PRINTED C ISPL=PRINT SKIP FOR I C JSPL=PRINT SKIP FOR J C SCALE=DIVISOR FOR VALUES OF A C SCALE=SCALA IF (SCALE.GT.ZERO) GO TO 150 AMX=ZERO DO 140 K=1,KB DO 140 J=1,JM,JSKP DO 140 I=1,IM,ISKP AMX=MAX(ABS(A(I,J,K)),AMX) 140 CONTINUE IF(AMX.EQ.0.) THEN SCALEI=0. GOTO 165 ENDIF SCALE=10.E0**(INT( LOG10(AMX)+1.E2)-103) 150 CONTINUE SCALEI=1./SCALE 165 CONTINUE WRITE(10,9980) LABEL 9980 FORMAT(' ',A30) WRITE(10,9981) TIME,SCALE 9981 FORMAT(' TIME = ',F8.2,' MULTIPLY ALL VALUES BY',E8.2) DO 10 IPP=1,2 I=IP(IPP) WRITE(10,30) I 30 FORMAT(3X,/' SECTION I =',I3) JB=1 JE=JB+23*JSKP 50 CONTINUE IF(JE.GT.JM) JE=JM DO 100 J=JB,JE,JSKP 100 NUM(J)=J WRITE(10,9990) (NUM(J),J=JB,JE,JSKP) 9990 FORMAT(/,' J = ',24I5,/) WRITE(10,9991) (DT(I,J),J=JB,JE,JSKP) 9991 FORMAT(7X,'D =',24F5.0,/,' ZZ ') DO 200 K=1,KB DO 190 J=JB,JE,JSKP 190 PLINE(J)=SCALEI*A(I,J,K) WRITE(10,9966) K,ZZ(K),(PLINE(J),J=JB,JE,JSKP) 9966 FORMAT(1X,I2,2X,F6.3,24(F5.1)) 200 CONTINUE WRITE(10,9984) IF(JE.EQ.JM) GO TO 10 JB=JB+24*JSKP JE=JB+23*JSKP GO TO 50 9984 FORMAT(//) 10 CONTINUE RETURN END c SUBROUTINE COUNT(HH,DD,MM,YY,DTI) C C COUNTER OF TIMING - UPDATE HOUR,DAY,MONTH,YEAR C THIS ROUTINE REPLACES ORIGINAL *COUNT* FROM C THE POMCFS, FIXING ITS ORIGINAL LEAP YEAR BUG C AND EXTENDING ITS VALIDITY BEYOND DEC. 31, 2003. C ALSO, A CHANGE WAS MADE TO FIT CONVENTIONS: C HOUR 24:00 IS CODED AS 0:00 THE NEXT DAY. C EVERY TIME THIS ROUTINE IS CALLED, TIME INCREMENT DTI C (SECONDS) IS ADDED, AND THEN DATE/TIME COMPONENTS C - HOUR OF THE DAY HH, DAY OF THE MONTH DD, MONTH OF C THE YEAR MM AND YEAR ARE ADJUSTED. C C PROGRAMMED AND TESTED: MARCH 28, 1996, LECH LOBOCKI C U.S. NATIONAL CENTERS FOR ENVIRONMENTAL PREDICTION C WD22LL@SUN1.WWB.NOAA.GOV TEL (301) 763-8133 C REAL MM C--------------------------------- advance clock HH=HH+DTI/3600. C--------------------------------- advance date when needed IF(HH.GE.24.) THEN HH=HH-24. CALL W3FS26 ( JDN(NINT(YY),NINT(MM),NINT(DD))+1, & IYEAR,MONTH,IDAY,IDAYWK,IDAYYR) YY=IYEAR MM=MONTH DD=IDAY ENDIF RETURN END C************************************************************* SUBROUTINE W3FS26(JLDAYN,IYEAR,MONTH,IDAY,IDAYWK,IDAYYR) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: W3FS26 YEAR, MONTH, DAY FROM JULIAN DAY NUMBER C AUTHOR: JONES,R.E. ORG: W342 DATE: 87-03-29 C C ABSTRACT: COMPUTES YEAR (4 DIGITS), MONTH, DAY, DAY OF WEEK, DAY C OF YEAR FROM JULIAN DAY NUMBER. THIS SUBROUTINE WILL WORK C FROM 1583 A.D. TO 3300 A.D. C C PROGRAM HISTORY LOG: C 87-03-29 R.E.JONES C 89-10-25 R.E.JONES CONVERT TO CRAY CFT77 FORTRAN C C USAGE: CALL W3FS26(JLDAYN,IYEAR,MONTH,IDAY,IDAYWK,IDAYYR) C C INPUT VARIABLES: C NAMES INTERFACE DESCRIPTION OF VARIABLES AND TYPES C ------ --------- ----------------------------------------------- C JLDAYN ARG LIST INTEGER JULIAN DAY NUMBER C C OUTPUT VARIABLES: C NAMES INTERFACE DESCRIPTION OF VARIABLES AND TYPES C ------ --------- ----------------------------------------------- C IYEAR ARG LIST INTEGER YEAR (4 DIGITS) C MONTH ARG LIST INTEGER MONTH C IDAY ARG LIST INTEGER DAY C IDAYWK ARG LIST INTEGER DAY OF WEEK (1 IS SUNDAY, 7 IS SAT) C IDAYYR ARG LIST INTEGER DAY OF YEAR (1 TO 366) C C REMARKS: A JULIAN DAY NUMBER CAN BE COMPUTED BY USING ONE OF THE C FOLLOWING STATEMENT FUNCTIONS. A DAY OF WEEK CAN BE COMPUTED C FROM THE JULIAN DAY NUMBER. A DAY OF YEAR CAN BE COMPUTED FROM C A JULIAN DAY NUMBER AND YEAR. C C IYEAR (4 DIGITS) C C JDN(IYEAR,MONTH,IDAY) = IDAY - 32075 C & + 1461 * (IYEAR + 4800 + (MONTH - 14) / 12) / 4 C & + 367 * (MONTH - 2 - (MONTH -14) / 12 * 12) / 12 C & - 3 * ((IYEAR + 4900 + (MONTH - 14) / 12) / 100) / 4 C C IYR (4 DIGITS) , IDYR(1-366) DAY OF YEAR C C JULIAN(IYR,IDYR) = -31739 + 1461 * (IYR + 4799) / 4 C & -3 * ((IYR + 4899) / 100) / 4 + IDYR C C DAY OF WEEK FROM JULIAN DAY NUMBER, 1 IS SUNDAY, 7 IS SATURDAY. C C JDAYWK(JLDAYN) = MOD((JLDAYN + 1),7) + 1 C C DAY OF YEAR FROM JULIAN DAY NUMBER AND 4 DIGIT YEAR. C C JDAYYR(JLDAYN,IYEAR) = JLDAYN - C & (-31739+1461*(IYEAR+4799)/4-3*((IYEAR+4899)/100)/4) C C THE FIRST FUNCTION WAS IN A LETTER TO THE EDITOR COMMUNICATIONS C OF THE ACM VOLUME 11 / NUMBER 10 / OCTOBER, 1968. THE 2ND C FUNCTION WAS DERIVED FROM THE FIRST. THIS SUBROUTINE WAS ALSO C INCLUDED IN THE SAME LETTER. JULIAN DAY NUMBER 1 IS C JAN 1,4713 B.C. A JULIAN DAY NUMBER CAN BE USED TO REPLACE A C DAY OF CENTURY, THIS WILL TAKE CARE OF THE DATE PROBLEM IN C THE YEAR 2000, OR REDUCE PROGRAM CHANGES TO ONE LINE CHANGE C OF 1900 TO 2000. JULIAN DAY NUMBERS CAN BE USED FOR FINDING C RECORD NUMBERS IN AN ARCHIVE OR DAY OF WEEK, OR DAY OF YEAR. C C ATTRIBUTES: C LANGUAGE: CRAY CFT77 FORTRAN C MACHINE: CRAY Y-MP8/864 C C$$$ C L = JLDAYN + 68569 N = 4 * L / 146097 L = L - (146097 * N + 3) / 4 I = 4000 * (L + 1) / 1461001 L = L - 1461 * I / 4 + 31 J = 80 * L / 2447 IDAY = L - 2447 * J / 80 L = J / 11 MONTH = J + 2 - 12 * L IYEAR = 100 * (N - 49) + I + L IDAYWK = MOD((JLDAYN + 1),7) + 1 IDAYYR = JLDAYN - & (-31739 +1461 * (IYEAR+4799) / 4 - 3 * ((IYEAR+4899)/100)/4) RETURN END C************************************************************* INTEGER FUNCTION JDN(IYEAR,MONTH,IDAY) JDN = IDAY - 32075 & + 1461 * (IYEAR + 4800 + (MONTH - 14) / 12) / 4 & + 367 * (MONTH - 2 - (MONTH -14) / 12 * 12) / 12 & - 3 * ((IYEAR + 4900 + (MONTH - 14) / 12) / 100) / 4 RETURN END C************************************************************* FUNCTION XXTIME(hour,day,rmonth,year) C THIS ROUTINE REPLACES ORIGINAL *XXTIME* FROM C THE POMCFS, BRINGING A CONSISTENT USE OF W3FS26. C PROGRAMMED AND TESTED: MARCH 28, 1996, LECH LOBOCKI C U.S. NATIONAL CENTERS FOR ENVIRONMENTAL PREDICTION C WD22LL@SUN1.WWB.NOAA.GOV TEL (301) 763-8133 ndays=JDN(NINT(year),NINT(rmonth),NINT(day))-JDN(0,1,1) XXTIME=ndays*24+hour return end C************************************************************* C+++++++++++++++++++++++ 8/2/1998 SUBROUTINE SAVELOLD(NUNIT,hh,dd,mm,yy,EL,IM,JM) C C SAVE SEA LEVEL AT SELECTED STATIONS. C REAL MM PARAMETER(NSTA=22) DIMENSION IST(NSTA),JST(NSTA),DUM(NSTA),ELSTA(NSTA) DIMENSION EL(IM,JM) CHARACTER*36 STATION(NSTA) C ORIGINAL SETUP AS RECEIVED FROM PRINCETON DATA (STATION(I),IST(I),JST(I),DUM(I),DUM(I),I=1,NSTA) / C STATION I J SLAT SLON 1"FORT PULASKI, SAVANNAH RIVER, GA ", 154,157,-80.655,32.053, 2"SPRINGMAID PIER, SC ", 181,170,-78.731,33.656, 3"WILMINGTON, NC ", 193,172,-77.752,34.237, 4"DUCK, FRF PIER, NC ", 222,186,-75.689,36.169, 5"CHESAPEAKE BAY BRIDGE TUNNEL, VA ", 226,200,-75.849,36.945, 6"LEWES, FT. MILES, DE ", 246,221,-74.969,38.769, 7"ATLANTIC CITY, ATLANTIC OCEAN, NJ ", 255,223,-74.316,39.338, 8"SANDY HOOK, NJ ", 267,235,-73.876,40.472, 9"MONTAUK, FORT POND BAY, NY ", 284,225,-71.950,41.072, &"NEWPORT, NARRAGANSETT BAY, RI ", 291,225,-71.318,41.486, 1"BOSTON, BOSTON HARBOR, MA ", 300,232,-70.863,42.357, 2"PORTLAND, CASCO BAY, ME ", 313,241,-70.165,43.656, 3"CLEARWATER BEACH, FL ", 112,110,-83.184,27.987, 4"ST. PETERSBURG, FL ", 114,106,-82.902,27.721, 5"NAPLES, GULF OF MEXICO, FL ", 119, 87,-81.999,26.139, 6"KEY WEST, FL ", 119, 72,-81.710,24.526, 7"FREEPORT, DOW BARGE CANAL, TX ", 18,179,-95.087,28.942, 8"CORPUS CHRISTI, GULF OF MEXICO, TX", 5,143,-97.054,27.588, 9"GALVESTON PLEASURE PIER, TX ", 22,189,-94.574,29.295, &"SABINE PASS NORTH, TX ", 30,194,-93.518,29.623, 1"BURMUDA ", 303,107,-64.666,32.364, 2"PANAMA CITY ", 93,157,-85.957,30.204/ CLHC for new grid C STATION I J H SLAT SLON c 1"Lake Worth, FL ", 2, 2, 433.00, 26.61, -80.03, c 2"Settlement Pt. Bahama", 6, 1, 350.00, 26.68, -79.00, c 3"St. Augustine Beach, ", 2, 40, 23.00, 29.86, -81.26, c 4"Fort Pulaski, GA ", 8, 63, 12.00, 32.03, -80.90, c 5"Bermuda, ESSO ", 91, 21, 297.00, 32.37, -64.70, c 6"Springmaid Pier, SC ", 23, 74, 10.00, 33.65, -78.92, c 7"Wilmington, NC ", 30, 74, 13.00, 34.23, -77.95, c 8"Hatteras, NC ", 44, 73, 11.00, 35.22, -75.63, c 9"Duck, NC ", 47, 81, 16.00, 36.18, -75.75, c *"Bay Bridge Tunnel, V ", 50, 90, 15.00, 36.97, -76.11, c 1"Lewes, DE ", 62, 96, 19.00, 38.78, -75.13, c 2"Atlantic City, NJ ", 67, 95, 19.00, 39.35, -74.42, c 3"Sandy Hook, NJ ", 73, 99, 11.00, 40.47, -74.00, c 4"Montauk, NY ", 81, 95, 44.00, 41.05, -71.96, c 5"Newport, RI ", 84, 94, 16.00, 41.50, -71.33, c 6"Boston, MA ", 89, 97, 28.00, 42.35, -71.05, c 7"Portland, ME ", 95, 100, 87.00, 43.66, -70.25, c 8"Halifax ", 121, 94, 54.00, 44.67, -63.58, c 9"Eastport, ME ", 109, 100, 55.00, 44.90, -66.98, c *"St. John's Newfoundl ", 168, 100, 99.00, 47.57, -52.70/ C LL_960327 ** WARNING: NOTE THAT Boston, MA has i=87 now C LL_960510 ** THE LATTER WAS WRONG, IT HAS TO BE i=89 C OLD BATHYMETRY ACCORDING TO EUGENE WEI C DATA (STATION(I),IST(I),JST(I),DUM(I),DUM(I),DUM(I),I=1,NSTA) / C STATION I J H SLAT SLON C 1 "Lake Worth, FL ", 2, 2, 337.70, 26.61, -80.03, C 2 "Settlement Pt. Bahamas ", 6, 1, 349.82, 26.68, -79.00, C 3 "St. Augustine Beach, FL", 2, 40, 11.76, 29.86, -81.26, C 4 "Fort Pulaski, GA ", 9, 61, 12.02, 32.03, -80.90, C 5 "Bermuda, ESSO ", 91, 21, 296.76, 32.37, -64.70, C 6 "Springmaid Pierr, SC ", 23, 70, 14.65, 33.65, -78.92, C 7 "Wilmington, NC ", 29, 71, 14.62, 34.23, -77.95, C 8 "Hatteras, NC ", 43, 72, 11.34, 35.22, -75.63, C 9 "Duck, NC ", 48, 82, 13.84, 36.18, -75.75, C & "Bay Bridge Tunnel, VA ", 51, 88, 12.12, 36.97, -76.11, C 1 "Lewes, DE ", 62, 94, 14.15, 38.78, -75.13, C 2 "Atlantic City, NJ ", 67, 94, 28.21, 39.35, -74.42, C 3 "Sandy Hook, NJ ", 73, 98, 25.99, 40.47, -74.00, C 4 "Montauk, NY ", 81, 95, 15.10, 41.05, -71.96, C 5 "Newport, RI ", 84, 94, 26.84, 41.50, -71.33, C 6 "Boston, MA ", 89, 97, 22.16, 42.35, -71.05, C 7 "Portland, ME ", 95, 100, 96.99, 43.66, -70.25, C 8 "Halifax ",121, 94, 54.43, 44.67, -63.58, C 9 "Eastport, ME ",109, 100, 54.82, 44.90, -66.98, C & "St. John's Newfoundland",169, 100, 120.17, 47.57, -52.70 / C NEW (NOS 15) BATHYMETRY ACCORDING TO EUGENE WEI C DATA (STATION(I),IST(I),JST(I),DUM(I),DUM(I),DUM(I),I=1,NSTA) / C STATION I J H SLAT SLON C 1 "Lake Worth, FL ", 2, 2, 337.70, 26.61, -80.03, C 2 "Settlement Pt. Bahamas ", 6, 1, 349.82, 26.68, -79.00, C 3 "St. Augustine Beach, FL", 2, 39, 11.76, 29.86, -81.26, C 4 "Fort Pulaski, GA ", 8, 63, 12.02, 32.03, -80.90, C 5 "Bermuda, ESSO ", 91, 20, 296.76, 32.37, -64.70, C 6 "Springmaid Pierr, SC ", 22, 73, 14.65, 33.65, -78.92, C 7 "Wilmington, NC ", 29, 73, 14.62, 34.23, -77.95, C 8 "Hatteras, NC ", 43, 72, 11.34, 35.22, -75.63, C 9 "Duck, NC ", 47, 81, 13.84, 36.18, -75.75, C & "Bay Bridge Tunnel, VA ", 50, 88, 12.12, 36.97, -76.11, C 1 "Lewes, DE ", 62, 95, 14.15, 38.78, -75.13, C 2 "Atlantic City, NJ ", 66, 94, 28.21, 39.35, -74.42, C 3 "Sandy Hook, NJ ", 72, 98, 25.99, 40.47, -74.00, C 4 "Montauk, NY ", 80, 94, 15.10, 41.05, -71.96, C 5 "Newport, RI ", 84, 94, 26.84, 41.50, -71.33, C 6 "Boston, MA ", 89, 97, 22.16, 42.35, -71.05, C 7 "Portland, ME ", 94, 100, 96.99, 43.66, -70.25, C 8 "Halifax ",121, 94, 54.43, 44.67, -63.58, C 9 "Eastport, ME ",108, 100, 54.82, 44.90, -66.98, C & "St. John's Newfoundland",168, 100, 120.17, 47.57, -52.70/ DO I=1,NSTA ELSTA(I)=EL(IST(I),JST(I)) ENDDO C WRITE(NUNIT,701) hh,dd,mm,yy,(ELSTA(I),I=1,NSTA) RETURN c 701 FORMAT (4f5.0,22(1X,F7.4)) 701 FORMAT (f7.3,3f5.0,22(1X,F7.4)) END C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE SAVEL(NUNIT,hh,dd,mm,yy,EL,IM,JM) C C SAVE SEA LEVEL AT SELECTED STATIONS. C REAL MM PARAMETER(NSTA=21) DIMENSION IST(NSTA),JST(NSTA),DUM(NSTA),ELSTA(NSTA) DIMENSION EL(IM,JM) CHARACTER*36 STATION(NSTA) C ORIGINAL SETUP AS RECEIVED FROM PRINCETON DATA (STATION(I),IST(I),JST(I),DUM(I),DUM(I),I=1,NSTA) / C STATION I J SLAT SLON 1"BURMUDA ", 303,107,-64.666,32.364, 2"PORTLAND, CASCO BAY, ME ", 313,241,-70.165,43.656, 3"BOSTON, BOSTON HARBOR, MA ", 300,232,-70.863,42.357, 4"NEWPORT, NARRAGANSETT BAY, RI ", 291,225,-71.318,41.486, 5"MONTAUK, FORT POND BAY, NY ", 284,225,-71.950,41.072, 6"SANDY HOOK, NJ ", 267,235,-73.876,40.472, 7"ATLANTIC CITY, ATLANTIC OCEAN, NJ ", 255,223,-74.316,39.338, 8"LEWES, FT. MILES, DE ", 246,221,-74.969,38.769, 9"CHESAPEAKE BAY BRIDGE TUNNEL, VA ", 226,200,-75.849,36.945, &"DUCK, FRF PIER, NC ", 222,186,-75.689,36.169, 1"WILMINGTON, NC ", 193,172,-77.752,34.237, 2"SPRINGMAID PIER, SC ", 181,170,-78.731,33.656, 3"FORT PULASKI, SAVANNAH RIVER, GA ", 154,157,-80.655,32.053, 4"PORT CANAVERAL, FL ", 141,107,-80.593,28.415, 5"VIRGINIA KEY, FL ", 138, 80,-80.162,25.737, 6"KEY WEST, FL ", 119, 72,-81.710,24.526, 7"NAPLES, GULF OF MEXICO, FL ", 119, 87,-81.999,26.139, 8"ST. PETERSBURG, FL ", 114,106,-82.902,27.721, 9"PANAMA CITY ", 93,157,-85.957,30.204, &"GALVESTON PLEASURE PIER, TX ", 22,189,-94.574,29.295, 1"CORPUS CHRISTI, GULF OF MEXICO, TX", 5,143,-97.054,27.588/ CLHC for new grid C STATION I J H SLAT SLON c 1"Lake Worth, FL ", 2, 2, 433.00, 26.61, -80.03, c 2"Settlement Pt. Bahama", 6, 1, 350.00, 26.68, -79.00, c 3"St. Augustine Beach, ", 2, 40, 23.00, 29.86, -81.26, c 4"Fort Pulaski, GA ", 8, 63, 12.00, 32.03, -80.90, c 5"Bermuda, ESSO ", 91, 21, 297.00, 32.37, -64.70, c 6"Springmaid Pier, SC ", 23, 74, 10.00, 33.65, -78.92, c 7"Wilmington, NC ", 30, 74, 13.00, 34.23, -77.95, c 8"Hatteras, NC ", 44, 73, 11.00, 35.22, -75.63, c 9"Duck, NC ", 47, 81, 16.00, 36.18, -75.75, c *"Bay Bridge Tunnel, V ", 50, 90, 15.00, 36.97, -76.11, c 1"Lewes, DE ", 62, 96, 19.00, 38.78, -75.13, c 2"Atlantic City, NJ ", 67, 95, 19.00, 39.35, -74.42, c 3"Sandy Hook, NJ ", 73, 99, 11.00, 40.47, -74.00, c 4"Montauk, NY ", 81, 95, 44.00, 41.05, -71.96, c 5"Newport, RI ", 84, 94, 16.00, 41.50, -71.33, c 6"Boston, MA ", 89, 97, 28.00, 42.35, -71.05, c 7"Portland, ME ", 95, 100, 87.00, 43.66, -70.25, c 8"Halifax ", 121, 94, 54.00, 44.67, -63.58, c 9"Eastport, ME ", 109, 100, 55.00, 44.90, -66.98, c *"St. John's Newfoundl ", 168, 100, 99.00, 47.57, -52.70/ C LL_960327 ** WARNING: NOTE THAT Boston, MA has i=87 now C LL_960510 ** THE LATTER WAS WRONG, IT HAS TO BE i=89 C OLD BATHYMETRY ACCORDING TO EUGENE WEI C DATA (STATION(I),IST(I),JST(I),DUM(I),DUM(I),DUM(I),I=1,NSTA) / C STATION I J H SLAT SLON C 1 "Lake Worth, FL ", 2, 2, 337.70, 26.61, -80.03, C 2 "Settlement Pt. Bahamas ", 6, 1, 349.82, 26.68, -79.00, C 3 "St. Augustine Beach, FL", 2, 40, 11.76, 29.86, -81.26, C 4 "Fort Pulaski, GA ", 9, 61, 12.02, 32.03, -80.90, C 5 "Bermuda, ESSO ", 91, 21, 296.76, 32.37, -64.70, C 6 "Springmaid Pierr, SC ", 23, 70, 14.65, 33.65, -78.92, C 7 "Wilmington, NC ", 29, 71, 14.62, 34.23, -77.95, C 8 "Hatteras, NC ", 43, 72, 11.34, 35.22, -75.63, C 9 "Duck, NC ", 48, 82, 13.84, 36.18, -75.75, C & "Bay Bridge Tunnel, VA ", 51, 88, 12.12, 36.97, -76.11, C 1 "Lewes, DE ", 62, 94, 14.15, 38.78, -75.13, C 2 "Atlantic City, NJ ", 67, 94, 28.21, 39.35, -74.42, C 3 "Sandy Hook, NJ ", 73, 98, 25.99, 40.47, -74.00, C 4 "Montauk, NY ", 81, 95, 15.10, 41.05, -71.96, C 5 "Newport, RI ", 84, 94, 26.84, 41.50, -71.33, C 6 "Boston, MA ", 89, 97, 22.16, 42.35, -71.05, C 7 "Portland, ME ", 95, 100, 96.99, 43.66, -70.25, C 8 "Halifax ",121, 94, 54.43, 44.67, -63.58, C 9 "Eastport, ME ",109, 100, 54.82, 44.90, -66.98, C & "St. John's Newfoundland",169, 100, 120.17, 47.57, -52.70 / C NEW (NOS 15) BATHYMETRY ACCORDING TO EUGENE WEI C DATA (STATION(I),IST(I),JST(I),DUM(I),DUM(I),DUM(I),I=1,NSTA) / C STATION I J H SLAT SLON C 1 "Lake Worth, FL ", 2, 2, 337.70, 26.61, -80.03, C 2 "Settlement Pt. Bahamas ", 6, 1, 349.82, 26.68, -79.00, C 3 "St. Augustine Beach, FL", 2, 39, 11.76, 29.86, -81.26, C 4 "Fort Pulaski, GA ", 8, 63, 12.02, 32.03, -80.90, C 5 "Bermuda, ESSO ", 91, 20, 296.76, 32.37, -64.70, C 6 "Springmaid Pierr, SC ", 22, 73, 14.65, 33.65, -78.92, C 7 "Wilmington, NC ", 29, 73, 14.62, 34.23, -77.95, C 8 "Hatteras, NC ", 43, 72, 11.34, 35.22, -75.63, C 9 "Duck, NC ", 47, 81, 13.84, 36.18, -75.75, C & "Bay Bridge Tunnel, VA ", 50, 88, 12.12, 36.97, -76.11, C 1 "Lewes, DE ", 62, 95, 14.15, 38.78, -75.13, C 2 "Atlantic City, NJ ", 66, 94, 28.21, 39.35, -74.42, C 3 "Sandy Hook, NJ ", 72, 98, 25.99, 40.47, -74.00, C 4 "Montauk, NY ", 80, 94, 15.10, 41.05, -71.96, C 5 "Newport, RI ", 84, 94, 26.84, 41.50, -71.33, C 6 "Boston, MA ", 89, 97, 22.16, 42.35, -71.05, C 7 "Portland, ME ", 94, 100, 96.99, 43.66, -70.25, C 8 "Halifax ",121, 94, 54.43, 44.67, -63.58, C 9 "Eastport, ME ",108, 100, 54.82, 44.90, -66.98, C & "St. John's Newfoundland",168, 100, 120.17, 47.57, -52.70/ DO I=1,NSTA ELSTA(I)=EL(IST(I),JST(I)) ENDDO C WRITE(NUNIT,701) hh,dd,mm,yy,(ELSTA(I),I=1,NSTA) RETURN c 701 FORMAT (4f5.0,21(1X,F7.4)) 701 FORMAT (f7.3,3f5.0,21(1X,F7.4)) END C c C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE PRXY (LABEL,TIME,A,IM,ISKP,JM,JSKP,SCALA) C >>> C THIS WRITES A 2-D FIELD C TIME=TIME IN DAYS C A = ARRAY(IM,JM) TO BE PRINTED C ISKP=PRINT SKIP FOR I C JSKP=PRINT SKIP FOR J C SCALE=DIVISOR FOR VALUES OF A C C IMPLICIT HALF PRECISION (A-H,O-Z) DIMENSION A(IM,JM),NUM(IM),LINE(IM) CHARACTER LABEL*(*) DATA ZERO /1.E-12/ C SCALE=SCALA IF (SCALE.GT.ZERO) GO TO 160 AMX=ZERO DO 150 J=1,JM,JSKP DO 150 I=1,IM,ISKP AMX=MAX(ABS(A(I,J)),AMX) 150 CONTINUE IF(AMX.EQ.0.) THEN SCALEI=0. GOTO 165 ENDIF SCALE=10.E0**(INT(LOG10(AMX)+1.E2)-103) 160 CONTINUE SCALEI=1.E0/SCALE 165 CONTINUE WRITE(10,170) LABEL 170 FORMAT(1X,A40) WRITE(10,180) TIME,SCALE 180 FORMAT(' TIME =',F9.4,' DAYS MULTIPLY ALL VALUES BY',1PE9.2) DO 190 I=1,IM 190 NUM(I)=I IB=1 C 200 CONTINUE IE=IB+23*ISKP IF(IE.GT.IM) IE=IM WRITE(10,210) (NUM(I),I=IB,IE,ISKP) 210 FORMAT(/,2X,24I5,/) DO 260 J=1,JM,JSKP JWR=JM+1-J DO 220 I=IB,IE,ISKP 220 LINE(I)=INT(SCALEI*A(I,JWR)) WRITE(10,240) JWR,(LINE(I),I=IB,IE,ISKP) 240 FORMAT(1X,I3,24I5) 260 CONTINUE WRITE(10,280) 280 FORMAT(//) IF(IE.GE.IM) RETURN IB=IB+24*ISKP GO TO 200 END SUBROUTINE TIDE_BODY(hh,dd,mm,yy,DTE,IEXT,EEQ,ELF,ALAT,ALON) C C astronomical tidal force. C PARAMETER (IM=401,JM=251) REAL EEQ(IM,JM),ELF(IM,JM),ALAT(IM,JM),ALON(IM,JM), > A_EEQ(IM,JM,6),B_EEQ(IM,JM,6) REAL MASK_T COMMON/TIDE1/AMP(6),OME(6),CHI(6),MASK_T(6) ! initiated in TIDE_SETUP REAL dd0,mm0,yy0,hh,dd,mm,yy DATA RAD/0.0174533/,PI/3.14159/ DATA NCALL/0/ SAVE C NCALL=NCALL+1 C C CALL DMN2JD (DD, MM, YY, D) C TIS=D+hh/24.0+IEXT*DTE/86400. D=(XXTIME(hh,dd,mm,yy)-XXTIME(0.,1.,1.,yy))/24.0 ! julian day - 1 c TIS=D+IEXT*DTE/86400. TIS=(float(IINT*isplit+IEXT)*dte)/86400.0 C WRITE(6,'(/" TIDE_BODY CALL:",I5," AT TIME:",1pe12.5)') NCALL,TIS IF (YY .NE. YY0) THEN ! CORRECTION, update CHI yearly C-----parameters ------------------------------- iday=1 ! CORRECTION iyear=yy inty=(iyear-73)/4 db=iday+365*(iyear-75)+inty t=(27392.500528+1.000000035*db)/36525 h0=279.69668+36000.768930485*t+3.03e-4*t**2 s0=270.434358+481267.88314137*t-0.001133*t**2+1.9e-6*t**3 p0=334.329653+4069.0340329575*t-0.010325*t**2-1.2e-5*t**3 C-----astronomical arguments-------------------- CHI(1)=(2*h0-2*s0)*rad ! M2 CHI(2)=0. ! S2 CHI(3)=(2*h0-3*s0+p0)*rad ! N2 CHI(4)=(h0+90)*rad ! K1 CHI(5)=(h0-2*s0-90)*rad ! O1 CHI(6)=-(h0+90)*rad ! P1 ! CORRECTION C*********************** do n=1,6 chi(n)=0.0 enddo C*********************** DO N=1,3 A_EEQ(:,:,N)= 1 AMP(N)*(COS(RAD*ALAT(:,:)))**2*COS(RAD*2.*ALON(:,:)) B_EEQ(:,:,N)= 1 AMP(N)*(COS(RAD*ALAT(:,:)))**2*SIN(RAD*2.*ALON(:,:)) M=N+3 A_EEQ(:,:,M)= 1 AMP(M)*SIN(2.*RAD*ALAT(:,:))*COS(RAD*ALON(:,:)) B_EEQ(:,:,M)= 1 AMP(M)*SIN(2.*RAD*ALAT(:,:))*SIN(RAD*ALON(:,:)) ENDDO dd0=dd mm0=mm yy0=yy c WRITE(6,'(" astronomical arguments updated!")') c WRITE(6,'(" CHI : ",6E12.4)') CHI c WRITE(6,'(" AMP : ",6F8.4)') AMP c WRITE(6,'(" OME : ",6F8.4)') OME c WRITE(6,'(" MASK: ",6F8.4)') MASK_T ENDIF EEQ(:,:)=0. DO N=1,6 EEQ(:,:)=EEQ(:,:)+ 1 MASK_T(N)* ( A_EEQ(:,:,N)*COS(OME(N)*TIS+CHI(N)) 1 -B_EEQ(:,:,N)*SIN(OME(N)*TIS+CHI(N)) ) ENDDO EEQ=0.69*EEQ C+++++++++++++ 5/26/98 C EEQ=EEQ-0.10*ELF EEQ=EEQ+0.10*ELF C++++++++++++++++ RETURN END SUBROUTINE TIDE_SETUP PARAMETER (IM=401,JM=251,NC=6, LJN=JM*NC,LIN=IM*NC) REAL MASK_T COMMON/TIDE1/AMP(6),OME(6),CHI(6),MASK_T(6) COMMON/TIDE2/EAMPE(JM,NC),EPHAE(JM,NC), 1 UAMPE(JM,NC),UPHAE(JM,NC), 2 VAMPE(JM,NC),VPHAE(JM,NC) DATA EAMPE/LJN*0.0/,UAMPE/LJN*0.0/,VAMPE/LJN*0.0/, 1 EPHAE/LJN*0.0/,UPHAE/LJN*0.0/,VPHAE/LJN*0.0/ C PC 07/10/96 C AMP_CO and PHA_CO are used TEMPORARILY to fix diurnal tidal amplitude and phase; DATA AMP_CO/0.45/,PHA_CO/-100.0/ c c open(33,file='/t90/lhc/Atl/Code/Tide/eamphe.dat' c open(33,file='/t90/lhc/Atl/Code/Tide/tidebc2.dat' c $ ,form='FORMATTED') c C---------------- M2 S2 N2 K1 O1 P1--------------------- DATA MASK_T/1., 0., 0., 0., 0., 0./ ! =1, tide included DATA AMP/0.242334,0.112841,0.046398, ! unit: m > 0.141565,0.100514,0.046843/, > OME/1.40519e-4,1.45444E-4,1.37880E-4, ! unit: 1/s > 0.72921e-4,0.67598e-4,0.72523e-4/ OME(:)=OME(:)*86400.0 ! unit: 1/day DO N=1,NC C------------------ At first test all zero, No letteral forcing do j=1,jm EAMPE(j,n)=0.0 EPHAE(j,n)=0.0 UAMPE(j,n)=0.0 UPHAE(j,n)=0.0 VAMPE(j,n)=0.0 VPHAE(j,n)=0.0 C c READ(7,'(i5,6f10.4)') (ii,EAMPE(j,n),EPHAE(j,n), c > UAMPE(j,n),UPHAE(j,n),VAMPE(j,n),VPHAE(j,n),j=1,jm) enddo ENDDO C--------------------------------------------------------------------- do j=1,JM read(33,'(i5,2f10.3)') ii,eampe(j,1),ephae(j,1) eampe(j,1)=eampe(j,1)*0.01 enddo C--------------------------------------------------------------------- EAMPE(:,4:6)=AMP_CO*EAMPE(:,4:6) UAMPE(:,4:6)=AMP_CO*UAMPE(:,4:6) VAMPE(:,4:6)=AMP_CO*VAMPE(:,4:6) EPHAE(:,4:6)=EPHAE(:,4:6)+PHA_CO UPHAE(:,4:6)=UPHAE(:,4:6)+PHA_CO VPHAE(:,4:6)=VPHAE(:,4:6)+PHA_CO c DO N=1,NC c WRITE(6,'(/" TIDAL CONSTITUENT: ",I2)') N c WRITE(6,'(i5,6f10.4)') (i,EAMPS(i,n),EPHAS(i,n), c > UAMPS(i,n),UPHAS(i,n),VAMPS(i,n),VPHAS(i,n),i=1,im) c WRITE(6,'(i5,6f10.4)') (j,EAMPE(j,n),EPHAE(j,n), c > UAMPE(j,n),UPHAE(j,n),VAMPE(j,n),VPHAE(j,n),j=1,jm) c ENDDO RETURN END SUBROUTINE PRXYM(LABEL,TIME,A,IM,IP,JM,JP,SCALA) C >>> C THIS WRITES A 2-D FIELD C TIME=TIME IN DAYS C A = ARRAY(IM,JM) TO BE PRINTED C ISKP=PRINT SKIP FOR I C JSKP=PRINT SKIP FOR J C SCALE=DIVISOR FOR VALUES OF A C C IMPLICIT HALF PRECISION (A-H,O-Z) DIMENSION A(IM,JM),NUM(IM),LINE(IM) CHARACTER LABEL*(*) DATA ZERO /1.E-12/ C SCALE=SCALA IF (SCALE.GT.ZERO) GO TO 160 AMX=ZERO DO 150 J=JP,JM DO 150 I=IP,IM AMX=MAX(ABS(A(I,J)),AMX) 150 CONTINUE IF(AMX.EQ.0.) THEN SCALEI=0. GOTO 165 ENDIF SCALE=10.E0**(INT(LOG10(AMX)+1.E2)-103) 160 CONTINUE SCALEI=1.E0/SCALE 165 CONTINUE WRITE(10,170) LABEL 170 FORMAT(1X,A40) WRITE(10,180) TIME,SCALE 180 FORMAT(' TIME =',F9.4,' DAYS MULTIPLY ALL VALUES BY',1PE9.2) DO 190 I=1,IM 190 NUM(I)=I c iskp=1 jskp=1 ib=ip c IB=1 C 200 CONTINUE IE=IB+23*ISKP IF(IE.GT.IM) IE=IM WRITE(10,210) (NUM(I),I=IB,IE,ISKP) 210 FORMAT(/,2X,24I5,/) DO 260 J=jp,JM,JSKP JWR=JM+1-J DO 220 I=IB,IE,ISKP 220 LINE(I)=INT(SCALEI*A(I,JWR)) WRITE(10,240) JWR,(LINE(I),I=IB,IE,ISKP) 240 FORMAT(1X,I3,24I5) 260 CONTINUE WRITE(10,280) 280 FORMAT(//) IF(IE.GE.IM) RETURN IB=IB+24*ISKP GO TO 200 END