;+
;
; @file_comments
; compute the gradient of a variable
;
; @categories
; Calculation
;
; @param FIELD
; The field for which we want to compute the gradient. A 2D (xy),
; 3D (xyz or yt) or 4D (xyzt) array or a structure readable by litchamp
; and containing a 2D (xy), 3D (xyz or yt) or 4D (xyzt) array.
; note that the dimension of the array must suit the domain dimension.
;
; @param DIREC {type=scalar string}
; the gradient direction: 'x', 'y', 'z'
;
; @returns
; the gradient of the input data with the same size 2D, 3D or 4D array
;
; @uses
; cm_4cal
; cm_4data
; cm_4mmesh
;
; @restrictions
; - Works only for Arakawa C-grid.
; - When computing the gradient, the result is not on the same grid point
;   than the input data. In consequence, we update, vargrid and the grid position
;   parameters (firstx[tuvf], lastx[tuvf], nx[tuvf], firsty[tuvf], lasty[tuvf],
;   ny[tuvf], firstz[tw], lastz[tw], nz[tw]).
; - points that cannot be computed (domain boundaries, coastline) are set to NaN
; - the common variable jpt is used to differ xyz (jpt=1) and xyt (jpt\=1) cases.
;
; @examples
; IDL> \@tst_initorca2
; IDL> plt, grad({arr:gphit,g:'T'}, 'x')
; IDL> plt, grad({arr:gphit,g:'T'}, 'y')
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;
; @version
; $Id$
;
;-
;
FUNCTION grad, field, direc
;
  compile_opt idl2, strictarrsubs
;
@cm_4cal   ; for jpt
@cm_4data  ; for varname, vargrid, vardate, varunit, valmask
@cm_4mesh
;------------------------------------------------------------
;
  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)
  szres = size(res)
  grille, mask, glam, gphi, gdep, nx, ny, nz $
          , firstx, firsty, firstz, lastx, lasty, lastz
;
  if n_elements(valmask) EQ 0 then valmask = 1.e20
  varname = 'grad of '+varname
  varunit = varunit+'/m'
  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'
          firstxu = firstxt & lastxu = lastxt & nxu = nxt
          firstyu = firstyt & lastyu = lastyt & nyu = nyt
        END
        'y':BEGIN
          divi = e2v[firstx:lastx, firsty:lasty]
          newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
          vargrid = 'V'
          firstxv = firstxt & lastxv = lastxt & nxv = nxt
          firstyv = firstyt & lastyv = lastyt & nyv = nyt
        END
        'z':BEGIN
          divi = e3w[firstz:lastz]
          newmask = mask
          vargrid = 'W'
          firstzw = firstzt & lastzw = lastzt & nzw = nzt
        END
        ELSE:return, report('Bad definition of direction argument')
      ENDCASE
    END
    'W':BEGIN
      case direc of
        'x':BEGIN
          divi = e1u[firstx:lastx, firsty:lasty]
          newmask = (umask())[firstx:lastx, firsty:lasty, firstz:lastz]
          vargrid = 'U'
          firstxu = firstxt & lastxu = lastxt & nxu = nxt
          firstyu = firstyt & lastyu = lastyt & nyu = nyt
        END
        'y':BEGIN
          divi = e2v[firstx:lastx, firsty:lasty]
          newmask = (vmask())[firstx:lastx, firsty:lasty, firstz:lastz]
          vargrid = 'V'
          firstxv = firstxt & lastxv = lastxt & nxv = nxt
          firstyv = firstyt & lastyv = lastyt & nyv = nyt
        END
        'z':BEGIN
          divi = e3t[firstz:lastz]
          newmask = mask
          vargrid = 'T'
          firstzt = firstzw & lastzt = lastzw & nzt = nzw
        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'
          firstxt = firstxu & lastxt = lastxu & nxt = nxu
          firstyt = firstyu & lastyt = lastyu & nyt = nyu
        END
        'y':BEGIN
          divi = e2f[firstx:lastx, firsty:lasty]
          newmask = (fmask())[firstx:lastx, firsty:lasty, firstz:lastz]
          vargrid = 'F'
          firstxf = firstxu & lastxf = lastxu & nxf = nxu
          firstyf = firstyu & lastyf = lastyu & nyf = nyu
        END
        'z':BEGIN
          divi = e3w[firstz:lastz]
          newmask = mask
          vargrid = 'W'
          firstzw = firstzt & lastzw = lastzt & nzw = nzt
        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'
          firstxf = firstxv & lastxf = lastxv & nxf = nxv
          firstyf = firstyv & lastyf = lastyv & nyf = nyv
        END
        'y':BEGIN
          divi = (shift(e2t, 0, -1))[firstx:lastx, firsty:lasty]
          newmask = tmask[firstx:lastx, firsty:lasty, firstz:lastz]
          vargrid = 'T'
          firstxt = firstxv & lastxt = lastxv & nxt = nxv
          firstyt = firstyv & lastyt = lastyv & nyt = nyv
        END
        'z':BEGIN
          divi = e3w[firstz:lastz]
          newmask = mask
          vargrid = 'W'
          firstzw = firstzt & lastzw = lastzt & nzw = nzt
        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
      return, report('F grid: case not coded, please contact SAXO team...')
    END
    ELSE:return, report('Bad definition of vargrid')
  ENDCASE
  res = fitintobox(temporary(res))
  IF n_elements(res) EQ 1 AND res[0] EQ -1 THEN return, res
  case 1 of
