Skip to content

Update the kpp test to include running with StokesMOST #102

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

Merged
mnlevy1981 merged 7 commits into from
Jan 24, 2025
Merged

Update the kpp test to include running with StokesMOST #102

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
1d19506
Update build for intel
mnlevy1981 Jan 17, 2025
c601b04
Update build with intel
mnlevy1981 Jan 17, 2025
94e2c38
Add StokesMOST test to reg_tests/kpp
mnlevy1981 Jan 17, 2025
2f2952a
Add call to cvmix_kpp_compute_StokesXi() in test
mnlevy1981 Jan 17, 2025
b419733
Fix invalid argument to compute_StokesXI
mnlevy1981 Jan 17, 2025
369a285
Need to compute denom inside k loop
mnlevy1981 Jan 17, 2025
529b922
Bug fixes
mnlevy1981 Jan 23, 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 bld/CompileFlags.mak
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ ifeq ($(FC),ftn)
FC_TMP = pgf90
endif
ifeq ($(PE_ENV),INTEL)
FC_TMP = ifort
FC_TMP = ifx
endif
ifeq ($(PE_ENV),PATHSCALE)
FC_TMP = pathf95
Expand All @@ -27,8 +27,8 @@ ifeq ($(FC_TMP),pgf90)
FCFLAGS = -O2 -Mfree -module $(OBJ_DIR)
endif

ifeq ($(FC_TMP),ifort)
FCFLAGS = -O2 -free -module $(OBJ_DIR) -cpp -warn all -diag-error warn -nogen-interface -fp-model source
ifeq ($(FC_TMP),ifx)
FCFLAGS = -O0 -g -fpe0 -free -module $(OBJ_DIR) -cpp -warn all -diag-error warn -nogen-interface -fp-model source
endif

ifeq ($(FC_TMP),xlf90)
Expand Down
10 changes: 10 additions & 0 deletions reg_tests/kpp/input.nl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,3 +52,13 @@ mix_type = 'kpp'
Cp0 = 3992.0d0
OBL_depth6 = 6000.0d0
/

! Test 7 params
&kpp_col7_nml
ltest7 = .true.
nlev7 = 10
layer_thick7 = 5.0d0
hmix7 = 17.0d0
! Parameter settings to match LMD94 (linear interp, average Nsqr)
interp_type_t7 = "linear"
/
9 changes: 7 additions & 2 deletions src/Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,12 @@ ifeq ($(USE_NETCDF),TRUE)
LINKING_FLAGS += -L$(NETCDF_LIB) -lnetcdff -lnetcdf
else
# get flibs from netcdf.pc
FLIBS = $(subst flibs=,,$(shell grep flibs $(NETCDF_PC)))
# older versions of netCDF have "flibs=" line
FLIBS = $(shell grep flibs $(NETCDF_PC) | cut -d '=' -f 2-)
ifeq ($(wildcard $(FLIBS)),)
# newer versions of netCDF have Libs: line (and use {libdir} in output
FLIBS = $(shell grep Libs $(NETCDF_PC) | grep -v private | cut -d '}' -f 2-)
endif
# Workaround for yellowstone, I need to figure out where this
# @NC_FLIBS@ comes from
ifeq ($(FLIBS),@NC_FLIBS@)
Expand Down Expand Up @@ -145,7 +150,7 @@ endif

all: exe

$(EXE): $(addprefix $(OBJ_DIR)/,$(OBJS) $(DRIVER_OBJS)) \
$(EXE): $(addprefix $(OBJ_DIR)/,$(DRIVER_OBJS) $(OBJS)) \
$(LIB_DIR)/libcvmix.a
$(FC) -o $(EXE) $(addprefix $(OBJ_DIR)/,$(OBJS) $(DRIVER_OBJS)) $(LINKING_FLAGS)

Expand Down
6 changes: 6 additions & 0 deletions src/cvmix_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,12 @@ Program cvmix_driver
real(kind=cvmix_r8) :: ocn_depth
character(len=cvmix_strlen) :: mix_type

external cvmix_BL_driver
external cvmix_shear_driver
external cvmix_tidal_driver
external cvmix_ddiff_driver
external cvmix_kpp_driver

namelist/cvmix_nml/mix_type, nlev, max_nlev, ocn_depth

