Skip to content

Commit

Permalink
+Rescaled variables in coord_adapt
Browse files Browse the repository at this point in the history 
  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.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Apr 17, 2020
1 parent 9b53927 commit 760dc39
Show file tree
Hide file tree
Showing 2 changed files with 42 additions and 38 deletions.
13 changes: 7 additions & 6 deletions src/ALE/MOM_regridding.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)
Expand Down
67 changes: 35 additions & 32 deletions src/ALE/coord_adapt.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -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

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)) + &
Expand All @@ -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)) + &
Expand All @@ -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)) + &
Expand All @@ -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)) + &
Expand All @@ -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
Expand All @@ -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
Expand Down

0 comments on commit 760dc39

Please sign in to comment.