Skip to content

Commit

Permalink
Merge pull request #274 from UMEP-dev/sunt05/issue271
Browse files Browse the repository at this point in the history 
Sunt05/issue271
  • Loading branch information
@sunt05
sunt05 authored Jul 4, 2024
2 parents 1f83e94 + 03e9884 commit f2212f3
Show file tree
Hide file tree
Showing 9 changed files with 142 additions and 149 deletions.
15 changes: 7 additions & 8 deletions .github/workflows/build-publish_to_pypi.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,13 +23,12 @@ jobs:
strategy:
matrix:
buildplat:
- [ubuntu-latest, manylinux, x86_64]
- [macos-11, macosx, x86_64]
- [windows-2019, win, AMD64]
- [macos-14, macosx, arm64]

python: ["cp39", "cp310", "cp311", "cp312"]
- [ubuntu-latest, manylinux, x86_64]
- [macos-11, macosx, x86_64]
- [macos-14, macosx, arm64]
- [windows-2019, win, AMD64]

python: ["cp310", "cp311", "cp312"]

fail-fast: false

Expand All @@ -45,7 +44,7 @@ jobs:
- name: Set up Python
uses: actions/setup-python@v5
with:
python-version: "3.12"
python-version: 3.12

- name: Build wheels
uses: pypa/[email protected]
Expand All @@ -54,7 +53,7 @@ jobs:
CIBW_ARCHS: ${{ matrix.buildplat[2] }}
CIBW_ENVIRONMENT_PASS_LINUX: RUNNER_OS
CIBW_BEFORE_ALL_MACOS: "brew install gfortran && brew unlink gfortran && brew link gfortran"
CIBW_BEFORE_BUILD_WINDOWS: "pip install delvewheel" # Use delvewheel on windows
CIBW_BEFORE_BUILD_WINDOWS: "pip install delvewheel" # Use delvewheel on windows
CIBW_REPAIR_WHEEL_COMMAND_WINDOWS: "delvewheel repair -w {dest_dir} {wheel}"
CIBW_TEST_REQUIRES: pytest
CIBW_TEST_COMMAND_MACOS: "python -m pytest '{project}/test'"
Expand Down
5 changes: 5 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,3 +36,8 @@

