BLOCK DATA C JELESNIANSKI/CHEN JULY 1988 TDL FORTRAN 77 VERSION C C PURPOSE C DATA INITIALIZATION C NONE C C VARIABLES C IP(4) JP(4) = SHIFT SUBSCRIPTS T 4 ADJACENT MOMENTUM POINTS C IH(4) JH(4) = SHIFT SUBSCRIPTS TO 4 ADJACENT SURGE POINTS ABOUT C A MOMUMTUM POINT C IIH(4) JJH(4) = SHIFT SUBSCRIPTS TO 4 ADJACENT SURGE POINTS ABOUT C A SURGE POINT C IOPERL = 2 AS DEFAULT IN OPERATIONAL MODE C MAPIN(HR) = SPECIFIES THE PRINT HOUR, BEFORE AND AFTER THE C NEAREST APPROACH OF STORM, FOR SNAPSHOT OUTPUT C G = GRAVITY C C7 = SLIP COEFFICIENT C C17 = AIR DENSITY C C19 = FRICTIONAL COEFFICIENT C C25 = EDDY VISCOSITY C MONTH(12) = ACCUMULATED DAYS AT MONTHLY INCREMENTS C AMON(12) = FIRST 3 CHARACTERS FOR EACH MONTH C C GENERAL COMMENTS C THIS MEMBER RESIDES IN THE ROOT OVERLAY C COMMON /DUMB8/ IP(4),JP(4),IH(4),JH(4) COMMON /DUM88/ IIH(4),JJH(4),IHH(4),JHH(4) COMMON /SWTCH/ IOPERL(5) COMMON /DELH11/ MAPIN(10),NUMBER COMMON /EGTH/ DELS,DELT,G,COR COMMON /THRD/ C7,C17,C19,C25 COMMON /MTH/ MONTH(12),AMON(12) CHARACTER*4 AMON C DATA IP/0,1,1,0/,JP/0,0,1,1/,IH/-1,0,-1,0/,JH/-1,-1,0,0/ DATA IIH/0,1,0,-1/,JJH/-1,0,1,0/,IHH/-1,1,0,0/,JHH/0,0,-1,1/ DATA IOPERL/2,2,2,2,2/ DATA MAPIN/-3,0,2,7*100/ DATA DELS,G/5280.,32.2/ DATA C7,C17,C19,C25 / .006 , .00115 , 3.E-6,.25 / DATA MONTH /-1,30,58,89,119,150,180,211,242,272,303,333/ DATA AMON/'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP', 1 'OCT','NOV','DEC' / C SET DEFAULTS TO NONARCHIVING MODE ( NOARCH='X'). NO TIME HISTRY C KHSPT=0. COMMON /ARCH/ NOARCH COMMON /NFILE/ KNPT,NDATA(3) COMMON /CESAV/ IPN(250),IPL(250),KHSPT,IPN2(250),IPL2(250), $ KHSPT2 CHARACTER*1 NOARCH DATA KNPT/0/,KHSPT/0/,NDATA/13,14,15/,NOARCH/'X'/ C THIS BLOCK DATA IS TO BE USED FOR LAKE PONTCHARTRAIN C BASIN ONLY. c INCLUDE 'parm.for' C INCLUDE 'parmmsy.for' C DATA ZSUBCE/-1./,JSUB,JSUB1/38,50/ DATA LDTMG/5*33,8*35,5*34, 1 30,29,29,28,27,26,26,24,22,20,21,22,22,21, 2 68*1,899*1/ DATA NODRY/74/ DATA IDRY/21,22,22,23,23,23,23,24,24,24,24,25,25,25,25,26,26,26, 1 27,27,27,28,28,28,29,29,29,29,30,30,30,30,30,30,31,31, 2 31,31,28,28,28,28,28,29,29,29,29,30,30,30,30,30,31,31, 3 31,32,32,32,32,32,32,32,33,33,33,33,33,34,34,34,34,34, 4 34,34,5926*1/ DATA JDRY/28,28,29,27,28,29,30,27,28,29,30,26,27,28,29,26,27,28, 1 24,25,26,23,24,25,22,23,24,25,20,21,22,23,24,25,19,20, 2 25,21,26,27,28,29,30,27,28,29,30,26,27,28,29,30,26,27, 3 28,23,24,25,26,27,28,29,23,26,27,28,29,20,21,22,23,27, 4 28,29,5926*1/ C END SUBROUTINE INTDLT C JELESNIANSKI AUGUST 1980 TDL IBM 360/195 C C PURPOSE C THIS SUBROUTINE (INTDLT) ASCERTAINS TIME STEP 'DELT' FOR C COMPUTATIONS, INITIALIZES STORM AT SUBSCRIPT 'IBGNT' AND C SETS 'ITEND' FOR ENDING COMPUTATIONS AFTER NEAREST STORM C APPROACH. C C DATA SET USE C NONE C C VARIABLES C DELT = TIME INTERVAL FOR EXPLICIT FINITE DIFFERENCING C DLTB = THE SMALLEST REDUCTION OF DELT IN THE RUN C RLNGTH = DISTANCE OF STORM FROM NEAREST APPROACH C JHR = SUBSCRIPT FOR NEAREST APPROACH TO BASIN CENTER C P(100) = HOURLY PRESSURE DROPS C X(100) Y(100) = HOURLY COMPONENTS FOR STORM POSITIONS C R(100) = HOURLY STORM SIZES C XMOUTH = X-DISTANCE FROM POLE TO BASIN CENTER C YMOUTH = Y-DISTANCE (USUALLY EQUAL TO 0) C IBGNT = TABULATED SUBSCRIPT FOR STORM INITIALIZATION C ITEND = TABULATED SUBSCRIPT FOR STORM TERMINATION C IBGS = NUMBER OF HOURS FOR COMPUTATIONS BEFORE C NEAREST STORM APPROACH C ITMS = NUMBER OF REMAINING HOURS, AFTER NEAREST APPROACH, C DLTIN(3) = TIME STEPS FOR 3 CATAGORIES OF STORMS FROM DATA C C GENERAL COMMENTS C THIS SUBROUTINE LIES IN OVERLAY 'STRMPR'. IT IS CALLED IN C BY SUBROUTINE 'STRMPR'. C COMMON /EGTH/ DELS,DELT,G,COR COMMON /DLTDTA/ DLTIN(3),DLTB COMMON /BSN/ PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C COMMON /SPLN/ ALT(15),ALN(15),AX(15),AY(15),RL(15),ANGD(15), 1 PT(15),RT(15),XLAT(999),YLONG(999),RLNGTH COMMON /STRMPS/ X(999),Y(999),P(999),R(999),DIR(999),SP(999) C STIME interferes with C code, so switched to STIME2 COMMON /STIME2/ ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND COMMON /ORIGN/ XMOUTH,YMOUTH COMMON /INZSTM/ IBGS,ITMS COMMON /OPT/ NOPT COMMON /FLES/ FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1 CHARACTER*256 FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1 C C SET DLTB FOR REDUCTION OF DELT. DLTB=DLTIN(2) C C RESET A SMALLER DLTB IF STORM IS STRONG OR STORM IS CLOSE. !------------------------------------------------------------------ !rewrite codes to remove GO TO statments Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------------ ! IF(RLNGTH.GT.90.) GO TO 100 ! MAXP=AMAX1(P(JHR),P(JHR-3)) ! IF (MAXP.GT.60) DLTB=DLTIN(3) C ! 100 CONTINUE C SET INITIAL DELT. ! DELT=DLTIN(1) C ! IF (DELT.GT.500..OR.DELT.LT.5.) GO TO 130 ! GO TO 150 ! 130 CONTINUE CC WRITE(*,140) DELT,(DLTIN(I),I=1,3) ! 140 FORMAT(' TIME STEP=',F10.3,' IS IMPROPERLY CHOSEN FROM BASIN', ! 1' DATA ',3F15.6) c STOP ! 150 CONTINUE C ! ! IF (RLNGTH.LE.90.) THEN MAXP=AMAX1(P(JHR),P(JHR-3)) IF (MAXP.GT.60) THEN DLTB=DLTIN(3) END IF END IF DELT=DLTIN(1) IF (DELT.GT.500..OR.DELT.LT.5.) THEN ! WRITE(*,140) DELT,(DLTIN(I),I=1,3) 140 FORMAT(' TIME STEP=',F10.3,' IS IMPROPERLY CHOSEN FROM BASIN', 1' DATA ',3F15.6) ! STOP END IF ! ! ! ! IBGNT=MAX0(JHR-IBGS,1) RETURN END SUBROUTINE INTVAL C JELESNIANSKI AUGUST 1980 TDL IBM 360/195 C PURPOSE C THIS SUBROUTINE (INTVAL) SETS INITIAL STORM PARAMETERS, C SOME FIXED CONSTANTS FOR INITIALIZATION, AND OUTPUT TIMES. C C DATA SET USE C NONE C C VARIABLES C DELT = TIME INTERVAL FOR EXPLICIT FINITE DIFFERENCING C NDLTHR NDLTH2 = NUMBER OF TIME STEPS IN ONE AND 1/2 HOUR C IPRTSV = AT THIS TIME PRINT HISTORICAL SURGES (10 POINTS) C IPRTAD = ADVANCE 'IPRTSV' BY 'IPRTAD' C IBGNT = SUBSCRIPT FOR STORM INITIALIZATION C X(100) Y(100) = HOURLY COMPONENTS FOR STORM POSITIONS C SP(100) = HOURLY FORWARD STORM SPEEDS C DIR(100) = HOURLY STORM AZIMUTHS C R(100) = HOURLY STORM SIZES C IBGS = NO OF HRS FOR COMPUTATIONS BEFORE NEAREST APPROACH C ITMS = NO OF HRS FOR COMPUTATIONS AFTER NEAREST APPROACH C MHALT = number of time steps for total compuations C ITEND = TABULATED SUBSCRIPT FOR STORM TERMINATION C IC19 = GROWTH TIME TO MATURE STORM C ALTO = LATITUDE OF BASIN CENTER C MAPIN(10) = PRINT OUT TIME FOR ENVELOPE SNAPSHOT C IOPERL = 2 IF IN OPERATIONAL MODE C ISTMAP = AT THIS TIME, PRINT A SNAPSHOT OF SURGES C NUMBER = NEAREST HOUR TO NEAREST APPROACH C C GENERAL COMMENTS C USE ONLY TIME STEPS DIVISIBLE BY 600 S(10 MIN),OR 900 S( C 15 MIN). C THIS SUBROUTINE LIES IN OVERLAY STRMPR. IT IS CALLED IN BY C SUBROUTINE STRMPR. C COMMON /EGTH/ DELS,DELT,G,COR COMMON /DELTH/ NDLTHR,NDLTH2,DLTHR,DLTH2 COMMON /FRTH/ IC12,IC19,BT COMMON /PRTSV/ IPRTSV,IPRTAD C STIME interferes with C code, so switched to STIME2 COMMON /STIME2/ ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND COMMON /STRMPS/ X(999),Y(999),P(999),R(999),DIR(999),SP(999) COMMON /INZSTM/ IBGS,ITMS COMMON /FFTH/ ITIME,MHALT COMMON /BSN/ PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C COMMON /DELH11/ MAPIN(10),NUMBER COMMON /SWTCH/ IOPERL(5) COMMON /MAPTME/ ISTMAP,IADMAP C C NUMBER OF TIME STEPS IN AN HOUR/HALF-HOUR NDLTHR=3600./(DELT-.001) NDLTH2=NDLTHR/2 DLTHR =NDLTHR DLTH2 =NDLTH2 C C NUMBER OF TIME STEPS BETWEEN TIME-HISTORY PRINT-OUT. IPRTAD=NDLTHR/4 IF(MOD(NDLTHR,6).EQ.0) IPRTAD=NDLTHR/6 IPRTSV=IPRTAD C C C C TOTAL RUN TIME T=ITEND-IBGNT C write (*,*) "Mhalt", MHALT, (T*3600.+1.)/DELT MHALT=(T*3600.+1.)/DELT CC write (*,*)' total hours and time steps',t,mhalt C C PRE-SET STORM GROWTH TO 10 HOURS IC19=10*NDLTHR C C COMPUTE STORM GROWTH TIME TO MATURITY ! IF (SP(IBGNT).LT.1.E-5) GO TO 100 C C NO OF TIME STEPS FOR STORM TO MOVE 1 DIAMETER IN LENGTH ! IGRWTH=(2.*R(IBGNT)/SP(IBGNT))*DLTHR ! GO TO 110 ! 100 IGRWTH=4*NDLTHR+4 ! 110 CONTINUE ! IC19=MAX0(IC19,MIN0(4*NDLTHR+4,IGRWTH)) C ! DO 120 I=1,5 ! J=IBGNT+I-1 ! IF(ABS(X(J)-1.).LT.(.25*R(IBGNT))) IC19=4*NDLTHR ! 120 CONTINUE C CC write (*,*) ' spin-up hours = ',(ic19+1)/ndlthr C SET MAPIN(10) AS SELECTED HOURS, FOR SNAPSHOT OUTPUTS C IF(IOPERL(1).NE.2) GO TO 125 C MAPIN(1)=0 C DO 122 I=2,10 C 122 MAPIN(I)=200 C 125 CONTINUE ! DO 130 I=1,10 ! 130 MAPIN(I)=JHR+MAPIN(I) C ! DO 140 I=1,10 ! IF(MAPIN(I).GT.IBGNT+1) GO TO 150 ! 140 CONTINUE C ! 150 NUMBER=I C SET PRINT OUT TIME OF STORM AT NEAREST APPROACH C INITIALIZE OUTPUT TIME FOR SNAPSHOT ENVELOPE OF SURGES ! ISTMAP=MAPIN(I) RETURN END SUBROUTINE CNSTNT C JELESNIANSKI SEPTEMBER 1980 TDL IBM 360/195 C C PURPOSE C THIS SUBROUTINE (CNSTNT) FIXES SOME PROGRAM CONSTANTS. C C DATA SET USE C NONE C C VARIABLES C X1(50) = STORED ARRAY OF CONSTANTS C DELT = TIME INCRMNT FOR FINITE DIFFRNCNG, SET IN 'INTDLT' C COR = CORIOLIS (1/SEC), (SET IN 'INTVAL') C IMXB JMXB = NO. OF GRID PTS ALONG A RAY/CIRCLE, IN 'CRDRD1' C IMXB1 JMXB1 = IMXB LESS 1, JMXB LESS 1 C IMXB2 JMXB2 = IMXB LESS 2, JMXB LESS 2 C DEGREE = ANGLE BETWEEN TWO RAYS, READ IN 'CRDRRD1' C DELA = VALUE OF 'DEGREE' IN RADIANS C RAD = RADIAN MEASURE C X1B(10) = STORED ARRAY OF CONSTANTS C G = GRAVITY = 32.2 FT/SEC-2 C AMAN->AGB13 = FIXED CONSTANT FOR SUBROUTINE 'CHANNL' C C GENERAL COMMENTS C THIS SUBROUTINE IS CALLED BY SUBROUTINES INITLZ AND STMVAL. C C COMMON /SLAT/ FSOUTH COMMON /FRST/ X1(50),X12(50) COMMON /EGTH/ DELS,DELT,G,COR COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /POLAR/ DEGREE,RMOUTH,AZMTH,DELA COMMON /GAMAF/ GAMA,GAMA1,GAMFP(4) COMMON /DUMB4/ X1B(10) COMMON /DUMB4I/ CORL(90) COMMON /CHNLNO/ AMAN,BMAN,GB,AGB,AGB13 X1B(1)=G*DELT/DELA X1B(2)=DELT/DELA X1B(4)=DELT*COR DO 10 L=1,90 CORL(L)=2.*(7.292116E-5)*SIN((L-1)*1.74532925E-2)/COR 10 CONTINUE X1B(5)=.5*X1B(2) X1B(6)=.5*X1B(1) C WEIR COEFFICIENT X1B(7)=.6*DELT*SQRT(2.*G/3.)/3. C SET SMOOTHING CONSTSNTS FOR POLAR COORDINATES EXPDLA=EXP(DELA) COSHLA=.5*(EXPDLA+1./EXPDLA) GAMA =EXPDLA/COSHLA-1. GAMA1=-GAMA GAMFP(1)=0. GAMFP(3)=0. GAMFP(2)=GAMA GAMFP(4)=GAMA1 C COMPUTING FIXED CONSTANTS FOR SUBROUTINE 'CHANNL' CMAN =1.5/.03 AMAN =3.*(CMAN**2)/13. BMAN =3.*(CMAN**2)/(4.*G) GB =G*BMAN AGB =AMAN-G*BMAN AGB13=13.*AGB/3. RETURN END SUBROUTINE BTMSTR(ZLATO) C JELESNIANSKI SEPTEMBER 1980 TDL IBM 360/195 C C PURPOSE C THIS SUBROUTINE (BTMSTR) COMPUTES BOTTOM STRESS COEFFCNTS C AT 1 FOOT INTERVALS, FROM 1 TO 300 FEET. C THE EQUATIONS ARE IN REFERENCE: JELESNIANSKI, C "NUMERICAL COMPUTATIONS OF STORM SURGES WITH BOTTOM C STRESS", MONTHLY WEATHER REVIEW, VOL 95, NO. 11, C NOV 1967, PP 740-756. C C DATA SET USE C NONE C C VARIABLES C ZLATO = LATITUDE IN DEGREES C COR = CORIOLIS PARAMTER C E = EKMAN PARAMETER, DEPTH*SQRT(COR/2*C25) C C25 = EDDY VISCOSITY COEFFICIENT, .25 FT**2/SEC C C7 = SLIP COEFFICIENT, .006 FT/SEC C AR(300) AI(300) = BOTTOM FRICTION COEFFS FOR CORIOLIS TERMS C BR(300) BI(300) = BOTTOM FRICTION COEFFS FOR SFC GRAD TERMS C CR(300) CI(300) = BOTTOM FRICTION COEFFS FOR SURFACE STRESS TERMS C C GENERAL COMMENTS C THIS MEMBER 'BTMSTR' RESIDES IN OVERLAY 'INITLZ'. IT IS C CALLED IN BY SUBROUTINE 'INITLZ'. C PROGRAM REWRITTEN (1980), REPLACING SINH BY SIN AND COSH C BY COS. C C parameter (ndp=600) COMMON /SCND/ AR(ndp),AI(ndp),BR(ndp),BI(ndp),CR(ndp),CI(ndp) COMPLEX*16 SIGMAI,CO1,CO2,CO3,CO4,CO5,CO7,CO17 COMPLEX*16 CO16,CO10,CP REAL*8 E DATA C25,C7/.25,.006/ C CORIOLIS PARAMETER C FSOUTH=1. C IF (ZLATO.LT.0.) FSOUTH=-1. COR=2.*(7.292116E-5)*SIN(ABS(ZLATO)*1.74532925E-2) CC PRINT 5 5 FORMAT(//' CONSTANTS FOR SLIP CONDITION AT BOTTOM'// 1' E AR AI BR BI CR CI 2 DEPTH') NN=1 NND=1 C DO 100 N=1,ndp A=N E=A*SQRT(COR/(2.*C25)) COR1=COR*A/C7 C SIGMAI=DCMPLX(-E,E) CO1=ZCOS(SIGMAI) CO2=SIGMAI/ZSIN(SIGMAI) CO1=CO1*CO2 C CO1=SIGMA/TANH(SIGMA) , CO2=SIGMA/SINH(SIGMA) C CO1=SIGMAI/TAN(SIGMAI), CO2=SIGMAI/SIN(SIGMAI), SIGMAI=I*SIGMA CO3=1./(CO1+DCMPLX(-1.,COR1)) CO4=CO3*CO3 CO5=CO2*CO2+CO1 CP=1.-CO2 CO7=(.5*CO5-1.)*CO4 CO17=1./(1.+CO7) CO16=(1.+CO3)*CO17 AR(N)=DREAL(CO16) AI(N)=DIMAG(CO16) BR(N)=DREAL(CO17) BI(N)=DIMAG(CO17) CO10=(1.+CP*CO3)*CO17 CR(N)=DREAL(CO10) CI(N)=DIMAG(CO10) !-------------------------------------------------- ! rewrite codes to remove GO TO statments ! Huiqing.Liu /MDL Feb. 2015 !-------------------------------------------------- IF (N.NE.NN) THEN CYCLE END IF IF (N.GT.10) NND=15 NN=NN+NND C PRINT 120,E,AR(N),AI(N),BR(N),BI(N),CR(N),CI(N),N 120 FORMAT(F6.2,6F10.5,I10) 100 CONTINUE RETURN END SUBROUTINE CNSDLT C J CHEN AUG 1982 IBM 360/195 C C PURPOSE C THIS SUBROUTINE (CNSDLT) FIXES COEFFICIENTS DEPENDING ON C TIME STEP FOR VARIATIONS IN THE COURSE OF COMPUTATION.. C C DATA SET USE C NONE C C VARIABLES C X1(50) = STORED ARRAY OF CONSTANTS C DELT = TIME INCRMNT FOR FINITE DIFFRNCNG, SET IN 'INTDLT' C COR = CORIOLIS (1/SEC), (SET IN 'INTVAL') C DEGREE = ANGLE BETWEEN TWO RAYS, READ IN 'CRDRD1' C DELA = VALUE OF 'DEGREE' IN RADIANS C RAD = RADIAN MEASURE C X1B(10) = STORED ARRAY OF CONSTANTS C G = GRAVITY = 32.2 FT/SEC-2 C C GENERAL COMMENTS C THIS SUBROUTINE IS CALLED BY SUBROUTINES INITLZ AND STMVAL. C C COMMON /SLAT/ FSOUTH COMMON /EGTH/ DELS,DELT,G,COR COMMON /FRST/ X1(50),X12(50) COMMON /POLAR/ DEGREE,RMOUTH,AZMTH,DELA COMMON /DUMB4/ X1B(10) C X1B(1)=G*DELT/DELA X1B(2)=DELT/DELA c X1B(4)=DELT*COR C C FSOHTH=-1. FOR SOUTHERN HEMISPHERE, NEGATIVE CORIOLIS PARAMETER C = 1. FOR NORTHERN HEMISPHERE, c C X1B(4)=DELT*COR*FSOUTH C X1B(5)=.5*X1B(2) X1B(6)=.5*X1B(1) C WEIR COEFFICIENT X1B(7)=.6*DELT*SQRT(2.*G/3.)/3. RETURN END SUBROUTINE CRDRD1(IVER,IPRJ) C JYE CHEN MAY 1990 TDL IBM 360/195 C PURPOSE C THIS SUBROUTINE (CRDRD1) READS IN BASIN PROJECTION DATA, C AS INPUT DATA TO THE SLOSH PROGRAM. EACH BASIN DATA IS ON C A SEPARATE FILE, AND CALLED IN BY NAME. C C DATA SET USE C FT09F001 C C VARIABLES C STA = ALPHAMERIC NAME OF A BASIN (10 LETTERS) C EBSN = '$' INDICATES TYPE I ELLIPTIC COORDINATES, C '+' INDICATES TYPE II ELLIPTIC COORDINATES, C OTHERWISE, POLAR COORDINATES C DOLLAR = '$' USED FOR ISLAND, OTHERWISE FOR REGULAR BASIN C AZMTH = SLANT OF X-AXIS FROM NORTH C ALTO ALNO = LATITUDE/LONGITUDE OF BASIN CENTER OR ORIGIN C DEGREE = ANGLE BETWEEN TWO RAYS OF A POLAR GRID C DELA = RADIAN MEASURE OF DEGREE C IRB IRE = TOTAL GRID POINTS IN -+P DIRECTION FROM TANGENT PT C NDB NDE = TOTAL GRID POINTS IN -+Q DIRECTION FROM TANGENT PT C XMOUTH = X-COORDINATE FROM TANGENT POINT TO (X,Y)=(P,Q)=(0, C YMOUTH = Y-COORDINATE FROM TANGENT POINT TO (X,Y)=(P,Q)=(0, C PHI = SLANT FROM X-AXIS TO N/S-AXIS C IMXB JMXB = TOTAL GRID SPACINGS IN P/Q DIRECTIONS C IMXB1 JMXB1 = IMXB LESS 1, JMXB LESS 1 C C GENERAL COMMENTS C common /opts/ nofld,nof1d character*1 nofld,nof1d COMMON /SLAT/ FSOUTH CHARACTER*1 ISOUTH COMMON /XOKE/ XOKE CHARACTER*1 XOKE COMMON /GPRT/ STA COMMON /GPRT1/ DOLLAR,EBSN COMMON /GPRT2/ EBSN1,EBSN2 COMMON /POLAR/ DEGREE,RMOUTH,AZMTH,DELA COMMON /INDX/ IRB,IRE,NDB,NDE COMMON /BSN/ PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C COMMON /DLTDTA/ DLTIN(3),DLTB COMMON /ORIGN/ XMOUTH,YMOUTH COMMON /ELLIP/ AAXIS,BAXIS,ABQAB COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /EGTH/ DELS,DELT,G,COR COMMON /HTERAIN/ HTER CHARACTER*1 ITERRAIN COMMON /SMTH/ ISMTH CHARACTER*1 ISMOOTH C COMMON /FLES/ FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1 CHARACTER*256 FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1 CHARACTER*16 STA CHARACTER*2 DOLLAR CHARACTER*1 EBSN,EBSN1,EBSN2 CHARACTER*12 ELLIPSOID C ADDED FOR ELLIPSOID CHOICE BETWEEN GRS80 AND CLARKE D.Y 8/2014 INTEGER IVER,IPRJ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DTA OPTIONS FOR DTA VERSION AND ELLIPSOID TYPE IF EXIST C IVER = 199201(ORIGINAL DTA FORMAT) / 201408(NEW DTA FORMAT) C IPRJ = 0(CLARKE ELLIPSOID) / 1(GRS80 ELLIPSOID) C D.Y 8/2014 C COMMON /DTAOPT/ IVER IPRJ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ---------------------------------------------------------------------- C WRITE(*,*)' ENTER BASIN PROJECTION FILE ' C READ (*,'(A)') FLE9 OPEN(9,FILE=FLE9) C--------------------------------------------------------------------------- C AT COLUMN 12, '$' INDICATES TYPE I ELLIPTIC COORDINATES, C '+' INDICATES TYPE II ELLIPTIC COORDINATES, C OTHERWISE, POLAR COORDINATES C AT COLUMN 13, '$ ' INDICATES FOR CLOSED ISLAND (PERIODIC B.C.) C '2$' INDICATES FOR MSY BASIN C AT COLUMN 15, '&' INDICATES NO CORNER SMOOTHING C '+' INDICATES GLOBAL SMOOTHING EVERY DELT C AT COLUMN 16, '+' INDICATES OUTPUT OPTION (J-I). DEFAULT (I-J) c 17, '&' SOUTHERN HEMISPHERE C '$' XOKE = 'X' c 18, '+' noflooding and overtopping of barriers c 19, '+' no 1d flow allowed C 20, '+' allow between 35 and 56 feet to flood. C 21, '+' x-filter smoothing method. (HCRT=0.5) C '-' x filter smoothing method. (HCRT=0.1) C '@' x filter smoothing method. (HCRT=0.2) C '#' x filter smoothing method. (HCRT=0.3) C '$' x filter smoothing method. (HCRT=0.4) C '*' xfilter exclusion of cells by Epsilon. (HCRT=0.1) C '*' redefined as (HCRT=Epsilon) C '%' xfilter exclusion of cells by Epsilon. (HCRT=0.2) C '^' xfilter exclusion of cells by Epsilon. (HCRT=0.3) C '&' xfilter exclusion of cells by Epsilon. (HCRT=0.4) C '=' xfilter exclusion of cells by Epsilon. (HCRT=0.5) C--------------------------------------------------------------------------- READ (9,'(A10,1X,A1,A2,8A1)') STA,EBSN,DOLLAR,EBSN1,EBSN2,ISOUTH, 1 nofld,nof1d,iterrain,ismooth C READ IN DTA VERSION AND ELLIPSOID TYPE FROM THE FIRST TWO LINES IF EXIST C D.Y 8/2014 IF (STA .EQ. 'Ver2014-08') THEN IVER = 201408 C WRITE(*,*) 'READING NEW DTA FORMAT.' READ (9,'(A11)') ELLIPSOID IF (TRIM(ELLIPSOID) .EQ. 'Proj=GRS80') THEN IPRJ = 1 C WRITE(*,*) 'ELLIPSOID = GRS80.' ELSEIF (TRIM(ELLIPSOID) .EQ. 'Proj=Clarke') THEN WRITE(*,*) 'ELLIPSOID = CLARKE' ELSE WRITE(*,*) 'UNRECOGNIZED ELLIPSOID INFO IN FILE: ',ELLIPSOID STOP ENDIF READ (9,'(A10,1X,A1,A2,8A1)')STA,EBSN,DOLLAR,EBSN1,EBSN2,ISOUTH, 1 nofld,nof1d,iterrain,ismooth ENDIF ! Huiqing.Liu /MDL basin file changed from 9 to 19 IF (ISOUTH.EQ.'$') XOKE='X' FSOUTH=1. IF (ISOUTH.EQ.'&') FSOUTH=-1. HTER=35. IF (ITERRAIN.EQ.'+') HTER=56. ISMTH=1 IF (ISMOOTH.EQ.'+') ISMTH=2 IF (ISMOOTH.EQ.'-') ISMTH=3 IF (ISMOOTH.EQ.'@') ISMTH=4 IF (ISMOOTH.EQ.'#') ISMTH=5 IF (ISMOOTH.EQ.'$') ISMTH=6 IF (ISMOOTH.EQ.'*') ISMTH=7 IF (ISMOOTH.EQ.'%') ISMTH=8 IF (ISMOOTH.EQ.'^') ISMTH=9 IF (ISMOOTH.EQ.'&') ISMTH=10 IF (ISMOOTH.EQ.'=') ISMTH=11 C READ (9,580) AZMTH C TLAT,TLONG NOT BEING USED. USE XMOUTH,YMOUTH, INSTEAD. IF (EBSN.EQ.'$'.OR.EBSN.EQ.'+') READ (9,580) TLAT,TLONG C LAT/LONG OF BASIN CENTER AND TANGENT POINT ON EARTH READ (9,580) ALTO,ALNO C C ------ FOR POLAR GRID SYSTEM, READ IN DEGREE OF SPREAD OF RAYS --- IF (EBSN.NE.'$'.AND.EBSN.NE.'+') THEN READ (9,580) DEGREE DELA=DEGREE*1.74532925E-2 ENDIF C ------ READ IN DISTANCE FROM BASIN CENTER TO 1ST HYPERBOLA, C ANGLE OF 1ST ASYMPTOTE FOR ELLIPTIC GRIDS I ------ IF (EBSN.EQ.'$') READ (9,580) DSTNT,ASMPT C ------ READ IN AAXIS,BAXIS FOR ELLIPTIC GRIDS II --- IF (EBSN.EQ.'+') READ (9,580) AAXIS,BAXIS,GSIZE C READ (9,'(5I5)') IRB,IRE,NDB,NDE,NDB1 C DEFINE DELTA THETA FOR ELLIPTIC GRIDS I AND II. IF (EBSN.EQ.'$') THEN IF (NDB.NE.0) THEN NDBZ=NDB ELSE NDBZ=NDE ENDIF C SPECIAL USE OF ELLIPTIC BASIN I (CARTESIAN APPROX.) C -- Arthur -- doesn't occur (3/30/2011) in any basin? --- IF (EBSN1.EQ.'$') NDBZ=NDB1 DEGREE=(90.-ASMPT)/NDBZ DELA=DEGREE*1.74532925E-2 ASMPT=ASMPT*1.74532925E-2 AAXIS=0. BAXIS=DSTNT/COS(ASMPT) ENDIF IF (EBSN.EQ.'+') THEN C DEGREE=360./(NDB+NDE) C DELA=DEGREE*1.74532925E-2 DELA=2.*GSIZE/(AAXIS+BAXIS) DEGREE=DELA/1.74532925E-2 XMOUTH=AAXIS YMOUTH=0 ENDIF C PHI =270.+AZMTH PHI1=1.74532925E-2*PHI ALT1=1.74532925E-2*ALTO C C FORMULATE CORIOLIS, SEC-1 COR=2.*(7.292116E-5)*SIN(ALT1) ALT1C=1.74532925E-2*(90.-ALTO) ALN1=1.74532925E-2*ALNO C C READ IN DELTS FOR 3 CATAGORIES OF STORM READ (9,100) (DLTIN(I),I=1,3) 100 FORMAT(3F10.6) C IF (EBSN.NE.'+') READ (9,580) XMOUTH,YMOUTH C DEFINE AAXIS=BAXIS= DISATANCE FROM POLE TO TANGENT POINT C FOR POLAR GRIDS IF (EBSN.NE.'$'.AND.EBSN.NE.'+') THEN AAXIS=XMOUTH BAXIS=AAXIS ENDIF C C DEFINE CONSTANTS FOR GENERAL CODES. RMOUTH=.5*(AAXIS+BAXIS) ABQAB=(AAXIS-BAXIS)/(AAXIS+BAXIS) C IMXB =IRB+IRE+1 JMXB =NDB+NDE+1 C CC WRITE (*,*) ' ---- GRID TRANSFORMATION DATA ----' CC WRITE (*,'(A)') ' AAXIS BAXIS DELA(DEG) IMXB JMXB (A-B CC 1)/(A+B)' CC WRITE (*,111) AAXIS,BAXIS,DELA/1.74532925E-2,IMXB,JMXB,ABQAB 111 FORMAT(2F10.4,F10.6,2I6,F10.5) C IMXB1 =IMXB-1 JMXB1 =JMXB-1 IMXB2=IMXB1-1 JMXB2=JMXB1-1 C 580 FORMAT (3F10.6) RETURN END SUBROUTINE TTIMER C JELESNIANSKI/CHEN JULY 1989 TDL IBM 360/195 C C PURPOSE C IT KEEPS TRACK OF TIME AND ACCUMULATES TIME STEPS. C IT UPDATES THE STORM POSITION AT EACH TIME STEP ON A C CARTESIAN GRID. IT COMPUTES A GROWTH FACTOR TO UPDATE C A STORM TO MATURITY. C C DATA SET USE C NONE C C VARIABLES C IC19 = STORM GROWTH TIME TO MATURITY, SET IN 'INTVAL' C C Y ( C Y *(X,Y) C Y ) C Y C PXXXXX0XXXXXXXXXXXXXXXXXXXXXXXX C 0 ORIGIN=TANGENT POINT C L Y ON EARTH C E Y C Y C C ITIME = INCREASES BY ONE FOR EACH TIME STEP C BT = STORM GROWTH COEFFICIENT, 0= 10 feet above datum so it is dry. C EOKE (v2) NorthWest Channel wasn't completely wet. ELSE HB(I,J)=-ZB(I,J) ENDIF 10 CONTINUE C Start added Arthur for Ok3... DO L=1,NODRY I=IDRY(L) J=JDRY(L) HB(I,J)=-ZB(I,J) ENDDO C Finished added Arthur for Ok3... C RESET CANNEL WATER LEVEL AS LAKE LEVEL C SUPPRESS WIND FOR ALL CANALS WITH 'TREE' OPTION IN 1D FLOW. DO 12 L=1,NSQRWC IF (KTREE(L).EQ.'T') THEN I=ISQR(L) J=JSQR(L) HB(I,J)=AMAX1(DTMLAK,-ZB(I,J)) K=ISIDE(L) II=I+IHH(K) JJ=J+JHH(K) HB(II,JJ)=AMAX1(DTMLAK,-ZB(II,JJ)) ENDIF 12 CONTINUE C GOTO 1112 C C ------------------------ FOR MSY ONLY ----------------- C 1211 CONTINUE C INITIALIZE WATER HEIGHTS ACCORDING TO LAKE DATUM DO 190 J=1,JMXB1 IFN=LDTMG(J) DO 190 I=1,IMXB1 IF (I.LE.IFN) THEN HB(I,J)=AMAX1(-ZB(I,J),AMAX1(DTMLAK,-ZB(I,J))) ELSE IF (ITREE(I,J).EQ.'2'.OR.ITREE(I,J+1).EQ.'2'.OR. 1 ITREE(I+1,J).EQ.'2'.OR.ITREE(I+1,J+1).EQ.'2') THEN ! Test Initial water level 0.5ft by H.Liu 8/2018 ! SEADTM = 0.5 HB(I,J)=.25*(PP(I,J)+PP(I+1,J)+PP(I,J+1)+PP(I+1,J+1))+SEADTM ! HB(I,J) = SEADTM HB(I,J)=AMAX1(HB(I,J),-ZB(I,J)) ELSE ! DTMLAK = 0.5 ZZ=.25*(PP(I,J)+PP(I+1,J)+PP(I,J+1)+PP(I+1,J+1)) +DTMLAK ! ZZ = DTMLAK HB(I,J)=AMAX1(-ZB(I,J),AMAX1(ZZ,-ZB(I,J))) ENDIF 190 CONTINUE C C DRY OUT LAND SQUARES BELOW SEA LEVEL. DO 195 L=1,NODRY I=IDRY(L) J=JDRY(L) HB(I,J)=-ZB(I,J) 195 CONTINUE c C 1112 CONTINUE C INITIALIZE U/V/HMX C DO 320 I=1,IMXB1 320 HB(I,JMXB)=HB(I,1) DO 200 J=1,JMXB DO 200 I=1,IMXB UB(I,J)=0. 200 VB(I,J)=0. DO 210 J=1,JMXB1 DO 210 I=1,IMXB1 HMX(I,J)=HB(I,J) 210 CONTINUE C C C INITLIZE FLWSQR TO BE ZERO. C DO 220 L=1,NSQRWC 220 FLWSQR(L)=0. C REDUCE PRINTOUT IF IN OPERATIONAL MODE. ! IF(IOPERL(1).EQ.2) GO TO 300 CC WRITE(*,140) C PRINT OUT INITIAL HEIGHT VALUES CC WRITE(*,150) ! NGRP=(JMXB1-1)/25+1 C 1130 CONTINUE C 140 FORMAT(1H1) 150 FORMAT(/' INITIAL HEIGHT VALUES') 160 FORMAT(1H ,I2,26F5.0) 170 FORMAT(1H ,I2,26F5.1) 180 FORMAT(/I6,26I5) 300 CONTINUE RETURN END SUBROUTINE RDLTLG INCLUDE 'parm.for' C COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /FLES/ FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1 CHARACTER*256 FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1 CHARACTER*256 FILNAM CHARACTER*80 DUM REAL *4 XX C CC WRITE (*,*)' READ IN LAT/LONG OF THE GRIDS' C READ (*,'(A)') FILNAM FILNAM = FLE99 CALL TOPEN(25,FILNAM,0,2) CALL TREAD(25,IDUM,1) CALL TREAD(25,JDUM,1) CALL TREADC(25,DUM,80) CALL TREADC(25,DUM,80) IF (IDUM.NE.IMXB.OR.JDUM.NE.JMXB) THEN WRITE (*,*) "LLX ERROR:", IDUM, JDUM, IMXB, JMXB ENDIF DO L=1,JMXB DO K=1,IMXB CALL TREADF(25,XX,1) YLT(K,L) = XX END DO END DO DO L=1,JMXB DO K=1,IMXB CALL TREADF(25,XX,1) YLT(K,L) = XX END DO END DO CALL TCLOSE(25) C OPEN (25,FILE=FILNAM,ACCESS='TRANSPARENT') C READ (25,end=26) IDUM,JDUM C READ (25) DUM,DUM C READ (25) ((YLT(K,L),K=1,IMXB),L=1,JMXB) C READ (25) ((YLG(K,L),K=1,IMXB),L=1,JMXB) C CLOSE (25) 26 continue C RETURN END SUBROUTINE CRDRD2 INCLUDE 'parm.for' C INCLUDE 'parmmsy.for' C COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /BSN/ PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C COMMON /SPLN/ ALT(15),ALN(15),AX(15),AY(15),RL(15),ANGD(15), 1 PT(15),RT(15),XLAT(999),YLONG(999),RLNGTH COMMON /DUMBSV/ SAVEH(10),IPN(10),JPN(10),STATS(10) COMMON /PRTOVL/ NSHEET,JPGS(10),JPGF(10),IPGS(10),IPGF(10) COMMON /FLWCPT/ NSQRS,NSQRW,NSQRWC,NPSS,NCUT COMMON /DLTDTA/ DLTIN(3),DLTB COMMON /CHNLD/ CHLNH(5),CHWTH(5),CHDPH(5) COMMON /CHNL/ IPT0(2,5),JPT0(2,5),IPTL(2,5),JPTL(2,5),ENTEXT(5) COMMON /IPRX/ IX60 COMMON /GPRT1/ DOLLAR,EBSN CHARACTER*1 EBSN !------------------------------------------------------- ! Added Basin Name Variable by Huiqing.Liu/MDL June/2015 ! for Modifying Cd wind drag coefficient parameter only ! in New Bering Sea extra-tropical basin (eno3) !------------------------------------------------------- COMMON /GPRT/ STA CHARACTER*16 STA CHARACTER*255 FN_DEP_NEP logical :: ThereCSV integer :: SizeCSV C CHARACTER*16 STATS DIMENSION IANG(L_) CHARACTER*16 AAA DIMENSION INW(7),JNW(7),HW(7) DIMENSION IBNK(NBK_),JBNK(NBK_),ISBNK(NBK_),BKI(NBK_), 1 BKE(NBK_),WCHN(NBK_) CHARACTER*2 AC,BC,IX60,XA1,YA1,ZA1,DOLLAR,SIG_ DATA XA1/'X'/,YA1/'Y'/,ZA1/'Z'/,SIG_/'$'/ C C READ BOUNDARY TYPES ALONG LEFT,RIGHT,TOP AND BOTTOM BOUNDARIES C ALONG MOMENTUM POINTS. C .****.****.****.****.****.****.****.****. TOP BNDRY C * * C * * C . . . . . . . . . C LEFT BNDY * * RIGHT BNDY C * + + + + * C * * C .****.****.****.****.****.****.****.****. BOTTOM BNDRY C 2 8 8 8 9 9 9 1 1 C DEEP WATER-) (-INTERMEDIATE----) (-SHALLOW----) (--LAND---- C IF (IMXB.GT.M_.OR.JMXB.GT.N_) THEN WRITE (*,*) ' BASIN DIMENSIONS ',IMXB,JMXB WRITE (*,*) ' HAS EXCEEDS THE COMPILE DIMENSIONS',M_,N_ STOP endif !------------------------------------------------------------------- ! Huiqing.Liu /MDL read real depth from another file to overcome ! the 4 digit limit in dta files !------------------------------------------------------------------- if (STA(1:7) == 'AK GULF'.or.STA(1:10) == 'NEP AK GUL') then FN_DEP_NEP='grid_depth_nep.txt' inquire(file=trim(FN_DEP_NEP),exist=ThereCSV) inquire(file=trim(FN_DEP_NEP),size=SizeCSV) if (ThereCSV.and.SizeCSV.gt.0) then ! Added by Huiqing.Liu /MDL in July 2020 to make sure the parm is not empty open(101,file=trim(FN_DEP_NEP)) do i=1,IMXB-1 do j=1,JMXB-1 read(101,'(2I5,F11.2,I2)')i_dex,j_dex,zb_tmp,i_boun(i,j) enddo enddo close(101) write(*,*) 'Depth_nep is ok' else i_boun=1 write(*,*) 'Depth_nep is not ok' endif endif !-------------------------------------------------------------------- C INITIALIZE ZBM(I,J) FOR ENTIRE GRIDS DO 120 J=1,JMXB DO 120 I=1,IMXB 120 ZBM(I,J)=-300. C C READING IN BARRIER HEIGHTS AT SELECTED GRID POINTS READ (9,520) INWP,AC ! Huiqing.Liu /MDL basin file changed from 9 to 19 IF(AC.EQ.'X') INWP=INWP+1000 IF(AC.EQ.'Y') INWP=INWP+2000 IF(AC.EQ.'Z') INWP=INWP+3000 IF(AC.EQ.'W') INWP=INWP+4000 IF(AC.EQ.'V') INWP=INWP+5000 IF(AC.EQ.'U') INWP=INWP+6000 C WRITE (*,520) INWP !--------------------------------------------- ! rewrite codes to remove GO TO statments ! Huiqing.Liu /MDL Feb. 2015 !--------------------------------------------- ! IF (INWP.EQ.0) GO TO 220 IF (INWP.NE.0) THEN ! ! INWX=(INWP-1)/7+1 I7X=INWP-(INWX-1)*7 I7=7 DO 210 IP=1,INWX IF(IP.EQ.INWX) I7=I7X READ (9,550) (INW(I),JNW(I),HW(I),I=1,I7) DO 200 N=1,I7 I=INW(N) J=JNW(N) !---------------------------------------------------- ! Huiqing.Liu/MDL for using mask in NEP basin ! In order to avoid warning messages, only assgin ! cut, barrier height if cell is not blocked by mask ! if(i_boun(i,j).ne.1) zbm(i,j)=hw(n) 200 continue ! 200 ZBM(I,J)=HW(N) !----------------------------------------------------- 210 CONTINUE ! ! END IF ! ! C C COMPLETE DEFINING ZBM( , ) IN SUBROUTINE 'DEPSFC' C C READ IN DATA FOR 1-DIM FLOW, (I,J) POINTS, SIDES 1-4 220 CONTINUE C C INITIALIZE DATA FOR REGULAR 1D FLOW DO 230 I=1,L_ BANK(I,1)=0. BANK(I,2)=0. HWEIR(I)=0. DELCUT(I)=1. 230 CONTINUE C AC='X' IF 1D-FLOW POINTS EXCEED 1000 . READ (9,520) NSQRS,AC,BC C WRITE (*,520) NSQRS,AC,BC IF (AC.EQ.'X') NSQRS=1000+NSQRS IF (AC.EQ.'Y') NSQRS=2000+NSQRS IF (AC.EQ.'Z') NSQRS=3000+NSQRS IF (AC.EQ.'W') NSQRS=4000+NSQRS IF (AC.EQ.'V') NSQRS=5000+NSQRS IF (AC.EQ.'U') NSQRS=6000+NSQRS IF (AC.EQ.'T') NSQRS=7000+NSQRS IF (AC.EQ.'S') NSQRS=8000+NSQRS IF (AC.EQ.'R') NSQRS=9000+NSQRS IF (AC.EQ.'Q') NSQRS=10000+NSQRS IF (AC.EQ.'P') NSQRS=11000+NSQRS IF (AC.EQ.'O') NSQRS=12000+NSQRS IF (AC.EQ.'N') NSQRS=13000+NSQRS IF (AC.EQ.'M') NSQRS=14000+NSQRS IF (AC.EQ.'L') NSQRS=15000+NSQRS IF (AC.EQ.'K') NSQRS=16000+NSQRS IF (AC.EQ.'J') NSQRS=17000+NSQRS IF (AC.EQ.'I') NSQRS=18000+NSQRS IF (AC.EQ.'H') NSQRS=19000+NSQRS IF (AC.EQ.'G') NSQRS=20000+NSQRS IF (AC.EQ.'F') NSQRS=21000+NSQRS IF (AC.EQ.'E') NSQRS=22000+NSQRS IF (AC.EQ.'D') NSQRS=23000+NSQRS IF (AC.EQ.'C') NSQRS=24000+NSQRS IF (AC.EQ.'B') NSQRS=25000+NSQRS IF (AC.EQ.'A') NSQRS=26000+NSQRS IF (AC.EQ.'z') NSQRS=27000+NSQRS IF (AC.EQ.'y') NSQRS=28000+NSQRS IF (AC.EQ.'x') NSQRS=29000+NSQRS IF (AC.EQ.'w') NSQRS=30000+NSQRS IF (AC.EQ.'v') NSQRS=31000+NSQRS IF (AC.EQ.'u') NSQRS=32000+NSQRS IF (AC.EQ.'t') NSQRS=33000+NSQRS IF (AC.EQ.'s') NSQRS=34000+NSQRS IF (AC.EQ.'r') NSQRS=35000+NSQRS IF (AC.EQ.'q') NSQRS=36000+NSQRS C IF (AC.EQ.'p') NSQRS=37000+NSQRS C IF (NSQRS.GT.L_) THEN WRITE(*,*)' NO. OF 1D FLOW POINTS ',NSQRS WRITE(*,*)' EXCEEDS THE MAXIMUM ALLOWED. MAX= ',L_ STOP ENDIF !--------------------------------------------- ! rewrite codes to remove GO TO statments ! Huiqing.Liu /MDL Feb. 2015 !--------------------------------------------- IF (NSQRS.NE.0) THEN READ (9,540) (ISQR(I),JSQR(I),KTREE(I),ISIDE(I),IANG(I), 1 I=1,NSQRS) C C FOR BC='X' ENTER SPECIAL 1D FLOW SQUARES,BE SURE OF THE SAME C ORDER AS IN 1D FLOW IF (BC.EQ.'X') THEN C C READ IN BANK DATA READ (9,520) NBNKS,AC IF (AC.EQ.'X') NBNKS=1000+NBNKS IF (AC.EQ.'Y') NBNKS=2000+NBNKS IF (AC.EQ.'Z') NBNKS=3000+NBNKS IF (AC.EQ.'W') NBNKS=4000+NBNKS IF (AC.EQ.'V') NBNKS=5000+NBNKS IF (AC.EQ.'U') NBNKS=6000+NBNKS IF (AC.EQ.'T') NBNKS=7000+NBNKS IF (AC.EQ.'S') NBNKS=8000+NBNKS IF (AC.EQ.'R') NBNKS=9000+NBNKS IF (AC.EQ.'Q') NBNKS=10000+NBNKS C WRITE(*,520) NBNKS READ (9,250)(IBNK(I),JBNK(I),AC,ISBNK(I),BKI(I),BKE(I), 1 WCHN(I),I=1,NBNKS) 250 FORMAT (4(2I3,A1,I1,2F3.0,F4.2)) C C REPOSITION BANKS WITH RESPECT TO 1D FLOW SQUARES C LIST OF BANK POINTS MUST BE IN THE SAME ORDER AS 1D FLOW PT. LLL=1 DO I=1,NBNKS II=IBNK(I) JJ=JBNK(I) K=ISBNK(I) LL=LLL DO L=LL,NSQRS ISET=ISQR(L) IF (ISET.EQ.II) THEN ! 260 JSET=JSQR(L) IF (JSET.EQ.JJ) THEN ! 260 KSET=ISIDE(L) IF (KSET.EQ.K) THEN ! 260 LLL=L+1 BANK(L,1)=BKI(I) BANK(L,2)=BKE(I) DELCUT(L)=WCHN(I) EXIT ! 290 END IF END IF END IF END DO END DO END IF END IF ! ! ! ! IF (NSQRS.EQ.0) GO TO 300 ! READ (9,540) (ISQR(I),JSQR(I),KTREE(I),ISIDE(I),IANG(I), ! 1 I=1,NSQRS) !C !C FOR BC='X' ENTER SPECIAL 1D FLOW SQUARES,BE SURE OF THE SAME !C ORDER AS IN 1D FLOW ! IF (BC.NE.'X') GO TO 300 !C !C READ IN BANK DATA ! READ (9,520) NBNKS,AC ! IF (AC.EQ.'X') NBNKS=1000+NBNKS ! IF (AC.EQ.'Y') NBNKS=2000+NBNKS ! IF (AC.EQ.'Z') NBNKS=3000+NBNKS ! IF (AC.EQ.'W') NBNKS=4000+NBNKS ! IF (AC.EQ.'V') NBNKS=5000+NBNKS ! IF (AC.EQ.'U') NBNKS=6000+NBNKS ! IF (AC.EQ.'T') NBNKS=7000+NBNKS ! IF (AC.EQ.'S') NBNKS=8000+NBNKS ! IF (AC.EQ.'R') NBNKS=9000+NBNKS ! IF (AC.EQ.'Q') NBNKS=10000+NBNKS !C WRITE(*,520) NBNKS ! READ (9,250)(IBNK(I),JBNK(I),AC,ISBNK(I),BKI(I),BKE(I), ! 1 WCHN(I),I=1,NBNKS) ! 250 FORMAT (4(2I3,A1,I1,2F3.0,F4.2)) C C REPOSITION BANKS WITH RESPECT TO 1D FLOW SQUARES C LIST OF BANK POINTS MUST BE IN THE SAME ORDER AS 1D FLOW PT. ! LLL=1 ! DO 290 I=1,NBNKS ! II=IBNK(I) ! JJ=JBNK(I) ! K=ISBNK(I) ! LL=LLL ! DO 280 L=LL,NSQRS ! ISET=ISQR(L) ! IF (ISET.NE.II) GO TO 260 ! JSET=JSQR(L) ! IF (JSET.NE.JJ) GO TO 260 ! KSET=ISIDE(L) ! IF (KSET.NE.K) GO TO 260 ! LLL=L+1 ! BANK(L,1)=BKI(I) ! BANK(L,2)=BKE(I) ! DELCUT(L)=WCHN(I) ! GO TO 290 ! 260 IF (L.NE.NSQRS) GO TO 280 !C WRITE(*,270) IBNK(I),JBNK(I),ISBNK(I) ! 270 FORMAT(' THE SPECIAL FLOW POINT FOR BANKS AT I=',I3,' J=', ! 1I3,' SIDE=',I3,' DOES NOT EXIST IN THE 1D FLOW REGION') ! 280 CONTINUE ! 290 CONTINUE C C READ IN DATA FOR RAISED WEIRS FOR 1-DIM FLOW 300 READ (9,520) NWEIRS,AC IF (AC.EQ.'X') NWEIRS=NWEIRS+1000 IF (AC.EQ.'Y') NWEIRS=NWEIRS+2000 IF (AC.EQ.'Z') NWEIRS=NWEIRS+3000 ! ! ! IF (NWEIRS.EQ.0) GO TO 320 ! READ (9,530)(ISQR(NSQRS+I),JSQR(NSQRS+I),KTREE(NSQRS+I), ! 1 ISIDE(NSQRS+I),HWEIR(NSQRS+I),I=1,NWEIRS) ! DO 310 I=1,NWEIRS ! 310 IANG(I+NSQRS)=0 ! 320 NSQRW=NSQRS+NWEIRS IF (NWEIRS.NE.0) THEN READ (9,530)(ISQR(NSQRS+I),JSQR(NSQRS+I),KTREE(NSQRS+I), 1 ISIDE(NSQRS+I),HWEIR(NSQRS+I),I=1,NWEIRS) DO I=1,NWEIRS IANG(I+NSQRS)=0 END DO END IF NSQRW=NSQRS+NWEIRS ! ! ! C C C READ IN DATA FOR 10 SELECTED POINTS (SURGE HISTORY) READ (9,630) (STATS(I),I=1,10) C WRITE(*,630) (STATS(I),I=1,10) READ (9,610) (IPN(I) ,I=1,10) READ (9,610) (JPN(I) ,I=1,10) C C READ IN DATA FOR CHANNELS READ (9,520) NPSS ! ! ! IF (NPSS.EQ.0) GO TO 340 IF (NPSS.NE.0) THEN READ (9,650) (CHLNH(I),I=1,NPSS) READ (9,680) (CHWTH(I),I=1,NPSS) READ (9,690) (CHDPH(I),I=1,NPSS) READ (9,700) ((IPT0(I,J),I=1,2),J=1,NPSS) READ (9,700) ((JPT0(I,J),I=1,2),J=1,NPSS) READ (9,700) ((IPTL(I,J),I=1,2),J=1,NPSS) READ (9,700) ((JPTL(I,J),I=1,2),J=1,NPSS) END IF ! 340 CONTINUE ! C C READ IN DATA FOR FLOW ACROSS CUTS READ (9,520) NCUT,AC IF (AC.EQ.'X') NCUT=1000+NCUT IF (AC.EQ.'Y') NCUT=2000+NCUT NSQRWC=NSQRW+NCUT IF (NSQRWC.GT.L_) THEN WRITE(*,*)' NO. OF 1D FLOW POINTS + CUTS',NSQRWC WRITE(*,*)' EXCEEDS THE MAXIMUM ALLOWED. MAX= ',L_ STOP ENDIF C ! ! ! ! IF (NCUT.EQ.0) GO TO 390 !C ! DO 380 I=1,NCUT ! M=I+NSQRW ! READ (9,370)ISQR(M),JSQR(M),KTREE(M),ISIDE(M),BANK(M,1), ! 1 BANK(M,2),CUTL(I),CUTLI(I),CUTLE(I),IANG(M),HWEIR(M) ! 370 FORMAT(2I3,A1,I1,2F4.0,3F5.2,I3,F5.1) ! IF (CUTL(I).EQ.0.) CUTL(I)=1. ! IF (CUTLI(I).EQ.0.) CUTLI(I)=1. ! IF (CUTLE(I).EQ.0.) CUTLE(I)=1. ! 380 CONTINUE ! 390 CONTINUE IF (NCUT.NE.0) THEN DO I=1,NCUT M=I+NSQRW READ (9,370)ISQR(M),JSQR(M),KTREE(M),ISIDE(M),BANK(M,1), 1 BANK(M,2),CUTL(I),CUTLI(I),CUTLE(I),IANG(M),HWEIR(M) 370 FORMAT(2I3,A1,I1,2F4.0,3F5.2,I3,F5.1) IF (CUTL(I).EQ.0.) CUTL(I)=1. IF (CUTLI(I).EQ.0.) CUTLI(I)=1. IF (CUTLE(I).EQ.0.) CUTLE(I)=1. END DO END IF ! ! C C COMPUTE SINE AND COSINE OF THE ANGLE THAT IS USED TO C ROTATE WIND STRESS DIRECTION FROM ACTUAL TO MODEL'S 1D FLOW C DIRECTIONS. ! DO 410 N=1,NSQRWC ! IF (IANG(N).EQ.0) GO TO 400 ! ANG=IANG(N)*1.74532925E-2 ! COS1D(N)=COS(ANG) ! SIN1D(N)=SIN(ANG) ! GO TO 410 ! 400 COS1D(N)=1. ! SIN1D(N)=0. ! 410 CONTINUE ! 420 CONTINUE ! DO N=1,NSQRWC IF (IANG(N).EQ.0) THEN COS1D(N)=1. SIN1D(N)=0. ELSE ANG=IANG(N)*1.74532925E-2 COS1D(N)=COS(ANG) SIN1D(N)=SIN(ANG) END IF END DO ! ! C READ IN DATA FOR MANGROVE POINTS OR TREE LOCATIONS (MGPT=-1) READ (9,520) MGPT ! ! ! IF (MGPT.EQ.0) GOTO 470 !C READ IN LOCATIONS OF TREE POINTS (0 OR 1) ! DO 460 J=1,JMXB ! READ (9,450) (ITREE(I,J),I=1,IMXB) ! 450 FORMAT(60A1) ! 460 CONTINUE !C SET UP BOUNDARY VALUES TO ZBM( , ) TO BE CHECKED FOR OPEN !C BOUNDARY CALCULATIONS. ! IF (DOLLAR.NE.SIG_) THEN ! DO 160 J=1,JMXB,JMXB1 ! DO 160 I=1,IMXB ! ZBM(I,J)=100. ! 160 CONTINUE ! ENDIF ! DO 180 I=1,IMXB,IMXB1 ! DO 180 J=1,JMXB ! ZBM(I,J)=100. ! 180 CONTINUE ! ! IF (MGPT.NE.0) THEN C READ IN LOCATIONS OF TREE POINTS (0 OR 1) DO J=1,JMXB READ (9,450) (ITREE(I,J),I=1,IMXB) 450 FORMAT(60A1) END DO C SET UP BOUNDARY VALUES TO ZBM( , ) TO BE CHECKED FOR OPEN C BOUNDARY CALCULATIONS. IF (DOLLAR.NE.SIG_) THEN DO J=1,JMXB,JMXB1 DO I=1,IMXB ZBM(I,J)=100. END DO END DO ENDIF DO I=1,IMXB,IMXB1 DO J=1,JMXB ZBM(I,J)=100. END DO END DO END IF ! ! C C READ IN DATA FOR PRINT OUT OF SHEETS OF MAX SURGES C OPTIONAL FORMATS TO ACCOMMODATE DIFFERENT BASINS DATA C FORMAT-- I3,A1 C IX60--' ' FOR DISPLAYING I=1,IXB,WHERE IXB < 60 C 'X' FOR DISPLAYING I=1,IXB,WHERE IXB >=60 C 'Y' FOR I=1,IPGF(I) C 'Z' FOR I=IPGS(I),IPGF(I) MOST GENERAL 470 READ (9,520) NSHEET,IX60 NSHEET=MIN(NSHEET,6) C WRITE(*,520) NSHEET,IX60 READ (9,610) (JPGS(I),I=1,NSHEET) C WRITE(*,610) (JPGS(I),I=1,NSHEET) READ (9,610) (JPGF(I),I=1,NSHEET) C WRITE(*,610) (JPGF(I),I=1,NSHEET) DO 480 N=1,NSHEET IPGS(N)=1 IPGF(N)=MIN0(59,IMXB1) IF (IX60.EQ.'X') IPGF(N)=IMXB1 480 CONTINUE IF (IX60.EQ.'Z') THEN READ (9,610) (IPGS(I),I=1,NSHEET) C WRITE (*,610) (IPGS(I),I=1,NSHEET) ENDIF IF (IX60.EQ.YA1.OR.IX60.EQ.ZA1) THEN READ (9,610) (IPGF(I),I=1,NSHEET) C WRITE (*,610) (IPGF(I),I=1,NSHEET) ENDIF C C READING IN DEPTH VALUES IMX=(IMXB-1)/15+1 DO 500 K=1,IMX IST=(K-1)*15+1 IFN=MIN0(IST+14,IMXB) DO 490 J=1,JMXB READ (9,710) (ZB(I,J),I=IST,IFN) CC IF (J.EQ.JMXB) WRITE(*,710) (ZB(I,J),I=IST,IFN) ! Huiqing.Liu /MDL Remove the depth limit ! DO 405 I=IST,IFN ! 405 ZB(I,J)=MIN(600.,ZB(I,J)) ! 405 ZB(I,J)=MIN(1500.,ZB(I,J)) ! 490 CONTINUE 500 CONTINUE c------------ FOR BIX BASIN ONLY IF (DOLLAR.EQ.'1$') THEN READ (9,'(i4)') NODRY IF(NODRY.GT.ND_) THEN WRITE(*,*) "Too many Dry cells" STOP ENDIF IF(NODRY.GT.1)THEN READ (9,'(9(I3,1X,I3,1X))')(IDRY(K),JDRY(K),K=1,NODRY) ENDIF ENDIF C IF (DOLLAR.EQ.'$') THEN DO 510 I=1,IMXB 510 ZB(I,JMXB)=ZB(I,1) ELSE DO 512 I=1,IMXB 512 ZB(I,JMXB)=ZB(I,JMXB1) ENDIF C CLOSE (9) 520 FORMAT (I3,2A1) 530 FORMAT (5(2I3,A1,I1,F5.1)) 540 FORMAT (6(2I3,A1,I1,I3)) 550 FORMAT (7(I3,I3,F3.0)) 560 FORMAT (A10,2X,A1) 570 FORMAT (F10.6) 580 FORMAT (2F10.6) 590 FORMAT (4I5) 600 FORMAT (22I3) 610 FORMAT (15I3) 620 FORMAT (45I1) 630 FORMAT (6A10) 640 FORMAT (15F3.0) 650 FORMAT (5F5.1) 660 FORMAT (13F5.0) 670 FORMAT (A4,3I3,1X,3I3) 680 FORMAT (5F6.0) 690 FORMAT (5F5.0) 700 FORMAT (10I3) 710 FORMAT (15F4.0,A16) RETURN END SUBROUTINE CRDRD2_N INCLUDE 'parm.for' C INCLUDE 'parmmsy.for' C COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /BSN/ PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C PARAMETER (NT_=160) COMMON /SPLN/ ALT(15),ALN(15),AX(15),AY(15),RL(15),ANGD(15), 1 PT(15),RT(15),XLAT(NT_),YLONG(NT_),RLNGTH COMMON /DUMBSV/ SAVEH(10),IPN(10),JPN(10),STATS(10) COMMON /PRTOVL/ NSHEET,JPGS(10),JPGF(10),IPGS(10),IPGF(10) COMMON /FLWCPT/ NSQRS,NSQRW,NSQRWC,NPSS,NCUT COMMON /DLTDTA/ DLTIN(3),DLTB COMMON /CHNLD/ CHLNH(5),CHWTH(5),CHDPH(5) COMMON /CHNL/ IPT0(2,5),JPT0(2,5),IPTL(2,5),JPTL(2,5),ENTEXT(5) COMMON /IPRX/ IX60 COMMON /GPRT1/ DOLLAR,EBSN CHARACTER*1 EBSN !------------------------------------------------------- ! Added Basin Name Variable by Huiqing.Liu/MDL June/2015 ! for Modifying Cd wind drag coefficient parameter only ! in New Bering Sea extra-tropical basin (eno3) !------------------------------------------------------- COMMON /GPRT/ STA CHARACTER*16 STA C CHARACTER*16 STATS DIMENSION IANG(L_) CHARACTER*16 AAA DIMENSION INW(7),JNW(7),HW(7) DIMENSION IBNK(NBK_),JBNK(NBK_),ISBNK(NBK_),BKI(NBK_), 1 BKE(NBK_),WCHN(NBK_) CHARACTER*2 AC,BC,IX60,XA1,YA1,ZA1,DOLLAR,SIG_ DATA XA1/'X'/,YA1/'Y'/,ZA1/'Z'/,SIG_/'$'/ C C READ BOUNDARY TYPES ALONG LEFT,RIGHT,TOP AND BOTTOM BOUNDARIES C ALONG MOMENTUM POINTS. C .****.****.****.****.****.****.****.****. TOP BNDRY C * * C * * C . . . . . . . . . C LEFT BNDY * * RIGHT BNDY C * + + + + * C * * C .****.****.****.****.****.****.****.****. BOTTOM BNDRY C 2 8 8 8 9 9 9 1 1 C DEEP WATER-) (-INTERMEDIATE----) (-SHALLOW----) (--LAND---- C IF (IMXB.GT.M_.OR.JMXB.GT.N_) THEN WRITE (*,*) ' BASIN DIMENSIONS ',IMXB,JMXB WRITE (*,*) ' HAS EXCEEDS THE COMPILE DIMENSIONS',M_,N_ STOP endif C INITIALIZE ZBM(I,J) FOR ENTIRE GRIDS DO 120 J=1,JMXB DO 120 I=1,IMXB 120 ZBM(I,J)=-300. C C READING IN BARRIER HEIGHTS AT SELECTED GRID POINTS READ (9,520) INWP,AC IF(AC.EQ.'X') INWP=INWP+1000 IF(AC.EQ.'Y') INWP=INWP+2000 IF(AC.EQ.'Z') INWP=INWP+3000 IF(AC.EQ.'W') INWP=INWP+4000 IF(AC.EQ.'V') INWP=INWP+5000 IF(AC.EQ.'U') INWP=INWP+6000 IF(AC.EQ.'T') INWP=INWP+7000 ! WRITE (*,*) INWP,AC !--------------------------------------------- ! rewrite codes to remove GO TO statments ! Huiqing.Liu /MDL July. 2020 !--------------------------------------------- ! IF (INWP.EQ.0) GO TO 220 IF (INWP.NE.0) THEN INWX=(INWP-1)/7+1 I7X=INWP-(INWX-1)*7 I7=7 DO 210 IP=1,INWX IF(IP.EQ.INWX) I7=I7X READ (9,550) (INW(I),JNW(I),HW(I),I=1,I7) DO 200 N=1,I7 I=INW(N) J=JNW(N) 200 ZBM(I,J)=HW(N) 210 CONTINUE ENDIF C C COMPLETE DEFINING ZBM( , ) IN SUBROUTINE 'DEPSFC' C C READ IN DATA FOR 1-DIM FLOW, (I,J) POINTS, SIDES 1-4 220 CONTINUE C C INITIALIZE DATA FOR REGULAR 1D FLOW DO 230 I=1,L_ BANK(I,1)=0. BANK(I,2)=0. HWEIR(I)=0. DELCUT(I)=1. 230 CONTINUE C AC='X' IF 1D-FLOW POINTS EXCEED 1000 . READ (9,520) NSQRS,AC,BC ! WRITE (*,*) NSQRS,AC,BC IF (AC.EQ.'X') NSQRS=1000+NSQRS IF (AC.EQ.'Y') NSQRS=2000+NSQRS IF (AC.EQ.'Z') NSQRS=3000+NSQRS IF (AC.EQ.'W') NSQRS=4000+NSQRS IF (AC.EQ.'V') NSQRS=5000+NSQRS IF (AC.EQ.'U') NSQRS=6000+NSQRS IF (AC.EQ.'T') NSQRS=7000+NSQRS IF (AC.EQ.'S') NSQRS=8000+NSQRS IF (AC.EQ.'R') NSQRS=9000+NSQRS IF (AC.EQ.'Q') NSQRS=10000+NSQRS IF (AC.EQ.'P') NSQRS=11000+NSQRS IF (AC.EQ.'O') NSQRS=12000+NSQRS IF (AC.EQ.'N') NSQRS=13000+NSQRS IF (AC.EQ.'M') NSQRS=14000+NSQRS IF (AC.EQ.'L') NSQRS=15000+NSQRS IF (AC.EQ.'K') NSQRS=16000+NSQRS IF (AC.EQ.'J') NSQRS=17000+NSQRS IF (AC.EQ.'I') NSQRS=18000+NSQRS IF (AC.EQ.'H') NSQRS=19000+NSQRS IF (AC.EQ.'G') NSQRS=20000+NSQRS IF (AC.EQ.'F') NSQRS=21000+NSQRS IF (AC.EQ.'E') NSQRS=22000+NSQRS IF (AC.EQ.'D') NSQRS=23000+NSQRS IF (AC.EQ.'C') NSQRS=24000+NSQRS IF (AC.EQ.'B') NSQRS=25000+NSQRS IF (AC.EQ.'A') NSQRS=26000+NSQRS IF (AC.EQ.'z') NSQRS=27000+NSQRS IF (AC.EQ.'y') NSQRS=28000+NSQRS IF (AC.EQ.'x') NSQRS=29000+NSQRS IF (AC.EQ.'w') NSQRS=30000+NSQRS IF (AC.EQ.'v') NSQRS=31000+NSQRS IF (AC.EQ.'u') NSQRS=32000+NSQRS IF (AC.EQ.'t') NSQRS=33000+NSQRS IF (AC.EQ.'s') NSQRS=34000+NSQRS IF (AC.EQ.'r') NSQRS=35000+NSQRS IF (AC.EQ.'q') NSQRS=36000+NSQRS IF (AC.EQ.'p') NSQRS=37000+NSQRS IF (AC.EQ.'o') NSQRS=38000+NSQRS IF (AC.EQ.'n') NSQRS=39000+NSQRS C IF (NSQRS.GT.L_) THEN WRITE(*,*)' NO. OF 1D FLOW POINTS ',NSQRS WRITE(*,*)' EXCEEDS THE MAXIMUM ALLOWED. MAX= ',L_ STOP ENDIF ! IF (NSQRS.EQ.0) GO TO 300 IF (NSQRS.NE.0) THEN READ (9,540) (ISQR(I),JSQR(I),KTREE(I),ISIDE(I),IANG(I), 1 I=1,NSQRS) C C FOR BC='X' ENTER SPECIAL 1D FLOW SQUARES,BE SURE OF THE SAME C ORDER AS IN 1D FLOW C IF (BC.EQ.'X') THEN C READ IN BANK DATA READ(9,520) NBNKS,AC IF (AC.EQ.'X') NBNKS=1000+NBNKS IF (AC.EQ.'Y') NBNKS=2000+NBNKS IF (AC.EQ.'Z') NBNKS=3000+NBNKS IF (AC.EQ.'W') NBNKS=4000+NBNKS IF (AC.EQ.'V') NBNKS=5000+NBNKS IF (AC.EQ.'U') NBNKS=6000+NBNKS IF (AC.EQ.'T') NBNKS=7000+NBNKS IF (AC.EQ.'S') NBNKS=8000+NBNKS IF (AC.EQ.'R') NBNKS=9000+NBNKS IF (AC.EQ.'Q') NBNKS=10000+NBNKS C WRITE(*,520) NBNKS READ(9,250)(IBNK(I),JBNK(I),AC,ISBNK(I),BKI(I),BKE(I), 1 WCHN(I),I=1,NBNKS) 250 FORMAT (4(2I4,A1,I1,2F3.0,F4.2)) C C REPOSITION BANKS WITH RESPECT TO 1D FLOW SQUARES C LIST OF BANK POINTS MUST BE IN THE SAME ORDER AS 1D FLOW PT. LLL=1 DO I=1,NBNKS II=IBNK(I) JJ=JBNK(I) K=ISBNK(I) LL=LLL DO L=LL,NSQRS ISET=ISQR(L) IF (ISET.EQ.II) THEN ! 260 JSET=JSQR(L) IF (JSET.EQ.JJ) THEN ! 260 KSET=ISIDE(L) IF (KSET.EQ.K) THEN ! 260 LLL=L+1 BANK(L,1)=BKI(I) BANK(L,2)=BKE(I) DELCUT(L)=WCHN(I) EXIT ! 290 END IF END IF END IF END DO END DO END IF END IF C C READ IN DATA FOR RAISED WEIRS FOR 1-DIM FLOW 300 READ (9,520) NWEIRS,AC IF (AC.EQ.'X') NWEIRS=NWEIRS+1000 IF (AC.EQ.'Y') NWEIRS=NWEIRS+2000 IF (AC.EQ.'Z') NWEIRS=NWEIRS+3000 IF (NWEIRS.NE.0) THEN READ (9,530)(ISQR(NSQRS+I),JSQR(NSQRS+I),KTREE(NSQRS+I), 1 ISIDE(NSQRS+I),HWEIR(NSQRS+I),I=1,NWEIRS) DO I=1,NWEIRS IANG(I+NSQRS)=0 END DO END IF NSQRW=NSQRS+NWEIRS C C C READ IN DATA FOR 10 SELECTED POINTS (SURGE HISTORY) READ (9,630) (STATS(I),I=1,10) C WRITE(*,630) (STATS(I),I=1,10) READ (9,610) (IPN(I) ,I=1,10) READ (9,610) (JPN(I) ,I=1,10) C C READ IN DATA FOR CHANNELS READ (9,520) NPSS IF (NPSS.NE.0) THEN READ (9,650) (CHLNH(I),I=1,NPSS) READ (9,680) (CHWTH(I),I=1,NPSS) READ (9,690) (CHDPH(I),I=1,NPSS) READ (9,700) ((IPT0(I,J),I=1,2),J=1,NPSS) READ (9,700) ((JPT0(I,J),I=1,2),J=1,NPSS) READ (9,700) ((IPTL(I,J),I=1,2),J=1,NPSS) READ (9,700) ((JPTL(I,J),I=1,2),J=1,NPSS) END IF C C READ IN DATA FOR FLOW ACROSS CUTS READ (9,520) NCUT,AC IF (AC.EQ.'X') NCUT=1000+NCUT IF (AC.EQ.'Y') NCUT=2000+NCUT NSQRWC=NSQRW+NCUT IF (NSQRWC.GT.L_) THEN WRITE(*,*)' NO. OF 1D FLOW POINTS + CUTS',NSQRWC WRITE(*,*)' EXCEEDS THE MAXIMUM ALLOWED. MAX= ',L_ STOP ENDIF C C IF (NCUT.NE.0) THEN DO 380 I=1,NCUT M=I+NSQRW READ (9,370)ISQR(M),JSQR(M),KTREE(M),ISIDE(M),BANK(M,1),BANK(M,2), 1 CUTL(I),CUTLI(I),CUTLE(I),IANG(M),HWEIR(M) 370 FORMAT(2I4,A1,I1,2F4.0,3F5.2,I3,F5.1) IF (CUTL(I).EQ.0.) CUTL(I)=1. IF (CUTLI(I).EQ.0.) CUTLI(I)=1. IF (CUTLE(I).EQ.0.) CUTLE(I)=1. 380 CONTINUE ENDIF C COMPUTE SINE AND COSINE OF THE ANGLE THAT IS USED TO C ROTATE WIND STRESS DIRECTION FROM ACTUAL TO MODEL'S 1D FLOW C DIRECTIONS. DO N=1,NSQRWC IF (IANG(N).EQ.0) THEN COS1D(N)=1. SIN1D(N)=0. ELSE ANG=IANG(N)*1.74532925E-2 COS1D(N)=COS(ANG) SIN1D(N)=SIN(ANG) END IF END DO C READ IN DATA FOR MANGROVE POINTS OR TREE LOCATIONS (MGPT=-1) READ (9,520) MGPT IF (MGPT.NE.0) THEN C READ IN LOCATIONS OF TREE POINTS (0 OR 1) DO 460 J=1,JMXB READ (9,450) (ITREE(I,J),I=1,IMXB) 450 FORMAT(60A1) 460 CONTINUE C SET UP BOUNDARY VALUES TO ZBM( , ) TO BE CHECKED FOR OPEN C BOUNDARY CALCULATIONS. IF (DOLLAR.NE.SIG_) THEN DO 160 J=1,JMXB,JMXB1 DO 160 I=1,IMXB ZBM(I,J)=100. 160 CONTINUE ENDIF DO 180 I=1,IMXB,IMXB1 DO 180 J=1,JMXB ZBM(I,J)=100. 180 CONTINUE ENDIF C C READ IN DATA FOR PRINT OUT OF SHEETS OF MAX SURGES C OPTIONAL FORMATS TO ACCOMMODATE DIFFERENT BASINS DATA C FORMAT-- I3,A1 C IX60--' ' FOR DISPLAYING I=1,IXB,WHERE IXB < 60 C 'X' FOR DISPLAYING I=1,IXB,WHERE IXB >=60 C 'Y' FOR I=1,IPGF(I) C 'Z' FOR I=IPGS(I),IPGF(I) MOST GENERAL 470 READ (9,520) NSHEET,IX60 NSHEET=MIN(NSHEET,6) C WRITE(*,520) NSHEET,IX60 READ (9,610) (JPGS(I),I=1,NSHEET) C WRITE(*,610) (JPGS(I),I=1,NSHEET) READ (9,610) (JPGF(I),I=1,NSHEET) C WRITE(*,610) (JPGF(I),I=1,NSHEET) DO 480 N=1,NSHEET IPGS(N)=1 IPGF(N)=MIN0(59,IMXB1) IF (IX60.EQ.'X') IPGF(N)=IMXB1 480 CONTINUE IF (IX60.EQ.'Z') THEN READ (9,610) (IPGS(I),I=1,NSHEET) C WRITE (*,610) (IPGS(I),I=1,NSHEET) ENDIF IF (IX60.EQ.YA1.OR.IX60.EQ.ZA1) THEN READ (9,610) (IPGF(I),I=1,NSHEET) C WRITE (*,610) (IPGF(I),I=1,NSHEET) ENDIF C C READING IN DEPTH VALUES IMX=(IMXB-1)/15+1 DO 500 K=1,IMX IST=(K-1)*15+1 IFN=MIN0(IST+14,IMXB) DO 490 J=1,JMXB READ (9,710) (ZB(I,J),I=IST,IFN) CC IF (J.EQ.JMXB) WRITE(*,710) (ZB(I,J),I=IST,IFN) ! DO 405 I=IST,IFN ! 405 ZB(I,J)=MIN(600.,ZB(I,J)) 490 CONTINUE 500 CONTINUE IF (DOLLAR.EQ.'1$') THEN READ (9,'(i4)') NODRY IF(NODRY.GT.ND_) THEN WRITE(*,*) "Too many Dry cells" STOP ENDIF IF(NODRY.GT.1)THEN READ (9,'(9(I4,1X,I4,1X))')(IDRY(K),JDRY(K),K=1,NODRY) ENDIF ENDIF C IF (DOLLAR.EQ.'$') THEN DO 510 I=1,IMXB 510 ZB(I,JMXB)=ZB(I,1) ELSE DO 512 I=1,IMXB 512 ZB(I,JMXB)=ZB(I,JMXB1) ENDIF C CLOSE (9) 520 FORMAT (I3,2A1) 530 FORMAT (5(2I3,A1,I1,F5.1)) 540 FORMAT (6(2I4,A1,I1,I3)) 550 FORMAT (7(I4,I4,F3.0)) 560 FORMAT (A10,2X,A1) 570 FORMAT (F10.6) 580 FORMAT (2F10.6) 590 FORMAT (4I5) 600 FORMAT (22I3) 610 FORMAT (15I3) 620 FORMAT (45I1) 630 FORMAT (6A10) 640 FORMAT (15F3.0) 650 FORMAT (5F5.1) 660 FORMAT (13F5.0) 670 FORMAT (A4,3I3,1X,3I3) 680 FORMAT (5F6.0) 690 FORMAT (5F5.0) 700 FORMAT (10I3) 710 FORMAT (15F5.0,A16) RETURN END SUBROUTINE DEPSFC C JELESNIANSKI SEPTEMBER 1980 TDL IBM 360/195 C C PURPOSE C THIS SUBROUTINE (DEPSFC) INITIALIZES INLAND SURFACE WATERS C TO LAKE DATUM (SEA WATER TO SEA DATUM IS DONE IN SUBROU- C TINE 'INILHT'). BARRIERS AND THEIR HEIGHTS ARE PLACED AT C SELECTED MOMENTUM GRID POINTS. FOR 1-DIM FLOW, DEPTHS AND C LATERAL BOUNDARIES ARE SET. BOUNDARIES ARE TEMPORARILY C SET AS A 100 FOOT WALL. PRINT OUTS ARE MADE OF THE DEPTH C FIELD, AND HIGHEST LAND VALUE AT MOMENTUM GRID POINTS C ARE SET. C C DATA SET USE C NONE C C VARIABLES C IMXB JMXB = TOTAL GRID SPACINGS, SET IN 'CRDRD1' C IMXB1 JMXB1 = IMXB LESS 1, JMXB LESS 1 C DOLLAR = '$' FOR ISLAND, OTHERWISE FOR REGULAR BASIN C ZB( , ) = BATHY/TOPO HEIGHTS, RELATIVE TO NGVD C ZBM( , ) = HIGHEST ELEVATION AT A MOMENTUM POINT C OR, VALUES INDICATING BOUNDARY TYPES C NSQRW = NUMBER OF SQUARES WITH 1-DIM FLOW C ISQR JSQR( ) = SUBSCRIPTS FOR SQUARES WITH 1-DIM FLOW C ISIDE( ) = THE SIDE OF A SQUARE TO PASS 1-DIM FLOW C IIH(4) JJH(4) = SHIFT SUBSCRIPTS TO 4 MOMENTUM POINTS ON A SQUARE C HWEIR( ) = LOWESTT BARRIER ELEVATION AT ISIDE( ) C IZ1(4) IZ2(4) = SET P SUBSCRIPTS ON 4 SIDES OF ISIDE( ) C JZ1(4) JZ2(4) = SET Q SUBSCRIPTS ON 4 SIDES OF ISIDE( ) C ZBMIN( ) = LOWEST DEPTH OF TWO SQUARES BOUNDING ISIDE(30 C MGPT = TOTAL NUMBER OF MANGROVE PTS C IMG( ) JMG( ) = I/J COORDINATES OF MANGROVE PTS C ITREE( , ) = 0 FOR NO TREE, 1 FOR TREE ON A GRID POINT C C GENERAL COMMENTS C THIS SUBROUTINE RESIDES IN OVERLAY 'INITLZ'. IT IS C CALLED IN SUBROUTINE 'INITLZ'. C INCLUDE 'parm.for' C PARAMETER (NBCPTS=12000) common /opts/ nofld,nof1d character*1 nofld,nof1d COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /DUMB8/ IP(4),JP(4),IH(4),JH(4) COMMON /DUM88/ IIH(4),JJH(4),IHH(4),JHH(4) COMMON /BSN/ PHI,ALTO,ALNO,PHI1,ALT1,ALN1,ALT1C COMMON /SPLN/ ALT(15),ALN(15),AX(15),AY(15),RL(15),ANGD(15), 1 PT(15),RT(15),XLAT(999),YLONG(999),RLNGTH COMMON /BCPTS/ TIDESH(NBCPTS),BOUN_NEST(NBCPTS), !Added by Huiqing Liu /MDL Nov. 2013 * NBCPT,ISH(NBCPTS),JSH(NBCPTS) COMMON /FLWCPT/ NSQRS,NSQRW,NSQRWC,NPSS,NCUT COMMON /SWTCH/ IOPERL(5) COMMON /GPRT1/ DOLLAR,EBSN COMMON /HTERAIN/ HTER COMMON /WETDRY/ IFLD C New Common Block for controling wet&dry Huiqing Liu /MDL Jan. 2014 !------------------------------------------------------- ! Added Basin Name Variable by Huiqing.Liu/MDL June/2015 ! for Modifying Cd wind drag coefficient parameter only ! in New Bering Sea extra-tropical basin (eno3) !------------------------------------------------------- COMMON /GPRT/ STA CHARACTER*16 STA CHARACTER*2 DOLLAR CHARACTER*1 EBSN DIMENSION IZ1(4),IZ2(4),JZ1(4),JZ2(4) DATA IZ1/0,1,0,0/, IZ2/0,1,1,1/, JZ1/1,1,0,1/, JZ2/0,0,0,1/ C DO 401 J=2,JMXB DO 401 I=2,IMXB ZMX=AMAX1(-ZB(I-1,J-1),-ZB(I,J-1),-ZB(I-1,J),-ZB(I,J)) cd Huiqing Liu /MDL need to change ZMX .LT.-4999. ! IF (ZMX.LT.-1499.) ITREE(I,J)='6' ! IF (ZMX.LT.-9999.) ITREE(I,J)='6' ! IF (ZMX.LT.-4900.) ITREE(I,J)='6' if (STA(1:7) == 'AK GULF'.or.STA(1:10) == 'NEP AK GUL')then if (i_boun(i,j).eq.1) ITREE(I,J)='6' elseif (STA(2:10) == 'EXTRATROP'.or.STA(1:10) =='GULF of M')then IF (ZMX.LT.-2999.) ITREE(I,J)='6' else IF (ZMX.LT.-9999.) ITREE(I,J)='6' endif ! Huiqing.Liu /MDL added new codes to setup boundary line ! No flooding is allowed ! IF (ZMX.GT.0.) ITREE(I,J)='4' 401 CONTINUE ! Huiqing.Liu /MDL added new codes to setup boundary line ! do i=470,imxb ! do j=1,jmxb ! ITREE(I,J)='6' ! enddo ! enddo ! ! C IF (DOLLAR.EQ.'$') THEN C DO I=2,IMXB C ZMX=AMAX1(-ZB(I-1,JMXB),-ZB(I,JMXB),-ZB(I-1,1),-ZB(I,1)) C IF (ZMX.LE.-289.) ITREE(I,1)='6' C ENDDO C ENDIF C C DETERMINE THE ACTIVE/NONACTIVE SQUARES FOR HEIGHT COMPUTATIONS DO 400 J=1,JMXB1 DO 400 I=1,IMXB1 C CHANGED BY NSM TO -HTER FROM -35 11/20/2010 : Accepted 4/4/2011 !------------------------------------------------------------ ! rewrite the following codes to remove the computed GOTO ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! old codes !--------------------------------------------------------------- ! IF (ZB(I,J).LE.-HTER) GOTO 1417 ! KSKP(I,J)='1' ! KCT=1 !C SPECIAL TRUNCATED SQUARES IN DEEP OCEAN ! IF (ITREE(I,J).EQ.'6') KCT=KCT+1 ! IF (ITREE(I+1,J).EQ.'6') KCT=KCT+1 ! IF (ITREE(I,J+1).EQ.'6') KCT=KCT+1 ! IF (ITREE(I+1,J+1).EQ.'6') KCT=KCT+1 ! GOTO (400,415,415,415,417),KCT ! 415 KSKP(I,J)='2' ! GOTO 400 ! 417 ZB(I,J)=-67. ! 1417 KSKP(I,J)='0' ! 400 CONTINUE !-------------------------------------------------------------- ! New codes !--------------------------------------------------------------- IF (ZB(I,J) < -HTER) then KSKP(I,J) = '0' ELSE KSKP(I,J) = '1' KCT = 1 !C SPECIAL TRUNCATED SQUARES IN DEEP OCEAN IF (ITREE(I,J).EQ.'6') KCT=KCT+1 IF (ITREE(I+1,J).EQ.'6') KCT=KCT+1 IF (ITREE(I,J+1).EQ.'6') KCT=KCT+1 IF (ITREE(I+1,J+1).EQ.'6') KCT=KCT+1 IF (KCT == 2 .OR. KCT == 3 .OR. KCT == 4) THEN KSKP(I,J) = '2' ELSEIF (KCT == 5) THEN ZB(I,J) = -67. KSKP(I,J) = '0' ENDIF ENDIF ! ! 400 CONTINUE !------------------------------------------------------------------ DO 421 J=1,JMXB 421 KSKP(IMXB,J)='0' DO 420 I=1,IMXB 420 KSKP(I,JMXB)=KSKP(I,1) C C C ISOLATE HIGHEST LAND SQUARE ABOUT A TRANSPORT POINT C FOR ISLAND ZB( ,JMXB)=ZB( ,1) J2=JMXB1 IF (DOLLAR.EQ.'$') J2=JMXB ! ! ! DO 130 I=2,IMXB ! DO 130 J=2,J2 ! ZMX=AMAX1(-ZB(I-1,J-1),-ZB(I,J-1),-ZB(I-1,J),-ZB(I,J)) ! IF (ZBM(I,J).GE.ZMX) GO TO 130 ! IF (ZBM(I,J).NE.-300.) GO TO 110 ! ZBM(I,J)=ZMX ! GO TO 130 ! 110 WRITE(*,120)I,J,ZBM(I,J) ! 120 FORMAT (' ERROR IN BARRIER HEIGHT AT POSITION I=',I5,' J=',I5, ! 1 ' HW=',F5.0) ! ZBM(I,J)=ZMX ! 130 CONTINUE ! DO I=2,IMXB DO J=2,J2 ZMX=AMAX1(-ZB(I-1,J-1),-ZB(I,J-1),-ZB(I-1,J),-ZB(I,J)) IF (ZBM(I,J).LT.ZMX) THEN IF (ZBM(I,J).NE.-300.) THEN WRITE(*,120)I,J,ZBM(I,J) 120 FORMAT (' WARNING IN BARRIER HEIGHT AT POSITION I=', 1 I5,' J=',I5,' HW=',F5.0) END IF ZBM(I,J)=ZMX END IF END DO END DO ! ! ! IF (DOLLAR.EQ.'$') THEN DO 132 I=1,IMXB ZBM(I,1)=ZBM(I,JMXB) 132 CONTINUE ENDIF C DO 1401 J=2,JMXB DO 1401 I=2,IMXB IF (ITREE(I,J).EQ.'5'.AND.ZBM(I,J).GE.-15.) ITREE(I,J)='3' C ADDED BY AAT 4/8/2011 : This should cause what used to be high C terrain (> 35 feet) that was marked with a '4' to have friction C winds and be able to be inundated. IF (HTER.GT.35.) THEN IF (ITREE(I,J).EQ.'4'.AND.ZBM(I,J).LT.HTER) ITREE(I,J)='1' ENDIF 1401 CONTINUE C C IF ZBM GREATER THAN 35 FT IN INTERIOR POINTS, MOMENTURM PTS C ARE NONACTIVE. DO 40 J=2,J2 DO 40 I=2,IMXB1 C CHANGED BY NSM TO HTER FROM 35 11/20/2010 : Accepted 4/4/2011 IF (ZBM(I,J).GE.HTER) ITREE(I,J)='4' IF (ITREE(I,J).EQ.'5') ITREE(I,J)='2' 40 CONTINUE C C REDEFINE MS(J). FOR EACH J, FIRST I THE COMPUTATION C STARTS IN MOMNTM (ZBM LESS THAN 35 FT). !------------------------------------------------------------ ! rewrite the following codes to remove the GOTO statments ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! old codes !--------------------------------------------------------------- ! DO 50 J=2,JMXB ! MS(J)=1 ! DO 60 I=1,IMXB !C CHANGED BY NSM TO HTER FROM 35 11/20/2010 : Accepted 4/4/2011 ! IF (ZBM(I,J).LT.HTER) GOTO 55 ! 60 CONTINUE ! I=I-1 ! 55 MS(J)=I ! 50 CONTINUE !-------------------------------------------------------------- ! New codes !--------------------------------------------------------------- DO J=2,JMXB MS(J)=1 DO I=1,IMXB C CHANGED BY NSM TO HTER FROM 35 11/20/2010 : Accepted 4/4/2011 IF (ZBM(I,J).LT.HTER) THEN MS(J)=I EXIT ELSE IF (I == IMXB) THEN MS(J) = I-1 EXIT END IF END IF END DO END DO ! ! ! IF (DOLLAR.EQ.'$') THEN MS(1)=MS(JMXB) ELSE MS(1)=MS(2) ENDIF c nxx=(jmxb-1)/15+1 do n=1,nxx j1=(n-1)*15+1 j2=MIN(JMXB,j1+14) CC write (*,'(15i3)') (ms(j),j=j1,j2) enddo c !------------------------------------------------------------ ! rewrite the following codes to remove the GOTO statments ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! old codes !--------------------------------------------------------------- C DEFINE IS(J). FOR EACH J, FIRST I THE COMPUTATION C STARTS IN CONTINUITY. ! DO 10 J=1,JMXB1 ! IS(J)=1 ! DO 20 I=1,IMXB1 ! IF (KSKP(I,J).EQ.'1') GOTO 15 ! 20 CONTINUE ! I=I-1 ! 15 IS(J)=I ! 10 CONTINUE !-------------------------------------------------------------- ! New codes !--------------------------------------------------------------- DO J=1,JMXB1 IS(J)=1 DO I=1,IMXB1 IF (KSKP(I,J).EQ.'1') THEN IS(J)=I EXIT ELSE IF (I == IMXB1) THEN IS(J) = I-1 EXIT END IF END IF END DO END DO ! ! ! IF (DOLLAR.EQ.'$')THEN IS(JMXB)=IS(1) ELSE IS(JMXB)=IS(JMXB1) ENDIF !------------------------------------------------------------ ! rewrite the following codes to remove the GOTO statments ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! New codes !--------------------------------------------------------------- DO J=2,JMXB1 ME(J)=IMXB DO II=2,IMXB1 I=IMXB-II+1 C CHANGED BY NSM TO HTER FROM 35 11/20/2010 : Accepted 4/4/2011 IF (ZBM(I,J).GE.HTER) THEN CYCLE END IF IF (ITREE(I,J).NE.'6') THEN ME(J)=I EXIT ELSE IF (II == IMXB1) THEN ME(J)=I+1 EXIT END IF END IF END DO END DO ! ! ME(1)=IMXB1 ME(JMXB)=IMXB1 !------------------------------------------------------------ ! rewrite the following codes to remove the GOTO statments ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! old codes !--------------------------------------------------------------- C DEFINE LOWER LIMIT OF I FOR CONTTY POINTS C ! DO 610 J=1,JMXB1 ! IE(J)=IMXB1 ! DO 620 II=1,IMXB1 ! I=IMXB1-II+1 ! IF (KSKP(I,J).EQ.'1') GOTO 615 ! 620 CONTINUE ! I=I+1 ! 615 IE(J)=I ! 610 CONTINUE ! IE(JMXB)=IE(1) !C !C ! DO 42 J=2,JMXB ! MF(J)=IMXB ! DO 45 II=1,IMXB1 ! I=IMXB-II ! IF (ITREE(I,J).EQ.'1'.OR.ITREE(I,J).EQ.'3') GOTO 48 ! 45 CONTINUE ! 48 MF(J)=I ! 42 CONTINUE !-------------------------------------------------------------- ! New codes !--------------------------------------------------------------- DO J=1,JMXB1 IE(J)=IMXB1 DO II=1,IMXB1 I=IMXB1-II+1 IF (KSKP(I,J).EQ.'1') THEN IE(J)=I EXIT ELSE IF (II == IMXB1) THEN IE(J) = I+1 EXIT END IF END IF END DO END DO IE(JMXB)=IE(1) DO J=2,JMXB MF(J)=IMXB DO II=1,IMXB1 I=IMXB-II IF (ITREE(I,J).EQ.'1'.OR.ITREE(I,J).EQ.'3') THEN MF(J)=I EXIT END IF END DO END DO ! ! ! IF (DOLLAR.EQ.'$') THEN MF(1)=MF(JMXB) ELSE MF(1)=MF(2) ENDIF ! ! !------------------------------------------------------------ ! rewrite the following codes to remove the computed GOTO ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! old codes !--------------------------------------------------------------- C IF NO BANKS, THEN SET BANK HEIGHTS TO BARRIER HEIGHT ! IF (NSQRWC.EQ.0) GO TO 250 ! DO 200 L=1,NSQRWC ! IF (BANK(L,1).NE.0.) GO TO 190 ! I=ISQR(L) ! J=JSQR(L) ! ISS=1 ! 180 BANK(L,ISS)=AMIN1(ZBM(I,J),ZBM(I,J+1),ZBM(I+1,J),ZBM(I+1,J+1)) ! GO TO (190,200),ISS ! 190 IF (BANK(L,2).NE.0.) GO TO 200 ! K=ISIDE(L) ! I=ISQR(L)+IHH(K) ! J=JSQR(L)+JHH(K) ! ISS=2 ! GO TO 180 ! 200 CONTINUE !-------------------------------------------------------------- ! New codes !--------------------------------------------------------------- ! IF NO BANKS, THEN SET BANK HEIGHTS TO BARRIER HEIGHT ! IF (NSQRWC.EQ.0) THEN ! GO TO 250 ! END IF IF (NSQRWC /= 0) THEN ! 250 DO L=1,NSQRWC IF (BANK(L,1).NE.0.) THEN IF (BANK(L,2).EQ.0.) THEN K=ISIDE(L) I=ISQR(L)+IHH(K) J=JSQR(L)+JHH(K) ISS=2 BANK(L,ISS)=AMIN1(ZBM(I,J),ZBM(I,J+1), $ ZBM(I+1,J),ZBM(I+1,J+1)) END IF ELSE I=ISQR(L) J=JSQR(L) ISS=1 BANK(L,ISS)=AMIN1(ZBM(I,J),ZBM(I,J+1), $ ZBM(I+1,J),ZBM(I+1,J+1)) IF (ISS == 1) THEN IF (BANK(L,2).EQ.0.) THEN K=ISIDE(L) I=ISQR(L)+IHH(K) J=JSQR(L)+JHH(K) ISS=2 BANK(L,ISS)=AMIN1(ZBM(I,J),ZBM(I,J+1), $ ZBM(I+1,J),ZBM(I+1,J+1)) END IF END IF END IF END DO ! ! C C SET UP DEPTHS AND SIDE HEIGHTS FOR 1-DIM FLOW AND RIVER BANKS C SET UP SPECIAL RAISED WEIRS OR SILLS FOR 1-DIM FLOW C AND ZBMIN FOR FLOW CUTS IN SUBROUTINE CRDRD2. ! ! ! IF (NSQRW.EQ.NSQRS) GOTO 2333 IF (NSQRW /= NSQRS) THEN ! N1=NSQRS+1 ! DO 232 L=N1,NSQRW DO L=N1,NSQRW ! I=ISQR(L) J=JSQR(L) K=ISIDE(L) I1=I+IZ1(K) I2=I+IZ2(K) J1=J+JZ1(K) J2=J+JZ2(K) IF (HWEIR(L).EQ.0.) HWEIR(L)=AMIN1(ZBM(I1,J1),ZBM(I2,J2)) ! 232 CONTINUE END DO ! ! DO 1232 N=N1,NSQRW DO N=N1,NSQRW ! I=ISQR(N) J=JSQR(N) K=ISIDE(N) I1=I+IZ1(K) I2=I+IZ2(K) J1=J+JZ1(K) J2=J+JZ2(K) ITREE(I1,J1)='Q' ITREE(I2,J2)='Q' ! 1232 CONTINUE END DO ! END IF ! 2333 CONTINUE ! DO 3222 J=2,JMXB1 ! DO 3222 I=2,IMXB1 DO J=2,JMXB1 DO I=2,IMXB1 IF (ITREE(I,J).EQ.'Q') THEN C COMMENTED BY NSM 11/20/2010 : Should it be HTER or HTER + 1? C ZBM(I,J)=36. ZBM(I,J)=HTER+1 ITREE(I,J)='4' ENDIF END DO END DO ! 3222 CONTINUE ! DO 240 L=1,NSQRWC DO L=1,NSQRWC I=ISQR(L) J=JSQR(L) ZINT=-ZB(I,J) K=ISIDE(L) II=I+IHH(K) JJ=J+JHH(K) ZEXT=-ZB(II,JJ) C FOR WEIR HEIGHTS NOT A SILL ! ! IF (HWEIR(L).NE.0.) GO TO 230 ! 220 HWEIR(L)=AMAX1(ZINT,ZEXT) ! IF (HWEIR(L) == 0.) THEN HWEIR(L)=AMAX1(ZINT,ZEXT) END IF ! ! I1=I+IZ1(K) ! 230 I1=I+IZ1(K) ! ! I2=I+IZ2(K) J1=J+JZ1(K) J2=J+JZ2(K) IZBX=AMIN1(ZBM(I1,J1),ZBM(I2,J2)) ZBMIN(L)=IZBX END DO ! 240 CONTINUE C ! IF (NSQRW.EQ.NSQRS) GOTO 250 ! DO 2221 N=N1,NSQRW ! I=ISQR(N) ! J=JSQR(N) ! K=ISIDE(N) ! I1=I+IZ1(K) ! I2=I+IZ2(K) ! J1=J+JZ1(K) ! J2=J+JZ2(K) !C WRITE(*,2222) I,J,K,HWEIR(N),ZBMIN(N),ZBM(I1,J1),ZBM(I2,J2) !C 2222 FORMAT(3I5,4F6.1) IF (NSQRW /= NSQRS) THEN DO N=N1,NSQRW I=ISQR(N) J=JSQR(N) K=ISIDE(N) I1=I+IZ1(K) I2=I+IZ2(K) J1=J+JZ1(K) J2=J+JZ2(K) END DO END IF ! END IF !IF (NSQRWC /= 0) THEN ! ! ! C 250 CONTINUE C DEFINE GRID COORDINATES FORSTATIC BOUNDARY CONDITION. N=0 ! DO 1200 I=1,IMXB1 ! DO 1200 J=1,JMXB1 ! IF (KSKP(I,J).NE.'2') GOTO 1200 ! N=N+1 ! IF (N.GT.NBCPTS) THEN ! WRITE (*,*) 'Number of boundary points exceeded', N, NBCPTS ! STOP ! ENDIF ! ISH(N)=I ! JSH(N)=J ! KSKP(I,J)='0' ! 1200 CONTINUE ! ! DO I=1,IMXB1 DO J=1,JMXB1 IF (KSKP(I,J) == '2') THEN N=N+1 IF (N.GT.NBCPTS) THEN WRITE (*,*) 'Number of boundary points exceeded', N, NBCPTS STOP END IF ISH(N)=I JSH(N)=J KSKP(I,J)='0' END IF END DO END DO ! ! ! NBCPT=N C WRITE(*,225) NBCPT C 225 FORMAT(I5,' TOTAL POINTS FOR STATIC BOUNDARY CONDITION') C NNN=(NBCPT-1)/15+1 C DO 223 K=1,NNN C N1=1+(K-1)*15 C N2=MIN(N1+14,NBCPT) C WRITE(*,222) K,(ISH(N),N=N1,N2) C WRITE(*,224) K,(JSH(N),N=N1,N2) C 222 FORMAT(T62,'I-B.C.',I2,T1,15I4) C 224 FORMAT(T62,'J-B.C.',I2,T1,15I4) C 223 CONTINUE c noflooding option cd Debug Huiqing Liu/MDL cd if (nofld.eq.'+') then IF (NOFLD == '+'.OR.IFLD == 0) THEN DO J=2,JMXB1 DO I=2,IMXB1 IF (ZBM(I,J) > 0.) THEN ITREE(I,J) = '4' ZBM(I,J) = HTER+1 END IF END DO END DO END IF C RETURN END SUBROUTINE CONTTY C JELESNIANSKI SEPTEMBER 198TDL IBM 360/195 C C PURPOSE C THIS SUBROUTINE (CONTTY) IS AN ALGORITHM FOR THE CONTINUITY C EQUATION, FOR AN IMAGE PLANE; POLAR COORDINATES ARE C TRANSFORMED ONTO THE IMAGE PLANE. THE SYSTEM AXIS ARE C (P,Q), BUT THE SUBROUTINE USES (I,J) FOR INTEGER SBSCRPTNG C TWO LEVELS IN TIME (FORWARD IN CONTTY PLUS BACKWARD IN MOMN) C ARE USED. THE GRID (SEE DIAGRAM BELOW) IS STAGGERED IN SPACE C A 5-POINT SCHEME (4 CORNERS OF A SQUARE) ARE USED FOR C THE MOMENTUM SPACE GRADIENT. C C .****.****. WHEN WATER DESCENDS BELOW SQUARE C *I-1 * I * (I,J), THE MOMENTUM ON 4-SURROUNDNG C * + * + * CORNERS ARE REDUCED TO ALLOW THE C * J * J * SQUARE TO BECOME DRY. THE C .****.****.****. PREVIOUSLY COMPUTED SURGES ON THE C *I-1 * I *I+1 * 4-SURROUNDING SQUARES ARE UPDATED C * + * + * + * ACCORDING TO THE REVISED MOMENTUMS C *J-1 *J-1 *J-1 * ON THE 4-CORNER POINTS OF (I,J). C .****.****.****. SEE DIAGRAM BELOW FOR SUBSCRIPTING. C C DATA SET USE C C VARIABLES C IMXB1 JMXB1 = MAX SUBSCRIPTS FOR SURGE SQUARES C HB( , ) = SURGE HEIGHTS C HSUB( , ) = STORAGE FOR A SECOND COPY OF HB C IS( ) = START LOOP, EXCLUDE HI TERRAIN, SET IN 'CRDRD2' C IB IFN = BEGINNING/ENDING A 'DO LOOP' IN I-DIRECTION C ISKIP = 1,FOR J=1,CHECK ZBM AT 2 BNDRY CORNERS FOR SKIPPIN C STAATIC CONDITION. C 2,AND I=IMXB1 FOR BOTTOM BOUNDARY C 3,FOR J=JMXB1 FOR RIGHT BOUNDARY C JUMP = CONDITIONAL INTEGER FOR STANDARD LOOP, C OR RETRIEVAL OF PREVIOUSLY COMPUTED SURGES, C SEE DIAGRAM ABOVE FOR AFFECTED SQUARES C LEAP = CONDITIONAL 'GO TO' TO DRY AN (I,J) SQUARE, C SEE DIAGRAM ABOVE FOR AN (I,J) SQUARE C IFST ISND = BEGIN/END I SUBSCRIPT ON UPDATED (J-1) SQUARES C IPASS = CONDITIONAL INTEGER, TO SET A SQUARE DRY C ZB( , ) = DEPTHS AT CENTER OF SURGE SQUARES C IP(4) JP(4) = SHIFT SUBSCRIPTS TO 4 MOMENTUM POINTS C UB VB( , ) = COMPONENTS OF MOMENTUM FIELDS C HPST = DUMMY VALUE FOR SUEGE VALUEC C DUX DVY = COMPONENTS OF MOMENTUM GRADIENT C FCT = FACTOR TO MULTIPLY MOMENTUM GRADIENT C CHNG = SURGE PLUS DEPTHS (TO TEST FOR DRY SQUARE) C DENM = CHANGE IN SURGE, ONE TIME ITERATION C ANUM = CHANGE IN SURGE HT, TO TO DRY A SQUARE C RFCT = RATIO, <1, OF ANUM TO DENM, REDUCTION FACTOR C ISET = I-SUBSCRIPT FOR DRIED SQUARE C C GENERAL COMMENTS C THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS CALLED C IN BY SUBROUTINE 'CMPUTE' C THE SYSTEM ROUTINE 'ERRSET' IS USED TO ISOLATE VARIABLES C IF AN UNDER/OVERFLOW OCCURS; VARIABLE PRINT OUT IS C AIDED VIA 'HELP' AND THE 'COMMON/CHELP/'. C C THE SUBSCRIPT FOR COMPUTATIONS OF SURGES ARE OFFSET (.5,.5) C FROM THE MOMENTUM POINTS. SEE BELOW: C C (I,J) (I,J+1) C . . C .=MOMENTUM POINTS (UB,VB) C +=SURGE POINTS (HB) C (I,J) C + C .=BARRIER POINTS (ZBM) C +=DEPTH POINTS (ZB) C C . . IP(4)=0,1,1,0 JP(4)=0,0,1,1 C (I+1,J) (I+1,J+1) C INCLUDE 'parm.for' C COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /DUMB8/ IP(4),JP(4),IH(4),JH(4) COMMON /DUM88/ IIH(4),JJH(4),IHH(4),JHH(4) COMMON /DUMB4/ X1B(10) COMMON /CHELP/ ANUM,DENM,RFCT,I,J,HBJ,HPST,UP(4),VP(4) ! added by Huiqing.Liu /MDL hold for previous time step HB(I,J) COMMON /DUMBOLD/ HB_OLD(M_,N_) LOGICAL :: Jump_n C C PRODUCE A DUPLICATE COPY OF HB(,) IN CASE WATER DECEND BELOW C GROUND. !------------------------------------------------------------ ! rewrite the following codes to remove the computed GOTO ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! old codes !--------------------------------------------------------------- ! DO 100 J=1,JMXB ! DO 100 I=1,IMXB1 ! 100 HSUB(I,J)=HB(I,J) !C SET UB,VB TO LOWER BOUNDS OR ZEROS TO AVOID UNDERFLOW. ! DO 110 J=2,JMXB1 ! DO 110 I=2,IMXB1 ! IF(ABS(UB(I,J)).GT.1.E-10) GO TO 110 ! UB(I,J)=0. ! VB(I,J)=0. ! 110 CONTINUE !-------------------------------------------------------------- ! New codes !--------------------------------------------------------------- DO J=1,JMXB DO I=1,IMXB1 HB_OLD(I,J)=HB(I,J)!store the previous time step HB(I,J), which will be used in Tide nudging HSUB(I,J)=HB(I,J) END DO END DO C SET UB,VB TO LOWER BOUNDS OR ZEROS TO AVOID UNDERFLOW. DO J=2,JMXB1 DO I=2,IMXB1 IF(ABS(UB(I,J)) <= 1.E-10) THEN UB(I,J)=0. VB(I,J)=0. END IF END DO END DO ! ! ! C C SET PARAMETERS FOR J=1 C J=1 IB=IS(J) IFN=IE(J) JUMP=4 LEAP=1 GO TO 130 C C BEGINING LOOP WHEN J DECENDS BY ONE FOR DRYING OF A SQUARE C ! DO ! 120 IB=IFST ! IFN=ISND C C BEGINING OF REGULAR J LOOP FROM 1 THRU JMXB1 C 130 CONTINUE ! DO IF (IFN.LT.IB) GOTO 330 ! IF (IFN.LT.IB) EXIT ! Jump_n = .false. DO 320 I=IB,IFN Jump_n = .false. C IF (KSKP(I,J).EQ.'0') CYCLE ! C LEAP=1, REGILAR CASE. LEAP=2, RETURN TO THE SQUARE BEING DRIED. C !------------------------------------------------------------ ! rewrite the following codes to remove the computed GOTO ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! New codes !--------------------------------------------------------------- ! IF (LEAP == 2) THEN IF (I.EQ.IPASS) THEN HB(I,J)=-ZB(I,J) LEAP=1 CYCLE END IF END IF ! ! ! C C BYPASS CHECKINGS FOR DEPTH >30 FT, OR TERRAIN > 35 FT. C DUX=-UB(I,J)+UB(I+1,J)+UB(I+1,J+1)-UB(I,J+1) C IF (ZB(I,J) <= 30.) THEN IF (DUX == 0.) THEN CYCLE END IF END IF ! CONTINUITY EQUATION. HPST=HSUB(I,J) ! ! Z=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J) FCT=X1B(5)/Z DVY=-VB(I,J)-VB(I+1,J)+VB(I+1,J+1)+VB(I,J+1) HBJ=HPST-FCT*(DUX+DVY) C C TEST IF SURGE IS BELOW DRY LAND, IF SO, SET SURGE TO LAND HGT. C REVISE SURROUNDING TRANSPORTS TO GIVE ZERO TOTAL DEPTH. C CHNG=HBJ+ZB(I,J) ! ! C**********END OF LOOP FOR REGULAR RETURN******************************* C C*********************************************************************** C C READJUST TRANSPORTS TO ACCOUNT SUQARES DECENDING BELOW TERRAIN. C ! ! IF (CHNG >= 0.) THEN HB(I,J) = HBJ CYCLE ELSE !CCCCCC MODIFICATION 10/27/99 RE-TESTED by Pro Fortran 09/20/01 IF (-CHNG < 1E-3) THEN HB(I,J) = -ZB(I,J) CYCLE ENDIF END IF ! ! 230 CONTINUE ! ! ! ! IF (HSUB(I,J).EQ.-ZB(I,J)) THEN CYCLE END IF DENM=HBJ-HPST DENM=SIGN(AMAX1(1.E-5,ABS(DENM)),DENM) ANUM=-ZB(I,J)-HPST RFCT=ANUM/DENM IF (RFCT >= 0.) THEN IF (ABS(RFCT).LT.1.) THEN IF(ABS(RFCT).LT.1.E-5) RFCT=0. END IF ELSE IF(ABS(RFCT).LT.1.E-5) RFCT=0. END IF ! ! DO 300 K=1,4 IA=I+IP(K) JA=J+JP(K) UB(IA,JA)=RFCT*UB(IA,JA) VB(IA,JA)=RFCT*VB(IA,JA) 300 CONTINUE C RECOMPUTE 3 SQUARES ON J-1 LINE J=J-1 JUMP=1 ISET=I LEAP=1 ! ! IF (J.EQ.0) GO TO 330 ! IF (KSKP(I,J).EQ.'0') GOTO 330 ! 310 IF (IS(J).GT.I) GO TO 330 IF (J.EQ.0) THEN Jump_n = .true. EXIT END IF IF (KSKP(I,J).EQ.'0') THEN Jump_n = .true. EXIT ENDIF IF (IS(J).GT.I) THEN Jump_n = .true. EXIT ENDIF ! IF (J.EQ.0.OR.KSKP(I,J).EQ.'0'.OR.IS(J).GT.I) THEN ! Jump_n = .true. ! EXIT ! END IF ! ELSE ! ! IFST=MAX0(I-1,IS(J)) ISND=MIN0(I+1,IE(J)) ! GO TO 120 IB=IFST IFN=ISND GO TO 130 ! ! 320 CONTINUE 330 CONTINUE ! IF (Jump_n) EXIT ! END DO !DO WHILE end Huiqing.Liu/MDL C INCREMENT J, OR CONTINUE ON OLD J LINE JUMP=JUMP+1 LEAP=1 C TEST FOR INCREMENTED OR OLD J ! ! IF (JUMP.NE.2) GO TO 340 C OLD J LINE, WITH 2 OLD SQUARES TO BE UPDATED ! J=J+1 ! LEAP=2 ! IB=MAX0(ISET-1,IS(J)) ! IFN=IE(J) ! IPASS=ISET ! GO TO 130 c INCREMENT J FOR NEW COMPUTATION LINE ! 340 J=J+1 ! IF (J.GT.JMXB1) GO TO 350 ! IB=IS(J) ! IFN=IE(J) ! GO TO 130 C C*********** RETURN WITH A NEW J-LINE ********************** C ! IF (JUMP.NE.2) THEN J=J+1 ! IF (J.GT.JMXB1) GO TO 350 IF (J.LE.JMXB1) THEN IB=IS(J) IFN=IE(J) GO TO 130 END IF ELSE C OLD J LINE, WITH 2 OLD SQUARES TO BE UPDATED J=J+1 LEAP=2 IB=MAX0(ISET-1,IS(J)) IFN=IE(J) IPASS=ISET GO TO 130 c INCREMENT J FOR NEW COMPUTATION LINE END IF ! 350 CONTINUE DO 360 I=1,IMXB1 360 HB(I,JMXB)=HB(I,1) C RETURN END SUBROUTINE BDRYHT C JELESNIANSKI SEPTEMBER 1980 TDL IBM 360/195 C C PURPOSE C THIS SUBROUTINE (BDRYHT) COMPUTES STATIC HEIGHTS ON DEEP C WATER BOUNDARIES ONLY, SEE DIAGRAM BELOW: C C ITREE( , ) C = C B.9 . C O C U + C N .*********. C D.8 . * * C | A* * * C / R* + * * C / Y* * + * COSL,SINL C / .6 . * * C DEEP / / * * C WATERS / / * * C / / / .*********. C / / / | C / .6 .----. . C / * C / * + + C / * C V .****.----.****.BOUNDARY----> C 6 6 8 9 4= ITREE(I,J) C DEEP C <----WATERS----> C C C DEEP WATER CAN BE ON ONE, TWO OR NEITHER OF THE TWO SIDE C BOUNDARIES, BUT IT WILL EXIST AT LEAST ON A SEGMENT OF C THE CIRCLE WITH LARGEST RADIUS OF THE POLAR GRID. C C DATA SET USE C NONE C C VARIABLES C AX AY = COMPS OF TOTAL STRM MOTION, ADVANCED IN 'STMVAL' C C1 C2 = INITIAL COMPS OF STORM, SET IN 'INTVAL' C COSL( ) SINL( ) = CO-SINE OF ANGLE, RAYS TO X-AXIS, (HEIGHT POINTS C IMXB1 = MAX I-SUBSCRIPT FOR HEIGHT POINTS C SEADTM = INITIAL HEIGHT OF THE SEA (NO STATIC HEIGHTS) C DELP(800) = STATIC HEIGHTS AT MILE INTERVALS FROM STORM CENTER C HB( , ) = SURGE HEIGHTS C ITREE( , ) = A, ON AT LEAST ONE BOUNDARY CORNER AS STATIC C HEIGHT IS USED C C GENERAL COMMENTS C THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS CALLED C IN BY SUBROUTINE 'CMPUTE'. C THIS SUBROUTINE IS CALLED IN BY SUBROUTINE 'CONTTY' C WHICH RESIDES IN OVERLAY 'CMPUTE' C INCLUDE 'parm.for' C PARAMETER (NBCPTS=12000) COMMON /BCPTS/ TIDESH(NBCPTS),BOUN_NEST(NBCPTS), ! Added by Huiqing Liu/MDL Nov. 2013 * NBCPT,ISH(NBCPTS),JSH(NBCPTS) COMMON /DUMMY4/ S(800),C(800),P(800),DELP(800) COMMON /STRMSB/ C1,C2,C21,C22,AX,AY,PTENCY,RTENCY C STIME interferes with C code, so switched to STIME2 COMMON /STIME2/ ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /FRST/ X1(50),X12(50) COMMON /DATUM/ SEADTM,DTMLAK C C C BOUNDARY SQUARES DETERMINED IN SUBROUTINE DEPSFC C DO 100 N=1,NBCPT I=ISH(N) J=JSH(N) XR=ELPCL(I)*COSL(J) YR=ELPDL(I)*SINL(J) X=XR-C1-AX Y=YR-C2-AY RSQ=X*X+Y*Y R1=SQRT(RSQ)/5280.+1. K=R1 R2=K DR=R1-R2 K=MIN0(K,790) HB(I,J)=DELP(K)+DR*(DELP(K+1)-DELP(K))+SEADTM+TIDESH(N) 100 CONTINUE C C STATIC HEIGHTS ON SIDE BOUNDARIES C RETURN END SUBROUTINE CHANNL C JELESNIANSKI SEPTEMBER 1980 TDL IBM 360/195 C C PURPOSE C THIS SUBROUTINE (CHANNL) COMPUTES TOTAL FLOW THROUGH A C PRISM, INVARIANT IN SPACE. THE PRISM BED IS HORIZONTAL. C THE HEAD BETWEEN THE TWO ENDS IS THE DRIVING FORCE TO PASS C WATER THRU THE CHANNEL. SEE GEOMETRY BELOW: C C XXXXXXXXXXXXXXXXXXXX C X' XX C X ' X X YL,Y0 C X ' X X ^ C X ' X X Y0,YL * C - XXXXXXXXXXXXXXXXXXXX X ^ * C . X ' X X * * C . X ' X X * * C C X ' X X * * C H X ' X X * * C D X ''''''''''''''X''''X XXXXXXXXXXXXXXXXXXXXX C P X ' X X <-----CHLNH---------> C H X ' X X H C . X ' X X T C . X' XX W C - XXXXXXXXXXXXXXXXXXXX H C C C <.....CHLNH........> C C THE EQUATION USED FOR THE CHANNEL FLOW IS; C C Q**2=(C1(YL**13/3-Y0**13/3))/(C2(YL**4/3-Y0**4/3)+CHLNH) C C DATA SET USE C NONE C C VARIABLES C NPSS = NUMBER OF PASSES C CHLNH(5) = LENGTH OF CHANNELS C G = GRAVITY C Y0 YL = SURGE AT ENDS OF THE CHANNEL C IPTO JPTO(2,5) = POSITIONS FOR 2-SQUARES AT ONE CHANNEL END C HB( , ) = SURGE HEIGHT C CHDPH(5) = CHANNEL DEPTHS C IPTL JPTL(2,5) = POSITIONS OF 2-SQUARES AT OTHER CHANNEL END C JPTL(2,5) = J-SUBSCRIPT FOR 2-SQUARES AT OTHER CHANNEL END C SIG = SIGN , + OR - C CDELY = CUT-OFF HEIGHT OF .1 FOOT FOR THE HEAD C AMAN-AGB13 = FIXED CONSTANTS C Q(5) = FLOW C DELT = UNIT TIME STEP C CHWTH(5) = WIDTH OF CHANNELS C EXTEXT(5) = TOTAL FLOW THROUGH CHANNEL C C GENERAL COMMENTS C THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS CALLED C IN BY SUBROUTINE 'CMPUTE'. C INCLUDE 'parm.for' C COMMON /FLWCPT/ NSQRS,NSQRW,NSQRWC,NPSS,NCUT COMMON /CHNLD/ CHLNH(5),CHWTH(5),CHDPH(5) COMMON /EGTH/ DELS,DELT,G,COR COMMON /CHNL/ IPT0(2,5),JPT0(2,5),IPTL(2,5),JPTL(2,5),ENTEXT(5) COMMON /CHNARA/ ARESQ0(5),ARESQL(5) COMMON /CHNLNO/ AMAN,BMAN,GB,AGB,AGB13 DIMENSION Q(5) DATA Q/5*0./ DATA CDELY/.1/ C C UPDATE HEIGHTS ON BOTH ENDS OF CHANNELS.(CONTINUITY) C DO 110 NPASS=1,NPSS QWDT=Q(NPASS)*DELT DO 100 K=1,2 I=IPT0(K,NPASS) J=JPT0(K,NPASS) II=IPTL(K,NPASS) JJ=JPTL(K,NPASS) HB(I,J)=HB(I,J)-.5*QWDT/ARESQ0(NPASS) HB(II,JJ)=HB(II,JJ)+.5*QWDT/ARESQL(NPASS) HB(I,J)=AMAX1(HB(I,J),-ZB(I,J)) HB(II,JJ)=AMAX1(HB(II,JJ),-ZB(II,JJ)) 100 CONTINUE 110 CONTINUE C C COMPUTE FLOW FROM HEAD BETWEEN THE CHANNEL ENDS. C DO 210 NPASS=1,NPSS CLNGTH=CHLNH(NPASS)*5280. GL=G*CLNGTH Y0=0. DO 120 K=1,2 I=IPT0(K,NPASS) J=JPT0(K,NPASS) 120 Y0=Y0+.5*HB(I,J) Y0=Y0+CHDPH(NPASS) YL=0. DO 130 K=1,2 II=IPTL(K,NPASS) JJ=JPTL(K,NPASS) 130 YL=YL+.5*HB(II,JJ) YL=YL+CHDPH(NPASS) SIG=1. ! IF(Y0.GT.YL) GO TO 140 IF (Y0.LE.YL) THEN Z=Y0 Y0=YL YL=Z SIG =-1. END IF 140 CONTINUE DELY=Y0-YL C C DO NOT CONSIDER HEAD LESS THAN .1 FT, IF(DELY.LT.CDELY) YL=Y0-CDELY C INITIAL GUESS FOR NEWTON'METHOD. TEST FOR SUPER-CRITICAL FLOW. HIT=YL Y=Y0 J=0 AF=GL+GB*(Y**(4./3.)) AF3=3.*AF AE=AMAN*(Y**(13./3.)) ! ! DO 150 HITO3=HIT**(1./3.) HIT2=HIT*HIT HIT3=HIT*HIT2 HIT10=HIT3*HITO3 HIT13=HIT*HIT10 CNUM=AF*HIT3+AGB*HIT13-AE CDNM=AF3*HIT2+AGB13*HIT10 HCRIT=HIT-CNUM/CDNM ! IF (ABS(HCRIT-HIT).LT.(.001)) GO TO 180 IF (ABS(HCRIT-HIT).LT.(.001)) EXIT J=J+1 ! IF (J.GT.10) GO TO 160 IF (J.GT.10) EXIT HIT=HCRIT ! GO TO 150 ! 160 CONTINUE END DO !!! ! ! CC WRITE(*,170)HCRIT,HIT 170 FORMAT (10X,'TOO MAY ITERATIONS. HCRIT=',F10.3,' HIT=',F10.3) 180 HIT=HCRIT ! ! ! C C DO NOT PERMIT SUPER-CRITICAL FLOW. ! IF(YL.GT.HCRIT) GO TO 190 IF (YL.GT.HCRIT) THEN Y0TRD=Y0**(1./3.) YLTRD=YL**(1./3.) Y04=Y0*Y0TRD YL4=YL*YLTRD Z=Y0*Y0 Y013=Z*Z*Y0TRD Z=YL*YL YL13=Z*Z*YLTRD QFLWSQ=AMAN*(Y013-YL13)/(CLNGTH+BMAN*(Y04-YL4)) QFLW=SQRT(ABS(QFLWSQ)) ELSE QFLW= SQRT(G*HCRIT)*HCRIT END IF ! QFLW= SQRT(G*HCRIT)*HCRIT ! GO TO 200 ! 190 Y0TRD=Y0**(1./3.) ! YLTRD=YL**(1./3.) ! Y04=Y0*Y0TRD ! YL4=YL*YLTRD ! Z=Y0*Y0 ! Y013=Z*Z*Y0TRD ! Z=YL*YL ! YL13=Z*Z*YLTRD ! QFLWSQ=AMAN*(Y013-YL13)/(CLNGTH+BMAN*(Y04-YL4)) ! QFLW=SQRT(ABS(QFLWSQ)) 200 CONTINUE Q(NPASS)=QFLW*SIG Q(NPASS)=Q(NPASS)*CHWTH(NPASS) C FOR HEAD < .1 FT,INTERPOLATE LINEARLY BETWEEN .1 AND 0. IF(DELY.LT.CDELY) Q(NPASS)=Q(NPASS)*DELY/CDELY 210 CONTINUE RETURN END SUBROUTINE FLW1DM C JELESNIANSKI SEPTEMBER 1980 TDL IBM 360/195 C C PURPOSE C THIS SUBROUTINE (FLW1DM) COMPUTES 1-DIM FLOW ACROSS CHANNEL C SIDES OF INTERLOCKING SQUARES. WATER PASSES ONLY THRU 'OPEN' C SIDES OF A 1-DIM FLOW CHANNEL, SEE DIAGRAM BELOW: C C .****.****.****.****. -FLOW *NO FLOW C - - - - * -THRU *THRU C - + - + - + - + * -CHANNEL *CHANNEL C - - - - * -SIDE *SIDE C -****.****.****.----. C * * (OPEN) (CLOSED) C * + * C * * C .----. C THE 'LOOP' OF THE CONTINUITY EQUATION FILLS/DEPLETES C TWO ADJOINING SQUARES WITH FLOW THRU AN 'OPEN' SIDE, C ONE SIDE AT A TIME. IF WATER DROPS BELOW A SQUARE C FROM THE ACTION OF AN 'OPEN' SIDE, THEN THE FLOW THRU C THE SIDE IS REDUCED TO EXHAUST WATER ONLY TO A DRY C SQUARE. C C EACH 'OPEN' SIDE HAS AN INTERIOR AND AN EXTERIOR SQUARE. C SUBROUTINE 'CRDRD2' READS THE SUBSCRIPTS FOR INTERIOR C SQUARES AND SETS UP PROCEDURES TO ISOLATE THE C ADJOINING EXTERIOR SQUARES AND 'OPEN' SIDES VIA C 'COMMON/DUMB18/'. THE PRESENT SUBROUTINE ISOLATES THE TWO C INTERIOR/EXTERIOR SQUARES AND THE ADJOINING SIDE FOR C 1-DIM FLOW. A 1-DIM FLOW SQUARE CAN BE EITHER AN INT/EXT C SQUARE, OR BOTH, 1-4 TIMES, DEPENDING ON THE AMOUNT OF C 'OPEN' SIDES BOUNDING A SQUARE. C C I,J SHIFTS ABOUT AN (I,J) INTERIOR SQUARE FOR 4 POSSIBLE C EXTERIOR SQUARES ARE: C C .////. C IHH=0/ C / 4+ / C JHH=1/ C .////.////.////. C IHH=-1/ I /IHH=1 C / 1+ / + / 2+ / C JHH=0 / J /JHH=0 C .////.////.////. C IHH=0/ C / 3+ / C JHH=-1 C .////. C C DATA SET USE C NONE C C VARIABLES C AI(800) = BOTTOM FRICTION, AT ONE FOOT INTERVALS C NSQRW = NO OF SQUARES WITH 1-DIM FLOW POSSIBLE C ISQR(L) JSQRL(L) = I/J-SUBSCRIPT FOR SQUARE WITH 1-DIM FLOW C HB( , ) = SURGE FIELD C ISIDE(300) = SIDE OF FLOW RELATIVE TO INTERIOR 1-DIM FLOW SQRE C IHH(4) JHH(4) = SHIFT I/J-SUBSCRIPT FOR CONTIGUOUS SQUARE C HINT HEXT = INTERIOR/EXTERIOR SURGE FOR SFC GRADIENT C HINT1 HEXT1 = INTERIOR/EXTERIOR SURGE HEIGHTS ON SQUARES C ZBMIN(300) = LOWEST BARRIER VALUE AT END OF AN 'OPEN' SIDE C ZL = DUMMY VALUE FOR ZBMIN(300) C ZB( , ) = DEPTH FIELD C HWEIR(300) = HIGHEST DEPTH OF INT/EXT SQUARES, OR WEIR HEIGHT C HOVT = DUMMY VALUE FOR HWEIR(300) C FLWSQR(300 ) = TRANSPORT ACROSS AN 'OPEN' SIDE, TWO TIME LEVELS C GRIDR2( ) = JACOBEAN FOR POLAR COORDINATES C SIGA(4) = SIGN FOR FLOW ACROSS SIDES OF SQUARES C DPH1D(2) = DEPTH+SURGE AT INT/EXT SQUARES, PRESENT TIME C CHNGHI CHNGHE = SURGE INCREASE ON INT/EXT SQRES, FROM 'OPEN' SIDE C DEPLI DEPLE = TOTAL DEPTH OF INT/EXT SQRES, FUTURE TIME C CII = ACCUMULATED WATER FROM 'OPEN' SIDES, INTERIOR SQR C SUBE = SURGE FROM EXT SQR TO PREVENT WATER BLW INT SQR C ) CEE = ACCUMULATED WATER FROM 'OPEN' SIDES, EXTERIOR SQR C SUBI = SURGE FROM INT SQR TO PREVENT WATER BLW EXT SQR C ID DPH = MEAN TOTAL DEPTH OF INT/EXT SQUARES, PRESENT TIME C DZL = HEIGHT OF HWEIR ABOVE MEAN DEPTHS C FXB FYB = COMPONENTS OF SURFACE STRESS C X1B(N) = FIXED COEFFICIENTS FOR FINITE DIFFERENCING C BR(300) BI(300) = FRICTION COEFFICIENTS AT 1-FOOT INTERVALS C JUMP = 2 FOR STRESS COMPUTATIONS, 1 FOR NO STRESS C C GENERAL COMMENTS C THIS IS A VERSION FOR 2-STEP SCHEME (JULY 1981). C THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS CALLED C IN BY SUBROUTINE 'CMPUTE'. SUBROUTINE 'FLW1DM' CALLS IN C C SOME (NOT ALL) ELEVTNS ON INT/EXT SQRES AT PRESENT TIME ARE: C C T=ZL=ZBMIN(300) C I C I C I C I C I***********>HINT=HB(I,JC) C HEXT<***********I / C / I /DPH1D(1)=HINT+ZB C / I / C DPH1D(2)/ I----------->-ZB(I,J) C / I / C / I / C -ZB<----------I / C / I / C / I / C / I / C / I / C DEPTH=ZB=V I V=ZB=DEPTH C <----------------I----------------->NGVD DATUM C (EXTERIOR) (INTERIOR) C INCLUDE 'parm.for' C parameter (ndp=600) COMMON /EGTH/ DELS,DELT,G,COR COMMON /SCND/ AR(ndp),AI(ndp),BR(ndp),BI(ndp),CR(ndp),CI(ndp) COMMON /DUMB4/ X1B(10) COMMON /DUMB8/ IP(4),JP(4),IH(4),JH(4) COMMON /DUM88/ IIH(4),JJH(4),IHH(4),JHH(4) COMMON /MF/ DPH,HPD(4),ID COMMON /FLWCPT/ NSQRS,NSQRW,NSQRWC,NPSS,NCUT C DIMENSION SIGA(4),DPH1D(2) ! Huiqing.Liu /MDL Feb. 2015 add one more logical variable Skip_h to ! help remvoe computed GO TO statments LOGICAL :: Skip_h C DATA SIGA/-1.,1.,-1.,1./ ! Initialize Skip_h to .false. Skip_h = .false. C DO 100 L=1,NSQRWC I=ISQR(L) J=JSQR(L) K=ISIDE(L) II=I+IHH(K) JJ=J+JHH(K) HSUB(I,J)=0. HSUB(II,JJ)=0. 100 F1DACT(L)='F' C CXXXXXXXXX BEGIN 1ST 'LOOP'. TREAT 2-DIM FLOW (FROM 'CONTTY') ONTO CXXXXXXXXX A 1-DIM SQUARE AS A SOURCE. C DO 510 L=1,NSQRWC C ! IF (FLWSQR(L).EQ.0.) CYCLE ! ! C ISOLATE A PARTICULAR SQUARE AND ITS SURFACE HEIGHT I=ISQR(L) J=JSQR(L) C C DETERMINE WHICH OF 4 SIDES HAS FLOW AND ISOLATE ADJACENT SQUARE C AND ITS SURFACE HEIGHT K=ISIDE(L) II=I+IHH(K) JJ=J+JHH(K) C C START CONTINUITY EQUATION FOR ONE SIDE OF A SQUARE C COMPUTE TOTAL FLOW IN/OUT OF INTERIOR SQUARE 110 Z=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J) Z1=ELPDL2(II)+(ELPCL2(II)-ELPDL2(II))*SINL2(JJ) CHNGHI=-SIGA(K)*X1B(5)*FLWSQR(L)/Z CHNGHE=-CHNGHI*Z/Z1 C C DETERMINE 1-DIM, OR EXPANSION, FLOW KCUT=1 IF (L.GT.NSQRW) KCUT=2 !------------------------------------------------------------ ! rewrite the following codes to remove the computed GOTO ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! Old codes !--------------------------------------------------------------- ! GO TO (120,130), KCUT C C CHECK FOR EXISTENCE OF BANKS ! 120 IF (DELCUT(L).GE.1.) GO TO 470 ! GO TO 140 ! 130 LL=L-NSQRW ! IF (CUTLI(LL).GE.1..AND.CUTLE(LL).GE.1.) GO TO 420 !C !C TWO CASES EXIST: OUTFLOW VERSUS INFLOW !C CASE 1: CHNGHI NEG (OUTFLOW), CHNGHE POS (INFLOW), LEAP=1 !C CASE 2: CHNBHI POS (INFLOW), CHNGHE NEG (OUTFLOW), LEAP=2 ! 140 LEAP=1 ! IF (CHNGHI.GT.0.) LEAP=2 ! GO TO (192,194),LEAP ! 192 CHNG=CHNGHI ! CSTRT=CHNGHI ! HT=HB(I,J) ! BK=BANK(L,1) ! IF(KCUT.EQ.2) CUTW=CUTLI(LL) ! GO TO (200,300),LEAP ! 194 CHNG=CHNGHE ! CSTRT=CHNGHE ! HT=HB(II,JJ) ! BK=BANK(L,2) ! IF(KCUT.EQ.2) CUTW=CUTLE(LL) ! GO TO (300,200),LEAP !C******LOGIC FOR NEGATIVE WATER CHANGES (OUTFLOW) ! 200 IF (HT.GT.BK) GO TO 230 ! GO TO (210,220),KCUT !C INITIAL WATER BELOW BANK AND DESCENDING ! 210 CHNG=CHNG ! GO TO 280 ! 220 CHNG=CHNG*CUTL(LL)/CUTW ! GO TO 272 !C INITIAL WATER ABOVE BANK AND DESCENDING ! 230 DELB=HT-BK ! GO TO (240,260),KCUT !C WATER DOES NOT DESCEND BELOW BANK, PURE 1-D FLOW ! 240 IF(-CHNG*DELCUT(L).GT.DELB) GO TO 250 ! CHNG=CHNG*DELCUT(L) ! GO TO 272 !C WATER DESCENDS BELOW BANK, PURE 1-D FLOW ! 250 CHNG=CHNG-DELB+DELB/DELCUT(L) ! GO TO 272 !C WATER DOES NOT DESCEND BELOW BANK, CUT-TYPE FLOW ! 260 IF (-CHNG*CUTL(LL).GT.DELB) GO TO 270 ! CHNG=CHNG*CUTL(LL) ! GO TO 272 !C WATER DESCENDS BELOW BANK, CUT-TYPE FLOW ! 270 CHNG=(CHNG*CUTL(LL)+DELB)/CUTW-DELB ! 272 CONTINUE ! 280 GO TO (290,295),LEAP ! 290 CHNGHI=CHNG ! HSUB(I,J)=HSUB(I,J)+ABS(CSTRT-CHNG) ! GO TO 194 ! 295 CHNGHE=CHNG ! HSUB(II,JJ)=HSUB(II,JJ)+ABS(CSTRT-CHNG) ! GO TO 192 !C$$$$$$LOGIC FOR POSITIVE WATER CHANGES (INFLOW) ! 300 IF (HT.LT.BK) GO TO 330 ! GO TO (310,320),KCUT !C INITIAL WATER ABOVE BANK AND ASCENDING ! 310 CHNG=CHNG*DELCUT(L) ! GO TO 372 ! 320 CHNG=CHNG*CUTL(LL) ! GO TO 372 !C INITIAL WATER BELOW BANK AND ASCENDING ! 330 DELB=BK-HT ! GO TO (340,360),KCUT !C WATER DOES NOT ASCEND ABOVE BANK, PURE 1-D FLOW ! 340 IF (CHNG.GT.DELB) GO TO 350 ! CHNG=CHNG ! GO TO 380 C WATER ASCENDS ABOVE BANK, PURE 1-D FLOW ! 350 CHNG=DELB+(CHNG-DELB)*DELCUT(L) ! GO TO 372 C WATER DOES NOT ASCEND ABOVE BANK, CUT-TYPE FLOW ! 360 IF (CHNG*CUTL(LL).GT.DELB*CUTW) GO TO 370 ! CHNG=CHNG*CUTL(LL)/CUTW ! GO TO 372 C WATER ASCENDS ABOVE BANK, CUT-TYPE FLOW ! 370 CHNG=CHNG*CUTL(LL)-DELB*CUTW+DELB ! 372 CONTINUE ! 380 GO TO (390,395),LEAP ! 390 CHNGHE=CHNG ! HSUB(II,JJ)=HSUB(II,JJ)+ABS(CSTRT-CHNG) ! GO TO 470 ! 395 CHNGHI=CHNG ! HSUB(I,J)=HSUB(I,J)+ABS(CSTRT-CHNG) ! GO TO 470 ! 420 CHNGHI=CHNGHI*CUTL(LL)/CUTLI(LL) ! CHNGHE=CHNGHE*CUTL(LL)/CUTLE(LL) ! 470 HBJ=HB(I,J)+CHNGHI ! HBJJ=HB(II,JJ)+CHNGHE !-------------------------------------------------------------- ! New codes !--------------------------------------------------------------- IF (KCUT == 2) THEN LL=L-NSQRW IF (CUTLI(LL).GE.1..AND.CUTLE(LL).GE.1.) THEN ! GO TO 420 CHNGHI=CHNGHI*CUTL(LL)/CUTLI(LL) CHNGHE=CHNGHE*CUTL(LL)/CUTLE(LL) ! 470 HBJ=HB(I,J)+CHNGHI HBJJ=HB(II,JJ)+CHNGHE ! Set Skip_h .true. to let code skip the following IF BLOCK to meet the ! logic in the old codes Skip_h = .true. ELSE LEAP = 1 Skip_h = .false. END IF ELSE IF (DELCUT(L).GE.1.) THEN ! GO TO 470 HBJ=HB(I,J)+CHNGHI HBJJ=HB(II,JJ)+CHNGHE ! Set Skip_h .true. to let code skip the following IF BLOCK to meet the ! logic in the old codes Skip_h = .true. ELSE LEAP = 1 Skip_h = .false. END IF END IF ! ! IF (.not.Skip_h) THEN IF (CHNGHI.GT.0.) LEAP=2 ! IF (LEAP == 1) THEN IF (LEAP /= 2) THEN CHNG=CHNGHI !192 CSTRT=CHNGHI HT=HB(I,J) BK=BANK(L,1) IF(KCUT.EQ.2) CUTW=CUTLI(LL) ! 200 IF (HT.GT.BK) THEN ! 200 DELB=HT-BK !230 ! IF (KCUT == 1) THEN !230 IF (KCUT /= 2) THEN ! For safe consideration IF (-CHNG*DELCUT(L).GT.DELB) THEN !240 CHNG=CHNG-DELB+DELB/DELCUT(L) ELSE CHNG=CHNG*DELCUT(L) END IF ! 272 CHNGHI=CHNG !290 HSUB(I,J)=HSUB(I,J)+ABS(CSTRT-CHNG) CHNG=CHNGHE ! 194 CSTRT=CHNGHE HT=HB(II,JJ) BK=BANK(L,2) IF(KCUT.EQ.2) CUTW=CUTLE(LL) !300 IF (HT.LT.BK) THEN !The value of HT is changed DELB=BK-HT !330 IF (CHNG.GT.DELB) THEN !340 CHNG=DELB+(CHNG-DELB)*DELCUT(L) !350 ! 372/380 ELSE CHNG=CHNG ! 372/380 END IF ELSE CHNG=CHNG*DELCUT(L) !310 ! 372/380 END IF ELSEIF (KCUT == 2) THEN !230 IF (-CHNG*CUTL(LL).GT.DELB) THEN ! 260 CHNG=(CHNG*CUTL(LL)+DELB)/CUTW-DELB ELSE CHNG=CHNG*CUTL(LL) END IF !272 CHNGHI=CHNG !290 HSUB(I,J)=HSUB(I,J)+ABS(CSTRT-CHNG) CHNG=CHNGHE ! 194 CSTRT=CHNGHE HT=HB(II,JJ) BK=BANK(L,2) IF(KCUT.EQ.2) CUTW=CUTLE(LL) !300 IF (HT.LT.BK) THEN !300 DELB=BK-HT !330 IF (CHNG*CUTL(LL).GT.DELB*CUTW) THEN !360 CHNG=CHNG*CUTL(LL)-DELB*CUTW+DELB !370 !372 ELSE CHNG=CHNG*CUTL(LL)/CUTW !372 END IF ELSE CHNG=CHNG*CUTL(LL) ! 320 !372 END IF END IF ELSE !HT <= BK !200 ! IF (KCUT == 1) THEN !200 IF (KCUT /= 2) THEN !200 CHNG=CHNG ! 210 CHNGHI=CHNG ! 290 HSUB(I,J)=HSUB(I,J)+ABS(CSTRT-CHNG) CHNG=CHNGHE ! 194 CSTRT=CHNGHE HT=HB(II,JJ) BK=BANK(L,2) IF (KCUT.EQ.2) CUTW=CUTLE(LL) ! 300 IF (HT.LT.BK) THEN !300 DELB=BK-HT !330 IF (CHNG.GT.DELB) THEN !340 CHNG=DELB+(CHNG-DELB)*DELCUT(L) !350 ELSE CHNG=CHNG END IF ! 372 !380 ELSE CHNG=CHNG*DELCUT(L) !310 !372 END IF ELSEIF (KCUT == 2) THEN !200 CHNG=CHNG*CUTL(LL)/CUTW ! 220 ! 272 CHNGHI=CHNG ! 290 HSUB(I,J)=HSUB(I,J)+ABS(CSTRT-CHNG) CHNG=CHNGHE ! 194 CSTRT=CHNGHE HT=HB(II,JJ) BK=BANK(L,2) IF (KCUT.EQ.2) CUTW=CUTLE(LL) IF (HT.LT.BK) THEN !300 DELB=BK-HT !330 IF (CHNG*CUTL(LL).GT.DELB*CUTW) THEN !360 CHNG=CHNG*CUTL(LL)-DELB*CUTW+DELB !370 ELSE CHNG=CHNG*CUTL(LL)/CUTW END IF !372 ELSE CHNG=CHNG*CUTL(LL) !320 !372 END IF END IF END IF CHNGHE=CHNG !390 HSUB(II,JJ)=HSUB(II,JJ)+ABS(CSTRT-CHNG) ! GO TO 470 HBJ=HB(I,J)+CHNGHI HBJJ=HB(II,JJ)+CHNGHE ELSEIF (LEAP == 2) THEN CHNG=CHNGHE ! 194 CSTRT=CHNGHE HT=HB(II,JJ) BK=BANK(L,2) IF(KCUT.EQ.2) CUTW=CUTLE(LL) IF (HT.GT.BK) THEN ! 200 DELB=HT-BK !230 ! IF (KCUT == 1) THEN !230 IF (KCUT /= 2) THEN ! For safe consideration IF (-CHNG*DELCUT(L).GT.DELB) THEN !240 CHNG=CHNG-DELB+DELB/DELCUT(L) !250 ELSE CHNG=CHNG*DELCUT(L) END IF !272 CHNGHE=CHNG ! 295 HSUB(II,JJ)=HSUB(II,JJ)+ABS(CSTRT-CHNG) CHNG=CHNGHI ! 192 CSTRT=CHNGHI HT=HB(I,J) BK=BANK(L,1) IF(KCUT.EQ.2) CUTW=CUTLI(LL) ! 300 IF (HT.LT.BK) THEN !300 DELB=BK-HT !330 IF (CHNG.GT.DELB) THEN !340 CHNG=DELB+(CHNG-DELB)*DELCUT(L) !350 !372 ELSE CHNG=CHNG !380 END IF ELSE CHNG=CHNG*DELCUT(L) !310 END IF ! 372 ELSEIF (KCUT == 2) THEN !230 IF (-CHNG*CUTL(LL).GT.DELB) THEN ! 260 CHNG=(CHNG*CUTL(LL)+DELB)/CUTW-DELB !270 ELSE CHNG=CHNG*CUTL(LL) END IF !272 CHNGHE=CHNG ! 295 HSUB(II,JJ)=HSUB(II,JJ)+ABS(CSTRT-CHNG) CHNG=CHNGHI ! 192 CSTRT=CHNGHI HT=HB(I,J) BK=BANK(L,1) IF(KCUT.EQ.2) CUTW=CUTLI(LL) ! 300 IF (HT.LT.BK) THEN !300 DELB=BK-HT !330 IF (CHNG*CUTL(LL).GT.DELB*CUTW) THEN !360 CHNG=CHNG*CUTL(LL)-DELB*CUTW+DELB !370 !372 ELSE CHNG=CHNG*CUTL(LL)/CUTW !372 END IF ELSE CHNG=CHNG*CUTL(LL) !320 !372 END IF END IF ELSE !200 ! IF (KCUT == 1) THEN IF (KCUT /= 2) THEN CHNG=CHNG !210 !272 CHNGHE=CHNG !295 HSUB(II,JJ)=HSUB(II,JJ)+ABS(CSTRT-CHNG) CHNG=CHNGHI !192 CSTRT=CHNGHI HT=HB(I,J) BK=BANK(L,1) IF(KCUT.EQ.2) CUTW=CUTLI(LL) IF (HT.LT.BK) THEN ! 300 DELB=BK-HT !330 IF (CHNG.GT.DELB) THEN !340 CHNG=DELB+(CHNG-DELB)*DELCUT(L) !350 ELSE CHNG=CHNG END IF ELSE CHNG=CHNG*DELCUT(L) !310 END IF ELSEIF (KCUT == 2) THEN !200 CHNG=CHNG*CUTL(LL)/CUTW !220 !280 CHNGHE=CHNG !295 HSUB(II,JJ)=HSUB(II,JJ)+ABS(CSTRT-CHNG) CHNG=CHNGHI !192 CSTRT=CHNGHI HT=HB(I,J) BK=BANK(L,1) IF(KCUT.EQ.2) CUTW=CUTLI(LL) IF (HT.LT.BK) THEN ! 300 DELB=BK-HT !330 IF (CHNG*CUTL(LL).GT.DELB*CUTW) THEN !360 CHNG=CHNG*CUTL(LL)-DELB*CUTW+DELB !370 ELSE CHNG=CHNG*CUTL(LL)/CUTW END IF ELSE CHNG=CHNG*CUTL(LL) !320 END IF END IF END IF CHNGHI=CHNG ! 395 HSUB(I,J)=HSUB(I,J)+ABS(CSTRT-CHNG) ! GO TO 470 HBJ=HB(I,J)+CHNGHI HBJJ=HB(II,JJ)+CHNGHE END IF END IF ! Skip_h ! Reset Skip_h back to initial value .false. Skip_h = .false. !------------------------------------------------------------ C FINISH CONTINUITY EQUATION FOR ONE SIDE OF A SQUARE BY C TESTING IF WATER LIES BELOW SQUARES 480 DEPLI=HBJ+ZB(I,J) DEPLE=HBJJ+ZB(II,JJ) ! ! ! ! IF (DEPLI.GE.0.) THEN IF (DEPLE.GE.0.) THEN HB(I,J)=HBJ HB(II,JJ)=HBJJ C F1DACT(L)='T' ELSE CEE=CHNGHE-DEPLE SUBI=-CEE*Z1/Z HB(I,J)=AMAX1(HB(I,J)-CHNGHI+SUBI,-ZB(I,J)) HB(II,JJ)=-ZB(II,JJ) FLWSQR(L)=FLWSQR(L)*CEE/CHNGHE ! CYCLE END IF ELSE CII=CHNGHI-DEPLI SUBE=-CII*Z/Z1 HB(I,J)=-ZB(I,J) HB(II,JJ)=AMAX1(HB(II,JJ)-CHNGHE+SUBE,-ZB(II,JJ)) FLWSQR(L)=FLWSQR(L)*CII/CHNGHI ! CYCLE END IF 510 CONTINUE C*********************************************************************** C END OF HEIGHT UPDATES BY COMPUTED FLOW ACCROSS SIDES C*********************************************************************** C C USE UPDATED HEIGHTS COMPUTE FLOW ACCROSS SIDES. C C FRICTION CHANGE BETWEEN 1-2 FEET OF TOTAL DEPTH, USED C FOR EXTRAPOLATION FROM 1 TO 0 FT. C SET PARAMETERS FOR FRCPNT ARGUMENTS, LAKE WINDS, WITH STATIC C PRESSURE FORCE. SLPA=AI(2)-AI(1) NCATG=1 NPLS=1 C DO 630 L=1,NSQRWC C ISOLATE A PARTICULAR SQUARE AND ITS SURFACE HEIGHT I=ISQR(L) J=JSQR(L) HINT=HB(I,J) C DETERMINE WHICH OF 4 SIDES HAS FLOW AND ISOLATE ADJACENT SQUARE C AND ITS SURFACE HEIGHT K=ISIDE(L) II=I+IHH(K) JJ=J+JHH(K) HEXT=HB(II,JJ) C ISOLATE MINIMUM ZBM ON SIDE CONNECTING TWO SQUARES ZL=ZBMIN(L) C TEST IF SFC LIES ABOVE ZL, E.G., 2-DIM FLOW IN BOTH SQUARES ! ! ! IF (AMIN1(HINT,HEXT).GE.ZL) THEN FLWSQR(L)=0. HSUB(I,J)=0. HSUB(II,JJ)=0. CYCLE ELSE IF (HINT.EQ.-ZB(I,J).AND.HEXT.EQ.-ZB(II,JJ)) THEN FLWSQR(L)=0. HSUB(I,J)=0. HSUB(II,JJ)=0. CYCLE ELSE HOVT=HWEIR(L) IF (AMAX1(HINT,HEXT).LE.HOVT) THEN FLWSQR(L)=0. HSUB(I,J)=0. HSUB(II,JJ)=0. CYCLE END IF END IF END IF ! ! 530 CONTINUE C COMPUTE MOMENTUM ON ONE SIDE OF A SQUARE HINT1=HB(I,J) HEXT1=HB(II,JJ) HINT=AMIN1(AMAX1(HINT1,HOVT),ZL) HEXT=AMIN1(AMAX1(HEXT1,HOVT),ZL) JUMP=2 IF (AMAX1(HINT,HEXT).GE.ZL) JUMP=1 HEAD=HEXT-HINT HEAD=SIGA(K)*HEAD DPH1D(1)=HINT+ZB(I,J) DPH1D(2)=HEXT+ZB(II,JJ) DPH=.5*(DPH1D(1)+DPH1D(2)) CCCCCCCCCCC Arthur fix of dph close to 0. if (dph.lt.0.00000001) then FLWSQR(L)=0. HSUB(I,J)=0. HSUB(II,JJ)=0. CYCLE ! end if CCCCCCCCCCC Arthur end fix of dph close to 0. ID=DPH ID=AMIN1(ndp-1.,DPH) ID=MAX0(ID,1) BIOBR=BI(ID)/BR(ID) ! IF (NSQRW.NE.NSQRS) THEN IF (L.GT.NSQRS.AND.L.LE.NSQRW) BIOBR=0. END IF ! 1492 CONTINUE HCOFF=BR(ID)+BI(ID)*BIOBR FLWT=-X1B(1)*HCOFF*DPH*HEAD ! ! ! ! IF (JUMP.NE.1) THEN C C AVERAGE SHELTERING COEFF IF WATER BECOME LESS THAN 1 FT ABOVE C BED,OR 1 FT BELOW CHANNEL BANK,OR BARRIER HEIGHT. IF (KTREE(L).EQ.'T') THEN CSHLTR=(DPH/50.)**2 ELSE CSHINT=1.-DIM(1.,DIM(HINT1,HOVT)) CSHEXT=1.-DIM(1.,DIM(HEXT1,HOVT)) CSH2I=1.-DIM(1.,DIM(ZL,HINT1)) CSH2E=1.-DIM(1.,DIM(ZL,HEXT1)) CSHLTR=.5*AMIN1(CSHINT+CSHEXT,CSH2I+CSH2E) ENDIF C III=I JJJ=J !------------------------------------------------------------ ! rewrite the following codes to remove the computed GOTO ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- !-------------------------------------------------------------- ! New codes !--------------------------------------------------------------- IF (K == 1) THEN XR0=COSL(JJJ) YR0=SINL(JJJ) ELPCZ=ELPCT(III) ELPDZ=ELPDT(III) ELSEIF (K == 2) THEN III=I+1 XR0=COSL(JJJ) YR0=SINL(JJJ) ELPCZ=ELPCT(III) ELPDZ=ELPDT(III) ELSEIF (k == 3) THEN XR0=COST(JJJ) YR0=SINT(JJJ) ELPCZ=ELPCL(III) ELPDZ=ELPDL(III) ELSEIF (k == 4) THEN JJJ=J+1 XR0=COST(JJJ) YR0=SINT(JJJ) ELPCZ=ELPCL(III) ELPDZ=ELPDL(III) ENDIF ! ! 115 CONTINUE XR=XR0*ELPCZ YR=YR0*ELPDZ C ! ! Huiqing.Liu /MDL modified to ETSS which need read GFS wind instead of ! using parmametric winds June 2014 ! ! CALL FRCPNT(XR,YR,FXB,FYB,NCATG,CSHLTR,NPLS) !------------------------------------------------------------ ! rewrite the following codes to remove the computed GOTO ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! New codes !--------------------------------------------------------------- IF (k == 1) THEN III1=I JJJ1=J+1 ELSEIF (K == 2) THEN III=III+1 III1=III+1 JJJ1=JJJ+1 ELSEIF (K == 3) THEN III1=III+1 JJJ1=JJJ ELSEIF (K == 4) THEN JJJ=JJJ+1 III1=III+1 JJJ1=JJJ+1 ENDIF ! ! ! 1121 CONTINUE CALL FRCPNT1(III,JJJ,FXB,FYB,NCATG,CSHLTR,NPLS) CALL FRCPNT1(III1,JJJ1,FXB1,FYB1,NCATG,CSHLTR,NPLS) ! ! Huiqing.Liu /MDL ! C CHANGE FORCING FROM (X,Y) SYSTEM TO IMAGE PLANE XR=XR0*ELPDZ YR=YR0*ELPCZ FXBP=XR*FXB+YR*FYB FYBP=XR*FYB-YR*FXB C CORRECT STRESS DIRECTION ACCORDING TO ACTURAL FLOW IF(COS1D(L).NE.1.) THEN FXB=COS1D(L)*FXBP-SIN1D(L)*FYBP FYB=SIN1D(L)*FXBP+COS1D(L)*FYBP ELSE FXB=FXBP FYB=FYBP ENDIF C FRC1=.75*FXB FRC2=.75*FYB IF (K.EQ.3.OR.K.EQ.4) THEN Z=FRC2 FRC2=-FRC1 FRC1=Z ENDIF C FRC12T= DELT*(FRC1+BIOBR*FRC2) C ENDED STRESS COMPUTATIONS. ELSE ! if jump .eq. 1 C C INITIALIZE FLOW WITH GRAVITY AND WITHOUT SURFACE STRESS FRC12T=0. END IF ! if jump.ne.1 ! ! C 590 CONTINUE C TAKE AVERAGE OF 'AI' ON TWO SQUARES AINTRP=0. DO 620 KK=1,2 IDD=AMIN1(299.,DPH1D(KK)) IF (IDD.GT.1) THEN SLPAI=AI(IDD+1)-AI(IDD) DPHI=IDD AITR=AI(IDD)+SLPAI*(DPH1D(KK)-DPHI) ELSE AITR=AI(1)+SLPA*(DPH1D(KK)-1.) END IF ! ! 610 AINTRP=AINTRP+AITR 620 CONTINUE AINTRP=.5*AINTRP FLWCOF=X1B(4)*(AINTRP-AR(ID)*BIOBR)+1. FLW=FLWT+FLWCOF*FLWSQR(L)+FRC12T FWD=HSUB(II,JJ)-HSUB(I,J) ! IF (FWD.EQ.0.) THEN FLWSQR(L)=FLW ELSE FD2=-SIGA(K)*FWD*X1B(2)/(DPH*DPH) FD2=ABS(FD2) Z=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J) Z1=ELPDL2(II)+(ELPCL2(II)-ELPDL2(II))*SINL2(JJ) XY=.5*(1./Z+1./Z1) FD2=FD2*XY FLW1=FLW DO ITR=1,20 FLW1=FLW/(1.+FD2*ABS(FLW1)) END DO FLWSQR(L)=FLW1 END IF 630 CONTINUE C++++++++END MOMENTUM COMPUTATIONS FOR 1-DIM FLOW+++++++++++++++++++++++ C C--------END MAIN 'LOOP'------------------------------------------------ RETURN END SUBROUTINE FRCPNT(X,Y,FXB,FYB,NCATG,CSHLTR,NPLS) C JELESNIANSKI SEPTEMBER 1980 TDL IBM 360/195 C C PURPOSE C THIS SUBROUTINE (FRCPNT) COMPUTES COMPONENTS OF SURFACE C DRIVING FORCES AT A MOMENTUM GRID POINT. ENTRY VALUES FOR C COMPUTATIONS ARE THRU 'COMMON/MF/'. THE LAKE WINDS ARE C DISTORTED ACCORDING TO FORWARD SPEED ALONG TRACK; OCEAN C WINDS ARE NOT DISTORTED. C DATA SET USE C NONE C C VARIABLES C CSHLTR = EXTINCTION COEF OF STRESS FOR WATER LESS THAN 1 FT C AX AY = COMPONENTS OF TOTAL STORM TRAVERSE C C1 C2 = INITIAL (X,Y) STORM POSITION, SET IN 'INTVAL' C X12(4) = RADIUS OR MAX WIND IN FEET C X12(15) = SQUARE OF X12(4) C X12(9) X12(10) = COMPONENTS OF FORWARD SPEED, 1ST HOUR C C19 = FRIC CFFCNT (3X10-6) C W1218(2) = TWICE THE MAXIMUM WINDS, SET IN 'CMPUTE' C CSHLTR = STRESS SHELTERING COEFFICIENT, FOR DEPTHS<1 FT C CW(800) SW(800) = INFLOW ANGLE, MILE INCRMNT, LAKE WINDS C C(800) S(800) = INFLOW ANGLE, MILE INCRMNT, OCEAN WINDS C BI BR CI CR(800) = FRICTION FACTORS C P(800) = PRESSURE GRADIENT, MILE INTERVALS C FXX FYY = COMPONENTS OF STRESS C PXX PYY = COMPONENTS OF SURFACE PRESSURE GRADIENT FORCE C FXBP FYBP = COMPONENTS OF SURFACE DRIVING FORCES IN (X,Y) C FXB FYB = COMPONENTS OF SURFACE FORCES ON IMAGE PLANE C C GENERAL COMMENTS C THIS SUBROUTINE LIES IN OVERLAY 'CMPUTE'. IT IS CALLED IN C BY SUBROUTINE 'MOMNTM'. SUBROUTINE 'MOMNTM' LIES IN C SUBROUTINE 'CMPUTE'. C C COMMON /SLAT/ FSOUTH parameter (ndp=600) COMMON /FRST/ X1(50),X12(50) COMMON /FRTH/ IC12,IC19,BT COMMON /THRD/ C7,C17,C19,C25 COMMON /STRMSB/ C1,C2,C21,C22,AX,AY,PTENCY,RTENCY COMMON /DUMY44/ SW(800),CW(800),W1218(2),WMAX(2) COMMON /DUMMY4/ S(800),C(800),P(800),DELP(800) COMMON /SCND/ AR(ndp),AI(ndp),BR(ndp),BI(ndp),CR(ndp),CI(ndp) COMMON /MF/ DPH,HPD(4),ID C C COMPUTING FORCE AT A GRID POINT XP=X-C1-AX YP=Y-C2-AY RSQ=XP*XP+YP*YP RS=SQRT(RSQ) R1=RS/5280.+1. K=R1 R2=K DR=R1-R2 C C VECTOR WIND IS CONSTANT FOR DISTANCE >790 MLS FRM STM CENTER K=MIN0(K,790) CCN=X12(15)+RSQ X1218=W1218(NCATG) IF (NCATG.EQ.1) THEN C COMPUTE LAKE WINDS CK1=CW(K)+DR*(CW(K+1)-CW(K)) SK1=SW(K)+DR*(SW(K+1)-SW(K)) ELSE C COMPUTE OCEAN WINDS CSHLTR=1. CK1=C(K)+DR*(C(K+1)-C(K)) SK1=S(K)+DR*(S(K+1)-S(K)) ENDIF C A=RS*X12(9) -X1218*(YP*CK1+XP*SK1) B=RS*X12(10)+X1218*(XP*CK1-YP*SK1) C C FSOHTH=-1. FOR SOUTHERN HEMISPHERE, CLOCKWISE CIRCLULATION; C = 1. FOR NORTHERN HEMISPHERE, COUNTER-CLOCKWISE CIRCULATION. C c A=RS*X12(9) +X1218*(-FSOUTH*YP*CK1-XP*SK1) c B=RS*X12(10)+X1218*( FSOUTH*XP*CK1-YP*SK1) C COMPUTE WIND DISTORTION, FROM LAND EFFECTS, FOR LAKE WINDS IF (NCATG.EQ.1) THEN RHOL=ABS(XP*X12(20)+YP*X12(21)) CCC=CCN*RHOL/(X12(15)+RHOL*RHOL) A=A+CCC*X12(23) B=B+CCC*X12(24) ENDIF C WIND COMPONENTS (A,B), FT/SEC, IN X,Y DIRECTIONS A=A*X12(4)/CCN B=B*X12(4)/CCN C STRESS =DRAG COEFF* SPEED * WIND VECTOR C STRESS REDUCTION FOR WATER LESS THAN 1 FT, AND SPIN-UP PERIOD FCTT=C19*SQRT(A*A+B*B) FCTT=FCTT*CSHLTR*BT FXX=FCTT*A FYY=FCTT*B C FXB=CR(ID)*FXX-CI(ID)*FYY FYB=CR(ID)*FYY+CI(ID)*FXX C PRESSURE TERM IN THE FORCING FUNCTIONS C SKIPS IF INTERMEDIATE WATER B. C IS USED. IF (NPLS.NE.0) THEN AA=-DPH*(P(K)+DR*(P(K+1)-P(K))) PXX=AA*XP PYY=AA*YP FXB=FXB+BR(ID)*PXX-BI(ID)*PYY FYB=FYB+BR(ID)*PYY+BI(ID)*PXX ENDIF C RETURN END SUBROUTINE MOMNTM (ITIDE,INEST,IREX) C JELESNIANSKI,CHEN SEPTEMBER 1980 TDL IBM 360/195 C C PURPOSE C THIS SUBROUTINE (MOMNTM) COMPUTES FUTURE MOMENTUM COMPONENTS C ON ALL INTERIOR GRID POINTS, AFTER FIRST COMPUTING FUTURE C SURGE BY SUBROUTINE 'CONTTY'. IT CHECKS IF MOMNTM CAN EXIST C AT A GRID POINT; E.G., OVERTOPPING THE GRID POINT . C FUTURE SURGES UPDATED FROM SUB. 'CONTTY' ARE USED TO COMPUTE C HEIGHT GRADIENTS AND TO CHECK FOR ELIGIBLE MOMENTUM POINTS-- C BACKWARD SCHEME---. CHECKS ARE MADE IF MOMENTUM IS TO BE C COMPUTED AT OPEN BOUNDARY POINTS. C C DATA SET USE C NONE C C VARIABLES C NCATG = 1 FOR LAKE WINDS, 2 FOR OCEAN WINDS C NPLS = 0 REMOVE STATIC PRESSURE TERM IN FORCING(INT.B.C.) C = 1 OTHERWISE C JMXB1 = MAX NO. OF J-POINTS, LESS A BOUNDARY POINT C JUMP = 1, FULLY FLOODED; 2,PARTIALLY FLOODED (NO STRESS) C MS( ) = BEGIN I-SUBSCRIPT ON A J-LINE C ITREE( , ) = INDICATOR ARRAY FOR LAKE/OCEAN, TREE/NO TREE C A,B,C,D FOR BNOUDARY TYPES C ZBM( , ) = MAXIMUM BARRIER HT AT A SURGE POINT C IH(4) JH(4) = SHIFT I/J-SUBSCRIPTS TO 4 SURRNDNG SRGE PTS C IP(4) JP(4) = SHIFTS I/J-SUBPTS FROM SURGE PTS TO 4 MOM. CORNERS C HB( , ) = SURGE HEIGHTS C ZB( , ) = DEPTHS C HF(4) = STORAGE FOR 4 SURRNDNG SURGE PTS, C HPD(4) = STORAGE FOR 4 SURRNDNG TOTAL DEPTH POINTS C UB VB( , ) = COMPONENTS OF TRANSPORT C DHX DHY = COMPONENTS OF SURFACE GRADIENT C ID DPH = MEAN TOTAL DEPTH C AI(300) AR(300) = FRICTION COEFFICIENTS ON TRANSPORTS C FXB FYB = COMPONENTS OF DRIVING FORCES C X1B(4) = DEL*CORIOLIS PARM C CSHLTR = EXTINCTION COEF OF STRESS FOR WATER LESS THAN 1 FT C ISKIP = 1, COMPUTE MOMENTUM ON J=1 BOUNDARY C 2, COMPUTE MOMENTUM ON I=2 BOUNDARY C 3, COMPUTE MOMENTUM ON J=JMXB BOUNDARY C 4, COMPUTE MOMENTUM ON I=IMXB BOUNDARY C IJUMP = 1, SHALLOW WATER BOUNDARY COMPUTATIONS C 2, INTER DEPTH BOUNDARY COMPUTATIONS C ISBS JSBS(4) = SHIFT SUBSCRIPTS FROM INTERIOR TO BOUNDARY PT C ISHF(3,4) = SUBSCIPTS SHIFTS FROM INTERIOR TO 3 CORNER C MOM-POINTS FOR ALL 4 CORNERS C C GENERAL COMMENTS C THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS CALLED C IN BY SUBROUTINE 'CMPUTE'. THE PRESENT SUBROUTINE CALLS C IN 2 SUBROUTINES, 'FRCPNT' AND 'MNTMBD'. C C INCLUDE 'parm.for' C parameter (ndp=600) COMMON /FRTH/ IC12,IC19,BT COMMON /MF/ DPH,HPD(4),ID COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /EGTH/ DELS,DELT,G,COR COMMON /SCND/ AR(ndp),AI(ndp),BR(ndp),BI(ndp),CR(ndp),CI(ndp) COMMON /DUMB8/ IP(4),JP(4),IH(4),JH(4) COMMON /DUM88/ IIH(4),JJH(4),IHH(4),JHH(4) COMMON /DUMB4/ X1B(10) COMMON /DUMB4I/ CORL(90) COMMON /GPRT1/ DOLLAR,EBSN Cbeg current output by Jindong/Brian COMMON /POLAR/ DEGREE,RMOUTH,AZMTH,DELA COMMON /INDX/ IRB,IRE,NDB,NDE COMMON /CURRENT/ CurrentMag(M_,N_),U_Current(M_,N_), 1 V_Current(M_,N_) Cend current output by Jindong/Brian DIMENSION SLPAI(ndp),HFX(4),ISBS(4),JSBS(4),ISHF(3,4) CHARACTER*2 DOLLAR,SIG_ CHARACTER*1 EBSN DATA ISBS/0,-1,0,1/, JSBS/-1,0,1,0/ DATA SIG_/'$'/,SHGT/50./ DATA ISHF/1,2,4, 1,2,3, 1,3,4, 2,3,4/ C PRECALCULATED SLOPES OF AI(ID). MAY BE PERMANENTLY STORED. DO 10 I=1,ndp-1 10 SLPAI(I)=AI(I+1)-AI(I) SLPAI(ndp)=SLPAI(ndp-1) Cd Debug Huiqing Liu/MDL Jan. 2014 for initializing the current magnitude to 99.9 DO I=1,IMXB DO J=1,JMXB CurrentMag(I,J) = 99.9 ENDDO ENDDO cd Debug Huiqing Liu /MDL Jan. 2014 JEND=JMXB1 IF (DOLLAR.EQ.SIG_) JEND=JMXB C DO 570 J=2,JEND C C ISKIP SET 3 RANGES FOR J C ISKIP=2 C FOR CILCULAR ISLAND, NO SIDE-BOUNDARIES (1 OR 3) TO BE TESTED. C FOR J=1, U/V REPEAT THOSE AT J=JMXB (JEND). IF (DOLLAR.NE.SIG_) THEN IF(J.EQ.2) ISKIP=1 IF(J.EQ.JEND) ISKIP=3 ENDIF C IST=MS(J) IFN=ME(J) C SKIP A J-COLUMN IF MS(J) IS SET TO IMXB OR GREATER. IF ((IFN-IST).LT.0) GO TO 570 C DO 560 I=IST,IFN JUMP=1 C NO TERRAIN HIGHER THAN 35 FT CAN BE FLOODED. C CHANGED BY NSM 11/20/2010 TO ALLOW TERRAIN HIGHER THAN 35FT TO FLOOD C DID NOT Accept. The '4' is now correctly shifted to 56 and above C and if it was a '4' for 35..56, it is now a '1'. IF (ITREE(I,J).EQ.'4'.OR.ITREE(I,J).EQ.'6') GO TO 560 C IF (ITREE(I,J).EQ.'6') GO TO 560 C ITREE IS SET TO 1 OR 3 FOR LAKE WINDS, 2 OR 5 FOR OCEAN WINDS IF (ITREE(I,J).EQ.'2'.OR.ITREE(I,J).EQ.'5') THEN NCATG=2 ELSE NCATG=1 ENDIF IF (NCATG.EQ.1) GOTO 140 IF (ITREE(I,J).EQ.'2'.AND.ZBM(I,J).LT.-50.) GOTO 160 C C TEST IF SURRNDNG SQRS ARE DRY, PART/FULL FLOOD . C 140 CONTINUE c HMN=AMIN1(HB(I,J),HB(I-1,J),HB(I,J-1),HB(I-1,J-1)) IF (HMN.GT.ZBM(I,J)) GO TO 200 HMXX=AMAX1(HB(I,J),HB(I-1,J),HB(I,J-1),HB(I-1,J-1)) IF (HMXX.GT.ZBM(I,J)) GO TO 300 C NO TRANSPORTS ARE ALLOWED, GO TO CHECK IF NEXT TO BOUNDARY. GO TO 220 C ALL SURROUNDING SQUARES ARE FLOODED, CALCULATE GRADIENT TREM C WITH CENTER DIFFERENCES 200 JUMP=1 DHX=-HB(I-1,J-1)+HB(I,J-1)-HB(I-1,J)+HB(I,J) DHY=-HB(I-1,J-1)-HB(I,J-1)+HB(I-1,J)+HB(I,J) GO TO 340 C C AT LEAST ONE SQUARE IS NOT OVERTOPPED. USE HYDRULIC HEAD FOR C SLOPE ON THAT SQUARE (1/4 EFFECT) AS APPROXMATION 300 JUMP=2 C DEFINE A MODIFIED SURFACE GRADIENT FOR A PARTIALLY FLOODED POINT HF1=AMAX1(HB(I-1,J-1),ZBM(I,J)) HF2=AMAX1(HB(I ,J-1),ZBM(I,J)) HF3=AMAX1(HB(I-1,J ),ZBM(I,J)) HF4=AMAX1(HB(I ,J ),ZBM(I,J)) DHX=-HF1+HF2-HF3+HF4 DHY=-HF1-HF2+HF3+HF4 C C TAKE AVERAGE (D+H), OR TOTAL DEPTH, OF 4 SURROUNDING SQUARES C USED TO DEFINE BOTTOM STRESS COEFFICIENTS 340 CONTINUE HPD(1)=HB(I-1,J-1)+ZB(I-1,J-1) HPD(2)=HB(I ,J-1)+ZB(I ,J-1) HPD(3)=HB(I-1,J )+ZB(I-1,J ) HPD(4)=HB(I ,J )+ZB(I ,J ) DPH=.25*(HPD(1)+HPD(2)+HPD(3)+HPD(4)) DPH=MAX(1.,DPH) DPHZ=DPH ID=AMIN1(ndp-1.0,DPHZ) C C FOR VERY THIN SHEET OF WATER ON BARE TERRAIN, BOTTOM FRICTIONS C ARE ARTIFICIALLY INCREASED WITH A CUBIC FORMULA IF (DPH.GT.5.) GO TO 405 C IT IS UNNECESSARY IF PARTIALLY FLOODED. IF (JUMP.EQ.2) GO TO 400 C C REVISE DPH AND ID DO 380 K=1,4 IF(HPD(K).GE.5.) GO TO 370 C USE ALTERNATIVE X-1+(1-X/5)**3 = X*(.4+.12*X-.008*X**2) HPD(K)=HPD(K)*(.4+.12*HPD(K)-.008*HPD(K)*HPD(K)) GO TO 380 370 HPD(K)=HPD(K)-1. 380 CONTINUE C DPH=.25*(HPD(1)+HPD(2)+HPD(3)+HPD(4)) DPHZ=DPH ID=AMIN1(ndp-1.,DPHZ) 400 ID=MAX0(1,ID) C C+++++++ TAKE AVERAGE 'AI(TOTAL DEPTH)' ON 4 SQUARES++++++++++++++++++++ C EXTRAPOLATION IF DPH IS LESS THAN 1 FT 405 AINTRP=0. DO 430 K=1,4 DPHZ=HPD(K) IDD=MIN(ndp-1.,DPHZ) IDD=MAX0(IDD,1) C EXTRAPOLATION RESULTS IF DEPTH IS LESS THAN 1 FT. DPHI=IDD AITR=AI(IDD)+SLPAI(IDD)*(DPHZ-DPHI) 420 AINTRP=AINTRP+AITR 430 CONTINUE AINTRP=.25*AINTRP GOTO 450 C C ALL SURROUNDING SQUARES ARE FLOODED 160 CONTINUE JUMP=1 DHX=-HB(I-1,J-1)+HB(I,J-1)-HB(I-1,J)+HB(I,J) DHY=-HB(I-1,J-1)-HB(I,J-1)+HB(I-1,J)+HB(I,J) DPH=.25*(HB(I-1,J-1)+HB(I,J-1)+HB(I-1,J)+HB(I,J)+ 1 ZB(I-1,J-1)+ZB(I,J-1)+ZB(I-1,J)+ZB(I,J)) DPH=MAX(1.,DPH) DPHZ=MIN(ndp-1.,DPH) ID=MIN(ndp-1.,DPHZ) DPHI=ID AINTRP=AI(ID)+SLPAI(ID)*(DPHZ-DPHI) C 450 AII=1.+X1B(4)*AINTRP C C CORIOLIS PARAMETER DEPENDS ON THE THE LATITUDE OF THE GRID. LLAT=YLT(I,J) ZCOR=CORL(LLAT+1)+(YLT(I,J)-LLAT)*(CORL(LLAT+2)-CORL(LLAT+1)) ARJ=AR(ID)*X1B(4)*ZCOR C C AIU=AII*UB(I,J) AIV=AII*VB(I,J) FCT=X1B(6)*DPH C C IF PARTIALLY FOODED, NO STRESS IS APPLIED FOR CURRENT TIME. IF (JUMP.EQ.2) GOTO 490 C ----------------------------------------------------------------------- IF (NCATG.EQ.2) THEN CSHLTR=1. ELSE HFX(1)=HB(I-1,J-1)-ZBM(I,J) HFX(2)=HB(I ,J-1)-ZBM(I,J) HFX(3)=HB(I-1,J )-ZBM(I,J) HFX(4)=HB(I ,J )-ZBM(I,J) C C SHELTERING WIND STRESS LINEARLY IF WATER IS LESS THAN 1 FT C ABOVE TERRAIN OR BARRIER. IF (ITREE(I,J).NE.'1') THEN CSH1=AMIN1(HFX(1),1.) CSH2=AMIN1(HFX(2),1.) CSH3=AMIN1(HFX(3),1.) CSH4=AMIN1(HFX(4),1.) CSHLTR=.25*(CSH1+CSH2+CSH3+CSH4) ELSE C C STRONGER CSHLTERING FOR MANGROVES OR TREES. SHGT = HEIGHT OF C TREE TOP. CSHLTR=0. DO 484 KX=1,4 HFXX=1. IF(HFX(KX).LE.0.) THEN HFXX=0. ELSE IF(HFX(KX).LT.SHGT) THEN HFXX=(HFX(KX)/SHGT)**2 ENDIF CSHLTR=CSHLTR+HFXX 484 CONTINUE CSHLTR=.25*CSHLTR ENDIF ENDIF C C LOCATE GRID (I,J) ON PHYSICAT PLANE Z=(XR,YR)=ZETA XR=ELPCT(I)*COST(J) YR=ELPDT(I)*SINT(J) NPLS=1 c ! ! Huiqing.Liu /MDL modified to ETSS which need read GFS wind instead of ! using parmametric winds June 2014 ! ! CALL FRCPNT(XR,YR,FXBP,FYBP,NCATG,CSHLTR,NPLS) CALL FRCPNT1(I,J,FXBP,FYBP,NCATG,CSHLTR,NPLS) ! ! Huiqing.Liu /MDL ! C C CHANGE FORCING FROM (X,Y) SYSTEM BACK TO IMAGE PLANE, C CONJ(DZ/DZETA)*(FXBP,FYBP), WHERE DZ/DZETA=(XR,YR) XR=ELPDT(I)*COST(J) YR=ELPCT(I)*SINT(J) FXB=(XR*FXBP+YR*FYBP) FYB=(XR*FYBP-YR*FXBP) C XTTX=DELT*FXB XTTY=DELT*FYB GO TO 480 C ----------------------------------------------------------------------- C NO FORCINGS 490 XTTX=0. XTTY=0. C 480 CONTINUE C IF (IREX.EQ.1) THEN ! only activate this part of code when it is necessary to save ! computing time Huiqing Liu/MDL CBeg output current by Jindong/Brian c ! Calculate the grid cell distance in feet r_ft=SQRT(ELPDT(I)**2+(ELPCT(I)**2-ELPDT(I)**2)*SINT(J)**2) c ! Determine the angle rotations AngU_North=(AZMTH*(3.14159265358979/180.))-((J-NDB)*DELA) AngU_East=AngU_North-(3.14159265358979/2.) AngV_North=AngU_East AngV_East=AngV_North-(3.14159265358979/2.) c ! Calcualte the transport components U,V U_Current(I,J)=((UB(I,J)*COS(AngU_East)) 1 +(VB(I,J)*COS(AngV_East)))/r_ft V_Current(I,J)=((UB(I,J)*COS(AngU_North)) 1 +(VB(I,J)*COS(AngV_North)))/r_ft c ! Calculate the currents U,V,Magnitude U_Current(I,J)=U_Current(I,J)/(HB(I,J)+ZB(I,J)) V_Current(I,J)=V_Current(I,J)/(HB(I,J)+ZB(I,J)) CurrentMag(I,J)=SQRT(U_Current(I,J)**2+V_Current(I,J)**2) Cend output current by Jindong/Brian ENDIF C TEMPORARILY STORE PRESENT U/V TO BE USED FOR BOUNDARY POINTS UBJ=UB(I,J) VBJ=VB(I,J) C COMPUTE MOMENTUM ON INTERIOR POINT, SIELECKI'S SCHEME FOR C CORIOLIS TERM IS USED. UB(I,J)=AIU-FCT*(BR(ID)*DHX-BI(ID)*DHY)+ARJ*VB(I,J)+XTTX VB(I,J)=AIV-FCT*(BR(ID)*DHY+BI(ID)*DHX)-ARJ*UB(I,J)+XTTY C GO TO 500 C C TRANSPORT POINT NOT OVERTOPPED IN FUTURE 220 UB(I,J)=0. VB(I,J)=0. C ------------------------------------------------------------------ C CHECK IF MOMENTUM COMPUTATIONS REQUIRED ON BOUNDARY C ISKIP=1, LEFT BOUNDARY(J=1) C =2, AND I=2, NEXT TO TOP BOUNDARY C =3, RIGHT BOUNDARY C =4, AND I=IMXB1,NEXT TO BOTTOM BOUNDARY C IJUMP=1, INTERMEDIATE WATER; =2, SHALLOW WATER C 500 IJUMP=1 ! !------------------------------------------------------------ ! rewrite the following codes to remove the computed GOTO ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! Old codes !--------------------------------------------------------------- ! GO TO (530,520,540),ISKIP C C CHECK FOR 4 CORNER POINTS AND SET ICOR. C ICOR=1,LEFT-TOP; 2,LEFT-BOTTOM;3,RIGHT-TOP;4,RIGHT-BOTTOM ! 530 IF(I.NE.ME(J)) GO TO 512 ! ICOR=2 ! 512 IF(I.NE.MS(J)) GO TO 510 ! ICOR=1 ! GO TO 600 ! 540 IF(I.NE.ME(J)) GO TO 514 ! ICOR=4 ! 514 IF(I.NE.MS(J)) GO TO 510 ! ICOR=3 ! GO TO 600 C C TEST FOR BOTTOM BOUNDARY ! 520 IF(I.EQ.MS(J)) GO TO 510 ! IF(I.EQ.ME(J)) GO TO 516 ! GO TO 560 ! 516 ISKIP=4 !-------------------------------------------------------------- ! New codes !--------------------------------------------------------------- IF (ISKIP == 1) THEN IF (I.EQ.ME(J)) THEN ICOR = 2 END IF IF (I.NE.MS(J)) THEN ! 510 IA=ISBS(ISKIP) JA=JSBS(ISKIP) ISB=I+IA JSB=J+JA !C IF AT (I,J),TRANSPORTS ARE ZERO,SET TRANSPORTS AT BOUNDY ZEROS. IF (UB(I,J).EQ.0.) THEN UB(ISB,JSB)=0. VB(ISB,JSB)=0. CYCLE ENDIF !C PRE-SET ITREEALONG BOUNDARIES; =4 FOR LAND (SKIP COMP.) !C =9 FOR SHALLOW WATER !C =8 FOR INTERMEDIATE DEPTHS !C =6 FOR DEEP WATER (SKIP COMP.) IF (ITREE(ISB,JSB).EQ.'6'.OR.ITREE(ISB,JSB).EQ.'4') THEN CYCLE END IF IF (ITREE(ISB,JSB).EQ.'8' ) THEN IJUMP=2 END IF !C WRITE(*,*) ' MNTMBD ISB,JSB',ISB,JSB,IJUMP IF (INEST.NE.2.AND.ITIDE.NE.31) THEN ! ALL boundary set to water grid Added by Huiqing Liu/MDL CALL MNTMBD(JUMP,IJUMP,ISB,JSB,ARJ,AII,I,J,UBJ,VBJ, 1 FXB,FYB,NCATG,CSHLTR,NPLS) ENDIF CYCLE ELSE ICOR = 1 DO KK=1,3 K=ISHF(KK,ICOR) ISB=I+IH(ICOR)+IP(K) JSB=J+JH(ICOR)+JP(K) !C IF AT (I,J),TRANSPORTS ARE ZERO,SET TRANSPORTS AT BOUNDY ZEROS. IF (UB(I,J).EQ.0.) THEN UB(ISB,JSB)=0. VB(ISB,JSB)=0. CYCLE END IF C IF (ITREE(ISB,JSB).EQ.'6'.OR.ITREE(ISB,JSB).EQ.'4') THEN EXIT END IF IF (ITREE(ISB,JSB).EQ.'8' ) THEN IJUMP=2 END IF IF (INEST.NE.2.AND.ITIDE.NE.31) THEN ! ALL boundary set to water grid Added by Huiqing Liu/MDL CALL MNTMBD(JUMP,IJUMP,ISB,JSB,ARJ,AII,I,J,UBJ,VBJ, 1 FXB,FYB,NCATG,CSHLTR,NPLS) END IF END DO END IF ELSEIF (ISKIP == 2) THEN IF (I.EQ.MS(J).OR.I.EQ.ME(J)) THEN IF (I.EQ.ME(J)) THEN ISKIP = 4 END IF IA=ISBS(ISKIP) JA=JSBS(ISKIP) ISB=I+IA JSB=J+JA !C IF AT (I,J),TRANSPORTS ARE ZERO,SET TRANSPORTS AT BOUNDY ZEROS. IF (UB(I,J).EQ.0.) THEN UB(ISB,JSB)=0. VB(ISB,JSB)=0. CYCLE ENDIF !C PRE-SET ITREEALONG BOUNDARIES; =4 FOR LAND (SKIP COMP.) !C =9 FOR SHALLOW WATER !C =8 FOR INTERMEDIATE DEPTHS !C =6 FOR DEEP WATER (SKIP COMP.) IF (ITREE(ISB,JSB).EQ.'6'.OR.ITREE(ISB,JSB).EQ.'4') THEN CYCLE END IF IF (ITREE(ISB,JSB).EQ.'8' ) THEN IJUMP=2 END IF !C WRITE(*,*) ' MNTMBD ISB,JSB',ISB,JSB,IJUMP IF (INEST.NE.2.AND.ITIDE.NE.31) THEN ! ALL boundary set to water grid Added by Huiqing Liu/MDL CALL MNTMBD(JUMP,IJUMP,ISB,JSB,ARJ,AII,I,J,UBJ,VBJ, 1 FXB,FYB,NCATG,CSHLTR,NPLS) ENDIF CYCLE ELSE CYCLE END IF ELSEIF (ISKIP == 3) THEN IF (I.EQ.ME(J)) THEN ICOR = 4 END IF IF(I.NE.MS(J)) THEN IA=ISBS(ISKIP) JA=JSBS(ISKIP) ISB=I+IA JSB=J+JA !C IF AT (I,J),TRANSPORTS ARE ZERO,SET TRANSPORTS AT BOUNDY ZEROS. IF (UB(I,J).EQ.0.) THEN UB(ISB,JSB)=0. VB(ISB,JSB)=0. CYCLE ENDIF !C PRE-SET ITREEALONG BOUNDARIES; =4 FOR LAND (SKIP COMP.) !C =9 FOR SHALLOW WATER !C =8 FOR INTERMEDIATE DEPTHS !C =6 FOR DEEP WATER (SKIP COMP.) IF (ITREE(ISB,JSB).EQ.'6'.OR.ITREE(ISB,JSB).EQ.'4') THEN CYCLE END IF IF (ITREE(ISB,JSB).EQ.'8' ) THEN IJUMP=2 END IF !C WRITE(*,*) ' MNTMBD ISB,JSB',ISB,JSB,IJUMP IF (INEST.NE.2.AND.ITIDE.NE.31) THEN ! ALL boundary set to water grid Added by Huiqing Liu/MDL CALL MNTMBD(JUMP,IJUMP,ISB,JSB,ARJ,AII,I,J,UBJ,VBJ, 1 FXB,FYB,NCATG,CSHLTR,NPLS) ENDIF CYCLE ELSE ICOR=3 DO KK=1,3 K=ISHF(KK,ICOR) ISB=I+IH(ICOR)+IP(K) JSB=J+JH(ICOR)+JP(K) !C IF AT (I,J),TRANSPORTS ARE ZERO,SET TRANSPORTS AT BOUNDY !ZEROS. IF (UB(I,J).EQ.0.) THEN UB(ISB,JSB)=0. VB(ISB,JSB)=0. CYCLE END IF C IF (ITREE(ISB,JSB).EQ.'6'.OR.ITREE(ISB,JSB).EQ.'4') THEN EXIT END IF IF (ITREE(ISB,JSB).EQ.'8' ) THEN IJUMP=2 END IF IF (INEST.NE.2.AND.ITIDE.NE.31) THEN ! ALL boundary set to water grid Added by Huiqing Liu/MDL CALL MNTMBD(JUMP,IJUMP,ISB,JSB,ARJ,AII,I,J,UBJ,VBJ, 1 FXB,FYB,NCATG,CSHLTR,NPLS) END IF END DO END IF END IF ! ! ! 560 CONTINUE 570 CONTINUE 580 CONTINUE C FOR CIRCULAR ISLAND, PERIODIC BOUNDARY CONDITIONS ARE USED. IF (DOLLAR.EQ.'$') THEN DO 700 I=1,IMXB UB(I,1)=UB(I,JMXB) 700 VB(I,1)=VB(I,JMXB) ENDIF RETURN END SUBROUTINE MNTMBD(JUMP,IJUMP,ISB,JSB,ARJ,AII,I,J,UBJ,VBJ, 1 FXB,FYB,NCATG,CSHLTR,NPLS) C JELESNIANSKI OCTOBER 1980 TDL IBM 360/195 C C PURPOSE C THIS SUBROUTINE (MNTMBD) COMPUTES TRANSPORTS ON MOMENTUM C BOUNDARY POINTS, BETWEEN LAND AND DEEP WATERS. IN SHALLOW C WATERS THE BOUNDARY CONDITION EQUATES SURFACE GRADIENT AT C BOUNDARY AND NEAREST INTERIOR POINT. IN INTERMEDIATE C WATERS THE BOUNDARY CONDITION IGNORES TIME VARIATIONS AND C EQUATES SURGE GRADIENT AS STATIC HEIGHT GRADIENT AT A C BOUNDARY POINT. C BELOW IS ARRANGEMENT OF MOMENTUM BOUNDARY LINE C AND NEAREST INTERIOR LINE. C C . .(I,J) . .(I,J) . . INTERIOR LINE C * * C + + * + * C * * C . .<---IJUMP=2--->. .< IJUMP=1 >.****. BOUNDARY LINE C (ISB,JSB) (ISB,JSB) C DEEP---> <---LAND----- C C DATA SET USE C NONE C C VARIABLES C ISKIP = 1 TO COMPUTE MOMENTUM ON J=1 BOUNDARY C 2 TO COMPUTE MOMENTUM ON I=IMXB BOUNDARY C 3 TO COMPUTE MOMENTUM ON J=JMXB BOUNDARY C ISB JSB = REVISED SUBSCRIPTS FOR A BOUNDARY POINT C UB VB( , ) = TRANSPORTS C JUMP = 1 FOR SURFACE FORCING,2 FOR NO SURFACE FORCING C IJUMP = 1 FOR SHALLOW WATER BOUNDARY COMPUTATIONS C 2 FOR INTER DEPTH BOUNDARY COMPUTATIONS C FXB FYB = COMPONENTS OF DRIVING FORCES C AIT AII AIU = PRE-COMPUTED VALUES IN 'MOMNTM' C AIV ARJ XTT = DITTO C AINTRP = DITTO C X1B(10) = FIXED CONSTANTS C UBJ VBJ = TRANSPORTS AT (I,J) AT PRESENT TIME,SAVED BEFORE C UPDATED C C GENERAL COMMENTS C THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS CALLED C IN BY SUBROUTINE 'MOMNTM'. SUBROUTINE 'MOMNTM' IS CALLED C IN BY SUBROUTINE 'CMPUTE'. THIS SUBROUTINE ALSO CALLS IN C SUBROUTINE 'STRESS'. C INCLUDE 'parm.for' C parameter (ndp=600) COMMON /EGTH/ DELS,DELT,G,COR COMMON /DUMB4/ X1B(10) COMMON /SCND/ AR(ndp),AI(ndp),BR(ndp),BI(ndp),CR(ndp),CI(ndp) C AIUL=AII*(UB(ISB,JSB)-UBJ) AIVB=AII*(VB(ISB,JSB)-VBJ) C TEST FOR PARTIAL FLOODING (NO FORCING). IF (JUMP.EQ.2) THEN XTTX=0. XTTY=0. GO TO 140 ENDIF C TEST FOR SHALLOW WATER BOUNDARY AND INTERMEDIATE WATER BRY. IF (IJUMP.EQ.2) GOTO 160 C C SHALLOW WATER COMPUTATIONS, FORCING COMPUTED ON BOUNDRY PTS C FXD=FXB FYD=FYB C GRID TRANSFORMATION FROM (I,J) TO (X,Y) XR=ELPCT(ISB)*COST(JSB) YR=ELPDT(ISB)*SINT(JSB) NPLS=1 ! ! Huiqing.Liu /MDL modified to ETSS which need read GFS wind instead of ! using parmametric winds June 2014 ! ! CALL FRCPNT(XR,YR,FXBP,FYBP,NCATG,CSHLTR,NPLS) CALL FRCPNT1(ISB,JSB,FXBP,FYBP,NCATG,CSHLTR,NPLS) ! ! Huiqing.Liu /MDL ! C CHANGE FORCING FROM (X,Y) SYSTEM BACK TO IMAGE PLANE XR=ELPDT(ISB)*COST(JSB) YR=ELPCT(ISB)*SINT(JSB) FXB=XR*FXBP+YR*FYBP FYB=XR*FYBP-YR*FXBP C STRESS DIFFERENCE BETWEEN BOUNDARY POINT AND INTERIOR POINT XTTX=DELT*(FXB-FXD) XTTY=DELT*(FYB-FYD) C C COMPUTE MOMENTUM ON SHALLOW DEPTH BOUNDARIES. SURFACE GRADIENT C TERMS CANCELED BY ASSUMED B.C. C NOTE: USE SIELECKI'S SCHEME FOR CORIOLIS TERM. C VBJ = V-COMPONENT AT PRESENT TIME AT I,J C UB(I,J),VB(I,J) = U,V IN FUTURE TIME 140 UB(ISB,JSB)=UB(I,J)+AIUL+ARJ*(VB(ISB,JSB)-VBJ) +XTTX VB(ISB,JSB)=VB(I,J)+AIVB-ARJ*(UB(ISB,JSB)-UB(I,J))+XTTY GO TO 170 C C NO FORCIING WHEN JUMP=2 C C COMPUTE MOMENTUM ON INTERMEDIATE DEPTH BOUNDARIES C ELIMINATE STATIC PRESSURE TERM IN FRCPNT(NPLS=0). 160 NPLS=0 CSHLTR=1. C GRID TRANSFORMATION FROM (I,J) TO (X,Y) XR=ELPCT(ISB)*COST(JSB) YR=ELPDT(ISB)*SINT(JSB) C ! ! Huiqing.Liu /MDL modified to ETSS which need read GFS wind instead of ! using parmametric winds June 2014 ! ! CALL FRCPNT(XR,YR,FXBP,FYBP,NCATG,CSHLTR,NPLS) CALL FRCPNT1(ISB,JSB,FXBP,FYBP,NCATG,CSHLTR,NPLS) ! ! Huiqing.Liu /MDL ! C CHANGE FROM (X,Y) SYSTEM TO IMAGE PLANE XR=ELPDT(ISB)*COST(JSB) YR=ELPCT(ISB)*SINT(JSB) FXB=XR*FXBP+YR*FYBP FYB=XR*FYBP-YR*FXBP C 166 UB(ISB,JSB)=AII*UB(ISB,JSB)+ARJ*VB(ISB,JSB)+ DELT*FXB VB(ISB,JSB)=AII*VB(ISB,JSB)-ARJ*UB(ISB,JSB)+ DELT*FYB 170 RETURN END SUBROUTINE FHMXSV C CONVERT SURGE ENVELOP DATA WITH 999 FOR DRY SQUARES SO THAT ZB() C CAN BE FREE FROM THE MEMORY. INCLUDE 'parm.for' C COMMON /GPRT1/ DOLLAR,EBSN COMMON /SWTCH/ IOPERL(5) COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 CHARACTER*2 DOLLAR CHARACTER*1 EBSN C REPLACE HMX BY 999 IF THE SQUARE IS DRY. ! DO 1 J=1,JMXB ! DO 1 I=1,IMXB DO J=1,JMXB DO I=1,IMXB IF (HMX(I,J)+ZB(I,J).EQ.0) THEN IHMX(I,J)=999 ELSE IHMX(I,J)=HMX(I,J)*10.+.5 ENDIF ! 1 CONTINUE END DO END DO !------------------------------------------------------------ ! rewrite the following codes to remove the computed GOTO ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! Old codes !--------------------------------------------------------------- ! DO 20 J=1,JMXB ! DO 20 I=1,IMXB ! IF (IHMX(I,J).EQ.999) GOTO 20 ! IF (DOLLAR.EQ.'2$') IHMX(I,J)=(HMX(I,J)-ZSUB(J))*10.+.5 ! 20 CONTINUE DO J=1,JMXB DO I=1,IMXB IF (IHMX(I,J).EQ.999) cycle IF (DOLLAR.EQ.'2$') IHMX(I,J)=(HMX(I,J)-ZSUB(J))*10.+.5 END DO END DO IF(IOPERL(1).NE.2) THEN NGRP=(JMXB1-1)/25+1 DO 1130 NN=1,NGRP J1=1+(NN-1)*25 J2=J1+MIN0(24,JMXB1-J1 ) CC WRITE(*,180)(J,J=J1,J2) 180 FORMAT(/I6,26I5) C DO 1120 I=1,IMXB1 CC WRITE(*,'(1H ,I2,26I5)') I,(IHMX(I,J),J=J1,J2) 1120 CONTINUE 1130 CONTINUE ENDIF C RETURN END SUBROUTINE ARCHIV C ARCHIV SURGE ENVELOP DATA FOR DISPLAY INCLUDE 'parm.for' COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /IDENT/ AIDENT(40),DACLOK(7) COMMON /FLES/ FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1 CHARACTER*256 FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1 C ARCHIVE SURGE ENVELOP TO A FILE C OPEN (91,FILE=FLE91,FORM='UNFORMATTED',STATUS='UNKNOWN') CALL TOPEN (91,FLE91,1,2) C WRITE(91) IMXB,JMXB C WRITE(91) AIDENT CALL TWRITE (91,IMXB,1) CALL TWRITE (91,JMXB,1) CALL TWRITEC (91,AIDENT,160) C SAVE HMX*10 (10TH A FT ACCURACY) IN INTEGER FORMAT C WRITE(91) ((IHMX(I,J),I=1,IMXB),J=1,JMXB) DO J=1,JMXB DO I=1,IMXB CALL TWRITES (91,IHMX(I,J),1) END DO END DO C CLOSE (91) CALL TCLOSE (91) C RETURN END SUBROUTINE FILTRP C JELESNIANSKI OCTOBER 1980 TDL IBM 360/195 C C PURPOSE C THIS SUBROUTINE (FILTRP) SMOOTHS THE COMPUTED SURGE ENVELOPE C OF HIGHEST WATERS, IN SPACE. THE PROCEDURES ARE THE SILLILER C IN 'FILTER' AND 'HMXSV', 2 PASSES OF 3-POINT SMOOTHER. C IT IS FOR COSMATIC APPEARANCE. NO SMOOTHING ACROSS THE NON- C OVERTOPPED BARRIERS IS ALLOWED. C DATA SET USE C NONE C C VARIABLES C IMAX JMAX = MAXIMUM GRID POINTS IN THE BASIN C IMXB1 JMXB1 = IMAX LESS 1, JMXB LESS 1 C HMX( , ) = SURGE ENVELOPE C HZ( ,3) = DUMMY STORAGE FOR 3 LINES C ZB( , ) = DEPTH FIELD C IIH(4) JJH(4) = SHIFT SUBSCRIPTS TO 4 ADJOINING MOMENTUM POINTS C IP(4) JP(4) = SHIFT SUBSCRIPTS TO 4 ADJOINING SURGE POINTS C ZBM( , ) = MAX BARIER HEIGHTS AT MOMENTUM POINTS C C GENERAL COMMENTS C THIS SUBROUTINE RESIDES IN OVERLAY 'PLOTS'. IT IS CALLED C IN BY SUBROUTINE 'PLOTS'. C INCLUDE 'parm.for' C COMMON /FFTH/ ITIME,MHALT COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /DUMB8/ IP(4),JP(4),IH(4),JH(4) COMMON /DUM88/ IIH(4),JJH(4),IHH(4),JHH(4) COMMON /HTERAIN/ HTER DIMENSION HZ(M_,3) DIMENSION GAMF(4) DATA HT/0./ C C STORE 2 OF 3 INITIAL LINES FOR SMOOTHING MIDDLE LINE DO 300 IPS=1,2 KR=IPS C !------------------------------------------------------------ ! rewrite the following codes to remove the computed GOTO ! Huiqing.Liu /MDL Feb. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! Old codes !--------------------------------------------------------------- ! GO TO (104,102),IPS ! 102 I1=2 ! I2=IMXB2 ! J1=1 ! J2=JMXB1 ! DO 100 J=1,2 ! DO 100 I=1,IMXB1 ! 100 HZ(I,J+1)=HMX(I,J) ! GO TO 105 ! 104 I1=1 ! I2=IMXB1 ! J1=2 ! J2=JMXB2 ! DO 112 J=1,3 ! DO 112 I=1,IMXB1 ! 112 HZ(I,J)=HMX(I,J) !-------------------------------------------------------------- ! New codes !--------------------------------------------------------------- IF (IPS == 1) THEN I1=1 I2=IMXB1 J1=2 J2=JMXB2 DO J=1,3 DO I=1,IMXB1 HZ(I,J)=HMX(I,J) END DO END DO ELSEIF (IPS == 2) THEN I1=2 I2=IMXB2 J1=1 J2=JMXB1 DO J=1,2 DO I=1,IMXB1 HZ(I,J+1)=HMX(I,J) END DO END DO ENDIF ! ! ! 105 CONTINUE !------------------------------------------------------------ ! rewrite the following codes to remove the GOTO statments ! Huiqing.Liu /MDL March. 2015 !------------------------------------------------------------ !-------------------------------------------------------------- ! Old codes !--------------------------------------------------------------- DO 220 J=J1,J2 DO 190 I=I1,I2 C CHANGED BY NSM TO HTER FROM 35 11/20/2010 : Accepted 4/4/2011 ! IF(ZB(I,J).LE.(-HTER)) GO TO 190 IF(ZB(I,J).LE.(-HTER)) cycle HST=HMX(I,J)+ZB(I,J) ! IF(HST.LT.HT) GO TO 190 IF(HST.LT.HT) cycle HSM=0. HDIF=0. GAM0=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J) C C**************SMOOTHING ROUTINE**************************************** DO 150 K=KR,4,2 II=I+IIH(K) JK=J+JJH(K) ! IF (II.EQ.0.OR.JK.EQ.0) THEN ! GAMF(K)=0. ! GOTO 160 ! ENDIF ! GAMFK=ELPDL2(II)+(ELPCL2(II)-ELPDL2(II))*SINL2(JK) !c GAMF(K)=.5*(1.+GAMFK/GAM0)-1. ! Z=GAMFK/GAM0 ! GAMF(K)=2.*Z/(1.+Z)-1. ! IF (II.EQ.IMXB.OR.JK.EQ.JMXB) GO TO 160 ! JJ=JJH(K)+2 ! IF (HZ(II,JJ)+ZB(II,JK).LT.HT) GO TO 160 ! KK=MOD(K,4)+1 ! IA=I+IP(K) ! IB=I+IP(KK) ! JA=J+JP(K) ! JB=J+JP(KK) ! Z =ZBM(IA,JA) ! ZZ=ZBM(IB,JB) ! ZZZ=AMIN1(Z,ZZ)+HT ! IF( HZ(I,2).LT.ZZZ.OR.HZ(II,JJ).LT.ZZZ) THEN ! HSM=HSM+HZ(I,2) ! HZZ=HZ(I,2) ! ELSE ! HSM=HSM+HZ(II,JJ) ! HZZ=HZ(II,JJ) ! ENDIF ! GOTO 210 ! 160 CONTINUE ! HSM=HSM+HZ(I,2) ! HZZ=HZ(I,2) IF (II.EQ.0.OR.JK.EQ.0) THEN GAMF(K)=0. ! 160 CONTINUE HSM=HSM+HZ(I,2) HZZ=HZ(I,2) ELSE GAMFK=ELPDL2(II)+(ELPCL2(II)-ELPDL2(II))*SINL2(JK) Z=GAMFK/GAM0 GAMF(K)=2.*Z/(1.+Z)-1. IF (II.EQ.IMXB.OR.JK.EQ.JMXB) THEN ! 160 CONTINUE HSM=HSM+HZ(I,2) HZZ=HZ(I,2) ELSE JJ=JJH(K)+2 IF (HZ(II,JJ)+ZB(II,JK).LT.HT) THEN ! 160 CONTINUE HSM=HSM+HZ(I,2) HZZ=HZ(I,2) ELSE KK=MOD(K,4)+1 IA=I+IP(K) IB=I+IP(KK) JA=J+JP(K) JB=J+JP(KK) Z =ZBM(IA,JA) ZZ=ZBM(IB,JB) ZZZ=AMIN1(Z,ZZ)+HT IF( HZ(I,2).LT.ZZZ.OR.HZ(II,JJ).LT.ZZZ) THEN HSM=HSM+HZ(I,2) HZZ=HZ(I,2) ELSE HSM=HSM+HZ(II,JJ) HZZ=HZ(II,JJ) ENDIF END IF END IF END IF ! 210 CONTINUE C WEIGHTING ADJUSTMENTS IN I-DIRECTION DUE TO POLAR GRIDS. HDIF=HDIF+HZZ*GAMF(K) 150 CONTINUE C**************END SMOOTHING ROUTINE************************************ C HMX(I,J)=(2.*HMX(I,J)+HSM)/4. C 180 HMX(I,J)=HMX(I,J)+HDIF/4. HMX(I,J)=AMAX1(HMX(I,J),-ZB(I,J)) 190 CONTINUE C C UPDATE 3 LINES FOR SMOOTHING ! IF(J.GE.J2) GO TO 220 IF(J.GE.J2) cycle DO 200 II=1,IMXB1 HZ(II,1)=HZ(II,2) HZ(II,2)=HZ(II,3) 200 HZ(II,3)=HMX(II,J+2) 220 CONTINUE 300 CONTINUE RETURN END SUBROUTINE HMXSV C SEPTEMBER 1980 JYE CHEN TDL IBM 360/195 C C PURPOSE C THIS SUBROUTINE (HMXSV) STORES MAXIMUM SURGE AS ADVANCING C IN TIME AT EACH GRID POINT, INTO ARRAY HMX( , ). C THE 'SMOOTHED' SURGE VALUE IS COMPUTED AND COMPARED LOCALLY C AGAINST PREVIOUSLY SAVED HMX. THE REASON IS TO AVOID C TEMPORARY NOISE ENTERING THE FINAL ENVELOP. C C DATA SET USE C C VARIABLES C HMX( , ) = MAXIMUM SURGE AT EACH GRID PT C IP(4) JP(4) = SHIFTS OF I J FROM HEIGHT PT TO 4 MOMNTN CORNERS C IBB(6,M) = RANGES AND INCREMENTS FOR 3 DO-LOOPS OF I,J C M=1 INTERIOR,M=2 HORZNTL BNDY,M=3 VERTICAL BNDY C RANGES ON THREE BASIN SEGMENTS ARE: C SEE SUBROUTINE FILTER C HB( , ) = SURGE FIELD C ZB( , ) = DEPTH FIELD C IIH(4) JJH(4) = I/J-SHIFT SUBSCRIPTS FOR SURGE POINTS C C GENERAL COMMENTS C THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. IT IS CALLED C BY SUBROUTINE 'CMPUTE'. C A 5-POINT SMOOTHING OPERATOR SIMILLER TO ONE C IN SUBROUTINE FILTER IS USED, SEE 'FILTER'. C INCLUDE 'parm.for' COMMON /FFTH/ ITIME,MHALT COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /EGTH/ DELS,DELT,G,COR COMMON /DUMB8/ IP(4),JP(4),IH(4),JH(4) COMMON /DUM88/ IIH(4),JJH(4),IHH(4),JHH(4) COMMON /GPRT2/ EBSN1,EBSN2 CHARACTER*1 EBSN,EBSN1,EBSN2 C DIMENSION GAMF(4) DIMENSION IBB(6,3) C DATA IIH/0,1,0,-1/,JJH/-1,0,1,0/ DATA IBB(1,1),IBB(1,2),IBB(1,3)/2,1,2/, 1 IBB(3,1), IBB(3,3)/1, 1/, 2 IBB(4,1),IBB(4,2),IBB(4,3)/2,1,1/, 3 IBB(6,1),IBB(6,2) /1,1 / IBB(2,1)=JMXB2 IBB(2,3)=JMXB2 IBB(3,2)=JMXB2 IBB(2,2)=JMXB1 IBB(5,1)=IMXB2 IBB(6,3)=IMXB2 IBB(5,2)=IMXB1 IBB(5,3)=IMXB1 C DO 200 M=1,3 JL =IBB(1,M) JM =IBB(2,M) JNCR=IBB(3,M) IL =IBB(4,M) IM =IBB(5,M) INCR=IBB(6,M) C DO 190 J=JL,JM,JNCR C DO 180 I=IL,IM,INCR C NO UPDATING IF SURGE IS DECREASING. HX=HMX(I,J) IF(HX.GE.HB(I,J)) GOTO 180 C SMOOTHING IS UNNECESSARY FOR WATER DEPTH GREATER THAN 40 FT. IF (ZB(I,J).GT.40..OR.EBSN1.EQ.'&') GOTO 175 HST=HB(I,J)+ZB(I,J) C TEST FOR LAND WETTED LESS THAN 1 FT. IF(HST.LT.1.) GO TO 175 HSM=0. HDIF=0. GAM0=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J) C C******** 'LOOP' FOR SMOOTHING****************************************** DO 150 K=1,4 II=I+IIH(K) JJ=J+JJH(K) IF (M.EQ.1) GOTO 158 IF (II.EQ.0.OR.JJ.EQ.0) THEN GAMF(K)=0. GO TO 140 ENDIF 158 CONTINUE GAMFK=ELPDL2(II)+(ELPCL2(II)-ELPDL2(II))*SINL2(JJ) c GAMF(K)=.5*(1.+GAMFK/GAM0)-1. Z=GAMFK/GAM0 GAMF(K)=2.*Z/(1.+Z)-1. IF (JJ.EQ.JMXB.OR.II.EQ.IMXB) GO TO 140 C C THE NEIGHBORING SQUARE IS DRY. IF (HB(II,JJ)+ZB(II,JJ).LT.1.) GO TO 140 KK=MOD(K,4)+1 IA=I+IP(K) IB=I+IP(KK) JA=J+JP(K) JB=J+JP(KK) Z =ZBM(IA,JA) ZZ=ZBM(IB,JB) ZZZ=AMIN1(Z,ZZ)+1. C C WATER ON EITHER SIDE MUST EXCEED THE BARRIER(LOWER) BY AT C LEAST 1 FT. IF( HB(I,J).LT.ZZZ.OR.HB(II,JJ).LT.ZZZ) THEN HSM=HSM+HB(I,J) HZZ=HB(I,J) ELSE HSM=HSM+HB(II,JJ) HZZ=HB(II,JJ) ENDIF GOTO 212 140 CONTINUE HSM=HSM+HB(I,J) HZZ=HB(I,J) 212 CONTINUE C WEIGHTING ADJUSTMENTS IN I-DIRECTION DUE TO POLAR GRIDS. HDIF=HDIF+HZZ*GAMF(K) 150 CONTINUE C******* END 'LOOP' FOR SMOOTHING ************************************** C HTEMP=(4.*HB(I,J)+HSM+HDIF)/8. HMX(I,J)=AMAX1(HX,HTEMP) GOTO 180 175 HMX(I,J)=HB(I,J) 180 CONTINUE 190 CONTINUE 200 CONTINUE C write (*,*) 'hmxsv:E',HMX(17,18),HB(17,18),HMX(64,21),HB(64,21) RETURN END SUBROUTINE SMPT2G C J. CHEN , AUG 22 1984 TDL IBM 360/195 C PURPOSE C TO COLLECT MOMENTUM GRID POINTS AJACENT TO ANY INTERIOR C BOUNDARIES (U=V=0) THEN CONVERT TO HEIGHT POINTS. C REPETITION IS AVOIDED. SPECIAL POINTS, ONE OF TH 4 IS C DRY, ARE TAGGED AS 888. FOR HEAVIER SMOOTHING. C CALL SUBROUTINE 'FLT2G' FOR SMOOTHING. C DATA SET USE C NONE C VARIABLES C HSUB( , ) = TEMPORARY STORAGE SPACE FOR HEIGHT FIELD C INCLUDE 'parm.for' C MCT_ set to 8000 06/09/10 PARAMETER (MCT_=8000) C PARAMETER (MCT_=1500) C COMMON /FLWCPT/ NSQRS,NSQRW,NSQRWC,NPSS,NCUT COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /SM2G/ MM2G,MCMX,I2G(M2G_),J2G(M2G_) COMMON /DUMB8/ IP(4),JP(4),IH(4),JH(4) COMMON /DUM88/ IIH(4),JJH(4),IHH(4),JHH(4) COMMON /FFTH/ ITIME,MHALT INTEGER*2 MM2G,MCMX,I2G,J2G INTEGER*2 MCT(MCT_),I2GG(MCT_),J2GG(MCT_) M=0 DO 10 J=2,JMXB2 DO 10 I=2,IMXB2 IF (UB(I,J).EQ.0..OR.ZBM(I,J).LE.-30.) GOTO 10 IF (UB(I-1,J ).EQ.0.) GOTO 20 IF (UB(I ,J-1).EQ.0.) GOTO 20 IF (UB(I ,J+1).EQ.0.) GOTO 20 IF (UB(I+1,J ).EQ.0.) GOTO 20 GOTO 10 20 IF (M+1.GE.M2G_) GOTO 22 M=M+1 I2G(M)=I J2G(M)=J 10 CONTINUE GOTO 23 C 22 CONTINUE CC WRITE(*,998)I2G(M),J2G(M) 998 FORMAT(' TOTAL MOMN. POINTS FOR SMOOTHING REACHES 3000.', 1 ' LAST I,J =',2I8) STOP C SAVE THE TOTAL COUNT OF CORNER MOMN. POINTS 23 MM2G=M C write (*,*) "num momn points for smooth", mm2g C DO 100 J=1,JMXB1 DO 100 I=1,IMXB1 C C COLLECT MOMNTURMUMTURM POINTS NEXT TO THE INTERIOR BOUNDARY OF C COMPUTATIONS 100 HSUB(I,J)=999. C MC=0 DO 120 M=1,MM2G II=I2G(M) IF(II.LE.1.OR.II.GE.IMXB) GO TO 120 JJ=J2G(M) IF(JJ.LE.1.OR.JJ.GE.JMXB) GO TO 120 DO 110 K=1,4 I=II+IH(K) J=JJ+JH(K) IF(HB(I,J)+ZB(I,J).EQ.0.) GO TO 112 HSUB(I,J)=HB(I,J) GO TO 110 112 IF(MC+1.GE.MCT_) GO TO 113 MC=MC+1 MCT(MC)=M 110 CONTINUE 120 CONTINUE GOTO 125 113 CONTINUE CC WRITE(*,997) I2G(M),J2G(M) 997 FORMAT(' TOTAL SPECIAL CORNER POINTS (MOMN) REACHES 1500.', 1 ' LAST I,J =',2I8) STOP 125 MCMX=MC C write (*,*) "num special corner points smooth", mcmx C IF(MCMX.EQ.0) GO TO 115 DO 114 M=1,MCMX MM=MCT(M) II=I2G(MM) JJ=J2G(MM) DO 116 K=1,4 I=II+IH(K) J=JJ+JH(K) IF(HB(I,J)+ZB(I,J).EQ.0.) GO TO 116 HSUB(I,J)=888. 116 CONTINUE 114 CONTINUE 115 CONTINUE C IF(NSQRW.EQ.0) GO TO 124 DO 122 L=1,NSQRW IF(F1DACT(L).EQ.'T') GO TO 122 I=ISQR(L) IF(I.LE.1.OR.I.GE.IMXB1) GO TO 122 J=JSQR(L) IF(J.LE.1.OR.J.GE.JMXB1) GO TO 122 IF(HSUB(I,J).EQ.999.) GO TO 122 K=ISIDE(L) II=I+IHH(K) JJ=J+JHH(K) IF(HSUB(II,JJ).EQ.999.) GO TO 122 IF(HSUB(I,J).LE.(ZBMIN(L)+1.)) GO TO 123 IF(HSUB(II,JJ).LE.(ZBMIN(L)+1.)) GO TO 123 GO TO 122 123 HSUB(I,J)=888. HSUB(II,JJ)=888. 122 CONTINUE 124 CONTINUE M=0 MM=0 DO 130 J=1,JMXB1 DO 130 I=1,IMXB1 IF(HSUB(I,J).EQ.999.) GO TO 130 IF(HSUB(I,J).EQ.888.) GO TO 132 IF(M+MM.GE.M2G_) GO TO 140 M=M+1 I2G(M)=I J2G(M)=J GO TO 130 132 IF(MM.GE.MCT_) GO TO 140 MM=MM+1 I2GG(MM)=I J2GG(MM)=J 130 CONTINUE 140 MM2G=M MCMX=MM IF (MM2G.GE.M2G_) THEN WRITE(*,1998) M2G_,I2G(M),J2G(M) 1998 FORMAT(' TOTAL HGHT. POINTS FOR SMOOTHING REACHES MAX.',I8, 1 ' LAST I,J =',2I8) STOP ENDIF IF (MCMX.GE.MCT_) THEN WRITE(*,1997) MCT_,I2G(M),J2G(M) 1997 FORMAT(' TOTAL SPECIAL CORNER POINTS (HGHTS) REACHES MAX.',I8, 1 ' LAST I,J =',2I8) STOP ENDIF C IF(MCMX.EQ.0) GO TO 400 DO 310 M=1,MCMX I2G(M+MM2G)=I2GG(M) J2G(M+MM2G)=J2GG(M) 310 CONTINUE 400 CONTINUE C C write (*,*)' itime = ',ITIME C WRITE(*,150)MM2G,MCMX C 150 FORMAT(' TOTAL NO OF SQUARES FOR SMOOTHING = ',I10,' + ',I10) C WRITE(*,332)MCMX C 332 FORMAT(' TOTAL SPECIAL CORNER POINTS=',I7,' I AND J =') c PRINT 333,(I2GG(M),M=1,MCMX) c PRINT 333,(J2GG(M),M=1,MCMX) C 333 FORMAT (40I3) c CALL FLT2G RETURN END SUBROUTINE FLT2G C AUG 22 1984 J. CHEN TDL IBM 360/195 C PURPOSE C TO ELIMINATE 2-GRID COMPUTATIONAL NOISE NEARBY LAND-WATER C BOUNDARIES. GRIG POINTS NEARBY INTERIOR BOUNDARIES ARE C PICKED OUT FROM TESTING IN SUBROUTINE 'SMPT2G'. SMOOTHING C DOES CONSERVE MASS BY ASSUMING OUSIDE OF SELECTED POINTS IS C NON-INTERACTING WITH THE POINTS TO BE SMOOTHED. C DATA SET USE C NONE C RETURNS C NO C VARIABLES C MM2G = TOTAL NO. OF POINTS TO BE SMOOTHED C MCMX = SPECIAL SMOOTHING POINTS C I2G(),J2G() = I,J COORD'S OF GRID POINTS FOR SMOOTHING C HB( , ) = SURGE HEIGHT FILED C ZB( , ) = TERRAIN HEIGHTS C ZBM( , ) = BARRIER HEIGHTS C HCT = CUTOFF DEPTH ,1 FT EXCEPT THE SPECIAL POINTS C WHICH 0 FT IS USED C COORDINATES ARE PICKED OUT FROM SUBROUTINE 'CONTTY)'. C INCLUDE 'parm.for' C COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /GAMAF/ GAMA,GAMA1,GAMFP(4) COMMON /GPRT1/ DOLLAR,EBSN CHARACTER*2 DOLLAR CHARACTER*1 EBSN COMMON /SM2G/ MM2G,MCMX,I2G(M2G_),J2G(M2G_) INTEGER*2 MM2G,MCMX,I2G,J2G COMMON /DUMB8/ IP(4),JP(4),IH(4),JH(4) COMMON /DUM88/ IIH(4),JJH(4),IHH(4),JHH(4) DIMENSION GAMF(4) C C IF (MCMX.EQ.0) GO TO 10 DO 1010 M=1,MCMX I=I2G(MM2G+M) J=J2G(MM2G+M) 1010 HSUB(I,J)=999. 10 CONTINUE C HCT=1. M1=1 MMM=MM2G NNN=0 3000 CONTINUE NNN=NNN+1 IF (NNN.GT.2) GO TO 2000 DO 190 M=M1,MMM I=I2G(M) J=J2G(M) IF(HSUB(I,J)+ZB(I,J).LE.HCT) GO TO 190 C HSM=0. HDIF=0. GAM0=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J) C DO 170 K=1,4 II=I+IIH(K) JJ=J+JJH(K) IF (II.EQ.0.OR.JJ.EQ.0) THEN GAMF(K)=0. GOTO 160 ENDIF IF (EBSN.EQ.'$'.OR.EBSN.EQ.'+') THEN GAMFK=ELPDL2(II)+(ELPCL2(II)-ELPDL2(II))*SINL2(JJ) c GAMF(K)=.5*(1.+GAMFK/GAM0)-1. Z=GAMFK/GAM0 GAMF(K)=2.*Z/(1.+Z)-1. ELSE GAMF(K)=GAMFP(K) ENDIF c IF (II.EQ.IMXB.OR.JJ.EQ.JMXB) GOTO 160 IF (HSUB(II,JJ).EQ.999.) GO TO 160 IF (HSUB(II,JJ)+ZB(II,JJ).LE.HCT) GO TO 160 KK=MOD(K,4)+1 IA=I+IP(K) IB=I+IP(KK) JA=J+JP(K) JB=J+JP(KK) Z =ZBM(IA,JA) ZZ=ZBM(IB,JB) ZZZ=AMIN1(Z,ZZ)+HCT C WATER ON EITHER SIDE MUST EXCEED THE BARRIER(LOWER) BY AT IF(NNN.EQ.2) GO TO 200 IF( HSUB(I,J).LT.ZZZ.OR.HSUB(II,JJ).LT.ZZZ) THEN HSM=HSM+HSUB(I,J) HZZ=HSUB(I,J) ELSE HSM=HSM+HSUB(II,JJ) HZZ=HSUB(II,JJ) ENDIF GOTO 210 200 IF( HSUB(I,J).LT.ZZZ.AND.HSUB(II,JJ).LT.ZZZ) THEN HSM=HSM+HSUB(I,J) HZZ=HSUB(I,J) ELSE HSM=HSM+HSUB(II,JJ) HZZ=HSUB(II,JJ) ENDIF GOTO 210 160 CONTINUE HSM=HSM+HSUB(I,J) HZZ=HSUB(I,J) 210 CONTINUE C WEIGHTING ADJUSTMENTS IN I-DIRECTION DUE TO POLAR GRIDS. HDIF=HDIF+HZZ*GAMF(K) 170 CONTINUE C HBTEMP=(1./3.)*HSUB(I,J)+(HSM+HDIF)/6. C HBTEMP=.5*HSUB(I,J)+.125*(HSM+HDIF) HB(I,J)=AMAX1(HBTEMP,-ZB(I,J)) 190 CONTINUE IF(MCMX.EQ.0) GO TO 2000 DO 1020 M=1,MM2G I=I2G(M) J=J2G(M) 1020 HSUB(I,J)=999. HCT=0. M1=MM2G+1 MMM=MM2G+MCMX DO 1030 M=M1,MMM I=I2G(M) J=J2G(M) 1030 HSUB(I,J)=HB(I,J) GO TO 3000 2000 CONTINUE RETURN END SUBROUTINE FILTER C SEPTEMBER 1980 JYE CHEN TDL IBM 360/195 C PURPOSE C TO SMOOTH THE SURGE HEIGHTS ONCE AN HOUR TO ELIMINATE C 2-INTERVAL NOISES IN THE FIELDS. C C DATA SET USE C NONE C C VARIABLES C HSUB( , ) = SCRATCH SPACES C IP(4) JP(4) = SHIFTS OF I J FROM HEIGHT PT TO 4 MOMNTN CORNERS C IBB(6,M) = RANGES AND INCREMENTS FOR 3 DO-LOOPS OF I,J C M=1 INTERIOR,M=2 HORZNTL BNDY,M=3 VERTICAL BNDY C RANGES ON THREE BASIN SEGMENTS ARE: C C 222222222........2222222 C 311111111........1111113 C 311111111........1111113 C . . J C . . J C 311111111........1111113 J C 311111111........1111113 J C 222222222........2222222 ---->I C C J2(J) JC = SHIFT J-SUBSCRIPT TO PRESENT TIME OF SURGE FIELD C HB( , ) = SURGE FIELD C ZB( , ) = DEPTH FIELD C IIH(4) JJH(4) = I/J-SHIFT SUBSCRIPTS FOR SURGE POINTS C IP(4) JP(4) = I/J-SHIFT SUBSCRPTS FOR MOMNTM PTS, IN 'BLCKDT' C C THE (I,J) SHIFTS TO MOMENTUM POINTS, VIA IP(4) AND JP(4), C ON CORNER POINTS K=1 TO 4 ARE: C C I+0,J+1 I+1,J+1 C K=4.-----.K=3 C I I C I I*J I C I I C K=1.-----.K=2 C I+0,J+0 I+1,J+0 C C ZBM( , ) = MAX BARRIER HEIGHTS AT MOMENTUM POINTS C GENERAL COMMENTS C THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. C A 5-POINT SMOOTHING OPERATOR IS USED, NEIGHBORING GRIDS C WEIGHED CONDITIONALLY IF THEY ARE DRY OR BLOCKED C BY BARRIERS. SHIFTS OF (I,J) TO ADJACENT SQUARES ARE SET C VIA IIH(4) AND JJH(4). IIH/0,1,0,-1/,JJH/-1,0,1,0/ C C I-1 I I+1 C .-----. +HMX,HB,ZB HEIGHT PTS C I+IIH I I+0 I .ZBM BARRIER PTS C J+1 + I 3+ I C J+JJH I J+1 I EXAMPLE 0.INTERIOR POINTS C .-----.-----.-----. I 1 I I 0 I C I I-1 I I I+1 I I I I I C J I 4+ I + I 2+ I (1/8) I1 4 1I + (1/8)IA 0 BI C I J+0 I(I,J)I J+0 I I I I I C .-----.-----.-----. I 1 I I 0 I C I I+0 I C J-1 I 1+ I C I J-1 I A=GAMA1,B=GAMA C .-----. C I-1 I I+1 C C .-----. C I I C J+1 I 3+ I EXAMPLE 1. SUQARE 2 EXCLUDED C I I C .-----.-----. I 1 I I 0 I C I I I I I I I C J I 4+ I + I (1/8) I1 5 0I + (1/8)IA B 0I C I I(I,J)I I I I I C .-----.-----. I 1 I I 0 I C I I C J-1 I 1+ I C I I C .-----. C C EXAMPLE 2. SQUARES 2 AND 3 C ARE EXCLUDED C .-----.-----. I 0 I I 0 I C I I I I I I I C J I 4+ I + I (1/8) I1 6 0I + (1/8)IA B 0I C I I(I,J)I I I I I C .-----.-----. I 1 I I 0 I C I I C J-1 I 1+ I C I I C .-----. C INCLUDE 'parm.for' C COMMON /FFTH/ ITIME,MHALT COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /EGTH/ DELS,DELT,G,COR COMMON /DUMB8/ IP(4),JP(4),IH(4),JH(4) COMMON /DUM88/ IIH(4),JJH(4),IHH(4),JHH(4) COMMON /GAMAF/ GAMA,GAMA1,GAMFP(4) COMMON /GPRT1/ DOLLAR,EBSN COMMON /HTERAIN/ HTER DIMENSION GAMF(4) DIMENSION IBB(6,3) CHARACTER*2 DOLLAR CHARACTER*1 EBSN DATA HCRT/0.5/ DATA IBB(1,1),IBB(1,2),IBB(1,3)/2,1,2/, 1 IBB(3,1), IBB(3,3)/1, 1/, 2 IBB(4,1),IBB(4,2),IBB(4,3)/2,1,1/, 3 IBB(6,1),IBB(6,2) /1,1 / C write (*,*) "Inside filter" IBB(2,1)=JMXB2 IBB(2,3)=JMXB2 IBB(3,2)=JMXB2 IBB(2,2)=JMXB1 IBB(5,1)=IMXB2 IBB(6,3)=IMXB2 IBB(5,2)=IMXB1 IBB(5,3)=IMXB1 DO 100 J=1,JMXB1 DO 100 I=1,IMXB1 100 HSUB(I,J)=HB(I,J) DO 210 M=1,3 JL =IBB(1,M) JM =IBB(2,M) JNCR=IBB(3,M) IL =IBB(4,M) IM =IBB(5,M) INCR=IBB(6,M) DO 200 J=JL,JM,JNCR DO 190 I=IL,IM,INCR C SKIP COMPUTATIONS FOR TERRAIN HIGHER THAN 35 FT. C CHANGED BY NSM TO HTER FROM 35 11/20/2010 : Accepted 4/4/2011 IF(ZB(I,J).LE.-HTER) GO TO 190 HST=HSUB(I,J)+ZB(I,J) C TEST FOR LAND WETTED LESS THAN 1 FT. IF(HST.LT.HCRT) GO TO 190 HSM=0. HDIF=0. C C******** 'LOOP' FOR SMOOTHING****************************************** GAM0=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J) DO 170 K=1,4 II=I+IIH(K) JJ=J+JJH(K) IF (M.EQ.1) GOTO 158 IF (II.EQ.0.OR.JJ.EQ.0) THEN GAMF(K)=0. GO TO 160 ENDIF 158 CONTINUE IF (EBSN.EQ.'$'.OR.EBSN.EQ.'+') THEN GAMFK=ELPDL2(II)+(ELPCL2(II)-ELPDL2(II))*SINL2(JJ) C GAMF(K)=.5*(1.+GAMFK/GAM0)-1. Z=GAMFK/GAM0 GAMF(K)=2.*Z/(1.+Z)-1. ELSE GAMF(K)=GAMFP(K) ENDIF IF (II.EQ.IMXB.OR.JJ.EQ.JMXB) GOTO 160 C THE NEIGHBORING SQUARE IS MERELY WET(LESS THAN 1 FT). IF (HSUB(II,JJ)+ZB(II,JJ).LT.HCRT) GO TO 160 KK=MOD(K,4)+1 IA=I+IP(K) IB=I+IP(KK) JA=J+JP(K) JB=J+JP(KK) Z =ZBM(IA,JA) ZZ=ZBM(IB,JB) ZZZ=AMIN1(Z,ZZ)+1. !Huiqing.Liu on June 2014 !C AAT Modified on 4/5/2011: So it is 0.5 feet (ie what HCRT is) ! ZZZ=AMIN1(Z,ZZ)+HCRT C WATER ON EITHER SIDE MUST EXCEED THE BARRIER(LOWER) BY AT C LEAST 1 FT. IF( HSUB(I,J).LT.ZZZ.OR.HSUB(II,JJ).LT.ZZZ) THEN HSM=HSM+HSUB(I,J) HZZ=HSUB(I,J) ELSE HSM=HSM+HSUB(II,JJ) HZZ=HSUB(II,JJ) ENDIF GOTO 212 160 CONTINUE HSM=HSM+HSUB(I,J) HZZ=HSUB(I,J) 212 CONTINUE C WEIGHTING ADJUSTMENTS IN I-DIRECTION DUE TO POLAR GRIDS. HDIF=HDIF+HZZ*GAMF(K) C 170 CONTINUE HTEMP=.5*HSUB(I,J)+(HSM+HDIF)/8. C******* END 'LOOP' FOR SMOOTHING ************************************** C 180 HB(I,J)=AMAX1(-ZB(I,J),HTEMP) 190 CONTINUE 200 CONTINUE 210 CONTINUE RETURN END SUBROUTINE FLTER2 INCLUDE 'parm.for' C COMMON /FFTH/ ITIME,MHALT COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /EGTH/ DELS,DELT,G,COR COMMON /DUMB8/ IP(4),JP(4),IH(4),JH(4) COMMON /DUM88/ IIH(4),JJH(4),IHH(4),JHH(4) COMMON /GAMAF/ GAMA,GAMA1,GAMFP(4) COMMON /GPRT1/ DOLLAR,EBSN COMMON /HTERAIN/ HTER DIMENSION GAMF(4) CHARACTER*2 DOLLAR CHARACTER*1 EBSN DATA HCRT/0.5/ C DO 100 J=1,JMXB DO 100 I=1,IMXB 100 HSUB(I,J)=HB(I,J) DO 210 M=1,2 C TOP AND BOTTOM BOUNDARY CONDITIONS M=2 C INTERIOR POINTS M=1 INCLUDING SIDE BOUNDARIES. IF (M.EQ.1) THEN IL =2 IM =IMXB2 INCR=1 ELSE IL =1 IM =IMXB1 INCR=IMXB2 ENDIF DO 200 J=1,JMXB1 DO 190 I=IL,IM,INCR C SKIP COMPUTATIONS FOR TERRAIN HIGHER THAN 35 FT. C CHANGED BY NSM TO HTER FROM 35 11/20/2010 : Accepted 4/4/2011 IF(ZB(I,J).LE.-HTER) GO TO 190 HST=HSUB(I,J)+ZB(I,J) C TEST FOR LAND WETTED LESS THAN 1 FT. IF(HST.LT.HCRT) GO TO 190 HSM=0. HDIF=0. C C******** 'LOOP' FOR SMOOTHING****************************************** GAM0=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J) DO 170 K=1,4 II=I+IIH(K) JJ=J+JJH(K) IF (JJ.EQ.0) JJ=JMXB1 IF (II.EQ.0) THEN GAMF(K)=0. GO TO 160 ENDIF 158 CONTINUE IF (EBSN.EQ.'$'.OR.EBSN.EQ.'+') THEN GAMFK=ELPDL2(II)+(ELPCL2(II)-ELPDL2(II))*SINL2(JJ) c GAMF(K)=.5*(1.+GAMFK/GAM0)-1. Z=GAMFK/GAM0 GAMF(K)=2.*Z/(1.+Z)-1. ELSE GAMF(K)=GAMFP(K) ENDIF IF (II.EQ.IMXB) GOTO 160 C THE NEIGHBORING SQUARE IS MERELY WET(LESS THAN 1 FT). IF (HSUB(II,JJ)+ZB(II,JJ).LT.HCRT) GO TO 160 KK=MOD(K,4)+1 IA=I+IP(K) IB=I+IP(KK) JA=J+JP(K) JB=J+JP(KK) Z =ZBM(IA,JA) ZZ=ZBM(IB,JB) C ZZZ=AMIN1(Z,ZZ)+1. ZZZ=AMIN1(Z,ZZ)+HCRT C WATER ON EITHER SIDE MUST EXCEED THE BARRIER(LOWER) BY AT C LEAST 1 FT. IF( HSUB(I,J).LT.ZZZ.OR.HSUB(II,JJ).LT.ZZZ) THEN HSM=HSM+HSUB(I,J) HZZ=HSUB(I,J) ELSE HSM=HSM+HSUB(II,JJ) HZZ=HSUB(II,JJ) ENDIF GOTO 212 160 CONTINUE HSM=HSM+HSUB(I,J) HZZ=HSUB(I,J) 212 CONTINUE C WEIGHTING ADJUSTMENTS IN I-DIRECTION DUE TO POLAR GRIDS. HDIF=HDIF+HZZ*GAMF(K) C 170 CONTINUE HTEMP=.5*HSUB(I,J)+(HSM+HDIF)/8. C******* END 'LOOP' FOR SMOOTHING ************************************** C 180 HB(I,J)=AMAX1(-ZB(I,J),HTEMP) 190 CONTINUE 200 CONTINUE 210 CONTINUE DO 300 I=1,IMXB 300 HB(I,JMXB)=HB(I,1) RETURN END SUBROUTINE SNAPSH C CONVERT SURGE ENVELOP DATA WITH 999 FOR DRY SQUARES SO THAT ZB() C CAN BE FREE FROM THE MEMORY. INCLUDE 'parm.for' C COMMON /IDENT/ AIDENT(40),DACLOK(7) COMMON /GPRT1/ DOLLAR,EBSN COMMON /SWTCH/ IOPERL(5) COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 CHARACTER*2 DOLLAR CHARACTER*1 EBSN C REPLACE HMX BY 999 IF THE SQUARE IS DRY. DO 1 J=1,JMXB DO 1 I=1,IMXB IF (HB(I,J)+ZB(I,J).EQ.0) THEN IHMX(I,J)=999 ELSE IHMX(I,J)=HB(I,J)*10.+.5 ENDIF 1 CONTINUE C C WRITE(77) IMXB,JMXB C WRITE(77) AIDENT C WRITE (77) ((IHMX(I,J),I=1,IMXB),J=1,JMXB) C IF (DOLLAR.EQ.'2$') THEN DO 20 J=1,JMXB DO 20 I=1,IMXB IF (IHMX(I,J).EQ.999) GOTO 20 IHMX(I,J)=(HMX(I,J)-ZSUB(J))*10.+.5 20 CONTINUE ENDIF C IF(IOPERL(1).NE.2) THEN NGRP=(JMXB1-1)/25+1 DO 1130 NN=1,NGRP J1=1+(NN-1)*25 J2=J1+MIN0(24,JMXB1-J1 ) CC WRITE(*,180)(J,J=J1,J2) 180 FORMAT(/I6,26I5) C DO 1120 I=1,IMXB1 CC WRITE(*,'(1H ,I2,26I5)') I,(IHMX(I,J),J=J1,J2) 1120 CONTINUE 1130 CONTINUE ENDIF C RETURN END !Huiqing Liu/MDL implment a new subroutine xfilter written by Craig A. Mattocks SUBROUTINE XFILTER(VERSION) C SEPTEMBER 1980 JYE CHEN TDL IBM 360/195 C PURPOSE C TO SMOOTH THE SURGE HEIGHTS ONCE AN HOUR TO ELIMINATE C 2-INTERVAL NOISES IN THE FIELDS. C C VERSION=1 means use original version for backward compat. C VERSION=2 means use current correct version. C C DATA SET USE C NONE C C VARIABLES C HSUB( , ) = SCRATCH SPACES C IP(4) JP(4) = SHIFTS OF I J FROM HEIGHT PT TO 4 MOMNTN CORNERS C IBB(6,M) = RANGES AND INCREMENTS FOR 3 DO-LOOPS OF I,J C M = 1 INTERIOR,M = 2 HORZNTL BNDY,M = 3 VERTICAL BNDY C RANGES ON THREE BASIN SEGMENTS ARE: C C 222222222........2222222 C 311111111........1111113 C 311111111........1111113 C . . J C . . J C 311111111........1111113 J C 311111111........1111113 J C 222222222........2222222 ---->I C C J2(J) JC = SHIFT J-SUBSCRIPT TO PRESENT TIME OF SURGE FIELD C HB( , ) = SURGE FIELD C ZB( , ) = DEPTH FIELD C IIH(4) JJH(4) = I/J-SHIFT SUBSCRIPTS FOR SURGE POINTS C IP(4) JP(4) = I/J-SHIFT SUBSCRPTS FOR MOMNTM PTS, IN 'BLCKDT' C C THE (I,J) SHIFTS TO MOMENTUM POINTS, VIA IP(4) AND JP(4), C ON CORNER POINTS K = 1 TO 4 ARE: C C I+0,J+1 I+1,J+1 C K = 4.-----.K = 3 C I I C I I*J I C I I C K = 1.-----.K = 2 C I+0,J+0 I+1,J+0 C C ZBM( , ) = MAX BARRIER HEIGHTS AT MOMENTUM POINTS C GENERAL COMMENTS C THIS SUBROUTINE RESIDES IN OVERLAY 'CMPUTE'. C A 5-POINT SMOOTHING OPERATOR IS USED, NEIGHBORING GRIDS C WEIGHED CONDITIONALLY IF THEY ARE DRY OR BLOCKED C BY BARRIERS. SHIFTS OF (I,J) TO ADJACENT SQUARES ARE SET C VIA IIH(4) AND JJH(4). IIH/0,1,0,-1/,JJH/-1,0,1,0/ C C I-1 I I+1 C .-----. +HMX,HB,ZB HEIGHT PTS C I+IIH I I+0 I .ZBM BARRIER PTS C J+1 + I 3+ I C J+JJH I J+1 I EXAMPLE 0.INTERIOR POINTS C .-----.-----.-----. I 1 I I 0 I C I I-1 I I I+1 I I I I I C J I 4+ I + I 2+ I (1/8) I1 4 1I + (1/8)IA 0 BI C I J+0 I(I,J)I J+0 I I I I I C .-----.-----.-----. I 1 I I 0 I C I I+0 I C J-1 I 1+ I C I J-1 I A = GAMA1,B = GAMA C .-----. C I-1 I I+1 C C .-----. C I I C J+1 I 3+ I EXAMPLE 1. SUQARE 2 EXCLUDED C I I C .-----.-----. I 1 I I 0 I C I I I I I I I C J I 4+ I + I (1/8) I1 5 0I + (1/8)IA B 0I C I I(I,J)I I I I I C .-----.-----. I 1 I I 0 I C I I C J-1 I 1+ I C I I C .-----. C C EXAMPLE 2. SQUARES 2 AND 3 C ARE EXCLUDED C .-----.-----. I 0 I I 0 I C I I I I I I I C J I 4+ I + I (1/8) I1 6 0I + (1/8)IA B 0I C I I(I,J)I I I I I C .-----.-----. I 1 I I 0 I C I I C J-1 I 1+ I C I I C .-----. C C ====================================================================== C C Development History C ------------------- C C Author Date Purpose C ------ ----- --------- C Jye Chen/TDL September 1980 Wrote the SLOSH filter code C A. Taylor and C. Forbes April 2011 Modified the SLOSH filter code C to include the cross Laplacian C filter C C. Forbes June 2011 Incorporated masks according C to Killworth (1991) and C Deleersnijder (1995), C implemented masks for coast- C lines, barriers, boundaries, C and included averaging. C Recoded the DO loop and GO TO C logic to make code more efficient. C A. Taylor August 2011 Allowed user to choose different C HCRT values. (0.5 or 0.1) C C References C ---------- C C Killworth, P.D., Stainforth, D., Webb, D.J., Paerson, S.M., 1991. C The development of a free-surface Bryan-Cox-Semtner ocean model. C Journal of Physical Oceanography, 21, 1333-1348. C C Deleersnijder, E., Campin, J.M., 1995. On the computation of the C barotropic mode of a free-surface world ocean model. C Annales Geophysicae, 13, 675-688. C C ====================================================================== INCLUDE 'parm.for' C IMPLICIT REAL(A-H, O-Z), INTEGER(I-N) PARAMETER (NP = 4) COMMON /FFTH/ ITIME, MHALT COMMON /DUMB3/ IMXB, JMXB, IMXB1, JMXB1, IMXB2, JMXB2 COMMON /EGTH/ DELS, DELT, G, COR COMMON /DUMB8/ IP(NP), JP(NP), IH(NP), JH(NP) COMMON /DUM88/ IIH(NP), JJH(NP), IHH(NP), JHH(NP) COMMON /GAMAF/ GAMA, GAMA1, GAMFP(NP) COMMON /GPRT1/ DOLLAR, EBSN COMMON /HTERAIN/ HTER COMMON /SMTH/ ISMTH C COMMON /MASS/ AATMAX, AATMIN, AATSUM, AATCNT REAL GAMF(NP) REAL GAMFX(NP) INTEGER IBB(NP+2, NP-1) INTEGER IIX(NP), JJX(NP) INTEGER IP1(NP), JP1(NP), IP2(NP), JP2(NP) INTEGER MASKII(NP), MASKIX(NP), MBARRX(NP), MBARRI(NP) INTEGER MBOUNI(NP), MBOUNX(NP) REAL EXCLUDE CHARACTER*2 DOLLAR CHARACTER*1 EBSN C DATA HCRT /0.1/ ! Water surface height difference that determines C DATA HCRT /0.5/ ! Water surface height difference that determines C DATA HCRT /1.0/ ! whether water spills over a barrier DATA IBB(1,1), IBB(1,2), IBB(1,3) /2, 1, 2/, 1 IBB(3,1), IBB(3,3) /1, 1/, 2 IBB(4,1), IBB(4,2), IBB(4,3) /2, 1, 1/, 3 IBB(6,1), IBB(6,2) /1, 1 / DATA IIX /-1, 1, 1, -1/ DATA JJX /-1, -1, 1, 1/ DATA IP1 /-1, 1, 2, 0/ DATA JP1 / 0, -1, 1, 2/ DATA IP2 / 0, 2, 1, -1/ DATA JP2 /-1, 0, 2, 1/ IF (ISMTH.EQ.2) THEN HCRT = 0.5 EXCLUDE = HCRT ELSE IF (ISMTH.EQ.3) THEN HCRT = 0.1 EXCLUDE = HCRT ELSE IF (ISMTH.EQ.4) THEN HCRT = 0.2 EXCLUDE = HCRT ELSE IF (ISMTH.EQ.5) THEN HCRT = 0.3 EXCLUDE = HCRT ELSE IF (ISMTH.EQ.6) THEN HCRT = 0.4 EXCLUDE = HCRT ELSE IF (ISMTH.EQ.7) THEN C HCRT = 0.1 HCRT = 1.E-5 EXCLUDE = 1.E-5 ELSE IF (ISMTH.EQ.8) THEN HCRT = 0.2 EXCLUDE = 1.E-5 ELSE IF (ISMTH.EQ.9) THEN HCRT = 0.3 EXCLUDE = 1.E-5 ELSE IF (ISMTH.EQ.10) THEN HCRT = 0.4 EXCLUDE = 1.E-5 ELSE IF (ISMTH.EQ.11) THEN HCRT = 0.5 EXCLUDE = 1.E-5 ENDIF IBB(2,1) = JMXB2 IBB(2,3) = JMXB2 IBB(3,2) = JMXB2 IBB(2,2) = JMXB1 IBB(5,1) = IMXB2 IBB(6,3) = IMXB2 IBB(5,2) = IMXB1 IBB(5,3) = IMXB1 C -------------------------- C Check on mass conservation C -------------------------- HSUMHB = 0. DO J = 1,JMXB1 DO I = 1,IMXB1 HSUB(I,J) = HB(I,J) HSUMHB = HSUMHB + HB(I,J) END DO END DO C -------------- C Filter weights C -------------- ALFA = 1./8. ! Weight of the combination Laplacian filter BETA = 1. ! Weight of the 'x' Laplacian filter C ---------------------------------------- C Start Main Loop for all (i,j) grid cells C M = 1 interior cells C M = 2,3 boundary cells {AT??} C ---------------------------------------- DO 100 M = 1, 3 JL = IBB(1,M) JM = IBB(2,M) JNCR = IBB(3,M) IL = IBB(4,M) IM = IBB(5,M) INCR = IBB(6,M) DO 90 J = JL, JM, JNCR DO 80 I = IL, IM, INCR C ------------------------------------------ C Skip computations for terrain > 35 ft. C Changed from a hard-wired constant to HTER C by NSM on 11/20/2010, accepted 4/4/2011 C ------------------------------------------ IF (ZB(I,J) .LE. -HTER) GO TO 80 C -------------------------------- C Test for land wetted, but < EXCLUDE C -------------------------------- HST = HSUB(I,J) + ZB(I,J) IF (HST .LT. EXCLUDE) GO TO 80 C --------------------------- C Initialize filter variables C --------------------------- HSM2 = 0. ! Sum of all + cells for + filter XSM2 = 0. ! Sum of all x cells for x filter MASKC = 1 ! Mask for center cell C HDIF = 0. C XDIF = 0. C GAM0 = ELPDL2(I) + (ELPCL2(I) - ELPDL2(I)) * SINL2(J) C ---------------------------------------------------------- C Start smoothing loop around the center cell C counterclockwise from: C south , east , north and west for + filter C and from: C southwest, southeast, northeast and northwest for x filter C ---------------------------------------------------------- DO 30 K = 1, 4 C |---------|---------|---------| C | | | | C | IX | II | IX | C | JX | JJ | JX | C | K = 4 | K = 3 | K = 3 | C |---------|---------|---------| C | | | | C | II | I | II | C | JJ | J | JJ | C | K = 4 | | K = 2 | C |---------|---------|---------| C | | | | C | IX | II | IX | C | JX | JJ | JX | C | K = 1 | K = 1 | K = 2 | C |---------|---------|---------| C ------------------------- C Initialize mask variables C ------------------------- MASKII(K) = 0 ! Mask for II cells according to height (+ filter) (0 no flow, 1 flow) MBARRI(K) = 0 ! Mask for II cells according to barrier (+ filter) (1 no flow, 0 flow) MBOUNI(K) = 0 ! Mask for II cells according to boundary (+ filter) (1 no flow, 0 flow) MASKIX(K) = 0 ! Mask for IX cells according to height (x filter) (0 no flow, 1 flow) MBARRX(K) = 0 ! Mask for IX cells according to barrier (x filter) (1 no flow, 0 flow) MBOUNX(K) = 0 ! Mask for IX cells according to boundary (x filter) (1 no flow, 0 flow) C -------------------------- C Start '+' Laplacian filter C -------------------------- II = I + IIH(K) JJ = J + JJH(K) IX = I + IIX(K) JX = J + JJX(K) C ---------------------------------------------- C If any variable from filter hits the boundary, C do not do anything C ---------------------------------------------- IF (II .EQ. 0 .OR. JJ .EQ. 0 .OR. & IX .EQ. 0 .OR. JX .EQ. 0 .OR. & I .EQ. 0 .OR. J .EQ. 0 ) GO TO 30 CC IF (M .NE. 1) THEN CC IF (II .EQ. 0 .OR. JJ .EQ. 0) THEN CC MASKII(K) = 0 CC MBARRI(K) = 0 CC MBOUNI(K) = 1 CC GO TO 10 CC END IF CC END IF C -------------------------------------------------------- C If the cell from + filter hits the IMXB or JMXB boundary C or the height < critical height (HCRT) set masks values C accordingly C -------------------------------------------------------- IF (II .EQ. IMXB .OR. JJ .EQ. JMXB .OR. & HSUB(II,JJ) + ZB(II,JJ) .LT. HCRT) THEN CCC already set to 0 MASKII(K) = 0 CCC already set to 0 MBARRI(K) = 0 MBOUNI(K) = 1 ELSE C -------------------------------------------- C Water on either side must exceed the (lower) C barrier by at least HCRT C -------------------------------------------- KK = MOD(K,4) + 1 IA = I + IP(K) IB = I + IP(KK) JA = J + JP(K) JB = J + JP(KK) Z = ZBM(IA,JA) ZZ = ZBM(IB,JB) ZZZ = AMIN1(Z,ZZ) + HCRT C ?? IF (HSUB(I,J) .LT. ZZZ .OR. HSUB(II,JJ) .LT. ZZZ) IF (HSUB(I,J) .LT. ZZZ .AND. HSUB(II,JJ) .LT. ZZZ) & THEN CCC already set to 0 MASKII(K) = 0 MBARRI(K) = 1 ELSE MASKII(K) = 1 CCC already set to 0 MBARRI(K) = 0 END IF END IF 10 CONTINUE C -------------------------- C Start 'x' Laplacian filter C -------------------------- KK = MOD(K, 4) + 1 KM1 = MOD(K + 2, 4) + 1 IA = I + IP(K) IB = I + IP(KK) JA = J + JP(K) JB = J + JP(KK) IA2 = I + IP2(K) JA2 = J + JP2(K) IA1 = I + IP1(K) JA1 = J + JP1(K) IAM1 = I + IP(KM1) JAM1 = J + JP(KM1) C ------------------------------------------------------ C If any variable from the 'x' filter hits the boundary, C do not do anything C ------------------------------------------------------ CC IF (M .NE. 1) THEN CC IF (IX .EQ. 0 .OR. JX .EQ. 0) THEN CC MASKIX(K) = 0 CC MBARRX(K) = 0 CC MBOUNX(K) = 1 CC GO TO 20 CC END IF CC END IF C -------------------------------------------------------- C If the cell from x filter hits the IMXB or JMXB boundary C or the height < critical height (HCRT) set masks values C accordingly C -------------------------------------------------------- IF (IX .EQ. IMXB .OR. JX .EQ. JMXB .OR. & HSUB(IX,JX) + ZB(IX,JX) .LT. HCRT) THEN CCC already set to 0 MASKIX(K) = 0 CCC already set to 0 MBARRX(K) = 0 MBOUNX(K) = 1 ELSE C ------------------------------------------------- C Water on either side must exceed the (lower) C barrier by at least HCRT. C Water should be able to flow from IJ to II C and from II to IX, clockwise or counter-clockwise C ------------------------------------------------- Z = ZBM(IA,JA) ZZ = ZBM(IB,JB) ZZZ = AMIN1(Z,ZZ) + HCRT YY = ZBM(IA2,JA2) YYY = AMIN1(Z,YY) + HCRT XX = ZBM(IA1,JA1) XXX = AMIN1(Z,XX) + HCRT WW = ZBM(IAM1,JAM1) WWW = AMIN1(Z,WW) + HCRT IF (((HSUB(I ,J ) .GE. ZZZ .AND. & HSUB(II,JJ) .GE. YYY) .OR. & (HSUB(IX,JX) .GE. YYY .AND. & HSUB(II,JJ) .GE. ZZZ)) & .OR. & ((HSUB(I ,J ) .GE. WWW .AND. & HSUB(IA1,JA1) .GE. XXX) .OR. & (HSUB(IX ,JX ) .GE. XXX .AND. & HSUB(IA1,JA1) .GE. WWW)) ) & THEN CCC un-used variable XSM = XSM + HSUB(IX,JX) CCC un-used variable XZZ = HSUB(IX,JX) MASKIX(K) = 1 CCC already set to 0 MBARRX(K) = 0 ELSE CCC un-used variable XSM = XSM + HSUB(I,J) CCC un-used variable XZZ = HSUB(I,J) CCC already set to 0 MASKIX(K) = 0 MBARRX(K) = 1 END IF END IF 20 CONTINUE 30 CONTINUE C -------------------------------------------------------- C Start combination of '+' and 'x' Laplacian filters: C H = H + (alpha * (plusFilter - beta * crossFilter)) * dt C -------------------------------------------------------- DO 40 KM = 1, 4 KKM = MOD(KM , 4) + 1 KM1M = MOD(KM + 2, 4) + 1 IX = I + IIX(KM) JX = J + JJX(KM) II = I + IIH(KM) JJ = J + JJH(KM) CC IF (M .EQ. 1) GOTO 158 CC IF (II.EQ.0.OR.JJ.EQ.0) THEN CC GAMF(K)=0. CC GO TO 160 CC ENDIF CC 158 CONTINUE CC IF (EBSN.EQ.'$'.OR.EBSN.EQ.'+') THEN CC GAMFK=ELPDL2(II)+(ELPCL2(II)-ELPDL2(II))*SINL2(JJ) CC GAMFKX=ELPDL2(IX)+(ELPCL2(IX)-ELPDL2(IX))*SINL2(JX) C GAMF(K)=.5*(1.+GAMFK/GAM0)-1. CC Z=GAMFK/GAM0 CC GAMF(K)=2.*Z/(1.+Z)-1. CC ZX=GAMFKX/GAM0 CC GAMFX(K)=2.*ZX/(1.+ZX)-1. CC ELSE CC GAMF(K)=GAMFP(K) CC GAMFX(K)=GAMFP(K) CC ENDIF IX2 = I + IIX(KKM) JX2 = J + JJX(KKM) II2 = I + IIH(KKM) JJ2 = J + JJH(KKM) IX3 = I + IIX(KM1M) JX3 = J + JJX(KM1M) II3 = I + IIH(KM1M) JJ3 = J + JJH(KM1M) MSUMII = MASKIX(KM ) + MASKIX(KKM) + MASKC MSUMIX = MASKIX(KM ) + MASKIX(KM1M) + MASKC C --------------------------------------------- C If any cell of the filter is at the boundary, C do not do any combination or averaging C --------------------------------------------- IF (II .EQ. 0 .OR. JJ .EQ. 0 .OR. & IX .EQ. 0 .OR. JX .EQ. 0 .OR. & I .EQ. 0 .OR. J .EQ. 0 .OR. & IX2 .EQ. 0 .OR. JX2 .EQ. 0 .OR. & IX3 .EQ. 0 .OR. JX3 .EQ. 0 .OR. & II2 .EQ. 0 .OR. JJ2 .EQ. 0 .OR. & II3 .EQ. 0 .OR. JJ3 .EQ. 0 ) GO TO 40 C ---------------------------------------------------- C If a cell in the '+' the filter is dry (e.g., II is C coastline and IJ, IX, and IX2 are water cells), then C average the two closest cells (IX and IX2) and C the center cell IJ, multiplying by their masks and C weighting them with the addition of all masks, C so that if any of the 3 are dry, they do not get C averaged. C C |---------|--------|---------| C | | | | C |---------|--------|---------| C | | IJ | | C |---------|--------|---------| C | IX | II | IX2 | C | JX | JJ | JX2 | C |---------|--------|---------| C C ---------------------------------------------------- IF (MASKII(KM) .EQ. 0 .AND. MSUMII .NE. 0 .AND. & MBARRI(KM) .EQ. 0 .AND. MBOUNX(KM) .EQ. 0 .AND. & MBOUNI(KM) .EQ. 0) THEN write (*,*) "TROUBLE: HSUB is supposed to be read only here" write (*,*) "because it is needed for calc in other cells." HSUB(II,JJ) = (MASKIX(KM )*HSUB(IX ,JX) + & MASKIX(KKM)*HSUB(IX2,JX2) + & MASKC *HSUB(I ,J) ) / & (MASKIX(KM ) + MASKIX(KKM) + MASKC) END IF C ---------------------------------------------------- C If a cell in the 'x' the filter is dry (e.g., IX is C coastline and IJ, II, and II3 are water cells), then C average the two closest cells (II3 and II) and C the center cell IJ, multiplying by their masks and C weighting them with the addition of all masks, C so that if any of the 3 are dry, they do not get C averaged. C C |---------|--------|---------| C | | | | C | | | | C |---------|--------|---------| C | II3 | I | | C | JJ3 | J | | C |---------|--------|---------| C | IX | II | | C | JX | JJ | | C |---------|--------|---------| C C ---------------------------------------------------- IF (MASKIX(KM) .EQ. 0 .AND. MSUMIX .NE. 0 .AND. & MBARRX(KM) .EQ. 0 .AND. MBOUNX(KM) .EQ. 0 .AND. & MBOUNI(KM) .EQ. 0. ) & THEN write (*,*) "TROUBLE: HSUB is supposed to be read only here" write (*,*) "because it is needed for calc in other cells." HSUB(IX,JX) = (MASKII(KM )*HSUB(II ,JJ ) + & MASKIX(KM1M)*HSUB(II3,JJ3) + & MASKC *HSUB(I ,J ) ) / & (MASKIX(KM ) + MASKIX(KM1M) + MASKC) END IF C ------------------------------------------------------ C Calculate the two terms of the combination '+' and 'x' C filter: C HSM2 = + filter: SUM (II,JJ) from K = 1, 4 C XSM2 = x filter: SUM (IX,JX) from K = 1, 4 C C |---------|---------|---------| C | | | | C | IX | II | IX | C | JX | JJ | JX | C | K = 4 | K = 3 | K = 3 | C |---------|---------|---------| C | | | | C | II | I | II | C | JJ | J | JJ | C | K = 4 | | K = 2 | C |---------|---------|---------| C | | | | C | IX | II | IX | C | JX | JJ | JX | C | K = 1 | K = 1 | K = 2 | C |---------|---------|---------| C C ------------------------------------------------------ C ??? MASKC is always 1. HSM2 = HSM2 + & MASKC * MASKII(KM) * HSUB(II,JJ) - & MASKC * MASKII(KM) * HSUB( I, J) C ------------------------------------------------------ C A B C Capital letters are interior, lower are exterior. C a D e Killsworth (91) recommended that: C f g h (1) If a cell 'a' is exterior for adjacent calculations to C cell 'D' it be replaced with the average of any C interior cells that it 'a' was adjacent to (ie 'D' and C diagonals to 'D' that touch 'a' (i.e. 'A')) C (2) If a cell 'f' is exterior and diagonal to 'D' it is C ignored. C C A B C Capital letters are interior, lower are exterior. C a D e adj= ALFA * (HSM2 - BETA*(XSM2)) (assume Beta = 1) C f g h adj=ALFA ((B-D)+(a-D)+(g-D)+(e-D) C -1/2[(f-D)+(h-D)+(C-D)+(A-D)] C Now (f-D), and (h-D) can be ignored since they are diagonals C (g-D) can be replaced with (D-D) = 0 C => adj=ALFA ((B-D)+ (a-D) + (e-D) -1/2 [(C-D) + (A-D)] C (a-D) becomes ((A + D)/2 - D) = (A - D)/2 C (e-D) becomes ((C + D)/2 - D) = (C - D)/2 C => adj=ALFA ((B-D)+ (A-D)/2 + (C-D)/2 - (C-D)/2 - (A-D)/2 C => adj=ALFA (B-D) C C In looking at this, the A and C terms are dropped from the C XSM component. However if 'e' had been 'E', we would have C kept the C term in the XSM. So a diagonal cell is included C in XSM if it is interior (connected to point in question) C and the adjacent cells to it are interior. C C One permutation that Killsworth didn't prepare for is: C A B C Capital letters are interior, lower are exterior. C E D F adj= ALFA * (HSM2 - BETA*(XSM2)) (assume Beta = 1) C G a H adj=ALFA ((B-D)+(E-D)+(F-D)+(a-D) C -1/2[(A-D)+(C-D)+(G-D)+(H-D)] C (a-D) becomes ((G + D + H)/3 - D) = (G - D)/3 + (H - D)/3 C => adj=ALFA ((B-D)+(E-D)+(F-D)+(G-D)/3+(H-D)/3 C -1/2[(A-D)+(C-D)+(G-D)+(H-D)] C In this case 0 .ne. (G-D)/3 -(G-D)/2 +(H-D)/3 -(H-D)/2 C I would argue that they should cancel, so XSM doesn't need C to take this into consideration. C ------------------------------------------------------ if (VERSION.eq.1) then C Use original bad formulation for backward compatibility with hch2. XSM2 = XSM2 + 0.5 * & ((MASKIX(KM) * MASKII(KM) * MASKC * MASKIX(KKM)) & * HSUB(IX,JX) - & (MASKIX(KM) * MASKII(KM) * MASKC * MASKIX(KKM)) & * HSUB( I, J)) else XSM2 = XSM2 + 0.5 * & MASKC * MASKIX(KM)*MASKII(KM)*MASKII(KM1M)* & (HSUB(IX,JX) - HSUB(I, J)) endif cc if (MASKII(KM).eq.0) then cc XSM2 = XSM2 + 1/6 * cc & MASKC * MASKIX(KM)*MASKIX(KKM)* cc & (HSUB(IX,JX) - HSUB(I, J)) cc endif cc if (MASKII(KM1M).eq.0) then cc XSM2 = XSM2 + 1/6 * cc & MASKC * MASKIX(KM)*MASKIX(KM1M)* cc & (HSUB(IX,JX) - HSUB(I, J)) cc endif ccc if ((i.eq.148).and.(j.eq.69)) then ccc if (MASKIX(KM)*MASKII(KM)*MASKII(KM1M).eq.1) then ccc write (*,*) "148,69,x",KM,HSUB(IX,JX),ZB(IX,JX) ccc write (*,*) XSM2, HSM2, HSUB(I,J) ccc end if ccc if (MASKII(KM).eq.1) then ccc write (*,*) "148,69,p",KM,HSUB(II,JJ),ZB(II,JJ) ccc write (*,*) XSM2, HSM2, HSUB(I,J) ccc end if ccc end if 40 CONTINUE C -------------------------------------------------------- C AT: Original plus filter was: C HTEMP=.5*HSUB(I,J)+(HSM+HDIF)/8. C =HSUB(I,J)+ 1/8*(HSM+HDIF - 4*HSUB(I,J)) C =HSUB(I,J)+ ALFA*(HSM+HDIF - 4*HSUB(I,J)) C C Introducing cross lapacian resulted in: C HTEMP=HSUB(I,J)+ ALFA*((HSM+HDIF-4*HSUB(I,J)) - C 1/2*((XSM+XDIF)-4*HSUB(I,J))) C C Now HSM2 and XSM2 are calculated in such a way that C HSM2 is the sum of HSM - 4*HSUB C XSM2 is the sum of 0.5 * (XSM - 4*HSUB) C So... C C Calculate the combination '+' and 'x' Laplacian filter: C H = H + (alpha * (plusFilter - beta * crossFilter)) * dt C -------------------------------------------------------- HTEMP = HSUB(I,J) + ALFA * (HSM2 - BETA*(XSM2)) C --------------------------------------- C Assign the value to the height variable C --------------------------------------- HB(I,J) = AMAX1(-ZB(I,J), HTEMP) 80 CONTINUE 90 CONTINUE 100 CONTINUE cc HSUMHB2 = 0. cc DO J = 1,JMXB1 cc DO I = 1,IMXB1 cc HSUMHB2 = HSUMHB2 + HB(I,J) cc END DO cc END DO cc AAT = HSUMHB - HSUMHB2 cc if (AATMAX.lt.AAT) then cc AATMAX=AAT cc end if cc if (AATMIN.gt.AAT) then cc AATMIN=AAT cc end if cc AATSUM = AATSUM + abs(AAT) cc AATCNT = AATCNT + 1 cc write (*,*) "mass?", AATCNT, AATMAX, AATMIN, AATSUM, AAT, HSUMHB, cc & IMXB1 * JMXB1 RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE INPUTN2 ! AUTHOR: BRIAN ZACHRY, JAN 2014 ! ! PURPOSE: TO ALLOW FOR A NEW TRACK FILE FORMAT THAT ALLOWS ! FOR MORE OPTIONS AND TRACK LENGTH UP TO 999 HOURS ! ! IT IS COMPATIBLE WITH THE ORIGINAL TRACK FILE ! FORMAT WHICH IS THE DEFAULT FLAG COMMON /STRMPS/ X(999),Y(999),P(999),R(999),DIR(999),SP(999) COMMON /SPLN/ ALT(15),ALN(15),AX(15),AY(15),RL(15),ANGD(15), 1 PT(15),RT(15),XLAT(999),YLONG(999),RLNGTH COMMON /IDENT/ AIDENT(40),DACLOK(7) C STIME interferes with C code, so switched to STIME2 COMMON /STIME2/ ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND COMMON /DATUM/ SEADTM,DTMLAK COMMON /DATUM1/ DTMCHN COMMON /XOKE/ XOKE CHARACTER*1 XOKE COMMON /FLES/ FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1 CHARACTER*256 FLE5,FLE9,FLE8,FLE91,FLE99,FLE10,FLE20,FLE30,FLE1 C CHARACTER*80 dummy COMMON /LANDFL/ LFTIME CHARACTER*80 LFTIME ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !New track file Declarations !Declare integer parms INTEGER i,jj,iError,iLine,Flag_TrkFile,iAdjustHour INTEGER iDateHr,iLandfallHr,iBegHr,iEndHr,iTrkLen INTEGER iInterp,TrkFileLen,TrkHeadLen,iHour,interpTimes INTEGER Flag_Press,iNumDataHead !Declare character parms CHARACTER(LEN=100) chLine,chDum,chStorm,chDateTime,chFileVer CHARACTER(LEN=100) chAuthor,chExtraInfo CHARACTER(LEN=100) Flag_Interp !Declare character arrays CHARACTER(LEN=100) TrkFileStr(1100) !Declare real parms REAL fOceanDat,fLakeDat,fCanalDat,fWindAdj,fRMWAdj REAL fLat,fLon,fDeltaP,fRMW,fPress,fAmbPress REAL yLat,yLon,yDeltaP,yRMW !Declare integer arrays INTEGER iarrHour(999) INTEGER hrlyHour(999) !Declare real arrays REAL hrlyLat(999),hrlyLon(999),hrlyRMW(999),hrlyDeltaP(999) REAL farrLat(999),farrLon(999),farrDeltaP(999),farrRMW(999) REAL y2Lat(999),y2Lon(999),y2DeltaP(999),y2RMW(999) !New common blocks with track file info COMMON /TRKHRS/ ITRACKLEN COMMON /TRKSTR/ TRKFILESTR,TRKHEADLEN,TRKFILELEN cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Code to read in the new 2013 track file format ! This format provides more flexability for the user ! Up to 999 total hours (included interpolated) are allowed ! ! Documentation for the new format can be found ... ! ! Date: 03/22/2012 ! Author: Brian Zachry !Open the track file OPEN(5,FILE=FLE5,STATUS='OLD',ACTION='READ',IOSTAT=iError) IF (iError /= 0) THEN WRITE(*,*) 'ERROR: Cannot Find the Track File' WRITE(*,'(A,A)') ' Track File Name: ',FLE5 STOP END IF !Set default track type flag to old track file format !--New track format = 2013 !--Old track format = 1992 Flag_TrkFile = 1992 !Check first line of track file for version number (year) READ(5,'(A100)',IOSTAT=iError) chLine IF (chLine(1:12) == 'FileVersion=') THEN READ(chLine(13:len_trim(chLine)),*) Flag_TrkFile ELSE IF (chLine(1:13) == ' FileVersion=') THEN READ(chLine(14:len_trim(chLine)),*) Flag_TrkFile ELSE IF (chLine(1:14) == ' FileVersion =') THEN READ(chLine(15:len_trim(chLine)),*) Flag_TrkFile END IF !Write out track file version to screen WRITE(*,'(A,I0)') 'track file ver: ', Flag_TrkFile ! CALL FLUSH (6) ! Huiqing Liu /MDL flush the screen output !Rewind track file after reading first line REWIND(UNIT=5) ! ---------------------- !Main code block of how to processes the track file !Process new SLOSH track file format (2013) IF (Flag_TrkFile == 2013) THEN !Initialize the arrays iarrHour(:) = 0 hrlyHour(:) = 0 hrlyLat(:) = 0.0 hrlyLon(:) = 0.0 hrlyRMW(:) = 0.0 hrlyDeltaP(:)= 0.0 farrLat(:) = 0.0 farrLon(:) = 0.0 farrDeltaP(:)= 0.0 farrRMW(:) = 0.0 y2Lat(:) = 0.0 y2Lon(:) = 0.0 y2DeltaP(:) = 0.0 y2RMW(:) = 0.0 !Set defaults for mandatory variables iBegHr = -999 iEndHr = -999 iLandfallHr = -999 iDateHr = -999 fOceanDat = -999.9 fLakeDat = -999.9 chDateTime = ' ' !***NEED TO ADD CHECK FOR THIS *** !Set defaults for optional variables chAuthor = ' ' chExtraInfo = ' ' chStorm = ' ' fWindAdj = 1.0 fRMWAdj = 1.0 Flag_Interp = 'Linear' Flag_Press = 0 fAmbPress = -999.9 fCanalDat = -999.9 ! ---------------------- !Code block to find length of header block and store !the entire track file as a string (A100) !--Store track file as string in array TrkFileStr !--Store header file lengh in TrkHeadLen !--Store track file length in TrkFileLen TrkHeadLen = -999 iLine = 0 iNumDataHead = 0 DO READ(5,'(A100)',IOSTAT=iError) chLine IF (iError /= 0) EXIT iLine = iLine + 1 TrkFileStr(iLine) = chLine !Various string triggers for the end of the header block IF (chLine(1:31) == 'Hour,Lat,Lon,DeltaP(mb),RMW(mi)' .OR. 1 chLine(1:33) == 'Hour,Lat,Lon,Pressure(mb),RMW(mi)' .OR. 2 chLine(1:5) == 'Hour,' .OR. 3 chLine(1:6) == ' Hour,' .OR. 4 chLine(1:6) == 'Hour ,' .OR. 5 chLine(1:7) == ' Hour ,') THEN iNumDataHead = iNumDataHead + 1 !Count for header TrkHeadLen = iLine !Set the length of header block END IF END DO TrkFileLen = iLine !Set the length of the track file !Check to make sure the data block (one header) was found IF (TrkHeadLen <= 0) THEN WRITE(*,*) 'ERROR: Data Block Not Found In Track File' WRITE(*,*) ' Insert Data Block Header: "Hour,Lat,Lon,..."' ELSE IF (iNumDataHead > 1) THEN WRITE(*,*) 'ERROR: More Than One Data Block Header Found' WRITE(*,*) ' Triggered By Variations On "Hour,"' STOP END IF !Rewind file to go back through header block and parse data REWIND(UNIT=5) ! ---------------------- !Code block to read and parse the header block DO jj = 1,TrkHeadLen READ(5,'(A100)',IOSTAT=iError) chLine IF (iError /= 0) EXIT !Ignore comments, denoted by '#', provide check for extra spaces IF (chLine(1:1) == '#' .OR. 1 chLine(1:2) == ' #' .OR. 2 chLine(1:3) == ' #' .OR. 3 chLine(1:4) == ' #' .OR. 4 chLine(1:5) == ' #') THEN CYCLE END IF !Required - Version number (year) of track file IF (chLine(1:12) == 'FileVersion=') THEN chFileVer = chLine(13:100) WRITE(*,*) 'FileVersion=',chFileVer(1:len_trim(chFileVer)) !Optional - Author of track file ELSE IF (chLine(1:7) == 'Author=') THEN chAuthor = chLine(8:100) WRITE(*,*) 'Author=',chAuthor(1:len_trim(chAuthor)) !Optional - Any extra info author wants to add ELSE IF (chLine(1:10) == 'ExtraInfo=') THEN chExtraInfo = chLine(11:100) WRITE(*,*) 'ExtraInfo=',chExtraInfo(1:len_trim(chExtraInfo)) !Optional - Name of storm in track file ELSE IF (chLine(1:10) == 'StormName=') THEN chStorm = chLine(11:100) WRITE(*,*) 'StormName=',chStorm(1:len_trim(chStorm)) !Required - Date and time of reference point in track file ELSE IF (chLine(1:9) == 'DateTime=') THEN chDateTime = chLine(10:100) WRITE(*,*) 'DateTime=',chDateTime(1:len_trim(chDateTime)) !Required - Hour of track file associated with DateTime variable ELSE IF (chLine(1:9) == 'DateHour=') THEN READ(chLine(10:len_trim(chLine)),*) iDateHr WRITE(*,'(A,I0)') ' DateHour=',iDateHr !Required - Landfall hour or nearest approach NAP ELSE IF (chLine(1:13) == 'LandfallHour=') THEN READ(chLine(14:len_trim(chLine)),*) iLandfallHr WRITE(*,'(A,I0)') ' LandfallHour=',iLandfallHr !Required - Begin hour of simulation ELSE IF (chLine(1:10) == 'BeginHour=') THEN READ(chLine(11:len_trim(chLine)),*) iBegHr WRITE(*,'(A,I0)') ' BeginHour=',iBegHr !Required - End hour of simulation ELSE IF (chLine(1:8) == 'EndHour=') THEN READ(chLine(9:len_trim(chLine)),*) iEndHr WRITE(*,'(A,I0)') ' EndHour=',iEndHr !Required (for PressureType=1) - Environmental pressure ELSE IF (chLine(1:16) == 'AmbientPressure=') THEN READ(chLine(17:len_trim(chLine)),*) fAmbPress WRITE(*,'(A,F0.1)') ' AmbientPressure=',fAmbPress !Required - Ocean datum ELSE IF (chLine(1:11) == 'OceanDatum=') THEN READ(chLine(12:len_trim(chLine)),*) fOceanDat WRITE(*,'(A,F0.1)') ' OceanDatum=',fOceanDat !Required - Lake datum ELSE IF (chLine(1:10) == 'LakeDatum=') THEN READ(chLine(11:len_trim(chLine)),*) fLakeDat WRITE(*,'(A,F0.1)') ' LakeDatum=',fLakeDat !Optional (for XOKE) - Canal datum for OKE basin ELSE IF (chLine(1:11) == 'CanalDatum=') THEN READ(chLine(12:len_trim(chLine)),*) fCanalDat WRITE(*,'(A,F0.1)') ' CanalDatum=',fCanalDat !Optional - Wind (pressure) adjustment percentage ELSE IF (chLine(1:11) == 'WindFactor=') THEN READ(chLine(12:len_trim(chLine)),*) fWindAdj WRITE(*,'(A,F0.1)') ' WindFactor=',fWindAdj !Optional - RMW adjustment percentage ELSE IF (chLine(1:10) == 'RMWFactor=') THEN READ(chLine(11:len_trim(chLine)),*) fRMWAdj WRITE(*,'(A,F0.1)') ' RMWFactor=',fRMWAdj !Optional - Interpolation type to use on track file data ELSE IF (chLine(1:18) == 'InterpolationType=') THEN Flag_Interp = chLine(19:100) WRITE(*,*) 'InterpolationType=', 1 Flag_Interp(1:len_trim(Flag_Interp)) !Optional - Format of input pressure data ELSE IF (chLine(1:13) == 'PressureType=') THEN READ(chLine(14:len_trim(chLine)),*) Flag_Press WRITE(*,'(A,I0)') ' PressureType=',Flag_Press END IF END DO !Loop/read of the header block !Perform rough sanity checks for the key variables !Check begin hour of simulation (0-999) IF (iBegHr < 0 .OR. iBegHr > 999) THEN WRITE(*,*) 'ERROR: Begin Hour Out Of Bounds' WRITE(*,'(A,I0)') ' Begin Hour = ',iBegHr STOP !Check end hour of simulation (0-999) ELSE IF (iEndHr < 0 .OR. iEndHr > 999) THEN WRITE(*,*) 'ERROR: End Hour Out Of Bounds' WRITE(*,'(A,I0)') ' End Hour = ',iEndHr STOP !Check hour of reference date in track file (0-999) ELSE IF (iDateHr < 0 .OR. iDateHr > 999) THEN WRITE(*,*) 'ERROR: Date Hour Out Of Bounds' WRITE(*,'(A,I0)') ' Date Hour = ',iDateHr STOP !Check landfall hour or nearest approach hour NAP (0-999) ELSE IF (iLandfallHr < 0 .OR. iLandfallHr > 999) THEN WRITE(*,*) 'ERROR: Landfall Hour Out Of Bounds' WRITE(*,'(A,I0)') ' Landfall Hour = ',iLandfallHr STOP !Check ambient pressure, must be provided for Flag_Press=1 ELSE IF (Flag_Press == 1 .AND. 1 (fAmbPress < 900. .OR. fAmbPress > 1100.)) THEN WRITE(*,*) 'ERROR: Ambient Pressure Out Of Bounds' WRITE(*,'(A,F0.3)') ' Ambient Pressure = ',fAmbPress WRITE(*,*) ' "AmbientPressure=" Required For "PressureType=1"' STOP !Check ocean datum (-50-50) ELSE IF (fOceanDat < -50. .OR. fOceanDat > 50.) THEN WRITE(*,*) 'ERROR: Ocean Datum Out Of Bounds' WRITE(*,'(A,F0.3)') ' Ocean Datum = ',fOceanDat STOP !Check lake datum (-50-50) ELSE IF (fLakeDat < -50. .OR. fLakeDat > 50.) THEN WRITE(*,*) 'ERROR: Lake Datum Out Of Bounds' WRITE(*,'(A,F0.3)') ' Lake Datum = ',fLakeDat STOP !Check simulation length, begin and end hours ELSE IF (iBegHr >= iEndHr) THEN WRITE(*,*) 'ERROR: Begin Hour >= End Hour' WRITE(*,'(A,I0)') ' Begin Hour = ',iBegHr WRITE(*,'(A,I0)') ' End Hour = ',iEndHr STOP !Check landfall hour compared to begin and end hours ELSE IF (iLandfallHr < iBegHr .OR. iLandfallHr > iEndHr) THEN WRITE(*,*) 'ERROR: Landfall Hour Outside Simulation Hours' WRITE(*,'(A,I0)') ' Landfall Hour = ',iLandfallHr WRITE(*,'(A,I0)') ' Begin Hour = ',iBegHr WRITE(*,'(A,I0)') ' End Hour = ',iEndHr STOP !Check reference hour compared to begin and end hours ELSE IF (iDateHr < iBegHr .OR. iDateHr > iEndHr) THEN WRITE(*,*) 'ERROR: Date Hour Outside Simulation Hours' WRITE(*,'(A,I0)') ' Date Hour = ',iDateHr WRITE(*,'(A,I0)') ' Begin Hour = ',iBegHr WRITE(*,'(A,I0)') ' End Hour = ',iEndHr STOP END IF ! ---------------------- !Code block to read data block and store values in array iTrkLen = 0 iAdjustHour = 0 DO !Read for the typical DeltaP format (PressureType=0) IF (Flag_Press == 0) THEN READ(5,*,IOSTAT=iError) iHour,fLat,fLon,fDeltaP,fRMW !Read for central pressure format (PressureType=1) ELSE IF (Flag_Press == 1) THEN READ(5,*,IOSTAT=iError) iHour,fLat,fLon,fPress,fRMW fDeltaP = fAmbPress - fPress !Read for central pressure format with ambient (PressureType=2) ELSE IF (Flag_Press == 2) THEN READ(5,*,IOSTAT=iError) iHour,fLat,fLon,fPress,fRMW,fAmbPress fDeltaP = fAmbPress - fPress !Error for user has provided a value not supported ELSE WRITE(*,*) 'ERROR: Flag_Press Not Valid (0,1,2)' WRITE(*,'(A,I0)') ' Flag_Press: ',Flag_Press WRITE(*,*) ' 0: Hour,Lat,Lon,DeltaP,RMW' WRITE(*,*) ' 1: Hour,Lat,Lon,Press,RMW' WRITE(*,*) ' *1*: Requirement In Header: AmbientPressure=' WRITE(*,*) ' 2: Hour,Lat,Lon,Press,RMW,AmbPress' STOP END IF IF (iError /= 0) EXIT !Perform rough sanity checks !Check to make sure track hour is in bounds (0-999) IF (iHour < 0 .OR. iHour > 999) THEN WRITE(*,*) 'ERROR: Track Hour Out Of Bounds (0-999)' WRITE(*,'(A,I0)') ' Hour: ',iHour STOP !Check to make sure latitude is in bounds (0-90 N) ELSE IF (fLat < 0. .OR. fLat > 90.) THEN WRITE(*,*) 'ERROR: Latitude Out Of Bounds (0-90 N)' WRITE(*,'(A,f0.3)') ' Latitude: ',fLat WRITE(*,'(A,I0)') ' At Hour: ',iHour STOP !Check to make sure longitude is in bounds (0-180 W) ELSE IF (fLon > 0. .OR. fLon < -180. ) THEN WRITE(*,*) 'ERROR: Longitude Out Of Bounds (0-180 W)' WRITE(*,'(A,f0.3)') ' Longitude: ',fLon WRITE(*,'(A,I0)') ' At Hour: ',iHour STOP !Check to make sure RMW is in bounds (1-500) ELSE IF (fRMW < 1. .OR. fRMW > 500.) THEN WRITE(*,*) 'ERROR: RMW Out Of Bounds (1-500)' WRITE(*,'(A,f0.3)') ' RMW: ',fRMW WRITE(*,'(A,I0)') ' At Hour: ',iHour STOP !Check to make sure delta pressure is in bounds (0.1-213) ELSE IF (fDeltaP > 213. .OR. fDeltaP < 0.1) THEN WRITE(*,*) 'ERROR: Delta Pressure Out Of Bounds (0.1-213)' WRITE(*,'(A,f0.3)') ' Delta Pressure: ',fDeltaP WRITE(*,'(A,I0)') ' At Hour: ',iHour STOP !Check to make sure pressure is in bounds (800 - 1100) ELSE IF ((Flag_Press == 1 .OR. Flag_Press == 2) .AND. 1 (fPress < 800. .OR. fPress > 1100.)) THEN WRITE(*,*) 'ERROR: Pressure Out Of Bounds (800-1100)' WRITE(*,'(A,f0.3)') ' Pressure: ',fPress WRITE(*,'(A,I0)') ' At Hour: ',iHour STOP !Check to make sure ambient pressure is in bounds (900 - 1100) ELSE IF ((Flag_Press == 2) .AND. 1 (fAmbPress < 900. .OR. fAmbPress > 1100.)) THEN WRITE(*,*) 'ERROR: Ambient Pressure Out Of Bounds (900-1100)' WRITE(*,'(A,f0.3)') ' Ambient Pressure: ',fAmbPress WRITE(*,'(A,I0)') ' At Hour: ',iHour STOP END IF !Keep track of the line number in the track file data block iTrkLen = iTrkLen + 1 !Compute adjustment to track hours to start at hour 1 !This is needed for the SLOSH arrays to work correctly IF (iTrkLen == 1) THEN iAdjustHour = 1 - iHour END IF !Adjust each hour based on offset from hour 1 iHour = iHour + iAdjustHour !Make the track file arrays iarrHour(iTrkLen) = iHour farrLat(iTrkLen) = fLat farrLon(iTrkLen) = fLon * (-1.0) farrDeltaP(iTrkLen) = fDeltaP farrRMW(iTrkLen) = fRMW !Perform simple checks for the hour increments in file IF (iTrkLen > 1) THEN !Check to make sure hours are increasing IF (iarrHour(iTrkLen) < iarrHour(iTrkLen-1)) THEN WRITE(*,*) 'ERROR: Track Hours Are In Decreasing Order' WRITE(*,*) ' Previous Hour: ', 1 iarrHour(iTrkLen-1) - iAdjustHour WRITE(*,*) ' Current Hour: ', 1 iarrHour(iTrkLen) - iAdjustHour STOP !Check for duplicate hours in track file ELSE IF (iarrHour(iTrkLen) - iarrHour(iTrkLen-1) == 0) THEN WRITE(*,*) 'ERROR: Duplicate Hours In Track File' WRITE(*,*) ' Previous Hour: ', 1 iarrHour(iTrkLen-1) - iAdjustHour WRITE(*,*) ' Current Hour: ', 1 iarrHour(iTrkLen) - iAdjustHour STOP END IF END IF END DO !Read of the data block ! ---------------------- !Code block to set the houly data arrays needed to run SLOSH !Interpolate track file data to hourly (if necessary) !Currently supports 'Linear' and 'Spline' interpolations !Perform linear interpolation IF (Flag_Interp == 'Linear' .OR. Flag_Interp == 'linear') THEN !Loop through the track file and interpolate where necessary DO i = 1,iTrkLen !Set the index for the array based on the hour iHour = iarrHour(i) !Hourly data (line) exists in the track file hrlyHour(iHour) = iarrHour(i) hrlyLat(iHour) = farrLat(i) hrlyLon(iHour) = farrLon(i) hrlyDeltaP(iHour) = farrDeltaP(i) hrlyRMW(iHour) = farrRMW(i) !Interpolation is necessary IF ((iarrHour(i+1) - iarrHour(i)) > 1 .AND. i < iTrkLen) THEN !Find the number of interpolations to hourly interpTimes = iarrHour(i+1) - iarrHour(i) !Interpolate to hourly DO jj = 1,interpTimes-1 iHour = iarrHour(i) + jj fFrac = REAL(jj) / REAL(interpTimes) !Set hourly data based on interpolation hrlyHour(iHour) = iarrHour(i) + 1 (fFrac*(iarrHour(i+1)-iarrHour(i))) hrlyLat(iHour) = farrLat(i) + 1 (fFrac*(farrLat(i+1)-farrLat(i))) hrlyLon(iHour) = farrLon(i) + 1 (fFrac*(farrLon(i+1)-farrLon(i))) hrlyDeltaP(iHour)= farrDeltaP(i) + 1 (fFrac*(farrDeltaP(i+1)-farrDeltaP(i))) hrlyRMW(iHour) = farrRMW(i) + 1 (fFrac*(farrRMW(i+1)-farrRMW(i))) END DO END IF END DO !Perform cubic spline interpolation ELSEIF (Flag_Interp == 'Spline' .OR. Flag_Interp == 'spline') THEN !Calculate cubic functions CALL spline(REAL(iarrHour(1:iTrkLen)),farrLat(1:iTrkLen), 1 iTrklen,1.0E+30,1.0E+30,y2Lat) CALL spline(REAL(iarrHour(1:iTrkLen)),farrLon(1:iTrkLen), 1 iTrkLen,1.0E+30,1.0E+30,y2Lon) CALL spline(REAL(iarrHour(1:iTrkLen)),farrDeltaP(1:iTrkLen), 1 iTrkLen,1.0E+30,1.0E+30,y2DeltaP) CALL spline(REAL(iarrHour(1:iTrkLen)),farrRMW(1:iTrkLen), 1 iTrkLen,1.0E+30,1.0E+30,y2RMW) !Loop through the track file and interpolate where necessary DO i = 1,iTrkLen !Set the index for the array based on the hour iHour = iarrHour(i) !Hourly data (line) exists in the track file hrlyHour(iHour) = iarrHour(i) hrlyLat(iHour) = farrLat(i) hrlyLon(iHour) = farrLon(i) hrlyDeltaP(iHour) = farrDeltaP(i) hrlyRMW(iHour) = farrRMW(i) !Interpolation is necessary IF ((iarrHour(i+1) - iarrHour(i)) > 1 .AND. i < iTrkLen) THEN !Find the number of interpolations to hourly interpTimes = iarrHour(i+1) - iarrHour(i) !Interpolate to hourly DO jj = 1,interpTimes-1 iHour = iarrHour(i) + jj !Set hourly data based on interpolation hrlyHour(iHour) = iHour !Latitude CALL splint(REAL(iarrHour(1:iTrkLen)),farrLat(1:iTrkLen), 1 y2Lat(1:iTrkLen),iTrkLen,REAL(iHour),yLat) hrlyLat(iHour) = yLat !Longtiude CALL splint(REAL(iarrHour(1:iTrkLen)),farrLon(1:iTrkLen), 1 y2Lon(1:iTrkLen),iTrkLen,REAL(iHour),yLon) hrlyLon(iHour) = yLon !DeltaP CALL splint(REAL(iarrHour(1:iTrkLen)), 1 farrDeltaP(1:iTrkLen),y2DeltaP(1:iTrkLen), 2 iTrkLen,REAL(iHour),yDeltaP) hrlyDeltaP(iHour) = yDeltaP !RMW CALL splint(REAL(iarrHour(1:iTrkLen)),farrRMW(1:iTrkLen), 1 y2RMW(1:iTrkLen),iTrkLen,REAL(iHour),yRMW) hrlyRMW(iHour) = yRMW END DO END IF END DO !Error user has provided an interpolation not supported ELSE WRITE(*,*) 'ERROR: Interpolation Type Not Specified' WRITE(*,*) ' Flag_Interp:',Flag_Interp(1:len_trim(Flag_Interp)) WRITE(*,*) ' "Linear": Linear Interpolation, Default' WRITE(*,*) ' "Spline": Cubic Spline Interpolation' STOP END IF !Special canal datum for OKE, must be > 0.0 to implement XOKE flag IF (fCanalDat > 0.) THEN XOKE = 'X' DTMCHN = fCanalDat END IF !Set the ocean and lake datums SEADTM = fOceanDat DTMLAK = fLakeDat !Set the start, end, and landfall hour IBGNT = iBegHr + iAdjustHour ITEND = iEndHr + iAdjustHour JHR = iLandfallHr + iAdjustHour LFTIME = 'HR1200 27 AUG 2013' !***NEEDS TO COME FROM HEADER*** !Set the track file length for common block iTrackLen = iHour !Set the arrays needed for SLOSH to run !Perform any factor adjustments here DO i = 1,iTrackLen XLAT(i) = hrlyLat(i) YLONG(i) = hrlyLon(i) P(i) = hrlyDeltaP(i) * fWindAdj R(i) = hrlyRMW(i) * fRMWAdj !WRITE(*,*) i,hrlyHour(i),XLAT(i),YLONG(i),P(i),R(i) END DO !Close the track file CLOSE(5) ! ---------------------------------------------------------- !Read the old track file format ELSE IF (Flag_TrkFile == 1992) THEN !Set the track file length for common block iTrackLen = 100 C READ IN 2 TITLE CARDS READ (5,'(20a4)') AIDENT C READ 100 STRM PSTNS IN LAT AND LONG, MM PRESSURE DROPS, C RADII OF MAX WINDS IN ST MILES, ALL 1 HOURS APART DO 110 I=1,100 READ (5,300) ITM,XLAT(I),YLONG(I),SPeed,DIRr,P(I),R(I) 300 FORMAT(15X,I5,8F8.2) 110 CONTINUE READ (5,'(3I3)') IBGNT,ITEND,JHR READ (5,'(A80)') LFTIME READ (5,'(2f5.1,A1,F5.1)') SEADTM,DTMLAK,XOKE,DTMCHN CLOSE (5) !Track file flag is not supported (user input error) ELSE WRITE(*,*) 'ERROR: Track File Version Not Supported' WRITE(*,*) ' User Input FileVersion: ',Flag_TrkFile WRITE(*,*) ' Supported Versions: 1992 (Default) or 2013' WRITE(*,*) ' 1992 - Original Track File Format' WRITE(*,*) ' 2013 - New Track File Format' CLOSE (5) STOP END IF C ht1 < 99.9 implies init water was for tide + anomaly (so tide) C ht1 = 99.9 implies init water was missing (so surge) C 150 < ht1 or -250 > ht1 implies anomaly can be found by C mod ((ht1 + 50), 100) - 50) C The exception would be 999.9, but that shouldn't be used anymore C and I don't believe that was in any .trk files. C IF (INT(SEADTM * 10 + .5) == 999) THEN C ht1 = 99.9 implies init water was missing (so surge) SEADTM = 0 DTMLAK = 0 ENDIF IF ((SEADTM.GE.150).OR.(SEADTM.LE.-250)) THEN SEADTM = MOD ((SEADTM + 50), 100.) - 50 C Then round to nearest 10th of a foot. C SEADTM = INT(SEADTM * 10 + .5)/10. ENDIF c write (*,*) ' sea datum, and lake datum =', seadtm, dtmlak c read (*,*) seadtm,dtmlak ! call flush (6) ! flush the screen output Huiqing Liu/MDL call tflush ! flush the screen output Huiqing Liu/MDL RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE spline(X,Y,N,YP1,YPN,Y2) c SPLINE use: given an 1D array of X data and an array of Y data, c both of length N, this routine computes the 2nd derivatives, Y2 at c each X data point. The user needs to specify the values of YP1 c and YP2, which flags how the Y2 are computed at the edges. For c natural spline fitting (recommended), set YP1 and YPN to numbers c greater than 1.0E+30. c this routine called once, prior to using routine SPLINT, as a set c up for using routine SPLINT, which performs the actual c interpolation c IMPORTANT NOTE: the X data values in array X must be in ascending c order or the interpolation will fail !IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NMAX=999) DIMENSION X(N),Y(N),Y2(N),U(NMAX) c if YP1>1.0E+30 use natural spline, otherwise estimate Y2 at the c first point IF (YP1.GT..99D30) THEN Y2(1)=0. U(1)=0. ELSE Y2(1)=-0.5 U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1) ENDIF c store intermediate values of terms in the expansion series DO 11 I=2,N-1 SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1)) P=SIG*Y2(I-1)+2. Y2(I)=(SIG-1.)/P U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) * /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 11 CONTINUE c if YPN>1.0E+30 use natural spline, otherwise estimate Y2 at the c last point point IF (YPN.GT..99D30) THEN QN=0. UN=0. ELSE QN=0.5 UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1))) ENDIF Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.) c compute the Y2 from the 2nd order expansion series DO 12 K=N-1,1,-1 Y2(K)=Y2(K)*Y2(K+1)+U(K) 12 CONTINUE RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE splint(XA,YA,Y2A,N,X,Y) c SPLINT use: given an 1D array of XA data, an array of YA data, and c an array of the 2nd derivatives Y2A, all of length N, this routine c performs cubic spline interpolation, returning the interpolated c value Y at the user input value X. The Y2A are computed in c routine SPLINE, which is called once before calling SPLINT. c IMPORTANT NOTE: the X data values in array X must be in ascending c order or the interpolation will fail ! IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XA(N),YA(N),Y2A(N) KLO=1 KHI=N c determine the indices of array XA that bracket the input X value 1 IF (KHI-KLO.GT.1) THEN K=(KHI+KLO)/2 IF(XA(K).GT.X)THEN KHI=K ELSE KLO=K ENDIF GOTO 1 ENDIF c determine the finite difference along the X dimension H=XA(KHI)-XA(KLO) IF (H.EQ.0.) STOP 'Bad XA input in routine SPLINE.' c interpolate A=(XA(KHI)-X)/H B=(X-XA(KLO))/H Y=A*YA(KLO)+B*YA(KHI)+ * ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6. RETURN END !*****************************************************************************/ ! * PCGWR --- Huiqing Liu / MDL ! * ! * DESCRIPTION ! * ! * Added brand new open boundary other than the static height boundary (Clamped) ! * The new open boundary called Partial Clamped Gravity Wave Radiation boundary ! * ! * ! * HISTORY ! * 01/20/2014 Written by --- Huiqing Liu / MDL ! * ! * NOTES ! * ! * A.F. Blumberg and L.H. Kantha, Open Boundary Condition for Circulation Models ! * Journal of Hydraulic Engineering, ASCE, (111)237-255,1985. ! * ! * ! *****************************************************************************/ SUBROUTINE PCGWR (TF_HR) C INCLUDE 'parm.for' C PARAMETER (NBCPTS=12000) COMMON /BCPTS/ TIDESH(NBCPTS),BOUN_NEST(NBCPTS), !Added by Huiqing Liu /MDL Nov. 2013 * NBCPT,ISH(NBCPTS),JSH(NBCPTS) COMMON /DUMMY4/ S(800),C(800),P(800),DELP(800) COMMON /STRMSB/ C1,C2,C21,C22,AX,AY,PTENCY,RTENCY C STIME interferes with C code, so switched to STIME2 COMMON /STIME2/ ISTM,JHR,ITMADV,NHRAD,IBGNT,ITEND COMMON /DUMB3/ IMXB,JMXB,IMXB1,JMXB1,IMXB2,JMXB2 COMMON /FRST/ X1(50),X12(50) COMMON /DATUM/ SEADTM,DTMLAK cc Added New Common Block COMMON /DUMB4/ X1B(10) COMMON /EGTH/ DELS,DELT,G,COR COMMON /POLAR/ DEGREE,RMOUTH,AZMTH,DELA COMMON /INDX/ IRB,IRE,NDB,NDE C CC Tf=TF_HR*3600. IF (TF_HR.EQ.0)Tf=1.*3600. DO 100 N=1,NBCPT I=ISH(N) J=JSH(N) XR=ELPCL(I)*COSL(J) YR=ELPDL(I)*SINL(J) X=XR-C1-AX Y=YR-C2-AY RSQ=X*X+Y*Y R1=SQRT(RSQ)/5280.+1. K=R1 R2=K DR=R1-R2 K=MIN0(K,790) HB_KN=DELP(K)+DR*(DELP(K+1)-DELP(K))+SEADTM+TIDESH(N) c HB_KN=TIDESH(N) c_opn=sqrt(Zb(i,j)*G) IF(I.NE.IMXB-1.AND.J.EQ.JMXB-1)THEN Z=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J) delr=sqrt(Z)*DELA FCT1=1.-c_opn*DELT/delr-DELT/Tf IF(FCT1.LT.0)FCT1=0. FCT2=c_opn*delt/delr FCT3=DELT/Tf HB(I,J)=FCT1*HB(I,J)+FCT2*HB(I,J-1)+FCT3*HB_KN ELSEIF(I.NE.1.AND.J.EQ.1)THEN Z=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J) delr=sqrt(Z)*DELA FCT1=1.-c_opn*DELT/delr-DELT/Tf IF(FCT1.LT.0)FCT1=0. FCT2=c_opn*delt/delr FCT3=DELT/Tf HB(I,J)=FCT1*HB(I,J)+FCT2*HB(I,J+1)+FCT3*HB_KN c HB(I,J)=HB_KN ELSEIF(J.NE.JMXB-1.AND.I.EQ.IMXB-1)THEN Z=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J) Z1=ELPDL2(I-1)+(ELPCL2(I-1)-ELPDL2(I-1))*SINL2(J) delr=abs(sqrt(Z)-sqrt(Z1)) FCT1=1.-c_opn*DELT/delr-DELT/Tf IF(FCT1.LT.0)FCT1=0. FCT2=c_opn*delt/delr FCT3=DELT/Tf HB(I,J)=FCT1*HB(I,J)+FCT2*HB(I-1,J)+FCT3*HB_KN ELSEIF(J.EQ.JMXB-1.AND.I.EQ.IMXB-1)THEN Z=ELPDL2(I)+(ELPCL2(I)-ELPDL2(I))*SINL2(J) Z1=ELPDL2(I-1)+(ELPCL2(I-1)-ELPDL2(I-1))*SINL2(J) delr=abs(sqrt(Z)-sqrt(Z1)) FCT1=1.-c_opn*DELT/delr-DELT/Tf IF(FCT1.LT.0)FCT1=0. FCT2=c_opn*delt/delr FCT3=DELT/Tf HB(I,J)=FCT1*HB(I,J)+FCT2*HB(I-1,J)+FCT3*HB_KN ELSE HB(I,J) = HB_KN ENDIF 100 CONTINUE C C RETURN END