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

Last change on this file since 5278 was 5272, checked in by abarral, 2 days ago

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