- 31 May 2024:
- [feature] Added `dict_debug` an optional output of `run_supy` to help debug the model (for developers: add a `debug` flag to `df_state` to activate this feature) (#233)

- 04 Jul 2024:
- [bugfix] Fixed a bug causing an abrupt change in results due to a less smooth transition in `z0` from surfaces without roughness elements to those with them. (#271)
- [bugfix] Improved the discretisation of the vertical levels in the RSL scheme for better interpolation of surface diagnostics (e.g. `T2`) (#271)
- [maintenance] Added support for NumPy 2.0 (#271)
19 changes: 3 additions & 16 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,22 +6,6 @@ This is a public repo for SUEWS source code and documentation.

---

- [SUEWS](#suews)
- [Documentation](#documentation)
- [Compilation](#compilation)
- [Sample Run](#sample-run)
- [Developer Note](#developer-note)
- [Branch](#branch)
- [`master` branch](#master-branch)
- [Manual](#manual)
- [Test](#test)
- [Tests and purposes](#tests-and-purposes)
- [Workflow](#workflow)
- [Preparation of tests](#preparation-of-tests)
- [Debugging with GDB](#debugging-with-gdb)
- [GDB on macOS](#gdb-on-macos)
- [debugging with GDB](#debugging-with-gdb-1)
- [Questions](#questions)


## Documentation
Expand Down Expand Up @@ -98,6 +82,9 @@ pip show supy


## Developer Note

:note: **the follow is deprecated and will be updated**

- When doing `pip install -e supy-driver` using WSL in VS Code on Windows 10 I got the error "[Errno 13] Permission denied: 'build/bdist.linux-x86_64/wheel/supy_driver-2021a2.dist-info'". The solution was in the Windows file explorer to right-click the project directory (SUEWS) -> properties -> security -> edit -> everyone -> tick allow -> apply.

### Branch
Expand Down
11 changes: 4 additions & 7 deletions pyproject.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,8 @@
requires = [
"wheel",
"pytest",
"f90wrap==0.2.14",
"oldest-supported-numpy",
"f90wrap==0.2.15",
"numpy>=2.0",
"meson-python>=0.12.0",
]
build-backend = 'mesonpy'
Expand All @@ -22,13 +22,11 @@ dynamic = ["version"]
license = { text = "GPL-V3.0" }

dependencies = [
"pandas< 1.5; python_version <= '3.9'", # to fix scipy dependency issue in UMEP under QGIS3 wtih python 3.9
"pandas; python_version > '3.9'",
"importlib_resources; python_version < '3.9'", # to fix importlib issue in UMEP under QGIS3
"pandas",
"matplotlib",
"chardet",
"scipy",
"f90wrap==0.2.14", # f90wrap is required for f2py-based supy driver
"f90wrap==0.2.15", # f90wrap is required for f2py-based supy driver
"dask", # needs dask for parallel tasks
"f90nml", # utility for namelist files
"seaborn", # stat plotting
Expand All @@ -44,7 +42,6 @@ dependencies = [
]
classifiers = [
"Programming Language :: Python :: 3 :: Only",
"Programming Language :: Python :: 3.9",
"Programming Language :: Python :: 3.10",
"Programming Language :: Python :: 3.11",
"Programming Language :: Python :: 3.12",
Expand Down
98 changes: 61 additions & 37 deletions src/suews/src/suews_phys_resist.f95
View file
Original file line number Diff line number Diff line change
Expand Up @@ -461,23 +461,33 @@ SUBROUTINE SUEWS_cal_RoughnessParameters( &

INTEGER, PARAMETER :: notUsedI = -55
REAL(KIND(1D0)), PARAMETER :: notUsed = -55.5
REAL(KIND(1D0)) :: z0m4Paved, z0m4Grass, z0m4BSoil, z0m4Water !Default values for roughness lengths [m]
REAL(KIND(1D0)) :: z0m_Paved, z0m_Grass, z0m_BSoil, z0m_Water !Default values for roughness lengths [m]

! calculated values of FAI
REAL(KIND(1D0)) :: FAIBldg_use
REAL(KIND(1D0)) :: FAIEveTree_use
REAL(KIND(1D0)) :: FAIDecTree_use
REAL(KIND(1D0)) :: z0m_zh ! z0m for roughness elements (i.e. zh>0)
REAL(KIND(1D0)) :: zdm_zh ! zdm for roughness elements (i.e. zh>0)
REAL(KIND(1D0)) :: z0m_zh0 ! z0m for non-roughness elements (i.e. zh=0)
REAL(KIND(1D0)) :: zdm_zh0 ! zdm for non-roughness elements (i.e. zh=0)

!Total area of buildings and trees
! areaZh = (sfr_surf(BldgSurf) + sfr_surf(ConifSurf) + sfr_surf(DecidSurf))
! TS 19 Jun 2022: take porosity of trees into account; to be consistent with PAI calculation in RSL
PAI = DOT_PRODUCT(sfr_surf([BldgSurf, ConifSurf, DecidSurf]), [1D0, 1 - porosity_evetr, 1 - porosity_dectr])

! Set default values (using Moene & van Dam 2013, Atmos-Veg-Soil Interactions, Table 3.3)
Z0m4Paved = 0.003 !estimate
Z0m4Grass = 0.02
Z0m4BSoil = 0.002
Z0m4Water = 0.0005
z0m_Paved = 0.003 !estimate
z0m_Grass = 0.02
z0m_BSoil = 0.002
z0m_Water = 0.0005
! z0m for non-roughness elements (i.e. zh=0)
z0m_zh0 = (z0m_Paved*sfr_surf(PavSurf) &
+ z0m_Grass*sfr_surf(GrassSurf) &
+ z0m_BSoil*sfr_surf(BSoilSurf) &
+ z0m_Water*sfr_surf(WaterSurf))/(1 - PAI)
zdm_zh0 = 0

!------------------------------------------------------------------------------
!If total area of buildings and trees is larger than zero, use tree heights and building heights to calculate zH and FAI
Expand Down Expand Up @@ -518,29 +528,30 @@ SUBROUTINE SUEWS_cal_RoughnessParameters( &
IF (Zh /= 0) THEN
!Calculate z0m and zdm depending on the Z0 method
IF (RoughLenMomMethod == 2) THEN !Rule of thumb (G&O 1999)
z0m = 0.1*Zh
zdm = 0.7*Zh
z0m_zh = 0.1*Zh
zdm_zh = 0.7*Zh
ELSEIF (RoughLenMomMethod == 3) THEN !MacDonald 1998
! zdm = (1 + 4.43**(-sfr_surf(BldgSurf))*(sfr_surf(BldgSurf) - 1))*Zh
zdm = (1 + 4.43**(-PAI)*(PAI - 1))*Zh ! changed sfr_surf(BldgSurf) to PAI to be consistent with z0m calculation, TS 04 Jun 2023
z0m = ((1 - zdm/Zh)*EXP(-(0.5*1.0*1.2/0.4**2*(1 - zdm/Zh)*FAI)**(-0.5)))*Zh
zdm_zh = (1 + 4.43**(-PAI)*(PAI - 1))*Zh ! changed sfr_surf(BldgSurf) to PAI to be consistent with z0m calculation, TS 04 Jun 2023
z0m_zh = ((1 - zdm/Zh)*EXP(-(0.5*1.0*1.2/0.4**2*(1 - zdm/Zh)*FAI)**(-0.5)))*Zh
ELSEIF (RoughLenMomMethod == 4) THEN ! lambdaP dependent as in Fig.1a of G&O (1999)
! these are derived using digitalised points
zdm = (-0.182 + 0.722*sigmoid(-1.16 + 3.89*PAI) + 0.493*sigmoid(-5.17 + 32.7*PAI))*Zh
z0m = (0.00208 + &
0.0165*MIN(PAI, .7) + 2.52*MIN(PAI, .7)**2 + &
3.21*MIN(PAI, .7)**3 - 43.6*MIN(PAI, .7)**4 + &
76.5*MIN(PAI, .7)**5 - 40.*MIN(PAI, .7)**6)*Zh
zdm_zh = (-0.182 + 0.722*sigmoid(-1.16 + 3.89*PAI) + 0.493*sigmoid(-5.17 + 32.7*PAI))*Zh
z0m_zh = (0.00208 + &
0.0165*MIN(PAI, .7) + 2.52*MIN(PAI, .7)**2 + &
3.21*MIN(PAI, .7)**3 - 43.6*MIN(PAI, .7)**4 + &
76.5*MIN(PAI, .7)**5 - 40.*MIN(PAI, .7)**6)*Zh
END IF
! #271: to smooth the z0m (and zdm) values when other non-rough surfaces are present
z0m = DOT_PRODUCT([z0m_zh, z0m_zh0], [PAI, 1 - PAI])
zdm = DOT_PRODUCT([zdm_zh, zdm_zh0], [PAI, 1 - PAI])

ELSEIF (Zh == 0) THEN !If zh calculated to be zero, set default roughness length and displacement height
IF (PAI /= 0) CALL ErrorHint(15, 'In SUEWS_RoughnessParameters.f95, zh = 0 m but areaZh > 0', zh, PAI, notUsedI)
!Estimate z0 and zd using default values and surfaces that do not contribute to areaZh
IF (PAI /= 1) THEN
z0m = (z0m4Paved*sfr_surf(PavSurf) &
+ z0m4Grass*sfr_surf(GrassSurf) &
+ z0m4BSoil*sfr_surf(BSoilSurf) &
+ z0m4Water*sfr_surf(WaterSurf))/(1 - PAI)
zdm = 0
z0m = z0m_zh0
zdm = zdm_zh0
CALL ErrorHint(15, 'Setting z0m and zdm using default values', z0m, zdm, notUsedI)
ELSEIF (PAI == 1) THEN !If, for some reason, Zh = 0 and areaZh == 1, assume height of 10 m and use rule-of-thumb
z0m = 1
Expand Down Expand Up @@ -642,10 +653,15 @@ SUBROUTINE SUEWS_cal_RoughnessParameters_DTS( &

!Default values for roughness lengths [m]
! Set default values (using Moene & van Dam 2013, Atmos-Veg-Soil Interactions, Table 3.3)
REAL(KIND(1D0)), PARAMETER :: z0m4Paved = 0.003 !estimate
REAL(KIND(1D0)), PARAMETER :: z0m4Grass = 0.02
REAL(KIND(1D0)), PARAMETER :: z0m4BSoil = 0.002
REAL(KIND(1D0)), PARAMETER :: z0m4Water = 0.0005
REAL(KIND(1D0)), PARAMETER :: z0m_Paved = 0.003 !estimate
REAL(KIND(1D0)), PARAMETER :: z0m_Grass = 0.02
REAL(KIND(1D0)), PARAMETER :: z0m_BSoil = 0.002
REAL(KIND(1D0)), PARAMETER :: z0m_Water = 0.0005

REAL(KIND(1D0)) :: z0m_zh ! z0m for roughness elements (i.e. zh>0)
REAL(KIND(1D0)) :: zdm_zh ! zdm for roughness elements (i.e. zh>0)
REAL(KIND(1D0)) :: z0m_zh0 ! z0m for non-roughness elements (i.e. zh=0)
REAL(KIND(1D0)) :: zdm_zh0 ! zdm for non-roughness elements (i.e. zh=0)

! calculated values of FAI
! REAL(KIND(1D0)), INTENT(out) :: FAIBldg_use
Expand Down Expand Up @@ -712,6 +728,13 @@ SUBROUTINE SUEWS_cal_RoughnessParameters_DTS( &
! TS 19 Jun 2022: take porosity of trees into account; to be consistent with PAI calculation in RSL
PAI = DOT_PRODUCT(sfr_surf([BldgSurf, ConifSurf, DecidSurf]), [1D0, 1 - porosity_evetr, 1 - porosity_dectr])

! z0m for non-roughness elements (i.e. zh=0)
z0m_zh0 = (z0m_Paved*sfr_surf(PavSurf) &
+ z0m_Grass*sfr_surf(GrassSurf) &
+ z0m_BSoil*sfr_surf(BSoilSurf) &
+ z0m_Water*sfr_surf(WaterSurf))/(1 - PAI)
zdm_zh0 = 0

!------------------------------------------------------------------------------
!If total area of buildings and trees is larger than zero, use tree heights and building heights to calculate zH and FAI
IF (PAI /= 0) THEN
Expand Down Expand Up @@ -752,28 +775,29 @@ SUBROUTINE SUEWS_cal_RoughnessParameters_DTS( &
IF (Zh /= 0) THEN
!Calculate z0m and zdm depending on the Z0 method
IF (RoughLenMomMethod == 2) THEN !Rule of thumb (G&O 1999)
z0m = 0.1*Zh
zdm = 0.7*Zh
z0m_zh = 0.1*Zh
zdm_zh = 0.7*Zh
ELSEIF (RoughLenMomMethod == 3) THEN !MacDonald 1998
zdm = (1 + 4.43**(-sfr_surf(BldgSurf))*(sfr_surf(BldgSurf) - 1))*Zh
z0m = ((1 - zdm/Zh)*EXP(-(0.5*1.0*1.2/0.4**2*(1 - zdm/Zh)*FAI)**(-0.5)))*Zh
zdm_zh = (1 + 4.43**(-sfr_surf(BldgSurf))*(sfr_surf(BldgSurf) - 1))*Zh
z0m_zh = ((1 - zdm/Zh)*EXP(-(0.5*1.0*1.2/0.4**2*(1 - zdm/Zh)*FAI)**(-0.5)))*Zh
ELSEIF (RoughLenMomMethod == 4) THEN ! lambdaP dependent as in Fig.1a of G&O (1999)
! these are derived using digitalised points
zdm = (-0.182 + 0.722*sigmoid(-1.16 + 3.89*PAI) + 0.493*sigmoid(-5.17 + 32.7*PAI))*Zh
z0m = (0.00208 + &
0.0165*MIN(PAI, .7) + 2.52*MIN(PAI, .7)**2 + &
3.21*MIN(PAI, .7)**3 - 43.6*MIN(PAI, .7)**4 + &
76.5*MIN(PAI, .7)**5 - 40.*MIN(PAI, .7)**6)*Zh
zdm_zh = (-0.182 + 0.722*sigmoid(-1.16 + 3.89*PAI) + 0.493*sigmoid(-5.17 + 32.7*PAI))*Zh
z0m_zh = (0.00208 + &
0.0165*MIN(PAI, .7) + 2.52*MIN(PAI, .7)**2 + &
3.21*MIN(PAI, .7)**3 - 43.6*MIN(PAI, .7)**4 + &
76.5*MIN(PAI, .7)**5 - 40.*MIN(PAI, .7)**6)*Zh
END IF
! #271: to smooth the z0m (and zdm) values when other non-rough surfaces are present
z0m = DOT_PRODUCT([z0m_zh, z0m_zh0], [PAI, 1 - PAI])
zdm = DOT_PRODUCT([zdm_zh, zdm_zh0], [PAI, 1 - PAI])

ELSEIF (Zh == 0) THEN !If zh calculated to be zero, set default roughness length and displacement height
IF (PAI /= 0) CALL ErrorHint(15, 'In SUEWS_RoughnessParameters.f95, zh = 0 m but areaZh > 0', zh, PAI, notUsedI)
!Estimate z0 and zd using default values and surfaces that do not contribute to areaZh
IF (PAI /= 1) THEN
z0m = (z0m4Paved*sfr_surf(PavSurf) &
+ z0m4Grass*sfr_surf(GrassSurf) &
+ z0m4BSoil*sfr_surf(BSoilSurf) &
+ z0m4Water*sfr_surf(WaterSurf))/(1 - PAI)
zdm = 0
z0m = z0m_zh0
zdm = zdm_zh0
CALL ErrorHint(15, 'Setting z0m and zdm using default values', z0m, zdm, notUsedI)
ELSEIF (PAI == 1) THEN !If, for some reason, Zh = 0 and areaZh == 1, assume height of 10 m and use rule-of-thumb
z0m = 1
Expand Down
Loading

0 comments on commit f2212f3

Please sign in to comment.