Skip to content

Commit

Permalink
Merge branch 'sbrus89/ocn/wave_drag' (PR #6310)
Browse files Browse the repository at this point in the history 
Add SAL updates and topographic wave drag schemes for barotropic tides

This PR incorporates updates to the global tidal forcing in MPAS-Ocean.
It combines two contributions developed under the ICoM project:
* @dengwirda's modifications to
  a) include surface pressure in the self attraction and loading (SAL)
     calculations
  b) use a depth rather than distance-based criteria in the SAL coastal
     smoothing factor, and
  c) update the harmonic analysis AM to account for initial SSH offsets.
* @knbarton's refactoring of the existing Jayne and St. Laurent
  topographic wave drag (TWD) parameterization and the addition of the
  Zaron and Egbert and Local Generation Formula TWD schemes.

[NML]
[BFB] - mpas-ocean standalone
  • Loading branch information
@jonbob
jonbob committed May 20, 2024
2 parents 41a0ccf + 6279764 commit b25c065
Show file tree
Hide file tree
Showing 17 changed files with 464 additions and 157 deletions.
5 changes: 4 additions & 1 deletion components/mpas-ocean/bld/build-namelist
View file
Original file line number Diff line number Diff line change
Expand Up @@ -720,7 +720,7 @@ add_default($nl, 'config_enable_shortwave_energy_fixer');
###########################################

add_default($nl, 'config_use_self_attraction_loading');
add_default($nl, 'config_self_attraction_loading_smoothing_width');
add_default($nl, 'config_self_attraction_loading_depth_cutoff');
add_default($nl, 'config_mpas_to_grid_weights_file');
add_default($nl, 'config_grid_to_mpas_weights_file');
add_default($nl, 'config_self_attraction_loading_compute_interval');
Expand Down Expand Up @@ -820,8 +820,11 @@ add_default($nl, 'config_loglaw_bottom_drag_min');
add_default($nl, 'config_loglaw_bottom_drag_max');
add_default($nl, 'config_explicit_bottom_drag_coeff');
add_default($nl, 'config_use_topographic_wave_drag');
add_default($nl, 'config_topographic_wave_drag_scheme');
add_default($nl, 'config_topographic_wave_drag_coeff');
add_default($nl, 'config_thickness_drag_type');
add_default($nl, 'config_topographic_wave_drag_cutoff_depth');
add_default($nl, 'config_topographic_wave_drag_cutoff_width');

####################################
# Namelist group: Rayleigh_damping #
Expand Down
5 changes: 4 additions & 1 deletion components/mpas-ocean/bld/build-namelist-section
View file
Original file line number Diff line number Diff line change
Expand Up @@ -245,7 +245,7 @@ add_default($nl, 'config_enable_shortwave_energy_fixer');
###########################################

add_default($nl, 'config_use_self_attraction_loading');
add_default($nl, 'config_self_attraction_loading_smoothing_width');
add_default($nl, 'config_self_attraction_loading_depth_cutoff');
add_default($nl, 'config_mpas_to_grid_weights_file');
add_default($nl, 'config_grid_to_mpas_weights_file');
add_default($nl, 'config_self_attraction_loading_compute_interval');
Expand Down Expand Up @@ -337,8 +337,11 @@ add_default($nl, 'config_loglaw_bottom_drag_min');
add_default($nl, 'config_loglaw_bottom_drag_max');
add_default($nl, 'config_explicit_bottom_drag_coeff');
add_default($nl, 'config_use_topographic_wave_drag');
add_default($nl, 'config_topographic_wave_drag_scheme');
add_default($nl, 'config_topographic_wave_drag_coeff');
add_default($nl, 'config_thickness_drag_type');
add_default($nl, 'config_topographic_wave_drag_cutoff_depth');
add_default($nl, 'config_topographic_wave_drag_cutoff_width');

####################################
# Namelist group: Rayleigh_damping #
Expand Down
5 changes: 4 additions & 1 deletion components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -347,7 +347,7 @@

