      MODULE LOADR_N
      
      CONTAINS
      
      SUBROUTINE LOADER_N (HMAX, KMAX, MMAX, KMQ_MAX, MQ_MAX, MINE, DYSTAR, SCAPSUM, 
     & MINSTATE, MAXSTATE, MMI1, MMI2, MMI3, MSC, IYMIN, IYMAX, IY, WT, Y, NR, 
     & NPATTERNS, PATTERN, NRUNS, RUN, DERIVS, AM, CM, COL_LBL, ROW_LBL, VOCAL)
      USE CODE ! CONTAINS "ENCODEC" & "DECODEC"
      USE CONSTANTS
      USE DERIVED_TYPES
      USE EM_IMSL ! CONTAINS EMULATIONS OF IMSL "ANORDF" AND "ANORPR" FUNCTIONS
      USE FILES
      USE GAUSS ! CONTAINS "GPDF", "GPMF", AND "Z"
      USE PATTERN_DATA
      USE UTILS ! CONTAINS "PRN_RMAT" AND "PRN_IMAT"
      IMPLICIT NONE
C *** 
C *** PURPOSE: THIS IS THE PROTOTYPE "LOADER" ROUTINE FOR "COGNIT_N".
C ***
C *** METHODOLOGY: LOADS PMF VALUES INTO THE "A" AND "C" SUB-MATRICES OF THE MARKOV CHAIN
C *** TRANSITION MATRIX.  THE APPROACH USED HERE IS INTENDED TO FACILITATE CALCULATION OF
C *** PARTIALS OF "HY" AND "HFY" W/RT OTHER VARIABLES ($y_1$, $A$, $F$, $%SIG$, ETC.).
C *** NB:  FOR STEADY-STATE FORMULATION OF TRANSITION MATRICES, SEE $\citep{MCCU1998}$.
C ***
C *** MODIFICATION HISTORY:
C *** 1) 7/9/18: REVISED TO CORRECT AND EXTEND LOGIC RE SUPRESSION OF DERIVATIVES
C *** CALCULATION (DERIVS_NE_0 LOGIC) FOR UNIFORM SCDs.
C *** 2) 6/16/19: REVISED TO AVOID "CRASH POINT" ASSOCIATED WITH CALCULATION OF "HY".
C *** SEE PUA/PLA LOGIC BELOW.
C *** 3) 8/2/19: UNDID MODIFICATION 2 ABOVE, AND STREAMLINED CALCULATION OF BOTH "HY"
C *** AND "HFY".  CONDUCTED MINOR CODE CLEANUP.
C *** 4) 2/7/20: STARTED REVISIONS TO CONVERT "HY/HFY" CALCULATIONS TO USING MMI3 DATA,
C *** VICE MMI2
C ***
C *** TODO:  SEE "TODO" RELATED TO TEST "HY > HY_MIN"
C ***
C *** ------------------------------------------------------------------
C ***
C *** ARGUMENTS:
C ***
      INTEGER, INTENT(IN) :: HMAX     ! NUMBER OF CHANNEL SEGMENTS
      INTEGER, INTENT(IN) :: KMAX     ! NUMBER OF SWEEPER TYPES
      INTEGER, INTENT(IN) :: MMAX     ! NUMBER OF MINE TYPES
      INTEGER, INTENT(IN) :: KMQ_MAX  ! KMAX*MMAX*QMAX, FOR USE BY "ENCODEC"/"DECODEC"
      INTEGER, INTENT(IN) :: MQ_MAX   ! MMAX*QMAX, FOR USE BY "ENCODEC"/"DECODEC"
      TYPE (MINES), INTENT(IN) :: MINE(HMAX,MMAX) ! MINE CHARACTERISTICS
      REAL, INTENT(IN) :: DYSTAR      ! INCREMENT USED IN NUMERICAL INTEGRATION
      INTEGER, INTENT(IN) :: SCAPSUM(HMAX,0:MMAX) ! TABLE OF CUMULATIVE "SCAP" VALUES
      INTEGER, INTENT(IN) :: MINSTATE ! MINIMUM VALUE OF STATE INDEX $s$
      INTEGER, INTENT(IN) :: MAXSTATE ! MAXIMUM VALUE OF STATE INDEX $s$
      TYPE (MMI1S), INTENT(IN) :: MMI1(HMAX,KMAX,MMAX) ! GENERAL MINESWEEPER-MINE
C                                     ! INTERACTION DATA
C !!!    NOTE THAT MMI2 AND MMI3 ARE BOTH HARD-WIRED FOR MAXIACTMAX = 1 !
      TYPE (MMI2S), INTENT(IN) :: MMI2(HMAX,KMAX,MMAX,0:1,1) ! DETAILED MINESWEEPER-MINE
C                                     ! INTERACTION DATA IN RECTANGULAR (A/B) FORMAT
      TYPE (MMI3S), INTENT(IN) :: MMI3(HMAX,KMAX,MMAX,0:1,1,3) ! DETAILED
                                      ! MINESWEEPER-MINE INTERACTION DATA IN PIECEWISE
                                      ! LINEAR FORMAT
      TYPE (MSCS), INTENT(IN) :: MSC(HMAX,KMAX) ! MINESWEEPER CHARACTERISTICS
      INTEGER, INTENT(IN) :: IYMIN    ! LOWER INDEX ON LATERAL LOCATIONS
      INTEGER, INTENT(IN) :: IYMAX    ! UPPER INDEX ON LATERAL LOCATIONS
      INTEGER, INTENT(IN) :: IY       ! CURRENT LATERAL LOCATION INDEX
      REAL, INTENT(IN) :: WT(IYMIN:IYMAX,HMAX,MMAX,0:3) ! ARRAY OF WEIGHTS USED IN 
C                                                       ! NUMERICAL INTEGRATION
      REAL, INTENT(IN) :: Y(IYMIN:IYMAX) ! VECTOR OF y-COORDINATES FOR PMF CALCULATIONS,
                                         ! yd
      INTEGER, INTENT(IN) :: NR        ! CURRENT RUN NUMBER
      INTEGER, INTENT(IN) :: NPATTERNS ! NUMBER OF PATTERNS COMPRISING THIS TACTIC
      TYPE(PATTERNS), INTENT(IN) :: PATTERN(NPATTERNS) ! PATTERN DATA
      INTEGER, INTENT(IN) :: NRUNS    ! NUMBER OF RUNS
      TYPE (RUNS), INTENT(IN) :: RUN(0:NRUNS) ! ARRAY OF RUN DATA
      LOGICAL, INTENT(IN) :: DERIVS   ! WHEN .TRUE., SPECIFIES CALCULATION OF DERIVATIVES
C        NOTE:  AM AND CM CONTAIN DERIVATIVES W/RT $y_i^*$
      REAL, INTENT(OUT) :: AM(0:2,MAXSTATE,MINSTATE:0) ! "A" SUB-MATRIX, W/ DERIVS
      REAL, INTENT(OUT) :: CM(0:2,MAXSTATE,MAXSTATE) ! "C" SUBMATRIX, W/ DERIVS
      CHARACTER*32, INTENT(IN) :: COL_LBL(MINSTATE:1) ! COLUMN LABELS
      CHARACTER*32, INTENT(IN) :: ROW_LBL(0:MAXSTATE) ! ROW LABELS
      LOGICAL, OPTIONAL, INTENT(IN) :: VOCAL ! INVOKES PRINTED OUTPUT WHEN PRESENT AND TRUE.
                                             ! ALSO ENABLES DEBUG POINTS.
C ***
C *** ------------------------------------------------------------------
C ***
C ***  LOCAL VARIABLE DEFINITIONS
C ***
      REAL AO2                         ! HALF OF ACTUATION WIDTH ($A$), yd
      REAL B                           ! CHARACTERISTIC ACTUATION PROBABILITY
      REAL BF                          ! CHARACTERISTIC DAMAGE PROBABILITY
      LOGICAL CALC_DERIVS              ! INVOKES CALCULATION OF DERIVATIVES
      REAL COEFF                       ! COEFFICIENT ($m_\ell y_i^* + b_\ell$)
                                       ! USED IN CALCULATING h(y) & DERIVS.
      CHARACTER*10 CTMP                ! TEMPORARY(SCRATCH) CHARACTER VARIABLE
      CHARACTER*10 CTMP2               ! DITTO
      REAL DELTA_BIGG                  ! FACTOR USED TO CALCULATE h(y) AND h_F(y)
      LOGICAL DERIVS_NE_0              ! INDICATES MATRIX DERIVATIVES CAN BE $\neq$ 0
