Context Navigation

filtreg_mod.F90 @ 3994

Last change on this file since 3994 was 1601, checked in by Ehouarn Millour, 13 years ago

Updates and upgrades of filter:

add minor corrections (r1591 of trunk) in filtreg_mod.F90: some arrays were oversized and the computation of the indexes from which the filter is applied could go wrong in extreme cases.
adapt filtreg_p.F (r1597 of trunk) so that BLAS routine DGEMM is only used if 'BLAS' preprocessing flag is set.
update FFT filter routines to match current 'trunk' versions (mostly cosmetic changes, exept for the implementation of use of FFTW in mod_fft_fftw.F90).
update "arch-PW6_VARGAS.fcm" to enable use of FFTW

NB: implemented FFTs assume working precision to be double precision ; trying to use them when working precision is single precision will clearly end in tragedy.
EM

File size: 14.7 KB

Rev	Line
[1601]	1	!
	2	! $Id $
	3	!
[1108]	4	MODULE filtreg_mod
[1086]	5
	6	REAL, DIMENSION(:,:,:), ALLOCATABLE :: matriceun,matriceus,matricevn
	7	REAL, DIMENSION(:,:,:), ALLOCATABLE :: matricevs,matrinvn,matrinvs
	8
	9	CONTAINS
	10
	11	SUBROUTINE inifilr
[1279]	12	USE mod_filtre_fft
[1086]	13	!
	14	! ... H. Upadhyaya, O.Sharma ...
	15	!
	16	IMPLICIT NONE
	17	!
	18	! version 3 .....
	19
	20	! Correction le 28/10/97 P. Le Van .
	21	! -------------------------------------------------------------------
	22	#include "dimensions.h"
	23	#include "paramet.h"
	24	! -------------------------------------------------------------------
	25	#include "comgeom.h"
	26	#include "coefils.h"
	27	#include "logic.h"
	28	#include "serre.h"
	29
	30	REAL dlonu(iim),dlatu(jjm)
	31	REAL rlamda( iim ), eignvl( iim )
	32	!
	33
	34	REAL lamdamax,pi,cof
	35	INTEGER i,j,modemax,imx,k,kf,ii
	36	REAL dymin,dxmin,colat0
	37	REAL eignft(iim,iim), coff
	38
	39	LOGICAL, SAVE :: first_call_inifilr = .TRUE.
	40
	41	#ifdef CRAY
	42	INTEGER ISMIN
	43	EXTERNAL ISMIN
	44	INTEGER iymin
	45	INTEGER ixmineq
	46	#endif
	47	!
	48	! ------------------------------------------------------------
	49	! This routine computes the eigenfunctions of the laplacien
	50	! on the stretched grid, and the filtering coefficients
	51	!
	52	! We designate:
	53	! eignfn eigenfunctions of the discrete laplacien
	54	! eigenvl eigenvalues
	55	! jfiltn indexof the last scalar line filtered in NH
	56	! jfilts index of the first line filtered in SH
	57	! modfrst index of the mode from WHERE modes are filtered
	58	! modemax maximum number of modes ( im )
	59	! coefil filtering coefficients ( lamda_max*COS(rlat)/lamda )
	60	! sdd SQRT( dx )
	61	!
	62	! the modes are filtered from modfrst to modemax
	63	!
	64	!-----------------------------------------------------------
	65	!
	66
	67	pi = 2. * ASIN( 1. )
	68
	69	DO i = 1,iim
	70	dlonu(i) = xprimu( i )
	71	ENDDO
	72	!
	73	CALL inifgn(eignvl)
	74	!
[1601]	75	PRINT *,'inifilr: EIGNVL '
[1086]	76	PRINT 250,eignvl
[1601]	77	250 FORMAT( 1x,5e14.6)
[1086]	78	!
	79	! compute eigenvalues and eigenfunctions
	80	!
	81	!
	82	!.................................................................
	83	!
	84	! compute the filtering coefficients for scalar lines and
	85	! meridional wind v-lines
	86	!
	87	! we filter all those latitude lines WHERE coefil < 1
	88	! NO FILTERING AT POLES
	89	!
	90	! colat0 is to be used when alpha (stretching coefficient)
	91	! is set equal to zero for the regular grid CASE
	92	!
	93	! ....... Calcul de colat0 .........
	94	! ..... colat0 = minimum de ( 0.5, min dy/ min dx ) ...
	95	!
	96	!
	97	DO j = 1,jjm
	98	dlatu( j ) = rlatu( j ) - rlatu( j+1 )
	99	ENDDO
	100	!
	101	#ifdef CRAY
	102	iymin = ISMIN( jjm, dlatu, 1 )
	103	ixmineq = ISMIN( iim, dlonu, 1 )
	104	dymin = dlatu( iymin )
	105	dxmin = dlonu( ixmineq )
	106	#else
	107	dxmin = dlonu(1)
	108	DO i = 2, iim
	109	dxmin = MIN( dxmin,dlonu(i) )
	110	ENDDO
	111	dymin = dlatu(1)
	112	DO j = 2, jjm
	113	dymin = MIN( dymin,dlatu(j) )
	114	ENDDO
	115	#endif
	116	!
[1601]	117	! For a regular grid, we want the filter to start at latitudes
	118	! corresponding to lengths dx of the same size as dy (in terms
	119	! of angles: dx=2*dy) => at colat0=0.5 (i.e. colatitude=30 degrees
	120	! <=> latitude=60 degrees).
	121	! Same idea for the zoomed grid: start filtering polewards as soon
	122	! as length dx becomes of the same size as dy
[1086]	123	!
	124	colat0 = MIN( 0.5, dymin/dxmin )
	125	!
	126	IF( .NOT.fxyhypb.AND.ysinus ) THEN
	127	colat0 = 0.6
	128	! ...... a revoir pour ysinus ! .......
	129	alphax = 0.
	130	ENDIF
	131	!
	132	PRINT 50, colat0,alphax
	133	50 FORMAT(/15x,' Inifilr colat0 alphax ',2e16.7)
	134	!
	135	IF(alphax.EQ.1. ) THEN
	136	PRINT *,' Inifilr alphax doit etre < a 1. Corriger '
	137	STOP
	138	ENDIF
	139	!
	140	lamdamax = iim / ( pi * colat0 * ( 1. - alphax ) )
	141
	142	! ... Correction le 28/10/97 ( P.Le Van ) ..
	143	!
	144	DO i = 2,iim
	145	rlamda( i ) = lamdamax/ SQRT( ABS( eignvl(i) ) )
	146	ENDDO
	147	!
	148
	149	DO j = 1,jjm
	150	DO i = 1,iim
	151	coefilu( i,j ) = 0.0
	152	coefilv( i,j ) = 0.0
	153	coefilu2( i,j ) = 0.0
	154	coefilv2( i,j ) = 0.0
	155	ENDDO
	156	ENDDO
	157
	158	!
	159	! ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....
	160	! .........................................................
	161	!
	162	modemax = iim
	163
	164	!!!! imx = modemax - 4 * (modemax/iim)
	165
	166	imx = iim
	167	!
[1601]	168	PRINT *,'inifilr: TRUNCATION AT ',imx
[1086]	169	!
[1601]	170	! Ehouarn: set up some defaults
	171	jfiltnu=2 ! avoid north pole
	172	jfiltsu=jjm ! avoid south pole (which is at jjm+1)
	173	jfiltnv=1 ! NB: no poles on the V grid
	174	jfiltsv=jjm
	175
[1086]	176	DO j = 2, jjm/2+1
	177	cof = COS( rlatu(j) )/ colat0
	178	IF ( cof .LT. 1. ) THEN
[1601]	179	IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) THEN
	180	jfiltnu= j
	181	ENDIF
[1086]	182	ENDIF
	183
	184	cof = COS( rlatu(jjp1-j+1) )/ colat0
	185	IF ( cof .LT. 1. ) THEN
[1601]	186	IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) THEN
[1086]	187	jfiltsu= jjp1-j+1
[1601]	188	ENDIF
[1086]	189	ENDIF
	190	ENDDO
	191	!
	192	DO j = 1, jjm/2
	193	cof = COS( rlatv(j) )/ colat0
	194	IF ( cof .LT. 1. ) THEN
[1601]	195	IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) THEN
	196	jfiltnv= j
	197	ENDIF
[1086]	198	ENDIF
	199
	200	cof = COS( rlatv(jjm-j+1) )/ colat0
	201	IF ( cof .LT. 1. ) THEN
[1601]	202	IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) THEN
[1086]	203	jfiltsv= jjm-j+1
[1601]	204	ENDIF
[1086]	205	ENDIF
	206	ENDDO
	207	!
	208
	209	IF( jfiltnu.GT. jjm/2 +1 ) THEN
	210	PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
	211	STOP
	212	ENDIF
	213
	214	IF( jfiltsu.GT. jjm +1 ) THEN
	215	PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
	216	STOP
	217	ENDIF
	218
	219	IF( jfiltnv.GT. jjm/2 ) THEN
	220	PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
	221	STOP
	222	ENDIF
	223
	224	IF( jfiltsv.GT. jjm ) THEN
	225	PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
	226	STOP
	227	ENDIF
	228
