diff --git a/src/core_atmosphere/Externals.cfg b/src/core_atmosphere/Externals.cfg
index 84dc47d1d8..393940f5f0 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=KIM_GWDO
required = True
[GSL_UGWP]
diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index 4281c40bba..971a0a5f5b 100644
--- a/src/core_atmosphere/Registry.xml
+++ b/src/core_atmosphere/Registry.xml
@@ -521,6 +521,7 @@
+
@@ -2278,7 +2279,7 @@
+ possible_values="`suite',`bl_kim_gwdo',`bl_ugwp_gwdo',`off'"/>
+
+
+
+
@@ -3782,6 +3799,8 @@
+
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_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_vars.F b/src/core_atmosphere/physics/mpas_atmphys_vars.F
index 134009f537..084fc9f0e0 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_vars.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_vars.F
@@ -469,7 +469,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/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_sf_sfclayrev.F b/src/core_atmosphere/physics/physics_wrf/module_sf_sfclayrev.F
index ac70882989..e3fcd527a9 100644
--- a/src/core_atmosphere/physics/physics_wrf/module_sf_sfclayrev.F
+++ b/src/core_atmosphere/physics/physics_wrf/module_sf_sfclayrev.F
@@ -19,10 +19,11 @@ subroutine sfclayrev(u3d,v3d,t3d,qv3d,p3d,dz8w, &
znt,ust,pblh,mavail,zol,mol,regime,psim,psih, &
fm,fh, &
xland,hfx,qfx,lh,tsk,flhc,flqc,qgh,qsfc,rmol, &
- u10,v10,th2,t2,q2, &
+ var2d,u10,v10,th2,t2,q2, &
gz1oz0,wspd,br,isfflx,dx, &
svp1,svp2,svp3,svpt0,ep1,ep2, &
karman,p1000mb,lakemask, &
+ if_kim_tofd,tofd_factor, &
ids,ide,jds,jde,kds,kde, &
ims,ime,jms,jme,kms,kme, &
its,ite,jts,jte,kts,kte, &
@@ -45,9 +46,12 @@ subroutine sfclayrev(u3d,v3d,t3d,qv3d,p3d,dz8w, &
real(kind=kind_phys),intent(in):: ep1,ep2,karman
real(kind=kind_phys),intent(in):: p1000mb
real(kind=kind_phys),intent(in):: cp,g,rovcp,r,xlv
+ real(kind=kind_phys),intent(in):: tofd_factor
+ logical ,intent(in):: if_kim_tofd
real(kind=kind_phys),intent(in),dimension(ims:ime,jms:jme):: &
dx, &
+ var2d, &
mavail, &
pblh, &
psfc, &
@@ -125,7 +129,7 @@ subroutine sfclayrev(u3d,v3d,t3d,qv3d,p3d,dz8w, &
dz_hv,u_hv,v_hv,qv_hv,p_hv,t_hv
real(kind=kind_phys),dimension(its:ite):: &
- lh_hv,u10_hv,v10_hv,th2_hv,t2_hv,q2_hv
+ lh_hv,var2d_hv,u10_hv,v10_hv,th2_hv,t2_hv,q2_hv
real(kind=kind_phys),dimension(its:ite):: &
ck_hv,cka_hv,cd_hv,cda_hv
@@ -157,6 +161,7 @@ subroutine sfclayrev(u3d,v3d,t3d,qv3d,p3d,dz8w, &
xland_hv(i) = xland(i,j)
lakemask_hv(i) = lakemask(i,j)
water_depth_hv(i) = water_depth(i,j)
+ var2d_hv(i) = var2d(i,j)
do k = kts,kte
dz_hv(i,k) = dz8w(i,k,j)
@@ -209,8 +214,9 @@ subroutine sfclayrev(u3d,v3d,t3d,qv3d,p3d,dz8w, &
rmol=rmol_hv,znt=znt_hv,ust=ust_hv,mavail=mavail_hv, &
zol=zol_hv,mol=mol_hv,regime=regime_hv,psim=psim_hv, &
psih=psih_hv,fm=fm_hv,fh=fh_hv,xland=xland_hv,lakemask=lakemask_hv, &
- hfx=hfx_hv,qfx=qfx_hv,tsk=tsk_hv,u10=u10_hv, &
+ hfx=hfx_hv,qfx=qfx_hv,tsk=tsk_hv,varf=var2d_hv,u10=u10_hv, &
v10=v10_hv,th2=th2_hv,t2=t2_hv,q2=q2_hv,flhc=flhc_hv, &
+ if_kim_tofd=if_kim_tofd,tofd_factor=tofd_factor, &
flqc=flqc_hv,qgh=qgh_hv,qsfc=qsfc_hv,lh=lh_hv, &
gz1oz0=gz1oz0_hv,wspd=wspd_hv,br=br_hv,isfflx=l_isfflx,dx=dx_hv, &
svp1=svp1,svp2=svp2,svp3=svp3,svpt0=svpt0,ep1=ep1,ep2=ep2,karman=karman, &
diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml
index cf4934a81b..e2e3be78dd 100644
--- a/src/core_init_atmosphere/Registry.xml
+++ b/src/core_init_atmosphere/Registry.xml
@@ -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"/>
-
@@ -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