-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathprofiles_T_S_rho.py
executable file
·104 lines (84 loc) · 3.45 KB
/
profiles_T_S_rho.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
#!/usr/bin/env python
'''
Script to compare some scalar values from different runs of Thwaites melt variability experiment.
'''
import sys
import os
import netCDF4
import datetime
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal
from matplotlib import cm
pii=3.14159
fmesh=netCDF4.Dataset('/project/projectdirs/acme/inputdata/ocn/mpas-o/oEC60to30v3wLI/oEC60to30v3wLI60lev.171031.nc')
latCell = fmesh.variables['latCell'][:]
lonCell = fmesh.variables['lonCell'][:]
xCell = fmesh.variables['xCell'][:]
yCell = fmesh.variables['yCell'][:]
depths = fmesh.variables['refBottomDepth'][:]
#idx=67250-1
#idx = np.argmin( (latCell-(-1.323514))**2 + (lonCell-5.672896)**2) #122901-1
idx=198673-1
#idx = np.argmin( (latCell-(-77.75/180.0*pii))**2 + (lonCell-((360.0-36.15)/180.0*pii ))**2) # Darelius 2016 Msouth
#idx = np.argmin( (latCell-(-77.0083/180.0*pii))**2 + (lonCell-((360.0-34.05)/180.0*pii ))**2) # Darelius 2016 Mnorth
#idx=210384-1 # filchner sill
print "idx=",idx
maxLevelCell=fmesh.variables['maxLevelCell'][idx]
nz = maxLevelCell
bottomDepth = fmesh.variables['bottomDepth'][idx]
layerThickness = fmesh.variables['layerThickness'][0, idx, :nz]
thicknessSum = layerThickness.sum()
thicknessCumSum = layerThickness.cumsum()
zSurf = bottomDepth - thicknessSum
zLayerBot = zSurf - thicknessCumSum
z = zLayerBot + 0.5 * layerThickness
path='/global/cscratch1/sd/dcomeau/acme_scratch/cori-knl/20190225.GMPAS-DIB-IAF.T62_oEC60to30v3wLI.cori-knl/archive/ocn/hist'
#path='/global/cscratch1/sd/hoffman2/acme_scratch/edison/20190423.GMPAS-DIB-IAF-ISMF.T62_oEC60to30v3wLI.edison.restrictedMelt/run'
#path='/global/cscratch1/sd/dcomeau/acme_scratch/cori-knl/20190225.GMPAS-DIB-IAF-ISMF.T62_oEC60to30v3wLI.cori-knl/run'
#path='/global/cscratch1/sd/hoffman2/acme_scratch/edison/archive/20190423.GMPAS-DIB-IAF-ISMF.T62_oEC60to30v3wLI.edison.restrictedMelt/ocn/hist'
#years = np.arange(70,105,3)
years = np.arange(5,85,10)
mo=1
colors = [ cm.jet(x) for x in np.linspace(0.0, 1.0, len(years))]
fig = plt.figure(1, facecolor='w')
nrow=1
ncol=4
axT = fig.add_subplot(nrow, ncol, 1)
plt.xlabel('temperature (deg. C)')
plt.ylabel('depth (m)')
plt.grid()
axS = fig.add_subplot(nrow, ncol, 2, sharey=axT)
plt.xlabel('salinity (psu)')
plt.ylabel('depth (m)')
plt.grid()
axrho = fig.add_subplot(nrow, ncol, 3, sharey=axT)
plt.xlabel('potential density (kg/m^3)')
plt.ylabel('depth (m)')
plt.grid()
axvel = fig.add_subplot(nrow, ncol, 4, sharey=axT)
plt.xlabel('velo (m/s)')
plt.ylabel('depth (m)')
plt.grid()
for i in range(len(years)):
yr = years[i]
c = colors[i]
print "yr=",yr
f=netCDF4.Dataset(path+'/'+'mpaso.hist.am.timeSeriesStatsMonthly.{0:04d}-{1:02d}-01.nc'.format(yr, mo),'r')
T=f.variables['timeMonthly_avg_activeTracers_temperature']
S=f.variables['timeMonthly_avg_activeTracers_salinity']
rho = f.variables['timeMonthly_avg_potentialDensity']
v=f.variables['timeMonthly_avg_velocityMeridional']
u=f.variables['timeMonthly_avg_velocityZonal']
axT.plot(T[0,idx,:nz], z, label="yr{0:04d}".format(yr), color=c)
axT.plot(-0.0575*S[0,idx,:nz]+0.0901-7.61e-4*z, z, 'k:')
axS.plot(S[0,idx,:nz], z, label="yr{0:04d}".format(yr), color=c)
axrho.plot(rho[0,idx,:nz], z, label="yr{0:04d}".format(yr), color=c)
axvel.plot((v[0,idx,:nz]**2+u[0,idx,:nz]**2)**0.5, z, label="yr{0:04d}".format(yr), color=c)
#axT.set_xlim([-10,20])
#axS.set_xlim([28,36])
#axrho.set_xlim([1000,1040])
#axT.legend()
axS.legend()
#axT.plot([-1.8, -1.8], [-4000.0, 0.0], 'k:')
plt.show()