1 | MODULE FTRAFO |
---|
2 | |
---|
3 | CONTAINS |
---|
4 | |
---|
5 | ! Berechnung des Gradienten eines Skalars aus dem Feld des |
---|
6 | ! Skalars XMN im Phasenraum. Zurueckgegeben werden die Felder der |
---|
7 | ! Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter. |
---|
8 | ! GWSAVE ist ein Hilfsfeld fuer die FFT |
---|
9 | ! P enthaelt die assoziierten Legendrepolynome, H deren Ableitung |
---|
10 | ! MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis |
---|
11 | ! MNAUF gibt die spektrale Aufloesung an, |
---|
12 | ! NI = Anzahl der Gauss'schen Gitterpunkte, |
---|
13 | ! NJ = Anzahl der Gauss'schen Breiten, |
---|
14 | ! NK = Anzahl der Niveaus |
---|
15 | |
---|
16 | SUBROUTINE PHGRAD(XMN,XLAM,XPHI,GWSAVE,IFAX,P,H,MLAT, & |
---|
17 | MNAUF,NI,NJ,NK) |
---|
18 | |
---|
19 | USE PHTOGR |
---|
20 | IMPLICIT NONE |
---|
21 | INTEGER J,K,M,N,NI,NJ,NK,MNAUF,GGIND,LL,LLP,LLH,LLS,LLPS,LLHS |
---|
22 | INTEGER MLAT(NJ),IFAX(10,NJ) |
---|
23 | REAL UFOUC(0:MAXAUF),MUFOUC(0:MAXAUF) |
---|
24 | REAL VFOUC(0:MAXAUF),MVFOUC(0:MAXAUF) |
---|
25 | REAL XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK) |
---|
26 | REAL P(0:(MNAUF+3)*(MNAUF+4)/2,NJ/2) |
---|
27 | REAL H(0:(MNAUF+2)*(MNAUF+3)/2) |
---|
28 | REAL XLAM(NI,NK),XPHI(NI,NK) |
---|
29 | REAL GWSAVE(8*NJ+15,NJ/2) |
---|
30 | REAL ERAD |
---|
31 | REAL SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI,RT,IT |
---|
32 | |
---|
33 | ERAD = 6367470.0 |
---|
34 | |
---|
35 | GGIND=0 |
---|
36 | DO J = 1,NJ/2 |
---|
37 | CALL DPLGND(MNAUF,P(0,J),H) |
---|
38 | DO K = 1,NK |
---|
39 | LL=0 |
---|
40 | LLP=0 |
---|
41 | LLH=0 |
---|
42 | DO M = 0,MNAUF |
---|
43 | SCR=0.D0 |
---|
44 | SCI=0.D0 |
---|
45 | ACR=0.D0 |
---|
46 | ACI=0.D0 |
---|
47 | MUSCR=0.D0 |
---|
48 | MUSCI=0.D0 |
---|
49 | MUACR=0.D0 |
---|
50 | MUACI=0.D0 |
---|
51 | LLS=LL |
---|
52 | LLPS=LLP |
---|
53 | LLHS=LLH |
---|
54 | IF(2*M+1.LT.MLAT(J)) THEN |
---|
55 | DO N = M,MNAUF,2 |
---|
56 | RT=XMN(2*LL,K) |
---|
57 | IT=XMN(2*LL+1,K) |
---|
58 | SCR =SCR+ RT*P(LLP,J) |
---|
59 | SCI =SCI+ IT*P(LLP,J) |
---|
60 | MUACR =MUACR+RT*H(LLH) |
---|
61 | MUACI =MUACI+ IT*H(LLH) |
---|
62 | LL=LL+2 |
---|
63 | LLP=LLP+2 |
---|
64 | LLH=LLH+2 |
---|
65 | ENDDO |
---|
66 | LL=LLS+1 |
---|
67 | LLP=LLPS+1 |
---|
68 | LLH=LLHS+1 |
---|
69 | DO N = M+1,MNAUF,2 |
---|
70 | RT=XMN(2*LL,K) |
---|
71 | IT=XMN(2*LL+1,K) |
---|
72 | ACR =ACR+ RT*P(LLP,J) |
---|
73 | ACI =ACI+ IT*P(LLP,J) |
---|
74 | MUSCR =MUSCR+ RT*H(LLH) |
---|
75 | MUSCI =MUSCI+ IT*H(LLH) |
---|
76 | LL=LL+2 |
---|
77 | LLP=LLP+2 |
---|
78 | LLH=LLH+2 |
---|
79 | ENDDO |
---|
80 | ENDIF |
---|
81 | LL=LLS+(MNAUF-M+1) |
---|
82 | LLP=LLPS+(MNAUF-M+3) |
---|
83 | LLH=LLHS+(MNAUF-M+2) |
---|
84 | |
---|
85 | UFOUC(2*M)=-M*(SCI-ACI)/ERAD |
---|
86 | UFOUC(2*M+1)=M*(SCR-ACR)/ERAD |
---|
87 | VFOUC(2*M)=-M*(SCI+ACI)/ERAD |
---|
88 | VFOUC(2*M+1)=M*(SCR+ACR)/ERAD |
---|
89 | |
---|
90 | MUFOUC(2*M)=-(MUSCR-MUACR)/ERAD |
---|
91 | MUFOUC(2*M+1)=-(MUSCI-MUACI)/ERAD |
---|
92 | MVFOUC(2*M)=-(MUSCR+MUACR)/ERAD |
---|
93 | MVFOUC(2*M+1)=-(MUSCI+MUACI)/ERAD |
---|
94 | ENDDO |
---|
95 | |
---|
96 | CALL RFOURTR(VFOUC,GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J)) |
---|
97 | XLAM(GGIND+1:GGIND+MLAT(J),K)=VFOUC(0:MLAT(J)-1) |
---|
98 | CALL RFOURTR(UFOUC,GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J)) |
---|
99 | XLAM(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=UFOUC(0:MLAT(J)-1) |
---|
100 | |
---|
101 | CALL RFOURTR(MVFOUC, GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J)) |
---|
102 | XPHI(GGIND+1:GGIND+MLAT(J),K)=MVFOUC(0:MLAT(J)-1) |
---|
103 | CALL RFOURTR(MUFOUC,GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J)) |
---|
104 | XPHI(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=MUFOUC(0:MLAT(J)-1) |
---|
105 | |
---|
106 | ENDDO |
---|
107 | GGIND=GGIND+MLAT(J) |
---|
108 | ENDDO |
---|
109 | |
---|
110 | END SUBROUTINE PHGRAD |
---|
111 | |
---|
112 | ! Berechnung der Divergenz aus dem Windfeld (U,V) |
---|
113 | ! im Phasenraum. Zurueckgegeben werden die Felder der |
---|
114 | ! Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter. |
---|
115 | ! GWSAVE ist ein Hilfsfeld fuer die FFT |
---|
116 | ! P enthaelt die assoziierten Legendrepolynome, H deren Ableitung |
---|
117 | ! MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis |
---|
118 | ! MNAUF gibt die spektrale Aufloesung an, |
---|
119 | ! NI = Anzahl der Gauss'schen Gitterpunkte, |
---|
120 | ! NJ = Anzahl der Gauss'schen Breiten, |
---|
121 | ! NK = Anzahl der Niveaus |
---|
122 | ! Beachte, dass das Windfeld eine um 1 erhoehte Aufloesung in mu-Richtung hat. |
---|
123 | |
---|
124 | SUBROUTINE CONTGL(PS,DPSDL,DPSDM,DIV,U,V,BREITE,ETA,MLAT,A,B,NI,NJ,NK) |
---|
125 | |
---|
126 | IMPLICIT NONE |
---|
127 | |
---|
128 | INTEGER NI,NJ,NK,I,J,K,MLAT(NJ),L |
---|
129 | |
---|
130 | REAL A(NK+1),B(NK+1) |
---|
131 | REAL PS(NI),DPSDL(NI),DPSDM(NI) |
---|
132 | REAL DIV(NI,NK),U(NI,NK),V(NI,NK),ETA(NI,NK) |
---|
133 | REAL BREITE(NJ) |
---|
134 | |
---|
135 | REAL DIVT1,DIVT2,POB,PUN,DPSDT,COSB |
---|
136 | |
---|
137 | L=0 |
---|
138 | DO J=1,NJ |
---|
139 | COSB=(1.0-BREITE(J)*BREITE(J)) |
---|
140 | DO I=1,MLAT(J) |
---|
141 | L=L+1 |
---|
142 | DIVT1=0.0 |
---|
143 | DIVT2=0.0 |
---|
144 | DO K=1,NK |
---|
145 | POB=A(K)+B(K)*PS(L) |
---|
146 | PUN=A(K+1)+B(K+1)*PS(L) |
---|
147 | |
---|
148 | DIVT1=DIVT1+DIV(L,K)*(PUN-POB) |
---|
149 | if(cosb .gt. 0.) then |
---|
150 | DIVT2=DIVT2+(B(K+1)-B(K))*PS(L)* & |
---|
151 | (U(L,K)*DPSDL(L)+V(L,K)*DPSDM(L))/COSB |
---|
152 | endif |
---|
153 | |
---|
154 | ETA(L,K)=-DIVT1-DIVT2 |
---|
155 | ENDDO |
---|
156 | |
---|
157 | DPSDT=(-DIVT1-DIVT2)/PS(L) |
---|
158 | |
---|
159 | DO K=1,NK |
---|
160 | ETA(L,K)=ETA(L,K)-DPSDT*B(K+1)*PS(L) |
---|
161 | ENDDO |
---|
162 | PS(L)=DPSDT*PS(L) |
---|
163 | ENDDO |
---|
164 | ENDDO |
---|
165 | END SUBROUTINE CONTGL |
---|
166 | |
---|
167 | END MODULE FTRAFO |
---|