Skip to content

Commit 65ef5c7

Browse files
mnlevy1981mnlevy1981
authored
Merge pull request #102 from mnlevy1981/test_StokesMOST_edgecases
Update the kpp test to include running with StokesMOST
2 parents ff9357e + 529b922 commit 65ef5c7

6 files changed

Lines changed: 222 additions & 26 deletions

File tree

‎bld/CompileFlags.mak‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ ifeq ($(FC),ftn)
66
FC_TMP = pgf90
77
endif
88
ifeq ($(PE_ENV),INTEL)
9-
FC_TMP = ifort
9+
FC_TMP = ifx
1010
endif
1111
ifeq ($(PE_ENV),PATHSCALE)
1212
FC_TMP = pathf95
@@ -27,8 +27,8 @@ ifeq ($(FC_TMP),pgf90)
2727
FCFLAGS = -O2 -Mfree -module $(OBJ_DIR)
2828
endif
2929

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

3434
ifeq ($(FC_TMP),xlf90)

‎reg_tests/kpp/input.nl‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,3 +52,13 @@ mix_type = 'kpp'
5252
Cp0 = 3992.0d0
5353
OBL_depth6 = 6000.0d0
5454
/
55+
56+
! Test 7 params
57+
&kpp_col7_nml
58+
ltest7 = .true.
59+
nlev7 = 10
60+
layer_thick7 = 5.0d0
61+
hmix7 = 17.0d0
62+
! Parameter settings to match LMD94 (linear interp, average Nsqr)
63+
interp_type_t7 = "linear"
64+
/

‎src/Makefile‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -90,7 +90,12 @@ ifeq ($(USE_NETCDF),TRUE)
9090
LINKING_FLAGS += -L$(NETCDF_LIB) -lnetcdff -lnetcdf
9191
else
9292
# get flibs from netcdf.pc
93-
FLIBS = $(subst flibs=,,$(shell grep flibs $(NETCDF_PC)))
93+
# older versions of netCDF have "flibs=" line
94+
FLIBS = $(shell grep flibs $(NETCDF_PC) | cut -d '=' -f 2-)
95+
ifeq ($(wildcard $(FLIBS)),)
96+
# newer versions of netCDF have Libs: line (and use {libdir} in output
97+
FLIBS = $(shell grep Libs $(NETCDF_PC) | grep -v private | cut -d '}' -f 2-)
98+
endif
9499
# Workaround for yellowstone, I need to figure out where this
95100
# @NC_FLIBS@ comes from
96101
ifeq ($(FLIBS),@NC_FLIBS@)
@@ -145,7 +150,7 @@ endif
145150

146151
all: exe
147152

148-
$(EXE): $(addprefix $(OBJ_DIR)/,$(OBJS) $(DRIVER_OBJS)) \
153+
$(EXE): $(addprefix $(OBJ_DIR)/,$(DRIVER_OBJS) $(OBJS)) \
149154
$(LIB_DIR)/libcvmix.a
150155
$(FC) -o $(EXE) $(addprefix $(OBJ_DIR)/,$(OBJS) $(DRIVER_OBJS)) $(LINKING_FLAGS)
151156

‎src/cvmix_driver.F90‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,12 @@ Program cvmix_driver
3636
real(kind=cvmix_r8) :: ocn_depth
3737
character(len=cvmix_strlen) :: mix_type
3838

39+
external cvmix_BL_driver
40+
external cvmix_shear_driver
41+
external cvmix_tidal_driver
42+
external cvmix_ddiff_driver
43+
external cvmix_kpp_driver
44+
3945
namelist/cvmix_nml/mix_type, nlev, max_nlev, ocn_depth
4046

4147
mix_type = 'unset'

‎src/drivers/cvmix_kpp_drv.F90‎

Lines changed: 179 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ Subroutine cvmix_kpp_driver()
2424
cvmix_kpp_compute_unresolved_shear, &
2525
cvmix_kpp_compute_turbulent_scales, &
2626
cvmix_kpp_compute_shape_function_coeffs, &
27+
cvmix_kpp_compute_StokesXi, &
2728
cvmix_coeffs_kpp
2829
use cvmix_put_get, only : cvmix_put
2930
use cvmix_io, only : cvmix_io_open, &
@@ -41,7 +42,7 @@ Subroutine cvmix_kpp_driver()
4142
real(cvmix_r8), parameter :: third = cvmix_one / real(3,cvmix_r8)
4243

