Addition of TRiSK++ to ocean dynamical core

components/mpas-ocean/bld/build-namelist

@@ -511,6 +511,7 @@ add_default($nl, 'config_transport_tests_flow_id');
 
 add_default($nl, 'config_dt');
 add_default($nl, 'config_time_integrator');
+add_default($nl, 'config_triskCV');
 
 ########################
 # Namelist group: hmix #

components/mpas-ocean/bld/build-namelist-section

@@ -54,6 +54,7 @@ add_default($nl, 'config_transport_tests_flow_id');
 
 add_default($nl, 'config_dt');
 add_default($nl, 'config_time_integrator');
+add_default($nl, 'config_triskCV');
 
 ########################
 # Namelist group: hmix #

components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml

@@ -53,6 +53,7 @@
 <config_dt ocn_grid="SOwISC12to60E2r4">'00:10:00'</config_dt>
 <config_dt ocn_grid="ECwISC30to60E2r1">'00:30:00'</config_dt>
 <config_time_integrator>'split_explicit'</config_time_integrator>
+<config_triskCV>.true.</config_triskCV>
 
 <!-- hmix -->
 <config_hmix_scaleWithMesh ocn_grid="oQU480">.false.</config_hmix_scaleWithMesh>

components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml

@@ -235,6 +235,13 @@ Valid values: 'split_explicit', 'RK4', 'unsplit_explicit', 'split_implicit'
 Default: Defined in namelist_defaults.xml
 </entry>
 
+<entry id="config_triskCV" type="logical"
+	category="time_integration" group="time_integration">
+Turns on the computation of the kinetic energy in the TRiSK-CV way.
+
+Valid values: .true. or .false.
+Default: Defined in namelist_defaults.xml
+</entry>
 
 <!-- hmix -->
 

components/mpas-ocean/src/Registry.xml

@@ -208,6 +208,10 @@
 					description="Time integration method."
 					possible_values="'split_explicit', 'RK4', 'unsplit_explicit', 'split_implicit'"
 		/>
+		<nml_option name="config_triskCV" type="logical" default_value=".true." units="unitless"
+                                        description="Turns on the computation of the kinetic energy in the TRiSK-CV way."
+                                        possible_values=".true. or .false."
+		/>
 	</nml_record>
 	<nml_record name="hmix" mode="forward">
 		<nml_option name="config_hmix_scaleWithMesh" type="logical" default_value=".false." units="unitless"
@@ -271,6 +275,28 @@
 					possible_values="any positive real"
 		/>
 	</nml_record>
+        <nml_record name="hmix_aust" mode="forward">
+                <nml_option name="config_aust_scale_factor" type="real" default_value="0.0" units="unitless"
+                                        description="Anticipated potential vorticity (APV) method scale factor, $c_{apv}$. When zero, APV is off."
+                                        possible_values="Any non-negative number, typically between zero and one."
+                />
+                <nml_option name="config_use_mom_div2" type="logical" default_value=".false." units="unitless"
+                                        description="If true, AUST horizontal mixing is used on the momentum equation."
+                                        possible_values=".true. or .false."
+                />
+                <nml_option name="config_mom_div2" type="real" default_value="1.0e5" units="m^2 s^{-1}"
+                                        description="Horizontal viscosity, $\nu_{2h}$. If config_hmix_use_ref_cell_width = .true. then $\nu_h$ = config_mom_del2*(cellWidth / config_hmix_use_ref_cell_width)."
+                                        possible_values="any positive real"
+                />
+                <nml_option name="config_use_mom_div4" type="logical" default_value=".false." units="unitless"
+                                        description="If true, biharmonic horizontal mixing is used on the momentum equation."
+                                        possible_values=".true. or .false."
+                />
+                <nml_option name="config_mom_div4" type="real" default_value="1.0e12" units="m^4 s^{-1}"
+                                        description="Coefficient for horizontal biharmonic operator on momentum.  If config_hmix_use_ref_cell_width = .true. then $\nu_{4h}$ = config_mom_del4*(cellWidth / config_hmix_use_ref_cell_width)^3."
+                                        possible_values="any positive real"
+                />
+        </nml_record>
 	<nml_record name="hmix_Leith" mode="forward">
 		<nml_option name="config_use_Leith_del2" type="logical" default_value=".false." units="unitless"
 					description="If true, the Leith enstrophy-cascade closure is turned on"
@@ -1802,8 +1828,12 @@
 			<var name="refLayerThickness"/>
 			<var name="xtime"/>
 			<var name="zMid"/>
-			<var name="zTop"/>
+			<var name="zTop"/>  
+                        <var name="weightsOnEdge_lsqr"/>
 			<var name="kineticEnergyCell"/>
