source: flexpart.git/preproc/src/phgrreal.f90 @ 5763793

10.4.1_peseiFPv9.3.1FPv9.3.1b_testingFPv9.3.2GFS_025bugfixes+enhancementsdevfp9.3.1-20161214-nc4grib2nc4_repairrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 5763793 was 5763793, checked in by Anne Fouilloux <annefou@…>, 9 years ago

Add pre-processing programs and examples. This is still under development!

  • Property mode set to 100644
File size: 14.8 KB
Line 
1MODULE PHTOGR
2
3  INTEGER, PARAMETER :: MAXAUF=36000
4
5CONTAINS
6
7  SUBROUTINE PHGR213(CXMN,FELD,WSAVE,IFAX,Z,MLAT,MNAUF, &
8    &MAXL,MAXB,MLEVEL)
9
10    !     DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
11    !     FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
12    !     RAUM AUF DAS REDUZIERTE GAUSS'SCHE GITTER DURCH
13    !
14    !     CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
15    !              CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
16    !     FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
17    !   WSAVE  = Working Array fuer Fouriertransformation
18    !     Z      = LEGENDREFUNKTIONSWERTE
19    !
20    !     MNAUF    ANZAHL DER FOURIERKOEFFIZIENTEN
21    !     MAXL     ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
22    !     MAXB     ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
23    !     MLEVEL   ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
24    !
25    IMPLICIT NONE
26
27    !                   Anzahl der Gitterpunkte auf jedem Breitenkreis
28    INTEGER MLAT(MAXB/2)
29    INTEGER K,MAXL,MAXB,MLEVEL,MNAUF
30    INTEGER IND(MAXB)
31
32
33    !    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
34    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB/2)
35
36    REAL, INTENT(IN) ::  CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
37    REAL FELD(MAXL,MLEVEL)
38    REAL WSAVE(8*MAXB+15,MAXB/2)
39    INTEGER :: IFAX(10,MAXB)
40
41    IND(1)=0
42    DO  K=2,MAXB/2
43       IND(K)=IND(K-1)+MLAT(K-1)
44    ENDDO
45
46    !$OMP PARALLEL DO SCHEDULE(DYNAMIC)
47    DO  K=1,MAXB/2
48       CALL PHSYM(K,IND,CXMN,FELD,Z,WSAVE,IFAX,MLAT, &
49       &MNAUF,MAXL,MAXB,MLEVEL)
50
51    ENDDO
52    !$OMP END PARALLEL DO
53
54  END SUBROUTINE PHGR213
55  !
56  !
57  SUBROUTINE PHSYM(K,IND,CXMN,FELD,Z,WSAVE,IFAX,MLAT, &
58       &MNAUF,MAXL,MAXB,MLEVEL)
59
60    IMPLICIT NONE
61
62    INTEGER MLAT(MAXB/2)
63    INTEGER K,L,I,J,LLS,LLPS,LL,LLP,MAXL,MAXB,MLEVEL,MNAUF
64    INTEGER IND(MAXB)
65    INTEGER :: IFAX(10,MAXB)
66
67
68    !    FELD DER FOURIERKOEFFIZIENTEN
69    REAL :: CXMS(0:MAXAUF-1),CXMA(0:MAXAUF-1)
70
71    !    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
72    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB/2)
73    REAL ACR,ACI,SCR,SCI
74
75    REAL, INTENT(IN) ::  CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
76    REAL FELD(MAXL,MLEVEL)
77    REAL WSAVE(8*MAXB+15,MAXB/2)
78
79    DO  L=1,MLEVEL
80       LL=0
81       LLP=0
82       DO  I=0,MNAUF
83          SCR=0.D0
84          SCI=0.D0
85          ACR=0.D0
86          ACI=0.D0
87          LLS=LL
88          LLPS=LLP
89          IF(2*I+1.LT.MLAT(K)) THEN
90             !  Innerste Schleife aufgespalten um if-Abfrage zu sparen
91             DO J=I,MNAUF,2
92                SCR=SCR+Z(LLP,K)*CXMN(2*LL,L)
93                SCI=SCI+Z(LLP,K)*CXMN(2*LL+1,L)
94                LL=LL+2           
95                LLP=LLP+2           
96             ENDDO
97             LL=LLS+1
98             LLP=LLPS+1
99             DO  J=I+1,MNAUF,2
100                ACR=ACR+Z(LLP,K)*CXMN(2*LL,L)
101                ACI=ACI+Z(LLP,K)*CXMN(2*LL+1,L)
102                LL=LL+2           
103                LLP=LLP+2           
104             ENDDO
105          ENDIF
106          LL=LLS+(MNAUF-I+1)
107          LLP=LLPS+(MNAUF-I+3)
108          CXMS(2*I)=SCR+ACR
109          CXMS(2*I+1)=SCI+ACI
110          CXMA(2*I)=SCR-ACR
111          CXMA(2*I+1)=SCI-ACI
112       ENDDO
113       CALL RFOURTR(CXMS,WSAVE(:,K),IFAX(:,K),MNAUF,&
114            &MLAT(K))
115       FELD(IND(k)+1:IND(K)+MLAT(K),L)=CXMS(0:MLAT(K)-1)
116       CALL RFOURTR(CXMA,&
117            &WSAVE(:,K),IFAX(:,K),MNAUF,MLAT(K))
118       FELD(MAXL-IND(k)-MLAT(K)+1:MAXL-IND(k),L)=CXMA(0:MLAT(K)-1)
119
120    ENDDO
121
122  END SUBROUTINE PHSYM
123
124  SUBROUTINE PHGCUT(CXMN,FELD,WSAVE,IFAX,Z,&
125       &                  MNAUF,MMAX,MAUF,MANF,MAXL,MAXB,MLEVEL)
126
127    !     DIE ROUTINE FUEHRT EINE TRANSFORMATION EINER
128    !     FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
129    !     RAUM AUF KUGELKOORDINATEN DURCH. Es kann ein Teilausschnitt
130    !                   Der Erde angegeben werden. Diese Routine ist langsamer als
131    !                   phgrph
132    !
133    !     CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
134    !              CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
135    !     FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
136    !     Z      = SINUS DER GEOGRAFISCHEN BREITEN
137    !
138    !     MNAUF    ANZAHL DER FOURIERKOEFFIZIENTEN
139    !     MAUF     ANZAHL DER LAENGEN UND DER FOURIERKOEFFIZIENTEN
140    !     MANF     ANFANG DES LAENGENBEREICHS FUER DAS GITTER,
141    !              AUF DAS INTERPOLIERT WERDEN SOLL
142    !     MAXL     ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
143    !     MAXB     ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
144    !     MLEVEL   ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
145    !
146    IMPLICIT NONE
147    INTEGER, INTENT(IN)  :: MNAUF, MMAX, MAUF, MANF, MAXL, MAXB, MLEVEL
148    REAL, INTENT(OUT)    :: FELD(MAXL,MAXB,MLEVEL)
149    INTEGER, INTENT(INOUT) :: IFAX(10)
150
151    !    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
152    REAL, INTENT(IN)     :: Z(0:((MMAX+3)*(MMAX+4))/2,MAXB)
153    REAL, INTENT(IN)     :: CXMN(0:(MMAX+1)*(MMAX+2)-1,MLEVEL)
154    REAL, INTENT(INOUT)  :: WSAVE(4*MAUF+15)  ! work array previously
155                                              ! initialized with set99
156                                              ! (will be used in fft99)
157    INTEGER              :: J
158    LOGICAL              :: SYM
159
160
161    IF(MAUF.LE.MNAUF) WRITE(*,*) 'TOO COARSE LONGITUDE RESOLUTION'
162    IF((MANF.LT.1).OR.(MAXL.LT.1).OR. &
163         &       (MANF.GT.MAUF).OR.(MAXL.GT.MAUF)) THEN
164       WRITE(*,*) 'WRONG LONGITUDE RANGE',MANF,MAXL
165       STOP
166    ENDIF
167
168    ! Pruefe, ob Ausgabegitter symmetrisch zum Aequator ist
169    ! Wenn ja soll Symmetrie der Legendrepolynome ausgenutzt werden
170    ! Check if the output grid is symmetrical to the equator
171    ! If yes we exploit it
172
173    IF(MAXB .GT. 4) THEN
174       SYM=.TRUE.
175       DO  J=5,5
176          IF(ABS(ABS(Z(100,J))-ABS(Z(100,MAXB+1-J))).GT.1E-11) THEN
177             SYM=.FALSE.
178          ENDIF
179       ENDDO
180       WRITE(*,*) 'Symmetrisch: ',SYM
181    ELSE
182       SYM=.FALSE.
183    ENDIF
184
185
186    IF(SYM) THEN
187       !$OMP PARALLEL DO
188       DO J=1,(MAXB+1)/2
189          CALL PHSYMCUT(J,CXMN,FELD,Z,WSAVE,IFAX, &
190               &MAUF,MNAUF,MAXL,MAXB,MLEVEL,MANF)
191
192       ENDDO
193       !$OMP END PARALLEL DO
194    ELSE
195       !$OMP PARALLEL DO
196       DO J=1,MAXB
197          CALL PHGPNS(CXMN,FELD,Z,WSAVE,IFAX, &
198               &J,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
199       ENDDO
200       !$OMP END PARALLEL DO
201
202    ENDIF
203
204
205    RETURN
206  END SUBROUTINE PHGCUT
207
208  SUBROUTINE PHSYMCUT(J,CXMN,FELD,Z,WSAVE,IFAX,&
209       &MAUF,MNAUF,MAXL,MAXB,MLEVEL,MANF)
210
211    IMPLICIT NONE
212
213    INTEGER, INTENT(IN)     :: J
214    REAL,    INTENT(IN)     :: CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
215    REAL,    INTENT(OUT)    :: FELD(MAXL,MAXB,MLEVEL)
216    REAL,    INTENT(IN)     :: Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB)
217    REAL,    INTENT(INOUT)  :: WSAVE(4*MAUF+15)
218    INTEGER, INTENT(IN)     :: IFAX(10)
219    INTEGER, INTENT(IN)     :: MAUF,MNAUF,MAXL,MAXB,MLEVEL,MANF
220   
221    !    FELD DER FOURIERKOEFFIZIENTEN
222
223    REAL :: CXM(0:MAXAUF-1),CXMA(0:MAXAUF-1)
224
225    !    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
226    REAL SCR,SCI,ACR,ACI
227
228
229    INTEGER            :: L, LL, LLP, LLS, LLPS, I, K
230
231    DO L=1,MLEVEL
232       LL=0
233       LLP=0
234       DO  I=0,MNAUF
235          SCR=0.D0
236          SCI=0.D0
237          ACR=0.D0
238          ACI=0.D0
239          LLS=LL
240          LLPS=LLP
241          !     Innerste Schleife aufgespalten um if-Abfrage zu sparen
242          DO  K=I,MNAUF,2
243             SCR=SCR+Z(LLP,J)*CXMN(2*LL,L)
244             SCI=SCI+Z(LLP,J)*CXMN(2*LL+1,L)
245             LL=LL+2           
246             LLP=LLP+2           
247          ENDDO
248          LL=LLS+1
249          LLP=LLPS+1
250          DO  K=I+1,MNAUF,2
251             ACR=ACR+Z(LLP,J)*CXMN(2*LL,L)
252             ACI=ACI +Z(LLP,J)*CXMN(2*LL+1,L)
253             LL=LL+2           
254             LLP=LLP+2           
255          ENDDO
256          LL=LLS+MNAUF-I+1
257          LLP=LLPS+MNAUF-I+3
258          CXM(2*I)=SCR+ACR
259          CXM(2*I+1)=SCI+ACI
260          CXMA(2*I)=SCR-ACR
261          CXMA(2*I+1)=SCI-ACI
262       ENDDO
263
264       CALL RFOURTR(CXM,WSAVE,IFAX,MNAUF,MAUF)
265       DO  I=0,MAXL-1
266          IF(MANF+I.LE.MAUF) THEN
267             FELD(I+1,J,L)=CXM(MANF+I-1)
268          ELSE
269             FELD(I+1,J,L)=CXM(MANF-MAUF+I-1)
270          ENDIF
271       ENDDO
272       CALL RFOURTR(CXMA,WSAVE,IFAX,MNAUF,MAUF)
273       DO  I=0,MAXL-1
274          IF(MANF+I.LE.MAUF) THEN
275             FELD(I+1,MAXB+1-J,L)=CXMA(MANF+I-1)
276          ELSE
277             FELD(I+1,MAXB+1-J,L)=CXMA(MANF-MAUF+I-1)
278          ENDIF
279       ENDDO
280    ENDDO
281
282  END SUBROUTINE PHSYMCUT
283
284  SUBROUTINE PHGPNS(CXMN,FELD,Z,WSAVE,IFAX, &
285       J,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
286
287    IMPLICIT NONE
288    REAL,    intent(in)    :: CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
289    REAL,    intent(out)   :: FELD(MAXL,MAXB,MLEVEL)
290    REAL,    intent(in)    :: Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB)
291    REAL,    intent(inout) :: WSAVE(4*MAUF+15)
292    INTEGER, intent(in)    :: IFAX(10)
293    INTEGER, intent(in)    :: MNAUF,MAUF,MANF,J,MAXL,MAXB,MLEVEL
294
295    REAL :: CXM(0:MAXAUF-1)
296
297    INTEGER I,L
298
299    DO L=1,MLEVEL
300       CALL LEGTR(CXMN(:,L),CXM,Z(:,J),MNAUF,MAUF)
301       CALL RFOURTR(CXM,WSAVE,IFAX,MNAUF,MAUF)
302
303       DO I=0,MAXL-1
304          IF(MANF+I.LE.MAUF) THEN
305             FELD(I+1,J,L)=CXM(MANF+I-1)
306          ELSE
307             FELD(I+1,J,L)=CXM(MANF-MAUF+I-1)
308          ENDIF
309       ENDDO
310    ENDDO
311  END SUBROUTINE PHGPNS
312  !
313  SUBROUTINE LEGTR(CXMN,CXM,Z,MNAUF,MAUF)
314    IMPLICIT NONE
315    INTEGER MNAUF,MAUF,LL,LLP,I,J
316    REAL CXM(0:MAXAUF-1)
317    REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1)
318    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2)
319    REAL CI,CR
320    !
321    !     DIESE ROUTINE BERECHNET DIE FOURIERKOEFFIZIENTEN CXM
322    !
323    LL=0
324    LLP=0
325    DO  I=0,MNAUF
326       CR=0.D0
327       CI=0.D0
328       DO  J=I,MNAUF
329          CR=CR+Z(LLP)*CXMN(2*LL)
330          CI=CI+Z(LLP)*CXMN(2*LL+1)
331          LL=LL+1
332          LLP=LLP+1
333       ENDDO
334       LLP=LLP+2
335       CXM(2*I)=CR
336       CXM(2*I+1)=CI
337    ENDDO
338
339  END SUBROUTINE LEGTR
340  !
341  !     
342  SUBROUTINE RFOURTR(CXM,TRIGS,IFAX,MNAUF,MAXL)
343    !     BERECHNET DIE FOURIERSUMME MIT EINEM FFT-ALGORITHMUS
344    IMPLICIT NONE
345
346    INTEGER, INTENT(IN) :: MNAUF, MAXL
347    INTEGER,intent(IN)  :: IFAX(10)
348    REAL, intent(inout) :: CXM(0:MAXAUF-1)
349    REAL                :: WSAVE(2*MAXL),TRIGS(2*MAXL)
350
351    INTEGER :: I
352
353    DO I=MNAUF+1,MAXL-1
354       CXM(2*I)=0.0
355       CXM(2*I+1)=0.0
356    ENDDO
357    CALL FFT99(CXM,WSAVE,TRIGS,IFAX,1,1,MAXL,1,1)
358    DO I=0,MAXL-1
359       CXM(I)=CXM(I+1)
360    ENDDO
361
362    RETURN
363  END SUBROUTINE RFOURTR
364  !     
365  !   
366  SUBROUTINE GAULEG(X1,X2,X,W,N)
367  ! From numerical recipes
368  ! Given the lower and upper limits of integration X1 and X2,
369  ! this routine returns arrayx X and W of length N containing the
370  ! abscissa and weights of the Gauss-Legendre N-point quadrature
371  ! formula. The parameter EPS is the relative precision. Note that
372  ! internal computations are done in double precision.
373    IMPLICIT NONE
374
375    INTEGER, INTENT(IN)             :: N
376    REAL, INTENT(IN)                :: X1, X2
377    REAL, DIMENSION(N), INTENT(OUT) :: X, W
378
379    REAL                            :: PI
380    REAL, PARAMETER                 :: EPS=3.D-14
381    INTEGER                         :: I,J,M
382    REAL                            :: P1, P2, P3
383    REAL                            :: PP
384    REAL                            :: Z1, Z
385    REAL                            :: XM, XL
386
387    PI=ACOS(-1.D0)
388   
389    M=(N+1)/2
390    XM=0.5D0*(X2+X1)
391    XL=0.5D0*(X2-X1)
392    DO  I=1,M
393       Z=DCOS(3.141592654D0*(I-.25D0)/(N+.5D0))
394       DO
395          P1=1.D0
396          P2=0.D0
397          DO  J=1,N
398             P3=P2
399             P2=P1
400             P1=((2.D0*J-1.D0)*Z*P2-(J-1.D0)*P3)/J
401          ENDDO
402          PP=N*(Z*P1-P2)/(Z*Z-1.D0)
403          Z1=Z
404          Z=Z1-P1/PP
405          IF (ABS(Z-Z1).LE.EPS) EXIT
406       ENDDO
407
408       X(I)=XM-XL*Z
409       X(N+1-I)=XM+XL*Z
410       W(I)=2.D0*XL/((1.D0-Z*Z)*PP*PP)
411       W(N+1-I)=W(I)
412    ENDDO
413  END SUBROUTINE GAULEG
414  !
415  !
416  SUBROUTINE PLGNFA(LL,X,Z)
417    !
418    ! PLGNFA is associated with normalization of
419    ! legendre functions OF P00(X) to PLL(X)
420    ! and write result in Z
421    ! The polynomials are indexed as the ECMWF, i.e.
422    ! P00, P10, P11, P20, P21, P22, ...
423    ! This routine is analogous to PLGNDN
424    ! X is the cosinus of the zenith angle or
425    ! sinus of latitude
426    !
427    IMPLICIT NONE
428
429    INTEGER, INTENT(IN)                                :: LL
430    REAL, INTENT(IN)                                   :: X
431    REAL, DIMENSION (0:((LL+3)*(LL+4))/2), INTENT(OUT) :: Z
432
433    REAL                 :: FACT, POT, SOMX2
434    REAL                 :: DJ, DK, DDK
435    INTEGER              :: I, J, K, L
436
437    !
438    L=LL+2
439    I=1
440    Z(0)=1.D0
441    FACT=1.D0
442    POT=1.D0
443    SOMX2=DSQRT(1.D0-X*X)
444    DO  J=0,L
445       DJ=DBLE(J)
446       IF(J.GT.0) THEN
447          FACT=FACT*(2.D0*DJ-1.D0)/(2.D0*DJ)
448          POT=POT*SOMX2
449          Z(I)=DSQRT((2.D0*DJ+1.D0)*FACT)*POT
450          I=I+1
451       ENDIF
452       IF(J.LT.L) THEN
453          Z(I)=X* &
454               &DSQRT((4.D0*DJ*DJ+8.D0*DJ+3.D0)/(2.D0*DJ+1.D0))*Z(I-1)
455          I=I+1
456       ENDIF
457       DK=DJ+2.D0
458       DO  K=J+2,L
459          DDK=(DK*DK-DJ*DJ)
460          Z(I)=X*DSQRT((4.D0*DK*DK-1.D0)/DDK)*Z(I-1)-       &
461               &    DSQRT(((2.D0*DK+1.D0)*(DK-DJ-1.D0)*(DK+DJ-1.D0))/ &
462               &    ((2.D0*DK-3.D0)*DDK))*Z(I-2)
463          DK=DK+1.D0
464          I=I+1
465       ENDDO
466    ENDDO
467    RETURN
468  END SUBROUTINE PLGNFA
469
470
471  SUBROUTINE DPLGND(MNAUF,Z,DZ)
472    !
473    ! DPLGND BERECHNET DIE ABLEITUNG DER NORMIERTEN ASSOZIIERTEN
474    ! LEGENDREFUNKTIONEN VON P00(X) BIS PLL(X)
475    ! UND SCHREIBT SIE IN DAS FELD DZ
476    ! DIE REIHENFOLGE IST
477    ! P00(X),P01(X),P11(X),P02(X),P12(X),P22(X),..PLL(X)
478    !
479    IMPLICIT NONE
480    INTEGER, INTENT(IN) :: MNAUF
481    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2)
482    REAL DZ(0:((MNAUF+2)*(MNAUF+3))/2)
483
484    INTEGER   :: I, J, LLP, LLH
485    REAL      :: WURZELA, WURZELB
486    !
487    IF(Z(0).NE.1.D0) THEN
488       WRITE(*,*) 'DPLGND: Z(0) must be 1.0'
489       STOP
490    ENDIF
491
492    LLP=0
493    LLH=0
494    DO  I=0,MNAUF+1
495       DO  J=I,MNAUF+1
496          IF(I.EQ.J) THEN
497             WURZELA=DSQRT(DBLE((J+1)*(J+1)-I*I)/DBLE(4*(J+1)*(J+1)-1))
498             DZ(LLH)=DBLE(J)*WURZELA*Z(LLP+1)
499          ELSE
500             WURZELB=DSQRT(DBLE((J+1)*(J+1)-I*I)/DBLE(4*(J+1)*(J+1)-1))
501             DZ(LLH)=DBLE(J)*WURZELB*Z(LLP+1)-DBLE(J+1)*WURZELA*Z(LLP-1)
502             WURZELA=WURZELB
503          ENDIF
504          LLH=LLH+1             
505          LLP=LLP+1
506       ENDDO
507       LLP=LLP+1
508    ENDDO
509  END SUBROUTINE DPLGND
510
511
512  ! Spectral Filter of Sardeshmukh and Hoskins (1984, MWR)
513  ! MM=Spectral truncation of field
514  ! MMAX= Spectral truncation of filter
515  !
516  SUBROUTINE SPFILTER(FELDMN,MM,MMAX)
517
518    IMPLICIT NONE
519
520    INTEGER, INTENT(IN) :: MM,MMAX
521    REAL, INTENT(INOUT) :: FELDMN(0:(MM+1)*(MM+2)-1)
522! local variables
523    INTEGER             :: I,J,K,L
524    REAL                :: KMAX,SMAX,FAK
525
526    SMAX=0.1
527    KMAX=-ALOG(SMAX)
528    KMAX=KMAX/(float(MMAX)*float(MMAX+1))**2
529    l=0
530    do i=0,MM
531       do j=i,MM
532          if(j .le. MMAX) then
533             fak=1.0
534             feldmn(2*l)=feldmn(2*l)*fak
535             feldmn(2*l+1)=feldmn(2*l+1)*fak
536          else
537             feldmn(2*l)=0.
538             feldmn(2*l+1)=0.
539          endif
540          l=l+1
541       enddo
542    enddo
543  END SUBROUTINE SPFILTER
544
545END MODULE PHTOGR
546   
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG