source: flex_extract.git/python/pythontest/TestInstallTar/test_untar/src/ftrafo.f @ 2fb99de

ctbtodev
Last change on this file since 2fb99de was 2fb99de, checked in by Anne Philipp <anne.philipp@…>, 6 years ago

introduced config with path definitions and changed py files accordingly; Installation works; some tests were added for tarball making; Problems in submission to ecgate

  • Property mode set to 100644
File size: 12.4 KB
Line 
1        MODULE FTRAFO
2
3        CONTAINS
4
5
6
7C
8C       Implementierung der spektralen Transformationsmethode unter Verwendung
9C       des reduzierten Gauss'schen Gitters
10C       
11C Berechnung der scale winds aus Vorticity und Divergenz
12C uebergibt man in XMN die Divergenz, so wird der divergente Anteil des
13C Windes (XPHI=Ud,XPHI=Vd) zurueckgegeben, uebergibt man die Vorticity, so
14C erhaelt man den rotationellen Wind (XLAM=Vrot,XPHI=-Urot).
15C Summiert man beide, erhaelt man den gesamten Scale wind
16C       GWSAVE ist ein Hilfsfeld fuer die FFT
17C       P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
18C       MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
19C       MNAUF gibt die spektrale Aufloesung an,
20C       NI = Anzahl der Gauss'schen Gitterpunkte pro Flaeche
21C       NJ = Anzahl der Gauss'schen Breiten,
22C       NK = Anzahl der Niveaus
23 
24        SUBROUTINE VDTOUV(XMN,XLAM,XPHI,GWSAVE,IFAX,P,MLAT,MNAUF,NI,NJ,NK)
25
26
27        USE PHTOGR
28
29        IMPLICIT NONE
30        INTEGER   J,N,NI,NJ,NK,MNAUF,GGIND(NJ/2)
31        INTEGER         MLAT(NJ),IFAX(10,NJ)
32        REAL                    XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
33        REAL            P(0:(MNAUF+3)*(MNAUF+4)/2,NJ/2)
34        REAL            H(0:(MNAUF+2)*(MNAUF+3)/2)
35        REAL                    XLAM(NI,NK),XPHI(NI,NK)
36        REAL                    GWSAVE(8*NJ+15,NJ/2)
37        REAL            SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI
38        REAL      RT,IT
39
40        GGIND(1)=0
41        DO 4 J = 2,NJ/2           
42          GGIND(J)=GGIND(J-1)+MLAT(J-1)
434       CONTINUE
44!$OMP PARALLEL DO SCHEDULE(DYNAMIC)
45        DO 5 J = 1,NJ/2
46           CALL VDUVSUB(J,XMN,XLAM,XPHI,GWSAVE,IFAX,P,GGIND(J),
47     *MLAT,MNAUF,NI,NJ,NK)   
48 5      CONTINUE
49!$OMP END PARALLEL DO
50        RETURN
51        END SUBROUTINE VDTOUV
52
53        SUBROUTINE VDUVSUB(J,XMN,XLAM,XPHI,GWSAVE,IFAX,P,
54     *GGIND,MLAT,MNAUF,NI,NJ,NK)
55
56        USE PHTOGR
57
58        IMPLICIT NONE
59        INTEGER   J,K,M,N,NI,NJ,NK,MNAUF,GGIND,LL,LLP,LLH,LLS,LLPS,LLHS
60        INTEGER         MLAT(NJ),IFAX(10,NJ)
61        REAL            UFOUC(0:MAXAUF),MUFOUC(0:MAXAUF)
62        REAL                    VFOUC(0:MAXAUF),MVFOUC(0:MAXAUF)
63        REAL                    XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
64        REAL            P(0:(MNAUF+3)*(MNAUF+4)/2,NJ/2)
65        REAL            H(0:(MNAUF+2)*(MNAUF+3)/2)
66        REAL                    XLAM(NI,NK),XPHI(NI,NK)
67        REAL                    GWSAVE(8*NJ+15,NJ/2)
68        REAL            ERAD,SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI
69        REAL      FAC(0:MNAUF),RT,IT
70 
71       
72        ERAD = 6367470.D0
73       
74        FAC(0)=0.D0
75        DO 12 N=1,MNAUF
76          FAC(N)=-ERAD/DBLE(N)/DBLE(N+1)
7712      CONTINUE
78
79         CALL DPLGND(MNAUF,P(0,J),H)
80         DO 3 K = 1,NK
81          LL=0
82          LLP=0
83          LLH=0
84          DO 2 M = 0,MNAUF
85           SCR=0.D0
86           SCI=0.D0
87           ACR=0.D0
88           ACI=0.D0
89           MUSCR=0.D0
90           MUSCI=0.D0
91           MUACR=0.D0
92           MUACI=0.D0
93           LLS=LL
94           LLPS=LLP
95           LLHS=LLH
96           IF(2*M+1.LT.MLAT(J)) THEN
97              DO 1 N = M,MNAUF,2
98               RT=XMN(2*LL,K)*FAC(N)
99               IT=XMN(2*LL+1,K)*FAC(N)
100               SCR =SCR+ RT*P(LLP,J)
101               SCI =SCI+ IT*P(LLP,J)
102               MUACR =MUACR+ RT*H(LLH)
103               MUACI =MUACI+ IT*H(LLH)
104               LL=LL+2
105               LLP=LLP+2
106               LLH=LLH+2
107 1            CONTINUE
108              LL=LLS+1
109              LLP=LLPS+1
110              LLH=LLHS+1
111              DO 11 N = M+1,MNAUF,2
112               RT=XMN(2*LL,K)*FAC(N)
113               IT=XMN(2*LL+1,K)*FAC(N)
114               ACR =ACR+ RT*P(LLP,J)
115               ACI =ACI+ IT*P(LLP,J)
116               MUSCR =MUSCR+ RT*H(LLH)
117               MUSCI =MUSCI+ IT*H(LLH)
118               LL=LL+2
119               LLP=LLP+2
120               LLH=LLH+2
121 11           CONTINUE
122           ENDIF
123           LL=LLS+(MNAUF-M+1)
124           LLP=LLPS+(MNAUF-M+3)
125           LLH=LLHS+(MNAUF-M+2)
126
127           UFOUC(2*M)=-M*(SCI-ACI)
128           UFOUC(2*M+1)=M*(SCR-ACR)
129           VFOUC(2*M)=-M*(SCI+ACI)
130           VFOUC(2*M+1)=M*(SCR+ACR)
131           
132           MUFOUC(2*M)=-(MUSCR-MUACR)
133           MUFOUC(2*M+1)=-(MUSCI-MUACI)
134           MVFOUC(2*M)=-(MUSCR+MUACR)
135           MVFOUC(2*M+1)=-(MUSCI+MUACI)
136 2        CONTINUE
137                       
138         CALL RFOURTR(VFOUC,
139     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
140         XLAM(GGIND+1:GGIND+MLAT(J),K)=VFOUC(0:MLAT(J)-1)
141         CALL RFOURTR(UFOUC,
142     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
143         XLAM(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=UFOUC(0:MLAT(J)-1)
144                       
145         CALL RFOURTR(MVFOUC,
146     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
147         XPHI(GGIND+1:GGIND+MLAT(J),K)=MVFOUC(0:MLAT(J)-1)
148         CALL RFOURTR(MUFOUC,
149     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
150         XPHI(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=MUFOUC(0:MLAT(J)-1)
151
1523               CONTINUE
153
154      RETURN
155      END SUBROUTINE VDUVSUB
156
157C Berechnung des Gradienten eines Skalars aus dem Feld des
158C       Skalars XMN im Phasenraum. Zurueckgegeben werden die Felder der
159C       Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter.
160C       GWSAVE ist ein Hilfsfeld fuer die FFT
161C       P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
162C       MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
163C       MNAUF gibt die spektrale Aufloesung an, 
164C       NI = Anzahl der Gauss'schen Gitterpunkte,
165C       NJ = Anzahl der Gauss'schen Breiten,
166C       NK = Anzahl der Niveaus
167 
168        SUBROUTINE PHGRAD(XMN,XLAM,XPHI,GWSAVE,IFAX,P,H,MLAT,
169     *MNAUF,NI,NJ,NK)
170
171        USE PHTOGR
172        IMPLICIT NONE
173        INTEGER   J,K,M,N,NI,NJ,NK,MNAUF,GGIND,LL,LLP,LLH,LLS,LLPS,LLHS
174        INTEGER         MLAT(NJ),IFAX(10,NJ)
175        REAL            UFOUC(0:MAXAUF),MUFOUC(0:MAXAUF)
176        REAL                    VFOUC(0:MAXAUF),MVFOUC(0:MAXAUF)
177        REAL                    XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
178        REAL            P(0:(MNAUF+3)*(MNAUF+4)/2,NJ/2)
179        REAL            H(0:(MNAUF+2)*(MNAUF+3)/2)
180        REAL                    XLAM(NI,NK),XPHI(NI,NK)
181        REAL                    GWSAVE(8*NJ+15,NJ/2)
182        REAL      ERAD
183        REAL            SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI,RT,IT
184       
185        ERAD = 6367470.0
186       
187        GGIND=0
188        DO 4 J = 1,NJ/2
189         CALL DPLGND(MNAUF,P(0,J),H)
190         DO 3 K = 1,NK
191          LL=0
192          LLP=0
193          LLH=0
194          DO 2 M = 0,MNAUF
195           SCR=0.D0
196           SCI=0.D0
197           ACR=0.D0
198           ACI=0.D0
199           MUSCR=0.D0
200           MUSCI=0.D0
201           MUACR=0.D0
202           MUACI=0.D0
203           LLS=LL
204           LLPS=LLP
205           LLHS=LLH
206           IF(2*M+1.LT.MLAT(J)) THEN
207              DO 1 N = M,MNAUF,2
208               RT=XMN(2*LL,K)
209               IT=XMN(2*LL+1,K)
210               SCR =SCR+ RT*P(LLP,J)
211               SCI =SCI+ IT*P(LLP,J)
212               MUACR =MUACR+RT*H(LLH)
213               MUACI =MUACI+ IT*H(LLH)
214               LL=LL+2
215               LLP=LLP+2
216               LLH=LLH+2
217 1            CONTINUE
218              LL=LLS+1
219              LLP=LLPS+1
220              LLH=LLHS+1
221              DO 11 N = M+1,MNAUF,2
222               RT=XMN(2*LL,K)
223               IT=XMN(2*LL+1,K)
224               ACR =ACR+ RT*P(LLP,J)
225               ACI =ACI+ IT*P(LLP,J)
226               MUSCR =MUSCR+ RT*H(LLH)
227               MUSCI =MUSCI+ IT*H(LLH)
228               LL=LL+2
229               LLP=LLP+2
230               LLH=LLH+2
231 11           CONTINUE
232           ENDIF
233          LL=LLS+(MNAUF-M+1)
234          LLP=LLPS+(MNAUF-M+3)
235          LLH=LLHS+(MNAUF-M+2)
236
237                UFOUC(2*M)=-M*(SCI-ACI)/ERAD
238                UFOUC(2*M+1)=M*(SCR-ACR)/ERAD
239                VFOUC(2*M)=-M*(SCI+ACI)/ERAD
240                VFOUC(2*M+1)=M*(SCR+ACR)/ERAD
241
242                MUFOUC(2*M)=-(MUSCR-MUACR)/ERAD
243                MUFOUC(2*M+1)=-(MUSCI-MUACI)/ERAD
244                MVFOUC(2*M)=-(MUSCR+MUACR)/ERAD
245                MVFOUC(2*M+1)=-(MUSCI+MUACI)/ERAD
2462       CONTINUE
247
248         CALL RFOURTR(VFOUC,
249     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
250         XLAM(GGIND+1:GGIND+MLAT(J),K)=VFOUC(0:MLAT(J)-1)
251         CALL RFOURTR(UFOUC,
252     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
253         XLAM(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=UFOUC(0:MLAT(J)-1)
254                       
255         CALL RFOURTR(MVFOUC,
256     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
257         XPHI(GGIND+1:GGIND+MLAT(J),K)=MVFOUC(0:MLAT(J)-1)
258         CALL RFOURTR(MUFOUC,
259     *GWSAVE(:,J),IFAX(:,J),MNAUF,MLAT(J),1)
260         XPHI(NI-GGIND-MLAT(J)+1:NI-GGIND,K)=MUFOUC(0:MLAT(J)-1)
261                       
2623               CONTINUE
263                GGIND=GGIND+MLAT(J)
2644       CONTINUE
265
266
267        RETURN
268        END SUBROUTINE PHGRAD
269
270C Berechnung des Gradienten eines Skalars aus dem Feld des
271C       Skalars XMN im Phasenraum. Zurueckgegeben werden die Felder der
272C       Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter.
273C       GWSAVE ist ein Hilfsfeld fuer die FFT
274C       P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
275C       MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
276C       MNAUF gibt die spektrale Aufloesung an,
277C       NI = Anzahl der Gauss'schen Gitterpunkte,
278C       NJ = Anzahl der Gauss'schen Breiten,
279C       NK = Anzahl der Niveaus
280 
281        SUBROUTINE PHGRACUT(XMN,XLAM,XPHI,GWSAVE,IFAX,P,H,MAUF,
282     *MNAUF,NI,NJ,MANF,NK)
283
284        USE PHTOGR
285        IMPLICIT NONE
286        INTEGER   J,K,M,N,NI,NJ,NK,MNAUF,GGIND,LL,LLP,LLH,LLS,LLPS,LLHS
287        INTEGER         MAUF,MANF,I,IFAX(10)
288        REAL            UFOUC(0:MAXAUF),MUFOUC(0:MAXAUF)
289        REAL                    VFOUC(0:MAXAUF),MVFOUC(0:MAXAUF)
290        REAL                    XMN(0:(MNAUF+1)*(MNAUF+2)-1,NK)
291        REAL            P(0:(MNAUF+3)*(MNAUF+4)/2,NJ)
292        REAL            H(0:(MNAUF+2)*(MNAUF+3)/2)
293        REAL                    XLAM(NI,NJ,NK),XPHI(NI,NJ,NK)
294        REAL                    HLAM(MAXAUF,2),HPHI(MAXAUF,2)
295        REAL                    GWSAVE(4*MAUF+15)
296        REAL      ERAD
297        REAL            SCR,SCI,ACR,ACI,MUSCR,MUSCI,MUACR,MUACI,RT,IT
298       
299        ERAD = 6367470.0
300       
301        GGIND=0
302        DO 4 J = 1,NJ
303         CALL DPLGND(MNAUF,P(0,J),H)
304         DO 3 K = 1,NK
305          LL=0
306          LLP=0
307          LLH=0
308          DO 2 M = 0,MNAUF
309           SCR=0.D0
310           SCI=0.D0
311           ACR=0.D0
312           ACI=0.D0
313           MUSCR=0.D0
314           MUSCI=0.D0
315           MUACR=0.D0
316           MUACI=0.D0
317           LLS=LL
318           LLPS=LLP
319           LLHS=LLH
320           IF(2*M+1.LT.MAUF) THEN
321              DO 1 N = M,MNAUF,2
322               RT=XMN(2*LL,K)
323               IT=XMN(2*LL+1,K)
324               SCR =SCR+ RT*P(LLP,J)
325               SCI =SCI+ IT*P(LLP,J)
326               MUACR =MUACR+RT*H(LLH)
327               MUACI =MUACI+ IT*H(LLH)
328               LL=LL+2
329               LLP=LLP+2
330               LLH=LLH+2
331 1            CONTINUE
332              LL=LLS+1
333              LLP=LLPS+1
334              LLH=LLHS+1
335              DO 11 N = M+1,MNAUF,2
336               RT=XMN(2*LL,K)
337               IT=XMN(2*LL+1,K)
338               ACR =ACR+ RT*P(LLP,J)
339               ACI =ACI+ IT*P(LLP,J)
340               MUSCR =MUSCR+ RT*H(LLH)
341               MUSCI =MUSCI+ IT*H(LLH)
342               LL=LL+2
343               LLP=LLP+2
344               LLH=LLH+2
345 11           CONTINUE
346           ENDIF
347          LL=LLS+(MNAUF-M+1)
348          LLP=LLPS+(MNAUF-M+3)
349          LLH=LLHS+(MNAUF-M+2)
350
351                UFOUC(2*M)=-M*(SCI-ACI)/ERAD
352                UFOUC(2*M+1)=M*(SCR-ACR)/ERAD
353                VFOUC(2*M)=-M*(SCI+ACI)/ERAD
354                VFOUC(2*M+1)=M*(SCR+ACR)/ERAD
355
356                MUFOUC(2*M)=-(MUSCR-MUACR)/ERAD
357                MUFOUC(2*M+1)=-(MUSCI-MUACI)/ERAD
358                MVFOUC(2*M)=-(MUSCR+MUACR)/ERAD
359                MVFOUC(2*M+1)=-(MUSCI+MUACI)/ERAD
3602                       CONTINUE
361                       
362                CALL RFOURTR(VFOUC,
363     *GWSAVE,IFAX,MNAUF,MAUF,1)
364                       
365                CALL RFOURTR(MVFOUC,
366     *GWSAVE,IFAX,MNAUF,MAUF,1)
367                       
368            DO 6 I=0,NI-1
369                  IF(MANF+I.LE. MAUF) THEN
370                    XLAM(I+1,J,K)=VFOUC(MANF+I-1)
371                    XPHI(I+1,J,K)=MVFOUC(MANF+I-1)
372                  ELSE
373                    XLAM(I+1,J,K)=VFOUC(MANF-MAUF+I-1)
374                    XPHI(I+1,J,K)=MVFOUC(MANF-MAUF+I-1)
375                  ENDIF
376    6   CONTINUE
3773               CONTINUE
378                GGIND=GGIND+MAUF
3794       CONTINUE
380
381        RETURN
382        END SUBROUTINE PHGRACUT
383
384C Berechnung der Divergenz aus dem Windfeld (U,V)
385C       im Phasenraum. Zurueckgegeben werden die Felder der
386C       Komponenten des horizontalen Gradienten XLAM,XPHI auf dem Gauss'schen Gitter.
387C       GWSAVE ist ein Hilfsfeld fuer die FFT
388C       P enthaelt die assoziierten Legendrepolynome, H deren Ableitung
389C       MLAT enthaelt die Anzahl der Gitterpunkte pro Breitenkreis
390C       MNAUF gibt die spektrale Aufloesung an, 
391C       NI = Anzahl der Gauss'schen Gitterpunkte,
392C       NJ = Anzahl der Gauss'schen Breiten,
393C       NK = Anzahl der Niveaus
394C Beachte, dass das Windfeld eine um 1 erhoehte Aufloesung in mu-Richtung hat.
395 
396        SUBROUTINE CONTGL(PS,DPSDL,DPSDM,DIV,U,V,BREITE,ETA,
397     *MLAT,A,B,NI,NJ,NK)
398                               
399        IMPLICIT NONE
400       
401        INTEGER NI,NJ,NK,I,J,K,MLAT(NJ),L
402       
403        REAL A(NK+1),B(NK+1)
404        REAL PS(NI),DPSDL(NI),DPSDM(NI)
405        REAL DIV(NI,NK),U(NI,NK),V(NI,NK),ETA(NI,NK)
406        REAL BREITE(NJ)
407       
408        REAL DIVT1,DIVT2,POB,PUN,DPSDT,COSB
409       
410        L=0
411        DO 4 J=1,NJ
412                COSB=(1.0-BREITE(J)*BREITE(J))
413                DO 3 I=1,MLAT(J)
414                        L=L+1
415                        DIVT1=0.0
416                        DIVT2=0.0
417                        DO 1 K=1,NK
418                                POB=A(K)+B(K)*PS(L)
419                                PUN=A(K+1)+B(K+1)*PS(L)
420                               
421                                DIVT1=DIVT1+DIV(L,K)*(PUN-POB)
422                                if(cosb .gt. 0.) then
423                                  DIVT2=DIVT2+(B(K+1)-B(K))*PS(L)*
424     *(U(L,K)*DPSDL(L)+V(L,K)*DPSDM(L))/COSB
425                                endif
426     
427                ETA(L,K)=-DIVT1-DIVT2
4281                       CONTINUE
429
430                        DPSDT=(-DIVT1-DIVT2)/PS(L)
431                       
432                        DO 2 K=1,NK
433                                ETA(L,K)=ETA(L,K)-DPSDT*B(K+1)*PS(L)
4342                       CONTINUE
435                        PS(L)=DPSDT*PS(L)
4363               CONTINUE
4374       CONTINUE
438        RETURN
439        END SUBROUTINE CONTGL
440       
441C       OMEGA berechnet omega im Hybridkoordinatensystem
442C       PS ist der Bodendruck,
443C       DPSDL,DPSDM sind die Komponenten des Gradienten des Logarithmus des
444C       Bodendrucks
445C       DIV,U,V sind die horizontale Divergenz und das horizontale Windfeld
446C       BREITE ist das Feld der Gauss'schen Breiten
447C       E ist omega,
448
449        SUBROUTINE OMEGA(PS,DPSDL,DPSDM,DIV,U,V,BREITE,E,MLAT,A,B,NGI
450     *   ,NGJ,MKK) 
451                               
452        IMPLICIT NONE
453               
454        INTEGER I,J,K,L,NGI,NGJ,MKK,MLAT(NGJ)
455       
456        REAL PS(NGI),DPSDL(NGI),DPSDM(NGI),A(MKK+1),B(MKK+1)
457        REAL DIV(NGI,MKK),U(NGI,MKK),V(NGI,MKK),E(NGI,MKK)
458        REAL BREITE(NGJ)
459       
460        REAL DIVT1,DIVT2,POB,PUN,DP,X,Y,COSB
461        REAL DIVT3(MKK+2)
462       
463        L=0
464        DO 4 J=1,NGJ
465                COSB=(1.0-BREITE(J)*BREITE(J))
466                DO 3 I=1,MLAT(J)
467                        L=L+1
468                        DIVT1=0.0
469                        DIVT2=0.0
470                        DIVT3(1)=0.0
471                        DO 1 K=1,MKK
472                                POB=A(K)+B(K)*PS(L)
473                                PUN=A(K+1)+B(K+1)*PS(L)
474                                DP=PUN-POB
475                               
476                Y=PS(L)*(U(L,K)*DPSDL(L)+V(L,K)*DPSDM(L))/COSB
477                IF(K.LT.3) THEN
478                X=0.0
479                ELSE
480                                X=(B(K+1)-B(K))*Y
481                                ENDIF
482
483                                DIVT1=DIVT1+DIV(L,K)*DP                 
484                                DIVT2=DIVT2+X
485     
486                DIVT3(K+1)=-DIVT1-DIVT2
487               
488                IF(K.GT.1) THEN
489                                E(L,K) = 0.5*(POB+PUN)/DP*Y*
490     *((B(K+1)-B(K))+(A(K+1)*B(K)-A(K)*B(K+1))/
491     *DP*LOG(PUN/POB))
492                                ELSE
493                                        E(L,K) = 0.0
494                                ENDIF
495                                                               
496                        E(L,K) = E(L,K)+0.5*(DIVT3(K)+DIVT3(K+1))
497
4981                       CONTINUE
4993               CONTINUE
5004       CONTINUE
501        RETURN
502        END SUBROUTINE OMEGA
503       
504        END MODULE FTRAFO
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG