Context Navigation

module_generic.f90 @ 2462

Last change on this file since 2462 was 2359, checked in by lfita, 6 years ago
Adding `Flike' in different subroutines to take into account indices from non-Fortran origins
File size: 34.0 KB

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

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

Download in other formats:

Original Format