source: branches/ORCHIDEE_2_2/ORCHIDEE/src_global/interregxy.f90 @ 7796

Last change on this file since 7796 was 7597, checked in by josefine.ghattas, 3 years ago

Integrated bug correction done ORCHIDEE-ROUTING [3608] by Jan Polcher:
This interpolation did not consider the mask provided to aggregate. Needed for interpolation of Safran forcing for exemple.
File size: 19.9 KB
Line 
1MODULE interregxy
2
3  ! Modules used :
4
5  USE constantes
6  USE mod_orchidee_para
7  USE grid
8  USE module_llxy
9  USE polygones
10
11  IMPLICIT NONE
12
13  PRIVATE
14  PUBLIC interregxy_aggr2d, interregxy_aggrve
15
16CONTAINS
17
18  SUBROUTINE interregxy_aggr2d(nbpt, lalo, neighbours, resolution, contfrac, &
19       &                     iml, jml, lon_rel, lat_rel, mask, callsign, &
20       &                     incmax, indinc, areaoverlap, ok)
21    !
22    ! INPUT
23    !
24    INTEGER(i_std), INTENT(in)   :: nbpt                 ! Number of points for which the data needs to be interpolated
25    REAL(r_std), INTENT(in)      :: lalo(nbpt,2)         ! Vector of latitude and longitudes (beware of the order !)
26    INTEGER(i_std), INTENT(in)   :: neighbours(nbpt,NbNeighb)! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W)
27    REAL(r_std), INTENT(in)      :: resolution(nbpt,2)   ! The size in km of each grid-box in X and Y
28    REAL(r_std), INTENT(in)      :: contfrac(nbpt)       ! Fraction of land in each grid box.
29    INTEGER(i_std), INTENT(in)   :: iml, jml             ! Size of the source grid
30    REAL(r_std), INTENT(in)      :: lon_rel(iml, jml)    ! Longitudes for the source grid
31    REAL(r_std), INTENT(in)      :: lat_rel(iml, jml)    ! Latitudes for the souce grid
32    INTEGER(i_std), INTENT(in)   :: mask(iml, jml)       ! Mask which retains only the significative points
33    ! of the source grid.
34    CHARACTER(LEN=*), INTENT(in) :: callsign             ! Allows to specify which variable is beeing treated
35    INTEGER(i_std), INTENT(in)   :: incmax               ! Maximum point of the source grid we can store.
36    !
37    ! Output
38    !
39    INTEGER(i_std), INTENT(out)    :: indinc(nbpt,incmax,2)
40    REAL(r_std), INTENT(out)       :: areaoverlap(nbpt,incmax)
41    LOGICAL, OPTIONAL, INTENT(out) :: ok                 ! return code
42    !
43    ! Local Variables
44    !
45    ! The number of sub-division we wish to have on each side of the rectangle.
46    ! sbd=1 means that on each side we take only the middle point.
47    !
48    INTEGER(i_std), PARAMETER  :: sbd=5
49    INTEGER(i_std), PARAMETER  :: idom=1
50    !
51    INTEGER(i_std)                                :: i, is, js, in
52    INTEGER(i_std)                                :: ilow, iup, jlow, jup
53    REAL(r_std)                                   :: dle, dlw, dls, dln
54    !
55    REAL(r_std), DIMENSION((sbd+1)*4+1)           :: lons, lats, ri, rj
56    !
57    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:)   :: nbov
58    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sourcei, sourcej
59    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)    :: overlap
60    LOGICAL                                       :: checkprint=.FALSE.
61    !
62    !
63    REAL(r_std), DIMENSION(2)   :: BBlon, BBlat
64    !
65    WRITE(*,*) "interregxy_aggr2d working on ", callsign
66    !
67    ALLOCATE(sourcei(iim_g, jjm_g, incmax))
68    ALLOCATE(sourcej(iim_g, jjm_g, incmax))
69    ALLOCATE(nbov(iim_g, jjm_g))
70    ALLOCATE(overlap(iim_g, jjm_g, incmax))
71    nbov(:,:) = 0
72    !
73    BBlon(1) = MINVAL(lalo(:,2))
74    BBlon(2) = MAXVAL(lalo(:,2))
75    BBlat(1) = MINVAL(lalo(:,1))
76    BBlat(2) = MAXVAL(lalo(:,1))
77    !
78    ! Go through the entire source grid and try to place each grid box in a target grid box.
79    !
80    DO is=1,iml
81       DO js=1,jml
82          !
83          IF ( mask(is,js) == 1 ) THEN
84             !
85             ! 1.0 Get the center of the box and their corresponding grid point
86             !
87             lons(1) = lon_rel(is,js)
88             lats(1) = lat_rel(is,js)
89             !
90             ! Assume all grids go from West to Est
91             dle=(lon_rel(MIN(is+1,iml),js)-lon_rel(is,js))/2.0
92             dlw=(lon_rel(is,js)-lon_rel(MAX(is-1,1),js))/2.0
93             !
94             IF ( lat_rel(1,1) < lat_rel(iml,jml)) THEN
95                ! Grid in SN direction
96                dls=(lat_rel(is,js)-lat_rel(is,MAX(js-1,1)))/2.0
97                dln=(lat_rel(is,MIN(js+1,jml))-lat_rel(is,js))/2.0
98             ELSE
99                dls=(lat_rel(is,js)-lat_rel(is,MIN(js+1,jml)))/2.0
100                dln=(lat_rel(is,MAX(js-1,1))-lat_rel(is,js))/2.0
101             ENDIF
102             ! Treat the border cases
103             IF ( dlw < dle/2. ) dlw = dle
104             IF ( dle < dlw/2. ) dle = dlw
105             !
106             IF ( dls < dln/2. ) dls = dln
107             IF ( dln < dls/2. ) dln = dls
108             !
109             DO i=0,sbd
110                ! Put the points in the order : SW, SE, NE, NW
111                lons(2+i) = lons(1) - dlw + i*((dlw+dle)/FLOAT(sbd+1))
112                lats(2+i) = lats(1) - dls
113                !
114                lons(2+(sbd+1)+i) = lons(1) + dle
115                lats(2+(sbd+1)+i) = lats(1) - dls + i*((dls+dln)/FLOAT(sbd+1))
116                !
117                lons(2+2*(sbd+1)+i) = lons(1) + dle - i*((dlw+dle)/FLOAT(sbd+1))
118                lats(2+2*(sbd+1)+i) = lats(1) + dln
119                !
120                lons(2+3*(sbd+1)+i) = lons(1) - dlw
121                lats(2+3*(sbd+1)+i) = lats(1) + dln - i*((dls+dln)/FLOAT(sbd+1))
122             ENDDO
123             !
124             ! Test that the polygone overlaps with the Bounding box. We might have projections
125             ! which are not valid outside of the RCM domain.
126             !
127             IF ( ((MINVAL(lons) > BBlon(1) .AND. MINVAL(lons) < BBlon(2)) .OR.   &
128                  &  (MAXVAL(lons) > BBlon(1) .AND. MAXVAL(lons) < BBlon(2))) .AND. & 
129                  & ((MINVAL(lats) > BBlat(1) .AND. MINVAL(lats) < BBlat(2)) .OR.   &
130                  &  (MAXVAL(lats) > BBlat(1) .AND. MAXVAL(lats) < BBlat(2)))) THEN
131                !
132                !
133                CALL Grid_Toij (lons, lats, ri, rj)
134                !
135                !
136                ! Find the points of the WRF grid boxes we need to scan to place this box (is, js)
137                ! of the source grid.
138                !
139                ilow = MAX(FLOOR(MINVAL(ri)),1)
140                iup  = MIN(CEILING(MAXVAL(ri)), iim_g)
141                jlow = MAX(FLOOR(MINVAL(rj)),1)
142                jup  = MIN(CEILING(MAXVAL(rj)), jjm_g)
143                !
144                !
145                IF ( ilow < iup .OR. jlow < jup ) THEN
146                   !
147                   ! Find the grid boxes in WRF and compute the overlap area.
148                   !
149                   CALL interregxy_findpoints(iim_g, jjm_g, lon_rel(is,js), lon_rel(is,js), is, js, ri, rj, incmax, &
150                        &    checkprint, nbov, sourcei, sourcej, overlap)
151
152                ENDIF
153                !
154                !
155             ENDIF
156          ENDIF
157       ENDDO
158    ENDDO
159    !
160    ! Gather the land points from the global grid
161    !
162    areaoverlap(:,:) = zero
163    DO in=1,nbpt
164       !
165       is = ilandindex(in)
166       js = jlandindex(in)
167       !
168       indinc(in,1:nbov(is,js),1) = sourcei(is,js,1:nbov(is,js))
169       indinc(in,1:nbov(is,js),2) = sourcej(is,js,1:nbov(is,js))
170       areaoverlap(in,1:nbov(is,js)) = overlap(is,js,1:nbov(is,js))
171    ENDDO
172    !
173    ! The model will stop if incmax is not sufficient. So if we are here all is OK.
174    !
175    ok=.TRUE.
176    !
177  END SUBROUTINE interregxy_aggr2d
178!
179! ===================================================================================
180!
181  SUBROUTINE interregxy_aggrve(nbpt, lalo, neighbours, resolution, contfrac, &
182       &                iml, lon_rel, lat_rel, resol_lon, resol_lat, callsign, &
183       &                incmax, indinc, areaoverlap, ok)
184    !
185    ! INPUT
186    !
187    INTEGER(i_std), INTENT(in)   :: nbpt                 ! Number of points for which the data needs to be interpolated
188    REAL(r_std), INTENT(in)      :: lalo(nbpt,2)         ! Vector of latitude and longitudes (beware of the order !)
189    INTEGER(i_std), INTENT(in)   :: neighbours(nbpt,NbNeighb)   ! Vector of neighbours for each grid point (1=N, 2=E, 3=S, 4=W)
190    REAL(r_std), INTENT(in)      :: resolution(nbpt,2)   ! The size in km of each grid-box in X and Y
191    REAL(r_std), INTENT(in)      :: contfrac(nbpt)       ! Fraction of land in each grid box.
192    INTEGER(i_std), INTENT(in)   :: iml                  ! Size of the source grid
193    REAL(r_std), INTENT(in)      :: lon_rel(iml)         ! Longitudes for the source grid
194    REAL(r_std), INTENT(in)      :: lat_rel(iml)         ! Latitudes for the source grid
195    REAL(r_std), INTENT(in)      :: resol_lon, resol_lat ! Resolution in meters of the source grid
196    CHARACTER(LEN=*), INTENT(in) :: callsign             ! Allows to specify which variable is beeing treated
197    INTEGER(i_std), INTENT(in)   :: incmax               ! Maximum point of the source grid we can store.
198    !
199    ! Output
200    !
201    INTEGER(i_std), INTENT(out)    :: indinc(nbpt,incmax)
202    REAL(r_std), INTENT(out)       :: areaoverlap(nbpt,incmax)
203    LOGICAL, OPTIONAL, INTENT(out) :: ok                 ! return code
204    !
205    !
206    ! Local Variables
207    !
208    ! The number of sub-division we wish to have on each side of the rectangle.
209    ! sbd=1 means that on each side we take only the middle point.
210    !
211    INTEGER(i_std), PARAMETER  :: sbd=5
212    INTEGER(i_std), PARAMETER  :: idom=1
213    !
214    INTEGER(i_std)                                :: i, is, js, in
215    INTEGER(i_std)                                :: ilow, iup, jlow, jup
216    REAL(r_std)                                   :: dle, dlw, dls, dln, coslat
217    !
218    REAL(r_std), DIMENSION((sbd+1)*4+1)           :: lons, lats, ri, rj
219    !
220    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:)   :: nbov
221    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sourcei, sourcej
222    REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)    :: overlap
223    !
224    !
225    REAL(r_std), DIMENSION(2)   :: BBlon, BBlat
226    REAL(r_std)                 :: half_resol_lon, half_resol_lat, trtmp
227    LOGICAL                     :: checkprint=.FALSE.
228    !
229    WRITE(*,*) "interregxy_aggrve working on ", callsign
230    !
231    ALLOCATE(sourcei(iim_g, jjm_g, incmax))
232    ALLOCATE(sourcej(iim_g, jjm_g, incmax))
233    ALLOCATE(nbov(iim_g, jjm_g))
234    ALLOCATE(overlap(iim_g, jjm_g, incmax))
235    nbov(:,:) = 0
236    !
237    BBlon(1) = MINVAL(lalo(:,2))
238    BBlon(2) = MAXVAL(lalo(:,2))
239    BBlat(1) = MINVAL(lalo(:,1))
240    BBlat(2) = MAXVAL(lalo(:,1))
241    !
242    ! Go through the entire source grid and try to place each grid box in a target grid box.
243    !
244    DO is=1,iml
245       !
246       ! 1.0 Get the center of the box and their corresponding grid point
247       !
248       lons(1) = lon_rel(is)
249       lats(1) = lat_rel(is)
250       ! Some pre-calculations
251       trtmp = 180./(2.*pi)
252       half_resol_lon = resol_lon/2.0
253       half_resol_lat = resol_lat/2.0
254       !
255       dln = half_resol_lat*trtmp/R_Earth
256       dls = half_resol_lat*trtmp/R_Earth
257       !
258       DO i=0,sbd
259          ! Put the points in the order : SW, SE, NE, NW
260          lats(2+i) = lats(1) - dls
261          coslat = COS( lats(2+i) * pi/180. )
262          dlw = half_resol_lon*trtmp/R_Earth/coslat
263          lons(2+i) = lons(1) - dlw + i*((2.*dlw)/FLOAT(sbd+1))
264          !
265          lats(2+(sbd+1)+i) = lats(1) - dls + i*((dls+dln)/FLOAT(sbd+1))
266          coslat = COS( lats(2+(sbd+1)+i) * pi/180. )
267          dle = half_resol_lon*trtmp/R_Earth/coslat
268          lons(2+(sbd+1)+i) = lons(1) + dle
269          !
270          lats(2+2*(sbd+1)+i) = lats(1) + dln
271          coslat = COS( lats(2+2*(sbd+1)+i) * pi/180. )
272          dle = half_resol_lon*trtmp/R_Earth/coslat
273          lons(2+2*(sbd+1)+i) = lons(1) + dle - i*((2*dle)/FLOAT(sbd+1))
274          !
275          lats(2+3*(sbd+1)+i) = lats(1) + dln - i*((dls+dln)/FLOAT(sbd+1))
276          coslat = COS( lats(2+3*(sbd+1)+i) * pi/180. )
277          dlw = half_resol_lon*trtmp/R_Earth/coslat
278          lons(2+3*(sbd+1)+i) = lons(1) - dlw
279       ENDDO
280       !
281       ! Test that the polygone overlaps with the Bounding box. We might have projections
282       ! which are not valid outside of the RCM domain.
283       !
284       IF ( ((MINVAL(lons) > BBlon(1) .AND. MINVAL(lons) < BBlon(2)) .OR.   &
285            &  (MAXVAL(lons) > BBlon(1) .AND. MAXVAL(lons) < BBlon(2))) .AND. & 
286            & ((MINVAL(lats) > BBlat(1) .AND. MINVAL(lats) < BBlat(2)) .OR.   &
287            &  (MAXVAL(lats) > BBlat(1) .AND. MAXVAL(lats) < BBlat(2)))) THEN
288          !
289          CALL Grid_Toij (lons, lats, ri, rj)
290          !
291          !
292          ! Find the points of the WRF grid boxes we need to scan to place this box (is, js)
293          ! of the source grid.
294          !
295          ilow = MAX(FLOOR(MINVAL(ri)),1)
296          iup  = MIN(CEILING(MAXVAL(ri)), iim_g)
297          jlow = MAX(FLOOR(MINVAL(rj)),1)
298          jup  = MIN(CEILING(MAXVAL(rj)), jjm_g)
299          !
300          !
301          IF ( ilow < iup .OR. jlow < jup ) THEN
302             !
303             ! Find the grid boxes in WRF and compute the overlap area.
304             !
305             CALL interregxy_findpoints(iim_g, jjm_g, lon_rel(is), lon_rel(is), is, 1, ri, rj, incmax, &
306                  &    checkprint, nbov, sourcei, sourcej, overlap)
307             
308          ENDIF
309          !
310          !
311       ENDIF
312    ENDDO
313    !
314    ! Gather the land points from the global grid
315    !
316    areaoverlap(:,:) = zero
317    DO in=1,nbpt
318       !
319       is = ilandindex(in)
320       js = jlandindex(in)
321       !
322       indinc(in,1:nbov(is,js)) = sourcei(is,js,1:nbov(is,js))
323       areaoverlap(in,1:nbov(is,js)) = overlap(is,js,1:nbov(is,js))
324    ENDDO
325    !
326    ! The model will stop if incmax is not sufficient. So if we are here all is OK.
327    !
328    ok=.TRUE.
329    !
330  END SUBROUTINE interregxy_aggrve
331!
332! ===================================================================================
333!
334  SUBROUTINE interregxy_findpoints(nb_we, nb_sn, lon_cen, lat_cen, is, js, ri, rj, &
335       &                           nbpt_max, checkprint, nbpt, sourcei, sourcej, overlap_area)
336    !
337    !
338    !   This Subroutine finds the overlap of grid box ri,rj with the WRF grid.  (ri,rj) are the WRF indexes of the (is,js) point
339    !   of the source grid.
340    !   Each time we find an overlap of (is,js) with a WRF grid box (ix,jx) we add that in the (sourcei,sourcej) fields.
341    !
342    !
343    ! Arguments
344    !
345    INTEGER(i_std), INTENT(in)                      :: nb_we, nb_sn
346    REAL(r_std), INTENT(in)                         :: lon_cen, lat_cen
347    INTEGER(i_std), INTENT(in)                      :: is, js
348    REAL(r_std), DIMENSION(:), INTENT(in)           :: ri, rj
349    INTEGER(i_std), INTENT(in)                      :: nbpt_max
350    LOGICAL, INTENT(in)                             :: checkprint
351    INTEGER(i_std), DIMENSION(:,:), INTENT(out)     :: nbpt
352    INTEGER(i_std), DIMENSION(:,:,:), INTENT(out)   :: sourcei, sourcej
353    REAL(r_std), DIMENSION(:,:,:), INTENT(out)      :: overlap_area
354    !
355    ! Local
356    !
357    INTEGER(i_std)                                  :: ilow, iup, jlow, jup, nbr
358    INTEGER(i_std)                                  :: ix, jx, i
359
360    INTEGER(i_std)                                  :: nbint, nbcross, nbsorted, nbfinal, nbtmpa, nbtmpb
361    !
362    INTEGER(i_std)                                  :: nbdots=5, dimpoly = 200
363    REAL(r_std), DIMENSION(4,2)                     :: poly_wrf
364    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)        :: poly_in, poly_tmpa, poly_tmpb, poly_int, poly_cross
365    !
366    REAL(r_std)                                     :: ax, ay, overlaparea
367    !
368    INTEGER(i_std), PARAMETER                       :: icheck=50, jcheck=46
369    !
370    !
371    !
372    nbr = SIZE(ri)
373    !
374    IF ( dimpoly < 4*(nbr-1) ) THEN
375       dimpoly =  4*nbr
376    ENDIF
377    !
378    ALLOCATE(poly_in(nbr-1,2))
379    ALLOCATE(poly_tmpa(dimpoly,2))
380    ALLOCATE(poly_tmpb(dimpoly,2))
381    ALLOCATE(poly_int(dimpoly,2))
382    ALLOCATE(poly_cross(dimpoly,2))
383    !
384    ! Scan all Determine the range of indicis of the WRF grid we need to scan
385    ! within the box of the source grid.
386    !
387    ilow = MAX(FLOOR(MINVAL(ri)),1)
388    iup  = MIN(CEILING(MAXVAL(ri)), nb_we)
389    jlow = MAX(FLOOR(MINVAL(rj)),1)
390    jup  = MIN(CEILING(MAXVAL(rj)), nb_sn)
391    !
392    !
393    poly_wrf(:,1) = (/ -0.5, 0.5, 0.5, -0.5 /)
394    poly_wrf(:,2) = (/ -0.5, -0.5, 0.5, 0.5 /)
395    !
396    DO ix = ilow,iup
397       DO jx = jlow,jup
398          !
399          ! Bring the points of the source grid into the +-0.5 range of the WRF grid ...
400          ! and remove the center point.
401          !
402          poly_in(1:nbr-1,1) = ri(2:nbr) - ix
403          poly_in(1:nbr-1,2) = rj(2:nbr) - jx
404          !
405          IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN
406             WRITE(*,*) "X = ", poly_in(1:nbr-1,1)
407             WRITE(*,*) "Y = ", poly_in(1:nbr-1,2)
408          ENDIF
409          !
410          CALL polygones_extend(nbr-1, poly_in, nbdots, nbtmpa, poly_tmpa)
411          CALL polygones_extend(4, poly_wrf, nbdots, nbtmpb, poly_tmpb)
412          !
413          IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN
414             WRITE(*,*) "X_extend = ", poly_tmpa(1:nbtmpa,1)
415             WRITE(*,*) "Y_extend = ", poly_tmpa(1:nbtmpa,2)
416          ENDIF
417          !
418          CALL polygones_intersection(nbtmpa, poly_tmpa, nbtmpb, poly_tmpb, nbint, poly_int) 
419          CALL polygones_crossing(nbr-1, poly_in,  4, poly_wrf, nbcross, poly_cross)
420          !
421          IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN
422              WRITE(*,*) "nbint = nbcross = ", nbint, nbcross
423           ENDIF
424          !
425          IF ( nbint+nbcross > dimpoly ) THEN
426             WRITE(*,*) "The intersection of polygones is larger than the memory foressen : ",  nbint+nbcross, dimpoly
427             STOP "interregxy_compoverlap"
428          ENDIF
429          !
430          !
431          IF ( nbint+nbcross > 2 ) THEN
432
433             poly_tmpa(1:nbint,:) = poly_int(1:nbint,:)
434             poly_tmpa(nbint+1:nbint+nbcross,:) =  poly_cross(1:nbcross,:)
435             
436             IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN
437                WRITE(*,*) "X_overlap = ", poly_tmpa(1:(nbint+nbcross),1)
438                WRITE(*,*) "Y_overlap = ", poly_tmpa(1:(nbint+nbcross),2)
439             ENDIF
440
441             CALL polygones_cleanup(nbint+nbcross, poly_tmpa, nbfinal, poly_tmpb)
442             IF ( nbint+nbcross < 3 ) THEN
443                WRITE(*,*) "poly_tmpa X : ", poly_tmpa(1:(nbint+nbcross),1)
444                WRITE(*,*) "poly_tmpa Y : ", poly_tmpa(1:(nbint+nbcross),2)
445                WRITE(*,*) "Cleaning goes too far : ", nbint+nbcross, " to ", nbfinal
446                WRITE(*,*) "poly_tmpb X : ", poly_tmpb(1:nbfinal,1)
447                WRITE(*,*) "poly_tmpb Y : ", poly_tmpb(1:nbfinal,2)
448             ENDIF
449
450             IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN
451                WRITE(*,*) "X_final = ", poly_tmpb(1:nbfinal,1)
452                WRITE(*,*) "Y_final = ", poly_tmpb(1:nbfinal,2)
453             ENDIF
454             IF ( nbfinal > 2 ) THEN
455                CALL polygones_convexhull(nbfinal, poly_tmpb, nbsorted, poly_tmpa)
456                IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN
457                   WRITE(*,*) "X_sorted = ", poly_tmpa(1:nbsorted,1)
458                   WRITE(*,*) "Y_sorted = ", poly_tmpa(1:nbsorted,2)
459                ENDIF
460                CALL polygones_area(nbsorted, poly_tmpa, dxWRF(ix,jx), dyWRF(ix,jx), overlaparea)
461             ELSE
462                overlaparea = 0.0
463             ENDIF
464
465             IF ( checkprint .AND. ix == icheck .AND. jx == jcheck ) THEN
466                WRITE(*,*) ix, jx, "overlaparea = ", overlaparea
467                WRITE(*,*) "============================================================================="
468             ENDIF
469
470             IF ( overlaparea > 0.0 ) THEN
471                nbpt(ix,jx) = nbpt(ix,jx)+1
472                IF ( nbpt(ix,jx) > nbpt_max ) THEN
473                   WRITE(*,*) "ERROR in interregxy_findpoints. "
474                   WRITE(*,*) "Consider your algorithm to estimate nbpt_max. "
475                   WRITE(*,*) "You should change it as it does not give enough points."
476                   WRITE(*,*) "nbpt(ix,jx) = ", nbpt(ix,jx)
477                   WRITE(*,*) "nbpt_max = ", nbpt_max
478                   STOP "interregxy_findpoints"
479                ENDIF
480                sourcei(ix,jx,nbpt(ix,jx)) = is
481                sourcej(ix,jx,nbpt(ix,jx)) = js
482                overlap_area(ix,jx,nbpt(ix,jx)) =  overlaparea
483                !
484             ENDIF
485          ENDIF
486       ENDDO
487    ENDDO
488    !
489    DEALLOCATE(poly_in, poly_tmpa, poly_tmpb, poly_int, poly_cross)
490    !
491  END SUBROUTINE interregxy_findpoints
492
493END MODULE interregxy
Note: See TracBrowser for help on using the repository browser.