Context Navigation

lmdz_filtreg.F90 @ 5116

Last change on this file since 5116 was 5116, checked in by abarral, 2 months ago
rename modules properly lmdz_* move ismin, ismax, minmax into new lmdz_libmath.f90 (lint) uppercase fortran keywords
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 28.4 KB

Line
1	! $Id$
2
3	MODULE lmdz_filtreg
4	IMPLICIT NONE; PRIVATE
5	PUBLIC matriceun, matriceus, matricevn, matricevs, matrinvn, matrinvs, &
6	inifilr, filtreg
7
8	REAL, DIMENSION(:, :, :), ALLOCATABLE :: matriceun, matriceus, matricevn
9	REAL, DIMENSION(:, :, :), ALLOCATABLE :: matricevs, matrinvn, matrinvs
10
11	CONTAINS
12
13	SUBROUTINE filtreg(champ, nlat, nbniv, ifiltre, iaire, &
14	griscal, iter)
15	USE lmdz_coefils, ONLY: jfiltnu, jfiltnv, jfiltsu, jfiltsv, sddu, sddv, unsddu, unsddv, modfrstv, modfrstu
16
17	!=======================================================================
18	!
19	! Auteur: P. Le Van 07/10/97
20	! ------
21	!
22	! Objet: filtre matriciel longitudinal ,avec les matrices precalculees
23	! pour l'operateur Filtre .
24	! ------
25	!
26	! Arguments:
27	! ----------
28	!
29	! nblat nombre de latitudes a filtrer
30	! nbniv nombre de niveaux verticaux a filtrer
31	! champ(iip1,nblat,nbniv) en entree : champ a filtrer
32	! en sortie : champ filtre
33	! ifiltre +1 Transformee directe
34	! -1 Transformee inverse
35	! +2 Filtre directe
36	! -2 Filtre inverse
37	!
38	! iaire 1 si champ intensif
39	! 2 si champ extensif (pondere par les aires)
40	!
41	! iter 1 filtre simple
42	!
43	!=======================================================================
44	!
45	!
46	! Variable Intensive
47	! ifiltre = 1 filtre directe
48	! ifiltre =-1 filtre inverse
49	!
50	! Variable Extensive
51	! ifiltre = 2 filtre directe
52	! ifiltre =-2 filtre inverse
53	!
54	!
55	INCLUDE "dimensions.h"
56	INCLUDE "paramet.h"
57
58	INTEGER :: nlat, nbniv, ifiltre, iter
59	INTEGER :: i, j, l, k
60	INTEGER :: iim2, immjm
61	INTEGER :: jdfil1, jdfil2, jffil1, jffil2, jdfil, jffil
62
63	REAL :: champ(iip1, nlat, nbniv)
64
65	REAL :: eignq(iim, nlat, nbniv), sdd1(iim), sdd2(iim)
66	LOGICAL :: griscal
67	INTEGER :: hemisph, iaire
68
69	LOGICAL, SAVE :: first = .TRUE.
70
71	REAL, SAVE :: sdd12(iim, 4)
72
73	INTEGER, PARAMETER :: type_sddu = 1
74	INTEGER, PARAMETER :: type_sddv = 2
75	INTEGER, PARAMETER :: type_unsddu = 3
76	INTEGER, PARAMETER :: type_unsddv = 4
77
78	INTEGER :: sdd1_type, sdd2_type
79
80	if (iim == 1) return ! no filtre in 2D y-z
81
82	IF (first) THEN
83	sdd12(1:iim, type_sddu) = sddu(1:iim)
84	sdd12(1:iim, type_sddv) = sddv(1:iim)
85	sdd12(1:iim, type_unsddu) = unsddu(1:iim)
86	sdd12(1:iim, type_unsddv) = unsddv(1:iim)
87
88	first = .FALSE.
89	ENDIF
90
91	IF(ifiltre==1.or.ifiltre==-1) &
92	stop 'Pas de transformee simple dans cette version'
93
94	IF(iter== 2) THEN
95	PRINT *, ' Pas d iteration du filtre dans cette version !'&
96	&, ' Utiliser old_filtreg et repasser !'
97	STOP
98	ENDIF
99
100	IF(ifiltre== -2 .AND..NOT.griscal) THEN
101	PRINT *, ' Cette routine ne calcule le filtre inverse que ' &
102	, ' sur la grille des scalaires !'
103	STOP
104	ENDIF
105
106	IF(ifiltre/=2 .AND.ifiltre/= - 2) THEN
107	PRINT *, ' Probleme dans filtreg car ifiltre NE 2 et NE -2' &
108	, ' corriger et repasser !'
109	STOP
110	ENDIF
111
112	iim2 = iim * iim
113	immjm = iim * jjm
114
115	IF(griscal) THEN
116	IF(nlat /= jjp1) THEN
117	PRINT 1111
118	STOP
119	ELSE
120
121	IF(iaire==1) THEN
122	sdd1_type = type_sddv
123	sdd2_type = type_unsddv
124	ELSE
125	sdd1_type = type_unsddv
126	sdd2_type = type_sddv
127	ENDIF
128
129	! IF( iaire.EQ.1 ) THEN
130	! CALL SCOPY( iim, sddv, 1, sdd1, 1 )
131	! CALL SCOPY( iim, unsddv, 1, sdd2, 1 )
132	! ELSE
133	! CALL SCOPY( iim, unsddv, 1, sdd1, 1 )
134	! CALL SCOPY( iim, sddv, 1, sdd2, 1 )
135	! END IF
136
137	jdfil1 = 2
138	jffil1 = jfiltnu
139	jdfil2 = jfiltsu
140	jffil2 = jjm
141	END IF
142	ELSE
143	IF(nlat/=jjm) THEN
144	PRINT 2222
145	STOP
146	ELSE
147
148	IF(iaire==1) THEN
149	sdd1_type = type_sddu
150	sdd2_type = type_unsddu
151	ELSE
152	sdd1_type = type_unsddu
153	sdd2_type = type_sddu
154	ENDIF
155
156	! IF( iaire.EQ.1 ) THEN
157	! CALL SCOPY( iim, sddu, 1, sdd1, 1 )
158	! CALL SCOPY( iim, unsddu, 1, sdd2, 1 )
159	! ELSE
160	! CALL SCOPY( iim, unsddu, 1, sdd1, 1 )
161	! CALL SCOPY( iim, sddu, 1, sdd2, 1 )
162	! END IF
163
164	jdfil1 = 1
165	jffil1 = jfiltnv
166	jdfil2 = jfiltsv
167	jffil2 = jjm
168	END IF
169	END IF
170
171	DO hemisph = 1, 2
172
173	IF (hemisph==1) THEN
174	jdfil = jdfil1
175	jffil = jffil1
176	ELSE
177	jdfil = jdfil2
178	jffil = jffil2
179	END IF
180
181	DO l = 1, nbniv
182	DO j = jdfil, jffil
183	DO i = 1, iim
184	champ(i, j, l) = champ(i, j, l) * sdd12(i, sdd1_type) ! sdd1(i)
185	END DO
186	END DO
187	END DO
188
189	IF(hemisph == 1) THEN
190
191	IF(ifiltre == -2) THEN
192
193	DO j = jdfil, jffil
194	#ifdef BLAS
195	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
196	matrinvn(1,1,j), &
197	iim, champ(1,j,1), iip1*nlat, 0.0, &
198	eignq(1,j-jdfil+1,1), iim*nlat)
199	#else
200	eignq(:, j - jdfil + 1, :) &
201	= matmul(matrinvn(:, :, j), champ(:iim, j, :))
202	#endif
203	END DO
204
205	ELSE IF (griscal) THEN
206
207	DO j = jdfil, jffil
208	#ifdef BLAS
209	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
210	matriceun(1,1,j), &
211	iim, champ(1,j,1), iip1*nlat, 0.0, &
212	eignq(1,j-jdfil+1,1), iim*nlat)
213	#else
214	eignq(:, j - jdfil + 1, :) &
215	= matmul(matriceun(:, :, j), champ(:iim, j, :))
216	#endif
217	END DO
218
219	ELSE
220
221	DO j = jdfil, jffil
222	#ifdef BLAS
223	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
224	matricevn(1,1,j), &
225	iim, champ(1,j,1), iip1*nlat, 0.0, &
226	eignq(1,j-jdfil+1,1), iim*nlat)
227	#else
228	eignq(:, j - jdfil + 1, :) &
229	= matmul(matricevn(:, :, j), champ(:iim, j, :))
230	#endif
231	END DO
232
233	ENDIF
234
235	ELSE
236
237	IF(ifiltre == -2) THEN
238
239	DO j = jdfil, jffil
240	#ifdef BLAS
241	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
242	matrinvs(1,1,j-jfiltsu+1), &
243	iim, champ(1,j,1), iip1*nlat, 0.0, &
244	eignq(1,j-jdfil+1,1), iim*nlat)
245	#else
246	eignq(:, j - jdfil + 1, :) &
247	= matmul(matrinvs(:, :, j - jfiltsu + 1), &
248	champ(:iim, j, :))
249	#endif
250	END DO
251
252	ELSE IF (griscal) THEN
253
254	DO j = jdfil, jffil
255	#ifdef BLAS
256	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
257	matriceus(1,1,j-jfiltsu+1), &
258	iim, champ(1,j,1), iip1*nlat, 0.0, &
259	eignq(1,j-jdfil+1,1), iim*nlat)
260	#else
261	eignq(:, j - jdfil + 1, :) &
262	= matmul(matriceus(:, :, j - jfiltsu + 1), &
263	champ(:iim, j, :))
264	#endif
265	END DO
266
267	ELSE
268
269	DO j = jdfil, jffil
270	#ifdef BLAS
271	CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
272	matricevs(1,1,j-jfiltsv+1), &
273	iim, champ(1,j,1), iip1*nlat, 0.0, &
274	eignq(1,j-jdfil+1,1), iim*nlat)
275	#else
276	eignq(:, j - jdfil + 1, :) &
277	= matmul(matricevs(:, :, j - jfiltsv + 1), &
278	champ(:iim, j, :))
279	#endif
280	END DO
281
282	ENDIF
283
284	ENDIF
285
286	IF(ifiltre== 2) THEN
287
288	DO l = 1, nbniv
289	DO j = jdfil, jffil
290	DO i = 1, iim
291	champ(i, j, l) = &
292	(champ(i, j, l) + eignq(i, j - jdfil + 1, l)) &
293	* sdd12(i, sdd2_type) ! sdd2(i)
294	END DO
295	END DO
296	END DO
297
298	ELSE
299
300	DO l = 1, nbniv
301	DO j = jdfil, jffil
302	DO i = 1, iim
303	champ(i, j, l) = &
304	(champ(i, j, l) - eignq(i, j - jdfil + 1, l)) &
305	* sdd12(i, sdd2_type) ! sdd2(i)
306	END DO
307	END DO
308	END DO
309
310	ENDIF
311
312	DO l = 1, nbniv
313	DO j = jdfil, jffil
314	champ(iip1, j, l) = champ(1, j, l)
315	END DO
316	END DO
317
318	ENDDO
319
320	1111 FORMAT(//20x, 'ERREUR dans le dimensionnement du tableau CHAMP a&
321	& filtrer, sur la grille des scalaires'/)
322	2222 FORMAT(//20x, 'ERREUR dans le dimensionnement du tableau CHAMP a fi&
323	& ltrer, sur la grille de V ou de Z'/)
324	RETURN
325	END SUBROUTINE filtreg
326
327	SUBROUTINE inifgn(dv)
328	!
329	! ... H.Upadyaya , O.Sharma ...
330	!
331	USE lmdz_coefils, ONLY: sddv, sddu, unsddu, unsddv, eignfnv, eignfnu
332	IMPLICIT NONE
333	!
334	include "dimensions.h"
335	include "paramet.h"
336	include "comgeom.h"
337	!
338	REAL :: vec(iim, iim), vec1(iim, iim)
339	REAL :: dlonu(iim), dlonv(iim)
340	REAL :: du(iim), dv(iim), d(iim)
341	REAL :: pi
342	INTEGER :: i, j, k, imm1, nrot
343	EXTERNAL SSUM
344	REAL :: SSUM
345	!
346
347	imm1 = iim - 1
348	pi = 2. * ASIN(1.)
349	!
350	DO i = 1, iim
351	dlonu(i) = xprimu(i)
352	dlonv(i) = xprimv(i)
353	END DO
354
355	DO i = 1, iim
356	sddv(i) = SQRT(dlonv(i))
357	sddu(i) = SQRT(dlonu(i))
358	unsddu(i) = 1. / sddu(i)
359	unsddv(i) = 1. / sddv(i)
360	END DO
361	!
362	DO j = 1, iim
363	DO i = 1, iim
364	vec(i, j) = 0.
365	vec1(i, j) = 0.
366	eignfnv(i, j) = 0.
367	eignfnu(i, j) = 0.
368	END DO
369	END DO
370	!
371	!
372	eignfnv(1, 1) = -1.
373	eignfnv(iim, 1) = 1.
374	DO i = 1, imm1
375	eignfnv(i + 1, i + 1) = -1.
376	eignfnv(i, i + 1) = 1.
377	END DO
378	DO j = 1, iim
379	DO i = 1, iim
380	eignfnv(i, j) = eignfnv(i, j) / (sddu(i) * sddv(j))
381	END DO
382	END DO
383	DO j = 1, iim
384	DO i = 1, iim
385	eignfnu(i, j) = -eignfnv(j, i)
386	END DO
387	END DO
388	!
389	DO j = 1, iim
390	DO i = 1, iim
391	vec (i, j) = 0.0
392	vec1(i, j) = 0.0
393	DO k = 1, iim
394	vec (i, j) = vec(i, j) + eignfnu(i, k) * eignfnv(k, j)
395	vec1(i, j) = vec1(i, j) + eignfnv(i, k) * eignfnu(k, j)
396	ENDDO
397	ENDDO
398	ENDDO
399
400	!
401	CALL jacobi(vec, iim, iim, dv, eignfnv, nrot)
402	CALL acc(eignfnv, d, iim)
403	CALL eigen_sort(dv, eignfnv, iim, iim)
404	!
405	CALL jacobi(vec1, iim, iim, du, eignfnu, nrot)
406	CALL acc(eignfnu, d, iim)
407	CALL eigen_sort(du, eignfnu, iim, iim)
408
409	!c ancienne version avec appels IMSL
410	!
411	! CALL MXM(eignfnu,iim,eignfnv,iim,vec,iim)
412	! CALL MXM(eignfnv,iim,eignfnu,iim,vec1,iim)
413	! CALL EVCSF(iim,vec,iim,dv,eignfnv,iim)
414	! CALL acc(eignfnv,d,iim)
415	! CALL eigen(eignfnv,dv)
416	!
417	! CALL EVCSF(iim,vec1,iim,du,eignfnu,iim)
418	! CALL acc(eignfnu,d,iim)
419	! CALL eigen(eignfnu,du)
420
421	RETURN
422	END SUBROUTINE inifgn
423
424	SUBROUTINE JACOBI(A, N, NP, D, V, NROT)
425	IMPLICIT NONE
426	! Arguments:
427	integer, intent(in) :: N
428	integer, intent(in) :: NP
429	integer, intent(out) :: NROT
430	real, intent(inout) :: A(NP, NP)
431	real, intent(out) :: D(NP)
432	real, intent(out) :: V(NP, NP)
433
434	! local variables:
435	INTEGER :: IP, IQ, I, J
436	REAL :: SM, TRESH, G, H, T, THETA, C, S, TAU
437	REAL :: B(N)
438	REAL :: Z(N)
439
440	DO IP = 1, N
441	DO IQ = 1, N
442	V(IP, IQ) = 0.
443	ENDDO
444	V(IP, IP) = 1.
445	ENDDO
446	DO IP = 1, N
447	B(IP) = A(IP, IP)
448	D(IP) = B(IP)
449	Z(IP) = 0.
450	ENDDO
451	NROT = 0
452	DO I = 1, 50 ! 50? I suspect this should be NP
453	! but convergence is fast enough anyway
454	SM = 0.
455	DO IP = 1, N - 1
456	DO IQ = IP + 1, N
457	SM = SM + ABS(A(IP, IQ))
458	ENDDO
459	ENDDO
460	IF(SM==0.)RETURN
461	IF(I<4)THEN
462	TRESH = 0.2 * SM / N**2
463	ELSE
464	TRESH = 0.
465	ENDIF
466	DO IP = 1, N - 1
467	DO IQ = IP + 1, N
468	G = 100. * ABS(A(IP, IQ))
469	IF((I>4).AND.(ABS(D(IP)) + G==ABS(D(IP))) &
470	.AND.(ABS(D(IQ)) + G==ABS(D(IQ))))THEN
471	A(IP, IQ) = 0.
472	ELSE IF(ABS(A(IP, IQ))>TRESH)THEN
473	H = D(IQ) - D(IP)
474	IF(ABS(H) + G==ABS(H))THEN
475	T = A(IP, IQ) / H
476	ELSE
477	THETA = 0.5 * H / A(IP, IQ)
478	T = 1. / (ABS(THETA) + SQRT(1. + THETA**2))
479	IF(THETA<0.)T = -T
480	ENDIF
481	C = 1. / SQRT(1 + T**2)
482	S = T * C
483	TAU = S / (1. + C)
484	H = T * A(IP, IQ)
485	Z(IP) = Z(IP) - H
486	Z(IQ) = Z(IQ) + H
487	D(IP) = D(IP) - H
488	D(IQ) = D(IQ) + H
489	A(IP, IQ) = 0.
490	DO J = 1, IP - 1
491	G = A(J, IP)
492	H = A(J, IQ)
493	A(J, IP) = G - S * (H + G * TAU)
494	A(J, IQ) = H + S * (G - H * TAU)
495	ENDDO
496	DO J = IP + 1, IQ - 1
497	G = A(IP, J)
498	H = A(J, IQ)
499	A(IP, J) = G - S * (H + G * TAU)
500	A(J, IQ) = H + S * (G - H * TAU)
501	ENDDO
502	DO J = IQ + 1, N
503	G = A(IP, J)
504	H = A(IQ, J)
505	A(IP, J) = G - S * (H + G * TAU)
506	A(IQ, J) = H + S * (G - H * TAU)
507	ENDDO
508	DO J = 1, N
509	G = V(J, IP)
510	H = V(J, IQ)
511	V(J, IP) = G - S * (H + G * TAU)
512	V(J, IQ) = H + S * (G - H * TAU)
513	ENDDO
514	NROT = NROT + 1
515	ENDIF
516	ENDDO
517	ENDDO
518	DO IP = 1, N
519	B(IP) = B(IP) + Z(IP)
520	D(IP) = B(IP)
521	Z(IP) = 0.
522	ENDDO
523	ENDDO ! of DO I=1,50
524	STOP 'Jacobi: 50 iterations should never happen'
525
526	END SUBROUTINE JACOBI
527
528	SUBROUTINE eigen_sort(d, v, n, np)
529	INTEGER :: n, np
530	REAL :: d(np), v(np, np)
531	INTEGER :: i, j, k
532	REAL :: p
533
534	DO i = 1, n - 1
535	k = i
536	p = d(i)
537	DO j = i + 1, n
538	IF(d(j)>=p) THEN
539	k = j
540	p = d(j)
541	ENDIF
542	ENDDO
543
544	IF(k/=i) THEN
545	d(k) = d(i)
546	d(i) = p
547	DO j = 1, n
548	p = v(j, i)
549	v(j, i) = v(j, k)
550	v(j, k) = p
551	ENDDO
552	ENDIF
553	ENDDO
554
555	RETURN
556	END SUBROUTINE eigen_sort
557
558	SUBROUTINE acc(vec, d, im)
559	IMPLICIT NONE
560	INTEGER :: im
561	REAL :: vec(im, im), d(im)
562	INTEGER :: i, j
563	REAL :: sum
564	real, external :: ssum
565	do j = 1, im
566	do i = 1, im
567	d(i) = vec(i, j) * vec(i, j)
568	enddo
569	sum = ssum(im, d, 1)
570	sum = sqrt(sum)
571	do i = 1, im
572	vec(i, j) = vec(i, j) / sum
573	enddo
574	enddo
575	RETURN
576	END SUBROUTINE acc
577
578
579	SUBROUTINE inifilr
580	#ifdef CPP_PARA
581	USE mod_filtre_fft, ONLY: use_filtre_fft,Init_filtre_fft
582	USE mod_filtre_fft_loc, ONLY: Init_filtre_fft_loc=>Init_filtre_fft !
583	#endif
584	USE serre_mod, ONLY: alphax
585	USE logic_mod, ONLY: fxyhypb, ysinus
586	USE comconst_mod, ONLY: maxlatfilter
587	USE lmdz_coefils, ONLY: modfrstv, modfrstu, jfiltnu, jfiltnv, coefilu, coefilv, &
588	coefilu2, coefilv2, eignfnv, eignfnu, jfiltsu, jfiltsv
589
590	! ... H. Upfiltreg_modadhyaya, O.Sharma ...
591
592	! version 3 .....
593
594	! Correction le 28/10/97 P. Le Van .
595	! -------------------------------------------------------------------
596	include "dimensions.h"
597	include "paramet.h"
598	include "comgeom.h"
599
600	REAL dlonu(iim), dlatu(jjm)
601	REAL rlamda(iim), eignvl(iim)
602
603	REAL lamdamax, pi, cof
604	INTEGER i, j, modemax, imx, k, kf, ii
605	REAL dymin, dxmin, colat0
606	REAL eignft(iim, iim), coff
607
608	LOGICAL, SAVE :: first_call_inifilr = .TRUE.
609
610	INTEGER ISMIN
611	EXTERNAL ISMIN
612	INTEGER iymin
613	INTEGER ixmineq
614
615	! ------------------------------------------------------------
616	! This routine computes the eigenfunctions of the laplacien
617	! on the stretched grid, and the filtering coefficients
618
619	! We designate:
620	! eignfn eigenfunctions of the discrete laplacien
621	! eigenvl eigenvalues
622	! jfiltn indexof the last scalar line filtered in NH
623	! jfilts index of the first line filtered in SH
624	! modfrst index of the mode from WHERE modes are filtered
625	! modemax maximum number of modes ( im )
626	! coefil filtering coefficients ( lamda_max*COS(rlat)/lamda )
627	! sdd SQRT( dx )
628
629	! the modes are filtered from modfrst to modemax
630
631	!-----------------------------------------------------------
632
633	if (iim == 1) return ! No filtre in 2D y-z
634
635	pi = 2. * ASIN(1.)
636
637	DO i = 1, iim
638	dlonu(i) = xprimu(i)
639	ENDDO
640
641	CALL inifgn(eignvl)
642
643	PRINT *, 'inifilr: EIGNVL '
644	PRINT 250, eignvl
645	250 FORMAT(1x, 5e14.6)
646
647	! compute eigenvalues and eigenfunctions
648
649
650	!.................................................................
651
652	! compute the filtering coefficients for scalar lines and
653	! meridional wind v-lines
654
655	! we filter all those latitude lines WHERE coefil < 1
656	! NO FILTERING AT POLES
657
658	! colat0 is to be used when alpha (stretching coefficient)
659	! is set equal to zero for the regular grid CASE
660
661	! ....... Calcul de colat0 .........
662	! ..... colat0 = minimum de ( 0.5, min dy/ min dx ) ...
663
664	DO j = 1, jjm
665	dlatu(j) = rlatu(j) - rlatu(j + 1)
666	ENDDO
667
668	dxmin = dlonu(1)
669	DO i = 2, iim
670	dxmin = MIN(dxmin, dlonu(i))
671	ENDDO
672	dymin = dlatu(1)
673	DO j = 2, jjm
674	dymin = MIN(dymin, dlatu(j))
675	ENDDO
676
677	! For a regular grid, we want the filter to start at latitudes
678	! corresponding to lengths dx of the same size as dy (in terms
679	! of angles: dx=2*dy) => at colat0=0.5 (i.e. colatitude=30 degrees
680	! <=> latitude=60 degrees).
681	! Same idea for the zoomed grid: start filtering polewards as soon
682	! as length dx becomes of the same size as dy
683
684	! if maxlatfilter >0, prescribe the colat0 value from the .def files
685
686	IF (maxlatfilter < 0.) THEN
687
688	colat0 = MIN(0.5, dymin / dxmin)
689	! colat0 = 1.
690
691	IF(.NOT.fxyhypb.AND.ysinus) THEN
692	colat0 = 0.6
693	! ...... a revoir pour ysinus ! .......
694	alphax = 0.
695	ENDIF
696
697	ELSE
698
699	colat0 = (90.0 - maxlatfilter) / 180.0 * pi
700
701	ENDIF
702
703	PRINT 50, colat0, alphax
704	50 FORMAT(/15x, ' Inifilr colat0 alphax ', 2e16.7)
705
706	IF(alphax==1.) THEN
707	PRINT *, ' Inifilr alphax doit etre < a 1. Corriger '
708	STOP
709	ENDIF
710
711	lamdamax = iim / (pi * colat0 * (1. - alphax))
712
713	! ... Correction le 28/10/97 ( P.Le Van ) ..
714
715	DO i = 2, iim
716	rlamda(i) = lamdamax / SQRT(ABS(eignvl(i)))
717	ENDDO
718
719	DO j = 1, jjm
720	DO i = 1, iim
721	coefilu(i, j) = 0.0
722	coefilv(i, j) = 0.0
723	coefilu2(i, j) = 0.0
724	coefilv2(i, j) = 0.0
725	ENDDO
726	ENDDO
727
728	! ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....
729	! .........................................................
730
731	modemax = iim
732
733	!!!! imx = modemax - 4 * (modemax/iim)
734
735	imx = iim
736
737	PRINT *, 'inifilr: TRUNCATION AT ', imx
738
739	! Ehouarn: set up some defaults
740	jfiltnu = 2 ! avoid north pole
741	jfiltsu = jjm ! avoid south pole (which is at jjm+1)
742	jfiltnv = 1 ! NB: no poles on the V grid
743	jfiltsv = jjm
744
745	DO j = 2, jjm / 2 + 1
746	cof = COS(rlatu(j)) / colat0
747	IF (cof < 1.) THEN
748	IF(rlamda(imx) * COS(rlatu(j))<1.) THEN
749	jfiltnu = j
750	ENDIF
751	ENDIF
752
753	cof = COS(rlatu(jjp1 - j + 1)) / colat0
754	IF (cof < 1.) THEN
755	IF(rlamda(imx) * COS(rlatu(jjp1 - j + 1))<1.) THEN
756	jfiltsu = jjp1 - j + 1
757	ENDIF
758	ENDIF
759	ENDDO
760
761	DO j = 1, jjm / 2
762	cof = COS(rlatv(j)) / colat0
763	IF (cof < 1.) THEN
764	IF(rlamda(imx) * COS(rlatv(j))<1.) THEN
765	jfiltnv = j
766	ENDIF
767	ENDIF
768
769	cof = COS(rlatv(jjm - j + 1)) / colat0
770	IF (cof < 1.) THEN
771	IF(rlamda(imx) * COS(rlatv(jjm - j + 1))<1.) THEN
772	jfiltsv = jjm - j + 1
773	ENDIF
774	ENDIF
775	ENDDO
776
777	IF(jfiltnu> jjm / 2 + 1) THEN
778	PRINT *, ' jfiltnu en dehors des valeurs acceptables ', jfiltnu
779	STOP
780	ENDIF
781
782	IF(jfiltsu> jjm + 1) THEN
783	PRINT *, ' jfiltsu en dehors des valeurs acceptables ', jfiltsu
784	STOP
785	ENDIF
786
787	IF(jfiltnv> jjm / 2) THEN
788	PRINT *, ' jfiltnv en dehors des valeurs acceptables ', jfiltnv
789	STOP
790	ENDIF
791
792	IF(jfiltsv> jjm) THEN
793	PRINT *, ' jfiltsv en dehors des valeurs acceptables ', jfiltsv
794	STOP
795	ENDIF
796
797	PRINT *, 'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ', &
798	jfiltnv, jfiltsv, jfiltnu, jfiltsu
799
800	IF(first_call_inifilr) THEN
801	ALLOCATE(matriceun(iim, iim, jfiltnu))
802	ALLOCATE(matriceus(iim, iim, jjm - jfiltsu + 1))
803	ALLOCATE(matricevn(iim, iim, jfiltnv))
804	ALLOCATE(matricevs(iim, iim, jjm - jfiltsv + 1))
805	ALLOCATE(matrinvn(iim, iim, jfiltnu))
806	ALLOCATE(matrinvs(iim, iim, jjm - jfiltsu + 1))
807	first_call_inifilr = .FALSE.
808	ENDIF
809
810	! ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
811	!................................................................
812
813	DO j = 1, jjm
814	!default initialization: all modes are retained (i.e. no filtering)
815	modfrstu(j) = iim
816	modfrstv(j) = iim
817	ENDDO
818
819	DO j = 2, jfiltnu
820	DO k = 2, modemax
821	cof = rlamda(k) * COS(rlatu(j))
822	IF (cof < 1.) GOTO 82
823	ENDDO
824	GOTO 84
825	82 modfrstu(j) = k
826
827	kf = modfrstu(j)
828	DO k = kf, modemax
829	cof = rlamda(k) * COS(rlatu(j))
830	coefilu(k, j) = cof - 1.
831	coefilu2(k, j) = cof * cof - 1.
832	ENDDO
833	84 CONTINUE
834	ENDDO
835
836	DO j = 1, jfiltnv
837
838	DO k = 2, modemax
839	cof = rlamda(k) * COS(rlatv(j))
840	IF (cof < 1.) GOTO 87
841	ENDDO
842	GOTO 89
843	87 modfrstv(j) = k
844
845	kf = modfrstv(j)
846	DO k = kf, modemax
847	cof = rlamda(k) * COS(rlatv(j))
848	coefilv(k, j) = cof - 1.
849	coefilv2(k, j) = cof * cof - 1.
850	ENDDO
851	89 CONTINUE
852	ENDDO
853
854	DO j = jfiltsu, jjm
855	DO k = 2, modemax
856	cof = rlamda(k) * COS(rlatu(j))
857	IF (cof < 1.) GOTO 92
858	ENDDO
859	GOTO 94
860	92 modfrstu(j) = k
861
862	kf = modfrstu(j)
863	DO k = kf, modemax
864	cof = rlamda(k) * COS(rlatu(j))
865	coefilu(k, j) = cof - 1.
866	coefilu2(k, j) = cof * cof - 1.
867	ENDDO
868	94 CONTINUE
869	ENDDO
870
871	DO j = jfiltsv, jjm
872	DO k = 2, modemax
873	cof = rlamda(k) * COS(rlatv(j))
874	IF (cof < 1.) GOTO 97
875	ENDDO
876	GOTO 99
877	97 modfrstv(j) = k
878
879	kf = modfrstv(j)
880	DO k = kf, modemax
881	cof = rlamda(k) * COS(rlatv(j))
882	coefilv(k, j) = cof - 1.
883	coefilv2(k, j) = cof * cof - 1.
884	ENDDO
885	99 CONTINUE
886	ENDDO
887
888	IF(jfiltnv>=jjm / 2 .OR. jfiltnu>=jjm / 2)THEN
889	! Ehouarn: and what are these for??? Trying to handle a limit case
890	! where filters extend to and meet at the equator?
891	IF(jfiltnv==jfiltsv)jfiltsv = 1 + jfiltnv
892	IF(jfiltnu==jfiltsu)jfiltsu = 1 + jfiltnu
893
894	PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', &
895	jfiltnv, jfiltsv, jfiltnu, jfiltsu
896	ENDIF
897
898	PRINT *, ' Modes premiers v '
899	PRINT 334, modfrstv
900	PRINT *, ' Modes premiers u '
901	PRINT 334, modfrstu
902
903	! ...................................................................
904
905	! ... Calcul de la matrice filtre 'matriceu' pour les champs situes
906	! sur la grille scalaire ........
907	! ...................................................................
908
909	DO j = 2, jfiltnu
910
911	DO i = 1, iim
912	coff = coefilu(i, j)
913	IF(i<modfrstu(j)) coff = 0.
914	DO k = 1, iim
915	eignft(i, k) = eignfnv(k, i) * coff
916	ENDDO
917	ENDDO ! of DO i=1,iim
918
919	#ifdef BLAS
920	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
921	eignfnv, iim, eignft, iim, 0.0, matriceun(1,1,j), iim)
922	#else
923	DO k = 1, iim
924	DO i = 1, iim
925	matriceun(i, k, j) = 0.0
926	DO ii = 1, iim
927	matriceun(i, k, j) = matriceun(i, k, j) &
928	+ eignfnv(i, ii) * eignft(ii, k)
929	ENDDO
930	ENDDO
931	ENDDO ! of DO k = 1, iim
932	#endif
933
934	ENDDO ! of DO j = 2, jfiltnu
935
936	DO j = jfiltsu, jjm
937
938	DO i = 1, iim
939	coff = coefilu(i, j)
940	IF(i<modfrstu(j)) coff = 0.
941	DO k = 1, iim
942	eignft(i, k) = eignfnv(k, i) * coff
943	ENDDO
944	ENDDO ! of DO i=1,iim
945	#ifdef BLAS
946	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
947	eignfnv, iim, eignft, iim, 0.0, &
948	matriceus(1,1,j-jfiltsu+1), iim)
949	#else
950	DO k = 1, iim
951	DO i = 1, iim
952	matriceus(i, k, j - jfiltsu + 1) = 0.0
953	DO ii = 1, iim
954	matriceus(i, k, j - jfiltsu + 1) = matriceus(i, k, j - jfiltsu + 1) &
955	+ eignfnv(i, ii) * eignft(ii, k)
956	ENDDO
957	ENDDO
958	ENDDO ! of DO k = 1, iim
959	#endif
960
961	ENDDO ! of DO j = jfiltsu, jjm
962
963	! ...................................................................
964
965	! ... Calcul de la matrice filtre 'matricev' pour les champs situes
966	! sur la grille de V ou de Z ........
967	! ...................................................................
968
969	DO j = 1, jfiltnv
970
971	DO i = 1, iim
972	coff = coefilv(i, j)
973	IF(i<modfrstv(j)) coff = 0.
974	DO k = 1, iim
975	eignft(i, k) = eignfnu(k, i) * coff
976	ENDDO
977	ENDDO
978
979	#ifdef BLAS
980	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
981	eignfnu, iim, eignft, iim, 0.0, matricevn(1,1,j), iim)
982	#else
983	DO k = 1, iim
984	DO i = 1, iim
985	matricevn(i, k, j) = 0.0
986	DO ii = 1, iim
987	matricevn(i, k, j) = matricevn(i, k, j) &
988	+ eignfnu(i, ii) * eignft(ii, k)
989	ENDDO
990	ENDDO
991	ENDDO
992	#endif
993
994	ENDDO ! of DO j = 1, jfiltnv
995
996	DO j = jfiltsv, jjm
997
998	DO i = 1, iim
999	coff = coefilv(i, j)
1000	IF(i<modfrstv(j)) coff = 0.
1001	DO k = 1, iim
1002	eignft(i, k) = eignfnu(k, i) * coff
1003	ENDDO
1004	ENDDO
1005
1006	#ifdef BLAS
1007	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
1008	eignfnu, iim, eignft, iim, 0.0, &
1009	matricevs(1,1,j-jfiltsv+1), iim)
1010	#else
1011	DO k = 1, iim
1012	DO i = 1, iim
1013	matricevs(i, k, j - jfiltsv + 1) = 0.0
1014	DO ii = 1, iim
1015	matricevs(i, k, j - jfiltsv + 1) = matricevs(i, k, j - jfiltsv + 1) &
1016	+ eignfnu(i, ii) * eignft(ii, k)
1017	ENDDO
1018	ENDDO
1019	ENDDO
1020	#endif
1021
1022	ENDDO ! of DO j = jfiltsv, jjm
1023
1024	! ...................................................................
1025
1026	! ... Calcul de la matrice filtre 'matrinv' pour les champs situes
1027	! sur la grille scalaire , pour le filtre inverse ........
1028	! ...................................................................
1029
1030	DO j = 2, jfiltnu
1031
1032	DO i = 1, iim
1033	coff = coefilu(i, j) / (1. + coefilu(i, j))
1034	IF(i<modfrstu(j)) coff = 0.
1035	DO k = 1, iim
1036	eignft(i, k) = eignfnv(k, i) * coff
1037	ENDDO
1038	ENDDO
1039
1040	#ifdef BLAS
1041	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
1042	eignfnv, iim, eignft, iim, 0.0, matrinvn(1,1,j), iim)
1043	#else
1044	DO k = 1, iim
1045	DO i = 1, iim
1046	matrinvn(i, k, j) = 0.0
1047	DO ii = 1, iim
1048	matrinvn(i, k, j) = matrinvn(i, k, j) &
1049	+ eignfnv(i, ii) * eignft(ii, k)
1050	ENDDO
1051	ENDDO
1052	ENDDO
1053	#endif
1054
1055	ENDDO ! of DO j = 2, jfiltnu
1056
1057	DO j = jfiltsu, jjm
1058
1059	DO i = 1, iim
1060	coff = coefilu(i, j) / (1. + coefilu(i, j))
1061	IF(i<modfrstu(j)) coff = 0.
1062	DO k = 1, iim
1063	eignft(i, k) = eignfnv(k, i) * coff
1064	ENDDO
1065	ENDDO
1066	#ifdef BLAS
1067	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
1068	eignfnv, iim, eignft, iim, 0.0, matrinvs(1,1,j-jfiltsu+1), iim)
1069	#else
1070	DO k = 1, iim
1071	DO i = 1, iim
1072	matrinvs(i, k, j - jfiltsu + 1) = 0.0
1073	DO ii = 1, iim
1074	matrinvs(i, k, j - jfiltsu + 1) = matrinvs(i, k, j - jfiltsu + 1) &
1075	+ eignfnv(i, ii) * eignft(ii, k)
1076	ENDDO
1077	ENDDO
1078	ENDDO
1079	#endif
1080
1081	ENDDO ! of DO j = jfiltsu, jjm
1082
1083	#ifdef CPP_PARA
1084	IF (use_filtre_fft) THEN
1085	CALL Init_filtre_fft(coefilu,modfrstu,jfiltnu,jfiltsu, &
1086	coefilv,modfrstv,jfiltnv,jfiltsv)
1087	CALL Init_filtre_fft_loc(coefilu,modfrstu,jfiltnu,jfiltsu, &
1088	coefilv,modfrstv,jfiltnv,jfiltsv)
1089	ENDIF
1090	#endif
1091	! ...................................................................
1092
1093	334 FORMAT(1x, 24i3)
1094
1095	END SUBROUTINE inifilr
1096
1097	END MODULE lmdz_filtreg

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

Download in other formats:

Original Format