C        "DERIVS" AND "DERIVS_NE_0" REPRESENT DIFFERENT CONDITIONS.  DERIVS = .TRUE.
C        IMPLIES THAT CALCULATION OF DERIVATIVES IS DESIRED.  DERIVS_NE_0 = .TRUE.
C        IMPLIES THAT A CONDITION UNDER WHICH THE DERIVATIVES ARE $\equiv$ 0 DOES NOT
C        APPLY.  BOTH NEED TO BE .TRUE. FOR DERIVATIVES TO BE CALCULATED.
      INTEGER ELL                      ! LOOP INDEX
      INTEGER ELLMAX                   ! CURRENT VALUE OF MMI3()%ELLMAX
      REAL EM                          ! CURRENT VALUE OF MMI3()%M(ELL)
      REAL F                           ! DANGEROUS FRONT ($F$), yd, FOR I'TH RUN
      LOGICAL FIRST_DNZ /.TRUE./       ! INDICATES WHETHER DERIVS_NE_0 IS EVER FALSE
      REAL FO2                         ! HALF OF DANGEROUS FRONT ($F$), yd
      CHARACTER*80 FORMAT1             ! OUTPUT FORMAT
      CHARACTER*80 FORMAT2             ! DITTO
      REAL, ALLOCATABLE :: G(:,:)      ! GAUSSIAN PDF VALUES, W/ DERIVS
      REAL GLA                         ! GAUSSIAN PDF VALUE AT ZLA
      REAL GLF                         ! GAUSSIAN PDF VALUE AT ZLF
      REAL GUA                         ! GAUSSIAN PDF VALUE AT ZUA
      REAL GUF                         ! GAUSSIAN PDF VALUE AT ZUF
      INTEGER H                        ! CHANNEL SEGMENT NUMBER $h$
      REAL HFOH(0:2)                   ! $h_F(y*)/h(y*)$, W/ DERIVS.
      INTEGER HTMP                     ! CHANNEL SEGMENT INDEX FOR I'TH RUN
      REAL HY                          ! ACTUATION PROBABILITY $h(y)$
      REAL HFY                         ! DAMAGE PROBABILITY
      REAL HNY                         ! NEUTRALIZATION PROBABILITY (HUNTING)
      INTEGER I                        ! RUN INDEX
      INTEGER IACTSUM                  ! SUM OF MAXIMUM NUMBERS OF ACTUATIONS
      INTEGER IORD                     ! ORDER OF PARTIAL DERIVATIVE W/RT $y_i^*$
                                       ! (0 TO 2)
      INTEGER IOS                      ! I/O STATUS VALUE, FOR USE IN ERROR HANDLING
      INTEGER IPAT                     ! PATTERN NUMBER
      INTEGER K                        ! SWEEPER TYPE $k$
      INTEGER KTMP                     ! SWEEPER TYPE USED ON RUN I
      CHARACTER*80 LABEL               ! ALPHANUMERIC LABEL USED IN PRINTED OUTPUTS
      LOGICAL LTMP                     ! TEMPORARY (SCRATCH) LOGICAL VARIABLE
      LOGICAL LTMP2                    ! DITTO
      INTEGER M                        ! MINE TYPE $m$
      REAL MAXABSZ                     ! $|z_\ell|_\max$
      REAL, ALLOCATABLE :: P(:,:)      ! GAUSSIAN PMF VALUES, W/ DERIVS
      REAL PHFY_PYSTAR                 ! FIRST PARTIAL OF $h_F(y_i^*)$ W/RT $y_i^*$
      REAL P2HFY_PYSTAR2               ! SECOND PARTIAL 
      REAL PHNY_PYSTAR                 ! FIRST PARTIAL OF $h_N(y_i^*)$ W/RT $y_i^*$
      REAL P2HNY_PYSTAR2               ! SECOND PARTIAL 
      REAL PHY_PYSTAR                  ! FIRST PARTIAL OF $h(y_i^*)$ W/RT $y_i^*$
      REAL P2HY_PYSTAR2                ! SECOND PARTIAL
      REAL PUA                         ! PMF VALUE AT ZUA
      REAL PLA                         ! PMF VALUE AT ZLA
      REAL U                           ! NUMERATOR IN STEADY-STATE CALCULATIONS
      INTEGER Q                        ! SWEEPER DAMAGE INDEX (1 = DAMAGED, 2 = UNDAMAGED)
      REAL RTMP                        ! TEMPORARY (SCRATCH) REAL VARIABLE
      REAL RTMP1                       ! DITTO
      REAL RTMP2                       ! DITTO
      INTEGER S                        ! STATE INDEX $s$
      INTEGER S2                       ! DITTO
      INTEGER SCAP                     ! SCAP = SMAX + 1
      INTEGER SCOL                     ! STATE INDEX $s$; COLUMN NUMBER IN "A" MATRIX
      INTEGER SCOL2                    ! STATE INDEX; COLUMN NUMBER WHEN Q=2 (NO DAMAGE TO SWEEPER)
      INTEGER SCOL3                    ! STATE INDEX; COLUMN NUMBER WHEN Q=3 (HUNTER NEUTRALIZES MINE)
      REAL SIG                         ! STANDARD DEVIATION OF NAVIGATION ERROR (SDNE) $\sigma$, yd
      REAL SIGSQ                       ! $\sigma^2$, yd^2
      LOGICAL SINGLE                   ! .TRUE. IF THIS RUN IS A SINGLE RUN THROUGH THE CHANNEL;
                                       ! .FALSE., IF IT REPRESENTS AN INFINITE NUMBER OF
                                       ! RUNS ON THIS TRACK (I.E., "STEADY STATE")
      INTEGER SROW                     ! ROW INDEX
      INTEGER SUBPAT                   ! SUB-PATTERN NUMBER FOR THIS RUN
      CHARACTER*140 TWEET              ! ALPHANUMERIC MESSAGE FOR TRANSFER TO LOG FILE
      REAL YSTAR                       ! DISTANCE OF Y-COORDINATE FROM INTENDED TRACK, yd ($y_i^*$)
      REAL, ALLOCATABLE :: Z_(:)       ! STANDARIZED u_\ELL VALUES
      REAL ZL                          ! LEFT LIMIT OF INTEGRATION
      REAL ZLA                         ! STANDARDIZED NORMAL DEVIATE AT $y_i^* - A/2$
      REAL ZR                          ! RIGHT LIMIT OF INTEGRATION
      REAL ZUA                         ! STANDARDIZED NORMAL DEVIATE AT $y_i^* + A/2$
      REAL ZLF                         ! STANDARDIZED NORMAL DEVIATE AT $y_i^* - F/2$
      REAL ZUF                         ! STANDARDIZED NORMAL DEVIATE AT $y_i^* + F/2$
C ***
C *** ------------------------------------------------------------------
C ***
C ***
C *** EXECUTABLE CODE
C ***
C      WRITE (*,'("LOADER_N: IY =",I10)') IY
      IF (IY == 99) THEN
         CONTINUE ! DEBUG POINT
      ENDIF ! IY = 99
      SINGLE = RUN(NR)%SINGLE
      IF (NR == NRUNS) THEN
         CONTINUE ! DEBUG POINT
      ENDIF ! NR = NRUNS
C        INITIALIZE SUB-MATRICES
      AM = 0.
      CM = 0.
      IF (SINGLE) THEN
         DO S = 1, MAXSTATE
            CM(0,S,S) = 1. ! CM(1,*,*) AND CM(2,*,*) $\equiv$ 0.
         ENDDO ! S
      ELSE ! "STEADY-STATE", SO CM $\equiv$ 0
         CONTINUE
      ENDIF ! SINGLE
      IPAT = RUN(NR)%IPAT
      H = PATTERN(IPAT)%H
      K = PATTERN(IPAT)%K
      SUBPAT = RUN(NR)%SUBPAT
      IF (SUBPAT == 0) THEN
         CONTINUE ! DEBUG POINT
      ENDIF ! SUBPAT=0
      DO M = 1, MMAX
C           FOLLOWING LOGIC DETERMINES WHETHER THE SCD IS UNIFORM, PRECLUDING
C           NON-ZERO DERIVATIVES FOR SOME VALUES OF "NR".  SEE NOTES OF 11/24/15.
C           THE LOGIC BELOW IS SIMILAR TO THAT IN "IS_BREAKTHROUGH".  HOWEVER, THE
C           PENALTY FUNCTION OPTIMIZATION LOGIC ALLOWS THE SWEEPER AND TRANSITOR
C           NAVIGATION LIMITS TO BE VIOLATED.  ACCORDINGLY, THE LOGIC BELOW USES
C           RUN LOCATIONS (RUN(I)%YI) RATHER THAN NAVIGATION LIMITS.  SEE NOTES OF
C           7/9/18.
         LTMP = .FALSE. ! DEFAULT IS TO CALCULATE DERIVATIVES
         IACTSUM = 0
         I_LOOP: DO I = 1, NR
            HTMP = PATTERN(RUN(I)%IPAT)%H
            KTMP = PATTERN(RUN(I)%IPAT)%K
            IF (HTMP == H) THEN ! THIS RUN IS IN THE SAME SEGMENT AS RUN "NR"
C                 THE FOLLOWING TESTS WHETHER THE SCD IS UNIFORM UP TO S/C %SMAX
               LTMP = (MINE(H,M)%BETA == 1.) .AND. (MINE(H,M)%IDIST <= 2)
               LTMP = LTMP .AND. MINE(H,M)%RHO > 0. ! ADDS REQ'T THAT DENSITY BE > 0
               IF (LTMP) THEN
C                    THE FOLLOWING STMT, AND OTHERS IN THIS CODE BLOCK, IS
C                    HARD-WIRED FOR SUBPAT = 1 AND IACT = 1
                  F = MMI2(H,KTMP,M,1,1)%F
C                    TODO: CONSIDER ADDING "& MMI2(H,KTMP,M,1,1)%BF > 0."
                  LTMP = LTMP .AND. F > 0.
                  RTMP = 0.5*F + ZSSC*MSC(H,KTMP)%SIG
                  LTMP = LTMP .AND. RUN(I)%YI - RTMP >= MINE(H,M)%YFL
     &                        .AND. RUN(I)%YI + RTMP <= MINE(H,M)%YFR
C                    THE FOLLOWING DEPENDS ON THE "BUSINESS RULE" THAT %SINGLE IS
C                    ALWAYS TRUE WHEN %SWPR IS, I.E., TRUE SWEEPERS ONLY MAKE SINGLE RUNS
                  IF (MSC(H,KTMP)%SWPR .OR. MSC(H,KTMP)%SINGLE) THEN
                     IACTSUM = IACTSUM + MMI1(H,KTMP,M)%IACTMAX
                  ELSE ! NOT MSC(H,KTMP)%SWPR & MSC(H,KTMP)%SINGLE, I.E.,
                       ! THE "INFINITE HOMOGENEOUS TRANSITOR CASE"
                     IACTSUM = HUGE(IACTSUM)
                  ENDIF ! %SWPR OR %SINGLE
               ENDIF ! LTMP
            ENDIF ! HTMP = H
            IF (IACTSUM .GT. MINE(H,M)%SMAX) THEN ! TODO: DOUBLE-CHECK THIS LOGIC
               LTMP = .FALSE.
               EXIT I_LOOP
            ENDIF ! IACTSUM > %SMAX
         ENDDO I_LOOP
         DERIVS_NE_0 = .NOT. LTMP
         CALC_DERIVS = DERIVS .AND. DERIVS_NE_0
         IF (FIRST_DNZ .AND. .NOT. DERIVS_NE_0) THEN
            FIRST_DNZ = .FALSE.
            TWEET = "LOADER_N: DERIVS_NE_0 IS FALSE"
            CALL TWITTER (TWEET)
         ENDIF ! FIRST_DNZ & ~DERIVS_NE_0
         SCAP = MINE(H,M)%SCAP
