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 ! coincident_gridsin2D: Subroutine to determine which lists of 2D gridsin points of an A list are also found in a B list ! coincident_list_2Dcoords: Subroutine to determine which 2D points of an A list are also found in a B list ! clean_polygons: Subroutine to clean polygons from non-used paths, polygons only left as path since they are inner path of a hole ! coincident_polygon: Subroutine to provide the intersection polygon between two polygons ! crossingpoints_polys: Subroutine to provide the crossing points between two polygons ! distanceRK: Function to provide the distance between two points ! FindMinimumR_K*: Function returns the location of the minimum in the section between Start and End. ! fill3DI_2Dvec: Subroutine to fill a 3D integer matrix from a series of indices from a 2D matrix ! fill3DR_2Dvec: Subroutine to fill a 3D float matrix from a series of indices from a 2D matrix ! get_multindices_lonlat: Subroutine to provide the nearest grid point of a series of longitudes and latitudes values ! grid_within_polygon: Subroutine to determine which grid cells from a matrix lay inside a polygon ! grid_spacepercen: Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another ! series of grid-cells (A) ! grid_spacepercen_providing_polys: Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another ! series of grid-cells (A) providing coincident polygons ! grid_spacepercen_within_reg: Subroutine to compute the percentage of grid space of a series of grid cells which are encompassed ! by a polygon ! grid_spacepercen_within_reg_providing_polys: Subroutine to compute the percentage of grid space of a series of grid ! cells which are encompassed by a polygon providing coordinates of the resultant polygons ! gridpoints_InsidePolygon: Subroutine to determine if a series of grid points are inside a polygon ! following ray casting algorithm ! Heron_area_triangle: Subroutine to compute area of a triangle using Heron's formula ! HistogramRK: Subroutine to provide the histogram of a series of RK values ! Histogram2D1_RK: Subroutine to provide the histogram of a 2D matrix RK values along the first dimension ! intersectfaces: Subroutine to provide if two faces of two polygons intersect ! intersection_2Dlines: Subroutine to provide the intersection point between two lines on the plane using Cramer's method ! join_polygon: Subroutine to join two polygons ! look_clockwise_borders: Subroutine to look clock-wise for a next point within a collection of borders ! (limits of a region) ! multi_index_mat2DI: Subroutine to provide the indices of the different locations of a value inside a 2D integer matrix ! multi_index_mat3DI: Subroutine to provide the indices of the different locations of a value inside a 3D integer matrix ! multi_index_mat4DI: Subroutine to provide the indices of the different locations of a value inside a 4D integer matrix ! multi_index_mat2DRK: Subroutine to provide the indices of the different locations of a value inside a 2D RK matrix ! multi_index_mat3DRK: Subroutine to provide the indices of the different locations of a value inside a 3D RK matrix ! multi_index_mat4DRK: Subroutine to provide the indices of the different locations of a value inside a 4D RK matrix ! multi_spaceweightcount_in3DRK3_slc3v3: Subroutine to compute an spatial percentage of coverture of ! categories from a 3D RK matrix using 3rd dimension as running one into a matrix of 3-variables ! slices of rank 3 using spatial weights ! multi_spaceweightstats_in1DRKno_slc3v3: Subroutine to compute an spatial statistics value from a 1D RK matrix without ! running one into a matrix of 3-variables slices of rank 3 using spatial weights ! multi_spaceweightstats_in2DRKno_slc3v3: Subroutine to compute an spatial statistics value from a 2D RK matrix without ! running one into a matrix of 3-variables slices of rank 3 using spatial weights ! multi_spaceweightstats_in2DRKno_slc2v2: Subroutine to compute an spatial statistics value from a 2D RK matrix without ! running one into a matrix of 2-variables slices of rank 2 using spatial weights ! multi_spaceweightstats_in3DRK3_slc2v2: Subroutine to compute an spatial statistics value from a 3D RK matrix using ! 3rd dimension as running one into a matrix of 2-variables slices of rank 2 using spatial weights ! multi_spaceweightstats_in3DRK3_slc3v3: Subroutine to compute an spatial statistics value from a 3D RK matrix using ! 3rd dimension as running one into a matrix of 3-variables slices of rank 3 using spatial weights ! multi_spaceweightstats_in3DRK3_slc3v4: Subroutine to compute an spatial statistics value from a 3D RK matrix using ! 3rd dimension as running one into a matrix of 3-variables slices of rank 4 using spatial weights ! multi_spaceweightstats_in4DRK3_4_slc3v3: Subroutine to compute an spatial statistics value from a 4D RK matrix using ! 3rd and 4th dimensions as running ones into a matrix of 3-variables slices of rank 3 using spatial weights ! NcountR: Subroutine to count real values ! paths_border: Subroutine to search the paths of a border field. ! path_properties: Subroutine to determine the properties of a path ! percentiles_R_K*D: Subroutine to compute the percentiles of a *D R_K array along given set of axis ! point_in_face: Function to determine if a given point is on a face of a polygon ! point_inside: Function to determine if a given point is inside a polygon providing its sorted vertices ! poly_has_point: Function to determine if a polygon has already a given point as one of its vertex ! polygon_properties: Subroutine to determine the properties of a polygon (as .TRUE. matrix) ! polygons: Subroutine to search the polygons of a border field. FORTRAN based. 1st = 1! ! polygons_t: Subroutine to search the polygons of a temporal series of boolean fields. FORTRAN based. 1st = 1! ! poly_overlap_tracks: Subroutine to determine tracks of a series of consecutive 2D field with polygons using ! maximum overlaping/coincidence ! PrintQuantilesR_K: Subroutine to print the quantiles of values REAL(r_k) ! poly_overlap_tracks_area: Subroutine to determine tracks of a series of consecutive 2D field with polygons using ! maximum overlaping/coincidence filtrered by a minimal area ! poly_overlap_tracks_area_ascii: Subroutine to determine tracks of a series of consecutive 2D field with polygons using maximum ! overlaping/coincidence filtrered by a minimal area writting theoutput on an ASCII file (memory limitations) ! quantilesR_K: Subroutine to provide the quantiles of a given set of values of type real 'r_k' ! rand_sample: Subroutine to randomly sample a range of indices ! read_finaltrack_ascii: Subroutine to read the final trajectories from an ASCII file ! read_overlap_single_track_ascii: Subroutine to read the values for a given trajectory ! read_overlap_polys_ascii: Subroutine to read from an ASCII file the associated polygons at a given time-step ! read_overlap_tracks_ascii: Subroutine to write to an ASCII the polygons associated to a trajectory at a given time step ! reconstruct_matrix: Subroutine to reconstruct a 2D matrix from a pair of syncronized vectors with the positions on x ! and y coordinates ! runmean_F1D: Subroutine fo computing the running mean of a given set of float 1D values ! shoelace_area_polygon: Computing the area of a polygon using sholace formula ! sort_polygon: Subroutine to sort a polygon using its center as average of the coordinates and remove duplicates ! SortR_K*: Subroutine receives an array x() r_K and sorts it into ascending order. ! spacepercen: Subroutine to compute the space-percentages of a series of grid cells (B) into another series of grid-cells (A) ! spacepercen_within_reg: Subroutine to compute the percentage of a series of grid cells which are encompassed by a polygon ! spaceweight_Icount: Subroutine to compute the percentage of a series of categories using weights ! spaceweightstats: Subroutine to compute an spatial statistics value from a matrix B into a matrix A using weights ! 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. ! unique_matrixRK2D: Subroutine to provide the unique values within a 2D RK matrix ! 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 !! Vectorial calculus ! curl2D_1o: Subroutine to compute the first order curl of a 2D vectorial field ! curl2D_c1o: Subroutine to compute the first order centered curl of a 2D vectorial field ! deformation3D: Subroutine to compute the deformation of a 3D field ! deltat3D: Subroutine to compute the temporal derivative of a 3D field ! divergence2D_1o: Subroutine to compute the first order divergence of a 2D vectorial field ! divergence2D_c1o: Subroutine to compute the first order centered divergence of a 2D vectorial field ! gradient2D_1o: Subroutine to compute the first order gradient of a 2D field ! gradient3D_1o: Subroutine to compute the first order gradient of a 3D field ! gradient2D_c1o: Subroutine to compute the first order centered gradient of a 2D field ! lap2D_1o: Subroutine to compute the first order laplacian of a 2D vectorial field ! lap2D_c1o: Subroutine to compute the first order centered laplacian of a 2D vectorial field ! matmodule3D: Subroutine to compute the module of a 3D matrix with 3 components ! tilting3D: Subroutine to compute the tilting of a 3D field ! vecmodule3D: Function to compute the module of a 3D vector !!! *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, k, ip, ipp, Nppt INTEGER :: ierr INTEGER, DIMENSION(:,:), ALLOCATABLE :: borders LOGICAL, DIMENSION(dx,dy) :: isborder, isbordery, borderp 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 CHARACTER(len=1000) :: boundsS !!!!!!! 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/100 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) boundsS = vectorI_S(3, (/Nppt, Nppt, 2/)) msg = "Problems allocating matrix 'paths' shape: " // TRIM(boundsS) // " try to reduce Nppt " // & "and recompile" 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(dy,2), STAT=ierr) msg = "Problems allocating matrix 'points'" CALL ErrMsg(msg, fname, ierr) CALL borders_matrixL(dbg, 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. borderp = .FALSE. DO j=1,Nptpaths(ip) borderp(paths(ip,j,1),paths(ip,j,2)) = .TRUE. END DO 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(dbg, dx, dy, isbordery, Nptpaths(ip), paths(ip,1:Nptpaths(ip),:), & ! Nvertx, xtrx, xtry, vertxs, Npts, points, isin) ! We only want to localize that points 'inside' Npts = 1 DO i=xtrx(1), xtrx(2) DO j=1, dy IF (boolmat(i,j)) THEN points(Npts,1) = i points(Npts,2) = j Npts = Npts + 1 END IF END DO CALL gridpoints_InsidePolygon_ray(dbg,dy,borderp(i,:), 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 *,' path boolmat isborder isborderp polygon (', xtrx(1), ',', xtry(1), ')x(', xtrx(2), & ',', xtry(2), ') _______' , Npts PRINT *,ip,'<>', i,':',boolmat(i,xtry(1):xtry(2)), ' border ', isborder(i,xtry(1):xtry(2)), & ' isbordery ', borderp(i,xtry(1):xtry(2)), ' polygon ', polys(i,xtry(1):xtry(2)) END IF END DO END DO PRINT *,' PRE-clean ' ! Cleaning polygons matrix of not-used and paths around holes !CALL clean_polygons(dx, dy, boolmat, polys, Npoly, dbg) IF (ALLOCATED(borders)) DEALLOCATE (borders) IF (ALLOCATED(Nptpaths)) DEALLOCATE (Nptpaths) IF (ALLOCATED(paths)) DEALLOCATE (paths) IF (ALLOCATED(vertxs)) DEALLOCATE (vertxs) IF (ALLOCATED(points)) DEALLOCATE (points) IF (ALLOCATED(isin)) 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 (NnotPol == Npols) THEN PRINT *,' ' // TRIM(fname) // ": avoiding to remove all polygons !!" NnotPol = 0 END IF 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(dbg, 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, INTENT(in) :: dbg 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. IF (dbg) PRINT *,'Border _______' DO i=1,dx halo_brdr(i+1,2:dy+1) = isbrdr(i,:) IF (dbg) PRINT *,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 IF (dbg) PRINT *,' ',Nbrbrdr,' Consec brdrs:', halo_brdr(ix+1,j+1), halo_brdr(ix+1,j+2), & '(', ix,j,';', ix,j+1,')', isbrdr(ix,j), isbrdr(ix,j+1) ELSE ! Will remove that consecutive borders above 2 IF (Nbrbrdr /= 0) THEN IF (dbg) PRINT *, ix,',',iy,';', Nintersecs, ' amount of consecutive borders:', Nbrbrdr, & ' removing:', MAX(Nbrbrdr-1, 0) Nintersecs = Nintersecs - MAX(Nbrbrdr-1, 0) Nbrbrdr = 0 END IF END IF END DO IF (MOD(Nintersecs,2) /= 0) inside(ip) = .TRUE. IF (dbg) PRINT *,ip,' point:', ix, iy, 'isbrdr:', isbrdr(ix,1:iy-1), 'y-ray:', halo_brdr(ix+1,1:iy), 'inside:', inside(ip) END IF END DO RETURN END SUBROUTINE gridpoints_InsidePolygon SUBROUTINE gridpoints_InsidePolygon_ray(dbg, 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) :: dy,Npath,Nvrtx,Npts LOGICAL, INTENT(in) :: dbg LOGICAL, DIMENSION(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(dy,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 ! dy: y-axis 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_ray' ! 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(dy+2), STAT=ierr) msg = "Problems allocating matrix 'halo_brdr'" CALL ErrMsg(msg, fname, ierr) halo_brdr = .FALSE. IF (dbg) PRINT *,'Border _______' halo_brdr(2:dy+1) = isbrdr(:) IF (dbg) PRINT *,isbrdr(:) inside = .FALSE. DO ip=1,dy 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(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(j+1) .AND. (ispath /= -1) .AND. (isvertex == -1) ) Nintersecs = Nintersecs + 1 IF (halo_brdr(j+1) .AND. (ispath /= -1) .AND. (halo_brdr(j+1) .EQV. halo_brdr(j+2))) THEN Nbrbrdr = Nbrbrdr + 1 IF (dbg) PRINT *,' ',Nbrbrdr,' Consec brdrs:', halo_brdr(j+1), halo_brdr(j+2), '(', & ix,j,';', ix,j+1,')', '(', ix,j,';', ix,j+1,')', isbrdr(j), isbrdr(j+1) ELSE ! Will remove that consecutive borders above 2 IF (Nbrbrdr /= 0) THEN IF (dbg) PRINT *, ix,',',iy,';', Nintersecs, ' amount of consecutive borders:', Nbrbrdr, & ' removing:', MAX(Nbrbrdr-1, 0) Nintersecs = Nintersecs - MAX(Nbrbrdr-1, 0) Nbrbrdr = 0 END IF END IF END DO IF (MOD(Nintersecs,2) /= 0) inside(ip) = .TRUE. IF (dbg) PRINT *,ip,' point:', ix, iy, 'isbrdr:', isbrdr(1:iy-1), 'y-ray:', halo_brdr(1:iy), & 'inside:', inside(ip) END IF END DO RETURN END SUBROUTINE gridpoints_InsidePolygon_ray 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(dbg,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, INTENT(in) :: dbg 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 IF (dbg) THEN PRINT *,' BORDERS _______ x y' DO i=1,dx PRINT *,isbrdr(i,:), ' ', isbrdry(i,:) END DO END IF 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' IF (dbg) PRINT *, TRIM(fname) // ' ...' ! 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 *,' isborder ______' DO i=1,dx PRINT *,isborder(i,:) END DO 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,dbg,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, dbg, iif, & jjf,iff) IF (iif /= -1 .AND. iff /= -1) THEN IF (dbg) PRINT *,' re-take path from point:', iif,',',jjf,' n-path:', iff found = .TRUE. iipth = index_list_coordsI(Nppt, borders, (/iip,ijp/)) EXIT END IF END DO IF (.NOT.found) THEN ! Looking for the next available border point for the new path DO i=1,Nbrdr IF (.NOT.gotbrdr(i)) THEN iipth = i EXIT END IF END DO IF (dbg) PRINT *,' Looking for next path starting at:', iipth, ' point:', & borders(iipth,:) ipath=ipath+1 EXIT END IF END IF ELSE IF (dbg) PRINT *,' looking for next point...' END IF IF (finishedstep) EXIT END DO END IF newpath END DO Npath = ipath Nptpaths(ipath) = Nptspath DEALLOCATE (gotbrdr) DEALLOCATE (emptygotbrdr) RETURN END SUBROUTINE paths_border SUBROUTINE rand_sample(Nvals, Nsample, sample) ! Subroutine to randomly sample a range of indices IMPLICIT NONE INTEGER, INTENT(in) :: Nvals, Nsample INTEGER, DIMENSION(Nsample), INTENT(out) :: sample ! Local INTEGER :: i, ind, jmax REAL, DIMENSION(Nsample) :: randv CHARACTER(len=50) :: fname LOGICAL :: found LOGICAL, DIMENSION(Nvals) :: issampled CHARACTER(len=256) :: msg CHARACTER(len=10) :: IS1, IS2 !!!!!!! Variables ! Nvals: number of values ! Nsamples: number of samples ! sample: samnple fname = 'rand_sample' IF (Nsample > Nvals) THEN WRITE(IS1,'(I10)')Nvals WRITE(IS2,'(I10)')Nsample msg = 'Sampling of ' // TRIM(IS1) // ' is too big for ' // TRIM(IS1) // 'values' CALL ErrMsg(msg, fname, -1) END IF ! Generation of random numbers always the same series during the whole program! CALL RANDOM_NUMBER(randv) ! Making sure that we do not repeat any value issampled = .FALSE. DO i=1, Nsample ! Generation of the index from the random numbers ind = MAX(INT(randv(i)*Nvals), 1) IF (.NOT.issampled(ind)) THEN sample(i) = ind issampled(ind) = .TRUE. ELSE ! Looking around the given index !PRINT *,' Index :', ind, ' already sampled!', issampled(ind) found = .FALSE. DO jmax=1, Nvals ind = MIN(ind+jmax, Nvals) IF (.NOT.issampled(ind)) THEN sample(i) = ind issampled(ind) = .TRUE. found = .TRUE. EXIT END IF ind = MAX(1, ind-jmax) IF (.NOT.issampled(ind)) THEN sample(i) = ind issampled(ind) = .TRUE. found = .TRUE. EXIT END IF END DO IF (.NOT.found) THEN msg = 'sampling could not be finished due to absence of available value!!' CALL ErrMsg(msg, fname, -1) END IF END IF END DO RETURN END SUBROUTINE rand_sample SUBROUTINE PrintQuantilesR_K(Nvals, vals, Nquants, qtvs, bspc) ! Subroutine to print the quantiles of values REAL(r_k) IMPLICIT NONE INTEGER, INTENT(in) :: Nvals, Nquants REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals REAL(r_k), DIMENSION(Nquants), INTENT(in) :: qtvs CHARACTER(len=1000), OPTIONAL :: bspc ! Local INTEGER :: iq LOGICAL, DIMENSION(Nvals) :: search1, search2, search CHARACTER(len=6) :: RS1 CHARACTER(len=50) :: fname CHARACTER(len=1000) :: bspcS !!!!!!! Variables ! vals: series of values ! qtvs: values of the quantiles ! bspc: base quantity of spaces fname = 'PrintQuantilesR_K' IF (PRESENT(bspc)) THEN bspcS = bspc ELSE bspcS = ' ' END IF DO iq=1, Nquants-1 WHERE (vals >= qtvs(iq)) search1 = .TRUE. ELSEWHERE search1 = .FALSE. END WHERE WHERE (vals < qtvs(iq+1)) search2 = .TRUE. ELSEWHERE search2 = .FALSE. END WHERE WHERE (search1 .AND. search2) search = .TRUE. ELSEWHERE search = .FALSE. END WHERE WRITE(RS1, '(F6.2)')(iq)*100./(Nquants-1) PRINT *, TRIM(bspcS) // '[',iq,']', TRIM(RS1) // ' %:', qtvs(iq), 'N:', COUNT(search) END DO RETURN END SUBROUTINE PrintQuantilesR_K INTEGER FUNCTION FindMinimumR_K(x, dsize, Startv, Endv) ! Function returns the location of the minimum in the section between Start and End. IMPLICIT NONE INTEGER, INTENT(in) :: dsize REAL(r_k), DIMENSION(dsize), INTENT(in) :: x INTEGER, INTENT(in) :: Startv, Endv ! Local REAL(r_k) :: Minimum INTEGER :: Location INTEGER :: i Minimum = x(Startv) ! assume the first is the min Location = Startv ! record its position DO i = Startv+1, Endv ! start with next elements IF (x(i) < Minimum) THEN ! if x(i) less than the min? Minimum = x(i) ! Yes, a new minimum found Location = i ! record its position END IF END DO FindMinimumR_K = Location ! return the position END FUNCTION FindMinimumR_K SUBROUTINE SwapR_K(a, b) ! Subroutine swaps the values of its two formal arguments. IMPLICIT NONE REAL(r_k), INTENT(INOUT) :: a, b ! Local REAL(r_k) :: Temp Temp = a a = b b = Temp END SUBROUTINE SwapR_K SUBROUTINE SortR_K(x, Nx) ! Subroutine receives an array x() r_K and sorts it into ascending order. IMPLICIT NONE INTEGER, INTENT(IN) :: Nx REAL(r_k), DIMENSION(Nx), INTENT(INOUT) :: x ! Local INTEGER :: i INTEGER :: Location DO i = 1, Nx-1 ! except for the last Location = FindMinimumR_K(x, Nx-i+1, i, Nx) ! find min from this to last CALL SwapR_K(x(i), x(Location)) ! swap this and the minimum END DO END SUBROUTINE SortR_K SUBROUTINE quantilesR_K(Nvals, vals, Nquants, quants) ! Subroutine to provide the quantiles of a given set of values of type real 'r_k' IMPLICIT NONE INTEGER, INTENT(in) :: Nvals, Nquants REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals REAL(r_k), DIMENSION(Nquants), INTENT(out) :: quants ! Local INTEGER :: i REAL(r_k) :: minv, maxv REAL(r_k), DIMENSION(Nvals) :: sortedvals !!!!!!! Variables ! Nvals: number of values ! Rk: kind of real of the values ! vals: values ! Nquants: number of quants ! quants: values at which the quantile start minv = MINVAL(vals) maxv = MAXVAL(vals) sortedvals = vals ! Using from: http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/chap08/sorting.f90 CALL SortR_K(sortedvals, Nvals) quants(1) = minv DO i=2, Nquants quants(i) = sortedvals(INT((i-1)*Nvals/Nquants)) END DO END SUBROUTINE quantilesR_K SUBROUTINE StatsR_K(Nvals, vals, minv, maxv, mean, mean2, stdev) ! Subroutine to provide the minmum, maximum, mean, the quadratic mean, and the standard deviation of a ! series of r_k numbers IMPLICIT NONE INTEGER, INTENT(in) :: Nvals REAL(r_k), DIMENSION(Nvals), INTENT(in) :: vals REAL(r_k), INTENT(out) :: minv, maxv, mean, mean2, stdev !!!!!!! Variables ! Nvals: number of values ! vals: values ! minv: minimum value of values ! maxv: maximum value of values ! mean: mean value of values ! mean2: quadratic mean value of values ! stdev: standard deviation of values minv = MINVAL(vals) maxv = MAXVAL(vals) mean=SUM(vals) mean2=SUM(vals*vals) mean=mean/Nvals mean2=mean2/Nvals stdev=SQRT(mean2 - mean*mean) RETURN END SUBROUTINE StatsR_k SUBROUTINE NcountR(values, d1, Ndiffvals, counts) ! Subroutine to count real values IMPLICIT NONE INTEGER, INTENT(in) :: d1 REAL(r_k), DIMENSION(d1), INTENT(in) :: values INTEGER, INTENT(out) :: Ndiffvals REAL(r_k), DIMENSION(d1,2), INTENT(out) :: counts ! Local INTEGER :: i, ival REAL(r_k), DIMENSION(d1) :: diffv !!!!!!! Variables ! values: values to count ! counts: counts of time for each value fname = 'NcountR' counts = -1. counts(1,1) = values(1) counts(1,2) = 1 Ndiffvals = 1 DO i=2,d1 diffv(1:Ndiffvals) = counts(1:Ndiffvals,1) - values(i) IF (ANY(diffv(1:Ndiffvals) == 0)) THEN ival = Index1DArrayR(counts(1:Ndiffvals,1), Ndiffvals, values(i)) counts(ival,2) = counts(ival,2) + 1 ELSE Ndiffvals = Ndiffvals + 1 counts(Ndiffvals,1) = values(i) counts(Ndiffvals,2) = 1 END IF END DO END SUBROUTINE NcountR SUBROUTINE runmean_F1D(d1, values, Nmean, headertail, runmean) ! Subroutine fo computing the running mean of a given set of float 1D values IMPLICIT NONE INTEGER, INTENT(in) :: d1, Nmean REAL(r_k), DIMENSION(d1), INTENT(in) :: values CHARACTER(len=*), INTENT(in) :: headertail REAL(r_k), DIMENSION(d1), INTENT(out) :: runmean ! Local INTEGER :: i, j, Nmean2 CHARACTER(len=5) :: NmeanS !!!!!!! Variables ! values: values to compute the running mean ! Nmean: number of odd points to use for the running mean ! headertail: How to proceed for the grid points at the beginning of the values which are not ! encompassed by the Nmean ! 'miss': set as missing values (1.d20) ! 'original': keep the original values ! 'progressfill': mean the values as a progressive running filter (e.g. for Nmean=5): ! runmean[values(1)] = values(1) ! runmean[values(2)] = MEAN(values(1:3)) ! runmean[values(3)] = MEAN(values(1:5)) ! runmean[values(4)] = MEAN(values(2:6)) ! (...) ! runmean[values(d1-2)] = MEAN(values(d1-5:d1)) ! runmean[values(d1-1)] = MEAN(values(d1-2:d1)) ! runmean[values(d1)] = MEAN(values(dd1)) ! 'zero': set as zero values ! runmean: runnig mean values fname = 'runmean_F1D' IF (MOD(Nmean,2) == 0) THEN WRITE(NmeanS,'(I5)')Nmean msg="Nmean has to be odd!! value provided: "// NmeanS CALL ErrMsg(msg, fname, -1) END IF Nmean2 = Nmean/2 SELECT CASE (TRIM(headertail)) CASE ('missing') runmean = fillval64 CASE ('original') runmean = values CASE ('progressfill') DO i=1, Nmean2 runmean(i) = SUM(values(1:2*(i-1)+1))/(2*(i-1)+1) END DO runmean(d1) = values(d1) DO i=2, Nmean2 j = d1-(2*(i-1)) runmean(d1-(i-1)) = SUM(values(j:d1))/(2*(i-1)+1) END DO CASE ('zero') runmean = zeroRK CASE DEFAULT msg = "'" // TRIM(headertail) // "' not available !!" //CHAR(44) // " available ones: " // & "'missing', 'original', 'progressfill', 'zero'" CALL ErrMsg(msg, fname, -1) END SELECT DO i= 1+Nmean2, d1 - Nmean2 runmean(i) = SUM(values(i-Nmean2:i+Nmean2))/Nmean END DO END SUBROUTINE runmean_F1D SUBROUTINE percentiles_R_K2D(values, axisS, Npercen, d1, d2, percentiles) ! Subroutine to compute the percentiles of a 2D R_K array along given set of axis IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, Npercen REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: values CHARACTER(LEN=*), INTENT(in) :: axisS REAL(r_k), DIMENSION(d1, d2, Npercen), INTENT(out) :: percentiles ! Local INTEGER :: i INTEGER :: Lstring, LaxisS, iichar CHARACTER(LEN=1000) :: splitaxis INTEGER, DIMENSION(1) :: axis1 CHARACTER(LEN=200), DIMENSION(2) :: axis2S INTEGER, DIMENSION(2) :: axis2 CHARACTER(LEN=1) :: Naxs !!!!!!! Variables ! d1,d2: length of the 2D dimensions ! values: values to use to compute the percentiles ! axisS: ':' separated list of axis to use to compute the percentiles ('all' for all axes) ! Npercen: number of percentiles ! percentiles: percentiles of the daata fname = 'percentiles_R_K2D' LaxisS = LEN_TRIM(axisS) iichar = numberTimes(axisS(1:LaxisS), ':') splitaxis = '' splitaxis(1:LaxisS) = axisS(1:LaxisS) percentiles = 0. IF (iichar == 0) THEN READ(axisS,'(I1)')axis1(1) ELSE IF (iichar == 1) THEN CALL split(splitaxis, ':', 2, axis2S) ELSE WRITE(Naxs,'(A1)')iichar msg = "' rank 2 values can not compute percentiles using " // Naxs // "' number of axis !!" CALL ErrMsg(msg, fname, -1) END IF IF (TRIM(axisS) == 'all') iichar = 2 IF (iichar == 0) THEN ! Might be a better way, but today I can't think it !! IF (axis1(1) == 1) THEN DO i=1, d2 CALL quantilesR_K(d1, values(:,i), Npercen, percentiles(1,i,:)) END DO ELSE IF (axis1(1) == 2) THEN DO i=1, d1 CALL quantilesR_K(d2, values(i,:), Npercen, percentiles(i,1,:)) END DO ELSE WRITE(Naxs,'(A1)')axis1(1) msg = "' rank 2 values can not compute percentiles using axis " // Naxs // "' !!" CALL ErrMsg(msg, fname, -1) END IF ELSE CALL quantilesR_K(d1*d2, RESHAPE(values, (/d1*d2/)), Npercen, percentiles(1,1,:)) END IF END SUBROUTINE percentiles_R_K2D SUBROUTINE percentiles_R_K3D(values, axisS, Npercen, d1, d2, d3, percentiles) ! Subroutine to compute the percentiles of a 3D R_K array along given set of axis IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, Npercen REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: values CHARACTER(LEN=*), INTENT(in) :: axisS REAL(r_k), DIMENSION(d1, d2, d3, Npercen), INTENT(out) :: percentiles ! Local INTEGER :: i, j INTEGER :: Lstring, LaxisS, iichar CHARACTER(LEN=1000) :: splitaxis INTEGER, DIMENSION(1) :: axis1 CHARACTER(LEN=200), DIMENSION(2) :: axis2S INTEGER, DIMENSION(2) :: axis2 CHARACTER(LEN=200), DIMENSION(3) :: axis3S INTEGER, DIMENSION(3) :: axis3 CHARACTER(LEN=1) :: Naxs1, Naxs2 !!!!!!! Variables ! d1,d2: length of the 2D dimensions ! values: values to use to compute the percentiles ! axisS: ':' separated list of axis to use to compute the percentiles ('all' for all axes) ! Npercen: number of percentiles ! percentiles: percentiles of the daata fname = 'percentiles_R_K3D' LaxisS = LEN_TRIM(axisS) iichar = numberTimes(axisS(1:LaxisS), ':') splitaxis = '' splitaxis(1:LaxisS) = axisS(1:LaxisS) percentiles = 0. IF (iichar == 0) THEN READ(axisS,'(I1)')axis1(1) ELSE IF (iichar == 1) THEN CALL split(splitaxis, ':', 2, axis2S) DO i=1,2 READ(axis2S(i), '(I1)')axis2(i) END DO ELSE IF (iichar == 2) THEN CALL split(splitaxis, ':', 3, axis3S) ELSE READ(Naxs1,'(A1)')iichar msg = "' rank 3 values can not compute percentiles using " // Naxs1 // "' number of axis !!" CALL ErrMsg(msg, fname, -1) END IF IF (TRIM(axisS) == 'all') iichar = 3 IF (iichar == 0) THEN ! Might be a better way, but today I can't think it !! IF (axis1(1) == 1) THEN DO i=1, d2 DO j=1, d3 CALL quantilesR_K(d1, values(:,i,j), Npercen, percentiles(1,i,j,:)) END DO END DO ELSE IF (axis1(1) == 2) THEN DO i=1, d1 DO j=1, d3 CALL quantilesR_K(d2, values(i,:,j), Npercen, percentiles(i,1,j,:)) END DO END DO ELSE IF (axis1(1) == 3) THEN DO i=1, d1 DO j=1, d2 CALL quantilesR_K(d3, values(i,j,:), Npercen, percentiles(i,j,1,:)) END DO END DO ELSE WRITE(Naxs1,'(A1)')axis1(1) msg = "' rank 3 values can not compute percentiles using axis " // Naxs1 // "' !!" CALL ErrMsg(msg, fname, -1) END IF ELSE IF (iichar == 1) THEN ! Might be a better way, but today I can't think it !! IF (axis2(1) == 1 .AND. axis2(2) == 2) THEN DO i=1, d3 CALL quantilesR_K(d1*d2, RESHAPE(values(:,:,i), (/d1*d2/)), Npercen, percentiles(1,1,i,:)) END DO ELSE IF (axis2(1) == 1 .AND. axis2(2) == 3) THEN DO i=1, d2 CALL quantilesR_K(d1*d3, RESHAPE(values(:,i,:), (/d1*d3/)), Npercen, percentiles(1,i,1,:)) END DO ELSE IF (axis2(1) == 2 .AND. axis2(2) == 3) THEN DO i=1, d1 CALL quantilesR_K(d2*d3, RESHAPE(values(i,:,:), (/d2*d3/)), Npercen, percentiles(i,1,1,:)) END DO ELSE WRITE(Naxs1,'(A1)')axis2(1) WRITE(Naxs2,'(A1)')axis2(2) msg="' rank 3 values can not compute percentiles using axis "//Naxs1// ', ' // Naxs2 // "' !!" CALL ErrMsg(msg, fname, -1) END IF ELSE CALL quantilesR_K(d1*d2*d3, RESHAPE(values, (/d1*d2*d3/)), Npercen, percentiles(1,1,1,:)) END IF END SUBROUTINE percentiles_R_K3D SUBROUTINE percentiles_R_K4D(values, axisS, Npercen, d1, d2, d3, d4, percentiles) ! Subroutine to compute the percentiles of a 4D R_K array along given set of axis IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d4, Npercen REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: values CHARACTER(LEN=*), INTENT(in) :: axisS REAL(r_k), DIMENSION(d1,d2,d3,d4,Npercen), INTENT(out) :: percentiles ! Local INTEGER :: i, j, k INTEGER :: Lstring, LaxisS, iichar CHARACTER(LEN=1000) :: splitaxis INTEGER, DIMENSION(1) :: axis1 CHARACTER(LEN=200), DIMENSION(2) :: axis2S INTEGER, DIMENSION(2) :: axis2 CHARACTER(LEN=200), DIMENSION(3) :: axis3S INTEGER, DIMENSION(3) :: axis3 CHARACTER(LEN=200), DIMENSION(4) :: axis4S CHARACTER(LEN=1) :: Naxs1, Naxs2, Naxs3 !!!!!!! Variables ! d1,d2: length of the 2D dimensions ! values: values to use to compute the percentiles ! axisS: ':' separated list of axis to use to compute the percentiles ('all' for all axes) ! Npercen: number of percentiles ! percentiles: percentiles of the daata fname = 'percentiles_R_K3D' LaxisS = LEN_TRIM(axisS) iichar = numberTimes(axisS(1:LaxisS), ':') splitaxis = '' splitaxis(1:LaxisS) = axisS(1:LaxisS) percentiles = 0. PRINT *,'iichar:', iichar, axisS(1:LaxisS) IF (iichar == 0) THEN READ(axisS,'(I1)')axis1(1) ELSE IF (iichar == 1) THEN CALL split(splitaxis, ':', 2, axis2S) DO i=1,2 READ(axis2S(i), '(I1)')axis2(i) END DO ELSE IF (iichar == 2) THEN CALL split(splitaxis, ':', 3, axis3S) DO i=1,3 READ(axis3S(i), '(I1)')axis3(i) END DO ELSE IF (iichar == 3) THEN CALL split(splitaxis, ':', 4, axis4S) ELSE READ(Naxs1,'(A1)')iichar msg = "' rank 4 values can not compute percentiles using " // Naxs1 // "' number of axis !!" CALL ErrMsg(msg, fname, -1) END IF IF (TRIM(axisS) == 'all') iichar = 4 IF (iichar == 0) THEN ! Might be a better way, but today I can't think it !! IF (axis1(1) == 1) THEN DO i=1, d2 DO j=1, d3 DO k=1, d4 CALL quantilesR_K(d1, values(:,i,j,k), Npercen, percentiles(1,i,j,k,:)) END DO END DO END DO ELSE IF (axis1(1) == 2) THEN DO i=1, d1 DO j=1, d3 DO k=1, d4 CALL quantilesR_K(d2, values(i,:,j,k), Npercen, percentiles(i,1,j,k,:)) END DO END DO END DO ELSE IF (axis1(1) == 3) THEN DO i=1, d1 DO j=1, d2 DO k=1, d4 CALL quantilesR_K(d3, values(i,j,:,k), Npercen, percentiles(i,j,1,k,:)) END DO END DO END DO ELSE IF (axis1(1) == 4) THEN DO i=1, d1 DO j=1, d2 DO k=1, d3 CALL quantilesR_K(d4, values(i,j,k,:), Npercen, percentiles(i,j,k,1,:)) END DO END DO END DO ELSE WRITE(Naxs1,'(A1)')axis1(1) msg = "' rank 3 values can not compute percentiles using axis " // Naxs1 // "' !!" CALL ErrMsg(msg, fname, -1) END IF ELSE IF (iichar == 1) THEN ! Might be a better way, but today I can't think it !! IF (axis2(1) == 1 .AND. axis2(2) == 2) THEN DO i=1, d3 DO j=1, d4 CALL quantilesR_K(d1*d2, RESHAPE(values(:,:,i,j), (/d1*d2/)), Npercen, & percentiles(1,1,i,j,:)) END DO END DO ELSE IF (axis2(1) == 1 .AND. axis2(2) == 3) THEN DO i=1, d2 DO j=1, d4 CALL quantilesR_K(d1*d3, RESHAPE(values(:,i,:,j), (/d1*d3/)), Npercen, & percentiles(1,i,1,j,:)) END DO END DO ELSE IF (axis2(1) == 1 .AND. axis2(2) == 4) THEN DO i=1, d2 DO j=1, d3 CALL quantilesR_K(d1*d4, RESHAPE(values(:,i,j,:), (/d1*d4/)), Npercen, & percentiles(1,i,j,1,:)) END DO END DO ELSE IF (axis2(1) == 2 .AND. axis2(2) == 3) THEN DO i=1, d1 DO j=1, d4 CALL quantilesR_K(d2*d3, RESHAPE(values(i,:,:,j), (/d2*d3/)), Npercen, & percentiles(i,1,1,j,:)) END DO END DO ELSE IF (axis2(1) == 2 .AND. axis2(2) == 4) THEN DO i=1, d1 DO j=1, d3 CALL quantilesR_K(d2*d4, RESHAPE(values(i,:,j,:), (/d2*d4/)), Npercen, & percentiles(i,1,j,1,:)) END DO END DO ELSE IF (axis2(1) == 3 .AND. axis2(2) == 4) THEN DO i=1, d1 DO j=1, d2 CALL quantilesR_K(d3*d4, RESHAPE(values(i,j,:,:), (/d3*d4/)), Npercen, & percentiles(i,j,1,1,:)) END DO END DO ELSE WRITE(Naxs1,'(A1)')axis2(1) WRITE(Naxs2,'(A1)')axis2(2) msg="' rank 4 values can not compute percentiles using axis "//Naxs1// ', ' // Naxs2 // "' !!" CALL ErrMsg(msg, fname, -1) END IF ELSE IF (iichar == 2) THEN IF (axis2(1) == 1 .AND. axis2(2) == 2 .AND. axis3(3) == 3) THEN DO i=1, d4 CALL quantilesR_K(d1*d2*d3, RESHAPE(values(:,:,:,i), (/d1*d2*d3/)), Npercen, & percentiles(1,1,1,i,:)) END DO ELSE IF (axis2(1) == 1 .AND. axis2(2) == 2 .AND. axis3(3) == 4) THEN DO i=1, d3 CALL quantilesR_K(d1*d2*d4, RESHAPE(values(:,:,i,:), (/d1*d2*d4/)), Npercen, & percentiles(1,1,i,1,:)) END DO ELSE IF (axis2(1) == 1 .AND. axis2(2) == 3 .AND. axis3(3) == 4) THEN DO i=1, d2 CALL quantilesR_K(d1*d3*d4, RESHAPE(values(:,i,:,:), (/d1*d3*d4/)), Npercen, & percentiles(1,i,1,1,:)) END DO ELSE IF (axis2(1) == 2 .AND. axis2(2) == 3 .AND. axis3(3) == 4) THEN DO i=1, d1 CALL quantilesR_K(d2*d3*d4, RESHAPE(values(i,:,:,:), (/d2*d3*d4/)), Npercen, & percentiles(i,1,1,1,:)) END DO ELSE WRITE(Naxs1,'(A1)')axis3(1) WRITE(Naxs2,'(A1)')axis3(2) WRITE(Naxs3,'(A1)')axis3(2) msg="' rank 4 values can not compute percentiles using axis "// Naxs1 // ', ' // Naxs2 // & ', ' // Naxs3 //"' !!" CALL ErrMsg(msg, fname, -1) END IF ELSE CALL quantilesR_K(d1*d2*d3*d4, RESHAPE(values, (/d1*d2*d3*d4/)), Npercen, percentiles(1,1,1,1,:)) END IF END SUBROUTINE percentiles_R_K4D REAL(r_k) FUNCTION distanceRK(pointA, pointB) ! Function to provide the distance between two points IMPLICIT NONE REAL(r_k), DIMENSION(2), INTENT(in) :: pointA, pointB !!!!!!! Variables ! pointA, B: couple of points to compute the distance between them fname = 'distanceRK' distanceRK = SQRT( (pointB(1)-pointA(1))**2 + (pointB(2)-pointA(2))**2 ) END FUNCTION distanceRK REAL(r_k) FUNCTION shoelace_area_polygon(Nvertex, poly) ! Computing the area of a polygon using sholace formula ! FROM: https://en.wikipedia.org/wiki/Shoelace_formula IMPLICIT NONE INTEGER, INTENT(in) :: Nvertex REAL(r_k), DIMENSION(Nvertex,2), INTENT(in) :: poly ! Local INTEGER :: i REAL(r_k) :: areapos, areaneg !!!!!!! Variables ! Nvertex: number of vertices of the polygon ! poly: coordinates of the vertex of the polygon (sorted) fname = 'shoelace_area_polygon' areapos = 0. areaneg = 0. DO i=1, Nvertex-1 areapos = areapos + poly(i,1)*poly(i+1,2) areaneg = areaneg + poly(i+1,1)*poly(i,2) END DO areapos = areapos + poly(Nvertex,1)*poly(1,2) areaneg = areaneg + poly(1,1)*poly(Nvertex,2) shoelace_area_polygon = 0.5*(areapos - areaneg) END FUNCTION shoelace_area_polygon SUBROUTINE intersection_2Dlines(lineA, lineB, intersect, ptintersect) ! Subroutine to provide the intersection point between two lines on the plane using Cramer's method IMPLICIT NONE REAL(r_k), DIMENSION(2,2), INTENT(in) :: lineA, lineB LOGICAL, INTENT(out) :: intersect REAL(r_k), DIMENSION(2), INTENT(out) :: ptintersect ! Local REAL(r_k), DIMENSION(2) :: segmentA, segmentB REAL(r_k) :: a11, a12, a21, a22, z1, z2 REAL(r_k) :: det, detX, detY LOGICAL :: axisAx, axisBx, axisAy, axisBy !!!!!!! Variables ! lineA: couple of coordinates for the line A ! lineB: couple of coordinates for the line B ! intersect: whether two lines intersect ! ptintersect: point of intersection [(0,0) if they do not intersect] fname = 'intersection_2Dlines' axisAx = .FALSE. axisAy = .FALSE. axisBx = .FALSE. axisBy = .FALSE. ! Setting segment parameters y = A + B*x IF (lineA(2,1) /= lineA(1,1)) THEN segmentA(2) = (lineA(2,2)-lineA(1,2))/(lineA(2,1)-lineA(1,1)) ! This might be to ask too much... ? !IF ( (lineA(1,1)*segmentA(2) - lineA(1,2)) /= (lineA(2,1)*segmentA(2) - lineA(2,2)) ) THEN ! PRINT *,'A = y1 - x2*B = ', lineA(1,2) - lineA(1,1)*segmentA(2) ! PRINT *,'A = y2 - x2*B = ', lineA(2,2) - lineA(2,1)*segmentA(2) ! msg = 'Wrong calculation of parameter A, for lineA' ! CALL ErrMSg(msg, fname, -1) !END IF segmentA(1) = lineA(1,2) - lineA(1,1)*segmentA(2) a11 = segmentA(2) a12 = -oneRK z1 = -segmentA(1) IF (lineA(2,2) == lineA(1,2)) axisAx = .TRUE. ELSE ! lineA || y-axis axisAy = .TRUE. END IF IF (lineB(2,1) /= lineB(1,1)) THEN segmentB(2) = (lineB(2,2)-lineB(1,2))/(lineB(2,1)-lineB(1,1)) ! This might be to ask too much... ? !IF ( (lineB(1,1)*segmentB(2) - lineB(1,2)) /= (lineB(2,1)*segmentB(2) - lineB(2,2)) ) THEN ! PRINT *,'A = x1*B -y1 = ', lineB(1,1)*segmentB(2) - lineB(1,2) ! PRINT *,'A = x2*B -y2 = ', lineB(2,1)*segmentB(2) - lineB(2,2) ! msg = 'Wrong calculation of parameter A, for lineB' ! CALL ErrMSg(msg, fname, -1) !END IF segmentB(1) = lineB(1,2) - lineB(1,1)*segmentB(2) a21 = segmentB(2) a22 = -oneRK z2 = -segmentB(1) IF (lineB(2,2) == lineB(1,2)) axisBx = .TRUE. ELSE ! lineB || y-axis axisBy = .TRUE. END IF ! Cramer's method ! a11 = B1; a12 = -1 ! a21 = B2; a22 = -1 ! z1 = -A1 ! z2 = -A2 ! (a11 a12)(x) (z1) ! (a21 a22)(y) (z2) ! -------- ------ ----- ---- --- -- - ! det = (a11*a22-a12*a21) ! detX = (z1*a22-z2*a21) ! detY = (a11*z1-a12*z2) ! ptintercept = (detX/det, detY/det) ! Cases when some of the lines are parallel to any given axis ! PRINT *,' axisAx', axisAx, 'axisAy', axisAy, 'axisBx', axisBx, 'axisBy', axisBy IF (axisAx .OR. axisAy .OR. axisBx .OR. axisBy) THEN IF (axisAx) THEN IF (axisBy) THEN intersect = .TRUE. ptintersect(1) = lineB(1,1) ptintersect(2) = lineA(1,2) ELSE intersect = .TRUE. ptintersect(1) = (lineA(1,2)-segmentB(1))/segmentB(2) ptintersect(2) = lineA(1,2) END IF END IF IF (axisAy) THEN IF (axisBy) THEN intersect = .TRUE. ptintersect(1) = lineA(1,1) ptintersect(2) = lineB(1,2) ELSE intersect = .TRUE. ptintersect(1) = lineA(1,1) ptintersect(2) = segmentB(1) + lineA(1,1)*segmentB(2) END IF END IF IF (axisBx) THEN IF (axisAy) THEN intersect = .TRUE. ptintersect(1) = lineA(1,1) ptintersect(2) = lineB(1,2) ELSE intersect = .TRUE. ptintersect(1) = (lineB(1,2)-segmentA(1))/segmentA(2) ptintersect(2) = lineB(1,2) END IF END IF IF (axisBy) THEN IF (axisAx) THEN intersect = .TRUE. ptintersect(1) = lineB(1,1) ptintersect(2) = lineA(1,2) ELSE intersect = .TRUE. ptintersect(1) = lineB(1,1) ptintersect(2) = segmentA(1) + lineB(1,1)*segmentA(2) END IF END IF ELSE det = (a11*a22-a12*a21) ! Parallel lines ! IF (det == zeroRK) THEN intersect = .FALSE. ptintersect = zeroRK ELSE intersect = .TRUE. detX = (z1*a22-z2*a12) detY = (a11*z2-a21*z1) ptintersect(1) = detX/det ptintersect(2) = detY/det END IF END IF END SUBROUTINE intersection_2Dlines !refs: !https://www.mathopenref.com/heronsformula.html !https://math.stackexchange.com/questions/1406340/intersect-area-of-two-polygons-in-cartesian-plan !http://www.cap-lore.com/MathPhys/IP/ !http://www.cap-lore.com/MathPhys/IP/IL.html !https://www.element84.com/blog/determining-the-winding-of-a-polygon-given-as-a-set-of-ordered-points !https://stackoverflow.com/questions/1165647/how-to-determine-if-a-list-of-polygon-points-are-in-clockwise-order !https://en.wikipedia.org/wiki/Shoelace_formula !https://en.wikipedia.org/wiki/Winding_number !https://en.wikipedia.org/wiki/Simple_polygon !https://en.wikipedia.org/wiki/Polygon#Properties !https://en.wikipedia.org/wiki/Convex_polygon !https://en.wikipedia.org/wiki/Jordan_curve_theorem !https://www.sangakoo.com/ca/temes/metode-de-cramer !https://www.geogebra.org/m/pw4QHFYT SUBROUTINE intersectfaces(faceA, faceB, intersect, intersectpt) ! Subroutine to provide if two faces of two polygons intersect ! AFTER: http://www.cap-lore.com/MathPhys/IP/IL.html ! A: faceA(1,:) ! B: faceA(2,:) ! C: faceB(1,:) ! D: faceB(2,:) IMPLICIT NONE REAL(r_k), DIMENSION(2,2), INTENT(in) :: faceA, faceB INTEGER, INTENT(out) :: intersect REAL(r_k), DIMENSION(2), INTENT(out) :: intersectpt ! Local REAL(r_k) :: Axmin, Aymin, Axmax, Aymax REAL(r_k) :: Bxmin, Bymin, Bxmax, Bymax REAL(r_k) :: areaABD, areaACD, areaBDC, areaDAB REAL(r_k), DIMENSION(3,2) :: triangle LOGICAL :: Lintersect !!!!!!! Variables ! faceA/B: coordinates of faces A and B to determine if they intersect ! intersect: integer to say if they intersect (0, no-intersect, +/-1 intersect) ! intersectpt: point where faces intersect [(0,0) otherwise] fname = 'intersectfaces' ! PRINT *,' ' // TRIM(fname) // ' ________' ! PRINT *,' faceA:', faceA(1,:), ';',faceA(2,:) ! PRINT *,' faceB:', faceB(1,:), ';',faceB(2,:) Axmin = MINVAL(faceA(:,1)) Axmax = MAXVAL(faceA(:,1)) Aymin = MINVAL(faceA(:,2)) Aymax = MAXVAL(faceA(:,2)) Bxmin = MINVAL(faceB(:,1)) Bxmax = MAXVAL(faceB(:,1)) Bymin = MINVAL(faceB(:,2)) Bymax = MAXVAL(faceB(:,2)) ! No intersection IF ( (Axmax <= Bxmin) .OR. (Axmin >= Bxmax) .OR. (Aymax <= Bymin) .OR. (Aymin >= Bymax) ) THEN intersect = 0 intersectpt = zeroRK ELSE ! Triangle ABD triangle(1,:) = faceA(1,:) triangle(2,:) = faceA(2,:) triangle(3,:) = faceB(2,:) areaABD = shoelace_area_polygon(3, triangle) ! Triangle ACD triangle(1,:) = faceA(1,:) triangle(2,:) = faceB(1,:) triangle(3,:) = faceB(2,:) areaACD = shoelace_area_polygon(3, triangle) ! Triangle BDC triangle(1,:) = faceA(2,:) triangle(2,:) = faceB(2,:) triangle(3,:) = faceB(1,:) areaBDC = shoelace_area_polygon(3, triangle) ! Triangle DAB triangle(1,:) = faceB(2,:) triangle(2,:) = faceA(1,:) triangle(3,:) = faceA(2,:) areaDAB = shoelace_area_polygon(3, triangle) IF (areaABD>zeroRK .AND. areaACD>zeroRK .AND. areaBDC>zeroRK .AND. areaDAB>zeroRK) THEN intersect = INT(ABS(areaABD)/areaABD) CALL intersection_2Dlines(faceA, faceB, Lintersect, intersectpt) ELSE IF (areaABD zeroRK, 'I', intersect END IF END SUBROUTINE intersectfaces LOGICAL FUNCTION poly_has_point(Nvertex, polygon, point) ! Function to determine if a polygon has already a given point as one of its vertex IMPLICIT NONE INTEGER, INTENT(in) :: Nvertex REAL(r_k), DIMENSION(Nvertex,2), INTENT(in) :: polygon REAL(r_k), DIMENSION(2), INTENT(in) :: point ! Local INTEGER :: iv REAL(r_k), DIMENSION(2) :: diff !!!!!!! Vertrex ! Nvertex: number of vertexs of the polygon ! polygon: vertexs of the polygon ! point: point to look for its ownership into the polygon fname = 'poly_has_point' poly_has_point = .FALSE. DO iv=1, Nvertex diff = polygon(iv,:)-point IF ( (diff(1) == zeroRK) .AND. (diff(2) == zeroRK)) THEN poly_has_point = .TRUE. EXIT END IF END DO END FUNCTION poly_has_point SUBROUTINE join_polygon(NvertexA, NvertexB, NvertexAB, polyA, polyB, Ncoinvertex, coinpoly) ! Subroutine to join two polygons ! AFTER: http://www.cap-lore.com/MathPhys/IP/ and http://www.cap-lore.com/MathPhys/IP/IL.html IMPLICIT NONE INTEGER, INTENT(in) :: NvertexA, NvertexB, NvertexAB REAL(r_k), DIMENSION(NvertexA,2), INTENT(in) :: polyA REAL(r_k), DIMENSION(NvertexB,2), INTENT(in) :: polyB INTEGER, INTENT(out) :: Ncoinvertex REAL(r_k), DIMENSION(NvertexAB,2), INTENT(out) :: coinpoly ! Local INTEGER :: iA, iB, icoin, ii REAL(r_k), DIMENSION(2,2) :: face1, face2 INTEGER :: intersct REAL(r_k), DIMENSION(2) :: ptintersct !!!!!!! variables ! NvertexA: number of vertexs polygon A ! NvertexB: number of vertexs polygon B ! polyA: pairs of coordinates for the polygon A (clockwise) ! polyB: pairs of coordinates for the polygon B (clockwise) ! Ncoinvertex: number of vertexes for the coincident polygon ! coinpoly: pairs of coordinates for the coincident polygon (clockwise) fname = 'join_polygon' icoin = 0 coinpoly = 0. ! First, include that vertex which do not lay within any polygon DO iA=1, NvertexA !PRINT *, ' iA:', iA, ':', polyA(iA,:), point_inside(polyA(iA,:), NvertexB, polyB) IF (.NOT. point_inside(polyA(iA,:), NvertexB, polyB)) THEN icoin = icoin + 1 coinpoly(icoin,:) = polyA(iA,:) END IF END DO DO iB=1, NvertexB !PRINT *, ' iB:', iB, ':', polyB(iB,:), point_inside(polyB(iB,:), NvertexA, polyA) IF (.NOT. point_inside(polyB(iB,:), NvertexA, polyA)) THEN icoin = icoin + 1 coinpoly(icoin,:) = polyB(iB,:) END IF END DO DO iA=1, NvertexA ! Getting couple of vertexs from polyA and polyB IF (iA /= NvertexA) THEN face1(1,:) = polyA(iA,:) face1(2,:) = polyA(iA+1,:) ELSE face1(1,:) = polyA(iA,:) face1(2,:) = polyA(1,:) END IF DO iB=1, NvertexB IF (iB /= NvertexB) THEN face2(1,:) = polyB(iB,:) face2(2,:) = polyB(iB+1,:) ELSE face2(1,:) = polyB(iB,:) face2(2,:) = polyB(1,:) END IF ! Compute areas of the four possible triangles. Introduce the coincident vertexs not included CALL intersectfaces(face1, face2, intersct, ptintersct) !PRINT *,iA,':',face1(1,:),';',face1(2,:), '=', iB, face2(1,:),';',face2(2,:), '<>', intersct,':', ptintersct IF (intersct == 1) THEN IF (.NOT.poly_has_point(icoin,coinpoly(1:icoin,:),ptintersct) ) THEN icoin = icoin + 1 coinpoly(icoin,:) = ptintersct END IF ELSE IF (intersct == -1) THEN IF (.NOT.poly_has_point(icoin,coinpoly(1:icoin,:),ptintersct) ) THEN icoin = icoin + 1 coinpoly(icoin,:) = ptintersct END IF END IF END DO END DO Ncoinvertex = icoin END SUBROUTINE join_polygon SUBROUTINE sort_polygon(Nvertex, polygon, sense, Nnewvertex, newpoly) ! Subroutine to sort a polygon using its center as average of the coordinates and remove duplicates ! Should be used the centroid instead, but by now let do it simple ! https://en.wikipedia.org/wiki/Centroid IMPLICIT NONE INTEGER, INTENT(in) :: Nvertex, sense REAL(r_k), DIMENSION(Nvertex,2), INTENT(in) :: polygon INTEGER, INTENT(out) :: Nnewvertex REAL(r_k), DIMENSION(Nvertex,2), INTENT(out) :: newpoly ! Local INTEGER :: iv, j REAL(r_k) :: vang REAL(r_k), DIMENSION(2) :: center REAL(r_k), DIMENSION(Nvertex) :: angles REAL(r_k), DIMENSION(Nvertex,2) :: sortpoly !!!!!!! Variables ! Nvertex: number of vertices ! polygon: coordinates of the vertices of the polygon ! sense: sens of sorting thepolygon (1: clockwise, -1: anti-clockwise) ! sortpoly: sorted polygon ! Nnewvertex: number of vertices new polygon ! newpoly: sorted and duplicate removed polygon fname = 'sort_polygon' ! To be substituted by centroid calculation (which requires already sorted vetexs...) center(1) = SUM(polygon(:,1))/Nvertex center(2) = SUM(polygon(:,2))/Nvertex DO iv=1, Nvertex angles(iv) = ATAN2(polygon(iv,2)-center(2),polygon(iv,1)-center(1)) END DO CALL sortR_K(angles, Nvertex) sortpoly = zeroRK DO iv=1, Nvertex DO j=1, Nvertex vang = ATAN2(polygon(j,2)-center(2), polygon(j,1)-center(1)) IF (angles(iv) == vang) THEN IF (sense == -1) THEN sortpoly(iv,:) = polygon(j,:) ELSE sortpoly(Nvertex-iv+1,:) = polygon(j,:) END IF EXIT END IF END DO END DO newpoly(1,:) = sortpoly(1,:) j = 1 DO iv=2, Nvertex IF (.NOT.poly_has_point(j,newpoly(1:j,:),sortpoly(iv,:)) ) THEN j = j+1 newpoly(j,:) = sortpoly(iv,:) END IF END DO Nnewvertex = j END SUBROUTINE sort_polygon LOGICAL FUNCTION point_inside(point, Nvertex, polygon) ! Function to determine if a given point is inside a polygon providing its sorted vertices ! FROM: https://en.wikipedia.org/wiki/Point_in_polygon IMPLICIT NONE REAL(r_k), DIMENSION(2), INTENT(in) :: point INTEGER, INTENT(in) :: Nvertex REAL(r_k), DIMENSION(Nvertex,2), INTENT(in) :: polygon ! Local INTEGER :: iv, Nintersect INTEGER :: cross REAL(r_k) :: xmin REAL(r_k), DIMENSION(2) :: crosspoint REAL(r_k), DIMENSION(2,2) :: face1, face2 REAL(r_k), DIMENSION(Nvertex) :: abovebelow !!!!!!! Variables ! point: point to look for ! Nvertrex: number of vertices of a polygon ! polygon: vertices of a polygon fname = 'point_inside' xmin = MINVAL(polygon(:,1)) ! Looking for the intersection with the ray Nintersect = 0 face1(1,:) = (/ xmin-0.5, point(2) /) face1(2,:) = (/ point(1), point(2) /) DO iv = 1, Nvertex IF (iv /= Nvertex) THEN face2(1,:) = polygon(iv,:) face2(2,:) = polygon(iv+1,:) ELSE face2(1,:) = polygon(iv,:) face2(2,:) = polygon(1,:) END IF CALL intersectfaces(face1, face2, cross, crosspoint) IF (cross /= 0) THEN Nintersect = Nintersect + 1 abovebelow(Nintersect) = iv END IF END DO IF (MOD(Nintersect,2) == 0) THEN point_inside = .FALSE. ELSE point_inside = .TRUE. END IF END FUNCTION point_inside LOGICAL FUNCTION point_in_face(pt, Nvertex, poly) ! Function to determine if a given point is on a face of a polygon IMPLICIT NONE REAL(r_k), DIMENSION(2), INTENT(in) :: pt INTEGER, INTENT(in) :: Nvertex REAL(r_k), DIMENSION(Nvertex,2), INTENT(in) :: poly ! Local INTEGER :: iv REAL(r_k) :: ix, ex, iy, ey, tmpv REAL(r_k) :: dx, dy, A, B !!!!!!! Variables ! pt: point to look for ! Nvertex: Number of vertices of the polygon ! poly: polygon fname = 'point_in_face' point_in_face = .FALSE. DO iv=1, Nvertex IF (iv < Nvertex) THEN ix = poly(iv,1) ex = poly(iv+1,1) iy = poly(iv,2) ey = poly(iv+1,2) ELSE ix = poly(iv,1) ex = poly(1,1) iy = poly(iv,2) ey = poly(1,2) END IF dx = ex - ix dy = ey - iy IF (dx == zeroRK) THEN IF (pt(1) == ix) THEN IF ( (iy < ey) .AND. (pt(2) >= iy) .AND. pt(2) <= ey) THEN point_in_face = .TRUE. EXIT ELSE IF ( (iy > ey) .AND. (pt(2) >= ey) .AND. pt(2) <= iy) THEN point_in_face = .TRUE. EXIT END IF END IF ELSE IF (dy == zeroRK) THEN IF (pt(2) == iy) THEN IF ((ix < ex) .AND. (pt(1) >= ix) .AND. pt(1) <= ex) THEN point_in_face = .TRUE. EXIT ELSE IF ((ix > ex) .AND. (pt(1) >= ex) .AND. pt(1) <= ix) THEN point_in_face = .TRUE. EXIT END IF END IF ELSE A = iy B = (ey-iy)/(ex-ix) IF (A+B*(pt(1)-ix) == pt(2)) THEN point_in_face = .TRUE. EXIT END IF END IF END IF END DO END FUNCTION point_in_face SUBROUTINE coincident_polygon(NvertexA, NvertexB, NvertexAB, polyA, polyB, Ncoinvertex, coinpoly) ! Subroutine to provide the intersection polygon between two polygons ! AFTER: http://www.cap-lore.com/MathPhys/IP/ and http://www.cap-lore.com/MathPhys/IP/IL.html IMPLICIT NONE INTEGER, INTENT(in) :: NvertexA, NvertexB, NvertexAB REAL(r_k), DIMENSION(NvertexA,2), INTENT(in) :: polyA REAL(r_k), DIMENSION(NvertexB,2), INTENT(in) :: polyB INTEGER, INTENT(out) :: Ncoinvertex REAL(r_k), DIMENSION(NvertexAB,2), INTENT(out) :: coinpoly ! Local INTEGER :: iA, iB, icoin, ii REAL(r_k), DIMENSION(2,2) :: face1, face2 INTEGER :: intersct REAL(r_k), DIMENSION(2) :: ptintersct !!!!!!! variables ! NvertexA: number of vertexs polygon A ! NvertexB: number of vertexs polygon B ! polyA: pairs of coordinates for the polygon A (clockwise) ! polyB: pairs of coordinates for the polygon B (clockwise) ! Ncoinvertex: number of vertexes for the coincident polygon ! coinpoly: pairs of coordinates for the coincident polygon (clockwise) fname = 'coincident_polygon' icoin = 0 coinpoly = 0. ! First, include that vertex which lay within any polygon DO iA=1, NvertexA IF (point_inside(polyA(iA,:), NvertexB, polyB)) THEN icoin = icoin + 1 coinpoly(icoin,:) = polyA(iA,:) END IF IF (point_in_face(polyA(iA,:), NvertexB, polyB)) THEN icoin = icoin + 1 coinpoly(icoin,:) = polyA(iA,:) END IF END DO DO iB=1, NvertexB IF (point_inside(polyB(iB,:), NvertexA, polyA)) THEN icoin = icoin + 1 coinpoly(icoin,:) = polyB(iB,:) END IF IF (point_in_face(polyB(iB,:), NvertexA, polyA)) THEN icoin = icoin + 1 coinpoly(icoin,:) = polyB(iB,:) END IF END DO ! Look interesections DO iA=1, NvertexA ! Getting couple of vertexs from polyA and polyB IF (iA /= NvertexA) THEN face1(1,:) = polyA(iA,:) face1(2,:) = polyA(iA+1,:) ELSE face1(1,:) = polyA(iA,:) face1(2,:) = polyA(1,:) END IF DO iB=1, NvertexB IF (iB /= NvertexB) THEN face2(1,:) = polyB(iB,:) face2(2,:) = polyB(iB+1,:) ELSE face2(1,:) = polyB(iB,:) face2(2,:) = polyB(1,:) END IF ! Compute areas of the four possible triangles. Introduce the coincident vertexs not included CALL intersectfaces(face1, face2, intersct, ptintersct) !PRINT *,iA,':',face1(1,:),';',face1(2,:), '=', iB, face2(1,:),';',face2(2,:), '<>', intersct,':', ptintersct IF ((intersct /= 0) .AND. (.NOT.poly_has_point(icoin,coinpoly(1:icoin,:),ptintersct)) ) THEN icoin = icoin + 1 coinpoly(icoin,:) = ptintersct END IF END DO END DO Ncoinvertex = icoin END SUBROUTINE coincident_polygon SUBROUTINE grid_within_polygon(NvertexA, polygonA, dx, dy, dxy, xCvals, yCvals, Nvertexmax, xBvals, & yBvals, Ngridsin, gridsin) ! Subroutine to determine which grid cells from a matrix lay inside a polygon IMPLICIT NONE INTEGER, INTENT(in) :: NvertexA, dx, dy, dxy, Nvertexmax REAL(r_k), DIMENSION(NvertexA,2), INTENT(in) :: polygonA REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: xCvals, yCvals REAL(r_k), DIMENSION(dx,dy,Nvertexmax), INTENT(in) :: xBvals, yBvals INTEGER, INTENT(out) :: Ngridsin INTEGER, DIMENSION(dxy,2), INTENT(out) :: gridsin ! Local INTEGER :: ix, iy, iv, ivA INTEGER :: cross REAL(r_k), DIMENSION(2) :: centergrid, vertex, crosspoint REAL(r_k), DIMENSION(2,2) :: face, faceA LOGICAL, DIMENSION(dx,dy) :: within REAL(r_k), DIMENSION(:,:), ALLOCATABLE :: gridpoly !!!!!!! Variables ! NvertexA: Number of vertices of the polygon to find the grids ! polygonA: ordered vertices of the polygon ! dx, dy: shape of the matrix with the grid points ! xCvals, yCvals: coordinates of the center of the grid cells ! Nvertexmax: Maximum number of vertices of the grid cells ! xBvals, yBvals: coordinates of th vertices of the grid cells (-99999 for no vertex) ! Ngridsin: number of grids with some extension within the polygon ! gridsin: grids within the polygin ! percentages: percentages of area of each of the grids within the polygon fname = 'grid_within_polygon' Ngridsin = 0 gridsin = 0 within = .FALSE. DO ix = 1, dx DO iy = 1, dy IF (.NOT.within(ix,iy)) THEN centergrid = (/ xCvals(ix,iy), yCvals(ix,iy) /) ! By grid center IF (point_inside(centergrid, NvertexA, polygonA)) THEN Ngridsin = Ngridsin + 1 ! Getting coordinates gridsin(Ngridsin,1) = ix gridsin(Ngridsin,2) = iy within(ix,iy) = .TRUE. CYCLE END IF ! Getting grid vertices DO iv=1, Nvertexmax IF (.NOT.within(ix,iy)) THEN IF (xBvals(ix,iy,iv) /= fillvalI) THEN vertex = (/ xBvals(ix,iy,iv), yBvals(ix,iy,iv) /) IF (point_inside(vertex, NvertexA, polygonA)) THEN Ngridsin = Ngridsin + 1 ! Getting coordinates gridsin(Ngridsin,1) = ix gridsin(Ngridsin,2) = iy within(ix,iy) = .TRUE. CYCLE END IF END IF END IF END DO ! Getting grid faces DO iv=1, Nvertexmax-1 IF (.NOT.within(ix,iy)) THEN IF (xBvals(ix,iy,iv) /= fillvalI) THEN face(1,:) = (/ xBvals(ix,iy,iv), yBvals(ix,iy,iv) /) face(2,:) = (/ xBvals(ix,iy,iv+1), yBvals(ix,iy,iv+1) /) DO ivA=1, NvertexA-1 faceA(1,:) = (/ polygonA(ivA,1), polygonA(ivA,2) /) faceA(2,:) = (/polygonA(ivA+1,1), polygonA(ivA+1,2) /) CALL intersectfaces(face, faceA, cross, crosspoint) IF (cross /= 0) THEN Ngridsin = Ngridsin + 1 ! Getting coordinates gridsin(Ngridsin,1) = ix gridsin(Ngridsin,2) = iy within(ix,iy) = .TRUE. CYCLE END IF END DO IF (.NOT.within(ix,iy)) THEN faceA(1,:) = (/ polygonA(NvertexA,1), polygonA(NvertexA,2) /) faceA(2,:) = (/polygonA(1,1), polygonA(1,2) /) CALL intersectfaces(face, faceA, cross, crosspoint) IF (cross /= 0) THEN Ngridsin = Ngridsin + 1 ! Getting coordinates gridsin(Ngridsin,1) = ix gridsin(Ngridsin,2) = iy within(ix,iy) = .TRUE. END IF END IF END IF END IF END DO iv = Nvertexmax IF (.NOT.within(ix,iy)) THEN IF (xBvals(ix,iy,iv) /= fillvalI) THEN face(1,:) = (/ xBvals(ix,iy,iv), yBvals(ix,iy,iv) /) face(2,:) = (/ xBvals(ix,iy,1), yBvals(ix,iy,1) /) DO ivA=1, NvertexA-1 faceA(1,:) = (/ polygonA(ivA,1), polygonA(ivA,2) /) faceA(2,:) = (/polygonA(ivA+1,1), polygonA(ivA+1,2) /) CALL intersectfaces(face, faceA, cross, crosspoint) IF (cross /= 0) THEN Ngridsin = Ngridsin + 1 ! Getting coordinates gridsin(Ngridsin,1) = ix gridsin(Ngridsin,2) = iy within(ix,iy) = .TRUE. CYCLE END IF END DO IF (.NOT.within(ix,iy)) THEN faceA(1,:) = (/ polygonA(NvertexA,1), polygonA(NvertexA,2) /) faceA(2,:) = (/polygonA(1,1), polygonA(1,2) /) CALL intersectfaces(face, faceA, cross, crosspoint) IF (cross /= 0) THEN Ngridsin = Ngridsin + 1 ! Getting coordinates gridsin(Ngridsin,1) = ix gridsin(Ngridsin,2) = iy within(ix,iy) = .TRUE. END IF END IF END IF END IF ! Inverting process, look polygon within grid IF (.NOT.within(ix,iy)) THEN IF (ALLOCATED(gridpoly)) DEALLOCATE(gridpoly) ALLOCATE(gridpoly(Nvertexmax,2)) DO iv = 1, Nvertexmax gridpoly(iv,1) = xBvals(ix,iy,iv) gridpoly(iv,2) = yBvals(ix,iy,iv) END DO DO ivA = 1, NvertexA IF (.NOT.within(ix,iy)) THEN vertex = polygonA(ivA,:) IF (point_inside(vertex, Nvertexmax, gridpoly)) THEN Ngridsin = Ngridsin + 1 ! Getting coordinates gridsin(Ngridsin,1) = ix gridsin(Ngridsin,2) = iy within(ix,iy) = .TRUE. CYCLE END IF END IF END DO END IF END IF END DO END DO IF (ALLOCATED(gridpoly)) DEALLOCATE(gridpoly) END SUBROUTINE grid_within_polygon SUBROUTINE spacepercen_within_reg(NvertexA, polygonA, dx, dy, Nvertexmax, xBvals, yBvals, & Ngridsin, gridsin, strict, percentages) ! Subroutine to compute the percentage of a series of grid cells which are encompassed by a polygon ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x ! and y axes. IMPLICIT NONE INTEGER, INTENT(in) :: NvertexA, dx, dy, Nvertexmax REAL(r_k), DIMENSION(NvertexA,2), INTENT(in) :: polygonA REAL(r_k), DIMENSION(dx,dy,Nvertexmax), INTENT(in) :: xBvals, yBvals INTEGER, INTENT(in) :: Ngridsin INTEGER, DIMENSION(Ngridsin,2), INTENT(in) :: gridsin LOGICAL, INTENT(in) :: strict REAL(r_k), DIMENSION(Ngridsin), INTENT(out) :: percentages ! Local INTEGER :: ig, iv, ix, iy INTEGER :: Nvertex, NvertexAgrid, Ncoin, Nsort CHARACTER(len=20) :: DS REAL(r_k) :: areapoly, areagpoly, totarea, totpercent REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: vertexgrid, icoinpoly, coinpoly, & sortpoly, poly !!!!!!! Variables ! NvertexA: Number of vertices of the polygin to find the grids ! polygonA: ordered vertices of the polygon ! dx, dy: shape of the matrix with the grid points ! xCvals, yCvals: coordinates of the center of the grid cells ! Nvertexmax: Maximum number of vertices of the grid cells ! xBvals, yBvals: coordinates of th vertices of the grid cells (-99999 for no vertex) ! Ngridsin: number of grids with some extension within the polygon ! gridsin: grids within the polygon ! strict: give an error if the area of the polygon is not fully covered ! percentages: percentages of area of each of the grids within the polygon fname = 'spacepercen_within_reg' percentages = zeroRK totpercent = zeroRK totarea = zeroRK areapoly = shoelace_area_polygon(NvertexA, polygonA) DO ig = 1, Ngridsin ix = gridsin(ig,1) iy = gridsin(ig,2) ! Getting grid vertices Nvertex = 0 DO iv=1, Nvertexmax IF (xBvals(ix,iy,iv) /= fillvalI) THEN Nvertex = Nvertex + 1 END IF END DO IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) ALLOCATE(vertexgrid(Nvertex,2)) vertexgrid(:,1) = xBvals(ix,iy,1:Nvertex) vertexgrid(:,2) = yBvals(ix,iy,1:Nvertex) ! Getting common vertices NvertexAgrid = NvertexA*Nvertex*2 IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly) ALLOCATE(icoinpoly(NvertexAgrid,2)) CALL coincident_polygon(NvertexA, Nvertex, NvertexAgrid, polygonA, vertexgrid, Ncoin, icoinpoly) IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly) ALLOCATE(coinpoly(Ncoin,2)) DO iv=1, Ncoin coinpoly(iv,:) = icoinpoly(iv,:) END DO IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly) ALLOCATE(sortpoly(Ncoin,2)) CALL sort_polygon(Ncoin, coinpoly, 1, Nsort, sortpoly) IF (ALLOCATED(poly)) DEALLOCATE(poly) ALLOCATE(poly(Nsort,2)) DO iv=1, Nsort poly(iv,:) = sortpoly(iv,:) END DO areagpoly = shoelace_area_polygon(Nsort, poly) IF (INT(LOG10(EPSILON(totpercent))) < 12) THEN totarea = totarea + ABS(areagpoly) percentages(ig) = ABS(areagpoly / areapoly) ! f2py does not like it! ! ELSE ! totarea = totarea + DABS(areagpoly) ! percentages(ig) = DABS(areagpoly / areapoly) END IF totpercent = totpercent + percentages(ig) END DO IF (INT(LOG10(EPSILON(totpercent))) < 12) THEN IF (strict .AND. ABS(totpercent - oneRK) > epsilonRK) THEN PRINT *, 'totarea:', totarea, ' area polygon:', areapoly PRINT *, 'totpercent:', totpercent, ' oneRK:', oneRK, ' diff:', totpercent - oneRK WRITE(DS,'(F20.8)')ABS(totpercent - oneRK) msg = 'sum of all grid space percentages does not cover (' // TRIM(DS) // ') all polygon' CALL ErrMsg(msg, fname, -1) END IF ELSE ! f2py does not like it! ! IF (strict .AND. ABS(totpercent - oneRK) > epsilonRK) THEN ! PRINT *, 'totarea:', totarea, ' area polygon:', areapoly ! PRINT *, 'totpercent:', totpercent, ' oneRK:', oneRK, ' diff:', totpercent - oneRK ! WRITE(DS,'(F20.16)')ABS(totpercent - oneRK) ! msg = 'sum of all grid space percentages does not cover (' // TRIM(DS) // ') all polygon' ! CALL ErrMsg(msg, fname, -1) ! END IF END IF IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly) IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly) IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly) IF (ALLOCATED(poly)) DEALLOCATE(poly) END SUBROUTINE spacepercen_within_reg SUBROUTINE grid_spacepercen_within_reg(NvertexA, polygonA, dx, dy, Nvertexmax, xBvals, yBvals, & Ngridsin, gridsin, strict, gridspace, percentages) ! Subroutine to compute the percentage of grid space of a series of grid cells which are encompassed ! by a polygon ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x ! and y axes. IMPLICIT NONE INTEGER, INTENT(in) :: NvertexA, dx, dy, Nvertexmax REAL(r_k), DIMENSION(NvertexA,2), INTENT(in) :: polygonA REAL(r_k), DIMENSION(dx,dy,Nvertexmax), INTENT(in) :: xBvals, yBvals INTEGER, INTENT(in) :: Ngridsin INTEGER, DIMENSION(Ngridsin,2), INTENT(in) :: gridsin LOGICAL, INTENT(in) :: strict REAL(r_k), DIMENSION(Ngridsin), INTENT(out) :: gridspace, percentages ! Local INTEGER :: ig, iv, ix, iy INTEGER :: Nvertex, NvertexAgrid, Ncoin, Nsort CHARACTER(len=20) :: DS REAL(r_k) :: areapoly, areagpoly REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: vertexgrid, icoinpoly, coinpoly, & sortpoly, poly !!!!!!! Variables ! NvertexA: Number of vertices of the polygon to find the grids ! polygonA: ordered vertices of the polygon ! dx, dy: shape of the matrix with the grid points ! xCvals, yCvals: coordinates of the center of the grid cells ! Nvertexmax: Maximum number of vertices of the grid cells ! xBvals, yBvals: coordinates of th vertices of the grid cells (-99999 for no vertex) ! Ngridsin: number of grids with some extension within the polygon ! gridsin: grids within the polygon ! strict: give an error if the area of the polygon is not fully covered ! gridspace: area of each of the grids ! percentages: percentages of grid area of each of the grids within the polygon fname = 'grid_spacepercen_within_reg' gridspace = zeroRK percentages = zeroRK DO ig = 1, Ngridsin ix = gridsin(ig,1) iy = gridsin(ig,2) ! Getting grid vertices Nvertex = 0 DO iv=1, Nvertexmax IF (xBvals(ix,iy,iv) /= fillvalI) THEN Nvertex = Nvertex + 1 END IF END DO IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) ALLOCATE(vertexgrid(Nvertex,2)) vertexgrid(:,1) = xBvals(ix,iy,1:Nvertex) vertexgrid(:,2) = yBvals(ix,iy,1:Nvertex) areapoly = shoelace_area_polygon(Nvertex, vertexgrid) ! Getting common vertices NvertexAgrid = NvertexA*Nvertex*2 IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly) ALLOCATE(icoinpoly(NvertexAgrid,2)) CALL coincident_polygon(NvertexA, Nvertex, NvertexAgrid, polygonA, vertexgrid, Ncoin, icoinpoly) IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly) ALLOCATE(coinpoly(Ncoin,2)) DO iv=1, Ncoin coinpoly(iv,:) = icoinpoly(iv,:) END DO IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly) ALLOCATE(sortpoly(Ncoin,2)) CALL sort_polygon(Ncoin, coinpoly, 1, Nsort, sortpoly) IF (ALLOCATED(poly)) DEALLOCATE(poly) ALLOCATE(poly(Nsort,2)) DO iv=1, Nsort poly(iv,:) = sortpoly(iv,:) END DO areagpoly = shoelace_area_polygon(Nsort, poly) gridspace(ig) = ABS(areapoly) percentages(ig) = ABS(areagpoly / areapoly) END DO IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly) IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly) IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly) IF (ALLOCATED(poly)) DEALLOCATE(poly) END SUBROUTINE grid_spacepercen_within_reg SUBROUTINE grid_spacepercen_within_reg_providing_polys(NvertexA, polygonA, dx, dy, Nvertexmax, & xBvals, yBvals, Ngridsin, gridsin, strict, Nmaxver2, Ncoinpoly, ccoinpoly, gridspace, percentages) ! Subroutine to compute the percentage of grid space of a series of grid cells which are encompassed ! by a polygon providing coordinates of the resultant polygons ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x ! and y axes. IMPLICIT NONE INTEGER, INTENT(in) :: NvertexA, dx, dy, Nvertexmax, Nmaxver2 REAL(r_k), DIMENSION(NvertexA,2), INTENT(in) :: polygonA REAL(r_k), DIMENSION(dx,dy,Nvertexmax), INTENT(in) :: xBvals, yBvals INTEGER, INTENT(in) :: Ngridsin INTEGER, DIMENSION(Ngridsin,2), INTENT(in) :: gridsin LOGICAL, INTENT(in) :: strict INTEGER, DIMENSION(Ngridsin), INTENT(out) :: Ncoinpoly REAL(r_k), DIMENSION(Ngridsin,Nmaxver2,2), & INTENT(out) :: ccoinpoly REAL(r_k), DIMENSION(Ngridsin), INTENT(out) :: gridspace, percentages ! Local INTEGER :: ig, iv, ix, iy INTEGER :: Nvertex, NvertexAgrid, Ncoin, Nsort CHARACTER(len=20) :: DS REAL(r_k) :: areapoly, areagpoly REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: vertexgrid, icoinpoly, coinpoly, & sortpoly, poly !!!!!!! Variables ! NvertexA: Number of vertices of the polygon to find the grids ! polygonA: ordered vertices of the polygon ! dx, dy: shape of the matrix with the grid points ! xCvals, yCvals: coordinates of the center of the grid cells ! Nvertexmax: Maximum number of vertices of the grid cells ! xBvals, yBvals: coordinates of th vertices of the grid cells (-99999 for no vertex) ! Ngridsin: number of grids with some extension within the polygon ! gridsin: grids within the polygon ! strict: give an error if the area of the polygon is not fully covered ! Nmaxver2: maximum possible number of vertices of the coincident polygon ! Ncoinpoly: number of vertices of the coincident polygon ! coinpoly: coordinates of the vertices of the coincident polygon ! gridspace: area of each of the grids ! percentages: percentages of grid area of each of the grids within the polygon fname = 'grid_spacepercen_within_reg_providing_polys' gridspace = zeroRK percentages = zeroRK DO ig = 1, Ngridsin ix = gridsin(ig,1) iy = gridsin(ig,2) ! Getting grid vertices Nvertex = 0 DO iv=1, Nvertexmax IF (xBvals(ix,iy,iv) /= fillvalI) THEN Nvertex = Nvertex + 1 END IF END DO IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) ALLOCATE(vertexgrid(Nvertex,2)) vertexgrid(:,1) = xBvals(ix,iy,1:Nvertex) vertexgrid(:,2) = yBvals(ix,iy,1:Nvertex) areapoly = shoelace_area_polygon(Nvertex, vertexgrid) ! Getting common vertices NvertexAgrid = NvertexA*Nvertex*2 IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly) ALLOCATE(icoinpoly(NvertexAgrid,2)) CALL coincident_polygon(NvertexA, Nvertex, NvertexAgrid, polygonA, vertexgrid, Ncoin, icoinpoly) IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly) ALLOCATE(coinpoly(Ncoin,2)) DO iv=1, Ncoin coinpoly(iv,:) = icoinpoly(iv,:) END DO IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly) ALLOCATE(sortpoly(Ncoin,2)) CALL sort_polygon(Ncoin, coinpoly, 1, Nsort, sortpoly) IF (ALLOCATED(poly)) DEALLOCATE(poly) ALLOCATE(poly(Nsort,2)) DO iv=1, Nsort poly(iv,:) = sortpoly(iv,:) END DO areagpoly = shoelace_area_polygon(Nsort, poly) Ncoinpoly(ig)= Nsort ccoinpoly(ig,:,:) = poly(:,:) gridspace(ig) = ABS(areapoly) percentages(ig) = ABS(areagpoly / areapoly) END DO IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) IF (ALLOCATED(icoinpoly)) DEALLOCATE(icoinpoly) IF (ALLOCATED(coinpoly)) DEALLOCATE(coinpoly) IF (ALLOCATED(sortpoly)) DEALLOCATE(sortpoly) IF (ALLOCATED(poly)) DEALLOCATE(poly) END SUBROUTINE grid_spacepercen_within_reg_providing_polys SUBROUTINE spacepercen(xCAvals, yCAvals, xBAvals, yBAvals, xCBvals, yCBvals, xBBvals, yBBvals, & dxA, dyA, NAvertexmax, dxB, dyB, dxyB, NBvertexmax, strict, Ngridsin, gridsin, areas, percentages) ! Subroutine to compute the space-percentages of a series of grid cells (B) into another series of ! grid-cells (A) ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x ! and y axes. IMPLICIT NONE INTEGER, INTENT(in) :: dxA, dyA, NAvertexmax INTEGER, INTENT(in) :: dxB, dyB, NBvertexmax, dxyB REAL(r_k), DIMENSION(dxA,dyA), INTENT(in) :: xCAvals, yCAvals REAL(r_k), DIMENSION(dxB,dyB), INTENT(in) :: xCBvals, yCBvals REAL(r_k), DIMENSION(dxA,dyA,NAvertexmax), INTENT(in):: xBAvals, yBAvals REAL(r_k), DIMENSION(dxB,dyB,NBvertexmax), INTENT(in):: xBBvals, yBBvals LOGICAL, INTENT(in) :: strict INTEGER, DIMENSION(dxA,dyA), INTENT(out) :: Ngridsin INTEGER, DIMENSION(dxA,dyA,dxyB,2), INTENT(out) :: gridsin REAL(r_k), DIMENSION(dxA,dyA), INTENT(out) :: areas REAL(r_k), DIMENSION(dxA,dyA,dxyB), INTENT(out) :: percentages ! Local INTEGER :: iv, ix, iy INTEGER :: Nvertex INTEGER, ALLOCATABLE, DIMENSION(:,:) :: poinsin CHARACTER(len=20) :: IS REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: vertexgrid !!!!!!! Variables ! dxA, dyA: shape of the matrix with the grid points A ! xCAvals, yCAvals: coordinates of the center of the grid cells A ! NAvertexmax: Maximum number of vertices of the grid cells A ! xBAvals, yBAvals: coordinates of th vertices of the grid cells A (-99999 for no vertex) ! dxB, dyB: shape of the matrix with the grid points B ! xCBvals, yCBvals: coordinates of the center of the grid cells B ! NBvertexmax: Maximum number of vertices of the grid cells B ! xBBvals, yBBvals: coordinates of th vertices of the grid cells B (-99999 for no vertex) ! strict: give an error if the area of the polygon is not fully covered ! Ngridsin: number of grids from grid B with some extension within the grid cell A ! gridsin: indices of B grids within the grids of A ! areas: areas of the polygons ! percentages: percentages of area of cells B of each of the grids within the grid cell A fname = 'spacepercen' DO ix = 1, dxA DO iy = 1, dyA ! Getting grid vertices Nvertex = 0 DO iv=1, NAvertexmax IF (xBAvals(ix,iy,iv) /= fillval64) THEN Nvertex = Nvertex + 1 END IF END DO IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) ALLOCATE(vertexgrid(Nvertex,2)) vertexgrid(:,1) = xBAvals(ix,iy,1:Nvertex) vertexgrid(:,2) = yBAvals(ix,iy,1:Nvertex) CALL grid_within_polygon(Nvertex, vertexgrid, dxB, dyB, dxB*dyB, xCBvals, yCBvals, & NBvertexmax, xBBvals, yBBvals, Ngridsin(ix,iy), gridsin(ix,iy,:,:)) IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin) ALLOCATE(poinsin(Ngridsin(ix,iy),2)) DO iv=1, Ngridsin(ix,iy) poinsin(iv,1) = gridsin(ix,iy,iv,1) poinsin(iv,2) = gridsin(ix,iy,iv,2) END DO areas(ix,iy) = shoelace_area_polygon(Nvertex, vertexgrid) CALL spacepercen_within_reg(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, xBBvals, yBBvals, & Ngridsin(ix,iy), poinsin, strict, percentages(ix,iy,:)) END DO END DO IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin) END SUBROUTINE spacepercen SUBROUTINE grid_spacepercen(xCAvals, yCAvals, xBAvals, yBAvals, xCBvals, yCBvals, xBBvals, yBBvals, & dxA, dyA, NAvertexmax, dxB, dyB, dxyB, NBvertexmax, strict, Ngridsin, gridsin, areas2D, areas, & percentages) ! Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another ! series of grid-cells (A) porviding coincident polygons ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x ! and y axes. IMPLICIT NONE INTEGER, INTENT(in) :: dxA, dyA, NAvertexmax INTEGER, INTENT(in) :: dxB, dyB, NBvertexmax, dxyB REAL(r_k), DIMENSION(dxA,dyA), INTENT(in) :: xCAvals, yCAvals REAL(r_k), DIMENSION(dxB,dyB), INTENT(in) :: xCBvals, yCBvals REAL(r_k), DIMENSION(dxA,dyA,NAvertexmax), INTENT(in):: xBAvals, yBAvals REAL(r_k), DIMENSION(dxB,dyB,NBvertexmax), INTENT(in):: xBBvals, yBBvals LOGICAL, INTENT(in) :: strict INTEGER, DIMENSION(dxA,dyA), INTENT(out) :: Ngridsin INTEGER, DIMENSION(dxA,dyA,dxyB,2), INTENT(out) :: gridsin REAL(r_k), DIMENSION(dxB,dyB), INTENT(out) :: areas2D REAL(r_k), DIMENSION(dxA,dyA,dxyB), INTENT(out) :: areas,percentages ! Local INTEGER :: iv, ix, iy INTEGER :: Nvertex, Nptin INTEGER, ALLOCATABLE, DIMENSION(:,:) :: poinsin CHARACTER(len=20) :: IS REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: vertexgrid !!!!!!! Variables ! dxA, dyA: shape of the matrix with the grid points A ! xCAvals, yCAvals: coordinates of the center of the grid cells A ! NAvertexmax: Maximum number of vertices of the grid cells A ! xBAvals, yBAvals: coordinates of th vertices of the grid cells A (-99999 for no vertex) ! dxB, dyB: shape of the matrix with the grid points B ! xCBvals, yCBvals: coordinates of the center of the grid cells B ! NBvertexmax: Maximum number of vertices of the grid cells B ! xBBvals, yBBvals: coordinates of th vertices of the grid cells B (-99999 for no vertex) ! strict: give an error if the area of the polygon is not fully covered ! Ngridsin: number of grids from grid B with some extension within the grid cell A ! gridsin: indices of B grids within the grids of A ! areas2D: areas of the grids as 2D matrix in the original shape ! areas: areas of cells B of each of the grids inside the grid cell A ! percentages: percentages of area of cells B of each of the grids inside the grid cell A fname = 'grid_spacepercen' areas2D = zeroRK areas = zeroRK percentages = zeroRK DO ix = 1, dxA DO iy = 1, dyA ! Getting grid vertices Nvertex = 0 DO iv=1, NAvertexmax IF (xBAvals(ix,iy,iv) /= fillval64) THEN Nvertex = Nvertex + 1 END IF END DO IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) ALLOCATE(vertexgrid(Nvertex,2)) vertexgrid(:,1) = xBAvals(ix,iy,1:Nvertex) vertexgrid(:,2) = yBAvals(ix,iy,1:Nvertex) CALL grid_within_polygon(Nvertex, vertexgrid, dxB, dyB, dxB*dyB, xCBvals, yCBvals, & NBvertexmax, xBBvals, yBBvals, Ngridsin(ix,iy), gridsin(ix,iy,1:dxyB,:)) IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin) ALLOCATE(poinsin(Ngridsin(ix,iy),2)) DO iv=1, Ngridsin(ix,iy) poinsin(iv,1) = gridsin(ix,iy,iv,1) poinsin(iv,2) = gridsin(ix,iy,iv,2) END DO Nptin = Ngridsin(ix,iy) CALL grid_spacepercen_within_reg(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, xBBvals, & yBBvals, Ngridsin(ix,iy), poinsin, strict, areas(ix,iy,1:Nptin), percentages(ix,iy,1:Nptin)) ! Filling areas DO iv = 1, Ngridsin(ix,iy) IF (areas2D(poinsin(iv,1), poinsin(iv,2)) == zeroRK) THEN areas2D(poinsin(iv,1), poinsin(iv,2)) = areas(ix,iy,iv) END IF END DO END DO END DO IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin) END SUBROUTINE grid_spacepercen SUBROUTINE grid_spacepercen_providing_polys(xCAvals, yCAvals, xBAvals, yBAvals, xCBvals, yCBvals, & xBBvals, yBBvals, dxA, dyA, NAvertexmax, dxB, dyB, dxyB, NBvertexmax, strict, Nmaxvercoin, & Nvercoinpolys, vercoinpolys, Ngridsin, gridsin, areas, percentages) ! Subroutine to compute the space-percentages of a series of grid cells (B) which lay inside another ! series of grid-cells (A) providing coincident polygons ! NOTE: Assuming coordinates on the plane with rectilinar, distance preserved and perpendicular x ! and y axes. IMPLICIT NONE INTEGER, INTENT(in) :: dxA, dyA, NAvertexmax INTEGER, INTENT(in) :: dxB, dyB, NBvertexmax, dxyB INTEGER, INTENT(in) :: Nmaxvercoin REAL(r_k), DIMENSION(dxA,dyA), INTENT(in) :: xCAvals, yCAvals REAL(r_k), DIMENSION(dxB,dyB), INTENT(in) :: xCBvals, yCBvals REAL(r_k), DIMENSION(dxA,dyA,NAvertexmax), INTENT(in):: xBAvals, yBAvals REAL(r_k), DIMENSION(dxB,dyB,NBvertexmax), INTENT(in):: xBBvals, yBBvals LOGICAL, INTENT(in) :: strict INTEGER, DIMENSION(dxA,dyA,dxyB,Nmaxvercoin), & INTENT(out) :: Nvercoinpolys REAL(r_k), DIMENSION(dxA,dyA,dxyB,Nmaxvercoin,2), & INTENT(out) :: vercoinpolys INTEGER, DIMENSION(dxA,dyA), INTENT(out) :: Ngridsin INTEGER, DIMENSION(dxA,dyA,dxyB,2), INTENT(out) :: gridsin REAL(r_k), DIMENSION(dxB,dyB), INTENT(out) :: areas REAL(r_k), DIMENSION(dxA,dyA,dxyB), INTENT(out) :: percentages ! Local INTEGER :: iv, ix, iy INTEGER :: Nvertex INTEGER, ALLOCATABLE, DIMENSION(:,:) :: poinsin CHARACTER(len=20) :: IS REAL(r_k), ALLOCATABLE, DIMENSION(:) :: pareas REAL(r_k), ALLOCATABLE, DIMENSION(:,:) :: vertexgrid !!!!!!! Variables ! dxA, dyA: shape of the matrix with the grid points A ! xCAvals, yCAvals: coordinates of the center of the grid cells A ! NAvertexmax: Maximum number of vertices of the grid cells A ! xBAvals, yBAvals: coordinates of th vertices of the grid cells A (-99999 for no vertex) ! dxB, dyB: shape of the matrix with the grid points B ! xCBvals, yCBvals: coordinates of the center of the grid cells B ! NBvertexmax: Maximum number of vertices of the grid cells B ! xBBvals, yBBvals: coordinates of th vertices of the grid cells B (-99999 for no vertex) ! strict: give an error if the area of the polygon is not fully covered ! Nvercoinpolys: number of vertices of the coincident polygon of each grid ! coinpolys: of vertices of the coincident polygon of each grid ! Ngridsin: number of grids from grid B with some extension within the grid cell A ! gridsin: indices of B grids within the grids of A ! areas: areas of the grids ! percentages: percentages of area of cells B of each of the grids inside the grid cell A fname = 'grid_spacepercen_providing_polys' areas = zeroRK DO ix = 1, dxA DO iy = 1, dyA ! Getting grid vertices Nvertex = 0 DO iv=1, NAvertexmax IF (xBAvals(ix,iy,iv) /= fillval64) THEN Nvertex = Nvertex + 1 END IF END DO IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) ALLOCATE(vertexgrid(Nvertex,2)) vertexgrid(:,1) = xBAvals(ix,iy,1:Nvertex) vertexgrid(:,2) = yBAvals(ix,iy,1:Nvertex) CALL grid_within_polygon(Nvertex, vertexgrid, dxB, dyB, dxB*dyB, xCBvals, yCBvals, & NBvertexmax, xBBvals, yBBvals, Ngridsin(ix,iy), gridsin(ix,iy,:,:)) IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin) ALLOCATE(poinsin(Ngridsin(ix,iy),2)) IF (ALLOCATED(pareas)) DEALLOCATE(pareas) ALLOCATE(pareas(Ngridsin(ix,iy))) DO iv=1, Ngridsin(ix,iy) poinsin(iv,1) = gridsin(ix,iy,iv,1) poinsin(iv,2) = gridsin(ix,iy,iv,2) END DO CALL grid_spacepercen_within_reg_providing_polys(Nvertex, vertexgrid, dxB, dyB, NBvertexmax, & xBBvals, yBBvals, Ngridsin(ix,iy), poinsin, strict, Nmaxvercoin, Nvercoinpolys(ix,iy,:,:), & vercoinpolys(ix,iy,:,:,:), pareas, percentages(ix,iy,:)) ! Filling areas DO iv = 1, Ngridsin(ix,iy) IF (areas(poinsin(iv,1), poinsin(iv,2)) == zeroRK) THEN areas(poinsin(iv,1), poinsin(iv,2)) = pareas(iv) END IF END DO END DO END DO IF (ALLOCATED(vertexgrid)) DEALLOCATE(vertexgrid) IF (ALLOCATED(pareas)) DEALLOCATE(pareas) IF (ALLOCATED(poinsin)) DEALLOCATE(poinsin) END SUBROUTINE grid_spacepercen_providing_polys SUBROUTINE unique_matrixRK2D(dx, dy, dxy, matrix2D, Nunique, unique) ! Subroutine to provide the unique values within a 2D RK matrix IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, dxy REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: matrix2D INTEGER, INTENT(out) :: Nunique REAL(r_k), DIMENSION(dxy), INTENT(out) :: unique ! Local INTEGER :: ix, iy, iu, minvalv LOGICAL :: single REAL(r_k), ALLOCATABLE, DIMENSION(:) :: uniques !!!!!!! Variables ! dx, dy: dimensions of the matrix ! dxy: dx*dy, maximum possible amount of different values ! matrix2D: matgrix of values ! Nunique: amount of unique values ! unique: sorted from minimum to maximum vector with the unique values fname = 'unique_matrixRK2D' minvalv = MINVAL(matrix2D) Nunique = 1 unique(1) = minvalv DO ix= 1, dx DO iy= 1, dy single = .TRUE. DO iu = 1, Nunique IF (matrix2D(ix,iy) == unique(iu)) THEN single = .FALSE. EXIT END IF END DO IF (single) THEN Nunique = Nunique + 1 unique(Nunique) = matrix2D(ix,iy) END IF END DO END DO IF (ALLOCATED(uniques)) DEALLOCATE(uniques) ALLOCATE(uniques(Nunique)) uniques(1:Nunique) = unique(1:Nunique) CALL sortR_K(uniques(1:Nunique), Nunique) unique(1:Nunique) = uniques(1:Nunique) END SUBROUTINE unique_matrixRK2D SUBROUTINE spaceweightstats(varin, Ngridsin, gridsin, percentages, stats, varout, dxA, dyA, dxB, & dyB, maxNgridsin, Lstats) ! Subroutine to compute an spatial statistics value from a matrix B into a matrix A using weights IMPLICIT NONE INTEGER, INTENT(in) :: dxA, dyA, dxB, dyB, maxNgridsin, Lstats CHARACTER(len=*), INTENT(in) :: stats INTEGER, DIMENSION(dxA,dyA), INTENT(in) :: Ngridsin INTEGER, DIMENSION(dxA,dyA,maxNgridsin,2), INTENT(in):: gridsin REAL(r_k), DIMENSION(dxB,dyB), INTENT(in) :: varin REAL(r_k), DIMENSION(dxA,dyA,maxNgridsin), INTENT(in):: percentages REAL(r_k), DIMENSION(dxA,dyA,Lstats), INTENT(out) :: varout ! Local INTEGER :: ix, iy, iv, ic, iu, ii, jj INTEGER :: Ncounts CHARACTER(len=3) :: val1S, val2S CHARACTER(len=30) :: val3S REAL(r_k) :: val1, val2 REAL(r_k), DIMENSION(Lstats) :: icounts !!!!!!! Variables ! dxA, dyA: length of dimensions of matrix A ! dxB, dyB: length of dimensions of matrix B ! maxNgridsin: maximum number of grid points from B to be used to compute into a grid of matrix A ! Lstats: length of the dimension of the statistics ! varin: variable from matrix B to be used ! Ngridsin: number of grids from matrix B for each grid of matrix A ! gridsin: coordinates of grids of matrix B for each grid of matrix A ! percentages: weights as percentages of space of grid in matrix A covered by grid of matrix B ! stats: name of the spatial statistics to compute inside each grid of matrix A using values from ! matrix B. Available ones: ! 'min': minimum value ! 'max': maximum value ! 'mean': space weighted mean value ! 'mean2': space weighted quadratic mean value ! 'stddev': space weighted standard deviation value ! 'count': percentage of the space of matrix A covered by each different value of matrix B ! varout: output statistical variable fname = 'spaceweightstats' ! Let's be efficvient? statn: SELECT CASE(TRIM(stats)) CASE('min') varout = fillVal64 DO ix=1, dxA DO iy=1, dyA DO iv=1, Ngridsin(ix,iy) ii = gridsin(ix,iy,iv,1) jj = gridsin(ix,iy,iv,2) IF (varin(ii,jj) < varout(ix,iy,Lstats)) varout(ix,iy,1) = varin(ii,jj) END DO END DO END DO CASE('max') varout = -fillVal64 DO ix=1, dxA DO iy=1, dyA DO iv=1, Ngridsin(ix,iy) ii = gridsin(ix,iy,iv,1) jj = gridsin(ix,iy,iv,2) IF (varin(ii,jj) > varout(ix,iy,Lstats)) varout(ix,iy,1) = varin(ii,jj) END DO END DO END DO CASE('mean') varout = zeroRK DO ix=1, dxA DO iy=1, dyA DO iv=1, Ngridsin(ix,iy) ii = gridsin(ix,iy,iv,1) jj = gridsin(ix,iy,iv,2) varout(ix,iy,1) = varout(ix,iy,1) + varin(ii,jj)*percentages(ix,iy,iv) END DO END DO END DO CASE('mean2') varout = zeroRK DO ix=1, dxA DO iy=1, dyA DO iv=1, Ngridsin(ix,iy) ii = gridsin(ix,iy,iv,1) jj = gridsin(ix,iy,iv,2) varout(ix,iy,1) = varout(ix,iy,1) + percentages(ix,iy,iv)*(varin(ii,jj))**2 END DO varout(ix,iy,1) = varout(ix,iy,1) / Ngridsin(ix,iy) END DO END DO CASE('stddev') varout = zeroRK DO ix=1, dxA DO iy=1, dyA val1 = zeroRK val2 = zeroRK DO iv=1, Ngridsin(ix,iy) ii = gridsin(ix,iy,iv,1) jj = gridsin(ix,iy,iv,2) val1 = val1 + varin(ii,jj)*percentages(ix,iy,iv) val2 = val2 + percentages(ix,iy,iv)*(varin(ii,jj))**2 END DO varout(ix,iy,1) = SQRT(val2 - val1**2) END DO END DO CASE('count') CALL unique_matrixRK2D(dxB, dyB, dxB*dyB, varin, Ncounts, icounts) IF (Lstats /= Ncounts) THEN PRINT *,' ' // TRIM(fname) // 'provided:', Lstats PRINT *,' ' // TRIM(fname) // 'found:', Ncounts, ' :', icounts WRITE(val1S,'(I3)')Lstats WRITE(val2S,'(I3)')Ncounts msg = "for 'count' different amount of passed categories: " // TRIM(val1S) // & ' and found ' // TRIM(val2S) CALL ErrMsg(msg, fname, -1) END IF varout = zeroRK DO ix=1, dxA DO iy=1, dyA DO iv=1, Ngridsin(ix,iy) ii = gridsin(ix,iy,iv,1) jj = gridsin(ix,iy,iv,2) ic = Index1DArrayR_K(icounts, Ncounts, varin(ii,jj)) IF (ic == -1) THEN WRITE(val3S,'(f30.20)')varin(ii,jj) msg = "value '" // val3S // "' for 'count' not found" CALL ErrMSg(msg, fname, -1) ELSE varout(ix,iy,ic) = varout(ix,iy,ic) + percentages(ix,iy,iv) END IF END DO END DO END DO CASE DEFAULT msg = "statisitcs '" // TRIM(stats) // "' not ready !!" // CHAR(44) // " available ones: " // & "'min', 'max', 'mean', 'mean2', 'stddev', 'count'" CALL ErrMsg(msg, fname, -1) END SELECT statn END SUBROUTINE spaceweightstats SUBROUTINE spaceweight_Icount(Ncat, cats, varin, Ngridsin, percentages, varout) ! Subroutine to compute the percentage of a series of categories using weights IMPLICIT NONE INTEGER, INTENT(in) :: Ncat, Ngridsin INTEGER, DIMENSION(Ncat), INTENT(in) :: cats INTEGER, DIMENSION(Ngridsin), INTENT(in) :: varin REAL(r_k), DIMENSION(Ngridsin), INTENT(in) :: percentages REAL(r_k), DIMENSION(Ncat), INTENT(out) :: varout ! Local INTEGER :: iv, ic !!!!!!! Variables ! Ncat: number of categories ! cats: categories to use ! varin: variable with the integer values ! Ngridsin: number of grids ! percentages: weights as percentages of the grids to use ! varout: output statistical variable fname = 'spaceweight_Icount' varout = zeroRK DO iv=1, Ngridsin DO ic=1, Ncat IF (varin(iv) == cats(ic)) THEN varout(ic) = varout(ic) + percentages(iv) END IF END DO END DO END SUBROUTINE spaceweight_Icount SUBROUTINE multi_spaceweightstats_in1DRKno_slc3v3(varin, idv, Ngridsin, gridsin, percentages, & varout, di1, ds1, ds2, ds3, maxNgridsin) ! Subroutine to compute an spatial statistics value from a 1D RK matrix without running one into a ! matrix of 3-variables slices of rank 3 using spatial weights IMPLICIT NONE INTEGER, INTENT(in) :: di1, idv, ds1, ds2, ds3 INTEGER, INTENT(in) :: maxNgridsin INTEGER, DIMENSION(ds1,ds2,ds3), INTENT(in) :: Ngridsin INTEGER, DIMENSION(ds1,ds2,ds3,maxNgridsin,2), & INTENT(in) :: gridsin REAL(r_k), DIMENSION(di1), INTENT(in) :: varin REAL(r_k), INTENT(in), & DIMENSION(ds1,ds2,ds3,maxNgridsin) :: percentages REAL(r_k), DIMENSION(ds1,ds2,ds3,7), INTENT(out) :: varout ! Local INTEGER :: i1, i2, i3, s1, s2, s3, iv INTEGER :: ii3, ss1, ss2, ss3 INTEGER :: Ncounts, Nin INTEGER, DIMENSION(1) :: dmaxvarin CHARACTER(len=3) :: val1S, val2S CHARACTER(len=30) :: val3S REAL(r_k) :: minv, maxv, meanv, mean2v, stdv, medv REAL(r_k), DIMENSION(:), ALLOCATABLE :: pin INTEGER, DIMENSION(:,:), ALLOCATABLE :: gin REAL(r_k), DIMENSION(:), ALLOCATABLE :: svin REAL(r_k), DIMENSION(:), ALLOCATABLE :: vin !!!!!!! Variables ! di1: length of dimensions of the 1D matrix of values ! ds[1-3]: length of dimensions of matrix with the slices ! maxNgridsin: maximum number of grid points from the 3D matrix in any slice ! varin: 1D RK variable to be used ! idv: which dimension of the sliced grids coincide with the dimension of 1D varin ! Ngridsin: number of grids from 3D RK matrix for each slice ! gridsin: coordinates of grids of the 3D RK matrix B to matrix of slices ! percentages: weights as percentages of space of 3D RK matrix for each slice !!!!! ! Available spatial statistics to compute inside each slice using values from 3D RK matrix ! 'min': minimum value ! 'max': maximum value ! 'mean': space weighted mean value ! 'mean2': space weighted quadratic mean value ! 'stddev': space weighted standard deviation value ! 'median': median value ! 'count': percentage of the space of matrix A covered by each different value of matrix B ! varout: output statistical variable fname = 'multi_spaceweightstats_in1DRKno_slc3v3' ss1 = 8 + 1 ss2 = 5 + 1 ss3 = 3 + 1 ii3 = 1 + 1 dmaxvarin = UBOUND(varin) ! Let's be efficient? varout = fillVal64 DO s1 =1, ds1 DO s2 =1, ds2 DO s3 =1, ds3 Nin = Ngridsin(s1,s2,s3) ! Computing along d3 IF (Nin > 1) THEN IF (ALLOCATED(gin)) DEALLOCATE(gin) ALLOCATE(gin(Nin,2)) IF (ALLOCATED(pin)) DEALLOCATE(pin) ALLOCATE(pin(Nin)) IF (ALLOCATED(vin)) DEALLOCATE(vin) ALLOCATE(vin(Nin)) IF (ALLOCATED(svin)) DEALLOCATE(svin) ALLOCATE(svin(Nin)) gin = gridsin(s1,s2,s3,1:Nin,:) pin = percentages(s1,s2,s3,1:Nin) ! Getting the values DO iv=1, Nin i1 = gin(iv,idv) vin(iv) = varin(i1) END DO minv = fillVal64 maxv = -fillVal64 meanv = zeroRK mean2v = zeroRK stdv = zeroRK minv = MINVAL(vin) maxv = MAXVAL(vin) meanv = SUM(vin*pin) mean2v = SUM(vin**2*pin) DO iv=1,Nin stdv = stdv + ( (meanv - vin(iv))*pin(iv) )**2 END DO stdv = SQRT(stdv) svin = vin(:) CALL SortR_K(svin, Nin) medv = svin(INT(Nin/2)) varout(s1,s2,s3,1) = minv varout(s1,s2,s3,2) = maxv varout(s1,s2,s3,3) = meanv varout(s1,s2,s3,4) = mean2v varout(s1,s2,s3,5) = stdv varout(s1,s2,s3,6) = medv varout(s1,s2,s3,7) = Nin*1. ELSE IF (Nin == 1) THEN i1 = gridsin(s1,s2,s3,1,idv) IF (i1 > 0 .AND. i1 <= dmaxvarin(1)) THEN varout(s1,s2,s3,1) = varin(i1) varout(s1,s2,s3,2) = varin(i1) varout(s1,s2,s3,3) = varin(i1) varout(s1,s2,s3,4) = varin(i1)*varin(i1) varout(s1,s2,s3,5) = zeroRK varout(s1,s2,s3,6) = varin(i1) varout(s1,s2,s3,7) = Nin*1. ELSE varout(s1,s2,s3,1) = fillval64 varout(s1,s2,s3,2) = fillval64 varout(s1,s2,s3,3) = fillval64 varout(s1,s2,s3,4) = fillval64 varout(s1,s2,s3,5) = fillval64 varout(s1,s2,s3,6) = fillval64 varout(s1,s2,s3,7) = zeroRK END IF ELSE varout(s1,s2,s3,1) = fillval64 varout(s1,s2,s3,2) = fillval64 varout(s1,s2,s3,3) = fillval64 varout(s1,s2,s3,4) = fillval64 varout(s1,s2,s3,5) = fillval64 varout(s1,s2,s3,6) = fillval64 varout(s1,s2,s3,7) = zeroRK END IF END DO END DO END DO IF (ALLOCATED(gin)) DEALLOCATE(gin) IF (ALLOCATED(pin)) DEALLOCATE(pin) IF (ALLOCATED(vin)) DEALLOCATE(vin) IF (ALLOCATED(svin)) DEALLOCATE(svin) RETURN END SUBROUTINE multi_spaceweightstats_in1DRKno_slc3v3 SUBROUTINE multi_spaceweightstats_in2DRKno_slc3v3(varin, Ngridsin, gridsin, percentages, varout, & di1, di2, ds1, ds2, ds3, maxNgridsin) ! Subroutine to compute an spatial statistics value from a 2D RK matrix without running one into a ! matrix of 3-variables slices of rank 3 using spatial weights IMPLICIT NONE INTEGER, INTENT(in) :: di1, di2, ds1, ds2, ds3 INTEGER, INTENT(in) :: maxNgridsin INTEGER, DIMENSION(ds1,ds2,ds3), INTENT(in) :: Ngridsin INTEGER, DIMENSION(ds1,ds2,ds3,maxNgridsin,2), & INTENT(in) :: gridsin REAL(r_k), DIMENSION(di1,di2), INTENT(in) :: varin REAL(r_k), INTENT(in), & DIMENSION(ds1,ds2,ds3,maxNgridsin) :: percentages REAL(r_k), DIMENSION(ds1,ds2,ds3,7), INTENT(out) :: varout ! Local INTEGER :: i1, i2, i3, s1, s2, s3, iv INTEGER :: ii3, ss1, ss2, ss3 INTEGER :: Ncounts, Nin CHARACTER(len=3) :: val1S, val2S CHARACTER(len=30) :: val3S REAL(r_k) :: minv, maxv, meanv, mean2v, stdv, medv REAL(r_k), DIMENSION(:), ALLOCATABLE :: pin INTEGER, DIMENSION(:,:), ALLOCATABLE :: gin REAL(r_k), DIMENSION(:), ALLOCATABLE :: svin REAL(r_k), DIMENSION(:), ALLOCATABLE :: vin !!!!!!! Variables ! di1, di2: length of dimensions of the 2D matrix of values ! ds[1-3]: length of dimensions of matrix with the slices ! maxNgridsin: maximum number of grid points from the 3D matrix in any slice ! varin: 2D RK variable to be used ! Ngridsin: number of grids from 3D RK matrix for each slice ! gridsin: coordinates of grids of the 3D RK matrix B to matrix of slices ! percentages: weights as percentages of space of 3D RK matrix for each slice !!!!! ! Available spatial statistics to compute inside each slice using values from 3D RK matrix ! 'min': minimum value ! 'max': maximum value ! 'mean': space weighted mean value ! 'mean2': space weighted quadratic mean value ! 'stddev': space weighted standard deviation value ! 'median': median value ! 'count': percentage of the space of matrix A covered by each different value of matrix B ! varout: output statistical variable fname = 'multi_spaceweightstats_in2DRKno_slc3v3' ss1 = 8 + 1 ss2 = 5 + 1 ss3 = 3 + 1 ii3 = 1 + 1 ! Let's be efficient? varout = fillVal64 DO s1 =1, ds1 DO s2 =1, ds2 DO s3 =1, ds3 Nin = Ngridsin(s1,s2,s3) ! Computing along d3 IF (Nin > 1) THEN IF (ALLOCATED(gin)) DEALLOCATE(gin) ALLOCATE(gin(Nin,2)) IF (ALLOCATED(pin)) DEALLOCATE(pin) ALLOCATE(pin(Nin)) IF (ALLOCATED(vin)) DEALLOCATE(vin) ALLOCATE(vin(Nin)) IF (ALLOCATED(svin)) DEALLOCATE(svin) ALLOCATE(svin(Nin)) gin = gridsin(s1,s2,s3,1:Nin,:) pin = percentages(s1,s2,s3,1:Nin) ! Getting the values DO iv=1, Nin i1 = gin(iv,1) i2 = gin(iv,2) vin(iv) = varin(i1,i2) END DO minv = fillVal64 maxv = -fillVal64 meanv = zeroRK mean2v = zeroRK stdv = zeroRK minv = MINVAL(vin) maxv = MAXVAL(vin) meanv = SUM(vin*pin) mean2v = SUM(vin**2*pin) DO iv=1,Nin stdv = stdv + ( (meanv - vin(iv))*pin(iv) )**2 END DO stdv = SQRT(stdv) svin = vin(:) CALL SortR_K(svin, Nin) medv = svin(INT(Nin/2)) varout(s1,s2,s3,1) = minv varout(s1,s2,s3,2) = maxv varout(s1,s2,s3,3) = meanv varout(s1,s2,s3,4) = mean2v varout(s1,s2,s3,5) = stdv varout(s1,s2,s3,6) = medv varout(s1,s2,s3,7) = Nin*1. ELSE IF (Nin == 1) THEN i1 = gridsin(s1,s2,s3,1,1) i2 = gridsin(s1,s2,s3,1,2) varout(s1,s2,s3,1) = varin(i1,i2) varout(s1,s2,s3,2) = varin(i1,i2) varout(s1,s2,s3,3) = varin(i1,i2) varout(s1,s2,s3,4) = varin(i1,i2)*varin(i1,i2) varout(s1,s2,s3,5) = zeroRK varout(s1,s2,s3,6) = varin(i1,i2) varout(s1,s2,s3,7) = Nin*1. ELSE varout(s1,s2,s3,1) = fillval64 varout(s1,s2,s3,2) = fillval64 varout(s1,s2,s3,3) = fillval64 varout(s1,s2,s3,4) = fillval64 varout(s1,s2,s3,5) = fillval64 varout(s1,s2,s3,6) = fillval64 varout(s1,s2,s3,7) = zeroRK END IF END DO END DO END DO IF (ALLOCATED(gin)) DEALLOCATE(gin) IF (ALLOCATED(pin)) DEALLOCATE(pin) IF (ALLOCATED(vin)) DEALLOCATE(vin) IF (ALLOCATED(svin)) DEALLOCATE(svin) RETURN END SUBROUTINE multi_spaceweightstats_in2DRKno_slc3v3 SUBROUTINE multi_spaceweightstats_in2DRKno_slc2v2(varin, Ngridsin, gridsin, percentages, varout, & di1, di2, ds1, ds2, maxNgridsin) ! Subroutine to compute an spatial statistics value from a 2D RK matrix without running one into a ! matrix of 2-variables slices of rank 2 using spatial weights IMPLICIT NONE INTEGER, INTENT(in) :: di1, di2, ds1, ds2 INTEGER, INTENT(in) :: maxNgridsin INTEGER, DIMENSION(ds1,ds2), INTENT(in) :: Ngridsin INTEGER, DIMENSION(ds1,ds2,maxNgridsin,2), & INTENT(in) :: gridsin REAL(r_k), DIMENSION(di1,di2), INTENT(in) :: varin REAL(r_k), INTENT(in), & DIMENSION(ds1,ds2,maxNgridsin) :: percentages REAL(r_k), DIMENSION(ds1,ds2,7), INTENT(out) :: varout ! Local INTEGER :: i1, i2, i3, s1, s2, iv INTEGER :: ii3, ss1, ss2 INTEGER :: Ncounts, Nin CHARACTER(len=3) :: val1S, val2S CHARACTER(len=30) :: val3S REAL(r_k) :: minv, maxv, meanv, mean2v, stdv, medv REAL(r_k), DIMENSION(:), ALLOCATABLE :: pin INTEGER, DIMENSION(:,:), ALLOCATABLE :: gin REAL(r_k), DIMENSION(:), ALLOCATABLE :: svin REAL(r_k), DIMENSION(:), ALLOCATABLE :: vin !!!!!!! Variables ! di1, di2: length of dimensions of the 2D matrix of values ! ds[1-2]: length of dimensions of matrix with the slices ! maxNgridsin: maximum number of grid points from the 2D matrix in any slice ! varin: 2D RK variable to be used ! Ngridsin: number of grids from 2D RK matrix for each slice ! gridsin: coordinates of grids of the 2D RK matrix B to matrix of slices ! percentages: weights as percentages of space of 2D RK matrix for each slice !!!!! ! Available spatial statistics to compute inside each slice using values from 2D RK matrix ! 'min': minimum value ! 'max': maximum value ! 'mean': space weighted mean value ! 'mean2': space weighted quadratic mean value ! 'stddev': space weighted standard deviation value ! 'median': median value ! 'count': percentage of the space of matrix A covered by each different value of matrix B ! varout: output statistical variable fname = 'multi_spaceweightstats_in2DRKno_slc2v2' ss1 = 8 + 1 ss2 = 5 + 1 ii3 = 1 + 1 ! Let's be efficient? varout = fillVal64 DO s1 =1, ds1 DO s2 =1, ds2 Nin = Ngridsin(s1,s2) ! Computing along d2 IF (Nin > 1) THEN IF (ALLOCATED(gin)) DEALLOCATE(gin) ALLOCATE(gin(Nin,2)) IF (ALLOCATED(pin)) DEALLOCATE(pin) ALLOCATE(pin(Nin)) IF (ALLOCATED(vin)) DEALLOCATE(vin) ALLOCATE(vin(Nin)) IF (ALLOCATED(svin)) DEALLOCATE(svin) ALLOCATE(svin(Nin)) gin = gridsin(s1,s2,1:Nin,:) pin = percentages(s1,s2,1:Nin) ! Getting the values DO iv=1, Nin i1 = gin(iv,1) i2 = gin(iv,2) vin(iv) = varin(i1,i2) END DO minv = fillVal64 maxv = -fillVal64 meanv = zeroRK mean2v = zeroRK stdv = zeroRK minv = MINVAL(vin) maxv = MAXVAL(vin) meanv = SUM(vin*pin) mean2v = SUM(vin**2*pin) DO iv=1,Nin stdv = stdv + ( (meanv - vin(iv))*pin(iv) )**2 END DO stdv = SQRT(stdv) svin = vin(:) CALL SortR_K(svin, Nin) medv = svin(INT(Nin/2)) varout(s1,s2,1) = minv varout(s1,s2,2) = maxv varout(s1,s2,3) = meanv varout(s1,s2,4) = mean2v varout(s1,s2,5) = stdv varout(s1,s2,6) = medv varout(s1,s2,7) = Nin*1. ELSE IF (Nin == 1) THEN i1 = gridsin(s1,s2,1,1) i2 = gridsin(s1,s2,1,2) varout(s1,s2,1) = varin(i1,i2) varout(s1,s2,2) = varin(i1,i2) varout(s1,s2,3) = varin(i1,i2) varout(s1,s2,4) = varin(i1,i2)*varin(i1,i2) varout(s1,s2,5) = zeroRK varout(s1,s2,6) = varin(i1,i2) varout(s1,s2,7) = Nin*1. ELSE varout(s1,s2,1) = fillval64 varout(s1,s2,2) = fillval64 varout(s1,s2,3) = fillval64 varout(s1,s2,4) = fillval64 varout(s1,s2,5) = fillval64 varout(s1,s2,6) = fillval64 varout(s1,s2,7) = zeroRK END IF END DO END DO IF (ALLOCATED(gin)) DEALLOCATE(gin) IF (ALLOCATED(pin)) DEALLOCATE(pin) IF (ALLOCATED(vin)) DEALLOCATE(vin) IF (ALLOCATED(svin)) DEALLOCATE(svin) RETURN END SUBROUTINE multi_spaceweightstats_in2DRKno_slc2v2 SUBROUTINE multi_spaceweightstats_in3DRK3_slc2v2(varin, Ngridsin, gridsin, percentages, varout, & di1, di2, di3, ds1, ds2, maxNgridsin) ! Subroutine to compute an spatial statistics value from a 3D RK matrix using 3rd dimension as ! running one into a matrix of 2-variables slices of rank 2 using spatial weights IMPLICIT NONE INTEGER, INTENT(in) :: di1, di2, di3, ds1, ds2 INTEGER, INTENT(in) :: maxNgridsin INTEGER, DIMENSION(ds1,ds2), INTENT(in) :: Ngridsin INTEGER, DIMENSION(ds1,ds2,maxNgridsin,2), & INTENT(in) :: gridsin REAL(r_k), DIMENSION(di1,di2,di3), INTENT(in) :: varin REAL(r_k), INTENT(in), & DIMENSION(ds1,ds2,maxNgridsin) :: percentages REAL(r_k), DIMENSION(ds1,ds2,di3,7), INTENT(out) :: varout ! Local INTEGER :: i1, i2, i3, s1, s2, iv INTEGER :: ii3, ss1, ss2 INTEGER :: Ncounts, Nin CHARACTER(len=3) :: val1S, val2S CHARACTER(len=30) :: val3S REAL(r_k) :: minv, maxv, meanv, mean2v, stdv, medv REAL(r_k), DIMENSION(:), ALLOCATABLE :: pin INTEGER, DIMENSION(:,:), ALLOCATABLE :: gin REAL(r_k), DIMENSION(:), ALLOCATABLE :: svin REAL(r_k), DIMENSION(:,:), ALLOCATABLE :: vin !!!!!!! Variables ! di1, di2, di3: length of dimensions of the 3D matrix of values ! ds[1-2]: length of dimensions of matrix with the slices ! maxNgridsin: maximum number of grid points from the 3D matrix in any slice ! varin: 3D RK variable to be used ! Ngridsin: number of grids from 3D RK matrix for each slice ! gridsin: coordinates of grids of the 3D RK matrix B to matrix of slices ! percentages: weights as percentages of space of 3D RK matrix for each slice !!!!! ! Available spatial statistics to compute inside each slice using values from 3D RK matrix ! 'min': minimum value ! 'max': maximum value ! 'mean': space weighted mean value ! 'mean2': space weighted quadratic mean value ! 'stddev': space weighted standard deviation value ! 'median': median value ! 'count': percentage of the space of matrix A covered by each different value of matrix B ! varout: output statistical variable fname = 'multi_spaceweightstats_in3DRK3_slc2v2' ss1 = 8 + 1 ss2 = 5 + 1 ii3 = 1 + 1 ! Let's be efficient? varout = fillVal64 DO s1 =1, ds1 DO s2 =1, ds2 Nin = Ngridsin(s1,s2) ! Computing along d3 IF (Nin > 1) THEN IF (ALLOCATED(gin)) DEALLOCATE(gin) ALLOCATE(gin(Nin,2)) IF (ALLOCATED(pin)) DEALLOCATE(pin) ALLOCATE(pin(Nin)) IF (ALLOCATED(vin)) DEALLOCATE(vin) ALLOCATE(vin(Nin,di3)) IF (ALLOCATED(svin)) DEALLOCATE(svin) ALLOCATE(svin(Nin)) gin = gridsin(s1,s2,1:Nin,:) pin = percentages(s1,s2,1:Nin) ! Getting the values DO iv=1, Nin i1 = gin(iv,1) i2 = gin(iv,2) vin(iv,:) = varin(i1,i2,:) END DO DO i3=1, di3 minv = fillVal64 maxv = -fillVal64 meanv = zeroRK mean2v = zeroRK stdv = zeroRK minv = MINVAL(vin(:,i3)) maxv = MAXVAL(vin(:,i3)) meanv = SUM(vin(:,i3)*pin) mean2v = SUM(vin(:,i3)**2*pin) DO iv=1,Nin stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2 END DO stdv = SQRT(stdv) svin = vin(:,i3) CALL SortR_K(svin, Nin) medv = svin(INT(Nin/2)) varout(s1,s2,i3,1) = minv varout(s1,s2,i3,2) = maxv varout(s1,s2,i3,3) = meanv varout(s1,s2,i3,4) = mean2v varout(s1,s2,i3,5) = stdv varout(s1,s2,i3,6) = medv varout(s1,s2,i3,7) = Nin*1. END DO ELSE IF (Nin == 1) THEN i1 = gridsin(s1,s2,1,1) i2 = gridsin(s1,s2,1,2) varout(s1,s2,:,1) = varin(i1,i2,:) varout(s1,s2,:,2) = varin(i1,i2,:) varout(s1,s2,:,3) = varin(i1,i2,:) varout(s1,s2,:,4) = varin(i1,i2,:)*varin(i1,i2,:) varout(s1,s2,:,5) = zeroRK varout(s1,s2,:,6) = varin(i1,i2,:) varout(s1,s2,:,7) = Nin*1. ELSE varout(s1,s2,:,1) = fillval64 varout(s1,s2,:,2) = fillval64 varout(s1,s2,:,3) = fillval64 varout(s1,s2,:,4) = fillval64 varout(s1,s2,:,5) = fillval64 varout(s1,s2,:,6) = fillval64 varout(s1,s2,:,7) = zeroRK END IF END DO END DO IF (ALLOCATED(gin)) DEALLOCATE(gin) IF (ALLOCATED(pin)) DEALLOCATE(pin) IF (ALLOCATED(vin)) DEALLOCATE(vin) IF (ALLOCATED(svin)) DEALLOCATE(svin) RETURN END SUBROUTINE multi_spaceweightstats_in3DRK3_slc2v2 SUBROUTINE multi_spaceweightstats_in3DRK3_slc3v3(varin, Ngridsin, gridsin, percentages, varout, & di1, di2, di3, ds1, ds2, ds3, maxNgridsin) ! Subroutine to compute an spatial statistics value from a 3D RK matrix using 3rd dimension as ! running one into a matrix of 3-variables slices of rank 3 using spatial weights IMPLICIT NONE INTEGER, INTENT(in) :: di1, di2, di3, ds1, ds2, ds3 INTEGER, INTENT(in) :: maxNgridsin INTEGER, DIMENSION(ds1,ds2,ds3), INTENT(in) :: Ngridsin INTEGER, DIMENSION(ds1,ds2,ds3,maxNgridsin,2), & INTENT(in) :: gridsin REAL(r_k), DIMENSION(di1,di2,di3), INTENT(in) :: varin REAL(r_k), INTENT(in), & DIMENSION(ds1,ds2,ds3,maxNgridsin) :: percentages REAL(r_k), DIMENSION(ds1,ds2,ds3,di3,7), INTENT(out) :: varout ! Local INTEGER :: i1, i2, i3, s1, s2, s3, iv INTEGER :: ii3, ss1, ss2, ss3 INTEGER :: Ncounts, Nin CHARACTER(len=3) :: val1S, val2S CHARACTER(len=30) :: val3S REAL(r_k) :: minv, maxv, meanv, mean2v, stdv, medv REAL(r_k), DIMENSION(:), ALLOCATABLE :: pin INTEGER, DIMENSION(:,:), ALLOCATABLE :: gin REAL(r_k), DIMENSION(:), ALLOCATABLE :: svin REAL(r_k), DIMENSION(:,:), ALLOCATABLE :: vin !!!!!!! Variables ! di1, di2, di3: length of dimensions of the 3D matrix of values ! ds[1-3]: length of dimensions of matrix with the slices ! maxNgridsin: maximum number of grid points from the 3D matrix in any slice ! varin: 3D RK variable to be used ! Ngridsin: number of grids from 3D RK matrix for each slice ! gridsin: coordinates of grids of the 3D RK matrix B to matrix of slices ! percentages: weights as percentages of space of 3D RK matrix for each slice !!!!! ! Available spatial statistics to compute inside each slice using values from 3D RK matrix ! 'min': minimum value ! 'max': maximum value ! 'mean': space weighted mean value ! 'mean2': space weighted quadratic mean value ! 'stddev': space weighted standard deviation value ! 'median': median value ! 'count': percentage of the space of matrix A covered by each different value of matrix B ! varout: output statistical variable fname = 'multi_spaceweightstats_in3DRK3_slc3v3' ss1 = 8 + 1 ss2 = 5 + 1 ss3 = 3 + 1 ii3 = 1 + 1 ! Let's be efficient? varout = fillVal64 DO s1 =1, ds1 DO s2 =1, ds2 DO s3 =1, ds3 Nin = Ngridsin(s1,s2,s3) ! Computing along d3 IF (Nin > 1) THEN IF (ALLOCATED(gin)) DEALLOCATE(gin) ALLOCATE(gin(Nin,2)) IF (ALLOCATED(pin)) DEALLOCATE(pin) ALLOCATE(pin(Nin)) IF (ALLOCATED(vin)) DEALLOCATE(vin) ALLOCATE(vin(Nin,di3)) IF (ALLOCATED(svin)) DEALLOCATE(svin) ALLOCATE(svin(Nin)) gin = gridsin(s1,s2,s3,1:Nin,:) pin = percentages(s1,s2,s3,1:Nin) ! Getting the values DO iv=1, Nin i1 = gin(iv,1) i2 = gin(iv,2) vin(iv,:) = varin(i1,i2,:) END DO DO i3=1, di3 minv = fillVal64 maxv = -fillVal64 meanv = zeroRK mean2v = zeroRK stdv = zeroRK minv = MINVAL(vin(:,i3)) maxv = MAXVAL(vin(:,i3)) meanv = SUM(vin(:,i3)*pin) mean2v = SUM(vin(:,i3)**2*pin) DO iv=1,Nin stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2 END DO stdv = SQRT(stdv) svin = vin(:,i3) CALL SortR_K(svin, Nin) medv = svin(INT(Nin/2)) varout(s1,s2,s3,i3,1) = minv varout(s1,s2,s3,i3,2) = maxv varout(s1,s2,s3,i3,3) = meanv varout(s1,s2,s3,i3,4) = mean2v varout(s1,s2,s3,i3,5) = stdv varout(s1,s2,s3,i3,6) = medv varout(s1,s2,s3,i3,7) = Nin*1. END DO ELSE IF (Nin == 1) THEN i1 = gridsin(s1,s2,s3,1,1) i2 = gridsin(s1,s2,s3,1,2) varout(s1,s2,s3,:,1) = varin(i1,i2,:) varout(s1,s2,s3,:,2) = varin(i1,i2,:) varout(s1,s2,s3,:,3) = varin(i1,i2,:) varout(s1,s2,s3,:,4) = varin(i1,i2,:)*varin(i1,i2,:) varout(s1,s2,s3,:,5) = zeroRK varout(s1,s2,s3,:,6) = varin(i1,i2,:) varout(s1,s2,s3,:,7) = Nin*1. ELSE varout(s1,s2,s3,:,1) = fillval64 varout(s1,s2,s3,:,2) = fillval64 varout(s1,s2,s3,:,3) = fillval64 varout(s1,s2,s3,:,4) = fillval64 varout(s1,s2,s3,:,5) = fillval64 varout(s1,s2,s3,:,6) = fillval64 varout(s1,s2,s3,:,7) = zeroRK END IF END DO END DO END DO IF (ALLOCATED(gin)) DEALLOCATE(gin) IF (ALLOCATED(pin)) DEALLOCATE(pin) IF (ALLOCATED(vin)) DEALLOCATE(vin) IF (ALLOCATED(svin)) DEALLOCATE(svin) RETURN END SUBROUTINE multi_spaceweightstats_in3DRK3_slc3v3 SUBROUTINE multi_spaceweightstats_in4DRK3_4_slc3v3(varin, Ngridsin, gridsin, percentages, varout, & di1, di2, di3, di4, ds1, ds2, ds3, maxNgridsin) ! Subroutine to compute an spatial statistics value from a 4D RK matrix using 3rd and 4th dimensions ! as running ones into a matrix of 3-variables slices of rank 3 using spatial weights IMPLICIT NONE INTEGER, INTENT(in) :: di1, di2, di3, di4, ds1, ds2, ds3 INTEGER, INTENT(in) :: maxNgridsin INTEGER, DIMENSION(ds1,ds2,ds3), INTENT(in) :: Ngridsin INTEGER, DIMENSION(ds1,ds2,ds3,maxNgridsin,2), & INTENT(in) :: gridsin REAL(r_k), DIMENSION(di1,di2,di3,di4), INTENT(in) :: varin REAL(r_k), INTENT(in), & DIMENSION(ds1,ds2,ds3,maxNgridsin) :: percentages REAL(r_k), DIMENSION(ds1,ds2,ds3,di3,di4,7), & INTENT(out) :: varout ! Local INTEGER :: i1, i2, i3, i4, s1, s2, s3, iv INTEGER :: ii3, ss1, ss2, ss3 INTEGER :: Ncounts, Nin CHARACTER(len=3) :: val1S, val2S CHARACTER(len=30) :: val3S REAL(r_k) :: minv, maxv, meanv, mean2v, stdv, medv REAL(r_k), DIMENSION(:), ALLOCATABLE :: pin INTEGER, DIMENSION(:,:), ALLOCATABLE :: gin REAL(r_k), DIMENSION(:), ALLOCATABLE :: svin REAL(r_k), DIMENSION(:,:,:), ALLOCATABLE :: vin !!!!!!! Variables ! di1, di2, di3, di4: length of dimensions of the 4D matrix of values ! ds[1-3]: length of dimensions of matrix with the slices ! maxNgridsin: maximum number of grid points from the 3D matrix in any slice ! varin: 3D RK variable to be used ! Ngridsin: number of grids from 3D RK matrix for each slice ! gridsin: coordinates of grids of the 3D RK matrix B to matrix of slices ! percentages: weights as percentages of space of 3D RK matrix for each slice !!!!! ! Available spatial statistics to compute inside each slice using values from 3D RK matrix ! 'min': minimum value ! 'max': maximum value ! 'mean': space weighted mean value ! 'mean2': space weighted quadratic mean value ! 'stddev': space weighted standard deviation value ! 'median': median value ! 'count': percentage of the space of matrix A covered by each different value of matrix B ! varout: output statistical variable fname = 'multi_spaceweightstats_in4DRK3_4_slc3v3' ! Let's be efficient? varout = fillVal64 DO s1 =1, ds1 DO s2 =1, ds2 DO s3 =1, ds3 Nin = Ngridsin(s1,s2,s3) ! Computing along d3 and d4 IF (Nin > 1) THEN IF (ALLOCATED(gin)) DEALLOCATE(gin) ALLOCATE(gin(Nin,2)) IF (ALLOCATED(pin)) DEALLOCATE(pin) ALLOCATE(pin(Nin)) IF (ALLOCATED(vin)) DEALLOCATE(vin) ALLOCATE(vin(Nin,di3,di4)) IF (ALLOCATED(svin)) DEALLOCATE(svin) ALLOCATE(svin(Nin)) gin = gridsin(s1,s2,s3,1:Nin,:) pin = percentages(s1,s2,s3,1:Nin) ! Getting the values DO iv=1, Nin i1 = gin(iv,1) i2 = gin(iv,2) vin(iv,:,:) = varin(i1,i2,:,:) END DO DO i3=1, di3 DO i4=1, di4 minv = fillVal64 maxv = -fillVal64 meanv = zeroRK mean2v = zeroRK stdv = zeroRK minv = MINVAL(vin(:,i3,i4)) maxv = MAXVAL(vin(:,i3,i4)) meanv = SUM(vin(:,i3,i4)*pin) mean2v = SUM(vin(:,i3,i4)**2*pin) DO iv=1,Nin stdv = stdv + ( (meanv - vin(iv,i3,i4))*pin(iv) )**2 END DO stdv = SQRT(stdv) svin = vin(:,i3,i4) CALL SortR_K(svin, Nin) medv = svin(INT(Nin/2)) varout(s1,s2,s3,i3,i4,1) = minv varout(s1,s2,s3,i3,i4,2) = maxv varout(s1,s2,s3,i3,i4,3) = meanv varout(s1,s2,s3,i3,i4,4) = mean2v varout(s1,s2,s3,i3,i4,5) = stdv varout(s1,s2,s3,i3,i4,6) = medv varout(s1,s2,s3,i3,i4,7) = Nin*1. END DO END DO ELSE IF (Nin == 1) THEN i1 = gridsin(s1,s2,s3,1,1) i2 = gridsin(s1,s2,s3,1,2) varout(s1,s2,s3,:,:,1) = varin(i1,i2,:,:) varout(s1,s2,s3,:,:,2) = varin(i1,i2,:,:) varout(s1,s2,s3,:,:,3) = varin(i1,i2,:,:) varout(s1,s2,s3,:,:,4) = varin(i1,i2,:,:)*varin(i1,i2,:,:) varout(s1,s2,s3,:,:,5) = zeroRK varout(s1,s2,s3,:,:,6) = varin(i1,i2,:,:) varout(s1,s2,s3,:,:,7) = Nin*1. ELSE varout(s1,s2,s3,:,:,1) = fillval64 varout(s1,s2,s3,:,:,2) = fillval64 varout(s1,s2,s3,:,:,3) = fillval64 varout(s1,s2,s3,:,:,4) = fillval64 varout(s1,s2,s3,:,:,5) = fillval64 varout(s1,s2,s3,:,:,6) = fillval64 varout(s1,s2,s3,:,:,7) = zeroRK END IF END DO END DO END DO IF (ALLOCATED(gin)) DEALLOCATE(gin) IF (ALLOCATED(pin)) DEALLOCATE(pin) IF (ALLOCATED(vin)) DEALLOCATE(vin) IF (ALLOCATED(svin)) DEALLOCATE(svin) RETURN END SUBROUTINE multi_spaceweightstats_in4DRK3_4_slc3v3 SUBROUTINE multi_spaceweightstats_in3DRK3_slc3v4(varin, Ngridsin, gridsin, percentages, varout, & di1, di2, di3, ds1, ds2, ds3, ds4, maxNgridsin) ! Subroutine to compute an spatial statistics value from a 3D RK matrix using 3rd dimension as ! running one into a matrix of 3-variables slices of rank 4 using spatial weights IMPLICIT NONE INTEGER, INTENT(in) :: di1, di2, di3, ds1, ds2, ds3, ds4 INTEGER, INTENT(in) :: maxNgridsin INTEGER, DIMENSION(ds1,ds2,ds3,ds4), INTENT(in) :: Ngridsin INTEGER, INTENT(in), & DIMENSION(ds1,ds2,ds3,ds4,maxNgridsin,2) :: gridsin REAL(r_k), DIMENSION(di1,di2,di3), INTENT(in) :: varin REAL(r_k), INTENT(in), & DIMENSION(ds1,ds2,ds3,ds4,maxNgridsin) :: percentages REAL(r_k), DIMENSION(ds1,ds2,ds3,ds4,di3,7), & INTENT(out) :: varout ! Local INTEGER :: i1, i2, i3, s1, s2, s3, s4, iv INTEGER :: Ncounts, Nin CHARACTER(len=3) :: val1S, val2S CHARACTER(len=30) :: val3S REAL(r_k) :: minv, maxv, meanv, mean2v, stdv, medv REAL(r_k), DIMENSION(:), ALLOCATABLE :: pin INTEGER, DIMENSION(:,:), ALLOCATABLE :: gin REAL(r_k), DIMENSION(:), ALLOCATABLE :: svin REAL(r_k), DIMENSION(:,:), ALLOCATABLE :: vin !!!!!!! Variables ! di1, di2, di3: length of dimensions of the 3D matrix of values ! ds[1-4]: length of dimensions of matrix with the slices ! maxNgridsin: maximum number of grid points from the 3D matrix in any slice ! varin: 3D RK variable to be used ! Ngridsin: number of grids from 3D RK matrix for each slice ! gridsin: coordinates of grids of the 3D RK matrix B to matrix of slices ! percentages: weights as percentages of space of 3D RK matrix for each slice !!!!! ! Available spatial statistics to compute inside each slice using values from 3D RK matrix ! 'min': minimum value ! 'max': maximum value ! 'mean': space weighted mean value ! 'mean2': space weighted quadratic mean value ! 'stddev': space weighted standard deviation value ! 'median': median value ! 'count': percentage of the space of matrix A covered by each different value of matrix B ! varout: output statistical variable fname = 'multi_spaceweightstats_in3DRK3_slc3v4' ! Let's be efficient? varout = fillVal64 DO s1 =1, ds1 DO s2 =1, ds2 DO s3 =1, ds3 DO s4 =1, ds4 Nin = Ngridsin(s1,s2,s3,s4) IF (Nin > 1) THEN IF (ALLOCATED(gin)) DEALLOCATE(gin) ALLOCATE(gin(Nin,2)) IF (ALLOCATED(pin)) DEALLOCATE(pin) ALLOCATE(pin(Nin)) IF (ALLOCATED(vin)) DEALLOCATE(vin) ALLOCATE(vin(Nin,di3)) IF (ALLOCATED(svin)) DEALLOCATE(svin) ALLOCATE(svin(Nin)) gin = gridsin(s1,s2,s3,s4,1:Nin,:) pin = percentages(s1,s2,s3,s4,1:Nin) ! Getting the values DO iv=1, Nin i1 = gin(iv,1) i2 = gin(iv,2) vin(iv,:) = varin(i1,i2,:) END DO ! Computing along d3 DO i3=1, di3 minv = fillVal64 maxv = -fillVal64 meanv = zeroRK mean2v = zeroRK stdv = zeroRK minv = MINVAL(vin(:,i3)) maxv = MAXVAL(vin(:,i3)) meanv = SUM(vin(:,i3)*pin) mean2v = SUM(vin(:,i3)**2*pin) DO iv=1,Nin stdv = stdv + ( (meanv - vin(iv,i3))*pin(iv) )**2 END DO stdv = SQRT(stdv) svin = vin(:,i3) CALL SortR_K(svin, Nin) medv = svin(INT(Nin/2)) varout(s1,s2,s3,s4,i3,1) = minv varout(s1,s2,s3,s4,i3,2) = maxv varout(s1,s2,s3,s4,i3,3) = meanv varout(s1,s2,s3,s4,i3,4) = mean2v varout(s1,s2,s3,s4,i3,5) = stdv varout(s1,s2,s3,s4,i3,6) = medv varout(s1,s2,s3,s4,i3,7) = Nin*1. END DO ELSE IF (Nin == 1) THEN i1 = gridsin(s1,s2,s3,s4,1,1) i2 = gridsin(s1,s2,s3,s4,1,2) varout(s1,s2,s3,s4,:,1) = varin(i1,i2,:) varout(s1,s2,s3,s4,:,2) = varin(i1,i2,:) varout(s1,s2,s3,s4,:,3) = varin(i1,i2,:) varout(s1,s2,s3,s4,:,4) = varin(i1,i2,:)*varin(i1,i2,:) varout(s1,s2,s3,s4,:,5) = zeroRK varout(s1,s2,s3,s4,:,6) = varin(i1,i2,:) varout(s1,s2,s3,s4,:,7) = Nin*1. ELSE varout(s1,s2,s3,s4,:,1) = fillval64 varout(s1,s2,s3,s4,:,2) = fillval64 varout(s1,s2,s3,s4,:,3) = fillval64 varout(s1,s2,s3,s4,:,4) = fillval64 varout(s1,s2,s3,s4,:,5) = fillval64 varout(s1,s2,s3,s4,:,6) = fillval64 varout(s1,s2,s3,s4,:,7) = zeroRK END IF END DO END DO END DO END DO IF (ALLOCATED(gin)) DEALLOCATE(gin) IF (ALLOCATED(pin)) DEALLOCATE(pin) IF (ALLOCATED(vin)) DEALLOCATE(vin) IF (ALLOCATED(svin)) DEALLOCATE(svin) RETURN END SUBROUTINE multi_spaceweightstats_in3DRK3_slc3v4 SUBROUTINE multi_index_mat2DI(d1, d2, d12, mat, value, Nindices, indices) ! Subroutine to provide the indices of the different locations of a value inside a 2D integer matrix IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d12 INTEGER, DIMENSION(d1,d2), INTENT(in) :: mat INTEGER,INTENT(in) :: value INTEGER, INTENT(out) :: Nindices INTEGER, DIMENSION(2,d12), INTENT(out) :: indices ! Local INTEGER :: i,j INTEGER :: Ncounts1D, icount1D INTEGER, DIMENSION(d2) :: diffmat1D !!!!!!! Variables ! d1, d2: shape of the 2D matrix ! mat: 2D matrix ! value: value to be looking for ! Nindices: number of times value found within matrix ! indices: indices of the found values fname = 'multi_index_mat2DI' Nindices = 0 indices = 0 DO i=1, d1 diffmat1D = mat(i,:) - value IF (ANY(diffmat1D == 0)) THEN Ncounts1D = COUNT(diffmat1D == 0) icount1D = 0 DO j=1, d2 IF (diffmat1D(j) == 0) THEN Nindices = Nindices + 1 indices(1,Nindices) = i indices(2,Nindices) = j icount1D = icount1D + 1 IF (icount1D == Ncounts1D) EXIT END IF END DO END IF END DO END SUBROUTINE multi_index_mat2DI SUBROUTINE multi_index_mat3DI(d1, d2, d3, d123, mat, value, Nindices, indices) ! Subroutine to provide the indices of the different locations of a value inside a 3D integer matrix IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d123 INTEGER, DIMENSION(d1,d2,d3), INTENT(in) :: mat INTEGER, INTENT(in) :: value INTEGER, INTENT(out) :: Nindices INTEGER, DIMENSION(3,d123), INTENT(out) :: indices ! Local INTEGER :: i,j,k INTEGER :: Ncounts1D, icount1D INTEGER :: Ncounts2D, icount2D INTEGER, DIMENSION(d2,d3) :: diffmat2D INTEGER, DIMENSION(d3) :: diffmat1D !!!!!!! Variables ! d1, d2, d3: shape of the 3D matrix ! mat: 3D matrix ! value: value to be looking for ! Nindices: number of times value found within matrix ! indices: indices of the found values fname = 'multi_index_mat3DI' Nindices = 0 indices = 0 DO i=1, d1 diffmat2D = mat(i,:,:) - value IF (ANY(diffmat2D == 0)) THEN Ncounts2D = COUNT(diffmat2D == 0) icount2D = 0 DO j=1, d2 diffmat1D = mat(i,j,:) - value IF (ANY(diffmat1D == 0)) THEN Ncounts1D = COUNT(diffmat1D == 0) icount1D = 0 DO k=1, d3 IF (diffmat1D(k) == 0) THEN Nindices = Nindices + 1 indices(1,Nindices) = i indices(2,Nindices) = j indices(3,Nindices) = k icount1D = icount1D + 1 IF (icount1D == Ncounts1D) EXIT END IF END DO icount2D = icount2D + icount1D IF (icount2D == Ncounts2D) EXIT END IF END DO END IF END DO END SUBROUTINE multi_index_mat3DI SUBROUTINE multi_index_mat4DI(d1, d2, d3, d4, d1234, mat, value, Nindices, indices) ! Subroutine to provide the indices of the different locations of a value inside a 4D integer matrix IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d4, d1234 INTEGER, DIMENSION(d1,d2,d3,d4), INTENT(in) :: mat INTEGER, INTENT(in) :: value INTEGER, INTENT(out) :: Nindices INTEGER, DIMENSION(4,d1234), INTENT(out) :: indices ! Local INTEGER :: i,j,k,l INTEGER :: Ncounts1D, icount1D INTEGER :: Ncounts2D, icount2D INTEGER :: Ncounts3D, icount3D INTEGER, DIMENSION(d2,d3,d4) :: diffmat3D INTEGER, DIMENSION(d3,d4) :: diffmat2D INTEGER, DIMENSION(d4) :: diffmat1D !!!!!!! Variables ! d1, d2, d3, d4: shape of the 4D matrix ! mat: 4D matrix ! value: value to be looking for ! Nindices: number of times value found within matrix ! indices: indices of the found values fname = 'multi_index_mat4DI' Nindices = 0 indices = 0 DO i=1, d1 diffmat3D = mat(i,:,:,:) - value IF (ANY(diffmat3D == 0)) THEN Ncounts3D = COUNT(diffmat3D == 0) icount3D = 0 DO j=1, d2 diffmat2D = mat(i,j,:,:) - value IF (ANY(diffmat2D == 0)) THEN Ncounts2D = COUNT(diffmat2D == 0) icount2D = 0 DO k=1, d3 diffmat1D = mat(i,j,k,:) - value IF (ANY(diffmat1D == 0)) THEN Ncounts1D = COUNT(diffmat1D == 0) icount1D = 0 DO l=1, d4 IF (diffmat1D(l) == 0) THEN Nindices = Nindices + 1 indices(1,Nindices) = i indices(2,Nindices) = j indices(3,Nindices) = k indices(4,Nindices) = l icount1D = icount1D + 1 IF (icount1D == Ncounts1D) EXIT END IF END DO icount2D = icount2D + icount1D IF (icount2D == Ncounts2D) EXIT END IF END DO icount3D = icount3D + icount1D IF (icount3D == Ncounts3D) EXIT END IF END DO END IF END DO END SUBROUTINE multi_index_mat4DI SUBROUTINE multi_index_mat2DRK(d1, d2, d12, mat, value, Nindices, indices) ! Subroutine to provide the indices of the different locations of a value inside a 2D RK matrix IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d12 REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: mat REAL(r_k),INTENT(in) :: value INTEGER, INTENT(out) :: Nindices INTEGER, DIMENSION(2,d12), INTENT(out) :: indices ! Local INTEGER :: i,j INTEGER :: Ncounts1D, icount1D REAL(r_k), DIMENSION(d2) :: diffmat1D !!!!!!! Variables ! d1, d2: shape of the 2D matrix ! mat: 2D matrix ! value: value to be looking for ! Nindices: number of times value found within matrix ! indices: indices of the found values fname = 'multi_index_mat2DRK' Nindices = 0 indices = 0 DO i=1, d1 diffmat1D = mat(i,:) - value IF (ANY(diffmat1D == zeroRK)) THEN Ncounts1D = COUNT(diffmat1D == zeroRK) icount1D = 0 DO j=1, d2 IF (diffmat1D(j) == zeroRK) THEN Nindices = Nindices + 1 indices(1,Nindices) = i indices(2,Nindices) = j icount1D = icount1D + 1 IF (icount1D == Ncounts1D) EXIT END IF END DO END IF END DO END SUBROUTINE multi_index_mat2DRK SUBROUTINE multi_index_mat3DRK(d1, d2, d3, d123, mat, value, Nindices, indices) ! Subroutine to provide the indices of the different locations of a value inside a 3D RK matrix IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d123 REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in) :: mat REAL(r_k),INTENT(in) :: value INTEGER, INTENT(out) :: Nindices INTEGER, DIMENSION(3,d123), INTENT(out) :: indices ! Local INTEGER :: i,j,k INTEGER :: Ncounts1D, icount1D INTEGER :: Ncounts2D, icount2D REAL(r_k), DIMENSION(d2,d3) :: diffmat2D REAL(r_k), DIMENSION(d3) :: diffmat1D !!!!!!! Variables ! d1, d2, d3: shape of the 3D matrix ! mat: 3D matrix ! value: value to be looking for ! Nindices: number of times value found within matrix ! indices: indices of the found values fname = 'multi_index_mat3DRK' Nindices = 0 indices = 0 DO i=1, d1 diffmat2D = mat(i,:,:) - value IF (ANY(diffmat2D == zeroRK)) THEN Ncounts2D = COUNT(diffmat2D == zeroRK) icount2D = 0 DO j=1, d2 diffmat1D = mat(i,j,:) - value IF (ANY(diffmat1D == zeroRK)) THEN Ncounts1D = COUNT(diffmat1D == zeroRK) icount1D = 0 DO k=1, d3 IF (diffmat1D(k) == zeroRK) THEN Nindices = Nindices + 1 indices(1,Nindices) = i indices(2,Nindices) = j indices(3,Nindices) = k icount1D = icount1D + 1 IF (icount1D == Ncounts1D) EXIT END IF END DO icount2D = icount2D + icount1D IF (icount2D == Ncounts2D) EXIT END IF END DO END IF END DO END SUBROUTINE multi_index_mat3DRK SUBROUTINE multi_index_mat4DRK(d1, d2, d3, d4, d1234, mat, value, Nindices, indices) ! Subroutine to provide the indices of the different locations of a value inside a 4D RK matrix IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, d3, d4, d1234 REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in) :: mat REAL(r_k),INTENT(in) :: value INTEGER, INTENT(out) :: Nindices INTEGER, DIMENSION(4,d1234), INTENT(out) :: indices ! Local INTEGER :: i,j,k,l INTEGER :: Ncounts1D, icount1D INTEGER :: Ncounts2D, icount2D INTEGER :: Ncounts3D, icount3D REAL(r_k), DIMENSION(d2,d3,d4) :: diffmat3D REAL(r_k), DIMENSION(d3,d4) :: diffmat2D REAL(r_k), DIMENSION(d4) :: diffmat1D !!!!!!! Variables ! d1, d2, d3, d4: shape of the 4D matrix ! mat: 4D matrix ! value: value to be looking for ! Nindices: number of times value found within matrix ! indices: indices of the found values fname = 'multi_index_mat4DRK' Nindices = 0 indices = 0 DO i=1, d1 diffmat3D = mat(i,:,:,:) - value IF (ANY(diffmat3D == zeroRK)) THEN Ncounts3D = COUNT(diffmat3D == zeroRK) icount3D = 0 DO j=1, d2 diffmat2D = mat(i,j,:,:) - value IF (ANY(diffmat2D == zeroRK)) THEN Ncounts2D = COUNT(diffmat2D == zeroRK) icount2D = 0 DO k=1, d3 diffmat1D = mat(i,j,k,:) - value IF (ANY(diffmat1D == zeroRK)) THEN Ncounts1D = COUNT(diffmat1D == zeroRK) icount1D = 0 DO l=1, d4 IF (diffmat1D(l) == zeroRK) THEN Nindices = Nindices + 1 indices(1,Nindices) = i indices(2,Nindices) = j indices(3,Nindices) = k indices(4,Nindices) = l icount1D = icount1D + 1 IF (icount1D == Ncounts1D) EXIT END IF END DO icount2D = icount2D + icount1D IF (icount2D == Ncounts2D) EXIT END IF END DO icount3D = icount3D + icount1D IF (icount3D == Ncounts3D) EXIT END IF END DO END IF END DO END SUBROUTINE multi_index_mat4DRK SUBROUTINE coincident_list_2Dcoords(NpointsA, pointsA, NpointsB, pointsB, Npoints, points, inpA, & inpB) ! Subroutine to determine which 2D points of an A list are also found in a B list IMPLICIT NONE INTEGER, INTENT(in) :: NpointsA, NpointsB INTEGER, DIMENSION(NpointsA,2), INTENT(in) :: pointsA INTEGER, DIMENSION(NpointsB,2), INTENT(in) :: pointsB INTEGER, INTENT(out) :: Npoints INTEGER, DIMENSION(NpointsA,2), INTENT(out) :: points INTEGER, DIMENSION(NpointsA), INTENT(out) :: inpA, inpB ! Local INTEGER :: iA, iB !!!!!!! Variables ! NpointsA: Number of points of the list A ! pointsA: points of the list A ! NpointsB: Number of points of the list B ! pointsB: points of the list B ! Npoints: Number of coincident points ! points: coincident points ! inpA: coincident points list A ! inpB: coincident points list B fname = 'coincident_list_2Dcoords' Npoints = 0 points = 0 inpA = 0 inpB = 0 DO iA = 1, NpointsA DO iB = 1, NpointsB IF ( (pointsA(iA,1) == pointsB(iB,1)) .AND. (pointsA(iA,2) == pointsB(iB,2)) ) THEN Npoints = Npoints + 1 points(Npoints,1) = pointsA(iA,1) points(Npoints,2) = pointsA(iA,2) inpA(Npoints) = iA inpB(Npoints) = iB EXIT END IF END DO END DO END SUBROUTINE coincident_list_2Dcoords SUBROUTINE coincident_gridsin2D_old(dxA, dyA, dxyA, NpointsA, pointsA, dxB, dyB, dxyB, NpointsB, & pointsB, Npoints, points, inpointsA, inpointsB) ! Subroutine to determine which lists of 2D gridsin points of an A list are also found in a B list IMPLICIT NONE INTEGER, INTENT(in) :: dxA, dyA, dxyA INTEGER, INTENT(in) :: dxB, dyB, dxyB INTEGER, DIMENSION(dxA, dyA), INTENT(in) :: NpointsA INTEGER, DIMENSION(dxB, dyB), INTENT(in) :: NpointsB INTEGER, DIMENSION(dxA, dyA, dxyA, 2), INTENT(in) :: pointsA INTEGER, DIMENSION(dxB, dyB, dxyB, 2), INTENT(in) :: pointsB INTEGER, DIMENSION(dxA, dyA, dxB, dyB), INTENT(out) :: Npoints INTEGER, DIMENSION(dxA, dyA, dxB, dyB, dxyA, 2), & INTENT(out) :: points INTEGER, DIMENSION(dxA, dyA, dxB, dyB, dxyA), & INTENT(out) :: inpointsA INTEGER, DIMENSION(dxA, dyA, dxB, dyB, dxyA), & INTENT(out) :: inpointsB ! Local INTEGER :: ixA, iyA, ixB, iyB, iv, ii INTEGER :: NA, NB INTEGER, DIMENSION(dxyA) :: ptsA, ptsB INTEGER, DIMENSION(dxyA, 2) :: pts !!!!!!! Variables ! dxA, dyA: 2D shape of the list A ! NpointsA: 2D Number of points of the list A ! pointsA: points of the list A ! dxB, dyB: 2D shape of the list B ! NpointsB: 2D Number of points of the list B ! pointsB: points of the list B ! Npoints: Number of coincident points ! points: coincident points ! inpointsA: coincident points list A ! inpointsB: coincident points list B fname = 'coincident_gridsin2D_old' Npoints = 0 points = 0 inpointsA = 0 inpointsB = 0 DO ixA=1, dxA DO iyA=1, dyA NA = NpointsA(ixA,iyA) DO ixB=1, dxB DO iyB=1, dyB NB = NpointsB(ixB,iyB) pts = -1 CALL coincident_list_2Dcoords(NA, pointsA(ixA,iyA,1:NA,:), NB, pointsB(ixB,iyB,1:NB,:), & Npoints(ixA,iyA,ixB,iyB), pts(1:NA,:), ptsA, ptsB) DO iv = 1, Npoints(ixA,iyA,ixB,iyB) points(ixA,iyA,ixB,iyB,iv,1) = pts(iv,1) points(ixA,iyA,ixB,iyB,iv,2) = pts(iv,2) inpointsA(ixA,iyA,ixB,iyB,iv) = ptsA(iv) inpointsB(ixA,iyA,ixB,iyB,iv) = ptsB(iv) END DO END DO END DO END DO END DO END SUBROUTINE coincident_gridsin2D_old SUBROUTINE coincident_gridsin2D(dxA, dyA, dxyA, NpointsA, pointsA, dxB, dyB, dxyB, NpointsB, & pointsB, Npoints, points, inpointsA, inpointsB) ! Subroutine to determine which lists of 2D gridsin points of an A list are also found in a B list IMPLICIT NONE INTEGER, INTENT(in) :: dxA, dyA, dxyA INTEGER, INTENT(in) :: dxB, dyB, dxyB INTEGER, DIMENSION(dxA, dyA), INTENT(in) :: NpointsA INTEGER, DIMENSION(dxB, dyB), INTENT(in) :: NpointsB INTEGER, DIMENSION(dxA, dyA, dxyA, 2), INTENT(in) :: pointsA INTEGER, DIMENSION(dxB, dyB, dxyB, 2), INTENT(in) :: pointsB INTEGER, DIMENSION(dxA, dyA, dxB, dyB), INTENT(out) :: Npoints INTEGER, DIMENSION(dxA, dyA, dxB, dyB, dxyA, 2), & INTENT(out) :: points INTEGER, DIMENSION(dxA, dyA, dxB, dyB, dxyA), & INTENT(out) :: inpointsA INTEGER, DIMENSION(dxA, dyA, dxB, dyB, dxyA), & INTENT(out) :: inpointsB ! Local INTEGER :: ixA, iyA, ixB, iyB, iv, iv1, iv2, iiv INTEGER :: NA, NB INTEGER, DIMENSION(dxyA) :: ptsA, ptsB INTEGER, DIMENSION(dxyA, 2) :: pts LOGICAL :: found !!!!!!! Variables ! dxA, dyA: 2D shape of the list A ! NpointsA: 2D Number of points of the list A ! pointsA: points of the list A ! dxB, dyB: 2D shape of the list B ! NpointsB: 2D Number of points of the list B ! pointsB: points of the list B ! Npoints: Number of coincident points ! points: coincident points ! inpointsA: coincident points list A ! inpointsB: coincident points list B fname = 'coincident_gridsin2D' Npoints = 0 points = 0 inpointsA = 0 inpointsB = 0 DO ixA=1, dxA DO iyA=1, dyA NA = NpointsA(ixA,iyA) DO ixB=1, dxB DO iyB=1, dyB NB = NpointsB(ixB,iyB) iv = 0 DO iv1=1, NA DO iv2=1, NB IF ( (pointsA(ixA,iyA,iv1,1) == pointsB(ixB,iyB,iv2,1)) .AND. & (pointsA(ixA,iyA,iv1,2) == pointsB(ixB,iyB,iv2,2)) ) THEN ! need to avoid double count (in ix[X] == 1) found = .FALSE. DO iiv = 1, iv IF ( (pointsA(ixA,iyA,iv1,1) == points(ixA,iyB,ixB,iyB,iiv,1)) .AND. & (pointsA(ixA,iyA,iv1,2) == points(ixA,iyB,ixB,iyB,iiv,2)) ) THEN found = .TRUE. END IF END DO IF (found) CYCLE iv = iv + 1 points(ixA,iyA,ixB,iyB,iv,1) = pointsA(ixA,iyA,iv1,1) points(ixA,iyA,ixB,iyB,iv,2) = pointsA(ixA,iyA,iv1,2) inpointsA(ixA,iyA,ixB,iyB,iv) = iv1 inpointsB(ixA,iyA,ixB,iyB,iv) = iv2 END IF END DO END DO Npoints(ixA,iyA,ixB,iyB) = iv END DO END DO END DO END DO END SUBROUTINE coincident_gridsin2D SUBROUTINE multi_spaceweightcount_in3DRK3_slc3v3(varin, Ngridsin, gridsin, percentages, Ncat, & categories, varout, di1, di2, di3, ds1, ds2, ds3, maxNgridsin) ! Subroutine to compute an spatial percentage of coverture of categories from a 3D RK matrix using ! 3rd dimension as running one into a matrix of 3-variables slices of rank 3 using spatial weights IMPLICIT NONE INTEGER, INTENT(in) :: di1, di2, di3, ds1, ds2, ds3 INTEGER, INTENT(in) :: Ncat, maxNgridsin INTEGER, DIMENSION(ds1,ds2,ds3), INTENT(in) :: Ngridsin INTEGER, DIMENSION(ds1,ds2,ds3,maxNgridsin,2), & INTENT(in) :: gridsin INTEGER, DIMENSION(di1,di2,di3), INTENT(in) :: varin REAL(r_k), INTENT(in), & DIMENSION(ds1,ds2,ds3,maxNgridsin) :: percentages INTEGER, DIMENSION(Ncat), INTENT(in) :: categories REAL(r_k), DIMENSION(ds1,ds2,ds3,di3,Ncat), & INTENT(out) :: varout ! Local INTEGER :: i1, i2, i3, s1, s2, s3, iv INTEGER :: ii3, ss1, ss2, ss3, Nin REAL(r_k), DIMENSION(:), ALLOCATABLE :: pin INTEGER, DIMENSION(:,:), ALLOCATABLE :: gin REAL(r_k), DIMENSION(:), ALLOCATABLE :: svin INTEGER, DIMENSION(:,:), ALLOCATABLE :: vin !!!!!!! Variables ! di1, di2, di3: length of dimensions of the 3D matrix of values ! ds[1-3]: length of dimensions of matrix with the slices ! maxNgridsin: maximum number of grid points from the 3D matrix in any slice ! varin: 3D RK variable to be used ! Ngridsin: number of grids from 3D RK matrix for each slice ! gridsin: coordinates of grids of the 3D RK matrix B to matrix of slices ! percentages: weights as percentages of space of 3D RK matrix for each slice ! Ncat: number of categories ! categories: category values to classify with ! varout: percentage of each category at each point along d3 fname = 'multi_spaceweightcount_in3DRK3_slc3v3' ! Let's be efficient? varout = fillVal64 DO s1 =1, ds1 DO s2 =1, ds2 DO s3 =1, ds3 Nin = Ngridsin(s1,s2,s3) ! Computing along d3 IF (Nin >= 1) THEN IF (ALLOCATED(gin)) DEALLOCATE(gin) ALLOCATE(gin(Nin,2)) IF (ALLOCATED(pin)) DEALLOCATE(pin) ALLOCATE(pin(Nin)) IF (ALLOCATED(vin)) DEALLOCATE(vin) ALLOCATE(vin(Nin,di3)) IF (ALLOCATED(svin)) DEALLOCATE(svin) ALLOCATE(svin(Nin)) gin = gridsin(s1,s2,s3,1:Nin,:) pin = percentages(s1,s2,s3,1:Nin) ! Getting the values DO iv=1, Nin i1 = gin(iv,1) i2 = gin(iv,2) vin(iv,:) = varin(i1,i2,:) END DO DO i3=1, di3 CALL spaceweight_Icount(Ncat, categories, vin(:,i3), Nin, pin, varout(s1,s2,s3,i3,:)) END DO ELSE varout(s1,s2,s3,:,1) = fillval64 END IF END DO END DO END DO IF (ALLOCATED(gin)) DEALLOCATE(gin) IF (ALLOCATED(pin)) DEALLOCATE(pin) IF (ALLOCATED(vin)) DEALLOCATE(vin) IF (ALLOCATED(svin)) DEALLOCATE(svin) RETURN END SUBROUTINE multi_spaceweightcount_in3DRK3_slc3v3 SUBROUTINE crossingpoints_polys(NvertexA, NvertexB, NvertexAB, polyA, polyB, Ncross, crossA, crossB,& crosspoints) ! Subroutine to provide the crossing points between two polygons IMPLICIT NONE INTEGER, INTENT(in) :: NvertexA, NvertexB, NvertexAB REAL(r_k), DIMENSION(NvertexA,2), INTENT(in) :: polyA REAL(r_k), DIMENSION(NvertexB,2), INTENT(in) :: polyB INTEGER, INTENT(out) :: Ncross INTEGER, DIMENSION(NvertexAB), INTENT(out) :: crossA, crossB REAL(r_k), DIMENSION(NvertexAB,2), INTENT(out) :: crosspoints ! Local INTEGER :: iA, iB, icoin, ii REAL(r_k), DIMENSION(2,2) :: face1, face2 INTEGER :: intersect REAL(r_k), DIMENSION(2) :: intersectpt !!!!!!! variables ! NvertexA: number of vertexs polygon A ! NvertexB: number of vertexs polygon B ! polyA: pairs of coordinates for the polygon A (clockwise) ! polyB: pairs of coordinates for the polygon B (clockwise) ! Ncross: number of crossing vertices ! crossA: vertex of polygon A where there is a crossing ! crossB: vertex of polygon B where there is a crossing ! crosspoints: coordinates of the crossing points fname = 'crossingpoints_polys' Ncross = 0 crossA = 0 crossB = 0 crosspoints = 0. ! Looping DO iA=1, NvertexA-1 face1(1,:) = polyA(iA,:) face1(2,:) = polyA(iA+1,:) DO iB=1, NvertexB-1 face2(1,:) = polyB(iB,:) face2(2,:) = polyB(iB+1,:) CALL intersectfaces(face1, face2, intersect, intersectpt) IF (intersect /= 0) THEN Ncross = Ncross + 1 crossA(Ncross) = iA crossB(Ncross) = iB crosspoints(Ncross,:) = intersectpt END IF END DO ! Last face B iB = NvertexB face2(1,:) = polyB(iB,:) face2(2,:) = polyB(1,:) CALL intersectfaces(face1, face2, intersect, intersectpt) IF (intersect /= 0) THEN Ncross = Ncross + 1 crossA(Ncross) = iA crossB(Ncross) = iB crosspoints(Ncross,:) = intersectpt END IF END DO ! Last face A iA = NvertexA face1(1,:) = polyA(iA,:) face1(2,:) = polyA(1,:) DO iB=1, NvertexB-1 face2(1,:) = polyB(iB,:) face2(2,:) = polyB(iB+1,:) CALL intersectfaces(face1, face2, intersect, intersectpt) IF (intersect /= 0) THEN Ncross = Ncross + 1 crossA(Ncross) = iA crossB(Ncross) = iB crosspoints(Ncross,:) = intersectpt END IF END DO ! Last face B iB = NvertexB face2(1,:) = polyB(iB,:) face2(2,:) = polyB(1,:) CALL intersectfaces(face1, face2, intersect, intersectpt) IF (intersect /= 0) THEN Ncross = Ncross + 1 crossA(Ncross) = iA crossB(Ncross) = iB crosspoints(Ncross,:) = intersectpt END IF END SUBROUTINE crossingpoints_polys SUBROUTINE gradient2D_1o(dx, dy, var, dsx, dsy, grad) ! Subroutine to compute the first order gradient of a 2D field ! FROM: From: https://en.wikipedia.org/wiki/Finite_difference_coefficient IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: var, dsx, dsy REAL(r_k), DIMENSION(dx,dy,2), INTENT(out) :: grad ! Local INTEGER :: i, j !!!!!!! Variables ! dx, dy: shape of the 2D field ! var: variable ! dsx, dsy: matrices of distances betweeen grid points along each axis ! grad: gradient fname = 'gradient2D_1o' grad = fillval64 DO i=1, dx-1 DO j=1, dy-1 grad(i,j,:) = (/ (var(i+1,j)-var(i,j))/dsx(i,j), (var(i,j+1)-var(i,j))/dsy(i,j) /) END DO END DO END SUBROUTINE gradient2D_1o SUBROUTINE gradient3D_1o(dx, dy, dz, var, dsx, dsy, dsz, grad) ! Subroutine to compute the first order gradient of a 3D field ! FROM: From: https://en.wikipedia.org/wiki/Finite_difference_coefficient IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, dz REAL(r_k), DIMENSION(dx,dy,dz), INTENT(in) :: var REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: dsx, dsy REAL(r_k), DIMENSION(dz), INTENT(in) :: dsz REAL(r_k), DIMENSION(dx,dy,dz,3), INTENT(out) :: grad ! Local INTEGER :: i, j, k !!!!!!! Variables ! dx, dy, dz: shape of the 3D field ! var: variable ! dsx, dsy, dsz: matrices of distances betweeen grid points along each axis ! grad: gradient fname = 'gradient3D_1o' grad = fillval64 DO i=1, dx-1 DO j=1, dy-1 DO k=1, dz-1 grad(i,j,k,:) = (/ (var(i+1,j,k)-var(i,j,k))/dsx(i,j), (var(i,j+1,k)-var(i,j,k))/dsy(i,j), & (var(i,j,k+1)-var(i,j,k))/dsz(k) /) END DO END DO END DO END SUBROUTINE gradient3D_1o SUBROUTINE gradient2D_c1o(dx, dy, var, dsx, dsy, grad) ! Subroutine to compute the first order centerd gradient of a 2D field IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: var, dsx, dsy REAL(r_k), DIMENSION(dx,dy,2), INTENT(out) :: grad ! Local INTEGER :: i, j REAL(r_k) :: d1x, d1y !!!!!!! Variables ! dx, dy: shape of the 2D field ! var: variable ! dsx, dsy: matrices of distances betweeen grid points along each axis ! grad: gradient fname = 'gradient2D_c1o' grad = fillval64 DO i=2, dx-1 DO j=2, dy-1 d1x = (0.5*var(i+1,j)-0.5*var(i-1,j))/dsx(i,j) d1y = (0.5*var(i,j+1)-0.5*var(i,j-1))/dsy(i,j) grad(i,j,:) = (/ d1x, d1y /) END DO END DO END SUBROUTINE gradient2D_c1o SUBROUTINE divergence2D_1o(dx, dy, dsx, dsy, xvar, yvar, diver) ! Subroutine to compute the first order divergence of a 2D vectorial field IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: xvar, yvar, dsx, dsy REAL(r_k), DIMENSION(dx,dy), INTENT(out) :: diver ! Local INTEGER :: i, j !!!!!!! Variables ! dx, dy: shape of the 2D field ! xvar, yvar: x-y components of vectorial variable ! dsx, dsy: matrices of distances betweeen grid points along each axis ! diver: divergence fname = 'divergence2D_1o' diver = fillval64 DO i=1, dx-1 DO j=1, dy-1 diver(i,j) = (xvar(i+1,j)-xvar(i,j))/dsx(i,j) + (yvar(i,j+1)-yvar(i,j))/dsy(i,j) END DO END DO END SUBROUTINE divergence2D_1o SUBROUTINE divergence2D_c1o(dx, dy, dsx, dsy, xvar, yvar, diver) ! Subroutine to compute the first order centerd divergence of a 2D vectorial field IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: xvar, yvar, dsx, dsy REAL(r_k), DIMENSION(dx,dy), INTENT(out) :: diver ! Local INTEGER :: i, j REAL(r_k) :: d1x, d1y !!!!!!! Variables ! dx, dy: shape of the 2D field ! xvar, yvar: x-y components of vectorial variable ! dsx, dsy: matrices of distances betweeen grid points along each axis ! diver: divergence fname = 'divergence2D_c1o' diver = fillval64 DO i=2, dx-1 DO j=2, dy-1 d1x = 0.5*(xvar(i+1,j)-0.5*xvar(i-1,j))/dsx(i,j) d1y = 0.5*(yvar(i,j+1)-0.5*yvar(i,j-1))/dsy(i,j) diver(i,j) = d1x + d1y END DO END DO END SUBROUTINE divergence2D_c1o SUBROUTINE curl2D_1o(dx, dy, dsx, dsy, xvar, yvar, curl) ! Subroutine to compute the first order curl of a 2D vectorial field IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: xvar, yvar, dsx, dsy REAL(r_k), DIMENSION(dx,dy), INTENT(out) :: curl ! Local INTEGER :: i, j !!!!!!! Variables ! dx, dy: shape of the 2D field ! xvar, yvar: x-y components of vectorial variable ! dsx, dsy: matrices of distances betweeen grid points along each axis ! curl: curl fname = 'curl2D_1o' curl = fillval64 DO i=1, dx-1 DO j=1, dy-1 curl(i,j) = (yvar(i+1,j)-yvar(i,j))/dsx(i,j) - (xvar(i,j+1)-xvar(i,j))/dsy(i,j) END DO END DO END SUBROUTINE curl2D_1o SUBROUTINE curl2D_c1o(dx, dy, dsx, dsy, xvar, yvar, curl) ! Subroutine to compute the first order centered curl of a 2D vectorial field IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: xvar, yvar, dsx, dsy REAL(r_k), DIMENSION(dx,dy), INTENT(out) :: curl ! Local INTEGER :: i, j REAL(r_k) :: d1x, d1y !!!!!!! Variables ! dx, dy: shape of the 2D field ! xvar, yvar: x-y components of vectorial variable ! dsx, dsy: matrices of distances betweeen grid points along each axis ! curl: curl fname = 'curl2D_c1o' curl = fillval64 DO i=1, dx-1 DO j=1, dy-1 d1x = 0.5*(yvar(i+1,j)-0.5*yvar(i-1,j))/dsx(i,j) d1y = 0.5*(xvar(i,j+1)-0.5*xvar(i,j-1))/dsy(i,j) curl(i,j) = d1x - d1y END DO END DO END SUBROUTINE curl2D_c1o SUBROUTINE lap2D_1o(dx, dy, dsx, dsy, var, lap) ! Subroutine to compute the first order laplacian of a 2D vectorial field IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: var, dsx, dsy REAL(r_k), DIMENSION(dx,dy), INTENT(out) :: lap ! Local INTEGER :: i, j REAL(r_k) :: d2x, d2y !!!!!!! Variables ! dx, dy: shape of the 2D field ! xvar, yvar: x-y components of vectorial variable ! dsx, dsy: matrices of distances betweeen grid points along each axis ! lap: laplacian fname = 'lap2D_1o' lap = fillval64 DO i=1, dx-2 DO j=1, dy-2 d2x = (1.*var(i,j)-2.*var(i+1,j)+1.*var(i+2,j))/(dsx(i,j)**2) d2y = (1.*var(i,j)-2.*var(i,j+1)+1.*var(i,j+2))/(dsy(i,j)**2) lap(i,j) = d2x+d2y END DO END DO END SUBROUTINE lap2D_1o SUBROUTINE lap2D_c1o(dx, dy, dsx, dsy, var, lap) ! Subroutine to compute the first order centered laplacian of a 2D vectorial field IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: var, dsx, dsy REAL(r_k), DIMENSION(dx,dy), INTENT(out) :: lap ! Local INTEGER :: i, j REAL(r_k) :: d2x, d2y !!!!!!! Variables ! dx, dy: shape of the 2D field ! xvar, yvar: x-y components of vectorial variable ! dsx, dsy: matrices of distances betweeen grid points along each axis ! lap: laplacian fname = 'lap2D_c1o' lap = fillval64 DO i=2, dx-1 DO j=2, dy-1 d2x = (1.*var(i-1,j)-2.*var(i,j)+1.*var(i+1,j))/(dsx(i,j)**2) d2y = (1.*var(i,j-1)-2.*var(i,j)+1.*var(i,j+1))/(dsy(i,j)**2) lap(i,j) = d2x+d2y END DO END DO END SUBROUTINE lap2D_c1o REAL(r_k) FUNCTION vecmodule3D(vec) ! Function to compute the module of a 3D vector IMPLICIT NONE REAL(r_k), DIMENSION(3), INTENT(in) :: vec fname = 'vecmodule3D' vecmodule3D = SQRT(vec(1)*vec(1) + vec(2)*vec(2) + vec(3)*vec(3)) END FUNCTION vecmodule3D SUBROUTINE matmodule3D(d1, d2, d3, mat, matmodule) ! Subroutine to compute the module of a 3D matrix with 3 components IMPLICIT NONE INTEGER, INTENT(in) :: d1,d2,d3 REAL(r_k), DIMENSION(d1,d2,d3,3), INTENT(in) :: mat REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out) :: matmodule fname = 'matmodule3D' matmodule = SQRT(mat(:,:,:,1)*mat(:,:,:,1) + mat(:,:,:,2)*mat(:,:,:,2) + mat(:,:,:,3)*mat(:,:,:,3)) END SUBROUTINE matmodule3D SUBROUTINE deformation3D(dx, dy, dz, vargrad, uagrad, vagrad, def) ! Subroutine to compute the deformation of a 3D field IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, dz REAL(r_k), DIMENSION(dx,dy,dz,3), INTENT(in) :: vargrad, uagrad, vagrad REAL(r_k), DIMENSION(dx,dy,dz,3), INTENT(out) :: def ! Local INTEGER :: i, j, k !!!!!!! Variables ! dx, dy, dz: shape of the 3D field ! vargrad: gradient of the variable ! uagrad, vagrad: gradient of x/y components of the wind field ! def: deformation fname = 'deformation3D' def = fillval64 def(:,:,:,1) = -uagrad(:,:,:,1)*vargrad(:,:,:,1) -vagrad(:,:,:,1)*vargrad(:,:,:,2) def(:,:,:,2) = -uagrad(:,:,:,2)*vargrad(:,:,:,1) -vagrad(:,:,:,2)*vargrad(:,:,:,2) def(:,:,:,3) = -uagrad(:,:,:,3)*vargrad(:,:,:,1) -vagrad(:,:,:,3)*vargrad(:,:,:,2) END SUBROUTINE deformation3D SUBROUTINE tilting3D(dx, dy, dz, vargrad, wagrad, tilt) ! Subroutine to compute the tilting of a 3D field IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, dz REAL(r_k), DIMENSION(dx,dy,dz,3), INTENT(in) :: vargrad, wagrad REAL(r_k), DIMENSION(dx,dy,dz,3), INTENT(out) :: tilt ! Local INTEGER :: i, j, k !!!!!!! Variables ! dx, dy, dz: shape of the 3D field ! vargrad: gradient of the variable ! wagrad: gradient of z component of the wind field ! tilt: tilting fname = 'tilting3D' tilt = fillval64 tilt(:,:,:,1) = -wagrad(:,:,:,1)*vargrad(:,:,:,3) tilt(:,:,:,2) = -wagrad(:,:,:,2)*vargrad(:,:,:,3) tilt(:,:,:,3) = -wagrad(:,:,:,3)*vargrad(:,:,:,3) END SUBROUTINE tilting3D SUBROUTINE deltat3D(dx, dy, dz, dt, var, ddt, dtkind, vardt) ! Subroutine to compute the temporal derivative of a 3D field IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, dz, dt REAL(r_k), INTENT(in) :: ddt REAL(r_k), DIMENSION(dx,dy,dz,dt), INTENT(in) :: var CHARACTER(LEN=*), INTENT(in) :: dtkind REAL(r_k), DIMENSION(dx,dy,dz,dt), INTENT(out) :: vardt ! Local INTEGER :: i, j, k, it !!!!!!! Variables ! dx, dy, dz: shape of the 3D field ! dt: number of time-steps ! var: variable [X] ! ddt: delta t [s] ! dtkind: kind of delta t to compute ! 'backward': [X(it) - X(it-1)] / ddt ! 'centered': [X(it+1) - 2*X(it) + X(it-1)] / 2*ddt ! 'forward': [X(it+1) - X(it)] / ddt ! vardt: temporal derivative of var [X s-1] fname = 'deltat3D' vardt = fillval64 IF (TRIM(dtkind) == 'backward') THEN DO it=2, dt vardt(:,:,:,it) = (var(:,:,:,it) - var(:,:,:,it-1))/ddt END DO ELSE IF (TRIM(dtkind) == 'centered') THEN DO it=2, dt-1 vardt(:,:,:,it) = (var(:,:,:,it+1) - 2*var(:,:,:,it) + var(:,:,:,it-1))/(2.*ddt) END DO ELSE IF (TRIM(dtkind) == 'forward') THEN DO it=1, dt-1 vardt(:,:,:,it) = (var(:,:,:,it+1) - var(:,:,:,it))/ddt END DO ELSE msg = "kind of delta t '" // TRIM(dtkind) // "' not ready !!" CALL ErrMsg(msg, fname, -1) END IF END SUBROUTINE deltat3D SUBROUTINE HistogramRK(Nvals, vals, missval, Nbins, Nbins1, bins, hist) ! Subroutine to provide the histogram of a series of RK values ! Histogram will be provided with Nbins+1 ! 1st hist: couns(vals= bins(ib-1)) .AND. (vals(i) < bins(ib)) ) THEN hist(ib) = hist(ib) + 1 CYCLE END IF END DO IF (vals(i) >= bins(Nbins)) hist(Nbins+1) = hist(Nbins+1) + 1 END DO RETURN END SUBROUTINE HistogramRK SUBROUTINE Histogram2D1_RK(d1, d2, vals, missval, Nbins, Nbins1, bins, hist) ! Subroutine to provide the histogram of a 2D matrix RK values along the first dimension IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, Nbins, Nbins1 REAL(r_k),INTENT(in) :: missval REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: vals REAL(r_k), DIMENSION(Nbins), INTENT(in) :: bins INTEGER, DIMENSION(d2, Nbins1), INTENT(out) :: hist ! Local INTEGER :: i, j !!!!!!! Variables ! d1, d2: dimension of values ! vals: values ! missval: value for the missing value ! Nbins: number of bins ! bins: bins of the histogram ! hist: histogram fname = 'Histogram2D1_RK' hist = zeroRK DO i=1, d2 CALL HistogramRK(d1, vals(:,i), missval, Nbins, Nbins1, bins, hist(i,:)) END DO RETURN END SUBROUTINE Histogram2D1_RK SUBROUTINE get_multindices_lonlat(dx, dy, lons, lats, dvals, lonvals, latvals, indices) ! Subroutine to provide the nearest grid point of a series of longitudes and latitudes values IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, dvals REAL(r_k), DIMENSION(dx, dy), INTENT(in) :: lons, lats REAL(r_k), DIMENSION(dvals), INTENT(in) :: lonvals, latvals INTEGER, DIMENSION(2,dvals), INTENT(out) :: indices ! Local INTEGER :: iv, ip REAL(r_k) :: mindiff REAL(r_k), DIMENSION(dx,dy) :: diff !!!!!!! Variables ! dx, dy: shape of the longitude and latitude matrices where to find the indices ! lons, lats: matrices of longitudes and latitudes ! dvals: number of longitudes and latitudes values ! lonvals, latvals: longitudes and latitudes values for which search the indices ! indices: indices found fname = 'get_multindices_lonlat' indices = zeroRK DO iv=1, dvals diff = SQRT((lons-lonvals(iv))**2+(lats-latvals(iv))**2) mindiff = MINVAL(diff) CALL multi_index_mat2DRK(dx, dy, 1, diff, mindiff, ip, indices(:,iv)) END DO END SUBROUTINE get_multindices_lonlat END MODULE module_scientific