text:
c******************************************************************************
c
c      Manifold valve routine
c
c******************************************************************************
c
       SUBROUTINE FDRPXX(RW,PIN,Y,Z,PC,YC,AD,AV,CFI,CFO,FEGR,TREGR,
     ;                   CD,PD,VD,YD,VVA,FLOW,IERR)
c
c      Determines the flow from a duct through a restriction to a plenum
c      given the duct variables Y,Z and the cylinder pressure PC and entropy
c      parameter Y=P/rho^k. Handles flow in either direction and choking. The
c      restriction can be a valve, a throttle, or junction blockage. The
c      plenum can be the cylinder, the junction, or the atmosphere.
c------------------------------------------------------------------------------
c      This was written as an intake valve program. Therefore, below the
c      restriction is called the "valve" and the plenum the "cylinder".
c
c      The analysis assumes PIN = .TRUE. and changes signs at exit if .FALSE.
c------------------------------------------------------------------------------
c      Model:
c
c      Inflow to the plenum is isentropic to an area = AV*CFI, where the
c      pressure is PC. Outflow from the plenum is isentropic to an area
c      AV*CFO, where the pressure is PV. The pressure recovery is based
c      in a momentum analysis
c
c      If the flow at the restriction is choked, the flow pressure PV
c      can be greater than the back pressure (PC or PD).
c
c      The analysis allows for EGR, which is a set fraction of the flow
c      into the cylinder and enters just upstream of the restriction at a
c      specified temperature. Alteration of the gas specific heats by the
c      recirculated flow is neglected. Specific heats are assumed constant.
c
c      Should be called with AV < AD. If not, the duct area is treated as
c      the valve area.
c
c      With backflow from the valve into the duct, the model assumes that the
c      pressure in the duct is the valve lip pressure (no pressure recovery).
c      The analysis carefully takes into account the entropy discontinuity
c      between that of the massin the duct and that of the gas coming out of
c      the cylinder into the duct.  This is a non-trivial gasdunamic analysis.
c      FDRPT.F is the test program that has been used to test this routine
c      and confirm that it s correct under all circumstances.
c
c      The Z parameter must be corrected for the entropy change between launch
c      of the characteristic and its reception at the valve. This correction
c      depends on the entropy at the valve, which in the case of outlow from
c      the valve is part of the solution. Hence the correction for entropy
c      change along the characteristic is handled as part of the valve flow
c      solution.
c------------------------------------------------------------------------------
c      Input:
c        RW    real work array for the manifold
c        PIN   if .TRUE.,  Z=[2/(k-1)]c+V, V is positive into the plenum
c        PIN   if .FALSE., Z=[2/(k-1)]c-V, V is positive out of the plenum
c        Y     P/rho^k, the entropy parameter, at launch of Z
c        Z     acoustic parameter, m/s
c        PC    cylinder pressure, Pa
c        YC    Y in cylinder
c        AD    duct flow area, m^2
c        AV    valve (restriction) flow area, m^2
c        CFI   valve flow coefficient for flow into cylinder
c        CFO   valve flow coefficient for flow out of cylinder
c        FEGR  fraction of cylinder flow recirculated into the manifold
c        TREGR ratio of stagnation T of the EGR to that of the duct flow
c              (input the value from the last time step)
c
c      Output:
c        CD    sound speed in duct, m/s
c        PD    pressure in duct, Pa
c        VD    velocity in duct, m/s
c        YD    entropy parameter Y at the duct (=Y for inflow)
c        FLOW  mass flow rate into the plenum, kg/s
c        VVA   absolute value of valve velocity, m/s
c
c      IERR = 0 if ok, >0 if problems
c------------------------------------------------------------------------------
       PARAMETER           (NMPTR = 96)
       PARAMETER           (NDRW = 21024)
       DIMENSION           MPTR(NMPTR)
       DIMENSION           RW(NDRW)
       COMMON  /MPTRXX/    MPTR
c------------------------------------------------------------------------------
       EQUIVALENCE (IoK   ,MPTR(1))
       EQUIVALENCE (IoKM1 ,MPTR(2))
       EQUIVALENCE (IoKP1 ,MPTR(3))
       EQUIVALENCE (IoKMKP,MPTR(4))
       EQUIVALENCE (IoKM12,MPTR(6))
       EQUIVALENCE (IoKM3M,MPTR(7))
       EQUIVALENCE (IoKKM1,MPTR(8))
       EQUIVALENCE (IoKPO2,MPTR(10))
       EQUIVALENCE (IoRK  ,MPTR(11))
       EQUIVALENCE (IoR2K ,MPTR(12))
       EQUIVALENCE (IoRKM1,MPTR(13))
       EQUIVALENCE (IoTK  ,MPTR(15))
       EQUIVALENCE (IoTRKM,MPTR(16))
       EQUIVALENCE (IoTRKP,MPTR(17))
       EQUIVALENCE (IoHKPM,MPTR(22))
       EQUIVALENCE (IoHKMK,MPTR(23))
       EQUIVALENCE (IoC0CS,MPTR(24))
       EQUIVALENCE (IoPRPC,MPTR(27))
       EQUIVALENCE (IoCRCC,MPTR(28))
       EQUIVALENCE (IoRCHM,MPTR(29))
       EQUIVALENCE (IoSQTK,MPTR(30))
       EQUIVALENCE (IoK2M1,MPTR(34))
       EQUIVALENCE (IoHKP1,MPTR(35))
       EQUIVALENCE (IoSFMC,MPTR(36))
       EQUIVALENCE (IoCMD2,MPTR(37))
c------------------------------------------------------------------------------
       REAL    MZ,MD,ME,MV,MDX,MMVX,MV2,MD2,ME2,MDC,MX,MXN,MEN,
     ;         MDZC,MDZC2,MVZC,MVZC2,MDVC,MDVC2,MZVC2,
     ;         MEMIN,MEMAX,MVMIN,MVMAX,MDLAST,MVLAST,
     ;         MXMIN,MXMAX
       LOGICAL PIN,MINED,MAXED,VCHOKE,ZCHOKE
c------------------------------------------------------------------------------
c    diagnostic monitor
       LOGICAL OUT
       COMMON /MONMXX/IUMONM,MON,OUT
c------------------------------------------------------------------------------
       DATA  KTRMAX,ERRH,DMOK /20,1.E-5,1E-5/
c------------------------------------------------------------------------------
c     restart point and diagnostic header
2      IF (OUT) WRITE (IUMONM,3)  PIN,Y,Z,PC,YC,YD,
     ;                            AD,AV,CFI,CFO,FEGR,TREGR
3      FORMAT (/' FDRPXX call for:  PIN = ',L2/
     ;   '  Y,Z,PC,YC,YD/AD,AV,CFI,CFO/FEGR,TREGR'/5(1PE14.6)/4E14.6/
     ;      2e14.6)
c
c    check for valve closed
       IF (AV.LE.0) GOTO 80
c
c    check for all FEGR
       IF (FEGR.GE.1)  THEN
           WRITE (IUMONM,4)
4          FORMAT (' Excessive EGR')
           CALL WARNZZ('@','@Excessive EGR; run cancelled.@@')
           IERR = 1
           RETURN
           ENDIF
c
c    determine the flow direction
       BETA = Z/CPYXX(RW,PC,Y)
       GAMMA = RW(IoTRKM)
       IF (OUT) WRITE (IUMONM,8) BETA,GAMMA
8      FORMAT (' BETA,GAMMA',2(1PE14.6))
       IF (BETA.GT.1.00001*GAMMA) THEN
c            inflow
               GOTO 20
           ELSEIF (BETA.LT.0.99999*GAMMA)  THEN
c            outflow
               GOTO 50
           ELSE
c            no flow
               GOTO 80
           ENDIF
c
c ** inflow
c
c    check valve area
20     IF (FEGR.GT.0)  THEN
               PHI = (1 - FEGR) + TREGR*FEGR
               C0VOC0 = SQRT(PHI)
               TERM = AD*C0VOC0/AV
           ELSE
               C0VOC0 = 1
               TERM = AD/AV
               PHI = 1
           ENDIF
       IF (TERM.LT.1) THEN
c            reduce effective valve area so valve chokes first
               AVX = AV*TERM
               TERM = 1
           ELSE
c            valve area ok
               AVX = AV
           ENDIF
       IF (OUT) WRITE (IUMONM,21) AD,AV,AVX
21     FORMAT (' AD,AV,AVX=',3(1PE14.6))
c    get parameters
       CPSI = (1 - FEGR)/PHI
       CPSI2 = CPSI*CPSI
       ZETA = AD/(AVX*CFI)
c
c    set to compute ME for MV = 1 to check for choking
       MV =  1
       TV = RW(IoHKP1)
       STV = RW(IoC0CS)
       MODE = 1
c    trial value based on small ME approximation with MV=1
       ME = RW(IoRCHM)/ZETA
       IF (ME.GT.0.5) ME = 0.5