4344
! CVMix datatypes
44-
type(cvmix_data_type) :: CVmix_vars1, CVmix_vars4, CVmix_vars5
45+
type(cvmix_data_type) :: CVmix_vars1, CVmix_vars4, CVmix_vars5, CVmix_vars7
4546

4647
real(cvmix_r8), dimension(:), allocatable, target :: Mdiff, Tdiff, Sdiff
4748
real(cvmix_r8), dimension(:), allocatable, target :: zt, zw_iface, &
@@ -52,19 +53,25 @@ Subroutine cvmix_kpp_driver()
5253
delta_vel_sqr, &
5354
buoy_freq_iface
5455
real(cvmix_r8), dimension(:,:), allocatable, target :: hor_vel
56+
real(cvmix_r8), dimension(:), allocatable :: uS, vS, uSbar, vSbar
57+
real(cvmix_r8), dimension(:), allocatable :: zeros, ones
5558
real(cvmix_r8), dimension(2) :: ref_vel
5659
real(cvmix_r8), dimension(4) :: shape_coeffs
57-
integer :: i, fid, kt, kw, nlev1, nlev3, nlev4, max_nlev4, OBL_levid4, nlev5
58-
real(cvmix_r8) :: hmix1, hmix5, ri_crit, layer_thick1, layer_thick4, &
59-
layer_thick5, OBL_depth4, OBL_depth5, N, Nsqr
60+
integer :: i, fid, kt, kw, nlev1, nlev3, nlev4, max_nlev4, OBL_levid4, &
61+
nlev5, nlev7
62+
real(cvmix_r8) :: hmix1, hmix5, hmix7, ri_crit, layer_thick1, &
63+
layer_thick4, layer_thick5, layer_thick7, OBL_depth4, &
64+
OBL_depth5, OBL_depth7, N, Nsqr
65+
real(cvmix_r8) :: StokesXI
6066
real(cvmix_r8) :: kOBL_depth, Bslope, Vslope
6167
real(cvmix_r8) :: sigma6, OBL_depth6, surf_buoy_force6, surf_fric_vel6, &
6268
vonkarman6, tao, rho0, grav, alpha, Qnonpen, Cp0, &
6369
w_m6, w_s6, wm6_true, ws6_true
64-
character(len=cvmix_strlen) :: interp_type_t1, interp_type_t4, interp_type_t5
70+
character(len=cvmix_strlen) :: interp_type_t1, interp_type_t4, &
71+
interp_type_t5, interp_type_t7
6572
character(len=cvmix_strlen) :: filename
6673
! True => run specified test
67-
logical :: ltest1, ltest2, ltest3, ltest4, ltest5, ltest6
74+
logical :: ltest1, ltest2, ltest3, ltest4, ltest5, ltest6, ltest7
6875
logical :: lnoDGat1 ! True => G'(1) = 0 (in test 4)
6976

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

7987
! Read namelists
8088

@@ -117,12 +125,20 @@ Subroutine cvmix_kpp_driver()
117125
Cp0 = real(3992, cvmix_r8)
118126
OBL_depth6 = real(6000, cvmix_r8)
119127

128+
! Defaults for test 7
129+
ltest7 = .false.
130+
nlev7 = 10
131+
layer_thick7 = real(5, cvmix_r8)
132+
hmix7 = real(17, cvmix_r8)
133+
interp_type_t7 = "linear"
134+
120135
read(*, nml=kpp_col1_nml)
121136
read(*, nml=kpp_col2_nml)
122137
read(*, nml=kpp_col3_nml)
123138
read(*, nml=kpp_col4_nml)
124139
read(*, nml=kpp_col5_nml)
125140
read(*, nml=kpp_col6_nml)
141+
read(*, nml=kpp_col7_nml)
126142

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

591607
end if ! ltest6
592608

609+
! Test 7: Test 5, but lStokesMOST = .true.
610+
if (ltest7) then
611+
print*, ""
612+
print*, "Test 7: Computing Bulk Richardson number (with StokesMOST)"
613+
print*, "----------"
614+
615+
! using linear interpolation, averaging Nsqr, and setting Cv = 1.5 to
616+
! match LMD result
617+
call cvmix_init_kpp(Cv = 1.5_cvmix_r8, interp_type=interp_type_t7, &
618+
lStokesMOST=.true., lMonOb=.true.)
619+
620+
! Set up vertical levels (centers and interfaces) and compute bulk
621+
! Richardson number
622+
allocate(zt(nlev7), zw_iface(nlev7+1))
623+
do kw=1,nlev7+1
624+
zw_iface(kw) = -layer_thick7*real(kw-1, cvmix_r8)
625+
end do
626+
do kt=1,nlev7
627+
zt(kt) = 0.5_cvmix_r8 * (zw_iface(kt)+zw_iface(kt+1))
628+
end do
629+
630+
! Compute Br-B(d), |Vr-V(d)|^2, and Vt^2
631+
allocate(buoyancy(nlev7), delta_vel_sqr(nlev7), hor_vel(nlev7,2), &
632+
shear_sqr(nlev7), w_s(nlev7), Ri_bulk(nlev7), Ri_bulk2(nlev7), &
633+
buoy_freq_iface(nlev7+1))
634+
635+
ref_vel(1) = 0.1_cvmix_r8
636+
ref_vel(2) = cvmix_zero
637+
N = 0.01_cvmix_r8
638+
Nsqr = N*N
639+
Bslope = -Nsqr
640+
Vslope = -0.1_cvmix_r8 / (real(nlev7,cvmix_r8)*layer_thick7-hmix7)
641+
do kt=1,nlev7
642+
if ((zt(kt).ge.-hmix7).or.(kt.eq.1)) then
643+
buoyancy(kt) = Nsqr
644+
hor_vel(kt,1) = 0.1_cvmix_r8
645+
buoy_freq_iface(kt) = cvmix_zero
646+
else
647+
if (zw_iface(kt).ge.-hmix7) then
648+
! derivatives of buoyancy and horizontal velocity component are
649+
! discontinuous in this layer (no change -> non-zero linear change)
650+
! so we compute area-average of analytic function over layer
651+
buoyancy(kt) = Bslope*(-zw_iface(kt+1)-real(hmix7,cvmix_r8))**2 / &
652+
(real(2,cvmix_r8)*layer_thick7) + Nsqr
653+
hor_vel(kt,1) = Vslope*(-zw_iface(kt+1)-real(hmix7,cvmix_r8))**2 / &
654+
(real(2,cvmix_r8)*layer_thick7) + 0.1_cvmix_r8
655+
else
656+
buoyancy(kt) = Nsqr+Bslope*(-zt(kt)-real(hmix7,cvmix_r8))
657+
hor_vel(kt,1) = 0.1_cvmix_r8 + Vslope * &
658+
(-zt(kt)-real(hmix7,cvmix_r8))
659+
end if
660+
buoy_freq_iface(kt) = sqrt(-(buoyancy(kt)-buoyancy(kt-1)) / &
661+
layer_thick7)
662+
end if
663+
! Compute w_s with zeta=0 per LMD page 393
664+
! => w_s = von Karman * surf_fric_vel = 0.4*0.01 = 4e-3
665+
call cvmix_kpp_compute_turbulent_scales(cvmix_zero, -zt(kt), &
666+
buoyancy(1), 0.01_cvmix_r8, &
667+
w_s=w_s(kt))
668+
hor_vel(kt,2) = cvmix_zero
669+
delta_vel_sqr(kt) = (ref_vel(1)-hor_vel(kt,1))**2 + &
670+
(ref_vel(2)-hor_vel(kt,2))**2
671+
end do
672+
buoy_freq_iface(nlev7+1) = N
673+
! MNL: test both uses of compute_bulk_Richardson
674+
Ri_bulk = cvmix_kpp_compute_bulk_Richardson(zt, (buoyancy(1)-buoyancy), &
675+
delta_vel_sqr, &
676+
Nsqr_iface = buoy_freq_iface**2, &
677+
ws_cntr = w_s)
678+
679+
shear_sqr = cvmix_kpp_compute_unresolved_shear(zt, w_s, &
680+
Nsqr_iface = buoy_freq_iface**2)
681+
! Note that Vt_shear_sqr is the fourth argument in compute_bulk_Richardson
682+
! so it does not need to declared explicitly (even though it is optional)
683+
Ri_bulk2 = cvmix_kpp_compute_bulk_Richardson(zt, (buoyancy(1)-buoyancy), &
684+
delta_vel_sqr, shear_sqr)
685+
allocate(zeros(nlev7), source=cvmix_zero)
686+
allocate(ones(nlev7), &
687+
uS(nlev7), &
688+
vS(nlev7), &
689+
uSbar(nlev7), &
690+
vSbar(nlev7), &
691+
source=cvmix_one)
692+
StokesXI = cvmix_one
693+
call cvmix_kpp_compute_StokesXi(zi=zt, &
694+
zk=zw_iface, &
695+
kSL=(nlev7+1)/2, &
696+
SLDepth=cvmix_zero, &
697+
surf_buoy_force=cvmix_one, &
698+
surf_fric_vel=cvmix_one, &
699+
omega_w2x=cvmix_one, &
700+
uE=ones, &
701+
vE=ones, &
702+
uS=uS, &
703+
vS=vS, &
704+
uSbar=uSbar, &
705+
vSbar=vSbar, &
706+
uS_SLD=cvmix_one, &
707+
vS_SLD=cvmix_one, &
708+
uSbar_SLD=cvmix_one, &
709+
vSbar_SLD=cvmix_one, &
710+
StokesXI=StokesXI)
711+
call cvmix_kpp_compute_OBL_depth(Ri_bulk, zw_iface, OBL_depth7, &
712+
kOBL_depth, zt, surf_fric=cvmix_one, &
713+
surf_buoy=ones, Xi=ones)
714+
deallocate(zeros, ones, uS, vS, uSbar, vSbar)
715+
do kt=1,nlev7
716+
if (abs(Ri_bulk(kt)-Ri_bulk2(kt)).gt.1e-12_cvmix_r8) then
717+
print*, "WARNING: two Ri_bulk computations did not match!"
718+
print*, zt(kt), Ri_bulk(kt), Ri_bulk2(kt)
719+
else
720+
print*, zt(kt), Ri_bulk(kt)
721+
end if
722+
end do
723+
print*, "OBL has depth of ", OBL_depth7
724+
#ifdef _NETCDF
725+
print*, "Done! Data is stored in test7.nc, run plot_bulk_Rich.ncl ", &
726+
"to see output."
727+
#else
728+
print*, "Done! Data is stored in test7.out, run plot_bulk_Rich.ncl ", &
729+
"to see output."
730+
#endif
731+
732+
CVmix_vars7%nlev = nlev7
733+
CVmix_vars7%BoundaryLayerDepth = OBL_depth7
734+
CVmix_vars7%zt_cntr => zt
735+
CVmix_vars7%BulkRichardson_cntr => Ri_bulk
736+
CVmix_vars7%Vx_cntr => hor_vel(:,1)
737+
#ifdef _NETCDF
738+
call cvmix_io_open(fid, "test7.nc", "nc")
739+
#else
740+
call cvmix_io_open(fid, "test7.out", "ascii")
741+
#endif
742+
call cvmix_output_write(fid, CVmix_vars7, (/"zt ", &
743+
"Ri_bulk ", &
744+
"Vx ", &
745+
"buoyancy"/), buoyancy)
746+
#ifdef _NETCDF
747+
call cvmix_output_write_att(fid, "OBL_depth", &
748+
CVmix_vars7%BoundaryLayerDepth)
749+
call cvmix_output_write_att(fid, "longname", "ocean height (cell center)",&
750+
"zt")
751+
call cvmix_output_write_att(fid, "units", "m", "zt")
752+
call cvmix_output_write_att(fid, "longname", "horizontal velocity", "U")
753+
call cvmix_output_write_att(fid, "units", "m/s", "U")
754+
call cvmix_output_write_att(fid, "units", "m/s^2", "buoyancy")
755+
call cvmix_output_write_att(fid, "longname", "Bulk Richardson number", &
756+
"BulkRichardson")
757+
call cvmix_output_write_att(fid, "units", "unitless", "BulkRichardson")
758+
#endif
759+
call cvmix_io_close(fid)
760+
761+
deallocate(zt, zw_iface)
762+
deallocate(buoyancy, delta_vel_sqr, hor_vel, shear_sqr, w_s, Ri_bulk, &
763+
Ri_bulk2, buoy_freq_iface)
764+
end if ! ltest7
765+
593766
!EOC
594767

595768
End Subroutine cvmix_kpp_driver

‎src/shared/cvmix_kpp.F90‎

Lines changed: 17 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1576,7 +1576,7 @@ subroutine cvmix_kpp_compute_OBL_depth_low(Ri_bulk, zw_iface, OBL_depth, &
15761576
! Local variables
15771577
real(kind=cvmix_r8), dimension(:), pointer :: depth
15781578
real(kind=cvmix_r8), dimension(4) :: coeffs
1579-
real(kind=cvmix_r8) :: Ekman, MoninObukhov, OBL_Limit
1579+
real(kind=cvmix_r8) :: Ekman, MoninObukhov, OBL_Limit, numer, denom
15801580
integer :: nlev, k, kRi
15811581
logical :: lstable
15821582

@@ -1671,16 +1671,18 @@ subroutine cvmix_kpp_compute_OBL_depth_low(Ri_bulk, zw_iface, OBL_depth, &
16711671
if (CVmix_kpp_params_in%lMonOb ) then
16721672
if ( present(Xi) .and. present(surf_buoy) ) then
16731673
MoninObukhov = OBL_limit
1674+
numer = surf_fric**3
16741675
do k = 0, kRi-1
1675-
if (surf_buoy(k+1) .gt. cvmix_zero) MoninObukhov = &
1676-
surf_fric**3 / (surf_buoy(k+1) * (cvmix_one-Xi(k+1)))
1676+
denom = surf_buoy(k+1) * (cvmix_one-Xi(k+1))
1677+
if ( denom*OBL_limit .gt. numer ) MoninObukhov = &
1678+
numer / denom
16771679
if ( MoninObukhov .lt. abs(zt_cntr(k+1)) ) &
16781680
exit
16791681
end do
16801682
if (k.eq.0) then
16811683
OBL_limit = abs(zt_cntr(1))
16821684
elseif (k.lt.kRi) then
1683-
OBL_limit = min( OBL_limit, abs(zw_iface(k)) )
1685+
OBL_limit = min( OBL_limit, abs(zw_iface(k+1)) )
16841686
end if
16851687

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

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

16961697
else ! not Stokes_MOST
16971698

@@ -1763,13 +1764,13 @@ subroutine cvmix_kpp_compute_OBL_depth_low(Ri_bulk, zw_iface, OBL_depth, &
17631764
! because we know OBL_depth will equal OBL_limit?
17641765
OBL_depth = min(OBL_depth, OBL_limit)
17651766
end if
1767+
end if ! lStokesMOST
17661768

1767-
OBL_depth = max(OBL_depth, CVmix_kpp_params_in%minOBLdepth)
1768-
if (CVmix_kpp_params_in%maxOBLdepth.gt.cvmix_zero) &
1769-
OBL_depth = min(OBL_depth, CVmix_kpp_params_in%maxOBLdepth)
1770-
kOBL_depth = cvmix_kpp_compute_kOBL_depth(zw_iface, zt_cntr, OBL_depth)
1769+
OBL_depth = max(OBL_depth, CVmix_kpp_params_in%minOBLdepth)
1770+
if (CVmix_kpp_params_in%maxOBLdepth.gt.cvmix_zero) &
1771+
OBL_depth = min(OBL_depth, CVmix_kpp_params_in%maxOBLdepth)
1772+
kOBL_depth = cvmix_kpp_compute_kOBL_depth(zw_iface, zt_cntr, OBL_depth)
17711773

1772-
end if ! lStokesMOST
17731774

17741775
!EOC
17751776

@@ -3344,12 +3345,13 @@ subroutine cvmix_kpp_composite_Gshape(sigma , Gat1, Gsig, dGdsig)
33443345
!EOP
33453346
!BOC
33463347

3347-
real(cvmix_r8) :: a2Gsig, a3Gsig, bGsig, sig_m, G_m, G_1, sig
3348+
real(cvmix_r8), parameter :: a2Gsig = -2.1637_cvmix_r8, &
3349+
a3Gsig = 0.5831_cvmix_r8, &
3350+
sig_m = 0.35_cvmix_r8, &
3351+
G_m = 0.11_cvmix_r8
3352+
3353+
real(cvmix_r8) :: bGsig, G_1, sig
33483354

3349-
a2Gsig = -2.1637_cvmix_r8
3350-
a3Gsig = 0.5831_cvmix_r8
3351-
sig_m = 0.35_cvmix_r8
3352-
G_m = 0.11_cvmix_r8 ! sig_m + sig_m * sig_m * (a2Gsig + a3Gsig * sig_m)
33533355
G_1 = MAX( cvmix_zero , MIN( Gat1 , G_m ) )
33543356

33553357
if (sigma .lT. sig_m) then

0 commit comments

Comments
 (0)