Index: preproc/src/phgrreal.f90
===================================================================
--- preproc/src/phgrreal.f90	(revision 57637935b8c736ecb0a376df5e8f6a044bbd05af)
+++ preproc/src/phgrreal.f90	(revision 57637935b8c736ecb0a376df5e8f6a044bbd05af)
@@ -0,0 +1,580 @@
+MODULE PHTOGR
+
+  INTEGER, PARAMETER :: MAXAUF=36000
+
+CONTAINS
+
+  SUBROUTINE PHGR213(CXMN,FELD,WSAVE,IFAX,Z,MLAT,MNAUF, &
+    &MAXL,MAXB,MLEVEL)
+
+    !     DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
+    !     FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
+    !     RAUM AUF DAS REDUZIERTE GAUSS'SCHE GITTER DURCH
+    !
+    !     CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
+    !              CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
+    !     FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
+    !	WSAVE  = Working Array fuer Fouriertransformation
+    !     Z 	 = LEGENDREFUNKTIONSWERTE
+    !
+    !     MNAUF    ANZAHL DER FOURIERKOEFFIZIENTEN
+    !     MAXL     ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
+    !     MAXB     ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
+    !     MLEVEL   ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
+    !
+    IMPLICIT NONE
+
+    !			Anzahl der Gitterpunkte auf jedem Breitenkreis
+    INTEGER MLAT(MAXB/2)
+    INTEGER K,MAXL,MAXB,MLEVEL,MNAUF
+    INTEGER IND(MAXB)
+
+
+    !    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
+    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB/2)
+
+    REAL, INTENT(IN) ::  CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
+    REAL FELD(MAXL,MLEVEL)
+    REAL WSAVE(8*MAXB+15,MAXB/2)
+    INTEGER :: IFAX(10,MAXB)
+
+    IND(1)=0
+    DO  K=2,MAXB/2
+       IND(K)=IND(K-1)+MLAT(K-1)
+    ENDDO
+
+    !$OMP PARALLEL DO SCHEDULE(DYNAMIC)
+    DO  K=1,MAXB/2
+       CALL PHSYM(K,IND,CXMN,FELD,Z,WSAVE,IFAX,MLAT, &
+       &MNAUF,MAXL,MAXB,MLEVEL)
+
+    ENDDO
+    !$OMP END PARALLEL DO
+
+  END SUBROUTINE PHGR213
+  !
+  !
+  SUBROUTINE PHSYM(K,IND,CXMN,FELD,Z,WSAVE,IFAX,MLAT, &
+       &MNAUF,MAXL,MAXB,MLEVEL)
+
+    IMPLICIT NONE
+
+    INTEGER MLAT(MAXB/2)
+    INTEGER K,L,I,J,LLS,LLPS,LL,LLP,MAXL,MAXB,MLEVEL,MNAUF
+    INTEGER IND(MAXB)
+    INTEGER :: IFAX(10,MAXB)
+
+
+    !    FELD DER FOURIERKOEFFIZIENTEN
+    REAL :: CXMS(0:MAXAUF-1),CXMA(0:MAXAUF-1)
+
+    !    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
+    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB/2)
+    REAL ACR,ACI,SCR,SCI
+
+    REAL, INTENT(IN) ::  CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
+    REAL FELD(MAXL,MLEVEL)
+    REAL WSAVE(8*MAXB+15,MAXB/2)
+
+    DO  L=1,MLEVEL
+       LL=0
+       LLP=0
+       DO  I=0,MNAUF
+          SCR=0.D0
+          SCI=0.D0
+          ACR=0.D0
+          ACI=0.D0
+          LLS=LL
+          LLPS=LLP
+          IF(2*I+1.LT.MLAT(K)) THEN
+             !	Innerste Schleife aufgespalten um if-Abfrage zu sparen
+             DO J=I,MNAUF,2
+                SCR=SCR+Z(LLP,K)*CXMN(2*LL,L)
+                SCI=SCI+Z(LLP,K)*CXMN(2*LL+1,L)
+                LL=LL+2            
+                LLP=LLP+2            
+             ENDDO
+             LL=LLS+1
+             LLP=LLPS+1
+             DO  J=I+1,MNAUF,2
+                ACR=ACR+Z(LLP,K)*CXMN(2*LL,L)
+                ACI=ACI+Z(LLP,K)*CXMN(2*LL+1,L)
+                LL=LL+2            
+                LLP=LLP+2            
+             ENDDO
+          ENDIF
+          LL=LLS+(MNAUF-I+1)
+          LLP=LLPS+(MNAUF-I+3)
+          CXMS(2*I)=SCR+ACR
+          CXMS(2*I+1)=SCI+ACI
+          CXMA(2*I)=SCR-ACR
+          CXMA(2*I+1)=SCI-ACI
+       ENDDO
+       CALL RFOURTR(CXMS,WSAVE(:,K),IFAX(:,K),MNAUF,&
+            &MLAT(K))
+       FELD(IND(k)+1:IND(K)+MLAT(K),L)=CXMS(0:MLAT(K)-1)
+       CALL RFOURTR(CXMA,&
+            &WSAVE(:,K),IFAX(:,K),MNAUF,MLAT(K))
+       FELD(MAXL-IND(k)-MLAT(K)+1:MAXL-IND(k),L)=CXMA(0:MLAT(K)-1)
+
+    ENDDO
+
+  END SUBROUTINE PHSYM
+
+  SUBROUTINE PHGCUT(CXMN,FELD,WSAVE,IFAX,Z,&
+       &                  MNAUF,MMAX,MAUF,MANF,MAXL,MAXB,MLEVEL)
+
+    !     DIE ROUTINE FUEHRT EINE TRANSFORMATION EINER
+    !     FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
+    !     RAUM AUF KUGELKOORDINATEN DURCH. Es kann ein Teilausschnitt
+    !			Der Erde angegeben werden. Diese Routine ist langsamer als
+    !			phgrph
+    !
+    !     CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
+    !              CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
+    !     FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
+    !     Z      = SINUS DER GEOGRAFISCHEN BREITEN
+    !
+    !     MNAUF    ANZAHL DER FOURIERKOEFFIZIENTEN
+    !     MAUF     ANZAHL DER LAENGEN UND DER FOURIERKOEFFIZIENTEN
+    !     MANF     ANFANG DES LAENGENBEREICHS FUER DAS GITTER,
+    !              AUF DAS INTERPOLIERT WERDEN SOLL
+    !     MAXL     ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
+    !     MAXB     ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
+    !     MLEVEL   ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
+    !
+    IMPLICIT NONE
+    INTEGER, INTENT(IN)  :: MNAUF, MMAX, MAUF, MANF, MAXL, MAXB, MLEVEL
+    REAL, INTENT(OUT)    :: FELD(MAXL,MAXB,MLEVEL)
+    INTEGER, INTENT(INOUT) :: IFAX(10)
+
+    !    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
+    REAL, INTENT(IN)     :: Z(0:((MMAX+3)*(MMAX+4))/2,MAXB)
+    REAL, INTENT(IN)     :: CXMN(0:(MMAX+1)*(MMAX+2)-1,MLEVEL)
+    REAL, INTENT(INOUT)  :: WSAVE(4*MAUF+15)  ! work array previously
+                                              ! initialized with set99 
+                                              ! (will be used in fft99)
+    INTEGER              :: J
+    LOGICAL              :: SYM
+
+
+    IF(MAUF.LE.MNAUF) WRITE(*,*) 'TOO COARSE LONGITUDE RESOLUTION'
+    IF((MANF.LT.1).OR.(MAXL.LT.1).OR. &
+         &       (MANF.GT.MAUF).OR.(MAXL.GT.MAUF)) THEN
+       WRITE(*,*) 'WRONG LONGITUDE RANGE',MANF,MAXL
+       STOP
+    ENDIF
+
+    ! Pruefe, ob Ausgabegitter symmetrisch zum Aequator ist
+    ! Wenn ja soll Symmetrie der Legendrepolynome ausgenutzt werden
+    ! Check if the output grid is symmetrical to the equator
+    ! If yes we exploit it
+
+    IF(MAXB .GT. 4) THEN
+       SYM=.TRUE.
+       DO  J=5,5
+          IF(ABS(ABS(Z(100,J))-ABS(Z(100,MAXB+1-J))).GT.1E-11) THEN
+             SYM=.FALSE.
+          ENDIF
+       ENDDO
+       WRITE(*,*) 'Symmetrisch: ',SYM
+    ELSE
+       SYM=.FALSE.
+    ENDIF
+
+
+    IF(SYM) THEN
+       !$OMP PARALLEL DO 
+       DO J=1,(MAXB+1)/2
+          CALL PHSYMCUT(J,CXMN,FELD,Z,WSAVE,IFAX, &
+               &MAUF,MNAUF,MAXL,MAXB,MLEVEL,MANF)
+
+       ENDDO
+       !$OMP END PARALLEL DO
+    ELSE
+       !$OMP PARALLEL DO 
+       DO J=1,MAXB
+          CALL PHGPNS(CXMN,FELD,Z,WSAVE,IFAX, &
+               &J,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
+       ENDDO
+       !$OMP END PARALLEL DO
+
+    ENDIF
+
+
+    RETURN
+  END SUBROUTINE PHGCUT
+
+  SUBROUTINE PHSYMCUT(J,CXMN,FELD,Z,WSAVE,IFAX,&
+       &MAUF,MNAUF,MAXL,MAXB,MLEVEL,MANF)
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN)     :: J
+    REAL,    INTENT(IN)     :: CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
+    REAL,    INTENT(OUT)    :: FELD(MAXL,MAXB,MLEVEL)
+    REAL,    INTENT(IN)     :: Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB)
+    REAL,    INTENT(INOUT)  :: WSAVE(4*MAUF+15)
+    INTEGER, INTENT(IN)     :: IFAX(10)
+    INTEGER, INTENT(IN)     :: MAUF,MNAUF,MAXL,MAXB,MLEVEL,MANF
+    
+    !    FELD DER FOURIERKOEFFIZIENTEN
+
+    REAL :: CXM(0:MAXAUF-1),CXMA(0:MAXAUF-1)
+
+    !    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
+    REAL SCR,SCI,ACR,ACI
+
+
+    INTEGER            :: L, LL, LLP, LLS, LLPS, I, K
+
+    DO L=1,MLEVEL
+       LL=0
+       LLP=0
+       DO  I=0,MNAUF
+          SCR=0.D0
+          SCI=0.D0
+          ACR=0.D0
+          ACI=0.D0
+          LLS=LL
+          LLPS=LLP
+          !	Innerste Schleife aufgespalten um if-Abfrage zu sparen
+          DO  K=I,MNAUF,2
+             SCR=SCR+Z(LLP,J)*CXMN(2*LL,L)
+             SCI=SCI+Z(LLP,J)*CXMN(2*LL+1,L)
+             LL=LL+2            
+             LLP=LLP+2            
+          ENDDO
+          LL=LLS+1
+          LLP=LLPS+1
+          DO  K=I+1,MNAUF,2
+             ACR=ACR+Z(LLP,J)*CXMN(2*LL,L)
+             ACI=ACI +Z(LLP,J)*CXMN(2*LL+1,L)
+             LL=LL+2            
+             LLP=LLP+2            
+          ENDDO
+          LL=LLS+MNAUF-I+1
+          LLP=LLPS+MNAUF-I+3
+          CXM(2*I)=SCR+ACR
+          CXM(2*I+1)=SCI+ACI
+          CXMA(2*I)=SCR-ACR
+          CXMA(2*I+1)=SCI-ACI
+       ENDDO
+
+       CALL RFOURTR(CXM,WSAVE,IFAX,MNAUF,MAUF)
+       DO  I=0,MAXL-1
+          IF(MANF+I.LE.MAUF) THEN
+             FELD(I+1,J,L)=CXM(MANF+I-1)
+          ELSE
+             FELD(I+1,J,L)=CXM(MANF-MAUF+I-1)
+          ENDIF
+       ENDDO
+       CALL RFOURTR(CXMA,WSAVE,IFAX,MNAUF,MAUF)
+       DO  I=0,MAXL-1
+          IF(MANF+I.LE.MAUF) THEN
+             FELD(I+1,MAXB+1-J,L)=CXMA(MANF+I-1)
+          ELSE
+             FELD(I+1,MAXB+1-J,L)=CXMA(MANF-MAUF+I-1)
+          ENDIF
+       ENDDO
+    ENDDO
+
+  END SUBROUTINE PHSYMCUT
+
+  SUBROUTINE PHGPNS(CXMN,FELD,Z,WSAVE,IFAX, &
+       J,MNAUF,MAUF,MANF,MAXL,MAXB,MLEVEL)
+
+    IMPLICIT NONE
+    REAL,    intent(in)    :: CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
+    REAL,    intent(out)   :: FELD(MAXL,MAXB,MLEVEL)
+    REAL,    intent(in)    :: Z(0:((MNAUF+3)*(MNAUF+4))/2,MAXB)
+    REAL,    intent(inout) :: WSAVE(4*MAUF+15)
+    INTEGER, intent(in)    :: IFAX(10)
+    INTEGER, intent(in)    :: MNAUF,MAUF,MANF,J,MAXL,MAXB,MLEVEL
+
+    REAL :: CXM(0:MAXAUF-1)
+
+    INTEGER I,L
+
+    DO L=1,MLEVEL
+       CALL LEGTR(CXMN(:,L),CXM,Z(:,J),MNAUF,MAUF)
+       CALL RFOURTR(CXM,WSAVE,IFAX,MNAUF,MAUF)
+
+       DO I=0,MAXL-1
+          IF(MANF+I.LE.MAUF) THEN
+             FELD(I+1,J,L)=CXM(MANF+I-1)
+          ELSE
+             FELD(I+1,J,L)=CXM(MANF-MAUF+I-1)
+          ENDIF
+       ENDDO
+    ENDDO
+  END SUBROUTINE PHGPNS
+  !
+  SUBROUTINE LEGTR(CXMN,CXM,Z,MNAUF,MAUF)
+    IMPLICIT NONE
+    INTEGER MNAUF,MAUF,LL,LLP,I,J
+    REAL CXM(0:MAXAUF-1)
+    REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1)
+    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2)
+    REAL CI,CR
+    !
+    !     DIESE ROUTINE BERECHNET DIE FOURIERKOEFFIZIENTEN CXM
+    !
+    LL=0
+    LLP=0
+    DO  I=0,MNAUF
+       CR=0.D0
+       CI=0.D0
+       DO  J=I,MNAUF
+          CR=CR+Z(LLP)*CXMN(2*LL)
+          CI=CI+Z(LLP)*CXMN(2*LL+1)
+          LL=LL+1
+          LLP=LLP+1
+       ENDDO
+       LLP=LLP+2
+       CXM(2*I)=CR
+       CXM(2*I+1)=CI
+    ENDDO
+
+  END SUBROUTINE LEGTR
+  !
+  !     
+  SUBROUTINE RFOURTR(CXM,TRIGS,IFAX,MNAUF,MAXL)
+    !     BERECHNET DIE FOURIERSUMME MIT EINEM FFT-ALGORITHMUS
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN) :: MNAUF, MAXL
+    INTEGER,intent(IN)  :: IFAX(10)
+    REAL, intent(inout) :: CXM(0:MAXAUF-1)
+    REAL                :: WSAVE(2*MAXL),TRIGS(2*MAXL)
+
+    INTEGER :: I
+
+    DO I=MNAUF+1,MAXL-1
+       CXM(2*I)=0.0
+       CXM(2*I+1)=0.0
+    ENDDO
+    CALL FFT99(CXM,WSAVE,TRIGS,IFAX,1,1,MAXL,1,1)
+    DO I=0,MAXL-1
+       CXM(I)=CXM(I+1)
+    ENDDO
+
+    RETURN
+  END SUBROUTINE RFOURTR
+  !     
+  !    
+  SUBROUTINE GAULEG(X1,X2,X,W,N)
+  ! From numerical recipes
+  ! Given the lower and upper limits of integration X1 and X2, 
+  ! this routine returns arrayx X and W of length N containing the
+  ! abscissa and weights of the Gauss-Legendre N-point quadrature
+  ! formula. The parameter EPS is the relative precision. Note that
+  ! internal computations are done in double precision.
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN)             :: N
+    REAL, INTENT(IN)                :: X1, X2
+    REAL, DIMENSION(N), INTENT(OUT) :: X, W
+
+    REAL                            :: PI
+    REAL, PARAMETER                 :: EPS=3.D-14
+    INTEGER                         :: I,J,M 
+    REAL                            :: P1, P2, P3
+    REAL                            :: PP
+    REAL                            :: Z1, Z
+    REAL                            :: XM, XL
+
+    PI=ACOS(-1.D0)
+   
+    M=(N+1)/2
+    XM=0.5D0*(X2+X1)
+    XL=0.5D0*(X2-X1)
+    DO  I=1,M
+       Z=DCOS(3.141592654D0*(I-.25D0)/(N+.5D0))
+       DO 
+          P1=1.D0
+          P2=0.D0
+          DO  J=1,N
+             P3=P2
+             P2=P1
+             P1=((2.D0*J-1.D0)*Z*P2-(J-1.D0)*P3)/J
+          ENDDO
+          PP=N*(Z*P1-P2)/(Z*Z-1.D0)
+          Z1=Z
+          Z=Z1-P1/PP
+          IF (ABS(Z-Z1).LE.EPS) EXIT
+       ENDDO
+
+       X(I)=XM-XL*Z
+       X(N+1-I)=XM+XL*Z
+       W(I)=2.D0*XL/((1.D0-Z*Z)*PP*PP)
+       W(N+1-I)=W(I)
+    ENDDO
+  END SUBROUTINE GAULEG
+  !
+  !
+  SUBROUTINE PLGNFA(LL,X,Z)
+    !
+    ! PLGNFA is associated with normalization of 
+    ! legendre functions OF P00(X) to PLL(X)
+    ! and write result in Z
+    ! The polynomials are indexed as the ECMWF, i.e.
+    ! P00, P10, P11, P20, P21, P22, ... 
+    ! This routine is analogous to PLGNDN
+    ! X is the cosinus of the zenith angle or
+    ! sinus of latitude
+    !
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN)                                :: LL
+    REAL, INTENT(IN)                                   :: X
+    REAL, DIMENSION (0:((LL+3)*(LL+4))/2), INTENT(OUT) :: Z
+
+    REAL                 :: FACT, POT, SOMX2
+    REAL                 :: DJ, DK, DDK
+    INTEGER              :: I, J, K, L 
+
+    !
+    L=LL+2
+    I=1
+    Z(0)=1.D0
+    FACT=1.D0
+    POT=1.D0
+    SOMX2=DSQRT(1.D0-X*X)
+    DO  J=0,L
+       DJ=DBLE(J)
+       IF(J.GT.0) THEN
+          FACT=FACT*(2.D0*DJ-1.D0)/(2.D0*DJ)
+          POT=POT*SOMX2
+          Z(I)=DSQRT((2.D0*DJ+1.D0)*FACT)*POT
+          I=I+1
+       ENDIF
+       IF(J.LT.L) THEN
+          Z(I)=X* &
+               &DSQRT((4.D0*DJ*DJ+8.D0*DJ+3.D0)/(2.D0*DJ+1.D0))*Z(I-1)
+          I=I+1
+       ENDIF
+       DK=DJ+2.D0
+       DO  K=J+2,L
+          DDK=(DK*DK-DJ*DJ)
+          Z(I)=X*DSQRT((4.D0*DK*DK-1.D0)/DDK)*Z(I-1)-       &
+               &    DSQRT(((2.D0*DK+1.D0)*(DK-DJ-1.D0)*(DK+DJ-1.D0))/ &
+               &    ((2.D0*DK-3.D0)*DDK))*Z(I-2)
+          DK=DK+1.D0
+          I=I+1
+       ENDDO
+    ENDDO
+    RETURN
+  END SUBROUTINE PLGNFA
+
+
+  SUBROUTINE DPLGND(MNAUF,Z,DZ)
+    !
+    ! DPLGND BERECHNET DIE ABLEITUNG DER NORMIERTEN ASSOZIIERTEN
+    ! LEGENDREFUNKTIONEN VON P00(X) BIS PLL(X)
+    ! UND SCHREIBT SIE IN DAS FELD DZ
+    ! DIE REIHENFOLGE IST
+    ! P00(X),P01(X),P11(X),P02(X),P12(X),P22(X),..PLL(X)
+    !
+    IMPLICIT NONE
+    INTEGER, INTENT(IN) :: MNAUF
+    REAL Z(0:((MNAUF+3)*(MNAUF+4))/2)
+    REAL DZ(0:((MNAUF+2)*(MNAUF+3))/2)
+
+    INTEGER   :: I, J, LLP, LLH
+    REAL      :: WURZELA, WURZELB
+    !
+    IF(Z(0).NE.1.D0) THEN
+       WRITE(*,*) 'DPLGND: Z(0) must be 1.0'
+       STOP
+    ENDIF
+
+    LLP=0
+    LLH=0
+    DO  I=0,MNAUF+1
+       DO  J=I,MNAUF+1
+          IF(I.EQ.J) THEN
+             WURZELA=DSQRT(DBLE((J+1)*(J+1)-I*I)/DBLE(4*(J+1)*(J+1)-1))
+             DZ(LLH)=DBLE(J)*WURZELA*Z(LLP+1)
+          ELSE
+             WURZELB=DSQRT(DBLE((J+1)*(J+1)-I*I)/DBLE(4*(J+1)*(J+1)-1))
+             DZ(LLH)=DBLE(J)*WURZELB*Z(LLP+1)-DBLE(J+1)*WURZELA*Z(LLP-1)
+             WURZELA=WURZELB
+          ENDIF
+          LLH=LLH+1              
+          LLP=LLP+1
+       ENDDO
+       LLP=LLP+1
+    ENDDO
+  END SUBROUTINE DPLGND
+
+
+  ! Spectral Filter of Sardeshmukh and Hoskins (1984, MWR)
+  ! MM=Spectral truncation of field
+  ! MMAX= Spectral truncation of filter
+  !
+  SUBROUTINE SPFILTER(FELDMN,MM,MMAX)
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN) :: MM,MMAX
+    REAL, INTENT(INOUT) :: FELDMN(0:(MM+1)*(MM+2)-1)
+! local variables
+    INTEGER             :: I,J,K,L
+    REAL                :: KMAX,SMAX,FAK
+
+    SMAX=0.1
+    KMAX=-ALOG(SMAX)
+    KMAX=KMAX/(float(MMAX)*float(MMAX+1))**2
+    l=0
+    do i=0,MM
+       do j=i,MM
+          if(j .le. MMAX) then
+             fak=1.0
+             feldmn(2*l)=feldmn(2*l)*fak
+             feldmn(2*l+1)=feldmn(2*l+1)*fak
+          else
+             feldmn(2*l)=0.
+             feldmn(2*l+1)=0.
+          endif
+          l=l+1
+       enddo
+    enddo
+  END SUBROUTINE SPFILTER
+
+END MODULE PHTOGR
+    
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
