Context Navigation

module_generic.f90 @ 2736

Last change on this file since 2736 was 2727, checked in by lfita, 5 years ago
Fixing zero indices analysis
File size: 43.2 KB

Rev	Line
[715]	1	MODULE module_generic
	2	! Module with generic functions
[1608]	3
	4	!!!!!!! Subroutines/Functions
[2340]	5	! continguos_homogene_zones: Subroutine to look for contiguous zones by looking by continuous grid points
[1655]	6	! freeunit: provides the number of a free unit in which open a file
[2724]	7	! fill_matrix2DRK_winmat2D_list1D: Subroutine to fill a 2D RK matrix using a list of 1D indices from
	8	! another given 2D matrix
[2357]	9	! from_coordlist_2DRKmatrix: Subroutine to construct a 2D RK matrix from a list of values accompaigned
	10	! by a list of coordinates to find i,j grid-point coordinates by minimum distance
[2356]	11	! from_ptlist_2DRKmatrix: Subroutine to construct a 2D RK matrix from a list of values accompaigned
[2357]	12	! by a list of grid-point coordinates
	13	! from_ptlist_2DRKNmatrix: Subroutine to construct N 2D RK matrix from a list of values accompaigned
	14	! by a list of grid-point coordinates
[1655]	15	! GetInNamelist: Subroutine to get a paramter from a namelistfile
[2475]	16	! GDATE: Subroutine to compute the gregorian calendar date (year,month,day) given the julian date (JD)
[2340]	17	! get_xyconlimits: Subroutine for getting the limits of contiguous values from a given point in a 2D matrix
[1612]	18	! index_list_coordsI: Function to provide the index of a given coordinate within a list of integer coordinates
[1615]	19	! Index1DArrayI: Function to provide the first index of a given value inside a 1D integer array
[1313]	20	! Index1DArrayR: Function to provide the first index of a given value inside a 1D real array
[1608]	21	! Index1DArrayR_K: Function to provide the first index of a given value inside a 1D real(r_k) array
[2340]	22	! Index1DArrayL: Function to provide the first index of a given value inside a 1D boolean array
[1313]	23	! Index2DArrayR: Function to provide the first index of a given value inside a 2D real array
[1653]	24	! Index2DArrayR_K: Function to provide the first index of a given value inside a 2D real(r_k) array
[2724]	25	! index_samevals1D_RK: Subroutine to search for the indices of the same values between 2 1D RK series
	26	! of values allowing repetitions
	27	! index_samevals2D_RK: Subroutine to search for the indices of the same values between 2 2D RK series
	28	! of values allowing repetitions
[2475]	29	! JD: Fucntion to compute the julian date (JD) given a gregorian calendar
[1660]	30	! Nvalues_2DArrayI: Number of different values of a 2D integer array
[1619]	31	! mat2DPosition: Function to provide the i, j indices of a given value inside a 2D matrix
[2340]	32	! multi_Index1DArrayL: Subroutine to provide the indices of a given value inside a 1D boolean array
[2235]	33	! numberTimes: Function to provide the number of times that a given set of characters happen within a string
[1608]	34	! RangeI: Function to provide a range of d1 values from 'iniv' to 'endv', of integer values in a vector
	35	! RangeR: Function to provide a range of d1 values from 'iniv' to 'endv', of real values in a vector
	36	! RangeR_K: Function to provide a range of d1 from 'iniv' to 'endv', of real(r_k) values in a vector
[2618]	37	! rm_values_vecRK: Subroutine to remove a given value (e.g. fill_Value) from a RK vector
[1655]	38	! stoprun: Subroutine to stop running and print a message
[2331]	39	! vectorI_S: Function to transform a vector of integers to a string of characters
	40	! vectorR_S: Function to transform a vector of reals to a string of characters
[2330]	41	! xzones_homogenization: Subroutine to homogenize 2D contiguous zones along x-axis. Zones might be
	42	! contiguous, but with different number assigned !
	43	! yzones_homogenization: Subroutine to homogenize 2D contiguous zones along y-axis. Zones might be
	44	! contiguous, but with different number assigned !
[715]	45
[1608]	46	USE module_definitions
[1655]	47	USE module_basic
[1184]	48
[715]	49	CONTAINS
	50
[1660]	51	SUBROUTINE Nvalues_2DArrayI(dx, dy, dxy, mat2DI, Nvals, vals)
	52	! Subroutine to give the number of different values of a 2D integer array
	53
	54	IMPLICIT NONE
	55
	56	INTEGER, INTENT(in) :: dx, dy, dxy
	57	INTEGER, DIMENSION(dx,dy), INTENT(in) :: mat2DI
	58	INTEGER, INTENT(out) :: Nvals
	59	INTEGER, DIMENSION(dxy), INTENT(out) :: vals
	60
	61	! Local
	62	INTEGER :: i, j, ij
	63
	64	!!!!!!! Variables
	65	! dx, dy: size of the 2D space
	66	! mat2DI: 2D integer matrix
	67	! Nvals: number of different values
	68	! vals: vector with the different values
	69
	70	fname = 'Nvalues_2DArrayI'
	71
	72	vals = 0
	73
	74	Nvals = 1
	75	vals(1) = mat2DI(1,1)
	76	DO i=1,dx
	77	DO j=1,dy
	78	IF (Index1DArrayI(vals, Nvals, mat2DI(i,j)) == -1) THEN
	79	Nvals = Nvals + 1
	80	vals(Nvals) = mat2DI(i,j)
	81	END IF
	82	END DO
	83	END DO
	84
	85	RETURN
	86
	87	END SUBROUTINE Nvalues_2DArrayI
	88
[1612]	89	INTEGER FUNCTION index_list_coordsI(Ncoords, coords, icoord)
	90	! Function to provide the index of a given coordinate within a list of integer coordinates
	91
	92	IMPLICIT NONE
	93
	94	INTEGER, INTENT(in) :: Ncoords
	95	INTEGER, DIMENSION(Ncoords,2), INTENT(in) :: coords
	96	INTEGER, DIMENSION(2), INTENT(in) :: icoord
	97
	98	! Local
	99	INTEGER, DIMENSION(Ncoords) :: dist
	100	INTEGER :: i,mindist
	101	INTEGER, DIMENSION(1) :: iloc
	102
	103	!!!!!!! Variables
	104	! Ncoords: number of coordinates in the list
	105	! coords: list of coordinates
	106	! icoord: coordinate to find
	107
	108	fname = 'index_list_coordsI'
	109
	110	dist = (coords(:,1)-icoord(1))2+(coords(:,2)-icoord(2))2
	111
	112	IF (ANY(dist == 0)) THEN
	113	iloc = MINLOC(dist)
	114	index_list_coordsI = iloc(1)
	115	ELSE
	116	index_list_coordsI = -1
	117	END IF
	118
	119	END FUNCTION index_list_coordsI
	120
[1615]	121	INTEGER FUNCTION Index1DArrayI(array1D, d1, val)
	122	! Function to provide the first index of a given value inside a 1D integer array
	123
	124	IMPLICIT NONE
	125
	126	INTEGER, INTENT(in) :: d1
	127	INTEGER, INTENT(in) :: val
	128	INTEGER, DIMENSION(d1), INTENT(in) :: array1D
	129
	130	! Local
	131	INTEGER :: i
	132
	133	fname = 'Index1DArrayI'
	134
	135	Index1DArrayI = -1
	136
	137	DO i=1,d1
	138	IF (array1d(i) == val) THEN
	139	Index1DArrayI = i
	140	EXIT
	141	END IF
	142	END DO
	143
	144	END FUNCTION Index1DArrayI
	145
