source: trunk/extrap.pro

Last change on this file was 48, checked in by pinsard, 10 years ago

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