1 | ;+ |
---|
2 | ; |
---|
3 | ; @file_comments |
---|
4 | ; compute the grid parameters (cm_4mesh) common |
---|
5 | ; |
---|
6 | ; domains sizes: |
---|
7 | ; --------------- |
---|
8 | ; jpi, jpj, jpk, jpiglo, jpjglo, jpkglo, jpidta, jpjdta, jpkdta |
---|
9 | ; |
---|
10 | ; domains positions regarding to the original grid: |
---|
11 | ; -------------------------------------------------- |
---|
12 | ; ixminmesh, ixmaxmesh, iyminmesh, iymaxmesh, izminmesh, izmaxmesh |
---|
13 | ; ixmindta, ixmaxdta, iymindta, iymaxdta, izmindta, izmaxdta |
---|
14 | ; |
---|
15 | ; horizontal parameters: |
---|
16 | ; ---------------------- |
---|
17 | ; glamt, glamf, gphit, gphit, e1t, e2t |
---|
18 | ; |
---|
19 | ; additional horizontal parameters if FULLCGRID keyword is defined: |
---|
20 | ; ----------------------------------------------------------------- |
---|
21 | ; glamu, glamv, gphiu, gphiv, e1u, e1v, e1f, e2u, e2v, e2f |
---|
22 | ; |
---|
23 | ; verticals parameters: |
---|
24 | ; --------------------- |
---|
25 | ; gdept, gdepw, e3t, e3w |
---|
26 | ; |
---|
27 | ; masks: |
---|
28 | ; ------ |
---|
29 | ; tmask |
---|
30 | ; |
---|
31 | ; additional masks if FULLCGRID keyword is defined: |
---|
32 | ; ------------------------------------------------- |
---|
33 | ; umaskred, vmaskred, fmaskredx, fmaskredy |
---|
34 | ; |
---|
35 | ; triangles_list: |
---|
36 | ; --------------- |
---|
37 | ; triangulation |
---|
38 | ; |
---|
39 | ; key_* variables: |
---|
40 | ; ---------------- |
---|
41 | ; key_onearth, key_periodic, key_shift, key_stride, key_partialstep, |
---|
42 | ; key_yreverse, key_zreverse, key_gridtype |
---|
43 | ; |
---|
44 | ; xxx related variables: |
---|
45 | ; ---------------------- |
---|
46 | ; ccmeshparameters, ccreadparameters |
---|
47 | ; |
---|
48 | ; @categories |
---|
49 | ; Grid |
---|
50 | ; |
---|
51 | ; @param startx {in}{optional}{type=scalar} |
---|
52 | ; x starting point, optional if [XY]AXIS keyword is used |
---|
53 | ; |
---|
54 | ; @param starty {in}{optional}{type=scalar} |
---|
55 | ; y starting point, optional if [XY]AXIS keyword is used |
---|
56 | ; |
---|
57 | ; @param stepxin {in}{optional}{type=scalar or vector} |
---|
58 | ; x direction step, optional if [XY]AXIS keyword is used, must be > 0 |
---|
59 | ; if stepxin is a vector nx is not used |
---|
60 | ; |
---|
61 | ; @param stepyin {in}{optional}{type=scalar or vector} |
---|
62 | ; y direction step, optional if [XY]AXIS keyword is used, |
---|
63 | ; could be > 0 (south to north) or < 0 (north to south) |
---|
64 | ; if stepyin is a vector ny is not used |
---|
65 | ; |
---|
66 | ; @param nxin {in}{optional}{type=scalar} |
---|
67 | ; number of points in x direction, |
---|
68 | ; optional if [XY]AXIS keyword is used or stepxin is a vector |
---|
69 | ; |
---|
70 | ; @param nyin {in}{optional}{type=scalar} |
---|
71 | ; number of points in y direction, |
---|
72 | ; optional if [XY]AXIS keyword is used or stepyin is a vector |
---|
73 | ; |
---|
74 | ; @keyword FULLCGRID {default=0}{type=scalar: 0 or 1} |
---|
75 | ; Activate to specify that you want to compute all the C grid parameters: |
---|
76 | ; definition of glam[uv], gphi[uv], e1[uvf], e2[uvf], [uv]maskred and |
---|
77 | ; fmaskred[xy] will be add to the default computations |
---|
78 | ; |
---|
79 | ; @keyword GLAMBOUNDARY {default=those defined in the file}{type=2 elements vector} |
---|
80 | ; Longitude boundaries that should be used to visualize the data. |
---|
81 | ; lon2 > lon1 |
---|
82 | ; lon2 - lon1 le 360 |
---|
83 | ; By default, the common (cm_4mesh) variable key_shift will be automatically |
---|
84 | ; defined according to GLAMBOUNDARY. |
---|
85 | ; |
---|
86 | ; @keyword MASK {default=array of 1}{type=2D or 3D array} |
---|
87 | ; Specify the land(0)/sea(1) mask |
---|
88 | ; |
---|
89 | ; @keyword ONEARTH {default=1}{type=scalar: 0 or 1} |
---|
90 | ; Force the manual definition of data localization on the earth or not |
---|
91 | ; 0) if the data are not on the earth |
---|
92 | ; 1) if the data are on earth (in that case we can for example use |
---|
93 | ; the labels 'longitude', 'latitude' in plots). |
---|
94 | ; The resulting value will be stored in the common (cm_4mesh) variable key_onearth |
---|
95 | ; ONEARTH=0 forces PERIODIC=0, SHIFT=0 and is cancelling GLAMBOUNDARY |
---|
96 | ; |
---|
97 | ; @keyword PERIODIC {default=computed by using the first line of glamt}{type=scalar: 0 or 1} |
---|
98 | ; Force the manual definition of the grid zonal periodicity. |
---|
99 | ; The resulting value will be stored in the common (cm_4mesh) variable key_periodic |
---|
100 | ; PERIODIC=0 forces SHIFT=0 |
---|
101 | ; |
---|
102 | ; @keyword PLAIN {default=0}{type=scalar: 0 or 1} |
---|
103 | ; Force YREVERSE=0, ZREVERSE=0, PERIODIC=0, SHIFT=0, STRIDE=[1, 1, 1] and |
---|
104 | ; suppress the automatic redefinition of the domain in case of x periodicity overlap, |
---|
105 | ; y periodicity overlap (ORCA type only) and mask border to 0. |
---|
106 | ; |
---|
107 | ; @keyword SHIFT {default=computed according to glamboundary}{type=scalar} |
---|
108 | ; Force the manual definition of the zonal shift that must be apply to the data. |
---|
109 | ; The resulting value will be stored in the common (cm_4mesh) variable key_shift |
---|
110 | ; Note that if key_periodic=0 then in any case key_shift=0. |
---|
111 | ; |
---|
112 | ; @keyword STRCALLING {type=string} |
---|
113 | ; a string containing the calling command used to |
---|
114 | ; call computegrid (this is used by <pro>xxx</pro>) |
---|
115 | ; |
---|
116 | ; @keyword STRIDE {default=[1, 1, 1]}{type=3 elements vector} |
---|
117 | ; Specify the stride in x, y and z direction. The resulting |
---|
118 | ; value will be stored in the common (cm_4mesh) variable key_stride |
---|
119 | ; |
---|
120 | ; @keyword XAXIS {type=1D or 2D array} |
---|
121 | ; Specify longitudes in this case startx, stepx and nx are not used but |
---|
122 | ; could be necessary if the y axis is not defined with yaxis. It must be |
---|
123 | ; possible to sort the first line of xaxis in the increasing order by |
---|
124 | ; shifting its elements. |
---|
125 | ; |
---|
126 | ; @keyword YAXIS {type=1D or 2D array} |
---|
127 | ; Specify latitudes in this case starty, stepy and ny are not used but |
---|
128 | ; starty and stepy could be necessary if the x axis is not defined with |
---|
129 | ; xaxis. It must be sorted in the increasing or decreasing order (along each column if 2d array). |
---|
130 | ; |
---|
131 | ; @keyword XYINDEX activate to specify that the horizontal grid should |
---|
132 | ; be simply defined by using the index of the points |
---|
133 | ; (xaxis = findgen(nx) and yaxis = findgen(ny)) |
---|
134 | ; using this keyword forces key_onearth=0 |
---|
135 | ; |
---|
136 | ; @keyword XMINMESH {default=0L}{type=scalar} |
---|
137 | ; Define common (cm_4mesh) variables ixminmesh used to define the localization |
---|
138 | ; of the first point of the grid along the x direction in a zoom of the original grid |
---|
139 | ; |
---|
140 | ; @keyword YMINMESH {default=0L}{type=scalar} |
---|
141 | ; Define common (cm_4mesh) variables iyminmesh used to define the localization |
---|
142 | ; of the first point of the grid along the y direction in a zoom of the original grid |
---|
143 | ; |
---|
144 | ; @keyword ZMINMESH {default=0L}{type=scalar} |
---|
145 | ; Define common (cm_4mesh) variables izminmesh used to define the localization |
---|
146 | ; of the first point of the grid along the z direction in a zoom of the original grid |
---|
147 | ; |
---|
148 | ; @keyword XMAXMESH {default=jpiglo-1}{type=scalar} |
---|
149 | ; Define common (cm_4mesh) variables ixmaxmesh used to define the localization |
---|
150 | ; of the last point of the grid along the x direction in a zoom of the original grid |
---|
151 | ; Note that if XMAXMESH < 0 then ixmaxmesh is defined as ixmaxmesh = jpiglo -1 + xmaxmesh |
---|
152 | ; |
---|
153 | ; @keyword YMAXMESH {default=jpjglo-1}{type=scalar} |
---|
154 | ; Define common (cm_4mesh) variables iymaxmesh used to define the localization |
---|
155 | ; of the last point of the grid along the y direction in a zoom of the original grid |
---|
156 | ; Note that if YMAXMESH < 0 then iymaxmesh is defined as iymaxmesh = jpjglo -1 + ymaxmesh |
---|
157 | ; |
---|
158 | ; @keyword ZMAXMESH {default=jpkglo-1}{type=scalar} |
---|
159 | ; Define common (cm_4mesh) variables izmaxmesh used to define the localization |
---|
160 | ; of the last point of the grid along the z direction in a zoom of the original grid |
---|
161 | ; Note that if ZMAXMESH < 0 then izmaxmesh is defined as izmaxmesh = jpkglo -1 + maxmesh |
---|
162 | ; |
---|
163 | ; @keyword FBASE2TBASE |
---|
164 | ; Activate when the model is a C grid based on a F point |
---|
165 | ; (with a F point at the bottom-left corner and a T point at the |
---|
166 | ; upper-right corner). In this case, we ignore |
---|
167 | ; - the first line of F and V points |
---|
168 | ; - the last line of T and U points |
---|
169 | ; - if the grid is not x-periodic, the first column of F and U points |
---|
170 | ; - if the grid is not x-periodic, the last column of T and V points. |
---|
171 | ; => we are back to a C grid based on T point as for OPA model. |
---|
172 | ; Note that in that case, key_gridtype = 'c_f' and not 'c' (-> used in read_ncdf) |
---|
173 | ; Note that activate FBASE2TBASE forces FULLCGRID=1 |
---|
174 | ; |
---|
175 | ; @keyword UBASE2TBASE |
---|
176 | ; Activate when the model is a C grid based on a U point |
---|
177 | ; (with a U point at the bottom-left corner and a T point at the |
---|
178 | ; upper-right corner). In this case, we ignore |
---|
179 | ; - if the grid is not x-periodic, the first column of F and U points |
---|
180 | ; - if the grid is not x-periodic, the last column of T and V points. |
---|
181 | ; => we are back to a C grid based on T point as for OPA model. |
---|
182 | ; Note that in that case, key_gridtype = 'c_u' and not 'c' (-> used in read_ncdf) |
---|
183 | ; Note that activate UBASE2TBASE forces FULLCGRID=1 |
---|
184 | ; |
---|
185 | ; @keyword VBASE2TBASE |
---|
186 | ; Activate when the model is a C grid based on a V point |
---|
187 | ; (with a V point at the bottom-left corner and a T point at the |
---|
188 | ; upper-right corner). In this case, we ignore |
---|
189 | ; - the first line of F and V points |
---|
190 | ; - the last line of T and U points |
---|
191 | ; => we are back to a C grid based on T point as for OPA model. |
---|
192 | ; Note that in that case, key_gridtype = 'c_v' and not 'c' (-> used in read_ncdf) |
---|
193 | ; Note that activate VBASE2TBASE forces FULLCGRID=1 |
---|
194 | ; |
---|
195 | ; @keyword ROMSH {type=2D array} |
---|
196 | ; This array is the final bathymetry at RHO-points. It is stored in the common |
---|
197 | ; variable (cm_4mesh) romszinfos.h |
---|
198 | ; Used when the model is a ROMS C-grid with one more point |
---|
199 | ; in longitude for T and V grid and one more point in latitude |
---|
200 | ; for T and U grid. In this case, we ignore |
---|
201 | ; - the last line of T and U points |
---|
202 | ; - the last column of T and V points. |
---|
203 | ; => we are back to a C grid based on T point as for OPA model. |
---|
204 | ; Note that activate ROMSH forces FULLCGRID=1 |
---|
205 | ; |
---|
206 | ; @keyword STRCALLING {type=scalar string} |
---|
207 | ; Used by xxx... |
---|
208 | ; |
---|
209 | ; @keyword YREVERSE {default=computed according to gphit[0, 1:jpj-1] LT gphit[0, 0:jpj-2]}{type=scalar} |
---|
210 | ; Force the manual definition of the y reverse that must be apply to the data. |
---|
211 | ; The resulting value will be stored in the common (cm_4mesh) variable key_yreverse |
---|
212 | ; |
---|
213 | ; @keyword ZAXIS {type=1D} |
---|
214 | ; Specify the vertical axis. Must be sorted in the increasing or decreasing order |
---|
215 | ; |
---|
216 | ; @keyword ZREVERSE {default=computed according to gdept[0] GT gdept[1]}{type=scalar} |
---|
217 | ; Force the manual definition of the z reverse that must be apply to the data. |
---|
218 | ; The resulting value will be stored in the common (cm_4mesh) variable key_zreverse |
---|
219 | ; |
---|
220 | ; @keyword _EXTRA |
---|
221 | ; not used in the present case ... |
---|
222 | ; |
---|
223 | ; @uses |
---|
224 | ; cm_4mesh |
---|
225 | ; cm_4data |
---|
226 | ; cm_4cal |
---|
227 | ; |
---|
228 | |
---|
229 | ; @restrictions |
---|
230 | ; if the grid has x/y periodicity overlap and/or if |
---|
231 | ; the mask has 0 everywhere at the border (like a closed sea) and |
---|
232 | ; if (we did not activate /plain and xminmesh, xmaxmesh, yminmesh, |
---|
233 | ; ymaxmesh keywords are defined to their default values), we redefine |
---|
234 | ; xminmesh, xmaxmesh, yminmesh, ymaxmesh in order to reove the |
---|
235 | ; overlapping part and/or to open the domain (avoid it be forced |
---|
236 | ; to use cell_fill = 1). |
---|
237 | ; |
---|
238 | ; FUV points definition is not exact if the grid is irregular |
---|
239 | ; |
---|
240 | ; @history |
---|
241 | ; Sebastien Masson (smasson\@lodyc.jussieu.fr) |
---|
242 | ; 2000-04-20 |
---|
243 | ; Sept 2004, several bug fix to suit C grid type... |
---|
244 | ; Aug 2005, rewritte almost everything... |
---|
245 | ; |
---|
246 | ; @version |
---|
247 | ; $Id$ |
---|
248 | ; |
---|
249 | ;- |
---|
250 | PRO computegrid, startx, starty, stepxin, stepyin, nxin, nyin $ |
---|
251 | , XAXIS=xaxis, YAXIS=yaxis, ZAXIS=zaxis $ |
---|
252 | , MASK=mask, GLAMBOUNDARY=glamboundary $ |
---|
253 | , XMINMESH=xminmesh, XMAXMESH=xmaxmesh $ |
---|
254 | , YMINMESH=yminmesh, YMAXMESH=ymaxmesh $ |
---|
255 | , ZMINMESH=zminmesh, ZMAXMESH=zmaxmesh $ |
---|
256 | , ONEARTH=onearth, PERIODIC=periodic $ |
---|
257 | , PLAIN=plain, SHIFT=shift, STRIDE=stride $ |
---|
258 | , YREVERSE=yreverse, ZREVERSE=zreverse $ |
---|
259 | , FULLCGRID=fullcgrid, XYINDEX=xyindex $ |
---|
260 | , UBASE2TBASE=ubase2tbase, VBASE2TBASE=vbase2tbase $ |
---|
261 | , FBASE2TBASE=fbase2tbase $ |
---|
262 | , STRCALLING=strcalling, ROMSH=romsh, _EXTRA=ex |
---|
263 | ; |
---|
264 | compile_opt idl2, strictarrsubs |
---|
265 | ; |
---|
266 | @cm_4mesh |
---|
267 | @cm_4data |
---|
268 | @cm_4cal |
---|
269 | IF NOT keyword_set(key_forgetold) THEN BEGIN |
---|
270 | @updatenew |
---|
271 | @updatekwd |
---|
272 | ENDIF |
---|
273 | ;--------------------------------------------------------- |
---|
274 | ;------------------------------------------------------------ |
---|
275 | time1 = systime(1) ; for key_performance |
---|
276 | ;------------------------------------------------------------ |
---|
277 | ; |
---|
278 | ;==================================================== |
---|
279 | ; Check input parameters |
---|
280 | ;==================================================== |
---|
281 | ; |
---|
282 | ; xaxis related parameters |
---|
283 | ; |
---|
284 | if n_elements(xaxis) NE 0 then BEGIN |
---|
285 | CASE (size(xaxis))[0] OF |
---|
286 | 0:nx = 1L |
---|
287 | 1:nx = (size(xaxis))[1] |
---|
288 | 2:nx = (size(xaxis))[1] |
---|
289 | ENDCASE |
---|
290 | ENDIF ELSE BEGIN |
---|
291 | IF n_elements(startx) EQ 0 THEN BEGIN |
---|
292 | dummy = report('If xaxis is not given, startx must be defined') |
---|
293 | return |
---|
294 | ENDIF |
---|
295 | CASE n_elements(stepxin) OF |
---|
296 | 0:BEGIN |
---|
297 | dummy = report('If xaxis is not given, stepxin must be defined') |
---|
298 | return |
---|
299 | END |
---|
300 | 1:BEGIN |
---|
301 | IF n_elements(nxin) EQ 0 THEN BEGIN |
---|
302 | dummy = report('If xaxis is not given and stepxin has only one element, nx must be defined') |
---|
303 | return |
---|
304 | ENDIF ELSE nx = nxin |
---|
305 | END |
---|
306 | ELSE:nx = n_elements(stepxin) |
---|
307 | ENDCASE |
---|
308 | ENDELSE |
---|
309 | ; |
---|
310 | ; yaxis related parameters |
---|
311 | ; |
---|
312 | if n_elements(yaxis) NE 0 then BEGIN |
---|
313 | CASE (size(yaxis))[0] OF |
---|
314 | 0:ny = 1L |
---|
315 | 1:ny = (size(yaxis))[1] |
---|
316 | 2:ny = (size(yaxis))[2] |
---|
317 | ENDCASE |
---|
318 | ENDIF ELSE BEGIN |
---|
319 | IF n_elements(starty) EQ 0 THEN BEGIN |
---|
320 | dummy = report('If yaxis is not given, starty must be defined') |
---|
321 | return |
---|
322 | ENDIF |
---|
323 | CASE n_elements(stepyin) OF |
---|
324 | 0:BEGIN |
---|
325 | dummy = report('If yaxis is not given, stepyin must be defined') |
---|
326 | return |
---|
327 | END |
---|
328 | 1:BEGIN |
---|
329 | IF n_elements(nyin) EQ 0 THEN BEGIN |
---|
330 | dummy = report('If yaxis is not given and stepyin has only one element, ny must be defined') |
---|
331 | return |
---|
332 | ENDIF ELSE ny = nyin |
---|
333 | END |
---|
334 | ELSE:ny = n_elements(stepyin) |
---|
335 | ENDCASE |
---|
336 | ENDELSE |
---|
337 | ; |
---|
338 | ; zaxis related parameters |
---|
339 | ; |
---|
340 | if n_elements(zaxis) NE 0 then BEGIN |
---|
341 | CASE (size(zaxis))[0] OF |
---|
342 | 0:nz = 1L |
---|
343 | 1:nz = (size(zaxis))[1] |
---|
344 | ELSE:BEGIN |
---|
345 | ras = report( 'not coded') |
---|
346 | stop |
---|
347 | END |
---|
348 | ENDCASE |
---|
349 | ENDIF ELSE nz = 1L |
---|
350 | ; |
---|
351 | ;==================================================== |
---|
352 | ; Others automatic definitions... |
---|
353 | ;==================================================== |
---|
354 | ; |
---|
355 | jpiglo = long(nx) |
---|
356 | jpjglo = long(ny) |
---|
357 | jpkglo = long(nz) |
---|
358 | IF keyword_set(romsh) THEN BEGIN |
---|
359 | jpiglo = jpiglo - 1 |
---|
360 | jpjglo = jpjglo - 1 |
---|
361 | fullcgrid = 1 |
---|
362 | ENDIF |
---|
363 | ; |
---|
364 | ; impact of plain keyword: |
---|
365 | ; |
---|
366 | IF keyword_set(plain) THEN BEGIN |
---|
367 | yreverse = 0 |
---|
368 | zreverse = 0 |
---|
369 | periodic = 0 |
---|
370 | shift = 0 |
---|
371 | stride = [1, 1, 1] |
---|
372 | ENDIF |
---|
373 | ; |
---|
374 | IF n_elements(xminmesh) NE 0 THEN ixminmesh = long(xminmesh[0]) ELSE ixminmesh = 0l |
---|
375 | IF n_elements(xmaxmesh) NE 0 THEN ixmaxmesh = long(xmaxmesh[0]) ELSE ixmaxmesh = jpiglo-1 |
---|
376 | IF n_elements(yminmesh) NE 0 THEN iyminmesh = long(yminmesh[0]) ELSE iyminmesh = 0l |
---|
377 | IF n_elements(ymaxmesh) NE 0 THEN iymaxmesh = long(ymaxmesh[0]) ELSE iymaxmesh = jpjglo-1 |
---|
378 | IF n_elements(zminmesh) NE 0 THEN izminmesh = long(zminmesh[0]) ELSE izminmesh = 0l |
---|
379 | IF n_elements(zmaxmesh) NE 0 THEN izmaxmesh = long(zmaxmesh[0]) ELSE izmaxmesh = jpkglo-1 |
---|
380 | ; |
---|
381 | CASE 1 OF |
---|
382 | keyword_set(fbase2tbase):key_gridtype = 'c_f' |
---|
383 | keyword_set(ubase2tbase):key_gridtype = 'c_u' |
---|
384 | keyword_set(vbase2tbase):key_gridtype = 'c_v' |
---|
385 | else:key_gridtype = 'c' |
---|
386 | ENDCASE |
---|
387 | IF key_gridtype EQ 'c_v' OR key_gridtype EQ 'c_f' THEN BEGIN |
---|
388 | iymaxmesh = iymaxmesh-1 |
---|
389 | ENDIF |
---|
390 | IF strlen(key_gridtype) EQ 3 THEN fullcgrid = 1 |
---|
391 | ; |
---|
392 | IF ixmaxmesh LT 0 THEN ixmaxmesh = jpiglo -1 + ixmaxmesh |
---|
393 | IF iymaxmesh LT 0 THEN iymaxmesh = jpjglo -1 + iymaxmesh |
---|
394 | IF izmaxmesh LT 0 THEN izmaxmesh = jpkglo -1 + izmaxmesh |
---|
395 | ; avoid basics errors... |
---|
396 | ixmaxmesh = 0 > ixmaxmesh < (jpiglo-1) |
---|
397 | ixminmesh = 0 > ixminmesh < ixmaxmesh |
---|
398 | iymaxmesh = 0 > iymaxmesh < (jpjglo-1) |
---|
399 | iyminmesh = 0 > iyminmesh < iymaxmesh |
---|
400 | izmaxmesh = 0 > izmaxmesh < (jpkglo-1) |
---|
401 | izminmesh = 0 > izminmesh < izmaxmesh |
---|
402 | ; |
---|
403 | jpi = ixmaxmesh-ixminmesh+1 |
---|
404 | jpj = iymaxmesh-iyminmesh+1 |
---|
405 | jpk = izmaxmesh-izminmesh+1 |
---|
406 | ; |
---|
407 | jpidta = jpiglo |
---|
408 | jpjdta = jpjglo |
---|
409 | jpkdta = jpkglo |
---|
410 | ixmindta = 0 |
---|
411 | ixmaxdta = jpidta-1 |
---|
412 | iymindta = 0 |
---|
413 | iymaxdta = jpjdta-1 |
---|
414 | izmindta = 0 |
---|
415 | izmaxdta = jpkdta-1 |
---|
416 | ; |
---|
417 | key_partialstep = 0 |
---|
418 | if n_elements(stride) eq 3 then key_stride = stride $ |
---|
419 | ELSE key_stride = [1, 1, 1] |
---|
420 | ; |
---|
421 | ; check xyindex and its consequences |
---|
422 | ; |
---|
423 | if keyword_set(xyindex) then onearth = 0 |
---|
424 | ; |
---|
425 | ; check onearth and its consequences |
---|
426 | ; |
---|
427 | IF n_elements(onearth) EQ 0 THEN key_onearth = 1b $ |
---|
428 | ELSE key_onearth = keyword_set(onearth) |
---|
429 | IF NOT key_onearth THEN BEGIN |
---|
430 | periodic = 0 |
---|
431 | shift = 0 |
---|
432 | ENDIF |
---|
433 | |
---|
434 | r = 6371000. |
---|
435 | ; |
---|
436 | ;==================================================== |
---|
437 | ; X direction : glamt |
---|
438 | ;==================================================== |
---|
439 | ; |
---|
440 | ; def of glamt |
---|
441 | ; |
---|
442 | if n_elements(xaxis) NE 0 then BEGIN |
---|
443 | if keyword_set(xyindex) THEN glamt = findgen(jpiglo) ELSE glamt = xaxis |
---|
444 | ENDIF ELSE BEGIN |
---|
445 | if keyword_set(xyindex) THEN stepx = 1. ELSE stepx = stepxin |
---|
446 | CASE 1 OF |
---|
447 | n_elements(stepx):glamt = startx + findgen(jpiglo)*stepx |
---|
448 | size(stepx, /n_dimensions):glamt = startx + total(stepx, /cumulative) |
---|
449 | ELSE:BEGIN |
---|
450 | dummy = report('Wrong definition of stepx...') |
---|
451 | return |
---|
452 | END |
---|
453 | ENDCASE |
---|
454 | ENDELSE |
---|
455 | ; |
---|
456 | ; apply glamboundary |
---|
457 | ; |
---|
458 | IF keyword_set(glamboundary) AND key_onearth THEN BEGIN |
---|
459 | IF glamboundary[0] GE glamboundary[1] THEN stop |
---|
460 | IF glamboundary[1]-glamboundary[0] GT 360 THEN stop |
---|
461 | glamt = glamt MOD 360 |
---|
462 | smaller = where(glamt LT glamboundary[0]) |
---|
463 | if smaller[0] NE -1 then glamt[smaller] = glamt[smaller]+360. |
---|
464 | bigger = where(glamt GE glamboundary[1]) |
---|
465 | if bigger[0] NE -1 then glamt[bigger] = glamt[bigger]-360. |
---|
466 | ENDIF |
---|
467 | ; |
---|
468 | ; force glamt to have 2 dimensions |
---|
469 | ; |
---|
470 | CASE size(reform(glamt), /n_dimensions) OF |
---|
471 | 0:glamt = replicate(glamt, jpi, jpj) |
---|
472 | 1:glamt = glamt[ixminmesh:ixmaxmesh]#replicate(1, jpj) |
---|
473 | 2:glamt = glamt[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh] |
---|
474 | ENDCASE |
---|
475 | ; keep 2d array even with degenerated dimension |
---|
476 | IF jpj EQ 1 THEN glamt = reform(glamt, jpi, jpj, /over) |
---|
477 | ; |
---|
478 | ;==================================================== |
---|
479 | ; Y direction : gphit |
---|
480 | ;==================================================== |
---|
481 | ; |
---|
482 | ; def of gphit |
---|
483 | ; |
---|
484 | if n_elements(yaxis) NE 0 THEN BEGIN |
---|
485 | if keyword_set(xyindex) THEN gphit = findgen(jpjglo) ELSE gphit = yaxis |
---|
486 | ENDIF ELSE BEGIN |
---|
487 | if keyword_set(xyindex) THEN stepy = 1. ELSE stepy = stepyin |
---|
488 | CASE 1 OF |
---|
489 | n_elements(stepy):gphit = starty + findgen(jpjglo)*stepy |
---|
490 | size(stepy, /n_dimensions):gphit = starty + total(stepy, /cumulative) |
---|
491 | ELSE:BEGIN |
---|
492 | dummy = report('Wrong definition of stepy...') |
---|
493 | return |
---|
494 | END |
---|
495 | ENDCASE |
---|
496 | ENDELSE |
---|
497 | ; |
---|
498 | ; force gphit to have 2 dimensions |
---|
499 | ; |
---|
500 | CASE size(reform(gphit), /n_dimensions) OF |
---|
501 | 0:gphit = replicate(gphit, jpi, jpj) |
---|
502 | 1:gphit = replicate(1, jpi)#gphit[iyminmesh:iymaxmesh] |
---|
503 | 2:gphit = gphit[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh] |
---|
504 | ENDCASE |
---|
505 | ; keep 2d array even with degenerated dimension |
---|
506 | IF jpj EQ 1 THEN gphit = reform(gphit, jpi, jpj, /over) |
---|
507 | ; |
---|
508 | ;==================================================== |
---|
509 | ; check y periodicity... Only according to ORCA grid |
---|
510 | ;==================================================== |
---|
511 | ; check the periodicity if iyminmesh and iymaxmesh have the default definitions... |
---|
512 | IF NOT keyword_set(plain) AND key_onearth EQ 1 $ |
---|
513 | AND iyminmesh EQ 0l AND iymaxmesh eq jpjglo-1 AND jpj GE 3 AND jpi GE 2 THEN BEGIN |
---|
514 | |
---|
515 | CASE 1 OF |
---|
516 | ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $ |
---|
517 | AND array_equal(gphit[1:*, jpj-1], reverse(gphit[1:*, jpj-3])) EQ 1:BEGIN |
---|
518 | ; T pivot |
---|
519 | ymaxmesh = -1 |
---|
520 | recall = 1 |
---|
521 | END |
---|
522 | ixminmesh EQ 1l AND ixmaxmesh eq jpiglo-2 $ |
---|
523 | AND array_equal(gphit[*, jpj-1], reverse(shift(gphit[*, jpj-3], -1))) EQ 1:BEGIN |
---|
524 | ; T pivot |
---|
525 | ymaxmesh = -1 |
---|
526 | recall = 1 |
---|
527 | END |
---|
528 | ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $ |
---|
529 | AND array_equal(gphit[*, jpj-1], reverse(gphit[*, jpj-2])) EQ 1:BEGIN |
---|
530 | ; F pivot |
---|
531 | ymaxmesh = -1 |
---|
532 | recall = 1 |
---|
533 | END |
---|
534 | ixminmesh EQ 1l AND ixmaxmesh eq jpiglo-2 $ |
---|
535 | AND array_equal(gphit[*, jpj-1], reverse(gphit[*, jpj-2])) EQ 1:BEGIN |
---|
536 | ; F pivot |
---|
537 | ymaxmesh = -1 |
---|
538 | recall = 1 |
---|
539 | END |
---|
540 | ELSE: |
---|
541 | ENDCASE |
---|
542 | ENDIF |
---|
543 | ; |
---|
544 | ;==================================================== |
---|
545 | ; check x periodicity... |
---|
546 | ;==================================================== |
---|
547 | IF n_elements(periodic) NE 0 THEN forcenoperio = 1 - keyword_set(periodic) |
---|
548 | ; check the periodicity if ixminmesh and ixmaxmesh have the default definitions... |
---|
549 | IF NOT keyword_set(plain) AND NOT keyword_set(forcenoperio) AND key_onearth EQ 1 $ |
---|
550 | AND ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 AND jpi GE 3 THEN BEGIN |
---|
551 | CASE 0 OF |
---|
552 | total((glamt[0, *] - glamt[jpi-2, *]) MOD 360) $ |
---|
553 | + total((glamt[1, *] - glamt[jpi-1, *]) MOD 360):BEGIN |
---|
554 | xminmesh = 1 |
---|
555 | xmaxmesh = -1 |
---|
556 | recall = 1 |
---|
557 | END |
---|
558 | total((glamt[0, *] - glamt[jpi-2, *]) MOD 360):BEGIN |
---|
559 | xminmesh = 1 |
---|
560 | recall = 1 |
---|
561 | END |
---|
562 | total((glamt[1, *] - glamt[jpi-1, *]) MOD 360):BEGIN |
---|
563 | xmaxmesh = -1 |
---|
564 | recall = 1 |
---|
565 | END |
---|
566 | ELSE: |
---|
567 | ENDCASE |
---|
568 | ENDIF |
---|
569 | ;==================================================== |
---|
570 | ; recall computegrid if needed... |
---|
571 | ;==================================================== |
---|
572 | IF keyword_set(recall) THEN BEGIN |
---|
573 | computegrid, XAXIS = glamt, YAXIS = gphit, ZAXIS = zaxis $ |
---|
574 | , MASK = mask, GLAMBOUNDARY = glamboundary $ |
---|
575 | , XMINMESH = xminmesh, XMAXMESH = xmaxmesh $ |
---|
576 | , YMINMESH = yminmesh, YMAXMESH = ymaxmesh $ |
---|
577 | , ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $ |
---|
578 | , PERIODIC = periodic, SHIFT = shift, STRIDE = stride $ |
---|
579 | , FULLCGRID = fullcgrid, XYINDEX = xyindex $ |
---|
580 | , STRCALLING = strcalling $ |
---|
581 | , ROMSH = romsh, _extra = ex |
---|
582 | return |
---|
583 | ENDIF |
---|
584 | ; |
---|
585 | ;==================================================== |
---|
586 | ; def key_yreverse |
---|
587 | ;==================================================== |
---|
588 | ; |
---|
589 | IF n_elements(yreverse) EQ 0 THEN BEGIN |
---|
590 | IF jpj GT 1 THEN BEGIN |
---|
591 | IF total(gphit[0, 1:jpj-1] LT gphit[0, 0:jpj-2]) GT jpj/2 THEN key_yreverse = 1 ELSE key_yreverse = 0 |
---|
592 | ENDIF ELSE key_yreverse = 0 |
---|
593 | ENDIF ELSE key_yreverse = yreverse |
---|
594 | IF keyword_set(key_yreverse) THEN BEGIN |
---|
595 | gphit = reverse(gphit, 2) |
---|
596 | glamt = reverse(glamt, 2) |
---|
597 | ENDIF |
---|
598 | ; |
---|
599 | ;==================================================== |
---|
600 | ; def of key_shift |
---|
601 | ;==================================================== |
---|
602 | ; |
---|
603 | ; definition of key_shift by shifting the array to have the min |
---|
604 | ; values of glamt[*, 0] in glamt[0, 0] |
---|
605 | ; |
---|
606 | IF n_elements(shift) EQ 0 THEN BEGIN |
---|
607 | IF jpi GT 1 then BEGIN |
---|
608 | xtest = glamt[*, 0] |
---|
609 | key_shift = (where(xtest EQ min(xtest)))[0] |
---|
610 | IF key_shift NE 0 THEN key_shift = jpi - key_shift |
---|
611 | ENDIF ELSE key_shift = 0 |
---|
612 | ENDIF ELSE key_shift = shift |
---|
613 | ; |
---|
614 | ;==================================================== |
---|
615 | ; def of key_periodic |
---|
616 | ;==================================================== |
---|
617 | ; |
---|
618 | IF n_elements(periodic) EQ 0 THEN BEGIN |
---|
619 | IF jpi GT 1 THEN BEGIN |
---|
620 | xtest = shift(glamt[*, 0], key_shift) |
---|
621 | ; check that xtest is now sorted in the increasing order |
---|
622 | IF array_equal(sort(xtest), lindgen(jpi)) EQ 0 THEN BEGIN |
---|
623 | ras = report(['WARNING: we cannot sort the xaxis with a simple shift...', $ |
---|
624 | 'we force key_periodic = 0 and key_shift = 0', $ |
---|
625 | 'only horizontal plot may be ok...']) |
---|
626 | key_periodic = 0 |
---|
627 | xnotsorted = 1 |
---|
628 | ENDIF ELSE BEGIN |
---|
629 | key_periodic = (xtest[jpi-1]+2*(xtest[jpi-1]-xtest[jpi-2])) $ |
---|
630 | GE (xtest[0]+360) |
---|
631 | ENDELSE |
---|
632 | ENDIF ELSE key_periodic = 0 |
---|
633 | ENDIF ELSE key_periodic = keyword_set(periodic) |
---|
634 | ; |
---|
635 | ; update key_shift |
---|
636 | ; |
---|
637 | key_shift = key_shift * (key_periodic EQ 1) |
---|
638 | ; |
---|
639 | IF (key_gridtype EQ 'c_u' OR key_gridtype EQ 'c_f') AND NOT keyword_set(key_periodic) THEN BEGIN |
---|
640 | ixmaxmesh = ixmaxmesh-1 |
---|
641 | jpi = jpi-1 |
---|
642 | ENDIF |
---|
643 | ; |
---|
644 | ;==================================================== |
---|
645 | ; apply key_shift |
---|
646 | ;==================================================== |
---|
647 | ; |
---|
648 | if keyword_set(key_shift) then BEGIN |
---|
649 | glamt = shift(glamt, key_shift, 0) |
---|
650 | gphit = shift(gphit, key_shift, 0) |
---|
651 | IF jpj EQ 1 THEN BEGIN |
---|
652 | glamt = reform(glamt, jpi, jpj, /over) |
---|
653 | gphit = reform(gphit, jpi, jpj, /over) |
---|
654 | ENDIF |
---|
655 | ENDIF |
---|
656 | ; |
---|
657 | ;==================================================== |
---|
658 | ; Are we using a "regular" grid (that can be described |
---|
659 | ; with x vector and y vector)? |
---|
660 | ;==================================================== |
---|
661 | ; |
---|
662 | ; to get faster, we first test the most basic cases before |
---|
663 | ; testing the full array. |
---|
664 | ; |
---|
665 | CASE 1 OF |
---|
666 | keyword_set(xyindex):key_irregular = 0b |
---|
667 | jpi EQ 1 OR jpj EQ 1:key_irregular = 0b |
---|
668 | n_elements(xaxis) EQ 0 AND n_elements(yaxis) EQ 0:key_irregular = 0b |
---|
669 | size(reform(xaxis), /n_dimensions) EQ 1 AND size(reform(xaxis), /n_dimensions) EQ 1:key_irregular = 0b |
---|
670 | n_elements(xaxis) EQ 0 AND size(reform(yaxis), /n_dimensions) EQ 1:key_irregular = 0b |
---|
671 | n_elements(yaxis) EQ 0 AND size(reform(xaxis), /n_dimensions) EQ 1:key_irregular = 0b |
---|
672 | array_equal(glamt[*, 0], glamt[*, jpj-1]) EQ 0:key_irregular = 1b |
---|
673 | array_equal(gphit[0, *], gphit[jpi-1, *]) EQ 0:key_irregular = 1b |
---|
674 | array_equal(glamt, glamt[*, 0]#replicate(1, jpj)) EQ 0:key_irregular = 1b |
---|
675 | array_equal(gphit, replicate(1, jpi)#(gphit[0, *])[*]) EQ 0:key_irregular = 1b |
---|
676 | ELSE:key_irregular = 0b |
---|
677 | ENDCASE |
---|
678 | ; |
---|
679 | ;==================================================== |
---|
680 | ; def of glamf: defined as the middle of T(i,j) T(i+1,j+1) |
---|
681 | ;==================================================== |
---|
682 | ; |
---|
683 | IF jpi GT 1 THEN BEGIN |
---|
684 | ; we must compute stepxf: x distance between T(i,j) T(i+1,j+1) |
---|
685 | CASE 1 OF |
---|
686 | n_elements(stepx):stepxf = stepx |
---|
687 | size(stepx, /n_dimensions):stepxf = stepx#replicate(1, jpj) |
---|
688 | ELSE:BEGIN |
---|
689 | if (keyword_set(key_onearth) AND keyword_set(xnotsorted)) $ |
---|
690 | OR (keyword_set(key_periodic) AND key_irregular) then BEGIN |
---|
691 | stepxf = (glamt + 720) MOD 360 |
---|
692 | IF jpj EQ 1 THEN stepxf = reform(stepxf, jpi, jpj, /over) |
---|
693 | stepxf = shift(stepxf, -1, -1) - stepxf |
---|
694 | stepxf = [ [[stepxf]], [[stepxf + 360]], [[stepxf - 360]] ] |
---|
695 | stepxf = min(abs(stepxf), dimension = 3) |
---|
696 | IF NOT keyword_set(key_periodic) THEN $ |
---|
697 | stepxf[jpi-1, *] = stepxf[jpi-2, *] |
---|
698 | ENDIF ELSE BEGIN |
---|
699 | stepxf = shift(glamt, -1, -1) - glamt |
---|
700 | IF keyword_set(key_periodic) THEN $ |
---|
701 | stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $ |
---|
702 | ELSE stepxf[jpi-1, *] = stepxf[jpi-2, *] |
---|
703 | ENDELSE |
---|
704 | IF jpj GT 1 THEN BEGIN |
---|
705 | stepxf[*, jpj-1] = stepxf[*, jpj-2] |
---|
706 | stepxf[jpi-1, jpj-1] = stepxf[jpi-2, jpj-2] |
---|
707 | ENDIF |
---|
708 | END |
---|
709 | ENDCASE |
---|
710 | glamf = glamt + 0.5 * stepxf |
---|
711 | ENDIF ELSE glamf = glamt + 0.5 |
---|
712 | ; |
---|
713 | IF keyword_set(key_periodic) AND (max(glamf)-min(glamt)) GE 360 THEN BEGIN |
---|
714 | IF NOT keyword_set(glamboundary) THEN BEGIN |
---|
715 | bigger = where(glamf GE min(glamt)+360) |
---|
716 | glamf[bigger] = glamf[bigger]-360. |
---|
717 | ENDIF ELSE glamf = glamboundary[0] > temporary(glamf) < glamboundary[1] |
---|
718 | ENDIF |
---|
719 | ; |
---|
720 | IF jpj EQ 1 THEN glamf = reform(glamf, jpi, jpj, /over) |
---|
721 | ; |
---|
722 | ;==================================================== |
---|
723 | ; def of gphif: defined as the middle of T(i,j) T(i+1,j+1) |
---|
724 | ;==================================================== |
---|
725 | ; |
---|
726 | IF jpj GT 1 THEN BEGIN |
---|
727 | ; we must compute stepyf: y distance between T(i,j) T(i+1,j+1) |
---|
728 | CASE 1 OF |
---|
729 | n_elements(stepy):stepyf = stepy |
---|
730 | size(stepy, /n_dimensions):stepyf = replicate(1, jpi)#stepy |
---|
731 | ELSE:BEGIN |
---|
732 | stepyf = shift(gphit, -1, -1) - gphit |
---|
733 | stepyf[*, jpj-1] = stepyf[*, jpj-2] |
---|
734 | IF jpi GT 1 THEN BEGIN |
---|
735 | if NOT keyword_set(key_periodic) THEN $ |
---|
736 | stepyf[jpi-1, *] = stepyf[jpi-2, *] |
---|
737 | stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2] |
---|
738 | ENDIF |
---|
739 | END |
---|
740 | ENDCASE |
---|
741 | gphif = gphit + 0.5 * stepyf |
---|
742 | ENDIF ELSE gphif = gphit + 0.5 |
---|
743 | IF key_onearth THEN gphif = -90. > gphif < 90. |
---|
744 | ; |
---|
745 | IF jpj EQ 1 THEN gphif = reform(gphif, jpi, jpj, /over) |
---|
746 | ; |
---|
747 | ;==================================================== |
---|
748 | ; e1t: x distance between U(i-1,j) and U(i,j) |
---|
749 | ;==================================================== |
---|
750 | ; |
---|
751 | ; *-|-*---|---*---| |
---|
752 | ; |
---|
753 | IF jpi GT 1 THEN BEGIN |
---|
754 | IF n_elements(stepx) NE 1 THEN BEGIN |
---|
755 | IF keyword_set(irregular) THEN BEGIN |
---|
756 | ; we must compute stepxu: x distance between T(i,j) T(i+1,j) |
---|
757 | IF keyword_set(key_periodic) THEN BEGIN |
---|
758 | stepxu = (glamt + 720) MOD 360 |
---|
759 | stepxu = shift(stepxu, -1, 0) - stepxu |
---|
760 | stepxu = [ [[stepxu]], [[stepxu + 360]], [[stepxu - 360]] ] |
---|
761 | stepxu = min(abs(stepxu), dimension = 3) |
---|
762 | ENDIF ELSE BEGIN |
---|
763 | stepxu = shift(glamt, -1, 0) - glamt |
---|
764 | stepxu[jpi-1, *] = stepxf[jpi-2, *] |
---|
765 | ENDELSE |
---|
766 | ENDIF ELSE stepxu = stepxf |
---|
767 | IF jpj EQ 1 THEN stepxu = reform(stepxu, jpi, jpj, /over) |
---|
768 | e1t = 0.5*(stepxu+shift(stepxu, 1, 0)) |
---|
769 | IF NOT keyword_set(key_periodic) THEN $ |
---|
770 | e1t[0, *] = e1t[1, *] |
---|
771 | ENDIF ELSE e1t = replicate(stepx, jpi, jpj) |
---|
772 | ENDIF ELSE e1t = replicate(1b, jpi, jpj) |
---|
773 | ; |
---|
774 | IF jpj EQ 1 THEN e1t = reform(e1t, jpi, jpj, /over) |
---|
775 | ; |
---|
776 | ;==================================================== |
---|
777 | ; e2t: y distance between V(i,j-1) and V(i,j) |
---|
778 | ;==================================================== |
---|
779 | ; |
---|
780 | IF jpj GT 1 THEN BEGIN |
---|
781 | ; we must compute stepyv: y distance between T(i,j) T(i,j+1) |
---|
782 | IF n_elements(stepy) NE 1 THEN BEGIN |
---|
783 | IF keyword_set(key_irregular) THEN BEGIN |
---|
784 | stepyv = shift(gphit, 0, -1) - gphit |
---|
785 | stepyv[*, jpj-1] = stepyv[*, jpj-2] |
---|
786 | ENDIF ELSE stepyv = stepyf |
---|
787 | e2t = 0.5*(stepyv+shift(stepyv, 0, 1)) |
---|
788 | e2t[*, 0] = e2t[*, 1] |
---|
789 | ENDIF ELSE e2t = replicate(stepy, jpi, jpj) |
---|
790 | ENDIF ELSE e2t = replicate(1b, jpi, jpj) |
---|
791 | ; |
---|
792 | IF key_onearth THEN e2t = r * !pi/180. * temporary(e2t) |
---|
793 | ; |
---|
794 | IF jpj EQ 1 THEN e2t = reform(e2t, jpi, jpj, /over) |
---|
795 | ; |
---|
796 | ;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
797 | IF keyword_set(fullcgrid) THEN BEGIN |
---|
798 | ; |
---|
799 | ;==================================================== |
---|
800 | ; def of glamu: defined as the middle of T(i,j) T(i+1,j) |
---|
801 | ;==================================================== |
---|
802 | ; |
---|
803 | IF keyword_set(irregular) THEN BEGIN |
---|
804 | glamu = glamt + 0.5 * stepxu |
---|
805 | IF keyword_set(glamboundary) AND key_onearth THEN $ |
---|
806 | glamu = glamboundary[0] > temporary(glamu) < glamboundary[1] |
---|
807 | ENDIF ELSE glamu = glamf |
---|
808 | ; |
---|
809 | IF jpj EQ 1 THEN glamu = reform(glamu, jpi, jpj, /over) |
---|
810 | ; |
---|
811 | ;==================================================== |
---|
812 | ; def of gphiu: defined as the middle of T(i,j) T(i+1,j) |
---|
813 | ;==================================================== |
---|
814 | ; |
---|
815 | IF jpi GT 1 THEN BEGIN |
---|
816 | ; we must compute stepyu: y distance between T(i+1,j) T(i,j) |
---|
817 | IF keyword_set(key_irregular) THEN BEGIN |
---|
818 | stepyu = shift(gphit, -1, 0) - gphit |
---|
819 | IF NOT keyword_set(key_periodic) THEN $ |
---|
820 | stepyu[jpi-1, *] = stepyu[jpi-2, *] |
---|
821 | gphiu = gphit + 0.5 * stepyu |
---|
822 | ENDIF ELSE gphiu = gphit |
---|
823 | ENDIF ELSE gphiu = gphit |
---|
824 | IF key_onearth THEN gphiu = -90. > gphiu < 90. |
---|
825 | ; |
---|
826 | IF jpj EQ 1 THEN gphiu = reform(gphiu, jpi, jpj, /over) |
---|
827 | ; |
---|
828 | ;==================================================== |
---|
829 | ; def of glamv: defined as the middle of T(i,j) T(i,j+1) |
---|
830 | ;==================================================== |
---|
831 | ; |
---|
832 | IF jpj GT 1 THEN BEGIN |
---|
833 | ; we must compute stepxv: x distance between T(i,j) T(i,j+1) |
---|
834 | IF keyword_set(irregular) THEN BEGIN |
---|
835 | IF keyword_set(key_periodic) THEN BEGIN |
---|
836 | stepxv = (glamt + 720) MOD 360 |
---|
837 | stepxv = shift(stepxv, 0, -1) - stepxv |
---|
838 | stepxv = [ [[stepxv]], [[stepxv + 360]], [[stepxv - 360]] ] |
---|
839 | stepxv = min(abs(stepxv), dimension = 3) |
---|
840 | ENDIF ELSE stepxv = shift(glamt, 0, -1) - glamt |
---|
841 | stepxv[*, jpj-1] = stepxv[*, jpj-2] |
---|
842 | glamv = glamt + 0.5 * stepxv |
---|
843 | IF keyword_set(glamboundary) AND key_onearth THEN $ |
---|
844 | glamv = glamboundary[0] > temporary(glamv) < glamboundary[1] |
---|
845 | ENDIF ELSE glamv = glamt |
---|
846 | ENDIF ELSE glamv = glamt |
---|
847 | ; |
---|
848 | ;==================================================== |
---|
849 | ; def of gphiv: defined as the middle of T(i,j) T(i,j+1) |
---|
850 | ;==================================================== |
---|
851 | ; |
---|
852 | IF keyword_set(key_irregular) THEN $ |
---|
853 | gphiv = gphit + 0.5 * stepyv $ |
---|
854 | ELSE gphiv = gphif |
---|
855 | IF key_onearth THEN gphiv = -90. > gphiv < 90. |
---|
856 | ; |
---|
857 | IF jpj EQ 1 THEN gphiv = reform(gphiv, jpi, jpj, /over) |
---|
858 | ; |
---|
859 | ;==================================================== |
---|
860 | ; e1u: x distance between T(i,j) and T(i+1,j) |
---|
861 | ;==================================================== |
---|
862 | ; |
---|
863 | IF jpi GT 1 AND n_elements(stepx) NE 1 THEN $ |
---|
864 | e1u = stepxu ELSE e1u = e1t |
---|
865 | ; |
---|
866 | ;==================================================== |
---|
867 | ; e2u: y distance between F(i,j-1) and F(i,j) |
---|
868 | ;==================================================== |
---|
869 | ; |
---|
870 | IF keyword_set(key_irregular) THEN BEGIN |
---|
871 | e2u = gphif - shift(gphif, 0, 1) |
---|
872 | e2u[*, 0] = e2u[*, 1] |
---|
873 | IF key_onearth THEN e2u = r * !pi/180. * temporary(e2u) |
---|
874 | ENDIF ELSE e2u = e2t |
---|
875 | ; |
---|
876 | IF jpj EQ 1 THEN e2u = reform(e2u, jpi, jpj, /over) |
---|
877 | ; |
---|
878 | ;==================================================== |
---|
879 | ; e1v: x distance between F(i-1,j) and F(i,j) |
---|
880 | ;==================================================== |
---|
881 | ; |
---|
882 | IF keyword_set(irregular) THEN BEGIN |
---|
883 | IF keyword_set(key_periodic) THEN BEGIN |
---|
884 | e1v = (glamf + 720) MOD 360 |
---|
885 | e1v = e1v - shift(e1v, 1, 0) |
---|
886 | e1v = [ [[e1v]], [[e1v + 360]], [[e1v - 360]] ] |
---|
887 | e1v = min(abs(e1v), dimension = 3) |
---|
888 | ENDIF ELSE BEGIN |
---|
889 | e1v = glamf - shift(glamf, 1, 0) |
---|
890 | e1v[0, *] = stepxf[1, *] |
---|
891 | ENDELSE |
---|
892 | ENDIF ELSE e1v = e1t |
---|
893 | ; |
---|
894 | IF jpj EQ 1 THEN e1v = reform(e1v, jpi, jpj, /over) |
---|
895 | ; |
---|
896 | ;==================================================== |
---|
897 | ; e2v: y distance between T(i,j) and T(i+1,j) |
---|
898 | ;==================================================== |
---|
899 | ; |
---|
900 | IF jpj GT 1 and n_elements(stepy) NE 1 THEN BEGIN |
---|
901 | e2v = stepyv |
---|
902 | IF key_onearth THEN e2v = r * !pi/180. * temporary(e2v) |
---|
903 | ENDIF ELSE e2v = e2t |
---|
904 | ; |
---|
905 | ;==================================================== |
---|
906 | ; e1f: x distance between V(i,j) and V(i+1,j) |
---|
907 | ;==================================================== |
---|
908 | ; |
---|
909 | IF keyword_set(irregular) THEN BEGIN |
---|
910 | IF keyword_set(key_periodic) THEN BEGIN |
---|
911 | e1f = (glamv + 720) MOD 360 |
---|
912 | e1f = shift(e1f, -1, 0) - e1f |
---|
913 | e1f = [ [[e1f]], [[e1f + 360]], [[e1f - 360]] ] |
---|
914 | e1f = min(abs(e1f), dimension = 3) |
---|
915 | ENDIF ELSE BEGIN |
---|
916 | e1f = shift(glamv, -1, 0) - glamt |
---|
917 | e1f[jpi-1, *] = stepxf[jpi-2, *] |
---|
918 | ENDELSE |
---|
919 | ENDIF ELSE e1f = e1u |
---|
920 | ; |
---|
921 | IF jpj EQ 1 THEN e1f = reform(e1f, jpi, jpj, /over) |
---|
922 | ; |
---|
923 | ;==================================================== |
---|
924 | ; e2f: y distance between U(i,j) and U(i,j+1) |
---|
925 | ;==================================================== |
---|
926 | ; |
---|
927 | IF keyword_set(key_irregular) THEN BEGIN |
---|
928 | e2f = shift(gphiu, 0, -1) - gphiu |
---|
929 | e2f[*, jpj-1] = e2f[*, jpj-2] |
---|
930 | IF key_onearth THEN e2f = r * !pi/180. * temporary(e2f) |
---|
931 | ENDIF ELSE e2f = e2v |
---|
932 | ; |
---|
933 | IF jpj EQ 1 THEN e2f = reform(e2f, jpi, jpj, /over) |
---|
934 | ; |
---|
935 | ENDIF |
---|
936 | ;~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
937 | ; |
---|
938 | ; |
---|
939 | ;==================================================== |
---|
940 | ; e1[tuvf] from degree to meters |
---|
941 | ;==================================================== |
---|
942 | ; |
---|
943 | IF keyword_set(key_onearth) THEN BEGIN |
---|
944 | e1t = r * !pi/180. * temporary(e1t) * cos(!pi/180.*gphit) |
---|
945 | IF keyword_set(fullcgrid) THEN BEGIN |
---|
946 | e1u = r * !pi/180. * temporary(e1u) * cos(!pi/180.*gphiu) |
---|
947 | e1v = r * !pi/180. * temporary(e1v) * cos(!pi/180.*gphiv) |
---|
948 | e1f = r * !pi/180. * temporary(e1f) * cos(!pi/180.*gphif) |
---|
949 | ENDIF |
---|
950 | ENDIF |
---|
951 | ; |
---|
952 | IF jpj EQ 1 THEN BEGIN |
---|
953 | e1t = reform(e1t, jpi, jpj, /over) |
---|
954 | IF keyword_set(fullcgrid) THEN BEGIN |
---|
955 | e1u = reform(e1u, jpi, jpj, /over) |
---|
956 | e1v = reform(e1v, jpi, jpj, /over) |
---|
957 | e1f = reform(e1f, jpi, jpj, /over) |
---|
958 | ENDIF |
---|
959 | ENDIF |
---|
960 | ; |
---|
961 | ;==================================================== |
---|
962 | ; if not fullcgrid: make sure we don't use glam[uv], gphi[uv], e[12][uvf] |
---|
963 | ;==================================================== |
---|
964 | ; |
---|
965 | IF NOT keyword_set(fullcgrid) THEN BEGIN |
---|
966 | glamu = !values.f_nan & glamv = !values.f_nan |
---|
967 | gphiu = !values.f_nan & gphiv = !values.f_nan |
---|
968 | e1u = !values.f_nan & e1v = !values.f_nan & e1f = !values.f_nan |
---|
969 | e2u = !values.f_nan & e2v = !values.f_nan & e2f = !values.f_nan |
---|
970 | firstxu = !values.f_nan & lastxu = !values.f_nan & nxu = !values.f_nan |
---|
971 | firstyu = !values.f_nan & lastyu = !values.f_nan & nyu = !values.f_nan |
---|
972 | firstxv = !values.f_nan & lastxv = !values.f_nan & nxv = !values.f_nan |
---|
973 | firstyv = !values.f_nan & lastyv = !values.f_nan & nyv = !values.f_nan |
---|
974 | ENDIF |
---|
975 | ; |
---|
976 | ;==================================================== |
---|
977 | ; Z direction |
---|
978 | ;==================================================== |
---|
979 | ; |
---|
980 | ; z axis |
---|
981 | ; |
---|
982 | CASE n_elements(zaxis) OF |
---|
983 | 0:BEGIN |
---|
984 | gdept = 0. |
---|
985 | key_zreverse = 0 |
---|
986 | END |
---|
987 | 1:BEGIN |
---|
988 | gdept = zaxis |
---|
989 | key_zreverse = 0 |
---|
990 | END |
---|
991 | ELSE:BEGIN |
---|
992 | gdept = zaxis[izminmesh:izmaxmesh] |
---|
993 | IF n_elements(zreverse) EQ 0 THEN BEGIN |
---|
994 | IF jpk GT 1 THEN BEGIN |
---|
995 | if gdept[0] GT gdept[1] then key_zreverse = 1 ELSE key_zreverse = 0 |
---|
996 | ENDIF ELSE key_zreverse = 0 |
---|
997 | ENDIF ELSE key_zreverse = zreverse |
---|
998 | IF keyword_set(key_zreverse) THEN gdept = reverse(gdept) |
---|
999 | END |
---|
1000 | ENDCASE |
---|
1001 | ; |
---|
1002 | if n_elements(gdept) GT 1 then BEGIN |
---|
1003 | stepz = shift(gdept, -1)-gdept |
---|
1004 | stepz[jpk-1] = stepz[jpk-2] |
---|
1005 | gdepw = 0. > (gdept-stepz/2.) |
---|
1006 | ENDIF ELSE BEGIN |
---|
1007 | stepz = 1. |
---|
1008 | gdepw = gdept |
---|
1009 | ENDELSE |
---|
1010 | IF keyword_set(romsh) THEN gdepw = gdept |
---|
1011 | ; |
---|
1012 | ;==================================================== |
---|
1013 | ; e3[tw]: |
---|
1014 | ;==================================================== |
---|
1015 | ; |
---|
1016 | e3t = stepz |
---|
1017 | IF n_elements(stepz) GT 1 THEN BEGIN |
---|
1018 | e3w = 0.5*(stepz+shift(stepz, 1)) |
---|
1019 | e3w[0] = 0.5*e3t[0] |
---|
1020 | ENDIF ELSE e3w = e3t |
---|
1021 | ; |
---|
1022 | ;==================================================== |
---|
1023 | ; Mask |
---|
1024 | ;==================================================== |
---|
1025 | ; |
---|
1026 | ; default mask eq 1 |
---|
1027 | if NOT keyword_set(mask) then tmask = -1 ELSE tmask = mask |
---|
1028 | ; |
---|
1029 | if tmask[0] NE -1 then BEGIN |
---|
1030 | IF keyword_set(romsh) THEN tmask = tmask[0:jpiglo-1, 0:jpjglo-1] |
---|
1031 | IF n_elements(mask) EQ jpiglo*jpjglo AND jpkglo GT 1 THEN BEGIN |
---|
1032 | tmask = tmask[*]#replicate(1, jpkglo) |
---|
1033 | tmask = reform(tmask, jpiglo, jpjglo, jpkglo, /overwrite) |
---|
1034 | ENDIF |
---|
1035 | IF jpiglo EQ 1 OR jpjglo EQ 1 THEN tmask = reform(tmask, jpiglo, jpjglo, jpkglo, /overwrite) |
---|
1036 | tmask = byte(tmask[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh, izminmesh:izmaxmesh]) |
---|
1037 | tmask = reform(tmask, jpi, jpj, jpk, /over) |
---|
1038 | if key_shift NE 0 then tmask = shift(tmask, key_shift, 0, 0) |
---|
1039 | ; because tmask = reverse(tmask, 2) is not working if the 3rd |
---|
1040 | ; dimension of tmask = 1, we call reform. |
---|
1041 | IF jpk EQ 1 THEN tmask = reform(tmask, /over) |
---|
1042 | IF key_yreverse EQ 1 THEN tmask = reverse(tmask, 2) |
---|
1043 | IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over) |
---|
1044 | IF key_zreverse EQ 1 THEN tmask = reverse(tmask, 3) |
---|
1045 | IF jpk EQ 1 THEN tmask = reform(tmask, jpi, jpj, jpk, /over) |
---|
1046 | IF keyword_set(fullcgrid) THEN BEGIN |
---|
1047 | IF keyword_set(key_periodic) THEN BEGIN |
---|
1048 | msk = tmask*shift(tmask, -1, 0, 0) |
---|
1049 | umaskred = msk[jpi-1, *, *] |
---|
1050 | ENDIF ELSE umaskred = tmask[jpi-1, *, *] |
---|
1051 | vmaskred = tmask[*, jpj-1, *] |
---|
1052 | ENDIF |
---|
1053 | fmaskredy = tmask[jpi-1, *, *] |
---|
1054 | fmaskredx = tmask[*, jpj-1, *] |
---|
1055 | ENDIF ELSE BEGIN |
---|
1056 | tmask = replicate(1b, jpi, jpj, jpk) |
---|
1057 | IF keyword_set(fullcgrid) THEN BEGIN |
---|
1058 | umaskred = replicate(1b, jpj, jpk) |
---|
1059 | vmaskred = replicate(1b, jpi, jpk) |
---|
1060 | ENDIF |
---|
1061 | fmaskredy = replicate(1b, jpj, jpk) |
---|
1062 | fmaskredx = replicate(1b, jpi, jpk) |
---|
1063 | ENDELSE |
---|
1064 | ; |
---|
1065 | IF jpi GT 2 AND jpj GT 2 AND NOT keyword_set(plain) $ |
---|
1066 | AND ixminmesh EQ 0l AND ixmaxmesh eq jpiglo-1 $ |
---|
1067 | AND iyminmesh EQ 0l AND iymaxmesh eq jpjglo-1 $ |
---|
1068 | AND total(tmask[*, 0, *]) EQ 0 AND total(tmask[*, jpj-1, *]) EQ 0 $ |
---|
1069 | AND total(tmask[0, *, *]) EQ 0 AND total(tmask[jpi-1, *, *]) EQ 0 THEN BEGIN |
---|
1070 | xminmesh = 1 |
---|
1071 | xmaxmesh = -1 |
---|
1072 | yminmesh = 1 |
---|
1073 | ymaxmesh = -1 |
---|
1074 | computegrid, XAXIS = glamt, YAXIS = gphit, ZAXIS = zaxis $ |
---|
1075 | , MASK = mask, GLAMBOUNDARY = glamboundary $ |
---|
1076 | , XMINMESH = xminmesh, XMAXMESH = xmaxmesh $ |
---|
1077 | , YMINMESH = yminmesh, YMAXMESH = ymaxmesh $ |
---|
1078 | , ZMINMESH = zminmesh, ZMAXMESH = zmaxmesh $ |
---|
1079 | , ONEARTH = onearth, PERIODIC = periodic $ |
---|
1080 | , PLAIN = plain, SHIFT = shift, STRIDE = stride $ |
---|
1081 | , FULLCGRID = fullcgrid, XYINDEX = xyindex $ |
---|
1082 | , UBASE2TBASE = ubase2tbase, VBASE2TBASE = vbase2tbase $ |
---|
1083 | , FBASE2TBASE = fbase2tbase, STRCALLING = strcalling $ |
---|
1084 | , ROMSH = romsh, _extra = ex |
---|
1085 | return |
---|
1086 | ENDIF |
---|
1087 | ; |
---|
1088 | IF NOT keyword_set(fullcgrid) THEN BEGIN |
---|
1089 | umaskred = !values.f_nan |
---|
1090 | vmaskred = !values.f_nan |
---|
1091 | ENDIF |
---|
1092 | ; |
---|
1093 | ;==================================================== |
---|
1094 | ; stride... |
---|
1095 | ;==================================================== |
---|
1096 | ; |
---|
1097 | IF total(key_stride) GT 3 THEN BEGIN |
---|
1098 | IF key_shift NE 0 THEN BEGIN |
---|
1099 | ; for explanation, see header of read_ncdf_varget.pro |
---|
1100 | jpiright = key_shift |
---|
1101 | jpileft = jpi - key_shift - ( (key_stride[0]-1)-((key_shift-1) MOD key_stride[0]) ) |
---|
1102 | jpi = ((jpiright-1)/key_stride[0]+1) + ((jpileft-1)/key_stride[0]+1) |
---|
1103 | ENDIF ELSE jpi = (jpi-1)/key_stride[0]+1 |
---|
1104 | jpj = (jpj-1)/key_stride[1]+1 |
---|
1105 | jpk = (jpk-1)/key_stride[2]+1 |
---|
1106 | ; |
---|
1107 | glamt = (temporary(glamt))[0:*:stride[0], 0:*:stride[1]] |
---|
1108 | gphit = (temporary(gphit))[0:*:stride[0], 0:*:stride[1]] |
---|
1109 | e1t = (temporary(e1t))[0:*:stride[0], 0:*:stride[1]] |
---|
1110 | e2t = (temporary(e2t))[0:*:stride[0], 0:*:stride[1]] |
---|
1111 | tmask = (temporary(tmask))[0:*:stride[0], 0:*:stride[1], 0:*:stride[2]] |
---|
1112 | gdept = gdept[0:*:stride[2]] |
---|
1113 | gdepw = gdepw[0:*:stride[2]] |
---|
1114 | e3t = e3t[0:*:stride[2]] |
---|
1115 | e3w = e3w[0:*:stride[2]] |
---|
1116 | ; we must recompute glamf and gphif... |
---|
1117 | IF jpi GT 1 THEN BEGIN |
---|
1118 | if (keyword_set(key_onearth) AND keyword_set(xnotsorted)) $ |
---|
1119 | OR (keyword_set(key_periodic) AND key_irregular) then BEGIN |
---|
1120 | stepxf = (glamt + 720) MOD 360 |
---|
1121 | stepxf = shift(stepxf, -1, -1) - stepxf |
---|
1122 | stepxf = [ [[stepxf]], [[stepxf + 360]], [[stepxf - 360]] ] |
---|
1123 | stepxf = min(abs(stepxf), dimension = 3) |
---|
1124 | IF NOT keyword_set(key_periodic) THEN $ |
---|
1125 | stepxf[jpi-1, *] = stepxf[jpi-2, *] |
---|
1126 | ENDIF ELSE BEGIN |
---|
1127 | stepxf = shift(glamt, -1, -1) - glamt |
---|
1128 | IF keyword_set(key_periodic) THEN $ |
---|
1129 | stepxf[jpi-1, *] = 360 + stepxf[jpi-1, *] $ |
---|
1130 | ELSE stepxf[jpi-1, *] = stepxf[jpi-2, *] |
---|
1131 | ENDELSE |
---|
1132 | IF jpj GT 1 THEN BEGIN |
---|
1133 | stepxf[*, jpj-1] = stepxf[*, jpj-2] |
---|
1134 | stepxf[jpi-1, jpj-1] = stepxf[jpi-2, jpj-2] |
---|
1135 | ENDIF |
---|
1136 | glamf = glamt + 0.5 * stepxf |
---|
1137 | ENDIF ELSE glamf = glamt + 0.5 |
---|
1138 | IF keyword_set(key_periodic) AND (max(glamf)-min(glamt)) GE 360 THEN BEGIN |
---|
1139 | IF NOT keyword_set(glamboundary) THEN BEGIN |
---|
1140 | bigger = where(glamf GE min(glamt)+360) |
---|
1141 | glamf[bigger] = glamf[bigger]-360. |
---|
1142 | ENDIF ELSE glamf = glamboundary[0] > temporary(glamf) < glamboundary[1] |
---|
1143 | ENDIF |
---|
1144 | IF jpj GT 1 THEN BEGIN |
---|
1145 | ; we must compute stepyf: y distance between T(i,j) T(i+1,j+1) |
---|
1146 | stepyf = shift(gphit, -1, -1) - gphit |
---|
1147 | stepyf[*, jpj-1] = stepyf[*, jpj-2] |
---|
1148 | IF jpi GT 1 THEN BEGIN |
---|
1149 | if NOT keyword_set(key_periodic) THEN $ |
---|
1150 | stepyf[jpi-1, *] = stepyf[jpi-2, *] |
---|
1151 | stepyf[jpi-1, jpj-1] = stepyf[jpi-2, jpj-2] |
---|
1152 | ENDIF |
---|
1153 | gphif = gphit + 0.5 * stepyf |
---|
1154 | ENDIF ELSE gphif = gphit + 0.5 |
---|
1155 | ; |
---|
1156 | IF jpj EQ 1 THEN BEGIN |
---|
1157 | glamt = reform(glamt, jpi, jpj, /over) |
---|
1158 | gphit = reform(gphit, jpi, jpj, /over) |
---|
1159 | glamf = reform(glamf, jpi, jpj, /over) |
---|
1160 | gphif = reform(gphif, jpi, jpj, /over) |
---|
1161 | e1t = reform(e1t, jpi, jpj, /over) |
---|
1162 | e2t = reform(e2t, jpi, jpj, /over) |
---|
1163 | ENDIF |
---|
1164 | ; |
---|
1165 | IF keyword_set(fullcgrid) THEN BEGIN |
---|
1166 | glamu = (temporary(glamu))[0:*:stride[0], 0:*:stride[1]] |
---|
1167 | gphiu = (temporary(gphiu))[0:*:stride[0], 0:*:stride[1]] |
---|
1168 | e1u = (temporary(e1u))[0:*:stride[0], 0:*:stride[1]] |
---|
1169 | e2u = (temporary(e2u))[0:*:stride[0], 0:*:stride[1]] |
---|
1170 | glamv = (temporary(glamv))[0:*:stride[0], 0:*:stride[1]] |
---|
1171 | gphiv = (temporary(gphiv))[0:*:stride[0], 0:*:stride[1]] |
---|
1172 | e1v = (temporary(e1v))[0:*:stride[0], 0:*:stride[1]] |
---|
1173 | e2v = (temporary(e2v))[0:*:stride[0], 0:*:stride[1]] |
---|
1174 | e1f = (temporary(e1f))[0:*:stride[0], 0:*:stride[1]] |
---|
1175 | e2f = (temporary(e2f))[0:*:stride[0], 0:*:stride[1]] |
---|
1176 | umaskred = (temporary(umaskred))[0, 0:*:stride[1], 0:*:stride[2]] |
---|
1177 | vmaskred = (temporary(vmaskred))[0:*:stride[0], 0, 0:*:stride[2]] |
---|
1178 | fmaskredy = (temporary(fmaskredy))[0, 0:*:stride[1], 0:*:stride[2]] |
---|
1179 | fmaskredx = (temporary(fmaskredx))[0:*:stride[0], 0, 0:*:stride[2]] |
---|
1180 | IF jpj EQ 1 THEN BEGIN |
---|
1181 | glamu = reform(glamu, jpi, jpj, /over) |
---|
1182 | gphiu = reform(gphiu, jpi, jpj, /over) |
---|
1183 | e1u = reform(e1u, jpi, jpj, /over) |
---|
1184 | e2u = reform(e2u, jpi, jpj, /over) |
---|
1185 | glamv = reform(glamv, jpi, jpj, /over) |
---|
1186 | gphiv = reform(gphiv, jpi, jpj, /over) |
---|
1187 | e1v = reform(e1v, jpi, jpj, /over) |
---|
1188 | e2v = reform(e2v, jpi, jpj, /over) |
---|
1189 | e1f = reform(e1f, jpi, jpj, /over) |
---|
1190 | e2f = reform(e2f, jpi, jpj, /over) |
---|
1191 | ENDIF |
---|
1192 | ENDIF |
---|
1193 | ENDIF |
---|
1194 | ; |
---|
1195 | ;==================================================== |
---|
1196 | ; apply all the grid parameters |
---|
1197 | ;==================================================== |
---|
1198 | ; |
---|
1199 | @updateold |
---|
1200 | domdef |
---|
1201 | ; |
---|
1202 | ;==================================================== |
---|
1203 | ; Triangulation |
---|
1204 | ;==================================================== |
---|
1205 | ; |
---|
1206 | IF total(tmask) EQ jpi*jpj*jpk $ |
---|
1207 | AND NOT keyword_set(key_irregular) THEN triangles_list = -1 $ |
---|
1208 | ELSE BEGIN |
---|
1209 | ; are we using ORCA2 ? |
---|
1210 | IF jpiglo EQ 182 AND jpi EQ 180 AND jpjglo EQ 149 AND jpj EQ 148 THEN $ |
---|
1211 | triangles_list = triangule() ELSE triangles_list = triangule(/keep_cont) |
---|
1212 | ENDELSE |
---|
1213 | ; |
---|
1214 | ;==================================================== |
---|
1215 | ; time axis (default definition) |
---|
1216 | ;==================================================== |
---|
1217 | ; |
---|
1218 | IF n_elements(time) EQ 0 OR n_elements(jpt) EQ 0 THEN BEGIN |
---|
1219 | jpt = 1 |
---|
1220 | time = 0 |
---|
1221 | ENDIF |
---|
1222 | ; |
---|
1223 | IF NOT keyword_set(key_forgetold) THEN BEGIN |
---|
1224 | @updateold |
---|
1225 | ENDIF |
---|
1226 | ;==================================================== |
---|
1227 | ; grid parameters used by xxx |
---|
1228 | ;==================================================== |
---|
1229 | ; |
---|
1230 | IF NOT keyword_set(strcalling) THEN BEGIN |
---|
1231 | IF n_elements(ccmeshparameters) EQ 0 THEN strcalling = 'computegrid' $ |
---|
1232 | ELSE strcalling = ccmeshparameters.filename |
---|
1233 | ENDIF |
---|
1234 | IF n_elements(glamt) GE 2 THEN BEGIN |
---|
1235 | glaminfo = moment(glamt) |
---|
1236 | IF finite(glaminfo[2]) EQ 0 THEN glaminfo = glaminfo[0:1] |
---|
1237 | gphiinfo = moment(gphit) |
---|
1238 | IF finite(gphiinfo[2]) EQ 0 THEN gphiinfo = gphiinfo[0:1] |
---|
1239 | ENDIF ELSE BEGIN |
---|
1240 | glaminfo = glamt |
---|
1241 | gphiinfo = gphit |
---|
1242 | ENDELSE |
---|
1243 | IF keyword_set(romsh) THEN $ |
---|
1244 | romszinfos = {h:romsh[ixminmesh:ixmaxmesh, iyminmesh:iymaxmesh], zeta:replicate(0., jpi, jpj), theta_s:-1, theta_b:-1, hc:-1} $ |
---|
1245 | ELSE romszinfos = {h:-1, zeta:-1, theta_s:-1, theta_b:-1, hc:-1} |
---|
1246 | |
---|
1247 | ccmeshparameters = {filename:strcalling $ |
---|
1248 | , glaminfo:float(string(glaminfo, format = '(E11.4)')) $ |
---|
1249 | , gphiinfo:float(string(gphiinfo, format = '(E11.4)')) $ |
---|
1250 | , jpiglo:jpiglo, jpjglo:jpjglo, jpkglo:jpkglo $ |
---|
1251 | , jpi:jpi, jpj:jpj, jpk:jpk $ |
---|
1252 | , ixminmesh:ixminmesh, ixmaxmesh:ixmaxmesh $ |
---|
1253 | , iyminmesh:iyminmesh, iymaxmesh:iymaxmesh $ |
---|
1254 | , izminmesh:izminmesh, izmaxmesh:izmaxmesh $ |
---|
1255 | , key_shift:key_shift, key_periodic:key_periodic $ |
---|
1256 | , key_stride:key_stride, key_gridtype:key_gridtype $ |
---|
1257 | , key_yreverse:key_yreverse, key_zreverse:key_zreverse $ |
---|
1258 | , key_partialstep:key_partialstep, key_onearth:key_onearth} |
---|
1259 | |
---|
1260 | ccreadparameters = {funclec_name:'read_ncdf' $ |
---|
1261 | , jpidta:jpidta, jpjdta:jpjdta, jpkdta:jpkdta $ |
---|
1262 | , ixmindta:ixmindta, ixmaxdta:ixmaxdta $ |
---|
1263 | , iymindta:iymindta, iymaxdta:iymaxdta $ |
---|
1264 | , izmindta:izmindta, izmaxdta:izmaxdta} |
---|
1265 | ;------------------------------------------------------------ |
---|
1266 | IF keyword_set(key_performance) EQ 1 THEN $ |
---|
1267 | print, 'time computegrid', systime(1)-time1 |
---|
1268 | ;------------------------------------------------------------ |
---|
1269 | return |
---|
1270 | end |
---|