Context Navigation

module_generic.f90 @ 2515

Last change on this file since 2515 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

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	! stoprun: Subroutine to stop running and print a message
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
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 !
38
39	USE module_definitions
40	USE module_basic
41
42	CONTAINS
43
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
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
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
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
145	REAL, INTENT(in) :: val
146	REAL, DIMENSION(d1), INTENT(in) :: array1D
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
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
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
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
257	REAL, INTENT(in) :: val
258	REAL, DIMENSION(d1,d2), INTENT(in) :: array2D
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
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
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
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
401	INTEGER FUNCTION freeunit()
402	! provides the number of a free unit in which open a file
403
404	IMPLICIT NONE
405
406	LOGICAL :: is_used
407
408	is_used = .true.
409	DO freeunit=10,100
410	INQUIRE(unit=freeunit, opened=is_used)
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)
516	STOP
517
518	END SUBROUTINE stoprun
519
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 !
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
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
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
649	fname = 'yzones_homogenization'
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
718	END SUBROUTINE yzones_homogenization
719
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
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
995	SUBROUTINE from_ptlist_2DRKmatrix(Npts, pts, Flike, vals, dx, dy, NOval, matrix)
996	! Subroutine to construct a 2D RK matrix from a list of values accompaigned by a list of grid-point
997	! coordinates
998
999	IMPLICIT NONE
1000
1001	INTEGER, INTENT(in) :: Npts, dx, dy
1002	LOGICAL, INTENT(in) :: Flike
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
1014	! Flike: whether input is Fortran-like (1: x, 2: y) or not (0-based; 1: y, 2: x)
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
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
1037
1038	RETURN
1039
1040	END SUBROUTINE from_ptlist_2DRKmatrix
1041
1042	SUBROUTINE from_ptlist_2DRKNmatrix(Npts, pts, Flike, Nvals, vals, dx, dy, NOval, Nmatrix)
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
1049	LOGICAL, INTENT(in) :: Flike
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
1061	! Flike: whether input is Fortran-like (1: x, 2: y) or not (0-based; 1: y, 2: x)
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
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
1085
1086	RETURN
1087
1088	END SUBROUTINE from_ptlist_2DRKNmatrix
1089
1090	SUBROUTINE from_coordlist_2DRKmatrix(Npts, coords, Flike, xcoord, ycoord, vals, dx, dy, NOval, &
1091	matrix)
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
1098	LOGICAL, INTENT(in) :: Flike
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
1114	! Flike: whether input is Fortran-like (1: x, 2: y) or not (0-based; 1: y, 2: x)
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
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
1142
1143	RETURN
1144
1145	END SUBROUTINE from_coordlist_2DRKmatrix
1146
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
1195	END MODULE module_generic

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

Download in other formats:

Original Format