Merge branch 'jonbob/ocn/add-akt-tilt-coupling' (PR #1278)

Modifications needed for time filtered ssh gradient coupling

This PR adds changes needed to the ocean core to allow filtering of high
frequency modes from sea surface tilt. It includes modifications to the ocn
driver as well as the mpas-o model itself:
* Added a new namelist option specifying the filtering timescale
* Added a new set of variables to registry that are the filtered ssh gradient
* Includes new ICG files for the oEC60to30v3 grid
* Also brings in a new version of the MPAS framework
It incorporates MPAS PR 1245. This PR will change answers for all compsets with
active ocean.

Tested to be BFB when the new namelist option, config_ssh_grad_relax_timescale,
is set to 0.0.  But also tested to have the expected impact in ongoing
T62-oRRS18to6 GMPAS simulation.

Passed automated testing in next

Performed a two year A_WCYCL2000S.ne30_oECv3_ICG test on anvil

Fixes #1160

[NML]
[non-BFB]

components/mpas-cice/cime_config/buildnml

@@ -62,8 +62,8 @@ if ( $ICE_GRID eq 'oEC60to30' ) {
 	$decomp_date .= '161222';
 	$decomp_prefix .= 'mpas-cice.graph.info.';
 } elsif ( $ICE_GRID eq 'oEC60to30v3_ICG' ) {
-        $grid_date .= '170213';
-	$grid_prefix .= 'seaice.EC60to30v3.a16q_smRi_yr102';
+        $grid_date .= '170301';
+	$grid_prefix .= 'seaice.EC60to30v3.restartFrom_anvil0221';
 	$decomp_date .= '161222';
 	$decomp_prefix .= 'mpas-cice.graph.info.';
 } elsif ( $ICE_GRID eq 'oEC60to30wLI' ) {
@@ -276,8 +276,17 @@ if ( -e "$CASEROOT/SourceMods/src.mpascice/$STREAM_NAME" ) {
 	print $stream_file '        type="output"' . "\n";
 	print $stream_file '        filename_template="mpascice.hist.$Y-$M-$D_$S.nc"' . "\n";
 	print $stream_file '        filename_interval="01-00-00_00:00:00"' . "\n";
-	print $stream_file '        clobber_mode="replace_files"' . "\n";
-	print $stream_file '        output_interval="none" >' . "\n";
+	print $stream_file '        clobber_mode="truncate"' . "\n";
+	print $stream_file '        output_interval="00-01-00_00:00:00">' . "\n";
+	print $stream_file '' . "\n";
+	print $stream_file '    <stream name="mesh"/>' . "\n";
+	print $stream_file '    <var name="xtime"/>' . "\n";
+	print $stream_file '    <var name="daysSinceStartOfSim"/>' . "\n";
+	print $stream_file '    <var name="iceAreaCell"/>' . "\n";
+	print $stream_file '    <var name="iceVolumeCell"/>' . "\n";
+	print $stream_file '    <var name="snowVolumeCell"/>' . "\n";
+	print $stream_file '    <var name="uVelocityGeo"/>' . "\n";
+	print $stream_file '    <var name="vVelocityGeo"/>' . "\n";
 	print $stream_file '' . "\n";
 	print $stream_file '</stream>' . "\n";
 	print $stream_file '' . "\n";

components/mpas-o/bld/build-namelist

@@ -671,6 +671,12 @@ add_default($nl, 'config_use_bulk_thickness_flux');
 add_default($nl, 'config_flux_attenuation_coefficient');
 add_default($nl, 'config_flux_attenuation_coefficient_runoff');
 
+############################
+# Namelist group: coupling #
+############################
+
+add_default($nl, 'config_ssh_grad_relax_timescale');
+
 ######################################
 # Namelist group: shortwaveRadiation #
 ######################################
@@ -1794,6 +1800,7 @@ my @groups = qw(run_modes
                 vmix_tanh
                 cvmix
                 forcing
+                coupling
                 shortwaveradiation
                 frazil_ice
                 land_ice_fluxes

components/mpas-o/bld/build-namelist-group-list

@@ -23,6 +23,7 @@ my @groups = qw(run_modes
                 vmix_tanh
                 cvmix
                 forcing
+                coupling
                 shortwaveradiation
                 frazil_ice
                 land_ice_fluxes

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

@@ -258,6 +258,12 @@ add_default($nl, 'config_use_bulk_thickness_flux');
 add_default($nl, 'config_flux_attenuation_coefficient');
 add_default($nl, 'config_flux_attenuation_coefficient_runoff');
 
+############################
+# Namelist group: coupling #
+############################
+
+add_default($nl, 'config_ssh_grad_relax_timescale');
+
 ######################################
 # Namelist group: shortwaveRadiation #
 ######################################

components/mpas-o/bld/namelist_files/namelist_defaults_mpas-o.xml

@@ -227,6 +227,9 @@
 <config_flux_attenuation_coefficient>0.001</config_flux_attenuation_coefficient>
 <config_flux_attenuation_coefficient_runoff>10.0</config_flux_attenuation_coefficient_runoff>
 
+<!-- coupling -->
+<config_ssh_grad_relax_timescale>86400.0</config_ssh_grad_relax_timescale>
+
 <!-- shortwaveRadiation -->
 <config_sw_absorption_type>'jerlov'</config_sw_absorption_type>
 <config_jerlov_water_type>3</config_jerlov_water_type>

components/mpas-o/bld/namelist_files/namelist_definition_mpas-o.xml

@@ -1157,6 +1157,17 @@ Default: Defined in namelist_defaults.xml
 </entry>
 
 
+<!-- coupling -->
+
+<entry id="config_ssh_grad_relax_timescale" type="real"
+	category="coupling" group="coupling">
+Timescale for relaxation of the ssh gradient for coupling. A value of 0.0 removes any relaxation and gives instantaneous response.
+
+Valid values: Any positive real number.
+Default: Defined in namelist_defaults.xml
+</entry>
+
+
 <!-- shortwaveRadiation -->
 
 <entry id="config_sw_absorption_type" type="char*1024"

components/mpas-o/cime_config/buildnml

@@ -64,8 +64,8 @@ if ( $OCN_GRID eq 'oEC60to30' ) {
 } elsif ( $OCN_GRID eq 'oEC60to30v3_ICG' ) {
         $grid_date .= '161222';
         $grid_prefix .= 'oEC60to30v3';
-        $ic_date .= '170213';
-        $ic_prefix .= 'oEC60to30v3.a16q_smRi_yr102';
+        $ic_date .= '170301';
+        $ic_prefix .= 'oEC60to30v3.restartFrom_anvil0221';
         $decomp_prefix .= 'mpas-o.graph.info.';
 } elsif ( $OCN_GRID eq 'oEC60to30wLI' ) {
         $grid_date .= '160830';
@@ -312,7 +312,8 @@ if ( -e "$CASEROOT/SourceMods/src.mpaso/$STREAM_NAME" ) {
 	print $stream_file '	<var name="indexSurfaceLayerDepth"/>' . "\n";
 	print $stream_file '	<var name="surfaceFrictionVelocity"/>' . "\n";
 	print $stream_file '	<var name="surfaceBuoyancyForcing"/>' . "\n";
-	print $stream_file '	<var name="seaSurfacePressure"/>' . "\n";
+	print $stream_file '	<var name="atmosphericPressure"/>' . "\n";
+	print $stream_file '	<var name="seaIcePressure"/>' . "\n";
 	print $stream_file '	<var name="frazilLayerThicknessTendency"/>' . "\n";
 	print $stream_file '	<var_struct name="tracersSurfaceFlux"/>' . "\n";
 	print $stream_file '	<var name="surfaceStressMagnitude"/>' . "\n";

components/mpas-o/driver/ocn_comp_mct.F

@@ -194,6 +194,11 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename )!{{{
     real (kind=RKIND) :: dt
     type (MPAS_timeInterval_type) :: timeStep
 
+    ! Added for coupling interval initialization
+    integer, pointer :: index_avgZonalSSHGradient, index_avgMeridionalSSHGradient
+    real (kind=RKIND), dimension(:),   pointer :: filteredSSHGradientZonal, filteredSSHGradientMeridional
+    real (kind=RKIND), dimension(:,:), pointer :: avgSSHGradient
+
     interface
        subroutine xml_stream_parser(xmlname, mgr_p, comm, ierr) bind(c)
           use iso_c_binding, only : c_char, c_ptr, c_int
@@ -628,6 +633,28 @@ end subroutine xml_stream_get_attributes
        end do
     end if
 
+!-----------------------------------------------------------------------
+!
+!   initialize coupling variables
+!   NOTE: could be moved to subroutine
+!
+!-----------------------------------------------------------------------
+
+    if ( runtype == "initial") then
+       block_ptr => domain % blocklist
+       do while(associated(block_ptr))
+          call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool)
+          call mpas_pool_get_dimension(forcingPool, 'index_avgSSHGradientZonal', index_avgZonalSSHGradient)
+          call mpas_pool_get_dimension(forcingPool, 'index_avgSSHGradientMeridional', index_avgMeridionalSSHGradient)
+          call mpas_pool_get_array(forcingPool, 'avgSSHGradient', avgSSHGradient)
+          call mpas_pool_get_array(forcingPool, 'filteredSSHGradientZonal', filteredSSHGradientZonal)
+          call mpas_pool_get_array(forcingPool, 'filteredSSHGradientMeridional', filteredSSHGradientMeridional)
+          filteredSSHGradientZonal      = avgSSHGradient(index_avgZonalSSHGradient,     :)
+          filteredSSHGradientMeridional = avgSSHGradient(index_avgMeridionalSSHGradient, :)
+          block_ptr => block_ptr % next
+       end do
+    endif
+
 !-----------------------------------------------------------------------
 !
 !   send intial state to driver
@@ -768,6 +795,13 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{
       logical, pointer :: config_write_output_on_startup
       character (len=StrKIND), pointer :: config_restart_timestamp_name
 
+      ! Added for coupling interval initialization
+      integer, pointer :: index_avgZonalSSHGradient, index_avgMeridionalSSHGradient
+      real (kind=RKIND), dimension(:),   pointer :: filteredSSHGradientZonal, filteredSSHGradientMeridional
+      real (kind=RKIND), dimension(:,:), pointer :: avgSSHGradient
+      real (kind=RKIND), pointer :: config_ssh_grad_relax_timescale
+      real (kind=RKIND) :: timeFilterFactor
+
       iam = domain % dminfo % my_proc_id
       call mpas_io_units_set_stdout(ocnStdOutUnit)
       call mpas_io_units_set_stderr(ocnStdErrUnit)
@@ -903,6 +937,31 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{
          write(stderrUnit,*) 'Completed timestep ', trim(timeStamp)
       end do
 
+      ! update coupled variables that get calculated on coupling intervals
+      ! NOTE: could be moved to subroutine
+      ! time filter ssh gradient
+      block_ptr => domain_ptr % blocklist
+      do while(associated(block_ptr))
+         call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool)
+         call mpas_pool_get_dimension(forcingPool, 'index_avgSSHGradientZonal', index_avgZonalSSHGradient)
+         call mpas_pool_get_dimension(forcingPool, 'index_avgSSHGradientMeridional', index_avgMeridionalSSHGradient)
+         call mpas_pool_get_array(forcingPool, 'avgSSHGradient', avgSSHGradient)
+         call mpas_pool_get_array(forcingPool, 'filteredSSHGradientZonal', filteredSSHGradientZonal)
+         call mpas_pool_get_array(forcingPool, 'filteredSSHGradientMeridional', filteredSSHGradientMeridional)
+         call mpas_pool_get_config(domain % configs, 'config_ssh_grad_relax_timescale', &
+                                                      config_ssh_grad_relax_timescale)
+         if (config_ssh_grad_relax_timescale < real(ocn_cpl_dt,RKIND)) then
+            timeFilterFactor = 1.0_RKIND
+         else
+            timeFilterFactor = real(ocn_cpl_dt,RKIND) / config_ssh_grad_relax_timescale
+         endif
+         filteredSSHGradientZonal = filteredSSHGradientZonal * (1.0_RKIND - timeFilterFactor) + &
+              avgSSHGradient(index_avgZonalSSHGradient, :) * timeFilterFactor
+         filteredSSHGradientMeridional = filteredSSHGradientMeridional * (1.0_RKIND - timeFilterFactor) + &
+              avgSSHGradient(index_avgMeridionalSSHGradient, :) * timeFilterFactor
+         block_ptr => block_ptr % next
+      end do
+
       ! Check if coupler wants us to write a restart file.
       ! We only write restart files at the end of a coupling interval
       if (seq_timemgr_RestartAlarmIsOn(EClock)) then
@@ -1391,8 +1450,8 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{
                                   seaIceSalinityFluxField, &
                                   riverRunoffFluxField, iceRunoffFluxField, &
                                   shortWaveHeatFluxField, rainFluxField, &
-                                  seaSurfacePressureField, iceFractionField, &
-                                  windSpeedSquared10mField, &
+                                  atmosphericPressureField, iceFractionField, &
+                                  seaIcePressureField, windSpeedSquared10mField, &
                                   iceFluxDICField, &
                                   iceFluxDONField, &
                                   iceFluxNO3Field, &
@@ -1417,8 +1476,8 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{
                                                seaIceSalinityFlux, &
                                                riverRunoffFlux, iceRunoffFlux, &
                                                shortWaveHeatFlux, rainFlux, &
-                                               seaSurfacePressure, iceFraction, &
-                                               windSpeedSquared10m, &
+                                               atmosphericPressure, iceFraction, &
+                                               seaIcePressure, windSpeedSquared10m, &
                                                iceFluxDIC,       &
                                                iceFluxDON, &
                                                iceFluxNO3, &
@@ -1485,7 +1544,8 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{
       call mpas_pool_get_field(forcingPool, 'iceRunoffFlux', iceRunoffFluxField)
       call mpas_pool_get_field(forcingPool, 'shortWaveHeatFlux', shortWaveHeatFluxField)
       call mpas_pool_get_field(forcingPool, 'rainFlux', rainFluxField)
-      call mpas_pool_get_field(forcingPool, 'seaSurfacePressure', seaSurfacePressureField)
+      call mpas_pool_get_field(forcingPool, 'atmosphericPressure', atmosphericPressureField)
+      call mpas_pool_get_field(forcingPool, 'seaIcePressure', seaIcePressureField)
       call mpas_pool_get_field(forcingPool, 'iceFraction', iceFractionField)
       call mpas_pool_get_field(forcingPool, 'iceRunoffFlux', iceRunoffFluxField)
 
@@ -1504,7 +1564,8 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{
       iceRunoffFlux => iceRunoffFluxField % array
       shortWaveHeatFlux => shortWaveHeatFluxField % array
       rainFlux => rainFluxField % array
-      seaSurfacePressure => seaSurfacePressureField % array
+      atmosphericPressure => atmosphericPressureField % array
+      seaIcePressure => seaIcePressureField % array
       iceFraction => iceFractionField % array
       iceRunoffFlux => iceRunoffFluxField % array
 
@@ -1607,11 +1668,12 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{
         if ( rainFluxField % isActive ) then
            rainFlux(i) = x2o_o % rAttr(index_x2o_Faxa_rain, n)
         end if
-        if ( seaSurfacePressureField % isActive ) then
-           ! Set seaSurfacePressure to be the sum of atmospheric bottom pressure 
-           ! and a limited version of sea ice pressure where it limits it to 5m of pressure
-           seaSurfacePressure(i) = x2o_o % rAttr(index_x2o_Sa_pbot,   n) &
-                                 + min( x2o_o % rAttr(index_x2o_Si_bpress, n), config_density0 * gravity * 5.0_RKIND )
+        if ( atmosphericPressureField % isActive ) then
+           atmosphericPressure(i) = x2o_o % rAttr(index_x2o_Sa_pbot,   n)
+        end if
+        if ( seaIcePressureField % isActive ) then
+           ! Set seaIcePressure to be limited to 5m of pressure
+           seaIcePressure(i) =  min( x2o_o % rAttr(index_x2o_Si_bpress, n), config_density0 * gravity * 5.0_RKIND )
         end if
         if ( iceFractionField % isActive ) then
            iceFraction(i) = x2o_o % rAttr(index_x2o_Si_ifrac, n)
@@ -1694,7 +1756,8 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{
    call mpas_pool_get_field(forcingPool, 'iceRunoffFlux', iceRunoffFluxField)
    call mpas_pool_get_field(forcingPool, 'shortWaveHeatFlux', shortWaveHeatFluxField)
    call mpas_pool_get_field(forcingPool, 'rainFlux', rainFluxField)
-   call mpas_pool_get_field(forcingPool, 'seaSurfacePressure', seaSurfacePressureField)
+   call mpas_pool_get_field(forcingPool, 'atmosphericPressure', atmosphericPressureField)
+   call mpas_pool_get_field(forcingPool, 'seaIcePressure', seaIcePressureField)
    call mpas_pool_get_field(forcingPool, 'iceFraction', iceFractionField)
    call mpas_pool_get_field(forcingPool, 'iceRunoffFlux', iceRunoffFluxField)
 
@@ -1772,8 +1835,11 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{
    if ( rainFluxField % isActive ) then
       call mpas_dmpar_exch_halo_field(rainFluxField)
    end if
-   if ( seaSurfacePressureField % isActive ) then
-      call mpas_dmpar_exch_halo_field(seaSurfacePressureField)
+   if ( atmosphericPressureField % isActive ) then
+      call mpas_dmpar_exch_halo_field(atmosphericPressureField)
+   end if
+   if ( seaIcePressureField % isActive ) then
+      call mpas_dmpar_exch_halo_field(seaIcePressureField)
    end if
    if ( iceFractionField % isActive ) then
       call mpas_dmpar_exch_halo_field(iceFractionField)
@@ -1868,8 +1934,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
 
    integer :: i, n
    integer, pointer :: nCellsSolve, index_temperatureSurfaceValue, index_salinitySurfaceValue, &
-                       index_avgZonalSurfaceVelocity, index_avgMeridionalSurfaceVelocity,      &
-                       index_avgZonalSSHGradient, index_avgMeridionalSSHGradient
+                       index_avgZonalSurfaceVelocity, index_avgMeridionalSurfaceVelocity
 
    type (block_type), pointer :: block_ptr
 
@@ -1884,6 +1949,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
    integer, dimension(:), pointer :: landIceMask
 
    real (kind=RKIND), dimension(:), pointer :: seaIceEnergy, accumulatedFrazilIceMass, frazilSurfacePressure, &
+                                               filteredSSHGradientZonal, filteredSSHGradientMeridional, &
                                                CO2Flux, DMSFlux, surfaceUpwardCO2Flux, &
                                                avgOceanSurfaceDIC, &
                                                avgOceanSurfaceDON, &
@@ -1896,7 +1962,7 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
                                                avgOceanSurfaceFeParticulate, &
                                                avgOceanSurfaceFeDissolved
 
-   real (kind=RKIND), dimension(:,:), pointer :: avgTracersSurfaceValue, avgSurfaceVelocity, avgSSHGradient, &
+   real (kind=RKIND), dimension(:,:), pointer :: avgTracersSurfaceValue, avgSurfaceVelocity, &
                                                  avgOceanSurfacePhytoC, &
                                                  avgOceanSurfaceDOC, layerThickness
 
@@ -1939,15 +2005,14 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
      call mpas_pool_get_dimension(forcingPool, 'index_avgSalinitySurfaceValue', index_salinitySurfaceValue)
      call mpas_pool_get_dimension(forcingPool, 'index_avgSurfaceVelocityZonal', index_avgZonalSurfaceVelocity)
      call mpas_pool_get_dimension(forcingPool, 'index_avgSurfaceVelocityMeridional', index_avgMeridionalSurfaceVelocity)
-     call mpas_pool_get_dimension(forcingPool, 'index_avgSSHGradientZonal', index_avgZonalSSHGradient)
-     call mpas_pool_get_dimension(forcingPool, 'index_avgSSHGradientMeridional', index_avgMeridionalSSHGradient)
 
      call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1)
 
      call mpas_pool_get_array(forcingPool, 'landIceMask', landIceMask)
      call mpas_pool_get_array(forcingPool, 'avgTracersSurfaceValue', avgTracersSurfaceValue)
      call mpas_pool_get_array(forcingPool, 'avgSurfaceVelocity', avgSurfaceVelocity)
-     call mpas_pool_get_array(forcingPool, 'avgSSHGradient', avgSSHGradient)
+     call mpas_pool_get_array(forcingPool, 'filteredSSHGradientZonal', filteredSSHGradientZonal)
+     call mpas_pool_get_array(forcingPool, 'filteredSSHGradientMeridional', filteredSSHGradientMeridional)
      if ( frazilIceActive ) then
         call mpas_pool_get_array(forcingPool, 'seaIceEnergy', seaIceEnergy)
         call mpas_pool_get_array(forcingPool, 'frazilSurfacePressure', frazilSurfacePressure)
@@ -1991,8 +2056,8 @@ subroutine ocn_export_mct(o2x_o, errorCode)   !{{{
        o2x_o % rAttr(index_o2x_So_u, n) = avgSurfaceVelocity(index_avgZonalSurfaceVelocity, i)
        o2x_o % rAttr(index_o2x_So_v, n) = avgSurfaceVelocity(index_avgMeridionalSurfaceVelocity, i)
 
-       o2x_o % rAttr(index_o2x_So_dhdx, n) = avgSSHGradient(index_avgZonalSSHGradient, i)
-       o2x_o % rAttr(index_o2x_So_dhdy, n) = avgSSHGradient(index_avgMeridionalSSHGradient, i)
+       o2x_o % rAttr(index_o2x_So_dhdx, n) = filteredSSHGradientZonal(i)
+       o2x_o % rAttr(index_o2x_So_dhdy, n) = filteredSSHGradientMeridional(i)
        if ( frazilIceActive ) then
           ! negative when frazil ice can be melted
           keepFrazil = .true.