C           NOTE USE OF ZERO AS LAST INDEX OF "WT"; ENSURES THAT ALL MINES CAN AFFECT EMS,
C           E(S)C, AND ETC
C           USE OF 1 AS LAST INDEXES OF "MMI3" HARD-WIRES CODE FOR IACT = 1 & INDEX = 1
C           (A)
         IF (WT(IY,H,M,0) .GT. 0. .AND. MMI3(H,K,M,SUBPAT,1,1)%WIDTH > 0.) THEN
C              PMF VALUES FOR THIS y-COORDINATE RECEIVE NON-ZERO WEIGHTS IN NUMERICAL
C              INTEGRATION
            YSTAR = Y(IY) - RUN(NR)%YI
            SIG = MSC(H,K)%SIG
            SIGSQ = SIG*SIG
C              CALL "PIECEWISE" WITH INDEX = 1 (ACTUATION PROBABILITIES)
            LTMP = YSTAR == -1200.
            CALL PIECEWISE (HMAX, KMAX, MMAX, MMI3, H, K, M, SUBPAT, 1,
     &       CALC_DERIVS, YSTAR, SIG, HY, PHY_PYSTAR, P2HY_PYSTAR2, VOCAL = LTMP)
C            ELLMAX = MMI3(H,K,M,SUBPAT,1,1)%ELLMAX
C            ALLOCATE ( P(ELLMAX,0:2), G(ELLMAX,0:2), Z_(ELLMAX) )
C            DO ELL = 1, ELLMAX
C               Z_(ELL) = Z(YSTAR-MMI3(H,K,M,SUBPAT,1,1)%U(ELL),0.,SIG)
C               P(ELL,0) = ANORDF(Z_(ELL))
C               G(ELL,0) = ANORPR(Z_(ELL))
C            ENDDO ! ELL
C            HY = 0.
C            PHY_PYSTAR = 0.
C            P2HY_PYSTAR2 = 0.
C            DO ELL = 1, ELLMAX-1
C               COEFF = MMI3(H,K,M,SUBPAT,1,1)%M(ELL)*YSTAR + MMI3(H,K,M,SUBPAT,1,1)%B(ELL)
C               MAXABSZ = MAX(ABS(Z_(ELL)), ABS(Z_(ELL+1)))
C               IF (MAXABSZ <= 10.) THEN ! FORMERLY "<= ZMAX"
C                  DELTA_BIGG = P(ELL,0) - P(ELL+1,0)
C               ELSE ! MAXABSZ > ZMAX - USE "ALTERNATIVE FORM" (MLM 5/13/20)
C                    ! NB: EXACT FORM ALWAYS USED FOR DERIVATIVES.
C                  DELTA_BIGG = 0.5 * (G(ELL,0)+G(ELL+1,0)) * (Z_(ELL)-Z_(ELL+1))
C               ENDIF ! MAXABSZ <= 10
C               RTMP = COEFF * DELTA_BIGG
C               IF (MMI3(H,K,M,SUBPAT,1,1)%M(ELL) /= 0.) THEN
C                  RTMP = RTMP + MMI3(H,K,M,SUBPAT,1,1)%M(ELL)*SIG*(G(ELL,0)-G(ELL+1,0))
C               ENDIF ! %M(ELL) /= 0
C               HY = HY + RTMP
C               IF (CALC_DERIVS) THEN
C                  ZR = Z_(ELL)
C                  ZL = Z_(ELL+1)
C                  RTMP1 = MMI3(H,K,M,SUBPAT,1,1)%M(ELL) * DELTA_BIGG
C     &                    + COEFF * (G(ELL,0) - G(ELL+1,0)) / SIG
C                  RTMP2 = 2.*MMI3(H,K,M,SUBPAT,1,1)%M(ELL) * (G(ELL,0) - G(ELL+1,0)) / SIG
C     &                    - COEFF * (ZR*G(ELL,0) - ZL*G(ELL+1,0)) / SIGSQ
C                  IF (MMI3(H,K,M,SUBPAT,1,1)%M(ELL) /= 0.) THEN
C                     RTMP1 = RTMP1 - MMI3(H,K,M,SUBPAT,1,1)%M(ELL) *(ZR*G(ELL,0) - ZL*G(ELL+1,0))
C                     RTMP2 = RTMP2 - MMI3(H,K,M,SUBPAT,1,1)%M(ELL) * ((1.-ZR)*G(ELL,0) - (1.-ZL)*G(ELL+1,0)) / SIG
C                  ENDIF ! %M /= 0
C                  PHY_PYSTAR = PHY_PYSTAR + RTMP1
C                  P2HY_PYSTAR2 = P2HY_PYSTAR2 + RTMP2
C              ENDIF ! CALC_DERIVS
C            ENDDO ! ELL
C            DEALLOCATE ( P, G, Z_ )

CC              TODO:  MODIFY THE FOLLOWING TO ALLOW "HY_MIN = 0" IF ONLY "ZEROTH-
CC              DERIVATIVE" DATA ARE DESIRED.  THIS WOULD ALLOW CALCULATION AND OUTPUT OF
CC              DETONATION PROBABILITIES THAT WOULD OTHERWISE NOT BE CALCULATED TO PRECLUDE
CC              UNDERFLOW OF $h(y)^2$ AND $h(y)^3$.
            IF (ISNAN(HY)) THEN
               TWEET = "LOADER_N: WARNING - HY IS A NaN!"
               CALL TWITTER (TWEET, PAWS = .TRUE.)
            ENDIF ! ISNAN(HY)
            IF (HY < 0. .OR. HY > 1.) THEN ! ERROR TRAP
               TWEET = "LOADER_N: WARNING - HY < 0 OR > 1"
C               CALL TWITTER (TWEET, PAWS = .FALSE.)
            ENDIF ! HY < 0 OR > 1
            IF (HY .GT. HY_MIN) THEN ! REVISE VALUES IN "A" AND "C"
C               GUA = ANORPR(ZUA)
C               GLA = ANORPR(ZLA)
C               PHY_PYSTAR = B * (GUA - GLA) / SIG ! $\partial h(y_i^*) / \partial y_i^*$
C               P2HY_PYSTAR2 = -B * (ZUA*GUA - ZLA*GLA) / SIGSQ ! SECOND PARTIAL
               SROW = SCAPSUM(H,M-1) + 1
C               IF (M .EQ. 2) THEN
C                  CONTINUE ! DEBUG POINT
C               ENDIF ! M = 2
               IF (SINGLE) THEN
                  CM(0,SROW,SROW) = 1. - HY ! $q_i(y)$
                  CM(0,SROW+1,SROW) = HY
C                    THE FOLLOWING STATEMENTS USE THE "CHAIN RULE" TO CALCULATE FIRST AND
C                    SECOND PARTIALS W/RT $D$
                  IF (CALC_DERIVS) THEN
                     CM(1,SROW,SROW) = -PHY_PYSTAR
                     CM(1,SROW+1,SROW) = PHY_PYSTAR
                     CM(2,SROW,SROW) = -P2HY_PYSTAR2
                     CM(2,SROW+1,SROW) = P2HY_PYSTAR2
                  ENDIF ! CALC_DERIVS
               ELSE ! STEADY-STATE
                  CONTINUE ! NO ACTION REQUIRED
               ENDIF ! SINGLE       

C                 INDEX=2 INDICATES DAMAGE DATA BEING USED
               CALL PIECEWISE (HMAX, KMAX, MMAX, MMI3, H, K, M, SUBPAT, 2,
     &          CALC_DERIVS, YSTAR, SIG, HFY, PHFY_PYSTAR, P2HFY_PYSTAR2, VOCAL = YSTAR == -1200.)
