text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   EVB-QMDFF - RPMD molecular dynamics and rate constant calculations on
!               black-box generated potential energy surfaces
!
!   Copyright (c) 2021 by Julien Steffen (steffen@pctc.uni-kiel.de)
!                         Stefan Grimme (grimme@thch.uni-bonn.de) (QMDFF code)
!
!   Permission is hereby granted, free of charge, to any person obtaining a
!   copy of this software and associated documentation files (the "Software"),
!   to deal in the Software without restriction, including without limitation
!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!   and/or sell copies of the Software, and to permit persons to whom the
!   Software is furnished to do so, subject to the following conditions:
!
!   The above copyright notice and this permission notice shall be included in
!   all copies or substantial portions of the Software.
!
!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
!   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
!   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!   DEALINGS IN THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      subroutine egrad_h3(q,Natoms,Nbeads,V,dVdq,info)
      integer, intent(in) :: Natoms
      integer, intent(in) :: Nbeads
      double precision, intent(in) :: q(3,Natoms,Nbeads)
      double precision, intent(out) :: V(Nbeads)
      double precision, intent(out) :: dVdq(3,Natoms,Nbeads)

      double precision :: R(3), dVdr(3)
      double precision :: xAB, yAB, zAB, rAB
      double precision :: xAC, yAC, zAC, rAC
      double precision :: xBC, yBC, zBC, rBC
      integer k, info
	info = 0
      do k = 1, Nbeads

        xAB = q(1,2,k) - q(1,1,k)
        yAB = q(2,2,k) - q(2,1,k)
        zAB = q(3,2,k) - q(3,1,k)
        rAB = sqrt(xAB * xAB + yAB * yAB + zAB * zAB)
        R(1) = rAB

        xAC = q(1,1,k) - q(1,3,k)
        yAC = q(2,1,k) - q(2,3,k)
        zAC = q(3,1,k) - q(3,3,k)
        rAC = sqrt(xAC * xAC + yAC * yAC + zAC * zAC)
        R(2) = rAC

        xBC = q(1,3,k) - q(1,2,k)
        yBC = q(2,3,k) - q(2,2,k)
        zBC = q(3,3,k) - q(3,2,k)
        rBC = sqrt(xBC * xBC + yBC * yBC + zBC * zBC)
        R(3) = rBC

        call pote(R, V(k), dVdr)

        dVdq(1,1,k) = dVdr(2) * xAC / rAC - dVdr(1) * xAB / rAB
        dVdq(2,1,k) = dVdr(2) * yAC / rAC - dVdr(1) * yAB / rAB
        dVdq(3,1,k) = dVdr(2) * zAC / rAC - dVdr(1) * zAB / rAB

        dVdq(1,2,k) = dVdr(1) * xAB / rAB - dVdr(3) * xBC / rBC
        dVdq(2,2,k) = dVdr(1) * yAB / rAB - dVdr(3) * yBC / rBC
        dVdq(3,2,k) = dVdr(1) * zAB / rAB - dVdr(3) * zBC / rBC

        dVdq(1,3,k) = dVdr(3) * xBC / rBC - dVdr(2) * xAC / rAC
        dVdq(2,3,k) = dVdr(3) * yBC / rBC - dVdr(2) * yAC / rAC
        dVdq(3,3,k) = dVdr(3) * zBC / rBC - dVdr(2) * zAC / rAC
      end do

      end subroutine egrad_h3

      SUBROUTINE pote (R,pe,dpe)
C    *************************
C    *  B K M P 2     P E S  *   Implemented in carlon by L.B.
C    *************************
C
C CALCULATE TOTAL H3 POTENTIAL FROM ALL OF ITS PARTS
C IF (ID.GT.0) ALSO CALCULATE THE DV/DR DERIVATIVES
C ALL DISTANCES ARE IN BOHRS AND ALL ENERGIES ARE IN HARTREES
C
C FOR A DISCUSSION OF THIS SURFACE, SEE:
C    A.I.BOOTHROYD, W.J.KEOGH, P.G.MARTIN, M.R.PETERSON
C    JOURNAL OF CHEMICAL PHYSICS 95 PP4343-4359 (SEPT15/91)
C    AND JCP 104 PP 7139-7152 (MAY8/96)
C
C NOTE: THIS FILE CONTAINS PARAMETERS FOR A SURFACE REFITTED
C       ON JUNE21/95 TO A SET OF SEVERAL THOUSAND AB INITIO
C       POINTS.  THE '706' SURFACE PARAMETERS HAVE BEEN
C       COMMENTED OUT.  THE ROUTINE NAMES HAVE BEEN MODIFIED
C       SLIGHTLY (USUALLY A '95' APPENDED) SO THAT A PROGRAMME
C       COULD EASILY CALL AND COMPARE BOTH OF OUR SURFACES.
C
C NOTE: THE SURFACE PARAMETER VALUES AS PUBLISHED LEAD TO AN
C       ANOMOLOUSLY DEEP VAN DER WAALS WELL FOR A VERY COMPACT
C       H2 MOLECULE (SAY R=0.8).  AFTER THAT PAPER WAS SUBMITTED,
C       THIS PROBLEM WAS FIXED AND THE CORRECTED CBEND COEFFICIENTS
C       ARE USED IN THIS VERSION OF THE SURFACE (VERSION 706).
C       THE OLD COEFFICIENTS ARE STILL IN SUBR.VBCB BUT HAVE BEEN
C       COMMENTED OUT.
C
C ANY QUESTIONS/PROBLEMS/COMMENTS CONCERNING THIS PROGRAMME CAN BE
C ADDRESSED TO :    WKEOGH@ALCHEMY.CHEM.UTORONTO.CA
C (IF NECESSARY, VIA   BOOTHROY@CITA.UTORONTO.CA    BOOTHROY@UTORDOP.BITNET
C                 OR   PGMARTIN@CITA.UTORONTO.CA    PGMARTIN@UTORDOP.BITNET )
C
C VERSION:
C APR12/95 ... PARAMETERS FOR SURFACE850308 ADDED
C JUL27/91 ... SURF706D.OUT CBEND VALUES PUT IN
C----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      dimension R(3),DPE(3)
c     DIMENSION R(3)
      DIMENSION DVLON(3),DVAS(3),DVBNDA(3),DVBNDB(3),
     .          DCAL(3), DCAS(3),DCBNDA(3),DCBNDB(3)
      DIMENSION VB(2,25),CB(2,25)
C     DIMENSION VBP(2,25),CBP(2,25)
C     DIMENSION DB1A(3),DB1B(3)
      DIMENSION DT(3),T(3)
C     DIMENSION DB2(3),DB3(6),EXPASS(4)
      IPR = 0
cC
c      R(1)=RAB
c      R(2)=RBC
c      R(3)=RAC
c     print*,r(1),r(2),r(3),rab,rbc,rac
C
c      print*,r(1),r(2),r(3)
      CALL CHGEOM(R,IVALID)
      IF(IVALID.LT.1)THEN
         VTOT   = 99.0
         DVTOT1 = 0.D0
         DVTOT2 = 0.D0
         DVTOT3 = 0.D0
         WRITE(6,6999) VTOT,DVTOT1,DVTOT2,DVTOT3
         RETURN
      END IF
 6999 FORMAT('INVALID GEOMETRY, V,DV=',4F8.2)
C
C
C  ZERO EVERYTHING TO AVOID ANY 'FUNNY' VALUES:
      VLON  = 0.D0
      VAS   = 0.D0
      VBNDA = 0.D0
      VBNDB = 0.D0
      CAL   = 0.D0
      CAS   = 0.D0
      CBNDA = 0.D0
      CBNDB = 0.D0
      DO I=1,3
         DVLON(I) = 0.D0
         DVAS(I)  = 0.D0
         DVBNDA(I)= 0.D0
         DVBNDB(I)= 0.D0
         DCAL(I)  = 0.D0
         DCAS(I)  = 0.D0
         DCBNDA(I)= 0.D0
         DCBNDB(I)= 0.D0
      END DO
C  ZERO THE VB AND CB ARRAYS (USED ONLY FOR NUMERICAL DERIVATIVES)
      DO I=1,25
         VB(1,I) = 0.D0
         VB(2,I) = 0.D0
         CB(1,I) = 0.D0
         CB(2,I) = 0.D0
      END DO
      CALL H3LOND95( R, VLON, DVLON )
      CALL VASCAL95( R, VAS, DVAS )
C    NOW DO ANY CORRECTIONS REQUIRED FOR COMPACT GEOMETRIES:
      CALL COMPAC95(R,ICOMPC,T,DT)
C    OCT.3/90 COMPACT ROUTINES ONLY CALLED FOR COMPACT GEOMETRIES:
      IF(ICOMPC.GE.1)THEN
         CALL CSYM95  ( R, CAL, DCAL )
         CALL CASYM95 ( R, CAS, DCAS,  1 ,T,DT)
      END IF
      CALL VBCB95(R,ICOMPC,T,DT,IPR,
     .          VBNDA,VBNDB,DVBNDA,DVBNDB,
     .          CBNDA,CBNDB,DCBNDA,DCBNDB)
C    ADD UP THE VARIOUS PARTS OF THE POTENTIAL:
      VTOT = VLON + VAS + VBNDA + VBNDB + CAL + CAS + CBNDA + CBNDB
      VTOT=VTOT
	PE=VTOT
C    ADD UP THE VARIOUS PARTS OF THE DERIVATIVE:
      DVTOT1 =   DVLON(1) + DVAS(1) + DVBNDA(1) + DVBNDB(1)
     .           +  DCAL(1) + DCAS(1) + DCBNDA(1) + DCBNDB(1)
      DVTOT2 =   DVLON(2) + DVAS(2) + DVBNDA(2) + DVBNDB(2)
     .           +  DCAL(2) + DCAS(2) + DCBNDA(2) + DCBNDB(2)
      DVTOT3 =   DVLON(3) + DVAS(3) + DVBNDA(3) + DVBNDB(3)
     .           +  DCAL(3) + DCAS(3) + DCBNDA(3) + DCBNDB(3)

	DPE(1)=DVTOT1
	DPE(2)=DVTOT2
	DPE(3)=DVTOT3
C
      IF(IPR.GT.0) THEN
          WRITE(7,*)
          WRITE(7,*) 'H3TOT ENTER --------------------------------'
          WRITE(7,7000) R,ICOMPC,VTOT,
     .             VLON,VAS,CAL,CAS,VBNDA,VBNDB,CBNDA,CBNDB
          IF(IPR.GT.1)THEN
                   WRITE(7,7100) DVTOT1,DVTOT2,DVTOT3,
     1                           DVLON,DVAS,DVBNDA,DVBNDB,
     .                           DCAL,DCAS,DCBNDA,DCBNDB
          END IF
          WRITE(7,7400) (VB(1,I),I=1,5),(VB(2,I),I=1,5)
          WRITE(7,7500) (CB(1,I),I=1,9),(CB(2,I),I=1,9)
          WRITE(7,7999) VBNDA,VBNDB,CBNDA,CBNDB,
     .                  DVBNDA,DVBNDB,DCBNDA,DCBNDB
          WRITE(7,*) 'EXITING SUBR.H3TOT'
          WRITE(7,*) 'H3TOT EXIT ---------------------------------'
      END IF
      RETURN
 7400 FORMAT('VBA VALUES: ',5(1X,G12.6),/,'VBB VALUES: ',5(1X,G12.6))
 7500 FORMAT('CBA VALUES: ',3(1X,F16.8),/,
     .       '            ',3(1X,F16.8),/,
     .       '            ',3(1X,F16.8),/,
     .       'CBB VALUES: ',3(1X,F16.8),/,
     .       '            ',3(1X,F16.8),/,
     .       '            ',3(1X,F16.8))
 7000 FORMAT(5X,'    R =',3(1X,F16.10),' ICOMPAC = ',I1,/,
     .       5X,'VTOT  = ',F18.12,/,
     .       5X,'VLON  = ',F18.12,'       VAS   = ',G18.12,/,
     .       5X,'CAL   = ',G18.12,'       CAS   = ',G18.12,/,
     .       5X,'VBNDA = ',G18.12,'       VBNDB = ',G18.12,/,
     .       5X,'CBNDA = ',G18.12,'       CBNDB = ',G18.12)
 7100 FORMAT(' DVTOT   = ',3(1X,G18.12),/,
     .       ' DVLON   = ',3(1X,G18.12),/,
     .       ' DVAS    = ',3(1X,G18.12),/,
     .       ' DVBNDA  = ',3(1X,G18.12),/,
     .       ' DVBNDB  = ',3(1X,G18.12),/,
     .       ' DCAL    = ',3(1X,G18.12),/,
     .       ' DCAS    = ',3(1X,G18.12),/,
     .       ' DCBNDA  = ',3(1X,G18.12),/,
     .       ' DCBNDB  = ',3(1X,G18.12))
 7750 FORMAT(3(1X,G16.10))
 7751 FORMAT(/,3(1X,G16.10))
 7999 FORMAT('   VBNDA        VBNDB        CBNDA        CBNDB ',/,
     .        4(E12.6,2X),/,
     .       'DVBNDA: ',3(E16.10,1X),/,
     .       'DVBNDB: ',3(E16.10,1X),/,
     .       'DCBNDA: ',3(E16.10,1X),/,
     .       'DCBNDB: ',3(E16.10,1X))
      END
C
      SUBROUTINE TRIPLET95(R,E3,ID)
C--------------------------------------------------------------------
C APR12/95 SURFACE950308 VALUES ADDED
C OCT04/90 SURFACE626 VALUES ADDED
C  H2 TRIPLET CURVE AND DERIVATIVES:
C  CALCULATES TRIPLET POTENTIAL AND FIRST DERIVATIVE
C  USES TRUHLAR HOROWITZ EQUATION WITH OUR EXTENSION
C  USES THE JOHNSON CORRECTION AT SHORT DISTANCES (R < RR)
C     IF R .GE. RR         USE MODIFIED T/H TRIPLET EQUATION
C     IF R .LE. RL         USE THE JOHNSON CORRECTION
C     IN BETWEEN           USE THE TRANSITION EQUATION
C--------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION E3(3),E1(3)
      PARAMETER( RL=0.95D0, RR=1.15D0 )
      PARAMETER( Z1=1.D0, Z2=2.D0, Z4=4.D0)
C TRIPLET AND JOHNSON VALUES MAR08/95
C     PARAMETER(
C    .  A1= -0.0460730469,A2=-16.3948753611,A3=-27.6026451962,
C    .  A4=  2.0383359014,A5= -6.8887433326,A6=  1.6216005938,
C    .  C1= -0.4110341110310669,C2= -0.0767398928636804,
C    .  C3=  0.4302169435606523)
C TRIPLET VALUES JUN21
      PARAMETER(
     .  A1= -0.0298546962,A2=-23.9604445036,A3=-42.5185569474,
     .  A4=  2.0382390988,A5=-11.5214861455,A6=  1.5309487826,
     .  C1= -0.4106358351531854,C2= -0.0770355790707090,
     .  C3=  0.4303193846943223) !JUN21 FIT
C  TRIPLET AND JOHNSON VALUES FROM FIT621:
C     PARAMETER(A1=-0.0253496194,A2=-29.2302126444,
C    .          A3=-50.7225015503,A4= 2.0452676876 ,
C    .          A5=-12.2408908509,A6= 1.6733157383 )
C     PARAMETER( C1=-0.4170298146519658,
C    .           C2=-0.0746027774843370,C3= 0.4297899952237434 )
C SURFACE PARAMETERS FROM FIT601.OUT
C     PARAMETER(A1=-0.66129429,A2=-1.99434198,A3=-2.37604328,
C    .          A4= 2.08107802,A5=-0.0313032510,A6=3.76546699,
C    .          C1=-0.4222590135447196,C2=-0.0731117796738824,
C    .          C3= 0.4295918082189010 )
      IPR = 0
      E3(3) = 0.D0
      IF(R.GE.RR )THEN
C       MODIFIED TRUHLAR/HOROWITZ TRIPLET EQUATION:
         EXDR = DEXP( -A4*R )
         RSQ  = R*R
         RA6  = R**(-A6)
         E3(1) = A1* ( A2 + R + A3*RSQ + A5*RA6 )*EXDR
C       FIRST DERIVATIVE OF TRIPLET CURVE:
         RA61 = R**(-A6-Z1)
         E3(2) = A1*EXDR*
     .   ( Z1 -A2*A4 +(Z2*A3-A4)*R -A3*A4*RSQ -A5*A6*RA61 -A4*A5*RA6 )
      END IF
      IF(R.LT.RR )THEN
          DR = R- RL
          CALL VH2OPT95(R,E1,2)
          IF( R.LE.RL ) THEN
C          JOHNSON TRIPLET EQUATION:
            E3(1) = E1(1) + C2*DR + C3
            E3(2) = E1(2) + C2
          ELSE
C          TRANSITION EQUATION:
            E3(1) = E1(1) + C1*DR*DR*DR + C2*DR + C3
            E3(2) = E1(2) + 3.D0*C1*DR*DR + C2
          END IF
      END IF
      IF(IPR.GT.0) WRITE(7,7100) E3
 7100 FORMAT('  TRIPLET:  E3 = ',3(1X,F12.8))
      RETURN
      END
C
C
      SUBROUTINE VH2OPT95(R,E,IDERIV)
C-----------------------------------------------------------------
C JUL09/90 ... SUPER DUPER SPEEDY VERSION
C SELF-CONTAINED VERSION OF SCHWENKE'S H2 POTENTIAL
C ALL DISTANCES IN BOHRS AND ALL ENERGIES IN HARTREES
C (1ST DERIV ADDED ON MAY 2 1989; 2ND DERIV ON MAY 28 1989)
C-----------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION E(3)
      PARAMETER(A0=0.03537359271649620,    A1= 2.013977588700072    ,
     .  A2= -2.827452449964767        ,    A3= 2.713257715593500    ,
     .  A4= -2.792039234205731        ,    A5= 2.166542078766724    ,
     .  A6= -1.272679684173909        ,    A7= 0.5630423099212294   ,
     .  A8= -0.1879397372273814       ,    A9= 0.04719891893374140  ,
     . A10= -0.008851622656489644     ,   A11= 0.001224998776243630 ,
     . A12= -1.227820520228028D-04    ,   A13= 8.638783190083473D-06,
     . A14= -4.036967926499151D-07    ,   A15= 1.123286608335365D-08,
     . A16= -1.406619156782167D-10 )
      PARAMETER( R0=3.5284882D0,  DD=0.160979391D0,
     .           C6=6.499027D0,   C8=124.3991D0,    C10=3285.828D0)
C     EDISS = 0.174445D0
      E(1)  = -999.D0
      E(2)  = -999.D0
      E(3)  = -999.D0
      R2  = R  *R
      R3  = R2 *R
      R4  = R3 *R
      R5  = R4 *R
      R6  = R5 *R
      R7  = R6 *R
      R8  = R7 *R
      R9  = R8 *R
      R10 = R9 *R
      R11 = R10*R
      R12 = R11*R
      R13 = R12*R
      R14 = R13*R
      R15 = R14*R
      R02 = R0*R0
      R04 = R02*R02
      R06 = R04*R02
      RR2  = R2 + R02
      RR4  = R4 + R04
      RR6  = R6 + R06
      RR25 = RR2*RR2*RR2*RR2*RR2
C     GENERAL TERM:  A(I)*R(I-1),  I=0,16
      ALPHAR =  A0/R + A1
     .            + A2 *R   + A3 *R2  + A4 *R3  + A5 *R4  + A6 *R5
     .            + A7 *R6  + A8 *R7  + A9 *R8  + A10*R9  + A11*R10
     .            + A12*R11 + A13*R12 + A14*R13 + A15*R14 + A16*R15
      EXALPH = DEXP(ALPHAR)
      VSR = DD*(EXALPH-1.D0)*(EXALPH-1.D0) - DD
      VLR =  -C6/RR6   -C8/(RR4*RR4)   -C10/RR25
      E(1)=  VSR +  VLR
C
C  CALCULATE FIRST DERIVATIVE IF REQUIRED:
      IF(IDERIV.GE.1)THEN
         R3 = R2*R
         R5 = R4*R
         RR26 = RR25*RR2
         RR43 = RR4*RR4*RR4
C        GENERAL TERM:  (I-1)*A(I)*R**(I-2)  ,  I=0,16
         DALPHR = -A0/R2 + A2 + 2.D0*A3*R
     .          +  3.D0 *A4 *R2  +  4.D0 *A5 *R3  +  5.D0 *A6 *R4
     .          +  6.D0 *A7 *R5  +  7.D0 *A8 *R6  +  8.D0 *A9 *R7
     .          +  9.D0 *A10*R8  + 10.D0 *A11*R9  + 11.D0 *A12*R10
     .          + 12.D0 *A13*R11 + 13.D0 *A14*R12 + 14.D0 *A15*R13
     .          + 15.D0 *A16*R14
         DVSR = 2.D0 *DD *(EXALPH-1.D0) *EXALPH *DALPHR
         DVLR =    6.D0*C6*R5 / (RR6*RR6)
     .          +  8.D0*C8*R3 / RR43  + 10.D0*C10*R / RR26
         E(2) = DVSR + DVLR
      END IF
C
C  CALCULATE SECOND DERIVATIVE IF REQUIRED:
      IF(IDERIV.GE.2)THEN
         R10  = R6*R4
         RR27 = RR26*RR2
         RR44 = RR43*RR4
         RR62 = RR6*RR6
         RR63 = RR62*RR6
C        GENERAL TERM: (I-1)*(I-2)*A(I)*R**(I-3),    I=0,16
         DDALPH = 2.D0*A0/R3 +  2.D0*A3      +  6.D0*A4 *R
     .     + 12.D0*A5 *R2    + 20.D0*A6 *R3  + 30.D0*A7 *R4
     .     + 42.D0*A8 *R5    + 56.D0*A9 *R6  + 72.D0*A10*R7
     .     + 90.D0*A11*R8    +110.D0*A12*R9  +132.D0*A13*R10
     .     +156.D0*A14*R11   +182.D0*A15*R12 +210.D0*A16*R13
         DDVSR = 2.D0 *DD *EXALPH
     .   *( (2.D0*EXALPH-1.D0)*DALPHR*DALPHR + (EXALPH-1.D0)*DDALPH )
         DDVLR =- 72.D0* C6 *R10 / RR63 -96.D0 *C8 *R6 / RR44
     .          -120.D0*C10 *R2  / RR27 +30.D0 *C6 *R4 / RR62
     .          + 24.D0* C8 *R2  / RR43 +10.D0 *C10    / RR26
         E(3) = DDVSR + DDVLR
      END IF
      RETURN
      END
C
      SUBROUTINE H3LOND95(R,VLON,DVLON)
C----------------------------------------------------------------------
C VERSION OF MAY 12/90  ... DERIVATIVES CORRECTED (0.5 CHANGED TO 0.25)
C CALCULATES THE H3 LONDON TERMS AND DERIVATIVES
C MODIFIED OCT 7/89 TO INCLUDE EPS**2 TERM WHICH ROUNDS OFF THE
C CUSP IN THE H3 POTENTIAL WHICH OCCURS AT EQUILATERAL TRIANGLE
C CONFIGURATIONS
C----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL*8 Q(3),J(3),JT
      DIMENSION R(3),E1(3),E3(3),ESING(3),ETRIP(3),DVLON(3)
      DIMENSION DE1(3),DE3(3)
C     DIMENSION DJ1(3),DJ2(3),DJ3(3)
      PARAMETER( HALF=0.5D0, TWO=2.D0 , EPS2=1.D-12 )
      IPR = 0
      DO I=1,3
         CALL VH2OPT95(R(I),ESING,2)
          E1(I) = ESING(1)
         DE1(I) = ESING(2)
         CALL TRIPLET95(R(I),ETRIP,2)
         IF(IPR.GT.0) WRITE(7,7400) I,ESING,ETRIP
          E3(I) = ETRIP(1)
         DE3(I) = ETRIP(2)
         Q(I)  = HALF*(E1(I) + E3(I))
         J(I)  = HALF*(E1(I) - E3(I))
      END DO
      SUMQ  =   Q(1) + Q(2) + Q(3)
      SUMJ  =   DABS( J(2)-J(1) )**2
     .        + DABS( J(3)-J(2) )**2
     .        + DABS( J(3)-J(1) )**2
      JT     = HALF*SUMJ + EPS2
      ROOTJT = DSQRT(JT)
      VLON   = SUMQ - ROOTJT
      IF(IPR.GT.0) THEN
         WRITE(7,7410) SUMQ,SUMJ
         WRITE(7,7420) VLON,ROOTJT
      END IF
C  CALCULATE THE DERIVATIVES WITH RESPECT TO R(I):
      DVLON(1) = HALF*(DE1(1)+DE3(1))
     .         - 0.25D0*(TWO*J(1)-J(2)-J(3))*(DE1(1)-DE3(1))/ROOTJT
      DVLON(2) = HALF*(DE1(2)+DE3(2))
     .         - 0.25D0*(TWO*J(2)-J(3)-J(1))*(DE1(2)-DE3(2))/ROOTJT
      DVLON(3) = HALF*(DE1(3)+DE3(3))
     .         - 0.25D0*(TWO*J(3)-J(1)-J(2))*(DE1(3)-DE3(3))/ROOTJT
      IF(IPR.GT.0) THEN
         WRITE(7,7000) R,E1,E3,VLON
         WRITE(7,7100) Q,J
         WRITE(7,7200) DVLON
      END IF
 7000 FORMAT('             R = ',3(1X,F12.6),/,
     .       '            E1 = ',3(1X,F12.8),/,
     .       '            E3 = ',3(1X,F12.8),/,
     .       '          VLON = ',1X,F12.8)
 7100 FORMAT(13X,'Q = ',3(1X,E12.6),/,13X,'J = ',3(1X,E12.6))
 7200 FORMAT('         DVLON = ',3(1X,G12.6))
 7400 FORMAT('FROM SUBR.LONDON: ',/,
     .       '  USING R',I1,':   ESINGLET=',3(1X,F12.8),/,
     .       '         ',1X,'    ETRIPLET=',3(1X,F12.8))
 7410 FORMAT('         SUMQ = ',G12.6,'        SUMJ = ',G12.6)
 7420 FORMAT('         VLON = ',F12.8,'      ROOTJT = ',G12.6)
      RETURN
      END
C
      SUBROUTINE VASCAL95(RPASS,VAS,DVAS)
C------------------------------------------------------------------
C VERSION OF APR12/95  950308 VALUES
C VERSION OF OCT11/90  FIT632C.OUT VALUES
C VERSION OF OCT5/90   SURF626 VALUES
C  CALCULATE THE ASYMMETRIC CORRECTION TERM AND ITS DERIVATIVES
C  SEE EQUATIONS [14] TO [16] OF TRUHLAR/HOROWITZ 1978 PAPER
C------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION DVAS(3),DA(3),DS(3),RPASS(3)
C VASYM VALUES JUN21
      PARAMETER(
     . AA1=0.3788951192E-02, AA2=0.1478100901E-02,
     . AA3=-.1848513849E-03, AA4=0.9230803609E-05,
     . AA5=-.1293180255E-06, AA6=0.5237179303E+00,
     . AA7=-.1112326215E-02) !JUN21 FIT
C--- VASYM VALUES MAR08.95
C     PARAMETER(
C    . AA1=0.3759731624E-02, AA2=0.1476254095E-02,
C    . AA3=-.1866759453E-03, AA4=0.9218646237E-05,
C    . AA5=-.1287906069E-06, AA6=0.5201790843E+00,
C    . AA7=-.1062909514E-02)
C--- VASYM VALUES FROM FIT632C ----------------------
C     PARAMETER(
C    . AA1=0.3438222224E-02,AA2=0.1398145763E-02,
C    . AA3=-.1923999449E-03,AA4=0.9712737075E-05,
C    . AA5=-.1263794562E-06,AA6=0.5181432712E+00,
C    . AA7=-.9487002995E-03 )
      IPR = 0
      R1 = RPASS(1)
      R2 = RPASS(2)
      R3 = RPASS(3)
      R  = R1 + R2 + R3
      RSQ = R*R
      RCU = RSQ*R
C   CALCULATE THE VAS TERM FIRST (EQ.14 OF TRUHLAR/HOROWITZ)
      CALL ACALC95(R1,R2,R3,A,DA)
        A2 = A *A
        A3 = A2*A
        A4 = A3*A
        A5 = A4*A
        EXP1 = DEXP(-AA1*RCU)
        EXP6 = DEXP(-AA6*R)
        S = AA2*A2 + AA3*A3 + AA4*A4 + AA5*A5
      VAS = S*EXP1 +  AA7*A2 *EXP6 / R
         DS(1)  = ( 2.D0*AA2*A  +3.D0*AA3*A2
     .             +4.D0*AA4*A3 +5.D0*AA5*A4) * DA(1)
         DVAS(1) =  -3.D0*AA1*RSQ*S*EXP1 + DS(1)*EXP1
     .             -AA7*A2*EXP6/RSQ  + 2.D0*AA7*A*DA(1)*EXP6/R
     .             -AA6*AA7*A2*EXP6/R
         DS(2)  = ( 2.D0*AA2*A  +3.D0*AA3*A2
     .             +4.D0*AA4*A3 +5.D0*AA5*A4) * DA(2)
         DVAS(2) =  -3.D0*AA1*RSQ*S*EXP1 + DS(2)*EXP1
     .             -AA7*A2*EXP6/RSQ  + 2.D0*AA7*A*DA(2)*EXP6/R
     .             -AA6*AA7*A2*EXP6/R
         DS(3)  = ( 2.D0*AA2*A  +3.D0*AA3*A2
     .             +4.D0*AA4*A3 +5.D0*AA5*A4) * DA(3)
         DVAS(3) =  -3.D0*AA1*RSQ*S*EXP1 + DS(3)*EXP1
     .             -AA7*A2*EXP6/RSQ  + 2.D0*AA7*A*DA(3)*EXP6/R
     .             -AA6*AA7*A2*EXP6/R
      IF(IPR.GT.0) WRITE(7,7100) VAS,DS,DVAS
 7100 FORMAT('    FROM SUBR.VASCAL:   VAS  = ',  1X,G12.6, /,
     .       '                        DS   = ',3(1X,G12.6),/,
     .       '                        DVAS = ',3(1X,G12.6))
      RETURN
      END
C
      SUBROUTINE ACALC95(R1,R2,R3,A,DA)
C---------------------------------------------------------------------
C  VERSION OF MAY 1, 1990
C  BASED ON EQUATIONS FROM NOTES DATE APRIL 4, 1990
C  CALCULATE THE A VALUE AND ITS DERIVATIVES
C           A  = DABS[ (R1-R2)*(R2-R3)*(R3-R1) ]
C---------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION DA(3)
      IPR = 0
      A = (R1-R2)*(R2-R3)*(R3-R1)
      DA(1) = ( -2.D0*R1 + R2 + R3 )*(R2-R3)
      DA(2) = ( -2.D0*R2 + R3 + R1 )*(R3-R1)
      DA(3) = ( -2.D0*R3 + R1 + R2 )*(R1-R2)
      IF(A.LT.0.D0)THEN
         A = -A
         DA(1) = -DA(1)
         DA(2) = -DA(2)
         DA(3) = -DA(3)
      END IF
      IF(IPR.GT.0) WRITE(7,7100) A,DA
 7100 FORMAT('ACALC ENTER ---------------------------------------',/,
     .       '     A = ',F15.10,'   DA = ',3(1X,G12.6),/,
     .       'ACALC EXIT ----------------------------------------')
      RETURN
      END
C
      SUBROUTINE COMPAC95(R,ICOMPC,T,DT)
C---------------------------------------------------------------
C  VERSION OF MAY 14/90 ... CALCULATES T AND DT VALUES ALSO
C  DECIDE WHETHER OR NOT THIS PARTICULAR GEOMETRY IS COMPACT,
C  THAT IS, ARE ANY OF THE THREE DISTANCES SMALLER THAN THE
C  RR VALUE FROM THE JOHNSON CORRECTION.
C-----------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION R(3),T(3),DT(3)
      PARAMETER( RR = 1.15D0, RP = 1.25D0 )
C  FIRST SEE IF THIS IS A COMPACT GEOMETRY:
      DO I=1,3
         T(I) = 0.D0
         DT(I)= 0.D0
      END DO
      ICOMPC = 0
      IF(R(1).LT.RR) ICOMPC = ICOMPC + 1
      IF(R(2).LT.RR) ICOMPC = ICOMPC + 1
      IF(R(3).LT.RR) ICOMPC = ICOMPC + 1
      IF(ICOMPC.EQ.0) RETURN
C CALCULATE THE T(I) VALUES:
      DO I=1,3
      IF(R(I).LT.RR) THEN
         TOP = RR-R(I)
         BOT = RP-R(I)
         TOP2 = TOP  * TOP
         TOP3 = TOP2 * TOP
         BOT2 = BOT  * BOT
         T(I)  = TOP3/BOT
         DT(I) = -3.D0*TOP2/BOT + TOP3/BOT2
      END IF
      END DO
      RETURN
      END
C
      SUBROUTINE CASYM95( R,CAS,DCAS,ID,T,DT )
C---------------------------------------------------------------
C APR12/95 ... SURFACE950308 PARAMETERS ADDED
C VERSION OF SEPT14/90
C  THE COMPACT ASYMMETRIC CORRECTION TERM AND DERIVATIVES
C---------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION R(3),T(3)
      DIMENSION DPR(3),DT(3),DSUMT(3),DTERM1(3),DETERM(3),DSERIES(3),
     .          DCAS(3),DA(3)
C CASYM VALUES JUN21
      PARAMETER(
     . U1=0.2210243144E+00, U2=0.4367417579E+00, U3=0.6994985432E-02,
     . U4=0.1491096501E+01, U5=0.1602896673E+01, U6=-.2821747323E+01,
     . U7=0.4948310833E+00, U8=-.3540394679E-01, U9=-.3305809954E+01,
     .U10=0.3644382172E+01,U11=-.9997570970E+00,U12=0.7989919534E-01,
     .U13=-.1075807322E-02) !JUN21 FIT
