source: LMDZ6/branches/contrails/libf/filtrez/filtreg_mod.F90 @ 5445

Last change on this file since 5445 was 5297, checked in by abarral, 8 weeks ago

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