;+
;
; @file_comments
; north stereographic polar projection
;
; @categories
; Interpolation
;
; @param plam {in}{required}
; longitude position
;
; @param pphi {in}{required}
; latitude position
;
; @keyword DOUBLE {default=0}
; use double precision (default is float)
;
; @returns
; structure: {x:x, y:y} containing the point position in north stereographic polar projection
;
; @hidden
;
;-
;
FUNCTION fsnspp, plam, pphi, DOUBLE = double
;
  compile_opt idl2, strictarrsubs
;
  IF keyword_set(double) THEN BEGIN
    a = 2.d * tan( !dpi/4.d - !dpi/180.d*pphi/2.d )
    x = cos( !dpi/180.d*plam ) * a
    y = sin( !dpi/180.d*plam ) * a
  ENDIF ELSE BEGIN
    a = 2. * tan( !pi/4. - !pi/180.*float(pphi)/2. )
    x = cos( !pi/180.*float(plam) ) * a
    y = sin( !pi/180.*float(plam) ) * a
  ENDELSE
  RETURN, {x:x, y:y}
END
;+
; @file_comments Compute angles between grid lines and direction of the North pole
;(fom angle.F,v 2.2 in OPA8.2)
;
; @categories
; Interpolation
;
; @param fileocemesh {in}{required}{type=scalar string}
; a netcdf file that contains (at least) the following variables:
;        glamu, gphiu: longitudes and latitudes at U-points
;        glamv, gphiv: longitudes and latitudes at V-points
;        glamf, gphif: longitudes and latitudes at F-points
;
; @param gcosu {out}{type=2d array}
; cosinus of the angle between grid lines at U points and direction of the North pole
;
; @param gsinu {out}{type=2d array}
; sinus of the angle between grid lines at U points and direction of the North pole
;
; @param gcosv {out}{type=2d array}
; cosinus of the angle between grid lines at V points and direction of the North pole
;
; @param gsinv {out}{type=2d array}
; sinus of the angle between grid lines at V points and direction of the North pole
;
; @param gcost {out}{type=2d array}
; cosinus of the angle between grid lines at T points and direction of the North pole
;
; @param gsint {out}{type=2d array}
; sinus of the angle between grid lines at T points and direction of the North pole
;
; @keyword IODIRECTORY {type=scalar string}{default=''}
; the directory path where is located fileocemesh
;
; @keyword DOUBLE {type=1 ou 2}{default=0}
; put 1 to use double precision (default is float)
;
; @restrictions
; to compute the lateral boundary conditions, we assume that:
;     (1) the first line is similar to the second line
;       =>    gcosu[*, 0] = gcosu[*, 1]
;       =>    gsinu[*, 0] = gsinu[*, 1]
;     (2) the grid follows OPA x periodicity rule, first column is
;     equal to the next to last column
;       =>    gcosv[0, *] = gcosv[jpj-2, *]
;       =>    gsinv[0, *] = gsinv[jpj-2, *]
;
; @history
;       Original :  96-07 (O. Marti)
;                   98-06 (G. Madec)
;       Feb 2005: IDL adaptation S. Masson
;
; @version
; $Id$
;-
;
PRO angle, fileocemesh, gcosu, gsinu, gcosv, gsinv, gcost, gsint $
           , IODIRECTORY = iodirectory, DOUBLE = double
;
; 0. read oceanic grid parameters
; ================================
;
;
  compile_opt idl2, strictarrsubs
;
  IF keyword_set(IODIRECTORY) THEN BEGIN
    IF  strpos(iodirectory,'/',/reverse_search) NE (strlen(iodirectory)-1) THEN $
      iodirectory = iodirectory+'/'
  ENDIF ELSE iodirectory = ''
  fileoce = iodirectory+fileocemesh
;
  fileoce = findfile(fileoce, count = okfile)
  IF okfile NE 1 THEN BEGIN
    print, 'the file '+fileoce+' is not found... we stop'
    stop
  ENDIF
;
  cdfido = ncdf_open(fileoce[0])
  ncdf_varget, cdfido, 'glamt', glamt
  ncdf_varget, cdfido, 'glamu', glamu
  ncdf_varget, cdfido, 'glamv', glamv
  ncdf_varget, cdfido, 'glamf', glamf
  ncdf_varget, cdfido, 'gphit', gphit
  ncdf_varget, cdfido, 'gphiu', gphiu
  ncdf_varget, cdfido, 'gphiv', gphiv
  ncdf_varget, cdfido, 'gphif', gphif
  ncdf_close, cdfido
;
  glamt = reform(glamt, /over)
  glamu = reform(glamu, /over)
  glamv = reform(glamv, /over)
  glamf = reform(glamf, /over)
  gphit = reform(gphit, /over)
  gphiu = reform(gphiu, /over)
  gphiv = reform(gphiv, /over)
  gphif = reform(gphif, /over)
  jpj = (size(glamf, /dimension))[1]
