source: trunk/tools/density_binning/density_bin_IDL_gm/msfs.pro @ 2

function msfs, msfin, sigma_values, mask=mask, orca=orca
@common
msf=msfin
;-------------------------
jpk=(size(msf))[3]
gdept=[min(sigma_values)-0.1,sigma_values]
gdepw=[min(sigma_values)-0.1,sigma_values]+0.05
e3t=replicate(1,jpk)
e3w=replicate(1,jpk)
tmask=replicate(1b,jpi,jpj,jpk)
domdef,0,jpi-1,0,jpj-1,0,jpk-1,/index


if keyword_set(mask) then begin
      print, 'subdomain: use of mask argument'
;       print, mask
     if (size(mask))[0] eq 2 then msf=msf*reform(mask[*]#replicate(1,jpk),jpi,jpj,jpk)
endif 
;
msf=total(msf*(e1v[*]#replicate(1,jpk)),1,/nan)
msf=1e-6*msf

mmsk=msf
indm=where(mmsk NE 0 )
mmsk(indm)=1

msf=total(msf,2,/cumulative)
;msf=total(msf,2,/cumulative,/nan)
msf=shift(msf,0,1)
; if sub basin, substract the bottom value (equivalent to upward integration)
if keyword_set(mask) then begin
     for k=0, jpk-1 do begin
       for j=0,jpj-1 do begin
            msf(j,k) = msf(j,k) - msf(j,jpk-1)
          endfor
        endfor
endif

; mask the msf
msf=msf*mmsk

varname='Meridional Stream Function'
varunit='Sv'
vargrid='W'


if keyword_set(orca) then begin
        case orca of
                2:ixref=51-1
                .5:ixref=204-1
                .25:ixref=408-1
                else:ixref=0
        endcase
        ixref=ixref-ixminmesh
        ixref=ixref+key_shift
        ixref=ixref mod jpi
        if ixref lt 0 then ixref=ixref+jpi
endif else ixref=0
domdef,ixref,ixref,0,jpj-1,0,jpk-1,/index,grille=['T','W']
;terre=where(total(vmask()*mask,1) eq 0)
;if n_elements(valmask) eq 0 then valmask=1e20
;if terre[0] ne -1 then msf[terre]=valmask

;------------------------
return, msf
end