source: perso/abdelouhab.djerrah/ORCHIDEE_ORIG/src_global/interpol_help.f90 @ 938

Last change on this file since 938 was 41, checked in by mmaipsl, 14 years ago

MM, MasaK : Group all definitions of R_Earth in the whole ORCHIDEE code : use a global

definition in constantes.f90

MM, NVui, DS: inverse compilation order of src_parallel and src_parameters

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