+                        <var name="kineticEnergyEdge"/>
+                        <var name="u_perp"/>
+                        <var name="divergence"/>
 			<var name="relativeVorticityCell"/>
 			<var name="areaCellGlobal"/>
 			<var name="areaEdgeGlobal"/>
@@ -2260,6 +2290,9 @@
 		<var name="weightsOnEdge" type="real" dimensions="maxEdges2 nEdges" units="unitless"
 			 description="Reconstruction weights associated with each of the edgesOnEdge, used to reconstruct the tangentialVelocity from normalVelocities on neighboring edges."
 		/>
+                <var name="weightsOnEdge_lsqr" type="real" dimensions="maxEdges2 nEdges" units="unitless"
+			 desription="Weights associated with each of the edgesOnEdge needed for the tangential velocity in the kinetic energy for TRiSK-CV."
+                />   
 		<var name="dvEdge" type="real" dimensions="nEdges" units="m"
 			 description="Length of each edge, computed as the distance between verticesOnEdge."
 		/>
@@ -2304,8 +2337,8 @@
 		/>
 		<var name="kiteAreasOnVertex" type="real" dimensions="vertexDegree nVertices" units="m^2"
 			 description="Area of the portions of each dual cell that are part of each cellsOnVertex."
-		/>
-		<var name="fEdge" type="real" dimensions="nEdges" units="radians s^{-1}"
+		/> 
+		<var name="fEdge" type="real" dimensions="nEdges" units="s^{-1}"
 			 description="Coriolis parameter at edges."
 		/>
 		<var name="fVertex" type="real" dimensions="nVertices" units="radians s^{-1}"
@@ -2656,6 +2689,14 @@
 			 description="kinetic energy of horizontal velocity on cells"
 			 packages="forwardMode;analysisMode"
 		/>
+                <var name="kineticEnergyEdge" type="real" dimensions="nVertLevels nEdges Time" units="m^2 s^{-2}"
+                         description="kinetic energy of horizontal velocity on edges"
+                         packages="forwardMode;analysisMode"
+                />
+                <var name="u_perp" type="real" dimensions="nVertLevels nEdges Time" units="m^2 s^{-2}"
+                         description="tangential velocity on edges needed for the computation of kineticEnergyEdge in TRSK++"
+                         packages="forwardMode;analysisMode"
+                />
 		<var name="viscosity" type="real" dimensions="nVertLevels nEdges Time" units="m^2 s^{-1}"
 			 description="horizontal viscosity"
 			 packages="forwardMode;analysisMode"

components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F

@@ -89,6 +89,8 @@ module ocn_forward_mode
    use ocn_constants
    use ocn_config
 
+   use ocn_ke_triskcv
+
    implicit none
    private
 
@@ -171,6 +173,8 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{
       !-----------------------------------------------------------------
       ierr = 0
 
+      call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool)
+
       ! Setup ocean config pool
       call ocn_config_init(domain % configs)
 
@@ -301,6 +305,12 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{
 
       call ocn_init_routines_block(block, dt, ierr)
 
+      if (config_triskCV) then 
+         ! Compute weightsOnEdge_lsqr needed in TRiSK-CV for the computation of
+         ! the tangential velocity used in the new kinetic energy
+         call Wlsqr(domain)
+      end if
+
       if(ierr.eq.1) then
           call mpas_log_write('An error was encountered in ocn_init_routines_block', MPAS_LOG_CRIT)
       endif

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F

@@ -703,7 +703,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
 
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          ! BEGIN baroclinic iterations on linear Coriolis term
-         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
 
          do j=1,numClinicIterations(splitExplicitStep)
 
@@ -785,15 +785,15 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                               'normalBaroclinicVelocity', timeLevel=2)
             call mpas_timer_stop("se halo normalBaroclinicVelocity")
 
-            call mpas_timer_stop('bcl iters on linear Coriolis')
+            call mpas_timer_stop('bcl iters on linear Coriolis') 
 
          end do  ! do j=1,config_n_bcl_iter
 
          call mpas_timer_start('se halo barotropicForcing')
          call mpas_dmpar_field_halo_exch(domain, 'barotropicForcing')
          call mpas_timer_stop('se halo barotropicForcing')
 
-         call mpas_timer_stop("se bcl vel")
+         call mpas_timer_stop("se bcl vel")       
 
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          ! END baroclinic iterations on linear Coriolis term
@@ -1890,7 +1890,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
          ! time n+1
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-         elseif (splitExplicitStep == numTSIterations) then
+         elseif (splitExplicitStep == numTSIterations) then 
 
             !$omp parallel
             !$omp do schedule(runtime) private(k)
@@ -2107,6 +2107,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
 
             !$omp parallel
             !$omp do schedule(runtime) private(k)
+
             do iEdge = 1, nEdgesAll
             do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
                normalVelocityNew(k,iEdge) = &
