darwin3 |
|
|
|
|||
File indexing completed on 2024-12-17 18:37:27 UTC
view on githubraw file Latest commit 87dd4f7d on 2024-01-17 18:17:24 UTC87dd4f7d5f Oliv*0001 #include "OASIM_OPTIONS.h" 0002 0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| 0004 CBOP 0005 C !ROUTINE: OASIM_SUN_VECTOR 0006 0007 C !INTERFACE: 0008 SUBROUTINE OASIM_SUN_VECTOR( 0009 I iday, iyr, gmt, 0010 O sunvec, rs ) 0011 0012 C !DESCRIPTION: 0013 C compute sun vector in earth-fixed coordinate system 0014 C 0015 C Together with a local zenith-pointing vector, this can be used 0016 C to compute the solar zenith angle. 0017 C 0018 C Original comment: 0019 C Given year, day of year, time in hours (GMT) and latitude and 0020 C longitude, returns an accurate solar zenith and azimuth angle. 0021 C Based on IAU 1976 Earth ellipsoid. Method for computing solar 0022 C vector and local vertical from Patt and Gregg, 1994, Int. J. 0023 C Remote Sensing. Only outputs solar zenith angle. This version 0024 C utilizes a pre-calculation of the local up, north, and east 0025 C vectors, since the locations where the solar zenith angle are 0026 C calculated in the model are fixed. 0027 C 0028 C !USES: 0029 IMPLICIT NONE 0030 #include "SIZE.h" 0031 #include "OASIM_SIZE.h" 0032 #include "OASIM_INTERNAL.h" 0033 0034 C !INPUT PARAMETERS: 0035 C iday :: yearday (1..366) 0036 C iyr :: year (e.g., 1970) 0037 C gmt :: time of day in hours (e.g., 12.5) 0038 INTEGER iday, iyr 0039 _RL gmt 0040 0041 C !OUTPUT PARAMETERS: 0042 C sunvec(3) :: sun vector in earth-fixed coordinate system 0043 C rs :: correction to earth-sun distance 0044 _RL sunvec(3), rs 0045 CEOP 0046 0047 #ifdef ALLOW_OASIM 0048 0049 C !FUNCTIONS: 0050 EXTERNAL OASIM_JULIAN_DAY 0051 INTEGER OASIM_JULIAN_DAY 0052 0053 C !LOCAL VARIABLES: 0054 INTEGER imon, jday, nutime, nt, iiday 0055 _RL rjd, t, sec, day, suni(3), fday, gha, ghar 0056 _RL xls, gs, xlm, omega, dpsi, eps 0057 #ifdef OASIM_DAILY_NUTATE 0058 COMMON/oasim_eph_parms/ dpsi, eps, nutime 0059 DATA nutime/-99999/ 0060 #endif 0061 DATA imon/1/ 0062 0063 sec = gmt*3600.0 _d 0 0064 c Compute floating point days since Jan 1.5, 2000 0065 c Note that the Julian day starts at noon on the specified date 0066 rjd = FLOAT(OASIM_JULIAN_DAY(iyr,imon,iday)) 0067 t = rjd - 2451545.0 _d 0 + (sec-43200.0 _d 0)/86400.0 _d 0 0068 CALL OASIM_EPHPARMS(t,xls,gs,xlm,omega) 0069 #ifdef OASIM_DAILY_NUTATE 0070 nt = INT(t) 0071 IF (nt.NE.nutime) THEN 0072 nutime = nt 0073 CALL OASIM_NUTATE(rad,t,xls,gs,xlm,omega,dpsi,eps) 0074 ENDIF 0075 #else 0076 CALL OASIM_NUTATE(rad,t,xls,gs,xlm,omega,dpsi,eps) 0077 #endif 0078 0079 c Compute sun vector 0080 c Compute unit sun vector in geocentric inertial coordinates 0081 CALL OASIM_SUN2000(rad,t,xls,gs,xlm,omega,dpsi,eps,suni,rs) 0082 0083 c Get Greenwich mean sidereal angle 0084 #ifdef OASIM_DAILY_NUTATE 0085 c Compute days since J2000 0086 day = FLOAT(iday) + sec/86400.0 _d 0 0087 iiday = INT(day) 0088 fday = day - iiday 0089 #else 0090 fday = sec/86400.0 _d 0 0091 #endif 0092 CALL OASIM_GHA2000(rad, t, fday, dpsi, eps, gha) 0093 ghar = gha/rad 0094 0095 c Transform Sun vector into geocentric rotating frame 0096 sunvec(1) = suni(1)*COS(ghar) + suni(2)*SIN(ghar) 0097 sunvec(2) = suni(2)*COS(ghar) - suni(1)*SIN(ghar) 0098 sunvec(3) = suni(3) 0099 0100 #endif /* ALLOW_OASIM */ 0101 0102 RETURN 0103 END 0104 0105 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| 0106 CBOP 0 0107 C !ROUTINE: OASIM_SUN2000 0108 0109 C !INTERFACE: 0110 SUBROUTINE OASIM_SUN2000( 0111 I radeg, t, xls , gs, xlm, omega, dpsi, eps, 0112 O sunvec, rs ) 0113 0114 C !DESCRIPTION: 0115 C This subroutine computes the Sun vector in geocentric inertial 0116 C (equatorial) coodinates. It uses the model referenced in The 0117 C Astronomical Almanac for 1984, Section S (Supplement) and documented 0118 C in Exact closed-form geolocation algorithm for Earth survey 0119 C sensors, by F.S. Patt and W.W. Gregg, Int. Journal of Remote 0120 C Sensing, 1993. The accuracy of the Sun vector is approximately 0.1 0121 C arcminute. 0122 C 0123 C Coded by: Frederick S. Patt, GSC, November 2, 1992 0124 C Modified to include Earth constants subroutine by W. Gregg, 0125 C May 11, 1993. 0126 0127 C !USES: 0128 IMPLICIT NONE 0129 0130 C !INPUT PARAMETERS: 0131 C radeg :: conversion factor for radians to degrees 0132 C t :: floating point days since Jan 1.5, 2000 0133 C xls :: Mean solar longitude (degrees) 0134 C gs :: Mean solar anomaly (degrees) 0135 C xlm :: Mean lunar longitude (degrees) 0136 C omega :: Ascending node of mean lunar orbit (degrees) 0137 C dpsi :: Nutation in longitude (degrees) 0138 C eps :: Obliquity of the Ecliptic (degrees) 0139 _RL radeg, t 0140 _RL xls, gs, xlm, omega, dpsi, eps 0141 0142 C !OUTPUT PARAMETERS: 0143 C sunvec :: Unit sun vector in geocentric inertial coords of date 0144 C rs :: Magnitude of the Sun vector (AU) 0145 _RL sunvec(3), rs 0146 CEOP 0147 0148 #ifdef ALLOW_OASIM 0149 0150 C !LOCAL VARIABLES: 0151 INTEGER nt 0152 _RL g2, g4, g5, dls, xlsg, xlsa 0153 0154 _RL xk 0155 DATA xk/0.0056932 _d 0/ !Constant of aberration 0156 0157 c Compute planet mean anomalies 0158 c Venus Mean Anomaly 0159 g2 = 50.40828 _d 0 + 1.60213022 _d 0*t 0160 g2 = MOD(g2,360.0 _d 0) 0161 0162 c Mars Mean Anomaly 0163 g4 = 19.38816 _d 0 + 0.52402078 _d 0*t 0164 g4 = MOD(g4,360.0 _d 0) 0165 0166 c Jupiter Mean Anomaly 0167 g5 = 20.35116 _d 0 + 0.08309121 _d 0*t 0168 g5 = MOD(g5,360.0 _d 0) 0169 0170 c Compute solar distance (AU) 0171 rs = 1.00014 _d 0 - 0.01671 _d 0*COS(gs/radeg) 0172 & - 0.00014 _d 0*COS(2.0 _d 0*gs/radeg) 0173 0174 c Compute Geometric Solar Longitude 0175 dls = (6893.0 _d 0 - 4.6543463 _d -4*t)*SIN(gs/radeg) 0176 & + 72.0 _d 0*SIN(2.0 _d 0*gs/radeg) 0177 & - 7.0 _d 0*COS((gs - g5)/radeg) 0178 & + 6.0 _d 0*SIN((xlm - xls)/radeg) 0179 & + 5.0 _d 0*SIN((4 _d 0*gs - 8 _d 0*g4 + 3 _d 0*g5)/radeg) 0180 & - 5.0 _d 0*COS((2 _d 0*gs - 2 _d 0*g2)/radeg) 0181 & - 4.0 _d 0*SIN((gs - g2)/radeg) 0182 & + 4.0 _d 0*COS((4 _d 0*gs - 8 _d 0*g4 + 3 _d 0*g5)/radeg) 0183 & + 3.0 _d 0*SIN((2 _d 0*gs - 2 _d 0*g2)/radeg) 0184 & - 3.0 _d 0*SIN(g5/radeg) 0185 & - 3.0 _d 0*SIN((2.0 _d 0*gs - 2.0 _d 0*g5)/radeg) !arcseconds 0186 0187 xlsg = xls + dls/3600.0 _d 0 0188 0189 c Compute Apparent Solar Longitude; includes corrections for nutation 0190 c in longitude and velocity aberration 0191 xlsa = xlsg + dpsi - xk/rs 0192 0193 c Compute unit Sun vector 0194 sunvec(1) = COS(xlsa/radeg) 0195 sunvec(2) = SIN(xlsa/radeg)*COS(eps/radeg) 0196 sunvec(3) = sin(xlsa/radeg)*SIN(eps/radeg) 0197 c type *,' Sunlon = ',xlsg,xlsa,eps 0198 0199 #endif /* ALLOW_OASIM */ 0200 0201 RETURN 0202 END 0203 0204 0205 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| 0206 CBOP 0 0207 C !ROUTINE: OASIM_GHA2000 0208 0209 C !INTERFACE: 0210 SUBROUTINE OASIM_GHA2000( radeg,t,fday,dpsi,eps,gha ) 0211 0212 C !DESCRIPTION: 0213 C This subroutine computes the Greenwich hour angle in degrees for the 0214 C input time. It uses the model referenced in The Astronomical Almanac 0215 C for 1984, Section S (Supplement) and documented in Exact 0216 C closed-form geolocation algorithm for Earth survey sensors, by 0217 C F.S. Patt and W.W. Gregg, Int. Journal of Remote Sensing, 1993. 0218 C It includes the correction to mean sideral time for nutation 0219 C as well as precession. 0220 C 0221 C Program written by: Frederick S. Patt 0222 C General Sciences Corporation 0223 C November 2, 1992 0224 0225 C !USES: 0226 IMPLICIT NONE 0227 0228 C !INPUT PARAMETERS: 0229 C t :: floating point days since Jan 1.5, 2000 0230 C fday :: fractional day 0231 C dpsi :: Nutation in longitude (degrees) 0232 C eps :: Obliquity of the Ecliptic (degrees) 0233 _RL radeg, t, fday, dpsi, eps 0234 0235 C !OUTPUT PARAMETERS: 0236 C gha :: Greenwich hour angle (degrees) 0237 _RL gha 0238 CEOP 0239 0240 #ifdef ALLOW_OASIM 0241 0242 C !LOCAL VARIABLES: 0243 _RL gmst 0244 0245 c Compute Greenwich Mean Sidereal Time (degrees) 0246 gmst = 100.4606184 _d 0 + 0.9856473663 _d 0*t + 2.908 _d -13*t*t 0247 0248 c Include apparent time correction and time-of-day 0249 gha = gmst + dpsi*COS(eps/radeg) + fday*360.0 _d 0 0250 gha = MOD(gha,360.0 _d 0) 0251 IF (gha.LT.0.0 _d 0) THEN 0252 gha = gha + 360.0 _d 0 0253 ENDIF 0254 0255 #endif /* ALLOW_OASIM */ 0256 0257 RETURN 0258 END 0259 0260 0261 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| 0262 CBOP 0 0263 C !ROUTINE: OASIM_EPHPARMS 0264 0265 C !INTERFACE: 0266 SUBROUTINE OASIM_EPHPARMS( t, 0267 O xls, gs, xlm, omega ) 0268 0269 C !DESCRIPTION: 0270 C This subroutine computes ephemeris parameters used by other Mission 0271 C Operations routines: the solar mean longitude and mean anomaly, and 0272 C the lunar mean longitude and mean ascending node. It uses the model 0273 C referenced in The Astronomical Almanac for 1984, Section S 0274 C (Supplement) and documented and documented in Exact closed-form 0275 C geolocation algorithm for Earth survey sensors, by F.S. Patt and 0276 C W.W. Gregg, Int. Journal of Remote Sensing, 1993. These parameters 0277 C are used to compute the solar longitude and the nutation in 0278 C longitude and obliquity. 0279 C 0280 C Program written by: Frederick S. Patt 0281 C General Sciences Corporation 0282 C November 2, 1992 0283 0284 C !USES: 0285 IMPLICIT NONE 0286 0287 C !INPUT PARAMETERS: 0288 C t :: time in days since January 1, 2000 at 12 hours UT 0289 _RL t 0290 0291 C !OUTPUT PARAMETERS: 0292 C xls :: Mean solar longitude (degrees) 0293 C gs :: Mean solar anomaly (degrees) 0294 C xlm :: Mean lunar longitude (degrees) 0295 C omega :: Ascending node of mean lunar orbit (degrees) 0296 _RL xls, gs, xlm, omega 0297 CEOP 0298 0299 #ifdef ALLOW_OASIM 0300 0301 c Sun Mean Longitude 0302 xls = 280.46592 _d 0 + 0.9856473516 _d 0*t 0303 xls = MOD(xls, 360.0 _d 0) 0304 0305 c Sun Mean Anomaly 0306 gs = 357.52772 _d 0 + 0.9856002831 _d 0*t 0307 gs = MOD(gs, 360.0 _d 0) 0308 0309 c Moon Mean Longitude 0310 xlm = 218.31643 _d 0 + 13.17639648 _d 0*t 0311 xlm = MOD(xlm, 360.0 _d 0) 0312 0313 c Ascending Node of Moons Mean Orbit 0314 omega = 125.04452 _d 0 - 0.0529537648 _d 0*t 0315 omega = MOD(omega, 360.0 _d 0) 0316 0317 #endif /* ALLOW_OASIM */ 0318 0319 RETURN 0320 END 0321 0322 0323 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| 0324 CBOP 0 0325 C !ROUTINE: OASIM_NUTATE 0326 0327 C !INTERFACE: 0328 SUBROUTINE OASIM_NUTATE( radeg,t,xls,gs,xlm,omega, 0329 O dpsi,eps ) 0330 0331 C !DESCRIPTION: 0332 C This subroutine computes the nutation in longitude and the obliquity 0333 C of the ecliptic corrected for nutation. It uses the model referenced 0334 C in The Astronomical Almanac for 1984, Section S (Supplement) and 0335 C documented in Exact closed-form geolocation algorithm for Earth 0336 C survey sensors, by F.S. Patt and W.W. Gregg, Int. Journal of 0337 C Remote Sensing, 1993. These parameters are used to compute the 0338 C apparent time correction to the Greenwich Hour Angle and for the 0339 C calculation of the geocentric Sun vector. The input ephemeris 0340 C parameters are computed using subroutine ephparms. Terms are 0341 C included to 0.1 arcsecond. 0342 C 0343 C Program written by: Frederick S. Patt 0344 C General Sciences Corporation 0345 C October 21, 1992 0346 0347 C !USES: 0348 IMPLICIT NONE 0349 0350 C !INPUT PARAMETERS: 0351 C radeg :: conversion factor for radians to degrees 0352 C t :: time in days since January 1, 2000 at 12 hours UT 0353 C xls :: Mean solar longitude (degrees) 0354 C gs :: Mean solar anomaly (degrees) 0355 C xlm :: Mean lunar longitude (degrees) 0356 C omega :: Ascending node of mean lunar orbit (degrees) 0357 _RL radeg, t, xls, gs, xlm, omega 0358 0359 C !OUTPUT PARAMETERS: 0360 C dpsi :: Nutation in longitude (degrees) 0361 C eps :: Obliquity of the Ecliptic (degrees) 0362 C (includes nutation in obliquity) 0363 _RL dpsi, eps 0364 CEOP 0365 0366 #ifdef ALLOW_OASIM 0367 0368 C !LOCAL VARIABLES: 0369 _RL epsm, deps 0370 0371 c Nutation in Longitude 0372 dpsi = - 17.1996 _d 0*SIN(omega/radeg) 0373 & + 0.2062 _d 0*SIN(2.0 _d 0*omega/radeg) 0374 & - 1.3187 _d 0*SIN(2.0 _d 0*xls/radeg) 0375 & + 0.1426 _d 0*SIN(gs/radeg) 0376 & - 0.2274 _d 0*SIN(2.0 _d 0*xlm/radeg) 0377 0378 c Mean Obliquity of the Ecliptic 0379 epsm = 23.439291 _d 0 - 3.560 _d -7*t 0380 0381 c Nutation in Obliquity 0382 deps = 9.2025 _d 0*COS(omega/radeg) 0383 & + 0.5736 _d 0*COS(2.0 _d 0*xls/radeg) 0384 0385 c True Obliquity of the Ecliptic 0386 eps = epsm + deps/3600.0 _d 0 0387 0388 dpsi = dpsi/3600.0 _d 0 0389 0390 #endif /* ALLOW_OASIM */ 0391 0392 RETURN 0393 END 0394 0395 0396 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| 0397 CBOP 0398 C !ROUTINE: OASIM_JULIAN_DAY 0399 0400 C !INTERFACE: 0401 FUNCTION OASIM_JULIAN_DAY( i, j, k ) 0402 0403 C !DESCRIPTION: 0404 C This function converts a calendar date to the corresponding Julian 0405 C day starting at noon on the calendar date. The algorithm used is 0406 C from Van Flandern and Pulkkinen, Ap. J. Supplement Series 41, 0407 C November 1979, p. 400. 0408 C 0409 C Written by Frederick S. Patt, GSC, November 4, 1992 0410 0411 C !USES: 0412 IMPLICIT NONE 0413 0414 C !INPUT PARAMETERS: 0415 C i :: Year - e.g. 1970 0416 C j :: Month - (1-12) 0417 C k :: Day - (1-31) 0418 INTEGER i, j, k 0419 0420 C !OUTPUT PARAMETERS: 0421 C oasim_julian_day :: Julian day 0422 INTEGER OASIM_JULIAN_DAY 0423 CEOP 0424 0425 #ifdef ALLOW_OASIM 0426 0427 OASIM_JULIAN_DAY= 367*i - 7*(i+(j+9)/12)/4 + 275*j/9 + k + 1721014 0428 0429 c This additional calculation is needed only for dates outside of the 0430 c period March 1, 1900 to February 28, 2100 0431 c oasim_julian_day = oasim_julian_day + 15 - 3*((i+(j-9)/7)/100+1)/4 0432 0433 #endif /* ALLOW_OASIM */ 0434 0435 RETURN 0436 END
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated from https://github.com/darwinproject/darwin3 by the 2.3.7-MITgcm-0.1 LXR engine. The LXR team |
|