Updates to Redi scheme #117

vanroekel · 2025-01-06T16:14:29Z

This addresses the need to transition the Redi isopycnal mixing scheme from isopycnal at depth to pure horizontal mixing in the surface boundary layer (see, e.g. - https://journals.ametsoc.org/view/journals/clim/21/6/2007jcli1508.1.xml). Currently the default behavior of master is to taper redi to zero in the BLD. To accomplish this, the RediKappa array is now 3d to allow for horizontal variability of the boundary layer and mixed layer depths. The transition is linear between the density threshold mixed layer depth and the boundary layer depth. There are places when the MLD is shallower than the BLD. In these cases, the MLD is set to BLD + a config flag option (config_redi_transition_layer_mld_offset).

The temperature and salinity gradients are also converted from model index space to z-coordinates. This addresses #113 - this part can easily be separated out in a final E3SM pull request.

Redi is set to GM and the taper is fixed to be consistent with Danabasoglu and Marshall 2007

There are cases where the density threshold MLD < BLD. In these cases, the MLD for the transition layer in the Redi scheme is set to BLD + the config flag value

vanroekel · 2025-01-06T16:16:28Z

Initial results from a G-case

This PR - https://web2.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.vanroekel/icos_rediks_yr_11_20/
Baseline - https://web2.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.vanroekel/icos_base_yr_11_20/

vanroekel · 2025-01-06T16:16:55Z

@jonbob - could you please run the namelist scripts for the new config option on this branch when you have a minute?

vanroekel · 2025-01-06T16:18:22Z

also pinging @maltrud, @irenavankova, and @proteanplanet as an FYI.

jonbob · 2025-01-06T17:18:11Z

OK, I had to fix some text in Registry that caused the xml parser to fail -- so please make sure you're OK with it. Otherwise I ran the scripts to make the bld files consistent and pushed the updates to your branch

vanroekel · 2025-01-06T18:58:16Z

Thanks @jonbob - I think the changes are okay, but may change lt to shallower than in another commit

jonbob · 2025-01-06T19:04:59Z

Thanks @vanroekel -- that makes sense. We'll have to change it to match in the namelist definitions file

irenavankova · 2025-01-08T20:20:39Z

I did some tests with Luke's changes when I was looking into GM tuning for polar SORRM configuration, focusing on Amundsen and Amery regions. Because the changes are quite impactful I want to bring them up here.

What I found is that there are some quite significant changes in ice-shelf cavities and around the continental shelves of Antarctica. Specifically mixed layer deepens a lot in the cavities and also around the coast. With the N2_dependent option there is a very high increase in melt rates, especially at the large cold cavities. I think these high melt rates are not driven by intrusions at depth (if anything sea-floor continental shelf temperatures cool slightly), but by a decrease in locally formed cold dense water. As a result we essentially stop producing Ice Shelf Water as indicated by the absence of basal freezing. I can recover the lower ice shelf melt rates, but for that need very low GM/Redi kappas and then the Antarctic slope front looks unrealistic compared to observations.

I just made a confluence page on those simulations here, and posted links to analysis etc there:

https://acme-climate.atlassian.net/wiki/spaces/PSC/pages/4879155201/Redi+scheme+changes+tested+in+SORRM+G-case

I wanted to point out that this is with a V2 mesh, and a bit older code base - I was using a reference case from August for comparisons that I had sufficiently spun up. I merged Luke's changes with the old branch and there were no conflicts.

proteanplanet · 2025-01-08T23:18:16Z

I suggested to @irenavankova that she confirm this problem in the SORRM V3 alfred1 B-case configuration.

Tagging @darincomeau and @cbegeman

vanroekel · 2025-01-09T04:26:22Z

@irenavankova thanks for the analysis, could I ask a few questions?

if you do rerun with the new config, can you please use this updated branch? The previous code had a number of issues, including ERS failures.
In this new branch, there is a new parameter that sets how thick the transition layer is, previously I hard coded it to 100m, it is a parameter now. I've not tested sensitivity to this. It could impact ML deepening.
Have you ever looked at boundary layer depth vs mixed layer depth in the cavities? I wonder if the boundary layer depth from KPP is an appropriate quantity in cavities. Perhaps this transition is not good for cavities and we should think on alternatives.
in your analysis are you concluding that constant vs. N2_dependent have the same melt rate increases and MLD deepening?
noting your strong horizontal variability in MLD in the reference case, it is possible there is significant horizontal entrainment and overall deepening of MLD across the cavity. Another possibility to consider - we could horizontally smooth MLD before computing transition layer depths, or perhaps better we have a second surface taper option that fixes a sfc depth and transition layer, which may be better in cavities.

Thinking about it there may need to be a different rotation of the flux to dampen the cross mixed layer transport. I'd be happy to talk about this further with you or any others on the polar team

irenavankova · 2025-01-09T16:02:27Z

I will discuss with Luke tomorrow and coordinate more tests afterwards.

vanroekel · 2025-01-09T16:39:31Z

In pondering this further, I want to make some changes to this branch before more testing commences. I'll try do that today. @irenavankova I'll show you the changes tomorrow.

1. allows redi to be constant when gm=N2_dependent 2. improves the new boundary layer taper for less entrainment 3. reinstates original N2_dependent as a new option 'N2_dependent_orig'

irenavankova · 2025-01-15T23:50:46Z

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

+           !$omp do schedule(runtime) private(k)
+           do iEdge = 1,nEdges
+              do k = minLevelEdgeBot(iEdge),maxLevelEdgeTop(iEdge)
+                 RediKappa(k,iEdge) = gmKappaScaling(k,iEdge)*RediKappa(k,iEdge)


I am probably missing something here, but wanted to clarify anyway. As of now it seems that RediKappa is not getting multiplied by RediHorizontalTaper when config_Redi_closure == 'N2_dependent'. Also, it is not clear to me why RediKappa needs to be multiplied by gmKappaScaling in the config_Redi_closure == 'N2_dependent' case - they were both calculated above and are very similar, I thought that calculation (when local_config_GM_kappa_lat_depth_variable) already introduced the N2 dependency into RediKappa.

Good catch. Thanks. The lines above where it is calculated need to be removed. With two versions of N2_dependent I didn't want to calculate the Redi twice, so in this way, RediKappa gets the original N2_dependent scaling or the new one depending on the choices, but I forgot to delete the one above

irenavankova · 2025-01-16T00:02:44Z

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

-               RediKappa(iEdge) = gmHorizontalTaper(iEdge) * gmBolusKappa(iEdge)
+               do k = minLevelEdgeBot(iEdge),maxLevelEdgeTop(iEdge)
+                  RediKappa(k,iEdge) = gmHorizontalTaper(iEdge) * gmBolusKappa(iEdge)
+               end do
            end do
            !$omp end do
            !$omp end parallel
         else if (config_Redi_closure == 'constant'.and. &
                   config_Redi_horizontal_taper == 'RossbyRadius') then


Why is this done only for RossbyRadius and not for any RediHorizontalTaper?

sorry I'm not sure what you are pointing at here. Can you explain what you mean by 'this'? Maybe you mean why there is no 'ramp' here? If so, that's applied in mpas_ocn_tracer_hmix_redi.F toward the bottom (in the init routine). The ramp taper is constant so it's applied once. However, if we move to a separate N2 dependent that varies in time we do need to apply the taper every time and can't assume it's already in RediKappa. I'll fix that. If that's not what you mean here, let me know

Let me know what you think of the fix I've proposed.

Yes, that is what I was refering to, I am not seeing where RediKappa is multiplied by RediHorizontalTaper in a situation where RediHorizontalTaper = 'ramp'. But yes, I see that in the constant case neither field is evolving so that case can be only done in init. I was also thinking more generally though, when we introduce another taper - like the one from the Dettling paper (slope dependent) that this could be handled in the same place for all tapers, but maybe that is not the most efficient way computationally.

thanks for the clarification. I think what I changed will do what you need. See line 765. But it is just if N2_dependent is chosen, maybe better is

config_Redi_closure .ne. 'constant' .or. config_Redi_closure .ne. 'equalGM'

what do you think?

vanroekel · 2025-01-16T04:30:22Z

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

+                  limiterUp(k,iCellSelf,iEdge) = slopeTaperUp !1.0_RKIND !(1.0_RKIND - RediKappaSfcTaper(k,iEdge)* &
+                                   ! (1.0_RKIND - slopeTaperUp))
+                  limiterDown(k,iCellSelf,iEdge) = slopeTaperDown! 1.0_RKIND !(1.0_RKIND - RediKappaSfcTaper(k,iEdge))!* &
+                                   ! (1.0_RKIND - slopeTaperDown))


