grphreal.f90 Source File

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