Merge branch 'ocean/coastal' into ocean/develop

This PR is an accumulation of PRs into the ocean/coastal branch. We are bringing them in at once for efficiency: #285, #289, #284, #295, #310, #311, #312, #335, #354, #356, #358, #359, #365, #371
E3SM-Project · Oct 19, 2019 · c9ef38e · c9ef38e
2 parents cb866ed + abc3267
commit c9ef38e
Show file tree

Hide file tree

Showing 22 changed files with 1,891 additions and 314 deletions.
diff --git a/Registry.xml b/Registry.xml
@@ -95,6 +95,10 @@
 			description="Number of forcing groups"
 			definition="4"
 		/>
+		<dim name="maxTidalConstituents"
+			description="Max number of tidal potential constituents"
+			definition="37"
+		/>
 	</dims>
 
 	<!--*************************************-->
@@ -484,6 +488,10 @@
 					description="Inverse-time coefficient for the Rayleigh damping term, $c_R$, applied at every depth level."
 					possible_values="Any positive real value."
 		/>
+		<nml_option name="config_Rayleigh_damping_depth_variable" type="logical" default_value=".false." units="unitless"
+					description="If true applies r h^-1 instead of just r."
+					possible_values=".true. or .false."
+		/>
 		<nml_option name="config_Rayleigh_bottom_friction" type="logical" default_value=".false." units="unitless"
 					description="If true, Rayleigh friction is only applied to the bottom"
 					possible_values=".true. or .false."
@@ -845,6 +853,10 @@
 					possible_values="Real non-zero value."
 		/>
 		<nml_option name="config_tidal_forcing_type" type="character" default_value="off" units="unitless"
+					description="Selects the way tidal forcing is applied."
+					possible_values="'thickness_source','direct'"
+		/>
+		<nml_option name="config_tidal_forcing_model" type="character" default_value="off" units="unitless"
 					description="Selects the mode in which tidal forcing is computed."
 					possible_values="'off','monochromatic'"
 		/>
@@ -861,6 +873,56 @@
 					possible_values="Any positive real number."
 		/>
 	</nml_record>
+	<nml_record name="tidal_potential_forcing" mode="forward">
+		<nml_option name="config_use_tidal_potential_forcing" type="logical" default_value=".false." units="unitless"
+					description="Controls if tidal potential forcing is used."
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_tidal_potential_reference_time" type="character" default_value="0001-01-01_00:00:00" units="unitless"
+					description="Timestamp describing the time used to initialize nodal factors."
+					possible_values="'YYYY-MM-DD_HH:MM:SS'"
+		/>
+		<nml_option name="config_use_tidal_potential_forcing_M2" type="logical" default_value=".true." units="unitless"
+					description="Controls if tidal potential forcing for the M2 constituent is used."
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_use_tidal_potential_forcing_S2" type="logical" default_value=".true." units="unitless"
+					description="Controls if tidal potential forcing for the S2 constituent is used."
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_use_tidal_potential_forcing_N2" type="logical" default_value=".true." units="unitless"
+					description="Controls if tidal potential forcing for the N2 constituent is used."
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_use_tidal_potential_forcing_K2" type="logical" default_value=".true." units="unitless"
+					description="Controls if tidal potential forcing for the K2 constituent is used."
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_use_tidal_potential_forcing_K1" type="logical" default_value=".true." units="unitless"
+					description="Controls if tidal potential forcing for the K1 constituent is used."
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_use_tidal_potential_forcing_O1" type="logical" default_value=".true." units="unitless"
+					description="Controls if tidal potential forcing for the O1 constituent is used."
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_use_tidal_potential_forcing_Q1" type="logical" default_value=".true." units="unitless"
+					description="Controls if tidal potential forcing for the Q1 constituent is used."
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_use_tidal_potential_forcing_P1" type="logical" default_value=".true." units="unitless"
+					description="Controls if tidal potential forcing for the P1 constituent is used."
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_tidal_potential_ramp" type="real" default_value="10.0" units="days"
+					description="Number of days over which the tidal potential forcing is ramped"
+					possible_values="Any positive real number"
+		/>
+		<nml_option name="config_self_attraction_and_loading_beta" type="real" default_value="0.09" units="unitless"
+					description="Coefficient for SAL scalar approximation"
+					possible_values="0.0 to turn off, 0.09 is typical value to use for scalar approximation"
+		/>
+	</nml_record>
 	<nml_record name="frazil_ice" mode="forward">
 		<nml_option name="config_use_frazil_ice_formation" type="logical" default_value=".false." units="unitless"
 					description="Controls if fluxes related to frazil ice process are computed."
@@ -1031,6 +1093,10 @@
 					description="If 'upwind', use upwind to evaluate the edge-value for layerThickness, i.e., layerThicknessEdge.  The standard MPAS-O approach is 'centered'."
 					possible_values="'upwind', 'centered'"
 		/>
+		<nml_option name="config_drying_safety_height" type="real" default_value="1.0e-10" units="m"
+					description="Safety factor on minimum cell height to ensure the minimum height is not violated due to round-off."
+					possible_values="Real value greater or equal to 0."
+		/>
 	</nml_record>
 	<nml_record name="ocean_constants">
 		<nml_option name="config_density0" type="real" default_value="1026.0" units="kg m^{-3}"
@@ -1040,8 +1106,8 @@
 	</nml_record>
 	<nml_record name="pressure_gradient" mode="forward">
 		<nml_option name="config_pressure_gradient_type" type="character" default_value="pressure_and_zmid" units="unitless"
-					description="Form of pressure gradient terms in momentum equation. For most applications, the gradient of pressure and layer mid-depth are appropriate.  For isopycnal coordinates, one may use the gradient of the Montgomery potential."
-					possible_values="'pressure_and_zmid' or 'Jacobian_from_density' or 'Jacobian_from_TS' or 'MontgomeryPotential'"
+					description="Form of pressure gradient terms in momentum equation. For most applications, the gradient of pressure and layer mid-depth are appropriate.  For isopycnal coordinates, one may use the gradient of the Montgomery potential.  The sea surface height gradient (ssh_gradient) option is for barotropic, depth-averaged pressure."
+					possible_values="'ssh_gradient', 'pressure_and_zmid' or 'Jacobian_from_density' or 'Jacobian_from_TS' or 'MontgomeryPotential'"
 		/>
 		<nml_option name="config_common_level_weight" type="real" default_value="0.5" units="unitless"
 					description="The weight between standard Jacobian and weighted Jacobian, $\gamma$."
@@ -1377,6 +1443,7 @@
 		<package name="cullCells" description="This package controls variables that provide information on what cells should be culled."/>
 		<package name="tracerBudget" description="This package controls variables that record the tracer budget."/>
 		<package name="gm" description="This package is for GM variables, which are only needed at lower resolutions."/>
+		<package name="tidalPotentialForcingPKG" description="This package controls variables required for tidal potential forcing"/>
 	</packages>
 
 	<streams>
@@ -3152,6 +3219,10 @@
 			 description="Input mask for application of tidal forcing where 1 is applied tidal forcing"
 			 packages="tidalForcing"
 		/>
+		<var name="tidalBCValue" type="real" dimensions="nCells Time" units="m"
+			 description="Value of ssh height in cell for tidal boundary condition"
+			 packages="tidalForcing"
+		/>
 
 		<!-- Output fields for forcing due to tides -->
 		<var name="tidalLayerThicknessTendency" type="real" dimensions="nVertLevels nCells Time" units="m s^{-1}"
@@ -3178,6 +3249,60 @@
 			 description="surface pressure forcing due to weight of frazil ice"
 			 packages="frazilIce"
 		/>
+
+		<!-- Fields for tidal potential forcing -->
+		<var name="tidalPotentialEta" type="real" dimensions="nCells Time" units="m"
+			description="Equilibrium tidal potential"
+			packages="forwardMode"
+		/>
+		<var name="nTidalPotentialConstituents" type="integer" dimensions="Time" units="unitless"
+			description="Number of tidal constituents"
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="tidalPotentialConstituentFrequency" type="real" dimensions="maxTidalConstituents Time" units="s^{-1}"
+			description="Frequency of tidal constituents"
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="tidalPotentialConstituentAmplitude" type="real" dimensions="maxTidalConstituents Time" units="m"
+			description="Amplitude of tidal constituents"
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="tidalPotentialConstituentLoveNumbers" type="real" dimensions="maxTidalConstituents Time" units="unitless"
+			description="Love number combinations for tidal constituents"
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="tidalPotentialConstituentNodalAmplitude" type="real" dimensions="maxTidalConstituents Time" units="m"
+			description="Amplitude nodal factor for tidal constituents"
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="tidalPotentialConstituentNodalPhase" type="real" dimensions="maxTidalConstituents Time" units="s^{-1}"
+			description="Phase nodal factor for tidal constituents"
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="tidalPotentialConstituentAstronomical" type="real" dimensions="maxTidalConstituents Time" units="s^{-1}"
+			description="Astronomical argument for tidal constituents"
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="tidalPotentialConstituentType" type="integer" dimensions="maxTidalConstituents Time" units="unitless"
+			description="Spieces code for tidal constituents: long-period = 0, diurnal = 1, semi-diurnal = 2"
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="tidalPotentialLatitudeFunction" type="real" dimensions="nCells R3 Time" units="unitless"
+			description="Latitude function for tidal constituents: long-period = 3\sin^2(\phi)-1, diurnal = \sin(2\phi), semi-diurnal = \cos^2(\phi)"
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="tidalPotentialZMid" type="real" dimensions="nVertLevels nCells Time" units="m"
+			description="zMid - Equilibrium tidal potential in RK4"
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="sshSubcycleCurWithTides" type="real" dimensions="nCells Time" units="m"
+			description="SSH - tidal potential in split explicit "
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="sshSubcycleNewWithTides" type="real" dimensions="nCells Time" units="m"
+			description="SSH - tidal potential in split explicit "
+			packages="tidalPotentialForcingPKG"
+		/>
 	</var_struct>
 	<var_struct name="timeVaryingForcing" time_levs="1">
 		<var name="windSpeedU" type="real" dimensions="nCells Time" units="m/s"

diff --git a/driver/mpas_ocn_core_interface.F b/driver/mpas_ocn_core_interface.F
@@ -120,6 +120,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
       logical, pointer :: thicknessBulkPKGActive
       logical, pointer :: frazilIceActive
       logical, pointer :: tidalForcingActive
+      logical, pointer :: tidalPotentialForcingPKGActive
       logical, pointer :: inSituEOSActive
       logical, pointer :: variableShortwaveActive
       logical, pointer :: gmActive
@@ -148,6 +149,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
       logical, pointer :: config_use_freq_filtered_thickness
       logical, pointer :: config_use_frazil_ice_formation
       logical, pointer :: config_use_tidal_forcing
+      logical, pointer :: config_use_tidal_potential_forcing
       logical, pointer :: config_use_standardGM
       logical, pointer :: config_use_time_varying_atmospheric_forcing
       logical, pointer :: config_use_time_varying_land_ice_forcing
@@ -285,6 +287,15 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
          tidalForcingActive = .true.
       end if
 
+      !
+      ! test for use of tidal potential forcing, tidalPotentialForcingPKGActive
+      !
+      call mpas_pool_get_package(packagePool, 'tidalPotentialForcingPKGActive', tidalPotentialForcingPKGActive)
+      call mpas_pool_get_config(configPool, 'config_use_tidal_potential_forcing', config_use_tidal_potential_forcing)
+      if (config_use_tidal_potential_forcing) then
+         tidalPotentialForcingPKGActive = .true.
+      end if
+
       !
       ! test for form of pressure gradient computation
       !

diff --git a/mode_forward/mpas_ocn_forward_mode.F b/mode_forward/mpas_ocn_forward_mode.F
@@ -55,6 +55,7 @@ module ocn_forward_mode
    use ocn_surface_land_ice_fluxes
    use ocn_frazil_forcing
    use ocn_tidal_forcing
+   use ocn_tidal_potential_forcing
 
    use ocn_tracer_hmix
    use ocn_tracer_hmix_redi