C---- CASYM VALUES MAR08.95
C     PARAMETER(
C    . U1=0.5120831287E+00, U2=0.1002277242E+01, U3=0.6850007419E-02,
C    . U4=-.2038751706E+01, U5=0.7027811909E+01, U6=-.4881767278E+01,
C    . U7=0.8801769106E+00, U8=-.6296419648E-01, U9=-.8125516783E+01,
C    .U10=0.6073964424E+01,U11=-.1451402523E+01,U12=0.1165084183E+00,
C    .U13=-.1176579871E-02)
C---- CASYM VALUES FROM FIT633A
C     PARAMETER(
C    . U1=0.1537481166E+00, U2=0.2745950036E+00, U3=0.7501206780E-02,
C    . U4=0.3119136023E+01, U5=0.9969170798E+00, U6=-.3373682823E+01,
C    . U7=0.6807215913E+00, U8=-.4920325491E-01, U9=-.3919467989E+01,
C    .U10=0.5085532326E+01,U11=-.1415264778E+01,U12=0.1138681785E+00,
C    .U13=-.1525367566E-02)
      IPR = 0
      CAS = 0.D0
      DCAS(1) = 0.D0
      DCAS(2) = 0.D0
      DCAS(3) = 0.D0
      CALL ACALC95(R(1),R(2),R(3),A,DA)
      A2 = A*A
C     IF(A.EQ.0.D0) RETURN
      SUMT = T(1) + T(2) + T(3)
      SR   = R(1) + R(2) + R(3)
      PR   = R(1) * R(2) * R(3)
      SR2  = SR*SR
      SR3  = SR2*SR
      PR2  = PR*PR
      PR3  = PR2*PR
C    WRITE OUT THE SERIES EXPLICITLY:
      SERIES = 1.D0 + U4/PR2 +  U5/PR + U6 +   U7*PR +  U8*PR2
     .           + A*(U9/PR2 + U10/PR + U11 + U12*PR + U13*PR2)
      TERM1  = U1/PR**U2
      ETERM  = DEXP(-U3*SR3)
      CAS    = SUMT * A2 * TERM1 * SERIES * ETERM
      IF(IPR.GT.0) THEN
            WRITE(7,*)
            WRITE(7,*) 'CASYM: ------------------'
            WRITE(7,*) ' T(I) = ',T
            WRITE(7,*) 'DT(I) = ',DT
            WRITE(7,*) ' ID = ',ID
            WRITE(7,7400) SERIES,TERM1,ETERM,CAS
            WRITE(7,*) ' CAS = ',CAS
      END IF
      IF(ID.GT.0)THEN
            DPR(1) = R(2)*R(3)
            DPR(2) = R(3)*R(1)
            DPR(3) = R(1)*R(2)
         DO I=1,3
            DSUMT(I) = DT(I)
            DTERM1(I)=  -1.D0*U1*U2*PR**(-U2-1.D0)*DPR(I)
            DETERM(I)= ETERM *(-3.D0*U3*SR2)
            DSERIES(I)=
     .        DPR(I)*( -2.D0*U4/PR3  -U5/PR2 +U7  +2.D0*U8*PR    )
     .        +DA(I)*( U9/PR2 +U10/PR +U11 +U12*PR +U13*PR2    )
     .        +A*DPR(I)*( -2.D0*U9/PR3 -U10/PR2 +U12 + 2.D0*U13*PR )
            DCAS(I) =
     .        DSUMT(I)      * A2    * TERM1 * SERIES * ETERM
     .      + 2.D0*A*DA(I)  * SUMT  * TERM1 * SERIES * ETERM
     .      + DTERM1(I)     * SUMT  * A2    * SERIES * ETERM
     .      + DSERIES(I)    * SUMT  * A2    * TERM1  * ETERM
     .      + DETERM(I)     * SUMT  * A2    * TERM1  * SERIES
         END DO
         IF(IPR.GT.0) THEN
            WRITE(7,7500) DPR,DT,DTERM1,DETERM,DSERIES,DCAS
            WRITE(7,*) ' DCAS = ',DCAS
         END IF
      END IF
      RETURN
 7400 FORMAT('SUBR.CAS:   SERIES = ',G12.6,'     TERM1 = ',G12.6,/,
     .       '             ETERM = ',G12.6,'      CAS  = ',G12.6)
 7500 FORMAT('    DERIVATIVES: DPR = ',3(1X,G12.6),/,
     .       '                 DT  = ',3(1X,G12.6),/,
     .       '              DTERM1 = ',3(1X,G12.6),/,
     .       '              DETERM = ',3(1X,G12.6),/,
     .       '             DSERIES = ',3(1X,G12.6),/,
     .       '              DCAS   = ',3(1X,G12.6))
      END
C
      SUBROUTINE CSYM95(R,CAL,DCAL)
C---------------------------------------------------------------
C VERSION OF APR12/95  950308 VALUES
C VERSION OF OCT12/90  SURF636 VALUES
C  CALCULATE THE 'COMPACT ALL' CORRECTION TERM AND DERIVATIVES
C  A CORRECTION TERM (ADDED SEPT 11/89), WHICH ADDS A SMALL
C  CORRECTION TO ALL COMPACT GEOMETRIES
C---------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION R(3),DCAL(3),G(3)
C     DIMENSION T(3)
      DIMENSION DG(3),SUMV(3)
      PARAMETER( RR=1.15D0, RP=1.25D0 )
C CSYM VALUES JUN21
      PARAMETER(
     . V1=-.2071708868E+00, V2=-.5672350377E+00, V3=0.9058780367E-02) !JUN21 FIT
C---- CSYM VALUES MAR08/95
C     PARAMETER(
C    . V1=-.2210049400E+00, V2=-.7054469608E+00, V3=0.4088898394E-02)
C---- CSYM VALUES FROM FIT633A
C     PARAMETER(
C    . V1=-.2070776049E+00, V2=-.5350898737E+00, V3=0.1011861942E-01)
      CAL   = 0.D0
      IPR   = 0
      SR    = R(1)+R(2)+R(3)
      SR2   = SR*SR
      SR3   = SR*SR2
      EXP3  = DEXP( -V3*SR3 )
      DEXP3 = -3.D0*V3*SR2*EXP3
      DO I=1,3
         RI      = R(I)
         RRRI    = RR-RI
         RRRI2   = RRRI*RRRI
         RRRI3   = RRRI*RRRI2
         RPRI    = RP-RI
         RPRI2   = RPRI*RPRI
         G(I)    = 0.D0
         DG(I)   = 0.D0
         SUMV(I) = V1+V1*V2*RI
         IF(RI.LT.RR) THEN
            G(I)  = (RRRI3/RPRI) *SUMV(I)
            DG(I) = (RRRI3/RPRI2)*SUMV(I)
     .                -3.D0*(RRRI2/RPRI)*SUMV(I)
     .                + (RRRI3/RPRI)*V1*V2
         END IF
      END DO
      SUMG = G(1) + G(2) + G(3)
      CAL  = SUMG*EXP3
      DCAL(1) = DG(1)*EXP3 + SUMG*DEXP3
      DCAL(2) = DG(2)*EXP3 + SUMG*DEXP3
      DCAL(3) = DG(3)*EXP3 + SUMG*DEXP3
      IF(IPR.GT.0)THEN
         WRITE(7,*)
         WRITE(7,7100) CAL
         WRITE(7,7000) G,DG,SUMV,DCAL
         WRITE(7,*) 'EXITING SUBR.CSYM'
      END IF
 7100 FORMAT('CSYM ENTER ----------------------------------------',/,
     .       '    CAL = ',1X,G20.14)
 7000 FORMAT('      G = ',3(1X,G12.6),/,
     .       '     DG = ',3(1X,G12.6),/,
     .       '   SUMV = ',3(1X,G12.6),/,
     .       '   DCAL = ',3(1X,G16.10),/,
     .       'CSYM EXIT -----------------------------------------')
      RETURN
      END
C
      SUBROUTINE VBCB95(RPASS,ICOMPC,T,DT,IPR,
     .                  VBNDA,VBNDB,DVBNDA,DVBNDB,
     .                  CBNDA,CBNDB,DCBNDA,DCBNDB)
C---------------------------------------------------------------
C APR12/95 SURFACE950308 VALUES ADDED
C JUL24/91 FIT705 CBEND VALUES ADDED
C FEB27/91 MODIFIED TO MATCH EQUATION IN H3 PAPER MORE CLOSELY
C NOV 4/90 VBEND COEFFICIENTS NOW A,G  CBEND COEFF'S STILL C,D
C IN THIS VERSION, THE DERIVATIVES ARE ALWAYS CALCULATED
C SUBROUTINES VBEND, CBEND AND B1AB ALL COMBINED INTO THIS ONE
C MODULE IN ORDER TO IMPROVE EFFICIENCY BY NOT HAVING TO PASS
C AROUND B1A,B1B,B2,B3A,B3B FUNCTIONS AND DERIVATIVES
C  B1A = 1 - SUM OF [ P1(COS(THETA(I))) ]
C  B1B = 1 - SUM OF [ P3(COS(THETA(I))) ]
C--------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C     ARRAYS REQUIRED FOR B1,B2,B3 CALCULATIONS:
      DIMENSION RPASS(3),DB1A(3),DB1B(3),DB2(3),DB3(6)
C     DIMENSION TH(3)
C     ARRAYS REQUIRED FOR VBEND CALCULATIONS:
      DIMENSION DVBNDA(3),DVBNDB(3)
      DIMENSION DVB1(3),DVB2(3),DVB3(3),DVB4(3),DVB5(3)
C     DIMENSION DVBND(3)
      DIMENSION VBND(2)
      DIMENSION VB(2,25)
C     ARRAYS REQUIRED FOR CBEND CALCULATIONS:
      DIMENSION CBNDS(2),DCBNDS(2,3)
      DIMENSION DCB1(3),DCB2(3),DCB3(3),DCB4(3),DCB5(3),DCB6(3),
     .          DCB7(3),DCB8(3)
      DIMENSION T(3),DP(3)
      DIMENSION DT(3),DCBNDA(3),DCBNDB(3)
      DIMENSION CB(2,25)
      PARAMETER( Z58=0.625D0, Z38=0.375D0 )
C     PARAMETERS REQUIRED FOR VBEND CALCULATIONS:
      PARAMETER(BETA1=0.52D0, BETA2=0.052D0, BETA3=0.79D0 )
C VBENDA TERMS JUN21         EXPTERMS= 0.52D0, 0.052D0, 0.79D0
      PARAMETER(
     .A11=-.1838073394E+03,A12=0.1334593242E+02,A13=-.2358129537E+00,
     .A21=-.4668193478E+01,A22=0.7197506670E+01,A23=0.2162004275E+02,
     .A24=0.2106294028E+02,A31=0.4242962586E+01,A32=0.4453505045E+01,
     .A41=-.1456918088E+00,A42=-.1692657366E-01,A43=0.1279520698E+01,
     .A44=-.4898940075E+00,A51=0.1742295219E+03,A52=0.3142175348E+02,
     .A53=0.5152903406E+01) !JUN21 FIT
C VBENDB TERMS JUN21
      PARAMETER(
     .G11=-.4765732725E+02,G12=0.3648933563E+01,G13=-.7141145244E-01,
     .G21=0.1002349176E-01,G22=0.9989856329E-02,G23=-.4161953634E-02,
     .G24=0.9075807910E-03,G31=-.2693628729E+00,G32=-.1399065763E-01,
     .G41=-.1417634346E-01,G42=-.4870024792E-03,G43=0.1312231847E+00,
     .G44=-.4409850519E-01,G51=0.5382970863E+02,G52=0.4587102824E+01,
     .G53=0.1768550515E+01) !JUN21 FIT
C    VBENDA MAR08/95   EXPTERMS= 0.52D0, 0.052D0, 0.79D0
C     PARAMETER(
C    .A11=-.1779469255E+03,A12=0.1292511538E+02,A13=-.2284450520E+00,
C    .A21=-.4671094689E+01,A22=0.7810423901E+01,A23=0.2359968959E+02,
C    .A24=0.2293876132E+02,A31=0.4191795902E+01,A32=0.4456025797E+01,
C    .A41=-.1419862958E+00,A42=-.1667760578E-01,A43=0.1246997125E+01,
C    .A44=-.4777517290E+00,A51=0.1681125268E+03,A52=0.3067567046E+02,
C    .A53=0.4942041527E+01)
C    VBENDB MAR08/95
C     PARAMETER(
C    .G11=-.4677295200E+02,G12=0.3587873970E+01,G13=-.7051720560E-01,
C    .G21=0.9549632701E-02,G22=0.9971623818E-02,G23=-.4040392363E-02,
C    .G24=0.9045746595E-03,G31=-.2746842691E+00,G32=-.1412945503E-01,
C    .G41=-.1337907240E-01,G42=-.6200741002E-03,G43=0.1257251291E+00,
C    .G44=-.4195584422E-01,G51=0.5287171694E+02,G52=0.4490366264E+01,
C    .G53=0.1735452866E+01)
C
C CBENDA TERMS JUN21/95
      PARAMETER(
     .C11=0.1860299931E+04,C12=-.6134458037E+03,C13=0.7337207161E+02,
     .C14=-.2676717625E+04,C15=0.1344099415E+04,C21=0.1538913137E+03,
     .C22=0.4348007369E+02,C23=0.1719720677E+03,C24=0.2115963042E+03,
     .C31=-.7026089414E+02,C32=-.1300938992E+03,C41=0.1310273564E+01,
     .C42=-.6175149574E+00,C43=-.2679089358E+02,C44=0.5577477171E+01,
     .C51=-.3543353539E+04,C52=-.3740709591E+03,C53=0.7979303144E+02,
     .C61=-.1104230585E+04,C62=0.4603572025E+04,C63=-.5593496634E+04,
     .C71=-.1069406434E+02,C72=0.1021807153E+01,C73=0.6669828341E-01,
     .C74=0.4168542348E+02,C75=0.1751608567E+02,C81=0.9486883238E+02,
     .C82=-.1519334221E+02,C83=0.4024697252E+04,C84=-.2225159395E+02) !JUN21
C           (LAST 4 PARAMETERS RENUMBERED)
C CBENDB TERMS JUN21/95
      PARAMETER(
     .D11=0.4203543357E+03,D12=-.4922474096E+02,D13=0.3362942544E+00,
     .D14=-.3827423082E+03,D15=0.1746726001E+03,D21=0.1699995737E-01,
     .D22=0.1513036778E-01,D23=0.2659119354E-01,D24=-.5760387483E-02,
     .D31=0.1020622621E+02,D32=0.1050536271E-01,D41=0.6836172780E+00,
     .D42=-.1627858240E+00,D43=-.6925485045E+01,D44=0.1632567385E+01,
     .D51=0.1083595009E+04,D52=0.4641431791E+01,D53=-.8233144461E+00,
     .D61=-.6157225942E+02,D62=0.3094361471E+03,D63=-.3299631143E+03,
     .D71=0.8866227120E+01,D72=-.1382126854E+01,D73=0.7620770145E-01,
     .D74=-.5145757859E+02,D75=0.2046097265E+01,D81=0.2540775558E+01,
     .D82=-.4889246569E+00,D83=-.1127439280E+04,D84=-.2269932295E+01) !JUN21
C           (LAST 4 PARAMETERS RENUMBERED)
C
C    CBENDA MAR08/95
C     PARAMETER(
C    .C11=0.1983833377E+04,C12=-.7161674985E+03,C13=0.9480354622E+02,
C    .C14=-.2860199829E+04,C15=0.1424701828E+04,C21=0.1544592401E+03,
C    .C22=0.3684848293E+02,C23=0.1717816399E+03,C24=0.2060270139E+03,
C    .C31=-.1215015279E+03,C32=-.1370574974E+03,C41=0.8903858200E+00,
C    .C42=-.6496543267E+00,C43=-.2367464822E+02,C44=0.5098769966E+01,
C    .C51=0.6728307258E+03,C52=-.3825078986E+03,C53=0.8128185587E+02,
C    .C61=-.1141379097E+04,C62=0.4771127549E+04,C63=-.5713570141E+04,
C    .C71=-.7394045711E+01,C72=0.9888694698E+00,C73=0.7493009017E-01,
C    .C74=0.9590644736E+02,C75=0.1925781563E+02,C81=0.1001295468E+03,
C    .C82=-.1533150643E+02,C83=-.1850297943E+03,C84=-.2420454635E+02)
C     (LAST 4 PARAMETERS RENUMBERED)
C    CBENDB TERMS MAR08/95
C     PARAMETER(
C    .D11=0.4441478933E+03,D12=-.5345816836E+02,D13=0.1443665553E+01,
C    .D14=-.4155134339E+03,D15=0.1912235377E+03,D21=-.3697316400E-02,
C    .D22=0.1688497284E-01,D23=0.2330055461E-01,D24=-.4670024282E-02,
C    .D31=0.6492099426E+01,D32=0.1597054504E-01,D41=0.6701003237E+00,
C    .D42=-.1868909419E+00,D43=-.6709596498E+01,D44=0.1713120242E+01,
C    .D51=-.5022064388E+04,D52=0.2322997491E+01,D53=-.2710031511E+00,
C    .D61=-.7265968242E+02,D62=0.3593588680E+03,D63=-.3933599161E+03,
C    .D71=0.9312252325E+01,D72=-.1380396384E+01,D73=0.7373026488E-01,
C    .D74=-.5070954743E+02,D75=0.7931122062E+00,D81=0.3599502276E+01,
C    .D82=-.5609868931E+00,D83=0.4980220561E+04,D84=-.1087179983E+01)
C     (LAST 4 PARAMETERS RENUMBERED)
C
C VBENDA VALUES FROM FIT634B (USED FOR SURF706)
C     PARAMETER(
C    .A11=-.2918252280E+03,A12=0.2164569141E+02,A13=-.4005497543E+00,
C    .A21=-.2890774947E+01,A22=0.1032032542E+02,A23=0.2681748056E+02,
C    .A24=0.2633751055E+02,A31=0.6180641351E+01,A32=0.5037667041E+01,
C    .A41=-.1125570079E+00,A42=-.3176529304E-02,A43=0.9068915355E+00,
C    .A44=-.7228418516E+00,A51=0.2785898232E+03,A52=0.4764442446E+02,
C    .A53=0.8621545662E+01)
C   VBENDB VALUES FROM FIT634B (USED FOR SURF706)
C     PARAMETER(
C    .G11=-.4241912309E+02,G12=0.2951365281E+01,G13=-.4840201514E-01,
C    .G21=-.1159168549E-03,G22=0.8688003567E-02,G23=-.3486923900E-02,
C    .G24=0.8312212839E-03,G31=0.5621810473E-01,G32=-.9776930747E-02,
C    .G41=-.1178456251E-01,G42=0.3491086729E-02,G43=0.7430516993E-01,
C    .G44=-.9643636957E-01,G51=0.4735533782E+02,G52=0.3001038808E+01,
C    .G53=0.1896630453E+01)
C
C PARAMETERS REQUIRED FOR CBEND CALCULATIONS:
C
C-- CBENDA VALUES FROM FIT706D.OUT
C     PARAMETER(
C    .C11=0.7107064647E+04,C12=-.3728421267E+04,C13=0.1757654624E+03,
C    .C14=-.9725998132E+04,C15=0.4665086074E+04,C21=-.4435165986E+02,
C    .C22=0.1604477309E+03,C23=0.5805142925E+03,C24=0.6892349445E+03,
C    .C31=0.6581730442E+03,C32=0.6078389739E+02,C41=0.2885182566E+01,
C    .C42=-.1728169916E+01,C43=-.1119535503E+03,C44=0.4052536250E+02,
C    .C51=0.2540505673E+03,C52=-.5762083627E+03,C53=0.1295901320E+03,
C    .C61=-.2131706075E+04,C62=0.9084452020E+04,C63=-.1138253963E+05,
C    .C71=-.3964298833E+02,C72=-.5019979693E+01,C73=0.2906541488E+00,
C    .C74=0.1212943686E+04,C75=0.4140463415E+02,C81=0.1752855549E+03,
C    .C82=-.2496320107E+02,C83=0.3765413052E+03,C84=-.5480488130E+02)
C-- CBENDB VALUES FROM FIT706D.OUT
C     PARAMETER(
C    .D11=0.1917166552E+04,D12=-.6542563392E+03,D13=0.6793758367E+02,
C    .D14=-.1694968847E+04,D15=0.6866649703E+03,D21=-.2137567948E+00,
C    .D22=0.4975938228E-01,D23=0.9364998295E-01,D24=-.2444320779E-01,
C    .D31=-.2863126914E+02,D32=0.5443219625E-01,D41=0.9673956120E+00,
C    .D42=-.1160159706E+01,D43=-.2424199759E+02,D44=0.8569424490E+01,
C    .D51=-.6517635862E+04,D52=0.1518098147E+03,D53=-.2706514366E+02,
C    .D61=0.4308392956E+02,D62=-.1234851732E+03,D63=0.2320626055E+03,
C    .D71=0.1049541418E+02,D72=-.2424169341E+01,D73=0.1745646946E+00,
C    .D74=0.2603615561E+02,D75=0.1345799970E+02,D81=-.6653710513E+01,
C    .D82=-.2576854447E+00,D83=0.6172608425E+04,D84=-.1328142473E+02)
C END OF PARAMETERS
C
C   FEB27/91 NEW C51 = C51+C83;   NEW D51 = D51+D83
      CX1 = C51 + C83
      DX1 = D51 + D83
      R1 = RPASS(1)
      R2 = RPASS(2)
      R3 = RPASS(3)
      T1 = R1*R1 -R2*R2 -R3*R3
      T2 = R2*R2 -R3*R3 -R1*R1
      T3 = R3*R3 -R1*R1 -R2*R2
C  CALCULATE THE COSINES OF THE THREE INTERNAL ANGLES:
      C1 = T1 / (-2.D0*R2*R3)
      C2 = T2 / (-2.D0*R3*R1)
      C3 = T3 / (-2.D0*R1*R2)
      SUM = C1 + C2 + C3
      B1A = 1.D0 - SUM
      C1CUBE = C1*C1*C1
      C2CUBE = C2*C2*C2
      C3CUBE = C3*C3*C3
      COS3T1 = 4.D0*C1CUBE - 3.D0*C1
      COS3T2 = 4.D0*C2CUBE - 3.D0*C2
      COS3T3 = 4.D0*C3CUBE - 3.D0*C3
      SUMB   = COS3T1 + COS3T2 + COS3T3
      B1B    = 1.D0 - ( Z58*SUMB + Z38*SUM )
C  CALCULATE DERIVATIVES IF DESIRED
         DC1DR1 = -R1/(R2*R3)
         DC2DR2 = -R2/(R1*R3)
         DC3DR3 = -R3/(R1*R2)
         DC1DR2 = ( T1/(R2*R2) + 2.D0 )/(2.D0*R3)
         DC1DR3 = ( T1/(R3*R3) + 2.D0 )/(2.D0*R2)
         DC2DR1 = ( T2/(R1*R1) + 2.D0 )/(2.D0*R3)
         DC2DR3 = ( T2/(R3*R3) + 2.D0 )/(2.D0*R1)
         DC3DR1 = ( T3/(R1*R1) + 2.D0 )/(2.D0*R2)
         DC3DR2 = ( T3/(R2*R2) + 2.D0 )/(2.D0*R1)
         DB1A(1) = -1.D0*( DC1DR1 + DC2DR1 + DC3DR1 )
         DB1A(2) = -1.D0*( DC1DR2 + DC2DR2 + DC3DR2 )
         DB1A(3) = -1.D0*( DC1DR3 + DC2DR3 + DC3DR3 )
         D1 = 12.D0*C1*C1 - 3.D0
         D2 = 12.D0*C2*C2 - 3.D0
         D3 = 12.D0*C3*C3 - 3.D0
         DB1B(1) = -Z58*( D1*DC1DR1 + D2*DC2DR1 + D3*DC3DR1 )
     .             -Z38*(    DC1DR1 +    DC2DR1 +    DC3DR1 )
         DB1B(2) = -Z58*( D1*DC1DR2 + D2*DC2DR2 + D3*DC3DR2 )
     .             -Z38*(    DC1DR2 +    DC2DR2 +    DC3DR2 )
         DB1B(3) = -Z58*( D1*DC1DR3 + D2*DC2DR3 + D3*DC3DR3 )
     .             -Z38*(    DC1DR3 +    DC2DR3 +    DC3DR3 )
CD    IF(IPR.GT.1)THEN
CD        WRITE(7,*)
CD        WRITE(7,*) ' FROM SUBR.VBCB -----------------------------'
CD        WRITE(7,6000) RPASS
CD        WRITE(7,6010) B1A,B1B
CD        WRITE(7,6020) DB1A,DB1B
CD        IF(IPR.GT.2)THEN
CD           WRITE(7,7000) DC1DR1,DC1DR2,DC1DR3,
CD   .                     DC2DR1,DC2DR2,DC2DR3,
CD   .                     DC3DR1,DC3DR2,DC3DR3
CD           WRITE(7,7010) T1,T2,T3,C1,C2,C3,D1,D2,D3
CD        END IF
CD    END IF
C
C  CALCULATE THE QUANTITIES USED BY BOTH VBEND AND CBEND:
C
      R    = R1 + R2 + R3
      RSQ  = R*R
      B2   = 1.D0/R1 + 1.D0/R2 + 1.D0/R3
      B3   = (R2-R1)*(R2-R1) + (R3-R2)*(R3-R2) + (R1-R3)*(R1-R3)
      EPS2 = 1.D-12
      B3B  = DSQRT( B3 + EPS2 )
         DB2(1) = -1.D0/(R1*R1)
         DB2(2) = -1.D0/(R2*R2)
         DB2(3) = -1.D0/(R3*R3)
         DB3(1) = 4.D0*R1 -2.D0*R2 -2.D0*R3
         DB3(2) = 4.D0*R2 -2.D0*R3 -2.D0*R1
         DB3(3) = 4.D0*R3 -2.D0*R1 -2.D0*R2
         DB3(4) = 0.5D0*DB3(1)/B3B
         DB3(5) = 0.5D0*DB3(2)/B3B
         DB3(6) = 0.5D0*DB3(3)/B3B
         EXP1   = DEXP(-BETA1*R)
         EXP2   = DEXP(-BETA2*RSQ)
         EXP7   = DEXP(-BETA3*R)
         DEXP1  = -BETA1*EXP1
         DEXP2  = -2.D0*BETA2*R*EXP2
         DEXP7  = -BETA3*EXP7
C   DO THE VBNDA CALCULATIONS:
         B1    = B1A
         B12   = B1*B1
         B13   = B12*B1
         B14   = B13*B1
         B15   = B14*B1
         ASUM  = A11 +  A12*R +  A13*RSQ
         BSUM  = A21*B12 +  A22*B13 +  A23*B14 +  A24*B15
         CSUM  = A31*B1*EXP1 +  A32*B12*EXP2
         DSUM1 = A41*EXP1 +  A42*EXP2
         DSUM2 = A43*EXP1 +  A44*EXP2
         FSUM  = A51 +  A52*R +  A53*RSQ
         VB(1,1) = B1*ASUM*EXP1
         VB(1,2) = BSUM * EXP2
         VB(1,3) = B2*CSUM
         VB(1,4) = B1*B3*DSUM1 + B1*B3B*DSUM2
         VB(1,5) = B1* FSUM  *EXP7
         VBND(1) = VB(1,1)+VB(1,2)+VB(1,3)+VB(1,4)+VB(1,5)
C
         DASUM   = A12+2.D0*A13*R
         DBSUM   =  2.D0* A21*B1  + 3.D0* A22*B12
     .            + 4.D0* A23*B13 + 5.D0* A24*B14
         DDSUM1  =  A41*DEXP1 +  A42*DEXP2
         DDSUM2  =  A43*DEXP1 +  A44*DEXP2
         DFSUM   =  A52 +2.D0*A53*R
         DVB1(1) = DB1A(1)*ASUM*EXP1 + B1*DASUM*EXP1 +B1*ASUM*DEXP1
         DVB1(2) = DB1A(2)*ASUM*EXP1 + B1*DASUM*EXP1 +B1*ASUM*DEXP1
         DVB1(3) = DB1A(3)*ASUM*EXP1 + B1*DASUM*EXP1 +B1*ASUM*DEXP1
C
         DVB2(1) = DBSUM*DB1A(1)*EXP2 + BSUM*DEXP2
         DVB2(2) = DBSUM*DB1A(2)*EXP2 + BSUM*DEXP2
         DVB2(3) = DBSUM*DB1A(3)*EXP2 + BSUM*DEXP2
C
C CALCULATE THE VB3 DERIVATIVES:
         DVB3(1)=DB2(1)*CSUM + B2*(  A31*DB1A(1)*EXP1 +  A31*B1*DEXP1
     .                  +2.D0* A32*B1*DB1A(1)*EXP2 +  A32*B12*DEXP2)
         DVB3(2)=DB2(2)*CSUM + B2*(  A31*DB1A(2)*EXP1 +  A31*B1*DEXP1
     .                  +2.D0* A32*B1*DB1A(2)*EXP2 +  A32*B12*DEXP2)
         DVB3(3)=DB2(3)*CSUM + B2*(  A31*DB1A(3)*EXP1 +  A31*B1*DEXP1
     .                  +2.D0* A32*B1*DB1A(3)*EXP2 +  A32*B12*DEXP2)
C
C CALCULATE THE VB4 DERIVATIVES (MAY 27/90):
         DVB4(1)= DB1A(1)*B3 *DSUM1+  B1*DB3(1)*DSUM1 + B1*B3 *DDSUM1
     .          + DB1A(1)*B3B*DSUM2 + B1*DB3(4)*DSUM2 + B1*B3B*DDSUM2
         DVB4(2)= DB1A(2)*B3 *DSUM1+  B1*DB3(2)*DSUM1 + B1*B3 *DDSUM1
     .          + DB1A(2)*B3B*DSUM2 + B1*DB3(5)*DSUM2 + B1*B3B*DDSUM2
         DVB4(3)= DB1A(3)*B3 *DSUM1+  B1*DB3(3)*DSUM1 + B1*B3 *DDSUM1
     .          + DB1A(3)*B3B*DSUM2 + B1*DB3(6)*DSUM2 + B1*B3B*DDSUM2
C
C CALCULATE THE VB5 DERIVATIVES:
         DVB5(1) = DB1A(1)*FSUM*EXP7 + B1*DFSUM*EXP7 + B1*FSUM*DEXP7
         DVB5(2) = DB1A(2)*FSUM*EXP7 + B1*DFSUM*EXP7 + B1*FSUM*DEXP7
         DVB5(3) = DB1A(3)*FSUM*EXP7 + B1*DFSUM*EXP7 + B1*FSUM*DEXP7
C
C CALCULATE THE OVERALL DERIVATIVES:
         DVBNDA(1) = DVB1(1)+DVB2(1)+DVB3(1)+DVB4(1)+DVB5(1)
         DVBNDA(2) = DVB1(2)+DVB2(2)+DVB3(2)+DVB4(2)+DVB5(2)
         DVBNDA(3) = DVB1(3)+DVB2(3)+DVB3(3)+DVB4(3)+DVB5(3)
         VBNDA = VBND(1)
C
C---- NOW DO THE VBENDB CALCULATIONS --------
         B1    = B1B
         B12   = B1*B1
         B13   = B12*B1
         B14   = B13*B1
         B15   = B14*B1
         ASUM  = G11 +  G12*R +  G13*RSQ
         BSUM  = G21*B12 +  G22*B13 +  G23*B14 +  G24*B15
         CSUM  = G31*B1*EXP1 +  G32*B12*EXP2
         DSUM1 = G41*EXP1 +  G42*EXP2
         DSUM2 = G43*EXP1 +  G44*EXP2
         FSUM  = G51 +  G52*R +  G53*RSQ
         VB(2,1) = B1*ASUM*EXP1
         VB(2,2) = BSUM * EXP2
         VB(2,3) = B2*CSUM
         VB(2,4) = B1*B3*DSUM1 + B1*B3B*DSUM2
         VB(2,5) = B1* FSUM  *EXP7
         VBND(2) = VB(2,1)+VB(2,2)+VB(2,3)+VB(2,4)+VB(2,5)
C
         DASUM   = G12+2.D0*G13*R
         DBSUM   =  2.D0* G21*B1  + 3.D0* G22*B12
     .             +4.D0* G23*B13 + 5.D0* G24*B14
         DDSUM1  =  G41*DEXP1 +  G42*DEXP2
         DDSUM2  =  G43*DEXP1 +  G44*DEXP2
         DFSUM   =  G52 +2.D0*G53*R
         DVB1(1) = DB1B(1)*ASUM*EXP1 + B1*DASUM*EXP1 +B1*ASUM*DEXP1
         DVB1(2) = DB1B(2)*ASUM*EXP1 + B1*DASUM*EXP1 +B1*ASUM*DEXP1
         DVB1(3) = DB1B(3)*ASUM*EXP1 + B1*DASUM*EXP1 +B1*ASUM*DEXP1
C
         DVB2(1) = DBSUM*DB1B(1)*EXP2 + BSUM*DEXP2
         DVB2(2) = DBSUM*DB1B(2)*EXP2 + BSUM*DEXP2
         DVB2(3) = DBSUM*DB1B(3)*EXP2 + BSUM*DEXP2
C
C CALCULATE THE VB3 DERIVATIVES:
         DVB3(1)=DB2(1)*CSUM + B2*(  G31*DB1B(1)*EXP1 +  G31*B1*DEXP1
     .                  +2.D0* G32*B1*DB1B(1)*EXP2 +  G32*B12*DEXP2)
         DVB3(2)=DB2(2)*CSUM + B2*(  G31*DB1B(2)*EXP1 +  G31*B1*DEXP1
     .                  +2.D0* G32*B1*DB1B(2)*EXP2 +  G32*B12*DEXP2)
         DVB3(3)=DB2(3)*CSUM + B2*(  G31*DB1B(3)*EXP1 +  G31*B1*DEXP1
     .                  +2.D0* G32*B1*DB1B(3)*EXP2 +  G32*B12*DEXP2)
