Merge pull request #61 from NOAA-GFDL/dev/gfdl

Sync with NOAA-GFDL dev/gfdl

.testing/tc2/MOM_input

@@ -297,6 +297,10 @@ BOUND_CORIOLIS = True           !   [Boolean] default = False
                                 ! v-points, and similarly at v-points.  This option would
                                 ! have no effect on the SADOURNY Coriolis scheme if it
                                 ! were possible to use centered difference thickness fluxes.
+PGF_STANLEY_T2_DET_COEFF = 0.5  !   [nondim] default = -1.0
+                                ! The coefficient correlating SGS temperature variance with the mean temperature
+                                ! gradient in the deterministic part of the Stanley form of the Brankart
+                                ! correction. Negative values disable the scheme.
 
 ! === module MOM_hor_visc ===
 LAPLACIAN = True
@@ -426,6 +430,10 @@ KHTH = 1.0                      !   [m2 s-1] default = 0.0
                                 ! The background horizontal thickness diffusivity.
 KHTH_MAX = 900.0                !   [m2 s-1] default = 0.0
                                 ! The maximum horizontal thickness diffusivity.
+STANLEY_PRM_DET_COEFF = 0.5     !   [nondim] default = -1.0
+                                ! The coefficient correlating SGS temperature variance with the mean temperature
+                                ! gradient in the deterministic part of the Stanley parameterization. Negative
+                                ! values disable the scheme.
 
 ! === module MOM_mixed_layer_restrat ===
 FOX_KEMPER_ML_RESTRAT_COEF = 5.0 !   [nondim] default = 0.0

src/ALE/MOM_ALE.F90

@@ -52,6 +52,7 @@ module MOM_ALE
 use regrid_consts,        only : coordinateUnits, coordinateMode, state_dependent
 use regrid_edge_values,   only : edge_values_implicit_h4
 use PLM_functions,        only : PLM_reconstruction, PLM_boundary_extrapolation
+use PLM_functions,        only : PLM_extrapolate_slope, PLM_monotonized_slope, PLM_slope_wa
 use PPM_functions,        only : PPM_reconstruction, PPM_boundary_extrapolation
 
 implicit none ; private
@@ -110,8 +111,9 @@ module MOM_ALE
 public ALE_build_grid
 public ALE_regrid_accelerated
 public ALE_remap_scalar
-public pressure_gradient_plm
-public pressure_gradient_ppm
+public ALE_PLM_edge_values
+public TS_PLM_edge_values
+public TS_PPM_edge_values
 public adjustGridForIntegrity
 public ALE_initRegridding
 public ALE_getCoordinate
@@ -1006,12 +1008,9 @@ subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_c
 end subroutine ALE_remap_scalar
 
 
-!> Use plm reconstruction for pressure gradient (determine edge values)
-!! By using a PLM (limited piecewise linear method) reconstruction, this
-!! routine determines the edge values for the salinity and temperature
-!! within each layer. These edge values are returned and are used to compute
-!! the pressure gradient (by computing the densities).
-subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )
+!> Calculate edge values (top and bottom of layer) for T and S consistent with a PLM reconstruction
+!! in the vertical direction. Boundary reconstructions are PCM unless bdry_extrap is true.
+subroutine TS_PLM_edge_values( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )
   type(ocean_grid_type),   intent(in)    :: G    !< ocean grid structure
   type(verticalGrid_type), intent(in)    :: GV   !< Ocean vertical grid structure
   type(ALE_CS),            intent(inout) :: CS   !< module control structure
@@ -1029,12 +1028,31 @@ subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_ext
   logical,                 intent(in)    :: bdry_extrap !< If true, use high-order boundary
                                                  !! extrapolation within boundary cells
 
+  call ALE_PLM_edge_values( CS, G, GV, h, tv%S, bdry_extrap, S_t, S_b )
+  call ALE_PLM_edge_values( CS, G, GV, h, tv%T, bdry_extrap, T_t, T_b )
+
+end subroutine TS_PLM_edge_values
+
+!> Calculate edge values (top and bottom of layer) 3d scalar array.
+!! Boundary reconstructions are PCM unless bdry_extrap is true.
+subroutine ALE_PLM_edge_values( CS, G, GV, h, Q, bdry_extrap, Q_t, Q_b )
+  type(ALE_CS),            intent(in)    :: CS   !< module control structure
+  type(ocean_grid_type),   intent(in)    :: G    !< ocean grid structure
+  type(verticalGrid_type), intent(in)    :: GV   !< Ocean vertical grid structure
+  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),SZK_(GV)), &
+                           intent(in)    :: Q    !< 3d scalar array
+  logical,                 intent(in)    :: bdry_extrap !< If true, use high-order boundary
+                                                 !! extrapolation within boundary cells
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                           intent(inout) :: Q_t  !< Scalar at the top edge of each layer
+  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
+                           intent(inout) :: Q_b  !< Scalar at the bottom edge of each layer
   ! Local variables
   integer :: i, j, k
-  real    :: hTmp(GV%ke)
-  real    :: tmp(GV%ke)
-  real, dimension(CS%nk,2) :: ppol_E     !Edge value of polynomial
-  real, dimension(CS%nk,2) :: ppol_coefs !Coefficients of polynomial
+  real :: slp(GV%ke)
+  real :: mslp
   real :: h_neglect
 
   if (.not.CS%answers_2018) then
@@ -1045,48 +1063,40 @@ subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_ext
     h_neglect = GV%kg_m2_to_H*1.0e-30
   endif
 
-  ! Determine reconstruction within each column
-  !$OMP parallel do default(shared) private(hTmp,ppol_E,ppol_coefs,tmp)
+  !$OMP parallel do default(shared) private(slp,mslp)
   do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
-    ! Build current grid
-    hTmp(:) = h(i,j,:)
-    tmp(:) = tv%S(i,j,:)
-    ! Reconstruct salinity profile
-    ppol_E(:,:) = 0.0
-    ppol_coefs(:,:) = 0.0
-    call PLM_reconstruction( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )
-    if (bdry_extrap) &
-      call PLM_boundary_extrapolation( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )
-
-    do k = 1,GV%ke
-      S_t(i,j,k) = ppol_E(k,1)
-      S_b(i,j,k) = ppol_E(k,2)
+    slp(1) = 0.
+    do k = 2, GV%ke-1
+      slp(k) = PLM_slope_wa(h(i,j,k-1), h(i,j,k), h(i,j,k+1), h_neglect, Q(i,j,k-1), Q(i,j,k), Q(i,j,k+1))
     enddo
+    slp(GV%ke) = 0.
 
-    ! Reconstruct temperature profile
-    ppol_E(:,:) = 0.0
-    ppol_coefs(:,:) = 0.0
-    tmp(:) = tv%T(i,j,:)
-    call PLM_reconstruction( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )
-    if (bdry_extrap) &
-      call PLM_boundary_extrapolation( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )
-
-    do k = 1,GV%ke
-      T_t(i,j,k) = ppol_E(k,1)
-      T_b(i,j,k) = ppol_E(k,2)
+    do k = 2, GV%ke-1
+      mslp = PLM_monotonized_slope(Q(i,j,k-1), Q(i,j,k), Q(i,j,k+1), slp(k-1), slp(k), slp(k+1))
+      Q_t(i,j,k) = Q(i,j,k) - 0.5 * mslp
+      Q_b(i,j,k) = Q(i,j,k) + 0.5 * mslp
     enddo
+    if (bdry_extrap) then
+      mslp = - PLM_extrapolate_slope(h(i,j,2), h(i,j,1), h_neglect, Q(i,j,2), Q(i,j,1))
+      Q_t(i,j,1) = Q(i,j,1) - 0.5 * mslp
+      Q_b(i,j,1) = Q(i,j,1) + 0.5 * mslp
+      mslp = PLM_extrapolate_slope(h(i,j,GV%ke-1), h(i,j,GV%ke), h_neglect, Q(i,j,GV%ke-1), Q(i,j,GV%ke))
+      Q_t(i,j,GV%ke) = Q(i,j,GV%ke) - 0.5 * mslp
+      Q_b(i,j,GV%ke) = Q(i,j,GV%ke) + 0.5 * mslp
+    else
+      Q_t(i,j,1) = Q(i,j,1)
+      Q_b(i,j,1) = Q(i,j,1)
+      Q_t(i,j,GV%ke) = Q(i,j,GV%ke)
+      Q_b(i,j,GV%ke) = Q(i,j,GV%ke)
+    endif
 
   enddo ; enddo
 
-end subroutine pressure_gradient_plm
-
+end subroutine ALE_PLM_edge_values
 
