source: tags/ORCHIDEE/src_global/interpol_help.f90 @ 6

Last change on this file since 6 was 6, checked in by orchidee, 14 years ago

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