Modifications needed for time filtered ssh gradient coupling

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.