Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add SAL updates and topographic wave drag schemes for barotropic tides #6310

Merged
merged 19 commits into from
May 20, 2024
Merged

Add SAL updates and topographic wave drag schemes for barotropic tides #6310

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
67d16df
Update SAL and AMHarmonicAnalysis to respect initial SSH offsets
dengwirda Apr 25, 2023
14e19d3
Update SAL and AMHarmonicAnalysis to respect initial SSH offsets
dengwirda Apr 25, 2023
992db9b
Update SAL and AMHarmonicAnalysis to respect initial SSH offsets
dengwirda Apr 25, 2023
3c32bd4
Update SAL and AMHarmonicAnalysis to respect initial SSH offsets
dengwirda Apr 25, 2023
75ebe3d
Added two new topographic wave drag schemes.
knbarton Nov 18, 2022
d02d740
Added local generation formula topograhpic wave drag scheme
knbarton Feb 13, 2023
687c935
Remove double multiplication by topographic wave drag coefficient
knbarton Feb 21, 2023
f5eaa4b
Fix issue where topographic wave drag was being multiplied by tendanc…
knbarton Feb 27, 2023
987126d
Fix initialization of bottom depth in init mode
knbarton Feb 27, 2023
d0d376d
Adds 1/H term into the JSL and ZAE schemes in the topographic wave drag
knbarton Feb 27, 2023
b3dda60
Changed 1/H variable from layerThickness to bottomDepthEdge
knbarton Feb 28, 2023
395f060
Fixed sign of bottom depth in topographic wave drag routine
knbarton Mar 3, 2023
bf001b1
Minor changes to ZAE topographic wave drag scheme
knbarton Mar 20, 2023
9c588a8
updates to topographic wave drag routine to prevent NaNs
knbarton May 2, 2023
0857b8c
Removed references to unused variable bottomDepthEdgeTWD
knbarton May 2, 2023
617cef2
Fix rebase issues with ab2
sbrus89 Mar 27, 2024
d67e580
Topographic wave drag cleanup
sbrus89 Apr 1, 2024
ea270ea
Make ZAE the default wave drag scheme
sbrus89 Apr 4, 2024
6279764
Update bld files to match mpaso Registry changes
jonbob May 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion components/mpas-ocean/bld/build-namelist
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -718,7 +718,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 @@ -818,8 +818,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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -243,7 +243,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 @@ -335,8 +335,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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -295,7 +295,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 @@ -399,8 +399,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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1229,9 +1229,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 @@ -1780,11 +1780,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 @@ -1796,6 +1804,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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -829,8 +829,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 @@ -1108,14 +1108,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"
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the implementation/tuning for LGF well established? I know that JSL and especially ZAE are working well, but am less confident re: LGF.
Should ZAE be the default if it's typically the best performing?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Having ZAE as the default sounds good to me. I'm going to do a full tidal run with all three and I'll report back on the errors. I believe that LGF needs to be re-evaluated given the indexing issue I found above.

/>
<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 @@ -1962,7 +1974,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 @@ -2262,7 +2274,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 @@ -3233,10 +3251,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 @@ -3419,12 +3461,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 @@ -3436,6 +3478,22 @@
description="Coriolis term at the previous time step used in the split-explicit AB2 time stepping"
packages="splitAB2TimeIntegrator"
/>
<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
Open in desktop
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
Open in desktop
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

Comment on lines +369 to +374
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the change that makes the harmonic analysis work on the ssh anomaly.

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