;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
;
; @file_comments
; Calculate the vertical component of the curl of a field of horizontal vectors
;
; @categories
; Calculation on matrixes
; 
; @param UU
; Matrix representing coordinates of a field of vectors
;
; @param VV 
; Matrix representing coordinates of a field of vectors
;
; @returns RES
; A 2d matrix
;
; @uses
; common.pro
;
; @restrictions
; U and V matrixes can be 2 or 4d.
; Beware, to discern differents configuration of U and V (xy, xyz, xyt, xyzt), 
; we look at the variable of the common 
;        -time which contain the calendar in IDL julian days to which U and 
; V refered to, in the same way as the variable 
;        -jpt which is the number of time's step to consider in time.
; U and V arrays ae cut in the same geographic domain. Because of the gap of 
; T, U, V and F grids, it is possible that these two arrays hase not the same 
; size and refered to different indexes. In this case, arrays are recut on 
; common indexesand the domain is redifined to match with these common 
; indexes. To avoid these recuts, use the keyword /memeindice in domdef.pro
;
;
; Points on the drawing edge are at !values.f_nan 
;
; @history 
; Guillaume Roullet (grlod@ipsl.jussieu.fr)
; 
;                      Sebastien Masson (smasson@lodyc.jussieu.fr)
;                      adaptation pour marcher avec un domaine reduit
;
;                      21/5/1999: valeurs manquantes a !values.f_nan
;periodicite 
; 
; @version
; $Id$
;
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION curl, uu, vv
;
  compile_opt idl2, strictarrsubs
;
@common
   tempsun = systime(1)         ; To key_performance
;
   IF finite(glamu[0])*finite(gphiu[0])*finite(glamv[0])*finite(gphiv[0]) EQ 0 THEN $
     return, report(['This version of curl is based on Arakawa C-grid.' $
                     , 'U and V grids must therefore be defined'])
;
   u = litchamp(uu)
   v = litchamp(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

;------------------------------------------------------------
; We find common points between U and V
;------------------------------------------------------------
   indicexu = (lindgen(jpi))[firstxu:firstxu+nxu-1]
   indicexv = (lindgen(jpi))[firstxv:firstxv+nxv-1]
   indicex = inter(indicexu, indicexv)
   indiceyu = (lindgen(jpj))[firstyu:firstyu+nyu-1]
   indiceyv = (lindgen(jpj))[firstyv:firstyv+nyv-1]
   indicey = inter(indiceyu, indiceyv)
   nx = n_elements(indicex) 
   ny = n_elements(indicey)
   case 1 of
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyz
;----------------------------------------------------------------------------
      (size(u))[0] EQ 3 AND date1 EQ date2 :BEGIN 
         indice2d = lindgen(jpi, jpj)
         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
;------------------------------------------------------------
; extraction of U and V on the appropriated domain
;------------------------------------------------------------
         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 firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *]
               nxv NE nx:if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
               nyu NE ny:if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *]
               nyv NE ny:if indicey[0] EQ firstyv 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
