source: trunk/SRC/Computation/norm.pro @ 325

;+
;
; @file_comments
; calculate the norm of vectors field located on Arakawa C-grid
;
; @categories
; Calculation
;
; @param UU {in}{required}
; Matrix representing the zonal coordinates (at U/V point) of a field of vectors
; A 2D (xy), 3D (xyz or yt), 4D (xyzt) or a structure readable by
; <pro>litchamp</pro> and containing a 2D (xy), 3D (xyz or yt), 4D (xyzt) array.
; Note that the dimension of the array must suit the domain dimension.
;
; @param VV {in}{required}
; Matrix representing the meridional coordinates (at V/U point) of a field of vectors
; A 2D (xy), 3D (xyz or yt), 4D (xyzt) or a structure readable by
; <pro>litchamp</pro> and containing a 2D (xy), 3D (xyz or yt), 4D (xyzt) array.
; Note that the dimension of the array must suit the domain dimension.
;
; @keyword DIREC
; 't' 'x' 'y' 'z' 'xys' 'xz' 'yz' 'xyz' 'xt' 'yt' 'zt' 'xyt'
;       'xzt' 'yzt' 'xyzt' Direction on which do averages
;
; @keyword _EXTRA
; Used to declare that this routine accepts inherited keywords
;
; @returns
; A 2D (xy), 3D (xyz or yt), 4D (xyzt) Array
;
; @uses
; cm_4mesh
; cm_4data
; cm_4cal
;
; @restrictions
; The norm is calculated on points T. To do this calculation, we average
; field U and V on points T before calculate the norm. At the edge of
; coast and of domain, we can not calculate fields U and V at points T,
; that is why these points are at value !values.f_nan.
;
; When we calculate on a reduce geographic domain, field U and V have not
; necessarily the same number of point. In this case, we recut U and V to
; keep only common points. We profit of this to redo a <pro>domdef</pro> which redefine
; a geographic domain on which fields U and V are extracted on same points
;
; To know what type of array we work with, we  test its size and dates
; gave by time[0] and time[jpt-1] to know if thee is a temporal dimension.
; Before to start norm, make sure that time and jpt are defined how
; they have to!
;
; @examples
; To calculate the average of the norm of streams on all the domain
; between 0 and 50:
;      IDL> domdef, 0, 50
;      IDL> res = norm(un, vn, dir = 'xyz')
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                       9/6/1999
; @version
; $Id$
;
;-
FUNCTION norm, uu, vv, DIREC = direc, _EXTRA = ex
;
  compile_opt idl2, strictarrsubs
;
@cm_4mesh
@cm_4data
@cm_4cal
  IF NOT keyword_set(key_forgetold) THEN BEGIN
@updatenew
@updatekwd
  ENDIF
;---------------------------------------------------------
  time1 = 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 norm is based on Arakawa C-grid.' $
                     , 'U and V grids must therefore be defined'])
;------------------------------------------------------------
  if NOT keyword_set(direc) then direc = 0
; construction of u and v at points T
  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
;
  grilleu = litchamp(uu, /grid)
  if grilleu EQ '' then grilleu = 'U'
  grillev = litchamp(vv, /grid)
  if grillev EQ '' then grillev = 'V'
  IF grilleu EQ 'V' AND grillev EQ 'U' THEN inverse = 1
  IF grilleu EQ 'T' AND grillev EQ 'T' THEN BEGIN
    interpolle  = 0
    return, report('Case not coded, but easy to do...!')
  ENDIF ELSE interpolle = 1
  if keyword_set(inverse) then begin
    tmp = u
    u = temporary(v)
    v = temporary(tmp)
  endif
;------------------------------------------------------------
; 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)
;----------------------------------------------------------------------------
  vargrid = 'T'
  varname = 'norm'
  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
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
    (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]
      indice3d = lindgen(jpi, jpj, jpk)
      indice3d = indice3d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
