;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:div
;
; PURPOSE:calcule la divergence d'un champ 2D
;
; CATEGORY:calcule sur les matrices
;
; CALLING SEQUENCE:res=div(u,v)
;
; INPUTS:
;       u et v deux matrices representant les coordonnes d''un
;       champ de vecteur
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:res: une matrice 2d
;
; COMMON BLOCKS:
;	common.pro
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
; les matrices u et v peuvent de 2 a 4 dimensions.
; attention pour distinger les differents configurations de u et v
; (xy, xyz, xyt, xyzt), on regarde la variable du common 
;        -time qui contient le calendrier en jour julien d''IDL auquel
;        se rapportent u et v ansi que la variable 
;        -jpt qui est le nombre de pas de temps a considerer ds time.
; les tableaux u et v sont decoupes sur le meme domaine
; geographique. A cause du decalage des grilles T, U, V et F il est
; possiible que ces 2 tableaux n''aient pas la meme taille et se
; repportent a des indices differents. Si tel est le cas les tableaux
; sont redecoupes sur les indices qu'ils ont en commun et le dommaine
; est redefinit pour qu'il colle a ces indices communs.
; pour eviter ces redecoupes utiliser le mot cles /memeindice ds
; domdef.pro
;
; les points sur le bord du dessin sont mis a !values.f_nan
;
; EXAMPLE:
;
; MODIFICATION HISTORY:Guillaume Roullet (grlod@ipsl.jussieu.fr)
;                      Creation : printemps 1998
;                      Sebastien Masson (smasson@lodyc.jussieu.fr)
;                      adaptation pour marcher avec un domaine reduit
;                      12/1/2000
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION divfred, uu, vv
   tempsun = systime(1)         ; pour key_performance
@common
   u = uu
   v = vv

   date1 = time[0]
   if n_elements(jpt) EQ 0 then date2 = date1 ELSE date2 = time[jpt-1]
   if (size(u))[0] NE (size(v))[0] then return,  -1

;------------------------------------------------------------
; on trouve les points que u et v ont en communs
;------------------------------------------------------------
   indicexu = (lindgen(jpi))[premierxu:premierxu+nxu-1]
   indicexv = (lindgen(jpi))[premierxv:premierxv+nxv-1]
   indicex = inter(indicexu, indicexv)
   indiceyu = (lindgen(jpj))[premieryu:premieryu+nyu-1]
   indiceyv = (lindgen(jpj))[premieryv:premieryv+nyv-1]
   indicey = inter(indiceyu, indiceyv)
   nx = n_elements(indicex) 
   ny = n_elements(indicey)
   indice2d = lindgen(jpi, jpj)
   indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
   case 1 of
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyz
;----------------------------------------------------------------------------
      (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN 
;------------------------------------------------------------
; extraction de u et v sur le domaine qui convient
;------------------------------------------------------------
         case 1 of
            (size(u))[0] NE 3 OR (size(v))[0] NE 3: return,  -1
            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
               case 1 of
                  nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
                  nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
                  nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
                  nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
                  ELSE :
               endcase
            END
            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
            END
            ELSE:BEGIN 
               zdiv = -1
               GOTO, sortie
            end
            
         endcase
;------------------------------------------------------------
; calcul de la divergence
;------------------------------------------------------------
         zu = (e2u[indice2d])[*]#replicate(1, nzt)
         zu = reform(zu, nx, ny, nzt, /over)
         zu = temporary(u)*temporary(zu) $
          *(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]
         terreu = where(zu EQ 0)
;;;         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
;
         zv = (e1v[indice2d])[*]#replicate(1, nzt)
         zv = reform(zv, nx, ny, nzt, /over)
         zv = temporary(v)*temporary(zv) $
          *(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]
         terrev = where(zv EQ 0)
;;;         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
;
         zdiv = 1e6/(e1t[indice2d]*e2t[indice2d])
         zdiv = (zdiv)[*]#replicate(1, nzt)
         zdiv = reform(zdiv, nx, ny, nzt, /over)
         zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0)) $
          *tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt]
;------------------------------------------------------------
; mise a !values.f_nan de la bordure
;------------------------------------------------------------
         if  NOT keyword_set(key_periodique) OR nx NE jpi then begin
            zdiv(0, *, *) = !values.f_nan
            zdiv(nx-1, *, *) = !values.f_nan
         endif
;;         zdiv(*, 0, *) = !values.f_nan
;;         zdiv(*, ny-1, *) = !values.f_nan
;
         zdiv = temporary(zdiv)
;
         if n_elements(valmask) EQ 0 THEN valmask = 1e20
         terre =  where(tmask[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, premierzt:dernierzt] EQ 0)
         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
;------------------------------------------------------------
; pour le trace graphique
;------------------------------------------------------------
         vargrid = 'T'
         varname = 'div'
         varunits = '1e6*s-1'
         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t']
         if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan, boite = boite)
      END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyt
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
      date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN 
