-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathplot_cross-section_om4.py
More file actions
156 lines (143 loc) · 4.94 KB
/
Copy pathplot_cross-section_om4.py
File metadata and controls
156 lines (143 loc) · 4.94 KB
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
#!/home/cod/bin/python
##################################################################
## ##
## - Objective: ##
## This script plots sections of constant longitude for ##
## temperature and salinity, and isopycnal surfaces. ##
## ##
## - History: Aug. 2018: original C.O. Dufour ##
## ##
##################################################################
# import modules
from netCDF4 import Dataset
from eos import *
from plottools import *
from utiltools import *
import matplotlib.pyplot as plt
from om4_models import *
import numpy as np
from PyRaf import *
import cmocean
##### to be set #####
# var name
yrdeb = 1988
yrend = 2007
# model
#modelist = [OM4p25,OM4p5,OM4p5_noeddy]
modelist = [OM4p5_noeddy]
# paths
path2grid = '/work/Carolina.Dufour/Mosaic'
path2figs = '/work/Carolina.Dufour/Figures/comp_model2model/OM4'
# domain
latmax = -50
region = 'Weddell'
#####################
## domain
if region=='Weddell': lon0 = -42
## time period
yrs = str(yrdeb).zfill(4)+'-'+str(yrend).zfill(4)
freq = yrend - yrdeb + 1
############### Main ###############
for k in range(len(modelist)):
arch = modelist[k]
print arch.__name__
### paths
path2data = arch.path2ocean+'/av/annual_'+str(freq)+'yr'
### grid
filegrid = arch.path2ocean + '/ocean_annual_z.static.nc'
fid = Dataset(filegrid)
lati = fid.variables['yh'][:]
jmax = find_j(lati,latmax)
lat = lati[:jmax+1]
loni = fid.variables['xq'][:]
i0 = find_j(loni,lon0)
lon = loni[i0]
area = fid.variables['areacello'][:jmax+1,i0]
fid.close()
### load time period
filein = path2data + '/ocean_annual_z.'+str(yrdeb)+'-'+str(yrend)+'.ann.nc'
fid = Dataset(filein)
temp = fid.variables['thetao'][...,:jmax+1,i0].squeeze()
salt = fid.variables['so'][...,:jmax+1,i0].squeeze()
dep = fid.variables['z_l'][:].squeeze()
vol = fid.variables['volcello'][...,:jmax+1,i0].squeeze()
fid.close()
# Compute density
sig0 = sigma_n(temp,salt,0)
sig2 = sigma_n(temp,salt,2)
### Load initial time period
filein = path2data + '/ocean_annual_z.1708-1727.ann.nc'
fid = Dataset(filein)
temp0 = fid.variables['thetao'][...,:jmax+1,i0].squeeze()
salt0 = fid.variables['so'][...,:jmax+1,i0].squeeze()
fid.close()
# Compute density
sig00 = sigma_n(temp0,salt0,0)
sig20 = sigma_n(temp0,salt0,2)
### Deal with partial cell
dz = vol/area
kmax = dz.shape[0]
zl = 0.0 * dz[:]
zl[0] = 0.5 * dz[0]
k=1
while k < kmax:
zl[k] = zl[k-1] + 0.5 * (dz[k]+dz[k-1])
k = k+1
zl_nm = np.array(zl)
######## Plot
## temp
lats,deps = np.meshgrid(lat,dep)
#vmin =-1.5 ; vmax = 4
vmin =-2.5 ; vmax = 3
tab = limit_range(temp,vmin,vmax)
fsize=12
plt.figure(figsize=(12,6))
plt.subplot(111,axisbg='gray')
plt.contourf(lats,zl_nm,tab,30,cmap=plt.cm.magma,vmin=vmin,vmax=vmax)
plt.colorbar(orientation='vertical')
## sigma_0 contours
cl0 = plt.contour(lats,zl_nm,sig0,[27.6,27.7,27.8],colors='k')
cl00 = plt.contour(lats,zl_nm,sig00,[27.6,27.7,27.8],colors='w')
## sigma_2 contours
cl2 = plt.contour(lats,zl_nm,sig2,[37.18,37.2],colors='k')
cl20 = plt.contour(lats,zl_nm,sig20,[37.18,37.2],colors='w')
#
plt.gca().invert_yaxis()
#
plt.clabel(cl0,inline=True,fmt='%.2f',manual=True)
plt.clabel(cl00,inline=True,fmt='%.2f',manual=True)
plt.clabel(cl2,inline=True,fmt='%.2f',manual=True)
plt.clabel(cl20,inline=True,fmt='%.2f',manual=True)
#
plt.xlabel('Latitude',size=fsize)
plt.ylabel('Depth (m)',size=fsize)
plt.xticks(size=fsize) ; plt.yticks(size=fsize)
plt.savefig(path2figs+'/'+arch.__name__+'_'+str(yrdeb).zfill(4)+'-'+str(yrend).zfill(4)+'_temp_sig0_section_'+region+'.png')
plt.show()
## salt
vmin =33.8 ; vmax = 34.8
tab = limit_range(salt,vmin,vmax)
fsize=12
plt.figure(figsize=(12,6))
plt.subplot(111,axisbg='gray')
plt.contourf(lats,zl_nm,tab,30,cmap=cmocean.cm.haline,vmin=vmin,vmax=vmax)
plt.colorbar(orientation='vertical')
## sigma_0 contours
cl0 = plt.contour(lats,zl_nm,sig0,[27.6,27.7,27.8],colors='k')
cl00 = plt.contour(lats,zl_nm,sig00,[27.6,27.7,27.8],colors='w')
## sigma_2 contours
cl2 = plt.contour(lats,zl_nm,sig2,[37.18,37.2],colors='k')
cl20 = plt.contour(lats,zl_nm,sig20,[37.18,37.2],colors='w')
#
plt.gca().invert_yaxis()
#
plt.clabel(cl0,inline=True,fmt='%.2f',manual=True)
plt.clabel(cl00,inline=True,fmt='%.2f',manual=True)
plt.clabel(cl2,inline=True,fmt='%.2f',manual=True)
plt.clabel(cl20,inline=True,fmt='%.2f',manual=True)
#
plt.xlabel('Latitude',size=fsize)
plt.ylabel('Depth (m)',size=fsize)
plt.xticks(size=fsize) ; plt.yticks(size=fsize)
plt.savefig(path2figs+'/'+arch.__name__+'_'+str(yrdeb).zfill(4)+'-'+str(yrend).zfill(4)+'_salt_sig0_section_'+region+'.png')
plt.show()