mix_type = 'unset'
Expand Down
185 changes: 179 additions & 6 deletions src/drivers/cvmix_kpp_drv.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ Subroutine cvmix_kpp_driver()
cvmix_kpp_compute_unresolved_shear, &
cvmix_kpp_compute_turbulent_scales, &
cvmix_kpp_compute_shape_function_coeffs, &
cvmix_kpp_compute_StokesXi, &
cvmix_coeffs_kpp
use cvmix_put_get, only : cvmix_put
use cvmix_io, only : cvmix_io_open, &
Expand All @@ -41,7 +42,7 @@ Subroutine cvmix_kpp_driver()
real(cvmix_r8), parameter :: third = cvmix_one / real(3,cvmix_r8)

! CVMix datatypes
type(cvmix_data_type) :: CVmix_vars1, CVmix_vars4, CVmix_vars5
type(cvmix_data_type) :: CVmix_vars1, CVmix_vars4, CVmix_vars5, CVmix_vars7

real(cvmix_r8), dimension(:), allocatable, target :: Mdiff, Tdiff, Sdiff
real(cvmix_r8), dimension(:), allocatable, target :: zt, zw_iface, &
Expand All @@ -52,19 +53,25 @@ Subroutine cvmix_kpp_driver()
delta_vel_sqr, &
buoy_freq_iface
real(cvmix_r8), dimension(:,:), allocatable, target :: hor_vel
real(cvmix_r8), dimension(:), allocatable :: uS, vS, uSbar, vSbar
real(cvmix_r8), dimension(:), allocatable :: zeros, ones
real(cvmix_r8), dimension(2) :: ref_vel
real(cvmix_r8), dimension(4) :: shape_coeffs
integer :: i, fid, kt, kw, nlev1, nlev3, nlev4, max_nlev4, OBL_levid4, nlev5
real(cvmix_r8) :: hmix1, hmix5, ri_crit, layer_thick1, layer_thick4, &
layer_thick5, OBL_depth4, OBL_depth5, N, Nsqr
integer :: i, fid, kt, kw, nlev1, nlev3, nlev4, max_nlev4, OBL_levid4, &
nlev5, nlev7
real(cvmix_r8) :: hmix1, hmix5, hmix7, ri_crit, layer_thick1, &
layer_thick4, layer_thick5, layer_thick7, OBL_depth4, &
OBL_depth5, OBL_depth7, N, Nsqr
real(cvmix_r8) :: StokesXI
real(cvmix_r8) :: kOBL_depth, Bslope, Vslope
real(cvmix_r8) :: sigma6, OBL_depth6, surf_buoy_force6, surf_fric_vel6, &
vonkarman6, tao, rho0, grav, alpha, Qnonpen, Cp0, &
w_m6, w_s6, wm6_true, ws6_true
character(len=cvmix_strlen) :: interp_type_t1, interp_type_t4, interp_type_t5
character(len=cvmix_strlen) :: interp_type_t1, interp_type_t4, &
interp_type_t5, interp_type_t7
character(len=cvmix_strlen) :: filename
! True => run specified test
logical :: ltest1, ltest2, ltest3, ltest4, ltest5, ltest6
logical :: ltest1, ltest2, ltest3, ltest4, ltest5, ltest6, ltest7
logical :: lnoDGat1 ! True => G'(1) = 0 (in test 4)

namelist/kpp_col1_nml/ltest1, nlev1, layer_thick1, interp_type_t1, hmix1, &
Expand All @@ -75,6 +82,7 @@ Subroutine cvmix_kpp_driver()
namelist/kpp_col5_nml/ltest5, nlev5, layer_thick5, hmix5, interp_type_t5
namelist/kpp_col6_nml/ltest6, vonkarman6, tao, rho0, grav, alpha, Qnonpen, &
Cp0, OBL_depth6
namelist/kpp_col7_nml/ltest7, nlev7, layer_thick7, hmix7, interp_type_t7

! Read namelists

Expand Down Expand Up @@ -117,12 +125,20 @@ Subroutine cvmix_kpp_driver()
Cp0 = real(3992, cvmix_r8)
OBL_depth6 = real(6000, cvmix_r8)

! Defaults for test 7
ltest7 = .false.
nlev7 = 10
layer_thick7 = real(5, cvmix_r8)
hmix7 = real(17, cvmix_r8)
interp_type_t7 = "linear"

