Context Navigation

module_generic.f90 @ 2640

Last change on this file since 2640 was 2618, checked in by lfita, 6 years ago

Adding:

`rm_values_vecRK': Subroutine to remove a given value (e.g. fill_Value) from a RK vector

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

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

Download in other formats:

Original Format