source: trunk/SRC/ToBeReviewed/HOPE/computehopegrid.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.6 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;
36  compile_opt idl2, strictarrsubs
37;
38@cm_4mesh
39@cm_4data
40  IF NOT keyword_set(key_forgetold) THEN BEGIN
41@updatenew
42  ENDIF
43;---------------------------------------------------------
44;------------------------------------------------------------
45   if not (keyword_set(scalar)*keyword_set(vector)) then scalar=1
46;------------------------------------------------------------
47   jpiglo = n_elements(xaxis)
48   jpjglo = n_elements(yaxis)
49   jpkglo = n_elements(zaxis)
50;
51   jpi = jpiglo
52   jpj = jpjglo
53   jpk = jpkglo
54;
55   if NOT keyword_set(firsts) then firsts = [0, 0, 0]
56   if NOT keyword_set(lasts) then lasts = [jpi-1, jpj-1, jpk-1]
57;
58; depermination of the grid type and of the point type
59;
60   if keyword_set(pttype) then vargrid = pttype
61   if linetype EQ 'odd-even' then key_gridtype = 'e' ELSE key_gridtype = 'c'
62;
63; computation of the horizontal grid
64;
65   if key_gridtype EQ 'e' then begin
66      if vargrid EQ 'T' then begin
67         glamt=xaxis
68         glamt = glamt#replicate(1, jpj)
69         xstep=(glamt[1,0]-glamt[0,0])/2.
70         glamt[*,2*lindgen((jpj)/2)]=glamt[*,2*lindgen(jpj/2)]-xstep
71         glamu=glamt+xstep
72      ENDIF ELSE BEGIN
73         glamu=xaxis
74         glamu = glamu#replicate(1, jpj)
75         xstep=(glamu[1,0]-glamu[0,0])/2.
76         glamu[*,2*lindgen((jpj)/2)]=glamu[*,2*lindgen(jpj/2)]-xstep
77         glamt=glamu-xstep
78      ENDELSE
79;zoom
80      glamt = glamt[firsts[0]:lasts[0], firsts[1]:lasts[1]]
81      glamu = glamu[firsts[0]:lasts[0], firsts[1]:lasts[1]]
82      jpiglo = lasts[0]-firsts[0]+1
83      jpi = jpiglo
84      jpjglo = lasts[1]-firsts[1]+1
85      jpj = jpjglo
86      glamv = glamu
87      glamf = glamu
88      gphit = yaxis[firsts[1]:lasts[1]]
89      gphit = replicate(1, jpi)#gphit
90      gphif = gphit
91      gphiu = gphit
92      gphiv = gphif
93   ENDIF ELSE BEGIN
94      if vargrid eq 'T' then begin
95         glamt=xaxis
96         glamt = glamt#replicate(1, jpj)
97         glamu=glamt+(glamt[1,0]-glamt[0,0])/2.
98      ENDIF ELSE BEGIN
99         glamu=xaxis
100         glamu = glamu#replicate(1, jpj)
101         xstep=(glamu[1,0]-glamu[0,0])/2.
102         glamt=glamu-(glamu[1,0]-glamu[0,0])/2.
103      ENDELSE
104;zoom
105      glamt = glamt[firsts[0]:lasts[0], firsts[1]:lasts[1]]
106      glamu = glamu[firsts[0]:lasts[0], firsts[1]:lasts[1]]
107      jpiglo = lasts[0]-firsts[0]+1
108      jpi = jpiglo
109      jpjglo = lasts[1]-firsts[1]+1
110      jpj = jpjglo
111      glamv = glamt
112      glamf = glamu
113      gphit = yaxis[firsts[1]:lasts[1]]
114      gphit = replicate(1, jpi)#gphit
115      gphiu = gphit
116      if jpj GT 1 then begin
117         gphiv = gphit+(gphit[0, 1]-gphit[0, 0])/2.
118         gphif = gphit+(gphit[0, 1]-gphit[0, 0])/2.
119      ENDIF ELSE BEGIN
120         gphiv = gphit
121         gphif = gphit
122      ENDELSE
123   ENDELSE
124;
125; computation of the vertical grid
126;
127   if keyword_set(wpoint) then begin
128      gdepw = zaxis
129      if jpk ne 1 then begin
130         gdept=(shift(gdepw, -1)+gdepw)/2.
131         gdept[jpk-1]=gdepw[jpk-1]+.5*(gdepw[jpk-1]-gdepw[jpk-2])
132      endif else gdept=zaxis
133   endif else begin
134      gdept = zaxis
135      if jpk ne 1 then begin
136         gdepw=(shift(gdept, 1)+gdept)/2.
137         gdepw[0]=0
138      endif else gdepw = zaxis
139   endelse
140;
141; computation of the vertical scale factors
142;
143   if jpk ne 1 then begin
144      e3t = abs(shift(gdepw,-1)-gdepw)
145      e3t[jpk-1] = abs(gdept[jpk-1]-gdepw[jpk-1])
146      e3w = abs(gdept-shift(gdept,1))
147      e3w[0] = abs(gdept[0]-gdepw[0])
148   endif else begin
149      e3t=1
150      e3w=1
151   endelse
152; zoom
153   gdept = gdept[firsts[2]:lasts[2]]
154   gdepw = gdepw[firsts[2]:lasts[2]]
155   e3t = e3t[firsts[2]:lasts[2]]
156   e3w = e3w[firsts[2]:lasts[2]]
157   jpkglo = lasts[2]-firsts[2]+1
158   jpk = jpkglo
159;
160; computation of the horizontal scale factors
161;
162   e1t = replicate(1b,jpi,jpj)
163   e2t = replicate(1b,jpi,jpj)
164   e1u = e1t
165   e2u = e2t
166   e1v = e1t
167   e2v = e2t
168   e1f = e1t
169   e2f = e2t
170;
171; mask
172;
173   tmask = replicate(1b, jpi, jpj, jpk)
174   if keyword_set(forthemask) then BEGIN
175      land=where(forthemask ge valmask/10)
176      if land[0] ne -1 then tmask[land] = 0b
177   endif
178   umaskred = replicate(1, jpj, jpk)
179   vmaskred = replicate(1, jpi, jpk)
180   fmaskredy = replicate(1, jpj, jpk)
181   fmaskredx = replicate(1, jpi, jpk)
182   @updateold
183   domdef
184   if keyword_set(firsts) AND keyword_set(lasts) then BEGIN
185      if vargrid EQ 'T' then BEGIN
186         if jpj GT 1 then begin
187            lon1 = min(glamt[0, [0, 1]])
188            lon2 = max(glamt[jpi-1, [0, 1]])
189         endif ELSE BEGIN
190            lon1 = min(glamt[0, 0])
191            lon2 = max(glamt[jpi-1, [0]])
192         ENDELSE
193      ENDIF ELSE BEGIN
194         if jpj GT 1 then begin
195            lon1 = min(glamu[0, [0, 1]])
196            lon2 = max(glamu[jpi-1, [0, 1]])
197         endif ELSE BEGIN
198            lon1 = min(glamu[0, 0])
199            lon2 = max(glamu[jpi-1, 0])
200         ENDELSE
201      ENDELSE
202      lat1 = min([gphit[0, 0], gphit[0, jpj-1]])
203      lat2 = max([gphit[0, 0], gphit[0, jpj-1]])
204      domdef, lon1, lon2, lat1, lat2, gdepw[0], gdept[jpk-1], gridtype = vargrid
205   ENDIF
206;
207   ixminmesh  =0l
208   ixmaxmesh  =long(jpi-1)
209   iyminmesh  =0l
210   iymaxmesh  =long(jpj-1)
211   izminmesh  =0l
212   izmaxmesh  =long(jpk-1)
213;----------------------------------------------------------
214; for the triangulation
215;----------------------------------------------------------
216   key_periodic = glamt[0] EQ ((glamt[jpi-1]+glamt[+1]-glamt[0]) MOD 360)
217   if jpi gt 4 AND jpj GT 4 then begin
218      triangles_list = triangule(shifted = glamt[0, 0] LT glamt[0, 1])
219      twin_corners_up=-1
220      twin_corners_dn=-1
221   ENDIF ELSE BEGIN
222      triangles_list=-1
223      twin_corners_up=-1
224      twin_corners_dn=-1
225   ENDELSE
226;------------------------------------------------------------
227  IF NOT keyword_set(key_forgetold) THEN BEGIN
228   @updateold
229  ENDIF
230;------------------------------------------------------------
231   return
232end
233
234
Note: See TracBrowser for help on using the repository browser.