1 | #!/usr/bin/env python |
---|
2 | # -*- coding: utf-8 -*- |
---|
3 | |
---|
4 | """ Plot mechanical energy evolution |
---|
5 | in the Marsaleix et al (2008) test case (IW1) |
---|
6 | """ |
---|
7 | |
---|
8 | from netCDF4 import Dataset |
---|
9 | import matplotlib.pyplot as plt |
---|
10 | import numpy as np |
---|
11 | |
---|
12 | rau0 = 1026. |
---|
13 | grav = 9.81 |
---|
14 | nexp = 5 |
---|
15 | |
---|
16 | def comp_ke(filedatat, filedatau ): |
---|
17 | global rau0 |
---|
18 | global grav |
---|
19 | |
---|
20 | # Open file at T-points: |
---|
21 | ncid = Dataset(filedatat) |
---|
22 | time = ncid.variables['time_counter'][:] # Time [s] |
---|
23 | lon = ncid.variables['nav_lon'][1, :] # Longitude [km] |
---|
24 | ssh = ncid.variables['ssh_inst'][:, 1, :] # Read sea level anomaly [m] |
---|
25 | e3t = ncid.variables['e3t_inst'][:, :, 1, :] # thicknesses [m] |
---|
26 | temp = ncid.variables['thetao_inst'][:, :, 1, :] # temperature [°C] |
---|
27 | e3t = np.where((temp == 0.), 0., e3t) |
---|
28 | ncid.close() |
---|
29 | |
---|
30 | # Open file at U-points: |
---|
31 | ncid = Dataset(filedatau) |
---|
32 | # ubar = ncid.variables['ubar_inst'][:, 1, :] # Read barotropic velocity [m/s] |
---|
33 | e3u = ncid.variables['e3u_inst'][:, :, 1, :] # U thicknesses [m] |
---|
34 | uu = ncid.variables['uo_inst'][:, :, 1, :] # velocity [m/s] |
---|
35 | ubar = np.zeros(np.shape(ssh)) |
---|
36 | hu = np.zeros(np.shape(ssh)) |
---|
37 | e3u = np.where((uu == 0.), 0., e3u) |
---|
38 | hu = np.sum(e3u, axis=(1)) |
---|
39 | hu = np.where((hu == 0.), 1., hu) |
---|
40 | ubar = np.sum(uu*e3u, axis=(1))/hu |
---|
41 | ncid.close() |
---|
42 | |
---|
43 | # Set velocity at T-points: |
---|
44 | (jpt, jpi) = np.shape(ssh) |
---|
45 | work = np.zeros(np.shape(ssh)) |
---|
46 | work[:, 1:jpi - 1] = 0.5 * (ubar[:, 1:jpi - 1] + ubar[:, 0:jpi - 2]) |
---|
47 | |
---|
48 | # mechanical energy: |
---|
49 | dx = (lon[1]-lon[0])*1.e3 # Horizontal resolution in m |
---|
50 | dt = (time[1]-time[0]) |
---|
51 | dl = (jpi-2)*dx/1.e3 |
---|
52 | h = np.zeros((jpt,1)) # Total depth |
---|
53 | h = np.sum(e3t, axis=(1)) |
---|
54 | ke = np.zeros((jpt, 1)) # energy |
---|
55 | ke = np.sum(0.5*rau0*dx*dx* (h*work**2 + grav*ssh**2), axis=(1)) |
---|
56 | ssh0 = lon / (0.5*dl) |
---|
57 | # ssh0 = np.sin(np.pi*lon/dl) |
---|
58 | ssh0[0] = 0. |
---|
59 | ssh0[jpi-1] = 0. |
---|
60 | ke0 = np.sum(0.5*rau0*dx*dx*grav*ssh0**2, axis=(0)) |
---|
61 | ke = ke/ke0 |
---|
62 | # Set time from time increment: |
---|
63 | time = np.arange(dt,(jpt+1)*dt,dt) |
---|
64 | |
---|
65 | return time, ke |
---|
66 | |
---|
67 | fig = plt.figure(figsize=(7, 5)) |
---|
68 | for exp in np.arange(1, nexp+1,1): |
---|
69 | filet='SEICHE_mar_rk3_dt_%i_grid_T.nc' %exp |
---|
70 | fileu='SEICHE_mar_rk3_dt_%i_grid_U.nc' %exp |
---|
71 | print(filet,fileu) |
---|
72 | (time, ke) = comp_ke(filet, fileu) |
---|
73 | plt.plot(time/86400., ke) |
---|
74 | |
---|
75 | plt.ylim((0.6, 1.05)) |
---|
76 | plt.xlim((0., 12.)) |
---|
77 | plt.xlabel('Time [days]') |
---|
78 | plt.ylabel(r'$\frac{E(t)}{E(0)}$',rotation=0) |
---|
79 | plt.grid() |
---|
80 | plt.title('DM19 test case - Barotropic energy vs baroclinic time step - RK3') |
---|
81 | plt.legend(('960s','480s','240s','120s','60s'),loc='lower left') |
---|
82 | plt.savefig('SEICHE_mar_rk3_dt.png') |
---|
83 | plt.show() |
---|