;----------------------------------------------------------------------------
;xy
;----------------------------------------------------------------------------
    szres[0] EQ 2:BEGIN
      land = where((temporary(mask))[*, *, firstz] EQ 0)
      if land[0] NE -1 then res[temporary(land)] = !values.f_nan
      case direc of
        'x':BEGIN
          res = (shift(res, -1, 0)-res)/temporary(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(temporary(res), 1, 0)
        END
        'y':BEGIN
          res = (shift(res, 0, -1)-res)/temporary(divi)
          res[*, ny-1] = !values.f_nan
          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1)
        END
        ELSE:return,  report('Bad definition of direction argument for the type of array')
      ENDCASE
      land = where((temporary(newmask))[*, *, firstz] EQ 0)
      if land[0] NE -1 then res[temporary(land)] = valmask
    END
;----------------------------------------------------------------------------
;xyt
;----------------------------------------------------------------------------
    szres[0] EQ 3 AND jpt NE 1:BEGIN
      land = where((temporary(mask))[*, *, firstz] EQ 0, cnt)
      if land[0] NE -1 then BEGIN
        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt))
        res[temporary(land)] = !values.f_nan
      ENDIF
      divi = (temporary(divi))[*]#replicate(1., jpt)
      case direc of
        'x':BEGIN
          res = (shift(res, -1, 0, 0)-res)/temporary(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(temporary(res), 1, 0, 0)
        END
        'y':BEGIN
          res = (shift(res, 0, -1, 0)-res)/temporary(divi)
          res[*, ny-1, *] = !values.f_nan
          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0)
        END
        ELSE:return,  report('Bad definition of direction argument for the type of array')
      ENDCASE
      land = where((temporary(newmask))[*, *, firstz] EQ 0, cnt)
      if land[0] NE -1 then BEGIN
        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*lindgen(jpt))
        res[temporary(land)] = valmask
      ENDIF
    END
;----------------------------------------------------------------------------
;xyz
;----------------------------------------------------------------------------
    szres[0] EQ 3 AND jpt EQ 1:BEGIN
      land = where(mask EQ 0)
      if land[0] NE -1 then res[temporary(land)] = !values.f_nan
      case direc OF
        'x':BEGIN
          divi = (temporary(divi))[*]#replicate(1., nz)
          res = (shift(res, -1, 0, 0)-res)/temporary(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(temporary(res), 1, 0, 0)
        END
        'y':BEGIN
          divi = (temporary(divi))[*]#replicate(1., nz)
          res = (shift(res, 0, -1, 0)-res)/temporary(divi)
          res[*, ny-1, *] = !values.f_nan
          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0)
        END
        'z':BEGIN
          divi = replicate(1., nx*ny)#(temporary(divi))[*]
          if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, /overwrite)
          if vargrid EQ 'W' THEN BEGIN
            res = (shift(res, 0, 0, 1)-res)/temporary(divi)
            res[*, *, 0] = !values.f_nan
          ENDIF ELSE BEGIN
            res = (res-shift(res, 0, 0, -1))/temporary(divi)
            res[*, *, nz-1] = !values.f_nan
          ENDELSE
        END
      ENDCASE
      land = where(temporary(newmask) EQ 0)
      if land[0] NE -1 then res[temporary(land)] = valmask
    END
;----------------------------------------------------------------------------
;xyzt
;----------------------------------------------------------------------------
    szres[0] EQ 4:BEGIN
      land = where(temporary(mask) EQ 0, cnt)
      if land[0] NE -1 then BEGIN
        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*nz*lindgen(jpt))
        res[temporary(land)] = !values.f_nan
      ENDIF
      case direc OF
        'x':BEGIN
          divi = (temporary(divi))[*]#replicate(1., nz*jpt)
          res = (shift(res, -1, 0, 0, 0)-res)/temporary(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(temporary(res), 1, 0, 0, 0)
        END
        'y':BEGIN
          divi = (temporary(divi))[*]#replicate(1., nz*jpt)
          res = (shift(res, 0, -1, 0, 0)-res)/temporary(divi)
          res[*, ny-1, *, *] = !values.f_nan
          if vargrid EQ 'T' OR vargrid EQ 'U' then res = shift(temporary(res), 0, 1, 0, 0)
        END
        'z':BEGIN
          divi = replicate(1., nx*ny)#(temporary(divi))[*]
          divi = (temporary(divi))[*]#replicate(1L, jpt)
          if nx EQ 1 OR ny EQ 1 then res = reform(res, nx, ny, nz, jpt, /overwrite)
          if vargrid EQ 'W' THEN BEGIN
            res = (shift(res, 0, 0, 1, 0)-res)/temporary(divi)
            res[*, *, 0, *] = !values.f_nan
          ENDIF ELSE BEGIN
            res = (res-shift(res, 0, 0, -1, 0))/temporary(divi)
            res[*, *, nz-1, *] = !values.f_nan
          ENDELSE
        END
      ENDCASE
      land = where(newmask EQ 0, cnt)
      if land[0] NE -1 then BEGIN
        land = (temporary(land))#replicate(1L, jpt) + replicate(1L, cnt)#(nx*ny*nz*lindgen(jpt))
        res[temporary(land)] = valmask
      ENDIF
    END
;------------------------------------------------------------
;------------------------------------------------------------
    ELSE:return, report('input array must have 2, 3 or 4 dimensions')
  ENDCASE
;------------------------------------------------------------


;------------------------------------------------------------
  return, res
END