Skip to content

Commit

Permalink
Merge branch 'dev/gfdl' into restore-ZB20
Browse files Browse the repository at this point in the history
  • Loading branch information
@marshallward
marshallward committed Aug 21, 2024
2 parents d86630e + 2024f6f commit 855119e
Show file tree
Hide file tree
Showing 3 changed files with 65 additions and 31 deletions.
54 changes: 40 additions & 14 deletions src/core/MOM_barotropic.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ module MOM_barotropic
use MOM_grid, only : ocean_grid_type
use MOM_harmonic_analysis, only : HA_accum_FtSSH, harmonic_analysis_CS
use MOM_hor_index, only : hor_index_type
use MOM_io, only : vardesc, var_desc, MOM_read_data, slasher
use MOM_io, only : vardesc, var_desc, MOM_read_data, slasher, NORTH_FACE, EAST_FACE
use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, open_boundary_query
use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W
use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_segment_type
Expand Down Expand Up @@ -4459,6 +4459,10 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
! drag piston velocity.
character(len=80) :: wave_drag_var ! The wave drag piston velocity variable
! name in wave_drag_file.
character(len=80) :: wave_drag_u ! The wave drag piston velocity variable
! name in wave_drag_file.
character(len=80) :: wave_drag_v ! The wave drag piston velocity variable
! name in wave_drag_file.
real :: mean_SL ! The mean sea level that is used along with the bathymetry to estimate the
! geometry when LINEARIZED_BT_CORIOLIS is true or BT_NONLIN_STRESS is false [Z ~> m].
real :: Z_to_H ! A local unit conversion factor [H Z-1 ~> nondim or kg m-3]
Expand Down Expand Up @@ -4703,8 +4707,17 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
"piston velocities.", default="", do_not_log=.not.CS%linear_wave_drag)
call get_param(param_file, mdl, "BT_WAVE_DRAG_VAR", wave_drag_var, &
"The name of the variable in BT_WAVE_DRAG_FILE with the "//&
"barotropic linear wave drag piston velocities at h points.", &
"barotropic linear wave drag piston velocities at h points. "//&
"It will not be used if both BT_WAVE_DRAG_U and BT_WAVE_DRAG_V are defined.", &
default="rH", do_not_log=.not.CS%linear_wave_drag)
call get_param(param_file, mdl, "BT_WAVE_DRAG_U", wave_drag_u, &
"The name of the variable in BT_WAVE_DRAG_FILE with the "//&
"barotropic linear wave drag piston velocities at u points.", &
default="", do_not_log=.not.CS%linear_wave_drag)
call get_param(param_file, mdl, "BT_WAVE_DRAG_V", wave_drag_v, &
"The name of the variable in BT_WAVE_DRAG_FILE with the "//&
"barotropic linear wave drag piston velocities at v points.", &
default="", do_not_log=.not.CS%linear_wave_drag)
call get_param(param_file, mdl, "BT_WAVE_DRAG_SCALE", wave_drag_scale, &
"A scaling factor for the barotropic linear wave drag "//&
"piston velocities.", default=1.0, units="nondim", &
Expand Down Expand Up @@ -4924,19 +4937,32 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
wave_drag_file = trim(slasher(inputdir))//trim(wave_drag_file)
call log_param(param_file, mdl, "INPUTDIR/BT_WAVE_DRAG_FILE", wave_drag_file)

allocate(lin_drag_h(isd:ied,jsd:jed), source=0.0)
if (len_trim(wave_drag_u) > 0 .and. len_trim(wave_drag_v) > 0) then
call MOM_read_data(wave_drag_file, wave_drag_u, CS%lin_drag_u, G%Domain, &
position=EAST_FACE, scale=GV%m_to_H*US%T_to_s)
call pass_var(CS%lin_drag_u, G%Domain)
CS%lin_drag_u(:,:) = wave_drag_scale * CS%lin_drag_u(:,:)

call MOM_read_data(wave_drag_file, wave_drag_var, lin_drag_h, G%Domain, scale=GV%m_to_H*US%T_to_s)
call pass_var(lin_drag_h, G%Domain)
do j=js,je ; do I=is-1,ie
CS%lin_drag_u(I,j) = wave_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i+1,j))
enddo ; enddo
do J=js-1,je ; do i=is,ie
CS%lin_drag_v(i,J) = wave_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i,j+1))
enddo ; enddo
deallocate(lin_drag_h)
endif
endif
call MOM_read_data(wave_drag_file, wave_drag_v, CS%lin_drag_v, G%Domain, &
position=NORTH_FACE, scale=GV%m_to_H*US%T_to_s)
call pass_var(CS%lin_drag_v, G%Domain)
CS%lin_drag_v(:,:) = wave_drag_scale * CS%lin_drag_v(:,:)