[764]	146	INTEGER FUNCTION Index1DArrayR(array1D, d1, val)
	147	! Function to provide the first index of a given value inside a 1D real array
	148
	149	IMPLICIT NONE
	150
	151	INTEGER, INTENT(in) :: d1
[1653]	152	REAL, INTENT(in) :: val
	153	REAL, DIMENSION(d1), INTENT(in) :: array1D
[764]	154
	155	! Local
	156	INTEGER :: i
	157
	158	fname = 'Index1DArrayR'
	159
	160	Index1DArrayR = -1
	161
	162	DO i=1,d1
	163	IF (array1d(i) == val) THEN
	164	Index1DArrayR = i
	165	EXIT
	166	END IF
	167	END DO
	168
	169	END FUNCTION Index1DArrayR
	170
[1608]	171	INTEGER FUNCTION Index1DArrayR_K(array1D, d1, val)
	172	! Function to provide the first index of a given value inside a 1D real(r_k) array
	173
	174	IMPLICIT NONE
	175
	176	INTEGER, INTENT(in) :: d1
	177	REAL(r_k), INTENT(in) :: val
	178	REAL(r_k), DIMENSION(d1), INTENT(in) :: array1D
	179
	180	! Local
	181	INTEGER :: i
	182
	183	fname = 'Index1DArrayR_K'
	184
	185	Index1DArrayR_K = -1
	186
	187	DO i=1,d1
	188	IF (array1d(i) == val) THEN
	189	Index1DArrayR_K = i
	190	EXIT
	191	END IF
	192	END DO
	193
	194	END FUNCTION Index1DArrayR_K
	195
[2340]	196	INTEGER FUNCTION Index1DArrayL(array1D, d1, val)
	197	! Function to provide the first index of a given value inside a 1D boolean array
	198
	199	IMPLICIT NONE
	200
	201	INTEGER, INTENT(in) :: d1
	202	LOGICAL, INTENT(in) :: val
	203	LOGICAL, DIMENSION(d1), INTENT(in) :: array1D
	204
	205	! Local
	206	INTEGER :: i
	207	CHARACTER(LEN=50) :: fname
	208
	209	fname = 'Index1DArrayL'
	210
	211	Index1DArrayL = -1
	212
	213	DO i=1,d1
	214	IF (array1d(i) .EQV. val) THEN
	215	Index1DArrayL = i
	216	EXIT
	217	END IF
	218	END DO
	219
	220	END FUNCTION Index1DArrayL
	221
	222	SUBROUTINE multi_Index1DArrayL(array1D, d1, val, Ntimes, pos)
	223	! Subroutine to provide the indices of a given value inside a 1D boolean array
	224
	225	IMPLICIT NONE
	226
	227	INTEGER, INTENT(in) :: d1
	228	LOGICAL, INTENT(in) :: val
	229	LOGICAL, DIMENSION(d1), INTENT(in) :: array1D
	230	INTEGER, INTENT(out) :: Ntimes
	231	INTEGER, DIMENSION(d1), INTENT(out) :: pos
	232
	233	! Local
	234	INTEGER :: i
	235	CHARACTER(LEN=50) :: fname
	236
	237	!!!!!!! Variables
	238	! array1D: 1D Array of values
	239	! d1: length of array
	240	! val: Value to look for
	241	! Ntimes: Number of times val is found within array1D
	242	! pos: positions of the values of val
	243
	244	fname = 'multi_Index1DArrayL'
	245
	246	Ntimes = 0
	247	pos = -1
	248
	249	DO i=1,d1
	250	IF (array1d(i) .EQV. val) THEN
	251	Ntimes = Ntimes + 1
	252	pos(Ntimes) = i
	253	END IF
	254	END DO
	255
	256	END SUBROUTINE multi_Index1DArrayL
	257
[715]	258	FUNCTION Index2DArrayR(array2D, d1, d2, val)
	259	! Function to provide the first index of a given value inside a 2D real array
	260
	261	IMPLICIT NONE
	262
	263	INTEGER, INTENT(in) :: d1, d2
[1653]	264	REAL, INTENT(in) :: val
	265	REAL, DIMENSION(d1,d2), INTENT(in) :: array2D
[715]	266	INTEGER, DIMENSION(2) :: Index2DArrayR
	267
	268	! Local
	269	INTEGER :: i, j
	270
	271	fname = 'Index2DArrayR'
	272
	273	Index2DArrayR = -1
	274
	275	DO i=1,d1
	276	DO j=1,d2
	277	IF (array2d(i,j) == val) THEN
	278	Index2DArrayR(1) = i
	279	Index2DArrayR(2) = j
	280	EXIT
	281	END IF
	282	END DO
	283	END DO
	284
	285	END FUNCTION Index2DArrayR
	286
[1653]	287	FUNCTION Index2DArrayR_K(array2D, d1, d2, val)
	288	! Function to provide the first index of a given value inside a 2D real array
	289
	290	IMPLICIT NONE
	291
	292	INTEGER, INTENT(in) :: d1, d2
	293	REAL(r_k), INTENT(in) :: val
	294	REAL(r_k), DIMENSION(d1,d2), INTENT(in) :: array2D
	295	INTEGER, DIMENSION(2) :: Index2DArrayR_K
	296
	297	! Local
	298	INTEGER :: i, j
	299
	300	fname = 'Index2DArrayR_K'
	301
	302	Index2DArrayR_K = -1
	303
	304	DO i=1,d1
	305	DO j=1,d2
	306	IF (array2d(i,j) == val) THEN
	307	Index2DArrayR_K(1) = i
	308	Index2DArrayR_K(2) = j
	309	EXIT
	310	END IF
	311	END DO
	312	END DO
	313
	314	END FUNCTION Index2DArrayR_K
	315
[2235]	316	INTEGER FUNCTION numberTimes(String, charv)
	317	! Function to provide the number of times that a given set of characters happen within a string
	318
	319	IMPLICIT NONE
	320
	321	CHARACTER(LEN=*), INTENT(IN) :: String, charv
	322
	323	! Local
	324	INTEGER :: i, Lstring, Lcharv
	325
	326	numberTimes = 0
	327
	328	Lstring = LEN_TRIM(String)
	329	Lcharv = LEN_TRIM(charv)
	330
	331	DO i=1,Lstring - Lcharv
	332	IF (String(i:i+Lcharv-1) == TRIM(charv)) numberTimes = numberTimes + 1
	333	END DO
	334
	335	END FUNCTION numberTimes
	336
	337
