source: LMDZ6/branches/Amaury_dev/libf/filtrez/lmdz_filtreg.F90 @ 5133

Last change on this file since 5133 was 5123, checked in by abarral, 5 months ago

Correct various minor mistakes from previous commits

  • 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: 28.4 KB
Line 
1! $Id$
2
3MODULE lmdz_filtreg
4  IMPLICIT NONE; PRIVATE
5  PUBLIC matriceun, matriceus, matricevn, matricevs, matrinvn, matrinvs, &
6          inifilr, filtreg
7
8  REAL, DIMENSION(:, :, :), ALLOCATABLE :: matriceun, matriceus, matricevn
9  REAL, DIMENSION(:, :, :), ALLOCATABLE :: matricevs, matrinvn, matrinvs
10
11CONTAINS
12
13  SUBROUTINE filtreg(champ, nlat, nbniv, ifiltre, iaire, &
14          griscal, iter)
15    USE lmdz_coefils, ONLY: jfiltnu, jfiltnv, jfiltsu, jfiltsv, sddu, sddv, unsddu, unsddv, modfrstv, modfrstu
16
17    !=======================================================================
18    !
19    !   Auteur: P. Le Van        07/10/97
20    !   ------
21    !
22    !   Objet: filtre matriciel longitudinal ,avec les matrices precalculees
23    !                 pour l'operateur  Filtre    .
24    !   ------
25    !
26    !   Arguments:
27    !   ----------
28    !
29    !  nblat                 nombre de latitudes a filtrer
30    !  nbniv                 nombre de niveaux verticaux a filtrer
31    !  champ(iip1,nblat,nbniv)  en entree : champ a filtrer
32    !                        en sortie : champ filtre
33    !  ifiltre               +1  Transformee directe
34    !                        -1  Transformee inverse
35    !                        +2  Filtre directe
36    !                        -2  Filtre inverse
37    !
38    !  iaire                 1   si champ intensif
39    !                        2   si champ extensif (pondere par les aires)
40    !
41    !  iter                  1   filtre simple
42    !
43    !=======================================================================
44    !
45    !
46    !                  Variable Intensive
47    !            ifiltre = 1     filtre directe
48    !            ifiltre =-1     filtre inverse
49    !
50    !                  Variable Extensive
51    !            ifiltre = 2     filtre directe
52    !            ifiltre =-2     filtre inverse
53    !
54    !
55    INCLUDE "dimensions.h"
56    INCLUDE "paramet.h"
57
58    INTEGER :: nlat, nbniv, ifiltre, iter
59    INTEGER :: i, j, l, k
60    INTEGER :: iim2, immjm
61    INTEGER :: jdfil1, jdfil2, jffil1, jffil2, jdfil, jffil
62
63    REAL :: champ(iip1, nlat, nbniv)
64
65    REAL :: eignq(iim, nlat, nbniv), sdd1(iim), sdd2(iim)
66    LOGICAL :: griscal
67    INTEGER :: hemisph, iaire
68
69    LOGICAL, SAVE :: first = .TRUE.
70
71    REAL, SAVE :: sdd12(iim, 4)
72
73    INTEGER, PARAMETER :: type_sddu = 1
74    INTEGER, PARAMETER :: type_sddv = 2
75    INTEGER, PARAMETER :: type_unsddu = 3
76    INTEGER, PARAMETER :: type_unsddv = 4
77
78    INTEGER :: sdd1_type, sdd2_type
79
80    IF (iim == 1) return ! no filtre in 2D y-z
81
82    IF (first) THEN
83      sdd12(1:iim, type_sddu) = sddu(1:iim)
84      sdd12(1:iim, type_sddv) = sddv(1:iim)
85      sdd12(1:iim, type_unsddu) = unsddu(1:iim)
86      sdd12(1:iim, type_unsddv) = unsddv(1:iim)
87
88      first = .FALSE.
89    ENDIF
90
91    IF(ifiltre==1.OR.ifiltre==-1) &
92            stop 'Pas de transformee simple dans cette version'
93
94    IF(iter== 2)  THEN
95      PRINT *, ' Pas d iteration du filtre dans cette version !'&
96              &, ' Utiliser old_filtreg et repasser !'
97      STOP
98    ENDIF
99
100    IF(ifiltre== -2 .AND..NOT.griscal)     THEN
101      PRINT *, ' Cette routine ne calcule le filtre inverse que ' &
102              , ' sur la grille des scalaires !'
103      STOP
104    ENDIF
105
106    IF(ifiltre/=2 .AND.ifiltre/= - 2)  THEN
107      PRINT *, ' Probleme dans filtreg car ifiltre NE 2 et NE -2' &
108              , ' corriger et repasser !'
109      STOP
110    ENDIF
111
112    iim2 = iim * iim
113    immjm = iim * jjm
114
115    IF(griscal)   THEN
116      IF(nlat /= jjp1)  THEN
117        PRINT  1111
118        STOP
119      ELSE
120
121        IF(iaire==1)  THEN
122          sdd1_type = type_sddv
123          sdd2_type = type_unsddv
124        ELSE
125          sdd1_type = type_unsddv
126          sdd2_type = type_sddv
127        ENDIF
128
129        ! IF( iaire.EQ.1 )  THEN
130        !    CALL SCOPY(  iim,    sddv, 1,  sdd1, 1 )
131        !    CALL SCOPY(  iim,  unsddv, 1,  sdd2, 1 )
132        ! ELSE
133        !    CALL SCOPY(  iim,  unsddv, 1,  sdd1, 1 )
134        !    CALL SCOPY(  iim,    sddv, 1,  sdd2, 1 )
135        ! END IF
136
137        jdfil1 = 2
138        jffil1 = jfiltnu
139        jdfil2 = jfiltsu
140        jffil2 = jjm
141      END IF
142    ELSE
143      IF(nlat/=jjm)  THEN
144        PRINT  2222
145        STOP
146      ELSE
147
148        IF(iaire==1)  THEN
149          sdd1_type = type_sddu
150          sdd2_type = type_unsddu
151        ELSE
152          sdd1_type = type_unsddu
153          sdd2_type = type_sddu
154        ENDIF
155
156        ! IF( iaire.EQ.1 )  THEN
157        !    CALL SCOPY(  iim,    sddu, 1,  sdd1, 1 )
158        !    CALL SCOPY(  iim,  unsddu, 1,  sdd2, 1 )
159        ! ELSE
160        !    CALL SCOPY(  iim,  unsddu, 1,  sdd1, 1 )
161        !    CALL SCOPY(  iim,    sddu, 1,  sdd2, 1 )
162        ! END IF
163
164        jdfil1 = 1
165        jffil1 = jfiltnv
166        jdfil2 = jfiltsv
167        jffil2 = jjm
168      END IF
169    END IF
170
171    DO hemisph = 1, 2
172
173      IF (hemisph==1)  THEN
174        jdfil = jdfil1
175        jffil = jffil1
176      ELSE
177        jdfil = jdfil2
178        jffil = jffil2
179      END IF
180
181      DO l = 1, nbniv
182        DO j = jdfil, jffil
183          DO i = 1, iim
184            champ(i, j, l) = champ(i, j, l) * sdd12(i, sdd1_type) ! sdd1(i)
185          END DO
186        END DO
187      END DO
188
189      IF(hemisph == 1)      THEN
190
191        IF(ifiltre == -2)   THEN
192
193          DO j = jdfil, jffil
194#ifdef BLAS
195              CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
196                    matrinvn(1,1,j), &
197                    iim, champ(1,j,1), iip1*nlat, 0.0, &
198                    eignq(1,j-jdfil+1,1), iim*nlat)
199#else
200            eignq(:, j - jdfil + 1, :) &
201                    = matmul(matrinvn(:, :, j), champ(:iim, j, :))
202#endif
203          END DO
204
205        ELSE IF (griscal)     THEN
206
207          DO j = jdfil, jffil
208#ifdef BLAS
209              CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
210                    matriceun(1,1,j), &
211                    iim, champ(1,j,1), iip1*nlat, 0.0, &
212                    eignq(1,j-jdfil+1,1), iim*nlat)
213#else
214            eignq(:, j - jdfil + 1, :) &
215                    = matmul(matriceun(:, :, j), champ(:iim, j, :))
216#endif
217          END DO
218
219        ELSE
220
221          DO j = jdfil, jffil
222#ifdef BLAS
223              CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
224                    matricevn(1,1,j), &
225                    iim, champ(1,j,1), iip1*nlat, 0.0, &
226                    eignq(1,j-jdfil+1,1), iim*nlat)
227#else
228            eignq(:, j - jdfil + 1, :) &
229                    = matmul(matricevn(:, :, j), champ(:iim, j, :))
230#endif
231          END DO
232
233        ENDIF
234
235      ELSE
236
237        IF(ifiltre == -2)   THEN
238
239          DO j = jdfil, jffil
240#ifdef BLAS
241              CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
242                    matrinvs(1,1,j-jfiltsu+1), &
243                    iim, champ(1,j,1), iip1*nlat, 0.0, &
244                    eignq(1,j-jdfil+1,1), iim*nlat)
245#else
246            eignq(:, j - jdfil + 1, :) &
247                    = matmul(matrinvs(:, :, j - jfiltsu + 1), &
248                    champ(:iim, j, :))
249#endif
250          END DO
251
252        ELSE IF (griscal)     THEN
253
254          DO j = jdfil, jffil
255#ifdef BLAS
256              CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
257                    matriceus(1,1,j-jfiltsu+1), &
258                    iim, champ(1,j,1), iip1*nlat, 0.0, &
259                    eignq(1,j-jdfil+1,1), iim*nlat)
260#else
261            eignq(:, j - jdfil + 1, :) &
262                    = matmul(matriceus(:, :, j - jfiltsu + 1), &
263                    champ(:iim, j, :))
264#endif
265          END DO
266
267        ELSE
268
269          DO j = jdfil, jffil
270#ifdef BLAS
271              CALL SGEMM("N", "N", iim, nbniv, iim, 1.0, &
272                    matricevs(1,1,j-jfiltsv+1), &
273                    iim, champ(1,j,1), iip1*nlat, 0.0, &
274                    eignq(1,j-jdfil+1,1), iim*nlat)
275#else
276            eignq(:, j - jdfil + 1, :) &
277                    = matmul(matricevs(:, :, j - jfiltsv + 1), &
278                    champ(:iim, j, :))
279#endif
280          END DO
281
282        ENDIF
283
284      ENDIF
285
286      IF(ifiltre== 2)  THEN
287
288        DO l = 1, nbniv
289          DO j = jdfil, jffil
290            DO i = 1, iim
291              champ(i, j, l) = &
292                      (champ(i, j, l) + eignq(i, j - jdfil + 1, l)) &
293                              * sdd12(i, sdd2_type) ! sdd2(i)
294            END DO
295          END DO
296        END DO
297
298      ELSE
299
300        DO l = 1, nbniv
301          DO j = jdfil, jffil
302            DO i = 1, iim
303              champ(i, j, l) = &
304                      (champ(i, j, l) - eignq(i, j - jdfil + 1, l)) &
305                              * sdd12(i, sdd2_type) ! sdd2(i)
306            END DO
307          END DO
308        END DO
309
310      ENDIF
311
312      DO l = 1, nbniv
313        DO j = jdfil, jffil
314          champ(iip1, j, l) = champ(1, j, l)
315        END DO
316      END DO
317
318    ENDDO
319
320    1111   FORMAT(//20x, 'ERREUR dans le dimensionnement du tableau  CHAMP a&
321            &     filtrer, sur la grille des scalaires'/)
322    2222   FORMAT(//20x, 'ERREUR dans le dimensionnement du tableau CHAMP a fi&
323            &     ltrer, sur la grille de V ou de Z'/)
324    RETURN
325  END SUBROUTINE filtreg
326
327  SUBROUTINE inifgn(dv)
328    !
329    !    ...  H.Upadyaya , O.Sharma  ...
330    !
331    USE lmdz_coefils, ONLY: sddv, sddu, unsddu, unsddv, eignfnv, eignfnu
332    USE lmdz_ssum_scopy, ONLY: ssum
333    IMPLICIT NONE
334    !
335    include "dimensions.h"
336    include "paramet.h"
337    include "comgeom.h"
338    !
339    REAL :: vec(iim, iim), vec1(iim, iim)
340    REAL :: dlonu(iim), dlonv(iim)
341    REAL :: du(iim), dv(iim), d(iim)
342    REAL :: pi
343    INTEGER :: i, j, k, imm1, nrot
344    !
345
346    imm1 = iim - 1
347    pi = 2. * ASIN(1.)
348    !
349    DO i = 1, iim
350      dlonu(i) = xprimu(i)
351      dlonv(i) = xprimv(i)
352    END DO
353
354    DO i = 1, iim
355      sddv(i) = SQRT(dlonv(i))
356      sddu(i) = SQRT(dlonu(i))
357      unsddu(i) = 1. / sddu(i)
358      unsddv(i) = 1. / sddv(i)
359    END DO
360    !
361    DO j = 1, iim
362      DO i = 1, iim
363        vec(i, j) = 0.
364        vec1(i, j) = 0.
365        eignfnv(i, j) = 0.
366        eignfnu(i, j) = 0.
367      END DO
368    END DO
369    !
370    !
371    eignfnv(1, 1) = -1.
372    eignfnv(iim, 1) = 1.
373    DO i = 1, imm1
374      eignfnv(i + 1, i + 1) = -1.
375      eignfnv(i, i + 1) = 1.
376    END DO
377    DO j = 1, iim
378      DO i = 1, iim
379        eignfnv(i, j) = eignfnv(i, j) / (sddu(i) * sddv(j))
380      END DO
381    END DO
382    DO j = 1, iim
383      DO i = 1, iim
384        eignfnu(i, j) = -eignfnv(j, i)
385      END DO
386    END DO
387    !
388    DO j = 1, iim
389      DO i = 1, iim
390        vec (i, j) = 0.0
391        vec1(i, j) = 0.0
392        DO k = 1, iim
393          vec (i, j) = vec(i, j) + eignfnu(i, k) * eignfnv(k, j)
394          vec1(i, j) = vec1(i, j) + eignfnv(i, k) * eignfnu(k, j)
395        ENDDO
396      ENDDO
397    ENDDO
398
399    !
400    CALL jacobi(vec, iim, iim, dv, eignfnv, nrot)
401    CALL acc(eignfnv, d, iim)
402    CALL eigen_sort(dv, eignfnv, iim, iim)
403    !
404    CALL jacobi(vec1, iim, iim, du, eignfnu, nrot)
405    CALL acc(eignfnu, d, iim)
406    CALL eigen_sort(du, eignfnu, iim, iim)
407
408    !c   ancienne version avec appels IMSL
409    !
410    ! CALL MXM(eignfnu,iim,eignfnv,iim,vec,iim)
411    ! CALL MXM(eignfnv,iim,eignfnu,iim,vec1,iim)
412    !     CALL EVCSF(iim,vec,iim,dv,eignfnv,iim)
413    ! CALL acc(eignfnv,d,iim)
414    ! CALL eigen(eignfnv,dv)
415    !
416    ! CALL EVCSF(iim,vec1,iim,du,eignfnu,iim)
417    ! CALL acc(eignfnu,d,iim)
418    ! CALL eigen(eignfnu,du)
419
420    RETURN
421  END SUBROUTINE inifgn
422
423  SUBROUTINE JACOBI(A, N, NP, D, V, NROT)
424    IMPLICIT NONE
425    ! Arguments:
426    INTEGER, INTENT(IN) :: N
427    INTEGER, INTENT(IN) :: NP
428    INTEGER, INTENT(OUT) :: NROT
429    REAL, INTENT(INOUT) :: A(NP, NP)
430    REAL, INTENT(OUT) :: D(NP)
431    REAL, INTENT(OUT) :: V(NP, NP)
432
433    ! local variables:
434    INTEGER :: IP, IQ, I, J
435    REAL :: SM, TRESH, G, H, T, THETA, C, S, TAU
436    REAL :: B(N)
437    REAL :: Z(N)
438
439    DO IP = 1, N
440      DO IQ = 1, N
441        V(IP, IQ) = 0.
442      ENDDO
443      V(IP, IP) = 1.
444    ENDDO
445    DO IP = 1, N
446      B(IP) = A(IP, IP)
447      D(IP) = B(IP)
448      Z(IP) = 0.
449    ENDDO
450    NROT = 0
451    DO I = 1, 50 ! 50? I suspect this should be NP
452      !     but convergence is fast enough anyway
453      SM = 0.
454      DO IP = 1, N - 1
455        DO IQ = IP + 1, N
456          SM = SM + ABS(A(IP, IQ))
457        ENDDO
458      ENDDO
459      IF(SM==0.)RETURN
460      IF(I<4)THEN
461        TRESH = 0.2 * SM / N**2
462      ELSE
463        TRESH = 0.
464      ENDIF
465      DO IP = 1, N - 1
466        DO IQ = IP + 1, N
467          G = 100. * ABS(A(IP, IQ))
468          IF((I>4).AND.(ABS(D(IP)) + G==ABS(D(IP))) &
469                  .AND.(ABS(D(IQ)) + G==ABS(D(IQ))))THEN
470            A(IP, IQ) = 0.
471          ELSE IF(ABS(A(IP, IQ))>TRESH)THEN
472            H = D(IQ) - D(IP)
473            IF(ABS(H) + G==ABS(H))THEN
474              T = A(IP, IQ) / H
475            ELSE
476              THETA = 0.5 * H / A(IP, IQ)
477              T = 1. / (ABS(THETA) + SQRT(1. + THETA**2))
478              IF(THETA<0.)T = -T
479            ENDIF
480            C = 1. / SQRT(1 + T**2)
481            S = T * C
482            TAU = S / (1. + C)
483            H = T * A(IP, IQ)
484            Z(IP) = Z(IP) - H
485            Z(IQ) = Z(IQ) + H
486            D(IP) = D(IP) - H
487            D(IQ) = D(IQ) + H
488            A(IP, IQ) = 0.
489            DO J = 1, IP - 1
490              G = A(J, IP)
491              H = A(J, IQ)
492              A(J, IP) = G - S * (H + G * TAU)
493              A(J, IQ) = H + S * (G - H * TAU)
494            ENDDO
495            DO J = IP + 1, IQ - 1
496              G = A(IP, J)
497              H = A(J, IQ)
498              A(IP, J) = G - S * (H + G * TAU)
499              A(J, IQ) = H + S * (G - H * TAU)
500            ENDDO
501            DO J = IQ + 1, N
502              G = A(IP, J)
503              H = A(IQ, J)
504              A(IP, J) = G - S * (H + G * TAU)
505              A(IQ, J) = H + S * (G - H * TAU)
506            ENDDO
507            DO J = 1, N
508              G = V(J, IP)
509              H = V(J, IQ)
510              V(J, IP) = G - S * (H + G * TAU)
511              V(J, IQ) = H + S * (G - H * TAU)
512            ENDDO
513            NROT = NROT + 1
514          ENDIF
515        ENDDO
516      ENDDO
517      DO IP = 1, N
518        B(IP) = B(IP) + Z(IP)
519        D(IP) = B(IP)
520        Z(IP) = 0.
521      ENDDO
522    ENDDO ! of DO I=1,50
523    STOP 'Jacobi: 50 iterations should never happen'
524
525  END SUBROUTINE JACOBI
526
527  SUBROUTINE eigen_sort(d, v, n, np)
528    INTEGER :: n, np
529    REAL :: d(np), v(np, np)
530    INTEGER :: i, j, k
531    REAL :: p
532
533    DO i = 1, n - 1
534      k = i
535      p = d(i)
536      DO j = i + 1, n
537        IF(d(j)>=p) THEN
538          k = j
539          p = d(j)
540        ENDIF
541      ENDDO
542
543      IF(k/=i) THEN
544        d(k) = d(i)
545        d(i) = p
546        DO j = 1, n
547          p = v(j, i)
548          v(j, i) = v(j, k)
549          v(j, k) = p
550        ENDDO
551      ENDIF
552    ENDDO
553
554    RETURN
555  END SUBROUTINE eigen_sort
556
557  SUBROUTINE acc(vec, d, im)
558    USE lmdz_ssum_scopy, ONLY: ssum
559    IMPLICIT NONE
560    INTEGER :: im
561    REAL :: vec(im, im), d(im)
562    INTEGER :: i, j
563    REAL :: sum
564    do j = 1, im
565      do i = 1, im
566        d(i) = vec(i, j) * vec(i, j)
567      enddo
568      sum = ssum(im, d, 1)
569      sum = sqrt(sum)
570      do i = 1, im
571        vec(i, j) = vec(i, j) / sum
572      enddo
573    enddo
574    RETURN
575  END SUBROUTINE acc
576
577
578  SUBROUTINE inifilr
579#ifdef CPP_PARA
580  USE lmdz_filtre_fft, ONLY: use_filtre_fft,Init_filtre_fft
581  USE lmdz_filtre_fft_loc, ONLY: Init_filtre_fft_loc=>Init_filtre_fft    !
582#endif
583    USE serre_mod, ONLY: alphax
584    USE logic_mod, ONLY: fxyhypb, ysinus
585    USE comconst_mod, ONLY: maxlatfilter
586    USE lmdz_coefils, ONLY: modfrstv, modfrstu, jfiltnu, jfiltnv, coefilu, coefilv, &
587            coefilu2, coefilv2, eignfnv, eignfnu, jfiltsu, jfiltsv
588
589    !    ... H. Upfiltreg_modadhyaya, O.Sharma   ...
590
591    !     version 3 .....
592
593    !     Correction  le 28/10/97    P. Le Van .
594    !  -------------------------------------------------------------------
595    include "dimensions.h"
596    include "paramet.h"
597    include "comgeom.h"
598
599    REAL  dlonu(iim), dlatu(jjm)
600    REAL  rlamda(iim), eignvl(iim)
601
602    REAL    lamdamax, pi, cof
603    INTEGER i, j, modemax, imx, k, kf, ii
604    REAL dymin, dxmin, colat0
605    REAL eignft(iim, iim), coff
606
607    LOGICAL, SAVE :: first_call_inifilr = .TRUE.
608
609    INTEGER   ISMIN
610    EXTERNAL  ISMIN
611    INTEGER iymin
612    INTEGER ixmineq
613
614    ! ------------------------------------------------------------
615    !   This routine computes the eigenfunctions of the laplacien
616    !   on the stretched grid, and the filtering coefficients
617
618    !  We designate:
619    !   eignfn   eigenfunctions of the discrete laplacien
620    !   eigenvl  eigenvalues
621    !   jfiltn   indexof the last scalar line filtered in NH
622    !   jfilts   index of the first line filtered in SH
623    !   modfrst  index of the mode from WHERE modes are filtered
624    !   modemax  maximum number of modes ( im )
625    !   coefil   filtering coefficients ( lamda_max*COS(rlat)/lamda )
626    !   sdd      SQRT( dx )
627
628    !     the modes are filtered from modfrst to modemax
629
630    !-----------------------------------------------------------
631
632    IF (iim == 1) return ! No filtre in 2D y-z
633
634    pi = 2. * ASIN(1.)
635
636    DO i = 1, iim
637      dlonu(i) = xprimu(i)
638    ENDDO
639
640    CALL inifgn(eignvl)
641
642    PRINT *, 'inifilr: EIGNVL '
643    PRINT 250, eignvl
644    250 FORMAT(1x, 5e14.6)
645
646    ! compute eigenvalues and eigenfunctions
647
648
649    !.................................................................
650
651    !  compute the filtering coefficients for scalar lines and
652    !  meridional wind v-lines
653
654    !  we filter all those latitude lines WHERE coefil < 1
655    !  NO FILTERING AT POLES
656
657    !  colat0 is to be used  when alpha (stretching coefficient)
658    !  is set equal to zero for the regular grid CASE
659
660    !    .......   Calcul  de  colat0   .........
661    !     .....  colat0 = minimum de ( 0.5, min dy/ min dx )   ...
662
663    DO j = 1, jjm
664      dlatu(j) = rlatu(j) - rlatu(j + 1)
665    ENDDO
666
667    dxmin = dlonu(1)
668    DO  i = 2, iim
669      dxmin = MIN(dxmin, dlonu(i))
670    ENDDO
671    dymin = dlatu(1)
672    DO j = 2, jjm
673      dymin = MIN(dymin, dlatu(j))
674    ENDDO
675
676    ! For a regular grid, we want the filter to start at latitudes
677    ! corresponding to lengths dx of the same size as dy (in terms
678    ! of angles: dx=2*dy) => at colat0=0.5 (i.e. colatitude=30 degrees
679    !  <=> latitude=60 degrees).
680    ! Same idea for the zoomed grid: start filtering polewards as soon
681    ! as length dx becomes of the same size as dy
682
683    ! if maxlatfilter >0, prescribe the colat0 value from the .def files
684
685    IF (maxlatfilter < 0.) THEN
686
687      colat0 = MIN(0.5, dymin / dxmin)
688      ! colat0  =  1.
689
690      IF(.NOT.fxyhypb.AND.ysinus)  THEN
691        colat0 = 0.6
692        !         ...... a revoir  pour  ysinus !   .......
693        alphax = 0.
694      ENDIF
695
696    ELSE
697
698      colat0 = (90.0 - maxlatfilter) / 180.0 * pi
699
700    ENDIF
701
702    PRINT 50, colat0, alphax
703    50  FORMAT(/15x, ' Inifilr colat0 alphax ', 2e16.7)
704
705    IF(alphax==1.)  THEN
706      PRINT *, ' Inifilr  alphax doit etre  <  a 1.  Corriger '
707      STOP
708    ENDIF
709
710    lamdamax = iim / (pi * colat0 * (1. - alphax))
711
712    !                        ... Correction  le 28/10/97  ( P.Le Van ) ..
713
714    DO i = 2, iim
715      rlamda(i) = lamdamax / SQRT(ABS(eignvl(i)))
716    ENDDO
717
718    DO j = 1, jjm
719      DO i = 1, iim
720        coefilu(i, j) = 0.0
721        coefilv(i, j) = 0.0
722        coefilu2(i, j) = 0.0
723        coefilv2(i, j) = 0.0
724      ENDDO
725    ENDDO
726
727    !    ... Determination de jfiltnu,jfiltnv,jfiltsu,jfiltsv ....
728    !    .........................................................
729
730    modemax = iim
731
732    !!!!    imx = modemax - 4 * (modemax/iim)
733
734    imx = iim
735
736    PRINT *, 'inifilr: TRUNCATION AT ', imx
737
738    ! Ehouarn: set up some defaults
739    jfiltnu = 2 ! avoid north pole
740    jfiltsu = jjm ! avoid south pole (which is at jjm+1)
741    jfiltnv = 1 ! NB: no poles on the V grid
742    jfiltsv = jjm
743
744    DO j = 2, jjm / 2 + 1
745      cof = COS(rlatu(j)) / colat0
746      IF (cof < 1.) THEN
747        IF(rlamda(imx) * COS(rlatu(j))<1.) THEN
748          jfiltnu = j
749        ENDIF
750      ENDIF
751
752      cof = COS(rlatu(jjp1 - j + 1)) / colat0
753      IF (cof < 1.) THEN
754        IF(rlamda(imx) * COS(rlatu(jjp1 - j + 1))<1.) THEN
755          jfiltsu = jjp1 - j + 1
756        ENDIF
757      ENDIF
758    ENDDO
759
760    DO j = 1, jjm / 2
761      cof = COS(rlatv(j)) / colat0
762      IF (cof < 1.) THEN
763        IF(rlamda(imx) * COS(rlatv(j))<1.) THEN
764          jfiltnv = j
765        ENDIF
766      ENDIF
767
768      cof = COS(rlatv(jjm - j + 1)) / colat0
769      IF (cof < 1.) THEN
770        IF(rlamda(imx) * COS(rlatv(jjm - j + 1))<1.) THEN
771          jfiltsv = jjm - j + 1
772        ENDIF
773      ENDIF
774    ENDDO
775
776    IF(jfiltnu> jjm / 2 + 1)  THEN
777      PRINT *, ' jfiltnu en dehors des valeurs acceptables ', jfiltnu
778      STOP
779    ENDIF
780
781    IF(jfiltsu>  jjm + 1)  THEN
782      PRINT *, ' jfiltsu en dehors des valeurs acceptables ', jfiltsu
783      STOP
784    ENDIF
785
786    IF(jfiltnv> jjm / 2)  THEN
787      PRINT *, ' jfiltnv en dehors des valeurs acceptables ', jfiltnv
788      STOP
789    ENDIF
790
791    IF(jfiltsv>     jjm)  THEN
792      PRINT *, ' jfiltsv en dehors des valeurs acceptables ', jfiltsv
793      STOP
794    ENDIF
795
796    PRINT *, 'inifilr: jfiltnv jfiltsv jfiltnu jfiltsu ', &
797            jfiltnv, jfiltsv, jfiltnu, jfiltsu
798
799    IF(first_call_inifilr) THEN
800      ALLOCATE(matriceun(iim, iim, jfiltnu))
801      ALLOCATE(matriceus(iim, iim, jjm - jfiltsu + 1))
802      ALLOCATE(matricevn(iim, iim, jfiltnv))
803      ALLOCATE(matricevs(iim, iim, jjm - jfiltsv + 1))
804      ALLOCATE(matrinvn(iim, iim, jfiltnu))
805      ALLOCATE(matrinvs(iim, iim, jjm - jfiltsu + 1))
806      first_call_inifilr = .FALSE.
807    ENDIF
808
809    !   ... Determination de coefilu,coefilv,n=modfrstu,modfrstv ....
810    !................................................................
811
812    DO j = 1, jjm
813      !default initialization: all modes are retained (i.e. no filtering)
814      modfrstu(j) = iim
815      modfrstv(j) = iim
816    ENDDO
817
818    DO j = 2, jfiltnu
819      DO k = 2, modemax
820        cof = rlamda(k) * COS(rlatu(j))
821        IF (cof < 1.) GOTO 82
822      ENDDO
823      GOTO 84
824      82     modfrstu(j) = k
825
826      kf = modfrstu(j)
827      DO k = kf, modemax
828        cof = rlamda(k) * COS(rlatu(j))
829        coefilu(k, j) = cof - 1.
830        coefilu2(k, j) = cof * cof - 1.
831      ENDDO
832      84     CONTINUE
833    ENDDO
834
835    DO j = 1, jfiltnv
836
837      DO k = 2, modemax
838        cof = rlamda(k) * COS(rlatv(j))
839        IF (cof < 1.) GOTO 87
840      ENDDO
841      GOTO 89
842      87     modfrstv(j) = k
843
844      kf = modfrstv(j)
845      DO k = kf, modemax
846        cof = rlamda(k) * COS(rlatv(j))
847        coefilv(k, j) = cof - 1.
848        coefilv2(k, j) = cof * cof - 1.
849      ENDDO
850      89     CONTINUE
851    ENDDO
852
853    DO j = jfiltsu, jjm
854      DO k = 2, modemax
855        cof = rlamda(k) * COS(rlatu(j))
856        IF (cof < 1.) GOTO 92
857      ENDDO
858      GOTO 94
859      92     modfrstu(j) = k
860
861      kf = modfrstu(j)
862      DO k = kf, modemax
863        cof = rlamda(k) * COS(rlatu(j))
864        coefilu(k, j) = cof - 1.
865        coefilu2(k, j) = cof * cof - 1.
866      ENDDO
867      94     CONTINUE
868    ENDDO
869
870    DO j = jfiltsv, jjm
871      DO k = 2, modemax
872        cof = rlamda(k) * COS(rlatv(j))
873        IF (cof < 1.) GOTO 97
874      ENDDO
875      GOTO 99
876      97     modfrstv(j) = k
877
878      kf = modfrstv(j)
879      DO k = kf, modemax
880        cof = rlamda(k) * COS(rlatv(j))
881        coefilv(k, j) = cof - 1.
882        coefilv2(k, j) = cof * cof - 1.
883      ENDDO
884      99     CONTINUE
885    ENDDO
886
887    IF(jfiltnv>=jjm / 2 .OR. jfiltnu>=jjm / 2)THEN
888      ! Ehouarn: and what are these for??? Trying to handle a limit case
889      !          where filters extend to and meet at the equator?
890      IF(jfiltnv==jfiltsv)jfiltsv = 1 + jfiltnv
891      IF(jfiltnu==jfiltsu)jfiltsu = 1 + jfiltnu
892
893      PRINT *, 'jfiltnv jfiltsv jfiltnu jfiltsu', &
894              jfiltnv, jfiltsv, jfiltnu, jfiltsu
895    ENDIF
896
897    PRINT *, '   Modes premiers  v  '
898    PRINT 334, modfrstv
899    PRINT *, '   Modes premiers  u  '
900    PRINT 334, modfrstu
901
902    !   ...................................................................
903
904    !   ... Calcul de la matrice filtre 'matriceu'  pour les champs situes
905    !                       sur la grille scalaire                 ........
906    !   ...................................................................
907
908    DO j = 2, jfiltnu
909
910      DO i = 1, iim
911        coff = coefilu(i, j)
912        IF(i<modfrstu(j)) coff = 0.
913        DO k = 1, iim
914          eignft(i, k) = eignfnv(k, i) * coff
915        ENDDO
916      ENDDO ! of DO i=1,iim
917
918#ifdef BLAS
919       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
920            eignfnv, iim, eignft, iim, 0.0, matriceun(1,1,j), iim)
921#else
922      DO k = 1, iim
923        DO i = 1, iim
924          matriceun(i, k, j) = 0.0
925          DO ii = 1, iim
926            matriceun(i, k, j) = matriceun(i, k, j) &
927                    + eignfnv(i, ii) * eignft(ii, k)
928          ENDDO
929        ENDDO
930      ENDDO ! of DO k = 1, iim
931#endif
932
933    ENDDO ! of DO j = 2, jfiltnu
934
935    DO j = jfiltsu, jjm
936
937      DO i = 1, iim
938        coff = coefilu(i, j)
939        IF(i<modfrstu(j)) coff = 0.
940        DO k = 1, iim
941          eignft(i, k) = eignfnv(k, i) * coff
942        ENDDO
943      ENDDO ! of DO i=1,iim
944#ifdef BLAS
945       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
946            eignfnv, iim, eignft, iim, 0.0, &
947            matriceus(1,1,j-jfiltsu+1), iim)
948#else
949      DO k = 1, iim
950        DO i = 1, iim
951          matriceus(i, k, j - jfiltsu + 1) = 0.0
952          DO ii = 1, iim
953            matriceus(i, k, j - jfiltsu + 1) = matriceus(i, k, j - jfiltsu + 1) &
954                    + eignfnv(i, ii) * eignft(ii, k)
955          ENDDO
956        ENDDO
957      ENDDO ! of DO k = 1, iim
958#endif
959
960    ENDDO ! of DO j = jfiltsu, jjm
961
962    !   ...................................................................
963
964    !   ... Calcul de la matrice filtre 'matricev'  pour les champs situes
965    !                       sur la grille   de V ou de Z           ........
966    !   ...................................................................
967
968    DO j = 1, jfiltnv
969
970      DO i = 1, iim
971        coff = coefilv(i, j)
972        IF(i<modfrstv(j)) coff = 0.
973        DO k = 1, iim
974          eignft(i, k) = eignfnu(k, i) * coff
975        ENDDO
976      ENDDO
977
978#ifdef BLAS
979       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
980            eignfnu, iim, eignft, iim, 0.0, matricevn(1,1,j), iim)
981#else
982      DO k = 1, iim
983        DO i = 1, iim
984          matricevn(i, k, j) = 0.0
985          DO ii = 1, iim
986            matricevn(i, k, j) = matricevn(i, k, j) &
987                    + eignfnu(i, ii) * eignft(ii, k)
988          ENDDO
989        ENDDO
990      ENDDO
991#endif
992
993    ENDDO ! of DO j = 1, jfiltnv
994
995    DO j = jfiltsv, jjm
996
997      DO i = 1, iim
998        coff = coefilv(i, j)
999        IF(i<modfrstv(j)) coff = 0.
1000        DO k = 1, iim
1001          eignft(i, k) = eignfnu(k, i) * coff
1002        ENDDO
1003      ENDDO
1004
1005#ifdef BLAS
1006       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
1007            eignfnu, iim, eignft, iim, 0.0, &
1008            matricevs(1,1,j-jfiltsv+1), iim)
1009#else
1010      DO k = 1, iim
1011        DO i = 1, iim
1012          matricevs(i, k, j - jfiltsv + 1) = 0.0
1013          DO ii = 1, iim
1014            matricevs(i, k, j - jfiltsv + 1) = matricevs(i, k, j - jfiltsv + 1) &
1015                    + eignfnu(i, ii) * eignft(ii, k)
1016          ENDDO
1017        ENDDO
1018      ENDDO
1019#endif
1020
1021    ENDDO ! of DO j = jfiltsv, jjm
1022
1023    !   ...................................................................
1024
1025    !   ... Calcul de la matrice filtre 'matrinv'  pour les champs situes
1026    !              sur la grille scalaire , pour le filtre inverse ........
1027    !   ...................................................................
1028
1029    DO j = 2, jfiltnu
1030
1031      DO i = 1, iim
1032        coff = coefilu(i, j) / (1. + coefilu(i, j))
1033        IF(i<modfrstu(j)) coff = 0.
1034        DO k = 1, iim
1035          eignft(i, k) = eignfnv(k, i) * coff
1036        ENDDO
1037      ENDDO
1038
1039#ifdef BLAS
1040       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
1041            eignfnv, iim, eignft, iim, 0.0, matrinvn(1,1,j), iim)
1042#else
1043      DO k = 1, iim
1044        DO i = 1, iim
1045          matrinvn(i, k, j) = 0.0
1046          DO ii = 1, iim
1047            matrinvn(i, k, j) = matrinvn(i, k, j) &
1048                    + eignfnv(i, ii) * eignft(ii, k)
1049          ENDDO
1050        ENDDO
1051      ENDDO
1052#endif
1053
1054    ENDDO ! of DO j = 2, jfiltnu
1055
1056    DO j = jfiltsu, jjm
1057
1058      DO i = 1, iim
1059        coff = coefilu(i, j) / (1. + coefilu(i, j))
1060        IF(i<modfrstu(j)) coff = 0.
1061        DO k = 1, iim
1062          eignft(i, k) = eignfnv(k, i) * coff
1063        ENDDO
1064      ENDDO
1065#ifdef BLAS
1066       CALL SGEMM ('N', 'N', iim, iim, iim, 1.0, &
1067            eignfnv, iim, eignft, iim, 0.0, matrinvs(1,1,j-jfiltsu+1), iim)
1068#else
1069      DO k = 1, iim
1070        DO i = 1, iim
1071          matrinvs(i, k, j - jfiltsu + 1) = 0.0
1072          DO ii = 1, iim
1073            matrinvs(i, k, j - jfiltsu + 1) = matrinvs(i, k, j - jfiltsu + 1) &
1074                    + eignfnv(i, ii) * eignft(ii, k)
1075          ENDDO
1076        ENDDO
1077      ENDDO
1078#endif
1079
1080    ENDDO ! of DO j = jfiltsu, jjm
1081
1082#ifdef CPP_PARA
1083    IF (use_filtre_fft) THEN
1084       CALL Init_filtre_fft(coefilu,modfrstu,jfiltnu,jfiltsu,  &
1085                           coefilv,modfrstv,jfiltnv,jfiltsv)
1086       CALL Init_filtre_fft_loc(coefilu,modfrstu,jfiltnu,jfiltsu,  &
1087                           coefilv,modfrstv,jfiltnv,jfiltsv)
1088    ENDIF
1089#endif
1090    !   ...................................................................
1091
1092    334 FORMAT(1x, 24i3)
1093
1094  END SUBROUTINE inifilr
1095
1096END MODULE lmdz_filtreg
Note: See TracBrowser for help on using the repository browser.