;------------------------------------------------------------
; extraction de u et v sur le domaine qui convient
;------------------------------------------------------------
         case 1 of
            (size(u))[0] NE 3 OR (size(v))[0] NE 3: return,  -1
            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
               case 1 of
                  nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
                  nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
                  nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
                  nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
                  ELSE :
               endcase
            END
            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
               u = u[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
               v = v[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, *]
            END
            ELSE:return,  -1
         endcase
;------------------------------------------------------------
; calcul de la divergence
;------------------------------------------------------------
         zu = e2u[indice2d]*(umask())[indice2d+jpi*jpj*premierzt]
         terreu = where(zu EQ 0)
;;;         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
         zu = (zu)[*]#replicate(1, jpt)
         zu = reform(zu, nx, ny, jpt, /over)
         zu = temporary(u)*temporary(zu)
;
         zv = e1v[indice2d]*(vmask())[indice2d+jpi*jpj*premierzt]
         terrev = where(zv EQ 0)
;;;         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
         zv = (zv)[*]#replicate(1, jpt)
         zv = reform(zv, nx, ny, jpt, /over)
         zv = temporary(v)*temporary(zv)
;
         zdiv = 1e6*tmask[indice2d+jpi*jpj*premierzt]/(e1t[indice2d]*e2t[indice2d])
         zdiv = (zdiv)[*]#replicate(1, jpt)
         zdiv = reform(zdiv, nx, ny, jpt, /over)
         terre =  where(zdiv EQ 0)
         zdiv = temporary(zdiv)*(zu-shift(zu, 1, 0, 0)+zv-shift(zv, 0, 1, 0))
;------------------------------------------------------------
; mise a !values.f_nan de la bordure
;------------------------------------------------------------
         if  NOT keyword_set(key_periodique) OR nx NE jpi then begin
            zdiv(0, *, *) = !values.f_nan
            zdiv(nx-1, *, *) = !values.f_nan
         endif
;;         zdiv(*, 0, *) = !values.f_nan
;;         zdiv(*, ny-1, *) = !values.f_nan
;
         if n_elements(valmask) EQ 0 THEN valmask = 1e20
         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
;------------------------------------------------------------
; pour le trace graphique
;------------------------------------------------------------
         vargrid = 'T'
         varname = 'div'
         varunits = '1e6*s-1'
         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t']
         if keyword_set(direc) then  zdiv = grossemoyenne(zdiv,direc,/nan, boite = boite)
      END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyzt
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
      date1 NE date2 AND (size(u))[0] EQ 4:BEGIN 
         return, report('non code!')
      END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xy
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
      ELSE:BEGIN                ;xy
         indice3d = lindgen(jpi, jpj, jpk)
         indice3d = indice3d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, niveau -1]
;------------------------------------------------------------
; extraction de u et v sur le domaine qui convient
;------------------------------------------------------------
         case 1 of
            (size(u))[0] NE 2 OR (size(v))[0] NE 2: BEGIN 
               zdiv = -1
               GOTO, sortie
            end
            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
               case 1 of
                  nxu NE nx:if indicex[0] EQ premierxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *]
                  nxv NE nx:if indicex[0] EQ premierxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
                  nyu NE ny:if indicey[0] EQ premieryu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny]
                  nyv NE ny:if indicey[0] EQ premieryv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
                  ELSE :
               endcase
            END
            (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND $
             (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj:BEGIN 
               u = u[indice2d]
               v = v[indice2d]
            END
            ELSE:return,  -1
         endcase
;------------------------------------------------------------
; calcul de la divergence
;------------------------------------------------------------
         zu = temporary(u)*e2u[indice2d]*(umask())[indice3d]
         terreu = where(zu EQ 0)
;;;         if terreu[0] NE -1 then zu[temporary(terreu)] = !values.f_nan
         zv = temporary(v)*e1v[indice2d]*(vmask())[indice3d]
         terrev = where(zv EQ 0)
;;;         if terrev[0] NE -1 then zv[temporary(terrev)] = !values.f_nan
         zdiv = zu-shift(zu, 1, 0)+zv-shift(zv, 0, 1)
         zdiv = temporary(zdiv)*tmask[indice3d]/(e1t[indice2d]*e2t[indice2d])
;------------------------------------------------------------
; mise a !values.f_nan de la bordure
;------------------------------------------------------------
         if  NOT keyword_set(key_periodique) OR nx NE jpi then begin
            zdiv(0, *) = !values.f_nan
            zdiv(nx-1, *) = !values.f_nan
         endif
;;         zdiv(*, 0) = !values.f_nan
;;         zdiv(*, ny-1) = !values.f_nan
;
         zdiv = temporary(zdiv)*1e6
;
         if n_elements(valmask) EQ 0 THEN valmask = 1e20
         terre =  where(tmask[indice3d] EQ 0)
         if terre[0] NE -1 then zdiv[temporary(terre)] = valmask
;------------------------------------------------------------
; pour le trace graphique
;------------------------------------------------------------
         vargrid = 'T'
         varname = 'div'
         varunits = '1e6*s-1'
         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], prof1, prof2, grille = ['t']
         if keyword_set(direc) then  zdiv = moyenne(zdiv,direc,/nan, boite = boite)

      END
   endcase
   
;------------------------------------------------------------
sortie:
   if keyword_set(key_performance) THEN print, 'temps div', systime(1)-tempsun 
   
   return, zdiv
end