c    set limits
       MEMAX = 1
       MEMIN = 0
c
c    prepare to iterate
28     MINED = .FALSE.
       MAXED = .FALSE.
       KTR = 0
       DME = 1
c
c - Loop point for NR iterations
c
30     KTR = KTR + 1
       ME2 = ME*ME
       TE = 1 + RW(IoKM12)*ME2
       STE = SQRT(TE)
c
c    compute MD for trial ME
       PSI2 = CPSI2*ME2*TE
       TERM = 2*RW(IoKM1)*PSI2
       IF (TERM.LT.0.01)  THEN
               TERM2 = TERM*TERM
               QN = 0.5*TERM - 0.125*TERM2 + 0.0625*TERM*TERM2
           ELSE
               QN = SQRT(1 + TERM) - 1
           ENDIF
       IF (QN.LT.0)  THEN
           WRITE (IUMONM,31) QN
31         FORMAT (' FDRPXX error; QN =',E14.6)
           GOTO 900
           ENDIF
       MD2 = QN*RW(IoRKM1)
       MD = SQRT(MD2)
       TD = 1 + RW(IoKM12)*MD2
       STD = SQRT(TD)
c
       IF (MODE.EQ.2)  THEN
c        compute MV for trial MD and ME
           TERM = BETA/(MD + GAMMA)
           MV2 = (TERM*TERM*TE - 1)*RW(IoTRKM)
c        if zero, set for no flow
           IF (MV2.LE.0) GOTO 80
           MV = SQRT(MV2)
           TV = 1 + RW(IoKM12)*MV2
           STV = SQRT(TV)
           ENDIF

c    compute H
       H = ZETA*(ME/MV)*(TV/TE)**RW(IoHKPM)
       IF (OUT) WRITE (IUMONM,34) KTR,MD,ME,MV,H
34     FORMAT (' KTR,MD,ME,MV,H',I3,3F10.6,1PE14.6)
       HERR = H - 1
c    check convergence
       IF (((ABS(HERR).LT.ERRH).OR.(ABS(MEMAX-MEMIN).LT.DMOK)).OR.
     ;      (ABS(DME).LT.DMOK))  GOTO 40

c    check count
       IF (KTR.GT.KTRMAX) GOTO 900
c
c    compute H derivatives
       DHDME = H*(1/ME - RW(IoHKPM)*RW(IoKM1)*ME/TE)
       DMDDME = CPSI*(STE + RW(IoKM12)*ME2/STE)/
     ;               (STD + RW(IoKM12)*MD2/STD)
       IF (MODE.EQ.2)  THEN
               DHDMV = H*( - 1/MV + RW(IoHKPM)*RW(IoKM1)*MV/TV)
               DMVDME = (TV/MV)*(- DMDDME*STE/(BETA*STV*RW(IoKM12))
     ;                     + ME/TE)
               DHDMET = DHDME + DHDMV*DMVDME
           ELSE
               DHDMET = DHDME
           ENDIF
c    reset limits
       IF (DHDMET.GT.0)  THEN
               IF (HERR.LT.0)  THEN
                       MEMIN = ME
                       HMIN = H
                       MINED = .TRUE.
                   ELSE
                       MEMAX = ME
                       HMAX = H
                       MAXED = .TRUE.
                   ENDIF
           ELSE
               IF (HERR.LT.0) THEN
                       MEMAX = ME
                       HMAX = H
                       MAXED = .TRUE.
                   ELSE
                       MEMIN = ME
                       HMIN = H
                       MINED = .TRUE.
                   ENDIF
           ENDIF
c    compute change
       DME = - HERR/DHDMET
       MEN = ME + DME
c    limit change
       IF (MAXED.AND.MINED)  THEN
c            solution bracketed
               IF (ABS(DME).GT.0.5*(MEMAX - MEMIN)) THEN
c                oscillating; split bracket
                   MEN = 0.5*(MEMAX + MEMIN)
                   ENDIF
           ELSE
c            not yet bracketed; limit change
               IF (ABS(DME).GT.0.2) THEN
                   DME = 0.2*DME/ABS(DME)
                   MEN = ME + DME
                   ENDIF
c            stay in bounds
               IF (MEN.GE.MEMAX)  THEN
                       MEN = 0.5*(ME + MEMAX)
                   ELSEIF (MEN.LE.MEMIN)  THEN
                       MEN = 0.5*(ME + MEMIN)
                   ENDIF
           ENDIF
       ME = MEN
       GOTO 30
c
c    H converged
40     IF (MODE.EQ.1)  THEN
c        check for choking
           BETAC = (MD + GAMMA)*STV/STE
           IF (BETA.LT.BETAC) THEN
c            not choked; set to find solution
               MEMAX = ME
               MEMIN = 0
               MODE = 2
c            trial value based on small M approximation
               BZ2 = BETA*ZETA*ZETA
               TERM = (BETA - GAMMA)*BZ2*RW(IoKM1)/CPSI2
               IF (TERM.LT.0.001)  THEN
                       QN = CPSI*0.5*TERM
                   ELSE
                       QN = CPSI*(SQRT(1+TERM) - 1)
                   ENDIF
               ME = QN/(RW(IoKM12)*BZ2)
               IF (ME.GT.0.5) ME = 0.5
               GOTO 28
               ENDIF
           ENDIF
c
c    converged on solution
       CD = Z/(MD + GAMMA)
       VD = MD*CD
c
c    final inflow calculations
       YD = Y
       PD = PCYXX(RW,CD,YD)
       RHOD = RW(IoK)*PD/(CD*CD)
       FLOW = AD*RHOD*VD/(1 - FEGR)
       VVA = MV*CD*C0VOC0*SQRT(TD/TV)
       GOTO 90
c
c ** outflow; EGR is assumed to be outflow when the main flow is outflow
c
c      State V is at the valve area AV
c      State D is in the duct between the valve and entropy discontinuity
c      State Z is in the duct where Z and Y are specified.
c
c    set effective valve area if duct smaller
50     AVX = MIN(AD,AV)
       ALPHA = AD/((1 - FEGR)*CFO*AVX)
       IF (OUT) WRITE (IUMONM,51) ALPHA,AD,AV,AVX
51     FORMAT (' ALPHA,AD,AV,AVX=',4(1PE14.6))
c
c    set stagnation sound speed
       C0 = CPYXX(RW,PC,YC)
c    calculate constants
       YCOYZP = (YC/Y)**RW(IoR2K)
       XT = Z/C0
       XT2 = XT*XT
       ALPHA2 = ALPHA*ALPHA
       B1 = ALPHA
       B2 = RW(IoK)
       B3= RW(IoK)/(1 - FEGR)
       C1 = RW(IoTRKM)/(YCOYZP*ALPHA**RW(IoHKMK))
c    assume Z not choled
       ZCHOKE = .FALSE.
c
c *  re-entry point if Z is found to be choked
c
c    set to compute MD for MV = -1 to check for choking
58     MV = - 1
       KTR  = - 1
       VCHOKE = .FALSE.
c    set limits
       MXMAX = 0
       MXMIN = - 1
c    prepare to iterate
       MINED = .FALSE.
       MAXED = .FALSE.
       DMX = 1
c
c  - loop point for NR iteration
c
60     KTR = KTR + 1
c
       IF (.NOT.VCHOKE)  THEN
c        calculate MD for trial MV from the momentum balance with PB = PV
           MV2 = MV*MV
           TV = 1 + RW(IoKM12)*MV2
           STV = SQRT(TV)
           PHI = (B1 + B3*MV2)/(MV*STV)
           PHI2 = PHI*PHI
           EPS = PHI2*RW(IoHKMK)*RW(IoRK) - 1
c        quadratic equation solution
           IF (ABS(EPS).LT.0.05) THEN
c                series expansion near the singularity
                   MD2 = RW(IoCMD2+1)
                   EPSN = 1
                   DO 61 N=2,6
                       EPSN = EPSN*EPS
                       MD2 = MD2 + RW(IoCMD2+N)*EPSN
61                     CONTINUE
                   TERM = - 1
               ELSE
c                full evaluation with accurate square root
                   RPHI2 = 1/PHI2
                   TERM =2*RW(IoKP1)*RPHI2
                   IF (TERM.GT.1.0) THEN
                       IF (TERM.GT.1.00001)  THEN
                           IF (OUT) WRITE (IUMONM,62) MV,TERM
62                         FORMAT (' Error: MV,TERM',2(1PE14.6))
                           IERR = 1
                           GOTO 900
                           ENDIF
                       TERM = 1
                       ENDIF
                   IF (TERM.LT.0.1) THEN
                           TERM2 = TERM*TERM
                           QN = RW(IoTK)*RPHI2
     ;                      - 0.5*TERM - 0.125*TERM2 - 0.0625*TERM2*TERM
     ;                      -  0.03906*TERM2*TERM2
                       ELSE
                           QN = RW(IoTK)*RPHI2 - 1 + SQRT(1 - TERM)

                       ENDIF
                   MD2 = QN/(2*(RW(IoKM12) - B2*B2*RPHI2))
               ENDIF
           IF (MD2.LE.0) THEN
               IF (OUT) WRITE (IUMONM,63) MV,MD2
