Merge pull request #1913 from GEOS-ESM/feature/pnorris/#1882-add-sola…

…r-noon-diagnostic

Feature/pnorris/#1882 add local solar hour angle diagnostic

CHANGELOG.md

@@ -9,6 +9,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- Added subroutine `MAPL_SunGetLocalSolarHourAngle()` in `base/MAPL_sun_uc.F90`. This provides a
+  convenient local solar hour angle diagnostic which will be used to detect local
+  solar noon via the EXAMPLE OF USE in the subroutine header. See DESCRIPTION in code for more
+  details. Provides the TRUE local solar hour angle (i.e., with equation of time included), but
+  can also provide the MEAN value (without EOT) via FORCE_MLSHA=.TRUE. optional argument.
+
 ### Changed
 
 ### Fixed

base/MAPL_sun_uc.F90

@@ -41,6 +41,7 @@ module MAPL_SunMod
   public MAPL_SunGetSolarConstant
   public MAPL_SunGetDaylightDuration
   public MAPL_SunGetDaylightDurationMax
+  public MAPL_SunGetLocalSolarHourAngle
 
 ! !PUBLIC TYPES:
 
@@ -2979,4 +2980,157 @@ end subroutine MAPL_SunGetDaylightDurationMax
 
 !==========================================================================
 
+!BOPI
+
+! !IROUTINE:  MAPL_SunGetLocalSolarHourAngle
+
+! !DESCRIPTION:
+
+! 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.
+
+! !INTERFACE:
+
+   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
+
+! !EXAMPLE OF USE:
+!
+!     ! detecting noon within the current timestep
+!     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.)
+
+!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