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

Last change on this file since 5285 was 5285, checked in by abarral, 4 days ago

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