text: MODULE GRTOPH
  
  !! Module for transformation from phase space into physical space
  
    USE PHTOGR
  
  CONTAINS
  
    SUBROUTINE GRPH213(CXMN,FELD,WSAVE,IFAX,Z,W,MLAT,MNAUF,MAXL,MAXB,MLEVEL)
  
  ! WRONG>>> 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 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 W(MAXB)
      DIMENSION IND(MAXB)
  
      IND(1)=0
      DO 6 J=2,MAXB/2
        IND(j)=IND(J-1)+MLAT(J-1)
  6   CONTINUE
  !$OMP PARALLEL DO SCHEDULE(DYNAMIC)
      DO 16 L=1,MLEVEL
        CALL GRPHSUB(L,IND,CXMN,FELD,WSAVE,IFAX,Z,W,MLAT,MNAUF,MAXL,MAXB,MLEVEL)
  16  CONTINUE
  !$omp end parallel do
  
      RETURN
  
    END SUBROUTINE GRPH213
  
    SUBROUTINE GRPHSUB(L,IND,CXMN,FELD,WSAVE,IFAX,Z,W,MLAT,MNAUF,MAXL,MAXB,MLEVEL)
  
  !! Carries out the transformation from the phase space into the physical
  !! space on spherical coordinates
  
  ! 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 5 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 4 I=1,2*(MNAUF+1)
          CXM(I,J)=CXMS(I)+CXMA(I)
          CXM(I,MAXB+1-J)=CXMS(I)-CXMA(I)
  4     CONTINUE
  5   CONTINUE
      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)
  
  !! Calculates the KFFKs (Fourier coefficients) CXMN
  
      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
  
  
      LL=0
      LLP=0
      DO 1 I=0,MNAUF
        KM=0
  9     KM=KM+1
        IF (MLAT(KM) .LE. 2*I) THEN
          GOTO 9
        END IF
        DO 2 J=I,MNAUF
          CR=0
          CI=0
          EVEN=MOD(I+J,2) .EQ. 0
          IF (EVEN) THEN
            DO 3 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
  3         CONTINUE
          ELSE
            DO 4 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
  4         CONTINUE
          END IF
  5       CXMN(2*LL)=CR
          CXMN(2*LL+1)=CI
          LL=LL+1
          LLP=LLP+1
  2     CONTINUE
        LLP=LLP+2
  1   CONTINUE
      RETURN
      
    END SUBROUTINE LGTR213
  
    SUBROUTINE RFOUFTR(CXM,TRIGS,IFAX,MNAUF,MAXL,ISIGN)
  
  !! Calculates the Fourier sum with an FFT algorithm
  
      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 CFFTF(MAXL,CXM,WSAVE)
      CALL FFT99(CXM,WSAVE,TRIGS,IFAX,1,1,MAXL,1,-1)
      RETURN
    END SUBROUTINE RFOUFTR
  
  END MODULE GRTOPH
  