Add GM option based on prognostic MEKE

components/mpas-ocean/meke_histograms.py

@@ -0,0 +1,106 @@
+import numpy
+import xarray
+import matplotlib.pyplot as plt 
+
+comparison=False
+if comparison:
+    master_path = '/compyfs/bege567/mpaso_scratch/compass-v1.2.0-E3SM-master/ocean/global_ocean/QU240/PHC/performance_test'
+#meke_path = '/compyfs/bege567/mpaso_scratch/compass-v1.2.0-E3SM-meke-off/ocean/global_ocean/QU240/PHC/performance_test'
+#meke_path = '/compyfs/bege567/mpaso_scratch/compass-v1.2.0-E3SM-meke/ocean/global_ocean/QU240/PHC/performance_test'
+meke_path = '/compyfs/bege567/mpaso_scratch/compass-v1.2.0-E3SM-meke-rk4/ocean/baroclinic_channel/10km/rpe_test/'
+#meke_path = '/compyfs/bege567/mpaso_scratch/compass-v1.2.0-E3SM-meke-rk4/ocean/baroclinic_channel/10km/default'
+#subdir='forward'
+subdir='rpe_test_5_nu_200'
+if comparison:
+    master_data = xarray.open_dataset(f'{master_path}/{subdir}/output.nc')
+meke_data = xarray.open_dataset(f'{meke_path}/{subdir}/output.nc')
+
+nt = meke_data.sizes['Time']
+
+if comparison:
+    gmBolusKappa1 = master_data.gmBolusKappa
+    gmCellsAvg1 = numpy.mean(gmBolusKappa1, axis=1) # nt
+    gmCellsMin1 = numpy.amin(gmBolusKappa1, axis=1) # nt
+    gmCellsMax1 = numpy.amax(gmBolusKappa1, axis=1) # nt
+
+gmBolusKappa = meke_data.gmBolusKappa
+gmCellsAvg = numpy.mean(gmBolusKappa, axis=1) # nt
+gmCellsMin = numpy.amin(gmBolusKappa, axis=1) # nt
+gmCellsMax = numpy.amax(gmBolusKappa, axis=1) # nt
+
+
+mekeSource = meke_data.mekeSourceTopOfEdge.values
+mekeSink = meke_data.mekeSink.values
+meke = meke_data.meke.values
+
+mekeCellsAvg = numpy.mean(meke, axis=(1,2)) # nt, nVertLevels
+mekeCellsMin = numpy.amin(meke, axis=(1,2)) # nt
+mekeCellsMax = numpy.amax(meke, axis=(1,2)) # nt
+
+mekeSourceEdgesAvg = numpy.mean(mekeSource, axis=(1,2)) # nt, nVertLevels
+mekeSourceEdgesMin = numpy.amin(mekeSource, axis=(1,2)) # nt
+mekeSourceEdgesMax = numpy.amax(mekeSource, axis=(1,2)) # nt
+
+mekeSinkCellsAvg = numpy.mean(mekeSink, axis=(1,2)) # nt, nVertLevels
+mekeSinkCellsMin = numpy.amin(mekeSink, axis=(1,2)) # nt
+mekeSinkCellsMax = numpy.amax(mekeSink, axis=(1,2)) # nt
+
+tidx = 0
+#tidx = nt-1
+
+print('plotting gm')
+fig = plt.figure()
+if comparison:
+    plt.hist(gmBolusKappa1[tidx,:],label='master',bins=100)
+plt.hist(gmBolusKappa[tidx,:],label='meke',bins=100)
+plt.legend()
+plt.savefig(f'{meke_path}/gm_histogram.png')
+plt.close(fig)
+
+fig,ax = plt.subplots(2)
+print('plotting meke terms')
+ax[0].hist(mekeSource[tidx,:].flatten(),label='source',bins=100)
+ax[0].hist(mekeSink[tidx,:].flatten(),label='sink',bins=100)
+ax[0].set_xlabel('MEKE terms')
+ax[0].set_xscale('log')
+ax[0].legend()
+
+print('plotting meke')
+ax[1].hist(meke[tidx,:].flatten(),bins=100)
+ax[1].set_xlabel('MEKE')
+ax[1].set_xscale('log')
+ax[1].legend()
+fig.savefig(f'{meke_path}/meke_histogram.png')
+plt.close(fig)
+
+fig,ax = plt.subplots(3,1)
+ax[0].plot(mekeSourceEdgesMin, '--b')
+ax[0].plot(mekeSourceEdgesAvg, '-b', label='Source')
+ax[0].plot(mekeSourceEdgesMax, ':b')
+ax[0].plot(mekeSinkCellsMin, '--k')
+ax[0].plot(mekeSinkCellsAvg, '-k', label='Sink')
+ax[0].plot(mekeSinkCellsMax, ':k')
+ax[0].set_yscale('log')
+ax[0].legend()
+ax[1].plot(mekeCellsMin, '--k')
+ax[1].plot(mekeCellsAvg, '-k', label='MEKE')
+ax[1].plot(mekeCellsMax, ':k')
+ax[1].set_yscale('log')
+ax[1].legend()
+ax[2].plot(gmCellsMin, '--k')
+ax[2].plot(gmCellsAvg, '-k', label='GM MEKE')
+ax[2].plot(gmCellsMax, ':k')
+if comparison:
+    ax[2].plot(gmCellsMin1, '--b')
+    ax[2].plot(gmCellsAvg1, '-b', label='GM master')
+    ax[2].plot(gmCellsMax1, ':b')
+    ax[2].legend()
+fig.savefig(f'{meke_path}/meke_timeseries.png')
+plt.close(fig)
+
+print('mekeSource: min, 10th, 50th, 90th, max')
+print(f'{numpy.min(mekeSource[tidx,:])},{numpy.percentile(mekeSource[tidx,:],[10,50,90])},{numpy.max(mekeSource[tidx,:])}')
+print('mekeSink: min, 10th, 50th, 90th, max')
+print(f'{numpy.min(mekeSink[tidx,:])},{numpy.percentile(mekeSink[tidx,:],[10,50,90])},{numpy.max(mekeSink[tidx,:])}')
+print('meke: min, 10th, 50th, 90th, max')
+print(f'{numpy.min(meke[tidx,:])},{numpy.percentile(meke[tidx,:],[10,50,90])},{numpy.max(meke[tidx,:])}')

