source: LMDZ4/trunk/libf/filtrez/filtreg_mod.F90 @ 1152

Last change on this file since 1152 was 1146, checked in by Laurent Fairhead, 16 years ago

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