MODULE module_scientific ! Module of the scientific function/subroutines !!!!!!! Functions & subroutines ! all_polygons_properties: Subroutine to determine the properties of all polygons in a 2D field: ! borders_matrixL: Subroutine to provide the borders of a logical array (interested in .TRUE.) ! coincidence_all_polys: Subtourine to determine which is the coincident polygon when a boolean polygon is provided ! to a map of integer polygons ! coincidence_all_polys_area: Subtourine to determine which is the coincident polygon when a boolean polygon is provided ! to a map of integer polygons filtrered by a minimal area ! coincidence_poly: Subtourine to determine which is the coincident polygon when a boolean polygon is provided ! to a map of integer polygons ! coincidence_poly_area: Subtourine to determine which is the coincident polygon when a boolean polygon is provided ! to a map of integer polygons filtrered by a minimal area ! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole ! coincident_polygon: Subroutine to provide the intersection polygon between two polygons ! distanceRK: Function to provide the distance between two points ! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End. ! fill3DI_2Dvec: Subroutine to fill a 3D integer matrix from a series of indices from a 2D matrix ! fill3DR_2Dvec: Subroutine to fill a 3D float matrix from a series of indices from a 2D matrix ! grid_within_polygon: Subroutine to determine which grid cells from a matrix lay inside a polygon ! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon ! following ray casting algorithm ! Heron_area_triangle: Subroutine to compute area of a triangle using Heron's formula ! intersectfaces: Subroutine to provide if two faces of two polygons intersect ! intersection_2Dlines: Subroutine to provide the intersection point between two lines on the plane using Cramer's method ! join_polygon: Subroutine to join two polygons ! look_clockwise_borders: Subroutine to look clock-wise for a next point within a collection of borders ! (limits of a region) ! NcountR: Subroutine to count real values ! paths_border: Subroutine to search the paths of a border field. ! path_properties: Subroutine to determine the properties of a path ! percentiles_R_K*D: Subroutine to compute the percentiles of a *D R_K array along given set of axis ! point_in_face: Function to determine if a given point is on a face of a polygon ! point_inside: Function to determine if a given point is inside a polygon providing its sorted vertices ! poly_has_point: Function to determine if a polygon has already a given point as one of its vertex ! polygon_properties: Subroutine to determine the properties of a polygon (as .TRUE. matrix) ! polygons: Subroutine to search the polygons of a border field. FORTRAN based. 1st = 1! ! polygons_t: Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1! ! poly_overlap_tracks: Subroutine to determine tracks of a series of consecutive 2D field with polygons using ! maximum overlaping/coincidence ! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k) ! poly_overlap_tracks_area: Subroutine to determine tracks of a series of consecutive 2D field with polygons using ! maximum overlaping/coincidence filtrered by a minimal area ! poly_overlap_tracks_area_ascii: Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum ! overlaping/coincidence filtrered by a minimal area writting theoutput on an ASCII file (memory limitations) ! quantilesR_K: Subroutine to provide the quantiles of a given set of values of type real 'r_k' ! rand_sample: Subroutine to randomly sample a range of indices ! read_finaltrack_ascii: Subroutine to read the final trajectories from an ASCII file ! read_overlap_single_track_ascii: Subroutine to read the values for a given trajectory ! read_overlap_polys_ascii: Subroutine to read from an ASCII file the associated polygons at a given time-step ! read_overlap_tracks_ascii: Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step ! reconstruct_matrix: Subroutine to reconstruct a 2D matrix from a pair of syncronized vectors with the positions on x ! and y coordinates ! runmean_F1D: Subroutine fo computing the running mean of a given set of float 1D values ! shoelace_area_polygon: Computing the area of a polygon using sholace formula ! sort_polygon: Subroutine to sort a polygon using its center as average of the coordinates and remove duplicates ! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order. ! StatsR_K: Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a ! series of r_k numbers ! SwapR_K*: Subroutine swaps the values of its two formal arguments. ! write_finaltrack_ascii: Subroutine to read the final trajectories into an ASCII file ! write_overlap_polys_ascii: Subroutine to write to an ASCII file the associated polygons at a given time-step ! write_overlap_tracks_ascii: Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step !!! *Functions/Subroutines to sort values adpated. The method used is usually referred to as "selection" method. ! from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90 USE module_definitions USE module_generic CONTAINS SUBROUTINE fill3DI_2Dvec(matind, inmat, id1, id2, od1, od2, od3, outmat) ! Subroutine to fill a 3D integer matrix from a series of indices from a 2D matrix IMPLICIT NONE INTEGER, INTENT(in) :: id1, id2, od1, od2, od3 INTEGER, DIMENSION(id1,id2), INTENT(in) :: inmat INTEGER, DIMENSION(od1,od2), INTENT(in) :: matind INTEGER, DIMENSION(od1,od2,od3), INTENT(out) :: outmat ! Local INTEGER :: i, j !!!!!!! Variables ! matind: equivalence on the 2D space of the location in inmat ! inmat: values of the input matrix (1Dvec, time) ! id1/2: dimensions of the input matrix ! od1/3: dimensions of the output matrix ! outmat: output matrix ! NOTE: id2 == od3 fname = 'fill3DR_2Dvec' DO i=1, od1 DO j=1, od2 IF (matind(i,j) /= -1) THEN outmat(i,j,:) = inmat(matind(i,j),:) ELSE outmat(i,j,:) = fillvalI END IF END DO END DO END SUBROUTINE fill3DI_2Dvec SUBROUTINE fill3DR_2Dvec(matind, inmat, id1, id2, od1, od2, od3, outmat) ! Subroutine to fill a 3D float matrix from a series of indices from a 2D matrix IMPLICIT NONE INTEGER, INTENT(in) :: id1, id2, od1, od2, od3 REAL(r_k), DIMENSION(id1,id2), INTENT(in) :: inmat INTEGER, DIMENSION(od1,od2), INTENT(in) :: matind REAL(r_k), DIMENSION(od1,od2,od3), INTENT(out) :: outmat ! Local INTEGER :: i, j !!!!!!! Variables ! matind: equivalence on the 2D space of the location in inmat ! inmat: values of the input matrix (1Dvec, time) ! id1/2: dimensions of the input matrix ! od1/3: dimensions of the output matrix ! outmat: output matrix ! NOTE: id2 == od3 fname = 'fill3DR_2Dvec' DO i=1, od1 DO j=1, od2 IF (matind(i,j) /= -1) THEN outmat(i,j,:) = inmat(matind(i,j),:) ELSE outmat(i,j,:) = fillval64 END IF END DO END DO END SUBROUTINE fill3DR_2Dvec SUBROUTINE reconstruct_matrix(vectorXpos, vectorYpos, dvec, Xmin, Xmax, Ymin, Ymax, dmatx, dmaty, & matProj, maxdiff, matind, matX, matY, matdiff) ! Subroutine to reconstruct a 2D matrix from a pair of syncronized vectors with the positions on x ! and y coordinates IMPLICIT NONE INTEGER, INTENT(in) :: dvec, dmatx, dmaty REAL(r_k), DIMENSION(dvec), INTENT(in) :: vectorXpos, vectorYpos REAL(r_k), INTENT(in) :: Xmin, Xmax, Ymin, Ymax, maxdiff CHARACTER(len=50), INTENT(in) :: matPRoj INTEGER, DIMENSION(dmatx, dmaty), INTENT(out) :: matind REAL(r_k), DIMENSION(dmatx, dmaty), INTENT(out) :: matX, matY, matdiff ! Local INTEGER :: i,j,iv, idiff, jdiff REAL(r_k) :: ddx, ddy, Xpos, Ypos, mindiff REAL(r_k), DIMENSION(dmatx,dmaty) :: diff REAL(r_k) :: nXpos, xXpos, nYpos, xYpos INTEGER, DIMENSION(2) :: mindiffloc CHARACTER(LEN=50) :: RSa, RSb !!!!!!! Variables ! vectorXpos, vectorYpos: syncronized vectors with the x,y coordinates from which the matrix will be reconstructed ! dvec: quantitiy of coordinates to use ! Xmin,Xmax,Ymin,Ymax: Location range of the matrix to be reconstructed ! maxdiff: Authorized maximum distance between matrix location and vector location ! dmatx, dmaty: Size in grid points of the matrix to be reconstructed ! matProj: Assumed projection of the values of the matrix ! 'latlon': regular lat/lon projection ! matind: matrix with the correspondant indiced from the vector of positions ! matX, matY: matrix with the coordinates of the elements of the matrix ! matdiff: matrix with the differences at each grid point fname = 'reconstruct_matrix' matind = -oneRK matdiff = -oneRK nXpos = MINVAL(vectorXpos) xXpos = MAXVAL(vectorXpos) nYpos = MINVAL(vectorYpos) xYpos = MAXVAL(vectorYpos) PRINT *, 'Limits of the positions to localize in a matrix: longitudes:', nXpos, & ' ,',xXpos, ' latitudes:', nYpos, ' ,', xYpos Projection: SELECT CASE(TRIM(matProj)) CASE('latlon') ! Resolution along matrix axes ddx = (Xmax - Xmin)/(dmatx-1) ddy = (Ymax - Ymin)/(dmaty-1) ! Filling matrix positions DO i=1,dmatx Xpos = Xmin + ddx*(i-1) DO j=1,dmaty Ypos = Ymin + ddy*(j-1) matX(i,j) = Xpos matY(i,j) = Ypos END DO END DO CASE DEFAULT msg = "Projection of matrix '" // TRIM(matProj) // "' not ready " // CHAR(10) // & " Available ones: 'latlon'" CALL ErrMsg(msg, fname, -1) END SELECT Projection ! Let's do it properly... DO iv=1,dvec diff = SQRT((matX - vectorXpos(iv))**2 + (matY - vectorYpos(iv))**2) mindiffloc = MINLOC(diff) mindiff = MINVAL(diff) IF (mindiff > maxdiff) THEN PRINT *,' vectorXpos:', vectorXpos(iv), 'vectorYpos:', vectorYpos(iv) !PRINT *,' Xpos:', Xpos, 'Ypos:', Ypos WRITE(RSa, '(f20.10)')mindiff WRITE(RSb, '(f20.10)')maxdiff msg = 'Difference: ' // TRIM(RSa) // ' overpasses the maximum authorized distance: ' // & TRIM(RSb) CALL ErrMsg(msg, fname, -1) ELSE i = mindiffloc(1) j = mindiffloc(2) matind(i,j) = iv matdiff(i,j) = mindiff END IF END DO RETURN END SUBROUTINE reconstruct_matrix SUBROUTINE read_finaltrack_ascii(funit, dt, itrack, ftrack) ! Subroutine to read the final trajectories from an ASCII file IMPLICIT NONE INTEGER, INTENT(in) :: funit, dt, itrack REAL(r_k), DIMENSION(5,dt), INTENT(out) :: ftrack ! Local INTEGER :: i, j, it LOGICAL :: found !!!!!!! Variables ! funit: unit where to write the trajectory ! dt: number of time-steps ! itrack: trajectory to read the values ! ftrack: values of the trajectory fname = 'read_finaltrack_ascii' ftrack = 0. REWIND(funit) it = 1 DO WHILE (.NOT.found) READ(funit,10)ftrack(1,1), Str1, ((ftrack(i,j),Str1,i=2,5),Str1,j=1,dt) IF (INT(ftrack(1,1)) == itrack) THEN ftrack(1,2:dt) = ftrack(1,1) found = .TRUE. END IF ! Just in case IF (it >= dt) found = .TRUE. END DO RETURN 10 FORMAT(I10000000,1x,A1,1x,10000000(4(F20.10,A1),A1)) END SUBROUTINE read_finaltrack_ascii SUBROUTINE write_finaltrack_ascii(funit, dt, ftrack) ! Subroutine to write the final trajectories into an ASCII file IMPLICIT NONE INTEGER, INTENT(in) :: funit, dt REAL(r_k), DIMENSION(5,dt), INTENT(in) :: ftrack ! Local INTEGER :: i, j !!!!!!! Variables ! funit: unit where to write the trajectory ! dt: number of time-steps ! ftrack: values of the trajectory fname = 'write_finaltrack_ascii' WRITE(funit,10)INT(ftrack(1,1)), ';', ((ftrack(i,j), ',', i=2,5), ':', j=1,dt) RETURN 10 FORMAT(I10,1x,A1,1x,10000000(4(F20.10,A1),A1)) END SUBROUTINE write_finaltrack_ascii SUBROUTINE read_overlap_single_track_ascii(funit, dt, Nxp, Nxtr, itrack, strack) ! Subroutine to read the values for a given trajectory IMPLICIT NONE INTEGER, INTENT(in) :: funit, dt, Nxp, Nxtr, itrack REAL(r_k), DIMENSION(5,Nxp,dt), INTENT(out) :: strack ! Local INTEGER :: i,j,k,l INTEGER :: read_it, itt, it, Ntrcks INTEGER, DIMENSION(Nxp) :: Npindep LOGICAL :: looking REAL(r_k), DIMENSION(5,Nxp,Nxtr) :: trcks !!!!!!! Variables ! funit: unit from where retrieve the values of the trajectory ! dt: time-dimension ! Nxp: maximum allowed number of polygons per time-step ! Nxp: maximum allowed number of trajectories ! itrack: trajectory Id to look for ! strack: Values for the given trajectory fname = 'read_overlap_single_track_ascii' strack = 0. REWIND(funit) looking = .TRUE. itt = 0 it = 1 DO WHILE (looking) READ(funit,5,END=100)Str10, read_it READ(funit,*)Ntrcks DO i=1, Ntrcks READ(funit,10)l, Str1, Npindep(i), Str1, ((trcks(k,j,i),Str1,k=1,5),Str1,j=1,Npindep(i)) END DO ! There is the desired trajectory at this time-step? IF (ANY(INT(trcks(1,1,:)) == itrack)) THEN itt = itt + 1 DO i=1, Ntrcks IF (INT(trcks(1,1,i)) == itrack) THEN DO j=1, Npindep(i) strack(:,j,it) = trcks(:,j,i) END DO END IF END DO ELSE ! It trajectory has already been initialized this is the end IF (itt > 0) looking = .FALSE. END IF ! Just in case... ;) IF (read_it >= dt) looking = .FALSE. it = it + 1 IF (it > dt) looking = .FALSE. END DO 100 CONTINUE RETURN 5 FORMAT(A10,1x,I4) 10 FORMAT(I4,1x,A1,I4,1x,A1,1x,1000000(5(F20.10,A1),A1)) END SUBROUTINE read_overlap_single_track_ascii SUBROUTINE read_overlap_tracks_ascii(funit, tstep, Nxp, Ntrcks, trcks) ! Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step IMPLICIT NONE INTEGER, INTENT(in) :: funit, tstep, Nxp INTEGER, INTENT(out) :: Ntrcks REAL(r_k), DIMENSION(5,Nxp,Nxp), INTENT(out) :: trcks ! Local INTEGER :: i, j, k, l, Npindep INTEGER :: read_it !!!!!!! Variables ! funit: unit where to write the file ! tstep: time-step to write the trajectories ! Nxp: maximum number of polygons per time-step ! Nrtcks: Number of trajectories of the given time-step ! trcks: trajectories fname = 'read_overlap_tracks_ascii' Ntrcks = 0 trcks = 0 READ(funit,5)Str10, read_it IF (read_it /= tstep) THEN WRITE(numSa,'(I4)')read_it WRITE(numSb,'(I4)')tstep msg = 'File time-step;' // TRIM(numSa) // ' does not coincide with the one from program:' // & TRIM(numSb) END IF READ(funit,*)Ntrcks DO i=1, Ntrcks READ(funit,10)l, Str1, Npindep, Str1, ((trcks(k,j,i),Str1,k=1,5),Str1,j=1,Npindep) END DO RETURN 5 FORMAT(A10,1x,I4) 10 FORMAT(I4,1x,A1,I4,1x,A1,1x,1000000(5(F20.10,A1),A1)) END SUBROUTINE read_overlap_tracks_ascii SUBROUTINE write_overlap_tracks_ascii(funit, tstep, Nxp, Ntrcks, trcks) ! Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step IMPLICIT NONE INTEGER, INTENT(in) :: funit, tstep, Nxp, Ntrcks REAL(r_k), DIMENSION(5,Nxp,Ntrcks) :: trcks ! Local INTEGER :: i, j, k, ii, Npindep, Nrealtracks !!!!!!! Variables ! funit: unit where to write the file ! tstep: time-step to write the trajectories ! Nxp: maximum number of polygons per time-step ! Nrtcks: Number of trajectories of the given time-step ! trcks: trajectories fname = 'write_overlap_tracks_ascii' WRITE(funit,5)'time-step:', tstep ! Looking for the non-zero trajectories Nrealtracks = 0 DO i=1, Ntrcks Npindep = COUNT(trcks(2,:,i) /= zeroRK) IF (Npindep /= 0) Nrealtracks = Nrealtracks + 1 END DO WRITE(funit,*)Nrealtracks ! Only writting the trajectories with values ii = 1 DO i=1, Ntrcks Npindep = COUNT(trcks(2,:,i) /= zeroRK) IF (Npindep /= 0) THEN WRITE(funit,10)ii,';', Npindep, ';', ((trcks(k,j,i),',',k=1,5),':',j=1,Npindep) ii = ii + 1 END IF END DO RETURN 5 FORMAT(A10,1x,I4) 10 FORMAT(I4,1x,A1,I4,1x,A1,1x,1000000(5(F20.10,A1),A1)) END SUBROUTINE write_overlap_tracks_ascii SUBROUTINE read_overlap_polys_ascii(funit, tstep, Nxp, Nindep, SpIndep, NpIndep, pIndep) ! Subroutine to read from an ASCII file the associated polygons at a given time-step IMPLICIT NONE INTEGER, INTENT(in) :: funit, tstep, Nxp INTEGER, INTENT(out) :: Nindep INTEGER, DIMENSION(Nxp), INTENT(out) :: SpIndep, NpIndep INTEGER, DIMENSION(Nxp,Nxp), INTENT(out) :: pIndep ! Local INTEGER :: i, j, k INTEGER :: read_it !!!!!!! Variables ! funit: unit associated to the file ! tstep: time-step of the values ! Nxp: allowed maximum numbe of polygons per time-step ! Nindpe: Number of independent polygons at this time-step ! SpIndep: Associated polygon to the independent one from the previous time-step ! NpIndep: Number of associated polygons to the independent time-step ! pIndep: polygons associated to a given independent polygon fname = 'read_overlap_polys_ascii' Nindep = 0 SpIndep = 0 NpIndep = 0 READ(funit,5)Str10, read_it IF (read_it /= tstep) THEN WRITE(numSa,'(I4)')read_it WRITE(numSb,'(I4)')tstep msg = 'File time-step;' // TRIM(numSa) // ' does not coincide with the one from program:' // & TRIM(numSb) END IF READ(funit,*)Nindep DO i=1, Nindep READ(funit,10) k, Str1, SpIndep(i), Str1, NpIndep(i), Str1, (pIndep(i,j), Str1, j=1,NpIndep(i)) END DO RETURN 5 FORMAT(A10,1x,I4) 10 FORMAT(I4,1x,A1,1x,I4,1x,A1,1x,I4,A1,1x,100000(I4,A1)) END SUBROUTINE read_overlap_polys_ascii SUBROUTINE write_overlap_polys_ascii(funit, tstep, Nxp, Nindep, SpIndep, NpIndep, pIndep) ! Subroutine to write into an ASCII file the associated polygons at a given time-step IMPLICIT NONE INTEGER, INTENT(in) :: funit, tstep, Nxp, Nindep INTEGER, DIMENSION(Nindep), INTENT(in) :: SpIndep, NpIndep INTEGER, DIMENSION(Nindep,Nxp), INTENT(in) :: pIndep ! Local INTEGER :: i, j !!!!!!! Variables ! funit: unit associated to the file ! tstep: time-step of the values ! Nxp: allowed maximum numbe of polygons per time-step ! Nindpe: Number of independent polygons at this time-step ! SpIndep: Associated polygon to the independent one from the previous time-step ! NpIndep: Number of associated polygons to the independent time-step ! pIndep: polygons associated to a given independent polygon fname = 'write_overlap_polys_ascii' WRITE(funit,5)'time-step:', tstep WRITE(funit,*)Nindep, ' ! Number of independent polygons' DO i=1, Nindep WRITE(funit,10) i, ';', SpIndep(i), ';', NpIndep(i), ':', (pIndep(i,j), ',', j=1,NpIndep(i)) END DO RETURN 5 FORMAT(A10,1x,I4) 10 FORMAT(I4,1x,A1,1x,I4,1x,A1,1x,I4,A1,1x,100000(I4,A1)) END SUBROUTINE write_overlap_polys_ascii SUBROUTINE poly_overlap_tracks_area_ascii(dbg, compute, dx, dy, dt, minarea, inNallpolys, allpolys, & ctrpolys, areapolys, Nmaxpoly, Nmaxtracks, methodmulti) ! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum ! overlaping/coincidence filtrered by a minimal area writting theoutput on an ASCII file (memory limitations) IMPLICIT NONE LOGICAL, INTENT(in) :: dbg CHARACTER(LEN=*), INTENT(in) :: compute, methodmulti INTEGER, INTENT(in) :: dx, dy, dt, Nmaxpoly, Nmaxtracks INTEGER, DIMENSION(dt), INTENT(in) :: inNallpolys INTEGER, DIMENSION(dx,dy,dt), INTENT(in) :: allpolys REAL(r_k), INTENT(in) :: minarea REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in) :: ctrpolys REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in) :: areapolys ! Local INTEGER :: i, j, ip, it, iip, itt, iit INTEGER :: fprevunit, ftrackunit, ftrunit, ierr, ios LOGICAL :: file_exist, dooverlap, dotracks, doftracks REAL(r_k), DIMENSION(Nmaxpoly) :: Aprevpolys, Acurrpolys REAL(r_k), DIMENSION(2,Nmaxpoly) :: Cprevpolys, Ccurrpolys INTEGER, DIMENSION(dx,dy) :: meetpolys, searchpolys INTEGER, DIMENSION(Nmaxpoly) :: coincidencies INTEGER, DIMENSION(Nmaxpoly) :: prevID, currID REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,2) :: tracks REAL(r_k), DIMENSION(5,dt) :: finaltracks INTEGER, DIMENSION(:), ALLOCATABLE :: coins INTEGER, DIMENSION(:,:), ALLOCATABLE :: coinsNpts INTEGER :: Nmeet, Nsearch, Nindep INTEGER, DIMENSION(2) :: Nindeppolys, Npolystime CHARACTER(len=5) :: NcoinS INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,2) :: polysIndep INTEGER, DIMENSION(Nmaxpoly,2) :: NpolysIndep INTEGER, DIMENSION(Nmaxpoly,2) :: SpolysIndep INTEGER :: iindep, iiprev INTEGER :: Nprev, NNprev, Ntprev LOGICAL :: Indeppolychained INTEGER :: itrack, ictrack INTEGER :: ixp, iyp, ttrack INTEGER, DIMENSION(2) :: Ntracks INTEGER :: idtrack, maxtrack REAL(r_k), DIMENSION(5,Nmaxpoly,dt) :: singletrack REAL(r_k) :: totArea, dist, mindist, maxarea, areai !!!!!!! Variables ! dx,dy,dt: space/time dimensions ! compute: how to copmute ! 'scratch': everything from the beginning ! 'continue': skipt that parts which already have the ascii file written ! inNallpolys: Vector with the original number of polygons at each time-step ! allpolys: Series of 2D field with the polygons ! minarea: minimal area (in same units as areapolys) to perform the tracking ! ctrpolys: center of the polygons ! areapolys: area of the polygons ! Nmaxpoly: Maximum possible number of polygons ! Nmaxtracks: maximum number of tracks ! methodmulti: methodology to follow when multiple polygons are given for the same track ! 'mean': get coordinates from the areal-weighted mean of the centers of the given polygons and their areas ! 'largest': get the coorindates of the largest polygon ! 'closest': get the coordinates of the closest polygon fname = 'poly_overlap_tracks_area_ascii' IF (dbg) PRINT *,TRIM(fname) SELECT CASE (TRIM(compute)) CASE ('scratch') dooverlap = .TRUE. dotracks = .TRUE. doftracks = .TRUE. CASE ('continue') INQUIRE(file='polygons_overlap.dat', exist=file_exist) IF (.NOT.file_exist) THEN dooverlap = .TRUE. ELSE IF (dbg) THEN PRINT *, TRIM(warnmsg) PRINT *," "//TRIM(fname) // ": File 'polygons_overlap.dat' already exists, skipping it !!" END IF dooverlap = .FALSE. END IF INQUIRE(file='trajectories_overlap.dat', exist=file_exist) IF (.NOT.file_exist) THEN dotracks = .TRUE. ELSE IF (dbg) THEN PRINT *, TRIM(warnmsg) PRINT *, " " // TRIM(fname) // ": File 'trajectories_overlap.dat' already exists, " // & "skipping it !!" END IF dotracks = .FALSE. END IF INQUIRE(file='trajectories.dat', exist=file_exist) IF (.NOT.file_exist) THEN doftracks = .TRUE. ELSE IF (dbg) THEN PRINT *, TRIM(warnmsg) PRINT *," "//TRIM(fname) // ": File 'trajectories.dat' already exists, skipping it !!" END IF doftracks = .FALSE. END IF CASE DEFAULT msg = "compute case: '" // TRIM(compute) // "' not ready !!" CALL ErrMsg(msg, fname, -1) END SELECT ! Checking multi-polygon methodology IF ( (TRIM(methodmulti) /= 'mean') .AND. (TRIM(methodmulti) /= 'largest') .AND. & (TRIM(methodmulti) /= 'closest')) THEN msg= "methodology for multiple-polygons: '"//TRIM(methodmulti)//"' not ready" // NEW_LINE('a')//& " available ones: 'mean', 'largest', 'closest'" CALL ErrMsg(msg, fname, -1) END IF IF (dooverlap) THEN ! ASCII file for all the polygons and their previous associated one fprevunit = freeunit() OPEN(fprevunit, file='polygons_overlap.dat', status='new', form='formatted', iostat=ios) msg = "Problems opening file: 'polygons_overlap.dat'" IF (ios == 17) PRINT *," Careful: 'polygons_overlap.dat' already exists!!" CALL ErrMsg(msg, fname, ios) ! Number of independent polygons by time step Nindeppolys = 0 ! Number of polygons attached to each independent polygons by time step NpolysIndep = 0 ! ID of searching polygon attached to each independent polygons by time step SpolysIndep = 0 ! ID of polygons attached to each independent polygons by time step polysIndep = 0 ! ID of polygons from previous time-step prevID = 0 ! ID of polygons from current time-step currID = 0 ! First time-step all are independent polygons it = 1 Nmeet = inNallpolys(it) Nindeppolys(it) = Nmeet ip = 0 meetpolys = allpolys(:,:,it) DO i=1, Nmeet IF (areapolys(i,it) >= minarea) THEN ip = ip + 1 SpolysIndep(ip,it) = i currID(ip) = i Acurrpolys(ip) = areapolys(i,it) Ccurrpolys(1,ip) = ctrpolys(1,i,it) Ccurrpolys(2,ip) = ctrpolys(2,i,it) NpolysIndep(ip,it) = 1 polysIndep(ip,1,it) = i ELSE WHERE (meetpolys == i) meetpolys = 0 END WHERE END IF END DO Nindeppolys(1) = ip Npolystime(1) = ip ! Starting step it = 0 IF (dbg) THEN PRINT *,' time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',0 PRINT *,' number of independent polygons:', Nindeppolys(it+1) PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______' DO i=1, Nindeppolys(it+1) PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', & polysIndep(i,1:NpolysIndep(i,it+1),it+1) END DO END IF ! Writting to the ASCII file Independent polygons and their associated CALL write_overlap_polys_ascii(fprevunit,it+1, Nmaxpoly, Nindeppolys(it+1), & SpolysIndep(1:Nindeppolys(it+1),it+1), NpolysIndep(1:Nindeppolys(it+1),it+1), & polysIndep(1:Nindeppolys(it+1),:,it+1)) it = 1 ! Looking for the coincidencies at each time step DO iit=1, dt-1 ! Number of times that a polygon has a coincidence coincidencies = 0 ! Preparing for next time-step searchpolys = meetpolys prevID = 0 prevID = currID Aprevpolys = Acurrpolys Cprevpolys = Ccurrpolys Nmeet = inNallpolys(iit+1) meetpolys = allpolys(:,:,iit+1) ip = 0 DO i=1, Nmeet IF (areapolys(i,iit+1) >= minarea) THEN ip = ip + 1 currID(ip) = i Acurrpolys(ip) = areapolys(i,iit+1) Acurrpolys(ip) = areapolys(i,iit+1) Ccurrpolys(1,ip) = ctrpolys(1,i,iit+1) Ccurrpolys(2,ip) = ctrpolys(2,i,iit+1) ELSE WHERE (meetpolys == i) meetpolys = 0 END WHERE END IF END DO Nindeppolys(it+1) = ip Npolystime(it+1) = ip ! Looking throughout the independent polygons Nmeet = Nindeppolys(it+1) !Nsearch = Nindeppolys(it) ! Previous space might have more polygons that their number of independent ones Nsearch = Npolystime(it) IF (ALLOCATED(coins)) DEALLOCATE(coins) ALLOCATE(coins(Nmeet), STAT=ierr) msg="Problems allocating 'coins'" CALL ErrMsg(msg,fname,ierr) IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts) ALLOCATE(coinsNpts(Nmeet, Nsearch), STAT=ierr) msg="Problems allocating 'coinsNpts'" CALL ErrMsg(msg,fname,ierr) CALL coincidence_all_polys_area(dbg, dx, dy, Nmeet, currID, meetpolys, Ccurrpolys(:,1:Nmeet), & Nsearch, prevID, searchpolys, Cprevpolys(:,1:Nsearch), Aprevpolys(1:Nsearch), coins, & coinsNpts) ! Counting the number of times that a polygon has a coincidency IF (dbg) THEN PRINT *,' Coincidencies for the given time-step:', iit+1,' _______' DO i=1, Nmeet PRINT *,currID(i), coins(i),' N search pts:', coinsNpts(i,:) END DO END IF ! Looking for the same equivalencies Nindep = 0 DO i=1, Nmeet IF (coins(i) == -1) THEN Nindep = Nindep + 1 SpolysIndep(Nindep,it+1) = -1 NpolysIndep(Nindep,it+1) = NpolysIndep(Nindep,it+1) + 1 polysIndep(Nindep,NpolysIndep(Nindep,it+1),it+1) = currID(i) ELSE IF (coins(i) == -9) THEN WRITE(NcoinS,'(I5)')coins(i) msg="coins= "//TRIM(NcoinS)//" This is an error. One should have always only one " // & "coincidence of polygon" CALL ErrMsg(msg, fname, -1) ELSE ! Looking for coincidences with previous independent polygons DO ip=1, Nsearch ! Looking into the polygons associated NNprev = NpolysIndep(ip,it) DO j=1, NNprev IF (coins(i) == polysIndep(ip,j,it)) THEN ! Which index corresponds to this coincidence? iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i)) IF (iindep == -1) THEN Nindep = Nindep + 1 SpolysIndep(Nindep,it+1) = coins(i) END IF iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i)) IF (iindep < 0) THEN PRINT *,' Looking for:', coins(i) PRINT *,' Within:', SpolysIndep(1:Nindep,it+1) PRINT *,' Might content:', polysIndep(ip,1:NNprev,it) PRINT *,' From an initial list:', coins(1:Nmeet) msg = 'Wrong index! There must be an index here' CALL ErrMsg(msg,fname,iindep) END IF coincidencies(ip) = coincidencies(ip) + 1 NpolysIndep(iindep,it+1) = NpolysIndep(iindep,it+1) + 1 polysIndep(iindep,NpolysIndep(iindep,it+1),it+1) = currID(i) EXIT END IF END DO END DO END IF END DO Nindeppolys(it+1) = Nindep IF (dbg) THEN PRINT *,' time step:',iit+1,' number to look polygons:', Nmeet,' searching polygons:',Nsearch PRINT *,' number of independent polygons:', Nindeppolys(it+1) PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______' DO i=1, Nindeppolys(it+1) PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', & polysIndep(i,1:NpolysIndep(i,it+1),it+1) END DO END IF ! Writting to the ASCII file Independent polygons and their associated CALL write_overlap_polys_ascii(fprevunit, iit+1, Nmaxpoly, Nindeppolys(it+1), & SpolysIndep(1:Nindeppolys(it+1),it+1), NpolysIndep(1:Nindeppolys(it+1),it+1), & polysIndep(1:Nindeppolys(it+1),:,it+1)) ! Preparing for the next time-step SpolysIndep(:,it) = 0 NpolysIndep(:,it) = 0 polysIndep(:,:,it) = 0 Nindeppolys(it) = Nindeppolys(it+1) SpolysIndep(1:Nindeppolys(it),it) = SpolysIndep(1:Nindeppolys(it+1),it+1) NpolysIndep(1:Nindeppolys(it),it) = NpolysIndep(1:Nindeppolys(it+1),it+1) Npolystime(it) = Npolystime(it+1) DO ip=1, Nindeppolys(it) polysIndep(ip,1,it) = polysIndep(ip,1,it+1) polysIndep(ip,2,it) = polysIndep(ip,2,it+1) END DO SpolysIndep(:,it+1) = 0 NpolysIndep(:,it+1) = 0 polysIndep(:,:,it+1) = 0 END DO CLOSE(fprevunit) IF (dbg) PRINT *," Succesful writting of ASCII chain of polygons 'polygons_overlap.dat' !!" END IF ! ASCII file for all the polygons and their previous associated one fprevunit = freeunit() OPEN(fprevunit, file='polygons_overlap.dat', status='old', form='formatted', iostat=ios) msg = "Problems opening file: 'polygons_overlap.dat'" CALL ErrMsg(msg, fname, ios) it = 1 IF (dbg) THEN PRINT *, 'Coincidencies to connect _______' DO iit=1, dt ! Reading from the ASCII file Independent polygons and their associated CALL read_overlap_polys_ascii(fprevunit, iit, Nmaxpoly, Nindeppolys(it), SpolysIndep(:,it), & NpolysIndep(:,it), polysIndep(:,:,it)) PRINT *,' it:', iit, ' Nindep:', Nindeppolys(it) PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs' DO ip=1, Nindeppolys(it) PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':', & polysIndep(ip,1:NpolysIndep(ip,it),it) END DO END DO END IF REWIND(fprevunit) ! Trajectories ! It should be done following the number of 'independent' polygons ! One would concatenate that independent polygons which share IDs from one step to another IF (dotracks) THEN ! ASCII file for the trajectories ftrackunit = freeunit() OPEN(ftrackunit, file='trajectories_overlap.dat', status='new', form='formatted', iostat=ios) msg = "Problems opening file: 'trajectories_overlap.dat'" IF (ios == 17) PRINT *," Careful: 'trajectories_overlap.dat' already exists!!" CALL ErrMsg(msg,fname,ios) ! First time-step. Take all polygons itrack = 0 tracks = zeroRK Ntracks = 0 it = 1 iit = 1 CALL read_overlap_polys_ascii(fprevunit, iit, Nmaxpoly, Nindeppolys(it), SpolysIndep(:,it), & NpolysIndep(:,it), polysIndep(:,:,it)) DO ip=1, Nindeppolys(1) itrack = itrack + 1 tracks(1,1,itrack,1) = itrack*1. tracks(2,1,itrack,1) = SpolysIndep(ip,1) tracks(3,1,itrack,1) = ctrpolys(1,ip,1) tracks(4,1,itrack,1) = ctrpolys(2,ip,1) tracks(5,1,itrack,1) = 1 Ntracks(1) = Ntracks(1) + 1 END DO ! Writting first time-step trajectories to the intermediate file CALL write_overlap_tracks_ascii(ftrackunit,iit,Nmaxpoly, Ntracks(it), tracks(:,:,1:Ntracks(it),it)) ! Looping allover already assigned tracks it = 2 maxtrack = Ntracks(1) timesteps: DO iit=2, dt CALL read_overlap_polys_ascii(fprevunit, iit, Nmaxpoly, Nindeppolys(it), SpolysIndep(:,it), & NpolysIndep(:,it), polysIndep(:,:,it)) IF (dbg) PRINT *,'track-timestep:', iit, 'N indep polys:', Nindeppolys(it) ! Indep polygons current time-step current_poly: DO i=1, Nindeppolys(it) IF (dbg) PRINT *,' curent poly:', i, 'Prev poly:', SpolysIndep(i,it), ' N ass. polygons:', & NpolysIndep(i,it), 'ass. poly:', polysIndep(i,1:NpolysIndep(i,it),it) Indeppolychained = .FALSE. ! Number of tracks previous time-step ! Looping overall it1_tracks: DO itt=1, Ntracks(it-1) itrack = tracks(1,1,itt,it-1) ! Number polygons ID assigned Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0) IF (dbg) PRINT *,itt,' track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1) ! Looking for coincidencies DO iip=1, Ntprev IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN Indeppolychained = .TRUE. IF (dbg) PRINT *,' coincidence found by polygon:', tracks(2,iip,itt,it-1) EXIT END IF END DO IF (Indeppolychained) THEN Ntracks(it) = Ntracks(it) + 1 ictrack = Ntracks(it) ! Assigning all the IDs to the next step of the track DO iip=1, NpolysIndep(i,it) iiprev = polysIndep(i,iip,it) tracks(1,iip,ictrack,it) = itrack tracks(2,iip,ictrack,it) = iiprev tracks(3,iip,ictrack,it) = ctrpolys(1,iiprev,iit) tracks(4,iip,ictrack,it) = ctrpolys(2,iiprev,iit) tracks(5,iip,ictrack,it) = iit END DO EXIT END IF IF (Indeppolychained) EXIT END DO it1_tracks ! Creation of a new track IF (.NOT.Indeppolychained) THEN Ntracks(it) = Ntracks(it) + 1 ictrack = Ntracks(it) ! ID of new track maxtrack = maxtrack + 1 IF (dbg) PRINT *,' New track!', maxtrack ! Assigning all the IDs to the next step of the track DO j=1, NpolysIndep(i,it) iiprev = polysIndep(i,j,it) tracks(1,j,ictrack,it) = maxtrack tracks(2,j,ictrack,it) = iiprev tracks(3,j,ictrack,it) = ctrpolys(1,iiprev,iit) tracks(4,j,ictrack,it) = ctrpolys(2,iiprev,iit) tracks(5,j,ictrack,it) = iit END DO END IF END DO current_poly IF (dbg) THEN PRINT *,' At this time-step:', iit, ' N trajectories:', Ntracks(it) DO i=1, Ntracks(it) Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0) PRINT *,i ,'ID tracks:', tracks(1,1,i,it), 'ID polygon:', tracks(2,1:Nprev,i,it) END DO END IF CALL write_overlap_tracks_ascii(ftrackunit,iit,Nmaxpoly,Ntracks(it),tracks(:,:,1:Ntracks(it),it)) ! Re-initializing for the next time-step tracks(:,:,:,it-1) = zeroRK Ntracks(it-1) = Ntracks(it) tracks(:,:,1:Ntracks(it-1),it-1) = tracks(:,:,1:Ntracks(it),it) Ntracks(it) = 0 tracks(:,:,:,it) = zeroRK END DO timesteps CLOSE(ftrackunit) IF (dbg) PRINT *," Succesful writting of ASCII chain of polygons 'trajectories_overlap.dat' !!" CLOSE(fprevunit) END IF ! Summarizing trajectories ! When multiple polygons are available, the mean of their central positions determines the position IF (doftracks) THEN ! ASCII file for the trajectories ftrackunit = freeunit() OPEN(ftrackunit, file='trajectories_overlap.dat', status='old', form='formatted', iostat=ios) msg = "Problems opening file: 'trajectories_overlap.dat'" CALL ErrMsg(msg,fname,ios) ! ASCII file for the final trajectories ftrunit = freeunit() OPEN(ftrunit, file='trajectories.dat', status='new', form='formatted', iostat=ios) msg = "Problems opening file: 'trajectories.dat'" IF (ios == 17) PRINT *," Careful: 'trajectories.dat' already exists!!" CALL ErrMsg(msg,fname,ios) finaltracks = zeroRK DO itt=1, Nmaxtracks CALL read_overlap_single_track_ascii(ftrackunit, dt, Nmaxpoly, Nmaxtracks, itt, singletrack) ! It might reach the las trajectory IF (ALL(singletrack == zeroRK)) EXIT itrack = INT(MAXVAL(singletrack(1,1,:))) IF (dbg) THEN PRINT *,' Trajectory:', itt, '_______', itrack DO it=1, dt IF (singletrack(2,1,it) /= zeroRK) THEN j = COUNT(singletrack(2,:,it) /= zeroRK) PRINT *,it,':',(singletrack(3,i,it),',',singletrack(4,i,it),' ; ',i=1,j) END IF END DO END IF finaltracks = zeroRK finaltracks(1,:) = itrack*oneRK DO it =1, dt Nprev = COUNT(INT(singletrack(2,:,it)) /= zeroRK) IF (Nprev /= 0) THEN finaltracks(5,it) = it*oneRK IF (TRIM(methodmulti) == 'largest') THEN maxarea = -10.*oneRK DO ip=1, Nprev IF (areapolys(singletrack(2,ip,it),it) > maxarea) THEN maxarea = areapolys(singletrack(2,ip,it),it) i = ip END IF END DO IF (dbg) THEN PRINT *,' Determine the trajectory coordinates to the largest polygon:', i, & ' area:', maxarea END IF finaltracks(2,it) = singletrack(2,i,it)*oneRK finaltracks(3,it) = singletrack(3,i,it) finaltracks(4,it) = singletrack(4,i,it) ELSE IF (TRIM(methodmulti) == 'closest') THEN IF (it > 1) THEN mindist = 10000000.*oneRK DO ip=1, Nprev dist = SQRT((singletrack(3,ip,it)-finaltracks(3,it-1))**2 + & (singletrack(4,ip,it)-finaltracks(4,it-1))**2 ) IF (dist < mindist) THEN mindist = dist i = ip END IF END DO finaltracks(2,it) = singletrack(3,i,it)*oneRK finaltracks(3,it) = singletrack(3,i,it) finaltracks(4,it) = singletrack(4,i,it) IF (dbg) THEN PRINT *,' Determine the trajectory coordinates to the closest previous polygon:',i,& ' distance:', mindist END IF ELSE maxarea = -10.*oneRK DO ip=1, Nprev IF (areapolys(singletrack(2,ip,it),it) > maxarea) THEN maxarea = areapolys(singletrack(2,ip,it),it) i = ip END IF END DO IF (dbg) THEN PRINT *, ' Determine the trajectory coordinates to the largest polygon:', i, & ' area:', maxarea, ' at the first time-step then to the closest' END IF finaltracks(2,it) = i*oneRK finaltracks(3,it) = singletrack(3,i,it) finaltracks(4,it) = singletrack(4,i,it) END IF ELSE totArea = zeroRK finaltracks(2,it) = -oneRK finaltracks(3,it) = zeroRK finaltracks(4,it) = zeroRK DO ip=1, Nprev areai = areapolys(singletrack(2,ip,it),it) totArea = totArea + areai finaltracks(3,it) = finaltracks(3,it) + singletrack(3,ip,it)*areai finaltracks(4,it) = finaltracks(4,it) + singletrack(4,ip,it)*areai END DO finaltracks(3,it) = finaltracks(3,it)/totArea finaltracks(4,it) = finaltracks(4,it)/totArea IF (dbg) THEN PRINT *,' Determine the trajectory coordinates to the area-averaged polygon ' // & ' total area:', totArea END IF END IF END IF END DO ! Writting the final track into the ASCII file CALL write_finaltrack_ascii(ftrunit, dt, finaltracks) END DO CLOSE(ftrackunit) IF (dbg) PRINT *," Succesful writting of ASCII trajectories 'trajectories.dat' !!" CLOSE(ftrunit) END IF IF (ALLOCATED(coins)) DEALLOCATE(coins) IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts) RETURN END SUBROUTINE poly_overlap_tracks_area_ascii SUBROUTINE poly_overlap_tracks_area(dbg, dx, dy, dt, minarea, inNallpolys, allpolys, ctrpolys, & areapolys, Nmaxpoly, Nmaxtracks, tracks, finaltracks) ! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum ! overlaping/coincidence filtrered by a minimal area IMPLICIT NONE LOGICAL, INTENT(in) :: dbg INTEGER, INTENT(in) :: dx, dy, dt, Nmaxpoly, Nmaxtracks INTEGER, DIMENSION(dt), INTENT(in) :: inNallpolys INTEGER, DIMENSION(dx,dy,dt), INTENT(in) :: allpolys REAL(r_k), INTENT(in) :: minarea REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in) :: ctrpolys REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in) :: areapolys REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,dt), & INTENT(out) :: tracks REAL(r_k), DIMENSION(4,Nmaxtracks,dt), INTENT(out) :: finaltracks ! Local INTEGER :: i, j, ip, it, iip, itt INTEGER :: ierr REAL(r_k), DIMENSION(Nmaxpoly) :: Aprevpolys, Acurrpolys REAL(r_k), DIMENSION(2,Nmaxpoly) :: Cprevpolys, Ccurrpolys INTEGER, DIMENSION(dt) :: Nallpolys INTEGER, DIMENSION(dx,dy) :: meetpolys, searchpolys INTEGER, DIMENSION(Nmaxpoly) :: coincidencies INTEGER, DIMENSION(Nmaxpoly) :: prevID, currID INTEGER, DIMENSION(:), ALLOCATABLE :: coins INTEGER, DIMENSION(:,:), ALLOCATABLE :: coinsNpts INTEGER :: Nmeet, Nsearch, Nindep INTEGER, DIMENSION(dt) :: Nindeppolys CHARACTER(len=5) :: NcoinS INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt) :: polysIndep INTEGER, DIMENSION(Nmaxpoly,dt) :: NpolysIndep INTEGER, DIMENSION(Nmaxpoly,dt) :: SpolysIndep INTEGER :: iindep, iiprev INTEGER :: Nprev, NNprev, Ntprev LOGICAL :: Indeppolychained INTEGER :: itrack, ictrack REAL(r_k) :: ixp, iyp INTEGER :: ttrack INTEGER, DIMENSION(dt) :: Ntracks INTEGER :: idtrack, maxtrack !!!!!!! Variables ! dx,dy,dt: space/time dimensions ! Nallpolys: Vector with the number of polygons at each time-step ! allpolys: Series of 2D field with the polygons ! minarea: minimal area (in same units as areapolys) to perform the tracking ! ctrpolys: center of the polygons ! areapolys: area of the polygons ! Nmaxpoly: Maximum possible number of polygons ! Nmaxtracks: maximum number of tracks ! tracks: series of consecutive polygons ! trackperiod: period of the track in time-steps fname = 'poly_overlap_tracks_area' IF (dbg) PRINT *,TRIM(fname) ! Number of independent polygons by time step Nindeppolys = 0 ! Number of polygons attached to each independent polygons by time step NpolysIndep = 0 ! ID of searching polygon attached to each independent polygons by time step SpolysIndep = 0 ! ID of polygons attached to each independent polygons by time step polysIndep = 0 ! ID of polygons from previous time-step prevID = 0 ! ID of polygons from current time-step currID = 0 ! First time-step all are independent polygons it = 1 Nmeet = inNallpolys(it) Nindeppolys(it) = Nmeet ip = 0 meetpolys = allpolys(:,:,it) DO i=1, Nmeet IF (areapolys(i,it) >= minarea) THEN ip = ip + 1 SpolysIndep(ip,it) = i currID(ip) = i Acurrpolys(ip) = areapolys(i,it) Ccurrpolys(1,ip) = ctrpolys(1,i,it) Ccurrpolys(2,ip) = ctrpolys(2,i,it) NpolysIndep(ip,it) = 1 polysIndep(ip,1,it) = i ELSE WHERE (meetpolys == i) meetpolys = 0 END WHERE END IF END DO Nallpolys(1) = ip Nindeppolys(1) = ip ! Starting step it = 0 IF (dbg) THEN PRINT *,' time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',0 PRINT *,' number of independent polygons:', Nindeppolys(it+1) PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______' DO i=1, Nindeppolys(it+1) PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', & polysIndep(i,1:NpolysIndep(i,it+1),it+1) END DO END IF ! Looking for the coincidencies at each time step DO it=1, dt-1 ! Number of times that a polygon has a coincidence coincidencies = 0 Nmeet = inNallpolys(it+1) searchpolys = meetpolys meetpolys = allpolys(:,:,it+1) prevID = 0 prevID = currID Aprevpolys = Acurrpolys Cprevpolys = Ccurrpolys ip = 0 DO i=1, Nmeet IF (areapolys(i,it+1) >= minarea) THEN ip = ip + 1 currID(ip) = i Acurrpolys(ip) = areapolys(i,it+1) Acurrpolys(ip) = areapolys(i,it+1) Ccurrpolys(1,ip) = ctrpolys(1,i,it+1) Ccurrpolys(2,ip) = ctrpolys(2,i,it+1) ELSE WHERE (meetpolys == i) meetpolys = 0 END WHERE END IF END DO Nallpolys(it+1) = ip Nindeppolys(it+1) = ip Nmeet = Nallpolys(it+1) ! Looking throughout the independent polygons Nsearch = Nindeppolys(it) IF (ALLOCATED(coins)) DEALLOCATE(coins) ALLOCATE(coins(Nmeet), STAT=ierr) msg="Problems allocating 'coins'" CALL ErrMsg(msg,fname,ierr) IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts) ALLOCATE(coinsNpts(Nmeet, Nsearch), STAT=ierr) msg="Problems allocating 'coinsNpts'" CALL ErrMsg(msg,fname,ierr) CALL coincidence_all_polys_area(dbg, dx,dy, Nmeet, currID, meetpolys, Acurrpolys(1:Nmeet), & Nsearch, prevID, searchpolys, Cprevpolys(:,1:Nsearch), Aprevpolys(1:Nsearch), coins, & coinsNpts) ! Counting the number of times that a polygon has a coincidency IF (dbg) THEN PRINT *,' Coincidencies for the given time-step:', it+1,' _______' DO i=1, Nmeet PRINT *,currID(i), coins(i),' N search pts:', coinsNpts(i,:) END DO END IF ! Looking for the same equivalencies Nindep = 0 DO i=1, Nmeet IF (coins(i) == -1) THEN Nindep = Nindep + 1 SpolysIndep(Nindep,it+1) = -1 NpolysIndep(Nindep,it+1) = NpolysIndep(Nindep,it+1) + 1 polysIndep(Nindep,NpolysIndep(Nindep,it+1),it+1) = currID(i) ELSE IF (coins(i) == -9) THEN WRITE(NcoinS,'(I5)')coins(i) msg="coins= "//TRIM(NcoinS)//" This is an error. One should have always only one " // & "coincidence of polygon" CALL ErrMsg(msg, fname, -1) ELSE ! Looking for coincidences with previous independent polygons DO ip=1, Nsearch ! Looking into the polygons associated NNprev = NpolysIndep(ip,it) DO j=1, NNprev IF (coins(i) == polysIndep(ip,j,it)) THEN ! Which index corresponds to this coincidence? iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i)) IF (iindep == -1) THEN Nindep = Nindep + 1 SpolysIndep(Nindep,it+1) = coins(i) END IF iindep = Index1DArrayI(SpolysIndep(1:Nindep,it+1), Nindep, coins(i)) IF (iindep < 0) THEN PRINT *,' Looking for:', coins(i) PRINT *,' Within:', SpolysIndep(1:Nindep,it+1) PRINT *,' Might content:', polysIndep(ip,1:NNprev,it) PRINT *,' From an initial list:', coins(1:Nmeet) msg = 'Wrong index! There must be an index here' CALL ErrMsg(msg,fname,iindep) END IF coincidencies(ip) = coincidencies(ip) + 1 NpolysIndep(iindep,it+1) = NpolysIndep(iindep,it+1) + 1 polysIndep(iindep,NpolysIndep(iindep,it+1),it+1) = currID(i) EXIT END IF END DO END DO END IF END DO Nindeppolys(it+1) = Nindep IF (dbg) THEN PRINT *,' time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',Nsearch PRINT *,' number of independent polygons:', Nindeppolys(it+1) PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______' DO i=1, Nindeppolys(it+1) PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', & polysIndep(i,1:NpolysIndep(i,it+1),it+1) END DO END IF END DO IF (dbg) THEN PRINT *, 'Coincidencies to connect _______' DO it=1, dt PRINT *,' it:', it, ' Nindep:', Nindeppolys(it) PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs' DO ip=1, Nindeppolys(it) PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':', & polysIndep(ip,1:NpolysIndep(ip,it),it) END DO END DO END IF ! Trajectories ! It should be done following the number of 'independent' polygons ! One would concatenate that independent polygons which share IDs from one step to another ! First time-step. Take all polygons itrack = 0 tracks = 0. Ntracks = 0 DO ip=1, Nindeppolys(1) itrack = itrack + 1 tracks(1,1,itrack,1) = itrack*1. tracks(2,1,itrack,1) = SpolysIndep(ip,1) tracks(3,1,itrack,1) = ctrpolys(1,ip,1) tracks(4,1,itrack,1) = ctrpolys(2,ip,1) tracks(5,1,itrack,1) = 1 Ntracks(1) = Ntracks(1) + 1 END DO ! Looping allover already assigned tracks timesteps: DO it=2, dt IF (dbg) PRINT *,'track-timestep:', it, 'N indep polys:', Nindeppolys(it) ! Indep polygons current time-step current_poly: DO i=1, Nindeppolys(it) IF (dbg) PRINT *,' curent poly:', i, 'Prev poly:', SpolysIndep(i,it), ' N ass. polygons:', & NpolysIndep(i,it), 'ass. poly:', polysIndep(i,1:NpolysIndep(i,it),it) Indeppolychained = .FALSE. ! Number of tracks previous time-step ! Looping overall it1_tracks: DO itt=1, Ntracks(it-1) itrack = tracks(1,1,itt,it-1) ! Number polygons ID assigned Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0) IF (dbg) PRINT *,itt,' track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1) ! Looking for coincidencies DO iip=1, Ntprev IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN Indeppolychained = .TRUE. IF (dbg) PRINT *,' coincidence found by polygon:', tracks(2,iip,itt,it-1) EXIT END IF END DO IF (Indeppolychained) THEN Ntracks(it) = Ntracks(it) + 1 ictrack = Ntracks(it) ! Assigning all the IDs to the next step of the track DO iip=1, NpolysIndep(i,it) iiprev = polysIndep(i,iip,it) tracks(1,iip,ictrack,it) = itrack tracks(2,iip,ictrack,it) = iiprev ixp = ctrpolys(1,iiprev,it) iyp = ctrpolys(2,iiprev,it) tracks(3,iip,ictrack,it) = ixp tracks(4,iip,ictrack,it) = iyp tracks(5,iip,ictrack,it) = it END DO EXIT END IF IF (Indeppolychained) EXIT END DO it1_tracks ! Creation of a new track IF (.NOT.Indeppolychained) THEN Ntracks(it) = Ntracks(it) + 1 ictrack = Ntracks(it) ! ID of new track maxtrack = INT(MAXVAL(tracks(1,:,:,:)*1.)) IF (dbg) PRINT *,' New track!', maxtrack+1 ! Assigning all the IDs to the next step of the track DO j=1, NpolysIndep(i,it) iiprev = polysIndep(i,j,it) tracks(1,j,ictrack,it) = maxtrack+1 tracks(2,j,ictrack,it) = iiprev ixp = ctrpolys(1,iiprev,it) iyp = ctrpolys(2,iiprev,it) tracks(3,j,ictrack,it) = ixp tracks(4,j,ictrack,it) = iyp tracks(5,j,ictrack,it) = it END DO END IF END DO current_poly IF (dbg) THEN PRINT *,' At this time-step:', it, ' N trajectories:', Ntracks(it) DO i=1, Ntracks(it) Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0) PRINT *,i ,'ID tracks:', tracks(1,1,i,it), 'ID polygon:', tracks(2,1:Nprev,i,it) END DO END IF END DO timesteps ! Summarizing trajectories ! When multiple polygons are available, the mean of their central positions determines the position finaltracks = 0. maxtrack = MAXVAL(tracks(1,:,:,:)) DO it=1, dt DO itt=1, Ntracks(it) itrack = INT(tracks(1,1,itt,it)) Nprev = COUNT(INT(tracks(2,:,itt,it)) /= 0) finaltracks(1,itrack,it) = itrack*1. finaltracks(2,itrack,it) = SUM(tracks(3,:,itt,it))/Nprev*1. finaltracks(3,itrack,it) = SUM(tracks(4,:,itt,it))/Nprev*1. finaltracks(4,itrack,it) = it*1. END DO END DO DEALLOCATE(coins) DEALLOCATE(coinsNpts) RETURN END SUBROUTINE poly_overlap_tracks_area SUBROUTINE coincidence_all_polys_area(dbg, dx, dy, Nallpoly, IDallpoly, allpoly, icpolys, Npoly, & IDpolys, polys, cpolys, apolys, polycoins, coinNptss) ! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons ! In case of multiple coincidencies, the closest and then the largest is taken filtrered by a minimal area ! Here the difference is that the index does not coincide with its ID IMPLICIT NONE LOGICAL, INTENT(in) :: dbg INTEGER, INTENT(in) :: dx, dy, Nallpoly, Npoly INTEGER, DIMENSION(dx,dy), INTENT(in) :: allpoly, polys INTEGER, DIMENSION(Nallpoly), INTENT(in) :: IDallpoly INTEGER, DIMENSION(Npoly), INTENT(in) :: IDpolys REAL(r_k), DIMENSION(2,Nallpoly), INTENT(in) :: icpolys REAL(r_k), DIMENSION(2,Npoly), INTENT(in) :: cpolys REAL(r_k), DIMENSION(Npoly), INTENT(in) :: apolys INTEGER, DIMENSION(Nallpoly), INTENT(out) :: polycoins INTEGER, DIMENSION(Nallpoly,Npoly), INTENT(out) :: coinNptss ! Local INTEGER :: i, j, ip INTEGER :: maxcorr INTEGER :: Nmaxcorr LOGICAL, DIMENSION(dx,dy) :: boolpoly INTEGER :: maxcoin REAL :: dist, maxcoindist, maxcoinarea INTEGER, DIMENSION(Npoly) :: IDcoins !!!!!!! Variables ! dx,dy: dimension of the space ! Nallpoly: Number of polygons to find coincidence ! allpoly: space with the polygons to meet ! IDallpoly: ID of the polygon to find coincidence ! icpolys: center of the polygons to look for the coincidence ! Npoly: number of polygons on the 2D space ! polys: 2D field of polygons identified by their integer number (0 for no polygon) ! IDpolys: ID of the polygon to search for coincidences ! cpolys: center of the polygons ! apolys: area of the polygons ! polycoins: coincident polyogn ! -1: no-coincidence ! 1 < Npoly: single coincidence with a given polygon ! -9: coincidence with more than one polygon ! coinNptss: number of points coincident with each polygon fname = 'coincidence_all_polys_area' IF (dbg) PRINT *,TRIM(fname) DO ip=1, Nallpoly boolpoly = allpoly == IDallpoly(ip) CALL coincidence_poly_area(dbg, dx, dy, boolpoly, Npoly, polys, polycoins(ip), IDcoins, & coinNptss(ip,:)) IF (dbg) PRINT *,' polygon', IDallpoly(ip), ' coincidence with:', polycoins(ip), 'IDpolys:', IDpolys(1:Npoly) ! Coincidence with more than one polygon IF (polycoins(ip) == -9) THEN maxcoindist = -10. maxcoinarea = -10. maxcoin = MAXVAL(coinNptss(ip,:)) DO j=1, Npoly IF (coinNptss(ip,j) == maxcoin) THEN dist = SQRT( (icpolys(1,ip)*1.-cpolys(1,j)*1.)**2 + (icpolys(2,ip)*1.-cpolys(2,j)*1.)**2 ) IF ( dist > maxcoindist) THEN maxcoindist = dist maxcoinarea = apolys(j) polycoins(ip) = IDcoins(j) ELSE IF ( dist == maxcoindist) THEN IF (apolys(j) > maxcoinarea) THEN polycoins(ip) = IDcoins(j) maxcoinarea = apolys(j) END IF END IF END IF END DO END IF END DO RETURN END SUBROUTINE coincidence_all_polys_area SUBROUTINE coincidence_poly_area(dbg, dx, dy, poly, Npoly, polys, polycoin, IDpoly, coinNpts) ! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons ! Here the difference is that the index does not coincide with its ID IMPLICIT NONE LOGICAL, INTENT(in) :: dbg INTEGER, INTENT(in) :: dx, dy, Npoly LOGICAL, DIMENSION(dx,dy), INTENT(in) :: poly INTEGER, DIMENSION(dx,dy), INTENT(in) :: polys INTEGER, INTENT(out) :: polycoin INTEGER, DIMENSION(Npoly), INTENT(out) :: IDpoly, coinNpts ! Local INTEGER :: i, j, ip INTEGER :: maxcorr INTEGER :: Nmaxcorr ! Lluis INTEGER :: Ndiffvals INTEGER, DIMENSION(:), ALLOCATABLE :: diffvals !!!!!!! Variables ! dx,dy: dimension of the space ! poly: bolean polygon to meet ! Npoly: number of polygons on the 2D space ! polys: 2D field of polygons identified by their integer number (0 for no polygon) ! polycoin: coincident polyogn ! -1: no-coincidence ! 1 < Npoly: single coincidence with a given polygon ! -9: coincidence with more than one polygon ! IDpoly: ID of the found polygon ! coinNpts: number of points coincident with each polygon fname = 'coincidence_poly_area' IF (dbg) PRINT *,TRIM(fname) IF (dbg) THEN PRINT *,' Boolean polygon to search coincidences ...' DO i=1,dx PRINT *,poly(i,:) END DO PRINT *,' 2D polygons space ...' DO i=1,dx PRINT '(1000(I7,1x))',polys(i,:) END DO END IF IF (ALLOCATED(diffvals)) DEALLOCATE(diffvals) ALLOCATE(diffvals(dx*dy)) ! Checking for consistency on number of polygons and real content (except 0 value) CALL Nvalues_2DArrayI(dx, dy, dx*dy, polys, Ndiffvals, diffvals) IF (Ndiffvals -1 /= Npoly) THEN PRINT *,TRIM(emsg) PRINT *,' number of different values:', Ndiffvals-1, ' theoretical Npoly:', Npoly PRINT *,' Different values:', diffvals(1:Ndiffvals) msg = 'Number of different values and Npoly must coincide' CALL ErrMsg(msg, fname, -1) END IF ! Looking for coincient points for the polygon coinNpts = 0 IDpoly = 0 ip = 0 DO i=1,dx DO j=1,dy IF (poly(i,j) .AND. polys(i,j) .NE. 0) THEN IF (.NOT.ANY(IDpoly == polys(i,j))) THEN ip = ip + 1 IDpoly(ip) = polys(i,j) ELSE ip = Index1DarrayI(IDpoly, Npoly, polys(i,j)) END IF coinNpts(ip) = coinNpts(ip) + 1 END IF END DO END DO maxcorr = 0 maxcorr = INT(MAXVAL(coinNpts*1.)) IF (dbg) PRINT *,' Maximum coincidence:', maxcorr IF (maxcorr == 0) THEN polycoin = -1 ELSE Nmaxcorr = 0 DO ip=1, Npoly IF (coinNpts(ip) == maxcorr) THEN Nmaxcorr = Nmaxcorr+1 polycoin = IDpoly(ip) END IF END DO IF (Nmaxcorr > 1) polycoin = -9 END IF IF (dbg) THEN PRINT *,' Coincidences for each polygon _______', Npoly DO ip=1, Npoly PRINT *,' ',ip, ' ID:', IDpoly(ip),': ', coinNpts(ip) END DO END IF RETURN END SUBROUTINE coincidence_poly_area SUBROUTINE poly_overlap_tracks(dbg, dx, dy, dt, minarea, Nallpolys, allpolys, ctrpolys, & areapolys, Nmaxpoly, Nmaxtracks, tracks, finaltracks) ! Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum overlaping/coincidence IMPLICIT NONE LOGICAL, INTENT(in) :: dbg INTEGER, INTENT(in) :: dx, dy, dt, Nmaxpoly, Nmaxtracks INTEGER, DIMENSION(dt), INTENT(in) :: Nallpolys INTEGER, DIMENSION(dx,dy,dt), INTENT(in) :: allpolys REAL(r_k), INTENT(in) :: minarea REAL(r_k), DIMENSION(2,Nmaxpoly,dt), INTENT(in) :: ctrpolys REAL(r_k), DIMENSION(Nmaxpoly,dt), INTENT(in) :: areapolys REAL(r_k), DIMENSION(5,Nmaxpoly,Nmaxtracks,dt), & INTENT(out) :: tracks REAL(r_k), DIMENSION(4,Nmaxtracks,dt), INTENT(out) :: finaltracks ! Local INTEGER :: i, j, ip, it, iip, itt INTEGER :: ierr INTEGER, DIMENSION(Nmaxpoly,dt) :: coincidencies, NOcoincidencies INTEGER, DIMENSION(:), ALLOCATABLE :: coins INTEGER, DIMENSION(:,:), ALLOCATABLE :: coinsNpts INTEGER, DIMENSION(Nmaxpoly,dt) :: polycoincidencies INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt) :: coincidenciesNpts INTEGER :: Nmeet, Nsearch, Nindep INTEGER, DIMENSION(dt) :: Nindeppolys CHARACTER(len=5) :: NcoinS INTEGER, DIMENSION(Nmaxpoly,Nmaxpoly,dt) :: polysIndep INTEGER, DIMENSION(Nmaxpoly,dt) :: NpolysIndep INTEGER, DIMENSION(Nmaxpoly,dt) :: SpolysIndep INTEGER :: iindep, iiprev INTEGER :: Nprev, NNprev, Ntprev LOGICAL :: Indeppolychained INTEGER :: itrack, ictrack INTEGER :: ixp, iyp, ttrack INTEGER, DIMENSION(dt) :: Ntracks INTEGER :: idtrack, maxtrack !!!!!!! Variables ! dx,dy,dt: space/time dimensions ! Nallpolys: Vector with the number of polygons at each time-step ! allpolys: Series of 2D field with the polygons ! minarea: minimal area (in same units as areapolys) to perform the tracking ! ctrpolys: center of the polygons ! areapolys: area of the polygons ! Nmaxpoly: Maximum possible number of polygons ! Nmaxtracks: maximum number of tracks ! tracks: series of consecutive polygons ! trackperiod: period of the track in time-steps fname = 'poly_overlap_tracks' IF (dbg) PRINT *,TRIM(fname) polycoincidencies = fillvalI coincidenciesNpts = fillvalI ! Number of times that a polygon has a coincidence coincidencies = 0 ! Polygons without a coincidence NOcoincidencies = 0 ! Number of independent polygons by time step Nindeppolys = 0 ! Number of polygons attached to each independent polygons by time step NpolysIndep = 0 ! ID of searching polygon attached to each independent polygons by time step SpolysIndep = 0 ! ID of polygons attached to each independent polygons by time step polysIndep = 0 ! First time-step all are independent polygons it = 1 Nmeet = Nallpolys(it) Nindeppolys(it) = Nmeet DO i=1, Nmeet SpolysIndep(i,it) = i NpolysIndep(1:Nmeet,it) = 1 polysIndep(1,i,it) = i END DO ! Looking for the coincidencies at each time step DO it=1, dt-1 Nmeet = Nallpolys(it+1) Nsearch = Nallpolys(it) IF (ALLOCATED(coins)) DEALLOCATE(coins) ALLOCATE(coins(Nmeet), STAT=ierr) msg="Problems allocating 'coins'" CALL ErrMsg(msg,fname,ierr) IF (ALLOCATED(coinsNpts)) DEALLOCATE(coinsNpts) ALLOCATE(coinsNpts(Nmeet, Nsearch), STAT=ierr) msg="Problems allocating 'coinsNpts'" CALL ErrMsg(msg,fname,ierr) CALL coincidence_all_polys(dbg, dx, dy, Nmeet, allpolys(:,:,it+1), ctrpolys(:,1:Nmeet,it+1), & Nsearch, allpolys(:,:,it), ctrpolys(:,1:Nsearch,it), areapolys(1:Nsearch,it), coins, coinsNpts) polycoincidencies(1:Nmeet,it+1) = coins coincidenciesNpts(1:Nmeet,1:Nsearch,it+1) = coinsNpts ! Counting the number of times that a polygon has a coincidency IF (dbg) THEN PRINT *,' Coincidencies for the given time-step:', it+1,' _______' DO i=1, Nmeet PRINT *,coins(i),' N search pts:', coinsNpts(i,:) END DO END IF Nindep = 0 DO i=1, Nmeet IF (coins(i) == -1) THEN Nindep = Nindep + 1 NOcoincidencies(i,it+1) = 1 SpolysIndep(Nindep,it+1) = -1 NpolysIndep(Nindep,it+1) = NpolysIndep(Nindep,it+1) + 1 polysIndep(Nindep,NpolysIndep(Nindep,it+1),it+1) = i ELSE IF (coins(i) == -9) THEN WRITE(NcoinS,'(I5)')coins(i) msg="coins= "//TRIM(NcoinS)//" This is an error. One should have always only one " // & "coincidence of polygon" CALL ErrMsg(msg, fname, -1) ELSE DO ip=1, Nsearch IF (coins(i) == ip) THEN IF (coincidencies(ip,it+1) == 0) THEN Nindep = Nindep + 1 SpolysIndep(Nindep,it+1) = ip END IF coincidencies(ip,it+1) = coincidencies(ip,it+1) + 1 DO iindep=1, Nindep IF (SpolysIndep(iindep,it+1) == coins(i)) THEN NpolysIndep(iindep,it+1) = NpolysIndep(iindep,it+1) + 1 polysIndep(iindep,NpolysIndep(iindep,it+1),it+1) = i END IF END DO END IF END DO END IF END DO Nindeppolys(it+1) = Nindep IF (dbg) THEN PRINT *,' time step:',it+1,' number to look polygons:', Nmeet,' searching polygons:',Nsearch PRINT *,' number of independent polygons:', Nindeppolys(it+1) PRINT *,' indep_polygon prev_step_polygon Nassociated_polygons curr_ass_polygons _______' DO i=1, Nindeppolys(it+1) PRINT *,i, SpolysIndep(i,it+1), NpolysIndep(i,it+1), ':', & polysIndep(i,1:NpolysIndep(i,it+1),it+1) END DO END IF END DO IF (dbg) THEN PRINT *, 'Coincidencies to connect _______' DO it=1, dt PRINT *,' it:', it, ' Nindep:', Nindeppolys(it) PRINT '(4x,3(A6,1x))','Nindep', 'PrevID', 'IDs' DO ip=1, Nindeppolys(it) PRINT '(4x,I6,A1,I6,A1,100(I6))', ip, ',', SpolysIndep(ip,it), ':', & polysIndep(ip,1:NpolysIndep(ip,it),it) END DO END DO END IF ! Trajectories ! It should be done following the number of 'independent' polygons ! One would concatenate that independent polygons which share IDs from one step to another ! First time-step. Take all polygons itrack = 0 tracks = 0. Ntracks = 0 DO ip=1, Nindeppolys(1) itrack = itrack + 1 tracks(1,1,itrack,1) = itrack*1. tracks(2,1,itrack,1) = SpolysIndep(ip,1) tracks(3,1,itrack,1) = ctrpolys(1,ip,1) tracks(4,1,itrack,1) = ctrpolys(2,ip,1) tracks(5,1,itrack,1) = 1 Ntracks(1) = Ntracks(1) + 1 END DO ! Looping allover already assigned tracks timesteps: DO it=2, dt IF (dbg) PRINT *,'timestep:', it, 'N indep polys:', Nindeppolys(it) ! Indep polygons current time-step current_poly: DO i=1, Nindeppolys(it) IF (dbg) PRINT *,' curent poly:', i, 'Prev poly:', SpolysIndep(i,it), ' N ass. polygons:', & NpolysIndep(i,it), 'ass. poly:', polysIndep(i,1:NpolysIndep(i,it),it) Indeppolychained = .FALSE. ! Number of tracks previous time-step ! Looping overall it1_tracks: DO itt=1, Ntracks(it-1) itrack = tracks(1,1,itt,it-1) ! Number polygons ID assigned Ntprev = COUNT(tracks(2,:,itt,it-1) /= 0) IF (dbg) PRINT *,itt,' track:', itrack, 'assigned:', tracks(2,1:Ntprev,itt,it-1) ! Looking for coincidencies DO iip=1, Ntprev IF (tracks(2,iip,itt,it-1) == SpolysIndep(i,it)) THEN Indeppolychained = .TRUE. IF (dbg) PRINT *,' coincidence found by polygon:', tracks(2,iip,itt,it-1) EXIT END IF END DO IF (Indeppolychained) THEN Ntracks(it) = Ntracks(it) + 1 ictrack = Ntracks(it) ! Assigning all the IDs to the next step of the track DO iip=1, NpolysIndep(i,it) iiprev = polysIndep(i,iip,it) tracks(1,iip,ictrack,it) = itrack tracks(2,iip,ictrack,it) = iiprev ixp = ctrpolys(1,iiprev,it) iyp = ctrpolys(2,iiprev,it) tracks(3,iip,ictrack,it) = ixp tracks(4,iip,ictrack,it) = iyp tracks(5,iip,ictrack,it) = it END DO EXIT END IF END DO it1_tracks ! Creation of a new track IF (.NOT.Indeppolychained) THEN Ntracks(it) = Ntracks(it) + 1 ictrack = Ntracks(it) ! ID of new track maxtrack = INT(MAXVAL(tracks(1,:,:,:)*1.)) IF (dbg) PRINT *,' New track!', maxtrack+1 ! Assigning all the IDs to the next step of the track DO j=1, NpolysIndep(i,it) iiprev = polysIndep(i,j,it) tracks(1,j,ictrack,it) = maxtrack+1 tracks(2,j,ictrack,it) = iiprev ixp = ctrpolys(1,iiprev,it) iyp = ctrpolys(2,iiprev,it) tracks(3,j,ictrack,it) = ixp tracks(4,j,ictrack,it) = iyp tracks(5,j,ictrack,it) = it END DO END IF END DO current_poly IF (dbg) THEN PRINT *,' At this time-step:', it, ' N trajectories:', Ntracks(it) DO i=1, Ntracks(it) Nprev = COUNT(INT(tracks(2,:,i,it)) /= 0) PRINT *,i,'tracks:', tracks(2,1:Nprev,i,it) END DO END IF END DO timesteps ! Summarizing trajectories ! When multiple polygons are available, the mean of their central positions determines the position finaltracks = 0. maxtrack = MAXVAL(tracks(1,:,:,:)) DO it=1, dt DO itt=1, Ntracks(it) itrack = INT(tracks(1,1,itt,it)) Nprev = COUNT(INT(tracks(2,:,itt,it)) /= 0) PRINT *,'it:', it,'itrack:', itrack, 'Nprev:', Nprev finaltracks(1,itrack,it) = itrack*1. finaltracks(2,itrack,it) = SUM(tracks(3,:,itt,it))/Nprev finaltracks(3,itrack,it) = SUM(tracks(4,:,itt,it))/Nprev finaltracks(4,itrack,it) = it*1. PRINT *,' finaltrack:', finaltracks(:,itrack,it) END DO END DO RETURN END SUBROUTINE poly_overlap_tracks SUBROUTINE coincidence_all_polys(dbg, dx, dy, Nallpoly, allpoly, icpolys, Npoly, polys, cpolys, & apolys, polycoins, coinNptss) ! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons ! In case of multiple coincidencies, the closest and then the largest is taken IMPLICIT NONE LOGICAL, INTENT(in) :: dbg INTEGER, INTENT(in) :: dx, dy, Nallpoly, Npoly INTEGER, DIMENSION(dx,dy), INTENT(in) :: allpoly, polys REAL(r_k), DIMENSION(2,Nallpoly), INTENT(in) :: icpolys REAL(r_k), DIMENSION(2,Npoly), INTENT(in) :: cpolys REAL(r_k), DIMENSION(Npoly), INTENT(in) :: apolys INTEGER, DIMENSION(Nallpoly), INTENT(out) :: polycoins INTEGER, DIMENSION(Nallpoly,Npoly), INTENT(out) :: coinNptss ! Local INTEGER :: i, j, ip INTEGER :: maxcorr INTEGER :: Nmaxcorr LOGICAL, DIMENSION(dx,dy) :: boolpoly INTEGER :: maxcoin REAL :: dist, maxcoindist, maxcoinarea !!!!!!! Variables ! dx,dy: dimension of the space ! Nallpoly: Number of polygons to find coincidence ! allpoly: space with the polygons to meet ! icpolys: center of the polygons to look for the coincidence ! Npoly: number of polygons on the 2D space ! polys: 2D field of polygons identified by their integer number (0 for no polygon) ! cpolys: center of the polygons ! apolys: area of the polygons ! polycoins: coincident polyogn ! -1: no-coincidence ! 1 < Npoly: single coincidence with a given polygon ! -9: coincidence with more than one polygon ! coinNptss: number of points coincident with each polygon fname = 'coincidence_all_polys' IF (dbg) PRINT *,TRIM(fname) DO ip=1, Nallpoly boolpoly = allpoly == ip CALL coincidence_poly(dbg, dx, dy, boolpoly, Npoly, polys, polycoins(ip), coinNptss(ip,:)) IF (dbg) PRINT *,' polygon', ip, ' coincidence with:', polycoins(ip) ! Coincidence with more than one polygon IF (polycoins(ip) == -9) THEN maxcoindist = -10. maxcoinarea = -10. maxcoin = MAXVAL(coinNptss(ip,:)) DO j=1, Npoly IF (coinNptss(ip,j) == maxcoin) THEN dist = SQRT( (icpolys(1,ip)*1.-cpolys(1,j)*1.)**2 + (icpolys(2,ip)*1.-cpolys(2,j)*1.)**2 ) IF ( dist > maxcoindist) THEN maxcoindist = dist maxcoinarea = apolys(j) polycoins(ip) = j ELSE IF ( dist == maxcoindist) THEN IF (apolys(j) > maxcoinarea) THEN polycoins(ip) = j maxcoinarea = apolys(j) END IF END IF END IF END DO END IF END DO RETURN END SUBROUTINE coincidence_all_polys SUBROUTINE coincidence_poly(dbg, dx, dy, poly, Npoly, polys, polycoin, coinNpts) ! Subtourine to determine which is the coincident polygon when a boolean polygon is provided to a map of integer polygons IMPLICIT NONE LOGICAL, INTENT(in) :: dbg INTEGER, INTENT(in) :: dx, dy, Npoly LOGICAL, DIMENSION(dx,dy), INTENT(in) :: poly INTEGER, DIMENSION(dx,dy), INTENT(in) :: polys INTEGER, INTENT(out) :: polycoin INTEGER, DIMENSION(Npoly), INTENT(out) :: coinNpts ! Local INTEGER :: i, j, ip INTEGER :: maxcorr INTEGER :: Nmaxcorr !!!!!!! Variables ! dx,dy: dimension of the space ! poly: bolean polygon to meet ! Npoly: number of polygons on the 2D space ! polys: 2D field of polygons identified by their integer number (0 for no polygon) ! polycoin: coincident polyogn ! -1: no-coincidence ! 1 < Npoly: single coincidence with a given polygon ! -9: coincidence with more than one polygon ! coinNpts: number of points coincident with each polygon fname = 'coincidence_poly' IF (dbg) PRINT *,TRIM(fname) IF (dbg) THEN PRINT *,' Boolean polygon to search coincidences ...' DO i=1,dx PRINT *,poly(i,:) END DO PRINT *,' 2D polygons space ...' DO i=1,dx PRINT '(1000(I7,1x))',polys(i,:) END DO END IF ! Looking for coincient points for the polygon coinNpts = 0 DO i=1,dx DO j=1,dy IF (poly(i,j) .AND. polys(i,j) .NE. 0) coinNpts(polys(i,j)) = coinNpts(polys(i,j)) + 1 END DO END DO maxcorr = 0 maxcorr = INT(MAXVAL(coinNpts*1.)) IF (dbg) PRINT *,' Maximum coincidence:', maxcorr IF (maxcorr == 0) THEN polycoin = -1 ELSE Nmaxcorr = 0 DO ip=1, Npoly IF (coinNpts(ip) == maxcorr) THEN Nmaxcorr=Nmaxcorr+1 polycoin = ip END IF END DO IF (Nmaxcorr > 1) polycoin = -9 END IF IF (dbg) THEN PRINT *,' Coincidences for each polygon _______' DO ip=1, Npoly PRINT *,' ', ip,': ', coinNpts(ip) END DO END IF RETURN END SUBROUTINE coincidence_poly SUBROUTINE all_polygons_properties(dbg, dx, dy, Npoly, polys, lon, lat, values, xres, yres, projN, & Npolyptss, xxtrms, yxtrms, meanctrs, meanwctrs, areas, nvals, xvals, mvals, m2vals, stdvals, & Nquant, quants, nvcoords, xvcoords, meanvnctrs, meanvxctrs) ! Subroutine to determine the properties of all polygons in a 2D field: ! Number of grid points ! grid-point coordinates of the minimum and maximum of the path along x,y axes ! grid coordinates of center from the mean of the coordinates of the poygon locations ! lon, lat center from the area weighted mean of the coordinates of the polygon locations ! area of the polygon (km2) ! minimum and maximum of the values within the polygon ! mean of the values within the polygon ! quadratic mean of the values within the polygon ! standard deviation of the values within the polygon ! number of quantiles ! quantiles of the values within the polygon ! grid coordinates of the minimum, maximum value within the polygon ! lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon IMPLICIT NONE LOGICAL, INTENT(in) :: dbg INTEGER, INTENT(in) :: dx, dy, Npoly, Nquant INTEGER, DIMENSION(dx,dy), INTENT(in) :: polys REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: lon, lat, values REAL(r_k), INTENT(in) :: xres, yres CHARACTER(len=1000), INTENT(in) :: projN INTEGER, DIMENSION(Npoly), INTENT(out) :: Npolyptss INTEGER, DIMENSION(Npoly,2), INTENT(out) :: xxtrms, yxtrms, meanctrs REAL(r_k), DIMENSION(Npoly), INTENT(out) :: areas REAL(r_k), DIMENSION(Npoly), INTENT(out) :: nvals, xvals, mvals, m2vals, stdvals REAL(r_k), DIMENSION(Npoly, Nquant), INTENT(out) :: quants INTEGER, DIMENSION(Npoly,2), INTENT(out) :: nvcoords, xvcoords REAL(r_k), DIMENSION(Npoly,2), INTENT(out) :: meanwctrs, meanvnctrs, meanvxctrs ! Local INTEGER :: ip LOGICAL, DIMENSION(dx,dy) :: boolpoly !!!!!!! Variables ! dx,dy: size of the space ! Npoly: number of polygons ! polys: polygon matrix with all polygons (as integer number per polygon) ! lon, lat: geographical coordinates of the grid-points of the matrix ! values: values of the 2D field to use ! [x/y]res resolution along the x and y axis ! projN: name of the projection ! 'ctsarea': Constant Area ! 'lon/lat': for regular longitude-latitude ! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR) ! Npolyptss: number of points of the polygons ! [x/y]xtrms: grid-point coordinates of minimum and maximum coordinates of the polygons ! meanctrs: center from the mean of the coordinates of the polygons ! meanwctrs: lon, lat coordinates of the center from the spatial-weighted mean of the polygons ! areas: area of the polygons [km] ! [n/x]vals: minimum and maximum of the values within the polygons ! mvals: mean of the values within the polygons ! m2vals: quadratic mean of the values within the polygons ! stdvals: standard deviation of the values within the polygons ! Nquant: number of quantiles ! quants: quantiles of the values within the polygons ! [n/x]vcoords: grid coordinates of the minimum/maximum of the values within the polygons ! meanv[n/x]ctrs: lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygons fname = 'all_polygons_properties' ! Initializing matrices Npolyptss = -1 xxtrms = fillval64 yxtrms = fillval64 meanctrs = fillval64 meanwctrs = fillval64 areas = fillval64 nvals = fillvalI xvals = fillval64 mvals = fillval64 m2vals = fillval64 stdvals = fillval64 quants = fillval64 nvcoords = fillvalI xvcoords = fillvalI meanvnctrs = fillval64 meanvxctrs = fillval64 DO ip=1, Npoly boolpoly = polys == ip CALL polygon_properties(dbg, dx, dy, boolpoly, lon, lat, values, xres, yres, projN, Npolyptss(ip),& xxtrms(ip,:), yxtrms(ip,:), meanctrs(ip,:), meanwctrs(ip,:), areas(ip), nvals(ip), xvals(ip), & mvals(ip), m2vals(ip), stdvals(ip), Nquant, quants(ip,:), nvcoords(ip,:), xvcoords(ip,:), & meanvnctrs(ip,:), meanvxctrs(ip,:)) END DO RETURN END SUBROUTINE all_polygons_properties SUBROUTINE polygon_properties(dbg, dx, dy, poly, lon, lat, values, xres, yres, projN, Npolypts, & xxtrm, yxtrm, meanctr, meanwctr, area, nval, xval, mval, m2val, stdval, Nquant, quant, nvcoord, & xvcoord, meanvnctr, meanvxctr) ! Subroutine to determine the properties of a polygon (as .TRUE. matrix) ! Number of grid points ! grid-point coordinates of the minimum and maximum of the path along x,y axes ! grid coordinates of center from the mean of the coordinates of the poygon locations ! lon, lat center from the area weighted mean of the coordinates of the polygon locations ! area of the polygon (km2) ! minimum and maximum of the values within the polygon ! mean of the values within the polygon ! quadratic mean of the values within the polygon ! standard deviation of the values within the polygon ! number of quantiles ! quantiles of the values within the polygon ! grid coordinates of the minimum, maximum value within the polygon ! lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon IMPLICIT NONE LOGICAL, INTENT(in) :: dbg INTEGER, INTENT(in) :: dx, dy, Nquant LOGICAL, DIMENSION(dx,dy), INTENT(in) :: poly REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: lon, lat, values REAL(r_k), INTENT(in) :: xres, yres CHARACTER(len=1000), INTENT(in) :: projN INTEGER, INTENT(out) :: Npolypts INTEGER, DIMENSION(2), INTENT(out) :: xxtrm, yxtrm, meanctr INTEGER, DIMENSION(2), INTENT(out) :: nvcoord, xvcoord REAL(r_k), DIMENSION(2), INTENT(out) :: meanwctr, meanvnctr, meanvxctr REAL(r_k), INTENT(out) :: area REAL(r_k), INTENT(out) :: nval, xval, mval, m2val, stdval REAL(r_k), DIMENSION(Nquant), INTENT(out) :: quant ! Local INTEGER :: i, j, ip INTEGER :: ierr INTEGER, DIMENSION(:,:), ALLOCATABLE :: polypts REAL(r_k), DIMENSION(:), ALLOCATABLE :: polyvals, distvn, distvx REAL(r_k) :: lonNADIR, latNADIR REAL(r_k) :: sumRESx, sumRESy REAL(r_k), DIMENSION(dx,dy) :: xcorr, ycorr CHARACTER(len=200), DIMENSION(3) :: satSvals CHARACTER(len=50) :: projS REAL(r_k) :: sumDISTnlon, sumDISTnlat, sumDISTxlon, & sumDISTxlat !!!!!!! Variables ! dx,dy: size of the space ! poly: polygon matrix (boolean) ! lon, lat: geographical coordinates of the grid-points of the matrix ! values: values of the 2D field to use ! [x/y]res resolution along the x and y axis ! projN: name of the projection ! 'ctsarea': Constant Area ! 'lon/lat': for regular longitude-latitude ! 'nadir-sat,[lonNADIR],[latNADIR]': for satellite data with the resolution given at nadir (lonNADIR, latNADIR) ! Npolypts: number of points of the polygon ! [x/y]xtrm: grid-point coordinates of minimum and maximum coordinates of the polygon ! meanctr: center from the mean of the coordinates of the polygon ! meanwctr: lon, lat coordinates of the center from the spatial-weighted mean of the polygon ! area: area of the polygon [km] ! [n/x]val: minimum and maximum of the values within the polygon ! mval: mean of the values within the polygon ! m2val: quadratic mean of the values within the polygon ! stdval: standard deviation of the values within the polygon ! Nquant: number of quantiles ! quant: quantiles of the values within the polygon ! [n/x]vcoord: grid coordinates of the minimum/maximum of the values within the polygon ! meanv[n/x]ctr: lon, lat coordinates of the area center weighted and also by distance to the lowest or highest values of the polygon fname = 'polygon_properties' IF (dbg) PRINT *," '" // TRIM(fname) // "' ..." ! Getting grid-point coordinates of the polygon Npolypts = COUNT(poly) IF (ALLOCATED(polypts)) DEALLOCATE(polypts) ALLOCATE(polypts(Npolypts,2), STAT=ierr) msg = "Problems allocating 'polypts'" CALL ErrMsg(msg, fname, ierr) IF (ALLOCATED(polyvals)) DEALLOCATE(polyvals) ALLOCATE(polyvals(Npolypts), STAT=ierr) msg = "Problems allocating 'polyvals'" CALL ErrMsg(msg, fname, ierr) IF (ALLOCATED(distvn)) DEALLOCATE(distvn) ALLOCATE(distvn(Npolypts), STAT=ierr) msg = "Problems allocating 'distvn'" CALL ErrMsg(msg, fname, ierr) IF (ALLOCATED(distvx)) DEALLOCATE(distvx) ALLOCATE(distvx(Npolypts), STAT=ierr) msg = "Problems allocating 'distvx'" CALL ErrMsg(msg, fname, ierr) IF (projN(1:7) == 'lon/lat') THEN projS = projN ELSE IF (projN(1:7) == 'ctsarea') THEN projS = projN ELSE IF (projN(1:9) == 'nadir-sat') THEN projS = 'nadir-sat' CALL split(projN, ',', 3, satSvals) READ(satSvals(2),'(F200.100)')lonNadir READ(satSvals(3),'(F200.100)')latNadir IF (dbg) PRINT *," 'nadir-geostationary-satellite' based projection of data with nadir " // & "location at:", lonNadir, latNadir ELSE msg = "Projection '" // TRIM(projN) // "' not ready" // CHAR(10) // " available ones: " // & "'ctsarea', 'lon/lat', 'nadir-sat'" CALL ErrMsg(msg,fname,-1) END IF area = 0. sumRESx = 0. sumRESy = 0. meanwctr = 0. meanvnctr = 0. meanvxctr = 0. xcorr = 0. ycorr = 0. nval = fillval64 xval = -fillval64 ip = 1 DO i=1,dx DO j=1,dy IF (poly(i,j)) THEN polypts(ip,1) = i polypts(ip,2) = j polyvals(ip) = values(i,j) SELECT CASE (TRIM(projS)) CASE ('ctsarea') ! Constant Area xcorr(i,j) = 1. ycorr(i,j) = 1. CASE ('lon/lat') ! Area as fixed yres and sinus-lat varying for xres ! IF (KIND(xcorr(i,j)) == KIND(1.d0)) THEN ! xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad)) ! ELSE xcorr(i,j) = ABS(SIN(lon(i,j)*DegRad)) ! END IF ycorr(i,j) = 1. CASE ('nadir-sat') ! Area from nadir resolution and degrading as we get far from satellite's nadir ! GOES-E: 0 N, 75 W ! IF (KIND(xcorr(i,j)) == KIND(1.d0)) THEN ! xcorr(i,j) = DABS(DSIN(lon(i,j)*DegRad)) ! ELSE xcorr(i,j) = ABS(SIN(lon(i,j)*DegRad)) ! END IF ycorr(i,j) = 1. END SELECT area = area + xres*xcorr(i,j)*yres*ycorr(i,j) meanwctr(1) = meanwctr(1) + lon(i,j)*xres*xcorr(i,j) meanwctr(2) = meanwctr(2) + lat(i,j)*yres*ycorr(i,j) IF (nval > values(i,j)) THEN nvcoord(1) = i nvcoord(2) = j nval = values(i,j) END IF IF (xval < values(i,j)) THEN xvcoord(1) = i xvcoord(2) = j xval = values(i,j) END IF ip = ip + 1 END IF END DO END DO IF (dbg) THEN PRINT *,' grid_coord lon lat value _______ ' DO ip=1, Npolypts PRINT *, polypts(ip,:), ';', lon(polypts(ip,1),polypts(ip,2)), lat(polypts(ip,1),polypts(ip,2)),& ':', polyvals(ip) END DO END IF sumRESx = xres*SUM(xcorr) sumRESy = yres*SUM(ycorr) xxtrm = (/ MINVAL(polypts(:,1)), MAXVAL(polypts(:,1)) /) yxtrm = (/ MINVAL(polypts(:,2)), MAXVAL(polypts(:,2)) /) meanctr = (/ SUM(polypts(:,1))/Npolypts, SUM(polypts(:,2))/Npolypts /) meanwctr = (/ meanwctr(1)/sumRESx, meanwctr(2)/sumRESy /) IF (dbg) THEN PRINT *,' mean grid center: ', meanctr, ' weighted mean center: ', meanwctr END IF ! Statistics of the values within the polygon CALL StatsR_K(Npolypts, polyvals, nval, xval, mval, m2val, stdval) IF (dbg) THEN PRINT *,' minimum value: ', nval, ' maximum:', xval, ' mean:', mval PRINT *,' coor. minimum: ', nvcoord PRINT *,' coor. maximum: ', xvcoord END IF ! Mean center weighted to minimum and maximum values ! IF (KIND(polyvals(1)) == KIND(1.d0)) THEN ! distvn = DABS(polyvals - nval) ! distvx = DABS(polyvals - xval) ! ELSE distvn = ABS(polyvals - nval) distvx = ABS(polyvals - xval) ! END IF meanvnctr = 0. meanvxctr = 0. sumDISTnlon = 0. sumDISTnlat = 0. sumDISTxlon = 0. sumDISTxlat = 0. ip = 1 DO i=1,dx DO j=1,dy IF (poly(i,j)) THEN meanvnctr(1) = meanvnctr(1)+lon(i,j)*distvn(ip)*xres*xcorr(i,j) meanvnctr(2) = meanvnctr(2)+lat(i,j)*distvn(ip)*yres*ycorr(i,j) meanvxctr(1) = meanvxctr(1)+lon(i,j)*distvx(ip)*xres*xcorr(i,j) meanvxctr(2) = meanvxctr(2)+lat(i,j)*distvx(ip)*yres*ycorr(i,j) sumDISTnlon = sumDISTnlon + distvn(ip)*xres*xcorr(i,j) sumDISTnlat = sumDISTnlat + distvn(ip)*yres*ycorr(i,j) sumDISTxlon = sumDISTxlon + distvx(ip)*xres*xcorr(i,j) sumDISTxlat = sumDISTxlat + distvx(ip)*yres*ycorr(i,j) ip = ip + 1 END IF END DO END DO meanvnctr = (/ meanvnctr(1)/sumDISTnlon, meanvnctr(2)/sumDISTnlat /) meanvxctr = (/ meanvxctr(1)/sumDISTxlon, meanvxctr(2)/sumDISTxlat /) IF (dbg) THEN PRINT *,' mean center to minimum: ', meanvnctr, ' to maximum: ', meanvxctr END IF ! Quantiles of the values within the polygon quant = -9999.d0 IF (Npolypts > Nquant) THEN CALL quantilesR_K(Npolypts, polyvals, Nquant, quant) ELSE CALL SortR_K(polyvals, Npolypts) END IF DEALLOCATE (polypts) DEALLOCATE (polyvals) RETURN END SUBROUTINE polygon_properties SUBROUTINE polygons_t(dbg, dx, dy, dt, boolmatt, polys, Npoly) ! Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1! IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, dt LOGICAL, DIMENSION(dx,dy,dt), INTENT(in) :: boolmatt LOGICAL, INTENT(in) :: dbg INTEGER, DIMENSION(dt), INTENT(out) :: Npoly INTEGER, DIMENSION(dx,dy,dt), INTENT(out) :: polys ! Local INTEGER :: i,it !!!!!!! Variables ! dx,dy: spatial dimensions of the space ! boolmatt: boolean matrix tolook for the polygons (.TRUE. based) ! polys: found polygons ! Npoly: number of polygons found fname = 'polygons' IF (dbg) PRINT *,TRIM(fname) polys = -1 Npoly = 0 DO it=1,dt IF (ANY(boolmatt(:,:,it))) THEN IF (dbg) THEN PRINT *,' it:', it, ' num. TRUE:', COUNT(boolmatt(:,:,it)), 'bool _______' DO i=1,dx PRINT *,boolmatt(i,:,it) END DO END IF CALL polygons(dbg, dx, dy, boolmatt(:,:,it), polys(:,:,it), Npoly(it)) ELSE IF (dbg) THEN PRINT *,' it:', it, " without '.TRUE.' values skipiing it!!" END IF END IF END DO END SUBROUTINE polygons_t SUBROUTINE polygons(dbg, dx, dy, boolmat, polys, Npoly) ! Subroutine to search the polygons of a boolean field. FORTRAN based. 1st = 1! IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy LOGICAL, DIMENSION(dx,dy), INTENT(in) :: boolmat LOGICAL, INTENT(in) :: dbg INTEGER, INTENT(out) :: Npoly INTEGER, DIMENSION(dx,dy), INTENT(out) :: polys ! Local INTEGER :: i, j, ip, ipp, Nppt INTEGER :: ierr INTEGER, DIMENSION(:,:), ALLOCATABLE :: borders LOGICAL, DIMENSION(dx,dy) :: isborder, isbordery INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: paths INTEGER :: Npath INTEGER, DIMENSION(:), ALLOCATABLE :: Nptpaths INTEGER, DIMENSION(2) :: xtrx, xtry, meanpth INTEGER :: Nvertx, Npts INTEGER, DIMENSION(:,:), ALLOCATABLE :: vertxs, points LOGICAL, DIMENSION(:), ALLOCATABLE :: isin !!!!!!! Variables ! dx,dy: spatial dimensions of the space ! boolmat: boolean matrix tolook for the polygons (.TRUE. based) ! polys: found polygons ! Npoly: number of polygons found fname = 'polygons' polys = -1 ! The mathematical maximum woiuld be dx*dy/4, but let's be optimistic... (sorry Jero) Nppt = dx*dy/10 IF (ALLOCATED(borders)) DEALLOCATE(borders) ALLOCATE(borders(Nppt,2), STAT=ierr) msg = "Problems allocating matrix 'borders'" CALL ErrMsg(msg, fname, ierr) IF (ALLOCATED(paths)) DEALLOCATE(paths) ALLOCATE(paths(Nppt,Nppt,2), STAT=ierr) msg = "Problems allocating matrix 'paths'" CALL ErrMsg(msg, fname, ierr) IF (ALLOCATED(Nptpaths)) DEALLOCATE(Nptpaths) ALLOCATE(Nptpaths(Nppt), STAT=ierr) msg = "Problems allocating matrix 'Nptpaths'" CALL ErrMsg(msg, fname, ierr) ! Filling with the points of all the space with .TRUE. Npts = COUNT(boolmat) IF (ALLOCATED(points)) DEALLOCATE(points) ALLOCATE(points(Npts,2), STAT=ierr) msg = "Problems allocating matrix 'points'" CALL ErrMsg(msg, fname, ierr) ! We only want to localize that points 'inside' ip = 1 DO i=1, dx DO j=1, dy IF (boolmat(i,j)) THEN points(ip,1) = i points(ip,2) = j ip = ip + 1 END IF END DO END DO CALL borders_matrixL(dx, dy, Nppt, boolmat, borders, isborder, isbordery) CALL paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths) Npoly = Npath DO ip=1, Npath IF (ALLOCATED(vertxs)) DEALLOCATE(vertxs) ALLOCATE(vertxs(Nptpaths(ip),2)) msg = "Problems allocating matrix 'vertxs'" CALL ErrMsg(msg, fname, ierr) IF (ALLOCATED(isin)) DEALLOCATE(isin) ALLOCATE(isin(Npts), STAT=ierr) msg = "Problems allocating matrix 'isin'" CALL ErrMsg(msg, fname, ierr) isin = .FALSE. IF (dbg) PRINT *, ' path:', ip, ' N pts:', Nptpaths(ip) CALL path_properties(dx, dy, boolmat, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), xtrx, xtry, & meanpth, 'y', Nvertx, vertxs) IF (dbg) THEN PRINT *, ' properties _______' PRINT *, ' x-extremes:', xtrx PRINT *, ' y-extremes:', xtry PRINT *, ' center mean:', meanpth PRINT *, ' y-vertexs:', Nvertx,' ________' DO i=1, Nvertx PRINT *,' ',i,':',vertxs(i,:) END DO END IF CALL gridpoints_InsidePolygon(dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), Nvertx,& xtrx, xtry, vertxs, Npts, points, isin) ! Filling polygons DO ipp=1, Npts IF (isin(ipp)) polys(points(ipp,1),points(ipp,2)) = ip END DO IF (dbg) THEN PRINT *,' boolmat isborder isbordery polygon (',xtrx(1),',',xtry(1),')x(',xtrx(2),',',xtry(2), & ') _______' DO i=xtrx(1), xtrx(2) PRINT *,i,':',boolmat(i,xtry(1):xtry(2)), ' border ', isborder(i,xtry(1):xtry(2)), & ' isbordery ', isbordery(i,xtry(1):xtry(2)), ' polygon ', polys(i,xtry(1):xtry(2)) END DO END IF END DO ! Cleaning polygons matrix of not-used and paths around holes CALL clean_polygons(dx, dy, boolmat, polys, Npoly, dbg) DEALLOCATE (borders) DEALLOCATE (Nptpaths) DEALLOCATE (paths) DEALLOCATE (vertxs) DEALLOCATE (points) DEALLOCATE (isin) RETURN END SUBROUTINE polygons SUBROUTINE clean_polygons(dx, dy, Lmat, pols, Npols, dbg) ! Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat INTEGER, INTENT(inout) :: Npols INTEGER, DIMENSION(dx,dy), INTENT(inout) :: pols LOGICAL, INTENT(in) :: dbg ! Local INTEGER :: i,j,ip,iprm INTEGER, DIMENSION(Npols) :: origPol, NotPol, neigPol INTEGER :: ispol, NnotPol CHARACTER(len=4) :: ISa !!!!!!! Variables ! dx, dy: size of the space ! Lmat: original bolean matrix from which the polygons come from ! Npols: original number of polygons ! pols: polygons space fname = 'clean_polygons' IF (dbg) PRINT *," At '" // TRIM(fname) // "' ..." origPol = -1 ! Looking for polygons already in space NnotPol = 0 DO ip=1, Npols ispol = COUNT(pols-ip == 0) IF (ispol > 0) THEN origPol(ip) = ip ELSE NnotPol = NnotPol + 1 NotPol(NnotPol) = ip neigPol(NnotPol) = -1 END IF END DO IF (dbg) THEN PRINT *,' It should be:', Npols, ' polygons, but already there are:', Npols - NnotPol PRINT *,' Polygons to remove:', NotPol(1:NnotPol) END IF ! Looking for the hole border of a polygon. This is identify as such polygon point which along ! y-axis has NpolygonA, Npolygon, .FALSE. DO i=1,dx DO j=2,dy-1 IF ( (pols(i,j-1) /= pols(i,j) .AND. pols(i,j+1) == -1) .AND. (COUNT(NotPol-pols(i,j)==0)==0) & .AND. (pols(i,j) /= -1) .AND. (pols(i,j-1) /= -1)) THEN IF (dbg) PRINT *,' Polygon:', pols(i,j), ' to be removed at point (',i,',',j,'); j-1:', & pols(i,j-1), ' j:', pols(i,j), ' j+1:', pols(i,j+1) NnotPol = NnotPol + 1 NotPol(NnotPol) = pols(i,j) neigPol(NnotPol) = pols(i,j-1) END IF END DO END DO IF (dbg) THEN PRINT *,' It should be:', Npols, ' polygons, but already there are:', Npols - NnotPol PRINT *,' Polygons to remove after looking for fake border-of-hole polygons _______' DO i=1, NnotPol PRINT *, ' Polygon:', NotPol(i), ' to be replaced by:', neigPol(i) END DO END IF ! Removing polygons DO iprm=1, NnotPol IF (neigPol(iprm) == -1) THEN WHERE (pols == NotPol(iprm)) pols = -1 END WHERE IF (dbg) THEN PRINT *,' removing polygon:', NotPol(iprm) END IF ELSE WHERE (pols == NotPol(iprm)) pols = neigPol(iprm) END WHERE IF (dbg) THEN PRINT *,' replacing polygon:', NotPol(iprm), ' by:', neigPol(iprm) END IF END IF END DO ! Re-numbering (descending values) DO i = 1, NnotPol iprm = MAXVAL(NotPol(1:NnotPol)) WHERE(pols > iprm) pols = pols - 1 END WHERE j = Index1DArrayI(NotPol, NnotPol, iprm) NotPol(j) = -9 END DO Npols = Npols - NnotPol RETURN END SUBROUTINE clean_polygons SUBROUTINE path_properties(dx, dy, Lmat, Nptspth, pth, xxtrm, yxtrm, meanctr, axs, Nvrtx, vrtxs) ! Subroutine to determine the properties of a path: ! extremes: minimum and maximum of the path along x,y axes ! meancenter: center from the mean of the coordinates of the paths locations ! vertexs: path point, without neighbours along a given axis IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, Nptspth LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat INTEGER, DIMENSION(Nptspth,2), INTENT(in) :: pth CHARACTER, INTENT(in) :: axs INTEGER, DIMENSION(2), INTENT(out) :: meanctr, xxtrm, yxtrm INTEGER, INTENT(out) :: Nvrtx INTEGER, DIMENSION(Nptspth,2), INTENT(out) :: vrtxs ! Local INTEGER :: i, ip, jp INTEGER :: neig1, neig2 !!!!!!! Variables ! dx,dy: size of the space ! Lmat: original matrix of logical values for the path ! Nptspth: number of points of the path ! pth: path coordinates (clockwise) ! axs: axis of finding the vertex ! [x/y]xtrm: minimum and maximum coordinates of the path ! meanctr: center from the mean of the coordinates of the path ! Nvrtx: Number of vertexs of the path ! vrtxs: coordinates of the vertexs fname = 'path_properties' vrtxs = -1 Nvrtx = 0 xxtrm = (/ MINVAL(pth(:,1)), MAXVAL(pth(:,1)) /) yxtrm = (/ MINVAL(pth(:,2)), MAXVAL(pth(:,2)) /) meanctr = (/ SUM(pth(:,1))/Nptspth, SUM(pth(:,2))/Nptspth /) IF (axs == 'x' .OR. axs == 'X') THEN ! Looking vertexs along x-axis DO i=1, Nptspth ip = pth(i,1) jp = pth(i,2) neig1 = 0 neig2 = 0 ! W-point IF (ip == 1) THEN neig1 = -1 ELSE IF (.NOT.Lmat(ip-1,jp)) neig1 = -1 END IF ! E-point IF (ip == dx) THEN neig2 = -1 ELSE IF (.NOT.Lmat(ip+1,jp)) neig2 = -1 END IF IF (neig1 == -1 .AND. neig2 == -1) THEN Nvrtx = Nvrtx + 1 vrtxs(Nvrtx,:) = (/ip,jp/) END IF END DO ELSE IF (axs == 'y' .OR. axs == 'Y') THEN ! Looking vertexs along x-axis DO i=1, Nptspth ip = pth(i,1) jp = pth(i,2) neig1 = 0 neig2 = 0 ! S-point IF (jp == 1) THEN neig1 = -1 ELSE IF (.NOT.Lmat(ip,jp-1)) neig1 = -1 END IF ! N-point IF (jp == dy) THEN neig2 = -1 ELSE IF (.NOT.Lmat(ip,jp+1)) neig2 = -1 END IF IF (neig1 == -1 .AND. neig2 == -1) THEN Nvrtx = Nvrtx + 1 vrtxs(Nvrtx,:) = (/ ip, jp /) END IF END DO ELSE msg = "Axis '" // axs // "' not available" // CHAR(10) // " Available ones: 'x', 'X', 'y, 'Y'" CALL ErrMsg(msg, fname, -1) END IF RETURN END SUBROUTINE path_properties SUBROUTINE gridpoints_InsidePolygon(dx, dy, isbrdr, Npath, path, Nvrtx, xpathxtrm, ypathxtrm, & vrtxs, Npts, pts, inside) ! Subroutine to determine if a series of grid points are inside a polygon following ray casting algorithm ! FROM: https://en.wikipedia.org/wiki/Point_in_polygon IMPLICIT NONE INTEGER, INTENT(in) :: dx,dy,Npath,Nvrtx,Npts LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isbrdr INTEGER, DIMENSION(Npath,2), INTENT(in) :: path INTEGER, DIMENSION(2), INTENT(in) :: xpathxtrm, ypathxtrm INTEGER, DIMENSION(Npath,2) :: vrtxs INTEGER, DIMENSION(Npts,2), INTENT(in) :: pts LOGICAL, DIMENSION(Npts), INTENT(out) :: inside ! Local INTEGER :: i,j,ip,ix,iy INTEGER :: Nintersecs, isvertex, ispath INTEGER :: ierr LOGICAL, DIMENSION(:,:), ALLOCATABLE :: halo_brdr INTEGER :: Nbrbrdr !!!!!!! Variables ! dx,dy: space size ! Npath: number of points of the path of the polygon ! path: path of the polygon ! isbrdr: boolean matrix of the space wqith .T. on polygon border ! Nvrtx: number of vertexs of the path ! [x/y]pathxtrm extremes of the path ! vrtxs: vertexs of the path along y-axis ! Npts: number of points ! pts: points to look for ! inside: vector wether point is inside or not (coincident to a border is inside) fname = 'gridpoints_InsidePolygon' ! Creation of a 1-grid point larger matrix to deal with points reaching the limits IF (ALLOCATED(halo_brdr)) DEALLOCATE(halo_brdr) ALLOCATE(halo_brdr(dx+2,dy+2), STAT=ierr) msg = "Problems allocating matrix 'halo_brdr'" CALL ErrMsg(msg, fname, ierr) halo_brdr = .FALSE. DO i=1,dx halo_brdr(i+1,2:dy+1) = isbrdr(i,:) END DO inside = .FALSE. DO ip=1,Npts Nintersecs = 0 ix = pts(ip,1) iy = pts(ip,2) ! Point might be outside path range... IF (ix >= xpathxtrm(1) .AND. ix <= xpathxtrm(2) .AND. iy >= ypathxtrm(1) .AND. & iy <= ypathxtrm(2)) THEN ! It is a border point? ispath = index_list_coordsI(Npath, path, (/ix,iy/)) IF (isbrdr(ix,iy) .AND. (ispath /= -1)) THEN inside(ip) = .TRUE. CYCLE END IF ! Looking along y-axis ! Accounting for consecutives borders Nbrbrdr = 0 DO j=MAX(1,ypathxtrm(1)-1),iy-1 ! Only counting that borders that are not vertexs ispath = index_list_coordsI(Npath, path, (/ix,j/)) isvertex = index_list_coordsI(Npath, vrtxs, (/ix,j/)) IF (halo_brdr(ix+1,j+1) .AND. (ispath /= -1) .AND. (isvertex == -1) ) Nintersecs = Nintersecs + 1 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 Nbrbrdr = Nbrbrdr + 1 ELSE ! Will remove that consecutive borders above 2 IF (Nbrbrdr /= 0) THEN Nintersecs = Nintersecs - MAX(Nbrbrdr-1, 0) Nbrbrdr = 0 END IF END IF END DO IF (MOD(Nintersecs,2) /= 0) inside(ip) = .TRUE. END IF END DO RETURN END SUBROUTINE gridpoints_InsidePolygon SUBROUTINE look_clockwise_borders(dx,dy,Nbrdrs,brdrs,gbrdr,isbrdr,ix,iy,dbg,xf,yf,iff) ! Subroutine to look clock-wise for a next point within a collection of borders (limits of a region) IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, Nbrdrs, ix, iy INTEGER, DIMENSION(Nbrdrs,2), INTENT(in) :: brdrs LOGICAL, DIMENSION(Nbrdrs), INTENT(in) :: gbrdr LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isbrdr LOGICAL, INTENT(in) :: dbg INTEGER, INTENT(out) :: xf, yf, iff ! Local INTEGER :: isch CHARACTER(len=2), DIMENSION(8) :: Lclock INTEGER, DIMENSION(8,2) :: spt INTEGER :: iif, jjf !!!!!!! Variables ! dx, dy: 2D shape ot the space ! Nbrdrs: number of brdrs found in this 2D space ! brdrs: list of coordinates of the borders ! gbrdr: accounts for the use if the given border point ! isbrdr: accounts for the matrix of the point is a border or not ! ix,iy: coordinates of the point to start to find for ! xf,yf: coordinates of the found point ! iff: position of the border found within the list of borders fname = 'look_clockwise_borders' ! Looking clock-wise assuming that one starts from the westernmost point ! Label of the search lclock = (/ 'W ', 'NW', 'N ', 'NE', 'E ', 'SE', 'S ', 'SW' /) ! Transformation to apply !spt = (/ (/-1,0/), (/-1,1/), (/0,1/), (/1,1/), (/1,0/), (/1,-1/), (/0,-1/), (/-1,-1/) /) spt(:,1) = (/ -1, -1, 0, 1, 1, 1, 0, -1 /) spt(:,2) = (/ 0, 1, 1, 1, 0, -1, -1, -1 /) xf = -1 yf = -1 DO isch=1, 8 ! clock-wise search IF (spt(isch,1) >= 0) THEN iif = MIN(dx,ix+spt(isch,1)) ELSE iif = MAX(1,ix+spt(isch,1)) END IF IF (spt(isch,2) >= 0) THEN jjf = MIN(dy,iy+spt(isch,2)) ELSE jjf = MAX(1,iy+spt(isch,2)) END IF iff = index_list_coordsI(Nbrdrs, brdrs,(/iif,jjf/)) IF (iff > 0) THEN IF (dbg) PRINT *,' ' // lclock(isch) // '-point:', iif,jjf, ':', iff, 'is',isbrdr(iif,jjf), & 'got',gbrdr(iff) IF (isbrdr(iif,jjf) .AND. .NOT.gbrdr(iff)) THEN xf = iif yf = jjf EXIT END IF END IF END DO RETURN END SUBROUTINE look_clockwise_borders SUBROUTINE borders_matrixL(dx,dy,dxy,Lmat,brdrs,isbrdr,isbrdry) ! Subroutine to provide the borders of a logical array (interested in .TRUE.) IMPLICIT NONE INTEGER, INTENT(in) :: dx,dy,dxy LOGICAL, DIMENSION(dx,dy), INTENT(in) :: Lmat INTEGER, DIMENSION(dxy,2), INTENT(out) :: brdrs LOGICAL, DIMENSION(dx,dy), INTENT(out) :: isbrdr, isbrdry ! Local INTEGER :: i,j,ib !!!!!!! Variables ! dx,dy: size of the space ! dxy: maximum number of border points ! Lmat: Matrix to look for the borders ! brdrs: list of coordinates of the borders ! isbrdr: matrix with .T./.F. wether the given matrix point is a border or not ! isbrdry: matrix with .T./.F. wether the given matrix point is a border or not only along y-axis fname = 'borders_matrixL' isbrdr = .FALSE. brdrs = -1 ib = 1 ! Starting with the borders. If a given point is TRUE it is a path-vertex ! Along y-axis DO i=1, dx IF (Lmat(i,1) .AND. .NOT.isbrdr(i,1)) THEN brdrs(ib,1) = i brdrs(ib,2) = 1 isbrdr(i,1) = .TRUE. ib=ib+1 END IF IF (Lmat(i,dy) .AND. .NOT.isbrdr(i,dy)) THEN brdrs(ib,1) = i brdrs(ib,2) = dy isbrdr(i,dy) = .TRUE. ib=ib+1 END IF END DO ! Along x-axis DO j=1, dy IF (Lmat(1,j) .AND. .NOT.isbrdr(1,j)) THEN brdrs(ib,1) = 1 brdrs(ib,2) = j isbrdr(1,j) = .TRUE. ib=ib+1 END IF IF (Lmat(dx,j) .AND. .NOT.isbrdr(dx,j)) THEN brdrs(ib,1) = dx brdrs(ib,2) = j isbrdr(dx,j) = .TRUE. ib=ib+1 END IF END DO isbrdry = isbrdr ! Border as that when looking on x-axis points with Lmat(i) /= Lmat(i+1) DO i=1, dx-1 DO j=1, dy-1 IF ( Lmat(i,j) .NEQV. Lmat(i+1,j) ) THEN IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN brdrs(ib,1) = i brdrs(ib,2) = j isbrdr(i,j) = .TRUE. ib=ib+1 ELSE IF (Lmat(i+1,j) .AND. .NOT.isbrdr(i+1,j)) THEN brdrs(ib,1) = i+1 brdrs(ib,2) = j isbrdr(i+1,j) = .TRUE. ib=ib+1 END IF END IF ! y-axis IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN IF (Lmat(i,j) .AND. .NOT.isbrdr(i,j)) THEN brdrs(ib,1) = i brdrs(ib,2) = j isbrdr(i,j) = .TRUE. isbrdry(i,j) = .TRUE. ib=ib+1 ELSE IF (Lmat(i,j+1) .AND. .NOT.isbrdr(i,j+1)) THEN brdrs(ib,1) = i brdrs(ib,2) = j+1 isbrdr(i,j+1) = .TRUE. isbrdry(i,j+1) = .TRUE. ib=ib+1 END IF END IF END DO END DO DO i=1, dx-1 DO j=1, dy-1 ! y-axis IF ( Lmat(i,j) .NEQV. Lmat(i,j+1) ) THEN IF (Lmat(i,j)) THEN isbrdry(i,j) = .TRUE. ELSE IF (Lmat(i,j+1)) THEN isbrdry(i,j+1) = .TRUE. END IF END IF END DO END DO ! only y-axis adding bands of 2 grid points DO i=1, dx-1 DO j=2, dy-2 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 IF (Lmat(i,j)) THEN isbrdry(i,j) = .TRUE. isbrdry(i,j+1) = .TRUE. END IF END IF END DO END DO RETURN END SUBROUTINE borders_matrixL SUBROUTINE paths_border(dbg, dx, dy, isborder, Nppt, borders, paths, Npath, Nptpaths) ! Subroutine to search the paths of a border field. IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, Nppt LOGICAL, INTENT(in) :: dbg LOGICAL, DIMENSION(dx,dy), INTENT(in) :: isborder INTEGER, DIMENSION(Nppt,2), INTENT(in) :: borders INTEGER, DIMENSION(Nppt,Nppt,2), INTENT(out) :: paths INTEGER, INTENT(out) :: Npath INTEGER, DIMENSION(Nppt), INTENT(out) :: Nptpaths ! Local INTEGER :: i,j,k,ib INTEGER :: ierr INTEGER :: Nbrdr LOGICAL, DIMENSION(:), ALLOCATABLE :: gotbrdr, emptygotbrdr INTEGER :: iipth, ipath, ip, Nptspath INTEGER :: iib, jjb, iip, ijp, iif, jjf, iff LOGICAL :: found, finishedstep !!!!!!! Variables ! dx,dy: spatial dimensions of the space ! Nppt: possible number of paths and points that the paths can have ! isborder: boolean matrix which provide the borders of the polygon ! borders: coordinates of the borders of the polygon ! paths: coordinates of each found path ! Npath: number of paths found ! Nptpaths: number of points per path fname = 'paths_border' ! Sarting matrix paths = -1 Npath = 0 Nptspath = 0 Nptpaths = -1 ib=1 finishedstep = .FALSE. ! Number of border points DO ib=1, Nppt IF (borders(ib,1) == -1 ) EXIT END DO Nbrdr = ib-1 IF (dbg) THEN PRINT *,' borders _______' DO i=1,Nbrdr PRINT *,' ',i,':',borders(i,:) END DO END IF ! Matrix which keeps track if a border point has been located IF (ALLOCATED(gotbrdr)) DEALLOCATE(gotbrdr) ALLOCATE(gotbrdr(Nbrdr), STAT=ierr) msg = "Problems allocating matrix 'gotbrdr'" CALL ErrMsg(msg, fname, ierr) IF (ALLOCATED(emptygotbrdr)) DEALLOCATE(emptygotbrdr) ALLOCATE(emptygotbrdr(Nbrdr), STAT=ierr) msg = "Problems allocating matrix 'emptygotbrdr'" CALL ErrMsg(msg, fname, ierr) gotbrdr = .FALSE. emptygotbrdr = .FALSE. ! Starting the fun... ! Looking along the lines and when a border is found, starting from there in a clock-wise way iipth = 1 ipath = 1 DO ib=1,Nbrdr iib = borders(iipth,1) jjb = borders(iipth,2) ! Starting new path newpath: IF (.NOT.gotbrdr(iipth)) THEN ip = 1 Nptspath = 1 paths(ipath,ip,:) = borders(iipth,:) gotbrdr(iipth) = .TRUE. ! Looking for following clock-wise search ! Not looking for W, because search starts from the W iip = iib ijp = jjb DO k=1,Nbrdr IF (dbg) PRINT *,ipath,'iip jip:', iip, ijp found = .FALSE. CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE.,iif,jjf,iff) IF (iif /= -1) THEN ip=ip+1 paths(ipath,ip,:) = (/ iif,jjf /) found = .TRUE. gotbrdr(iff) = .TRUE. iip = iif ijp = jjf Nptspath = Nptspath + 1 END IF IF (dbg) THEN PRINT *,iib,jjb,' end of this round path:', ipath, '_____', gotbrdr DO i=1, Nptspath PRINT *,' ',i,':',paths(ipath,i,:) END DO END IF ! If it is not found a next point, might be because it is a non-polygon related value IF (.NOT.found) THEN IF (dbg) PRINT *,'NOT FOUND !!!', gotbrdr ! Are still there available borders? IF (ALL(gotbrdr) .EQV. .TRUE.) THEN finishedstep = .TRUE. Npath = ipath Nptpaths(ipath) = Nptspath EXIT ELSE Nptpaths(ipath) = Nptspath ! Let's have a look if the previous points in the path have already some 'non-located' neighbourgs DO i=Nptspath,1,-1 iip = paths(ipath,i,1) ijp = paths(ipath,i,2) CALL look_clockwise_borders(dx,dy,Nppt,borders,gotbrdr,isborder,iip,ijp,.FALSE., iif, & jjf,iff) IF (iif /= -1 .AND. iff /= -1) THEN IF (dbg) PRINT *,' re-take path from point:', iif,',',jjf,' n-path:', iff found = .TRUE. iipth = index_list_coordsI(Nppt, borders, (/iip,ijp/)) EXIT END IF END DO IF (.NOT.found) THEN ! Looking for the next available border point for the new path DO i=1,Nbrdr IF (.NOT.gotbrdr(i)) THEN iipth = i EXIT END IF END DO IF (dbg) PRINT *,' Looking for next path starting at:', iipth, ' point:', & borders(iipth,:) ipath=ipath+1 EXIT END IF END IF ELSE IF (dbg) PRINT *,' looking for next point...' END IF IF (finishedstep) EXIT END DO END IF newpath END DO Npath = ipath Nptpaths(ipath) = Nptspath DEALLOCATE (gotbrdr) DEALLOCATE (emptygotbrdr) RETURN END SUBROUTINE paths_border SUBROUTINE rand_sample(Nvals, Nsample, sample) ! Subroutine to randomly sample a range of indices IMPLICIT NONE INTEGER, INTENT(in) :: Nvals, Nsample INTEGER, DIMENSION(Nsample), INTENT(out) :: sample ! Local INTEGER :: i, ind, jmax REAL, DIMENSION(Nsample) :: randv CHARACTER(len=50) :: fname LOGICAL :: found LOGICAL, DIMENSION(Nvals) :: issampled CHARACTER(len=256) :: msg CHARACTER(len=10) :: IS1, IS2 !!!!!!! Variables ! Nvals: number of values ! Nsamples: number of samples ! sample: samnple fname = 'rand_sample' IF (Nsample > Nvals) THEN WRITE(IS1,'(I10)')Nvals WRITE(IS2,'(I10)')Nsample msg = 'Sampling of ' // TRIM(IS1) // ' is too big for ' // TRIM(IS1) // 'values' CALL ErrMsg(msg, fname, -1) END IF ! Generation of random numbers always the same series during the whole program! CALL RANDOM_NUMBER(randv) ! Making sure that we do not repeat any value issampled = .FALSE. DO i=1, Nsample ! Generation of the index from the random numbers ind = MAX(INT(randv(i)*Nvals), 1) IF (.NOT.issampled(ind)) THEN sample(i) = ind issampled(ind) = .TRUE. ELSE ! Looking around the given index !PRINT *,' Index :', ind, ' already sampled!', issampled(ind) found = .FALSE. DO jmax=1, Nvals ind = MIN(ind+jmax, Nvals) IF (.NOT.issampled(ind)) THEN sample(i) = ind issampled(ind) = .TRUE. found = .TRUE. EXIT END IF ind = MAX(1, ind-jmax) IF (.NOT.issampled(ind)) THEN sample(i) = ind issampled(ind) = .TRUE. found = .TRUE. EXIT END IF END DO IF (.NOT.found) THEN msg = 'sampling could not be finished due to absence of available value!!' CALL ErrMsg(msg, fname, -1) END IF END IF END DO RETURN END SUBROUTINE rand_sample SUBROUTINE PrintQuantilesR_K(Nvals, vals, Nquants, qtvs, bspc) ! Subroutine to print the quantiles of values REAL(r_k) IMPLICIT NONE INTEGER, INTENT(in) :: Nvals, Nquants REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals REAL(r_k), DIMENSION(Nquants), INTENT(in) :: qtvs CHARACTER(len=1000), OPTIONAL :: bspc ! Local INTEGER :: iq LOGICAL, DIMENSION(Nvals) :: search1, search2, search CHARACTER(len=6) :: RS1 CHARACTER(len=50) :: fname CHARACTER(len=1000) :: bspcS !!!!!!! Variables ! vals: series of values ! qtvs: values of the quantiles ! bspc: base quantity of spaces fname = 'PrintQuantilesR_K' IF (PRESENT(bspc)) THEN bspcS = bspc ELSE bspcS = ' ' END IF DO iq=1, Nquants-1 WHERE (vals >= qtvs(iq)) search1 = .TRUE. ELSEWHERE search1 = .FALSE. END WHERE WHERE (vals < qtvs(iq+1)) search2 = .TRUE. ELSEWHERE search2 = .FALSE. END WHERE WHERE (search1 .AND. search2) search = .TRUE. ELSEWHERE search = .FALSE. END WHERE WRITE(RS1, '(F6.2)')(iq)*100./(Nquants-1) PRINT *, TRIM(bspcS) // '[',iq,']', TRIM(RS1) // ' %:', qtvs(iq), 'N:', COUNT(search) END DO RETURN END SUBROUTINE PrintQuantilesR_K INTEGER FUNCTION FindMinimumR_K(x, dsize, Startv, Endv) ! Function returns the location of the minimum in the section between Start and End. IMPLICIT NONE INTEGER, INTENT(in) :: dsize REAL(r_k), DIMENSION(dsize), INTENT(in) :: x INTEGER, INTENT(in) :: Startv, Endv ! Local REAL(r_k) :: Minimum INTEGER :: Location INTEGER :: i Minimum = x(Startv) ! assume the first is the min Location = Startv ! record its position DO i = Startv+1, Endv ! start with next elements IF (x(i) < Minimum) THEN ! if x(i) less than the min? Minimum = x(i) ! Yes, a new minimum found Location = i ! record its position END IF END DO FindMinimumR_K = Location ! return the position END FUNCTION FindMinimumR_K SUBROUTINE SwapR_K(a, b) ! Subroutine swaps the values of its two formal arguments. IMPLICIT NONE REAL(r_k), INTENT(INOUT) :: a, b ! Local REAL(r_k) :: Temp Temp = a a = b b = Temp END SUBROUTINE SwapR_K SUBROUTINE SortR_K(x, Nx) ! Subroutine receives an array x() r_K and sorts it into ascending order. IMPLICIT NONE INTEGER, INTENT(IN) :: Nx REAL(r_k), DIMENSION(Nx), INTENT(INOUT) :: x ! Local INTEGER :: i INTEGER :: Location DO i = 1, Nx-1 ! except for the last Location = FindMinimumR_K(x, Nx-i+1, i, Nx) ! find min from this to last CALL SwapR_K(x(i), x(Location)) ! swap this and the minimum END DO END SUBROUTINE SortR_K SUBROUTINE quantilesR_K(Nvals, vals, Nquants, quants) ! Subroutine to provide the quantiles of a given set of values of type real 'r_k' IMPLICIT NONE INTEGER, INTENT(in) :: Nvals, Nquants REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals REAL(r_k), DIMENSION(Nquants), INTENT(out) :: quants ! Local INTEGER :: i REAL(r_k) :: minv, maxv REAL(r_k), DIMENSION(Nvals) :: sortedvals !!!!!!! Variables ! Nvals: number of values ! Rk: kind of real of the values ! vals: values ! Nquants: number of quants ! quants: values at which the quantile start minv = MINVAL(vals) maxv = MAXVAL(vals) sortedvals = vals ! Using from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90 CALL SortR_K(sortedvals, Nvals) quants(1) = minv DO i=2, Nquants quants(i) = sortedvals(INT((i-1)*Nvals/Nquants)) END DO END SUBROUTINE quantilesR_K SUBROUTINE StatsR_K(Nvals, vals, minv, maxv, mean, mean2, stdev) ! Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a ! series of r_k numbers IMPLICIT NONE INTEGER, INTENT(in) :: Nvals REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals REAL(r_k), INTENT(out) :: minv, maxv, mean, mean2, stdev !!!!!!! Variables ! Nvals: number of values ! vals: values ! minv: minimum value of values ! maxv: maximum value of values ! mean: mean value of values ! mean2: quadratic mean value of values ! stdev: standard deviation of values minv = MINVAL(vals) maxv = MAXVAL(vals) mean=SUM(vals) mean2=SUM(vals*vals) mean=mean/Nvals mean2=mean2/Nvals stdev=SQRT(mean2 - mean*mean) RETURN END SUBROUTINE StatsR_k SUBROUTINE NcountR(values, d1, Ndiffvals, counts) ! Subroutine to count real values IMPLICIT NONE INTEGER, INTENT(in) :: d1 REAL(r_k), DIMENSION(d1), INTENT(in) :: values INTEGER, INTENT(out) :: Ndiffvals REAL(r_k), DIMENSION(d1,2), INTENT(out) :: counts ! Local INTEGER :: i, ival REAL(r_k), DIMENSION(d1) :: diffv !!!!!!! Variables ! values: values to count ! counts: counts of time for each value fname = 'NcountR' counts = -1. counts(1,1) = values(1) counts(1,2) = 1 Ndiffvals = 1 DO i=2,d1 diffv(1:Ndiffvals) = counts(1:Ndiffvals,1) - values(i) IF (ANY(diffv(1:Ndiffvals) == 0)) THEN ival = Index1DArrayR(counts(1:Ndiffvals,1), Ndiffvals, values(i)) counts(ival,2) = counts(ival,2) + 1 ELSE Ndiffvals = Ndiffvals + 1 counts(Ndiffvals,1) = values(i) counts(Ndiffvals,2) = 1 END IF END DO END SUBROUTINE NcountR SUBROUTINE runmean_F1D(d1, values, Nmean, headertail, runmean) ! Subroutine fo computing the running mean of a given set of float 1D values IMPLICIT NONE INTEGER, INTENT(in) :: d1, Nmean REAL(r_k), DIMENSION(d1), INTENT(in) :: values CHARACTER(len=*), INTENT(in) :: headertail REAL(r_k), DIMENSION(d1), INTENT(out) :: runmean ! Local INTEGER :: i, j, Nmean2 CHARACTER(len=5) :: NmeanS !!!!!!! Variables ! values: values to compute the running mean ! Nmean: number of odd points to use for the running mean ! headertail: How to proceed for the grid points at the beginning of the values which are not ! encompassed by the Nmean ! 'miss': set as missing values (1.d20) ! 'original': keep the original values ! 'progressfill': mean the values as a progressive running filter (e.g. for Nmean=5): ! runmean[values(1)] = values(1) ! runmean[values(2)] = MEAN(values(1:3)) ! runmean[values(3)] = MEAN(values(1:5)) ! runmean[values(4)] = MEAN(values(2:6)) ! (...) ! runmean[values(d1-2)] = MEAN(values(d1-5:d1)) ! runmean[values(d1-1)] = MEAN(values(d1-2:d1)) ! runmean[values(d1)] = MEAN(values(dd1)) ! 'zero': set as zero values ! runmean: runnig mean values fname = 'runmean_F1D' IF (MOD(Nmean,2) == 0) THEN WRITE(NmeanS,'(I5)')Nmean msg="Nmean has to be odd!! value provided: "// NmeanS CALL ErrMsg(msg, fname, -1) END IF Nmean2 = Nmean/2 SELECT CASE (TRIM(headertail)) CASE ('missing') runmean = fillval64 CASE ('original') runmean = values CASE ('progressfill') DO i=1, Nmean2 runmean(i) = SUM(values(1:2*(i-1)+1))/(2*(i-1)+1) END DO runmean(d1) = values(d1) DO i=2, Nmean2 j = d1-(2*(i-1)) runmean(d1-(i-1)) = SUM(values(j:d1))/(2*(i-1)+1) END DO CASE ('zero') runmean = zeroRK CASE DEFAULT msg = "'" // TRIM(headertail) // "' not available !!" //CHAR(44) // " available ones: " // & "'missing', 'original', 'progressfill', 'zero'" CALL ErrMsg(msg, fname, -1) END SELECT DO i= 1+Nmean2, d1 - Nmean2 runmean(i) = SUM(values(i-Nmean2:i+Nmean2))/Nmean END DO END SUBROUTINE runmean_F1D SUBROUTINE percentiles_R_K2D(values, axisS, Npercen, d1, d2, percentiles) ! Subroutine to compute the percentiles of a 2D R_K array along given set of axis IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, Npercen REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: values CHARACTER(LEN=*), INTENT(in) :: axisS REAL(r_k), DIMENSION(d1, d2, Npercen), INTENT(out) :: percentiles ! Local INTEGER :: i INTEGER :: Lstring, LaxisS, iichar CHARACTER(LEN=1000) :: splitaxis INTEGER, DIMENSION(1) :: axis1 CHARACTER(LEN=200), DIMENSION(2) :: axis2S INTEGER, DIMENSION(2) :: axis2 CHARACTER(LEN=1) :: Naxs !!!!!!! Variables ! d1,d2: length of the 2D dimensions ! values: values to use to compute the percentiles ! axisS: ':' separated list of axis to use to compute the percentiles ('all' for all axes) ! Npercen: number of percentiles ! percentiles: percentiles of the daata fname = 'percentiles_R_K2D' LaxisS = LEN_TRIM(axisS) iichar = numberTimes(axisS(1:LaxisS), ':') splitaxis = '' splitaxis(1:LaxisS) = axisS(1:LaxisS) percentiles = 0. IF (iichar == 0) THEN READ(axisS,'(I1)')axis1(1) ELSE IF (iichar == 1) THEN CALL split(splitaxis, ':', 2, axis2S) ELSE WRITE(Naxs,'(A1)')iichar msg = "' rank 2 values can not compute percentiles using " // Naxs // "' number of axis !!" CALL ErrMsg(msg, fname, -1) END IF IF (TRIM(axisS) == 'all') iichar = 2 IF (iichar == 0) THEN ! Might be a better way, but today I can't think it !! IF (axis1(1) == 1) THEN DO i=1, d2 CALL quantilesR_K(d1, values(:,i), Npercen, percentiles(1,i,:)) END DO ELSE IF (axis1(1) == 2) THEN DO i=1, d1 CALL quantilesR_K(d2, values(i,:), Npercen, percentiles(i,1,:)) END DO ELSE WRITE(Naxs,'(A1)')axis1(1) msg = "' rank 2 values can not compute percentiles using axis " // Naxs // "' !!" CALL ErrMsg(msg, fname, -1) END IF ELSE CALL quantilesR_K(d1*d2, RESHAPE(values, (/d1*d2/)), Npercen, percentiles(1,1,:)) END IF END SUBROUTINE percentiles_R_K2D SUBROUTINE percentiles_R_K3D(values, axisS, Npercen, d1, d2, d3, percentiles) ! Subroutine to compute the percentiles of a 3D R_K array along given set of axis IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, Npercen REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: values CHARACTER(LEN=*), INTENT(in) :: axisS REAL(r_k), DIMENSION(d1, d2, d3, Npercen), INTENT(out) :: percentiles ! Local INTEGER :: i, j INTEGER :: Lstring, LaxisS, iichar CHARACTER(LEN=1000) :: splitaxis INTEGER, DIMENSION(1) :: axis1 CHARACTER(LEN=200), DIMENSION(2) :: axis2S INTEGER, DIMENSION(2) :: axis2 CHARACTER(LEN=200), DIMENSION(3) :: axis3S INTEGER, DIMENSION(3) :: axis3 CHARACTER(LEN=1) :: Naxs1, Naxs2 !!!!!!! Variables ! d1,d2: length of the 2D dimensions ! values: values to use to compute the percentiles ! axisS: ':' separated list of axis to use to compute the percentiles ('all' for all axes) ! Npercen: number of percentiles ! percentiles: percentiles of the daata fname = 'percentiles_R_K3D' LaxisS = LEN_TRIM(axisS) iichar = numberTimes(axisS(1:LaxisS), ':') splitaxis = '' splitaxis(1:LaxisS) = axisS(1:LaxisS) percentiles = 0. IF (iichar == 0) THEN READ(axisS,'(I1)')axis1(1) ELSE IF (iichar == 1) THEN CALL split(splitaxis, ':', 2, axis2S) DO i=1,2 READ(axis2S(i), '(I1)')axis2(i) END DO ELSE IF (iichar == 2) THEN CALL split(splitaxis, ':', 3, axis3S) ELSE READ(Naxs1,'(A1)')iichar msg = "' rank 3 values can not compute percentiles using " // Naxs1 // "' number of axis !!" CALL ErrMsg(msg, fname, -1) END IF IF (TRIM(axisS) == 'all') iichar = 3 IF (iichar == 0) THEN ! Might be a better way, but today I can't think it !! IF (axis1(1) == 1) THEN DO i=1, d2 DO j=1, d3 CALL quantilesR_K(d1, values(:,i,j), Npercen, percentiles(1,i,j,:)) END DO END DO ELSE IF (axis1(1) == 2) THEN DO i=1, d1 DO j=1, d3 CALL quantilesR_K(d2, values(i,:,j), Npercen, percentiles(i,1,j,:)) END DO END DO ELSE IF (axis1(1) == 3) THEN DO i=1, d1 DO j=1, d2 CALL quantilesR_K(d3, values(i,j,:), Npercen, percentiles(i,j,1,:)) END DO END DO ELSE WRITE(Naxs1,'(A1)')axis1(1) msg = "' rank 3 values can not compute percentiles using axis " // Naxs1 // "' !!" CALL ErrMsg(msg, fname, -1) END IF ELSE IF (iichar == 1) THEN ! Might be a better way, but today I can't think it !! IF (axis2(1) == 1 .AND. axis2(2) == 2) THEN DO i=1, d3 CALL quantilesR_K(d1*d2, RESHAPE(values(:,:,i), (/d1*d2/)), Npercen, percentiles(1,1,i,:)) END DO ELSE IF (axis2(1) == 1 .AND. axis2(2) == 3) THEN DO i=1, d2 CALL quantilesR_K(d1*d3, RESHAPE(values(:,i,:), (/d1*d3/)), Npercen, percentiles(1,i,1,:)) END DO ELSE IF (axis2(1) == 2 .AND. axis2(2) == 3) THEN DO i=1, d1 CALL quantilesR_K(d2*d3, RESHAPE(values(i,:,:), (/d2*d3/)), Npercen, percentiles(i,1,1,:)) END DO ELSE WRITE(Naxs1,'(A1)')axis2(1) WRITE(Naxs2,'(A1)')axis2(2) msg="' rank 3 values can not compute percentiles using axis "//Naxs1// ', ' // Naxs2 // "' !!" CALL ErrMsg(msg, fname, -1) END IF ELSE CALL quantilesR_K(d1*d2*d3, RESHAPE(values, (/d1*d2*d3/)), Npercen, percentiles(1,1,1,:)) END IF END SUBROUTINE percentiles_R_K3D SUBROUTINE percentiles_R_K4D(values, axisS, Npercen, d1, d2, d3, d4, percentiles) ! Subroutine to compute the percentiles of a 4D R_K array along given set of axis IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d4, Npercen REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: values CHARACTER(LEN=*), INTENT(in) :: axisS REAL(r_k), DIMENSION(d1,d2,d3,d4,Npercen), INTENT(out) :: percentiles ! Local INTEGER :: i, j, k INTEGER :: Lstring, LaxisS, iichar CHARACTER(LEN=1000) :: splitaxis INTEGER, DIMENSION(1) :: axis1 CHARACTER(LEN=200), DIMENSION(2) :: axis2S INTEGER, DIMENSION(2) :: axis2 CHARACTER(LEN=200), DIMENSION(3) :: axis3S INTEGER, DIMENSION(3) :: axis3 CHARACTER(LEN=200), DIMENSION(4) :: axis4S CHARACTER(LEN=1) :: Naxs1, Naxs2, Naxs3 !!!!!!! Variables ! d1,d2: length of the 2D dimensions ! values: values to use to compute the percentiles ! axisS: ':' separated list of axis to use to compute the percentiles ('all' for all axes) ! Npercen: number of percentiles ! percentiles: percentiles of the daata fname = 'percentiles_R_K3D' LaxisS = LEN_TRIM(axisS) iichar = numberTimes(axisS(1:LaxisS), ':') splitaxis = '' splitaxis(1:LaxisS) = axisS(1:LaxisS) percentiles = 0. PRINT *,'iichar:', iichar, axisS(1:LaxisS) IF (iichar == 0) THEN READ(axisS,'(I1)')axis1(1) ELSE IF (iichar == 1) THEN CALL split(splitaxis, ':', 2, axis2S) DO i=1,2 READ(axis2S(i), '(I1)')axis2(i) END DO ELSE IF (iichar == 2) THEN CALL split(splitaxis, ':', 3, axis3S) DO i=1,3 READ(axis3S(i), '(I1)')axis3(i) END DO ELSE IF (iichar == 3) THEN CALL split(splitaxis, ':', 4, axis4S) ELSE READ(Naxs1,'(A1)')iichar msg = "' rank 4 values can not compute percentiles using " // Naxs1 // "' number of axis !!" CALL ErrMsg(msg, fname, -1) END IF IF (TRIM(axisS) == 'all') iichar = 4 IF (iichar == 0) THEN ! Might be a better way, but today I can't think it !! IF (axis1(1) == 1) THEN DO i=1, d2 DO j=1, d3 DO k=1, d4 CALL quantilesR_K(d1, values(:,i,j,k), Npercen, percentiles(1,i,j,k,:)) END DO END DO END DO ELSE IF (axis1(1) == 2) THEN DO i=1, d1 DO j=1, d3 DO k=1, d4 CALL quantilesR_K(d2, values(i,:,j,k), Npercen, percentiles(i,1,j,k,:)) END DO END DO END DO ELSE IF (axis1(1) == 3) THEN DO i=1, d1 DO j=1, d2 DO k=1, d4 CALL quantilesR_K(d3, values(i,j,:,k), Npercen, percentiles(i,j,1,k,:)) END DO END DO END DO ELSE IF (axis1(1) == 4) THEN DO i=1, d1 DO j=1, d2 DO k=1, d3 CALL quantilesR_K(d4, values(i,j,k,:), Npercen, percentiles(i,j,k,1,:)) END DO END DO END DO ELSE WRITE(Naxs1,'(A1)')axis1(1) msg = "' rank 3 values can not compute percentiles using axis " // Naxs1 // "' !!" CALL ErrMsg(msg, fname, -1) END IF ELSE IF (iichar == 1) THEN ! Might be a better way, but today I can't think it !! IF (axis2(1) == 1 .AND. axis2(2) == 2) THEN DO i=1, d3 DO j=1, d4 CALL quantilesR_K(d1*d2, RESHAPE(values(:,:,i,j), (/d1*d2/)), Npercen, & percentiles(1,1,i,j,:)) END DO END DO ELSE IF (axis2(1) == 1 .AND. axis2(2) == 3) THEN DO i=1, d2 DO j=1, d4 CALL quantilesR_K(d1*d3, RESHAPE(values(:,i,:,j), (/d1*d3/)), Npercen, & percentiles(1,i,1,j,:)) END DO END DO ELSE IF (axis2(1) == 1 .AND. axis2(2) == 4) THEN DO i=1, d2 DO j=1, d3 CALL quantilesR_K(d1*d4, RESHAPE(values(:,i,j,:), (/d1*d4/)), Npercen, & percentiles(1,i,j,1,:)) END DO END DO ELSE IF (axis2(1) == 2 .AND. axis2(2) == 3) THEN DO i=1, d1 DO j=1, d4 CALL quantilesR_K(d2*d3, RESHAPE(values(i,:,:,j), (/d2*d3/)), Npercen, & percentiles(i,1,1,j,:)) END DO END DO ELSE IF (axis2(1) == 2 .AND. axis2(2) == 4) THEN DO i=1, d1 DO j=1, d3 CALL quantilesR_K(d2*d4, RESHAPE(values(i,:,j,:), (/d2*d4/)), Npercen, & percentiles(i,1,j,1,:)) END DO END DO ELSE IF (axis2(1) == 3 .AND. axis2(2) == 4) THEN DO i=1, d1 DO j=1, d2 CALL quantilesR_K(d3*d4, RESHAPE(values(i,j,:,:), (/d3*d4/)), Npercen, & percentiles(i,j,1,1,:)) END DO END DO ELSE WRITE(Naxs1,'(A1)')axis2(1) WRITE(Naxs2,'(A1)')axis2(2) msg="' rank 4 values can not compute percentiles using axis "//Naxs1// ', ' // Naxs2 // "' !!" CALL ErrMsg(msg, fname, -1) END IF ELSE IF (iichar == 2) THEN IF (axis2(1) == 1 .AND. axis2(2) == 2 .AND. axis3(3) == 3) THEN DO i=1, d4 CALL quantilesR_K(d1*d2*d3, RESHAPE(values(:,:,:,i), (/d1*d2*d3/)), Npercen, & percentiles(1,1,1,i,:)) END DO ELSE IF (axis2(1) == 1 .AND. axis2(2) == 2 .AND. axis3(3) == 4) THEN DO i=1, d3 CALL quantilesR_K(d1*d2*d4, RESHAPE(values(:,:,i,:), (/d1*d2*d4/)), Npercen, & percentiles(1,1,i,1,:)) END DO ELSE IF (axis2(1) == 1 .AND. axis2(2) == 3 .AND. axis3(3) == 4) THEN DO i=1, d2 CALL quantilesR_K(d1*d3*d4, RESHAPE(values(:,i,:,:), (/d1*d3*d4/)), Npercen, & percentiles(1,i,1,1,:)) END DO ELSE IF (axis2(1) == 2 .AND. axis2(2) == 3 .AND. axis3(3) == 4) THEN DO i=1, d1 CALL quantilesR_K(d2*d3*d4, RESHAPE(values(i,:,:,:), (/d2*d3*d4/)), Npercen, & percentiles(i,1,1,1,:)) END DO ELSE WRITE(Naxs1,'(A1)')axis3(1) WRITE(Naxs2,'(A1)')axis3(2) WRITE(Naxs3,'(A1)')axis3(2) msg="' rank 4 values can not compute percentiles using axis "// Naxs1 // ', ' // Naxs2 // & ', ' // Naxs3 //"' !!" CALL ErrMsg(msg, fname, -1) END IF ELSE CALL quantilesR_K(d1*d2*d3*d4, RESHAPE(values, (/d1*d2*d3*d4/)), Npercen, percentiles(1,1,1,1,:)) END IF END SUBROUTINE percentiles_R_K4D REAL(r_k) FUNCTION distanceRK(pointA, pointB) ! Function to provide the distance between two points IMPLICIT NONE REAL(r_k), DIMENSION(2), INTENT(in) :: pointA, pointB !!!!!!! Variables ! pointA, B: couple of points to compute the distance between them fname = 'distanceRK' distanceRK = SQRT( (pointB(1)-pointA(1))**2 + (pointB(2)-pointA(2))**2 ) END FUNCTION distanceRK REAL(r_k) FUNCTION shoelace_area_polygon(Nvertex, poly) ! Computing the area of a polygon using sholace formula ! FROM: https://en.wikipedia.org/wiki/Shoelace_formula IMPLICIT NONE INTEGER, INTENT(in) :: Nvertex REAL(r_k), DIMENSION(Nvertex,2), INTENT(in) :: poly ! Local INTEGER :: i REAL(r_k) :: areapos, areaneg !!!!!!! Variables ! Nvertex: number of vertices of the polygon ! poly: coordinates of the vertex of the polygon (sorted) fname = 'shoelace_area_polygon' areapos = 0. areaneg = 0. DO i=1, Nvertex-1 areapos = areapos + poly(i,1)*poly(i+1,2) areaneg = areaneg + poly(i+1,1)*poly(i,2) END DO areapos = areapos + poly(Nvertex,1)*poly(1,2) areaneg = areaneg + poly(1,1)*poly(Nvertex,2) shoelace_area_polygon = 0.5*(areapos - areaneg) END FUNCTION shoelace_area_polygon SUBROUTINE intersection_2Dlines(lineA, lineB, intersect, ptintersect) ! Subroutine to provide the intersection point between two lines on the plane using Cramer's method IMPLICIT NONE REAL(r_k), DIMENSION(2,2), INTENT(in) :: lineA, lineB LOGICAL, INTENT(out) :: intersect REAL(r_k), DIMENSION(2), INTENT(out) :: ptintersect ! Local REAL(r_k), DIMENSION(2) :: segmentA, segmentB REAL(r_k) :: a11, a12, a21, a22, z1, z2 REAL(r_k) :: det, detX, detY LOGICAL :: axisAx, axisBx, axisAy, axisBy !!!!!!! Variables ! lineA: couple of coordinates for the line A ! lineB: couple of coordinates for the line B ! intersect: whether two lines intersect ! ptintersect: point of intersection [(0,0) if they do not intersect] fname = 'intersection_2Dlines' axisAx = .FALSE. axisAy = .FALSE. axisBx = .FALSE. axisBy = .FALSE. ! Setting segment parameters y = A + B*x IF (lineA(2,1) /= lineA(1,1)) THEN segmentA(2) = (lineA(2,2)-lineA(1,2))/(lineA(2,1)-lineA(1,1)) IF ( (lineA(1,1)*segmentA(2) - lineA(1,2)) /= (lineA(2,1)*segmentA(2) - lineA(2,2)) ) THEN PRINT *,'A = y1 - x2*B = ', lineA(1,2) - lineA(1,1)*segmentA(2) PRINT *,'A = y2 - x2*B = ', lineA(2,2) - lineA(2,1)*segmentA(2) msg = 'Wrong calculation of parameter A, for lineA' CALL ErrMSg(msg, fname, -1) END IF segmentA(1) = lineA(1,2) - lineA(1,1)*segmentA(2) a11 = segmentA(2) a12 = -oneRK z1 = -segmentA(1) IF (lineA(2,2) == lineA(1,2)) axisAx = .TRUE. ELSE ! lineA || y-axis axisAy = .TRUE. END IF IF (lineB(2,1) /= lineB(1,1)) THEN segmentB(2) = (lineB(2,2)-lineB(1,2))/(lineB(2,1)-lineB(1,1)) IF ( (lineB(1,1)*segmentB(2) - lineB(1,2)) /= (lineB(2,1)*segmentB(2) - lineB(2,2)) ) THEN PRINT *,'A = x1*B -y1 = ', lineB(1,1)*segmentB(2) - lineB(1,2) PRINT *,'A = x2*B -y2 = ', lineB(2,1)*segmentB(2) - lineB(2,2) msg = 'Wrong calculation of parameter A, for lineB' CALL ErrMSg(msg, fname, -1) END IF segmentB(1) = lineB(1,2) - lineB(1,1)*segmentB(2) a21 = segmentB(2) a22 = -oneRK z2 = -segmentB(1) IF (lineB(2,2) == lineB(1,2)) axisBx = .TRUE. ELSE ! lineB || y-axis axisBy = .TRUE. END IF ! Cramer's method ! a11 = B1; a12 = -1 ! a21 = B2; a22 = -1 ! z1 = -A1 ! z2 = -A2 ! (a11 a12)(x) (z1) ! (a21 a22)(y) (z2) ! -------- ------ ----- ---- --- -- - ! det = (a11*a22-a12*a21) ! detX = (z1*a22-z2*a21) ! detY = (a11*z1-a12*z2) ! ptintercept = (detX/det, detY/det) ! Cases when some of the lines are parallel to any given axis ! PRINT *,' axisAx', axisAx, 'axisAy', axisAy, 'axisBx', axisBx, 'axisBy', axisBy IF (axisAx .OR. axisAy .OR. axisBx .OR. axisBy) THEN IF (axisAx) THEN IF (axisBy) THEN intersect = .TRUE. ptintersect(1) = lineB(1,1) ptintersect(2) = lineA(1,2) ELSE intersect = .TRUE. ptintersect(1) = (lineA(1,2)-segmentB(1))/segmentB(2) ptintersect(2) = lineA(1,2) END IF END IF IF (axisAy) THEN IF (axisBy) THEN intersect = .TRUE. ptintersect(1) = lineA(1,1) ptintersect(2) = lineB(1,2) ELSE intersect = .TRUE. ptintersect(1) = lineA(1,1) ptintersect(2) = segmentB(1) + lineA(1,1)*segmentB(2) END IF END IF IF (axisBx) THEN IF (axisAy) THEN intersect = .TRUE. ptintersect(1) = lineA(1,1) ptintersect(2) = lineB(1,2) ELSE intersect = .TRUE. ptintersect(1) = (lineB(1,2)-segmentA(1))/segmentA(2) ptintersect(2) = lineB(1,2) END IF END IF IF (axisBy) THEN IF (axisAx) THEN intersect = .TRUE. ptintersect(1) = lineB(1,1) ptintersect(2) = lineA(1,2) ELSE intersect = .TRUE. ptintersect(1) = lineB(1,1) ptintersect(2) = segmentA(1) + lineB(1,1)*segmentA(2) END IF END IF ELSE det = (a11*a22-a12*a21) ! Parallel lines ! IF (det == zeroRK) THEN intersect = .FALSE. ptintersect = zeroRK ELSE intersect = .TRUE. detX = (z1*a22-z2*a12) detY = (a11*z2-a21*z1) ptintersect(1) = detX/det ptintersect(2) = detY/det END IF END IF END SUBROUTINE intersection_2Dlines !refs: !https://www.mathopenref.com/heronsformula.html !https://math.stackexchange.com/questions/1406340/intersect-area-of-two-polygons-in-cartesian-plan !http://www.cap-lore.com/MathPhys/IP/ !http://www.cap-lore.com/MathPhys/IP/IL.html !https://www.element84.com/blog/determining-the-winding-of-a-polygon-given-as-a-set-of-ordered-points !https://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order !https://en.wikipedia.org/wiki/Shoelace_formula !https://en.wikipedia.org/wiki/Winding_number !https://en.wikipedia.org/wiki/Simple_polygon !https://en.wikipedia.org/wiki/Polygon#Properties !https://en.wikipedia.org/wiki/Convex_polygon !https://en.wikipedia.org/wiki/Jordan_curve_theorem !https://www.sangakoo.com/ca/temes/metode-de-cramer !https://www.geogebra.org/m/pw4QHFYT SUBROUTINE intersectfaces(faceA, faceB, intersect, intersectpt) ! Subroutine to provide if two faces of two polygons intersect ! AFTER: http://www.cap-lore.com/MathPhys/IP/IL.html ! A: faceA(1,:) ! B: faceA(2,:) ! C: faceB(1,:) ! D: faceB(2,:) IMPLICIT NONE REAL(r_k), DIMENSION(2,2), INTENT(in) :: faceA, faceB INTEGER, INTENT(out) :: intersect REAL(r_k), DIMENSION(2), INTENT(out) :: intersectpt ! Local REAL(r_k) :: Axmin, Aymin, Axmax, Aymax REAL(r_k) :: Bxmin, Bymin, Bxmax, Bymax REAL(r_k) :: areaABD, areaACD, areaBDC, areaDAB REAL(r_k), DIMENSION(3,2) :: triangle LOGICAL :: Lintersect !!!!!!! Variables ! faceA/B: coordinates of faces A and B to determine if they intersect ! intersect: integer to say if they intersect (0, no-intersect, +/-1 intersect) ! intersectpt: point where faces intersect [(0,0) otherwise] fname = 'intersectfaces' ! PRINT *,' ' // TRIM(fname) // ' ________' ! PRINT *,' faceA:', faceA(1,:), ';',faceA(2,:) ! PRINT *,' faceB:', faceB(1,:), ';',faceB(2,:) Axmin = MINVAL(faceA(:,1)) Axmax = MAXVAL(faceA(:,1)) Aymin = MINVAL(faceA(:,2)) Aymax = MAXVAL(faceA(:,2)) Bxmin = MINVAL(faceB(:,1)) Bxmax = MAXVAL(faceB(:,1)) Bymin = MINVAL(faceB(:,2)) Bymax = MAXVAL(faceB(:,2)) ! No intersection IF ( (Axmax <= Bxmin) .OR. (Axmin >= Bxmax) .OR. (Aymax <= Bymin) .OR. (Aymin >= Bymax) ) THEN intersect = 0 intersectpt = zeroRK ELSE ! Triangle ABD triangle(1,:) = faceA(1,:) triangle(2,:) = faceA(2,:) triangle(3,:) = faceB(2,:) areaABD = shoelace_area_polygon(3, triangle) ! Triangle ACD triangle(1,:) = faceA(1,:) triangle(2,:) = faceB(1,:) triangle(3,:) = faceB(2,:) areaACD = shoelace_area_polygon(3, triangle) ! Triangle BDC triangle(1,:) = faceA(2,:) triangle(2,:) = faceB(2,:) triangle(3,:) = faceB(1,:) areaBDC = shoelace_area_polygon(3, triangle) ! Triangle DAB triangle(1,:) = faceB(2,:) triangle(2,:) = faceA(1,:) triangle(3,:) = faceA(2,:) areaDAB = shoelace_area_polygon(3, triangle) IF (areaABD>zeroRK .AND. areaACD>zeroRK .AND. areaBDC>zeroRK .AND. areaDAB>zeroRK) THEN intersect = INT(ABS(areaABD)/areaABD) CALL intersection_2Dlines(faceA, faceB, Lintersect, intersectpt) ELSE IF (areaABD zeroRK, 'I', intersect END IF END SUBROUTINE intersectfaces LOGICAL FUNCTION poly_has_point(Nvertex, polygon, point) ! Function to determine if a polygon has already a given point as one of its vertex IMPLICIT NONE INTEGER, INTENT(in) :: Nvertex REAL(r_k), DIMENSION(Nvertex,2), INTENT(in) :: polygon REAL(r_k), DIMENSION(2), INTENT(in) :: point ! Local INTEGER :: iv REAL(r_k), DIMENSION(2) :: diff !!!!!!! Vertrex ! Nvertex: number of vertexs of the polygon ! polygon: vertexs of the polygon ! point: point to look for its ownership into the polygon fname = 'poly_has_point' poly_has_point = .FALSE. DO iv=1, Nvertex diff = polygon(iv,:)-point IF ( (diff(1) == zeroRK) .AND. (diff(2) == zeroRK)) THEN poly_has_point = .TRUE. EXIT END IF END DO END FUNCTION poly_has_point SUBROUTINE join_polygon(NvertexA, NvertexB, NvertexAB, polyA, polyB, Ncoinvertex, coinpoly) ! Subroutine to join two polygons ! AFTER: http://www.cap-lore.com/MathPhys/IP/ and http://www.cap-lore.com/MathPhys/IP/IL.html IMPLICIT NONE INTEGER, INTENT(in) :: NvertexA, NvertexB, NvertexAB REAL(r_k), DIMENSION(NvertexA,2), INTENT(in) :: polyA REAL(r_k), DIMENSION(NvertexB,2), INTENT(in) :: polyB INTEGER, INTENT(out) :: Ncoinvertex REAL(r_k), DIMENSION(NvertexAB,2), INTENT(out) :: coinpoly ! Local INTEGER :: iA, iB, icoin, ii REAL(r_k), DIMENSION(2,2) :: face1, face2 INTEGER :: intersct REAL(r_k), DIMENSION(2) :: ptintersct !!!!!!! variables ! NvertexA: number of vertexs polygon A ! NvertexB: number of vertexs polygon B ! polyA: pairs of coordinates for the polygon A (clockwise) ! polyB: pairs of coordinates for the polygon B (clockwise) ! Ncoinvertex: number of vertexes for the coincident polygon ! coinpoly: pairs of coordinates for the coincident polygon (clockwise) fname = 'join_polygon' icoin = 0 coinpoly = 0. ! First, include that vertex which do not lay within any polygon DO iA=1, NvertexA !PRINT *, ' iA:', iA, ':', polyA(iA,:), point_inside(polyA(iA,:), NvertexB, polyB) IF (.NOT. point_inside(polyA(iA,:), NvertexB, polyB)) THEN icoin = icoin + 1 coinpoly(icoin,:) = polyA(iA,:) END IF END DO DO iB=1, NvertexB !PRINT *, ' iB:', iB, ':', polyB(iB,:), point_inside(polyB(iB,:), NvertexA, polyA) IF (.NOT. point_inside(polyB(iB,:), NvertexA, polyA)) THEN icoin = icoin + 1 coinpoly(icoin,:) = polyB(iB,:) END IF END DO DO iA=1, NvertexA ! Getting couple of vertexs from polyA and polyB IF (iA /= NvertexA) THEN face1(1,:) = polyA(iA,:) face1(2,:) = polyA(iA+1,:) ELSE face1(1,:) = polyA(iA,:) face1(2,:) = polyA(1,:) END IF DO iB=1, NvertexB IF (iB /= NvertexB) THEN face2(1,:) = polyB(iB,:) face2(2,:) = polyB(iB+1,:) ELSE face2(1,:) = polyB(iB,:) face2(2,:) = polyB(1,:) END IF ! Compute areas of the four possible triangles. Introduce the coincident vertexs not included CALL intersectfaces(face1, face2, intersct, ptintersct) !PRINT *,iA,':',face1(1,:),';',face1(2,:), '=', iB, face2(1,:),';',face2(2,:), '<>', intersct,':', ptintersct IF (intersct == 1) THEN IF (.NOT.poly_has_point(icoin,coinpoly(1:icoin,:),ptintersct) ) THEN icoin = icoin + 1 coinpoly(icoin,:) = ptintersct END IF ELSE IF (intersct == -1) THEN IF (.NOT.poly_has_point(icoin,coinpoly(1:icoin,:),ptintersct) ) THEN icoin = icoin + 1 coinpoly(icoin,:) = ptintersct END IF END IF END DO END DO Ncoinvertex = icoin END SUBROUTINE join_polygon SUBROUTINE sort_polygon(Nvertex, polygon, sense, Nnewvertex, newpoly) ! Subroutine to sort a polygon using its center as average of the coordinates and remove duplicates ! Should be used the centroid instead, but by now let do it simple ! https://en.wikipedia.org/wiki/Centroid IMPLICIT NONE INTEGER, INTENT(in) :: Nvertex, sense REAL(r_k), DIMENSION(Nvertex,2), INTENT(in) :: polygon INTEGER, INTENT(out) :: Nnewvertex REAL(r_k), DIMENSION(Nvertex,2), INTENT(out) :: newpoly ! Local INTEGER :: iv, j REAL(r_k) :: vang REAL(r_k), DIMENSION(2) :: center REAL(r_k), DIMENSION(Nvertex) :: angles REAL(r_k), DIMENSION(Nvertex,2) :: sortpoly !!!!!!! Variables ! Nvertex: number of vertices ! polygon: coordinates of the vertices of the polygon ! sense: sens of sorting thepolygon (1: clockwise, -1: anti-clockwise) ! sortpoly: sorted polygon ! Nnewvertex: number of vertices new polygon ! newpoly: sorted and duplicate removed polygon fname = 'sort_polygon' ! To be substituted by centroid calculation (which requires already sorted vetexs...) center(1) = SUM(polygon(:,1))/Nvertex center(2) = SUM(polygon(:,2))/Nvertex DO iv=1, Nvertex angles(iv) = ATAN2(polygon(iv,2)-center(2),polygon(iv,1)-center(1)) END DO CALL sortR_K(angles, Nvertex) sortpoly = zeroRK DO iv=1, Nvertex DO j=1, Nvertex vang = ATAN2(polygon(j,2)-center(2), polygon(j,1)-center(1)) IF (angles(iv) == vang) THEN IF (sense == -1) THEN sortpoly(iv,:) = polygon(j,:) ELSE sortpoly(Nvertex-iv+1,:) = polygon(j,:) END IF EXIT END IF END DO END DO newpoly(1,:) = sortpoly(1,:) j = 1 DO iv=2, Nvertex IF (.NOT.poly_has_point(j,newpoly(1:j,:),sortpoly(iv,:)) ) THEN j = j+1 newpoly(j,:) = sortpoly(iv,:) END IF END DO Nnewvertex = j END SUBROUTINE sort_polygon LOGICAL FUNCTION point_inside(point, Nvertex, polygon) ! Function to determine if a given point is inside a polygon providing its sorted vertices ! FROM: https://en.wikipedia.org/wiki/Point_in_polygon IMPLICIT NONE REAL(r_k), DIMENSION(2), INTENT(in) :: point INTEGER, INTENT(in) :: Nvertex REAL(r_k), DIMENSION(Nvertex,2), INTENT(in) :: polygon ! Local INTEGER :: iv, Nintersect INTEGER :: cross REAL(r_k) :: xmin REAL(r_k), DIMENSION(2) :: crosspoint REAL(r_k), DIMENSION(2,2) :: face1, face2 REAL(r_k), DIMENSION(Nvertex) :: abovebelow !!!!!!! Variables ! point: point to look for ! Nvertrex: number of vertices of a polygon ! polygon: vertices of a polygon fname = 'point_inside' xmin = MINVAL(polygon(:,1)) ! Looking for the intersection with the ray Nintersect = 0 face1(1,:) = (/ xmin-0.5, point(2) /) face1(2,:) = (/ point(1), point(2) /) DO iv = 1, Nvertex IF (iv /= Nvertex) THEN face2(1,:) = polygon(iv,:) face2(2,:) = polygon(iv+1,:) ELSE face2(1,:) = polygon(iv,:) face2(2,:) = polygon(1,:) END IF CALL intersectfaces(face1, face2, cross, crosspoint) IF (cross /= 0) THEN Nintersect = Nintersect + 1 abovebelow(Nintersect) = iv END IF END DO IF (MOD(Nintersect,2) == 0) THEN point_inside = .FALSE. ELSE point_inside = .TRUE. END IF END FUNCTION point_inside LOGICAL FUNCTION point_in_face(pt, Nvertex, poly) ! Function to determine if a given point is on a face of a polygon IMPLICIT NONE REAL(r_k), DIMENSION(2), INTENT(in) :: pt INTEGER, INTENT(in) :: Nvertex REAL(r_k), DIMENSION(Nvertex,2), INTENT(in) :: poly ! Local INTEGER :: iv REAL(r_k) :: ix, ex, iy, ey, tmpv REAL(r_k) :: dx, dy, A, B !!!!!!! Variables ! pt: point to look for ! Nvertex: Number of vertices of the polygon ! poly: polygon fname = 'point_in_face' point_in_face = .FALSE. DO iv=1, Nvertex IF (iv < Nvertex) THEN ix = poly(iv,1) ex = poly(iv+1,1) iy = poly(iv,2) ey = poly(iv+1,2) ELSE ix = poly(iv,1) ex = poly(1,1) iy = poly(iv,2) ey = poly(1,2) END IF dx = ex - ix dy = ey - iy IF (dx == zeroRK) THEN IF (pt(1) == ix) THEN IF ( (iy < ey) .AND. (pt(2) >= iy) .AND. pt(2) <= ey) THEN point_in_face = .TRUE. EXIT ELSE IF ( (iy > ey) .AND. (pt(2) >= ey) .AND. pt(2) <= iy) THEN point_in_face = .TRUE. EXIT END IF END IF ELSE IF (dy == zeroRK) THEN IF (pt(2) == iy) THEN IF ((ix < ex) .AND. (pt(1) >= ix) .AND. pt(1) <= ex) THEN point_in_face = .TRUE. EXIT ELSE IF ((ix > ex) .AND. (pt(1) >= ex) .AND. pt(1) <= ix) THEN point_in_face = .TRUE. EXIT END IF END IF ELSE A = iy B = (ey-iy)/(ex-ix) IF (A+B*(pt(1)-ix) == pt(2)) THEN point_in_face = .TRUE. EXIT END IF END IF END IF END DO END FUNCTION point_in_face SUBROUTINE coincident_polygon(NvertexA, NvertexB, NvertexAB, polyA, polyB, Ncoinvertex, coinpoly) ! Subroutine to provide the intersection polygon between two polygons ! AFTER: http://www.cap-lore.com/MathPhys/IP/ and http://www.cap-lore.com/MathPhys/IP/IL.html IMPLICIT NONE INTEGER, INTENT(in) :: NvertexA, NvertexB, NvertexAB REAL(r_k), DIMENSION(NvertexA,2), INTENT(in) :: polyA REAL(r_k), DIMENSION(NvertexB,2), INTENT(in) :: polyB INTEGER, INTENT(out) :: Ncoinvertex REAL(r_k), DIMENSION(NvertexAB,2), INTENT(out) :: coinpoly ! Local INTEGER :: iA, iB, icoin, ii REAL(r_k), DIMENSION(2,2) :: face1, face2 INTEGER :: intersct REAL(r_k), DIMENSION(2) :: ptintersct !!!!!!! variables ! NvertexA: number of vertexs polygon A ! NvertexB: number of vertexs polygon B ! polyA: pairs of coordinates for the polygon A (clockwise) ! polyB: pairs of coordinates for the polygon B (clockwise) ! Ncoinvertex: number of vertexes for the coincident polygon ! coinpoly: pairs of coordinates for the coincident polygon (clockwise) fname = 'coincident_polygon' icoin = 0 coinpoly = 0. ! First, include that vertex which lay within any polygon DO iA=1, NvertexA IF (point_inside(polyA(iA,:), NvertexB, polyB)) THEN icoin = icoin + 1 coinpoly(icoin,:) = polyA(iA,:) END IF IF (point_in_face(polyA(iA,:), NvertexB, polyB)) THEN icoin = icoin + 1 coinpoly(icoin,:) = polyA(iA,:) END IF END DO DO iB=1, NvertexB IF (point_inside(polyB(iB,:), NvertexA, polyA)) THEN icoin = icoin + 1 coinpoly(icoin,:) = polyB(iB,:) END IF IF (point_in_face(polyB(iB,:), NvertexA, polyA)) THEN icoin = icoin + 1 coinpoly(icoin,:) = polyB(iB,:) END IF END DO ! Look interesections DO iA=1, NvertexA ! Getting couple of vertexs from polyA and polyB IF (iA /= NvertexA) THEN face1(1,:) = polyA(iA,:) face1(2,:) = polyA(iA+1,:) ELSE face1(1,:) = polyA(iA,:) face1(2,:) = polyA(1,:) END IF DO iB=1, NvertexB IF (iB /= NvertexB) THEN face2(1,:) = polyB(iB,:) face2(2,:) = polyB(iB+1,:) ELSE face2(1,:) = polyB(iB,:) face2(2,:) = polyB(1,:) END IF ! Compute areas of the four possible triangles. Introduce the coincident vertexs not included CALL intersectfaces(face1, face2, intersct, ptintersct) !PRINT *,iA,':',face1(1,:),';',face1(2,:), '=', iB, face2(1,:),';',face2(2,:), '<>', intersct,':', ptintersct IF ((intersct /= 0) .AND. (.NOT.poly_has_point(icoin,coinpoly(1:icoin,:),ptintersct)) ) THEN icoin = icoin + 1 coinpoly(icoin,:) = ptintersct END IF END DO END DO Ncoinvertex = icoin END SUBROUTINE coincident_polygon SUBROUTINE grid_within_polygon(NvertexA, polygonA, dx, dy, dxy, xCvals, yCvals, Nvertexmax, xBvals, & yBvals, Ngridsin, gridsin) ! Subroutine to determine which grid cells from a matrix lay inside a polygon IMPLICIT NONE INTEGER, INTENT(in) :: NvertexA, dx, dy, dxy, Nvertexmax REAL(r_k), DIMENSION(NvertexA,2), INTENT(in) :: polygonA REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: xCvals, yCvals REAL(r_k), DIMENSION(dx,dy,Nvertexmax), INTENT(in) :: xBvals, yBvals INTEGER, INTENT(out) :: Ngridsin INTEGER, DIMENSION(dxy,2), INTENT(out) :: gridsin ! Local INTEGER :: ix, iy, iv REAL(r_k), DIMENSION(2) :: centergrid, vertex LOGICAL, DIMENSION(dx,dy) :: within !!!!!!! Variables ! NvertexA: Number of vertices of the polygin to find the grids ! polygonA: ordered vertices of the polygon ! dx, dy: shape of the matrix with the grid points ! xCvals, yCvals: coordinates of the center of the grid cells ! Nvertexmax: Maximum number of vertices of the grid cells ! xBvals, yBvals: coordinates of th vertices of the grid cells (-99999 for no vertex) ! Ngridsin: number of grids with some extension within the polygon ! gridsin: grids within the polygin ! percentages: percentages of area of each of the grids within the polygon fname = 'spacepercen_within_reg' Ngridsin = 0 gridsin = 0 within = .FALSE. DO ix = 1, dx DO iy = 1, dy IF (.NOT.within(ix,iy)) THEN centergrid = (/ xCvals(ix,iy), yCvals(ix,iy) /) ! By grid center IF (point_inside(centergrid, NvertexA, polygonA)) THEN Ngridsin = Ngridsin + 1 ! Getting coordinates gridsin(Ngridsin,1) = ix gridsin(Ngridsin,2) = iy within(ix,iy) = .TRUE. CYCLE END IF ! Getting grid vertices DO iv=1, Nvertexmax IF (.NOT.within(ix,iy)) THEN IF (xBvals(ix,iy,iv) /= fillvalI) THEN vertex = (/ xBvals(ix,iy,iv), yBvals(ix,iy,iv) /) IF (point_inside(vertex, NvertexA, polygonA)) THEN Ngridsin = Ngridsin + 1 ! Getting coordinates gridsin(Ngridsin,1) = ix gridsin(Ngridsin,2) = iy within(ix,iy) = .TRUE. CYCLE END IF END IF END IF END DO END IF END DO END DO END SUBROUTINE grid_within_polygon SUBROUTINE spacepercen_within_reg(NvertexA, polygonA, dx, dy, Nvertexmax, xBvals, yBvals, & Ngridsin, gridsin, strict, percentages) ! Subroutine to compute the percentage of a series of grid cells which are encompassed by a polygon ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x ! and y axes. IMPLICIT NONE INTEGER, INTENT(in) :: NvertexA, dx, dy, Nvertexmax REAL(r_k), DIMENSION(NvertexA,2), INTENT(in) :: polygonA REAL(r_k), DIMENSION(dx,dy,Nvertexmax), INTENT(in) :: xBvals, yBvals INTEGER, INTENT(in) :: Ngridsin INTEGER, DIMENSION(Ngridsin,2), INTENT(in) :: gridsin LOGICAL, INTENT(in) :: strict REAL(r_k), DIMENSION(Ngridsin), INTENT(out) :: percentages ! Local INTEGER :: ig, iv, ix, iy INTEGER :: Nvertex, NvertexAgrid, Ncoin, Nsort CHARACTER(len=20) :: DS REAL(r_k) :: areapoly, areagpoly, totarea, totpercent REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: vertexgrid, icoinpoly, coinpoly, & sortpoly, poly !!!!!!! Variables ! NvertexA: Number of vertices of the polygin to find the grids ! polygonA: ordered vertices of the polygon ! dx, dy: shape of the matrix with the grid points ! xCvals, yCvals: coordinates of the center of the grid cells ! Nvertexmax: Maximum number of vertices of the grid cells ! xBvals, yBvals: coordinates of th vertices of the grid cells (-99999 for no vertex) ! Ngridsin: number of grids with some extension within the polygon ! gridsin: grids within the polygon ! strict: give an error if the area of the polygon is not fully covered ! percentages: percentages of area of each of the grids within the polygon fname = 'spacepercen_within_reg' percentages = zeroRK totpercent = zeroRK totarea = zeroRK areapoly = shoelace_area_polygon(NvertexA, polygonA) DO ig = 1, Ngridsin ix = gridsin(ig,1) iy = gridsin(ig,2) ! Getting grid vertices Nvertex = 0 DO iv=1, Nvertexmax IF (xBvals(ix,iy,iv) /= fillvalI) THEN Nvertex = Nvertex + 1 END IF END DO IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) ALLOCATE(vertexgrid(Nvertex,2)) vertexgrid(:,1) = xBvals(ix,iy,1:Nvertex) vertexgrid(:,2) = yBvals(ix,iy,1:Nvertex) ! Getting common vertices NvertexAgrid = NvertexA*Nvertex IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly) ALLOCATE(icoinpoly(NvertexAgrid,2)) CALL coincident_polygon(NvertexA, Nvertex, NvertexAgrid, polygonA, vertexgrid, Ncoin, icoinpoly) IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly) ALLOCATE(coinpoly(Ncoin,2)) DO iv=1, Ncoin coinpoly(iv,:) = icoinpoly(iv,:) END DO IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly) ALLOCATE(sortpoly(Ncoin,2)) CALL sort_polygon(Ncoin, coinpoly, 1, Nsort, sortpoly) IF (ALLOCATED(poly)) DEALLOCATE(poly) ALLOCATE(poly(Nsort,2)) DO iv=1, Nsort poly(iv,:) = sortpoly(iv,:) END DO areagpoly = shoelace_area_polygon(Nsort, poly) IF (INT(LOG10(EPSILON(totpercent))) < 12) THEN totarea = totarea + ABS(areagpoly) percentages(ig) = ABS(areagpoly / areapoly) ! f2py does not like it! ! ELSE ! totarea = totarea + DABS(areagpoly) ! percentages(ig) = DABS(areagpoly / areapoly) END IF totpercent = totpercent + percentages(ig) END DO IF (INT(LOG10(EPSILON(totpercent))) < 12) THEN IF (strict .AND. ABS(totpercent - oneRK) > epsilonRK) THEN PRINT *, 'totarea:', totarea, ' area polygon:', areapoly PRINT *, 'totpercent:', totpercent, ' oneRK:', oneRK, ' diff:', totpercent - oneRK WRITE(DS,'(F20.8)')ABS(totpercent - oneRK) msg = 'sum of all grid space percentages does not cover (' // TRIM(DS) // ') all polygon' CALL ErrMsg(msg, fname, -1) END IF ELSE ! f2py does not like it! ! IF (strict .AND. ABS(totpercent - oneRK) > epsilonRK) THEN ! PRINT *, 'totarea:', totarea, ' area polygon:', areapoly ! PRINT *, 'totpercent:', totpercent, ' oneRK:', oneRK, ' diff:', totpercent - oneRK ! WRITE(DS,'(F20.16)')ABS(totpercent - oneRK) ! msg = 'sum of all grid space percentages does not cover (' // TRIM(DS) // ') all polygon' ! CALL ErrMsg(msg, fname, -1) ! END IF END IF END SUBROUTINE spacepercen_within_reg SUBROUTINE spacepercen(dxA, dyA, xCAvals, yCAvals, NAvertexmax, xBAvals, yBAvals, dxB, dyB, xCBvals,& yCBvals, NBvertexmax, xBBvals, yBBvals, dxyB, strict, Ngridsin, gridsin, percentages) ! Subroutine to compute the space-percentages of a series of grid cells (B) into another series of ! grid-cells (A) ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x ! and y axes. IMPLICIT NONE INTEGER, INTENT(in) :: dxA, dyA, NAvertexAma, dxB, dyB, & NBvertexmax, dxyB REAL(r_k), DIMENSION(dxA,dyA), INTENT(in) :: xCAvals, yCAvals REAL(r_k), DIMENSION(dxB,dyB), INTENT(in) :: xCBvals, yCBvals REAL(r_k), DIMENSION(dxA,dyA,NAvertexmax), INTENT(in):: xBAvals, yBAvals REAL(r_k), DIMENSION(dxB,dyB,NBvertexmax), INTENT(in):: xBBvals, yBBvals LOGICAL, INTENT(in) :: strict INTEGER, DIMENSION(dxA,dyA), INTENT(out) :: Ngridsin INTEGER, DIMENSION(dxA,dyA,dxyB,2), INTENT(out) :: gridsin REAL(r_k), DIMENSION(dxA,dyA,Ngridsin), INTENT(out) :: percentages ! Local INTEGER :: iv, ix, iy INTEGER :: Nvertex, NvertexAgrid, Ncoin, Nsort CHARACTER(len=20) :: DS REAL(r_k) :: areapoly, areagpoly, totarea, totpercent REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: vertexgrid, poinsin !!!!!!! Variables ! dxA, dyA: shape of the matrix with the grid points A ! xCAvals, yCAvals: coordinates of the center of the grid cells A ! NAvertexmax: Maximum number of vertices of the grid cells A ! xBAvals, yBAvals: coordinates of th vertices of the grid cells A (-99999 for no vertex) ! dxB, dyB: shape of the matrix with the grid points B ! xCBvals, yCBvals: coordinates of the center of the grid cells B ! NBvertexmax: Maximum number of vertices of the grid cells B ! xBBvals, yBBvals: coordinates of th vertices of the grid cells B (-99999 for no vertex) ! strict: give an error if the area of the polygon is not fully covered ! Ngridsin: number of grids from grid B with some extension within the grid cell A ! gridsin: indices of B grids within the grids of A ! percentages: percentages of area of cells B of each of the grids within the grid cell A fname = 'spacepercen' DO ix = 1, dimxA DO iy = 1, dimyA ! Getting grid vertices Nvertex = 0 DO iv=1, NAvertexmax IF (xBAvals(ix,iy,iv) /= fillvalI) THEN Nvertex = Nvertex + 1 END IF END DO IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) ALLOCATE(vertexgrid(Nvertex,2)) vertexgrid(:,1) = xBAvals(ix,iy,1:Nvertex) vertexgrid(:,2) = yBAvals(ix,iy,1:Nvertex) CALL grid_within_polygon(Nvertex, vertexgrid, dxB, dyB, dxB*dyB, xBCvals, yBCvals, NBvertexmax, & xBBvals, yBBvals, Ngridsin(ix,iy), gridsin(ix,iy,:,:)) IF (ALLOCATE(poinsin)) DEALLOCATE(poinsin) ALLOCATE(poinsin(Ngridsin(ix,iy),2)) DO iv=1, Ngridsin(ix,iy) poinsin(i,1) = gridsin(ix,iy,iv,1) poinsin(i,2) = gridsin(ix,iy,iv,2) END DO CALL spacepercen_within_reg(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, xBBvals, yBBvals, & Ngridsin(ix,iy), poinsin, strict, percentages(ix,iy,:)) END DO END DO END MODULE module_scientific