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.

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

Keep glc-derived runoff as separate fields to ocn if possible #547

Open
wants to merge 1 commit into
base: main
Choose a base branch
from

Conversation

billsacks
Copy link
Member

@billsacks billsacks commented Mar 20, 2025

Description of changes

When @mvertens reworked the GLC -> OCN runoff mapping, she put in place hooks to keep the glc-derived runoff separate from the lnd-derived runoff. However, @gustavo-marques found that these are still being combined into a single merged flux. (That may have been the intent at the time, until ocean models were ready to accept this new flux. I know we had extensive conversations about this, but don't remember all of the details.)

This PR separates these two fluxes so that the ocean can separately account for the lnd vs glc-derived runoff and the associated heat fluxes. I have done this in a way so that, if we're coupled to an ocean model that doesn't expect this separate glc-derived runoff (Forr_rofl_glc and Forr_rofi_glc), then it will continue to be merged into the Foxx_rofl and Foxx_rofi fluxes. (MOM6 is set up to handle the glc-specific fluxes; I'm not sure about BLOM in NorESM.)

Specific notes

Contributors other than yourself, if any: none

CMEPS Issues Fixed (include github issue #): none

Are changes expected to change answers? YES - substantial diffs in a few fields, and expected to be roundoff-level diffs system-wide.

There are large differences in the fields ocnExp_Foxx_rofi, ocnExp_Forr_rofi_glc, ocnExp_Foxx_hrofi, and ocnExp_Foxx_hrofi_glc due to redistribution of fluxes between these different fields. Beyond that, there will be roundoff-level differences if the separation of these runoff terms is just for diagnostic purposes and the ocean model treats them identically. (The differences then just arise from order of operations differences due to having two separate fluxes in place of a single, combined flux.) In theory an ocean model could do different things with these fluxes, and that would lead to more substantial differences, but I don't think that's currently the case, at least for CESM / MOM.

Any User Interface Changes (namelist or namelist defaults changes)? No

Testing performed

Tested in the context of cesm3_0_alpha06c.

I ran a few permutations of the test SMS_Ld5.ne30pg3_t232.BLT1850.derecho_intel.allactive-decstart.

I confirmed that, with the changes here, ocnExp_Forr_rofi_glc is non-zero, and ocnExp_Foxx_rofi is 0 around Greenland, as expected.

I added hourly coupler history output by making this change:

./xmlchange HIST_OPTION=nsteps,HIST_N=1
./xmlchange --file env_test.xml HIST_OPTION=nsteps,HIST_N=1
./xmlchange ROF_NCPL=24,GLC_NCPL=24,OCN_NCPL=24

Comparing with a baseline run, the first diffs appear at time 2002-01-01-07200.

There are diffs in the following fields, as expected:

 RMS ocnExp_Forr_rofi_glc             1.2583E-07            NORMALIZED  2.0450E+01
 RMS ocnExp_Foxx_hrofi                3.8657E-04            NORMALIZED  4.2437E-01
 RMS ocnExp_Foxx_hrofi_glc            3.8657E-04            NORMALIZED  2.7026E+01
 RMS ocnExp_Foxx_rofi                 1.2583E-07            NORMALIZED  7.2900E-01

I summed ocnExp_Forr_rofi_glc+ocnExp_Foxx_rofi in the new run and the baseline and verified that, in this first time with differences, the sums are identical in the new and baseline. Similarly for ocnExp_Foxx_hrofi+ocnExp_Foxx_hrofi_glc, although there are roundoff-level differences in that sum.

In addition, there are diffs in these fields from the ocean:

 RMS ocnImp_Fioo_q                    7.8187E-13            NORMALIZED  7.9219E-18
 RMS ocnImp_So_bldepth                2.4854E-14            NORMALIZED  1.1975E-15
 RMS ocnImp_So_s                      6.7742E-16            NORMALIZED  3.3799E-17
 RMS ocnImp_So_t                      2.7349E-16            NORMALIZED  1.6198E-18

Note that those are all roundoff-level different.

After that time, differences grow, but the above analysis makes me comfortable that the changes in this PR only introduce roundoff-level diffs in the system.

I also verified that, if MOM doesn't advertise Forr_rofl_glc and Forr_rofi_glc, the behavior reverts to what it was before the changes in this PR. To do this, I needed to change some code in CMEPS related to the calculation of heat fluxes in addition to the changes in MOM. Full diffs were:

diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90
index bb7572ed7..2265a0a65 100644
--- a/config_src/drivers/nuopc_cap/mom_cap.F90
+++ b/config_src/drivers/nuopc_cap/mom_cap.F90
@@ -810,10 +810,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc)
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Sa_pslv"        , "will provide")
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofl"      , "will provide") !-> liquid runoff
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofi"      , "will provide") !-> ice runoff
-  if (cesm_coupled) then
-    call fld_list_add(fldsToOcn_num, fldsToOcn, "Forr_rofl_glc"  , "will provide") !-> liquid glc runoff
-    call fld_list_add(fldsToOcn_num, fldsToOcn, "Forr_rofi_glc"  , "will provide") !-> frozen glc runoff
-  endif
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Si_ifrac"       , "will provide") !-> ice fraction
   call fld_list_add(fldsToOcn_num, fldsToOcn, "So_duu10n"      , "will provide") !-> wind^2 at 10m
   call fld_list_add(fldsToOcn_num, fldsToOcn, "Fioi_meltw"     , "will provide")
diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90
index 180202c7e..fe243f7bf 100644
--- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90
+++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90
@@ -218,17 +218,11 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
   ! liquid glc runoff
   if ( associated(ice_ocean_boundary%lrunoff_glc) ) then
     ice_ocean_boundary%lrunoff_glc (:,:) = 0._ESMF_KIND_R8
-    call state_getimport(importState, 'Forr_rofl_glc',  &
-         isc, iec, jsc, jec, ice_ocean_boundary%lrunoff_glc, areacor=med2mod_areacor, rc=rc)
-    if (ChkErr(rc,__LINE__,u_FILE_u)) return
   endif
 
   ! frozen glc runoff
   if ( associated(ice_ocean_boundary%frunoff_glc) ) then
     ice_ocean_boundary%frunoff_glc (:,:) = 0._ESMF_KIND_R8
-    call state_getimport(importState, 'Forr_rofi_glc',  &
-         isc, iec, jsc, jec, ice_ocean_boundary%frunoff_glc, areacor=med2mod_areacor, rc=rc)
-    if (ChkErr(rc,__LINE__,u_FILE_u)) return
   endif
 
   !----

and

diff --git a/mediator/med_phases_prep_ocn_mod.F90 b/mediator/med_phases_prep_ocn_mod.F90
index e30c4ada..e5222288 100644
--- a/mediator/med_phases_prep_ocn_mod.F90
+++ b/mediator/med_phases_prep_ocn_mod.F90
@@ -171,9 +171,7 @@ contains
          FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hrofl'     , rc=rc) .and. &
          FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_rofi'      , rc=rc) .and. &
          FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hrofi'     , rc=rc) .and. &
-         FB_fldchk(is_local%wrap%FBExp(compocn), 'Forr_rofl_glc'  , rc=rc) .and. &
          FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hrofl_glc' , rc=rc) .and. &
-         FB_fldchk(is_local%wrap%FBExp(compocn), 'Forr_rofi_glc'  , rc=rc) .and. &
          FB_fldchk(is_local%wrap%FBExp(compocn), 'Foxx_hrofi_glc' , rc=rc)) then
 
        call FB_GetFldPtr(is_local%wrap%FBImp(compocn,compocn), 'So_t', tocn, rc=rc)
@@ -206,13 +204,9 @@ contains
        call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Foxx_hrofi', hrofi, rc=rc)
        if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
-       call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Forr_rofl_glc' , rofl_glc, rc=rc)
-       if (ChkErr(rc,__LINE__,u_FILE_u)) return
        call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Foxx_hrofl_glc', hrofl_glc, rc=rc)
        if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
-       call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Forr_rofi_glc' , rofi_glc, rc=rc)
-       if (ChkErr(rc,__LINE__,u_FILE_u)) return
        call FB_GetFldPtr(is_local%wrap%FBExp(compocn), 'Foxx_hrofi_glc', hrofi_glc, rc=rc)
        if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
@@ -224,8 +218,8 @@ contains
           hcond(n)  = max((tocn(n) - shr_const_tkfrz), 0._r8) * max(evap(n), 0._r8)  * shr_const_cpsw
           hrofl(n)  = max((tocn(n) - shr_const_tkfrz), 0._r8) * rofl(n)  * shr_const_cpsw
           hrofi(n)  = min((tocn(n) - shr_const_tkfrz), 0._r8) * rofi(n)  * shr_const_cpsw
-          hrofl_glc(n) = max((tocn(n) - shr_const_tkfrz), 0._r8) * rofl_glc(n)  * shr_const_cpsw
-          hrofi_glc(n) = min((tocn(n) - shr_const_tkfrz), 0._r8) * rofi_glc(n)  * shr_const_cpsw
+          hrofl_glc(n) = max((tocn(n) - shr_const_tkfrz), 0._r8) * 0._r8  * shr_const_cpsw
+          hrofi_glc(n) = min((tocn(n) - shr_const_tkfrz), 0._r8) * 0._r8  * shr_const_cpsw
        end do
 
        ! Determine enthalpy correction factor that will be added to the sensible heat flux sent to the atm

Copy link
Collaborator

@jedwards4b jedwards4b left a comment

Choose a reason for hiding this comment

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

LGTM thank you!

@billsacks
Copy link
Member Author

I went over this with @mvertens and she approves of these changes.

@jedwards4b
Copy link
Collaborator

Do you still want a review from @gustavo-marques ?

@billsacks
Copy link
Member Author

Do you still want a review from @gustavo-marques ?

Yes. I've asked Gustavo to test with these changes to confirm that this resolves the issue he was seeing. I think it makes sense to wait to merge this until he gives the go-ahead.

@gustavo-marques
Copy link
Collaborator

Thank you for this PR @billsacks
Yes, I will do some tests in G and B simulations.

@gustavo-marques
Copy link
Collaborator

Hi @billsacks,

I am trying to test this PR using the compset we have been using for the CESM3 fully-coupled development runs. This is how I created the case:

/glade/work/gmarques/cesm.sandboxes/cesm3_0_alpha06c_testing/cime/scripts/create_newcase --case b.e30_a06c.cmeps547 --compset B1850C_LTso --res ne30pg3_t232_wg37 --run-unsupported

In my sandbox (cesm3_0_alpha06c_testing), cmeps has been updated to billsacks/separate_glc_rof_fields.

I am keeping the original NCPL's:

Results in group run_coupling
	ATM_NCPL: 48
	GLC_NCPL: 1
	ICE_NCPL: 48
	LND_NCPL: 48
	NCPL_BASE_PERIOD: day
	OCN_NCPL: 24
	ROF_NCPL: 8
	WAV_NCPL: 48 

Writing CPL history files every day

./xmlchange HIST_OPTION=nsteps,HIST_N=24

However, both ocnExp_Foxx_hrofi_glc and ocnExp_Foxx_hrofl_glc are zero. Do I need to do something else to enable these fluxes?

@gustavo-marques
Copy link
Collaborator

Also, it would be great to check the "budget tables" generated by the coupled (diags.log.*). Is it possible to get these tables in short runs (5 days)?

@billsacks
Copy link
Member Author

Great questions, @gustavo-marques .

I also noticed in testing that the fluxes from GLC were 0 until the first crossing of the year boundary. To perform testing in a short test case (5 days), I set the model to start near the end of the year. These changes came from the decstart testmod that I used for testing, which does the following:

./xmlchange RUN_STARTDATE=2001-12-30
./xmlchange CLM_BLDNML_OPTS=-ignore_warnings --append

I have actually never created or looked at budget tables myself, so I'm probably not the best person to ask your question about producing budget tables. But if you remain stuck with that one, I can help look into it.

@jedwards4b
Copy link
Collaborator

@Gustavo you can control the frequency of budget output using the namelist control variables:

budget_ann = 1
 budget_daily = 0
 budget_inst = 0
 budget_ltann = 1
  budget_ltend = 0
budget_month = 1
budget_table_version = v1
do_budgets = .true.

Set in user_nl_cpl, so you want to turn on budget_daily or perhaps even budget_inst.

@gustavo-marques
Copy link
Collaborator

Thanks, @billsacks and @jedwards4b!

Bill, I couldn't get MOM6 to work with RUN_STARTDATE=2001-12-30, so I used RUN_STARTDATE=0002-01-01 instead. However, I am still seeing zeros for hrofl_glc and hrofi_glc. Do you have any ideas?

@billsacks
Copy link
Member Author

Bill, I couldn't get MOM6 to work with RUN_STARTDATE=2001-12-30, so I used RUN_STARTDATE=0002-01-01 instead. However, I am still seeing zeros for hrofl_glc and hrofi_glc. Do you have any ideas?

Just to be clear: are you seeing zeroes for all of the _glc terms - i.e., rofi_glc in addition to hrofi_glc? That's what I'm imagining... if you're seeing non-zero values for rofi_glc but zeros for hrofi_glc, that's a different issue....

I have a feeling that you need to cross the year boundary to start seeing non-zeros. So, it's not enough to just start in year 2: the run seems to need to actually cross the Dec 31 – Jan 1 boundary.

You might be able to solve this with ./xmlchange GLC_AVG_PERIOD=glc_coupling_period. We do that in short tests with CISM. I don't remember (or maybe never knew) all of the relevant details of DGLC, so I'm not sure if the same thing applies there, but you could try it.

I'm puzzled, though, why I was able to do a test with MOM starting Dec 30 but you're not. So you could also try just reproducing exactly what I did: running ./create_test SMS_Ld5.ne30pg3_t232.BLT1850.derecho_intel.allactive-decstart from cesm3_0_alpha06c.

@jedwards4b
Copy link
Collaborator

@gustavo-marques @billsacks What is the status of this PR?

@gustavo-marques
Copy link
Collaborator

By setting ./xmlchange GLC_AVG_PERIOD=glc_coupling_period, I successfully generated non-zero GLC frozen runoff within the first five days of simulation. To validate the consistency of water, heat, and salt fluxes at the ocean surface with changes in the ocean interior after incorporating this PR, I conducted short tests using G and B compsets. Results are documented in the following notebooks: g.e30_a06c.cmeps547 (G compset) and b.e30_a06c.cmeps547 (B compset).

These tests indicate that the GLC runoff fluxes are correctly transferred to the ocean model, and the ocean interior changes are consistent with the surface fluxes. In the G compset, GLC fluxes are zero. In the B compset, only frozen runoff (primarily around Greenland) is non-zero. While I am confident in approving this PR, I would like confirmation from the land ice developers (@Katetc and @gunterl) that only frozen GLC runoff should be non-zero in the B compset. Additionally, are there any scenarios where liquid GLC runoff is expected to be non-zero?

@billsacks
Copy link
Member Author

Thank you very much for doing that verification @gustavo-marques ! That's great to hear!

I would like confirmation from the land ice developers (@Katetc and @gunterl) that only frozen GLC runoff should be non-zero in the B compset. Additionally, are there any scenarios where liquid GLC runoff is expected to be non-zero?

I can answer this: For a B compset with DGLC, yes, only frozen GLC runoff should be non-zero: DGLC does not send any liquid runoff. I believe that compsets with CISM could send some (small) liquid runoff in the GLC terms.

@jedwards4b - given Gustavo's approval above, are you okay with me merging this?

@gustavo-marques
Copy link
Collaborator

Thanks, @billsacks!
For completeness, I’d like to test a case where CISM sends liquid runoff to the ocean. Could you let me know how I can set that up?

@billsacks
Copy link
Member Author

Thanks, @gustavo-marques - that seems like a great idea, and I appreciate your doing that additional testing.

It's been a while since I've set up a case with active CISM, so I may be forgetting something and/or it's possible that this won't work, but if you want to test this in a short run, I think you can do this by:

  1. In the compset long name, replace DGLC%NOEVOLVE with CISM2%GRIS-EVOLVE
  2. In user_nl_cism, set test_coupling = .true. (this changes the ice dynamics time step to 1 day instead of 1 year)
  3. In user_nl_clm, set for_testing_allow_non_annual_changes = .true. and for_testing_zero_dynbal_fluxes = .true. (these are needed when setting CISM's test_coupling = .true..

(and keep the ./xmlchange GLC_AVG_PERIOD=glc_coupling_period that you were using before).

But it's quite possible that won't work. If it doesn't, I'd ask Kate or Gunter for help, though feel free to loop me back in as well.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants