	PROGRAM BGMITT

C***********************************************************************
C******   Bilanzgerechte Mittelung inhomogener Strömungsfelder   *******
C******               VDI Richtlinie 4675, Teil 2                *******
C******                    Demo - Programm                       *******
C***********************************************************************
C                                                                      *
C        BEARBEITER:           LUECKING, JADE HOCHSCHULE               *
C        STAND:                MAI 2016                                *
C                                                                      *
C***********************************************************************

      INCLUDE 'mkl.fi'

	PARAMETER (IX=200,JX=100)

	REAL KT,KT1,KTM,KM,KD,KS,KSS,KSI,KSF

      DIMENSION X   (IX,JX),Y (IX,JX),Z (IX,JX)
      DIMENSION U   (IX,JX),V (IX,JX),W (IX,JX)
      DIMENSION P   (IX,JX),T (IX,JX)
	DIMENSION ETAT(IX,JX),KT(IX,JX)

      DIMENSION X1 (3,IX*JX),U1 (3,IX*JX),P1(IX*JX),T1(IX*JX)
      DIMENSION ETAT1(IX*JX),KT1  (IX*JX)

      DIMENSION UM (3,IX*JX),PM  (IX*JX),TM(IX*JX)
      DIMENSION ETATM(IX*JX),KTM (IX*JX)
      DIMENSION KM   (IX*JX),RHOM(IX*JX),HM(IX*JX),SM(IX*JX)

      DIMENSION IQ(4,IX*JX),IT(3,IX*JX*2)
      DIMENSION A(3,3),B(3,3),IPIV(3),WORK(3)
      DIMENSION RCCT(3,3),RHCT(3),RKCT(3),RSCT(3)

	DIMENSION DITP(3)
	DIMENSION UD(3,IX*JX),PD(IX*JX),TD(IX*JX)
      DIMENSION KS(  IX*JX),PS(IX*JX),TS(IX*JX)
      DIMENSION CHIK(IX*JX),XIK(IX*JX),ZETAK(IX*JX)

	DIMENSION UDS(3,IX*JX),PDS(IX*JX),TDS(IX*JX)
      DIMENSION KSS(  IX*JX),PSS(IX*JX),TSS(IX*JX)
      DIMENSION CHS(  IX*JX),XIS(IX*JX),ZES(IX*JX)
      DIMENSION ISUM( IX*JX)
      
      DIMENSION AREA(3),SIP(3),SPS(3),SITP(3),UDF(3)
      DIMENSION AZ3(3),UDZN(3)

      CHARACTER DATEI_B*20,DATEI_C*20,DATEI_D*20,DATEI_E*20,DATEI_X*20
      CHARACTER TEXTM1*40,TEXTM2*40


C **** ABFRAGEN ZU INPUT UND ABLAUF DES PROGRAMMS

	WRITE (*,*) 'Programm BGMITT'

      WRITE (*,*) 'Quereinstieg ueber Datei "D" ja/nein?'
      WRITE (*,*) 'Fuer ja 1, fuer nein 0 eingeben!'
      READ  (*,*) IQUE
      IF(IQUE.NE.0.AND.IQUE.NE.1) STOP 'IQUE'
      
      WRITE (*,*) 'Zylinderkoordinaten ja/nein?'
      WRITE (*,*) 'Fuer ja 1, fuer nein 0 eingeben!'
      READ  (*,*) IZYL
      IF(IZYL.NE.0.AND.IZYL.NE.1) STOP 'IZYL'
     
      IF(IQUE.EQ.0) THEN
          
        WRITE (*,*) 'Datei "B" angeben!'
        READ  (*,*) DATEI_B
	  DATEI_C = DATEI_B
	  DATEI_D = DATEI_B
	  DATEI_E = DATEI_B
        
        WRITE (*,*) 'Format der "B"-Datei angeben!'
        WRITE (*,*) '0 fuer quadrilaterale Elemente   strukturiert'
        WRITE (*,*) '1 fuer quadrilaterale Elemente unstrukturiert'
        WRITE (*,*) '2 fuer trilaterale    Elemente unstrukturiert'
        READ  (*,*) ITYP
        IF(ITYP.NE.0.AND.ITYP.NE.1.AND.ITYP.NE.2) STOP 'ITYP'
        
        WRITE (*,*) 'Bilanzgerechte Kurzzeitmittelung ja/nein?'
        WRITE (*,*) 'Fuer ja 1, fuer nein 0 eingeben!'
        READ  (*,*) IKZM
        IF(IKZM.NE.0.AND.IKZM.NE.1) STOP 'IKZM'
        
      ELSE IF(IQUE.EQ.1) THEN
          
        WRITE (*,*) 'Datei "D" angeben!'
        READ  (*,*) DATEI_D
	  DATEI_E = DATEI_D
	  DATEI_B = DATEI_D
	  DATEI_C = DATEI_D

      ENDIF

	DO N=20,1,-1
	  IF(DATEI_B(N:N).EQ.'.') THEN
	    WRITE (DATEI_B(N-1:N+3),'(A5)') 'B.plt'
	    WRITE (DATEI_C(N-1:N+3),'(A5)') 'C.plt'
	    WRITE (DATEI_D(N-1:N+3),'(A5)') 'D.dat'
	    WRITE (DATEI_E(N-1:N+3),'(A5)') 'E.dat'
	    GOTO 10
	  ENDIF
	ENDDO

   10 CONTINUE


C **** EINLESEN UND AUFBEREITUNG DER DATEI "B" ODER "D"

      IF(IQUE.EQ.0) THEN                      ! Datei "B"

        OPEN  (8,FILE=DATEI_B,STATUS='OLD')

        READ  (8,*)
        READ  (8,*)
        READ  (8,*)
        READ  (8,'(23X,A40)') TEXTM1
        READ  (8,'(23X,A40)') TEXTM2
        READ  (8,'(38X,1E15.0)') RGAS
        READ  (8,'(38X,5E15.0)') A0,A1,A2,A3,A4
        READ  (8,'(38X,4E15.0)') PSTD,TSTD,HSTD,SSTD
        IF(IKZM.EQ.1) READ  (8,'(38X,2E15.0)') PRT,SIGK
        
        WRITE (*,*) 'TEXTM1: ',TEXTM1
        WRITE (*,*) 'TEXTM2: ',TEXTM2

        IF(ITYP.EQ.0) THEN                    ! QUAD strukturiert

          READ  (8,'(12X,2(3X,I5))') NI,NJ

          WRITE (*,*) 'NI,NJ: ',NI,NJ

          IF(NI.GT.IX) THEN
            WRITE (*,*) 'PARAMETER IX ZU KLEIN'
            WRITE (*,*) 'Zum Beenden Taste druecken!'
            READ  (*,*)
            STOP
          ELSE IF(NJ.GT.JX) THEN
            WRITE (*,*) 'PARAMETER JX ZU KLEIN'
            WRITE (*,*) 'Zum Beenden Taste druecken!'
            READ  (*,*)
            STOP
          ENDIF

          IF(IKZM.EQ.0) THEN

	      DO J=1,NJ
	      DO I=1,NI
	        READ (8,*) X   (I,J),Y (I,J),Z (I,J),
     +                   U   (I,J),V (I,J),W (I,J),
     +                   P   (I,J),T (I,J)
            ENDDO
            ENDDO

          ELSE IF(IKZM.EQ.1) THEN

	      DO J=1,NJ
	      DO I=1,NI
	        READ (8,*) X   (I,J),Y (I,J),Z (I,J),
     +                   U   (I,J),V (I,J),W (I,J),
     +                   P   (I,J),T (I,J),
     +                   ETAT(I,J),KT(I,J)
            ENDDO
            ENDDO

          ENDIF
          
          IF(IZYL.EQ.1) THEN                  ! Umrechnung zyl. -> kart.
              
	      DO J=1,NJ
	      DO I=1,NI
              RAD  = Y(I,J)
              PHI  = Z(I,J)
              CRAD = V(I,J)
              CPHI = W(I,J)
              Y(I,J) = RAD *COS(PHI)
              Z(I,J) = RAD *SIN(PHI)
              V(I,J) = CRAD*COS(PHI) - CPHI*SIN(PHI)
              W(I,J) = CRAD*SIN(PHI) + CPHI*COS(PHI)
            ENDDO
            ENDDO

          ENDIF

          DO J=1,NJ                           ! 1D Anordnung
          DO I=1,NI
            NRKN = (J-1)*NI + I
            X1   (1,NRKN) = X   (I,J)
            X1   (2,NRKN) = Y   (I,J)
            X1   (3,NRKN) = Z   (I,J)
            U1   (1,NRKN) = U   (I,J)
            U1   (2,NRKN) = V   (I,J)
            U1   (3,NRKN) = W   (I,J)
            P1   (  NRKN) = P   (I,J)
            T1   (  NRKN) = T   (I,J)
            ETAT1(  NRKN) = ETAT(I,J)
            KT1  (  NRKN) = KT  (I,J)
          ENDDO
          ENDDO
          
          NKNOT = NI*NJ

                                              ! Dreieckszerlegung

          NQUAD = (NI-1)*(NJ-1)               ! QUAD Connectivity
      
          DO J=1,NJ-1
          DO I=1,NI-1
            NRKN = (J-1)*NI + I
            NREL = NRKN - (J-1)
            IQ(1,NREL) = NRKN
            IQ(2,NREL) = NRKN + 1
            IQ(3,NREL) = NRKN + 1 + NI
            IQ(4,NREL) = NRKN + NI
          ENDDO
          ENDDO

          NTRIA = NQUAD * 2                   ! TRIA Connectivity
      
          DO N=1,NQUAD
            NREL = N*2 - 1
            IT(1,NREL) = IQ(1,N)
            IT(2,NREL) = IQ(2,N)
            IT(3,NREL) = IQ(3,N)
            NREL = N*2
            IT(1,NREL) = IQ(3,N)
            IT(2,NREL) = IQ(4,N)
            IT(3,NREL) = IQ(1,N)
          ENDDO

        ELSE IF(ITYP.EQ.1) THEN               ! QUAD unstrukturiert

          READ  (8,'(31X,2(3X,I5))') NKNOT,NQUAD

          WRITE (*,*) 'NKNOT,NQUAD: ',NKNOT,NQUAD

          IF(NKNOT.GT.IX*JX) THEN
            WRITE (*,*) 'PARAMETER IX*JX ZU KLEIN'
            WRITE (*,*) 'Zum Beenden Taste druecken!'
            READ  (*,*)
            STOP
          ELSE IF(NQUAD.GT.IX*JX) THEN
            WRITE (*,*) 'PARAMETER IX*JX ZU KLEIN'
            WRITE (*,*) 'Zum Beenden Taste druecken!'
            READ  (*,*)
            STOP
          ENDIF

          IF(IKZM.EQ.0) THEN

	      DO N=1,NKNOT
	        READ (8,*) X1(1,N),X1(2,N),X1(3,N),
     +                   U1(1,N),U1(2,N),U1(3,N),P1(N),T1(N)
            ENDDO

          ELSE IF(IKZM.EQ.1) THEN

	      DO N=1,NKNOT
	        READ (8,*) X1(1,N),X1(2,N),X1(3,N),
     +                   U1(1,N),U1(2,N),U1(3,N),P1(N),T1(N),
     +                   ETAT1(N),KT1(N)
            ENDDO

          ENDIF
          
          DO N=1,NQUAD
            READ (8,*) IQ(1,N),IQ(2,N),IQ(3,N),IQ(4,N)
          ENDDO
          
          IF(IZYL.EQ.1) THEN                  ! Umrechnung zyl. -> kart.
              
	      DO N=1,NKNOT
              RAD  = X1(2,N)
              PHI  = X1(3,N)
              CRAD = U1(2,N)
              CPHI = U1(3,N)
              X1(2,N) = RAD *COS(PHI)
              X1(3,N) = RAD *SIN(PHI)
              U1(2,N) = CRAD*COS(PHI) - CPHI*SIN(PHI)
              U1(3,N) = CRAD*SIN(PHI) + CPHI*COS(PHI)
            ENDDO

          ENDIF

                                              ! Dreieckszerlegung

          NTRIA = NQUAD * 2                   ! TRIA Connectivity
      
          DO N=1,NQUAD
            NREL = N*2 - 1
            IT(1,NREL) = IQ(1,N)
            IT(2,NREL) = IQ(2,N)
            IT(3,NREL) = IQ(3,N)
            NREL = N*2
            IT(1,NREL) = IQ(3,N)
            IT(2,NREL) = IQ(4,N)
            IT(3,NREL) = IQ(1,N)
          ENDDO

        ELSE IF(ITYP.EQ.2) THEN               ! TRIA unstrukturiert

          READ  (8,'(26X,2(3X,I5))') NKNOT,NTRIA

          WRITE (*,*) 'NKNOT,NTRIA: ',NKNOT,NTRIA

          IF(NKNOT.GT.IX*JX  ) THEN
            WRITE (*,*) 'PARAMETER IX*JX   ZU KLEIN'
            WRITE (*,*) 'Zum Beenden Taste druecken!'
            READ  (*,*)
            STOP
          ELSE IF(NTRIA.GT.IX*JX*2) THEN
            WRITE (*,*) 'PARAMETER IX*JX*2 ZU KLEIN'
            WRITE (*,*) 'Zum Beenden Taste druecken!'
            READ  (*,*)
            STOP
          ENDIF

          IF(IKZM.EQ.0) THEN

	      DO N=1,NKNOT
	        READ (8,*) X1(1,N),X1(2,N),X1(3,N),
     +                   U1(1,N),U1(2,N),U1(3,N),P1(N),T1(N)
            ENDDO

          ELSE IF(IKZM.EQ.1) THEN

	      DO N=1,NKNOT
	        READ (8,*) X1(1,N),X1(2,N),X1(3,N),
     +                   U1(1,N),U1(2,N),U1(3,N),P1(N),T1(N),
     +                   ETAT1(N),KT1(N)
            ENDDO

          ENDIF
          
          DO N=1,NTRIA
            READ (8,*) IT(1,N),IT(2,N),IT(3,N)
          ENDDO
          
          IF(IZYL.EQ.1) THEN                  ! Umrechnung zyl. -> kart.
              
	      DO N=1,NKNOT
              RAD  = X1(2,N)
              PHI  = X1(3,N)
              CRAD = U1(2,N)
              CPHI = U1(3,N)
              X1(2,N) = RAD *COS(PHI)
              X1(3,N) = RAD *SIN(PHI)
              U1(2,N) = CRAD*COS(PHI) - CPHI*SIN(PHI)
              U1(3,N) = CRAD*SIN(PHI) + CPHI*COS(PHI)
            ENDDO

          ENDIF

        ENDIF

      ELSE IF(IQUE.EQ.1) THEN                 ! Datei "D"

        OPEN (8,FILE=DATEI_D,STATUS='OLD')

	  READ (8,*)
	  READ (8,*)
	  READ (8,'(32X,A40)') TEXTM1
	  READ (8,'(32X,A40)') TEXTM2
	  READ (8,*)
	  READ (8,*)
	  READ (8,'(32X,1E15.0)') RGAS
	  READ (8,'(32X,5E15.0)') A0,A1,A2,A3,A4
	  READ (8,*)
	  READ (8,*)
	  READ (8,'(32X,4E15.0)') PSTD,TSTD,HSTD,SSTD
	  READ (8,*)
	  READ (8,'(32X,3E15.0)') AREA(1),AREA(2),AREA(3)
	  READ (8,'(32X,1E15.0)') SMP
	  READ (8,*)
	  READ (8,*)
	  READ (8,*)
	  READ (8,'(32X,1E15.0)') PFM
        SPS(1) = PFM * AREA(1)
        SPS(2) = PFM * AREA(2)
        SPS(3) = PFM * AREA(3)
	  READ (8,*)
	  READ (8,*)
	  READ (8,*)
	  READ (8,'(32X,3E15.0)') SIP(1),SIP(2),SIP(3)
        SIP(1) = SIP(1) * SMP
        SIP(2) = SIP(2) * SMP
        SIP(3) = SIP(3) * SMP
	  READ (8,'(32X,1E15.0)') SHDP
        SHDP = SHDP * SMP
	  READ (8,'(32X,1E15.0)') SKDP
        SKDP = SKDP * SMP
	  READ (8,*)
	  READ (8,*)
	  READ (8,*)
	  READ (8,'(32X,1E15.0)') SHSP
        SHSP = SHSP * SMP
	  READ (8,'(32X,1E15.0)') SKSP
        SKSP = SKSP * SMP
	  READ (8,'(32X,1E15.0)') SSSP
        SSSP = SSSP * SMP
        
        TDM = TFUN(SHDP/SMP,HSTD,A0,A1,A2,TSTD)
        RHOM(1) = PFM / RGAS / TDM
        
        GOTO 500

      ENDIF


C **** FLUSSINTEGRATION ZUR FLÄCHENMITTELUNG

      SMP = 0.0
      DO I=1,3
        AREA(I) = 0.0
        SIP (I) = 0.0
        SPS (I) = 0.0
      ENDDO
      SHDP = 0.0
      SKDP = 0.0
      SHSP = 0.0
      SKSP = 0.0
      SSSP = 0.0

      DO N=1,NTRIA

C **** Metrik: Tangentenvektoren

        DO NZ=1,3
          A(NZ,1) = X1(NZ,IT(2,N)) - X1(NZ,IT(1,N))
          A(NZ,2) = X1(NZ,IT(3,N)) - X1(NZ,IT(1,N))
        ENDDO
        A(1,3) = A(2,1)*A(3,2) - A(3,1)*A(2,2)
        A(2,3) = A(3,1)*A(1,2) - A(1,1)*A(3,2)
        A(3,3) = A(1,1)*A(2,2) - A(2,1)*A(1,2)

C **** Mittelwerte im Element

        UM (1,N) = ( U1 (1,IT(1,N))+U1 (1,IT(2,N))+U1 (1,IT(3,N)) ) / 3.
        UM (2,N) = ( U1 (2,IT(1,N))+U1 (2,IT(2,N))+U1 (2,IT(3,N)) ) / 3.
        UM (3,N) = ( U1 (3,IT(1,N))+U1 (3,IT(2,N))+U1 (3,IT(3,N)) ) / 3.
        PM   (N) = ( P1   (IT(1,N))+P1   (IT(2,N))+P1   (IT(3,N)) ) / 3.
        TM   (N) = ( T1   (IT(1,N))+T1   (IT(2,N))+T1   (IT(3,N)) ) / 3.
        ETATM(N) = ( ETAT1(IT(1,N))+ETAT1(IT(2,N))+ETAT1(IT(3,N)) ) / 3.
        KTM  (N) = ( KT1  (IT(1,N))+KT1  (IT(2,N))+KT1  (IT(3,N)) ) / 3.
        
        KM   (N) = (UM(1,N)**2+UM(2,N)**2+UM(3,N)**2) / 2.
        RHOM (N) = PM(N) / RGAS / TM(N)
        HM   (N) = HFUN(TM(N),TSTD,A0,A1,A2,HSTD)
	  SM   (N) = SFUN(TM(N),PM(N),TSTD,PSTD,A0,A1,A2,RGAS,SSTD)

