source: LMDZ6/branches/Amaury_dev/libf/filtrez/filtreg_mod.F90 @ 5099

Last change on this file since 5099 was 5099, checked in by abarral, 2 months ago

Replace most uses of CPP_DUST by the corresponding logical defined in lmdz_cppkeys_wrapper.F90
Convert several files from .F to .f90 to allow Dust to compile w/o rrtm/ecrad
Create lmdz_yoerad.f90
(lint) Remove "!" on otherwise empty line
  • 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: 14.1 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    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    INTEGER   ISMIN
45    EXTERNAL  ISMIN
46    INTEGER iymin
47    INTEGER ixmineq
48
49    ! ------------------------------------------------------------
50    !   This routine computes the eigenfunctions of the laplacien
51    !   on the stretched grid, and the filtering coefficients
52
53    !  We designate:
54    !   eignfn   eigenfunctions of the discrete laplacien
55    !   eigenvl  eigenvalues
56    !   jfiltn   indexof the last scalar line filtered in NH
57    !   jfilts   index of the first line filtered in SH
58    !   modfrst  index of the mode from WHERE modes are filtered
59    !   modemax  maximum number of modes ( im )
60    !   coefil   filtering coefficients ( lamda_max*COS(rlat)/lamda )
61    !   sdd      SQRT( dx )
62
63    !     the modes are filtered from modfrst to modemax
64
65    !-----------------------------------------------------------
66
67    if ( iim == 1 ) return ! No filtre in 2D y-z
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
77    PRINT *,'inifilr: EIGNVL '
78    PRINT 250,eignvl
79250 FORMAT( 1x,5e14.6)
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    dxmin   =  dlonu(1)
104    DO  i  = 2, iim
105       dxmin = MIN( dxmin,dlonu(i) )
106    ENDDO
107    dymin  = dlatu(1)
108    DO j  = 2, jjm
109       dymin = MIN( dymin,dlatu(j) )
110    ENDDO
111
112    ! For a regular grid, we want the filter to start at latitudes
113    ! corresponding to lengths dx of the same size as dy (in terms
114    ! of angles: dx=2*dy) => at colat0=0.5 (i.e. colatitude=30 degrees
115    !  <=> latitude=60 degrees).
116    ! Same idea for the zoomed grid: start filtering polewards as soon
117    ! as length dx becomes of the same size as dy
118
119    ! if maxlatfilter >0, prescribe the colat0 value from the .def files
120   
121    IF (maxlatfilter < 0.) THEN
122
123    colat0  =  MIN( 0.5, dymin/dxmin )
124    ! colat0  =  1.
125
126    IF( .NOT.fxyhypb.AND.ysinus )  THEN
127       colat0 = 0.6
128       !         ...... a revoir  pour  ysinus !   .......
129       alphax = 0.
130    ENDIF
131
132    ELSE
133
134    colat0=(90.0-maxlatfilter)/180.0*pi       
135
136    ENDIF
137
138    PRINT 50, colat0,alphax
13950  FORMAT(/15x,' Inifilr colat0 alphax ',2e16.7)
140
141    IF(alphax==1. )  THEN
142       PRINT *,' Inifilr  alphax doit etre  <  a 1.  Corriger '
143       STOP
144    ENDIF
145
146    lamdamax = iim / ( pi * colat0 * ( 1. - alphax ) )
147
148    !                        ... Correction  le 28/10/97  ( P.Le Van ) ..
149
150    DO i = 2,iim
151       rlamda( i ) = lamdamax/ SQRT( ABS( eignvl(i) ) )
152    ENDDO
153
154    DO j = 1,jjm
155       DO i = 1,iim
156          coefilu( i,j )  = 0.0
157          coefilv( i,j )  = 0.0
158          coefilu2( i,j ) = 0.0
159          coefilv2( i,j ) = 0.0
160       ENDDO
161    ENDDO
162
163    !    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....
164    !    .........................................................
165
166    modemax = iim
167
168!!!!    imx = modemax - 4 * (modemax/iim)
169
170    imx  = iim
171
172    PRINT *,'inifilr: TRUNCATION AT ',imx
173
174! Ehouarn: set up some defaults
175    jfiltnu=2 ! avoid north pole
176    jfiltsu=jjm ! avoid south pole (which is at jjm+1)
177    jfiltnv=1 ! NB: no poles on the V grid
178    jfiltsv=jjm
179
180    DO j = 2, jjm/2+1
181       cof = COS( rlatu(j) )/ colat0
182       IF ( cof < 1. ) THEN
183          IF( rlamda(imx) * COS(rlatu(j) )<1. ) THEN
184            jfiltnu= j
185          ENDIF
186       ENDIF
187
188       cof = COS( rlatu(jjp1-j+1) )/ colat0
189       IF ( cof < 1. ) THEN
190          IF( rlamda(imx) * COS(rlatu(jjp1-j+1) )<1. ) THEN
191               jfiltsu= jjp1-j+1
192          ENDIF
193       ENDIF
194    ENDDO
195
196    DO j = 1, jjm/2
197       cof = COS( rlatv(j) )/ colat0
198       IF ( cof < 1. ) THEN
199          IF( rlamda(imx) * COS(rlatv(j) )<1. ) THEN
200            jfiltnv= j
201          ENDIF
202       ENDIF
203
204       cof = COS( rlatv(jjm-j+1) )/ colat0
205       IF ( cof < 1. ) THEN
206          IF( rlamda(imx) * COS(rlatv(jjm-j+1) )<1. ) THEN
207               jfiltsv= jjm-j+1
208          ENDIF
209       ENDIF
210    ENDDO
211
212    IF( jfiltnu> jjm/2 +1 )  THEN
213       PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
214       STOP
215    ENDIF
216
217    IF( jfiltsu>  jjm  +1 )  THEN
218       PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
219       STOP
220    ENDIF
221
222    IF( jfiltnv> jjm/2    )  THEN
223       PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
224       STOP
225    ENDIF
226
227    IF( jfiltsv>     jjm  )  THEN
228       PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
229       STOP
230    ENDIF
231
232    PRINT *,'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ' , &
233         jfiltnv,jfiltsv,jfiltnu,jfiltsu
234
235    IF(first_call_inifilr) THEN
236       ALLOCATE(matriceun(iim,iim,jfiltnu))
237       ALLOCATE(matriceus(iim,iim,jjm-jfiltsu+1))
238       ALLOCATE(matricevn(iim,iim,jfiltnv))
239       ALLOCATE(matricevs(iim,iim,jjm-jfiltsv+1))
240       ALLOCATE( matrinvn(iim,iim,jfiltnu))
241       ALLOCATE( matrinvs(iim,iim,jjm-jfiltsu+1))
242       first_call_inifilr = .FALSE.
243    ENDIF
244
245    !   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
246    !................................................................
247
248
249    DO j = 1,jjm
250    !default initialization: all modes are retained (i.e. no filtering)
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 < 1. ) GOTO 82
259       ENDDO
260       GOTO 84
26182     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
26984     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 < 1. ) GOTO 87
278       ENDDO
279       GOTO 89
28087     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
28889     CONTINUE
289    ENDDO
290
291    DO j = jfiltsu,jjm
292       DO k = 2,modemax
293          cof = rlamda(k) * COS( rlatu(j) )
294          IF ( cof < 1. ) GOTO 92
295       ENDDO
296       GOTO 94
29792     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
30594     CONTINUE
306    ENDDO
307
308    DO j = jfiltsv,jjm
309       DO k = 2,modemax
310          cof = rlamda(k) * COS( rlatv(j) )
311          IF ( cof < 1. ) GOTO 97
312       ENDDO
313       GOTO 99
31497     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
32299     CONTINUE
323    ENDDO
324
325    IF(jfiltnv>=jjm/2 .OR. jfiltnu>=jjm/2)THEN
326! Ehouarn: and what are these for??? Trying to handle a limit case
327!          where filters extend to and meet at the equator?
328       IF(jfiltnv==jfiltsv)jfiltsv=1+jfiltnv
329       IF(jfiltnu==jfiltsu)jfiltsu=1+jfiltnu
330
331       PRINT *,'jfiltnv jfiltsv jfiltnu jfiltsu' , &
332            jfiltnv,jfiltsv,jfiltnu,jfiltsu
333    ENDIF
334
335    PRINT *,'   Modes premiers  v  '
336    PRINT 334,modfrstv
337    PRINT *,'   Modes premiers  u  '
338    PRINT 334,modfrstu
339
340    !   ...................................................................
341
342    !   ... Calcul de la matrice filtre 'matriceu'  pour les champs situes
343    !                       sur la grille scalaire                 ........
344    !   ...................................................................
345
346    DO j = 2, jfiltnu
347
348       DO i=1,iim
349          coff = coefilu(i,j)
350          IF( i<modfrstu(j) ) coff = 0.
351          DO k=1,iim
352             eignft(i,k) = eignfnv(k,i) * coff
353          ENDDO
354       ENDDO ! of DO i=1,iim
355
356#ifdef BLAS
357       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
358            eignfnv, iim, eignft, iim, 0.0, matriceun(1,1,j), iim)
359#else
360       DO k = 1, iim
361          DO i = 1, iim
362             matriceun(i,k,j) = 0.0
363             DO ii = 1, iim
364                matriceun(i,k,j) = matriceun(i,k,j) &
365                     + eignfnv(i,ii)*eignft(ii,k)
366             ENDDO
367          ENDDO
368       ENDDO ! of DO k = 1, iim
369#endif
370
371    ENDDO ! of DO j = 2, jfiltnu
372
373    DO j = jfiltsu, jjm
374
375       DO i=1,iim
376          coff = coefilu(i,j)
377          IF( i<modfrstu(j) ) coff = 0.
378          DO k=1,iim
379             eignft(i,k) = eignfnv(k,i) * coff
380          ENDDO
381       ENDDO ! of DO i=1,iim
382#ifdef BLAS
383       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
384            eignfnv, iim, eignft, iim, 0.0, &
385            matriceus(1,1,j-jfiltsu+1), iim)
386#else
387       DO k = 1, iim
388          DO i = 1, iim
389             matriceus(i,k,j-jfiltsu+1) = 0.0
390             DO ii = 1, iim
391                matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) &
392                     + eignfnv(i,ii)*eignft(ii,k)
393             ENDDO
394          ENDDO
395       ENDDO ! of DO k = 1, iim
396#endif
397
398    ENDDO ! of DO j = jfiltsu, jjm
399
400    !   ...................................................................
401
402    !   ... Calcul de la matrice filtre 'matricev'  pour les champs situes
403    !                       sur la grille   de V ou de Z           ........
404    !   ...................................................................
405
406    DO j = 1, jfiltnv
407
408       DO i = 1, iim
409          coff = coefilv(i,j)
410          IF( i<modfrstv(j) ) coff = 0.
411          DO k = 1, iim
412             eignft(i,k) = eignfnu(k,i) * coff
413          ENDDO
414       ENDDO
415
416#ifdef BLAS
417       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
418            eignfnu, iim, eignft, iim, 0.0, matricevn(1,1,j), iim)
419#else
420       DO k = 1, iim
421          DO i = 1, iim
422             matricevn(i,k,j) = 0.0
423             DO ii = 1, iim
424                matricevn(i,k,j) = matricevn(i,k,j) &
425                     + eignfnu(i,ii)*eignft(ii,k)
426             ENDDO
427          ENDDO
428       ENDDO
429#endif
430
431    ENDDO ! of DO j = 1, jfiltnv
432
433    DO j = jfiltsv, jjm
434
435       DO i = 1, iim
436          coff = coefilv(i,j)
437          IF( i<modfrstv(j) ) coff = 0.
438          DO k = 1, iim
439             eignft(i,k) = eignfnu(k,i) * coff
440          ENDDO
441       ENDDO
442
443#ifdef BLAS
444       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
445            eignfnu, iim, eignft, iim, 0.0, &
446            matricevs(1,1,j-jfiltsv+1), iim)
447#else
448       DO k = 1, iim
449          DO i = 1, iim
450             matricevs(i,k,j-jfiltsv+1) = 0.0
451             DO ii = 1, iim
452                matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) &
453                     + eignfnu(i,ii)*eignft(ii,k)
454             ENDDO
455          ENDDO
456       ENDDO
457#endif
458
459    ENDDO ! of DO j = jfiltsv, jjm
460
461    !   ...................................................................
462
463    !   ... Calcul de la matrice filtre 'matrinv'  pour les champs situes
464    !              sur la grille scalaire , pour le filtre inverse ........
465    !   ...................................................................
466
467    DO j = 2, jfiltnu
468
469       DO i = 1,iim
470          coff = coefilu(i,j)/ ( 1. + coefilu(i,j) )
471          IF( i<modfrstu(j) ) coff = 0.
472          DO k=1,iim
473             eignft(i,k) = eignfnv(k,i) * coff
474          ENDDO
475       ENDDO
476
477#ifdef BLAS
478       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
479            eignfnv, iim, eignft, iim, 0.0, matrinvn(1,1,j), iim)
480#else
481       DO k = 1, iim
482          DO i = 1, iim
483             matrinvn(i,k,j) = 0.0
484             DO ii = 1, iim
485                matrinvn(i,k,j) = matrinvn(i,k,j) &
486                     + eignfnv(i,ii)*eignft(ii,k)
487             ENDDO
488          ENDDO
489       ENDDO
490#endif
491
492    ENDDO ! of DO j = 2, jfiltnu
493
494    DO j = jfiltsu, jjm
495
496       DO i = 1,iim
497          coff = coefilu(i,j) / ( 1. + coefilu(i,j) )
498          IF( i<modfrstu(j) ) coff = 0.
499          DO k=1,iim
500             eignft(i,k) = eignfnv(k,i) * coff
501          ENDDO
502       ENDDO
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
518    ENDDO ! of DO j = jfiltsu, jjm
519
520#ifdef CPP_PARA
521    IF (use_filtre_fft) THEN
522       CALL Init_filtre_fft(coefilu,modfrstu,jfiltnu,jfiltsu,  &
523                           coefilv,modfrstv,jfiltnv,jfiltsv)
524       CALL Init_filtre_fft_loc(coefilu,modfrstu,jfiltnu,jfiltsu,  &
525                           coefilv,modfrstv,jfiltnv,jfiltsv)
526    ENDIF
527#endif
528    !   ...................................................................
529
530334 FORMAT(1x,24i3)
531
532  END SUBROUTINE inifilr
533
534END MODULE filtreg_mod
Note: See TracBrowser for help on using the repository browser.