1 | ! |
---|
2 | ! Aggregation routines. These routines allow to interpolate from the finer grid on which the |
---|
3 | ! surface parameter is available to the coarser one of the model. |
---|
4 | ! |
---|
5 | ! The routines work for the fine data on a regular lat/lon grid. This grid can come in as either |
---|
6 | ! a rank2 array or a vector. Two procedure exist which require slightly different input fields. |
---|
7 | ! |
---|
8 | ! |
---|
9 | !< $HeadURL$ |
---|
10 | !< $Date$ |
---|
11 | !< $Author$ |
---|
12 | !< $Revision$ |
---|
13 | ! |
---|
14 | ! |
---|
15 | MODULE interpol_help |
---|
16 | |
---|
17 | ! Modules used : |
---|
18 | |
---|
19 | USE constantes |
---|
20 | USE mod_orchidee_para |
---|
21 | |
---|
22 | IMPLICIT NONE |
---|
23 | |
---|
24 | PRIVATE |
---|
25 | PUBLIC aggregate, aggregate_p |
---|
26 | ! |
---|
27 | INTERFACE aggregate |
---|
28 | MODULE PROCEDURE aggregate_2d, aggregate_vec |
---|
29 | END INTERFACE |
---|
30 | ! |
---|
31 | INTERFACE aggregate_p |
---|
32 | MODULE PROCEDURE aggregate_2d_p, aggregate_vec_p |
---|
33 | END INTERFACE |
---|
34 | ! |
---|
35 | LOGICAL, PARAMETER :: check_grid=.FALSE. |
---|
36 | ! |
---|
37 | CONTAINS |
---|
38 | ! |
---|
39 | ! This routing will get for each point of the coarse grid the |
---|
40 | ! indexes of the finer grid and the area of overlap. |
---|
41 | ! This routine is designed for a fine grid which is regular in lat/lon. |
---|
42 | ! |
---|
43 | SUBROUTINE aggregate_2d (nbpt, lalo, neighbours, resolution, contfrac, & |
---|
44 | & iml, jml, lon_rel, lat_rel, mask, callsign, & |
---|
45 | & incmax, indinc, areaoverlap, ok) |
---|
46 | |
---|
47 | USE grid, ONLY : global |
---|
48 | |
---|
49 | ! |
---|
50 | ! INPUT |
---|
51 | ! |
---|
52 | INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated |
---|
53 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) |
---|
54 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W) |
---|
55 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y |
---|
56 | REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box. |
---|
57 | INTEGER(i_std), INTENT(in) :: iml, jml ! Size of the finer grid |
---|
58 | REAL(r_std), INTENT(in) :: lon_rel(iml, jml) ! Longitudes for the finer grid |
---|
59 | REAL(r_std), INTENT(in) :: lat_rel(iml, jml) ! Latitudes for the finer grid |
---|
60 | INTEGER(i_std), INTENT(in) :: mask(iml, jml) ! Mask which retains only the significative points |
---|
61 | ! of the fine grid. |
---|
62 | CHARACTER(LEN=*), INTENT(in) :: callsign ! Allows to specify which variable is beeing treated |
---|
63 | INTEGER(i_std), INTENT(in) :: incmax ! Maximum point of the fine grid we can store. |
---|
64 | ! |
---|
65 | ! Output |
---|
66 | ! |
---|
67 | INTEGER(i_std), INTENT(out) :: indinc(nbpt,incmax,2) |
---|
68 | REAL(r_std), INTENT(out) :: areaoverlap(nbpt,incmax) |
---|
69 | LOGICAL, OPTIONAL, INTENT(out) :: ok ! return code |
---|
70 | ! |
---|
71 | ! Local Variables |
---|
72 | ! |
---|
73 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_ful, lon_ful |
---|
74 | REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: loup_rel, lolow_rel, laup_rel, lalow_rel |
---|
75 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: searchind |
---|
76 | REAL(r_std) :: lon_up, lon_low, lat_up, lat_low |
---|
77 | REAL(r_std) :: coslat, ax, ay, sgn, lonrel, lolowrel, louprel |
---|
78 | INTEGER(i_std) :: fopt, fopt_max, ip, jp, ib, i, itmp, iprog, nbind |
---|
79 | REAL(r_std) :: domain_minlon,domain_maxlon,domain_minlat,domain_maxlat |
---|
80 | INTEGER(i_std) :: minLon(1), maxLon(1) |
---|
81 | |
---|
82 | INTEGER :: ALLOC_ERR |
---|
83 | LOGICAL :: err_fopt |
---|
84 | err_fopt = .FALSE. |
---|
85 | ! |
---|
86 | ! Some inital assignmens |
---|
87 | ! |
---|
88 | areaoverlap(:,:) = moins_un |
---|
89 | indinc(:,:,:) = zero |
---|
90 | |
---|
91 | ALLOCATE (laup_rel(iml,jml), STAT=ALLOC_ERR) |
---|
92 | IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of laup_rel','','') |
---|
93 | |
---|
94 | ALLOCATE (loup_rel(iml,jml), STAT=ALLOC_ERR) |
---|
95 | IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of loup_rel','','') |
---|
96 | |
---|
97 | ALLOCATE (lalow_rel(iml,jml), STAT=ALLOC_ERR) |
---|
98 | IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of lalow_rel','','') |
---|
99 | |
---|
100 | ALLOCATE (lolow_rel(iml,jml), STAT=ALLOC_ERR) |
---|
101 | IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of lolow_rel','','') |
---|
102 | |
---|
103 | ALLOCATE (lat_ful(iml+2,jml+2), STAT=ALLOC_ERR) |
---|
104 | IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of lat_ful','','') |
---|
105 | |
---|
106 | ALLOCATE (lon_ful(iml+2,jml+2), STAT=ALLOC_ERR) |
---|
107 | IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of lon_ful','','') |
---|
108 | |
---|
109 | ALLOCATE (searchind(iml*jml,2), STAT=ALLOC_ERR) |
---|
110 | IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of searchind','','') |
---|
111 | |
---|
112 | IF (PRESENT(ok)) ok = .TRUE. |
---|
113 | ! |
---|
114 | ! Duplicate the border assuming we have a global grid going from west to east |
---|
115 | ! |
---|
116 | lon_ful(2:iml+1,2:jml+1) = lon_rel(1:iml,1:jml) |
---|
117 | lat_ful(2:iml+1,2:jml+1) = lat_rel(1:iml,1:jml) |
---|
118 | ! |
---|
119 | IF ( lon_rel(iml,1) .LT. lon_ful(2,2)) THEN |
---|
120 | lon_ful(1,2:jml+1) = lon_rel(iml,1:jml) |
---|
121 | lat_ful(1,2:jml+1) = lat_rel(iml,1:jml) |
---|
122 | ELSE |
---|
123 | lon_ful(1,2:jml+1) = lon_rel(iml,1:jml)-360 |
---|
124 | lat_ful(1,2:jml+1) = lat_rel(iml,1:jml) |
---|
125 | ENDIF |
---|
126 | |
---|
127 | IF ( lon_rel(1,1) .GT. lon_ful(iml+1,2)) THEN |
---|
128 | lon_ful(iml+2,2:jml+1) = lon_rel(1,1:jml) |
---|
129 | lat_ful(iml+2,2:jml+1) = lat_rel(1,1:jml) |
---|
130 | ELSE |
---|
131 | lon_ful(iml+2,2:jml+1) = lon_rel(1,1:jml)+360 |
---|
132 | lat_ful(iml+2,2:jml+1) = lat_rel(1,1:jml) |
---|
133 | ENDIF |
---|
134 | ! |
---|
135 | sgn = lat_rel(1,1)/ABS(lat_rel(1,1)) |
---|
136 | lat_ful(2:iml+1,1) = sgn*180 - lat_rel(1:iml,1) |
---|
137 | sgn = lat_rel(1,jml)/ABS(lat_rel(1,jml)) |
---|
138 | lat_ful(2:iml+1,jml+2) = sgn*180 - lat_rel(1:iml,jml) |
---|
139 | lat_ful(1,1) = lat_ful(iml+1,1) |
---|
140 | lat_ful(iml+2,1) = lat_ful(2,1) |
---|
141 | lat_ful(1,jml+2) = lat_ful(iml+1,jml+2) |
---|
142 | lat_ful(iml+2,jml+2) = lat_ful(2,jml+2) |
---|
143 | ! |
---|
144 | ! Add the longitude lines to the top and bottom |
---|
145 | ! |
---|
146 | lon_ful(:,1) = lon_ful(:,2) |
---|
147 | lon_ful(:,jml+2) = lon_ful(:,jml+1) |
---|
148 | ! |
---|
149 | ! Get the upper and lower limits of each grid box |
---|
150 | ! |
---|
151 | DO ip=1,iml |
---|
152 | DO jp=1,jml |
---|
153 | ! |
---|
154 | loup_rel(ip,jp) =MAX(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)),& |
---|
155 | & 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1))) |
---|
156 | lolow_rel(ip,jp) =MIN(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)),& |
---|
157 | & 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1))) |
---|
158 | laup_rel(ip,jp) =MAX(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)),& |
---|
159 | & 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2))) |
---|
160 | lalow_rel(ip,jp) =MIN(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)),& |
---|
161 | & 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2))) |
---|
162 | ! |
---|
163 | ENDDO |
---|
164 | ENDDO |
---|
165 | IF (check_grid) THEN |
---|
166 | WRITE(numout,*) "================================" |
---|
167 | WRITE(numout,*) "interpol_aggregate_2d : " |
---|
168 | WRITE(numout,*) "lalo(:,1) :",lalo(:,1) |
---|
169 | WRITE(numout,*) "lalo(:,2) :",lalo(:,2) |
---|
170 | WRITE(numout,*) "Map meshes : " |
---|
171 | WRITE(numout,*) "lat read(1,:) :",lat_rel(1,:) |
---|
172 | WRITE(numout,*) "lat_ful(1,:) :",lat_ful(1,:) |
---|
173 | WRITE(numout,*) "lat_ful(2,:) :",lat_ful(2,:) |
---|
174 | WRITE(numout,*) "lalow_rel(1,:) :",lalow_rel(1,:) |
---|
175 | WRITE(numout,*) "laup_rel(1,:) :",laup_rel(1,:) |
---|
176 | WRITE(numout,*) "================================" |
---|
177 | WRITE(numout,*) "lon read(:,1) :",lon_rel(:,1) |
---|
178 | WRITE(numout,*) "lon_ful(:,1) :",lon_ful(:,1) |
---|
179 | WRITE(numout,*) "lon_ful(:,2) :",lon_ful(:,2) |
---|
180 | WRITE(numout,*) "lolow_rel(:,1) :",lolow_rel(:,1) |
---|
181 | WRITE(numout,*) "loup_rel(:,1) :",loup_rel(:,1) |
---|
182 | WRITE(numout,*) "================================" |
---|
183 | ENDIF |
---|
184 | ! |
---|
185 | ! |
---|
186 | ! To speedup calculations we will get the limits of the domain of the |
---|
187 | ! coarse grid and select all the points of the fine grid which are potentialy |
---|
188 | ! in this domain. |
---|
189 | ! |
---|
190 | ! |
---|
191 | minLon = MINLOC(lalo(1:nbpt,2)) |
---|
192 | coslat = MAX(COS(lalo(minLon(1),1) * pi/180. ), mincos )*pi/180. * R_Earth |
---|
193 | domain_minlon = lalo(minLon(1),2) - resolution(minLon(1),1)/(2.0*coslat) |
---|
194 | maxLon = MAXLOC(lalo(1:nbpt,2)) |
---|
195 | coslat = MAX(COS(lalo(maxLon(1),1) * pi/180. ), mincos )*pi/180. * R_Earth |
---|
196 | domain_maxlon = lalo(maxLon(1),2) + resolution(maxLon(1),1)/(2.0*coslat) |
---|
197 | ! |
---|
198 | coslat = pi/180. * R_Earth |
---|
199 | domain_minlat = MINVAL(lalo(1:nbpt,1)) - resolution(maxLon(1),2)/(2.0*coslat) |
---|
200 | domain_maxlat = MAXVAL(lalo(1:nbpt,1)) + resolution(maxLon(1),2)/(2.0*coslat) |
---|
201 | ! |
---|
202 | IF (check_grid) THEN |
---|
203 | WRITE(numout,*) "indices min/max of longitude :",minLon,maxLon, & |
---|
204 | & "; longitude min/max : ",lalo(minLon,1),lalo(maxLon,1) |
---|
205 | WRITE(numout,*) "Domain for coarse grid :" |
---|
206 | WRITE(numout,*) '(',domain_minlat,',',domain_minlon,')',& |
---|
207 | & '(',domain_maxlat,',',domain_maxlon,')' |
---|
208 | WRITE(numout,*) "================================" |
---|
209 | ENDIF |
---|
210 | ! |
---|
211 | ! we list a first approximation of all point we will need to |
---|
212 | ! scan to fill our coarse grid. |
---|
213 | ! |
---|
214 | IF ( global ) THEN |
---|
215 | ! Here we do the entire globe |
---|
216 | WRITE(numout,*) 'In aggregate_p : do interpolation to global model domain' |
---|
217 | nbind=0 |
---|
218 | DO ip=1,iml |
---|
219 | DO jp=1,jml |
---|
220 | IF (mask(ip,jp) == 1 ) THEN |
---|
221 | nbind = nbind + 1 |
---|
222 | searchind(nbind,1) = ip |
---|
223 | searchind(nbind,2) = jp |
---|
224 | ENDIF |
---|
225 | ENDDO |
---|
226 | ENDDO |
---|
227 | ! |
---|
228 | ELSE |
---|
229 | ! Now we get a limited number of points |
---|
230 | WRITE(numout,*) 'In aggregate_p : do interpolation to regional model domain' |
---|
231 | nbind=0 |
---|
232 | DO ip=1,iml |
---|
233 | DO jp=1,jml |
---|
234 | IF ( loup_rel(ip,jp) >= domain_minlon .AND. lolow_rel(ip,jp) <= domain_maxlon .AND.& |
---|
235 | & laup_rel(ip,jp) >= domain_minlat .AND. lalow_rel(ip,jp) <= domain_maxlat ) THEN |
---|
236 | IF (mask(ip,jp) == 1 ) THEN |
---|
237 | nbind = nbind + 1 |
---|
238 | searchind(nbind,1) = ip |
---|
239 | searchind(nbind,2) = jp |
---|
240 | ENDIF |
---|
241 | ENDIF |
---|
242 | ENDDO |
---|
243 | ENDDO |
---|
244 | ENDIF |
---|
245 | ! |
---|
246 | WRITE(numout,*) 'We will work with ', nbind, ' points of the fine grid' |
---|
247 | ! |
---|
248 | WRITE(numout,*) 'Aggregate_2d : ', callsign |
---|
249 | #ifdef INTERPOL_ADVANCE |
---|
250 | WRITE(numout,'(2a40)')'0%--------------------------------------', & |
---|
251 | & '------------------------------------100%' |
---|
252 | #endif |
---|
253 | ! |
---|
254 | ! Now we take each grid point and find out which values from the forcing we need to average |
---|
255 | ! |
---|
256 | fopt_max = -1 |
---|
257 | DO ib =1, nbpt |
---|
258 | ! |
---|
259 | ! Give a progress meter |
---|
260 | ! |
---|
261 | #ifdef INTERPOL_ADVANCE |
---|
262 | iprog = NINT(REAL(ib,r_std)/REAL(nbpt,r_std)*79.) - NINT(REAL(ib-1,r_std)/REAL(nbpt,r_std)*79.) |
---|
263 | IF ( iprog .NE. 0 ) THEN |
---|
264 | WRITE(numout,'(a1,$)') 'x' |
---|
265 | ENDIF |
---|
266 | #endif |
---|
267 | ! |
---|
268 | ! We find the 4 limits of the grid-box. As we transform the resolution of the model |
---|
269 | ! into longitudes and latitudes we do not have the problem of periodicity. |
---|
270 | ! coslat is a help variable here ! |
---|
271 | ! |
---|
272 | coslat = MAX(COS(lalo(ib,1) * pi/180. ), mincos )*pi/180. * R_Earth |
---|
273 | ! |
---|
274 | lon_up = lalo(ib,2) + resolution(ib,1)/(2.0*coslat) |
---|
275 | lon_low =lalo(ib,2) - resolution(ib,1)/(2.0*coslat) |
---|
276 | ! |
---|
277 | coslat = pi/180. * R_Earth |
---|
278 | ! |
---|
279 | lat_up =lalo(ib,1) + resolution(ib,2)/(2.0*coslat) |
---|
280 | lat_low =lalo(ib,1) - resolution(ib,2)/(2.0*coslat) |
---|
281 | ! |
---|
282 | ! Find the grid boxes from the data that go into the model's boxes |
---|
283 | ! We still work as if we had a regular grid ! Well it needs to be localy regular so |
---|
284 | ! so that the longitude at the latitude of the last found point is close to the one |
---|
285 | ! of the next point. |
---|
286 | ! |
---|
287 | fopt = zero |
---|
288 | ! |
---|
289 | DO i=1,nbind |
---|
290 | ! |
---|
291 | ip = searchind(i,1) |
---|
292 | jp = searchind(i,2) |
---|
293 | ! |
---|
294 | ! Either the center of the data grid point is in the interval of the model grid or |
---|
295 | ! the East and West limits of the data grid point are on either sides of the border of |
---|
296 | ! the data grid. |
---|
297 | ! |
---|
298 | ! To do that correctly we have to check if the grid box sits on the date-line. |
---|
299 | ! |
---|
300 | IF ( lon_low < -180.0 ) THEN |
---|
301 | ! -179 -> -179 |
---|
302 | ! 179 -> -181 |
---|
303 | lonrel = MOD( lon_rel(ip,jp) - 360.0, 360.0) |
---|
304 | lolowrel = MOD( lolow_rel(ip,jp) - 360.0, 360.0) |
---|
305 | louprel = MOD( loup_rel(ip,jp) - 360.0, 360.0) |
---|
306 | ! |
---|
307 | ELSE IF ( lon_up > 180.0 ) THEN |
---|
308 | ! -179 -> 181 |
---|
309 | ! 179 -> 179 |
---|
310 | lonrel = MOD( lon_rel(ip,jp) + 360., 360.0) |
---|
311 | lolowrel = MOD( lolow_rel(ip,jp) + 360., 360.0) |
---|
312 | louprel = MOD( loup_rel(ip,jp) + 360., 360.0) |
---|
313 | ELSE |
---|
314 | lonrel = lon_rel(ip,jp) |
---|
315 | lolowrel = lolow_rel(ip,jp) |
---|
316 | louprel = loup_rel(ip,jp) |
---|
317 | ENDIF |
---|
318 | ! |
---|
319 | ! |
---|
320 | ! |
---|
321 | IF ( lonrel > lon_low .AND. lonrel < lon_up .OR. & |
---|
322 | & lolowrel < lon_low .AND. louprel > lon_low .OR. & |
---|
323 | & lolowrel < lon_up .AND. louprel > lon_up ) THEN |
---|
324 | ! |
---|
325 | ! Now that we have the longitude let us find the latitude |
---|
326 | ! |
---|
327 | IF ( lat_rel(ip,jp) > lat_low .AND. lat_rel(ip,jp) < lat_up .OR. & |
---|
328 | & lalow_rel(ip,jp) < lat_low .AND. laup_rel(ip,jp) > lat_low .OR.& |
---|
329 | & lalow_rel(ip,jp) < lat_up .AND. laup_rel(ip,jp) > lat_up) THEN |
---|
330 | ! |
---|
331 | fopt = fopt + 1 |
---|
332 | IF ( fopt > incmax) THEN |
---|
333 | err_fopt=.TRUE. |
---|
334 | EXIT |
---|
335 | ELSE |
---|
336 | ! |
---|
337 | ! If we sit on the date line we need to do the same transformations as above. |
---|
338 | ! |
---|
339 | IF ( lon_low < -180.0 ) THEN |
---|
340 | lolowrel = MOD( lolow_rel(ip,jp) - 360.0, 360.0) |
---|
341 | louprel = MOD( loup_rel(ip,jp) - 360.0, 360.0) |
---|
342 | ! |
---|
343 | ELSE IF ( lon_up > 180.0 ) THEN |
---|
344 | lolowrel = MOD( lolow_rel(ip,jp) + 360., 360.0) |
---|
345 | louprel = MOD( loup_rel(ip,jp) + 360., 360.0) |
---|
346 | ELSE |
---|
347 | lolowrel = lolow_rel(ip,jp) |
---|
348 | louprel = loup_rel(ip,jp) |
---|
349 | ENDIF |
---|
350 | ! |
---|
351 | ! Get the area of the fine grid in the model grid |
---|
352 | ! |
---|
353 | coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), mincos ) |
---|
354 | ax = (MIN(lon_up,louprel)-MAX(lon_low, lolowrel))*pi/180. * R_Earth * coslat |
---|
355 | ay = (MIN(lat_up, laup_rel(ip,jp))-MAX(lat_low,lalow_rel(ip,jp)))*pi/180. * R_Earth |
---|
356 | ! |
---|
357 | areaoverlap(ib, fopt) = ax*ay |
---|
358 | indinc(ib, fopt, 1) = ip |
---|
359 | indinc(ib, fopt, 2) = jp |
---|
360 | ! |
---|
361 | ! If this point was 100% within the grid then we can de-select it from our |
---|
362 | ! list as it can not be in another mesh of the coarse grid. |
---|
363 | ! |
---|
364 | IF ( louprel < lon_up .AND. lolowrel > lon_low .AND. & |
---|
365 | & laup_rel(ip,jp) < lat_up .AND. lalow_rel(ip,jp) > lat_low ) THEN |
---|
366 | searchind(i,1) = 0 |
---|
367 | searchind(i,2) = 0 |
---|
368 | ENDIF |
---|
369 | ! |
---|
370 | ENDIF |
---|
371 | ENDIF ! IF lat |
---|
372 | ENDIF ! IF lon |
---|
373 | ENDDO |
---|
374 | |
---|
375 | IF (err_fopt) THEN |
---|
376 | WRITE(numout,*) 'Working on variable :', callsign |
---|
377 | WRITE(numout,*) 'Reached value ', fopt,' for fopt on point', ib, lalo(ib,2), lalo(ib,1) |
---|
378 | CALL ipslerr_p(2,'aggregate_2d', & |
---|
379 | 'Working on variable :'//callsign, & |
---|
380 | 'Reached incmax value for fopt.',& |
---|
381 | 'Please increase incmax in subroutine calling aggregate') |
---|
382 | IF (PRESENT(ok)) THEN |
---|
383 | ok = .FALSE. |
---|
384 | RETURN |
---|
385 | ELSE |
---|
386 | CALL ipslerr_p(3,'aggregate_2d','Stop now','','') |
---|
387 | ENDIF |
---|
388 | ENDIF |
---|
389 | fopt_max = MAX ( fopt, fopt_max ) |
---|
390 | ! |
---|
391 | ! De-select the marked points |
---|
392 | ! |
---|
393 | itmp = nbind |
---|
394 | nbind = 0 |
---|
395 | DO i=1,itmp |
---|
396 | IF ( searchind(i,1) > 0 .AND. searchind(i,2) > 0 ) THEN |
---|
397 | nbind = nbind + 1 |
---|
398 | searchind(nbind,1) = searchind(i,1) |
---|
399 | searchind(nbind,2) = searchind(i,2) |
---|
400 | ENDIF |
---|
401 | ENDDO |
---|
402 | ! |
---|
403 | ENDDO |
---|
404 | ! |
---|
405 | DO ib=1,nbpt |
---|
406 | DO fopt=1,incmax |
---|
407 | IF (( indinc(ib,fopt,1) == 0 .AND. indinc(ib,fopt,2) > 0) .OR.& |
---|
408 | & ( indinc(ib,fopt,2) == 0 .AND. indinc(ib,fopt,1) > 0) ) THEN |
---|
409 | WRITE(*,*) "aggregate_2d PROBLEM : point =",ib, fopt," Indicies = ", & |
---|
410 | & indinc(ib,fopt,1), indinc(ib,fopt,2), areaoverlap(ib,fopt) |
---|
411 | ENDIF |
---|
412 | ENDDO |
---|
413 | ENDDO |
---|
414 | |
---|
415 | |
---|
416 | WRITE(numout,*) "" |
---|
417 | WRITE(numout,*) "aggregate_2D nbvmax = ",incmax, "max used = ",fopt_max |
---|
418 | ! |
---|
419 | ! Do some memory management. |
---|
420 | ! |
---|
421 | DEALLOCATE (laup_rel) |
---|
422 | DEALLOCATE (loup_rel) |
---|
423 | DEALLOCATE (lalow_rel) |
---|
424 | DEALLOCATE (lolow_rel) |
---|
425 | DEALLOCATE (lat_ful) |
---|
426 | DEALLOCATE (lon_ful) |
---|
427 | DEALLOCATE (searchind) |
---|
428 | ! |
---|
429 | ! Close the progress meter |
---|
430 | ! |
---|
431 | WRITE(numout,*) ' ' |
---|
432 | ! |
---|
433 | END SUBROUTINE aggregate_2d |
---|
434 | |
---|
435 | ! |
---|
436 | ! This routing will get for each point of the coarse grid the |
---|
437 | ! indexes of the finer grid and the area of overlap. |
---|
438 | ! This routine is designed for a fine grid which is regular in meters along lat lon axes. |
---|
439 | ! |
---|
440 | SUBROUTINE aggregate_vec (nbpt, lalo, neighbours, resolution, contfrac, & |
---|
441 | & iml, lon_rel, lat_rel, resol_lon, resol_lat, callsign, & |
---|
442 | & incmax, indinc, areaoverlap, ok) |
---|
443 | ! |
---|
444 | ! INPUT |
---|
445 | ! |
---|
446 | INTEGER(i_std), INTENT(in) :: nbpt ! Number of points for which the data needs to be interpolated |
---|
447 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) ! Vector of latitude and longitudes (beware of the order !) |
---|
448 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W) |
---|
449 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) ! The size in km of each grid-box in X and Y |
---|
450 | REAL(r_std), INTENT(in) :: contfrac(nbpt) ! Fraction of land in each grid box. |
---|
451 | INTEGER(i_std), INTENT(in) :: iml ! Size of the finer grid |
---|
452 | REAL(r_std), INTENT(in) :: lon_rel(iml) ! Longitudes for the finer grid |
---|
453 | REAL(r_std), INTENT(in) :: lat_rel(iml) ! Latitudes for the finer grid |
---|
454 | REAL(r_std), INTENT(in) :: resol_lon, resol_lat ! Resolution in meters of the fine grid |
---|
455 | CHARACTER(LEN=*), INTENT(in) :: callsign ! Allows to specify which variable is beeing treated |
---|
456 | INTEGER(i_std), INTENT(in) :: incmax ! Maximum point of the fine grid we can store. |
---|
457 | ! |
---|
458 | ! Output |
---|
459 | ! |
---|
460 | INTEGER(i_std), INTENT(out) :: indinc(nbpt,incmax) |
---|
461 | REAL(r_std), INTENT(out) :: areaoverlap(nbpt,incmax) |
---|
462 | LOGICAL, OPTIONAL, INTENT(out) :: ok ! return code |
---|
463 | ! |
---|
464 | ! Local Variables |
---|
465 | ! |
---|
466 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: searchind |
---|
467 | REAL(r_std) :: lon_up, lon_low, lat_up, lat_low |
---|
468 | REAL(r_std) :: coslat, ax, ay, lonrel, lolowrel, louprel |
---|
469 | REAL(r_std) :: latrel, lauprel, lalowrel |
---|
470 | INTEGER(i_std), DIMENSION(nbpt) :: fopt |
---|
471 | INTEGER(i_std) :: fopt_max, not_found_fopt |
---|
472 | INTEGER(i_std) :: ip, ib, i, j, itmp, iprog, nbind, pp, ipp |
---|
473 | REAL(r_std) :: domain_minlon,domain_maxlon,domain_minlat,domain_maxlat |
---|
474 | REAL(r_std) :: minlon, minlat, mini |
---|
475 | INTEGER(i_std) :: ff(1), incp |
---|
476 | INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: fine_ind |
---|
477 | INTEGER(i_std) :: pos_pnt(5) |
---|
478 | INTEGER :: ALLOC_ERR |
---|
479 | ! |
---|
480 | LOGICAL :: err_fopt |
---|
481 | err_fopt = .FALSE. |
---|
482 | ! |
---|
483 | ! Some inital assignmens |
---|
484 | ! |
---|
485 | areaoverlap(:,:) = moins_un |
---|
486 | indinc(:,:) = zero |
---|
487 | |
---|
488 | ALLOCATE (searchind(iml), STAT=ALLOC_ERR) |
---|
489 | IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_vec','ERROR IN ALLOCATION of searchbind','','') |
---|
490 | |
---|
491 | IF (PRESENT(ok)) ok = .TRUE. |
---|
492 | ! |
---|
493 | ! To speedup calculations we will get the limits of the domain of the |
---|
494 | ! coarse grid and select all the points of the fine grid which are potentialy |
---|
495 | ! in this domain. |
---|
496 | ! |
---|
497 | ! |
---|
498 | ff = MINLOC(lalo(:,2)) |
---|
499 | coslat = MAX(COS(lalo(ff(1),1) * pi/180. ), mincos )*pi/180. * R_Earth |
---|
500 | domain_minlon = lalo(ff(1),2) - resolution(ff(1),1)/(2.0*coslat) |
---|
501 | ff = MAXLOC(lalo(:,2)) |
---|
502 | coslat = MAX(COS(lalo(ff(1),1) * pi/180. ), mincos )*pi/180. * R_Earth |
---|
503 | domain_maxlon = lalo(ff(1),2) + resolution(ff(1),1)/(2.0*coslat) |
---|
504 | ! |
---|
505 | coslat = pi/180. * R_Earth |
---|
506 | domain_minlat = MINVAL(lalo(:,1)) - resolution(ff(1),2)/(2.0*coslat) |
---|
507 | domain_maxlat = MAXVAL(lalo(:,1)) + resolution(ff(1),2)/(2.0*coslat) |
---|
508 | ! |
---|
509 | ! Find appropriate resolution for index table |
---|
510 | ! |
---|
511 | ff=MINLOC(resolution(:,1)) |
---|
512 | coslat = MAX(COS(lalo(ff(1),1) * pi/180. ), mincos )*pi/180. * R_Earth |
---|
513 | minlon=resolution(ff(1),1)/(2.0*coslat) |
---|
514 | ff=MINLOC(resolution(:,2)) |
---|
515 | coslat = pi/180. * R_Earth |
---|
516 | minlat=resolution(ff(1),2)/(2.0*coslat) |
---|
517 | mini=MIN(minlon, minlat) |
---|
518 | ! |
---|
519 | ! This interpolation only works if the model grid is coarser than the data grid |
---|
520 | ! |
---|
521 | IF (MINVAL(resolution(:,1)) < resol_lon .OR. MINVAL(resolution(:,2)) < resol_lat) THEN |
---|
522 | WRITE(numout,*) " === WARNING == " |
---|
523 | WRITE(numout,*) "Resolution minima of the model (lon, lat) : ", & |
---|
524 | & MINVAL(resolution(:,1)), MINVAL(resolution(:,2)) |
---|
525 | WRITE(numout,*) "Resolution of the file to be interpolated (fine grid) : ", resol_lon, resol_lat |
---|
526 | WRITE(numout,*) "This interpolation assumes that we aggregate from a fine grid to a coarser grid" |
---|
527 | WRITE(numout,*) "In the data submitted it apears that the model is runing on a finer grid than the data" |
---|
528 | ENDIF |
---|
529 | ! |
---|
530 | incp = 10 |
---|
531 | IF (mini < 0.1) THEN |
---|
532 | incp=100 |
---|
533 | ELSE IF (mini < 0.01) THEN |
---|
534 | incp = 1000 |
---|
535 | ENDIF |
---|
536 | ! |
---|
537 | ! Allocate the needed memory for fine_ind |
---|
538 | ! |
---|
539 | ALLOCATE (fine_ind(NINT(domain_minlon*incp)-2:NINT(domain_maxlon*incp)+2, & |
---|
540 | & NINT(domain_minlat*incp)-2:NINT(domain_maxlat*incp)+2), STAT=ALLOC_ERR) |
---|
541 | IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_vec','ERROR IN ALLOCATION of find_ind','','') |
---|
542 | ! |
---|
543 | ! Generate a quick access table for the coarse grid |
---|
544 | ! |
---|
545 | fine_ind(:,:) = zero |
---|
546 | ! |
---|
547 | DO ib=1,nbpt |
---|
548 | coslat = MAX(COS(lalo(ib,1) * pi/180. ), mincos )*pi/180. * R_Earth |
---|
549 | ! |
---|
550 | lon_up = lalo(ib,2) + resolution(ib,1)/(2.0*coslat) |
---|
551 | lon_low =lalo(ib,2) - resolution(ib,1)/(2.0*coslat) |
---|
552 | ! |
---|
553 | coslat = pi/180. * R_Earth |
---|
554 | ! |
---|
555 | lat_up =lalo(ib,1) + resolution(ib,2)/(2.0*coslat) |
---|
556 | lat_low =lalo(ib,1) - resolution(ib,2)/(2.0*coslat) |
---|
557 | ! |
---|
558 | fine_ind(NINT(lon_low*incp):NINT(lon_up*incp),NINT(lat_low*incp):NINT(lat_up*incp))=ib |
---|
559 | ! |
---|
560 | ENDDO |
---|
561 | ! |
---|
562 | WRITE(numout,*) 'Domaine LON range : ', domain_minlon, domain_maxlon |
---|
563 | WRITE(numout,*) 'Domaine LAT range : ', domain_minlat, domain_maxlat |
---|
564 | ! |
---|
565 | ! we list a first approximation of all point we will need to |
---|
566 | ! scan to fill our coarse grid. |
---|
567 | ! |
---|
568 | IF ( domain_minlon <= -179.5 .AND. domain_maxlon >= 179.5 .AND. & |
---|
569 | & domain_minlat <= -89.5 .AND. domain_maxlat >= 89.5 ) THEN |
---|
570 | ! Here we do the entire globe |
---|
571 | nbind=0 |
---|
572 | DO ip=1,iml |
---|
573 | nbind = nbind + 1 |
---|
574 | searchind(nbind) = ip |
---|
575 | ENDDO |
---|
576 | ! |
---|
577 | ELSE |
---|
578 | ! Now we get a limited number of points |
---|
579 | nbind=0 |
---|
580 | DO ip=1,iml |
---|
581 | ! Compute the limits of the meshes of the fine grid |
---|
582 | coslat = MAX(COS(lat_rel(ip) * pi/180. ), mincos )*pi/180. * R_Earth |
---|
583 | louprel = MIN(lon_rel(ip) + resol_lon/(2.0*coslat), 180.) |
---|
584 | lolowrel = MAX(lon_rel(ip) - resol_lon/(2.0*coslat), -180.) |
---|
585 | coslat = pi/180. * R_Earth |
---|
586 | lauprel = MIN(lat_rel(ip) + resol_lat/(2.0*coslat), 90.) |
---|
587 | lalowrel = MAX(lat_rel(ip) - resol_lat/(2.0*coslat), -90.) |
---|
588 | ! |
---|
589 | IF ( louprel >= domain_minlon .AND. lolowrel <= domain_maxlon .AND.& |
---|
590 | & lauprel >= domain_minlat .AND. lalowrel <= domain_maxlat ) THEN |
---|
591 | nbind = nbind + 1 |
---|
592 | searchind(nbind) = ip |
---|
593 | ENDIF |
---|
594 | ENDDO |
---|
595 | ENDIF |
---|
596 | ! |
---|
597 | WRITE(numout,*) 'We will work with ', nbind, ' points of the fine grid and ', nbpt, 'for the coarse grid' |
---|
598 | ! |
---|
599 | WRITE(numout,*) 'Aggregate_vec : ', callsign |
---|
600 | ! |
---|
601 | ! Now we take each grid point and find out which values from the forcing we need to average |
---|
602 | ! |
---|
603 | fopt(:) = zero |
---|
604 | fopt_max = -1 |
---|
605 | ! |
---|
606 | ! |
---|
607 | ! |
---|
608 | loopnbind : DO i=1,nbind |
---|
609 | ! |
---|
610 | ! |
---|
611 | ip = searchind(i) |
---|
612 | ! |
---|
613 | ! Either the center of the data grid point is in the interval of the model grid or |
---|
614 | ! the East and West limits of the data grid point are on either sides of the border of |
---|
615 | ! the data grid. |
---|
616 | ! |
---|
617 | lonrel = lon_rel(ip) |
---|
618 | coslat = MAX(COS(lat_rel(ip) * pi/180. ), mincos )*pi/180. * R_Earth |
---|
619 | louprel = MIN(lon_rel(ip) + resol_lon/(2.0*coslat), domain_maxlon) |
---|
620 | lolowrel = MAX(lon_rel(ip) - resol_lon/(2.0*coslat), domain_minlon) |
---|
621 | ! |
---|
622 | latrel = lat_rel(ip) |
---|
623 | coslat = pi/180. * R_Earth |
---|
624 | lauprel = MIN(lat_rel(ip) + resol_lat/(2.0*coslat), domain_maxlat) |
---|
625 | lalowrel = MAX(lat_rel(ip) - resol_lat/(2.0*coslat), domain_minlat) |
---|
626 | ! |
---|
627 | ! |
---|
628 | pos_pnt(:) = zero |
---|
629 | ipp = zero |
---|
630 | pp = fine_ind(NINT(lonrel*incp),NINT(latrel*incp)) |
---|
631 | ! |
---|
632 | IF (COUNT(pos_pnt(:) == pp) == zero ) THEN |
---|
633 | pos_pnt(ipp+1) = pp |
---|
634 | ipp = ipp + 1 |
---|
635 | ENDIF |
---|
636 | pp = fine_ind(NINT(louprel*incp),NINT(lauprel*incp)) |
---|
637 | ! |
---|
638 | IF (COUNT(pos_pnt(:) == pp) == zero ) THEN |
---|
639 | pos_pnt(ipp+1) = pp |
---|
640 | ipp = ipp + 1 |
---|
641 | ENDIF |
---|
642 | pp = fine_ind(NINT(louprel*incp),NINT(lalowrel*incp)) |
---|
643 | ! |
---|
644 | IF (COUNT(pos_pnt(:) == pp) == zero ) THEN |
---|
645 | pos_pnt(ipp+1) = pp |
---|
646 | ipp = ipp + 1 |
---|
647 | ENDIF |
---|
648 | pp = fine_ind(NINT(lolowrel*incp),NINT(lauprel*incp)) |
---|
649 | ! |
---|
650 | IF (COUNT(pos_pnt(:) == pp) == zero ) THEN |
---|
651 | pos_pnt(ipp+1) = pp |
---|
652 | ipp = ipp + 1 |
---|
653 | ENDIF |
---|
654 | pp = fine_ind(NINT(lolowrel*incp),NINT(lalowrel*incp)) |
---|
655 | ! |
---|
656 | IF (COUNT(pos_pnt(:) == pp) == zero ) THEN |
---|
657 | pos_pnt(ipp+1) = pp |
---|
658 | ipp = ipp + 1 |
---|
659 | ENDIF |
---|
660 | ! |
---|
661 | ! |
---|
662 | IF ( ipp > zero ) THEN |
---|
663 | ! |
---|
664 | DO pp=1,ipp |
---|
665 | ib = pos_pnt(pp) |
---|
666 | ! |
---|
667 | ! We find the 4 limits of the grid-box. As we transform the resolution of the model |
---|
668 | ! into longitudes and latitudes we do not have the problem of periodicity. |
---|
669 | ! coslat is a help variable here ! |
---|
670 | ! |
---|
671 | coslat = MAX(COS(lalo(ib,1) * pi/180. ), mincos )*pi/180. * R_Earth |
---|
672 | ! |
---|
673 | lon_up = lalo(ib,2) + resolution(ib,1)/(2.0*coslat) |
---|
674 | lon_low =lalo(ib,2) - resolution(ib,1)/(2.0*coslat) |
---|
675 | ! |
---|
676 | coslat = pi/180. * R_Earth |
---|
677 | ! |
---|
678 | lat_up =lalo(ib,1) + resolution(ib,2)/(2.0*coslat) |
---|
679 | lat_low =lalo(ib,1) - resolution(ib,2)/(2.0*coslat) |
---|
680 | ! |
---|
681 | IF ( lonrel > lon_low .AND. lonrel < lon_up .OR. & |
---|
682 | & lolowrel < lon_low .AND. louprel > lon_low .OR. & |
---|
683 | & lolowrel < lon_up .AND. louprel > lon_up ) THEN |
---|
684 | ! |
---|
685 | ! Now that we have the longitude let us find the latitude |
---|
686 | ! |
---|
687 | IF ( latrel > lat_low .AND. latrel < lat_up .OR. & |
---|
688 | & lalowrel < lat_low .AND. lauprel > lat_low .OR.& |
---|
689 | & lalowrel < lat_up .AND. lauprel > lat_up) THEN |
---|
690 | ! |
---|
691 | fopt(ib) = fopt(ib) + 1 |
---|
692 | fopt_max = MAX ( fopt(ib), fopt_max ) |
---|
693 | ! |
---|
694 | IF ( fopt(ib) > incmax) THEN |
---|
695 | err_fopt=.TRUE. |
---|
696 | EXIT loopnbind |
---|
697 | ELSE |
---|
698 | ! |
---|
699 | ! Get the area of the fine grid in the model grid |
---|
700 | ! |
---|
701 | coslat = MAX( COS( lat_rel(ip) * pi/180. ), mincos ) |
---|
702 | ax = (MIN(lon_up,louprel)-MAX(lon_low,lolowrel))*pi/180. * R_Earth * coslat |
---|
703 | ay = (MIN(lat_up,lauprel)-MAX(lat_low,lalowrel))*pi/180. * R_Earth |
---|
704 | ! |
---|
705 | areaoverlap(ib, fopt(ib)) = ax*ay |
---|
706 | indinc(ib, fopt(ib)) = ip |
---|
707 | ! |
---|
708 | ENDIF |
---|
709 | ENDIF |
---|
710 | ENDIF |
---|
711 | ENDDO |
---|
712 | ENDIF |
---|
713 | ENDDO loopnbind |
---|
714 | ! |
---|
715 | IF (err_fopt) THEN |
---|
716 | WRITE(numout,*) 'Reached value ', fopt(ib),' for fopt on point', ib |
---|
717 | CALL ipslerr_p(2,'aggregate_vec (nbpt < nbind)', & |
---|
718 | 'Working on variable :'//callsign, & |
---|
719 | 'Reached incmax value for fopt.',& |
---|
720 | 'Please increase incmax in subroutine calling aggregate') |
---|
721 | IF (PRESENT(ok)) THEN |
---|
722 | ok = .FALSE. |
---|
723 | RETURN |
---|
724 | ELSE |
---|
725 | CALL ipslerr_p(3,'aggregate_vec','Stop now','','') |
---|
726 | ENDIF |
---|
727 | ENDIF |
---|
728 | ! |
---|
729 | WRITE(numout,*) |
---|
730 | not_found_fopt = COUNT(fopt(:) .EQ. zero) |
---|
731 | WRITE(numout,*) "aggregate_vec : ",not_found_fopt, & |
---|
732 | & "did not find any corresponding data in the input file." |
---|
733 | WRITE(numout,*) "aggregate_vec : This is ", not_found_fopt/FLOAT(nbpt)*100., & |
---|
734 | & " % of the grid" |
---|
735 | WRITE(numout,*) "aggregate_vec : nbvmax = ",incmax, "max used = ",fopt_max |
---|
736 | ! |
---|
737 | ! Do some memory management. |
---|
738 | ! |
---|
739 | DEALLOCATE (searchind) |
---|
740 | DEALLOCATE (fine_ind) |
---|
741 | ! |
---|
742 | ! Close the progress meter |
---|
743 | ! |
---|
744 | WRITE(numout,*) ' ' |
---|
745 | ! |
---|
746 | END SUBROUTINE aggregate_vec |
---|
747 | ! |
---|
748 | ! |
---|
749 | |
---|
750 | SUBROUTINE aggregate_vec_p(nbpt, lalo, neighbours, resolution, contfrac, & |
---|
751 | & iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign, & |
---|
752 | & nbvmax, sub_index, sub_area, ok) |
---|
753 | |
---|
754 | IMPLICIT NONE |
---|
755 | |
---|
756 | INTEGER(i_std), INTENT(in) :: nbpt |
---|
757 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) |
---|
758 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) |
---|
759 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) |
---|
760 | REAL(r_std), INTENT(in) :: contfrac(nbpt) |
---|
761 | INTEGER(i_std), INTENT(in) :: iml |
---|
762 | REAL(r_std), INTENT(in) :: lon_ful(iml) |
---|
763 | REAL(r_std), INTENT(in) :: lat_ful(iml) |
---|
764 | REAL(r_std), INTENT(in) :: resol_lon, resol_lat |
---|
765 | CHARACTER(LEN=*), INTENT(in) :: callsign |
---|
766 | INTEGER(i_std), INTENT(in) :: nbvmax |
---|
767 | INTEGER(i_std), INTENT(out) :: sub_index(nbpt,nbvmax) |
---|
768 | REAL(r_std), INTENT(out) :: sub_area(nbpt,nbvmax) |
---|
769 | LOGICAL, OPTIONAL, INTENT(out) :: ok ! return code |
---|
770 | |
---|
771 | INTEGER(i_std) :: sub_index_g(nbp_glo,nbvmax) |
---|
772 | REAL(r_std) :: sub_area_g(nbp_glo,nbvmax) |
---|
773 | |
---|
774 | IF (is_root_prc) CALL aggregate(nbp_glo, lalo_g, neighbours_g, resolution_g, contfrac_g, & |
---|
775 | & iml, lon_ful, lat_ful, resol_lon, resol_lat, callsign, & |
---|
776 | & nbvmax, sub_index_g, sub_area_g, ok) |
---|
777 | |
---|
778 | CALL BCAST(ok) |
---|
779 | CALL scatter(sub_index_g,sub_index) |
---|
780 | CALL scatter(sub_area_g,sub_area) |
---|
781 | |
---|
782 | |
---|
783 | END SUBROUTINE aggregate_vec_p |
---|
784 | |
---|
785 | SUBROUTINE aggregate_2d_p(nbpt, lalo, neighbours, resolution, contfrac, & |
---|
786 | & iml, jml, lon_ful, lat_ful, mask, callsign, & |
---|
787 | & nbvmax, sub_index, sub_area, ok) |
---|
788 | |
---|
789 | IMPLICIT NONE |
---|
790 | |
---|
791 | INTEGER(i_std), INTENT(in) :: nbpt |
---|
792 | REAL(r_std), INTENT(in) :: lalo(nbpt,2) |
---|
793 | INTEGER(i_std), INTENT(in) :: neighbours(nbpt,8) |
---|
794 | REAL(r_std), INTENT(in) :: resolution(nbpt,2) |
---|
795 | REAL(r_std), INTENT(in) :: contfrac(nbpt) |
---|
796 | INTEGER(i_std), INTENT(in) :: iml,jml |
---|
797 | REAL(r_std), INTENT(in) :: lon_ful(iml,jml) |
---|
798 | REAL(r_std), INTENT(in) :: lat_ful(iml,jml) |
---|
799 | INTEGER(i_std), INTENT(in) :: mask(iml, jml) |
---|
800 | CHARACTER(LEN=*), INTENT(in) :: callsign |
---|
801 | INTEGER(i_std), INTENT(in) :: nbvmax |
---|
802 | INTEGER(i_std), INTENT(out) :: sub_index(nbpt,nbvmax,2) |
---|
803 | REAL(r_std), INTENT(out) :: sub_area(nbpt,nbvmax) |
---|
804 | LOGICAL, OPTIONAL, INTENT(out) :: ok ! return code |
---|
805 | |
---|
806 | INTEGER(i_std) :: sub_index_g(nbp_glo,nbvmax,2) |
---|
807 | REAL(r_std) :: sub_area_g(nbp_glo,nbvmax) |
---|
808 | |
---|
809 | IF (is_root_prc) CALL aggregate_2d(nbp_glo, lalo_g, neighbours_g, resolution_g, contfrac_g, & |
---|
810 | & iml, jml, lon_ful, lat_ful, mask, callsign, & |
---|
811 | & nbvmax, sub_index_g, sub_area_g, ok) |
---|
812 | CALL BCAST(ok) |
---|
813 | CALL scatter(sub_index_g,sub_index) |
---|
814 | CALL scatter(sub_area_g,sub_area) |
---|
815 | |
---|
816 | END SUBROUTINE aggregate_2d_p |
---|
817 | ! |
---|
818 | END MODULE interpol_help |
---|