source: trunk/SRC/Computation/curl.pro @ 248

;+
;
; @file_comments
; Calculate the vertical component of the curl of a field of horizontal vectors
;
; @categories
; Calculation
;
; @param UU
; Matrix representing the zonal coordinates (U point) of a field of vectors
; A 2D (xy), 3D (xyz or yt) or a structure readable by litchamp and containing
; a 2D (xy), 3D (xyz or yt) array (4D case is not coded yet).
; Note that the dimension of the array must suit the domain dimension.
;
; @param VV
; Matrix representing the meridional coordinates (V point) of a field of vectors
; A 2D (xy), 3D (xyz or yt) or a structure readable by litchamp and containing
; a 2D (xy), 3D (xyz or yt) array (4D case is not coded yet).
; Note that the dimension of the array must suit the domain dimension.
;
; @keyword DIREC {type=scalar string}
; Use if you want to call <pro>moyenne</pro> or
; <pro>grossemoyenne</pro> after the div computation
; with a mean done in the DIREC direction
;
; @returns
; the vertical component of the curl of the input data (with the same size)
;
; @uses
; cm_4cal
; cm_4data
; cm_4mmesh
;
; @restrictions
;
; - Works only for Arakawa C-grid.
; - UU must be on U grid, VV must be on V grid
; - 4D case is not coded yet
; - the common variable jpt is used to differ xyz (jpt=1) and xyt (jpt\=1) cases.
; - U and V arrays are cut in the same geographic domain. Because of the shift between
;   T, U, V and F grids, it is possible that these two arrays do not have the same
;   size and refer to different indexes. In this case, arrays are re-cut on
;   common indexes. To avoid these re-cuts, use the keyword /memeindice in
; <pro>domdef</pro>
; - When computing the divergence, we update, vargrid, varname, varunits and the
;   grid position parameters (firstxf, lastxf, nxf, firstyf, lastyf, nyf).
; - points that cannot be computed (domain boundaries, coastline) are set to NaN
;
; @examples
; IDL> \@tst_initorca2
; IDL> plt, curl(dist(jpi,jpj), dist(jpi,jpj))
;
; @history
; Guillaume Roullet (grlod\@ipsl.jussieu.fr)
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
; adaptation to work with a reduce domain
; 21/5/1999: missing values at !values.f_nan
;
; @version
; $Id$
;
; @todo
; code the 4D case
;
;-
;
FUNCTION curl, uu, vv, DIREC = DIREC
;
  compile_opt idl2, strictarrsubs
;
@cm_4cal                        ; for jpt
@cm_4data                       ; for varname, vargrid, vardate, varunit, valmask
@cm_4mesh
;
  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)

  szu = size(u)
  szv = size(v)

  if szu[0] NE szv[0] then return, report('U and V input data must have the same number of dimensions')

;------------------------------------------------------------
; 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)
  indice2d = lindgen(jpi, jpj)
  indice2d = indice2d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1]
;----------------------------------------------------------------------------
  vargrid = 'F'
  varname = 'vorticity'
  varunits = 's-1'
  if n_elements(valmask) EQ 0 THEN valmask = 1e20
  firstxt = indicex[0] & lastxt = indicex[0]+nx-1 & nxt = nx
  firstyt = indicey[0] & lastyt = indicey[0]+ny-1 & nyt = ny
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
  case 1 of
;----------------------------------------------------------------------------
;xyz
;----------------------------------------------------------------------------
    szu[0] EQ 3 AND jpt EQ 1:BEGIN
