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