Context Navigation

lmdz_filtreg.F90 @ 5456

Last change on this file since 5456 was 5159, checked in by abarral, 6 months ago
Put dimensions.h and paramet.h into modules
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

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

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

Download in other formats:

Original Format