Skip to content
Snippets Groups Projects
Commit e18f4b5f authored by skomo's avatar skomo
Browse files

initial git repo of version 7.0.3

parents
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
src/ftrafo.f 0 → 100644
+ 504
− 0
View file @ e18f4b5f
MODULE FTRAFO
CONTAINS
C
C Implementierung der spektralen Transformationsmethode unter Verwendung
C des reduzierten Gauss'schen Gitters
C
C Berechnung der scale winds aus Vorticity und Divergenz
C uebergibt man in XMN die Divergenz, so wird der divergente Anteil des
C Windes (XPHI=Ud,XPHI=Vd) zurueckgegeben, uebergibt man die Vorticity, so
C erhaelt man den rotationellen Wind (XLAM=Vrot,XPHI=-Urot).
C Summiert man beide, erhaelt man den gesamten Scale wind
C GWSAVE ist ein Hilfsfeld fuer die FFT
C P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
C MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
C MNAUF gibt die spektrale Aufloesung an,
C NI = Anzahl der Gauss'schen Gitterpunkte pro Flaeche
C NJ = Anzahl der Gauss'schen Breiten,
C NK = Anzahl der Niveaus
SUBROUTINE VDTOUV(XMN,XLAM,XPHI,GWSAVE,IFAX,P,MLAT,MNAUF,NI,NJ,NK)
USE PHTOGR
IMPLICIT NONE
INTEGER J,N,NI,NJ,NK,MNAUF,GGIND(NJ/2)
INTEGER MLAT(NJ),IFAX(10,NJ)
REAL XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
REAL P(0:(MNAUF+3)*(MNAUF+4)/2,NJ/2)
REAL H(0:(MNAUF+2)*(MNAUF+3)/2)
REAL XLAM(NI,NK),XPHI(NI,NK)
REAL GWSAVE(8*NJ+15,NJ/2)
REAL SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI
REAL RT,IT
GGIND(1)=0
DO 4 J = 2,NJ/2
GGIND(J)=GGIND(J-1)+MLAT(J-1)
4 CONTINUE
!$OMP PARALLEL DO SCHEDULE(DYNAMIC)
DO 5 J = 1,NJ/2
CALL VDUVSUB(J,XMN,XLAM,XPHI,GWSAVE,IFAX,P,GGIND(J),
*MLAT,MNAUF,NI,NJ,NK)
5 CONTINUE
!$OMP END PARALLEL DO
RETURN
END SUBROUTINE VDTOUV
SUBROUTINE VDUVSUB(J,XMN,XLAM,XPHI,GWSAVE,IFAX,P,
*GGIND,MLAT,MNAUF,NI,NJ,NK)
USE PHTOGR
IMPLICIT NONE
INTEGER J,K,M,N,NI,NJ,NK,MNAUF,GGIND,LL,LLP,LLH,LLS,LLPS,LLHS
INTEGER MLAT(NJ),IFAX(10,NJ)
REAL UFOUC(0:MAXAUF),MUFOUC(0:MAXAUF)
REAL VFOUC(0:MAXAUF),MVFOUC(0:MAXAUF)
REAL XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
REAL P(0:(MNAUF+3)*(MNAUF+4)/2,NJ/2)
REAL H(0:(MNAUF+2)*(MNAUF+3)/2)
REAL XLAM(NI,NK),XPHI(NI,NK)
REAL GWSAVE(8*NJ+15,NJ/2)
REAL ERAD,SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI
REAL FAC(0:MNAUF),RT,IT
ERAD = 6367470.D0
FAC(0)=0.D0
DO 12 N=1,MNAUF
FAC(N)=-ERAD/DBLE(N)/DBLE(N+1)
12 CONTINUE
CALL DPLGND(MNAUF,P(0,J),H)
DO 3 K = 1,NK
LL=0
LLP=0
LLH=0
DO 2 M = 0,MNAUF
SCR=0.D0
SCI=0.D0
ACR=0.D0
ACI=0.D0
MUSCR=0.D0
MUSCI=0.D0
MUACR=0.D0
MUACI=0.D0
LLS=LL
LLPS=LLP
LLHS=LLH
IF(2*M+1.LT.MLAT(J)) THEN
DO 1 N = M,MNAUF,2
RT=XMN(2*LL,K)*FAC(N)
IT=XMN(2*LL+1,K)*FAC(N)
SCR =SCR+ RT*P(LLP,J)
SCI =SCI+ IT*P(LLP,J)
MUACR =MUACR+ RT*H(LLH)
MUACI =MUACI+ IT*H(LLH)
LL=LL+2
LLP=LLP+2
LLH=LLH+2
1 CONTINUE
LL=LLS+1
LLP=LLPS+1
LLH=LLHS+1
DO 11 N = M+1,MNAUF,2
RT=XMN(2*LL,K)*FAC(N)
IT=XMN(2*LL+1,K)*FAC(N)
ACR =ACR+ RT*P(LLP,J)
ACI =ACI+ IT*P(LLP,J)
MUSCR =MUSCR+ RT*H(LLH)
MUSCI =MUSCI+ IT*H(LLH)
LL=LL+2
LLP=LLP+2
LLH=LLH+2
11 CONTINUE
ENDIF
LL=LLS+(MNAUF-M+1)
LLP=LLPS+(MNAUF-M+3)
LLH=LLHS+(MNAUF-M+2)
UFOUC(2*M)=-M*(SCI-ACI)
UFOUC(2*M+1)=M*(SCR-ACR)
VFOUC(2*M)=-M*(SCI+ACI)
VFOUC(2*M+1)=M*(SCR+ACR)
MUFOUC(2*M)=-(MUSCR-MUACR)
MUFOUC(2*M+1)=-(MUSCI-MUACI)
MVFOUC(2*M)=-(MUSCR+MUACR)
MVFOUC(2*M+1)=-(MUSCI+MUACI)
2 CONTINUE
CALL RFOURTR(VFOUC,
*GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
XLAM(GGIND+1:GGIND+MLAT(J),K)=VFOUC(0:MLAT(J)-1)
CALL RFOURTR(UFOUC,
*GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
XLAM(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=UFOUC(0:MLAT(J)-1)
CALL RFOURTR(MVFOUC,
*GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
XPHI(GGIND+1:GGIND+MLAT(J),K)=MVFOUC(0:MLAT(J)-1)
CALL RFOURTR(MUFOUC,
*GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
XPHI(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=MUFOUC(0:MLAT(J)-1)
3 CONTINUE
RETURN
END SUBROUTINE VDUVSUB
C Berechnung des Gradienten eines Skalars aus dem Feld des
C Skalars XMN im Phasenraum. Zurueckgegeben werden die Felder der
C Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter.
C GWSAVE ist ein Hilfsfeld fuer die FFT
C P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
C MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
C MNAUF gibt die spektrale Aufloesung an,
C NI = Anzahl der Gauss'schen Gitterpunkte,
C NJ = Anzahl der Gauss'schen Breiten,
C NK = Anzahl der Niveaus
SUBROUTINE PHGRAD(XMN,XLAM,XPHI,GWSAVE,IFAX,P,H,MLAT,
*MNAUF,NI,NJ,NK)
USE PHTOGR
IMPLICIT NONE
INTEGER J,K,M,N,NI,NJ,NK,MNAUF,GGIND,LL,LLP,LLH,LLS,LLPS,LLHS
INTEGER MLAT(NJ),IFAX(10,NJ)
REAL UFOUC(0:MAXAUF),MUFOUC(0:MAXAUF)
REAL VFOUC(0:MAXAUF),MVFOUC(0:MAXAUF)
REAL XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
REAL P(0:(MNAUF+3)*(MNAUF+4)/2,NJ/2)
REAL H(0:(MNAUF+2)*(MNAUF+3)/2)
REAL XLAM(NI,NK),XPHI(NI,NK)
REAL GWSAVE(8*NJ+15,NJ/2)
REAL ERAD
REAL SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI,RT,IT
ERAD = 6367470.0
GGIND=0
DO 4 J = 1,NJ/2
CALL DPLGND(MNAUF,P(0,J),H)
DO 3 K = 1,NK
LL=0
LLP=0
LLH=0
DO 2 M = 0,MNAUF
SCR=0.D0
SCI=0.D0
ACR=0.D0
ACI=0.D0
MUSCR=0.D0
MUSCI=0.D0
MUACR=0.D0
MUACI=0.D0
LLS=LL
LLPS=LLP
LLHS=LLH
IF(2*M+1.LT.MLAT(J)) THEN
DO 1 N = M,MNAUF,2
RT=XMN(2*LL,K)
IT=XMN(2*LL+1,K)
SCR =SCR+ RT*P(LLP,J)
SCI =SCI+ IT*P(LLP,J)
MUACR =MUACR+RT*H(LLH)
MUACI =MUACI+ IT*H(LLH)
LL=LL+2
LLP=LLP+2
LLH=LLH+2
1 CONTINUE
LL=LLS+1
LLP=LLPS+1
LLH=LLHS+1
DO 11 N = M+1,MNAUF,2
RT=XMN(2*LL,K)
IT=XMN(2*LL+1,K)
ACR =ACR+ RT*P(LLP,J)
ACI =ACI+ IT*P(LLP,J)
MUSCR =MUSCR+ RT*H(LLH)
MUSCI =MUSCI+ IT*H(LLH)
LL=LL+2
LLP=LLP+2
LLH=LLH+2
11 CONTINUE
ENDIF
LL=LLS+(MNAUF-M+1)
LLP=LLPS+(MNAUF-M+3)
LLH=LLHS+(MNAUF-M+2)
UFOUC(2*M)=-M*(SCI-ACI)/ERAD
UFOUC(2*M+1)=M*(SCR-ACR)/ERAD
VFOUC(2*M)=-M*(SCI+ACI)/ERAD
VFOUC(2*M+1)=M*(SCR+ACR)/ERAD
MUFOUC(2*M)=-(MUSCR-MUACR)/ERAD
MUFOUC(2*M+1)=-(MUSCI-MUACI)/ERAD
MVFOUC(2*M)=-(MUSCR+MUACR)/ERAD
MVFOUC(2*M+1)=-(MUSCI+MUACI)/ERAD
2 CONTINUE
CALL RFOURTR(VFOUC,
*GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
XLAM(GGIND+1:GGIND+MLAT(J),K)=VFOUC(0:MLAT(J)-1)
CALL RFOURTR(UFOUC,
*GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
XLAM(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=UFOUC(0:MLAT(J)-1)
CALL RFOURTR(MVFOUC,
*GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
XPHI(GGIND+1:GGIND+MLAT(J),K)=MVFOUC(0:MLAT(J)-1)
CALL RFOURTR(MUFOUC,
*GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
XPHI(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=MUFOUC(0:MLAT(J)-1)
3 CONTINUE
GGIND=GGIND+MLAT(J)
4 CONTINUE
RETURN
END SUBROUTINE PHGRAD
C Berechnung des Gradienten eines Skalars aus dem Feld des
C Skalars XMN im Phasenraum. Zurueckgegeben werden die Felder der
C Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter.
C GWSAVE ist ein Hilfsfeld fuer die FFT
C P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
C MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
C MNAUF gibt die spektrale Aufloesung an,
C NI = Anzahl der Gauss'schen Gitterpunkte,
C NJ = Anzahl der Gauss'schen Breiten,
C NK = Anzahl der Niveaus
SUBROUTINE PHGRACUT(XMN,XLAM,XPHI,GWSAVE,IFAX,P,H,MAUF,
*MNAUF,NI,NJ,MANF,NK)
USE PHTOGR
IMPLICIT NONE
INTEGER J,K,M,N,NI,NJ,NK,MNAUF,GGIND,LL,LLP,LLH,LLS,LLPS,LLHS
INTEGER MAUF,MANF,I,IFAX(10)
REAL UFOUC(0:MAXAUF),MUFOUC(0:MAXAUF)
REAL VFOUC(0:MAXAUF),MVFOUC(0:MAXAUF)
REAL XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
REAL P(0:(MNAUF+3)*(MNAUF+4)/2,NJ)
REAL H(0:(MNAUF+2)*(MNAUF+3)/2)
REAL XLAM(NI,NJ,NK),XPHI(NI,NJ,NK)
REAL HLAM(MAXAUF,2),HPHI(MAXAUF,2)
REAL GWSAVE(4*MAUF+15)
REAL ERAD
REAL SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI,RT,IT
ERAD = 6367470.0
GGIND=0
DO 4 J = 1,NJ
CALL DPLGND(MNAUF,P(0,J),H)
DO 3 K = 1,NK
LL=0
LLP=0
LLH=0
DO 2 M = 0,MNAUF
SCR=0.D0
SCI=0.D0
ACR=0.D0
ACI=0.D0
MUSCR=0.D0
MUSCI=0.D0
MUACR=0.D0
MUACI=0.D0
LLS=LL
LLPS=LLP
LLHS=LLH
IF(2*M+1.LT.MAUF) THEN
DO 1 N = M,MNAUF,2
RT=XMN(2*LL,K)
IT=XMN(2*LL+1,K)
SCR =SCR+ RT*P(LLP,J)
SCI =SCI+ IT*P(LLP,J)
MUACR =MUACR+RT*H(LLH)
MUACI =MUACI+ IT*H(LLH)
LL=LL+2
LLP=LLP+2
LLH=LLH+2
1 CONTINUE
LL=LLS+1
LLP=LLPS+1
LLH=LLHS+1
DO 11 N = M+1,MNAUF,2
RT=XMN(2*LL,K)
IT=XMN(2*LL+1,K)
ACR =ACR+ RT*P(LLP,J)
ACI =ACI+ IT*P(LLP,J)
MUSCR =MUSCR+ RT*H(LLH)
MUSCI =MUSCI+ IT*H(LLH)
LL=LL+2
LLP=LLP+2
LLH=LLH+2
11 CONTINUE
ENDIF
LL=LLS+(MNAUF-M+1)
LLP=LLPS+(MNAUF-M+3)
LLH=LLHS+(MNAUF-M+2)
UFOUC(2*M)=-M*(SCI-ACI)/ERAD
UFOUC(2*M+1)=M*(SCR-ACR)/ERAD
VFOUC(2*M)=-M*(SCI+ACI)/ERAD
VFOUC(2*M+1)=M*(SCR+ACR)/ERAD
MUFOUC(2*M)=-(MUSCR-MUACR)/ERAD
MUFOUC(2*M+1)=-(MUSCI-MUACI)/ERAD
MVFOUC(2*M)=-(MUSCR+MUACR)/ERAD
MVFOUC(2*M+1)=-(MUSCI+MUACI)/ERAD
2 CONTINUE
CALL RFOURTR(VFOUC,
*GWSAVE,IFAX,MNAUF,MAUF,1)
CALL RFOURTR(MVFOUC,
*GWSAVE,IFAX,MNAUF,MAUF,1)
DO 6 I=0,NI-1
IF(MANF+I.LE. MAUF) THEN
XLAM(I+1,J,K)=VFOUC(MANF+I-1)
XPHI(I+1,J,K)=MVFOUC(MANF+I-1)
ELSE
XLAM(I+1,J,K)=VFOUC(MANF-MAUF+I-1)
XPHI(I+1,J,K)=MVFOUC(MANF-MAUF+I-1)
ENDIF
6 CONTINUE
3 CONTINUE
GGIND=GGIND+MAUF
4 CONTINUE
RETURN
END SUBROUTINE PHGRACUT
C Berechnung der Divergenz aus dem Windfeld (U,V)
C im Phasenraum. Zurueckgegeben werden die Felder der
C Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter.
C GWSAVE ist ein Hilfsfeld fuer die FFT
C P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
C MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
C MNAUF gibt die spektrale Aufloesung an,
C NI = Anzahl der Gauss'schen Gitterpunkte,
C NJ = Anzahl der Gauss'schen Breiten,
C NK = Anzahl der Niveaus
C Beachte, dass das Windfeld eine um 1 erhoehte Aufloesung in mu-Richtung hat.
SUBROUTINE CONTGL(PS,DPSDL,DPSDM,DIV,U,V,BREITE,ETA,
*MLAT,A,B,NI,NJ,NK)
IMPLICIT NONE
INTEGER NI,NJ,NK,I,J,K,MLAT(NJ),L
REAL A(NK+1),B(NK+1)
REAL PS(NI),DPSDL(NI),DPSDM(NI)
REAL DIV(NI,NK),U(NI,NK),V(NI,NK),ETA(NI,NK)
REAL BREITE(NJ)
REAL DIVT1,DIVT2,POB,PUN,DPSDT,COSB
L=0
DO 4 J=1,NJ
COSB=(1.0-BREITE(J)*BREITE(J))
DO 3 I=1,MLAT(J)
L=L+1
DIVT1=0.0
DIVT2=0.0
DO 1 K=1,NK
POB=A(K)+B(K)*PS(L)
PUN=A(K+1)+B(K+1)*PS(L)
DIVT1=DIVT1+DIV(L,K)*(PUN-POB)
if(cosb .gt. 0.) then
DIVT2=DIVT2+(B(K+1)-B(K))*PS(L)*
*(U(L,K)*DPSDL(L)+V(L,K)*DPSDM(L))/COSB
endif
ETA(L,K)=-DIVT1-DIVT2
1 CONTINUE
DPSDT=(-DIVT1-DIVT2)/PS(L)
DO 2 K=1,NK
ETA(L,K)=ETA(L,K)-DPSDT*B(K+1)*PS(L)
2 CONTINUE
PS(L)=DPSDT*PS(L)
3 CONTINUE
4 CONTINUE
RETURN
END SUBROUTINE CONTGL
C OMEGA berechnet omega im Hybridkoordinatensystem
C PS ist der Bodendruck,
C DPSDL,DPSDM sind die Komponenten des Gradienten des Logarithmus des
C Bodendrucks
C DIV,U,V sind die horizontale Divergenz und das horizontale Windfeld
C BREITE ist das Feld der Gauss'schen Breiten
C E ist omega,
SUBROUTINE OMEGA(PS,DPSDL,DPSDM,DIV,U,V,BREITE,E,MLAT,A,B,NGI
* ,NGJ,MKK)
IMPLICIT NONE
INTEGER I,J,K,L,NGI,NGJ,MKK,MLAT(NGJ)
REAL PS(NGI),DPSDL(NGI),DPSDM(NGI),A(MKK+1),B(MKK+1)
REAL DIV(NGI,MKK),U(NGI,MKK),V(NGI,MKK),E(NGI,MKK)
REAL BREITE(NGJ)
REAL DIVT1,DIVT2,POB,PUN,DP,X,Y,COSB
REAL DIVT3(MKK+2)
L=0
DO 4 J=1,NGJ
COSB=(1.0-BREITE(J)*BREITE(J))
DO 3 I=1,MLAT(J)
L=L+1
DIVT1=0.0
DIVT2=0.0
DIVT3(1)=0.0
DO 1 K=1,MKK
POB=A(K)+B(K)*PS(L)
PUN=A(K+1)+B(K+1)*PS(L)
DP=PUN-POB
Y=PS(L)*(U(L,K)*DPSDL(L)+V(L,K)*DPSDM(L))/COSB
IF(K.LT.3) THEN
X=0.0
ELSE
X=(B(K+1)-B(K))*Y
ENDIF
DIVT1=DIVT1+DIV(L,K)*DP
DIVT2=DIVT2+X
DIVT3(K+1)=-DIVT1-DIVT2
IF(K.GT.1) THEN
E(L,K) = 0.5*(POB+PUN)/DP*Y*
*((B(K+1)-B(K))+(A(K+1)*B(K)-A(K)*B(K+1))/
*DP*LOG(PUN/POB))
ELSE
E(L,K) = 0.0
ENDIF
E(L,K) = E(L,K)+0.5*(DIVT3(K)+DIVT3(K+1))
1 CONTINUE
3 CONTINUE
4 CONTINUE
RETURN
END SUBROUTINE OMEGA
END MODULE FTRAFO
src/grphreal.f 0 → 100644
+ 188
− 0
View file @ e18f4b5f
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
src/jparams.h 0 → 100644
+ 34
− 0
View file @ e18f4b5f
C
C Parameters
C
INTEGER JP32, JPLONO, J2NFFT, JPFFT, JPLOOK, JPMAX, JPMAXITER
INTEGER JPMXTRY, JPTRNC, JPK, JPTRP1
PARAMETER ( JP32 = 32 )
C
C The following value for JPLONO (2560) will handle regular grids
C from N1 to N720 derived from spectral truncations from T1 to
C T639.
C
Cjdc PARAMETER ( JPLONO = 2560 , J2NFFT = 2 + JPLONO, JPFFT = 12000)
PARAMETER ( JPLONO = 6000 , J2NFFT = 2 + JPLONO, JPFFT = 12000)
PARAMETER ( JPLOOK = 50)
PARAMETER ( JPMAX = 2048 )
PARAMETER ( JPMAXITER = 10)
PARAMETER ( JPMXTRY = 3 )
PARAMETER ( JPTRNC = 2047, JPK = (JPTRNC + 1)*(JPTRNC + 4) )
PARAMETER ( JPTRP1 = (JPTRNC + 1) )
C
REAL PPEPSA, PPQUART, PPHALF, PPTWO, PP90
PARAMETER ( PPEPSA = 1.0E-6)
PARAMETER ( PPQUART = 0.25E0)
PARAMETER ( PPHALF = 0.5E0)
PARAMETER ( PPTWO = 2.0E0)
PARAMETER ( PP90 = 90.0E0)
C
REAL PPI
PARAMETER ( PPI = 3.14159265358979 )
C
C Debug parameters
C
INTEGER NDBG, NDBGLP
COMMON /JDCNDBG/ NDBG, NDBGLP
src/phgrreal.f 0 → 100644
+ 0
− 0
View file @ e18f4b5f
This diff is collapsed.
src/posnam.f 0 → 100644
+ 25
− 0
View file @ e18f4b5f
SUBROUTINE POSNAM(KULNAM,CDNAML)
!-------------------------------------
!--- position in namelist file.
! author: Mats Hamrud, ECMWF
INTEGER, INTENT(IN) :: KULNAM
CHARACTER*(*), INTENT(IN) :: CDNAML
CHARACTER*120 CLINE
CHARACTER*1 CLTEST
REWIND(KULNAM)
ILEN=LEN(CDNAML)
102 CONTINUE
CLINE=' '
READ(KULNAM,'(A)') CLINE
IND1=INDEX(CLINE,'&'//CDNAML)
IF(IND1.EQ.0) GO TO 102
CLTEST=CLINE(IND1+ILEN+1:IND1+ILEN+1)
IF((LGE(CLTEST,'0').AND.LLE(CLTEST,'9')).OR.
& (LGE(CLTEST,'A').AND.LLE(CLTEST,'Z'))) GO TO 102
BACKSPACE(KULNAM)
RETURN
END SUBROUTINE POSNAM
src/preconvert.f90 0 → 100644
+ 0
− 0
View file @ e18f4b5f
This diff is collapsed.
src/rwGRIB2.f90 0 → 100644
+ 0
− 0
View file @ e18f4b5f
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment