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

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

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