Context Navigation

module_generic.f90 @ 2828

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

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

Download in other formats:

Original Format