C
C CALCULATE THE VB4 DERIVATIVES (MAY 27/90):
         DVB4(1)= DB1B(1)*B3 *DSUM1+  B1*DB3(1)*DSUM1 + B1*B3 *DDSUM1
     .          + DB1B(1)*B3B*DSUM2 + B1*DB3(4)*DSUM2 + B1*B3B*DDSUM2
         DVB4(2)= DB1B(2)*B3 *DSUM1+  B1*DB3(2)*DSUM1 + B1*B3 *DDSUM1
     .          + DB1B(2)*B3B*DSUM2 + B1*DB3(5)*DSUM2 + B1*B3B*DDSUM2
         DVB4(3)= DB1B(3)*B3 *DSUM1+  B1*DB3(3)*DSUM1 + B1*B3 *DDSUM1
     .          + DB1B(3)*B3B*DSUM2 + B1*DB3(6)*DSUM2 + B1*B3B*DDSUM2
C
C CALCULATE THE VB5 DERIVATIVES:
         DVB5(1) = DB1B(1)*FSUM*EXP7 + B1*DFSUM*EXP7 + B1*FSUM*DEXP7
         DVB5(2) = DB1B(2)*FSUM*EXP7 + B1*DFSUM*EXP7 + B1*FSUM*DEXP7
         DVB5(3) = DB1B(3)*FSUM*EXP7 + B1*DFSUM*EXP7 + B1*FSUM*DEXP7
C
C CALCULATE THE OVERALL DERIVATIVES:
         DVBNDB(1) = DVB1(1)+DVB2(1)+DVB3(1)+DVB4(1)+DVB5(1)
         DVBNDB(2) = DVB1(2)+DVB2(2)+DVB3(2)+DVB4(2)+DVB5(2)
         DVBNDB(3) = DVB1(3)+DVB2(3)+DVB3(3)+DVB4(3)+DVB5(3)
         VBNDB = VBND(2)
C
C  NOW CALCULATE THE CBEND CORRECTION TERMS:
C  CBNDA USES THE B1A FORMULA AND CBNDB USES THE B1B FORMULA
      IF(ICOMPC.EQ.0) RETURN
      SUMT  = T(1) + T(2) + T(3)
      RCU   = RSQ*R
      P     = R1*R2*R3
      PSQ   = P*P
      PCU   = PSQ*P
      DP(1) = R2*R3
      DP(2) = R3*R1
      DP(3) = R1*R2
      CBNDS(1)  = 0.D0
      CBNDS(2)  = 0.D0
      DCBNDA(1) = 0.D0
      DCBNDA(2) = 0.D0
      DCBNDA(3) = 0.D0
      DCBNDB(1) = 0.D0
      DCBNDB(2) = 0.D0
      DCBNDB(3) = 0.D0
C  CALCULATE EXPONENTIALS AND DERIVATIVES (COMMON TO CBNDA AND CBNDB):
            EXP7  = DEXP( -BETA2*RCU )
            DEXP7 = -3.D0*BETA2*RSQ*EXP7
C ----- CALCULATE THE CBNDA CORRECTION TERM ----------------
       B1    = B1A
       B12   = B1*B1
       B13   = B12*B1
       B14   = B13*B1
       B15   = B14*B1
       ASUM  = C11 + C12*R + C13*RSQ + C14/R + C15/RSQ
       BSUM  = C21*B12 + C22*B13 + C23*B14 + C24*B15
       CSUM  = C31*B1 *EXP1  + C32*B12*EXP2
       DSUM1 = C41*EXP1 + C42*EXP2
       DSUM2 = C43*EXP1 + C44*EXP2
       FSUM  = CX1 + C52*R + C53*RSQ
       GSUM  = C61 + C62/R + C63/RSQ
       AASUM = C71 + C72*P + C73*PSQ + C74/P + C75/PSQ
       FFSUM =       C81*P + C82*PSQ + C84/PSQ
       DASUM = C12 + 2.D0*C13*R -C14/RSQ -2.D0*C15/RCU
       DBSUM = 2.D0*C21*B1 +3.D0*C22*B12 +4.D0*C23*B13 +5.D0*C24*B14
       DDSUM1= C41*DEXP1 + C42*DEXP2
       DDSUM2= C43*DEXP1 + C44*DEXP2
       DFSUM = C52 + 2.D0*C53*R
       DGSUM = -C62/RSQ - 2.D0*C63/RCU
       DAASUM = C72 +2.D0*C73*P -C74/PSQ -2.D0*C75/PCU
       DFFSUM = C81 + 2.D0*C82*P -2.D0*C84/PCU
         CB(1,1)  = B1*ASUM * EXP1  /P
         CB(1,2)  = BSUM * EXP2
         CB(1,3)  = B2*CSUM
         CB(1,4)  = B1*B3*DSUM1 + B1*B3B*DSUM2
         CB(1,5)  = B1 * FSUM * EXP7 /P
         CB(1,6)  = B1 * GSUM * EXP7
         CB(1,7)  = B1 * AASUM * EXP2
         CB(1,8)  = B1 * FFSUM * EXP7
         CBNDS(1) = CB(1,1)+CB(1,2)+CB(1,3)+CB(1,4)+CB(1,5)
     .                     +CB(1,6)+CB(1,7)+CB(1,8)
C       CALCULATE THE DERIVATIVES:
         DCB1(1) = DB1A(1)*ASUM*EXP1/P + B1*DASUM*EXP1/P
     .                +B1*ASUM*DEXP1/P - B1*ASUM*EXP1*DP(1)/PSQ
         DCB1(2) = DB1A(2)*ASUM*EXP1/P + B1*DASUM*EXP1/P
     .                +B1*ASUM*DEXP1/P - B1*ASUM*EXP1*DP(2)/PSQ
         DCB1(3) = DB1A(3)*ASUM*EXP1/P + B1*DASUM*EXP1/P
     .                +B1*ASUM*DEXP1/P - B1*ASUM*EXP1*DP(3)/PSQ
C
         DCB2(1) = DBSUM*DB1A(1)*EXP2 + BSUM*DEXP2
         DCB2(2) = DBSUM*DB1A(2)*EXP2 + BSUM*DEXP2
         DCB2(3) = DBSUM*DB1A(3)*EXP2 + BSUM*DEXP2
C
         DCB3(1) = DB2(1)*CSUM + B2*( C31*DB1A(1)*EXP1 +C31*B1*DEXP1
     .            +C32 *2.D0*B1*DB1A(1)*EXP2 + C32 *B12*DEXP2 )
         DCB3(2) = DB2(2)*CSUM + B2*( C31 *DB1A(2)*EXP1 +C31 *B1*DEXP1
     .            +C32 *2.D0*B1*DB1A(2)*EXP2 + C32 *B12*DEXP2 )
         DCB3(3) = DB2(3)*CSUM + B2*( C31 *DB1A(3)*EXP1 +C31 *B1*DEXP1
     .            +C32 *2.D0*B1*DB1A(3)*EXP2 + C32 *B12*DEXP2 )
C          DCB4 EQUATIONS MAY 27/90:
         DCB4(1) = DB1A(1)*B3 *DSUM1 + B1*DB3(1)*DSUM1 + B1*B3 *DDSUM1
     .           + DB1A(1)*B3B*DSUM2 + B1*DB3(4)*DSUM2 + B1*B3B*DDSUM2
         DCB4(2) = DB1A(2)*B3 *DSUM1 + B1*DB3(2)*DSUM1 + B1*B3 *DDSUM1
     .           + DB1A(2)*B3B*DSUM2 + B1*DB3(5)*DSUM2 + B1*B3B*DDSUM2
         DCB4(3) = DB1A(3)*B3 *DSUM1 + B1*DB3(3)*DSUM1 + B1*B3 *DDSUM1
     .           + DB1A(3)*B3B*DSUM2 + B1*DB3(6)*DSUM2 + B1*B3B*DDSUM2
C          DCB5 EQUATIONS MAY 27/90: (CORRECTED JUN01)
         DCB5(1) = DB1A(1)*FSUM*EXP7/P + B1*DFSUM*EXP7/P
     .               + B1*FSUM*DEXP7/P - B1*FSUM*EXP7*DP(1)/PSQ
         DCB5(2) = DB1A(2)*FSUM*EXP7/P + B1*DFSUM*EXP7/P
     .               + B1*FSUM*DEXP7/P - B1*FSUM*EXP7*DP(2)/PSQ
         DCB5(3) = DB1A(3)*FSUM*EXP7/P + B1*DFSUM*EXP7/P
     .               + B1*FSUM*DEXP7/P - B1*FSUM*EXP7*DP(3)/PSQ
C
         DCB6(1) = DB1A(1)*GSUM*EXP7 + B1*DGSUM*EXP7 + B1*GSUM*DEXP7
         DCB6(2) = DB1A(2)*GSUM*EXP7 + B1*DGSUM*EXP7 + B1*GSUM*DEXP7
         DCB6(3) = DB1A(3)*GSUM*EXP7 + B1*DGSUM*EXP7 + B1*GSUM*DEXP7
C
         DCB7(1) = DB1A(1)*AASUM*EXP2 + B1*DAASUM*DP(1)*EXP2
     .                               + B1* AASUM*DEXP2
         DCB7(2) = DB1A(2)*AASUM*EXP2 + B1*DAASUM*DP(2)*EXP2
     .                               + B1* AASUM*DEXP2
         DCB7(3) = DB1A(3)*AASUM*EXP2 + B1*DAASUM*DP(3)*EXP2
     .                               + B1* AASUM*DEXP2
         DCB8(1) = DB1A(1)*FFSUM*EXP7  +B1*DFFSUM*DP(1)*EXP7
     .                                +B1* FFSUM*DEXP7
         DCB8(2) = DB1A(2)*FFSUM*EXP7  +B1*DFFSUM*DP(2)*EXP7
     .                                +B1* FFSUM*DEXP7
         DCB8(3) = DB1A(3)*FFSUM*EXP7  +B1*DFFSUM*DP(3)*EXP7
     .                                +B1* FFSUM*DEXP7
           DCBNDS(1,1) = DCB1(1) +DCB2(1) +DCB3(1) +DCB4(1) +DCB5(1)
     .                  +DCB6(1) +DCB7(1) +DCB8(1)
           DCBNDS(1,2) = DCB1(2) +DCB2(2) +DCB3(2) +DCB4(2) +DCB5(2)
     .                  +DCB6(2) +DCB7(2) +DCB8(2)
           DCBNDS(1,3) = DCB1(3) +DCB2(3) +DCB3(3) +DCB4(3) +DCB5(3)
     .                  +DCB6(3) +DCB7(3) +DCB8(3)
C
C ----- CALCULATE THE CBNDB CORRECTION TERM ----------------
       B1    = B1B
       B12   = B1*B1
       B13   = B12*B1
       B14   = B13*B1
       B15   = B14*B1
       ASUM  = D11 + D12*R + D13*RSQ + D14/R + D15/RSQ
       BSUM  = D21*B12 + D22*B13 + D23*B14 + D24*B15
       CSUM  = D31*B1 *EXP1  + D32*B12*EXP2
       DSUM1 = D41*EXP1 + D42*EXP2
       DSUM2 = D43*EXP1 + D44*EXP2
       FSUM  = DX1 + D52*R + D53*RSQ
       GSUM  = D61 + D62/R + D63/RSQ
       AASUM = D71 + D72*P + D73*PSQ + D74/P + D75/PSQ
       FFSUM =       D81*P + D82*PSQ + D84/PSQ
       DASUM = D12 + 2.D0*D13*R -D14/RSQ -2.D0*D15/RCU
       DBSUM = 2.D0*D21*B1 +3.D0*D22*B12 +4.D0*D23*B13 +5.D0*D24*B14
       DDSUM1= D41*DEXP1 + D42*DEXP2
       DDSUM2= D43*DEXP1 + D44*DEXP2
       DFSUM = D52 + 2.D0*D53*R
       DGSUM = -D62/RSQ - 2.D0*D63/RCU
       DAASUM = D72 +2.D0*D73*P -D74/PSQ -2.D0*D75/PCU
       DFFSUM = D81 + 2.D0*D82*P -2.D0*D84/PCU
         CB(2,1)  = B1*ASUM * EXP1  /P
         CB(2,2)  = BSUM * EXP2
         CB(2,3)  = B2*CSUM
         CB(2,4)  = B1*B3*DSUM1 + B1*B3B*DSUM2
         CB(2,5)  = B1 * FSUM * EXP7 /P
         CB(2,6)  = B1 * GSUM * EXP7
         CB(2,7)  = B1 * AASUM * EXP2
         CB(2,8)  = B1 * FFSUM * EXP7
         CBNDS(2) = CB(2,1)+CB(2,2)+CB(2,3)+CB(2,4)+CB(2,5)
     .                     +CB(2,6)+CB(2,7)+CB(2,8)
C       CALCULATE THE DERIVATIVES:
         DCB1(1) = DB1B(1)*ASUM*EXP1/P + B1*DASUM*EXP1/P
     .                +B1*ASUM*DEXP1/P - B1*ASUM*EXP1*DP(1)/PSQ
         DCB1(2) = DB1B(2)*ASUM*EXP1/P + B1*DASUM*EXP1/P
     .                +B1*ASUM*DEXP1/P - B1*ASUM*EXP1*DP(2)/PSQ
         DCB1(3) = DB1B(3)*ASUM*EXP1/P + B1*DASUM*EXP1/P
     .                +B1*ASUM*DEXP1/P - B1*ASUM*EXP1*DP(3)/PSQ
C
         DCB2(1) = DBSUM*DB1B(1)*EXP2 + BSUM*DEXP2
         DCB2(2) = DBSUM*DB1B(2)*EXP2 + BSUM*DEXP2
         DCB2(3) = DBSUM*DB1B(3)*EXP2 + BSUM*DEXP2
C
         DCB3(1)= DB2(1)*CSUM + B2*( D31 *DB1B(1)*EXP1 +D31 *B1*DEXP1
     .           +D32 *2.D0*B1*DB1B(1)*EXP2 + D32 *B12*DEXP2 )
         DCB3(2)= DB2(2)*CSUM + B2*( D31 *DB1B(2)*EXP1 +D31 *B1*DEXP1
     .           +D32 *2.D0*B1*DB1B(2)*EXP2 + D32 *B12*DEXP2 )
         DCB3(3)= DB2(3)*CSUM + B2*( D31 *DB1B(3)*EXP1 +D31 *B1*DEXP1
     .           +D32 *2.D0*B1*DB1B(3)*EXP2 + D32 *B12*DEXP2 )
C          DCB4 EQUATIONS MAY 27/90:
         DCB4(1) = DB1B(1)*B3 *DSUM1 + B1*DB3(1)*DSUM1 + B1*B3 *DDSUM1
     .           + DB1B(1)*B3B*DSUM2 + B1*DB3(4)*DSUM2 + B1*B3B*DDSUM2
         DCB4(2) = DB1B(2)*B3 *DSUM1 + B1*DB3(2)*DSUM1 + B1*B3 *DDSUM1
     .           + DB1B(2)*B3B*DSUM2 + B1*DB3(5)*DSUM2 + B1*B3B*DDSUM2
         DCB4(3) = DB1B(3)*B3 *DSUM1 + B1*DB3(3)*DSUM1 + B1*B3 *DDSUM1
     .           + DB1B(3)*B3B*DSUM2 + B1*DB3(6)*DSUM2 + B1*B3B*DDSUM2
C          DCB5 EQUATIONS MAY 27/90:
         DCB5(1) = DB1B(1)*FSUM*EXP7/P + B1*DFSUM*EXP7/P
     .               + B1*FSUM*DEXP7/P - B1*FSUM*EXP7*DP(1)/PSQ
         DCB5(2) = DB1B(2)*FSUM*EXP7/P + B1*DFSUM*EXP7/P
     .               + B1*FSUM*DEXP7/P - B1*FSUM*EXP7*DP(2)/PSQ
         DCB5(3) = DB1B(3)*FSUM*EXP7/P + B1*DFSUM*EXP7/P
     .               + B1*FSUM*DEXP7/P - B1*FSUM*EXP7*DP(3)/PSQ
C
         DCB6(1) = DB1B(1)*GSUM*EXP7 + B1*DGSUM*EXP7 + B1*GSUM*DEXP7
         DCB6(2) = DB1B(2)*GSUM*EXP7 + B1*DGSUM*EXP7 + B1*GSUM*DEXP7
         DCB6(3) = DB1B(3)*GSUM*EXP7 + B1*DGSUM*EXP7 + B1*GSUM*DEXP7
C
         DCB7(1) = DB1B(1)*AASUM*EXP2 + B1*DAASUM*DP(1)*EXP2
     .                               + B1* AASUM*DEXP2
         DCB7(2) = DB1B(2)*AASUM*EXP2 + B1*DAASUM*DP(2)*EXP2
     .                               + B1* AASUM*DEXP2
         DCB7(3) = DB1B(3)*AASUM*EXP2 + B1*DAASUM*DP(3)*EXP2
     .                               + B1* AASUM*DEXP2
         DCB8(1) = DB1B(1)*FFSUM*EXP7  +B1*DFFSUM*DP(1)*EXP7
     .                                +B1* FFSUM*DEXP7
         DCB8(2) = DB1B(2)*FFSUM*EXP7  +B1*DFFSUM*DP(2)*EXP7
     .                                +B1* FFSUM*DEXP7
         DCB8(3) = DB1B(3)*FFSUM*EXP7  +B1*DFFSUM*DP(3)*EXP7
     .                                +B1* FFSUM*DEXP7
           DCBNDS(2,1) = DCB1(1) +DCB2(1) +DCB3(1) +DCB4(1) +DCB5(1)
     .                  +DCB6(1) +DCB7(1) +DCB8(1)
           DCBNDS(2,2) = DCB1(2) +DCB2(2) +DCB3(2) +DCB4(2) +DCB5(2)
     .                  +DCB6(2) +DCB7(2) +DCB8(2)
           DCBNDS(2,3) = DCB1(3) +DCB2(3) +DCB3(3) +DCB4(3) +DCB5(3)
     .                  +DCB6(3) +DCB7(3) +DCB8(3)
C
C  CALCULATE THE TOTAL DERIVATIVE FROM THE PIECES:
         DCBNDA(1) = DT(1) * CBNDS(1)    + SUMT * DCBNDS(1,1)
         DCBNDB(1) = DT(1) * CBNDS(2)    + SUMT * DCBNDS(2,1)
         DCBNDA(2) = DT(2) * CBNDS(1)    + SUMT * DCBNDS(1,2)
         DCBNDB(2) = DT(2) * CBNDS(2)    + SUMT * DCBNDS(2,2)
         DCBNDA(3) = DT(3) * CBNDS(1)    + SUMT * DCBNDS(1,3)
         DCBNDB(3) = DT(3) * CBNDS(2)    + SUMT * DCBNDS(2,3)
      CBNDA = SUMT*CBNDS(1)
      CBNDB = SUMT*CBNDS(2)
C  SOME DEBUGGING PRINT STATEMENTS:
CD    IF(IPR.GT.1)THEN
CD       WRITE(7,7650) B2,B3
CD       WRITE(7,7700) DB2,DB3
CD       WRITE(7,*)
CD       WRITE(7,*) 'VBEND VALUES :'
CD       K=1
CD       WRITE(7,7150) K,DVB1,DVB2,DVB3,DVB4,DVB5,DVBNDA
CD       WRITE(7,7155) (VB(1,I),I=1,5)
CD       K=2
CD       WRITE(7,7150) K,DVB1,DVB2,DVB3,DVB4,DVB5,DVBNDB
CD       WRITE(7,7155) (VB(2,I),I=1,5)
CD       WRITE(7,7100) VBNDA,VBNDB
CD       WRITE(7,*)
CD       WRITE(7,*) 'CBEND VALUES:'
CD       WRITE(7,*) 'C81,C81,C82=',C81,C81,C82
CD       WRITE(7,*) 'C91,C83,C84=',C91,C83,C84
CD       WRITE(7,*) 'D81,D81,D82=',D81,D81,D82
CD       WRITE(7,*) 'D91,D83,D84=',D91,D83,D84
CD       K=1
CD          WRITE(7,2020) K,(CB(1,J),J=1,9)
CD          WRITE(7,2033) SUMT,CBNDS(1)
CD       K=2
CD          WRITE(7,2020) K,(CB(K,J),J=1,9)
CD          WRITE(7,2033) SUMT,CBNDS(2)
CD       K=1
CD       WRITE(7,7152) K,DCB1,DCB2,DCB3,DCB4,DCB5,DCB6,
CD   .               DCB7,DCB8,DCB9,(DCBNDS(1,I),I=1,3)
CD       K=2
CD       WRITE(7,7152) K,DCB1,DCB2,DCB3,DCB4,DCB5,DCB6,
CD   .                 DCB7,DCB8,DCB9,(DCBNDS(K,I),I=1,3)
CD       WRITE(7,7200) DT,DCBNDA,DCBNDB
CD       WRITE(7,*) 'VBCB  EXIT -------------------------'
CD       WRITE(7,*)
CD    END IF
      RETURN
 2000 FORMAT(5X,'R1, R2, R3 = ',3(1X,F12.8))
 2010 FORMAT(5X,'B1A, B1B, B2, B3 = ',4(1X,G12.6))
 7100 FORMAT('    VBNDA = ',G12.6,'     VBNDB = ',G12.6)
 7650 FORMAT('    B2 =    ',G12.6,'       B3 =  ',G12.6)
 7700 FORMAT('    DB2 =   ',3(1X,G12.6),/,
     .       '    DB3 =   ',6(1X,G12.6))
 7150 FORMAT(' K=',I1,'     DERIVATIVES: ',/,
     .       ' DVB1  = ',3(1X,G16.10),/,
     .       ' DVB2  = ',3(1X,G16.10),/,
     .       ' DVB3  = ',3(1X,G16.10),/,
     .       ' DVB4  = ',3(1X,G16.10),/,
     .       ' DVB5  = ',3(1X,G16.10),/,
     .       ' DVBND = ',3(1X,G16.10))
 7155 FORMAT('  VB1 = ',G16.10,'  VB2 = ',G16.10,/,
     .       '  VB3 = ',G16.10,'  VB4 = ',G16.10,/,
     .       '  VB5 = ',G16.10)
 2020 FORMAT(5X,'K=',I1,' CB1,CB2,CB3= ',3(1X,G16.10),
     .     /,8X,' CB4,CB5,CB6= ',3(1X,G16.10),
     .     /,8X,' CB7,CB8,CB9= ',3(1X,G16.10))
 2033 FORMAT(5X,'SUMT = ',G12.6 ,'  CBNDS(K) = ',G12.6)
 7152 FORMAT('SUBR.CBEND:   K=',I1,'     DERIVATIVES: ',/,
     .       '  DCB1  = ',3(1X,G18.12),/,
     .       '  DCB2  = ',3(1X,G18.12),/,
     .       '  DCB3  = ',3(1X,G18.12),/,
     .       '  DCB4  = ',3(1X,G18.12),/,
     .       '  DCB5  = ',3(1X,G18.12),/,
     .       '  DCB6  = ',3(1X,G18.12),/,
     .       '  DCB7  = ',3(1X,G18.12),/,
     .       '  DCB8  = ',3(1X,G18.12),/,
     .       '  DCB9  = ',3(1X,G18.12),/,
     .       '  DCBND = ',3(1X,G18.12))
 7200 FORMAT('  DT    = ',3(1X,G18.12),/,
     .       '  DCBNDA= ',3(1X,G18.12),/,
     .       '  DCBNDB= ',3(1X,G18.12))
 6000 FORMAT(' FROM SUBR.B1AB:  R=',3(1X,F9.6))
 6010 FORMAT('    B1A = ',F16.10,'     B1B = ',F16.10)
 6020 FORMAT('    DB1A: ',3(1X,E12.6),/,'    DB1B: ',3(1X,E12.6))
 7000 FORMAT('    DC(J)/DR(I) DERIVATIVES: ',/,
     .    8X,3(1X,E12.6),/,8X,3(1X,E12.6),/,8X,3(1X,E12.6))
 7010 FORMAT('    RI2-RJ2-RK2: ',3(1X,E12.6),/,
     .       '       COSINES : ',3(1X,E12.6),/,
     .       '    12C**2 - 3 : ',3(1X,E12.6))
      END
C
      SUBROUTINE CHGEOM(R,IVALID)
C----------------------------------------------------------------------C
C APR21/95
C FIRST, CHECK THAT IT IS A VALID H3 GEOMETRY, WITHIN TOLERANCE "DERR":
C THIS TEST CORRECTED ON NOV23/93
C     GENERALLY, RLO + RMED > RHI
C     SOMETIMES, RLO + RMED = RHI (LINEAR GEOMETRY)
C     BUT IF     RLO + RMED < RHI THERE'S A PROBLEM
C
C ERROR1 ... CHECK THAT RLO+RMID>RHI
C ISTOP1 ... [0] WARNING [1]STOP CALCULATIONS
C
C ERROR2 ... CHECK THAT ALL R(I)'S ARE > 0.2 BOHR
C ISTOP2 ... [0] WARNING [1]STOP CALCULATIONS
C----------------------------------------------------------------------C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(DERR = 1.D-5)
      PARAMETER( ISTOP1 = 1 , ISTOP2 = 1 )
      DIMENSION R(3)
c       print*,r
      IVALID = 1
      RHI = DMAX1(R(1),R(2))
      IF(R(3).GE.RHI)THEN
         RMID = RHI
         RHI = R(3)
      ELSE
         RMID = DMAX1( R(3) , DMIN1(R(1),R(2)) )
      ENDIF
      RLO = DMIN1(R(1),R(2),R(3))
CX    IF(RLO+RMID.GT.RHI+DERR)THEN  {OLD CONDITION}
      IF(RLO+RMID+DERR.LT.RHI)THEN
         WRITE(6,*) ' RLONGEST          = ',RHI+DERR
         WRITE(6,*) ' RMIDDLE,RSHORTEST = ',RMID,RLO
         WRITE(6,*) 'WARNING: INVALID GEOMETRY R',R
c         IF(ISTOP1.GT.0) STOP ' STOP -- GEOMETRY ERROR '
C         IVALID = -1
      END IF
      IF(RLO.LT.0.2D0)THEN
CD       WRITE(6,*) ' SHORTEST DISTANCE: ',RLO
         WRITE(6,*) 'WARNING: INVALID GEOMETRY RLO',R
c         IF(ISTOP2.GT.0) STOP ' STOP -- GEOMETRY ERROR '
C         IVALID = -2
      END IF
      RETURN
      END

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense


text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

program id_set_def_stdin

use     types_mod, only : r8, i8
use  location_mod, only : location_type
use utilities_mod, only : get_unit
use  obs_kind_mod, only : QTY_SURFACE_PRESSURE
use     model_mod, only : static_init_model, get_model_size, &
                          get_state_meta_data

implicit none

! version controlled file description for error handling, do not edit
character(len=256), parameter :: source   = &
   "$URL$"
character(len=32 ), parameter :: revision = "$Revision$"
character(len=128), parameter :: revdate  = "$Date$"

integer :: i, model_size, var_type, iunit
type(location_type) :: location

! Write to file
iunit  = get_unit()
open(unit = iunit, file = 'id_set_def_stdin.out')

! Get the model size
call static_init_model()
model_size = get_model_size()

! Set the number of state variables, all observed
write(iunit, *) model_size

! No values or qc
write(iunit, *) 0
write(iunit, *) 0

! Loop through all the state variables, set the obs variance accordingly
do i = 1, model_size
   ! There are more obs
   write(iunit, *) 0

   ! Identity obs
   write(iunit, *) -1 * i

   ! Time is 0 days 0 seconds for create obs sequence
   write(iunit, *) 0, 0

   call get_state_meta_data(int(i,i8), location, var_type)

   ! Output the appropriate observational error variance
   if(var_type == QTY_SURFACE_PRESSURE) then
   write(iunit, *)  10000.0_r8
   else
      write(iunit, *) 1.0_r8
   endif

end do

! Output the default set_def.out file name
write(iunit, *) 'set_def.out'

end program id_set_def_stdin

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
!! This tests for issue ROSE-2567
!
!
module memory_pool_subs
! contains memory pool-related routines

contains

  subroutine free_memory(cptr)
    use iso_c_binding
    implicit none

    ! incoming

    type(c_ptr)               :: cptr

  end subroutine free_memory

end module memory_pool_subs

max_stars_repo_name: ouankou/rose

max_issues_repo_name: WildeGeist/rose

max_forks_repo_name: WildeGeist/rose


text:
!   Copyright 2014 College of William and Mary
!
!   Licensed under the Apache License, Version 2.0 (the "License");
!   you may not use this file except in compliance with the License.
!   You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
!   Unless required by applicable law or agreed to in writing, software
!   distributed under the License is distributed on an "AS IS" BASIS,
!   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!   See the License for the specific language governing permissions and
!   limitations under the License.

! Generate *_nu.nc from gridded data (nc file)
! The section on nc read needs to be modified as appropriate- search for
! 'start11' and 'end11'
! Beware the order of vertical levels in the nc file!!!
! Assume elev=0 in interpolation of 3D variables
! The code will do all it can to extrapolate: below bottom/above surface.
! If a pt in hgrid.ll is outside the background nc grid, const. values will be filled (from gen_hot_3Dth_from_nc.in)
!
! ifort -O2 -mcmodel=medium -CB -o gen_nu ~yinglong/git/schism/src/Utility/UtilLib/compute_zcor.f90 gen_nu.f90 -I$NETCDF/include -I$NETCDF_FORTRAN/include -L$NETCDF_FORTRAN/lib -L$NETCDF/lib -lnetcdf -lnetcdff

!   Input:
!     (1) hgrid.gr3;
!     (2) hgrid.ll;
!     (3) vgrid.in (SELFE R1703 and up);
!     (4) include.gr3: if depth=0, skip the interpolation to speed up.
!                      Should be larger than the non-0 regions in *_nudge.gr3
!     (5) gen_nu.in:
!                     1st line: T,S values for pts outside bg grid in nc
!                     2nd line: time step in .nc in sec; output stride
!                     3rd line: start yr, month, day (assume to be in PST)
!                     4th line: # of days
!     (6) ROMS nc files. The extent the ROMS files cover needs to be
!                      larger than the region specified in include.gr3
!   Output: TEM_nu.in
!   Debug outputs: fort.11 (fatal errors); fort.2[0-9], fort.9[5-9], fort.100

      program gen_hot
!      use typeSizes
      use netcdf
      use compute_zcor

!      implicit none

!      integer, parameter :: debug=1
      integer, parameter :: mnp=300000
      integer, parameter :: mne=400000
      integer, parameter :: mns=600000
      integer, parameter :: mnv=50
      integer, parameter :: mnope=50 !max # of open bnd segments
      integer, parameter :: mnond=1000 !max # of open bnd nodes in each segment
      integer, parameter :: mnei=20 !neighbor
      integer, parameter :: nbyte=4
      real, parameter :: small1=1.e-2 !used to check area ratios

!     netcdf related variables
      integer :: sid, ncids(100),one_dim ! Netcdf file IDs
      integer :: latvid, lonvid, zmvid, hvid ! positional variables
      integer :: xdid, xvid ! longitude index
      integer :: ydid, yvid ! latitude index
      integer :: ldid, lvid ! vertical level, 1 is top
      integer :: svid, tvid ! salt & temp variable IDs
      integer :: uvid, vvid ! vel variable IDs
      integer :: evid ! SSH variable IDs
      integer :: var1d_dims(1),var4d_dims(4)
      integer, dimension(nf90_max_var_dims) :: dids

!     Local variables for data
      real (kind = 4), allocatable :: xind(:), yind(:), lind(:)
!     Lat, lon, bathymetry
      real (kind = 4), allocatable :: lat(:,:), lon(:,:), hnc(:)
!     Vertical postion, salinity, and temperature
      real (kind = 4), allocatable :: zm(:,:,:),salt(:,:,:),temp(:,:,:), &
     &uvel(:,:,:),vvel(:,:,:),ssh(:,:)
      integer, allocatable :: kbp(:,:),ihope(:,:)
!     File names for netcdf files
      character(len=1024) :: ncfile1
!     Command line arguments
      character(len=1024) :: s, yr, md
      character(len=4) :: iyear_char
      character(len=1) :: char1,char2
      character(len=2) :: char3,char4
      character(len=2) :: hr_char(4)
!     external function for number of command line arguments
      integer :: iargc

      integer :: status ! netcdf local status variable
      integer :: ier ! allocate error return.
      integer :: ixlen, iylen, ilen ! sampled lengths in each coordinate direction
      integer :: ixlen1, iylen1,ixlen2, iylen2 !reduced indices for CORIE grid to speed up interpolation
!     integer :: i,j,k,i1,i2,i3,j1,j2,j3
      integer :: elnode(4,mne)
      integer :: elside(4,mne)
      integer :: isdel(2,mns)
      integer :: indel(mnei,mnp)
      dimension xl(mnp),yl(mnp),dp(mnp),i34(mne),include2(mnp),imap(mnp)
      dimension ztot(0:mnv),sigma(mnv),cs(mnv),iest(mnp),ixy(mnp,3),arco(3,mnp)
      dimension wild(100),wild2(100,2)
      dimension tempout(mnv,mnp), saltout(mnv,mnp),month_day(12)
      dimension uout(mnp,mnv),vout(mnp,mnv),eout(mnp),su2(mns,mnv),sv2(mns,mnv)
      dimension tsd(mns,mnv),ssd(mns,mnv),tsel(mnv,mne,2)
      dimension nne(mnp),ic3(4,mne),nx(4,4,3),isidenode(2,mns)
      dimension xcj(mns),ycj(mns),nond(mnope),iond(mnope,mnond),iob(mnope),iond2(mnope*mnond)
      dimension ndays_mon(12)
      allocatable :: z(:,:),sigma_lcl(:,:),kbp2(:)
      real*8 :: aa1(1)

!     First statement
      ndays_mon=(/31,28,31,30,31,30,31,31,30,31,30,31/)  !# of days in each month for non-leap yr
      hr_char=(/'03','09','15','21'/) !each day has 4 starting hours in ROMS

      open(10,file='gen_nu.in',status='old')
      read(10,*) tem_outside,sal_outside !T,S values for pts outside bg grid in nc
      read(10,*) dtout,nt_out !time step in .nc [sec], output stride
!      read(10,*) nob,iob(1:nob) !# of open bnds that need *3D.th; list of IDs
      read(10,*) istart_year,istart_mon,istart_day !in PST
      read(10,*) nndays !# of days to process
      close(10)
!      if(tem_es<0.or.sal_es<0.or.tem_outside<0.or.sal_outside<0) &
!     &stop 'Invalid T,S constants'

