Changeset 24
- Timestamp:
- 11/28/07 16:13:23 (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/procs/macros/make_energetics.pro
r2 r24 6 6 @common 7 7 @com_eg 8 8 IF debug_w THEN print, 'cmd_wrk.grid =', cmd_wrk.grid 9 9 CASE cmd_wrk.grid OF 10 10 'T': source_model = 'opa' 11 'T#': source_model = 'opa#' 11 12 ELSE: source_model = 'ipcc' 12 13 ENDCASE … … 22 23 ; Read T, S, U, V, W, taux, tauy 23 24 ; 25 print, ' Data size to read (M)= ', nxt*nyt*jpt*(5*nzt + 2)*1.e-6 26 24 27 CASE source_model OF 25 28 'opa': BEGIN … … 43 46 file_temp = file_name 44 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 45 71 'ipcc': BEGIN 46 72 base_file_name_grd = base_file_name+base_suffix … … 57 83 ENDCASE 58 84 59 60 t = tn.data 61 s = sn.data 62 u = un.data 63 v = vn.data 64 w = wn.data 65 tx = tauxn.data 66 ty = tauyn.data 85 t_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 67 96 68 97 rg = 9.81 98 99 print, ' read ok' 69 100 70 101 ; rearrange data depending on source … … 72 103 CASE source_model OF 73 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 74 120 ; transform W fields onto T grid 75 121 maskw = w LT valmask/10. 76 122 w_T = 0.5*( w*maskw + shift(w, 0, 0, -1, 0)*shift(maskw, 0, 0, -1, 0) ) 77 123 w_T(*, *, (size(w))(3)-1, *) = w_T(*, *, (size(w))(3)-2, *) 124 w = 0 125 maskw = 0 ; free memory 126 78 127 END 79 128 'ipcc': BEGIN 80 w_T = w 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 81 141 idx_2d = where (u(*, *, 0, 0) GT valmask/10.) 82 tx(idx_2d) = valmask142 IF idx_2d(0) NE -1 THEN tx(idx_2d) = valmask 83 143 idx_2d = where (v(*, *, 0, 0) GT valmask/10.) 84 ty(idx_2d) = valmask144 IF idx_2d(0) NE -1 THEN ty(idx_2d) = valmask 85 145 idx = where (t LT valmask/10.) 86 t(idx) = t(idx)-273.15 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 87 165 END 88 166 ENDCASE … … 91 169 ; compute potential density rho 92 170 93 idxt=where(t GT valmask/10.) 94 idxs=where(s GT valmask/10.) 95 IF idxt[0] NE -1 THEN t(idxt)=0. 96 IF idxs[0] NE -1 THEN s(idxs)=0. 171 print, ' compute rho...' 172 97 173 98 174 sr=sqrt(abs(s)) … … 102 178 r3=(-1.6546E-6*t+1.0227E-4)*t-5.72466E-3 103 179 rhop = ( ( 4.8314E-4*s + r3*sr +r2)*s +r1) 180 print, ' rhop min/max', min(rhop), max(rhop) 104 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 105 183 106 184 ; compute mean profiles on T grid 107 185 108 186 vargrid = 'T' 109 rho_s = grossemoyenne(rhop, 'xyt', boite = zbox, NaN = valmask) 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) 110 190 111 191 rho_s4d = replicate(1, nxt*nyt)#rho_s … … 116 196 117 197 rho_diff = (rho_s-shift(rho_s,-1))/shift(e3w, -1) 118 rho_diff = shift( rho_diff, 1)198 rho_diff = shift(temporary(rho_diff), 1) 119 199 rho_diff(0) = 0. 120 200 … … 126 206 stab_inv = ABS(1./rho_diff_T) 127 207 208 rho_diff_T = 0 ; free memory 209 128 210 ; remove first 2 levels (MXL too unstable) 129 211 … … 142 224 IF idxt[0] NE -1 THEN int_val2(idxt) = 0. 143 225 226 print, ' compute APE...' 227 144 228 ape = 0.5*rg*grossemoyenne(int_val2, 'xyz', /integration) 229 230 int_val2 = 0 ; free memory 145 231 146 232 ape_wr = ape … … 149 235 ; compute buoyancy forcing bfx = int[(rho-rho_s).w]dxdydz 150 236 ; 237 print, ' compute BF...' 238 151 239 int_val = (rhop-rho_s4d)*(w_T) 152 240 ; remove first 2 levels (MXL too unstable) … … 156 244 bfx = rg*grossemoyenne(int_val, 'xyz', /integration) 157 245 246 int_val = 0 ; free memory 247 158 248 bfx_wr = bfx 159 249 bfx_b = bfx … … 161 251 162 252 ; compute wind work = int(tau.um)dx.dy where um=u(over 30 m) 253 254 print, ' compute WW...' 163 255 164 256 umean=grossemoyenne(u,'z',boite=[0,30]) … … 170 262 idyv = where(vmean GT valmask/10.) 171 263 172 tx(idx) = 0.173 ty(idy) = 0.174 umean(idxu) = 0.175 vmean(idyv) = 0.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. 176 268 177 269 dot_prodx = tx*umean … … 213 305 IF ps EQ 1 THEN openps 214 306 215 pltt, ape, 't', petit = [2, 4, 1], landscape = 1, /rempli, /BASICMARGES, title = 'APE (full) '216 pltt, ww, 't', petit = [2, 4, 2], min = - 1, max = 5, /noerase, /rempli, /BASICMARGES, title = 'Wind work (full)'217 pltt, bfx, 't', petit = [2, 4, 8], min = - 1, max = 5, color = 4, /noerase, /rempli, /BASICMARGES, title = 'B (full)'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)' 218 310 219 311 ape_1mm = trends(ape, 412, 't') … … 243 335 ; compute and plot sst in nino 3 244 336 245 domdef, [210., 270., -5., 5., 0., 10.] 337 maxdepth = 10 338 IF gdept(0) GT maxdepth THEN maxdepth = gdept(0) 339 340 domdef, [210., 280., -5., 5., 0., maxdepth] 246 341 tn = nc_read(file_temp,var_temp, ncdf_db, TIME_1 = time_1, TIME_2 = time_2) 247 342 st = tn.data … … 299 394 print, ' B/W efficiency (SC) = ', efficiency_sc 300 395 301 stop396 ; stop 302 397 303 398 field = {name: '', data: rhop, legend: '', units: '', origin: '', direc: '', dim: 0}
Note: See TracChangeset
for help on using the changeset viewer.