Attachment 'eady2.py'
Download 1 import sys; sys.path.append('../py_src')
2 from pyOM_gui import pyOM_gui as pyOM
3 from numpy import *
4 import lin_stab
5
6 #H0 = 200 # water depth
7 #N0 = 0.001 # stability freq.
8 #EPS = 1e-6
9 #U0 = 0.15 # velocity mag.
10
11
12 H0 = 2000 # water depth
13 N0 = 0.004 # stability freq.
14 U0 = 0.5 # velocity mag.
15 EPS = 1e-5
16
17 f0 = 1e-4 # coriolis freq.
18 beta= 0e-11 # beta
19 HY = 0.0e-3 # slope of topography
20
21 KX_INITAL = 1.0 # initial with wave with wavenumber k = KX_INITAL * kmax
22 # where kmax is wavenumber of fastest growing mode
23
24 class eady2(pyOM):
25 """ Eady (1941) solution
26 """
27 def set_parameter(self):
28 """set main parameter
29 """
30 M=self.fortran.main_module
31 (M.nx,M.ny,M.nz) = (34,34,10)
32
33 M.congr_epsilon = 1e-12
34 M.congr_max_iterations = 5000
35 M.enable_streamfunction = 1
36
37 M.enable_superbee_advection = 1
38 M.enable_explicit_vert_friction = 1
39 M.enable_hor_friction = 1
40 M.kappam_0 = 1.e-4
41 M.kappah_0 = 1.e-4
42
43 M.enable_hydrostatic = 1
44 M.enable_cyclic_x = 1
45 M.enable_conserve_energy = 0
46 M.coord_degree = 0
47 M.eq_of_state_type = 1
48 M.enable_tempsalt_sources = 1
49 return
50
51 def set_grid(self):
52 M=self.fortran.main_module
53 # print some numbers
54 Lr = N0*H0/f0 # Rossby radius
55 delta = H0/Lr # aspect ratio
56 Ro = U0/(f0*Lr) # Rossby number
57 Ri = N0**2*H0**2/U0**2 # Richardson number
58 print
59 print ' L = %f km'%(Lr/1e3)
60 print ' Ro = %f '%Ro
61 print ' Ri = %f '%Ri
62 print ' delta = %f '%delta
63 print ' ell = %f '%(Lr/6400e3)
64 print
65 # solve linear stability problem first
66 ksx=linspace(0,3.2,40); kx=ksx/Lr
67 ky = array( [0./Lr] )
68 M.dzt[:] = H0/M.nz
69 print ' Delta z = ',M.dzt[0]
70 zw=arange(M.nz)*M.dzt[0]+ M.dzt[0]
71 zt=zw-M.dzt[0]/2.0
72 U=U0/2+U0*zt/H0
73 V=U*0
74 B=N0**2*zt
75 if 1:
76 om_max,om,kmax,lmax,u,v,w,b,p=lin_stab.qg(U,V,B,M.dzt[0],kx,ky,beta,f0,0,HY)
77 print ' Max. growth rate QG %f 1/days ' % (-imag(om)*86400)
78 print ' k_max Lr = %f , l_max L_r = %f ' % (kmax*Lr/pi,lmax*Lr/pi)
79 if 0:
80 om_max,om,kmax,lmax,u,v,w,b,p=lin_stab.pe(U,V,B,M.dzt[0],kx,ky,0.,beta,0.,f0,0.,HY)
81 print ' Max. growth rate PE %f 1/days ' % (-imag(om)*86400)
82 print ' k_max = %f Lr , l_max = %f Lr' % (kmax*Lr,lmax*Lr)
83 print
84 self.lin_stab_kmax = kmax
85 self.lin_stab_b = b
86 self.lin_stab_p = p
87 self.lin_stab_u = u
88 self.lin_stab_v = v
89 self.lin_stab_w = w
90 L = 2*2*pi/kmax # two waves
91 M.dxt[:] = L/M.nx
92 M.dyt[:] = L/M.nx
93 M.dt_tracer = M.dxt[0]*1200/20e3 # c = 20e3/1200.0, dt =dx/c
94 M.dt_mom = M.dxt[0]*1200/20e3 # c = 20e3/1200.0, dt =dx/c
95 print "dx=%f km, dt= %f s "%(M.dxt[0]/1e3,M.dt_tracer)
96 #M.a_h = (M.dxt[0])**3*2e-11
97 return
98
99
100
101 def set_coriolis(self):
102 M=self.fortran.main_module
103 for j in range( M.yt.shape[0] ): M.coriolis_t[:,j] = f0+beta*M.yt[j]
104 return
105
106 def set_forcing(self):
107 M=self.fortran.main_module
108
109 # update density, etc of last time step
110 M.temp[:,:,:,M.tau-1] = M.temp[:,:,:,M.tau-1] + self.t0
111 self.fortran.calc_eq_of_state(M.tau)
112 M.temp[:,:,:,M.tau-1] = M.temp[:,:,:,M.tau-1] - self.t0
113
114 # advection of background temperature
115 self.fortran.advect_tracer(M.is_pe-M.onx,M.ie_pe+M.onx,M.js_pe-M.onx,M.je_pe+M.onx,self.t0,self.dt0[...,M.tau-1],M.nz)
116 M.temp_source[:] = (1.5+M.ab_eps)*self.dt0[...,M.tau-1] - ( 0.5+M.ab_eps)*self.dt0[...,M.taum1-1]
117 return
118
119
120 def set_initial_conditions(self):
121 """ setup all initial conditions
122 """
123 M=self.fortran.main_module
124 alpha = self.fortran.linear_eq_of_state.linear_eq_of_state_drhodt()
125 grav = 9.81; rho0 = 1024.0
126
127 # background velocity
128 for k in range(M.nz):
129 M.u[:,:,k,M.tau-1]=(U0/2+U0*M.zt[k]/H0)*M.masku[:,:,k]
130
131 # background buoyancy
132 self.t0 = zeros( (M.i_blk+2*M.onx,M.j_blk+2*M.onx,M.nz) , 'd', order='F')
133 self.dt0 = zeros( (M.i_blk+2*M.onx,M.j_blk+2*M.onx,M.nz,3) , 'd', order='F')
134
135 for k in range(M.nz):
136 self.t0[:,:,k]=-N0**2*M.zt[k]/grav/alpha*rho0*M.maskt[:,:,k]
137
138 for k in range(M.nz):
139 for j in range(1,M.ny+1):
140 jj = self.jf2py(j)
141 f = (M.coriolis_t[0,jj]+M.coriolis_t[0,jj+1])/2.0
142 uz = U0/M.ht[:,jj]
143 self.t0[:,jj+1,k]=(self.t0[:,jj,k]+M.dyt[jj]*uz*f/grav/alpha*rho0)*M.maskt[:,jj,k]
144
145 # perturbation buoyancy, etc
146 for k in range(M.nz):
147 for j in range(M.yt.shape[0]):
148 ky=pi/(M.ny*M.dyt[0])
149 kx = KX_INITAL*self.lin_stab_kmax
150 phase = cos(kx*M.xt) +complex(0,1)*sin(kx*M.xt)
151 phase = phase*sin(ky*M.yt[j])
152 amp = 0.2
153 M.temp[:,j,k,M.tau-1] = amp*real(phase*self.lin_stab_b[k])*M.maskt[:,j,k]*rho0/(grav*alpha)
154 M.u[:,j,k,M.tau-1] = M.u[:,j,k,M.tau-1] + amp*real(phase*self.lin_stab_u[k])*M.masku[:,j,k]
155 M.v[:,j,k,M.tau-1] = amp*real(phase*self.lin_stab_v[k])*M.maskv[:,j,k]
156 M.w[:,j,k,M.tau-1] = amp*real(phase*self.lin_stab_w[k])*M.maskw[:,j,k]
157
158 for d in (M.temp,M.u,M.v,M.w): d[...,M.taum1-1] = d[...,M.tau-1]
159
160 return
161
162
163 def topography(self):
164 """ setup topography
165 """
166 M=self.fortran.main_module
167 #for j in range(M.ny):
168 # z0 = M.yt[j]*HY
169 # for k in range(M.nz):
170 # if z0 > M.zt[k]-M.zt[0] :
171 # M.maskt[:,j,k]=0
172 M.kbot[:]=0
173 M.kbot[:,2:-2]=1
174 return
175
176
177
178 def make_plot(self):
179 """ make a plot using methods of self.figure
180 """
181 if hasattr(self,'figure'):
182 M=self.fortran.main_module # fortran module with model variables
183 k=M.nz*3/4
184 k2=M.nz/4
185 i=int(M.nx/2)
186 j=int(M.ny/2)
187 x=M.xt[2:-2]/1e3
188 y=M.yt[2:-2]/1e3
189 z=M.zt[:]
190
191 self.figure.clf()
192 ax=self.figure.add_subplot(221)
193 a=M.temp[2:-2,2:-2,:,M.tau-1] + self.t0[2:-2,2:-2,:]
194 b=M.maskt[2:-2,2:-2,:]
195 a=sum(a*b,axis=0)/sum(b,axis=0)
196 co=ax.contourf(y,z,a.transpose())
197 a=M.u[2:-2,2:-2,:,M.tau-1]
198 b=M.masku[2:-2,2:-2,:]
199 a=sum(a*b,axis=0)/sum(b,axis=0)
200 ax.contour(y,z,a.transpose(),10,colors='black')
201 ax.set_ylabel('z [km]')
202
203 self.figure.colorbar(co)
204 ax.set_title("total temp. at x=%4.1f km"%x[i]);
205 ax.axis('tight')
206
207 ax=self.figure.add_subplot(222)
208 co=ax.contourf(x,y,M.temp[2:-2,2:-2,k,M.tau-1].transpose())
209 a=M.temp[2:-2,2:-2,k,M.tau-1] + self.t0[2:-2,2:-2,k]
210 ax.contour(x,y,a.transpose(),10,colors='black')
211 ax.set_title("temp. perturb. at z=%4.1f m"%z[k]);
212 self.figure.colorbar(co)
213 ax.axis('tight')
214
215 ax=self.figure.add_subplot(223)
216 co=ax.contourf(x,y,M.temp[2:-2,2:-2,k2,M.tau-1].transpose())
217 a=M.temp[2:-2,2:-2,k2,M.tau-1] + self.t0[2:-2,2:-2,k2]
218 ax.contour(x,y,a.transpose(),10,colors='black')
219 ax.set_title("temp. perturb. at z=%4.1f m"%z[k2]);
220 a=M.u[2:-2:2,2:-2:2,k,M.tau-1]
221 b=M.v[2:-2:2,2:-2:2,k,M.tau-1]
222 ax.quiver(x[::2],y[::2],a.transpose(),b.transpose() )
223 self.figure.colorbar(co)
224 ax.axis('tight')
225
226
227 ax=self.figure.add_subplot(224)
228 a=M.temp[2:-2,j,:,M.tau-1]
229 co=ax.contourf(x,z,a.transpose())
230 ax.set_title("temp. perturb. at y=%4.1f km"%y[j]);
231 ax.set_xlabel('x [km]');
232 self.figure.colorbar(co)
233 ax.axis('tight')
234 return
235
236 if __name__ == "__main__": eady2().mainloop( snapint = 86400.0)
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.