Skip to content

INIT_GWDO: Revisions in computing mountain statistics #1408

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Songyou184 wants to merge 2 commits into
base: develop
Choose a base branch
from
+72 −32
Open

INIT_GWDO: Revisions in computing mountain statistics #1408

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
aa5e8d5
modified: src/core_init_atmosphere/Registry.xml
Songyou184 Feb 6, 2026
4189550
remove unused critical height variable (hc)
jihyeonjang Jun 18, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 7 additions & 1 deletion src/core_init_atmosphere/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -288,7 +288,7 @@
description="Whether to recompute subgrid-scale orography statistics for the GSL drag directly on the native MPAS mesh (case 7 only)"
possible_values="true or false"/>

<nml_option name="config_gwd_cell_scaling" type="real" default_value="1.0" in_defaults="false"
<nml_option name="config_gwd_cell_scaling" type="real" default_value="2.0" in_defaults="false"
units="-"
description="Scaling factor for the effective grid cell diameter used in computation of GWD static fields"
possible_values="Positive real values"/>
Expand Down Expand Up @@ -463,6 +463,7 @@
<var name="ol2" packages="vertical_stage_in;met_stage_in"/>
<var name="ol3" packages="vertical_stage_in;met_stage_in"/>
<var name="ol4" packages="vertical_stage_in;met_stage_in"/>
<var name="elvmax" packages="vertical_stage_in;met_stage_in"/>
<var name="deriv_two" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
<var name="defc_a" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
<var name="defc_b" packages="gwd_stage_in;vertical_stage_in;met_stage_in"/>
Expand Down Expand Up @@ -568,6 +569,7 @@
<var name="ol2" packages="gwd_stage_out;vertical_stage_out;met_stage_out"/>
<var name="ol3" packages="gwd_stage_out;vertical_stage_out;met_stage_out"/>
<var name="ol4" packages="gwd_stage_out;vertical_stage_out;met_stage_out"/>
<var name="elvmax" packages="gwd_stage_out;vertical_stage_out;met_stage_out"/>
<var name="deriv_two"/>
<var name="defc_a"/>
<var name="defc_b"/>
Expand Down Expand Up @@ -932,6 +934,10 @@
<var name="ol4" type="real" dimensions="nCells" units="unitless"
description="effective orographic length for north-westerly flow"
packages="gwd_stage_out;vertical_stage_out;met_stage_out"/>
<var name="elvmax" type="real" dimensions="nCells" units="m"
description="elevation maximum over a grid cell"
packages="gwd_stage_out;vertical_stage_out;met_stage_out"/>


<!-- GSL GWDO fields - meso-scale -->
<var name="var2dls" type="real" dimensions="nCells" units="m"
Expand Down
92 changes: 63 additions & 29 deletions src/core_init_atmosphere/mpas_init_atm_gwd.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,9 @@ end subroutine read_geogrid
!> con
!> ol{1,2,3,4}
!> oa{1,2,3,4}
!>
!> Added elvmax for blocking, revised oa and ol computations
!> date : 24 August 2025, Songyou Hong (hong@ucar.edu)
!
!-----------------------------------------------------------------------
function compute_gwd_fields(domain) result(iErr)
Expand All @@ -113,7 +116,7 @@ function compute_gwd_fields(domain) result(iErr)

type (mpas_pool_type), pointer :: mesh, state
integer :: iCell, i
real (kind=RKIND) :: dc
real (kind=RKIND) :: dc, hratio
real (kind=RKIND), pointer :: config_gwd_cell_scaling
integer, pointer :: nCells
integer, pointer :: nEdges
Expand Down Expand Up @@ -210,7 +213,7 @@ function compute_gwd_fields(domain) result(iErr)
call mpas_pool_get_array(mesh, 'oa2', oa2)
call mpas_pool_get_array(mesh, 'oa3', oa3)
call mpas_pool_get_array(mesh, 'oa4', oa4)
! call mpas_pool_get_array(mesh, 'elvmax', elvmax)
call mpas_pool_get_array(mesh, 'elvmax', elvmax)
! call mpas_pool_get_array(mesh, 'theta', htheta)
! call mpas_pool_get_array(mesh, 'gamma', hgamma)
! call mpas_pool_get_array(mesh, 'sigma', hsigma)
Expand Down Expand Up @@ -247,6 +250,7 @@ function compute_gwd_fields(domain) result(iErr)
! cell by averaging the distances to each of the neighboring cells
!
dc = 0.0
hratio = 0.0
do i=1,nEdgesOnCell(iCell)
dc = dc + dcEdge(edgesOnCell(i,iCell))
end do
Expand All @@ -255,7 +259,13 @@ function compute_gwd_fields(domain) result(iErr)
dc = dc * sphere_radius
end if
dc = dc * config_gwd_cell_scaling

