source: branches/ORCHIDEE_EXT/ORCHIDEE/src_global/interpol_help.f90 @ 64

Last change on this file since 64 was 64, checked in by didier.solyga, 13 years ago

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