!     Read in hgrid and vgrid
!      open(17,file='estuary.gr3',status='old')
      open(16,file='hgrid.ll',status='old')
      open(14,file='hgrid.gr3',status='old') !only need depth info and connectivity
      open(15,file='include.gr3',status='old')
      open(19,file='vgrid.in',status='old')
      open(11,file='fort.11',status='replace')
      read(14,*); read(14,*)ne,np
      if(np.gt.mnp.or.ne.gt.mne) then
        write(*,*)'Increase mnp/mne'
        stop
      endif
      read(15,*); read(15,*)
      read(16,*)
      read(16,*)
!      read(17,*)
!      read(17,*)
      do i=1,np
        read(14,*)j,xtmp,ytmp,dp(i)
        read(16,*)j,xl(i),yl(i) !,dp(i)
        read(15,*)j,xtmp,ytmp,tmp
        include2(i)=ceiling(tmp)
      enddo !i
      do i=1,ne
        read(14,*)j,i34(i),(elnode(l,i),l=1,i34(i))
      enddo !i
!     Open bnds
      read(14,*) nope
      read(14,*) neta
      ntot=0
      if(nope>mnope) stop 'Increase mnope (2)'
      do k=1,nope
        read(14,*) nond(k)
        if(nond(k)>mnond) stop 'Increase mnond'
        do i=1,nond(k)
          read(14,*) iond(k,i)
        enddo
      enddo

!      nond0=0
!      do i=1,nob
!        ibnd=iob(i)
!        do j=1,nond(ibnd)
!          nond0=nond0+1
!          iond2(nond0)=iond(ibnd,j)
!        enddo !j
!      enddo !i

      close(14)
      close(16)

!     V-grid
      read(19,*)ivcor
      read(19,*)nvrt
      rewind(19)
      allocate(sigma_lcl(nvrt,np),z(np,nvrt),kbp2(np))
      call get_vgrid_single('vgrid.in',np,nvrt,ivcor,kz,h_s,h_c,theta_b,theta_f,ztot,sigma,sigma_lcl,kbp2)

!     Compute z-coord.
      do i=1,np
        if(ivcor==2) then
          call zcor_SZ_single(max(0.11,dp(i)),0.,0.1,h_s,h_c,theta_b,theta_f,kz,nvrt,ztot,sigma,z(i,:),idry,kbp2(i))
        else if(ivcor==1) then
          z(i,kbp2(i):nvrt)=max(0.11,dp(i))*sigma_lcl(kbp2(i):nvrt,i)
        else
          write(11,*)'Unknown ivcor:',ivcor
          stop
        endif

        !Extend below bottom for interpolation later
        z(i,1:kbp2(i)-1)=z(i,kbp2(i))
      enddo !i

!     Compute geometry
      do k=3,4
        do i=1,k
          do j=1,k-1
            nx(k,i,j)=i+j
            if(nx(k,i,j)>k) nx(k,i,j)=nx(k,i,j)-k
            if(nx(k,i,j)<1.or.nx(k,i,j)>k) then
              write(*,*)'nx wrong',i,j,k,nx(k,i,j)
              stop
            endif
          enddo !j
        enddo !i
      enddo !k

      do i=1,np
        nne(i)=0
      enddo

      do i=1,ne
        do j=1,i34(i)
          nd=elnode(j,i)
          nne(nd)=nne(nd)+1
          if(nne(nd)>mnei) then
            write(11,*)'Too many neighbors',nd
            stop
          endif
          indel(nne(nd),nd)=i
        enddo
      enddo

!     Compute ball info; this won't be affected by re-arrangement below
      do i=1,ne
        do j=1,i34(i)
          ic3(j,i)=0 !index for bnd sides
          nd1=elnode(nx(i34(i),j,1),i)
          nd2=elnode(nx(i34(i),j,2),i)
          do k=1,nne(nd1)
            ie=indel(k,nd1)
            if(ie/=i.and.(elnode(1,ie)==nd2.or.elnode(2,ie)==nd2.or.elnode(3,ie)==nd2.or.(i34(ie)==4.and.elnode(4,ie)==nd2))) ic3(j,i)=ie
          enddo !k
        enddo !j
      enddo !i

      ns=0 !# of sides
      do i=1,ne
        do j=1,i34(i)
          nd1=elnode(nx(i34(i),j,1),i)
          nd2=elnode(nx(i34(i),j,2),i)
          if(ic3(j,i)==0.or.i<ic3(j,i)) then !new sides
            ns=ns+1
            if(ns>mns) then
              write(11,*)'Too many sides'
              stop
            endif
            elside(j,i)=ns
            isdel(1,ns)=i
            isidenode(1,ns)=nd1
            isidenode(2,ns)=nd2
            xcj(ns)=(xl(nd1)+xl(nd2))/2
            ycj(ns)=(yl(nd1)+yl(nd2))/2
!            dps(ns)=(dp(nd1)+dp(nd2))/2
!            distj(ns)=dsqrt((x(nd2)-x(nd1))**2+(y(nd2)-y(nd1))**2)
!            if(distj(ns)==0) then
!              write(11,*)'Zero side',ns
!              stop
!            endif
!            thetan=datan2(x(nd1)-x(nd2),y(nd2)-y(nd1))
!            snx(ns)=dcos(thetan)
!            sny(ns)=dsin(thetan)

            isdel(2,ns)=ic3(j,i) !bnd element => bnd side
!           Corresponding side in element ic3(j,i)
            if(ic3(j,i)/=0) then !old internal side
              iel=ic3(j,i)
              index=0
              do k=1,i34(iel)
                if(ic3(k,iel)==i) then
                  index=k
                  exit
                endif
              enddo !k
              if(index==0) then
                write(11,*)'Wrong ball info',i,j
                stop
              endif
              elside(index,iel)=ns
            endif !ic3(j,i).ne.0
          endif !ic3(j,i)==0.or.i<ic3(j,i)
        enddo !j=1,i34
      enddo !i=1,ne

      if(ns<ne.or.ns<np) then
        write(11,*)'Weird grid with ns < ne or ns < np',np,ne,ns
        stop
      endif
     print*, 'Done computing geometry...'

!     open(54,file='TEM_nu.in',form='unformatted',status='replace')

!start11
!     Define limits for variables for sanity checks
      tempmin=0
      tempmax=30
      saltmin=0
      saltmax=37
      vmag_max=10 !max. |u| or |v|
      ssh_max=8 !max. of |SSH|

!     Assuming elev, u,v, S,T have same dimensions (and time step) and grids do not change over time
      timeout=-dtout !sec
      irecout=0
      irecout2=0 !final time record # in the output
      iyr_now=istart_year
      imon_now=istart_mon
      iday_now=istart_day-1

      do ifile=1,nndays !days; each day has 4 files
        !Construct the nc file name
        if(mod(iyr_now,4)==0) then
          ndays_mon(2)=29
        else
          ndays_mon(2)=28
        endif

        iday_now=iday_now+1
        if(iday_now>ndays_mon(imon_now)) then
          iday_now=iday_now-ndays_mon(imon_now)
          imon_now=imon_now+1
          if(imon_now>12) then
            imon_now=imon_now-12
            iyr_now=iyr_now+1
          endif
        endif

        write(iyear_char,'(i4)')iyr_now
        write(char3,'(i2.2)')imon_now
        write(char4,'(i2.2)')iday_now
        !print*, 'imon_now=',char3,char4,iyear_char

        do isub=1,4 !# of files each day
!--------------------------------------------------------------------
!     Open nc file
      !ncfile1='~yinglong/vims20/CA_DWR/ROMS/ca_subSFB_fcst_'//iyear_char//char3//char4//hr_char(isub)//'.nc'
      ncfile1='/sanstore/schism_archive/roms/ca_subSFB_fcst_'//iyear_char//char3//char4//hr_char(isub)//'.nc'
      print*, 'doing ',trim(ncfile1)
      status = nf90_open(trim(ncfile1),nf90_nowrite, sid)
      !Define junk value for sid; the test is abs()>rjunk
      rjunk=9998
      call check(status)

!     Get dimnesions from first file, assumed same for all variables
      if(ifile==1.and.isub==1) then
        status = nf90_inq_varid(sid, "temp", svid)
        print*, 'Done reading variable ID'

!       WARNING: indices reversed from ncdump!
        status = nf90_Inquire_Variable(sid, svid, dimids = dids)
        status = nf90_Inquire_Dimension(sid, dids(1), len = ixlen)
        status = nf90_Inquire_Dimension(sid, dids(2), len = iylen)
        status = nf90_Inquire_Dimension(sid, dids(3), len = ilen)
        status = nf90_Inquire_Dimension(sid, dids(4), len = ntime)
        print*, 'ixlen,iylen,ilen,ntime= ',ixlen,iylen,ilen,ntime

!       allocate arrays
        allocate(xind(ixlen),stat=ier)
        allocate(yind(iylen),stat=ier)
        allocate(lind(ilen),stat=ier)
        allocate(lat(ixlen,iylen))
        allocate(lon(ixlen,iylen))
        allocate(zm(ixlen, iylen, ilen))
!        allocate(hnc(ixlen,iylen))
        allocate(hnc(ilen))
        allocate(kbp(ixlen,iylen))
        allocate(ihope(ixlen,iylen))
        allocate(uvel(ixlen,iylen,ilen),stat=ier)
        allocate(vvel(ixlen,iylen,ilen),stat=ier)
        allocate(salt(ixlen,iylen,ilen),stat=ier)
        allocate(temp(ixlen,iylen,ilen),stat=ier)
        allocate(ssh(ixlen,iylen),stat=ier)
        uvel=0; vvel=0; ssh=0

!       get static info (lat/lon grids etc)
        status = nf90_inq_varid(sid, "lon", xvid)
        status = nf90_get_var(sid, xvid, xind)
        status = nf90_inq_varid(sid, "lat", yvid)
        status = nf90_get_var(sid, yvid, yind)
        status = nf90_inq_varid(sid, "depth", hvid)
        status = nf90_get_var(sid, hvid, hnc)

!       processing static info
        do i=1,ixlen
          lon(i,:)=xind(i)
        enddo !i
        do j=1,iylen
          lat(:,j)=yind(j)
        enddo !j
        lon=lon-360e0 !convert to our long.

!       Compute z-coord. (assuming eta=0)
!       WARNING: In zm(), 1 is bottom; ilen is surface (SELFE convention)
        do i=1,ixlen
          do j=1,iylen
            do k=1,ilen
              zm(i,j,k)=-hnc(1+ilen-k)
            enddo !k
          enddo !j
        enddo !i

!       Get rest of varid's
        status = nf90_inq_varid(sid, "temp", tvid)
!        status = nf90_inq_varid(sid, "zeta", evid)
!        status = nf90_inq_varid(sid, "u", uvid)
!        status = nf90_inq_varid(sid, "v", vvid)

!       Arrays no longer used after this: hnc,lind
        deallocate(hnc,lind)
      endif !ifile==1 && isub==1

!     Vertical convention follows SELFE from now on; i.e., 1 is at bottom
!     Compute bottom indices
      if(ifile==1.and.isub==1) then
        ilo=5 !5th record is 00:00 PST
      else
        ilo=2 !there is an overlap between files
      endif !ifile

      do it2=ilo,ntime
        irecout=irecout+1
        timeout=timeout+dtout !sec
        print*, 'Time out (days)=',timeout/86400

        !Make sure it2=ilo is output (some init. work below)
        if(mod(irecout-1,nt_out)/=0) cycle
        irecout2=irecout2+1

!       Read T,S,u,v,SSH
!       WARNING! Make sure the order of vertical indices is 1
!                at bottom, ilen at surface; revert if necessary!
        status = nf90_get_var(sid,svid,salt(:,:,ilen:1:-1),start=(/1,1,1,it2/),count=(/ixlen,iylen,ilen,1/))
        status = nf90_get_var(sid,tvid,temp(:,:,ilen:1:-1),start=(/1,1,1,it2/),count=(/ixlen,iylen,ilen,1/))

        kbp=1  !init.
        do i=1,ixlen
          do j=1,iylen
            if(abs(salt(i,j,ilen))>rjunk) kbp(i,j)=-1 !dry

            if(kbp(i,j)==1) then !wet
              !Extend near bottom
              klev0=-1 !flag
              do k=1,ilen
                if(abs(salt(i,j,k))<=rjunk) then
                  klev0=k; exit
                endif
              enddo !k
              if(klev0<=0) then
                write(11,*)'Impossible (1):',i,j,salt(i,j,ilen)
                stop
              endif !klev0
              salt(i,j,1:klev0-1)=salt(i,j,klev0)
              temp(i,j,1:klev0-1)=temp(i,j,klev0)
!              uvel(i,j,1:klev0-1)=uvel(i,j,klev0)
!              vvel(i,j,1:klev0-1)=vvel(i,j,klev0)

              !Check
              do k=1,ilen
                if(salt(i,j,k)<saltmin.or.salt(i,j,k)>saltmax.or. &
                   temp(i,j,k)<tempmin.or.temp(i,j,k)>tempmax) then
!                   abs(uvel(i,j,k))>vmag_max.or.abs(vvel(i,j,k))>vmag_max.or. &
!                   abs(ssh(i,j))>ssh_max) then
                  write(11,*)'Fatal: no valid values:',it2,i,j,k,salt(i,j,k), &
     &temp(i,j,k) !,uvel(i,j,k),vvel(i,j,k),ssh(i,j)
                  stop
                endif
              enddo !k
            endif !kbp=1
          enddo !j
        enddo !i

!       Compute S,T etc@ invalid pts based on nearest neighbor
!       Search around neighborhood of a pt
        do i=1,ixlen
          do j=1,iylen
            if(kbp(i,j)==-1) then !invalid pts
              !Compute max possible tier # for neighborhood
              mmax=max(i-1,ixlen-i,j-1,iylen-j)

              m=0 !tier #
              loop6: do
                m=m+1
                do ii=max(-m,1-i),min(m,ixlen-i)
                  i3=max(1,min(ixlen,i+ii))
                  do jj=max(-m,1-j),min(m,iylen-j)
                    j3=max(1,min(iylen,j+jj))
                    if(kbp(i3,j3)==1) then !found
                      i1=i3; j1=j3
                      exit loop6
                    endif
                  enddo !jj
                enddo !ii

                if(m==mmax) then
                  write(11,*)'Max. exhausted:',i,j,mmax
                  write(11,*)'kbp'
                  do ii=1,ixlen
                    do jj=1,iylen
                      write(11,*)ii,jj,kbp(ii,jj)
                    enddo !jj
                  enddo !ii
                  stop
                endif
              end do loop6

              salt(i,j,1:ilen)=salt(i1,j1,1:ilen)
              temp(i,j,1:ilen)=temp(i1,j1,1:ilen)
!              uvel(i,j,1:ilen)=uvel(i1,j1,1:ilen)
!              vvel(i,j,1:ilen)=vvel(i1,j1,1:ilen)
!              ssh(i,j)=ssh(i1,j1)

              !Check
              do k=1,ilen
                if(salt(i,j,k)<saltmin.or.salt(i,j,k)>saltmax.or. &
                   temp(i,j,k)<tempmin.or.temp(i,j,k)>tempmax) then
!                   abs(uvel(i,j,k))>vmag_max.or.abs(vvel(i,j,k))>vmag_max.or. &
!                   abs(ssh(i,j))>ssh_max) then
                  write(11,*)'Fatal: no valid values after searching:',it2,i,j,k,salt(i,j,k), &
     &temp(i,j,k) !,uvel(i,j,k),vvel(i,j,k),ssh(i,j)
                  stop
                endif
              enddo !k

            endif !kbp(i,j)==-1
          enddo !j=iylen1,iylen2
        enddo !i=ixlen1,ixlen2

!       Do something for 1st step only
        if(ifile==1.and.isub==1.and.it2==ilo) then
!         Test outputs
          icount=0
          do i=1,ixlen
            do j=1,iylen
              icount=icount+1
!              write(95,*)icount,lon(i,j),lat(i,j),ssh(i,j)
!              write(96,*)icount,lon(i,j),lat(i,j),uvel(i,j,ilen)
!              write(97,*)icount,lon(i,j),lat(i,j),vvel(i,j,ilen)
              write(98,*)icount,lon(i,j),lat(i,j),salt(i,j,1) !,i,j
              write(99,*)icount,lon(i,j),lat(i,j),temp(i,j,ilen)
              write(100,*)icount,lon(i,j),lat(i,j),temp(i,j,1)
            enddo !j
          enddo !i
          print*, 'done outputting test outputs for nc'

!         Find parent elements for hgrid.ll
          ixy=0
          loop4: do i=1,np
            if(include2(i)==0) cycle loop4

            do ix=1,ixlen-1
              do iy=1,iylen-1
                x1=lon(ix,iy); x2=lon(ix+1,iy); x3=lon(ix+1,iy+1); x4=lon(ix,iy+1)
                y1=lat(ix,iy); y2=lat(ix+1,iy); y3=lat(ix+1,iy+1); y4=lat(ix,iy+1)
                a1=abs(signa(xl(i),x1,x2,yl(i),y1,y2))
                a2=abs(signa(xl(i),x2,x3,yl(i),y2,y3))
                a3=abs(signa(xl(i),x3,x4,yl(i),y3,y4))
                a4=abs(signa(xl(i),x4,x1,yl(i),y4,y1))
                b1=abs(signa(x1,x2,x3,y1,y2,y3))
                b2=abs(signa(x1,x3,x4,y1,y3,y4))
                rat=abs(a1+a2+a3+a4-b1-b2)/(b1+b2)
                if(rat<small1) then
                  ixy(i,1)=ix; ixy(i,2)=iy
