forked from knaughten/roms_tools
-
Notifications
You must be signed in to change notification settings - Fork 0
/
adv_sss.py
138 lines (121 loc) · 5.37 KB
/
adv_sss.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
from netCDF4 import *
from numpy import *
from matplotlib.pyplot import *
#import colormaps as cmaps
# Create a 3x2 plot of sea surface salinity as seen by CICE on 23 August (the
# sea ice area max) for each advection experiment: the absolute values for
# U3_LIM, and the anomalies from U3_LIM for the other 5 experiments.
def adv_sss ():
# Paths to simulation directories
paths = ['/short/m68/kaa561/advection/u3_lim/', '/short/m68/kaa561/advection/u3/', '/short/m68/kaa561/advection/c4_l/', '/short/m68/kaa561/advection/c4_h/', '/short/m68/kaa561/advection/a4_l/', '/short/m68/kaa561/advection/a4_h/']
# Titles for plotting
labels = ['a) U3_LIM', 'b) U3 - U3_LIM', 'c) C4_LD - U3_LIM', 'd) C4_HD - U3_LIM', 'e) A4_LD - U3_LIM', 'f) A4_HD - U3_LIM']
# File name: daily average for 23 August
file_tail = 'cice/rundir/history/iceh.1992-08-23.nc'
# Bounds and ticks for colour scales
min_abs = 33
max_abs = 35
tick_abs = 0.5
max_anom = 1.0 #2.0
tick_anom = 0.5 #1.0
# Degrees to radians conversion factor
deg2rad = pi/180.
# Centre of missing circle in grid
lon_c = 50
lat_c = -83
# Radius of missing circle
radius = 10.5
# Boundary of regular grid to embed circle in
circle_bdry = -70+90
lon_ticks = array([-120, -60, 60, 120, 180])
lat_ticks = array([-44, -42, -42, -44, -41])
lon_labels = [r'120$^{\circ}$W', r'60$^{\circ}$W', r'60$^{\circ}$E', r'120$^{\circ}$E', r'180$^{\circ}$']
lon_rot = [-60, 60, -60, 60, 0]
# Read salinity data from U3_LIM simulation; also grid and mask variables
id = Dataset(paths[0]+file_tail, 'r')
data_tmp = id.variables['sss'][0,:350,:]
lon_tmp = id.variables['TLON'][:350,:]
lat_tmp = id.variables['TLAT'][:350,:]
mask_tmp = id.variables['tmask'][:350,:]
id.close()
# Wrap periodic boundary so there isn't a gap in the plot
lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1])
lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1])
mask = ma.empty([size(mask_tmp,0), size(mask_tmp,1)+1])
data0 = ma.empty([size(data_tmp,0), size(data_tmp,1)+1])
lon[:,:-1] = lon_tmp
lon[:,-1] = lon_tmp[:,0]
lat[:,:-1] = lat_tmp
lat[:,-1] = lat_tmp[:,0]
mask[:,:-1] = mask_tmp
mask[:,-1] = mask_tmp[:,0]
data0[:,:-1] = data_tmp
data0[:,-1] = data_tmp[:,0]
# Land mask
land = ma.masked_where(mask==1, mask)
# Circumpolar x and y coordinates for plotting
x = -(lat+90)*cos(lon*deg2rad+pi/2)
y = (lat+90)*sin(lon*deg2rad+pi/2)
# Coordinates of centre of missing circle
x_c = -(lat_c+90)*cos(lon_c*deg2rad+pi/2)
y_c = (lat_c+90)*sin(lon_c*deg2rad+pi/2)
# Longitude labels
x_ticks = -(lat_ticks+90)*cos(lon_ticks*deg2rad+pi/2)
y_ticks = (lat_ticks+90)*sin(lon_ticks*deg2rad+pi/2)
# Regular grid to embed missing circle in
x_reg, y_reg = meshgrid(linspace(-circle_bdry, circle_bdry, num=100), linspace(-circle_bdry, circle_bdry, num=100))
# Mask everything except the circle out of the regular grid
land_circle = zeros(shape(x_reg))
land_circle = ma.masked_where(sqrt((x_reg-x_c)**2 + (y_reg-y_c)**2) > radius, land_circle)
# Set up figure
fig = figure(figsize=(9,15))
ax = fig.add_subplot(3, 2, 1, aspect='equal')
# First shade land
contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6')))
# Fill in missing circle
contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6')))
# Shade the salinity data (pcolor not contourf so we don't misrepresent the
# model grid)
img0 = pcolor(x, y, data0, vmin=min_abs, vmax=max_abs, cmap='jet') #cmaps.viridis)
# Add longitude labels
for i in range(size(x_ticks)):
text(x_ticks[i], y_ticks[i], lon_labels[i], ha='center', rotation=lon_rot[i])
axis('off')
# Add title
title(labels[0], fontsize=20)
# Add colorbar
cbaxes0 = fig.add_axes([0.025, 0.7, 0.02, 0.2])
cbar0 = colorbar(img0, ticks=arange(min_abs, max_abs+tick_abs, tick_abs), cax=cbaxes0, extend='both')
cbar0.ax.tick_params(labelsize=16)
# Loop over the other simulations
for sim in range(1, len(paths)):
# Read the salinity data
id = Dataset(paths[sim]+file_tail, 'r')
data_tmp = id.variables['sss'][0,:350,:]
id.close()
# Wrap the periodic boundary
data = ma.empty([size(data_tmp,0), size(data_tmp,1)+1])
data[:,:-1] = data_tmp
data[:,-1] = data_tmp[:,0]
# Calculate anomaly from U3_LIM
data = data - data0
# Add to plot, same as before
ax = fig.add_subplot(3, 2, sim+1, aspect='equal')
contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6')))
contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6')))
img = pcolor(x, y, data, vmin=-max_anom, vmax=max_anom, cmap='RdBu_r')
axis('off')
title(labels[sim], fontsize=20)
if sim == 3:
# Only add an anomaly colorbar for one of the simulations
cbaxes = fig.add_axes([0.025, 0.4, 0.02, 0.2])
cbar = colorbar(img, ticks=arange(-max_anom, max_anom+tick_anom, tick_anom), cax=cbaxes, extend='both')
cbar.ax.tick_params(labelsize=16)
# Main title
suptitle('Sea surface salinity (psu) on 23 August', fontsize=28)
subplots_adjust(wspace=0.025,hspace=0.15)
#fig.show()
fig.savefig('adv_sss.png')
# Command-line interface
if __name__ == '__main__':
adv_sss()