<!-- self_attraction_loading -->
<config_use_self_attraction_loading>.false.</config_use_self_attraction_loading>
<config_self_attraction_loading_smoothing_width>1.0</config_self_attraction_loading_smoothing_width>
<config_self_attraction_loading_depth_cutoff>10.0</config_self_attraction_loading_depth_cutoff>
<config_mpas_to_grid_weights_file>'mpas_to_grid.nc'</config_mpas_to_grid_weights_file>
<config_grid_to_mpas_weights_file>'grid_to_mpas.nc'</config_grid_to_mpas_weights_file>
<config_self_attraction_loading_compute_interval>'0000-00-00_00:30:00'</config_self_attraction_loading_compute_interval>
Expand Down Expand Up @@ -476,8 +476,11 @@
<config_loglaw_bottom_drag_max>1.0e-1</config_loglaw_bottom_drag_max>
<config_explicit_bottom_drag_coeff>1.0e-3</config_explicit_bottom_drag_coeff>
<config_use_topographic_wave_drag>.false.</config_use_topographic_wave_drag>
<config_topographic_wave_drag_scheme>'ZAE'</config_topographic_wave_drag_scheme>
<config_topographic_wave_drag_coeff>5.0e-4</config_topographic_wave_drag_coeff>
<config_thickness_drag_type>'centered'</config_thickness_drag_type>
<config_topographic_wave_drag_cutoff_depth>500</config_topographic_wave_drag_cutoff_depth>
<config_topographic_wave_drag_cutoff_width>10</config_topographic_wave_drag_cutoff_width>

<!-- Rayleigh_damping -->
<config_Rayleigh_damping_coeff>0.0</config_Rayleigh_damping_coeff>
Expand Down
30 changes: 27 additions & 3 deletions components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1245,9 +1245,9 @@ Valid values: .true. or .false.
Default: Defined in namelist_defaults.xml
</entry>

<entry id="config_self_attraction_loading_smoothing_width" type="real"
<entry id="config_self_attraction_loading_depth_cutoff" type="real"
category="self_attraction_loading" group="self_attraction_loading">
Defines region over which ssh is smoothed to zero at coasts for SAL calculation.
Defines depths over which ssh is smoothed to zero for SAL calculation.

Valid values: Any positive real number.
Default: Defined in namelist_defaults.xml
Expand Down Expand Up @@ -1796,11 +1796,19 @@ Valid values: .true. or .false.
Default: Defined in namelist_defaults.xml
</entry>

<entry id="config_topographic_wave_drag_scheme" type="char*1024"
category="bottom_drag" group="bottom_drag">
Which of the following three wave drag schemes to use: JSL: Jayne and St. Laurent; ZAE: Zaron and Egbert; or LGF: Local generation formula

Valid values: JSL or ZAE or LGF
Default: Defined in namelist_defaults.xml
</entry>

<entry id="config_topographic_wave_drag_coeff" type="real"
category="bottom_drag" group="bottom_drag">
Dimensionless topographic wave drag coefficient, $c_{topo}$.

Valid values: any positive real, typically 5.0e-4
Valid values: any positive real
Default: Defined in namelist_defaults.xml
</entry>

Expand All @@ -1812,6 +1820,22 @@ Valid values: 'harmonic-mean', 'centered'
Default: Defined in namelist_defaults.xml
</entry>

<entry id="config_topographic_wave_drag_cutoff_depth" type="real"
category="bottom_drag" group="bottom_drag">
Specifies minimum depth at which topographic wave drag is applied

Valid values: any positive real
Default: Defined in namelist_defaults.xml
</entry>

<entry id="config_topographic_wave_drag_cutoff_width" type="real"
category="bottom_drag" group="bottom_drag">
Width over which the topographic wave drag is progressively turned off

Valid values: any positive real
Default: Defined in namelist_defaults.xml
</entry>


<!-- Rayleigh_damping -->

