source: LMDZ6/trunk/libf/filtrez/filtreg_mod.F90

Last change on this file was 4519, checked in by evignon, 14 months ago

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