;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:
;
; PURPOSE:
;
; CATEGORY:
;
; CALLING SEQUENCE:
; 
; INPUTS:
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:
;
; COMMON BLOCKS:common.pro
;
; SIDE EFFECTS:
;
; RESTRICTIONS:
;
; EXAMPLE:
;
; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
;
; @todo seb: remplir les truc!
;
;
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION grad, field, direc
;
  compile_opt idl2, strictarrsubs
;
@common
;------------------------------------------------------------
;
   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
     return, report(['This version of grad is based on Arakawa C-grid.' $
                     , 'U and V grids must therefore be defined'])
;
   res = litchamp(field)
   taille=size(res)
   grille, mask, glam, gphi, gdep, nx, ny,nz,firstx,firsty,firstz,lastx, lasty, lastz     
   if n_elements(valmask) EQ 0 then valmask = 1e20
   case strupcase(vargrid) of
      'T':BEGIN
         case direc of
            'x':BEGIN 
               divi = e1u[firstx:lastx, firsty:lasty]
               newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
               vargrid = 'U'
               domdef, glamt[firstx, 0], glamu[lastx, 0] $
                , gphit[0, firsty], gphiu[0, lasty], gridtype = ['T','U']
            END
            'y':BEGIN
               divi = e2v[firstx:lastx, firsty:lasty]
               newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
               vargrid = 'V'
               domdef, glamt[firstx, 0], glamv[lastx, 0] $
                , gphit[0, firsty], gphiv[0, lasty], gridtype = ['T','V']
            END
            'z':BEGIN
               divi = e3w[firstz:lastz]
               newmask = mask
               vargrid = 'W'
            END
            ELSE:return, report('Bad definition of direction argument')
         ENDCASE
      END
      'W':BEGIN
         case direc of
            'x':divi = e1u[firstx:lastx, firsty:lasty]
            'y':divi = e2v[firstx:lastx, firsty:lasty]
            'z':BEGIN
               divi = e3t[firstz:lastz]
               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))[firstx:lastx, firsty:lasty]
               newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
               vargrid = 'T'
               domdef, glamt[firstx, 0], glamu[lastx] $
                , gphit[0, firsty], gphiu[0, lasty], gridtype = ['T','U']
            END
            'y':BEGIN
               divi = e2f[firstx:lastx, firsty:lasty]
               newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
               vargrid = 'F'
               domdef, glamu[firstx, 0], glamf[lastx, 0] $
                , gphiu[0, firsty], gphif[0, lasty], gridtype = ['U','F']
            END
            'z':BEGIN
               divi = e3w[firstz:lastz]
               newmask = mask
               vargrid = 'W'
            END
            ELSE:return, report('Bad definition of direction argument')
         endcase
      END
      'V':BEGIN
         case direc of
            'x':BEGIN
               divi = e1f[firstx:lastx, firsty:lasty]
               newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
               vargrid = 'F'
               domdef, glamv[firstx, 0], glamf[lastx, 0] $
                , gphiv[0, firsty], gphif[0, lasty], gridtype = ['V','F']
            END
            'y':BEGIN
               divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty]
               newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
               vargrid = 'T'
               domdef, glamt[firstx, 0], glamv[lastx, 0] $
                , gphit[0, firsty], gphiv[0, lasty], gridtype = ['T','V']
            END
            'z':BEGIN
               divi = e3w[firstz:lastz]
               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))[firstx:lastx, firsty:lasty]
;             'y':divi = (shift(e2u, 0, -1))[firstx:lastx, firsty:lasty]
;             'z':divi = e3w[firstz:lastz]
;             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[*, *, firstz] 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_periodic 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[*, *, firstz] 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[*, *, firstz] 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_periodic 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[*, *, firstz] EQ 0)
         if earth[0] NE -1 THEN res[earth] = valmask
      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_periodic 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
            'x':BEGIN 
               divi = divi[*]#replicate(1, nz*jpt)
               res = (shift(res, -1, 0, 0, 0)-res)/divi
               if key_periodic 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, 0)
            END
            'y':BEGIN 
               divi = divi[*]#replicate(1, nz*jpt)
               res = (shift(res, 0, -1, 0, 0)-res)/divi
               res[*, ny-1, *, *] = !values.f_nan
               if vargrid EQ 'T' OR vargrid EQ 'U'  then res = shift(res, 0, 1, 0, 0)        
            END
            '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
            END
         ENDCASE
         if earth[0] NE -1 then res[earth] = valmask
      END
;------------------------------------------------------------
;------------------------------------------------------------
   endcase
   varname = 'grad of '+varname
   varunit = varunit+'/m'


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