read(*, nml=kpp_col1_nml)
read(*, nml=kpp_col2_nml)
read(*, nml=kpp_col3_nml)
read(*, nml=kpp_col4_nml)
read(*, nml=kpp_col5_nml)
read(*, nml=kpp_col6_nml)
read(*, nml=kpp_col7_nml)

! Test 1: user sets up levels via namelist (constant thickness) and specifies
! critical Richardson number as well as depth parameter hmix1. The
Expand Down Expand Up @@ -590,6 +606,163 @@ Subroutine cvmix_kpp_driver()

end if ! ltest6

! Test 7: Test 5, but lStokesMOST = .true.
if (ltest7) then
print*, ""
print*, "Test 7: Computing Bulk Richardson number (with StokesMOST)"
print*, "----------"

! using linear interpolation, averaging Nsqr, and setting Cv = 1.5 to
! match LMD result
call cvmix_init_kpp(Cv = 1.5_cvmix_r8, interp_type=interp_type_t7, &
lStokesMOST=.true., lMonOb=.true.)

! Set up vertical levels (centers and interfaces) and compute bulk
! Richardson number
allocate(zt(nlev7), zw_iface(nlev7+1))
do kw=1,nlev7+1
zw_iface(kw) = -layer_thick7*real(kw-1, cvmix_r8)
end do
do kt=1,nlev7
zt(kt) = 0.5_cvmix_r8 * (zw_iface(kt)+zw_iface(kt+1))
end do

! Compute Br-B(d), |Vr-V(d)|^2, and Vt^2
allocate(buoyancy(nlev7), delta_vel_sqr(nlev7), hor_vel(nlev7,2), &
shear_sqr(nlev7), w_s(nlev7), Ri_bulk(nlev7), Ri_bulk2(nlev7), &
buoy_freq_iface(nlev7+1))

ref_vel(1) = 0.1_cvmix_r8
ref_vel(2) = cvmix_zero
N = 0.01_cvmix_r8
Nsqr = N*N
Bslope = -Nsqr
Vslope = -0.1_cvmix_r8 / (real(nlev7,cvmix_r8)*layer_thick7-hmix7)
do kt=1,nlev7
if ((zt(kt).ge.-hmix7).or.(kt.eq.1)) then
buoyancy(kt) = Nsqr
hor_vel(kt,1) = 0.1_cvmix_r8
buoy_freq_iface(kt) = cvmix_zero
else
if (zw_iface(kt).ge.-hmix7) then
! derivatives of buoyancy and horizontal velocity component are
! discontinuous in this layer (no change -> non-zero linear change)
! so we compute area-average of analytic function over layer
buoyancy(kt) = Bslope*(-zw_iface(kt+1)-real(hmix7,cvmix_r8))**2 / &
(real(2,cvmix_r8)*layer_thick7) + Nsqr
hor_vel(kt,1) = Vslope*(-zw_iface(kt+1)-real(hmix7,cvmix_r8))**2 / &
(real(2,cvmix_r8)*layer_thick7) + 0.1_cvmix_r8
else
buoyancy(kt) = Nsqr+Bslope*(-zt(kt)-real(hmix7,cvmix_r8))
hor_vel(kt,1) = 0.1_cvmix_r8 + Vslope * &
(-zt(kt)-real(hmix7,cvmix_r8))
end if
buoy_freq_iface(kt) = sqrt(-(buoyancy(kt)-buoyancy(kt-1)) / &
layer_thick7)
end if
! Compute w_s with zeta=0 per LMD page 393
! => w_s = von Karman * surf_fric_vel = 0.4*0.01 = 4e-3
call cvmix_kpp_compute_turbulent_scales(cvmix_zero, -zt(kt), &
buoyancy(1), 0.01_cvmix_r8, &
w_s=w_s(kt))
hor_vel(kt,2) = cvmix_zero
delta_vel_sqr(kt) = (ref_vel(1)-hor_vel(kt,1))**2 + &
(ref_vel(2)-hor_vel(kt,2))**2
end do
buoy_freq_iface(nlev7+1) = N
! MNL: test both uses of compute_bulk_Richardson
Ri_bulk = cvmix_kpp_compute_bulk_Richardson(zt, (buoyancy(1)-buoyancy), &
delta_vel_sqr, &
Nsqr_iface = buoy_freq_iface**2, &
ws_cntr = w_s)

