GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: filtrez/filtreg_mod.F90 Lines: 185 202 91.6 %
Date: 2023-06-30 12:56:34 Branches: 170 210 81.0 %

Line Branch Exec Source
1
!
2
! $Id$
3
!
4
MODULE filtreg_mod
5
6
  REAL, DIMENSION(:,:,:), ALLOCATABLE :: matriceun,matriceus,matricevn
7
  REAL, DIMENSION(:,:,:), ALLOCATABLE :: matricevs,matrinvn,matrinvs
8
9
CONTAINS
10
11
1
  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
33
    DO i = 1,iim
75
33
       dlonu(i) = xprimu( i )
76
    ENDDO
77
    !
78
1
    CALL inifgn(eignvl)
79
    !
80
1
    PRINT *,'inifilr: EIGNVL '
81
1
    PRINT 250,eignvl
82
250 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
33
    DO j = 1,jjm
103
33
       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
1
    dxmin   =  dlonu(1)
113
32
    DO  i  = 2, iim
114
32
       dxmin = MIN( dxmin,dlonu(i) )
115
    ENDDO
116
1
    dymin  = dlatu(1)
117
32
    DO j  = 2, jjm
118
32
       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
1
    IF (maxlatfilter .LT. 0.) THEN
132
133
1
    colat0  =  MIN( 0.5, dymin/dxmin )
134
    ! colat0  =  1.
135
    !
136

1
    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
1
    PRINT 50, colat0,alphax
152
50  FORMAT(/15x,' Inifilr colat0 alphax ',2e16.7)
153
    !
154
1
    IF(alphax.EQ.1. )  THEN
155
       PRINT *,' Inifilr  alphax doit etre  <  a 1.  Corriger '
156
       STOP
157
    ENDIF
158
    !
159
1
    lamdamax = iim / ( pi * colat0 * ( 1. - alphax ) )
160
161
    !                        ... Correction  le 28/10/97  ( P.Le Van ) ..
162
    !
163
32
    DO i = 2,iim
164
32
       rlamda( i ) = lamdamax/ SQRT( ABS( eignvl(i) ) )
165
    ENDDO
166
    !
167
168
33
    DO j = 1,jjm
169
1057
       DO i = 1,iim
170
1024
          coefilu( i,j )  = 0.0
171
1024
          coefilv( i,j )  = 0.0
172
1024
          coefilu2( i,j ) = 0.0
173
1056
          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
1
    imx  = iim
186
    !
187
1
    PRINT *,'inifilr: TRUNCATION AT ',imx
188
    !
189
! Ehouarn: set up some defaults
190
1
    jfiltnu=2 ! avoid north pole
191
1
    jfiltsu=jjm ! avoid south pole (which is at jjm+1)
192
1
    jfiltnv=1 ! NB: no poles on the V grid
193
1
    jfiltsv=jjm
194
195
17
    DO j = 2, jjm/2+1
196
16
       cof = COS( rlatu(j) )/ colat0
197
16
       IF ( cof .LT. 1. ) THEN
198
5
          IF( rlamda(imx) * COS(rlatu(j) ).LT.1. ) THEN
199
5
            jfiltnu= j
200
          ENDIF
201
       ENDIF
202
203
16
       cof = COS( rlatu(jjp1-j+1) )/ colat0
204
17
       IF ( cof .LT. 1. ) THEN
205
5
          IF( rlamda(imx) * COS(rlatu(jjp1-j+1) ).LT.1. ) THEN
206
5
               jfiltsu= jjp1-j+1
207
          ENDIF
208
       ENDIF
209
    ENDDO
210
    !
211
17
    DO j = 1, jjm/2
212
16
       cof = COS( rlatv(j) )/ colat0
213
16
       IF ( cof .LT. 1. ) THEN
214
5
          IF( rlamda(imx) * COS(rlatv(j) ).LT.1. ) THEN
215
5
            jfiltnv= j
216
          ENDIF
217
       ENDIF
218
219
16
       cof = COS( rlatv(jjm-j+1) )/ colat0
220
17
       IF ( cof .LT. 1. ) THEN
221
5
          IF( rlamda(imx) * COS(rlatv(jjm-j+1) ).LT.1. ) THEN
222
5
               jfiltsv= jjm-j+1
223
          ENDIF
224
       ENDIF
225
    ENDDO
226
    !
227
228
1
    IF( jfiltnu.GT. jjm/2 +1 )  THEN
229
       PRINT *,' jfiltnu en dehors des valeurs acceptables ' ,jfiltnu
230
       STOP
231
    ENDIF
232
233
1
    IF( jfiltsu.GT.  jjm  +1 )  THEN
234
       PRINT *,' jfiltsu en dehors des valeurs acceptables ' ,jfiltsu
235
       STOP
236
    ENDIF
237
238
1
    IF( jfiltnv.GT. jjm/2    )  THEN
239
       PRINT *,' jfiltnv en dehors des valeurs acceptables ' ,jfiltnv
240
       STOP
241
    ENDIF
242
243
1
    IF( jfiltsv.GT.     jjm  )  THEN
244
       PRINT *,' jfiltsv en dehors des valeurs acceptables ' ,jfiltsv
245
       STOP
246
    ENDIF
247
248
1
    PRINT *,'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ' , &
249
2
         jfiltnv,jfiltsv,jfiltnu,jfiltsu
250
251
1
    IF(first_call_inifilr) THEN
252

1
       ALLOCATE(matriceun(iim,iim,jfiltnu))
253

1
       ALLOCATE(matriceus(iim,iim,jjm-jfiltsu+1))
254

1
       ALLOCATE(matricevn(iim,iim,jfiltnv))
255

1
       ALLOCATE(matricevs(iim,iim,jjm-jfiltsv+1))
256

1
       ALLOCATE( matrinvn(iim,iim,jfiltnu))
257

1
       ALLOCATE( matrinvs(iim,iim,jjm-jfiltsu+1))
258
1
       first_call_inifilr = .FALSE.
259
    ENDIF
260
261
    !
262
    !   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
263
    !................................................................
264
    !
265
    !
266
33
    DO j = 1,jjm
267
    !default initialization: all modes are retained (i.e. no filtering)
268
32
       modfrstu( j ) = iim
269
33
       modfrstv( j ) = iim
270
    ENDDO
271
    !
272
6
    DO j = 2,jfiltnu
273
69
       DO k = 2,modemax
274
69
          cof = rlamda(k) * COS( rlatu(j) )
275
69
          IF ( cof .LT. 1. ) GOTO 82
276
       ENDDO
277
5
       GOTO 84
278
5
82     modfrstu( j ) = k
279
       !
280
       kf = modfrstu( j )
281
96
       DO k = kf , modemax
282
91
          cof = rlamda(k) * COS( rlatu(j) )
283
91
          coefilu(k,j) = cof - 1.
284
96
          coefilu2(k,j) = cof*cof - 1.
285
       ENDDO
286
1
84     CONTINUE
287
    ENDDO
288
    !
289
    !
290
6
    DO j = 1,jfiltnv
291
       !
292
57
       DO k = 2,modemax
293
57
          cof = rlamda(k) * COS( rlatv(j) )
294
57
          IF ( cof .LT. 1. ) GOTO 87
295
       ENDDO
296
5
       GOTO 89
297
5
87     modfrstv( j ) = k
298
       !
299
       kf = modfrstv( j )
300
108
       DO k = kf , modemax
301
103
          cof = rlamda(k) * COS( rlatv(j) )
302
103
          coefilv(k,j) = cof - 1.
303
108
          coefilv2(k,j) = cof*cof - 1.
304
       ENDDO
305
1
89     CONTINUE
306
    ENDDO
307
    !
308
6
    DO j = jfiltsu,jjm
309
69
       DO k = 2,modemax
310
69
          cof = rlamda(k) * COS( rlatu(j) )
311
69
          IF ( cof .LT. 1. ) GOTO 92
312
       ENDDO
313
5
       GOTO 94
314
5
92     modfrstu( j ) = k
315
       !
316
       kf = modfrstu( j )
317
96
       DO k = kf , modemax
318
91
          cof = rlamda(k) * COS( rlatu(j) )
319
91
          coefilu(k,j) = cof - 1.
320
96
          coefilu2(k,j) = cof*cof - 1.
321
       ENDDO
322
1
94     CONTINUE
323
    ENDDO
324
    !
325
6
    DO j = jfiltsv,jjm
326
57
       DO k = 2,modemax
327
57
          cof = rlamda(k) * COS( rlatv(j) )
328
57
          IF ( cof .LT. 1. ) GOTO 97
329
       ENDDO
330
5
       GOTO 99