@irenavankova these lines are important to pick up in your runs.

vanroekel · 2025-01-16T04:35:56Z

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

+           !$omp do schedule(runtime) private(k)
+           do iEdge = 1,nEdges
+              do k = minLevelEdgeBot(iEdge),maxLevelEdgeTop(iEdge)
+                 RediKappa(k,iEdge) = RediHorizontalTaper(iEdge)*gmKappaScaling(k,iEdge)*RediKappa(k,iEdge)


I added the horizontal taper here as it will be the ramp, one, or rossby radius depending on parameters

irenavankova · 2025-01-16T18:05:32Z

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

+            cell2 = cellsOnEdge(2, iEdge)
+
+            bldEdge = 0.5_RKIND*(boundaryLayerDepth(cell1) + boundaryLayerDepth(cell2))
+            mldEdge = config_redi_transition_layer_mld_offset*bldEdge


I don't think mldEdge is used in the code anymore, so can be deleted?
But if it were - then with the default value of 100 m it would be very large - is that intended?

Yes you are right, I rewrote to not use it, but I think we should keep it for now. We may want to consider having this approach in cavities and a different one globally. And yes 100m is very large that was the default I tried, but we'd probably want something closer to a couple of grid cells (20-50m)

ok, so can we rename config_redi_transition_layer_mld_offset in Registry to something like: The fraction of boundary layer depth where RediKappaSfcTaper begins to transition linearly from 0 to 1, where 1 is reached at boundary layer depth.
And also change units to unitless and default to something between 0 and 1?

