program csrtide_test c...to compute example output from csrtide c...coded by: Richard Eanes, UT/CSR, eanes@csr.utexas.edu c...version 2.0 - 03 Nov 94 c...version 2.1 - 13 Nov 94 - Added definition of xmjd in first call to c csrtide. Fixed bug in 2.0 by changing dmjd c in second csrtide call to xmjd. c...version 3.0 - 15 Mar 95 - Changed orthoweight file name implicit double precision (a-h,o-z) double precision jd logical isdata,pseudo,radiation dimension mdyhms(6), u(3,2), v(3,2) dimension x(100),y(100),h1(100),h2(100),amp(100),pha(100) c...xmjd0=48988 is the modified julian date for 1 Jan 1993 00:00:00 in days c modified julian date is julian date - 2400000.5 c See subroutine JULDAY below to convert from calendar date to c julian date, and KALDAY to convert from julian date to calendar date data xmjd0 /48988.d0/, dmjd /0.125d0/, nmjd /9/ c...open the orthoweight file open (51, file='ortho_csr_3.0') c...read in and echo the desired location 20 write (6,*) write (6,*) 'Enter the location, lat,lon in deg (99,0) to stop' read (5,*) dlat,dlon write (6,*) ' Latitude = ' , dlat, ' (deg) ', . ' Longitude = ', dlon, ' (deg)' write (6,*) if (dlat.gt.91.) stop c...illustrate use of subroutine csrtide and subroutine admit2 c...admit2 is used to compute harmonic constants from the u,v orthoweights c see the documentation in the source code for the ordering of the c major tide harmonics in the output arrays x,y,h1,h2,amp,pha xmjd = xmjd0 call csrtide (xmjd,dlat,dlon,tide,tlp,tds,isdata, . u, v, pseudo, radiation) call admit2 (pseudo,radiation,u,v,x,y,h1,h2,amp,pha) write (6,*) write (6,*) ' Latitude = ' , dlat, ' (deg) ', . ' Longitude = ', dlon, ' (deg)' if (isdata) then write (6,102) x(25),y(25) ! admittance at M2 write (6,103) h1(25),h2(25) ! in-phase and quadrature components write (6,104) amp(25),pha(25) ! amplitude and phase else write (6,*) ' Selected location not in tide model domain' endif write (6,*) c...illustrate use of csrtide to synthesize the tide at a given time,lat,lon c tide is the total tide in mm, tlp is the long period and tds is the sum c of the diurnal and semidiurnal tides in mm. tide=tlp+tds do i=1,nmjd xmjd = xmjd0 + (i-1)*dmjd jd = xmjd + 2400000.5d0 ! Julian date from modified julian date call kalday (jd,mdyhms) ! Convert to calendar date call csrtide (xmjd,dlat,dlon,tide,tlp,tds,isdata, . u,v,pseudo,radiation) if (isdata) then if (i.eq.1) write (6,101) write (6,100) xmjd,mdyhms,dlat,dlon,tide,tlp,tds else write (6,*) ' Selected location not in tide model domain' go to 20 endif enddo write (6,*) go to 20 100 format (1x,f10.3,1x,i2.2,'/',i2.2,'/',i2.2,1x,i2.2,':',i2.2, . ':',i2.2,3x,2f10.2,3f13.0) 101 format (/, ' MJD (days) Calender Date', . ' Lat (deg) Lon (deg) Total Tide (mm) ', . ' Long Per. ', ' Diur.+Semid.') 102 format (' Admittance at M2 (Real, Imag) = ', 2f8.3) 103 format (' In-phase and Quadrature Components at M2 (mm) = ', . 2f8.0) 104 format (' Amplitude (mm) and Phase (deg) at M2 = ', f8.0, f8.2) end 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...Note: An "orthoweight" file can contain geocentric tide, ocean tide, c load tide, or tide model difference orthoweights. If the c orthoweight file in use contains differences or loads then c the user should beware of using the "tide" return variable c and should use "tds" instead. 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 Version 3.0 - March 1995, new version of tptide allows compact integer c format for orthoweight file 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 using c subroutine admit2. 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 is based upon 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, -90 to 90) * * dlon d i east longitude (in degrees, 0 to 360). * * 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 user supplied orthoweight file, 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 iunit - global array of orthoweights defining tidal * constants at all valid locations. (see parameter statement) * * 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 compiler option to ignore double * precision should be selected. * * 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) * 3.0 1/24/95 r. eanes parameter statements for dimensions and unit * number. compact format option for orthoweight file. * increased dimensions to allow 0.5 x 0.5 deg grid. * changed wmin to 0.2 to make rsc94b work at gauges. * implicit double precision (a-h,o-z) parameter (nxmax=720,nymax=361) parameter (iunit=51) character*80 title1,title2 real uv(2,3,2,nxmax,nymax), undef double precision latmin,latmax,lonmin,lonmax dimension u(3,2),v(3,2),iuv(12) save logical init,isdata,pseudo,radiation data init/.true./, undef/99999./, eps /1.d-10/, wmin /0.2/ * on first call, read t/p -derived ortho-weights * ------------------------------------------------ if (init) then init = .false. do 4 j=1,nymax do 4 i=1,nxmax uv(1,1,1,i,j) = undef 4 continue read (iunit,14) title1,title2 write (6,14) title1,title2 read (iunit,11) nx,ny,latmin,latmax,lonmin,lonmax,idiff, . itype,irad,ifmt c c...orthoweight file format options c idiff - if .eq. 1 then zero tide returned for points outside of domain c and isdata set to true c itype - if .eq. 1 then cartwright&ray 1991 pseudo-orthoweights are assumed c irad - if .eq. 1 then the cartwright&ray radiation potential modification is used c ifmt - if .eq. 1 then compact integer format used for orthoweights c 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.nxmax .or. ny.gt.nymax) 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 if (ifmt.eq.0) then ! original c&r format read(iunit,12,end=10) i,j,u1,v1,u2,v2,u3,v3, . 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 endif if (ifmt.eq.1) then ! compact integer format read(iunit,*,end=10) i,j,iuv uv(1,1,1,i,j) = iuv(1)/1000.d0 uv(2,1,1,i,j) = iuv(2)/1000.d0 uv(1,2,1,i,j) = iuv(3)/1000.d0 uv(2,2,1,i,j) = iuv(4)/1000.d0 uv(1,3,1,i,j) = iuv(5)/1000.d0 uv(2,3,1,i,j) = iuv(6)/1000.d0 uv(1,1,2,i,j) = iuv(7)/1000.d0 uv(2,1,2,i,j) = iuv(8)/1000.d0 uv(1,2,2,i,j) = iuv(9)/1000.d0 uv(2,2,2,i,j) = iuv(10)/1000.d0 uv(1,3,2,i,j) = iuv(11)/1000.d0 uv(2,3,2,i,j) = iuv(12)/1000.d0 endif go to 8 10 continue 11 format(2i4,4f8.3,4i2) 12 format(2i4,16x,6f8.3,/,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 if (ilon1.gt.nx) ilon1 = ilon1 - 1 ! avoid problem for dlon=360. ilon2 = ilon1 + 1 if (ilon2.gt.nx) ilon2 = 1 ! longitude grid must span globe if (ilon1.lt.1 .or. ilon1.gt.nx) then write (6,*) 'error in tptide: incorrect longitude input' write (6,*) 'mjd,dlon,xlon,ilon1,ilon2', . time/86400.d0,dlon,xlon,ilon1,ilon2 stop 301 endif 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 subroutine kalday (utc,mdyhms) c implicit double precision (a-h,o-z) dimension mdyhms(6),iptime(5) c c...compute the hour, minute, and second. c t = utc 5 jd = t xtime = mod((t-jd)*24.0d0+12.0d0,24.0d0) iptime(4) = xtime xtime = (xtime-iptime(4))*60.0d0 iptime(5) = xtime ptimes = (xtime-iptime(5))*60.0d0 c c...compute the year, month, and day. c if (iptime(4).lt.12) jd = jd+1 lx = jd+68569 nx = 4*lx/146097 lx = lx-(146097*nx+3)/4 iptime(3) = 4000*(lx+1)/1461001 lx = lx-1461*iptime(3)/4+31 iptime(1) = 80*lx/2447 iptime(2) = lx-2447*iptime(1)/80 lx = iptime(1)/11 iptime(1) = iptime(1)+2-12*lx iptime(3) = 100*(nx-49)+iptime(3)+lx do 10 i = 1,5 mdyhms(i) = iptime(i) 10 continue mdyhms(3) = mdyhms(3)-1900 mdyhms(6) = ptimes+0.5d0 c if (mdyhms(4).eq.23 .and. mdyhms(5).eq.59 .and. . mdyhms(6).eq.60) then t = t + 1.0d-6 go to 5 endif c return c end subroutine julday (itimes,stime,utc) c c...purpose: to convert between calendar day and julian date (utc). c c output: c utc double precision julian date c c input: c itimes integer array containing month, day, year, hour, minute c stime floating point seconds c implicit double precision (a-h,o-z) dimension itimes(5) c jd = itimes(2)-32075+1461*(itimes(3)+4800+(itimes(1)-14)/12)/4+367 * *(itimes(1)-2-(itimes(1)-14)/12*12)/12-3*((itimes(3)+4900+ * (itimes(1)-14)/12)/100)/4 c utc = dble(float(jd))-0.5d+0+dble(float(itimes(4)))/24.0d+0+ * dble(float(itimes(5)))/1440.0d+0+dble(stime)/86400.0d+0 c return c c.....end julday c end