SUBROUTINE SCHWID (XMJD,RLAT,RLON,TSC,TSCLP,TSCSD,ISDAT) c c Ordering of the Schwiderski Constituents c c (1) = M2 (PRINCIPAL LUNAR, SEMIDIURNAL) c (2) = S2 (PRINCIPAL SOLAR, " ) c (3) = K1 (DECLINATION LUNI-SOLAR, DIURNAL) c (4) = O1 (PRINCIPAL LUNAR, DIURNAL) c (5) = N2 (ELLIPTICAL LUNAR, SEMIDIURNAL) c (6) = P1 (PRINCIPAL SOLAR, DIURNAL) c (7) = K2 (DECLINATION LUNI-SOLAR, SEMIDIURNAL) c (8) = Q1 (ELLIPTICAL LUNAR, DIURNAL) c (9) = MF (FORTNIGHTLY LUNAR, LONG PERIOD) c (10) = MM (MONTHLY LUNAR, LONG PERIOD) c (11) = SSA (SEMIANNUAL SOLAR, LONG PERIOD) c LOGICAL ISDAT,FIRST DOUBLE PRECISION UTC,JAN1 INTEGER YEAR COMMON /TEST/ YEAR,FIRST,FACT(11),TICON(11) DIMENSION MDY(5),XCOSP(11),XSINP(11) SAVE DATA FIRST /.TRUE./ C UTC = DBLE(XMJD) + 2400000.5D0 IF (FIRST) THEN CALL KALDAY (MDY,SECS,UTC) YEAR = MDY(3) MDY(1) = 1 MDY(2) = 1 MDY(4) = 0 MDY(5) = 0 CALL JULDAY (MDY,0.,JAN1) ENDIF C ISEC = (UTC-JAN1)*86400.D0 CALL SCTIDE (RLAT,RLON,ISEC,Y,ISDAT,XCOSP,XSINP) IF (ISDAT) THEN TSC = Y/10. TSCLP = (TICON(9) + TICON(10) + TICON(11))/10. TSCSD = TSC - TSCLP ELSE TSCLP = 0. TSCSD = 0. TSC = 0. ENDIF C RETURN END SUBROUTINE SCTIDE(X,Y,TIME,TIDEC,ISDATA,XCOSP,XSINP) * * PURPOSE: COMPUTE OCEAN TIDE (IN MM.) FOR POINT X,Y AT INPUT TIME * * SOURCE: SCHWIDERSKI M2 TIDE TAPE * * PROGRAMMER: EDWARD TIMMERMAN * * SUBROUTINES CALLED: LINTERP * * SYSTEM DEPENDENT CODE: USE OF VIRTUAL MEMORY AREA AND EXTENDED * MEMORY AREA (VMA/EMA) FEATURE OF HP1000 * RTE-A TO STORE 1,330,560 TWO BYTE WORDS * IN TIDEM VECTOR. OTHERWISE, ALL CODE * IS STANDARD FORTRAN 77. * * DATE: JANUARY 13, 1986. * * *============================================================================ * MODIFICATION HISTORY *============================================================================ * * REVISED 28MAY86 ADDING FUNCTION IZSAME TO CORRECT ERROR * IN HANDLING "SAME QUAD" DATA. ALSO CORRECTED LOGIC BUG * IN TIDAL PREDICTIONS NEAR THE EQUATOR. * * REVISED 22OCT86 ADDING TEST FOR NO DATA AVAILABLE FROM * THE SCHWIDERSKI MODELS. THE VMA BACKING STORE FILE WAS * MODIFIED TO CONTAIN 9999'S AT COORDINATES THAT HAVE NO DATA. * * REVISED 15JAN87 REPLACING THE JUNKINS INTERPOLATING ROUTINE * CALLED "INTERP" WITH THE GENERAL PURPOSE BI-LINEAR (AND LINEAR, * IF FEWER THAN 4 POINTS ARE OVER OCEAN) INTERPOLATING ROUTINE * CALLED "LINTERP". CLEANED UP SOME LOOSE ENDS AND REMOVED THE * CHECKING FOR "NODATA" FROM THE TIDES ROUTINE AND INCORPORATED * "NODATA" CHECKING IN THE INTERPOLATING ROUTINE. PREVIOUSLY * USED SUBROUTINE "IZSAME" WAS REMOVED SINCE THE NEW INTERPOLATOR * CAN BE IMMEDIATELY CALLED IF THE CURRENT X,Y IS IN THE SAME * 1 BY 1 DEGREE QUAD AS THE PREVIOUS CALL TO TIDES. "IZSAME" * USED TO DETERMINE IF IN THE SAME 30 BY 30 MINUTE QUAD AS THE * LAST CALL. * * OCT 15, 1987...FIXED TWO SET OF ERRORS: * * (1) TWO PARAMETERS WERE MISTYPED... * HX3 WAS 3.30D-4, IT SHOULD BE 3.03D-4; * PX2 WAS 4069.034032957D-0, IT SHOULD END IN ...575 * * (2) THE WEIGHTS FOR LONGITUDNAL VALUES IN THE LINTERP SUBROUTINE * WERE INCORRECT. A WEIGHT OF 1 WAS USED WHEN COS(LAT) SHOULD * HAVE BEEN IN CASES WHERE THERE WERE FEWER THAN 4 CORNERS OVER * WATER. * * * JAN 8, 1988...FIXED Q1 MODEL ERROR AS PER LETTER FROM P. L. * WOODWORTH. * THE LINE: * CHI(8) = H0 - 3.D0*S0 - HALFPI * * WAS CHANGED TO: * CHI(8) = H0 - 3.D0*S0 - HALFPI + P0 * * * MARCH 17, 1988...ADDED CODE FOR NODAL VARIATION CORRECTIONS AS * DESCRIBED BY DAVID E. CARTWRIGHT IN HIS LETTER * TO LAURY MILLER, 3/8/88. * * * * * AUG 15, 1988...SUBTLE ERROR AT TEST OF WHETHER WE CAN SKIP RECOMPUTING * THE DAY DEPENDENT TERMS. IF LASTDATA = .FALSE., * SAMEQUAD = .FALSE., AND .NOT. FIRST, THEN THE TERMS ARE * SKIPPED. BUT WHAT IF THIS IS THE FIRST POINT OVER * WATER AND A NEW QUAD? ==> THE ASTRONOMICAL TERMS WILL * NOT GET COMPUTED FOR THAT ENTIRE DAY!!!!!! * * CHANGED LINE OF CODE: * * IF (.NOT. FIRST ) THEN * * TO...... * * IF (.NOT. FIRST .AND. LASTDATA) THEN * *********************************************************************** ****************** FORMAL VARIABLES ************************** *********************************************************************** * * * VARIABLE I/O TYPE PURPOSE * * -------- ----- ---- ------- * * * * X INPUT REAL*8 NORTH LATITUDE, IN * * DECIMAL DEGREES. * * (-78 < X < 89) * * * * Y INPUT REAL*8 EAST LONGITUDE, IN * * DECIMAL DEGREES. * * * * TIME INPUT INTEGER*4 UNIVERSAL TIME OF OBS IN * * SECONDS SINCE BEGINNING * * OF YEAR. * * * * TIDEC OUTPUT REAL*8 OCEAN TIDE FOR POINT X,Y * * AT INPUT TIME. * * * * ISDATA OUTPUT LOGICAL ISDATA = .FALSE. IF NO * * MODEL INFORMATION IS * * AVAILABLE FOR POINT X,Y. * * OTHERWISE, ISDATA=.TRUE. * * * *********************************************************************** ***************** COMMON BLOCKS ***************************** *********************************************************************** * COMMON BLOCK /TEST/ CONTAINS.... * * * * VARIABLE TYPE PURPOSE * * -------- ---- ------- * * YEAR INTEGER*2 YEAR OF OBSERVATION. (I.E., 1985) * * * * FIRST LOGICAL .TRUE. IF FIRST CALL TO SUBROUTINE * * TIDES. FIRST IS RESET TO .FALSE. * * BY THIS SUBROUTINE AFTER THE * * INITIAL CALL. THIS VARIABLE * * ALLOWS BYPASSING COSTLY COMPUTER * * COMPUTATIONS ON SUBSEQUENT CALLS. * *CSR* .TRUE. TRIGGERS READING THE LARGE * *CSR* FILE TO FILL /BIG/ TIDEM(1330560) * *CSR* FACT(11) REAL FACTOR TO MULTIPLY EACH CONSTITUENT * *CSR* BY FOR SUM RETURNED IN TIDEC * *CSR* TICON(11) REAL TIDE HEIGHT FOR EACH CONSTITUENT * * * * COMMON BLOCK /BIG/ CONTAINS THE TIDEM ARRAY, 1330560 INTEGER*2 * * ELEMENTS IN VMA/EMA. * * * *********************************************************************** *********************************************************************** IMPLICIT REAL*8 (A-H,O-Z) REAL*8 FREQ(11),XSINP(11),XCOSP(11),CHI(11),X,Y,TIDEC REAL*8 PF(11),PU(11),QF(11),QU(11) INTEGER*2 ACOSP(4,11),ASINP(4,11),TIDEM(1330560),YEAR, . VCOSP(44),VSINP(44),LAT,LONG,INDX,COLAT INTEGER*4 ILAT1,ILAT2,TIME,TIMELAST,OFFSET LOGICAL NODAL LOGICAL FIRST,ISDATA,LASTDATA,NODATA,SAMEQUAD,DEBUG PARAMETER (DEBUG = .FALSE.) PARAMETER (PI = 3.141592653589793D0) PARAMETER (DEGRAD = PI/180.D0) PARAMETER (HALFPI = PI/2.D0) * * DEFINE CONSTANT COEFFICIENTS USED TO COMPUTE H0, S0, AND P0, * THE ASTRONOMICAL ARGUMENTS IN THE TIDE COMPUTATION. * * (H0,S0,P0)= MEAN LONGITUDES OF SUN, MOON, AND LUNAR PERIGEE * AT GREENWICH MIDNIGHT. * * (ALL ARE IN RADIANS.) * PARAMETER (HX1=279.69668D0*DEGRAD) PARAMETER (HX2=36000.768930485D0*DEGRAD) PARAMETER (HX3=3.03D-4*DEGRAD) PARAMETER (SX1=270.434358D0*DEGRAD) PARAMETER (SX2=481267.88314137D0*DEGRAD) PARAMETER (SX3=1133.D-6*DEGRAD) PARAMETER (SX4=1.9D-6*DEGRAD) PARAMETER (PX1=334.329653D0*DEGRAD) PARAMETER (PX2=4069.0340329575D0*DEGRAD) PARAMETER (PX3=10325.D-6*DEGRAD) PARAMETER (PX4=1.2D-5*DEGRAD) * * QX1, QX2, QX3 VALUES ADDED 3/17/88, PROVIDED BY D. E. CARTWRIGHT. * PARAMETER (QX1=259.1833D0*DEGRAD) PARAMETER (QX2=1934.14201D0*DEGRAD) PARAMETER (QX3=2078.0D-6*DEGRAD) * * EQUIVALENCE, FOR CONVIENENCE, THE ASINP & ACOSP ARRAYS * WITH VECTORS VSINP & VCOSP. THESE VECTORS ARE USED * IN CALLS TO THE INTERPOLATING SUBROUTINE, LINTERP. * EQUIVALENCE (ASINP,VSINP),(ACOSP,VCOSP) * * "SAVE" ALL LOCAL VARIABLES TO ALLOW SAME DAY & SAME QUAD * COMPUTING SHORTCUTS ON SUBSEQUENT CALLS! * SAVE * * COMMON BLOCK /BIG/ IN VMA/EMA. * * COMMON/BIG/TIDEM * * COMMON BLOCK /TEST/ CONTAINS THE YEAR OF THE OBSERVATION * AND LOGICAL VARIABLE FIRST. THIS COMMON BLOCK MUST BE * INITIALIZED IN THE CALLING PROGRAM WITH THE OBSERVATION * YEAR AND FIRST=.TRUE. PRIOR TO THE FIRST CALL TO TIDES. * COMMON/TEST/YEAR,FIRST,FACT(11),TICON(11) * * DEFINE FREQUENCIES OF PARTIAL TIDES FOR THE 11 MODELS. * (UNIT IS 1/SEC.) * DATA (FREQ(I),I=1,11)/1.40519D-4,1.45444D-4,0.72921D-4, . 0.67598D-4,1.37880D-4,0.72523D-4,1.45842D-4,0.64959D-4, . 0.053234D-4,0.026392D-4,0.003982D-4/ * * THESE VALUES PROVIDED 3/17/88 BY D.E CARTWRIGHT FOR NODAL VARIATION * CORRECTIONS. * DATA QU/-2.14D0, 0.D0, -8.86D0, 10.80D0, -2.14D0, 0.D0, -17.74D0, . 10.80D0, -23.74D0, 0.D0, 0.D0/ DATA QF/-37.D-3, 0.D0, 115.D-3, 187.D-3, -37.D-3, 0.D0, 286.D-3, . 187.D-3, 413.D-3, -130.D-3, 0.D0/ DATA NODAL/.TRUE./ * * INITIALIZE TIDEC VARIABLE (RETURNED OCEAN TIDE). * SET LOGICAL VARIABLE "LASTDATA" ON FIRST CALL TO TIDES. * TIDEC = 0.D0 *CSR IF (FIRST) CALL RDTIDE *CSR IF (FIRST) LASTDATA = .TRUE. * * TRANSLATE LAT, LONG TO COORDINATE SYSTEM COMPATABLE WITH * M2 TIDE TAPE. TIDE TAPE VALUES REFER TO "CENTERS" OF * (COLAT,LONG) PAIRS, THE CENTER BEING AT (COLAT-.5,LONG-.5). * THUS SHIFT X,Y HALF A DEGREE TO THE SOUTH AND TO THE EAST. * XTRANS = X - .5D0 YTRANS = Y + .5D0 * * CHECK THAT THE INPUT LAT AND LONG (X,Y) ARE LEGAL. I.E., -78 < LAT < 89, * AND 0 <= LONG <= 360. IF NOT, SET ISDATA FLAG TO .FALSE. AND RETURN * IMMEDIATELY TO CALLING ROUTINE. SET OTHER VARIABLES FOR SUBSEQUENT CALLS. * IF ( (X .LT. -78.0D0) .OR. (X .GT. 89.0D0) .OR. . (Y .LT. 0.0D0) .OR. (Y .GT. 360.0D0) ) THEN ISDATA = .FALSE. LASTDATA = .FALSE. FIRST = .FALSE. TIMELAST = TIME XLAST = XTRANS YLAST = YTRANS RETURN ENDIF * ACCOUNT FOR TRANSLATING PASS GREENWICH. IF(YTRANS .GT. 360.D0) YTRANS = YTRANS - 360.D0 * CHECK IF THE FIRST CALL TO TIDES. IF NOT..... * * * IF SAME 1 DEGREE QUAD AS LAST CALL TO TIDES, JUST INTERPOLATE * FOR NEW X,Y POINT, BYPASSING ACCESSING THE VMA ARRAY AND * THE COMPUTATION OF POSITION-DEPENDENT TERMS (I.E. GO TO 3 OR 4). * * IF SAME DAY AS LAST CALL TO TIDES, BYPASS RE-COMPUTING * THE DAY DEPENDENT TERMS -- THE ASTRONOMICAL ARGUMENTS (I.E. GO TO 4). * * IF LASTDATA IS FALSE THERE IS NO DATA FOR THIS QUAD, SET ISDATA TO * FALSE AND RETURN WITHOUT FURTHER PROCESSING. * IF (.NOT.FIRST) THEN SAMEQUAD = (IDINT(XLAST) .EQ. IDINT(XTRANS) .AND. . IDINT(YLAST) .EQ. IDINT(YTRANS)) * TEST THAT THE SIGNS ARE THE SAME, EITHER BOTH POSITIVE OR * BOTH NEGATIVE. (USE THE FORTRAN '.EQV.' LOGICAL OPERATOR.) SAMEQUAD = SAMEQUAD .AND. . ( (XLAST .GT. 0.0D0) .EQV. (XTRANS .GT. 0.0D0) ) IF (SAMEQUAD .AND. LASTDATA) THEN DO 1 J= 1,11 INDX = (J-1) * 4 + 1 CALL LINTERP(VCOSP(INDX),XTRANS,YTRANS,DATA,NODATA) XCOSP(J) = DATA CALL LINTERP(VSINP(INDX),XTRANS,YTRANS,DATA,NODATA) XSINP(J) = DATA 1 CONTINUE IF( (TIMELAST / 86400) .EQ. (TIME / 86400) ) GO TO 4 GO TO 3 ELSEIF (.NOT. LASTDATA .AND. SAMEQUAD) THEN ISDATA = .FALSE. RETURN ENDIF ENDIF ***************************************************************** **** P O S I T I O N D E P E N D E N T T E R M S **** ***************************************************************** * * ASINP AND ACOSP STORE THE VALUES AMPLITUDE * COS(PHASE) AND * AMPLITUDE * SIN(PHASE), RESPECTIVELY, AT THE 4 CORNERS OF A 1 DEGREE * SQUARE CONTAINING THE OBSERVATION'S GEODETIC POSITION. FOR EXAMPLE: * * LAT + 1 * ACOSP(4,J)---> X*********X <---ACOSP(1,J) * * * * LONG * * LONG + 1 * * * * ACOSP(3,J)---> X*********X <---ACOSP(2,J) * LAT * * ACOSP & ASINP VALUES ARE IN MILLIMETERS. * * *************************************************************** * THE TIDEM VECTOR CONTAINS 11 TIDE MODELS STORED AS FOLLOWS: * THE DATA ARE ARRANGED BY COLATITUDES (E.G. 90-LAT) * (N IN INTEGER DEGREES: N = 1,2,...,168) IN 11 CONSECUTIVE PAIRS * OF BLOCKS, EACH BLOCK CONTAINING 720 WORDS ARRANGED BY LONGITUDES * (M IN INTEGER DEGREES: M = 1,2,...,360) WITH THE FIRST * 360 WORDS CONTAINING AMPLITUDE * COS( PHASE ) AND THE * SECOND 360 WORDS CONTAINING AMPLITUDE * SIN( PHASE ). * * COLAT 1 COLAT2 * |-----------^-------------|-----------^-------------|------ * MODEL1 MODEL2 ... MODEL11 MODEL1 MODEL2 ... MODEL11 ..ETC. * 720WDS 720WDS 720WDS 720WDS 720WDS 720WDS ..ETC. * * THUS THERE ARE 11 * 720 WORDS FOR EACH COLATITUDE, * AND 168 * 11 * 720 = 1,330,560 WORDS IN THE TIDEM VECTOR. * * M2 TIDAL CONSTANTS ARE REPRESENTATIVE FOR THE ONE-BY-ONE * DEGREE SQUARE WITH GEOGRAPHICAL CENTER POINT (An,Bm) WHERE: * * An = n - 0.5 = COLATITUDE IN DEGREES (n= 1,2,...,168) * Bm = m - 0.5 = LONGITUDE IN DEGREES (m= 1,2,...,360). * * LAT = IDINT(XTRANS) LONG = IDINT(YTRANS) COLAT = 90 - LAT * * IF XTRANS < 0, THEN THE POINT IS IN THE SOUTHERN * HEMISPHERE. ADD ONE TO COLAT TO PRESERVE THE ORDERING * SCHEME OF THE ACOSP & ASINP VECTORS. * NOTE M2 TIDE TAPE NUMBERS LONGITUDES 1 TO 360 RATHER THAN * 0 TO 359. * IF(XTRANS .LT. 0.D0) COLAT = COLAT + 1 IF(LONG .EQ. 0) LONG = 360 ILAT1 = (COLAT - 1) * 7920 ILAT2 = (COLAT - 2) * 7920 DO 2 J=1,11 OFFSET = (J-1) * 720 ACOSP(1,J) = TIDEM(ILAT2 + MOD(LONG,360)+1 + OFFSET) TIDELM = ILAT2 + MOD(LONG,360)+1 + OFFSET ACOSP(2,J) = TIDEM(ILAT1 + MOD(LONG,360)+1 + OFFSET) ACOSP(3,J) = TIDEM(ILAT1 + LONG + OFFSET) ACOSP(4,J) = TIDEM(ILAT2 + LONG + OFFSET) * * CALL INTERPOLATING ROUTINE TO GET THE VALUE OF ACOSP AT X,Y. * STORE RESULT IN XCOSP(J). * INDX = (J-1) * 4 + 1 **DEBUG IF (DEBUG) THEN WRITE(7,*) 'VCOSP,INDX ',(VCOSP(KK),KK=INDX,INDX+3),INDX WRITE(7,*) 'TID ARRAY ELEMENT 1 = ',TIDELM WRITE(7,*) ' OFFSET = ',OFFSET ENDIF CALL LINTERP(VCOSP(INDX),XTRANS,YTRANS,DATA,NODATA) * IF NODATA IS .TRUE. , THEN ALL 4 CORNERS ARE OVER LAND. SET * VARAIBLES FOR NEXT CALL TO TIDES AND RETURN. IF ( NODATA ) THEN ISDATA = .FALSE. LASTDATA = .FALSE. FIRST = .FALSE. TIMELAST = TIME XLAST = XTRANS YLAST = YTRANS RETURN ENDIF XCOSP(J) = DATA OFFSET = OFFSET + 360 ASINP(1,J) = TIDEM(ILAT2 + MOD(LONG,360)+1 + OFFSET) ASINP(2,J) = TIDEM(ILAT1 + MOD(LONG,360)+1 + OFFSET) ASINP(3,J) = TIDEM(ILAT1 + LONG + OFFSET) ASINP(4,J) = TIDEM(ILAT2 + LONG + OFFSET) * * INTERPOLATE, STORING RESULT IN XSINP(J) * **DEBUG IF (DEBUG) THEN WRITE(7,*) 'VSINP,INDX ',(VSINP(KK),KK=INDX,INDX+3),INDX ENDIF CALL LINTERP(VSINP(INDX),XTRANS,YTRANS,DATA,NODATA) XSINP(J) = DATA 2 CONTINUE ************************************************************* **** D A Y D E P E N D E N T T E R M S *** ************************************************************* * * IF A NEW DAY, COMPUTE DAY OF THE YEAR THE OBSERVATION WAS TAKEN. * IF SAME DAY AS LAST CALL TO TIDES, BYPASS RE-COMPUTING THE DAY * DEPENDENT CONTRIBUTIONS TO THE TIDE MODEL (GO TO 4). NOTE HAVE * TO CHECK THAT THERE WAS DATA ON THE PREVIOUS CALL. * IF( .NOT. FIRST .AND. LASTDATA) THEN IF( (TIMELAST / 86400) .EQ. (TIME / 86400) ) GO TO 4 ENDIF 3 IDAY = TIME/86400 + 1 * * COMPUTE NUMBER OF DAYS SINCE BEGINNING OF 1975. * 'YEAR' IS IN COMMON BLOCK /TEST/. * JDAY = IDAY + 365*(YEAR-1975) + (YEAR-1973)/4 * * COMPUTE UNIVERSAL DAY (JULIAN CENTURIES.) * T = (27392.500528D0 + 1.0000000356D0*JDAY) / 36525.D0 * * COMPUTE MEAN LONGITUDES OF SUN, MOON, & LUNAR PERIGEE AT * GREENWICH MIDNIGHT. * H0 = HX1 + (HX2 + HX3*T)*T S0 = SX1 + (SX2 + (SX4*T - SX3)*T)*T P0 = PX1 + ((-PX4*T - PX3)*T + PX2)*T * * CHI IS THE ASTRONOMICAL ARGUMENT OF THE PARTIAL TIDE. * * CHI(1) = M2 (PRINCIPAL LUNAR, SEMIDIURNAL) * CHI(2) = S2 (PRINCIPAL SOLAR, " ) * CHI(3) = K1 (DECLINATION LUNI-SOLAR, DIURNAL) * CHI(4) = O1 (PRINCIPAL LUNAR, DIURNAL) * CHI(5) = N2 (ELLIPTICAL LUNAR, SEMIDIURNAL) * CHI(6) = P1 (PRINCIPAL SOLAR, DIURNAL) * CHI(7) = K2 (DECLINATION LUNI-SOLAR, SEMIDIURNAL) * CHI(8) = Q1 (ELLIPTICAL LUNAR, DIURNAL) * CHI(9) = MF (FORTNIGHTLY LUNAR, LONG PERIOD) * CHI(10) = MM (MONTHLY LUNAR, LONG PERIOD) * CHI(11) = SSA (SEMIANNUAL SOLAR, LONG PERIOD) * * CHI(1) = 2.D0*(H0-S0) CHI(2) = 0.D0 CHI(3) = H0 + HALFPI CHI(4) = H0 - 2.D0*S0 - HALFPI CHI(5) = 2.D0*H0 - 3.D0*S0 + P0 CHI(6) = -H0 - HALFPI CHI(7) = 2.D0*H0 CHI(8) = H0 - 3.D0*S0 - HALFPI + P0 CHI(9) = 2.D0*S0 CHI(10) = S0 - P0 CHI(11) = 2.D0*H0 * * COMPUTE UNIVERSAL TIME OF DAY IN SECONDS. * 4 UTS = MOD(TIME,86400) * * COMPUTE OCEAN TIDE. ACCUMULATE CONTRIBUTION FROM EACH MODEL * IN TIDEC. NOTE NODAL FACTORS PF, PU ARE COMPUTED ONLY FOR THE * FIRST DAY, THEN HELD CONSTANT FOR THE REMAINDER OF THE SEQUENCE. * THIS IS ACCURATE ENOUGH FOR A SEQUENCE SPANNING 6 MONTHS OR LESS * ---SO SEZ D. E. CARTWRIGHT 3/17/88. * * NODAL VARIATION FACTORS PF, PU COMPUTED AS PER D. E. CARTWRIGHT 3/17/88. * IF (NODAL) THEN Q0 = QX1 - (QX2 - QX3 * T) * T CNODE = DCOS(Q0) SNODE = DSIN(Q0) DO 1003 I = 1,11 PU(I) = QU(I) * SNODE * DEGRAD 1003 PF(I) = 1.D0 + QF(I) * CNODE NODAL = .FALSE. ENDIF * COMPUTE TIDE. DO 5 I = 1,11 TARG = FREQ(I) * UTS + CHI(I) + PU(I) TICON(I) = (XCOSP(I) * DCOS(TARG) . + XSINP(I) * DSIN(TARG)) * PF(I) * FACT(I) TIDEC = TIDEC + TICON(I) c TIDEC = TIDEC + (XCOSP(I) * DCOS(TARG) c . + XSINP(I) * DSIN(TARG)) * PF(I) * FACT(I) C WRITE (*,*) '----- IN TIDF1 -----' C WRITE(*,88)TIDEC,TARG 88 FORMAT('TIDEC, TARG: ',2E16.9) 5 CONTINUE * * SET VARIABLES FOR NEXT CALL TO TIDES. * XLAST = XTRANS YLAST = YTRANS TIMELAST = TIME ISDATA = .TRUE. LASTDATA = .TRUE. FIRST = .FALSE. RETURN END subroutine linterp(corner_data, lat, long, mean, nodata) * * PURPOSE * * This routine computes a weighted mean of the four integer values * in parameter array corner_data. The weights are computed as in a * bi-linear interpolation. The form of the two dimensional linear * interpolation was taken from subroutine CSINT by Laury Miller. * * -- Ed Timmerman 14Jan87. implicit real*8 (a-h,o-z) integer*2 corner_data(4) real*8 lat, long, mean, dlat, dlong, degrad real*8 cos_lat, X1Weight, X2Weight, Y1Weight, Y2Weight real*8 X1, X2, Y1, Y2, W(8), WeightSum logical X1Flag, X2Flag, Y1Flag, Y2Flag, isdata(4), nodata parameter (degrad = 3.1415926535898d0 / 180.0d0) * Compute DLAT and DLONG, the "Delta" latitude and longitude. * This is a "fast & dirty" approximation of the distances from * the point (lat,long) to the northeastern corner of the 1-degree * quadrangle containing it. These "distances" are used for weights. cos_lat = DCOS(lat*degrad) if (lat .lt. 0.d0) then dlat = (IDINT(lat) - lat) else dlat = (IDINT(lat+1) - lat) endif dlong = (IDINT(long+1) - long) * cos_lat * Set logical flags, depending if over land (9999) or not. isdata(1) = (corner_data(1) .ne. 9999) isdata(2) = (corner_data(2) .ne. 9999) isdata(3) = (corner_data(3) .ne. 9999) isdata(4) = (corner_data(4) .ne. 9999) X1Flag = ( isdata(1) .or. isdata(4) ) X2Flag = ( isdata(2) .or. isdata(3) ) Y1Flag = ( isdata(1) .or. isdata(2) ) Y2Flag = ( isdata(3) .or. isdata(4) ) * Test if all 4 corners are over land. If so, set the nodata flag to * .true. and return. Otherwise, nodata = .false. indicating that there * is valid data. nodata = .false. if ( .not. (Y1Flag .or. Y2Flag)) then nodata = .true. return endif * Set weights for the X1, X2, Y1, and Y2 components. Note that at least * one X and Y must be nonzero if there is at least one corner not over * land. * Set weights for X1 and X2 components. if (X1Flag) then if (X2Flag) then X1Weight = 1.0d0 - dlat X2Weight = dlat else X1Weight = 1.0d0 X2Weight = 0.0d0 endif else * ASSERTION: X2FLAG must be .true. X1Weight = 0.0d0 X2Weight = 1.0d0 endif * Set weights for Y1 and Y2 components. if (Y1Flag) then if (Y2Flag) then Y1Weight = cos_lat - dlong Y2Weight = dlong else Y1Weight = cos_lat Y2Weight = 0.0d0 endif else * ASSERTION: Y2FLAG must be .true. Y1Weight = 0.0d0 Y2Weight = cos_lat endif * Initialize weights for the corner components of X1, X2, Y1, Y2 * to zero. do 1 i = 1,8 1 W(i) = 0.0d0 * Set weights for the corner components. This could have been done * in the above sections where the X's and Y's were set, but the code * is clearer (?) if separated out. if ( isdata(1) .and. isdata(4) ) then W(1) = cos_lat - dlong W(2) = dlong elseif ( (.not. isdata(4)) .and. isdata(1) ) then W(1) = cos_lat elseif ( (.not. isdata(1)) .and. isdata(4) ) then W(2) = cos_lat endif if ( isdata(2) .and. isdata(3) ) then W(3) = cos_lat - dlong W(4) = dlong elseif ( (.not. isdata(3)) .and. isdata(2) ) then W(3) = cos_lat elseif ( (.not. isdata(2)) .and. isdata(3) ) then W(4) = cos_lat endif if ( isdata(1) .and. isdata(2) ) then W(5) = 1.0d0 - dlat W(6) = dlat elseif ( (.not. isdata(2)) .and. isdata(1) ) then W(5) = 1.0d0 elseif ( (.not. isdata(1)) .and. isdata(2) ) then W(6) = 1.0d0 endif if ( isdata(3) .and. isdata(4) ) then W(7) = dlat W(8) = 1.0d0 - dlat elseif ( (.not. isdata(4)) .and. isdata(3) ) then W(7) = 1.0d0 elseif ( (.not. isdata(3)) .and. isdata(4) ) then W(8) = 1.0d0 endif * Compute weighted mean. (Bi-linear interpolation.) if (X1Flag) then X1 = ( W(1) * corner_data(1) + . W(2) * corner_data(4)) / (W(1)+W(2)) * X1Weight else X1 = 0.0d0 endif if (X2Flag) then X2 = ( W(3) * corner_data(2) + . W(4) * corner_data(3)) / (W(3)+W(4)) * X2Weight else X2 = 0.0d0 endif Y1 = ( W(5) * corner_data(1) + . W(6) * corner_data(2)) * Y1Weight Y2 = ( W(7) * corner_data(3) + . W(8) * corner_data(4)) * Y2Weight WeightSum = X1Weight + X2Weight + Y1Weight + Y2Weight mean = ( (X1+X2) + (Y1+Y2) ) / WeightSum return end SUBROUTINE RDTIDE C C...TO READ THE SCHWIDERSKI DATAQ BASE FILE USED IN SUBROUTINE TIDES C INTEGER TIDEM COMMON /BIG/ TIDEM(1330560) C c OPEN (10,STATUS='OLD',ERR=998) OPEN (10,ERR=998) C READ (10,100,ERR=999) (TIDEM(J),J=1,1330560) 100 FORMAT (16I5) WRITE (6,*) ' SCHWIDERSKI TIDE FILE SUCCESSFULLY READ' RETURN C C...ERRORS 998 WRITE (6,*) ' *** ERROR OPENING tidf2.tide ' CALL EXIT (1) 999 WRITE (6,*) ' *** ERROR READING tidf2.tide ' CALL EXIT (2) C END