C **** Massenstrom

	  DMP  = RHOM(N) * (UM(1,N)*A(1,3)+UM(2,N)*A(2,3)+UM(3,N)*A(3,3))


C **** KURZZEITMITTELUNG

        IF(IKZM.EQ.1) THEN

C **** Metrik: Inverse

          DO NZ=1,3
          DO NS=1,3
            B(NZ,NS) = A(NZ,NS)
          ENDDO
          ENDDO

          CALL SGETRF(3,3,B,3,IPIV,INFO)

          CALL SGETRI(3,B,3,IPIV,WORK,3,INFO)

C **** Ableitungen

          DUDX  = (U1 (1,IT(2,N)) - U1 (1,IT(1,N))) * B(1,1) +
     F            (U1 (1,IT(3,N)) - U1 (1,IT(1,N))) * B(2,1)
          DUDY  = (U1 (1,IT(2,N)) - U1 (1,IT(1,N))) * B(1,2) +
     F            (U1 (1,IT(3,N)) - U1 (1,IT(1,N))) * B(2,2)
          DUDZ  = (U1 (1,IT(2,N)) - U1 (1,IT(1,N))) * B(1,3) +
     F            (U1 (1,IT(3,N)) - U1 (1,IT(1,N))) * B(2,3)

          DVDX  = (U1 (2,IT(2,N)) - U1 (2,IT(1,N))) * B(1,1) +
     F            (U1 (2,IT(3,N)) - U1 (2,IT(1,N))) * B(2,1)
          DVDY  = (U1 (2,IT(2,N)) - U1 (2,IT(1,N))) * B(1,2) +
     F            (U1 (2,IT(3,N)) - U1 (2,IT(1,N))) * B(2,2)
          DVDZ  = (U1 (2,IT(2,N)) - U1 (2,IT(1,N))) * B(1,3) +
     F            (U1 (2,IT(3,N)) - U1 (2,IT(1,N))) * B(2,3)

          DWDX  = (U1 (3,IT(2,N)) - U1 (3,IT(1,N))) * B(1,1) +
     F            (U1 (3,IT(3,N)) - U1 (3,IT(1,N))) * B(2,1)
          DWDY  = (U1 (3,IT(2,N)) - U1 (3,IT(1,N))) * B(1,2) +
     F            (U1 (3,IT(3,N)) - U1 (3,IT(1,N))) * B(2,2)
          DWDZ  = (U1 (3,IT(2,N)) - U1 (3,IT(1,N))) * B(1,3) +
     F            (U1 (3,IT(3,N)) - U1 (3,IT(1,N))) * B(2,3)

          DTDX  = (T1 (  IT(2,N)) - T1 (  IT(1,N))) * B(1,1) +
     F            (T1 (  IT(3,N)) - T1 (  IT(1,N))) * B(2,1)
          DTDY  = (T1 (  IT(2,N)) - T1 (  IT(1,N))) * B(1,2) +
     F            (T1 (  IT(3,N)) - T1 (  IT(1,N))) * B(2,2)
          DTDZ  = (T1 (  IT(2,N)) - T1 (  IT(1,N))) * B(1,3) +
     F            (T1 (  IT(3,N)) - T1 (  IT(1,N))) * B(2,3)

          DKTDX = (KT1(  IT(2,N)) - KT1(  IT(1,N))) * B(1,1) +
     F            (KT1(  IT(3,N)) - KT1(  IT(1,N))) * B(2,1)
          DKTDY = (KT1(  IT(2,N)) - KT1(  IT(1,N))) * B(1,2) +
     F            (KT1(  IT(3,N)) - KT1(  IT(1,N))) * B(2,2)
          DKTDZ = (KT1(  IT(2,N)) - KT1(  IT(1,N))) * B(1,3) +
     F            (KT1(  IT(3,N)) - KT1(  IT(1,N))) * B(2,3)

C **** Reynolds-Spannungen etc.

          RCCT(1,1) = -ETATM(N) * (DUDX+DUDX)
	    RCCT(1,2) = -ETATM(N) * (DUDY+DVDX)
	    RCCT(1,3) = -ETATM(N) * (DUDZ+DWDX)

	    RCCT(2,1) = RCCT(1,2)
          RCCT(2,2) = -ETATM(N) * (DVDY+DVDY)
	    RCCT(2,3) = -ETATM(N) * (DVDZ+DWDY)

	    RCCT(3,1) = RCCT(1,3)
	    RCCT(3,2) = RCCT(2,3)
          RCCT(3,3) = -ETATM(N) * (DWDZ+DWDZ)
          
          CPM       = A0 + A1*TM(N) + A2*TM(N)**2

	    RHCT(1  ) = -ETATM(N)*CPM/PRT * DTDX
	    RHCT(2  ) = -ETATM(N)*CPM/PRT * DTDY
	    RHCT(3  ) = -ETATM(N)*CPM/PRT * DTDZ

	    RKCT(1  ) = -ETATM(N)/SIGK * DKTDX
	    RKCT(2  ) = -ETATM(N)/SIGK * DKTDY
	    RKCT(3  ) = -ETATM(N)/SIGK * DKTDZ

	    RSCT(1  ) = -ETATM(N)*(CPM-RGAS)/TM(N) * DTDX
	    RSCT(2  ) = -ETATM(N)*(CPM-RGAS)/TM(N) * DTDY
	    RSCT(3  ) = -ETATM(N)*(CPM-RGAS)/TM(N) * DTDZ

C **** Flüsse

        DO I=1,3
          DITP(I) = UM(I,N)*DMP + RCCT(I,1)*A(1,3) +
     +                            RCCT(I,2)*A(2,3) +
     +                            RCCT(I,3)*A(3,3) + PM(N)*A(I,3)
        ENDDO
	  DHP  = HM(N)*DMP + RHCT(1)*A(1,3)+RHCT(2)*A(2,3)+RHCT(3)*A(3,3)
	  DKP  = KM(N)*DMP + KTM(N)*DMP
     +                   + RKCT(1)*A(1,3)+RKCT(2)*A(2,3)+RKCT(3)*A(3,3)
	  DSP  = SM(N)*DMP + RSCT(1)*A(1,3)+RSCT(2)*A(2,3)+RSCT(3)*A(3,3)
	  DHTP = DHP + DKP

C **** Irreversible Kurzzeit-Mittelwerte, Unterschall-Lösung

	  RHOD = RHOM(N)

	  IZAE = 0
  100   IZAE = IZAE + 1
	  IF(IZAE.GT.1000) THEN
	    WRITE (*,*) 'KZM IZAE,N',IZAE,N
	    STOP
	  ENDIF

	  PD(N)   = ( DITP(1)*A(1,3)+DITP(2)*A(2,3)+DITP(3)*A(3,3)-
     +              DMP**2/RHOD ) / ( A(1,3)**2+A(2,3)**2+A(3,3)**2 )
	  UD(1,N) = ( DITP(1)-PD(N)*A(1,3) ) / DMP
	  UD(2,N) = ( DITP(2)-PD(N)*A(2,3) ) / DMP
	  UD(3,N) = ( DITP(3)-PD(N)*A(3,3) ) / DMP
	  KD      = (UD(1,N)**2+UD(2,N)**2+UD(3,N)**2)/2.
	  HD      = DHTP/DMP - KD
	  TD(N)   = TFUN(HD,HSTD,A0,A1,A2,TSTD)
        SD      = SFUN(TD(N),PD(N),TSTD,PSTD,A0,A1,A2,RGAS,SSTD)

	  RHODA   = RHOD
	  RHOD    = PD(N)/RGAS/TD(N)
	  IF(ABS(RHOD-RHODA)/RHOD.GT.1.E-6) GOTO 100

C **** Reversible Kurzzeit-Mittelwerte

        HS      = DHP/DMP
        KS(N)   = DKP/DMP
        SS      = DSP/DMP
        TS(N)   = TFUN(HS,HSTD,A0,A1,A2,TSTD)
        PS(N)   = PFUN(SS,TS(N),SSTD,TSTD,A0,A1,A2,RGAS,PSTD)

C **** Kenngrößen

	  YPS  = YFUN(HS,HD,SS,SD,TS(N),TD(N),A0,A1,A2)

	  CHIK (N) = YPS / KS(N)
	  XIK  (N) = KD  / KS(N)
	  ZETAK(N) = 1.-CHIK(N)-XIK(N)

        ELSE IF(IKZM.EQ.0) THEN

          RHOD    = RHOM(N)
          PD(N)   = PM(N)
          UD(1,N) = UM(1,N)
          UD(2,N) = UM(2,N)
          UD(3,N) = UM(3,N)
          KD      = KM(N)
          HD      = HM(N)
          TD(N)   = TM(N)
          HS      = HM(N)
          KS(N)   = KM(N)
          SS      = SM(N)
          TS(N)   = TM(N)
          PS(N)   = PM(N)

        ENDIF

C **** Integration der Flüsse etc. zur Flächenmittelung

        SMP = SMP + DMP/2.
        
        IF(IZYL.EQ.0) THEN
          DO I=1,3
            AREA(I) = AREA(I) + A(I,3)/2.
            SIP (I) = SIP (I) + UD(I,N)*DMP/2.
            SPS (I) = SPS (I) + PD(  N)*A(I,3)/2.
          ENDDO
        ELSE IF(IZYL.EQ.1) THEN               ! Umrechnung kart. -> zyl.
          YM = ( X1(2,IT(1,N)) + X1(2,IT(2,N)) + X1(2,IT(3,N)) ) / 3.
          ZM = ( X1(3,IT(1,N)) + X1(3,IT(2,N)) + X1(3,IT(3,N)) ) / 3.
          PHI = ATAN2(ZM,YM)
          AZ3 (1) = A (1,3)
          AZ3 (2) = A (3,3)*SIN(PHI) + A (2,3)*COS(PHI)
          AZ3 (3) = A (3,3)*COS(PHI) - A (2,3)*SIN(PHI)
          UDZN(1) = UD(1,N)
          UDZN(2) = UD(3,N)*SIN(PHI) + UD(2,N)*COS(PHI)
          UDZN(3) = UD(3,N)*COS(PHI) - UD(2,N)*SIN(PHI)
          DO I=1,3
            AREA(I) = AREA(I) + AZ3(I)/2.
            SIP (I) = SIP (I) + UDZN(I)*DMP/2.
            SPS (I) = SPS (I) + PD(  N)*AZ3(I)/2.
          ENDDO
        ENDIF
        
        SHDP = SHDP + HD   *DMP/2.
        SKDP = SKDP + KD   *DMP/2.
        SHSP = SHSP + HS   *DMP/2.
        SKSP = SKSP + KS(N)*DMP/2.
        SSSP = SSSP + SS   *DMP/2.

      ENDDO


C **** AUSGABE DER DATEI "C"

      IF(IKZM.EQ.1) THEN                      ! bei Kurzzeitmittelung

C **** Interpolation der Knotenwerte: Summe der Elementeinflüsse

	DO N=1,NKNOT
        UDS(1,N) = 0.0
        UDS(2,N) = 0.0
        UDS(3,N) = 0.0
        PDS(  N) = 0.0
        TDS(  N) = 0.0
        KSS(  N) = 0.0
        PSS(  N) = 0.0
        TSS(  N) = 0.0
        CHS(  N) = 0.0
        XIS(  N) = 0.0
        ZES(  N) = 0.0
        ISUM( N) = 0
      ENDDO

      DO N=1,NTRIA
        DO I=1,3
          UDS(1,IT(I,N)) = UDS(1,IT(I,N)) + UD(1,N)
          UDS(2,IT(I,N)) = UDS(2,IT(I,N)) + UD(2,N)
          UDS(3,IT(I,N)) = UDS(3,IT(I,N)) + UD(3,N)
          PDS(  IT(I,N)) = PDS(  IT(I,N)) + PD(  N)
          TDS(  IT(I,N)) = TDS(  IT(I,N)) + TD(  N)
          KSS(  IT(I,N)) = KSS(  IT(I,N)) + KS(  N)
          PSS(  IT(I,N)) = PSS(  IT(I,N)) + PS(  N)
          TSS(  IT(I,N)) = TSS(  IT(I,N)) + TS(  N)
          CHS(  IT(I,N)) = CHS(  IT(I,N)) + CHIK (N)
          XIS(  IT(I,N)) = XIS(  IT(I,N)) + XIK  (N)
          ZES(  IT(I,N)) = ZES(  IT(I,N)) + ZETAK(N)
          ISUM( IT(I,N)) = ISUM( IT(I,N)) + 1
        ENDDO
      ENDDO

C **** Ausgabe im TECPLOT-Format "finite-element data"

      OPEN  (10,FILE=DATEI_C,STATUS='UNKNOWN')

	WRITE (10,'(A,A20,A)') 'TITLE = "',DATEI_C,'"'

	WRITE (10,'(3(A))') 'VARIABLES = "X" "Y" "Z"',
     F                    ' "UD" "VD" "WD" "PD" "TD"',
     F                    ' "KS" "PS" "TS" "CHI" "XI" "ZETA"'

	WRITE (10,'(A,A20,A)') 'TEXT X=5 Y=95 H=2 T = "',DATEI_C,'"'
	WRITE (10,'(A,A40,A)') 'TEXT X=5 Y=92 H=2 T = "',TEXTM1 ,'"'
	WRITE (10,'(A,A40,A)') 'TEXT X=5 Y=89 H=2 T = "',TEXTM2 ,'"'

	WRITE (10,1001) 'TEXT X=5 Y=86 H=1 T = "RGAS         = ',RGAS,'"'
	WRITE (10,1002) 'TEXT X=5 Y=84 H=1 T = "CPPOL        = ',
     F                                               A0,A1,A2,A3,A4,'"'
	WRITE (10,1003) 'TEXT X=5 Y=82 H=1 T = "P,T,H,S STD  = ',
     F                                          PSTD,TSTD,HSTD,SSTD,'"'
	WRITE (10,1004) 'TEXT X=5 Y=80 H=1 T = "PRTUR,SIGMAK = ',
     F                                                     PRT,SIGK,'"'
 1001 FORMAT(A,1P1E15.7,A)
 1002 FORMAT(A,1P5E15.7,A)
 1003 FORMAT(A,1P4E15.7,A)
 1004 FORMAT(A,1P2E15.7,A)

	WRITE (10,'(A,2(A,I5))') 'ZONE F=FEPOINT ET=TRIANGLE',
     F                         ' N=',NKNOT,' E=',NTRIA

	DO N=1,NKNOT                            ! arithmetische Mittelwerte
        UDI  = UDS(1,N) / ISUM(N)
        VDI  = UDS(2,N) / ISUM(N)
        WDI  = UDS(3,N) / ISUM(N)
        PDI  = PDS(  N) / ISUM(N)
        TDI  = TDS(  N) / ISUM(N)
        KSI  = KSS(  N) / ISUM(N)
        PSI  = PSS(  N) / ISUM(N)
        TSI  = TSS(  N) / ISUM(N)
        CHI  = CHS(  N) / ISUM(N)
        XII  = XIS(  N) / ISUM(N)
        ZEI  = ZES(  N) / ISUM(N)
        IF(IZYL.EQ.0) THEN
          WRITE (10,1005) X1(1,N),X1(2,N),X1(3,N),
     F                    UDI,VDI,WDI,PDI,TDI,
     F                    KSI,PSI,TSI,CHI,XII,ZEI
        ELSE IF(IZYL.EQ.1) THEN               ! Umrechnung kart. -> zyl.
          RAD = SQRT(X1(2,N)**2+X1(3,N)**2)
          PHI = ATAN2(X1(3,N),X1(2,N))
          CRAD = WDI*SIN(PHI) + VDI*COS(PHI)
          CPHI = WDI*COS(PHI) - VDI*SIN(PHI)
          WRITE (10,1005) X1(1,N),RAD,PHI,
     F                    UDI,CRAD,CPHI,PDI,TDI,
     F                    KSI,PSI,TSI,CHI,XII,ZEI
        ENDIF
      ENDDO
 1005 FORMAT(1P14E15.7)
      
	DO N=1,NTRIA
        WRITE (10,1010) IT(1,N),IT(2,N),IT(3,N)
      ENDDO
 1010 FORMAT(3I5)

      ENDIF


C **** AUSGABE DER DATEI "D"

      OPEN  (12,FILE=DATEI_D,STATUS='UNKNOWN')

	WRITE (12,*) 'Gewichtete Mittelwerte für die bilanzgerechte
     F Flächenmittelung'
	WRITE (12,1200) 'Datei Name:                     ',DATEI_D
	WRITE (12,1201) '1. Überschrift:                 ',TEXTM1
	WRITE (12,1201) '2. Überschrift:                 ',TEXTM2

      WRITE (12,*)
      WRITE (12,*) 'Gas Eigenschaften'
	WRITE (12,1202) 'Spezielle Gaskonstante:         ',RGAS
	WRITE (12,1203) 'cp-Polynom Koeffizienten:       ',A0,A1,A2,A3,A4
      
      WRITE (12,*)
      WRITE (12,*) 'Referenzwerte'
	WRITE (12,1204) 'p,T,h,s ref:                    ',PSTD,TSTD,
     F                                                   HSTD,SSTD

      WRITE (12,*)
	WRITE (12,1205) 'Fläche in Komponenten:          ',AREA(1),
     F                                                   AREA(2),
     F                                                   AREA(3)
	WRITE (12,1202) 'Massenstrom:                    ',SMP
      
      WRITE (12,*)
      WRITE (12,*) 'Flächen gewichtete Mittelwerte'
      WRITE (12,*) 'ggf. bilanzgerecht irreversibel kurzzeitgemittelt'
      PFM = SQRT( (SPS (1)**2+SPS (2)**2+SPS (3)**2) /
     F            (AREA(1)**2+AREA(2)**2+AREA(3)**2) )
	WRITE (12,1202) 'Statischer Druck:               ',PFM
      
      WRITE (12,*)
      WRITE (12,*) 'Massenstrom gewichtete Mittelwerte'
      WRITE (12,*) 'ggf. bilanzgerecht irreversibel kurzzeitgemittelt'
	WRITE (12,1205) 'Geschwindigkeit in Komponenten: ',SIP (1)/SMP,
     F                                                   SIP (2)/SMP,
     F                                                   SIP (3)/SMP
	WRITE (12,1202) 'Spezifische Enthalpie:          ',SHDP   /SMP
	WRITE (12,1202) 'Spezifische kinetische Energie: ',SKDP   /SMP
      
      WRITE (12,*)
      WRITE (12,*) 'Massenstrom gewichtete Mittelwerte'
      WRITE (12,*) 'ggf. bilanzgerecht   reversibel kurzzeitgemittelt'
	WRITE (12,1202) 'Spezifische Enthalpie:          ',SHSP   /SMP
	WRITE (12,1202) 'Spezifische kinetische Energie: ',SKSP   /SMP
	WRITE (12,1202) 'Spezifische Entropie:           ',SSSP   /SMP

 1200 FORMAT(A,A20)
 1201 FORMAT(A,A40)
 1202 FORMAT(A,1P1E15.7)
 1203 FORMAT(A,1P5E15.7)
 1204 FORMAT(A,1P4E15.7)
 1205 FORMAT(A,1P3E15.7)


C **** FLÄCHENMITTELUNG

C **** Irreversible Mittelwerte, Unterschall-Lösung

  500 CONTINUE
     
      DO I=1,3
        SITP(I) = SIP(I) + SPS(I)
      ENDDO
      SHTDP = SHDP + SKDP

      RHOD = RHOM(1)

      IZAE = 0
  200 IZAE = IZAE + 1
      IF(IZAE.GT.1000) THEN
        WRITE (*,*) 'FLM IZAE ',IZAE
	  STOP
      ENDIF

	PDF    = ( SITP(1)*AREA(1)+SITP(2)*AREA(2)+SITP(3)*AREA(3)-
     +           SMP**2/RHOD ) / ( AREA(1)**2+AREA(2)**2+AREA(3)**2 )
      UDF(1) = ( SITP(1)-PDF*AREA(1) ) / SMP
      UDF(2) = ( SITP(2)-PDF*AREA(2) ) / SMP
      UDF(3) = ( SITP(3)-PDF*AREA(3) ) / SMP
      KD     = (UDF(1)**2+UDF(2)**2+UDF(3)**2)/2.
      HD     = SHTDP/SMP - KD
      WRITE (*,*) 'HD',HD
      TDF    = TFUN(HD,HSTD,A0,A1,A2,TSTD)
      SD     = SFUN(TDF,PDF,TSTD,PSTD,A0,A1,A2,RGAS,SSTD)

      RHODA  = RHOD
      RHOD   = PDF/RGAS/TDF
      IF(ABS(RHOD-RHODA)/RHOD.GT.1.E-6) GOTO 200

