source: trunk/procs/macros/make_energetics.pro @ 24

Last change on this file since 24 was 24, checked in by kolasinski, 17 years ago

Add make_energetics in the macros directory
File size: 14.3 KB
Line 
1;
2; make energetics computations
3;
4FUNCTION make_energetics, file_name, ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data, ZMTYP = zmtyp
5
6@common
7@com_eg
8IF debug_w THEN print, 'cmd_wrk.grid =',  cmd_wrk.grid
9CASE cmd_wrk.grid OF
10   'T': source_model = 'opa'
11   'T#': source_model = 'opa#'
12   ELSE: source_model = 'ipcc'
13ENDCASE
14
15;; full vertical domain
16;; imposes vert_type = '0' in plt_def
17vert_switch = 0
18IF debug_w THEN BEGIN
19   print, 'base_file_name:', base_file_name
20   print, 'file_name:', file_name
21ENDIF
22;
23; Read T, S, U, V, W, taux, tauy
24;
25   print, '    Data size to read (M)= ', nxt*nyt*jpt*(5*nzt + 2)*1.e-6
26
27CASE source_model OF
28   'opa': BEGIN
29      tn = nc_read(file_name,'votemper', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
30      sn = nc_read(file_name,'vosaline', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
31      IF data_domain EQ 'pacific' THEN BEGIN
32         file_namu = strmid(file_name, 0, strlen(file_name)-8)+'U_pac.nc'
33         file_namv = strmid(file_name, 0, strlen(file_name)-8)+'V_pac.nc'
34         file_namw = strmid(file_name, 0, strlen(file_name)-8)+'W_pac.nc'
35      ENDIF ELSE BEGIN
36         file_namu = strmid(file_name, 0, strlen(file_name)-4)+'U.nc'
37         file_namv = strmid(file_name, 0, strlen(file_name)-4)+'V.nc'
38         file_namw = strmid(file_name, 0, strlen(file_name)-4)+'W.nc'
39      ENDELSE     
40      un = nc_read(file_namu,'vozocrtx', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
41      vn = nc_read(file_namv,'vomecrty', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
42      wn = nc_read(file_namw,'vovecrtz', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
43      tauxn = nc_read(file_namu,'sozotaux', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
44      tauyn = nc_read(file_namv,'sometauy', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
45      var_temp = 'votemper'
46      file_temp = file_name
47   END
48   'opa#': BEGIN
49      tn = nc_read(file_name,'votemper', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
50      IF data_domain EQ 'pacific' THEN BEGIN
51         file_nams = strmid(file_name, 0, strlen(file_name)-17)+'T_pac_vosaline.nc'
52         file_namu = strmid(file_name, 0, strlen(file_name)-17)+'U_pac_vozocrtx.nc'
53         file_namv = strmid(file_name, 0, strlen(file_name)-17)+'V_pac_vomecrty.nc'
54         file_namw = strmid(file_name, 0, strlen(file_name)-17)+'W_pac_vovecrtz.nc'
55      ENDIF ELSE BEGIN
56         file_nams = strmid(file_name, 0, strlen(file_name)-13)+'T_vosaline.nc'
57         file_namu = strmid(file_name, 0, strlen(file_name)-13)+'U_vozocrtx.nc'
58         file_namv = strmid(file_name, 0, strlen(file_name)-13)+'V_vomecrty.nc'
59         file_namw = strmid(file_name, 0, strlen(file_name)-13)+'W_vovecrtz.nc'
60      ENDELSE     
61      sn = nc_read(file_nams,'vosaline', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
62      un = nc_read(file_namu,'vozocrtx', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
63      vn = nc_read(file_namv,'vomecrty', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
64      wn = nc_read(file_namw,'vovecrtz', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
65      tauxn = nc_read(file_namu,'sozotaux', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
66      tauyn = nc_read(file_namv,'sometauy', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
67      var_temp = 'votemper'
68      file_temp = file_name
69      source_model = 'opa'
70   END
71   'ipcc': BEGIN
72      base_file_name_grd = base_file_name+base_suffix
73      tn = nc_read(base_file_name_grd+'_thetao.nc','thetao', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
74      sn = nc_read(base_file_name_grd+'_so.nc','so', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
75      un = nc_read(base_file_name_grd+'_uo.nc','uo', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
76      vn = nc_read(base_file_name_grd+'_vo.nc','vo', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
77      wn = nc_read(base_file_name_grd+'_wo.nc','wo', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
78      tauxn = nc_read(base_file_name_grd+'_tauu.nc','tauu', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
79      tauyn = nc_read(base_file_name_grd+'_tauv.nc','tauv', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
80      var_temp = 'thetao'
81      file_temp = base_file_name_grd+'_thetao.nc'
82   END
83ENDCASE
84
85t_unit = tn.units
86
87; assign and free memory
88
89   t = tn.data & tn = 0
90   s = sn.data & sn = 0
91   u = un.data & un = 0
92   v = vn.data & vn = 0
93   w = wn.data & wn = 0
94   tx = tauxn.data & tauxn = 0
95   ty = tauyn.data & tauyn = 0
96
97   rg = 9.81
98
99   print, '     read ok'
100
101; rearrange data depending on source
102
103CASE source_model OF
104   'opa': BEGIN 
105      ; build mask
106      IF valmask EQ 0. THEN BEGIN
107         tmaskr = abs(s) gt 1.e-6
108         valmask = 1.e20
109         w(where (tmaskr EQ 0)) = valmask
110         idxt = where (tmaskr EQ 0)
111         t (where (tmaskr EQ 0)) = 0.
112         s (where (tmaskr EQ 0)) = 0.
113         tmaskr = 0 ; free memory
114      ENDIF ELSE BEGIN
115         idxt = where (t EQ valmask)
116         IF idxt(0) NE -1 THEN t (idxt) = 0.
117         IF idxt(0) NE -1 THEN s (idxt) = 0.
118      ENDELSE
119
120      ; transform W fields onto T grid
121      maskw = w LT valmask/10.
122      w_T = 0.5*( w*maskw + shift(w, 0, 0, -1, 0)*shift(maskw, 0, 0, -1, 0) )
123      w_T(*, *, (size(w))(3)-1, *) = w_T(*, *, (size(w))(3)-2, *)
124      w = 0
125      maskw = 0  ; free memory
126     
127   END
128   'ipcc': BEGIN
129      idx_t = where (t GT valmask/10)
130      ; special case CCCMA
131      IF strpos(cmd_wrk.exp, 'CCCMA') NE -1 THEN BEGIN
132         ; transform W fields onto T grid
133         maskw = w LT valmask/10.
134         w_T = 0.5*( w*maskw + shift(w, 0, 0, -1, 0)*shift(maskw, 0, 0, -1, 0) )
135         w_T(*, *, (size(w))(3)-1, *) = w_T(*, *, (size(w))(3)-2, *)
136         w = 0
137         w_T(idx_t) = valmask
138         maskw = 0 ; free memory
139      ENDIF ELSE w_T = w
140     
141      idx_2d = where (u(*, *, 0, 0) GT valmask/10.)
142      IF idx_2d(0) NE -1 THEN tx(idx_2d) = valmask
143      idx_2d = where (v(*, *, 0, 0) GT valmask/10.)
144      IF idx_2d(0) NE -1 THEN ty(idx_2d) = valmask
145      idx = where (t LT valmask/10.)
146      IF t_unit NE "C" THEN BEGIN
147         IF idx(0) NE -1 THEN t(idx) = t(idx)-273.15
148      ENDIF
149      idx1 = n_elements(idx_t)
150      idx_w = where (w_T GT valmask/10)
151      idx2 = n_elements(idx_w)
152      IF idx1 NE idx2 THEN BEGIN
153         print, ' *** WARNING wo and thetao do not have same masks !',  idx1,  idx2, min(idx_t -idx_w), max(idx_t- idx_w)
154         print, ' *** Setting wo to zero on all masked points'
155         w_T(idx_w) = 0.
156      ENDIF ELSE print, '       check these are 0 (W,T): ', min(idx_t -idx_w), max(idx_t- idx_w)
157      w = 0 & idx = 0 & idx_t = 0 & idx_w = 0 ; free memory
158      idxt=where(t GT valmask/10.)
159      idxs=where(s GT valmask/10.)
160      IF source_model EQ 'ipcc' THEN print, '       check these are 0 (T,S): ', min(idxt -idxs), max(idxt-idxs)
161      IF idxt[0] NE -1 THEN t(idxt)=0.
162      IF idxs[0] NE -1 THEN s(idxs)=0.
163      idxs = 0                  ; free memory
164
165   END
166ENDCASE
167     
168
169; compute potential density rho
170
171   print, '     compute rho...'
172
173   
174   sr=sqrt(abs(s))
175   r1=((((6.536332E-9*t-1.120083E-6)*t+1.001685E-4)*t $
176        -9.095290E-3)*t+6.793952E-2)*t+999.842594
177   r2=(((5.3875E-9*t-8.2467E-7)*t+7.6438E-5)*t-4.0899E-3)*t+8.24493E-1
178   r3=(-1.6546E-6*t+1.0227E-4)*t-5.72466E-3
179   rhop = ( ( 4.8314E-4*s + r3*sr +r2)*s +r1) 
180   print, '      rhop min/max', min(rhop), max(rhop)
181   IF idxt[0] NE -1 THEN rhop(idxt) = valmask
182   sr = 0 & r1 = 0 &  r2 = 0 &  r3 = 0 & t = 0 & s = 0 ; free memory
183
184; compute mean profiles on T grid
185
186   vargrid = 'T'
187   IF valmask NE 0 THEN BEGIN
188      rho_s = grossemoyenne(rhop, 'xyt', boite = zbox, NaN = valmask)
189   ENDIF ELSE rho_s = grossemoyenne(rhop, 'xyt', boite = zbox)
190
191   rho_s4d = replicate(1, nxt*nyt)#rho_s
192   rho_s4d = reform(rho_s4d[*]#replicate(1, jpt), nxt, nyt, nzt, jpt, /overwrite)
193
194
195; compute mean stability = d(rho_s)/dz (on W grid)
196
197   rho_diff = (rho_s-shift(rho_s,-1))/shift(e3w, -1)
198   rho_diff = shift(temporary(rho_diff), 1)
199   rho_diff(0) = 0.
200
201; transform onto T grid
202
203   rho_diff_T = 0.5*(rho_diff+shift(rho_diff, -1))
204   rho_diff_T((size(rho_diff))(1)-1) = rho_diff((size(rho_diff))(1)-2)
205
206   stab_inv = ABS(1./rho_diff_T)
207
208   rho_diff_T = 0 ; free memory
209
210; remove first 2 levels (MXL too unstable)
211
212   stab_inv[0:1] = 0.
213
214; test: remove only top level
215
216;   stab_inv[0:0] = 0.
217
218; compute [(rho-rho_s)**2]/stability
219
220   stab_inv = replicate(1, nxt*nyt)#stab_inv
221   stab_inv = reform(stab_inv[*]#replicate(1, jpt), nxt, nyt, nzt, jpt, /overwrite)
222
223   int_val2 = ((rhop-rho_s4d)^2)*stab_inv
224   IF idxt[0] NE -1 THEN int_val2(idxt) = 0.
225
226   print, '     compute APE...'
227
228   ape = 0.5*rg*grossemoyenne(int_val2, 'xyz', /integration)
229
230   int_val2 = 0 ; free memory
231
232   ape_wr = ape
233   ape = ape*1.e-18
234;
235; compute buoyancy forcing bfx = int[(rho-rho_s).w]dxdydz
236;
237   print, '     compute BF...'
238
239   int_val = (rhop-rho_s4d)*(w_T)
240; remove first 2 levels (MXL too unstable)
241   IF idxt[0] NE -1 THEN int_val(idxt) = 0.
242   int_val[*, *, 0:1, *] = 0.
243
244   bfx = rg*grossemoyenne(int_val, 'xyz', /integration)
245
246   int_val = 0 ; free memory
247
248   bfx_wr = bfx
249   bfx_b = bfx
250   bfx = bfx*1.e-11
251
252; compute wind work = int(tau.um)dx.dy where um=u(over 30 m)
253
254   print, '     compute WW...'
255
256   umean=grossemoyenne(u,'z',boite=[0,30])
257   vmean=grossemoyenne(v,'z',boite=[0,30])
258
259   idx = where(tx GT valmask/10.)
260   idy = where(ty GT valmask/10.)
261   idxu = where(umean GT valmask/10.)
262   idyv = where(vmean GT valmask/10.)
263
264   IF idx[0] NE -1 THEN tx(idx) = 0.
265   IF idy[0] NE -1 THEN ty(idy) = 0.
266   IF idxu[0] NE -1 THEN umean(idxu) = 0.
267   IF idyv[0] NE -1 THEN vmean(idyv) = 0.
268
269   dot_prodx = tx*umean
270   dot_prody = ty*vmean
271
272   wwx = grossemoyenne(dot_prodx, 'xy', /integration)
273   wwy = grossemoyenne(dot_prody, 'xy', /integration)
274   ww = wwx + wwy
275
276   ww_wr = ww
277   ww_b = ww
278   ww = ww*1.e-11
279   wwx = wwx*1.e-11
280   wwy = wwy*1.e-11
281
282; compute forcing efficiency: stddev(B)/stddev(W)
283
284   bfx_1mm = trends(bfx_b, 412, 't')
285   bfx_sc = mean_sc
286   ww_1mm = trends(ww_b, 412, 't')
287   ww_sc = mean_sc
288   efficiency = sqrt((moment(bfx_1mm))[1])/sqrt((moment(ww_1mm))[1])
289   efficiency_sc = sqrt((moment(bfx_sc[0:11]))[1])/sqrt((moment(ww_sc[0:11]))[1])
290
291; plotting stuff
292
293   ps = 0
294
295   red =   [0, 255,   0,   0, 0, 255]
296   green = [0,   0, 255,   0, 0,   0]
297   blue =  [0,   0,   0, 255, 0, 255]
298   red = [0, red, red, red, red, red, red, red ]
299   green = [0, green, green, green, green, green, green, green]
300   blue = [0, blue, blue, blue, blue, blue, blue, blue ]
301   tvlct, red, green, blue
302
303   IF cmd_wrk.out EQ 'ps' THEN ps = 1
304
305   IF ps EQ 1 THEN openps
306
307   pltt, ape, 't',  petit = [2, 4, 1],  landscape = 1, /rempli, /BASICMARGES,  title = 'APE (full) '+cmd_wrk.exp
308   pltt, ww, 't',  petit = [2, 4, 2],  min = -2, max = 5, /noerase, /rempli, /BASICMARGES,  title = 'Wind work (full) '+cmd_wrk.exp
309   pltt, bfx, 't',  petit = [2, 4, 8],  min = -2, max = 5, color = 4, /noerase, /rempli, /BASICMARGES,  title = 'B (full)'
310
311   ape_1mm = trends(ape, 412, 't')
312   pltt, ape_1mm, 't',  petit = [2, 4, 3],  /noerase, /rempli, /BASICMARGES,  title = 'APE (inter)'
313   jpt_b = jpt
314   jpt = 24
315   pltt, mean_sc[0:23], 't',  petit = [2, 4, 5],  /noerase, /rempli, /BASICMARGES,  title = 'APE (seasonal cycle x 2)'
316
317   jpt = jpt_b
318
319   ww_1mm = trends(ww, 412, 't')
320   tmp = mean_sc
321   wwx_1mm = trends(wwx, 412, 't')
322
323   pltt, ww_1mm, 't',  petit = [2, 4, 4], color = 2, /noerase, /rempli, /BASICMARGES,  title = 'Interannual W (red) B (blue) [efficiency = '+string(strcompress(efficiency))+']'
324   pltt, bfx_1mm*1.e-11, 't',  petit = [2, 4, 4],  /ov1d, color = 4, thick = 2, /noerase, /rempli, /BASICMARGES
325
326   jpt_b = jpt
327   jpt = 24
328   pltt, tmp[0:23], 't',  petit = [2, 4, 6],  min = -1,  max = 3.5,  /noerase, /rempli, /BASICMARGES,  title = 'Seasonal Cycle x 2 (W total: black, B: blue, Wx/y: red/green)'
329   pltt, mean_sc[0:23], 't',  petit = [2, 4, 6], /ov1d, color = 2, thick = 2, /noerase, /rempli, /BASICMARGES
330   wwy_1mm = trends(wwy, 412, 't')
331   pltt, mean_sc[0:23], 't',  petit = [2, 4, 6], /ov1d, color = 3, thick = 2, /noerase, /rempli, /BASICMARGES
332   pltt, bfx_sc[0:23]*1.e-11, 't',  petit = [2, 4, 6], /ov1d, color = 4, thick = 1, /noerase, /rempli, /BASICMARGES
333   jpt = jpt_b
334
335; compute and plot sst in nino 3   
336
337   maxdepth = 10
338   IF gdept(0) GT maxdepth THEN maxdepth = gdept(0)
339
340   domdef, [210., 280., -5., 5., 0., maxdepth]
341   tn = nc_read(file_temp,var_temp, ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
342   st = tn.data
343   sst = grossemoyenne(st, 'xyz')
344   sst_wr = sst
345   sst = trends(sst, 412, 't')
346
347   pltt, sst, 't',  petit = [2, 4, 7],  /noerase, /rempli, /BASICMARGES,  title = 'Nino3 SSTA'
348
349   correlation = C_CORRELATE(ape, sst, [0])
350
351   IF ps EQ 1 THEN BEGIN
352      closeps
353      printps
354   ENDIF
355
356; write to ascii file
357
358   get_lun, nuldat
359   filename = cmd_wrk.exp+'_'+cmd_wrk.date1+'_'+cmd_wrk.spec+'_'+cmd_wrk.plt+'_sst.asc'
360   openw, nuldat, asciidir+filename
361   print,  '     -> writing nino 3 sst data to ',  asciidir+filename & print,  ' '
362   printf, nuldat, sst_wr, format = '(f8.3)'
363   free_lun, nuldat & close, nuldat
364
365   get_lun, nuldat
366   filename = cmd_wrk.exp+'_'+cmd_wrk.date1+'_'+cmd_wrk.spec+'_'+cmd_wrk.plt+'_ape.asc'
367   openw, nuldat, asciidir+filename
368   print,  '     -> writing ape data to ',  asciidir+filename & print,  ' '
369   printf, nuldat, ape_wr, format = '(g10.4)'
370   free_lun, nuldat & close, nuldat
371
372   get_lun, nuldat
373   filename = cmd_wrk.exp+'_'+cmd_wrk.date1+'_'+cmd_wrk.spec+'_'+cmd_wrk.plt+'_ww.asc'
374   openw, nuldat, asciidir+filename
375   print,  '     -> writing ww data to ',  asciidir+filename & print,  ' '
376   printf, nuldat, ww_wr, format = '(g10.4)'
377   free_lun, nuldat & close, nuldat
378
379   get_lun, nuldat
380   filename = cmd_wrk.exp+'_'+cmd_wrk.date1+'_'+cmd_wrk.spec+'_'+cmd_wrk.plt+'_bf.asc'
381   openw, nuldat, asciidir+filename
382   print,  '     -> writing bf data to ',  asciidir+filename & print,  ' '
383   printf, nuldat, bfx_wr, format = '(g10.4)'
384   free_lun, nuldat & close, nuldat
385
386; check that d(APE)/dt ~ ww
387   
388   dapedt = (ape-shift(ape, 1))/(86400.*30.)
389;   pltt, dapedt-(ww*1.e11),'t',petit=[1,2,1],/rempli,/portrait         
390;   pltt, dapedt/(ww*1.e11),'t',petit=[1,2,2],/rempli,/portrait,/noerase
391   print, ' d(APE)/dt / wind work correlation', C_CORRELATE(dapedt, ww, [0])
392   print, ' APE/nino3 sst correlation=', correlation
393   print, ' B/W efficiency (interannual) = ', efficiency
394   print, ' B/W efficiency (SC) = ', efficiency_sc
395
396;   stop
397
398   field = {name: '', data: rhop, legend: '', units: '', origin: '', direc: '', dim: 0}
399   
400   field.origin = tn.origin
401   field.dim = tn.dim - 1
402
403   field.direc = 'xyzt'
404
405   return, field
406END
Note: See TracBrowser for help on using the repository browser.