text:
! TODO:
! - Study convergence pattern, it's not as reliable as I hoped.
! - Compare loss per step and loss per unit time with old version.
! - Implement custom block matrix multiplication routines.
! - Training logs.
! - Save model and load model (during training).
! - Well-spaced validation data (multiply distances in input and output).
! - Validation-based training (and greedy saving).
! - Embedding support, for integer valued inputs.
! - Embedding support for vector of assigned neighbors for each point.
! - Train only the output, internal, or input layers.
! - Test larger networks on image data.
!
! MAYBE TODO:
! - Use ALLOCATE instead of automatic, because of large N seg faults.
! - Implement similar code in C, compare speeds.

! A piecewise linear regression model.
MODULE PLRM
  USE ISO_FORTRAN_ENV, RT => REAL32
  IMPLICIT NONE
  ! Value used for rectification, default is 0.
  REAL(KIND=RT), PARAMETER :: INITIAL_SHIFT_RANGE = 2.0_RT
  REAL(KIND=RT), PARAMETER :: INITIAL_FLEX = 0.01_RT
  REAL(KIND=RT), PARAMETER :: DISCONTINUITY = 0.0_RT
  REAL(KIND=RT), PARAMETER :: INITIAL_OUTPUT_SCALE = 0.001_RT
  REAL(KIND=RT), PARAMETER :: DEFAULT_INITIAL_STEP = 0.001_RT
  REAL(KIND=RT), PARAMETER :: STEP_COMPOUND_RATE = 1.001_RT
  REAL(KIND=RT), PARAMETER :: INITIAL_STEP_MEAN_CHANGE = 0.1_RT
  REAL(KIND=RT), PARAMETER :: INITIAL_STEP_VAR_CHANGE = 0.001_RT
  INTEGER,       PARAMETER :: MIN_STEPS_TO_STABILITY = 0
  INTEGER,       PARAMETER :: DEFAULT_STEPS = 100

  ! Value for bias base (in terms of standard deviation for unit
  !  variance data), default value is 2.

  ! MDI = Model dimension -- input
  ! MDS = Model dimension -- states (internal)
  ! MNS = Model number -- states (internal)
  ! MDO = Model dimension -- output
  INTEGER :: MDI, MDS, MNS, MDO
  ! Model parameters, for evaluation.
  REAL(KIND=RT), DIMENSION(:,:),   ALLOCATABLE :: INPUT_VECS
  REAL(KIND=RT), DIMENSION(:),     ALLOCATABLE :: INPUT_SHIFT
  REAL(KIND=RT), DIMENSION(:),     ALLOCATABLE :: INPUT_FLEX
  REAL(KIND=RT), DIMENSION(:,:,:), ALLOCATABLE :: INTERNAL_VECS
  REAL(KIND=RT), DIMENSION(:,:),   ALLOCATABLE :: INTERNAL_SHIFT
  REAL(KIND=RT), DIMENSION(:,:),   ALLOCATABLE :: INTERNAL_FLEX
  REAL(KIND=RT), DIMENSION(:,:),   ALLOCATABLE :: OUTPUT_VECS
  REAL(KIND=RT), DIMENSION(:),     ALLOCATABLE :: OUTPUT_SHIFT

  PUBLIC :: NEW_MODEL, INIT_MODEL, MINIMIZE_MSE, EVALUATE_ONE
  PRIVATE :: SSE_GRADIENT

  ! Externally provided matrix-matrix multiplication routine.
  EXTERNAL :: SGEMM
  ! EXTERNAL :: DGEMM

