Add cryo terms to coupler budget #74

darincomeau · 2024-01-29T16:54:59Z

Adds five new o2x coupling fields and terms to coupler budget for polar configurations:

wism - water from ice shelf basal melting (either prognostic or data)
wrrof - water removed from Antarctica liquid runoff
wriof - water removed from Antarctica solid runoff

hism - heat from ice shelf basal melting
hriof - heat from removed ice runoff

darincomeau · 2024-01-29T16:59:50Z

Results from the last month of one year tests.

ne30pg2_ECwISC30to60E2r1.WCYCL1850:

(seq_diag_print_mct) NET WATER BUDGET (kg/m2s*1e6): period = all_time: date =     20101     0
                           atm            lnd            rof            ocn         ice nh         ice sh            glc        *SUM*
         wfreeze     0.00000000     0.00000000     0.00000000    -0.14447987     0.07808967     0.06639020     0.00000000     0.00000000
           wmelt     0.00000000     0.00000000     0.00000000     0.56578336    -0.56105289    -0.00474712     0.00000000    -0.00001666
           wrain   -32.31931275     6.58372094     0.00000000    25.67277765     0.04991955     0.01230397     0.00000000    -0.00059064
           wsnow    -2.13925165     0.83594827     0.00000000     0.95016920     0.13399754     0.21911914     0.00000000    -0.00001750
           wevap    34.46006843    -4.38902620     0.00000000   -30.04343738    -0.01192597    -0.01577898     0.00000000    -0.00010010
         wrunoff     0.00000000    -0.60896813     0.06726314     0.54133451     0.00000000     0.00000000     0.00000000    -0.00037048
         wfrzrof     0.00000000    -0.20327073     0.00000000     0.20325413     0.00000000     0.00000000     0.00000000    -0.00001660
           wberg     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000
            wism     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000
           wrrof     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000
           wriof     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000
          wirrig     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000
           *SUM*     0.00150403     2.21840414     0.06726314    -2.25459840    -0.31097210     0.27728721     0.00000000    -0.00111198

ne30pg2_ECwISC30to60E2r1.CRYO1850-DISMF:

(seq_diag_print_mct) NET WATER BUDGET (kg/m2s*1e6): period = all_time: date =     20101     0
                           atm            lnd            rof            ocn         ice nh         ice sh            glc        *SUM*
         wfreeze     0.00000000     0.00000000     0.00000000    -0.10856025     0.07964831     0.02891194     0.00000000    -0.00000000
           wmelt     0.00000000     0.00000000     0.00000000     0.56418650    -0.54820973    -0.01596599     0.00000000     0.00001078
           wrain   -32.10275677     6.38852996     0.00000000    25.65772076     0.04453218     0.01147925     0.00000000    -0.00049462
           wsnow    -2.13062436     0.83462355     0.00000000     0.95319482     0.13564092     0.20717841     0.00000000     0.00001335
           wevap    34.22371427    -4.36568440     0.00000000   -29.83287799    -0.01057680    -0.01467655     0.00000000    -0.00010146
         wrunoff     0.00000000    -0.59464619     0.06129152     0.53293445     0.00000000     0.00000000     0.00000000    -0.00042022
         wfrzrof     0.00000000    -0.21128857     0.00000000     0.21127484     0.00000000     0.00000000     0.00000000    -0.00001372
           wberg     0.00000000     0.00000000     0.00000000     0.07937360     0.00000000     0.00000000     0.00000000     0.07937360
            wism     0.00000000     0.00000000     0.00000000     0.03377146     0.00000000     0.00000000     0.00000000     0.03377146
           wrrof     0.00000000     0.00000000     0.00000000    -0.01411851     0.00000000     0.00000000     0.00000000    -0.01411851
           wriof     0.00000000     0.00000000     0.00000000    -0.17140508     0.00000000     0.00000000     0.00000000    -0.17140508
          wirrig     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000
           *SUM*    -0.00966686     2.05153436     0.06129152    -2.09450538    -0.29896512     0.21692706     0.00000000    -0.07338442