331
5
97     modfrstv( j ) = k
332
       !
333
       kf = modfrstv( j )
334
108
       DO k = kf , modemax
335
103
          cof = rlamda(k) * COS( rlatv(j) )
336
103
          coefilv(k,j) = cof - 1.
337
108
          coefilv2(k,j) = cof*cof - 1.
338
       ENDDO
339
1
99     CONTINUE
340
    ENDDO
341
    !
342
343

1
    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
1
    PRINT *,'   Modes premiers  v  '
354
1
    PRINT 334,modfrstv
355
1
    PRINT *,'   Modes premiers  u  '
356
1
    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
6
    DO j = 2, jfiltnu
366
367
165
       DO i=1,iim
368
160
          coff = coefilu(i,j)
369
160
          IF( i.LT.modfrstu(j) ) coff = 0.
370
5285
          DO k=1,iim
371
5280
             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
166
       DO k = 1, iim
382
5285
          DO i = 1, iim
383
5120
             matriceun(i,k,j) = 0.0
384
169120
             DO ii = 1, iim
385
                matriceun(i,k,j) = matriceun(i,k,j) &
386
168960
                     + 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
6
    DO j = jfiltsu, jjm
396
397
165
       DO i=1,iim
398
160
          coff = coefilu(i,j)
399
160
          IF( i.LT.modfrstu(j) ) coff = 0.
400
5285
          DO k=1,iim
401
5280
             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
166
       DO k = 1, iim
413
5285
          DO i = 1, iim
414
5120
             matriceus(i,k,j-jfiltsu+1) = 0.0
415
169120
             DO ii = 1, iim
416
                matriceus(i,k,j-jfiltsu+1) = matriceus(i,k,j-jfiltsu+1) &
417
168960
                     + 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
6
    DO j = 1, jfiltnv
433
434
165
       DO i = 1, iim
435
160
          coff = coefilv(i,j)
436
160
          IF( i.LT.modfrstv(j) ) coff = 0.
437
5285
          DO k = 1, iim
438
5280
             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
166
       DO k = 1, iim
449
5285
          DO i = 1, iim
450
5120
             matricevn(i,k,j) = 0.0
451
169120
             DO ii = 1, iim
452
                matricevn(i,k,j) = matricevn(i,k,j) &
453
168960
                     + 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
6
    DO j = jfiltsv, jjm
463
464
165
       DO i = 1, iim
465
160
          coff = coefilv(i,j)
466
160
          IF( i.LT.modfrstv(j) ) coff = 0.
467
5285
          DO k = 1, iim
468
5280
             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
166
       DO k = 1, iim
480
5285
          DO i = 1, iim
481
5120
             matricevs(i,k,j-jfiltsv+1) = 0.0
482
169120
             DO ii = 1, iim
483
                matricevs(i,k,j-jfiltsv+1) = matricevs(i,k,j-jfiltsv+1) &
484
168960
                     + 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
6
    DO j = 2, jfiltnu
500
501
165
       DO i = 1,iim
502
160
          coff = coefilu(i,j)/ ( 1. + coefilu(i,j) )
503
160
          IF( i.LT.modfrstu(j) ) coff = 0.
504
5285
          DO k=1,iim
505
5280
             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
166
       DO k = 1, iim
516
5285
          DO i = 1, iim
517
5120
             matrinvn(i,k,j) = 0.0
518
169120
             DO ii = 1, iim
519
                matrinvn(i,k,j) = matrinvn(i,k,j) &
520
168960
                     + 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
6
    DO j = jfiltsu, jjm
530
531
165
       DO i = 1,iim
532
160
          coff = coefilu(i,j) / ( 1. + coefilu(i,j) )
533
160
          IF( i.LT.modfrstu(j) ) coff = 0.
534
5285
          DO k=1,iim
535
5280
             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
166
       DO k = 1, iim
546
5285
          DO i = 1, iim
547
5120
             matrinvs(i,k,j-jfiltsu+1) = 0.0
548
169120
             DO ii = 1, iim
549
                matrinvs(i,k,j-jfiltsu+1) = matrinvs(i,k,j-jfiltsu+1) &
550
168960
                     + 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
    !
570
334 FORMAT(1x,24i3)
571
755 FORMAT(1x,6f10.3,i3)
572
573
1
    RETURN
574
  END SUBROUTINE inifilr
575
576
END MODULE filtreg_mod