;
; I. Compute the cosinus and sinus
; ================================
; (computation done on the north stereographic polar plan
;
;   ... north pole direction & modulous (at t-point)
  znpt = fsnspp( glamt, gphit, DOUBLE = double )
  glamt = -1 & gphit = -1; free memory
  znpt.x = - znpt.x
  znpt.y = - znpt.y
  znnpt = znpt.x*znpt.x + znpt.y*znpt.y
;   ... north pole direction & modulous (at u-point)
  znpu = fsnspp( glamu, gphiu, DOUBLE = double )
  glamu = -1 & gphiu = -1; free memory
  znpu.x = - znpu.x
  znpu.y = - znpu.y
  znnpu = znpu.x*znpu.x + znpu.y*znpu.y
;   ... north pole direction & modulous (at v-point)
  znpv = fsnspp( glamv, gphiv, DOUBLE = double )
  znpv00 = znpv
  znpv01 = fsnspp( shift(glamv, 0, 1), shift(gphiv, 0, 1), DOUBLE = double )
  glamv = -1 & gphiv = -1; free memory
  znpv.x = - znpv.x
  znpv.y = - znpv.y
  znnpv = znpv.x*znpv.x + znpv.y*znpv.y
;   ... f-point
  znpf00 = fsnspp( glamf, gphif, DOUBLE = double )
  znpf01 = fsnspp( shift(glamf, 0, 1), shift(gphif, 0, 1), DOUBLE = double )
  znpf10 = fsnspp( shift(glamf, 1, 0), shift(gphif, 1, 0), DOUBLE = double )
  glamf = -1 & gphif = -1; free memory
;   ... j-direction: v-point segment direction (t-point)
  zxvvt = znpv00.x - znpv01.x
  zyvvt = znpv00.y - znpv01.y
  zmnpvt = sqrt ( temporary(znnpt) * ( zxvvt*zxvvt + zyvvt*zyvvt )  )
  znpv00 = -1; free memory
  znpv01 = -1; free memory
  IF keyword_set(double) THEN zmnpvt = 1.e-14 > zmnpvt $
  ELSE zmnpvt = 1.e-6 > zmnpvt
;   ... j-direction: f-point segment direction (u-point)
  zxffu = znpf00.x - znpf01.x
  zyffu = znpf00.y - znpf01.y
  zmnpfu = sqrt ( temporary(znnpu) * ( zxffu*zxffu + zyffu*zyffu )  )
  znpf01 = -1; free memory
  IF keyword_set(double) THEN zmnpfu = 1.e-14 > zmnpfu $
  ELSE zmnpfu = 1.e-6 > zmnpfu
;   ... i-direction: f-point segment direction (v-point)
  zxffv = znpf00.x - znpf10.x
  zyffv = znpf00.y - znpf10.y
  znpf00 = -1 &  znpf10 = -1; free memory
  zmnpfv = sqrt ( temporary(znnpv) * ( zxffv*zxffv + zyffv*zyffv )  )
  IF keyword_set(double) THEN zmnpfv = 1.e-14 > zmnpfv $
  ELSE zmnpfv = 1.e-6 > zmnpfv
;   ... cosinus and sinus using scalar and vectorial products
  gsint = ( znpt.x*zyvvt - znpt.y*zxvvt ) / zmnpvt
  gcost = ( znpt.x*zxvvt + znpt.y*zyvvt ) / zmnpvt
;   ... cosinus and sinus using scalar and vectorial products
  gsinu = ( znpu.x*zyffu - znpu.y*zxffu ) / zmnpfu
  gcosu = ( znpu.x*zxffu + znpu.y*zyffu ) / zmnpfu
;   ... cosinus and sinus using scalar and vectorial products
;       (caution, rotation of 90 degres)
  gsinv =  ( znpv.x*zxffv + znpv.y*zyffv ) / zmnpfv
  gcosv = -( znpv.x*zyffv - znpv.y*zxffv ) / zmnpfv
;
; II. Geographic mesh
; ===================
;
;       bad = where(abs(glamf-shift(glamf, 0, 1)) LT 1.e-8)
;       IF bad[0] NE -1 THEN BEGIN
;         gcosu[bad] = 1.
;         gsinu[bad] = 0.
;       ENDIF
;       bad = where(abs(gphif-shift(gphif, 1, 0)) LT 1.e-8)
;       IF bad[0] NE -1 THEN BEGIN
;         gcosv[bad] = 1.
;         gsinv[bad] = 0.
;       ENDIF
;
; III. Lateral boundary conditions
; ================================
;
  gcost[*, 0] = gcost[*, 1]
  gsint[*, 0] = gsint[*, 1]
  gcosu[*, 0] = gcosu[*, 1]
  gsinu[*, 0] = gsinu[*, 1]
  gcosv[0, *] = gcosv[jpj-2, *]
  gsinv[0, *] = gsinv[jpj-2, *]
;
  RETURN
END