subroutine BiasTide(doodson,bias) real*8::doodson,bias integer i * ------------------------------------- bias=0.d0 !(n-m)*90 if(dabs(doodson- 55565)<1.d-3)bias=180 ! omiga1 if(dabs(doodson-125755)<1.d-3)bias=-90 ! 2Q1 if(dabs(doodson-127555)<1.d-3)bias=-90 ! sigma1 if(dabs(doodson-135655)<1.d-3)bias=-90 ! Q1 if(dabs(doodson-137455)<1.d-3)bias=-90 ! rou1 if(dabs(doodson-145545)<1.d-3)bias=-90 ! rou1 if(dabs(doodson-145555)<1.d-3)bias=-90 ! O1 if(dabs(doodson-155555)<1.d-3)bias=-90 ! M1 if(dabs(doodson-155655)<1.d-3)bias=90 ! M1 if(dabs(doodson-157455)<1.d-3)bias=90 ! kama1 if(dabs(doodson-162556)<1.d-3)bias=-90 ! pi1 if(dabs(doodson-163555)<1.d-3)bias=-90 ! P1 if(dabs(doodson-163655)<1.d-3)bias=90 ! theta1 if(dabs(doodson-164556)<1.d-3)bias=90 ! S1 if(dabs(doodson-165555)<1.d-3)bias=90 ! K1 if(dabs(doodson-166554)<1.d-3)bias=90 ! psi1 if(dabs(doodson-167555)<1.d-3)bias=90 ! phi1 if(dabs(doodson-175455)<1.d-3)bias=90 ! J1 if(dabs(doodson-185555)<1.d-3)bias=90 ! OO1 if(dabs(doodson-253755)<1.d-3)bias=180 ! gama2 if(dabs(doodson-255545)<1.d-3)bias=180 ! alfa2 if(dabs(doodson-265455)<1.d-3)bias=180 ! L2 if(dabs(doodson-263655)<1.d-3)bias=180 ! namta2 if(dabs(doodson-275554)<1.d-3)bias=180 ! R2 if(dabs(doodson-355555)<1.d-3)bias=180 ! M3 end ! !********************************************************* ! subroutine CalcTidefu(mjd,doodson,df,du) real*8::mjd,doodson,df,du,pi,RAD real*8::shpn(5),td,f(22),u(22),ddsn(22),SINN,SIN2N,COSN,COS2N,OMEGA integer i * ------------------------------------- pi=datan(1.d0)*4.d0;RAD=pi/180.d0;df=1.d0;du=0.d0 CALL ASTRO5(mjd, shpn ) OMEGA=shpn(4) ddsn(1)=135655! Q1 ddsn(2)=145555! O1 ddsn(3)=163555! P1 ddsn(4)=165555! K1 ddsn(5)=245655! N2 ddsn(6)=255555! M2 ddsn(7)=273555! S2 ddsn(8)=275555! K2 ddsn(9)=125755! 2Q1 ddsn(10)=127555! sigma1 ddsn(11)=137455! rho1 ddsn(12)=155655! M1 ddsn(13)=175455! J1 ddsn(14)=185555! OO1 ddsn(15)=235755! 2N2 ddsn(16)=237555! mu2 ddsn(17)=247455! nu2 ddsn(18)=265455! L2 ddsn(19)=164556! S1 (Doodson's phase) ddsn(20)=355555! M3 ddsn(21)=455555! M4 ddsn(22)=655555! M6 * * determine nodal corrections f and u * Note: Update this code next iteration of model -RDR * ----------------------------------- SINN = DSIN(OMEGA*RAD) COSN = DCOS(OMEGA*RAD) SIN2N = DSIN(2.d0*OMEGA*RAD) COS2N = DCOS(2.d0*OMEGA*RAD) F(1) = 1.009 + 0.187*COSN - 0.015*COS2N F(2) = F(1) F(4) = 1.006 + 0.115*COSN - 0.009*COS2N F(3) = F(4) F(5) = 1.000 - 0.037*COSN F(6) = F(5) F(8) = 1.024 + 0.286*COSN + 0.008*COS2N F(7) = F(8) F(9) = DSQRT((1.0 + 0.189*COSN - 0.0058*COS2N)**2 + * (0.189*SINN - 0.0058*SIN2N)**2) F(10) = F(9) F(11) = F(9) F(12) = DSQRT((1.0 + 0.185*COSN)**2 + (0.185*SINN)**2) F(13) = DSQRT((1.0 + 0.198*COSN)**2 + (0.198*SINN)**2) F(14) = DSQRT((1.0 + 0.640*COSN + 0.134*COS2N)**2 + * (0.640*SINN + 0.134*SIN2N)**2 ) F(15) = DSQRT((1.0 - 0.0373*COSN)**2 + (0.0373*SINN)**2) F(16) = F(15) F(17) = F(15) F(18) = F(15) F(19) = F(3) F(20) = F(6)*dsqrt(F(6)) F(21) = F(6)*F(6) F(22) = F(6)*F(6)*F(6) U(1) = 10.8*SINN - 1.3*SIN2N U(2) = U(1) U(4) = -8.9*SINN + 0.7*SIN2N U(3) = U(4) U(5) = -2.1*SINN U(6) = U(5) U(8) = -17.7*SINN + 0.7*SIN2N U(7) = U(8) U(9) = DATAN2(0.189*SINN - 0.0058*SIN2N, * 1.0 + 0.189*COSN - 0.0058*SIN2N)/RAD U(10) = U(9) U(11) = U(9) U(12) = DATAN2( 0.185*SINN, 1.0 + 0.185*COSN)/RAD U(13) = DATAN2(-0.198*SINN, 1.0 + 0.198*COSN)/RAD U(14) = DATAN2(-0.640*SINN - 0.134*SIN2N, * 1.0 + 0.640*COSN + 0.134*COS2N)/RAD U(15) = DATAN2(-0.0373*SINN, 1.0 - 0.0373*COSN)/RAD U(16) = U(15) U(17) = U(15) U(18) = U(15) U(19) = U(3) U(20) = 1.5d0*U(6) U(21) = 2.d0*U(6) U(22) = 3.d0*U(6) do i=1,22 if(dabs(doodson-ddsn(i))<0.1)then df=f(i);du=u(i);return endif enddo return end ! !********************************************************************************* ! SUBROUTINE ASTRO5( TIME, SHPNP ) * *--------------------------------------------------------------------- * Computes the 5 basic astronomical mean longitudes s, h, p, N, p'. * * Note N is not N', i.e. N is decreasing with time. * * TIME is UTC in decimal Modified Julian Day (MJD). * All longitudes returned in degrees. * * R. D. Ray, NASA/GSFC August 2003 * * Most of the formulae for mean longitudes are extracted from * Jean Meeus, Astronomical Algorithms, 2nd ed., 1998. * Page numbers below refer to this book. * * Note: This routine uses TIME in UT and does not distinguish between * the subtle differences of UTC, UT1, etc. This is more than adequate * for the calculation of these arguments, especially in tidal studies. *--------------------------------------------------------------------- * IMPLICIT NONE DOUBLE PRECISION TIME, SHPNP(5) DOUBLE PRECISION TJD,T,CIRCLE,DEL,TJLAST,D PARAMETER (CIRCLE=360.0D0) * Convert to Julian Day and to Ephemeris Time * ------------------------------------------- TJD = TIME + 2400000.5D0 * Compute time argument in centuries relative to J2000 * ---------------------------------------------------- T = ( TJD - 2451545.d0 )/36525.d0 * * * mean longitude of moon (p.338) * ------------------------------ SHPNP(1) = (((-1.53388d-8*T + 1.855835d-6)*T - 1.5786d-3)*T + * 481267.88123421d0)*T + 218.3164477d0 * * mean elongation of moon (p.338) * ------------------------------- D = (((-8.8445d-9*T + 1.83195d-6)*T - 1.8819d-3)*T + * 445267.1114034d0)*T + 297.8501921d0 * * mean longitude of sun * --------------------- SHPNP(2) = SHPNP(1) - D * * mean longitude of lunar perigee (p.343) * --------------------------------------- SHPNP(3) = ((-1.249172d-5*T - 1.032d-2)*T + 4069.0137287d0)*T + * 83.3532465d0 * * mean longitude of ascending lunar node (p.144) * ---------------------------------------------- SHPNP(4) = ((2.22222d-6*T + 2.0708d-3)*T - 1934.136261d0)*T + * 125.04452d0 * * mean longitude of solar perigee (Simon et al., 1994) * ---------------------------------------------------- SHPNP(5) = 282.94d0 + 1.7192d0 * T SHPNP(1) = DMOD(SHPNP(1),CIRCLE) SHPNP(2) = DMOD(SHPNP(2),CIRCLE) SHPNP(3) = DMOD(SHPNP(3),CIRCLE) SHPNP(4) = DMOD(SHPNP(4),CIRCLE) IF (SHPNP(1).LT.0.D0) SHPNP(1) = SHPNP(1) + CIRCLE IF (SHPNP(2).LT.0.D0) SHPNP(2) = SHPNP(2) + CIRCLE IF (SHPNP(3).LT.0.D0) SHPNP(3) = SHPNP(3) + CIRCLE IF (SHPNP(4).LT.0.D0) SHPNP(4) = SHPNP(4) + CIRCLE RETURN END