Skip to content

Enable nonhydrostatic option for FV3 #75

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
jtruesdal wants to merge 6 commits into
base: stormspeed
Choose a base branch
from
Open

Enable nonhydrostatic option for FV3 #75

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
cb9c729
enable nonhydrostatic fv3
jtruesdal Dec 19, 2025
f43f27c
Fix variable check for W in initial_conditions
jtruesdal Dec 19, 2025
fb153e4
Fix nlev assignment for W variable in ic_baroclinic.F90
jtruesdal Dec 19, 2025
210db2c
Apply suggestion from @Copilot
jtruesdal Dec 19, 2025
ab080bb
Merge remote-tracking branch 'origin/stormspeed' into fv3nh to update…
jtruesdal Dec 19, 2025
34c8b4b
Merge remote-tracking branch 'origin/fv3nh' with copilot review mods …
jtruesdal Dec 19, 2025
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
6 changes: 3 additions & 3 deletions .gitmodules
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,10 +43,10 @@

[submodule "fv3"]
path = src/dynamics/fv3
url = https://github.com/ESCOMP/CAM_FV3_interface.git
url = https://github.com/jtruesdal/CAM_FV3_interface.git
fxrequired = AlwaysRequired
fxtag = fv3int_061924
fxDONOTUSEurl = https://github.com/ESCOMP/CAM_FV3_interface.git
fxtag = fv3int_nh_121825
fxDONOTUSEurl = https://github.com/jtruesdal/CAM_FV3_interface.git

