diff --git a/.gitmodules b/.gitmodules index 01c7283aa319..0610f98d3e22 100644 --- a/.gitmodules +++ b/.gitmodules @@ -61,3 +61,6 @@ [submodule "components/mpas-ocean/src/FFTW"] path = components/mpas-ocean/src/FFTW url = https://github.com/knbarton/fftw-3.3.8.git +[submodule "components/mpas-ocean/src/ppr"] + path = components/mpas-ocean/src/ppr + url = https://github.com/dengwirda/PPR.git diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist index 5fee738461f8..64a922e9781d 100755 --- a/components/mpas-ocean/bld/build-namelist +++ b/components/mpas-ocean/bld/build-namelist @@ -768,11 +768,14 @@ add_default($nl, 'config_land_ice_flux_jenkins_salt_transfer_coefficient'); # Namelist group: advection # ############################# -add_default($nl, 'config_vert_tracer_adv'); -add_default($nl, 'config_vert_tracer_adv_order'); +add_default($nl, 'config_vert_advection_method'); +add_default($nl, 'config_vert_remap_order'); +add_default($nl, 'config_vert_remap_interval'); +add_default($nl, 'config_vert_tracer_adv_flux_order'); add_default($nl, 'config_horiz_tracer_adv_order'); add_default($nl, 'config_coef_3rd_order'); -add_default($nl, 'config_monotonic'); +add_default($nl, 'config_flux_limiter'); +add_default($nl, 'config_remap_limiter'); ############################### # Namelist group: bottom_drag # diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section index d7d4efe740f0..0a9d5459111c 100644 --- a/components/mpas-ocean/bld/build-namelist-section +++ b/components/mpas-ocean/bld/build-namelist-section @@ -301,11 +301,14 @@ add_default($nl, 'config_land_ice_flux_jenkins_salt_transfer_coefficient'); # Namelist group: advection # ############################# -add_default($nl, 'config_vert_tracer_adv'); -add_default($nl, 'config_vert_tracer_adv_order'); +add_default($nl, 'config_vert_advection_method'); +add_default($nl, 'config_vert_remap_order'); +add_default($nl, 'config_vert_remap_interval'); +add_default($nl, 'config_vert_tracer_adv_flux_order'); add_default($nl, 'config_horiz_tracer_adv_order'); add_default($nl, 'config_coef_3rd_order'); -add_default($nl, 'config_monotonic'); +add_default($nl, 'config_flux_limiter'); +add_default($nl, 'config_remap_limiter'); ############################### # Namelist group: bottom_drag # diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml index 789acecc6e54..a1c31bf79ab2 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -342,11 +342,14 @@ 8.42e-5 -'stencil' -3 +'flux-form' +3 +0 +3 3 0.25 -.true. +'monotonic' +'monotonic' .true. diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml index 6e0153ead725..841de7c4e7fc 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -1582,17 +1582,33 @@ Default: Defined in namelist_defaults.xml - -Method for interpolating tracer values from layer centers to layer edges +Method for advecting tracers, momentum, and thickness vertically. -Valid values: 'spline' and 'stencil' +Valid values: 'flux-form' and 'remap' Default: Defined in namelist_defaults.xml - -Order of polynomial used for tracer reconstruction at layer edges +Order of remapping method used for momentum and tracer advection + +Valid values: 1, 2, 3 and 5 +Default: Defined in namelist_defaults.xml + + + +Number of timesteps between each remapping. If 0, remapping occurs every timestep + +Valid values: Any integer greater than or equal to 0 +Default: Defined in namelist_defaults.xml + + + +Order of polynomial used for tracer reconstruction at layer edges for flux-form method Valid values: 2, 3 and 4 Default: Defined in namelist_defaults.xml @@ -1614,11 +1630,19 @@ Valid values: any real between 0 and 1 Default: Defined in namelist_defaults.xml - -If .true. then fluxes are limited to produce a monotonic advection scheme +Slope limiter for the flux-form advection scheme. -Valid values: .true. and .false. +Valid values: 'none','monotonic' +Default: Defined in namelist_defaults.xml + + + +Slope limiter for the vertical remap advection scheme. + +Valid values: 'none','monotonic','weno' Default: Defined in namelist_defaults.xml @@ -2188,7 +2212,7 @@ Default: Defined in namelist_defaults.xml -Enables a change on tracer monotonicity at the end of the monotonic advection routine. Only used if config_monotonic is set to .true. +Enables a change on tracer monotonicity at the end of the monotonic advection routine. Only used if config_flux_limiter is set to monotonic Valid values: .true. or .false. Default: Defined in namelist_defaults.xml diff --git a/components/mpas-ocean/docs/design_docs/index.rst b/components/mpas-ocean/docs/design_docs/index.rst index 703a9d4efabb..190686a4a264 100644 --- a/components/mpas-ocean/docs/design_docs/index.rst +++ b/components/mpas-ocean/docs/design_docs/index.rst @@ -8,3 +8,4 @@ Design document describing new capabilities added to MPAS-Ocean. time-varying-wind prognostic-sediment-transport + vert-lagrang-remap/vert-lagrang-remap diff --git a/components/mpas-ocean/docs/design_docs/vert-lagrang-remap/vert-lagrang-remap.rst b/components/mpas-ocean/docs/design_docs/vert-lagrang-remap/vert-lagrang-remap.rst new file mode 100644 index 000000000000..a2328e75be2c --- /dev/null +++ b/components/mpas-ocean/docs/design_docs/vert-lagrang-remap/vert-lagrang-remap.rst @@ -0,0 +1,582 @@ + +Vertical Lagrangian-remap method +================================ + +date: 2020/12/21 +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + + + +Summary +------- + +We propose to implement the vertical Lagrangian-remap method to support the +choice of hybrid vertical coordinates, which may be required for evolving ice- +shelf cavity geometries. It can be seen as an extension of the Arbitrary +Lagrangian-Eulerian method which is already implemented in MPAS-Ocean, a key +difference being that the solution of diasurface transport is coordinate-free. +This has the advantage of removing the vertical CFL constraint, which is of +particular concern for an evolving ice shelf geometry in terrain-following +coordinate schemes. This development may also reduce numerical diffusion in the +vertical, but we do not require this result for this implementation to be +successful. + +Generalization of the vertical coordinate is not in the scope of this development. +Additionally, the infrastructure for specifying a target grid that differs from +the existing z-star and z-tilde options is reserved for future development. +The minimum criteria for success are not degrading the accuracy of the solution +or computational performance. + + +Requirements +------------ + +Requirement: Ability to perform vertical Lagrangian-remapping +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +The core of the proposed development is the ability to perform vertical +Lagrangian-remapping (VLR). When VLR is performed, the momentum, volume and +tracer conservation equations are first solved in a Lagrangian sense; that is, +layer thicknesses evolve such that there is no vertical transport across layer +interfaces. A target grid is then specified. For this stage of development, +we only require that the algorithm be able to use the existing grid +specifications. However, the implementation should be general enough that +a different grid could be used for remapping. This regridding step is followed +by a conservative remapping of velocity and scalars to the target grid. + +Requirement: Grid is sufficiently conditioned for accuracy and stability +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +A severely deformed grid can lead to issues with the accuracy and stability of +the numerical methods used for the dynamics and thermodynamics. If layers are +allowed to become very thin they can also degrade the accuracy and stability of +the solution. Thus, a requirement of the VLR implementation is that there are +safeguards to prevent layers from becoming too thin or deformed. + +Requirement: Ability to specify remapping method and its options +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +There are several choices to make for the remapping operation to balance +accuracy, stability and computational cost. To the extent that it is feasible, +the user should have the ability to control the order of accuracy of the +remapping method, and other remapping options (e.g., via namelist options). + +Requirement: Support for existing time stepping schemes +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Vertical Lagrangian remapping should be supported (eventually) in both the RK4 +and the split-explicit time stepping schemes currently implemented in MPAS- +Ocean. The scheme will also have a straightforward interface so it can be easily +added to new time stepping schemes. + +Requirement: not climate changing +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +For the same target grid, time-stepping scheme and problem, the solution for +layer thicknesses and key prognostic variables using VLR is within a non- +climate-changing threshold value of the solution when VLR is not performed. +Tests should span from idealized tests of simplified dynamics to full climate +simulations. Acceptable thresholds specified for each test could be similar to +those from changing time stepping schemes. + +Requirement: tracer and volume conservation +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2020/12/04 + +Contributors: Darren Engwirda, Xylar Asay-Davis, Carolyn Begeman, Mark Petersen + +The column-wise tracer content should be conserved over the remapping step to a +specified threshold, and using a specified method to integrate the tracer. The +column-integrated volume must remain unchanged over the remapping step, as +expected. + +Requirement: performance +^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2020/12/04 + +Contributors: Darren Engwirda, Xylar Asay-Davis, Carolyn Begeman, Mark Petersen + +Performance is not degraded below an acceptable level. Compute time for +stand-alone ocean, without i/o, on full primitive equations, is expected to +remain very similar - within 5% of compute time when using vertical +advection method but likely less because remapping is a column-wise, in-cache +operation. + +Requirement: no interference +^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Remapping does not introduce errors in the computation of prognostic and +diagnostic variables, such as through the use of an Eulerian vertical velocity +that is not consistent with the Lagrangian solution. + +Requirement: modularity +^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2020/12/21 + +Contributors: Darren Engwirda, Xylar Asay-Davis, Carolyn Begeman + +To the degree possible, the code for determining the target grid and performing +remapping should be kept in its own Fortran module(s) for better readability. +These modules should be called by both timestepping routines to maximize code +reuse. + + +Algorithm Design (optional) +--------------------------- + +Algorithm Design: Ability to perform vertical Lagrangian-remapping +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2020/12/15 + +Contributors: Darren Engwirda, Xylar Asay-Davis, Carolyn Begeman + +The conservation of momentum, volume, and tracer equations in MPAS-Ocean ( +`Ringler et al. 2013 `_; +`Petersen et al. 2014 `_) +are: + +.. math:: + + \frac{\partial u_k}{\partial t} + q_k h_k u_k^{normal} + \overline{w^t \delta z^t u} = -\frac{1}{\rho_0} \nabla p_k - \frac{\rho_k g}{\rho_0} \nabla z_k - \nabla K_k + [D_h^u]_k + [D_{\nu}^u]_k + F_k^u + + \frac{\partial h_k}{\partial t} + \nabla \cdot (h_k \mathbf{u}_k) + w_k^t - w_{k+1}^t = 0 + + \frac{\partial (h_k \phi_k)}{\partial t} + \nabla \cdot (h_k \mathbf{u}_k \phi_k) + \overline{\phi}_k^t w_k^t - \overline{\phi}_{k+1}^t w_{k+1}^t = [D_h^{\phi}]_k + [D_v^{\phi}]_k + F_k^{\phi} + +For the Lagrangian step, the vertical velocity through the top of the cell, +:math:`w_k^t`, is set to zero in all of the above equations. Thus, these +equations simplify to: + +.. math:: + + \frac{\partial u_k}{\partial t} + q_k h_k u_k^{\perp} = -\frac{1}{\rho_0} \nabla p_k - \frac{\rho_k g}{\rho_0} \nabla z_k - \nabla K_k + [D_h^u]_k + [D_v^u]_k + F_k^u + + \frac{\partial h_k}{\partial t} + \nabla \cdot (h_k \mathbf{u}_k) = 0 + + \frac{\partial (h_k \phi_k)}{\partial t} + \nabla \cdot (h_k \mathbf{u}_k \phi_k) = [D_h^{\phi}]_k + [D_v^{\phi}]_k + F_k^{\phi} + +The time-stepping algorithm (RK4 or split-explicit) advances the prognostic +variables and layer thickness from their values at time n +:math:`u_k^{n},\phi_k^{n},h_k^{n}`, to their values after the Lagrangian step, +designated by the superscript *lg*, :math:`u_k^{lg},h_k^{lg},\phi_k^{lg}`. + +Note that the vertical mixing terms :math:`D_v^h, D_v^{\phi}` +are retained here. We opt to compute these terms prior to remapping as this +allow for future development in which the dynamics are subcycled relative to +the thermodynamics and remapping is scheduled on the thermodynamic timestep. +This computation of vertical mixing terms prior to remapping is similar to +both MOM6 and HYCOM. We anticipate that there could be a trade-off between (a) +loss of accuracy of vertical mixing terms when their computation precedes +remapping due to grid deformation and (b) loss of accuracy when their +computation follows remapping due to remapping errors in vertical gradients of +prognostic variables. We do not intend to test this at this time. + +The target grid needs to be determined after the solution for prognostic +variables so that the vertical Lagrangian-remapping method is general enough to +be used with coordinate systems that depend on the ocean state (this includes +the z-star coordinate system in which SSH perturbations are vertically +distributed between layers). We do not present an algorithmic design for +regridding to coordinate systems not already supported in MPAS-Ocean, as this +will be the subject of future development. For now, the target grid is based on a +constant set of z-star levels that are specified at initialization. + +For the grid selection step, the target grid, :math:`h_k^{target}`, is +determined, conserving volume: + +.. math:: + + \sum_{k=1}^{kmax}h_k^{target} = \sum_{k=1}^{kmax}h_k^{lg} + + +For scalar remapping, layer thicknesses at the next timestep, +:math:`h_k^{n+1}` are set according to the target grid and scalars are remapped +to the target grid using the remapping operations represented by the function +:math:`G`: + +.. math:: + + h_k^{n+1} = h_k^{target} + + hEdge_k^{n+1} = 0.5 (h_{k,cell1}^{n+1} + h_{k,cell2}^{n+1}) + + u_k^{n+1} = G(u_k^{lg},hEdge_k^{lg},hEdge_k^{n+1}) + + \phi_k^{n+1} = G(\phi_k^{lg},\phi_k^{lg},h_k^{n+1}) + +For velocity remapping, we solve for layer thicknesses at edges after the +lagrangian step and the regridded thickness. In this development, we only plan +to support centered edge layer thicknesses consistent with the centered +advection scheme. There does not appear to be a precedent among ocean models +(HYCOM, MOM6) at this time for using upwinded layer thickness in the remapping +operation. We touch on a few of the implementation challenges with using +upwinded layer thicknesses for remapping in the corresponding implementation +section. If VLR is run with an upwinded thickness flux, the horizontal momentum +flux will not be conserved as :math:`hEdge^{n+1}` will be reassigned to the +upwinded layer thickness (errors will likely increase as horizontal gradients +in layer thickness increase). Otherwise, the remapping operation :math:`G` +conserves volume flux and scalar concentration. + +.. math:: + + hEdge_k^{lg} = 0.5 (h_{k,cell1}^{lg} + h_{k,cell2}^{lg}) + + hEdge_k^{n+1} = 0.5 (h_{k,cell1}^{n+1} + h_{k,cell2}^{n+1}) + + u_k^{n+1} = G(u_k^{lg},hEdge_k^{lg},hEdge_k^{n+1}) + + +.. math:: + + \sum_{k=1}^{kmax} u_k^{n+1} h_k^{n+1} = \sum_{k=1}^{kmax} u_k^{lg} h_k^{lg} + + \sum_{k=1}^{kmax} \phi_k^{n+1} h_k^{n+1} = \sum_{k=1}^{kmax} \phi_k^{lg} h_k^{lg} + +The vertical velocity across layer interfaces may be computed anytime after +regridding. It can be computed as + +.. math:: + + w = - \nabla \cdot (h_k \mathbf{u}_k) - (h_k^{t+1} - h_k^t)/dt + +or + +.. math:: + + w = (h_k^{t+1} - h_k^{lg})/dt + +The choice between the two is discussed in the Implementation section. + + +Implementation +-------------- + +Implementation: Ability to perform vertical Lagrangian-remapping +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Namelist options: + +- To turn VLR on/off: + :code:`advection, config_vert_advection_method = 'remap' or 'flux-form'` + + +Lagrangian step: + +The solution for prognostic variables in RK4 and split-explicit remains +largely the same. The main difference is that the vertical velocity through +the top of layers is set to zero in the routine +:code:`ocn_vert_transport_velocity_top`. This is the same as what is done when +:code:`config_vert_coord_movement` is :code:`impermeable_interfaces`. + + +Grid selection step: + +#. :math:`z_k^{target}`, the depth of the top of the layer, is determined based + on an analytical expression for the grid. + The simplest case is constant z-levels, :math:`z_k^{target} = z_k^{init}`. + Since :math:`z_k^{target}` can be a function of the ocean state (e.g., + :math:`\rho` for isopycnal coordinates, grid selection doesn't occur until + after the solution for prognostic variables. +#. Superimpose SSH perturbations according to one of the existing depth- + dependent functions, :math:`z_k^{target} = z_k^{target} + c(z) \: \eta`. As + in :code:`ocn_ALE_thickness`, layer thicknesses are adjusted from the + seafloor upwards. Currently :code:`layerThicknessTarget` is computed in a + new routine :code:`ocn_vert_regrid` following other instances in which a z- + star grid is recomputed. +#. Apply conditioning steps outlined in the following section. (Not + implemented until coordinates other than z-star are implemented.) +#. In preparation for remapping, compute :math:`z_k^{target}` from + :code:`layerThicknessTarget`. + +All of the grid selection steps will be performed from a separate module called +:code:`ocn_vert_regrid`. This topic is further addressed in section +Implementation: modularity. + + +Remapping step: + +There is a new remapping routine :code:`ocn_remap_vert_state` which updates +all state variables with depth coordinates at :code:`timeLevel=2`. On input, +:code:`layerThickness(tlev=2)` reflects the Lagrangian layer thickness +determined by :code:`ocn_tend_thick`. On output, it is equal to +:code:`layerThicknessNew` as determined by the regridding routine +:code:`ocn_vert_regrid`. + +Members of :code:`statePool` that will be remapped: + +- All members of :code:`tracerPool` unless :code:`activeTracersOnly`, in which + case only the :code:`activeTracers` will be remapped +- :code:`normalVelocity` +- :code:`highFreqThickness` +- :code:`lowFreqDivergence` +- :code:`normalBarotropicVelocity`, only for split-explicit time-stepping + +In preparation for remapping, we compute the scratch variables +:code:`layerThickEdgeNew` as the mean of :code:`layerThicknessNew` at +neighboring cells. We do these locally rather than through a call to +:code:`ocn_diagnostic_solve_layerThicknessEdge`. For the Lagrangian +:code:`layerThickEdgeMean` we use the existing solution from +:code:`ocn_diagnostic_solve_layerThicknessEdge`. This approach requires that +:code:`config_thickness_flux_type` is :code:`'centered'`. At initialization, +we throw an error but do not terminate the run if +:code:`config_thickness_flux_type` is not :code:`'centered'` and VLR is active. + + +A note about difficulties of implementing upwinded :code:`layerThicknessEdge` +for remapping: + +Currently, the PPR library assumes that the total height is the same before and +after remapping (i.e., :code:`sum_k(layerThicknessOld)` equals +:code:`sum_k(layerThicknessTarget)`. Over the course of remapping, the upwind +cell could change for one or more layers and thus the total column height could +change. PPR would have to be carefully adapted to deal with this condition in +order to preserve the total volume flux as well as the vertical distribution of +momentum during remapping. + +An alternative to modifying the remapping library is to use centered edge layer +thicknesses for remapping and correct :code:`uNormal` prior to remapping such +that :code:`uNormalCorr(k) * layerThicknessEdgeCntr(k) = uNormal(k) * layerThicknessEdgeUpwind(k)`. +When edge layer thicknesses are upwinded based on the remapped :code:`uNormal`, +:code:`uNormal` must be corrected again to preserve layerwise fluxes. There +will still be some error in the vertical distribution of volume flux with this +approach. Given the complexity of either of these implementation options, we +leave this issue for future development. + + +After determining the layer thicknesses to remap to, the vertical layer +interface locations are determined, :code:`heightCellNow`, +:code:`heightCellNew`, :code:`heightEdgeNow`, :code:`heightEdgeNew`. +These variables and the Lagrangian state variable are the inputs to the PPR +routine :code:`rmap1d`, and the output is the remapped state variable. + + +Some implementation considerations for PPR: + +- Error-checking in PPR: make consistent with MPAS errors, consider additional + error checks. (Not implemented) + +After remapping, :code:`ocn_diagnostic_solve` is called. This is needed to +compute the density and pressure fields based on the remapped ocean state and +the diagnostic field :code:`vertVelocityTop` which is the vertical velocity +through the top of the layer. This is only used as a diagnostic variable for +computing the MOC streamfunction and tracer budgets. None of the mixing +parameterizations require a vertical velocity (Eulerian or diasurface velocity). + +Note: if `vertVelocityTop` is computed between regridding and remapping then it +can be computed as + +.. code:: + + vertVelocityTop(k) = vertVelocityTop(k+1) - div_hu(k) - + (layerThickness(k,tlev=2) - layerThickness(k,tlev=1))/dt + +In our implementation, `vertVelocityTop` is computed after remapping and +:code:`div_hu` is no longer appropriate as it has been remapped. In this case, +the Lagrangian layer thickness is stored in the state variable +:code:`layerThicknessLag` and then the vertical velocity through the top of the +layer is computed as: + +.. code:: + + layerThicknessLag = layerThickness(tlev=2) + + layerThickness(tlev=2) = layerThicknessTarget + + vertVelocityTop = (layerThickness(tlev=2) - layerThicknessLag)/dt + +The function that performs this computation is +:code:`ocn_diagnostic_solve_vertVel_remap`. + +If :code:`normalGMBolusVelocity` is computed based on the remapped ocean state +then the computation of :code:`vertTransportVelocityTop` and +:code:`vertGMBolusVelocitytop` is unchanged as these fields represent Eulerian +velocities. + +However, the current implementation will not compute +:code:`normalGMBolusVelocity` based on the remapped ocean state before +:code:`ocn_diagnostic_solve` is called. Thus, :code:`vertGMBolusVelocityTop` and +:code:`vertTransportVelocityTop` will be inaccurate. This will only pose an +issue when the number of vertical levels changes during regridding, a +capability which isn't included in this development scope. + +Tracer tendencies that are computed as diagnostics will also be inaccurate +after regridding, as they will not be remapped. Remapping these variables does +not make physical sense without also computing vertical tracer fluxes, which +would be overly burdensome. Analysis members that currently use these diagnostic +variables are :code:`mpas_ocn_layer_volume_weighted_averages` and +:code:`mpas_ocn_mixed_layer_heat_budget`. + +Implementation: Grid is sufficiently conditioned for accuracy and stability +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Since the grid is remapped back to z-star in the current implementation, we +assume that the grid is sufficiently regular. + +For future coordinate choices, consdider applying the following after +determining the target grid: + +#. Assign :math:`h_k^{t+1}` to :math:`h_k^{lg}` if :math:`h_k^{t+1} - h_k^{lg}` + is less than a minimum thickness alteration. This is motivated by accuracy + considerations, as each remapping may introduce errors. (Not implemented) +#. Apply minimum layer thickness criterion. (Not implemented) +#. Smooth layers in space and time. (Not implemented) + +Namelist options: + +- Minimum layer thickness +- Minimum thickness change for remapping to occur + + +Implementation: Ability to specify remapping method and its options +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Namelist options: + +- Order of remapping: + :code:`advection, config_vert_remap_order = 1, 2, 3, or 5` +- Interval number of timesteps for remapping: + :code:`advection, config_vert_remap_interval = integer >= 0` +- Slope limiter: + :code:`advection, config_remap_limiter = 'none', 'monotonic', or 'weno'` + +Implementation: Support for existing time stepping schemes +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Currently, only the split-explicit time integration scheme is supported. We +explored the RK4 implementation but it was not found to be stable. + +Implementation: performance +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +OMP calls were implemented for all cell loops in the remapping routine. Our +main strategy for improving performance is remapping an interval of time steps +rather than every time step. For EC60to30 global configurations, an interval of +3-4 timesteps was found to not significantly degrade the solution. + +Options for improving performance (not implemented): + +- Splitting the scalar and momentum timesteps +- Only remapping when the change in thickness exceeds given threshold +- Optimizing/parallelizing PPR + +Implementation: no interference +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2020/12/15 + +Contributors: Darren Engwirda, Xylar Asay-Davis, Carolyn Begeman + +Ensure that no calculations are made outside of remapping itself using +prognostic (or diagnostic) variables both before and after remapping to avoid +introducing additional errors. + +Look for places in the code where prognostic variables are used at the previous +timestep. + +Implementation: modularity +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Remapping operations (PPR) are performed in a separate routine, +:code:`ocn_vert_remap_state`. + +Target grid levels are determined in a separate routine, +:code:`ocn_vert_regrid`. + + +Testing +------- + +Testing and Validation: Ability to perform vertical Lagrangian-remapping +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Ability to handle strong vertical velocities: + +- Baroclinic channel test case + +Evaluating spurious mixing due to remapping: Compare with and without VLR + +- Internal wave test case +- Baroclinic channel test case + + +Testing and Validation: not climate changing +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +Global ocean case :code:`GMPAS-JRA1p4.TL319_EC30to60E2r2` 50 year run comparison +between VLR branch with VLR on and master without remapping: +`Link to analysis `_ + +Testing and Validation: performance +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Date last modified: 2022/05/21 + +Contributors: Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis + +For the same global ocean case above we show that remapping increases run time +by 20% for remapping every time step and that it scales almost perfectly with as +the remapping interval increases, i.e., 5% increase for remapping every 4 time +steps. We also show that 90% of the remapping time is spent in the PPR routine +:code:`rmap1d`, indicating that this would be the target for performance +improvements. diff --git a/components/mpas-ocean/src/Makefile b/components/mpas-ocean/src/Makefile index 6ace38223a57..822ed12d5409 100644 --- a/components/mpas-ocean/src/Makefile +++ b/components/mpas-ocean/src/Makefile @@ -8,7 +8,7 @@ endif FRAMEWORK_DIR = ../../../mpas-framework/src OCEAN_SHARED_INCLUDES = -I$(FRAMEWORK_DIR)/framework -I$(FRAMEWORK_DIR)/external/esmf_time_f90 -I$(FRAMEWORK_DIR)/operators OCEAN_SHARED_INCLUDES += -I$(PWD)/shared -I$(PWD)/analysis_members -I$(PWD)/mode_forward -I$(PWD)/mode_analysis -OCEAN_SHARED_INCLUDES += -I$(PWD)/BGC -I$(PWD)/MARBL/include -I$(PWD)/cvmix/src/shared -I$(PWD)/gotm/build -I$(PWD)/gotm/build/modules -I$(PWD)/SHTNS +OCEAN_SHARED_INCLUDES += -I$(PWD)/BGC -I$(PWD)/MARBL/include -I$(PWD)/cvmix/src/shared -I$(PWD)/gotm/build -I$(PWD)/gotm/build/modules -I$(PWD)/SHTNS -I$(PWD)/ppr/src ifneq "$(EXCLUDE_INIT_MODE)" "true" OCEAN_SHARED_INCLUDES += -I$(PWD)/mode_init endif @@ -22,7 +22,7 @@ ifneq "$(EXCLUDE_INIT_MODE)" "true" OCEAN_SHARED_INCLUDES += -I$(PWD)/mode_init endif -all: shared libcvmix analysis_members libBGC libMARBL libgotm libfftw libshtns +all: shared libcvmix analysis_members libBGC libMARBL libgotm libfftw libshtns libppr (cd mode_forward; $(MAKE) FCINCLUDES="$(FCINCLUDES) $(OCEAN_SHARED_INCLUDES)" all ) (cd mode_analysis; $(MAKE) FCINCLUDES="$(FCINCLUDES) $(OCEAN_SHARED_INCLUDES)" all ) ifneq "$(EXCLUDE_INIT_MODE)" "true" @@ -131,7 +131,14 @@ libfftw: fi \ fi -shared: libcvmix libBGC libMARBL libgotm libshtns libfftw +libppr: + if [ -e ppr/.git ]; then \ + (cd ppr/src; $(CPP) $(CPPFLAGS) $(CPPINCLUDES) ppr_1d.F90 > ppr_1d.F; $(FC) $(FFLAGS) -c ppr_1d.F -o ppr_1d.o) \ + else \ + (pwd ; echo "Missing core_ocean/ppr/.git, did you forget to 'git submodule update --init --recursive' ?"; exit 1) \ + fi + +shared: libcvmix libBGC libMARBL libgotm libshtns libfftw libppr (cd shared; $(MAKE) FCINCLUDES="$(FCINCLUDES) $(OCEAN_SHARED_INCLUDES)") analysis_members: libcvmix shared libgotm diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 9e9128f39dda..fda2d17e651c 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1000,12 +1000,20 @@ - - + + - + @@ -1340,7 +1352,7 @@ possible_values=".true. or .false." /> + @@ -1687,6 +1700,25 @@ + + + + + + + + + + + + + \brief Computes diagnostic variables for velocity +!> \author Carolyn Begeman +!> \date October 2021 +!> \details +!> This routine computes the Eulerian vertical velocity valid for +!> Vertical Lagrangian Remapping +! +!----------------------------------------------------------------------- + subroutine ocn_diagnostic_solve_vertVel_remap(layerThickness, & + layerThicknessLag, dt, vertVelocityTop)!{{{ + + use ocn_mesh + + implicit none + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness ! layerThickness after remapping at timeLevelIn + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThicknessLag ! layerThickness before remapping at timeLevelIn + real (kind=RKIND), intent(in) :: dt + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + vertVelocityTop + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iCell, k + real (kind=RKIND), dimension(:), allocatable :: & + heightCell, heightCellLag + + allocate(heightCell(nVertLevels+1)) + allocate(heightCellLag(nVertLevels+1)) + + !$omp parallel + !$omp do schedule(runtime) private(k, heightCell, heightCellLag) + do iCell = 1, nCellsAll + + heightCell(:) = -1.0_RKIND * bottomDepth(iCell) + heightCellLag(:) = -1.0_RKIND * bottomDepth(iCell) + + do k = maxLevelCell(iCell), minLevelCell(iCell), -1 + ! reconstruct heights + heightCellLag(k) = heightCellLag(k + 1) + layerThicknessLag(k, iCell) + heightCell(k) = heightCell(k + 1) + layerThickness(k, iCell) + end do + + vertVelocityTop(:, iCell) = 0.0_RKIND + + do k = minLevelCell(iCell), maxLevelCell(iCell) + ! compute vertical velocity at the top of each Cell + vertVelocityTop(k,iCell) = (heightCellLag(k) - heightCell(k))/dt + end do + end do + !$omp end do + !$omp end parallel + + deallocate(heightCell) + deallocate(heightCellLag) + + end subroutine ocn_diagnostic_solve_vertVel_remap!}}} + !*********************************************************************** ! ! routine ocn_diagnostic_solve_surface_pressure @@ -2209,7 +2335,8 @@ subroutine ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, oldLayerT err = 0 - if (config_vert_coord_movement.eq.'impermeable_interfaces') then + if (config_vert_coord_movement == 'impermeable_interfaces' .or. & + configVertAdvMethod == vertAdvRemap) then vertAleTransportTop=0.0_RKIND !$acc update device(vertAleTransportTop) return diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection.F index 8efe2433c048..b2fe794421af 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection.F @@ -143,7 +143,14 @@ subroutine ocn_tracer_advection_init(err)!{{{ ! set some basic flags for options tracerAdvOff = config_disable_tr_adv - monotonicOn = config_monotonic + if (config_flux_limiter == 'monotonic') then + monotonicOn = .true. + elseif (config_flux_limiter == 'none') then + monotonicOn = .false. + else + CALL mpas_log_write('Advection limiter not supported', & + MPAS_LOG_CRIT) + endif ! set all other options from submodule initialization routines call ocn_tracer_advection_shared_init(err1) diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_vert.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_vert.F index 2bb1eb147730..acf90b4cd5db 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_vert.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_advection_vert.F @@ -297,7 +297,7 @@ subroutine ocn_tracer_advection_vert_init(err) !{{{ err = 0 ! set error code to success ! set choice of vertical advection order - select case (config_vert_tracer_adv_order) + select case (config_vert_tracer_adv_flux_order) case (2) vertOrder = vertOrder2 coef3rdOrder = 0.0_RKIND diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_advection.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_advection.F new file mode 100644 index 000000000000..a390d6529ef5 --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_advection.F @@ -0,0 +1,98 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_vertical_advection +! +!> \brief MPAS ocean vertical advection module +!> \author Carolyn Begeman +!> \date February 2022 +!> \details +!> This module contains the initialization for vertical advection +!> schemes. +! +!----------------------------------------------------------------------- + +module ocn_vertical_advection + + use ocn_config + use mpas_log + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + integer, public :: & + configVertAdvMethod ! choice of vertical advection method + integer, parameter, public :: &! supported vertical advection methods + vertAdvFluxForm = 1, &! flux form + vertAdvRemap = 2 ! remapping + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + public :: ocn_vertical_advection_init + +contains + +!*********************************************************************** + +!*********************************************************************** +! +! routine ocn_vertical_advection_init +! +!> \brief Initializes vertical advection scheme properties +!> \author Carolyn Begeman +!> \date February 2022 +!> \details +!> This routine initializes quantities related to +!> the vertical advection scheme. +! +!---------------------------------------------------------------------- + + subroutine ocn_vertical_advection_init(err)!{{{ + + !-------------------------------------------------------------------- + + integer, intent(out) :: err + + + err = 0 + + select case (trim(config_vert_advection_method)) + + case ('flux-form') + + configVertAdvMethod = vertAdvFluxForm + + case ('remap') + + configVertAdvMethod = vertAdvRemap + + case default + + call mpas_log_write( & + 'Invalid choice for config_vert_advection_method. Choices are: flux-form, remap') + err = 1 + + end select + + end subroutine ocn_vertical_advection_init!}}} + +!*********************************************************************** + +end module ocn_vertical_advection + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_regrid.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_regrid.F new file mode 100644 index 000000000000..369a7919f31e --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_regrid.F @@ -0,0 +1,131 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_vertical_regrid +! +!> \brief MPAS ocean vertical regridding +!> \author Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis +!> \date July 2021 +!> \details +!> This module contains the vertical regridding routine, used for +!> vertical Lagrangian remapping. +! +!----------------------------------------------------------------------- + +module ocn_vertical_regrid + + use mpas_derived_types + use mpas_pool_routines + use mpas_constants + use mpas_dmpar + use mpas_threading + use mpas_vector_reconstruction + use mpas_spline_interpolation + use mpas_timer + use mpas_log + + use ocn_constants + use ocn_config + use ocn_mesh + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_vert_regrid + public :: ocn_vert_regrid_init + + contains + +!*********************************************************************** +! +! routine ocn_vert_regrid +! +!> \brief Regridding for vertical Lagrangian remapping +!> \author Carolyn Begeman +!> \date 2021 +!> \details +!> This routine determines the layerThickness to remap to subject to +!> constraints +! +!----------------------------------------------------------------------- + + subroutine ocn_vert_regrid(restingThickness, layerThicknessLag, & + layerThicknessTarget, err) + + real (kind=RKIND), dimension(:, :), intent(in) :: & + layerThicknessLag, & ! layerThickness after the lagrangian step + restingThickness + real (kind=RKIND), dimension(:, :), intent(out) :: & + layerThicknessTarget ! adjusted target layerThickness for remapping + integer, intent(out) :: err !< Output: Error flag + + integer :: iCell, k, kmin, kmax + + real (kind=RKIND) :: totalThickness + + err = 0 + + ! Calculate z-star layer locations + layerThicknessTarget = 0.0_RKIND + + !$omp parallel + !$omp do schedule(runtime) private(k,totalThickness) + do iCell = 1, nCellsAll + totalThickness = 0.0_RKIND + do k = maxLevelCell(iCell), minLevelCell(iCell), -1 + totalThickness = totalThickness + layerThicknessLag(k,iCell) + end do + do k = maxLevelCell(iCell), minLevelCell(iCell), -1 + layerThicknessTarget(k, iCell) = restingThickness(k, iCell) * & + totalThickness / bottomDepth(iCell) + end do + end do + !$omp end do + !$omp end parallel + + + end subroutine ocn_vert_regrid + + +!*********************************************************************** +! +! routine ocn_vert_regrid_init +! +!> \brief Initializes ocean vertical regridding +!> \author Carolyn Begeman +!> \date July 2021 +!> \details +!> This routine initializes parameters required for vertical Lagrangian +!> regridding +! +!----------------------------------------------------------------------- + + subroutine ocn_vert_regrid_init(err) + + integer, intent(out) :: err !< Output: Error flag + + err = 0 + + end subroutine ocn_vert_regrid_init + +end module ocn_vertical_regrid +! vim: foldmethod=marker diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F new file mode 100644 index 000000000000..7bde91d922b0 --- /dev/null +++ b/components/mpas-ocean/src/shared/mpas_ocn_vertical_remap.F @@ -0,0 +1,457 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_vertical_remap +! +!> \brief MPAS ocean vertical Lagrangian remapping +!> \author Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis +!> \date July 2021 +!> \details +!> This module contains the vertical remapping routine. +! +!----------------------------------------------------------------------- + +module ocn_vertical_remap + + use mpas_derived_types + use mpas_pool_routines + use mpas_constants + use mpas_dmpar + use mpas_threading + use mpas_vector_reconstruction + use mpas_spline_interpolation + use mpas_timer + + use ocn_constants + use ocn_config + use ocn_diagnostics_variables + use ocn_mesh + use ocn_vertical_regrid + + use ppr_1d + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + + type(rmap_opts) :: opts + integer :: bc_upper, bc_lower + integer :: itimestepLastRemap + + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_remap_vert_state + public :: ocn_vertical_remap_init + + contains + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_remap_vert_state +! +!> \brief MPAS ocean vertical Lagrangian remapping +!> \author Carolyn Begeman, Darren Engwirda, Xylar Asay-Davis +!> \date July 2021 +!> \details +!> This routine remaps state variables from current layerThickness to +!> a new layerThickness +! +!----------------------------------------------------------------------- + + subroutine ocn_remap_vert_state(block, err) + + type (block_type), intent(in) :: block + integer, intent(out) :: err + + integer :: nCells, nEdges + integer :: nLayers, nLevels, nVars, nDoFs, nTracers + integer :: cell1, cell2, iCell, jCell, iEdge, iTrac + integer :: k, kmin, kmax, kTop, kBot + + type (mpas_pool_type), pointer :: verticalMeshPool + type (mpas_pool_type), pointer :: statePool + type (mpas_pool_type), pointer :: tracersPool + type (mpas_pool_iterator_type) :: groupItr + + type(rmap_work) :: work ! Internal workspace for PPR + type(rcon_ends), dimension(:), allocatable :: bcUpper, bcLower + + real (kind=RKIND), dimension(:, :, :), allocatable :: dataNow, dataNew + real (kind=RKIND), dimension(:, :), allocatable :: & + layerThicknessNew, & ! layerThickness (new target grid) + layerThickEdgeNew, & ! layerThickness at edges (new target grid) + heightCellNow, heightEdgeNow, & ! depth at layer interfaces (lagrangian grid) + heightCellNew, heightEdgeNew ! depth at layer interfaces (new target grid) + + real (kind=RKIND), dimension(:, :), pointer :: layerThickness, restingThickness + real (kind=RKIND), dimension(:, :), pointer :: normalVelocity + real (kind=RKIND), dimension(:, :), pointer :: highFreqThickness,lowFreqDivergence + real (kind=RKIND), dimension(:, :, :), pointer :: tracersGroup + + err = 0 + + itimestepLastRemap = itimestepLastRemap + 1 + if (itimestepLastRemap < config_vert_remap_interval) return + itimestepLastRemap = 0 + + ! Remapping currently only supports one block for the whole domain rather + ! than sub-blocks + call mpas_pool_get_subpool(block % structs, 'verticalMesh', verticalMeshPool) + call mpas_pool_get_subpool(block % structs, 'state', statePool) + + call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) + + call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness) + + call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 2) + call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, 2) + if (config_use_freq_filtered_thickness) then + call mpas_pool_get_array(statePool, 'highFreqThickness', highFreqThickness, 2) + call mpas_pool_get_array(statePool, 'lowFreqDivergence', lowFreqDivergence, 2) + end if + + nCells = nCellsAll + nEdges = nEdgesAll + + ! SETUP VERTICAL GRIDS: + + allocate(heightCellNow(nVertLevels + 1, nCells)) + allocate(heightCellNew(nVertLevels + 1, nCells)) + allocate(heightEdgeNow(nVertLevels + 1, nEdges)) + allocate(heightEdgeNew(nVertLevels + 1, nEdges)) + allocate(layerThicknessNew(nVertLevels, nCells + 1)) + allocate(layerThickEdgeNew(nVertLevels, nEdges)) + + ! Compute new layer thicknesses based on vertical coordinate choice + ! Currently, zstar coordinate is hard-coded + call ocn_vert_regrid(restingThickness, layerThickness, & + layerThicknessNew, err) + + ! Compute the layer interface locations between bottomDepth and -ssh + ! (rather than between -bottomDepth and ssh due to PPR limitations) + + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCells + heightCellNow(:, iCell) = bottomDepth(iCell) + heightCellNew(:, iCell) = bottomDepth(iCell) + + ! reconstruct "now" heights + do k = maxLevelCell(iCell), minLevelCell(iCell), -1 + heightCellNow(k, iCell) = & + heightCellNow(k + 1, iCell) - layerThickness(k, iCell) + end do + + ! reconstruct "new" heights + heightCellNew(minLevelCell(iCell), iCell) = heightCellNow(minLevelCell(iCell), iCell) + do k = maxLevelCell(iCell), minLevelCell(iCell)+1, -1 + heightCellNew(k, iCell) = & + heightCellNew(k + 1, iCell) - layerThicknessNew(k, iCell) + end do + end do + !$omp end do + !$omp end parallel + +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(layerThickness, layerThickEdgeMean, & + !$acc minLevelEdgeBot, maxLevelEdgeTop, cellsOnEdge) & + !$acc private(k, kmin, kmax, cell1, cell2) +#else + !$omp parallel + !$omp do schedule(runtime) private(k, kmin, kmax, cell1, cell2) +#endif + do iEdge = 1, nEdgesAll + kmin = minLevelEdgeBot(iEdge) + kmax = maxLevelEdgeTop(iEdge) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1, nVertLevels + ! initialize layerThicknessEdgeMean to avoid divide by + ! zero and NaN problems. + layerThickEdgeNew(k,iEdge) = -1.0e34_RKIND + end do + do k = kmin, kmax + ! central differenced + layerThickEdgeNew(k,iEdge) = 0.5_RKIND * & + (layerThicknessNew(k,cell1) + & + layerThicknessNew(k,cell2)) + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + !$omp parallel + !$omp do schedule(runtime) private(k, kTop, kBot, cell1, cell2) + do iEdge = 1, nEdges + + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + kTop = minLevelEdgeBot(iEdge) + kBot = maxLevelEdgeTop(iEdge) + + ! We only need to define the kBot+1 depth for one active cell + ! as the absolute kBot+1 depth has no effect on remapping + if (cell1 <= nCellsAll .and. maxLevelCell(cell1) > 0) then + heightEdgeNow(:, iEdge) = heightCellNow(kBot+1, cell1) + heightEdgeNew(:, iEdge) = heightCellNew(kBot+1, cell1) + elseif (cell2 <= nCellsAll .and. maxLevelCell(cell2) > 0) then + heightEdgeNow(:, iEdge) = heightCellNow(kBot+1, cell2) + heightEdgeNew(:, iEdge) = heightCellNew(kBot+1, cell2) + else + call mpas_log_write('No valid cell for edge remapping', MPAS_LOG_WARN) + heightEdgeNow(:,iEdge) = 0.0_RKIND + cycle + endif + + do k = kBot, kTop, -1 + heightEdgeNow(k, iEdge) = heightEdgeNow(k+1, iEdge) - layerThickEdgeMean(k, iEdge) + heightEdgeNew(k, iEdge) = heightEdgeNew(k+1, iEdge) - layerThickEdgeNew(k, iEdge) + end do + + end do + !$omp end do + !$omp end parallel + + ! Assign new layerThicknesses + layerThickness = layerThicknessNew + + !---------------------------------------------------------------------------- + ! ACTUAL REMAPPING FROM HERE + + nDoFs = 1 ! always 1 for MPAS, it'd be >1 for DG methods (with more DoF per layer) + + ! how many variables to be remapped concurrently per column + if (config_use_freq_filtered_thickness) then + nVars = 3 + else + nVars = 1 + endif + + allocate(dataNow(nDoFs, nVars, nVertLevels+1)) + allocate(dataNew(nDoFs, nVars, nVertLevels+1)) + allocate(bcUpper(nVars)) + allocate(bcLower(nVars)) + + ! set boundary conditions + bcUpper%bcopt = bc_upper + bcLower%bcopt = bc_lower + + dataNow = 0.0_RKIND + dataNew = 0.0_RKIND + + call work%init(nVertLevels + 1, nVars, opts) + + ! Remap all edge-centred variables from heightEdgeNow to heightEdgeNew + if (config_use_freq_filtered_thickness) then + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(kTop, kBot, nLayers, nLevels, dataNow, dataNew, work) + do iEdge = 1, nEdges + + kTop = minLevelEdgeBot(iEdge) + kBot = maxLevelEdgeTop(iEdge) + + nLayers = kBot - kTop + 1 + nLevels = nLayers + 1 + + if (nLevels .lt. 2) cycle + + dataNow(1, 1, 1:nLayers) = normalVelocity (kTop:kBot, iEdge) + dataNow(1, 2, 1:nLayers) = highFreqThickness(kTop:kBot, iEdge) + dataNow(1, 3, 1:nLayers) = lowFreqDivergence(kTop:kBot, iEdge) + + call rmap1d(nLevels, nLevels, nVars ,nDoFs, & + heightEdgeNow(kTop:kBot+1, iEdge), & + heightEdgeNew(kTop:kBot+1, iEdge), & + dataNow, dataNew, & + bcUpper, bcLower, work, opts) + + normalVelocity (kTop:kBot, iEdge) = dataNew(1, 1, 1:nLayers) + highFreqThickness(kTop:kBot, iEdge) = dataNew(1, 2, 1:nLayers) + lowFreqDivergence(kTop:kBot, iEdge) = dataNew(1, 3, 1:nLayers) + end do + !$omp end do + !$omp end parallel + else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(kTop, kBot, nLayers, nLevels, dataNow, dataNew, work) + do iEdge = 1, nEdges + + kTop = minLevelEdgeBot(iEdge) + kBot = maxLevelEdgeTop(iEdge) + + nLayers = kBot - kTop + 1 + nLevels = nLayers + 1 + + if (nLevels .lt. 2) cycle + + dataNow(1, 1, 1:nLayers) = normalVelocity(kTop:kBot, iEdge) + + call rmap1d(nLevels, nLevels, nVars, nDoFs, & + heightEdgeNow(kTop:kBot+1, iEdge), & + heightEdgeNew(kTop:kBot+1, iEdge), & + dataNow, dataNew, & + bcUpper, bcLower, work, opts) + + normalVelocity (kTop:kBot, iEdge) = dataNew(1, 1, 1:nLayers) + + end do + !$omp end do + !$omp end parallel + end if + deallocate(dataNow,dataNew) + deallocate(bcUpper,bcLower) + call work%free() + + ! Remap all cell-centred variables from heightCellNow to heightCellNew + + call mpas_pool_begin_iteration(tracersPool) + do while ( mpas_pool_get_next_member(tracersPool, groupItr) ) + if ( groupItr % memberType == MPAS_POOL_FIELD ) then + call mpas_pool_get_array(tracersPool, groupItr % memberName, tracersGroup, 2) + if ( associated(tracersGroup) ) then + + nTracers = size(tracersGroup, dim=1) + + allocate(dataNow(nDoFs, nTracers, nVertLevels)) + allocate(dataNew(nDoFs, nTracers, nVertLevels)) + + call work%init(nVertLevels + 1, nTracers, opts) + + allocate(bcUpper(nTracers)) + allocate(bcLower(nTracers)) + bcUpper%bcopt = bc_upper + bcLower%bcopt = bc_lower + + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(iTrac, kTop, kBot, nLayers, nLevels, dataNow, dataNew, work) + do iCell = 1, nCells + + kTop = minLevelCell(iCell) + kBot = maxLevelCell(iCell) + + nLayers = kBot - kTop + 1 + nLevels = nLayers + 1 + + if (nLevels .lt. 2) cycle + + do iTrac = 1, nTracers + + dataNow(1, iTrac, 1:nLayers) = tracersGroup(iTrac, kTop:kBot, iCell) + + end do + call rmap1d(nLevels, nLevels, nTracers,nDoFs, & + heightCellNow(kTop:kBot+1, iCell), & + heightCellNew(kTop:kBot+1, iCell), & + dataNow, dataNew, & + bcUpper, bcLower, work, opts) + + do iTrac = 1, nTracers + tracersGroup(iTrac, kTop:kBot, iCell) = dataNew(1, iTrac, 1:nLayers) + end do + + end do + !$omp end do + !$omp end parallel + + deallocate(dataNow,dataNew) + deallocate(bcUpper,bcLower) + call work%free() + + end if + end if + end do + + deallocate(heightCellNow) + deallocate(heightEdgeNow) + deallocate(heightCellNew) + deallocate(heightEdgeNew) + + end subroutine ocn_remap_vert_state + +!*********************************************************************** +! +! routine ocn_vertical_remap_init +! +!> \brief Initializes ocean vertical remapping +!> \author Carolyn Begeman +!> \date July 2021 +!> \details +!> This routine initializes parameters required for vertical Lagrangian +!> remapping +! +!----------------------------------------------------------------------- + + subroutine ocn_vertical_remap_init(err) + + integer, intent(out) :: err !< Output: Error flag + + err = 0 + + call ocn_vert_regrid_init(err) + + if ( config_time_integrator == 'RK4' ) then + CALL mpas_log_write('Vertical remap not supported for RK4', & + MPAS_LOG_CRIT) + endif + + ! Options for remapping + if ( config_vert_remap_order == 1 ) then + opts%cell_meth = pcm_method ! PCM method + elseif ( config_vert_remap_order == 2 ) then + opts%cell_meth = plm_method ! PLM method + elseif ( config_vert_remap_order == 3 ) then + opts%cell_meth = ppm_method ! PPM method + opts%edge_meth = p3e_method ! 3rd-order edge interp. + elseif ( config_vert_remap_order == 5 ) then + opts%cell_meth = pqm_method ! PPM method + opts%edge_meth = p5e_method ! 3rd-order edge interp. + else + CALL mpas_log_write('Vertical remap order $i not supported', & + MPAS_LOG_CRIT, intArgs=(/config_vert_remap_order/)) + endif + + if ( config_remap_limiter == 'monotonic' ) then + opts%cell_lims = mono_limit ! monotone slope limits + elseif ( config_remap_limiter == 'weno' ) then + opts%cell_lims = weno_limit ! WENO slope limits + opts%wall_lims = mono_limit ! monotone slope limits + elseif ( config_remap_limiter == 'none' ) then + opts%cell_lims = null_limit ! no slope limits + else + CALL mpas_log_write('Remap limiter not supported', & + MPAS_LOG_CRIT) + endif + + bc_upper = bcon_loose + bc_lower = bcon_loose + + itimestepLastRemap = 0 + + end subroutine ocn_vertical_remap_init + +end module ocn_vertical_remap +! vim: foldmethod=marker