[1601]	229	PRINT *,'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ' , &
[1086]	230	jfiltnv,jfiltsv,jfiltnu,jfiltsu
	231
	232	IF(first_call_inifilr) THEN
	233	ALLOCATE(matriceun(iim,iim,jfiltnu))
[1601]	234	ALLOCATE(matriceus(iim,iim,jjm-jfiltsu+1))
[1086]	235	ALLOCATE(matricevn(iim,iim,jfiltnv))
[1601]	236	ALLOCATE(matricevs(iim,iim,jjm-jfiltsv+1))
[1086]	237	ALLOCATE( matrinvn(iim,iim,jfiltnu))
[1601]	238	ALLOCATE( matrinvs(iim,iim,jjm-jfiltsu+1))
[1086]	239	first_call_inifilr = .FALSE.
	240	ENDIF
	241
	242	!
	243	! ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
	244	!................................................................
	245	!
	246	!
	247	DO j = 1,jjm
[1601]	248	!default initialization: all modes are retained (i.e. no filtering)
[1086]	249	modfrstu( j ) = iim
	250	modfrstv( j ) = iim
	251	ENDDO
	252	!
	253	DO j = 2,jfiltnu
	254	DO k = 2,modemax
	255	cof = rlamda(k) * COS( rlatu(j) )
	256	IF ( cof .LT. 1. ) GOTO 82
	257	ENDDO
	258	GOTO 84
	259	82 modfrstu( j ) = k
	260	!
	261	kf = modfrstu( j )
	262	DO k = kf , modemax
	263	cof = rlamda(k) * COS( rlatu(j) )
	264	coefilu(k,j) = cof - 1.
	265	coefilu2(k,j) = cof*cof - 1.
	266	ENDDO
	267	84 CONTINUE
	268	ENDDO
	269	!
	270	!
	271	DO j = 1,jfiltnv
	272	!
	273	DO k = 2,modemax
	274	cof = rlamda(k) * COS( rlatv(j) )
	275	IF ( cof .LT. 1. ) GOTO 87
	276	ENDDO
	277	GOTO 89
	278	87 modfrstv( j ) = k
	279	!
	280	kf = modfrstv( j )
	281	DO k = kf , modemax
	282	cof = rlamda(k) * COS( rlatv(j) )
	283	coefilv(k,j) = cof - 1.
	284	coefilv2(k,j) = cof*cof - 1.
	285	ENDDO
	286	89 CONTINUE
	287	ENDDO
	288	!
	289	DO j = jfiltsu,jjm
	290	DO k = 2,modemax
	291	cof = rlamda(k) * COS( rlatu(j) )
	292	IF ( cof .LT. 1. ) GOTO 92
	293	ENDDO
	294	GOTO 94
	295	92 modfrstu( j ) = k
	296	!
	297	kf = modfrstu( j )
	298	DO k = kf , modemax
	299	cof = rlamda(k) * COS( rlatu(j) )
	300	coefilu(k,j) = cof - 1.
	301	coefilu2(k,j) = cof*cof - 1.
	302	ENDDO
	303	94 CONTINUE
	304	ENDDO
	305	!
	306	DO j = jfiltsv,jjm
	307	DO k = 2,modemax
	308	cof = rlamda(k) * COS( rlatv(j) )
	309	IF ( cof .LT. 1. ) GOTO 97
	310	ENDDO
	311	GOTO 99
	312	97 modfrstv( j ) = k
	313	!
	314	kf = modfrstv( j )
	315	DO k = kf , modemax
	316	cof = rlamda(k) * COS( rlatv(j) )
	317	coefilv(k,j) = cof - 1.
	318	coefilv2(k,j) = cof*cof - 1.
	319	ENDDO
	320	99 CONTINUE
	321	ENDDO
	322	!
	323
	324	IF(jfiltnv.GE.jjm/2 .OR. jfiltnu.GE.jjm/2)THEN
[1601]	325	! Ehouarn: and what are these for??? Trying to handle a limit case
	326	! where filters extend to and meet at the equator?