@@ -221,6 +222,8 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{
       ierr = ior(ierr, err_tmp)
       call ocn_tidal_forcing_init(err_tmp)
       ierr = ior(ierr, err_tmp)
+      call ocn_tidal_potential_forcing_init(domain,err_tmp)
+      ierr = ior(ierr, err_tmp)
 
       call ocn_tracer_hmix_init(err_tmp)
       ierr = ior(ierr, err_tmp)

diff --git a/mode_forward/mpas_ocn_time_integration_rk4.F b/mode_forward/mpas_ocn_time_integration_rk4.F
@@ -144,6 +144,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
 
       ! Mesh array pointers
       integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop
+      real (kind=RKIND), dimension(:), pointer :: bottomDepth
 
       ! Provis Array Pointers
       real (kind=RKIND), dimension(:,:), pointer :: normalVelocityProvis, layerThicknessProvis
@@ -190,6 +191,13 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
       character (len=StrKIND) :: modifiedGroupName
       character (len=StrKIND) :: configName
 
+      ! Tidal boundary condition
+      logical, pointer :: config_use_tidal_forcing
+      character (len=StrKIND), pointer :: config_tidal_forcing_type
+      real (kind=RKIND), dimension(:), pointer :: tidalInputMask, tidalBCValue
+      real (kind=RKIND), dimension(:,:), pointer :: restingThickness
+      real (kind=RKIND) :: totalDepth
+
       ! Get config options
       call mpas_pool_get_config(domain % configs, 'config_mom_del4', config_mom_del4)
       call mpas_pool_get_config(domain % configs, 'config_filter_btr_mode', config_filter_btr_mode)
@@ -207,6 +215,8 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
       call mpas_pool_get_config(domain % configs, 'config_verify_not_dry', config_verify_not_dry)
       call mpas_pool_get_config(domain % configs, 'config_drying_min_cell_height', config_drying_min_cell_height)
       call mpas_pool_get_config(domain % configs, 'config_zero_drying_velocity', config_zero_drying_velocity)
+      call mpas_pool_get_config(domain % configs, 'config_use_tidal_forcing', config_use_tidal_forcing)
+      call mpas_pool_get_config(domain % configs, 'config_tidal_forcing_type', config_tidal_forcing_type)
 
       !
       ! Initialize time_levs(2) with state at current time
@@ -628,9 +638,11 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
          call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
          call mpas_pool_get_subpool(block % structs, 'diagnostics', diagnosticsPool)
          call mpas_pool_get_subpool(block % structs, 'scratch', scratchPool)
+         call mpas_pool_get_subpool(block % structs, 'verticalMesh', verticalMeshPool)
 
          call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
          call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges)
+         call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth)
 
          call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocityCur, 1)
          call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocityNew, 2)
@@ -657,6 +669,9 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
          call mpas_pool_get_array(diagnosticsPool, 'gradSSHMeridional', gradSSHMeridional)
          call mpas_pool_get_array(diagnosticsPool, 'surfaceVelocity', surfaceVelocity)
          call mpas_pool_get_array(diagnosticsPool, 'SSHGradient', SSHGradient)
+         call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness)
+         call mpas_pool_get_array(forcingPool, 'tidalInputMask', tidalInputMask)
+         call mpas_pool_get_array(forcingPool, 'tidalBCValue', tidalBCValue)
 
 
          if (config_prescribe_velocity) then
@@ -675,6 +690,26 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
             !$omp end do
          end if
 
+         ! direct application of tidal boundary condition
+         if (config_use_tidal_forcing .and. trim(config_tidal_forcing_type) == 'direct') then
+           do iCell=1, nCells
+             ! artificially assumes boolean mask for now, could generalize to tappered sponge layer
+             if (tidalInputMask(iCell) == 1.0_RKIND) then
+               ! compute total depth for relative thickness contribution
+               totalDepth = 0.0_RKIND
+               do k = 1, maxLevelCell(iCell)
+                 totalDepth = totalDepth + restingThickness(k,iCell)
+               end do
+
+               ! only modify layer thicknesses on tidal boundary
+               do k = 1, maxLevelCell(iCell)
+                 layerThicknessNew(k, iCell) = tidalInputMask(iCell)*(tidalBCValue(iCell) + bottomDepth(iCell))*(restingThickness(k,iCell)/totalDepth)
+                 !(1.0_RKIND - tidalInputMask(iCell))*layerThicknessNew(k, iCell)  ! generalized tappered assumption code
+               end do
+             end if
+           end do
+         end if
+
          call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, diagnosticsPool, scratchPool, tracersPool, 2)
          call mpas_threading_barrier()