C**** TIDES.FOR 2004/10/18 C**** Gravitational perturbations due to Moon and Sun at given time C**** IMPLICIT REAL*8 (A-H,M,O-Z) PARAMETER (IM=90,JM=60, * TWOPI = 6.283185307179586477d0, * GUNIV = 6.672d-11, ! universal gravity coef. (m^3/kg*s^2) * GRAV = 9.807d0, ! mean gravity acceleration due to Earth * RADIUS = 6371000, ! mean radius of Earth (m) * ASEMI = 1.4960d11, ! semi-major axis of Earth's orbit (m) * MEARTH = 5.976d24, ! mass of Earth (kg) * MMOON = 7.35028d22, ! mass of Moon (kg) * MSUN = 1.98941d30, ! mass of Sun (kg) * ECCEN = .0167d0, ! eccentricity of Earth's orbit for 2000 * OBLIQD = 23.44d0, ! obliquity of Earth's orgit for 2000 * OMEGVP = 4.938d0) ! longitude of perihelion for year 2000 CHARACTER*80 ARGMNT, TITLE(3) INTEGER*4 YMDHtoT REAL*4 UG(IM,JM),VG(IM,JM),WG(IM,JM), * UM(IM,JM),VM(IM,JM),WM(IM,JM), * US(IM,JM),VS(IM,JM),WS(IM,JM) COMMON /GEOMCB/ SINI(IM),COSI(IM), SINJ(JM),COSJ(JM) C**** DATA TITLE / 1'EASTWARD TIDES (10^-6 m/s^2) Moon: , Sun: ,', 2'NORTHWARD TIDES (10^-6 m/s^2) Moon: , Sun: ,', 3'VERTICAL TIDES (10^-6 m/s^2) Moon: , Sun: ,'/ C**** C**** X, Y and Z use Earth's rectangular geocentric coordinates C**** U, V and W use Earth's local topographic coordinates C**** NARGS = IARGC() IF(NARGS < 1) GO TO 800 CALL SGEOM C**** C**** Calculate desired DAY from command line arguments C**** JMONT = 1 JDATE = 1 JHOUR = 0 JMINU = 0 CALL GETARG (1,ARGMNT) READ (ARGMNT,*,ERR=801) JYEAR IF(JYEAR < -240000 .or. JYEAR > 240000) GO TO 801 IF(NARGS==1) GO TO 10 CALL GETARG (2,ARGMNT) READ (ARGMNT,*,ERR=802) JMONT IF(JMONT < 1 .or. JMONT > 12) GO TO 802 IF(NARGS==2) GO TO 10 CALL GETARG (3,ARGMNT) READ (ARGMNT,*,ERR=803) JDATE IF(JDATE < 1 .or. JDATE > 31) GO TO 803 IF(NARGS==3) GO TO 10 CALL GETARG (4,ARGMNT) READ (ARGMNT,*,ERR=804) JHOUR IF(JHOUR < 0 .or. JHOUR > 23) GO TO 804 IF(NARGS==4) GO TO 10 CALL GETARG (5,ARGMNT) READ (ARGMNT,*,ERR=805) JMINU IF(JMINU < 0 .or. JMINU > 59) GO TO 805 10 JTIME = YMDHtoT (JYEAR,JMONT,JDATE,JHOUR) DAY = (JTIME + JMINU/60d0) / 24d0 C**** C**** Locate Moon and Sun for given day C**** CALL MOON (DAY, MOOLON,MOOLAT,MOODIS) OBLIQ = OBLIQD*TWOPI/360 ! convert obliquity to radians CALL ORBIT (ECCEN,OBLIQ,OMEGVP, DAY, * DIST,SIND,COSD,SUNLON,SUNLAT,EQTIME) WRITE (TITLE(1)(37:40),901) NINT(MOOLON*360/TWOPI) WRITE (TITLE(1)(42:44),902) NINT(MOOLAT*360/TWOPI) WRITE (TITLE(1)(52:55),901) NINT(SUNLON*360/TWOPI) WRITE (TITLE(1)(57:59),902) NINT(SUNLAT*360/TWOPI) C**** Calculate distances between centers (m) DMCmEC = RADIUS*MOODIS ! distance between Moon and Earth DSCmEC = ASEMI*DIST ! distance between Sun and Earth C**** Calculate location of Moon and Sun XMC = DMCmEC*COS(MOOLAT)*COS(MOOLON) YMC = DMCmEC*COS(MOOLAT)*SIN(MOOLON) ZMC = DMCmEC*SIN(MOOLAT) XSC = DSCmEC*COS(SUNLAT)*COS(SUNLON) YSC = DSCmEC*COS(SUNLAT)*SIN(SUNLON) ZSC = DSCmEC*SIN(SUNLAT) C**** C**** Loop over Model grid cells C**** DO 20 J=1,JM DO 20 I=1,IM C**** Calculate location of cell on Earth's surface XES = RADIUS*COSJ(J)*COSI(I) YES = RADIUS*COSJ(J)*SINI(I) ZES = RADIUS*SINJ(J) C**** Calculate spatial perturbation of gravity due to Moon and Sun DMCmES = SQRT ((XMC-XES)**2 + (YMC-YES)**2 + (ZMC-ZES)**2) XGM = GUNIV*MMOON*((XMC-XES)/DMCmES**3 - XMC/DMCmEC**3) YGM = GUNIV*MMOON*((YMC-YES)/DMCmES**3 - YMC/DMCmEC**3) ZGM = GUNIV*MMOON*((ZMC-ZES)/DMCmES**3 - ZMC/DMCmEC**3) DSCmES = SQRT ((XSC-XES)**2 + (YSC-YES)**2 + (ZSC-ZES)**2) XGS = GUNIV*MSUN*((XSC-XES)/DSCmES**3 - XSC/DSCmEC**3) YGS = GUNIV*MSUN*((YSC-YES)/DSCmES**3 - YSC/DSCmEC**3) ZGS = GUNIV*MSUN*((ZSC-ZES)/DSCmES**3 - ZSC/DSCmEC**3) C**** Calculate unit vectors of local topographic coordinates in C**** the rectangular geocentric coordinate system XU = - SINI(I) YU = COSI(I) ZU = 0 XV = - SINJ(J)*COSI(I) YV = - SINJ(J)*SINI(I) ZV = COSJ(J) XW = COSJ(J)*COSI(I) YW = COSJ(J)*SINI(I) ZW = SINJ(J) C**** Project the gravity pertubations onto the U, V and W unit vectors UM(I,J) = (XU*XGM + YU*YGM + ZU*ZGM) * 1d6 VM(I,J) = (XV*XGM + YV*YGM + ZV*ZGM) * 1d6 WM(I,J) = (XW*XGM + YW*YGM + ZW*ZGM) * 1d6 US(I,J) = (XU*XGS + YU*YGS + ZU*ZGS) * 1d6 VS(I,J) = (XV*XGS + YV*YGS + ZV*ZGS) * 1d6 WS(I,J) = (XW*XGS + YW*YGS + ZW*ZGS) * 1d6 UG(I,J) = (XU*(XGM+XGS) + YU*(YGM+YGS) + ZU*(ZGM+ZGS)) * 1d6 VG(I,J) = (XV*(XGM+XGS) + YV*(YGM+YGS) + ZV*(ZGM+ZGS)) * 1d6 20 WG(I,J) = (XW*(XGM+XGS) + YW*(YGM+YGS) + ZW*(ZGM+ZGS)) * 1d6 C**** C**** Write gravity perturbations (tides) to disk C**** OPEN (2,FILE='TIDES.O',FORM='UNFORMATTED') WRITE (TITLE(1)(61:80),903) JYEAR,JMONT,JDATE,JHOUR,JMINU TITLE(2)(32:80) = TITLE(1)(32:80) TITLE(3)(32:80) = TITLE(1)(32:80) WRITE (2) TITLE(1),UG WRITE (2) TITLE(2),VG WRITE (2) TITLE(3),WG ARGMNT = TITLE(1) TITLE(1)(46:60) = ' ' TITLE(2)(46:60) = ' ' TITLE(3)(46:60) = ' ' WRITE (2) TITLE(1),UM WRITE (2) TITLE(2),VM WRITE (2) TITLE(3),WM TITLE(1)(31:45) = ARGMNT(47:61) TITLE(2)(31:45) = ARGMNT(47:61) TITLE(3)(31:45) = ARGMNT(47:61) WRITE (2) TITLE(1),US WRITE (2) TITLE(2),VS WRITE (2) TITLE(3),WS CLOSE (2) GO TO 999 C**** 800 WRITE (6,*) *'Usage: TIDES Year [Month [Date [Hour [Minute]]]] 2004/10/18' WRITE (6,*) *' Tidal accelerations for local topographic coordinates.' GO TO 999 801 WRITE (0,*) 'Integer year between -240000 and 240000 is' // * ' required, not: ',ARGMNT STOP 801 802 WRITE (0,*) 'Integer month between 1 and 12 is required, not: ', * ARGMNT STOP 802 803 WRITE (0,*) 'Integer date between 1 and 31 is required, not: ', * ARGMNT STOP 803 804 WRITE (0,*) 'Integer hour between 0 and 23 is required, not: ', * ARGMNT STOP 804 805 WRITE (0,*) 'Integer minute between 0 and 59 is required, not: ', * ARGMNT STOP 805 C**** 901 FORMAT (I4) 902 FORMAT (I3) 903 FORMAT (I8,3('/',I2.2),':',I2.2) 999 END SUBROUTINE SGEOM C**** C**** SGEOM calculates the Spherical GEOMetry parameters for the C-grid C**** and fills in the GEOMCB common block C**** IMPLICIT REAL*8 (A-H,M,O-Z) PARAMETER (IM=90,JM=60, * TWOPI = 6.283185307179586477d0, * RADIUS = 6371000) ! mean radius of Earth (m) COMMON /GEOMCB/ SINI(IM),COSI(IM), SINJ(JM),COSJ(JM) C**** DO 10 I=1,IM SINI(I) = SIN ((I-.5-IM/2)*TWOPI/IM) 10 COSI(I) = COS ((I-.5-IM/2)*TWOPI/IM) FJEQ = .5*(1+JM) DO 20 J=1,JM SINJ(J) = SIN ((J-FJEQ)*TWOPI*.5/JM) 20 COSJ(J) = COS ((J-FJEQ)*TWOPI*.5/JM) RETURN END FUNCTION YMDHtoT (JYEAR,JMONT,JDATE,JHOUR) C**** C**** For a given JYEAR (A.D.), JMONT, JDATE, and JHOUR, calculate the C**** number of HOURs measured from 2000 January 1, hour 0 using the C**** Gregorian calendar. C**** INTEGER*4, PARAMETER :: * JDofMN(0:12)=(/0,31,59,90,120,151,181,212,243,273,304,334,365/), * JDofML(0:12)=(/0,31,60,91,121,152,182,213,244,274,305,335,366/) INTEGER*4 YMDHtoT PARAMETER (KDAY4C = 365*400 + 97, ! number of days in 4 centuries * KDAY1C = 365*100 + 24, ! number of days in 1 century * KDAY4Y = 365* 4 + 1, ! number of days in 4 years * KDAY1Y = 365) ! number of days in 1 year C**** IF(ABS(JYEAR).gt.240000) GO TO 800 N4CENT = FLOOR((JYEAR-2000.d0)/400.) IYR4C = JYEAR-2000 - N4CENT*400 N1CENT = IYR4C/100 IYR1C = IYR4C - N1CENT*100 N4YEAR = IYR1C/4 IYR4Y = IYR1C - N4YEAR*4 N1YEAR = IYR4Y IDAY = N4CENT*KDAY4C IF(N1CENT.gt.0) GO TO 10 C**** First of every fourth century: 16??, 20??, 24??, etc. IDAY = IDAY + N4YEAR*KDAY4Y IF(N1YEAR.gt.0) GO TO 200 GO TO 100 C**** Second to fourth of every fourth century: 21??, 22??, 23??, etc. 10 IDAY = IDAY + KDAY1C+1 + (N1CENT-1)*KDAY1C IF(N4YEAR.gt.0) GO TO 20 C**** First four years of every second to fourth century when there is C**** no leap year during the four years: 2100-2103, 2200-2203, etc. IDAY = IDAY + N1YEAR*KDAY1Y GO TO 210 C**** Subsequent four years of every second to fourth century when C**** there is a leap year: 2104-2107, 2108-2111 ... 2204-2207, etc. 20 IDAY = IDAY + KDAY4Y-1 + (N4YEAR-1)*KDAY4Y IF(N1YEAR.gt.0) GO TO 200 C**** C**** Current year is a leap year C**** 100 IDAY = IDAY + JDofML(JMONT-1) + JDATE-1 YMDHtoT = IDAY*24 + JHOUR RETURN C**** C**** Current year is not a leap year C**** 200 IDAY = IDAY + KDAY1Y+1 + (N1YEAR-1)*KDAY1Y 210 IDAY = IDAY + JDofMN(JMONT-1) + JDATE-1 YMDHtoT = IDAY*24 + JHOUR RETURN C**** 800 WRITE (6,*) 'JYEAR exceeds output limits of YMDHtoT. JYEAR =', * JYEAR STOP 800 END SUBROUTINE MOON (DAY, MOOLON,MOOLAT,MOODIS) C**** C**** Input: DAY = days measured since 2000 January 1, hour 0 C**** C**** Output: MOOLON = Earth's LONgitude directly beneath the MOOn C**** MOOLAT = Earth's LATitude directly beneath the MOOn C**** MOODIS = DIStance from MOOn to Earth in Earth radii units C**** C**** Constants: OBLIQ = dihedral angle between Earth's equatorial C**** plane and orbital plane of year 2000 C**** TAofVE = true anomaly of vernal equinox of year 2000 C**** MAofVE = mean anomaly of vernal equinox of year 2000 C**** EDAYzY = tropical year = Earth days per year = 365.2425 C**** VE200D = days from 2000 January 1, hour 0 till vernal C**** equinox of year 2000 = 31 + 29 + 19 + 7.5/24 C**** C**** Intermediate quantities: C**** RA = right ascension of the Moon = C**** = longitude of Moon in geocentric equatorial coordinates C**** VEQLON = longitude of Greenwich Meridion in geocentric C**** equatorial coordinates at vernal equinox of year 2000 C**** ROTATE = change in longitude in geocentric nonrotating C**** equatorial coordinates from a point's location on C**** vernal equinox to is current location where the point C**** is fixed on rotating Earth C**** C**** Angles, longitude and latitude are measured in radians unless C**** the variable is terminated by a capital D in which case it is C**** measured in degrees C**** IMPLICIT REAL*8 (A-H,M,O-Z) PARAMETER (TWOPI=6.283185307179586477d0, PIz180=TWOPI/360., * OBLIQ=.409101d0, TAofVE=1.345728d0, MAofVE=1.313255d0, * EDAYzY=365.2425d0, VE200D=79.3125) INTEGER*4 RADEG,RAMIN, MLDEG,MLMIN REAL*8 ECLONN(3,0:6),ECLATN(3,4),HOPARN(3,0:4) C**** DATA ECLONN /218.32d0, 0. , 481267.883d0, * 6.29d0, 134.9d0, 477198.85d0, * -1.27d0, 259.2d0, -413335.38d0, * .66d0, 235.7d0, 890534.23d0, * .21d0, 269.9d0, 954397.70d0, * -.19d0, 357.5d0, 35999.05d0, * -.11d0, 186.6d0, 966404.05d0/ DATA ECLATN / 5.13d0, 93.3d0, 483202.03d0, * .28d0, 228.2d0, 960400.87d0, * -.28d0, 318.3d0, 6003.18d0, * -.17d0, 217.6d0, -407332.20d0/ DATA HOPARN / .9508d0, 0. , 0. , * .0518d0, 134.9d0, 477198.85d0, * .0095d0, 259.2d0, -413335.38d0, * .0078d0, 235.7d0, 890534.23d0, * .0028d0, 269.9d0, 954397.70d0/ C**** Calculate time in Julian centuries TJC = (DAY-.5) / (100.*365.25) C**** C**** Geocentric Ecliptic coordinates: C**** EcLat = angle between an object and its perpendicular projection C**** onto the Ecliptic plane with Earth center as angle vertex C**** EcLon = angle from vernal equinox vector (ray from Earth to Sun C**** at vernal equinox) to projection onto the ecliptic plane C**** of vector of current Earth to object C**** C**** Calculate ecliptic longitude ECLOND in degrees and radians ECLON ECLOND = ECLONN(1,0) + ECLONN(3,0)*TJC DO 10 I=1,6 10 ECLOND = ECLOND + + ECLONN(1,I)*SIN(PIz180*(ECLONN(2,I)+ECLONN(3,I)*TJC)) ECLOND = MOD (ECLOND, 360.d0) ECLON = ECLOND*PIz180 C**** Calculate ecliptic latitude ECLATD in degrees and radians ECLAT ECLATD = 0. DO 20 I=1,4 20 ECLATD = ECLATD + + ECLATN(1,I)*SIN(PIz180*(ECLATN(2,I)+ECLATN(3,I)*TJC)) ECLAT = ECLATD*PIz180 C**** Calculate horizontal parallax in degrees HOPARD and radians HOPAR HOPARD = HOPARN(1,0) DO 30 I=1,4 30 HOPARD = HOPARD + + HOPARN(1,I)*COS(PIz180*(HOPARN(2,I)+HOPARN(3,I)*TJC)) HOPAR = HOPARD*PIz180 C**** Calculate MOODIS in Earth radii MOODIS = 1. / SIN(HOPAR) C**** Calculate geocentric ecliptic rectangular coordinates XEC = COS(ECLAT)*COS(ECLON) YEC = COS(ECLAT)*SIN(ECLON) ZEC = SIN(ECLAT) C**** C**** Geocentric nonrotating Equatorial coordinates: C**** Declination = angle between an object and its perpendicular C**** projection onto the Equatorial plane with Earth C**** center as angle vertex C**** RightAscention = angle from the vernal equinox vector (ray from C**** Earth to Sun at vernal equinox) to projection C**** onto equatorial plane of vector of current Earth C**** to object C**** C**** Calculate geocentric equatorial rectangular coordinates XEQ = XEC YEQ = YEC*COS(OBLIQ) - ZEC*SIN(OBLIQ) ZEQ = ZEC*COS(OBLIQ) + YEC*SIN(OBLIQ) C**** Calculate right ascension and MOOLON RA = ATAN2 (YEQ,XEQ) VEQLON = TWOPI*VE200D - TWOPI/2. + MAofVE - TAofVE ! modulo 2*Pi ROTATE = TWOPI*(DAY-VE200D)*(EDAYzY+1.)/EDAYzY MOOLON = MODULO (RA-ROTATE-VEQLON, TWOPI) IF(MOOLON.gt.TWOPI/2.) MOOLON = MOOLON - TWOPI C**** Calculate declination = MLAT in radians MOOLAT = ASIN (ZEQ) RETURN END SUBROUTINE ORBIT (ECCEN,OBLIQ,OMEGVP, DAY, * DIST,SIND,COSD,SUNLON,SUNLAT,EQTIME) C**** C**** ORBIT receives orbital parameters and time of year, and returns C**** distance from Sun and Sun's position. C**** Reference for following caculations is: V.M.Blanco and C**** S.W.McCuskey, 1961, "Basic Physics of the Solar System", pages C**** 135 - 151. Existence of Moon and heavenly bodies other than C**** Earth and Sun are ignored. Earth is assumed to be spherical. C**** C**** Program author: Gary L. Russell 2002/09/25 C**** Angles, longitude and latitude are measured in radians. C**** C**** Input: ECCEN = eccentricity of the orbital ellipse C**** OBLIQ = latitude of Tropic of Cancer C**** OMEGVP = longitude of perihelion (sometimes Pi is added) = C**** = spatial angle from vernal equinox to perihelion C**** with Sun as angle vertex C**** DAY = days measured since 2000 January 1, hour 0 C**** C**** Constants: EDAYzY = tropical year = Earth days per year = 365.2425 C**** VE200D = days from 2000 January 1, hour 0 till vernal C**** equinox of year 2000 = 31 + 29 + 19 + 7.5/24 C**** C**** Intermediate quantities: C**** BSEMI = semi minor axis in units of semi major axis C**** PERIHE = perihelion in days since 2000 January 1, hour 0 C**** in its annual revolution about Sun C**** TA = true anomaly = spatial angle from perihelion to C**** current location with Sun as angle vertex C**** EA = eccentric anomaly = spatial angle measured along C**** eccentric circle (that circumscribes Earth's orbit) C**** from perihelion to point above (or below) Earth's C**** absisca (where absisca is directed from center of C**** eccentric circle to perihelion) C**** MA = mean anomaly = temporal angle from perihelion to C**** current time in units of 2*Pi per tropical year C**** TAofVE = TA(VE) = true anomaly of vernal equinox = - OMEGVP C**** EAofVE = EA(VE) = eccentric anomaly of vernal equinox C**** MAofVE = MA(VE) = mean anomaly of vernal equinox C**** SLNORO = longitude of Sun in Earth's nonrotating reference frame C**** VEQLON = longitude of Greenwich Meridion in Earth's nonrotating C**** reference frame at vernal equinox C**** ROTATE = change in longitude in Earth's nonrotating reference C**** frame from point's location on vernal equinox to its C**** current location where point is fixed on rotating Earth C**** SLMEAN = longitude of fictitious mean Sun in Earth's rotating C**** reference frame (normal longitude and latitude) C**** C**** Output: DIST = distance to Sun in units of semi major axis C**** SIND = sine of declination angle = sin(SUNLAT) C**** COSD = cosine of the declination angle = cos(SUNLAT) C**** SUNLON = longitude of point on Earth directly beneath Sun C**** SUNLAT = latitude of point on Earth directly beneath Sun C**** EQTIME = Equation of Time = C**** = longitude of fictitious mean Sun minus SUNLON C**** C**** From the above reference: C**** (4-54): [1 - ECCEN*cos(EA)]*[1 + ECCEN*cos(TA)] = (1 - ECCEN^2) C**** (4-55): tan(TA/2) = sqrt[(1+ECCEN)/(1-ECCEN)]*tan(EA/2) C**** Yield: tan(EA) = sin(TA)*sqrt(1-ECCEN^2) / [cos(TA) + ECCEN] C**** or: tan(TA) = sin(EA)*sqrt(1-ECCEN^2) / [cos(EA) - ECCEN] C**** IMPLICIT REAL*8 (A-H,M-Z) PARAMETER (TWOPI=6.283185307179586477d0, * EDAYzY=365.2425d0, VE200D=79.3125) C**** C**** Determine EAofVE from geometry: tan(EA) = b*sin(TA) / [e+cos(TA)] C**** Determine MAofVE from Kepler's equation: MA = EA - e*sin(EA) C**** Determine MA knowing time from vernal equinox to current day C**** BSEMI = SQRT (1. - ECCEN*ECCEN) TAofVE = - OMEGVP EAofVE = ATAN2 (BSEMI*SIN(TAofVE), ECCEN+COS(TAofVE)) MAofVE = EAofVE - ECCEN*SIN(EAofVE) C PERIHE = VE200D - MAofVE*EDAYzY/TWOPI MA = MODULO (TWOPI*(DAY-VE200D)/EDAYzY + MAofVE, TWOPI) C**** C**** Numerically invert Kepler's equation: MA = EA - e*sin(EA) C**** EA = MA + ECCEN*(SIN(MA) + ECCEN*SIN(2.*MA)/2.) 10 dEA = (MA - EA + ECCEN*SIN(EA)) / (1. - ECCEN*COS(EA)) EA = EA + dEA IF(ABS(dEA).gt.1.d-10) GO TO 10 C**** C**** Calculate distance to Sun and true anomaly C**** DIST = 1. - ECCEN*COS(EA) TA = ATAN2 (BSEMI*SIN(EA), COS(EA)-ECCEN) C**** C**** Change reference frame to be nonrotating reference frame, angles C**** fixed according to stars, with Earth at center and positive x C**** axis be ray from Earth to Sun were Earth at vernal equinox, and C**** x-y plane be Earth's equatorial plane. Distance from current Sun C**** to this x axis is DIST sin(TA-TAofVE). At vernal equinox, Sun is C**** located at (DIST,0,0). At other times, Sun is located at: C**** C**** SUN = (DIST cos(TA-TAofVE), C**** DIST sin(TA-TAofVE) cos(OBLIQ), C**** DIST sin(TA-TAofVE) sin(OBLIQ)) C**** SIND = SIN(TA-TAofVE) * SIN(OBLIQ) COSD = SQRT (1. - SIND*SIND) SUNX = COS(TA-TAofVE) SUNY = SIN(TA-TAofVE) * COS(OBLIQ) SLNORO = ATAN2 (SUNY,SUNX) C**** C**** Determine Sun location in Earth's rotating reference frame C**** (normal longitude and latitude) C**** VEQLON = TWOPI*VE200D - TWOPI/2. + MAofVE - TAofVE ! modulo 2*Pi ROTATE = TWOPI*(DAY-VE200D)*(EDAYzY+1.)/EDAYzY SUNLON = MODULO (SLNORO-ROTATE-VEQLON, TWOPI) IF(SUNLON.gt.TWOPI/2.) SUNLON = SUNLON - TWOPI SUNLAT = ASIN (SIN(TA-TAofVE)*SIN(OBLIQ)) C**** C**** Determine longitude of fictitious mean Sun C**** Calculate Equation of Time C**** SLMEAN = TWOPI/2. - TWOPI*(DAY-FLOOR(DAY)) EQTIME = MODULO (SLMEAN-SUNLON, TWOPI) IF(EQTIME.gt.TWOPI/2.) EQTIME = EQTIME - TWOPI C**** RETURN END