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