source: trunk/SRC/ToBeReviewed/CALCULS/remplit.pro @ 285

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

corrections of some headers and parameters and keywords case. change of pro2href to replace proidl
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 8.2 KB
Line 
1;+
2; @file_comments
3;
4; @categories
5;
6; @param ZINPUT
7;
8; @keyword NAN
9;
10; @keyword NITER
11;
12; @keyword BASIQUE
13;
14; @keyword MASK
15;
16; @keyword FILLXDIR
17;
18; @keyword FILLYDIR
19;
20; @keyword FILLVAL
21;
22; @keyword _EXTRA
23; Used to pass keywords
24;
25; @returns
26;
27; @uses
28;
29; @restrictions
30;
31; @examples
32;
33; @history
34;
35; @version
36; $Id$
37;
38;;
39;; Extrapole zinout[jpi,jpj] sur les continents en utilisant les 4
40;; plus proches voisins masques oceaniquement et construit un nouveau masque
41;; contenant l'ancien masque oceanique PLUS les points extrapoles.
42;; Reitere le processus niter fois.
43;; C'est pas clair, essayez !
44;;
45;;
46;
47;    /Nan: to fill the point which have the value
48;    !values.f_nan. Without this keyword, these point are not filling
49;    and stays at !values.f_nan.
50;
51;
52; @todo seb
53;
54;-
55;
56FUNCTION remplit, zinput, NAN = nan, NITER = niter, BASIQUE = basique, MASK = mask, FILLXDIR = fillxdir, FILLYDIR = fillydir, FILLVAL = fillval, _EXTRA = ex
57;
58  compile_opt idl2, strictarrsubs
59;
60@common
61  tempsun = systime(1)          ; pour key_performance
62; les points non remplis sont masques a valmask
63  IF n_elements(niter) EQ 0 THEN niter = 1
64  IF niter EQ 0 THEN return, zinput
65  z = zinput
66  if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c'
67  oldkey_gridtype = key_gridtype
68; keep only the first letter of the grid type
69  key_gridtype = strmid(key_gridtype, 0, 1)
70  if keyword_set(basique) then begin
71    key_gridtype = 'c'
72    nx = (size(zinput))[1]
73    ny = (size(zinput))[2]
74    if NOT keyword_set(mask) then mmmask = basique ELSE mmmask = mask
75  ENDIF ELSE grille, mmmask, glam, gphi, gdep, nx, ny, nz, _extra = ex
76  if keyword_set(mask) then mmmask = mask
77;---------------------------------------------------------------
78  if (size(mmmask))[0] EQ 3 THEN mmmask = mmmask[*, *, 0]
79;
80  if n_elements(mmmask) EQ 1 then mmmask = replicate(1b, nx, ny)
81  if keyword_set(nan) then begin
82    nanpoint = where(finite(z) EQ 0)
83    if nanpoint[0] NE -1 then begin
84      mmmask[nanpoint] = 0b
85      z[nanpoint] = 0
86    endif
87  ENDIF
88  mmmask = byte(mmmask)
89;---------------------------------------------------------------
90; on ajoute un cadre de zero a z, mask, e1, e2
91; comme ca apres on peut faire des shifts ds tous les sens sans se
92; soucier des bords du domaine!
93;---------------------------------------------------------------
94  tempdeux = systime(1)         ; pour key_performance =2
95  nx2 = nx+2
96  case key_gridtype of
97    'c':BEGIN
98      ztmp = bytarr(nx+2, ny+2)
99      ztmp[1:nx, 1:ny] = mmmask
100      mmmask = temporary(ztmp)
101      ztmp = fltarr(nx+2, ny+2)
102      ztmp[1:nx, 1:ny] = z
103      if keyword_set(key_periodic) AND nx EQ jpi then begin
104        ztmp[0, 1:ny] = z[jpi-1, *]
105        ztmp[nx+1, 1:ny] = z[0, *]
106      endif
107      z = temporary(ztmp)
108    END
109    'e':BEGIN
110      ztmp = bytarr(nx+2, ny+4)
111      ztmp[1:nx, 2:ny+1] = mmmask
112      mmmask = temporary(ztmp)
113      ztmp = fltarr(nx+2, ny+4)
114      ztmp[1:nx, 2:ny+1] = z
115      if keyword_set(key_periodic) AND nx EQ jpi then begin
116        ztmp[0, 2:ny+1] = z[jpi-1, *]
117        ztmp[nx+1, 2:ny+1] = z[0, *]
118      endif
119      z = temporary(ztmp)
120    END
121  endcase
122  IF testvar(var = key_performance) EQ 2 THEN $
123    print, 'temps remplit: on ajoute un cadre de zero ', systime(1)-tempdeux
124;---------------------------------------------------------------
125;---------------------------------------------------------------
126; iteration
127;---------------------------------------------------------------
128;---------------------------------------------------------------
129  FOR n = 1, niter DO BEGIN
130; on trouve les points coast
131    tempdeux = systime(1)       ; pour key_performance =2
132; les points du bord du cadre ne doivent pas etre selectionnes comme
133; la coast
134    case key_gridtype of
135      'c':BEGIN
136        mmmask[0, *] = 1b
137        mmmask[nx+1, *] = 1b
138        mmmask[*, 0] = 1b
139        mmmask[*, ny+1] = 1b
140      END
141      'e':BEGIN
142        mmmask[0, *] = 1b
143        mmmask[nx+1, *] = 1b
144        mmmask[*, 0:1] = 1b
145        mmmask[*, ny+2:ny+3] = 1b
146      END
147    endcase
148; liste des points terre restant
149    IF keyword_set(fillxdir) THEN BEGIN
150; we stop if all the lines, that contains data, have been filled
151      test = total(mmmask[1:nx, *], 1)
152      IF total((test EQ 0)+(test EQ nx)) EQ ny+2 THEN GOTO, fini
153    ENDIF
154    IF keyword_set(fillydir) THEN BEGIN
155; we stop if all the columns, that contains data, have been filled
156      test = total(mmmask[*, 1:ny], 2)
157      IF total((test EQ 0)+(test EQ ny)) EQ nx+2 THEN GOTO, fini
158    ENDIF
159    land = where(mmmask EQ 0)
160    if land[0] EQ -1 then GOTO, fini
161; les points du bord du cadre doivent maintenant etre dans la terre
162    case key_gridtype of
163      'c':BEGIN
164        mmmask[0, *] = 0b
165        mmmask[nx+1, *] = 0b
166        mmmask[*, 0] = 0b
167        mmmask[*, ny+1] = 0b
168      END
169      'e':BEGIN
170        mmmask[0, *] = 0b
171        mmmask[nx+1, *] = 0b
172        mmmask[*, 0:1] = 0b
173        mmmask[*, ny+2:ny+3] = 0b
174      END
175    endcase
176    if keyword_set(key_periodic) AND nx EQ jpi then begin
177      mmmask[0, *] = mmmask[nx, *]
178      mmmask[nx+1, *] = mmmask[1, *]
179    endif
180; liste des voisins mer
181    case key_gridtype of
182      'c':BEGIN
183        CASE 1 OF
184          keyword_set(fillxdir):weight = mmmask[1+land]+mmmask[-1+land]
185          keyword_set(fillydir):weight = mmmask[nx2+land]+mmmask[-nx2+land]
186          ELSE:weight = mmmask[1+land]+mmmask[-1+land]+mmmask[nx2+land]+mmmask[-nx2+land] $
187            +1./sqrt(2)*(mmmask[nx2+1+land]+mmmask[nx2-1+land] $
188                         +mmmask[-nx2+1+land]+mmmask[-nx2-1+land])
189        ENDCASE
190      END
191      'e':BEGIN
192        shifted = glam[0, 0] LT glam[0, 1]
193        oddeven = (land/nx2+1-shifted) MOD 2
194        weight = mmmask[1+land]+mmmask[-1+land] $
195          +mmmask[2*nx2+land]+mmmask[-2*nx2+land] $
196          +sqrt(2)*(mmmask[-nx2+oddeven+land]+mmmask[-nx2-1+oddeven+land] $
197                    +mmmask[nx2+oddeven+land]+mmmask[nx2-1+oddeven+land])
198      END
199    endcase
200
201    ok = where(weight GT 0)
202    weight = weight[ok]
203    coast = land[temporary(ok)]
204;
205    IF testvar(var = key_performance) EQ 2 THEN $
206      print, 'temps remplit: trouver la coast ', systime(1)-tempdeux
207;---------------------------------------------------------------
208; remplissage des points coast
209;---------------------------------------------------------------
210    tempdeux = systime(1)       ; pour key_performance =2
211; on masque z
212    z = temporary(z)*mmmask
213;
214    case key_gridtype of
215      'c':BEGIN
216        CASE 1 OF
217          keyword_set(fillxdir):zcoast = z[1+coast]+z[-1+coast]
218          keyword_set(fillydir):zcoast = z[nx2+coast]+z[-nx2+coast]
219          ELSE:zcoast = z[1+coast]+z[-1+coast]+z[nx2+coast]+z[-nx2+coast] $
220            +1./sqrt(2)*(z[nx2+1+coast]+z[nx2-1+coast] $
221                         +z[-nx2+1+coast]+z[-nx2-1+coast])
222        ENDCASE
223      END
224      'e':BEGIN
225        oddeven = (coast/nx2+1-shifted) MOD 2
226        zcoast = z[1+coast]+z[-1+coast]+z[2*nx2+coast]+z[-2*nx2+coast] $
227          +sqrt(2)*(z[-nx2+oddeven+coast]+z[-nx2-1+oddeven+coast] $
228                    +z[nx2+oddeven+coast]+z[nx2-1+oddeven+coast])
229      END
230    endcase
231;
232    z[coast] =  temporary(zcoast)/ temporary(weight)
233; we update the boundary conditions of z
234    if keyword_set(key_periodic) AND nx EQ jpi then begin
235      z[0, *] = z[nx, *]
236      z[nx+1, *] = z[1, *]
237    endif
238;---------------------------------------------------------------
239; IV) on reduit le masque
240;---------------------------------------------------------------
241    mmmask[ temporary(coast)] = 1
242;
243    IF testvar(var = key_performance) EQ 2 THEN $
244      print, 'temps remplit: une iteration ', systime(1)-tempdeux
245  ENDFOR
246fini:
247;
248; on masque les valeurs sur les lands restantes
249;
250  IF n_elements(valmask) EQ 0 then valmask = 1e20
251  IF n_elements(fillval) EQ 0 THEN fillval = valmask
252  z = temporary(z)*mmmask + fillval*(1b-mmmask)
253;---------------------------------------------------------------
254; on redecoupe le tableau pour retirer le cadre!
255;---------------------------------------------------------------
256  case key_gridtype of
257    'c':BEGIN
258      z = z[1:nx, 1:ny]
259    END
260    'e':BEGIN
261      z = z[1:nx, 2:ny+1]
262    END
263  endcase
264;
265  key_gridtype = oldkey_gridtype
266;---------------------------------------------------------------
267  if keyword_set(key_performance) THEN print, 'temps remplit', systime(1)-tempsun
268  return, z
269END
270
Note: See TracBrowser for help on using the repository browser.