source: trunk/SRC/ToBeReviewed/HOPE/computehopegrid.pro @ 231

Last change on this file since 231 was 231, checked in by pinsard, 17 years ago

improvements/corrections of some *.pro headers
  • Property svn:keywords set to Id
File size: 6.3 KB
Line 
1;+
2;
3; @file_comments
4;
5; @categories
6; Grid
7;
8; @param XAXIS
9;
10; @param YAXIS
11;
12; @param ZAXIS
13;
14; @param LINETYPE
15;
16; @keyword FORTHEMASK
17;
18; @keyword WPOINT
19;
20; @keyword FIRSTS
21;
22; @keyword LASTS
23;
24; @keyword PTTYPE
25;
26; @uses
27; common.pro
28;
29; @history
30; Sebastien Masson (smasson\@lodyc.jussieu.fr)
31;                      2001-06
32;
33; @version
34; $Id$
35;
36; @todo seb
37;
38;-
39;
40PRO computehopegrid, xaxis, yaxis, zaxis, linetype, FORTHEMASK = forthemask, WPOINT = wpoint, FIRSTS = firsts, LASTS = lasts, PTTYPE = pttype
41;---------------------------------------------------------
42;
43  compile_opt idl2, strictarrsubs
44;
45@cm_4mesh
46@cm_4data
47  IF NOT keyword_set(key_forgetold) THEN BEGIN
48@updatenew
49  ENDIF
50;---------------------------------------------------------
51;------------------------------------------------------------
52   if not (keyword_set(scalar)*keyword_set(vector)) then scalar=1
53;------------------------------------------------------------
54   jpiglo = n_elements(xaxis)
55   jpjglo = n_elements(yaxis)
56   jpkglo = n_elements(zaxis)
57;
58   jpi = jpiglo
59   jpj = jpjglo
60   jpk = jpkglo
61;
62   if NOT keyword_set(firsts) then firsts = [0, 0, 0]
63   if NOT keyword_set(lasts) then lasts = [jpi-1, jpj-1, jpk-1]
64;
65; determination of the grid type and of the point type
66;
67   if keyword_set(pttype) then vargrid = pttype
68   if linetype EQ 'odd-even' then key_gridtype = 'e' ELSE key_gridtype = 'c'
69;
70; computation of the horizontal grid
71;
72   if key_gridtype EQ 'e' then begin
73      if vargrid EQ 'T' then begin
74         glamt=xaxis
75         glamt = glamt#replicate(1, jpj)
76         xstep=(glamt[1,0]-glamt[0,0])/2.
77         glamt[*,2*lindgen((jpj)/2)]=glamt[*,2*lindgen(jpj/2)]-xstep
78         glamu=glamt+xstep
79      ENDIF ELSE BEGIN
80         glamu=xaxis
81         glamu = glamu#replicate(1, jpj)
82         xstep=(glamu[1,0]-glamu[0,0])/2.
83         glamu[*,2*lindgen((jpj)/2)]=glamu[*,2*lindgen(jpj/2)]-xstep
84         glamt=glamu-xstep
85      ENDELSE
86;zoom
87      glamt = glamt[firsts[0]:lasts[0], firsts[1]:lasts[1]]
88      glamu = glamu[firsts[0]:lasts[0], firsts[1]:lasts[1]]
89      jpiglo = lasts[0]-firsts[0]+1
90      jpi = jpiglo
91      jpjglo = lasts[1]-firsts[1]+1
92      jpj = jpjglo
93      glamv = glamu
94      glamf = glamu
95      gphit = yaxis[firsts[1]:lasts[1]]
96      gphit = replicate(1, jpi)#gphit
97      gphif = gphit
98      gphiu = gphit
99      gphiv = gphif
100   ENDIF ELSE BEGIN
101      if vargrid eq 'T' then begin
102         glamt=xaxis
103         glamt = glamt#replicate(1, jpj)
104         glamu=glamt+(glamt[1,0]-glamt[0,0])/2.
105      ENDIF ELSE BEGIN
106         glamu=xaxis
107         glamu = glamu#replicate(1, jpj)
108         xstep=(glamu[1,0]-glamu[0,0])/2.
109         glamt=glamu-(glamu[1,0]-glamu[0,0])/2.
110      ENDELSE
111;zoom
112      glamt = glamt[firsts[0]:lasts[0], firsts[1]:lasts[1]]
113      glamu = glamu[firsts[0]:lasts[0], firsts[1]:lasts[1]]
114      jpiglo = lasts[0]-firsts[0]+1
115      jpi = jpiglo
116      jpjglo = lasts[1]-firsts[1]+1
117      jpj = jpjglo
118      glamv = glamt
119      glamf = glamu
120      gphit = yaxis[firsts[1]:lasts[1]]
121      gphit = replicate(1, jpi)#gphit
122      gphiu = gphit
123      if jpj GT 1 then begin
124         gphiv = gphit+(gphit[0, 1]-gphit[0, 0])/2.
125         gphif = gphit+(gphit[0, 1]-gphit[0, 0])/2.
126      ENDIF ELSE BEGIN
127         gphiv = gphit
128         gphif = gphit
129      ENDELSE
130   ENDELSE
131;
132; computation of the vertical grid
133;
134   if keyword_set(wpoint) then begin
135      gdepw = zaxis
136      if jpk ne 1 then begin
137         gdept=(shift(gdepw, -1)+gdepw)/2.
138         gdept[jpk-1]=gdepw[jpk-1]+.5*(gdepw[jpk-1]-gdepw[jpk-2])
139      endif else gdept=zaxis
140   endif else begin
141      gdept = zaxis
142      if jpk ne 1 then begin
143         gdepw=(shift(gdept, 1)+gdept)/2.
144         gdepw[0]=0
145      endif else gdepw = zaxis
146   endelse
147;
148; computation of the vertical scale factors
149;
150   if jpk ne 1 then begin
151      e3t = abs(shift(gdepw,-1)-gdepw)
152      e3t[jpk-1] = abs(gdept[jpk-1]-gdepw[jpk-1])
153      e3w = abs(gdept-shift(gdept,1))
154      e3w[0] = abs(gdept[0]-gdepw[0])
155   endif else begin
156      e3t=1
157      e3w=1
158   endelse
159; zoom
160   gdept = gdept[firsts[2]:lasts[2]]
161   gdepw = gdepw[firsts[2]:lasts[2]]
162   e3t = e3t[firsts[2]:lasts[2]]
163   e3w = e3w[firsts[2]:lasts[2]]
164   jpkglo = lasts[2]-firsts[2]+1
165   jpk = jpkglo
166;
167; computation of the horizontal scale factors
168;
169   e1t = replicate(1b,jpi,jpj)
170   e2t = replicate(1b,jpi,jpj)
171   e1u = e1t
172   e2u = e2t
173   e1v = e1t
174   e2v = e2t
175   e1f = e1t
176   e2f = e2t
177;
178; mask
179;
180   tmask = replicate(1b, jpi, jpj, jpk)
181   if keyword_set(forthemask) then BEGIN
182      land=where(forthemask ge valmask/10)
183      if land[0] ne -1 then tmask[land] = 0b
184   endif
185   umaskred = replicate(1, jpj, jpk)
186   vmaskred = replicate(1, jpi, jpk)
187   fmaskredy = replicate(1, jpj, jpk)
188   fmaskredx = replicate(1, jpi, jpk)
189   @updateold
190   domdef
191   if keyword_set(firsts) AND keyword_set(lasts) then BEGIN
192      if vargrid EQ 'T' then BEGIN
193         if jpj GT 1 then begin
194            lon1 = min(glamt[0, [0, 1]])
195            lon2 = max(glamt[jpi-1, [0, 1]])
196         endif ELSE BEGIN
197            lon1 = min(glamt[0, 0])
198            lon2 = max(glamt[jpi-1, [0]])
199         ENDELSE
200      ENDIF ELSE BEGIN
201         if jpj GT 1 then begin
202            lon1 = min(glamu[0, [0, 1]])
203            lon2 = max(glamu[jpi-1, [0, 1]])
204         endif ELSE BEGIN
205            lon1 = min(glamu[0, 0])
206            lon2 = max(glamu[jpi-1, 0])
207         ENDELSE
208      ENDELSE
209      lat1 = min([gphit[0, 0], gphit[0, jpj-1]])
210      lat2 = max([gphit[0, 0], gphit[0, jpj-1]])
211      domdef, lon1, lon2, lat1, lat2, gdepw[0], gdept[jpk-1], gridtype = vargrid
212   ENDIF
213;
214   ixminmesh  =0l
215   ixmaxmesh  =long(jpi-1)
216   iyminmesh  =0l
217   iymaxmesh  =long(jpj-1)
218   izminmesh  =0l
219   izmaxmesh  =long(jpk-1)
220;----------------------------------------------------------
221; for the triangulation
222;----------------------------------------------------------
223   key_periodic = glamt[0] EQ ((glamt[jpi-1]+glamt[+1]-glamt[0]) MOD 360)
224   if jpi gt 4 AND jpj GT 4 then begin
225      triangles_list = triangule(shifted = glamt[0, 0] LT glamt[0, 1])
226      twin_corners_up=-1
227      twin_corners_dn=-1
228   ENDIF ELSE BEGIN
229      triangles_list=-1
230      twin_corners_up=-1
231      twin_corners_dn=-1
232   ENDELSE
233;------------------------------------------------------------
234  IF NOT keyword_set(key_forgetold) THEN BEGIN
235   @updateold
236  ENDIF
237;------------------------------------------------------------
238   return
239end
240
241
Note: See TracBrowser for help on using the repository browser.