!                 Find a triangle
                  in=0 !flag
                  do l=1,2 !split quad
                    ap=abs(signa(xl(i),x1,x3,yl(i),y1,y3))
                    if(l==1) then !nodes 1,2,3
                      bb=abs(signa(x1,x2,x3,y1,y2,y3))
                      wild(l)=abs(a1+a2+ap-bb)/bb
                      if(wild(l)<small1*5) then
                        in=1
                        arco(1,i)=max(0.,min(1.,a2/bb))
                        arco(2,i)=max(0.,min(1.,ap/bb))
                        exit
                      endif
                    else !nodes 1,3,4
                      bb=abs(signa(x1,x3,x4,y1,y3,y4))
                      wild(l)=abs(a3+a4+ap-bb)/bb
                      if(wild(l)<small1*5) then
                        in=2
                        arco(1,i)=max(0.,min(1.,a3/bb))
                        arco(2,i)=max(0.,min(1.,a4/bb))
                        exit
                      endif
                    endif
                  enddo !l=1,2
                  if(in==0) then
                    write(11,*)'Cannot find a triangle:',(wild(l),l=1,2)
                    stop
                  endif
                  ixy(i,3)=in
                  arco(3,i)=max(0.,min(1.,1-arco(1,i)-arco(2,i)))
                  cycle loop4
                endif !rat<small1
              enddo !iy=iylen1,iylen2-1
            enddo !ix=ixlen1,ixlen2-1
          end do loop4 !i=1,np
        endif !ifile==1.and.it2==

!       Do interpolation: all at 1st  step, bnd only for the rest
!        if(ifile==1.and.it2==ilo) then
!          iup=np
!        else
!          iup=nond0
!        endif

!        print*, 'Time out (days)=',timeout/86400,it2,ilo,ntime
        tempout=-99; saltout=-99
        npout=0 !# of _valid_ output points
        do ii=1,np
          i=ii !iond2(ii) !node #

          if(ixy(i,1)==0.or.ixy(i,2)==0) then
            ! write(*,*)'Cannot find a parent element:',i
            ! if(include2(i)/=0.) then
              !npout=npout+1
              !imap(npout)=i
              tempout(:,i)=tem_outside
              saltout(:,i)=sal_outside
            ! endif
            if (include2(i)/=0.) then
                print*,'This node is not covered by the data', i
            endif
          else !found parent
            npout=npout+1
            imap(npout)=i
            ix=ixy(i,1); iy=ixy(i,2); in=ixy(i,3)
            !Find vertical level
            tempout(:,i)=-999
            do k=1,nvrt
              if(kbp(ix,iy)==-1) then
                lev=ilen-1; vrat=1
              else if(z(i,k)<=zm(ix,iy,kbp(ix,iy))) then
                lev=kbp(ix,iy); vrat=0
              else if(z(i,k)>=zm(ix,iy,ilen)) then !above f.s.
                lev=ilen-1; vrat=1
              else
                lev=-99 !flag
                do kk=1,ilen-1
                  if(z(i,k)>=zm(ix,iy,kk).and.z(i,k)<=zm(ix,iy,kk+1)) then
                    lev=kk
                    vrat=(zm(ix,iy,kk)-z(i,k))/(zm(ix,iy,kk)-zm(ix,iy,kk+1))
                    exit
                  endif
                enddo !kk
                if(lev==-99) then
                  write(11,*)'Cannot find a level:',i,k,z(i,k),(zm(ix,iy,l),l=1,kbp(ix,iy))
                  stop
                endif
              endif

!              write(18,*)i,k,ix,iy,lev,vrat,kbp(ix,iy)
              wild2(1,1)=temp(ix,iy,lev)*(1-vrat)+temp(ix,iy,lev+1)*vrat
              wild2(1,2)=salt(ix,iy,lev)*(1-vrat)+salt(ix,iy,lev+1)*vrat
              wild2(2,1)=temp(ix+1,iy,lev)*(1-vrat)+temp(ix+1,iy,lev+1)*vrat
              wild2(2,2)=salt(ix+1,iy,lev)*(1-vrat)+salt(ix+1,iy,lev+1)*vrat
              wild2(3,1)=temp(ix+1,iy+1,lev)*(1-vrat)+temp(ix+1,iy+1,lev+1)*vrat
              wild2(3,2)=salt(ix+1,iy+1,lev)*(1-vrat)+salt(ix+1,iy+1,lev+1)*vrat
              wild2(4,1)=temp(ix,iy+1,lev)*(1-vrat)+temp(ix,iy+1,lev+1)*vrat
              wild2(4,2)=salt(ix,iy+1,lev)*(1-vrat)+salt(ix,iy+1,lev+1)*vrat

!              wild2(5,1)=uvel(ix,iy,lev)*(1-vrat)+uvel(ix,iy,lev+1)*vrat
!              wild2(5,2)=vvel(ix,iy,lev)*(1-vrat)+vvel(ix,iy,lev+1)*vrat
!              wild2(6,1)=uvel(ix+1,iy,lev)*(1-vrat)+uvel(ix+1,iy,lev+1)*vrat
!              wild2(6,2)=vvel(ix+1,iy,lev)*(1-vrat)+vvel(ix+1,iy,lev+1)*vrat
!              wild2(7,1)=uvel(ix+1,iy+1,lev)*(1-vrat)+uvel(ix+1,iy+1,lev+1)*vrat
!              wild2(7,2)=vvel(ix+1,iy+1,lev)*(1-vrat)+vvel(ix+1,iy+1,lev+1)*vrat
!              wild2(8,1)=uvel(ix,iy+1,lev)*(1-vrat)+uvel(ix,iy+1,lev+1)*vrat
!              wild2(8,2)=vvel(ix,iy+1,lev)*(1-vrat)+vvel(ix,iy+1,lev+1)*vrat
              if(in==1) then
                tempout(k,i)=wild2(1,1)*arco(1,i)+wild2(2,1)*arco(2,i)+wild2(3,1)*arco(3,i)
                saltout(k,i)=wild2(1,2)*arco(1,i)+wild2(2,2)*arco(2,i)+wild2(3,2)*arco(3,i)
              else
                tempout(k,i)=wild2(1,1)*arco(1,i)+wild2(3,1)*arco(2,i)+wild2(4,1)*arco(3,i)
                saltout(k,i)=wild2(1,2)*arco(1,i)+wild2(3,2)*arco(2,i)+wild2(4,2)*arco(3,i)
              endif

              !Check
              if(tempout(k,i)<tempmin.or.tempout(k,i)>tempmax.or. &
                saltout(k,i)<saltmin.or.saltout(k,i)>saltmax) then
                write(11,*)'Interpolated values invalid:',i,k,tempout(k,i),saltout(k,i)
                stop
              endif

              !Correct near surface T bias
              if(kbp(ix,iy)/=-1.and.z(i,k)>-10) tempout(k,i)=tempout(k,i)-1

!             Enforce lower bound for temp. for eqstate
              tempout(k,i)=max(0.,tempout(k,i))
            enddo !k=1,nvrt
          endif !ixy(i,1)==0.or.
        enddo !ii

!       Output
!        print*, 'Time out (days)=',timeout/86400,it2,ilo,ntime
        print*, 'outputting at day ',timeout/86400,npout
!        write(54)timeout
!        do i=1,np
!          write(54)tempout(1:nvrt,i)
!        enddo !i
        if(ifile==1.and.isub==1.and.it2==ilo) then
          iret=nf90_create('TEM_nu.nc',OR(NF90_CLOBBER,NF90_NETCDF4),ncid_T)
          iret=nf90_def_dim(ncid_T,'node',npout,node_dim)
          iret=nf90_def_dim(ncid_T,'nLevels',nvrt,nv_dim)
          iret=nf90_def_dim(ncid_T,'one',1,one_dim)
          iret=nf90_def_dim(ncid_T,'time', NF90_UNLIMITED,itime_dim)

          var1d_dims(1)=itime_dim
          iret=nf90_def_var(ncid_T,'time',NF90_DOUBLE,var1d_dims,itime_id_T)
          var1d_dims(1)=node_dim
          iret=nf90_def_var(ncid_T,'map_to_global_node',NF90_INT,var1d_dims,id_map_T)
          var4d_dims(1)=one_dim; var4d_dims(2)=nv_dim;
          var4d_dims(3)=node_dim; var4d_dims(4)=itime_dim
          iret=nf90_def_var(ncid_T,'tracer_concentration',NF90_FLOAT,var4d_dims,ivar_T)
          iret=nf90_enddef(ncid_T)

          iret=nf90_create('SAL_nu.nc',OR(NF90_CLOBBER,NF90_NETCDF4),ncid_S)
          iret=nf90_def_dim(ncid_S,'node',npout,node_dim)
          iret=nf90_def_dim(ncid_S,'nLevels',nvrt,nv_dim)
          iret=nf90_def_dim(ncid_S,'one',1,one_dim)
          iret=nf90_def_dim(ncid_S,'time', NF90_UNLIMITED,itime_dim)

          var1d_dims(1)=itime_dim
          iret=nf90_def_var(ncid_S,'time',NF90_DOUBLE,var1d_dims,itime_id_S)
          var1d_dims(1)=node_dim
          iret=nf90_def_var(ncid_S,'map_to_global_node',NF90_INT,var1d_dims,id_map_S)
          iret=nf90_def_var(ncid_S,'tracer_concentration',NF90_FLOAT,var4d_dims,ivar_S)
          iret=nf90_enddef(ncid_S)
        endif !ifile

        aa1(1)=timeout/86400
        iret=nf90_put_var(ncid_T,itime_id_T,aa1,(/irecout2/),(/1/))
        iret=nf90_put_var(ncid_S,itime_id_S,aa1,(/irecout2/),(/1/))
        iret=nf90_put_var(ncid_T,id_map_T,imap(1:npout),(/1/),(/npout/))
        iret=nf90_put_var(ncid_S,id_map_S,imap(1:npout),(/1/),(/npout/))
        iret=nf90_put_var(ncid_T,ivar_T,tempout(:,imap(1:npout)),(/1,1,1,irecout2/),(/1,nvrt,npout,1/))
        iret=nf90_put_var(ncid_S,ivar_S,saltout(:,imap(1:npout)),(/1,1,1,irecout2/),(/1,nvrt,npout,1/))
        ! iret=nf90_put_var(ncid_T,ivar_T,tempout(:,:npout),(/1,1,1,irecout2/),(/1,nvrt,npout,1/))
        ! iret=nf90_put_var(ncid_S,ivar_S,saltout(:,:npout),(/1,1,1,irecout2/),(/1,nvrt,npout,1/))
      enddo !it2=ilo,ntime

      status = nf90_close(sid)
!end11
      print*, 'done reading nc for file ',ifile
!--------------------------------------------------------------------
        enddo !isub=1,4
      enddo !ifile=1,

      iret=nf90_close(ncid_T)
      iret=nf90_close(ncid_S)

!      deallocate(lat,lon,zm,h,kbp,ihope,xind,yind,lind,salt,temp)

      print*, 'Finished'
!     End of main

!     Subroutines
      contains
!     Internal subroutine - checks error status after each netcdf,
!     prints out text message each time an error code is returned.
      subroutine check(status)
      integer, intent ( in) :: status

      if(status /= nf90_noerr) then
        print *, trim(nf90_strerror(status))
        print *, 'failed to open nc files'
        stop
      end if
      end subroutine check

      end program gen_hot

      function signa(x1,x2,x3,y1,y2,y3)

      signa=((x1-x3)*(y2-y3)-(x2-x3)*(y1-y3))/2

      return
      end

max_stars_repo_name: CADWRDeltaModeling/schimpy

max_issues_repo_name: CADWRDeltaModeling/schimpy

max_forks_repo_name: CADWRDeltaModeling/schimpy


text:
subroutine SHSCouplingMatrix(kij, galpha, lmax, nmax, exitstatus)
!------------------------------------------------------------------------------
!
!   SHSCouplingMatrix returns the spherical harmonic coupling matrix that
!   relates the expectation of the function expressed in a subset of the best-
!   localized Slepian functions to the expectation of the global power spectrum
!   (assumed to be stationary). The Slepian functions are determined by a call
!   to either (1) `SHReturnTapers` and then `SHRotateTapers`, or
!   (2) `SHReturnTapersMap`. Each row of galpha contains the (lmax+1)**2
!   spherical harmonic coefficients of a Slepian function that can be unpacked
!   using `SHVectorToCilm`. The Slepian functions must be normalized to have
!   unit power (that is the sum of the coefficients squared is 1).
!
!   The relationship between the global and localized power spectra is given by
!
!       < S_{\tilde{f}}(l) > = \sum_{l'=0}^lmax K_{ll'} S_{f}(l')
!
!   where S_{\tilde{f}} is the expectation of the power spectrum at degree l of
!   the function expressed in Slepian functions, S_{f}(l') is the expectation
!   of the global power spectrum, and < ... > is the expectation operator.
!   The coupling matrix is given explicitly by
!
!   K_{ij} = \frac{1}{2l'+1} Sum_{m=-l}^l Sum_{m'=-l'}^l'
!            ( Sum_{alpha=1}^nmax g_{l'm'}(alpha) g_{lm}(alpha) )**2
!
!   Calling Parameters
!
!       IN
!           galpha  Matrix of dimension ( (LMAX+1)**2, nmax) containing the
!                   spherical harmonic coefficients of the Slepian functions.
!                   Each column corresponds to a Slepian function ordered from
!                   best to worst concentrated.
!           lmax    Maximum spherical harmonic degree of the Slepian functions.
!           nmax    Maximum number of Slepian coefficients used to construct
!                   the function.
!
!       OUT
!           kij     Spherical harmonic coupling matrix of dimension
!                   ((LMAX+1)**2, nmax).
!
!       OPTIONAL (OUT)
!           exitstatus  If present, instead of executing a STOP when an error
!                       is encountered, the variable exitstatus will be
!                       returned describing the error.
!                       0 = No errors;
!                       1 = Improper dimensions of input array;
!                       2 = Improper bounds for input variable;
!                       3 = Error allocating memory;
!                       4 = File IO error.
!
!   Copyright (c) 2005-2019, SHTOOLS
!   All rights reserved.
!
!------------------------------------------------------------------------------
    use ftypes

    implicit none

    real(dp), intent(out) :: kij(:,:)
    real(dp), intent(in) :: galpha(:,:)
    integer, intent(in) :: lmax, nmax
    integer, intent(out), optional :: exitstatus
    integer :: l, lp, m, mp, alpha, ind, indp
    real(dp) :: temp

    if (present(exitstatus)) exitstatus = 0

    if (size(kij(:,1)) < (lmax+1) .or. size(kij(1,:)) < (lmax+1)) then
        print*, "Error --- SHSCouplingMatrix"
        print*, "KIJ must be dimensioned as (LMAX+1, LMAX+1)."
        print*, "LMAX = ", lmax
        print*, "Dimension of KIJ = ", size(kij(:,1)), size(kij(1,:))
        if (present(exitstatus)) then
            exitstatus = 1
            return
        else
            stop
        end if

    else if (size(galpha(:,1)) < (lmax+1)**2 .or. &
            size(galpha(1,:)) < nmax) then
        print*, "Error --- SHSCouplingMatrix"
        print*, "GALPHA must be dimensioned as ( (LMAX+1)**2, NMAX )."
        print*, "LMAX = ", lmax
        print*, "NMAX = ", nmax
        print*, "Input array is dimensioned as ", size(galpha(:,1)), &
                size(galpha(1,:))
        if (present(exitstatus)) then
            exitstatus = 1
            return
        else
            stop
        end if

    end if

    kij = 0.0_dp

    ! Calculate k(l, lp) * (2lp +1 ), which is symmetric
    do l = 0, lmax
        do lp = l, lmax, 1
            do m = -l, l, 1

                if (m >= 0) then
                    ind = l**2+m+1
                else
                    ind = l**2+l+abs(m)+1
                end if

                do mp = -lp, lp, 1

                    if (mp >= 0) then
                        indp = lp**2+mp+1
                    else
                        indp = lp**2+lp+abs(mp)+1
                    end if

                    temp = 0.0_dp
                    do alpha = 1, nmax, 1

                       temp = temp + galpha(indp, alpha) * galpha(ind, alpha)

                    end do

                    kij(l+1, lp+1) = kij(l+1, lp+1) + temp**2

                end do

            end do

        end do

    end do

    ! Fill lower half of matrix
    do l = 1, lmax, 1

        do lp = 0, l-1, 1

            kij(l+1, lp+1) = kij(lp+1, l+1)

        end do

    end do

    ! Multiply by 1/(2lp+1)
    do lp = 0, lmax

        kij(1:lmax+1, lp+1) = kij(1:lmax+1, lp+1) / dble(2*lp+1)

    end do

end subroutine SHSCouplingMatrix

max_stars_repo_name: dilkins/SHTOOLS

max_issues_repo_name: dilkins/SHTOOLS

max_forks_repo_name: dilkins/SHTOOLS