ne30pg2_ECwISC30to60E2r1.CRYO1850:

(seq_diag_print_mct) NET WATER BUDGET (kg/m2s*1e6): period = all_time: date =     20101     0
                           atm            lnd            rof            ocn         ice nh         ice sh            glc        *SUM*
         wfreeze     0.00000000     0.00000000     0.00000000    -0.11310092     0.08161947     0.03148145     0.00000000    -0.00000000
           wmelt     0.00000000     0.00000000     0.00000000     0.61626948    -0.58662693    -0.02965021     0.00000000    -0.00000766
           wrain   -32.14911941     6.56787763     0.00000000    25.51340667     0.05426512     0.01304402     0.00000000    -0.00052597
           wsnow    -2.16887503     0.84515116     0.00000000     0.96718653     0.12969461     0.22683606     0.00000000    -0.00000666
           wevap    34.32869776    -4.37683964     0.00000000   -29.92584677    -0.01067165    -0.01544155     0.00000000    -0.00010186
         wrunoff     0.00000000    -0.60071222     0.06399359     0.53630126     0.00000000     0.00000000     0.00000000    -0.00041738
         wfrzrof     0.00000000    -0.20471793     0.00000000     0.20465379     0.00000000     0.00000000     0.00000000    -0.00006414
           wberg     0.00000000     0.00000000     0.00000000     0.07937360     0.00000000     0.00000000     0.00000000     0.07937360
            wism     0.00000000     0.00000000     0.00000000     0.07718944     0.00000000     0.00000000     0.00000000     0.07718944
           wrrof     0.00000000     0.00000000     0.00000000    -0.01619780     0.00000000     0.00000000     0.00000000    -0.01619780
           wriof     0.00000000     0.00000000     0.00000000    -0.16214656     0.00000000     0.00000000     0.00000000    -0.16214656
          wirrig     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000
           *SUM*     0.01070332     2.23075900     0.06399359    -2.22291127    -0.33171938     0.22626976     0.00000000    -0.02290498

darincomeau · 2024-01-29T17:03:09Z

driver-mct/main/seq_diag_mct.F90

       end if

       lSize = mct_avect_lSize(o2x_o)
       ic = c_ocn_or
       do n=1,lSize
          ca_o =  dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ko,n)
          ca_i =  dom_o%data%rAttr(kArea,n) * frac_o%rAttr(ki,n)
+          ca_c =  dom_o%data%rAttr(kArea,n) !DC area including ice-shelf cavities


Adding this as a scaling factor was a bit of guesswork, but seems to work. Without it, I was getting 0 for ice-shelf melt in the coupler table, I think because the coupler is seeing the ice-shelf cavity grid cells as land, so 0 area for ocean and ice.

@xylar here is where I got around the coupler masking out ocean cells under ice shelves.

Got it. Thanks!

darincomeau · 2024-01-29T17:04:33Z

driver-mct/main/seq_diag_mct.F90

@@ -1652,7 +1673,8 @@ subroutine seq_diag_ice_mct( ice, frac_i, infodata, do_i2x, do_x2i)
          nf = f_hsen  ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_i*i2x_i%rAttr(index_i2x_Faii_sen,n)
          nf = f_hberg ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - (ca_o+ca_i)*i2x_i%rAttr(index_i2x_Fioi_bergh,n)
          nf = f_wmelt ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_i*i2x_i%rAttr(index_i2x_Fioi_meltw,n)
-          nf = f_wberg ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - (ca_o+ca_i)*i2x_i%rAttr(index_i2x_Fioi_bergw,n)
+!DC propose removing wberg from ice entry


Note I've removed the wberg calculation for ice. We don't want this term to be cancelled out by itself across ice/ocn since it's an additional forcing to the system. Rather we want to see how closely it's balanced by the removed runoff term(s).

darincomeau · 2024-01-29T17:06:26Z

