grphreal.f Source File

      MODULE GRTOPH

      USE PHTOGR

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

C     DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
C     FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
C     RAUM AUF KUGELKOORDINATEN DURCH
C
C     CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
C              CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
C			CXM		 = FOURIERKOEFFIZIENTEN - nur ein Hilfsfeld
C     FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
C			WSAVE  = Working Array fuer Fouriertransformation
C     Z 	 = LEGENDREFUNKTIONSWERTE
C
C     MNAUF    ANZAHL DER FOURIERKOEFFIZIENTEN
C     MAXL     ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
C     MAXB     ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
C     MLEVEL   ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
C
      IMPLICIT REAL (A-H,O-Z)


C			Anzahl der Gitterpunkte pro Breitenkreis des reduzierten 
C			Gauss'schen Gitters
	INTEGER MLAT(MAXB),ISIZE,IFAX(10,MAXB)
			
C    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
      REAL*8 Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)
      
C      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 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
C
      SUBROUTINE GRPHSUB(L,IND,CXMN,FELD,WSAVE,IFAX,Z,W,MLAT,
     *MNAUF,MAXL,MAXB,MLEVEL)

C     DIE ROUTINE F]HRT EINE TRANSFORMATION EINER
C     FELDVARIABLEN VOM PHASENRAUM IN  DEN PHYSIKALISCHEN
C     RAUM AUF KUGELKOORDINATEN DURCH
C
C     CXMN   = SPEKTRALKOEFFIZIENTEN IN DER REIHENFOLGE
C              CX00,CX01,CX11,CX02,....CXMNAUFMNAUF
C			CXM		 = FOURIERKOEFFIZIENTEN - nur ein Hilfsfeld
C     FELD   = FELD DER METEOROLOGISCHEN VARIABLEN
C			WSAVE  = Working Array fuer Fouriertransformation
C     Z 	 = LEGENDREFUNKTIONSWERTE
C
C     MNAUF    ANZAHL DER FOURIERKOEFFIZIENTEN
C     MAXL     ANZAHL DER FUER DAS GITTER BENUTZTEN LAENGEN
C     MAXB     ANZAHL DER FUER DAS GITTER BENOETIGTEN BREITEN
C     MLEVEL   ANZAHL DER LEVELS, DIE TRANSFORMIERT WERDEN
C
      IMPLICIT REAL (A-H,O-Z)

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

C			Anzahl der Gitterpunkte pro Breitenkreis des reduzierten 
C			Gauss'schen Gitters
			INTEGER MLAT(MAXB),ISIZE
			
C    FELD DER LEGENDREPOLYNOME FUER EINE BREITE
      REAL Z(MAXB/2,0:((MNAUF+3)*(MNAUF+4))/2)
      
C      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
C
        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
C
C     DIESE ROUTINE BERECHNET DIE KFFKs CXMN
C
			LL=0
			LLP=0
      DO 1 I=0,MNAUF
        KM=0
    9   KM=KM+1
        IF(MLAT(KM).LE.2*I) THEN
           GOTO 9
        ENDIF
        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
          ENDIF
    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
C

C
      SUBROUTINE RFOUFTR(CXM,TRIGS,IFAX,MNAUF,MAXL,ISIGN)
C     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)


C 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