Closed trhille closed 1 year ago
This gives the same results as shown in PR 490, but also gets rid of the divide-by-zero warning.
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])
@matthewhoffman, that's true in the sense of the total run, but t=0 in the flux output files is t=1 (using zero-based indexing for both, even though one is python and one is fortran) in the model run because the diagnostic solve is not written to them. So unless we were able to append the diagnostic solve to the beginning of the flux output files, I think we want this. Since the state output files don't have all the necessary fields, that wouldn't work unless we created a bunch of all-zeros fields.
Do not allow t-1 to access last time level. Use t=0 in that case instead. This was addressed in one place during code cleanup in PR 490, but I missed it in two other places. Testing indicates it makes no difference, but I think it is technically correct.