Expand Down
80 changes: 69 additions & 11 deletions components/mpas-ocean/src/Registry.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -840,8 +840,8 @@
description="Controls if self-attraction and loading is applied to ssh"
possible_values=".true. or .false."
/>
<nml_option name="config_self_attraction_loading_smoothing_width" type="real" default_value="1.0" units="km"
description="Defines region over which ssh is smoothed to zero at coasts for SAL calculation."
<nml_option name="config_self_attraction_loading_depth_cutoff" type="real" default_value="10.0" units="m"
description="Defines depths over which ssh is smoothed to zero for SAL calculation."
possible_values="Any positive real number."
/>
<nml_option name="config_mpas_to_grid_weights_file" type="character" default_value="mpas_to_grid.nc"
Expand Down Expand Up @@ -1119,14 +1119,26 @@
description="If true, topographic wave drag is used on the momentum equation."
possible_values=".true. or .false."
/>
<nml_option name="config_topographic_wave_drag_coeff" type="real" default_value="5.0e-4"
<nml_option name="config_topographic_wave_drag_scheme" type="character" default_value="ZAE" units="unitless"
description="Which of the following three wave drag schemes to use: JSL: Jayne and St. Laurent; ZAE: Zaron and Egbert; or LGF: Local generation formula"
possible_values="JSL or ZAE or LGF"
/>
<nml_option name="config_topographic_wave_drag_coeff" type="real" default_value="0.3" units="unitless"
description="Dimensionless topographic wave drag coefficient, $c_{topo}$."
possible_values="any positive real, typically 5.0e-4"
possible_values="any positive real"
/>
<nml_option name="config_thickness_drag_type" type="character" default_value="centered"
description="The type of layerThickness averaging to use on the drag term. The standard MPAS-O approach is 'centered'."
possible_values="'harmonic-mean', 'centered'"
/>
<nml_option name="config_topographic_wave_drag_cutoff_depth" type="real" default_value="500" units="m"
description="Specifies minimum depth at which topographic wave drag is applied"
possible_values="any positive real"
/>
<nml_option name="config_topographic_wave_drag_cutoff_width" type="real" default_value="10" units="m"
description="Width over which the topographic wave drag is progressively turned off"
possible_values="any positive real"
/>
</nml_record>
<nml_record name="Rayleigh_damping" mode="forward">
<nml_option name="config_Rayleigh_damping_coeff" type="real" default_value="1.0e-4" units="s^-1"
Expand Down Expand Up @@ -1988,7 +2000,7 @@
<var name="gotmEpsbTopOfCell"/>
<var name="gotmDissTopOfCell"/>
<var name="gotmLengthTopOfCell"/>
<var name="ssh_sal"/>
<var name="pgf_sal"/>
</stream>

<stream name="output"
Expand Down Expand Up @@ -2288,7 +2300,13 @@
filename_template="topographic_wave_drag.nc"
runtime_format="single_file"
mode="forward;analysis">
<var name="topographic_wave_drag"/>
<var name="topo_rinv"/>
<var name="topo_buoyancy_N1V"/>
<var name="topo_buoyancy_N1B"/>
<var name="bathy_stddev"/>
<var name="bed_slope_edges"/>
<var name="lonGradEdge"/>
<var name="latGradEdge"/>
</stream>

</streams>
Expand Down Expand Up @@ -3324,10 +3342,34 @@
<var name="surfaceFluxAttenuationCoefficientRunoff" type="real" dimensions="nCells Time" units="m"
description="The spatially-dependent length scale of exponential decay of river runoff. Fluxes are multiplied by $e^{z/\gamma}$, where this coefficient is $\gamma$."
/>
<var name="topographic_wave_drag" type="real" dimensions="nEdges" units="1"
<var name="topo_rinv" type="real" dimensions="nEdges" units="unitless"
description="wave drag coefficient or 1/(rinv) where rinv is the e-folding time used in HyCOM"
packages="topographicWaveDragPKG"
/>
<var name="topo_buoyancy_N1V" type="real" dimensions="nEdges" units="unitless"
description="Vertically averaged buoyancy frequencies from WOA 2013 data"
packages="topographicWaveDragPKG"
/>
<var name="topo_buoyancy_N1B" type="real" dimensions="nEdges" units="unitless"
description="Bottom buoyancy frequencies from WOA 2013 data"
packages="topographicWaveDragPKG"
/>
<var name="bathy_stddev" type="real" dimensions="nEdges" units="unitless"
description="Standard deviation of the bathymetry within a given cell"
packages="topographicWaveDragPKG"
/>
<var name="bed_slope_edges" type="real" dimensions="nEdges" units="unitless"
description="The integral gradient average per cell"
packages="topographicWaveDragPKG"
/>
<var name="lonGradEdge" type="real" dimensions="nEdges" units="unitless"
description="Meridional gradient of bathymetry"
packages="topographicWaveDragPKG"
/>
<var name="latGradEdge" type="real" dimensions="nEdges" units="unitless"
description="Zonal gradient of bathymetry"
packages="topographicWaveDragPKG"
/>
<!-- diagnostic fields for land-ice fluxes -->
<var name="landIceFrictionVelocity" type="real" dimensions="nCells Time" units="m s^-1"
description="The friction velocity $u_*$ under land ice"
Expand Down Expand Up @@ -3510,12 +3552,12 @@
description="A vector used in the split-implicit barotropic mode solver"
packages="semiImplicitTimePKG"
/>
<var name="ssh_sal" type="real" dimensions="nCells Time" units="m"
description="sea surface height perturbation from self-attraction and loading"
<var name="pgf_sal" type="real" dimensions="nCells Time" units="m"
description="Barotropic pressure perturbation due to self-attraction and loading"
packages="tidalPotentialForcingPKG"
/>
<var name="ssh_sal_grad" type="real" dimensions="nEdges Time" units=""
description="gradient of sea surface height perturbation from self-attraction and loading"
<var name="pgf_sal_grad" type="real" dimensions="nEdges Time" units=""
description="Gradient of barotropic pressure perturbation due to self-attraction and loading"
packages="tidalPotentialForcingPKG"
/>
<!-- FIELD split_explicit AB2 time stepping -->
Expand All @@ -3535,6 +3577,22 @@
description="intermediate estimate ssh for wetting/drying purposes"
packages="forwardMode"
/>
<var name="temp_twd" type="real" dimensions="nEdges Time" units=""
description="tendency perturbation due to topographic wave drag"
packages="tidalPotentialForcingPKG"
/>
<var name="normalCoeffTWD" type="real" dimensions="nEdges" units=""
description="Coefficient multiplied by the normal velocity in the topographic wave drag"
packages="tidalPotentialForcingPKG"
/>
<var name="tangentialCoeffTWD" type="real" dimensions="nEdges" units=""
description="Coefficient multiplied by the tangential velocity in the topographic wave drag"
packages="tidalPotentialForcingPKG"
/>
<var name="topographicWaveDrag" type="real" dimensions="nEdges" units=""
description="topographic wave drag multiplicative factor"
packages="tidalPotentialForcingPKG"
/>
</var_struct>
<var_struct name="shortwave" time_levs="1">
<!-- **********************************************************************
Expand Down
4 changes: 4 additions & 0 deletions components/mpas-ocean/src/analysis_members/Registry_harmonic_analysis.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,9 @@
<var name="analysisConstituentNodalPhase" type="real" dimensions="maxTidalConstituents Time" units="rad/s"
description="Frequency of each constituent"
/>
<var name="sshInit" type="real" dimensions="nCells" units="m"
description="Copy of the initial ssh distribution. Used to compute analysis with respect to ssh anomalies, ensuring cells with initial ssh offsets are handled correctly"
/>
<var name="leastSquaresLHSMatrix" type="real" dimensions="maxTidalConstituentsX2 maxTidalConstituentsX2 Time" units="rad/s"
description="Frequency of each constituent"
/>
Expand Down Expand Up @@ -154,6 +157,7 @@
output_interval="stream:restart:output_interval"
immutable="false"
mode="forward;analysis">
<var name="sshInit"/>
<var name="leastSquaresLHSMatrix"/>
<var name="leastSquaresRHSVector"/>
</stream>
Expand Down
36 changes: 27 additions & 9 deletions components/mpas-ocean/src/analysis_members/mpas_ocn_harmonic_analysis.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,11 +117,12 @@ subroutine ocn_init_harmonic_analysis(domain, err)!{{{

type (dm_info) :: dminfo
type (block_type), pointer :: block
type (mpas_pool_type), pointer :: harmonicAnalysisAMPool
type (mpas_pool_type), pointer :: statePool, harmonicAnalysisAMPool

integer, pointer :: nAnalysisConstituents
integer, pointer :: nAnalysisConstituents, nCellsSolve
real (kind=RKIND), pointer :: harmonicAnalysisStart
real (kind=RKIND), pointer :: harmonicAnalysisEnd
real (kind=RKIND), dimension(:), pointer :: ssh, sshInit
real (kind=RKIND), dimension(:,:), pointer :: leastSquaresLHSMatrix
real (kind=RKIND), dimension(:,:), pointer :: leastSquaresRHSVector
real (kind=RKIND), dimension(:,:), pointer :: decomposedConstituentAmplitude
Expand All @@ -135,7 +136,7 @@ subroutine ocn_init_harmonic_analysis(domain, err)!{{{
integer, dimension(:), allocatable :: tidalConstituentType

type (MPAS_Time_Type) :: refTime
integer :: iCon
integer :: iCon, iCell

err = 0

Expand All @@ -145,8 +146,13 @@ subroutine ocn_init_harmonic_analysis(domain, err)!{{{

block => domain % blocklist
do while (associated(block))
call mpas_pool_get_subpool(block % structs, 'state', statePool)
call mpas_pool_get_subpool(block % structs, 'harmonicAnalysisAM', harmonicAnalysisAMPool)

call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', nCellsSolve)
call mpas_pool_get_array(statePool, 'ssh', ssh, 1)

call mpas_pool_get_array(harmonicAnalysisAMPool, 'sshInit', sshInit)
call mpas_pool_get_array(harmonicAnalysisAMPool, 'nAnalysisConstituents', nAnalysisConstituents)
call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresLHSMatrix', leastSquaresLHSMatrix)
call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresRHSVector', leastSquaresRHSVector)
Expand Down Expand Up @@ -214,11 +220,14 @@ subroutine ocn_init_harmonic_analysis(domain, err)!{{{
tidalConstituentType, &
err)

if (.not. config_do_restart) then
leastSquaresRHSVector = 0.0_RKIND
leastSquaresLHSMatrix = 0.0_RKIND
end if
if (.not. config_do_restart) then
leastSquaresRHSVector = 0.0_RKIND
leastSquaresLHSMatrix = 0.0_RKIND
end if

do iCell = 1,nCellSSolve
sshInit(iCell) = ssh(iCell)
end do

do iCon = 1,nAnalysisConstituents
call mpas_log_write('Constituent '//constituentList(iCon)%constituent)
Expand Down Expand Up @@ -291,7 +300,8 @@ subroutine ocn_compute_harmonic_analysis(domain, timeLevel, err)!{{{

integer, pointer :: nCellsSolve
integer, pointer :: nAnalysisConstituents
real (kind=RKIND), dimension(:), pointer :: ssh
real (kind=RKIND), dimension(:), pointer :: ssh, sshInit
real (kind=RKIND), dimension(:), allocatable :: sshDiff
real (kind=RKIND), dimension(:,:), pointer :: leastSquaresLHSMatrix
real (kind=RKIND), dimension(:,:), pointer :: leastSquaresRHSVector
real (kind=RKIND), dimension(:), pointer :: analysisConstituentFrequency
Expand Down Expand Up @@ -334,6 +344,7 @@ subroutine ocn_compute_harmonic_analysis(domain, timeLevel, err)!{{{

call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', nCellsSolve)
call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevel)
call mpas_pool_get_array(harmonicAnalysisAMPool, 'sshInit', sshInit)
call mpas_pool_get_array(harmonicAnalysisAMPool, 'nAnalysisConstituents', nAnalysisConstituents)
call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresLHSMatrix', leastSquaresLHSMatrix)
call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresRHSVector', leastSquaresRHSVector)
Expand All @@ -355,6 +366,12 @@ subroutine ocn_compute_harmonic_analysis(domain, timeLevel, err)!{{{

! update if within harmonic analysis period
if (daysSinceStartOfSim .le. harmonicAnalysisEnd) then
allocate(sshDiff(nCellsSolve))

do iCell = 1,nCellsSolve
sshDiff(iCell) = ssh(iCell) - sshInit(iCell)
end do

CALL update_least_squares_LHS_matrix(nAnalysisConstituents, &
time, &
analysisConstituentFrequency, &
Expand All @@ -364,8 +381,9 @@ subroutine ocn_compute_harmonic_analysis(domain, timeLevel, err)!{{{
time, &
nCellsSolve, &
analysisConstituentFrequency, &
ssh, &
sshDiff, &
leastSquaresRHSVector)
deallocate(sshDiff)
call mpas_log_write('harmonicAnalysisAM update')
end if

Expand Down
Loading

0 comments on commit b25c065

Please sign in to comment.