CONTAINS

  ! Allocate storage for a new model given paramters:
  !   DI -- Dimension of input vector.
  !   DS -- Dimension of internal vectors spaces.
  !   NS -- Number of sequential internal vector spaces.
  !   DO -- Dimension of output vector.
  SUBROUTINE NEW_MODEL(DI, DS, NS, DO)
    INTEGER, INTENT(IN) :: DI, DS, NS, DO
    ! Store the integer sizes related to this model.
    MDI = DI
    MDS = DS
    MNS = NS
    MDO = DO
    ! Deallocate any of the internal arrays if they have already been allocated.
    IF (ALLOCATED(INPUT_VECS))     DEALLOCATE(INPUT_VECS)
    IF (ALLOCATED(INPUT_SHIFT))    DEALLOCATE(INPUT_SHIFT)
    IF (ALLOCATED(INPUT_FLEX))     DEALLOCATE(INPUT_FLEX)
    IF (ALLOCATED(INTERNAL_VECS))  DEALLOCATE(INTERNAL_VECS)
    IF (ALLOCATED(INTERNAL_SHIFT)) DEALLOCATE(INTERNAL_SHIFT)
    IF (ALLOCATED(INTERNAL_FLEX))  DEALLOCATE(INTERNAL_FLEX)
    IF (ALLOCATED(OUTPUT_VECS))    DEALLOCATE(OUTPUT_VECS)
    IF (ALLOCATED(OUTPUT_SHIFT))   DEALLOCATE(OUTPUT_SHIFT)
    ! Allocate space for all model parameters.
    ALLOCATE( &
         INPUT_VECS(DI,DS), &
         INPUT_SHIFT(DS), &
         INPUT_FLEX(DS), &
         INTERNAL_VECS(DS, DS, NS-1), &
         INTERNAL_SHIFT(DS, NS-1), &
         INTERNAL_FLEX(DS, NS-1), &
         OUTPUT_VECS(DS, DO), &
         OUTPUT_SHIFT(DO) &
         )
  END SUBROUTINE NEW_MODEL

  ! Initialize the weights for a model, optionally provide a random seed.
  SUBROUTINE INIT_MODEL(SEED)
    INTEGER, INTENT(IN), OPTIONAL :: SEED
    !  Storage for seeding the random number generator (for repeatability).
    INTEGER, DIMENSION(:), ALLOCATABLE :: SEED_ARRAY
    ! Local iterator.
    INTEGER :: I
    ! Set a random seed, if one was provided (otherwise leave default).
    IF (PRESENT(SEED)) THEN
       CALL RANDOM_SEED(SIZE=I)
       ALLOCATE(SEED_ARRAY(I))
       SEED_ARRAY(:) = SEED
       CALL RANDOM_SEED(PUT=SEED_ARRAY(:))
    END IF
    ! Generate well spaced random unit-length vectors (no scaling biases)
    !  for all initial parameters in the input, internal, and output.
    CALL RANDOM_UNIT_VECTORS(INPUT_VECS(:,:))
    DO I = 1, MNS-1
       CALL RANDOM_UNIT_VECTORS(INTERNAL_VECS(:,:,I))
    END DO
    CALL RANDOM_UNIT_VECTORS(OUTPUT_VECS(:,:))
    ! Make the output vectors have very small magnitude initially.
    OUTPUT_VECS(:,:) = OUTPUT_VECS(:,:) * INITIAL_OUTPUT_SCALE
    ! Generate random biases for inputs and internal layers, zero
    !  bias for the output layer (first two will be rescaled).
    CALL RANDOM_NUMBER(INPUT_SHIFT(:))
    CALL RANDOM_NUMBER(INTERNAL_SHIFT(:,:))
    OUTPUT_SHIFT(:) = 0.0_RT
    ! Values are in range [-2,2], assuming no scaling changes,
    !  this will capture 2 standard deviations in either direction.
    INPUT_SHIFT(:) = INITIAL_SHIFT_RANGE * (INPUT_SHIFT(:) * 2.0_RT - 1.0_RT)
    INTERNAL_SHIFT(:,:) = INITIAL_SHIFT_RANGE * (INTERNAL_SHIFT(:,:) * 2.0_RT - 1.0_RT)
    ! Set initial *_FLEX to be equal to default small value.
    INPUT_FLEX(:) = INITIAL_FLEX
    INTERNAL_FLEX(:,:) = INITIAL_FLEX

  CONTAINS

    ! Generate randomly distributed vectors on the N-sphere.
    SUBROUTINE RANDOM_UNIT_VECTORS(COLUMN_VECTORS)
      REAL(KIND=RT), INTENT(OUT), DIMENSION(:,:) :: COLUMN_VECTORS
      INTEGER :: I
      ! Generate random numbers in the range [0,1].
      CALL RANDOM_NUMBER(COLUMN_VECTORS(:,:))
      ! Map those random numbers to the range [-1,1].
      COLUMN_VECTORS(:,:) = 2.0_RT * COLUMN_VECTORS(:,:) - 1.0_RT
      ! Make the vectors uniformly distributed on the unit ball (for dimension > 1).
      IF (SIZE(COLUMN_VECTORS,1) .GT. 1) THEN
         ! Compute the inverse hyperbolic tangent of all values.
         COLUMN_VECTORS(:,:) = ATANH(COLUMN_VECTORS(:,:))
         ! Normalize all vectors to have unit length.
         DO I = 1, SIZE(COLUMN_VECTORS,2)
            COLUMN_VECTORS(:,I) = COLUMN_VECTORS(:,I) / NORM2(COLUMN_VECTORS(:,I))
         END DO
      END IF
      ! Orthogonalize the first components of the column
      !  vectors to ensure those are well spaced.
      I = MIN(SIZE(COLUMN_VECTORS,1), SIZE(COLUMN_VECTORS,2))
      CALL ORTHOGONALIZE(COLUMN_VECTORS(:,1:I))
    END SUBROUTINE RANDOM_UNIT_VECTORS

    ! Orthogonalize A with Householder reflections.
    SUBROUTINE ORTHOGONALIZE(A)
      REAL(KIND=RT), INTENT(INOUT), DIMENSION(:,:) :: A
      ! Local variables (perform orthogonalization by computing
      !  Householder reflectors that make A upper diagonal).
      REAL(KIND=RT), DIMENSION(SIZE(A,1), SIZE(A,2)) :: Q
      REAL(KIND=RT), DIMENSION(SIZE(A,1), SIZE(A,1)) :: QI
      INTEGER :: I, J, M, N, L
      REAL(KIND=RT) :: R
      M = SIZE(A,1)
      N = SIZE(A,2)
      L = MIN(M, N)
      ! Set Q to have ones on the diagonal.
      Q(:,:) = 0.0_RT
      FORALL (I=1:L) Q(I,I) = 1.0_RT
      ! Set Qi to the identity.
      QI(:,:) = 0.0_RT
      FORALL (I=1:M) QI(I,I) = 1.0_RT
      ! Compute Householder reflectors.
      DO I = 1, L-1
         ! Copy the column vector from A and compute its length.
         QI(I:M,I) = A(I:M,I)
         QI(I,I) = NORM2(QI(I:M,I))
         ! Make the norm have the opposite sign of the first component.
         IF (A(I,I) .LT. 0.0_RT) QI(I,I) = -QI(I,I)
         ! Subtract this computed value from the first component and
         !  normalize the vector to have unit length.
         QI(I,I) = A(I,I) - QI(I,I)
         QI(I:M,I) = QI(I:M,I) / NORM2(QI(I:M,I))
         ! Perform the outer product between the new vector and itself.
         FORALL (J=I+1:M)
            QI(I:M,J) = QI(J,I) * QI(I:M,I)
         END FORALL
         QI(I:M,I) = QI(I,I) * QI(I:M,I)
         ! Subtract twice this value all from the identity matrix.
         QI(I:M,I:M) = -2.0_RT * QI(I:M,I:M)
         FORALL (J=I:M) QI(J,J) = QI(J,J) + 1.0_RT
         ! Update the matrix as well as Q.
         !  TODO: Use SGEMM instead of builtin MATMUL.
         A(:,:) = MATMUL(QI(:,:), A(:,:))
         Q(:,:) = MATMUL(QI(:,:), Q(:,:))
         ! Reset the current row and column of QI to be from the identity.
         QI(I,I) = 1.0_RT
         QI(I+1:M,I) = 0.0_RT
         QI(I,I+1:M) = 0.0_RT
      END DO
      ! Overwrite A with the orthognal Q.
      A(1:M,1:L) = Q(1:M,1:L)
      ! Randomly flip the sign of some of the vectors in A.
      DO I = 1, L
         CALL RANDOM_NUMBER(R)
         IF (R .LT. 0.5_RT) THEN
            A(:,I) = -A(:,I)
         END IF
      END DO
    END SUBROUTINE ORTHOGONALIZE

  END SUBROUTINE INIT_MODEL

  ! Evaluate the piecewise linear regression model.
  SUBROUTINE EVALUATE_ONE(INPUT, OUTPUT)
    REAL(KIND=RT), INTENT(IN), DIMENSION(:) :: INPUT
    ! Record of values after the last transformation.
    REAL(KIND=RT), INTENT(OUT), DIMENSION(:) :: OUTPUT
    ! Internal values.
    REAL(KIND=RT), DIMENSION(MDS) :: INTERNAL_VALUES
    ! Compute the number of nodes (to know where "bias" value is).
    INTEGER :: I, IP1
    ! Make sure that this model has been initialized.
    IF (.NOT. ALLOCATED(INPUT_VECS)) THEN
       OUTPUT(:) = 0.0
       RETURN
    END IF
    ! Compute the input layer.
    INTERNAL_VALUES(:) = INPUT_SHIFT(:) + MATMUL(INPUT(:), INPUT_VECS(:,:))
    WHERE (INTERNAL_VALUES(:) .LT. DISCONTINUITY)
       INTERNAL_VALUES(:) = INTERNAL_VALUES(:) * INPUT_FLEX(:)
    END WHERE
    ! Compute the next set of internal values with a rectified activation.
    DO I = 1, MNS-1
       IP1 = I+1
       INTERNAL_VALUES(:) = INTERNAL_SHIFT(:,I) + &
            MATMUL(INTERNAL_VALUES(:), INTERNAL_VECS(:,:,I))
       WHERE (INTERNAL_VALUES(:) .LT. DISCONTINUITY)
          INTERNAL_VALUES(:) = INTERNAL_VALUES(:) * INTERNAL_FLEX(:,I)
       END WHERE
    END DO
    ! Compute the output.
    OUTPUT(:) = OUTPUT_SHIFT(:) + &
         MATMUL(INTERNAL_VALUES(:), OUTPUT_VECS(:,:))
  END SUBROUTINE EVALUATE_ONE

  ! Evaluate the piecewise linear regression model.
  SUBROUTINE EVALUATE(INPUTS, OUTPUTS)
    REAL(KIND=RT), INTENT(IN), DIMENSION(:,:) :: INPUTS
    REAL(KIND=RT), INTENT(OUT), DIMENSION(:,:) :: OUTPUTS
    ! Internal values.
    REAL(KIND=RT), DIMENSION(:,:), ALLOCATABLE :: TEMP_VALUES, INTERNAL_VALUES
    ! Compute the number of nodes (to know where "bias" value is).
    INTEGER :: I, D, L, N
    ! Make sure that this model has been initialized.
    IF (.NOT. ALLOCATED(INPUT_VECS)) THEN
       OUTPUTS(:,:) = 0.0
       RETURN
    END IF
    ! Get the number of points.
    N = SIZE(INPUTS,2)
    ! Allocate storage for the internal values.
    ALLOCATE(INTERNAL_VALUES(1:N,1:MDS), TEMP_VALUES(1:N,1:MDS))
    ! Compute the input transformation.
    CALL SGEMM('T', 'N', N, MDS, MDI, 1.0_RT, &
         INPUTS(:,:), SIZE(INPUTS,1), &
         INPUT_VECS(:,:), SIZE(INPUT_VECS,1), &
         0.0_RT, INTERNAL_VALUES(:,:), SIZE(INTERNAL_VALUES,1))
    DO D = 1, MDS
       INTERNAL_VALUES(:,D) = INPUT_SHIFT(D) + INTERNAL_VALUES(:,D)
       WHERE (INTERNAL_VALUES(:,D) .LT. DISCONTINUITY)
          INTERNAL_VALUES(:,D) = INTERNAL_VALUES(:,D) * INPUT_FLEX(D)
       END WHERE
    END DO
    ! Compute the next set of internal values with a rectified activation.
    DO L = 1, MNS-1
       CALL SGEMM('N', 'N', N, MDS, MDS, 1.0_RT, &
            INTERNAL_VALUES(:,:), SIZE(INTERNAL_VALUES,1), &
            INTERNAL_VECS(:,:,L), SIZE(INTERNAL_VECS,1), &
            0.0_RT, TEMP_VALUES(:,:), SIZE(TEMP_VALUES,1))
       INTERNAL_VALUES(:,:) = TEMP_VALUES(:,:)
       DO D = 1, MDS
          INTERNAL_VALUES(:,D) =  INTERNAL_SHIFT(D,L) + INTERNAL_VALUES(:,D)
          WHERE (INTERNAL_VALUES(:,D) .LT. DISCONTINUITY)
             INTERNAL_VALUES(:,D) = INTERNAL_VALUES(:,D) * INTERNAL_FLEX(D,L)
          END WHERE
       END DO
    END DO
    ! Return the final embedded layer if the outputs have the size of the embedding.
    IF ((MDS .NE. MDO) .AND. (SIZE(OUTPUTS,1) .EQ. MDS)) THEN
       ! Do the necessary transpose operation one dimension at a time, produce embedding.
       DO L = 1, MDS
          OUTPUTS(L,1:N) = INTERNAL_VALUES(1:N,L)
       END DO
    ! Otherwise, assume regular output computation.
    ELSE
       CALL SGEMM('T', 'T', MDO, N, MDS, 1.0_RT, &
            OUTPUT_VECS(:,:), SIZE(OUTPUT_VECS,1), &
            INTERNAL_VALUES(:,:), SIZE(INTERNAL_VALUES,1), &
            0.0_RT, OUTPUTS(:,:), SIZE(OUTPUTS,1))
       DO D = 1, MDO
          OUTPUTS(:,D) = OUTPUT_SHIFT(D) + OUTPUTS(:,D)
       END DO
    END IF
    ! Deallocate temporary variables.
    DEALLOCATE(INTERNAL_VALUES, TEMP_VALUES)
  END SUBROUTINE EVALUATE

  ! Compute the gradient of the sum of squared error of this regression
  ! model with respect to its parameters given input and output pairs.
  SUBROUTINE SSE_GRADIENT(INPUTS, OUTPUTS, SUM_SQUARED_ERROR, &
       INPUT_VECS_GRADIENT, INPUT_SHIFT_GRADIENT, INPUT_FLEX_GRADIENT,&
       INTERNAL_VECS_GRADIENT, INTERNAL_SHIFT_GRADIENT, INTERNAL_FLEX_GRADIENT,&
       OUTPUT_VECS_GRADIENT, OUTPUT_SHIFT_GRADIENT)
    ! Data values stored with contiguous points: shape = (MDI,N)
    REAL(KIND=RT), INTENT(IN), DIMENSION(:,:) :: INPUTS
    REAL(KIND=RT), INTENT(IN), DIMENSION(:,:) :: OUTPUTS
    ! Sum (over all data) squared error (summed over dimensions).
    REAL(KIND=RT), INTENT(INOUT) :: SUM_SQUARED_ERROR
    ! Model parameter gradients.
    REAL(KIND=RT), INTENT(INOUT), DIMENSION(:,:)   :: INPUT_VECS_GRADIENT
    REAL(KIND=RT), INTENT(INOUT), DIMENSION(:)     :: INPUT_SHIFT_GRADIENT
    REAL(KIND=RT), INTENT(INOUT), DIMENSION(:)     :: INPUT_FLEX_GRADIENT
    REAL(KIND=RT), INTENT(INOUT), DIMENSION(:,:,:) :: INTERNAL_VECS_GRADIENT
    REAL(KIND=RT), INTENT(INOUT), DIMENSION(:,:)   :: INTERNAL_SHIFT_GRADIENT
    REAL(KIND=RT), INTENT(INOUT), DIMENSION(:,:)   :: INTERNAL_FLEX_GRADIENT
    REAL(KIND=RT), INTENT(INOUT), DIMENSION(:,:)   :: OUTPUT_VECS_GRADIENT
    REAL(KIND=RT), INTENT(INOUT), DIMENSION(:)     :: OUTPUT_SHIFT_GRADIENT
    ! Internal values and lists of which ones are nonzero.
    REAL(KIND=RT), DIMENSION(SIZE(INPUTS,2),MDO)     :: OUTPUT_GRADIENT
    REAL(KIND=RT), DIMENSION(SIZE(INPUTS,2),MDS,MNS) :: INTERNAL_VALUES
    LOGICAL,       DIMENSION(SIZE(INPUTS,2),MDS,MNS) :: LT_DISCONTINUITY
    REAL(KIND=RT), DIMENSION(SIZE(INPUTS,2),MDS)     :: INTERNAL_TEMP
    ! Number of points.
    REAL(KIND=RT) :: N
    INTEGER :: D
    ! Evaluate the model at all data points, store outputs in "OUTPUT_GRADIENT"
    CALL EVALUATE_BATCH()
    ! Compute the new gradient.
    DO D = 1, MDO
       OUTPUT_GRADIENT(:,D) = OUTPUT_GRADIENT(:,D) - OUTPUTS(D,:)
    END DO
    ! Get the number of points.
    N = REAL(SIZE(INPUTS,2), RT)
    ! Compute the current sum squared error.
    SUM_SQUARED_ERROR = SUM_SQUARED_ERROR + SUM(OUTPUT_GRADIENT(:,:)**2)
    ! Compute the gradient of parameters with respect to the output gradient.
    CALL GRADIENT_BATCH()

  CONTAINS

    ! Evaluate the piecewise linear regression model.
    SUBROUTINE EVALUATE_BATCH()
      ! D   - dimension index
      ! L   - layer index
      ! LP1 - layer index "plus 1" -> "P1"
      INTEGER :: D, L, LP1
      ! Compute the input transformation.
      CALL SGEMM('T', 'N', SIZE(INPUTS,2), MDS, MDI, 1.0_RT, &
           INPUTS(:,:), SIZE(INPUTS,1), &
           INPUT_VECS(:,:), SIZE(INPUT_VECS,1), &
           0.0_RT, INTERNAL_VALUES(:,:,1), SIZE(INTERNAL_VALUES,1))
      DO D = 1, MDS
         INTERNAL_VALUES(:,D,1) = INTERNAL_VALUES(:,D,1) + INPUT_SHIFT(D)
         LT_DISCONTINUITY(:,D,1) = INTERNAL_VALUES(:,D,1) .LT. DISCONTINUITY
         WHERE (LT_DISCONTINUITY(:,D,1))
            INTERNAL_VALUES(:,D,1) = INTERNAL_VALUES(:,D,1) * INPUT_FLEX(D)
         END WHERE
      END DO
      ! Compute the next set of internal values with a rectified activation.
      DO L = 1, MNS-1
         LP1 = L+1
         CALL SGEMM('N', 'N', SIZE(INPUTS,2), MDS, MDS, 1.0_RT, &
              INTERNAL_VALUES(:,:,L), SIZE(INTERNAL_VALUES,1), &
              INTERNAL_VECS(:,:,L), SIZE(INTERNAL_VECS,1), &
              0.0_RT, INTERNAL_VALUES(:,:,LP1), SIZE(INTERNAL_VALUES,1))
         DO D = 1, MDS
            INTERNAL_VALUES(:,D,LP1) = INTERNAL_VALUES(:,D,LP1) + INTERNAL_SHIFT(D,L)
            LT_DISCONTINUITY(:,D,LP1) = INTERNAL_VALUES(:,D,LP1) .LT. DISCONTINUITY
            WHERE (LT_DISCONTINUITY(:,D,LP1))
               INTERNAL_VALUES(:,D,LP1) = INTERNAL_VALUES(:,D,LP1) * INTERNAL_FLEX(D,L)
            END WHERE
         END DO
      END DO
      ! Compute the output.
      CALL SGEMM('N', 'N', SIZE(INPUTS,2), MDO, MDS, 1.0_RT, &
           INTERNAL_VALUES(:,:,L), SIZE(INTERNAL_VALUES,1), &
           OUTPUT_VECS(:,:), SIZE(OUTPUT_VECS,1), &
           0.0_RT, OUTPUT_GRADIENT(:,:), SIZE(OUTPUT_GRADIENT,1))
      DO D = 1, MDO
         OUTPUT_GRADIENT(:,D) = OUTPUT_SHIFT(D) + OUTPUT_GRADIENT(:,D)
      END DO
    END SUBROUTINE EVALUATE_BATCH

    SUBROUTINE GRADIENT_BATCH()
      ! D   - dimension index
      ! L   - layer index
      ! LP1 - layer index "plus 1" -> "P1"
      INTEGER :: D, L, LP1
      ! Compute the average output gradient.
      OUTPUT_SHIFT_GRADIENT(:) = SUM(OUTPUT_GRADIENT(:,:), 1) / N &
           + OUTPUT_SHIFT_GRADIENT(:)
      ! Compute output parameter gradients.
      CALL SGEMM('T', 'N', MDS, MDO, SIZE(INPUTS,2), 1.0_RT / N, &
           INTERNAL_VALUES(:,:,MNS), SIZE(INTERNAL_VALUES,1), &
           OUTPUT_GRADIENT(:,:), SIZE(OUTPUT_GRADIENT,1), &
           1.0_RT, OUTPUT_VECS_GRADIENT(:,:), SIZE(OUTPUT_VECS_GRADIENT,1))
      ! Propogate the gradient back to the last internal vector space.
      CALL SGEMM('N', 'T', SIZE(INPUTS,2), MDS, MDO, 1.0_RT, &
           OUTPUT_GRADIENT(:,:), SIZE(OUTPUT_GRADIENT,1), &
           OUTPUT_VECS(:,:), SIZE(OUTPUT_VECS,1), &
           0.0_RT, INTERNAL_TEMP(:,:), SIZE(INTERNAL_TEMP,1))
      ! Cycle over all internal layers.
      INTERNAL_REPRESENTATIONS : DO L = MNS-1, 1, -1
         LP1 = L+1
         DO D = 1, MDS
            ! Compute the flexure gradient.
            IF (INTERNAL_FLEX(D,L) .NE. 0.0_RT) THEN
               INTERNAL_FLEX_GRADIENT(D,L) = SUM(&
                    INTERNAL_TEMP(:,D) * INTERNAL_VALUES(:,D,LP1), 1, &
                    LT_DISCONTINUITY(:,D,LP1)) / INTERNAL_FLEX(D,L)
            ELSE
               INTERNAL_FLEX_GRADIENT(D,L) = 0.0_RT
            END IF
            ! Propogate the error gradient back to the preceding vectors.
            WHERE (LT_DISCONTINUITY(:,D,LP1))
               INTERNAL_VALUES(:,D,LP1) = INTERNAL_TEMP(:,D) * INTERNAL_FLEX(D,L)
            ELSEWHERE
               INTERNAL_VALUES(:,D,LP1) = INTERNAL_TEMP(:,D)
            END WHERE
         END DO
         ! Compute the shift gradient.
         INTERNAL_SHIFT_GRADIENT(:,L) = SUM(INTERNAL_VALUES(:,:,LP1), 1) / N &
              + INTERNAL_SHIFT_GRADIENT(:,L)
         ! Compute the gradient with respect to each output and all inputs.
         CALL SGEMM('T', 'N', MDS, MDS, SIZE(INPUTS,2), 1.0_RT / N, &
              INTERNAL_VALUES(:,:,L), SIZE(INTERNAL_VALUES,1), &
              INTERNAL_VALUES(:,:,LP1), SIZE(INTERNAL_VALUES,1), &
              1.0_RT, INTERNAL_VECS_GRADIENT(:,:,L), SIZE(INTERNAL_VECS_GRADIENT,1))
         ! Propogate the gradient to the immediately preceding layer's flexure.
         CALL SGEMM('N', 'T', SIZE(INPUTS,2), MDS, MDS, 1.0_RT, &
              INTERNAL_VALUES(:,:,LP1), SIZE(INTERNAL_VALUES,1), &
              INTERNAL_VECS(:,:,L), SIZE(INTERNAL_VECS,1), &
              0.0_RT, INTERNAL_TEMP(:,:), SIZE(INTERNAL_TEMP,1))
      END DO INTERNAL_REPRESENTATIONS
      ! Compute the gradients going into the first layer.
      DO D = 1, MDS
         ! Compute the flexure gradient.
         IF (INPUT_FLEX(D) .NE. 0.0_RT) THEN
            INPUT_FLEX_GRADIENT(D) = SUM(&
                 INTERNAL_TEMP(:,D) * INTERNAL_VALUES(:,D,1), 1, &
                 LT_DISCONTINUITY(:,D,1)) / INPUT_FLEX(D)
         ELSE
            INPUT_FLEX_GRADIENT(D) = 0.0_RT
         END IF
         ! Propogate the error gradient back to the preceding vectors.
         WHERE (LT_DISCONTINUITY(:,D,1))
            INTERNAL_VALUES(:,D,1) = INTERNAL_TEMP(:,D) * INPUT_FLEX(D)
         ELSEWHERE
            INTERNAL_VALUES(:,D,1) = INTERNAL_TEMP(:,D)
         END WHERE
      END DO
      ! Compute the input parameter gradients.
      INPUT_SHIFT_GRADIENT(:) = SUM(INTERNAL_VALUES(:,:,1), 1) / N &
           + INPUT_SHIFT_GRADIENT(:)
      ! Compute the gradient of all input parameters.
      !   [the INPUTS are transposed already, shape = (MDI,N)]
      CALL SGEMM('N', 'N', MDI, MDS, SIZE(INPUTS,2), 1.0_RT / N, &
           INPUTS(:,:), SIZE(INPUTS,1), &
           INTERNAL_VALUES(:,:,1), SIZE(INTERNAL_VALUES,1), &
           1.0_RT, INPUT_VECS_GRADIENT(:,:), SIZE(INPUT_VECS_GRADIENT,1))
    END SUBROUTINE GRADIENT_BATCH

  END SUBROUTINE SSE_GRADIENT

  ! Fit input / output pairs by minimizing mean squared error.
  SUBROUTINE MINIMIZE_MSE(X, Y, STEPS, BATCH_SIZE, NUM_THREADS, &
       STEP_SIZE, KEEP_BEST, MEAN_SQUARED_ERROR, RECORD)
    REAL(KIND=RT), INTENT(IN), DIMENSION(:,:) :: X
    REAL(KIND=RT), INTENT(IN), DIMENSION(:,:) :: Y
    INTEGER,       INTENT(IN) :: STEPS
    INTEGER,       INTENT(IN), OPTIONAL :: BATCH_SIZE, NUM_THREADS
    REAL(KIND=RT), INTENT(IN), OPTIONAL :: STEP_SIZE
    LOGICAL,       INTENT(IN), OPTIONAL :: KEEP_BEST
    REAL(KIND=RT), INTENT(OUT) :: MEAN_SQUARED_ERROR
    REAL(KIND=RT), INTENT(OUT), DIMENSION(STEPS), OPTIONAL :: RECORD
    !  gradient step arrays, 4 copies of internal model (total of 5).
    REAL(KIND=RT), DIMENSION(MDI,MDS)       :: V1_GRAD, V1_MEAN, V1_VAR, BEST_V1
    REAL(KIND=RT), DIMENSION(MDS)           :: S1_GRAD, S1_MEAN, S1_VAR, BEST_S1
    REAL(KIND=RT), DIMENSION(MDS)           :: F1_GRAD, F1_MEAN, F1_VAR, BEST_F1
    REAL(KIND=RT), DIMENSION(MDS,MDS,MNS-1) :: V2_GRAD, V2_MEAN, V2_VAR, BEST_V2
    REAL(KIND=RT), DIMENSION(MDS,MNS-1)     :: S2_GRAD, S2_MEAN, S2_VAR, BEST_S2
    REAL(KIND=RT), DIMENSION(MDS,MNS-1)     :: F2_GRAD, F2_MEAN, F2_VAR, BEST_F2
    REAL(KIND=RT), DIMENSION(MDS,MDO)       :: V3_GRAD, V3_MEAN, V3_VAR, BEST_V3
    REAL(KIND=RT), DIMENSION(MDO)           :: S3_GRAD, S3_MEAN, S3_VAR, BEST_S3
    !  local integers
    LOGICAL :: REVERT_TO_BEST
    INTEGER :: BS, I, J, L, NX, NT, BATCH_START, BATCH_END
    REAL(KIND=RT) :: RNX, BATCHES, PREV_MSE, BEST_MSE, SCALAR
    REAL(KIND=RT) :: STEP_FACTOR, STEP_MEAN_CHANGE, STEP_MEAN_REMAIN, &
         STEP_VAR_CHANGE, STEP_VAR_REMAIN
    ! Function that is defined by OpenMP.
    INTERFACE
       FUNCTION OMP_GET_MAX_THREADS()
         INTEGER :: OMP_GET_MAX_THREADS
       END FUNCTION OMP_GET_MAX_THREADS
    END INTERFACE
    ! Number of points.
    NX = SIZE(X,2)
    RNX = REAL(NX, RT)
    ! Set the number of threads.
    IF (PRESENT(NUM_THREADS)) THEN
       NT = NUM_THREADS
    ELSE
       NT = OMP_GET_MAX_THREADS()
    END IF
    ! Set the batch size.
    IF (PRESENT(BATCH_SIZE)) THEN
       BS = MIN(BATCH_SIZE, NX)
    ELSE
       BS = (NX + NT - 1) / NT
    END IF
    ! Set the step factor.
    IF (PRESENT(STEP_SIZE)) THEN
       STEP_FACTOR = STEP_SIZE
    ELSE
       STEP_FACTOR = DEFAULT_INITIAL_STEP
    END IF
    ! Set the "keep best" boolean.
    IF (PRESENT(KEEP_BEST)) THEN
       REVERT_TO_BEST = KEEP_BEST
    ELSE
       REVERT_TO_BEST = (STEPS .GE. MIN_STEPS_TO_STABILITY)
    END IF
    ! Initial rates of change of mean and variance values.
    STEP_MEAN_CHANGE = INITIAL_STEP_MEAN_CHANGE
    STEP_MEAN_REMAIN = 1.0_RT - STEP_MEAN_CHANGE
    STEP_VAR_CHANGE = INITIAL_STEP_VAR_CHANGE
    STEP_VAR_REMAIN = 1.0_RT - STEP_VAR_CHANGE
    ! Initial mean squared error is "max float value".
    PREV_MSE = HUGE(PREV_MSE)
    BEST_MSE = HUGE(BEST_MSE)
    ! Set the average step sizes.
    V1_MEAN(:,:)   = 0.0_RT
    S1_MEAN(:)     = 0.0_RT
    F1_MEAN(:)     = 0.0_RT
    V2_MEAN(:,:,:) = 0.0_RT
    S2_MEAN(:,:)   = 0.0_RT
    F2_MEAN(:,:)   = 0.0_RT
    V3_MEAN(:,:)   = 0.0_RT
    S3_MEAN(:)     = 0.0_RT
    ! Set the average step size variances.
    V1_VAR(:,:)   = 0.0_RT
    S1_VAR(:)     = 0.0_RT
    F1_VAR(:)     = 0.0_RT
    V2_VAR(:,:,:) = 0.0_RT
    S2_VAR(:,:)   = 0.0_RT
    F2_VAR(:,:)   = 0.0_RT
    V3_VAR(:,:)   = 0.0_RT
    S3_VAR(:)     = 0.0_RT
    ! Iterate, taking steps with the average gradient over all data.
    fit_loop : DO I = 1, STEPS
       ! Compute the average gradient over all points.
       MEAN_SQUARED_ERROR = 0.0_RT
       ! Set gradients to zero initially.
       V1_GRAD(:,:)   = 0.0_RT
       S1_GRAD(:)     = 0.0_RT
       F1_GRAD(:)     = 0.0_RT
       V2_GRAD(:,:,:) = 0.0_RT
       S2_GRAD(:,:)   = 0.0_RT
       F2_GRAD(:,:)   = 0.0_RT
       V3_GRAD(:,:)   = 0.0_RT
       S3_GRAD(:)     = 0.0_RT
       ! Count the number of batches.
       BATCHES = 0.0_RT
       !$OMP PARALLEL DO NUM_THREADS(NT) PRIVATE(BATCH_END) &
       !$OMP&  REDUCTION(+: BATCHES, MEAN_SQUARED_ERROR, &
       !$OMP&  V1_GRAD, S1_GRAD, F1_GRAD, V2_GRAD, S2_GRAD, F2_GRAD, &
       !$OMP&  V3_GRAD, S3_GRAD)
       DO BATCH_START = 1, NX, BS
          BATCHES = BATCHES + 1.0_RT
          BATCH_END = MIN(NX, BATCH_START+BS-1)
          ! Compute the gradient from all data.
          CALL SSE_GRADIENT(X(:,BATCH_START:BATCH_END), &
               Y(:,BATCH_START:BATCH_END), MEAN_SQUARED_ERROR, &
               V1_GRAD, S1_GRAD, F1_GRAD, V2_GRAD, S2_GRAD, F2_GRAD, &
               V3_GRAD, S3_GRAD)
       END DO
       !$OMP END PARALLEL DO
       ! Convert the sum of squared errors into the mean squared error.
       MEAN_SQUARED_ERROR = MEAN_SQUARED_ERROR / RNX
       ! Update the saved "best" model based on error (only when dropout is disabled).
       IF (MEAN_SQUARED_ERROR .LT. BEST_MSE) THEN
          BEST_MSE = MEAN_SQUARED_ERROR
          BEST_V1 = INPUT_VECS
          BEST_S1 = INPUT_SHIFT
          BEST_F1 = INPUT_FLEX
          BEST_V2 = INTERNAL_VECS
          BEST_S2 = INTERNAL_SHIFT
          BEST_F2 = INTERNAL_SHIFT
          BEST_V3 = OUTPUT_VECS
          BEST_S3 = OUTPUT_SHIFT
       END IF
       ! Update the step factor based on model improvement.
       IF (MEAN_SQUARED_ERROR .LT. PREV_MSE) THEN
          STEP_FACTOR = STEP_FACTOR * STEP_COMPOUND_RATE
          STEP_MEAN_CHANGE = STEP_MEAN_CHANGE * STEP_COMPOUND_RATE
          STEP_MEAN_REMAIN = 1.0_RT - STEP_MEAN_CHANGE
          STEP_VAR_CHANGE = STEP_VAR_CHANGE * STEP_COMPOUND_RATE
          STEP_VAR_REMAIN = 1.0_RT - STEP_VAR_CHANGE
       ELSE
          STEP_FACTOR = STEP_FACTOR / STEP_COMPOUND_RATE
          STEP_MEAN_CHANGE = STEP_MEAN_CHANGE / STEP_COMPOUND_RATE
          STEP_MEAN_REMAIN = 1.0_RT - STEP_MEAN_CHANGE
          STEP_VAR_CHANGE = STEP_VAR_CHANGE / STEP_COMPOUND_RATE
          STEP_VAR_REMAIN = 1.0_RT - STEP_VAR_CHANGE
       END IF
       ! Store this error if a record is being kept.
       IF (PRESENT(RECORD)) RECORD(I) = MEAN_SQUARED_ERROR
       ! Store the previous error for tracking the best-so-far.
       PREV_MSE = MEAN_SQUARED_ERROR
       ! Convert the average gradients.
       V1_GRAD = V1_GRAD / BATCHES
       S1_GRAD = S1_GRAD / BATCHES
       F1_GRAD = F1_GRAD / BATCHES
       V2_GRAD = V2_GRAD / BATCHES
       S2_GRAD = S2_GRAD / BATCHES
       F2_GRAD = F2_GRAD / BATCHES
       V3_GRAD = V3_GRAD / BATCHES
       S3_GRAD = S3_GRAD / BATCHES
       ! Update the steps for all different parameters.
       CALL COMPUTE_STEP(SIZE(V1_GRAD(:,:)), V1_GRAD(:,:), V1_MEAN(:,:), V1_VAR(:,:))
       CALL COMPUTE_STEP(SIZE(S1_GRAD(:)), S1_GRAD(:), S1_MEAN(:), S1_VAR(:))
       CALL COMPUTE_STEP(SIZE(F1_GRAD(:)), F1_GRAD(:), F1_MEAN(:), F1_VAR(:))
       DO L = 1, MNS-1
          CALL COMPUTE_STEP(SIZE(V2_GRAD(:,:,L)), V2_GRAD(:,:,L), V2_MEAN(:,:,L), V2_VAR(:,:,L))
          CALL COMPUTE_STEP(SIZE(S2_GRAD(:,L)), S2_GRAD(:,L), S2_MEAN(:,L), S2_VAR(:,L))
          CALL COMPUTE_STEP(SIZE(F2_GRAD(:,L)), F2_GRAD(:,L), F2_MEAN(:,L), F2_VAR(:,L))
       END DO
       CALL COMPUTE_STEP(SIZE(V3_GRAD(:,:)), V3_GRAD(:,:), V3_MEAN(:,:), V3_VAR(:,:))
       CALL COMPUTE_STEP(SIZE(S3_GRAD), S3_GRAD, S3_MEAN, S3_VAR)
       ! Take the gradient steps (based on the computed "step" above).
       INPUT_VECS(:,:)      = INPUT_VECS(:,:)      - V1_GRAD(:,:)
       INPUT_SHIFT(:)       = INPUT_SHIFT(:)       - S1_GRAD(:)
       INPUT_FLEX(:)        = INPUT_FLEX(:)        - F1_GRAD(:)
       INTERNAL_VECS(:,:,:) = INTERNAL_VECS(:,:,:) - V2_GRAD(:,:,:)
       INTERNAL_SHIFT(:,:)  = INTERNAL_SHIFT(:,:)  - S2_GRAD(:,:)
       INTERNAL_FLEX(:,:)   = INTERNAL_FLEX(:,:)   - F2_GRAD(:,:)
       OUTPUT_VECS(:,:)     = OUTPUT_VECS(:,:)     - V3_GRAD(:,:)
       OUTPUT_SHIFT(:)      = OUTPUT_SHIFT(:)      - S3_GRAD(:)

       ! Maintain a constant max-norm across the magnitue of internal vectors.
       DO L = 1, MNS-1
          SCALAR = SQRT(MAXVAL(SUM(INTERNAL_VECS(:,:,L)**2, 1)))
          INTERNAL_VECS(:,:,L) = INTERNAL_VECS(:,:,L) / SCALAR
       END DO

    END DO fit_loop
    ! Restore the best model seen so far (if a reasonable number of steps were taken).
    IF (REVERT_TO_BEST) THEN
       MEAN_SQUARED_ERROR = BEST_MSE
       INPUT_VECS = BEST_V1
       INPUT_SHIFT = BEST_S1
       INTERNAL_VECS = BEST_V2
       INTERNAL_SHIFT = BEST_S2
       OUTPUT_VECS = BEST_V3
       OUTPUT_SHIFT = BEST_S3
    END IF
  CONTAINS

    ! Compute the current step factor as the running average step,
    ! inversely scaled by the average magnitude of the step.
    SUBROUTINE COMPUTE_STEP(N, CURR_STEP, STEP_MEAN, STEP_VAR)
      INTEGER, INTENT(IN) :: N
      REAL(KIND=RT), DIMENSION(1:N), INTENT(INOUT) :: CURR_STEP
      REAL(KIND=RT), DIMENSION(1:N), INTENT(INOUT) :: STEP_MEAN
      REAL(KIND=RT), DIMENSION(1:N), INTENT(INOUT) :: STEP_VAR
      ! Update sliding window mean and variance calculations.
      STEP_MEAN = STEP_MEAN_CHANGE * CURR_STEP                + STEP_MEAN_REMAIN * STEP_MEAN
      STEP_VAR = STEP_VAR_CHANGE * (STEP_MEAN - CURR_STEP)**2 + STEP_VAR_REMAIN * STEP_VAR
      STEP_VAR = MAX(STEP_VAR, EPSILON(STEP_FACTOR))
      ! Take averaged gradient descent step.
      CURR_STEP = STEP_FACTOR * STEP_MEAN
      ! Start scaling by step magnitude by variance once enough data is collected.
      IF (I .GT. MIN_STEPS_TO_STABILITY) THEN
         CURR_STEP = CURR_STEP / SQRT(STEP_VAR)
      END IF
    END SUBROUTINE COMPUTE_STEP

  END SUBROUTINE MINIMIZE_MSE