!
! for effective grid spacing (dc) being too small to compute mountain statistics,
! apply a scale awareness with the ratio of dc / mtn data resolution (sg_delta )
! minimum number of points in a box is assumed at 4 x 4
! it apples to standard deviation (var)
hratio = dc / (sg_delta)
if(hratio < 4.0) dc = sg_delta * 4.0
!
! Cut out a rectangular piece of the global 30-arc-second topography
! data that is centered at the lat/lon (in radians) of the current cell being
Expand All @@ -279,6 +289,9 @@ function compute_gwd_fields(domain) result(iErr)
! subroutines to compute each sub-grid orography statistic
!
var2d(iCell) = get_var(box, box_mean, nx, ny)
if(hratio < 4.0) then
var2d(iCell) = var2d(iCell) * hratio * 0.25
endif
con(iCell) = get_con(box, box_landuse, box_mean, nx, ny)
oa1(iCell) = get_oa1(box, box_mean, nx, ny)
oa2(iCell) = get_oa2(box, box_mean, nx, ny)
Expand All @@ -290,14 +303,14 @@ function compute_gwd_fields(domain) result(iErr)
! in a General Circulation Model. J. Climate, 9, 2698-2717.
hc = 1116.2_RKIND - 0.878_RKIND * var2d(iCell)

ol1(iCell) = get_ol1(box, nx, ny)
ol2(iCell) = get_ol2(box, nx, ny)
ol3(iCell) = get_ol3(box, nx, ny)
ol4(iCell) = get_ol4(box, nx, ny)
ol1(iCell) = get_ol1(box, box_mean, nx, ny)
ol2(iCell) = get_ol2(box, box_mean, nx, ny)
ol3(iCell) = get_ol3(box, box_mean, nx, ny)
ol4(iCell) = get_ol4(box, box_mean, nx, ny)

hlanduse(iCell) = get_dom_landmask(box_landuse, nx, ny) ! get dominant land mask in cell

! elvmax(iCell) = get_elvmax(box, nx, ny)
elvmax(iCell) = get_elvmax(box, nx, ny)
! htheta(iCell) = get_htheta(box, dxm, nx, ny)
! hgamma(iCell) = get_hgamma(box, dxm, nx, ny)
! hsigma(iCell) = get_hsigma(box, dxm, nx, ny)
Expand Down Expand Up @@ -945,9 +958,16 @@ real (kind=RKIND) function get_oa1(box, box_mean, nx, ny)
nu = 0
nd = 0
do j=1,ny
do i=1,nx/2
if (box(i,j) > box_mean) nu = nu + 1
end do
if(mod(nx,2).eq.0.) then
do i=1,nx/2
if (box(i,j) > box_mean) nu = nu + 1
end do
else
do i=1,nx/2 + 1
if (box(i,j) > box_mean) nu = nu + 1
end do
endif

do i=nx/2+1,nx
if (box(i,j) > box_mean) nd = nd + 1
end do
Expand Down Expand Up @@ -987,11 +1007,19 @@ real (kind=RKIND) function get_oa2(box, box_mean, nx, ny)