;------------------------------------------------------------
; 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
          case (size(u))[3] OF
            nzt: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
            jpk:BEGIN
              if nxu NE nx then $
                 if indicex[0] EQ firstxu then u = u[0:nx-1, *, firstzt:lastzt] ELSE u = u[1: nx, *, firstzt:lastzt]
              IF nxv NE nx THEN $
                 if indicex[0] EQ firstxv then v = v[0:nx-1, *, firstzt:lastzt] ELSE v = v[1: nx, *, firstzt:lastzt]
              IF nyu NE ny THEN $
                 if indicey[0] EQ firstyu then u = u[*, 0:ny-1, firstzt:lastzt] ELSE u = u[*, 1: ny, firstzt:lastzt]
              IF nyv NE ny THEN $
                 if indicey[0] EQ firstyv then v = v[*, 0:ny-1, firstzt:lastzt] ELSE v = v[*, 1: ny, firstzt:lastzt]
            end
            ELSE: return, report(['the third dimension of u (' + strtrim((size(u))[3], 1) $
                                  +') must be equal to nzt (' + strtrim(nzt, 1) + ') or jpk ('+strtrim(jpk, 1)+')'])
          endcase
        END
        (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $
           (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN
          u = u[indice3d]
          v = v[indice3d]
        END
        ELSE: return $
           , report(['We cannot find re-cut on common indexes (see restriction un the header.).' $
                     , 'To avoid this problem, read the full domain' $
                     , 'or use the keyword /memeindice in domdef when defining the zoom area'])
      endcase
;------------------------------------------------------------------
; We reshape u and v to make sure that no dimension has been erased
;------------------------------------------------------------------
      if nzt EQ 1 then begin
        u = reform(u, nx, ny, nzt, /over)
        v = reform(v, nx, ny, nzt, /over)
      endif
;------------------------------------------------------------------
; construction of u and v at points T
;-----------------------------------------------------------
      a = u[0, *, *]
      u = (u+shift(u, 1, 0, 0))/2.
      if NOT keyword_set(key_periodic) OR nx NE jpi then u[0, *, *] = a
      a = v[*, 0, *]
      v = (v+shift(v, 0, 1, 0))/2.
      if NOT keyword_set(key_periodic) OR nx NE jpi then v[*, 0, *] = a
;----------------------------------------------------------------------------
; attribution of the mask and of longitude and latitude arrays
;----------------------------------------------------------------------------
      mask = tmask[indice3d]
      if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over)
;-----------------------------------------------------------
      if n_elements(valmask) EQ 0 THEN valmask = 1e20
      landu = where(u GE valmask/10)
      if landu[0] NE -1 then u[landu] = 0
      landv = where(v GE valmask/10)
      if landv[0] NE -1 then v[landv] = 0
      res = sqrt(u^2+v^2)
      if NOT keyword_set(key_periodic) OR nx NE jpi then res[0, *, *] = !values.f_nan
      res[*, 0, *] = !values.f_nan
      mask = where(mask eq 0)
      IF mask[0] NE -1 THEN res[mask] = valmask
; All kind of average
      domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0], (gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
      if keyword_set(direc) then res = moyenne(res, direc, /nan, boxzoom = boxzoom, /nodomdef)
;-----------------------------------------------------------
    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:return $
           , report(['We cannot find re-cut on common indexes (see restriction un the header.).' $
                     , 'To avoid this problem, read the full domain' $
                     , 'or use the keyword /memeindice in domdef when defining the zoom area'])
      endcase
;------------------------------------------------------------------
; construction of u and v at points T
;-----------------------------------------------------------
      a = u[0, *, *]
      u = (u+shift(u, 1, 0, 0))/2.
      if NOT keyword_set(key_periodic) OR nx NE jpi then u[0, *, *] = a
      a = v[*, 0, *]
      v = (v+shift(v, 0, 1, 0))/2.
      if NOT keyword_set(key_periodic) OR nx NE jpi then v[*, 0, *] = a
;----------------------------------------------------------------------------
; attribution of the mask and of longitude and latitude arrays.
; We recover the complete grid to establish a big mask extent in the four
; direction to cover pointsfor which a land point has been
; considerated (make a small drawing)
;----------------------------------------------------------------------------
      mask = tmask[indice2d+jpi*jpj*firstzt]
      if ny EQ 1 then mask = reform(mask, nx, ny, /over)
;-----------------------------------------------------------
; construction of land containing all points to mask
;-----------------------------------------------------------
      if n_elements(valmask) EQ 0 THEN valmask = 1e20
      landu = where(u GE valmask/10)
      if landu[0] NE -1 then u[landu] = 0
      landv = where(v GE valmask/10)
      if landv[0] NE -1 then v[landv] = 0
      res = sqrt(u^2+v^2)
      if NOT keyword_set(key_periodic) OR nx NE jpi then res[0, *, *] = !values.f_nan
      res[*, 0, *] = !values.f_nan
      mask = where(mask eq 0)
      IF mask[0] NE -1 THEN BEGIN
        coeftps = lindgen(jpt)*nx*ny
        coeftps = replicate(1, n_elements(mask))#coeftps
        mask = (temporary(mask))[*]#replicate(1, jpt)
        mask = temporary(mask[*]) + temporary(coeftps[*])
        res[temporary(mask)] = valmask
      ENDIF
; moyennes en tous genres
      domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0], (gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
      if keyword_set(direc) then res = grossemoyenne(res, direc, /nan, boxzoom = boxzoom, /nodomdef)
    END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xyzt
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
    date1 NE date2 AND (size(u))[0] EQ 4:BEGIN
      indice2d = lindgen(jpi, jpj)
      indice2d = indice2d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1]
      indice3d = lindgen(jpi, jpj, jpk)
      indice3d = indice3d[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt]
