;+
;;
;;  Calcul de la Fonction de Courant a travers les lignes J
;;
;;
;;  creation : 12/02/99 M Imbard
;;  Modification : 12/11/98 G. Roullet 
;;                 18/11/98 G. Roullet
;-
PRO msf_mult, v, msf, msfmsk, MSK = msk
@common
;COMMON moy, zonmsk, mermsk, jpl
IF n_elements(jpl) EQ 0 THEN jpl = jpj
;
;  jpl bandes de latitudes
;
fm     = fltarr(jpl, jpk, jpt)
msf    = fltarr(jpl, jpk, jpt)
msfmsk = fltarr(jpl, jpk, jpt)
zv = v
z  = zv
ze1v   = zv
ze3v   = zv
vert = replicate(1, jpk, jpt)
;
;  Masque de u pour un calcul de msf sur le sous domaine msk(jpi,jpj)
;
IF n_elements(msk) NE 0 THEN BEGIN 
; msk est souvent vide
; donc cette procedure est ineffective.
      zvmask = boundperio( (msk +shift(msk, 0, -1) ) < 1 )
      zvmask = zvmask(*)#vert
      zv = zv*zvmask
ENDIF 
;
;  Ecriture "tricky" optimisee... si vous avez compris comment
;  on calculait un flux a partir de champ 2D, vous avez fait le
;  plus dur : ici on generalise a 3D d''ou l''utilisation de #vert
; 
;  calcul du flux
;
FOR i = 0, jpi-1 DO BEGIN for j= 0, jpj-1 do begin ze1v(i,j,*,*) = replicate(e1v(i,j),jpk*jpt) & endfor &endfor
FOR k = 0, jpk-1 DO BEGIN ze3v(*,*,k,*)=replicate(e3t(k),jpi*jpj*jpt) & endfor
;
z= -v*ze1v*ze3v
;la multiplication ne fait plus disparaitre la dimension temporelle
;maintenant que ze1v et ze3v ont ete modifie depuis msf_simple.pro 
; integration zonale du flux !
;
fm = total(z, 1) 
;
;  calcul de la msf en integrant depuis le fond
;
FOR k = jpk-2, 0, -1 DO begin msf[*, k,*] = msf[*, k+1,*]-fm(*, k,*) & endfor
;
;  msfmsk est le masque associe a msf (utilise pour les graphiques)
;
msfmsk = fm NE 0.
;
END