nu = 0
nd = 0
do j=1,ny/2
do i=1,nx
if(mod(ny,2).eq.0.) then
do j=1,ny/2
do i=1,nx
if (box(i,j) > box_mean) nu = nu + 1
end do
end do
end do
end do
else
do j=1,ny/2 + 1
do i=1,nx
if (box(i,j) > box_mean) nu = nu + 1
end do
end do
endif
do j=ny/2+1,ny
do i=1,nx
if (box(i,j) > box_mean) nd = nd + 1
Expand Down Expand Up @@ -1036,9 +1064,10 @@ real (kind=RKIND) function get_oa3(box, box_mean, nx, ny)
ratio = real(ny,RKIND)/real(nx,RKIND)
do j=1,ny
do i=1,nx
if (nint(real(i,RKIND) * ratio) < (ny - j)) then
if (nint(real(i,RKIND) * ratio) <= (ny - j + 1)) then
if (box(i,j) > box_mean) nu = nu + 1
else
endif
if (nint(real(i,RKIND) * ratio) >= (ny - j + 1)) then
if (box(i,j) > box_mean) nd = nd + 1
end if
end do
Expand Down Expand Up @@ -1082,9 +1111,10 @@ real (kind=RKIND) function get_oa4(box, box_mean, nx, ny)
ratio = real(ny,RKIND)/real(nx,RKIND)
do j=1,ny
do i=1,nx
if (nint(real(i,RKIND) * ratio) < j) then
if (nint(real(i,RKIND) * ratio) <= j) then
if (box(i,j) > box_mean) nu = nu + 1
else
endif
if (nint(real(i,RKIND) * ratio) >= j) then
if (box(i,j) > box_mean) nd = nd + 1
end if
end do
Expand All @@ -1109,10 +1139,11 @@ end function get_oa4
!> \details
!
!-----------------------------------------------------------------------
real (kind=RKIND) function get_ol1(box, nx, ny)
real (kind=RKIND) function get_ol1(box, box_mean, nx, ny)

implicit none

real (kind=RKIND), intent(in) :: box_mean
real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell
integer, intent(in) :: nx, ny

Expand All @@ -1125,7 +1156,7 @@ real (kind=RKIND) function get_ol1(box, nx, ny)

do j=ny/4,3*ny/4
do i=1,nx
if (box(i,j) > hc) nw = nw + 1
if (box(i,j) > box_mean) nw = nw + 1
nt = nt + 1
end do
end do
Expand All @@ -1145,10 +1176,11 @@ end function get_ol1
!> \details
!
!-----------------------------------------------------------------------
real (kind=RKIND) function get_ol2(box, nx, ny)
real (kind=RKIND) function get_ol2(box, box_mean, nx, ny)

implicit none

real (kind=RKIND), intent(in) :: box_mean
real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell
integer, intent(in) :: nx, ny

Expand All @@ -1161,7 +1193,7 @@ real (kind=RKIND) function get_ol2(box, nx, ny)

do j=1,ny
do i=nx/4,3*nx/4
if (box(i,j) > hc) nw = nw + 1
if (box(i,j) > box_mean) nw = nw + 1
nt = nt + 1
end do
end do
Expand All @@ -1181,10 +1213,11 @@ end function get_ol2
!> \details
!
!-----------------------------------------------------------------------
real (kind=RKIND) function get_ol3(box, nx, ny)
real (kind=RKIND) function get_ol3(box, box_mean, nx, ny)

implicit none

real (kind=RKIND), intent(in) :: box_mean
real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell
integer, intent(in) :: nx, ny

Expand All @@ -1197,13 +1230,13 @@ real (kind=RKIND) function get_ol3(box, nx, ny)

do j=1,ny/2
do i=1,nx/2
if (box(i,j) > hc) nw = nw + 1
if (box(i,j) > box_mean) nw = nw + 1
nt = nt + 1
end do
end do
do j=ny/2+1,ny
do i=nx/2+1,nx
if (box(i,j) > hc) nw = nw + 1
if (box(i,j) > box_mean) nw = nw + 1
nt = nt + 1
end do
end do
Expand All @@ -1223,10 +1256,11 @@ end function get_ol3
!> \details
!
!-----------------------------------------------------------------------
real (kind=RKIND) function get_ol4(box, nx, ny)
real (kind=RKIND) function get_ol4(box, box_mean, nx, ny)

implicit none

real (kind=RKIND), intent(in) :: box_mean
real (kind=RKIND), dimension(:,:), pointer, intent(in) :: box ! Subset of topography covering a grid cell
integer, intent(in) :: nx, ny

Expand All @@ -1239,13 +1273,13 @@ real (kind=RKIND) function get_ol4(box, nx, ny)

do j=ny/2+1,ny
do i=1,nx/2
if (box(i,j) > hc) nw = nw + 1
if (box(i,j) > box_mean) nw = nw + 1
nt = nt + 1
end do
end do
do j=1,ny/2
do i=nx/2+1,nx
if (box(i,j) > hc) nw = nw + 1
if (box(i,j) > box_mean) nw = nw + 1
nt = nt + 1
end do
end do
Expand Down