Changeset 24

Timestamp:

11/28/07 16:13:23 (16 years ago)

Author:

kolasinski

Message:

Add make_energetics in the macros directory

File:

• trunk/procs/macros/make_energetics.pro (modified) (17 diffs)

trunk/procs/macros/make_energetics.pro

@@ -6 +6 @@
 @common
 @com_eg
-
+IF debug_w THEN print, 'cmd_wrk.grid =',  cmd_wrk.grid
 CASE cmd_wrk.grid OF
    'T': source_model = 'opa'
+   'T#': source_model = 'opa#'
    ELSE: source_model = 'ipcc'
 ENDCASE 
@@ -22 +23 @@
 ; Read T, S, U, V, W, taux, tauy
 ;
+   print, '    Data size to read (M)= ', nxt*nyt*jpt*(5*nzt + 2)*1.e-6
+
 CASE source_model OF
    'opa': BEGIN 
@@ -43 +46 @@
       file_temp = file_name
    END
+   'opa#': BEGIN 
+      tn = nc_read(file_name,'votemper', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
+      IF data_domain EQ 'pacific' THEN BEGIN 
+         file_nams = strmid(file_name, 0, strlen(file_name)-17)+'T_pac_vosaline.nc'
+         file_namu = strmid(file_name, 0, strlen(file_name)-17)+'U_pac_vozocrtx.nc'
+         file_namv = strmid(file_name, 0, strlen(file_name)-17)+'V_pac_vomecrty.nc'
+         file_namw = strmid(file_name, 0, strlen(file_name)-17)+'W_pac_vovecrtz.nc'
+      ENDIF ELSE BEGIN 
+         file_nams = strmid(file_name, 0, strlen(file_name)-13)+'T_vosaline.nc'
+         file_namu = strmid(file_name, 0, strlen(file_name)-13)+'U_vozocrtx.nc'
+         file_namv = strmid(file_name, 0, strlen(file_name)-13)+'V_vomecrty.nc'
+         file_namw = strmid(file_name, 0, strlen(file_name)-13)+'W_vovecrtz.nc'
+      ENDELSE      
+      sn = nc_read(file_nams,'vosaline', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
+      un = nc_read(file_namu,'vozocrtx', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
+      vn = nc_read(file_namv,'vomecrty', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
+      wn = nc_read(file_namw,'vovecrtz', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
+      tauxn = nc_read(file_namu,'sozotaux', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
+      tauyn = nc_read(file_namv,'sometauy', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
+      var_temp = 'votemper'
+      file_temp = file_name
+      source_model = 'opa'
+   END
    'ipcc': BEGIN
       base_file_name_grd = base_file_name+base_suffix
@@ -57 +83 @@
 ENDCASE 
 
-
-   t = tn.data
-   s = sn.data
-   u = un.data
-   v = vn.data
-   w = wn.data
-   tx = tauxn.data
-   ty = tauyn.data
+t_unit = tn.units
+
+; assign and free memory
+
+   t = tn.data & tn = 0
+   s = sn.data & sn = 0
+   u = un.data & un = 0
+   v = vn.data & vn = 0
+   w = wn.data & wn = 0
+   tx = tauxn.data & tauxn = 0
+   ty = tauyn.data & tauyn = 0
 
    rg = 9.81
+
+   print, '     read ok'
 
 ; rearrange data depending on source
@@ -72 +103 @@
 CASE source_model OF
    'opa': BEGIN  
+      ; build mask
+      IF valmask EQ 0. THEN BEGIN 
+         tmaskr = abs(s) gt 1.e-6
+         valmask = 1.e20
+         w(where (tmaskr EQ 0)) = valmask
+         idxt = where (tmaskr EQ 0)
+         t (where (tmaskr EQ 0)) = 0.
+         s (where (tmaskr EQ 0)) = 0.
+         tmaskr = 0 ; free memory
+      ENDIF ELSE BEGIN 
+         idxt = where (t EQ valmask)
+         IF idxt(0) NE -1 THEN t (idxt) = 0.
+         IF idxt(0) NE -1 THEN s (idxt) = 0.
+      ENDELSE 
+
       ; transform W fields onto T grid
       maskw = w LT valmask/10.
       w_T = 0.5*( w*maskw + shift(w, 0, 0, -1, 0)*shift(maskw, 0, 0, -1, 0) )
       w_T(*, *, (size(w))(3)-1, *) = w_T(*, *, (size(w))(3)-2, *)
+      w = 0 
+      maskw = 0  ; free memory
+      
    END 
    'ipcc': BEGIN
-      w_T = w
+      idx_t = where (t GT valmask/10) 
+      ; special case CCCMA
+      IF strpos(cmd_wrk.exp, 'CCCMA') NE -1 THEN BEGIN 
+         ; transform W fields onto T grid
+         maskw = w LT valmask/10.
+         w_T = 0.5*( w*maskw + shift(w, 0, 0, -1, 0)*shift(maskw, 0, 0, -1, 0) )
+         w_T(*, *, (size(w))(3)-1, *) = w_T(*, *, (size(w))(3)-2, *)
+         w = 0 
+         w_T(idx_t) = valmask
+         maskw = 0 ; free memory
+      ENDIF ELSE w_T = w
+      
       idx_2d = where (u(*, *, 0, 0) GT valmask/10.)
-      tx(idx_2d) = valmask
+      IF idx_2d(0) NE -1 THEN tx(idx_2d) = valmask
       idx_2d = where (v(*, *, 0, 0) GT valmask/10.)
-      ty(idx_2d) = valmask
+      IF idx_2d(0) NE -1 THEN ty(idx_2d) = valmask
       idx = where (t LT valmask/10.)
-      t(idx) = t(idx)-273.15
+      IF t_unit NE "C" THEN BEGIN 
+         IF idx(0) NE -1 THEN t(idx) = t(idx)-273.15
+      ENDIF 
+      idx1 = n_elements(idx_t)
+      idx_w = where (w_T GT valmask/10)
+      idx2 = n_elements(idx_w)
+      IF idx1 NE idx2 THEN BEGIN
+         print, ' *** WARNING wo and thetao do not have same masks !',  idx1,  idx2, min(idx_t -idx_w), max(idx_t- idx_w)
+         print, ' *** Setting wo to zero on all masked points'
+         w_T(idx_w) = 0.
+      ENDIF ELSE print, '       check these are 0 (W,T): ', min(idx_t -idx_w), max(idx_t- idx_w)
+      w = 0 & idx = 0 & idx_t = 0 & idx_w = 0 ; free memory
+      idxt=where(t GT valmask/10.)
+      idxs=where(s GT valmask/10.)
+      IF source_model EQ 'ipcc' THEN print, '       check these are 0 (T,S): ', min(idxt -idxs), max(idxt-idxs)
+      IF idxt[0] NE -1 THEN t(idxt)=0.
+      IF idxs[0] NE -1 THEN s(idxs)=0.
+      idxs = 0                  ; free memory
+
    END 
 ENDCASE 
@@ -91 +169 @@
 ; compute potential density rho
 
-   idxt=where(t GT valmask/10.)
-   idxs=where(s GT valmask/10.)
-   IF idxt[0] NE -1 THEN t(idxt)=0.
-   IF idxs[0] NE -1 THEN s(idxs)=0.
+   print, '     compute rho...'
+
    
    sr=sqrt(abs(s))
@@ -102 +178 @@
    r3=(-1.6546E-6*t+1.0227E-4)*t-5.72466E-3
    rhop = ( ( 4.8314E-4*s + r3*sr +r2)*s +r1)  
+   print, '      rhop min/max', min(rhop), max(rhop)
    IF idxt[0] NE -1 THEN rhop(idxt) = valmask
+   sr = 0 & r1 = 0 &  r2 = 0 &  r3 = 0 & t = 0 & s = 0 ; free memory
 
 ; compute mean profiles on T grid
 
    vargrid = 'T'
-   rho_s = grossemoyenne(rhop, 'xyt', boite = zbox, NaN = valmask)
+   IF valmask NE 0 THEN BEGIN 
+      rho_s = grossemoyenne(rhop, 'xyt', boite = zbox, NaN = valmask)
+   ENDIF ELSE rho_s = grossemoyenne(rhop, 'xyt', boite = zbox)
 
    rho_s4d = replicate(1, nxt*nyt)#rho_s
@@ -116 +196 @@
 
    rho_diff = (rho_s-shift(rho_s,-1))/shift(e3w, -1)
-   rho_diff = shift(rho_diff, 1)
+   rho_diff = shift(temporary(rho_diff), 1)
    rho_diff(0) = 0.
 
@@ -126 +206 @@
    stab_inv = ABS(1./rho_diff_T)
 
+   rho_diff_T = 0 ; free memory
+
 ; remove first 2 levels (MXL too unstable)
 
@@ -142 +224 @@
    IF idxt[0] NE -1 THEN int_val2(idxt) = 0.
 
+   print, '     compute APE...'
+
    ape = 0.5*rg*grossemoyenne(int_val2, 'xyz', /integration)
+
+   int_val2 = 0 ; free memory
 
    ape_wr = ape
@@ -149 +235 @@
 ; compute buoyancy forcing bfx = int[(rho-rho_s).w]dxdydz
 ;
+   print, '     compute BF...'
+
    int_val = (rhop-rho_s4d)*(w_T)
 ; remove first 2 levels (MXL too unstable)
@@ -156 +244 @@
    bfx = rg*grossemoyenne(int_val, 'xyz', /integration)
 
+   int_val = 0 ; free memory
+
    bfx_wr = bfx
    bfx_b = bfx
@@ -161 +251 @@
 
 ; compute wind work = int(tau.um)dx.dy where um=u(over 30 m)
+
+   print, '     compute WW...'
 
    umean=grossemoyenne(u,'z',boite=[0,30])
@@ -170 +262 @@
    idyv = where(vmean GT valmask/10.)
 
-   tx(idx) = 0.
-   ty(idy) = 0.
-   umean(idxu) = 0.
-   vmean(idyv) = 0.
+   IF idx[0] NE -1 THEN tx(idx) = 0.
+   IF idy[0] NE -1 THEN ty(idy) = 0.
+   IF idxu[0] NE -1 THEN umean(idxu) = 0.
+   IF idyv[0] NE -1 THEN vmean(idyv) = 0.
 
    dot_prodx = tx*umean
@@ -213 +305 @@
    IF ps EQ 1 THEN openps
 
-   pltt, ape, 't',  petit = [2, 4, 1],  landscape = 1, /rempli, /BASICMARGES,  title = 'APE (full)'
-   pltt, ww, 't',  petit = [2, 4, 2],  min = -1, max = 5, /noerase, /rempli, /BASICMARGES,  title = 'Wind work (full)'
-   pltt, bfx, 't',  petit = [2, 4, 8],  min = -1, max = 5, color = 4, /noerase, /rempli, /BASICMARGES,  title = 'B (full)'
+   pltt, ape, 't',  petit = [2, 4, 1],  landscape = 1, /rempli, /BASICMARGES,  title = 'APE (full) '+cmd_wrk.exp
+   pltt, ww, 't',  petit = [2, 4, 2],  min = -2, max = 5, /noerase, /rempli, /BASICMARGES,  title = 'Wind work (full) '+cmd_wrk.exp
+   pltt, bfx, 't',  petit = [2, 4, 8],  min = -2, max = 5, color = 4, /noerase, /rempli, /BASICMARGES,  title = 'B (full)'
 
    ape_1mm = trends(ape, 412, 't')
@@ -243 +335 @@
 ; compute and plot sst in nino 3   
 
-   domdef, [210., 270., -5., 5., 0., 10.]
+   maxdepth = 10
+   IF gdept(0) GT maxdepth THEN maxdepth = gdept(0)
+
+   domdef, [210., 280., -5., 5., 0., maxdepth]
    tn = nc_read(file_temp,var_temp, ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
    st = tn.data
@@ -299 +394 @@
    print, ' B/W efficiency (SC) = ', efficiency_sc
 
-   stop
+;   stop
 
    field = {name: '', data: rhop, legend: '', units: '', origin: '', direc: '', dim: 0}