;------------------------------------------------------------
; 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
          case (size(u))[3] OF
            nzt: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
            jpk:BEGIN
              if nxu NE nx then $
                 if indicex[0] EQ firstxu then u = u[0:nx-1, *, firstzt:lastzt, *] ELSE u = u[1: nx, *, firstzt:lastzt, *]
              IF nxv NE nx THEN $
                 if indicex[0] EQ firstxv then v = v[0:nx-1, *, firstzt:lastzt, *] ELSE v = v[1: nx, *, firstzt:lastzt, *]
              IF nyu NE ny THEN $
                 if indicey[0] EQ firstyu then u = u[*, 0:ny-1, firstzt:lastzt, *] ELSE u = u[*, 1: ny, firstzt:lastzt, *]
              IF nyv NE ny THEN $
                 if indicey[0] EQ firstyv then v = v[*, 0:ny-1, firstzt:lastzt, *] ELSE v = v[*, 1: ny, firstzt:lastzt, *]
            end
            ELSE: report, (['the third dimension of u (' + strtrim((size(u))[3], 1) $
                            +') must be equal to nzt (' + strtrim(nzt, 1) + ') or jpk ('+strtrim(jpk, 1)+')'])
          endcase
        END
        (size(u))[1] EQ jpi AND (size(u))[2] EQ jpj AND (size(u))[3] EQ jpk AND $
           (size(v))[1] EQ jpi AND (size(v))[2] EQ jpj AND (size(u))[3] EQ jpk :BEGIN
          u = u[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt, *]
          v = v[indicex[0]:indicex[0]+nx-1, indicey[0]:indicey[0]+ny-1, firstzt:lastzt, *]
        END
        ELSE: return $
           , report(['We cannot find re-cut on common indexes (see restriction un the header.).' $
                     , 'To avoid this problem, read the full domain' $
                     , 'or use the keyword /memeindice in domdef when defining the zoom area'])
      endcase
;------------------------------------------------------------------
; construction of u and v at points T
;-----------------------------------------------------------
      a = u[0, *, *, *]
      u = (u+shift(u, 1, 0, 0, 0))/2.
      if NOT keyword_set(key_periodic) OR nx NE jpi then u[0, *, *, *] = a
      a = v[*, 0, *, *]
      v = (v+shift(v, 0, 1, 0, 0))/2.
      if NOT keyword_set(key_periodic) OR nx NE jpi then v[*, 0, *, *] = a