[1608]	338	FUNCTION RangeI(d1, iniv, endv)
	339	! Function to provide a range of d1 values from 'iniv' to 'endv', of integer values in a vector
	340
	341	IMPLICIT NONE
	342
	343	INTEGER, INTENT(in) :: d1, iniv, endv
	344	INTEGER, DIMENSION(d1) :: RangeI
	345
	346	! Local
	347	INTEGER :: i, intv
	348
	349	fname = 'RangeI'
	350
	351	intv = (endv - iniv) / (d1*1 - 1)
	352
	353	RangeI(1) = iniv
	354	DO i=2,d1
	355	RangeI(i) = RangeI(i-1) + intv
	356	END DO
	357
	358	END FUNCTION RangeI
	359
	360	FUNCTION RangeR(d1, iniv, endv)
	361	! Function to provide a range of d1 from 'iniv' to 'endv', of real values in a vector
	362
	363	IMPLICIT NONE
	364
	365	INTEGER, INTENT(in) :: d1
	366	REAL, INTENT(in) :: iniv, endv
	367	REAL, DIMENSION(d1) :: RangeR
	368
	369	! Local
	370	INTEGER :: i
	371	REAL :: intv
	372
	373	fname = 'RangeR'
	374
	375	intv = (endv - iniv) / (d1*1. - 1.)
	376
	377	RangeR(1) = iniv
	378	DO i=2,d1
	379	RangeR(i) = RangeR(i-1) + intv
	380	END DO
	381
	382	END FUNCTION RangeR
	383
	384	FUNCTION RangeR_K(d1, iniv, endv)
	385	! Function to provide a range of d1 from 'iniv' to 'endv', of real(r_k) values in a vector
	386
	387	IMPLICIT NONE
	388
	389	INTEGER, INTENT(in) :: d1
	390	REAL(r_k), INTENT(in) :: iniv, endv
	391	REAL(r_k), DIMENSION(d1) :: RangeR_K
	392
	393	! Local
	394	INTEGER :: i
	395	REAL(r_k) :: intv
	396
	397	fname = 'RangeR_K'
	398
	399	intv = (endv - iniv) / (d1*oneRK-oneRK)
	400
	401	RangeR_K(1) = iniv
	402	DO i=2,d1
	403	RangeR_K(i) = RangeR_K(i-1) + intv
	404	END DO
	405
	406	END FUNCTION RangeR_K
	407
[1655]	408	INTEGER FUNCTION freeunit()
	409	! provides the number of a free unit in which open a file
[1608]	410
[1619]	411	IMPLICIT NONE
	412
[1655]	413	LOGICAL :: is_used
[1619]	414
[1655]	415	is_used = .true.
	416	DO freeunit=10,100
[1660]	417	INQUIRE(unit=freeunit, opened=is_used)
[1655]	418	IF (.not. is_used) EXIT
	419	END DO
	420
	421	RETURN
	422
	423	END FUNCTION freeunit
	424
	425	SUBROUTINE GetInNamelist(namelistfile, param, kindparam, Ival, Rval, Lval, Sval)
	426	! Subroutine to get a paramter from a namelistfile
	427
	428	IMPLICIT NONE
	429
	430	CHARACTER(LEN=*), INTENT(IN) :: namelistfile, param
	431	CHARACTER(LEN=1), INTENT(IN) :: kindparam
	432	INTEGER, OPTIONAL, INTENT(OUT) :: Ival
	433	REAL, OPTIONAL, INTENT(OUT) :: Rval
	434	LOGICAL, OPTIONAL, INTENT(OUT) :: Lval
	435	CHARACTER(LEN=200), OPTIONAL, INTENT(OUT) :: Sval
	436
	437	! Local
	438	INTEGER :: i, funit, ios
	439	INTEGER :: Lparam, posparam
	440	LOGICAL :: is_used
	441	CHARACTER(LEN=1000) :: line, message
	442	CHARACTER(LEN=200), DIMENSION(2) :: lvals
	443	CHARACTER(LEN=200) :: pval
	444
	445	!!!!!!! Variables
	446	! namelistfile: name of the namelist file
	447	! param: parameter to get
	448	! paramkind: kind of the parameter (I: Integer, L: boolean, R: Real, S: String)
	449
	450	fname = 'GetInNamelist'
	451
	452	! Reading dimensions file and defining dimensions
	453	is_used = .true.
	454	DO funit=10,100
	455	INQUIRE(unit=funit, opened=is_used)
	456	IF (.not. is_used) EXIT
	457	END DO
	458
	459	OPEN(funit, FILE=TRIM(namelistfile), STATUS='old', FORM='formatted', IOSTAT=ios)
	460	IF ( ios /= 0 ) CALL stoprun(message, fname)
	461
	462	Lparam = LEN_TRIM(param)
	463
	464	DO i=1,10000
	465	READ(funit,"(A200)",END=100)line
	466	posparam = INDEX(TRIM(line), TRIM(param))
	467	IF (posparam /= 0) EXIT
	468
	469	END DO
	470	100 CONTINUE
	471
	472	IF (posparam == 0) THEN
	473	message = "namelist '" // TRIM(namelistfile) // "' does not have parameter '" // TRIM(param) // &
	474	"' !!"
	475	CALL stoprun(message, fname)
	476	END IF
	477
	478	CLOSE(UNIT=funit)
	479
	480	CALL split(line, '=', 2, lvals)
	481	IF (kindparam /= 'S') THEN
	482	CALL RemoveNonNum(lvals(2), pval)
	483	END IF
	484
	485	! L. Fita, LMD. October 2015
	486	! Up to now, only getting scalar values
	487	kparam: SELECT CASE (kindparam)
	488	CASE ('I')
	489	Ival = StoI(pval)
	490	! PRINT *,TRIM(param),'= ', Ival
	491	CASE ('L')
	492	Lval = StoL(pval)
	493	! PRINT *,TRIM(param),'= ', Lval
	494	CASE ('R')
	495	Rval = StoR(pval)
	496	! PRINT *,TRIM(param),'= ', Rval
	497	CASE ('S')
	498	Sval = lvals(2)
	499
	500	CASE DEFAULT
	501	message = "type of parameter '" // kindparam // "' not ready !!"
	502	CALL stoprun(message, fname)
	503
	504	END SELECT kparam
	505
	506	END SUBROUTINE GetInNamelist
	507
	508	SUBROUTINE stoprun(msg, fname)
	509	! Subroutine to stop running and print a message
	510
	511	IMPLICIT NONE
	512
	513	CHARACTER(LEN=*), INTENT(IN) :: fname
	514	CHARACTER(LEN=*), INTENT(IN) :: msg
	515
	516	! local
	517	CHARACTER(LEN=50) :: errmsg, warnmsg
	518
	519	errmsg = 'ERROR -- error -- ERROR -- error'
	520
	521	PRINT *, TRIM(errmsg)
	522	PRINT *, ' ' // TRIM(fname) // ': ' // TRIM(msg)
[1683]	523	STOP
[1655]	524
	525	END SUBROUTINE stoprun
	526
[2330]	527	SUBROUTINE xzones_homogenization(dx, dy, inzones, outzones)
	528	! Subroutine to homogenize 2D contiguous zones along x-axis, zones might be contiguous, but with
	529	! different number assigned !
[2328]	530	! Here we have a 2D matrix of integers, with contiguous integer filled zones, zero outside any zone
	531	! It might be, that within the same zone, might be lines which do not share the same integer
	532	! 0 0 0 0 0 0 0 0 0 0
	533	! 0 1 1 0 0 0 1 1 0 0
	534	! 0 2 0 0 1 == > 0 1 0 0 2
	535	! 0 1 1 0 0 0 1 1 0 0
	536
	537	IMPLICIT NONE
	538
	539	INTEGER, INTENT(in) :: dx, dy
	540	INTEGER, DIMENSION(dx,dy), INTENT(in) :: inzones
	541	INTEGER, DIMENSION(dx,dy), INTENT(out) :: outzones
	542
	543	! Local
	544	INTEGER :: i,j,k
	545	INTEGER :: Nmaxzones, TOTzones
	546	LOGICAL :: assigned
