!------------------------------------------------------------------------------ ! Global Modeling and Assimilation Office (GMAO) ! ! Goddard Earth Observing System (GEOS) ! ! MAPL Component ! !------------------------------------------------------------------------------ ! #include "MAPL_ErrLog.h" ! !> !### MODULE: `MAPL_SunMod` ! ! Author: GMAO SI-Team ! ! This class is intended to manage the sun`s position and provide ! the insolation at the top of the atmosphere. The main method ! is `GEOS_SunGetInsolation`, which depends on an Orbit object. ! The Orbit object defines this class and has public opaque type `GEOS_SunOrbit`. ! Methods are provided for creating it, destroying it, and making various queries. ! module MAPL_SunMod ! !USES: use ESMF use MAPL_Constants use MAPL_BaseMod use MAPL_IOMod use MAPL_CommsMod use MAPL_ExceptionHandling use netcdf use, intrinsic :: iso_fortran_env, only: REAL64 use pflogger, only: logging, Logger implicit none private ! !PUBLIC MEMBER FUNCTIONS: public MAPL_SunOrbitCreate public MAPL_SunOrbitCreateFromConfig public MAPL_SunOrbitCreated public MAPL_SunOrbitDestroy public MAPL_SunOrbitQuery public MAPL_SunGetInsolation public MAPL_SunGetSolarConstant public MAPL_SunGetDaylightDuration public MAPL_SunGetDaylightDurationMax public MAPL_SunGetLocalSolarHourAngle ! !PUBLIC TYPES: public MAPL_SunOrbit !EOP integer, public, parameter :: MAPL_SunAutumnalEquinox = 1 integer, public, parameter :: MAPL_SunWinterSolstice = 2 integer, public, parameter :: MAPL_SunVernalEquinox = 3 integer, public, parameter :: MAPL_SunSummerSolstice = 4 integer, public, parameter :: MAPL_SunDailyMean = 5 integer, public, parameter :: MAPL_SunAnnualMean = 6 ! Default solar orbital system parameters (private). ! Dont change these unless you know what you are doing. ! They are appropriate for the current modern epoch circa 2000. ! ------------------------------------------------------------- ! Parameters of old orbital system (tabularized intercalation cycle) ! ------------------------------------------------------------------ real, parameter :: DEFAULT_ORBIT_ECCENTRICITY = 0.0167 real, parameter :: DEFAULT_ORBIT_OBLIQUITY = 23.45 ! degrees real, parameter :: DEFAULT_ORBIT_PERIHELION = 102.0 ! degrees integer, parameter :: DEFAULT_ORBIT_EQUINOX = 80 ! days ! Parameters of new orbital system (analytic two-body), which allows some ! time-varying behavior, namely, linear variation in LAMBDAP, ECC, and OBQ. ! ------------------------------------------------------------------------- ! Fixed anomalistic year length in mean solar days real, parameter :: DEFAULT_ORB2B_YEARLEN = 365.2596 ! Reference date and time for orbital parameters ! (defaults to J2000 = 01Jan2000 12:00:00 TT = 11:58:56 UTC) integer, parameter :: DEFAULT_ORB2B_REF_YYYYMMDD = 20000101 integer, parameter :: DEFAULT_ORB2B_REF_HHMMSS = 115856 ! Orbital eccentricity at reference date real, parameter :: DEFAULT_ORB2B_ECC_REF = 0.016710 ! Rate of change of orbital eccentricity per Julian century real, parameter :: DEFAULT_ORB2B_ECC_RATE = -4.2e-5 ! Earth's obliquity (axial tilt) at reference date [degrees] real, parameter :: DEFAULT_ORB2B_OBQ_REF = 23.44 ! Rate of change of obliquity [degrees per Julian century] real, parameter :: DEFAULT_ORB2B_OBQ_RATE = -1.3e-2 ! Longitude of perihelion at reference date [degrees] ! (from March equinox to perihelion in direction of earth's motion) real, parameter :: DEFAULT_ORB2B_LAMBDAP_REF = 282.947 ! Rate of change of LAMBDAP [degrees per Julian century] ! (Combines both equatorial and ecliptic precession) real, parameter :: DEFAULT_ORB2B_LAMBDAP_RATE = 1.7195 ! March Equinox date and time ! (defaults to March 20, 2000 at 07:35:00 UTC) integer, parameter :: DEFAULT_ORB2B_EQUINOX_YYYYMMDD = 20000320 integer, parameter :: DEFAULT_ORB2B_EQUINOX_HHMMSS = 73500 ! ------------------------------------------------------------- interface MAPL_SunGetInsolation module procedure SOLAR_1D module procedure SOLAR_2D module procedure SOLAR_ARR_INT end interface interface MAPL_SunGetSolarConstant module procedure MAPL_SunGetSolarConstantByYearDoY module procedure MAPL_SunGetSolarConstantByTime module procedure MAPL_SunGetSolarConstantFromNetcdfFile module procedure MAPL_SunGetSolarConstantFromNRLFile end interface type MAPL_SunOrbit private logical :: CREATED = .FALSE. type(ESMF_Clock) :: CLOCK real :: OB, ECC, PER, YEARLEN integer :: EQNX, YEARS_PER_CYCLE, DAYS_PER_CYCLE real, pointer, dimension(:) :: ZC => null() real, pointer, dimension(:) :: ZS => null() real, pointer, dimension(:) :: PP => null() real, pointer, dimension(:) :: TH => null() real, pointer, dimension(:) :: ET => null() logical :: EOT logical :: FIX_SUN logical :: ANAL2B real :: ORB2B_YEARLEN type(ESMF_Time) :: ORB2B_TIME_REF real :: ORB2B_ECC_REF real :: ORB2B_ECC_RATE real :: ORB2B_OBQ_REF real :: ORB2B_OBQ_RATE real :: ORB2B_LAMBDAP_REF real :: ORB2B_LAMBDAP_RATE type(ESMF_Time) :: ORB2B_TIME_EQUINOX real :: ORB2B_OMG0 type(ESMF_Time) :: ORB2B_TIME_PERI end type MAPL_SunOrbit contains !========================================================================== ! utlity functions ! rectify to [-pi,+pi) function RECT_PMPI(X) real :: X, RECT_PMPI RECT_PMPI = MODULO( X + MAPL_PI, 2*MAPL_PI ) - MAPL_PI end function ! true anomaly from eccentric anomaly function calcTAfromEA(EA,EAFAC) result(TA) real :: EA, EAFAC, TA TA = 2. * atan( tan(EA / 2.) / EAFAC) end function ! eccentric anomaly from true anomaly function calcEAfromTA(TA,EAFAC) result(EA) real :: TA, EAFAC, EA EA = 2. * atan( EAFAC * tan(TA / 2.)) end function ! mean anomaly from eccentric anomaly (Kepler's equation) function calcMAfromEA(EA,ECC) result(MA) real :: EA, ECC, MA MA = EA - ECC * sin(EA) end function ! eccentric anomaly from mean anomaly ! (invert Kepler's equation by Newton-Raphson) subroutine invert_Keplers_Newton( & MA, ECC, & EA, dE, nits, & tol, max_its) real, intent(in ) :: MA real, intent(in ) :: ECC real, intent(out) :: EA real, intent(out) :: dE integer, intent(out) :: nits real, optional, intent(in ) :: tol integer, optional, intent(in ) :: max_its real :: f, df, tol_ integer :: max_its_ if (present(tol)) then tol_ = tol else tol_ = 1.e-10 endif if (present(max_its)) then max_its_ = max_its else max_its_ = 10 endif EA = MA do nits = 1, max_its_ f = EA - ECC * sin(EA) - MA df = 1. - ECC * cos(EA) dE = -f / df EA = EA + dE if (abs(dE) < tol_) exit enddo end subroutine ! Earth-Sun distance from true anomaly function calcRadfromTA(TA,ECC,OMSQECC) result(Rad) real :: TA, ECC, OMSQECC, Rad Rad = OMSQECC / (1. + ECC * cos(TA)) end function !========================================================================== !> ! Integrates the earth`s orbit and stores the necessary ! parameters to easily compute the earth`s position for each day ! of the full (usually 4-year) intercalation cycle. ! The orbital parameters are passed as arguments. ! The full calendar intercalation cycle is obtained from the ! ESMF clock passed as an argument. This becomes the orbit`s ! attached clock. Currently we assume a single intercalation. ! ! A good introduction to celestial mechanics for understanding ! this code can be found in Blanco & McCuskey, 1961: "Basic ! Physics of the Solar System", hereafter BM. ! !``` ! CLOCK: The orbit will depend on the calendar in this clock ! This is used for the length of year, to set intercalation cycle. ! ! ECCENTRICITY: Eccentricity of the Earth`s orbit. ! ! PERIHELION: Longitude of perihelion, measured in degrees from ! autumnal equinox in the direction of the Earth`s motion. ! ! OBLIQUITY: Tilt of the Earth`s rotation axis from a ! normal to the plane of the orbit. In degrees. ! ! EQUINOX: Day of year of vernal equinox. ! Equinox is assumed to occur at 0Z on this day on the ! first year of the cycle. ! ! EOT: Apply Equation of Time correction? ! ! ORBIT_ANAL2B: New orbital system (analytic two-body) allows some ! time-varying behavior, namely, linear time variation in LAMBDAP, ! ECC, and OBQ. If .TRUE., the following ORB2B parameters are used ! and only CLOCK and EOT above are used, i.e., the ECCENTRICITY, ! OBLIQUITY, PERIHELION and EQUINOX above are NOT used and are ! replaced by the relevant ORB2B parameters below. ! ! ORB2B_YEARLEN: Fixed anomalistic year length in mean solar days. ! ! ORB2B_REF_YYYYMMDD: Reference date for orbital parameters. ! ! ORB2B_REF_HHMMSS: Reference time for orbital parameters. ! ! ORB2B_ECC_REF: Orbital eccentricity at reference date. ! ! ORB2B_ECC_RATE: Rate of change of orbital eccentricity per Julian century. ! ! ORB2B_OBQ_REF: Earth's obliquity (axial tilt) at reference date [degrees]. ! ! ORB2B_OBQ_RATE: Rate of change of obliquity [degrees per Julian century]. ! ! ORB2B_LAMBDAP_REF: Longitude of perihelion at reference date [degrees] ! from March equinox to perihelion in direction of earth's motion). ! ! ORB2B_LAMBDAP_RATE: Rate of change of LAMBDAP [degrees per Julian century] ! (Combines both equatorial and ecliptic precession). ! ! ORB2B_EQUINOX_YYYYMMDD: March equinox date. ! ! ORB2B_EQUINOX_HHMMSS: March equinox time. !``` ! type(MAPL_SunOrbit) function MAPL_SunOrbitCreate(CLOCK, & ECCENTRICITY, & OBLIQUITY, & PERIHELION, & EQUINOX, & EOT, & ORBIT_ANAL2B, & ORB2B_YEARLEN, & ORB2B_REF_YYYYMMDD, & ORB2B_REF_HHMMSS, & ORB2B_ECC_REF, & ORB2B_ECC_RATE, & ORB2B_OBQ_REF, & ORB2B_OBQ_RATE, & ORB2B_LAMBDAP_REF, & ORB2B_LAMBDAP_RATE, & ORB2B_EQUINOX_YYYYMMDD, & ORB2B_EQUINOX_HHMMSS, & FIX_SUN, & RC ) ! !ARGUMENTS: type(ESMF_Clock) , intent(IN ) :: CLOCK real , intent(IN ) :: ECCENTRICITY real , intent(IN ) :: OBLIQUITY real , intent(IN ) :: PERIHELION integer , intent(IN ) :: EQUINOX logical , intent(IN ) :: EOT logical , intent(IN ) :: ORBIT_ANAL2B real , intent(IN ) :: ORB2B_YEARLEN integer , intent(IN ) :: ORB2B_REF_YYYYMMDD integer , intent(IN ) :: ORB2B_REF_HHMMSS real , intent(IN ) :: ORB2B_ECC_REF real , intent(IN ) :: ORB2B_ECC_RATE real , intent(IN ) :: ORB2B_OBQ_REF real , intent(IN ) :: ORB2B_OBQ_RATE real , intent(IN ) :: ORB2B_LAMBDAP_REF real , intent(IN ) :: ORB2B_LAMBDAP_RATE integer , intent(IN ) :: ORB2B_EQUINOX_YYYYMMDD integer , intent(IN ) :: ORB2B_EQUINOX_HHMMSS logical, optional , intent(IN ) :: FIX_SUN integer, optional , intent(OUT) :: RC !EOPI ! ======================================= ! PMN Dec 2019: Notes on Equation of Time ! ======================================= ! (Part of a more complete analysis available in orbit.pdf) ! ! @ Introduction: ! ! The Earth rotates on its axis with a period T_S called the sidereal ! day (after the Latin for "star", since it is the rotation period of ! the Earth with respect to distant stars). T_S is slightly shorter ! than the so-called mean solar day, or clock day, of duration T_M = ! 86400 seconds. This is because the Earth is a prograde planet, i.e., ! it rotates about its axis in the same sense (counterclockwise look- ! ing down on the North Pole) as it orbits the sun. Specifically, say ! the sun crosses the meridian of some location at a particular time. ! And imagine there is a distant star directly behind the sun at that ! moment. After one sidereal day the location will rotate 360 degrees ! about the earth's axis and the distant star will again cross its ! merdian. But during that period the earth will have moved a small ! counterclockwise distance around its orbit and so it will take a ! small additional rotation of the earth for the sun to also cross ! the meridian and thereby complete a solar day. ! ! Put another way, a solar day is slightly longer than a sidereal day ! because the sun appears to move slowly eastward across the celestial ! sphere with respect to distant stars as the year passes. The path ! of this motion is called the ecliptic. And clearly what governs the ! length of a solar day is the apparent velocity of the sun along the ! ecliptic, or, more particularly, the equatorial component of that ! velocity. But both the magnitude and equatorial component of the solar ! ecliptic velocity change during the year, the former because the earth's ! orbit is elliptical, not circular, and the latter because the earth's ! axis of rotation is tilted with respect to the orbital (ecliptic) plane. ! Thus the length of a solar day changes during the year. While these ! factors cause only a small perturbation to the length of the solar ! day (less than 30 seconds), the perturbations accumulate so that at ! different times of the year apparent solar time ("sundial time") and ! mean solar time ("clock time") can differ by as much as ~15 minutes. ! This difference is called the Equation of Time. ! ! To be more rigorous, consider a fictitious "Mean Sun" that moves at ! constant eastward speed around the celestial equator, completing a ! full revolution in a year, namely in the period Y * T_M, where Y is ! the number of mean solar days in a year (e.g., 365.25). Thus, in one ! mean solar day, T_M, the mean sun has moved an angle 2*PI/Y eastward. ! Hence, beyond one full earth revolution, period T_S, an additional ! earth rotation of (T_M-T_S) * 2*PI/T_S = 2*PI/Y is required to "catch ! up with the moving sun", as described earlier. Hence T_M - T_S = T_S / Y ! and so ! ! T_M = T_S (Y+1)/Y,..... (1) ! ! a constant (near unity) multiple of the fixed sidereal day. T_M is ! the length of the solar day for the "mean sun", or the "mean solar ! day". Because it is invariant during the year, it is convenient for ! timekeeping, and forms the basis for "mean solar time", which at ! Greenwich is essentially UTC. By *definition*, T_M = 24h = 86400s. ! That is, what we know as "hours", "minutes" and "seconds", are just ! convenient integer fractions of the mean solar day. In these units, ! the sidereal day T_S is approximately 23h 56m 4s. ! ! The solar zenith angle calculation (in MAPL_SunGetInsolation) needs ! the true local solar hour angle, h_T, which is the angle, measured ! westward along the equator, from the local meridian to the true sun. ! This is just the Greenwich solar hour angle, H_T, plus the longitude, ! so we will henceforth work exclusively with Greenwich hour angles. ! We should use the hour angle of the *true* sun, H_T, but a common ! approximation replaces this with the hour angle of the mean sun ! ! H_M = 2*PI*(u-0.5),..... (2) ! ! where u is UTC time (in days) and the offset is needed because the mean ! solar hour angle is zero at "noon". If more accuracy is required, the ! hour angle of the true sun is typically obtained as a small correction ! to H_M called the Equation of time, EOT: ! ! H_T = H_M + EOT, where EOT = H_T - H_M. ! ! As discussed above, EOT corrects for two factors: ! (a) the variable speed of the earth in its elliptical orbit about ! the sun (e.g., moving fastest at perihelion), and ! (b) the tilt of the earth's axis of rotation wrt its orbital plane ! (the "obliquity"), which causes the equatorial projection of ! the sun's apparent ecliptic motion to vary with the season ! (e.g., being parallel to the equator at the solstices.) ! ! @ Derivation of Equation of Time: ! ! We can write ! ! H_T = H_1 - (H_1 - H_T) = H_1 - a_T, ! ! where H_1 is the Greenwich hour angle of the First Point of Aries (the ! location of the vernal equinox, denoted "1PoA"), and is also known as the ! Greenwich Sidereal hour angle, and where a_T is the right ascension of the ! true sun (since the right ascension of any object is just the difference ! between the hour angles of 1PoA and the object). Hence, ! ! EOT = H_1 - H_M - a_T...... (3) ! ! All three terms on the right of (3) are time variable: a_T changes slowly ! throughout the year, and is known from the earth-sun two-body elliptical ! orbit solution, while H_1 and H_M vary rapidly with earth's rotation. H_M ! has a period of one mean solar day, T_M, and H_1 has a period of one ! sidereal day, T_S. ! ! It may seem from from (2) that the mean sun and its hour angle are fully ! specified. That, in fact, is not yet the case: (2) is really just a def- ! inition of UTC, namely, that one UTC day is one mean solar day and that ! the time of culmination of the mean sun, what we call "noon", occurs at ! UTC 12h00m. What we are still at liberty to do is specify the phasing of ! the mean sun in its equatorial orbit, e.g., by specifying the time u_R ! at which the mean sun passes through 1PoA (both on the equator). At this ! time, H_1(u_R) = H_M(u_R), and so ! ! H_1(u) - H_M(u) = 2*PI*(u-u_R)*(Y+1)/Y - 2*PI*(u-u_R) ! = 2*PI * (u-u_R) / Y ! = MA(u) - MA(u_R),... (4) ! ! where MA(u) = 2*PI * (u-u_P) / Y is the so-called "mean anomaly", known ! from the earth-sun two-body orbital solution, and u_P is the time of ! perihelion. Thus, to fully determine EOT, through (3) and (4), we need ! only to specify MA(u_R). ! ! To understand the mean anomaly MA, consider the standard two-body earth- ! sun problem in which the earth E moves in an elliptical orbit about the ! sun S at one focus, all in the ecliptic plane. The point on the ellipse ! closest to S is called the perihelion P. Obviously, the center of the ! ellipse O, the focus S and the perihelion P are co-linear, the so-called ! major axis of the ellipse. Additionally, let C be the circumscribing circle ! of the ellipse, with center O and passing through P (and the corresponding ! aphelion A). By Kepler's Second Law, the sun-earth vector sweeps out equal ! areas in equal times, so the fractional area of the elliptical sector PSE ! is a linear function of time, being zero at perihelion and one a year later. ! Specifically, this fractional area is none other than the scaled mean anomaly ! MA(u)/(2*PI) = (u - u_P) / Y. Clearly MA(u) can also be interpreted as an ! angle, the angle POQ of a point Q orbiting on the circumcircle C at constant ! speed in the same direction as the earth, also with a yearly period, and ! passing through P at the same time u_P as the earth. Thus the point Q can ! be conceptualized as a sort of "mean earth" orbiting a "second mean sun" ! (different from M above) at O. Note that while the angle MA(u) = angle POQ ! of this mean earth at time u is a linear function of time, the corresponding ! angle of the real earth, namely TA(u) = angle PSE, called the true anomaly, ! is a non-linear function of time, since the real earth has a variable speed ! in its elliptical orbit, e.g., moving faster at perihelion, so that its areal ! fraction is linear in time. The relationship between MA(u) and TA(u) is known ! from the orbital solution. Finally, the ecliptic longitude of the earth, ! lambda = angle 1SE is the angle at the sun, measured in the same direction ! as the earth's motion, from the First Point of Aries to the earth. Then ! ! TA(u) = angle PSE(u) = angle PS1 + angle 1SE(u) = lambda(u) - lambda_P, ! ! where lambda_P = lambda(u_P) = angle 1SP = -angle PS1 is known as the ! longitude of perihelion, and is currently about 283 deg, or equivalently ! -77 deg. ! ! With this background, we can understand the quantity MA(u_R) we are trying ! to specify. If we *choose* ! ! MA(u_R) = -lambda_P = angle PS1... (5) ! <==> angle POQ(u_R) = angle PS1, ! ! then at u_R, viewed from the mean earth Q, the second (ecliptic) mean sun ! O is in direction of 1PoA. And at that same time, by definition of u_R, ! the first (equatorial) mean sun M, as seen from the real earth E, is also ! in direction of 1PoA. ! ! I (PMN) have verified that the choice (5) for MA(u_R) leads to zero mean ! EOT to first order in the eccentricity, e. I cannot say, at this point, ! that it is generally true for any order in e. I add below a final explicit ! enforcement of zero mean EOT to the code. [I found that the order e^2 term ! for the mean EOT was NOT zero, and was larger than the mean EOT produced ! by this code (which is valid for all orders in e) before any explicity ! correction of the mean to zero. This suggests: (a) that I made a mistake ! in my calculations, or (b) that higer order e terms provide some cancel- ! ation.] ! ! Hence, finally, ! ! EOT(u) = MA(u) + PRHV - a_T(u)... (6) ! ! where PRHV is the name for lamba_P in the code. ! =========================================================================== ! Locals character(len=ESMF_MAXSTR), parameter :: IAm = "SunOrbitCreate" real(kind=REAL64) :: YEARLEN integer :: K, KP, YEARS_PER_CYCLE, DAYS_PER_CYCLE real(kind=REAL64) :: TREL, T1, T2, T3, T4 real(kind=REAL64) :: SOB, COB, OMG0, OMG, PRH, PRHV real :: OMECC, OPECC, OMSQECC, EAFAC real(kind=REAL64) :: TA, EA, MA, TRRA, MNRA real :: meanEOT type(MAPL_SunOrbit) :: ORBIT integer :: STATUS integer :: year, month, day, hour, minute, second real(ESMF_KIND_R8) :: days real :: ECC_EQNX, LAMBDAP_EQNX, EAFAC_EQNX real :: TA_EQNX, EA_EQNX, MA_EQNX type(ESMF_TimeInterval) :: DT ! record inputs needed by both orbit methods ORBIT%CLOCK = CLOCK ORBIT%EOT = EOT ORBIT%ANAL2B = ORBIT_ANAL2B ! Analytic two-body orbit? if (ORBIT_ANAL2B) then ! record inputs in ORBIT type associate(D2R => MAPL_DEGREES_TO_RADIANS) ORBIT%ORB2B_YEARLEN = ORB2B_YEARLEN ORBIT%ORB2B_ECC_REF = ORB2B_ECC_REF ORBIT%ORB2B_OBQ_REF = ORB2B_OBQ_REF * D2R ! radians ORBIT%ORB2B_LAMBDAP_REF = ORB2B_LAMBDAP_REF * D2R ! radians ORBIT%ORB2B_ECC_RATE = ORB2B_ECC_RATE / 36525. ! per day ORBIT%ORB2B_OBQ_RATE = ORB2B_OBQ_RATE * D2R / 36525. ! radians per day ORBIT%ORB2B_LAMBDAP_RATE = ORB2B_LAMBDAP_RATE * D2R / 36525. ! radians per day end associate ! record MAPL Time object for REFerence time year = ORB2B_REF_YYYYMMDD / 10000 month = mod(ORB2B_REF_YYYYMMDD, 10000) / 100 day = mod(ORB2B_REF_YYYYMMDD, 100) hour = ORB2B_REF_HHMMSS / 10000 minute = mod(ORB2B_REF_HHMMSS, 10000) / 100 second = mod(ORB2B_REF_HHMMSS, 100) call ESMF_TimeSet(ORBIT%ORB2B_TIME_REF, & yy=year, mm=month, dd=day, h=hour, m=minute, s=second, rc=STATUS) _VERIFY(STATUS) ! record MAPL Time object for EQUINOX year = ORB2B_EQUINOX_YYYYMMDD / 10000 month = mod(ORB2B_EQUINOX_YYYYMMDD, 10000) / 100 day = mod(ORB2B_EQUINOX_YYYYMMDD, 100) hour = ORB2B_EQUINOX_HHMMSS / 10000 minute = mod(ORB2B_EQUINOX_HHMMSS, 10000) / 100 second = mod(ORB2B_EQUINOX_HHMMSS, 100) call ESMF_TimeSet(ORBIT%ORB2B_TIME_EQUINOX, & yy=year, mm=month, dd=day, h=hour, m=minute, s=second, rc=STATUS) _VERIFY(STATUS) ! time-invariant precalculations ORBIT%ORB2B_OMG0 = 2. * MAPL_PI / ORB2B_YEARLEN ! at provided ORB2B_TIME_EQUINOX LAMBDA=0 by definition call ESMF_TimeIntervalGet( & ORBIT%ORB2B_TIME_EQUINOX - ORBIT%ORB2B_TIME_REF, & d_r8=days, rc=STATUS) _VERIFY(STATUS) ECC_EQNX = ORBIT%ORB2B_ECC_REF + days * ORBIT%ORB2B_ECC_RATE LAMBDAP_EQNX = ORBIT%ORB2B_LAMBDAP_REF + days * ORBIT%ORB2B_LAMBDAP_RATE EAFAC_EQNX = sqrt((1.-ECC_EQNX)/(1.+ECC_EQNX)) TA_EQNX = -LAMBDAP_EQNX ! since LAMBDA=0 EA_EQNX = calcEAfromTA(TA_EQNX,EAFAC_EQNX) MA_EQNX = calcMAfromEA(EA_EQNX,ECC_EQNX) ! Time of one perihelion (all subsequent ORB2B_YEARLEN apart) call ESMF_TimeIntervalSet(DT, d_r8 = dble(MA_EQNX / ORBIT%ORB2B_OMG0)) ORBIT%ORB2B_TIME_PERI = ORBIT%ORB2B_TIME_EQUINOX - DT else ! ================================== ! Standard tabularized intercalation ! ================================== ! MJS: This needs to come from the calendar when the time manager works right. YEARLEN = 365.25 ! Factors involving the orbital parameters !----------------------------------------- OMECC = 1. - ECCENTRICITY OPECC = 1. + ECCENTRICITY OMSQECC = 1. - ECCENTRICITY**2 ! pmn: consider changing to line below when zero-diff not issue ! OMSQECC = OMECC * OPECC EAFAC = sqrt(OMECC/OPECC) associate(D2R => MAPL_DEGREES_TO_RADIANS) OMG0 = 2.*MAPL_PI/YEARLEN OMG = OMG0/sqrt(OMSQECC)**3 PRH = PERIHELION*D2R SOB = sin(OBLIQUITY*D2R) COB = cos(OBLIQUITY*D2R) end associate ! PRH is the ecliptic longitude of the perihelion, measured (at the Sun) ! from the autumnal equinox in the direction of the Earth`s orbital motion ! (counterclockwise as viewed from ecliptic north pole). ! For EOT calculations we will reference the perihelion wrt to the vernal ! equinox, PRHV. Of course, the difference is just PI. ! pmn: once the EOT code is established and zero-diff not an issue, ! consider removing original PRH and changing the original (non-EOT), ! code, which employs ! cos(Y \pm PI) = -COS(Y) ! to use PRHV, namely ! -cos(Y-PRH) = cos(Y-PRH-PI) = cos(Y-PRHV) PRHV = PRH + MAPL_PI_R8 ! Compute length of leap cycle ! ---------------------------- if(YEARLEN-int(YEARLEN) > 0.) then YEARS_PER_CYCLE = nint(1./(YEARLEN-int(YEARLEN))) else YEARS_PER_CYCLE = 1 endif DAYS_PER_CYCLE=nint(YEARLEN*YEARS_PER_CYCLE) ! save inputs and intercalculation details ! ---------------------------------------- ORBIT%OB = OBLIQUITY ORBIT%ECC = ECCENTRICITY ORBIT%PER = PERIHELION ORBIT%EQNX = EQUINOX ORBIT%YEARLEN = YEARLEN ORBIT%YEARS_PER_CYCLE = YEARS_PER_CYCLE ORBIT%DAYS_PER_CYCLE = DAYS_PER_CYCLE ! Allocate orbital cycle outputs ! ------------------------------ ! TH: Ecliptic longitude of the true Sun, TREL [radians] ! ZS: Sine of declination ! ZC: Cosine of declination ! PP: Inverse of square of earth-sun distance [1/(AU**2)] ! ET: Equation of time [radians] = ! True solar hour angle - Mean solar hour angle if(associated(ORBIT%TH)) deallocate(ORBIT%TH) allocate(ORBIT%TH(DAYS_PER_CYCLE), stat=status) _VERIFY(STATUS) if(associated(ORBIT%ZC)) deallocate(ORBIT%ZC) allocate(ORBIT%ZC(DAYS_PER_CYCLE), stat=status) _VERIFY(STATUS) if(associated(ORBIT%ZS)) deallocate(ORBIT%ZS) allocate(ORBIT%ZS(DAYS_PER_CYCLE), stat=status) _VERIFY(STATUS) if(associated(ORBIT%PP)) deallocate(ORBIT%PP) allocate(ORBIT%PP(DAYS_PER_CYCLE), stat=status) _VERIFY(STATUS) if (ORBIT%EOT) then if(associated(ORBIT%ET)) deallocate(ORBIT%ET) allocate(ORBIT%ET(DAYS_PER_CYCLE), stat=status) _VERIFY(STATUS) end if ! Begin integration at the vernal equinox (K=1, KP=EQUINOX), at ! which, by defn, the ecliptic longitude of the true sun is zero. ! Right ascension of true sun at EQUINOX is also zero by defn. ! -------------------------------------------------------------- KP = EQUINOX TREL = 0. ORBIT%ZS(KP) = sin(TREL)*SOB ORBIT%ZC(KP) = sqrt(1.-ORBIT%ZS(KP)**2) ORBIT%PP(KP) = ( (1.-ECCENTRICITY*cos(TREL-PRH)) / OMSQECC )**2 ORBIT%TH(KP) = TREL if (ORBIT%EOT) then ! Right ascension of "mean sun" = MA(u) + PRHV [eqn (6) above]. ! Calcn of True (TA), Eccentric (EA), and Mean Anomaly (MA). TA = TREL - PRHV ! by defn of TA and PRHV EA = 2.d0*atan(EAFAC*tan(TA/2.d0)) ! BM 4-55 MA = EA - ECCENTRICITY*sin(EA) ! Kepler's eqn (BM 4-49 ff.) MNRA = MA + PRHV ! see EOT notes above TRRA = 0. ! because TREL=0 at Equinox ! Equation of Time, ET [radians] ! True Solar hour angle = Mean Solar hour angle + ET ! (hour angle and right ascension are in reverse direction) ORBIT%ET(KP) = rect_pmpi(real(MNRA - TRRA)) end if ! Integrate orbit for entire leap cycle using Runge-Kutta ! Mean sun moves at constant speed around Celestial Equator ! --------------------------------------------------------- do K=2,DAYS_PER_CYCLE T1 = dTRELdDAY(TREL ,OMG,ECCENTRICITY,PRH) T2 = dTRELdDAY(TREL+T1*0.5,OMG,ECCENTRICITY,PRH) T3 = dTRELdDAY(TREL+T2*0.5,OMG,ECCENTRICITY,PRH) T4 = dTRELdDAY(TREL+T3 ,OMG,ECCENTRICITY,PRH) KP = mod(KP,DAYS_PER_CYCLE) + 1 TREL = TREL + (T1 + 2.0*(T2 + T3) + T4) / 6.0 ORBIT%ZS(KP) = sin(TREL)*SOB ORBIT%ZC(KP) = sqrt(1.-ORBIT%ZS(KP)**2) ORBIT%PP(KP) = ( (1.-ECCENTRICITY*cos(TREL-PRH)) / OMSQECC )**2 ORBIT%TH(KP) = TREL if (ORBIT%EOT) then ! From BM 1-15 and 1-16 with beta=0 (Sun is on ecliptic), ! and dividing through by common cos(dec) since it does not ! affect the ratio of sin(RA) to cos(RA). TRRA = atan2(sin(TREL)*COB,cos(TREL)) MNRA = MNRA + OMG0 ORBIT%ET(KP) = rect_pmpi(real(MNRA - TRRA)) end if enddo ! enforce zero mean EOT (just in case) if (ORBIT%EOT) then meanEOT = sum(ORBIT%ET)/DAYS_PER_CYCLE ORBIT%ET = ORBIT%ET - meanEOT end if end if if (present(FIX_SUN)) then ORBIT%FIX_SUN=FIX_SUN else ORBIT%FIX_SUN=.FALSE. end if ORBIT%CREATED = .TRUE. MAPL_SunOrbitCreate = ORBIT _RETURN(ESMF_SUCCESS) contains real(kind=REAL64) function dTRELdDAY(TREL,OMG,ECCENTRICITY,PRH) real(kind=REAL64), intent(in) :: TREL ! ecliptic longitude of true Sun real(kind=REAL64), intent(in) :: OMG real, intent(in) :: ECCENTRICITY real(kind=REAL64), intent(in) :: PRH dTRELdDAY = OMG*(1.0-ECCENTRICITY*cos(TREL-PRH))**2 end function dTRELdDAY end function MAPL_SunOrbitCreate !========================================================================== !> ! The function `MAPL_SunOrbitCreateFromConfig` works like `MAPL_SunOrbitCreate()` ! but gets orbital parameters from Config CF. ! function MAPL_SunOrbitCreateFromConfig ( & CF, CLOCK, FIX_SUN, RC) result (ORBIT) ! !ARGUMENTS: type (ESMF_Config), intent(INOUT) :: CF type (ESMF_Clock), intent(IN ) :: CLOCK logical, intent(IN ) :: FIX_SUN integer, optional, intent(OUT ) :: RC type (MAPL_SunOrbit) :: ORBIT character(len=ESMF_MAXSTR), parameter :: IAm = "SunOrbitCreateFromConfig" integer :: STATUS real :: ECC, OB, PER integer :: EQNX logical :: EOT, ORBIT_ANAL2B integer :: ORB2B_REF_YYYYMMDD, ORB2B_REF_HHMMSS, & ORB2B_EQUINOX_YYYYMMDD, ORB2B_EQUINOX_HHMMSS real :: ORB2B_YEARLEN, & ORB2B_ECC_REF, ORB2B_ECC_RATE, & ORB2B_OBQ_REF, ORB2B_OBQ_RATE, & ORB2B_LAMBDAP_REF, ORB2B_LAMBDAP_RATE ! pmn: There is one orbit is per STATE, so, for example, the MAPL states of the ! solar and land gridded components can potentially have independent solar orbits. ! Usually these "independent orbits" will be IDENTICAL because the configuration ! resources such as "ECCENTRICITY:" or "EOT:" will not be qualified by the name ! of the gridded component. But for example, if the resource file specifies ! "EOT: .FALSE." ! but ! "SOLAR_EOT: .TRUE." ! then only SOLAR will have an EOT correction. The same goes for the new orbital ! system choice ORBIT_ANAL2B. ! A state's orbit is actually created in this routine by requesting the ORBIT ! object. If its not already created then it will be made below. GridComps that ! don't needed an orbit and dont request one will not have one. ! Parameters of standard orbital system (tabularized intercalation cycle) ! ----------------------------------------------------------------------- call ESMF_ConfigGetAttribute (CF, & ECC, label="ECCENTRICITY:", & default=DEFAULT_ORBIT_ECCENTRICITY, _RC) call ESMF_ConfigGetAttribute (CF, & OB, label="OBLIQUITY:", & default=DEFAULT_ORBIT_OBLIQUITY, _RC) call ESMF_ConfigGetAttribute (CF, & PER, label="PERIHELION:", & default=DEFAULT_ORBIT_PERIHELION, _RC) call ESMF_ConfigGetAttribute (CF, & EQNX, label="EQUINOX:", & default=DEFAULT_ORBIT_EQUINOX, _RC) ! Apply Equation of Time correction? ! ---------------------------------- call ESMF_ConfigGetAttribute (CF, & EOT, label="EOT:", & default=.FALSE., _RC) ! New orbital system (analytic two-body) allows some time-varying ! behavior, namely, linear variation in LAMBDAP, ECC, and OBQ. ! --------------------------------------------------------------- call ESMF_ConfigGetAttribute (CF, & ORBIT_ANAL2B, label="ORBIT_ANAL2B:", & default=.FALSE., _RC) ! Fixed anomalistic year length in mean solar days call ESMF_ConfigGetAttribute (CF, & ORB2B_YEARLEN, label="ORB2B_YEARLEN:", & default=DEFAULT_ORB2B_YEARLEN, _RC) ! Reference date and time for orbital parameters call ESMF_ConfigGetAttribute (CF, & ORB2B_REF_YYYYMMDD, label="ORB2B_REF_YYYYMMDD:", & default=DEFAULT_ORB2B_REF_YYYYMMDD, _RC) call ESMF_ConfigGetAttribute (CF, & ORB2B_REF_HHMMSS, label="ORB2B_REF_HHMMSS:", & default=DEFAULT_ORB2B_REF_HHMMSS, _RC) ! Orbital eccentricity at reference date call ESMF_ConfigGetAttribute (CF, & ORB2B_ECC_REF, label="ORB2B_ECC_REF:", & default=DEFAULT_ORB2B_ECC_REF, _RC) ! Rate of change of orbital eccentricity per Julian century call ESMF_ConfigGetAttribute (CF, & ORB2B_ECC_RATE, label="ORB2B_ECC_RATE:", & default=DEFAULT_ORB2B_ECC_RATE, _RC) ! Earth's obliquity (axial tilt) at reference date [degrees] call ESMF_ConfigGetAttribute (CF, & ORB2B_OBQ_REF, label="ORB2B_OBQ_REF:", & default=DEFAULT_ORB2B_OBQ_REF, _RC) ! Rate of change of obliquity [degrees per Julian century] call ESMF_ConfigGetAttribute (CF, & ORB2B_OBQ_RATE, label="ORB2B_OBQ_RATE:", & default=DEFAULT_ORB2B_OBQ_RATE, _RC) ! Longitude of perihelion at reference date [degrees] ! (from March equinox to perihelion in direction of earth's motion) call ESMF_ConfigGetAttribute (CF, & ORB2B_LAMBDAP_REF, label="ORB2B_LAMBDAP_REF:", & default=DEFAULT_ORB2B_LAMBDAP_REF, _RC) ! Rate of change of LAMBDAP [degrees per Julian century] ! (Combines both equatorial and ecliptic precession) call ESMF_ConfigGetAttribute (CF, & ORB2B_LAMBDAP_RATE, label="ORB2B_LAMBDAP_RATE:", & default=DEFAULT_ORB2B_LAMBDAP_RATE, _RC) ! March Equinox date and time call ESMF_ConfigGetAttribute (CF, & ORB2B_EQUINOX_YYYYMMDD, label="ORB2B_EQUINOX_YYYYMMDD:", & default=DEFAULT_ORB2B_EQUINOX_YYYYMMDD, _RC) call ESMF_ConfigGetAttribute (CF, & ORB2B_EQUINOX_HHMMSS, label="ORB2B_EQUINOX_HHMMSS:", & default=DEFAULT_ORB2B_EQUINOX_HHMMSS, _RC) ! create the orbit object ORBIT = MAPL_SunOrbitCreate ( & CLOCK, ECC, OB, PER, EQNX, & EOT, ORBIT_ANAL2B, ORB2B_YEARLEN, & ORB2B_REF_YYYYMMDD, ORB2B_REF_HHMMSS, & ORB2B_ECC_REF, ORB2B_ECC_RATE, & ORB2B_OBQ_REF, ORB2B_OBQ_RATE, & ORB2B_LAMBDAP_REF, ORB2B_LAMBDAP_RATE, & ORB2B_EQUINOX_YYYYMMDD, ORB2B_EQUINOX_HHMMSS, & FIX_SUN=FIX_SUN,_RC) _RETURN(ESMF_SUCCESS) end function MAPL_SunOrbitCreateFromConfig !========================================================================== !> ! The routine `MAPL_SunOrbitDestroy` Destroys a *GEOS_SunOrbit* object, ! deallocating the space used to save the ephemeris. ! subroutine MAPL_SunOrbitDestroy(ORBIT, RC) ! !ARGUMENTS: type(MAPL_SunOrbit), intent(INOUT) :: ORBIT integer, optional, intent( OUT) :: RC character(len=ESMF_MAXSTR), parameter :: IAm = "SunOrbitDestroy" if(associated(ORBIT%TH)) deallocate(ORBIT%TH) if(associated(ORBIT%ZC)) deallocate(ORBIT%ZC) if(associated(ORBIT%ZS)) deallocate(ORBIT%ZS) if(associated(ORBIT%PP)) deallocate(ORBIT%PP) if(associated(ORBIT%ET)) deallocate(ORBIT%ET) ORBIT%CREATED = .FALSE. _RETURN(ESMF_SUCCESS) end subroutine MAPL_SunOrbitDestroy !========================================================================== !> ! The function `MAPL_SunOrbitCreated` returns `.true.` ! if the given orbit object has been initilized. ! logical function MAPL_SunOrbitCreated(ORBIT, RC) ! !ARGUMENTS: type(MAPL_SunOrbit), intent(IN ) :: ORBIT integer, optional, intent(OUT) :: RC character(len=ESMF_MAXSTR), parameter :: IAm = "SunOrbitCreated" MAPL_SunOrbitCreated = ORBIT%CREATED _RETURN(ESMF_SUCCESS) return end function MAPL_SunOrbitCreated !========================================================================== !> ! Query for quantities in an orbit object. ! Optionally returns the parameters of the orbit and its ! associated `ESMF_Clock`. It fails ! if the orbit has not been created. ! ! @bug ! Not updated for ORBIT_ANAL2B option, which does not precalc ! many of the above outputs. !@endbug ! subroutine MAPL_SunOrbitQuery(ORBIT, & ECCENTRICITY, & OBLIQUITY, & PERIHELION, & EQUINOX, & YEAR_LENGTH, & YEARS_PER_CYCLE, & DAYS_PER_CYCLE, & CLOCK, & ZS, & ZC, & TH, & PP, & ET, & RC ) ! !ARGUMENTS: type(MAPL_SunOrbit), intent(IN ) :: ORBIT real, optional, intent(OUT) :: OBLIQUITY real, optional, intent(OUT) :: ECCENTRICITY real, optional, intent(OUT) :: PERIHELION real, optional, intent(OUT) :: YEAR_LENGTH integer, optional, intent(OUT) :: EQUINOX integer, optional, intent(OUT) :: YEARS_PER_CYCLE integer, optional, intent(OUT) :: DAYS_PER_CYCLE type(ESMF_Clock ), optional, intent(OUT) :: CLOCK real, optional, pointer, dimension(:) :: ZS real, optional, pointer, dimension(:) :: ZC real, optional, pointer, dimension(:) :: TH real, optional, pointer, dimension(:) :: PP real, optional, pointer, dimension(:) :: ET integer, optional, intent(OUT) :: RC character(len=ESMF_MAXSTR), parameter :: IAm = "SunOrbitQuery" integer :: STATUS _ASSERT(MAPL_SunOrbitCreated(ORBIT,RC=STATUS),'MAPL_SunOrbit not yet created!') if(present(CLOCK )) CLOCK = ORBIT%CLOCK if(present(OBLIQUITY )) OBLIQUITY = ORBIT%OB if(present(ECCENTRICITY )) ECCENTRICITY = ORBIT%ECC if(present(PERIHELION )) PERIHELION = ORBIT%PER if(present(EQUINOX )) EQUINOX = ORBIT%EQNX if(present(YEAR_LENGTH )) YEAR_LENGTH = ORBIT%YEARLEN if(present(DAYS_PER_CYCLE )) DAYS_PER_CYCLE = ORBIT%DAYS_PER_CYCLE if(present(YEARS_PER_CYCLE)) YEARS_PER_CYCLE = ORBIT%YEARS_PER_CYCLE if(present(ZS )) ZS => ORBIT%ZS if(present(ZC )) ZC => ORBIT%ZC if(present(TH )) TH => ORBIT%TH if(present(PP )) PP => ORBIT%PP if(present(ET )) ET => ORBIT%ET _RETURN(ESMF_SUCCESS) end subroutine MAPL_SunOrbitQuery !========================================================================== !> ! `GEOS_SunGetInsolation` returns the cosine of the solar zenith angle and the ! insolation at the top of the atmosphere for the given reference time, latitudes, ! longitudes, and orbit. It is overloaded to accept either 1d or 2d ! FORTRAN arrays or ESMF arrays of lats and lons and to produce the ! corresponding outputs. ! ! The reference time is obtained as follows. If CurrTime is specified, it is used; ! otherwise, if the optional clock is specified, it is set to the time ! on that clock. If neither currTime nor clock are given, the time on the attached ! clock is used. ! ! If the optional time interval is specified, the return values are averages ! over that interval following the reference time. In this case, the cosine of ! the solar zenith angle is an insolation-weighted average. The straight average ! of the zenith angle and the average over the daylight part of the ! interval (ZTHB,ZTHD) and the values at the beginning and end of the interval ! (ZTH1,ZTHN) are also optionally available. The last two are supported only for ! first two (non-ESMF) overloads. ! PMN Jun 2020: Added optional ZTHP, like ZTHB a pure average of the cosine of the ! solar zenith angle, except that for ZTHP the value averaged is allowed to be ! negative (below the horizon). It will be used for Photolysis calculations. ! Not implemented for ESMF overload. ! ! If the interval is not specified, the values are instantaneous values valid at ! the reference time. ! ! The optional *TIME* argument is used to return some specialized ! insolations. For example, the orbit at any of four Equinox or Solstice ! positions. If *TIME* is present, only the time of day is used from the clock, ! and a time interval, if specified, must be less than 24 hours. It can also be ! used to return daily-mean insolation for the date on the clock, or the annual-mean ! insolation for the year on the clock. ! ! The *TIME* argument can be any of the following: !``` ! MAPL_SunAutumnalEquinox ! MAPL_SunWinterSolstice ! MAPL_SunVernalEquinox ! MAPL_SunSummerSolstice ! MAPL_SunDailyMean ! MAPL_SunAnnualMean !``` ! ! @note ! If ORBIT%EOT is .TRUE., an Equation of Time correction will be ! applied. This shifts the actual daylight period w.r.t. to mean solar ! noon, to account for small but cumulative eccentricity and obliquity ! effects on the actual length of the solar day. !@endnote ! #undef DIMENSIONS #define DIMENSIONS (:) #define THE_SIZE (size(LONS,1)) recursive subroutine SOLAR_1D(LONS, LATS, ORBIT,ZTH,SLR,INTV,CLOCK, & TIME,currTime,DIST,ZTHB,ZTHD,ZTH1,ZTHN,ZTHP,& STEPSIZE,RC) #include "sun.H" end subroutine SOLAR_1D !========================================================================== #undef DIMENSIONS #undef THE_SIZE #define DIMENSIONS (:,:) #define THE_SIZE (size(LONS,1),size(LONS,2)) recursive subroutine SOLAR_2D(LONS, LATS, ORBIT,ZTH,SLR,INTV,CLOCK, & TIME,currTime,DIST,ZTHB,ZTHD,ZTH1,ZTHN,ZTHP,& STEPSIZE,RC) #include "sun.H" end subroutine SOLAR_2D #undef DIMENSIONS #undef THE_SIZE !BOP !EOP !========================================================================== subroutine SOLAR_ARR_INT(LONS, LATS, ORBIT, ZTH, SLR, INTV, CLOCK, & TIME, currTime, DIST, ZTHB, ZTHD, RC) type (MAPL_SunOrbit), intent(IN ) :: ORBIT type (ESMF_Array), intent(IN ) :: LATS type (ESMF_Array), intent(IN ) :: LONS type (ESMF_Array), intent(OUT) :: ZTH type (ESMF_Array), intent(OUT) :: SLR type (ESMF_TimeInterval), optional, intent(INout) :: INTV type (ESMF_Clock), optional, intent(IN ) :: CLOCK type (ESMF_Time), optional, intent(IN ) :: currTime integer, optional, intent(IN ) :: TIME real, optional, intent(OUT) :: DIST type (ESMF_Array), optional, intent(OUT) :: ZTHB type (ESMF_Array), optional, intent(OUT) :: ZTHD integer, optional, intent(OUT) :: RC ! Locals !!$ character(len=ESMF_MAXSTR) :: IAm = "SunGetInsolationArr" integer :: STATUS real, pointer, dimension (: ) :: LONS1, LATS1, ZTH1, SLR1, ZTHB1, ZTHD1 real, pointer, dimension (:,:) :: LONS2, LATS2, ZTH2, SLR2, ZTHB2, ZTHD2 integer :: RANK ! Begin _FAIL(" pmn: this routine is not up to date, is it even used anywhere?") call ESMF_ArrayGet(LONS, RANK=RANK, RC=STATUS) _VERIFY(STATUS) select case(RANK) case(1) call ESMF_ArrayGet(LATS, localDE=0, farrayptr=LATS1, RC=STATUS) _VERIFY(STATUS) call ESMF_ArrayGet(LONS,localDE=0, farrayptr=LONS1, RC=STATUS) _VERIFY(STATUS) call ESMF_ArrayGet(ZTH ,localDE=0, farrayptr=ZTH1, RC=STATUS) _VERIFY(STATUS) call ESMF_ArrayGet(SLR ,localDE=0, farrayptr=SLR1, RC=STATUS) _VERIFY(STATUS) if(present(ZTHB) .and. present(ZTHD)) then call ESMF_ArrayGet(ZTHB ,localDE=0, farrayptr=ZTHB1 ,RC=STATUS) _VERIFY(STATUS) call ESMF_ArrayGet(ZTHD ,localDE=0, farrayptr=ZTHD1 ,RC=STATUS) _VERIFY(STATUS) call MAPL_SunGetInsolation(LONS1,LATS1,ORBIT,ZTH1,SLR1,INTV,CLOCK,& TIME,currTime,DIST=DIST,ZTHB=ZTHB1,ZTHD=ZTHD1,RC=STATUS) elseif(present(ZTHB)) then call ESMF_ArrayGet(ZTHB ,localDE=0, farrayptr=ZTHB1 ,RC=STATUS) _VERIFY(STATUS) call MAPL_SunGetInsolation(LONS1,LATS1,ORBIT,ZTH1,SLR1,INTV,CLOCK,& TIME,currTime,DIST=DIST,ZTHB=ZTHB1,RC=STATUS) elseif(present(ZTHD)) then call ESMF_ArrayGet(ZTHD ,localDE=0, farrayptr=ZTHD1 ,RC=STATUS) _VERIFY(STATUS) call MAPL_SunGetInsolation(LONS1,LATS1,ORBIT,ZTH1,SLR1,INTV,CLOCK,& TIME,currTime,DIST=DIST,ZTHD=ZTHD1,RC=STATUS) else call MAPL_SunGetInsolation(LONS1,LATS1,ORBIT,ZTH1,SLR1,INTV,CLOCK,& TIME,currTime,DIST=DIST,RC=STATUS) endif _VERIFY(STATUS) case(2) call ESMF_ArrayGet(LATS,localDE=0, farrayptr=LATS2,RC=STATUS) _VERIFY(STATUS) call ESMF_ArrayGet(LONS,localDE=0, farrayptr=LONS2,RC=STATUS) _VERIFY(STATUS) call ESMF_ArrayGet(ZTH ,localDE=0, farrayptr=ZTH2 ,RC=STATUS) _VERIFY(STATUS) call ESMF_ArrayGet(SLR ,localDE=0, farrayptr=SLR2 ,RC=STATUS) _VERIFY(STATUS) if(present(ZTHB) .and. present(ZTHD)) then call ESMF_ArrayGet(ZTHB ,localDE=0, farrayptr=ZTHB2 ,RC=STATUS) _VERIFY(STATUS) call ESMF_ArrayGet(ZTHD ,localDE=0, farrayptr=ZTHD2 ,RC=STATUS) _VERIFY(STATUS) call MAPL_SunGetInsolation(LONS2,LATS2,ORBIT,ZTH2,SLR2,INTV,CLOCK,& TIME,currTime,DIST=DIST,ZTHB=ZTHB2,ZTHD=ZTHD2,RC=STATUS) elseif(present(ZTHB)) then call ESMF_ArrayGet(ZTHB ,localDE=0, farrayptr=ZTHB2 ,RC=STATUS) _VERIFY(STATUS) call MAPL_SunGetInsolation(LONS2,LATS2,ORBIT,ZTH2,SLR2,INTV,CLOCK,& TIME,currTime,DIST=DIST,ZTHB=ZTHB2,RC=STATUS) elseif(present(ZTHD)) then call ESMF_ArrayGet(ZTHD ,localDE=0, farrayptr=ZTHD2 ,RC=STATUS) _VERIFY(STATUS) call MAPL_SunGetInsolation(LONS2,LATS2,ORBIT,ZTH2,SLR2,INTV,CLOCK,& TIME,currTime,DIST=DIST,ZTHD=ZTHD2,RC=STATUS) else call MAPL_SunGetInsolation(LONS2,LATS2,ORBIT,ZTH2,SLR2,INTV,CLOCK,& TIME,currTime,DIST=DIST,RC=STATUS) endif _VERIFY(STATUS) case default _RETURN(ESMF_FAILURE) end select _RETURN(ESMF_SUCCESS) end subroutine SOLAR_ARR_INT subroutine GETIDAY(IDAY,TIME,ORBIT,RC) integer, intent(OUT) :: IDAY integer, intent(IN ) :: TIME type(MAPL_SunORBIT), intent(IN ) :: ORBIT integer, optional, intent(OUT) :: RC real :: ANOMALY select case(TIME) case (MAPL_SunAutumnalEquinox) ANOMALY = MAPL_PI case (MAPL_SunWinterSolstice ) ANOMALY = (MAPL_PI*3.0)/2.0 case (MAPL_SunVernalEquinox ) ANOMALY = 0.0 case (MAPL_SunSummerSolstice ) ANOMALY = MAPL_PI/2.0 case default _RETURN(ESMF_FAILURE) end select do IDAY=1,ORBIT%DAYS_PER_CYCLE-1 if(ORBIT%TH(IDAY)<=ANOMALY .and. ORBIT%TH(IDAY+1)>ANOMALY) then _RETURN(ESMF_SUCCESS) end if end do _RETURN(ESMF_FAILURE) end subroutine GETIDAY !========================================================================== subroutine MAPL_SunGetSolarConstantByTime(Time,SC,HK,rc) type (ESMF_Time), intent(INOUT) :: Time real, intent(OUT) :: SC real, optional, intent(OUT) :: HK(:) integer, optional, intent(OUT) :: rc integer :: YY, DOY integer :: STATUS call ESMF_TimeGet (TIME, YY=YY, DayOfYear=DOY, RC=STATUS) _VERIFY(STATUS) call MAPL_SunGetSolarConstantByYearDoY(YY,DOY,SC,HK, RC=STATUS) _VERIFY(STATUS) _RETURN(ESMF_SUCCESS) end subroutine MAPL_SunGetSolarConstantByTime !------------------------------------------------------------------------- !NASA/GSFC, Global Modeling and Assimilation Office, Code 610.1, GEOS/DAS! !------------------------------------------------------------------------- !> ! Given the year and day-of-year, this function returns the CMIP5 solar constant ! and the 8 band fractions for the Chou solar code. These are based on annual ! values from 1610 to 2008 and a repeating cycle after 2008. For dates prior ! tp 1 January 1610, it returns the value at the start of 1610. ! ! DayOfYear is expected to have a value of 1.00 at 0:00 UTC Jan 1 and leap years ! are acounted for by repeating the 366th day on January first of the following year. ! ! The SC values have been multiplied by .9965 to calibrate to the SOURCE/TIM scale. ! The SC values include the "background" variation. ! ! @bug ! Band values for RRTMG not implemented !@endbug ! subroutine MAPL_SunGetSolarConstantByYearDoY(year,dayofyear,SC,HK, rc) integer, intent(IN) :: Year integer, intent(IN) :: DayOfYear real, intent(OUT) :: SC real, optional, intent(OUT) :: HK(:) integer, optional, intent(OUT) :: rc real :: F integer :: i1,i2,Current integer, parameter :: firstYear = 1610 integer, parameter :: finalYear = 2008 integer, parameter :: SolCycle = 12 ! Cycle 23 was 12.2 years real, parameter :: SO(firstYear:finalYear) = (/ & 1360.768, 1360.751, 1361.147, 1361.259, 1361.333,& 1361.051, 1360.636, 1360.510, 1360.500, 1360.590,& 1360.586, 1360.582, 1360.578, 1360.574, 1360.534,& 1360.751, 1360.640, 1360.563, 1360.573, 1360.514,& 1360.365, 1360.373, 1360.340, 1360.438, 1360.358,& 1360.363, 1360.330, 1360.327, 1360.792, 1360.844,& 1360.415, 1360.304, 1360.622, 1360.404, 1360.349,& 1360.268, 1360.267, 1360.267, 1360.265, 1360.262,& 1360.261, 1360.261, 1360.288, 1360.266, 1360.258,& 1360.248, 1360.247, 1360.242, 1360.237, 1360.235,& 1360.248, 1360.240, 1360.235, 1360.235, 1360.235,& 1360.235, 1360.235, 1360.235, 1360.235, 1360.235,& 1360.235, 1360.241, 1360.237, 1360.235, 1360.236,& 1360.234, 1360.246, 1360.236, 1360.235, 1360.234,& 1360.239, 1360.234, 1360.234, 1360.234, 1360.244,& 1360.234, 1360.238, 1360.235, 1360.238, 1360.237,& 1360.237, 1360.238, 1360.240, 1360.240, 1360.240,& 1360.241, 1360.240, 1360.240, 1360.240, 1360.241,& 1360.244, 1360.246, 1360.249, 1360.267, 1360.282,& 1360.297, 1360.286, 1360.302, 1360.288, 1360.284,& 1360.286, 1360.296, 1360.302, 1360.308, 1360.325,& 1360.352, 1360.391, 1360.451, 1360.393, 1360.564,& 1360.507, 1360.510, 1360.443, 1360.406, 1360.490,& 1360.479, 1360.648, 1360.654, 1360.844, 1360.569,& 1360.881, 1360.403, 1360.530, 1360.420, 1360.430,& 1360.560, 1360.769, 1360.604, 1360.558, 1360.800,& 1360.503, 1360.831, 1360.546, 1360.491, 1360.436,& 1360.438, 1360.444, 1360.452, 1360.875, 1360.895,& 1360.859, 1360.689, 1360.661, 1360.643, 1360.556,& 1360.543, 1360.571, 1360.694, 1360.807, 1360.870,& 1360.846, 1361.010, 1360.872, 1360.802, 1360.778,& 1360.628, 1360.597, 1360.802, 1361.054, 1361.238,& 1361.241, 1361.135, 1361.063, 1360.845, 1360.807,& 1360.670, 1360.725, 1360.886, 1361.115, 1361.173,& 1360.995, 1361.106, 1360.842, 1360.749, 1360.618,& 1360.677, 1360.988, 1361.153, 1361.100, 1361.079,& 1360.980, 1360.830, 1360.817, 1360.796, 1360.702,& 1360.586, 1360.556, 1360.503, 1360.478, 1360.485,& 1360.525, 1360.796, 1360.675, 1360.543, 1360.547,& 1360.557, 1360.460, 1360.393, 1360.368, 1360.340,& 1360.332, 1360.344, 1360.385, 1360.434, 1360.457,& 1360.512, 1360.608, 1360.581, 1360.533, 1360.518,& 1360.463, 1360.441, 1360.465, 1360.470, 1360.511,& 1360.584, 1360.687, 1360.798, 1360.885, 1360.903,& 1360.943, 1360.779, 1360.678, 1360.582, 1360.618,& 1360.878, 1361.258, 1361.331, 1361.119, 1361.037,& 1360.921, 1360.770, 1360.715, 1360.645, 1360.675,& 1360.803, 1360.900, 1361.001, 1361.197, 1361.188,& 1360.990, 1360.994, 1360.915, 1360.805, 1360.682,& 1360.604, 1360.606, 1360.715, 1360.908, 1361.125,& 1361.191, 1361.082, 1360.932, 1360.849, 1360.787,& 1360.696, 1360.633, 1360.581, 1360.740, 1360.976,& 1361.213, 1361.154, 1361.107, 1360.907, 1360.779,& 1360.631, 1360.581, 1360.581, 1360.552, 1360.575,& 1360.721, 1360.865, 1360.888, 1360.935, 1360.837,& 1360.728, 1360.563, 1360.520, 1360.484, 1360.456,& 1360.491, 1360.741, 1360.839, 1360.980, 1361.075,& 1360.959, 1360.779, 1360.634, 1360.611, 1360.559,& 1360.529, 1360.451, 1360.459, 1360.669, 1360.938,& 1360.750, 1360.946, 1360.830, 1360.895, 1360.785,& 1360.652, 1360.569, 1360.523, 1360.539, 1360.605,& 1360.909, 1361.118, 1361.267, 1361.229, 1361.000,& 1360.838, 1360.733, 1360.616, 1360.648, 1360.679,& 1360.882, 1360.983, 1361.166, 1361.044, 1361.003,& 1360.985, 1360.864, 1360.757, 1360.637, 1360.748,& 1360.964, 1361.352, 1361.286, 1361.222, 1361.206,& 1361.143, 1361.065, 1360.962, 1360.805, 1360.834,& 1361.100, 1361.198, 1361.437, 1361.567, 1361.474,& 1361.232, 1360.996, 1360.987, 1360.849, 1360.877,& 1360.997, 1361.529, 1361.885, 1361.850, 1361.601,& 1361.495, 1361.139, 1360.968, 1360.916, 1360.918,& 1360.954, 1361.137, 1361.333, 1361.383, 1361.466,& 1361.461, 1361.177, 1361.271, 1361.019, 1360.947,& 1360.755, 1360.866, 1361.053, 1361.493, 1361.852,& 1361.865, 1361.912, 1361.504, 1361.418, 1361.030,& 1360.862, 1360.858, 1361.010, 1361.301, 1361.865,& 1361.770, 1361.663, 1361.520, 1361.247, 1361.017,& 1360.920, 1360.832, 1360.960, 1361.321, 1361.603,& 1361.882, 1361.819, 1361.897, 1361.448, 1361.267,& 1361.074, 1361.030, 1360.944, 1360.912 /) real, parameter :: CHOUBAND1(firstYear:finalYear) = (/ & 0.00516, 0.00515, 0.00519, 0.00519, 0.00520,& 0.00518, 0.00515, 0.00514, 0.00514, 0.00514,& 0.00514, 0.00514, 0.00514, 0.00514, 0.00514,& 0.00516, 0.00515, 0.00514, 0.00514, 0.00514,& 0.00513, 0.00513, 0.00513, 0.00514, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00517, 0.00517,& 0.00514, 0.00513, 0.00515, 0.00514, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00514, 0.00513, 0.00515,& 0.00514, 0.00514, 0.00514, 0.00513, 0.00514,& 0.00514, 0.00515, 0.00515, 0.00517, 0.00514,& 0.00517, 0.00513, 0.00514, 0.00513, 0.00513,& 0.00514, 0.00516, 0.00515, 0.00514, 0.00516,& 0.00514, 0.00516, 0.00514, 0.00514, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00517, 0.00517,& 0.00517, 0.00515, 0.00515, 0.00515, 0.00514,& 0.00514, 0.00514, 0.00515, 0.00516, 0.00516,& 0.00516, 0.00517, 0.00516, 0.00516, 0.00515,& 0.00514, 0.00514, 0.00516, 0.00518, 0.00519,& 0.00519, 0.00518, 0.00518, 0.00516, 0.00515,& 0.00514, 0.00515, 0.00516, 0.00518, 0.00518,& 0.00517, 0.00518, 0.00516, 0.00515, 0.00514,& 0.00515, 0.00517, 0.00518, 0.00518, 0.00518,& 0.00517, 0.00516, 0.00516, 0.00516, 0.00515,& 0.00514, 0.00514, 0.00514, 0.00514, 0.00514,& 0.00514, 0.00516, 0.00515, 0.00514, 0.00514,& 0.00514, 0.00514, 0.00513, 0.00513, 0.00513,& 0.00513, 0.00513, 0.00513, 0.00514, 0.00514,& 0.00514, 0.00515, 0.00515, 0.00514, 0.00514,& 0.00514, 0.00513, 0.00514, 0.00513, 0.00514,& 0.00514, 0.00515, 0.00516, 0.00517, 0.00517,& 0.00517, 0.00516, 0.00515, 0.00514, 0.00514,& 0.00516, 0.00519, 0.00520, 0.00518, 0.00517,& 0.00517, 0.00515, 0.00515, 0.00514, 0.00515,& 0.00516, 0.00516, 0.00517, 0.00519, 0.00518,& 0.00517, 0.00517, 0.00517, 0.00516, 0.00515,& 0.00514, 0.00514, 0.00515, 0.00516, 0.00518,& 0.00519, 0.00518, 0.00517, 0.00516, 0.00516,& 0.00515, 0.00514, 0.00514, 0.00515, 0.00517,& 0.00519, 0.00518, 0.00518, 0.00517, 0.00516,& 0.00515, 0.00514, 0.00514, 0.00514, 0.00514,& 0.00515, 0.00516, 0.00516, 0.00517, 0.00517,& 0.00516, 0.00515, 0.00514, 0.00514, 0.00514,& 0.00514, 0.00516, 0.00517, 0.00518, 0.00518,& 0.00517, 0.00516, 0.00515, 0.00515, 0.00514,& 0.00514, 0.00514, 0.00514, 0.00515, 0.00516,& 0.00517, 0.00517, 0.00517, 0.00517, 0.00516,& 0.00515, 0.00514, 0.00514, 0.00514, 0.00515,& 0.00517, 0.00518, 0.00520, 0.00519, 0.00518,& 0.00516, 0.00516, 0.00515, 0.00515, 0.00515,& 0.00517, 0.00518, 0.00518, 0.00519, 0.00518,& 0.00517, 0.00516, 0.00515, 0.00515, 0.00515,& 0.00517, 0.00519, 0.00521, 0.00520, 0.00520,& 0.00518, 0.00517, 0.00517, 0.00516, 0.00516,& 0.00517, 0.00520, 0.00522, 0.00522, 0.00522,& 0.00519, 0.00518, 0.00517, 0.00516, 0.00516,& 0.00517, 0.00522, 0.00525, 0.00525, 0.00523,& 0.00521, 0.00517, 0.00517, 0.00516, 0.00516,& 0.00516, 0.00517, 0.00520, 0.00520, 0.00520,& 0.00521, 0.00518, 0.00519, 0.00517, 0.00516,& 0.00515, 0.00516, 0.00517, 0.00520, 0.00523,& 0.00523, 0.00523, 0.00522, 0.00519, 0.00517,& 0.00516, 0.00516, 0.00516, 0.00519, 0.00523,& 0.00522, 0.00522, 0.00520, 0.00518, 0.00517,& 0.00516, 0.00515, 0.00516, 0.00518, 0.00520,& 0.00522, 0.00522, 0.00522, 0.00519, 0.00518,& 0.00517, 0.00516, 0.00516, 0.00516 /) real, parameter :: CHOUBAND2(firstYear:finalYear) = (/ & 0.00567, 0.00567, 0.00569, 0.00570, 0.00570,& 0.00569, 0.00567, 0.00566, 0.00566, 0.00567,& 0.00567, 0.00567, 0.00567, 0.00567, 0.00566,& 0.00567, 0.00567, 0.00567, 0.00567, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00568, 0.00568,& 0.00566, 0.00566, 0.00567, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00567,& 0.00567, 0.00567, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00567, 0.00567, 0.00568, 0.00567,& 0.00568, 0.00566, 0.00567, 0.00566, 0.00566,& 0.00567, 0.00568, 0.00567, 0.00567, 0.00568,& 0.00566, 0.00568, 0.00567, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00568, 0.00568,& 0.00568, 0.00567, 0.00567, 0.00567, 0.00566,& 0.00566, 0.00566, 0.00567, 0.00568, 0.00568,& 0.00568, 0.00568, 0.00568, 0.00567, 0.00567,& 0.00567, 0.00566, 0.00567, 0.00569, 0.00569,& 0.00569, 0.00569, 0.00569, 0.00568, 0.00567,& 0.00567, 0.00567, 0.00568, 0.00569, 0.00569,& 0.00568, 0.00569, 0.00568, 0.00567, 0.00567,& 0.00567, 0.00568, 0.00569, 0.00569, 0.00569,& 0.00568, 0.00568, 0.00568, 0.00567, 0.00567,& 0.00567, 0.00567, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00568, 0.00567, 0.00567, 0.00567,& 0.00567, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00566, 0.00567, 0.00567, 0.00567, 0.00566,& 0.00566, 0.00566, 0.00566, 0.00566, 0.00566,& 0.00567, 0.00567, 0.00568, 0.00568, 0.00568,& 0.00568, 0.00567, 0.00567, 0.00566, 0.00567,& 0.00568, 0.00569, 0.00570, 0.00569, 0.00568,& 0.00568, 0.00567, 0.00567, 0.00567, 0.00567,& 0.00567, 0.00568, 0.00568, 0.00569, 0.00569,& 0.00568, 0.00568, 0.00568, 0.00567, 0.00567,& 0.00567, 0.00567, 0.00567, 0.00568, 0.00569,& 0.00569, 0.00569, 0.00568, 0.00568, 0.00567,& 0.00567, 0.00567, 0.00566, 0.00567, 0.00568,& 0.00569, 0.00569, 0.00569, 0.00568, 0.00567,& 0.00567, 0.00567, 0.00567, 0.00566, 0.00566,& 0.00567, 0.00568, 0.00568, 0.00568, 0.00568,& 0.00568, 0.00567, 0.00567, 0.00566, 0.00566,& 0.00567, 0.00568, 0.00568, 0.00569, 0.00569,& 0.00568, 0.00568, 0.00567, 0.00567, 0.00567,& 0.00567, 0.00566, 0.00566, 0.00567, 0.00568,& 0.00568, 0.00568, 0.00568, 0.00568, 0.00568,& 0.00567, 0.00567, 0.00567, 0.00567, 0.00567,& 0.00568, 0.00569, 0.00570, 0.00569, 0.00569,& 0.00568, 0.00567, 0.00567, 0.00567, 0.00567,& 0.00568, 0.00569, 0.00569, 0.00569, 0.00569,& 0.00568, 0.00568, 0.00567, 0.00567, 0.00567,& 0.00568, 0.00570, 0.00570, 0.00570, 0.00570,& 0.00569, 0.00568, 0.00568, 0.00568, 0.00567,& 0.00568, 0.00570, 0.00571, 0.00571, 0.00571,& 0.00569, 0.00569, 0.00568, 0.00568, 0.00567,& 0.00568, 0.00571, 0.00573, 0.00573, 0.00572,& 0.00570, 0.00569, 0.00568, 0.00568, 0.00568,& 0.00568, 0.00569, 0.00570, 0.00570, 0.00570,& 0.00570, 0.00569, 0.00569, 0.00568, 0.00568,& 0.00567, 0.00568, 0.00568, 0.00570, 0.00572,& 0.00572, 0.00572, 0.00571, 0.00570, 0.00568,& 0.00568, 0.00567, 0.00568, 0.00570, 0.00572,& 0.00571, 0.00571, 0.00570, 0.00569, 0.00568,& 0.00568, 0.00567, 0.00568, 0.00569, 0.00570,& 0.00571, 0.00571, 0.00572, 0.00570, 0.00569,& 0.00568, 0.00568, 0.00568, 0.00567 /) real, parameter :: CHOUBAND3(firstYear:finalYear) = (/ & 0.01202, 0.01202, 0.01203, 0.01203, 0.01203,& 0.01202, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01202, 0.01202, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01202, 0.01202,& 0.01201, 0.01201, 0.01202, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01202, 0.01202, 0.01202, 0.01201,& 0.01202, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01202, 0.01201, 0.01201, 0.01202,& 0.01201, 0.01202, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01202, 0.01202,& 0.01202, 0.01202, 0.01202, 0.01202, 0.01201,& 0.01201, 0.01201, 0.01202, 0.01202, 0.01202,& 0.01202, 0.01202, 0.01202, 0.01202, 0.01202,& 0.01201, 0.01201, 0.01202, 0.01202, 0.01203,& 0.01203, 0.01203, 0.01202, 0.01202, 0.01202,& 0.01201, 0.01202, 0.01202, 0.01203, 0.01203,& 0.01202, 0.01202, 0.01202, 0.01202, 0.01201,& 0.01202, 0.01202, 0.01203, 0.01203, 0.01202,& 0.01202, 0.01202, 0.01202, 0.01202, 0.01202,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01202, 0.01202, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01202, 0.01202, 0.01202, 0.01202,& 0.01202, 0.01202, 0.01202, 0.01201, 0.01201,& 0.01202, 0.01203, 0.01203, 0.01203, 0.01202,& 0.01202, 0.01202, 0.01202, 0.01201, 0.01202,& 0.01202, 0.01202, 0.01202, 0.01203, 0.01203,& 0.01202, 0.01202, 0.01202, 0.01202, 0.01202,& 0.01201, 0.01201, 0.01202, 0.01202, 0.01203,& 0.01203, 0.01202, 0.01202, 0.01202, 0.01202,& 0.01202, 0.01201, 0.01201, 0.01202, 0.01202,& 0.01203, 0.01203, 0.01203, 0.01202, 0.01202,& 0.01201, 0.01201, 0.01201, 0.01201, 0.01201,& 0.01202, 0.01202, 0.01202, 0.01202, 0.01202,& 0.01202, 0.01202, 0.01201, 0.01201, 0.01201,& 0.01201, 0.01202, 0.01202, 0.01203, 0.01203,& 0.01202, 0.01202, 0.01202, 0.01202, 0.01201,& 0.01201, 0.01201, 0.01201, 0.01202, 0.01202,& 0.01202, 0.01202, 0.01202, 0.01202, 0.01202,& 0.01202, 0.01201, 0.01201, 0.01201, 0.01202,& 0.01202, 0.01203, 0.01203, 0.01203, 0.01202,& 0.01202, 0.01202, 0.01202, 0.01202, 0.01202,& 0.01202, 0.01202, 0.01203, 0.01203, 0.01203,& 0.01202, 0.01202, 0.01202, 0.01202, 0.01202,& 0.01202, 0.01203, 0.01203, 0.01203, 0.01203,& 0.01203, 0.01202, 0.01202, 0.01202, 0.01202,& 0.01202, 0.01203, 0.01204, 0.01204, 0.01204,& 0.01203, 0.01202, 0.01202, 0.01202, 0.01202,& 0.01202, 0.01204, 0.01205, 0.01205, 0.01204,& 0.01203, 0.01203, 0.01202, 0.01202, 0.01202,& 0.01202, 0.01203, 0.01203, 0.01203, 0.01203,& 0.01203, 0.01203, 0.01203, 0.01202, 0.01202,& 0.01202, 0.01202, 0.01202, 0.01203, 0.01204,& 0.01204, 0.01204, 0.01204, 0.01203, 0.01202,& 0.01202, 0.01202, 0.01202, 0.01203, 0.01204,& 0.01204, 0.01204, 0.01203, 0.01203, 0.01202,& 0.01202, 0.01202, 0.01202, 0.01203, 0.01203,& 0.01204, 0.01204, 0.01204, 0.01203, 0.01203,& 0.01202, 0.01202, 0.01202, 0.01202 /) real, parameter :: CHOUBAND4(firstYear:finalYear) = (/ & 0.05687, 0.05687, 0.05690, 0.05690, 0.05691,& 0.05689, 0.05686, 0.05685, 0.05685, 0.05686,& 0.05686, 0.05686, 0.05686, 0.05686, 0.05685,& 0.05687, 0.05686, 0.05686, 0.05686, 0.05685,& 0.05684, 0.05685, 0.05684, 0.05685, 0.05684,& 0.05685, 0.05684, 0.05684, 0.05687, 0.05688,& 0.05685, 0.05684, 0.05686, 0.05685, 0.05685,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05684, 0.05684, 0.05684, 0.05684,& 0.05684, 0.05685, 0.05685, 0.05685, 0.05686,& 0.05685, 0.05685, 0.05685, 0.05685, 0.05685,& 0.05685, 0.05686, 0.05686, 0.05688, 0.05686,& 0.05688, 0.05685, 0.05686, 0.05685, 0.05685,& 0.05686, 0.05687, 0.05686, 0.05686, 0.05687,& 0.05685, 0.05688, 0.05686, 0.05685, 0.05685,& 0.05685, 0.05685, 0.05685, 0.05688, 0.05688,& 0.05688, 0.05687, 0.05686, 0.05686, 0.05686,& 0.05685, 0.05686, 0.05686, 0.05687, 0.05688,& 0.05687, 0.05689, 0.05688, 0.05687, 0.05687,& 0.05686, 0.05686, 0.05687, 0.05689, 0.05690,& 0.05690, 0.05689, 0.05689, 0.05687, 0.05687,& 0.05686, 0.05686, 0.05688, 0.05689, 0.05690,& 0.05688, 0.05689, 0.05687, 0.05687, 0.05686,& 0.05686, 0.05688, 0.05690, 0.05689, 0.05689,& 0.05688, 0.05687, 0.05687, 0.05687, 0.05687,& 0.05686, 0.05686, 0.05685, 0.05685, 0.05685,& 0.05685, 0.05687, 0.05686, 0.05686, 0.05686,& 0.05686, 0.05685, 0.05685, 0.05685, 0.05684,& 0.05684, 0.05684, 0.05685, 0.05685, 0.05685,& 0.05685, 0.05686, 0.05686, 0.05686, 0.05685,& 0.05685, 0.05685, 0.05685, 0.05685, 0.05685,& 0.05686, 0.05686, 0.05687, 0.05688, 0.05688,& 0.05688, 0.05687, 0.05686, 0.05686, 0.05686,& 0.05688, 0.05690, 0.05691, 0.05689, 0.05689,& 0.05688, 0.05687, 0.05686, 0.05686, 0.05686,& 0.05687, 0.05688, 0.05688, 0.05690, 0.05690,& 0.05688, 0.05688, 0.05688, 0.05687, 0.05686,& 0.05686, 0.05686, 0.05686, 0.05688, 0.05689,& 0.05690, 0.05689, 0.05688, 0.05687, 0.05687,& 0.05686, 0.05686, 0.05686, 0.05687, 0.05688,& 0.05690, 0.05690, 0.05689, 0.05688, 0.05687,& 0.05686, 0.05686, 0.05686, 0.05685, 0.05686,& 0.05687, 0.05688, 0.05688, 0.05688, 0.05688,& 0.05687, 0.05686, 0.05686, 0.05686, 0.05685,& 0.05686, 0.05687, 0.05688, 0.05689, 0.05690,& 0.05689, 0.05687, 0.05687, 0.05687, 0.05686,& 0.05686, 0.05685, 0.05686, 0.05687, 0.05688,& 0.05688, 0.05689, 0.05688, 0.05688, 0.05688,& 0.05687, 0.05686, 0.05686, 0.05686, 0.05686,& 0.05688, 0.05689, 0.05691, 0.05690, 0.05689,& 0.05688, 0.05687, 0.05686, 0.05686, 0.05687,& 0.05688, 0.05689, 0.05690, 0.05690, 0.05689,& 0.05688, 0.05688, 0.05687, 0.05686, 0.05687,& 0.05688, 0.05691, 0.05691, 0.05691, 0.05691,& 0.05690, 0.05689, 0.05688, 0.05687, 0.05687,& 0.05689, 0.05691, 0.05693, 0.05693, 0.05692,& 0.05690, 0.05689, 0.05688, 0.05688, 0.05688,& 0.05689, 0.05693, 0.05695, 0.05695, 0.05694,& 0.05692, 0.05689, 0.05688, 0.05688, 0.05688,& 0.05688, 0.05689, 0.05691, 0.05691, 0.05692,& 0.05692, 0.05690, 0.05690, 0.05689, 0.05688,& 0.05687, 0.05688, 0.05689, 0.05692, 0.05694,& 0.05694, 0.05695, 0.05693, 0.05691, 0.05689,& 0.05688, 0.05687, 0.05688, 0.05691, 0.05694,& 0.05694, 0.05693, 0.05692, 0.05690, 0.05688,& 0.05688, 0.05687, 0.05688, 0.05690, 0.05692,& 0.05694, 0.05694, 0.05694, 0.05691, 0.05690,& 0.05689, 0.05689, 0.05688, 0.05688 /) real, parameter :: CHOUBAND5(firstYear:finalYear) = (/ & 0.37831, 0.37831, 0.37831, 0.37832, 0.37832,& 0.37831, 0.37830, 0.37830, 0.37830, 0.37830,& 0.37830, 0.37830, 0.37830, 0.37830, 0.37830,& 0.37831, 0.37830, 0.37830, 0.37830, 0.37830,& 0.37829, 0.37829, 0.37829, 0.37830, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37830, 0.37830,& 0.37829, 0.37829, 0.37830, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37829, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37830, 0.37829, 0.37830,& 0.37830, 0.37830, 0.37830, 0.37830, 0.37830,& 0.37830, 0.37830, 0.37830, 0.37831, 0.37830,& 0.37831, 0.37830, 0.37830, 0.37830, 0.37830,& 0.37830, 0.37831, 0.37830, 0.37830, 0.37831,& 0.37830, 0.37831, 0.37830, 0.37830, 0.37830,& 0.37830, 0.37830, 0.37830, 0.37831, 0.37831,& 0.37831, 0.37830, 0.37830, 0.37830, 0.37830,& 0.37830, 0.37830, 0.37831, 0.37831, 0.37831,& 0.37831, 0.37831, 0.37831, 0.37831, 0.37831,& 0.37831, 0.37831, 0.37831, 0.37831, 0.37832,& 0.37832, 0.37832, 0.37832, 0.37831, 0.37831,& 0.37831, 0.37831, 0.37831, 0.37832, 0.37832,& 0.37831, 0.37832, 0.37831, 0.37831, 0.37831,& 0.37831, 0.37831, 0.37832, 0.37831, 0.37831,& 0.37831, 0.37831, 0.37831, 0.37831, 0.37831,& 0.37830, 0.37830, 0.37830, 0.37830, 0.37830,& 0.37830, 0.37831, 0.37830, 0.37830, 0.37830,& 0.37830, 0.37830, 0.37830, 0.37829, 0.37829,& 0.37829, 0.37829, 0.37830, 0.37830, 0.37830,& 0.37830, 0.37830, 0.37830, 0.37830, 0.37830,& 0.37830, 0.37830, 0.37830, 0.37830, 0.37830,& 0.37830, 0.37830, 0.37831, 0.37831, 0.37831,& 0.37831, 0.37831, 0.37831, 0.37830, 0.37831,& 0.37831, 0.37832, 0.37832, 0.37832, 0.37831,& 0.37831, 0.37831, 0.37831, 0.37831, 0.37831,& 0.37831, 0.37831, 0.37831, 0.37832, 0.37832,& 0.37831, 0.37831, 0.37831, 0.37831, 0.37831,& 0.37831, 0.37831, 0.37831, 0.37831, 0.37832,& 0.37832, 0.37832, 0.37831, 0.37831, 0.37831,& 0.37831, 0.37831, 0.37830, 0.37831, 0.37831,& 0.37832, 0.37832, 0.37832, 0.37831, 0.37831,& 0.37830, 0.37830, 0.37830, 0.37830, 0.37830,& 0.37831, 0.37831, 0.37831, 0.37831, 0.37831,& 0.37830, 0.37830, 0.37830, 0.37830, 0.37830,& 0.37830, 0.37831, 0.37831, 0.37831, 0.37832,& 0.37831, 0.37831, 0.37830, 0.37830, 0.37830,& 0.37830, 0.37830, 0.37830, 0.37831, 0.37832,& 0.37830, 0.37832, 0.37831, 0.37831, 0.37831,& 0.37831, 0.37830, 0.37830, 0.37830, 0.37831,& 0.37832, 0.37832, 0.37832, 0.37833, 0.37832,& 0.37831, 0.37831, 0.37830, 0.37831, 0.37831,& 0.37831, 0.37831, 0.37832, 0.37831, 0.37831,& 0.37832, 0.37832, 0.37831, 0.37831, 0.37831,& 0.37832, 0.37833, 0.37832, 0.37832, 0.37832,& 0.37832, 0.37832, 0.37832, 0.37831, 0.37832,& 0.37833, 0.37832, 0.37832, 0.37833, 0.37833,& 0.37832, 0.37831, 0.37832, 0.37832, 0.37832,& 0.37832, 0.37833, 0.37834, 0.37833, 0.37832,& 0.37833, 0.37833, 0.37832, 0.37832, 0.37832,& 0.37832, 0.37833, 0.37833, 0.37833, 0.37833,& 0.37833, 0.37832, 0.37833, 0.37832, 0.37832,& 0.37831, 0.37832, 0.37833, 0.37834, 0.37835,& 0.37835, 0.37835, 0.37833, 0.37834, 0.37832,& 0.37832, 0.37832, 0.37832, 0.37833, 0.37835,& 0.37834, 0.37833, 0.37834, 0.37833, 0.37832,& 0.37832, 0.37832, 0.37832, 0.37833, 0.37834,& 0.37835, 0.37835, 0.37835, 0.37833, 0.37833,& 0.37832, 0.37832, 0.37832, 0.37832 /) real, parameter :: CHOUBAND6(firstYear:finalYear) = (/ & 0.33338, 0.33338, 0.33334, 0.33332, 0.33332,& 0.33335, 0.33339, 0.33340, 0.33341, 0.33340,& 0.33340, 0.33340, 0.33340, 0.33340, 0.33340,& 0.33338, 0.33339, 0.33340, 0.33340, 0.33340,& 0.33342, 0.33341, 0.33342, 0.33341, 0.33342,& 0.33341, 0.33342, 0.33342, 0.33337, 0.33336,& 0.33341, 0.33342, 0.33339, 0.33341, 0.33341,& 0.33342, 0.33342, 0.33342, 0.33342, 0.33342,& 0.33342, 0.33342, 0.33342, 0.33342, 0.33342,& 0.33342, 0.33342, 0.33342, 0.33343, 0.33343,& 0.33342, 0.33342, 0.33343, 0.33343, 0.33343,& 0.33343, 0.33343, 0.33343, 0.33343, 0.33343,& 0.33343, 0.33342, 0.33342, 0.33343, 0.33343,& 0.33343, 0.33342, 0.33343, 0.33343, 0.33343,& 0.33342, 0.33343, 0.33343, 0.33343, 0.33342,& 0.33343, 0.33342, 0.33343, 0.33342, 0.33343,& 0.33343, 0.33342, 0.33342, 0.33342, 0.33342,& 0.33342, 0.33342, 0.33342, 0.33342, 0.33342,& 0.33342, 0.33342, 0.33342, 0.33342, 0.33342,& 0.33342, 0.33342, 0.33342, 0.33342, 0.33342,& 0.33342, 0.33342, 0.33342, 0.33342, 0.33342,& 0.33342, 0.33341, 0.33341, 0.33341, 0.33339,& 0.33340, 0.33340, 0.33341, 0.33341, 0.33340,& 0.33340, 0.33339, 0.33339, 0.33337, 0.33340,& 0.33336, 0.33341, 0.33340, 0.33341, 0.33341,& 0.33340, 0.33337, 0.33339, 0.33340, 0.33337,& 0.33340, 0.33337, 0.33340, 0.33340, 0.33341,& 0.33341, 0.33341, 0.33341, 0.33336, 0.33336,& 0.33337, 0.33338, 0.33339, 0.33339, 0.33340,& 0.33340, 0.33340, 0.33339, 0.33337, 0.33337,& 0.33337, 0.33335, 0.33337, 0.33338, 0.33338,& 0.33339, 0.33340, 0.33338, 0.33335, 0.33333,& 0.33333, 0.33334, 0.33335, 0.33337, 0.33338,& 0.33339, 0.33339, 0.33337, 0.33334, 0.33334,& 0.33336, 0.33335, 0.33337, 0.33338, 0.33340,& 0.33339, 0.33336, 0.33334, 0.33334, 0.33334,& 0.33336, 0.33337, 0.33337, 0.33338, 0.33338,& 0.33340, 0.33340, 0.33340, 0.33341, 0.33341,& 0.33340, 0.33337, 0.33338, 0.33340, 0.33340,& 0.33340, 0.33341, 0.33341, 0.33341, 0.33342,& 0.33342, 0.33342, 0.33341, 0.33341, 0.33341,& 0.33340, 0.33339, 0.33339, 0.33340, 0.33340,& 0.33341, 0.33341, 0.33341, 0.33341, 0.33340,& 0.33340, 0.33339, 0.33337, 0.33336, 0.33336,& 0.33336, 0.33338, 0.33339, 0.33340, 0.33340,& 0.33337, 0.33333, 0.33332, 0.33334, 0.33335,& 0.33336, 0.33338, 0.33339, 0.33339, 0.33339,& 0.33338, 0.33337, 0.33336, 0.33333, 0.33334,& 0.33336, 0.33336, 0.33336, 0.33338, 0.33339,& 0.33340, 0.33340, 0.33339, 0.33337, 0.33334,& 0.33334, 0.33335, 0.33336, 0.33337, 0.33338,& 0.33339, 0.33339, 0.33340, 0.33338, 0.33336,& 0.33333, 0.33334, 0.33334, 0.33336, 0.33338,& 0.33339, 0.33340, 0.33340, 0.33340, 0.33340,& 0.33338, 0.33337, 0.33337, 0.33336, 0.33336,& 0.33337, 0.33339, 0.33340, 0.33340, 0.33340,& 0.33340, 0.33337, 0.33336, 0.33334, 0.33334,& 0.33335, 0.33337, 0.33338, 0.33339, 0.33339,& 0.33340, 0.33340, 0.33340, 0.33338, 0.33336,& 0.33336, 0.33336, 0.33336, 0.33336, 0.33337,& 0.33338, 0.33339, 0.33340, 0.33340, 0.33339,& 0.33336, 0.33334, 0.33332, 0.33333, 0.33335,& 0.33337, 0.33338, 0.33339, 0.33339, 0.33338,& 0.33336, 0.33334, 0.33333, 0.33334, 0.33334,& 0.33336, 0.33337, 0.33338, 0.33339, 0.33338,& 0.33336, 0.33332, 0.33331, 0.33331, 0.33332,& 0.33334, 0.33335, 0.33336, 0.33337, 0.33337,& 0.33335, 0.33332, 0.33329, 0.33329, 0.33329,& 0.33333, 0.33334, 0.33336, 0.33337, 0.33337,& 0.33336, 0.33329, 0.33325, 0.33325, 0.33327,& 0.33330, 0.33334, 0.33336, 0.33337, 0.33337,& 0.33337, 0.33335, 0.33332, 0.33331, 0.33331,& 0.33330, 0.33333, 0.33333, 0.33336, 0.33336,& 0.33338, 0.33337, 0.33336, 0.33331, 0.33327,& 0.33327, 0.33326, 0.33329, 0.33332, 0.33335,& 0.33337, 0.33337, 0.33336, 0.33332, 0.33327,& 0.33328, 0.33328, 0.33331, 0.33334, 0.33336,& 0.33337, 0.33338, 0.33336, 0.33333, 0.33330,& 0.33328, 0.33328, 0.33327, 0.33332, 0.33334,& 0.33335, 0.33336, 0.33337, 0.33337 /) real, parameter :: CHOUBAND7(firstYear:finalYear) = (/ & 0.16512, 0.16512, 0.16508, 0.16507, 0.16506,& 0.16509, 0.16513, 0.16515, 0.16515, 0.16514,& 0.16514, 0.16514, 0.16514, 0.16514, 0.16515,& 0.16512, 0.16513, 0.16514, 0.16514, 0.16515,& 0.16516, 0.16516, 0.16516, 0.16515, 0.16516,& 0.16516, 0.16517, 0.16517, 0.16511, 0.16511,& 0.16516, 0.16517, 0.16513, 0.16516, 0.16516,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16517, 0.16517, 0.16517, 0.16517, 0.16517,& 0.16516, 0.16516, 0.16515, 0.16516, 0.16514,& 0.16515, 0.16515, 0.16515, 0.16516, 0.16515,& 0.16515, 0.16513, 0.16513, 0.16511, 0.16514,& 0.16511, 0.16516, 0.16514, 0.16516, 0.16516,& 0.16514, 0.16512, 0.16514, 0.16514, 0.16512,& 0.16515, 0.16511, 0.16514, 0.16515, 0.16516,& 0.16516, 0.16515, 0.16515, 0.16511, 0.16510,& 0.16511, 0.16513, 0.16513, 0.16513, 0.16514,& 0.16515, 0.16514, 0.16513, 0.16512, 0.16511,& 0.16511, 0.16509, 0.16511, 0.16512, 0.16512,& 0.16514, 0.16514, 0.16512, 0.16509, 0.16507,& 0.16507, 0.16508, 0.16509, 0.16511, 0.16512,& 0.16513, 0.16513, 0.16511, 0.16508, 0.16508,& 0.16510, 0.16508, 0.16511, 0.16512, 0.16514,& 0.16513, 0.16510, 0.16508, 0.16508, 0.16509,& 0.16510, 0.16511, 0.16512, 0.16512, 0.16513,& 0.16514, 0.16514, 0.16515, 0.16515, 0.16515,& 0.16515, 0.16512, 0.16513, 0.16514, 0.16514,& 0.16514, 0.16515, 0.16516, 0.16516, 0.16516,& 0.16517, 0.16516, 0.16516, 0.16515, 0.16515,& 0.16515, 0.16514, 0.16514, 0.16514, 0.16515,& 0.16515, 0.16515, 0.16515, 0.16515, 0.16515,& 0.16514, 0.16513, 0.16512, 0.16511, 0.16510,& 0.16510, 0.16512, 0.16513, 0.16514, 0.16514,& 0.16511, 0.16507, 0.16506, 0.16508, 0.16509,& 0.16510, 0.16512, 0.16513, 0.16514, 0.16513,& 0.16512, 0.16511, 0.16510, 0.16507, 0.16508,& 0.16510, 0.16510, 0.16511, 0.16512, 0.16513,& 0.16514, 0.16514, 0.16513, 0.16511, 0.16508,& 0.16507, 0.16509, 0.16510, 0.16511, 0.16512,& 0.16513, 0.16514, 0.16514, 0.16512, 0.16510,& 0.16507, 0.16508, 0.16508, 0.16511, 0.16512,& 0.16514, 0.16514, 0.16514, 0.16514, 0.16514,& 0.16513, 0.16511, 0.16511, 0.16510, 0.16510,& 0.16511, 0.16513, 0.16514, 0.16514, 0.16515,& 0.16514, 0.16512, 0.16510, 0.16508, 0.16508,& 0.16509, 0.16511, 0.16513, 0.16513, 0.16514,& 0.16514, 0.16515, 0.16515, 0.16512, 0.16510,& 0.16511, 0.16509, 0.16510, 0.16510, 0.16511,& 0.16513, 0.16514, 0.16514, 0.16514, 0.16513,& 0.16510, 0.16508, 0.16506, 0.16507, 0.16509,& 0.16511, 0.16512, 0.16513, 0.16513, 0.16512,& 0.16510, 0.16509, 0.16507, 0.16508, 0.16508,& 0.16509, 0.16511, 0.16512, 0.16513, 0.16512,& 0.16510, 0.16505, 0.16505, 0.16505, 0.16506,& 0.16507, 0.16509, 0.16510, 0.16511, 0.16511,& 0.16509, 0.16506, 0.16503, 0.16502, 0.16503,& 0.16506, 0.16509, 0.16510, 0.16511, 0.16511,& 0.16509, 0.16502, 0.16498, 0.16499, 0.16501,& 0.16503, 0.16508, 0.16510, 0.16510, 0.16511,& 0.16510, 0.16508, 0.16505, 0.16505, 0.16504,& 0.16504, 0.16507, 0.16506, 0.16509, 0.16510,& 0.16512, 0.16511, 0.16509, 0.16504, 0.16500,& 0.16500, 0.16499, 0.16503, 0.16505, 0.16509,& 0.16511, 0.16511, 0.16510, 0.16506, 0.16500,& 0.16501, 0.16501, 0.16504, 0.16507, 0.16509,& 0.16510, 0.16511, 0.16510, 0.16506, 0.16503,& 0.16500, 0.16501, 0.16500, 0.16505, 0.16507,& 0.16509, 0.16510, 0.16510, 0.16511 /) real, parameter :: CHOUBAND8(firstYear:finalYear) = (/ & 0.04348, 0.04348, 0.04347, 0.04347, 0.04347,& 0.04347, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04349, 0.04349, 0.04349, 0.04349, 0.04349,& 0.04349, 0.04349, 0.04349, 0.04349, 0.04349,& 0.04349, 0.04349, 0.04349, 0.04349, 0.04349,& 0.04349, 0.04349, 0.04349, 0.04349, 0.04349,& 0.04349, 0.04349, 0.04349, 0.04349, 0.04349,& 0.04349, 0.04349, 0.04349, 0.04349, 0.04349,& 0.04349, 0.04349, 0.04349, 0.04349, 0.04349,& 0.04349, 0.04349, 0.04349, 0.04349, 0.04349,& 0.04349, 0.04349, 0.04349, 0.04349, 0.04349,& 0.04349, 0.04349, 0.04349, 0.04349, 0.04349,& 0.04349, 0.04349, 0.04349, 0.04349, 0.04349,& 0.04349, 0.04349, 0.04349, 0.04349, 0.04349,& 0.04349, 0.04349, 0.04349, 0.04349, 0.04349,& 0.04349, 0.04349, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04347, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04347, 0.04347,& 0.04347, 0.04347, 0.04347, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04347, 0.04347,& 0.04347, 0.04347, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04347, 0.04347, 0.04347, 0.04347,& 0.04347, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04347, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04347, 0.04347, 0.04347, 0.04347,& 0.04347, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04347, 0.04347, 0.04347,& 0.04347, 0.04347, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04347,& 0.04347, 0.04347, 0.04347, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04347,& 0.04347, 0.04347, 0.04347, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04347, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04347, 0.04347,& 0.04347, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04347,& 0.04348, 0.04347, 0.04348, 0.04347, 0.04348,& 0.04348, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04347, 0.04347, 0.04347, 0.04347, 0.04347,& 0.04347, 0.04348, 0.04348, 0.04348, 0.04348,& 0.04347, 0.04347, 0.04347, 0.04347, 0.04347,& 0.04347, 0.04347, 0.04348, 0.04348, 0.04348,& 0.04347, 0.04347, 0.04347, 0.04347, 0.04347,& 0.04347, 0.04347, 0.04347, 0.04348, 0.04347,& 0.04347, 0.04347, 0.04347, 0.04346, 0.04346,& 0.04347, 0.04347, 0.04347, 0.04347, 0.04347,& 0.04347, 0.04346, 0.04346, 0.04346, 0.04346,& 0.04346, 0.04347, 0.04347, 0.04347, 0.04347,& 0.04347, 0.04347, 0.04347, 0.04347, 0.04346,& 0.04346, 0.04347, 0.04347, 0.04347, 0.04347,& 0.04348, 0.04347, 0.04347, 0.04346, 0.04346,& 0.04346, 0.04346, 0.04346, 0.04346, 0.04347,& 0.04347, 0.04347, 0.04347, 0.04347, 0.04346,& 0.04346, 0.04346, 0.04346, 0.04347, 0.04347,& 0.04347, 0.04347, 0.04347, 0.04347, 0.04346,& 0.04346, 0.04346, 0.04346, 0.04346, 0.04347,& 0.04347, 0.04347, 0.04347, 0.04347 /) ! Establish location in table for current, previous and next year. ! ---------------------------------------------------------------- if(Year>finalYear .or. (Year==finalYear .and. dayOfYear > 182) ) then Current = (finalYear - SolCycle) + mod(Year - finalYear, SolCycle) else Current = Year endif ! Divide the year into halves. ! --------------------------- if(dayOfYear <= 182 .or. Current<firstYear) then i1 = max(Current-1,firstYear) i2 = max(Current ,firstYear) F = (dayOfYear+183)/365.00 else i1 = Current i2 = Current+1 F = (dayOfYear-183)/365.00 end if ! Linear interpolation to the given day-of-year. ! ---------------------------------------------- SC = so (i1)*(1.-F) + so (i2)*F if(present(HK)) then if(size(HK)==8) then ! Chou Bands HK(1) = ChouBand1(i1)*(1.-F) + ChouBand1(i2)*F HK(2) = ChouBand2(i1)*(1.-F) + ChouBand2(i2)*F HK(3) = ChouBand3(i1)*(1.-F) + ChouBand3(i2)*F HK(4) = ChouBand4(i1)*(1.-F) + ChouBand4(i2)*F HK(5) = ChouBand5(i1)*(1.-F) + ChouBand5(i2)*F HK(6) = ChouBand6(i1)*(1.-F) + ChouBand6(i2)*F HK(7) = ChouBand7(i1)*(1.-F) + ChouBand7(i2)*F HK(8) = ChouBand8(i1)*(1.-F) + ChouBand8(i2)*F _ASSERT(abs(1.0-sum(HK))<1.e-4,'Chou Solar band weightings do not sum to unity!') else _FAIL('HK: Solar band weightings only available for Chou') endif end if _RETURN(ESMF_SUCCESS) end subroutine MAPL_SunGetSolarConstantByYearDoY !------------------------------------------------------------------------- !NASA/GSFC, Global Modeling and Assimilation Office, Code 610.1, GEOS/DAS! !------------------------------------------------------------------------- !> ! Acquire the solar constant and the eight band fractions for the Chou solar code ! from an external NetCDF file. The initial configuration assumes monthly values ! beginning in January, and does not interpolate beyond the range of available data. ! !@bug ! Band values for RRTMG not implemented !@endbug ! subroutine MAPL_SunGetSolarConstantFromNetcdfFile(CLOCK,fileName,SC,HK,MESOPHOT,JCALC4,rc) implicit none type(ESMF_Clock), intent(in) :: CLOCK character(len=*), intent(in) :: fileName real, intent(out) :: SC real, optional, intent(out) :: HK(:) real, optional, intent(out) :: MESOPHOT(:) real, optional, intent(out) :: JCALC4(:) integer, optional, intent(out) :: rc type(ESMF_Time) :: time integer :: N integer :: begYear, endYear integer :: INDX1, INDX2 integer :: MM, YY, DD, CCYY integer :: STATUS real :: FAC integer :: ncid integer :: ndate, dimid_ndate integer :: nbin_sorad, dimid_nbin_sorad integer :: nbin_meso_phot, dimid_nbin_meso_phot integer :: nbin_jcalc4, dimid_nbin_jcalc4 integer, dimension(: ), allocatable :: date_year integer :: varid_date_year integer, dimension(: ), allocatable :: date_month integer :: varid_date_month real, dimension(: ), allocatable :: tsi integer :: varid_tsi real, dimension(:,:), allocatable :: coef_sorad integer :: varid_coef_sorad real, dimension(:,:), allocatable :: coef_meso_phot integer :: varid_coef_meso_phot real, dimension(:,:), allocatable :: coef_jcalc4 integer :: varid_coef_jcalc4 character(len=ESMF_MAXSTR) :: IAm = "MAPL_SunGetSolarConstantFromNetcdfFile" ! Open the file ! ------------- STATUS = nf90_open(trim(fileName), NF90_NOWRITE, ncid) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error opening file ', trim(fileName), status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if ! Read in dimensions ! ------------------ status = nf90_inq_dimid(ncid, 'ndate', dimid_ndate) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting ndate dimid', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if status = nf90_inquire_dimension(ncid, dimid_ndate, len = ndate) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting ndate length', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if if (present(HK)) then status = nf90_inq_dimid(ncid, 'nbin_sorad', dimid_nbin_sorad) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting nbin_sorad dimid', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if status = nf90_inquire_dimension(ncid, dimid_nbin_sorad, len = nbin_sorad) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting nbin_sorad length', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if end if if (present(MESOPHOT)) then status = nf90_inq_dimid(ncid, 'nbin_meso_phot', dimid_nbin_meso_phot) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting nbin_meso_phot dimid', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if status = nf90_inquire_dimension(ncid, dimid_nbin_meso_phot, len = nbin_meso_phot) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting nbin_meso_phot length', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if end if if (present(JCALC4)) then status = nf90_inq_dimid(ncid, 'nbin_jcalc4', dimid_nbin_jcalc4) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting nbin_jcalc4 dimid', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if status = nf90_inquire_dimension(ncid, dimid_nbin_jcalc4, len = nbin_jcalc4) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting nbin_jcalc4 length', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if end if ! Allocate our arrays ! ------------------- allocate(date_year(ndate), source=0, stat=status) _VERIFY(STATUS) allocate(date_month(ndate), source=0, stat=status) _VERIFY(STATUS) allocate(tsi(ndate), source=0.0, stat=status) _VERIFY(STATUS) if (present(HK) ) then allocate(coef_sorad(nbin_sorad,ndate), source=0.0, stat=status) _VERIFY(STATUS) end if if (present(MESOPHOT) ) then allocate(coef_meso_phot(nbin_meso_phot,ndate), source=0.0, stat=status) _VERIFY(STATUS) end if if (present(JCALC4) ) then allocate(coef_jcalc4(nbin_jcalc4,ndate), source=0.0, stat=status) _VERIFY(STATUS) end if ! Read in date_year ! ----------------- status = nf90_inq_varid(ncid, 'date_year', varid_date_year) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting date_year varid', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if status = nf90_get_var(ncid, varid_date_year, date_year) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting date_year variable', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if ! Read in date_month ! ------------------ status = nf90_inq_varid(ncid, 'date_month', varid_date_month) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting date_month varid', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if status = nf90_get_var(ncid, varid_date_month, date_month) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting date_month variable', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if ! Read in tsi ! ----------- status = nf90_inq_varid(ncid, 'tsi', varid_tsi) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting tsi varid', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if status = nf90_get_var(ncid, varid_tsi, tsi) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting tsi variable', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if ! Read in coef_sorad ! ------------------ if (present(HK)) then status = nf90_inq_varid(ncid, 'coef_sorad', varid_coef_sorad) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting coef_sorad varid', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if status = nf90_get_var(ncid, varid_coef_sorad, coef_sorad) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting coef_sorad variable', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if end if ! Read in coef_meso_phot ! ---------------------- if (present(MESOPHOT)) then status = nf90_inq_varid(ncid, 'coef_meso_phot', varid_coef_meso_phot) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting coef_meso_phot varid', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if status = nf90_get_var(ncid, varid_coef_meso_phot, coef_meso_phot) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting coef_meso_phot variable', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if end if ! Read in coef_jcalc4 ! ------------------- if (present(JCALC4)) then status = nf90_inq_varid(ncid, 'coef_jcalc4', varid_coef_jcalc4) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting coef_jcalc4 varid', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if status = nf90_get_var(ncid, varid_coef_jcalc4, coef_jcalc4) if (STATUS /= NF90_NOERR) then write (*,*) trim(Iam)//': Error getting coef_jcalc4 variable', status write (*,*) nf90_strerror(status) _FAIL('needs informative message') end if end if ! Time interpolation parameters ! ----------------------------- call MAPL_ClimInterpFac(CLOCK, INDX1, INDX2, FAC, RC=STATUS) _VERIFY(STATUS) ! Time: Year and month ! -------------------- call ESMF_ClockGet(CLOCK, CURRTIME=time, RC=STATUS) _VERIFY(STATUS) call ESMF_TimeGet(time, YY=YY, MM=MM, DD=DD, RC=STATUS) _VERIFY(STATUS) ! This index search routine below is adaptedfrom PChemGridComp ! ! Read bracketing months, make sure INDX1 and INDX2 are in range. Annual ! cycle is preserved for years that precede and succeed the climatology. ! ---------------------------------------------------------------------- N = 12*ndate begYear = date_year(1) endYear = date_year(ndate) CCYY = YY IF(CCYY < begYear) CCYY = begYear IF(CCYY > endYear) CCYY = endYear INDX1 = INDX1+(CCYY-begYear)*12 INDX2 = INDX2+(CCYY-begYear)*12 IF(MM == 1 .AND. INDX1 > INDX2) INDX1 = INDX1-12 IF(MM == 12 .AND. INDX2 < INDX1) INDX2 = INDX2+12 IF(INDX1 == 0) INDX1 = 12 IF(YY < begYear .AND. INDX2 > 12) INDX2 = 1 IF(INDX2 == N+1) INDX2 = N-11 IF(YY > endYear .AND. INDX1 == N-12) INDX1 = N ! Linear Interpolation to the given day-of-month ! ---------------------------------------------- SC = tsi(INDX1)*FAC + tsi(INDX2)*(1.0-FAC) if (present(HK)) then HK = coef_sorad(:,INDX1)*FAC + coef_sorad(:,INDX2)*(1.0-FAC) end if if (present(MESOPHOT)) then HK = coef_meso_phot(:,INDX1)*FAC + coef_meso_phot(:,INDX2)*(1.0-FAC) end if if (present(JCALC4)) then HK = coef_jcalc4(:,INDX1)*FAC + coef_jcalc4(:,INDX2)*(1.0-FAC) end if ! Close the file ! -------------- STATUS = nf90_close(ncid) _VERIFY(STATUS) ! Bounds check ! ------------ if (PRESENT(HK)) then _ASSERT(ABS(1.0-SUM(HK)) < 1.E-4,'needs informative message') end if ! Deallocate our arrays ! --------------------- deallocate(date_year, stat=status) _VERIFY(STATUS) deallocate(date_month, stat=status) _VERIFY(STATUS) deallocate(tsi, stat=status) _VERIFY(STATUS) if (present(HK)) then deallocate(coef_sorad, stat=status) _VERIFY(STATUS) end if if (present(MESOPHOT)) then deallocate(coef_meso_phot, stat=status) _VERIFY(STATUS) end if if (present(JCALC4)) then deallocate(coef_jcalc4, stat=status) _VERIFY(STATUS) end if _RETURN(ESMF_SUCCESS) end subroutine MAPL_SunGetSolarConstantFromNetcdfFile !------------------------------------------------------------------------- !NASA/GSFC, Global Modeling and Assimilation Office, Code 610.1, GEOS/DAS! !------------------------------------------------------------------------- !> ! Acquire the TSI, Mg, and SB from file. ! ! @bug ! Band values for RRTMG not implemented !@endbug ! subroutine MAPL_SunGetSolarConstantFromNRLFile(CLOCK,filename_in,SC,MG,SB,PersistSolar,rc) implicit none type(ESMF_Clock), intent(in) :: CLOCK character(len=*), intent(in) :: filename_in real, intent(out) :: SC real, intent(out) :: MG real, intent(out) :: SB logical, optional, intent(in) :: PersistSolar integer, optional, intent(out) :: rc type(ESMF_VM) :: VM type(ESMF_Time) :: noonCurrentDay, prevNoon, nextNoon type(ESMF_TimeInterval) :: intToNextNoon, oneDayInterval type(ESMF_Time) :: currentTime, origTime type(ESMF_Time) :: timeBasedOnCycle24 type(ESMF_Time) :: startCycle24, startCycle25 type(ESMF_TimeInterval) :: timeSinceStartOfCycle25 integer :: currentYear, currentMon, currentDay, currentDOY integer :: prevDOY, nextDOY, prevNoonYear, nextNoonYear integer :: originalYear, originalMon, originalDay, origDOY real(ESMF_KIND_R8) :: days_r8 integer :: INDX1, INDX2, i integer :: STATUS, UNIT real :: FAC character(len=ESMF_MAXPATHLEN) :: FILENAME logical :: amIRoot integer :: deId, NPES logical :: outOfTable logical :: PersistSolar_ integer, parameter :: YEAR_NOT_FOUND = -99 character(len=256) :: line character(len=1), parameter :: commentChar = "#" integer :: stat logical, save :: TableCreated = .FALSE. integer, save, allocatable, dimension(:) :: yearTable, doyTable real, save, allocatable, dimension(:) :: tsi, mgindex, sbindex integer, save :: numlines class(Logger), pointer :: lgr lgr => logging%get_logger('MAPL.SUN') if (present(PersistSolar)) then PersistSolar_ = PersistSolar else PersistSolar_ = .TRUE. end if call ESMF_VMGetCurrent(vm, rc=status) call ESMF_VmGet(VM, localPet=deId, petCount=npes, rc=status) _VERIFY(STATUS) amIRoot = (deId == 0) CREATE_TABLE: if (.not. TableCreated) then ! First we open the file on root to get the ! number of lines so we can allocate our tables ! --------------------------------------------- if (amIRoot) then ! Open the file ! ------------- filename = trim(filename_in) open(newunit=unit, file=filename, form="formatted", status="old", iostat=status) _ASSERT(status==0,'Could not find NRL data file '// trim(filename )) ! Determine length of file ! ------------------------ call lgr%debug("Scanning the Solar Table to determine number of data points") numlines = num_lines_in_file(UNIT) call lgr%debug("Solar Table Data Points: %i0", numlines) end if ! Broadcast the number of lines ! ----------------------------- call MAPL_CommsBcast(vm, DATA=numlines, N=1, ROOT=0, _RC) ! Allocate our arrays on all processes ! ------------------------------------ allocate(yearTable(numlines), source=0, _STAT) allocate(doyTable(numlines), source=0, _STAT) allocate(tsi(numlines), source=0.0, _STAT) allocate(mgindex(numlines), source=0.0, _STAT) allocate(sbindex(numlines), source=0.0, _STAT) ! Back to root to read in the values ! ---------------------------------- if (amIRoot) then call lgr%debug("Reading the Solar Table") i = 1 do read(unit,'(A)',iostat=stat) line if (is_iostat_end(stat)) exit ! Skip lines starting with #-comments if (index (line, commentChar) /= 0) cycle read(line,*) yearTable(i), doyTable(i), tsi(i), mgindex(i), sbindex(i) i = i + 1 end do ! Belt and suspenders check that all data was read _ASSERT(size(yearTable) == numlines,"Inconsistency in NRL number of lines") close(unit, _IOSTAT) end if ! Broadcast the tables ! -------------------- call MAPL_CommsBcast(vm, DATA=yearTable, N=numlines, ROOT=0, _RC) call MAPL_CommsBcast(vm, DATA=doyTable, N=numlines, ROOT=0, _RC) call MAPL_CommsBcast(vm, DATA=tsi, N=numlines, ROOT=0, _RC) call MAPL_CommsBcast(vm, DATA=mgindex, N=numlines, ROOT=0, _RC) call MAPL_CommsBcast(vm, DATA=sbindex, N=numlines, ROOT=0, _RC) TableCreated = .TRUE. end if CREATE_TABLE ! Now we need to find the two bracketing days ! ------------------------------------------- ! Get current time ! ---------------- call ESMF_ClockGet(CLOCK, CURRTIME=currentTime, _RC) call ESMF_TimeGet( currentTime, YY = currentYear, & MM = currentMon, & DD = currentDay, & dayOfYear = currentDOY, _RC) ! Test if current time is outside our file ! ---------------------------------------- outOfTable = .FALSE. ! First is current year higher than last in file... if ( currentYear > yearTable(numlines) ) then outOfTable = .TRUE. ! ...or if a partial year, are we near the end else if ( currentYear == yearTable(numlines) .and. currentDOY >= doyTable(numlines)) then outOfTable = .TRUE. end if ! If we are out of the table... ! ----------------------------- OUT_OF_TABLE: if ( outOfTable ) then PERSIST_SOLAR: if ( PersistSolar_ ) then ! If we are outOfTable and we have the PersistSolar ! option we just use the last value in the table... ! ------------------------------------------------- SC = tsi(numlines) MG = mgindex(numlines) SB = sbindex(numlines) call lgr%debug('Off the end of table, persisting last values') call lgr%debug(' tsi at end of table: %F8.3', tsi(numlines)) call lgr%debug(' mgindex at end of table: %F8.6', mgindex(numlines)) call lgr%debug(' sbindex at end of table: %F9.4', sbindex(numlines)) _RETURN(ESMF_SUCCESS) else ! If we are out of the table and not persisting, we must ! recenter our day to be based on the last complete Solar Cycle ! ------------------------------------------------------------- ! Create an ESMF_Time at start of Cycle 24 ! ---------------------------------------- call ESMF_TimeSet( startCycle24, YY = 2008, & MM = 12, & DD = 1, & H = 12, & M = 00, & S = 00, _RC) ! Create an ESMF_Time at start of Cycle 25 ! ---------------------------------------- call ESMF_TimeSet( startCycle25, YY = 2019, & MM = 12, & DD = 1, & H = 12, & M = 00, & S = 00, _RC) ! Create TimeInterval based on interval ! from start of latest Cycle 25 ! ------------------------------------- timeSinceStartOfCycle25 = currentTime - startCycle25 ! Make a new time based on that ! interval past start of Cycle 24 ! ------------------------------- timeBasedOnCycle24 = startCycle24 + timeSinceStartOfCycle25 ! Store our original time just in case ! ------------------------------------ origTime = currentTime originalYear = currentYear originalMon = currentMon originalDay = currentDay origDOY = currentDOY ! Make our "current" time the one calculated above ! ------------------------------------------------ currentTime = timeBasedOnCycle24 ! Get new currentYear, currentMon, currentDay ! ------------------------------------------- call ESMF_TimeGet( currentTime, YY = currentYear, & MM = currentMon, & DD = currentDay, & dayOfYear = currentDOY, _RC) call lgr%debug('Off the end of table, moving into last complete cycle') call lgr%debug(' Original Year-Mon-Day to Find: %i0.4~-%i0.2~-%i0.2', originalYear,originalMon,originalDay) call lgr%debug(' Original Day of Year: %i0', origDOY) call lgr%debug(' New Year-Mon-Day to Find: %i0.4~-%i0.2~-%i0.2', currentYear,currentMon,currentDay) call lgr%debug(' New Day of Year: %i0', currentDOY) end if PERSIST_SOLAR end if OUT_OF_TABLE ! Create an ESMF_Time on noon of current day ! ------------------------------------------ call ESMF_TimeSet( noonCurrentDay, YY = currentYear, & MM = currentMon, & DD = currentDay, & H = 12, & M = 00, & S = 00, _RC) ! Figure out bracketing days for interpolation ! NOTE: nextNoon is mainly for debugging purposes ! ----------------------------------------------- call ESMF_TimeIntervalSet(oneDayInterval, D=1, _RC) if (currentTime <= noonCurrentDay) then prevNoon = noonCurrentDay - oneDayInterval nextNoon = noonCurrentDay else prevNoon = noonCurrentDay nextNoon = noonCurrentDay + oneDayInterval end if ! Get the DOYs ! ------------ call ESMF_TimeGet( prevNoon, YY = prevNoonYear, dayOfYear = prevDOY, _RC) call ESMF_TimeGet( nextNoon, YY = nextNoonYear, dayOfYear = nextDOY, _RC) ! Our interpolation factor is based of when we are compared to the next noon ! -------------------------------------------------------------------------- intToNextNoon = nextNoon-currentTime ! The FAC for interpolating is just the real version ! of the size of the timeinterval to the next noon ! -------------------------------------------------- call ESMF_TimeIntervalGet(intToNextNoon, d_r8=days_r8, _RC) FAC = real(days_r8) ! Use our find_file_index function to get the index for previous noon ! ------------------------------------------------------------------- INDX1 = find_file_index(numlines, yearTable, prevNoonYear, prevDOY) INDX2 = INDX1 + 1 ! Linear Interpolation to the given day-of-month ! ---------------------------------------------- SC = tsi(INDX1)*FAC + tsi(INDX2)*(1.0-FAC) MG = mgindex(INDX1)*FAC + mgindex(INDX2)*(1.0-FAC) SB = sbindex(INDX1)*FAC + sbindex(INDX2)*(1.0-FAC) call lgr%debug(' First DOY to Find: %i3', prevDOY) call lgr%debug(' file_index for date: %i6', INDX1) call lgr%debug(' yearTable(date): %i4', yearTable(INDX1)) call lgr%debug(' doyTable(date): %i3', doyTable(INDX1)) call lgr%debug(' tsi(date): %f8.3', tsi(INDX1)) call lgr%debug(' mgindex(date): %f8.6', mgindex(INDX1)) call lgr%debug(' sbindex(date): %f9.4', sbindex(INDX1)) call lgr%debug(' Second DOY to Find: %i3', nextDOY) call lgr%debug(' file_index for date: %i6', INDX2) call lgr%debug(' yearTable(date): %i4', yearTable(INDX2)) call lgr%debug(' doyTable(date): %i3', doyTable(INDX2)) call lgr%debug(' tsi(date): %f8.3', tsi(INDX2)) call lgr%debug(' mgindex(date): %f8.6', mgindex(INDX2)) call lgr%debug(' sbindex(date): %f9.4', sbindex(INDX2)) call lgr%debug(' Interpolation Factor: %f8.6', FAC) _RETURN(ESMF_SUCCESS) contains integer function num_lines_in_file(unit) result(count) implicit none integer, intent(in) :: unit character(len=256) :: line integer :: stat character(len=1), parameter :: commentChar = "#" count = 0 do read(unit,'(A)',iostat=stat) line if (is_iostat_end(stat)) exit ! Skip lines starting with #-comments if (index (line, commentChar) /= 0) cycle count = count + 1 end do rewind(unit) end function num_lines_in_file integer function find_file_index(numlines, year, year_to_find, doy_to_find) implicit none integer, intent(in) :: numlines integer, intent(in), dimension(numlines) :: year integer, intent(in) :: year_to_find, doy_to_find integer :: i, yr_index ! Find the index for the year ! --------------------------- yr_index = YEAR_NOT_FOUND do i = 1, numlines if (year(i) == year_to_find) then yr_index = i exit end if end do if (yr_index == YEAR_NOT_FOUND) then find_file_index = YEAR_NOT_FOUND else find_file_index = yr_index + doy_to_find - 1 end if end function find_file_index end subroutine MAPL_SunGetSolarConstantFromNRLFile !========================================================================== !> ! Return the daylight duration in seconds (i.e, the time between sunrise and sunset) for ! a day around the specified time. The routine is accurate enough for most purposes, but ! does not solve for precise sunrise/sunset times influenced by changes in the orbital ! parameters between those times. The time input does NOT need to be noon --- it is used ! simply to evaluate the solar declination needed for the daylight duration calculation. ! In order of preference, time is taken from currTime, if present, or else the currTime ! of CLOCK, if present, or else the currTime of the ORBIT's associated clock. ! subroutine MAPL_SunGetDaylightDuration(ORBIT,LATS,DAYL,currTime,CLOCK,RC) ! !ARGUMENTS: type(MAPL_SunOrbit), intent(IN ) :: ORBIT real, dimension(:) , intent(IN ) :: LATS real, dimension(:) , intent(OUT) :: DAYL type(ESMF_Time) , optional, intent(IN ) :: currTime type(ESMF_Clock) , optional, intent(IN ) :: CLOCK integer, optional, intent(OUT) :: RC ! Locals character(len=ESMF_MAXSTR), parameter :: IAm = "MAPL_SunGetDaylightDuration" integer :: STATUS type(ESMF_Time) :: CURRENTTIME integer :: YEAR, SEC_OF_DAY, DAY_OF_YEAR, IDAY, IDAYP1 real :: FAC, ZS, ZC real(ESMF_KIND_R8) :: days real :: ECC, OBQ, LAMBDAP real :: OMECC, OPECC, OMSQECC, EAFAC real :: MA, EA, dE, TA, LAMBDA integer :: nits ! Begin _ASSERT(MAPL_SunOrbitCreated(ORBIT),'MAPL_SunOrbit not yet created!') ! which current time? if (present(currTime)) then CURRENTTIME = CURRTIME else if (present(CLOCK)) then call ESMF_ClockGet( CLOCK, currTime=CURRENTTIME, RC=STATUS) else call ESMF_ClockGet(ORBIT%CLOCK, currTime=CURRENTTIME, RC=STATUS) end if _VERIFY(STATUS) end if if (ORBIT%ANAL2B) then ! include time variation in orbit from reference time call ESMF_TimeIntervalGet( & CURRENTTIME - ORBIT%ORB2B_TIME_REF, & d_r8=days, rc=STATUS) _VERIFY(STATUS) ECC = ORBIT%ORB2B_ECC_REF + days * ORBIT%ORB2B_ECC_RATE OBQ = ORBIT%ORB2B_OBQ_REF + days * ORBIT%ORB2B_OBQ_RATE LAMBDAP = ORBIT%ORB2B_LAMBDAP_REF + days * ORBIT%ORB2B_LAMBDAP_RATE ! derived quantities OMECC = 1. - ECC OPECC = 1. + ECC OMSQECC = OMECC * OPECC EAFAC = sqrt(OMECC/OPECC) ! time interval since perhelion in days call ESMF_TimeIntervalGet( & CURRENTTIME - ORBIT%ORB2B_TIME_PERI, & d_r8=days, rc=STATUS) _VERIFY(STATUS) ! mean anomaly MA = ORBIT%ORB2B_OMG0 * days ! eccentric anomaly call invert_Keplers_Newton(MA,ECC,EA,dE,nits) ! true anomaly TA = calcTAfromEA(EA,EAFAC) ! celestial longitude LAMBDA = TA + LAMBDAP ! sin and cos of solar declination ZS = sin(LAMBDA) * sin(OBQ) ZC = sqrt(1. - ZS**2) else call ESMF_TimeGet(CURRENTTIME, YY=YEAR, S=SEC_OF_DAY, & dayOfYear=DAY_OF_YEAR, RC=STATUS) _VERIFY(STATUS) YEAR = mod(YEAR-1,ORBIT%YEARS_PER_CYCLE) IDAY = YEAR*int(ORBIT%YEARLEN)+DAY_OF_YEAR IDAYP1 = mod(IDAY,ORBIT%DAYS_PER_CYCLE) + 1 FAC = real(SEC_OF_DAY)/86400. ZS = ORBIT%ZS(IDAYP1)*FAC + ORBIT%ZS(IDAY)*(1.-FAC) ZC = ORBIT%ZC(IDAYP1)*FAC + ORBIT%ZC(IDAY)*(1.-FAC) endif ! dayligt duration [secs] DAYL = (86400./MAPL_PI)*acos(min(1.,max(-1.,-tan(LATS)*ZS/ZC))) _RETURN(ESMF_SUCCESS) end subroutine MAPL_SunGetDaylightDuration !========================================================================== !> ! Return the daylight duration in seconds (i.e, the time between sunrise and sunset) for ! its MAXIMUM at the summer solstice. The routine is accurate enough for most purposes, ! but does not solve for precise sunrise/sunset times influenced by changes in the orbital ! parameters between those times. The time input does NOT need to be noon --- it is used ! simply to evaluate the obliquity needed for the maximum daylight duration calculation. ! In order of preference, time is taken from currTime, if present, or else the currTime ! of CLOCK, if present, or else the currTime of the ORBIT's associated clock. ! Note: Unless ORBIT_ANAL2B, the obliquity is fixed and the time is irrelevant. ! subroutine MAPL_SunGetDaylightDurationMax(ORBIT,LATS,DAYL,currTime,CLOCK,RC) ! !ARGUMENTS: type(MAPL_SunOrbit), intent(IN ) :: ORBIT real, dimension(:) , intent(IN ) :: LATS real, dimension(:) , intent(OUT) :: DAYL type(ESMF_Time) , optional, intent(IN ) :: currTime type(ESMF_Clock) , optional, intent(IN ) :: CLOCK integer, optional, intent(OUT) :: RC !EOPI character(len=ESMF_MAXSTR), parameter :: IAm = "MAPL_SunGetDaylightDurationMax" integer :: STATUS type(ESMF_Time) :: CURRENTTIME real(ESMF_KIND_R8) :: days real :: OBQ _ASSERT(MAPL_SunOrbitCreated(ORBIT),'MAPL_SunOrbit not yet created!') ! Which time? if (present(currTime)) then CURRENTTIME = CURRTIME else if (present(CLOCK)) then call ESMF_ClockGet( CLOCK, currTime=CURRENTTIME, RC=STATUS) else call ESMF_ClockGet(ORBIT%CLOCK, currTime=CURRENTTIME, RC=STATUS) end if _VERIFY(STATUS) end if if (ORBIT%ANAL2B) then ! time variation in obliquity from ref time call ESMF_TimeIntervalGet( & CURRENTTIME - ORBIT%ORB2B_TIME_REF, & d_r8=days, rc=STATUS) _VERIFY(STATUS) OBQ = ORBIT%ORB2B_OBQ_REF + days * ORBIT%ORB2B_OBQ_RATE else ! obliquity fixed in this case OBQ = ORBIT%OB * (MAPL_PI/180.) endif _ASSERT(OBQ >= 0, 'Obliquity less than 0 detected!') _ASSERT(OBQ < MAPL_PI, 'Obliquity greater than pi detected!') ! Maximum daylight duration at summer solstice [secs] ! (an even function of latitude) DAYL = (86400./MAPL_PI)*acos(min(1.,max(-1., & -tan(ABS(LATS))*tan(OBQ)))) _RETURN(ESMF_SUCCESS) end subroutine MAPL_SunGetDaylightDurationMax !========================================================================== !> ! Returns the local solar hour angle (in radians) at the single time and ! multiple longitudes specified. In order of preference, time is taken ! from TIME, if present, or else the CURRTIME of CLOCK, if present, or ! else the CURRTIME of the ORBIT's associated clock. ! ! NB: For accurate results, namely to receive the TRUE local solar hour ! angle, ensure the ORBIT has the EOT flag set true. Conversely, to get ! only the MEAN local solar hour angle, use the optional argument ! FORCE_MLSHA=.TRUE.. This will turn off the Equation of Time correction ! (for this LSHA calculation only) even if the ORBIT includes it. For ! example, in the local noon detection in the EXAMPLE below, this will ! give mean local noons that are exactly 24 hours apart at a particular ! location. But they will no longer exactly be the solar culmination ! times (the TRUE local noon) in that case. TRUE local noons are not ! exactly 24h apart because of orbital variations in length of day ! throughout the year, as described by the Equation of Time. ! ! @note ! `Example of use:`: ! ! type (ESMF_Time) :: NOW ! type (ESMF_TimeInterval) :: DELT ! real, dimension(size(LONS)) :: LSHA0, LSHA1 ! logical, dimension(size(LONS)) :: isNoon ! call ESMF_ClockGet (CLOCK, CURRTIME=NOW, TIMESTEP=DELT, __RC__) ! call MAPL_SunGetLocalSolarHourAngle (ORBIT, LONS, LSHA0, TIME=NOW, __RC__) ! call MAPL_SunGetLocalSolarHourAngle (ORBIT, LONS, LSHA1, TIME=NOW+DELT, __RC__) ! isnoon = (LSHA0 <= 0. .and. LSHA1 > 0.) !@endnote ! subroutine MAPL_SunGetLocalSolarHourAngle (ORBIT,LONS,LSHA, & TIME,CLOCK,FORCE_MLSHA,RC) ! !ARGUMENTS: type (MAPL_SunOrbit), intent(IN ) :: ORBIT real, dimension(:), intent(IN ) :: LONS ! [radians] real, dimension(:), intent(OUT) :: LSHA type (ESMF_Time), optional, intent(IN ) :: TIME type (ESMF_Clock), optional, intent(IN ) :: CLOCK logical, optional, intent(IN ) :: FORCE_MLSHA integer, optional, intent(OUT) :: RC !EOPI ! Locals character(len=ESMF_MAXSTR), parameter :: IAm = "MAPL_SunGetLocalSolarHourAngle" integer :: STATUS type (ESMF_Time) :: T real (ESMF_KIND_R8) :: days integer :: YEAR, SEC_OF_DAY, DAY_OF_YEAR, IDAY, IDAYP1 real :: DFRAC, GSHA real :: ECC, OBQ, LAMBDAP real :: OMECC, OPECC, OMSQECC, EAFAC real :: MA, EA, dE, TA, LAMBDA real :: RT, RM, ET integer :: i, nits logical :: do_EOT _ASSERT (MAPL_SunOrbitCreated(ORBIT),'MAPL_SunOrbit not yet created!') ! Which time? if (present(TIME)) then T = TIME else if (present(CLOCK)) then call ESMF_ClockGet ( CLOCK, CURRTIME=T, _RC) else call ESMF_ClockGet (ORBIT%CLOCK, CURRTIME=T, _RC) end if end if ! NB: include YY and dayOfYear here so that S is seconds WITHIN a day. ! YEAR and DAY_OF_YEAR are used within the non-ANAL2B branch anyway. call ESMF_TimeGet (T, YY=YEAR, dayOfYear=DAY_OF_YEAR, S=SEC_OF_DAY, RC=STATUS) _VERIFY(STATUS) ! fraction of day (0 at midnight) DFRAC = real(SEC_OF_DAY) / 86400. ! Greenwich MEAN solar hour angle (zero at noon) GSHA = 2. * MAPL_PI * (DFRAC - 0.5) ! Apply equation of time correction? do_EOT = ORBIT%EOT if (present(FORCE_MLSHA)) then if (FORCE_MLSHA) do_EOT = .FALSE. endif if (do_EOT) then if (ORBIT%ANAL2B) then ! include time variation in orbit from reference time call ESMF_TimeIntervalGet (T - ORBIT%ORB2B_TIME_REF, d_r8=days, _RC) ECC = ORBIT%ORB2B_ECC_REF + days * ORBIT%ORB2B_ECC_RATE OBQ = ORBIT%ORB2B_OBQ_REF + days * ORBIT%ORB2B_OBQ_RATE LAMBDAP = ORBIT%ORB2B_LAMBDAP_REF + days * ORBIT%ORB2B_LAMBDAP_RATE ! derived quantities OMECC = 1. - ECC OPECC = 1. + ECC OMSQECC = OMECC * OPECC EAFAC = sqrt(OMECC/OPECC) ! time interval since perhelion in days call ESMF_TimeIntervalGet (T - ORBIT%ORB2B_TIME_PERI, d_r8=days, _RC) ! mean anomaly MA = ORBIT%ORB2B_OMG0 * days ! eccentric anomaly call invert_Keplers_Newton (MA,ECC,EA,dE,nits) ! true anomaly TA = calcTAfromEA (EA,EAFAC) ! celestial longitude LAMBDA = TA + LAMBDAP ! solar right ascension (true and mean) RT = atan2(sin(LAMBDA)*cos(OBQ),cos(LAMBDA)) RM = MA + LAMBDAP ! equation of time ET = RECT_PMPI (RM - RT) else ! get equation of time by table interpolation YEAR = mod(YEAR-1,ORBIT%YEARS_PER_CYCLE) IDAY = YEAR*int(ORBIT%YEARLEN)+DAY_OF_YEAR IDAYP1 = mod(IDAY,ORBIT%DAYS_PER_CYCLE) + 1 ET = ORBIT%ET(IDAYP1)*DFRAC + ORBIT%ET(IDAY)*(1.-DFRAC) endif ! Gives Greenwich TRUE solar hour angle GSHA = GSHA + ET end if ! EOT correction ! LOCAL solar hour angle do i = 1, size(LONS) LSHA(i) = RECT_PMPI (GSHA + LONS(i)) end do _RETURN(ESMF_SUCCESS) end subroutine MAPL_SunGetLocalSolarHourAngle !========================================================================== end module MAPL_SunMod