text:
      SUBROUTINE DIFF_EQNS
     &           (STAGE,   NSD,
     &            NNW_SD,  NGP_SD,
     &            NW_IPAR, NW_RPAR,
     &            CRD,     PHI,     PHN,
     &            DCRD_DT, DPHI_DT, DPHN_DT,
     &            S0,      S1,      D0,      D1)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'bct_params.inc'
      INCLUDE 'chk_params.inc'
      INCLUDE 'tme_params.inc'
      INCLUDE 'tme_funcs.inc'
C
      INTEGER(KIND=IK) STAGE
      INTEGER(KIND=IK) NSD
      INTEGER(KIND=IK) NNW_SD (*)
      INTEGER(KIND=IK) NGP_SD (*)
      INTEGER(KIND=IK) NW_IPAR (N_IP, *)
      REAL   (KIND=RK) NW_RPAR (N_RP, *)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) DCRD_DT (4, *)
      REAL   (KIND=RK) DPHI_DT (2, *)
      REAL   (KIND=RK) DPHN_DT (2, *)
      REAL   (KIND=RK) S0 (*)
      REAL   (KIND=RK) S1 (*)
      REAL   (KIND=RK) D0 (*)
      REAL   (KIND=RK) D1 (*)
C     Solve Delta(Phi) = 0 (by means of Boundary Integral Equations)
      CALL SOLVE_LAPLACE
     &     (BCD_PHI, NSD,
     &      NNW_SD,  NGP_SD,
     &      NW_IPAR, NW_RPAR,
     &      CRD,     PHI,     PHN,
     &      S0,      S1,      D0,      D1)
C
      IF (PROB_TIME_DEPENDENT .AND. GET_TIME () .LE. GET_TEND ()) THEN
C        Compute the righthand sides the the time-dependent DE's
         CALL RIGHTHAND_SIDES
     &        (STAGE,   NSD,
     &         NNW_SD,  NGP_SD,
     &         NW_IPAR, NW_RPAR,
     &         CRD,     PHI,     PHN,
     &         DCRD_DT, DPHI_DT, DPHN_DT)
C        Solve Laplace equation for Phi_t
C        CALL SOLVE_LAPLACE
C    &     (BCD_PHI_T, NSD,
C    &      NNW_SD,    NGP_SD,
C    &      NW_IPAR,   NW_RPAR,
C    &      CRD,       DPHI_DT,   DPHN_DT,
C    &      S0,        S1,        D0,        D1)
      END IF
C
      IF (STAGE .EQ. 1 .AND. I_PLOT .EQ. 0) THEN
         CALL PLTOUT
     &        (NSD,     NNW_SD,  NW_IPAR,
     &         CRD,      PHI,     PHN,
     &         DCRD_DT, DPHI_DT, DPHN_DT)
         IF (CHK_CONTOUR) THEN
            CALL CONTOUR_INTEGRAL
     &           (NSD, NNW_SD, NGP_SD, NW_IPAR, CRD, PHN)
         END IF
      END IF
C
      RETURN
      END
      SUBROUTINE RIGHTHAND_SIDES
     &           (STAGE,   NSD,
     &            NNW_SD,  NGP_SD,
     &            NW_IPAR, NW_RPAR,
     &            CRD,     PHI,     PHN,
     &            DCRD_DT, DPHI_DT, DPHN_DT)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
C
      INTEGER(KIND=IK) STAGE
      INTEGER(KIND=IK) NSD
      INTEGER(KIND=IK) NNW_SD (*)
      INTEGER(KIND=IK) NGP_SD (*)
      INTEGER(KIND=IK) NW_IPAR (N_IP, *)
      REAL   (KIND=RK) NW_RPAR (N_RP, *)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) DCRD_DT (4, *)
      REAL   (KIND=RK) DPHI_DT (2, *)
      REAL   (KIND=RK) DPHN_DT (2, *)
C
      INTEGER(KIND=IK) ISD, INW, IGP, NNW, NGP
C
      INW = 1
      IGP = 1
C
      DO ISD = 1, NSD
         NNW = NNW_SD (ISD)
         NGP = NGP_SD (ISD)
         CALL RHS_SD
     &        (STAGE,             NNW,              NGP,
     &         NW_IPAR (1, INW),  NW_RPAR (1, INW),
     &         CRD     (1, IGP),  PHI     (1, IGP), PHN     (1, IGP),
     &         DCRD_DT (1, IGP),  DPHI_DT (1, IGP), DPHN_DT (1, IGP))
         INW = INW + NNW
         IGP = IGP + NGP
      END DO
C
      RETURN
      END
      SUBROUTINE RHS_SD
     &           (STAGE,   NNW,     NGP,
     &            NW_IPAR, NW_RPAR,
     &            CRD,     PHI,     PHN,
     &            DCRD_DT, DPHI_DT, DPHN_DT)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'tme_params.inc'
C
      INTEGER(KIND=IK) STAGE
      INTEGER(KIND=IK) NNW
      INTEGER(KIND=IK) NGP
      INTEGER(KIND=IK) NW_IPAR (N_IP, *)
      REAL   (KIND=RK) NW_RPAR (N_RP, *)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) DCRD_DT (4, *)
      REAL   (KIND=RK) DPHI_DT (2, *)
      REAL   (KIND=RK) DPHN_DT (2, *)
C
      CALL PARTIAL_DERIVS
     &     (STAGE,   NNW,
     &      NW_IPAR, NW_RPAR,
     &      CRD,     PHI,     PHN,
     &      DCRD_DT, DPHI_DT, DPHN_DT)
C
      IF (GRID_TIME_DEPENDENT) THEN
         CALL MATERIAL_DERIVS
     &        (STAGE,   NNW,     NGP,
     &         NW_IPAR, NW_RPAR,
     &         CRD,     PHI,     PHN,
     &         DCRD_DT, DPHI_DT, DPHN_DT)
      END IF
C
      RETURN
      END
      SUBROUTINE PARTIAL_DERIVS
     &           (STAGE,   NNW,
     &            NW_IPAR, NW_RPAR,
     &            CRD,     PHI,     PHN,
     &            CRD_T,   PHI_T,   PHN_T)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'bct_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'net_funcs.inc'
C
      INTEGER(KIND=IK) STAGE             ! Current RK stage (IN)
      INTEGER(KIND=IK) NNW               ! nr. of networks (IN)
      INTEGER(KIND=IK) NW_IPAR (N_IP, *) ! integer netw. params. (IN)
      REAL   (KIND=RK) NW_RPAR (N_RP, *) ! real netw. params. (IN)
      REAL   (KIND=RK) CRD (4, *)        ! grid point coords. (IN)
      REAL   (KIND=RK) PHI (2, *)        ! potential (IN)
      REAL   (KIND=RK) PHN (2, *)        ! normal deriv of Phi (IN)
      REAL   (KIND=RK) CRD_T (4, *)      ! partial d/dT of CRD (OUT)
      REAL   (KIND=RK) PHI_T (2, *)      ! partial d/dT of PHI (OUT)
      REAL   (KIND=RK) PHN_T (2, *)      ! partial d/dT of PHN (OUT)
C
      INTEGER(KIND=IK) INW, IGP, JGP, NGP, BCT
      REAL   (KIND=RK) X_XI, Z_XI, J_I, U_L (2), U_G (2)
C
      IGP = 1
C
      DO INW = 1, NNW
         NGP = GET_NGP (NW_IPAR (1, INW))
         BCT = GET_BCT (NW_IPAR (1, INW))
         DO JGP = IGP, IGP + NGP - 1
            CRD_T (1, JGP) = 0.0D0
            CRD_T (2, JGP) = 0.0D0
            PHI_T (1, JGP) = 0.0D0
            PHN_T (1, JGP) = 0.0D0
         END DO
         IF (BCT .EQ. -BCT_BERNOULLI) THEN
            CALL BERNOULLI
     &           (NGP,
     &            NW_IPAR (1, INW), NW_RPAR (1, INW),
     &            CRD     (1, IGP), PHI     (1, IGP), PHN (1, IGP),
     &            CRD_T   (1, IGP), PHI_T   (1, IGP))
         ELSE IF (BCT .EQ. -BCT_LINBERNOU) THEN
            CALL LINBERNOU
     &           (NGP,
     &            NW_IPAR (1, INW), NW_RPAR (1, INW),
     &            CRD     (1, IGP), PHI     (1, IGP), PHN (1, IGP),
     &            CRD_T   (1, IGP), PHI_T   (1, IGP))
         ELSE IF (BCT .EQ. -BCT_SOMMERFELD) THEN
            CALL SOMMERFELD
     &           (NGP,
     &            NW_IPAR (1, INW), NW_RPAR (1, INW),
     &            CRD     (1, IGP), PHI     (1, IGP), PHN (1, IGP),
     &            CRD_T   (1, IGP), PHI_T   (1, IGP))
         ELSE IF (BCT .EQ. -BCT_TIME_POLY) THEN
            CALL TIME_POLY
     &           (NGP,
     &            NW_IPAR (1, INW), NW_RPAR (1, INW),
     &            CRD     (1, IGP), PHI     (1, IGP), PHN (1, IGP),
     &            CRD_T   (1, IGP), PHI_T   (1, IGP))
         END IF
         DO JGP = IGP, IGP + NGP - 1
            X_XI = CRD (3, JGP)
            Z_XI = CRD (4, JGP)
            J_I  = 1.0D0 / SQRT (X_XI ** 2 + Z_XI ** 2)
            U_L (1) = PHI (2, JGP) * J_I
            U_L (2) = PHN (1, JGP)
            U_G (1) = (X_XI * U_L (1) - Z_XI * U_L (2)) * J_I
            U_G (2) = (X_XI * U_L (2) + Z_XI * U_L (1)) * J_I
            CALL LOCAL_TO_GLOBAL (1, 2, CRD (1, JGP), U_L, U_G)
C           WRITE (*, 1001) JGP, PHI_T (1, JGP), U_G
         END DO
         IGP = IGP + NGP
      END DO
C
      RETURN
      END
      SUBROUTINE BERNOULLI
     &           (NGP,
     &            NW_IPAR, NW_RPAR,
     &            CRD,     PHI,     PHN,
     &            CRD_T,   PHI_T)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'rle_params.inc'
C
      INTEGER(KIND=IK) NGP
      INTEGER(KIND=IK) NW_IPAR (*)
      REAL   (KIND=RK) NW_RPAR (*)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) CRD_T (4, *)
      REAL   (KIND=RK) PHI_T (2, *)
C
      INTEGER(KIND=IK) IGP
      REAL   (KIND=RK) X_XI, Z_XI, J_I, U_L (2), Q
C
      DO IGP = 1, NGP
         X_XI = CRD (3, IGP)
         Z_XI = CRD (4, IGP)
         J_I  = 1.0D0 / SQRT (X_XI ** 2 + Z_XI ** 2)
         U_L (1) = PHI (2, IGP) * J_I
         U_L (2) = PHN (1, IGP)
         Q       = U_L (1) ** 2 + U_L (2) ** 2
C        phi_t + 1/2 grad (phi) ** 2 + g z = 0
         PHI_T (1, IGP) =  -0.5D0 * Q - GRAV * CRD (2, IGP)
      END DO
C
      RETURN
      END
      SUBROUTINE LINBERNOU
     &           (NGP,
     &            NW_IPAR, NW_RPAR,
     &            CRD,     PHI,     PHN,
     &            CRD_T,   PHI_T)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'rle_params.inc'
      INCLUDE 'tme_funcs.inc'
C
      INTEGER(KIND=IK) NGP
      INTEGER(KIND=IK) NW_IPAR (*)
      REAL   (KIND=RK) NW_RPAR (*)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) CRD_T (4, *)
      REAL   (KIND=RK) PHI_T (2, *)
C
      INTEGER(KIND=IK) IGP
      REAL   (KIND=RK) X_XI, Z_XI, J_I, U_L (2), U_G (2)
      REAL   (KIND=RK) ETA (2, NGP)
      INTEGER(KIND=IK) BC (2)
      DATA             BC /0, 0/
C
      CALL DCOPY (NGP, CRD (2, 1), 4, ETA (1, 1), 2)
      CALL Spline (BC, NGP, 1, 2, ETA (1, 1), ETA (2, 1))
      DO IGP = 1, NGP
         X_XI = CRD (3, IGP)
         Z_XI = CRD (4, IGP)
         J_I  = 1.0D0 / SQRT (X_XI ** 2 + Z_XI ** 2)
         U_L (1) = PHI (2, IGP) * J_I
         U_L (2) = PHN (1, IGP)
         U_G (1) = (X_XI * U_L (1) - Z_XI * U_L (2)) * J_I
         U_G (2) = (X_XI * U_L (2) + Z_XI * U_L (1)) * J_I
C        phi_t = -g * eta
         PHI_T (1, IGP) =  -GRAV * CRD (2, IGP)
C        phi_xi_t = -g * eta
         PHI_T (2, IGP) =  -GRAV * ETA (2, IGP)
C        eta_t =  Phi_z
         CRD_T (2, IGP) = U_G (2)
C        WRITE (*, *) IGP, ':', PHI_T (1, IGP), PHI_T (2, IGP)
      END DO
C
      RETURN
      END
      SUBROUTINE SOMMERFELD
     &           (NGP,
     &            NW_IPAR, NW_RPAR,
     &            CRD,     PHI,     PHN,
     &            CRD_T,   PHI_T)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
C
      INTEGER(KIND=IK) NGP
      INTEGER(KIND=IK) NW_IPAR (*)
      REAL   (KIND=RK) NW_RPAR (*)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) CRD_T (4, *)
      REAL   (KIND=RK) PHI_T (2, *)
C
      INTEGER(KIND=IK) IGP
      REAL   (KIND=RK) PAR (N_BCP), C, B
C
      CALL GET_BCP (NW_RPAR, PAR)
C
      C = PAR (1)
      B = PAR (2)
C     B = g * (R- H) - 1/2 * C ** 2 (see page 48 of PdH's thesis)
      DO IGP = 1, NGP
         PHI_T (1, IGP) =  C * PHN (1, IGP) - B
      END DO
C
      RETURN
      END
      SUBROUTINE TIME_POLY
     &           (NGP,
     &            NW_IPAR, NW_RPAR,
     &            CRD,     PHI,     PHN,
     &            CRD_T,   PHI_T)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'anl_params.inc'
      INCLUDE 'tme_funcs.inc'
C
      INTEGER(KIND=IK) NGP
      INTEGER(KIND=IK) NW_IPAR (*)
      REAL   (KIND=RK) NW_RPAR (*)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) CRD_T (4, *)
      REAL   (KIND=RK) PHI_T (2, *)
C
      INTEGER(KIND=IK) IGP
      REAL   (KIND=RK) C, T
C
      T = GET_TIME ()
      C = ANL_PAR (5)
      DO IGP = 1, NGP
         PHI_T (1, IGP) = 3.0D0 * C * T ** 2
      END DO
C
      RETURN
      END
      SUBROUTINE MATERIAL_DERIVS
     &           (STAGE,   NNW,     NGP,
     &            NW_IPAR, NW_RPAR,
     &            CRD,     PHI,     PHN,
     &            DCRD_DT, DPHI_DT, DPHN_DT)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
C
      INTEGER(KIND=IK) STAGE
      INTEGER(KIND=IK) NNW
      INTEGER(KIND=IK) NGP
      INTEGER(KIND=IK) NW_IPAR (N_IP, *)
      REAL   (KIND=RK) NW_RPAR (N_RP, *)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) DCRD_DT (4, *)
      REAL   (KIND=RK) DPHI_DT (2, *)
      REAL   (KIND=RK) DPHN_DT (2, *)
C
      CALL GRID_VELOCITY
     &     (NNW, NGP, NW_IPAR, NW_RPAR, CRD, PHI, PHN, DCRD_DT)
C
      CALL TO_MATERIAL
     &     (NGP, CRD, PHI, PHN, DCRD_DT, DPHI_DT)
C
      RETURN
      END
      SUBROUTINE GRID_VELOCITY
     &           (NNW, NGP, NW_IPAR, NW_RPAR, CRD, PHI, PHN, V_GRID)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
C
      INTEGER(KIND=IK) NNW
      INTEGER(KIND=IK) NGP
      INTEGER(KIND=IK) NW_IPAR (N_IP, *)
      REAL   (KIND=RK) NW_RPAR (N_RP, *)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) V_GRID (4, *)
C
      CALL BASIC_VELOCITY
     &     (NGP, CRD, PHI, PHN, V_GRID)
C
      CALL CORRECT_TANGENTIAL
     &     (NNW, NW_IPAR, NW_RPAR, CRD, PHI, PHN, V_GRID)
C
      CALL ALIGN_NORMAL
     &     (NNW, NW_IPAR, NW_RPAR, CRD, PHI, PHN, V_GRID)
C
      CALL ALIGN_TANGENTIAL
     &     (NNW, NW_IPAR, NW_RPAR, CRD, PHI, PHN, V_GRID)
C
      RETURN
      END
      SUBROUTINE BASIC_VELOCITY
     &           (NGP, CRD, PHI, PHN, V_GRID)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
C
      INTEGER(KIND=IK) NGP
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) V_GRID (4, *)
C
      INTEGER(KIND=IK) JGP
      REAL   (KIND=RK) X_XI, Z_XI, J_INV
C
      DO JGP = 1, NGP
         X_XI  = CRD (3, JGP)
         Z_XI  = CRD (4, JGP)
         J_INV = 1.0D0 / SQRT (X_XI ** 2 + Z_XI ** 2)
C        d/dT (X) = Grad (Phi) = (phi_s, phi_n)
         V_GRID (1, JGP) = PHI (2, JGP) * J_INV
         V_GRID (2, JGP) = PHN (1, JGP)
      END DO
C
      RETURN
      END
      SUBROUTINE CORRECT_TANGENTIAL
     &           (NNW, NW_IPAR, NW_RPAR, CRD, PHI, PHN, V_GRID)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'bct_params.inc'
      INCLUDE 'net_funcs.inc'
      INCLUDE 'tme_funcs.inc'
C
      INTEGER(KIND=IK) NNW
      INTEGER(KIND=IK) NW_IPAR (N_IP, *)
      REAL   (KIND=RK) NW_RPAR (N_RP, *)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) V_GRID (4, *)
C
      INTEGER(KIND=IK) INW, IGP, JGP, NGP, BCT
      REAL   (KIND=RK) PAR (N_GMP)
      REAL   (KIND=RK) ALPHA, DT
      REAL   (KIND=RK) S_XI (100), S_XIXI (100)
      INTEGER(KIND=IK) BC (2)
      DATA             BC /0, 0/
C
      IGP = 1
C
      DO INW = 1, NNW
         NGP = GET_NGP (NW_IPAR (1, INW))
         BCT = GET_BCT (NW_IPAR (1, INW))
         IF (BCT .EQ. -BCT_BERNOULLI) THEN
            CALL GET_GMP (NW_RPAR (1, INW), PAR)
            DT    = GET_DELTA_T ()
            ALPHA = PAR (1) / DT
C           only correct the tangential velocities of networks that
C           are  bernoulli kind of networks
            IF (ALPHA .NE. 0.0D0) THEN
               DO JGP = IGP, IGP + NGP - 1
                  S_XI(JGP) = SQRT (CRD (3,JGP) ** 2 + CRD (4,JGP) ** 2)
               END DO
               CALL SPLINE (BC, NGP, 1, 1, S_XI, S_XIXI)
               DO JGP = IGP + 1, IGP + NGP - 2
                  V_GRID (1,JGP) = V_GRID (1,JGP) + ALPHA * S_XIXI (JGP)
               END DO
            END IF
         END IF
         IGP = IGP + NGP
      END DO
C
      RETURN
      END
      SUBROUTINE ALIGN_NORMAL
     &           (NNW, NW_IPAR, NW_RPAR, CRD, PHI, PHN, V_GRID)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'bct_params.inc'
      INCLUDE 'net_funcs.inc'
C
      INTEGER(KIND=IK) NNW
      INTEGER(KIND=IK) NW_IPAR (N_IP, *)
      REAL   (KIND=RK) NW_RPAR (N_RP, *)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) V_GRID (4, *)
C
      INTEGER(KIND=IK) INW, IGP, JGP, NGP, IED, ANW, AGP, ACT, BCT
      REAL   (KIND=RK) U_A (2)
C
      IGP = 1
C
      DO INW = 1, NNW
         NGP = GET_NGP (NW_IPAR (1, INW))
         BCT = GET_BCT (NW_IPAR (1, INW))
         DO IED = 1, 2
            ANW = GET_ANW (IED, NW_IPAR (1, INW))
            ACT = GET_BCT      (NW_IPAR (1, ANW))
            IF  (ACT .EQ. -BCT_BERNOULLI
     &     .AND. BCT .NE.  BCT_WAVEMAKER) THEN
C              only align the normal velocities of lateral networks that
C              are coupled to a network with a bernoulli kind of bc
               JGP = IGP + (IED - 1) * (NGP - 1)
               ANW = GET_ANW (IED,      NW_IPAR (1, INW))
               AGP = GET_AGP (IED, INW, NW_IPAR)
               CALL EDGE_VELOCITY
     &              (JGP, AGP, CRD, V_GRID, U_A)
C              Make the lateral network have a uniform normal velocity
               DO JGP = IGP, IGP + NGP - 1
                  V_GRID (2, JGP) = U_A (2)
C                 V_GRID (2, JGP) = PHN (1, IGP + (NGP - 1) * (IED - 1))
               END DO
            END IF
         END DO
         IGP = IGP + NGP
      END DO
C
      RETURN
      END
      SUBROUTINE ALIGN_TANGENTIAL
     &           (NNW, NW_IPAR, NW_RPAR, CRD, PHI, PHN, V_GRID)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'bct_params.inc'
      INCLUDE 'net_funcs.inc'
C
      INTEGER(KIND=IK) NNW
      INTEGER(KIND=IK) NW_IPAR (N_IP, *)
      REAL   (KIND=RK) NW_RPAR (N_RP, *)
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) V_GRID (4, *)
C
      INTEGER(KIND=IK) INW, IGP, JGP, NGP, IED, AGP, BCT
      REAL   (KIND=RK) U_A (2, 2), DVS
C
      IGP = 1
C
      DO INW = 1, NNW
         NGP = GET_NGP (NW_IPAR (1, INW))
         BCT = GET_BCT (NW_IPAR (1, INW))
         IF (BCT .NE. -BCT_BERNOULLI) THEN
C           only align the tangential velocities of networks that
C           are NOT bernoulli kind of networks
            DO IED = 1, 2
               JGP = IGP + (IED - 1) * (NGP - 1)
               AGP = GET_AGP (IED, INW, NW_IPAR)
               CALL EDGE_VELOCITY
     &              (JGP, AGP, CRD, V_GRID, U_A (1, IED))
            END DO
C           make the tangential velocity be linear, starting at
C           U_A (1, 1) and ending at U_A (1, 2)
            DVS = (U_A (1, 2) - U_A (1, 1)) / DBLE (NGP - 1)
            DO JGP = IGP, IGP + NGP - 1
               V_GRID (1, JGP) = U_A (1, 1) + DVS * DBLE (JGP - IGP)
            END DO
         END IF
         IGP = IGP + NGP
      END DO
C
      RETURN
      END
      SUBROUTINE EDGE_VELOCITY
     &           (IGP, AGP, CRD, VEL, U_A)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'net_funcs.inc'
C
      INTEGER(KIND=IK) IGP
      INTEGER(KIND=IK) AGP
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) VEL (4, *)
      REAL   (KIND=RK) U_A (2)
C
      REAL   (KIND=RK) X_XI, Z_XI, J_INV, U_G (2)
C     GET Velocities of adjacent network in global coordinates
      X_XI    = CRD (3, AGP)
      Z_XI    = CRD (4, AGP)
      J_INV   = 1.0D0 / SQRT (X_XI ** 2 + Z_XI ** 2)
      U_G (1) = (X_XI * VEL (1, AGP) - Z_XI * VEL (2, AGP)) * J_INV
      U_G (2) = (X_XI * VEL (2, AGP) + Z_XI * VEL (1, AGP)) * J_INV
C     And transform them to the local coordinates
      X_XI    = CRD (3, IGP)
      Z_XI    = CRD (4, IGP)
      J_INV   = 1.0D0 / SQRT (X_XI ** 2 + Z_XI ** 2)
      U_A (1) = (X_XI * U_G (1) + Z_XI * U_G (2)) * J_INV
      U_A (2) = (X_XI * U_G (2) - Z_XI * U_G (1)) * J_INV
C
      RETURN
      END
      SUBROUTINE TO_MATERIAL
     &           (NGP, CRD, PHI, PHN, V_GRID, DPHI_DT)
C ---------------------------------------------------------------------------
C     Converts the partial derivatives to the material derivatives and
C     converts the grid velocity from local to global coordinates.
C
C     D/Dt = d/dt + dx/dt * Grad
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
      INCLUDE 'net_params.inc'
      INCLUDE 'tme_funcs.inc'
