source: branches/publications/ORCHIDEE-ICE_SurfaceMassBalance/src_global/interpol_help.f90

Last change on this file was 4290, checked in by josefine.ghattas, 7 years ago

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