@@ -2133,6 +2134,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
 
       call mpas_timer_start("se implicit vert mix")
 
+      !call mpas_dmpar_field_halo_exch(domain, 'normalVelocity', timeLevel=2)
+
       ! Call ocean diagnostic solve in preparation for vertical mixing.
       ! Note it is called again after vertical mixing, because u and
       ! tracers change. For Richardson vertical mixing, only density,
@@ -2198,7 +2201,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
       !$acc exit data delete(layerThicknessNew, normalVelocityNew)
 #endif
 
-      call mpas_dmpar_field_halo_exch(domain, 'surfaceFrictionVelocity')
+      call mpas_dmpar_field_halo_exch(domain, 'surfaceFrictionVelocity') 
 
       ! Compute normalGMBolusVelocity; it will be added to the
       ! baroclinic modes in Stage 2 above.
@@ -2319,7 +2322,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
       end if
 
       call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, &
-                                scratchPool, tracersPool, 2)
+                                scratchPool, tracersPool, 2) 
 
       ! Update the effective desnity in land ice if we're coupling to
       ! land ice

components/mpas-ocean/src/ocean.cmake

@@ -99,6 +99,7 @@ list(APPEND RAW_SOURCES
   core_ocean/shared/mpas_ocn_time_varying_forcing.F
   core_ocean/shared/mpas_ocn_wetting_drying.F
   core_ocean/shared/mpas_ocn_vel_tidal_potential.F
+  core_ocean/shared/mpas_ocn_ke_triskCV.F
 )
 
 set(OCEAN_DRIVER

components/mpas-ocean/src/shared/Makefile

@@ -69,15 +69,16 @@ OBJS = mpas_ocn_init_routines.o \
 	   mpas_ocn_vel_tidal_potential.o \
 	   mpas_ocn_vel_forcing_topographic_wave_drag.o \
 	   mpas_ocn_transport_tests.o \
-	   mpas_ocn_vel_self_attraction_loading.o 
+	   mpas_ocn_vel_self_attraction_loading.o  \
+           mpas_ocn_ke_triskCV.o
 
 all: $(OBJS)
 
-mpas_ocn_init_routines.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics.o mpas_ocn_diagnostics_variables.o mpas_ocn_gm.o mpas_ocn_forcing.o mpas_ocn_surface_land_ice_fluxes.o
+mpas_ocn_init_routines.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics.o mpas_ocn_diagnostics_variables.o mpas_ocn_gm.o mpas_ocn_forcing.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_ke_triskCV.o
 
 mpas_ocn_tendency.o: mpas_ocn_high_freq_thickness_hmix_del2.o mpas_ocn_tracer_surface_restoring.o mpas_ocn_thick_surface_flux.o mpas_ocn_tracer_short_wave_absorption.o mpas_ocn_tracer_advection.o mpas_ocn_tracer_hmix.o mpas_ocn_tracer_nonlocalflux.o mpas_ocn_surface_bulk_forcing.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_tracer_surface_flux_to_tend.o mpas_ocn_tracer_interior_restoring.o mpas_ocn_tracer_exponential_decay.o mpas_ocn_tracer_ideal_age.o mpas_ocn_tracer_TTD.o mpas_ocn_vmix.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_frazil_forcing.o mpas_ocn_tidal_forcing.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_DMS.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_tracer_CFC.o mpas_ocn_diagnostics.o mpas_ocn_wetting_drying.o mpas_ocn_vel_self_attraction_loading.o  mpas_ocn_vel_tidal_potential.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_thick_hadv.o mpas_ocn_thick_vadv.o mpas_ocn_vel_hadv_coriolis.o mpas_ocn_vel_pressure_grad.o mpas_ocn_vel_vadv.o mpas_ocn_vel_hmix.o mpas_ocn_vel_forcing.o
 
-mpas_ocn_diagnostics.o: mpas_ocn_thick_ale.o mpas_ocn_equation_of_state.o mpas_ocn_gm.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_surface_land_ice_fluxes.o
+mpas_ocn_diagnostics.o: mpas_ocn_thick_ale.o mpas_ocn_equation_of_state.o mpas_ocn_gm.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_diagnostics_variables.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_ke_triskCV.o
 
 mpas_ocn_diagnostics_variables.o: mpas_ocn_config.o
 
@@ -217,6 +218,8 @@ mpas_ocn_vel_pressure_grad.o: mpas_ocn_constants.o
 
 mpas_ocn_vel_tidal_potential.o: mpas_ocn_diagnostics_variables.o
 
+mpas_ocn_ke_triskCV.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_diagnostics_variables.o mpas_ocn_mesh.o mpas_ocn_thick_ale.o
+
 clean:
 	$(RM) *.o *.i *.mod *.f90