components/mpas-ocean/driver/mpaso_cpl_indices.F

@@ -187,6 +189,10 @@ subroutine mpaso_cpl_indices_set( )
    index_o2x_Faoo_fco2_ocn = mct_avect_indexra(o2x,'Faoo_fco2_ocn',perrWith='quiet')
    index_o2x_Faoo_fdms_ocn = mct_avect_indexra(o2x,'Faoo_fdms_ocn',perrWith='quiet')
    index_o2x_So_ssh        = mct_avect_indexra(o2x,'So_ssh')
+!DC P needed since in ice-shelf ocean domain?


Noting I added 'P' for the ice-shelf melt term; I'm not sure if this matters with the way I did the scaling below. The removed runoff terms do not have P.

I just tested without the 'P' and get the same answers, so I'm removing this.

darincomeau · 2024-01-29T17:09:13Z

I've posted a single month of the water coupler budget, but for the polar runs, just noting that looking at a single month can be a bit misleading, since the ice-shelf/berg melt has the opposite seasonal cycle as the removed runoff terms.

Ideally we would like to see, over the course of the year, that (wism + wberg) \approx (wrrof + wriof).

EDIT: the tables above are the all-time sums over the 1 year simulations.

xylar · 2024-01-29T18:21:09Z

@darincomeau, I don't really speak coupler yet (though I'm keen to get on Duolingo or wherever and learn some!). Is the coupler budget output just the accumulated change in whatever field over the given period (e.g. a month?), as opposed to accumulated over the simulation?

xylar · 2024-01-29T18:21:57Z

By the way, thanks so much for getting this going! I'm realizing how far over my head this would have been.

darincomeau · 2024-01-29T18:40:09Z

Is the coupler budget output just the accumulated change in whatever field over the given period (e.g. a month?), as opposed to accumulated over the simulation?

Good question; I thought they were monthly sums, but given it says 'period = all_time:', perhaps they are accumulations over the entire simulation, which would be preferable.

jonbob · 2024-01-29T18:46:20Z

@xylar and @darincomeau -- by default, the coupler budgets are output monthly, annually, and all-time. But this is controlled by flags which also allow daily and instantaneous options

xylar · 2024-01-29T18:48:42Z

@jonbob, I guess my question was whether the field going in should be removedRiverRunoffFlux, etc. as @darincomeau has it or accumulatedRemovedRiverRunoffFlux, etc. from the analysis member that adds those up in MPAS-Ocean. It sounds like @darincomeau is probably doing things right and the coupler itself takes care of accumulating everything over the desired period?

darincomeau · 2024-01-29T18:51:26Z

It sounds like @darincomeau is probably doing things right and the coupler itself takes care of accumulating everything over the desired period?

I was also wondering about accumulation, in particular when the mpas-ocean timestep is shorter than the coupling interval.

xylar · 2024-01-29T18:53:45Z

Exactly!

xylar · 2024-01-29T18:54:28Z

@darincomeau, I think we probably need to at the very lease accumulate the field over the coupling interval like happens with other fields passed to the coupler.

darincomeau · 2024-01-29T18:55:37Z

Ok, I'll try to find an analog, but if you have one in mind, I'll take it as a starting point.

xylar · 2024-01-29T18:56:20Z

https://github.com/E3SM-Project/E3SM/blob/master/components/mpas-ocean/src/shared/mpas_ocn_time_average_coupled.F

jonbob · 2024-01-29T18:56:23Z

Sorry if I misunderstood. But you are correct -- the coupler will take care of accumulations in general but the ocean needs to do its own accumulation over the coupling interval

xylar · 2024-01-29T18:56:59Z

@jonbob, thanks, that clears things up.

darincomeau · 2024-01-29T18:58:06Z

https://github.com/E3SM-Project/E3SM/blob/master/components/mpas-ocean/src/shared/mpas_ocn_time_average_coupled.F

Thanks @xylar - looks like this is an important step I missed!

xylar · 2024-01-29T18:59:26Z

Maybe use avgTotalFreshWaterTemperatureFlux as your field to imitate?

https://github.com/E3SM-Project/E3SM/blob/3e727e0fbbae951added95895482e5a0da36069a/components/mpas-ocean/src/Registry.xml#L3689-L3691

darincomeau · 2024-02-01T02:02:29Z

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

+              call mpas_pool_get_array(forcingPool, 'landIceFreshwaterFlux', landIceFreshwaterFlux)
+           endif
+
+          ! load data melt rates if used


this didn't work, stays 0 when data ISMF is used.

You should be able to use landIceFreshwaterFlux for both standalone and data, since one is just copied to the other:
https://github.com/E3SM-Project/E3SM/blob/104a2027d9ad34fbd9198422641ad47f19db65c0/components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F#L300
But that doesn't explain why you're getting zero.

darincomeau · 2024-02-01T02:06:23Z

With that last commit, the results for ne30pg2_ECwISC30to60E2r1.CRYO1850 are exactly the same. I then realized for that resolution the ocean timestep is the same as the coupling interval. So in a way I guess that's good confirmation, but I'll test with a SORRM mesh for a case where the ocean timestep is less than the coupling interval.
Data ISMF are not being accounted for however, I'm not sure if the variable is not yet assigned when trying to be used in mpas_ocn_time_averaged.F, or something else.

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

darincomeau · 2024-02-08T02:24:55Z

@mark-petersen can you check the proposed unit changes in the conservation check analysis member's Registry here 285be91?

@cbegeman and I looked into this and see the flux variables being multiplied by m^2 here:
https://github.com/E3SM-Project/E3SM/blob/50485d0c929a873b5935ea3c2dc1da4d236319dd/components/mpas-ocean/src/analysis_members/mpas_ocn_conservation_check.F#L954
and not divided out again. Also when we look at the values that are written out to the conservation check analysis member, we can back out values in kg m^-2 s^-1 that appear in the ocn.log file from this analysis member by scaling by 1/area again.

The reason we're looking here is we want to use these accumulated*Flux quantities to produce a net AIS runoff time series plot for polar simulations to compare to global ssh anomaly, and to us everything is consistent if these quantities are actually in units kg s^-1.

Also, many thanks @mark-petersen for your work on this analysis member in the past. It's made verification of this coupler budget accounting much easier.

xylar · 2024-02-08T16:15:04Z

@darincomeau, thanks for the progress yesterday! I haven't managed to do any more debugging yet today. I'll see what I can manage this evening. @darincomeau, if I do a merge with E3SM-Project#6221, is this branch otherwise up-to-date with what you're testing with?

darincomeau · 2024-02-08T16:27:16Z

@xylar yes this is up to date. Would you like me to cherry-pick that sea ice IC commit over here, or were you going to do a local merge? I'm planning on doing a clean branch for the actual E3SM PR to clean up commit history, so it's fine to add changes to files here we otherwise won't have in the final PR.

I was planning to start the heat terms today, but I'll keep that local until it's complete.

xylar · 2024-02-08T23:12:27Z

@darincomeau, I'm still debugging, writing out the streams I recommended above.

One thought: could it be that the coupler is masking out the ice shelves when it sums up the fluxes? We certainly mask out ice shelves in some fields (e.g. those that get remapped to the atmosphere or those that get sent to sea ice). Of course, that wouldn't explain why we're not getting zeros for prognostic melt.

xylar · 2024-02-08T23:16:56Z

I can confirm that the variable avgLandIceFreshwaterFlux has non-zero values so the problem has to be on the coupler side of things.

darincomeau · 2024-02-08T23:23:09Z

One thought: could it be that the coupler is masking out the ice shelves when it sums up the fluxes? We certainly mask out ice shelves in some fields (e.g. those that get remapped to the atmosphere or those that get sent to sea ice). Of course, that wouldn't explain why we're not getting zeros for prognostic melt.

