;---------
;+
; @file_comments north stereographic polar projection
;
; @keyword    /DOUBLE use double precision (default is float)
;
; @returns
;       gsinu,gcosu : sinus and cosinus of the angle 
;       gsinv,gcosv   between north-south direction 
;       gsint,gcost   and the j-direction of the mesh
;
; @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 
;-
;---------
;
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
;(fom angle.F,v 2.2 in OPA8.2)
;
; @param fileocemesh {in}{required} a netcdf file that contains (at least):
;        glamu, gphiu: longitudes and latitudes at U-points
;        glamv, gphiv: longitudes and latitudes at V-points
;        glamf, gphif: longitudes and latitudes at F-points
;
; @keyword IODIRECTORY the directory path where is located fileocemesh
; @keyword    /DOUBLE use double precision (default is float)
;-
;---------
;
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