Skip to content

Commit

Permalink
Merge branch 'dev/gfdl' into add_check_scaling
Browse files Browse the repository at this point in the history
  • Loading branch information
@marshallward
marshallward authored Jan 29, 2022
2 parents 26216f6 + 65998cd commit 4b8645f
Show file tree
Hide file tree
Showing 8 changed files with 132 additions and 69 deletions.
2 changes: 1 addition & 1 deletion docs/conf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,7 +159,7 @@ def latexPassthru(name, rawtext, text, lineno, inliner, options={}, content=[]):

# General information about the project.
project = u'MOM6'
copyright = u'2017-2021, MOM6 developers'
copyright = u'2017-2022, MOM6 developers'

# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
Expand Down
Binary file added docs/images/channel_drag.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
56 changes: 30 additions & 26 deletions src/ALE/MOM_ALE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -378,7 +378,7 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h)
call diag_update_remap_grids(CS%diag)
endif
! Remap all variables from old grid h onto new grid h_new
call remap_all_state_vars( CS%remapCS, CS, G, GV, h, h_new, Reg, OBC, -dzRegrid, &
call remap_all_state_vars( CS%remapCS, CS, G, GV, h, h_new, Reg, OBC, dzRegrid, &
u, v, CS%show_call_tree, dt )

if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_main)")
Expand Down Expand Up @@ -732,10 +732,10 @@ end subroutine ALE_regrid_accelerated
!! new grids. When velocity components need to be remapped, thicknesses at
!! velocity points are taken to be arithmetic averages of tracer thicknesses.
!! This routine is called during initialization of the model at time=0, to
!! remap initiali conditions to the model grid. It is also called during a
!! remap initial conditions to the model grid. It is also called during a
!! time step to update the state.
subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg, OBC, &
dxInterface, u, v, debug, dt)
dzInterface, u, v, debug, dt)
type(remapping_CS), intent(in) :: CS_remapping !< Remapping control structure
type(ALE_CS), intent(in) :: CS_ALE !< ALE control structure
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
Expand All @@ -747,37 +747,42 @@ subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg,
type(tracer_registry_type), pointer :: Reg !< Tracer registry structure
type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
optional, intent(in) :: dxInterface !< Change in interface position
optional, intent(in) :: dzInterface !< Change in interface position
!! [H ~> m or kg m-2]
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
optional, intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
optional, intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
logical, optional, intent(in) :: debug !< If true, show the call tree
real, optional, intent(in) :: dt !< time step for diagnostics [T ~> s]

! Local variables
integer :: i, j, k, m
integer :: nz, ntr
real, dimension(GV%ke+1) :: dx
real, dimension(GV%ke) :: h1, u_column
real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: work_conc
real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: work_cont
real, dimension(SZI_(G), SZJ_(G)) :: work_2d
real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
real :: ppt2mks
real, dimension(GV%ke) :: h2
real :: h_neglect, h_neglect_edge
real, dimension(GV%ke+1) :: dz ! The change in interface heights interpolated to
! a velocity point [H ~> m or kg m-2]
real, dimension(GV%ke) :: h1 ! A column of initial thicknesses [H ~> m or kg m-2]
real, dimension(GV%ke) :: h2 ! A column of updated thicknesses [H ~> m or kg m-2]
real, dimension(GV%ke) :: u_column ! A column of properties, like tracer concentrations
! or velocities, being remapped [various units]
real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: work_conc ! The rate of change of concentrations [Conc T-1 ~> Conc s-1]
real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: work_cont ! The rate of change of cell-integrated tracer
! content [Conc H T-1 ~> Conc m s-1 or Conc kg m-2 s-1] or
! cell thickness [H T-1 ~> m s-1 or Conc kg m-2 s-1]
real, dimension(SZI_(G), SZJ_(G)) :: work_2d ! The rate of change of column-integrated tracer
! content [Conc H T-1 ~> Conc m s-1 or Conc kg m-2 s-1]
real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
real :: h_neglect, h_neglect_edge ! Tiny thicknesses used in remapping [H ~> m or kg m-2]
logical :: show_call_tree
type(tracer_type), pointer :: Tr => NULL()
integer :: i, j, k, m, nz, ntr

show_call_tree = .false.
if (present(debug)) show_call_tree = debug
if (show_call_tree) call callTree_enter("remap_all_state_vars(), MOM_ALE.F90")

