Context Navigation

module_generic.f90 @ 2495

Last change on this file since 2495 was 2476, checked in by lfita, 6 years ago

Adding:

`from360d_reg': Function to transform from 360d = 12 * 30d calendar to a regular one

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

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

Download in other formats:

Original Format