[1086]	327	IF(jfiltnv.EQ.jfiltsv)jfiltsv=1+jfiltnv
	328	IF(jfiltnu.EQ.jfiltsu)jfiltsu=1+jfiltnu
	329
	330	PRINT *,'jfiltnv jfiltsv jfiltnu jfiltsu' , &
	331	jfiltnv,jfiltsv,jfiltnu,jfiltsu
	332	ENDIF
	333
	334	PRINT *,' Modes premiers v '
	335	PRINT 334,modfrstv
	336	PRINT *,' Modes premiers u '
	337	PRINT 334,modfrstu
	338
	339	!
	340	! ...................................................................
	341	!
	342	! ... Calcul de la matrice filtre 'matriceu' pour les champs situes
	343	! sur la grille scalaire ........
	344	! ...................................................................
	345	!
	346	DO j = 2, jfiltnu
	347
	348	DO i=1,iim
	349	coff = coefilu(i,j)
	350	IF( i.LT.modfrstu(j) ) coff = 0.
	351	DO k=1,iim
	352	eignft(i,k) = eignfnv(k,i) * coff
	353	ENDDO
[1601]	354	ENDDO ! of DO i=1,iim
[1086]	355	#ifdef CRAY
	356	CALL MXM( eignfnv,iim,eignft,iim,matriceun(1,1,j),iim )
	357	#else
	358	#ifdef BLAS
	359	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
	360	eignfnv, iim, eignft, iim, 0.0, matriceun(1,1,j), iim)
	361	#else
	362	DO k = 1, iim
	363	DO i = 1, iim
	364	matriceun(i,k,j) = 0.0
	365	DO ii = 1, iim
	366	matriceun(i,k,j) = matriceun(i,k,j) &
	367	+ eignfnv(i,ii)*eignft(ii,k)
	368	ENDDO
	369	ENDDO
[1601]	370	ENDDO ! of DO k = 1, iim
[1086]	371	#endif
	372	#endif
	373
[1601]	374	ENDDO ! of DO j = 2, jfiltnu
[1086]	375
	376	DO j = jfiltsu, jjm
	377
	378	DO i=1,iim
	379	coff = coefilu(i,j)
	380	IF( i.LT.modfrstu(j) ) coff = 0.
	381	DO k=1,iim
	382	eignft(i,k) = eignfnv(k,i) * coff
	383	ENDDO
[1601]	384	ENDDO ! of DO i=1,iim
[1086]	385	#ifdef CRAY
	386	CALL MXM(eignfnv,iim,eignft,iim,matriceus(1,1,j-jfiltsu+1),iim)
	387	#else
	388	#ifdef BLAS
	389	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
	390	eignfnv, iim, eignft, iim, 0.0, &
	391	matriceus(1,1,j-jfiltsu+1), iim)
	392	#else
	393	DO k = 1, iim
	394	DO i = 1, iim
	395	matriceus(i,k,j-jfiltsu+1) = 0.0
	396	DO ii = 1, iim
	397	matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) &
	398	+ eignfnv(i,ii)*eignft(ii,k)
	399	ENDDO
	400	ENDDO
[1601]	401	ENDDO ! of DO k = 1, iim
[1086]	402	#endif
	403	#endif
	404
[1601]	405	ENDDO ! of DO j = jfiltsu, jjm
[1086]	406
	407	! ...................................................................
	408	!
	409	! ... Calcul de la matrice filtre 'matricev' pour les champs situes
	410	! sur la grille de V ou de Z ........
	411	! ...................................................................
	412	!
	413	DO j = 1, jfiltnv
	414
	415	DO i = 1, iim
	416	coff = coefilv(i,j)
	417	IF( i.LT.modfrstv(j) ) coff = 0.
	418	DO k = 1, iim
	419	eignft(i,k) = eignfnu(k,i) * coff
	420	ENDDO
	421	ENDDO
	422	#ifdef CRAY
	423	CALL MXM( eignfnu,iim,eignft,iim,matricevn(1,1,j),iim )
	424	#else
	425	#ifdef BLAS
	426	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
	427	eignfnu, iim, eignft, iim, 0.0, matricevn(1,1,j), iim)
	428	#else
	429	DO k = 1, iim
	430	DO i = 1, iim
	431	matricevn(i,k,j) = 0.0
	432	DO ii = 1, iim
	433	matricevn(i,k,j) = matricevn(i,k,j) &
	434	+ eignfnu(i,ii)*eignft(ii,k)
	435	ENDDO
	436	ENDDO
	437	ENDDO
	438	#endif
	439	#endif
	440