C               ELLMAX = MMI3(H,K,M,SUBPAT,1,2)%ELLMAX
C               ALLOCATE ( P(ELLMAX,0:2), G(ELLMAX,0:2), Z_(ELLMAX) )
C               P = 0.
C               G = 0.
C               Z_ = 0.
C               DO ELL = 1, ELLMAX
C                  Z_(ELL) = Z(YSTAR-MMI3(H,K,M,SUBPAT,1,2)%U(ELL), 0., SIG)
C                  P(ELL,0) = ANORDF(Z_(ELL))
C                  G(ELL,0) = ANORPR(Z_(ELL))
C               ENDDO ! ELL
C               HFY = 0.
C               PHFY_PYSTAR = 0.
C               P2HFY_PYSTAR2 = 0.
C               DO ELL = 1, ELLMAX-1
C                  EM = MMI3(H,K,M,SUBPAT,1,2)%M(ELL)
C                  MAXABSZ = MAX(ABS(Z_(ELL)), ABS(Z_(ELL+1)))
C                  IF (MAXABSZ <= 10.) THEN
C                     DELTA_BIGG = P(ELL,0) - P(ELL+1,0)
C                  ELSE ! MAXABSZ > 10
C                     DELTA_BIGG = 0.5 * (G(ELL,0)+G(ELL+1,0)) * (Z_(ELL)-Z_(ELL+1))
C                  ENDIF ! MAXABSZ <= 10
C                  COEFF = EM*YSTAR + MMI3(H,K,M,SUBPAT,1,2)%B(ELL)
C                  RTMP = COEFF * DELTA_BIGG
C                  RTMP1 = 0.
C                  RTMP2 = 0.
C                  IF (EM /= 0.) THEN
C                     RTMP = RTMP + EM*SIG*(G(ELL,0)-G(ELL+1,0))
C                  ENDIF ! EM /= 0
C                  HFY = HFY + RTMP
C                  IF (CALC_DERIVS) THEN
C                     ZR = Z_(ELL)
C                     ZL = Z_(ELL+1)
C                     RTMP1 = EM * DELTA_BIGG
C     &                       + COEFF * (G(ELL,0) - G(ELL+1,0))/SIG
C                     RTMP2 = 2. * EM * (G(ELL,0) - G(ELL+1,0)) / SIG
C     &                       - COEFF * (ZR*G(ELL,0) - ZL*G(ELL+1,0)) / SIGSQ
C                     IF (EM /= 0.) THEN
C                        RTMP1 = RTMP1 - EM*(ZR*G(ELL,0)-ZL*G(ELL+1,0))
C                        RTMP2 = RTMP2 - EM * ((1.-ZR)*G(ELL,0) - (1.-ZL)*G(ELL+1,0)) / SIG
C                     ENDIF ! EM /= 0
C                     PHFY_PYSTAR = PHFY_PYSTAR + RTMP1
C                     P2HFY_PYSTAR2 = P2HFY_PYSTAR2 + RTMP2
C                  ENDIF ! CALC_DERIVS
C               ENDDO ! ELL
C               DEALLOCATE ( P, G, Z_ )
C               IF (M == 1) WRITE (IOUT,'(E16.6,2(",",E16.6))') YSTAR, HY, HFY ! TEMPORARY CODE
               IF (ISNAN(HFY)) THEN
                  TWEET = "LOADER_N: WARNING - HFY IS A NaN!"
                  CALL TWITTER (TWEET, PAWS = .TRUE.)
               ENDIF ! ISNAN(HFY)
               IF (HFY < 0. .OR. HFY > 1.) THEN ! ERROR TRAP
                  TWEET = "LOADER_N: WARNING - HFY < 0 OR > 1"
                  CALL TWITTER (TWEET, PAWS = .FALSE.)
               ENDIF ! HFY < 0 OR > 1
C                 A DETONATION THAT DAMAGES THE SWEEPER CORRESPONDS TO Q = 1
               CALL ENCODEC (H, K, M, 1, KMQ_MAX, MQ_MAX, SCOL)
               CALL ENCODEC (H, K, M, 2, KMQ_MAX, MQ_MAX, SCOL2)
C               IF (SCOL .GT. 0) THEN
C                  CONTINUE ! DEBUG POINT
C               ENDIF ! SCOL > 0
               IF (SINGLE) THEN
                  AM(0,SROW,SCOL) = HFY
C                    SEE COMMENT ABOVE RE "CHAIN RULE"  
                  AM(0,SROW,SCOL2) = HY - HFY ! $s_i(y) = h(y_i^*) - h_F(y_i^*)$
                  IF (CALC_DERIVS) THEN
                     AM(1,SROW,SCOL) = PHFY_PYSTAR ! FIRST PARTIAL
                     AM(2,SROW,SCOL) = P2HFY_PYSTAR2 ! SECOND PARTIAL
                     AM(1,SROW,SCOL2) = PHY_PYSTAR - PHFY_PYSTAR
                     AM(2,SROW,SCOL2) = P2HY_PYSTAR2 - P2HFY_PYSTAR2
                  ENDIF ! CALC_DERIVS
               ELSE ! NOT SINGLE; $\equiv$ STEADY-STATE
C                  IF (PRESENT(VOCAL) .AND. VOCAL) THEN
C                     CONTINUE ! DEBUG POINT
C                  ENDIF ! PRESENT (VOCAL) & VOCAL
!                  IF (HFY .GT. HY .OR. HFY < EPS) THEN ! ROUND-OFF ERROR IS LIKELY RESPONSIBLE; SET HFY
                                        ! TO ZERO
                  IF (HFY .GT. HY) THEN
                     CALL CONVERT (NR, CTMP2)      
                     IF (IY .GE. 0) THEN
                        CALL CONVERT (IY, CTMP)
                        TWEET = "LOADER_N: HFY SET TO ZERO FOR IY = "//TRIM(ADJUSTL(CTMP))//
     &                     " AND RUN NO. "//TRIM(ADJUSTL(CTMP2))
                     ELSE ! IY < 0
                        CALL CONVERT (-IY, CTMP)
                        TWEET = "LOADER_N: HFY SET TO ZERO FOR IY = -"//TRIM(ADJUSTL(CTMP))//
     &                     " AND RUN NO. "//TRIM(ADJUSTL(CTMP2))
                     ENDIF ! IY >= 0
                     CALL TWITTER (TWEET, PAWS = .TRUE.)
                     HFY = 0.
                     PHFY_PYSTAR = 0.
                     P2HFY_PYSTAR2 = 0.
                  ENDIF ! HFY > HY
                  AM(0,SROW,SCOL) = HFY/HY
                  AM(0,SROW,SCOL2) = (HY-HFY)/HY
                  IF (CALC_DERIVS) THEN
                     RTMP = HY*HY
                     IF (RTMP == 0.) THEN
                        TWEET = "LOADER_N: HY**2 UNDERFLOW TO ZERO"
                        CALL TWITTER (TWEET, PAWS = .TRUE.)
                     ENDIF ! RTMP = 0
                     AM(1,SROW,SCOL) = (PHFY_PYSTAR*HY - PHY_PYSTAR*HFY) / RTMP
                     RTMP = RTMP*HY
                     IF (RTMP == 0.) THEN
                        TWEET = "LOADER_N: HY**3 UNDERFLOW TO ZERO"
                        CALL TWITTER (TWEET, PAWS = .TRUE.)
                     ENDIF ! RTMP = 0
                     RTMP2 = P2HFY_PYSTAR2*HY*HY - 2.*PHFY_PYSTAR*HY*PHY_PYSTAR
     &                       + 2.*HFY*(PHY_PYSTAR**2) - HFY*HY*P2HY_PYSTAR2
                     AM(2,SROW,SCOL) = RTMP2 / RTMP
                     AM(1,SROW,SCOL2) = -AM(1,SROW,SCOL)
                     AM(2,SROW,SCOL2) = -AM(2,SROW,SCOL)
                  ENDIF ! CALC_DERIVS
C                    FILL IN REMAINING ROWS OF "A" FOR THIS SEGMENT AND MINE TYPE
                  DO S = SROW+1, SCAPSUM(H,M)
                     DO S2 = SCOL2, SCOL ! SCOL2 < SCOL
                        DO IORD = 0, 2
                           AM(IORD,S,S2) = AM(IORD,SROW,S2)
                        ENDDO ! IORD
                     ENDDO ! S2
                  ENDDO ! S
               ENDIF ! SINGLE
            ELSE ! HY = 0
               SROW = SCAPSUM(H,M-1) + 1
C                 THE FOLLOWING ARE HARD-WIRED FOR Q = 1 & Q = 2, RESP'Y.
               CALL ENCODEC (H,K,M,1,KMQ_MAX,MQ_MAX,SCOL)
               CALL ENCODEC (H,K,M,2,KMQ_MAX,MQ_MAX,SCOL2)
               IF (SINGLE) THEN
                  CONTINUE ! NO ACTION REQUIRED
               ELSE ! ~SINGLE; INFINITE HOMOGENEOUS TRANSITOR CASE
C                    THE FOLLOWING IS HARD-WIRED FOR IACT = 1
                  IF (MMI2(H,K,M,RUN(NR)%SUBPAT,1)%F < MMI2(H,K,M,RUN(NR)%SUBPAT,1)%A) THEN
C                       MINES TEND TO BE SWEPT WITHOUT DAMAGE TO SWEEPER
                     AM(0,SROW,SCOL) = 0.
                     AM(0,SROW,SCOL2) = 1.
                  ELSE ! %F = %A
                     RTMP = MMI2(H,K,M,RUN(NR)%SUBPAT,1)%BF / MMI2(H,K,M,RUN(NR)%SUBPAT,1)%B
                     AM(0,SROW,SCOL) = 1. - RTMP
                     AM(0,SROW,SCOL2) = RTMP
                  ENDIF ! %F < %A
                  AM(1,SROW,SCOL) = 0.
                  AM(1,SROW,SCOL2) = 0.
                  AM(2,SROW,SCOL) = 0.
                  AM(2,SROW,SCOL2) = 0.
                  DO S = SROW+1, SCAPSUM(H,M)
                     DO S2 = SCOL2, SCOL ! SCOL2 < SCOL
                        DO IORD = 0, 2
                           AM(IORD,S,S2) = AM(IORD,SROW,S2)
                        ENDDO ! IORD
                     ENDDO ! S2
                  ENDDO ! S
               ENDIF ! SINGLE
            ENDIF ! HY > HY_MIN
            IF (MSC(H,K)%HUNTR) THEN
C                 CALCULATE NEUTRALIZATION PROBABILITIES (h_N(y*)); INDEX = 3
C               LTMP = IY == 20
               LTMP = .FALSE.
               CALL PIECEWISE (HMAX, KMAX, MMAX, MMI3, H, K, M, SUBPAT, 3,
     &          CALC_DERIVS, YSTAR, SIG, HNY, PHNY_PYSTAR, P2HNY_PYSTAR2,
     &          VOCAL = YSTAR == -1200.) ! FORMERLY LTMP
C               IF (IY == 20) WRITE (IOUT,'("LOADER_N: M;NR;YSTAR;HY;HFY;HNY =,",I10,",",I10,4(",",E16.6))')
C     &          M, NR, YSTAR, HY, HFY, HNY
               CONTINUE ! DEBUG POINT
               IF (HNY < 0.) THEN
                  HNY = 0.
                  PHNY_PYSTAR = 0.
                  P2HNY_PYSTAR2 = 0.
               ENDIF ! HNY < 0.