! If remap_uv_using_old_alg is .true. and u or v is requested, then we must have dxInterface. Otherwise,
! u and v can be remapped without dxInterface
if ( .not. present(dxInterface) .and. (CS_ALE%remap_uv_using_old_alg .and. (present(u) .or. present(v))) ) then
call MOM_error(FATAL, "remap_all_state_vars: dxInterface must be present if using old algorithm "// &
! If remap_uv_using_old_alg is .true. and u or v is requested, then we must have dzInterface. Otherwise,
! u and v can be remapped without dzInterface
if ( .not. present(dzInterface) .and. (CS_ALE%remap_uv_using_old_alg .and. (present(u) .or. present(v))) ) then
call MOM_error(FATAL, "remap_all_state_vars: dzInterface must be present if using old algorithm "// &
"and u/v are to be remapped")
endif

Expand All @@ -790,7 +795,6 @@ subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg,
endif

nz = GV%ke
ppt2mks = 0.001

ntr = 0 ; if (associated(Reg)) ntr = Reg%ntr

Expand Down Expand Up @@ -856,14 +860,14 @@ subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg,

! Remap u velocity component
if ( present(u) ) then
!$OMP parallel do default(shared) private(h1,h2,dx,u_column)
!$OMP parallel do default(shared) private(h1,h2,dz,u_column)
do j = G%jsc,G%jec ; do I = G%iscB,G%iecB ; if (G%mask2dCu(I,j)>0.) then
! Build the start and final grids
h1(:) = 0.5 * ( h_old(i,j,:) + h_old(i+1,j,:) )
if (CS_ALE%remap_uv_using_old_alg) then
dx(:) = 0.5 * ( dxInterface(i,j,:) + dxInterface(i+1,j,:) )
dz(:) = 0.5 * ( dzInterface(i,j,:) + dzInterface(i+1,j,:) )
do k = 1, nz
h2(k) = max( 0., h1(k) + ( dx(k+1) - dx(k) ) )
h2(k) = max( 0., h1(k) + ( dz(k) - dz(k+1) ) )
enddo
else
h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i+1,j,:) )
Expand All @@ -889,14 +893,14 @@ subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg,

! Remap v velocity component
if ( present(v) ) then
!$OMP parallel do default(shared) private(h1,h2,dx,u_column)
!$OMP parallel do default(shared) private(h1,h2,dz,u_column)
do J = G%jscB,G%jecB ; do i = G%isc,G%iec ; if (G%mask2dCv(i,j)>0.) then
! Build the start and final grids
h1(:) = 0.5 * ( h_old(i,j,:) + h_old(i,j+1,:) )
if (CS_ALE%remap_uv_using_old_alg) then
dx(:) = 0.5 * ( dxInterface(i,j,:) + dxInterface(i,j+1,:) )
dz(:) = 0.5 * ( dzInterface(i,j,:) + dzInterface(i,j+1,:) )
do k = 1, nz
h2(k) = max( 0., h1(k) + ( dx(k+1) - dx(k) ) )
h2(k) = max( 0., h1(k) + ( dz(k) - dz(k+1) ) )
enddo
else
h2(:) = 0.5 * ( h_new(i,j,:) + h_new(i,j+1,:) )
Expand Down
65 changes: 33 additions & 32 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1003,6 +1003,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
visc_bound_rem(i,j) = 0.0
Kh(i,j) = hrat_min(i,j) * CS%Kh_Max_xx(i,j)
else
! ### NOTE: The denominator could be zero here - AJA ###
visc_bound_rem(i,j) = 1.0 - Kh(i,j) / (hrat_min(i,j) * CS%Kh_Max_xx(i,j))
endif
enddo ; enddo
Expand Down Expand Up @@ -1194,7 +1195,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
if (CS%better_bound_Ah .or. CS%better_bound_Kh) then
do J=js-1,Jeq ; do I=is-1,Ieq
h_min = min(h_u(I,j), h_u(I,j+1), h_v(i,J), h_v(i+1,J))
hrat_min(i,j) = min(1.0, h_min / (hq(I,J) + h_neglect))
hrat_min(I,J) = min(1.0, h_min / (hq(I,J) + h_neglect))
enddo ; enddo

