source: lmdz_wrf/trunk/tools/module_scientific.f90 @ 1614

Last change on this file since 1614 was 1614, checked in by lfita, 8 years ago

Fixing last issues on polygons detection
Adding direct acces to the 'module_scientific.f90' as `module_ForSci'
File size: 31.0 KB
Line 
1MODULE module_scientific
2! Module of the scientific function/subroutines
3 
4!!!!!!! Functions & subroutines
5! borders_matrixL: Subroutine to provide the borders of a logical array (interested in .TRUE.)
6! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End.
7! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon
8!   following ray casting algorithm
9! look_clockwise_borders: Subroutine to look clock-wise for a next point within a collection of borders
10!   (limits of a region)
11! paths_border: Subroutine to search the paths of a border field.
12! path_properties: Subroutine to determine the properties of a path
13! polygons: Subroutine to search the polygons of a border field. FORTRAN based. 1st = 1!
14! polygons_t: Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!
15! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k)
16! quantilesR_K: Subroutine to provide the quantiles of a given set of values of type real 'r_k'
17! rand_sample: Subroutine to randomly sample a range of indices
18! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order.
19! SwapR_K*: Subroutine swaps the values of its two formal arguments.
20
21!!! *Functions/Subroutines to sort values adpated. The method used is usually referred to as "selection" method.
22! from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
23
24  USE module_definitions
25  USE module_generic
26
27  CONTAINS
28
29SUBROUTINE polygons_t(dbg, dx, dy, dt, boolmatt, polys, Npoly)
30! Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1!
31
32  IMPLICIT NONE
33
34  INTEGER, INTENT(in)                                    :: dx, dy, dt
35  LOGICAL, DIMENSION(dx,dy,dt), INTENT(in)               :: boolmatt
36  LOGICAL, INTENT(in)                                    :: dbg
37  INTEGER, DIMENSION(dt), INTENT(out)                    :: Npoly
38  INTEGER, DIMENSION(dx,dy,dt), INTENT(out)              :: polys
39
40! Local
41  INTEGER                                                :: i,it
42
43!!!!!!! Variables
44! dx,dy: spatial dimensions of the space
45! boolmatt: boolean matrix tolook for the polygons (.TRUE. based)
46! polys: found polygons
47! Npoly: number of polygons found
48
49  fname = 'polygons'
50
51  IF (dbg) PRINT *,TRIM(fname)
52
53  polys = -1
54  Npoly = 0
55
56  DO it=1,dt
57    IF (ANY(boolmatt(:,:,it))) THEN
58      IF (dbg) THEN
59        PRINT *,'  it:', it, ' num. TRUE:', COUNT(boolmatt(:,:,it)), 'bool _______'
60        DO i=1,dx
61          PRINT *,boolmatt(i,:,it)
62        END DO
63      END IF
64      CALL polygons(dbg, dx, dy, boolmatt(:,:,it), polys(:,:,it), Npoly(it))
65    ELSE
66      IF (dbg) THEN
67        PRINT *,'  it:', it, " without '.TRUE.' values skipiing it!!"
68      END IF
69    END IF
70  END DO
71
72END SUBROUTINE polygons_t
73
74SUBROUTINE polygons(dbg, dx, dy, boolmat, polys, Npoly)
75! Subroutine to search the polygons of a boolean field. FORTRAN based. 1st = 1!
76
77  IMPLICIT NONE
78
79  INTEGER, INTENT(in)                                    :: dx, dy
80  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: boolmat
81  LOGICAL, INTENT(in)                                    :: dbg
82  INTEGER, INTENT(out)                                   :: Npoly
83  INTEGER, DIMENSION(dx,dy), INTENT(out)                 :: polys
84
85! Local
86  INTEGER                                                :: i, j, ip, ipp, Nppt
87  INTEGER                                                :: ierr
88  INTEGER, DIMENSION(:,:), ALLOCATABLE                   :: borders
89  LOGICAL, DIMENSION(dx,dy)                              :: isborder, isbordery
90  INTEGER, DIMENSION(:,:,:), ALLOCATABLE                 :: paths
91  INTEGER                                                :: Npath
92  INTEGER, DIMENSION(:), ALLOCATABLE                     :: Nptpaths
93  INTEGER, DIMENSION(2)                                  :: xtrx, xtry, meanpth
94  INTEGER                                                :: Nvertx, Npts
95  INTEGER, DIMENSION(:,:), ALLOCATABLE                   :: vertxs, points
96  LOGICAL, DIMENSION(:), ALLOCATABLE                     :: isin
97
98!!!!!!! Variables
99! dx,dy: spatial dimensions of the space
100! boolmat: boolean matrix tolook for the polygons (.TRUE. based)
101! polys: found polygons
102! Npoly: number of polygons found
103
104  fname = 'polygons'
105
106  polys = -1
107
108  Nppt = dx*dy/10
109  Npts = dx*dy
110
111  IF (ALLOCATED(borders)) DEALLOCATE(borders)
112  ALLOCATE(borders(Nppt,2), STAT=ierr)
113  msg = "Problems allocating matrix 'borders'"
114  CALL ErrMsg(msg, fname, ierr)
115
116  IF (ALLOCATED(paths)) DEALLOCATE(paths)
117  ALLOCATE(paths(Nppt,Nppt,2), STAT=ierr)
118  msg = "Problems allocating matrix 'paths'"
119  CALL ErrMsg(msg, fname, ierr)
120
121  IF (ALLOCATED(Nptpaths)) DEALLOCATE(Nptpaths)
122  ALLOCATE(Nptpaths(Nppt), STAT=ierr)
123  msg = "Problems allocating matrix 'Nptpaths'"
124  CALL ErrMsg(msg, fname, ierr)
125
126  ! Filling with the points of all the space
127  IF (ALLOCATED(points)) DEALLOCATE(points)
128  ALLOCATE(points(Npts,2), STAT=ierr)
129  msg = "Problems allocating matrix 'points'"
130  CALL ErrMsg(msg, fname, ierr)
131
132  ip = 1
133  DO i=1, dx
134    DO j=1, dy
135      points(ip,1) = i
136      points(ip,2) = j
137      ip = ip + 1
138    END DO
139  END DO
140
141  CALL borders_matrixL(dx, dy, Nppt, boolmat, borders, isborder, isbordery)
142  CALL paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)
143
144  Npoly = Npath
145
146  DO ip=1, Npath
147    IF (ALLOCATED(vertxs)) DEALLOCATE(vertxs)
148    ALLOCATE(vertxs(Nptpaths(ip),2))
149    msg = "Problems allocating matrix 'vertxs'"
150    CALL ErrMsg(msg, fname, ierr)
151
152    IF (ALLOCATED(isin)) DEALLOCATE(isin)
153    ALLOCATE(isin(Npts), STAT=ierr)
154    msg = "Problems allocating matrix 'isin'"
155    CALL ErrMsg(msg, fname, ierr)
156
157    isin = .FALSE.
158
159    IF (dbg) PRINT *, '  path:', ip, ' N pts:', Nptpaths(ip)
160
161    CALL path_properties(dx, dy, boolmat, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), xtrx, xtry,       &
162      meanpth, 'y', Nvertx, vertxs)
163
164    IF (dbg) THEN
165      PRINT *, '    properties  _______'
166      PRINT *, '    x-extremes:', xtrx
167      PRINT *, '    y-extremes:', xtry
168      PRINT *, '    center mean:', meanpth
169      PRINT *, '    y-vertexs:', Nvertx,' ________'
170      DO i=1, Nvertx
171        PRINT *,'      ',i,':',vertxs(i,:)
172      END DO
173    END IF
174 
175    CALL gridpoints_InsidePolygon(dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), Nvertx,&
176      xtrx, xtry, vertxs, Npts, points, isin)
177
178    ! Filling polygons
179    ipp = 1
180    DO i=1, dx
181      DO j=1, dy
182        IF (isin(ipp)) polys(i,j) = ip
183        ipp = ipp + 1
184      END DO
185    END DO
186
187    IF (dbg) THEN
188      PRINT *,'  boolmat isborder isbordery polygon (',xtrx(1),',',xtry(1),')x(',xtrx(2),',',xtry(2), &
189        ') _______'
190      DO i=xtrx(1), xtrx(2)
191        PRINT *,i,':',boolmat(i,xtry(1):xtry(2)), ' border ', isborder(i,xtry(1):xtry(2)),            &
192          ' isbordery ', isbordery(i,xtry(1):xtry(2)), ' polygon ', polys(i,xtry(1):xtry(2))
193      END DO
194    END IF
195
196  END DO
197
198  DEALLOCATE (borders) 
199  DEALLOCATE (Nptpaths)
200  DEALLOCATE (paths)
201  DEALLOCATE (vertxs)
202  DEALLOCATE (points)
203  DEALLOCATE (isin)
204
205  RETURN
206
207END SUBROUTINE polygons
208
209  SUBROUTINE path_properties(dx, dy, Lmat, Nptspth, pth, xxtrm, yxtrm, meanctr, axs, Nvrtx, vrtxs)
210! Subroutine to determine the properties of a path:
211!   extremes: minimum and maximum of the path along x,y axes
212!   meancenter: center from the mean of the coordinates of the paths locations
213!   vertexs: path point, without neighbours along a given axis
214
215  IMPLICIT NONE
216
217  INTEGER, INTENT(in)                                    :: dx, dy, Nptspth
218  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: Lmat
219  INTEGER, DIMENSION(Nptspth,2), INTENT(in)              :: pth
220  CHARACTER, INTENT(in)                                  :: axs
221  INTEGER, DIMENSION(2), INTENT(out)                     :: meanctr, xxtrm, yxtrm
222  INTEGER, INTENT(out)                                   :: Nvrtx
223  INTEGER, DIMENSION(Nptspth,2), INTENT(out)             :: vrtxs
224
225! Local
226  INTEGER                                                :: i, ip, jp
227  INTEGER                                                :: neig1, neig2
228
229!!!!!!! Variables
230! dx,dy: size of the space
231! Lmat: original matrix of logical values for the path
232! Nptspth: number of points of the path
233! pth: path coordinates (clockwise)
234! axs: axis of finding the vertex
235! [x/y]xtrm: minimum and maximum coordinates of the path
236! meanctr: center from the mean of the coordinates of the path
237! Nvrtx: Number of vertexs of the path
238! vrtxs: coordinates of the vertexs
239
240  fname = 'path_properties'
241
242  vrtxs = -1
243  Nvrtx = 0
244
245  xxtrm = (/ MINVAL(pth(:,1)), MAXVAL(pth(:,1)) /)
246  yxtrm = (/ MINVAL(pth(:,2)), MAXVAL(pth(:,2)) /)
247  meanctr = (/ SUM(pth(:,1))/Nptspth, SUM(pth(:,2))/Nptspth /)
248
249  IF (axs == 'x' .OR. axs == 'X') THEN
250    ! Looking vertexs along x-axis
251    DO i=1, Nptspth
252      ip = pth(i,1)
253      jp = pth(i,2)
254      neig1 = 0
255      neig2 = 0
256      ! W-point
257      IF (ip == 1) THEN
258        neig1 = -1
259      ELSE
260        IF (.NOT.Lmat(ip-1,jp)) neig1 = -1
261      END IF
262      ! E-point
263      IF (ip == dx) THEN
264        neig2 = -1
265      ELSE
266        IF (.NOT.Lmat(ip+1,jp)) neig2 = -1
267      END IF
268   
269      IF (neig1 == -1 .AND. neig2 == -1) THEN
270        Nvrtx = Nvrtx + 1
271        vrtxs(Nvrtx,:) = (/ip,jp/)
272      END IF
273    END DO
274  ELSE IF (axs == 'y' .OR. axs == 'Y') THEN
275    ! Looking vertexs along x-axis
276    DO i=1, Nptspth
277      ip = pth(i,1)
278      jp = pth(i,2)
279
280      neig1 = 0
281      neig2 = 0
282      ! S-point
283      IF (jp == 1) THEN
284        neig1 = -1
285      ELSE
286        IF (.NOT.Lmat(ip,jp-1)) neig1 = -1
287      END IF
288      ! N-point
289      IF (jp == dy) THEN
290        neig2 = -1
291      ELSE
292        IF (.NOT.Lmat(ip,jp+1)) neig2 = -1
293      END IF
294
295      IF (neig1 == -1 .AND. neig2 == -1) THEN
296        Nvrtx = Nvrtx + 1
297        vrtxs(Nvrtx,:) = (/ ip, jp /)
298      END IF
299    END DO
300  ELSE
301    msg = "Axis '" // axs // "' not available" // CHAR(10) // "  Available ones: 'x', 'X', 'y, 'Y'"
302    CALL ErrMsg(msg, fname, -1)
303  END IF
304
305  RETURN
306
307  END SUBROUTINE path_properties
308
309  SUBROUTINE gridpoints_InsidePolygon(dx, dy, isbrdr, Npath, path, Nvrtx, xpathxtrm, ypathxtrm,       &
310    vrtxs, Npts, pts, inside)
311! Subroutine to determine if a series of grid points are inside a polygon following ray casting algorithm
312! FROM: https://en.wikipedia.org/wiki/Point_in_polygon
313
314  IMPLICIT NONE
315
316  INTEGER, INTENT(in)                                    :: dx,dy,Npath,Nvrtx,Npts
317  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: isbrdr
318  INTEGER, DIMENSION(Npath,2), INTENT(in)                :: path
319  INTEGER, DIMENSION(2), INTENT(in)                      :: xpathxtrm, ypathxtrm
320  INTEGER, DIMENSION(Npath,2)                            :: vrtxs
321  INTEGER, DIMENSION(Npts,2), INTENT(in)                 :: pts
322  LOGICAL, DIMENSION(Npts), INTENT(out)                  :: inside
323
324! Local
325  INTEGER                                                :: i,j,ip,ix,iy
326  INTEGER                                                :: Nintersecs, isvertex, ispath
327  INTEGER                                                :: ierr
328  LOGICAL, DIMENSION(:,:), ALLOCATABLE                   :: halo_brdr
329  INTEGER                                                :: Nbrbrdr
330
331!!!!!!! Variables
332! dx,dy: space size
333! Npath: number of points of the path of the polygon
334! path: path of the polygon
335! isbrdr: boolean matrix of the space wqith .T. on polygon border
336! Nvrtx: number of vertexs of the path
337! [x/y]pathxtrm extremes of the path
338! vrtxs: vertexs of the path along y-axis
339! Npts: number of points
340! pts: points to look for
341! inside: vector wether point is inside or not (coincident to a border is inside)
342
343  fname = 'gridpoints_InsidePolygon'
344
345  ! Creation of a 1-grid point larger matrix to deal with points reaching the limits
346  IF (ALLOCATED(halo_brdr)) DEALLOCATE(halo_brdr)
347  ALLOCATE(halo_brdr(dx+2,dy+2), STAT=ierr)
348  msg = "Problems allocating matrix 'halo_brdr'"
349  CALL ErrMsg(msg, fname, ierr)
350  halo_brdr = .FALSE.
351
352  DO i=1,dx
353    halo_brdr(i+1,2:dy+1) = isbrdr(i,:)
354  END DO
355
356  inside = .FALSE.
357
358  DO ip=1,Npts
359    Nintersecs = 0
360    ix = pts(ip,1)
361    iy = pts(ip,2)
362    ! Point might be outside path range...
363    IF (ix >= xpathxtrm(1) .AND. ix <= xpathxtrm(2) .AND. iy >= ypathxtrm(1) .AND.                    &
364      iy <= ypathxtrm(2)) THEN
365
366      ! It is a border point?
367      ispath = index_list_coordsI(Npath, path, (/ix,iy/))
368      IF (isbrdr(ix,iy) .AND. (ispath /= -1)) THEN
369        inside(ip) = .TRUE.
370        CYCLE
371      END IF
372
373      ! Looking along y-axis
374      ! Accounting for consecutives borders
375      Nbrbrdr = 0
376      DO j=MAX(1,ypathxtrm(1)-1),iy-1
377        ! Only counting that borders that are not vertexs
378        ispath = index_list_coordsI(Npath, path, (/ix,j/))
379        isvertex = index_list_coordsI(Npath, vrtxs, (/ix,j/))
380
381        IF (halo_brdr(ix+1,j+1) .AND. (ispath /= -1) .AND. (isvertex == -1) ) Nintersecs = Nintersecs + 1
382        IF (halo_brdr(ix+1,j+1) .AND. (ispath /= -1) .AND. (halo_brdr(ix+1,j+1) .EQV. halo_brdr(ix+1,j+2))) THEN
383          Nbrbrdr = Nbrbrdr + 1
384        ELSE
385          ! Will remove that consecutive borders above 2
386          IF (Nbrbrdr /= 0) THEN
387            Nintersecs = Nintersecs - MAX(Nbrbrdr-1, 0)
388            Nbrbrdr = 0
389          END IF
390        END IF
391      END DO
392      IF (MOD(Nintersecs,2) /= 0) inside(ip) = .TRUE.
393    END IF
394
395  END DO
396
397  RETURN
398
399END SUBROUTINE gridpoints_InsidePolygon
400
401SUBROUTINE look_clockwise_borders(dx,dy,Nbrdrs,brdrs,gbrdr,isbrdr,ix,iy,dbg,xf,yf,iff)
402! Subroutine to look clock-wise for a next point within a collection of borders (limits of a region)
403
404  IMPLICIT NONE
405
406  INTEGER, INTENT(in)                                    :: dx, dy, Nbrdrs, ix, iy
407  INTEGER, DIMENSION(Nbrdrs,2), INTENT(in)               :: brdrs
408  LOGICAL, DIMENSION(Nbrdrs), INTENT(in)                 :: gbrdr
409  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: isbrdr
410  LOGICAL, INTENT(in)                                    :: dbg
411  INTEGER, INTENT(out)                                   :: xf, yf, iff
412
413! Local
414  INTEGER                                                :: isch
415  CHARACTER(len=2), DIMENSION(8)                         :: Lclock
416  INTEGER, DIMENSION(8,2)                                :: spt
417  INTEGER                                                :: iif, jjf
418
419!!!!!!! Variables
420! dx, dy: 2D shape ot the space
421! Nbrdrs: number of brdrs found in this 2D space
422! brdrs: list of coordinates of the borders
423! gbrdr: accounts for the use if the given border point
424! isbrdr: accounts for the matrix of the point is a border or not
425! ix,iy: coordinates of the point to start to find for
426! xf,yf: coordinates of the found point
427! iff: position of the border found within the list of borders
428
429  fname = 'look_clockwise_borders'
430
431  ! Looking clock-wise assuming that one starts from the westernmost point
432
433  ! Label of the search
434  lclock = (/ 'W ', 'NW', 'N ', 'NE', 'E ', 'SE', 'S ', 'SW' /)
435  ! Transformation to apply
436  !spt = (/ (/-1,0/), (/-1,1/), (/0,1/), (/1,1/), (/1,0/), (/1,-1/), (/0,-1/), (/-1,-1/) /)
437  spt(:,1) = (/ -1, -1, 0, 1, 1, 1, 0, -1 /)
438  spt(:,2) = (/ 0, 1, 1, 1, 0, -1, -1, -1 /)
439
440  xf = -1
441  yf = -1
442  DO isch=1, 8
443    ! clock-wise search
444    IF (spt(isch,1) >= 0) THEN
445      iif = MIN(dx,ix+spt(isch,1))
446    ELSE
447      iif = MAX(1,ix+spt(isch,1))
448    END IF
449    IF (spt(isch,2) >= 0) THEN
450      jjf = MIN(dy,iy+spt(isch,2))
451    ELSE
452      jjf = MAX(1,iy+spt(isch,2))
453    END IF
454    iff = index_list_coordsI(Nbrdrs, brdrs,(/iif,jjf/))
455    IF (iff > 0) THEN
456      IF (dbg) PRINT *,'    ' // lclock(isch) // '-point:', iif,jjf, ':', iff, 'is',isbrdr(iif,jjf),  &
457        'got',gbrdr(iff)
458      IF (isbrdr(iif,jjf) .AND. .NOT.gbrdr(iff)) THEN
459        xf = iif
460        yf = jjf
461        EXIT
462      END IF
463    END IF
464  END DO
465
466  RETURN
467
468END SUBROUTINE look_clockwise_borders
469
470SUBROUTINE borders_matrixL(dx,dy,dxy,Lmat,brdrs,isbrdr,isbrdry)
471! Subroutine to provide the borders of a logical array (interested in .TRUE.)
472
473  IMPLICIT NONE
474
475  INTEGER, INTENT(in)                                    :: dx,dy,dxy
476  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: Lmat
477  INTEGER, DIMENSION(dxy,2), INTENT(out)                 :: brdrs
478  LOGICAL, DIMENSION(dx,dy), INTENT(out)                 :: isbrdr, isbrdry
479
480! Local
481  INTEGER                                                :: i,j,ib
482
483!!!!!!! Variables
484! dx,dy: size of the space
485! dxy: maximum number of border points
486! Lmat: Matrix to look for the borders
487! brdrs: list of coordinates of the borders
488! isbrdr: matrix with .T./.F. wether the given matrix point is a border or not
489! isbrdry: matrix with .T./.F. wether the given matrix point is a border or not only along y-axis
490
491  fname = 'borders_matrixL'
492
493  isbrdr = .FALSE.
494  brdrs = -1
495  ib = 1
496
497  ! Starting with the borders. If a given point is TRUE it is a path-vertex
498  ! Along y-axis
499  DO i=1, dx
500    IF (Lmat(i,1) .AND. .NOT.isbrdr(i,1)) THEN
501      brdrs(ib,1) = i
502      brdrs(ib,2) = 1
503      isbrdr(i,1) = .TRUE.
504      ib=ib+1
505    END IF
506    IF (Lmat(i,dy) .AND. .NOT.isbrdr(i,dy)) THEN
507      brdrs(ib,1) = i
508      brdrs(ib,2) = dy
509      isbrdr(i,dy) = .TRUE.
510      ib=ib+1
511    END IF
512  END DO
513  ! Along x-axis
514  DO j=1, dy
515    IF (Lmat(1,j) .AND. .NOT.isbrdr(1,j)) THEN
516      brdrs(ib,1) = 1
517      brdrs(ib,2) = j
518      isbrdr(1,j) = .TRUE.
519      ib=ib+1
520     END IF
521    IF (Lmat(dx,j) .AND. .NOT.isbrdr(dx,j)) THEN
522      brdrs(ib,1) = dx
523      brdrs(ib,2) = j
524      isbrdr(dx,j) = .TRUE.
525      ib=ib+1
526    END IF
527  END DO
528
529  isbrdry = isbrdr
530
531  ! Border as that when looking on x-axis points with Lmat(i) /= Lmat(i+1)
532  DO i=1, dx-1
533    DO j=1, dy-1
534      IF ( Lmat(i,j) .NEQV. Lmat(i+1,j) ) THEN
535        IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN
536          brdrs(ib,1) = i
537          brdrs(ib,2) = j
538          isbrdr(i,j) = .TRUE.
539          ib=ib+1
540        ELSE IF (Lmat(i+1,j) .AND. .NOT.isbrdr(i+1,j)) THEN
541          brdrs(ib,1) = i+1
542          brdrs(ib,2) = j
543          isbrdr(i+1,j) = .TRUE.
544          ib=ib+1
545        END IF
546      END IF
547      ! y-axis
548      IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN
549        IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN
550          brdrs(ib,1) = i
551          brdrs(ib,2) = j
552          isbrdr(i,j) = .TRUE.
553          isbrdry(i,j) = .TRUE.
554          ib=ib+1
555        ELSE IF (Lmat(i,j+1) .AND. .NOT.isbrdr(i,j+1)) THEN
556          brdrs(ib,1) = i
557          brdrs(ib,2) = j+1
558          isbrdr(i,j+1) = .TRUE.
559          isbrdry(i,j+1) = .TRUE.
560          ib=ib+1
561        END IF
562      END IF
563    END DO       
564  END DO
565
566  DO i=1, dx-1
567    DO j=1, dy-1
568      ! y-axis
569      IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN
570        IF (Lmat(i,j)) THEN
571          isbrdry(i,j) = .TRUE.
572        ELSE IF (Lmat(i,j+1)) THEN
573          isbrdry(i,j+1) = .TRUE.
574        END IF
575      END IF
576    END DO       
577  END DO
578  ! only y-axis adding bands of 2 grid points
579  DO i=1, dx-1
580    DO j=2, dy-2
581      IF ( (Lmat(i,j) .EQV. Lmat(i,j+1)) .AND. (Lmat(i,j).NEQV.Lmat(i,j-1)) .AND. (Lmat(i,j).NEQV.Lmat(i,j+2)) ) THEN
582        IF (Lmat(i,j)) THEN
583          isbrdry(i,j) = .TRUE.
584          isbrdry(i,j+1) = .TRUE.
585        END IF
586      END IF
587    END DO       
588  END DO
589
590  RETURN
591
592END SUBROUTINE borders_matrixL
593
594SUBROUTINE paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths)
595! Subroutine to search the paths of a border field.
596
597  IMPLICIT NONE
598
599  INTEGER, INTENT(in)                                    :: dx, dy, Nppt
600  LOGICAL, INTENT(in)                                    :: dbg
601  LOGICAL, DIMENSION(dx,dy), INTENT(in)                  :: isborder
602  INTEGER, DIMENSION(Nppt,2), INTENT(in)                 :: borders
603  INTEGER, DIMENSION(Nppt,Nppt,2), INTENT(out)           :: paths
604  INTEGER, INTENT(out)                                   :: Npath
605  INTEGER, DIMENSION(Nppt), INTENT(out)                  :: Nptpaths
606
607! Local
608  INTEGER                                                :: i,j,k,ib
609  INTEGER                                                :: ierr
610  INTEGER                                                :: Nbrdr
611  LOGICAL, DIMENSION(:), ALLOCATABLE                     :: gotbrdr, emptygotbrdr
612  INTEGER                                                :: iipth, ipath, ip, Nptspath
613  INTEGER                                                :: iib, jjb, iip, ijp, iif, jjf, iff
614  LOGICAL                                                :: found, finishedstep
615
616!!!!!!! Variables
617! dx,dy: spatial dimensions of the space
618! Nppt: possible number of paths and points that the paths can have
619! isborder: boolean matrix which provide the borders of the polygon
620! borders: coordinates of the borders of the polygon
621! paths: coordinates of each found path
622! Npath: number of paths found
623! Nptpaths: number of points per path
624
625  fname = 'paths_border'
626
627  ! Sarting matrix
628  paths = -1
629  Npath = 0
630  Nptspath = 0
631  Nptpaths = -1
632
633  ib=1
634  finishedstep = .FALSE.
635
636  ! Number of border points
637  DO ib=1, Nppt
638    IF (borders(ib,1) == -1 ) EXIT
639  END DO
640  Nbrdr = ib-1
641   
642  IF (dbg) THEN
643    PRINT *,'    borders _______'
644    DO i=1,Nbrdr
645      PRINT *,'    ',i,':',borders(i,:)
646    END DO
647  END IF
648
649  ! Matrix which keeps track if a border point has been located
650  IF (ALLOCATED(gotbrdr)) DEALLOCATE(gotbrdr)
651  ALLOCATE(gotbrdr(Nbrdr), STAT=ierr)
652  msg = "Problems allocating matrix 'gotbrdr'"
653  CALL ErrMsg(msg, fname, ierr)
654  IF (ALLOCATED(emptygotbrdr)) DEALLOCATE(emptygotbrdr)
655  ALLOCATE(emptygotbrdr(Nbrdr), STAT=ierr)
656  msg = "Problems allocating matrix 'emptygotbrdr'"
657  CALL ErrMsg(msg, fname, ierr)
658
659  gotbrdr = .FALSE.
660  emptygotbrdr = .FALSE.
661
662  ! Starting the fun...
663   
664  ! Looking along the lines and when a border is found, starting from there in a clock-wise way
665  iipth = 1
666  ipath = 1   
667  DO ib=1,Nbrdr
668    iib = borders(iipth,1)
669    jjb = borders(iipth,2)
670    ! Starting new path
671    newpath: IF (.NOT.gotbrdr(iipth)) THEN
672      ip = 1
673      Nptspath = 1
674      paths(ipath,ip,:) = borders(iipth,:)
675      gotbrdr(iipth) = .TRUE.
676      ! Looking for following clock-wise search
677      ! Not looking for W, because search starts from the W
678      iip = iib
679      ijp = jjb
680      DO k=1,Nbrdr
681        IF (dbg) PRINT *,ipath,'iip jip:', iip, ijp
682        found = .FALSE.
683        CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE.,iif,jjf,iff)
684        IF (iif /= -1) THEN
685          ip=ip+1
686          paths(ipath,ip,:) = (/ iif,jjf /)
687          found = .TRUE.
688          gotbrdr(iff) = .TRUE.
689          iip = iif
690          ijp = jjf
691          Nptspath = Nptspath + 1         
692        END IF
693
694        IF (dbg) THEN
695          PRINT *,iib,jjb,'    end of this round path:', ipath, '_____', gotbrdr
696          DO i=1, Nptspath
697            PRINT *,'      ',i,':',paths(ipath,i,:)
698          END DO
699        END IF
700        ! If it is not found a next point, might be because it is a non-polygon related value
701        IF (.NOT.found) THEN
702          IF (dbg) PRINT *,'NOT FOUND !!!', gotbrdr
703          ! Are still there available borders? 
704          IF (ALL(gotbrdr) .EQV. .TRUE.) THEN
705            finishedstep = .TRUE.
706            Npath = ipath
707            Nptpaths(ipath) = Nptspath
708            EXIT
709          ELSE
710            Nptpaths(ipath) = Nptspath
711            ! Let's have a look if the previous points in the path have already some 'non-located' neighbourgs
712            DO i=Nptspath,1,-1
713              iip = paths(ipath,i,1)
714              ijp = paths(ipath,i,2)
715              CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE., iif,   &
716                jjf,iff)
717              IF (iif /= -1 .AND. iff /= -1) THEN
718                IF (dbg) PRINT *,'    re-take path from point:', iif,',',jjf,' n-path:', iff
719                found = .TRUE.
720                iipth = index_list_coordsI(Nppt, borders, (/iip,ijp/))
721                EXIT
722              END IF
723            END DO
724            IF (.NOT.found) THEN
725              ! Looking for the next available border point for the new path
726              DO i=1,Nbrdr
727                IF (.NOT.gotbrdr(i)) THEN
728                  iipth = i
729                  EXIT
730                END IF
731              END DO
732              IF (dbg) PRINT *,'  Looking for next path starting at:', iipth, ' point:',              &
733                borders(iipth,:)
734              ipath=ipath+1
735              EXIT
736            END IF
737          END IF
738        ELSE
739          IF (dbg) PRINT *,'  looking for next point...'
740        END IF
741        IF (finishedstep) EXIT
742      END DO
743    END IF newpath
744  END DO
745  Npath = ipath
746  Nptpaths(ipath) = Nptspath
747   
748  DEALLOCATE (gotbrdr)
749  DEALLOCATE (emptygotbrdr)
750
751  RETURN
752
753END SUBROUTINE paths_border
754
755SUBROUTINE rand_sample(Nvals, Nsample, sample)
756! Subroutine to randomly sample a range of indices
757
758  IMPLICIT NONE
759
760  INTEGER, INTENT(in)                                    :: Nvals, Nsample
761  INTEGER, DIMENSION(Nsample), INTENT(out)               :: sample
762
763! Local
764  INTEGER                                                :: i, ind, jmax
765  REAL, DIMENSION(Nsample)                               :: randv
766  CHARACTER(len=50)                                      :: fname
767  LOGICAL                                                :: found
768  LOGICAL, DIMENSION(Nvals)                              :: issampled
769  CHARACTER(len=256)                                     :: msg
770  CHARACTER(len=10)                                      :: IS1, IS2
771
772!!!!!!! Variables
773! Nvals: number of values
774! Nsamples: number of samples
775! sample: samnple
776  fname = 'rand_sample'
777
778  IF (Nsample > Nvals) THEN
779    WRITE(IS1,'(I10)')Nvals
780    WRITE(IS2,'(I10)')Nsample
781    msg = 'Sampling of ' // TRIM(IS1) // ' is too big for ' // TRIM(IS1) // 'values'
782    CALL ErrMsg(msg, fname, -1)
783  END IF
784
785  ! Generation of random numbers always the same series during the whole program!
786  CALL RANDOM_NUMBER(randv)
787
788  ! Making sure that we do not repeat any value
789  issampled = .FALSE.
790
791  DO i=1, Nsample
792    ! Generation of the index from the random numbers
793    ind = MAX(INT(randv(i)*Nvals), 1)
794
795    IF (.NOT.issampled(ind)) THEN
796      sample(i) = ind
797      issampled(ind) = .TRUE.
798    ELSE
799      ! Looking around the given index
800      !PRINT *,' Index :', ind, ' already sampled!', issampled(ind)
801      found = .FALSE.
802      DO jmax=1, Nvals
803        ind = MIN(ind+jmax, Nvals)
804        IF (.NOT.issampled(ind)) THEN
805          sample(i) = ind
806          issampled(ind) = .TRUE.
807          found = .TRUE.
808          EXIT
809        END IF
810        ind = MAX(1, ind-jmax)
811        IF (.NOT.issampled(ind)) THEN
812          sample(i) = ind
813          issampled(ind) = .TRUE.
814          found = .TRUE.
815          EXIT
816        END IF
817      END DO
818      IF (.NOT.found) THEN
819        msg = 'sampling could not be finished due to absence of available value!!'
820        CALL ErrMsg(msg, fname, -1)
821      END IF
822    END IF
823
824  END DO
825
826  RETURN
827
828END SUBROUTINE rand_sample
829
830SUBROUTINE PrintQuantilesR_K(Nvals, vals, Nquants, qtvs, bspc)
831! Subroutine to print the quantiles of values REAL(r_k)
832
833  IMPLICIT NONE
834
835  INTEGER, INTENT(in)                                    :: Nvals, Nquants
836  REAL(r_k), DIMENSION(Nvals), INTENT(in)                :: vals
837  REAL(r_k), DIMENSION(Nquants), INTENT(in)              :: qtvs
838  CHARACTER(len=1000), OPTIONAL                          :: bspc
839
840! Local
841  INTEGER                                                :: iq
842  LOGICAL, DIMENSION(Nvals)                              :: search1, search2, search
843  CHARACTER(len=6)                                       :: RS1
844  CHARACTER(len=50)                                      :: fname
845  CHARACTER(len=1000)                                    :: bspcS
846
847!!!!!!! Variables
848! vals: series of values
849! qtvs: values of the quantiles
850! bspc: base quantity of spaces
851
852  fname = 'PrintQuantilesR_K'
853
854  IF (PRESENT(bspc)) THEN
855    bspcS = bspc
856  ELSE
857    bspcS = '      '
858  END IF
859
860  DO iq=1, Nquants-1
861
862    WHERE (vals >= qtvs(iq))
863      search1 = .TRUE.
864    ELSEWHERE
865      search1 = .FALSE.
866    END WHERE
867
868    WHERE (vals < qtvs(iq+1))
869      search2 = .TRUE.
870    ELSEWHERE
871      search2 = .FALSE.
872    END WHERE
873
874    WHERE (search1 .AND. search2)
875      search = .TRUE.
876    ELSEWHERE
877      search = .FALSE.
878    END WHERE
879
880    WRITE(RS1, '(F6.2)')(iq)*100./(Nquants-1)
881    PRINT *, TRIM(bspcS) // '[',iq,']', TRIM(RS1) // ' %:', qtvs(iq), 'N:', COUNT(search)
882
883  END DO
884
885  RETURN
886
887END SUBROUTINE PrintQuantilesR_K
888
889   INTEGER FUNCTION FindMinimumR_K(x, dsize, Startv, Endv)
890! Function returns the location of the minimum in the section between Start and End.
891
892      IMPLICIT NONE
893
894      INTEGER, INTENT(in)                                :: dsize
895      REAL(r_k), DIMENSION(dsize), INTENT(in)            :: x
896      INTEGER, INTENT(in)                                :: Startv, Endv
897
898! Local
899      REAL(r_k)                                          :: Minimum
900      INTEGER                                            :: Location
901      INTEGER                                            :: i
902
903      Minimum  = x(Startv)                               ! assume the first is the min
904      Location = Startv                                  ! record its position
905      DO i = Startv+1, Endv                              ! start with next elements
906         IF (x(i) < Minimum) THEN                        !   if x(i) less than the min?
907            Minimum  = x(i)                              !      Yes, a new minimum found
908            Location = i                                 !      record its position
909         END IF
910      END DO
911
912      FindMinimumR_K = Location                          ! return the position
913
914   END FUNCTION  FindMinimumR_K
915
916   SUBROUTINE SwapR_K(a, b)
917! Subroutine swaps the values of its two formal arguments.
918
919      IMPLICIT  NONE
920
921      REAL(r_k), INTENT(INOUT)                           :: a, b
922! Local
923      REAL(r_k)                                          :: Temp
924
925      Temp = a
926      a    = b
927      b    = Temp
928
929   END SUBROUTINE  SwapR_K
930
931   SUBROUTINE  SortR_K(x, Nx)
932! Subroutine receives an array x() r_K and sorts it into ascending order.
933
934      IMPLICIT NONE
935
936      INTEGER, INTENT(IN)                                :: Nx
937      REAL(r_k), DIMENSION(Nx), INTENT(INOUT)            :: x
938
939! Local
940      INTEGER                                            :: i
941      INTEGER                                            :: Location
942
943      DO i = 1, Nx-1                                     ! except for the last
944         Location = FindMinimumR_K(x, Nx-i+1, i, Nx)     ! find min from this to last
945         CALL  SwapR_K(x(i), x(Location))                ! swap this and the minimum
946      END DO
947
948   END SUBROUTINE  SortR_K
949
950SUBROUTINE quantilesR_K(Nvals, vals, Nquants, quants)
951! Subroutine to provide the quantiles of a given set of values of type real 'r_k'
952
953  IMPLICIT NONE
954
955  INTEGER, INTENT(in)                                    :: Nvals, Nquants
956  REAL(r_k), DIMENSION(Nvals), INTENT(in)                :: vals
957  REAL(r_k), DIMENSION(Nquants), INTENT(out)             :: quants
958
959! Local
960  INTEGER                                                :: i
961  REAL(r_k)                                              :: minv, maxv
962  REAL(r_k), DIMENSION(Nvals)                            :: sortedvals
963
964!!!!!!! Variables
965! Nvals: number of values
966! Rk: kind of real of the values
967! vals: values
968! Nquants: number of quants
969! quants: values at which the quantile start
970
971  minv = MINVAL(vals)
972  maxv = MAXVAL(vals)
973
974  sortedvals = vals
975  ! Using from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90
976  CALL SortR_K(sortedvals, Nvals)
977
978  quants(1) = minv
979  DO i=2, Nquants
980    quants(i) = sortedvals(INT((i-1)*Nvals/Nquants))
981  END DO
982
983END SUBROUTINE quantilesR_K
984
985END MODULE module_scientific
Note: See TracBrowser for help on using the repository browser.