C                 REVISE VALUES ALREADY ENTERED IN "A" & "C" MATRICES
               SROW = SCAPSUM(H,M-1) + 1 ! PROBABLY UNNECESSARY
               CM(0,SROW,SROW) = CM(0,SROW,SROW) * (1.-HNY) ! REVISED "q"
               CM(0,SROW+1,SROW) = CM(0,SROW+1,SROW) * (1.-HNY) ! REVISED "p"
               IF (CALC_DERIVS) THEN
                  TWEET = "LOADER_N: DERIVATIVES NOT LOADED IN A & C MATRICES"
                  CONTINUE ! DEBUG POINT
C                  CALL TWITTER (TWEET, PAWS = .FALSE.)
              ENDIF ! CALC_DERIVS
              CALL ENCODEC (H, K, M, 1, KMQ_MAX, MQ_MAX, SCOL) ! PROB. UNNECESSARY
              CALL ENCODEC (H, K, M, 2, KMQ_MAX, MQ_MAX, SCOL2) ! DITTO
              CALL ENCODEC (H, K, M, 3, KMQ_MAX, MQ_MAX, SCOL3) ! NECESSARY (NOT PREVIOUSLY DEFINED)
C             AM(0,SROW,SCOL) = HFY ! "r"; NO CHANGE
C             AM(0,SROW,SCOL2) = HY - HFY ! "s"; NO CHANGE
              AM(0,SROW,SCOL3) = (1.-HY) * HNY ! "u" - Pr(POISED MINE IS NEUTRALIZED) - SEE NOTES 5/25/20
              DO S = 2, MINE(H,M)%SCAP
                 SROW = SCAPSUM(H,M-1) + S
                 AM(0,SROW,SCOL3) = HNY ! "t", THE PROBABILITY THAT AN UNPOISED MINE IS NEUTRALIZED
              ENDDO ! S
            ENDIF ! MSC(H,K)%HUNTR
C              COPY VALUES TO REMAINDER OF "C" MATRIX
            DO S = 2, MINE(H,M)%SCAP
               SCOL = SCAPSUM(H,M-1) + S
               DO SROW = SCOL, SCAPSUM(H,M)
                  DO IORD = 0, 2
                     CM(IORD,SROW,SCOL) = CM(IORD,SROW-1,SCOL-1)
                  ENDDO ! IORD
               ENDDO ! SROW
            ENDDO ! S
         ENDIF ! WT(IY,H,M,0) > 0 & %A > 0
      ENDDO ! M
C        BEGIN DEBUGGING PRINTOUTS
      IF (VOCAL) THEN
         CALL CONVERT (NR, CTMP2)      
         IF (IY .GE. 0) THEN
            CALL CONVERT (IY, CTMP)
            LABEL = "LOADER_N: A & C MATRICES & DERIVS. W/RT y_i^* FOR IY = "//TRIM(ADJUSTL(CTMP))//
     &         " AND RUN NO. "//TRIM(ADJUSTL(CTMP2))
         ELSE ! IY < 0
            CALL CONVERT (-IY, CTMP)
            LABEL = "LOADER_N: A & C MATRICES & DERIVS. W/RT y_i^* FOR IY = -"//TRIM(ADJUSTL(CTMP))//
     &         " AND RUN NO. "//TRIM(ADJUSTL(CTMP2))
         ENDIF ! IY >= 0
         WRITE (IOUT,'(A80)') LABEL  
         CALL PRN_ACM (MINSTATE, MAXSTATE, AM, CM, COL_LBL, ROW_LBL)
      ENDIF ! VOCAL
C        END DEBUGGING PRINTOUTS

      RETURN
      END SUBROUTINE LOADER_N
      
      SUBROUTINE PRN_ACM (MINSTATE, MAXSTATE, AM, CM, COL_LBL, ROW_LBL)
      USE FILES
      USE UTILS ! CONTAINS "CONVERT"
      USE CODE ! CONTAINS "GEN_LBLS"
      IMPLICIT NONE
C *** 
C *** PURPOSE: PRINTS "A" AND "C" SUB-MATRICES
C ***
C *** METHODOLOGY: N/A
C ***
C *** MODIFICATION HISTORY: VERSION 1.0 - TBD (INITIAL RELEASE)
C ***
C *** ------------------------------------------------------------------
C ***
C *** ARGUMENTS:
C ***
      INTEGER, INTENT(IN) :: MINSTATE ! SMALLEST STATE INDEX
      INTEGER, INTENT(IN) :: MAXSTATE ! LARGEST STATE INDEX
      REAL, INTENT(IN) :: AM(0:2,MAXSTATE,MINSTATE:0) ! "A" SUB-MATRIX
      REAL, INTENT(IN) :: CM(0:2,MAXSTATE,MAXSTATE) ! "C" SUB-MATRIX
      CHARACTER*43, INTENT(IN) :: COL_LBL(MINSTATE:1) ! COLUMN LABELS
      CHARACTER*32, INTENT(IN) :: ROW_LBL(0:MAXSTATE) ! ROW LABELS
C ***
C *** ------------------------------------------------------------------
C ***
C ***  LOCAL VARIABLE DEFINITIONS
C ***
      CHARACTER*10 CTMP               ! TEMPORARY (SCRATCH) CHARACTER VARIABLE
      CHARACTER*10 CTMP2              ! DITTO
      CHARACTER*80 FORMAT1            ! RUN-TIME OUTPUT FORMAT
      CHARACTER*80 FORMAT2            ! DITTO
      CHARACTER*80 FORMAT3            ! DITTO
      INTEGER IORD                    ! ORDER OF DERIVATIVE (0 TO 2)
      INTEGER IOS                     ! I/O STATUS VARIABLE; USED FOR ERROR HANDLING
      CHARACTER*80 LABEL              ! OUTPUT LABEL
      INTEGER S                       ! STATE VARIABLE
      INTEGER SROW                    ! ROW NUMBER
C ***
C *** ------------------------------------------------------------------
C ***
C ***
C *** EXECUTABLE CODE
C ***
      CALL CONVERT (1-MINSTATE, CTMP)
      CALL CONVERT (MAXSTATE, CTMP2)
      FORMAT1 = '(10X,'//TRIM(ADJUSTL(CTMP))//'(",",A43),'//TRIM(ADJUSTL(CTMP2))//'(",",A32))'
      FORMAT3 = '(A32,",",A43,'//TRIM(ADJUSTL(CTMP))//'(","),A32)'
      CALL CONVERT (MAXSTATE-MINSTATE+1, CTMP2)
      FORMAT2 = '(A32,'//TRIM(ADJUSTL(CTMP2))//'(",",E16.6))'
      DO IORD = 0, 2
         CALL CONVERT (IORD, CTMP)
         LABEL = "IORD = "//TRIM(ADJUSTL(CTMP))
         WRITE (IOUT,'(/A80)') LABEL
         WRITE (IOUT,FORMAT3) ROW_LBL(0), COL_LBL(1), ROW_LBL(0)
         WRITE (IOUT,FORMAT1,IOSTAT=IOS) (COL_LBL(S), S = MINSTATE, 0), (ROW_LBL(S), S = 1, MAXSTATE)
         DO SROW = 1, MAXSTATE
            WRITE (IOUT,FORMAT2,IOSTAT=IOS) ROW_LBL(SROW),
     &         (AM(IORD,SROW,S), S = MINSTATE, 0), (CM(IORD,SROW,S), S = 1, MAXSTATE)   
         ENDDO ! SROW
      ENDDO ! IORD

      RETURN      
      END SUBROUTINE PRN_ACM
      
      SUBROUTINE PIECEWISE (HMAX, KMAX, MMAX, MMI3, H, K, M, SUBPAT, INDEX, 
     & CALC_DERIVS, YSTAR, SIG, HY, PHY_PYSTAR, P2HY_PYSTAR2, VOCAL)
      USE DERIVED_TYPES
      USE EM_IMSL
      USE FILES
      USE GAUSS
      USE UTILS
      IMPLICIT NONE
C *** 
C *** PURPOSE: CALCULATES HY/HFY/HNY, W/ DERIVS W/RT y*, USING PIECEWISE
C *** LINEAR DATA
C ***
C *** METHODOLOGY: ALGEBRA, PARTIAL DIFFERENTIATION
C ***
C *** PENDING MODIFICATIONS ("TODO" ITEMS):
C ***
C *** MODIFICATION HISTORY: VERSION 1.0 - TBD (INITIAL RELEASE)
C ***
C *** ------------------------------------------------------------------
C ***
C *** ARGUMENTS:
C ***                                 !
      INTEGER, INTENT(IN) :: HMAX     ! NUMBER OF CHANNEL SEGMENTS
      INTEGER, INTENT(IN) :: KMAX     ! NUMBER OF SWEEPER/HUNTER TYPES
      INTEGER, INTENT(IN) :: MMAX     ! NUMBER OF MINE TYPES
      TYPE(MMI3S), INTENT(IN) :: MMI3(HMAX,KMAX,MMAX,0:1,1,3) ! PIECEWISE DATA
      INTEGER, INTENT(IN) :: H        ! SEGMENT NUMBER
      INTEGER, INTENT(IN) :: K        ! SWEEPER/HUNTER INDEX
      INTEGER, INTENT(IN) :: M        ! MINE INDEX
      INTEGER, INTENT(IN) :: SUBPAT   ! SUB-PATTERN INDEX
      INTEGER, INTENT(IN) :: INDEX    ! DATA TYPE INDEX (ACTUATION/DAMAGE/
                                      ! NEUTRALIZATION)
      LOGICAL, INTENT(IN) :: CALC_DERIVS ! SPECIFIES CALCULATION OF DERIVATIVES
      REAL, INTENT(IN) :: YSTAR       ! RELATIVE LATERAL DISTANCE (y-y_i)
      REAL, INTENT(IN) :: SIG         ! STANDARD DEVIATION OF NAVIGATION ERROR, yd
      REAL, INTENT(OUT) :: HY         ! ACTUATION/DAMAGE/NEUTRALIZATION
                                      ! PROBABILITY
      REAL, INTENT(OUT) :: PHY_PYSTAR ! FIRST PARTIAL W/RT y*
      REAL, INTENT(OUT) :: P2HY_PYSTAR2 ! SECOND PARTIAL W/RT y*
      LOGICAL, OPTIONAL, INTENT(IN) :: VOCAL ! INVOKES PRINTED OUTPUT
