source: trunk/SRC/Interpolation/angle.pro @ 114

Last change on this file since 114 was 114, checked in by smasson, 18 years ago

new compilation options (compile_opt idl2, strictarrsubs) in each routine

  • Property svn:executable set to *
File size: 6.2 KB
Line 
1;---------
2;+
3; @file_comments north stereographic polar projection
4;
5; @keyword    /DOUBLE use double precision (default is float)
6;
7; @returns
8;       gsinu,gcosu : sinus and cosinus of the angle
9;       gsinv,gcosv   between north-south direction
10;       gsint,gcost   and the j-direction of the mesh
11;
12; @restrictions to compute the lateral boundary conditions, we assume
13; that:
14;     (1) the first line is similar to the second line
15;       =>    gcosu[*, 0] = gcosu[*, 1]
16;       =>    gsinu[*, 0] = gsinu[*, 1]
17;     (2) the grid follows OPA x periodicity rule, first column is
18;     equal to the next to last column
19;       =>    gcosv[0, *] = gcosv[jpj-2, *]
20;       =>    gsinv[0, *] = gsinv[jpj-2, *]
21;
22;
23; @history
24;   --------------
25;       Original :  96-07 (O. Marti)
26;                   98-06 (G. Madec)
27;       Feb 2005: IDL adaptation S. Masson
28;-
29;---------
30;
31FUNCTION fsnspp, plam, pphi, DOUBLE = double
32;
33  compile_opt idl2, strictarrsubs
34;
35  IF keyword_set(double) THEN BEGIN
36    a = 2.d * tan( !dpi/4.d - !dpi/180.d*pphi/2.d )
37    x = cos( !dpi/180.d*plam ) * a
38    y = sin( !dpi/180.d*plam ) * a
39  ENDIF ELSE BEGIN
40    a = 2. * tan( !pi/4. - !pi/180.*float(pphi)/2. )
41    x = cos( !pi/180.*float(plam) ) * a
42    y = sin( !pi/180.*float(plam) ) * a   
43  ENDELSE
44  RETURN, {x:x, y:y}
45END
46;---------
47;+
48; @file_comments Compute angles between grid lines and direction of the North
49;(fom angle.F,v 2.2 in OPA8.2)
50;
51; @param fileocemesh {in}{required} a netcdf file that contains (at least):
52;        glamu, gphiu: longitudes and latitudes at U-points
53;        glamv, gphiv: longitudes and latitudes at V-points
54;        glamf, gphif: longitudes and latitudes at F-points
55;
56; @keyword IODIRECTORY the directory path where is located fileocemesh
57; @keyword    /DOUBLE use double precision (default is float)
58;-
59;---------
60;
61PRO angle, fileocemesh, gcosu, gsinu, gcosv, gsinv, gcost, gsint $
62           , IODIRECTORY = iodirectory, DOUBLE = double
63;
64; 0. read oceanic grid parameters
65; ================================
66;
67;
68  compile_opt idl2, strictarrsubs
69;
70  IF keyword_set(IODIRECTORY) THEN BEGIN
71    IF  strpos(iodirectory,'/',/reverse_search) NE (strlen(iodirectory)-1) THEN $
72      iodirectory = iodirectory+'/'
73  ENDIF ELSE iodirectory = ''
74  fileoce = iodirectory+fileocemesh
75;
76  fileoce = findfile(fileoce, count = okfile)
77  IF okfile NE 1 THEN BEGIN
78    print, 'the file '+fileoce+' is not found... we stop'
79    stop
80  ENDIF
81;
82  cdfido = ncdf_open(fileoce[0])
83  ncdf_varget, cdfido, 'glamt', glamt
84  ncdf_varget, cdfido, 'glamu', glamu
85  ncdf_varget, cdfido, 'glamv', glamv
86  ncdf_varget, cdfido, 'glamf', glamf
87  ncdf_varget, cdfido, 'gphit', gphit
88  ncdf_varget, cdfido, 'gphiu', gphiu
89  ncdf_varget, cdfido, 'gphiv', gphiv
90  ncdf_varget, cdfido, 'gphif', gphif
91  ncdf_close, cdfido
92;
93  glamt = reform(glamt, /over)
94  glamu = reform(glamu, /over)
95  glamv = reform(glamv, /over)
96  glamf = reform(glamf, /over)
97  gphit = reform(gphit, /over)
98  gphiu = reform(gphiu, /over)
99  gphiv = reform(gphiv, /over)
100  gphif = reform(gphif, /over)
101  jpj = (size(glamf, /dimension))[1]
102;
103; I. Compute the cosinus and sinus
104; ================================
105; (computation done on the north stereographic polar plan
106;
107;   ... north pole direction & modulous (at t-point)
108  znpt = fsnspp( glamt, gphit, DOUBLE = double )
109  glamt = -1 & gphit = -1; free memory
110  znpt.x = - znpt.x
111  znpt.y = - znpt.y
112  znnpt = znpt.x*znpt.x + znpt.y*znpt.y
113;   ... north pole direction & modulous (at u-point)
114  znpu = fsnspp( glamu, gphiu, DOUBLE = double )
115  glamu = -1 & gphiu = -1; free memory
116  znpu.x = - znpu.x
117  znpu.y = - znpu.y
118  znnpu = znpu.x*znpu.x + znpu.y*znpu.y
119;   ... north pole direction & modulous (at v-point)
120  znpv = fsnspp( glamv, gphiv, DOUBLE = double )
121  znpv00 = znpv
122  znpv01 = fsnspp( shift(glamv, 0, 1), shift(gphiv, 0, 1), DOUBLE = double )
123  glamv = -1 & gphiv = -1; free memory
124  znpv.x = - znpv.x
125  znpv.y = - znpv.y
126  znnpv = znpv.x*znpv.x + znpv.y*znpv.y
127;   ... f-point
128  znpf00 = fsnspp( glamf, gphif, DOUBLE = double )
129  znpf01 = fsnspp( shift(glamf, 0, 1), shift(gphif, 0, 1), DOUBLE = double )
130  znpf10 = fsnspp( shift(glamf, 1, 0), shift(gphif, 1, 0), DOUBLE = double )
131  glamf = -1 & gphif = -1; free memory
132;   ... j-direction: v-point segment direction (t-point)
133  zxvvt = znpv00.x - znpv01.x
134  zyvvt = znpv00.y - znpv01.y
135  zmnpvt = sqrt ( temporary(znnpt) * ( zxvvt*zxvvt + zyvvt*zyvvt )  )
136  znpv00 = -1; free memory
137  znpv01 = -1; free memory
138  IF keyword_set(double) THEN zmnpvt = 1.e-14 > zmnpvt $
139  ELSE zmnpvt = 1.e-6 > zmnpvt
140;   ... j-direction: f-point segment direction (u-point)
141  zxffu = znpf00.x - znpf01.x
142  zyffu = znpf00.y - znpf01.y
143  zmnpfu = sqrt ( temporary(znnpu) * ( zxffu*zxffu + zyffu*zyffu )  )
144  znpf01 = -1; free memory
145  IF keyword_set(double) THEN zmnpfu = 1.e-14 > zmnpfu $
146  ELSE zmnpfu = 1.e-6 > zmnpfu
147;   ... i-direction: f-point segment direction (v-point)
148  zxffv = znpf00.x - znpf10.x
149  zyffv = znpf00.y - znpf10.y
150  znpf00 = -1 &  znpf10 = -1; free memory
151  zmnpfv = sqrt ( temporary(znnpv) * ( zxffv*zxffv + zyffv*zyffv )  )
152  IF keyword_set(double) THEN zmnpfv = 1.e-14 > zmnpfv $
153  ELSE zmnpfv = 1.e-6 > zmnpfv
154;   ... cosinus and sinus using scalar and vectorial products
155  gsint = ( znpt.x*zyvvt - znpt.y*zxvvt ) / zmnpvt
156  gcost = ( znpt.x*zxvvt + znpt.y*zyvvt ) / zmnpvt
157;   ... cosinus and sinus using scalar and vectorial products
158  gsinu = ( znpu.x*zyffu - znpu.y*zxffu ) / zmnpfu
159  gcosu = ( znpu.x*zxffu + znpu.y*zyffu ) / zmnpfu
160;   ... cosinus and sinus using scalar and vectorial products
161;       (caution, rotation of 90 degres)
162  gsinv =  ( znpv.x*zxffv + znpv.y*zyffv ) / zmnpfv
163  gcosv = -( znpv.x*zyffv - znpv.y*zxffv ) / zmnpfv
164;
165; II. Geographic mesh
166; ===================
167;
168;       bad = where(abs(glamf-shift(glamf, 0, 1)) LT 1.e-8)
169;       IF bad[0] NE -1 THEN BEGIN
170;         gcosu[bad] = 1.
171;         gsinu[bad] = 0.
172;       ENDIF
173;       bad = where(abs(gphif-shift(gphif, 1, 0)) LT 1.e-8)
174;       IF bad[0] NE -1 THEN BEGIN
175;         gcosv[bad] = 1.
176;         gsinv[bad] = 0.
177;       ENDIF
178;
179; III. Lateral boundary conditions
180; ================================
181;
182  gcost[*, 0] = gcost[*, 1]
183  gsint[*, 0] = gsint[*, 1]
184  gcosu[*, 0] = gcosu[*, 1]
185  gsinu[*, 0] = gsinu[*, 1]
186  gcosv[0, *] = gcosv[jpj-2, *]
187  gsinv[0, *] = gsinv[jpj-2, *]
188;
189  RETURN
190END
Note: See TracBrowser for help on using the repository browser.