components/mpas-ocean/src/Registry.xml

@@ -348,6 +348,14 @@
 					description="If true, the standard GM for the tracer advection and mixing is turned on."
 					possible_values=".true. or .false."
 		/>
+		<nml_option name="config_prognostic_meke" type="logical" default_value=".false." units="unitless"
+					description="If true, a prognostic meke is solved for at each time step"
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_use_prognostic_meke" type="logical" default_value=".false." units="unitless"
+					description="If true, the prognostic meke is used to compute gmBolusKappa"
+					possible_values=".true. or .false."
+		/>
 		<nml_option name="config_GM_closure" type="character" default_value="EdenGreatbatch" units="unitless"
 					description="Control what method used to compute GM $\kappa$. Both 'constant' and 'N2_dependent' use the method in Ferrari et al. 2010 (https://doi.org/10.1016/j.ocemod.2010.01.004). 'constant' uses a constant kappa in eqn 16a, while 'N2_dependent' varies kappa in the vertical according to Danabasoglu and Marshall 2007 (https://doi.org/10.1016/j.ocemod.2007.03.006). 'Visbeck' implements a horizontally varying diffusivity of Visbeck et al 1997. EdenGreatbatch implements a simplified form of the EKE scheme in Eden and Greatbatch (2008) Ocean modeling"
 					possible_values="'constant', 'N2_dependent', 'Visbeck', 'EdenGreatbatch'"
@@ -1472,6 +1480,7 @@
 		<package name="splitTimeIntegrator" description="This package includes variables required for either the split or unsplit explicit time integrators."/>
 		<package name="semiImplicitTimePKG" description="This package includes variables required for split-implicit time integrators."/>
 		<package name="thicknessFilter" description="This package includes variables required for frequency filtered thickness."/>
+		<package name="MEKE" description="This package includes variables required for prognostic mesoscale eddy kinetic energy."/>
 		<package name="windStressBulkPKG" description="This package includes varibles required for bulk wind stress forcing."/>
 		<package name="variableBottomDragPKG" description="This package includes varibles required for variable bottom drag."/>
 		<package name="thicknessBulkPKG" description="This package includes varibles required for bulk thickness forcing."/>
@@ -2129,6 +2138,12 @@
 			 packages="thicknessFilter"
 		/>
 
+		<!-- FIELDS FOR MEKE -->
+		<var name="meke" type="real" dimensions="nVertLevels nCells Time" units="m"
+			 description="Mesoscale eddy kinetic energy"
+			 packages="MEKE"
+		/>
+
 		<!-- FIELDS FOR FRAZIL ICE -->
 		<var name="accumulatedFrazilIceMass" type="real" dimensions="nCells Time" units="kg m^{-2}"
 			 description="Mass per unit area of frazil ice produced. Reset to zero at each coupling interval"
@@ -2482,6 +2497,10 @@
 			 description="time tendency of low frequency-filtered divergence"
 			 packages="thicknessFilter"
 		/>
+		<var name="tendMeke" type="real" dimensions="nVertLevels nCells Time" units="XX s^{-1}" name_in_code="meke"
+			 description="time tendency of mesoscale eddy kinetic energy"
+			 packages="MEKE"
+		/>
 	</var_struct>
 	<var_struct name="diagnostics" time_levs="1">
 		<var name="xtime" type="text" dimensions="Time" units="unitless"
@@ -2806,6 +2825,9 @@
 			description="phase speed for the bolus velocity calculation"
 			packages="forwardMode;analysisMode"
 		/>
+		<var name="LrEdge" type="real" dimensions="nEdges" units="m"
+			 description="Length scale for Eden and Greatbatch GM."
+		/>
 		<var name="betaEdge" type="real" dimensions="nEdges Time" units="m^{-1} s^{-1}"
 			description="meridional gradient of the coriolis parameter, used in the mesoscale eddy parameterization schemes"
 		/>