END MODULE PLRM

max_stars_repo_name: tchlux/util

max_issues_repo_name: tchlux/util

max_forks_repo_name: tchlux/util


module GwfHaloModule

  use KindModule,         only: DP, I4B
  use GwfModule,          only: GwfModelType
  use ConstantsModule,    only: LINELENGTH, DZERO, LENMODELNAME, LENSOLUTIONNAME
  use ConnectionsModule,  only: ConnectionsType
  use MpiExchangeModule,  only: MpiExchangeType
  use SimModule, only: store_error, store_error_unit, ustop

  implicit none

  private
  public :: gwfhalo_cr
  public :: GwfHaloModelType

  type, extends(GwfModelType) :: GwfHaloModelType
    integer(I4B), pointer :: nband
    type(MpiExchangeType), pointer :: MpiSol => null() !PAR
    type(GwfModelType), pointer :: gwf1 => null()
    type(GwfModelType), pointer :: gwf2 => null()
    character(len=LENMODELNAME), pointer :: m1name => null()        ! name of GWF Model 1
    character(len=LENMODELNAME), pointer :: m2name => null()        ! name of GWF Model 2
    integer(I4B), pointer :: nexg => null()
    integer(I4B), pointer :: m1ndim => null()
    integer(I4B), pointer :: m2ndim => null()
    integer(I4B), pointer :: m1irewet => null()
    integer(I4B), pointer :: m2irewet => null()
    integer(I4B), dimension(:), pointer :: nodem1 => null()
    integer(I4B), dimension(:), pointer :: nodem2 => null()
    type(ConnectionsType), pointer :: m1con => null()
    type(ConnectionsType), pointer :: m2con => null()
    integer(I4B), dimension(:), pointer, contiguous :: m1nbnod => null()
    integer(I4B), dimension(:), pointer, contiguous :: m2nbnod => null()
    integer(I4B), pointer :: offset => null()
    integer(I4B), dimension(:), pointer, contiguous :: imapnodem1tohalo => null()
    integer(I4B), dimension(:), pointer, contiguous :: imapnodem2tohalo => null()
    integer(I4B), dimension(:), pointer, contiguous :: imapm1tohalo => null()
    integer(I4B), dimension(:), pointer, contiguous :: imapm2tohalo => null()
    integer(I4B), dimension(:), pointer, contiguous :: m1nodes => null()
    integer(I4B), dimension(:), pointer, contiguous :: m2nodes => null()
  contains

    procedure :: gwfhalo_df1
    procedure :: gwfhalo_df2
    procedure :: gwfhalo_df3
    procedure :: gwfhalo_ar
    procedure :: gwfhalo_cf1
    procedure :: gwfhalo_cf2
    procedure :: gwfhalo_fc_calc
    procedure :: gwfhalo_fn_calc

  end type GwfHaloModelType

  contains

  subroutine gwfhalo_cr(this, id, modelname, m1name, m2name)
