source: trunk/INTERP2/newinterp.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: 6.8 KB
Line 
1;------------------------------------------------------------
2;------------------------------------------------------------
3;------------------------------------------------------------
4;
5; NAME: newinterp.pro
6;
7; PURPOSE: Uses weights to compute interpolated field
8;
9; CATEGORY : Subroutine
10;
11; CALLING SEQUENCE : newinterp
12;
13; INPUTS :
14;         ORCA_GEO grid file + ORCA_GEO weight file + Input field
15;         on ORCA grid
16; OUTPUTS :
17;         Output field on ORCA_GEO grid
18;
19; COMMON BLOCKS:
20;         None
21;
22; SIDE EFFECTS:
23;
24; RESTRICTIONS:
25;
26; EXAMPLE:
27;
28; MODIFICATION HISTORY: 08/2002 Robinson Hordoir
29;                       03/2003 Robinson Hordoir & Christian Ethe :
30;                       Added duplication stripes + extrapolation on
31;                       land points
32;
33;------------------------------------------------------------
34;------------------------------------------------------------
35;------------------------------------------------------------
36
37
38
39pro newinterp
40
41
42@initorca_reg
43@naminterp2
44@init_path2
45
46
47; Get the output mask
48
49mask_reg=read_ncdf('tmask',FILENAME=outputgrid,/NOSTRUCT,/cont_nofill)
50
51
52; Get the Pointers and weights...
53
54weight=replicate(0.,jpi,jpj,4)
55pointer_i=replicate(999,jpi,jpj,4)
56pointer_j=replicate(999,jpi,jpj,4)
57
58weight(*,*,0)=read_ncdf('weight_1', FILENAME=weightfile,/NOSTRUCT,/cont_nofill)
59weight(*,*,1)=read_ncdf('weight_2', FILENAME=weightfile,/NOSTRUCT,/cont_nofill)
60weight(*,*,2)=read_ncdf('weight_3', FILENAME=weightfile,/NOSTRUCT,/cont_nofill)
61weight(*,*,3)=read_ncdf('weight_4', FILENAME=weightfile,/NOSTRUCT,/cont_nofill)
62
63
64pointer_i(*,*,0)=read_ncdf('pointer_1_i', FILENAME=weightfile,/NOSTRUCT,/cont_nofill)
65pointer_i(*,*,1)=read_ncdf('pointer_2_i', FILENAME=weightfile,/NOSTRUCT,/cont_nofill)
66pointer_i(*,*,2)=read_ncdf('pointer_3_i', FILENAME=weightfile,/NOSTRUCT,/cont_nofill)
67pointer_i(*,*,3)=read_ncdf('pointer_4_i', FILENAME=weightfile,/NOSTRUCT,/cont_nofill)
68
69
70pointer_j(*,*,0)=read_ncdf('pointer_1_j', FILENAME=weightfile,/NOSTRUCT,/cont_nofill)
71pointer_j(*,*,1)=read_ncdf('pointer_2_j', FILENAME=weightfile,/NOSTRUCT,/cont_nofill)
72pointer_j(*,*,2)=read_ncdf('pointer_3_j', FILENAME=weightfile,/NOSTRUCT,/cont_nofill)
73pointer_j(*,*,3)=read_ncdf('pointer_4_j', FILENAME=weightfile,/NOSTRUCT,/cont_nofill)
74
75; Get the field to be interpolated on reg geogrid
76
77input_field=ncdf_lec(datafile,var=name_input_field)
78
79; ici mettre : continent=where(input_field GT 1.E13)
80;              input_field(continent)=0.
81
82input_mask =ncdf_lec(inputmesh,var='tmask')
83
84zmask=input_mask
85
86
87IF keyword_set(do_rempl) then begin
88
89; Extrapolate the input field on masked values
90for jt=0,ntimesteps-1 do begin
91   for jk=0, nlevel-1 do begin
92
93    not_extrap=input_field(*,*,jk,jt)
94
95    if (max(input_mask(*,*,jk)) GT 0 ) then input_field(*,*,jk,jt)=remplit(not_extrap,NITE=niterations,BASIQUE=input_mask(*,*,jk))
96
97   endfor
98endfor
99endif
100
101
102; Get the output field
103
104output_field=replicate(0.,jpi,jpj,nlevel,ntimesteps)
105
106
107for ji=0,jpi-1 do begin
108
109   for jj=0,jpj-1 do begin
110   
111   
112   if mask_reg(ji,jj,0) NE 0 then begin
113
114   output_field(ji,jj,*,*)=                  $
115                  weight(ji,jj,0)*input_field(pointer_i(ji,jj,0)+1,pointer_j(ji,jj,0),*,*) $
116                + weight(ji,jj,1)*input_field(pointer_i(ji,jj,1)+1,pointer_j(ji,jj,1),*,*) $
117                + weight(ji,jj,2)*input_field(pointer_i(ji,jj,2)+1,pointer_j(ji,jj,2),*,*) $
118                + weight(ji,jj,3)*input_field(pointer_i(ji,jj,3)+1,pointer_j(ji,jj,3),*,*)
119   
120   endif
121
122   endfor
123
124endfor
125
126; Initialisation of array for the northern boundary condition
127jplus = 1
128zdata_north = replicate(0.,jpi,jpj+jplus,nlevel,ntimesteps)
129zdata_north(*,0:jpj-1,*,*) = output_field
130zdata_inv = replicate(0.,jpi,jpj+jplus,nlevel,ntimesteps)
131; Data of the input grid is duplicated for the northern boundary condition
132
133zdata_inv(0:jpi/2 -1 ,*,*,*) = zdata_north(jpi/2:jpi-1,*,*,*)
134zdata_inv(jpi/2:jpi-1 ,*,*,*) = zdata_north(0:jpi/2-1,*,*,*)
135zdata_inv = reverse(zdata_inv,2)
136zdata_north(*,jpj:jpj+jplus-1,*,*) = zdata_inv(*,jpj:jpj+jplus-1,*,*)
137
138; Data of the input grid is duplicated for the West/east boundary condition
139final_output=replicate(0.,jpiglo,jpjglo,nlevel,ntimesteps)
140final_output(1:jpiglo-2,0:jpjglo-1,*,*) = zdata_north
141
142final_output(0,*,*,*) = final_output(jpiglo-2,*,*,*)
143final_output(jpiglo-1,*,*,*) = final_output(1,*,*,*)
144
145
146IF keyword_set(nlower_limit) THEN BEGIN
147    tab_min = WHERE( final_output LT min_value )
148    IF tab_min(0) NE -1 THEN final_output(tab_min) = min_value
149ENDIF
150
151IF keyword_set(nupper_limit) THEN BEGIN
152    tab_max = WHERE( final_output GT max_value )
153    IF tab_max(0) NE -1 THEN final_output(tab_max) = max_value
154ENDIF
155
156final_time=ncdf_lec(datafile,var=data_time_variable)
157
158nav_lon=ncdf_lec(outputgrid,var='nav_lon')
159nav_lat=ncdf_lec(outputgrid,var='nav_lat')
160nav_lev=ncdf_lec(outputgrid,var='nav_lev')
161
162vargrid = 'T'
163
164
165; Name
166   idout = NCDF_CREATE(outputfile,/clobber)
167   print, 'Creation du fichier Netcdf'
168   NCDF_CONTROL, idout, /nofill
169; Dimension
170
171   xidout = NCDF_DIMDEF(idout, 'x', jpiglo)
172   yidout = NCDF_DIMDEF(idout, 'y', jpjglo)
173   didout = NCDF_DIMDEF(idout, 'z', nlevel)
174   x_a    = NCDF_DIMDEF(idout, 'x_a',1)
175   y_a    = NCDF_DIMDEF(idout, 'y_a',1)
176   z_a    = NCDF_DIMDEF(idout, 'z_a',1)
177   tidout = NCDF_DIMDEF(idout, 'time_counter', /unlimited)
178
179; Attributes
180     
181   id0  = NCDF_VARDEF(idout, 'nav_lon'     , [xidout, yidout                ],/FLOAT )
182   id1  = NCDF_VARDEF(idout, 'nav_lat'     , [xidout, yidout                ],/FLOAT )
183   id2  = NCDF_VARDEF(idout, 'nav_lev'     , [                didout        ],/FLOAT )
184   id3  = NCDF_VARDEF(idout, 'time_counter'        , [                        tidout],/FLOAT )
185   id4  = NCDF_VARDEF(idout, 'time_steps'  , [                        tidout],/LONG  )
186   id5  = NCDF_VARDEF(idout,name_input_field,[xidout, yidout,didout,  tidout],/DOUBLE)
187
188
189; Variable 0
190   NCDF_ATTPUT, idout, id0, 'units', 'degrees_east'
191   NCDF_ATTPUT, idout, id0, 'long_name', 'Longitude'
192   NCDF_ATTPUT, idout, id0, 'nav_model', 'Default grid'
193; Variable 1
194   NCDF_ATTPUT, idout, id1, 'units', 'degrees_north'
195   NCDF_ATTPUT, idout, id1, 'long_name', 'Latitude'
196   NCDF_ATTPUT, idout, id1, 'nav_model', 'Default grid'
197; Variable 2
198   NCDF_ATTPUT, idout, id2, 'units','meters'
199   NCDF_ATTPUT, idout, id2, 'long_name','Depth'
200   NCDF_ATTPUT, idout, id2, 'nav_model','Default grid'
201; Variable 3
202   NCDF_ATTPUT, idout, id3, 'units', 'seconds since 0001-01-15 12:00:00 '
203   NCDF_ATTPUT, idout, id3, 'calendar','noleap'
204   NCDF_ATTPUT, idout, id3, 'title', 'Time'
205   NCDF_ATTPUT, idout, id3, 'long_name', 'Time axis'
206   NCDF_ATTPUT, idout, id3, 'time_origin','0000-DEC-15 00:00:00'
207
208
209 
210   NCDF_CONTROL, idout, /ENDEF
211
212; Writting
213
214   NCDF_VARPUT, idout, id0 , nav_lon
215   NCDF_VARPUT, idout, id1 , nav_lat
216   NCDF_VARPUT, idout, id2 , nav_lev(0:nlevel-1)
217   NCDF_VARPUT, idout, id3 , final_time
218   NCDF_VARPUT, idout, id4 , ntimesteps
219   NCDF_VARPUT, idout, id5 , final_output
220
221
222   NCDF_CLOSE, idout
223
224 
225   
226
227end
228
229
230
231
232
233
234
235
Note: See TracBrowser for help on using the repository browser.