From 3675c216975876b24d31353f560043d2eea58052 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 6 Mar 2022 08:58:20 -0500 Subject: [PATCH] +Add Hybgen remapping options This commit adds the ability to specify three remapping options derived from Hycom's hybgen code. - Adds the new module MOM_hybgen_remap, which contains the new subroutines hybgen_plm_coefs, hybgen_ppm_coefs, and hybgen_weno_coefs - Adds code to handle PLM_HYBGEN, PPM_HYBGEN and WENO_HYBGEN as valid entries for specifying the ALE remapping schemes. Also added descriptions of units to some internal remapping variables. - Adds the optional argument PCM_cell to regridding_main, remapping_core_h, and remap_all_state_vars to specify layers that should use piecewise constant remapping, regardless of the overall remapping scheme, to follow the approach used in Hycom. ALE_main uses these new optional arguments, although until the Hybgen regridding code is added, they will always be set to false. - Makes 7 character strings longer in 5 files to accommodate the new remapping options. All answers are bitwise identical, but there are changes to some entries in the MOM_parameter_doc files. --- src/ALE/MOM_ALE.F90 | 19 +- src/ALE/MOM_hybgen_remap.F90 | 391 ++++++++++++++++++ src/ALE/MOM_regridding.F90 | 10 +- src/ALE/MOM_remapping.F90 | 112 +++-- src/core/MOM_open_boundary.F90 | 4 +- .../MOM_state_initialization.F90 | 2 +- .../MOM_tracer_initialization_from_Z.F90 | 2 +- src/ocean_data_assim/MOM_oda_incupd.F90 | 2 +- .../vertical/MOM_ALE_sponge.F90 | 19 +- .../vertical/MOM_CVMix_KPP.F90 | 4 +- 10 files changed, 515 insertions(+), 50 deletions(-) create mode 100644 src/ALE/MOM_hybgen_remap.F90 diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 46669c20cb..d8eee55c14 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -380,6 +380,7 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h) real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: eta_preale real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_new ! New 3D grid obtained after last time step [H ~> m or kg m-2] + logical :: PCM_cell(SZI_(G),SZJ_(G),SZK_(GV)) !< If true, PCM remapping should be used in a cell. integer :: nk, i, j, k, isc, iec, jsc, jec, ntr nk = GV%ke; isc = G%isc; iec = G%iec; jsc = G%jsc; jec = G%jec @@ -405,7 +406,7 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h) ! Build new grid. The new grid is stored in h_new. The old grid is h. ! Both are needed for the subsequent remapping of variables. call regridding_main( CS%remapCS, CS%regridCS, G, GV, h, tv, h_new, dzRegrid, & - frac_shelf_h ) + frac_shelf_h, PCM_cell=PCM_cell) call check_grid( G, GV, h, 0. ) @@ -419,7 +420,7 @@ subroutine ALE_main( G, GV, US, h, u, v, tv, Reg, CS, OBC, dt, frac_shelf_h) ! Remap all variables from old grid h onto new grid h_new call remap_all_state_vars( CS, G, GV, h, h_new, Reg, OBC, dzRegrid, u, v, & - CS%show_call_tree, dt ) + CS%show_call_tree, dt, PCM_cell=PCM_cell ) if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_main)") @@ -776,7 +777,7 @@ end subroutine ALE_regrid_accelerated !! remap initial conditions to the model grid. It is also called during a !! time step to update the state. subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & - dzInterface, u, v, debug, dt ) + dzInterface, u, v, debug, dt, PCM_cell) type(ALE_CS), intent(in) :: CS !< ALE control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure @@ -795,6 +796,8 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & optional, intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1] logical, optional, intent(in) :: debug !< If true, show the call tree real, optional, intent(in) :: dt !< time step for diagnostics [T ~> s] + logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + optional, intent(in) :: PCM_cell !< Use PCM remapping in cells where true ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: h_tot ! The vertically summed thicknesses [H ~> m or kg m-2] @@ -861,8 +864,14 @@ subroutine remap_all_state_vars(CS, G, GV, h_old, h_new, Reg, OBC, & ! Build the start and final grids h1(:) = h_old(i,j,:) h2(:) = h_new(i,j,:) - call remapping_core_h(CS%remapCS, nz, h1, Tr%t(i,j,:), nz, h2, tr_column, & - h_neglect, h_neglect_edge) + if (present(PCM_cell)) then + PCM(:) = PCM_cell(i,j,:) + call remapping_core_h(CS%remapCS, nz, h1, Tr%t(i,j,:), nz, h2, tr_column, & + h_neglect, h_neglect_edge, PCM_cell=PCM) + else + call remapping_core_h(CS%remapCS, nz, h1, Tr%t(i,j,:), nz, h2, tr_column, & + h_neglect, h_neglect_edge) + endif ! Intermediate steps for tendency of tracer concentration and tracer content. if (present(dt)) then diff --git a/src/ALE/MOM_hybgen_remap.F90 b/src/ALE/MOM_hybgen_remap.F90 new file mode 100644 index 0000000000..539c992283 --- /dev/null +++ b/src/ALE/MOM_hybgen_remap.F90 @@ -0,0 +1,391 @@ +!> This module contains the hybgen remapping routines from HYCOM, with minor +!! modifications to follow the MOM6 coding conventions +module MOM_hybgen_remap + +! This file is part of MOM6. See LICENSE.md for the license. + +implicit none ; private + +public hybgen_plm_coefs, hybgen_ppm_coefs, hybgen_weno_coefs + +contains + +!> Set up the coefficients for PLM remapping of a set of scalars +subroutine hybgen_plm_coefs(si, dpi, slope, nk, ns, thin, PCM_lay) + integer, intent(in) :: nk !< The number of input layers + integer, intent(in) :: ns !< The number of scalar fields to work on + real, intent(in) :: si(nk,ns) !< The cell-averaged input scalar fields [A] + real, intent(in) :: dpi(nk) !< The input grid layer thicknesses [H ~> m or kg m-2] + real, intent(out) :: slope(nk,ns) !< The PLM slope times cell width [A] + real, intent(in) :: thin !< A negligible layer thickness that can be ignored [H ~> m or kg m-2] + logical, optional, intent(in) :: PCM_lay(nk) !< If true for a layer, use PCM remapping for that layer + +!----------------------------------------------------------------------- +! 1) coefficients for remapping from one set of vertical cells to another. +! method: piecewise linear across each input cell with +! monotonized central-difference limiter. +! +! van Leer, B., 1977, J. Comp. Phys., 23 276-299. +! +! 2) input arguments: +! si - initial scalar fields in pi-layer space +! dpi - initial layer thicknesses (dpi(k) = pi(k+1)-pi(k)) +! nk - number of layers +! ns - number of fields +! thin - layer thickness (>0) that can be ignored +! PCM_lay - use PCM for selected layers (optional) +! +! 3) output arguments: +! slope - coefficients for hybgen_plm_remap +! profile(y) = si+slope*(y-1), -0.5 <= y <= 0.5 +! +! 4) Tim Campbell, Mississippi State University, October 2002. +! Alan J. Wallcraft, Naval Research Laboratory, Aug. 2007. +!----------------------------------------------------------------------- +! + real :: qcen ! A layer's thickness divided by the distance between the centers + ! of the adjacent cells, usually ~0.5, but always <= 1 [nondim] + real :: zbot, zcen, ztop ! Tracer slopes times the layer thickness [A] + integer :: i, k + + do i=1,ns + slope(1, i) = 0.0 + slope(nk,i) = 0.0 + enddo !i + do k= 2,nk-1 + if (dpi(k) <= thin) then !use PCM + do i=1,ns ; slope(k,i) = 0.0 ; enddo + else +! --- use qcen in place of 0.5 to allow for non-uniform grid + qcen = dpi(k) / (dpi(k)+0.5*(dpi(k-1)+dpi(k+1))) !dpi(k)>thin + do i=1,ns +! --- PLM (non-zero slope, but no new extrema) +! --- layer value is si-0.5*slope at top interface, +! --- and si+0.5*slope at bottom interface. +! +! --- monotonized central-difference limiter (van Leer, 1977, +! --- JCP 23 pp 276-299). For a discussion of PLM limiters, see +! --- Finite Volume Methods for Hyperbolic Problems by R.J. Leveque. + ztop = 2.0*(si(k, i)-si(k-1,i)) + zbot = 2.0*(si(k+1,i)-si(k, i)) + zcen = qcen*(si(k+1,i)-si(k-1,i)) + if (ztop*zbot > 0.0) then !ztop,zbot are the same sign + slope(k,i) = sign(min(abs(zcen),abs(zbot),abs(ztop)), zbot) + else + slope(k,i) = 0.0 !local extrema, so no slope + endif + enddo !i + endif !PCM:PLM + enddo !k + + if (present(PCM_lay)) then + do k=1,nk ; if (PCM_lay(k)) then + do i=1,ns ; slope(k,i) = 0.0 ; enddo + endif ; enddo + endif + +end subroutine hybgen_plm_coefs + + +!> Set up the coefficients for PPM remapping of a set of scalars +subroutine hybgen_ppm_coefs(s, h_src, edges, nk, ns, thin, PCM_lay) + integer, intent(in) :: nk !< The number of input layers + integer, intent(in) :: ns !< The scalar fields to work on + real, intent(in) :: s(nk,ns) !< The input scalar fields [A] + real, intent(in) :: h_src(nk) !< The input grid layer thicknesses [H ~> m or kg m-2] + real, intent(out) :: edges(nk,2,ns) !< The PPM interpolation edge values of the scalar fields [A] + real, intent(in) :: thin !< A negligible layer thickness that can be ignored [H ~> m or kg m-2] + logical, optional, intent(in) :: PCM_lay(nk) !< If true for a layer, use PCM remapping for that layer + +!----------------------------------------------------------------------- +! 1) coefficients for remapping from one set of vertical cells to another. +! method: monotonic piecewise parabolic across each input cell +! +! Colella, P. & P.R. Woodward, 1984, J. Comp. Phys., 54, 174-201. +! +! 2) input arguments: +! s - initial scalar fields in pi-layer space +! h_src - initial layer thicknesses (>=0) +! nk - number of layers +! ns - number of fields +! thin - layer thickness (>0) that can be ignored +! PCM_lay - use PCM for selected layers (optional) +! +! 3) output arguments: +! edges - cell edge scalar values for the PPM reconstruction +! edges.1 is value at interface above +! edges.2 is value at interface below +! +! 4) Tim Campbell, Mississippi State University, October 2002. +! Alan J. Wallcraft, Naval Research Laboratory, Aug. 2007. +!----------------------------------------------------------------------- +! + real :: dp(nk) ! Input grid layer thicknesses, but with a minimum thickness given by thin [H ~> m or kg m-2] + logical :: PCM_layer(nk) ! True for layers that should use PCM remapping, either because they are + ! very thin, or because this is specified by PCM_lay. + real :: da ! Difference between the unlimited scalar edge value estimates [A] + real :: a6 ! Scalar field differences that are proportional to the curvature [A] + real :: slk, srk ! Differences between adjacent cell averages of scalars [A] + real :: sck ! Scalar differences across a cell. + real :: as(nk) ! Scalar field difference across each cell [A] + real :: al(nk), ar(nk) ! Scalar field at the left and right edges of a cell [A] + real :: h112(nk+1), h122(nk+1) ! Combinations of thicknesses [H ~> m or kg m-2] + real :: I_h12(nk+1) ! Inverses of combinations of thickesses [H-1 ~> m-1 or m2 kg-1] + real :: h2_h123(nk) ! A ratio of a layer thickness of the sum of 3 adjacent thicknesses [nondim] + real :: I_h0123(nk) ! Inverse of the sum of 4 adjacent thicknesses [H-1 ~> m-1 or m2 kg-1] + real :: h01_h112(nk+1) ! A ratio of sums of adjacent thicknesses [nondim], 2/3 in the limit of uniform thicknesses. + real :: h23_h122(nk+1) ! A ratio of sums of adjacent thicknesses [nondim], 2/3 in the limit of uniform thicknesses. + integer :: k, i + + ! This PPM remapper is not currently written to work with massless layers, so set + ! the thicknesses for very thin layers to some minimum value. + do k=1,nk ; dp(k) = max(h_src(k), thin) ; enddo + + ! Specify the layers that will use PCM remapping. + if (present(PCM_lay)) then + do k=1,nk ; PCM_layer(k) = (PCM_lay(k) .or. dp(k) <= thin) ; enddo + else + do k=1,nk ; PCM_layer(k) = (dp(k) <= thin) ; enddo + endif + + !compute grid metrics + do k=2,nk + h112(K) = 2.*dp(k-1) + dp(k) + h122(K) = dp(k-1) + 2.*dp(k) + I_h12(K) = 1.0 / (dp(k-1) + dp(k)) + enddo !k + do k=2,nk-1 + h2_h123(k) = dp(k) / (dp(k) + (dp(k-1)+dp(k+1))) + enddo + do K=3,nk-1 + I_h0123(K) = 1.0 / ((dp(k-2) + dp(k-1)) + (dp(k) + dp(k+1))) + + h01_h112(K) = (dp(k-2) + dp(k-1)) / (2.0*dp(k-1) + dp(k)) + h23_h122(K) = (dp(k) + dp(k+1)) / (dp(k-1) + 2.0*dp(k)) + enddo + + do i=1,ns + !Compute average slopes: Colella, Eq. (1.8) + as(1) = 0. + do k=2,nk-1 + if (PCM_layer(k)) then !use PCM + as(k) = 0.0 + else + slk = s(k, i)-s(k-1,i) + srk = s(k+1,i)-s(k, i) + if (slk*srk > 0.) then + sck = h2_h123(k)*( h112(K)*srk*I_h12(K+1) + h122(K+1)*slk*I_h12(K) ) + as(k) = sign(min(abs(2.0*slk), abs(sck), abs(2.0*srk)), sck) + else + as(k) = 0. + endif + endif !PCM:PPM + enddo !k + as(nk) = 0. + !Compute "first guess" edge values: Colella, Eq. (1.6) + al(1) = s(1,i) ! 1st layer PCM + ar(1) = s(1,i) ! 1st layer PCM + al(2) = s(1,i) ! 1st layer PCM + do K=3,nk-1 + ! This is a 4th order explicit edge value estimate. + al(k) = (dp(k)*s(k-1,i) + dp(k-1)*s(k,i)) * I_h12(K) & + + I_h0123(K)*( 2.*dp(k)*dp(k-1)*I_h12(K)*(s(k,i)-s(k-1,i)) * & + ( h01_h112(K) - h23_h122(K) ) & + + (dp(k)*as(k-1)*h23_h122(K) - dp(k-1)*as(k)*h01_h112(K)) ) + ar(k-1) = al(k) + enddo !k + ar(nk-1) = s(nk,i) ! last layer PCM + al(nk) = s(nk,i) ! last layer PCM + ar(nk) = s(nk,i) ! last layer PCM + !Impose monotonicity: Colella, Eq. (1.10) + do k=2,nk-1 + if ((PCM_layer(k)) .or. ((s(k+1,i)-s(k,i))*(s(k,i)-s(k-1,i)) <= 0.)) then !local extremum + al(k) = s(k,i) + ar(k) = s(k,i) + else + da = ar(k)-al(k) + a6 = 6.0*s(k,i) - 3.0*(al(k)+ar(k)) + if (da*a6 > da*da) then !peak in right half of zone + al(k) = 3.0*s(k,i) - 2.0*ar(k) + elseif (da*a6 < -da*da) then !peak in left half of zone + ar(k) = 3.0*s(k,i) - 2.0*al(k) + endif + endif + enddo !k + !Set coefficients + do k=1,nk + edges(k,1,i) = al(k) + edges(k,2,i) = ar(k) + enddo !k + enddo !i + +end subroutine hybgen_ppm_coefs + + +!> Set up the coefficients for PPM remapping of a set of scalars +subroutine hybgen_weno_coefs(s, h_src, edges, nk, ns, thin, PCM_lay) + integer, intent(in) :: nk !< The number of input layers + integer, intent(in) :: ns !< The number of scalar fields to work on + real, intent(in) :: s(nk,ns) !< The input scalar fields [A] + real, intent(in) :: h_src(nk) !< The input grid layer thicknesses [H ~> m or kg m-2] + real, intent(out) :: edges(nk,2,ns) !< The WENO interpolation edge values of the scalar fields [A] + real, intent(in) :: thin !< A negligible layer thickness that can be ignored [H ~> m or kg m-2] + logical, optional, intent(in) :: PCM_lay(nk) !< If true for a layer, use PCM remapping for that layer + +!----------------------------------------------------------------------- +! 1) coefficients for remapping from one set of vertical cells to another. +! method: monotonic WENO-like alternative to PPM across each input cell +! a second order polynomial approximation of the profiles +! using a WENO reconciliation of the slopes to compute the +! interfacial values +! +! This scheme might have ben developed by Shchepetkin. A.F., personal communication. +! See also Engwirda, D., and M. Kelley, A WENO-type slope-limiter for a family of piecewise +! polynomial methods, arXive:1606.08188v1, 27 June 2016. +! +! 2) input arguments: +! s - initial scalar fields in pi-layer space +! h_src - initial layer thicknesses (>=0) +! nk - number of layers +! ns - number of fields +! thin - layer thickness (>0) that can be ignored +! PCM_lay - use PCM for selected layers (optional) +! +! 3) output arguments: +! edges - cell edge scalar values for the WENO reconstruction +! edges.1 is value at interface above +! edges.2 is value at interface below +! +! 4) Laurent Debreu, Grenoble. +! Alan J. Wallcraft, Naval Research Laboratory, July 2008. +!----------------------------------------------------------------------- +! +! real, parameter :: dsmll=1.0e-8 ! This has units of [A2], and hence can not be a parameter. +! + real :: curv_cell ! An estimate of the tracer curvature centered on a cell times the grid + ! spacing [A H-1 ~> A m-1 or A kg m-2] + real :: seh1, seh2 ! Tracer slopes at the cell edges times the cell grid spacing [A] + real :: q01, q02 ! Various tracer differences between a cell average and the edge values [A] + real :: q001, q002 ! Tracer slopes at the cell edges times the cell grid spacing [A] + real :: ds2a, ds2b ! Squared tracer differences between a cell average and the edge values [A2] + logical :: PCM_layer(nk) ! True for layers that should use PCM remapping, either because they are + ! very thin, or because this is specified by PCM_lay. + real :: dp(nk) ! Input grid layer thicknesses, but with a minimum thickness given by thin [H ~> m or kg m-2] + real :: qdpkm(nk) ! Inverse of the sum of two adjacent thicknesses [H-1 ~> m-1 or m2 kg-1] + real :: qdpkmkp(nk) ! Inverse of the sum of three adjacent thicknesses [H-1 ~> m-1 or m2 kg-1] + real :: dpkm2kp(nk) ! Twice the distance between the centers of the layers two apart [H ~> m or kg m-2] + real :: zw(nk,2) ! Squared combinations of the differences between the the cell average tracer + ! concentrations and the left and right edges [A2] + real :: min_ratio ! The minimum ratio of the values of zw used to interpolate the edge values [nondim] + real :: wt1 ! The weight of the upper layer in the interpolated shared edge value [nondim] + real :: slope_edge(nk+1) ! Tracer slopes at the edges [A H-1 ~> A m-1 or A kg m-2] + real :: val_edge(nk+1) ! A weighted average edge concentration [A] + integer :: i, k + + min_ratio = 1.0e-8 + + ! The WENO remapper is not currently written to work with massless layers, so set + ! the thicknesses for very thin layers to some minimum value. + do k=1,nk ; dp(k) = max(h_src(k), thin) ; enddo + + ! Specify the layers that will use PCM remapping. + if (present(PCM_lay)) then + do k=1,nk ; PCM_layer(k) = (PCM_lay(k) .or. dp(k) <= thin) ; enddo + else + do k=1,nk ; PCM_layer(k) = (dp(k) <= thin) ; enddo + endif + + !compute grid metrics + do k=2,nk-1 + qdpkm( K) = 1.0 / (dp(k-1) + dp(k)) + qdpkmkp(k) = 1.0 / (dp(k-1) + dp(k) + dp(k+1)) + dpkm2kp(k) = dp(k-1) + 2.0*dp(k) + dp(k+1) + enddo !k + qdpkm(nk) = 1.0 / (dp(nk-1) + dp(nk)) + + do i=1,ns + do K=2,nk + slope_edge(K) = qdpkm(K) * (s(k,i)-s(k-1,i)) + enddo !k + k = 1 !PCM first layer + edges(k,1,i) = s(k,i) + edges(k,2,i) = s(k,i) + zw(k,1) = 0.0 + zw(k,2) = 0.0 + do k=2,nk-1 + if ((slope_edge(K)*slope_edge(K+1) < 0.0) .or. PCM_layer(k)) then !use PCM + edges(k,1,i) = s(k,i) + edges(k,2,i) = s(k,i) + zw(k,1) = 0.0 + zw(k,2) = 0.0 + else + seh1 = dp(k)*slope_edge(K+1) + seh2 = dp(k)*slope_edge(K) + q01 = dpkm2kp(k)*slope_edge(K+1) + q02 = dpkm2kp(k)*slope_edge(K) + if (abs(seh1) > abs(q02)) then + seh1 = q02 + endif + if (abs(seh2) > abs(q01)) then + seh2 = q01 + endif + curv_cell = (seh1 - seh2) * qdpkmkp(k) + q001 = seh1 - curv_cell*dp(k+1) + q002 = seh2 + curv_cell*dp(k-1) + ! q001 = (seh1 * (dp(k-1) + dp(k)) + seh2 * dp(k+1)) * qdpkmkp(k) + ! q002 = (seh2 * (dp(k+1) + dp(k)) + seh1 * dp(k-1)) * qdpkmkp(k) + + edges(k,2,i) = s(k,i) + q001 + edges(k,1,i) = s(k,i) - q002 + zw(k,1) = (2.0*q001 - q002)**2 + zw(k,2) = (2.0*q002 - q001)**2 + endif !PCM:WENO + enddo !k + k = nk !PCM last layer + edges(k,1,i) = s(k,i) + edges(k,2,i) = s(k,i) + zw(k, 1) = 0.0 + zw(k, 2) = 0.0 + + do k=2,nk + ! This was the original code based on that in Hycom, but because zw has + ! dimensions of [A2], it can not use a constant (hard coded) value of dsmll. + ! ds2a = max(zw(k-1,2), dsmll) + ! ds2b = max(zw(k, 1), dsmll) + ! val_edge(K) = (ds2b*edges(k-1,2,i)+ds2a*edges(k,1,i)) / (ds2b+ds2a) + ! Use a weighted average of the two layers' estimated edge values as the actual edge value. + if (zw(k,1) + zw(k-1,2) <= 0.0) then + wt1 = 0.5 + elseif (zw(k,1) <= min_ratio * (zw(k,1) + zw(k-1,2))) then + wt1 = min_ratio + elseif (zw(k-1,2) <= min_ratio * (zw(k,1) + zw(k-1,2))) then + wt1 = (1.0 - min_ratio) + else + wt1 = zw(k,1) / (zw(k,1) + zw(k-1,2)) + endif + val_edge(k) = wt1*edges(k-1,2,i) + (1.0-wt1)*edges(k,1,i) + enddo !k + val_edge( 1) = 2.0*s( 1,i)-val_edge( 2) !not used? + val_edge(nk+1) = 2.0*s(nk,i)-val_edge(nk) !not used? + + do k=2,nk-1 + if (.not.PCM_layer(k)) then !don't use PCM + q01 = val_edge(K+1) - s(k,i) + q02 = s(k,i) - val_edge(K) + if (q01*q02 < 0.0) then + q01 = 0.0 + q02 = 0.0 + elseif (abs(q01) > abs(2.0*q02)) then + q01 = 2.0*q02 + elseif (abs(q02) > abs(2.0*q01)) then + q02 = 2.0*q01 + endif + edges(k,1,i) = s(k,i) - q02 + edges(k,2,i) = s(k,i) + q01 + endif ! PCM:WENO + enddo !k + enddo !i + +end subroutine hybgen_weno_coefs + +end module MOM_hybgen_remap diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index dafd165245..008475b16d 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -753,7 +753,8 @@ end subroutine end_regridding !------------------------------------------------------------------------------ !> Dispatching regridding routine for orchestrating regridding & remapping -subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h, conv_adjust) +subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_shelf_h, & + conv_adjust, PCM_cell) !------------------------------------------------------------------------------ ! This routine takes care of (1) building a new grid and (2) remapping between ! the old grid and the new grid. The creation of the new grid can be based @@ -783,15 +784,16 @@ subroutine regridding_main( remapCS, CS, G, GV, h, tv, h_new, dzInterface, frac_ real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in position of each interface real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: frac_shelf_h !< Fractional ice shelf coverage logical, optional, intent(in ) :: conv_adjust !< If true, do convective adjustment + logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + optional, intent(out ) :: PCM_cell !< Use PCM remapping in cells where true + ! Local variables real :: trickGnuCompiler - logical :: use_ice_shelf logical :: do_convective_adjustment do_convective_adjustment = .true. if (present(conv_adjust)) do_convective_adjustment = conv_adjust - - use_ice_shelf = present(frac_shelf_h) + if (present(PCM_cell)) PCM_cell(:,:,:) = .false. select case ( CS%regridding_scheme ) diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 1b3c5884de..82087eea24 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -13,12 +13,12 @@ module MOM_remapping use PLM_functions, only : PLM_reconstruction, PLM_boundary_extrapolation use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation use PQM_functions, only : PQM_reconstruction, PQM_boundary_extrapolation_v1 +use MOM_hybgen_remap, only : hybgen_plm_coefs, hybgen_ppm_coefs, hybgen_weno_coefs + use MOM_io, only : stdout, stderr implicit none ; private -#include - !> Container for remapping parameters type, public :: remapping_CS private @@ -46,11 +46,14 @@ module MOM_remapping ! The following are private parameter constants integer, parameter :: REMAPPING_PCM = 0 !< O(h^1) remapping scheme -integer, parameter :: REMAPPING_PLM = 1 !< O(h^2) remapping scheme -integer, parameter :: REMAPPING_PPM_H4 = 2 !< O(h^3) remapping scheme -integer, parameter :: REMAPPING_PPM_IH4 = 3 !< O(h^3) remapping scheme -integer, parameter :: REMAPPING_PQM_IH4IH3 = 4 !< O(h^4) remapping scheme -integer, parameter :: REMAPPING_PQM_IH6IH5 = 5 !< O(h^5) remapping scheme +integer, parameter :: REMAPPING_PLM = 2 !< O(h^2) remapping scheme +integer, parameter :: REMAPPING_PLM_HYBGEN = 3 !< O(h^2) remapping scheme +integer, parameter :: REMAPPING_PPM_H4 = 4 !< O(h^3) remapping scheme +integer, parameter :: REMAPPING_PPM_IH4 = 5 !< O(h^3) remapping scheme +integer, parameter :: REMAPPING_PPM_HYBGEN = 6 !< O(h^3) remapping scheme +integer, parameter :: REMAPPING_WENO_HYBGEN= 7 !< O(h^3) remapping scheme +integer, parameter :: REMAPPING_PQM_IH4IH3 = 8 !< O(h^4) remapping scheme +integer, parameter :: REMAPPING_PQM_IH6IH5 = 9 !< O(h^5) remapping scheme integer, parameter :: INTEGRATION_PCM = 0 !< Piecewise Constant Method integer, parameter :: INTEGRATION_PLM = 1 !< Piecewise Linear Method @@ -60,11 +63,14 @@ module MOM_remapping character(len=40) :: mdl = "MOM_remapping" !< This module's name. !> Documentation for external callers -character(len=256), public :: remappingSchemesDoc = & +character(len=360), public :: remappingSchemesDoc = & "PCM (1st-order accurate)\n"//& "PLM (2nd-order accurate)\n"//& + "PLM_HYBGEN (2nd-order accurate)\n"//& "PPM_H4 (3rd-order accurate)\n"//& "PPM_IH4 (3rd-order accurate)\n"//& + "PPM_HYBGEN (3rd-order accurate)\n"//& + "WENO_HYBGEN (3rd-order accurate)\n"//& "PQM_IH4IH3 (4th-order accurate)\n"//& "PQM_IH6IH5 (5th-order accurate)\n" character(len=3), public :: remappingDefaultScheme = "PLM" !< Default remapping method @@ -185,41 +191,50 @@ function isPosSumErrSignificant(n1, sum1, n2, sum2) end function isPosSumErrSignificant !> Remaps column of values u0 on grid h0 to grid h1 assuming the top edge is aligned. -subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edge) +subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1, h_neglect, h_neglect_edge, PCM_cell) type(remapping_CS), intent(in) :: CS !< Remapping control structure 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), intent(in) :: h0 !< Cell widths on source grid [H] + real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid [A] integer, intent(in) :: n1 !< Number of cells on target grid - real, dimension(n1), intent(in) :: h1 !< Cell widths on target grid - real, dimension(n1), intent(out) :: u1 !< Cell averages on target grid + real, dimension(n1), intent(in) :: h1 !< Cell widths on target grid [H] + real, dimension(n1), intent(out) :: u1 !< Cell averages on target grid [A] real, optional, intent(in) :: h_neglect !< A negligibly small width for the !! purpose of cell reconstructions - !! in the same units as h0. + !! in the same units as h0 [H] real, optional, intent(in) :: h_neglect_edge !< A negligibly small width !! for the purpose of edge value - !! calculations in the same units as h0. + !! calculations in the same units as h0 [H] + logical, dimension(n0), optional, intent(in) :: PCM_cell !< If present, use PCM remapping for + !! cells in the source grid where this is true. + ! Local variables integer :: iMethod - real, dimension(n0,2) :: ppoly_r_E !Edge value of polynomial - real, dimension(n0,2) :: ppoly_r_S !Edge slope of polynomial - real, dimension(n0,CS%degree+1) :: ppoly_r_coefs !Coefficients of polynomial + real, dimension(n0,2) :: ppoly_r_E ! Edge value of polynomial + real, dimension(n0,2) :: ppoly_r_S ! Edge slope of polynomial + real, dimension(n0,CS%degree+1) :: ppoly_r_coefs ! Coefficients of polynomial + real :: edges(n0,2) ! Interpolation edge values [A] + real :: slope(n0) ! Interpolation slopes per cell width [A] + real :: h0tot, h0err ! Sum of source cell widths and round-off error in this sum [H] + real :: h1tot, h1err ! Sum of target cell widths and round-off error in this sum [H] + real :: u0tot, u0err ! Integrated values on the source grid and round-off error in this sum [H A] + real :: u1tot, u1err ! Integrated values on the target grid and round-off error in this sum [H A] + real :: u0min, u0max, u1min, u1max ! Extrema of values on the two grids [A] + real :: uh_err ! Difference in the total amounts on the two grids [H A] + real :: hNeglect, hNeglect_edge ! Negligibly small cell widths in the same units as h0 [H] integer :: k - real :: eps, h0tot, h0err, h1tot, h1err, u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err - real :: hNeglect, hNeglect_edge hNeglect = 1.0e-30 ; if (present(h_neglect)) hNeglect = h_neglect hNeglect_edge = 1.0e-10 ; if (present(h_neglect_edge)) hNeglect_edge = h_neglect_edge call build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, ppoly_r_E, ppoly_r_S, iMethod, & - hNeglect, hNeglect_edge ) + hNeglect, hNeglect_edge, PCM_cell ) if (CS%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, CS%degree, & - CS%boundary_extrapolation, ppoly_r_coefs, ppoly_r_E, ppoly_r_S) - + CS%boundary_extrapolation, ppoly_r_coefs, ppoly_r_E, ppoly_r_S) call remap_via_sub_cells( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, n1, h1, iMethod, & - CS%force_bounds_in_subcell, u1, uh_err ) + CS%force_bounds_in_subcell, u1, uh_err ) if (CS%check_remapping) then ! Check errors and bounds @@ -306,7 +321,7 @@ subroutine remapping_core_w( CS, n0, h0, u0, n1, dx, u1, h_neglect, h_neglect_ed endif enddo call remap_via_sub_cells( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, n1, h1, iMethod, & - CS%force_bounds_in_subcell,u1, uh_err ) + CS%force_bounds_in_subcell, u1, uh_err ) ! call remapByDeltaZ( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, n1, dx, iMethod, u1, hNeglect ) ! call remapByProjection( n0, h0, u0, CS%ppoly_r, n1, h1, iMethod, u1, hNeglect ) @@ -353,7 +368,7 @@ end subroutine remapping_core_w !> Creates polynomial reconstructions of u0 on the source grid h0. subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & ppoly_r_E, ppoly_r_S, iMethod, h_neglect, & - h_neglect_edge ) + h_neglect_edge, PCM_cell ) type(remapping_CS), intent(in) :: CS !< Remapping control structure integer, intent(in) :: n0 !< Number of cells on source grid real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid @@ -369,10 +384,14 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & real, optional, intent(in) :: h_neglect_edge !< A negligibly small width !! for the purpose of edge value !! calculations in the same units as h0. + logical, optional, intent(in) :: PCM_cell(n0) !< If present, use PCM remapping for + !! cells from the source grid where this is true. + ! Local variables integer :: local_remapping_scheme integer :: remapping_scheme !< Remapping scheme logical :: boundary_extrapolation !< Extrapolate at boundaries if true + integer :: k, n ! Reset polynomial ppoly_r_E(:,:) = 0.0 @@ -398,6 +417,16 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & call PLM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect) endif iMethod = INTEGRATION_PLM + case ( REMAPPING_PLM_HYBGEN ) + call hybgen_PLM_coefs(u0, h0, ppoly_r_coefs(:,2), n0, 1, h_neglect) + do k=1,n0 + ppoly_r_E(k,1) = u0(k) - 0.5 * ppoly_r_coefs(k,2) ! Left edge value of cell k + ppoly_r_E(k,2) = u0(k) + 0.5 * ppoly_r_coefs(k,2) ! Right edge value of cell k + ppoly_r_coefs(k,1) = ppoly_r_E(k,1) + enddo + if ( CS%boundary_extrapolation ) & + call PLM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) + iMethod = INTEGRATION_PLM case ( REMAPPING_PPM_H4 ) call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 ) call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect, answers_2018=CS%answers_2018 ) @@ -412,6 +441,18 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) endif iMethod = INTEGRATION_PPM + case ( REMAPPING_PPM_HYBGEN ) + call hybgen_PPM_coefs(u0, h0, ppoly_r_E, n0, 1, h_neglect) + call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect, answers_2018=.false. ) + if ( CS%boundary_extrapolation ) & + call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) + iMethod = INTEGRATION_PPM + case ( REMAPPING_WENO_HYBGEN ) + call hybgen_weno_coefs(u0, h0, ppoly_r_E, n0, 1, h_neglect) + call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect, answers_2018=.false. ) + if ( CS%boundary_extrapolation ) & + call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefs, h_neglect ) + iMethod = INTEGRATION_PPM case ( REMAPPING_PQM_IH4IH3 ) call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E, h_neglect_edge, answers_2018=CS%answers_2018 ) call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_S, h_neglect, answers_2018=CS%answers_2018 ) @@ -437,6 +478,16 @@ subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefs, & 'The selected remapping method is invalid' ) end select + if (present(PCM_cell)) then + ! Change the coefficients to those for the piecewise constant method in indicated cells. + do k=1,n0 ; if (PCM_cell(k)) then + ppoly_r_coefs(k,1) = u0(k) + ppoly_r_E(k,1:2) = u0(k) + ppoly_r_S(k,1:2) = 0.0 + do n=2,CS%degree+1 ; ppoly_r_coefs(k,n) = 0.0 ; enddo + endif ; enddo + endif + end subroutine build_reconstructions_1d !> Checks that edge values and reconstructions satisfy bounds @@ -1580,12 +1631,21 @@ subroutine setReconstructionType(string,CS) case ("PLM") CS%remapping_scheme = REMAPPING_PLM degree = 1 + case ("PLM_HYBGEN") + CS%remapping_scheme = REMAPPING_PLM_HYBGEN + degree = 1 case ("PPM_H4") CS%remapping_scheme = REMAPPING_PPM_H4 degree = 2 case ("PPM_IH4") CS%remapping_scheme = REMAPPING_PPM_IH4 degree = 2 + case ("PPM_HYBGEN") + CS%remapping_scheme = REMAPPING_PPM_HYBGEN + degree = 2 + case ("WENO_HYBGEN") + CS%remapping_scheme = REMAPPING_WENO_HYBGEN + degree = 2 case ("PQM_IH4IH3") CS%remapping_scheme = REMAPPING_PQM_IH4IH3 degree = 4 diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 41ba70f152..5fc168f5a4 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -360,8 +360,8 @@ subroutine open_boundary_config(G, US, param_file, OBC) character(len=128) :: inputdir logical :: answers_2018, default_2018_answers logical :: check_reconstruction, check_remapping, force_bounds_in_subcell - character(len=32) :: remappingScheme -! This include declares and sets the variable "version". + character(len=64) :: remappingScheme + ! This include declares and sets the variable "version". # include "version_variable.h" allocate(OBC) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 8378644ea9..45020c86f5 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2450,7 +2450,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just ! extrapolating the densities at the bottom of unstable profiles ! from data when finding the initial interface locations in ! layered mode from a dataset of T and S. - character(len=10) :: remappingScheme + character(len=64) :: remappingScheme real :: tempAvg, saltAvg integer :: nPoints, ans integer :: id_clock_routine, id_clock_read, id_clock_interp, id_clock_fill, id_clock_ALE diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index 8a67d71fe2..9b6800524b 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -58,7 +58,7 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ character(len=200) :: mesg real :: convert integer :: recnum - character(len=10) :: remapScheme + character(len=64) :: remapScheme logical :: homog,useALE ! This include declares and sets the variable "version". diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index ab3621296f..af9534e943 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -142,7 +142,7 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data, res real :: nhours_incupd, dt, dt_therm type(vardesc) :: vd character(len=256) :: mesg - character(len=10) :: remapScheme + character(len=64) :: remapScheme if (.not.associated(CS)) then call MOM_error(WARNING, "initialize_oda_incupd called without an associated "// & "control structure.") diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 472ee21e36..6b89a86b30 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -163,16 +163,18 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, ! time at v-points [T-1 ~> s-1]. -! This include declares and sets the variable "version". -#include "version_variable.h" + ! Local variables character(len=40) :: mdl = "MOM_sponge" ! This module's name. - logical :: use_sponge real, allocatable, dimension(:,:) :: Iresttime_u !< inverse of the restoring time at u points [T-1 ~> s-1] real, allocatable, dimension(:,:) :: Iresttime_v !< inverse of the restoring time at v points [T-1 ~> s-1] + ! This include declares and sets the variable "version". +# include "version_variable.h" + character(len=64) :: remapScheme + logical :: use_sponge logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries logical :: default_2018_answers integer :: i, j, k, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v - character(len=10) :: remapScheme + if (associated(CS)) then call MOM_error(WARNING, "initialize_ALE_sponge_fixed called with an associated "// & "control structure.") @@ -422,17 +424,18 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, param_file, CS, Irest real, dimension(SZI_(G),SZJB_(G)), intent(in), optional :: Iresttime_v_in !< The inverse of the restoring time !! for v [T-1 ~> s-1]. -! This include declares and sets the variable "version". -#include "version_variable.h" + ! Local variables character(len=40) :: mdl = "MOM_sponge" ! This module's name. - logical :: use_sponge real, allocatable, dimension(:,:) :: Iresttime_u !< inverse of the restoring time at u points [T-1 ~> s-1] real, allocatable, dimension(:,:) :: Iresttime_v !< inverse of the restoring time at v points [T-1 ~> s-1] + ! This include declares and sets the variable "version". +# include "version_variable.h" + character(len=64) :: remapScheme + logical :: use_sponge logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries logical :: default_2018_answers logical :: spongeDataOngrid = .false. integer :: i, j, k, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v - character(len=10) :: remapScheme if (associated(CS)) then call MOM_error(WARNING, "initialize_ALE_sponge_varying called with an associated "// & diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index d12d850a73..f7378bb37b 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -77,8 +77,8 @@ module MOM_CVMix_KPP real :: cs2 !< Parameter for multiplying by non-local term ! This is active for NLT_SHAPE_CUBIC_LMD only logical :: enhance_diffusion !< If True, add enhanced diffusivity at base of boundary layer. - character(len=10) :: interpType !< Type of interpolation to compute bulk Richardson number - character(len=10) :: interpType2 !< Type of interpolation to compute diff and visc at OBL_depth + character(len=32) :: interpType !< Type of interpolation to compute bulk Richardson number + character(len=32) :: interpType2 !< Type of interpolation to compute diff and visc at OBL_depth logical :: computeEkman !< If True, compute Ekman depth limit for OBLdepth logical :: computeMoninObukhov !< If True, compute Monin-Obukhov limit for OBLdepth logical :: passiveMode !< If True, makes KPP passive meaning it does NOT alter the diffusivity