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 ! 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 ! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon ! following ray casting algorithm ! 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 ! 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 ! 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 END MODULE module_scientific