From 6f1a191738937312c88ac227840175e1b1baf203 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 15 Aug 2021 04:54:02 -0400 Subject: [PATCH 01/11] +Add optional dZref argument to find_eta Added an optional dZref argument to the two find_eta routines so that the bathymetry can use a different reference height than is used for the interface heights. By default, they use the same reference height and all answers are bitwise identical. There is a new optional arguments (that is not yet being exercised with this commit, but has been tested and will be used in an upcoming commit) to two publicly visible interfaces. --- src/core/MOM_interface_heights.F90 | 91 +++++++++++++++++------------- 1 file changed, 53 insertions(+), 38 deletions(-) diff --git a/src/core/MOM_interface_heights.F90 b/src/core/MOM_interface_heights.F90 index ec7501c5f0..17729e586c 100644 --- a/src/core/MOM_interface_heights.F90 +++ b/src/core/MOM_interface_heights.F90 @@ -28,7 +28,7 @@ module MOM_interface_heights !! form for consistency with the calculation of the pressure gradient forces. !! Additionally, these height may be dilated for consistency with the !! corresponding time-average quantity from the barotropic calculation. -subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) +subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m, dZref) 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 @@ -37,14 +37,17 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) !! thermodynamic variables. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: eta !< layer interface heights !! [Z ~> m] or [1/eta_to_m m]. - real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic - !! variable that gives the "correct" free surface height (Boussinesq) or total water - !! column mass per unit area (non-Boussinesq). This is used to dilate the layer. - !! thicknesses when calculating interfaceheights [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic variable + !! that gives the "correct" free surface height (Boussinesq) or total water + !! column mass per unit area (non-Boussinesq). This is used to dilate the layer + !! thicknesses when calculating interface heights [H ~> m or kg m-2]. + !! In Boussinesq mode, eta_bt and G%bathyT use the same reference height. integer, optional, intent(in) :: halo_size !< width of halo points on !! which to calculate eta. real, optional, intent(in) :: eta_to_m !< The conversion factor from - !! the units of eta to m; by default this is US%Z_to_m. + !! the units of eta to m; by default this is US%Z_to_m. + real, optional, intent(in) :: dZref !< The difference in the + !! reference height between G%bathyT and eta [Z ~> m]. The default is 0. ! Local variables real :: p(SZI_(G),SZJ_(G),SZK_(GV)+1) ! Hydrostatic pressure at each interface [R L2 T-2 ~> Pa] @@ -55,6 +58,8 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) real :: I_gEarth ! The inverse of the gravitational acceleration times the ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1] real :: Z_to_eta, H_to_eta, H_to_rho_eta ! Unit conversion factors with obvious names. + real :: dZ_ref ! The difference in the reference height between G%bathyT and eta [Z ~> m]. + ! dZ_ref is 0 unless the optional argument dZref is present. integer i, j, k, isv, iev, jsv, jev, nz, halo halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) @@ -69,33 +74,35 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) H_to_eta = GV%H_to_Z * Z_to_eta H_to_rho_eta = GV%H_to_RZ * Z_to_eta I_gEarth = Z_to_eta / GV%g_Earth + dZ_ref = 0.0 ; if (present(dZref)) dZ_ref = dZref -!$OMP parallel default(shared) private(dilate,htot) -!$OMP do - do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -Z_to_eta*G%bathyT(i,j) ; enddo ; enddo + !$OMP parallel default(shared) private(dilate,htot) + !$OMP do + do j=jsv,jev ; do i=isv,iev ; eta(i,j,nz+1) = -Z_to_eta*(G%bathyT(i,j) + dZ_ref) ; enddo ; enddo if (GV%Boussinesq) then -!$OMP do + !$OMP do do j=jsv,jev ; do k=nz,1,-1 ; do i=isv,iev eta(i,j,K) = eta(i,j,K+1) + h(i,j,k)*H_to_eta enddo ; enddo ; enddo if (present(eta_bt)) then ! Dilate the water column to agree with the free surface height ! that is used for the dynamics. -!$OMP do + !$OMP do do j=jsv,jev do i=isv,iev dilate(i) = (eta_bt(i,j)*H_to_eta + Z_to_eta*G%bathyT(i,j)) / & - (eta(i,j,1) + Z_to_eta*G%bathyT(i,j)) + (eta(i,j,1) + Z_to_eta*(G%bathyT(i,j) + dZ_ref)) enddo do k=1,nz ; do i=isv,iev - eta(i,j,K) = dilate(i) * (eta(i,j,K) + Z_to_eta*G%bathyT(i,j)) - Z_to_eta*G%bathyT(i,j) + eta(i,j,K) = dilate(i) * (eta(i,j,K) + Z_to_eta*(G%bathyT(i,j) + dZ_ref)) - & + Z_to_eta*(G%bathyT(i,j) + dZ_ref) enddo ; enddo enddo endif else if (associated(tv%eqn_of_state)) then -!$OMP do + !$OMP do do j=jsv,jev if (associated(tv%p_surf)) then do i=isv,iev ; p(i,j,1) = tv%p_surf(i,j) ; enddo @@ -106,19 +113,19 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) p(i,j,K+1) = p(i,j,K) + GV%g_Earth*GV%H_to_RZ*h(i,j,k) enddo ; enddo enddo -!$OMP do + !$OMP do do k=1,nz call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,K), p(:,:,K+1), & 0.0, G%HI, tv%eqn_of_state, US, dz_geo(:,:,k), halo_size=halo) enddo -!$OMP do + !$OMP do do j=jsv,jev do k=nz,1,-1 ; do i=isv,iev eta(i,j,K) = eta(i,j,K+1) + I_gEarth * dz_geo(i,j,k) enddo ; enddo enddo else -!$OMP do + !$OMP do do j=jsv,jev ; do k=nz,1,-1 ; do i=isv,iev eta(i,j,K) = eta(i,j,K+1) + H_to_rho_eta*h(i,j,k) / GV%Rlay(k) enddo ; enddo ; enddo @@ -126,18 +133,19 @@ subroutine find_eta_3d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) if (present(eta_bt)) then ! Dilate the water column to agree with the free surface height ! from the time-averaged barotropic solution. -!$OMP do + !$OMP do do j=jsv,jev do i=isv,iev ; htot(i) = GV%H_subroundoff ; enddo do k=1,nz ; do i=isv,iev ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo do i=isv,iev ; dilate(i) = eta_bt(i,j) / htot(i) ; enddo do k=1,nz ; do i=isv,iev - eta(i,j,K) = dilate(i) * (eta(i,j,K) + Z_to_eta*G%bathyT(i,j)) - Z_to_eta*G%bathyT(i,j) + eta(i,j,K) = dilate(i) * (eta(i,j,K) + Z_to_eta*(G%bathyT(i,j) + dZ_ref)) - & + Z_to_eta*(G%bathyT(i,j) + dZ_ref) enddo ; enddo enddo endif endif -!$OMP end parallel + !$OMP end parallel end subroutine find_eta_3d @@ -145,7 +153,7 @@ end subroutine find_eta_3d !! with the calculation of the pressure gradient forces. Additionally, the sea !! surface height may be adjusted for consistency with the corresponding !! time-average quantity from the barotropic calculation. -subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) +subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m, dZref) 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 @@ -155,12 +163,16 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta !< free surface height relative to !! mean sea level (z=0) often [Z ~> m]. real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: eta_bt !< optional barotropic - !! variable that gives the "correct" free surface height (Boussinesq) or total - !! water column mass per unit area (non-Boussinesq) [H ~> m or kg m-2]. + !! variable that gives the "correct" free surface height (Boussinesq) or total + !! water column mass per unit area (non-Boussinesq) [H ~> m or kg m-2]. + !! In Boussinesq mode, eta_bt and G%bathyT use the same reference height. integer, optional, intent(in) :: halo_size !< width of halo points on !! which to calculate eta. real, optional, intent(in) :: eta_to_m !< The conversion factor from - !! the units of eta to m; by default this is US%Z_to_m. + !! the units of eta to m; by default this is US%Z_to_m. + real, optional, intent(in) :: dZref !< The difference in the + !! reference height between G%bathyT and eta [Z ~> m]. The default is 0. + ! Local variables real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & p ! Hydrostatic pressure at each interface [R L2 T-2 ~> Pa] @@ -170,6 +182,8 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) real :: I_gEarth ! The inverse of the gravitational acceleration times the ! rescaling factor derived from eta_to_m [T2 Z L-2 ~> s2 m-1] real :: Z_to_eta, H_to_eta, H_to_rho_eta ! Unit conversion factors with obvious names. + real :: dZ_ref ! The difference in the reference height between G%bathyT and eta [Z ~> m]. + ! dZ_ref is 0 unless the optional argument dZref is present. integer i, j, k, is, ie, js, je, nz, halo halo = 0 ; if (present(halo_size)) halo = max(0,halo_size) @@ -180,26 +194,27 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) H_to_eta = GV%H_to_Z * Z_to_eta H_to_rho_eta = GV%H_to_RZ * Z_to_eta I_gEarth = Z_to_eta / GV%g_Earth + dZ_ref = 0.0 ; if (present(dZref)) dZ_ref = dZref -!$OMP parallel default(shared) private(htot) -!$OMP do - do j=js,je ; do i=is,ie ; eta(i,j) = -Z_to_eta*G%bathyT(i,j) ; enddo ; enddo + !$OMP parallel default(shared) private(htot) + !$OMP do + do j=js,je ; do i=is,ie ; eta(i,j) = -Z_to_eta*(G%bathyT(i,j) + dZ_ref) ; enddo ; enddo if (GV%Boussinesq) then if (present(eta_bt)) then -!$OMP do + !$OMP do do j=js,je ; do i=is,ie - eta(i,j) = H_to_eta*eta_bt(i,j) + eta(i,j) = H_to_eta*eta_bt(i,j) - Z_to_eta*dZ_ref enddo ; enddo else -!$OMP do + !$OMP do do j=js,je ; do k=1,nz ; do i=is,ie eta(i,j) = eta(i,j) + h(i,j,k)*H_to_eta enddo ; enddo ; enddo endif else if (associated(tv%eqn_of_state)) then -!$OMP do + !$OMP do do j=js,je if (associated(tv%p_surf)) then do i=is,ie ; p(i,j,1) = tv%p_surf(i,j) ; enddo @@ -211,17 +226,17 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) p(i,j,k+1) = p(i,j,k) + GV%g_Earth*GV%H_to_RZ*h(i,j,k) enddo ; enddo enddo -!$OMP do + !$OMP do do k = 1, nz call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p(:,:,k), p(:,:,k+1), 0.0, & G%HI, tv%eqn_of_state, US, dz_geo(:,:,k), halo_size=halo) enddo -!$OMP do + !$OMP do do j=js,je ; do k=1,nz ; do i=is,ie eta(i,j) = eta(i,j) + I_gEarth * dz_geo(i,j,k) enddo ; enddo ; enddo else -!$OMP do + !$OMP do do j=js,je ; do k=1,nz ; do i=is,ie eta(i,j) = eta(i,j) + H_to_rho_eta*h(i,j,k) / GV%Rlay(k) enddo ; enddo ; enddo @@ -229,18 +244,18 @@ subroutine find_eta_2d(h, tv, G, GV, US, eta, eta_bt, halo_size, eta_to_m) if (present(eta_bt)) then ! Dilate the water column to agree with the time-averaged column ! mass from the barotropic solution. -!$OMP do + !$OMP do do j=js,je do i=is,ie ; htot(i) = GV%H_subroundoff ; enddo do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo do i=is,ie - eta(i,j) = (eta_bt(i,j) / htot(i)) * (eta(i,j) + Z_to_eta*G%bathyT(i,j)) - & - Z_to_eta*G%bathyT(i,j) + eta(i,j) = (eta_bt(i,j) / htot(i)) * (eta(i,j) + Z_to_eta*(G%bathyT(i,j) + dZ_ref)) - & + Z_to_eta*(G%bathyT(i,j) + dZ_ref) enddo enddo endif endif -!$OMP end parallel + !$OMP end parallel end subroutine find_eta_2d From a2df3b75e92d80b5751408d12cb5f1f3dd6f9f1c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 15 Aug 2021 04:58:27 -0400 Subject: [PATCH 02/11] Code cleanup in MOM_regridding Eliminated the depth argument to check_grid_column, along with some internal calculations that are never used. Also corrected some comments. All answers are bitwise identical. --- src/ALE/MOM_regridding.F90 | 38 +++++++++++++++++--------------------- 1 file changed, 17 insertions(+), 21 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 77a2b4ea8b..5b19a7549c 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -880,21 +880,20 @@ subroutine check_remapping_grid( G, GV, h, dzInterface, msg ) !$OMP parallel do default(shared) do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 - if (G%mask2dT(i,j)>0.) call check_grid_column( GV%ke, GV%Z_to_H*G%bathyT(i,j), h(i,j,:), dzInterface(i,j,:), msg ) + if (G%mask2dT(i,j)>0.) call check_grid_column( GV%ke, h(i,j,:), dzInterface(i,j,:), msg ) enddo ; enddo end subroutine check_remapping_grid !> Check that the total thickness of new and old grids are consistent -subroutine check_grid_column( nk, depth, h, dzInterface, msg ) +subroutine check_grid_column( nk, h, dzInterface, msg ) integer, intent(in) :: nk !< Number of cells - real, intent(in) :: depth !< Depth of bottom [Z ~> m] or arbitrary units real, dimension(nk), intent(in) :: h !< Cell thicknesses [Z ~> m] or arbitrary units real, dimension(nk+1), intent(in) :: dzInterface !< Change in interface positions (same units as h) character(len=*), intent(in) :: msg !< Message to append to errors ! Local variables integer :: k - real :: eps, total_h_old, total_h_new, h_new, z_old, z_new + real :: eps, total_h_old, total_h_new, h_new eps =1. ; eps = epsilon(eps) @@ -904,13 +903,8 @@ subroutine check_grid_column( nk, depth, h, dzInterface, msg ) total_h_old = total_h_old + h(k) enddo - ! Integrate upwards for the interfaces consistent with the rest of MOM6 - z_old = - depth - if (depth == 0.) z_old = - total_h_old total_h_new = 0. do k = nk,1,-1 - z_old = z_old + h(k) ! Old interface position above layer k - z_new = z_old + dzInterface(k) ! New interface position based on dzInterface h_new = h(k) + ( dzInterface(k) - dzInterface(k+1) ) ! New thickness if (h_new<0.) then write(0,*) 'k,h,hnew=',k,h(k),h_new @@ -1082,7 +1076,7 @@ subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g ) end subroutine filtered_grid_motion -!> Builds a z*-ccordinate grid with partial steps (Adcroft and Campin, 2004). +!> Builds a z*-coordinate grid with partial steps (Adcroft and Campin, 2004). !! z* is defined as !! z* = (z-eta)/(H+eta)*H s.t. z*=0 when z=eta and z*=-H when z=-H . subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) @@ -1118,7 +1112,7 @@ subroutine build_zstar_grid( CS, G, GV, h, dzInterface, frac_shelf_h) cycle endif - ! Local depth (G%bathyT is positive) + ! Local depth (G%bathyT is positive downward) nominalDepth = G%bathyT(i,j)*GV%Z_to_H ! Determine water column thickness @@ -1319,7 +1313,7 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS, frac_shel endif - ! Local depth (G%bathyT is positive) + ! Local depth (G%bathyT is positive downward) nominalDepth = G%bathyT(i,j)*GV%Z_to_H ! Determine total water column thickness @@ -1406,7 +1400,7 @@ end subroutine build_rho_grid !! density interpolated from the column profile and a clipping of depth for !! each interface to a fixed z* or p* grid. This should probably be (optionally?) !! changed to find the nearest location of the target density. -!! \remark { Based on Bleck, 2002: An oceanice general circulation model framed in +!! \remark { Based on Bleck, 2002: An ocean-ice general circulation model framed in !! hybrid isopycnic-Cartesian coordinates, Ocean Modelling 37, 55-88. !! http://dx.doi.org/10.1016/S1463-5003(01)00012-9 } subroutine build_grid_HyCOM1( G, GV, US, h, tv, h_new, dzInterface, CS, frac_shelf_h ) @@ -1575,7 +1569,9 @@ subroutine build_grid_SLight(G, GV, US, h, tv, dzInterface, CS) real, dimension(SZK_(GV)+1) :: z_col_new ! Interface positions relative to the surface [H ~> m or kg m-2] real, dimension(SZK_(GV)+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2] real, dimension(SZK_(GV)) :: p_col ! Layer center pressure [R L2 T-2 ~> Pa] - real :: depth + + ! Local variables + real :: depth ! Depth of the ocean relative to the mean sea surface height in thickness units [H ~> m or kg m-2] integer :: i, j, k, nz real :: h_neglect, h_neglect_edge @@ -1631,8 +1627,8 @@ end subroutine build_grid_SLight subroutine adjust_interface_motion( CS, nk, h_old, dz_int ) type(regridding_CS), intent(in) :: CS !< Regridding control structure integer, intent(in) :: nk !< Number of layers in h_old - real, dimension(nk), intent(in) :: h_old !< Minium allowed thickness of h [H ~> m or kg m-2] - real, dimension(CS%nk+1), intent(inout) :: dz_int !< Minium allowed thickness of h [H ~> m or kg m-2] + real, dimension(nk), intent(in) :: h_old !< Minimum allowed thickness of h [H ~> m or kg m-2] + real, dimension(CS%nk+1), intent(inout) :: dz_int !< Minimum allowed thickness of h [H ~> m or kg m-2] ! Local variables integer :: k real :: h_new, eps, h_total, h_err @@ -1710,8 +1706,8 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) real :: total_height real :: delta_h real :: max_depth - real :: eta ! local elevation - real :: local_depth + real :: eta ! local elevation [H ~> m or kg m-2] + real :: local_depth ! The local ocean depth relative to mean sea level in thickness units [H ~> m or kg m-2] real :: x1, y1, x2, y2 real :: x, t @@ -1769,7 +1765,7 @@ subroutine build_grid_arbitrary( G, GV, h, dzInterface, h_new, CS ) endif enddo - ! Chnage in interface position + ! Change in interface position x = 0. ! Left boundary at x=0 dzInterface(i,j,1) = 0. do k = 2,nz @@ -1797,7 +1793,7 @@ subroutine inflate_vanished_layers_old( CS, G, GV, h ) ! objective is to make sure all layers are at least as thick as the minimum ! thickness allowed for regridding purposes (this parameter is set in the ! MOM_input file or defaulted to 1.0e-3). When layers are too thin, they -! are inflated up to the minmum thickness. +! are inflated up to the minimum thickness. !------------------------------------------------------------------------------ ! Arguments @@ -1901,7 +1897,7 @@ function uniformResolution(nk,coordMode,maxDepth,rhoLight,rhoHeavy) ! Arguments integer, intent(in) :: nk !< Number of cells in source grid character(len=*), intent(in) :: coordMode !< A string indicating the coordinate mode. - !! See the documenttion for regrid_consts + !! See the documentation for regrid_consts !! for the recognized values. real, intent(in) :: maxDepth !< The range of the grid values in some modes real, intent(in) :: rhoLight !< The minimum value of the grid in RHO mode From db43b353fe507e277b7fa1ad98e2f5fd753611a2 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 15 Aug 2021 04:58:56 -0400 Subject: [PATCH 03/11] Longer error string in get_polynomial_coordinate Lengthened a message string in get_polynomial_coordinate so that it will give a valid fatal error message in some cases of failures rather than resulting a segmentation fault with no message output. All answers are bitwise identical. --- src/ALE/regrid_interp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ALE/regrid_interp.F90 b/src/ALE/regrid_interp.F90 index 0c758fadaf..87019d46cf 100644 --- a/src/ALE/regrid_interp.F90 +++ b/src/ALE/regrid_interp.F90 @@ -373,7 +373,7 @@ function get_polynomial_coordinate( N, h, x_g, edge_values, ppoly_coefs, & real :: grad ! gradient during N-R iterations [A] integer :: i, k, iter ! loop indices integer :: k_found ! index of target cell - character(len=200) :: mesg + character(len=320) :: mesg logical :: use_2018_answers ! If true use older, less acccurate expressions. eps = NR_OFFSET From 4d8d7afc7301188a97b542f0481e2cc683c3d4ac Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 15 Aug 2021 05:01:06 -0400 Subject: [PATCH 04/11] +Add optional Z_0p argument to int_density_dz Added an optional Z_0p argument to the various Boussinesq density integral routines as a new intercept for the linear expression between the (arbitrary) interface heights and pressure used for the equation of state routines. If omitted, this is equivalent to setting this intercept to 0, and all answers are bitwise identical. There are new optional arguments to several publicly visible interfaces (that are not yet being exercised with this commit, but have been tested and will be used in an upcoming commit). --- src/core/MOM_density_integrals.F90 | 43 +++++++++++++++--------- src/equation_of_state/MOM_EOS.F90 | 8 +++-- src/equation_of_state/MOM_EOS_Wright.F90 | 20 ++++++----- 3 files changed, 43 insertions(+), 28 deletions(-) diff --git a/src/core/MOM_density_integrals.F90 b/src/core/MOM_density_integrals.F90 index 302ba0a714..04e151d5a7 100644 --- a/src/core/MOM_density_integrals.F90 +++ b/src/core/MOM_density_integrals.F90 @@ -38,7 +38,7 @@ module MOM_density_integrals !! required for calculating the finite-volume form pressure accelerations in a !! Boussinesq model. subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p) type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the arrays real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T !< Potential temperature referenced to the surface [degC] @@ -78,13 +78,14 @@ subroutine int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, real, optional, intent(in) :: dz_neglect !< A minuscule thickness change [Z ~> m] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. + real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] if (EOS_quadrature(EOS)) then call int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, US, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p=Z_0p) else call analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p=Z_0p) endif end subroutine int_density_dz @@ -93,8 +94,8 @@ end subroutine int_density_dz !> Calculates (by numerical quadrature) integrals of pressure anomalies across layers, which !! are required for calculating the finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & - EOS, US, dpa, intz_dpa, intx_dpa, inty_dpa, & - bathyT, dz_neglect, useMassWghtInterp, use_inaccurate_form) + EOS, US, dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & + dz_neglect, useMassWghtInterp, use_inaccurate_form, Z_0p) type(hor_index_type), intent(in) :: HI !< Horizontal index type for input variables. real, dimension(SZI_(HI),SZJ_(HI)), & intent(in) :: T !< Potential temperature of the layer [degC] @@ -136,6 +137,8 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & !! interpolate T/S for top and bottom integrals. logical, optional, intent(in) :: use_inaccurate_form !< If true, uses an inaccurate form of !! density anomalies, as was used prior to March 2018. + real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] + ! Local variables real :: T5(5), S5(5) ! Temperatures and salinities at five quadrature points [degC] and [ppt] real :: p5(5) ! Pressures at five quadrature points, never rescaled from Pa [Pa] @@ -148,6 +151,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real :: rho_scale ! A scaling factor for densities from kg m-3 to R [R m3 kg-1 ~> 1] real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3] real :: dz ! The layer thickness [Z ~> m] + real :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: hWght ! A pressure-thickness below topography [Z ~> m] real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Z ~> m] real :: iDenom ! The inverse of the denominator in the weights [Z-2 ~> m-2] @@ -173,6 +177,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & GxRho = US%RL2_T2_to_Pa * G_e * rho_0 rho_ref_mks = rho_ref * US%R_to_kg_m3 I_Rho = 1.0 / rho_0 + z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p use_rho_ref = .true. if (present(use_inaccurate_form)) then if (use_inaccurate_form) use_rho_ref = .not. use_inaccurate_form @@ -191,7 +196,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dz = z_t(i,j) - z_b(i,j) do n=1,5 T5(n) = T(i,j) ; S5(n) = S(i,j) - p5(n) = -GxRho*(z_t(i,j) - 0.25*real(n-1)*dz) + p5(n) = -GxRho*((z_t(i,j) - z0pres) - 0.25*real(n-1)*dz) enddo if (use_rho_ref) then if (rho_scale /= 1.0) then @@ -245,7 +250,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i+1,j) - z_b(i+1,j)) T5(1) = wtT_L*T(i,j) + wtT_R*T(i+1,j) S5(1) = wtT_L*S(i,j) + wtT_R*S(i+1,j) - p5(1) = -GxRho*(wt_L*z_t(i,j) + wt_R*z_t(i+1,j)) + p5(1) = -GxRho*((wt_L*z_t(i,j) + wt_R*z_t(i+1,j)) - z0pres) do n=2,5 T5(n) = T5(1) ; S5(n) = S5(1) ; p5(n) = p5(n-1) + GxRho*0.25*dz enddo @@ -300,7 +305,7 @@ subroutine int_density_dz_generic_pcm(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i,j+1) - z_b(i,j+1)) T5(1) = wtT_L*T(i,j) + wtT_R*T(i,j+1) S5(1) = wtT_L*S(i,j) + wtT_R*S(i,j+1) - p5(1) = -GxRho*(wt_L*z_t(i,j) + wt_R*z_t(i,j+1)) + p5(1) = -GxRho*((wt_L*z_t(i,j) + wt_R*z_t(i,j+1)) - z0pres) do n=2,5 T5(n) = T5(1) ; S5(n) = S5(1) p5(n) = p5(n-1) + GxRho*0.25*dz @@ -336,7 +341,7 @@ end subroutine int_density_dz_generic_pcm subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, dpa, & intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp, & - use_inaccurate_form) + use_inaccurate_form, Z_0p) integer, intent(in) :: k !< Layer index to calculate integrals for type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the input arrays type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure @@ -379,6 +384,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & !! interpolate T/S for top and bottom integrals. logical, optional, intent(in) :: use_inaccurate_form !< If true, uses an inaccurate form of !! density anomalies, as was used prior to March 2018. + real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] ! This subroutine calculates (by numerical quadrature) integrals of ! pressure anomalies across layers, which are required for calculating the @@ -427,6 +433,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & real :: massWeightToggle ! A non-dimensional toggle factor (0 or 1) [nondim] real :: Ttl, Tbl, Ttr, Tbr ! Temperatures at the velocity cell corners [degC] real :: Stl, Sbl, Str, Sbr ! Salinities at the velocity cell corners [ppt] + real :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: hWght ! A topographically limited thicknes weight [Z ~> m] real :: hL, hR ! Thicknesses to the left and right [Z ~> m] real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2] @@ -443,6 +450,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & GxRho = US%RL2_T2_to_Pa * G_e * rho_0 rho_ref_mks = rho_ref * US%R_to_kg_m3 I_Rho = 1.0 / rho_0 + z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p massWeightToggle = 0. if (present(useMassWghtInterp)) then if (useMassWghtInterp) massWeightToggle = 1. @@ -473,7 +481,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & do i = Isq,Ieq+1 dz(i) = e(i,j,K) - e(i,j,K+1) do n=1,5 - p5(i*5+n) = -GxRho*(e(i,j,K) - 0.25*real(n-1)*dz(i)) + p5(i*5+n) = -GxRho*((e(i,j,K) - z0pres) - 0.25*real(n-1)*dz(i)) ! Salinity and temperature points are linearly interpolated S5(i*5+n) = wt_t(n) * S_t(i,j,k) + wt_b(n) * S_b(i,j,k) T5(i*5+n) = wt_t(n) * T_t(i,j,k) + wt_b(n) * T_b(i,j,k) @@ -581,7 +589,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & S15(pos+1) = w_left*Stl + w_right*Str S15(pos+5) = w_left*Sbl + w_right*Sbr - p15(pos+1) = -GxRho*(w_left*e(i,j,K) + w_right*e(i+1,j,K)) + p15(pos+1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i+1,j,K)) - z0pres) ! Pressure do n=2,5 @@ -692,7 +700,7 @@ subroutine int_density_dz_generic_plm(k, tv, T_t, T_b, S_t, S_b, e, rho_ref, & S15(pos+1) = w_left*Stl + w_right*Str S15(pos+5) = w_left*Sbl + w_right*Sbr - p15(pos+1) = -GxRho*(w_left*e(i,j,K) + w_right*e(i,j+1,K)) + p15(pos+1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i,j+1,K)) - z0pres) ! Pressure do n=2,5 @@ -775,7 +783,7 @@ end subroutine int_density_dz_generic_plm !! are parabolic profiles subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & rho_ref, rho_0, G_e, dz_subroundoff, bathyT, HI, GV, EOS, US, & - dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp) + dpa, intz_dpa, intx_dpa, inty_dpa, useMassWghtInterp, Z_0p) integer, intent(in) :: k !< Layer index to calculate integrals for type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structures for the input arrays type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure @@ -816,6 +824,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & !! divided by the y grid spacing [R L2 T-2 ~> Pa] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. + real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] ! This subroutine calculates (by numerical quadrature) integrals of ! pressure anomalies across layers, which are required for calculating the @@ -854,6 +863,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & real :: t6 ! PPM curvature coefficient for T [degC] real :: T_top, T_mn, T_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of T real :: S_top, S_mn, S_bot ! Left edge, cell mean and right edge values used in PPM reconstructions of S + real :: z0pres ! The height at which the pressure is zero [Z ~> m] real :: hWght ! A topographically limited thicknes weight [Z ~> m] real :: hL, hR ! Thicknesses to the left and right [Z ~> m] real :: iDenom ! The denominator of the thickness weight expressions [Z-2 ~> m-2] @@ -868,6 +878,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & GxRho = US%RL2_T2_to_Pa * G_e * rho_0 rho_ref_mks = rho_ref * US%R_to_kg_m3 I_Rho = 1.0 / rho_0 + z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p massWeightToggle = 0. if (present(useMassWghtInterp)) then if (useMassWghtInterp) massWeightToggle = 1. @@ -900,7 +911,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & endif dz = e(i,j,K) - e(i,j,K+1) do n=1,5 - p5(n) = -GxRho*(e(i,j,K) - 0.25*real(n-1)*dz) + p5(n) = -GxRho*((e(i,j,K) - z0pres) - 0.25*real(n-1)*dz) ! Salinity and temperature points are reconstructed with PPM S5(n) = wt_t(n) * S_t(i,j,k) + wt_b(n) * ( S_b(i,j,k) + s6 * wt_t(n) ) T5(n) = wt_t(n) * T_t(i,j,k) + wt_b(n) * ( T_b(i,j,k) + t6 * wt_t(n) ) @@ -978,7 +989,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & ! Pressure dz = w_left*(e(i,j,K) - e(i,j,K+1)) + w_right*(e(i+1,j,K) - e(i+1,j,K+1)) - p5(1) = -GxRho*(w_left*e(i,j,K) + w_right*e(i+1,j,K)) + p5(1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i+1,j,K)) - z0pres) do n=2,5 p5(n) = p5(n-1) + GxRho*0.25*dz enddo @@ -1066,7 +1077,7 @@ subroutine int_density_dz_generic_ppm(k, tv, T_t, T_b, S_t, S_b, e, & ! Pressure dz = w_left*(e(i,j,K) - e(i,j,K+1)) + w_right*(e(i,j+1,K) - e(i,j+1,K+1)) - p5(1) = -GxRho*(w_left*e(i,j,K) + w_right*e(i,j+1,K)) + p5(1) = -GxRho*((w_left*e(i,j,K) + w_right*e(i,j+1,K)) - z0pres) do n=2,5 p5(n) = p5(n-1) + GxRho*0.25*dz enddo diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 6e74c3ffa3..23f22d8a24 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -1253,7 +1253,7 @@ end subroutine analytic_int_specific_vol_dp !! pressure anomalies across layers, which are required for calculating the !! finite-volume form pressure accelerations in a Boussinesq model. subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, dpa, & - intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp) + intz_dpa, intx_dpa, inty_dpa, bathyT, dz_neglect, useMassWghtInterp, Z_0p) type(hor_index_type), intent(in) :: HI !< Ocean horizontal index structure real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature referenced to the surface [degC] @@ -1292,6 +1292,8 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, real, optional, intent(in) :: dz_neglect !< A miniscule thickness change [Z ~> m] logical, optional, intent(in) :: useMassWghtInterp !< If true, uses mass weighting to !! interpolate T/S for top and bottom integrals. + real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] + ! Local variables real :: rho_scale ! A multiplicative factor by which to scale density from kg m-3 to the ! desired units [R m3 kg-1 ~> 1] @@ -1322,11 +1324,11 @@ subroutine analytic_int_density_dz(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, EOS, if ((rho_scale /= 1.0) .or. (pres_scale /= 1.0)) then call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & - dz_neglect, useMassWghtInterp, rho_scale, pres_scale) + dz_neglect, useMassWghtInterp, rho_scale, pres_scale, Z_0p=Z_0p) else call int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, bathyT, & - dz_neglect, useMassWghtInterp) + dz_neglect, useMassWghtInterp, Z_0p=Z_0p) endif case default call MOM_error(FATAL, "No analytic integration option is available with this EOS!") diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index e5cc9555b7..730687fbf6 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -408,7 +408,7 @@ end subroutine calculate_compress_wright !! finite-volume form pressure accelerations in a Boussinesq model. subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & dpa, intz_dpa, intx_dpa, inty_dpa, & - bathyT, dz_neglect, useMassWghtInterp, rho_scale, pres_scale) + bathyT, dz_neglect, useMassWghtInterp, rho_scale, pres_scale, Z_0p) type(hor_index_type), intent(in) :: HI !< The horizontal index type for the arrays. real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed), & intent(in) :: T !< Potential temperature relative to the surface @@ -451,6 +451,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & !! from kg m-3 to the desired units [R m3 kg-1 ~> 1] real, optional, intent(in) :: pres_scale !< A multiplicative factor to convert pressure !! into Pa [Pa T2 R-1 L-2 ~> 1]. + real, optional, intent(in) :: Z_0p !< The height at which the pressure is 0 [Z ~> m] ! Local variables real, dimension(HI%isd:HI%ied,HI%jsd:HI%jed) :: al0_2d, p0_2d, lambda_2d @@ -461,7 +462,8 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real :: g_Earth ! The gravitational acceleration [m2 Z-1 s-2 ~> m s-2] real :: I_Rho ! The inverse of the Boussinesq density [m3 kg-1] real :: rho_ref_mks ! The reference density in MKS units, never rescaled from kg m-3 [kg m-3] - real :: p_ave, I_al0, I_Lzz + real :: p_ave ! The layer averaged pressure [Pa] + real :: I_al0, I_Lzz real :: dz ! The layer thickness [Z ~> m]. real :: hWght ! A pressure-thickness below topography [Z ~> m]. real :: hL, hR ! Pressure-thicknesses of the columns to the left and right [Z ~> m]. @@ -470,10 +472,11 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & real :: hWt_RL, hWt_RR ! hWt_RA is the weighted influence of A on the right column [nondim]. real :: wt_L, wt_R ! The linear weights of the left and right columns [nondim]. real :: wtT_L, wtT_R ! The weights for tracers from the left and right columns [nondim]. - real :: intz(5) ! The integrals of density with height at the - ! 5 sub-column locations [R L2 T-2 ~> Pa]. + real :: intz(5) ! The gravitational acceleration times the integrals of density + ! with height at the 5 sub-column locations [R L2 T-2 ~> Pa] or [Pa]. real :: Pa_to_RL2_T2 ! A conversion factor of pressures from Pa to the output units indicated by ! pres_scale [R L2 T-2 Pa-1 ~> 1] or [1]. + real :: z0pres ! The height at which the pressure is zero [Z ~> m] logical :: do_massWeight ! Indicates whether to do mass weighting. real, parameter :: C1_3 = 1.0/3.0, C1_7 = 1.0/7.0 ! Rational constants. real, parameter :: C1_9 = 1.0/9.0, C1_90 = 1.0/90.0 ! Rational constants. @@ -499,6 +502,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & else rho_ref_mks = rho_ref ; I_Rho = 1.0 / rho_0 endif + z0pres = 0.0 ; if (present(Z_0p)) z0pres = Z_0p do_massWeight = .false. if (present(useMassWghtInterp)) then ; if (useMassWghtInterp) then @@ -517,7 +521,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & al0 = al0_2d(i,j) ; p0 = p0_2d(i,j) ; lambda = lambda_2d(i,j) dz = z_t(i,j) - z_b(i,j) - p_ave = -0.5*GxRho*(z_t(i,j)+z_b(i,j)) + p_ave = -GxRho*(0.5*(z_t(i,j)+z_b(i,j)) - z0pres) I_al0 = 1.0 / al0 I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave) @@ -561,8 +565,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i+1,j) dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i+1,j) - z_b(i+1,j)) - p_ave = -0.5*GxRho*(wt_L*(z_t(i,j)+z_b(i,j)) + & - wt_R*(z_t(i+1,j)+z_b(i+1,j))) + p_ave = -GxRho*(0.5*(wt_L*(z_t(i,j)+z_b(i,j)) + wt_R*(z_t(i+1,j)+z_b(i+1,j))) - z0pres) I_al0 = 1.0 / al0 I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave) @@ -603,8 +606,7 @@ subroutine int_density_dz_wright(T, S, z_t, z_b, rho_ref, rho_0, G_e, HI, & lambda = wtT_L*lambda_2d(i,j) + wtT_R*lambda_2d(i,j+1) dz = wt_L*(z_t(i,j) - z_b(i,j)) + wt_R*(z_t(i,j+1) - z_b(i,j+1)) - p_ave = -0.5*GxRho*(wt_L*(z_t(i,j)+z_b(i,j)) + & - wt_R*(z_t(i,j+1)+z_b(i,j+1))) + p_ave = -GxRho*(0.5*(wt_L*(z_t(i,j)+z_b(i,j)) + wt_R*(z_t(i,j+1)+z_b(i,j+1))) - z0pres) I_al0 = 1.0 / al0 I_Lzz = 1.0 / (p0 + (lambda * I_al0) + p_ave) From e4ffa765454139d6fd4e9c1a5bed7097a286faa9 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 15 Aug 2021 05:01:53 -0400 Subject: [PATCH 05/11] Corrects the routine indicated in 2 error messages Corrected error messages to indicate the right subroutine in rescale_dyn_horgrid_bathymetry. All answers are bitwise identical. --- src/framework/MOM_dyn_horgrid.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90 index 2a9a381caa..89a59374a7 100644 --- a/src/framework/MOM_dyn_horgrid.F90 +++ b/src/framework/MOM_dyn_horgrid.F90 @@ -294,9 +294,9 @@ subroutine rescale_dyn_horgrid_bathymetry(G, m_in_new_units) if (m_in_new_units == 1.0) return if (m_in_new_units < 0.0) & - call MOM_error(FATAL, "rescale_grid_bathymetry: Negative depth units are not permitted.") + call MOM_error(FATAL, "rescale_dyn_horgrid_bathymetry: Negative depth units are not permitted.") if (m_in_new_units == 0.0) & - call MOM_error(FATAL, "rescale_grid_bathymetry: Zero depth units are not permitted.") + call MOM_error(FATAL, "rescale_dyn_horgrid_bathymetry: Zero depth units are not permitted.") rescale = 1.0 / m_in_new_units do j=jsd,jed ; do i=isd,ied From 7cc0582919c2d3f544e42dfcdf983cb0b334da3e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 15 Aug 2021 05:03:33 -0400 Subject: [PATCH 06/11] Prep initialization for reference height changes Added the optional argument dZ_ref_eta to adjustEtaToFitBathymetry so that the bathymetry can use a different reference height than is used in this routine. Also changed the name and reversed the sign of the depth (now Z_bot) argument to find_interfaces, and a new Z_bottom array in MOM_temp_salt_initialize_from_Z. An extra unit conversion factor was eliminated from MOM_state_init_tests. All answers are bitwise identical, and all of these changes are internal to this module. --- .../MOM_state_initialization.F90 | 72 +++++++++++-------- 1 file changed, 42 insertions(+), 30 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 0f8c772348..62b06a0b57 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -11,7 +11,7 @@ module MOM_state_initialization use MOM_cpu_clock, only : CLOCK_ROUTINE, CLOCK_LOOP use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast use MOM_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE, AGRID -use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, NOTE, WARNING, is_root_pe +use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, read_param, log_param, param_file_type use MOM_file_parser, only : log_version @@ -655,7 +655,7 @@ subroutine initialize_thickness_from_file(h, G, GV, US, param_file, file_has_thi !! only read parameters without changing h. ! Local variables - real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! Interface heights, in depth units. + real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! Interface heights, in depth units [Z ~> m]. integer :: inconsistent = 0 logical :: correct_thickness logical :: just_read ! If true, just read parameters but set nothing. @@ -696,7 +696,7 @@ subroutine initialize_thickness_from_file(h, G, GV, US, param_file, file_has_thi call MOM_read_data(filename, "eta", eta(:,:,:), G%Domain, scale=US%m_to_Z) if (correct_thickness) then - call adjustEtaToFitBathymetry(G, GV, US, eta, h) + call adjustEtaToFitBathymetry(G, GV, US, eta, h, dZ_ref_eta=0.0) else do k=nz,1,-1 ; do j=js,je ; do i=is,ie if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) then @@ -732,25 +732,30 @@ end subroutine initialize_thickness_from_file !! is dilated (expanded) to fill the void. !! @remark{There is a (hard-wired) "tolerance" parameter such that the !! criteria for adjustment must equal or exceed 10cm.} -subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h) +subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h, dZ_ref_eta) 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)+1), intent(inout) :: eta !< Interface heights [Z ~> m]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, optional, intent(in) :: dZ_ref_eta !< The difference between the + !! reference heights for bathyT and + !! eta [Z ~> m], 0 by default. ! Local variables integer :: i, j, k, is, ie, js, je, nz, contractions, dilations real :: hTolerance = 0.1 !< Tolerance to exceed adjustment criteria [Z ~> m] - real :: hTmp, eTmp, dilate + real :: dilate ! A factor by which the column is dilated [nondim] + real :: dZ_ref ! The difference in the reference heights for G%bathyT and eta [Z ~> m] character(len=100) :: mesg is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke hTolerance = 0.1*US%m_to_Z + dZ_ref = 0.0 ; if (present(dZ_ref_eta)) dZ_ref = dZ_ref_eta contractions = 0 do j=js,je ; do i=is,ie - if (-eta(i,j,nz+1) > G%bathyT(i,j) + hTolerance) then - eta(i,j,nz+1) = -G%bathyT(i,j) + if (-eta(i,j,nz+1) > (G%bathyT(i,j) + dZ_ref) + hTolerance) then + eta(i,j,nz+1) = -(G%bathyT(i,j) + dZ_ref) contractions = contractions + 1 endif enddo ; enddo @@ -779,12 +784,12 @@ subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h) ! The whole column is dilated to accommodate deeper topography than ! the bathymetry would indicate. ! This should be... if ((G%mask2dt(i,j)*(eta(i,j,1)-eta(i,j,nz+1)) > 0.0) .and. & - if (-eta(i,j,nz+1) < G%bathyT(i,j) - hTolerance) then + if (-eta(i,j,nz+1) < (G%bathyT(i,j) + dZ_ref) - hTolerance) then dilations = dilations + 1 if (eta(i,j,1) <= eta(i,j,nz+1)) then - do k=1,nz ; h(i,j,k) = (eta(i,j,1) + G%bathyT(i,j)) / real(nz) ; enddo + do k=1,nz ; h(i,j,k) = (eta(i,j,1) + (G%bathyT(i,j) + dZ_ref)) / real(nz) ; enddo else - dilate = (eta(i,j,1) + G%bathyT(i,j)) / (eta(i,j,1) - eta(i,j,nz+1)) + dilate = (eta(i,j,1) + (G%bathyT(i,j) + dZ_ref)) / (eta(i,j,1) - eta(i,j,nz+1)) do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo endif do k=nz,2,-1 ; eta(i,j,K) = eta(i,j,K+1) + h(i,j,k) ; enddo @@ -1746,11 +1751,11 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, param_f type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic !! variables. real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: u !< The zonal velocity that is being - !! initialized [L T-1 ~> m s-1] + intent(in) :: u !< The zonal velocity that is being + !! initialized [L T-1 ~> m s-1] real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - intent(in) :: v !< The meridional velocity that is being - !! initialized [L T-1 ~> m s-1] + intent(in) :: v !< The meridional velocity that is being + !! initialized [L T-1 ~> m s-1] type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. type(sponge_CS), pointer :: Layer_CSp !< A pointer that is set to point to the control !! structure for this module (in layered mode). @@ -2182,9 +2187,9 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p ! calculate increments if input are full fields if (oda_inc) then ! input are increments - if (is_root_pe()) call MOM_error(NOTE,"incupd using increments fields ") + if (is_root_pe()) call MOM_mesg("incupd using increments fields ") else ! inputs are full fields - if (is_root_pe()) call MOM_error(NOTE,"incupd using full fields ") + if (is_root_pe()) call MOM_mesg("incupd using full fields ") call calc_oda_increments(h, tv, u, v, G, GV, US, oda_incupd_CSp) if (save_inc) then call output_oda_incupd_inc(Time, G, GV, param_file, oda_incupd_CSp, US) @@ -2319,6 +2324,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param real, dimension(:,:,:), allocatable, target :: temp_z, salt_z, mask_z real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zi ! Interface heights [Z ~> m]. + real, dimension(SZI_(G),SZJ_(G)) :: Z_bottom ! The (usually negative) height of the seafloor + ! relative to the surface [Z ~> m]. integer, dimension(SZI_(G),SZJ_(G)) :: nlevs real, dimension(SZI_(G)) :: press ! Pressures [R L2 T-2 ~> Pa]. @@ -2513,6 +2520,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param call pass_var(mask_z,G%Domain) call pass_var(rho_z,G%Domain) + do j=js,je ; do i=is,ie + Z_bottom(i,j) = -G%bathyT(i,j) + enddo ; enddo + ! Done with horizontal interpolation. ! Now remap to model coordinates if (useALEremapping) then @@ -2530,11 +2541,11 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param tmp_mask_in(i,j,1:kd) = mask_z(i,j,:) do k = 1, nkd if (tmp_mask_in(i,j,k)>0. .and. k<=kd) then - zBottomOfCell = max( z_edges_in(k+1), -G%bathyT(i,j) ) + zBottomOfCell = max( z_edges_in(k+1), Z_bottom(i,j)) tmpT1dIn(i,j,k) = temp_z(i,j,k) tmpS1dIn(i,j,k) = salt_z(i,j,k) elseif (k>1) then - zBottomOfCell = -G%bathyT(i,j) + zBottomOfCell = Z_bottom(i,j) tmpT1dIn(i,j,k) = tmpT1dIn(i,j,k-1) tmpS1dIn(i,j,k) = tmpS1dIn(i,j,k-1) else ! This next block should only ever be reached over land @@ -2544,7 +2555,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param h1(i,j,k) = GV%Z_to_H * (zTopOfCell - zBottomOfCell) zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k enddo - h1(i,j,kd) = h1(i,j,kd) + GV%Z_to_H * max(0., zTopOfCell + G%bathyT(i,j) ) + h1(i,j,kd) = h1(i,j,kd) + GV%Z_to_H * max(0., zTopOfCell - Z_bottom(i,j) ) ! The max here is in case the data data is shallower than model endif ! mask2dT enddo ; enddo @@ -2567,7 +2578,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param ! Build the target grid combining hTarget and topography zTopOfCell = 0. ; zBottomOfCell = 0. do k = 1, nz - zBottomOfCell = max( zTopOfCell - hTarget(k), -G%bathyT(i,j) ) + zBottomOfCell = max( zTopOfCell - hTarget(k), Z_bottom(i,j)) h(i,j,k) = GV%Z_to_H * (zTopOfCell - zBottomOfCell) zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k enddo @@ -2624,11 +2635,11 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param nkml = 0 ; if (separate_mixed_layer) nkml = GV%nkml - call find_interfaces(rho_z, z_in, kd, Rb, G%bathyT, zi, G, GV, US, nlevs, nkml, & + call find_interfaces(rho_z, z_in, kd, Rb, Z_bottom, zi, G, GV, US, nlevs, nkml, & Hmix_depth, eps_z, eps_rho, density_extrap_bug) if (correct_thickness) then - call adjustEtaToFitBathymetry(G, GV, US, zi, h) + call adjustEtaToFitBathymetry(G, GV, US, zi, h, dZ_ref_eta=0.0) else do k=nz,1,-1 ; do j=js,je ; do i=is,ie if (zi(i,j,K) < (zi(i,j,K+1) + GV%Angstrom_Z)) then @@ -2640,7 +2651,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param enddo ; enddo ; enddo inconsistent=0 do j=js,je ; do i=is,ie - if (abs(zi(i,j,nz+1) + G%bathyT(i,j)) > 1.0*US%m_to_Z) & + if (abs(zi(i,j,nz+1) - Z_bottom(i,j)) > 1.0*US%m_to_Z) & inconsistent = inconsistent + 1 enddo ; enddo call sum_across_PEs(inconsistent) @@ -2710,7 +2721,7 @@ end subroutine MOM_temp_salt_initialize_from_Z !> Find interface positions corresponding to interpolated depths in a density profile -subroutine find_interfaces(rho, zin, nk_data, Rb, depth, zi, G, GV, US, nlevs, nkml, hml, & +subroutine find_interfaces(rho, zin, nk_data, Rb, Z_bot, zi, G, GV, US, nlevs, nkml, hml, & eps_z, eps_rho, density_extrap_bug) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure @@ -2720,7 +2731,8 @@ subroutine find_interfaces(rho, zin, nk_data, Rb, depth, zi, G, GV, US, nlevs, n real, dimension(nk_data), intent(in) :: zin !< Input data levels [Z ~> m]. real, dimension(SZK_(GV)+1), intent(in) :: Rb !< target interface densities [R ~> kg m-3] real, dimension(SZI_(G),SZJ_(G)), & - intent(in) :: depth !< ocean depth [Z ~> m]. + intent(in) :: Z_bot !< The (usually negative) height of the seafloor + !! relative to the surface [Z ~> m]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & intent(out) :: zi !< The returned interface heights [Z ~> m] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -2808,15 +2820,15 @@ subroutine find_interfaces(rho, zin, nk_data, Rb, depth, zi, G, GV, US, nlevs, n ! Linearly interpolate to find the depth, zi_, where Rb would be found. slope = (zin(k_int+1) - zin(k_int)) / max(rho_(k_int+1) - rho_(k_int), eps_rho) zi_(K) = -1.0*(zin(k_int) + slope*(Rb(K)-rho_(k_int))) - zi_(K) = min(max(zi_(K), -depth(i,j)), -1.0*hml) + zi_(K) = min(max(zi_(K), Z_bot(i,j)), -1.0*hml) enddo - zi_(nz+1) = -depth(i,j) + zi_(nz+1) = Z_bot(i,j) if (nkml > 0) then ; do K=2,nkml+1 - zi_(K) = max(hml*((1.0-real(K))/real(nkml)), -depth(i,j)) + zi_(K) = max(hml*((1.0-real(K))/real(nkml)), Z_bot(i,j)) enddo ; endif do K=nz,max(nkml+2,2),-1 if (zi_(K) < zi_(K+1) + eps_Z) zi_(K) = zi_(K+1) + eps_Z - if (zi_(K) > -1.0*hml) zi_(K) = max(-1.0*hml, -depth(i,j)) + if (zi_(K) > -1.0*hml) zi_(K) = max(-1.0*hml, Z_bot(i,j)) enddo do K=1,nz+1 @@ -2864,7 +2876,7 @@ subroutine MOM_state_init_tests(G, GV, US, tv) S_t(k) = 35. - (0. * I_z_scale)*e(k) S(k) = 35. + (0. * I_z_scale)*z(k) S_b(k) = 35. - (0. * I_z_scale)*e(k+1) - call calculate_density(0.5*(T_t(k)+T_b(k)), 0.5*(S_t(k)+S_b(k)), -GV%Rho0*GV%g_Earth*US%m_to_Z*z(k), & + call calculate_density(0.5*(T_t(k)+T_b(k)), 0.5*(S_t(k)+S_b(k)), -GV%Rho0*GV%g_Earth*z(k), & rho(k), tv%eqn_of_state) P_tot = P_tot + GV%g_Earth * rho(k) * GV%H_to_Z*h(k) enddo From 62446ec0ff60ea30de920cb2e321816355f4f468 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 15 Aug 2021 05:10:19 -0400 Subject: [PATCH 07/11] Cleanup of comments in MOM_ALE_sponge.F90 Extensive inconsequential cleanup in MOM_ALE_sponge.F90, including the elimination of a dozen unnecessary index-range elements of the ALE_sponge_CS and modifying a number of comments to be much more descriptive. This included the correction of numerous spelling errors and other typos in comments. All answers are bitwise identical. --- .../vertical/MOM_ALE_sponge.F90 | 148 ++++++++---------- 1 file changed, 66 insertions(+), 82 deletions(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index e122452368..419b012387 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -47,7 +47,7 @@ module MOM_ALE_sponge module procedure set_up_ALE_sponge_vel_field_varying end interface -!> Ddetermine the number of points which are within sponges in this computational domain. +!> Determine the number of points which are within sponges in this computational domain. !! !! Only points that have positive values of Iresttime and which mask2dT indicates are ocean !! points are included in the sponges. It also stores the target interface heights. @@ -90,31 +90,19 @@ module MOM_ALE_sponge !> ALE sponge control structure type, public :: ALE_sponge_CS ; private integer :: nz !< The total number of layers. - integer :: nz_data !< The total number of arbritary layers (used by older code). - integer :: isc !< The starting i-index of the computational domain at h. - integer :: iec !< The ending i-index of the computational domain at h. - integer :: jsc !< The starting j-index of the computational domain at h. - integer :: jec !< The ending j-index of the computational domain at h. - integer :: IscB !< The starting I-index of the computational domain at u/v. - integer :: IecB !< The ending I-index of the computational domain at u/v. - integer :: JscB !< The starting J-index of the computational domain at u/v. - integer :: JecB !< The ending J-index of the computational domain at h. - integer :: isd !< The starting i-index of the data domain at h. - integer :: ied !< The ending i-index of the data domain at h. - integer :: jsd !< The starting j-index of the data domain at h. - integer :: jed !< The ending j-index of the data domain at h. + integer :: nz_data !< The total number of arbitrary layers (used by older code). integer :: num_col !< The number of sponge tracer points within the computational domain. integer :: num_col_u !< The number of sponge u-points within the computational domain. integer :: num_col_v !< The number of sponge v-points within the computational domain. integer :: fldno = 0 !< The number of fields which have already been !! registered by calls to set_up_sponge_field logical :: sponge_uv !< Control whether u and v are included in sponge - integer, pointer :: col_i(:) => NULL() !< Array of the i-indicies of each tracer columns being damped. - integer, pointer :: col_j(:) => NULL() !< Array of the j-indicies of each tracer columns being damped. - integer, pointer :: col_i_u(:) => NULL() !< Array of the i-indicies of each u-columns being damped. - integer, pointer :: col_j_u(:) => NULL() !< Array of the j-indicies of each u-columns being damped. - integer, pointer :: col_i_v(:) => NULL() !< Array of the i-indicies of each v-columns being damped. - integer, pointer :: col_j_v(:) => NULL() !< Array of the j-indicies of each v-columns being damped. + integer, pointer :: col_i(:) => NULL() !< Array of the i-indices of each tracer column being damped. + integer, pointer :: col_j(:) => NULL() !< Array of the j-indices of each tracer column being damped. + integer, pointer :: col_i_u(:) => NULL() !< Array of the i-indices of each u-column being damped. + integer, pointer :: col_j_u(:) => NULL() !< Array of the j-indices of each u-column being damped. + integer, pointer :: col_i_v(:) => NULL() !< Array of the i-indices of each v-column being damped. + integer, pointer :: col_j_v(:) => NULL() !< Array of the j-indices of each v-column being damped. real, pointer :: Iresttime_col(:) => NULL() !< The inverse restoring time of each tracer column [T-1 ~> s-1]. real, pointer :: Iresttime_col_u(:) => NULL() !< The inverse restoring time of each u-column [T-1 ~> s-1]. @@ -124,8 +112,8 @@ module MOM_ALE_sponge type(p2d) :: Ref_val(MAX_FIELDS_) !< The values to which the fields are damped. type(p2d) :: Ref_val_u !< The values to which the u-velocities are damped. type(p2d) :: Ref_val_v !< The values to which the v-velocities are damped. - type(p3d) :: var_u !< Pointer to the u velocities. that are being damped. - type(p3d) :: var_v !< Pointer to the v velocities. that are being damped. + type(p3d) :: var_u !< Pointer to the u velocities that are being damped. + type(p3d) :: var_v !< Pointer to the v velocities that are being damped. type(p2d) :: Ref_h !< Grid on which reference data is provided (older code). type(p2d) :: Ref_hu !< u-point grid on which reference data is provided (older code). type(p2d) :: Ref_hv !< v-point grid on which reference data is provided (older code). @@ -137,7 +125,7 @@ module MOM_ALE_sponge logical :: remap_answers_2018 !< If true, use the order of arithmetic and expressions that !! recover the answers for remapping from the end of 2018. !! Otherwise, use more robust forms of the same expressions. - logical :: hor_regrid_answers_2018 !< If true, use the order of arithmetic for horizonal regridding + logical :: hor_regrid_answers_2018 !< If true, use the order of arithmetic for horizontal regridding !! that recovers the answers from the end of 2018. Otherwise, use !! rotationally symmetric forms of the same expressions. @@ -241,9 +229,6 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, CS%time_varying_sponges = .false. CS%nz = GV%ke - CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec - CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed - CS%iscB = G%iscB ; CS%iecB = G%iecB; CS%jscB = G%jscB ; CS%jecB = G%jecB ! number of columns to be restored CS%num_col = 0 ; CS%fldno = 0 @@ -265,7 +250,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, col = col +1 endif enddo ; enddo - ! same for total number of arbritary layers and correspondent data + ! same for total number of arbitrary layers and correspondent data CS%nz_data = nz_data allocate(CS%Ref_h%p(CS%nz_data,CS%num_col)) do col=1,CS%num_col ; do K=1,CS%nz_data @@ -295,11 +280,11 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, if (present(Iresttime_u_in)) then Iresttime_u(:,:) = Iresttime_u_in(:,:) else - do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec ; do I=G%iscB,G%iecB Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) enddo ; enddo endif - do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec ; do I=G%iscB,G%iecB if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) & CS%num_col_u = CS%num_col_u + 1 enddo ; enddo @@ -312,7 +297,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, ! Store the column indices and restoring rates in the CS structure col = 1 - do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec ; do I=G%iscB,G%iecB if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then CS%col_i_u(col) = I ; CS%col_j_u(col) = j CS%Iresttime_col_u(col) = Iresttime_u(I,j) @@ -320,7 +305,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, endif enddo ; enddo - ! same for total number of arbritary layers and correspondent data + ! same for total number of arbitrary layers and correspondent data allocate(CS%Ref_hu%p(CS%nz_data,CS%num_col_u)) do col=1,CS%num_col_u I = CS%col_i_u(col) ; j = CS%col_j_u(col) @@ -339,11 +324,11 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, if (present(Iresttime_v_in)) then Iresttime_v(:,:) = Iresttime_v_in(:,:) else - do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB; do i=G%isc,G%iec Iresttime_v(i,J) = 0.5 * (Iresttime(i,j) + Iresttime(i,j+1)) enddo ; enddo endif - do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB; do i=G%isc,G%iec if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) & CS%num_col_v = CS%num_col_v + 1 enddo ; enddo @@ -356,7 +341,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, ! pass indices, restoring time to the CS structure col = 1 - do J=CS%jscB,CS%jecB ; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB ; do i=G%isc,G%iec if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) then CS%col_i_v(col) = i ; CS%col_j_v(col) = j CS%Iresttime_col_v(col) = Iresttime_v(i,j) @@ -364,7 +349,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, endif enddo ; enddo - ! same for total number of arbritary layers and correspondent data + ! same for total number of arbitrary layers and correspondent data allocate(CS%Ref_hv%p(CS%nz_data,CS%num_col_v)) do col=1,CS%num_col_v i = CS%col_i_v(col) ; J = CS%col_j_v(col) @@ -430,7 +415,7 @@ subroutine get_ALE_sponge_thicknesses(G, data_h, sponge_mask, CS) end subroutine get_ALE_sponge_thicknesses -!> This subroutine determines the number of points which are to be restoref in the computational +!> This subroutine determines the number of points which are to be restored in the computational !! domain. Only points that have positive values of Iresttime and which mask2dT indicates are ocean !! points are included in the sponges. subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Iresttime_u_in, Iresttime_v_in) @@ -510,9 +495,6 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest CS%time_varying_sponges = .true. CS%nz = GV%ke - CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec - CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed - CS%iscB = G%iscB ; CS%iecB = G%iecB; CS%jscB = G%jscB ; CS%jecB = G%jecB ! number of columns to be restored CS%num_col = 0 ; CS%fldno = 0 @@ -551,12 +533,12 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest if (present(Iresttime_u_in)) then Iresttime_u(:,:) = Iresttime_u_in(:,:) else - do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec ; do I=G%iscB,G%iecB Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) enddo ; enddo endif CS%num_col_u = 0 ; - do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec; do I=G%iscB,G%iecB if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) & CS%num_col_u = CS%num_col_u + 1 enddo ; enddo @@ -566,14 +548,14 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest allocate(CS%col_j_u(CS%num_col_u)) ; CS%col_j_u = 0 ! pass indices, restoring time to the CS structure col = 1 - do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec ; do I=G%iscB,G%iecB if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then CS%col_i_u(col) = i ; CS%col_j_u(col) = j CS%Iresttime_col_u(col) = Iresttime_u(i,j) col = col + 1 endif enddo ; enddo - ! same for total number of arbritary layers and correspondent data + ! same for total number of arbitrary layers and correspondent data endif total_sponge_cols_u = CS%num_col_u call sum_across_PEs(total_sponge_cols_u) @@ -583,12 +565,12 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest if (present(Iresttime_v_in)) then Iresttime_v(:,:) = Iresttime_v_in(:,:) else - do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB; do i=G%isc,G%iec Iresttime_v(i,J) = 0.5 * (Iresttime(i,j) + Iresttime(i,j+1)) enddo ; enddo endif CS%num_col_v = 0 ; - do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB; do i=G%isc,G%iec if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) & CS%num_col_v = CS%num_col_v + 1 enddo ; enddo @@ -598,7 +580,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest allocate(CS%col_j_v(CS%num_col_v)) ; CS%col_j_v = 0 ! pass indices, restoring time to the CS structure col = 1 - do J=CS%jscB,CS%jecB ; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB ; do i=G%isc,G%iec if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) then CS%col_i_v(col) = i ; CS%col_j_v(col) = j CS%Iresttime_col_v(col) = Iresttime_v(i,j) @@ -630,16 +612,16 @@ subroutine init_ALE_sponge_diags(Time, G, diag, CS, US) CS%id_sp_tendency(1) = -1 CS%id_sp_tendency(1) = register_diag_field('ocean_model', 'sp_tendency_temp', diag%axesTL, Time, & - 'Time tendency due to temperature restoring', 'degC s-1',conversion=US%s_to_T) + 'Time tendency due to temperature restoring', 'degC s-1', conversion=US%s_to_T) CS%id_sp_tendency(2) = -1 CS%id_sp_tendency(2) = register_diag_field('ocean_model', 'sp_tendency_salt', diag%axesTL, Time, & - 'Time tendency due to salinity restoring', 'g kg-1 s-1',conversion=US%s_to_T) + 'Time tendency due to salinity restoring', 'g kg-1 s-1', conversion=US%s_to_T) CS%id_sp_u_tendency = -1 CS%id_sp_u_tendency = register_diag_field('ocean_model', 'sp_tendency_u', diag%axesCuL, Time, & - 'Zonal acceleration due to sponges', 'm s-2',conversion=US%L_T2_to_m_s2) + 'Zonal acceleration due to sponges', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_sp_v_tendency = -1 CS%id_sp_v_tendency = register_diag_field('ocean_model', 'sp_tendency_v', diag%axesCvL, Time, & - 'Meridional acceleration due to sponges', 'm s-2',conversion=US%L_T2_to_m_s2) + 'Meridional acceleration due to sponges', 'm s-2', conversion=US%L_T2_to_m_s2) end subroutine init_ALE_sponge_diags @@ -718,19 +700,18 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed CS%fldno = CS%fldno + 1 if (CS%fldno > MAX_FIELDS_) then - write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease & - &the number of fields to be damped in the call to & - &initialize_ALE_sponge." )') CS%fldno + write(mesg, '("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease "//& + "the number of fields to be damped in the call to initialize_ALE_sponge." )') CS%fldno call MOM_error(FATAL,"set_up_ALE_sponge_field: "//mesg) endif - ! get a unique time interp id for this field. If sponge data is ongrid, then setup + ! get a unique time interp id for this field. If sponge data is on-grid, then setup ! to only read on the computational domain if (CS%spongeDataOngrid) then CS%Ref_val(CS%fldno)%id = init_external_field(filename, fieldname, MOM_domain=G%Domain) else CS%Ref_val(CS%fldno)%id = init_external_field(filename, fieldname) endif - fld_sz(1:4)=-1 + fld_sz(1:4) = -1 call get_external_field_info(CS%Ref_val(CS%fldno)%id, size=fld_sz) nz_data = fld_sz(3) CS%Ref_val(CS%fldno)%nz_data = nz_data !< individual sponge fields may reside on a different vertical grid @@ -748,17 +729,19 @@ end subroutine set_up_ALE_sponge_field_varying !> This subroutine stores the reference profile at u and v points for the variable !! whose address is given by u_ptr and v_ptr. subroutine set_up_ALE_sponge_vel_field_fixed(u_val, v_val, G, GV, u_ptr, v_ptr, CS) - type(ocean_grid_type), intent(in) :: G !< Grid structure (in). + type(ocean_grid_type), intent(in) :: G !< Grid structure (in). type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - type(ALE_sponge_CS), pointer :: CS !< Sponge structure (in/out). + type(ALE_sponge_CS), pointer :: CS !< Sponge structure (in/out). real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: u_val !< u field to be used in the sponge, it has arbritary number of layers but - !! not to exceed the total number of model layers + intent(in) :: u_val !< u field to be used in the sponge [L T-1 ~> m s-1], + !! it is provided on its own vertical grid that may + !! have fewer layers than the model itself, but not more. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - intent(in) :: v_val !< v field to be used in the sponge, it has arbritary number of layers but - !! not to exceed the number of model layers - real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u pointer to the field to be damped - real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v pointer to the field to be damped + intent(in) :: v_val !< v field to be used in the sponge [L T-1 ~> m s-1], + !! it is provided on its own vertical grid that may + !! have fewer layers than the model itself, but not more. + real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u-field to be damped [L T-1 ~> m s-1] + real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v-field to be damped [L T-1 ~> m s-1] integer :: j, k, col, fld_sz(4) character(len=256) :: mesg ! String for error messages @@ -794,15 +777,16 @@ subroutine set_up_ALE_sponge_vel_field_varying(filename_u, fieldname_u, filename character(len=*), intent(in) :: filename_v !< File name for v field character(len=*), intent(in) :: fieldname_v !< Name of v variable in file type(time_type), intent(in) :: Time !< Model time - type(ocean_grid_type), intent(in) :: G !< Ocean grid (in) - type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(ALE_sponge_CS), pointer :: CS !< Sponge structure (in/out). - real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u pointer to the field to be damped (in). - real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v pointer to the field to be damped (in). + type(ocean_grid_type), intent(in) :: G !< Ocean grid (in) + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(ALE_sponge_CS), pointer :: CS !< Sponge structure (in/out). + real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u-field to be damped [L T-1 ~> m s-1] + real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v-field to be damped [L T-1 ~> m s-1] + ! Local variables - real, allocatable, dimension(:,:,:) :: u_val !< U field to be used in the sponge. - real, allocatable, dimension(:,:,:) :: v_val !< V field to be used in the sponge. + real, allocatable, dimension(:,:,:) :: u_val !< U field to be used in the sponge [L T-1 ~> m s-1]. + real, allocatable, dimension(:,:,:) :: v_val !< V field to be used in the sponge [L T-1 ~> m s-1]. real, allocatable, dimension(:), target :: z_in, z_edges_in real :: missing_value @@ -892,10 +876,13 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) real, dimension(:), allocatable :: tmpT1d integer :: c, m, nkmb, i, j, k, is, ie, js, je, nz, nz_data integer :: col, total_sponge_cols - real, allocatable, dimension(:), target :: z_in, z_edges_in - real :: missing_value, Idt - real :: h_neglect, h_neglect_edge - real :: zTopOfCell, zBottomOfCell ! Heights [Z ~> m]. + real, allocatable, dimension(:), target :: z_in ! The depths (positive downward) in the input file [Z ~> m] + real, allocatable, dimension(:), target :: z_edges_in ! The depths (positive downward) of the + ! edges in the input file [Z ~> m] + real :: missing_value + real :: Idt ! The inverse of the timestep [T-1 ~> s-1] + real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] + real :: zTopOfCell, zBottomOfCell ! Interface heights (positive upward) in the input dataset [Z ~> m]. integer :: nPoints is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -1020,7 +1007,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) call pass_var(sp_val, G%Domain) call pass_var(mask_z, G%Domain) - do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB + do j=G%jsc,G%jec; do I=G%iscB,G%iecB sp_val_u(I,j,1:nz_data) = 0.5*(sp_val(i,j,1:nz_data)+sp_val(i+1,j,1:nz_data)) mask_u(I,j,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i+1,j,1:nz_data)) enddo ; enddo @@ -1068,7 +1055,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) answers_2018=CS%hor_regrid_answers_2018) call pass_var(sp_val, G%Domain) call pass_var(mask_z, G%Domain) - do J=CS%jscB,CS%jecB; do i=CS%isc,CS%iec + do J=G%jscB,G%jecB; do i=G%isc,G%iec sp_val_v(i,J,1:nz_data) = 0.5*(sp_val(i,j,1:nz_data)+sp_val(i,j+1,1:nz_data)) mask_v(i,J,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i,j+1,1:nz_data)) enddo ; enddo @@ -1107,7 +1094,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) nz_data = CS%Ref_val_u%nz_data allocate(tmp_val2(nz_data)) if (CS%id_sp_u_tendency > 0) then - allocate(tmp_u(G%isdB:G%iedB,G%jsd:G%jed,nz));tmp_u(:,:,:)=0.0 + allocate(tmp_u(G%isdB:G%iedB,G%jsd:G%jed,nz)) ; tmp_u(:,:,:)=0.0 endif ! u points do c=1,CS%num_col_u @@ -1137,7 +1124,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) endif ! v points if (CS%id_sp_v_tendency > 0) then - allocate(tmp_v(G%isd:G%ied,G%jsdB:G%jedB,nz));tmp_v(:,:,:)=0.0 + allocate(tmp_v(G%isd:G%ied,G%jsdB:G%jedB,nz)) ; tmp_v(:,:,:)=0.0 endif nz_data = CS%Ref_val_v%nz_data allocate(tmp_val2(nz_data)) @@ -1170,9 +1157,6 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) deallocate(tmp_val2) endif - - - end subroutine apply_ALE_sponge !> Rotate the ALE sponge fields from the input to the model index map. From c72d441ead6bb722378cd36350c74c968b13d4eb Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 15 Aug 2021 05:12:04 -0400 Subject: [PATCH 08/11] (*)Fix dimensionally inconsistent MEKE beta calcs Corrected dimensional inconsistencies in the negligible thicknesses in the denominators of the expressions for the topographic betas in MEKE_equilibrium and MEKE_lengthScales. This could change answers in strange cases, but seems unlikely to do so (partly because it is in a max expression, and not added), and did not change any answers in the MOM6-examples test suite. --- src/parameterizations/lateral/MOM_MEKE.F90 | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 762b2edaea..9441ab7107 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -637,7 +637,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h end subroutine step_forward_MEKE -!> Calculates the equilibrium solutino where the source depends only on MEKE diffusivity +!> Calculates the equilibrium solution where the source depends only on MEKE diffusivity !! and there is no lateral diffusion of MEKE. !! Results is in MEKE%MEKE. subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass) @@ -667,6 +667,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m real :: resid, ResMin, ResMax ! Residuals [L2 T-3 ~> W kg-1] real :: FatH ! Coriolis parameter at h points; to compute topographic beta [T-1 ~> s-1] real :: beta_topo_x, beta_topo_y ! Topographic PV gradients in x and y [T-1 L-1 ~> s-1 m-1] + real :: dZ_neglect ! A negligible change in height [Z ~> m] integer :: i, j, is, ie, js, je, n1, n2 real :: tolerance ! Width of EKE bracket [L2 T-2 ~> m2 s-2]. logical :: useSecant, debugIteration @@ -680,6 +681,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m Ubg2 = CS%MEKE_Uscale**2 cd2 = CS%cdrag**2 tolerance = 1.0e-12*US%m_s_to_L_T**2 + dZ_neglect = GV%H_to_Z*GV%H_subroundoff !$OMP do do j=js,je ; do i=is,ie @@ -701,14 +703,14 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m ! of the water column thickness instead of the bathymetric depth. beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & (G%bathyT(i+1,j)-G%bathyT(i,j)) * G%IdxCu(I,j) & - / max(G%bathyT(i+1,j),G%bathyT(i,j), GV%H_subroundoff) & + / max(G%bathyT(i+1,j),G%bathyT(i,j), dZ_neglect) & + (G%bathyT(i,j)-G%bathyT(i-1,j)) * G%IdxCu(I-1,j) & - / max(G%bathyT(i,j),G%bathyT(i-1,j), GV%H_subroundoff) ) + / max(G%bathyT(i,j),G%bathyT(i-1,j), dZ_neglect) ) beta_topo_y = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & (G%bathyT(i,j+1)-G%bathyT(i,j)) * G%IdyCv(i,J) & - / max(G%bathyT(i,j+1),G%bathyT(i,j), GV%H_subroundoff) + & + / max(G%bathyT(i,j+1),G%bathyT(i,j), dZ_neglect) + & (G%bathyT(i,j)-G%bathyT(i,j-1)) * G%IdyCv(i,J-1) & - / max(G%bathyT(i,j),G%bathyT(i,j-1), GV%H_subroundoff) ) + / max(G%bathyT(i,j),G%bathyT(i,j-1), dZ_neglect) ) endif beta = sqrt((G%dF_dx(i,j) + beta_topo_x)**2 + & (G%dF_dy(i,j) + beta_topo_y)**2 ) @@ -853,9 +855,11 @@ subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, & real :: SN ! The local Eady growth rate [T-1 ~> s-1] real :: FatH ! Coriolis parameter at h points [T-1 ~> s-1] real :: beta_topo_x, beta_topo_y ! Topographic PV gradients in x and y [T-1 L-1 ~> s-1 m-1] + real :: dZ_neglect ! A negligible change in height [Z ~> m] integer :: i, j, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + dZ_neglect = GV%H_to_Z*GV%H_subroundoff !$OMP do do j=js,je ; do i=is,ie @@ -878,14 +882,14 @@ subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, & ! of the water column thickness instead of the bathymetric depth. beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & (G%bathyT(i+1,j)-G%bathyT(i,j)) * G%IdxCu(I,j) & - / max(G%bathyT(i+1,j),G%bathyT(i,j), GV%H_subroundoff) & + / max(G%bathyT(i+1,j),G%bathyT(i,j), dZ_neglect) & + (G%bathyT(i,j)-G%bathyT(i-1,j)) * G%IdxCu(I-1,j) & - / max(G%bathyT(i,j),G%bathyT(i-1,j), GV%H_subroundoff) ) + / max(G%bathyT(i,j),G%bathyT(i-1,j), dZ_neglect) ) beta_topo_y = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & (G%bathyT(i,j+1)-G%bathyT(i,j)) * G%IdyCv(i,J) & - / max(G%bathyT(i,j+1),G%bathyT(i,j), GV%H_subroundoff) + & + / max(G%bathyT(i,j+1),G%bathyT(i,j), dZ_neglect) + & (G%bathyT(i,j)-G%bathyT(i,j-1)) * G%IdyCv(i,J-1) & - / max(G%bathyT(i,j),G%bathyT(i,j-1), GV%H_subroundoff) ) + / max(G%bathyT(i,j),G%bathyT(i,j-1), dZ_neglect) ) endif beta = sqrt((G%dF_dx(i,j) + beta_topo_x)**2 + & (G%dF_dy(i,j) + beta_topo_y)**2 ) From 5017bf6c65702ef45ce843cdeed376da1c0cbecf Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 15 Aug 2021 05:41:48 -0400 Subject: [PATCH 09/11] Added comments highlighting bugs Added comments denoted with "###" indicating bugs or conceptual errors in update_OBC_segment_data, horizontal_viscosity, thickness_diffuse and wave_speed. Only comments and the case of some indices are changed in this commit, and all answers are bitwise identical. Actually correcting these bugs would probably change answers in some cases. --- src/core/MOM_open_boundary.F90 | 21 ++++++++++--------- src/diagnostics/MOM_wave_speed.F90 | 3 ++- .../lateral/MOM_hor_visc.F90 | 10 +++++---- .../lateral/MOM_thickness_diffuse.F90 | 7 ++++++- 4 files changed, 25 insertions(+), 16 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 61e20d14a6..318d10008c 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -3802,29 +3802,31 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) ishift=0;jshift=0 if (segment%is_E_or_W) then allocate(normal_trans_bt(segment%HI%IsdB:segment%HI%IedB,segment%HI%jsd:segment%HI%jed)) - normal_trans_bt(:,:)=0.0 + normal_trans_bt(:,:) = 0.0 if (segment%direction == OBC_DIRECTION_W) ishift=1 I=segment%HI%IsdB do j=segment%HI%jsd,segment%HI%jed - segment%Cg(I,j) = sqrt(GV%g_prime(1)*G%bathyT(i+ishift,j)) - segment%Htot(I,j)=0.0 + segment%Htot(I,j) = 0.0 do k=1,GV%ke segment%h(I,j,k) = h(i+ishift,j,k) - segment%Htot(I,j)=segment%Htot(I,j)+segment%h(I,j,k) + segment%Htot(I,j) = segment%Htot(I,j) + segment%h(I,j,k) enddo + segment%Cg(I,j) = sqrt(GV%g_prime(1)*G%bathyT(i+ishift,j)) + !### This should be: segment%Cg(I,j) = sqrt(GV%g_prime(1)*segment%Htot(I,j)*GV%H_to_Z) enddo else! (segment%direction == OBC_DIRECTION_N .or. segment%direction == OBC_DIRECTION_S) allocate(normal_trans_bt(segment%HI%isd:segment%HI%ied,segment%HI%JsdB:segment%HI%JedB)) - normal_trans_bt(:,:)=0.0 + normal_trans_bt(:,:) = 0.0 if (segment%direction == OBC_DIRECTION_S) jshift=1 J=segment%HI%JsdB do i=segment%HI%isd,segment%HI%ied - segment%Cg(i,J) = sqrt(GV%g_prime(1)*G%bathyT(i,j+jshift)) - segment%Htot(i,J)=0.0 + segment%Htot(i,J) = 0.0 do k=1,GV%ke segment%h(i,J,k) = h(i,j+jshift,k) - segment%Htot(i,J)=segment%Htot(i,J)+segment%h(i,J,k) + segment%Htot(i,J) = segment%Htot(i,J) + segment%h(i,J,k) enddo + segment%Cg(i,J) = sqrt(GV%g_prime(1)*G%bathyT(i,j+jshift)) + !### This should be: segment%Cg(i,J) = sqrt(GV%g_prime(1)*segment%Htot(i,J)*GV%H_to_Z) enddo endif @@ -4715,7 +4717,7 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) integer :: i, j integer :: l_seg logical :: fatal_error = .False. - real :: min_depth + real :: min_depth ! The minimum depth for ocean points [Z ~> m] integer, parameter :: cin = 3, cout = 4, cland = -1, cedge = -2 character(len=256) :: mesg ! Message for error messages. type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list @@ -4730,7 +4732,6 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) allocate(color(G%isd:G%ied, G%jsd:G%jed)) ; color = 0 allocate(color2(G%isd:G%ied, G%jsd:G%jed)) ; color2 = 0 - ! Paint a frame around the outside. do j=G%jsd,G%jed color(G%isd,j) = cedge diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 035386f92d..d363b185f8 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -444,7 +444,8 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_ hw = 0.5*(Hc(k-1)+Hc(k)) gp = gprime(K) if (l_mono_N2_column_fraction>0. .or. l_mono_N2_depth>=0.) then - if ( ((G%bathyT(i,j)-sum_hc < l_mono_N2_column_fraction*G%bathyT(i,j)) .or. & + !### Change to: if ( ((htot(i) - sum_hc < l_mono_N2_column_fraction*htot(i)) .or. & ) ) + if ( ((G%bathyT(i,j) - sum_hc < l_mono_N2_column_fraction*G%bathyT(i,j)) .or. & ((l_mono_N2_depth >= 0.) .and. (sum_hc > l_mono_N2_depth))) .and. & (L2_to_Z2*gp > N2min*hw) ) then ! Filters out regions where N2 increases with depth but only in a lower fraction diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index c588a1faa4..b4f857dec4 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -496,7 +496,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo do J=js-2,Jeq+1 ; do I=is-2,Ieq+1 - grad_vel_mag_bt_q(I,J) = boundary_mask_q(I,J) * (dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + & + grad_vel_mag_bt_q(I,J) = boundary_mask_q(I,J) * (dvdx_bt(I,J)**2 + dudy_bt(I,J)**2 + & (0.25*((dudx_bt(i,j)+dudx_bt(i+1,j+1))+(dudx_bt(i,j+1)+dudx_bt(i+1,j))))**2 + & (0.25*((dvdy_bt(i,j)+dvdy_bt(i+1,j+1))+(dvdy_bt(i,j+1)+dvdy_bt(i+1,j))))**2) enddo ; enddo @@ -1389,7 +1389,7 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 if (grad_vel_mag_bt_h(i,j)>0) then - GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * & + GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j) / CS%GME_h0, 1.0)**2) * & (0.25*(KH_u_GME(I,j,k)+KH_u_GME(I-1,j,k)+KH_v_GME(i,J,k)+KH_v_GME(i,J-1,k))) else GME_coeff = 0.0 @@ -1405,8 +1405,10 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, enddo ; enddo do J=js-1,Jeq ; do I=is-1,Ieq - if (grad_vel_mag_bt_q(i,j)>0) then - GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j)/CS%GME_h0,1.0)**2) * & + if (grad_vel_mag_bt_q(I,J)>0) then + !### This expression is not rotationally invariant - bathyT is to the SW of q points, + ! and it needs parentheses in the sum of the 4 diffusivities. + GME_coeff = CS%GME_efficiency * (MIN(G%bathyT(i,j) / CS%GME_h0, 1.0)**2) * & (0.25*(KH_u_GME(I,j,k)+KH_u_GME(I,j+1,k)+KH_v_GME(i,J,k)+KH_v_GME(i+1,J,k))) else GME_coeff = 0.0 diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index da62ffc6b7..3b3d72576c 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -481,14 +481,19 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp do j=js,je ; do I=is-1,ie hu(I,j) = 2.0*h(i,j,k)*h(i+1,j,k)/(h(i,j,k)+h(i+1,j,k)+h_neglect) if (hu(I,j) /= 0.0) hu(I,j) = 1.0 + !### The same result would be accomplished with the following without a division: + ! hu(I,j) = 0.0 ; if (h(i,j,k)*h(i+1,j,k) /= 0.0) hu(I,j) = 1.0 KH_u_lay(I,j) = 0.5*(KH_u(I,j,k)+KH_u(I,j,k+1)) enddo ; enddo do J=js-1,je ; do i=is,ie hv(i,J) = 2.0*h(i,j,k)*h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h_neglect) if (hv(i,J) /= 0.0) hv(i,J) = 1.0 + !### The same result would be accomplished with the following without a division: + ! hv(i,J) = 0.0 ; if (h(i,j,k)*h(i,j+1,k) /= 0.0) hv(i,J) = 1.0 KH_v_lay(i,J) = 0.5*(KH_v(i,J,k)+KH_v(i,J,k+1)) enddo ; enddo ! diagnose diffusivity at T-point + !### Because hu and hv are nondimensional here, the denominator is dimensionally inconsistent. do j=js,je ; do i=is,ie Kh_t(i,j,k) = ((hu(I-1,j)*KH_u_lay(i-1,j)+hu(I,j)*KH_u_lay(I,j)) & +(hv(i,J-1)*KH_v_lay(i,J-1)+hv(i,J)*KH_v_lay(i,J))) & @@ -505,7 +510,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp enddo do j=js,je ; do i=is,ie - MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(1.0,G%bathyT(i,j)) + MEKE%Kh_diff(i,j) = GV%H_to_Z * MEKE%Kh_diff(i,j) / MAX(1.0*US%m_to_Z, G%bathyT(i,j)) enddo ; enddo endif From 3cda6fd68d4105e623c0981077aca6588733a422 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 17 Aug 2021 07:19:01 -0400 Subject: [PATCH 10/11] +Pass depth_tot to initialization routines Added depth_tot arguments to explicitly pass the total depth to various initialization routines for the thicknesses, sponges, or temperatures. These routines had previously used G%bathT for this role, but this change prepares to allow the reference depth for the bathymetry to change without requiring these routines to be changed to accommodate it. All answers are bitwise identical, but there are new arguments to about two dozen routines. --- .../MOM_state_initialization.F90 | 109 ++++++++++-------- src/user/BFB_initialization.F90 | 6 +- src/user/DOME2d_initialization.F90 | 24 ++-- src/user/DOME_initialization.F90 | 38 +++--- src/user/ISOMIP_initialization.F90 | 30 +++-- src/user/Kelvin_initialization.F90 | 44 ++++--- src/user/Neverworld_initialization.F90 | 8 +- src/user/Phillips_initialization.F90 | 6 +- src/user/RGC_initialization.F90 | 15 ++- src/user/adjustment_initialization.F90 | 10 +- src/user/baroclinic_zone_initialization.F90 | 6 +- src/user/benchmark_initialization.F90 | 6 +- src/user/circle_obcs_initialization.F90 | 8 +- src/user/dense_water_initialization.F90 | 6 +- src/user/dumbbell_initialization.F90 | 16 ++- src/user/seamount_initialization.F90 | 10 +- src/user/sloshing_initialization.F90 | 6 +- src/user/soliton_initialization.F90 | 6 +- src/user/user_initialization.F90 | 2 +- 19 files changed, 212 insertions(+), 144 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 62b06a0b57..1907a75c74 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -149,6 +149,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & optional, intent(in) :: frac_shelf_h !< The fraction of the grid cell covered !! by a floating ice shelf [nondim]. ! Local variables + real :: depth_tot(SZI_(G),SZJ_(G)) ! The nominal total depth of the ocean [Z ~> m] character(len=200) :: filename ! The name of an input file. character(len=200) :: filename2 ! The name of an input files. character(len=200) :: inputdir ! The directory where NetCDF input files are. @@ -179,8 +180,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & logical :: debug_layers = .false. logical :: use_ice_shelf character(len=80) :: mesg -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB @@ -227,6 +228,13 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & !enddo endif + ! Set the nominal depth of the ocean, which might be different from the bathymetric + ! geopotential height, for use by the various initialization routines. G%bathyT has + ! already been initialized in previous calls. + do j=jsd,jed ; do i=isd,ied + depth_tot(i,j) = G%bathyT(i,j) + enddo ; enddo + ! The remaining initialization calls are done, regardless of whether the ! fields are actually initialized here (if just_read=.false.) or whether it ! is just to make sure that all valid parameters are read to enable the @@ -241,8 +249,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & if (.NOT.use_temperature) call MOM_error(FATAL,"MOM_initialize_state : "//& "use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true") - call MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_params=just_read,& - frac_shelf_h=frac_shelf_h) + call MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, & + just_read_params=just_read, frac_shelf_h=frac_shelf_h) else ! Initialize thickness, h. call get_param(PF, mdl, "THICKNESS_CONFIG", config, & @@ -275,9 +283,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & default="uniform", do_not_log=just_read) select case (trim(config)) case ("file") - call initialize_thickness_from_file(h, G, GV, US, PF, .false., just_read_params=just_read) + call initialize_thickness_from_file(h, depth_tot, G, GV, US, PF, .false., just_read_params=just_read) case ("thickness_file") - call initialize_thickness_from_file(h, G, GV, US, PF, .true., just_read_params=just_read) + call initialize_thickness_from_file(h, depth_tot, G, GV, US, PF, .true., just_read_params=just_read) case ("coord") if (new_sim .and. useALE) then call ALE_initThicknessToCoord( ALE_CSp, G, GV, h ) @@ -285,37 +293,37 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & call MOM_error(FATAL, "MOM_initialize_state: USE_REGRIDDING must be True "//& "for THICKNESS_CONFIG of 'coord'") endif - case ("uniform"); call initialize_thickness_uniform(h, G, GV, PF, & + case ("uniform"); call initialize_thickness_uniform(h, depth_tot, G, GV, PF, & just_read_params=just_read) - case ("list"); call initialize_thickness_list(h, G, GV, US, PF, & + case ("list"); call initialize_thickness_list(h, depth_tot, G, GV, US, PF, & just_read_params=just_read) - case ("DOME"); call DOME_initialize_thickness(h, G, GV, PF, & + case ("DOME"); call DOME_initialize_thickness(h, depth_tot, G, GV, PF, & just_read_params=just_read) - case ("ISOMIP"); call ISOMIP_initialize_thickness(h, G, GV, US, PF, tv, & + case ("ISOMIP"); call ISOMIP_initialize_thickness(h, depth_tot, G, GV, US, PF, tv, & just_read_params=just_read) - case ("benchmark"); call benchmark_initialize_thickness(h, G, GV, US, PF, & + case ("benchmark"); call benchmark_initialize_thickness(h, depth_tot, G, GV, US, PF, & tv%eqn_of_state, tv%P_Ref, just_read_params=just_read) - case ("Neverwoorld","Neverland"); call Neverworld_initialize_thickness(h, G, GV, US, PF, & + case ("Neverworld","Neverland"); call Neverworld_initialize_thickness(h, depth_tot, G, GV, US, PF, & tv%eqn_of_state, tv%P_Ref) - case ("search"); call initialize_thickness_search - case ("circle_obcs"); call circle_obcs_initialize_thickness(h, G, GV, PF, & + case ("search"); call initialize_thickness_search() + case ("circle_obcs"); call circle_obcs_initialize_thickness(h, depth_tot, G, GV, PF, & just_read_params=just_read) case ("lock_exchange"); call lock_exchange_initialize_thickness(h, G, GV, US, & PF, just_read_params=just_read) case ("external_gwave"); call external_gwave_initialize_thickness(h, G, GV, US, & PF, just_read_params=just_read) - case ("DOME2D"); call DOME2d_initialize_thickness(h, G, GV, US, PF, & + case ("DOME2D"); call DOME2d_initialize_thickness(h, depth_tot, G, GV, US, PF, & just_read_params=just_read) case ("adjustment2d"); call adjustment_initialize_thickness(h, G, GV, US, & PF, just_read_params=just_read) - case ("sloshing"); call sloshing_initialize_thickness(h, G, GV, US, PF, & + case ("sloshing"); call sloshing_initialize_thickness(h, depth_tot, G, GV, US, PF, & just_read_params=just_read) - case ("seamount"); call seamount_initialize_thickness(h, G, GV, US, PF, & + case ("seamount"); call seamount_initialize_thickness(h, depth_tot, G, GV, US, PF, & just_read_params=just_read) - case ("dumbbell"); call dumbbell_initialize_thickness(h, G, GV, US, PF, & + case ("dumbbell"); call dumbbell_initialize_thickness(h, depth_tot, G, GV, US, PF, & just_read_params=just_read) - case ("soliton"); call soliton_initialize_thickness(h, G, GV, US) - case ("phillips"); call Phillips_initialize_thickness(h, G, GV, US, PF, & + case ("soliton"); call soliton_initialize_thickness(h, depth_tot, G, GV, US) + case ("phillips"); call Phillips_initialize_thickness(h, depth_tot, G, GV, US, PF, & just_read_params=just_read) case ("rossby_front"); call Rossby_front_initialize_thickness(h, G, GV, US, & PF, just_read_params=just_read) @@ -363,11 +371,11 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & case ("DOME2D"); call DOME2d_initialize_temperature_salinity ( tv%T, & tv%S, h, G, GV, PF, eos, just_read_params=just_read) case ("ISOMIP"); call ISOMIP_initialize_temperature_salinity ( tv%T, & - tv%S, h, G, GV, US, PF, eos, just_read_params=just_read) + tv%S, h, depth_tot, G, GV, US, PF, eos, just_read_params=just_read) case ("adjustment2d"); call adjustment_initialize_temperature_salinity ( tv%T, & - tv%S, h, G, GV, PF, eos, just_read_params=just_read) + tv%S, h, depth_tot, G, GV, PF, eos, just_read_params=just_read) case ("baroclinic_zone"); call baroclinic_zone_init_temperature_salinity( tv%T, & - tv%S, h, G, GV, US, PF, just_read_params=just_read) + tv%S, h, depth_tot, G, GV, US, PF, just_read_params=just_read) case ("sloshing"); call sloshing_initialize_temperature_salinity(tv%T, & tv%S, h, G, GV, PF, eos, just_read_params=just_read) case ("seamount"); call seamount_initialize_temperature_salinity(tv%T, & @@ -547,22 +555,22 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & " \t\t for buoyancy-forced basin case.\n"//& " \t USER - call a user modified routine.", default="file") select case (trim(config)) - case ("DOME"); call DOME_initialize_sponges(G, GV, US, tv, PF, sponge_CSp) - case ("DOME2D"); call DOME2d_initialize_sponges(G, GV, US, tv, PF, useALE, & + case ("DOME"); call DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, sponge_CSp) + case ("DOME2D"); call DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, PF, useALE, & sponge_CSp, ALE_sponge_CSp) - case ("ISOMIP"); call ISOMIP_initialize_sponges(G, GV, US, tv, PF, useALE, & + case ("ISOMIP"); call ISOMIP_initialize_sponges(G, GV, US, tv, depth_tot, PF, useALE, & sponge_CSp, ALE_sponge_CSp) - case("RGC"); call RGC_initialize_sponges(G, GV, US, tv, u, v, PF, useALE, & + case("RGC"); call RGC_initialize_sponges(G, GV, US, tv, u, v, depth_tot, PF, useALE, & sponge_CSp, ALE_sponge_CSp) case ("USER"); call user_initialize_sponges(G, GV, use_temperature, tv, PF, sponge_CSp, h) - case ("BFB"); call BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, PF, & + case ("BFB"); call BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, depth_tot, PF, & sponge_CSp, h) - case ("DUMBBELL"); call dumbbell_initialize_sponges(G, GV, US, tv, PF, useALE, & + case ("DUMBBELL"); call dumbbell_initialize_sponges(G, GV, US, tv, depth_tot, PF, useALE, & sponge_CSp, ALE_sponge_CSp) case ("phillips"); call Phillips_initialize_sponges(G, GV, US, tv, PF, sponge_CSp, h) - case ("dense"); call dense_water_initialize_sponges(G, GV, US, tv, PF, useALE, & + case ("dense"); call dense_water_initialize_sponges(G, GV, US, tv, depth_tot, PF, useALE, & sponge_CSp, ALE_sponge_CSp) - case ("file"); call initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, PF, & + case ("file"); call initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_tot, PF, & sponge_CSp, ALE_sponge_CSp, Time) case default ; call MOM_error(FATAL, "MOM_initialize_state: "//& "Unrecognized sponge configuration "//trim(config)) @@ -639,14 +647,16 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & end subroutine MOM_initialize_state !> Reads the layer thicknesses or interface heights from a file. -subroutine initialize_thickness_from_file(h, G, GV, US, param_file, file_has_thickness, & +subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, file_has_thickness, & just_read_params) 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(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. - type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, intent(in) :: file_has_thickness !< If true, this file contains layer !! thicknesses; otherwise it contains @@ -696,7 +706,7 @@ subroutine initialize_thickness_from_file(h, G, GV, US, param_file, file_has_thi call MOM_read_data(filename, "eta", eta(:,:,:), G%Domain, scale=US%m_to_Z) if (correct_thickness) then - call adjustEtaToFitBathymetry(G, GV, US, eta, h, dZ_ref_eta=0.0) + call adjustEtaToFitBathymetry(G, GV, US, eta, h, dZ_ref_eta=0.0) else do k=nz,1,-1 ; do j=js,je ; do i=is,ie if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) then @@ -708,7 +718,7 @@ subroutine initialize_thickness_from_file(h, G, GV, US, param_file, file_has_thi enddo ; enddo ; enddo do j=js,je ; do i=is,ie - if (abs(eta(i,j,nz+1) + G%bathyT(i,j)) > 1.0*US%m_to_Z) & + if (abs(eta(i,j,nz+1) + depth_tot(i,j)) > 1.0*US%m_to_Z) & inconsistent = inconsistent + 1 enddo ; enddo call sum_across_PEs(inconsistent) @@ -811,11 +821,13 @@ subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h, dZ_ref_eta) end subroutine adjustEtaToFitBathymetry !> Initializes thickness to be uniform -subroutine initialize_thickness_uniform(h, G, GV, param_file, just_read_params) +subroutine initialize_thickness_uniform(h, depth_tot, G, GV, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, optional, intent(in) :: just_read_params !< If present and true, this call will @@ -850,7 +862,7 @@ subroutine initialize_thickness_uniform(h, G, GV, param_file, just_read_params) ! Angstrom thick, and 2. the interfaces are where they should be ! based on the resting depths and interface height perturbations, ! as long at this doesn't interfere with 1. - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(K) = e0(K) if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then @@ -866,12 +878,14 @@ subroutine initialize_thickness_uniform(h, G, GV, param_file, just_read_params) end subroutine initialize_thickness_uniform !> Initialize thickness from a 1D list -subroutine initialize_thickness_list(h, G, GV, US, param_file, just_read_params) +subroutine initialize_thickness_list(h, depth_tot, G, GV, US, param_file, just_read_params) 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(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, optional, intent(in) :: just_read_params !< If present and true, this call will @@ -926,7 +940,7 @@ subroutine initialize_thickness_list(h, G, GV, US, param_file, just_read_params) ! Angstrom thick, and 2. the interfaces are where they should be ! based on the resting depths and interface height perturbations, ! as long at this doesn't interfere with 1. - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(K) = e0(K) if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then @@ -1743,10 +1757,11 @@ end subroutine initialize_temp_salt_linear !! number of tracers should be restored within each sponge. The !! interface height is always subject to damping, and must always be !! the first registered field. -subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, param_file, Layer_CSp, ALE_CSp, Time) +subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_tot, param_file, & + Layer_CSp, ALE_CSp, Time) 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(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type logical, intent(in) :: use_temperature !< If true, T & S are state variables. type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic !! variables. @@ -1756,6 +1771,8 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, param_f real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & intent(in) :: v !< The meridional velocity that is being !! initialized [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. type(sponge_CS), pointer :: Layer_CSp !< A pointer that is set to point to the control !! structure for this module (in layered mode). @@ -1916,7 +1933,7 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, param_f call MOM_read_data(filename, eta_var, eta(:,:,:), G%Domain, scale=US%m_to_Z) do j=js,je ; do i=is,ie - eta(i,j,nz+1) = -G%bathyT(i,j) + eta(i,j,nz+1) = -depth_tot(i,j) enddo ; enddo do k=nz,1,-1 ; do j=js,je ; do i=is,ie if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) & @@ -1977,7 +1994,7 @@ subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, param_f allocate(h(isd:ied,jsd:jed,nz_data)) call MOM_read_data(filename, eta_var, eta(:,:,:), G%Domain, scale=US%m_to_Z) do j=js,je ; do i=is,ie - eta(i,j,nz+1) = -G%bathyT(i,j) + eta(i,j,nz+1) = -depth_tot(i,j) enddo ; enddo do k=nz,1,-1 ; do j=js,je ; do i=is,ie if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) & @@ -2256,13 +2273,15 @@ end subroutine set_velocity_depth_min !> This subroutine determines the isopycnal or other coordinate interfaces and !! layer potential temperatures and salinities directly from a z-space file on !! a latitude-longitude grid. -subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_params, frac_shelf_h) +subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just_read_params, frac_shelf_h) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(out) :: h !< Layer thicknesses being initialized [H ~> m or kg m-2] type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic !! variables including temperature and salinity + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: PF !< A structure indicating the open file !! to parse for model parameter values. @@ -2521,7 +2540,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param call pass_var(rho_z,G%Domain) do j=js,je ; do i=is,ie - Z_bottom(i,j) = -G%bathyT(i,j) + Z_bottom(i,j) = -depth_tot(i,j) enddo ; enddo ! Done with horizontal interpolation. diff --git a/src/user/BFB_initialization.F90 b/src/user/BFB_initialization.F90 index f632b95086..49c0a03235 100644 --- a/src/user/BFB_initialization.F90 +++ b/src/user/BFB_initialization.F90 @@ -76,13 +76,15 @@ end subroutine BFB_set_coord !> This subroutine sets up the sponges for the southern bouundary of the domain. Maximum damping occurs !! within 2 degrees lat of the boundary. The damping linearly decreases northward over the next 2 degrees. -subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, param_file, CSp, h) +subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, depth_tot, param_file, CSp, h) 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 logical, intent(in) :: use_temperature !< If true, temperature and salinity are used as !! state variables. type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters type(sponge_CS), pointer :: CSp !< A pointer to the sponge control structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & @@ -129,7 +131,7 @@ subroutine BFB_initialize_sponges_southonly(G, GV, US, use_temperature, tv, para max_damping = 1.0 / (86400.0*US%s_to_T) do j=js,je ; do i=is,ie - if (G%bathyT(i,j) <= min_depth) then ; Idamp(i,j) = 0.0 + if (depth_tot(i,j) <= min_depth) then ; Idamp(i,j) = 0.0 elseif (G%geoLatT(i,j) < slat+2.0) then ; Idamp(i,j) = max_damping elseif (G%geoLatT(i,j) < slat+4.0) then Idamp(i,j) = max_damping * (slat+4.0-G%geoLatT(i,j))/2.0 diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90 index 293d601757..f99f0b8d5c 100644 --- a/src/user/DOME2d_initialization.F90 +++ b/src/user/DOME2d_initialization.F90 @@ -90,12 +90,14 @@ subroutine DOME2d_initialize_topography( D, G, param_file, max_depth ) end subroutine DOME2d_initialize_topography !> Initialize thicknesses according to coordinate mode -subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_params ) +subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, just_read_params ) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, optional, intent(in) :: just_read_params !< If present and true, this call will @@ -150,7 +152,7 @@ subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_par case ( REGRIDDING_LAYER, REGRIDDING_RHO ) do j=js,je ; do i=is,ie - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(k) = e0(k) if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then @@ -172,7 +174,7 @@ subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_par ! case ( IC_RHO_C ) ! ! do j=js,je ; do i=is,ie - ! eta1D(nz+1) = -G%bathyT(i,j) + ! eta1D(nz+1) = -depth_tot(i,j) ! do k=nz,1,-1 ! eta1D(k) = e0(k) ! if (eta1D(k) < (eta1D(k+1) + min_thickness)) then @@ -194,7 +196,7 @@ subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_par case ( REGRIDDING_ZSTAR ) do j=js,je ; do i=is,ie - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(k) = e0(k) if (eta1D(k) < (eta1D(k+1) + min_thickness)) then @@ -208,7 +210,7 @@ subroutine DOME2d_initialize_thickness ( h, G, GV, US, param_file, just_read_par case ( REGRIDDING_SIGMA ) do j=js,je ; do i=is,ie - h(i,j,:) = GV%Z_to_H*G%bathyT(i,j) / nz + h(i,j,:) = GV%Z_to_H*depth_tot(i,j) / nz enddo ; enddo case default @@ -353,11 +355,13 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, param_file, end subroutine DOME2d_initialize_temperature_salinity !> Set up sponges in 2d DOME configuration -subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, ACSp) +subroutine DOME2d_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_ALE, CSp, ACSp) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< Parameter file structure logical, intent(in) :: use_ALE !< If true, indicates model is in ALE mode type(sponge_CS), pointer :: CSp !< Layer-mode sponge structure @@ -453,7 +457,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, AC enddo e0(nz+1) = -G%max_depth do j=js,je ; do i=is,ie - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(k) = e0(k) if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then @@ -470,7 +474,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, AC ! Construct temperature and salinity on the arbitrary grid T(:,:,:) = 0.0 ; S(:,:,:) = 0.0 do j=js,je ; do i=is,ie - z = -G%bathyT(i,j) + z = -depth_tot(i,j) do k = nz,1,-1 z = z + 0.5 * GV%H_to_Z * h(i,j,k) ! Position of the center of layer k S(i,j,k) = 34.0 - 1.0 * (z / (G%max_depth)) @@ -491,7 +495,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, AC ! Construct interface heights to restore toward do j=js,je ; do i=is,ie - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(K) = -G%max_depth * real(k-1) / real(nz) if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then @@ -508,7 +512,7 @@ subroutine DOME2d_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, AC d_eta(nz) = dome2d_depth_bay * G%max_depth - (nz-1) * GV%Angstrom_Z endif - eta(i,j,nz+1) = -G%bathyT(i,j) + eta(i,j,nz+1) = -depth_tot(i,j) do K=nz,1,-1 eta(i,j,K) = eta(i,j,K+1) + d_eta(k) enddo diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 81444704b3..1f3d24e1c9 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -87,11 +87,13 @@ end subroutine DOME_initialize_topography ! ----------------------------------------------------------------------------- !> This subroutine initializes layer thicknesses for the DOME experiment -subroutine DOME_initialize_thickness(h, G, GV, param_file, just_read_params) +subroutine DOME_initialize_thickness(h, depth_tot, G, GV, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, optional, intent(in) :: just_read_params !< If present and true, this call will @@ -124,7 +126,7 @@ subroutine DOME_initialize_thickness(h, G, GV, param_file, just_read_params) ! Angstrom thick, and 2. the interfaces are where they should be ! ! based on the resting depths and interface height perturbations, ! ! as long at this doesn't interfere with 1. ! - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(K) = e0(K) if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then @@ -145,17 +147,19 @@ end subroutine DOME_initialize_thickness !! number of tracers should be restored within each sponge. The ! !! interface height is always subject to damping, and must always be ! !! the first registered field. ! -subroutine DOME_initialize_sponges(G, GV, US, tv, PF, CSp) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. +subroutine DOME_initialize_sponges(G, GV, US, tv, depth_tot, PF, CSp) + 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 containing pointers to any available - !! thermodynamic fields, including potential temperature and - !! salinity or mixed layer density. Absent fields have NULL ptrs. - type(param_file_type), intent(in) :: PF !< A structure indicating the open file to - !! parse for model parameter values. - type(sponge_CS), pointer :: CSp !< A pointer that is set to point to the control - !! structure for this module. + type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available + !! thermodynamic fields, including potential temperature and + !! salinity or mixed layer density. Absent fields have NULL ptrs. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] + type(param_file_type), intent(in) :: PF !< A structure indicating the open file to + !! parse for model parameter values. + type(sponge_CS), pointer :: CSp !< A pointer that is set to point to the control + !! structure for this module. real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! A temporary array for interface heights [Z ~> m]. real :: temp(SZI_(G),SZJ_(G),SZK_(GV)) ! A temporary array for other variables. ! @@ -204,16 +208,16 @@ subroutine DOME_initialize_sponges(G, GV, US, tv, PF, CSp) ! depth space for Boussinesq or non-Boussinesq models. eta(i,j,1) = 0.0 do k=2,nz -! eta(i,j,K)=max(H0(k), -G%bathyT(i,j), GV%Angstrom_Z*(nz-k+1) - G%bathyT(i,j)) - e_dense = -G%bathyT(i,j) +! eta(i,j,K) = max(H0(k), -depth_tot(i,j), GV%Angstrom_Z*(nz-k+1) - depth_tot(i,j)) + e_dense = -depth_tot(i,j) if (e_dense >= H0(k)) then ; eta(i,j,K) = e_dense else ; eta(i,j,K) = H0(k) ; endif - if (eta(i,j,K) < GV%Angstrom_Z*(nz-k+1) - G%bathyT(i,j)) & - eta(i,j,K) = GV%Angstrom_Z*(nz-k+1) - G%bathyT(i,j) + if (eta(i,j,K) < GV%Angstrom_Z*(nz-k+1) - depth_tot(i,j)) & + eta(i,j,K) = GV%Angstrom_Z*(nz-k+1) - depth_tot(i,j) enddo - eta(i,j,nz+1) = -G%bathyT(i,j) + eta(i,j,nz+1) = -depth_tot(i,j) - if (G%bathyT(i,j) > min_depth) then + if (depth_tot(i,j) > min_depth) then Idamp(i,j) = damp / 86400.0 else ; Idamp(i,j) = 0.0 ; endif enddo ; enddo diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index aa1c6cdfe6..76f60d9b99 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -128,12 +128,14 @@ subroutine ISOMIP_initialize_topography(D, G, param_file, max_depth, US) end subroutine ISOMIP_initialize_topography !> Initialization of thicknesses -subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read_params) +subroutine ISOMIP_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv, just_read_params) 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(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any @@ -206,7 +208,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read ! Calculate thicknesses do j=js,je ; do i=is,ie - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(k) = e0(k) if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then @@ -221,7 +223,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR ) ! Initial thicknesses for z coordinates if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(k) = -G%max_depth * real(k-1) / real(nz) if (eta1D(k) < (eta1D(k+1) + min_thickness)) then @@ -236,7 +238,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie - h(i,j,:) = GV%Z_to_H * G%bathyT(i,j) / dfloat(nz) + h(i,j,:) = GV%Z_to_H * depth_tot(i,j) / dfloat(nz) enddo ; enddo case default @@ -248,7 +250,7 @@ subroutine ISOMIP_initialize_thickness ( h, G, GV, US, param_file, tv, just_read end subroutine ISOMIP_initialize_thickness !> Initial values for temperature and salinity -subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_file, & +subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, depth_tot, G, GV, US, param_file, & eqn_of_state, just_read_params) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure @@ -256,6 +258,8 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [ppt] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: depth_tot !< The nominal total bottom-to-top + !! depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< Parameter file structure type(EOS_type), pointer :: eqn_of_state !< Equation of state structure logical, optional, intent(in) :: just_read_params !< If present and true, this call will @@ -315,7 +319,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi dS_dz = (s_sur - s_bot) / G%max_depth dT_dz = (t_sur - t_bot) / G%max_depth do j=js,je ; do i=is,ie - xi0 = -G%bathyT(i,j) + xi0 = -depth_tot(i,j) do k = nz,1,-1 xi0 = xi0 + 0.5 * h(i,j,k) * GV%H_to_Z ! Depth in middle of layer S(i,j,k) = S_sur + dS_dz * xi0 @@ -420,7 +424,7 @@ end subroutine ISOMIP_initialize_temperature_salinity !> Sets up the the inverse restoration time (Idamp), and ! the values towards which the interface heights and an arbitrary ! number of tracers should be restored within each sponge. -subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp) +subroutine ISOMIP_initialize_sponges(G, GV, US, tv, depth_tot, PF, use_ALE, CSp, ACSp) 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 @@ -429,6 +433,8 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp) !! fields, potential temperature and !! salinity or mixed layer density. !! Absent fields have NULL ptrs. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: PF !< A structure indicating the !! open file to parse for model !! parameter values. @@ -508,7 +514,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp) ! and mask2dT is 1. do j=js,je ; do i=is,ie - if (G%bathyT(i,j) <= min_depth) then + if (depth_tot(i,j) <= min_depth) then Idamp(i,j) = 0.0 elseif (G%geoLonT(i,j) >= 790.0 .AND. G%geoLonT(i,j) <= 800.0) then dummy1 = (G%geoLonT(i,j)-790.0)/(800.0-790.0) @@ -549,7 +555,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp) ! Calculate thicknesses do j=js,je ; do i=is,ie - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(k) = e0(k) if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then @@ -563,7 +569,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp) case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR ) ! Initial thicknesses for z coordinates do j=js,je ; do i=is,ie - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(k) = -G%max_depth * real(k-1) / real(nz) if (eta1D(k) < (eta1D(k+1) + min_thickness)) then @@ -577,7 +583,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp) case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates do j=js,je ; do i=is,ie - h(i,j,:) = GV%Z_to_H * (G%bathyT(i,j) / dfloat(nz)) + h(i,j,:) = GV%Z_to_H * (depth_tot(i,j) / dfloat(nz)) enddo ; enddo case default @@ -593,7 +599,7 @@ subroutine ISOMIP_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp) dS_dz = (s_sur - s_bot) / G%max_depth dT_dz = (t_sur - t_bot) / G%max_depth do j=js,je ; do i=is,ie - xi0 = -G%bathyT(i,j) + xi0 = -depth_tot(i,j) do k = nz,1,-1 xi0 = xi0 + 0.5 * h(i,j,k) * GV%H_to_Z ! Depth in middle of layer S(i,j,k) = S_sur + dS_dz * xi0 diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90 index 4a136dd2db..ed944e5f0a 100644 --- a/src/user/Kelvin_initialization.F90 +++ b/src/user/Kelvin_initialization.F90 @@ -186,6 +186,7 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) real :: lambda ! Offshore decay scale [L-1 ~> m-1] real :: omega ! Wave frequency [T-1 ~> s-1] real :: PI + real :: depth_tot(SZI_(G),SZJ_(G)) ! The total depth of the ocean [Z ~> m] integer :: i, j, k, n, is, ie, js, je, isd, ied, jsd, jed, nz integer :: IsdB, IedB, JsdB, JedB real :: mag_SSH ! An overall magnitude of the external wave sea surface height at the coastline [Z ~> m] @@ -209,6 +210,17 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) PI = 4.0*atan(1.0) km_to_L_scale = 1000.0*US%m_to_L + do j=jsd,jed ; do i=isd,ied + depth_tot(i,j) = G%bathyT(i,j) + enddo ; enddo + !### Instead this should be: + ! do j=jsd,jed ; do i=isd,ied + ! depth_tot(i,j) = 0.0 + ! enddo ; enddo + ! do j=jsd,jed ; do i=isd,ied + ! depth_tot(i,j) = depth_tot(i,j) + GV%H_to_Z * h(i,j,k) + ! enddo ; enddo + if (CS%mode == 0) then mag_SSH = 1.0*US%m_to_Z omega = 2.0 * PI / (12.42 * 3600.0*US%s_to_T) ! M2 Tide period @@ -245,20 +257,17 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) y = -(x1 - CS%coast_offset1) * sina + y1 * cosa if (CS%mode == 0) then ! Use inside bathymetry - cff = sqrt(GV%g_Earth * G%bathyT(i+1,j) ) + cff = sqrt(GV%g_Earth * depth_tot(i+1,j) ) val2 = mag_SSH * exp(- CS%F_0 * y / cff) segment%eta(I,j) = GV%Z_to_H*val2 * cos(omega * time_sec) - segment%normal_vel_bt(I,j) = (val2 * (val1 * cff * cosa / & - (G%bathyT(i+1,j) )) ) + segment%normal_vel_bt(I,j) = val2 * (val1 * cff * cosa / depth_tot(i+1,j) ) if (segment%nudged) then do k=1,nz - segment%nudged_normal_vel(I,j,k) = (val2 * (val1 * cff * cosa / & - (G%bathyT(i+1,j))) ) + segment%nudged_normal_vel(I,j,k) = val2 * (val1 * cff * cosa / depth_tot(i+1,j) ) enddo elseif (segment%specified) then do k=1,nz - segment%normal_vel(I,j,k) = (val2 * (val1 * cff * cosa / & - (G%bathyT(i+1,j) )) ) + segment%normal_vel(I,j,k) = val2 * (val1 * cff * cosa / depth_tot(i+1,j) ) segment%normal_trans(I,j,k) = segment%normal_vel(I,j,k) * h(i+1,j,k) * G%dyCu(I,j) enddo endif @@ -288,11 +297,11 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) y1 = km_to_L_scale * G%geoLatBu(I,J) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = - (x1 - CS%coast_offset1) * sina + y1 * cosa - cff = sqrt(GV%g_Earth * G%bathyT(i+1,j) ) + cff = sqrt(GV%g_Earth * depth_tot(i+1,j) ) val2 = mag_SSH * exp(- CS%F_0 * y / cff) if (CS%mode == 0) then ; do k=1,nz segment%tangential_vel(I,J,k) = (val1 * val2 * cff * sina) / & - ( 0.5*(G%bathyT(i+1,j+1) + G%bathyT(i+1,j) ) ) + ( 0.5*(depth_tot(i+1,j+1) + depth_tot(i+1,j) ) ) enddo ; endif enddo ; enddo @@ -306,20 +315,17 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = - (x1 - CS%coast_offset1) * sina + y1 * cosa if (CS%mode == 0) then - cff = sqrt(GV%g_Earth * G%bathyT(i,j+1) ) + cff = sqrt(GV%g_Earth * depth_tot(i,j+1) ) val2 = mag_SSH * exp(- 0.5 * (G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J)) * y / cff) segment%eta(I,j) = GV%Z_to_H*val2 * cos(omega * time_sec) - segment%normal_vel_bt(I,j) = (val1 * cff * sina / & - (G%bathyT(i,j+1) )) * val2 + segment%normal_vel_bt(I,j) = (val1 * cff * sina / depth_tot(i,j+1) ) * val2 if (segment%nudged) then do k=1,nz - segment%nudged_normal_vel(I,j,k) = (val1 * cff * sina / & - (G%bathyT(i,j+1) )) * val2 + segment%nudged_normal_vel(I,j,k) = (val1 * cff * sina / depth_tot(i,j+1)) * val2 enddo elseif (segment%specified) then do k=1,nz - segment%normal_vel(I,j,k) = (val1 * cff * sina / & - (G%bathyT(i,j+1) )) * val2 + segment%normal_vel(I,j,k) = (val1 * cff * sina / depth_tot(i,j+1) ) * val2 segment%normal_trans(i,J,k) = segment%normal_vel(i,J,k) * h(i,j+1,k) * G%dxCv(i,J) enddo endif @@ -347,11 +353,11 @@ subroutine Kelvin_set_OBC_data(OBC, CS, G, GV, US, h, Time) y1 = km_to_L_scale * G%geoLatBu(I,J) x = (x1 - CS%coast_offset1) * cosa + y1 * sina y = - (x1 - CS%coast_offset1) * sina + y1 * cosa - cff = sqrt(GV%g_Earth * G%bathyT(i,j+1) ) + cff = sqrt(GV%g_Earth * depth_tot(i,j+1) ) val2 = mag_SSH * exp(- 0.5 * (G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J)) * y / cff) if (CS%mode == 0) then ; do k=1,nz - segment%tangential_vel(I,J,k) = ((val1 * val2 * cff * sina) / & - ( 0.5*((G%bathyT(i+1,j+1)) + G%bathyT(i,j+1))) ) + segment%tangential_vel(I,J,k) = (val1 * val2 * cff * sina) / & + ( 0.5*(depth_tot(i+1,j+1) + depth_tot(i,j+1)) ) enddo ; endif enddo ; enddo endif diff --git a/src/user/Neverworld_initialization.F90 b/src/user/Neverworld_initialization.F90 index 93a43e4a3e..3f5b8c8ab2 100644 --- a/src/user/Neverworld_initialization.F90 +++ b/src/user/Neverworld_initialization.F90 @@ -239,12 +239,14 @@ end function circ_ridge !! by finding the depths of interfaces in a specified latitude-dependent !! temperature profile with an exponentially decaying thermocline on top of a !! linear stratification. -subroutine Neverworld_initialize_thickness(h, G, GV, US, param_file, eqn_of_state, P_ref) +subroutine Neverworld_initialize_thickness(h, depth_tot, G, GV, US, param_file, eqn_of_state, P_ref) 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, intent(out), dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h !< The thickness that is being + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< The thickness that is being !! initialized [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open !! file to parse for model !! parameter values. @@ -283,7 +285,7 @@ subroutine Neverworld_initialize_thickness(h, G, GV, US, param_file, eqn_of_stat enddo do j=js,je ; do i=is,ie - e_interface = -G%bathyT(i,j) + e_interface = -depth_tot(i,j) do k=nz,2,-1 h(i,j,k) = GV%Z_to_H * (e0(k) - e_interface) ! Nominal thickness x=(G%geoLonT(i,j)-G%west_lon)/G%len_lon diff --git a/src/user/Phillips_initialization.F90 b/src/user/Phillips_initialization.F90 index dfa9c19460..448c86b5fb 100644 --- a/src/user/Phillips_initialization.F90 +++ b/src/user/Phillips_initialization.F90 @@ -35,12 +35,14 @@ module Phillips_initialization contains !> Initialize the thickness field for the Phillips model test case. -subroutine Phillips_initialize_thickness(h, G, GV, US, param_file, just_read_params) +subroutine Phillips_initialize_thickness(h, depth_tot, G, GV, US, param_file, just_read_params) 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(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, optional, intent(in) :: just_read_params !< If present and true, this call will @@ -113,7 +115,7 @@ subroutine Phillips_initialize_thickness(h, G, GV, US, param_file, just_read_par ! thicknesses are set to insure that: 1. each layer is at least an Angstrom thick, and ! 2. the interfaces are where they should be based on the resting depths and interface ! height perturbations, as long at this doesn't interfere with 1. - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(K) = eta_im(j,K) if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then diff --git a/src/user/RGC_initialization.F90 b/src/user/RGC_initialization.F90 index 4df728c22a..d051bccc6c 100644 --- a/src/user/RGC_initialization.F90 +++ b/src/user/RGC_initialization.F90 @@ -46,7 +46,7 @@ module RGC_initialization !> Sets up the the inverse restoration time, and the values towards which the interface heights, !! velocities and tracers should be restored within the sponges for the RGC test case. -subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp) +subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, depth_tot, PF, use_ALE, CSp, ACSp) 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 @@ -59,6 +59,8 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp) target, intent(in) :: u !< Array with the u velocity [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & target, intent(in) :: v !< Array with the v velocity [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: PF !< A structure indicating the !! open file to parse for model !! parameter values. @@ -93,10 +95,11 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed iscB = G%iscB ; iecB = G%iecB; jscB = G%jscB ; jecB = G%jecB - call get_param(PF,mod,"MIN_THICKNESS",min_thickness,'Minimum layer thickness',units='m',default=1.e-3) + call get_param(PF, mod,"MIN_THICKNESS", min_thickness, 'Minimum layer thickness', & + units='m', default=1.e-3) - call get_param(PF, mod, "RGC_TNUDG", TNUDG, 'Nudging time scale for sponge layers (days)', & - default=0.0, scale=86400.0*US%s_to_T) + call get_param(PF, mod, "RGC_TNUDG", TNUDG, 'Nudging time scale for sponge layers', & + units='days', default=0.0, scale=86400.0*US%s_to_T) call get_param(PF, mod, "LENLAT", lenlat, & "The latitudinal or y-direction length of the domain", & @@ -114,7 +117,7 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp) "Nudge velocities (u and v) towards zero in the sponge layer.", & default=.false., do_not_log=.true.) - T(:,:,:) = 0.0 ; S(:,:,:) = 0.0 ; Idamp(:,:) = 0.0; RHO(:,:,:) = 0.0 + T(:,:,:) = 0.0 ; S(:,:,:) = 0.0 ; Idamp(:,:) = 0.0 ; RHO(:,:,:) = 0.0 call get_param(PF, mod, "MINIMUM_DEPTH", min_depth, & "The minimum depth of the ocean.", units="m", default=0.0) @@ -130,7 +133,7 @@ subroutine RGC_initialize_sponges(G, GV, US, tv, u, v, PF, use_ALE, CSp, ACSp) ! and mask2dT is 1. do i=is,ie ; do j=js,je - if ((G%bathyT(i,j) <= min_depth) .or. (G%geoLonT(i,j) <= lensponge)) then + if ((depth_tot(i,j) <= min_depth) .or. (G%geoLonT(i,j) <= lensponge)) then Idamp(i,j) = 0.0 elseif (G%geoLonT(i,j) >= (lenlon - lensponge) .AND. G%geoLonT(i,j) <= lenlon) then dummy1 = (G%geoLonT(i,j)-(lenlon - lensponge))/(lensponge) diff --git a/src/user/adjustment_initialization.F90 b/src/user/adjustment_initialization.F90 index ad4eab33ff..b9f676dc55 100644 --- a/src/user/adjustment_initialization.F90 +++ b/src/user/adjustment_initialization.F90 @@ -57,8 +57,8 @@ subroutine adjustment_initialize_thickness ( h, G, GV, US, param_file, just_read real :: target_values(SZK_(GV)+1) ! Target densities or density anomalies [R ~> kg m-3] logical :: just_read ! If true, just read parameters but set nothing. character(len=20) :: verticalCoordinate -! This include declares and sets the variable "version". -#include "version_variable.h" + ! This include declares and sets the variable "version". +# include "version_variable.h" integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -193,13 +193,15 @@ subroutine adjustment_initialize_thickness ( h, G, GV, US, param_file, just_read end subroutine adjustment_initialize_thickness !> Initialization of temperature and salinity in the adjustment test case -subroutine adjustment_initialize_temperature_salinity(T, S, h, G, GV, param_file, & +subroutine adjustment_initialize_temperature_salinity(T, S, h, depth_tot, G, GV, param_file, & eqn_of_state, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< The temperature that is being initialized. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< The salinity that is being initialized. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< The model thicknesses [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to !! parse for model parameter values. type(EOS_type), pointer :: eqn_of_state !< Equation of state. @@ -260,7 +262,7 @@ subroutine adjustment_initialize_temperature_salinity(T, S, h, G, GV, param_file case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA ) dSdz = -delta_S_strat / G%max_depth do j=js,je ; do i=is,ie - eta1d(nz+1) = -G%bathyT(i,j) + eta1d(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1d(k) = eta1d(k+1) + h(i,j,k)*GV%H_to_Z enddo diff --git a/src/user/baroclinic_zone_initialization.F90 b/src/user/baroclinic_zone_initialization.F90 index b1c988e016..22f4d705a1 100644 --- a/src/user/baroclinic_zone_initialization.F90 +++ b/src/user/baroclinic_zone_initialization.F90 @@ -75,7 +75,7 @@ subroutine bcz_params(G, GV, US, param_file, S_ref, dSdz, delta_S, dSdx, T_ref, end subroutine bcz_params !> Initialization of temperature and salinity with the baroclinic zone initial conditions -subroutine baroclinic_zone_init_temperature_salinity(T, S, h, G, GV, US, param_file, & +subroutine baroclinic_zone_init_temperature_salinity(T, S, h, depth_tot, G, GV, US, param_file, & just_read_params) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -83,6 +83,8 @@ subroutine baroclinic_zone_init_temperature_salinity(T, S, h, G, GV, US, param_f real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [ppt] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< The model thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< Parameter file handle logical, optional, intent(in) :: just_read_params !< If present and true, this call will !! only read parameters without changing T & S. @@ -109,7 +111,7 @@ subroutine baroclinic_zone_init_temperature_salinity(T, S, h, G, GV, US, param_f PI = 4.*atan(1.) do j = G%jsc,G%jec ; do i = G%isc,G%iec - zi = -G%bathyT(i,j) + zi = -depth_tot(i,j) x = G%geoLonT(i,j) - (G%west_lon + 0.5*G%len_lon) ! Relative to center of domain xd = x / G%len_lon ! -1/2 < xd 1/2 y = G%geoLatT(i,j) - (G%south_lat + 0.5*G%len_lat) ! Relative to center of domain diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90 index ed0bbbf069..d077e0fa6f 100644 --- a/src/user/benchmark_initialization.F90 +++ b/src/user/benchmark_initialization.F90 @@ -82,13 +82,15 @@ end subroutine benchmark_initialize_topography !! by finding the depths of interfaces in a specified latitude-dependent !! temperature profile with an exponentially decaying thermocline on top of a !! linear stratification. -subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state, & +subroutine benchmark_initialize_thickness(h, depth_tot, G, GV, US, param_file, eqn_of_state, & P_Ref, just_read_params) 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(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. type(EOS_type), pointer :: eqn_of_state !< Equation of state structure @@ -181,7 +183,7 @@ subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state ! are set to insure that: 1. each layer is at least Gv%Angstrom_m thick, and ! 2. the interfaces are where they should be based on the resting depths and interface ! height perturbations, as long at this doesn't interfere with 1. - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,2,-1 T_int = 0.5*(T0(k) + T0(k-1)) diff --git a/src/user/circle_obcs_initialization.F90 b/src/user/circle_obcs_initialization.F90 index 4dd5a7c606..29fb6647b3 100644 --- a/src/user/circle_obcs_initialization.F90 +++ b/src/user/circle_obcs_initialization.F90 @@ -28,11 +28,13 @@ module circle_obcs_initialization contains !> This subroutine initializes layer thicknesses for the circle_obcs experiment. -subroutine circle_obcs_initialize_thickness(h, G, GV, param_file, just_read_params) +subroutine circle_obcs_initialize_thickness(h, depth_tot, G, GV, param_file, just_read_params) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, optional, intent(in) :: just_read_params !< If present and true, this call will @@ -79,8 +81,8 @@ subroutine circle_obcs_initialize_thickness(h, G, GV, param_file, just_read_para enddo ! Uniform thicknesses for base state - do j=js,je ; do i=is,ie ! - eta1D(nz+1) = -G%bathyT(i,j) + do j=js,je ; do i=is,ie + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(K) = e0(K) if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then diff --git a/src/user/dense_water_initialization.F90 b/src/user/dense_water_initialization.F90 index e8fe345bb0..c1eb4fa2e7 100644 --- a/src/user/dense_water_initialization.F90 +++ b/src/user/dense_water_initialization.F90 @@ -152,11 +152,13 @@ subroutine dense_water_initialize_TS(G, GV, param_file, eqn_of_state, T, S, h, j end subroutine dense_water_initialize_TS !> Initialize the restoring sponges for the dense water experiment -subroutine dense_water_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, ACSp) +subroutine dense_water_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_ALE, CSp, ACSp) type(ocean_grid_type), intent(in) :: G !< Horizontal grid control structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid control structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< Parameter file structure logical, intent(in) :: use_ALE !< ALE flag type(sponge_CS), pointer :: CSp !< Layered sponge control structure pointer @@ -234,7 +236,7 @@ subroutine dense_water_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CS do j = G%jsc,G%jec do i = G%isc,G%iec - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k = nz,1,-1 eta1D(k) = e0(k) diff --git a/src/user/dumbbell_initialization.F90 b/src/user/dumbbell_initialization.F90 index f7b647dd27..463fe018b0 100644 --- a/src/user/dumbbell_initialization.F90 +++ b/src/user/dumbbell_initialization.F90 @@ -90,12 +90,14 @@ subroutine dumbbell_initialize_topography( D, G, param_file, max_depth ) end subroutine dumbbell_initialize_topography !> Initializes the layer thicknesses to be uniform in the dumbbell test case -subroutine dumbbell_initialize_thickness ( h, G, GV, US, param_file, just_read_params) +subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file, just_read_params) 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(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, optional, intent(in) :: just_read_params !< If present and true, this call will @@ -169,7 +171,7 @@ subroutine dumbbell_initialize_thickness ( h, G, GV, US, param_file, just_read_p e0(K) = max(-G%max_depth, e0(K)) ! Bound by bottom enddo do j=js,je ; do i=is,ie - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(k) = e0(k) if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then @@ -184,7 +186,7 @@ subroutine dumbbell_initialize_thickness ( h, G, GV, US, param_file, just_read_p case ( REGRIDDING_ZSTAR ) ! Initial thicknesses for z coordinates if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(k) = -G%max_depth * real(k-1) / real(nz) if (eta1D(k) < (eta1D(k+1) + min_thickness)) then @@ -199,7 +201,7 @@ subroutine dumbbell_initialize_thickness ( h, G, GV, US, param_file, just_read_p case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie - h(i,j,:) = GV%Z_to_H * G%bathyT(i,j) / dfloat(nz) + h(i,j,:) = GV%Z_to_H * depth_tot(i,j) / dfloat(nz) enddo ; enddo end select @@ -284,11 +286,13 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, param_file end subroutine dumbbell_initialize_temperature_salinity !> Initialize the restoring sponges for the dumbbell test case -subroutine dumbbell_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, ACSp) +subroutine dumbbell_initialize_sponges(G, GV, US, tv, depth_tot, param_file, use_ALE, CSp, ACSp) type(ocean_grid_type), intent(in) :: G !< Horizontal grid control structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid control structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< Parameter file structure logical, intent(in) :: use_ALE !< ALE flag type(sponge_CS), pointer :: CSp !< Layered sponge control structure pointer @@ -354,7 +358,7 @@ subroutine dumbbell_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, if (use_ALE) then ! construct a uniform grid for the sponge do j=G%jsc,G%jec ; do i=G%isc,G%iec - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = depth_tot(i,j) do k=nz,1,-1 eta1D(k) = -G%max_depth * real(k-1) / real(nz) if (eta1D(k) < (eta1D(k+1) + min_thickness)) then diff --git a/src/user/seamount_initialization.F90 b/src/user/seamount_initialization.F90 index 9118133108..6bfaedc221 100644 --- a/src/user/seamount_initialization.F90 +++ b/src/user/seamount_initialization.F90 @@ -77,12 +77,14 @@ end subroutine seamount_initialize_topography !> Initialization of thicknesses. !! This subroutine initializes the layer thicknesses to be uniform. -subroutine seamount_initialize_thickness ( h, G, GV, US, param_file, just_read_params) +subroutine seamount_initialize_thickness (h, depth_tot, G, GV, US, param_file, just_read_params) 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(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, optional, intent(in) :: just_read_params !< If present and true, this call will @@ -152,7 +154,7 @@ subroutine seamount_initialize_thickness ( h, G, GV, US, param_file, just_read_p e0(K) = max(-G%max_depth, e0(K)) ! Bound by bottom enddo do j=js,je ; do i=is,ie - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(k) = e0(k) if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then @@ -167,7 +169,7 @@ subroutine seamount_initialize_thickness ( h, G, GV, US, param_file, just_read_p case ( REGRIDDING_ZSTAR ) ! Initial thicknesses for z coordinates if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie - eta1D(nz+1) = -G%bathyT(i,j) + eta1D(nz+1) = -depth_tot(i,j) do k=nz,1,-1 eta1D(k) = -G%max_depth * real(k-1) / real(nz) if (eta1D(k) < (eta1D(k+1) + min_thickness)) then @@ -182,7 +184,7 @@ subroutine seamount_initialize_thickness ( h, G, GV, US, param_file, just_read_p case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie - h(i,j,:) = GV%Z_to_H * G%bathyT(i,j) / dfloat(nz) + h(i,j,:) = GV%Z_to_H * depth_tot(i,j) / dfloat(nz) enddo ; enddo end select diff --git a/src/user/sloshing_initialization.F90 b/src/user/sloshing_initialization.F90 index e1c0a96d63..1c3334d8b0 100644 --- a/src/user/sloshing_initialization.F90 +++ b/src/user/sloshing_initialization.F90 @@ -53,12 +53,14 @@ end subroutine sloshing_initialize_topography !! same thickness but all interfaces (except bottom and sea surface) are !! displaced according to a half-period cosine, with maximum value on the !! left and minimum value on the right. This sets off a regular sloshing motion. -subroutine sloshing_initialize_thickness ( h, G, GV, US, param_file, just_read_params) +subroutine sloshing_initialize_thickness ( h, depth_tot, G, GV, US, param_file, just_read_params) 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(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, optional, intent(in) :: just_read_params !< If present and true, this call will @@ -152,7 +154,7 @@ subroutine sloshing_initialize_thickness ( h, G, GV, US, param_file, just_read_p enddo ! 3. The last interface must coincide with the seabed - z_inter(nz+1) = -G%bathyT(i,j) + z_inter(nz+1) = -depth_tot(i,j) ! Modify interface heights to make sure all thicknesses are strictly positive do k = nz,1,-1 if ( z_inter(k) < (z_inter(k+1) + GV%Angstrom_Z) ) then diff --git a/src/user/soliton_initialization.F90 b/src/user/soliton_initialization.F90 index ac6ec8c4bc..f62aa54f88 100644 --- a/src/user/soliton_initialization.F90 +++ b/src/user/soliton_initialization.F90 @@ -28,12 +28,14 @@ module soliton_initialization contains !> Initialization of thicknesses in Equatorial Rossby soliton test -subroutine soliton_initialize_thickness(h, G, GV, US) +subroutine soliton_initialize_thickness(h, depth_tot, G, GV, US) 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(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] integer :: i, j, k, is, ie, js, je, nz real :: x, y, x0, y0 @@ -55,7 +57,7 @@ subroutine soliton_initialize_thickness(h, G, GV, US) y = G%geoLatT(i,j)-y0 val3 = exp(-val1*x) val4 = val2 * ( 2.0*val3 / (1.0 + (val3*val3)) )**2 - h(i,j,k) = GV%Z_to_H * (0.25*val4*(6.0*y*y + 3.0) * exp(-0.5*y*y) + G%bathyT(i,j)) + h(i,j,k) = GV%Z_to_H * (0.25*val4*(6.0*y*y + 3.0) * exp(-0.5*y*y) + depth_tot(i,j)) enddo enddo ; enddo diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90 index 793b87f149..3338121d9e 100644 --- a/src/user/user_initialization.F90 +++ b/src/user/user_initialization.F90 @@ -247,7 +247,7 @@ end subroutine write_user_log !! - u - Zonal velocity [Z T-1 ~> m s-1]. !! - v - Meridional velocity [Z T-1 ~> m s-1]. !! - h - Layer thickness [H ~> m or kg m-2]. (Must be positive.) -!! - G%bathyT - Basin depth [Z ~> m]. (Must be positive.) +!! - G%bathyT - Basin depth [Z ~> m]. !! - G%CoriolisBu - The Coriolis parameter [T-1 ~> s-1]. !! - GV%g_prime - The reduced gravity at each interface [L2 Z-1 T-2 ~> m s-2]. !! - GV%Rlay - Layer potential density (coordinate variable) [R ~> kg m-3]. From 8a0ae944d181ec6d109a029568f27eb23d55e46f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 17 Aug 2021 07:19:29 -0400 Subject: [PATCH 11/11] +Pass depth_tot around within MOM_MEKE Added a new argument, depth_tot, to three routines within MOM_MEKE.F90 to replace places where G%bathyT had been used. This new depth_tot variable is currently set to G%bathyT, but there is commented-out code suggesting a more appropriate expression based on the temporally evolving total thickness of the water column. All answers are bitwise identical, but there are new arguments to 3 private routines. --- src/parameterizations/lateral/MOM_MEKE.F90 | 81 +++++++++++++--------- 1 file changed, 48 insertions(+), 33 deletions(-) diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 9441ab7107..8779968bcd 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -126,6 +126,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h real, dimension(SZI_(G),SZJ_(G)) :: & mass, & ! The total mass of the water column [R Z ~> kg m-2]. I_mass, & ! The inverse of mass [R-1 Z-1 ~> m2 kg-1]. + depth_tot, & ! The depth of the water column [Z ~> m]. src, & ! The sum of all MEKE sources [L2 T-3 ~> W kg-1] (= m2 s-3). MEKE_decay, & ! A diagnostic of the MEKE decay timescale [T-1 ~> s-1]. drag_rate_visc, & ! Near-bottom velocity contribution to bottom dratg [L T-1 ~> m s-1] @@ -161,7 +162,8 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h real :: advFac ! The product of the advection scaling factor and 1/dt [T-1 ~> s-1] real :: mass_neglect ! A negligible mass [R Z ~> kg m-2]. real :: ldamping ! The MEKE damping rate [T-1 ~> s-1]. - real :: Rho0 ! A density used to convert mass to distance [R ~> kg m-3]. + real :: Rho0 ! A density used to convert mass to distance [R ~> kg m-3] + real :: I_Rho0 ! The inverse of the density used to convert mass to distance [R-1 ~> m3 kg-1] real :: sdt ! dt to use locally [T ~> s] (could be scaled to accelerate) real :: sdt_damp ! dt for damping [T ~> s] (sdt could be split). logical :: use_drag_rate ! Flag to indicate drag_rate is finite @@ -204,6 +206,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h sdt = dt*CS%MEKE_dtScale ! Scaled dt to use for time-stepping Rho0 = GV%Rho0 + I_Rho0 = 1.0 / GV%Rho0 mass_neglect = GV%H_to_RZ * GV%H_subroundoff cdrag2 = CS%cdrag**2 @@ -280,13 +283,20 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo enddo + !$OMP parallel do default(shared) + do j=js-1,je+1 ; do i=is-1,ie+1 + depth_tot(i,j) = G%bathyT(i,j) + !### Try changing this to: + ! depth_tot(i,j) = mass(i,j) * I_Rho0 + enddo ; enddo + if (CS%initialize) then - call MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass) + call MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass, depth_tot) CS%initialize = .false. endif ! Calculates bottomFac2, barotrFac2 and LmixScale - call MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, MEKE%MEKE, bottomFac2, barotrFac2, LmixScale) + call MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, MEKE%MEKE, depth_tot, bottomFac2, barotrFac2, LmixScale) if (CS%debug) then if (CS%visc_drag) & call uvchksum("MEKE drag_vel_[uv]", drag_vel_u, drag_vel_v, G%HI, & @@ -323,7 +333,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h !$OMP parallel do default(shared) do j=js,je ; do i=is,ie src(i,j) = src(i,j) - CS%MEKE_GMcoeff*MEKE%GM_src(i,j) / & - (GV%Rho0 * MAX(1.0*US%m_to_Z, G%bathyT(i,j))) + (GV%Rho0 * MAX(1.0*US%m_to_Z, depth_tot(i,j))) enddo ; enddo else !$OMP parallel do default(shared) @@ -334,7 +344,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif if (CS%MEKE_equilibrium_restoring) then - call MEKE_equilibrium_restoring(CS, G, US, SN_u, SN_v) + call MEKE_equilibrium_restoring(CS, G, US, SN_u, SN_v, depth_tot) do j=js,je ; do i=is,ie src(i,j) = src(i,j) - CS%MEKE_restoring_rate*(MEKE%MEKE(i,j) - CS%equilibrium_value(i,j)) enddo ; enddo @@ -640,7 +650,7 @@ end subroutine step_forward_MEKE !> Calculates the equilibrium solution where the source depends only on MEKE diffusivity !! and there is no lateral diffusion of MEKE. !! Results is in MEKE%MEKE. -subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass) +subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_mass, depth_tot) type(ocean_grid_type), intent(inout) :: G !< Ocean grid. type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -651,6 +661,8 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m real, dimension(SZI_(G),SZJ_(G)), intent(in) :: drag_rate_visc !< Mean flow velocity contribution !! to the MEKE drag rate [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJ_(G)), intent(in) :: I_mass !< Inverse of column mass [R-1 Z-1 ~> m2 kg-1]. + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: depth_tot !< The depth of the water column [Z ~> m]. + ! Local variables real :: beta ! Combined topograpic and planetary vorticity gradient [T-1 L-1 ~> s-1 m-1] real :: SN ! The local Eady growth rate [T-1 ~> s-1] @@ -690,27 +702,27 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1)) if (CS%MEKE_equilibrium_alt) then - MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*G%bathyT(i,j))**2 / cd2 + MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*depth_tot(i,j))**2 / cd2 else FatH = 0.25*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))) ! Coriolis parameter at h points ! Since zero-bathymetry cells are masked, this avoids calculations on land - if (CS%MEKE_topographic_beta == 0. .or. G%bathyT(i,j) == 0.) then + if (CS%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then beta_topo_x = 0. ; beta_topo_y = 0. else !### Consider different combinations of these estimates of topographic beta, and the use ! of the water column thickness instead of the bathymetric depth. beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & - (G%bathyT(i+1,j)-G%bathyT(i,j)) * G%IdxCu(I,j) & - / max(G%bathyT(i+1,j),G%bathyT(i,j), dZ_neglect) & - + (G%bathyT(i,j)-G%bathyT(i-1,j)) * G%IdxCu(I-1,j) & - / max(G%bathyT(i,j),G%bathyT(i-1,j), dZ_neglect) ) + (depth_tot(i+1,j)-depth_tot(i,j)) * G%IdxCu(I,j) & + / max(depth_tot(i+1,j), depth_tot(i,j), dZ_neglect) & + + (depth_tot(i,j)-depth_tot(i-1,j)) * G%IdxCu(I-1,j) & + / max(depth_tot(i,j), depth_tot(i-1,j), dZ_neglect) ) beta_topo_y = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & - (G%bathyT(i,j+1)-G%bathyT(i,j)) * G%IdyCv(i,J) & - / max(G%bathyT(i,j+1),G%bathyT(i,j), dZ_neglect) + & - (G%bathyT(i,j)-G%bathyT(i,j-1)) * G%IdyCv(i,J-1) & - / max(G%bathyT(i,j),G%bathyT(i,j-1), dZ_neglect) ) + (depth_tot(i,j+1)-depth_tot(i,j)) * G%IdyCv(i,J) & + / max(depth_tot(i,j+1), depth_tot(i,j), dZ_neglect) + & + (depth_tot(i,j)-depth_tot(i,j-1)) * G%IdyCv(i,J-1) & + / max(depth_tot(i,j), depth_tot(i,j-1), dZ_neglect) ) endif beta = sqrt((G%dF_dx(i,j) + beta_topo_x)**2 + & (G%dF_dy(i,j) + beta_topo_y)**2 ) @@ -729,7 +741,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m do while (resid>0.) n1 = n1 + 1 EKE = EKEmax - call MEKE_lengthScales_0d(CS, US, G%areaT(i,j), beta, G%bathyT(i,j), & + call MEKE_lengthScales_0d(CS, US, G%areaT(i,j), beta, depth_tot(i,j), & MEKE%Rd_dx_h(i,j), SN, EKE, & bottomFac2, barotrFac2, LmixScale, LRhines, LEady) ! TODO: Should include resolution function in Kh @@ -804,12 +816,14 @@ end subroutine MEKE_equilibrium !< This subroutine calculates a new equilibrium value for MEKE at each time step. This is not copied into !! MEKE%MEKE; rather, it is used as a restoring term to nudge MEKE%MEKE back to an equilibrium value -subroutine MEKE_equilibrium_restoring(CS, G, US, SN_u, SN_v) +subroutine MEKE_equilibrium_restoring(CS, G, US, SN_u, SN_v, depth_tot) type(ocean_grid_type), intent(inout) :: G !< Ocean grid. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type. type(MEKE_CS), pointer :: CS !< MEKE control structure. real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1]. real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1]. + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: depth_tot !< The depth of the water column [Z ~> m]. + ! Local variables real :: SN ! The local Eady growth rate [T-1 ~> s-1] integer :: i, j, is, ie, js, je ! local indices @@ -826,7 +840,7 @@ subroutine MEKE_equilibrium_restoring(CS, G, US, SN_u, SN_v) ! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.) ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1)) - CS%equilibrium_value(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*G%bathyT(i,j))**2 / cd2 + CS%equilibrium_value(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_L*depth_tot(i,j))**2 / cd2 enddo ; enddo if (CS%id_MEKE_equilibrium>0) call post_data(CS%id_MEKE_equilibrium, CS%equilibrium_value, CS%diag) @@ -836,8 +850,8 @@ end subroutine MEKE_equilibrium_restoring !> Calculates the eddy mixing length scale and \f$\gamma_b\f$ and \f$\gamma_t\f$ !! functions that are ratios of either bottom or barotropic eddy energy to the !! column eddy energy, respectively. See \ref section_MEKE_equations. -subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, & - EKE, bottomFac2, barotrFac2, LmixScale) +subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, EKE, depth_tot, & + bottomFac2, barotrFac2, LmixScale) type(MEKE_CS), pointer :: CS !< MEKE control structure. type(MEKE_type), pointer :: MEKE !< MEKE data. type(ocean_grid_type), intent(inout) :: G !< Ocean grid. @@ -846,8 +860,9 @@ subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, & real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: SN_u !< Eady growth rate at u-points [T-1 ~> s-1]. real, dimension(SZI_(G),SZJB_(G)), intent(in) :: SN_v !< Eady growth rate at v-points [T-1 ~> s-1]. real, dimension(SZI_(G),SZJ_(G)), intent(in) :: EKE !< Eddy kinetic energy [L2 T-2 ~> m2 s-2]. - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: bottomFac2 !< gamma_b^2 - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: barotrFac2 !< gamma_t^2 + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: depth_tot !< The depth of the water column [Z ~> m]. + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: bottomFac2 !< gamma_b^2 [nondim] + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: barotrFac2 !< gamma_t^2 [nondim] real, dimension(SZI_(G),SZJ_(G)), intent(out) :: LmixScale !< Eddy mixing length [L ~> m]. ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: LRhines, LEady ! Possible mixing length scales [L ~> m] @@ -875,21 +890,21 @@ subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, & ! If bathyT is zero, then a division by zero FPE will be raised. In this ! case, we apply Adcroft's rule of reciprocals and set the term to zero. ! Since zero-bathymetry cells are masked, this should not affect values. - if (CS%MEKE_topographic_beta == 0. .or. G%bathyT(i,j) == 0.0) then + if (CS%MEKE_topographic_beta == 0. .or. (depth_tot(i,j) == 0.0)) then beta_topo_x = 0. ; beta_topo_y = 0. else !### Consider different combinations of these estimates of topographic beta, and the use ! of the water column thickness instead of the bathymetric depth. beta_topo_x = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & - (G%bathyT(i+1,j)-G%bathyT(i,j)) * G%IdxCu(I,j) & - / max(G%bathyT(i+1,j),G%bathyT(i,j), dZ_neglect) & - + (G%bathyT(i,j)-G%bathyT(i-1,j)) * G%IdxCu(I-1,j) & - / max(G%bathyT(i,j),G%bathyT(i-1,j), dZ_neglect) ) + (depth_tot(i+1,j)-depth_tot(i,j)) * G%IdxCu(I,j) & + / max(depth_tot(i+1,j), depth_tot(i,j), dZ_neglect) & + + (depth_tot(i,j)-depth_tot(i-1,j)) * G%IdxCu(I-1,j) & + / max(depth_tot(i,j), depth_tot(i-1,j), dZ_neglect) ) beta_topo_y = -CS%MEKE_topographic_beta * FatH * 0.5 * ( & - (G%bathyT(i,j+1)-G%bathyT(i,j)) * G%IdyCv(i,J) & - / max(G%bathyT(i,j+1),G%bathyT(i,j), dZ_neglect) + & - (G%bathyT(i,j)-G%bathyT(i,j-1)) * G%IdyCv(i,J-1) & - / max(G%bathyT(i,j),G%bathyT(i,j-1), dZ_neglect) ) + (depth_tot(i,j+1)-depth_tot(i,j)) * G%IdyCv(i,J) & + / max(depth_tot(i,j+1), depth_tot(i,j), dZ_neglect) + & + (depth_tot(i,j)-depth_tot(i,j-1)) * G%IdyCv(i,J-1) & + / max(depth_tot(i,j), depth_tot(i,j-1), dZ_neglect) ) endif beta = sqrt((G%dF_dx(i,j) + beta_topo_x)**2 + & (G%dF_dy(i,j) + beta_topo_y)**2 ) @@ -898,7 +913,7 @@ subroutine MEKE_lengthScales(CS, MEKE, G, GV, US, SN_u, SN_v, & beta = 0. endif ! Returns bottomFac2, barotrFac2 and LmixScale - call MEKE_lengthScales_0d(CS, US, G%areaT(i,j), beta, G%bathyT(i,j), & + call MEKE_lengthScales_0d(CS, US, G%areaT(i,j), beta, depth_tot(i,j), & MEKE%Rd_dx_h(i,j), SN, MEKE%MEKE(i,j), & bottomFac2(i,j), barotrFac2(i,j), LmixScale(i,j), & LRhines(i,j), LEady(i,j))