;------------------------------------------------------------
; extraction of U and V on the appropriated domain
;------------------------------------------------------------
      case 1 of
        szu[1] EQ nxu AND szu[2] EQ nyu AND $
           szv[1] EQ nxv AND szv[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
        szu[1] EQ jpi AND szu[2] EQ jpj AND $
           szv[1] EQ jpi AND szv[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
;------------------------------------------------------------
; curl computation
;------------------------------------------------------------
      coefu = ((e1u[indice2d])[*]#replicate(1., nzt)) $
              * (umask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
      landu = where(coefu EQ 0)
      if landu[0] NE -1 then coefu[temporary(landu)] = !values.f_nan

      coefv = ((e2v[indice2d])[*]#replicate(1., nzt)) $
              *(vmask())[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
      landv = where(coefv EQ 0)
      if landv[0] NE -1 then coefv[temporary(landv)] = !values.f_nan

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

      psi = (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
      psi = temporary(tabf) * temporary(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 landf[0] NE -1 then psi[temporary(landf)] = valmask
      if keyword_set(direc) then psi = moyenne(psi, direc, /nan)
    END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyt
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
    szu[0] EQ 3 AND jpt GT 1:BEGIN
;------------------------------------------------------------
; extraction of U and V on the appropriated domain
;------------------------------------------------------------
      case 1 of
        szu[1] EQ nxu AND szu[2] EQ nyu AND $
           szv[1] EQ nxv AND szv[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
        szu[1] EQ jpi AND szu[2] EQ jpj AND $
           szv[1] EQ jpi AND szv[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
;----------------------------------------------------------------------------
; curl computation
;----------------------------------------------------------------------------
      coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
      landu = where(coefu EQ 0)
      if landu[0] NE -1 then coefu[temporary(landu)] = !values.f_nan
      coefu = temporary(coefu[*])#replicate(1., jpt)
;
      coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
      landv = where(coefv EQ 0)
      if landv[0] NE -1 then coefv[temporary(landv)] = !values.f_nan
      coefv = temporary(coefv[*])#replicate(1., jpt)
;
      tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d])
      tabf = reform(temporary(tabf[*])#replicate(1., jpt), nx, ny, jpt, /overwrite)
      landf = where(tabf EQ 0)
;
      zu = temporary(u) * temporary(coefu)
      zv = temporary(v) * temporary(coefv)
;
      psi = (shift(zv, -1, 0, 0)-zv) + (zu-shift(zu, 0, -1, 0))
      psi = temporary(tabf) * temporary(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 landf[0] NE -1 then psi[temporary(landf)] = valmask
      if keyword_set(direc) then psi = grossemoyenne(psi, direc, /nan)
    END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyzt
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
    szu[0] EQ 4:BEGIN
      return, report('Case not coded contact saxo team or make a do loop!')
    END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xy
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
    szu[0] EQ 2:BEGIN
;------------------------------------------------------------
;------------------------------------------------------------
      case 1 of
        szu[1] EQ nxu AND szu[2] EQ nyu AND $
           szv[1] EQ nxv AND szv[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
        szu[1] EQ jpi AND szu[2] EQ jpj AND $
           szv[1] EQ jpi AND szv[2] EQ jpj:BEGIN
          u = u[indice2d]
          v = v[indice2d]
        END
        ELSE:return,  -1
      endcase
;------------------------------------------------------------
; curl computation
;------------------------------------------------------------
      coefu = e1u[indice2d]*(umask())[indice2d+jpi*jpj*firstzt]
      landu = where(coefu EQ 0)
      if landu[0] NE -1 then coefu[temporary(landu)] = !values.f_nan
      coefv = e2v[indice2d]*(vmask())[indice2d+jpi*jpj*firstzt]
      landv = where(coefv EQ 0)
      if landv[0] NE -1 then coefv[temporary(landv)] = !values.f_nan
      tabf = (fmask())[indice2d+jpi*jpj*firstzt]/(e1f[indice2d]*e2f[indice2d])
      landf =  where(tabf EQ 0)
;
      zu = temporary(u) * temporary(coefu)
      zv = temporary(v) * temporary(coefv)

      psi = (shift(zv, -1, 0)-zv) + (zu-shift(zu, 0, -1))
      psi = temporary(tabf) * temporary(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 landf[0] NE -1 then psi[temporary(landf)] = valmask
      if keyword_set(direc) then psi = moyenne(psi, direc, /nan)
    END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
    ELSE:return, report('U and V input arrays must have 2, 3 or 4 dimensions')
  ENDCASE
;------------------------------------------------------------
  if keyword_set(key_performance) THEN print, 'temps curl', systime(1)-tempsun

  return, psi
end