;
; make energetics computations
;
FUNCTION make_energetics, file_name, ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2, ALL_DATA = all_data, ZMTYP = zmtyp

@common
@com_eg

CASE cmd_wrk.grid OF
   'T': source_model = 'opa'
   ELSE: source_model = 'ipcc'
ENDCASE 

;; full vertical domain
;; imposes vert_type = '0' in plt_def
vert_switch = 0
IF debug_w THEN BEGIN
   print, 'base_file_name:', base_file_name 
   print, 'file_name:', file_name
ENDIF 
;
; Read T, S, U, V, W, taux, tauy
;
CASE source_model OF
   'opa': BEGIN 
      tn = nc_read(file_name,'votemper', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
      sn = nc_read(file_name,'vosaline', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
      IF data_domain EQ 'pacific' THEN BEGIN 
         file_namu = strmid(file_name, 0, strlen(file_name)-8)+'U_pac.nc'
         file_namv = strmid(file_name, 0, strlen(file_name)-8)+'V_pac.nc'
         file_namw = strmid(file_name, 0, strlen(file_name)-8)+'W_pac.nc'
      ENDIF ELSE BEGIN 
         file_namu = strmid(file_name, 0, strlen(file_name)-4)+'U.nc'
         file_namv = strmid(file_name, 0, strlen(file_name)-4)+'V.nc'
         file_namw = strmid(file_name, 0, strlen(file_name)-4)+'W.nc'
      ENDELSE      
      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
   END
   'ipcc': BEGIN
      base_file_name_grd = base_file_name+base_suffix
      tn = nc_read(base_file_name_grd+'_thetao.nc','thetao', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
      sn = nc_read(base_file_name_grd+'_so.nc','so', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
      un = nc_read(base_file_name_grd+'_uo.nc','uo', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
      vn = nc_read(base_file_name_grd+'_vo.nc','vo', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
      wn = nc_read(base_file_name_grd+'_wo.nc','wo', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
      tauxn = nc_read(base_file_name_grd+'_tauu.nc','tauu', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
      tauyn = nc_read(base_file_name_grd+'_tauv.nc','tauv', ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
      var_temp = 'thetao'
      file_temp = base_file_name_grd+'_thetao.nc'
   END
ENDCASE 


   t = tn.data
   s = sn.data
   u = un.data
   v = vn.data
   w = wn.data
   tx = tauxn.data
   ty = tauyn.data

   rg = 9.81

; rearrange data depending on source

CASE source_model OF
   'opa': 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, *)
   END 
   'ipcc': BEGIN
      w_T = w
      idx_2d = where (u(*, *, 0, 0) GT valmask/10.)
      tx(idx_2d) = valmask
      idx_2d = where (v(*, *, 0, 0) GT valmask/10.)
      ty(idx_2d) = valmask
      idx = where (t LT valmask/10.)
      t(idx) = t(idx)-273.15
   END 
ENDCASE 
      

; 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.
   
   sr=sqrt(abs(s))
   r1=((((6.536332E-9*t-1.120083E-6)*t+1.001685E-4)*t $
        -9.095290E-3)*t+6.793952E-2)*t+999.842594
   r2=(((5.3875E-9*t-8.2467E-7)*t+7.6438E-5)*t-4.0899E-3)*t+8.24493E-1
   r3=(-1.6546E-6*t+1.0227E-4)*t-5.72466E-3
   rhop = ( ( 4.8314E-4*s + r3*sr +r2)*s +r1)  
   IF idxt[0] NE -1 THEN rhop(idxt) = valmask

; compute mean profiles on T grid

   vargrid = 'T'
   rho_s = grossemoyenne(rhop, 'xyt', boite = zbox, NaN = valmask)

   rho_s4d = replicate(1, nxt*nyt)#rho_s
   rho_s4d = reform(rho_s4d[*]#replicate(1, jpt), nxt, nyt, nzt, jpt, /overwrite)


; compute mean stability = d(rho_s)/dz (on W grid)

   rho_diff = (rho_s-shift(rho_s,-1))/shift(e3w, -1)
   rho_diff = shift(rho_diff, 1)
   rho_diff(0) = 0.

; transform onto T grid

   rho_diff_T = 0.5*(rho_diff+shift(rho_diff, -1))
   rho_diff_T((size(rho_diff))(1)-1) = rho_diff((size(rho_diff))(1)-2)

   stab_inv = ABS(1./rho_diff_T)

; remove first 2 levels (MXL too unstable)

   stab_inv[0:1] = 0.

; test: remove only top level

;   stab_inv[0:0] = 0.

; compute [(rho-rho_s)**2]/stability

   stab_inv = replicate(1, nxt*nyt)#stab_inv
   stab_inv = reform(stab_inv[*]#replicate(1, jpt), nxt, nyt, nzt, jpt, /overwrite)

   int_val2 = ((rhop-rho_s4d)^2)*stab_inv
   IF idxt[0] NE -1 THEN int_val2(idxt) = 0.

   ape = 0.5*rg*grossemoyenne(int_val2, 'xyz', /integration)

   ape_wr = ape
   ape = ape*1.e-18
;
; compute buoyancy forcing bfx = int[(rho-rho_s).w]dxdydz
;
   int_val = (rhop-rho_s4d)*(w_T)
; remove first 2 levels (MXL too unstable)
   IF idxt[0] NE -1 THEN int_val(idxt) = 0.
   int_val[*, *, 0:1, *] = 0.

   bfx = rg*grossemoyenne(int_val, 'xyz', /integration)

   bfx_wr = bfx
   bfx_b = bfx
   bfx = bfx*1.e-11

; compute wind work = int(tau.um)dx.dy where um=u(over 30 m)

   umean=grossemoyenne(u,'z',boite=[0,30])
   vmean=grossemoyenne(v,'z',boite=[0,30])

   idx = where(tx GT valmask/10.)
   idy = where(ty GT valmask/10.)
   idxu = where(umean GT valmask/10.)
   idyv = where(vmean GT valmask/10.)

   tx(idx) = 0.
   ty(idy) = 0.
   umean(idxu) = 0.
   vmean(idyv) = 0.

   dot_prodx = tx*umean
   dot_prody = ty*vmean

   wwx = grossemoyenne(dot_prodx, 'xy', /integration)
   wwy = grossemoyenne(dot_prody, 'xy', /integration)
   ww = wwx + wwy

   ww_wr = ww
   ww_b = ww
   ww = ww*1.e-11
   wwx = wwx*1.e-11
   wwy = wwy*1.e-11

; compute forcing efficiency: stddev(B)/stddev(W)

   bfx_1mm = trends(bfx_b, 412, 't')
   bfx_sc = mean_sc
   ww_1mm = trends(ww_b, 412, 't')
   ww_sc = mean_sc
   efficiency = sqrt((moment(bfx_1mm))[1])/sqrt((moment(ww_1mm))[1])
   efficiency_sc = sqrt((moment(bfx_sc[0:11]))[1])/sqrt((moment(ww_sc[0:11]))[1])

; plotting stuff

   ps = 0

   red =   [0, 255,   0,   0, 0, 255]
   green = [0,   0, 255,   0, 0,   0]
   blue =  [0,   0,   0, 255, 0, 255]
   red = [0, red, red, red, red, red, red, red ]
   green = [0, green, green, green, green, green, green, green]
   blue = [0, blue, blue, blue, blue, blue, blue, blue ]
   tvlct, red, green, blue

   IF cmd_wrk.out EQ 'ps' THEN ps = 1

   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)'

   ape_1mm = trends(ape, 412, 't')
   pltt, ape_1mm, 't',  petit = [2, 4, 3],  /noerase, /rempli, /BASICMARGES,  title = 'APE (inter)'
   jpt_b = jpt
   jpt = 24
   pltt, mean_sc[0:23], 't',  petit = [2, 4, 5],  /noerase, /rempli, /BASICMARGES,  title = 'APE (seasonal cycle x 2)'

   jpt = jpt_b

   ww_1mm = trends(ww, 412, 't')
   tmp = mean_sc
   wwx_1mm = trends(wwx, 412, 't')

   pltt, ww_1mm, 't',  petit = [2, 4, 4], color = 2, /noerase, /rempli, /BASICMARGES,  title = 'Interannual W (red) B (blue) [efficiency = '+string(strcompress(efficiency))+']'
   pltt, bfx_1mm*1.e-11, 't',  petit = [2, 4, 4],  /ov1d, color = 4, thick = 2, /noerase, /rempli, /BASICMARGES

   jpt_b = jpt
   jpt = 24
   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)'
   pltt, mean_sc[0:23], 't',  petit = [2, 4, 6], /ov1d, color = 2, thick = 2, /noerase, /rempli, /BASICMARGES
   wwy_1mm = trends(wwy, 412, 't')
   pltt, mean_sc[0:23], 't',  petit = [2, 4, 6], /ov1d, color = 3, thick = 2, /noerase, /rempli, /BASICMARGES
   pltt, bfx_sc[0:23]*1.e-11, 't',  petit = [2, 4, 6], /ov1d, color = 4, thick = 1, /noerase, /rempli, /BASICMARGES
   jpt = jpt_b

; compute and plot sst in nino 3   

   domdef, [210., 270., -5., 5., 0., 10.]
   tn = nc_read(file_temp,var_temp, ncdf_db,  TIME_1 = time_1,  TIME_2 =  time_2)
   st = tn.data
   sst = grossemoyenne(st, 'xyz')
   sst_wr = sst
   sst = trends(sst, 412, 't')

   pltt, sst, 't',  petit = [2, 4, 7],  /noerase, /rempli, /BASICMARGES,  title = 'Nino3 SSTA'

   correlation = C_CORRELATE(ape, sst, [0])

   IF ps EQ 1 THEN BEGIN
      closeps
      printps
   ENDIF

; write to ascii file

   get_lun, nuldat
   filename = cmd_wrk.exp+'_'+cmd_wrk.date1+'_'+cmd_wrk.spec+'_'+cmd_wrk.plt+'_sst.asc'
   openw, nuldat, asciidir+filename
   print,  '     -> writing nino 3 sst data to ',  asciidir+filename & print,  ' '
   printf, nuldat, sst_wr, format = '(f8.3)'
   free_lun, nuldat & close, nuldat

   get_lun, nuldat
   filename = cmd_wrk.exp+'_'+cmd_wrk.date1+'_'+cmd_wrk.spec+'_'+cmd_wrk.plt+'_ape.asc'
   openw, nuldat, asciidir+filename
   print,  '     -> writing ape data to ',  asciidir+filename & print,  ' '
   printf, nuldat, ape_wr, format = '(g10.4)'
   free_lun, nuldat & close, nuldat

   get_lun, nuldat
   filename = cmd_wrk.exp+'_'+cmd_wrk.date1+'_'+cmd_wrk.spec+'_'+cmd_wrk.plt+'_ww.asc'
   openw, nuldat, asciidir+filename
   print,  '     -> writing ww data to ',  asciidir+filename & print,  ' '
   printf, nuldat, ww_wr, format = '(g10.4)'
   free_lun, nuldat & close, nuldat

   get_lun, nuldat
   filename = cmd_wrk.exp+'_'+cmd_wrk.date1+'_'+cmd_wrk.spec+'_'+cmd_wrk.plt+'_bf.asc'
   openw, nuldat, asciidir+filename
   print,  '     -> writing bf data to ',  asciidir+filename & print,  ' '
   printf, nuldat, bfx_wr, format = '(g10.4)'
   free_lun, nuldat & close, nuldat

; check that d(APE)/dt ~ ww
   
   dapedt = (ape-shift(ape, 1))/(86400.*30.)
;   pltt, dapedt-(ww*1.e11),'t',petit=[1,2,1],/rempli,/portrait         
;   pltt, dapedt/(ww*1.e11),'t',petit=[1,2,2],/rempli,/portrait,/noerase
   print, ' d(APE)/dt / wind work correlation', C_CORRELATE(dapedt, ww, [0]) 
   print, ' APE/nino3 sst correlation=', correlation 
   print, ' B/W efficiency (interannual) = ', efficiency
   print, ' B/W efficiency (SC) = ', efficiency_sc

   stop

   field = {name: '', data: rhop, legend: '', units: '', origin: '', direc: '', dim: 0}
   
   field.origin = tn.origin
   field.dim = tn.dim - 1

   field.direc = 'xyzt'

   return, field
END