! ******************************************************************************
! gwfhalo_cr -- Create a new groundwater flow model object for exchange
! Subroutine: (1) creates model object and add to exchange modellist
!             (2) assign values
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use GwfDisuModule,              only: disu_cr
    use GwfNpfModule,               only: npf_cr
    use Xt3dModule,                 only: xt3d_cr
    use GhostNodeModule,            only: gnc_cr
    use GwfMvrModule,               only: mvr_cr
    use GwfIcModule,                only: ic_cr
    use GwfObsModule,               only: gwf_obs_cr
    use MemoryHelperModule, only: create_mem_path
    ! -- dummy
    type(GwfHaloModelType), pointer   :: this
    integer(I4B), intent(in)      :: id
    character(len=*), intent(in)  :: modelname
    character(len=*), intent(in)  :: m1name
    character(len=*), intent(in)  :: m2name
    ! -- local
    integer(I4B), target :: iout_dum
    integer(I4B), pointer :: in_dum
    ! -- format
! ------------------------------------------------------------------------------
    !
    ! -- Allocate a new GWF Model (this)
    allocate(this)

    ! -- Set memory path before allocation in memory manager can be done
    this%memoryPath = create_mem_path(modelname)

    call this%allocate_scalars(modelname)
    !
    ! -- Assign values
    this%name = modelname
    this%macronym = 'GWF'
    this%id = id
    this%ishalo = .true.
    !
    ! -- TODO: to be replaced by optional arguments
    allocate(in_dum)
    in_dum = -1
    iout_dum = -1
    !
    ! -- Create discretization object
    call disu_cr(this%dis, this%name, in_dum, iout_dum)
    !
    ! -- Create packages that are tied directly to model
    call npf_cr(this%npf, this%name, in_dum, iout_dum)
    call xt3d_cr(this%xt3d, this%name, in_dum, iout_dum)
    call gnc_cr(this%gnc, this%name, in_dum, iout_dum)
    call ic_cr(this%ic, this%name, in_dum, iout_dum, this%dis)
    call mvr_cr(this%mvr, this%name, in_dum, iout_dum, this%npf%dis)
    call gwf_obs_cr(this%obs, in_dum)
    !
    ! -- Set model 2 name
    allocate(this%m1name, this%m2name)
    this%m1name = m1name
    this%m2name = m2name
    !
    ! -- return
    return
  end subroutine gwfhalo_cr

  subroutine gwfhalo_df1(this, m1id, m2id)
! ******************************************************************************
! gwfhalo_df -- define the halo model part 1.
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use MemoryManagerModule,  only: mem_allocate
    use ListsModule,          only: basemodellist
    use BaseModelModule,      only: GetBaseModelFromList, BaseModelType
    use MpiExchangeGwfModule, only: mpi_set_gwfhalo_world_var_int
    ! -- dummy
    class(GwfHaloModelType) :: this
    integer(I4B), intent(in) :: m1id
    integer(I4B), intent(in) :: m2id
    ! -- local
    class(BaseModelType), pointer :: mb
    type(GwfModelType), pointer :: gwf1
    type(GwfModelType), pointer :: gwf2
 ! ------------------------------------------------------------------------------
    !
    ! -- get gwf1
    gwf1 => null()
    mb => GetBaseModelFromList(basemodellist, m1id)
    select type (mb)
    type is (GwfModelType)
      gwf1 => mb
      this%gwf1 => gwf1
    end select
    !
    ! -- get gwf2
    gwf2 => null()
    if (m2id > 0) then
      mb => GetBaseModelFromList(basemodellist, m2id)
      select type (mb)
      type is (GwfModelType)
        gwf2 => mb
        this%gwf2 => gwf2
      end select
    endif
    !
    allocate(this%nband)
    this%nband = 1 ! TODO: read for input file
    !
    call mem_allocate(this%nexg, 'NEXG', this%name)
    call mem_allocate(this%m1ndim, 'M1NDIM', this%name)
    call mem_allocate(this%m2ndim, 'M2NDIM', this%name)
    this%m1ndim = gwf1%dis%ndim
    if (associated(this%gwf2)) then
      this%m2ndim = gwf2%dis%ndim
    else
      ! -- Get NDIM from the global reduction
       this%m2ndim = mpi_set_gwfhalo_world_var_int('DNDIM', this%m2name)
    endif
    !
    ! -- return
    return
  end subroutine gwfhalo_df1

  subroutine gwfhalo_df2(this, nodem1, inewton, m2_bympi)
! ******************************************************************************
! gwfhalo_df -- define the halo model by creating the ia/ja arrays and copying
!   information from each individual model into the halo model
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    ! -- dummy
    class(GwfHaloModelType) :: this
    integer(I4B), dimension(:), intent(in), target :: nodem1
    integer(I4B), intent(in) :: inewton
    logical, intent(in) :: m2_bympi
    ! -- local