Yes that's right, I'll do that when we settle on the final form

cbegeman · 2025-01-17T17:23:17Z

@vanroekel I wanted to double check with you about where in vertical space you evaluate dTdx and dSdx. I would think you want to evaluate this quantity at the top of cells since that's where the vertical diffusivity is defined.

My reading of the code is that you are doing:

dTdx|top = dTdx|center + dTdz|top * dzdx|center
dzdx|center = (zMid(k,cell2) - zMid(k,cell1)) / dx

And I'm wondering if it should be:

dTdx|top = dTdx|center + dTdz|top * dzdx|top
dzdx|top = dzdx|center - 0.5 * dhdx
dhdx = (layerThickness(k,cell2) - layerThickness(k,cell1)) / dx

In other words, dzdx would be the slope of the top of the layer.

This comment perhaps belong here #113.

vanroekel · 2025-01-18T04:23:06Z

Good question @cbegeman. I think what I have isn't quite right, but I think it also needs to be different from what you have as well. the dTdx lives on the edge in the middle of the cell (as defined by Griffies 1998 - https://journals.ametsoc.org/view/journals/phoc/28/5/1520-0485_1998_028_0805_idiazc_2.0.co_2.xml). I would cast this as

dTdx_z_center = dTdx_index_space_center + dTdz_center*dzdx_center

This is a correction from index space to physical depth space, so not a movement to the top. Diffusivity is at the top, but it's a combination of the slope which is formed by dTdx_z_center and dTdz_top (equation 34 of griffies 1998)

But right now, I have dTdz_top and dzdx_top (zmid(k) - zmid(k+1)) I think would live at the cell top, but those should be in the centers.

