source: branches/ORCHIDEE_2_2/ORCHIDEE/src_global/interpol_help.f90 @ 7259

Last change on this file since 7259 was 7259, checked in by agnes.ducharne, 3 years ago

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