Re-enable AnOHM with rolling buffer forcing

CHANGELOG.md

@@ -170,6 +170,12 @@
   - Same functionality, lighter dependencies
   - Consolidates all parallel processing to one library
 
+### 16 Nov 2025
+- [feature] Re-enabled AnOHM without requiring future forcing data
+  - Added rolling 24-hour forcing buffers to `OHM_STATE` and `SUEWS_cal_Qs` so StorageHeatMethod 3 now diagnoses coefficients from the most recently completed day
+  - Refactored `module_phys_anohm` to consume buffered shortwave, meteorology, and anthropogenic heat series directly (removed the legacy `MetForcingData_grid` dependency)
+  - Updated StorageHeatMethod documentation to clarify the trailing-day workflow and avoid confusion about forcing requirements
+
 ### 14 Nov 2025
 - [feature] Added `SUEWSSimulation.from_sample_data()` factory method and comprehensive OOP enhancements (#779)
   - New factory method for cleaner OOP workflow: `sim = SUEWSSimulation.from_sample_data()`

docs/source/inputs/tables/RunControl/csv-table/StorageHeatMethod.csv

@@ -1,5 +1,5 @@
 "Value","Comments"
 0,"Uses observed values of ΔQS supplied in meteorological forcing file."
 1,"ΔQS modelled using the objective hysteresis model (OHM) :cite:`G91` using parameters specified for each surface type."
-3,"ΔQS modelled using AnOHM :cite:`S17`. |NotRecmd|"
+3,"ΔQS modelled using AnOHM :cite:`S17` with coefficients diagnosed from the most recent completed day of forcing (no future data required). |NotRecmd|"
 4,"ΔQS modelled using the Element Surface Temperature Method (ESTM) :cite:`O05`. |NotRecmd|"

docs/source/parameterisations-and-sub-models.rst

@@ -99,6 +99,7 @@ Storage heat flux, ΔQ\ :sub:`S`
 
    -  **OHM** (Objective Hysteresis Model) :cite:`G91,GO99,GO02`. Storage heat heat flux is calculated using empirically-fitted relations with net all-wave radiation and the rate of change in net all-wave radiation.
    -  **AnOHM** (Analytical Objective Hysteresis Model) :cite:`S17`. OHM approach using analytically-derived coefficients. |NotRecmd|
+      -  Coefficients are now updated using the most recently completed 24 h of observed forcing, so coupling no longer requires access to future meteorological data.
    -  **ESTM** (Element Surface Temperature Method) :cite:`O05`. Heat transfer through urban facets (roof, wall, road, interior) is calculated from surface temperature measurements and knowledge of material properties. |NotRecmd|
    -  **EHC** (Explicit Heat Conduction). Calculates storage heat flux using explicit heat conduction through urban facets with separate roof/wall/ground temperature calculations. Provides detailed surface temperature outputs for each urban element.
    -  **DyOHM** (Dynamic Objective Hysteresis Model) (Liu et al., in prep.). Extends OHM by calculating coefficients dynamically based on material thermal properties (thermal conductivity) and meteorological conditions. Requires vertical wall layer configuration.

src/suews/src/suews_ctrl_driver.f95

@@ -24,14 +24,15 @@ MODULE SUEWS_Driver
                             ROUGHNESS_STATE, solar_State, atm_state, flag_STATE, &
                             SUEWS_STATE, SUEWS_DEBUG, STEBBS_STATE, BUILDING_ARCHETYPE_PRM, STEBBS_PRM, STEBBS_BLDG, &
                             SUEWS_STATE_BLOCK, &
-                            output_line, output_block
+                            output_line, output_block, &
+                            ANOHM_MAX_SAMPLES, ANOHM_MIN_VALID_SAMPLES
    USE meteo, ONLY: qsatf, RH2qa, qa2RH
    USE module_phys_atmmoiststab, ONLY: cal_AtmMoist, cal_Stab, stab_psi_heat, stab_psi_mom, SUEWS_update_atmState
    USE module_phys_narp, ONLY: NARP_cal_SunPosition, NARP_cal_SunPosition_DTS
    USE module_phys_atmmoiststab, ONLY: cal_AtmMoist, cal_Stab, stab_psi_heat, stab_psi_mom
    USE module_phys_narp, ONLY: NARP_cal_SunPosition
    USE module_phys_spartacus, ONLY: SPARTACUS
-   ! USE AnOHM_module, ONLY: AnOHM
+   USE module_phys_anohm, ONLY: AnOHM
    USE module_phys_resist, ONLY: AerodynamicResistance, BoundaryLayerResistance, SurfaceResistance, &
                             SUEWS_cal_RoughnessParameters
    USE module_phys_ohm, ONLY: OHM
@@ -1496,13 +1497,15 @@ SUBROUTINE SUEWS_cal_Qs( &
       REAL(KIND(1D0)), DIMENSION(nsurf) :: cpAnOHM ! heat capacity [J m-3 K-1]
       REAL(KIND(1D0)), DIMENSION(nsurf) :: kkAnOHM ! thermal conductivity [W m-1 K-1]
       REAL(KIND(1D0)), DIMENSION(nsurf) :: chAnOHM ! bulk transfer coef [J m-3 K-1]
+      REAL(KIND(1D0)), DIMENSION(nsurf) :: moist_surf ! normalised soil moisture state [-]
 
       REAL(KIND(1D0)), DIMENSION(27), INTENT(out) :: dataOutLineESTM !data output from ESTM
       REAL(KIND(1D0)), DIMENSION(ncolumnsDataOutSTEBBS - 5), INTENT(out) :: dataOutLineSTEBBS !data output from STEBBS
 
       ! internal use arrays
       REAL(KIND(1D0)) :: Tair_mav_5d ! Tair_mav_5d=HDD(id-1,4) HDD at the begining of today (id-1)
       REAL(KIND(1D0)) :: qn_use ! qn used in OHM calculations [W m-2]
+      INTEGER :: is
 
       ASSOCIATE ( &
          atmState => modState%atmState, &
@@ -1819,7 +1822,8 @@ SUBROUTINE SUEWS_cal_Qs( &
                IF (StorageHeatMethod == 0) THEN !Use observed QS
                   qs = qs_obs
 
-               ELSEIF (StorageHeatMethod == 1 .OR. StorageHeatMethod == 6 .OR. StorageHeatMethod == 7) THEN !Use OHM to calculate QS
+               ELSEIF (StorageHeatMethod == 1 .OR. StorageHeatMethod == 6 .OR. StorageHeatMethod == 7 .OR. &
+                      (StorageHeatMethod == 3 .AND. .NOT. ohmState%anohm_coeff_ready)) THEN !Use OHM to calculate QS or fallback for ANOHM spin-up
                   Tair_mav_5d = HDD_id(10)
                   IF (Diagnose == 1) WRITE (*, *) 'Calling OHM...'
                   CALL OHM(qn_use, ohmState%qn_av, ohmState%dqndt, &
@@ -1853,27 +1857,34 @@ SUBROUTINE SUEWS_cal_Qs( &
                   QS_roof = qs
                   QS_wall = qs
 
-                  ! use AnOHM to calculate QS, TS 14 Mar 2016
-                  ! disable AnOHM, TS 20 Jul 2023
-                  ! ELSEIF (StorageHeatMethod == 3) THEN !
-                  !    IF (Diagnose == 1) WRITE (*, *) 'Calling AnOHM...'
-                  !    ! CALL AnOHM(qn1_use,qn1_store_grid,qn1_av_store_grid,qf,&
-                  !    !      MetForcingData_grid,state_id/StoreDrainPrm(6,:),&
-                  !    !      alb, emis, cpAnOHM, kkAnOHM, chAnOHM,&
-                  !    !      sfr_surf,nsurf,nsh,EmissionsMethod,id,Gridiv,&
-                  !    !      a1,a2,a3,qs,deltaQi)
-                  !    moist_surf = state_id/StoreDrainPrm(6, :)
-                  !    ! CALL AnOHM( &
-                  !    !    tstep, dt_since_start, &
-                  !    !    qn_use, qn_av_prev, dqndt_prev, qf, &
-                  !    !    MetForcingData_grid, moist_surf, &
-                  !    !    alb, emis, cpAnOHM, kkAnOHM, chAnOHM, & ! input
-                  !    !    sfr_surf, nsurf, EmissionsMethod, id, Gridiv, &
-                  !    !    qn_av_next, dqndt_next, &
-                  !    !    a1, a2, a3, qs, deltaQi) ! output
-                  !    QS_surf = qs
-                  !    QS_roof = qs
-                  !    QS_wall = qs
+               ELSEIF (StorageHeatMethod == 3) THEN
+                  DO is = 1, nsurf
+                     IF (StoreDrainPrm(6, is) > 0.0D0) THEN
+                        moist_surf(is) = MAX(0.0D0, MIN(1.0D0, state_id(is)/StoreDrainPrm(6, is)))
+                     ELSE
+                        moist_surf(is) = 0.0D0
+                     END IF
+                  END DO
+
+                  IF (Diagnose == 1) WRITE (*, *) 'Calling AnOHM...'
+                  CALL AnOHM( &
+                     tstep, dt_since_start, &
+                     qn_use, ohmState%qn_av, ohmState%dqndt, &
+                     ohmState%anohm_coeff_sd(1:ohmState%anohm_coeff_count), &
+                     ohmState%anohm_coeff_ta(1:ohmState%anohm_coeff_count), &
+                     ohmState%anohm_coeff_rh(1:ohmState%anohm_coeff_count), &
+                     ohmState%anohm_coeff_pres(1:ohmState%anohm_coeff_count), &
+                     ohmState%anohm_coeff_ws(1:ohmState%anohm_coeff_count), &
+                     ohmState%anohm_coeff_ah(1:ohmState%anohm_coeff_count), &
+                     ohmState%anohm_coeff_tHr(1:ohmState%anohm_coeff_count), &
+                     moist_surf, &
+                     alb, emis, cpAnOHM, kkAnOHM, chAnOHM, &
+                     sfr_surf, nsurf, ohmState%anohm_coeff_day, Gridiv, &
+                     ohmState%qn_av, ohmState%dqndt, &
+                     a1, a2, a3, qs, deltaQi)
+                  QS_surf = qs
+                  QS_roof = qs
+                  QS_wall = qs
 
                   ! !Calculate QS using ESTM
                ELSEIF (StorageHeatMethod == 4 .OR. StorageHeatMethod == 14) THEN
@@ -1942,6 +1953,85 @@ SUBROUTINE SUEWS_cal_Qs( &
    END SUBROUTINE SUEWS_cal_Qs
 !=======================================================================
 
+   SUBROUTINE anohm_rollover_if_needed(timer, ohmState)
+      TYPE(SUEWS_TIMER), INTENT(IN) :: timer
+      TYPE(OHM_STATE), INTENT(INOUT) :: ohmState
+
+      IF (ohmState%anohm_working_day == -999) THEN
+         ohmState%anohm_working_day = timer%id
+      END IF
+
+      IF (timer%id /= ohmState%anohm_working_day) THEN
+         IF (ohmState%anohm_working_count > 0) THEN
+            ohmState%anohm_coeff_tHr = ohmState%anohm_working_tHr
+            ohmState%anohm_coeff_sd = ohmState%anohm_working_sd
+            ohmState%anohm_coeff_ta = ohmState%anohm_working_ta
+            ohmState%anohm_coeff_rh = ohmState%anohm_working_rh
+            ohmState%anohm_coeff_pres = ohmState%anohm_working_pres
+            ohmState%anohm_coeff_ws = ohmState%anohm_working_ws
+            ohmState%anohm_coeff_ah = ohmState%anohm_working_ah
+            ohmState%anohm_coeff_day = ohmState%anohm_working_day
+            ohmState%anohm_coeff_count = MIN(ohmState%anohm_working_count, ANOHM_MAX_SAMPLES)
+            ohmState%anohm_coeff_ready = (ohmState%anohm_coeff_count >= ANOHM_MIN_VALID_SAMPLES)
+         ELSE
+            ohmState%anohm_coeff_ready = .FALSE.
+            ohmState%anohm_coeff_count = 0
+            ohmState%anohm_coeff_day = ohmState%anohm_working_day
+            ohmState%anohm_coeff_tHr = -999.0D0
+            ohmState%anohm_coeff_sd = -999.0D0
+            ohmState%anohm_coeff_ta = -999.0D0
+            ohmState%anohm_coeff_rh = -999.0D0
+            ohmState%anohm_coeff_pres = -999.0D0
+            ohmState%anohm_coeff_ws = -999.0D0
+            ohmState%anohm_coeff_ah = -999.0D0
+         END IF
+
+         ohmState%anohm_working_day = timer%id
+         ohmState%anohm_working_count = 0
+         ohmState%anohm_working_tHr = -999.0D0
+         ohmState%anohm_working_sd = -999.0D0
+         ohmState%anohm_working_ta = -999.0D0
+         ohmState%anohm_working_rh = -999.0D0
+         ohmState%anohm_working_pres = -999.0D0
+         ohmState%anohm_working_ws = -999.0D0
+         ohmState%anohm_working_ah = -999.0D0
+      END IF
+
+   END SUBROUTINE anohm_rollover_if_needed
+
+   SUBROUTINE anohm_append_sample(timer, forcing, qf, ohmState)
+      TYPE(SUEWS_TIMER), INTENT(IN) :: timer
+      TYPE(SUEWS_FORCING), INTENT(IN) :: forcing
+      REAL(KIND(1D0)), INTENT(IN) :: qf
+      TYPE(OHM_STATE), INTENT(INOUT) :: ohmState
+
+      INTEGER :: idx
+
+      IF (ohmState%anohm_working_day == -999) THEN
+         ohmState%anohm_working_day = timer%id
+      END IF
+
+      IF (timer%id /= ohmState%anohm_working_day) THEN
+         CALL anohm_rollover_if_needed(timer, ohmState)
+      END IF
+
+      IF (timer%imin /= 0 .OR. timer%isec /= 0) RETURN
+
+      idx = MIN(timer%it + 1, ANOHM_MAX_SAMPLES)
+      ohmState%anohm_working_tHr(idx) = REAL(timer%it, KIND(1D0))
+      ohmState%anohm_working_sd(idx) = forcing%kdown
+      ohmState%anohm_working_ta(idx) = forcing%temp_c
+      ohmState%anohm_working_rh(idx) = forcing%RH
+      ohmState%anohm_working_pres(idx) = forcing%pres
+      ohmState%anohm_working_ws(idx) = forcing%U
+      ohmState%anohm_working_ah(idx) = qf
+
+      IF (idx > ohmState%anohm_working_count) THEN
+         ohmState%anohm_working_count = idx
+      END IF
+
+   END SUBROUTINE anohm_append_sample
+
 !==========================drainage and runoff================================
    SUBROUTINE SUEWS_cal_Water( &
       timer, config, forcing, siteInfo, & ! input
@@ -2005,6 +2095,10 @@ SUBROUTINE SUEWS_cal_Water( &
 
             state_id = hydroState%state_surf
             StoreDrainPrm = phenState%StoreDrainPrm
+            IF (config%StorageHeatMethod == 3) THEN
+               CALL anohm_rollover_if_needed(timer, modState%ohmState)
+               CALL anohm_append_sample(timer, forcing, modState%heatState%qf, modState%ohmState)
+            END IF
 
             ! sfr_surf = [pavedPrm%sfr, bldgPrm%sfr, evetrPrm%sfr, dectrPrm%sfr, grassPrm%sfr, bsoilPrm%sfr, waterPrm%sfr]
             WaterDist(1, 1) = pavedPrm%waterdist%to_paved

src/suews/src/suews_ctrl_type.f95

@@ -51,6 +51,8 @@ MODULE module_ctrl_type
       solar_State, atm_state
 
    IMPLICIT NONE
+   INTEGER, PARAMETER, PUBLIC :: ANOHM_MAX_SAMPLES = 48
+   INTEGER, PARAMETER, PUBLIC :: ANOHM_MIN_VALID_SAMPLES = 6
    ! in the following, the type definitions starting with `SUEWS_` are used in the main program
 
    ! ********** SUEWS_parameters schema (basic) **********