source: trunk/extrap.pro @ 2

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

initial import from /usr/work/fvi/OPA/geomag/
File size: 5.0 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;+
5; NAME: extrap.pro (based on remplit.pro)
6;
7; PURPOSE: Extrapolates ocean data values on land points
8;
9; CATEGORY: Subroutine
10;
11; CALLING SEQUENCE: extrap, zdata, zmask, it, val
12;
13; INPUTS:
14;          zdata : field to extrapolate
15;          zmask : field's mask
16;          it    : iteration
17;          val   : field value on masked points 
18;
19; KEYWORD PARAMETERS: None
20;
21; OUTPUTS:
22;          zdata : extrapolated field
23;          zmask : mask after extrapolation
24;
25; COMMON BLOCKS:
26;       common_interp.pro
27;
28; SIDE EFFECTS:
29;
30; RESTRICTIONS:
31;
32; EXAMPLE:
33;
34; MODIFICATION HISTORY: 19/11/99 Arnaud Jouzeau
35;                       25/02/00 Sebastien Masson (remplit.pro)
36;-
37;------------------------------------------------------------
38;------------------------------------------------------------
39;------------------------------------------------------------
40PRO extrap, z, mmmask, it, val
41@common_interp
42;
43; 1. Looking for land points
44; ==========================
45;
46; IF keyword_set(key_performance) THEN temp = systime(1)
47; znbvoisins=shift(zmask,1,0)+shift(zmask,-1,0)+shift(zmask,0,1)+shift(zmask,0,-1) $
48;           +shift(zmask,1,1)+shift(zmask,-1,1)+shift(zmask,1,-1)+shift(zmask,-1,-1)
49; IF keyword_set(key_performance) THEN print, systime(1)-temp
50;
51; 2. Extrapolating when needed
52; ============================
53;
54; IF keyword_set(key_performance) THEN temp = systime(1)
55; zdfm = zdata*zmask
56; zremplit =shift(zdfm,1,0)+shift(zdfm,-1,0)+shift(zdfm,0,1)+shift(zdfm, 0,-1) $
57;          +shift(zdfm,1,1)+shift(zdfm,-1,1)+shift(zdfm,1,-1)+shift(zdfm, -1,-1)
58; IF keyword_set(key_performance) THEN print, systime(1)-temp
59;
60;  ... coastal points
61;
62; IF keyword_set(key_performance) THEN temp = systime(1)
63; ocean = where(zmask NE 0.)
64; IF keyword_set(key_performance) THEN print, systime(1)-temp
65; IF keyword_set(key_performance) THEN temp = systime(1)
66; ind = where(znbvoisins NE 0. AND znbvoisins NE 8.)
67; IF keyword_set(key_performance) THEN print, systime(1)-temp
68; IF keyword_set(key_performance) THEN temp = systime(1)
69; znbvoisins(ind) = 1.d/temporary(znbvoisins(ind))
70; IF keyword_set(key_performance) THEN print, systime(1)-temp
71; IF keyword_set(key_performance) THEN temp = systime(1)
72; zremplit(ind) = temporary(zremplit(ind))*znbvoisins(ind)
73; IF keyword_set(key_performance) THEN print, systime(1)-temp
74; IF keyword_set(key_performance) THEN temp = systime(1)
75; zmask(ind) = 1.d
76; IF keyword_set(key_performance) THEN print, systime(1)-temp
77;
78;  ... inland points
79;
80; IF keyword_set(key_performance) THEN temp = systime(1)
81; ocean2 = where(zmask NE 0.)
82; IF keyword_set(key_performance) THEN print, systime(1)-temp
83; IF keyword_set(key_performance) THEN temp = systime(1)
84; IF n_elements(ocean2) NE jpiatm*(jpjatm+4) THEN BEGIN
85; land = where(zmask EQ 0.)
86; zremplit(land) = val
87; ENDIF
88; IF keyword_set(key_performance) THEN print, systime(1)-temp
89;
90;  ... filling with ocean values
91;
92; IF keyword_set(key_performance) THEN temp = systime(1)
93; zremplit(ocean)=zdata(ocean)
94; zdata = zremplit
95; IF keyword_set(key_performance) THEN print, systime(1)-temp
96; IF keyword_set(key_performance) THEN print, 'fin iteration', it
97;
98; on trouve les points cote
99   tempdeux = systime(1)        ; pour key_performance =2
100; les points du bord du cadre ne doivent pas etre dans la cote
101   mmmask[0, *] = 1
102   mmmask[jpiatm+1, *] = 1
103   mmmask[*, 0] = 1
104   mmmask[*, jpjatm+1] = 1
105; liste des points terre restant
106   terre = where(mmmask EQ 0)
107   if terre[0] EQ -1 then GOTO, fini
108; les points du bord du cadre doivent maintenan etre dans la terre
109   mmmask[0, *] = 0
110   mmmask[jpiatm+1, *] = 0
111   mmmask[*, 0] = 0
112   mmmask[*, jpjatm+1] = 0
113; liste des voisins mer
114   voisins = shift(mmmask,-1,0)+shift(mmmask,1,0)+shift(mmmask,0,-1)+shift(mmmask,0,1) $
115    +1./sqrt(2)*(shift(mmmask,-1,-1)+shift(mmmask,1,1)+shift(mmmask,1,-1)+shift(mmmask,-1,1))
116; liste des points cote
117   cote = where(voisins[terre] GT 0)
118   cote = terre[cote]
119;
120   IF testvar(var = key_performance) EQ 2 THEN $
121    print, 'temps remplit: trouver la cote ', systime(1)-tempdeux
122;---------------------------------------------------------------
123; remplissage des points cote
124;---------------------------------------------------------------
125   tempdeux = systime(1)        ; pour key_performance =2
126; on masque z
127   z = z*mmmask
128;
129   zcote = z[1+cote]+z[-1+cote]+z[(jpiatm+2)+cote]+z[-(jpiatm+2)+cote] $
130    +1./sqrt(2)*(z[(jpiatm+2)+1+cote]+z[(jpiatm+2)-1+cote]+z[-(jpiatm+2)+1+cote]+z[-(jpiatm+2)-1+cote])
131   poids = voisins[cote]
132     
133   z[cote] = zcote/poids
134;---------------------------------------------------------------
135; IV) on reduit le masque
136;---------------------------------------------------------------
137   mmmask[cote] = 1
138   z = z*mmmask + val*(1-mmmask)
139;   print, max(z)
140;
141   IF testvar(var = key_performance) EQ 2 THEN $
142    print, 'temps remplit: une iteration ', systime(1)-tempdeux
143fini:
144
145return
146END
Note: See TracBrowser for help on using the repository browser.