From 21423faa04d0ff854940620253a78e5ca8e9bf2d Mon Sep 17 00:00:00 2001 From: Ting Sun Date: Wed, 19 Nov 2025 14:10:45 +0000 Subject: [PATCH] feat: Re-enable 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 --- CHANGELOG.md | 6 + .../csv-table/StorageHeatMethod.csv | 2 +- .../parameterisations-and-sub-models.rst | 1 + src/suews/src/suews_ctrl_driver.f95 | 142 +++++++++-- src/suews/src/suews_ctrl_type.f95 | 2 + src/suews/src/suews_phys_anohm.f95 | 225 ++++-------------- src/suews/src/suews_type_surface.f95 | 24 ++ 7 files changed, 195 insertions(+), 207 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 48ceea6a0..7af716a15 100644 --- a/CHANGELOG.md +++ b/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()` diff --git a/docs/source/inputs/tables/RunControl/csv-table/StorageHeatMethod.csv b/docs/source/inputs/tables/RunControl/csv-table/StorageHeatMethod.csv index 2f6ef9b39..397f0e8c7 100644 --- a/docs/source/inputs/tables/RunControl/csv-table/StorageHeatMethod.csv +++ b/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|" diff --git a/docs/source/parameterisations-and-sub-models.rst b/docs/source/parameterisations-and-sub-models.rst index d2f03e714..bb68dbb4d 100644 --- a/docs/source/parameterisations-and-sub-models.rst +++ b/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. diff --git a/src/suews/src/suews_ctrl_driver.f95 b/src/suews/src/suews_ctrl_driver.f95 index bb9538966..5569c7879 100644 --- a/src/suews/src/suews_ctrl_driver.f95 +++ b/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,6 +1497,7 @@ 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 @@ -1503,6 +1505,7 @@ SUBROUTINE SUEWS_cal_Qs( & ! 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 diff --git a/src/suews/src/suews_ctrl_type.f95 b/src/suews/src/suews_ctrl_type.f95 index 743f14a2c..271247503 100644 --- a/src/suews/src/suews_ctrl_type.f95 +++ b/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) ********** diff --git a/src/suews/src/suews_phys_anohm.f95 b/src/suews/src/suews_phys_anohm.f95 index 4c561c30c..795092f20 100644 --- a/src/suews/src/suews_phys_anohm.f95 +++ b/src/suews/src/suews_phys_anohm.f95 @@ -29,10 +29,11 @@ MODULE module_phys_anohm !! -# grid ensemble OHM coefficients: a1, a2 and a3 SUBROUTINE AnOHM( & tstep, dt_since_start, & - qn1, qn_av_prev, dqndt_prev, qf, & - MetForcingData_grid, moist_surf, & + qn1, qn_av_prev, dqndt_prev, & + Sd_series, Ta_series, RH_series, pres_series, WS_series, AH_series, tHr_series, & + moist_surf, & alb, emis, cpAnOHM, kkAnOHM, chAnOHM, & ! input - sfr_surf, nsurf, EmissionsMethod, id, Gridiv, & + sfr_surf, nsurf, forcing_day, Gridiv, & qn_av_next, dqndt_next, & a1, a2, a3, qs, deltaQi) ! output @@ -40,10 +41,14 @@ SUBROUTINE AnOHM( & INTEGER, INTENT(in) :: tstep ! time step [s] INTEGER, INTENT(in) :: dt_since_start ! time since simulation starts [s] - REAL(KIND(1D0)), INTENT(in), DIMENSION(:, :) :: MetForcingData_grid !< met forcing array of grid - REAL(KIND(1D0)), INTENT(in) :: qn1 !< net all-wave radiation [W m-2] - REAL(KIND(1D0)), INTENT(in) :: qf !< anthropogenic heat flux [W m-2] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: Sd_series !< shortwave radiation time series [W m-2] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: Ta_series !< air temperature series [degC] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: RH_series !< relative humidity series [%] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: pres_series !< pressure series [hPa] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: WS_series !< wind speed series [m s-1] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: AH_series !< anthropogenic heat series [W m-2] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: tHr_series !< time of day series [h] REAL(KIND(1D0)), INTENT(in) :: sfr_surf(nsurf) !< surface fraction (0-1) [-] REAL(KIND(1D0)), INTENT(in) :: moist_surf(nsurf) !< non-dimensional surface wetness status (0-1) [-] @@ -53,9 +58,8 @@ SUBROUTINE AnOHM( & REAL(KIND(1D0)), INTENT(in), DIMENSION(:) :: kkAnOHM !< thermal conductivity [W m-1 K-1] REAL(KIND(1D0)), INTENT(in), DIMENSION(:) :: chAnOHM !< bulk transfer coef [J m-3 K-1] - INTEGER, INTENT(in) :: id !< day of year [-] + INTEGER, INTENT(in) :: forcing_day !< day of year that forcing series represents [-] INTEGER, INTENT(in) :: Gridiv !< grid id [-] - INTEGER, INTENT(in) :: EmissionsMethod !< AnthropHeat option [-] INTEGER, INTENT(in) :: nsurf !< number of surfaces [-] ! INTEGER,INTENT(in):: nsh !< number of timesteps in one hour [-] @@ -74,7 +78,7 @@ SUBROUTINE AnOHM( & REAL(KIND(1D0)), INTENT(out) :: deltaQi(nsurf) !< storage heat flux of snow surfaces INTEGER :: is, xid !< @var qn1 net all-wave radiation - INTEGER, SAVE :: id_save ! store index of the valid day with enough data ! TODO: Remove SAVE states from the model + INTEGER, SAVE :: id_save = -999 ! store index of the valid day with enough data ! TODO: Remove SAVE states from the model REAL(KIND(1D0)), PARAMETER :: NotUsed = -55.5 !< @var qn1 net all-wave radiation INTEGER, PARAMETER :: notUsedI = -55 !< @var qn1 net all-wave radiation LOGICAL :: idQ ! whether id contains enough data @@ -93,26 +97,28 @@ SUBROUTINE AnOHM( & ! to test if the current met block contains enough data for AnOHM ! TODO: more robust selection should be implemented ! daylight hours >= 6 - idQ = COUNT(MetForcingData_grid(:, 2) == id .AND. & ! day of year - MetForcingData_grid(:, 4) == 0 .AND. & ! minutes - MetForcingData_grid(:, 15) > 0) & ! Sd - >= 6 + idQ = COUNT(Sd_series > 5.0D0) >= 6 ! PRINT*, idQ IF (idQ) THEN ! given enough data, calculate coefficients of day `id` - xid = id - id_save = id ! store index of the valid day with enough data - ELSE - ! otherwise calculate coefficients of yesterday: `id-1` + xid = forcing_day + id_save = forcing_day ! store index of the valid day with enough data + ELSEIF (id_save /= -999) THEN + ! otherwise use coefficients stored for the last valid day xid = id_save + ELSE + xid = forcing_day END IF DO is = 1, nsurf ! call AnOHM to calculate the coefs. - CALL AnOHM_coef(is, xid, Gridiv, MetForcingData_grid, moist_surf, EmissionsMethod, qf, & !input - alb, emis, cpAnOHM, kkAnOHM, chAnOHM, & ! input - xa1(is), xa2(is), xa3(is)) ! output + CALL AnOHM_coef( & + is, xid, Gridiv, & + Sd_series, Ta_series, RH_series, pres_series, WS_series, AH_series, tHr_series, & + moist_surf, & + alb, emis, cpAnOHM, kkAnOHM, chAnOHM, & + xa1(is), xa2(is), xa3(is)) ! print*, 'AnOHM_coef are: ',xa1,xa2,xa3 END DO @@ -148,7 +154,8 @@ END SUBROUTINE AnOHM !! -# OHM coefficients of a given surface type: a1, a2 and a3 SUBROUTINE AnOHM_coef( & sfc_typ, xid, xgrid, & !input - MetForcingData_grid, moist, EmissionsMethod, qf, & !input + Sd_series, Ta_series, RH_series, pres_series, WS_series, AH_series, tHr_series, & !input + moist, & alb, emis, cpAnOHM, kkAnOHM, chAnOHM, & ! input xa1, xa2, xa3) ! output @@ -158,16 +165,20 @@ SUBROUTINE AnOHM_coef( & INTEGER, INTENT(in) :: sfc_typ !< surface type [-] INTEGER, INTENT(in) :: xid !< day of year [-] INTEGER, INTENT(in) :: xgrid !< grid id [-] - INTEGER, INTENT(in) :: EmissionsMethod !< AnthropHeat option [-] - REAL(KIND(1D0)), INTENT(in) :: qf !< anthropogenic heat flux [W m-2] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: Sd_series !< incoming solar radiation [W m-2] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: Ta_series !< air temperature [degC] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: RH_series !< relative humidity [%] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: pres_series !< Atmospheric pressure [hPa] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: WS_series ! wind speed [m s-1] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: AH_series ! anthropogenic heat [W m-2] + REAL(KIND(1D0)), DIMENSION(:), INTENT(in) :: tHr_series ! local time [hr] REAL(KIND(1D0)), INTENT(in), DIMENSION(:) :: alb !< albedo [-] REAL(KIND(1D0)), INTENT(in), DIMENSION(:) :: emis !< emissivity [-] REAL(KIND(1D0)), INTENT(in), DIMENSION(:) :: cpAnOHM !< heat capacity [J m-3 K-1] REAL(KIND(1D0)), INTENT(in), DIMENSION(:) :: kkAnOHM !< thermal conductivity [W m-1 K-1] REAL(KIND(1D0)), INTENT(in), DIMENSION(:) :: chAnOHM !< bulk transfer coef [J m-3 K-1] REAL(KIND(1D0)), INTENT(in), DIMENSION(:) :: moist !< surface wetness status [-] - REAL(KIND(1D0)), INTENT(in), DIMENSION(:, :) :: MetForcingData_grid !< met forcing array of grid ! output REAL(KIND(1D0)), INTENT(out) :: xa1 !< AnOHM coefficients of grid [-] @@ -186,14 +197,7 @@ SUBROUTINE AnOHM_coef( & REAL(KIND(1D0)) :: mWS, mWF, mAH !< mean values of WS, WF and AH ! forcings: - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: Sd ! incoming solar radiation [W m-2] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: Ta ! air temperature [degC] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: RH ! relative humidity [%] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: pres ! air pressure [hPa] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: WS ! wind speed [m s-1] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: WF ! water flux density [m3 s-1 m-2] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: AH ! anthropogenic heat [W m-2] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: tHr ! time [hr] + REAL(KIND(1D0)), DIMENSION(SIZE(Sd_series)) :: WF_series ! water flux density [m3 s-1 m-2] ! sfc. properties: REAL(KIND(1D0)) :: xalb ! albedo, @@ -209,7 +213,7 @@ SUBROUTINE AnOHM_coef( & ! locally saved variables: ! if coefficients have been calculated, just reload them ! otherwise, do the calculation - INTEGER, SAVE :: id_save, grid_save + INTEGER, SAVE :: id_save = -999, grid_save = -999 REAL(KIND(1D0)), SAVE :: coeff_grid_day(7, 3) = -999. ! PRINT*, 'xid,id_save',xid,id_save @@ -224,15 +228,11 @@ SUBROUTINE AnOHM_coef( & xa3 = coeff_grid_day(sfc_typ, 3) ELSE - ! load forcing characteristics: - CALL AnOHM_Fc( & - xid, MetForcingData_grid, EmissionsMethod, qf, & ! input - ASd, mSd, tSd, ATa, mTa, tTa, tau, mWS, mWF, mAH) ! output - - ! load forcing variables: - CALL AnOHM_FcLoad( & - xid, MetForcingData_grid, EmissionsMethod, qf, & ! input - Sd, Ta, RH, pres, WS, WF, AH, tHr) ! output + ! load forcing characteristics from provided series: + WF_series = 0.0D0 + CALL AnOHM_FcCal( & + Sd_series, Ta_series, WS_series, WF_series, AH_series, tHr_series, & + ASd, mSd, tSd, ATa, mTa, tTa, tau, mWS, mWF, mAH) ! load sfc. properties: xalb = alb(sfc_typ) @@ -246,7 +246,7 @@ SUBROUTINE AnOHM_coef( & ! calculate Bowen ratio: CALL AnOHM_Bo_cal( & sfc_typ, & - Sd, Ta, RH, pres, tHr, & ! input: forcing + Sd_series, Ta_series, RH_series, pres_series, tHr_series, & ! input: forcing ASd, mSd, ATa, mTa, tau, mWS, mWF, mAH, & ! input: forcing xalb, xemis, xcp, xk, xch, xmoist, & ! input: sfc properties tSd, & ! input: peaking time of Sd in hour @@ -663,145 +663,6 @@ SUBROUTINE AnOHM_coef_water_cal( & END SUBROUTINE AnOHM_coef_water_cal !======================================================================================== - !======================================================================================== - SUBROUTINE AnOHM_Fc( & - xid, MetForcingData_grid, EmissionsMethod, qf, & ! input - ASd, mSd, tSd, ATa, mTa, tTa, tau, mWS, mWF, mAH) ! output - - IMPLICIT NONE - - ! input - INTEGER, INTENT(in) :: xid - INTEGER, INTENT(in) :: EmissionsMethod - REAL(KIND(1D0)), INTENT(in), DIMENSION(:, :) :: MetForcingData_grid - REAL(KIND(1D0)), INTENT(in) :: qf !< anthropogenic heat flux [W m-2] - - ! output - REAL(KIND(1D0)), INTENT(out) :: ASd !< daily amplitude of solar radiation [W m-2] - REAL(KIND(1D0)), INTENT(out) :: mSd !< daily mean solar radiation [W m-2] - REAL(KIND(1D0)), INTENT(out) :: tSd !< local peaking time of solar radiation [hr] - REAL(KIND(1D0)), INTENT(out) :: ATa !< daily amplitude of air temperature [degC] - REAL(KIND(1D0)), INTENT(out) :: mTa !< daily mean air temperature [degC] - REAL(KIND(1D0)), INTENT(out) :: tTa !< local peaking time of air temperature [hour] - REAL(KIND(1D0)), INTENT(out) :: tau !< phase lag between Sd and Ta (Ta-Sd) [rad] - REAL(KIND(1D0)), INTENT(out) :: mWS !< daily mean wind speed [m s-1] - REAL(KIND(1D0)), INTENT(out) :: mWF !< daily mean underground moisture flux [m3 s-1 m-2] - REAL(KIND(1D0)), INTENT(out) :: mAH !< daily mean anthropogenic heat flux [W m-2] - - ! forcings: - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: Sd !< incoming solar radiation [W m-2] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: Ta !< air temperature [degC] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: RH !< relative humidity [%] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: pres !< Atmospheric pressure [hPa] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: WS ! wind speed [m s-1] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: WF ! water flux density [m3 s-1 m-2] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: AH ! anthropogenic heat [W m-2] - REAL(KIND(1D0)), DIMENSION(:), ALLOCATABLE :: tHr ! local time [hr] - - ! load forcing variables: - CALL AnOHM_FcLoad(xid, MetForcingData_grid, EmissionsMethod, qf, Sd, Ta, RH, pres, WS, WF, AH, tHr) - ! calculate forcing scales for AnOHM: - CALL AnOHM_FcCal(Sd, Ta, WS, WF, AH, tHr, ASd, mSd, tSd, ATa, mTa, tTa, tau, mWS, mWF, mAH) - - ! CALL r8vec_print(SIZE(sd, dim=1),sd,'Sd') - ! PRINT*, ASd,mSd,tSd - - ! CALL r8vec_print(SIZE(ta, dim=1),ta,'Ta') - ! PRINT*, ATa,mTa,tTa - - END SUBROUTINE AnOHM_Fc - !======================================================================================== - - !======================================================================================== - !> load forcing series for AnOHM_FcCal - SUBROUTINE AnOHM_FcLoad( & - xid, MetForcingData_grid, EmissionsMethod, qf, & ! input - Sd, Ta, RH, pres, WS, WF, AH, tHr) ! output - - IMPLICIT NONE - - ! input - INTEGER, INTENT(in) :: xid !< day of year - INTEGER, INTENT(in) :: EmissionsMethod !< AnthropHeat option - REAL(KIND(1D0)), INTENT(in), DIMENSION(:, :) :: MetForcingData_grid !< met forcing array of grid - REAL(KIND(1D0)), INTENT(in) :: qf !< anthropogenic heat flux [W m-2] - - ! output - REAL(KIND(1D0)), DIMENSION(:), INTENT(out), ALLOCATABLE :: Sd !< incoming solar radiation [W m-2] - REAL(KIND(1D0)), DIMENSION(:), INTENT(out), ALLOCATABLE :: Ta !< air temperature [degC] - REAL(KIND(1D0)), DIMENSION(:), INTENT(out), ALLOCATABLE :: RH !< relative humidity [%] - REAL(KIND(1D0)), DIMENSION(:), INTENT(out), ALLOCATABLE :: pres !< atmospheric pressure [mbar] - REAL(KIND(1D0)), DIMENSION(:), INTENT(out), ALLOCATABLE :: WS !< wind speed [m s-1] - REAL(KIND(1D0)), DIMENSION(:), INTENT(out), ALLOCATABLE :: WF !< water flux density [m3 s-1 m-2] - REAL(KIND(1D0)), DIMENSION(:), INTENT(out), ALLOCATABLE :: AH !< anthropogenic heat [W m-2] - REAL(KIND(1D0)), DIMENSION(:), INTENT(out), ALLOCATABLE :: tHr !< local time [hr] - - ! local variables: - REAL(KIND(1D0)), DIMENSION(:, :), ALLOCATABLE :: subMet ! subset array of daytime series - - INTEGER :: err - INTEGER :: lenMetData, nVar - - LOGICAL, ALLOCATABLE :: metMask(:) - - ! construct mask - IF (ALLOCATED(metMask)) DEALLOCATE (metMask, stat=err) - ALLOCATE (metMask(SIZE(MetForcingData_grid, dim=1))) - metMask = (MetForcingData_grid(:, 2) == xid & ! day=xid - .AND. MetForcingData_grid(:, 4) == 0) ! tmin=0 - - ! determine the length of subset - lenMetData = COUNT(metMask) - - ! construct array for time and met variables - nVar = 8 ! number of variables to retrieve - ! print*, 'good 1' - ! allocate subMet: - IF (ALLOCATED(subMet)) DEALLOCATE (subMet, stat=err) - ALLOCATE (subMet(lenMetData, nVar)) - subMet = RESHAPE(PACK(MetForcingData_grid(:, (/3, & !time: hour - 15, 12, 11, 13, 10, 12, 9/)), & ! met: Sd, Ta, RH, pres, WS, WF, AH - ! NB: WF not used: a placeholder - SPREAD(metMask, dim=2, ncopies=nVar)), & ! replicate mask vector to 2D array - (/lenMetData, nVar/)) ! convert to target shape - - ! re-allocate arrays as their sizes may change during passing - IF (ALLOCATED(tHr)) DEALLOCATE (tHr, stat=err) - ALLOCATE (tHr(lenMetData)) - IF (ALLOCATED(Sd)) DEALLOCATE (Sd, stat=err) - ALLOCATE (Sd(lenMetData)) - IF (ALLOCATED(Ta)) DEALLOCATE (Ta, stat=err) - ALLOCATE (Ta(lenMetData)) - IF (ALLOCATED(RH)) DEALLOCATE (RH, stat=err) - ALLOCATE (RH(lenMetData)) - IF (ALLOCATED(pres)) DEALLOCATE (pres, stat=err) - ALLOCATE (pres(lenMetData)) - IF (ALLOCATED(WS)) DEALLOCATE (WS, stat=err) - ALLOCATE (WS(lenMetData)) - IF (ALLOCATED(WF)) DEALLOCATE (WF, stat=err) - ALLOCATE (WF(lenMetData)) - IF (ALLOCATED(AH)) DEALLOCATE (AH, stat=err) - ALLOCATE (AH(lenMetData)) - - ! load the sublist into forcing variables - tHr = subMet(:, 1) ! time in hour - Sd = subMet(:, 2) - Ta = subMet(:, 3) - RH = subMet(:, 4) - pres = subMet(:, 5) - WS = subMet(:, 6) - WF = 0 ! set as 0 for the moment - IF (EmissionsMethod == 0) THEN - AH = subMet(:, 8) ! read in from MetForcingData_grid, - ELSE - ! AH = 0 ! temporarily change to zero; - AH = qf ! use modelled value - ! AH = mAH_grids(xid-1,xgrid) - END IF - - END SUBROUTINE AnOHM_FcLoad - !======================================================================================== - !======================================================================================== !> calculate the key parameters of a sinusoidal curve for AnOHM forcings !> i.e., a, b, c in a*Sin(Pi/12*t+b)+c diff --git a/src/suews/src/suews_type_surface.f95 b/src/suews/src/suews_type_surface.f95 index 983570048..90243bcb8 100644 --- a/src/suews/src/suews_type_surface.f95 +++ b/src/suews/src/suews_type_surface.f95 @@ -1,6 +1,8 @@ module module_type_surface + use module_ctrl_const_allocate, only: nsurf implicit none + INTEGER, PARAMETER, PUBLIC :: ANOHM_MAX_SAMPLES = 48 TYPE, PUBLIC :: LUMPS_PRM REAL(KIND(1D0)) :: raincover = 0.0D0 ! limit when surface totally covered with water for LUMPS [mm] @@ -59,6 +61,28 @@ module module_type_surface REAL(KIND(1D0)) :: a1_water = 0.0D0! Dynamic OHM coefficients of water REAL(KIND(1D0)) :: a2_water = 0.0D0! Dynamic OHM coefficients of water REAL(KIND(1D0)) :: a3_water = 0.0D0! Dynamic OHM coefficients of water + INTEGER :: anohm_working_day = -999 + INTEGER :: anohm_working_count = 0 + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_working_tHr = -999.0D0 + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_working_sd = -999.0D0 + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_working_ta = -999.0D0 + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_working_rh = -999.0D0 + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_working_pres = -999.0D0 + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_working_ws = -999.0D0 + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_working_ah = -999.0D0 + INTEGER :: anohm_coeff_day = -999 + INTEGER :: anohm_coeff_count = 0 + LOGICAL :: anohm_coeff_ready = .FALSE. + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_coeff_tHr = -999.0D0 + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_coeff_sd = -999.0D0 + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_coeff_ta = -999.0D0 + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_coeff_rh = -999.0D0 + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_coeff_pres = -999.0D0 + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_coeff_ws = -999.0D0 + REAL(KIND(1D0)), DIMENSION(ANOHM_MAX_SAMPLES) :: anohm_coeff_ah = -999.0D0 + REAL(KIND(1D0)), DIMENSION(nsurf) :: anohm_a1_surf = 0.0D0 + REAL(KIND(1D0)), DIMENSION(nsurf) :: anohm_a2_surf = 0.0D0 + REAL(KIND(1D0)), DIMENSION(nsurf) :: anohm_a3_surf = 0.0D0 ! flag for iteration safety - YES LOGICAL :: iter_safe = .TRUE. END TYPE OHM_STATE