From 760dc39609a47f7974066b67d1f590f4c7b5d784 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 17 Apr 2020 19:53:11 -0400 Subject: [PATCH] +Rescaled variables in coord_adapt Rescaled internal density and pressure variables in coord_adapt, as well as some input parameters. These changes include rescaling of the input and output variables associated with the calculate_density routines. One variable that was being reused with different units has been split into two, and there are new arguments to build_grid_adaptive, build_adapt_column, and init_coord_adapt. All answers in the MOM6-examples test suite are bitwise identical. --- src/ALE/MOM_regridding.F90 | 13 ++++---- src/ALE/coord_adapt.F90 | 67 ++++++++++++++++++++------------------ 2 files changed, 42 insertions(+), 38 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index f6791a3b73..8be7824193 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -599,7 +599,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m call set_regrid_params(CS, adaptTimeRatio=adaptTimeRatio, adaptZoom=adaptZoom, & adaptZoomCoeff=adaptZoomCoeff, adaptBuoyCoeff=adaptBuoyCoeff, adaptAlpha=adaptAlpha, & - adaptDoMin=tmpLogical, adaptDrho0=US%R_to_kg_m3*adaptDrho0) + adaptDoMin=tmpLogical, adaptDrho0=adaptDrho0) endif if (main_parameters .and. coord_is_state_dependent) then @@ -885,7 +885,7 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_ call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case ( REGRIDDING_ADAPTIVE ) - call build_grid_adaptive(G, GV, h, tv, dzInterface, remapCS, CS) + call build_grid_adaptive(G, GV, G%US, h, tv, dzInterface, remapCS, CS) call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new) case default @@ -1527,9 +1527,10 @@ end subroutine build_grid_HyCOM1 !> This subroutine builds an adaptive grid that follows density surfaces where !! possible, subject to constraints on the smoothness of interface heights. -subroutine build_grid_adaptive(G, GV, h, tv, dzInterface, remapCS, CS) +subroutine build_grid_adaptive(G, GV, US, h, tv, dzInterface, remapCS, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various !! thermodynamic variables @@ -1575,7 +1576,7 @@ subroutine build_grid_adaptive(G, GV, h, tv, dzInterface, remapCS, CS) cycle endif - call build_adapt_column(CS%adapt_CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) + call build_adapt_column(CS%adapt_CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNext) call filtered_grid_motion(CS, nz, zInt(i,j,:), zNext, dzInterface(i,j,:)) ! convert from depth to z @@ -1990,7 +1991,7 @@ subroutine initCoord(CS, GV, US, coord_mode) call init_coord_slight(CS%slight_CS, CS%nk, CS%ref_pressure, CS%target_density, & CS%interp_CS, GV%m_to_H) case (REGRIDDING_ADAPTIVE) - call init_coord_adapt(CS%adapt_CS, CS%nk, CS%coordinateResolution, GV%m_to_H) + call init_coord_adapt(CS%adapt_CS, CS%nk, CS%coordinateResolution, GV%m_to_H, US%kg_m3_to_R) end select end subroutine initCoord @@ -2272,7 +2273,7 @@ subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_gri !! preventing interfaces from being shallower than !! the depths specified by the regridding coordinate. real, optional, intent(in) :: adaptDrho0 !< Reference density difference for stratification-dependent - !! diffusion. [kg m-3] + !! diffusion. [R ~> kg m-3] if (present(interp_scheme)) call set_interp_scheme(CS%interp_CS, interp_scheme) if (present(boundary_extrapolation)) call set_interp_extrap(CS%interp_CS, boundary_extrapolation) diff --git a/src/ALE/coord_adapt.F90 b/src/ALE/coord_adapt.F90 index 3a083af2db..383bf6a055 100644 --- a/src/ALE/coord_adapt.F90 +++ b/src/ALE/coord_adapt.F90 @@ -5,6 +5,7 @@ module coord_adapt use MOM_EOS, only : calculate_density_derivs use MOM_error_handler, only : MOM_error, FATAL +use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : ocean_grid_type, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type @@ -36,7 +37,7 @@ module coord_adapt !> Stratification-dependent diffusion coefficient real :: adaptBuoyCoeff - !> Reference density difference for stratification-dependent diffusion [kg m-3] + !> Reference density difference for stratification-dependent diffusion [R ~> kg m-3] real :: adaptDrho0 !> If true, form a HYCOM1-like mixed layet by preventing interfaces @@ -49,31 +50,28 @@ module coord_adapt contains !> Initialise an adapt_CS with parameters -subroutine init_coord_adapt(CS, nk, coordinateResolution, m_to_H) +subroutine init_coord_adapt(CS, nk, coordinateResolution, m_to_H, kg_m3_to_R) type(adapt_CS), pointer :: CS !< Unassociated pointer to hold the control structure integer, intent(in) :: nk !< Number of layers in the grid real, dimension(:), intent(in) :: coordinateResolution !< Nominal near-surface resolution [m] or !! other units specified with m_to_H - real, optional, intent(in) :: m_to_H !< A conversion factor from m to the units of thicknesses - - real :: m_to_H_rescale ! A unit conversion factor. + real, intent(in) :: m_to_H !< A conversion factor from m to the units of thicknesses + real, intent(in) :: kg_m3_to_R !< A conversion factor from kg m-3 to the units of density if (associated(CS)) call MOM_error(FATAL, "init_coord_adapt: CS already associated") allocate(CS) allocate(CS%coordinateResolution(nk)) - m_to_H_rescale = 1.0 ; if (present(m_to_H)) m_to_H_rescale = m_to_H - CS%nk = nk CS%coordinateResolution(:) = coordinateResolution(:) ! Set real parameter default values CS%adaptTimeRatio = 1e-1 ! Nondim. CS%adaptAlpha = 1.0 ! Nondim. - CS%adaptZoom = 200.0 * m_to_H_rescale + CS%adaptZoom = 200.0 * m_to_H ! [H ~> m or kg m-2] CS%adaptZoomCoeff = 0.0 ! Nondim. CS%adaptBuoyCoeff = 0.0 ! Nondim. - CS%adaptDrho0 = 0.5 ! [kg m-3] + CS%adaptDrho0 = 0.5 * kg_m3_to_R ! [R ~> kg m-3] end subroutine init_coord_adapt @@ -98,7 +96,7 @@ subroutine set_adapt_params(CS, adaptTimeRatio, adaptAlpha, adaptZoom, adaptZoom real, optional, intent(in) :: adaptZoomCoeff !< Near-surface zooming coefficient real, optional, intent(in) :: adaptBuoyCoeff !< Stratification-dependent diffusion coefficient real, optional, intent(in) :: adaptDrho0 !< Reference density difference for - !! stratification-dependent diffusion + !! stratification-dependent diffusion [R ~> kg m-3] logical, optional, intent(in) :: adaptDoMin !< If true, form a HYCOM1-like mixed layer by !! preventing interfaces from becoming shallower than !! the depths set by coordinateResolution @@ -114,10 +112,11 @@ subroutine set_adapt_params(CS, adaptTimeRatio, adaptAlpha, adaptZoom, adaptZoom if (present(adaptDoMin)) CS%adaptDoMin = adaptDoMin end subroutine set_adapt_params -subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) +subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNext) type(adapt_CS), intent(in) :: CS !< The control structure for this module type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various !! thermodynamic variables integer, intent(in) :: i !< The i-index of the column to work on @@ -130,8 +129,12 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) ! Local variables integer :: k, nz - real :: h_up, b1, b_denom_1, d1, depth, drdz, nominal_z, stretching - real, dimension(SZK_(GV)+1) :: alpha, beta, del2sigma ! drho/dT and drho/dS + real :: h_up, b1, b_denom_1, d1, depth, nominal_z, stretching + real :: drdz ! The vertical density gradient [R H-1 ~> kg m-4 or m-1] + real, dimension(SZK_(GV)+1) :: alpha ! drho/dT [R degC-1 ~> kg m-3 degC-1] + real, dimension(SZK_(GV)+1) :: beta ! drho/dS [R ppt-1 ~> kg m-3 ppt-1] + real, dimension(SZK_(GV)+1) :: del2sigma ! Laplacian of in situ density times grid spacing [R ~> kg m-3] + real, dimension(SZK_(GV)+1) :: dh_d2s ! Thickness change in response to del2sigma [H ~> m or kg m-2] real, dimension(SZK_(GV)) :: kGrid, c1 ! grid diffusivity on layers, and tridiagonal work array nz = CS%nk @@ -143,8 +146,8 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) ! local depth for scaling diffusivity depth = G%bathyT(i,j) * GV%Z_to_H - ! initialize del2sigma to zero - del2sigma(:) = 0. + ! initialize del2sigma and the thickness change response to it zero + del2sigma(:) = 0.0 ; dh_d2s(:) = 0.0 ! calculate del-squared of neutral density by a ! stencilled finite difference @@ -155,8 +158,8 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) call calculate_density_derivs( & 0.5 * (tInt(i,j,2:nz) + tInt(i,j-1,2:nz)), & 0.5 * (sInt(i,j,2:nz) + sInt(i,j-1,2:nz)), & - 0.5 * (zInt(i,j,2:nz) + zInt(i,j-1,2:nz)) * GV%H_to_Pa, & - alpha, beta, tv%eqn_of_state, dom=(/2,nz/)) + 0.5 * (zInt(i,j,2:nz) + zInt(i,j-1,2:nz)) * (GV%H_to_RZ * GV%g_Earth), & + alpha, beta, tv%eqn_of_state, US, dom=(/2,nz/)) del2sigma(2:nz) = del2sigma(2:nz) + & (alpha(2:nz) * (tInt(i,j-1,2:nz) - tInt(i,j,2:nz)) + & @@ -167,8 +170,8 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) call calculate_density_derivs( & 0.5 * (tInt(i,j,2:nz) + tInt(i,j+1,2:nz)), & 0.5 * (sInt(i,j,2:nz) + sInt(i,j+1,2:nz)), & - 0.5 * (zInt(i,j,2:nz) + zInt(i,j+1,2:nz)) * GV%H_to_Pa, & - alpha, beta, tv%eqn_of_state, dom=(/2,nz/)) + 0.5 * (zInt(i,j,2:nz) + zInt(i,j+1,2:nz)) * (GV%H_to_RZ * GV%g_Earth), & + alpha, beta, tv%eqn_of_state, US, dom=(/2,nz/)) del2sigma(2:nz) = del2sigma(2:nz) + & (alpha(2:nz) * (tInt(i,j+1,2:nz) - tInt(i,j,2:nz)) + & @@ -179,8 +182,8 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) call calculate_density_derivs( & 0.5 * (tInt(i,j,2:nz) + tInt(i-1,j,2:nz)), & 0.5 * (sInt(i,j,2:nz) + sInt(i-1,j,2:nz)), & - 0.5 * (zInt(i,j,2:nz) + zInt(i-1,j,2:nz)) * GV%H_to_Pa, & - alpha, beta, tv%eqn_of_state, dom=(/2,nz/)) + 0.5 * (zInt(i,j,2:nz) + zInt(i-1,j,2:nz)) * (GV%H_to_RZ * GV%g_Earth), & + alpha, beta, tv%eqn_of_state, US, dom=(/2,nz/)) del2sigma(2:nz) = del2sigma(2:nz) + & (alpha(2:nz) * (tInt(i-1,j,2:nz) - tInt(i,j,2:nz)) + & @@ -191,8 +194,8 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) call calculate_density_derivs( & 0.5 * (tInt(i,j,2:nz) + tInt(i+1,j,2:nz)), & 0.5 * (sInt(i,j,2:nz) + sInt(i+1,j,2:nz)), & - 0.5 * (zInt(i,j,2:nz) + zInt(i+1,j,2:nz)) * GV%H_to_Pa, & - alpha, beta, tv%eqn_of_state, dom=(/2,nz/)) + 0.5 * (zInt(i,j,2:nz) + zInt(i+1,j,2:nz)) * (GV%H_to_RZ * GV%g_Earth), & + alpha, beta, tv%eqn_of_state, US, dom=(/2,nz/)) del2sigma(2:nz) = del2sigma(2:nz) + & (alpha(2:nz) * (tInt(i+1,j,2:nz) - tInt(i,j,2:nz)) + & @@ -205,23 +208,23 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) ! ! a positive curvature means we're too light relative to adjacent columns, ! so del2sigma needs to be positive too (push the interface deeper) - call calculate_density_derivs(tInt(i,j,:), sInt(i,j,:), zInt(i,j,:) * GV%H_to_Pa, & - alpha, beta, tv%eqn_of_state, dom=(/1,nz/)) !### This should be (/1,nz+1/) - see 25 lines below. + call calculate_density_derivs(tInt(i,j,:), sInt(i,j,:), zInt(i,j,:) * (GV%H_to_RZ * GV%g_Earth), & + alpha, beta, tv%eqn_of_state, US, dom=(/1,nz/)) !### This should be (/1,nz+1/) - see 25 lines below. do K = 2, nz ! TODO make lower bound here configurable - del2sigma(K) = del2sigma(K) * (0.5 * (h(i,j,k-1) + h(i,j,k))) / & + dh_d2s(K) = del2sigma(K) * (0.5 * (h(i,j,k-1) + h(i,j,k))) / & max(alpha(K) * (tv%T(i,j,k) - tv%T(i,j,k-1)) + & - beta(K) * (tv%S(i,j,k) - tv%S(i,j,k-1)), 1e-20) + beta(K) * (tv%S(i,j,k) - tv%S(i,j,k-1)), 1e-20*US%kg_m3_to_R) ! don't move the interface so far that it would tangle with another ! interface in the direction we're moving (or exceed a Nyquist limit ! that could cause oscillations of the interface) - h_up = merge(h(i,j,k), h(i,j,k-1), del2sigma(K) > 0.) - del2sigma(K) = 0.5 * CS%adaptAlpha * & - sign(min(abs(del2sigma(K)), 0.5 * h_up), del2sigma(K)) + h_up = merge(h(i,j,k), h(i,j,k-1), dh_d2s(K) > 0.) + dh_d2s(K) = 0.5 * CS%adaptAlpha * & + sign(min(abs(del2sigma(K)), 0.5 * h_up), dh_d2s(K)) ! update interface positions so we can diffuse them - zNext(K) = zInt(i,j,K) + del2sigma(K) + zNext(K) = zInt(i,j,K) + dh_d2s(K) enddo ! solve diffusivity equation to smooth grid @@ -233,7 +236,7 @@ subroutine build_adapt_column(CS, G, GV, tv, i, j, zInt, tInt, sInt, h, zNext) do k = 1, nz ! calculate the dr bit of drdz drdz = 0.5 * (alpha(K) + alpha(K+1)) * (tInt(i,j,K+1) - tInt(i,j,K)) + & - 0.5 * (beta(K) + beta(K+1)) * (sInt(i,j,K+1) - sInt(i,j,K)) + 0.5 * (beta(K) + beta(K+1)) * (sInt(i,j,K+1) - sInt(i,j,K)) ! divide by dz from the new interface positions drdz = drdz / (zNext(K) - zNext(K+1) + GV%H_subroundoff) ! don't do weird stuff in unstably-stratified regions