-!> Use ppm reconstruction for pressure gradient (determine edge values)
-!> By using a PPM (limited piecewise linear method) reconstruction, this
-!> routine determines the edge values for the salinity and temperature
-!> within each layer. These edge values are returned and are used to compute
-!> the pressure gradient (by computing the densities).
-subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )
+!> Calculate edge values (top and bottom of layer) for T and S consistent with a PPM reconstruction
+!! in the vertical direction. Boundary reconstructions are PCM unless bdry_extrap is true.
+subroutine TS_PPM_edge_values( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )
   type(ocean_grid_type),   intent(in)    :: G    !< ocean grid structure
   type(verticalGrid_type), intent(in)    :: GV   !< Ocean vertical grid structure
   type(ALE_CS),            intent(inout) :: CS   !< module control structure
@@ -1168,7 +1178,7 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_ext
 
   enddo ; enddo
 
-end subroutine pressure_gradient_ppm
+end subroutine TS_PPM_edge_values
 
 
 !> Initializes regridding for the main ALE algorithm

src/ALE/MOM_remapping.F90

@@ -1027,11 +1027,11 @@ real function average_value_ppoly( n0, u0, ppoly0_E, ppoly0_coefs, method, i0, x
 end function average_value_ppoly
 
 !> Measure totals and bounds on source grid
-subroutine measure_input_bounds( n0, h0, u0, ppoly_E, h0tot, h0err, u0tot, u0err, u0min, u0max )
+subroutine measure_input_bounds( n0, h0, u0, edge_values, h0tot, h0err, u0tot, u0err, u0min, u0max )
   integer,               intent(in)  :: n0 !< Number of cells on source grid
   real, dimension(n0),   intent(in)  :: h0 !< Cell widths on source grid
   real, dimension(n0),   intent(in)  :: u0 !< Cell averages on source grid
-  real, dimension(n0,2), intent(in)  :: ppoly_E !< Cell edge values on source grid
+  real, dimension(n0,2), intent(in)  :: edge_values !< Cell edge values on source grid
   real,                  intent(out) :: h0tot !< Sum of cell widths
   real,                  intent(out) :: h0err !< Magnitude of round-off error in h0tot
   real,                  intent(out) :: u0tot !< Sum of cell widths times values
@@ -1047,15 +1047,15 @@ subroutine measure_input_bounds( n0, h0, u0, ppoly_E, h0tot, h0err, u0tot, u0err
   h0err = 0.
   u0tot = h0(1) * u0(1)
   u0err = 0.
-  u0min = min( ppoly_E(1,1), ppoly_E(1,2) )
-  u0max = max( ppoly_E(1,1), ppoly_E(1,2) )
+  u0min = min( edge_values(1,1), edge_values(1,2) )
+  u0max = max( edge_values(1,1), edge_values(1,2) )
   do k = 2, n0
     h0tot = h0tot + h0(k)
     h0err = h0err + eps * max(h0tot, h0(k))
     u0tot = u0tot + h0(k) * u0(k)
     u0err = u0err + eps * max(abs(u0tot), abs(h0(k) * u0(k)))
-    u0min = min( u0min, ppoly_E(k,1), ppoly_E(k,2) )
-    u0max = max( u0max, ppoly_E(k,1), ppoly_E(k,2) )
+    u0min = min( u0min, edge_values(k,1), edge_values(k,2) )
+    u0max = max( u0max, edge_values(k,1), edge_values(k,2) )
   enddo
 
 end subroutine measure_input_bounds

src/ALE/P1M_functions.F90

@@ -24,11 +24,11 @@ module P1M_functions
 !!
 !! It is assumed that the size of the array 'u' is equal to the number of cells
 !! defining 'grid' and 'ppoly'. No consistency check is performed here.
-subroutine P1M_interpolation( N, h, u, ppoly_E, ppoly_coef, h_neglect, answers_2018 )
+subroutine P1M_interpolation( N, h, u, edge_values, ppoly_coef, h_neglect, answers_2018 )
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(:),   intent(in)    :: h !< cell widths (size N) [H]
   real, dimension(:),   intent(in)    :: u !< cell average properties (size N) [A]
-  real, dimension(:,:), intent(inout) :: ppoly_E !< Potentially modified edge values [A]
+  real, dimension(:,:), intent(inout) :: edge_values !< Potentially modified edge values [A]
   real, dimension(:,:), intent(inout) :: ppoly_coef !< Potentially modified
                                            !! piecewise polynomial coefficients, mainly [A]
   real,       optional, intent(in)    :: h_neglect !< A negligibly small width [H]
@@ -39,17 +39,17 @@ subroutine P1M_interpolation( N, h, u, ppoly_E, ppoly_coef, h_neglect, answers_2
   real      :: u0_l, u0_r   ! edge values (left and right)
 
   ! Bound edge values (routine found in 'edge_values.F90')
-  call bound_edge_values( N, h, u, ppoly_E, h_neglect, answers_2018 )
+  call bound_edge_values( N, h, u, edge_values, h_neglect, answers_2018 )
 
   ! Systematically average discontinuous edge values (routine found in
   ! 'edge_values.F90')
-  call average_discontinuous_edge_values( N, ppoly_E )
+  call average_discontinuous_edge_values( N, edge_values )
 
   ! Loop on interior cells to build interpolants
   do k = 1,N
 
-    u0_l = ppoly_E(k,1)
-    u0_r = ppoly_E(k,2)
+    u0_l = edge_values(k,1)
+    u0_r = edge_values(k,2)
 
     ppoly_coef(k,1) = u0_l
     ppoly_coef(k,2) = u0_r - u0_l
@@ -65,12 +65,12 @@ end subroutine P1M_interpolation
 !!
 !! It is assumed that the size of the array 'u' is equal to the number of cells
 !! defining 'grid' and 'ppoly'. No consistency check is performed here.
-subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef )
+subroutine P1M_boundary_extrapolation( N, h, u, edge_values, ppoly_coef )
   ! Arguments
   integer,              intent(in)    :: N !< Number of cells
   real, dimension(:),   intent(in)    :: h !< cell widths (size N) [H]
   real, dimension(:),   intent(in)    :: u !< cell averages (size N) [A]
-  real, dimension(:,:), intent(inout) :: ppoly_E !< edge values of piecewise polynomials [A]
+  real, dimension(:,:), intent(inout) :: edge_values !< edge values of piecewise polynomials [A]
   real, dimension(:,:), intent(inout) :: ppoly_coef !< coefficients of piecewise polynomials, mainly [A]
 
   ! Local variables
@@ -99,20 +99,20 @@ subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef )
   ! by using the edge value in the neighboring cell.
   u0_r = u0 + 0.5 * slope
 
-  if ( (u1 - u0) * (ppoly_E(2,1) - u0_r) < 0.0 ) then
-    slope = 2.0 * ( ppoly_E(2,1) - u0 )
+  if ( (u1 - u0) * (edge_values(2,1) - u0_r) < 0.0 ) then
+    slope = 2.0 * ( edge_values(2,1) - u0 )
   endif
 
   ! Using the limited slope, the left edge value is reevaluated and
   ! the interpolant coefficients recomputed
   if ( h0 /= 0.0 ) then
-    ppoly_E(1,1) = u0 - 0.5 * slope
+    edge_values(1,1) = u0 - 0.5 * slope
   else
-    ppoly_E(1,1) = u0
+    edge_values(1,1) = u0
   endif
 
-  ppoly_coef(1,1) = ppoly_E(1,1)
-  ppoly_coef(1,2) = ppoly_E(1,2) - ppoly_E(1,1)
+  ppoly_coef(1,1) = edge_values(1,1)
+  ppoly_coef(1,2) = edge_values(1,2) - edge_values(1,1)
 
   ! ------------------------------------------
   ! Right edge value in the left boundary cell
@@ -127,18 +127,18 @@ subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coef )
 
   u0_l = u1 - 0.5 * slope
 
-  if ( (u1 - u0) * (u0_l - ppoly_E(N-1,2)) < 0.0 ) then
-    slope = 2.0 * ( u1 - ppoly_E(N-1,2) )
+  if ( (u1 - u0) * (u0_l - edge_values(N-1,2)) < 0.0 ) then
+    slope = 2.0 * ( u1 - edge_values(N-1,2) )
   endif
 
   if ( h1 /= 0.0 ) then
-    ppoly_E(N,2) = u1 + 0.5 * slope
+    edge_values(N,2) = u1 + 0.5 * slope
   else
-    ppoly_E(N,2) = u1
+    edge_values(N,2) = u1
   endif
 
-  ppoly_coef(N,1) = ppoly_E(N,1)
-  ppoly_coef(N,2) = ppoly_E(N,2) - ppoly_E(N,1)
+  ppoly_coef(N,1) = edge_values(N,1)
+  ppoly_coef(N,2) = edge_values(N,2) - edge_values(N,1)
 
 end subroutine P1M_boundary_extrapolation