[2330]	547	INTEGER, DIMENSION(dx) :: prevline
	548	INTEGER, DIMENSION(dy) :: Nxzones
	549	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: zones
	550
	551	!!!!!!! Variables
	552	! dx, dy: Shape of the 2D space
	553	! inzones: zones to homogenize
	554	! outzones: zones homogenized
	555
	556	fname = 'xzones_homogenization'
	557
	558	! Maximum possible number of zones
	559	Nmaxzones = INT((dx/2)*(dy/2))
	560
	561	! Matrix with [i,j,Nzone,izone/ezone]
	562	IF (ALLOCATED(zones)) DEALLOCATE(zones)
	563	ALLOCATE(zones(dy,Nmaxzones,3))
	564
	565	zones = 0
	566	Nxzones = 0
	567	! Getting beginning/end of y-bands
	568	DO j=1, dy
	569	k = 0
	570	i = 1
	571	IF (inzones(i,j) /= 0) THEN
	572	k = k + 1
	573	zones(j,k,1) = i
	574	zones(j,k,3) = k
	575	END IF
	576	DO i=2, dx
	577	IF ( (inzones(i,j) /= 0) .AND. (inzones(i-1,j) == 0)) THEN
	578	k = k+1
	579	zones(j,k,1) = i
	580	zones(j,k,3) = k
	581	ELSE IF ( (inzones(i-1,j) /= 0) .AND. (inzones(i,j) == 0)) THEN
	582	zones(j,k,2) = i-1
	583	zones(j,k,3) = k
	584	END IF
	585	END DO
	586	IF (k > 0) THEN
	587	IF (zones(j,k,2) == 0) zones(j,k,2) = dx
	588	END IF
	589	Nxzones(j) = k
	590	END DO
	591
	592	! Homogenizing contigous zones
	593	outzones = 0
	594	TOTzones = 0
	595	j = 1
	596	DO k = 1, Nxzones(j)
	597	TOTzones = TOTzones + 1
	598	DO i=zones(j,k,1), zones(j,k,2)
	599	outzones(i,j) = TOTzones
	600	END DO
	601	END DO
	602
	603	DO j=2, dy
	604	prevline = outzones(:,j-1)
	605	DO k = 1, Nxzones(j)
	606	assigned = .FALSE.
	607	DO i=zones(j,k,1), zones(j,k,2)
	608	IF (prevline(i) /= 0) THEN
	609	outzones(zones(j,k,1):zones(j,k,2),j) = prevline(i)
	610	assigned = .TRUE.
	611	EXIT
	612	END IF
	613	END DO
	614	IF (.NOT.assigned) THEN
	615	TOTzones = TOTzones + 1
	616	DO i=zones(j,k,1), zones(j,k,2)
	617	outzones(i,j) = TOTzones
	618	END DO
	619	END IF
	620	END DO
	621	END DO
	622
	623	IF (ALLOCATED(zones)) DEALLOCATE(zones)
	624
	625	END SUBROUTINE xzones_homogenization
	626
	627	SUBROUTINE yzones_homogenization(dx, dy, inzones, outzones)
	628	! Subroutine to homogenize 2D contiguous zones along y-axis, zones might be contiguous, but with
	629	! different number assigned !
	630	! Here we have a 2D matrix of integers, with contiguous integer filled zones, zero outside any zone
	631	! It might be, that within the same zone, might be lines which do not share the same integer
	632	! 0 0 0 0 0 0 0 0 0 0
	633	! 0 1 1 0 0 0 1 1 0 0
	634	! 0 2 0 0 1 == > 0 1 0 0 2
	635	! 0 1 1 0 0 0 1 1 0 0
	636
	637	IMPLICIT NONE
	638
	639	INTEGER, INTENT(in) :: dx, dy
	640	INTEGER, DIMENSION(dx,dy), INTENT(in) :: inzones
	641	INTEGER, DIMENSION(dx,dy), INTENT(out) :: outzones
	642
	643	! Local
	644	INTEGER :: i,j,k
	645	INTEGER :: Nmaxzones, TOTzones
	646	LOGICAL :: assigned
[2328]	647	INTEGER, DIMENSION(dy) :: prevline
	648	INTEGER, DIMENSION(dx) :: Nyzones
	649	INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: zones
	650
	651	!!!!!!! Variables
	652	! dx, dy: Shape of the 2D space
	653	! inzones: zones to homogenize
	654	! outzones: zones homogenized
	655
[2330]	656	fname = 'yzones_homogenization'
[2328]	657
	658	! Maximum possible number of zones
	659	Nmaxzones = INT((dx/2)*(dy/2))
	660
	661	! Matrix with [i,j,Nzone,izone/ezone]
	662	IF (ALLOCATED(zones)) DEALLOCATE(zones)
	663	ALLOCATE(zones(dx,Nmaxzones,3))
	664
	665	zones = 0
	666	Nyzones = 0
	667	! Getting beginning/end of y-bands
	668	DO i=1, dx
	669	k = 0
	670	j = 1
	671	IF (inzones(i,j) /= 0) THEN
	672	k = k + 1
	673	zones(i,k,1) = j
	674	zones(i,k,3) = k
	675	END IF
	676	DO j=2, dy
	677	IF ( (inzones(i,j) /= 0) .AND. (inzones(i,j-1) == 0)) THEN
	678	k = k+1
	679	zones(i,k,1) = j
	680	zones(i,k,3) = k
	681	ELSE IF ( (inzones(i,j-1) /= 0) .AND. (inzones(i,j) == 0)) THEN
	682	zones(i,k,2) = j-1
	683	zones(i,k,3) = k
	684	END IF
	685	END DO
	686	IF (k > 0) THEN
	687	IF (zones(i,k,2) == 0) zones(i,k,2) = dy
	688	END IF
	689	Nyzones(i) = k
	690	END DO
	691
	692	! Homogenizing contigous zones
	693	outzones = 0
	694	TOTzones = 0
	695	i = 1
	696	DO k = 1, Nyzones(i)
	697	TOTzones = TOTzones + 1
	698	DO j=zones(i,k,1), zones(i,k,2)
	699	outzones(i,j) = TOTzones
	700	END DO
	701	END DO
	702
	703	DO i=2, dx
	704	prevline = outzones(i-1,:)
	705	DO k = 1, Nyzones(i)
	706	assigned = .FALSE.
	707	DO j=zones(i,k,1), zones(i,k,2)
	708	IF (prevline(j) /= 0) THEN
	709	outzones(i,zones(i,k,1):zones(i,k,2)) = prevline(j)
	710	assigned = .TRUE.
	711	EXIT
	712	END IF
	713	END DO
	714	IF (.NOT.assigned) THEN
	715	TOTzones = TOTzones + 1
	716	DO j=zones(i,k,1), zones(i,k,2)
	717	outzones(i,j) = TOTzones
	718	END DO
	719	END IF
	720	END DO
	721	END DO
	722
	723	IF (ALLOCATED(zones)) DEALLOCATE(zones)
	724
