-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathcompute_pointBudgets_history.py
330 lines (297 loc) · 16.4 KB
/
compute_pointBudgets_history.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
#!/usr/bin/env python
"""
Name: compute_pointBudgets.py
Author: Milena Veneziani
"""
# ensure plots are rendered on ICC
from __future__ import absolute_import, division, print_function, \
unicode_literals
import os
import matplotlib as mpl
mpl.use('Agg')
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import xarray as xr
import datetime
import netCDF4
from common_functions import plot_xtick_format, days_to_datetime
# Choose years
year1 = 69
year2 = 69
years = range(year1, year2+1)
referenceDate = '0001-01-01'
movingAverageMonths = 1
#movingAverageMonths = 12
# Settings for lcrc:
# NOTE: make sure to use the same mesh file that is in streams.ocean!
meshfile = '/lcrc/group/e3sm/public_html/inputdata/ocn/mpas-o/EC30to60E2r2/mpaso.EC30to60E2r2.rstFromG-anvil.201001.nc'
#casenameFull = 'v2_1.LR.historical_0101'
#casename = 'v2_1.LR.historical_0101'
#modeldir = f'/lcrc/group/e3sm/ac.golaz/E3SMv2_1/{casenameFull}/archive/ocn/hist'
casenameFull = '20230927b.branchfrom69.v3alpha04_bigrid.piControl.chrysalis'
casename = 'v3alpha04_bigrid.piControl'
modeldir = f'/lcrc/group/e3sm/ac.abarthel/E3SMv3_dev/{casenameFull}/run/debug_output/'
# Settings for nersc:
# NOTE: make sure to use the same mesh file that is in streams.ocean!
#maindir = '/global/cfs/projectdirs/e3sm'
#meshfile = f'{maindir}/inputdata/ocn/mpas-o/EC30to60E2r2/mpaso.EC30to60E2r2.rstFromG-anvil.201001.nc'
#casename = 'GM600_Redi600'
#casenameFull = 'GMPAS-JRA1p4_EC30to60E2r2_GM600_Redi600_perlmutter'
#modeldir = f'{maindir}/maltrud/archive/onHPSS/{casenameFull}/ocn/hist'
#meshfile = f'{maindir}/inputdata/ocn/mpas-o/ARRM10to60E2r1/mpaso.ARRM10to60E2r1.220730.nc'
#casenameFull = '20221201.WCYCL1950.arcticx4v1pg2_ARRM10to60E2r1.lat-dep-bd-submeso.cori-knl'
#casename = 'fullyRRM_lat-dep-bd-submeso'
#modeldir = f'/global/cscratch1/sd/milena/e3sm_scratch/cori-knl/{casenameFull}/run'
# Coordinates of point where to compute budgets
lonPoint = -10
latPoint = 67
m3ps_to_Sv = 1e-6 # m^3/s flux to Sverdrups
rho0 = 1027.0 # kg/m^3
earthRadius = 6367.44 # km
figdir = f'./volBudget/{casename}'
if not os.path.isdir(figdir):
os.makedirs(figdir)
nTime = 12*len(years)
# Read in relevant global mesh information
dsMesh = xr.open_dataset(meshfile)
nLevels = dsMesh.dims['nVertLevels']
nCells = dsMesh.dims['nCells']
# Identify index of selected ocean cell, by computing the minimum
# of the spherical distance between all points and lonPoint,latPoint
lonCell = dsMesh.lonCell
latCell = dsMesh.latCell
lonPoint = np.pi/180*lonPoint
latPoint = np.pi/180*latPoint
dlon = lonCell-lonPoint
dlat = latCell-latPoint
a = np.sin(dlat / 2) ** 2 + np.cos(latCell) * np.cos(latPoint) * np.sin(dlon / 2) ** 2
c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
spherDist = earthRadius * c
indices = xr.DataArray(data=np.arange(nCells).astype(int), dims='nCells')
iCell = indices.where(spherDist==np.min(spherDist), drop=True).values.astype(int)[0]
print(lonCell.values[iCell]*180/np.pi, latCell.values[iCell]*180/np.pi)
areaCell = dsMesh.areaCell.isel(nCells=iCell)
#
# edgeID of all edges bordering the chosen cell. If 0, edge is on land, so remove it.
edgesOnCell = dsMesh.edgesOnCell.isel(nCells=iCell).values
edgesOnCell = edgesOnCell[np.where(edgesOnCell>0)]
# for each ocean edge bordering the chosen cell, select IDs of straddling cells
cellsOnEdge = dsMesh.cellsOnEdge.isel(nEdges=edgesOnCell-1).values
coe0 = cellsOnEdge[:, 0] - 1
coe1 = cellsOnEdge[:, 1] - 1
dvEdge = dsMesh.dvEdge.isel(nEdges=edgesOnCell-1)
print(iCell)
print(coe0)
print(coe1)
# Compute volume budget
# Initialize volume budget terms
volNetLateralFlux = np.zeros(nTime)
evapFlux = np.zeros(nTime)
rainFlux = np.zeros(nTime)
snowFlux = np.zeros(nTime)
riverRunoffFlux = np.zeros(nTime)
iceRunoffFlux = np.zeros(nTime)
seaIceFreshWaterFlux = np.zeros(nTime)
layerThick = np.zeros(nTime)
frazilThick = np.zeros(nTime)
t = np.zeros(nTime)
ktime = 0
for year in years:
print(f'Year = {year:04d} out of {len(years)} years total')
#for month in range(1, 13):
for month in range(1, 2):
print(f' Month= {month:02d}')
modelfile = f'{modeldir}/{casenameFull}.mpaso.hist.am.highFrequencyOutput.{year:04d}-{month:02d}-21_00.00.00.nc'
ds = xr.open_dataset(modelfile, decode_times=False)
t[ktime] = ds.Time.isel(Time=0).values
# Compute net lateral fluxes:
if 'normalTransportVelocity' in ds.keys():
vel = ds.normalTransportVelocity.isel(Time=0, nEdges=edgesOnCell-1)
elif 'normalVelocity' in ds.keys():
vel = ds.normalVelocity.isel(Time=0, nEdges=edgesOnCell-1)
if 'normalGMBolusVelocity' in ds.keys():
vel = vel + ds.normalGMBolusVelocity.isel(Time=0, nEdges=edgesOnCell-1)
if 'normalMLEvelocity' in ds.keys():
vel = vel + ds.normalMLEvelocity.isel(Time=0, nEdges=edgesOnCell-1)
else:
raise KeyError('no appropriate normalVelocity variable found')
dzOnCells0 = ds.layerThickness.isel(Time=0, nCells=coe0)
dzOnCells1 = ds.layerThickness.isel(Time=0, nCells=coe1)
# Then, interpolate dz's onto edges, also considering the topomask
dzOnEdges = 0.5 * (dzOnCells0 + dzOnCells1)
dzOnEdges = dzOnEdges.rename({'nCells': 'nEdges'})
dArea = dvEdge * dzOnEdges
lateralFlux = (vel * dArea).sum(dim='nVertLevels', skipna=True).sum(dim='nEdges')
volNetLateralFlux[ktime] = lateralFlux.values
# Compute net surface fluxes:
if 'evaporationFlux' in ds.keys():
flux = ds.evaporationFlux.isel(Time=0, nCells=iCell)
evapFlux[ktime] = (flux * areaCell).values
else:
raise KeyError('no evaporation flux variable found')
if 'rainFlux' in ds.keys():
flux = ds.rainFlux.isel(Time=0, nCells=iCell)
rainFlux[ktime] = (flux * areaCell).values
else:
raise KeyError('no rain flux variable found')
if 'snowFlux' in ds.keys():
flux = ds.snowFlux.isel(Time=0, nCells=iCell)
snowFlux[ktime] = (flux * areaCell).values
else:
raise KeyError('no snow flux variable found')
if 'riverRunoffFlux' in ds.keys():
flux = ds.riverRunoffFlux.isel(Time=0, nCells=iCell)
riverRunoffFlux[ktime] = (flux * areaCell).values
else:
raise KeyError('no river runoff flux variable found')
if 'iceRunoffFlux' in ds.keys():
flux = ds.iceRunoffFlux.isel(Time=0, nCells=iCell)
iceRunoffFlux[ktime] = (flux * areaCell).values
else:
raise KeyError('no ice runoff flux variable found')
if 'seaIceFreshWaterFlux' in ds.keys():
flux = ds.seaIceFreshWaterFlux.isel(Time=0, nCells=iCell)
seaIceFreshWaterFlux[ktime] = (flux * areaCell).values
else:
raise KeyError('no sea ice freshwater flux variable found')
# Compute layer thickness tendencies:
if 'tendLayerThickness' in ds.keys():
layerThickTend = ds.tendLayerThickness.isel(Time=0, nCells=iCell)
layerThick[ktime] = (layerThickTend.sum(dim='nVertLevels', skipna=True) * areaCell).values
else:
raise KeyError('no layer thickness tendency variable found')
if 'frazilLayerThicknessTendency' in ds.keys():
frazilThickTend = ds.frazilLayerThicknessTendency.isel(Time=0, nCells=iCell)
frazilThick[ktime] = (frazilThickTend.sum(dim='nVertLevels', skipna=True) * areaCell).values
else:
raise KeyError('no frazil layer thickness tendency variable found')
ktime = ktime + 1
print('\nresidual')
print(m3ps_to_Sv*volNetLateralFlux[0]+1/rho0*m3ps_to_Sv*(evapFlux[0]+rainFlux[0]+snowFlux[0]+riverRunoffFlux[0]+iceRunoffFlux[0]+seaIceFreshWaterFlux[0])-m3ps_to_Sv*(layerThick[0]+frazilThick[0]))
print('\nnetlateral + evap + rain + snow + riverrunoff + icerunoff + seaiceflux, layerThickTend + frazilThickTend')
print(m3ps_to_Sv*volNetLateralFlux[0]+1/rho0*m3ps_to_Sv*(evapFlux[0]+rainFlux[0]+snowFlux[0]+riverRunoffFlux[0]+iceRunoffFlux[0]+seaIceFreshWaterFlux[0]), m3ps_to_Sv*(layerThick[0]+frazilThick[0]))
#print('\nnetlateral, evap + rain + snow, riverrunoff + icerunoff, seaiceflux, layerThickTend + frazilThickTend')
print(m3ps_to_Sv*volNetLateralFlux[0], 1/rho0*m3ps_to_Sv*(evapFlux[0]+rainFlux[0]+snowFlux[0]), 1/rho0*m3ps_to_Sv*(riverRunoffFlux[0]+iceRunoffFlux[0]), 1/rho0*m3ps_to_Sv*seaIceFreshWaterFlux[0], m3ps_to_Sv*(layerThick[0]+frazilThick[0]))
print('\nnetlateral, allSurfFluxes, layerThickTend + frazilThickTend')
print(m3ps_to_Sv*volNetLateralFlux[0], 1/rho0*m3ps_to_Sv*(evapFlux[0]+rainFlux[0]+snowFlux[0]+riverRunoffFlux[0]+iceRunoffFlux[0]+seaIceFreshWaterFlux[0]), m3ps_to_Sv*(layerThick[0]+frazilThick[0]))
boh
volNetLateralFlux = m3ps_to_Sv * volNetLateralFlux
evapFlux = 1/rho0 * m3ps_to_Sv * evapFlux
rainFlux = 1/rho0 * m3ps_to_Sv * rainFlux
snowFlux = 1/rho0 * m3ps_to_Sv * snowFlux
riverRunoffFlux = 1/rho0 * m3ps_to_Sv * riverRunoffFlux
iceRunoffFlux = 1/rho0 * m3ps_to_Sv * iceRunoffFlux
seaIceFreshWaterFlux = 1/rho0 * m3ps_to_Sv * seaIceFreshWaterFlux
thickTend = m3ps_to_Sv * (layerThick + frazilThick)
res = thickTend - (volNetLateralFlux + evapFlux + rainFlux + snowFlux + riverRunoffFlux + iceRunoffFlux + seaIceFreshWaterFlux)
figdpi = 300
figsize = (16, 16)
volNetLateralFlux_runavg = pd.Series.rolling(pd.DataFrame(volNetLateralFlux), movingAverageMonths, center=True).mean()
evapFlux_runavg = pd.Series.rolling(pd.DataFrame(evapFlux), movingAverageMonths, center=True).mean()
rainFlux_runavg = pd.Series.rolling(pd.DataFrame(rainFlux), movingAverageMonths, center=True).mean()
snowFlux_runavg = pd.Series.rolling(pd.DataFrame(snowFlux), movingAverageMonths, center=True).mean()
riverRunoffFlux_runavg = pd.Series.rolling(pd.DataFrame(riverRunoffFlux), movingAverageMonths, center=True).mean()
iceRunoffFlux_runavg = pd.Series.rolling(pd.DataFrame(iceRunoffFlux), movingAverageMonths, center=True).mean()
seaIceFreshWaterFlux_runavg = pd.Series.rolling(pd.DataFrame(seaIceFreshWaterFlux), movingAverageMonths, center=True).mean()
thickTend_runavg = pd.Series.rolling(pd.DataFrame(thickTend), movingAverageMonths, center=True).mean()
res_runavg = pd.Series.rolling(pd.DataFrame(res), movingAverageMonths, center=True).mean()
figfile = f'{figdir}/volBudget_pointIcell{iCell:d}_{casename}_years{year1:04d}-{year2:04d}.png'
fig, ax = plt.subplots(5, 2, figsize=figsize)
ax[0, 0].plot(t, volNetLateralFlux, 'k', alpha=0.5, linewidth=1.5)
ax[0, 1].plot(t, evapFlux, 'k', alpha=0.5, linewidth=1.5)
ax[1, 0].plot(t, rainFlux, 'k', alpha=0.5, linewidth=1.5)
ax[1, 1].plot(t, snowFlux, 'k', alpha=0.5, linewidth=1.5)
ax[2, 0].plot(t, riverRunoffFlux, 'k', alpha=0.5, linewidth=1.5)
ax[2, 1].plot(t, iceRunoffFlux, 'k', alpha=0.5, linewidth=1.5)
ax[3, 0].plot(t, seaIceFreshWaterFlux, 'k', alpha=0.5, linewidth=1.5)
ax[3, 1].plot(t, thickTend, 'k', alpha=0.5, linewidth=1.5)
ax[4, 0].plot(t, res, 'k', alpha=0.5, linewidth=1.5)
if movingAverageMonths==1:
ax[4, 1].plot(t, volNetLateralFlux, 'r', linewidth=2, label='netLateral')
ax[4, 1].plot(t, evapFlux+rainFlux+snowFlux, 'c', linewidth=2, label='E-P')
ax[4, 1].plot(t, riverRunoffFlux+iceRunoffFlux, 'g', linewidth=2, label='runoff')
ax[4, 1].plot(t, seaIceFreshWaterFlux, 'b', linewidth=2, label='seaiceFW')
ax[4, 1].plot(t, -thickTend, 'm', linewidth=2, label='thickTend')
ax[4, 1].plot(t, res, 'k', linewidth=2, label='res')
else:
ax[0, 0].plot(t, volNetLateralFlux_runavg, 'k', linewidth=3)
ax[0, 1].plot(t, evapFlux_runavg, 'k', linewidth=3)
ax[1, 0].plot(t, rainFlux_runavg, 'k', linewidth=3)
ax[1, 1].plot(t, snowFlux_runavg, 'k', linewidth=3)
ax[2, 0].plot(t, riverRunoffFlux_runavg, 'k', linewidth=3)
ax[2, 1].plot(t, iceRunoffFlux_runavg, 'k', linewidth=3)
ax[3, 0].plot(t, seaIceFreshWaterFlux_runavg, 'k', linewidth=3)
ax[3, 1].plot(t, thickTend_runavg, 'k', linewidth=3)
ax[4, 0].plot(t, res_runavg, 'k', linewidth=3)
ax[4, 1].plot(t, volNetLateralFlux_runavg, 'r', linewidth=2, label='netLateral')
ax[4, 1].plot(t, evapFlux_runavg+rainFlux_runavg+snowFlux_runavg, 'c', linewidth=2, label='E-P')
ax[4, 1].plot(t, riverRunoffFlux_runavg+iceRunoffFlux_runavg, 'g', linewidth=2, label='runoff')
ax[4, 1].plot(t, seaIceFreshWaterFlux_runavg, 'b', linewidth=2, label='seaiceFW')
ax[4, 1].plot(t, -thickTend_runavg, 'm', linewidth=2, label='thickTend')
ax[4, 1].plot(t, res_runavg, 'k', linewidth=2, label='res')
ax[4, 1].set_title(f'{movingAverageMonths}-month running averages', fontsize=16, fontweight='bold')
ax[4, 1].legend(loc='lower left')
#ax[0, 0].plot(t, np.zeros_like(t), 'k', linewidth=0.5)
#ax[0, 1].plot(t, np.zeros_like(t), 'k', linewidth=0.5)
#ax[1, 0].plot(t, np.zeros_like(t), 'k', linewidth=0.5)
#ax[1, 1].plot(t, np.zeros_like(t), 'k', linewidth=0.5)
#ax[2, 0].plot(t, np.zeros_like(t), 'k', linewidth=0.5)
#ax[2, 1].plot(t, np.zeros_like(t), 'k', linewidth=0.5)
#ax[3, 0].plot(t, np.zeros_like(t), 'k', linewidth=0.5)
#ax[3, 1].plot(t, np.zeros_like(t), 'k', linewidth=0.5)
#ax[4, 0].plot(t, np.zeros_like(t), 'k', linewidth=0.5)
#ax[4, 1].plot(t, np.zeros_like(t), 'k', linewidth=0.5)
ax[0, 0].autoscale(enable=True, axis='x', tight=True)
ax[0, 1].autoscale(enable=True, axis='x', tight=True)
ax[1, 0].autoscale(enable=True, axis='x', tight=True)
ax[1, 1].autoscale(enable=True, axis='x', tight=True)
ax[2, 0].autoscale(enable=True, axis='x', tight=True)
ax[2, 1].autoscale(enable=True, axis='x', tight=True)
ax[3, 0].autoscale(enable=True, axis='x', tight=True)
ax[3, 1].autoscale(enable=True, axis='x', tight=True)
ax[4, 0].autoscale(enable=True, axis='x', tight=True)
ax[4, 1].autoscale(enable=True, axis='x', tight=True)
plot_xtick_format('gregorian', np.min(t), np.max(t), maxXTicks=20)
ax[0, 0].grid(color='k', linestyle=':', linewidth = 0.5)
ax[0, 1].grid(color='k', linestyle=':', linewidth = 0.5)
ax[1, 0].grid(color='k', linestyle=':', linewidth = 0.5)
ax[1, 1].grid(color='k', linestyle=':', linewidth = 0.5)
ax[2, 0].grid(color='k', linestyle=':', linewidth = 0.5)
ax[2, 1].grid(color='k', linestyle=':', linewidth = 0.5)
ax[3, 0].grid(color='k', linestyle=':', linewidth = 0.5)
ax[3, 1].grid(color='k', linestyle=':', linewidth = 0.5)
ax[4, 0].grid(color='k', linestyle=':', linewidth = 0.5)
ax[4, 1].grid(color='k', linestyle=':', linewidth = 0.5)
ax[0, 0].set_title(f'mean={np.nanmean(volNetLateralFlux):.2e} $\pm$ {np.nanstd(volNetLateralFlux):.2e}', \
fontsize=16, fontweight='bold')
ax[0, 1].set_title(f'mean={np.nanmean(evapFlux):.2e} $\pm$ {np.nanstd(evapFlux):.2e}', \
fontsize=16, fontweight='bold')
ax[1, 0].set_title(f'mean={np.nanmean(rainFlux):.2e} $\pm$ {np.nanstd(rainFlux):.2e}', \
fontsize=16, fontweight='bold')
ax[1, 1].set_title(f'mean={np.nanmean(snowFlux):.2e} $\pm$ {np.nanstd(snowFlux):.2e}', \
fontsize=16, fontweight='bold')
ax[2, 0].set_title(f'mean={np.nanmean(riverRunoffFlux):.2e} $\pm$ {np.nanstd(riverRunoffFlux):.2e}', \
fontsize=16, fontweight='bold')
ax[2, 1].set_title(f'mean={np.nanmean(iceRunoffFlux):.2e} $\pm$ {np.nanstd(iceRunoffFlux):.2e}', \
fontsize=16, fontweight='bold')
ax[3, 0].set_title(f'mean={np.nanmean(seaIceFreshWaterFlux):.2e} $\pm$ {np.nanstd(seaIceFreshWaterFlux):.2e}', \
fontsize=16, fontweight='bold')
ax[3, 1].set_title(f'mean={np.nanmean(thickTend):.2e} $\pm$ {np.nanstd(thickTend):.2e}', \
fontsize=16, fontweight='bold')
ax[4, 0].set_title(f'mean={np.nanmean(res):.2e} $\pm$ {np.nanstd(res):.2e}', \
fontsize=16, fontweight='bold')
ax[4, 0].set_xlabel('Time (Days)', fontsize=12, fontweight='bold')
ax[4, 1].set_xlabel('Time (Years)', fontsize=12, fontweight='bold')
ax[0, 0].set_ylabel('Net lateral flux (Sv)', fontsize=12, fontweight='bold')
ax[0, 1].set_ylabel('Evap flux (Sv)', fontsize=12, fontweight='bold')
ax[1, 0].set_ylabel('Rain flux (Sv)', fontsize=12, fontweight='bold')
ax[1, 1].set_ylabel('Snow flux (Sv)', fontsize=12, fontweight='bold')
ax[2, 0].set_ylabel('River runoff flux (Sv)', fontsize=12, fontweight='bold')
ax[2, 1].set_ylabel('Ice runoff flux (Sv)', fontsize=12, fontweight='bold')
ax[3, 0].set_ylabel('Sea ice FW flux (Sv)', fontsize=12, fontweight='bold')
ax[3, 1].set_ylabel('Layer+frazil thickness tend (Sv)', fontsize=12, fontweight='bold')
ax[4, 0].set_ylabel('Sum of all terms (Sv)', fontsize=12, fontweight='bold')
ax[4, 1].set_ylabel('Sv', fontsize=12, fontweight='bold')
fig.tight_layout(pad=0.5)
fig.savefig(figfile, dpi=figdpi, bbox_inches='tight')