! ------------------------------------------------------------------------------
    !
    ! -- Set newton flag
    this%inewton = inewton
    this%npf%inewton = inewton
    !
    this%nexg = size(nodem1)
    !
    ! -- Set pointers
    this%nodem1 => nodem1
    !
    ! -- Set connectivity for M1
    call gwfhalo_df_con(this%m1con, this%gwf1, nodem1, this%m1nbnod, this%nband, &
      this%imapnodem1tohalo, this%imapm1tohalo, this%m1nodes, this%name, 'M1')
    !
    allocate(this%offset)
    this%offset = sum(this%m1nbnod)
    !
    ! -- Allocate for M2
    if (m2_bympi) then
      allocate(this%m2con)
      call this%m2con%allocate_scalars(trim(this%name)//'_M2')
    end if
    !
    ! -- return
    return
  end subroutine gwfhalo_df2

  subroutine gwfhalo_df3(this, nodem2, ihc, cl1, cl2,  &
    hwva)
! ******************************************************************************
! gwfhalo_df -- define the halo model by creating the ia/ja arrays and copying
!   information from each individual model into the halo model
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use ConnectionsModule, only: ConnectionsType, fillisym, filljas
    use SparseModule, only: sparsematrix
    use MemoryManagerModule, only: mem_allocate
    use MpiExchangeModule, only: MpiWorld
    ! -- dummy
    class(GwfHaloModelType) :: this
    integer(I4B), dimension(:), intent(in), target :: nodem2
    integer(I4B), dimension(:), intent(in) :: ihc
    real(DP), dimension(:), intent(in) :: cl1
    real(DP), dimension(:), intent(in) :: cl2
    real(DP), dimension(:), intent(in) :: hwva
    ! -- local
    integer(I4B) :: nexg
    integer(I4B) :: n, m, n1, n2
    integer(I4B) :: iexg, i, j, k
    integer(I4B) :: isym, isymh
    integer(I4B) :: ierror
    type(sparsematrix) :: sparse
    integer(I4B), dimension(:), allocatable :: rowmaxnnz
    !
    ! -- pointers
    type(ConnectionsType), pointer :: con
! ------------------------------------------------------------------------------
    !
    nexg = size(this%nodem1)
    !
    if (associated(this%gwf2)) then
      this%nodem2 => nodem2
      !
      call gwfhalo_df_con(this%m2con, this%gwf2, nodem2, this%m2nbnod, this%nband, &
        this%imapnodem2tohalo, this%imapm2tohalo, this%m2nodes, this%name, 'M2')
    endif
    !
    ! --- Initialize sparse matrix
    this%neq = sum(this%m1nbnod) + sum(this%m2nbnod)
    this%dis%nodes = this%neq
    this%dis%nodesuser = this%neq
    !
    ! -- Allocate the arrays of the halo model disu package
    call this%dis%allocate_arrays()
    !
    ! -- Copy the top, bot, and areas to halo
    call copydbltohalo(this%gwf1%dis%top, this%dis%top, this%imapm1tohalo, 0)
    call copydbltohalo(this%gwf1%dis%bot, this%dis%bot, this%imapm1tohalo, 0)
    call copydbltohalo(this%gwf1%dis%area, this%dis%area, this%imapm1tohalo, 0)
    !
    if (associated(this%gwf2)) then
      call copydbltohalo(this%gwf2%dis%top, this%dis%top, this%imapm2tohalo, this%offset)
      call copydbltohalo(this%gwf2%dis%bot, this%dis%bot, this%imapm2tohalo, this%offset)
      call copydbltohalo(this%gwf2%dis%area, this%dis%area, this%imapm2tohalo, this%offset)
    else
      call MpiWorld%mpicopydbltohalo('HALO_INIT_DIS', 'TOP', 'DIS', this%name, this%dis%top, this%offset)
      call MpiWorld%mpicopydbltohalo('HALO_INIT_DIS', 'BOT', 'DIS', this%name, this%dis%bot, this%offset)
      call MpiWorld%mpicopydbltohalo('HALO_INIT_DIS', 'AREA', 'DIS', this%name, this%dis%area, this%offset)
    endif
    !
    allocate(rowmaxnnz(this%neq))
    do n = 1, this%neq
      rowmaxnnz(n) = 6
    end do
    call sparse%init(this%neq, this%neq, rowmaxnnz)
    !
    ! -- Diagonals for all halo model cells
    do n = 1, this%neq
      call sparse%addconnection(n, n, 1)
    end do
    !
    ! -- Internal connections of M1
    con => this%m1con
    do i = 1, con%nodes
      n1 = i
      do j = con%ia(i) + 1, con%ia(i + 1) - 1
        n2 =  con%ja(j)
        call sparse%addconnection(n1, n2, 1)
      end do
    end do
    ! -- Internal connections of M2
    con => this%m2con
    do i = 1, con%nodes
      n1 = i + this%offset
      do j = con%ia(i) + 1, con%ia(i + 1) - 1
        n2 =  con%ja(j) + this%offset
        call sparse%addconnection(n1, n2, 1)
      end do
    end do
    ! -- Connections between M1 and M2
    do i = 1, size(this%nodem1)
      n1 = this%imapnodem1tohalo(i)
      n2 = this%imapnodem2tohalo(i) + this%offset
      call sparse%addconnection(n1, n2, 1)
      call sparse%addconnection(n2, n1, 1)
    end do
    !
    allocate(this%dis%con)
    call this%dis%con%allocate_scalars(this%name)
    this%dis%con%nodes = this%neq
    this%dis%con%nja = sparse%nnz
    this%dis%con%njas = (this%dis%con%nja - this%dis%con%nodes) / 2
    this%dis%njas = this%dis%con%njas
    !
    ! -- Allocate index arrays of size nja and symmetric arrays
    call this%dis%con%allocate_arrays()
    !
    ! -- Fill the IA and JA arrays from sparse, then destroy sparse
    call sparse%filliaja(this%dis%con%ia, this%dis%con%ja, ierror)
    call sparse%destroy()
    !
    ! -- Create the isym and jas arrays
    call fillisym(this%dis%con%nodes, this%dis%con%nja,                        &
      this%dis%con%ia, this%dis%con%ja, this%dis%con%isym)
    call filljas(this%dis%con%nodes, this%dis%con%nja,                         &
      this%dis%con%ia, this%dis%con%ja, this%dis%con%isym,                     &
      this%dis%con%jas)
    !
    ! -- Fill the ihc, cl1, cl2, and hwva arrays for connections of M1
    con => this%m1con
    do i = 1,con%nodes
      n1 = i
      do j = con%ia(i) + 1, con%ia(i + 1) - 1
        n2 = con%ja(j)
        k = con%getjaindex(n1, n2); isym = con%jas(k)
        k = this%dis%con%getjaindex(n1, n2); isymh = this%dis%con%jas(k)
        this%dis%con%ihc(isymh)  = con%ihc(isym)
        this%dis%con%cl1(isymh)  = con%cl1(isym)
        this%dis%con%cl2(isymh)  = con%cl2(isym)
        this%dis%con%hwva(isymh) = con%hwva(isym)
      end do
    end do
    ! -- Fill the ihc, cl1, cl2, and hwva arrays for connections of M2
    con => this%m2con
    do i = 1,con%nodes
      n1 = i
      do j = con%ia(i) + 1, con%ia(i + 1) - 1
        n2 = con%ja(j)
        k = con%getjaindex(n1, n2); isym = con%jas(k)
        k = this%dis%con%getjaindex(n1 + this%offset, n2 + this%offset)
        isymh = this%dis%con%jas(k)
        this%dis%con%ihc(isymh)  = con%ihc(isym)
        this%dis%con%cl1(isymh)  = con%cl1(isym)
        this%dis%con%cl2(isymh)  = con%cl2(isym)
        this%dis%con%hwva(isymh) = con%hwva(isym)
      end do
    end do
    !
    ! -- Fill the ihc, cl1, cl2, and hwva arrays for connections of M2
    do iexg = 1, nexg
      n = this%imapnodem1tohalo(iexg)
      m = this%imapnodem2tohalo(iexg) + this%offset
      i = this%dis%con%getjaindex(n, m)
      isym = this%dis%con%jas(i)
      this%dis%con%ihc(isym) = ihc(iexg)
      this%dis%con%cl1(isym) = cl1(iexg)
      this%dis%con%cl2(isym) = cl2(iexg)
      this%dis%con%hwva(isym) = hwva(iexg)
    enddo
    !
    ! -- Finish setting up the halo model
    this%nja = this%dis%con%nja
    call this%allocate_arrays()
    !allocate(this%amat(this%nja))
    !
    ! -- Store information needed for observations
    call this%obs%obs_df(this%iout, this%name, 'GWFH', this%dis)
    !
    ! -- return
    return
  end subroutine gwfhalo_df3

  subroutine gwfhalo_df_con(con, gwf, nodem, nbnodes, nband, imapnodemtohalo,   &
    imapmtohalo, mnodes, name, mstr)
! ******************************************************************************
! gwfhalo_df_con -- define the connection.
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use SparseModule, only: sparsematrix
    use MemoryManagerModule, only: mem_allocate
    use ConnectionsModule, only: ConnectionsType, fillisym, filljas

    ! -- dummy
    type(ConnectionsType), intent(out), pointer :: con !out
    type(GwfModelType), intent(in), target :: gwf
    integer(I4B), dimension(:), intent(in), target :: nodem
    integer(I4B), dimension(:), intent(out), pointer, contiguous :: nbnodes !out
    integer(I4B), intent(in) :: nband
    integer(I4B), dimension(:), intent(out), pointer, contiguous :: imapnodemtohalo !out
    integer(I4B), dimension(:), intent(out), pointer, contiguous :: imapmtohalo !out
    integer(I4B), dimension(:), intent(out), pointer, contiguous :: mnodes !out
    character(len=*), intent(in) :: name
    character(len=*), intent(in) :: mstr

    ! -- local
    integer(I4B) :: i, j, n, m, ipos, n1, n2, m1, m2, iband, neq, iact
    integer(I4B) :: isymn, isymm
    integer(I4B) :: ierror, nexg
    type(sparsematrix) :: sparse
    integer(I4B), dimension(:), allocatable :: rowmaxnnz
    integer(I4B), dimension(:), allocatable :: iwrk
! ------------------------------------------------------------------------------
    !
    !
    nexg = size(nodem)
    !
    ! -- Allocate
    call mem_allocate(nbnodes, nband, 'NBNODES', trim(name)//'_'//trim(mstr))
    nbnodes = 0
    call mem_allocate(imapnodemtohalo, nexg, 'IMAPNODEMTOHALO',                    &
       trim(name)//'_'//trim(mstr))
    !
    allocate(iwrk(gwf%neq))
    !
    do iact = 1, 2
      !
      iwrk = 0
      nbnodes = 0
      !
      ! -- Set first band (equals unique nodes in nodem1)
      m = 0
      do i = 1, nexg
        n = nodem(i)
        if (iwrk(n) == 0) then
          m = m + 1
          iwrk(n) = m
          if (iact == 2) then
            imapmtohalo(m) = n
            mnodes(m) = m
          end if
        end if
      end do
      nbnodes(1) = m
      !
      if (iact == 1) then
        do i = 1, nexg
          n = nodem(i)
          imapnodemtohalo(i) = iwrk(n)
        end do
      end if
      !
      ! -- Label the bands
      do iband = 2, nband
        do i= 1, gwf%neq
          n1 = iwrk(i)
          if ((n1 > 0) .and. (n1 <= gwf%neq)) then
            ! -- Label the source band nodes negative
            iwrk(i) = -iwrk(i)
            do ipos = gwf%dis%con%ia(i) + 1, gwf%dis%con%ia(i + 1) - 1
              n = gwf%dis%con%ja(ipos)
              if (iwrk(n) == 0) then
                nbnodes(iband) = nbnodes(iband) + 1
                m = m + 1; n2 = m
                ! -- Label the target band nodes larger than neq
                iwrk(n) = n2 + gwf%neq
                if (iact == 2) then
                  imapmtohalo(m) = n
                  mnodes(m) = m
                end if
              end if
            end do
          end if
        end do
        do i= 1, gwf%neq
          if (iwrk(i) > gwf%neq) then
            iwrk(i) = abs(iwrk(i))
            iwrk(i) = iwrk(i) - gwf%neq
          end if
        end do
      end do ! iband
      !
      do i= 1, gwf%neq
        iwrk(i) = abs(iwrk(i))
      end do
      !
      ! -- Set up the sparse matrix connections
      if (iact == 2) then
        n = 0
        do i = 1, size(imapmtohalo)
          n1 = imapmtohalo(i)
          do ipos = gwf%dis%con%ia(n1) + 1, gwf%dis%con%ia(n1 + 1) - 1
            n2 = gwf%dis%con%ja(ipos)
            if (iwrk(n2) > 0) then
              n = n + 1
              !write(s,'(5(i3.3,a))')&
              !  n, ' --> ',i,'(',n1,') ---> ',iwrk(n2),'(',n2,')'
              !write(*,'(a)') trim(s)
              call sparse%addconnection(i, iwrk(n2), 1)
            end if
          end do
        end do
      end if
      !
      if (iact == 1) then
        !
        ! -- Setup a sparse matrix object in order to determine halo connectivity
        neq = sum(nbnodes)
        allocate(rowmaxnnz(neq))
        do n = 1, neq
          rowmaxnnz(n) = 6
        end do
        call sparse%init(neq, neq, rowmaxnnz)
        !
        ! -- Diagonals for all halo model cells
        do n = 1, neq
          call sparse%addconnection(n, n, 1)
        end do
        !
        ! -- Allocate mapping array
        call mem_allocate(imapmtohalo, neq, 'IMAPMTOHALO', trim(name)//'_'//trim(mstr))
        call mem_allocate(mnodes, neq, 'MNODES', trim(name)//'_'//trim(mstr))
        !
      end if
    end do ! iact
    !
    ! --- Allocate connectivity data
    allocate(con)
    call con%allocate_scalars(trim(name)//'_'//trim(mstr))
    con%nodes = neq
    con%nja = sparse%nnz
    con%njas = (con%nja - con%nodes) / 2
    !
    ! -- Allocate index arrays of size nja and symmetric arrays
    call con%allocate_arrays()
    !
    ! -- Fill the IA and JA arrays from sparse, then destroy sparse
    call sparse%filliaja(con%ia, con%ja, ierror)
    call sparse%destroy()
    !
    ! -- Create the isym and jas arrays
    call fillisym(con%nodes, con%nja, con%ia, con%ja, con%isym)
    call filljas(con%nodes, con%nja, con%ia, con%ja, con%isym, con%jas)
    !
    ! -- Fill the ihc, cl1, cl2, and hwva arrays
    do i = 1, size(imapmtohalo)
      n1 = imapmtohalo(i) ! gwf node
      do ipos = gwf%dis%con%ia(n1) + 1, gwf%dis%con%ia(n1 + 1) - 1
        n2 = gwf%dis%con%ja(ipos) ! gwf node
        if (iwrk(n2) > 0) then
          m1 = i; m2 = iwrk(n2) ! halo node2
          j = gwf%dis%con%getjaindex(n1, n2); isymn = gwf%dis%con%jas(j)
          j = con%getjaindex(m1, m2); isymm = con%jas(j)
          con%ihc(isymm)  = gwf%dis%con%ihc(isymn)
          con%cl1(isymm)  = gwf%dis%con%cl1(isymn)
          con%cl2(isymm)  = gwf%dis%con%cl2(isymn)
          con%hwva(isymm) = gwf%dis%con%hwva(isymn)
        end if
      end do
    end do
    !
    ! -- Cleanup work array
    deallocate(iwrk)
    !
  end subroutine gwfhalo_df_con

  subroutine gwfhalo_ar(this)
! ******************************************************************************
! gwfhalo_ar -- allocate in the halo model, and copy the package information
!   from gwf1 and gwf2 into the halo model
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use MemoryManagerModule, only: mem_allocate, mem_reallocate
    use MpiExchangeGwfModule, only: mpi_set_gwfhalo_world_var_int
    use MpiExchangeModule, only: MpiExchangeType, MpiWorld
    use BndModule, only: BndType, GetBndFromList
    ! -- dummy
    class(GwfHaloModelType) :: this
    ! -- local
    class(BndType), pointer :: packobj
    type(MpiExchangeType), pointer :: mpi
    integer(I4B) :: n, m1flag, m2flag, ip
! ------------------------------------------------------------------------------
    !
    ! -- Create ibound and x arrays for halo model
    call mem_allocate(this%ibound, this%neq, 'IBOUND', this%name)
    call mem_allocate(this%x, this%neq, 'X', this%name)
    !
    ! -- Copy ibound from gwf models into this halo model
    call copyinttohalo(this%gwf1%ibound, this%ibound, this%imapm1tohalo, 0)
    if (associated(this%gwf2)) then
      call copyinttohalo(this%gwf2%ibound, this%ibound, this%imapm2tohalo,     &
        this%offset)
    else
      mpi => this%MpiSol
      call mpi%mpicopyinttohalo('X_IACTIVE', 'IACTIVE', '', this%name,         &
        this%ibound, this%offset)
    endif
    !
    mpi => MpiWorld
    !
    ! -- Set up ic package and copy strt from gwf models into halo model
    call this%ic%allocate_arrays(this%dis%nodes)
    call copydbltohalo(this%gwf1%ic%strt, this%ic%strt, this%imapm1tohalo, 0)
    if (associated(this%gwf2)) then
      call copydbltohalo(this%gwf2%ic%strt, this%ic%strt, this%imapm2tohalo,   &
        this%offset)
    else
      call mpi%mpicopydbltohalo('HALO_INIT_NPFIC', 'STRT', 'IC', this%name,    &
        this%ic%strt, this%offset)
    endif
    do n = 1, this%dis%nodes
      this%x(n) = this%ic%strt(n)
    enddo
    !
    ! -- Copy or point npf pointers
    this%npf%dis    => this%dis
    this%npf%ic     => this%ic
    this%npf%ibound => this%ibound
    this%npf%hnew   => this%x
    !
    ! -- Allocate NPF arrays
    call this%npf%allocate_arrays(this%dis%nodes, this%dis%njas)
    !
    ! -- Copy required NPF arrays into halo model
    call copyinttohalo(this%gwf1%npf%icelltype, this%npf%icelltype,            &
      this%imapm1tohalo, 0)
    if (associated(this%gwf2)) then
      call copyinttohalo(this%gwf2%npf%icelltype, this%npf%icelltype,          &
        this%imapm2tohalo, this%offset)
    else
      call mpi%mpicopyinttohalo('HALO_INIT_NPFIC', 'ICELLTYPE', 'NPF',         &
         this%name, this%npf%icelltype, this%offset)
    endif
    call copydbltohalo(this%gwf1%npf%k11, this%npf%k11, this%imapm1tohalo, 0)
    if (associated(this%gwf2)) then
      call copydbltohalo(this%gwf2%npf%k11, this%npf%k11, this%imapm2tohalo,   &
        this%offset)
    else
      call mpi%mpicopydbltohalo('HALO_INIT_NPFIC', 'K11', 'NPF', this%name,    &
        this%npf%k11, this%offset)
    endif
    !
    call copydbltohalo(this%gwf1%npf%sat, this%npf%sat, this%imapm1tohalo, 0)
    if (associated(this%gwf2)) then
      call copydbltohalo(this%gwf2%npf%sat, this%npf%sat, this%imapm2tohalo,   &
        this%offset)
    else
      call mpi%mpicopydbltohalo('HALO_INIT_NPFIC', 'SAT', 'NPF', this%name,    &
         this%npf%sat, this%offset)
    endif
    !
    ! -- Copy k22
    m1flag = this%gwf1%npf%ik22
    if (.not.associated(this%gwf2)) then
      m2flag = mpi_set_gwfhalo_world_var_int('IK22', this%m2name)
    else
      m2flag = this%gwf2%npf%ik22
    endif
    if ((m1flag /= 0).or.(m2flag /= 0)) then
      this%npf%ik22 = 1
      call mem_reallocate(this%npf%k22, this%dis%nodes, 'K22',                 &
        trim(this%npf%memoryPath))
      if (m1flag /= 0) then
        call copydbltohalo(this%gwf1%npf%k22, this%npf%k22, this%imapm1tohalo, &
          0)
      else
        call copydbltohalo(this%gwf1%npf%k11, this%npf%k22, this%imapm1tohalo, &
          0)
      endif
      if (associated(this%gwf2)) then
        if (m2flag /= 0) then
          call copydbltohalo(this%gwf2%npf%k22, this%npf%k22,                  &
            this%imapm2tohalo, this%offset)
        else
          call copydbltohalo(this%gwf2%npf%k11, this%npf%k22,                  &
            this%imapm2tohalo, this%offset)
        endif
      else
        if (m2flag /= 0) then
          call mpi%mpicopydbltohalo('HALO_INIT_NPFIC', 'K22', 'NPF', this%name,&
            this%npf%k22, this%offset)
        else
          call mpi%mpicopydbltohalo('HALO_INIT_NPFIC', 'K11', 'NPF', this%name,&
            this%npf%k22, this%offset)
        endif
      endif
    endif
    !
    ! -- Copy k33
    m1flag = this%gwf1%npf%ik33
    if (.not.associated(this%gwf2)) then
      m2flag = mpi_set_gwfhalo_world_var_int('IK33', this%m2name)
    else
      m2flag = this%gwf2%npf%ik33
    endif
    !
    if ((m1flag /= 0).or.(m2flag /= 0)) then
      this%npf%ik33 = 1
      call mem_reallocate(this%npf%k33, this%dis%nodes, 'K33',                 &
        trim(this%npf%memoryPath))
      if (m1flag /= 0) then
        call copydbltohalo(this%gwf1%npf%k33, this%npf%k33, this%imapm1tohalo, &
         0)
      else
        call copydbltohalo(this%gwf1%npf%k11, this%npf%k33, this%imapm1tohalo, &
          0)
      endif
      if (associated(this%gwf2)) then
        if (m2flag /= 0) then
          call copydbltohalo(this%gwf2%npf%k33, this%npf%k33,                  &
            this%imapm2tohalo, this%offset)
        else
          call copydbltohalo(this%gwf2%npf%k11, this%npf%k33,                  &
            this%imapm2tohalo, this%offset)
        endif
      else
        if (m2flag /= 0) then
          call mpi%mpicopydbltohalo('HALO_INIT_NPFIC', 'K33', 'NPF', this%name,&
            this%npf%k33, this%offset)
        else
          call mpi%mpicopydbltohalo('HALO_INIT_NPFIC', 'K11', 'NPF', this%name,&
            this%npf%k33, this%offset)
        endif
      endif
    endif
    !
    ! -- Copy angle1
    m1flag = this%gwf1%npf%iangle1
    if (.not.associated(this%gwf2)) then
      m2flag = mpi_set_gwfhalo_world_var_int('IANGLE1', this%m2name)
    else
      m2flag = this%gwf2%npf%iangle1
    endif
    if ((m1flag /= 0).or.(m2flag /= 0)) then
      this%npf%iangle1 = 1
      call mem_reallocate(this%npf%angle1, this%dis%nodes, 'ANGLE1',           &
        trim(this%npf%memoryPath))
      do n = 1, this%dis%nodes
        this%npf%angle1(n) = DZERO
      end do
      if (m1flag/= 0) then
        call copydbltohalo(this%gwf1%npf%angle1, this%npf%angle1,              &
          this%imapm1tohalo, 0)
      endif
      if (associated(this%gwf2)) then
        if (m2flag /= 0) then
          call copydbltohalo(this%gwf2%npf%angle1, this%npf%angle1,            &
            this%imapm2tohalo, this%offset)
        endif
      else
        if (m2flag /= 0) then
          call mpi%mpicopydbltohalo('HALO_INIT_NPFIC', 'ANGLE1', 'NPF',        &
            this%name, this%npf%angle1, this%offset)
        endif
      endif
    endif
    !
    ! -- Copy angle2
    m1flag = this%gwf1%npf%iangle2
    if (.not.associated(this%gwf2)) then
      m2flag = mpi_set_gwfhalo_world_var_int('IANGLE2', this%m2name)
    else
      m2flag = this%gwf2%npf%iangle2
    endif
    if ((m1flag /= 0).or.(m2flag /= 0)) then
      this%npf%iangle2 = 1
      call mem_reallocate(this%npf%angle2, this%dis%nodes, 'ANGLE2',           &
        trim(this%npf%memoryPath))
      do n = 1, this%dis%nodes
        this%npf%angle2(n) = DZERO
      end do
      if (m1flag/= 0) then
        call copydbltohalo(this%gwf1%npf%angle2, this%npf%angle2,              &
          this%imapm1tohalo, 0)
      endif
      if (associated(this%gwf2)) then
        if (m2flag /= 0) then
          call copydbltohalo(this%gwf2%npf%angle2, this%npf%angle2,            &
            this%imapm2tohalo, this%offset)
        endif
      else
        if (m2flag /= 0) then
          call mpi%mpicopydbltohalo('HALO_INIT_NPFIC', 'ANGLE2', 'NPF',        &
            this%name, this%npf%angle2, this%offset)
        endif
      endif
    endif
    !
   ! -- Copy angle3
    m1flag = this%gwf1%npf%iangle3
    if (.not.associated(this%gwf2)) then
      m2flag = mpi_set_gwfhalo_world_var_int('IANGLE3', this%m2name)
    else
      m2flag = this%gwf2%npf%iangle3
    endif
    if ((m1flag /= 0).or.(m2flag /= 0)) then
      this%npf%iangle3 = 1
      call mem_reallocate(this%npf%angle3, this%dis%nodes, 'ANGLE2',           &
        trim(this%npf%memoryPath))
      do n = 1, this%dis%nodes
        this%npf%angle3(n) = DZERO
      end do
      if (m1flag/= 0) then
        call copydbltohalo(this%gwf1%npf%angle3, this%npf%angle3,              &
          this%imapm1tohalo, 0)
      endif
      if (associated(this%gwf2)) then
        if (m2flag /= 0) then
          call copydbltohalo(this%gwf2%npf%angle3, this%npf%angle3,            &
            this%imapm2tohalo, this%offset)
        endif
      else
        if (m2flag /= 0) then
          call mpi%mpicopydbltohalo('HALO_INIT_NPFIC', 'ANGLE3', 'NPF',        &
            this%name, this%npf%angle3, this%offset)
        endif
      endif
    endif
    !
    ! -- prepcheck
    call this%npf%prepcheck()
    !
    ! -- Set rewet flags
    call mem_allocate(this%m1irewet, 'M1IREWET', this%name)
    this%m1irewet = this%gwf1%npf%irewet
    call mem_allocate(this%m2irewet, 'M2IREWET', this%name)
    if (.not.associated(this%gwf2)) then
      this%m2irewet = mpi_set_gwfhalo_world_var_int('IREWET', this%m2name)
    else
      this%m2irewet = this%gwf2%npf%irewet
    endif
    !
    ! -- Copy wetdry
    m1flag = this%m1irewet
    m2flag = this%m2irewet
   if ((m1flag /= 0).or.(m2flag /= 0)) then
      this%npf%irewet = 1
      call mem_reallocate(this%npf%wetdry, this%dis%nodes, 'WETDRY',           &
        trim(this%npf%memoryPath))
      if (m1flag/= 0) then
        call copydbltohalo(this%gwf1%npf%wetdry, this%npf%wetdry,              &
          this%imapm1tohalo, 0)
      endif
      if (associated(this%gwf2)) then
        if (m2flag /= 0) then
          call copydbltohalo(this%gwf2%npf%wetdry, this%npf%wetdry,            &
            this%imapm2tohalo, this%offset)
        endif
      else
        if (m2flag /= 0) then
          call mpi%mpicopydbltohalo('HALO_INIT_NPFIC', 'WETDRY', 'NPF',        &
            this%name, this%npf%wetdry, this%offset)
        endif
      endif
    endif
    !
    ! -- AR for mover packages
    if (.not.associated(this%gwf2)) then
      do ip=1,this%bndlist%Count()
        packobj => GetBndFromList(this%bndlist, ip)
        call packobj%bnd_ar()
      enddo
    endif
    !
    ! -- return
    return
  end subroutine gwfhalo_ar

  subroutine gwfhalo_cf1(this, kiter)
! ******************************************************************************
! gwfhalo_cf1 -- update heads and ibound in halo model
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use MpiExchangeModule, only: MpiExchangeType
    use BndModule, only: BndType,  GetBndFromList
    ! -- dummy
    class(GwfHaloModelType) :: this
    integer(I4B), intent(in) :: kiter
    ! -- local
    type(MpiExchangeType), pointer :: mpi
    class(BndType), pointer :: packobj
    integer(I4B) :: ip
! ------------------------------------------------------------------------------
    !
    ! -- Copy x from gwf models into this halo model x
    call copydbltohalo(this%gwf1%x, this%x, this%imapm1tohalo, 0)
    call copyinttohalo(this%gwf1%ibound, this%ibound, this%imapm1tohalo, 0)
    if (associated(this%gwf2)) then
      call copydbltohalo(this%gwf2%x, this%x, this%imapm2tohalo, this%offset)
      call copyinttohalo(this%gwf2%ibound, this%ibound, this%imapm2tohalo,      &
        this%offset)
    else
      mpi => this%MpiSol
      call mpi%mpicopydbltohalo('X_IACTIVE', 'X', '', this%name,                &
        this%x, this%offset)
      call mpi%mpicopyinttohalo('X_IACTIVE', 'IACTIVE', '', this%name,          &
        this%ibound, this%offset)
    endif
    !
    ! -- return
    return
  end subroutine gwfhalo_cf1

  subroutine gwfhalo_cf2(this, kiter)
! ******************************************************************************
! gwfhalo_cf2 -- rewetting calculations
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    use MpiExchangeModule, only: MpiExchangeType
    use BndModule, only: BndType,  GetBndFromList
    ! -- dummy
    class(GwfHaloModelType) :: this
    integer(I4B), intent(in) :: kiter
    ! -- local
    type(MpiExchangeType), pointer :: mpi
    class(BndType), pointer :: packobj
    integer(I4B) :: ip
! ------------------------------------------------------------------------------
    !
    ! -- rewet
    ! -- TODO: this is a problem.  Wetting and drying prints messages, which
    !    doesn't make sense for the halo model.  And the halo model should
    !    only check for wetting for the n-m connection, and not recheck for
    !    all the cells in the halo model.  This is not good.
    !call this%npf%npf_cf(kiter, this%dis%nodes, this%x)
    !
    ! -- Calculate saturation for convertible cells
    call this%npf%sat_calc(this%x)
    !
    ! -- CF for package movers
    if (.not.associated(this%gwf2)) then
      do ip=1,this%bndlist%Count()
        packobj => GetBndFromList(this%bndlist, ip)
        call packobj%bnd_cf()
      enddo
    endif
    !
    ! -- return
    return
  end subroutine gwfhalo_cf2

!  subroutine gwfhalo_fc(this, kiter, cond)
!! ******************************************************************************
!! gwfhalo_fc -- fill halo model conductance terms into separate amat and rhs
!!   vectors so that the gwf-gwf exchange can add them to the solution amat
!! ******************************************************************************
!!
!!    SPECIFICATIONS:
!! ------------------------------------------------------------------------------
!    ! -- modules
!    use SimModule, only: ustop, store_error, count_errors
!    ! -- dummy
!    class(GwfHaloModelType) :: this
!    integer(I4B), intent(in) :: kiter
!    real(DP), dimension(:), intent(inout) :: cond
!    ! -- local
!    character(len=LINELENGTH) :: errmsg
!    integer(I4B) :: n
!    integer(I4B) :: m
!    integer(I4B) :: nexg
!    integer(I4B) :: ipos
!    real(DP), dimension(:), allocatable :: rhs
!    integer(I4B), dimension(:), allocatable :: idxglo
!! ------------------------------------------------------------------------------
!    !
!    ! -- allocate
!    allocate(rhs(this%dis%nodes))
!    allocate(idxglo(this%nja))
!    !
!    ! -- Copy x from models into this halo model x
!    call copydbltohalo(this%gwf1%x, this%x, this%imapm1tohalo)
!    call copydbltohalo(this%gwf2%x, this%x, this%imapm2tohalo)
!    !
!    ! -- Fill halo amat so that the exchange object can add the terms
!    do n = 1, this%nja
!      this%amat(n) = DZERO
!      idxglo(n) = n
!    enddo
!    do n = 1, this%dis%nodes
!      rhs(n) = DZERO
!    enddo
!    call this%npf%npf_fc(kiter, this%dis%nodes, this%nja, this%nja, this%amat, &
!      idxglo, rhs, this%x)
!    !
!    ! -- TODO: for now just compare halo cond with gwf-gwf cond
!    do nexg = 1, size(this%nodem1)
!      n = this%imapm1tohalo(this%nodem1(nexg))
!      m = this%imapm2tohalo(this%nodem2(nexg))
!      ipos = this%dis%con%getjaindex(n, m)
!      if (this%amat(ipos) - cond(nexg) /= DZERO) then
!        write(errmsg, '(a, 1x, i0, 1x, i0, 2(1pg25.16))') 'HALO MODEL COND /= COND', n, m, &
!          this%amat(ipos), cond(nexg)
!        call store_error(errmsg)
!      endif
!    enddo
!    if (count_errors() > 0) call ustop()
!    !
!    ! -- return
!    return
!  end subroutine gwfhalo_fc

  subroutine gwfhalo_fc_calc(this, iexg, terms)
! ******************************************************************************
! gwfhalo_fc_calc -- calculate the amat and rhs terms for connected gwf cells
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    ! -- dummy
    class(GwfHaloModelType) :: this
    integer(I4B), intent(in) :: iexg
    real(DP), dimension(:, :), intent(inout) :: terms
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: m
    integer(I4B) :: ii
! ------------------------------------------------------------------------------
    n = this%imapnodem1tohalo(iexg)
    m = this%imapnodem2tohalo(iexg) + this%offset
    ii = this%dis%con%getjaindex(n, m)
    call this%npf%npf_fc_calc(n, m, ii, this%x, terms)
    !
    ! -- return
    return
  end subroutine gwfhalo_fc_calc

  subroutine gwfhalo_fn_calc(this, iexg, terms)
! ******************************************************************************
! gwfhalo_fn_calc -- calculate the newton amat and rhs terms for connected gwf cells
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    ! -- dummy
    class(GwfHaloModelType) :: this
    integer(I4B), intent(in) :: iexg
    real(DP), dimension(:, :), intent(inout) :: terms
    ! -- local
    integer(I4B) :: n
    integer(I4B) :: m
    integer(I4B) :: ii
! ------------------------------------------------------------------------------
    n = this%imapnodem1tohalo(iexg)
    m = this%imapnodem2tohalo(iexg) + this%offset
    ii = this%dis%con%getjaindex(n, m)
    call this%npf%npf_fn_calc(n, m, ii, this%x, terms)
    !
    ! -- return
    return
  end subroutine gwfhalo_fn_calc

  subroutine copydbltohalo(gwfarray, gwfhaloarray, imaptohalo, offset)
! ******************************************************************************
! copydbltohalo -- allocate in the halo model
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    ! -- dummy
    real(DP), dimension(:), intent(in) :: gwfarray
    real(DP), dimension(:), intent(inout) :: gwfhaloarray
    integer(I4B), dimension(:), intent(in) :: imaptohalo
    integer, intent(in) :: offset
    ! -- integer
    integer(I4B) :: n, i
! ------------------------------------------------------------------------------
    do i = 1, size(imaptohalo)
      n = imaptohalo(i)
      gwfhaloarray(i+offset) = gwfarray(n)
    enddo
    !
    ! -- return
    return
  end subroutine copydbltohalo

  subroutine copyinttohalo(gwfarray, gwfhaloarray, imaptohalo, offset)
! ******************************************************************************
! copyinttohalo -- allocate in the halo model
! ******************************************************************************
!
!    SPECIFICATIONS:
! ------------------------------------------------------------------------------
    ! -- modules
    ! -- dummy
    integer(I4B), dimension(:), intent(in) :: gwfarray
    integer(I4B), dimension(:), intent(inout) :: gwfhaloarray
    integer(I4B), dimension(:), intent(in) :: imaptohalo
    integer, intent(in) :: offset
    ! -- local
    integer(I4B) :: n, i
! ------------------------------------------------------------------------------
    do i = 1, size(imaptohalo)
      n = imaptohalo(i)
      gwfhaloarray(i+offset) = gwfarray(n)
    enddo
    !
    ! -- return
    return
  end subroutine copyinttohalo
!
end module GwfHaloModule


text:
The socalled Freshman Stripe is the strip of water/mud/grease which can sometimes be seen on the back of bicycle riders after riding over wet roads or through puddles. This road grime also soils backpacks, though its easier to clean off than one might expect. It is sonamed because naive freshmen tend to forget purchasing fenders when they buy their bikes and absentminded salespeople at bicycle shops forget to mention they need them. Used bicycles often lack fenders too. This is more an issue of firsttime bicycle riders than of freshmen.

In addition to the stripe on your back watch out for water thrown back by your front tire when you turn.

see: Bicycling in Winter.

 Avoiding the Stripe

  Buy a rear fender at a Bicycle Shops bike shop.
  Make a Recycling reuse fender out of a 2 liter bottle.
  Wear a waterproof jacket (works better for darker colors).
  Waterproof ponchos allow you to cover both you and your backpack or messenger bag.

Avoiding puddles is pretty easy to do, so I really havent had a problem with the stripe despite not having a fender.  However, every once in a while during the winter of 05 I encountered a big puddle and forgot to slow down.  Lol, safe to say my entire back became soaked with mud and water Users/PatrickSing

Walking into class with a large brown stain on the back of your jeans is not fun. Users/MattNagel

When I was a freshman, I found the biggest problem wasnt your own bike, it was everyone elses.  I had a fender that worked well to prevent my  own stripe....its the amount of people you bike past or around or even near that get you.  Either get side splashed, or whoevers in front tosses it behind them as well. Plays havoc with your pants. End up getting all around muddy, even with a rain jacket.  My back however, was always clean. Users/EdWins ES

That waterproof jacket solution would probably only work if you dont have a backpack on Users/AmyZimmerman

I think the best thing to do in the rain is just allow extra time, and walk with an umbrella.  i find it uncomfortable to bike when your forward movement causes rain to lang on your face.  But another thing to do is just not care, perhaps the stripe will be in next winter :)  Users/MattHh

20060205 16:33:47 nbsp if you go slow enough, you can avoid the stripe.  just plan ahead and give yourself some extra time to get to class. Users/PaulHarms

20060910 00:32:08 nbsp At least the stripe shows you rode your bike to campus instead of being lazy and driving. Users/JimEvans

20070429 20:06:49 nbsp I wear a complete rainsuit when I bike in the rain, stripe be damned. Also a couple of my bikes have winnar fenders. Users/StevenDaubert

20070505 18:27:26 nbsp Haha what stupid freshmen! Users/AlexNelson

20070510 14:03:16 nbsp If you need big ponchos (raincovers) that you can wear with a backpack, check out http://www.jeantex.com Its a german company with First Class bike (and sail) clothing and goodies. I have been testing the Jeantex Bergerac jacket on my bicycles (I am carless as well, but I live in NL) for over 6 months and it outperforms my North Face clothes by orders of magnitude. And jeantex is dirt cheap. Nice site! Users/JanVerhoeven

20071025 01:20:53 nbsp There should be instructions here on how to make a fender from a 2 liter bottle.  Users/BreeButler

20071025 04:27:09 nbsp well you just slit it down one side lengthwise. . .im not sure people NEED a diagram. but youre welcome to add one. Users/KellyCorcoran

20090219 17:49:27 nbsp Just a note on my edit of this page. First, thanks to the original author. Its a great page. But the previous description that bike store salespeople trick freshmen or beginning bikers to purchase bikes without fenders isnt just erroneous but foolhardy. Its erroneous because it just isnt true. Its foolhardy because if you actually walk into a bike store to purchase a bike youll find that sales associate are adamant you purchase fenders ... and for good reason! Its something else you purchase from them. Plus, the margins retailers make on accessories, like fenders, is around 4050%, whereas the margins they make on the bikes is almost nil. There are very few bike shops in town that will let you walk out of the store without encouraging you to purchase fenders. The freshman stripe is more often the result of someone purchasing a used bike, bringing a bike from a home, or refusing to pay for fenders after the sticker shock of the new bike. Users/condemned2bfree

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
subroutine postpr_hdt(nst, gdp)
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: postpr_hdt.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io/src/output/postpr_hdt.f90 $
!!--description-----------------------------------------------------------------
!
!    Function:
!              - Call POSTPR without having to define all arguments in the
!                main routine. This routine is called in TRISOL for half
!                time steps.
!
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use flow2d3d_timers
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer                             , pointer :: lundia
    integer                             , pointer :: lunprt
    real(fp)                            , pointer :: rhow
    integer                             , pointer :: ktemp
    real(fp)                            , pointer :: grdang
    character(19)                       , pointer :: prsmap
    character(21)                       , pointer :: selmap
    character(23)                       , pointer :: prshis
    character(23)                       , pointer :: selhis
    character(256)                      , pointer :: comfil  !!  Communication file name
    character(256)                      , pointer :: runid   !!  Run identification code for the current simulation (used to determine the names of the in- /output files used by the system)
    character(256)                      , pointer :: trifil  !!  File name for TRISULA NEFIS output files (tri"h/m"-"casl""labl".dat/def)
    character(5)                        , pointer :: versio  !!  Version nr. of the current package
    integer                             , pointer :: initi   ! Control parameter
    integer                             , pointer :: iphisc  ! Current time counter for printing history data
    integer                             , pointer :: itcomc  ! Current time counter for the communication file
    integer                             , pointer :: itcur   ! Current time counter for the communication file, where starting point depend on CYCLIC
    integer                             , pointer :: itdroc  ! Current time counter for the drogue data file
    integer                             , pointer :: ithisc  ! Current time counter for the history file
    integer                             , pointer :: itimc   ! Current time step counter for 2D system
    integer                             , pointer :: itmapc  ! Current time counter for the map file
    integer                             , pointer :: itrstc  ! Current time counter for the restart file. Start writing after first interval is passed. Last time will always be written to file for ITRSTI > 0
    integer                             , pointer :: npmap   ! Current array counter for printing map data
    integer                             , pointer :: ntcur   ! Total number of timesteps on comm. file (to write to)
    real(fp)                            , pointer :: dtsec   ! DT in seconds
!
! Call variables
!
    integer                                       :: nst           ! Current time step counter
!
! Local variables
!
    logical                                       :: error         ! Flag=TRUE if an error is encountered
!
!! executable statements -------------------------------------------------------
!
    lundia              => gdp%gdinout%lundia
    lunprt              => gdp%gdinout%lunprt
    rhow                => gdp%gdphysco%rhow
    ktemp               => gdp%gdtricom%ktemp
    grdang              => gdp%gdtricom%grdang
    prsmap              => gdp%gdtricom%prsmap
    selmap              => gdp%gdtricom%selmap
    prshis              => gdp%gdtricom%prshis
    selhis              => gdp%gdtricom%selhis
    comfil              => gdp%gdtricom%comfil
    runid               => gdp%runid
    trifil              => gdp%gdtricom%trifil
    versio              => gdp%gdtricom%versio
    initi               => gdp%gdtricom%initi
    iphisc              => gdp%gdtricom%iphisc
    itcomc              => gdp%gdtricom%itcomc
    itcur               => gdp%gdtricom%itcur
    itdroc              => gdp%gdtricom%itdroc
    ithisc              => gdp%gdtricom%ithisc
    itimc               => gdp%gdtricom%itimc
    itmapc              => gdp%gdtricom%itmapc
    itrstc              => gdp%gdtricom%itrstc
    npmap               => gdp%gdtricom%npmap
    ntcur               => gdp%gdtricom%ntcur
    dtsec               => gdp%gdtricom%dtsec
    !
    ! Call postpr
    !
    error = .false.
    !
    call psemnefis
    call timer_start(timer_postpr, gdp)
    call postpr(lundia    ,lunprt    ,error     ,versio    ,comfil    , &
              & trifil    ,runid     ,prsmap    ,prshis    ,selmap    , &
              & selhis    ,rhow      ,grdang    ,initi     ,dtsec     , &
              & nst+1     ,iphisc    ,npmap     ,itcomc    ,itimc     , &
              & itcur     ,ntcur     ,ithisc    ,itmapc    ,itdroc    , &
              & itrstc    ,ktemp     ,.true.    ,gdp       )
    call timer_stop(timer_postpr, gdp)
    call vsemnefis
    !
    !if (error) return
  end subroutine postpr_hdt

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
#
# test math operators and order of operations
#

print("Testing math operators and order of operations...")
errors = 0

try assert 1 + 2 * 3 == 7 rescue errors += 1
try assert (1 + 2) * 3 == 9 rescue errors += 1
try assert (1 - 2) * 3 == -3 rescue errors += 1
try assert 1 - 2 * 3 == -5 rescue errors += 1
try assert (3 * 3) ^ 2 == 81 rescue errors += 1
try assert 3 * 3 ^ 2 == 27 rescue errors += 1
try assert 4 * (3 + 3) == 24 rescue errors += 1
try assert 4 / 2 + 6 == 8 rescue errors += 1
try assert 4 / (2 + 6) == 0.5 rescue errors += 1
try assert 60 % 21 // 2 == 9 rescue errors += 1
try assert 63 % (21 // 2) == 3 rescue errors += 1

a = 8
assert (a += 1) == 9
assert (a -= 1) == 8
assert (a /= 2) == 4.0
assert (a //= 2) == 2
assert (a %= 2) == 0

print(if errors "Some tests failed ({errors}) errors\n" else "passed!\n")

max_stars_repo_name: cwells/franz

max_issues_repo_name: cwells/franz

max_forks_repo_name: cwells/franz


text:
#include "f90tw_defs_boost.h"

module assertions_boost

      use ISO_C_BINDING

      implicit none

#define HEADERLEVEL "boost_level.h"
#define HEADERLEVELWHAT "interface_boost_wrappers.h"
#define VERSIONSUB F90_F90TW_VER
#define BOOSTFUNCTION0 F90_BOOST_LEVEL
#define BOOSTFUNCTION1 F90_BOOST_LEVEL_MESSAGE
#define BOOSTFUNCTION2 F90_BOOST_LEVEL_CHECK1T2A
#define BOOSTFUNCTION3 F90_BOOST_LEVEL_CHECK1T3A
#define BOOSTSOMETHING F90_BOOST_SOMETHING
#define BOOSTVOID F90_BOOST_VOID

#define DOINTERFACE
#include "boost_what.h"

contains

#undef DOINTERFACE
#define CAST CAST_IMPLEMENT
#define DIM1 DIM1_IMPLEMENT
#define DIM2 DIM2_IMPLEMENT
#include "boost_what.h"

#undef BOOSTVOID
#undef BOOSTSOMETHING
#undef BOOSTFUNCTION3
#undef BOOSTFUNCTION2
#undef BOOSTFUNCTION1
#undef BOOSTFUNCTION0
#undef VERSIONSUB
#undef HEADERLEVELWHAT
#undef HEADERLEVEL

end module assertions_boost

ext: fpp


text:
*DECK DEXPRL
      DOUBLE PRECISION FUNCTION DEXPRL (X)
C***BEGIN PROLOGUE  DEXPRL
C***PURPOSE  Calculate the relative error exponential (EXP(X)-1)/X.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C4B
C***TYPE      DOUBLE PRECISION (EXPREL-S, DEXPRL-D, CEXPRL-C)
C***KEYWORDS  ELEMENTARY FUNCTIONS, EXPONENTIAL, FIRST ORDER, FNLIB
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Evaluate  EXPREL(X) = (EXP(X) - 1.0) / X.   For small ABS(X) the
C Taylor series is used.  If X is negative the reflection formula
C         EXPREL(X) = EXP(X) * EXPREL(ABS(X))
C may be used.  This reflection formula will be of use when the
C evaluation for small ABS(X) is done by Chebyshev series rather than
C Taylor series.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH
C***REVISION HISTORY  (YYMMDD)
C   770801  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   890911  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C***END PROLOGUE  DEXPRL
      DOUBLE PRECISION X, ABSX, ALNEPS, XBND, XLN, XN,  D1MACH
      LOGICAL FIRST
      SAVE NTERMS, XBND, FIRST
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DEXPRL
      IF (FIRST) THEN
         ALNEPS = LOG(D1MACH(3))
         XN = 3.72D0 - 0.3D0*ALNEPS
         XLN = LOG((XN+1.0D0)/1.36D0)
         NTERMS = XN - (XN*XLN+ALNEPS)/(XLN+1.36D0) + 1.5D0
         XBND = D1MACH(3)
      ENDIF
      FIRST = .FALSE.
C
      ABSX = ABS(X)
      IF (ABSX.GT.0.5D0) DEXPRL = (EXP(X)-1.0D0)/X
      IF (ABSX.GT.0.5D0) RETURN
C
      DEXPRL = 1.0D0
      IF (ABSX.LT.XBND) RETURN
C
      DEXPRL = 0.0D0
      DO 20 I=1,NTERMS
        DEXPRL = 1.0D0 + DEXPRL*X/(NTERMS+2-I)
 20   CONTINUE
C
      RETURN
      END


text:
        SUBROUTINE GETGIGSTATUS(ifail,
     &                          SENSITIVITY,
     &                          NTGRESULT,
     &                          GIGRESULT,
     &                          NRESULTS,
     &                          LEVEL,
     &                          PROD,
     &                          INGR,
     &                          NUMVIOLATED,
     &                          VIOLTYPE,
     &                          VIOLPROD,
     &                          VIOLNTRM,
     &                          VIOLVALUE)
*
*       +---------------------------------------------------+
*       | Copyright (c) FORMAT Int '03. All rights reserved |
*       +---------------------------------------------------+
*       | Function   : GETGIGSTATUS                         |
*       | Description: Set the STATUS for GIGs              |
*       | Revision   : November 2003                        |
*       +---------------------------------------------------+
*
*------------------------
*  Implicit definition  *
*------------------------
*
        IMPLICIT INTEGER (A-Z)
*
*-------------------
*  Include common  *
*-------------------
*
        INCLUDE 'SLPCOMM.INC'
        INCLUDE 'pusr.inc'
*
*------------------
*  Indirect data  *
*------------------
*
        DIMENSION        NTGRESULT(*),
     &                   VIOLTYPE (*),
     &                   VIOLPROD (*),
     &                   VIOLNTRM (*),
     &                   VIOLVALUE(*)
*
        DOUBLE PRECISION SENSITIVITY,NTGRESULT,LEVEL,VIOLVALUE
*
*---------------
*  Local data  *
*---------------
*
        INTEGER GIGINDEX
        DOUBLE PRECISION DLOWER,DUPPER,DIFFLEVEL
*
*---------------
*  Initialise  *
*---------------
*
        TEMPSTATUS=0
*
*--------------------------------------------
*  Get Min and Max constraints for the GIG  *
*--------------------------------------------
*
        NDX=GIGINDEX(GIGRESULT,2,nu,rm,mi,de,NUM_RATIO_RESULTS,NRESULTS)
        DLOWER=NTGRESULT(NDX)
*
        NDX=GIGINDEX(GIGRESULT,3,nu,rm,mi,de,NUM_RATIO_RESULTS,NRESULTS)
        DUPPER=NTGRESULT(NDX)
*
*----------------------------------
*  Correct the level for the GIG  *
*----------------------------------
*
        DIFFLEVEL=0D0
*
*  Check if the GIG is violated.
*
        DO I=1,NUMVIOLATED
          IF(VIOLTYPE(I).EQ.TYPE_NU.AND.
     &       VIOLPROD(I).EQ.INGR.AND.
     &       VIOLNTRM(I).EQ.PROD) THEN
            DIFFLEVEL=VIOLVALUE(I)
          END IF
        END DO
*
        LEVEL=LEVEL+DIFFLEVEL
*
*-----------------------------
*  Set STATUS when feasible  *
*-----------------------------
*
        IF(ifail.EQ.0) THEN
          IF(SENSITIVITY.EQ.0D0) THEN
*  Within limits
            TEMPSTATUS=1
          ELSE IF(SENSITIVITY.GT.0D0) THEN
*  On minimum limit
            TEMPSTATUS=2
          ELSE IF(SENSITIVITY.LT.0D0) THEN
*  On maximum limit
            TEMPSTATUS=3
          END IF
*
*  In some cases the above method is not correct.....
*
          CALL COMPAREFLOATINGNUMBERS(DLOWER,DUPPER,ISIGN)
*
*  If min/max constraints are different, correct where level equals min/max
*
          IF(ISIGN.NE.0D0) THEN
            CALL COMPAREFLOATINGNUMBERS(LEVEL,DLOWER,ISIGN)
            IF(ISIGN.EQ.0D0) THEN
              TEMPSTATUS=2
            END IF
            CALL COMPAREFLOATINGNUMBERS(LEVEL,DUPPER,ISIGN)
            IF(ISIGN.EQ.0D0) THEN
              TEMPSTATUS=3
            END IF
          END IF
        ELSE
*
*-----------------------------
*  Set STATUS if infeasible  *
*-----------------------------
*
          IF(SENSITIVITY.EQ.0D0) THEN
*  Within limits
            TEMPSTATUS=1
          ELSE
*
*  Put the level into RESULT
*
            NDX=GIGINDEX(GIGRESULT,1,nu,rm,mi,de,NUM_RATIO_RESULTS,
     &                   NRESULTS)
            NTGRESULT(NDX)=LEVEL
*
            IF(SENSITIVITY.GT.0D0) THEN
*  On minimum limit
              TEMPSTATUS=2
*
            ELSE IF(SENSITIVITY.LT.0D0) THEN
*  On maximum limit
              TEMPSTATUS=3
            END IF
*
*  In all cases
*
            CALL COMPAREFLOATINGNUMBERS(LEVEL,DLOWER,ISIGN)
*  Under minimum limit
            IF(ISIGN.LT.0) THEN
              TEMPSTATUS=4
            END IF
            CALL COMPAREFLOATINGNUMBERS(LEVEL,DUPPER,ISIGN)
*  Above maximum limit
            IF(ISIGN.GT.0) THEN
              TEMPSTATUS=5
            END IF
          END IF
        END IF
*
*  Put the status into RESULT
*
        IF(TEMPSTATUS.EQ.0) THEN
          write(*,*)'FATAL ERROR !!. STATUS not found for a GIG!'
          STOP
        END IF
*
        NDX=GIGINDEX(GIGRESULT,8,nu,rm,mi,de,NUM_RATIO_RESULTS,NRESULTS)
        NTGRESULT(NDX)=TEMPSTATUS
*
*  Put the level into RESULT
*
        NDX=GIGINDEX(GIGRESULT,1,nu,rm,mi,de,NUM_RATIO_RESULTS,NRESULTS)
        NTGRESULT(NDX)=LEVEL

*
        RETURN
*
        END

        SUBROUTINE GETRATIOSTATUS(ifail,
     &                            SENSITIVITY,
     &                            NTGRESULT,
     &                            RATIORESULT,
     &                            NRESULTS)
*
*       +---------------------------------------------------+
*       | Copyright (c) FORMAT Int '03. All rights reserved |
*       +---------------------------------------------------+
*       | Function   : GETRATIOSTATUS                       |
*       | Description: Set the STATUS for ratios            |
*       | Revision   : November 2003                        |
*       +---------------------------------------------------+
*
*------------------------
*  Implicit definition  *
*------------------------
*
        IMPLICIT INTEGER (A-Z)
*
*-------------------
*  Include common  *
*-------------------
*
        INCLUDE 'SLPCOMM.INC'
        INCLUDE 'pusr.inc'
*
*------------------
*  Indirect data  *
*------------------
*
        DIMENSION        NTGRESULT(*)
*
        DOUBLE PRECISION SENSITIVITY,NTGRESULT
*
*---------------
*  Local data  *
*---------------
*
        DOUBLE PRECISION LEVEL,DLOWER,DUPPER,LEVEL1,LEVEL2,DIFFLEVEL
*
*-----------------------------
*  Set STATUS when feasible  *
*-----------------------------
*
        IF(ifail.EQ.0) THEN
          IF(SENSITIVITY.EQ.0D0) THEN
*  Within limits
            TEMPSTATUS=1
          ELSE IF(SENSITIVITY.GT.0D0) THEN
*  On minimum limit
            TEMPSTATUS=2
          ELSE IF(SENSITIVITY.LT.0D0) THEN
*  On maximum limit
            TEMPSTATUS=3
          END IF
          NDX=RATIOINDEX(RATIORESULT,8,nu,rm,mi,de,NRESULTS)
          IF(TEMPSTATUS.GT.NTGRESULT(NDX)) THEN
            NTGRESULT(NDX)=TEMPSTATUS
          END IF
*  Set sensitivity
          NDX=RATIOINDEX(RATIORESULT,4,nu,rm,mi,de,NRESULTS)
          NTGRESULT(NDX)=SENSITIVITY
*
*-----------------------------
*  Set STATUS if infeasible  *
*-----------------------------
*
        ELSE
*  Set odd status and let Integra-Mix sort the status out
          NDX=RATIOINDEX(RATIORESULT,8,nu,rm,mi,de,NRESULTS)
          NTGRESULT(NDX)=9999D0
*  Set sensitivity
          NDX=RATIOINDEX(RATIORESULT,4,nu,rm,mi,de,NRESULTS)
          NTGRESULT(NDX)=SENSITIVITY
        END IF
        RETURN
*
        END

max_stars_repo_name: ERGO-Code/FortranFilterSLP

max_issues_repo_name: ERGO-Code/FortranFilterSLP

max_forks_repo_name: ERGO-Code/FortranFilterSLP


text:
MODULE EOS

  IMPLICIT NONE

  INTEGER(KIND=4), PARAMETER :: NUMGASINFOS  = 4
  INTEGER(KIND=4), PARAMETER :: NUMNONDIMENS = 16
  INTEGER(KIND=4), PARAMETER :: NUMGASPARAMS = 5

  INTERFACE

     SUBROUTINE GASDV(numDim,eosInfo,eosParams,eosNonDim,qIn,dvIn)
       IMPLICIT NONE
       INTEGER(4), INTENT(IN)  :: numDim
       INTEGER(8), INTENT(IN)  :: eosInfo
       REAL(8),    INTENT(IN)  :: eosParams
       REAL(8),    INTENT(IN)  :: eosNonDim
       REAL(8),    INTENT(IN)  :: qIn
       REAL(8),    INTENT(OUT) :: dvIn
     END SUBROUTINE GASDV

     SUBROUTINE GASDP(numDim,eosInfo,eosParams,eosNonDim,qIn,dvIn,dp)
       IMPLICIT NONE
       INTEGER(4), INTENT(IN)  :: numDim
       INTEGER(8), INTENT(IN)  :: eosInfo
       REAL(8),    INTENT(IN)  :: eosParams
       REAL(8),    INTENT(IN)  :: eosNonDim
       REAL(8),    INTENT(IN)  :: qIn
       REAL(8),    INTENT(IN)  :: dvIn
       REAL(8),    INTENT(OUT) :: dp
     END SUBROUTINE GASDP

     SUBROUTINE GASISOTHERMAL(numDim,eosInfo,eosParams,eosNonDim,isoT,q)
       IMPLICIT NONE
       INTEGER(4), INTENT(IN)    :: numDim
       INTEGER(8), INTENT(IN)    :: eosInfo
       REAL(8),    INTENT(IN)    :: eosParams
       REAL(8),    INTENT(IN)    :: eosNonDim
       REAL(8),    INTENT(IN)    :: isoT
       REAL(8),    INTENT(INOUT) :: q
     END SUBROUTINE GASISOTHERMAL

     SUBROUTINE GASEOSROE(numDim,eosInfo,eosParams,eosNonDim,qIn,dv,dp,c)
       IMPLICIT NONE
       INTEGER(4), INTENT(IN)  :: numDim
       INTEGER(8), INTENT(IN)  :: eosInfo
       REAL(8),    INTENT(IN)  :: eosParams
       REAL(8),    INTENT(IN)  :: eosNonDim
       REAL(8),    INTENT(IN)  :: qIn
       REAL(8),    INTENT(OUT) :: dv
       REAL(8),    INTENT(OUT) :: dp
       REAL(8),    INTENT(OUT) :: c
     END SUBROUTINE GASEOSROE

     SUBROUTINE GASSOUNDSPEED(numDim,eosInfo,eosParams,eosNonDim,qIn,dvIn,dpIn,c)
       IMPLICIT NONE
       INTEGER(4), INTENT(IN)  :: numDim
       INTEGER(8), INTENT(IN)  :: eosInfo
       REAL(8),    INTENT(IN)  :: eosParams
       REAL(8),    INTENT(IN)  :: eosNonDim
       REAL(8),    INTENT(IN)  :: qIn
       REAL(8),    INTENT(IN)  :: dvIn
       REAL(8),    INTENT(IN)  :: dpIn
       REAL(8),    INTENT(OUT) :: c
     END SUBROUTINE GASSOUNDSPEED

  END INTERFACE

  CONTAINS

    !> Calculate average material properties for mixture
    SUBROUTINE MIXPROPERTIES(numPointsBuffer,numSpecies,pointIndex,speciesMolWt, &
         speciesHf,speciesCp,univR,gasDensity,rhoY,tMix)

      IMPLICIT NONE

      INTEGER(KIND=8) :: numPointsbuffer,pointIndex
      INTEGER(KIND=4) :: numSpecies
      REAL(KIND=8)    :: speciesMolWt(numSpecies)
      REAL(KIND=8)    :: speciesHf(numSpecies)
      REAL(KIND=8)    :: speciesCp(numSpecies)
      REAL(KIND=8)    :: gasDensity(numPointsBuffer)
      REAL(KIND=8)    :: rhoY(numSpecies*numPointsBuffer)
      REAL(KIND=8)    :: tMix(3), univR

      ! Local
      REAL(KIND=8)    :: bufWtM1, baseR, baseHf, baseCp, rhoM1, specFrac
      INTEGER(KIND=8) :: specIndex
      INTEGER(KIND=4) :: iSpec

      bufWtM1 = 1.0_8/speciesMolWt(numSpecies)

      baseR   = univR*bufWtM1
      baseHf  = speciesHf(numSpecies)
      baseCp  = speciesCp(numSpecies)

      tMix(1) = baseR
      tMix(2) = baseHf
      tMix(3) = baseCp

      rhoM1   = 1.0_8/gasDensity(pointIndex)

      DO iSpec = 1, numSpecies-1
         specIndex   = (iSpec-1)*numPointsBuffer + pointIndex
         specFrac    = rhoY(specIndex)*rhoM1
         tMix(1) = tMix(1) + specFrac*(univR/speciesMolWt(iSpec) - baseR)
         tMix(2) = tMix(2) + specFrac*(speciesHf(iSpec) - baseHf)
         tMix(3) = tMix(3) + specFrac*(speciesCp(iSpec) - baseCp)
      END DO

    END SUBROUTINE MIXPROPERTIES

END MODULE EOS

max_stars_repo_name: MTCam/OKLib

max_issues_repo_name: MTCam/OKLib

max_forks_repo_name: MTCam/OKLib