shear_sqr = cvmix_kpp_compute_unresolved_shear(zt, w_s, &
Nsqr_iface = buoy_freq_iface**2)
! Note that Vt_shear_sqr is the fourth argument in compute_bulk_Richardson
! so it does not need to declared explicitly (even though it is optional)
Ri_bulk2 = cvmix_kpp_compute_bulk_Richardson(zt, (buoyancy(1)-buoyancy), &
delta_vel_sqr, shear_sqr)
allocate(zeros(nlev7), source=cvmix_zero)
allocate(ones(nlev7), &
uS(nlev7), &
vS(nlev7), &
uSbar(nlev7), &
vSbar(nlev7), &
source=cvmix_one)
StokesXI = cvmix_one
call cvmix_kpp_compute_StokesXi(zi=zt, &
zk=zw_iface, &
kSL=(nlev7+1)/2, &
SLDepth=cvmix_zero, &
surf_buoy_force=cvmix_one, &
surf_fric_vel=cvmix_one, &
omega_w2x=cvmix_one, &
uE=ones, &
vE=ones, &
uS=uS, &
vS=vS, &
uSbar=uSbar, &
vSbar=vSbar, &
uS_SLD=cvmix_one, &
vS_SLD=cvmix_one, &
uSbar_SLD=cvmix_one, &
vSbar_SLD=cvmix_one, &
StokesXI=StokesXI)
call cvmix_kpp_compute_OBL_depth(Ri_bulk, zw_iface, OBL_depth7, &
kOBL_depth, zt, surf_fric=cvmix_one, &
surf_buoy=ones, Xi=ones)
deallocate(zeros, ones, uS, vS, uSbar, vSbar)
do kt=1,nlev7
if (abs(Ri_bulk(kt)-Ri_bulk2(kt)).gt.1e-12_cvmix_r8) then
print*, "WARNING: two Ri_bulk computations did not match!"
print*, zt(kt), Ri_bulk(kt), Ri_bulk2(kt)
else
print*, zt(kt), Ri_bulk(kt)
end if
end do
print*, "OBL has depth of ", OBL_depth7
#ifdef _NETCDF
print*, "Done! Data is stored in test7.nc, run plot_bulk_Rich.ncl ", &
"to see output."
#else
print*, "Done! Data is stored in test7.out, run plot_bulk_Rich.ncl ", &
"to see output."
#endif

CVmix_vars7%nlev = nlev7
CVmix_vars7%BoundaryLayerDepth = OBL_depth7
CVmix_vars7%zt_cntr => zt
CVmix_vars7%BulkRichardson_cntr => Ri_bulk
CVmix_vars7%Vx_cntr => hor_vel(:,1)
#ifdef _NETCDF
call cvmix_io_open(fid, "test7.nc", "nc")
#else
call cvmix_io_open(fid, "test7.out", "ascii")
#endif
call cvmix_output_write(fid, CVmix_vars7, (/"zt ", &
"Ri_bulk ", &
"Vx ", &
"buoyancy"/), buoyancy)
#ifdef _NETCDF
call cvmix_output_write_att(fid, "OBL_depth", &
CVmix_vars7%BoundaryLayerDepth)
call cvmix_output_write_att(fid, "longname", "ocean height (cell center)",&
"zt")
call cvmix_output_write_att(fid, "units", "m", "zt")
call cvmix_output_write_att(fid, "longname", "horizontal velocity", "U")
call cvmix_output_write_att(fid, "units", "m/s", "U")
call cvmix_output_write_att(fid, "units", "m/s^2", "buoyancy")
call cvmix_output_write_att(fid, "longname", "Bulk Richardson number", &
"BulkRichardson")
call cvmix_output_write_att(fid, "units", "unitless", "BulkRichardson")
#endif
call cvmix_io_close(fid)

deallocate(zt, zw_iface)
deallocate(buoyancy, delta_vel_sqr, hor_vel, shear_sqr, w_s, Ri_bulk, &
Ri_bulk2, buoy_freq_iface)
end if ! ltest7

!EOC