C **** Reversible Mittelwerte

	HS  = SHSP / SMP
      WRITE (*,*) 'HS',HS
	KSF = SKSP / SMP
	SS  = SSSP / SMP
	TSF = TFUN(HS,HSTD,A0,A1,A2,TSTD)
	PSF = PFUN(SS,TSF,SSTD,TSTD,A0,A1,A2,RGAS,PSTD)

C **** Kenngrößen

	YPS  = YFUN(HS,HD,SS,SD,TSF,TDF,A0,A1,A2)

	CHI  = YPS / KSF
	XI   = KD  / KSF
	ZETA = 1.-CHI-XI


C **** AUSGABE DER DATEI "E"

      OPEN  (14,FILE=DATEI_E,STATUS='UNKNOWN')

	WRITE (14,*) 'Bilanzgerechte Flächenmittelwerte'
	WRITE (14,1400) 'Datei Name:                     ',DATEI_E
	WRITE (14,1401) '1. Überschrift:                 ',TEXTM1
	WRITE (14,1401) '2. Überschrift:                 ',TEXTM2

      WRITE (14,*)
      WRITE (14,*) 'Gas Eigenschaften'
	WRITE (14,1402) 'Spezielle Gaskonstante:         ',RGAS
	WRITE (14,1403) 'cp-Polynom Koeffizienten:       ',A0,A1,A2,A3,A4
      
      WRITE (14,*)
      WRITE (14,*) 'Referenzwerte'
	WRITE (14,1404) 'p,T,h,s ref:                    ',PSTD,TSTD,
     F                                                   HSTD,SSTD

      WRITE (14,*)
      WRITE (14,*) 'Irreversible Mittelwerte'
	WRITE (14,1405) 'Geschwindigkeit in Komponenten: ',UDF (1),
     F                                                   UDF (2),
     F                                                   UDF (3)
	WRITE (14,1402) 'Statischer Druck:               ',PDF
	WRITE (14,1402) 'Temperatur:                     ',TDF

      WRITE (14,*)
      WRITE (14,*) 'Reversible   Mittelwerte'
	WRITE (14,1402) 'Spezifische kinetische Energie: ',KSF
	WRITE (14,1402) 'Statischer Druck:               ',PSF
	WRITE (14,1402) 'Temperatur:                     ',TSF

      WRITE (14,*)
      WRITE (14,*) 'Kennzahlen für die Inhomogenität'
	WRITE (14,1402) 'Druckrückgewinn CHI:            ',CHI
	WRITE (14,1402) 'Kinetische Restenergie XI:      ',XI
	WRITE (14,1402) 'Dissipation ZETA:               ',ZETA

 1400 FORMAT(A,A20)
 1401 FORMAT(A,A40)
 1402 FORMAT(A,1P1E15.7)
 1403 FORMAT(A,1P5E15.7)
 1404 FORMAT(A,1P4E15.7)
 1405 FORMAT(A,1P3E15.7)


      WRITE (*,*) 'Zum Beenden Taste druecken!'
      READ  (*,*)
      
      END

      
      FUNCTION HFUN(T,TSTD,A0,A1,A2,HSTD)
      SCPDT = A0*(T-TSTD) + A1/2.*(T**2-TSTD**2) + A2/3.*(T**3-TSTD**3)
      HFUN  = HSTD + SCPDT
      END
      
      
      FUNCTION SFUN(T,P,TSTD,PSTD,A0,A1,A2,RGAS,SSTD)
      SCPDLNT = A0*LOG(T/TSTD) + A1*(T-TSTD) + A2/2.*(T**2-TSTD**2)
      SFUN    = SSTD + SCPDLNT - RGAS*LOG(P/PSTD)
      END

      
      FUNCTION TFUN(H,HSTD,A0,A1,A2,TSTD)
      TX  = TSTD
      HX  = HSTD
      ISTEU = 0
  100 ISTEU = ISTEU + 1
      IF(ISTEU.GT.1000) THEN
        WRITE (*,*) 'TFUN: H,HX,TX',H,HX,TX
        WRITE (*,*) 'Zum Beenden Taste druecken!'
        READ  (*,*)
        STOP 'TFUN'
      ENDIF
      CPM = A0 + A1*TX + A2*TX**2
      T   = TX + (H-HX)/CPM
      HX  = HFUN(T,TSTD,A0,A1,A2,HSTD)
      IF(ABS((HX-H)/H).GT.1.E-5) THEN
        TX = T
        GOTO 100
      ENDIF
      TFUN  = T
      END

      
      FUNCTION PFUN(S,T,SSTD,TSTD,A0,A1,A2,RGAS,PSTD)
      SCPDLNT = A0*LOG(T/TSTD) + A1*(T-TSTD) + A2/2.*(T**2-TSTD**2)
      PFUN    = PSTD * EXP((SCPDLNT-(S-SSTD))/RGAS)
      END

      
      FUNCTION YFUN(H1,H2,S1,S2,T1,T2,A0,A1,A2)
      SCPDLNT = A0*LOG(T2/T1) + A1*(T2-T1) + A2/2.*(T2**2-T1**2)
      YFUN    = (H2-H1) * (1.-(S2-S1)/SCPDLNT)
      END