C
      INTEGER(KIND=IK) NGP
      REAL   (KIND=RK) CRD (4, *)
      REAL   (KIND=RK) PHI (2, *)
      REAL   (KIND=RK) PHN (2, *)
      REAL   (KIND=RK) V_GRID  (4, *)
      REAL   (KIND=RK) DPHI_DT (2, *)
C
      INTEGER(KIND=IK) IGP
      REAL   (KIND=RK) X_XI, Z_XI, J_INV, PHI_S, PHI_N, V_S, V_N
C
      DO IGP = 1, NGP
         X_XI  = CRD (3, IGP)
         Z_XI  = CRD (4, IGP)
         J_INV = 1.0D0 / SQRT (X_XI ** 2 + Z_XI ** 2)
C        Grad Phi in local coordinates
         PHI_S = PHI (2, IGP) * J_INV
         PHI_N = PHN (1, IGP)
C        D/DT (Phi)= d/dT (Phi) + V_grid * Grad (Phi)
         DPHI_DT (1, IGP) = DPHI_DT (1, IGP)
     &                    + V_GRID  (1, IGP) * PHI_S
     &                    + V_GRID  (2, IGP) * PHI_N
C        local grid velocity
         V_S = V_GRID (1, IGP)
         V_N = V_GRID (2, IGP)
C        Convert to global velocity
         V_GRID (1, IGP)  = (X_XI * V_S - Z_XI * V_N) * J_INV
         V_GRID (2, IGP)  = (X_XI * V_N + Z_XI * V_S) * J_INV
C        WRITE (*, 1001) IGP, V_GRID (1, IGP), V_GRID (2, IGP)
      END DO
C
C1001 FORMAT (1X, I3, ':', 10E14.6)
      RETURN
      END

max_stars_repo_name: meggermo/waterwaves

max_issues_repo_name: meggermo/waterwaves

max_forks_repo_name: meggermo/waterwaves


text:
! wind_tm_event_def.for - data structures and declarations for working with
! WIND/WAVES telemetry events

! NOTE:  This file must be "included" after wind_tm_user_def.for in the
! wind_tm_lib.for source because of the definitions for max_channels
! and size_of_packet.
$IF ABSOFT_FORTRAN
!	integer*4       size_of_packet, max_channels, size_of_packet
!	parameter	(size_of_packet=431)
!	parameter	(max_channels=5)
!	parameter	(size_of_packet=431)

	integer*4       max_packets_per_event
	integer*4       perfect_event
	integer*4       truncated_event
	integer*4       gapped_event
	parameter	(max_packets_per_event=32)
	parameter	(perfect_event=1)
	parameter	(truncated_event=3)
	parameter	(gapped_event=5)
$ELSE
!	parameter	max_channels=5
!	parameter	size_of_packet=431

	parameter	max_packets_per_event=32
	parameter	perfect_event=1
	parameter	truncated_event=3
	parameter	gapped_event=5
$ENDIF

	! this is the event built by the event builder and operated on
	! by wind_tm_get_item

$IF ABSOFT_FORTRAN
	integer*4       max_len_h1, max_len_h2
	integer*4       max_len_h3, max_len_h4
	integer*4       max_len_hk
	parameter	(max_len_h1 = 4)
	parameter	(max_len_h2 = 72)
	parameter	(max_len_h3 = 72)
	parameter	(max_len_h4 = 32)
	parameter	(max_len_hk = 81+25+125+5)
$ELSE
	parameter	max_len_h1 = 4
	parameter	max_len_h2 = 72
	parameter	max_len_h3 = 72
	parameter	max_len_h4 = 32
	parameter	max_len_hk = 81+25+125+5
$ENDIF

	structure	/wind_packet_headers/
	   byte		h1(max_len_h1)			! primary header
	   byte		h2(max_len_h2)			! instrument header
	   byte		h3(max_len_h3)			! measurement header
	   byte		h4(max_len_h4)			! packet header
	   logical*4	got_1st
	   logical*4	got_2nd
	   logical*4	got_3rd
	   logical*4	got_4th
	   logical*4	got_data
	end structure
$IF ABSOFT_FORTRAN
	integer*4	sz_wind_packet_headers
	parameter	(sz_wind_packet_headers=max_len_h1+max_len_h2+
	1		 max_len_h3+max_len_h4+ (5*4))
	integer*4	sz_max_event_data
	parameter	(sz_max_event_data=max_packets_per_event *
	1                size_of_packet)
$ELSE
	parameter	sz_wind_packet_headers=max_len_h1+max_len_h2+
	1		max_len_h3+max_len_h4+ (5*4)
	parameter	sz_max_event_data=max_packets_per_event*size_of_packet
$ENDIF

	structure	/wind_event_buffer/
	   integer*4	event_completion_status
	   integer*4	alignment_dummy
	   union
	    map
	      integer*4	cpn	 		! current_packet_number
	      integer*4	num_packets_to_build_this_event
	      record	/wind_packet_headers/ ph(max_packets_per_event)
	      byte		data(sz_max_event_data)
	      integer*4	last_data_byte
	    end map
	    map
	      byte	major_frame(0:sizeof_major_frame-1)
	    end map
	    map
	      byte	minor_frame(0:sizeof_minor_frame-1)
	    end map
	    map
	      byte	packet(0:size_of_packet-1)
	    end map
	   end union
	   integer*4	hk_major		! MF of hk data assoc. w/ event
	   byte		hk(max_len_hk)		! array of hk data
	end structure
	! declare the event buffer, eb
	record		/wind_event_buffer/ eb(max_channels)
	common		/wind_user_blk1/ eb
$IF ABSOFT_FORTRAN
	integer*4	sizeof_eb
	parameter	(sizeof_eb= 4 + 4 +
	1		sizeof_major_frame +
	1		 4 + max_len_hk)
$ELSE
	parameter	sizeof_eb= 4 + 4 +
	1		sizeof_major_frame +
	1		 4 + max_len_hk

$ENDIF
	! this is the internal structure filled by wind_tm_get_event
	! in order to extract the event from the packet(s)

	structure	/wind_event_extract_info/
	   union
	   map
	   character*40	event_type		! event name
	   end map
	   map
	   byte		event_type_b
	   end map
	   map
	   real*8	alignment_dummy
	   end map
	   end union
	   integer*4	event_subtype
	   integer*4	event_status
	   integer*4	dbms_date_time(2)
	   integer*4	last_date_time(2)
	   integer*4	h1_start_bit
	   integer*4	h1_bit_size
	   integer*4	h2_bit_size
	   integer*4	h3_bit_size
	   integer*4	h4_bit_size

	   ! packet id information, constant for every packet in event
	   integer*4	packet_id_start_bit
	   integer*4	packet_id_bit_size
	   integer*4	packet_id

	   integer*4	h1_start_byte,   h1_byte_size,   h1_end_byte
	   integer*4	h2_start_byte,   h2_byte_size,   h2_end_byte
	   integer*4	h3_start_byte,   h3_byte_size,   h3_end_byte
	   integer*4	h4_start_byte,   h4_byte_size,   h4_end_byte
	   integer*4	data_start_byte, data_byte_size, data_end_byte
	   integer*4	tds_start_byte,  tds_byte_size,  tds_end_byte

	   ! flag bits
	   integer*4	first_packet_mask
	   integer*4	last_packet_mask

	   ! values calculated from first packet in event,
	   ! constant for every packet in event
	   integer*4	got_a_first_packet
	   integer*4	got_a_last_packet

	end structure
	record		/wind_event_extract_info/
	1		eei(max_channels)	! event_extract_info
	common		/wind_user_blk2/ eei
$IF ABSOFT_FORTRAN
	integer*4	sizeof_eei
	parameter	(sizeof_eei= 40 + (4*2) + (8*2) + (4*5) +
	1		(4*3) + (4*6*3) + (4*2) + (4*2))
$ELSE
	parameter	sizeof_eei= 40 + (4*2) + (8*2) + (4*5) +
	1		(4*3) + (4*6*3) + (4*2) + (4*2)
$ENDIF

ext: for


text:
program v
integer(kind=8) ::d = 9999999999_8
end program v

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
C    %W% %G%
      subroutine postbrn(record)

      include 'ipfinc/parametr.inc'

      include 'ipfinc/postdta.inc'
      include 'ipfinc/bscords.inc'
      include 'ipfinc/pasprm.inc'

      character record*(*), type*1, subtype*1
      character rcrdpar*300 ! record for parallel circuits
      character floa*25, flob*25, fmtflo*25, fmtflof*25, par*7
      character ratcod*1, ratcod1*1, ovldflg*3, fmtovld*20, ovldval*30
      character crcodam*1, crcodva*1
      character fmtcmp*15, pctcomp*7
      character begpath*30
      character parcord(11)*120 ! coordinate records for parallel circuits
      character in_buffer * (MAXBUFFER)
      character out_pbuf * (MAXBUFFER)
      character null*1, linefeed*1
      character fmttap*16, tapa1*10, tapa2*10  !tx taps
      character ckt*1             ! Circuit - used for separate parallels

      integer error, firstp

      data fmtflo /'(i9)'/
      data fmtflof /'(f9.2)'/
      data fmtovld/'(''('',a1,1x,i9,'')'')'/
      data fmtcmp /'(''('',i9,''%)'')'/
      data fmttap /'(''('',f9.2,'')'')'/
      data lenbrec /257/       ! length of branch record + null

      null = char(0)
      linefeed = char(10)
      lparsw = 0               ! default - processing parallels combined
********************* debug stuff **********************
*     print 9901, record
*9901 format(' entering postbrn record=:',a)
********************* debug stuff **********************
      read(record,101)
     & type, subtype, buspf1, basepf1, buspf2, basepf2, ckt, numpar,
     & pin, qin, pout, qout, ploss, qloss,
     & crfloam, crratam, crcodam, mxendam,
     & crflova, crratva, crcodva, mxendva,
     & ttpctam, ttfloam, ttpctva, ttflova, par
  101 format(bz, 2a1, t7, a8, f4.0, 1x, a8, f4.0, a1, i1, 6f15.7,
     & 2(e15.7, f8.1, a1, i1), 4e15.7, t250, a7)
********************* debug stuff **********************
*     print 9901, record
*     print 9901, nxtcord
********************* debug stuff **********************

      if (dgoptyp .eq. 5) then    ! If coordinate only, set bus powerflow
          type = nxtcord(1:1)     ! identification to coordinate id.
          read (nxtcord, 226)  buspf1, basepf1, buspf2, basepf2
      endif

      indcor2 = ksrchcor(buspf2, basepf2)
********************* debug stuff **********************
*     print 9931, buspf1, basepf1,buspf2, basepf2, indcor2
*9931 format(' buspf2, basepf2, indcor2 :', 2(a8,f6.1),i5)
*     if (subtype .eq. 'P') print 9933, record(234:249)
*9933 format (a)
********************* debug stuff **********************

*********** discard lines that are not of interest **************
      if (indcor2 .le. 0) go to 400
      if ( buspf1 .gt. buspf2 .or.
     &    (buspf1 .eq. buspf2 .and. basepf1 .gt. basepf2)) go to 1000
*********** lines that are not of interest have been discarded ********
      ibasepf1 = basepf1
      ibasepf2 = basepf2
      karosg = 1       ! section number that arrow is in

******* build the coordinate track for a line ******************
      indcor1 = ksrchcor(buspf1, basepf1)
********************* debug stuff **********************
*     print 9951, buspf1, basepf1, indcor1
*9951 format(' buspf1, basepf1, indcor1', a8,f6.1,i5)
********************* debug stuff **********************
      if (indcor1 .le. 0) go to 1000
      begpath = '[ '
      if (kprflg(indcor1) .eq. 0)
     & write(begpath(3:30), 201) xcord(indcor1), ycord(indcor1)
  201 format(2(f6.2, ' cmtr '))
      write (lpost,211) begpath
  211 format(a)

******* look for branch coordinate data

****** get next coordinate record *******************
  220 continue
********************* debug stuff **********************
*     print 9941, kcorflg
*9941 format('kcorflg = :', i2)
********************* debug stuff **********************

      if (kcorflg .eq. 0) then
          read (corfil, 221, end=250) nxtcord
  221     format(a)
          kcorflg = 1
      elseif (kcorflg .eq. 1) then
****** information is already in nxtcord *******************
          continue
      elseif (kcorflg .eq. 2) then
          go to 280
      endif
********************* debug stuff **********************
*     print 9961, nxtcord
*9961 format('nxtcord = :', a)
********************* debug stuff **********************

      if ( nxtcord(1:1) .eq. 'L' .or. nxtcord(1:1) .eq. 'T') then
         read (nxtcord, 226) buscd1, basecd1, buscd2, basecd2
  226    format(bz, t3, 2(a8, f4.0) )

         if     (buspf1 .gt. buscd1) then
             kcorflg = 0
             go to 220
         elseif (buspf1  .eq. buscd1 .and.
     1           basepf1 .gt. basecd1) then
              kcorflg = 0
              go to 220
         elseif (buspf1 .eq. buscd1 .and. basepf1 .eq. basecd1 .and.
     1           buspf2 .gt. buscd2) then
              kcorflg = 0
              go to 220
         elseif (buspf1 .eq. buscd1 .and. basepf1 .eq. basecd1 .and.
     &           buspf2 .eq. buscd2 .and. basepf2 .gt. basecd2) then
              kcorflg = 0
              go to 220

         elseif (buspf1 .eq. buscd1 .and. basepf1 .eq. basecd1 .and.
     &           buspf2 .eq. buscd2 .and. basepf2 .eq. basecd2) then
              kcorflg = 0

******* begin set up for plotting parallels separately *********
               if (broppar .eq. 1 .and.
     &            (numpar .gt. 1 .or. vaopdif .eq. 1)) then
               iparcor = 1
               parcord(iparcor) = nxtcord
  230          if (iparcor .lt. 11) iparcor = iparcor + 1
               read (corfil, 221, end=232) parcord(iparcor)
               if ( parcord(iparcor)(1:1) .ne. 'L' .and.
     &              parcord(iparcor)(1:1) .ne. 'T') then
                  kcorflg = 3
                  go to 234
               endif

               read (parcord(iparcor), 226)
     &          buscd1, basecd1, buscd2, basecd2
               if (buspf1 .eq. buscd1 .and. basepf1 .eq. basecd1 .and.
     &             buspf2  .eq. buscd2 .and. basepf2 .eq. basecd2)
     &           go to 230
               kcorflg = 3
               go to 234
  232          kcorflg = 2
  234          iparcor = iparcor - 1
***************** debug stuff **********************
*     print 9234, iparcor
*9234 format(' postbrn at 234-iparcor: ',i5)
***************** debug stuff **********************
               if (iparcor .gt. 1) then
                  write(in_buffer, 237) record(1:31)
  237             format(a)
                  lastp = lastch(in_buffer)
                  in_buffer(lastp+1:) = null
                  out_pbuf(1:) = null
************* debug stuff ***********************
*     print 9237, in_buffer(1:100)
*9237 format (' in_buffer at 237', a)
************* debug stuff ***********************
                  if (vaopdif .eq. 0) then
                     call gtparout(in_buffer, out_pbuf, error)

************* debug stuff ***********************
*      print 9235
* 9235 format (' Returned from gtparout')
************* debug stuff ***********************
                  else
                     call pp_odif(in_buffer, out_pbuf, error)
************* debug stuff ***********************
*      print 9236
* 9236 format (' Returned from pp_odif')
************* debug stuff ***********************
                  endif

                  jparcor = 0
                  firstp = 1
                  lastp = index(out_pbuf, null)
                  nparpf = (lastp - firstp + 1)/lenbrec

************* debug stuff ***********************
**** This debug block gave the run time error: 'output statement overflows record'
*      idbg1 = 1
*      do 9240 idbg = 1, nparpf
*      idbg2 = idbg1 + lenbrec - 3
*      idbg2 = idbg1 + 255
*      print 9239, lenbrec, lastp, nparpf, iparcor
* 9239 format(' lenbrec, lastp, nparpf, iparcor ',4i5)
*      print 9238, out_pbuf(idbg1:idbg2)
* 9238 format (' out_pbuf at 237+', a)
*      idbg1 = idbg1 + lenbrc
* 9240 continue
************* debug stuff ***********************

                  if (nparpf .le. iparcor) then
                     lparsw = 1
                     go to 260
                  endif
               endif
            endif
******* end   set up for plotting parallels separately *********
*********** debug stuff ******************
*     print 9241, nxtcord
*9241 format(' in postbrn before 241-nxtcord',a)
*********** debug stuff ******************

            read (nxtcord, 241) karosg
  241       format(bz, t29, i2)
            if (karosg .eq. 0) karosg = 1
*           kcorflg = 0
            ind1 = 31
            max = 90

            do while (ind1 .lt. max)
              if (nxtcord(ind1:ind1+11) .eq. ' ') go to 280
              read (nxtcord(ind1:ind1+11), 246) xcorda, ycorda
  246         format(bz, 2f6.2)
              write(lpost,201) xcorda, ycorda
              ind1 = ind1 + 12
            end do
         else
            kcorflg = 1
         endif

      endif
      go to 280
  250 kcorflg = 2
*************** begin retrieve coordinate data for each parallel ********
  260 jparcor = jparcor + 1
      if ( firstp .ge. lastp) then
         print 261
  261    format(' program error in plotting parallel circuits ')
         return
      else
         nextp = nxt_term(out_pbuf(firstp+1:)) + firstp
**************** debug stuff *********************
*      print 9261, firstp, nextp,  lastp, iparcor, jparcor
* 9261 format(' after 261 firstp, nextp,  lastp, iparcor, jparcor ', 5i5)
*      print 9263, out_pbuf(firstp:firstp+32)
* 9263 format(' out_pbuf(firstp:firstp+32)', a)
*      print 9264, out_pbuf(firstp+33:firstp+122)
* 9264 format(' out_pbuf(firstp+33:firstp+122)', a)
*
*      print 9265, out_pbuf(firstp+233:firstp+240)
* 9265 format(' after 261 out_pbuf(firstp+233:firstp+240)', a)
**************** debug stuff *********************
         rcrdpar = out_pbuf(firstp:nextp-1)
         read (out_pbuf(firstp:nextp-1),101)
     & type, subtype, buspf1, basepf1, buspf2, basepf2, ckt, numpar,
     & pin, qin, pout, qout, ploss, qloss,
     & crfloam, crratam, crcodam, mxendam,
     & crflova, crratva, crcodva, mxendva,
     & ttpctam, ttfloam, ttpctva, ttflova, par
         firstp = nextp
         if (out_pbuf(firstp:firstp) .eq. linefeed) firstp = firstp + 1

**************** debug stuff *********************
*     print 9262, parcord(jparcor)
*9262 format(' parcord(jparcor) ',a)
**************** debug stuff *********************
         read (parcord(jparcor), 241) karosg
         if (karosg .eq. 0) karosg = 1
*        kcorflg = 0
         ind1 = 31
         max = 90

         do while (ind1 .lt. max)
            if (parcord(jparcor)(ind1:ind1+11) .eq. ' ') go to 280
            read (parcord(jparcor)(ind1:ind1+11), 246) xcorda, ycorda
            write(lpost,201) xcorda, ycorda
            ind1 = ind1 + 12
         end do
      endif
*************** begin retrieve coordinate data for each parallel ********

  280 continue
      if (kprflg(indcor2) .eq. 0)
     1 write(lpost,201) xcord(indcor2), ycord(indcor2)

      write (lpost, 291)
  291 format (' ]')

****** build remainder of parameters for call to acLine *********
      ipinflg =  sign (1.0, pin)
****** build flo for PQ, MVA/I, LOSS *************
      if     (dgoptyp .eq. 1 .or. dgoptyp .eq. 4) then
****** build flo for PQ *************
****** or for Interchange if there is a branch in the interchange coord file
         floa = '()'
         flob = '()'
****** build flo for PQ -  Pin, Qin, Pout, Qout ***************
         if (ipinflg .lt. 0) then
            ptmp = pin
            qtmp = qin
            pin = - pout
            qin = - qout
            pout = - ptmp
            qout = - qtmp
         endif
************* begin formatting sending end P & Q *************
         if (broppsn .eq. 1 .or. bropqsn .eq. 1) then
            nxtflo = 3            ! next flow data starts here
         floa = '( '

            if (broppsn .eq. 1) then
               apin = abs(pin) + 0.5
               ipin = sign (apin, pin)
               lognum = alog10(amax1(apin, 1.0)) + 1.0
               if (ipin .lt. 0) lognum = lognum + 1
********************* debug stuff **********************
*     print 9971, lognum, fmtflo
*9971 format(' in postbrn at 9971, lognum, fmtflo = :',i5,1x,a)
********************* debug stuff **********************
               write (fmtflo(3:3), 311) lognum
  311          format(i1)
               lstflo =  nxtflo + lognum  - 1   ! last character of next flow data
********************* debug stuff **********************
*     print 9911, fmtflo
*9911 format(' in postbrn, fmtflo = :',a)
********************* debug stuff **********************
               write (floa(nxtflo:lstflo), fmtflo) ipin
********************* debug stuff **********************
*     print 9921, floa
*9921 format(' in postbrn, flo = :',a)
********************* debug stuff **********************
               nxtflo = lstflo + 1
               lstflo =  nxtflo
               floa(nxtflo:lstflo) = ' '
               nxtflo = nxtflo + 1
            endif

            if (bropqsn .eq. 1) then
               aqin = abs(qin) + 0.5
               iqin = sign (aqin, qin)
               lognum = alog10(amax1(aqin, 1.0)) + 1.0
               if (iqin .lt. 0) lognum = lognum + 1
               write (fmtflo(3:3), 311) lognum

               lstflo =  nxtflo
               floa(nxtflo:lstflo) = '('
               nxtflo = nxtflo + 1

               lstflo =  nxtflo + lognum  - 1   ! last character of next flow data
               write (floa(nxtflo:lstflo), fmtflo) iqin
               nxtflo = lstflo + 1
               lstflo =  nxtflo
               floa(nxtflo:lstflo) = ' '
               nxtflo = nxtflo + 1

               lstflo =  nxtflo
               floa(nxtflo:lstflo) = ')'
               nxtflo = nxtflo + 1
            endif

            floa(nxtflo:nxtflo+1) = ' )'
         endif

************* begin to format receiving end P & Q *********

         if (bropprc .eq. 1 .or. bropqrc .eq. 1) then
            nxtflo = 4            ! next flow data starts here
            flob = '( ['

            if (bropprc .eq. 1) then
               apout = abs(pout) + 0.5
               ipout = sign (apout, pout)
               lognum = alog10(amax1(apout, 1.0)) + 1.0
               if (ipout .lt. 0) lognum = lognum + 1
********************* debug stuff **********************
*     print 9971, lognum, fmtflo
********************* debug stuff **********************
               write (fmtflo(3:3), 311) lognum
               lstflo =  nxtflo + lognum  - 1   ! last character of next flow data
********************* debug stuff **********************
*     print 9911, fmtflo
********************* debug stuff **********************
               write (flob(nxtflo:lstflo), fmtflo) ipout
********************* debug stuff **********************
*     print 9921, flob
********************* debug stuff **********************
               nxtflo = lstflo + 1
               lstflo =  nxtflo
               flob(nxtflo:lstflo) = ' '
               nxtflo = nxtflo + 1
            endif

            if (bropqrc .eq. 1) then
               aqout = abs(qout) + 0.5
               iqout = sign (aqout, qout)
               lognum = alog10(amax1(aqout, 1.0)) + 1.0
               if (iqout .lt. 0) lognum = lognum + 1
********************* debug stuff **********************
*     print 9971, lognum, fmtflo
********************* debug stuff **********************
               write (fmtflo(3:3), 311) lognum
********************* debug stuff **********************
*     print 9911, fmtflo
********************* debug stuff **********************
               lstflo =  nxtflo
               flob(nxtflo:lstflo) = '('
               nxtflo = nxtflo + 1

               lstflo =  nxtflo + lognum  - 1   ! last character of next flow data
               write (flob(nxtflo:lstflo), fmtflo) iqout
********************* debug stuff **********************
*     print 9921, flob
********************* debug stuff **********************

               nxtflo = lstflo + 1
               lstflo =  nxtflo
               flob(nxtflo:lstflo) = ' '
               nxtflo = nxtflo + 1
               lstflo =  nxtflo
               flob(nxtflo:lstflo) = ')'
               nxtflo = nxtflo + 1
            endif

            flob(nxtflo:nxtflo+2) = '] )'
         endif

      elseif (dgoptyp .eq. 2) then
****** build flo for MVA and I
         flob = '()'
         nxtflo = 3            ! next flow data starts here
         floa = '(<'

*************** build current component of flow ******************
         if (mxendam .ne. 0) then
            pin = ttfloam
            apin = abs(pin) + 0.5
            ipin = sign (apin, pin)
            lognum = alog10(amax1(apin, 1.0)) + 1.0
            if (ipin .lt. 0) lognum = lognum + 1
********************* debug stuff **********************
*     print 9971, lognum, fmtflo
********************* debug stuff **********************
            write (fmtflo(3:3), 311) lognum
            lstflo =  nxtflo + lognum  - 1   ! last character of next flow data
********************* debug stuff **********************
*     print 9911, fmtflo
********************* debug stuff **********************
            write (floa(nxtflo:lstflo), fmtflo) ipin
********************* debug stuff **********************
*     print 9921, floa
*     print 9981, mxendva, ttflova
********************* debug stuff **********************
            nxtflo = lstflo + 1
            if (mxendva .ne. 0) then
              floa(nxtflo:nxtflo) = 'a'
              nxtflo = nxtflo + 1
            endif
         endif

*************** build MVA component of flow ***
         if (mxendva .ne. 0) then
           pin = ttflova
           apin = abs(pin) + 0.5
           ipin = sign (apin, pin)
           lognum = alog10(amax1(apin, 1.0)) + 1.0
           if (ipin .lt. 0) lognum = lognum + 1
********************* debug stuff **********************
*     print 9971, lognum, fmtflo
********************* debug stuff **********************
           write (fmtflo(3:3), 311) lognum
           lstflo =  nxtflo + lognum  - 1   ! last character of next flow data
********************* debug stuff **********************
*     print 9911, fmtflo
********************* debug stuff **********************
           write (floa(nxtflo:lstflo), fmtflo) ipin
********************* debug stuff **********************
*     print 9921, floa
*     print 9981, mxendva, ttflova
*9981 format('mxendva, ttflova = ',i3, f9.1)
********************* debug stuff **********************
           nxtflo = lstflo + 1
         endif
*** The line component overrides tx component in determining whether maxflow
*** is at the near or far end of mised L & T sections
         if (mxendam .ne. 0) then
           mvaend = mxendam
         else
           mvaend = mxendva
         endif
         if ((ipinflg .ge. 0 .and. mvaend .eq. 2) .or.
     1     (ipinflg .lt. 0 .and. mvaend .eq. 1)) then
            floa(nxtflo:nxtflo) = 'R'
            nxtflo = nxtflo + 1
         endif
         floa(nxtflo:nxtflo+1) = '>)'

      elseif (dgoptyp .eq. 3) then
****** build flo for Losses

         floa = '()'
         flob = '()'
         pin = ploss
         qin = qloss
************* begin formatting  P & Q losses *************
************* use broppsn and bropqsn  to select P and/or Q losses ****
         if (broppsn .eq. 1 .or. bropqsn .eq. 1) then
            nxtflo = 2            ! next flow data starts here
         floa = '('

            if (broppsn .eq. 1) then
               apin = abs(pin) + 0.5
               ipin = sign (apin, pin)
               lognum = alog10(amax1(apin, 1.0)) + 4.0
               if (ipin .lt. 0) lognum = lognum + 1
********************* debug stuff **********************
*     print 9971, lognum, fmtflof
********************* debug stuff **********************
               write (fmtflof(3:3), 311) lognum
               lstflo =  nxtflo + lognum  - 1   ! last character of next flow data
********************* debug stuff **********************
*     print 9911, fmtflof
********************* debug stuff **********************
               write (floa(nxtflo:lstflo), fmtflof) pin
********************* debug stuff **********************
*     print 9921, floa
********************* debug stuff **********************
               nxtflo = lstflo + 1
               lstflo =  nxtflo
               floa(nxtflo:lstflo) = ' '
               nxtflo = nxtflo + 1
            endif

            if (bropqsn .eq. 1) then
               aqin = abs(qin) + 0.5
               iqin = sign (aqin, qin)
               lognum = alog10(amax1(aqin, 1.0)) + 4.0
               if (iqin .lt. 0) lognum = lognum + 1
********************* debug stuff **********************
*     print 9971, lognum, fmtflof
********************* debug stuff **********************
               write (fmtflof(3:3), 311) lognum

               lstflo =  nxtflo
               floa(nxtflo:lstflo) = '('
               nxtflo = nxtflo + 1

               lstflo =  nxtflo + lognum  - 1   ! last character of next flow data
********************* debug stuff **********************
*     print 9911, fmtflof
********************* debug stuff **********************
               write (floa(nxtflo:lstflo), fmtflof) qin
********************* debug stuff **********************
*     print 9921, floa
********************* debug stuff **********************
               nxtflo = lstflo + 1
*              lstflo =  nxtflo
*              floa(nxtflo:lstflo) = ' '
*              nxtflo = nxtflo + 1

               lstflo =  nxtflo
               floa(nxtflo:lstflo) = ')'
               nxtflo = nxtflo + 1
            endif

            floa(nxtflo:nxtflo) = ')'
         endif

         continue

      elseif (dgoptyp .eq. 5) then
****** build flo for Coordinates only
         floa = '()'
         flob = '()'
********************* debug stuff **********************
*     print 9301, floa
*9301 format(' before 315 in postbrn floa = ',a)
********************* debug stuff **********************

      elseif (dgoptyp .eq. 6) then
****** build flo for % loading
         floa = '(    %)'
         flob = '()'
         if (mxendva .ne. 0) then
           olpct = ttpctva
           if (mxendam .ne. 0 .and. ttpctam .gt. ttpctva) then
             olpct = ttpctam
             floa = '(    %a)'
           endif
         else
           olpct = ttpctam
         endif
         if (abs(olpct) .le. 999.0) then
             aolpct = abs(olpct) + 0.5
             iolpct = sign (aolpct, olpct)
             write(floa(2:5),315) iolpct
  315        format(i4)
         else
             floa = '(****%)'
         endif
      endif

********** flo is complete ****************************
      if (vaopdif .eq. 0) then
         if (lparsw .eq. 1 )  then
             write (par, 320) ckt
  320        format('(#', a1, ')' )
         elseif (numpar .le. 1 )  then
*            if (numpar .le. 1 )  then
             par = '()'
         else
             write (par, 321) numpar
  321        format('([', i1, '])' )
         endif
      endif

********* build branch overload information *********
      ovldflg = '()'
      ovldval = '()'
*** check for overloads except for coords only of difference plots
      if (dgoptyp .ne. 5 .and. vaopdif .eq. 0) then
*** calculate % loading for worst parallel
        if (mxendam .ne. 0 .and. crratam .gt. 0.0) then
          crpctam = crfloam / crratam * 100.0
        else
          crpctam = 0.0
        endif
        if (mxendva .ne. 0 .and. crratva .gt. 0.0) then
          crpctva = crflova / crratva * 100.0
        else
          crpctva = 0.0
        endif

        if (mxendva .ne. 0) then
          if (mxendam .ne. 0) then
***         checking mixed L & T sections in circuit
            ratcod = '*'
            if (crpctva .gt. crpctam) then
***           Transformer is more overloaded
                ratcod1 = 'm'
                crflo = crflova
                crrat = crratva
                olpct = crpctva
            else
***           Line is more overloaded
                ratcod1 = 'a'
                crflo = crfloam
                crrat = crratam
                olpct = crpctam
            endif
          else       ! (mxendam .ne. 0) test is false
***         checking for parallel transformer circuits
            if (ttflova .gt. crflova) then
              ratcod = '*'
            else
              ratcod = crcodva
            endif
            ratcod1 = 'm'
            crflo = crflova
            crrat = crratva
            olpct = crpctva
          endif
        else         !(mxendva .ne. 0) test is false
***       checking for parallel line circuits
          if (ttfloam .gt. crfloam) then
            ratcod = '*'
          else
            ratcod = crcodam
          endif
          ratcod1 = 'a'
          crflo = crfloam
          crrat = crratam
          olpct = crpctam
        endif        !(mxendva .ne. 0) test is complete

        if (olpct .ge. 100.0) then
          ovldflg = '('//ratcod//')'
        endif
***********debug stuff ***************************
*     print 9311, olpct
*9311 format (' in postbrn at 321-olpct = ',f8.2)
***********debug stuff ***************************
******************************* temp debug stuff*************
        if (olpct .ge. 90.0) then
*       if (olpct .ge. 10.0) then
******************************* temp debug stuff*************

          ovldval = '('
          nxtflo = 2
          pin = crflo
          apin = abs(pin) + 0.5
          ipin = sign (apin, pin)
          lognum = alog10(amax1(apin, 1.0)) + 1.0
          if (ipin .lt. 0) lognum = lognum + 1
********************* debug stuff **********************
*     print 9971, lognum, fmtflo
********************* debug stuff **********************
          write (fmtflo(3:3), 311) lognum
          lstflo =  nxtflo + lognum  - 1   ! last character of next flow data
********************* debug stuff **********************
*     print 9911, fmtflo
********************* debug stuff **********************
          write (ovldval(nxtflo:lstflo), fmtflo) ipin
********************* debug stuff **********************
*     print 9921, ovldval
********************* debug stuff **********************
          nxtflo = lstflo + 1
          if (mxendva .ne. 0 .and. mxendam .ne. 0) then
            ovldval(nxtflo:nxtflo) = 'a'
            nxtflo = nxtflo + 1
          endif

          ovldval(nxtflo:nxtflo) = ratcod
          kolval = crrat + 0.5
          olvaln = kolval
          lognum = alog10 (amax1(olvaln, 1.0)) + 1
          nxtflo = nxtflo + 1
          lstflo = nxtflo + lognum - 1
          write (fmtflo(3:3), 311) lognum
          write (ovldval(nxtflo:lstflo),fmtflo) kolval
          nxtflo = lstflo + 1
          ovldval(nxtflo:nxtflo) = ')'
        endif    ! end (olpct .ge. 90.0) test
      endif      ! end (dgoptyp .ne. 5 .and. vaopdif .eq. 0) test

      if (type .eq. 'L' .or. type .eq. 'E') then
*************** check for compensation ***************
********** debug stuff ******************
*     print 9341, lparsw
*9341 format(' before 341 lparsw = ',i5)
********** debug stuff ******************
         if (lparsw .eq. 0) then
            read (record,341)  pctcmp
         else
            read (rcrdpar,341) pctcmp
         endif
  341    format(t234, f8.2)
         if (abs(pctcmp) .gt. 0.1) then
            apctcmp = abs(pctcmp) + 0.5
            ipctcmp = sign(apctcmp, pctcmp)
            lognum = alog10 (amax1(apctcmp, 1.0)) +1
            if (ipctcmp .lt. 0.0) lognum = lognum + 1
            if (lognum .gt. 4) then
               pctcomp = '(****%)'
            else
               write (fmtcmp(7:7), 311) lognum
               write (pctcomp, fmtcmp) ipctcmp
            endif
         else
            pctcomp = '()'
         endif

*************** write a postscript line record **********
         write (lpost, 351) ovldflg, ovldval, pctcomp, karosg,
     1   ipinflg, floa, flob, par, kotgflg, ibasepf1
  351    format(a3, 1x, a30, 1x, a7, 1x, i2, 1x, i8, 1x, 2(a25, 1x), a7,
     1   1x, i1, 1x, i5, ' acLine')
      elseif (type .eq. 'T') then
         if (lparsw .eq. 0) then
            read (record,355)  tap1, tap2
         else
            read (rcrdpar,355) tap1, tap2
         endif

  355    format(t234, 2f8.1)

         atap = abs(tap1) + 0.005
         itap = sign (atap, tap1)
         lognum = alog10(amax1(atap, 1.0)) + 4.0
         if (itap .lt. 0) lognum = lognum + 1
********************* debug stuff **********************
*     print 9871, lognum, fmttap
*9871 format(' in postbrn at 9871, lognum, fmttap = :',i5,1x,a)
********************* debug stuff **********************
         if (lognum .gt. 8) then
            tapa1 = '(*****.**)'
         else
            write (fmttap(7:7), 311) lognum
********************* debug stuff **********************
*     print 9811, fmttap
*9811 format(' in postbrn, fmttap = :',a)
********************* debug stuff **********************
            write (tapa1, fmttap) tap1
         endif
********************* debug stuff **********************
*     print 9821, tapa1
*9821 format(' in postbrn, tapa1 or tapa2 = :',a)
********************* debug stuff **********************
         if (subtype .eq. 'P') then
           tapa2 = '(DEG PS)'
         else
           atap = abs(tap2) + 0.005
           itap = sign (atap, tap2)
           lognum = alog10(amax1(atap, 1.0)) + 4.0
           if (itap .lt. 0) lognum = lognum + 1
           if (lognum .gt. 8) then
              tapa2 = '(*****.**)'
           else
              write (fmttap(7:7), 311) lognum
              write (tapa2, fmttap) tap2
           endif
         endif
         write (lpost, 361) ovldflg, ovldval, karosg, ipinflg,
     1   floa, flob, par, kotgflg
  361    format(a3, 1x, a30, 1x, i2, 1x, i8, 1x, 2(a25, 1x), a7, 1x
     1   , i1)
         write (lpost, 363) ibasepf1, ibasepf2, tapa1, tapa2
  363    format(1x, 2(1x, i5), 2(1x, a10), ' acTran')

      endif

****** begin write start of next path if processing circuits separately ********
      if (lparsw .eq. 1) then
************** debug stuff ***********************
*     print 9361, jparcor, nparpf
*9361 format(' after 361 jparcor, nparpf ', 2i5)
************** debug stuff ***********************
         if (jparcor .lt. nparpf) then
            write (lpost, 211) begpath
            go to 260
         else
*           if (kcorflg .eq. 3) then
*              kcorflg = 1
*              nxtcord = parcord(iparcor + 1)
*           endif
         endif
      endif
****** end   write start of next path if processing circuits separately ********

      if (kcorflg .eq. 3) then
         kcorflg = 1
         nxtcord = parcord(iparcor + 1)
********** debug stuff ***************************
*      print 9851, iparcor,nxtcord(1:30)
* 9851 format (' before 400: iparcor, nxtcord(1:10) ', i2, a)
********** debug stuff ***************************
      endif

      go to 1000

********* begin latent branch logic *******************
  400 if (bsopsum .eq. 1) then
         psumflo = psumflo + pin
         qsumflo = qsumflo + qin
      endif
********* end   latent branch logic *******************

 1000 continue
********************* debug stuff **********************
*     print 9999
*9999 format(' exiting postbrn')
********************* debug stuff **********************
      return
      end

ext: f

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
So you think youve been overcharged... you just got home, looked at your account online, and that that darn restaurant or gas station or rental company has taken money from your account that is over and above what the bill was, even including tip!  Screw em... time to fire up the trust Davis Wiki and chew them out about it...

Wait a second... you might want to learn a bit about how banking works first...

You see, there are different kinds of charges to accounts, and usually your online banking service doesnt show you the difference between a normal charge (the kind that actually takes money out), and an authorization hold, which (and this is the really important part) does not actually take money out of your account.  Yes, it is showing up on your bank statement right now, and yes, they did do a transaction in excess... but that is because they didnt know how much to charge you.  The extra charge will disappear in a day or two.  Really.  Thats how it is supposed to work.

Why do they do this?  Well, its pretty simple:  you just sat down at a Restaurants restaurant, ate a nice $40 meal, and the waiter hands you the bill.  You hand them your credit card, and they whisk it away.  When they run your card, you havent yet added a Tipping tip.  They run it, the transaction adds a small excess floating amount, and they return your card and the bill.  You then write in your tip, and sometime in the next day they do their paperwork and tell the credit card company the exact amount.  Thats why you dont have two charges on your account every time you tip somebody at a restaurant.

Even more obvious is why gas stations do this:  When you swipe your card at a pump, it is authorized.  At the time, they have no idea how much gas you are getting... are you putting a couple gallons into a scooter, or filling a recreation vehicle up with $250 worth of gas?  They toss an authorization hold on there, and you pump the gas.  Its usually not noticed in this case because the pump is fairly regularly updating the amounts, but you can catch it if you check your account right after pumping gas.  In this case, since they dont have any idea of the amount, they often show up as an undercharge.

Please note that when you pay by debit card without using a PIN, some gas stations put a hold of up to $75 on your card. That means that up to $75 of your checking account can be off limits (in addition to the price you paid for gas). The time of the hold can vary and depends on your bank. (i.e. 72 hours, 5 days, etc.) This could cause overdrafts fees, bounced checks, or insufficient funds on your account.

How can you protect yourself at the filling station?

    If you must use a debit card, pay inside where you can use your PIN number; PINbased transactions are registered immediately.

    If you pay at the pump, use a credit card.

    Oil companies proprietary charge cards often dont have any kind of block feature on them, but not always. Chevron, for example, briefly blocks out $18.

    Cash is still good  many gas stations, tired of paying high creditcard fees, now give a discount for cash. Look for such deals.

ARCO: $75 hold
Chevron: $1$18
Shell: $1

Rental places will sometimes put a hefty hold on there, reserving the right to charge you for expected room service (yeah, statistically, youre likely to use it) or repair or replacement of an item or room.  Hotels, car rentals and other places usually do this.  You also see extra dollar charges sometimes appear to verify accounts.

The cool thing is that the charge automatically goes away if the merchantFootnote(Merchant is banking speak for the restaurant or business) doesnt file within a certain time period.  It also often (depending on your card) doesnt affect your actual available balance.  That means you can have $100 available, go to two restaurants, have $45 meals at both and they charge you $55 each, but you wont have gone over $100 if you tipped less than $5 each time.  Oh, youre going to the fourth level of hell when you dieFootnote(4th Level of Hell: Being forced to work food service in perpetuity.) for skimping on the tip, but your account was never overdrawn.  Like most things involving finance and banking, it depends on your details of your account: debit/charge cards may well get overdrawn, but real credit cards wont.  By the way, ever wonder why rental places sometimes insist on a real credit card?

Authorization holds exist because a merchant can not check to see how much money you have available.  This would be a blatant violation of your privacy, and could be used by unscrupulous individuals to steal from you.  An authorization hold simply asks the credit card company if you have X amount available, and to hold that money for a short period of time.

In short, before writing a scathing review about how you were overcharged, you might want to wait 48 hours and see if that overcharge makes like a baker meeting a Boojum and softly and suddenly vanishes away.

Links for more info

  http://wiki.answers.com/Q/Why_would_a_restaurant_put_a_temporary_charge_on_your_credit_card Wiki Answers Why would a restaurant put a temporary charge on your credit card?
  wiki:wikipedia:Authorization hold Authorization Hold on Wikipedia

Well Yes and No

First off an Authorization hold or PreAuthorization is done to check to ensure the card is active or good. A $1 dollar charge is normally done and is considered one of the Best Practice for merchants accepting credit cards. There are of course different transaction types for different vertical markets (restaurants, hotels, car rental, gas station, etc) each having a little bit of a tweak on transaction types that best fit the needs of these different businesses.

As for the businesses there are merchants who prefer live transactions versus batch transactions. Having a system that charges against your card in real time is expensive and normally reserved for high end stores, boutiques and jewelers. Everyone else uses batch processing. A batch is a file of transactions kept on the POS (point of sale) terminal or the store/company server. At some regular interval, midnight of the business day, 2 days, 3 days, a week. The business transmits the information to THEIR bank and the money is applied against your account. You may have noticed at times that you bought something at a store it showed up as a preauth and then disappeared and the money went back into your account until a couple days later SURPRISE!! money is gone. This is the reason why.

And while the first commenter makes some good points Im in disagreement. Tonight on the way home from SFO I stopped in Vacaville at a Chevron. When I got home I decided to review my accounts since I had been on vacation for a week. I was shocked when I saw the amount of gas I pumped 11.86 and then a second transaction for $126.00 dollars. Yes your reading that correctly one hundred twenty six dollars. I called my credit card company straight away and challenged the transaction. While I understand this is a business, other business like hotels, car rental agencies explain in advance that they will make a preauth for the total amount that will linger for 23 days. Im not sure if this is required for them or just part of their best practices. Im outraged about this preauth and dont have the attitude that is ok or that we should just take it. You can fight city hall and win. Here are some of my suggestions:

                  Shop / Get your gas elsewhere, if dont like the service go to their competition
                  Write your congressman / congresswoman asking for preauth amounts to be stated before you make a purchase
                  Write a letter to parent company of the merchant in my case Chevron Oil
                  Start a web site featuring gas stations that employ this exploitive practice and drum up public outcry

Getting screwed by a merchant is not ok or right. If you feel you have been treated unfairly use your money and your voice to make your message heard. If a merchant losses enough business they will change or go out of business. Remember they have the products, but you have the money and the power. Use them both wisely. XG

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
           !------------------------!
           ! Cubpack User Interface !
           !------------------------!
Module CUI

USE Precision_Model
USE internal_types

Implicit NONE

PRIVATE

PUBLIC :: CUBATR, CUBPACK_INFO
!-----------------------------------------------------------------------
!***BEGIN PROLOGUE CUBATR
!***DATE WRITTEN   901114   (YYMMDD)
!***REVISION DATE  970620   (YYMMDD)
!***REVISION DATE  980406   (YYMMDD) (MDIV removed)
!***REVISION DATE  000809   (YYMMDD)
!***REVISION DATE  010719   (YYMMDD)
!***REVISION DATE  020715   (YYMMDD) (CUBPACK_INFO added)
!***AUTHOR
!          Ronald Cools, Dept. of Computer Science,
!          Katholieke Universiteit Leuven, Celestijnenlaan 200A,
!          B-3001 Heverlee, Belgium
!          Email:  Ronald.Cools@cs.kuleuven.ac.be
!***REVISION DATE  150729   (YYMMDD) (Made threadsafe / re-entrant)
!***MODIFIER
!          Pat Scott, Dept. of Physics, Imperial College London
!          Email: p.scott@imperial
!
!***PURPOSE  Computation of integrals over a collection of regions.
!
!***DESCRIPTION
!          CUBATR is the driver routine for CUBPACK and the only
!          routine that a user has to deal with (at the moment).
!
!-----------------------------------------------------------------------
PRIVATE :: CUBATR_X, CUBATR_1, CUBATR_CLEAR

!
! Module variables
!
INTEGER, SAVE, PRIVATE :: PreJob=0
INTEGER, PRIVATE                :: BOTTIH,BOTTRH
INTEGER,         DIMENSION(:), PRIVATE, ALLOCATABLE  :: IWork
REAL(kind=stnd), DIMENSION(:), PRIVATE, ALLOCATABLE  :: RWork
TYPE(EPSALG_MEM), PRIVATE              :: M
!$omp threadprivate (PreJob, BOTTIH, BOTTRH, IWork, RWork, M)

INTERFACE CUBATR
   MODULE PROCEDURE CUBATR_X, CUBATR_1, CUBATR_CLEAR
END INTERFACE

CONTAINS

SUBROUTINE CUBATR_X     &
     (DIMENS,NumFun,Integrand,NumRgn,Vertices,RgType,Value,AbsErr,    &
!   and optional parameters
IFAIL,Neval,EpsAbs,EpsRel,Restart,MinPts,MaxPts,Key,Job,Tune)
!-----------------------------------------------------------------------
!   Input parameters
!   ----------------
!
!   DIMENS Integer.
!          The dimension of the region of integration.
!
!   NumFun Integer.
!          Number of components of the integrand.
!
!   Integrand
!          Externally declared function for computing all components
!          of the integrand at the given evaluation point.
!          It must have input parameter X:
!              X(1)   The x-coordinate of the evaluation point.
!              X(2)   The y-coordinate of the evaluation point.
!              ...
!              X(DIMENS) The z-coordinate of the evaluation point.
!         and NumFun, the number of components of the integrand.
!         It must be compatible with the following interface:
!           INTERFACE
!              FUNCTION Integrand(NUMFUN,X) RESULT(Value)
!                USE Precision_Model
!                INTEGER, INTENT(IN) :: NUMFUN
!                REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
!                REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
!              END FUNCTION Integrand
!           END INTERFACE
!
!   NumRgn Integer.
!          The number of given regions.
!
!   Vertices
!          Real array of dimension (DIMENS,DIMENS+1,NumRgn).
!          Vertices(1:DIMENS,K,L) are the x, y, ... coordinates
!          respectively of vertex K of region L, where
!          K = 1,...,DIMENS+1 and L = 1,...,NumRgn.
!
!   RgType Integer array of dimension (NumRgn).
!          RgType(L) describes the type of region L.
!
!   Value  Real array of dimension NumFun.
!          Approximations to all components of the integral if
!          the procedure is restarted.
!
!   AbsErr Real array of dimension NumFun.
!          Estimates of absolute errors if the procedure is restarted.
!
!   IFAIL  Optional integer argument.
!          This follows the NAG convention:
!          IFAIL = 1 : soft silent error
!                      Control returned to calling program.
!          IFAIL = -1: soft noisy error
!                      Error message is printed.
!                      Control returned to calling program.
!          IFAIL = 0 : hard noisy error
!                      Error message is printed and program is stopped.
!          Default IFAIL = -1.
!
!   EpsAbs Optional real argument.
!          Requested absolute error.
!          Default  EpsAbs = 0.
!
!   EpsRel Optional real argument.
!          Requested relative error.
!          Default EpsRel = sqrt(machine precision).
!
!   Restart Optional boolean argument.
!          If Restart = FALSE, this is the first attempt to compute
!                              the integral.
!          If Restart = TRUE, then we restart a previous attempt.
!          In this case the only parameters for CUBATR that may
!          be changed (with respect to the previous call of CUBATR)
!          are MinPts, MaxPts, EpsAbs, EpsRel, Key and Restart.
!          Default Restart = FALSE.
!
!   MinPts Optional integer argument.
!          The minimum allowed number of integrand evaluations.
!          Default MinPts = 0.
!
!   MaxPts Optional integer argument.
!          The maximum allowed number of integrand evaluations.
!          Default MaxPts = enough to do 500 subdivisions.
!
!   Key    Optional integer argument.
!          Can be used by Rule_General to choose between several
!          local integration rules.
!          Default Key = 2 if Dimension=1 and extrapolation is used
!                                        (This corresponds to QAGS)
!          Default Key = 0 otherwise
!
!   Job    Optional integer argument.
!          If |Job| = 0, then nothing will be done except freeing all
!                        allocated memory.
!                        This is usefull after a call of CUBATR if no
!                        Restart will be done later and memory usage
!                        might become an issue later.
!                        Equivalently, one can call CUBATR()
!                        without any arguments.
!                   = 1, the global adaptive algorithm is called
!                   = 2, extrapolation using the epsilon algorithm is used.
!                   = 11, a region will be divided in 2**DIMENS subregions
!                        and the global adaptive algorithm is called.
!                        In combination with Key=0, this resembles DUCTRI and DCUTET.
!                   = 12, a region will be divided in 2 subregions
!                        and the global adaptive algorithm is called.
!                        In combination with Key=3 or 4, this resembles DCUHRE.
!          If Job < 0, then an overview of the Region Collection is dumped.
!          This will create the files tmp_integerstore and tmp_realstore.
!          Default Job = 1.
!
!   Tune   Optional real argument.
!          Can be used by Global_Adapt or the local error estimators
!          to influence the reliability. 0 <= Tune <= 1.
!          Tune = 1 is the most reliable available.
!          Default Tune = 1.
!          Note that this is an experimental and controversial parameter.
!          In this version, only Tune = 1 is supported for all regions.
!
!   Output parameters
!   -----------------
!
!   Value  Real array of dimension NumFun.
!          Approximations to all components of the integral
!
!   AbsErr Real array of dimension NumFun.
!          Estimates of absolute errors.
!
!   NEval  Optional Integer.
!          Number of integrand evaluations used by CUBATR for this call.
!
!   IFAIL  Optional Integer.
!          IFAIL = 0 for normal exit.
!
!            AbsErr(K) <=  EpsAbs or
!            AbsErr(K) <=  ABS(Value(K))*EpsRel with MaxPts or less
!            function evaluations for all values of K,
!            1 <= K <= NumFun .
!
!          IFAIL = 1 if MaxPts was too small to obtain the required
!            accuracy. In this case Global_Adapt returns values of
!            Value with estimated absolute errors AbsErr.
!
!          IFAIL > 1 in more serious case of trouble.
!-----------------------------------------------------------------------
! MODULES USED
USE Check_Input
USE Error_Handling
USE DS_ROUTINES, ONLY: DSCOPY, DSINIT, DSUSED, DSSTAT, DSSUM, DSPINT
USE CubatureRule_General, ONLY: Rule_Cost
USE Global_Adaptive_Algorithm
!***END PROLOGUE CUBATR
!-----------------------------------------------------------------------
!
! Global variables
!
INTERFACE
   FUNCTION Integrand(NUMFUN,X) RESULT(Value)
     USE Precision_Model
     INTEGER, INTENT(IN) :: NUMFUN
     REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
     REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
   END FUNCTION Integrand
END INTERFACE
INTEGER, INTENT(IN)               :: DIMENS,NumFun,NumRgn
INTEGER, DIMENSION(:), INTENT(IN) :: RgType
LOGICAL, INTENT(IN),    OPTIONAL  :: Restart
INTEGER, INTENT(OUT),   OPTIONAL  :: NEval
INTEGER, INTENT(IN),    OPTIONAL  :: Job,Key,MaxPts,MinPts
INTEGER, INTENT(IN OUT), OPTIONAL  :: IFAIL
REAL(kind=stnd), INTENT(IN), OPTIONAL :: Tune,EpsAbs,EpsRel
REAL(kind=stnd), INTENT(IN), DIMENSION(:,:,:) :: Vertices
REAL(kind=stnd), INTENT(IN OUT), DIMENSION(:)  :: AbsErr, Value
!
! Named constants
!
INTEGER, PARAMETER ::  NRINFO=1, NIINFO=5
!
!  Local variables
!
INTEGER :: BLOCK,i,Inform,Leval,LJob,LMaxPts,LMinPts,MinCost, &
           NRVERT,NrSub,NRVACA,MAXRGN,RULCLS,STATUS,Tmp
LOGICAL                                     :: EpsAlg,LRestart
REAL(kind=stnd)                             :: LEpsAbs, LEpsRel
REAL(kind=stnd), DIMENSION(:), ALLOCATABLE  :: TmpRWork
INTEGER,    DIMENSION(:), ALLOCATABLE       :: TmpIWork
TYPE(INTEGRATOR_INFO)                       :: CINFO
TYPE(USER_INFO)                             :: UINFO
!-----------------------------------------------------------------------
!
! Check array sizes
!      Array size mismatch results in hard error.
   Inform = 0
   IF (size(rgtype) < numrgn) THEN
       write(unit=*,fmt=*) "Error: size(rgtype) < numrgn"
       Inform = Inform + 1
   END IF
   IF (size(abserr) < numfun) THEN
       write(unit=*,fmt=*) "Error: size(abserr) < numfun"
       Inform = Inform + 1
   END IF
   IF (size(Value) < numfun) THEN
       write(unit=*,fmt=*) "Error: size(Value) < numfun"
       Inform = Inform + 1
   END IF
   IF ((size(vertices,1) /= dimens) .or. (size(vertices,2) /= dimens+1) &
          .or. (size(vertices,3) < numrgn)) THEN
       Inform = Inform + 1
       write(unit=*,fmt=*)"Error: size(vertices) /= (/dimens,dimens+1,numrgn/)"
   END IF
   IF (Inform /= 0) THEN
      WRITE(unit=*,fmt=*)  "Array size mismatch results in hard error."
      STOP   ! "Array size mismatch results in hard error."
   END IF

!-----------------------------------------------------------------------
IF (PRESENT(NEval)) THEN
    NEval = 0
END IF
!-----------------------------------------------------------------------
IF (PRESENT(Job)) THEN
    LJob = Job
    IF (Job == 0) THEN
       CALL CUBATR_CLEAR()
       RETURN
    END IF
 ELSE
    LJob = 1
END IF
!-----------------------------------------------------------------------
!
! Set optional arguments
!
IF ( PRESENT(Restart)) THEN
                              LRestart = Restart
                       ELSE
                              LRestart = .FALSE.
END IF
IF ( PRESENT(Key)) THEN
                        CINFO%Key = Key
                   ELSE
                       IF ((ABS(LJob) == 2) .AND. (DIMENS == 1)) THEN
                           CINFO%Key = 2   ! simulate QAGS
                       ELSE
                           CINFO%Key = 0
                       END IF
END IF
!-----------------------------------------------------------------------
!
!  Check input parameters
!
IF ( .NOT. LRestart) THEN
         CALL CHECK(DIMENS,NumFun,NumRgn,RgType,Inform, &
                    LJob,IFAIL,EpsAbs,EpsRel,MinPts,MaxPts,Tune)
ELSE IF ( ALLOCATED(IWork) ) THEN
         CALL CHECK(DIMENS,NumFun,BOTTRH,BOTTIH,IWork,  &
                    Inform,LJob,IFAIL,EpsAbs,EpsRel,MinPts,MaxPts,Tune)
ELSE
   Inform = 4096  ! There is nothing to restart from
END IF
IF (Inform /= 0) THEN
    CALL Handle_Error(Inform,IFAIL)
    RETURN
END IF
!-----------------------------------------------------------------------
RULCLS = Rule_Cost( DIMENS, RgType(1), CINFO%Key )
MinCost = RULCLS
DO i = 2,NumRgn
   Tmp = Rule_Cost( DIMENS, RgType(i), CINFO%Key)
   RULCLS = max(RULCLS,Tmp)
   MinCost = MinCost + Tmp
END DO
!-----------------------------------------------------------------------
!
! Set optional arguments
!

IF ( PRESENT(MinPts)) THEN
                              LMinPts = MinPts
                      ELSE
                              LMinPts = 0
END IF
IF ( PRESENT(MaxPts)) THEN
                              LMaxPts = MaxPts
                      ELSE
                              LMaxPts = 500*RULCLS
END IF
IF ( PRESENT(Tune)) THEN
                              CINFO%Tune = Tune
                    ELSE
                              CINFO%Tune = 1
END IF
IF ( PRESENT(EpsAbs)) THEN
                              LEpsAbs = EpsAbs
                      ELSE
                              LEpsAbs = 0
END IF
IF ( PRESENT(EpsRel)) THEN
                              LEpsRel = EpsRel
                      ELSE
                              LEpsRel = SQRT(EPSILON(LEpsRel))
END IF
!-----------------------------------------------------------------------
!
! Set other parameters of the Global Adaptive algorithm
!
!
!  NrSub is an upper limit for the number of subregions after subdivision.
!  This influence memory managment, so don't exagerate here.
IF (DIMENS <= 3) THEN
                              NrSub = 2**DIMENS
                 ELSE
                              NrSub = 4
END IF
!
EpsAlg = ( ABS(LJob) == 2 )
CINFO%UNIFORM_SUBDIV = EpsAlg
CINFO%NrSub = NrSub

IF (( ABS(LJob) == 11) .AND. (DIMENS <= 3)) THEN
   ! simulate dcutri and dcutet ; NrSub = 2**DIMENS ; EpsAlg = .FALSE.
   CINFO%UNIFORM_SUBDIV = .TRUE.
END IF

IF ( ABS(LJob) == 12 ) THEN
   ! simulate dcuhre; NrSub = 2 ; EpsAlg = .FALSE.
   CINFO%NrSub = 2
END IF

NRVERT = DIMENS + 1   ! Only cubes and simplices are implemented here.

UINFO%NumFun = NumFun
UINFO%NumRgn = NumRgn
UINFO%MinPts = LMinPts
UINFO%MaxPts = LMaxPts
UINFO%EpsAbs = LEpsAbs
UINFO%EpsRel = LEpsRel
!-----------------------------------------------------------------------
IF (LRestart) THEN
  ! This requires allocating larger arrays and copying
  ! the region collection.
  ALLOCATE(TmpRWork(SIZE(RWork)),STAT=status)
  IF (status /=  0) THEN
     WRITE(unit=*,fmt=*)  "Problem allocating real workspace."
     STOP  !  "Problem allocating real workspace."
  END IF
  ALLOCATE(TmpIWork(SIZE(IWork)),STAT=status)
  IF (status /=  0) THEN
     WRITE(unit=*,fmt=*) "Problem allocating integer workspace."
     STOP  ! "Problem allocating integer workspace."
  END IF
  MAXRGN = DSUSED(IWork)
  CALL DSCOPY(IWork,RWork,TmpIWork,TmpRWork)
ELSE
  MAXRGN = NumRgn
END IF
!-----------------------------------------------------------------------
! NRVACA is the number of regions the global adaptive algorithm
!        removes from the data structure for further processing.
!        In some routines for shared memory parallel machines
!        this is the variable MDIV
NRVACA = 1
! MAXRGN depends on the number of function evalutions
MAXRGN = MAXRGN + 1 + (NrSub-1)*(LMaxPts - RULCLS*NumRgn)/(RULCLS*NrSub)
!
! Compute length of workspace needed.
!
BOTTIH = MAXRGN*(1+NIINFO) + 15 + NRVACA
BLOCK = NRINFO+NRVERT*DIMENS+2*NumFun
IF (NumFun > 1) THEN
    BLOCK = BLOCK + 1
END IF
BOTTRH = MAXRGN*BLOCK
!
! Allocate space for the region collection
!
IF (ALLOCATED(RWork)) THEN
    DEALLOCATE(RWork)
END IF
ALLOCATE(RWork(BOTTRH),STAT=status)
IF (status /=  0) THEN
   WRITE(unit=*,fmt=*) "Problem allocating real workspace."
   STOP  ! "Problem allocating real workspace."
END IF
IF (ALLOCATED(IWork)) THEN
   DEALLOCATE(IWork)
END IF
ALLOCATE(IWork(BOTTIH),STAT=status)
IF (status  /=  0) THEN
   WRITE(unit=*,fmt=*) "Problem allocating integer workspace."
   STOP  ! "Problem allocating integer workspace."
END IF
!
! Initialise region collection
!
IF ( LRestart) THEN
    CALL DSCOPY(TmpIWork,TmpRWork,IWork,RWork)
    DEALLOCATE(TmpIWork,TmpRWork)
ELSE
   IF ( MinCost > LMaxPts ) THEN
       Inform = 128    ! Dit nummer werd al gebruikt !
   ELSE
       CALL DSINIT(DIMENS,NRVERT,NIINFO,NRINFO,NumFun,NRVACA,   &
                   BOTTIH,BOTTRH,IWork,Inform)
   END IF
   IF (Inform /= 0) THEN
       CALL Handle_Error(Inform,IFAIL)
       RETURN
   END IF
END IF
!-----------------------------------------------------------------------
!
!  Call integration routine
!
If (EpsAlg) THEN
   IF ( LRestart .AND. (PreJob /= ABS(LJob))) THEN
      Inform = 3
   ELSE
      ! Observe that only relevant array sections are passed !
      CALL Global_Adapt_Extrap(DIMENS,CINFO,UINFO,NRVERT,NIINFO,    &
                     NRINFO, Vertices(1:DIMENS,1:NRVERT,1:NUMRGN),  &
                     RgType(1:NUMRGN),Integrand,LRestart,           &
                     Value(1:NUMFUN),AbsErr(1:NUMFUN),LEval,Inform, &
                     RWork,IWork,M)
   END IF
ELSE
   IF ( LRestart .AND. (PreJob /= LJob)) THEN
       IF ( ASSOCIATED(M%RESLA)) THEN
           DEALLOCATE(M%RESLA,M%ERLARG,M%RESULT1,M%ABSERR1,M%RCOPY)
       END IF
       CALL DSPINT(IWork,RWork)
       CALL DSSUM(Value,Abserr,IWork,RWork,Inform)
   END IF
   ! Observe that only relevant array sections are passed !
   CALL Global_Adapt(DIMENS,CINFO,UINFO,NRVERT,NIINFO,NRINFO,    &
                  Vertices(1:DIMENS,1:NRVERT,1:NUMRGN),          &
                  RgType(1:NUMRGN),Integrand,LRestart,           &
                  Value(1:NUMFUN),AbsErr(1:NUMFUN),LEval,Inform, &
                  RWork,IWork)
END IF
IF (PRESENT(NEval)) THEN
   NEval = LEval
END IF
!-----------------------------------------------------------------------
 IF (LJob < 0) THEN
    WRITE(unit=*,fmt=*) "Debug mode: dumping region collection overview."
    CALL DSSTAT(IWork(:),RWork(:))          ! For debugging.
 END IF
!-----------------------------------------------------------------------
! IF ((Inform >= 8) .or. (Inform == 3)) THEN
!    Something went wrong but the data structure remains untouched
!    and so this call can be ignored.
IF ((Inform < 8) .AND. (Inform /= 3)) THEN
   PreJob = ABS(LJob)
END IF
CALL Handle_Error(Inform,IFAIL)
RETURN
END SUBROUTINE CUBATR_X

SUBROUTINE CUBATR_1                                              &
     (DIMENS,Integrand,SVertices,SRgType,SValue,SAbsErr,         &
!   and optional parameters                                      &
      IFAIL,Neval,EpsAbs,EpsRel,Restart,MaxPts,Key,Job)
!-----------------------------------------------------------------------
!   Input parameters
!   ----------------
!
!   DIMENS Integer.
!          The dimension of the region of integration.
!
!   Integrand
!          Externally declared function for computing all components
!          of the integrand at the given evaluation point.
!          It must have input parameter X:
!              X(1)   The x-coordinate of the evaluation point.
!              X(2)   The y-coordinate of the evaluation point.
!              ...
!              X(DIMENS) The z-coordinate of the evaluation point.
!         and NumFun, the number of components of the integrand.
!         It must be compatible with the following interface:
!           INTERFACE
!              FUNCTION Integrand(NUMFUN,X) RESULT(Value)
!                USE Precision_Model
!                INTEGER, INTENT(IN) :: NUMFUN
!                REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
!                REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
!              END FUNCTION Integrand
!           END INTERFACE
!
!   SVertices
!          Real array of dimension (DIMENS,DIMENS+1).
!          Vertices(1:DIMENS,K) are the x, y, ... coordinates
!          respectively of vertex K of the region, where
!          K = 1,...,DIMENS+1.
!
!   SRgType Integer.
!          RgType describes the type of region L.
!
!   SValue Real.
!          Approximation to the integral if the procedure is restarted.
!
!   SAbsErr Real.
!          Estimate of the absolute error if the procedure is restarted.
!
!   IFAIL  Optional integer argument.
!          This follows the NAG convention:
!          IFAIL = 1 : soft silent error
!                      Control returned to calling program.
!          IFAIL = -1: soft noisy error
!                      Error message is printed.
!                      Control returned to calling program.
!          IFAIL = 0 : hard noisy error
!                      Error message is printed and program is stopped.
!          Default IFAIL = -1.
!
!   EpsAbs Optional real argument.
!          Requested absolute error.
!          Default  EpsAbs = 0.
!
!   EpsRel Optional real argument.
!          Requested relative error.
!          Default EpsRel = sqrt(machine precision).
!
!   Restart Optional boolean argument.
!          If Restart = FALSE, this is the first attempt to compute
!                              the integral.
!          If Restart = TRUE, then we restart a previous attempt.
!          In this case the only parameters for CUBATR that may
!          be changed (with respect to the previous call of CUBATR)
!          are MinPts, MaxPts, EpsAbs, EpsRel, Key and Restart.
!          Default Restart = FALSE.
!
!   MaxPts Optional integer argument.
!          The maximum allowed number of integrand evaluations.
!          Default MaxPts = enough to do 500 subdivisions.
!
!   Key    Optional integer argument.
!          Can be used by Rule_General to choose between several
!          local integration rules.
!          Default Key = 2 if Dimension=1 and extrapolation is used
!                                        (This corresponds to QAGS)
!          Default Key = 0 otherwise
!
!   Job    Optional integer argument.
!          If |Job| = 0, then nothing will be done except freeing all
!                        allocated memory.
!                        This is usefull after a call of CUBATR if no
!                        Restart will be done later and memory usage
!                        might become an issue later.
!                        Equivalently, one can call CUBATR()
!                        without any arguments.
!                   = 1, the global adaptive algorithm is called
!                   = 2, extrapolation using the epsilon algorithm is used.
!                   = 11, a region will be divided in 2**DIMENS subregions
!                        and the global adaptive algorithm is called.
!                        In combination with Key=0, this resembles DUCTRI and DCUTET.
!                   = 12, a region will be divided in 2 subregions
!                        and the global adaptive algorithm is called.
!                        In combination with Key=3 or 4, this resembles DCUHRE.
!          If Job < 0, then an overview of the Region Collection is dumped.
!          This will create the files tmp_integerstore and tmp_realstore.
!          Default Job = 1.
!
!   Output parameters
!   -----------------
!
!   SValue Real.
!          Approximation to the integral
!
!   AbsErr Real.
!          Estimate of the absolute error.
!
!   NEval  Optional Integer.
!          Number of integrand evaluations used by CUBATR for this call.
!
!   IFAIL  Optional Integer.
!          IFAIL = 0 for normal exit.
!
!            AbsErr(K) <=  EpsAbs or
!            AbsErr(K) <=  ABS(Value(K))*EpsRel with MaxPts or less
!            function evaluations for all values of K,
!            1 <= K <= NumFun .
!
!          IFAIL = 1 if MaxPts was too small to obtain the required
!            accuracy. In this case Global_Adapt returns values of
!            Value with estimated absolute errors AbsErr.
!
!          IFAIL > 1 in more serious case of trouble.
!-----------------------------------------------------------------------
!
! Global variables
!
INTERFACE
   FUNCTION Integrand(NUMFUN,X) RESULT(Value)
      USE Precision_Model
      INTEGER, INTENT(IN) :: NUMFUN
      REAL(kind=stnd), DIMENSION(:), INTENT(IN) :: X
      REAL(kind=stnd), DIMENSION(NUMFUN) :: Value
   END FUNCTION Integrand
END INTERFACE
LOGICAL, OPTIONAL, INTENT(IN)                 :: Restart
INTEGER, INTENT(IN)                           :: DIMENS,SRgType
INTEGER, INTENT(OUT), OPTIONAL                :: NEval
INTEGER, INTENT(IN), OPTIONAL                 :: Key,MaxPts,Job
INTEGER, INTENT(IN OUT), OPTIONAL             :: IFAIL
REAL(kind=stnd), INTENT(IN), OPTIONAL         :: EpsAbs,EpsRel
REAL(kind=stnd), DIMENSION(:,:), INTENT(IN)   :: SVertices
REAL(kind=stnd), INTENT(IN OUT)               :: SValue,SAbsErr
!
! Local variables
!
INTEGER, DIMENSION(1)                         :: RgType
REAL(kind=stnd), DIMENSION(1)                 :: Value, AbsErr
REAL(kind=stnd), DIMENSION(DIMENS,DIMENS+1,1) :: Vertices
!-------------------
RgType(1) = SRgType
Vertices(:,:,1) = SVertices
IF (PRESENT(Restart)) THEN
    IF ( Restart ) THEN
       Value(1) = SValue
       AbsErr(1) = SAbsErr
    END IF
END IF
CALL CUBATR             &
  (DIMENS,1,Integrand,1,Vertices,RgType,Value,AbsErr,        &
   !   and optional parameters                               &
   ifail=IFAIL,neval=Neval,epsabs=EpsAbs,epsrel=EpsRel,      &
   restart=Restart,maxpts=MaxPts,key=key,job=Job)
SValue = Value(1)
SAbsErr = AbsErr(1)
RETURN
END SUBROUTINE CUBATR_1

SUBROUTINE CUBATR_CLEAR()
  IF ( ALLOCATED(Iwork) ) THEN
     DEALLOCATE(RWork,IWork)
  END IF
  IF ( ASSOCIATED(M%RESLA)) THEN
     DEALLOCATE(M%RESLA,M%ERLARG,M%RESULT1,M%ABSERR1,M%RCOPY)
  END IF
  PreJob = 0
  RETURN
END SUBROUTINE CUBATR_CLEAR

SUBROUTINE CUBPACK_INFO()

REAL(kind=stnd) :: x=1.0e-30  ! lowest accuracy of cubature formula constants

print *," ---------------------------------------------------------------"
print *,"              CUBPACK information"
print *,"              -------------------"
print *," The model for real numbers in the current installed version,"
print *," obtained with the declaration REAL(KIND=stnd), has the"
print *," following characteristics:"
print *,"                base =  ",radix(x)
print *," digits in this base = ", digits(x)
! print *," highest exponent = ", maxexponent(x)
! print *," lowest exponent = "minexponent(x), "(normalized numbers)
print *," This implies:"
print *,"            machine epsilon = ",epsilon(x)
print *,"        largest real number = ", huge(x)
print *," smallest normalized number = ", tiny(x)
print *
print *," The lowest relative error that may be obtained with this"
print "(""  version is about "",G8.2)",max(50*epsilon(x),x)
print *," Asking for lower error will push the routine to use the"
print *," maximal number of function evaluations it is allowed."
print *
print *
print *," This version accepts a collection of hyper-rectangles"
print *,"   (and parallelepipeds) and simplices as integration regions."
print *," Extrapolation using the epsilon-algorithm is available"
print *,"      for dimensions 1, 2 and 3."
print *," The following values of KEY give different integration rules:"
print *," - finite interval: KEY = 1, 2, 3, 4, 5."
print *,"      KEY < 1 defaults to 1; KEY > 5 defaults to 5."
print *," - n-cube:          KEY = 3, 4       uses rule of degree 2*KEY+1"
print *,"           otherwise, uses for a square a rule of degree 13"
print *,"                                 3-cube a rule of degree 11"
print *,"                                        a rule of degree  7"
print *," - n-simplex:       KEY = 1, 2, 3, 4 uses rule of degree 2*KEY+1"
print *,"         otherwise, uses for a triangle a rule of degree 13"
print *,"                            tetrahedron a rule of degree  8"
print *,"                                        a rule of degree  7"
print *," KEY = 0 corresponds to our preferred choice."
print *," ---------------------------------------------------------------"

RETURN
END SUBROUTINE CUBPACK_INFO

END Module CUI

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_name: bgin/MissileSimulation

max_forks_repo_name: bgin/MissileSimulation


text:
      subroutine newvec(v,p,t,nv,neq)
c     implicit double precision (a-h,o-z)                                    dp
      common/main_block/ a(1)
      dimension v(*),p(nv,1),t(*)

      k54=igtpnt(54)
      k16=igtpnt(16)
      do 20 i=1,nv
      call pzero (v,neq)
      do 10 j=1,nv
      l=neq*(j-1)
      call blkcpy (a(k54+l),t,neq)
      do 10 k=1,neq
   10 v(k)=v(k)+t(k)*p(j,i)
      l=neq*(i-1)
   20 call blkcpy (v,a(k16+l),neq)
      return
      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode