;+
;
; @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:
; glamt, gphit: longitudes and latitudes at T-points
; 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
;
; @param gcosf {out}{type=2d array}
; cosinus of the angle between grid lines at F points and direction of the North pole
;
; @param gsinf {out}{type=2d array}
; sinus of the angle between grid lines at F points and direction of the North pole
;
; @keyword IODIRECTORY {type=scalar string}{default=''}
; the directory path where is located fileocemesh
;
; @keyword DOUBLE {type=1 or 0}{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 and north pole periodicity rule (see lbcorca)
;
; @history
; Original : 96-07 (O. Marti)
; 98-06 (G. Madec)
; Feb 2005: IDL adaptation S. Masson
; May 2007: F points + call to lbcorca
;
; @version
; $Id$
;
;-
;
PRO angle, fileocemesh, gcosu, gsinu, gcosv, gsinv, gcost, gsint, gcosf, gsinf $
, IODIRECTORY = iodirectory, DOUBLE = double
;
compile_opt idl2, strictarrsubs
;
; 0. read oceanic grid parameters
; ================================
;
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
ras = report('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)
;
sz = size(glamf, /dimension)
jpi = sz[0]
jpj = sz[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 )
znpu00 = znpu
znpu0i = fsnspp( shift(glamu, 0, -1), shift(gphiu, 0, -1), 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
; ... north pole direction & modulous (at f-point)
znpf = fsnspp( glamf, gphif, DOUBLE = double )
znpf00 = znpf
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
znpf.x = - znpf.x
znpf.y = - znpf.y
znnpf = znpf.x*znpf.x + znpf.y*znpf.y
;
; ... 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 & 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
; ... j-direction: u-point segment direction (f-point)
zxuuf = znpu0i.x - znpu00.x
zyuuf = znpu0i.y - znpu00.y
zmnpuf = sqrt ( temporary(znnpf) * ( zxuuf*zxuuf + zyuuf*zyuuf ) )
znpu00 = -1 & znpu0i = -1; free memory
IF keyword_set(double) THEN zmnpuf = 1.e-14 > zmnpuf $
ELSE zmnpuf = 1.e-6 > zmnpuf
; ... 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
; ... cosinus and sinus using scalar and vectorial products
gsinf = ( znpf.x*zyuuf - znpf.y*zxuuf ) / zmnpuf
gcosf = ( znpf.x*zxuuf + znpf.y*zyuuf ) / zmnpuf
;
; 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]
;
IF keyword_set(double) THEN sgn = 1.d ELSE sgn = 1.
dummy = lbcorca(gcost, 'T', sgn, /correction)
dummy = lbcorca(gsint, 'T', -sgn, /correction)
dummy = lbcorca(gcosu, 'U', sgn, /correction)
dummy = lbcorca(gsinu, 'U', -sgn, /correction)
dummy = lbcorca(gcosv, 'V', sgn, /correction)
dummy = lbcorca(gsinv, 'V', -sgn, /correction)
dummy = lbcorca(gcosf, 'F', sgn, /correction)
dummy = lbcorca(gsinf, 'F', -sgn, /correction)
;
RETURN
END