MODULE module_generic ! Module with generic functions !!!!!!! Subroutines/Functions ! continguos_homogene_zones: Subroutine to look for contiguous zones by looking by continuous grid points ! freeunit: provides the number of a free unit in which open a file ! from_coordlist_2DRKmatrix: Subroutine to construct a 2D RK matrix from a list of values accompaigned ! by a list of coordinates to find i,j grid-point coordinates by minimum distance ! from_ptlist_2DRKmatrix: Subroutine to construct a 2D RK matrix from a list of values accompaigned ! by a list of grid-point coordinates ! from_ptlist_2DRKNmatrix: Subroutine to construct N 2D RK matrix from a list of values accompaigned ! by a list of grid-point coordinates ! GetInNamelist: Subroutine to get a paramter from a namelistfile ! get_xyconlimits: Subroutine for getting the limits of contiguous values from a given point in a 2D matrix ! index_list_coordsI: Function to provide the index of a given coordinate within a list of integer coordinates ! Index1DArrayI: Function to provide the first index of a given value inside a 1D integer array ! Index1DArrayR: Function to provide the first index of a given value inside a 1D real array ! Index1DArrayR_K: Function to provide the first index of a given value inside a 1D real(r_k) array ! Index1DArrayL: Function to provide the first index of a given value inside a 1D boolean array ! Index2DArrayR: Function to provide the first index of a given value inside a 2D real array ! Index2DArrayR_K: Function to provide the first index of a given value inside a 2D real(r_k) array ! Nvalues_2DArrayI: Number of different values of a 2D integer array ! mat2DPosition: Function to provide the i, j indices of a given value inside a 2D matrix ! multi_Index1DArrayL: Subroutine to provide the indices of a given value inside a 1D boolean array ! numberTimes: Function to provide the number of times that a given set of characters happen within a string ! RangeI: Function to provide a range of d1 values from 'iniv' to 'endv', of integer values in a vector ! RangeR: Function to provide a range of d1 values from 'iniv' to 'endv', of real values in a vector ! RangeR_K: Function to provide a range of d1 from 'iniv' to 'endv', of real(r_k) values in a vector ! stoprun: Subroutine to stop running and print a message ! vectorI_S: Function to transform a vector of integers to a string of characters ! vectorR_S: Function to transform a vector of reals to a string of characters ! xzones_homogenization: Subroutine to homogenize 2D contiguous zones along x-axis. Zones might be ! contiguous, but with different number assigned ! ! yzones_homogenization: Subroutine to homogenize 2D contiguous zones along y-axis. Zones might be ! contiguous, but with different number assigned ! USE module_definitions USE module_basic CONTAINS SUBROUTINE Nvalues_2DArrayI(dx, dy, dxy, mat2DI, Nvals, vals) ! Subroutine to give the number of different values of a 2D integer array IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy, dxy INTEGER, DIMENSION(dx,dy), INTENT(in) :: mat2DI INTEGER, INTENT(out) :: Nvals INTEGER, DIMENSION(dxy), INTENT(out) :: vals ! Local INTEGER :: i, j, ij !!!!!!! Variables ! dx, dy: size of the 2D space ! mat2DI: 2D integer matrix ! Nvals: number of different values ! vals: vector with the different values fname = 'Nvalues_2DArrayI' vals = 0 Nvals = 1 vals(1) = mat2DI(1,1) DO i=1,dx DO j=1,dy IF (Index1DArrayI(vals, Nvals, mat2DI(i,j)) == -1) THEN Nvals = Nvals + 1 vals(Nvals) = mat2DI(i,j) END IF END DO END DO RETURN END SUBROUTINE Nvalues_2DArrayI INTEGER FUNCTION index_list_coordsI(Ncoords, coords, icoord) ! Function to provide the index of a given coordinate within a list of integer coordinates IMPLICIT NONE INTEGER, INTENT(in) :: Ncoords INTEGER, DIMENSION(Ncoords,2), INTENT(in) :: coords INTEGER, DIMENSION(2), INTENT(in) :: icoord ! Local INTEGER, DIMENSION(Ncoords) :: dist INTEGER :: i,mindist INTEGER, DIMENSION(1) :: iloc !!!!!!! Variables ! Ncoords: number of coordinates in the list ! coords: list of coordinates ! icoord: coordinate to find fname = 'index_list_coordsI' dist = (coords(:,1)-icoord(1))**2+(coords(:,2)-icoord(2))**2 IF (ANY(dist == 0)) THEN iloc = MINLOC(dist) index_list_coordsI = iloc(1) ELSE index_list_coordsI = -1 END IF END FUNCTION index_list_coordsI INTEGER FUNCTION Index1DArrayI(array1D, d1, val) ! Function to provide the first index of a given value inside a 1D integer array IMPLICIT NONE INTEGER, INTENT(in) :: d1 INTEGER, INTENT(in) :: val INTEGER, DIMENSION(d1), INTENT(in) :: array1D ! Local INTEGER :: i fname = 'Index1DArrayI' Index1DArrayI = -1 DO i=1,d1 IF (array1d(i) == val) THEN Index1DArrayI = i EXIT END IF END DO END FUNCTION Index1DArrayI INTEGER FUNCTION Index1DArrayR(array1D, d1, val) ! Function to provide the first index of a given value inside a 1D real array IMPLICIT NONE INTEGER, INTENT(in) :: d1 REAL, INTENT(in) :: val REAL, DIMENSION(d1), INTENT(in) :: array1D ! Local INTEGER :: i fname = 'Index1DArrayR' Index1DArrayR = -1 DO i=1,d1 IF (array1d(i) == val) THEN Index1DArrayR = i EXIT END IF END DO END FUNCTION Index1DArrayR INTEGER FUNCTION Index1DArrayR_K(array1D, d1, val) ! Function to provide the first index of a given value inside a 1D real(r_k) array IMPLICIT NONE INTEGER, INTENT(in) :: d1 REAL(r_k), INTENT(in) :: val REAL(r_k), DIMENSION(d1), INTENT(in) :: array1D ! Local INTEGER :: i fname = 'Index1DArrayR_K' Index1DArrayR_K = -1 DO i=1,d1 IF (array1d(i) == val) THEN Index1DArrayR_K = i EXIT END IF END DO END FUNCTION Index1DArrayR_K INTEGER FUNCTION Index1DArrayL(array1D, d1, val) ! Function to provide the first index of a given value inside a 1D boolean array IMPLICIT NONE INTEGER, INTENT(in) :: d1 LOGICAL, INTENT(in) :: val LOGICAL, DIMENSION(d1), INTENT(in) :: array1D ! Local INTEGER :: i CHARACTER(LEN=50) :: fname fname = 'Index1DArrayL' Index1DArrayL = -1 DO i=1,d1 IF (array1d(i) .EQV. val) THEN Index1DArrayL = i EXIT END IF END DO END FUNCTION Index1DArrayL SUBROUTINE multi_Index1DArrayL(array1D, d1, val, Ntimes, pos) ! Subroutine to provide the indices of a given value inside a 1D boolean array IMPLICIT NONE INTEGER, INTENT(in) :: d1 LOGICAL, INTENT(in) :: val LOGICAL, DIMENSION(d1), INTENT(in) :: array1D INTEGER, INTENT(out) :: Ntimes INTEGER, DIMENSION(d1), INTENT(out) :: pos ! Local INTEGER :: i CHARACTER(LEN=50) :: fname !!!!!!! Variables ! array1D: 1D Array of values ! d1: length of array ! val: Value to look for ! Ntimes: Number of times val is found within array1D ! pos: positions of the values of val fname = 'multi_Index1DArrayL' Ntimes = 0 pos = -1 DO i=1,d1 IF (array1d(i) .EQV. val) THEN Ntimes = Ntimes + 1 pos(Ntimes) = i END IF END DO END SUBROUTINE multi_Index1DArrayL FUNCTION Index2DArrayR(array2D, d1, d2, val) ! Function to provide the first index of a given value inside a 2D real array IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2 REAL, INTENT(in) :: val REAL, DIMENSION(d1,d2), INTENT(in) :: array2D INTEGER, DIMENSION(2) :: Index2DArrayR ! Local INTEGER :: i, j fname = 'Index2DArrayR' Index2DArrayR = -1 DO i=1,d1 DO j=1,d2 IF (array2d(i,j) == val) THEN Index2DArrayR(1) = i Index2DArrayR(2) = j EXIT END IF END DO END DO END FUNCTION Index2DArrayR FUNCTION Index2DArrayR_K(array2D, d1, d2, val) ! Function to provide the first index of a given value inside a 2D real array IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2 REAL(r_k), INTENT(in) :: val REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: array2D INTEGER, DIMENSION(2) :: Index2DArrayR_K ! Local INTEGER :: i, j fname = 'Index2DArrayR_K' Index2DArrayR_K = -1 DO i=1,d1 DO j=1,d2 IF (array2d(i,j) == val) THEN Index2DArrayR_K(1) = i Index2DArrayR_K(2) = j EXIT END IF END DO END DO END FUNCTION Index2DArrayR_K INTEGER FUNCTION numberTimes(String, charv) ! Function to provide the number of times that a given set of characters happen within a string IMPLICIT NONE CHARACTER(LEN=*), INTENT(IN) :: String, charv ! Local INTEGER :: i, Lstring, Lcharv numberTimes = 0 Lstring = LEN_TRIM(String) Lcharv = LEN_TRIM(charv) DO i=1,Lstring - Lcharv IF (String(i:i+Lcharv-1) == TRIM(charv)) numberTimes = numberTimes + 1 END DO END FUNCTION numberTimes FUNCTION RangeI(d1, iniv, endv) ! Function to provide a range of d1 values from 'iniv' to 'endv', of integer values in a vector IMPLICIT NONE INTEGER, INTENT(in) :: d1, iniv, endv INTEGER, DIMENSION(d1) :: RangeI ! Local INTEGER :: i, intv fname = 'RangeI' intv = (endv - iniv) / (d1*1 - 1) RangeI(1) = iniv DO i=2,d1 RangeI(i) = RangeI(i-1) + intv END DO END FUNCTION RangeI FUNCTION RangeR(d1, iniv, endv) ! Function to provide a range of d1 from 'iniv' to 'endv', of real values in a vector IMPLICIT NONE INTEGER, INTENT(in) :: d1 REAL, INTENT(in) :: iniv, endv REAL, DIMENSION(d1) :: RangeR ! Local INTEGER :: i REAL :: intv fname = 'RangeR' intv = (endv - iniv) / (d1*1. - 1.) RangeR(1) = iniv DO i=2,d1 RangeR(i) = RangeR(i-1) + intv END DO END FUNCTION RangeR FUNCTION RangeR_K(d1, iniv, endv) ! Function to provide a range of d1 from 'iniv' to 'endv', of real(r_k) values in a vector IMPLICIT NONE INTEGER, INTENT(in) :: d1 REAL(r_k), INTENT(in) :: iniv, endv REAL(r_k), DIMENSION(d1) :: RangeR_K ! Local INTEGER :: i REAL(r_k) :: intv fname = 'RangeR_K' intv = (endv - iniv) / (d1*oneRK-oneRK) RangeR_K(1) = iniv DO i=2,d1 RangeR_K(i) = RangeR_K(i-1) + intv END DO END FUNCTION RangeR_K INTEGER FUNCTION freeunit() ! provides the number of a free unit in which open a file IMPLICIT NONE LOGICAL :: is_used is_used = .true. DO freeunit=10,100 INQUIRE(unit=freeunit, opened=is_used) IF (.not. is_used) EXIT END DO RETURN END FUNCTION freeunit SUBROUTINE GetInNamelist(namelistfile, param, kindparam, Ival, Rval, Lval, Sval) ! Subroutine to get a paramter from a namelistfile IMPLICIT NONE CHARACTER(LEN=*), INTENT(IN) :: namelistfile, param CHARACTER(LEN=1), INTENT(IN) :: kindparam INTEGER, OPTIONAL, INTENT(OUT) :: Ival REAL, OPTIONAL, INTENT(OUT) :: Rval LOGICAL, OPTIONAL, INTENT(OUT) :: Lval CHARACTER(LEN=200), OPTIONAL, INTENT(OUT) :: Sval ! Local INTEGER :: i, funit, ios INTEGER :: Lparam, posparam LOGICAL :: is_used CHARACTER(LEN=1000) :: line, message CHARACTER(LEN=200), DIMENSION(2) :: lvals CHARACTER(LEN=200) :: pval !!!!!!! Variables ! namelistfile: name of the namelist file ! param: parameter to get ! paramkind: kind of the parameter (I: Integer, L: boolean, R: Real, S: String) fname = 'GetInNamelist' ! Reading dimensions file and defining dimensions is_used = .true. DO funit=10,100 INQUIRE(unit=funit, opened=is_used) IF (.not. is_used) EXIT END DO OPEN(funit, FILE=TRIM(namelistfile), STATUS='old', FORM='formatted', IOSTAT=ios) IF ( ios /= 0 ) CALL stoprun(message, fname) Lparam = LEN_TRIM(param) DO i=1,10000 READ(funit,"(A200)",END=100)line posparam = INDEX(TRIM(line), TRIM(param)) IF (posparam /= 0) EXIT END DO 100 CONTINUE IF (posparam == 0) THEN message = "namelist '" // TRIM(namelistfile) // "' does not have parameter '" // TRIM(param) // & "' !!" CALL stoprun(message, fname) END IF CLOSE(UNIT=funit) CALL split(line, '=', 2, lvals) IF (kindparam /= 'S') THEN CALL RemoveNonNum(lvals(2), pval) END IF ! L. Fita, LMD. October 2015 ! Up to now, only getting scalar values kparam: SELECT CASE (kindparam) CASE ('I') Ival = StoI(pval) ! PRINT *,TRIM(param),'= ', Ival CASE ('L') Lval = StoL(pval) ! PRINT *,TRIM(param),'= ', Lval CASE ('R') Rval = StoR(pval) ! PRINT *,TRIM(param),'= ', Rval CASE ('S') Sval = lvals(2) CASE DEFAULT message = "type of parameter '" // kindparam // "' not ready !!" CALL stoprun(message, fname) END SELECT kparam END SUBROUTINE GetInNamelist SUBROUTINE stoprun(msg, fname) ! Subroutine to stop running and print a message IMPLICIT NONE CHARACTER(LEN=*), INTENT(IN) :: fname CHARACTER(LEN=*), INTENT(IN) :: msg ! local CHARACTER(LEN=50) :: errmsg, warnmsg errmsg = 'ERROR -- error -- ERROR -- error' PRINT *, TRIM(errmsg) PRINT *, ' ' // TRIM(fname) // ': ' // TRIM(msg) STOP END SUBROUTINE stoprun SUBROUTINE xzones_homogenization(dx, dy, inzones, outzones) ! Subroutine to homogenize 2D contiguous zones along x-axis, zones might be contiguous, but with ! different number assigned ! ! Here we have a 2D matrix of integers, with contiguous integer filled zones, zero outside any zone ! It might be, that within the same zone, might be lines which do not share the same integer ! 0 0 0 0 0 0 0 0 0 0 ! 0 1 1 0 0 0 1 1 0 0 ! 0 2 0 0 1 == > 0 1 0 0 2 ! 0 1 1 0 0 0 1 1 0 0 IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy INTEGER, DIMENSION(dx,dy), INTENT(in) :: inzones INTEGER, DIMENSION(dx,dy), INTENT(out) :: outzones ! Local INTEGER :: i,j,k INTEGER :: Nmaxzones, TOTzones LOGICAL :: assigned INTEGER, DIMENSION(dx) :: prevline INTEGER, DIMENSION(dy) :: Nxzones INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: zones !!!!!!! Variables ! dx, dy: Shape of the 2D space ! inzones: zones to homogenize ! outzones: zones homogenized fname = 'xzones_homogenization' ! Maximum possible number of zones Nmaxzones = INT((dx/2)*(dy/2)) ! Matrix with [i,j,Nzone,izone/ezone] IF (ALLOCATED(zones)) DEALLOCATE(zones) ALLOCATE(zones(dy,Nmaxzones,3)) zones = 0 Nxzones = 0 ! Getting beginning/end of y-bands DO j=1, dy k = 0 i = 1 IF (inzones(i,j) /= 0) THEN k = k + 1 zones(j,k,1) = i zones(j,k,3) = k END IF DO i=2, dx IF ( (inzones(i,j) /= 0) .AND. (inzones(i-1,j) == 0)) THEN k = k+1 zones(j,k,1) = i zones(j,k,3) = k ELSE IF ( (inzones(i-1,j) /= 0) .AND. (inzones(i,j) == 0)) THEN zones(j,k,2) = i-1 zones(j,k,3) = k END IF END DO IF (k > 0) THEN IF (zones(j,k,2) == 0) zones(j,k,2) = dx END IF Nxzones(j) = k END DO ! Homogenizing contigous zones outzones = 0 TOTzones = 0 j = 1 DO k = 1, Nxzones(j) TOTzones = TOTzones + 1 DO i=zones(j,k,1), zones(j,k,2) outzones(i,j) = TOTzones END DO END DO DO j=2, dy prevline = outzones(:,j-1) DO k = 1, Nxzones(j) assigned = .FALSE. DO i=zones(j,k,1), zones(j,k,2) IF (prevline(i) /= 0) THEN outzones(zones(j,k,1):zones(j,k,2),j) = prevline(i) assigned = .TRUE. EXIT END IF END DO IF (.NOT.assigned) THEN TOTzones = TOTzones + 1 DO i=zones(j,k,1), zones(j,k,2) outzones(i,j) = TOTzones END DO END IF END DO END DO IF (ALLOCATED(zones)) DEALLOCATE(zones) END SUBROUTINE xzones_homogenization SUBROUTINE yzones_homogenization(dx, dy, inzones, outzones) ! Subroutine to homogenize 2D contiguous zones along y-axis, zones might be contiguous, but with ! different number assigned ! ! Here we have a 2D matrix of integers, with contiguous integer filled zones, zero outside any zone ! It might be, that within the same zone, might be lines which do not share the same integer ! 0 0 0 0 0 0 0 0 0 0 ! 0 1 1 0 0 0 1 1 0 0 ! 0 2 0 0 1 == > 0 1 0 0 2 ! 0 1 1 0 0 0 1 1 0 0 IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy INTEGER, DIMENSION(dx,dy), INTENT(in) :: inzones INTEGER, DIMENSION(dx,dy), INTENT(out) :: outzones ! Local INTEGER :: i,j,k INTEGER :: Nmaxzones, TOTzones LOGICAL :: assigned INTEGER, DIMENSION(dy) :: prevline INTEGER, DIMENSION(dx) :: Nyzones INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: zones !!!!!!! Variables ! dx, dy: Shape of the 2D space ! inzones: zones to homogenize ! outzones: zones homogenized fname = 'yzones_homogenization' ! Maximum possible number of zones Nmaxzones = INT((dx/2)*(dy/2)) ! Matrix with [i,j,Nzone,izone/ezone] IF (ALLOCATED(zones)) DEALLOCATE(zones) ALLOCATE(zones(dx,Nmaxzones,3)) zones = 0 Nyzones = 0 ! Getting beginning/end of y-bands DO i=1, dx k = 0 j = 1 IF (inzones(i,j) /= 0) THEN k = k + 1 zones(i,k,1) = j zones(i,k,3) = k END IF DO j=2, dy IF ( (inzones(i,j) /= 0) .AND. (inzones(i,j-1) == 0)) THEN k = k+1 zones(i,k,1) = j zones(i,k,3) = k ELSE IF ( (inzones(i,j-1) /= 0) .AND. (inzones(i,j) == 0)) THEN zones(i,k,2) = j-1 zones(i,k,3) = k END IF END DO IF (k > 0) THEN IF (zones(i,k,2) == 0) zones(i,k,2) = dy END IF Nyzones(i) = k END DO ! Homogenizing contigous zones outzones = 0 TOTzones = 0 i = 1 DO k = 1, Nyzones(i) TOTzones = TOTzones + 1 DO j=zones(i,k,1), zones(i,k,2) outzones(i,j) = TOTzones END DO END DO DO i=2, dx prevline = outzones(i-1,:) DO k = 1, Nyzones(i) assigned = .FALSE. DO j=zones(i,k,1), zones(i,k,2) IF (prevline(j) /= 0) THEN outzones(i,zones(i,k,1):zones(i,k,2)) = prevline(j) assigned = .TRUE. EXIT END IF END DO IF (.NOT.assigned) THEN TOTzones = TOTzones + 1 DO j=zones(i,k,1), zones(i,k,2) outzones(i,j) = TOTzones END DO END IF END DO END DO IF (ALLOCATED(zones)) DEALLOCATE(zones) END SUBROUTINE yzones_homogenization CHARACTER(len=1000) FUNCTION vectorI_S(d1, vector) ! Function to transform a vector of integers to a string of characters IMPLICIT NONE INTEGER, INTENT(in) :: d1 INTEGER, DIMENSION(d1), INTENT(in) :: vector ! Local INTEGER :: iv CHARACTER(len=50) :: IS !!!!!!! Variables ! d1: length of the vector ! vector: values to transform fname = 'vectorI_S' vectorI_S = '' DO iv=1, d1 WRITE(IS, '(I50)')vector(iv) IF (iv == 1) THEN vectorI_S = TRIM(IS) ELSE vectorI_S = TRIM(vectorI_S) // ', ' // TRIM(IS) END IF END DO END FUNCTION vectorI_S CHARACTER(len=1000) FUNCTION vectorR_S(d1, vector) ! Function to transform a vector of reals to a string of characters IMPLICIT NONE INTEGER, INTENT(in) :: d1 REAL, DIMENSION(d1), INTENT(in) :: vector ! Local INTEGER :: iv CHARACTER(len=50) :: RS !!!!!!! Variables ! d1: length of the vector ! vector: values to transform fname = 'vectorR_S' vectorR_S = '' DO iv=1, d1 WRITE(RS, '(F50.25)')vector(iv) IF (iv == 1) THEN vectorR_S = TRIM(RS) ELSE vectorR_S = TRIM(vectorR_S) // ', ' // TRIM(RS) END IF END DO END FUNCTION vectorR_S SUBROUTINE continguos_homogene_zones(dx, dy, matvals, Nzones, contzones) ! Subroutine to look for contiguous zones by looking by continuous grid points IMPLICIT NONE INTEGER, INTENT(in) :: dx, dy INTEGER, DIMENSION(dx,dy), INTENT(in) :: matvals INTEGER, INTENT(out) :: Nzones INTEGER, DIMENSION(dx,dy), INTENT(out) :: contzones ! Local INTEGER :: i,j, k INTEGER :: ii, ei, ij, ej INTEGER :: Ncont, Nassigned, Ncontmin LOGICAL, DIMENSION(dx,dy) :: assigned, notdone INTEGER, DIMENSION(:), ALLOCATABLE :: pzones, allzones LOGICAL, DIMENSION(:), ALLOCATABLE :: passigns !!!!!!! Variables ! dx, dy: shape of the matrix ! matvals: matrix with the values ! contzones: homogeneous zones found fname = 'continguos_homogene_zones' ! Vector to keep track of all zone values IF (ALLOCATED(allzones)) DEALLOCATE(allzones) ALLOCATE(allzones(dx*dy/4)) allzones = 0 assigned = .FALSE. notdone = .TRUE. contzones = -1 Nzones = 0 DO i=1, dx DO j=1, dy ! First IF (matvals(i,j) /= 0 .AND. notdone(i,j)) THEN CALL get_xyconlimits(dx, dy, matvals, 0, i, j, ii, ei, ij, ej) ! Has any point of the rays already been assigned? ! Along x IF (ALLOCATED(pzones)) DEALLOCATE(pzones) ALLOCATE(pzones(ei-ii+1)) IF (ALLOCATED(passigns)) DEALLOCATE(passigns) ALLOCATE(passigns(ei-ii+1)) passigns = assigned(ii:ei,j) CALL multi_Index1DArrayL(passigns, ei-ii+1, .TRUE., Nassigned, pzones) IF (Nassigned /= 0) THEN Ncontmin = 10000000 DO k=1, Nassigned IF (contzones(ii+pzones(k)-1,j) < Ncontmin) Ncontmin = contzones(ii+pzones(k)-1,j) END DO ! If there is more than one assigned value change all values to the minimum one DO k=1, Nassigned IF (contzones(ii+pzones(k)-1,j) /= Ncontmin) THEN WHERE (contzones == contzones(ii+pzones(k)-1,j) .AND. contzones /= Ncontmin) contzones = Ncontmin END WHERE allzones(contzones(ii+pzones(k)-1,j)) = 0 END IF END DO Ncont = Ncontmin END IF ! Along y IF (ALLOCATED(pzones)) DEALLOCATE(pzones) ALLOCATE(pzones(ej-ij+1)) IF (ALLOCATED(passigns)) DEALLOCATE(passigns) ALLOCATE(passigns(ej-ij+1)) passigns = assigned(i,ij:ej) CALL multi_Index1DArrayL(passigns, ej-ij+1, .TRUE., Nassigned, pzones) IF (Nassigned /= 0) THEN Ncontmin = 10000000 DO k=1, Nassigned IF (contzones(i,ij+pzones(k)-1) < Ncontmin) Ncontmin = contzones(i,ij+pzones(k)-1) END DO ! If there is more than one assigned value change all values to the minimum one DO k=1, Nassigned IF (contzones(i,ij+pzones(k)-1) /= Ncontmin) THEN WHERE (contzones == contzones(i,ij+pzones(k)-1) .AND. contzones /= Ncontmin) contzones = Ncontmin END WHERE allzones(contzones(i,ij+pzones(k)-1)) = 0 END IF END DO Ncont = Ncontmin END IF IF (.NOT.assigned(i,j)) THEN Nzones = Nzones + 1 Ncont = Nzones allzones(Nzones) = 1 ELSE Ncont = contzones(i,j) END IF contzones(i,j) = Ncont contzones(ii:ei,j) = Ncont contzones(i,ij:ej) = Ncont notdone(i,j) = .FALSE. assigned(i,j) = .TRUE. assigned(ii:ei,j) = .TRUE. assigned(i,ij:ej) = .TRUE. END IF END DO END DO ! Using allzones to provide continuous assigned values Nzones = 0 DO k=1, dx*dy/4 IF (allzones(k) /= 0) THEN Nzones = Nzones + 1 IF (allzones(k) /= Nzones) THEN WHERE(contzones == allzones(k)) contzones = Nzones END WHERE END IF END IF END DO RETURN END SUBROUTINE continguos_homogene_zones SUBROUTINE get_xyconlimits(d1, d2, matv, NOval, i, j, ix, ex, iy, ey) ! Subroutine for getting the limits of contiguous values from a given point in a 2D matrix IMPLICIT NONE INTEGER, INTENT(in) :: d1, d2, i, j, NOval INTEGER, DIMENSION(d1,d2), INTENT(in) :: matv INTEGER, INTENT(out) :: ix, ex, iy, ey ! Local INTEGER :: i1, j1 LOGICAL :: found !!!!!!! Variables ! d1, d2: Shape of input 2D values ! i, j: grid point from which get the contiguous values ! NOval: value for grid points without data ! matv: 2D matrx values ! ix, ex, iy, ey: limits from a grid point fname = 'get_xyconlimits' ix = -1 ex = -1 iy = -1 ey = -1 ! Before i IF (i > 1) THEN ix = i found = .FALSE. DO i1=i-1,1,-1 IF (matv(i1,j) == NOval) THEN ix = i1 + 1 found = .TRUE. EXIT END IF END DO IF (.NOT.found) ix=1 ELSE ix = 1 END IF ! After i IF (i < d1) THEN ex = i found = .FALSE. DO i1=i+1,d1 IF (matv(i1,j) == NOval) THEN ex = i1 - 1 found = .TRUE. EXIT END IF END DO IF (.NOT.found) ex=d1 ELSE ex = d1 END IF ! Before j IF (j > 1) THEN iy = j found = .FALSE. DO j1=j-1,1,-1 IF (matv(i,j1) == NOval) THEN iy = j1 + 1 found = .TRUE. EXIT END IF END DO IF (.NOT.found) iy=1 ELSE iy = 1 END IF ! After j IF (j < d2) THEN ey = j found = .FALSE. DO j1=j+1,d2 IF (matv(i,j1) == NOval) THEN ey = j1 - 1 found = .TRUE. EXIT END IF END DO IF (.NOT.found) ey=d2 ELSE ey = d2 END IF RETURN END SUBROUTINE get_xyconlimits SUBROUTINE from_ptlist_2DRKmatrix(Npts, pts, Flike, vals, dx, dy, NOval, matrix) ! Subroutine to construct a 2D RK matrix from a list of values accompaigned by a list of grid-point ! coordinates IMPLICIT NONE INTEGER, INTENT(in) :: Npts, dx, dy LOGICAL, INTENT(in) :: Flike REAL(r_k), INTENT(in) :: NOval INTEGER, DIMENSION(Npts,2), INTENT(in) :: pts REAL(r_k), DIMENSION(Npts), INTENT(in) :: vals REAL(r_k), DIMENSION(dx,dy), INTENT(out) :: matrix ! Local INTEGER :: iv, i, j !!!!!!! Variables ! Npts: Number of values of the list ! pts: 2D matrix coordinates of the values ! Flike: whether input is Fortran-like (1: x, 2: y) or not (0-based; 1: y, 2: x) ! vals: list of values correspondant to the coordinates ! dx, dy: shape of the 2D matrix ! NOval: Value to assign when there is no coordinate ! matrix: resultant matrix fname = 'from_ptlist_2DRKmatrix' matrix = NOval IF (Flike) THEN DO iv=1, Npts i = pts(iv,1) j = pts(iv,2) matrix(i,j) = vals(iv) END DO ELSE DO iv=1, Npts i = pts(iv,2)+1 j = pts(iv,1)+1 matrix(i,j) = vals(iv) END DO END IF RETURN END SUBROUTINE from_ptlist_2DRKmatrix SUBROUTINE from_ptlist_2DRKNmatrix(Npts, pts, Flike, Nvals, vals, dx, dy, NOval, Nmatrix) ! Subroutine to construct N 2D RK matrix from a list of values accompaigned by a list of grid-point ! coordinates IMPLICIT NONE INTEGER, INTENT(in) :: Npts, dx, dy, Nvals LOGICAL, INTENT(in) :: Flike REAL(r_k), INTENT(in) :: NOval INTEGER, DIMENSION(Npts,2), INTENT(in) :: pts REAL(r_k), DIMENSION(Npts,Nvals), INTENT(in) :: vals REAL(r_k), DIMENSION(dx,dy,Nvals), INTENT(out) :: Nmatrix ! Local INTEGER :: iv, i, j, ic !!!!!!! Variables ! Npts: Number of values of the list ! pts: 2D matrix coordinates of the values ! Flike: whether input is Fortran-like (1: x, 2: y) or not (0-based; 1: y, 2: x) ! Nvals: number of different values ! vals: list of N-values correspondant to the coordinates ! dx, dy: shape of the 2D matrix ! NOval: Value to assign when there is no coordinate ! Nmatrix: resultant N-matrix fname = 'from_ptlist_2DRKNmatrix' Nmatrix = NOval IF (Flike) THEN DO iv=1, Npts i = pts(iv,1) j = pts(iv,2) Nmatrix(i,j,:) = vals(iv,:) END DO ELSE DO iv=1, Npts i = pts(iv,2)+1 j = pts(iv,1)+1 Nmatrix(i,j,:) = vals(iv,:) END DO END IF RETURN END SUBROUTINE from_ptlist_2DRKNmatrix SUBROUTINE from_coordlist_2DRKmatrix(Npts, coords, Flike, xcoord, ycoord, vals, dx, dy, NOval, & matrix) ! Subroutine to construct a 2D RK matrix from a list of values accompaigned by a list of coordinates ! to find i,j grid-point coordinates by minimum distance IMPLICIT NONE INTEGER, INTENT(in) :: Npts, dx, dy LOGICAL, INTENT(in) :: Flike REAL(r_k), INTENT(in) :: NOval REAL(r_k), DIMENSION(Npts,2), INTENT(in) :: coords REAL(r_k), DIMENSION(Npts), INTENT(in) :: vals REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: xcoord, ycoord REAL(r_k), DIMENSION(dx,dy), INTENT(out) :: matrix ! Local INTEGER :: iv, i, j REAL(r_k) :: xi, yi REAL(r_k), DIMENSION(dx,dy) :: diff INTEGER, DIMENSION(2) :: minpt !!!!!!! Variables ! Npts: Number of values of the list ! coords: 2D coordinates of the values ! Flike: whether input is Fortran-like (1: x, 2: y) or not (0-based; 1: y, 2: x) ! vals: list of values correspondant to the coordinates ! xcoord, ycoord: matrix of values correspondant to the each coordinate ! dx, dy: shape of the 2D matrix ! NOval: Value to assign when there is no coordinate ! matrix: resultant matrix fname = 'from_coordlist_2DRKmatrix' matrix = NOval IF (Flike) THEN DO iv=1, Npts xi = coords(iv,1) yi = coords(iv,2) diff = SQRT((xcoord-xi)**2+(ycoord-yi)**2) minpt = MINLOC(diff) matrix(minpt(1),minpt(2)) = vals(iv) END DO ELSE DO iv=1, Npts xi = coords(iv,2) yi = coords(iv,1) diff = SQRT((xcoord-xi)**2+(ycoord-yi)**2) minpt = MINLOC(diff) matrix(minpt(1),minpt(2)) = vals(iv) END DO END IF RETURN END SUBROUTINE from_coordlist_2DRKmatrix END MODULE module_generic