source: branches/publications/ORCHIDEE_CAN_r3069/src_global/interpol_help.f90 @ 7346

Last change on this file since 7346 was 2442, checked in by sebastiaan.luyssaert, 10 years ago

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