Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer, | intent(in) | :: | J | |||
integer, | intent(in) | :: | LFLAG | |||
real(kind=r64), | intent(in) | :: | PDROP | |||
integer, | intent(in) | :: | IL | |||
integer, | intent(in) | :: | N | |||
integer, | intent(in) | :: | M | |||
real(kind=r64), | intent(out) | :: | F(2) | |||
real(kind=r64), | intent(out) | :: | DF(2) | |||
integer, | intent(out) | :: | NF |
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
SUBROUTINE AFEDOP(J,LFLAG,PDROP,IL,N,M,F,DF,NF)
! SUBROUTINE INFORMATION:
! AUTHOR Lixing Gu
! DATE WRITTEN Oct. 2005
! MODIFIED na
! RE-ENGINEERED This subroutine is revised based on a vertical large opening subroutine from COMIS
! PURPOSE OF THIS SUBROUTINE:
! This subroutine simulates airflow and pressure of a detailed large opening component.
! METHODOLOGY EMPLOYED:
! Purpose: This routine calculates the massflow and its derivative
! through a large opening in both flow directions. As input
! the density profiles RhoProfF/T are required aswell as the
! effective pressure difference profile DpProfNew, which is the
! sum of the stack pressure difference profile DpProf and the
! difference of the actual pressures at reference height. The
! profiles are calculated in the routine PresProfile.
! The massflow and its derivative are calculated for each
! interval representing a step of the pressure difference
! profile. The total flow and derivative are obtained by
! summation over the whole opening.
! The calculation is split into different cases representing
! different situations of the opening:
! - closed opening (opening factor = 0): summation of top and
! bottom crack (crack length = lwmax) plus "integration" over
! a vertically distributed crack of length (2*lhmax+lextra).
! - type 1: normal rectangular opening: "integration" over NrInt
! openings with width actlw and height actlh/NrInt
! - type 2: horizontally pivoted window: flow direction assumed
! strictly perpendicular to the plane of the opening
! -> "integration" over normal rectangular openings at top
! and bottom of LO plus a rectangular opening in series with two
! triangular openings in the middle of the LO (most general
! situation). The geometry is defined by the input parameters
! actlw(=lwmax), actlh, axisheight, opening angle.
! Assuming the massflow perpendicular to the opening plane in all
! cases the ownheightfactor has no influence on the massflow.
! REFERENCES:
! Helmut E. Feustel and Alison Rayner-Hooson, "COMIS Fundamentals," LBL-28560,
! Lawrence Berkeley National Laboratory, Berkeley, CA, May 1990
! USE STATEMENTS:
! na
IMPLICIT NONE ! Enforce explicit typing of all variables in this routine
! SUBROUTINE ARGUMENT DEFINITIONS:
INTEGER, INTENT(IN) :: J ! Component number
INTEGER, INTENT(IN) :: LFLAG ! Initialization flag.If = 1, use laminar relationship
REAL(r64), INTENT(IN) :: PDROP ! Total pressure drop across a component (P1 - P2) [Pa]
INTEGER, INTENT(IN) :: IL ! Linkage number
INTEGER, INTENT(IN) :: N ! Node 1 number
INTEGER, INTENT(IN) :: M ! Node 2 number
INTEGER, INTENT(OUT) :: NF ! Number of flows, either 1 or 2
REAL(r64), INTENT(OUT) :: F(2) ! Airflow through the component [kg/s]
REAL(r64), INTENT(OUT) :: DF(2) ! Partial derivative: DF/DP
! SUBROUTINE PARAMETER DEFINITIONS:
REAL(r64), PARAMETER :: RealMax=0.1d+37
REAL(r64), PARAMETER :: RealMin=1d-37
! INTERFACE BLOCK SPECIFICATIONS
! na
! DERIVED TYPE DEFINITIONS
! na
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
INTEGER CompNum
REAL(r64) Width,Height
REAL(r64) fma12 ! massflow in direction "from-to" [kg/s]
REAL(r64) fma21 ! massflow in direction "to-from" [kg/s]
REAL(r64) dp1fma12 ! derivative d fma12 / d Dp [kg/s/Pa]
REAL(r64) dp1fma21 ! derivative d fma21 / d Dp [kg/s/Pa]
REAL(r64) DpProfNew(NrInt+2) ! Differential pressure profile for Large Openings, taking into account fixed
! pressures and actual zone pressures at reference height
REAL(r64) fact ! Actual opening factor
REAL(r64) DifLim ! Limit for the pressure difference where laminarization takes place [Pa]
REAL(r64) Cfact,FvVeloc
REAL(r64) ActLh,ActLw,Lextra,Axishght,ActCD,Cs,expn,Type
REAL(r64) Interval,fmasum,dfmasum,Prefact,EvalHghts(NrInt+2)
REAL(r64) h2,h4,alpha,rholink,c1,c2
REAL(r64) DpZeroOffset,area,WFact,HFact
INTEGER i,Loc,iNum
! FLOW:
! Get component properties
DifLim = 1.0d-4
CompNum = AirflowNetworkCompData(J)%TypeNum
Width = MultizoneSurfaceData(IL)%Width
Height = MultizoneSurfaceData(IL)%Height
fact = MultizoneSurfaceData(IL)%OpenFactor
Loc = (AirflowNetworkLinkageData(IL)%DetOpenNum-1)*(NrInt+2)
iNum = MultizoneCompDetOpeningData(CompNum)%NumFac
ActCD=0.0d0
if (iNum .eq. 2) then
if (fact .le. MultizoneCompDetOpeningData(CompNum)%OpenFac2) then
WFact = MultizoneCompDetOpeningData(CompNum)%WidthFac1+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
(MultizoneCompDetOpeningData(CompNum)%WidthFac2-MultizoneCompDetOpeningData(CompNum)%WidthFac1)
HFact = MultizoneCompDetOpeningData(CompNum)%HeightFac1+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
(MultizoneCompDetOpeningData(CompNum)%HeightFac2-MultizoneCompDetOpeningData(CompNum)%HeightFac1)
CFact = MultizoneCompDetOpeningData(CompNum)%DischCoeff1+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
(MultizoneCompDetOpeningData(CompNum)%DischCoeff2-MultizoneCompDetOpeningData(CompNum)%DischCoeff1)
Else
CALL ShowFatalError('Open Factor is above the maximum input range for opening factors in ' &
//'AirflowNetwork:MultiZone:Component:DetailedOpening = ' &
//TRIM(MultizoneCompDetOpeningData(CompNum)%Name))
End if
end if
if (iNum .eq. 3) then
if (fact .le. MultizoneCompDetOpeningData(CompNum)%OpenFac2) then
WFact = MultizoneCompDetOpeningData(CompNum)%WidthFac1+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
(MultizoneCompDetOpeningData(CompNum)%WidthFac2-MultizoneCompDetOpeningData(CompNum)%WidthFac1)
HFact = MultizoneCompDetOpeningData(CompNum)%HeightFac1+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
(MultizoneCompDetOpeningData(CompNum)%HeightFac2-MultizoneCompDetOpeningData(CompNum)%HeightFac1)
CFact = MultizoneCompDetOpeningData(CompNum)%DischCoeff1+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
(MultizoneCompDetOpeningData(CompNum)%DischCoeff2-MultizoneCompDetOpeningData(CompNum)%DischCoeff1)
else if (fact .le. MultizoneCompDetOpeningData(CompNum)%OpenFac3) then
WFact = MultizoneCompDetOpeningData(CompNum)%WidthFac2+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac2)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac3-MultizoneCompDetOpeningData(CompNum)%OpenFac2)* &
(MultizoneCompDetOpeningData(CompNum)%WidthFac3-MultizoneCompDetOpeningData(CompNum)%WidthFac2)
HFact = MultizoneCompDetOpeningData(CompNum)%HeightFac2+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac2)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac3-MultizoneCompDetOpeningData(CompNum)%OpenFac2)* &
(MultizoneCompDetOpeningData(CompNum)%HeightFac3-MultizoneCompDetOpeningData(CompNum)%HeightFac2)
CFact = MultizoneCompDetOpeningData(CompNum)%DischCoeff2+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac2)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac3-MultizoneCompDetOpeningData(CompNum)%OpenFac2)* &
(MultizoneCompDetOpeningData(CompNum)%DischCoeff3-MultizoneCompDetOpeningData(CompNum)%DischCoeff2)
Else
CALL ShowFatalError('Open Factor is above the maximum input range for opening factors in ' &
//'AirflowNetwork:MultiZone:Component:DetailedOpening = ' &
//TRIM(MultizoneCompDetOpeningData(CompNum)%Name))
end if
end if
if (iNum .eq. 4) then
if (fact .le. MultizoneCompDetOpeningData(CompNum)%OpenFac2) then
WFact = MultizoneCompDetOpeningData(CompNum)%WidthFac1+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
(MultizoneCompDetOpeningData(CompNum)%WidthFac2-MultizoneCompDetOpeningData(CompNum)%WidthFac1)
HFact = MultizoneCompDetOpeningData(CompNum)%HeightFac1+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
(MultizoneCompDetOpeningData(CompNum)%HeightFac2-MultizoneCompDetOpeningData(CompNum)%HeightFac1)
CFact = MultizoneCompDetOpeningData(CompNum)%DischCoeff1+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac1)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac2-MultizoneCompDetOpeningData(CompNum)%OpenFac1)* &
(MultizoneCompDetOpeningData(CompNum)%DischCoeff2-MultizoneCompDetOpeningData(CompNum)%DischCoeff1)
else if (fact .le. MultizoneCompDetOpeningData(CompNum)%OpenFac3) then
WFact = MultizoneCompDetOpeningData(CompNum)%WidthFac2+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac2)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac3-MultizoneCompDetOpeningData(CompNum)%OpenFac2)* &
(MultizoneCompDetOpeningData(CompNum)%WidthFac3-MultizoneCompDetOpeningData(CompNum)%WidthFac2)
HFact = MultizoneCompDetOpeningData(CompNum)%HeightFac2+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac2)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac3-MultizoneCompDetOpeningData(CompNum)%OpenFac2)* &
(MultizoneCompDetOpeningData(CompNum)%HeightFac3-MultizoneCompDetOpeningData(CompNum)%HeightFac2)
CFact = MultizoneCompDetOpeningData(CompNum)%DischCoeff2+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac2)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac3-MultizoneCompDetOpeningData(CompNum)%OpenFac2)* &
(MultizoneCompDetOpeningData(CompNum)%DischCoeff3-MultizoneCompDetOpeningData(CompNum)%DischCoeff2)
else if (fact .le. MultizoneCompDetOpeningData(CompNum)%OpenFac4) then
WFact = MultizoneCompDetOpeningData(CompNum)%WidthFac3+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac3)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac4-MultizoneCompDetOpeningData(CompNum)%OpenFac3)* &
(MultizoneCompDetOpeningData(CompNum)%WidthFac4-MultizoneCompDetOpeningData(CompNum)%WidthFac3)
HFact = MultizoneCompDetOpeningData(CompNum)%HeightFac3+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac3)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac4-MultizoneCompDetOpeningData(CompNum)%OpenFac3)* &
(MultizoneCompDetOpeningData(CompNum)%HeightFac4-MultizoneCompDetOpeningData(CompNum)%HeightFac3)
CFact = MultizoneCompDetOpeningData(CompNum)%DischCoeff3+ &
(fact-MultizoneCompDetOpeningData(CompNum)%OpenFac3)/ &
(MultizoneCompDetOpeningData(CompNum)%OpenFac4-MultizoneCompDetOpeningData(CompNum)%OpenFac3)* &
(MultizoneCompDetOpeningData(CompNum)%DischCoeff4-MultizoneCompDetOpeningData(CompNum)%DischCoeff3)
Else
CALL ShowFatalError('Open Factor is above the maximum input range for opening factors in ' &
//'AirflowNetwork:MultiZone:Component:DetailedOpening = ' &
//TRIM(MultizoneCompDetOpeningData(CompNum)%Name))
end if
end if
! calculate DpProfNew
Do i=1,NrInt+2
DpProfNew(i) = PDROP+DpProf(Loc+i)-DpL(1,iL)
End Do
! Get opening data based on the opening factor
if (fact .eq. 0) then
ActLw=MultizoneSurfaceData(IL)%Width
ActLh=MultizoneSurfaceData(IL)%Height
Cfact = 0.0d0
else
ActLw=MultizoneSurfaceData(IL)%Width*Wfact
ActLh=MultizoneSurfaceData(IL)%Height*HFact
ActCD=CFact
end if
Cs=MultizoneCompDetOpeningData(CompNum)%FlowCoef
Expn=MultizoneCompDetOpeningData(CompNum)%FlowExpo
Type=MultizoneCompDetOpeningData(CompNum)%LVOType
IF (Type.EQ.1) THEN
Lextra=MultizoneCompDetOpeningData(CompNum)%LVOValue
AxisHght=0.0d0
ELSE IF (Type.EQ.2) THEN
Lextra=0.0d0
AxisHght=MultizoneCompDetOpeningData(CompNum)%LVOValue
ActLw=MultizoneSurfaceData(IL)%Width
ActLh=MultizoneSurfaceData(IL)%Height
ENDIF
! Add window multiplier with window close
If (MultizoneSurfaceData(IL)%Multiplier > 1.0d0) Cs=Cs*MultizoneSurfaceData(IL)%Multiplier
! Add window multiplier with window open
If (fact .gt. 0.0d0) then
If (MultizoneSurfaceData(IL)%Multiplier > 1.0d0) ActLw=ActLw*MultizoneSurfaceData(IL)%Multiplier
End If
! Add recurring warnings
If (fact .gt. 0.0d0) then
If (ActLw .eq. 0.0d0) Then
MultizoneCompDetOpeningData(CompNum)%WidthErrCount = MultizoneCompDetOpeningData(CompNum)%WidthErrCount + 1
if (MultizoneCompDetOpeningData(CompNum)%WidthErrCount< 2) then
CALL ShowWarningError('The actual width of the AirflowNetwork:MultiZone:Component:DetailedOpening of '//&
TRIM(MultizoneCompDetOpeningData(CompNum)%Name)//' is 0.')
CALL ShowContinueError('The actual width is set to 1.0E-6 m.')
CALL ShowContinueErrorTimeStamp(' Occurrence info: ')
else
CALL ShowRecurringWarningErrorAtEnd('The actual width of the '// &
'AirflowNetwork:MultiZone:Component:DetailedOpening of ' &
//TRIM(MultizoneCompDetOpeningData(CompNum)%Name)//' is 0 error continues.', &
MultizoneCompDetOpeningData(CompNum)%WidthErrIndex, ActLw,ActLw)
end if
ActLw = 1.0d-6
End If
If (ActLh .eq. 0.0d0) Then
MultizoneCompDetOpeningData(CompNum)%HeightErrCount = MultizoneCompDetOpeningData(CompNum)%HeightErrCount + 1
if (MultizoneCompDetOpeningData(CompNum)%HeightErrCount< 2) then
CALL ShowWarningError('The actual height of the AirflowNetwork:MultiZone:Component:DetailedOpening of '//&
TRIM(MultizoneCompDetOpeningData(CompNum)%Name)//' is 0.')
CALL ShowContinueError('The actual height is set to 1.0E-6 m.')
CALL ShowContinueErrorTimeStamp(' Occurrence info: ')
else
CALL ShowRecurringWarningErrorAtEnd('The actual width of the '// &
'AirflowNetwork:MultiZone:Component:DetailedOpening of ' &
//TRIM(MultizoneCompDetOpeningData(CompNum)%Name)//' is 0 error continues.', &
MultizoneCompDetOpeningData(CompNum)%HeightErrIndex, ActLh,ActLh)
end if
ActLh = 1.0d-6
End If
End If
! Initialization:
NF = 1
Interval=ActLh/NrInt
fma12=0.0d0
fma21=0.0d0
dp1fma12=0.0d0
dp1fma21=0.0d0
! Closed LO
IF (Cfact.EQ.0) THEN
DpZeroOffset=DifLim
! bottom crack
IF (DpProfNew(1).GT.0) THEN
IF (Abs(DpProfNew(1)).LE.DpZeroOffset) THEN
dfmasum=Cs*ActLw*DpZeroOffset**expn/DpZeroOffset
fmasum=DpProfNew(1)*dfmasum
ELSE
fmasum=Cs*ActLw*(DpProfNew(1))**expn
dfmasum=fmasum*expn/DpProfNew(1)
ENDIF
fma12=fma12+fmasum
dp1fma12=dp1fma12+dfmasum
ELSE
IF (Abs(DpProfNew(1)).LE.DpZeroOffset) THEN
dfmasum=-Cs*ActLw*DpZeroOffset**expn/DpZeroOffset
fmasum=DpProfNew(1)*dfmasum
ELSE
fmasum=Cs*ActLw*(-DpProfNew(1))**expn
dfmasum=fmasum*expn/DpProfNew(1)
ENDIF
fma21=fma21+fmasum
dp1fma21=dp1fma21+dfmasum
ENDIF
! top crack
IF (DpProfNew(NrInt+2).GT.0) THEN
IF (Abs(DpProfNew(NrInt+2)).LE.DpZeroOffset) THEN
dfmasum=Cs*ActLw*DpZeroOffset**expn/DpZeroOffset
fmasum=DpProfNew(NrInt+2)*dfmasum
ELSE
fmasum=Cs*ActLw*(DpProfNew(NrInt+2))**expn
dfmasum=fmasum*expn/DpProfNew(NrInt+2)
ENDIF
fma12=fma12+fmasum
dp1fma12=dp1fma12+dfmasum
ELSE
IF (Abs(DpProfNew(NrInt+2)).LE.DpZeroOffset) THEN
dfmasum=-Cs*ActLw*DpZeroOffset**expn/DpZeroOffset
fmasum=DpProfNew(NrInt+2)*dfmasum
ELSE
fmasum=Cs*ActLw*(-DpProfNew(NrInt+2))**expn
dfmasum=fmasum*expn/DpProfNew(NrInt+2)
ENDIF
fma21=fma21+fmasum
dp1fma21=dp1fma21+dfmasum
ENDIF
! side and extra cracks
Prefact=Interval*(2+lextra/ActLh)*Cs
DO i=2,NrInt+1
IF (DpProfNew(i).GT.0) THEN
IF (Abs(DpProfNew(i)).LE.DpZeroOffset) THEN
dfmasum=Prefact*DpZeroOffset**expn/DpZeroOffset
fmasum=DpProfNew(i)*dfmasum
ELSE
fmasum=Prefact*(DpProfNew(i))**expn
dfmasum=fmasum*expn/DpProfNew(i)
ENDIF
fma12=fma12+fmasum
dp1fma12=dp1fma12+dfmasum
ELSE
IF (Abs(DpProfNew(i)).LE.DpZeroOffset) THEN
dfmasum=-Prefact*DpZeroOffset**expn/DpZeroOffset
fmasum=DpProfNew(i)*dfmasum
ELSE
fmasum=Prefact*(-DpProfNew(i))**expn
dfmasum=fmasum*expn/DpProfNew(i)
ENDIF
fma21=fma21+fmasum
dp1fma21=dp1fma21+dfmasum
ENDIF
ENDDO
ENDIF
! Open LO, type 1
IF ((Cfact.NE.0).AND.(type.EQ.1)) THEN
DpZeroOffset=DifLim*1d-3
Prefact=ActLw*ActCd*Interval*SQRT(2.0d0)
DO i=2,NrInt+1
IF (DpProfNew(i).GT.0) THEN
IF (Abs(DpProfNew(i)).LE.DpZeroOffset) THEN
dfmasum=SQRT(RhoProfF(Loc+i)*DpZeroOffset)/DpZeroOffset
fmasum=DpProfNew(i)*dfmasum
ELSE
fmasum=SQRT(RhoProfF(Loc+i)*DpProfNew(i))
dfmasum=0.5d0*fmasum/DpProfNew(i)
ENDIF
fma12=fma12+fmasum
dp1fma12=dp1fma12+dfmasum
ELSE
IF (Abs(DpProfNew(i)).LE.DpZeroOffset) THEN
dfmasum=-SQRT(RhoProfT(Loc+i)*DpZeroOffset)/DpZeroOffset
fmasum=DpProfNew(i)*dfmasum
ELSE
fmasum=SQRT(-RhoProfT(Loc+i)*DpProfNew(i))
dfmasum=0.5d0*fmasum/DpProfNew(i)
ENDIF
fma21=fma21+fmasum
dp1fma21=dp1fma21+dfmasum
ENDIF
ENDDO
fma12=prefact*fma12
fma21=prefact*fma21
dp1fma12=prefact*dp1fma12
dp1fma21=prefact*dp1fma21
ENDIF
! Open LO, type 2
IF ((Cfact.NE.0).AND.(type.EQ.2)) THEN
! Initialization
DpZeroOffset=DifLim*1d-3
! New definition for opening factors for LVO type 2: opening angle = 90 degrees --> opening factor = 1.0
! should be PIOvr2 in below?
alpha=fact*3.14159d0/2.d0
h2=Axishght*(1.0d0-COS(alpha))
h4=Axishght+(ActLh-Axishght)*COS(alpha)
EvalHghts(1)=0.0d0
EvalHghts(NrInt+2)=ActLh
! New definition for opening factors for LVO type 2: pening angle = 90 degrees --> opening factor = 1.0
IF (fact.EQ.1.0d0) THEN
h2=Axishght
h4=Axishght
ENDIF
DO i=2,NrInt+1
EvalHghts(i)=Interval*(i-1.5d0)
ENDDO
! Calculation of massflow and its derivative
DO i=2,NrInt+1
IF (DpProfNew(i).GT.0) THEN
rholink=RhoProfF(Loc+i)
ELSE
rholink=RhoProfT(Loc+i)
ENDIF
IF ((EvalHghts(i).LE.h2).OR.(EvalHghts(i).GE.h4)) THEN
IF (Abs(DpProfNew(i)).LE.DpZeroOffset) THEN
dfmasum=ActCd*ActLw*Interval*SQRT(2.0*rholink* &
DpZeroOffset)/DpZeroOffset*DSIGN(1d0,DpProfNew(i))
fmasum=DpProfNew(i)*dfmasum
ELSE
fmasum=ActCd*ActLw*Interval*SQRT(2.0*rholink*ABS(DpProfNew(i)))
dfmasum=0.5d0*fmasum/DpProfNew(i)
ENDIF
ELSE
! triangular opening at the side of LO
c1=ActCd*ActLw*Interval*SQRT(2.0*rholink)
c2=2*ActCd*ABS((Axishght-Evalhghts(i)))*TAN(alpha)* &
Interval*SQRT(2.0d0*rholink)
IF ((c1.NE.0).AND.(c2.NE.0)) THEN
IF (Abs(DpProfNew(i)).LE.DpZeroOffset) THEN
dfmasum=SQRT(DpZeroOffset/(1/c1/c1+1/c2/c2))/ &
DpZeroOffset*DSIGN(1d0,DpProfNew(i))
fmasum=DpProfNew(i)*dfmasum
ELSE
fmasum=SQRT(ABS(DpProfNew(i))/(1/c1/c1+1/c2/c2))
dfmasum=0.5d0*fmasum/DpProfNew(i)
ENDIF
ELSE
fmasum=0.0d0
dfmasum=0.0d0
ENDIF
ENDIF
IF (DpProfNew(i).GT.0) THEN
fma12=fma12+fmasum
dp1fma12=dp1fma12+dfmasum
ELSE
fma21=fma21+fmasum
dp1fma21=dp1fma21+dfmasum
ENDIF
ENDDO
ENDIF
! Calculate some velocity in the large opening
area=ActLh*ActLw*actCD
if (area.gt.(cs+RealMin)) then
if (area.gt.RealMin) then
FvVeloc=(fma21+fma12)/area
else
FvVeloc=0.0d0
end if
else
! here the average velocity over the full area, may blow half in half out.
! velocity= Fva/Nett area=Fma/Rho/(Cm/( (2**N)* sqrt(1.2) ) )
if (cs.gt.0.0d0) then
! get the average Rho for this closed window
DO i=2,NrInt+1
rholink=0.0d0
IF (DpProfNew(i).GT.0) THEN
rholink=RhoProfF(Loc+i)
ELSE
rholink=RhoProfT(Loc+i)
ENDIF
rholink=rholink/nrInt
rholink=1.2d0
END DO
FvVeloc=(Fma21+fma12)*(2.d0**ExpN)*SQRT(1.2d0)/(rhoLink*Cs)
else
FvVeloc=0.0d0
end if
end if
! Output mass flow rates and associated derivatives
F(1) = fma12-fma21
DF(1) = dp1fma12-dp1fma21
F(2) = 0.0d0
if (fma12 .NE. 0.0d0 .and. fma21 .NE. 0.0d0) then
F(2) = fma21
End if
DF(2) = 0.0d0
RETURN
END SUBROUTINE AFEDOP