text:
      SUBROUTINE CALPUV9 (ISTL_)
C
C CHANGE RECORD
C ** SUBROUTINE CALPUV9 CALCULATES THE EXTERNAL SOLUTION FOR P, UHDYE,
C ** AND VHDXE, FOR FREE SURFACE FLOWS WITH PROVISIONS FOR WETTING
C ** AND DRYING OF CELLS
C
      USE GLOBAL

      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::QSUMTMP

      IF(.NOT.ALLOCATED(QSUMTMP))THEN
        ALLOCATE(QSUMTMP(LCM))
        QSUMTMP=0.
      ENDIF
C
      IF(DEBUG)THEN
        IF(N.EQ.1.AND.ISDSOLV.EQ.1)THEN
          OPEN(1,FILE='FUV1.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
          OPEN(1,FILE='EQCOEF1.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
          OPEN(1,FILE='EQTERM1.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
          OPEN(1,FILE='FP1.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
        ENDIF
        IF(N.EQ.2.AND.ISDSOLV.EQ.1)THEN
          OPEN(1,FILE='FUV2.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
          OPEN(1,FILE='EQCOEF2.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
          OPEN(1,FILE='EQTERM2.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
          OPEN(1,FILE='FP2.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
        ENDIF
        IF(ISDSOLV.EQ.1)THEN
          OPEN(1,FILE='FUV.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
          OPEN(1,FILE='EQCOEF.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
          OPEN(1,FILE='EQTERM.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
          OPEN(1,FILE='FP.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
        ENDIF
      ENDIF
C
      DELT=DT2
      DELTD2=DT
      IF(ISTL_.EQ.2)THEN
        DELT=DT
        DELTD2=0.5*DT
      ENDIF
      DELTI=1./DELT
C
C **  CALCULATE EXTERNAL BUOYANCY INTEGRALS AT TIME LEVEL (N)
C
      IF(BSC.GT.1.E-6)THEN
        ! *** CALCULATE EXTERNAL BUOYANCY INTEGRALS AT TIME LEVEL (N)
        CALL CALEBI
C
        ! *** CALCULATE EXPLICIT EXTERNAL PRESSURE GRADIENTS
        DO L=2,LA
          !SBX(L)=0.5*SUB(L)*DYU(L)
          FPGXE(L)=ROLD*FPGXE(L)+RNEW*(
     &        -SBX(L)*HU(L)*GP*((BI2(L)+BI2(LWEST(L)))*(HP(L)-HP(LWEST(L)))
     &        +2.0*HU(L)*(BI1(L)-BI1(LWEST(L)))
     &        +(BE(L)+BE(LWEST(L)))*(BELV(L)-BELV(LWEST(L)))))
        ENDDO
        DO L=2,LA
          LS=LSC(L)
          !SBY(L)=0.5*SVB(L)*DXV(L)
          FPGYE(L)=ROLD*FPGYE(L)+RNEW*(
     &         -SBY(L)*HV(L)*GP*((BI2(L)+BI2(LS))*(HP(L)-HP(LS))
     &        +2.0*HV(L)*(BI1(L)-BI1(LS))
     &        +(BE(L)+BE(LS))*(BELV(L)-BELV(LS))))
        ENDDO
      ENDIF
C
C **  CALCULATE EXPLICIT EXTERNAL UHDYE AND VHDXE EQUATION TERMS
C **  HRU=SUB*HMU*DYU/DXU & HRV=SVB*HMV*DXV/DYV
C
      IF(ISTL_.EQ.2)THEN
        DO L=2,LA
          HUTMP(L)=0.5*(HU(L)+H1U(L))
          HVTMP(L)=0.5*(HV(L)+H1V(L))
        ENDDO
      ELSE
        DO L=2,LA
          HUTMP(L)=HU(L)
          HVTMP(L)=HV(L)
        ENDDO
      ENDIF
C
      DO L=2,LA
        LS=LSC(L)
        LN=LNC(L)
        !DXYU(L)=DXU(L)*DYU(L)
        !DXIU(L)=1./DXU(L)
        !HRUO(L)=SUBO(L)*DYU(L)*DXIU(L)
        FUHDYE(L)=UHDY1E(L)
     &      -DELTD2*SUB(L)*HRUO(L)*HUTMP(L)*(P1(L)-P1(LWEST(L)))
     &      +SUB(L)*DELT*DXIU(L)*(DXYU(L)*(TSX1(L)-RITB1*TBX1(L))
     &      +FCAXE(L)+FPGXE(L)-SNLT*FXE(L))
C
        !DXYV(L)=DXV(L)*DYV(L)
        !DYIV(L)=1./DYV(L)
        !HRVO(L)=SVBO(L)*DXV(L)*DYIV(L)
        FVHDXE(L)=VHDX1E(L)
     &      -DELTD2*SVB(L)*HRVO(L)*HVTMP(L)*(P1(L)-P1(LS ))
     &      +SVB(L)*DELT*DYIV(L)*(DXYV(L)*(TSY1(L)-RITB1*TBY1(L))
     &      -FCAYE(L)+FPGYE(L)-SNLT*FYE(L))
      ENDDO
C
      IF(ISDSOLV.GE.1.AND.DEBUG)THEN
        OPEN(1,FILE='FUV.OUT',POSITION='APPEND',STATUS='UNKNOWN')
        WRITE(1,1001)N,ISTL_
        DO L=2,LA
          WRITE(1,1001)IL(L),JL(L),UHDY1E(L),HRUO(L),HUTMP(L),P1(L),
     &        P1(LWEST(L)),TSX1(L),TBX1(L),FCAXE(L),FPGXE(L),FXE(L)
        ENDDO
        CLOSE(1)
        IF(N.EQ.1)THEN
          OPEN(1,FILE='FUV1.OUT',POSITION='APPEND',STATUS='UNKNOWN')
          WRITE(1,1001)N,ISTL_
          DO L=2,LA
            WRITE(1,1001)IL(L),JL(L),UHDY1E(L),HRUO(L),HUTMP(L),P1(L),
     &          P1(LWEST(L)),TSX1(L),TBX1(L),FCAXE(L),FPGXE(L),FXE(L)
          ENDDO
          CLOSE(1)
        ENDIF
        IF(N.EQ.2)THEN
          OPEN(1,FILE='FUV2.OUT',POSITION='APPEND',STATUS='UNKNOWN')
          WRITE(1,1001)N,ISTL_
          DO L=2,LA
            WRITE(1,1001)IL(L),JL(L),UHDY1E(L),HRUO(L),HUTMP(L),P1(L),
     &          P1(LWEST(L)),TSX1(L),TBX1(L),FCAXE(L),FPGXE(L),FXE(L)
          ENDDO
          CLOSE(1)
        ENDIF
      ENDIF
C
C **  SET IMPLICIT BOTTOM AND VEGETATION DRAG AS APPROPRIATE
C
      RCX(2:LA)=1./(1.+FBODYFXI(2:LA)) !Multiple layers no vegetation
      RCY(2:LA)=1./(1.+FBODYFYI(2:LA))

      RCX(1)=0.
      RCY(1)=0.
      RCX(LC)=0.
      RCY(LC)=0.
C
C * SINGLE LAYER NO VEGETATION
C
      IF(KC.EQ.1.AND.ISVEG.EQ.0.AND.RITB.GT.0.)THEN
        DO L=2,LA
          RCX(L)=1./( 1.+FBODYFXI(L)
     &        +RITB*DELT*HUI(L)*STBX(L)*SQRT(VU(L)*VU(L)+U(L,1)*U(L,1)))
          RCY(L)=1./( 1.+FBODYFYI(L)
     &        +RITB*DELT*HVI(L)*STBY(L)*SQRT(UV(L)*UV(L)+V(L,1)*V(L,1)))
          FUHDYE(L)=FUHDYE(L)*RCX(L)
          FVHDXE(L)=FVHDXE(L)*RCY(L)
        ENDDO
      ENDIF
C
C * SINGLE LAYER WITH VEGETATION
C
      IF(KC.EQ.1.AND.ISVEG.GE.1)THEN
        DO L=2,LA
          RCX(L)=1./( 1.+FBODYFXI(L)
     &        +RITB*DELT*HUI(L)*STBX(L)*SQRT(VU(L)*VU(L)+U(L,1)*U(L,1))
     &        +DELT*FXVEGE(L) )
          RCY(L)=1./( 1.+FBODYFYI(L)
     &        +RITB*DELT*HVI(L)*STBY(L)*SQRT(UV(L)*UV(L)+V(L,1)*V(L,1))
     &        +DELT*FYVEGE(L) )
          FUHDYE(L)=FUHDYE(L)*RCX(L)
          FVHDXE(L)=FVHDXE(L)*RCY(L)
        ENDDO
      ENDIF
C
C * MULTIPLE LAYERS WITH VEGETATION
C
      IF(KC.GT.1.AND.ISVEG.GE.1)THEN
        DO L=2,LA
          RCX(L)=1./( 1.+FBODYFXI(L)+DELT*FXVEGE(L) )
          RCY(L)=1./( 1.+FBODYFYI(L)+DELT*FYVEGE(L) )
          FUHDYE(L)=FUHDYE(L)*RCX(L)
          FVHDXE(L)=FVHDXE(L)*RCY(L)
        ENDDO
      ENDIF
C
C **  ADJUST VOLUME SOURCE AND SINKS
C
      IF(ISGWIE.EQ.0)THEN
        DO L=2,LA
          IF(QSUME(L).LE.0.)THEN
            IF(H1P(L).LE.HDRY)THEN
              QSUMTMP(L)=0.
            ELSE
              QSUMTMP(L)=-(H1P(L)-HDRY)*DXYP(L)*DELTI
              QSUMTMP(L)=MAX(QSUMTMP(L),QSUME(L))
            ENDIF
          ELSE
            QSUMTMP(L)=QSUME(L)
          ENDIF
        ENDDO
        DO L=2,LA
          DIFQVOL=QSUME(L)-QSUMTMP(L)
          DO K=1,KC
            QSUM(L,K)=QSUM(L,K)-DIFQVOL*DZC(K)
          ENDDO
          QSUME(L)=QSUMTMP(L)
        ENDDO
      ENDIF
C
C **  ADJUST SOURCES AND SINKS ESTIMATING SURFACE AND GROUNDWATER
C **  AVAILABLE FOR EVAPOTRANSPIRATON AND INFILTRATION
C
      IF(ISGWIE.GE.1)THEN
        DO L=2,LA
          RIFTR(L)=0.
          EVAPSW(L)=0.
          EVAPGW(L)=0.
          IF(H1P(L).GT.HDRY)THEN
C
C       APPLY MAXIMUM ET
C
            IF(EVAPCVT.LT.0.)THEN
              SVPW=(10.**((0.7859+0.03477*TEM(L,KC))/
     &            (1.+0.00412*TEM(L,KC))))
              EVAPT(L)=CLEVAP(L)*0.7464E-3*WINDST(L)*(SVPW-VPA(L))
     &            /PATMT(L)
            ENDIF
            EVAPSW(L)=EVAPT(L)*DXYP(L)
            RIFTR(L)=0.
C
C       CALCULATE DEPTH OF ACTIVE GROUNDWATER ELEV BELOW SURFACE
C
            DTAGW=BELV(L)-AGWELV(L)
            IF(DTAGW.GT.0.0)THEN
C
C         INFLITRATION CAN OCCUR, CALCULATE LIMITING RATE TO BRING
C         GW ELEV TO SOIL SURFACE
C
              RIFTRL=RNPOR*DTAGW*DELTI
C
C         SET RIFTRL TO MIN OF LIMITING RATE OR ACTUAL RATE
C
              RIFTRL=MIN(RIFTRM,RIFTRL)
C
C         ESTIMATE RATE BASED ON AVAILABLE SURFACE WATER
C
              RAVAIL=(H1P(L)-HDRY)*DELTI-EVAPT(L)
C
C         SET RIFTRL TO MIN OF AVAILABLE RATE OR LIMITING RATE
C
              RIFTRL=MIN(RAVAIL,RIFTRL)
C
C         CONVERT TO VOLUME FLOW UNITS
C
              RIFTR(L)=RIFTRL*DXYP(L)
            ENDIF
C
C       ADJUST VOLUME OUTFLOWS OF WET CELLS
C
            IF(QSUME(L).LT.0.0)THEN
              QSUMIET=RIFTR(L)+EVAPSW(L)
              QEAVAIL=DXYP(L)*(H1P(L)-HDRY)*DELTI-QSUMIET
              QEAVAIL=MAX(QEAVAIL,0.0)
              QEAVAIL=-QEAVAIL
              QSUMTMP(L)=MAX(QSUME(L),QEAVAIL)
            ELSE
              QSUMTMP(L)=QSUME(L)
            ENDIF
          ELSE
            RIFTR(L)=0.
            EVAPSW(L)=0.
            QSUMTMP(L)=MAX(QSUME(L),0.0)
          ENDIF
        ENDDO
        DO L=2,LA
          DIFQVOL=QSUME(L)-QSUMTMP(L)
          DO K=1,KC
            QSUM(L,K)=QSUM(L,K)-DIFQVOL*DZC(K)
          ENDDO
          QSUME(L)=QSUMTMP(L)
        ENDDO
      ENDIF
C
C **  SET OLD TIME LEVEL TERMS IN CONTINUITY EQUATION FOR
C **  NON BOUNDARY POINTS
C **  HRU=HMU*DYU/DXU & HRV=HMV*DXV/DYV
C **  DXYIP=1/(DXP*DYP)
C
      DO L=2,LA
        LN=LNC(L)
        FP(L)=DELTI*DXYP(L)*P1(L)-0.5*G*(UHDY1E(LEAST(L))-UHDY1E(L)
     &                                  +VHDX1E(LN )-VHDX1E(L))
      ENDDO
C
C **  ADVANCE EXTERNAL VARIABLES FOR THREE TIME LEVEL STEP
C
      IF(ISTL_.EQ.3)THEN
        DO L=2,LA
          UHDY2E(L)=UHDY1E(L)
          VHDX2E(L)=VHDX1E(L)
          UHDY1E(L)=UHDYE(L)
          VHDX1E(L)=VHDXE(L)
          U1V(L)=UV(L)
          V1U(L)=VU(L)
          P1(L)=P(L)
          H1U(L)=HU(L)
          H1V(L)=HV(L)
          H1UI(L)=HUI(L)
          H1VI(L)=HVI(L)
          H2P(L)=H1P(L)
          H1P(L)=HP(L)
        ENDDO
        IF(ISGWIE.GE.1)THEN
	      DO L=2,LA
	        AGWELV2(L)=AGWELV1(L)
	        AGWELV1(L)=AGWELV(L)
	      ENDDO
        ENDIF
      ENDIF
C
C **  SET NEW TIME LEVEL TERMS IN CONTINUITY EQUATION INCLUDING
C **  HOST-GUEST CHANNAL INTERACTION FOR NON BOUNDARY POINTS
C
      C1=0.5*G
      DO L=2,LA
        LN=LNC(L)
        ! ***  THE SUB & SVB SWITCHES ALREADY ACCOUNTED FOR
        FP(L)=FP(L)-C1*(FUHDYE(LEAST(L))-FUHDYE(L)
     &                 +FVHDXE(LN )-FVHDXE(L)
     &                -2.0*QSUME(L) )
      ENDDO
C
      IF(ISGWIE.GE.1)THEN
        DO L=2,LA
          FP(L)=FP(L)-G*SPB(L)*(RIFTR(L)+EVAPSW(L))
        ENDDO
      ENDIF
C
      IF(ISDSOLV.GE.1.AND.DEBUG)THEN
        OPEN(1,FILE='FP.OUT',POSITION='APPEND',STATUS='UNKNOWN')
        WRITE(1,1001)N,ISTL_
        DO L=2,LA
          WRITE(1,1001)IL(L),JL(L),FP1(L),FUHDYE(L),FUHDYE(LEAST(L)),
     &        FVHDXE(L),FVHDXE(LNC(L)),QSUME(L),RIFTR(L),EVAPSW(L)
        ENDDO
        CLOSE(1)
        IF(N.EQ.1)THEN
          OPEN(1,FILE='FP1.OUT',POSITION='APPEND',STATUS='UNKNOWN')
          WRITE(1,1001)N,ISTL_
          DO L=2,LA
            WRITE(1,1001)IL(L),JL(L),FP1(L),FUHDYE(L),FUHDYE(LEAST(L)),
     &          FVHDXE(L),FVHDXE(LNC(L)),QSUME(L),RIFTR(L),EVAPSW(L)
          ENDDO
          CLOSE(1)
        ENDIF
        IF(N.EQ.2)THEN
          OPEN(1,FILE='FP2.OUT',POSITION='APPEND',STATUS='UNKNOWN')
          WRITE(1,1001)N,ISTL_
          DO L=2,LA
            WRITE(1,1001)IL(L),JL(L),FP1(L),FUHDYE(L),FUHDYE(LEAST(L)),
     &          FVHDXE(L),FVHDXE(LNC(L)),QSUME(L),RIFTR(L),EVAPSW(L)
          ENDDO
          CLOSE(1)
        ENDIF
      ENDIF
C
      C1=-0.5*DELTD2*G
      DO L=2,LA
        CS(L)=C1*SVB(L  )*HRVO(L  )*RCY(L  )*HVTMP(L  )
        CW(L)=C1*SUB(L  )*HRUO(L  )*RCX(L  )*HUTMP(L  )
        CE(L)=C1*SUB(LEAST(L))*HRUO(LEAST(L))*RCX(LEAST(L))*HUTMP(LEAST(L))
      ENDDO
      DO L=2,LA
        LN=LNC(L)
        CN(L)=C1*SVB(LN )*HRVO(LN )*RCY(LN )*HVTMP(LN )
      ENDDO
C
C *** APPLY THE OPEN BOUNDARY CONDITIONS
C
      IF(NBCSOP.GT.0) CALL SETOPENBC(DELT,DELTD2,DELTI,HUTMP,HVTMP)
C
      ! *** SET THE CENTER
      DO L=2,LA
        CC(L)=DELTI*DXYP(L)-CS(L)-CW(L)-CE(L)-CN(L)
      ENDDO
C
      ! *** SCALE COEFFICIENTS IN EXTERNAL MODEL LINEAR EQUATION SYSTEM
      CCMNM=1.E+18
      DO L=2,LA
        CCMNM=MIN(CCMNM,CC(L))
        FPTMP(L)=FP(L)
      ENDDO
      CCMNMI=1./CCMNM
C
C *** APPLY THE OPEN BOUNDARY CONDITIONS FOR ADJACENT CELLS
C
      IF(NBCSOP.GT.0) CALL SETOPENBC2
C
      CC(1)=1.
      CC(LC)=1.
C
C **  SCALE BY MINIMUM DIAGONAL
C
      IF(IRVEC.EQ.9)THEN
        DO L=2,LA
          CCS(L)=CS(L)*CCMNMI
          CCW(L)=CW(L)*CCMNMI
          CCE(L)=CE(L)*CCMNMI
          CCN(L)=CN(L)*CCMNMI
          CCC(L)=CC(L)*CCMNMI
          FPTMP(L)=FPTMP(L)*CCMNMI
          CCCI(L)=1./CCC(L)
        ENDDO
      ENDIF
C
C **  SCALE TO NORMAL FORM
C
      IF(IRVEC.EQ.99)THEN
        DO L=2,LA
          CCS(L)=CS(L)/SQRT( CC(L)*CC(LSC(L)) )
          CCW(L)=CW(L)/SQRT( CC(L)*CC(LWEST(L)   ) )
          CCE(L)=CE(L)/SQRT( CC(L)*CC(LEAST(L)   ) )
          CCN(L)=CN(L)/SQRT( CC(L)*CC(LNC(L)) )
          CCC(L)=1.
          FPTMP(L)=FPTMP(L)/SQRT( CC(L) )
          P(L)=P(L)*SQRT( CC(L) )
          CCCI(L)=1.
        ENDDO
      ENDIF
C
      ! *** CALL THE PRECONDITIONED CONJUGATE GRADIENT SOLVER
      CALL CONGRAD (ISTL_)
C
      IF(IRVEC.EQ.99)THEN
        DO L=2,LA
          P(L)=P(L)/SQRT( CC(L) )
        ENDDO
      ENDIF
C
C ** DIAGNOSTICS
C
      IF(ISDSOLV.GE.1.AND.DEBUG)THEN
        OPEN(1,FILE='EQCOEF.OUT',POSITION='APPEND',STATUS='UNKNOWN')
        WRITE(1,1001)N,ISTL_
        DO L=2,LA
          SURFTMP=GI*P(L)
          WRITE(1,1001)IL(L),JL(L),CS(L),CW(L),CC(L),CE(L),CN(L),
     &        FP(L),SURFTMP
        ENDDO
        CLOSE(1)
        IF(N.EQ.1)THEN
          OPEN(1,FILE='EQCOEF1.OUT',POSITION='APPEND',STATUS='UNKNOWN')
          WRITE(1,1001)N,ISTL_
          DO L=2,LA
            SURFTMP=GI*P(L)
            WRITE(1,1001)IL(L),JL(L),CS(L),CW(L),CC(L),CE(L),CN(L),
     &          FP(L),SURFTMP
          ENDDO
          CLOSE(1)
        ENDIF
        IF(N.EQ.2)THEN
          OPEN(1,FILE='EQCOEF2.OUT',POSITION='APPEND',STATUS='UNKNOWN')
          WRITE(1,1001)N,ISTL_
          DO L=2,LA
            SURFTMP=GI*P(L)
            WRITE(1,1001)IL(L),JL(L),CS(L),CW(L),CC(L),CE(L),CN(L),
     &          FP(L),SURFTMP
          ENDDO
          CLOSE(1)
        ENDIF
      ENDIF
      IF(ISDSOLV.GE.1.AND.DEBUG)THEN
        OPEN(1,FILE='EQTERM.OUT',POSITION='APPEND',STATUS='UNKNOWN')
        WRITE(1,1001)N,ISTL_
        DO L=2,LA
          WRITE(1,1001)IL(L),JL(L),SUB(L),SVB(L),HRUO(L),
     &        HRVO(L),HUTMP(L),HVTMP(L)
        ENDDO
        CLOSE(1)
        IF(N.EQ.1)THEN
          OPEN(1,FILE='EQTERM1.OUT',POSITION='APPEND',STATUS='UNKNOWN')
          WRITE(1,1001)N,ISTL_
          DO L=2,LA
            WRITE(1,1001)IL(L),JL(L),SUB(L),SVB(L),HRUO(L),
     &          HRVO(L),HUTMP(L),HVTMP(L)
          ENDDO
          CLOSE(1)
        ENDIF
        IF(N.EQ.2)THEN
          OPEN(1,FILE='EQTERM2.OUT',POSITION='APPEND',STATUS='UNKNOWN')
          WRITE(1,1001)N,ISTL_
          DO L=2,LA
            WRITE(1,1001)IL(L),JL(L),SUB(L),SVB(L),HRUO(L),
     &          HRVO(L),HUTMP(L),HVTMP(L)
          ENDDO
          CLOSE(1)
        ENDIF
      ENDIF
 1001 FORMAT(2I5,10(1X,E12.4))
 1002 FORMAT(3I4,10(1X,E9.2))
C
C **  CALCULATE UHEX AND VHEX AND TOTAL DEPTHS AT TIME LEVEL (N+1)
C **  HRU=SUB*DYU/DXU & HRV=SVB*DXV/DYV
C
      DO L=2,LA
        LS=LSC(L)
        UHDYE(L)=SUB(L)*( FUHDYE(L)
     &      -DELTD2*HRUO(L)*RCX(L)*HUTMP(L)*(P(L)-P(LWEST(L))) )
        VHDXE(L)=SVB(L)*( FVHDXE(L)
     &      -DELTD2*HRVO(L)*RCY(L)*HVTMP(L)*(P(L)-P(LS )) )
      ENDDO
C
      DO L=2,LA
        UHE(L)=UHDYE(L)*DYIU(L)
        VHE(L)=VHDXE(L)*DXIV(L)
      ENDDO
C
C **  CALCULATE REVISED CELL DEPTHS BASED ON NEW HORIZONTAL
C **  TRANSPORTS AT (N+1)
C
      IF(ISTL_.EQ.3)THEN
        DO L=2,LA
          LN=LNC(L)
          HP(L)=H2P(L)+DELT*DXYIP(L)*(QSUME(L)
     &        -0.5*(UHDYE(LEAST(L))+UHDY2E(LEAST(L))-UHDYE(L)-UHDY2E(L)
     &             +VHDXE(LN )+VHDX2E(LN )-VHDXE(L)-VHDX2E(L)))
        ENDDO
      ELSE
        DO L=2,LA
          LN=LNC(L)
          HP(L)=H1P(L)+DELT*DXYIP(L)*(QSUME(L)
     &        -0.5*(UHDYE(LEAST(L))+UHDY1E(LEAST(L))-UHDYE(L)-UHDY1E(L)
     &             +VHDXE(LN )+VHDX1E(LN )-VHDXE(L)-VHDX1E(L)))
        ENDDO
      ENDIF
C
      IF(ISGWIE.GE.1)THEN
        DO L=2,LA
          HP(L)=HP(L)-DELT*DXYIP(L)*(RIFTR(L)+EVAPSW(L))
        ENDDO
      ENDIF
C
      ! *** APPLY OPEN BOUNDARYS
      DO LL=1,NBCSOP
        L=LOBCS(LL)
        HP(L)=GI*P(L)-BELV(L)
      ENDDO
C
C **  PERFORM FINAL UPDATES OF P,HU, AND HV
C
      DO L=2,LA
        P(L)=G*(HP(L)+BELV(L))
      ENDDO
      DO L=2,LA
        LS=LSC(L)
        HU(L)=0.5*(DXYP(L)*HP(L)+DXYP(LWEST(L))*HP(LWEST(L)))*DXYIU(L)
        HV(L)=0.5*(DXYP(L)*HP(L)+DXYP(LS )*HP(LS ))*DXYIV(L)
      ENDDO
      DO L=2,LA
        HPI(L)=1./HP(L)
        HUI(L)=1./HU(L)
        HVI(L)=1./HV(L)
      ENDDO
C
C **  PERFORM UPDATE ON GROUNDWATER ELEVATION
C
      IF(ISGWIE.GE.1)THEN
        DO L=2,LA
          QSUM(L,KC)=QSUM(L,KC)-EVAPSW(L)
          QSUM(L,1 )=QSUM(L,1 )-RIFTR(L)
        ENDDO
C
C       INFILTRATION STEP
C
        RNPORI=1./RNPOR
        IF(ISTL_.EQ.3)THEN
          DO L=2,LA
            AGWELV(L)=AGWELV2(L)+RNPORI*DELT*DXYIP(L)*RIFTR(L)
          ENDDO
        ELSE
          DO L=2,LA
            AGWELV(L)=AGWELV1(L)+RNPORI*DELT*DXYIP(L)*RIFTR(L)
          ENDDO
        ENDIF
        DO L=2,LA
          AGWELV(L)=MIN(AGWELV(L),BELV(L))
        ENDDO
C
C       ET STEP
C
        DO L=2,LA
          IF(EVAPCVT.LT.0.)THEN
            SVPW=(10.**((0.7859+0.03477*TEM(L,KC))/
     &          (1.+0.00412*TEM(L,KC))))
           EVAPT(L)=CLEVAP(L)*0.7464E-3*WINDST(L)*(SVPW-VPA(L))/PATMT(L)
          ENDIF
          ETGWTMP=EVAPT(L)-EVAPSW(L)*DXYIP(L)
          ETGWTMP=MAX(ETGWTMP,0.0)
          ETGWAVL=RNPOR*DELTI*(AGWELV(L)-BELAGW(L))
          ETGWAVL=MAX(ETGWAVL,0.0)
          ETGWTMP=MIN(ETGWTMP,ETGWAVL)
          EVAPGW(L)=ETGWTMP*DXYP(L)
        ENDDO
        DO L=2,LA
          AGWELV(L)=AGWELV(L)-RNPORI*DELT*DXYIP(L)*EVAPGW(L)
        ENDDO
        DO L=2,LA
          AGWELV(L)=MAX(AGWELV(L),BELAGW(L))
        ENDDO
      ENDIF
C
C **  CHECK FOR NEGATIVE DEPTHS
C
      IF(ISNEGH.GE.1)THEN
        INEGFLG=0
        DO L=2,LA
          IF(HP(L).LT.0.)THEN
            INEGFLG=1
            LN=LNC(L)
            WRITE (6,6060)IL(L),JL(L),HP(L),H1P(L),H2P(L)
            WRITE (6,6061)IL(L),JL(L),HU(L),H1U(L)
            WRITE (6,6062)IL(L),JL(L),HU(LEAST(L)),H1U(LEAST(L))
            WRITE (6,6063)IL(L),JL(L),HV(L),H1V(L)
            WRITE (6,6064)IL(L),JL(L),HV(LN),H1V(LN)
            WRITE (8,6060)IL(L),JL(L),HP(L),H1P(L),H2P(L)
            WRITE (8,6061)IL(L),JL(L),HU(L),H1U(L)
            WRITE (8,6062)IL(L),JL(L),HU(LEAST(L)),H1U(LEAST(L))
            WRITE (8,6063)IL(L),JL(L),HV(L),H1V(L)
            WRITE (8,6064)IL(L),JL(L),HV(LN),H1V(LN)
          ENDIF
        ENDDO
        DO L=2,LA
          IF(HU(L).LT.0.)THEN
            INEGFLG=1
            LN=LNC(L)
            WRITE (6,6060)IL(L),JL(L),HP(L),H1P(L),H2P(L)
            WRITE (6,6061)IL(L),JL(L),HU(L),H1U(L)
            WRITE (6,6062)IL(L),JL(L),HU(LEAST(L)),H1U(LEAST(L))
            WRITE (6,6063)IL(L),JL(L),HV(L),H1V(L)
            WRITE (6,6064)IL(L),JL(L),HV(LN),H1V(LN)
            WRITE (8,6060)IL(L),JL(L),HP(L),H1P(L),H2P(L)
            WRITE (8,6061)IL(L),JL(L),HU(L),H1U(L)
            WRITE (8,6062)IL(L),JL(L),HU(LEAST(L)),H1U(LEAST(L))
            WRITE (8,6063)IL(L),JL(L),HV(L),H1V(L)
            WRITE (8,6064)IL(L),JL(L),HV(LN),H1V(LN)
          ENDIF
        ENDDO
        DO L=2,LA
          IF(HV(L).LT.0.)THEN
            INEGFLG=1
            LN=LNC(L)
            WRITE (6,6060)IL(L),JL(L),HP(L),H1P(L),H2P(L)
            WRITE (6,6061)IL(L),JL(L),HU(L),H1U(L)
            WRITE (6,6062)IL(L),JL(L),HU(LEAST(L)),H1U(LEAST(L))
            WRITE (6,6063)IL(L),JL(L),HV(L),H1V(L)
            WRITE (6,6064)IL(L),JL(L),HV(LN),H1V(LN)
            WRITE (8,6060)IL(L),JL(L),HP(L),H1P(L),H2P(L)
            WRITE (8,6061)IL(L),JL(L),HU(L),H1U(L)
            WRITE (8,6062)IL(L),JL(L),HU(LEAST(L)),H1U(LEAST(L))
            WRITE (8,6063)IL(L),JL(L),HV(L),H1V(L)
            WRITE (8,6064)IL(L),JL(L),HV(LN),H1V(LN)
          ENDIF
        ENDDO
      ENDIF
      IF(ISNEGH.EQ.2)THEN
        IF(INEGFLG.EQ.1)THEN
          CALL RESTOUT(1)
          OPEN(1,FILE='EQCOEF.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
          OPEN(1,FILE='EQCOEF.OUT',POSITION='APPEND',STATUS='UNKNOWN')
          WRITE(1,1001)N,ISTL_
          DO L=2,LA
            SURFTMP=GI*P(L)
            WRITE(1,1001)IL(L),JL(L),CCS(L),CCW(L),CCC(L),CCE(L),CCN(L),
     &          FPTMP(L),SURFTMP
          ENDDO
          CLOSE(1)
          OPEN(1,FILE='EQTERM.OUT',STATUS='UNKNOWN')
          CLOSE(1,STATUS='DELETE')
          OPEN(1,FILE='EQTERM.OUT',POSITION='APPEND',STATUS='UNKNOWN')
          WRITE(1,1001)N,ISTL_
          DO L=2,LA
            WRITE(1,1002)IL(L),JL(L),SUB(L),SVB(L),HRUO(L),
     &          HRVO(L),HUTMP(L),HVTMP(L)
          ENDDO
          CLOSE(1)
          OPEN(1,FILE='CFLMAX.OUT')
          CLOSE(1,STATUS='DELETE')
          OPEN(1,FILE='CFLMAX.OUT')
          DO L=2,LA
            WRITE(1,1991)IL(L),JL(L),(CFLUUU(L,K),K=1,KC)
            WRITE(1,1992)(CFLVVV(L,K),K=1,KC)
            WRITE(1,1992)(CFLWWW(L,K),K=1,KC)
            WRITE(1,1992)(CFLCAC(L,K),K=1,KC)
          ENDDO
          CLOSE(1)
          STOP
        ENDIF
      ENDIF
 1991 FORMAT(2I5,12F7.2)
 1992 FORMAT(10X,12F7.2)
 6060 FORMAT('  NEG DEPTH AT I,J =',2I4,'  HP,H1P,H2P =',3(2X,E12.4))
 6061 FORMAT('  NEG DEPTH AT I,J =',2I4,'  HUW,H1UW =',2(2X,E12.4))
 6062 FORMAT('  NEG DEPTH AT I,J =',2I4,'  HUE,H1UE =',2(2X,E12.4))
 6063 FORMAT('  NEG DEPTH AT I,J =',2I4,'  HVS,H1VS =',2(2X,E12.4))
 6064 FORMAT('  NEG DEPTH AT I,J =',2I4,'  HVN,H1VN =',2(2X,E12.4))
 6910 FORMAT('  DRYING AT N,I,J =',I10,2I6,'  HP,H1P,H2P ='
     &    ,3(2X,E12.4))
 6911 FORMAT('  DRY W FACE N,I,J =',I10,2I6,' HU,H,H1 =',3(2X,E12.4))
 6912 FORMAT('  DRY E FACE N,I,J =',I10,2I6,' HU,H,H1 =',3(2X,E12.4))
 6913 FORMAT('  DRY S FACE N,I,J =',I10,2I6,' HV,H,H1 =',3(2X,E12.4))
 6914 FORMAT('  DRY N FACE N,I,J =',I10,2I6,' HV,H,H1 =',3(2X,E12.4))
 6920 FORMAT('  WETTING AT N,I,J =',I10,2I6,' HP,H1P,H2P ='
     &    ,3(2X,E12.4))
 6921 FORMAT('  WET S FACE N,I,J =',I10,2I6,' HV,H,H1 =',3(2X,E12.4))
 6922 FORMAT('  WET W FACE N,I,J =',I10,2I6,' HU,H,H1 =',3(2X,E12.4))
 6923 FORMAT('  WET E FACE N,I,J =',I10,2I6,' HU,H,H1 =',3(2X,E12.4))
 6924 FORMAT('  WET N FACE N,I,J =',I10,2I6,' HV,H,H1 =',3(2X,E12.4))
 6930 FORMAT('  WET BY VOL  N,I,J =',I10,2I6,' HP,H1P,H2P ='
     &    ,3(2X,E12.4))
 6940 FORMAT('  RESOLVE,  N,I,J =',I10,2I6,' HP,H1P,H2P ='
     &    ,3(2X,E12.4))
 6941 FORMAT('  RESOLVE,  N,I,J =',I10,2I6,' HUE,HP,H1P ='
     &    ,3(2X,E12.4))
 6942 FORMAT('  RESOLVE,  N,I,J =',I10,2I6,' HUW,HP,H1P ='
     &    ,3(2X,E12.4))
 6943 FORMAT('  RESOLVE,  N,I,J =',I10,2I6,' HVS,HP,H1P ='
     &    ,3(2X,E12.4))
 6944 FORMAT('  RESOLVE,  N,I,J =',I10,2I6,' HVN,HP,H1P ='
     &    ,3(2X,E12.4))
 6945 FORMAT('  RESOLVE NEG,  N,I,J =',I10,2I6,' HP,H1P,H2P ='
     &    ,3(2X,E12.4))
 6950 FORMAT('  RESOLVE, NEG DEP N,I,J =',I10,2I6,' HP,H1P,H2P ='
     &    ,3(2X,E12.4))
C
C **  CALCULATE THE EXTERNAL DIVERGENCE
C
      IF(ISDIVEX.EQ.1)THEN
        DIVEXMX=0.
        DIVEXMN=1000000.
        IF(ISTL_.EQ.3)THEN
          DO L=2,LA
            IF(SPB(L).NE.0)THEN
              LN=LNC(L)
              DIVEX=SPB(L)*(DXYP(L)*(HP(L)-H2P(L))*DELTI
     &            +0.5*(UHDYE(LEAST(L))+UHDY2E(LEAST(L))-UHDYE(L)-UHDY2E(L)
     &            +VHDXE(LN)+VHDX2E(LN)-VHDXE(L)-VHDX2E(L))-QSUME(L)
     &            +RIFTR(L)+EVAPSW(L))
              IF(DIVEX.GT.DIVEXMX)THEN
                DIVEXMX=DIVEX
                LMAX=L
              ENDIF
              IF(DIVEX.LT.DIVEXMN)THEN
                DIVEXMN=DIVEX
                LMIN=L
              ENDIF
            ENDIF
          ENDDO
        ELSE
          DO L=2,LA
            IF(SPB(L).NE.0)THEN
              LN=LNC(L)
              DIVEX=SPB(L)*(DXYP(L)*(HP(L)-H1P(L))*DELTI
     &            +0.5*(UHDYE(LEAST(L))+UHDY1E(LEAST(L))-UHDYE(L)-UHDY1E(L)
     &            +VHDXE(LN)+VHDX1E(LN)-VHDXE(L)-VHDX1E(L))-QSUME(L)
     &            +RIFTR(L)+EVAPSW(L))
              IF(DIVEX.GT.DIVEXMX)THEN
                DIVEXMX=DIVEX
                LMAX=L
              ENDIF
              IF(DIVEX.LT.DIVEXMN)THEN
                DIVEXMN=DIVEX
                LMIN=L
              ENDIF
            ENDIF
          ENDDO
        ENDIF
        IMAX=IL(LMAX)
        JMAX=JL(LMAX)
        IMIN=IL(LMIN)
        JMIN=JL(LMIN)
        WRITE(6,6628)DIVEXMX,IMAX,JMAX
        WRITE(6,6629)DIVEXMN,IMIN,JMIN
      ENDIF
  566 FORMAT('  I=',I5,3X,'J=',I5,3X,'HP=',F12.4)
 6628 FORMAT('  DIVEXMX=',E13.5,5X,2I10)
 6629 FORMAT('  DIVEXMN=',E13.5,5X,2I10)
C
C **  UPDATE ZERO DIMENSION VOLUME BALANCE
C
      IF(ISDRY.GE.1.AND.ISTL_.EQ.3)THEN
        VOLADD=0.
        DO L=2,LA
          IF(SPB(L).NE.0)THEN
            VOLADD=VOLADD+QSUME(L)-RIFTR(L)-EVAPSW(L)
          ENDIF
        ENDDO
        VOLADD=VOLADD*DT
        VOLZERD=VOLZERD+VOLADD
        VETZERD=VETZERD+VOLADD+DT*EVAPSW(L)
      ENDIF
 5303 FORMAT(2X,F10.4,2X,F10.5,3(2X,E13.5))
      RETURN
      END


text:
      SUBROUTINE STRAX1
C
C     THIS ROUTINE IS PHASE I OF STRESS DATA RECOVERY FOR THE AXI-
C     SYMMETRIC RING WITH TRIANGULAR CROSS SECTION RING
C
C     ECPT ( 1) = ELEMENT ID
C     ECPT ( 2) = SIL A                                   I
C     ECPT ( 3) = SIL B                                   I
C     ECPT ( 4) = SIL C                                   I
C     ECPT ( 5) = MATERIAL ORIENTATION ANGLE(DEGREES)     R
C     ECPT ( 7) = MATERIAL ID                             I
C     ECPT ( 8) TO ECPT(21) = STRESS PHASE ANG.           R
C     ECPT (22) = CORD. SYS. GRID POINT A (NOT USED)      I
C     ECPT (23) = R-CORD OF GRID A                        R
C     ECPT (24) = Z-CORD OF GRID A                        R
C     ECPT (25) = 0.0                                     R
C     ECPT (26) = CORD. SYS. GRID POINT B (NOT USED)      I
C     ECPT (27) = R-CORD OF GRID B                        R
C     ECPT (28) = Z-CORD OF GRID B                        R
C     ECPT (29) = 0.0                                     R
C     ECPT (30) = CORD. SYS. GRID POINT C (NOT USED)      I
C     ECPT (31) = R-CORD OF GRID C                        R
C     ECPT (32) = Z-CORD OF GRID C                        R
C     ECPT (33) = 0.0                                     R
C     ECPT (34) = EL. TEMPERATURE FOR MATERIAL PROP       R
C
C     ANY GROUP OF STATEMENTS PREFACED BY AN IF STATEMENT CONTAINING
C     ...KSYS78 OR LSYS78 ...  INDICATES CODING NECESSARY FOR THIS
C     ELEMENT*S PIEZOELECTRIC CAPABILITY
C
C     KSYS78 = 0   ELASTIC, NON-PIEZOELECTRIC MATERIAL
C     KSYS78 = 1   ELECTRICAL-ELASTIC COUPLED, PIEZOELETRIC MATERIAL
C     KSYS78 = 2   ELASTIC ONLY, PIEZOELECTRIC MATERIAL
C     LSYS78 = .TRUE. IF KSYS78 = 0, OR 2
C
      LOGICAL         PZMAT,LSYS78
      DIMENSION       IECPT(35),DELINT(12),TEO(45),ACURL(117),R(3),
     1                Z(3),ICS(3),D1(27),D2(9),ACURP1(27),ACURP2(9),
     2                GABABP(3,3),EE(63),WJP(3,3)
C     ECPT COMMON BLOCK
      COMMON /SDR2X5/ ECPT(34),DUM5(66),IDEL,IGP(3),TZ,SEL(54),TS(6),
     1                AK(81),PHI(14),SELP1(18),AKPH2(9),AKUPH(27),
     2                SELP2(27),SELP3(9)
      COMMON /SDR2X6/ GABABQ(9,9), D(81), WJ(6,9)
      COMMON /MATIN / MATIDC,MATFLG,ELTEMP,STRESS,SINTH,COSTH
      COMMON /MATOUT/ E(3),ANU(3),RHO,G(3),ALF(3),TZERO,GSUBE,MOSKP(9),
     1                SETMAT
      COMMON /MATPZ / PZOUT(51)
C     COMMON /MATPZ / CE11,CE12,CE13,CE14,CE15,CE16,CE22,CE23,CE24,CE25,
C                     CE26,CE33,CE34,CE35,CE36,CE44,CE45,CE46,CE55,CE56,
C                     CE66,E11,E12,E13,E14,E15,E16,E21,E22,E23,E24,E25,
C                     E26,E31,E32,E33,E34,E35,E36,EPS11,EPS12,EPS13,EPS2
C                     EPS23,EPS33,RHO,A1,A2,A12,TREF,GE
C
      COMMON /CONDAS/ CONSTS(5)
      COMMON /SYSTEM/ KSYSTM(77),KSYS78
      EQUIVALENCE     (IECPT(1),ECPT(1)), (Z(1),Z1), (Z(2),Z2),
     1                (R(1),R1), (R(2),R2), (R(3),R3), (Z(3),Z3),
     2                (CONSTS(1),PI), (CONSTS(4),DEGRAD),
     3                (ACURP1(1),ACURL(82)), (ACURP2(1),ACURL(109))
C
      LSYS78 = .FALSE.
      IF (KSYS78.EQ.0 .OR. KSYS78.EQ.2) LSYS78 = .TRUE.
C
C     START EXECUTION
C
C     STORE ECPT PARAMETERS IN LOCAL VARIABLES
C
      IDEL   = IECPT( 1)
      IGP(1) = IECPT( 2)
      IGP(2) = IECPT( 3)
      IGP(3) = IECPT( 4)
      MATID  = IECPT( 7)
      ICS(1) = IECPT(22)
      R(1)   = ECPT( 23)
      Z(1)   = ECPT( 24)
      D(1)   = ECPT( 25)
      ICS(2) = IECPT(26)
      R(2)   = ECPT( 27)
      Z(2)   = ECPT( 28)
      D(2)   = ECPT( 29)
      ICS(3) = IECPT(30)
      R(3)   = ECPT( 31)
      Z(3)   = ECPT( 32)
      D(3)   = ECPT( 33)
      DGAMA  = ECPT(  5)
      TEMPE  = ECPT( 34)
C
C     TEST THE VALIDITY OF THE GRID POINT COORDINATES
C
      DO 200 I = 1,3
      IF (R(I) .LE. 0.0) GO TO 910
      IF (D(I) .NE. 0.0) GO TO 910
  200 CONTINUE
C
C     COMPUTE THE ELEMENT COORDINATES
C
      ZMIN = AMIN1(Z1,Z2,Z3)
      Z1 = Z1 - ZMIN
      Z2 = Z2 - ZMIN
      Z3 = Z3 - ZMIN
C
C     FORM THE TRANSFORMATION MATRIX GABABQ (9X9) FROM FIELD COORDINATES
C     TO GRID POINT DEGREES OF FREEDOM
C
      DO 300 I = 1,81
  300 GABABQ(I,1) = 0.0
      AA = R2*Z3 + R1*Z2 + Z1*R3 - Z2*R3 - R1*Z3 - R2*Z1
      AA = 1.0/AA
      C1 = AA*(R2*Z3 - Z2*R3)
      C2 =-AA*(Z3 - Z2)
      C3 = AA*(R3 - R2)
      GABABQ(1,1) = C1
      GABABQ(1,2) = C2
      GABABQ(1,3) = C3
      GABABQ(2,4) = C1
      GABABQ(2,5) = C2
      GABABQ(2,6) = C3
      GABABQ(3,7) = C1
      GABABQ(3,8) = C2
      GABABQ(3,9) = C3
      IF (LSYS78) GO TO 302
      GABABP(1,1) = C1
      GABABP(1,2) = C2
      GABABP(1,3) = C3
  302 CONTINUE
      C1 =-AA*(R1*Z3 - Z1*R3)
      C2 = AA*(Z3 - Z1)
      C3 =-AA*(R3 - R1)
      GABABQ(4,1) = C1
      GABABQ(4,2) = C2
      GABABQ(4,3) = C3
      GABABQ(5,4) = C1
      GABABQ(5,5) = C2
      GABABQ(5,6) = C3
      GABABQ(6,7) = C1
      GABABQ(6,8) = C2
      GABABQ(6,9) = C3
      IF (LSYS78) GO TO 304
      GABABP(2,1) = C1
      GABABP(2,2) = C2
      GABABP(2,3) = C3
  304 CONTINUE
      C1 = AA*(R1*Z2 - Z1*R2)
      C2 =-AA*(Z2 - Z1)
      C3 = AA*(R2 - R1)
      GABABQ(7,1) = C1
      GABABQ(7,2) = C2
      GABABQ(7,3) = C3
      GABABQ(8,4) = C1
      GABABQ(8,5) = C2
      GABABQ(8,6) = C3
      GABABQ(9,7) = C1
      GABABQ(9,8) = C2
      GABABQ(9,9) = C3
      IF (LSYS78) GO TO 306
      GABABP(3,1) = C1
      GABABP(3,2) = C2
      GABABP(3,3) = C3
  306 CONTINUE
C
C     COMPUTE THE INTEGRAL VALUES IN ARRAY DELINT THE ORDER IS INDICATED
C     BY THE FOLLOWING TABLE
C
C       DELINT(01) = (-1,0)
C       DELINT(02) = (-1,1)
C       DELINT(03) = (-1,2)
C       DELINT(04) = (0, 0)
C       DELINT(05) = (0, 1)
C       DELINT(06) = (1, 0)
C
      RA  = (R1 + R2 + R3)/3.0
      ZA  = (Z1 + Z2 + Z3)/3.0
      RH  = AMIN1(R1,R2,R3)/10.0
      DR  = AMAX1(ABS(R1-R2),ABS(R2-R3),ABS(R3-R1))
      AREA= (R1*(Z2-Z3) + R2*(Z3-Z1) + R3*(Z1-Z2))/2.0
      I1  = 0
      DO 400 I = 1,2
      IP  = I - 2
      DO 350 J = 1,3
      IQ  = J  - 1
      I1  = I1 + 1
      IF (I1 .NE. 6) GO TO 310
      IP  = 1
      IQ  = 0
  310 IF (DR .GT. RH) GO TO 320
      DELINT(I1) = ((RA**IP)*(ZA**IQ))*AREA
      GO TO 330
  320 DELINT(I1) = AIS(3,IP,IQ,R,Z)
  330 DELINT(I1) = ABS(DELINT(I1))
  350 CONTINUE
  400 CONTINUE
C
C     LOCATE THE MATERIAL PROPERTIES IN THE MAT1 OR MAT3
C
      DGAMR  = DGAMA*DEGRAD
      COSG   = COS(DGAMR)
      SING   = SIN(DGAMR)
      SINTH  = SING
      COSTH  = COSG
      MATIDC = MATID
      MATFLG = 7
      IF (KSYS78 .GT. 0) MATFLG = 9
      ELTEMP = TEMPE
      CALL MAT (IDEL)
      PZMAT  = .FALSE.
      IF (SETMAT.EQ.4. .OR. SETMAT.EQ.5.) PZMAT = .TRUE.
      IF (PZMAT) GO TO 410
      KSAVE  = KSYS78
      KSYS78 = 0
      LSYS78 = .TRUE.
      GO TO 420
  410 RHO    = PZOUT(46)
      ALF(1) = PZOUT(47)
      ALF(2) = PZOUT(48)
      ALF(3) = PZOUT(49)
      TZERO  = PZOUT(50)
      GSUBE  = PZOUT(51)
  420 CONTINUE
      IF (SETMAT .EQ. 2.0) GO TO 920
      TZ = TZERO
      IF (KSYS78 .GT. 0) GO TO 500
C
C     SET MATERIAL PROPERTIES IN LOCAL VARIABLES (AGAIN)
C
      ER  = E(1)
      ET  = E(2)
      EZ  = E(3)
      VRO = ANU(1)
      VOZ = ANU(2)
      VZR = ANU(3)
      GOR = G(1)
      GZO = G(2)
      GRZ = G(3)
      VOR = VRO*ET/ER
      VZO = VOZ*EZ/ET
      VRZ = VZR*ER/EZ
      DEL = 1.0/(1.0-VRO*VOR-VOZ*VZO-VZR*VRZ-VRO*VOZ*VZR-VRZ*VOR*VZO)
C
C     COMPUTE ELASTIC CONSTANTS MATRIX FROM MATERIAL TO ELEMENT AXIS
C
  500 CONTINUE
      DO 510 I = 1,45
  510 TEO (I) = 0.0
      IF (KSYS78 .GT. 0) GO TO 512
      TEO ( 1) = ER*(1.0 - VOZ*VZO)*DEL
      TEO ( 2) = ER*(VZR + VZO*VOR)*DEL
      TEO ( 3) = EZ*(1.0 - VRO*VOR)*DEL
      TEO ( 4) = ER*(VOR + VZR*VOZ)*DEL
      TEO ( 5) = ET*(VZO + VRO*VZR)*DEL
      TEO ( 6) = ET*(1.0 - VRZ*VZR)*DEL
      TEO (10) = GRZ
      TEO (15) = GOR
      TEO (21) = GZO
      GO TO 514
  512 CONTINUE
C
C     PIEZOELECTRIC MATERIAL PROPERTIES STORED IN TEO(22-39)
C     DIELECTRIC MATERIAL PROPERTIES STORED IN TEO(40-45)
C     TEO(22-39) CONTAINS E-TRANSPOSE
C
      TEO( 1) = PZOUT( 1)
      TEO( 2) = PZOUT( 2)
      TEO( 3) = PZOUT( 7)
      TEO( 4) = PZOUT( 3)
      TEO( 5) = PZOUT( 8)
      TEO( 6) = PZOUT(12)
      TEO( 7) = PZOUT( 4)
      TEO( 8) = PZOUT( 9)
      TEO( 9) = PZOUT(13)
      TEO(10) = PZOUT(16)
      TEO(11) = PZOUT( 5)
      TEO(12) = PZOUT(10)
      TEO(13) = PZOUT(14)
      TEO(14) = PZOUT(17)
      TEO(15) = PZOUT(19)
      TEO(16) = PZOUT( 6)
      TEO(17) = PZOUT(11)
      TEO(18) = PZOUT(15)
      TEO(19) = PZOUT(18)
      TEO(20) = PZOUT(20)
      TEO(21) = PZOUT(21)
C
      IF (KSYS78 .EQ. 2) GO TO 514
      TEO(22) = PZOUT(22)
      TEO(23) = PZOUT(28)
      TEO(24) = PZOUT(34)
      TEO(25) = PZOUT(23)
      TEO(26) = PZOUT(29)
      TEO(27) = PZOUT(35)
      TEO(28) = PZOUT(24)
      TEO(29) = PZOUT(30)
      TEO(30) = PZOUT(36)
      TEO(31) = PZOUT(25)
      TEO(32) = PZOUT(31)
      TEO(33) = PZOUT(37)
      TEO(34) = PZOUT(26)
      TEO(35) = PZOUT(32)
      TEO(36) = PZOUT(38)
      TEO(37) = PZOUT(27)
      TEO(38) = PZOUT(33)
      TEO(39) = PZOUT(39)
      TEO(40) =-PZOUT(40)
      TEO(41) =-PZOUT(41)
      TEO(42) =-PZOUT(42)
      TEO(43) =-PZOUT(43)
      TEO(44) =-PZOUT(44)
      TEO(45) =-PZOUT(45)
  514 CONTINUE
C
      DO 520 I = 5,63
  520 EE(I) = 0.0
      C2  = COSG*COSG
      C4  = C2  *C2
      S2  = SING*SING
      S4  = S2  *S2
      C2S2= C2  *S2
      C3  = COSG*C2
      S3  = SING*S2
      CS2 = COSG*S2
      SC2 = SING*C2
      CS  = COSG*SING
C
      EE( 1) = TEO(1)*C4 + TEO(3)*S4 + 2.0*C2S2*(TEO(2) + 2.0*TEO(10))
      EE( 2) = TEO(2)*(C4+S4) + C2S2*(TEO(1) + TEO(3)-4.0*TEO(10))
      EE( 3) = TEO(4)*C2 + TEO(5)*S2
      EE( 4) = COSG*SING*S2*(TEO(2)-TEO(3) + 2.0*TEO(10))
     4       + SING*COSG*C2*(TEO(1)-TEO(2) - 2.0*TEO(10))
      EE( 7) = EE(2)
      EE( 8) = TEO(1)*S4 + 2.0*C2S2*(TEO(2) + 2.0*TEO(10)) + TEO(3)*C4
      EE( 9) = TEO(4)*S2 + TEO(5)*C2
      EE(10) = SING*COSG*C2 * (TEO(2)-TEO(3) + 2.0*TEO(10))
     O       + COSG*SING*S2 * (TEO(1)-TEO(2) - 2.0*TEO(10))
      EE(13) = EE(3)
      EE(14) = EE(9)
      EE(15) = TEO(6)
      EE(16) = SING*COSG*(TEO(4)-TEO(5))
      EE(19) = EE(4)
      EE(20) = EE(10)
      EE(21) = EE(16)
      EE(22) = C2S2*(TEO(1) - 2.0*TEO(2) + TEO(3)) + TEO(10)*(C2-S2)**2
      EE(29) = TEO(15)*C2 + TEO(21)*S2
      EE(30) = SING*COSG*(TEO(15)-TEO(21))
      EE(35) = EE(30)
      EE(36) = TEO(15)*S2 + TEO(21)*C2
      IF (LSYS78) GO TO 530
C
C     PIEZOELECTRIC MATERIAL PROPERTIES IN ELEMENT COORDINATES
C
      EE(37) = C3*TEO(22) - S3*TEO(26) + CS2*(TEO(25)+2.0*TEO(32))
     7       - SC2*(TEO(23)+2.0*TEO(31))
      EE(38) = C3*TEO(23) + S3*TEO(25) + CS2*(TEO(26)-2.0*TEO(31))
     8       + SC2*(TEO(22)-2.0*TEO(32))
      EE(39) = S2*TEO(27) + C2*TEO(24) - 2.0*CS*TEO(33)
      EE(40) = C3*TEO(25) - S3*TEO(23) + CS2*(TEO(22)-2.0*TEO(32))
     O       - SC2*(TEO(26)-2.0*TEO(31))
      EE(41) = C3*TEO(26) + S3*TEO(22) + CS2*(TEO(23)+2.0*TEO(31))
     1       + SC2*( TEO(25)+2.0*TEO(32))
      EE(42) = S2*TEO(24) + C2*TEO(27) + 2.0*CS*TEO(33)
      EE(43) = COSG*TEO(28) - SING*TEO(29)
      EE(44) = COSG*TEO(29) + SING*TEO(28)
      EE(45) = TEO(30)
      EE(46) = C3*TEO(31) + S3*TEO(32) - CS2*(TEO(23)-TEO(26)+TEO(31))
     6       + SC2*(-TEO(32)-TEO(25)+TEO(22))
      EE(47) = C3*TEO(32) - S3*TEO(31) - CS2*(TEO(25)-TEO(22)+TEO(32))
     7       + SC2*(TEO(23)+TEO(31)-TEO(26))
      EE(48) = (C2-S2)*TEO(33) + CS*(TEO(24)-TEO(27))
      EE(49) = C2*TEO(34) + S2*TEO(38) - CS*(TEO(35)+TEO(37))
      EE(50) = C2*TEO(35) - S2*TEO(37) + CS*(TEO(34)-TEO(38))
      EE(51) = COSG*TEO(36) - SING*TEO(39)
      EE(52) = C2*TEO(37) - S2*TEO(35) - CS*(TEO(38)-TEO(34))
      EE(53) = C2*TEO(38) + S2*TEO(34) + CS*(TEO(35)+TEO(37))
      EE(54) = COSG*TEO(39) + SING*TEO(36)
C
C     DIELECTRIC MATERIAL PROPERTIES IN ELEMENT COORDINTES
C
      EE(55) = S2*TEO(43) - 2.0*CS*TEO(41) + C2*TEO(40)
      EE(56) = (C2-S2)*TEO(41) - CS*(TEO(43)-TEO(40))
      EE(57) =-SING*TEO(44) + COSG*TEO(42)
      EE(59) = C2*TEO(43) + 2.0*CS*TEO(41) + S2*TEO(40)
      EE(60) = COSG*TEO(44) + SING*TEO(42)
      EE(63) = TEO(45)
  530 CONTINUE
C
      IECPT(1) = IECPT(1) - (IECPT(1)/1000)*1000 - 1
      AJHO  = IECPT(1)
      AJJHO = AJHO*AJHO
C
C     FORM THE ELEMENT STIFFNESS MATRIX IN FIELD SYSTEM
C
      ACURL( 1) = (EE(15) + AJJHO*EE(29))*DELINT(1)
      ACURL( 2) = (EE( 3) + EE(15) + AJJHO*EE(29))*DELINT(4)
      ACURL( 3) = (EE(15) + AJJHO*EE(29))*DELINT(2) + EE(16)*DELINT(4)
      ACURL( 4) = (EE(15) + EE(29))*AJHO*DELINT(1)
      ACURL( 5) = EE(15)*AJHO*DELINT(4)
      ACURL( 6) = (EE(15)+EE(29))*AJHO*DELINT(2) - EE(30)*AJHO*DELINT(4)
      ACURL( 7) = AJJHO*DELINT(1)*EE(35)
      ACURL( 8) = (EE(16) + AJJHO*EE(35))*DELINT(4)
      ACURL( 9) = EE(9)*DELINT(4) + AJJHO*DELINT(2)*EE(35)
      ACURL(11) = (EE(1) + 2.0*EE(3) + EE(15) + AJJHO*EE(29))*DELINT(6)
      ACURL(12) = (EE(3) + EE(15) + AJJHO*EE(29))*DELINT(5)
     1          + (EE(4) + EE (16))*DELINT(6)
      ACURL(13) = (EE(3) + EE(15) + EE(29))*AJHO*DELINT(4)
      ACURL(14) = (EE(3) + EE(15))*DELINT(6)*AJHO
      ACURL(15) = (EE(3) + EE(15) + EE(29))*AJHO*DELINT(5)
     1          - AJHO*EE(30)*DELINT(6)
      ACURL(16) = AJJHO*DELINT(4)*EE(35)
      ACURL(17) = (EE(4) + EE(16) + AJJHO*EE(35))*DELINT(6)
      ACURL(18) = (EE(2) + EE(9))*DELINT(6) + AJJHO*DELINT(5)*EE(35)
      ACURL(21) = (EE(15) + AJJHO*EE(29))*DELINT(3) + EE(22)*DELINT(6)
     1          + 2.0*EE(16)*DELINT(5)
      ACURL(22) = (EE(15) + EE(29))*AJHO*DELINT(2)
     1          + AJHO*DELINT(4)*EE(16)
      ACURL(23) = EE(15)*AJHO*DELINT(5) + AJHO*DELINT(6)*EE(16)
      ACURL(24) = (EE(15) + EE(29))*AJHO*DELINT(3)
     1          + (EE(16) - EE(30))*AJHO*DELINT(5)
      ACURL(25) = AJJHO*DELINT(2)*EE(35)
      ACURL(26) = EE(22)*DELINT(6) + (EE(21) + AJJHO*EE(35))*DELINT(5)
      ACURL(27) = EE(9)*DELINT(5) + EE(10)*DELINT(6)
     1          + AJJHO*DELINT(3)*EE(35)
      ACURL(31) = (EE(29) + AJJHO*EE(15))*DELINT(1)
      ACURL(32) = EE(15)*AJJHO*DELINT(4)
      ACURL(33) = (EE(29) + AJJHO*EE(15))*DELINT(2) - EE(30)*DELINT(4)
      ACURL(34) = AJHO*DELINT(1)*EE(35)
      ACURL(35) = AJHO*(EE(16) + EE(35))*DELINT(4)
      ACURL(36) = EE(9)*AJHO*DELINT(4) + AJHO*DELINT(2)*EE(35)
      ACURL(41) = AJJHO*DELINT(6)*EE(15)
      ACURL(42) = EE(15)*AJJHO*DELINT(5)
      ACURL(43) = 0.0
      ACURL(44) = AJHO*DELINT(6)*EE(16)
      ACURL(45) = EE(9)*AJHO*DELINT(6)
      ACURL(51) = (EE(29) + AJJHO*EE(15))*DELINT(3) + EE(36)*DELINT(6)
     1          - 2.0*DELINT(5)*EE(35)
      ACURL(52) = AJHO*(DELINT(2)*EE(30) - DELINT(4)*EE(36))
      ACURL(53) = -EE(36)*AJHO*DELINT(6) + AJHO*(EE(16)
     1          + EE(35))*DELINT(5)
      ACURL(54) = (EE(9) - EE(36))*AJHO*DELINT(5)
     1          + AJHO*DELINT(3)*EE(35)
      ACURL(61) = EE(36)*AJJHO*DELINT(1)
      ACURL(62) = EE(36)*AJJHO*DELINT(4)
      ACURL(63) = (EE(36))*AJJHO*DELINT(2)
      ACURL(71) = (EE(22) + AJJHO*EE(36))*DELINT(6)
      ACURL(72) = EE(36)*AJJHO*DELINT(5) + EE(20)*DELINT(6)
      ACURL(81) = EE(36)*AJJHO*DELINT(3) + EE(8)*DELINT(6)
C
      IF (LSYS78) GO TO 540
      ACURL(82) =-(EE(45)-AJHO*EE(51))*AJHO*DELINT(1)
      ACURL(83) = (EE(43)-AJHO*EE(45)-AJHO*EE(49)
     3          + AJJHO*EE(51))*DELINT(4)
      ACURL(84) = (EE(44)-AJHO*EE(50))*DELINT(4)-(EE(45)
     4          - AJHO*EE(51))*AJHO*DELINT(2)
      ACURL(85) =-(EE(39)+EE(45)-AJHO*EE(51))*AJHO*DELINT(4)
      ACURL(86) = (EE(37)+EE(43)-(EE(39)+EE(45)+EE(49)
     6          - AJHO*EE(51))*AJHO)*DELINT(6)
      ACURL(87) = (EE(38)+EE(44)-AJHO*EE(50))*DELINT(6)-(EE(39)+EE(45)
     7          - AJHO*EE(51))*AJHO*DELINT(5)
      ACURL(88) =-(EE(45)-AJHO*EE(51))*AJHO*DELINT(2)-EE(48)*AJHO*
     8            DELINT(4)
      ACURL(89) = (EE(43)-AJHO*EE(45)-AJHO*EE(49)+AJJHO*EE(51))*
     9            DELINT(5) +(EE(46)-EE(48)*AJHO)*DELINT(6)
      ACURL(90) = (EE(44)-AJHO*EE(48)-AJHO*EE(50))*DELINT(5)+EE(47)*
     O            DELINT(6) - (EE(45)-AJHO*EE(51))*AJHO*DELINT(3)
      ACURL(91) =-(EE(45)*AJHO-EE(51))*AJHO*DELINT(1)
      ACURL(92) = (AJHO*EE(43)-AJJHO*EE(45)-EE(49)+AJHO*EE(51))*
     2            DELINT(4)
      ACURL(93) = (EE(44)*AJHO-EE(50))*DELINT(4)-(EE(45)*AJHO-EE(51))*
     3            AJHO*DELINT(2)
      ACURL(94) =-EE(45)*AJJHO*DELINT(4)
      ACURL(95) = (EE(43)-AJHO*EE(45))*AJHO*DELINT(6)
      ACURL(96) = EE(44)*AJHO*DELINT(6)-EE(45)*AJJHO*DELINT(5)
      ACURL(97) =-(EE(45)*AJHO-EE(51))*AJHO*DELINT(2)-EE(54)*AJHO*
     7            DELINT(4)
      ACURL(98) = (EE(43)*AJHO-AJJHO*EE(45)-EE(49)+EE(51)*AJHO)*
     8            DELINT(5)+(EE(52)-AJHO*EE(54))*DELINT(6)
      ACURL(99) = (EE(44)*AJHO-EE(50)-EE(54)*AJHO)*DELINT(5)+EE(53)*
     9            DELINT(6)-(EE(45)*AJHO-EE(51))*AJHO*DELINT(3)
      ACURL(100) = EE(54)*AJJHO*DELINT(1)
      ACURL(101) =-(EE(52)-EE(54)*AJHO)*AJHO*DELINT(4)
      ACURL(102) =-(EE(53)*DELINT(4)-EE(54)*AJHO*DELINT(2))*AJHO
      ACURL(103) =-(EE(48)-EE(54)*AJHO)*AJHO*DELINT(4)
      ACURL(104) = (EE(46)-EE(48)*AJHO-EE(52)*AJHO+EE(54)*AJJHO)*
     4             DELINT(6)
      ACURL(105) = (EE(47)-EE(53)*AJHO)*DELINT(6)-(EE(48)-EE(54)*AJHO)*
     5             AJHO*DELINT(5)
      ACURL(106) = EE(54)*AJJHO*DELINT(2)-EE(42)*AJHO*DELINT(4)
      ACURL(107) = (EE(40)-EE(42)*AJHO)*DELINT(6)-(EE(52)-EE(54)*AJHO)*
     7             AJHO*DELINT(5)
      ACURL(108) = EE(41)*DELINT(6)+(-EE(42)-EE(53))*AJHO*DELINT(5)+
     8             EE(54)*AJJHO*DELINT(3)
      ACURL(109) = EE(63)*AJJHO*DELINT(1)
      ACURL(110) = (-EE(57)+EE(63)*AJHO)*AJHO*DELINT(4)
      ACURL(111) =-EE(60)*AJHO*DELINT(4)+EE(63)*AJJHO*DELINT(2)
      ACURL(112) = ACURL(110)
      ACURL(113) = (EE(55)-2.0*EE(57)*AJHO+EE(63)*AJJHO)*DELINT(6)
      ACURL(114) = (EE(56)-EE(60)*AJHO)*DELINT(6)+(-EE(57)+EE(63)*AJHO)*
     4             AJHO*DELINT(5)
      ACURL(115) = ACURL(111)
      ACURL(116) = ACURL(114)
      ACURL(117) = EE(59)*DELINT(6)-2.0*EE(60)*AJHO*DELINT(5)+EE(63)*
     7             AJJHO*DELINT(3)
  540 CONTINUE
C
C     EXPAND ACURL INTO (9X9)
C
      DO 610 IB = 2,9
      IC = 10*IB - 19
      I  = IC
      DO 605 J = IB,9
      IC = IC + 9
      I  = I + 1
  605 ACURL(IC) = ACURL(I)
  610 CONTINUE
      DGAMR = PI
      IF (AJHO .EQ. 0) DGAMR = 2.0*PI
      DO 630 I = 1,81
  630 ACURL(I) = DGAMR*ACURL(I)
C
      IF (LSYS78) GO TO 640
      DO 635 I = 82,117
  635 ACURL(I) = ACURL(I)*DGAMR
  640 CONTINUE
C
C     TRANSFORM THE ELEMENT STIFFNESS MATRIX FROM FIELD SYSTEM
C     TO GRID POINT DEGREES OF FREEDOM
C
      CALL GMMATS (GABABQ,9,9,1, ACURL,9,9,0,D )
      CALL GMMATS (     D,9,9,0,GABABQ,9,9,0,AK)
C
      IF (LSYS78) GO TO 650
      CALL GMMATS (GABABQ,9,9,1,ACURP1,9,3,0,D1   )
      CALL GMMATS (    D1,9,3,0,GABABP,3,3,0,AKUPH)
      CALL GMMATS (GABABP,3,3,1,ACURP2,3,3,0,D2   )
      CALL GMMATS (    D2,3,3,0,GABABP,3,3,0,AKPH2)
  650 CONTINUE
C
C
C     LOCATE THE TRANSFORMATION MATRICES FOR THE THREE GRID POINTS
C     **** COORDINATE SYS NOT POSSIBLE WITH RINGAX ********
C     **** THE FOLLOWING CODE COULD IMPLEMENT IT   ********
C     **   IF FOLLOWING CODE IS IMPLEMENTED MUST BE MODIFIED FOR
C          PIEZOELECTRIC
C.    DO 750 I = 1,81
C.750 AKI(I) = 0.00
C.    DO 800 I = 1,3
C.    CALL TRANSS (ICS(I),D(1))
C.    K = 30*(I-1) + 1
C.    DO 800 J = 1,3
C.    KK = K + 9*(J-1)
C.    JJ = 3*(J-1) + 1
C.    AKI(KK  ) = D(JJ  )
C.    AKI(KK+1) = D(JJ+1)
C.    AKI(KK+2) = D(JJ+2)
C.800 CONTINUE
C
C     TRANSFORM THE STIFFNESS MATRIX FROM BASIC TO LOCAL COORDINATES
C
C.    CALL GMMATS (AKI,9,9,1, AK,9,9,0,D )
C.    CALL GMMATS (  D,9,9,0,AKI,9,9,0,AK)
C
C     FORM WJ MATRIX
C
      DO 1000 I = 1,6
      DO 1000 J = 1,9
 1000 WJ (I,J) = 0.0
      RSUM = 0.0
      ZSUM = 0.0
      DO 1040 I = 1,3
      RSUM = RSUM + R(I)
 1040 ZSUM = ZSUM + Z(I)
      RSUM = RSUM/3.0
      ZSUM = ZSUM/3.0
      ZDR  = ZSUM/RSUM
      WJ (1,2) = 1.0
      WJ (2,9) = 1.0
      WJ (3,1) = 1.0/RSUM
      WJ (3,2) = 1.0
      WJ (3,3) = ZDR
      WJ (3,4) = AJHO/RSUM
      WJ (3,5) = AJHO
      WJ (3,6) = AJHO*ZDR
      WJ (4,3) = 1.0
      WJ (4,8) = 1.0
      WJ (5,1) =-AJHO/RSUM
      WJ (5,2) =-AJHO
      WJ (5,3) =-AJHO*ZDR
      WJ (5,4) =-1.0/RSUM
      WJ (5,6) =-ZDR
      WJ (6,6) = 1.0
      WJ (6,7) =-AJHO/RSUM
      WJ (6,8) =-AJHO
      WJ (6,9) =-AJHO*ZDR
C
      IF (LSYS78) GO TO 1060
C
C     FORM WJP MATRIX
C
      DO 1050 I = 1,3
      DO 1050 J = 1,3
 1050 WJP(I,J) = 0.0
C
      WJP(1,2) = 1.0
      WJP(2,3) = 1.0
      WJP(3,1) = -AJHO/RSUM
      WJP(3,2) = -AJHO
      WJP(3,3) = -AJHO*ZDR
 1060 CONTINUE
C
C     COMPUTE THE STRESS MATRIX
C
      CALL GMMATS (   WJ,9,6,1,GABABQ,9,9,0,D(1))
      CALL GMMATS (EE(1),6,6,0,  D(1),6,9,0,SEL )
      IF (LSYS78) GO TO 1070
      CALL GMMATS (   WJP,3,3,1,GABABP,3,3,0,D2(1))
      CALL GMMATS (EE(37),6,3,0, D2(1),3,3,0,SELP1)
      CALL GMMATS (EE(37),6,3,1,  D(1),6,9,0,SELP2)
      CALL GMMATS (EE(55),3,3,0, D2(1),3,3,0,SELP3)
 1070 CONTINUE
C
C    *** MORE CORD SYS REMOVAL.  FEL ABOVE IS WJ **********
C    **  IF FOLLOWING CODE IS IMPLEMENTED MUST BE MODIFIED FOR
C        PIEZOELECTRIC
C
C.    CALL GMMATS (WJ,6,9,0, AK,9,9,0, SEL)
C
C
C     TRANSFORM THE STRESS MATRIX FROM BASIC TO LOCAL COORDINATES
C
C     COMPUTE THE THE THERMAL STRESS
C
      TS(1) = EE( 1)*ALF(1) + EE( 2)*ALF(3) + EE( 3)*ALF(2)
      TS(2) = EE( 7)*ALF(1) + EE( 8)*ALF(3) + EE( 9)*ALF(2)
      TS(3) = EE(13)*ALF(1) + EE(14)*ALF(3) + EE(15)*ALF(2)
      TS(4) = EE(19)*ALF(1) + EE(20)*ALF(3) + EE(21)*ALF(2)
      TS(5) = 0.0
      TS(6) = 0.0
      DO 1100 IKI = 1,14
      PHI(IKI) = ECPT(7+IKI)
 1100 CONTINUE
      GO TO 940
C
C     SET FATAL ERROR FLAG AND ALLOWING ERROR MESSAGES TO ACCUMULATE
C
  910 I = 37
      GO TO 930
  920 I = 126
  930 CALL MESAGE (30,I,IDEL)
      NOGO = 1
  940 IF (.NOT.PZMAT) KSYS78 = KSAVE
      RETURN
      END

ext: f


text:
module EQS_KS
! Description:
!    Contains the subroutine for computing the RHS of the K-S equations in
!    double precision.
!
! Author:
!    Davide Amato
!    Space Dynamics Group - Technical University of Madrid
!    d.amato@upm.es
!
! ==============================================================================

use KINDS, only: ik,qk,dk
implicit none

contains

subroutine DKS_RHS(neq,s,u,udot)
! Description:
!    Computes the value of the right-hand side of the equations of motion of the
!    Kustaanheimo-Stiefel formulation.
!
! ==============================================================================
!                                  VARIABLES AND DECLARATIONS
! ==============================================================================

! MODULES
use AUXILIARIES,    only: inSoI
use TRANSFORM,      only: DKS2CART,DKSMAT
use THIRDBODIES,    only: DACC3B_ND
use CONSTANTS,      only: DU,TU,wEarth,smaEarth,muSun,muEarth

! VARIABLES
implicit none
! Arguments
integer(ik),intent(in)   ::  neq                  ! Number of equations.
real(dk),intent(in)      ::  s                    ! Value of fictitious time.
real(dk),intent(in)      ::  u(1:neq)             ! KS state vector.
real(dk),intent(out)     ::  udot(1:neq)          ! RHS of EoM's, ND.
! Local variables
! -- State variables
real(dk)    ::  x(1:4),xdot(1:4)    ! Radius and velocity in R^4, ND
real(dk)    ::  r                   ! Radius magnitude, ND
real(dk)    ::  t                   ! Physical time, ND
real(dk)    ::  L(1:4,1:4)          ! KS-matrix
! -- Perturbations
real(dk)    ::  rVpot               ! Perturbing potential term
real(dk)    ::  drVdu(1:4)          ! Derivatives of the perturbing potential term
real(dk)    ::  P(1:4)              ! Non-potential perturbing acceleration, ND
real(dk)    ::  r2(1:3)             ! Position vector of the third body [-]
real(dk)    ::  lon_Earth           ! Current Earth longitude [rad]
real(dk)    ::  Q(1:4)              ! Total acceleration, ND

! ==============================================================================

! STATE VECTOR DICTIONARY
! u(1:4)        u1,...,u4; KS-position, in R^4
! u(5:8)        u1',...,u4'; KS-velocity, in R^4
! u(9)          h; (-total energy) = (-Keplerian energy) + (-potential)
! u(10)         t; non-dimensional physical time

! ==============================================================================
! 01. COMPUTE CARTESIAN COORDINATES, ND
! ==============================================================================

call DKS2CART(u,x,xdot)
r = sqrt(dot_product(x,x))
t = u(10)

! ==============================================================================
! 02. COMPUTE PERTURBING POTENTIAL, ND
! ==============================================================================

! Initialize
rVpot = 0._dk; drVdu = 0._dk

! ==============================================================================
! 03. COMPUTE PERTURBING ACCELERATION, ND
! ==============================================================================

! Initialize
P = 0._dk
lon_Earth = wEarth*(t/TU)

if (inSoI) then
  ! Third body: Sun
  r2 = -(smaEarth/DU)*[cos(lon_Earth),sin(lon_Earth),0._dk]
  P(1:3) = DACC3B_ND(x(1:3),r2,muEarth,muSun)

else
  ! Third body: Earth
  r2 = (smaEarth/DU)*[cos(lon_Earth),sin(lon_Earth),0._dk]
  P(1:3) = DACC3B_ND(x(1:3),r2,muSun,muEarth)

end if

! Compute KS-matrix
L = DKSMAT(u)

! Total acceleration
Q = -drVdu/4._dk + .5_dk*r*matmul(transpose(L),P)

! ==============================================================================
! 04. COMPUTE RHS
! ==============================================================================

! Velocities
udot(1:4) = u(5:8)

! Accelerations
udot(5:8) = -.5_dk*u(9)*u(1:4) + Q

! Total energy
udot(9)   = -2._dk*dot_product(u(5:8),matmul(transpose(L),P))

! Time
udot(10)  = r

end subroutine DKS_RHS

end module EQS_KS

max_stars_repo_name: dkamato/NAPLES

max_issues_repo_name: dkamato/NAPLES

max_forks_repo_name: dkamato/NAPLES


text:
! ###################################################################
! Copyright (c) 2014-2019, Marc De Graef Research Group/Carnegie Mellon University
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
!     - Redistributions of source code must retain the above copyright notice, this list
!        of conditions and the following disclaimer.
!     - Redistributions in binary form must reproduce the above copyright notice, this
!        list of conditions and the following disclaimer in the documentation and/or
!        other materials provided with the distribution.
!     - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
!        of its contributors may be used to endorse or promote products derived from
!        this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################

!--------------------------------------------------------------------------
! EMsoft:files.f90
!--------------------------------------------------------------------------
!
! MODULE: files
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief everything that has to do with file-based input-output
!
!> @version
!
!> @date    1/ 5/99 MDG 1.0 original
!> @date    5/19/01 MDG 2.0 f90 version
!> @date   11/27/01 MDG 2.1 added kind support
!> @date   03/25/13 MDG 3.0 updated IO
!> @date   01/10/14 MDG 4.0 update after new cell type
!> @date   03/29/15 MDG 5.0 reformatted xtal files in HDF5 format; removed obsolete routines
!> @date   05/05/15 MDG 5.1 removed all getenv() calls; replaced with global path strings
!> @date   03/29/18 MDG 5.2 removed all stdout use
!--------------------------------------------------------------------------

module files

use local
use typedefs

contains

!--------------------------------------------------------------------------
!
! SUBROUTINE: DumpXtalInfo
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief Write a brief summary of the crystal structure on the screen
!
!> @date    1/ 5/99 MDG 1.0 original
!> @date    5/19/01 MDG 2.0 f90 version
!> @date   11/27/01 MDG 2.1 added kind support
!> @date   03/25/13 MDG 3.0 updated IO
!> @date   01/10/14 MDG 4.0 update after new cell type
!> @date   06/06/14 MDG 4.1 added cell pointer as argument, corrected Message routine
!> @date   07/31/18 MDG 4.2 added a few more output parameters
!--------------------------------------------------------------------------
recursive subroutine DumpXtalInfo(cell)
!DEC$ ATTRIBUTES DLLEXPORT :: DumpXtalInfo

use constants
use io
use symmetry

IMPLICIT NONE

type(unitcell), pointer                 :: cell

integer(kind=irg)                       :: i, j, oi_int(3)
real(kind=dbl)                          :: oi_real(5)

 call Message('', frm = "(A/)")
 call Message('Crystal Structure Information', frm = "('-->',A,'<--')")
 oi_real(1) = cell%a
 call WriteValue('  a [nm]             : ', oi_real, 1, "(F9.5)")
 oi_real(1) = cell%b
 call WriteValue('  b [nm]             : ', oi_real, 1, "(F9.5)")
 oi_real(1) = cell%c
 call WriteValue('  c [nm]             : ', oi_real, 1, "(F9.5)")
 oi_real(1) = cell%alpha
 call WriteValue('  alpha [deg]        : ', oi_real, 1, "(F9.5)")
 oi_real(1) = cell%beta
 call WriteValue('  beta  [deg]        : ', oi_real, 1, "(F9.5)")
 oi_real(1) = cell%gamma
 call WriteValue('  gamma [deg]        : ', oi_real, 1, "(F9.5)")
 oi_real(1) = cell%vol
 call WriteValue('  Volume [nm^3]      : ', oi_real, 1, "(F12.8)")
 oi_int(1) = cell%SYM_SGnum
 call WriteValue('  Space group #      : ', oi_int, 1, "(1x,I3)")
 call WriteValue('  Space group symbol : ', trim(SYM_SGname(cell%SYM_SGnum)) )
 call WriteValue('  Generator String   : ',  trim(SYM_GL(cell%SYM_SGnum)) )
 if ((cell%SYM_SGset.eq.2).AND.(cell%xtal_system.ne.5)) then
  call Message('   Using second origin setting', frm = "(A)")
 endif
 if ((cell%SYM_SGset.eq.2).AND.(cell%xtal_system.eq.5)) then
  call Message('   Using rhombohedral parameters', frm = "(A)")
 endif
  if (cell%SG%SYM_centrosym) then
    call Message('   Structure is centrosymmetric', frm = "(A)")
 else
   call Message('   Structure is non-centrosymmetric', frm = "(A)")
 end if

! space group and point group information
 oi_int(1) = cell%SG%SYM_GENnum
 call WriteValue('  # generators       : ', oi_int, 1, "(1x,I3)")
 oi_int(1) = cell%SG%SYM_MATnum
 call WriteValue('  # symmetry matrices: ', oi_int, 1, "(1x,I3)")
 oi_int(1) = cell%SG%SYM_NUMpt
 call WriteValue('  # point sym. matr. : ', oi_int, 1, "(1x,I3)")

! generate atom positions and dump output
 call Message('', frm = "(A/)")
 call CalcPositions(cell,'v')
 oi_int(1) = cell%ATOM_ntype
 call WriteValue('  Number of asymmetric atom positions ', oi_int, 1)
 do i=1,cell%ATOM_ntype
  oi_int(1:3) = (/i, cell%ATOM_type(i), cell%numat(i)/)
  call WriteValue('  General position / atomic number / multiplicity :', oi_int, 3,"(1x,I3,'/',I2,'/',I3,$)")
  call Message(' ('//ATOM_sym(cell%ATOM_type(i))//')', frm = "(A)")
  call Message('   Equivalent positions  (x y z  occ  DWF) ', frm = "(A)")
  do j=1,cell%numat(i)
    oi_real(1:5) = (/cell%apos(i, j,1:3),dble(cell%ATOM_pos(i,4:5))/)
    call WriteValue('         > ', oi_real, 5,"(2x,4(F9.5,','),F9.5)")
  end do
end do
call Message('', frm = "(A/)")

end subroutine DumpXtalInfo

!--------------------------------------------------------------------------
!
! SUBROUTINE: CopyTemplateFiles
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief  copy template files into local folder
!
!> @details In the resources folder, there is a text file called templatecodes.txt
!> in which each template file is given a unique ID number.  The present routine
!> receives the requested numbers and then looks into the file to figure out
!> which ones need to be copied.
!
!> @param nt number of files to copy
!> @param templatelist integer array identifying the template files to be copied
!
!> @date   06/26/13 MDG 1.0 first attempt
!> @date   06/08/14 MDG 2.0 added stdout argument
!> @date   10/06/14 MDG 2.1 corrected off-by-one error in templatelist
!> @date   05/05/15 MDG 2.2 removed getenv() call; replaced by global path string
!> @date   08/19/16 MDG 2.3 added handling of separate NameListTemplate folders for developers
!> @date   05/11/17 MDG 3.0 added support for json template files
!> @date   03/29/18 MDG 3.1 removed stdout argument
!--------------------------------------------------------------------------
recursive subroutine CopyTemplateFiles(nt,templatelist,json)
!DEC$ ATTRIBUTES DLLEXPORT :: CopyTemplateFiles

use io
use error

IMPLICIT NONE

integer(kind=irg),INTENT(IN)            :: nt
integer(kind=irg),INTENT(IN)            :: templatelist(*)
logical,INTENT(IN),OPTIONAL             :: json

integer(kind=irg),parameter     :: maxnumtemplates = 256
character(fnlen)                :: templates(maxnumtemplates)
character(fnlen)                :: input_name, output_name, tcf, tppath1, tppath2, tpl, tplextension
integer(kind=irg)               :: ios, i, j, ipos
character(255)                  :: line
logical                         :: fexist, develop

! first open and read the resources/templatecodes.txt file

tcf = trim(EMsoft_toNativePath(EMsoft_getTemplatecodefilename()))
inquire(file=trim(tcf),exist=fexist)
if (.not.fexist) then
  call FatalError('CopyTemplateFiles','template code file not found:'//tcf)
end if

open(UNIT=dataunit,FILE=trim(tcf), STATUS='old', FORM='formatted',ACCESS='sequential')

templates = ''
do
 read(dataunit,'(I3.3,A)',iostat=ios) j, line
! write (*,*) j,' ->',trim(line),'<-'
 if (ios.ne.0) then
  exit
 end if
 templates(j+1) = trim(line)
end do
CLOSE(UNIT=dataunit, STATUS='keep')

tplextension = '.template'
if (present(json)) then
  if (json.eqv..TRUE.) then
    tppath1 = trim(EMsoft_getTemplatepathname(json))
    tplextension = '.jtemplate'
  else
    tppath1 = trim(EMsoft_getTemplatepathname())
  end if
else
  tppath1 = trim(EMsoft_getTemplatepathname())
end if

! then, determine whether or not the user is working in develop mode by checking for the
! Develop keyword in the EMsoftconfig.json file... Regular users will only have a single
! NameListTemplates folder, but developers have two, so we need to make sure we check both
! locations.  The second location is the private folder...
develop = EMsoft_getEMdevelop()
tppath2 = ''
if (develop.eqv..TRUE.) then
  ipos = index(tppath1,'Public')
  do i=1,ipos-1
    tppath2(i:i) = tppath1(i:i)
  end do
  if (present(json)) then
    if (json.eqv..TRUE.) then
      tcf = 'Private/JSONTemplates/'
    else
      tcf = 'Private/NamelistTemplates/'
    end if
  else
    tcf = 'Private/NamelistTemplates/'
  endif
  do i=ipos,ipos+26
    j = i-ipos+1
    tppath2(i:i) = tcf(j:j)
  end do
end if

do i=1,nt
 tpl = trim(templates(templatelist(i)+1))
 input_name = trim(tppath1)//trim(tpl)//trim(tplextension)
 input_name = EMsoft_toNativePath(input_name)

 inquire(file=trim(input_name),exist=fexist)
 if (.not.fexist) then
  if (develop.eqv..TRUE.) then
   input_name = trim(tppath2)//trim(tpl)//trim(tplextension)
   input_name = EMsoft_toNativePath(input_name)
   inquire(file=trim(input_name),exist=fexist)
   if (.not.fexist) then
     call FatalError('CopyTemplateFiles','template file '//trim(templates(templatelist(i)+1))//trim(tplextension)// &
                    ' not found in either template folder')
   end if
  else
   call FatalError('CopyTemplateFiles','template file '//trim(templates(templatelist(i)+1))//trim(tplextension)//' not found')
  end if
 end if
 output_name = trim(templates(templatelist(i)+1))//trim(tplextension)
 output_name = EMsoft_toNativePath(output_name)
 open(UNIT=dataunit,FILE=trim(input_name), STATUS='old', FORM='formatted',ACCESS='sequential')
 open(UNIT=dataunit2,FILE=trim(output_name), STATUS='unknown', FORM='formatted',ACCESS='sequential')
 do
        read(dataunit,'(A)',iostat=ios) line
        if (ios.ne.0) then
          exit
        end if
        write(dataunit2,'(A)') trim(line)
  end do
 close(UNIT=dataunit, STATUS='keep')
 close(UNIT=dataunit2, STATUS='keep')
 call Message('  -> created template file '//trim(templates(templatelist(i)+1))//trim(tplextension), frm = "(A)")
end do

end subroutine CopyTemplateFiles

!--------------------------------------------------------------------------
!
! SUBROUTINE: Interpret_Program_Arguments
!
!> @author Marc De Graef, Carnegie Mellon University
!
!> @brief  interpret the command line arguments
!
!> @details In the resources folder, there is a text file called templatecodes.txt
!> in which each template file is given a unique ID number.  The present routine
!> receives the requested numbers and then looks into the file to figure out
!> which ones need to be copied.
!
!> @param nmldefault default nml file name on input, actual file name on output
!> @param numt number of files to (potentially) copy
!> @param templatelist integer array identifying the template files to be copied
!
!> @date   06/26/13 MDG 1.0 first attempt (this might belong elsewhere...)
!> @date   09/04/13 MDG 1.1 minor modifications to arguments
!> @date   06/08/14 MDG 2.0 added stdout argument
!> @date   05/11/17 MDG 3.0 added support for JSON-formatted template files
!> @date   03/29/18 MDG 3.1 removed stdout argument
!--------------------------------------------------------------------------
recursive subroutine Interpret_Program_Arguments(nmldefault,numt,templatelist,progname)
!DEC$ ATTRIBUTES DLLEXPORT :: Interpret_Program_Arguments

use io

IMPLICIT NONE

character(fnlen),INTENT(INOUT)          :: nmldefault
integer(kind=irg),INTENT(IN)            :: numt
integer(kind=irg),INTENT(IN)            :: templatelist(*)
character(fnlen),INTENT(IN)             :: progname

integer(kind=irg)                       :: numarg       !< number of command line arguments
integer(kind=irg)                       :: iargc        !< external function for command line
character(fnlen)                        :: arg          !< to be read from the command line
character(fnlen)                        :: nmlfile      !< nml file name
integer(kind=irg)                       :: i, io_int(1)
logical                                 :: haltprogram, json

json = .FALSE.

!numarg = iargc()
numarg = command_argument_count()
nmlfile = ''
nmlfile = trim(nmldefault)

if (numarg.gt.0) then
  io_int(1) = numarg
  call WriteValue('Number of command line arguments detected: ',io_int,1)
end if

haltprogram = .FALSE.
if (numarg.ge.1) haltprogram = .TRUE.

if (numarg.gt.0) then ! there is at least one argument
  do i=1,numarg
    call getarg(i,arg)
!    mess = 'Found the following argument: '//trim(arg); call Message("(/A/)")
! does the argument start with a '-' character?
    if (arg(1:1).eq.'-') then
        if (trim(arg).eq.'-h') then
         call Message(' Program should be called as follows: ', frm = "(/A)")
         call Message('        '//trim(progname)//' -h -t -j [nmlfile]', frm = "(A)")
         call Message(' where nmlfile is an optional file name for the namelist file;', frm = "(A/)")
         call Message(' If absent, the default name '''//trim(nmldefault)//''' will be used.', frm = "(A)")
         call Message(' To create templates of all possible input files, type '//trim(progname)//' -t', frm = "(A)")
         call Message(' To produce this message, type '//trim(progname)//' -h', frm = "(A)")
         call Message(' All program arguments can be combined in the same order;  ', frm = "(A)")
         call Message(' the argument without - will be interpreted as the input file name.', frm = "(A/)")
        end if
        if (trim(arg).eq.'-t') then
! with this option the program creates template namelist files in the current folder so that the
! user can edit them (file extension will be .template; should be changed by user to .nml)
                call Message('Creating program name list template files:', frm = "(/A)")
                call CopyTemplateFiles(numt,templatelist)
        end if
        if (trim(arg).eq.'-j') then
          json = .TRUE.
! with this option the program creates template JSON files in the current folder so that the
! user can edit them (file extension will be .jsontemplate; should be changed by user to .json)
!
! It should be noted that the template files contain comment lines starting with the "!" character;
! this is not standard JSON (which does not allow for comment lines).  The EMsoft JSON files will
! first be filtered to remove all the comment lines before being passed to the json parser routine.
                call Message('Creating program JSON template files:', frm = "(/A)")
                call CopyTemplateFiles(numt,templatelist,json=json)
        end if
    else
! no, the first character is not '-', so this argument must be the filename
! if it is present, but any of the other arguments were present as well, then
! we stop the program.
        nmlfile = arg
        if (numarg.eq.1) haltprogram = .FALSE.
    end if
  end do
end if

if (haltprogram) then
  call Message('To execute program, remove all flags except for nml/json input file name', frm = "(/A/)")
  stop
end if

nmldefault = nmlfile

end subroutine Interpret_Program_Arguments

end module files

max_stars_repo_name: ZachClayburn/EMsoft

max_issues_repo_name: ZachClayburn/EMsoft

max_forks_repo_name: ZachClayburn/EMsoft


text:
!--------------------------------------------------------------------------------------
!
!Copyright (C) 2014 TopCoder Inc., All Rights Reserved.
!
!--------------------------------------------------------------------------------------

!--------------------------------------------------------------------------------------
!
!   Copyright 2014 Ecole Centrale de Nantes, 1 rue de la No�, 44300 Nantes, France
!
!   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.
!
!   Contributors list:
!   - G. Delhommeau
!   - P. Guevel
!   - J.C. Daubisse
!   - J. Singh
!   - A. Babarit
!
!--------------------------------------------------------------------------------------

!   This module computes the Kochin function for problem already solved.
!
!
! Changes in version 1.2 (Implementation of Higher Order Panel Methods)
!       Added COMMON_TYPE module as dependency
!
!   @author yedtoss
!   @version 1.2

MODULE KOCHIN

    USE COMMON_TYPE
    !
    IMPLICIT NONE

CONTAINS
    SUBROUTINE COMPUTE_KOCHIN(kwave,Theta,HKochin, SolverVar)

        USE COM_VAR
        USE ELEMENTARY_FNS
        !
        !   Inputs / Outputs
        REAL :: kwave
        REAL :: THeta
        COMPLEX :: HKochin,HKochin1
        !   Locals
        COMPLEX,PARAMETER :: II=CMPLX(0.,1.)
        REAL :: PI
        INTEGER :: i,k,l
        REAL :: wbar
        COMPLEX,DIMENSION(NFA*2**NSYMY) :: ZS
        COMPLEX,DIMENSION(NP) :: CEP,CEM,ZJ
        INTEGER,DIMENSION(5) :: KK
        COMPLEX :: ZMS,ZNS,ZCS,ZC,ZHS,ZHD

        TYPE(TempVar), TARGET :: SolverVar
        REAL, POINTER :: T
        COMPLEX, DIMENSION(:), POINTER :: ZPB,ZPS
        COMPLEX, DIMENSION(:), POINTER :: ZIGB,ZIGS
        COMPLEX, DIMENSION(:, :), POINTER :: ZIJ
        REAL, POINTER :: FSP,FSM,VSXP,VSYP,VSZP,VSXM,VSYM,VSZM
        REAL, POINTER :: SP1,SM1,SP2,SM2
        REAL, POINTER :: VSXP1,VSXP2,VSYP1,VSYP2,VSZP1,VSZP2
        REAL, POINTER :: VSXM1,VSXM2,VSYM1,VSYM2,VSZM1,VSZM2
        INTEGER, POINTER:: NQ
        REAL, POINTER:: CQ(:),QQ(:),AMBDA(:),AR(:)
        T => SolverVar%T
        ZPB => SolverVar%ZPB
        ZPS => SolverVar%ZPS
        ZIGB => SolverVar%ZIGB
        ZIGS => SolverVar%ZIGS
        ZIJ => SolverVar%ZIJ
        FSP => SolverVar%FSP
        FSM => SolverVar%FSM
        VSXP => SolverVar%VSXP
        VSYP => SolverVar%VSYP
        VSZP => SolverVar%VSZP
        VSXM => SolverVar%VSXM
        VSYM => SolverVar%VSYM
        VSZM => SolverVar%VSZM
        SP1 => SolverVar%SP1
        SM1 => SolverVar%SM1
        SP2 => SolverVar%SP2
        SM2 => SolverVar%SM2
        VSXP1 => SolverVar%VSXP1
        VSXP2 => SolverVar%VSXP2
        VSYP1 => SolverVar%VSYP1
        VSYP2 => SolverVar%VSYP2
        VSZP1 => SolverVar%VSZP1
        VSZP2 => SolverVar%VSZP2
        VSXM1 => SolverVar%VSXM1
        VSXM2 => SolverVar%VSXM2
        VSYM1 => SolverVar%VSYM1
        VSYM2 => SolverVar%VSYM2
        VSZM1 => SolverVar%VSZM1
        VSZM2 => SolverVar%VSZM2
        NQ => SolverVar%NQ
        CQ => SolverVar%CQ
        QQ => SolverVar%QQ
        AMBDA => SolverVar%AMBDA
        AR => SolverVar%AR
        !
        PI=4.*ATAN(1.)
        !   Compute Kochin coefficients (integration using Gauss - 1 point)
        HKochin1=0.
        !    DO i=1,NFA
        !        wbar=(XG(i)-XEFF)*COS(Theta)+(YG(i)-YEFF)*SIN(Theta)
        !        HKochin1=HKochin1+ZIGB(i)*CIH(kwave,ZG(i),Depth)*CEXP(-II*kwave*wbar)*AIRE(i)
        !        IF (NSYMY.EQ.1) THEN
        !            wbar=(XG(i)-XEFF)*COS(Theta)+(-YG(i)-YEFF)*SIN(Theta)
        !            HKochin1=HKochin1+ZIGS(i)*CIH(kwave,ZG(i),Depth)*CEXP(-II*kwave*wbar)*AIRE(i)
        !        END IF
        !    END DO
        !   Compute Kochin coefficients (analytical integration)
        DO k=1,NP
            wbar=(X(k)-XEFF)*COS(Theta)+(Y(k)-YEFF)*SIN(Theta)
            ZJ(k)=kwave*(Z(k)-II*wbar)
            IF (REAL(ZJ(k)).LT.-18.) THEN
                CEP(k)=(0.,0.)
            ELSE
                CEP(k)=CEXP(ZJ(k))
            END IF
            IF ((REAL(-ZJ(k)-2.*kwave*Depth).LT.-18.).OR.(-kwave*Depth.GE.0.)) THEN
                CEM(k)=(0.,0.)
            ELSE
                CEM(k)=CEXP(-ZJ(k)-2.*kwave*Depth)
            END IF
        END DO
        DO i=1,NFA
            KK(1)=M1(i)
            KK(2)=M2(i)
            KK(3)=M3(i)
            KK(4)=M4(i)
            KK(5)=KK(1)
            ZMS=YN(i)-II*ZN(i)*SIN(Theta)
            ZNS=XN(i)-II*ZN(i)*COS(Theta)
            ZS(i)=(0.,0.)
            DO L=1,4
                ZCS=ZMS*(X(KK(l+1))-X(KK(l)))-ZNS*(Y(KK(l+1))-Y(KK(l)))
                ZC=ZJ(KK(L+1))-ZJ(KK(L))
                IF(ABS(AIMAG(ZC)).LT.1.E-04.AND.ABS(REAL(ZC)).LT.1.E-04)THEN
                    ZHS=0.5*(CEP(KK(L+1))+CEP(KK(L)))
                    ZHD=-0.5*(CEM(KK(L+1))+CEM(KK(L)))
                ELSE
                    ZHS=(CEP(KK(L+1))-CEP(KK(L)))/ZC
                    ZHD=(CEM(KK(L+1))-CEM(KK(L)))/ZC
                ENDIF
                ZS(i)=ZS(i)+ZCS*ZHS+CONJG(ZCS*ZHD)
            END DO
            IF ((kwave*Depth.GT.18.).OR.(kwave*Depth.LE.0.)) THEN
                ZS(i)=ZS(i)/kwave
            ELSE
                ZS(i)=0.5*ZS(i)*EXP(kwave*Depth)/SINH(kwave*Depth)/kwave
            END IF
        END DO
        IF (NSYMY.EQ.1) THEN
            DO k=1,NP
                wbar=(X(k)-XEFF)*COS(Theta)+(-Y(k)-YEFF)*SIN(Theta)
                ZJ(k)=kwave*(Z(k)-II*wbar)
                IF (REAL(ZJ(k)).LT.-18.) THEN
                    CEP(k)=(0.,0.)
                ELSE
                    CEP(k)=CEXP(ZJ(k))
                END IF
                IF ((REAL(-ZJ(k)-2.*kwave*Depth).LT.-18.).OR.(-kwave*Depth.GE.0.)) THEN
                    CEM(k)=(0.,0.)
                ELSE
                    CEM(k)=CEXP(-ZJ(k)-2.*kwave*Depth)
                END IF
            END DO
            DO i=1,NFA
                KK(1)=M1(i)
                KK(2)=M2(i)
                KK(3)=M3(i)
                KK(4)=M4(i)
                KK(5)=KK(1)
                ZMS=-YN(i)-II*ZN(i)*SIN(Theta)
                ZNS=XN(i)-II*ZN(i)*COS(Theta)
                !        Initialisation de l integrale sur la facette J
                ZS(i+NFA)=(0.,0.)
                DO L=1,4
                    ZCS=ZMS*(X(KK(l+1))-X(KK(l)))+ZNS*(Y(KK(l+1))-Y(KK(l)))
                    ZC=ZJ(KK(L+1))-ZJ(KK(L))
                    IF(ABS(AIMAG(ZC)).LT.1.E-04.AND.ABS(REAL(ZC)).LT.1.E-04)THEN
                        ZHS=0.5*(CEP(KK(L+1))+CEP(KK(L)))
                        ZHD=-0.5*(CEM(KK(L+1))+CEM(KK(L)))
                    ELSE
                        ZHS=(CEP(KK(L+1))-CEP(KK(L)))/ZC
                        ZHD=(CEM(KK(L+1))-CEM(KK(L)))/ZC
                    ENDIF
                    ZS(i+NFA)=ZS(i+NFA)+ZCS*ZHS+CONJG(ZCS*ZHD)
                END DO
                IF ((kwave*Depth.GT.18.).OR.(kwave*Depth.LE.0.)) THEN
                    ZS(i+NFA)=-ZS(i+NFA)/kwave
                ELSE
                    ZS(i+NFA)=-0.5*ZS(i+NFA)*EXP(kwave*Depth)/SINH(kwave*Depth)/kwave
                END IF
            END DO
        END IF
        HKochin=0.
        DO i=1,NFA
            HKochin=HKochin+ZIGB(i)*ZS(i)
            IF (NSYMY.EQ.1) THEN
                HKochin=HKochin+ZIGS(i)*ZS(i+NFA)
            END IF
        END DO
        HKochin=0.25/PI*HKochin
    !
    END SUBROUTINE COMPUTE_KOCHIN

END MODULE KOCHIN

max_stars_repo_name: NREL/OpenWARP

max_issues_repo_name: rhydar/Test

max_forks_repo_name: rhydar/Test


text:
Im a plant bio PhD student at UC Davis. Just moved to Davis September 2012

welcome to the wiki! ~Daubert Any idea as to what your dissertation will be about?

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text: This library is just north of Stevenson Bridge on the East side of the street.


text: #redirect Segundo Services Center


text:
program zero_size
implicit none
associate (r0 => [real ::], d0 => [real(kind(1.0d0)) ::], &
           i0 => [integer ::], c20 => [character (len=2) ::], &
           c50 => [character (len=5) ::])
print*,size(r0),kind(r0)  ! 0 4
print*,size(d0),kind(d0)  ! 0 8
print*,size(i0),kind(i0)  ! 0 4
print*,size(c20),len(c20) ! 0 2
print*,size(c50),len(c50) ! 0 5
end associate
end program zero_size

max_stars_repo_name: awvwgk/FortranTip

max_issues_repo_name: awvwgk/FortranTip

max_forks_repo_name: awvwgk/FortranTip


text:
c
c   Subroutine to find the JH2 yield stress, which is dependent on
c   pressure, strain rate and damage accumulated. All values are
c   normalized with respect to the HEL (stress wrt SHEL, pressure
c   wrt PHEL).
c
c   inputs
c   ======
c   epd = equivalent strain rate         (history var)
c   P = current hydrostatic pressure     (history var)
c   D = current damage                   (history var)
c   A = intact strength coefficient      (JH2 param.)
c   B = fractured strength coefficient   (JH2 param.)
c   c = strain rate coefficient          (JH2 param.)
c   n = intact strength exponent         (JH2 param.)
c   m = fractured strength exponent      (JH2 param.)
c   EPSI = quasi-static threshold s.rate (JH2 param.)
c   T = hydrostatic pressure to failure  (JH2 param.)
c   SFmax = maximum failure strength     (JH2 param.)
c   PHEL = Pressure comp. of HEL         (JH2 param.)
c   SHEL = Deviatoric comp. of HEL       (JH2 param.)
c
c   outputs
c   =======
c   syield = yield stress
c
        subroutine JH2_yield (
     &             epd,P,D,A,B,c,n,m,EPSI,T,SFmax,HEL,PHEL,syield)
c
         implicit none
         double precision, intent (in) :: epd,P,D,A,B,c,n,m,T,
     &                                    SFmax,HEL,PHEL,EPSI
         double precision, intent (inout) :: syield
         double precision :: si, sf, Ps, Ts, es, SHEL
c
c        Find normalized pressure and tensile strength constant (wrt PHEL)
c        =================================================================
         Ps = P / PHEL
         Ts = T / PHEL
         SHEL = 1.5_8 * (HEL - PHEL)
c
c        Find normalized strain rate (wrt EPSI) and find strain rate term
c        ================================================================
         es = 1.0_8+c*log(epd/EPSI)
c
c        Find normalized strength at D=0 (si) , D=1 (sf)
c        ===============================================
         si = A * (Ps+Ts) ** n * es
         sf = B * Ps ** m * es
         print *,"strain scale"
         print "(2e12.4)",es
c
c        Check if normalized fractured strength exceeds SFmax, if so sf=SFmax
c        ====================================================================
         if (sf.GT.SFmax) then
            sf = SFmax
         end if
c
c        Calculate flow stress based on si, sf and damage
c        ================================================
         syield = SHEL * (si - D*si + D*sf)
c
         return
c
        end subroutine JH2_yield
c
c  Subroutine to update damage in the JH2 model
c
c  inputs
c  ======
c  eps = effective plastic strain increment
c  P = current hydrostatic pressure
c  T = hydrostatic tensile stress to failure
c  PHEL = pressure component of HEL
c  D1 = damage coefficient
c  D2 = damage exponent
c
c  outputs
c  =======
c  D = damage parameter (to be updated in subroutine)
c
c
        subroutine JH2_damage (eps,P,T,PHEL,D1,D2,D)
c
         double precision, intent (in) :: eps,P,T,PHEL,D1,D2
         double precision, intent (inout) :: D
         double precision :: epf
c
         Ps = P / PHEL
         Ts = T / PHEL
c
         epf = D1 * (Ps + Ts) ** D2
         D = D + eps / epf
c
         return
c
        end subroutine JH2_damage
c
c  Subroutine to calculate equivalent stress as well as normals
c  required for the plastic loop
c
c  inputs
c  ======
c  A = input tensor
c  cflag = changes outputs of the function (note only cflag=2 is implemented to date)
c
c  outputs
c  =======
c  seq = von Mises stress as output
c  s1d = 1st derivatives of Seq wrt S (is the nomal tensor in case of von Mises)
c  s2d = 2nd derivatives of Seq	wrt S (apparently used for implicit formulations, not implemented here)
c
c
        subroutine vm_box (A, cflag, seq, s1d, s2d)
c
         implicit none
         double precision, intent (in) :: A(9)
         double precision, intent (out) :: seq, s1d(9), s2d(9,9)
         double precision :: S(9)
         integer, intent(in) :: cflag
c
c
         call VMS (A,seq)
c
         s1d = 0.0_8
         s2d = 0.0_8

         if (cflag.eq.2) then
            call A2S (A,S)
c            call t2print (S)
            s1d = (1.5_8/seq)*S
         else if (cflag.eq.3) then
         else
         end if
c
         return
c
        end subroutine vm_box
c
c
c
c  inputs
c  ======
c  P = hydrostatic pressure
c  D = damage
c  A = intact strength coefficent
c  B = damaged strength coefficient
c  M = damaged strength exponent
c  N = intact strength exponent
c  HEL = hugoniot elastic limit
c  PHEL = hydrostatic pressure at HEL
c  Lel = elastic tensor
c  N = normal vector (from vm_box)
c  D1 = damage accumulation coefficient
c  D2 = damage accumulation exponent
c
c  outputs
c  =======
c  h = hardening rate (ds_bar/dep)
c
        subroutine JH2_hardening
     &  (epd,P,D,A,B,c,M,N,T,HEL,PHEL,SFmax,EPSI,Lel,NOR,D1,D2,h)

c
          implicit none
          double precision, intent (in) :: epd,P,D,A,B,c,M,N,T
          double precision, intent (in) :: HEL,PHEL,SFmax,EPSI
          double precision, intent (in) :: Lel(9,9), NOR(9), D1, D2
          double precision, intent (inout) :: h
          double precision :: es, si, sf, ef, Ps, Ts, SHEL
          double precision :: alp, bet, kap, KIJ(9), LN(9), KLN
c
c
c        Find normalized pressure and tensile strength constant (wrt PHEL)
c        =================================================================
         Ps = P / PHEL
         Ts = T / PHEL
         SHEL = 1.5 * (HEL - PHEL)
c
c        Find normalized strain rate (wrt EPSI) and find strain rate term
c        ================================================================
         es = 1.0_8+c*log(epd/EPSI)
c
c        Find normalized strength at D=0 (si) , D=1 (sf)
c        ===============================================
         si = A * (Ps+Ts) ** n * es
         sf = B * Ps ** m * es
c
c        Check if normalized fractured strength exceeds SFmax, if so sf=SFmax
c        ====================================================================
         if (sf.GT.SFmax) then
            sf = SFmax
         end if
c
c        Find failure strain at given pressure
c        =====================================
         ef = D1 ** (Ps+Ts) ** D2
         if (ef.lt.1e-20) then
            ef = 1e-20
         end if
c
c        Find intermediates to find total derivative
c        ===========================================
         call KD (KIJ)
         call tc_4d2 (Lel,NOR,LN)
         call tc_2d2 (KIJ,LN,KLN)
         kap = Ps * KLN / PHEL / 3.0_8
         alp = A * N * (Ps-Ts)**(n-1.0_8)
         bet = B * m * (Ps)**(m-1.0_8)
c
         h = SHEL * ((1.0_8-D)*alp*kap+D*bet*kap+(sf-si)/ef)
c
         return
c
        end subroutine JH2_hardening
c
c
        subroutine umat_JH2 (cm, sig, deps, hisv, tt, dt)
c
          implicit none
          double precision, intent (in) :: cm (20), deps (9), tt, dt
          double precision, intent (inout) :: sig (9), hisv (12)
          double precision :: L4 (9,9), strial (9), dsig(9), E, v
          double precision :: seq, sbar, h, sflow, NOR(9), ph (9,9)
          double precision :: LN (9), NLN, pl_corr(9), SHEL, syield
          double precision :: MID,RO,G,A,B,C,M,N,EPSI,T,SFmax,HEL,PHEL
          double precision :: BETA, D1, D2, K1, K2, K3, FS
          double precision :: P,D, dlam, ddlam, dep, ef, Ps, Ts, sbar_t
          double precision :: tol, epd
          integer cnt
c
c         cm(1) = MID   , material id
c         cm(2) = RO    , density
c         cm(3) = G     , shear modulus
c         cm(4) = A     , intact material strength coefficient
c         cm(5) = B     , fractured material strength exponent
c         cm(6) = C     , strain rate coefficent
c         cm(7) = M     , fractured material strength exponent
c         cm(8) = N     , intact material strength coefficient
c         cm(9) = EPSI  , quasi-static strain rate threshold
c         cm(10) = T    , hydrostatic tensile strength
c         cm(11) = SFmax, maximum normalized fractured strength
c         cm(12) = HEL  , hugoniot elastic limit
c         cm(13) = PHEL , pressure component of HEL
c         cm(14) = BETA , fraction of elastic energy converted to pressure
c         cm(15) = D1   , damage accumulation coefficient
c         cm(16) = D2   , damage accumulation exponent
c         cm(17) = K1   , EOS term 1 (if K2,K3=0, K1 is the bulk mod.)
c         cm(18) = K2   , EOS term 2
c         cm(19) = K3   , EOS term 3
c         cm(20) = FS   , failure criteria
c
         MID = cm(1)
         RO = cm(2)
         G = cm(3)
         A = cm(4)
         B = cm(5)
         C = cm(6)
         M = cm(7)
         N = cm(8)
         EPSI = cm(9)
         T = cm(10)
         SFmax = cm(11)
         HEL = cm(12)
         PHEL = cm(13)
         BETA = cm(14)
         D1 = cm(15)
         D2 = cm(16)
         K1 = cm(17)
         K2 = cm(18)
         K3 = cm(19)
         FS = cm(20)
         D = hisv(12)
         call I1 (sig, P)
         P = P / 3.0_8
         syield = 0.0_8
         tol = G * 1.0e-6_8
         epd = 0.0_8
c
         if (ABS(epd).lt.1e-20) then
            epd = 1e-20
         end if
c
c ================================
c    1. Calculate Trial Stress
c ================================
c
c         Get elastic constants from G,K
c         ==============================
          E = 9.0_8*G*K1/(3.0_8*K1+G)
          v = (3.0_8*K1-2.0_8*G)/(2.0_8*(3.0_8*K1+G))
c          print "(2e12.4)",E,v
          SHEL = 1.5_8*(HEL-PHEL)
c
c         Get elastic fourth order tensor L4
c         ==================================
          call L4_EL (E,v,L4)
c          print *,
c          call t4print (L4)
c          print *,
c
c         Calculate stress increment (dsig = L:D)
c         =======================================
          call tc_4d2 (L4,deps,dsig)
c          print *,
c          call t2print (dsig)
c          print *,
c
c         Update strial with dsig (sn+1,trial = sn + dsig)
c         ================================================
          strial = sig + dsig
c
c ================================
c      2. Find Seq and Syield
c ================================
c
c         Find equivalent stress and normal using vm_box
c         ==============================================
          call vm_box (sig,2,seq,NOR,ph)
c
c         Find yield stress using JH2_yield
c         =================================
          call JH2_yield
     &    (epd,P,D,A,B,c,M,N,EPSI,T,SFmax,HEL,PHEL,syield)
          print *,"Yield Strength"
          print "(1e12.4)", syield
          print *,
c
c         Check for yielding
c         ==================
          sflow = seq - syield
c
          if (sflow.lt.0.0_8) then
             sig = strial
          else
c         =====================
c         Initiate PLASTIC LOOP
c         =====================
          call JH2_hardening
     &    (epd,P,D,A,B,c,M,N,T,HEL,PHEL,SFmax,EPSI,L4,NOR,D1,D2,h)
c
c  prime the loop, set all counters to 0
c  set sflow = 100 to ensure this runs once
c
            dlam = 0.0_8
            dep = 0.0_8
            cnt = 1
            sflow = 100.0_8
c
c           plastic loop time!
c
            do while (DABS(sflow).gt.tol.AND.cnt.lt.5)
c
c            +  use normal and elastic tensor to calculate
c            +  plastic correction
c
               call tc_4d2(L4,NOR,LN)
               pl_corr = dlam * LN
               sig = strial - pl_corr
c
c            +  recalculate seq and sflow to determine whether
c            +  material is still yielding
c
               call VMS (sig,seq)
               sbar_t = sbar + h*dep
               sflow = seq - sbar_t
c
c            +  update plastic multiplier increment, plastic multiplier
c            +  and plastic strain (the same in assoc. flow) to reflect
c            +  plastic correction
c
               call tc_2d2 (NOR,LN,NLN)
               ddlam = sflow / (NLN+h)
               dlam = dlam + ddlam
               dep = dlam
               cnt = cnt + 1
c
            end do
            hisv(1) = hisv(1) + dep
            hisv(11) = seq
c
c           +  update damage variable based on plastic strain increment
c
            Ps = P / PHEL
            Ts = T / PHEL
            ef = D1 ** (Ps+Ts) ** D2
            D = D + dep / ef
            hisv(12) = D
          end if
c         return
c
        end subroutine umat_JH2
c

max_stars_repo_name: kengwit/UMAT

max_issues_repo_name: kengwit/UMAT

max_forks_repo_name: kengwit/UMAT


text:
c     Subroutine to resolve qc info
c     AJ_Kettle, 14Dec2020
c     01Feb2021: slave designations for variables passed into subroutine

      SUBROUTINE resolve_qc_single_param(
     +  s_qckey_timescale_spec,l_qckey,
     +  s_qckey_sourceflag,s_qckey_c3sflag,s_qckey_timescale,
     +  l_timestamp_rgh,l_timestamp,
     +  s_prcp_qflag,

     +  s_prcp_cdmqc,s_prcp_cdmqcmethod,
     +  i_cnt_qc)

      IMPLICIT NONE
c************************************************************************
c     Variables passed into subroutine

      CHARACTER(LEN=300)  :: s_qckey_timescale_spec
      INTEGER             :: l_qckey
      CHARACTER(LEN=*)    :: s_qckey_sourceflag(100)             !2
      CHARACTER(LEN=*)    :: s_qckey_c3sflag(100)                !2
      CHARACTER(LEN=300)  :: s_qckey_timescale(100)

      INTEGER             :: l_timestamp_rgh
      INTEGER             :: l_timestamp

      CHARACTER(LEN=*)    :: s_prcp_qflag(l_timestamp_rgh)       !1

      CHARACTER(LEN=*)    :: s_prcp_cdmqc(l_timestamp_rgh)       !1
      CHARACTER(LEN=*)    :: s_prcp_cdmqcmethod(l_timestamp_rgh) !2

      INTEGER             :: i_cnt_qc
c******
c     Variables used in subroutine
      INTEGER             :: i,j,k,ii,jj,kk

      INTEGER             :: ilen
c************************************************************************
      ii=0

      DO i=1,l_timestamp

c      Initialize
       s_prcp_cdmqc(i)='0'
       s_prcp_cdmqcmethod(i)=''

       ilen=LEN_TRIM(s_prcp_qflag(i))

c      act if something in source qc field
       IF (ilen.GT.0) THEN
        s_prcp_cdmqc(i)='1'

c        print*,'problem identified',i,ilen,s_prcp_qflag(i)
c        print*,'s_qckey_timescale_spec=',TRIM(s_qckey_timescale_spec)

        DO j=1,l_qckey
c        Match time scales
         IF (TRIM(s_qckey_timescale(j)).EQ.TRIM(s_qckey_timescale_spec))
     +     THEN
c        comparing 1 & 2 element string
         IF (TRIM(s_prcp_qflag(i)).EQ.TRIM(s_qckey_sourceflag(j))) THEN
c          print*,'match found'
          ii=ii+1
          s_prcp_cdmqcmethod(i)=TRIM(ADJUSTL(s_qckey_c3sflag(j)))   !both length 2
          GOTO 10
         ENDIF
         ENDIF
        ENDDO

        print*,'problem: qc flag info not found',
     +    '='//TRIM(s_prcp_qflag(i))//'='

        DO j=1,l_qckey
         print*,'j=',j,'='//TRIM(s_qckey_sourceflag(j))//'='
        ENDDO

        STOP 'qc_converter'

 10     CONTINUE
       ENDIF

      ENDDO

      i_cnt_qc=ii

c      print*,'ii=',ii,l_timestamp

      RETURN
      END

max_stars_repo_name: glamod/glamod-nuim

max_issues_repo_name: glamod/glamod-nuim

max_forks_repo_name: glamod/glamod-nuim


! 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
!

! BEGIN DART PREPROCESS TYPE DEFINITIONS
!RAW_STATE_VARIABLE,           QTY_STATE_VARIABLE,        COMMON_CODE
!VELOCITY,                     QTY_VELOCITY,              COMMON_CODE
!TRACER_CONCENTRATION,         QTY_TRACER_CONCENTRATION,  COMMON_CODE
!TRACER_SOURCE,                QTY_TRACER_SOURCE,         COMMON_CODE
!MEAN_SOURCE,                  QTY_MEAN_SOURCE,           COMMON_CODE
!SOURCE_PHASE,                 QTY_SOURCE_PHASE,          COMMON_CODE
! END DART PREPROCESS TYPE DEFINITIONS


text:
! { dg-do run }
! Checks the fixes for PR34681 and PR34704, in which various mixtures
! of default initializer and allocatable array were not being handled
! correctly for derived types with allocatable components.
!
! Contributed by Paolo Giannozzi <p.giannozzi@fisica.uniud.it>
!
program boh
  integer :: c1, c2, c3, c4, c5
  !
  call mah (0, c1) ! These calls deal with PR34681
  call mah (1, c2)
  call mah (2, c3)
  !
  if (c1 /= c2) call abort
  if (c1 /= c3) call abort
  !
  call mah0 (c4) ! These calls deal with PR34704
  call mah1 (c5)
  !
  if (c4 /= c5) call abort
  !
end program boh
!
subroutine mah (i, c)
  !
  integer, intent(in) :: i
  integer, intent(OUT) :: c
  !
  type mix_type
     real(8), allocatable :: a(:)
     complex(8), allocatable :: b(:)
  end type mix_type
  type(mix_type), allocatable, save :: t(:)
  integer :: j, n=1024
  !
  if (i==0) then
     allocate (t(1))
     allocate (t(1)%a(n))
     allocate (t(1)%b(n))
     do j=1,n
        t(1)%a(j) = j
        t(1)%b(j) = n-j
     end do
  end if
  c = sum( t(1)%a(:) ) + sum( t(1)%b(:) )
  if ( i==2) then
     deallocate (t(1)%b)
     deallocate (t(1)%a)
     deallocate (t)
  end if
end subroutine mah

subroutine mah0 (c)
  !
  integer, intent(OUT) :: c
  type mix_type
     real(8), allocatable :: a(:)
     integer :: n=1023
  end type mix_type
  type(mix_type) :: t
  !
  allocate(t%a(1))
  t%a=3.1415926
  c = t%n
  deallocate(t%a)
  !
end subroutine mah0
!
subroutine mah1 (c)
  !
  integer, intent(OUT) :: c
  type mix_type
     real(8), allocatable :: a(:)
     integer :: n=1023
  end type mix_type
  type(mix_type), save :: t
  !
  allocate(t%a(1))
  t%a=3.1415926
  c = t%n
  deallocate(t%a)
  !
end subroutine mah1

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
The Davis Spinners Guild is a friendly, growing group of fiber and Crafts craft enthusiasts. We welcome members of all ages and skill levels. The guild typically meets at a members house on Sunday afternoons: we bring our projects, equipment, mystery fiber for identification, and sometimes baked goods or excess produce from our gardens.

Our members can often be found spinning at the Farmers Market, at the Square Tomatoes Craft Fair, on the The CoOp CoOp patio, or hanging out at the http://www.ravelry.com/groups/dfcknittingandcraftcircle Davis Food Coop Knitting and Craft Circle. If you see us, come say hello!

 Local Resources

 Yolo Wool Mill in Woodland does commercial fiber preparation and spinning. They also sell naturalcolored and blended dyed prepared fiber.
 http://www.meridianjacobs.com Meridian Jacobs is a farm on the Dixon/Vacaville line. They have a spinningoriented shop open on most Thursdays and by appointment. Robin also hosts a Spinners Night Out one Friday per month.
 http://www.fiberconfections.com Fiber Confections has both natural and dyed roving. Colleen is often at the Farmers Market.
 http://winddancerranch.us Wind Dancer Ranch in the Capay Valley raises NavajoChurro sheep and sells undyed shearings, batts, and roving.
 http://www.vreseis.com/ Sally Fox in Brooks grows organic naturallycolored cotton and raises Merino sheep. She sells sliver/roving from both.
 http://www.yarnyarnyarn.com/ Rumplestiltskin in Sacramento has dyed New Zealand roving, spinning wheels, and fiber prep tools.
 http://babettasyarnandgifts.com Babettas Yarns and Gifts in Fair Oaks carries a variety of plant and animal fibers.


text:
Caffino serves a variety of drive though hot/cold drinks, both caffeinated and decaffeinated as well as pastries. Caffinos prices are lower than either Starbucks or Silo Brenans, and the service is much faster, even when theres a line of cars in the morning. If youre a regular, they usually remember your drink, how you like it, and which side of the cup you sip from
 now thats what I call service!

Replaced by Dutch Bros Coffee.

20060220 16:36:10 nbsp The employees at Caffino are not only fast, efficient and courteous, they are the nicest group of people anywhere. Its the best of both worlds, Corporation services with community kindness...They Rock! KJS Users/KarenaSkinner

20060818 07:45:22 nbsp Caffino is one of my favorite coffee places in town! By the time I need caffeine enough to shell out bucks, I need it fast, I need it strong, and I need it good at a reasonable price. Caffino hits all three and the service is the best Ive ever received at any coffee bar. And they absolutely do remember how you like your coffee if you come in with any regularity. These guys actually deserve their tips! Users/AlphaDog

20071203 01:14:01 nbsp This place is great, and it is true that if you are a regular they do remember how you like your coffee. Staff is really nice.  Users/RDN

20080731 11:55:57 nbsp Fast and friendly service.They make the best classic caramel latte. Users/jasen

20090331 22:24:50 nbsp This place is super convenient and awesome. Im so sad they reduced the hours. Where am I to go to perk up for late night studyoutings? Users/Hyungmi

20090504 16:53:16 nbsp Who say theyre good? I gave them 2 chances and both times were bad. Their coffee are so sweet, containing a lot of sugar and syrup. I also tried the smoothie, cant taste or see the fruit at all, again, only syrup and sugar. Will never come back. Users/thtly

20090815 23:54:17 nbsp You cant say anything bad about Caffino. Fast, friendly service, and the best coffee I can think of for the price. Overly pretentious coffee houses have nothing on them. A triple iced espresso will set you back a dollar and eighty five.  Users/Darren22

20100917 08:03:50 nbsp I dont know about them remembering which side of the cup you sip from, but when I used to ride past here on my way to school in the mornings, when they saw me coming theyd start making my mocha! I even got one or two free on days I wasnt going to buy one. Much love guys.  Great service, and excellent teas as well.   Users/MasonMurray

20101230 14:22:02 nbsp Hands down the BEST coffe shop in town! They serve amazing Mighty Leaf Teas!!! The staff members are always smiling and they remember all the regular customers drinks! I look forward to driving through every morning!!!  Users/CorkyCredit

20120914 19:37:16 nbsp It looks like they are gone  replaced by Dutch Brothers... Users/RMS
This is true Daubert
 anybody know why they closed? apparently other caffinos shut down as well. the phone line is disconnected, though.