;+
;;
;; Calcul de la Fonction de Courant Barotrope 
;; a partir du champ de vitesse (u,v)
;;
;; creation : printemps 98 G. Roullet
;;
;-

PRO diag_bsf_mu, hu, hv, mu
@common
@com_eg
;;
;; calcul adapte a la grille ORCA
;;
   mu = fltarr(14)
; 0/ Antartic
;   mu(0) = total( (hu*e2u*umaskr(*, *, 0))[149, 5:35])
   mu(0) = total( (hu*e2u*umaskr(*, *, 0))[151, 14:48])
; 1/ America
;   mu(1) = -total( (hv*e1v*vmaskr(*, *, 0))[110:150, 60])
   mu(1) = -total( (hv*e1v*vmaskr(*, *, 0))[110:150, 60])
; 3/ Australia
;   mu(3) = total( (hu*e2u*umaskr(*, *, 0))[17, 40:49])
   mu(3) = total( (hu*e2u*umaskr(*, *, 0))[20, 50:61])
; 2/ New Zealand
;   mu(2) = mu(3)+ total( (hv*e1v*vmaskr(*, *, 0))[30:40, 25])
   mu(2) = mu(3)+ total( (hv*e1v*vmaskr(*, *, 0))[35:49, 44])
; 4/ Madagascar
;   mu(4) = total( (hv*e1v*vmaskr(*, *, 0))[150:159, 40])
   mu(4) = total( (hv*e1v*vmaskr(*, *, 0))[160:164, 56])
; 5/ New Guinea
;   mu(5) = - total( (hu*e2u*umaskr(*, *, 0))[16, 50:55])
   mu(5) = total( (hv*e2v*vmaskr(*, *, 0))[23:33, 61])
; 6/ Celebes
   mu(6) = total( (hv*e1v*vmaskr(*, *, 0))[13:22, 64])
; 7/ Borneo
;   mu(6) = total( (hv*e1v*vmaskr(*, *, 0))[7:12, 60])
   mu(7) = total( (hv*e1v*vmaskr(*, *, 0))[14:17, 67])
; 8/ Philippines
;   mu(7) = total( (hv*e1v*vmaskr(*, *, 0))[10:16, 90])
   mu(8) = total( (hv*e1v*vmaskr(*, *, 0))[14:22, 86])
; 9/ Tasmania
   mu(9) =  mu(3)+ total( (hu*e1u*umaskr(*, *, 0))[34, 41:45])
; 10/ Cuba
;   mu(8) = - total( (hv*e1v*vmaskr(*, *, 0))[101:127, 86])
   mu(10) = - total( (hv*e1v*vmaskr(*, *, 0))[105:134, 90])
; 11/ Island
;   mu(9) = - total( (hv*e1v*vmaskr(*, *, 0))[117:130, 124])
   mu(11) = - total( (hv*e1v*vmaskr(*, *, 0))[132:143, 120])
; 12/ Spitzberg
;   mu(10) = - total( (hv*e1v*vmaskr(*, *, 0))[110:120, 137])
   mu(12) = - total( (hv*e1v*vmaskr(*, *, 0))[140:160, 135])
; 13/ Japan
   mu(13) = total( (hv*e1v*vmaskr(*, *, 0))[25:131, 100])
end

PRO diag_bsf_lec, bsf
@common
@com_eg
   ;;
   ;; Lecture du fichier contenant les fonctions de courant associe a
   ;; chaque ile (pour ORCA)
   ;;
   ;; Creation : printemps 1998 / G. Roullet
   ;; modification : juin 99
   ;;
   jpisl=14
   bsf = fltarr(jpi, jpj, jpisl)
   z2d = fltarr(jpi, jpj)
   fichier = 'islands.nc'
   FOR k = 1, jpisl DO BEGIN 
      ile = 'island'+strtrim(string(k), 2)
      z2d = nc_get(fichier, ile, NCDF_DB =hom_idl+'grids/')
      bsf(*, *, k-1) = shift(z2d, 0, 0)
   ENDFOR 
end


FUNCTION diag_bsf, u, v, FLUX = flux, _EXTRA = extra
   ;;
   ;;  Routine principale
   ;;
@common
COMMON bsf_lec, bsfislands

   jpisl=14
   IF keyword_set(flux) EQ 0 THEN BEGIN 
;
;  integration sur la verticale
;
      hu = moyz(u, /int)
      hv = moyz(v, /int)
   ENDIF ELSE BEGIN 
;
;  on suppose que u et v sont deja integres verticalement
;
      hu = u
      hv = v
   ENDELSE 
   
 print, '     Inversion...'
   bsf = inversion( curl2(hu, hv), /f, _extra = extra)
;
;  calcul des potentiels entre chaque ile
;
 print, '     Potential...'
   diag_bsf_mu, hu, hv, mu
;
;  lecture des fonctions courants unite autour de chaque ile
;
   IF n_elements(bsfislands) EQ 0 THEN diag_bsf_lec, bsfislands
   
   FOR k = 0, jpisl-1 DO bsf = bsf+mu(k)*bsfislands(*, *, k)

   return, bsf

END