diff --git a/src/core_atmosphere/Externals.cfg b/src/core_atmosphere/Externals.cfg
index 84dc47d1d8..9154ac7d7e 100644
--- a/src/core_atmosphere/Externals.cfg
+++ b/src/core_atmosphere/Externals.cfg
@@ -1,8 +1,8 @@
[MMM-physics]
local_path = ./physics_mmm
protocol = git
-repo_url = https://github.com/NCAR/MMM-physics.git
-tag = 20250616-MPASv8.3
+repo_url = https://github.com/Songyou184/MMM-physics.git
+branch=GWDO_SH
required = True
[GSL_UGWP]
diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index 4281c40bba..8fe33fdaf5 100644
--- a/src/core_atmosphere/Registry.xml
+++ b/src/core_atmosphere/Registry.xml
@@ -406,6 +406,7 @@
+
@@ -471,6 +472,7 @@
+
@@ -521,6 +523,7 @@
+
@@ -1434,10 +1437,13 @@
-
-
+
+
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in;cu_ntiedtke_in;mp_kessler_in;mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in;cu_ntiedtke_in;mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in;cu_ntiedtke_in;mp_kessler_in;mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in;cu_ntiedtke_in;mp_thompson_in;mp_thompson_aers_in;mp_wsm6_in"/>
+
+
+
+
+
+
+ possible_values="`suite',`bl_ysu',`bl_shinhong',`bl_mynn',`off'"/>
+ possible_values="`suite',`bl_kim_gwdo',`bl_ugwp_gwdo',`off'"/>
+
+
+
+
@@ -2651,37 +2688,45 @@
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+
-
+
+
+
+
+ packages="bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_ysu_in;bl_shinhong_in"/>
@@ -2701,10 +2746,6 @@
description="liquid water - liquid water potential temperature covariance"
packages="bl_mynn_in"/>
-
-
@@ -2721,10 +2762,6 @@
description="liquid water potential temperature variance"
packages="bl_mynn_in"/>
-
-
@@ -2792,121 +2829,121 @@
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_ysu_in;bl_shinhong_in"/>
@@ -2928,19 +2965,19 @@
-
-
-
-
-
@@ -3549,27 +3586,27 @@
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+ packages="bl_mynn_in;bl_ysu_in;bl_shinhong_in"/>
+
-
@@ -3893,7 +3932,7 @@
units="kg kg^{-1}"
description="Rain water mixing ratio increment"/>
-
diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F
index f7d04a1f0c..7cfd10b1bf 100644
--- a/src/core_atmosphere/mpas_atm_core.F
+++ b/src/core_atmosphere/mpas_atm_core.F
@@ -372,6 +372,7 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
use mpas_vector_reconstruction
use mpas_stream_manager
use mpas_atm_boundaries, only : mpas_atm_setup_bdy_masks
+ use mpas_constants, only : omega
#ifdef DO_PHYSICS
! use mpas_atmphys_aquaplanet
use mpas_atmphys_control
@@ -405,6 +406,7 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
real (kind=RKIND), dimension(:), pointer :: dvEdge, invDvEdge
real (kind=RKIND), dimension(:), pointer :: dcEdge, invDcEdge
real (kind=RKIND), dimension(:), pointer :: areaTriangle, invAreaTriangle
+ real (kind=RKIND), dimension(:), pointer :: fcell, latcell
integer, pointer :: nCells_ptr, nEdges_ptr, nVertices_ptr, nVertLevels_ptr, nEdgesSolve
integer :: nCells, nEdges, nVertices, nVertLevels
integer :: thread
@@ -444,6 +446,8 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
call mpas_pool_get_array(mesh, 'areaTriangle', areaTriangle)
call mpas_pool_get_array(mesh, 'invAreaTriangle', invAreaTriangle)
+ call mpas_pool_get_array(mesh, 'fCell', fcell)
+ call mpas_pool_get_array(mesh, 'latCell', latcell)
call mpas_pool_get_dimension(mesh, 'nCells', nCells_ptr)
call mpas_pool_get_dimension(mesh, 'nEdges', nEdges_ptr)
@@ -469,6 +473,10 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
invAreaTriangle(iVertex) = 1.0_RKIND / areaTriangle(iVertex)
end do
+ do iCell=1,nCells
+ fcell(iCell) = 2. * omega * sin(latcell(iCell))
+ end do
+
!!!!! End compute inverses
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
diff --git a/src/core_atmosphere/physics/.a.swp b/src/core_atmosphere/physics/.a.swp
new file mode 100644
index 0000000000..d16aed2e43
Binary files /dev/null and b/src/core_atmosphere/physics/.a.swp differ
diff --git a/src/core_atmosphere/physics/mpas_atmphys_control.F b/src/core_atmosphere/physics/mpas_atmphys_control.F
index b3162019e5..933e28d282 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_control.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_control.F
@@ -83,6 +83,10 @@ module mpas_atmphys_control
! * in the mesoscale_reference suite, replaced the MM5 surface layer scheme with the MM5 revised surface layer
! scheme as the default option for config_sfclayer_scheme.
! Laura D. Fowler (laura@ucar.edu) / 2024-06-18.
+! * renamed "bl_ysu_gwdo" to "bl_kim_gwdo"
+! Songyou Hong (hong@ucar.edu) / 2025-08-24.
+! * added the option bl_shinhong
+! Songyou Hong (hong@ucar.edu) / 2025-11-27.
contains
@@ -133,7 +137,7 @@ subroutine physics_namelist_check(configs)
if (trim(config_microp_scheme) == 'suite') config_microp_scheme = 'mp_wsm6'
if (trim(config_convection_scheme) == 'suite') config_convection_scheme = 'cu_ntiedtke'
if (trim(config_pbl_scheme) == 'suite') config_pbl_scheme = 'bl_ysu'
- if (trim(config_gwdo_scheme) == 'suite') config_gwdo_scheme = 'bl_ysu_gwdo'
+ if (trim(config_gwdo_scheme) == 'suite') config_gwdo_scheme = 'bl_kim_gwdo'
if (trim(config_radt_lw_scheme) == 'suite') config_radt_lw_scheme = 'rrtmg_lw'
if (trim(config_radt_sw_scheme) == 'suite') config_radt_sw_scheme = 'rrtmg_sw'
if (trim(config_radt_cld_scheme) == 'suite') config_radt_cld_scheme = 'cld_fraction'
@@ -145,7 +149,7 @@ subroutine physics_namelist_check(configs)
if (trim(config_microp_scheme) == 'suite') config_microp_scheme = 'mp_thompson'
if (trim(config_convection_scheme) == 'suite') config_convection_scheme = 'cu_grell_freitas'
if (trim(config_pbl_scheme) == 'suite') config_pbl_scheme = 'bl_mynn'
- if (trim(config_gwdo_scheme) == 'suite') config_gwdo_scheme = 'bl_ysu_gwdo'
+ if (trim(config_gwdo_scheme) == 'suite') config_gwdo_scheme = 'bl_kim_gwdo'
if (trim(config_radt_lw_scheme) == 'suite') config_radt_lw_scheme = 'rrtmg_lw'
if (trim(config_radt_sw_scheme) == 'suite') config_radt_sw_scheme = 'rrtmg_sw'
if (trim(config_radt_cld_scheme) == 'suite') config_radt_cld_scheme = 'cld_fraction'
@@ -201,6 +205,7 @@ subroutine physics_namelist_check(configs)
!pbl scheme:
if(.not. (config_pbl_scheme .eq. 'off' .or. &
config_pbl_scheme .eq. 'bl_mynn' .or. &
+ config_pbl_scheme .eq. 'bl_shinhong' .or. &
config_pbl_scheme .eq. 'bl_ysu')) then
write(mpas_err_message,'(A,A20)') 'illegal value for pbl_scheme: ', &
@@ -211,7 +216,7 @@ subroutine physics_namelist_check(configs)
!gravity wave drag over orography scheme:
if(.not. (config_gwdo_scheme .eq. 'off' .or. &
- config_gwdo_scheme .eq. 'bl_ysu_gwdo' .or. &
+ config_gwdo_scheme .eq. 'bl_kim_gwdo' .or. &
config_gwdo_scheme .eq. 'bl_ugwp_gwdo')) then
write(mpas_err_message,'(A,A20)') 'illegal value for gwdo_scheme: ', &
@@ -279,10 +284,10 @@ subroutine physics_namelist_check(configs)
else
if(config_pbl_scheme == 'bl_mynn') then
config_sfclayer_scheme = 'sf_mynn'
- elseif(config_pbl_scheme == 'bl_ysu') then
+ elseif(config_pbl_scheme == 'bl_shinhong' .or. config_pbl_scheme == 'bl_ysu') then
if(config_sfclayer_scheme /= 'sf_monin_obukhov' .and. &
config_sfclayer_scheme /= 'sf_monin_obukhov_rev') then
- write(mpas_err_message,'(A,A20)') 'wrong choice for surface layer scheme with YSU PBL: ', &
+ write(mpas_err_message,'(A,A20)') 'wrong choice for surface layer scheme with SHINHONG/YSU PBL: ', &
trim(config_sfclayer_scheme)
call physics_error_fatal(mpas_err_message)
endif
@@ -481,7 +486,7 @@ subroutine physics_compatibility_check(dminfo, blockList, streamManager, ierr)
call mpas_pool_get_config(blocklist % configs, 'config_gwdo_scheme', gwdo_scheme)
- if (trim(gwdo_scheme) == 'bl_ysu_gwdo') then
+ if (trim(gwdo_scheme) == 'bl_kim_gwdo') then
maxvar2d_local = -huge(maxvar2d_local)
block => blockList
do while (associated(block))
diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver.F b/src/core_atmosphere/physics/mpas_atmphys_driver.F
index 8e31672657..134d8ee904 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_driver.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_driver.F
@@ -143,6 +143,7 @@ subroutine physics_driver(domain,itimestep,xtime_s)
config_sfclayer_scheme
logical, pointer:: config_oml1d
+ logical, pointer:: config_gwdo_nonhyd, config_kim_tofd
real(kind=RKIND),pointer:: config_bucket_radt
!local variables:
diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_gwdo.F b/src/core_atmosphere/physics/mpas_atmphys_driver_gwdo.F
index a96ba7bf2a..c95b4402a7 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_driver_gwdo.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_driver_gwdo.F
@@ -10,22 +10,17 @@ module mpas_atmphys_driver_gwdo
use mpas_kind_types
use mpas_pool_routines
use mpas_timer,only: mpas_timer_start,mpas_timer_stop
-
use mpas_atmphys_constants
use mpas_atmphys_vars
use mpas_atmphys_manager,only: curr_julday
-
!wrf physics:
use module_bl_gwdo
use module_bl_ugwp_gwdo
-
implicit none
private
public:: allocate_gwdo, &
deallocate_gwdo, &
driver_gwdo
-
-
!MPAS driver for parameterization of gravity wave drag over orography.
!Laura D. Fowler (send comments to laura@ucar.edu).
!2013-05-01.
@@ -66,29 +61,22 @@ module mpas_atmphys_driver_gwdo
! Laura D. Fowler (laura@ucar.edu) / 2023-05-15.
! * added the NOAA UFS unified gravity wave drag scheme
! Michael D. Toy (michael.toy@noaa.gov) / 2024-10-21
-
-
+! * Revised the ysu_gwdo follwing the studies (hong et al. 2025, koo and hong 2025), and renamed to kim_gwdo
+! Songyou Hong (hong@ucar.edu) / 2025-11-27
contains
-
-
!=================================================================================================================
subroutine allocate_gwdo(configs)
!=================================================================================================================
-
!input arguments:
type(mpas_pool_type),intent(in):: configs
-
!local variables:
character(len=StrKIND),pointer:: gwdo_scheme
logical,pointer:: ugwp_diags,ngw_scheme
-
call mpas_pool_get_config(configs,'config_gwdo_scheme',gwdo_scheme)
call mpas_pool_get_config(configs,'config_ugwp_diags',ugwp_diags)
call mpas_pool_get_config(configs,'config_ngw_scheme',ngw_scheme)
-
if(.not.allocated(cosa_p) ) allocate(cosa_p(ims:ime,jms:jme) )
if(.not.allocated(sina_p) ) allocate(sina_p(ims:ime,jms:jme) )
-
if(.not.allocated(dx_p) ) allocate(dx_p(ims:ime,jms:jme) )
if(.not.allocated(kpbl_p )) allocate(kpbl_p(ims:ime,jms:jme) )
if(.not.allocated(dusfcg_p)) allocate(dusfcg_p(ims:ime,jms:jme))
@@ -98,10 +86,8 @@ subroutine allocate_gwdo(configs)
if(.not.allocated(rublten_p)) allocate(rublten_p(ims:ime,kms:kme,jms:jme))
if(.not.allocated(rvblten_p)) allocate(rvblten_p(ims:ime,kms:kme,jms:jme))
if(.not.allocated(rthblten_p)) allocate(rthblten_p(ims:ime,kms:kme,jms:jme))
-
gwdo_select: select case (trim(gwdo_scheme))
-
- case("bl_ysu_gwdo")
+ case("bl_kim_gwdo")
if(.not.allocated(var2d_p) ) allocate(var2d_p(ims:ime,jms:jme) )
if(.not.allocated(con_p) ) allocate(con_p(ims:ime,jms:jme) )
if(.not.allocated(oa1_p) ) allocate(oa1_p(ims:ime,jms:jme) )
@@ -112,7 +98,7 @@ subroutine allocate_gwdo(configs)
if(.not.allocated(ol2_p) ) allocate(ol2_p(ims:ime,jms:jme) )
if(.not.allocated(ol3_p) ) allocate(ol3_p(ims:ime,jms:jme) )
if(.not.allocated(ol4_p) ) allocate(ol4_p(ims:ime,jms:jme) )
-
+ if(.not.allocated(elvmax_p)) allocate(elvmax_p(ims:ime,jms:jme))
case("bl_ugwp_gwdo")
if(.not.allocated(var2dls_p) ) allocate(var2dls_p(ims:ime,jms:jme) )
if(.not.allocated(conls_p) ) allocate(conls_p(ims:ime,jms:jme) )
@@ -169,31 +155,22 @@ subroutine allocate_gwdo(configs)
if(.not.allocated(ddy_j1tau_p)) allocate(ddy_j1tau_p(ims:ime,jms:jme))
if(.not.allocated(ddy_j2tau_p)) allocate(ddy_j2tau_p(ims:ime,jms:jme))
endif
-
case default
-
end select gwdo_select
-
end subroutine allocate_gwdo
-
!=================================================================================================================
subroutine deallocate_gwdo(configs)
!=================================================================================================================
-
!input arguments:
type(mpas_pool_type),intent(in):: configs
-
!local variables:
character(len=StrKIND),pointer:: gwdo_scheme
logical,pointer:: ugwp_diags,ngw_scheme
-
call mpas_pool_get_config(configs,'config_gwdo_scheme',gwdo_scheme)
call mpas_pool_get_config(configs,'config_ugwp_diags',ugwp_diags)
call mpas_pool_get_config(configs,'config_ngw_scheme',ngw_scheme)
-
if(allocated(cosa_p) ) deallocate(cosa_p )
if(allocated(sina_p) ) deallocate(sina_p )
-
if(allocated(dx_p) ) deallocate(dx_p )
if(allocated(kpbl_p) ) deallocate(kpbl_p )
if(allocated(dusfcg_p)) deallocate(dusfcg_p)
@@ -203,10 +180,8 @@ subroutine deallocate_gwdo(configs)
if(allocated(rublten_p)) deallocate(rublten_p)
if(allocated(rvblten_p)) deallocate(rvblten_p)
if(allocated(rthblten_p)) deallocate(rthblten_p)
-
gwdo_select: select case (trim(gwdo_scheme))
-
- case("bl_ysu_gwdo")
+ case("bl_kim_gwdo")
if(allocated(var2d_p) ) deallocate(var2d_p )
if(allocated(con_p) ) deallocate(con_p )
if(allocated(oa1_p) ) deallocate(oa1_p )
@@ -217,7 +192,7 @@ subroutine deallocate_gwdo(configs)
if(allocated(ol2_p) ) deallocate(ol2_p )
if(allocated(ol3_p) ) deallocate(ol3_p )
if(allocated(ol4_p) ) deallocate(ol4_p )
-
+ if(allocated(elvmax_p)) deallocate(elvmax_p)
case("bl_ugwp_gwdo")
if(allocated(var2dls_p) ) deallocate(var2dls_p )
if(allocated(conls_p) ) deallocate(conls_p )
@@ -274,17 +249,12 @@ subroutine deallocate_gwdo(configs)
if(allocated(ddy_j1tau_p)) deallocate(ddy_j1tau_p)
if(allocated(ddy_j2tau_p)) deallocate(ddy_j2tau_p)
endif
-
case default
-
end select gwdo_select
-
end subroutine deallocate_gwdo
-
!=================================================================================================================
subroutine gwdo_from_MPAS(configs,mesh,sfc_input,ngw_input,diag_physics,tend_physics,its,ite)
!=================================================================================================================
-
!input arguments:
type(mpas_pool_type),intent(in):: configs
type(mpas_pool_type),intent(in):: mesh
@@ -292,22 +262,19 @@ subroutine gwdo_from_MPAS(configs,mesh,sfc_input,ngw_input,diag_physics,tend_phy
type(mpas_pool_type),intent(in):: ngw_input
type(mpas_pool_type),intent(in):: diag_physics
type(mpas_pool_type),intent(in):: tend_physics
-
integer,intent(in):: its,ite
-
!local variables:
integer:: i,k,j
character(len=StrKIND),pointer:: gwdo_scheme
character(len=StrKIND),pointer:: convection_scheme,microp_scheme
logical,pointer:: ugwp_diags,ngw_scheme
real(kind=RKIND),parameter :: rad2deg = 180./3.1415926
-
!local pointers:
integer,dimension(:),pointer:: kpbl
integer,dimension(:),pointer:: jindx1_tau,jindx2_tau
real(kind=RKIND),pointer:: len_disp
real(kind=RKIND),dimension(:),pointer :: meshDensity
- real(kind=RKIND),dimension(:),pointer :: oa1,oa2,oa3,oa4,ol1,ol2,ol3,ol4,con,var2d
+ real(kind=RKIND),dimension(:),pointer :: oa1,oa2,oa3,oa4,ol1,ol2,ol3,ol4,con,var2d,elvmax
real(kind=RKIND),dimension(:),pointer :: oa1ls,oa2ls,oa3ls,oa4ls,ol1ls,ol2ls, &
ol3ls,ol4ls,conls,var2dls
real(kind=RKIND),dimension(:),pointer :: oa1ss,oa2ss,oa3ss,oa4ss,ol1ss,ol2ss, &
@@ -322,9 +289,7 @@ subroutine gwdo_from_MPAS(configs,mesh,sfc_input,ngw_input,diag_physics,tend_phy
real(kind=RKIND),dimension(:,:),pointer:: dtaux3d_ls,dtauy3d_ls,dtaux3d_bl,dtauy3d_bl, &
dtaux3d_ss,dtauy3d_ss,dtaux3d_fd,dtauy3d_fd
real(kind=RKIND),dimension(:,:),pointer:: dudt_ngw,dvdt_ngw,dtdt_ngw
-
!-----------------------------------------------------------------------------------------------------------------
-
call mpas_pool_get_config(configs,'config_len_disp',len_disp)
call mpas_pool_get_config(configs,'config_gwdo_scheme',gwdo_scheme)
call mpas_pool_get_config(configs,'config_ugwp_diags',ugwp_diags)
@@ -332,12 +297,10 @@ subroutine gwdo_from_MPAS(configs,mesh,sfc_input,ngw_input,diag_physics,tend_phy
call mpas_pool_get_config(configs,'config_convection_scheme',convection_scheme)
call mpas_pool_get_config(configs,'config_microp_scheme',microp_scheme)
call mpas_pool_get_array(mesh,'meshDensity',meshDensity)
-
-
gwdo_select: select case (trim(gwdo_scheme))
-
- case("bl_ysu_gwdo")
+ case("bl_kim_gwdo")
call mpas_pool_get_array(sfc_input,'var2d',var2d)
+ call mpas_pool_get_array(sfc_input,'elvmax',elvmax)
call mpas_pool_get_array(sfc_input,'con' ,con )
call mpas_pool_get_array(sfc_input,'oa1' ,oa1 )
call mpas_pool_get_array(sfc_input,'oa2' ,oa2 )
@@ -359,9 +322,9 @@ subroutine gwdo_from_MPAS(configs,mesh,sfc_input,ngw_input,diag_physics,tend_phy
ol2_p(i,j) = ol2(i)
ol3_p(i,j) = ol3(i)
ol4_p(i,j) = ol4(i)
+ elvmax_p(i,j)= elvmax(i)
enddo
enddo
-
case("bl_ugwp_gwdo")
call mpas_pool_get_array(sfc_input,'var2dls',var2dls)
call mpas_pool_get_array(sfc_input,'conls' ,conls )
@@ -509,14 +472,9 @@ subroutine gwdo_from_MPAS(configs,mesh,sfc_input,ngw_input,diag_physics,tend_phy
enddo
enddo
endif
-
endif
-
case default
-
end select gwdo_select
-
-
call mpas_pool_get_array(diag_physics,'kpbl' ,kpbl )
call mpas_pool_get_array(diag_physics,'dusfcg' ,dusfcg )
call mpas_pool_get_array(diag_physics,'dvsfcg' ,dvsfcg )
@@ -525,7 +483,6 @@ subroutine gwdo_from_MPAS(configs,mesh,sfc_input,ngw_input,diag_physics,tend_phy
call mpas_pool_get_array(tend_physics,'rublten' ,rublten )
call mpas_pool_get_array(tend_physics,'rvblten' ,rvblten )
call mpas_pool_get_array(tend_physics,'rthblten',rthblten)
-
do j = jts,jte
do i = its,ite
sina_p(i,j) = 0._RKIND
@@ -536,7 +493,6 @@ subroutine gwdo_from_MPAS(configs,mesh,sfc_input,ngw_input,diag_physics,tend_phy
dvsfcg_p(i,j) = dvsfcg(i)
enddo
enddo
-
do j = jts,jte
do k = kts,kte
do i = its,ite
@@ -548,31 +504,24 @@ subroutine gwdo_from_MPAS(configs,mesh,sfc_input,ngw_input,diag_physics,tend_phy
enddo
enddo
enddo
-
end subroutine gwdo_from_MPAS
-
!=================================================================================================================
subroutine gwdo_to_MPAS(configs,diag_physics,tend_physics,its,ite)
!=================================================================================================================
-
!input arguments:
integer,intent(in):: its,ite
type(mpas_pool_type),intent(in):: configs
-
!inout arguments:
type(mpas_pool_type),intent(inout):: diag_physics
type(mpas_pool_type),intent(inout):: tend_physics
-
!local variables:
integer:: i,k,j
character(len=StrKIND),pointer:: gwdo_scheme
logical,pointer:: ugwp_diags,ngw_scheme
-
!local pointers:
real(kind=RKIND),dimension(:),pointer :: dusfcg,dvsfcg
real(kind=RKIND),dimension(:,:),pointer:: dtaux3d,dtauy3d,rubldiff,rvbldiff,rublten,rvblten
real(kind=RKIND),dimension(:,:),pointer:: rthblten
-
real(kind=RKIND),dimension(:),pointer :: oa1ls,oa2ls,oa3ls,oa4ls,ol1ls,ol2ls, &
ol3ls,ol4ls,conls,var2dls
real(kind=RKIND),dimension(:),pointer :: oa1ss,oa2ss,oa3ss,oa4ss,ol1ss,ol2ss, &
@@ -582,9 +531,7 @@ subroutine gwdo_to_MPAS(configs,diag_physics,tend_physics,its,ite)
real(kind=RKIND),dimension(:,:),pointer:: dtaux3d_ls,dtauy3d_ls,dtaux3d_bl,dtauy3d_bl, &
dtaux3d_ss,dtauy3d_ss,dtaux3d_fd,dtauy3d_fd
real(kind=RKIND),dimension(:,:),pointer:: dudt_ngw,dvdt_ngw,dtdt_ngw
-
!-----------------------------------------------------------------------------------------------------------------
-
call mpas_pool_get_config(configs,'config_gwdo_scheme',gwdo_scheme)
call mpas_pool_get_config(configs,'config_ugwp_diags',ugwp_diags)
call mpas_pool_get_config(configs,'config_ngw_scheme',ngw_scheme)
@@ -597,10 +544,7 @@ subroutine gwdo_to_MPAS(configs,diag_physics,tend_physics,its,ite)
call mpas_pool_get_array(tend_physics,'rublten' ,rublten )
call mpas_pool_get_array(tend_physics,'rvblten' ,rvblten )
call mpas_pool_get_array(tend_physics,'rthblten',rthblten)
-
-
gwdo_select: select case (trim(gwdo_scheme))
-
case("bl_ugwp_gwdo")
if (ugwp_diags) then
call mpas_pool_get_array(diag_physics,'dusfc_ls' ,dusfc_ls )
@@ -660,18 +604,14 @@ subroutine gwdo_to_MPAS(configs,diag_physics,tend_physics,its,ite)
enddo
endif
endif
-
case default
-
end select gwdo_select
-
do j = jts,jte
do i = its,ite
dusfcg(i) = dusfcg_p(i,j)
dvsfcg(i) = dvsfcg_p(i,j)
enddo
enddo
-
do j = jts,jte
do k = kts,kte
do i = its,ite
@@ -685,55 +625,45 @@ subroutine gwdo_to_MPAS(configs,diag_physics,tend_physics,its,ite)
enddo
enddo
enddo
-
end subroutine gwdo_to_MPAS
-
!=================================================================================================================
subroutine driver_gwdo(itimestep,configs,mesh,sfc_input,ngw_input,diag_physics,tend_physics,its,ite)
!=================================================================================================================
-
!input arguments:
type(mpas_pool_type),intent(in):: configs
type(mpas_pool_type),intent(in):: mesh
type(mpas_pool_type),intent(in):: sfc_input
-
integer,intent(in):: its,ite
integer,intent(in):: itimestep
-
!inout arguments:
type(mpas_pool_type),intent(inout):: ngw_input
type(mpas_pool_type),intent(inout):: diag_physics
type(mpas_pool_type),intent(inout):: tend_physics
-
!local variables:
character(len=StrKIND),pointer:: gwdo_scheme
logical,pointer:: ugwp_diags,ngw_scheme
+ logical,pointer:: if_nonhyd
integer,pointer:: ntau_d1y_ptr,ntau_d2t_ptr
real(kind=RKIND),dimension(:),pointer :: days_limb_ptr
real(kind=RKIND),dimension(:,:),pointer:: tau_limb_ptr
integer:: ntau_d1y,ntau_d2t
real(kind=RKIND),dimension(:),allocatable:: days_limb
real(kind=RKIND),dimension(:,:),allocatable:: tau_limb
-
integer:: i
real(kind=RKIND),dimension(:),allocatable:: dx_max
-
+ real(kind=RKIND),pointer :: dx_factor
!CCPP-compliant flags:
character(len=StrKIND):: errmsg
integer:: errflg
-
!-----------------------------------------------------------------------------------------------------------------
!call mpas_log_write('')
!call mpas_log_write('--- enter subroutine driver_gwdo:')
-
!initialization of CCPP-compliant flags:
errmsg = ' '
errflg = 0
-
call mpas_pool_get_config(configs,'config_gwdo_scheme',gwdo_scheme)
call mpas_pool_get_config(configs,'config_ugwp_diags',ugwp_diags)
call mpas_pool_get_config(configs,'config_ngw_scheme',ngw_scheme)
-
! Call up variables needed for NGW scheme
if (ngw_scheme) then
call mpas_pool_get_dimension(mesh,'lat',ntau_d1y_ptr)
@@ -747,14 +677,12 @@ subroutine driver_gwdo(itimestep,configs,mesh,sfc_input,ngw_input,diag_physics,t
days_limb(:) = days_limb_ptr(:)
tau_limb (:,:) = tau_limb_ptr(:,:)
endif
-
-
!copy MPAS arrays to local arrays:
call gwdo_from_MPAS(configs,mesh,sfc_input,ngw_input,diag_physics,tend_physics,its,ite)
-
gwdo_select: select case (trim(gwdo_scheme))
-
- case("bl_ysu_gwdo")
+ case("bl_kim_gwdo")
+ call mpas_pool_get_config(configs,'config_gwdo_factor',dx_factor)
+ call mpas_pool_get_config(configs,'config_gwdo_nonhyd',if_nonhyd)
call mpas_timer_start('bl_gwdo')
call gwdo ( &
p3d = pres_hydd_p , p3di = pres2_hydd_p , pi3d = pi_p , &
@@ -768,14 +696,14 @@ subroutine driver_gwdo(itimestep,configs,mesh,sfc_input,ngw_input,diag_physics,t
var2d = var2d_p , oc12d = con_p , oa2d1 = oa1_p , &
oa2d2 = oa2_p , oa2d3 = oa3_p , oa2d4 = oa4_p , &
ol2d1 = ol1_p , ol2d2 = ol2_p , ol2d3 = ol3_p , &
- ol2d4 = ol4_p , sina = sina_p , cosa = cosa_p , &
- errmsg = errmsg , errflg = errflg , &
+ ol2d4 = ol4_p , elvmax = elvmax_p , sina = sina_p , &
+ cosa = cosa_p , errmsg = errmsg , errflg = errflg , &
+ dx_factor = dx_factor , if_nonhyd = if_nonhyd, &
ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde , &
ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme , &
its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte &
)
call mpas_timer_stop('bl_gwdo')
-
case("bl_ugwp_gwdo")
call mpas_timer_start('bl_ugwp_gwdo')
call gwdo_ugwp ( &
@@ -820,19 +748,13 @@ subroutine driver_gwdo(itimestep,configs,mesh,sfc_input,ngw_input,diag_physics,t
if(allocated(tau_limb) ) deallocate(tau_limb )
endif
call mpas_timer_stop('bl_ugwp_gwdo')
-
case default
-
end select gwdo_select
-
!copy local arrays to MPAS grid:
call gwdo_to_MPAS(configs,diag_physics,tend_physics,its,ite)
-
!call mpas_log_write('--- end subroutine driver_gwdo.')
!call mpas_log_write('')
-
end subroutine driver_gwdo
-
!=================================================================================================================
end module mpas_atmphys_driver_gwdo
!=================================================================================================================
diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_pbl.F b/src/core_atmosphere/physics/mpas_atmphys_driver_pbl.F
index 72a411aeba..c73c70eb2b 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_driver_pbl.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_driver_pbl.F
@@ -16,6 +16,7 @@ module mpas_atmphys_driver_pbl
use bl_mynn,only: bl_mynn_init
use module_bl_mynn,only: mynn_bl_driver
+ use module_bl_shinhong
use module_bl_ysu
implicit none
@@ -39,6 +40,7 @@ module mpas_atmphys_driver_pbl
!
! WRF physics called from driver_pbl:
! -----------------------------------
+! * module_bl_shinhong : SHINHONG PBL scheme.
! * module_bl_ysu : YSU PBL scheme.
!
! add-ons and modifications to sourcecode:
@@ -78,6 +80,8 @@ module mpas_atmphys_driver_pbl
! Laura D. Fowler (laura@ucar.edu) / 2024-02-14.
! * updated the MYNN PBL scheme to the sourcecode from WRF version 4.6.
! Laura D. Fowler (laura@ucar.edu) / 2024-02.15.
+! * added the option SHINHONG PBL scheme, incorporating revisions and additions from WRF version 4.7
+! Songyou Hong (hong@ucar.edu) / 2025-11.27.
contains
@@ -102,6 +106,7 @@ subroutine allocate_pbl(configs)
if(.not.allocated(ust_p) ) allocate(ust_p(ims:ime,jms:jme) )
if(.not.allocated(wspd_p) ) allocate(wspd_p(ims:ime,jms:jme) )
if(.not.allocated(xland_p)) allocate(xland_p(ims:ime,jms:jme))
+ if(.not.allocated(fcell_p)) allocate(fcell_p(ims:ime,jms:jme))
if(.not.allocated(hpbl_p) ) allocate(hpbl_p(ims:ime,jms:jme) )
if(.not.allocated(kpbl_p) ) allocate(kpbl_p(ims:ime,jms:jme) )
if(.not.allocated(znt_p) ) allocate(znt_p(ims:ime,jms:jme) )
@@ -125,6 +130,21 @@ subroutine allocate_pbl(configs)
pbl_select: select case (trim(pbl_scheme))
+ case("bl_shinhong")
+ !from surface-layer model:
+ if(.not.allocated(br_p) ) allocate(br_p(ims:ime,jms:jme) )
+ if(.not.allocated(dx_p) ) allocate(dx_p(ims:ime,jms:jme) )
+ if(.not.allocated(ctopo_p) ) allocate(ctopo_p(ims:ime,jms:jme) )
+ if(.not.allocated(ctopo2_p)) allocate(ctopo2_p(ims:ime,jms:jme) )
+ if(.not.allocated(delta_p) ) allocate(delta_p(ims:ime,jms:jme) )
+ if(.not.allocated(psih_p) ) allocate(psih_p(ims:ime,jms:jme) )
+ if(.not.allocated(psim_p) ) allocate(psim_p(ims:ime,jms:jme) )
+ if(.not.allocated(u10_p) ) allocate(u10_p(ims:ime,jms:jme) )
+ if(.not.allocated(v10_p) ) allocate(v10_p(ims:ime,jms:jme) )
+ if(.not.allocated(exch_p) ) allocate(exch_p(ims:ime,kms:kme,jms:jme))
+ if(.not.allocated(wstar_p) ) allocate(wstar_p(ims:ime,jms:jme) )
+ if(.not.allocated(elpbl_p) ) allocate(elpbl_p(ims:ime,kms:kme,jms:jme) )
+ if(.not.allocated(tkepbl_p) ) allocate(tkepbl_p(ims:ime,kms:kme,jms:jme) )
case("bl_ysu")
!from surface-layer model:
if(.not.allocated(br_p) ) allocate(br_p(ims:ime,jms:jme) )
@@ -211,6 +231,7 @@ subroutine deallocate_pbl(configs)
if(allocated(ust_p) ) deallocate(ust_p )
if(allocated(wspd_p) ) deallocate(wspd_p )
if(allocated(xland_p)) deallocate(xland_p)
+ if(allocated(fcell_p)) deallocate(fcell_p)
if(allocated(hpbl_p) ) deallocate(hpbl_p )
if(allocated(kpbl_p) ) deallocate(kpbl_p )
if(allocated(znt_p) ) deallocate(znt_p )
@@ -234,6 +255,21 @@ subroutine deallocate_pbl(configs)
pbl_select: select case (trim(pbl_scheme))
+ case("bl_shinhong")
+ !from surface-layer model:
+ if(allocated(br_p) ) deallocate(br_p )
+ if(allocated(dx_p) ) deallocate(dx_p )
+ if(allocated(ctopo_p) ) deallocate(ctopo_p )
+ if(allocated(ctopo2_p)) deallocate(ctopo2_p)
+ if(allocated(delta_p) ) deallocate(delta_p )
+ if(allocated(psih_p) ) deallocate(psih_p )
+ if(allocated(psim_p) ) deallocate(psim_p )
+ if(allocated(u10_p) ) deallocate(u10_p )
+ if(allocated(v10_p) ) deallocate(v10_p )
+ if(allocated(exch_p) ) deallocate(exch_p )
+ if(allocated(wstar_p) ) deallocate(wstar_p )
+ if(allocated(elpbl_p) ) deallocate(elpbl_p )
+ if(allocated(tkepbl_p) ) deallocate(tkepbl_p )
case("bl_ysu")
!from surface-layer model:
if(allocated(br_p) ) deallocate(br_p )
@@ -321,13 +357,17 @@ subroutine pbl_from_MPAS(configs,mesh,sfc_input,diag_physics,tend_physics,its,it
!local pointers:
character(len=StrKIND),pointer:: pbl_scheme
- real(kind=RKIND),dimension(:),pointer:: hfx,hpbl,qfx,ust,wspd,xland,znt
+ real(kind=RKIND),dimension(:),pointer:: hfx,hpbl,qfx,ust,wspd,xland,znt,fcell
real(kind=RKIND),dimension(:),pointer:: delta,wstar
!local pointers for YSU scheme:
logical,pointer:: config_ysu_pblmix
real(kind=RKIND),dimension(:),pointer:: br,fh,fm,u10,v10
real(kind=RKIND),dimension(:,:),pointer:: rthratenlw,rthratensw
+!local pointers for SHINHONG scheme:
+ logical,pointer:: config_shinhong_nonlocal_flux
+ logical,pointer:: config_shinhong_scu_mixing
+ logical,pointer:: config_shinhong_ke_dissipation
!local pointers for MYNN scheme:
real(kind=RKIND),pointer:: len_disp
@@ -353,6 +393,7 @@ subroutine pbl_from_MPAS(configs,mesh,sfc_input,diag_physics,tend_physics,its,it
call mpas_pool_get_array(tend_physics,'rthratensw',rthratensw)
call mpas_pool_get_array(sfc_input,'xland',xland)
+ call mpas_pool_get_array(mesh,'fCell',fcell)
do j = jts,jte
do i = its,ite
@@ -363,6 +404,7 @@ subroutine pbl_from_MPAS(configs,mesh,sfc_input,diag_physics,tend_physics,its,it
ust_p(i,j) = ust(i)
wspd_p(i,j) = wspd(i)
xland_p(i,j) = xland(i)
+ fcell_p(i,j) = fcell(i)
kpbl_p(i,j) = 1
znt_p(i,j) = znt(i)
!... ocean currents are set to zero:
@@ -378,6 +420,56 @@ subroutine pbl_from_MPAS(configs,mesh,sfc_input,diag_physics,tend_physics,its,it
pbl_select: select case (trim(pbl_scheme))
+ case("bl_shinhong")
+ call mpas_pool_get_config(configs,'config_shinhong_nonlocal_flux',config_shinhong_nonlocal_flux)
+ call mpas_pool_get_config(configs,'config_shinhong_scu_mixing',config_shinhong_scu_mixing)
+ call mpas_pool_get_config(configs,'config_shinhong_ke_dissipation',config_shinhong_ke_dissipation)
+ call mpas_pool_get_config(configs,'config_len_disp',len_disp)
+
+ call mpas_pool_get_array(mesh,'meshDensity',meshDensity)
+ call mpas_pool_get_array(diag_physics,'br' ,br )
+ call mpas_pool_get_array(diag_physics,'delta',delta)
+ call mpas_pool_get_array(diag_physics,'fm' ,fm )
+ call mpas_pool_get_array(diag_physics,'fh' ,fh )
+ call mpas_pool_get_array(diag_physics,'u10' ,u10 )
+ call mpas_pool_get_array(diag_physics,'v10' ,v10 )
+ call mpas_pool_get_array(diag_physics,'wstar',wstar)
+
+ call mpas_pool_get_array(diag_physics,'el_pbl' ,el_pbl )
+ call mpas_pool_get_array(diag_physics,'tke_pbl' ,tke_pbl )
+
+ do j = jts,jte
+ do i = its,ite
+ dx_p(i,j) = len_disp / meshDensity(i)**0.25
+ enddo
+ enddo
+
+ do j = jts,jte
+ do i = its,ite
+ !from surface-layer model:
+ br_p(i,j) = br(i)
+ psim_p(i,j) = fm(i)
+ psih_p(i,j) = fh(i)
+ u10_p(i,j) = u10(i)
+ v10_p(i,j) = v10(i)
+ delta_p(i,j) = delta(i)
+ wstar_p(i,j) = wstar(i)
+ !initialization for YSU/SHINHONG PBL scheme:
+ ctopo_p(i,j) = 1._RKIND
+ ctopo2_p(i,j) = 1._RKIND
+ enddo
+ enddo
+
+ do j = jts,jte
+ do k = kts,kte
+ do i = its,ite
+ exch_p(i,k,j) = 0._RKIND
+ tkepbl_p(i,k,j) = tke_pbl(k,i)
+ elpbl_p(i,k,j) = el_pbl(k,i)
+ enddo
+ enddo
+ enddo
+
case("bl_ysu")
call mpas_pool_get_config(configs,'config_ysu_pblmix',config_ysu_pblmix)
@@ -608,6 +700,27 @@ subroutine pbl_to_MPAS(configs,diag_physics,tend_physics,its,ite)
pbl_select: select case (trim(pbl_scheme))
+ case("bl_shinhong")
+ call mpas_pool_get_array(diag_physics,'delta',delta )
+ call mpas_pool_get_array(diag_physics,'wstar' ,wstar )
+ call mpas_pool_get_array(diag_physics,'exch_h',exch_h)
+ call mpas_pool_get_array(diag_physics,'tke_pbl',tke_pbl)
+ call mpas_pool_get_array(diag_physics,'el_pbl',el_pbl)
+
+ do j = jts,jte
+ do i = its,ite
+ delta(i) = delta_p(i,j)
+ wstar(i) = wstar_p(i,j)
+ enddo
+ do k = kts,kte
+ do i = its,ite
+ exch_h(k,i) = exch_p(i,k,j)
+ tke_pbl(k,i) = tkepbl_p(i,k,j)
+ el_pbl(k,i) = elpbl_p(i,k,j)
+ enddo
+ enddo
+ enddo
+
case("bl_ysu")
call mpas_pool_get_array(diag_physics,'delta',delta )
call mpas_pool_get_array(diag_physics,'wstar' ,wstar )
@@ -774,6 +887,10 @@ subroutine driver_pbl(itimestep,configs,mesh,sfc_input,diag_physics,tend_physics
config_do_restart, &
bl_mynn_tkeadvect
+ logical,pointer:: config_shinhong_nonlocal_flux, &
+ config_shinhong_scu_mixing, &
+ config_shinhong_ke_dissipation
+
character(len=StrKIND),pointer:: pbl_scheme
integer,pointer:: bl_mynn_cloudpdf, &
@@ -822,6 +939,43 @@ subroutine driver_pbl(itimestep,configs,mesh,sfc_input,diag_physics,tend_physics
pbl_select: select case (trim(pbl_scheme))
+ case("bl_shinhong")
+ call mpas_timer_start('bl_shinhong')
+ call mpas_pool_get_config(configs,'config_shinhong_nonlocal_flux',config_shinhong_nonlocal_flux)
+ call mpas_pool_get_config(configs,'config_shinhong_scu_mixing',config_shinhong_scu_mixing)
+ call mpas_pool_get_config(configs,'config_shinhong_ke_dissipation',config_shinhong_ke_dissipation)
+
+ call shinhong ( &
+ p3d = pres_hyd_p , p3di = pres2_hyd_p , psfc = psfc_p , &
+ t3d = t_p , dz8w = dz_p , pi3d = pi_p , &
+ u3d = u_p , v3d = v_p , qv3d = qv_p , &
+ qc3d = qc_p , qi3d = qi_p , rublten = rublten_p , &
+ rvblten = rvblten_p , rthblten = rthblten_p , rqvblten = rqvblten_p , &
+ rqcblten = rqcblten_p , rqiblten = rqiblten_p , flag_qc = f_qc , &
+ flag_qi = f_qi , cp = cp , g = gravity , &
+ rovcp = rcp , rd = R_d , rovg = rdg , &
+ ep1 = ep_1 , ep2 = ep_2 , karman = karman , &
+ xlv = xlv , rv = R_v , znt = znt_p , &
+ ust = ust_p , hpbl = hpbl_p , psim = psim_p , &
+ psih = psih_p , xland = xland_p , hfx = hfx_p , &
+ qfx = qfx_p , wspd = wspd_p , br = br_p , &
+ dt = dt_pbl , kpbl2d = kpbl_p , exch_h = kzh_p , &
+ exch_m = kzm_p , wstar = wstar_p , delta = delta_p , &
+ shinhong_nonlocal_flux = config_shinhong_nonlocal_flux , &
+ shinhong_dissi_heating = config_shinhong_ke_dissipation , &
+ shinhong_scu_mixing = config_shinhong_scu_mixing , &
+ tke = tkepbl_p , el= elpbl_p , corf=fcell_p , &
+ uoce = uoce_p , voce = voce_p , rthraten = rthraten_p , &
+ u10 = u10_p , v10 = v10_p , ctopo = ctopo_p , &
+ ctopo2 = ctopo2_p , flag_bep = flag_bep , idiff = idiff , &
+ dx = dx_p , &
+ errmsg = errmsg , errflg = errflg , &
+ ids = ids , ide = ide , jds = jds , jde = jde , kds = kds , kde = kde , &
+ ims = ims , ime = ime , jms = jms , jme = jme , kms = kms , kme = kme , &
+ its = its , ite = ite , jts = jts , jte = jte , kts = kts , kte = kte &
+ )
+ call mpas_timer_stop('bl_shinhong')
+
case("bl_ysu")
call mpas_timer_start('bl_ysu')
call ysu ( &
diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F b/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F
index afde4fa523..81d92d2f2e 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_driver_sfclayer.F
@@ -144,6 +144,7 @@ subroutine allocate_sfclayer(configs)
if(.not.allocated(v10_p) ) allocate(v10_p(ims:ime,jms:jme) )
if(.not.allocated(wspd_p) ) allocate(wspd_p(ims:ime,jms:jme) )
if(.not.allocated(xland_p) ) allocate(xland_p(ims:ime,jms:jme) )
+ if(.not.allocated(var2d_p) ) allocate(var2d_p(ims:ime,jms:jme) )
if(.not.allocated(zol_p) ) allocate(zol_p(ims:ime,jms:jme) )
if(.not.allocated(znt_p) ) allocate(znt_p(ims:ime,jms:jme) )
@@ -277,6 +278,7 @@ subroutine deallocate_sfclayer(configs)
if(allocated(v10_p) ) deallocate(v10_p )
if(allocated(wspd_p) ) deallocate(wspd_p )
if(allocated(xland_p) ) deallocate(xland_p )
+ if(allocated(var2d_p) ) deallocate(var2d_p )
if(allocated(zol_p) ) deallocate(zol_p )
if(allocated(znt_p) ) deallocate(znt_p )
@@ -369,7 +371,7 @@ subroutine sfclayer_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
real(kind=RKIND),pointer:: len_disp
real(kind=RKIND),dimension(:),pointer:: meshDensity
- real(kind=RKIND),dimension(:),pointer:: skintemp,sst,xice,xland
+ real(kind=RKIND),dimension(:),pointer:: skintemp,sst,xice,xland,var2d
real(kind=RKIND),dimension(:),pointer:: hpbl,mavail
real(kind=RKIND),dimension(:),pointer:: br,cpm,chs,chs2,cqs2,flhc,flqc,gz1oz0,hfx,qfx, &
qgh,qsfc,lh,mol,psim,psih,regime,rmol,ust,ustm, &
@@ -393,6 +395,7 @@ subroutine sfclayer_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
call mpas_pool_get_array(diag_physics,'mavail' ,mavail )
call mpas_pool_get_array(sfc_input ,'skintemp',skintemp)
call mpas_pool_get_array(sfc_input ,'xland' ,xland )
+ call mpas_pool_get_array(sfc_input ,'var2d' ,var2d )
!inout variables:
call mpas_pool_get_array(diag_physics,'br' ,br )
@@ -427,6 +430,7 @@ subroutine sfclayer_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite)
mavail_p(i,j) = mavail(i)
tsk_p(i,j) = skintemp(i)
xland_p(i,j) = xland(i)
+ var2d_p(i,j) = var2d(i)
!inout variables:
br_p(i,j) = br(i)
@@ -852,12 +856,14 @@ subroutine driver_sfclayer(itimestep,configs,mesh,diag_physics,sfc_input,its,ite
!local pointers:
logical,pointer:: config_do_restart,config_frac_seaice
+ logical,pointer:: config_kim_tofd
character(len=StrKIND),pointer:: sfclayer_scheme
real(kind=RKIND),dimension(:),pointer:: areaCell
!local variables:
integer:: initflag
real(kind=RKIND):: dx
+ real(kind=RKIND), pointer :: tofd_factor
!CCPP-compliant flags:
character(len=StrKIND):: errmsg
@@ -874,6 +880,8 @@ subroutine driver_sfclayer(itimestep,configs,mesh,diag_physics,sfc_input,its,ite
call mpas_pool_get_config(configs,'config_do_restart' ,config_do_restart )
call mpas_pool_get_config(configs,'config_frac_seaice' ,config_frac_seaice)
call mpas_pool_get_config(configs,'config_sfclayer_scheme',sfclayer_scheme )
+ call mpas_pool_get_config(configs,'config_kim_tofd' ,config_kim_tofd )
+ call mpas_pool_get_config(configs,'config_tofd_factor' ,tofd_factor)
call mpas_pool_get_array(mesh,'areaCell',areaCell)
@@ -963,7 +971,9 @@ subroutine driver_sfclayer(itimestep,configs,mesh,diag_physics,sfc_input,its,ite
xland = xland_p , hfx = hfx_p , qfx = qfx_p , &
lh = lh_p , tsk = tsk_p , flhc = flhc_p , &
flqc = flqc_p , qgh = qgh_p , qsfc = qsfc_p , &
- rmol = rmol_p , u10 = u10_p , v10 = v10_p , &
+ rmol = rmol_p , var2d = var2d_p , u10 = u10_p , &
+ if_kim_tofd = config_kim_tofd, tofd_factor = tofd_factor , &
+ v10 = v10_p , &
th2 = th2m_p , t2 = t2m_p , q2 = q2_p , &
gz1oz0 = gz1oz0_p , wspd = wspd_p , br = br_p , &
isfflx = isfflx , dx = dx_p , svp1 = svp1 , &
@@ -993,7 +1003,9 @@ subroutine driver_sfclayer(itimestep,configs,mesh,diag_physics,sfc_input,its,ite
xland = xland_sea , hfx = hfx_sea , qfx = qfx_sea , &
lh = lh_sea , tsk = tsk_sea , flhc = flhc_sea , &
flqc = flqc_sea , qgh = qgh_sea , qsfc = qsfc_sea , &
- rmol = rmol_sea , u10 = u10_sea , v10 = v10_sea , &
+ rmol = rmol_sea , var2d = var2d_p , u10 = u10_sea , &
+ if_kim_tofd = config_kim_tofd, tofd_factor = tofd_factor , &
+ v10 = v10_sea , &
th2 = th2m_sea , t2 = t2m_sea , q2 = q2_sea , &
gz1oz0 = gz1oz0_sea , wspd = wspd_sea , br = br_sea , &
isfflx = isfflx , dx = dx_p , svp1 = svp1 , &
diff --git a/src/core_atmosphere/physics/mpas_atmphys_packages.F b/src/core_atmosphere/physics/mpas_atmphys_packages.F
index 5d32cb297e..9fe858d96e 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_packages.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_packages.F
@@ -39,7 +39,7 @@ function atmphys_setup_packages(configs,packages,iocontext) result(ierr)
character(len=StrKIND),pointer:: config_lsm_scheme
logical,pointer:: mp_kessler_in,mp_thompson_in,mp_thompson_aers_in,mp_wsm6_in
logical,pointer:: cu_grell_freitas_in,cu_kain_fritsch_in,cu_ntiedtke_in
- logical,pointer:: bl_mynn_in,bl_ysu_in
+ logical,pointer:: bl_mynn_in,bl_ysu_in,bl_shinhong_in
logical,pointer:: sf_noahmp_in
integer :: ierr
@@ -150,8 +150,11 @@ function atmphys_setup_packages(configs,packages,iocontext) result(ierr)
nullify(bl_ysu_in)
call mpas_pool_get_package(packages,'bl_ysu_inActive',bl_ysu_in)
+ nullify(bl_shinhong_in)
+ call mpas_pool_get_package(packages,'bl_shinhong_inActive',bl_shinhong_in)
+
if(.not.associated(bl_mynn_in) .or. &
- .not.associated(bl_ysu_in)) then
+ .not.associated(bl_ysu_in) .or. .not.associated(bl_shinhong_in) ) then
call mpas_log_write('====================================================================================',messageType=MPAS_LOG_ERR)
call mpas_log_write('* Error while setting up packages for planetary layer options in atmosphere core.', messageType=MPAS_LOG_ERR)
call mpas_log_write('====================================================================================',messageType=MPAS_LOG_ERR)
@@ -161,15 +164,19 @@ function atmphys_setup_packages(configs,packages,iocontext) result(ierr)
bl_mynn_in = .false.
bl_ysu_in = .false.
+ bl_shinhong_in = .false.
if(config_pbl_scheme=='bl_mynn') then
bl_mynn_in = .true.
elseif(config_pbl_scheme == 'bl_ysu') then
bl_ysu_in = .true.
+ elseif(config_pbl_scheme == 'bl_shinhong') then
+ bl_shinhong_in = .true.
endif
call mpas_log_write(' bl_mynn_in = $l', logicArgs=(/bl_mynn_in/))
call mpas_log_write(' bl_ysu_in = $l', logicArgs=(/bl_ysu_in/))
+ call mpas_log_write(' bl_shinhong_in = $l', logicArgs=(/bl_shinhong_in/))
call mpas_log_write('')
!--- initialization of all packages for parameterizations of land surface processes:
diff --git a/src/core_atmosphere/physics/mpas_atmphys_vars.F b/src/core_atmosphere/physics/mpas_atmphys_vars.F
index 134009f537..5747847d40 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_vars.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_vars.F
@@ -374,11 +374,11 @@ module mpas_atmphys_vars
!=================================================================================================================
logical,parameter:: &
- flag_bep = .false. !flag to use BEP/BEP+BEM for use in the YSU PBL scheme (with urban physics). since we do
+ flag_bep = .false. !flag to use BEP/BEP+BEM for use in the SHINHONG/YSU PBL scheme (with urban physics). since we do
!not run urban physics, flag_bep is always set to false.
integer,parameter:: &
- idiff = 0 !BEP/BEM+BEM diffusion flag for use in the YSU PBL scheme (with urban physics). since we
+ idiff = 0 !BEP/BEM+BEM diffusion flag for use in the SHINHONG/YSU PBL scheme (with urban physics). since we
!do not run urban physics, idiff is set to zero.
integer:: ysu_pblmix
@@ -394,6 +394,7 @@ module mpas_atmphys_vars
hpbl_p, &!PBL height [m]
delta_p, &!
wstar_p, &!
+ fcell_p, &!
uoce_p, &!
voce_p !
@@ -469,7 +470,9 @@ module mpas_atmphys_vars
sina_p !sine of map rotation [-]
real(kind=RKIND),dimension(:,:),allocatable:: &
+ ter_p, &!orographic height [m]
var2d_p, &!orographic variance [m2]
+ elvmax_p, &!orographic maximum [m]
con_p, &!orographic convexity [m2]
oa1_p, &!orographic direction asymmetry function [-]
oa2_p, &!orographic direction asymmetry function [-]
diff --git a/src/core_atmosphere/physics/physics_wrf/Makefile b/src/core_atmosphere/physics/physics_wrf/Makefile
index 4495b74960..3058e079b3 100644
--- a/src/core_atmosphere/physics/physics_wrf/Makefile
+++ b/src/core_atmosphere/physics/physics_wrf/Makefile
@@ -14,6 +14,7 @@ OBJS = \
module_bl_ugwp_gwdo.o \
module_bl_mynn.o \
module_bl_ysu.o \
+ module_bl_shinhong.o \
module_cam_error_function.o \
module_cam_shr_kind_mod.o \
module_cam_support.o \
diff --git a/src/core_atmosphere/physics/physics_wrf/module_bl_gwdo.F b/src/core_atmosphere/physics/physics_wrf/module_bl_gwdo.F
index ae95ed5d62..b7ac54fef6 100644
--- a/src/core_atmosphere/physics/physics_wrf/module_bl_gwdo.F
+++ b/src/core_atmosphere/physics/physics_wrf/module_bl_gwdo.F
@@ -16,9 +16,10 @@ subroutine gwdo(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, &
rublten,rvblten, &
dtaux3d,dtauy3d,dusfcg,dvsfcg, &
var2d,oc12d,oa2d1,oa2d2,oa2d3,oa2d4,ol2d1,ol2d2,ol2d3,ol2d4, &
- sina,cosa,znu,znw,p_top, &
+ elvmax,sina,cosa,znu,znw,p_top, &
cp,g,rd,rv,ep1,pi, &
dt,dx,kpbl2d,itimestep, &
+ dx_factor,if_nonhyd, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
@@ -76,6 +77,8 @@ subroutine gwdo(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, &
integer,intent(in),dimension(ims:ime,jms:jme):: kpbl2d
real(kind=kind_phys),intent(in):: dt,cp,g,rd,rv,ep1,pi
+ real(kind=kind_phys),intent(in):: dx_factor
+ logical ,intent(in):: if_nonhyd
real(kind=kind_phys),intent(in),optional:: p_top
real(kind=kind_phys),intent(in),dimension(kms:kme),optional:: &
@@ -88,6 +91,7 @@ subroutine gwdo(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, &
oc12d, &
oa2d1,oa2d2,oa2d3,oa2d4, &
ol2d1,ol2d2,ol2d3,ol2d4, &
+ elvmax, &
sina,cosa
real(kind=kind_phys),intent(in),dimension(ims:ime,kms:kme,jms:jme):: &
@@ -124,7 +128,7 @@ subroutine gwdo(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, &
integer:: i,j,k
real(kind=kind_phys),dimension(its:ite):: &
- var2d_hv,oc12d_hv,dx_hv,sina_hv,cosa_hv
+ var2d_hv,oc12d_hv,dx_hv,sina_hv,cosa_hv,elvmax_hv
real(kind=kind_phys),dimension(its:ite):: &
oa2d1_hv,oa2d2_hv,oa2d3_hv,oa2d4_hv,ol2d1_hv,ol2d2_hv,ol2d3_hv,ol2d4_hv
real(kind=kind_phys),dimension(its:ite):: &
@@ -193,6 +197,7 @@ subroutine gwdo(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, &
ol2d2_hv(i) = ol2d2(i,j)
ol2d3_hv(i) = ol2d3(i,j)
ol2d4_hv(i) = ol2d4(i,j)
+ elvmax_hv(i) =elvmax(i,j)
enddo
call bl_gwdo_run(sina=sina_hv,cosa=cosa_hv &
@@ -209,6 +214,8 @@ subroutine gwdo(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, &
,oa2d3=oa2d3_hv, oa2d4=oa2d4_hv &
,ol2d1=ol2d1_hv, ol2d2=ol2d2_hv &
,ol2d3=ol2d3_hv, ol2d4=ol2d4_hv &
+ ,omax=elvmax_hv &
+ ,dx_factor=dx_factor,if_nonhyd=if_nonhyd &
,g_=g,cp_=cp,rd_=rd,rv_=rv,fv_=ep1,pi_=pi &
,dxmeter=dx_hv,deltim=dt &
,its=its,ite=ite,kte=kte,kme=kte+1 &
diff --git a/src/core_atmosphere/physics/physics_wrf/module_bl_shinhong.F b/src/core_atmosphere/physics/physics_wrf/module_bl_shinhong.F
new file mode 100644
index 0000000000..dbfe5f01f2
--- /dev/null
+++ b/src/core_atmosphere/physics/physics_wrf/module_bl_shinhong.F
@@ -0,0 +1,499 @@
+#define NEED_B4B_DURING_CCPP_TESTING 1
+!=================================================================================================================
+ module module_bl_shinhong
+ use mpas_kind_types,only: kind_phys => RKIND
+ use bl_shinhong
+
+
+ implicit none
+ private
+ public:: shinhong
+
+
+ contains
+
+
+!=================================================================================================================
+ subroutine shinhong(u3d,v3d,t3d,qv3d,qc3d,qi3d,p3d,p3di,pi3d, &
+ rublten,rvblten,rthblten, &
+ rqvblten,rqcblten,rqiblten,flag_qc,flag_qi, &
+ cp,g,rovcp,rd,rovg,ep1,ep2,karman,xlv,rv, &
+ dz8w,psfc, &
+ znt,ust,hpbl,psim,psih, &
+ xland,hfx,qfx,wspd,br, &
+ dt,kpbl2d, &
+ exch_h,exch_m, &
+ wstar,delta, &
+ shinhong_nonlocal_flux, &
+ tke,el,corf, &
+ u10,v10, &
+ uoce,voce, &
+ rthraten,shinhong_scu_mixing, &
+ shinhong_dissi_heating, &
+ ctopo,ctopo2,dx, &
+ idiff,flag_bep,frc_urb2d, &
+ a_u_bep,a_v_bep,a_t_bep, &
+ a_q_bep, &
+ a_e_bep,b_u_bep,b_v_bep, &
+ b_t_bep,b_q_bep, &
+ b_e_bep,dlg_bep, &
+ dl_u_bep,sf_bep,vl_bep, &
+ ids,ide, jds,jde, kds,kde, &
+ ims,ime, jms,jme, kms,kme, &
+ its,ite, jts,jte, kts,kte, &
+ errmsg,errflg &
+ )
+!-------------------------------------------------------------------------------
+ implicit none
+!-------------------------------------------------------------------------------
+!-- u3d 3d u-velocity interpolated to theta points (m/s)
+!-- v3d 3d v-velocity interpolated to theta points (m/s)
+!-- th3d 3d potential temperature (k)
+!-- t3d temperature (k)
+!-- qv3d 3d water vapor mixing ratio (kg/kg)
+!-- qc3d 3d cloud mixing ratio (kg/kg)
+!-- qi3d 3d ice mixing ratio (kg/kg)
+! (note: if P_QI
-
@@ -463,6 +463,7 @@
+
@@ -568,6 +569,7 @@
+
@@ -932,6 +934,10 @@
+
+
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)
@@ -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)
@@ -290,14 +293,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)
@@ -945,9 +948,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
@@ -987,11 +997,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
@@ -1036,9 +1054,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
@@ -1082,9 +1101,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
@@ -1109,10 +1129,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
@@ -1125,7 +1146,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
@@ -1145,10 +1166,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
@@ -1161,7 +1183,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
@@ -1181,10 +1203,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
@@ -1197,7 +1220,7 @@ 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
@@ -1223,10 +1246,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
@@ -1239,13 +1263,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