63             FORMAT (' Error: MV,MD2',2(1PE14.6))
               IERR = 1
               GOTO 900
               ENDIF
           MD = - SQRT(MD2)
           TD = 1 + RW(IoKM12)*MD2
           STD = SQRT(TD)
           ENDIF
c
c    compute test function
       IF (ZCHOKE)  THEN
c            from MZ= -1
               H = - MD*YCOYZP*((ALPHA*MD/MV)**RW(IoHKMK))
     ;                        *(STV/STD)**RW(IoKPO2)
           ELSE
c            from the acoustic equation
               H1 = XT*STD/MD
               H2 = -C1*((MV/MD)**RW(IoHKMK))*((STD/STV)**RW(IoKPO2))/MD
               H = H1 + H2
           ENDIF
       HERR = H - 1
c    diagnostic output
       IF (OUT) WRITE (IUMONM,65) KTR,MD,MV,H
65     FORMAT (' KTR,MD,MV,H=',I3,3(1PE14.6))
c
c    branch on task
       IF (KTR.EQ.0)  THEN
c        calculated MD for MV=-1
           IF ( ((.NOT.ZCHOKE).AND.(H.GE.1)).OR.
     ;           (ZCHOKE.AND.(H.LE.1)) ) THEN
c            choked at V
               VCHOKE = .TRUE.
               IF (MD2.GE.1)  THEN
c                also choked at D; revise MD (PB<PV)
                   MD = - 1
                   TD = RW(IoHKP1)
                   STD = RW(IoC0CS)
                   GOTO 68
                   ENDIF
c            V choked but D not choked; go iterate for MD at MV=-1
               GOTO 66
               ENDIF
c        V not choked; set trial MV
           IF (ZCHOKE)  THEN
                   MV = - ALPHA/YCOYZP
               ELSE
                   MV = ALPHA*(XT - RW(IoTRKM)/YCOYZP)
               ENDIF
           IF (MV.LT.-0.5) MV = - 0.5
           GOTO 60
           ENDIF
c
c    check convergence of H
66     IF (((ABS(HERR).LT.ERRH).OR.(ABS(MXMAX-MXMIN).LT.DMOK)).OR.
     ;   (ABS(DMX).LT.DMOK))  THEN
           IF (ABS(MV).LT.1E-5)  GOTO 80
           IF (ZCHOKE)  GOTO 70
c        converged with Z not choked
           GOTO 68
           ENDIF
c
c    check iteration count
       IF (KTR.GT.KTRMAX)  THEN
           IF (OUT)  WRITE (IUMONM,67)
67         FORMAT (' FDRPXX outflow iteration not convergent')
           GOTO 900
           ENDIF
c
c    N-R iteration
       IF (.NOT.VCHOKE)  THEN
c        calculate dMD/dMV from the momentum balance
           U1 = B1 + B3*MV2
           W1 = U1*STD
           U2 = B2*MD + 1/MD
           W2 = U2*MV*STV
           CDMV = W1*2*B3*MV/U1 - W2*(1/MV + RW(IoKM12)*MV/TV)
           CDMD = W1*RW(IoKM12)*MD/TD - W2*(B2 - 1/MD2)/U2
           DMDDMV = - CDMV/CDMD
           ENDIF
       IF (ZCHOKE)  THEN
c            iterating to find MD and MV for MZ = -1 (no acoustic condition)
               DHDMD = H*(1/MD + RW(IoHKMK)/MD
     ;                     - RW(IoKPO2)*RW(IoKM12)*MD/TD)
               IF (.NOT.VCHOKE)  THEN
                   DHDMV = H*( - RW(IoHKMK)/MV
     ;                     + RW(IoKPO2)*RW(IoKM12)*MV/TV)
                   ENDIF
           ELSE
c            iterating to find MD and MV to satisfy the acoustic condition
               DHDMD = H1*(RW(IoKM12)*MD/TD - 1/MD)
     ;           + H2*((- RW(IoHKMK) - 1)/MD
     ;                + RW(IoKPO2)*RW(IoKM12)*MD/TD)
               IF (.NOT.VCHOKE) THEN
                   DHDMV = H2*(RW(IoHKMK)/MV
     ;                    - RW(IoKPO2)*RW(IoKM12)*MV/TV)
                   ENDIF
           ENDIF
       HERR = H - 1
c    branch on mode
       IF (VCHOKE)  THEN
c            choked at V; adjusting MD keeping MV=-1
               MX = MD
               DHDMX = DHDMD
           ELSE
c            not choked at V; adjusting MV
               MX = MV
               DHDMX = DHDMV +  DHDMD*DMDDMV
           ENDIF
c    reset limits
       IF (DHDMX.GT.0)  THEN
               IF (HERR.LT.0)  THEN
                       MXMIN = MX
                       HMIN = H
                       MINED = .TRUE.
                   ELSE
                       MXMAX = MX
                       HMAX = H
                       MAXED = .TRUE.
                   ENDIF
           ELSE
               IF (HERR.LT.0) THEN
                       MXMAX = MX
                       HMAX = H
                       MAXED = .TRUE.
                   ELSE
                       MXMIN = MX
                       HMIN = H
                       MINED = .TRUE.
                   ENDIF
           ENDIF
c    compute change
       DMX = - HERR/DHDMX
       MXN = MX + DMX
c    limit change
       IF (MAXED.AND.MINED)  THEN
c            solution bracketed
               IF (ABS(DMX).GT.0.5*(MXMAX - MXMIN)) THEN
c                oscillating; split bracket
                   MXN = 0.5*(MXMAX + MXMIN)
                   ENDIF
           ELSE
c            not yet bracketed; limit change
               IF (ABS(DMX).GT.0.2) THEN
                   DMX = 0.2*DMX/ABS(DMX)
                   MXN = MX + DMX
                   ENDIF
c            stay in bounds
               IF (MXN.GE.MXMAX)  THEN
                       MXN = 0.5*(MX + MXMAX)
                   ELSEIF (MXN.LE.MXMIN)  THEN
                       MXN = 0.5*(MX + MXMIN)
                   ENDIF
           ENDIF
c    load the new value
       IF (VCHOKE)  THEN
c            set new trial MD
               MD = MXN
           ELSE
c            set new trial MV
               MV = MXN
           ENDIF
       GOTO 60

c    MD and MV are ok; check Z
68     MZ = MD*YCOYZP*((ALPHA*MD/MV)**RW(IoHKMK))*(STV/STD)**RW(IoKPO2)
       IF (MZ.LT.-1)  THEN
c        choked at Z; set to find MV and MD for MZ=-1
           ZCHOKE = .TRUE.
           GOTO 58
           ENDIF
c
c    final outflow calculations
70     CD = C0/STD
       VD = MD*CD
       CD2 = CD*CD
       PV = PC/TV**RW(IoKKM1)
       PD = PV*STV*MV/(ALPHA*STD*MD)
       YD = YC2PXX(RW,CD2,PD)
       RHOD = RW(IoK)*PD/CD2
       FLOW = AD*RHOD*VD/(1 - FEGR)
       VVA = ABS(MV)*C0/STV
       GOTO 90
c
c    no flow
80     VD = 0
       MD = 0
       MV = 0
       CD = Z*RW(IoKM12)
       YD = Y
       PD = PCYXX(RW,CD,YD)
       FLOW = 0
       VVA = 0
c
c    normal exit
90     IERR = 0

c    check flow direction
       IF (.NOT.PIN)  THEN
c        set postive out
           VD = - VD
           FLOW = - FLOW
           MD = - MD
           MV = - MV
           ENDIF
       RETURN
c
c *  error exit
c
c    check for rerun
900    IF (.NOT.OUT)  THEN
c        rerun with diagnostic output
           OUT = .TRUE.
           GOTO 2
           ENDIF
c
c    write the restart file
       OPEN (10,FILE='FINDFLOW.BAD',FORM='UNFORMATTED')
       WRITE (10) PIN,Y,Z,PC,YC,AD,AV,CFI,CFO,FEGR,TREGR,
     ;       (RW(I),I=1,70)
       CLOSE (10)
c    complete the monitor file and repoen a dummy
       WRITE (IUMONM,902)
902    FORMAT (' Test file FINDFLOW.BAD written')
       CLOSE (IUMONM)
       OPEN (IUMONM,FILE='MANTEMP.MON')
c    instruct user
       CALL WARNZZ('#','#The manifold entry/exit flow routine '//
     ;   'failed to converge.  Files MANIFOLD.MON and FINDFLOW.BAD '//
     ;   'contain data needed to diagnose and fix the problem.  '//
     ;   'If your input data seem ok, please quit ESP now and save '//
     ;   'these files and your setup file '//
     ;   'in a separate directory (so that they will '//
     ;   'not be overwritten by your next ESP run) and contact# #'//
     ;  '    wcr@thermo.stanford.edu# #'//
     ;  'to arrange for transfer of these files so that we can '//
     ;  'make ESP work better.  Do not email these files!##')
       IERR = 1
       RETURN
c
       END
c******************************************************************************
c
       SUBROUTINE CENDXX(RW,Y,Z,CD,PD,VD,YD)
c
c      Determines the duct state at a closed end.
c
c      If Z=[2/(k-1)]c+V then V is positive into the closure.
c      If Z=[2/(k-1)]c-V then V is positive out of the closure.
c------------------------------------------------------------------------------
c      Input:
c        RW    real work array for the manifold
c        Y     P/rho^k, the entropy parameter at the duct end
c        Z     acoustic parameter, m/s
c
c      Output:
c        CD    sound speed in duct, m/s
c        PD    pressure in duct, Pa
c        VD    velocity in duct, m/s  (zero at closed end)
c        YD    entropy parameter in the duct at the closed end
c------------------------------------------------------------------------------
       PARAMETER           (NMPTR = 96)
       PARAMETER           (NDRW = 21024)
       DIMENSION           MPTR(NMPTR)
       DIMENSION           RW(NDRW)
       COMMON  /MPTRXX/    MPTR
c------------------------------------------------------------------------------
       EQUIVALENCE (IoKM12,MPTR(6))
       EQUIVALENCE (IoR2K ,MPTR(12))
c------------------------------------------------------------------------------
       VD = 0
       CD = Z*RW(IoKM12)
       PD = PCYXX(RW,CD,YD)
       YD = Y
       RETURN
       END
c******************************************************************************

max_stars_repo_name: kjfergu/esp

max_issues_repo_name: fergu/esp

max_forks_repo_name: fergu/esp


text:
      subroutine wphase (nph, em, eref, lmax, ne, ph)

c     Writes phase data to file PHASExx.DAT for each shell

      implicit double precision (a-h, o-z)

      include 'dim.h'

      complex*16 eref(nex)
      complex*16 ph(nex,ltot+1,0:nphx)
      dimension em(nex)
      dimension lmax(0:nphx)
      character*30  fname

c     Dump phase data, eref and complex phase for each shell
      do 260  iph = 0, nph
c        prepare file for shell's phase data
         write(fname,242)  iph
  242    format('phase', i2.2, '.dat')
         open (unit=1, file=fname, status='unknown', iostat=ios)
         call chopen (ios, fname, 'wphase')
         call wthead  (1)
c        write out unique pot and lmax
         write(1,244)   iph, lmax(iph), ne
  244    format (1x, 3i4, '   unique pot,  lmax, ne ')
c        for each energy
c        ie, em, eref, p=sqrt(em-eref)
c        ph array to ltot+1, 5 values per line
         do 250  ie = 1, ne
            xp = sqrt(em(ie) - eref(ie))
            write(1,246)  ie, em(ie), eref(ie), sqrt(em(ie)-eref(ie))
  246       format ('   ie        energy      re(eref)',
     1              '      im(eref)',
     2              '         re(p)         im(p)', /,
     3              1x, i4, 1p, 5e14.6)
            write(1,248)  (ph(ie,ll,iph), ll=1,lmax(iph)+1)
  248       format (1x, 1p, 4e14.6)
  250    continue
         close(unit=1)
  260 continue

      return
      end

ext: f

max_stars_repo_name: keechul/ifeffit

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: bruceravel/ifeffit

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: bruceravel/ifeffit

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
program main
    integer, parameter :: k4 = 4

    ! equivalent: STyCharacter (CharLenInt 3)
    character(3)              :: s11 = "sup"
    character(len=3)          :: s12 = "sup"
    character(len=3 , kind=1) :: s13 = "sup"
    character(kind=1, len=3 ) :: s14 = "sup"

    ! equivalent: STyCharacter (CharLenInt 1)
    character                 :: s21 = "c"
    character(1)              :: s22 = "c"
    character(len=1)          :: s23 = "c"
    character(kind=1)         :: s24 = "c"
    character(len=1 , kind=1) :: s25 = "c"
    character(kind=1, len=1 ) :: s26 = "c"

    ! equivalent: STyCharacter (CharLenInt 3)
    character    :: s31*3 = "sup"
    character(3) :: s32   = "sup"
    character(3) :: s33*3 = "sup"   ! FS soft type error (...using decl length)
    character(2) :: s34*3 = "sup"   ! FS soft type error (...using decl length)
    character(9) :: s35*3 = "sup"   ! FS soft type error (...using decl length)

    ! type error: non-CHARACTER given length
    integer      :: i11*4 = 1234

    ! special RHS len expr
    character    :: s41*(*)   = "buh?"

    ! unimplemented: unsupported RHS len expr
    character    :: s43*k4    = "uhhh"
    character    :: s44*(2*2) = "huh."
end program main

max_stars_repo_name: envp/fortran-src

max_issues_repo_name: envp/fortran-src

max_forks_repo_name: envp/fortran-src


text:
The Arroyo Neighborhood Association is one of several Neighborhood Associations in Davis.

The Arroyo Neigborhood Association includes the Evergreen Forest, West Park and Evergreen Meadows South or #4 developments. Going clockwise it is roughly bounded by and including (where appropriate) Santa Rosa, Crystal Grove, Hampton, Shasta, El Capitan, White Sands, Pinnacles, and back out to Shasta.

the above was taken from the website

They even have a http://dcn.davis.ca.us/~arroyo/index.shtml website

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! { dg-do run }

program asyncwait
  integer, parameter :: N = 64
  real, allocatable :: a(:), b(:), c(:)
  integer i

  allocate (a(N))
  allocate (b(N))
  allocate (c(N))

  !$acc parallel async (0)
  !$acc loop
  do i = 1, N
    a(i) = 1
  end do
  !$acc end parallel

  !$acc parallel async (1)
  !$acc loop
  do i = 1, N
    b(i) = 1
  end do
  !$acc end parallel

  !$acc wait (0, 1)

  !$acc parallel
  !$acc loop
  do i = 1, N
    c(i) = a(i) + b(i)
  end do
  !$acc end parallel

  do i = 1, N
    if (c(i) .ne. 2.0) call abort
  end do

  !$acc kernels async (0)
  !$acc loop
  do i = 1, N
    a(i) = 1
  end do
  !$acc end kernels

  !$acc kernels async (1)
  !$acc loop
  do i = 1, N
    b(i) = 1
  end do
  !$acc end kernels

  !$acc wait (0, 1)

  !$acc kernels
  !$acc loop
  do i = 1, N
    c(i) = a(i) + b(i)
  end do
  !$acc end kernels

  do i = 1, N
    if (c(i) .ne. 2.0) call abort
  end do

  deallocate (a)
  deallocate (b)
  deallocate (c)
end program asyncwait


subroutine init_time
  use amr_commons
  use hydro_commons
  use pm_commons
  use cooling_module
#ifdef grackle
  use grackle_parameters
#endif
#ifdef RT
  use rt_cooling_module
#endif
  use mpi_mod
  implicit none
  integer::i,Nmodel
  real(kind=8)::T2_sim
#ifdef grackle
  real(dp)::scale_nH,scale_T2,scale_l,scale_d,scale_t,scale_v
  logical::file_exists
#ifndef WITHOUTMPI
  integer::info
#endif
#endif

  if(nrestart==0)then
     if(cosmo)then
        ! Get cosmological parameters from input files
        call init_cosmo
     else
        ! Get parameters from input files
        if(initfile(levelmin).ne.' '.and.filetype.eq.'grafic')then
           call init_file
        endif
        t=0
        aexp=1
     end if
  end if

  if(cosmo)then

     ! Allocate look-up tables
     n_frw=1000
     allocate(aexp_frw(0:n_frw),hexp_frw(0:n_frw))
     allocate(tau_frw(0:n_frw),t_frw(0:n_frw))

     ! Compute Friedman model look up table
     if(myid==1)write(*,*)'Computing Friedman model'
     call friedman(dble(omega_m),dble(omega_l),dble(omega_k), &
          & 1d-6,dble(aexp_ini), &
          & aexp_frw,hexp_frw,tau_frw,t_frw,n_frw)

     ! Compute initial conformal time
     ! Find neighboring expansion factors
     i=1
     do while(aexp_frw(i)>aexp.and.i<n_frw)
        i=i+1
     end do
     ! Interploate time
     if(nrestart==0)then
        t=tau_frw(i)*(aexp-aexp_frw(i-1))/(aexp_frw(i)-aexp_frw(i-1))+ &
             & tau_frw(i-1)*(aexp-aexp_frw(i))/(aexp_frw(i-1)-aexp_frw(i))
        aexp=aexp_frw(i)*(t-tau_frw(i-1))/(tau_frw(i)-tau_frw(i-1))+ &
             & aexp_frw(i-1)*(t-tau_frw(i))/(tau_frw(i-1)-tau_frw(i))
        hexp=hexp_frw(i)*(t-tau_frw(i-1))/(tau_frw(i)-tau_frw(i-1))+ &
             & hexp_frw(i-1)*(t-tau_frw(i))/(tau_frw(i-1)-tau_frw(i))
     end if
     texp=t_frw(i)*(t-tau_frw(i-1))/(tau_frw(i)-tau_frw(i-1))+ &
          & t_frw(i-1)*(t-tau_frw(i))/(tau_frw(i-1)-tau_frw(i))
  else
     texp=t
  end if

  ! Initialize cooling model
#ifdef grackle
  if(use_grackle==1)then
     if(myid==1)then
        write(*,'(A50)')"__________________________________________________"
        write(*,*)'Grackle - Computing cooling model'
        write(*,*)'Grackle - Loading ',TRIM(grackle_data_file)
        write(*,'(A50)')"__________________________________________________"
     endif
     INQUIRE(FILE=grackle_data_file,EXIST=file_exists)
     if(.not.file_exists) then
        if(myid==1) write(*,*) grackle_data_file," not found"
        call clean_stop
     endif

     iresult = set_default_chemistry_parameters(my_grackle_data)
     if(iresult.eq.0)then
         write(*,*) 'Grackle - error in initialize_chemistry_data'
#ifndef WITHOUTMPI
         call MPI_ABORT(MPI_COMM_WORLD,1,info)
#else
         stop
#endif
     endif
     my_grackle_data%use_grackle = use_grackle
     my_grackle_data%with_radiative_cooling = grackle_with_radiative_cooling
     my_grackle_data%primordial_chemistry = grackle_primordial_chemistry
     my_grackle_data%metal_cooling = grackle_metal_cooling
     my_grackle_data%UVbackground = grackle_UVbackground
     my_grackle_data%cmb_temperature_floor = grackle_cmb_temperature_floor
     my_grackle_data%h2_on_dust = grackle_h2_on_dust
     my_grackle_data%photoelectric_heating = grackle_photoelectric_heating
     my_grackle_data%use_volumetric_heating_rate = grackle_use_volumetric_heating_rate
     my_grackle_data%use_specific_heating_rate = grackle_use_specific_heating_rate
     my_grackle_data%three_body_rate = grackle_three_body_rate
     my_grackle_data%cie_cooling = grackle_cie_cooling
     my_grackle_data%h2_optical_depth_approximation = grackle_h2_optical_depth_approximation
     my_grackle_data%ih2co = grackle_ih2co
     my_grackle_data%ipiht = grackle_ipiht
     my_grackle_data%NumberOfTemperatureBins = grackle_NumberOfTemperatureBins
     my_grackle_data%CaseBRecombination = grackle_CaseBRecombination
     my_grackle_data%Compton_xray_heating = grackle_Compton_xray_heating
     my_grackle_data%LWbackground_sawtooth_suppression = grackle_LWbackground_sawtooth_suppression
     my_grackle_data%NumberOfDustTemperatureBins = grackle_NumberOfDustTemperatureBins
     my_grackle_data%use_radiative_transfer = grackle_use_radiative_transfer
     my_grackle_data%radiative_transfer_coupled_rate_solver = grackle_radiative_transfer_coupled_rate_solver
     my_grackle_data%radiative_transfer_intermediate_step = grackle_radiative_transfer_intermediate_step
     my_grackle_data%radiative_transfer_hydrogen_only = grackle_radiative_transfer_hydrogen_only
     my_grackle_data%self_shielding_method = grackle_self_shielding_method
     my_grackle_data%Gamma = grackle_Gamma
     my_grackle_data%photoelectric_heating_rate = grackle_photoelectric_heating_rate
     my_grackle_data%HydrogenFractionByMass = grackle_HydrogenFractionByMass
     my_grackle_data%DeuteriumToHydrogenRatio = grackle_DeuteriumToHydrogenRatio
     my_grackle_data%SolarMetalFractionByMass = grackle_SolarMetalFractionByMass
     my_grackle_data%TemperatureStart = grackle_TemperatureStart
     my_grackle_data%TemperatureEnd = grackle_TemperatureEnd
     my_grackle_data%DustTemperatureStart = grackle_DustTemperatureStart
     my_grackle_data%DustTemperatureEnd = grackle_DustTemperatureEnd
     my_grackle_data%LWbackground_intensity = grackle_LWbackground_intensity
     my_grackle_data%UVbackground_redshift_on = grackle_UVbackground_redshift_on
     my_grackle_data%UVbackground_redshift_off = grackle_UVbackground_redshift_off
     my_grackle_data%UVbackground_redshift_fullon = grackle_UVbackground_redshift_fullon
     my_grackle_data%UVbackground_redshift_drop = grackle_UVbackground_redshift_drop
     my_grackle_data%cloudy_electron_fraction_factor = grackle_cloudy_electron_fraction_factor
     grackle_data_file = TRIM(grackle_data_file)//C_NULL_CHAR
     my_grackle_data%grackle_data_file = C_LOC(grackle_data_file(1:1))

     ! Grackle units
     call units(scale_l,scale_t,scale_d,scale_v,scale_nH,scale_T2)
     my_grackle_units%comoving_coordinates = 0
     my_grackle_units%density_units = scale_d
     my_grackle_units%length_units = scale_l
     my_grackle_units%time_units = scale_t
     my_grackle_units%velocity_units = scale_v
     my_grackle_units%a_units = 1.0d0
     !Set initial expansion factor (for internal units).
     !Set expansion factor to 1 for non-cosmological simulation
     ! Safety for GRACKLE initialisation
     my_grackle_units%a_value = aexp_ini

     if(cosmo) then
        my_grackle_units%comoving_coordinates = 1
        ! Reonization redshift has to be later than starting redshift
        z_reion=min(1d0/(1.1d0*aexp_ini)-1d0,z_reion)
        ! Approximate initial temperature
        T2_start=1.356d-2/aexp_ini**2
        if(nrestart==0)then
           if(myid==1)write(*,*)'Starting with T/mu (K) = ',T2_start
        end if
     endif

     ! Initialize the Grackle data
     iresult = initialize_chemistry_data(my_grackle_units)
     ! Enforce UVbackground starting redshift in grackle
     my_grackle_data%UVbackground_redshift_on = grackle_UVbackground_redshift_on
     my_grackle_data%UVbackground_redshift_fullon = grackle_UVbackground_redshift_fullon
     if(iresult.eq.0)then
         write(*,*) 'Grackle - error in initialize_chemistry_data'
#ifndef WITHOUTMPI
         call MPI_ABORT(MPI_COMM_WORLD,1,info)
#else
         stop
#endif
     endif
     my_grackle_fields%grid_dimension = C_LOC(gr_dimension)
     my_grackle_fields%grid_start = C_LOC(gr_start)
     my_grackle_fields%grid_end = C_LOC(gr_end)
     ! Point to grackle fields
     my_grackle_fields%density = C_LOC(gr_density)
     my_grackle_fields%HI_density = C_LOC(gr_HI_density)
     my_grackle_fields%HII_density = C_LOC(gr_HII_density)
     my_grackle_fields%HM_density = C_LOC(gr_HM_density)
     my_grackle_fields%HeI_density = C_LOC(gr_HeI_density)
     my_grackle_fields%HeII_density = C_LOC(gr_HeII_density)
     my_grackle_fields%HeIII_density = C_LOC(gr_HeIII_density)
     my_grackle_fields%H2I_density = C_LOC(gr_H2I_density)
     my_grackle_fields%H2II_density = C_LOC(gr_H2II_density)
     my_grackle_fields%DI_density = C_LOC(gr_DI_density)
     my_grackle_fields%DII_density = C_LOC(gr_DII_density)
     my_grackle_fields%HDI_density = C_LOC(gr_HDI_density)
     my_grackle_fields%e_density = C_LOC(gr_e_density)
     my_grackle_fields%metal_density = C_LOC(gr_metal_density)
     my_grackle_fields%internal_energy = C_LOC(gr_energy)
     my_grackle_fields%x_velocity = C_LOC(gr_x_velocity)
     my_grackle_fields%y_velocity = C_LOC(gr_y_velocity)
     my_grackle_fields%z_velocity = C_LOC(gr_z_velocity)
     my_grackle_fields%volumetric_heating_rate = C_LOC(gr_volumetric_heating_rate)
     my_grackle_fields%specific_heating_rate = C_LOC(gr_specific_heating_rate)
     my_grackle_fields%RT_HI_ionization_rate = C_LOC(gr_RT_HI_ionization_rate)
     my_grackle_fields%RT_HeI_ionization_rate = C_LOC(gr_RT_HeI_ionization_rate)
     my_grackle_fields%RT_HeII_ionization_rate = C_LOC(gr_RT_HeII_ionization_rate)
     my_grackle_fields%RT_H2_dissociation_rate = C_LOC(gr_RT_H2_dissociation_rate)
     my_grackle_fields%RT_heating_rate = C_LOC(gr_RT_heating_rate)
     do i=1,nvector
        gr_x_velocity(i) = 0d0
        gr_y_velocity(i) = 0d0
        gr_z_velocity(i) = 0d0
        gr_HII_density(i) = 0d0
        gr_HM_density(i) = 0d0
        gr_HeII_density(i) = 0d0
        gr_HeIII_density(i) = 0d0
        gr_H2I_density(i) = 0d0
        gr_H2II_density(i) = 0d0
        gr_DII_density(i) = 0d0
        gr_HDI_density(i) = 0d0
        gr_e_density(i) = 0d0
        gr_volumetric_heating_rate(i) = 0d0
        gr_specific_heating_rate(i) = 0d0
        gr_RT_HI_ionization_rate(i) = 0d0
        gr_RT_HeI_ionization_rate(i) = 0d0
        gr_RT_HeII_ionization_rate(i) = 0d0
        gr_RT_H2_dissociation_rate(i) = 0d0
        gr_RT_heating_rate(i) = 0d0
     enddo
  else
     if(cooling.and..not.(neq_chem.or.rt))then
        if(myid==1)write(*,*)'Computing cooling model'
        Nmodel=-1
        if(.not. haardt_madau)then
           Nmodel=2
        endif
        if(cosmo)then
           ! Reonization redshift has to be later than starting redshift
           z_reion=min(1d0/(1.1d0*aexp_ini)-1d0,z_reion)
           call set_model(Nmodel,dble(J21*1d-21),-1.0d0,dble(a_spec),-1.0d0,dble(z_reion), &
                & -1,2, &
                & dble(h0/100),dble(omega_b),dble(omega_m),dble(omega_l), &
                & dble(aexp_ini),T2_sim)
           T2_start=T2_sim
           if(nrestart==0)then
              if(myid==1)write(*,*)'Starting with T/mu (K) = ',T2_start
           end if
        else
           call set_model(Nmodel,dble(J21*1d-21),-1.0d0,dble(a_spec),-1.0d0,dble(z_reion), &
                & -1,2, &
                & dble(70d0/100d0),0.04d0,0.3d0,0.7d0, &
                & dble(aexp_ini),T2_sim)
        endif
     end if
  endif
#else
  if(cooling.and..not.(neq_chem.or.rt))then
     if(myid==1)write(*,*)'Computing cooling model'
     Nmodel=-1
     if(.not. haardt_madau)then
        Nmodel=2
     endif
     if(cosmo)then
        ! Reonization redshift has to be later than starting redshift
        z_reion=min(1d0/(1.1d0*aexp_ini)-1d0,z_reion)
        call set_model(Nmodel,dble(J21*1d-21),-1.0d0,dble(a_spec),-1.0d0,dble(z_reion), &
             & -1,2, &
             & dble(h0/100.),dble(omega_b),dble(omega_m),dble(omega_l), &
             & dble(aexp_ini),T2_sim)
        T2_start=T2_sim
        if(nrestart==0)then
           if(myid==1)write(*,*)'Starting with T/mu (K) = ',T2_start
        end if
     else
        call set_model(Nmodel,dble(J21*1d-21),-1.0d0,dble(a_spec),-1.0d0,dble(z_reion), &
             & -1,2, &
             & dble(70./100.),dble(0.04),dble(0.3),dble(0.7), &
             & dble(aexp_ini),T2_sim)
     endif
  end if
#endif
#ifdef RT
  if(neq_chem.or.rt) then
     if(myid==1)write(*,*)'Computing thermochemistry model'
     Nmodel=-1
     if(.not. haardt_madau)then
        Nmodel=2
     endif
     if(cosmo)then
        ! Reonization redshift has to be later than starting redshift
        z_reion=min(1d0/(1.1d0*aexp_ini)-1d0,z_reion)
        call rt_set_model(dble(h0/100.),dble(omega_b),dble(omega_m),dble(omega_l), &
             & dble(aexp_ini),T2_sim)
        T2_start=T2_sim
        if(nrestart==0)then
           if(myid==1)write(*,*)'Starting with T/mu (K) = ',T2_start
        end if
     else
        call rt_set_model(dble(70./100.),dble(0.04),dble(0.3),dble(0.7), &
             & dble(aexp_ini),T2_sim)
     endif
  end if
#endif

end subroutine init_time

subroutine init_file
  use amr_commons
  use hydro_commons
  use pm_commons
  use mpi_mod
  implicit none
  !------------------------------------------------------
  ! Read geometrical parameters in the initial condition files.
  ! Initial conditions are supposed to be made by
  ! Bertschinger's grafic version 2.0 code.
  !------------------------------------------------------
  integer:: ilevel,nx_loc,ny_loc,nz_loc
  real(sp)::dxini0,xoff10,xoff20,xoff30,astart0,omega_m0,omega_l0,h00
  character(LEN=80)::filename
  logical::ok
  integer,parameter::tag=1116
#ifndef WITHOUTMPI
  integer::dummy_io,info2
#endif

  if(verbose)write(*,*)'Entering init_file'

  ! Reading initial conditions parameters only

  nlevelmax_part=levelmin-1
  do ilevel=levelmin,nlevelmax
     if(initfile(ilevel).ne.' ')then
        filename=TRIM(initfile(ilevel))//'/ic_d'

        ! Wait for the token
#ifndef WITHOUTMPI
        if(IOGROUPSIZE>0) then
           if (mod(myid-1,IOGROUPSIZE)/=0) then
              call MPI_RECV(dummy_io,1,MPI_INTEGER,myid-1-1,tag,&
                   & MPI_COMM_WORLD,MPI_STATUS_IGNORE,info2)
           end if
        endif
#endif
        INQUIRE(file=filename,exist=ok)
        if(.not.ok)then
           if(myid==1)then
              write(*,*)'File '//TRIM(filename)//' does not exist'
           end if
           call clean_stop
        end if
        open(10,file=filename,form='unformatted')
        if(myid==1)write(*,*)'Reading file '//TRIM(filename)
        rewind 10
        read(10)n1(ilevel),n2(ilevel),n3(ilevel),dxini0 &
             & ,xoff10,xoff20,xoff30 &
             & ,astart0,omega_m0,omega_l0,h00
        close(10)

        ! Send the token
#ifndef WITHOUTMPI
        if(IOGROUPSIZE>0) then
           if(mod(myid,IOGROUPSIZE)/=0 .and.(myid.lt.ncpu))then
              dummy_io=1
              call MPI_SEND(dummy_io,1,MPI_INTEGER,myid-1+1,tag, &
                   & MPI_COMM_WORLD,info2)
           end if
        endif
#endif

        dxini(ilevel)=dxini0
        xoff1(ilevel)=xoff10
        xoff2(ilevel)=xoff20
        xoff3(ilevel)=xoff30
        nlevelmax_part=nlevelmax_part+1
     endif
  end do

  ! Check compatibility with run parameters
  nx_loc=icoarse_max-icoarse_min+1
  ny_loc=jcoarse_max-jcoarse_min+1
  nz_loc=kcoarse_max-kcoarse_min+1
  if(         nx_loc.ne.n1(levelmin)/2**levelmin &
       & .or. ny_loc.ne.n2(levelmin)/2**levelmin &
       & .or. nz_loc.ne.n3(levelmin)/2**levelmin) then
     write(*,*)'coarser grid is not compatible with initial conditions file'
     write(*,*)'Found    n1=',n1(levelmin),&
          &            ' n2=',n2(levelmin),&
          &            ' n3=',n3(levelmin)
     write(*,*)'Expected n1=',nx_loc*2**levelmin &
          &           ,' n2=',ny_loc*2**levelmin &
          &           ,' n3=',nz_loc*2**levelmin
     call clean_stop
  end if

  ! Write initial conditions parameters
  if(myid==1)then
     do ilevel=levelmin,nlevelmax_part
        write(*,'(" Initial conditions for level =",I4)')ilevel
        write(*,'(" n1=",I4," n2=",I4," n3=",I4)') &
             & n1(ilevel),&
             & n2(ilevel),&
             & n3(ilevel)
        write(*,'(" dx=",1pe10.3)')dxini(ilevel)
        write(*,'(" xoff=",1pe10.3," yoff=",1pe10.3," zoff=",&
             & 1pe10.3)') &
             & xoff1(ilevel),&
             & xoff2(ilevel),&
             & xoff3(ilevel)
     end do
  end if

end subroutine init_file

subroutine init_cosmo
  use amr_commons
  use hydro_commons
  use pm_commons
  use gadgetreadfilemod
  use mpi_mod
  implicit none
  !------------------------------------------------------
  ! Read cosmological and geometrical parameters
  ! in the initial condition files.
  ! Initial conditions are supposed to be made by
  ! Bertschinger's grafic version 2.0 code.
  !------------------------------------------------------
  integer:: ilevel
  real(sp)::dxini0,xoff10,xoff20,xoff30,astart0,omega_m0,omega_l0,h00
  character(LEN=80)::filename
  character(LEN=5)::nchar
  logical::ok
  TYPE(gadgetheadertype) :: gadgetheader
  integer::i
  integer,parameter::tag=1117
#ifndef WITHOUTMPI
  integer::dummy_io,info2
#endif

  if(verbose)write(*,*)'Entering init_cosmo'

  if(initfile(levelmin)==' ')then
     write(*,*)'You need to specifiy at least one level of initial condition'
     call clean_stop
  end if

  SELECT CASE (filetype)
  case ('grafic', 'ascii')
     ! Reading initial conditions parameters only
     aexp=2
     nlevelmax_part=levelmin-1
     do ilevel=levelmin,nlevelmax
        if(initfile(ilevel).ne.' ')then
           if(multiple)then
              call title(myid,nchar)
              filename=TRIM(initfile(ilevel))//'/dir_deltab/ic_deltab.'//TRIM(nchar)
           else
              filename=TRIM(initfile(ilevel))//'/ic_deltab'
           endif

           ! Wait for the token
#ifndef WITHOUTMPI
           if(IOGROUPSIZE>0) then
              if (mod(myid-1,IOGROUPSIZE)/=0) then
                 call MPI_RECV(dummy_io,1,MPI_INTEGER,myid-1-1,tag,&
                      & MPI_COMM_WORLD,MPI_STATUS_IGNORE,info2)
              end if
           endif
#endif

           INQUIRE(file=filename,exist=ok)
           if(.not.ok)then
              if(myid==1)then
                 write(*,*)'File '//TRIM(filename)//' does not exist'
              end if
              call clean_stop
           end if
           open(10,file=filename,form='unformatted')
           if(myid==1)write(*,*)'Reading file '//TRIM(filename)
           rewind 10
           read(10)n1(ilevel),n2(ilevel),n3(ilevel),dxini0 &
                & ,xoff10,xoff20,xoff30 &
                & ,astart0,omega_m0,omega_l0,h00
           close(10)

           ! Send the token
#ifndef WITHOUTMPI
           if(IOGROUPSIZE>0) then
              if(mod(myid,IOGROUPSIZE)/=0 .and.(myid.lt.ncpu))then
                 dummy_io=1
                 call MPI_SEND(dummy_io,1,MPI_INTEGER,myid-1+1,tag, &
                      & MPI_COMM_WORLD,info2)
              end if
           endif