[2330]	725	END SUBROUTINE yzones_homogenization
[2328]	726
[2331]	727	CHARACTER(len=1000) FUNCTION vectorI_S(d1, vector)
	728	! Function to transform a vector of integers to a string of characters
	729
	730	IMPLICIT NONE
	731
	732	INTEGER, INTENT(in) :: d1
	733	INTEGER, DIMENSION(d1), INTENT(in) :: vector
	734
	735	! Local
	736	INTEGER :: iv
	737	CHARACTER(len=50) :: IS
	738
	739	!!!!!!! Variables
	740	! d1: length of the vector
	741	! vector: values to transform
	742
	743	fname = 'vectorI_S'
	744
	745	vectorI_S = ''
	746	DO iv=1, d1
	747	WRITE(IS, '(I50)')vector(iv)
	748	IF (iv == 1) THEN
	749	vectorI_S = TRIM(IS)
	750	ELSE
	751	vectorI_S = TRIM(vectorI_S) // ', ' // TRIM(IS)
	752	END IF
	753	END DO
	754
	755	END FUNCTION vectorI_S
	756
	757	CHARACTER(len=1000) FUNCTION vectorR_S(d1, vector)
	758	! Function to transform a vector of reals to a string of characters
	759
	760	IMPLICIT NONE
	761
	762	INTEGER, INTENT(in) :: d1
	763	REAL, DIMENSION(d1), INTENT(in) :: vector
	764
	765	! Local
	766	INTEGER :: iv
	767	CHARACTER(len=50) :: RS
	768
	769	!!!!!!! Variables
	770	! d1: length of the vector
	771	! vector: values to transform
	772
	773	fname = 'vectorR_S'
	774
	775	vectorR_S = ''
	776	DO iv=1, d1
	777	WRITE(RS, '(F50.25)')vector(iv)
	778	IF (iv == 1) THEN
	779	vectorR_S = TRIM(RS)
	780	ELSE
	781	vectorR_S = TRIM(vectorR_S) // ', ' // TRIM(RS)
	782	END IF
	783	END DO
	784
	785	END FUNCTION vectorR_S
	786
[2340]	787	SUBROUTINE continguos_homogene_zones(dx, dy, matvals, Nzones, contzones)
	788	! Subroutine to look for contiguous zones by looking by continuous grid points
	789
	790	IMPLICIT NONE
	791
	792	INTEGER, INTENT(in) :: dx, dy
	793	INTEGER, DIMENSION(dx,dy), INTENT(in) :: matvals
	794	INTEGER, INTENT(out) :: Nzones
	795	INTEGER, DIMENSION(dx,dy), INTENT(out) :: contzones
	796	! Local
	797	INTEGER :: i,j, k
	798	INTEGER :: ii, ei, ij, ej
	799	INTEGER :: Ncont, Nassigned, Ncontmin
	800	LOGICAL, DIMENSION(dx,dy) :: assigned, notdone
	801	INTEGER, DIMENSION(:), ALLOCATABLE :: pzones, allzones
	802	LOGICAL, DIMENSION(:), ALLOCATABLE :: passigns
	803
	804	!!!!!!! Variables
	805	! dx, dy: shape of the matrix
	806	! matvals: matrix with the values
	807	! contzones: homogeneous zones found
	808
	809	fname = 'continguos_homogene_zones'
	810
	811	! Vector to keep track of all zone values
	812	IF (ALLOCATED(allzones)) DEALLOCATE(allzones)
	813	ALLOCATE(allzones(dx*dy/4))
	814	allzones = 0
	815
	816	assigned = .FALSE.
	817	notdone = .TRUE.
	818	contzones = -1
	819	Nzones = 0
	820	DO i=1, dx
	821	DO j=1, dy
	822	! First
	823	IF (matvals(i,j) /= 0 .AND. notdone(i,j)) THEN
	824	CALL get_xyconlimits(dx, dy, matvals, 0, i, j, ii, ei, ij, ej)
	825
	826	! Has any point of the rays already been assigned?
	827	! Along x
	828	IF (ALLOCATED(pzones)) DEALLOCATE(pzones)
	829	ALLOCATE(pzones(ei-ii+1))
	830	IF (ALLOCATED(passigns)) DEALLOCATE(passigns)
	831	ALLOCATE(passigns(ei-ii+1))
	832	passigns = assigned(ii:ei,j)
	833	CALL multi_Index1DArrayL(passigns, ei-ii+1, .TRUE., Nassigned, pzones)
	834	IF (Nassigned /= 0) THEN
	835	Ncontmin = 10000000
	836	DO k=1, Nassigned
	837	IF (contzones(ii+pzones(k)-1,j) < Ncontmin) Ncontmin = contzones(ii+pzones(k)-1,j)
	838	END DO
	839	! If there is more than one assigned value change all values to the minimum one
	840	DO k=1, Nassigned
	841	IF (contzones(ii+pzones(k)-1,j) /= Ncontmin) THEN
	842	WHERE (contzones == contzones(ii+pzones(k)-1,j) .AND. contzones /= Ncontmin)
	843	contzones = Ncontmin
	844	END WHERE
	845	allzones(contzones(ii+pzones(k)-1,j)) = 0
	846	END IF
	847	END DO
	848	Ncont = Ncontmin
	849	END IF
	850	! Along y
	851	IF (ALLOCATED(pzones)) DEALLOCATE(pzones)
	852	ALLOCATE(pzones(ej-ij+1))
	853	IF (ALLOCATED(passigns)) DEALLOCATE(passigns)
	854	ALLOCATE(passigns(ej-ij+1))
	855	passigns = assigned(i,ij:ej)
	856	CALL multi_Index1DArrayL(passigns, ej-ij+1, .TRUE., Nassigned, pzones)
	857	IF (Nassigned /= 0) THEN
	858	Ncontmin = 10000000
	859	DO k=1, Nassigned
	860	IF (contzones(i,ij+pzones(k)-1) < Ncontmin) Ncontmin = contzones(i,ij+pzones(k)-1)
	861	END DO
	862	! If there is more than one assigned value change all values to the minimum one
	863	DO k=1, Nassigned
	864	IF (contzones(i,ij+pzones(k)-1) /= Ncontmin) THEN
	865	WHERE (contzones == contzones(i,ij+pzones(k)-1) .AND. contzones /= Ncontmin)
	866	contzones = Ncontmin
	867	END WHERE
	868	allzones(contzones(i,ij+pzones(k)-1)) = 0
	869	END IF
	870	END DO
	871	Ncont = Ncontmin
	872	END IF
	873
	874	IF (.NOT.assigned(i,j)) THEN
	875	Nzones = Nzones + 1
	876	Ncont = Nzones
	877	allzones(Nzones) = 1
	878	ELSE
	879	Ncont = contzones(i,j)
	880	END IF
	881	contzones(i,j) = Ncont
	882	contzones(ii:ei,j) = Ncont
	883	contzones(i,ij:ej) = Ncont
	884	notdone(i,j) = .FALSE.
	885	assigned(i,j) = .TRUE.
	886	assigned(ii:ei,j) = .TRUE.
	887	assigned(i,ij:ej) = .TRUE.
	888	END IF
	889	END DO
	890	END DO
	891
	892	! Using allzones to provide continuous assigned values
	893	Nzones = 0
	894	DO k=1, dx*dy/4
	895	IF (allzones(k) /= 0) THEN
	896	Nzones = Nzones + 1
	897	IF (allzones(k) /= Nzones) THEN
	898	WHERE(contzones == allzones(k))
	899	contzones = Nzones
	900	END WHERE
	901	END IF
	902	END IF
	903	END DO
	904
	905	RETURN
	906
	907	END SUBROUTINE continguos_homogene_zones
	908
	909	SUBROUTINE get_xyconlimits(d1, d2, matv, NOval, i, j, ix, ex, iy, ey)
	910	! Subroutine for getting the limits of contiguous values from a given point in a 2D matrix
	911
	912	IMPLICIT NONE
	913
	914	INTEGER, INTENT(in) :: d1, d2, i, j, NOval
	915	INTEGER, DIMENSION(d1,d2), INTENT(in) :: matv
	916	INTEGER, INTENT(out) :: ix, ex, iy, ey
	917
	918	! Local
	919	INTEGER :: i1, j1
	920	LOGICAL :: found
	921
	922	!!!!!!! Variables
	923	! d1, d2: Shape of input 2D values
	924	! i, j: grid point from which get the contiguous values
	925	! NOval: value for grid points without data
	926	! matv: 2D matrx values
	927	! ix, ex, iy, ey: limits from a grid point
	928
	929	fname = 'get_xyconlimits'
	930
	931	ix = -1
	932	ex = -1
	933	iy = -1
	934	ey = -1
	935
	936	! Before i
	937	IF (i > 1) THEN
	938	ix = i
	939	found = .FALSE.
	940	DO i1=i-1,1,-1
	941	IF (matv(i1,j) == NOval) THEN
	942	ix = i1 + 1
	943	found = .TRUE.
	944	EXIT
	945	END IF
	946	END DO
	947	IF (.NOT.found) ix=1
	948	ELSE
	949	ix = 1
	950	END IF
	951	! After i
	952	IF (i < d1) THEN
	953	ex = i
	954	found = .FALSE.
	955	DO i1=i+1,d1
	956	IF (matv(i1,j) == NOval) THEN
	957	ex = i1 - 1
	958	found = .TRUE.
	959	EXIT
	960	END IF
	961	END DO
	962	IF (.NOT.found) ex=d1
	963	ELSE
	964	ex = d1
	965	END IF
	966
	967	! Before j
	968	IF (j > 1) THEN
	969	iy = j
	970	found = .FALSE.
	971	DO j1=j-1,1,-1
	972	IF (matv(i,j1) == NOval) THEN
	973	iy = j1 + 1
	974	found = .TRUE.
	975	EXIT
	976	END IF
	977	END DO
	978	IF (.NOT.found) iy=1
	979	ELSE
	980	iy = 1
	981	END IF
	982	! After j
	983	IF (j < d2) THEN
	984	ey = j
	985	found = .FALSE.
	986	DO j1=j+1,d2
	987	IF (matv(i,j1) == NOval) THEN
	988	ey = j1 - 1
	989	found = .TRUE.
	990	EXIT
	991	END IF
	992	END DO
	993	IF (.NOT.found) ey=d2
	994	ELSE
	995	ey = d2
	996	END IF
	997
	998	RETURN
	999
	1000	END SUBROUTINE get_xyconlimits
	1001
