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

sbrus89 · 2024-03-27T22:25:20Z

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. See also: Self-attraction-and-loading and harmonic analysis for SSH != 0 E3SM-Ocean-Discussion/E3SM#52.
@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.

This PR has been tested using a new compass case, which performs a 125-day tidal run (with a 90 day analysis period) on the VR45to5 mesh using the Zaron and Egbert topographic wave drag option. This results in a 3.3 cm M2 deep water RMSE as compared to the TPXO9 tidal database, which is the most accurate MPAS-Ocean tidal result to date. Here is a link to the compass PR with testing results: MPAS-Dev/compass#802

[NML]
[BFB] - mpas-ocean standalone

- Include pressure forcing term in SAL computation to account for under isc and atm pressure cases. - Update SAL coastal smoothing factor to mask shallow cells where wet-dry effects may occur out of SAL computation. - Compute harmonic analysis AM wrt. SSH-SSH_init to account for initially non-zero SSH distributions.

- Add missing rho0gInv definition.

- Clean up descriptions, comments and use real-valued switches.

- Change sshSmoothed to btrPressure in SAL routines, since we're now working with surface pressures too.

…y instead of subtracted.

sbrus89 · 2024-03-28T21:27:32Z

components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F

-        ! Set up coastal ssh smoothing
-        allocate(sshSmoothed(nCellsAll))
-
-        if (config_self_attraction_loading_smoothing_width > 2.0_RKIND) then
-           trans_width = config_self_attraction_loading_smoothing_width*1000.0_RKIND
-           trans_start = 0.0_RKIND
-           do iCell = 1,nCellsAll
-               d = distanceToCoast(iCell)
-               coastalSmoothingFactor(iCell) = 0.5_RKIND*(tanh((d - trans_start - 0.5_RKIND * trans_width) / (0.2_RKIND *trans_width)) + 1.0_RKIND)
-           enddo
-        else
-           coastalSmoothingFactor(:) = 1.0_RKIND
-        endif
+        ! Set up coastal smoothing
+        allocate(btrPressure(nCellsAll))

+        ! Set-up smoothing mask for shallow regions where SAL should not
+        ! be computed across e.g. wet-dry cells
+        depth_tol = config_self_attraction_loading_depth_cutoff
+        do iCell = 1,nCellsAll
+           depth = ssh(iCell) + bottomDepth(iCell)
+           coastalSmoothingFactor(iCell) = 0.5_RKIND*(tanh( &
+              (depth - 0.5_RKIND*depth_tol) / (0.2_RKIND*depth_tol)) + 1.0_RKIND)
+        enddo
+


This is the switch from distance-based to depth-based SAL smoothing

sbrus89 · 2024-03-28T21:28:20Z

components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F

-          sshSmoothed(iCell) = coastalSmoothingFactor(iCell)*ssh(iCell)
+          ! compute SAL wrt. the full barotropic bottom pressure
+          btrPressure(iCell) = coastalSmoothingFactor(iCell) * ( &
+              ssh(iCell) + rho0gInv * surfacePressure(iCell) )


This is where the surface pressure is added to the SAL calculation

sbrus89 · 2024-03-28T21:29:57Z

components/mpas-ocean/src/analysis_members/mpas_ocn_harmonic_analysis.F

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


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

sbrus89 · 2024-03-28T22:14:10Z

components/mpas-ocean/src/Registry.xml

+            <var name="temp_twd"/>
+            <var name="normalCoeff"/>
+            <var name="tangentialCoeff"/>
+            <var name="topographicWaveDrag"/>


I don't believe these need to be in the restart stream, since they are calculated in an init routine

sbrus89 · 2024-03-28T22:30:56Z

components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F

+                      topographicWaveDrag = topographicWaveDragCoeff  &
+                                            * topoDragTurnOff         &
+                                            * (sqrt((topo_buoyancy_N1B(iEdge)**2 - omegaM2**2) * &
+                                                    (topo_buoyancy_N1V(iEdge)**2 - omegaM2**2))) &
+                                            / (4*pii*omegaM2)
+                      normalCoeff = (lonGradEdge(iEdge)**2) * &
+                                    (cos(angleEdge(iEdge))**2) + (latGradEdge(iEdge)**2)*(sin(angleEdge(iEdge))**2) - &
+                                        2*latGradEdge(iEdge)*lonGradEdge(iEdge)*sin(angleEdge(iEdge))*(cos(angleEdge(iEdge)))
+                      tangentialCoeff = sin(angleEdge(iEdge)) * &
+                                        cos(angleEdge(iEdge)) * &
+                                        (lonGradEdge(iEdge)**2 - latGradEdge(iEdge)**2) + &
+                                        latGradEdge(iEdge)*lonGradEdge(iEdge)*(cos(angleEdge(iEdge))**2 - sin(angleEdge(iEdge))**2)


These variables should be indexed, i.e., topographicWaveDrag(iEdge) = , normalCoeff(iEdge) =, tangentialCoeff(iEdge) =

I'm not sure what fortran will do in this case, but I feel like testing for LGF should be revisited...

I'm pretty sure it'll set the entire array to a single value, so this definitely needs to be retested.

sbrus89 · 2024-03-28T22:31:45Z

components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F

+                  topographicWaveDrag(iEdge) = 0.0_RKIND
+                  if (kmax > 0) then
+                      k =  maxLevelEdgeTop(iEdge)
+                      topoDragTurnOff = (TANH((layerThickEdgeMean(k,iEdge)-topoDragCutoffDepth)/topoDragCutoffWidth)+1)/2


I think I'm going to pull out the application of the topoDragTurnOff factor since it's common to all three drag schemes

sbrus89 · 2024-03-28T22:34:02Z

components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F

+            ! topographic_wave_drag = 1/rinv
+            if (k>0) then
+              temp_twd(iEdge) = topographicWaveDrag(iEdge) * normalVelocity(k,iEdge)
+              tend(k,iEdge) = tend(k,iEdge) - topographicWaveDrag(iEdge) * normalVelocity(k,iEdge)


We should just set tend(k,iEdge) = tend(k,iEdge) - temp_twd(iEdge)

sbrus89 · 2024-03-28T22:39:20Z

components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F

+                      topoDragTurnOff = (TANH((layerThickEdgeMean(k,iEdge)-topoDragCutoffDepth)/topoDragCutoffWidth)+1)/2
+                      topographicWaveDrag(iEdge) = topographicWaveDragCoeff*topoDragTurnOff*50.0_RKIND  * &
+                                (bed_slope_edges(iEdge)**2)                                             * &
+                                (5.24e-3_RKIND*1300.0_RKIND                                             * &


These should be defined as constant variables

dengwirda · 2024-03-31T17:41:56Z

components/mpas-ocean/src/Registry.xml

-		<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="JSL" units="unitless"
+					description="hich 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"


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?

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.

jonbob · 2024-05-02T21:59:45Z

@sbrus89 -- no problem. I'm just trying to keep PRs moving, so I also updated the bld files here

dengwirda · 2024-05-03T02:49:31Z

@jonbob, @sbrus89 and all, another option could be to just stick with the ZAE and JSL schemes --- I think they're both robust and well tested. If you're keen to merge LGF with the fixes to the array indexing as well though, I'd suggest to work through some tuning + benchmarking first.

sbrus89 · 2024-05-06T14:48:36Z

@dengwirda, I confirmed that with the LGF changes described above (e.g. MPAS-Dev/compass@ab0da31), the model achieves a 3.36 cm deep M2 RMSE with a wave drag coefficient of 2.5. That wave drag value is much more reasonable and is solidly in the rage used by ADCRIC (See figure 6 here). While I haven't done a full set of tuning runs to find an optimal value, this result gives me a lot more confidence that LGF is working correctly.

dengwirda · 2024-05-07T02:20:05Z

Looks good @sbrus89! Those array indexing issues must have been where the issues with LGF had been lurking in the past.
Thanks for chasing all of this down.

sbrus89 · 2024-05-07T13:28:20Z

@jonbob, this one is ready now when convenient.

sbrus89 · 2024-05-07T13:34:32Z

@dengwirda - no problem, I'm glad we got it sorted out. Thank you for the SAL/HA changes in this PR. I appreciate you taking the time to review as well!

@dengwirda

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

jonbob · 2024-05-16T21:32:43Z

Passes:

ERS_Ld5.T62_oQU120.CMPASO-NYF.chrysalis_intel
ERP_Ld3.ne30pg2_r05_IcoswISC30E3r5.WCYCL1850.chrysalis_intel.allactive-pioroot1

with expected NML DIFFs. Merged to next

jonbob · 2024-05-20T16:33:45Z

merged to master and expected NML DIFFs blessed

sbrus89 · 2024-05-20T17:56:36Z

Thanks very much, @jonbob!

xylar · 2024-05-29T19:39:10Z

A reminder that it is always a good idea to check the compass pr suite as part of testing changes that can affect standalone MPAS-Ocean. I'm seeing segfaults in standalone MPAS-Ocean forward runs starting with this branch.

sbrus89 · 2024-05-29T20:11:10Z

@xylar, I looked into this and wanted to note that the issue is in the namelist settings for the global_ocean cases which currently use topographic wave drag, but don't have the input stream data now required following this PR. So, the fix to the segfaults in the pr suite should be on the compass side and not in E3SM.

xylar · 2024-05-29T20:28:32Z

@sbrus89, as we discussed, that's great and quite a relief!

xylar · 2024-05-31T08:43:47Z

I'm afraid I think this needs to be fixed or reverted ASAP. I am seeing major changes in MPAS-Ocean init mode for global_ocean test cases. In particular, I am seeing zero layer thickness in the Icos240 test case for all layers >= 5.

Here is the deepest layer before this change:

And after:

xylar · 2024-05-31T09:27:58Z

components/mpas-ocean/src/mode_init/mpas_ocn_init_global_ocean.F

@@ -893,7 +893,8 @@ subroutine ocn_init_setup_global_ocean_create_model_topo(domain, iErr)!{{{
             end if
             if (isOcean) then
                ! Enforce minimum depth
-                bottomDepth(iCell) = max(bottomDepthObserved(iCell), refBottomDepth(minimum_levels))
+                bottomDepth(iCell) = max(bottomDepthObserved(iCell), config_global_ocean_minimum_depth)
+                bottomDepth(iCell) = min(bottomDepth(iCell), refBottomDepth(minimum_levels))


This should have been a max, not a min.

Thanks for fixing this @xylar! Sorry I overlooked this change and didn't catch it earlier.

A bug was introduced in E3SM-Project#6310 that made the ocean shallower rather than deeper than the minimum allowed depth. This merge fixes that bug.

… next (PR #6454) Fix bottom depth deepening in global_ocean init mode A bug was introduced in #6310 that made the ocean shallower rather than deeper than the minimum allowed depth. This merge fixes that bug. Fixes #6453 [BFB] -- mpas-ocean standalone only

…6454) Fix bottom depth deepening in global_ocean init mode A bug was introduced in #6310 that made the ocean shallower rather than deeper than the minimum allowed depth. This merge fixes that bug. Fixes #6453 [BFB] -- mpas-ocean standalone only

This merge updates the E3SM-Project submodule from [31e0924](https://github.com/E3SM-Project/E3SM/tree/31e0924) to [8939709](https://github.com/E3SM-Project/E3SM/tree/8939709). This update includes the following MPAS-Ocean and MPAS-Frameworks PRs (check mark indicates bit-for-bit with previous PR in the list): - [ ] (ocn) E3SM-Project/E3SM#6263 - [ ] (ocn) E3SM-Project/E3SM#6310 - [ ] (fwk) E3SM-Project/E3SM#6427 - [ ] (ocn) E3SM-Project/E3SM#6397 - [ ] (ocn) E3SM-Project/E3SM#6306 - [ ] (ocn) E3SM-Project/E3SM#6454

…(PR #6471) Fix OpenACC deletes for topographic wave drag Partially addresses #6470 The errors were introduced in #6310, when variables were introduced and renamed, including in the OpenACC create directives but corresponding changes were incomplete for the OpenACC directives for deleting them. [BFB]

Fix OpenACC deletes for topographic wave drag Partially addresses #6470 The errors were introduced in #6310, when variables were introduced and renamed, including in the OpenACC create directives but corresponding changes were incomplete for the OpenACC directives for deleting them. [BFB]

dengwirda and others added 16 commits March 26, 2024 12:35

Update SAL and AMHarmonicAnalysis to respect initial SSH offsets

14e19d3

- Add missing rho0gInv definition.

Update SAL and AMHarmonicAnalysis to respect initial SSH offsets

992db9b

- Clean up descriptions, comments and use real-valued switches.

Update SAL and AMHarmonicAnalysis to respect initial SSH offsets

3c32bd4

- Change sshSmoothed to btrPressure in SAL routines, since we're now working with surface pressures too.

Added two new topographic wave drag schemes.

75ebe3d

Added local generation formula topograhpic wave drag scheme

d02d740

Remove double multiplication by topographic wave drag coefficient

687c935

Fix issue where topographic wave drag was being multiplied by tendanc…

f5eaa4b

…y instead of subtracted.

Fix initialization of bottom depth in init mode

987126d

Adds 1/H term into the JSL and ZAE schemes in the topographic wave drag

d0d376d

Changed 1/H variable from layerThickness to bottomDepthEdge

b3dda60

Fixed sign of bottom depth in topographic wave drag routine

395f060

Minor changes to ZAE topographic wave drag scheme

bf001b1

updates to topographic wave drag routine to prevent NaNs

9c588a8

Removed references to unused variable bottomDepthEdgeTWD

0857b8c

Fix rebase issues with ab2

617cef2

sbrus89 added the MPAS-Ocean standalone Issues and features for standalone MPAS-Ocean code that dont impact E3SM. label Mar 27, 2024

sbrus89 requested review from mark-petersen and dengwirda March 27, 2024 22:25

sbrus89 assigned jonbob Mar 27, 2024

sbrus89 added mpas-ocean BFB PR leaves answers BFB labels Mar 27, 2024

sbrus89 mentioned this pull request Mar 27, 2024

Add a tides test case for the VR45to5 mesh MPAS-Dev/compass#802

Open

6 tasks

sbrus89 commented Mar 28, 2024

View reviewed changes

dengwirda reviewed Mar 31, 2024

View reviewed changes

jonbob added the NML label May 2, 2024

dengwirda approved these changes May 7, 2024

View reviewed changes

jonbob merged commit b25c065 into E3SM-Project:master May 20, 2024
13 checks passed

sbrus89 mentioned this pull request May 29, 2024

Turn off topographic wave drag in global ocean cases MPAS-Dev/compass#825

Merged

xylar reviewed May 31, 2024

View reviewed changes

xylar mentioned this pull request May 31, 2024

Bug in bottom depth in MPAS-Ocean global_ocean init mode #6453

Closed

xylar added a commit to xylar/E3SM that referenced this pull request May 31, 2024

Fix bottom depth deepening in global_ocean init mode

9b4046c

A bug was introduced in E3SM-Project#6310 that made the ocean shallower rather than deeper than the minimum allowed depth. This merge fixes that bug.

xylar mentioned this pull request May 31, 2024

Fix bottom depth deepening in global_ocean init mode #6454

Merged

xylar mentioned this pull request Jun 4, 2024

Update E3SM-Project submodule MPAS-Dev/compass#827

Merged

6 tasks

This was referenced Jun 13, 2024

Build fail in MPAS with nvidiagpu compiler #6470

Closed

Fix OpenACC deletes for topographic wave drag #6471

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

sbrus89 commented Mar 27, 2024 •

edited by jonbob

Loading

sbrus89 Mar 28, 2024

sbrus89 Mar 28, 2024

sbrus89 Mar 28, 2024

sbrus89 Mar 28, 2024

sbrus89 Mar 28, 2024

dengwirda Mar 31, 2024

sbrus89 Apr 1, 2024

sbrus89 Mar 28, 2024 •

edited

Loading

sbrus89 Mar 28, 2024

sbrus89 Mar 28, 2024

dengwirda Mar 31, 2024

sbrus89 Apr 1, 2024

jonbob commented May 2, 2024

dengwirda commented May 3, 2024

sbrus89 commented May 6, 2024 •

edited

Loading

dengwirda commented May 7, 2024 •

edited

Loading

sbrus89 commented May 7, 2024

sbrus89 commented May 7, 2024

jonbob commented May 16, 2024

jonbob commented May 20, 2024

sbrus89 commented May 20, 2024

xylar commented May 29, 2024

sbrus89 commented May 29, 2024

xylar commented May 29, 2024

xylar commented May 31, 2024

xylar May 31, 2024

sbrus89 May 31, 2024

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

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

Conversation

sbrus89 commented Mar 27, 2024 • edited by jonbob Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sbrus89 Mar 28, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jonbob commented May 2, 2024

dengwirda commented May 3, 2024

sbrus89 commented May 6, 2024 • edited Loading

dengwirda commented May 7, 2024 • edited Loading

sbrus89 commented May 7, 2024

sbrus89 commented May 7, 2024

jonbob commented May 16, 2024

jonbob commented May 20, 2024

sbrus89 commented May 20, 2024

xylar commented May 29, 2024

sbrus89 commented May 29, 2024

xylar commented May 29, 2024

xylar commented May 31, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sbrus89 commented Mar 27, 2024 •

edited by jonbob

Loading

sbrus89 Mar 28, 2024 •

edited

Loading

sbrus89 commented May 6, 2024 •

edited

Loading

dengwirda commented May 7, 2024 •

edited

Loading