[submodule "geoschem"]
path = src/chemistry/geoschem/geoschem_src
Expand Down
4 changes: 3 additions & 1 deletion bld/build-namelist
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4108,7 +4108,9 @@ if ($dyn eq 'fv') {

# FV3 dycore
if ( $dyn eq 'fv3') {

add_default($nl, 'fv3_a_imp');
add_default($nl, 'fv3_p_fac');
add_default($nl, 'fv3_use_logp');
add_default($nl, 'fv3_adjust_dry_mass');
add_default($nl, 'fv3_beta');
add_default($nl, 'fv3_clock_grain');
Expand Down
3 changes: 3 additions & 0 deletions bld/namelist_files/namelist_defaults_cam.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3114,6 +3114,9 @@

<fv3_print_memory_usage dyn="fv3">.true.</fv3_print_memory_usage>

<fv3_a_imp > 1.0D0 </fv3_a_imp>
<fv3_p_fac > 0.05D0 </fv3_p_fac>
<fv3_use_logp > .false. </fv3_use_logp>
<fv3_adjust_dry_mass > .false. </fv3_adjust_dry_mass>
<fv3_beta > 0 </fv3_beta>
<fv3_consv_am > .false. </fv3_consv_am>
Expand Down
36 changes: 36 additions & 0 deletions bld/namelist_files/namelist_definition.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9352,6 +9352,42 @@ Default: none

<!-- for fv3 dynamical core -->

<entry id="fv3_a_imp" type="real" category="dyn_fv3"
group="fv_core_nml" valid_values="" >
Real: Controls behavior of the non-hydrostatic solver. Values greater
than 0.5 enable the semi-implicit solver, in which the value of a_imp
controls the time-off-centering: use a_imp = 1.0 for a fully backward
time-stepping. For consistency, the sum of beta and a_imp should
be 1 when the semi-implicit solver is used. The semi-implicit algo-
rithm is substantially more efficient except at very high (km-scale)
resolutions with an acoustic time step of a few seconds or less. 0.75
by default. Proper values are 0, or between 0.5 and 1. This variable is
only used if hydrostatic =.false.
Default: 1.
</entry>

<entry id="fv3_p_fac" type="real" category="dyn_fv3"
group="fv_core_nml" valid_values="" >
Real: Safety factor for minimum nonhydrostatic pressures, which
will be limited so the full pressure is no less than p_fac times the
hydrostatic pressure. This is only of concern in mid-top or high-
top models with very low pressures near the model top, and has
no effect in most simulations. The pressure limiting activates only
when model is in danger of blowup due to unphysical negative
total pressures. 0.05 by default. Only used if hydrostatic =.false.
and the semi-implicit solver is used. Proper range is 0 to 0.25.
Default: 0.05
</entry>

<entry id="fv3_use_logp" type="logical" category="dyn_fv3"
group="fv_core_nml" valid_values="" >
Logical: Enables a variant of the Lin pressure-gradient force
algorithm, which uses the logarithm of pressure instead of the Exner
function (as in Lin 1997). This yields more accurate results for re-
gions that are nearly isothermal. Ignored if hydrostatic = true.
Default: FALSE
</entry>

<entry id="fv3_clock_grain" type="char*256" category="dyn_fv3"
group="fms_nml" valid_values="" >
The level of clock granularity used for performance timing sections
Expand Down
2 changes: 0 additions & 2 deletions cime_config/config_component.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -372,8 +372,6 @@
<value compset="_BGC%BDRD"> co2_cycle_rad_passive=.true. </value>
<!-- hemco (harmonized emissions component) option -->
<value compset="_CAM.*%HEMCO"> use_hemco=.true. </value>
<!-- rearth provided by FMS for FV3 and must remain fixed to that value to maintain consistency between phys/dyn. -->
<value grid="a%C[0-9]"> rearth = 6.37122D6 </value>
</values>
<group>run_component_cam</group>
<file>env_run.xml</file>
Expand Down
2 changes: 1 addition & 1 deletion src/dynamics/fv3
Open in desktop
Submodule fv3 updated 4 files
+9 −0 Externals.cfg
+27 −8 dyn_comp.F90
+2,160 −0 src_override/fv_arrays.F90
+1,931 −0 src_override/nh_utils.F90
54 changes: 45 additions & 9 deletions src/dynamics/tests/inic_analytic.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,8 @@ module inic_analytic
CONTAINS
!==============================================================================

subroutine dyn_set_inic_col(vcoord, latvals, lonvals, glob_ind, zint, U, V, T, &
PS, PHIS_IN, PHIS_OUT, Q, m_cnst, mask, verbose)
subroutine dyn_set_inic_col(vcoord, latvals, lonvals, glob_ind, zint, U, V, W, &
T, PS, PHIS_IN, PHIS_OUT, Q, m_cnst, mask, verbose)
use cam_initfiles, only: pertlim
#ifdef ANALYTIC_IC
use ic_held_suarez, only: hs94_set_ic
Expand All @@ -62,6 +62,7 @@ subroutine dyn_set_inic_col(vcoord, latvals, lonvals, glob_ind, zint, U, V, T, &
real(r8), optional, intent(in) :: zint(:,:) ! height at layer interfaces
real(r8), optional, intent(inout) :: U(:,:) ! zonal velocity
real(r8), optional, intent(inout) :: V(:,:) ! meridional velocity
real(r8), optional, intent(inout) :: W(:,:) ! vertical velocity (nonhydrostatic)
real(r8), optional, intent(inout) :: T(:,:) ! temperature
real(r8), optional, intent(inout) :: PS(:) ! surface pressure
real(r8), optional, intent(in) :: PHIS_IN(:) ! surface geopotential
Expand Down Expand Up @@ -116,6 +117,13 @@ subroutine dyn_set_inic_col(vcoord, latvals, lonvals, glob_ind, zint, U, V, T, &
return
end if
end if
if (present(W)) then
if (size(W) > 0) then
call check_array_size(W(:,1), 'W', latvals, subname)
else
return
end if
end if
if (present(T)) then
if (size(T) > 0) then
call check_array_size(T(:,1), 'T', latvals, subname)
Expand Down Expand Up @@ -160,20 +168,20 @@ subroutine dyn_set_inic_col(vcoord, latvals, lonvals, glob_ind, zint, U, V, T, &
end if
select case(trim(analytic_ic_type))
case('held_suarez_1994')
call hs94_set_ic(latvals, lonvals, U=U, V=V, T=T, PS=PS, PHIS=PHIS_OUT, &
call hs94_set_ic(latvals, lonvals, U=U, V=V, W=W, T=T, PS=PS, PHIS=PHIS_OUT, &
Q=Q, m_cnst=m_cnst, mask=mask_use, verbose=verbose_use)

case('moist_baroclinic_wave_dcmip2016', 'dry_baroclinic_wave_dcmip2016')
call bc_wav_set_ic(vcoord, latvals, lonvals, zint=zint, U=U, V=V, T=T, PS=PS, &
PHIS=PHIS_OUT, Q=Q, m_cnst=m_cnst, mask=mask_use, verbose=verbose_use)
call bc_wav_set_ic(vcoord, latvals, lonvals, zint=zint, U=U, V=V, W=W, T=T, &
PS=PS, PHIS=PHIS_OUT, Q=Q, m_cnst=m_cnst, mask=mask_use, verbose=verbose_use)

case('dry_baroclinic_wave_jw2006')
call bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U=U, V=V, T=T, PS=PS, &
call bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U=U, V=V, W=W, T=T, PS=PS, &
PHIS=PHIS_OUT, Q=Q, m_cnst=m_cnst, mask=mask_use, verbose=verbose_use)

case('us_standard_atmosphere')
call us_std_atm_set_ic(latvals, lonvals, zint=zint, U=U, V=V, T=T, PS=PS, PHIS_IN=PHIS_IN, &
PHIS_OUT=PHIS_OUT, Q=Q, m_cnst=m_cnst, mask=mask_use, verbose=verbose_use)
call us_std_atm_set_ic(latvals, lonvals, zint=zint, U=U, V=V, W=W, T=T, PS=PS, &
PHIS_IN=PHIS_IN, PHIS_OUT=PHIS_OUT, Q=Q, m_cnst=m_cnst, mask=mask_use, verbose=verbose_use)

case default
call endrun(subname//': Unknown analytic_ic_type, "'//trim(analytic_ic_type)//'"')
Expand Down Expand Up @@ -216,7 +224,7 @@ subroutine dyn_set_inic_col(vcoord, latvals, lonvals, glob_ind, zint, U, V, T, &

end subroutine dyn_set_inic_col

subroutine dyn_set_inic_cblock(vcoord,latvals, lonvals, glob_ind, U, V, T, &
subroutine dyn_set_inic_cblock(vcoord,latvals, lonvals, glob_ind, U, V, W, T, &
PS, PHIS_IN, PHIS_OUT, Q, m_cnst, mask)
!-----------------------------------------------------------------------
!
Expand All @@ -231,6 +239,7 @@ subroutine dyn_set_inic_cblock(vcoord,latvals, lonvals, glob_ind, U, V, T, &
integer, intent(in) :: glob_ind(:) ! global column index
real(r8), optional, intent(inout) :: U(:,:,:) ! zonal velocity
real(r8), optional, intent(inout) :: V(:,:,:) ! meridional velocity
real(r8), optional, intent(inout) :: W(:,:,:) ! vertical velocity (nonhydrostatic)
real(r8), optional, intent(inout) :: T(:,:,:) ! temperature
real(r8), optional, intent(inout) :: PS(:,:) ! surface pressure
real(r8), optional, intent(in) :: PHIS_IN(:,:)! surface geopotential
Expand Down Expand Up @@ -259,6 +268,9 @@ subroutine dyn_set_inic_cblock(vcoord,latvals, lonvals, glob_ind, U, V, T, &
if(present(V)) then
call get_input_shape(V, 'V', mname, size1, size2, size3, subname)
end if
if(present(W)) then
call get_input_shape(W, 'W', mname, size1, size2, size3, subname)
end if
if(present(T)) then
call get_input_shape(T, 'T', mname, size1, size2, size3, subname)
end if
Expand Down Expand Up @@ -301,6 +313,10 @@ subroutine dyn_set_inic_cblock(vcoord,latvals, lonvals, glob_ind, U, V, T, &
call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
glob_ind(bbeg:bend), V=V(:,:,i), mask=mask(bbeg:bend), verbose=verbose)
end if
if (present(W)) then
call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
glob_ind(bbeg:bend), W=W(:,:,i), mask=mask(bbeg:bend), verbose=verbose)
end if
if (present(PS).and.present(PHIS_IN).and.present(T)) then
call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
glob_ind(bbeg:bend), PS=PS(:,i), PHIS_IN=PHIS_IN(:,i), T=T(:,:,i), &
Expand Down Expand Up @@ -333,6 +349,10 @@ subroutine dyn_set_inic_cblock(vcoord,latvals, lonvals, glob_ind, U, V, T, &
call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
glob_ind(bbeg:bend), V=V(:,:,i), verbose=verbose)
end if
if (present(W)) then
call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
glob_ind(bbeg:bend), W=W(:,:,i), verbose=verbose)
end if
if (present(PS).and.present(PHIS_IN).and.present(T)) then
call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
glob_ind(bbeg:bend), PHIS_IN=PHIS_IN(:,i),PS=PS(:,i),T=T(:,:,i), &
Expand Down Expand Up @@ -382,6 +402,10 @@ subroutine dyn_set_inic_cblock(vcoord,latvals, lonvals, glob_ind, U, V, T, &
call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
glob_ind(bbeg:bend), V=V(:,i,:), mask=mask(bbeg:bend), verbose=verbose)
end if
if (present(W)) then
call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
glob_ind(bbeg:bend), W=W(:,i,:), mask=mask(bbeg:bend), verbose=verbose)
end if
if (present(PS).and.present(PHIS_IN).and.present(T)) then
call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
glob_ind(bbeg:bend), PHIS_IN=PHIS_IN(:,i),PS=PS(:,i),T=T(:,i,:), &
Expand Down Expand Up @@ -414,6 +438,10 @@ subroutine dyn_set_inic_cblock(vcoord,latvals, lonvals, glob_ind, U, V, T, &
call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
glob_ind(bbeg:bend), V=V(:,i,:), verbose=verbose)
end if
if (present(W)) then
call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
glob_ind(bbeg:bend), W=W(:,i,:), verbose=verbose)
end if
if (present(PS).and.present(PHIS_IN).and.present(T)) then
call dyn_set_inic_col(vcoord,latvals(bbeg:bend), lonvals(bbeg:bend), &
glob_ind(bbeg:bend), PHIS_IN=PHIS_IN(:,i),PS=PS(:,i),T=T(:,i,:), &
Expand Down Expand Up @@ -463,6 +491,10 @@ subroutine dyn_set_inic_cblock(vcoord,latvals, lonvals, glob_ind, U, V, T, &
call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
V=V(:,i,:), verbose=verbose)
end if
if (present(W)) then
call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
W=W(:,i,:), verbose=verbose)
end if
if (present(PS).and.present(PHIS_IN).and.present(T)) then
call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
PS=PS(:,i),T=T(:,i,:),PHIS_IN=PHIS_IN(:,i), verbose=verbose)
Expand Down Expand Up @@ -511,6 +543,10 @@ subroutine dyn_set_inic_cblock(vcoord,latvals, lonvals, glob_ind, U, V, T, &
call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
V=V(:,:,i), verbose=verbose)
end if
if (present(W)) then
call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
W=W(:,:,i), verbose=verbose)
end if
if (present(PS).and.present(PHIS_IN).and.present(T)) then
call dyn_set_inic_col(vcoord,lat_use, lonvals, glob_ind(bbeg:bend), &
T=T(:,:,i),PS=PS(:,i), PHIS_IN=PHIS_IN(:,i), verbose=verbose)
Expand Down
25 changes: 19 additions & 6 deletions src/dynamics/tests/initial_conditions/ic_baro_dry_jw06.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ module ic_baro_dry_jw06

contains

subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, &
subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, W, T, PS, PHIS, &
Q, m_cnst, mask, verbose)
use dyn_tests_utils, only: vc_moist_pressure, vc_dry_pressure, vc_height
use constituents, only: cnst_name
Expand All @@ -67,6 +67,7 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, &
! z_k for vccord 1)
real(r8), optional, intent(inout) :: U(:,:) ! zonal velocity
real(r8), optional, intent(inout) :: V(:,:) ! meridional velocity
real(r8), optional, intent(inout) :: W(:,:) ! vertical velocity (nonhydrostatic)
real(r8), optional, intent(inout) :: T(:,:) ! temperature
real(r8), optional, intent(inout) :: PS(:) ! surface pressure
real(r8), optional, intent(out) :: PHIS(:) ! surface geopotential
Expand All @@ -77,7 +78,7 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, &
! Local variables
logical, allocatable :: mask_use(:)
logical :: verbose_use
logical :: lu,lv,lt,lq,l3d_vars
logical :: lu,lv,lw,lt,lq,l3d_vars
integer :: i, k, m
integer :: ncol
integer :: nlev
Expand All @@ -91,8 +92,8 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, &
real(r8) :: phi_vertical
real(r8) :: u_wind(size(latvals))

a_omega = rearth*omega
exponent = rair*gamma/gravit
a_omega = rearth*omega
exponent = rair*gamma/gravit

allocate(mask_use(size(latvals)))
if (present(mask)) then
Expand Down Expand Up @@ -163,13 +164,15 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, &
!
lu = present(U)
lv = present(V)
lw = present(W)
lT = present(T)
lq = present(Q)
l3d_vars = lu .or. lv .or. lt .or.lq
l3d_vars = lu .or. lv .or. lw .or. lt .or.lq
nlev = -1
if (l3d_vars) then
if (lu) nlev = size(U, 2)
if (lv) nlev = size(V, 2)
if (lw) nlev = size(W, 2)
if (lt) nlev = size(T, 2)
if (lq) nlev = size(Q, 2)

Expand Down Expand Up @@ -201,6 +204,16 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, &
write(iulog,*) ' V initialized by "',subname,'"'
end if
end if
if (lw) then
do k = 1, nlev
where(mask_use)
W(:,k) = 0.0_r8
end where
end do
if(masterproc.and. verbose_use) then
write(iulog,*) ' W (nonhydrostatic) initialized by "',subname,'"'
end if
end if
if (lt) then
do k = 1, nlev
eta = hyam(k) + hybm(k)
Expand Down Expand Up @@ -232,7 +245,7 @@ subroutine bc_dry_jw06_set_ic(vcoord, latvals, lonvals, U, V, T, PS, PHIS, &
end where
end do
if(masterproc.and. verbose_use) then
write(iulog,*) ' ', trim(cnst_name(m_cnst(1))), ' initialized by "',subname,'"'
write(iulog,*) ' ', trim(cnst_name(m_cnst(1))), ' initialized by "',subname,'"'
end if
end if
end if
Expand Down
22 changes: 18 additions & 4 deletions src/dynamics/tests/initial_conditions/ic_baroclinic.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ module ic_baroclinic

contains

subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, T, PS, PHIS, &
subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, W, T, PS, PHIS, &
Q, m_cnst, mask, verbose)
use dyn_tests_utils, only: vc_moist_pressure, vc_dry_pressure, vc_height
use constituents, only: cnst_name
Expand All @@ -93,6 +93,7 @@ subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, T, PS, PHIS, &
real(r8), optional, intent(in) :: zint(:,:) ! interface height (ncol,ilev), ordered top to bottom
real(r8), optional, intent(inout) :: U(:,:) ! zonal velocity
real(r8), optional, intent(inout) :: V(:,:) ! meridional velocity
real(r8), optional, intent(inout) :: W(:,:) ! vertical velocity (nonhydrostatic)
real(r8), optional, intent(inout) :: T(:,:) ! temperature
real(r8), optional, intent(inout) :: PS(:) ! surface pressure
real(r8), optional, intent(out) :: PHIS(:) ! surface geopotential
Expand All @@ -113,7 +114,7 @@ subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, T, PS, PHIS, &
real(r8) :: uk,vk,Tvk,qk,pk !mid-level state
real(r8) :: psurface
real(r8) :: wvp,qdry
logical :: lU, lV, lT, lQ, l3d_vars
logical :: lu, lv, lw, lt, lq, l3d_vars
logical :: cnst1_is_moisture
real(r8), allocatable :: pdry_half(:), pwet_half(:),zdry_half(:),zk(:)
real(r8), allocatable :: zmid(:,:) ! layer midpoint heights for test tracer initialization
Expand Down Expand Up @@ -214,13 +215,15 @@ subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, T, PS, PHIS, &
!
lu = present(U)
lv = present(V)
lT = present(T)
lw = present(W)
lt = present(T)
lq = present(Q)
l3d_vars = lu .or. lv .or. lt .or. lq
l3d_vars = lu .or. lv .or. lw .or. lt .or. lq
nlev = -1
if (l3d_vars) then
if (lu) nlev = size(U, 2)
if (lv) nlev = size(V, 2)
if (lw) nlev = size(W, 2)
if (lt) nlev = size(T, 2)

if (lq) then
Expand Down Expand Up @@ -334,6 +337,17 @@ subroutine bc_wav_set_ic(vcoord,latvals, lonvals, zint, U, V, T, PS, PHIS, &
' ', trim(cnst_name(m_cnst(1))), ' initialized by "',subname,'"'
end if

if (lw) then
do k = 1, nlev
where(mask_use)
W(:,k) = 0.0_r8
end where
end do
if(masterproc.and. verbose_use) then
write(iulog,*) ' W (nonhydrostatic) initialized by "',subname,'"'
end if
end if

if (lq) then

ncnst = size(m_cnst, 1)
Expand Down
Loading