End Subroutine cvmix_kpp_driver
32 changes: 17 additions & 15 deletions src/shared/cvmix_kpp.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1576,7 +1576,7 @@ subroutine cvmix_kpp_compute_OBL_depth_low(Ri_bulk, zw_iface, OBL_depth, &
! Local variables
real(kind=cvmix_r8), dimension(:), pointer :: depth
real(kind=cvmix_r8), dimension(4) :: coeffs
real(kind=cvmix_r8) :: Ekman, MoninObukhov, OBL_Limit
real(kind=cvmix_r8) :: Ekman, MoninObukhov, OBL_Limit, numer, denom
integer :: nlev, k, kRi
logical :: lstable

Expand Down Expand Up @@ -1671,16 +1671,18 @@ subroutine cvmix_kpp_compute_OBL_depth_low(Ri_bulk, zw_iface, OBL_depth, &
if (CVmix_kpp_params_in%lMonOb ) then
if ( present(Xi) .and. present(surf_buoy) ) then
MoninObukhov = OBL_limit
numer = surf_fric**3
do k = 0, kRi-1
if (surf_buoy(k+1) .gt. cvmix_zero) MoninObukhov = &
surf_fric**3 / (surf_buoy(k+1) * (cvmix_one-Xi(k+1)))
denom = surf_buoy(k+1) * (cvmix_one-Xi(k+1))
if ( denom*OBL_limit .gt. numer ) MoninObukhov = &
numer / denom
if ( MoninObukhov .lt. abs(zt_cntr(k+1)) ) &
exit
end do
if (k.eq.0) then
OBL_limit = abs(zt_cntr(1))
elseif (k.lt.kRi) then
OBL_limit = min( OBL_limit, abs(zw_iface(k)) )
OBL_limit = min( OBL_limit, abs(zw_iface(k+1)) )
end if

else
Expand All @@ -1691,7 +1693,6 @@ subroutine cvmix_kpp_compute_OBL_depth_low(Ri_bulk, zw_iface, OBL_depth, &

! (4) OBL_depth must be at or above OBL_limit -zt_cntr(1) < OBL_depth < OBL_limit
OBL_depth = min(OBL_depth, OBL_limit )
kOBL_depth = cvmix_kpp_compute_kOBL_depth(zw_iface, zt_cntr, OBL_depth)

else ! not Stokes_MOST

Expand Down Expand Up @@ -1763,13 +1764,13 @@ subroutine cvmix_kpp_compute_OBL_depth_low(Ri_bulk, zw_iface, OBL_depth, &
! because we know OBL_depth will equal OBL_limit?
OBL_depth = min(OBL_depth, OBL_limit)
end if
end if ! lStokesMOST

OBL_depth = max(OBL_depth, CVmix_kpp_params_in%minOBLdepth)
if (CVmix_kpp_params_in%maxOBLdepth.gt.cvmix_zero) &
OBL_depth = min(OBL_depth, CVmix_kpp_params_in%maxOBLdepth)
kOBL_depth = cvmix_kpp_compute_kOBL_depth(zw_iface, zt_cntr, OBL_depth)
OBL_depth = max(OBL_depth, CVmix_kpp_params_in%minOBLdepth)
if (CVmix_kpp_params_in%maxOBLdepth.gt.cvmix_zero) &
OBL_depth = min(OBL_depth, CVmix_kpp_params_in%maxOBLdepth)
kOBL_depth = cvmix_kpp_compute_kOBL_depth(zw_iface, zt_cntr, OBL_depth)

end if ! lStokesMOST

!EOC

Expand Down Expand Up @@ -3344,12 +3345,13 @@ subroutine cvmix_kpp_composite_Gshape(sigma , Gat1, Gsig, dGdsig)
!EOP
!BOC

real(cvmix_r8) :: a2Gsig, a3Gsig, bGsig, sig_m, G_m, G_1, sig
real(cvmix_r8), parameter :: a2Gsig = -2.1637_cvmix_r8, &
a3Gsig = 0.5831_cvmix_r8, &
sig_m = 0.35_cvmix_r8, &
G_m = 0.11_cvmix_r8

real(cvmix_r8) :: bGsig, G_1, sig

a2Gsig = -2.1637_cvmix_r8
a3Gsig = 0.5831_cvmix_r8
sig_m = 0.35_cvmix_r8
G_m = 0.11_cvmix_r8 ! sig_m + sig_m * sig_m * (a2Gsig + a3Gsig * sig_m)
G_1 = MAX( cvmix_zero , MIN( Gat1 , G_m ) )

if (sigma .lT. sig_m) then
Expand Down