;+
;;
;; Resolution de l'equation A.X=b ou A est une matrice creuse
;; et (x,b) sont des champs horizontaux
;;
;; Creation : printemps 1998 / G. Roullet
;; 
;-
FUNCTION inversion, b, GUESS = x, F = f, TOL = tol, _extra = extra
@common
@com_eg
COMMON laplaciens, lapla_t, lapla_f

   IF (where(b NE 0))[0] EQ -1 THEN return, fltarr(jpi, jpj)

   IF n_elements(x) EQ 0 THEN x = fltarr(jpi, jpj)

   IF keyword_set(tol) EQ 0 THEN tol = 1e-3

   print, '          Solver precision:', tol

   IF n_elements(f) EQ 0 THEN BEGIN 
      IF n_elements(lapla_t) EQ 0 THEN build_laplacien_t   
      A = lapla_t 
   ENDIF ELSE BEGIN 
      IF n_elements(lapla_f) EQ 0 THEN build_laplacien_f   
      A = lapla_f
   ENDELSE 

   zb = double(b)
   zx = double(x)

   zr = linbcg(A, zb(*), zx(*), tol = tol, iter = iter, /double, _extra = extra)
   print, '          Number of iterations :', iter

   zr = boundperio(reform(zr, jpi, jpj), ff = f)
   
   return, zr
END