;----------------------------------------------------------------------------
; attribution of the mask and of longitude and latitude arrays
;----------------------------------------------------------------------------
      mask = tmask[indice3d]
      if nzt EQ 1 then mask = reform(mask, nx, ny, nzt, /over)
;-----------------------------------------------------------
      if n_elements(valmask) EQ 0 THEN valmask = 1e20
      landu = where(u GE valmask/10)
      if landu[0] NE -1 then u[landu] = 0
      landv = where(v GE valmask/10)
      if landv[0] NE -1 then v[landv] = 0
      res = sqrt(u^2+v^2)
      if NOT keyword_set(key_periodic) OR nx NE jpi then res[0, *, *, *] = !values.f_nan
      res[*, 0, *, *] = !values.f_nan
      mask = where(mask eq 0)
      IF mask[0] NE -1 THEN BEGIN
        coeftps = lindgen(jpt)*nx*ny*nzt
        coeftps = replicate(1, n_elements(mask))#coeftps
        mask = (temporary(mask))[*]#replicate(1, jpt)
        mask = temporary(mask[*]) + temporary(coeftps[*])
        res[temporary(mask)] = valmask
      ENDIF
; All kind of average
      domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0], (gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
      if keyword_set(direc) then res = grossemoyenne(res, direc, /nan, boxzoom = boxzoom, /nodomdef)
    END
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
;xy
;----------------------------------------------------------------------------
;----------------------------------------------------------------------------
    ELSE:BEGIN                  ;xy
      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[indice2d]
          v = v[indice2d]
        END
        ELSE:return $
           , report(['We cannot find re-cut on common indexes (see restriction un the header.).' $
                     , 'To avoid this problem, read the full domain' $
                     , 'or use the keyword /memeindice in domdef when defining the zoom area'])
      endcase
;------------------------------------------------------------------
; We reshape u and v to make sure that no dimension has been erased
;------------------------------------------------------------------
      if ny EQ 1 then begin
        u = reform(u, nx, ny, /over)
        v = reform(v, nx, ny, /over)
      endif
;------------------------------------------------------------------
; construction of u and v at points T
;-----------------------------------------------------------
      a = u[0, *]
      u = (u+shift(u, 1, 0))/2.
      if NOT keyword_set(key_periodic) OR nx NE jpi then u[0, *] = a
      a = v[*, 0]
      v = (v+shift(v, 0, 1))/2.
      if NOT keyword_set(key_periodic) OR nx NE jpi then v[*, 0] = a
;----------------------------------------------------------------------------
; attribution of the mask and of longitude and latitude arrays.
; We recover the complete grid to establish a big mask extent in the four
; direction to cover pointsfor which a land point has been
; considerated (make a small drawing)
;----------------------------------------------------------------------------
      mask = tmask[indice2d+jpi*jpj*firstzt]
      if nyt EQ 1 THEN mask = reform(mask, nx, ny, /over)
;-----------------------------------------------------------
; construction of land containing all points to mask
;-----------------------------------------------------------
      if n_elements(valmask) EQ 0 THEN valmask = 1e20
      landu = where(u GE valmask/10)
      if landu[0] NE -1 then u[landu] = 0
      landv = where(v GE valmask/10)
      if landv[0] NE -1 then v[landv] = 0
      res = sqrt(u^2+v^2)
      if NOT keyword_set(key_periodic) OR nx NE jpi then res[0, *] = !values.f_nan
      res[*, 0] = !values.f_nan
      mask = where(mask eq 0)
      IF mask[0] NE -1 THEN res[mask] = valmask
; All kind of average
      domdef, (glamt[indice2d])[0, 0], (glamu[indice2d])[nx-1, 0], (gphit[indice2d])[0, 0], (gphiv[indice2d])[0, ny-1], vert1, vert2, /meme
      if keyword_set(direc) then res = moyenne(res, direc, /nan, boxzoom = boxzoom, /nodomdef)
    END
;----------------------------------------------------------------------------
  endcase
;------------------------------------------------------------
  if keyword_set(key_performance) THEN print, 'time norm', systime(1)-time1
  return, res
end