MODULE GRTOPH

  USE PHTOGR

CONTAINS
  !
  SUBROUTINE GRPH213(CXMN,FELD,WSAVE,IFAX,Z,W,MLAT, &
       MNAUF,MAXL,MAXB,MLEVEL)

    !     DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
    !     FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
    !     RAUM AUF KUGELKOORDINATEN DURCH
    !
    !     CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
    !              CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
    !			CXM		 = FOURIERKOEFFIZIENTEN - nur ein Hilfsfeld
    !     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 REAL (A-H,O-Z)


    !			Anzahl der Gitterpunkte pro Breitenkreis des reduzierten 
    !			Gauss'schen Gitters
    INTEGER MLAT(MAXB),ISIZE,IFAX(10,MAXB)

    !    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
    REAL*8 Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)

    !      LOGICAL*1 USED(((216*217)/2+1)*160)

    DIMENSION CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
    REAL FELD(MAXL,MLEVEL)
    DIMENSION WSAVE(8*MAXB+15,MAXB/2)
    REAL*8 W(MAXB) 
    DIMENSION IND(MAXB)


    IND(1)=0
    DO J=2,MAXB/2
       IND(j)=IND(J-1)+MLAT(J-1)
    ENDDO
    !$OMP PARALLEL DO SCHEDULE(DYNAMIC)
    DO  L=1,MLEVEL
       CALL GRPHSUB(L,IND,CXMN,FELD,WSAVE,IFAX,Z,W,MLAT, &
            &MNAUF,MAXL,MAXB,MLEVEL)
    ENDDO
    !$omp end parallel do


    RETURN
  END SUBROUTINE GRPH213
  !
  SUBROUTINE GRPHSUB(L,IND,CXMN,FELD,WSAVE,IFAX,Z,W,MLAT, &
       &MNAUF,MAXL,MAXB,MLEVEL)

    !     DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
    !     FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
    !     RAUM AUF KUGELKOORDINATEN DURCH
    !
    !     CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
    !              CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
    !			CXM		 = FOURIERKOEFFIZIENTEN - nur ein Hilfsfeld
    !     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 REAL (A-H,O-Z)

    !    FELD DER FOURIERKOEFFIZIENTEN
    REAL CXMS(4*(MNAUF+1))
    REAL CXMA(4*(MNAUF+1))
    REAL,ALLOCATABLE :: CXM(:,:)

    !			Anzahl der Gitterpunkte pro Breitenkreis des reduzierten 
    !			Gauss'schen Gitters
    INTEGER MLAT(MAXB),ISIZE

    !    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
    REAL Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)

    !      LOGICAL*1 USED(((216*217)/2+1)*160)

    REAL CXMN(0:(MNAUF+1)*(MNAUF+2)-1,MLEVEL)
    REAL FELD(MAXL,MLEVEL)
    REAL WSAVE(8*MAXB+15,MAXB/2)
    INTEGER IFAX(10,MAXB)
    REAL W(MAXB)
    INTEGER IND(MAXB)

    ALLOCATE(CXM( 4*MAXB,MAXB))
    DO  J=1,MAXB/2
       CXMS(1:MLAT(J))=FELD(IND(J)+1:IND(J)+MLAT(J),L)
       CALL RFOUFTR(CXMS,WSAVE(1,J),IFAX(:,J),MNAUF,MLAT(J),1)
       CXMA(1:MLAT(J))=FELD(MAXL-IND(J)-MLAT(J)+1:MAXL-IND(J),L)
       CALL RFOUFTR(CXMA, &
            &WSAVE(1,J),IFAX(:,J),MNAUF,MLAT(J),1)
       DO  I=1,2*(MNAUF+1)
          CXM(I,J)=CXMS(I)+CXMA(I)
          CXM(I,MAXB+1-J)=CXMS(I)-CXMA(I)
       ENDDO
    ENDDO
    CALL LGTR213(CXMN(0,L),CXM,Z,W,MLAT,MNAUF,MAXB)

    DEALLOCATE(CXM)

    RETURN
  END SUBROUTINE GRPHSUB
  !
  SUBROUTINE LGTR213(CXMN,CXM,Z,W,MLAT,MNAUF,MAXB)
    IMPLICIT REAL (A-H,O-Z)
    INTEGER MLAT(MAXB)
    DIMENSION CXM(0:4*MAXB-1,MAXB)
    DIMENSION CXMN(0:2*(((MNAUF+1)*MNAUF)/2+MNAUF)+1)
    REAL*8 Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)
    REAL*8 W(MAXB),CR,CI,HILF
    LOGICAL EVEN
    !
    !     DIESE ROUTINE BERECHNET DIE KFFKs CXMN
    !
    LL=0
    LLP=0
    DO  I=0,MNAUF
       KM=0
       DO 
          KM=KM+1
          IF(MLAT(KM).GT.2*I) THEN
             EXIT
          ENDIF
       ENDDO
       DO  J=I,MNAUF
          CR=0
          CI=0
          EVEN=MOD(I+J,2).EQ.0
          IF(EVEN) THEN
             DO  K=KM,MAXB/2
                HILF=W(K)*Z(K,LLP)
                CR=CR+CXM(2*I,K)*HILF
                CI=CI+CXM(2*I+1,K)*HILF
             ENDDO
          ELSE
             DO  K=KM,MAXB/2
                HILF=W(K)*Z(K,LLP)
                CR=CR+CXM(2*I,MAXB+1-K)*HILF
                CI=CI+CXM(2*I+1,MAXB+1-K)*HILF
             ENDDO
          ENDIF
          CXMN(2*LL)=CR
          CXMN(2*LL+1)=CI
          LL=LL+1
          LLP=LLP+1
       ENDDO
       LLP=LLP+2
    ENDDO
    RETURN
  END SUBROUTINE LGTR213
  !

  !
  SUBROUTINE RFOUFTR(CXM,TRIGS,IFAX,MNAUF,MAXL,ISIGN)
    !     BERECHNET DIE FOURIERSUMME MIT EINEM FFT-ALGORITHMUS
    IMPLICIT REAL (A-H,O-Z)
    DIMENSION CXM(0:2*MAXL-1)
    DIMENSION FELD(MAXL),TRIGS(2*MAXL)
    DIMENSION WSAVE(MAXAUF)
    INTEGER IFAX(10)


    ! NORMIERUNG...
    WSAVE(1)=CXM(MAXL-1)

    CXM(1:MAXL)=CXM(0:MAXL-1)/2
    CXM(0)=WSAVE(1)/2
    CALL FFT99(CXM,WSAVE,TRIGS,IFAX,1,1,MAXL,1,-1)
    RETURN
  END SUBROUTINE RFOUFTR

END MODULE GRTOPH