#endif

           dxini(ilevel)=dxini0
           xoff1(ilevel)=xoff10
           xoff2(ilevel)=xoff20
           xoff3(ilevel)=xoff30
           astart(ilevel)=astart0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !BiP 15/02/2021 specific for the cas of Dirac-Milne : impose that the first species (particles is in same proportion than second species (hydro)
           omega_m=1.0!
	   omega_b=0.5
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
           omega_l=omega_l0
           h0=h00
           aexp=MIN(aexp,astart(ilevel))
           nlevelmax_part=nlevelmax_part+1
           ! Compute SPH equivalent mass (initial gas mass resolution)
           mass_sph=omega_b/omega_m*0.5d0**(ndim*ilevel)

        endif
     end do

     ! Compute initial expansion factor
     if(aexp_ini.lt.1.0)then
        aexp=aexp_ini
     else
        aexp_ini=aexp
     endif

     ! Check compatibility with run parameters
     if(.not. multiple) then
        if(         nx.ne.n1(levelmin)/2**levelmin &
             & .or. ny.ne.n2(levelmin)/2**levelmin &
             & .or. nz.ne.n3(levelmin)/2**levelmin) then
           write(*,*)'coarser grid is not compatible with initial conditions file'
           write(*,*)'Found    n1=',n1(levelmin),&
                &            ' n2=',n2(levelmin),&
                &            ' n3=',n3(levelmin)
           write(*,*)'Expected n1=',nx*2**levelmin &
                &           ,' n2=',ny*2**levelmin &
                &           ,' n3=',nz*2**levelmin
           call clean_stop
        endif
     end if

     ! Compute box length in the initial conditions in units of h-1 Mpc
     boxlen_ini=dble(nx)*2**levelmin*dxini(levelmin)*(h0/100)

  CASE ('gadget')
     if (verbose) write(*,*)'Reading in gadget format from '//TRIM(initfile(levelmin))
     call gadgetreadheader(TRIM(initfile(levelmin)), 0, gadgetheader, ok)
     if(.not.ok) call clean_stop
     do i=1,6
        if (i .ne. 2) then
           if (gadgetheader%nparttotal(i) .ne. 0) then
              write(*,*) 'Non DM particles present in bin ', i
              call clean_stop
           endif
        endif
     enddo
     if (gadgetheader%mass(2) == 0) then
        write(*,*) 'Particles have different masses, not supported'
        call clean_stop
     endif
     omega_m = gadgetheader%omega0
     omega_l = gadgetheader%omegalambda
     h0 = gadgetheader%hubbleparam * 100d0
     boxlen_ini = gadgetheader%boxsize
     aexp = gadgetheader%time
     aexp_ini = aexp
     ! Compute SPH equivalent mass (initial gas mass resolution)
     mass_sph=omega_b/omega_m*0.5d0**(ndim*levelmin)
     nlevelmax_part = levelmin
     astart(levelmin) = aexp
     xoff1(levelmin)=0
     xoff2(levelmin)=0
     xoff3(levelmin)=0
     dxini(levelmin) = boxlen_ini/(nx*2**levelmin*(h0/100))

  CASE DEFAULT
     write(*,*) 'Unsupported input format '//filetype
     call clean_stop
  END SELECT

  ! Write cosmological parameters
  if(myid==1)then
     write(*,'(" Cosmological parameters:")')
     write(*,'(" aexp=",1pe10.3," H0=",1pe10.3," km s-1 Mpc-1")')aexp,h0
     write(*,'(" omega_m=",F7.3," omega_l=",F7.3," omega_b=",F7.3)')omega_m,omega_l,omega_b
     write(*,'(" box size=",1pe10.3," h-1 Mpc")')boxlen_ini
  end if
  omega_k=1d0-omega_l-omega_m

  ! Compute linear scaling factor between aexp and astart(ilevel)
  do ilevel=levelmin,nlevelmax_part
     dfact(ilevel)=d1a(aexp)/d1a(astart(ilevel))
     vfact(ilevel)=astart(ilevel)*fpeebl(astart(ilevel)) & ! Same scale factor as in grafic1
          & *sqrt(omega_m/astart(ilevel)+omega_l*astart(ilevel)*astart(ilevel)+omega_k) &
          & /astart(ilevel)*h0
  end do

  ! Write initial conditions parameters
  do ilevel=levelmin,nlevelmax_part
     if(myid==1)then
        write(*,'(" Initial conditions for level =",I4)')ilevel
        write(*,'(" dx=",1pe10.3," h-1 Mpc")')dxini(ilevel)*h0/100.
     endif
     if(.not.multiple)then
        if(myid==1)then
           write(*,'(" n1=",I4," n2=",I4," n3=",I4)') &
                & n1(ilevel),&
                & n2(ilevel),&
                & n3(ilevel)
           write(*,'(" xoff=",1pe10.3," yoff=",1pe10.3," zoff=",&
                & 1pe10.3," h-1 Mpc")') &
                & xoff1(ilevel)*h0/100.,&
                & xoff2(ilevel)*h0/100.,&
                & xoff3(ilevel)*h0/100.
        endif
     else
        write(*,'(" myid=",I4," n1=",I4," n2=",I4," n3=",I4)') &
             & myid,n1(ilevel),n2(ilevel),n3(ilevel)
        write(*,'(" myid=",I4," xoff=",1pe10.3," yoff=",1pe10.3," zoff=",&
             & 1pe10.3," h-1 Mpc")') &
             & myid,xoff1(ilevel)*h0/100.,&
             & xoff2(ilevel)*h0/100.,&
             & xoff3(ilevel)*h0/100.
     endif
  end do

  ! Scale displacement in Mpc to code velocity (v=dx/dtau)
  ! in coarse cell units per conformal time
  vfact(1)=aexp*fpeebl(aexp)*sqrt(omega_m/aexp+omega_l*aexp*aexp+omega_k)
  ! This scale factor is different from vfact in grafic by h0/aexp

contains

  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  function fy(a)
    implicit none
    !      Computes the integrand
    real(dp)::fy
    real(dp)::y,a

    y=omega_m*(1d0/a-1d0) + omega_l*(a*a-1d0) + 1d0
    fy=1d0/y**1.5d0

    return
  end function fy
  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  function d1a(a)
    implicit none
    real(dp)::d1a
    !     Computes the linear growing mode D1 in a Friedmann-Robertson-Walker
    !     universe. See Peebles LSSU sections 11 and 14.
    real(dp)::a,y12,y,eps

    eps=1.0d-6
    if(a .le. 0.0d0)then
       write(*,*)'a=',a
       call clean_stop
    end if
    y=omega_m*(1d0/a-1d0) + omega_l*(a*a-1d0) + 1d0
    if(y .lt. 0.0D0)then
       write(*,*)'y=',y
       call clean_stop
    end if
    y12=y**0.5d0
    d1a=y12/a*rombint(eps,a,eps)

    return
  end function d1a
  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!!$  function ad1(d1)
!!$    implicit none
!!$    real(dp)::ad1
!!$    real(dp)::a,d1,da
!!$    integer::niter
!!$    ! Inverts the relation d1(a) given by function d1a(a) using
!!$    ! Newton-Raphson.
!!$    if (d1.eq.0.0) stop 'ad1 undefined for d1=0!'
!!$    ! Initial guess for Newton-Raphson iteration, good for Omega near 1.
!!$    a=1.e-7
!!$    niter=0
!!$10  niter=niter+1
!!$    da=(d1/d1a(a)-1d0)/fpeebl(a)*a
!!$    a=a+da
!!$    if (abs(da).gt.1.0e-8.and.niter.lt.10) go to 10
!!$    ad1=a
!!$    return
!!$  end function ad1
  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  function fpeebl(a)
    implicit none
    real(dp) :: fpeebl,a
    !     Computes the growth factor f=d\log D1/d\log a.
    real(dp) :: fact,y,eps

    eps=1.0d-6
    y=omega_m*(1d0/a-1d0) + omega_l*(a*a-1d0) + 1d0
    fact=rombint(eps,a,eps)
    fpeebl=(omega_l*a*a-0.5d0*omega_m/a)/y - 1d0 + a*fy(a)/fact
    return
  end function fpeebl
  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  function rombint(a,b,tol)
    implicit none
    real(dp)::rombint
    !
    !     Rombint returns the integral from a to b of f(x)dx using Romberg
    !     integration. The method converges provided that f(x) is continuous
    !     in (a,b). The function f must be double precision and must be
    !     declared external in the calling routine.
    !     tol indicates the desired relative accuracy in the integral.
    !
    integer::maxiter=16,maxj=5
    real(dp),dimension(100):: g
    real(dp)::a,b,tol,fourj
    real(dp)::h,error,gmax,g0,g1
    integer::nint,i,j,k,jmax

    h=0.5d0*(b-a)
    gmax=h*(fy(a)+fy(b))
    g(1)=gmax
    nint=1
    error=1.0d20
    i=0
10  i=i+1
    if(.not.  (i>maxiter.or.(i>5.and.abs(error)<tol)))then
       !  Calculate next trapezoidal rule approximation to integral.

       g0=0.0d0
       do k=1,nint
          g0=g0+fy(a+(k+k-1)*h)
       end do
       g0=0.5d0*g(1)+h*g0
       h=0.5d0*h
       nint=nint+nint
       jmax=min(i,maxj)
       fourj=1.0d0

       do j=1,jmax
          ! Use Richardson extrapolation.
          fourj=4.0d0*fourj
          g1=g0+(g0-g(j))/(fourj-1.0d0)
          g(j)=g0
          g0=g1
       enddo
       if (abs(g0).gt.tol) then
          error=1.0d0-gmax/g0
       else
          error=gmax
       end if
       gmax=g0
       g(jmax+1)=g0
       go to 10
    end if
    rombint=g0
    if (i>maxiter.and.abs(error)>tol) &
         &    write(*,*) 'Rombint failed to converge; integral, error=', &
         &    rombint,error
    return
  end function rombint

end subroutine init_cosmo

subroutine friedman(O_mat_0,O_vac_0,O_k_0,alpha,axp_min, &
     & axp_out,hexp_out,tau_out,t_out,ntable)
  use amr_parameters
  implicit none
  integer::ntable
  real(kind=8)::O_mat_0, O_vac_0, O_k_0
  real(kind=8)::alpha,axp_min
  real(dp),dimension(0:ntable)::axp_out,hexp_out,tau_out,t_out
  ! ######################################################!
  ! This subroutine assumes that axp = 1 at z = 0 (today) !
  ! and that t and tau = 0 at z = 0 (today).              !
  ! axp is the expansion factor, hexp the Hubble constant !
  ! defined as hexp=1/axp*daxp/dtau, tau the conformal    !
  ! time, and t the look-back time, both in unit of 1/H0. !
  ! alpha is the required accuracy and axp_min is the     !
  ! starting expansion factor of the look-up table.       !
  ! ntable is the required size of the look-up table.     !
  ! ######################################################!
  real(kind=8)::axp_tau, axp_t
  real(kind=8)::axp_tau_pre, axp_t_pre
  real(kind=8)::dadtau, dadt
  real(kind=8)::dtau,dt
  real(kind=8)::tau,t
  integer::nstep,nout,nskip

  if( (O_mat_0+O_vac_0+O_k_0) .ne. 1.0D0 )then
     write(*,*)'Error: non-physical cosmological constants'
     write(*,*)'O_mat_0,O_vac_0,O_k_0=',O_mat_0,O_vac_0,O_k_0
     write(*,*)'The sum must be equal to 1.0, but '
     write(*,*)'O_mat_0+O_vac_0+O_k_0=',O_mat_0+O_vac_0+O_k_0
     call clean_stop
  end if

  axp_tau = 1.0D0
  axp_t = 1.0D0
  tau = 0.0D0
  t = 0.0D0
  nstep = 0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!BiP: Specific to Dirac-Milne : impose Milne background

O_mat_0=0.0
O_vac_0=0.0d0
O_k_0=1.0d0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!BiP

  do while ( (axp_tau .ge. axp_min) .or. (axp_t .ge. axp_min) )

     nstep = nstep + 1
     dtau = alpha * axp_tau / dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)
     axp_tau_pre = axp_tau - dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)*dtau/2d0
     axp_tau = axp_tau - dadtau(axp_tau_pre,O_mat_0,O_vac_0,O_k_0)*dtau
     tau = tau - dtau

     dt = alpha * axp_t / dadt(axp_t,O_mat_0,O_vac_0,O_k_0)
     axp_t_pre = axp_t - dadt(axp_t,O_mat_0,O_vac_0,O_k_0)*dt/2d0
     axp_t = axp_t - dadt(axp_t_pre,O_mat_0,O_vac_0,O_k_0)*dt
     t = t - dt

  end do

  if(debug)then
     write(*,666)-t
  end if
  666 format(' Age of the Universe (in unit of 1/H0)=',1pe10.3)

  nskip=nstep/ntable

  axp_t = 1d0
  t = 0d0
  axp_tau = 1d0
  tau = 0d0
  nstep = 0
  nout=0
  t_out(nout)=t
  tau_out(nout)=tau
  axp_out(nout)=axp_tau
  hexp_out(nout)=dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)/axp_tau

  do while ( (axp_tau .ge. axp_min) .or. (axp_t .ge. axp_min) )

     nstep = nstep + 1
     dtau = alpha * axp_tau / dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)
     axp_tau_pre = axp_tau - dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)*dtau/2d0
     axp_tau = axp_tau - dadtau(axp_tau_pre,O_mat_0,O_vac_0,O_k_0)*dtau
     tau = tau - dtau

     dt = alpha * axp_t / dadt(axp_t,O_mat_0,O_vac_0,O_k_0)
     axp_t_pre = axp_t - dadt(axp_t,O_mat_0,O_vac_0,O_k_0)*dt/2d0
     axp_t = axp_t - dadt(axp_t_pre,O_mat_0,O_vac_0,O_k_0)*dt
     t = t - dt

     if(mod(nstep,nskip)==0)then
        nout=nout+1
        t_out(nout)=t
        tau_out(nout)=tau
        axp_out(nout)=axp_tau
        hexp_out(nout)=dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)/axp_tau
     end if

  end do
  t_out(ntable)=t
  tau_out(ntable)=tau
  axp_out(ntable)=axp_tau
  hexp_out(ntable)=dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)/axp_tau

end subroutine friedman

function dadtau(axp_tau,O_mat_0,O_vac_0,O_k_0)
  use amr_parameters
  real(kind=8)::dadtau,axp_tau,O_mat_0,O_vac_0,O_k_0
  dadtau = axp_tau*axp_tau*axp_tau *  &
       &   ( O_mat_0 + &
       &     O_vac_0 * axp_tau*axp_tau*axp_tau + &
       &     O_k_0   * axp_tau )
  dadtau = sqrt(dadtau)
  return
end function dadtau

function dadt(axp_t,O_mat_0,O_vac_0,O_k_0)
  use amr_parameters
  real(kind=8)::dadt,axp_t,O_mat_0,O_vac_0,O_k_0
  dadt   = (1.0D0/axp_t)* &
       &   ( O_mat_0 + &
       &     O_vac_0 * axp_t*axp_t*axp_t + &
       &     O_k_0   * axp_t )
  dadt = sqrt(dadt)
  return
end function dadt