Add facemelting to ismip6 processing script

[matthewhoffman]

This merge adds facemelting to the ismip6 output processing. If faceMeltSpeed is not in the output file, the processing will assume it was not needed. It also includes a modification requested by ISMIP6 to not mask the topg field by ice extent.

[matthewhoffman]

@trhille & @hollyhan , I added the missing code needed to include facemelting in the ISMIP6 post-processing. I have not tried it. I marked you both as reviewers, but I think a single review from whichever of you does the post-processing will suffice.

[trhille]

Testing on the first 19 time-steps in /global/cscratch1/sd/trhille/ISMIP6AE_4km_20230105/landice/ismip6_run_ais/ismip6AE/expAE01/output, we are now getting a reasonable ratio between face-melting and calving fluxes after the change in e3fc3a7:

 np.sum(faceMeltFluxArray, axis=1)
<xarray.DataArray 'calvingVelocity' (Time: 19)>
array([1.48544421, 1.51592272, 1.48098547, 1.49853168, 1.5134366 ,
       1.52476585, 1.48896772, 1.49788689, 1.49486955, 1.53907428,
       1.5355752 , 1.53331737, 1.5521225 , 1.5520087 , 1.54904589,
       1.51816453, 1.64156709, 1.56003856, 1.58241978])

np.sum(calvingFluxArray, axis=1)
<xarray.DataArray 'calvingVelocity' (Time: 19)>
array([54.0402839 , 54.17191084, 54.37512065, 54.45548776, 54.42315226,
       54.42744283, 54.1488912 , 54.34918999, 54.42410417, 53.99537311,
       53.71377185, 54.24357702, 53.97536237, 55.11158879, 54.63964945,
       54.57640848, 54.11411556, 54.55360337, 53.60329467])

 np.sum(faceMeltFluxArray, axis=1) / np.sum(calvingFluxArray, axis=1)
<xarray.DataArray 'calvingVelocity' (Time: 19)>
array([0.02748772, 0.02798356, 0.02723645, 0.02751847, 0.02780869,
       0.02801465, 0.02749766, 0.02756043, 0.02746705, 0.02850382,
       0.02858811, 0.02826726, 0.02875613, 0.0281612 , 0.02835022,
       0.02781723, 0.03033528, 0.02859643, 0.02952094])

It's a little hard to completely verify this because these fluxes are based on faceMeltSpeed and calvingVelocity, and don't explicitly account for varying thickness. However, these number are in the same ballpark as the ratio of the face-melting and calving volumes:

np.sum(faceMeltingThickness * areaCell, axis=1) / np.sum(calvingThickness * areaCell, axis=1)
array([0.03948997, 0.04677852, 0.04074631, 0.70849306, 0.04013701,
       0.04802554, 0.03989737, 0.33803736, 0.0401574 , 0.04250693,
       0.04779896, 0.04189592, 0.056934  , 0.04217208, 0.22731058,
       0.03994538, 0.05387756, 0.04491399, 0.05716555])

Note that there are large increases in the ratio of volumes at t=3, 7, and 14 (counting from 0), which also show up as increases in totalFaceMeltingFlux in globalStats.nc, while totalCalvingFlux remains fairly constant :

I don't know the precise reason for these increases, but it seems like face-melting is acting on higher thickness at those time-steps, as there is not a large increase in faceMeltFluxArray in those time-steps.

[trhille]

After further investigation, it turns out that the large spikes in face melting that we observed above were due to the deletion of stranded non-dynamic cells in Step 4 of li_apply_front_ablation_velocity: https://github.com/MALI-Dev/E3SM/blob/develop/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F#L2782-L2804

Here is a snapshot of an example area in Enderby Land from expAE01, with the white contour representing cellMask_dynamicIce immediately before face-melting is applied, and the black contour representing cellMask_ice immediately before face-melting is applied: So, while stranded non-dynamic cells like this are being removed in most time steps, when a large volume of non-dynamic ice removal corresponds with a small time step, we get a large spike in the reported totalFaceMeltingFlux, as shown in previous comments. This is technically correct, but was a bit misleading and confused me for an embarrassingly long time.

After discussion with @matthewhoffman, we decided that it makes more sense to put that mass loss into calvingThicknessFromThreshold in post-processing, since this is not physically really related to face-melting, but is just a necessary clean-up step.

After https://github.com/MPAS-Dev/MPAS-Tools/pull/490/commits/1fb836e29c902c8ce27ec681a3d3df710365ec38, the face-melting flux is consistent with expectations, with the variability from deletion of stranded non-dynamic cells in Step 4 of li_apply_front_ablation_velocity pushed into calvingThicknessFromThreshold.

np.sum(calvingThickness * areaCell / deltat_array, axis=1)
array([49228.153451  , 48798.73392818, 49032.44329166, 49214.95196977,
       49139.30843669, 49108.95946778, 48821.79868413, 48846.31918012,
       49401.00119255, 48382.79642358, 47904.21284329, 48979.80998666,
       48583.38359388, 49540.65512782, 49088.7300805 , 49470.88866793,
       48316.82709034, 48878.77000334, 47424.93068525])

np.sum(faceMeltingThicknessCleaned * areaCell / deltat_array, axis=1)
array([1944.01822659, 1873.30731363, 1887.97714191, 1920.20367756,
       1967.18496115, 1961.2663264 , 1918.65576488, 1927.41338517,
       1970.79965386, 1922.56628164, 1970.3633663 , 1973.88229605,
       1961.96986461, 2025.76840017, 1988.95688489, 1962.31699737,
       1950.14811708, 1964.74881235, 1958.0574656 ])

np.sum(calvingThicknessFromThreshold * areaCell / deltat_array, axis=1)
array([0.00000000e+00, 4.09425274e+02, 1.09914030e+02, 3.29482481e+04,
       5.12006266e+00, 3.97218190e+02, 2.92056028e+01, 1.45844674e+04,
       1.30159305e+01, 1.34037910e+02, 3.19408344e+02, 7.81719872e+01,
       8.04076717e+02, 6.34640333e+01, 9.16943069e+03, 1.38164811e+01,
       6.53044849e+02, 2.30591706e+02, 7.53014640e+02])

Also note that I had to remove all the conda initialization code from my ~/.bashrc file on Cori to get this to run, presumably due to an update in e3sm_unified.

[trhille]

Thanks for the review, @matthewhoffman. I've implemented all your suggested changes in https://github.com/MPAS-Dev/MPAS-Tools/pull/490/commits/e17b4b454b98185190f0ba94392cf2b8bc045941. Testing confirms that they do not change the answers I was getting before (see above comment):

python $mpas_tools/output_processing_li/ismip6_postprocessing/post_process_mali_to_ismip6.py -i_flux \
/global/cscratch1/sd/trhille/ISMIP6AE_4km_20230105/landice/ismip6_run_ais/ismip6AE/expAE01/output/test_ncrcat_flux.nc  \
-i_mesh /global/cscratch1/sd/trhille/ISMIP6AE_4km_20230105/landice/ismip6_run_ais/ismip6AE/expAE01/rst.2015-01-01.nc  \
--ismip6_grid_file /global/cscratch1/sd/hollyhan/ISMIP6-Projections-Forcing-2300/AIS/Ocean_forcing/ukesm1-0-ll_ssp126/1995-2300/UKESM1-0-LL_thermal_forcing_8km_x_60m_v2.nc \
--mapping_file map_AIS_4to20km_to_ismip6_8km_conserve.nc  \
-e expAE01 --mali_mesh_name AIS_4to20km -p $CSCRATCH/test_ismip6_processing/expAE01/

(Pdb) np.sum(calvingThickness * areaCell / deltat_array, axis=1)
array([49228.153451  , 48798.73392818, 49032.44329166, 49214.95196977,
       49139.30843669, 49108.95946778, 48821.79868413, 48846.31918012,
       49401.00119255, 48382.79642358, 47904.21284329, 48979.80998666,
       48583.38359388, 49540.65512782, 49088.7300805 , 49470.88866793,
       48316.82709034, 48878.77000334, 47424.93068525])
(Pdb) np.sum(faceMeltingThicknessCleaned * areaCell / deltat_array, axis=1)
array([1944.01822659, 1873.30731363, 1887.97714191, 1920.20367756,
       1967.18496115, 1961.2663264 , 1918.65576488, 1927.41338517,
       1970.79965386, 1922.56628164, 1970.3633663 , 1973.88229605,
       1961.96986461, 2025.76840017, 1988.95688489, 1962.31699737,
       1950.14811708, 1964.74881235, 1958.0574656 ])
(Pdb) np.sum(calvingThicknessFromThreshold * areaCell / deltat_array, axis=1)
array([0.00000000e+00, 4.09425274e+02, 1.09914030e+02, 3.29482481e+04,
       5.12006266e+00, 3.97218190e+02, 2.92056028e+01, 1.45844674e+04,
       1.30159305e+01, 1.34037910e+02, 3.19408344e+02, 7.81719872e+01,
       8.04076717e+02, 6.34640333e+01, 9.16943069e+03, 1.38164811e+01,
       6.53044849e+02, 2.30591706e+02, 7.53014640e+02])