Context Navigation

filtreg_mod.F90 @ 2069

Last change on this file since 2069 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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: 15.0 KB

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

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

Download in other formats:

Original Format