Does that make sense?

cbegeman · 2025-01-21T17:56:11Z

@vanroekel Yes, what you proposed looks consistent with evaluating dTdz at the centers. Thanks!

vanroekel added 7 commits December 16, 2024 23:29

Makes redi and GM consistent for N2_dependent option

a75665e

Redi is set to GM and the taper is fixed to be consistent with Danabasoglu and Marshall 2007

few fixes

becbd7c

Fixes equalGM interface

2c1b3d0

Removes unecessary features

43f4bed

fixes ERS and PEM failures

786d0fd

Fixes PET failure

d613cb4

Adds config flag for MLD offset

2c6bcc5

There are cases where the density threshold MLD < BLD. In these cases, the MLD for the transition layer in the Redi scheme is set to BLD + the config flag value

vanroekel added the enhancement New feature or request label Jan 6, 2025

vanroekel requested review from mark-petersen and alicebarthel January 6, 2025 16:14

jonbob added 2 commits January 6, 2025 11:15

Fix text that caused the xml parser to fail

854a99d

Make bld files consistent with Registry

86a7d3b

vanroekel added 2 commits January 10, 2025 16:12

New options

0c47b5c

1. allows redi to be constant when gm=N2_dependent 2. improves the new boundary layer taper for less entrainment 3. reinstates original N2_dependent as a new option 'N2_dependent_orig'

fixes limiter bug

a1b3cc4

irenavankova reviewed Jan 15, 2025

View reviewed changes

irenavankova reviewed Jan 16, 2025

View reviewed changes

Fixes to the redi taper for N2_dependent

0731a18

vanroekel commented Jan 16, 2025

View reviewed changes

irenavankova reviewed Jan 16, 2025

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Updates to Redi scheme #117

Updates to Redi scheme #117

vanroekel commented Jan 6, 2025 •

edited

Loading

vanroekel commented Jan 6, 2025

vanroekel commented Jan 6, 2025

vanroekel commented Jan 6, 2025

jonbob commented Jan 6, 2025

vanroekel commented Jan 6, 2025 •

edited

Loading

jonbob commented Jan 6, 2025

irenavankova commented Jan 8, 2025

proteanplanet commented Jan 8, 2025

vanroekel commented Jan 9, 2025

irenavankova commented Jan 9, 2025

vanroekel commented Jan 9, 2025

irenavankova Jan 15, 2025

vanroekel Jan 16, 2025

irenavankova Jan 16, 2025

vanroekel Jan 16, 2025

vanroekel Jan 16, 2025

irenavankova Jan 16, 2025

vanroekel Jan 16, 2025

vanroekel Jan 16, 2025

vanroekel Jan 16, 2025

irenavankova Jan 16, 2025

vanroekel Jan 16, 2025

irenavankova Jan 16, 2025

vanroekel Jan 16, 2025

cbegeman commented Jan 17, 2025

vanroekel commented Jan 18, 2025

cbegeman commented Jan 21, 2025 •

edited

Loading

Updates to Redi scheme #117

Are you sure you want to change the base?

Updates to Redi scheme #117

Conversation

vanroekel commented Jan 6, 2025 • edited Loading

vanroekel commented Jan 6, 2025

vanroekel commented Jan 6, 2025

vanroekel commented Jan 6, 2025

jonbob commented Jan 6, 2025

vanroekel commented Jan 6, 2025 • edited Loading

jonbob commented Jan 6, 2025

irenavankova commented Jan 8, 2025

proteanplanet commented Jan 8, 2025

vanroekel commented Jan 9, 2025

irenavankova commented Jan 9, 2025

vanroekel commented Jan 9, 2025

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

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

cbegeman commented Jan 17, 2025

vanroekel commented Jan 18, 2025

cbegeman commented Jan 21, 2025 • edited Loading

vanroekel commented Jan 6, 2025 •

edited

Loading

vanroekel commented Jan 6, 2025 •

edited

Loading

cbegeman commented Jan 21, 2025 •

edited

Loading