@xylar I don't think this is what's going on, though this issue did come up. I commented in the code, I'll ping you where I made a change to get the coupler to not mask out ice shelves. It seems to be doing the right thing, but it was a bit of guesswork.
Also, in the original implementation where I was sending the data ISMF directly from ocn_comp_mct.F, without averaging over coupling interval, DISMF did show up in the coupler table (in a table at the top of the PR).

darincomeau · 2024-02-08T23:25:40Z

I can confirm that the variable avgLandIceFreshwaterFlux has non-zero values

Right, I saw this too, it ends up non-zero in timeSeriesStatsMonthly, as well as in the ocean conservation analysis member. That led me to think it was something specific to the mpas_ocn_time_average.F routine.

darincomeau · 2024-02-08T23:27:07Z

@xylar it's late, get some rest! I haven't had a chance to do anything here today, maybe when I implement the heat terms this evening, something will magically pop out.

components/mpas-ocean/driver/ocn_comp_mct.F

darincomeau · 2024-02-09T04:24:14Z

(seq_diag_print_mct) NET HEAT BUDGET (W/m2): period =  monthly: date =     10201     0
                           atm            lnd            rof            ocn         ice nh         ice sh            glc        *SUM*
         hfreeze     0.00000000     0.00000000     0.00000000     0.06212584    -0.06212584     0.00000000     0.00000000     0.00000000
           hmelt     0.00000000     0.00000000     0.00000000    -0.41618844     0.33142624     0.08507855     0.00000000     0.00031635
          hnetsw  -172.11170145    35.69859350     0.00000000   136.29059592     0.06204593     0.06810322     0.00000000     0.00763712
           hlwdn  -332.39451972    80.08749536     0.00000000   247.09087671     5.01251824     0.20517240     0.00000000     0.00154300
           hlwup   388.84155761   -98.86623025     0.00000000  -283.46553532    -6.27045094    -0.23999653     0.00000000    -0.00065543
         hlatvap    83.00674983    -6.87900876     0.00000000   -76.09764783    -0.02550562    -0.00479259     0.00000000    -0.00020497
         hlatfus     0.67328674    -0.35937298     0.00000000    -0.24544797    -0.06723218    -0.00120155     0.00000000     0.00003206
          hiroff     0.00000000     0.04587322     0.00000000    -0.04559532     0.00000000     0.00000000     0.00000000     0.00027791
            hsen    19.75339851    -8.47541141     0.00000000   -11.26688638    -0.01226106     0.00067890     0.00000000    -0.00048143
           hberg     0.00000000     0.00000000     0.00000000    -0.06320892     0.00000000     0.00000000     0.00000000    -0.06320892
            hism     0.00000000     0.00000000     0.00000000     0.01125439     0.00000000     0.00000000     0.00000000     0.01125439
           hriof     0.00000000     0.00000000     0.00000000     0.03160610     0.00000000     0.00000000     0.00000000     0.03160610
        hh2otemp     0.42776435     0.00000000     0.00000000    -0.42799595     0.00000000     0.00000000     0.00000000    -0.00023160
           *SUM*   -11.80346412     1.25193868     0.00000000    11.45795284    -1.03158523     0.11304241     0.00000000    -0.01211542

(seq_diag_print_mct) NET WATER BUDGET (kg/m2s*1e6): period =  monthly: date =     10201     0
                           atm            lnd            rof            ocn         ice nh         ice sh            glc        *SUM*
         wfreeze     0.00000000     0.00000000     0.00000000    -0.18617274     0.18617274     0.00000000     0.00000000    -0.00000000
           wmelt     0.00000000     0.00000000     0.00000000    -2.45519385     2.78357026    -0.32692927     0.00000000     0.00144714
           wrain   -31.38352297     5.15651606     0.00000000    26.20419654     0.01291856     0.00159838     0.00000000    -0.00829343
           wsnow    -2.01764082     1.07693431     0.00000000     0.73553482     0.20147491     0.00360070     0.00000000    -0.00009608
           wevap    33.18177647    -2.74422912     0.00000000   -30.42688838    -0.00902837    -0.00171256     0.00000000    -0.00008195
         wrunoff     0.00000000    -0.36266438     0.17583865     0.18515025     0.00000000     0.00000000     0.00000000    -0.00167548
         wfrzrof     0.00000000    -0.13746845     0.00000000     0.13663565     0.00000000     0.00000000     0.00000000    -0.00083280
           wberg     0.00000000     0.00000000     0.00000000     0.18475440     0.00000000     0.00000000     0.00000000     0.18475440
            wism     0.00000000     0.00000000     0.00000000     0.03372607     0.00000000     0.00000000     0.00000000     0.03372607
           wrrof     0.00000000     0.00000000     0.00000000    -0.05490293     0.00000000     0.00000000     0.00000000    -0.05490293
           wriof     0.00000000     0.00000000     0.00000000    -0.09471410     0.00000000     0.00000000     0.00000000    -0.09471410
          wirrig     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000     0.00000000
           *SUM*    -0.21938732     2.98908842     0.17583865    -5.73787427     3.17510811    -0.32344276     0.00000000     0.05933083

Now with heat flux terms for ice shelf melt (hism) and removed ice runoff (hriof).

I left off heat flux due to removedRiverRunoff, which I'm not sure how to handle.

darincomeau · 2024-02-09T04:38:19Z

components/mpas-ocean/src/Registry.xml

+		<var name="avgLandIceFreshwaterFlux" type="real" dimensions="nCells Time" units="kg m^-2 s^-1"
+		         description="Time-averaged sum of freshwater fluxes associated with ice shelf basal melt fluxes cell centers sent to coupler. Positive into the ocean."
+		/>
+		<var name="avgLandIceHeatFlux" type="real" dimensions="nCells Time" units="J m^-2 s^-1"


need to check units here

darincomeau · 2024-02-09T04:40:42Z

driver-mct/shr/seq_flds_mod.F90

+    call seq_flds_add(o2x_fluxes,"Foxo_ismh")
+    longname = 'Heat flux due to basal melting of ice shelves'
+    stdname  = 'basal_iceshelf_heat_flux'
+    units    = 'J m-2 s-1'


same unit check

darincomeau · 2024-02-09T17:56:06Z

while the units were technically correct, changed to W/m^s

xylar

This looks ready for a "real" PR to me.

proteanplanet

Please proceed to a full PR with a single design document as discussed on the email thread between @stephenprice and @xylar. Please use the style of the PR E3SM-Project#5958 as a template for your work, including the link to https://acme-climate.atlassian.net/wiki/spaces/ICE/pages/3952476924/Overview+of+Icepack+physics+merge+into+MPAS-SI. This is an example of a major change to the fully coupled model that pretty much sailed through review. This may all seem like overkill, but it's important because this work could significantly impact coupler budgets across the project, so will be subject to a lot of scrutiny.

darincomeau · 2024-02-09T21:52:46Z

Thanks for pointing to that design document Andrew, I had been looking at an older version. We'll get started on that.

darincomeau · 2024-02-10T00:05:01Z

Closing in favor of E3SM-Project#6229

darincomeau added 4 commits January 29, 2024 10:50

Adding cryo coupling fields and budget terms (ismw broken)

efbd115

Changing scaling for ice-shelf cavities

10691dd

Calling in dataLandIceFreshwaterFlux also

4e16416

Commenting out wberg calculation in ice column

8295de2

darincomeau requested review from xylar, jonbob and proteanplanet January 29, 2024 16:55

darincomeau commented Jan 29, 2024

View reviewed changes

Accumulate added fluxes over coupling interval

d42b9cb

darincomeau commented Feb 1, 2024

View reviewed changes

xylar reviewed Feb 1, 2024

View reviewed changes

components/mpas-ocean/src/shared/mpas_ocn_time_average_coupled.F Outdated Show resolved Hide resolved

xylar reviewed Feb 8, 2024

View reviewed changes