source: branches/publications/ORCHIDEE_gmd_2018_MICT-LEAK/src_global/interregxy.f90 @ 6892

Last change on this file since 6892 was 3564, checked in by albert.jornet, 9 years ago

Merge: from [3313:3545/trunk/ORCHIDEE]
Clean: output subroutine variables compile warning messages are solved.

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