@@ -2853,6 +2875,14 @@
 			 description="gradient Richardson number defined at the edge (horizontally) and top (vertically)"
 			 packages="forwardMode;analysisMode"
 		/>
+		<var name="mekeSourceTopOfEdge" type="real" dimensions="nVertLevelsP1 nEdges Time" units="XX s^{-1}" name_in_code="mekeSourceTopOfEdge"
+			 description="source term tendency of mesoscale eddy kinetic energy at top of edges"
+			 packages="MEKE"
+		/>
+		<var name="mekeSink" type="real" dimensions="nVertLevels nCells Time" units="XX s^{-1}" name_in_code="mekeSink"
+			 description="source term tendency of mesoscale eddy kinetic energy at top of edges"
+			 packages="MEKE"
+		/>
 		<var name="vertViscTopOfEdge" type="real" dimensions="nVertLevelsP1 nEdges Time" units="m^2 s^{-1}"
 			 description="vertical viscosity defined at the edge (horizontally) and top (vertically)"
 			 packages="forwardMode;analysisMode"

components/mpas-ocean/src/driver/mpas_ocn_core_interface.F

@@ -134,6 +134,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
       logical, pointer :: timeVaryingLandIceForcingPKGActive
       logical, pointer :: gotmPKGActive
       logical, pointer :: verticalRemapPKGActive
+      logical, pointer :: MEKEActive
 
       type (mpas_pool_iterator_type) :: pkgItr
       logical, pointer :: packageActive
@@ -166,6 +167,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
       logical, pointer :: config_use_time_varying_atmospheric_forcing
       logical, pointer :: config_use_time_varying_land_ice_forcing
       logical, pointer :: config_use_gotm
+      logical, pointer :: config_prognostic_meke
 
       character (len=StrKIND), pointer :: config_time_integrator
       character (len=StrKIND), pointer :: config_ocean_run_mode
@@ -389,6 +391,16 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
          verticalRemapPKGActive = .true.
       end if
 
+      ! test for use of MEKE
+      !
+      call mpas_pool_get_package(packagePool, 'MEKEActive', MEKEActive)
+      call mpas_pool_get_config(configPool, 'config_prognostic_meke', config_prognostic_meke)
+      if ( forwardModeActive ) then
+         if (config_prognostic_meke) then
+            MEKEActive = .true.
+         end if
+      endif
+
       !
       ! call into analysis member driver to set analysis member packages
       !

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

@@ -307,6 +307,7 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{
       if(ierr.eq.1) then
           call mpas_log_write('An error was encountered in ocn_init_routines_block', MPAS_LOG_CRIT)
       endif
+      call mpas_log_write('Complete ocn_init_routines_block')
 
       xtime = startTimeStamp
 
@@ -409,6 +410,7 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{
               MPAS_LOG_ERR, masterOnly=.true., flushNow=.true.)
          return
       endif
+      call mpas_log_write('Complete vel tend modules')
 
       ! Initialize forcing modules
       call ocn_surface_bulk_forcing_init(err_tmp)
@@ -427,20 +429,27 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{
               MPAS_LOG_ERR, masterOnly=.true., flushNow=.true.)
          return
       endif
+      call mpas_log_write('Complete forcing modules')
 
       ! Initialize tracer tendencies
       call ocn_tracer_hmix_init(err_tmp)
       ierr = ior(ierr, err_tmp)
+      call mpas_log_write('End tracer hmix init')
       call ocn_tracer_hmix_redi_init(domain,err_tmp)
       ierr = ior(ierr, err_tmp)
+      call mpas_log_write('End tracer hmix redi init')
       call ocn_tracer_surface_flux_init(err_tmp)
       ierr = ior(ierr, err_tmp)
+      call mpas_log_write('End tracer surface flux init')
       call ocn_tracer_advection_init(err_tmp)
       ierr = ior(ierr,err_tmp)
+      call mpas_log_write('End tracer advection init')
       call ocn_tracer_short_wave_absorption_init(domain,err_tmp)
       ierr = ior(ierr,err_tmp)
+      call mpas_log_write('End tracer SW init')
       call ocn_gm_init(domain, err_tmp)
       ierr = ior(ierr,err_tmp)
+      call mpas_log_write('End gm init')
       call ocn_tracer_nonlocalflux_init(err_tmp)
       ierr = ior(ierr,err_tmp)
       call ocn_tracer_ecosys_init(domain, err_tmp)
@@ -459,6 +468,7 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{
               MPAS_LOG_ERR, masterOnly=.true., flushNow=.true.)
          return
       endif
+      call mpas_log_write('Complete tracer tend modules')
 
       ! Initialize other modules
       call ocn_vmix_init(domain, err_tmp)
@@ -476,6 +486,7 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{
               MPAS_LOG_ERR, masterOnly=.true., flushNow=.true.)
          return
       endif
+      call mpas_log_write('Complete ocn_forward_mode_init')
 
    !--------------------------------------------------------------------