[2359]	1002	SUBROUTINE from_ptlist_2DRKmatrix(Npts, pts, Flike, vals, dx, dy, NOval, matrix)
[2357]	1003	! Subroutine to construct a 2D RK matrix from a list of values accompaigned by a list of grid-point
	1004	! coordinates
[2356]	1005
	1006	IMPLICIT NONE
	1007
	1008	INTEGER, INTENT(in) :: Npts, dx, dy
[2359]	1009	LOGICAL, INTENT(in) :: Flike
[2356]	1010	REAL(r_k), INTENT(in) :: NOval
	1011	INTEGER, DIMENSION(Npts,2), INTENT(in) :: pts
	1012	REAL(r_k), DIMENSION(Npts), INTENT(in) :: vals
	1013	REAL(r_k), DIMENSION(dx,dy), INTENT(out) :: matrix
	1014
	1015	! Local
	1016	INTEGER :: iv, i, j
	1017
	1018	!!!!!!! Variables
	1019	! Npts: Number of values of the list
	1020	! pts: 2D matrix coordinates of the values
[2359]	1021	! Flike: whether input is Fortran-like (1: x, 2: y) or not (0-based; 1: y, 2: x)
[2356]	1022	! vals: list of values correspondant to the coordinates
	1023	! dx, dy: shape of the 2D matrix
	1024	! NOval: Value to assign when there is no coordinate
	1025	! matrix: resultant matrix
	1026
	1027	fname = 'from_ptlist_2DRKmatrix'
	1028
	1029	matrix = NOval
	1030
[2359]	1031	IF (Flike) THEN
	1032	DO iv=1, Npts
	1033	i = pts(iv,1)
	1034	j = pts(iv,2)
	1035	matrix(i,j) = vals(iv)
	1036	END DO
	1037	ELSE
	1038	DO iv=1, Npts
	1039	i = pts(iv,2)+1
	1040	j = pts(iv,1)+1
	1041	matrix(i,j) = vals(iv)
	1042	END DO
	1043	END IF
[2356]	1044
	1045	RETURN
	1046
	1047	END SUBROUTINE from_ptlist_2DRKmatrix
	1048
[2359]	1049	SUBROUTINE from_ptlist_2DRKNmatrix(Npts, pts, Flike, Nvals, vals, dx, dy, NOval, Nmatrix)
[2357]	1050	! Subroutine to construct N 2D RK matrix from a list of values accompaigned by a list of grid-point
	1051	! coordinates
	1052
	1053	IMPLICIT NONE
	1054
	1055	INTEGER, INTENT(in) :: Npts, dx, dy, Nvals
[2359]	1056	LOGICAL, INTENT(in) :: Flike
[2357]	1057	REAL(r_k), INTENT(in) :: NOval
	1058	INTEGER, DIMENSION(Npts,2), INTENT(in) :: pts
	1059	REAL(r_k), DIMENSION(Npts,Nvals), INTENT(in) :: vals
	1060	REAL(r_k), DIMENSION(dx,dy,Nvals), INTENT(out) :: Nmatrix
	1061
	1062	! Local
	1063	INTEGER :: iv, i, j, ic
	1064
	1065	!!!!!!! Variables
	1066	! Npts: Number of values of the list
	1067	! pts: 2D matrix coordinates of the values
[2359]	1068	! Flike: whether input is Fortran-like (1: x, 2: y) or not (0-based; 1: y, 2: x)
[2357]	1069	! Nvals: number of different values
	1070	! vals: list of N-values correspondant to the coordinates
	1071	! dx, dy: shape of the 2D matrix
	1072	! NOval: Value to assign when there is no coordinate
	1073	! Nmatrix: resultant N-matrix
	1074
	1075	fname = 'from_ptlist_2DRKNmatrix'
	1076
	1077	Nmatrix = NOval
	1078
[2359]	1079	IF (Flike) THEN
	1080	DO iv=1, Npts
	1081	i = pts(iv,1)
	1082	j = pts(iv,2)
	1083	Nmatrix(i,j,:) = vals(iv,:)
	1084	END DO
	1085	ELSE
	1086	DO iv=1, Npts
	1087	i = pts(iv,2)+1
	1088	j = pts(iv,1)+1
	1089	Nmatrix(i,j,:) = vals(iv,:)
	1090	END DO
	1091	END IF
[2357]	1092
	1093	RETURN
	1094
	1095	END SUBROUTINE from_ptlist_2DRKNmatrix
	1096
[2359]	1097	SUBROUTINE from_coordlist_2DRKmatrix(Npts, coords, Flike, xcoord, ycoord, vals, dx, dy, NOval, &
	1098	matrix)
[2357]	1099	! Subroutine to construct a 2D RK matrix from a list of values accompaigned by a list of coordinates
	1100	! to find i,j grid-point coordinates by minimum distance
	1101
	1102	IMPLICIT NONE
	1103
	1104	INTEGER, INTENT(in) :: Npts, dx, dy
