Attachment 'acc2.py'
Download 1 import sys; sys.path.append('../py_src')
2
3 from pyOM_gui import pyOM_gui as pyOM
4 #from pyOM_ave import pyOM_ave as pyOM
5
6 from numpy import *
7
8 yt_start = -39.0
9 yt_end = 43
10 yu_start = -40.0
11 yu_end = 42
12
13 class acc2(pyOM):
14 """ a simple global model with a Southern Ocean and Atlantic part
15 """
16 def set_parameter(self):
17 """set main parameter
18 """
19 M=self.fortran.main_module
20
21 (M.nx,M.ny,M.nz) = (30,42,15)
22 M.dt_mom = 4800
23 M.dt_tracer = 86400/2.0
24
25 M.coord_degree = 1
26 M.enable_cyclic_x = 1
27
28 M.congr_epsilon = 1e-12
29 M.congr_max_iterations = 5000
30 M.enable_streamfunction = 1
31
32 I=self.fortran.isoneutral_module
33 I.enable_neutral_diffusion = 1
34 I.k_iso_0 = 1000.0
35 I.k_iso_steep = 500.0
36 I.iso_dslope=0.005
37 I.iso_slopec=0.01
38 I.enable_skew_diffusion = 1
39
40 M.enable_hor_friction = 1
41 M.a_h = (2*M.degtom)**3*2e-11
42 M.enable_hor_friction_cos_scaling = 1
43 M.hor_friction_cospower=1
44
45 M.enable_bottom_friction = 1
46 M.r_bot = 1e-5
47
48 M.enable_implicit_vert_friction = 1
49 T=self.fortran.tke_module
50 T.enable_tke = 1
51 T.c_k = 0.1
52 T.c_eps = 0.7
53 T.alpha_tke = 30.0
54 T.mxl_min = 1e-8
55 T.tke_mxl_choice = 2
56 #T.enable_tke_superbee_advection = 1
57
58 M.k_gm_0 = 1000.0
59 E=self.fortran.eke_module
60 E.enable_eke = 1
61 E.eke_k_max = 1e4
62 E.eke_c_k = 0.4
63 E.eke_c_eps = 0.5
64 E.eke_cross = 2.
65 E.eke_crhin = 1.0
66 E.eke_lmin = 100.0
67 E.enable_eke_superbee_advection = 1
68 E.enable_eke_isopycnal_diffusion = 1
69
70 I=self.fortran.idemix_module
71 I.enable_idemix = 1
72 I.enable_idemix_hor_diffusion = 1
73 I.enable_eke_diss_surfbot = 1
74 I.eke_diss_surfbot_frac = 0.2
75 I.enable_idemix_superbee_advection = 1
76
77 M.eq_of_state_type = 3
78 return
79
80
81 def set_grid(self):
82 M=self.fortran.main_module
83 ddz = array([50.,70.,100.,140.,190.,240.,290.,340.,390.,440.,490.,540.,590.,640.,690.])
84 M.dxt[:] = 2.0
85 M.dyt[:] = 2.0
86 M.x_origin= 0.0
87 M.y_origin= -40.0
88 M.dzt[:] = ddz[::-1]/2.5
89 return
90
91
92 def set_coriolis(self):
93 M=self.fortran.main_module
94 for j in range( M.yt.shape[0] ): M.coriolis_t[:,j] = 2*M.omega*sin(M.yt[j]/180.*pi)
95 return
96
97 def set_topography(self):
98 """ setup topography
99 """
100 M=self.fortran.main_module
101 (X,Y)= meshgrid(M.xt,M.yt); X=X.transpose(); Y=Y.transpose()
102 M.kbot[:]=0
103 M.kbot[ X >1.0]=1
104 M.kbot[ Y <-20]=1
105 return
106
107 def set_initial_conditions(self):
108 """ setup initial conditions
109 """
110 M=self.fortran.main_module
111 (XT,YT)= meshgrid(M.xt,M.yt); XT=XT.transpose(); YT=YT.transpose()
112 (XU,YU)= meshgrid(M.xu,M.yu); XU=XU.transpose(); YU=YU.transpose()
113
114 # initial conditions
115 for k in range(M.nz):
116 M.temp[:,:,k,M.tau-1] = (1-M.zt[k]/M.zw[0])*15*M.maskt[:,:,k]
117 M.temp[:,:,k,M.taum1-1] = (1-M.zt[k]/M.zw[0])*15*M.maskt[:,:,k]
118 M.salt[:,:,:,M.tau-1] = 35.0*M.maskt[:]
119 M.salt[:,:,:,M.taum1-1] = 35.0*M.maskt[:]
120
121 # wind stress forcing
122 for j in range(M.js_pe,M.je_pe+1):
123 jj = self.jf2py(j)
124 taux=0.0
125 if M.yt[jj]<-20 : taux = .1e-3*sin(pi*(M.yu[jj]-yu_start)/(-20.0-yt_start))
126 if M.yt[jj]>10 : taux = .1e-3*(1-cos(2*pi*(M.yu[jj]-10.0)/(yu_end-10.0)))
127 M.surface_taux[:,jj] = taux*M.masku[:,jj,-1]
128
129 # surface heatflux forcing
130 self.t_rest = zeros( M.u[:,:,1,0].shape, order = 'F' )
131 self.t_star = zeros( M.u[:,:,1,0].shape, order = 'F' )
132
133 for j in range(M.js_pe,M.je_pe+1):
134 jj = self.jf2py(j)
135 t_star=15
136 if M.yt[jj]<-20.0: t_star=15*(M.yt[jj]-yt_start )/(-20.0-yt_start)
137 if M.yt[jj]> 20.0: t_star=15*(1-(M.yt[jj]-20)/(yt_end-20) )
138 self.t_star[:,jj] = t_star
139 self.t_rest[:,jj] = M.dzt[-1]/(30.*86400.)*M.maskt[:,jj,-1]
140
141
142 T=self.fortran.tke_module
143 if T.enable_tke:
144 for j in range(M.js_pe,M.je_pe+1):
145 for i in range(M.is_pe,M.ie_pe+1):
146 (ii,jj) = (self.if2py(i), self.jf2py(j) )
147 T.forc_tke_surface[ii,jj] = sqrt( (0.5*(M.surface_taux[ii,jj]+M.surface_taux[ii-1,jj]))**2 \
148 +(0.5*(M.surface_tauy[ii,jj]+M.surface_tauy[ii,jj-1]))**2 )**(3./2.)
149
150 I=self.fortran.idemix_module
151 if I.enable_idemix:
152 I.forc_iw_bottom[:] = 1.0e-6*M.maskw[:,:,-1]
153 I.forc_iw_surface[:] = 0.1e-6*M.maskw[:,:,-1]
154 return
155
156
157 def set_forcing(self):
158 M=self.fortran.main_module
159 M.forc_temp_surface[:]=self.t_rest*(self.t_star-M.temp[:,:,-1,M.tau-1])
160 return
161
162 def set_diagnostics(self):
163 M=self.fortran.main_module
164 self.register_average(name='temp',long_name='Temperature', units = 'deg C' , grid = 'TTT', var = M.temp)
165 self.register_average(name='salt',long_name='Salinity', units = 'g/kg' , grid = 'TTT', var = M.salt)
166 self.register_average(name='u', long_name='Zonal velocity', units = 'm/s' , grid = 'UTT', var = M.u)
167 self.register_average(name='v', long_name='Meridional velocity', units = 'm/s' , grid = 'TUT', var = M.v)
168 self.register_average(name='w', long_name='Vertical velocity', units = 'm/s' , grid = 'TTU', var = M.w)
169 self.register_average(name='taux',long_name='wind stress', units = 'm^2/s' , grid = 'UT', var = M.surface_taux)
170 self.register_average(name='tauy',long_name='wind stress', units = 'm^2/s' , grid = 'TU', var = M.surface_tauy)
171 self.register_average(name='psi' ,long_name='Streamfunction', units = 'm^3/s' , grid = 'UU', var = M.psi)
172 return
173
174 def user_defined_signal(self):
175 """ this routine must be called by all processors
176 """
177 M=self.fortran.main_module
178 a = zeros( (M.nx,M.ny), 'd', order = 'F')
179 a[M.is_pe-1:M.ie_pe,0] = M.xt[2:-2]
180 self.fortran.pe0_recv_2d(a)
181 self.xt_gl = a[:,0].copy()
182
183 a[0,M.js_pe-1:M.je_pe] = M.yt[2:-2]
184 self.fortran.pe0_recv_2d(a)
185 self.yt_gl = a[0,:].copy()
186
187 self.psi_gl = zeros( (M.nx,M.ny), 'd', order = 'F')
188 self.psi_gl[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe] = where( M.maskz[2:-2,2:-2,-1] >0, M.psi[2:-2,2:-2,M.tau-1] , NaN)
189 self.fortran.pe0_recv_2d(self.psi_gl)
190
191 self.temp_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
192 for k in range(M.nz):
193 a[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe] = where( M.maskt[2:-2,2:-2,k] >0, M.temp[2:-2,2:-2,k,M.tau-1] , NaN)
194 self.fortran.pe0_recv_2d(a)
195 self.temp_gl[:,:,k]=a.copy()
196
197 self.kappa_gl = zeros( (M.nx,M.ny,M.nz), 'd', order = 'F')
198 for k in range(M.nz):
199 a[M.is_pe-1:M.ie_pe,M.js_pe-1:M.je_pe] = where( M.maskw[2:-2,2:-2,k] >0, M.kappah[2:-2,2:-2,k] , NaN)
200 self.fortran.pe0_recv_2d(a)
201 self.kappa_gl[:,:,k]=a.copy()
202
203 return
204
205 def make_plot(self):
206 M=self.fortran.main_module
207
208 self.set_signal('user_defined') # following routine is called by all PEs
209 self.user_defined_signal()
210
211 self.figure.clf()
212 ax=self.figure.add_subplot(221)
213
214 co=ax.contourf(self.yt_gl,M.zt,self.temp_gl[M.nx/2-1,:,:].transpose())
215 self.figure.colorbar(co)
216 ax.set_title('temperature')
217 ax.set_ylabel('z [m]')
218 ax.axis('tight')
219
220 ax=self.figure.add_subplot(223)
221 try:
222 co=ax.contourf(self.yt_gl,M.zw,log10(self.kappa_gl[M.nx/2-1,:,:].transpose()) )
223 except:
224 pass
225 self.figure.colorbar(co)
226 ax.set_title('Diffusivity')
227 ax.set_xlabel('Latitude [deg N]')
228 ax.set_ylabel('z [m]')
229 ax.axis('tight')
230
231 ax=self.figure.add_subplot(122)
232 co=ax.contourf(self.xt_gl,self.yt_gl,self.psi_gl.transpose()*1e-6)
233 self.figure.colorbar(co)
234 ax.set_title('Streamfunction [Sv]')
235 ax.set_xlabel('Longitude [deg E]')
236 ax.axis('tight')
237
238 return
239
240 if __name__ == "__main__": acc2().run(snapint= 86400*10.0, runlen = 365*86400.*200)
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.