C ***
C *** ------------------------------------------------------------------
C ***
C ***  LOCAL VARIABLE DEFINITIONS
C ***                                 !
      REAL COEFF                      ! COEFFICIENT (my* + b)
      CHARACTER*10 CTMP               ! TEMPORARY (SCRATCH) CHARACTER VARIABLE
      REAL DELTA_BIGG                 ! DIFFERENCE IN GAUSSIAN PMF
      INTEGER ELL                     ! DATA SEGMENT INDEX
      INTEGER ELLMAX                  ! MAXIMUM VALUE OF ELL
      REAL EM                         ! DATA SEGMENT SLOPE
      CHARACTER*140 FORMAT1           ! RUN-TIME OUTPUT FORMAT
      CHARACTER*140 FORMAT2           ! DITTO
      REAL, ALLOCATABLE :: G(:,:)     ! GAUSSIAN PDF VALUES
      INTEGER IACT                    ! NUMBER OF ACTUATIONS
      INTEGER ITMP                    ! TEMPORARY (SCRATCH) INTEGER VARIABLE
      LOGICAL LTMP                    ! TEMPORARY (SCRATCH) LOGICAL VARIABLE
      REAL MAXZ                       ! MAXIMUM z_\ell
      REAL, ALLOCATABLE :: P(:,:)     ! GAUSSIAN PMF VALUES
      REAL RTMP                       ! TEMPORARY (SCRATCH) REAL VARIABLE
      REAL RTMP1                      ! DITTO
      REAL RTMP2                      ! DITTO
      REAL SIGSQ                      ! SIG**2
      CHARACTER*140 TWEET             ! MESSAGE FOR LOG FILE
      REAL, ALLOCATABLE :: Z_(:)      ! STANDARDIZED DEVIATE VALUES
      REAL ZL                         ! LEFT LIMIT OF INTEGRATION
      REAL ZR                         ! RIGHT LIMIT OF INTEGRATION
C ***
C *** ------------------------------------------------------------------
C ***
C *** EXECUTABLE CODE
C ***
      WRITE (*,'("PIECEWISE: ",5I5, L5, 2F10.0)') H,K,M,SUBPAT,INDEX,CALC_DERIVS,YSTAR,SIG
      IF (PRESENT(VOCAL) .AND. VOCAL) THEN
         IACT = 1
         WRITE (IOUT,'("PIECEWISE: YSTAR = ",E16.6)') YSTAR
         WRITE (IOUT,'("PIECEWISE: MMI3"/"H,K,M,SUBPAT,IACT,INDEX"/I10,5(",",I10))')
     &    H, K, M, SUBPAT, 1, INDEX
         FORMAT1 = '("h,k,m,subpat,iact,index,%HIST,%ELLMAX,%W,%SPREAD,%WIDTH,%HEIGHT,%NOTE"/
     &    I10,5(",",I10),",",L1,",",I10,4(",",E16.6),",",A20)'
         WRITE(IOUT,FORMAT1)
     &    H,K,M,SUBPAT,IACT,INDEX, MMI3(H,K,M,SUBPAT,IACT,INDEX)%HIST,
     &    MMI3(H,K,M,SUBPAT,IACT,INDEX)%ELLMAX, MMI3(H,K,M,SUBPAT,IACT,INDEX)%W,
     &    MMI3(H,K,M,SUBPAT,IACT,INDEX)%SPREAD,
     &    MMI3(H,K,M,SUBPAT,IACT,INDEX)%WIDTH, MMI3(H,K,M,SUBPAT,IACT,INDEX)%HEIGHT,
     &    MMI3(H,K,M,SUBPAT,IACT,INDEX)%NOTE
         IF (MMI3(H,K,M,SUBPAT,IACT,INDEX)%ELLMAX > 0) THEN
            CALL CONVERT (MMI3(H,K,M,SUBPAT,IACT,INDEX)%ELLMAX, CTMP)
            FORMAT2 = '(A6,'//TRIM(ADJUSTL(CTMP))//'(",",F10.2))'
            WRITE (IOUT,FORMAT2) '%U(yd)', (MMI3(H,K,M,SUBPAT,IACT,INDEX)%U(ELL),
     &       ELL = 1, MMI3(H,K,M,SUBPAT,IACT,INDEX)%ELLMAX)
            FORMAT2 = '(A6,'//TRIM(ADJUSTL(CTMP))//'(",",F10.3))'
            WRITE (IOUT,FORMAT2) '%P', (MMI3(H,K,M,SUBPAT,IACT,INDEX)%P(ELL),
     &       ELL = 1, MMI3(H,K,M,SUBPAT,IACT,INDEX)%ELLMAX)
            FORMAT2 = '(A6,'//TRIM(ADJUSTL(CTMP))//'(",",E16.6))'
            WRITE (IOUT,FORMAT2) '%M', (MMI3(H,K,M,SUBPAT,IACT,INDEX)%M(ELL),
     &       ELL = 1, MMI3(H,K,M,SUBPAT,IACT,INDEX)%ELLMAX)
            WRITE (IOUT,FORMAT2) '%B', (MMI3(H,K,M,SUBPAT,IACT,INDEX)%B(ELL),
     &       ELL = 1, MMI3(H,K,M,SUBPAT,IACT,INDEX)%ELLMAX)
         ENDIF ! %ELLMAX > 0
         CONTINUE ! DEBUG POINT
      ENDIF ! PRESENT(VOCAL) & VOCAL
      HY = 0.
      PHY_PYSTAR = 0.
      P2HY_PYSTAR2 = 0.
      SIGSQ = SIG * SIG
      ELLMAX = MMI3(H,K,M,SUBPAT,1,INDEX)%ELLMAX
      IF (ALLOCATED(P) .OR. ALLOCATED(G) .OR. ALLOCATED(Z_)) THEN
         CONTINUE ! DEBUG POINT
      ENDIF ! ALLOCATED(P)
      ALLOCATE ( P(ELLMAX,0:2), G(ELLMAX,0:2), Z_(ELLMAX) )
      DO ELL = 1, ELLMAX
         Z_(ELL) = Z(YSTAR-MMI3(H,K,M,SUBPAT,1,INDEX)%U(ELL),0.,SIG)
         P(ELL,0) = ANORDF(Z_(ELL))
         G(ELL,0) = ANORPR(Z_(ELL))
      ENDDO ! ELL
      DO ELL = 1, ELLMAX-1
         EM = MMI3(H,K,M,SUBPAT,1,INDEX)%M(ELL)
         COEFF = EM*YSTAR + MMI3(H,K,M,SUBPAT,1,INDEX)%B(ELL)
         MAXZ = MAX(Z_(ELL), Z_(ELL+1)) ! SHOULD ALWAYS = Z_(ELL)
         IF (MAXZ <= 10.) THEN
            DELTA_BIGG = P(ELL,0) - P(ELL+1,0)
         ELSE ! MAXZ > 10 - USE EQUIVALENT FORM
            DELTA_BIGG = ANORDF(-Z_(ELL+1)) - ANORDF(-Z_(ELL))            
         ENDIF ! MAXZ <= 10
         RTMP = COEFF * DELTA_BIGG
         IF (EM /= 0.) THEN
            RTMP = RTMP + EM*SIG*(G(ELL,0)-G(ELL+1,0))
         ENDIF ! EM /= 0
         IF (RTMP < 0.) THEN
            CONTINUE ! DEBUG POINT
         ENDIF ! RTMP < 0.
         HY = HY + RTMP
         IF (CALC_DERIVS) THEN
            ZR = Z_(ELL)
            ZL = Z_(ELL+1)
            RTMP1 = EM * DELTA_BIGG + COEFF * (G(ELL,0) - G(ELL+1,0)) / SIG
            RTMP2 = 2. * EM * (G(ELL,0) - G(ELL+1,0)) / SIG
     &              - COEFF * (ZR*G(ELL,0) - ZL*G(ELL+1,0)) / SIGSQ
            IF (EM /= 0.) THEN
               RTMP1 = RTMP1 - EM *(ZR*G(ELL,0) - ZL*G(ELL+1,0))
               RTMP2 = RTMP2 - EM * ((1.-ZR)*G(ELL,0) - (1.-ZL)*G(ELL+1,0)) / SIG
            ENDIF ! EM /= 0
            PHY_PYSTAR = PHY_PYSTAR + RTMP1
            P2HY_PYSTAR2 = P2HY_PYSTAR2 + RTMP2
         ENDIF ! CALC_DERIVS
      ENDDO ! ELL
      IF (HY < 0.) THEN
         CONTINUE ! DEBUG POINT
      ENDIF ! HY < 0
      IF (.NOT. ALLOCATED(P)) THEN
         TWEET = "PIECEWISE: P NOT ALLOCATED !"
         CALL TWITTER (TWEET, PAWS = .TRUE.)
      ELSEIF (.NOT. ALLOCATED(G)) THEN
         TWEET = "PIECEWISE: G NOT ALLOCATED !"
         CALL TWITTER (TWEET, PAWS = .TRUE.)
      ELSEIF (.NOT. ALLOCATED(Z_)) THEN
         TWEET = "PIECEWISE: Z_ NOT ALLOCATED !"
         CALL TWITTER (TWEET, PAWS = .TRUE.)
      ENDIF ! .NOT. ALLOCATED(P)
      DEALLOCATE ( P, G, Z_ )
      IF (ALLOCATED(P)) THEN
         CONTINUE ! DEBUG POINT
      ELSEIF (ALLOCATED(G)) THEN
         CONTINUE ! DEBUG POINT
      ELSEIF (ALLOCATED(Z_)) THEN
         CONTINUE ! DEBUG POINT
      ENDIF ! ALLOCATED(P)

      RETURN
      END SUBROUTINE PIECEWISE
      
      SUBROUTINE HFY_HY (HMAX, KMAX, MMAX, YSTAR, MSC, MMI3, M, NR, NRUNS, RUN,
     & NPAT, PAT, HFOH, VOCAL)
      USE CONSTANTS
      USE DERIVED_TYPES
      USE EM_IMSL ! CONTAINS "ANORDF" AND "ANORPR"
      USE FILES
      USE PATTERN_DATA
      USE UTILS ! CONTAINS "TWITTER"
      IMPLICIT NONE