[2359]	1105	LOGICAL, INTENT(in) :: Flike
[2357]	1106	REAL(r_k), INTENT(in) :: NOval
	1107	REAL(r_k), DIMENSION(Npts,2), INTENT(in) :: coords
	1108	REAL(r_k), DIMENSION(Npts), INTENT(in) :: vals
	1109	REAL(r_k), DIMENSION(dx,dy), INTENT(in) :: xcoord, ycoord
	1110	REAL(r_k), DIMENSION(dx,dy), INTENT(out) :: matrix
	1111
	1112	! Local
	1113	INTEGER :: iv, i, j
	1114	REAL(r_k) :: xi, yi
	1115	REAL(r_k), DIMENSION(dx,dy) :: diff
	1116	INTEGER, DIMENSION(2) :: minpt
	1117
	1118	!!!!!!! Variables
	1119	! Npts: Number of values of the list
	1120	! coords: 2D coordinates of the values
[2359]	1121	! Flike: whether input is Fortran-like (1: x, 2: y) or not (0-based; 1: y, 2: x)
[2357]	1122	! vals: list of values correspondant to the coordinates
	1123	! xcoord, ycoord: matrix of values correspondant to the each coordinate
	1124	! dx, dy: shape of the 2D matrix
	1125	! NOval: Value to assign when there is no coordinate
	1126	! matrix: resultant matrix
	1127
	1128	fname = 'from_coordlist_2DRKmatrix'
	1129
	1130	matrix = NOval
	1131
[2359]	1132	IF (Flike) THEN
	1133	DO iv=1, Npts
	1134	xi = coords(iv,1)
	1135	yi = coords(iv,2)
	1136	diff = SQRT((xcoord-xi)2+(ycoord-yi)2)
	1137	minpt = MINLOC(diff)
	1138	matrix(minpt(1),minpt(2)) = vals(iv)
	1139	END DO
	1140	ELSE
	1141	DO iv=1, Npts
	1142	xi = coords(iv,2)
	1143	yi = coords(iv,1)
	1144	diff = SQRT((xcoord-xi)2+(ycoord-yi)2)
	1145	minpt = MINLOC(diff)
	1146	matrix(minpt(1),minpt(2)) = vals(iv)
	1147	END DO
	1148	END IF
[2357]	1149
	1150	RETURN
	1151
	1152	END SUBROUTINE from_coordlist_2DRKmatrix
	1153
[2475]	1154	INTEGER FUNCTION JD (year,month,day)
	1155	! Fucntion to compute the julian date (JD) given a gregorian calendar
	1156	! FROM: https://aa.usno.navy.mil/faq/docs/JD_Formula.php
	1157
	1158	IMPLICIT NONE
	1159
	1160	INTEGER, INTENT(in) :: year,month,day
	1161
	1162	! Local
	1163	INTEGER :: I,J,K
	1164
	1165	I= year
	1166	J= month
	1167	K= day
	1168
	1169	JD= K-32075+1461(I+4800+(J-14)/12)/4+367(J-2-(J-14)/1212)/12-3((I+4900+(J-14)/12)/100)/4
	1170
	1171	END FUNCTION JD
	1172
	1173	SUBROUTINE GDATE(jd,year,month,day)
	1174	! Subroutine to compute the gregorian calendar date (year,month,day) given the julian date (JD)
	1175
	1176	IMPLICIT NONE
	1177
	1178	INTEGER, INTENT(in) :: jd
	1179	INTEGER, INTENT(out) :: year,month,day
	1180	! Local
	1181	INTEGER :: I,J,K,L,N
	1182
	1183	L= JD+68569
	1184	N= 4*L/146097
	1185	L= L-(146097*N+3)/4
	1186	I= 4000*(L+1)/1461001
	1187	L= L-1461*I/4+31
	1188	J= 80*L/2447
	1189	K= L-2447*J/80
	1190	L= J/11
	1191	J= J+2-12*L
	1192	I= 100*(N-49)+I+L
	1193
	1194	year= I
	1195	month= J
	1196	day= K
	1197
	1198	RETURN
	1199
	1200	END SUBROUTINE GDATE
	1201
[2618]	1202	SUBROUTINE rm_values_vecRK(d1, vector, rmvalue, rmd1, rmvector)
	1203	! Subroutine to remove a given value (e.g. fill_Value) from a RK vector
	1204
	1205	IMPLICIT NONE
	1206
	1207	INTEGER, INTENT(in) :: d1
	1208	REAL(r_k), INTENT(in) :: rmvalue
	1209	REAL(r_k), DIMENSION(d1), INTENT(in) :: vector
	1210	INTEGER, INTENT(out) :: rmd1
	1211	REAL(r_k), DIMENSION(d1), INTENT(out) :: rmvector
	1212
	1213	! Local
	1214	INTEGER :: i
	1215
	1216	!!!!!!! Variables
	1217	! d1: length of the vector
	1218	! vector: vector with values
	1219	! rmvalue: RK value to remove
	1220	! rmd1: length of the new vector
	1221	! rmvector: new vector
	1222
	1223	fname = 'rm_values_vecRK'
	1224
	1225	rmd1 = 0
	1226	rmvector = zeroRK
	1227
	1228	DO i=1, d1
	1229	IF (vector(i) /= rmvalue) THEN
	1230	rmd1 = rmd1 + 1
	1231	rmvector(rmd1) = vector(i)
	1232	END IF
	1233	END DO
	1234
	1235	END SUBROUTINE rm_values_vecRK
	1236
[2725]	1237	SUBROUTINE index_samevals1D_RK(d1r, refv, d1v, vals, ii, indices, missval, samevalues)
[2724]	1238	! Subroutine to search for the indices of the same values between 2 1D RK series of values allowing
	1239	! repetitions
	1240
	1241	IMPLICIT NONE
	1242
	1243	INTEGER, INTENT(in) :: d1r, d1v
[2725]	1244	REAL(r_k), INTENT(in) :: missval
[2724]	1245	REAL(r_k), DIMENSION(d1r), INTENT(in) :: refv
	1246	REAL(r_k), DIMENSION(d1v), INTENT(in) :: vals
	1247	INTEGER, INTENT(out) :: ii
	1248	INTEGER, DIMENSION(d1v,2), INTENT(out) :: indices
	1249	REAL(r_k), DIMENSION(d1v), INTENT(out) :: samevalues
	1250
	1251	! Local
	1252	INTEGER :: ir, iv, iiv
	1253
	1254	!!!!!!! Variables
	1255	! d1r: size of the reference data
	1256	! refv: reference values
	1257	! d1v: size of the values
	1258	! vals: values to look in
	1259	! ii: quantity of same values found
	1260	! indices: output
[2725]	1261	! missval: missing value
[2724]	1262	! samevalues: values where coincidence is found
	1263
	1264	fname = 'index_samevals1D_RK'
	1265
	1266	indices = 0
	1267	samevalues = zeroRK
	1268	ii = 0
	1269	DO ir=1, d1r
[2726]	1270	IF (refv(ir) /= missval) THEN
	1271	DO iv=1, d1v
	1272	IF (vals(iv) == refv(ir)) THEN
	1273	ii = ii + 1
	1274	indices(ii,:) = (/ ir, iv /)
	1275	samevalues(ii) = refv(ir)
	1276	END IF
	1277	END DO
	1278	END IF
[2724]	1279	END DO
	1280
	1281	RETURN
	1282
	1283	END SUBROUTINE index_samevals1D_RK
	1284
