subroutine csrtide (dmjd,dlat,dlon,tide,tlp,tds,isdata, . u, v, pseudo, radiation) c...To compute diurnal and semidiurnal ocean tides from an "orthoweight" c file and long period tides from equilibrium theory for a given c time, latitude and longitude. c...Coded by Richard Eanes, University of Texas Center for Space Research c eanes@csr.utexas.edu c...Version 2.0 - 03 Nov 94 c - 13 Nov 94, added better documentation on djmd. c...Inputs: (double precision) c dmjd Modified Julian Date in Days, number of days (including fraction) c since 17 Nov 1858 00:00. The standard astronomical julian date c (jd) is related to dmjd by jd = 2400000.5 + dmjd. The modified c julian date of 1 Jan 1958 00:00 is 36204. The modified julian c date of 1 Jan 1994 12:00 is 49353.5. c dlat North Latitude (in degrees) c dlon East Longitude (in degrees) c c...Outputs (double precision) c tds Diurnal and semidiurnal tide value in millimeters c tlp Long period tide value in millimeters (equilibrium theory) c tide Total tide value in millimeters c isdata Logical, true if location is within the domain of the c tide model, false otherwise. When false the c tide values are returned as zero. c The following outputs can be used for constructing maps of the tidal c fields at the major diurnal and semidiurnal constituents c c u,v Diurnal and semidiurnal orthoweights for input location c (can be used in subroutine admit to compute harmonic information) c pseudo Logical, switch passed to subroutine admit determining what c variety of orthoweights are on the tide file c radiation Logical, true if orthoweights modified for radiation potential implicit double precision (a-h,o-z) logical isdata,pseudo,radiation dimension u(6),v(6) smjd = dmjd*86400.d0 call lpeqmt (smjd,dlat,tlp) call tptide (dlat,dlon,smjd,tds,isdata,u,v,pseudo,radiation) if (isdata) then tlp = tlp*10.d0 ! convert from cm to mm tds = tds*10.d0 ! convert from cm to mm tide = tlp + tds else tlp = 0.d0 tds = 0.d0 tide = 0.d0 endif return end SUBROUTINE TPTIDE( DLAT,DLON,TIME,TIDE,ISDATA,U,V,PSEUDO, . RADIATION ) * * This code has only minor changes from the GSTIDE program used to * compute the Cartwright&Ray GeoSat tide solution. * * * Name - TPTIDE Name derivation - TOPEX/Poseidon derived tide * * Function - to compute the ocean tidal height at a given time * and location. * * Language - Fortran 77 * * Arguments - * name type I/O description * ---- ---- --- ----------- * DLAT D I north latitude (in degrees) for desired * location. * * DLON D I east longitude (in degrees). * * TIME D I desired time, in seconds of MJD. E.g., * Jan 1, 1986 00:00 is 4011638400. * * TIDE D O computed tidal height, in cm. * * ISDATA L O logical denoting whether tide data exist at * desired location. If FALSE, then TIDE is * not modified. * * U,V D O The orthoweights interpolated to DLAT,DLON * Dimensioned (3,2) * * PSEUDO L O True for C&R style "orthoweights" * This can be passed to admit.f for proper * interpretation of the orthoweights. * Set per integer on orthoweight file * * RADIATION L O True if the correction for radiation * potential is activated. * Set per integer on orthoweight file * * * Usage notes - * Using the input data as of Nov 19, 1990, this routine computes * the height of the (ocean + load) tide at the desired location. * This quantity is appropriate for use in satellite altimetry. * The computed tide is composed of the 30 largest spectral lines * within both the diurnal and semidiurnal bands, sufficient for * representing all major constituents including nodal modulations. * * Processing logic - * The tidal height is computed at four grid points surrounding * the desired location, with the final result being a bilinear * interpolation. * * File references - * The orthoweight data are read on initial call to routine: * File 51 - global array of orthoweights defining tidal * constants at all valid locations. * * Important local variables - * Array UV is loaded with orthoweights as follows: * UV(I,J,L,M,N), where * I = 1 for U; 2 for V * J = 1,2 or 3 (3 complex coeffs per species.) * L = tidal species (1 or 2) * M = 1 to NX, longitude index number * N = 1 to NY, latitude index * * For machines like the Cray, the DOUBLE PRECISION should be removed. * * Error processing - none. * * Technical references - * D. Cartwright and R. Ray, Oceanic tides from Geosat altimetry, * Journal of Geophysical Research, 95, 3069-3090, 1990. * (particularly Appendix A). * * History - * version date programmer change description * ------- ---- ---------- ------------------ * 1.0 8/17/90 Ray/Cartwright Initial version. * 1.1 12/18/90 Ray/Cartwright Enhanced documentation; * Use ascii input for portability. * 1.2 12/31/91 R. Ray Changed min W to 0.5 (old value of * 0.2 allowed too much extrapolation) * 2.0 2/25/94 R. Eanes Changed to use T/P derived orthoweights * 2.1 3/02/94 R. Eanes Generalized grid size up to 1 X 1 degree * and added data statements for the harmonic * amplitudes and phases (a la the JPL mods) * Interpolate orthoweights instead of tides * 2.2 3/03/94 R. Eanes Return orthoweights (u,v) for other use * Added ITYPE to switch to C&R style orthoweights * Change min W (WMIN) back to 0.2. * 2.3 4/18/94 R. Eanes Take out T/P kludge for latitudes > 66.2 * 2.4 5/20/94 R. Eanes Change treatment of extrapolating past the southern * and northernmost grid points. Set WMIN to 0.1. * 2.5 Added R.Ray's radiation tide correction * WMIN changed to 0.5 (again) * IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*80 TITLE1,TITLE2 REAL UV(2,3,2,360,180), UNDEF DOUBLE PRECISION LATMIN,LATMAX,LONMIN,LONMAX DIMENSION U(3,2),V(3,2) SAVE LOGICAL INIT,ISDATA,PSEUDO,RADIATION DATA INIT/.TRUE./, UNDEF/99999./, EPS /1.D-10/, WMIN /0.5/ * on first call, read T/P -derived ortho-weights * ------------------------------------------------ IF (INIT) THEN INIT = .FALSE. DO 4 J=1,180 DO 4 I=1,360 UV(1,1,1,I,J) = UNDEF 4 CONTINUE READ (51,14) TITLE1,TITLE2 WRITE (6,14) TITLE1,TITLE2 READ(51,11) NX,NY,LATMIN,LATMAX,LONMIN,LONMAX,IDIFF,ITYPE,IRAD IF (ITYPE.EQ.1) THEN PSEUDO = .TRUE. ! C&R style orthoweights ELSE PSEUDO = .FALSE. ! Corrected orthoweights ENDIF IF (IRAD.EQ.1) THEN RADIATION = .TRUE. ! Corrected for radiation potential ELSE RADIATION = .FALSE. ! Not corrected for radiation potential ENDIF IF (NX.GT.360 .OR. NY.GT.180) STOP 1 ! must change dimensions DX = (LONMAX-LONMIN)/(NX-1) DY = (LATMAX-LATMIN)/(NY-1) XLATS = LATMIN - DY ! southern limit of extrapolation XLATN = LATMAX + DY ! northern limit of extrapolation 8 CONTINUE READ(51,12,END=10) I,J,SLAT,SLON,U1,V1,U2,V2,U3,V3 READ(51,13) U4,V4,U5,V5,U6,V6 UV(1,1,1,I,J) = U1 UV(2,1,1,I,J) = V1 UV(1,2,1,I,J) = U2 UV(2,2,1,I,J) = V2 UV(1,3,1,I,J) = U3 UV(2,3,1,I,J) = V3 UV(1,1,2,I,J) = U4 UV(2,1,2,I,J) = V4 UV(1,2,2,I,J) = U5 UV(2,2,2,I,J) = V5 UV(1,3,2,I,J) = U6 UV(2,3,2,I,J) = V6 GO TO 8 10 CONTINUE 11 FORMAT(2I4,4F8.3,3I2) 12 FORMAT(2I4,2F8.3,6F8.3) 13 FORMAT(24X,6F8.3) 14 FORMAT(A80,/,A80) ENDIF ISDATA = .TRUE. TIDE = 0.0D0 DO I=1,3 DO J=1,2 U(I,J) = 0.D0 V(I,J) = 0.D0 ENDDO ENDDO * compute indices for desired position * ------------------------------------ XLAT = DLAT JLAT1 = INT((XLAT - LATMIN)/DY) + 1 IF (JLAT1.EQ.0 .AND. XLAT.GT.XLATS) JLAT1 = 1 ! allow extrapolation to XLATS IF (JLAT1.EQ.NY .AND. XLAT.LT.XLATN) JLAT1 = NY - 1 ! allow extrapolation to XLATN JLAT2 = JLAT1 + 1 IF (JLAT1.LT.1 .OR. JLAT2.GT.NY) THEN ISDATA = .FALSE. * if the orthoweight file repesents differences return zero * instead of saying there is no data * --------------------------------------------------------- IF (IDIFF.EQ.1) ISDATA = .TRUE. RETURN ENDIF XLON = DLON IF (XLON.LT.LONMIN) XLON = XLON + 360.D0 - EPS ! eps to avoid problem when xlon=lonmin ILON1 = INT((XLON - LONMIN)/DX) + 1 ILON2 = ILON1 + 1 IF (ILON2.GT.NX) ILON2 = 1 ! longitude grid must span globe IF (ILON1.LT.1 .OR. ILON1.GT.NX) STOP 301 ! should never happen. IF (UV(1,1,1,ILON1,JLAT1).EQ.UNDEF .AND. * UV(1,1,1,ILON2,JLAT1).EQ.UNDEF .AND. * UV(1,1,1,ILON1,JLAT2).EQ.UNDEF .AND. * UV(1,1,1,ILON2,JLAT2).EQ.UNDEF) THEN ISDATA = .FALSE. IF (IDIFF.EQ.1) ISDATA = .TRUE. RETURN ENDIF * Distances from corners in units of cell size * -------------------------------------------- WX1 = ( DX - (XLON - REAL(ILON1-1)*DX - LONMIN) )/DX WX2 = 1.D0 - WX1 WY1 = ( DY - (XLAT - REAL(JLAT1-1)*DY - LATMIN) )/DY WY2 = 1.D0 - WY1 IF (WY2.LT.0.) THEN ! north corner weights zero in southern extrapolation WY1 = WY1 + WY2 WY2 = 0. ENDIF IF (WY1.LT.0.) THEN ! south corner weights zero in northern extrapolation WY2 = WY1 + WY2 WY1 = 0. ENDIF * Interpolation weights: * W1,W2,W3,W4 are for northwest,northeast,southeast,southwest corners. * -------------------------------------------------------------------- W1 = WX1*WY2 W2 = WX2*WY2 W3 = WX2*WY1 W4 = WX1*WY1 * initialize sum of weights and orthoweights * ------------------------------------------ W = 0.D0 DO 50 I=1,3 DO 50 J=1,2 U(I,J) = 0.D0 50 V(I,J) = 0.D0 * interpolate the orthoweights * ---------------------------- IF (UV(1,1,1,ILON1,JLAT1).NE.UNDEF) THEN DO 100 I=1,3 U(I,1) = UV(1,I,1,ILON1,JLAT1)*W4 V(I,1) = UV(2,I,1,ILON1,JLAT1)*W4 U(I,2) = UV(1,I,2,ILON1,JLAT1)*W4 V(I,2) = UV(2,I,2,ILON1,JLAT1)*W4 100 CONTINUE W = W4 ENDIF IF (UV(1,1,1,ILON1,JLAT2).NE.UNDEF) THEN DO 200 I=1,3 U(I,1) = UV(1,I,1,ILON1,JLAT2)*W1 + U(I,1) V(I,1) = UV(2,I,1,ILON1,JLAT2)*W1 + V(I,1) U(I,2) = UV(1,I,2,ILON1,JLAT2)*W1 + U(I,2) V(I,2) = UV(2,I,2,ILON1,JLAT2)*W1 + V(I,2) 200 CONTINUE W = W + W1 ENDIF IF (UV(1,1,1,ILON2,JLAT2).NE.UNDEF) THEN DO 300 I=1,3 U(I,1) = UV(1,I,1,ILON2,JLAT2)*W2 + U(I,1) V(I,1) = UV(2,I,1,ILON2,JLAT2)*W2 + V(I,1) U(I,2) = UV(1,I,2,ILON2,JLAT2)*W2 + U(I,2) V(I,2) = UV(2,I,2,ILON2,JLAT2)*W2 + V(I,2) 300 CONTINUE W = W + W2 ENDIF IF (UV(1,1,1,ILON2,JLAT1).NE.UNDEF) THEN DO 400 I=1,3 U(I,1) = UV(1,I,1,ILON2,JLAT1)*W3 + U(I,1) V(I,1) = UV(2,I,1,ILON2,JLAT1)*W3 + V(I,1) U(I,2) = UV(1,I,2,ILON2,JLAT1)*W3 + U(I,2) V(I,2) = UV(2,I,2,ILON2,JLAT1)*W3 + V(I,2) 400 CONTINUE W = W + W3 ENDIF IF (W.GT.WMIN) THEN DO 500 I=1,3 DO 500 J=1,2 U(I,J) = U(I,J)/W 500 V(I,J) = V(I,J)/W CALL TPTID1( U,V,TIME,PSEUDO,RADIATION,TIDE ) ELSE ISDATA = .FALSE. IF (IDIFF.EQ.1) THEN ISDATA = .TRUE. TIDE = 0.0D0 DO I=1,3 DO J=1,2 U(I,J) = 0.D0 V(I,J) = 0.D0 ENDDO ENDDO ENDIF ENDIF RETURN END SUBROUTINE TPTID1( U,V,TS,PSEUDO,RADIATION,TIDE ) * Computes the tide at time TS from the orthoweights U,V. * This is a kernel routine, meant to be called only by TPTIDE. IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION INDX(30,2,4),AMP(30,2),FREQ(30,2), * PHA(30,2),CT(30,2),ST(30,2),AT(3),BT(3), * U00(2),U20(2),U21(2),U40(2),U41(2),V41(2), * P(3,2),Q(3,2),U(3,2),V(3,2), * SHPN(4),PHC(4),DPD(4) LOGICAL INIT,RADIATION,PSEUDO PARAMETER (RAD_FACTOR=0.97D0, RAD_PHASE=5.9D0) SAVE * DATA PHC/290.21D0, 280.12D0, 274.35D0, 343.51D0/, * DPD/13.1763965D0,0.9856473D0,0.1114041D0,0.0529539D0/, * U00/0.0298D0,0.0200D0/, U20/0.1408D0,0.0905D0/, * U21/0.0805D0,0.0638D0/, U40/0.6002D0,0.3476D0/, * U41/0.3025D0,0.1645D0/, V41/0.1517D0,0.0923D0/ DATA TC/40431744.D2/, TSLAST/-9.99E10/ DATA INIT/.TRUE./ DATA INDX !( L ,M,N) + /-3,-3,-2,-2,-2,-2,-1,-1,-1,-1,-1, 0, 0, 0, 0, !( 1:15,1,1) + 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, !(16:30,1,1) + -3,-3,-2,-2,-2,-1,-1,-1,-1,-1,-1, 0, 0, 0, 0, !( 1:15,2,1) + 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, !(16:30,2,1) + 0, 2, 0, 0, 2, 2, 0, 0, 0, 0, 2, 0, 0, 0, 2, !( 1:15,1,2) + -3,-2,-2, 0, 0, 0, 1, 2,-2, 0, 0,-2, 0, 0, 0, !(16:30,1,2) + 0, 2, 0, 2, 3,-1, 0, 0, 1, 2, 3,-2,-1, 0, 0, !( 1:15,2,2) + 1,-2, 0, 0, 0,-3,-2,-1, 0, 0, 0, 0,-2, 0, 0, !(16:30,2,2) + 2, 0, 1, 1,-1,-1, 0, 0, 0, 2, 0,-1, 1, 1,-1, !( 1:15,1,3) + 0, 0, 0, 0, 0, 0, 0, 0, 1,-1,-1, 0,-2, 0, 0, !(16:30,1,3) + 3, 1, 2, 0, 0, 1, 1, 1, 1,-1,-1, 2, 0, 0, 0, !( 1:15,2,3) + 0, 1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1,-1,-1, !(16:30,2,3) + 0, 0,-1, 0,-1, 0,-2,-1, 0, 0, 0, 0, 0, 1, 0, !( 1:15,1,4) + 0,-1, 0,-1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, !(16:30,1,4) + 0, 0, 0, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0,-1, 0, !( 1:15,2,4) + 0, 0, 0, 0, 1, 0, 0, 0,-1, 0, 1, 2, 0, 0, 1/ !(16:30,2,4) DATA AMP !( L ,M) + /0.00663,0.00802,0.00947,0.05019,0.0018 ,0.00954, !( 1: 6,1) + 0.00152,0.04946,0.26218,0.00171,0.00343,0.00741, !( 7:12,1) + 0.02062,0.00414,0.00394,0.00713,0.00137,0.122 , !(13:18,1) + 0.0073 ,0.36874,0.05002,0.00293,0.00524,0.00395, !(19:24,1) + 0.02061,0.00409,0.00342,0.00169,0.01128,0.00723, !(25:30,1) + 0.0018 ,0.00467,0.01601,0.01932,0.0013 ,0.00102, !( 1: 6,2) + 0.00451,0.121 ,0.00113,0.02298,0.00106,0.0019 , !( 7:12,2) + 0.00218,0.02358,0.63194,0.00193,0.00466,0.01786, !(13:18,2) + 0.00447,0.00197,0.01718,0.29402,0.00305,0.00102, !(19:24,2) + 0.07994,0.02382,0.00259,0.00086,0.00447,0.00195/ !(25:30,2) DATA PHA !( L ,M) + / 90., 90., 90., 90., 90., 90.,270., 90., 90.,270., !( 1:10,1) + 270.,270.,270.,270.,270., 13.,270., 90., 90.,270., !(11:20,1) + 270.,-13.,270.,270.,270.,270.,270.,270.,270.,270., !(21:30,1) + 0., 0., 0., 0., 77.,103.,180., 0., 77., 0., !( 1:10,2) + 77.,180.,103.,180., 0., 77.,180.,180., 0., 0., !(11:20,2) + 283., 0.,262.,180., 0., 0., 0., 0., 0., 0./ !(21:30,2) * IF (INIT) THEN INIT = .FALSE. RAD = DATAN(1.D0)/45.D0 DPLD = 36.D1 - DPD(1) + DPD(2) TWOPI = DATAN(1.D0)*8.D0 * FDPD = 0.D0 c WRITE(6,100) c 100 FORMAT(/' Semidiurnal & diurnal tidal potential terms:'/) DO 201 M=1,2 FDPD = FDPD + DPLD DO 101 L=1,30 c The data formerly read here is now in the data statements above c READ(9,*) (INDX(L,M,N),N=1,4), AMP(L,M), PHA(L,M) THETA = FDPD DO 1 N=1,4 THETA = THETA + REAL(INDX(L,M,N))*DPD(N) 1 CONTINUE FREQ(L,M) = THETA OMEGAT = 2.D0*THETA*RAD CT(L,M) = 2.D0*DCOS(OMEGAT) ST(L,M) = 2.D0*DSIN(OMEGAT) IF (RADIATION .AND. M.EQ.2 .AND. * (L.EQ.21.OR.L.EQ.22.OR.L.EQ.23)) THEN AMP(L,M) = AMP(L,M)*RAD_FACTOR PHA(L,M) = PHA(L,M) - RAD_PHASE WRITE(6,503) L,M,(INDX(L,M,N),N=1,4),FREQ(L,M),AMP(L,M), * PHA(L,M) ELSE c WRITE(6,502) L,M,(INDX(L,M,N),N=1,4),FREQ(L,M),AMP(L,M), c * PHA(L,M) ENDIF 101 CONTINUE c WRITE(6,*) ' ' 201 CONTINUE 502 FORMAT(1X,I3,I4,2X,4I2,F12.6,2(F10.5,F10.1,5X)) 503 FORMAT(1X,I3,I4,2X,4I2,F12.6,F10.5,F10.1,3X,'incl.rad.') ENDIF IF (TS.NE.TSLAST) THEN TD = (TS - TC)/864.D2 * Compute 4 principal mean longitudes for a given time TD DO 107 N=1,4 PH = PHC(N) + TD*DPD(N) SHPN(N) = DMOD(PH,36.D1) 107 CONTINUE FD = TD - DINT(TD) E = 36.D1*FD - SHPN(1) + SHPN(2) ENDIF * OTIDE = 0.D0 DO 300 M=1,2 IF (TS.NE.TSLAST) THEN * * Compute orthotides P(3),Q(3) for given time TS in terms * of potential amplitudes AT(3), BT(3) DO 207 K=1,3 AT(K) = 0.D0 BT(K) = 0.D0 207 CONTINUE DO 407 L=1,30 PH = M*E + PHA(L,M) DO 307 N=1,4 I = INDX(L,M,N) IF (I.EQ.0) GO TO 307 PH = PH + REAL(I)*SHPN(N) 307 CONTINUE THETA = PH*RAD A = AMP(L,M)*100.D0 CR = A*DCOS(THETA) SR = A*DSIN(THETA) CRCT = CR*CT(L,M) CRST = CR*ST(L,M) SRCT = SR*CT(L,M) SRST = SR*ST(L,M) AT(1) = AT(1) + CR BT(1) = BT(1) - SR AT(2) = AT(2) + CRCT BT(2) = BT(2) - SRCT AT(3) = AT(3) + CRST BT(3) = BT(3) - SRST 407 CONTINUE P(1,M) = U00(M)*AT(1) Q(1,M) = U00(M)*BT(1) P(2,M) = U20(M)*AT(1) - U21(M)*AT(2) Q(2,M) = U20(M)*BT(1) - U21(M)*BT(2) IF (PSEUDO) THEN * Cartwright and Ray pseudo-orthoweights P(3,M) = U40(M)*AT(1) - U41(M)*AT(2) - V41(M)*AT(3) Q(3,M) = U40(M)*BT(1) - U41(M)*BT(2) - V41(M)*BT(3) ELSE * Correct orthoweight definition P(3,M) = U40(M)*AT(1) - U41(M)*AT(2) + V41(M)*AT(3) Q(3,M) = U40(M)*BT(1) - U41(M)*BT(2) + V41(M)*BT(3) ENDIF ENDIF DO 30 K=1,3 OTIDE = OTIDE + U(K,M)*P(K,M) + V(K,M)*Q(K,M) 30 CONTINUE 300 CONTINUE TIDE = OTIDE TSLAST = TS RETURN END subroutine admit2 (pseudo,radiation,u,v,x,y,h1,h2,amp,pha) c c...to compute admittance and in-phase and quadrature components c from orthoweights c c...coded by Richard Eanes, UT/CSR, 1993 and 1994 c c Inputs: c pseudo -- Logical variable, true for CR91 style pseudo-orthoweights c Use output from tptide for this c radiation - Logical variable, true if the orthoweights are corrected c for the radiation potential at T2,S2,R2 (R2 not in this list) c u,v - The ortwhoweights for the diurnal and semidiurnal bands c c Outputs: c x,y - The tidal admittances at a list of frequencies spanning c the diurnal and semidiurnal bands (see data below) c h1,h2 - In-phase and quadrature amplitudes at the same frequencies (cm) c amp,pha - Amplitude and phase at the same list of frequencies. (cm,deg) c implicit double precision (a-h,o-z) parameter (nf=30) logical first,pseudo,radiation character*11 tname dimension u(3,2), v(3,2) dimension amp(100),pha(100),x(100),y(100),h1(100),h2(100) dimension tname(nf),freq(nf),per(nf),hsm(nf),chi(nf) dimension i1(0:2),i2(0:2),xp1(nf),xp2(nf),fxy(nf),sp(6,2) save c...the orthotide weight factors (signs with + were wrong in C&R 90 Table A1) data ((sp(i,m),i=1,6),m=1,2) / . 0.0298, 0.1408, +0.0805, 0.6002, +0.3025, 0.1517, . 0.0200, 0.0905, +0.0638, 0.3476, +0.1645, 0.0923/ c data (tname(i),freq(i),per(i),hsm(i),chi(i),i=1,nf) / c name f(cycle/d) period hs (m) phase c long period . 'LP 55.565', 0.000147094, 6798.382, 0.0279, 180.0, ! 1 . 'Sa 56.554', 0.002737779, 365.260, -0.0049, 0.0, ! 2 . 'Ssa 57.555', 0.005475819, 182.621, -0.0310, 0.0, ! 3 . 'TERa 58.554', 0.008213597, 121.749, -0.0018, 0.0, ! 4 . 'Mm 65.455', 0.036291647, 27.555, -0.0352, 0.0, ! 5 . 'Mf 75.555', 0.073202203, 13.661, -0.0666, 0.0, ! 6 . 'TERm 85.455', 0.109493850, 9.133, -0.0128, 0.0, ! 7 . '93a 93.555', 0.140928587, 7.096, -0.0020, 0.0, ! 8 c diurnal . '2Q1 125.755', 0.856952384, -6.904, -0.0066, 270.0, ! 9 . 'Q1 135.655', 0.893244048, -9.213, -0.0502, 270.0, !10 . 'O1 145.555', 0.929535695, -13.841, -0.2622, 270.0, !11 . 'M1 155.655', 0.966446233, -28.297, 0.0206, 90.0, !12 . 'P1 163.555', 0.997262079, -221.060, -0.1220, 270.0, !13 . 'S1 164.556', 1.000000000, -559.897, 0.0029, 90.0, !14 . 'K1 165.555', 1.002737898, 1050.246, 0.3688, 90.0, !15 . 'PHI1167.555', 1.008213699, 155.580, 0.0053, 90.0, !16 . 'J1 175.455', 1.039029527, 26.850, 0.0206, 90.0, !17 . 'OO1 185.555', 1.075940083, 13.485, 0.0113, 90.0, !18 . 'NU1 195.455', 1.112231730, 9.054, 0.0022, 90.0, !19 c semi-diurnal . '227 227.655', 1.828255527, -5.704, 0.0047, 0.0, !20 . '2N2 235.755', 1.859690264, -6.950, 0.0160, 0.0, !21 . 'MU2 237.555', 1.864547173, -7.193, 0.0193, 0.0, !22 . 'N2 245.655', 1.895981946, -9.295, 0.1210, 0.0, !23 . 'NU2 247.455', 1.900838820, -9.734, 0.0230, 0.0, !24 . 'M2 255.555', 1.932273593, -14.026, 0.6319, 0.0, !25 . 'L2 265.455', 1.968565204, -28.566, -0.0179, 180.0, !26 . 'T2 272.556', 1.997262163, -158.475, 0.0172, 0.0, !27 . 'S2 273.555', 2.000000000, -279.994, 0.2940, 0.0, !28 . 'K2 275.555', 2.005475795, 525.123, 0.0800, 0.0, !29 . '285 285.455', 2.041767407, 26.181, 0.0045, 0.0/ !30 data i1 / 1, 9, 20/, i2 / 8, 19, 30/ data dt /2.d0/ data rad_fac /0.97d0/, rad_phase /5.9d0/ ! Ray-Sanchez-Cartrwight radiation factors data first /.true./ c...compute constants depending only on frequency first = .true. if (first) then first = .false. twopi = 8.d0*atan(1.d0) deg = 360.d0/twopi do m = 1,2 mf = (-1)**m do i = i1(m),i2(m) theta = twopi*freq(i)*dt cw = 2.*cos(theta) sw = 2.*sin(theta) xp1(i) = sp(2,m) - sp(3,m)*cw if (pseudo) then * Cartwright&Ray pseudo-orthoweights xp2(i) = sp(4,m) - sp(5,m)*cw - sp(6,m)*sw else * True orthoweights xp2(i) = sp(4,m) - sp(5,m)*cw + sp(6,m)*sw endif fxy(i) = mf*abs(hsm(i))*1000.d0 ! scale to return values in mm c write (12,200) tname(i),fxy(i)*sp(1,m),fxy(i)*xp1(i), c . fxy(i)*xp2(i) enddo enddo endif do m = 1,2 do i = i1(m),i2(m) x(i) = sp(1,m)*u(1,m) + xp1(i)*u(2,m) + xp2(i)*u(3,m) y(i) = sp(1,m)*v(1,m) + xp1(i)*v(2,m) + xp2(i)*v(3,m) h1(i) = fxy(i)*x(i) h2(i) = -fxy(i)*y(i) amp(i) = sqrt(h1(i)**2 + h2(i)**2) if (amp(i).ne.0.) then pha(i) = atan2(h2(i),h1(i))*deg if (pha(i).lt.0.) pha(i) = pha(i) + 360. else pha(i) = 0. endif c...correct T2 and S2 for radiation potential if (radiation .and. (i.eq.27 .or. i.eq.28) ) then amp(i) = amp(i)*rad_fac pha(i) = pha(i) + rad_phase h1(i) = amp(i)*cos(pha(i)/deg) h2(i) = amp(i)*sin(pha(i)/deg) endif c write (6,100) i,m,tname(i),hsm(i),freq(i),x(i),y(i), c . h1(i),h2(i),amp(i),pha(i) enddo enddo 100 format (1x,2i2,1x,a11,1x,f8.4,f8.5,1x,4f9.4,1x,2f8.2) 200 format (a11,3x,3f15.6) return end SUBROUTINE LPEQMT( TS, DLAT, TLP ) * * Function - computes the long-period equilibrium ocean tides. * * Arguments - * name type I/O description * ---- ---- --- ----------- * TS D I modified Julian day, in seconds, denoting * time at which tide is to be computed. * * DLAT D I latitude in degrees (positive north) for the * position at which tide is computed. * * TLP D O computed long-period tide, in centimeters. * * * Processing logic - * Fifteen tidal spectral lines from the Cartwright-Tayler-Edden * tables are summed over to compute the long-period tide. * * Technical references - * Cartwright & Tayler, Geophys. J. R.A.S., 23, 45, 1971. * Cartwright & Edden, Geophys. J. R.A.S., 33, 253, 1973. * * History - * version date programmer change description * ------- ---- ---------- ------------------ * 1.0 11/27/90 D. Cartwright Documented by R. Ray * * IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION PHC(4),DPD(4),SHPN(4) SAVE * PARAMETER (RAD= 0.017453292519943 D0) PARAMETER (PSOL=283.D0*RAD) ! solar perigee held constant DATA PHC/290.21D0,280.12D0,274.35D0,343.51D0/, * DPD/13.1763965D0,0.9856473D0,0.1114041D0,0.0529539D0/, * TC/40431744.D2/ * * Compute 4 principal mean longitudes in radians at time TD * --------------------------------------------------------- TD = (TS - TC)/864.D2 ! days past 1 Jan 87 00:00 DO 1 N=1,4 PH = PHC(N) + TD*DPD(N) SHPN(N) = RAD*DMOD(PH,36.D1) 1 CONTINUE * * Assemble long-period tide potential from 15 CTE terms > 1 mm. * Nodal term is included but not the constant term. * ------------------------------------------------- ZLP = 2.79D0 * DCOS( SHPN(4) ) * -0.49D0 * DCOS( SHPN(2) - PSOL ) * -3.10D0 * DCOS( 2.D0*SHPN(2) ) PH = SHPN(1) ZLP = ZLP - 0.67D0 * DCOS( PH - 2.D0*SHPN(2) + SHPN(3) ) * -(3.52D0 - 0.46D0*DCOS(SHPN(4))) * DCOS( PH - SHPN(3) ) PH = PH + SHPN(1) ZLP = ZLP - 6.66D0 * DCOS( PH ) * - 2.76D0 * DCOS( PH + SHPN(4) ) * - 0.26D0 * DCOS( PH + 2.D0*SHPN(4) ) * - 0.58D0 * DCOS( PH - 2.D0*SHPN(2) ) * - 0.29D0 * DCOS( PH - 2.D0*SHPN(3) ) PH = PH + SHPN(1) ZLP = ZLP - 1.27D0 * DCOS( PH - SHPN(3) ) * - 0.53D0 * DCOS( PH - SHPN(3) + SHPN(4) ) * - 0.24D0 * DCOS( PH - 2.D0*SHPN(2) + SHPN(3) ) * * Multiply by gamma_2 * sqrt(5/4 pi) * P20(lat) * --------------------------------------------- S = DSIN(DLAT*RAD) TLP = 0.437D0*ZLP*(1.5D0*S*S - 0.5D0) RETURN END