[1601]	441	ENDDO ! of DO j = 1, jfiltnv
[1086]	442
	443	DO j = jfiltsv, jjm
	444
	445	DO i = 1, iim
	446	coff = coefilv(i,j)
	447	IF( i.LT.modfrstv(j) ) coff = 0.
	448	DO k = 1, iim
	449	eignft(i,k) = eignfnu(k,i) * coff
	450	ENDDO
	451	ENDDO
	452	#ifdef CRAY
	453	CALL MXM(eignfnu,iim,eignft,iim,matricevs(1,1,j-jfiltsv+1),iim)
	454	#else
	455	#ifdef BLAS
	456	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
	457	eignfnu, iim, eignft, iim, 0.0, &
	458	matricevs(1,1,j-jfiltsv+1), iim)
	459	#else
	460	DO k = 1, iim
	461	DO i = 1, iim
	462	matricevs(i,k,j-jfiltsv+1) = 0.0
	463	DO ii = 1, iim
	464	matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) &
	465	+ eignfnu(i,ii)*eignft(ii,k)
	466	ENDDO
	467	ENDDO
	468	ENDDO
	469	#endif
	470	#endif
	471
[1601]	472	ENDDO ! of DO j = jfiltsv, jjm
[1086]	473
	474	! ...................................................................
	475	!
	476	! ... Calcul de la matrice filtre 'matrinv' pour les champs situes
	477	! sur la grille scalaire , pour le filtre inverse ........
	478	! ...................................................................
	479	!
	480	DO j = 2, jfiltnu
	481
	482	DO i = 1,iim
	483	coff = coefilu(i,j)/ ( 1. + coefilu(i,j) )
	484	IF( i.LT.modfrstu(j) ) coff = 0.
	485	DO k=1,iim
	486	eignft(i,k) = eignfnv(k,i) * coff
	487	ENDDO
	488	ENDDO
	489	#ifdef CRAY
	490	CALL MXM( eignfnv,iim,eignft,iim,matrinvn(1,1,j),iim )
	491	#else
	492	#ifdef BLAS
	493	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
	494	eignfnv, iim, eignft, iim, 0.0, matrinvn(1,1,j), iim)
	495	#else
	496	DO k = 1, iim
	497	DO i = 1, iim
	498	matrinvn(i,k,j) = 0.0
	499	DO ii = 1, iim
	500	matrinvn(i,k,j) = matrinvn(i,k,j) &
	501	+ eignfnv(i,ii)*eignft(ii,k)
	502	ENDDO
	503	ENDDO
	504	ENDDO
	505	#endif
	506	#endif
	507
[1601]	508	ENDDO ! of DO j = 2, jfiltnu
[1086]	509
	510	DO j = jfiltsu, jjm
	511
	512	DO i = 1,iim
	513	coff = coefilu(i,j) / ( 1. + coefilu(i,j) )
	514	IF( i.LT.modfrstu(j) ) coff = 0.
	515	DO k=1,iim
	516	eignft(i,k) = eignfnv(k,i) * coff
	517	ENDDO
	518	ENDDO
	519	#ifdef CRAY
	520	CALL MXM(eignfnv,iim,eignft,iim,matrinvs(1,1,j-jfiltsu+1),iim)
	521	#else
	522	#ifdef BLAS
	523	CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
	524	eignfnv, iim, eignft, iim, 0.0, matrinvs(1,1,j-jfiltsu+1), iim)
	525	#else
	526	DO k = 1, iim
	527	DO i = 1, iim
	528	matrinvs(i,k,j-jfiltsu+1) = 0.0
	529	DO ii = 1, iim
	530	matrinvs(i,k,j-jfiltsu+1) = matrinvs(i,k,j-jfiltsu+1) &
	531	+ eignfnv(i,ii)*eignft(ii,k)
	532	ENDDO
	533	ENDDO
	534	ENDDO
	535	#endif
	536	#endif
	537
[1601]	538	ENDDO ! of DO j = jfiltsu, jjm
[1086]	539
[1279]	540	IF (use_filtre_fft) THEN
	541	CALL Init_filtre_fft(coefilu,modfrstu,jfiltnu,jfiltsu, &
	542	coefilv,modfrstv,jfiltnv,jfiltsv)
	543	ENDIF
	544
[1086]	545	! ...................................................................
	546
	547	!
	548	334 FORMAT(1x,24i3)
	549	755 FORMAT(1x,6f10.3,i3)
	550
	551	RETURN
	552	END SUBROUTINE inifilr
	553
[1108]	554	END MODULE filtreg_mod

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

Download in other formats:

Original Format