[2726]	1285	SUBROUTINE index_samevals2D_RK(d1r, d2r, refv, d1v, d2v, d12v, vals, ii, indices, missval, &
	1286	samevalues)
[2724]	1287	! Subroutine to search for the indices of the same values between 2 2D RK series of values allowing
	1288	! repetitions
	1289
	1290	IMPLICIT NONE
	1291
	1292	INTEGER, INTENT(in) :: d1r, d2r, d1v, d2v, d12v
[2726]	1293	REAL(r_k), INTENT(in) :: missval
[2724]	1294	REAL(r_k), DIMENSION(d1r,d2r), INTENT(in) :: refv
	1295	REAL(r_k), DIMENSION(d1v,d2v), INTENT(in) :: vals
	1296	INTEGER, INTENT(out) :: ii
	1297	INTEGER, DIMENSION(d12v,2,2), INTENT(out) :: indices
	1298	REAL(r_k), DIMENSION(d12v), INTENT(out) :: samevalues
	1299
	1300	! Local
	1301	INTEGER :: ir1, ir2, iv1, iv2
	1302
	1303	!!!!!!! Variables
	1304	! d1r, d2r: size of the reference data
	1305	! refv: reference values
	1306	! d1v, d2v: size of the values
	1307	! vals: values to look in
	1308	! ii: quantity of same values found
	1309	! indices: output
[2726]	1310	! missval: missing value
[2724]	1311	! samevalues: values where coincidence is found
	1312
	1313	fname = 'index_samevals2D_RK'
	1314
	1315	indices = 0
	1316	samevalues = zeroRK
	1317	ii = 0
	1318	DO ir1=1, d1r
	1319	DO ir2=1, d2r
[2726]	1320	IF (refv(ir1,ir2) /= missval) THEN
	1321	DO iv1=1, d1v
	1322	DO iv2=1, d2v
	1323	IF (vals(iv1,iv2) == refv(ir1,ir2)) THEN
	1324	ii = ii + 1
	1325	indices(ii,1,1) = ir1
	1326	indices(ii,1,2) = ir2
	1327	indices(ii,2,1) = iv1
	1328	indices(ii,2,2) = iv2
	1329	samevalues(ii) = refv(ir1,ir2)
	1330	END IF
	1331	END DO
[2724]	1332	END DO
[2726]	1333	END IF
[2724]	1334	END DO
	1335	END DO
	1336
	1337	RETURN
	1338
	1339	END SUBROUTINE index_samevals2D_RK
	1340
	1341	SUBROUTINE fill_matrix2DRK_winmat2D_list1D(d1i, d2i, inmatrix, ind, dlist, inlist, olist, missval, &
	1342	od, d1o, d2o, omat)
	1343	! Subroutine to fill a 2D RK matrix using a list of 1D indices from another given 2D matrix
	1344
	1345	IMPLICIT NONE
	1346
	1347	INTEGER, INTENT(in) :: d1i, d2i, ind, dlist, od, d1o, d2o
	1348	REAL(r_k), DIMENSION(d1i, d2i), INTENT(in) :: inmatrix
	1349	INTEGER, DIMENSION(dlist), INTENT(in) :: inlist, olist
	1350	REAL(r_k), INTENT(in) :: missval
	1351	REAL(r_k), DIMENSION(d1o, d2o), INTENT(out) :: omat
	1352
	1353	! Local
	1354	INTEGER :: ii, ij, oi, oj, il
	1355	INTEGER :: isame, osame, irun, orun
[2727]	1356	INTEGER :: ilx, olx, iln, oln
[2724]	1357	CHARACTER(len=3) :: isS, osS
	1358
	1359	!!!!!!! Variables
	1360	! d1i, d2i: size of the input matrix
	1361	! inmatrix: input matrix with the values to fill the output matrix
	1362	! ind: dimension of the input matrix that the list of indices refer to
	1363	! dlist: number of indices from the list
	1364	! inlist: list of indices from the input matrix
	1365	! olist: list of indices from the output matrix
	1366	! missval: missing value
	1367	! od: dimension of the output matix to which assign the indices of the list
	1368	! d1o, d2o: size of the output matrix
	1369	! omat: output matrix
	1370
	1371	fname = 'fill_matrix2DRK_winmat2D_list1D'
	1372
	1373	omat = missval
	1374
	1375	IF (ind == 1) THEN
	1376	isame = d2i
	1377	irun = d1i
	1378	ELSE
	1379	isame = d1i
	1380	irun = d2i
	1381	END IF
	1382
	1383	IF (od == 1) THEN
	1384	osame = d2o
	1385	orun = d1o
	1386	ELSE
	1387	osame = d1o
	1388	orun = d2o
	1389	END IF
	1390
	1391	IF (isame /= osame) THEN
	1392	WRITE(isS,'(I3)')isame
	1393	WRITE(osS,'(I3)')osame
	1394	msg = 'Resultant working size from input ' // isS // ' to output ' // osS // ' differ !!'
	1395	CALL ErrMsg(msg, fname, -1)
	1396	END IF
	1397
[2727]	1398	iln = MINVAL(inlist)
	1399	oln = MINVAL(olist)
	1400
	1401	IF (iln < 1) THEN
	1402	WRITE(isS,'(I3)')iln
	1403	msg = 'Incorrect minimum input index: ' // isS // ' Negative or zero indices not allowed !!'
	1404	CALL ErrMsg(msg, fname, -1)
	1405	END IF
	1406
	1407	IF (oln < 1) THEN
	1408	WRITE(isS,'(I3)')oln
	1409	msg = 'Incorrect minimum output index: ' // isS // ' Negative or zero indices not allowed !!'
	1410	CALL ErrMsg(msg, fname, -1)
	1411	END IF
	1412
[2724]	1413	ilx = MAXVAL(inlist)
	1414	olx = MAXVAL(olist)
	1415
	1416	IF (ilx > irun) THEN
	1417	WRITE(isS,'(I3)')ilx
	1418	WRITE(osS,'(I3)')irun
	1419	msg = 'Maximum value in input indices ' // isS // ' larger than assigned dimension ' // osS // &
	1420	' !!'
	1421	CALL ErrMsg(msg, fname, -1)
	1422	END IF
	1423
	1424	IF (olx > orun) THEN
	1425	WRITE(isS,'(I3)')olx
	1426	WRITE(osS,'(I3)')orun
	1427	msg = 'Maximum value in output indices ' // isS // ' larger than assigned dimension ' // osS // &
	1428	' !!'
	1429	CALL ErrMsg(msg, fname, -1)
	1430	END IF
	1431
	1432	IF (od == 1) THEN
	1433	IF (ind == 1) THEN
	1434	DO il=1, dlist
	1435	omat(olist(il),:) = inmatrix(inlist(il),:)
	1436	END DO
	1437	ELSE
	1438	DO il=1, dlist
	1439	omat(olist(il),:) = inmatrix(:,inlist(il))
	1440	END DO
	1441	END IF
	1442	ELSE
	1443	IF (ind == 1) THEN
	1444	DO il=1, dlist
	1445	omat(:,olist(il)) = inmatrix(inlist(il),:)
	1446	END DO
	1447	ELSE
	1448	DO il=1, dlist
	1449	omat(:,olist(il)) = inmatrix(:,inlist(il))
	1450	END DO
	1451	END IF
	1452	END IF
	1453
	1454	RETURN
	1455
	1456	END SUBROUTINE fill_matrix2DRK_winmat2D_list1D
	1457
	1458
[715]	1459	END MODULE module_generic

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format