;+
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
;-
FUNCTION grad, field, direc
;
@common
;
   res = litchamp(field)
   taille=size(res)
   grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz
   if n_elements(valmask) EQ 0 then valmask = 1e20
   case strupcase(vargrid) of
      'T':BEGIN
         case direc of
            'x':BEGIN
               divi = e1u[premierx:dernierx, premiery:derniery]
               newmask = (umask())[premierx:dernierx, premiery:derniery, premierz:dernierz]
               vargrid = 'U'
               domdef, glamt[premierx, 0], glamu[dernierx, 0] $
                , gphit[0, premiery], gphiu[0, derniery], grille = ['T','U']
            END
            'y':BEGIN
               divi = e2v[premierx:dernierx, premiery:derniery]
               newmask = (vmask())[premierx:dernierx, premiery:derniery, premierz:dernierz]
               vargrid = 'V'
               domdef, glamt[premierx, 0], glamv[dernierx, 0] $
                , gphit[0, premiery], gphiv[0, derniery], grille = ['T','V']
            END
            'z':BEGIN
               divi = e3w[premierz:dernierz]
               newmask = mask
               vargrid = 'W'
            END
            ELSE : return, report('Bad definition of direction argument')
         ENDCASE
      END
      'W':BEGIN
         case direc of
            'x':divi = e1u[premierx:dernierx, premiery:derniery]
            'y':divi = e2v[premierx:dernierx, premiery:derniery]
            'z':BEGIN
               divi = e3t[premierz:dernierz]
               newmask = mask
               vargrid = 'T'
            END
            ELSE : return, report('Bad definition of direction argument')
         endcase
      END
      'U':BEGIN
         case direc of
            'x':BEGIN
               divi = (shift(e1t, -1, 0))[premierx:dernierx, premiery:derniery]
               newmask = tmask[premierx:dernierx, premiery:derniery, premierz:dernierz]
               vargrid = 'T'
               domdef, glamt[premierx, 0], glamu[dernierx] $
                , gphit[0, premiery], gphiu[0, derniery], grille = ['T','U']
            END
            'y':BEGIN
               divi = e2f[premierx:dernierx, premiery:derniery]
               newmask = (fmask())[premierx:dernierx, premiery:derniery, premierz:dernierz]
               vargrid = 'F'
               domdef, glamu[premierx, 0], glamf[dernierx, 0] $
                , gphiu[0, premiery], gphif[0, derniery], grille = ['U','F']
            END
            'z':BEGIN
               divi = e3w[premierz:dernierz]
               newmask = mask
               vargrid = 'W'
            END
            ELSE : return, report('Bad definition of direction argument')
         endcase
      END
      'V':BEGIN
         case direc of
            'x':BEGIN
               divi = e1f[premierx:dernierx, premiery:derniery]
               newmask = (fmask())[premierx:dernierx, premiery:derniery, premierz:dernierz]
               vargrid = 'F'
               domdef, glamv[premierx, 0], glamf[dernierx, 0] $
                , gphiv[0, premiery], gphif[0, derniery], grille = ['V','F']
            END
            'y':BEGIN
               divi = (shift(e2t, 0, -1))[premierx:dernierx, premiery:derniery]
               newmask = tmask[premierx:dernierx, premiery:derniery, premierz:dernierz]
               vargrid = 'T'
               domdef, glamt[premierx, 0], glamv[dernierx, 0] $
                , gphit[0, premiery], gphiv[0, derniery], grille = ['T','V']
            END
            'z':BEGIN
               divi = e3w[premierz:dernierz]
               newmask = mask
               vargrid = 'W'
            END
            ELSE : return, report('Bad definition of direction argument')
         endcase
      END
;       'F':BEGIN
;          case direc of
;             'x':divi = (shift(e1v, -1, 0))[premierx:dernierx, premiery:derniery]
;             'y':divi = (shift(e2u, 0, -1))[premierx:dernierx, premiery:derniery]
;             'z':divi = e3w[premierz:dernierz]
;             ELSE : return, report('Bad definition of direction argument')
;          endcase
;       END
      ELSE : return, report('Bad definition of vargrid')
   ENDCASE
   res = fitintobox(res)
   case 1 of
;----------------------------------------------------------------------------
;xy
;----------------------------------------------------------------------------
      taille[0] EQ 2:BEGIN
         earth = where(mask[*, *, premierz] EQ 0)
         if earth[0] NE -1 then res[earth] = !values.f_nan
         case direc of
            'x':BEGIN
               res = (shift(res, -1, 0)-res)/divi
               if key_periodique EQ 0 OR nx NE jpi THEN res[nx-1, *] = !values.f_nan
               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0)
            END
            'y':BEGIN
               res = (shift(res, 0, -1)-res)/divi
               res[*, ny-1] = !values.f_nan
               if vargrid EQ 'T' OR vargrid EQ 'U' then res =   shift(res, 0, 1)
            END
            ELSE : return,  report('Bad definition of direction argument for the type of array')
         ENDCASE
         earth = where(newmask[*, *, premierz] EQ 0)
         if earth[0] NE -1 then res[earth] = valmask
      END
;----------------------------------------------------------------------------
;xyt
;----------------------------------------------------------------------------
      taille[0] EQ 3 AND jpt NE 1:BEGIN
         earth = where(mask[*, *, premierz] EQ 0)
         if earth[0] NE -1 then BEGIN
            earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt))
            res[earth[*]] = !values.f_nan
         ENDIF
         divi = divi[*]#replicate(1, jpt)
         case direc of
            'x':BEGIN
               res = (shift(res, -1, 0, 0)-res)/divi
               if key_periodique EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
               if vargrid EQ 'T' OR vargrid EQ 'V'  then res = shift(res, 1, 0, 0)
            END
            'y':BEGIN
               res = (shift(res, 0, -1, 0)-res)/divi
               res[*, ny-1, *] = !values.f_nan
               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0)
            END
            ELSE : return,  report('Bad definition of direction argument for the type of array')
         ENDCASE
         earth = where(newmask[*, *, premierz] EQ 0)
         if earth[0] NE -1 then BEGIN
            earth = earth#replicate(1, jpt)+replicate(1, n_elements(earth))#(nx*ny*lindgen(jpt))
            res[earth] = valmask
         ENDIF
      END
;----------------------------------------------------------------------------
;xyz
;----------------------------------------------------------------------------
      taille[0] EQ 3 AND jpt EQ 1:BEGIN
         earth = where(mask EQ 0)
         if earth[0] NE -1 then res[earth] = !values.f_nan
         case direc OF
            'x':BEGIN
               divi = divi[*]#replicate(1, nz)
               res = (shift(res, -1, 0, 0)-res)/divi
               if key_periodique EQ 0 OR nx NE jpi THEN res[nx-1, *, *] = !values.f_nan
               if vargrid EQ 'T' OR vargrid EQ 'V' then res = shift(res, 1, 0, 0)
            END
            'y':BEGIN
               divi = divi[*]#replicate(1, nz)
               res = (shift(res, 0, -1, 0)-res)/divi
               res[*, ny-1, *] = !values.f_nan
               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0)
            END
            'z':BEGIN
               divi = reform(replicate(1, nx*ny)#divi, nx, ny, nz)
               if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz)
               if vargrid EQ 'W' THEN BEGIN
                  res = (shift(res, 0, 0, 1)-res)/divi
                  res[*, *, 0] = !values.f_nan
               ENDIF ELSE BEGIN
                  res = (res-shift(res, 0, 0, -1))/divi
                  res[*, *, nz-1] = !values.f_nan
               ENDELSE
               if earth[0] NE -1 then res[earth] = valmask
            END
         ENDCASE
      END
;------------------------------------------------------------
;----------------------------------------------------------------------------
;xyzt
;----------------------------------------------------------------------------
      taille[0] EQ 4:BEGIN
         earth = where((mask[*]#replicate(1, jpt)) EQ 0)
         if earth[0] NE -1 then res[earth] = !values.f_nan
         case direc OF
            'z':BEGIN
               divi = replicate(1, nx*ny)#divi
               divi = reform(divi[*]#replicate(1, jpt), nx, ny, nz, jpt, /over)
               if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt)
               if vargrid EQ 'W' THEN BEGIN
                  res = (shift(res, 0, 0, 1, 0)-res)/divi
                  res[*, *, 0, *] = !values.f_nan
               ENDIF ELSE BEGIN
                  res = (res-shift(res, 0, 0, -1, 0))/divi
                  res[*, *, nz-1, *] = !values.f_nan
               ENDELSE
               if earth[0] NE -1 then res[earth] = valmask
            END
         ENDCASE
      END
;------------------------------------------------------------
;------------------------------------------------------------
   endcase
   varname = 'grad of '+varname
   varunit = varunit+'/m'


;------------------------------------------------------------
   return, res
end