C *** 
C *** PURPOSE: PERFORMS HIGH-RESOLUTION CALCULATION OF $h_F(y)/h(y)$
C ***
C *** METHODOLOGY: SERIES EXPANSION
C ***
C *** PENDING MODIFICATIONS ("TODO" ITEMS):
C ***
C *** MODIFICATION HISTORY: VERSION 1.0 - TBD (INITIAL RELEASE)
C ***
C *** ------------------------------------------------------------------
C ***
C *** ARGUMENTS:
C ***                                 !
      INTEGER, INTENT(IN) :: HMAX     ! NUMBER OF MINEFIELD SEGMENTS
      INTEGER, INTENT(IN) :: KMAX     ! NUMBER OF "SWEEPER" TYPES
      INTEGER, INTENT(IN) :: MMAX     ! NUMBER OF MINE TYPES
      REAL, INTENT(IN) :: YSTAR       ! LATERAL DISTANCE FROM TRACK, yd
      TYPE (MSCS), INTENT(IN) :: MSC(HMAX,KMAX) ! MINESWEEPER CHARACTERISTICS DATA
      TYPE (MMI3S), INTENT(IN) :: MMI3(HMAX,KMAX,MMAX,0:1,1,3)
                                      ! TYPE-3 MINESWEEPER-MINE INTERACTION
                                      ! DATA
      INTEGER, INTENT(IN) :: M        ! MINE TYPE
      INTEGER, INTENT(IN) :: NR       ! RUN NUMBER
      INTEGER, INTENT(IN) :: NRUNS    ! TOTAL NUMBER OF RUNS
      TYPE (RUNS), INTENT(IN) :: RUN(0:NRUNS) ! RUN DATA
      INTEGER, INTENT(IN) :: NPAT     ! NUMBER OF PATTERNS
      TYPE (PATTERNS), INTENT(IN) :: PAT(NPAT) ! PATTERN DATA
      REAL, INTENT(OUT) :: HFOH(0:2)  ! $h_F/h$, W/ DERIVS.
      LOGICAL, OPTIONAL, INTENT(IN) :: VOCAL ! INVOKES PRINTED OUTPUT
C ***
C *** ------------------------------------------------------------------
C ***
C ***  LOCAL VARIABLE DEFINITIONS
C ***                                 !
      REAL, ALLOCATABLE :: B(:,:)     ! INTERCEPT DATA
      REAL, ALLOCATABLE :: C(:,:)     ! COEFFICIENT DATA (SECOND TERM)
      REAL CONST                      ! CONSTANT FOR NU-TH TERM IN EXPANSION OF
                                      ! GUASSIAN pmf
      CHARACTER*10 CTMP               ! TEMPORARY (SCRATCH) CHARACTER VARIABLE
      REAL, ALLOCATABLE :: D(:,:,:)   ! COEFFICIENT DATA (FIRST TERM) - FORMERLY "a"
      LOGICAL DEBUG                   ! INVOKES DEBUGGING OUTPUT
      REAL DENOM                      ! DENOMINATOR OF $\nu^\mathrm{th}$ TERM
      REAL DHNU(2,0:2)                ! $\Delta h_\nu(y*)$, W/ DERIVS
      REAL DISCR                      ! DISCRIMINANT
      INTEGER ELL                     ! PIECEWISE DATA INDEX
      INTEGER ELLMAX(2)               ! MAXIMUM ELL VALUES FOR ACTUATION &
                                      ! DAMAGE DATA, RESP'Y
      REAL, ALLOCATABLE :: EM(:,:)    ! SLOPE DATA
      LOGICAL FIRST /.TRUE./          ! INDICATES FIRST CALL
      CHARACTER*80 FORMAT1            ! RUN-TIME OUTPUT FORMAT
      CHARACTER*80 FORMAT2            ! DITTO
C        NB: GPDF AND GPMF ARE \emph{ARRAYS} IN THE CONTEXT OF THIS SUBROUTINE
      REAL, ALLOCATABLE :: GPDF(:,:,:) ! ARRAY OF GAUSSIAN pdf VALUES, W/ DERIVS.
C      REAL, ALLOCATABLE :: GPMF(:,:)  ! ARRAY OF GAUSSIAN pmf VALUES
      INTEGER H                       ! SEGMENT INDEX
      REAL HY(2,0:2)                  ! $h_{(F)}(y)$, W/ DERIVS W/RT $y*$
      INTEGER INDEX                   ! ACTUATION/DAMAGE INDEX
      INTEGER IORD                    ! ORDER OF PARTIAL DERIVATIVE W/RT $y^*$
                                      ! (0 TO 2)
      INTEGER IPAT                    ! PATTERN INDEX
      INTEGER ITMP                    ! TEMPORARY (SCRATCH) INTEGER VARIABLE
      INTEGER K                       ! SWEEPER TYPE INDEX
      LOGICAL LTMP                    ! TEMPORARY (SCRATCH) LOGICAL VARIABLE
      REAL MAXABSZ                    ! MAXIMUM ABSOLUTE VALUE OF Z (ACTUATION
                                      ! DATA ONLY)
      INTEGER NU                      ! TERM NUMBER
      INTEGER NU_MIN                  ! MINIMUM VALUE OF "NU" REQUIRED FOR
                                      ! ITERATIVE CALCULATIONS
      REAL RTMP                       ! TEMPORARY (SCRATCH) REAL VARIABLE
      REAL RTMP2                      ! DITTO
      REAL SIG                        ! STANDARD DEVIATION OF NAVIGATION ERROR, yd
      REAL SIGSQ                      ! VARIANCE OF NAVIGATION ERROR $\sigma^2$, yd^2
      INTEGER SUBPAT                  ! SUB-PATTERN INDEX
      REAL, ALLOCATABLE :: TAU(:,:,:) ! $\tau$ VALUES
      CHARACTER*140 TWEET             ! MESSAGE FOR LOG FILE
      REAL, ALLOCATABLE :: Z_(:,:)    ! STANDARDIZED VARIABLES
      REAL, ALLOCATABLE :: ZSQ(:,:,:) ! ARRAY OF $z^2$ VALUES
C ***
C *** ------------------------------------------------------------------
C ***
C *** EXECUTABLE CODE
C ***
C      CALL PRN_PATTERN (PAT(RUN(NR)%IPAT), .TRUE.)
      WRITE (*,'("HFY_HY: YSTAR = ",E16.6)') YSTAR
      IF (FIRST) THEN
         TWEET = "HFY_HY: REMINDER - DERIVATIVES NOT YET IMPLEMENTED"
         CALL TWITTER (TWEET, PAWS = .TRUE.)
         FIRST = .FALSE.
      ENDIF ! FIRST
      DEBUG = PRESENT(VOCAL) .AND. VOCAL
      IF (DEBUG) WRITE (IOUT,'("HFY_HY: YSTAR = ",F10.2)') YSTAR
      HFOH = 0. ! TEMPORARY CODE
      HY = 0.
      IPAT = RUN(NR)%IPAT
      SUBPAT = RUN(NR)%SUBPAT
      H = PAT(IPAT)%H
      K = PAT(IPAT)%K
      SIG = MSC(H,K)%SIG
      SIGSQ = SIG*SIG
      MAXABSZ = -HUGE(RTMP)
      ALLOCATE ( EM(2,ELLCAP), B(2,ELLCAP), Z_(2,ELLCAP),
C     & GPMF(2,ELLCAP),
     & GPDF(2,ELLCAP,0:2), C(2,ELLCAP), D(2,ELLCAP,0:2), TAU(2,ELLCAP,0:2),
     & ZSQ(2,ELLCAP,0:2) )
      C = 0.
      D = 0.
      DO INDEX = 1, 2
         ELLMAX(INDEX) = MMI3(H,K,M,SUBPAT,1,INDEX)%ELLMAX
         DO ELL = 1, ELLCAP
            IF (ELL <= ELLMAX(INDEX)) THEN
               EM(INDEX,ELL) = MMI3(H,K,M,SUBPAT,1,INDEX)%M(ELL)
               B(INDEX,ELL) = MMI3(H,K,M,SUBPAT,1,INDEX)%B(ELL)
            ENDIF ! ELL <= ELLMAX(INDEX)
         ENDDO ! ELL
         DO ELL = 1, ELLMAX(INDEX)
            IF (YSTAR == -1125.) THEN
               CONTINUE ! DEBUG POINT
            ENDIF ! YSTAR = -1125
            Z_(INDEX,ELL) = (YSTAR - MMI3(H,K,M,SUBPAT,1,INDEX)%U(ELL)) / SIG
            IF (INDEX == 1) THEN
               RTMP = ABS(Z_(INDEX,ELL))
               IF (RTMP > MAXABSZ) MAXABSZ = RTMP
            ENDIF ! RTMP > MAXABSZ
C            GPMF(INDEX,ELL) = ANORDF(Z_(INDEX,ELL))
            GPDF(INDEX,ELL,0) = ANORPR(Z_(INDEX,ELL))
            GPDF(INDEX,ELL,1) = -GPDF(INDEX,ELL,0)*Z_(INDEX,ELL)/SIG
            GPDF(INDEX,ELL,2) = -GPDF(INDEX,ELL,1)*Z_(INDEX,ELL)/SIG - GPDF(INDEX,ELL,0)/SIGSQ
