Index: python/pythontest/TestInstallTar/test_untar/src/ftrafo.f
===================================================================
--- python/pythontest/TestInstallTar/test_untar/src/ftrafo.f	(revision 2fb99de9b96e65c8b5caa85c276ad8e0be540ed0)
+++ python/pythontest/TestInstallTar/test_untar/src/ftrafo.f	(revision 2fb99de9b96e65c8b5caa85c276ad8e0be540ed0)
@@ -0,0 +1,504 @@
+	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