else
allocate(lin_drag_h(isd:ied,jsd:jed), source=0.0)

call MOM_read_data(wave_drag_file, wave_drag_var, lin_drag_h, G%Domain, scale=GV%m_to_H*US%T_to_s)
call pass_var(lin_drag_h, G%Domain)
do j=js,je ; do I=is-1,ie
CS%lin_drag_u(I,j) = wave_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i+1,j))
enddo ; enddo
do J=js-1,je ; do i=is,ie
CS%lin_drag_v(i,J) = wave_drag_scale * 0.5 * (lin_drag_h(i,j) + lin_drag_h(i,j+1))
enddo ; enddo
deallocate(lin_drag_h)
endif ! len_trim(wave_drag_u) > 0 .and. len_trim(wave_drag_v) > 0
endif ! len_trim(wave_drag_file) > 0
endif ! CS%linear_wave_drag

CS%dtbt_fraction = 0.98 ; if (dtbt_input < 0.0) CS%dtbt_fraction = -dtbt_input

Expand Down
32 changes: 19 additions & 13 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,6 +123,7 @@ module MOM_hor_visc
real :: min_grid_Ah !< Minimun horizontal biharmonic viscosity used to
!! limit grid Reynolds number [L4 T-1 ~> m4 s-1]
logical :: use_cont_thick !< If true, thickness at velocity points adopts h[uv] in BT_cont from continuity solver.
logical :: use_cont_thick_bug !< If true, retain an answer-changing bug for thickness at velocity points.
type(ZB2020_CS) :: ZB2020 !< Zanna-Bolton 2020 control structure.
logical :: use_ZB2020 !< If true, use Zanna-Bolton 2020 parameterization.

Expand Down Expand Up @@ -282,9 +283,9 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
type(thickness_diffuse_CS), optional, intent(in) :: TD !< Thickness diffusion control structure
type(accel_diag_ptrs), optional, intent(in) :: ADp !< Acceleration diagnostics
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
optional, intent(in) :: hu_cont !< Layer thickness at u-points [H ~> m or kg m-2].
optional, intent(inout) :: hu_cont !< Layer thickness at u-points [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
optional, intent(in) :: hv_cont !< Layer thickness at v-points [H ~> m or kg m-2].
optional, intent(inout) :: hv_cont !< Layer thickness at v-points [H ~> m or kg m-2].

! Local variables
real, dimension(SZIB_(G),SZJ_(G)) :: &
Expand Down Expand Up @@ -477,6 +478,9 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
use_MEKE_Au = allocated(MEKE%Au)

use_cont_huv = CS%use_cont_thick .and. present(hu_cont) .and. present(hv_cont)
if (use_cont_huv .and. .not.CS%use_cont_thick_bug) then
call pass_vector(hu_cont, hv_cont, G%domain, To_All+Scalar_Pair, halo=2)
endif

rescale_Kh = .false.
if (VarMix%use_variable_mixing) then
Expand Down Expand Up @@ -709,7 +713,14 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
! in OBCs, which are not ordinarily be necessary, and might not be necessary
! even with OBCs if the accelerations are zeroed at OBC points, in which
! case the j-loop for h_u could collapse to j=js=1,je+1. -RWH
if (CS%use_land_mask) then
if (use_cont_huv) then
do j=js-2,je+2 ; do I=Isq-1,Ieq+1
h_u(I,j) = hu_cont(I,j,k)
enddo ; enddo
do J=Jsq-1,Jeq+1 ; do i=is-2,ie+2
h_v(i,J) = hv_cont(i,J,k)
enddo ; enddo
elseif (CS%use_land_mask) then
do j=js-2,je+2 ; do I=is-2,Ieq+1
h_u(I,j) = 0.5 * (G%mask2dT(i,j)*h(i,j,k) + G%mask2dT(i+1,j)*h(i+1,j,k))
enddo ; enddo
Expand All @@ -725,16 +736,6 @@ subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G,
enddo ; enddo
endif

! The following should obviously be combined with the previous block if adopted.
if (use_cont_huv) then
do j=js-2,je+2 ; do I=Isq-1,Ieq+1
h_u(I,j) = hu_cont(I,j,k)
enddo ; enddo
do J=Jsq-1,Jeq+1 ; do i=is-2,ie+2
h_v(i,J) = hv_cont(i,J,k)
enddo ; enddo
endif

! Adjust contributions to shearing strain and interpolated values of
! thicknesses on open boundaries.
if (apply_OBC) then ; do n=1,OBC%number_of_segments
Expand Down Expand Up @@ -2140,6 +2141,11 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
call get_param(param_file, mdl, "USE_CONT_THICKNESS", CS%use_cont_thick, &
"If true, use thickness at velocity points from continuity solver. This option "//&
"currently only works with split mode.", default=.false.)
call get_param(param_file, mdl, "USE_CONT_THICKNESS_BUG", CS%use_cont_thick_bug, &
"If true, retain an answer-changing halo update bug when "//&
"USE_CONT_THICKNESS=True. This is not recommended.", &
default=.false., do_not_log=.not.CS%use_cont_thick)

call get_param(param_file, mdl, "LAPLACIAN", CS%Laplacian, &
"If true, use a Laplacian horizontal viscosity.", &
default=.false.)
Expand Down
10 changes: 6 additions & 4 deletions src/parameterizations/lateral/MOM_mixed_layer_restrat.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1064,8 +1064,8 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
h_big = 0.5*( big_H(i,j) + big_H(i+1,j) ) ! H ~> m or kg m-3
grd_b = ( buoy_av(i+1,j) - buoy_av(i,j) ) * G%IdxCu(I,j) ! L H-1 T-2 ~> s-2 or m3 kg-1 s-2
r_wpup = 2. / ( wpup(i,j) + wpup(i+1,j) ) ! T2 L-1 H-1 ~> s2 m-2 or m s2 kg-1
psi_mag = ( ( ( CS%Cr_space(i,j) * grid_dsd ) * ( absf * h_sml ) ) & ! L2 H T-1 ~> m3 s-1 or kg s-1
* ( ( h_big**2 ) * grd_b ) ) * r_wpup
psi_mag = ( ( ( (0.5*(CS%Cr_space(i,j) + CS%Cr_space(i+1,j))) * grid_dsd ) & ! L2 H T-1 ~> m3 s-1 or kg s-1
* ( absf * h_sml ) ) * ( ( h_big**2 ) * grd_b ) ) * r_wpup
else ! There is no flux on land and no gradient at open boundary points.
psi_mag = 0.0
endif
Expand Down Expand Up @@ -1105,8 +1105,8 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d
h_big = 0.5*( big_H(i,j) + big_H(i,j+1) ) ! H ~> m or kg m-3
grd_b = ( buoy_av(i,j+1) - buoy_av(i,j) ) * G%IdyCv(I,j) ! L H-1 T-2 ~> s-2 or m3 kg-1 s-2
r_wpup = 2. / ( wpup(i,j) + wpup(i,j+1) ) ! T2 L-1 H-1 ~> s2 m-2 or m s2 kg-1
psi_mag = ( ( ( CS%Cr_space(i,j) * grid_dsd ) * ( absf * h_sml ) ) & ! L2 H T-1 ~> m3 s-1 or kg s-1
* ( ( h_big**2 ) * grd_b ) ) * r_wpup
psi_mag = ( ( ( (0.5*(CS%Cr_space(i,j) + CS%Cr_space(i,j+1))) * grid_dsd ) & ! L2 H T-1 ~> m3 s-1 or kg s-1
* ( absf * h_sml ) ) * ( ( h_big**2 ) * grd_b ) ) * r_wpup
else ! There is no flux on land and no gradient at open boundary points.
psi_mag = 0.0
endif
Expand Down Expand Up @@ -1670,6 +1670,7 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
filename = trim(inputdir) // "/" // trim(filename)
allocate(CS%MLD_Tfilt_space(G%isd:G%ied,G%jsd:G%jed), source=0.0)
call MOM_read_data(filename, varname, CS%MLD_Tfilt_space, G%domain, scale=US%s_to_T)
call pass_var(CS%MLD_Tfilt_space, G%domain)
endif
allocate(CS%Cr_space(G%isd:G%ied,G%jsd:G%jed), source=CS%Cr)
if (CS%Cr_grid) then
Expand All @@ -1681,6 +1682,7 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS,
default="Cr")
filename = trim(inputdir) // "/" // trim(filename)
call MOM_read_data(filename, varname, CS%Cr_space, G%domain)
call pass_var(CS%Cr_space, G%domain)
endif
call closeParameterBlock(param_file) ! The remaining parameters do not have MLE% prepended
call get_param(param_file, mdl, "MLE_USE_PBL_MLD", CS%MLE_use_PBL_MLD, &
Expand Down

0 comments on commit 855119e

Please sign in to comment.