source: trunk/ToBeReviewed/HOPE/computehopegrid.pro @ 53

Last change on this file since 53 was 53, checked in by pinsard, 18 years ago

upgrade of ToBeReviewed?/HOPE according to cerbere.lodyc.jussieu.fr: /usr/home/smasson/SAXO_RD/

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