PRO angle, glamu, gphiu, glamv, gphiv, glamf, gphif, gsinu, gsinv, gcosu, gcosv
;;;---------------------------------------------------------------------
;;;
;;;                       ROUTINE angle
;;;                     *******************
;;;
;;;  Purpose :
;;;  --------
;;;     Compute angles between grid lines and direction of the North
;;;
;;   Method :
;;   -------
;;
;;   Input :
;;   ------
;;      common
;;            /comcoh/ glamf,   : longitudes and latitudes at F-points
;;                     gphif
;;
;;   Output :
;;   -------
;;      common
;;            /comcoh/ gsinu,gcosu : sinus and cosinus of the angle
;;                     gsinv,gcosv   between north-south direction
;;                                   and the j-direction of the mesh
;;
;;   Modifications:
;;   --------------
;;       Original :  96-07 (O. Marti)
;;                   98-06 (G. Madec)
;;                   99-11 (A. Jouzeau) : Adaptation to IDL
;;----------------------------------------------------------------------
;;----------------------------------------------------------------------
;
;;;---------------------------------------------------------------------
;;;  OPA8, LODYC (1997)
;;;---------------------------------------------------------------------
;
; I. Compute the cosinus and sinus
; ================================
;
;     ... north pole direction & modulous (at u-point)
      zxnpu = 0. - fsxspp( glamu, gphiu )
      zynpu = 0. - fsyspp( glamu, gphiu )
      znnpu = zxnpu*zxnpu + zynpu*zynpu
;
;     ... north pole direction & modulous (at v-point)
      zxnpv = 0. - fsxspp( glamv, gphiv )
      zynpv = 0. - fsyspp( glamv, gphiv )
      znnpv = zxnpv*zxnpv + zynpv*zynpv
;
;     ... j-direction: f-point segment direction (u-point)
      zxffu=  fsxspp( glamf, gphif ) - fsxspp( shift(glamf, 0, 1), shift(gphif, 0, 1) )
      zyffu=  fsyspp( glamf, gphif ) - fsyspp( shift(glamf, 0, 1), shift(gphif, 0, 1) )
      zmnpfu= sqrt ( znnpu * ( zxffu*zxffu + zyffu*zyffu )  )


;
;     ... i-direction: f-point segment direction (v-point)
      zxffv=  fsxspp( glamf, gphif ) - fsxspp( shift(glamf, 1, 0), shift(gphif, 1, 0) )
      zyffv=  fsyspp( glamf, gphif ) - fsyspp( shift(glamf, 1, 0), shift(gphif, 1, 0) )
      zmnpfv= sqrt ( znnpv * ( zxffv*zxffv + zyffv*zyffv )  )

;
;     ... cosinus and sinus using scalar and vectorial products
      gsinu = ( zxnpu*zyffu - zynpu*zxffu ) / zmnpfu
      gcosu = ( zxnpu*zxffu + zynpu*zyffu ) / zmnpfu


;
;     ... cosinus and sinus using scalar and vectorial products
;     (caution, rotation of 90 degres)
      gsinv = ( zxnpv*zxffv + zynpv*zyffv ) / zmnpfv
      gcosv =-( zxnpv*zyffv - zynpv*zxffv ) / zmnpfv


;
;
;     II. Geographic mesh
;     ===================
;
;
      ind = where(abs(glamf-shift(glamf, 0, 1)) LT 1.e-8)
      gsinu(ind) = 0.d
      gcosu(ind) = 1.d
;
      ind = where(abs(gphif-shift(gphif, 1, 0)) LT 1.e-8)
      gsinv(ind) = 0.d
      gcosv(ind) = 1.d
;
return
END