if (CS%better_bound_Kh) then
Expand All @@ -1217,11 +1218,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
(G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) == 0.0) then
! Only one of hu and hv is nonzero, so just add them.
hq(I,J) = hu + hv
hrat_min(i,j) = 1.0
hrat_min(I,J) = 1.0
else
! Both hu and hv are nonzero, so take the harmonic mean.
hq(I,J) = 2.0 * (hu * hv) / ((hu + hv) + h_neglect)
hrat_min(i,j) = min(1.0, min(hu, hv) / (hq(I,J) + h_neglect) )
hrat_min(I,J) = min(1.0, min(hu, hv) / (hq(I,J) + h_neglect) )
endif
endif
endif
Expand All @@ -1234,11 +1235,11 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
do J=js-1,Jeq ; do I=is-1,Ieq
grad_vort = grad_vort_mag_q(I,J) + grad_div_mag_q(I,J)
grad_vort_qg = 3. * grad_vort_mag_q_2d(I,J)
vert_vort_mag(i,j) = min(grad_vort, grad_vort_qg)
vert_vort_mag(I,J) = min(grad_vort, grad_vort_qg)
enddo ; enddo
else
do J=js-1,Jeq ; do I=is-1,Ieq
vert_vort_mag(i,j) = grad_vort_mag_q(I,J) + grad_div_mag_q(I,J)
vert_vort_mag(I,J) = grad_vort_mag_q(I,J) + grad_div_mag_q(I,J)
enddo ; enddo
endif
endif
Expand All @@ -1254,23 +1255,23 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
if (CS%Smagorinsky_Kh) then
if (CS%add_LES_viscosity) then
do J=js-1,Jeq ; do I=is-1,Ieq
Kh(i,j) = Kh(i,j) + CS%Laplac2_const_xx(i,j) * Shear_mag(i,j)
Kh(I,J) = Kh(I,J) + CS%Laplac2_const_xx(i,j) * Shear_mag(i,j)
enddo ; enddo
else
do J=js-1,Jeq ; do I=is-1,Ieq
Kh(i,j) = max(Kh(i,j), CS%Laplac2_const_xy(I,J) * Shear_mag(i,j) )
Kh(I,J) = max(Kh(I,J), CS%Laplac2_const_xy(I,J) * Shear_mag(i,j) )
enddo ; enddo
endif
endif

if (CS%Leith_Kh) then
if (CS%add_LES_viscosity) then
do J=js-1,Jeq ; do I=is-1,Ieq
Kh(i,j) = Kh(i,j) + CS%Laplac3_const_xx(i,j) * vert_vort_mag(i,j) * inv_PI3
Kh(I,J) = Kh(I,J) + CS%Laplac3_const_xx(i,j) * vert_vort_mag(I,J) * inv_PI3 ! Is this right? -AJA
enddo ; enddo
else
do J=js-1,Jeq ; do I=is-1,Ieq
Kh(i,j) = max(Kh(i,j), CS%Laplac3_const_xy(I,J) * vert_vort_mag(i,j) * inv_PI3)
Kh(I,J) = max(Kh(I,J), CS%Laplac3_const_xy(I,J) * vert_vort_mag(I,J) * inv_PI3)
enddo ; enddo
endif
endif
Expand All @@ -1281,40 +1282,40 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
! stack size has been reduced.
do J=js-1,Jeq ; do I=is-1,Ieq
if (rescale_Kh) &
Kh(i,j) = VarMix%Res_fn_q(i,j) * Kh(i,j)
Kh(I,J) = VarMix%Res_fn_q(I,J) * Kh(I,J)

if (CS%res_scale_MEKE) &
meke_res_fn = VarMix%Res_fn_q(i,j)
meke_res_fn = VarMix%Res_fn_q(I,J)

! Older method of bounding for stability
if (legacy_bound) &
Kh(i,j) = min(Kh(i,j), CS%Kh_Max_xy(i,j))
Kh(I,J) = min(Kh(I,J), CS%Kh_Max_xy(I,J))

Kh(i,j) = max(Kh(i,j), CS%Kh_bg_min) ! Place a floor on the viscosity, if desired.
Kh(I,J) = max(Kh(I,J), CS%Kh_bg_min) ! Place a floor on the viscosity, if desired.

if (use_MEKE_Ku) then
! *Add* the MEKE contribution (might be negative)
Kh(i,j) = Kh(i,j) + 0.25*( (MEKE%Ku(i,j) + MEKE%Ku(i+1,j+1)) + &
Kh(I,J) = Kh(I,J) + 0.25*( (MEKE%Ku(i,j) + MEKE%Ku(i+1,j+1)) + &
(MEKE%Ku(i+1,j) + MEKE%Ku(i,j+1)) ) * meke_res_fn
endif

! Older method of bounding for stability
if (CS%anisotropic) &
! *Add* the shear component of anisotropic viscosity
Kh(i,j) = Kh(i,j) + CS%Kh_aniso * CS%n1n2_q(I,J)**2
Kh(I,J) = Kh(I,J) + CS%Kh_aniso * CS%n1n2_q(I,J)**2

! Newer method of bounding for stability
if (CS%better_bound_Kh) then
if (Kh(i,j) >= hrat_min(i,j) * CS%Kh_Max_xy(I,J)) then
if (Kh(i,j) >= hrat_min(I,J) * CS%Kh_Max_xy(I,J)) then
visc_bound_rem(i,j) = 0.0
Kh(i,j) = hrat_min(i,j) * CS%Kh_Max_xy(I,J)
elseif (CS%Kh_Max_xy(I,J)>0.) then
visc_bound_rem(i,j) = 1.0 - Kh(i,j) / (hrat_min(i,j) * CS%Kh_Max_xy(I,J))
Kh(i,j) = hrat_min(I,J) * CS%Kh_Max_xy(I,J)
elseif (hrat_min(I,J)*CS%Kh_Max_xy(I,J)>0.) then
visc_bound_rem(I,J) = 1.0 - Kh(I,J) / (hrat_min(I,J) * CS%Kh_Max_xy(I,J))
endif
endif

if (CS%id_Kh_q>0 .or. CS%debug) &
Kh_q(I,J,k) = Kh(i,j)
Kh_q(I,J,k) = Kh(I,J)

if (CS%id_vort_xy_q>0) &
vort_xy_q(I,J,k) = vort_xy(I,J)
Expand Down Expand Up @@ -1352,37 +1353,37 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
if (CS%Smagorinsky_Ah) then
if (CS%bound_Coriolis) then
do J=js-1,Jeq ; do I=is-1,Ieq
AhSm = Shear_mag(i,j) * (CS%Biharm_const_xy(I,J) &
+ CS%Biharm_const2_xy(I,J) * Shear_mag(i,j) &
AhSm = Shear_mag(I,J) * (CS%Biharm_const_xy(I,J) &
+ CS%Biharm_const2_xy(I,J) * Shear_mag(I,J) &
)
Ah(i,j) = max(Ah(I,J), AhSm)
enddo ; enddo
else
do J=js-1,Jeq ; do I=is-1,Ieq
AhSm = CS%Biharm_const_xy(I,J) * Shear_mag(i,j)
Ah(i,j) = max(Ah(I,J), AhSm)
AhSm = CS%Biharm_const_xy(I,J) * Shear_mag(I,J)
Ah(I,J) = max(Ah(I,J), AhSm)
enddo ; enddo
endif
endif

if (CS%Leith_Ah) then
do J=js-1,Jeq ; do I=is-1,Ieq
AhLth = CS%Biharm6_const_xy(I,J) * abs(Del2vort_q(I,J)) * inv_PI6
Ah(i,j) = max(Ah(I,J), AhLth)
Ah(I,J) = max(Ah(I,J), AhLth)
enddo ; enddo
endif

if (CS%bound_Ah .and. .not.CS%better_bound_Ah) then
do J=js-1,Jeq ; do I=is-1,Ieq
Ah(i,j) = min(Ah(i,j), CS%Ah_Max_xy(I,J))
Ah(I,J) = min(Ah(I,J), CS%Ah_Max_xy(I,J))
enddo ; enddo
endif
endif ! Smagorinsky_Ah or Leith_Ah

if (use_MEKE_Au) then
! *Add* the MEKE contribution
do J=js-1,Jeq ; do I=is-1,Ieq
Ah(i,j) = Ah(i,j) + 0.25 * ( &
Ah(I,J) = Ah(I,J) + 0.25 * ( &
(MEKE%Au(i,j) + MEKE%Au(i+1,j+1)) + (MEKE%Au(i+1,j) + MEKE%Au(i,j+1)) &
)
enddo ; enddo
Expand All @@ -1391,31 +1392,31 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
if (CS%Re_Ah > 0.0) then
do J=js-1,Jeq ; do I=is-1,Ieq
KE = 0.125 * ((u(I,j,k) + u(I,j+1,k))**2 + (v(i,J,k) + v(i+1,J,k))**2)
Ah(i,j) = sqrt(KE) * CS%Re_Ah_const_xy(i,j)
Ah(I,J) = sqrt(KE) * CS%Re_Ah_const_xy(i,j)
enddo ; enddo
endif

if (CS%better_bound_Ah) then
if (CS%better_bound_Kh) then
do J=js-1,Jeq ; do I=is-1,Ieq
Ah(i,j) = min(Ah(i,j), visc_bound_rem(i,j) * hrat_min(i,j) * CS%Ah_Max_xy(I,J))
Ah(I,J) = min(Ah(i,j), visc_bound_rem(i,j) * hrat_min(I,J) * CS%Ah_Max_xy(I,J))
enddo ; enddo
else
do J=js-1,Jeq ; do I=is-1,Ieq
Ah(i,j) = min(Ah(i,j), hrat_min(i,j) * CS%Ah_Max_xy(I,J))
Ah(I,J) = min(Ah(i,j), hrat_min(I,J) * CS%Ah_Max_xy(I,J))
enddo ; enddo
endif
endif

if (CS%id_Ah_q>0 .or. CS%debug) then
do J=js-1,Jeq ; do I=is-1,Ieq
Ah_q(I,J,k) = Ah(i,j)
Ah_q(I,J,k) = Ah(I,J)
enddo ; enddo
endif

! Again, need to initialize str_xy as if its biharmonic
do J=js-1,Jeq ; do I=is-1,Ieq
d_str = Ah(i,j) * (dDel2vdx(I,J) + dDel2udy(I,J))
d_str = Ah(I,J) * (dDel2vdx(I,J) + dDel2udy(I,J))

str_xy(I,J) = str_xy(I,J) + d_str

Expand Down
2 changes: 1 addition & 1 deletion src/parameterizations/vertical/MOM_CVMix_conv.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ logical function CVMix_conv_init(Time, G, GV, US, param_file, diag, CS)

call get_param(param_file, mdl, 'KD_CONV', CS%kd_conv_const, &
"Diffusivity used in convective regime. Corresponding viscosity "//&
"(KV_CONV) will be set to KD_CONV * PRANDTL_TURB.", &
"(KV_CONV) will be set to KD_CONV * PRANDTL_CONV.", &
units='m2/s', default=1.00)

call get_param(param_file, mdl, 'BV_SQR_CONV', CS%bv_sqr_conv, &
Expand Down
2 changes: 1 addition & 1 deletion src/parameterizations/vertical/MOM_set_viscosity.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -911,7 +911,7 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
else ; L(K) = L(K)*pbv%por_layer_widthV(i,J,K); endif

! Determine the drag contributing to the bottom boundary layer
! and the Raleigh drag that acts on each layer.
! and the Rayleigh drag that acts on each layer.
if (L(K) > L(K+1)) then
if (vol_below < bbl_thick) then
BBL_frac = (1.0-vol_below/bbl_thick)**2
Expand Down
2 changes: 1 addition & 1 deletion src/parameterizations/vertical/_V_diffusivity.dox
View file
Original file line number Diff line number Diff line change
Expand Up @@ -278,7 +278,7 @@ The original version concentrates buoyancy work in regions of strong stratificat
The shape of the \cite danabasoglu2012 background mixing has a uniform background value, with a dip
at the equator and a bump at \f$\pm 30^{\circ}\f$ degrees latitude. The form is shown in this figure

\image html background_varying.png "Form of the vertically uniform background mixing in \cite danabasoglu2012. The values are symmetric about the equator."
\image html background_varying.png "Form of the vertically uniform background mixing in Danabasoglu [2012]. The values are symmetric about the equator."
\imagelatex{background_varying.png,Form of the vertically uniform background mixing in \cite danabasoglu2012. The values are symmetric about the equator.,\includegraphics[width=\textwidth\,height=\textheight/2\,keepaspectratio=true]}

Some parameters of this curve are set in the input file, some are hard-coded in calculate_bkgnd_mixing.
Expand Down
Loading

0 comments on commit 4b8645f

Please sign in to comment.