;------------------------------------------------------------
; calculation of the curl
;------------------------------------------------------------
         coefu = (e1u[indice2d])[*]#replicate(1, nzt)
         coefu = reform(coefu, nx, ny, nzt, /over)
         coefu = coefu*(umask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
         terreu = where(coefu EQ 0)
         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
         
         coefv = (e2v[indice2d])[*]#replicate(1, nzt)
         coefv = reform(coefv, nx, ny, nzt, /over)
         coefv = coefv*(vmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
         terrev = where(coefv EQ 0)
         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan

         tabf = (fmask())[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
         div = (e1f[indice2d]*e2f[indice2d])[*]#replicate(1, nzt)
         div = reform(div, nx, ny, nzt, /over)
         tabf = tabf/div
;
         zu = u*temporary(coefu)
         zv = v*temporary(coefv)

         psi =  (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
         psi = tabf*psi
;------------------------------------------------------------
; Edging put at !values.f_nan
;------------------------------------------------------------
         if NOT keyword_set(key_periodic)  OR nx NE jpi then begin
            psi[0, *, *] = !values.f_nan
            psi[nx-1, *, *] = !values.f_nan
         endif
         psi[*, 0, *] = !values.f_nan
         psi[*, ny-1, *] = !values.f_nan
;
         if n_elements(valmask) EQ 0 THEN valmask = 1e20
         terref =  where(tabf EQ 0)
         if terref[0] NE -1 then psi[temporary(terref)] = valmask
;------------------------------------------------------------
; For the graphic drawing
;------------------------------------------------------------
         domdef, (glagmt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f']
         if keyword_set(direc) then psi = moyenne(psi,direc,/nan)

;----------------------------------------------------------------------------
      END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyt
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
      date1 NE date2 AND (size(u))[0] EQ 3 :BEGIN 
         indice2d = lindgen(jpi, jpj)
         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
;------------------------------------------------------------
; extraction of U and V on the appropriated domain
;------------------------------------------------------------
         case 1 of
            (size(u))[1] EQ nxu AND (size(u))[2] EQ nyu AND $
             (size(v))[1] EQ nxv AND (size(v))[2] EQ nyv:BEGIN 
               if nxu NE nx then $
                if indicex[0] EQ firstxu then u = u[0:nx-1, *, *] ELSE u = u[1: nx, *, *] 
               IF nxv NE nx THEN $
                if indicex[0] EQ firstxv then v = v[0:nx-1, *, *] ELSE v = v[1: nx, *, *]
               IF nyu NE ny THEN $
                if indicey[0] EQ firstyu then u = u[*, 0:ny-1, *] ELSE u = u[*, 1: ny, *] 
               IF nyv NE ny THEN $
                if indicey[0] EQ firstyv then v = v[*, 0:ny-1, *] ELSE v = v[*, 1: ny, *]
            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 
               print, 'problemes d''adequation entre la taille du domaine et la taille des matrices necessaires a tracer des vecteurs'
               return, -1
            end
         endcase
;----------------------------------------------------------------------------
; Calculation of the curl
;----------------------------------------------------------------------------
         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
         terreu = where(coefu EQ 0)
         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
         coefu = temporary(coefu[*])#replicate(1, jpt)
         coefu = reform(coefu, nx, ny, jpt, /over)
;
         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
         terrev = where(coefv EQ 0)
         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
         coefv = temporary(coefv[*])#replicate(1, jpt)
         coefv = reform(coefv, nx, ny, jpt, /over)
;
         tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d])
         tabf = temporary(tabf[*])#replicate(1, jpt)
         tabf = reform(tabf, nx, ny, jpt, /over)
;------------------------------------------------------------
; Calculation of the curl
;------------------------------------------------------------
         zu = u*temporary(coefu)
         zv = v*temporary(coefv)
;
         psi =  (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
         psi = tabf*psi
;------------------------------------------------------------
; extraction of U and V on the appropriated domain
;------------------------------------------------------------
         if NOT keyword_set(key_periodic) OR nx NE jpi then begin
            psi[0, *, *] = !values.f_nan
            psi[nx-1, *, *] = !values.f_nan
         endif
         psi[*, 0, *] = !values.f_nan
         psi[*, ny-1, *] = !values.f_nan
         if n_elements(valmask) EQ 0 THEN valmask = 1e20
         terref =  where(tabf EQ 0)
         if terref[0] NE -1 then psi[temporary(terref)] = valmask
;----------------------------------------------------------------------------
         domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f']
         if keyword_set(direc) then psi = grossemoyenne(psi,direc,/nan)
;----------------------------------------------------------------------------
      END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyzt
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
      date1 NE date2 AND (size(u))[0] EQ 4:BEGIN 
         return, report('non code!')
      END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xy
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
      ELSE:BEGIN                ;xy
         indice2d = lindgen(jpi, jpj)
         indice2d = indice2d[indicex[0]:indicex[0]+nx-1,indicey[0]:indicey[0]+ny-1]
;------------------------------------------------------------
;------------------------------------------------------------
         case 1 of
            (size(u))[0] NE 2 OR (size(v))[0] NE 2: 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 
               if nxu NE nx then $
                if indicex[0] EQ firstxu then u = u[0:nx-1, *] ELSE u = u[1: nx, *] 
               IF nxv NE nx THEN $
                if indicex[0] EQ firstxv then v = v[0:nx-1, *] ELSE v = v[1: nx, *]
               IF nyu NE ny THEN $
                if indicey[0] EQ firstyu then u = u[*, 0:ny-1] ELSE u = u[*, 1: ny] 
               IF nyv NE ny THEN $
                if indicey[0] EQ firstyv then v = v[*, 0:ny-1] ELSE v = v[*, 1: ny]
            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
;------------------------------------------------------------
; Calculation of the curl
;------------------------------------------------------------
         coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
         terreu = where(coefu EQ 0)
         if terreu[0] NE -1 then coefu[temporary(terreu)] = !values.f_nan
         coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
         terrev = where(coefv EQ 0)
         if terrev[0] NE -1 then coefv[temporary(terrev)] = !values.f_nan
         tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d])
;
         zu = u*temporary(coefu)
         zv = v*temporary(coefv)

         psi =  (shift(zv, -1, 0)-zv) + (zu-shift(zu, 0, -1))
         psi = tabf*psi

;------------------------------------------------------------
; Edging put at !values.f_nan
;------------------------------------------------------------
         if  NOT keyword_set(key_periodic) OR nx NE jpi then begin
            psi[0, *] = !values.f_nan
            psi[nx-1, *] = !values.f_nan
         endif
         psi[*, 0] = !values.f_nan
         psi[*, ny-1] = !values.f_nan
;
         if n_elements(valmask) EQ 0 THEN valmask = 1e20
         terref =  where(tabf EQ 0)
         if terref[0] NE -1 then psi[temporary(terref)] = valmask
;------------------------------------------------------------
; for the graphic drawing
;------------------------------------------------------------
         domdef, (glamt[indice2d])[0, 0], (glamf[indice2d])[nx-1, 0],(gphit[indice2d])[0, 0], (gphif[indice2d])[0, ny-1], vert1, vert2, gridtype = ['t', 'f']
         if keyword_set(direc) then psi = moyenne(psi,direc,/nan)

      END
;----------------------------------------------------------------------------
   endcase
;------------------------------------------------------------
   if keyword_set(key_performance) THEN print, 'temps curl', systime(1)-tempsun 
;------------------------------------------------------------
   vargrid = 'F'
   varname = 'vorticity'
   return, psi
end