source: trunk/procs/macros/make_transfo.pro @ 170

;+
;
; compute transformation
;
; @param todo {in}
;
; input denflx = density flux
;
; @uses
; <pro>common</pro>
; <propost_it>com_eg</propost_it>
;
; @history
; EG january 2001 (inspired from denflux.f)
;
; @version
; $Id$
;
;-
FUNCTION make_transfo, file_name, ncdf_db, todo $
         , BOXZOOM=boxzoom $
         , TIME_1=time_1 $
         , TIME_2=time_2 $
         , ALL_DATA=all_data $
         , ZMTYP=zmtyp
;
  compile_opt idl2, strictarrsubs
;
@common
@com_eg
;
;
; Read sst (C), sss (PSU), Q [W/m2], E-P [kg/m2/s]=[mm/s]
;
   sst = nc_read(file_name,'sosstsst', ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1, TIME_2 = time_2, ALL_DATA = all_data)
   sss = nc_read(file_name,'sosaline', ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1, TIME_2 = time_2, ALL_DATA = all_data)
   hef = nc_read(file_name,'sohefldo', ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1, TIME_2 = time_2, ALL_DATA = all_data)
   emp = nc_read(file_name,'sowaflup', ncdf_db, BOXZOOM = boxzoom, TIME_1 = time_1, TIME_2 = time_2, ALL_DATA = all_data)
   sst = sst.data
   sss = sss.data
   hef = hef.data
   emp = emp.data

   idxt=where(sst eq valmask)
   idxs=where(sss eq valmask)
   IF idxt[0] NE -1 THEN sst[idxt]=0.
   IF idxs[0] NE -1 THEN sss[idxs]=0.   

; Compute potential density (sigma_0)
                                                     
   rho0 = eos[sst, sss]-1000.
   IF idxs[0] NE -1 THEN rho0[idxt] = valmask   

; Define density grid

   n_sig = (sig_max - sig_min)/sig_del + 1
   z_sig = sig_min+findgen(n_sig)*sig_del

; other inits

   sizef = size(sst)
   jpi = sizef[1]
   jpj = sizef[2]
   jpt = 1
   IF sizef[0] EQ 3 THEN jpt = sizef[3]

   dt = 1./float(jpt)            ; integration delta time
   drho = sig_del          ; delta rho
   dsig = sig_del
   P = 0.                        ; pressure = 0 (potential density)

   t_area = 0.
   t_heat = 0.
   t_wafl = 0.
   fluxbar = 0.
   fluxsal = 0.

; Output arrays

   fheat = fltarr(jpi, jpj, jpt)
   fwafl = fltarr(jpi, jpj, jpt)
   dia = fltarr(n_sig+1)
   diawf = fltarr(n_sig+1)
   diat = fltarr( n_sig+1, jpt)
   diawft = fltarr(n_sig+1, jpt)
   hist = fltarr(n_sig+1, jpt)
   histwf = fltarr(n_sig+1, jpt)
   area = fltarr(n_sig+1, jpt)
   nval = lonarr(n_sig+1, jpt)

; Work arrays

   timwrk = fltarr(jpt) 

   trf_bin = fltarr(n_sig+1)
   trf_bin[*] = 0. 
;
;  Loop on 2D domain (and time as '*' dimension)
;
   FOR ji = 0, (jpi-1) DO BEGIN 
      FOR jj = 0, (jpj-1) DO BEGIN 
 
         IF rho0[ji, jj, 0] NE valmask THEN BEGIN 
;
;        input fields
            t = sst[ji, jj, *]
            s = sss[ji, jj, *]
            hf = hef[ji, jj, *]
            ep = emp[ji, jj, *]
            sig = rho0[ji, jj, *] 

;        compte alpha, beta, Cp
            al = alpha[t, s, P]
            be = betar[t, s, P]
            Cp = Cpsw[t, s, P]

;        find density bin
            ib = ROUND((sig-sig_min)/dsig)
            
            ib = ib > 0
            ib = ib < n_sig

;        compute area and integral of fluxes

            dA = e1t[ji, jj]*e2t[ji, jj]
        ; conv = m2 -> 1.e6 km2
            conv = 1.e-6*1.e-6
            t_area = t_area + dA*conv
        ; conv = W -> PW
            conv = 1.e-15
            t_heat = t_heat + total(hf)*dt*dA*conv  
        ; conv = mm -> m and m3/s to Sv
            conv = 1.e-3*1.e-6
            t_wafl = t_wafl + total(ep)*dt*dA*conv  

;        tally number of point in each bin 

            FOR jt = 0, (jpt-1) DO BEGIN 
               area[ib[jt], jt] = area[ib[jt], jt] + dA
               nval[ib[jt], jt] = nval[ib[jt], jt] + 1
            ENDFOR 

;        sum heat and fresh water fluxes as mass fluxes in kg/m2/s (all units SI)
;        convwf : kg/m2/s=mm/s -> m/s
            convwf = 1.e-3
            fheat[ji, jj, *] = ( -al/Cp)*hf
            fwafl[ji, jj, *] = (1000.+sig)*be*s*ep*convwf

;        volume transport per density class (m3/s)
;        period averaged density flux
            FOR jt = 0, (jpt-1) DO BEGIN 
               fht = fheat[ji, jj, jt]*dA
               fwf = fwafl[ji, jj, jt]*dA
               hist[ib[jt], jt] = hist[ib[jt], jt] + (fht+fwf)/drho
               histwf[ib[jt], jt] = histwf[ib[jt], jt] + fwf/drho
               fluxbar = fluxbar + (fht+fwf)*dt
               fluxsal = fluxsal + (fwf)*dt
            ENDFOR 

         ENDIF 
;     end of loop on ji,jj
      ENDFOR 
   ENDFOR 


   print, ' '
   print, ' ----------------------------------------------------------------'
   print, '                    Transformation statitics'
   print, ' ----------------------------------------------------------------'
   print, '   Net heat flux in ocean               (PW)  : ', t_heat
   print, '   Net water flux in ocean              (Sv)  : ', t_wafl
   print, '   Period averaged density flux       (kg/s)  : ', fluxbar
   print, '      - heat flux contribution        (kg/s)  : ', fluxbar-fluxsal
   print, '      - fresh water flux contribution (kg/s)  : ', fluxsal
   print, '   Integral area (glob.=357.208)  (10^6 km2)  : ', t_area
   print, ' '

; output

   CASE todo OF ;
      'denflx': result = (fheat+fwafl)           ; density flux
      'denflw': result = (fwafl)                 ; fresh water density flux
      ELSE: result = -1
   ENDCASE 
                                          
   return, result
END