!             NB:  SECOND PARTIALS OF "D" NOT ASSIGNED; $\equiv$ 0
            IF (ELL == 1) THEN
               RTMP = EM(INDEX,ELL)
               C(INDEX,ELL) = RTMP * SIG
               D(INDEX,ELL,0) = RTMP*YSTAR + B(INDEX,ELL)
               D(INDEX,ELL,1) = RTMP
            ELSEIF (ELL <= ELLMAX(INDEX)-1) THEN
               RTMP = EM(INDEX,ELL) - EM(INDEX,ELL-1)
               C(INDEX,ELL) = RTMP * SIG
               D(INDEX,ELL,0) = RTMP*YSTAR + B(INDEX,ELL) - B(INDEX,ELL-1)
               D(INDEX,ELL,1) = RTMP               
            ELSE ! ELL = ELLMAX
               RTMP = -EM(INDEX,ELL-1)
               C(INDEX,ELL) = RTMP * SIG
               D(INDEX,ELL,0) = RTMP*YSTAR - B(INDEX,ELL-1)
               D(INDEX,ELL,1) = RTMP
            ENDIF ! ELL = 1
            DO IORD = 0, 2
               HY(INDEX,IORD) = HY(INDEX,IORD) + C(INDEX,ELL)*GPDF(INDEX,ELL,IORD)
            ENDDO ! IORD
         ENDDO ! ELL
      ENDDO ! INDEX
      IF (DEBUG) WRITE (IOUT,'("HFY_HY: HY(1;0) =",E16.6)') HY(1,0)
      
      IF (MAXABSZ <= 10.) THEN ! FORMERLY "<= ZMAX"
         IF (DEBUG) THEN
              WRITE (IOUT,'("HFY_HY: ELLMAX =",2I10)') (ELLMAX(INDEX), INDEX = 1,2)
            DO INDEX = 1, 2
               WRITE (IOUT,'("HFY_HY: D")')
               CALL CONVERT (ELLMAX(INDEX), CTMP)
               FORMAT1 = '('//TRIM(ADJUSTL(CTMP))//'(E16.6,","))'
               WRITE (IOUT,FORMAT1) (D(INDEX,ELL,0), ELL = 1, ELLMAX(INDEX))
               WRITE (IOUT,'("HFY_HY: Z_")')
               WRITE (IOUT,FORMAT1) (Z_(INDEX,ELL), ELL = 1, ELLMAX(INDEX))
            ENDDO ! INDEX
            WRITE (IOUT,'("HFY_HY: NU, CONST, HY(1;0), HY(2;0)")')
         ENDIF ! DEBUG
C           FOR THE LOGIC BEHIND THE FOLLOWING, SEE MLM NOTES OF 5/3/20
         RTMP = MAXABSZ * MAXABSZ
         DISCR = 1. - 6.*RTMP + RTMP*RTMP
         IF (DISCR < 0.) THEN
            NU_MIN = 1
         ELSE ! DISCR >= 0
            NU_MIN  = CEILING (1. + 0.25*SQRT(DISCR))
         ENDIF ! DISCR < 0
         NU = 0
         DHNU(1,0) = HUGE(RTMP)
C           "DO UNTIL NU > NU_MIN & |DHNU(1,0)| <= 1.E-6 & HY > 0"
         DO WHILE (NU <= NU_MIN .OR. ABS(DHNU(1,0)) > 1.E-6 .OR. HY(1,0) <= 0.)
            CALL DELTA_HNU (NU, SIG, ELLMAX, D, Z_, ZSQ, TAU, DHNU)
            IF (DEBUG) WRITE (IOUT,'("HFY_HY: DHNU(*;0) =",2E16.6)') (DHNU(INDEX,0), INDEX = 1,2)
            HY = HY + DHNU ! MATRIX ADDITION
!            DO INDEX = 1, 2
!               DO IORD = 0, 2
!                  HY(INDEX,IORD) = HY(INDEX,IORD) + DHNU(INDEX,IORD)
!               ENDDO ! IORD
!            ENDDO ! INDEX
            IF (DEBUG) WRITE (IOUT,'(I10,3(",",E16.6))') NU, HY(1,0), HY(2,0)
            NU = NU + 1
         ENDDO ! WHILE NU <= NU_MIN
         
         HFOH(0) = HY(2,0) / HY(1,0)
         RTMP = HY(1,0) * HY(1,0)
         RTMP2 = RTMP * HY(1,0)
         HFOH(1) = (HY(2,1)*HY(1,0) - HY(2,0)*HY(1,1)) / RTMP2 ! /rtmp???
         HFOH(2) = (RTMP*HY(2,2) - 2.*HY(1,0)*HY(1,1)*HY(2,1)
     &              + 2.*HY(2,0)*HY(1,1)*HY(1,1) - HY(2,0)*HY(1,0)*HY(1,2))
     &              / RTMP2
      ENDIF ! MAXABSZ <= ZMAX
      DEALLOCATE ( EM, B, Z_, GPDF, C, D, TAU, ZSQ )
      RETURN
      END SUBROUTINE HFY_HY
      
      SUBROUTINE DELTA_HNU (NU, SIG, ELLMAX, D, Z_, ZSQ, TAU, DHNU)
      USE CONSTANTS
      USE UTILS
      IMPLICIT NONE
C *** 
C *** PURPOSE: CALCULATES THE $\nu^\mathrm{th}$ TERM OF h_{(F)}(y)
C ***
C *** METHODOLOGY: BASED ON POLYNOMIAL EXPANSION OF GAUSSIAN pmf.  SEE
C *** \citep{NBS1970} EQ. 26.2.10 & NOTES OF 4/28/20.
C ***
C *** PENDING MODIFICATIONS ("TODO" ITEMS):
C ***
C *** MODIFICATION HISTORY: VERSION 1.0 - TBD (INITIAL RELEASE)
C ***
C *** ------------------------------------------------------------------
C ***
C *** ARGUMENTS:
C ***                                 !
      INTEGER, INTENT(IN) :: NU       ! ITERATION NUMBER
      REAL, INTENT(IN) :: SIG         ! STANDARD DEVIATION OF NAVIGATION ERROR, yd
      INTEGER, INTENT(IN) :: ELLMAX(2) ! MAXIMUM \ell VALUE
      REAL, INTENT(IN) :: D(2,ELLCAP,0:2) ! COEFFICIENTS, W/ DERIVS W/RT y*
      REAL, INTENT(IN) :: Z_(2,ELLCAP) ! STANDARDIZED DEVIATE VALUES
      REAL, INTENT(INOUT) :: ZSQ(2,ELLCAP,0:2) ! $z^2$ VALUES
      REAL, INTENT(INOUT) :: TAU(2,ELLCAP,0:2) ! $\tau$ VALUES
      REAL, INTENT(OUT) :: DHNU(2,0:2) ! $\Delta h_\nu$, W/ DERIVS W/RT y* (ALSO
                                     ! USED AS ACCUMULATOR)
C ***
C *** ------------------------------------------------------------------
C ***
C ***  LOCAL VARIABLE DEFINITIONS
C ***                                 !
      CHARACTER*10 CTMP               ! TEMPORARY (SCRATCH) CHARACTER VARIABLE
      LOGICAL DEBUG                   ! INVOKES DEBUGGING PRINTOUT
      INTEGER ELL                     ! INDEX OVER DATA SEGMENTS (z VALUES)
      CHARACTER*80 FORMAT1            ! RUN-TIME OUTPUT FORMAT
      CHARACTER*80 FORMAT2            ! DITTO
      INTEGER INDEX                   ! A/F (h/h_F) INDEX
      INTEGER IORD                    ! ORDER OF DERIV W/RT $y*$
      INTEGER ITMP                    ! TEMPORARY (SCRATCH) INTEGER VARIABLE
      CHARACTER*80 LABEL              ! OUTPUT LABEL
      LOGICAL LTMP                    ! TEMPORARY (SCRATCH) LOGICAL VARIABLE
      REAL RTMP                       ! TEMPORARY (SCRATCH) REAL VARIABLE
      REAL RTMP2                      ! DITTO
      REAL SIGSQ                      ! $\sigma^2$
      CHARACTER*140 TWEET             ! MESSAGE FOR LOG FILE
      REAL ZTMP(0:2)                  ! $z_\ell^{2\nu+1}$, W/ DERIVS.
C ***
C *** ------------------------------------------------------------------
C ***
C *** EXECUTABLE CODE
C ***
      DHNU = 0.
      RTMP = FLOAT(2*NU + 1)
      RTMP2 = FLOAT(2*NU)
      SIGSQ = SIG * SIG
      DO INDEX = 1, 2
         IF (NU == 0) THEN
            DO ELL = 1, ELLMAX(INDEX)
               TAU(INDEX,ELL,0) = RRTWOPI * Z_(INDEX,ELL)
               ZSQ(INDEX,ELL,0) = Z_(INDEX,ELL) * Z_(INDEX,ELL)
               DHNU(INDEX,0) = DHNU(INDEX,0) + TAU(INDEX,ELL,0)
            ENDDO ! ELL
         ELSE ! NU > 0
            DO ELL = 1, ELLMAX(INDEX)
               TAU(INDEX,ELL,0) = -TAU(INDEX,ELL,0) * ZSQ(INDEX,ELL,0) / RTMP2
               DHNU(INDEX,0) = DHNU(INDEX,0) + TAU(INDEX,ELL,0)/RTMP
            ENDDO ! ELL
         ENDIF ! NU = 0
      ENDDO ! INDEX
C      LABEL = "DELTA_HNU: DHNU"
C      CALL PRN_RMAT(LABEL,DHNU,1,2,0,2)
      
      RETURN
      END SUBROUTINE DELTA_HNU
      
      END MODULE LOADR_N