module constants
    implicit none

    integer, parameter :: dp = kind(1.0d0)

    public dp

end module constants


text:
SUBROUTINE CALCONC (ISTL_,IS2TL_)

  ! *** SUBROUTINE CALCULATES THE CONCENTRATION OF DISSOLVED AND
  ! *** SUSPENDED CONSTITUTENTS, INCLUDING SALINITY, TEMPERATURE, DYE AND
  ! *** AND SUSPENDED SEDIMENT AT TIME LEVEL (N+1). THE VALUE OF ISTL
  ! *** INDICATES THE NUMBER OF TIME LEVELS IN THE STEP
  !
  ! CHANGE RECORD
  ! DATE MODIFIED     BY               DESCRIPTION
  !------------------------------------------------------------------------------------------------!
  !    2015-06       PAUL M. CRAIG     IMPLEMENTED SIGMA-Z (SGZ) IN EE7.3
  !    2015-01       PAUL M. CRAIG     Added fully coupled Ice Sub-model with Frazil Ice Transport
  !    2014-09       PAUL M. CRAIG     ADDED THE LWET BYPASS APPROACH
  !    2014-08       D H CHUNG         SET EXPLICIT PRECISIONS OF INTEGER & REAL
  !    2011-03       PAUL M. CRAIG     MERGED LATEST CODES AND RESTRUCTURED, ADDED OMP
  !    2002-05       John Hamrick      Modified calls to calbal and budget subroutines
  !                                     added calls to bal2t2, bal2t3
  !------------------------------------------------------------------------------------------------!

  USE GLOBAL
  USE OMP_LIB
  USE HEAT_MODULE, ONLY:CALHEAT

  IMPLICIT NONE

  INTEGER :: ISTL_,IS2TL_
  INTEGER :: K,L,LP,NT,NS,ND,IT,M
  INTEGER :: NTMP,LF,LL,LE,LN
  INTEGER, SAVE :: NICE

  REAL      :: RCDZKMK, RCDZKK
  REAL      :: DELTD2,CDYETMP,DAGE
  REAL,SAVE :: SEDTIME

  REAL(RKD), EXTERNAL :: DSTIME
  REAL(RKD)           :: CCUBTMP, CCMBTMP     ! VERTICAL DIFFUSION TEMPORARY VARIABLES
  REAL(RKD)           :: TTDS                 ! MODEL TIMING TEMPORARY VARIABLE

  ! *** VERTICAL DIFFUSION VARIABLES
  REAL(RKD),SAVE,ALLOCATABLE,DIMENSION(:)   :: CCLBTMP
  REAL(RKD),SAVE,ALLOCATABLE,DIMENSION(:)   :: EEB
  REAL(RKD),SAVE,ALLOCATABLE,DIMENSION(:,:) :: VCU

  REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: TOXASM
  REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: SEDASM
  REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: SNDASM

  IF(  .NOT. ALLOCATED(EEB) )THEN
    ALLOCATE(CCLBTMP(LCM))
    ALLOCATE(EEB(LCM))
    ALLOCATE(VCU(LCM,KCM))
    ALLOCATE(TOXASM(NTXM))
    ALLOCATE(SEDASM(NSCM))
    ALLOCATE(SNDASM(NSNM))
    EEB=0.0
    CCLBTMP=0.0
    VCU=0.
    TOXASM=0.0
    SEDASM=0.0
    SNDASM=0.0
    NICE=0
    SEDTIME=0.0
  ENDIF

  DELT=DT2
  IF( ISTL_ == 2 )THEN
    IF( ISDYNSTP == 0 )THEN
      DELT=DT
    ELSE
      DELT=DTDYN
    END IF
  ENDIF
  DELTD2=DELT

  ! *** MASS BALANCE
  IF( IS2TIM >= 1 )THEN
    IF( ISBAL >= 1 )THEN
      CALL BAL2T3A
    ENDIF
  ENDIF
  IT=1

  ! ***************************************************************************************
  ! *** 3D ADVECTI0N TRANSPORT CALCULATION, STANDARD TRANSPORT
  TTDS=DSTIME(0)

  !$OMP PARALLEL DEFAULT(SHARED)

  ! *** PRESPECIFY THE UPWIND CELLS FOR 3D ADVECTION
  !$OMP DO PRIVATE(ND,K,LP,L,LE,LN)
  DO ND=1,NDM
    DO K=1,KC
      DO LP=1,LLWET(K,ND)
        L=LKWET(LP,K,ND)
        IF( UHDY2(L,K) >= 0.0 )THEN
          LUPU(L,K)=LWC(L)
        ELSE
          LUPU(L,K)=L
        END IF
        IF( VHDX2(L,K) >= 0.0 )THEN
          LUPV(L,K)=LSC(L)
        ELSE
          LUPV(L,K)=L
        END IF
      ENDDO
    ENDDO

    IF( KC > 1 )THEN
      DO K=1,KS
        DO LP=1,LLWET(K,ND)
          L=LKWET(LP,K,ND)
          IF( W2(L,K) >= 0. )THEN
            KUPW(L,K)=K
          ELSE
            KUPW(L,K)=K+1
          END IF
        ENDDO
      ENDDO
    ENDIF
  ENDDO
  !$OMP END DO

  ! *** EE7.2 - MUST ZERO ALL THREADS FOR EVERY INSTANCE DO TO POTENTIAL OF DIFFERENT THREADS MAY NOT BE ZEROED FOR CURRENT TIME
  ! *** ZERO DRY CELL FLUXES
  IF( LADRY > 0 )THEN
    !$OMP SINGLE
    DO ND=1,NDM
      DO K=1,KC
        DO LP=1,LADRY
          L=LDRY(LP)
          FUHUD(L,K,ND)=0.
          FVHUD(L,K,ND)=0.
          FUHVD(L,K,ND)=0.
          FVHVD(L,K,ND)=0.
          UUUU(L,K,ND) =0.
          VVVV(L,K,ND) =0.
          DUU(L,K,ND)  =0.
          DVV(L,K,ND)  =0.
          POS(L,K,ND)  =0.
          WWWW(L,K,ND) =0.
          FWUU(L,K,ND) =0.

          LN=LNC(L)
          LE=LEC(L)
          FUHUD(LE,K,ND)=0.
          FVHUD(LN,K,ND)=0.
          FUHVD(LE,K,ND)=0.
          FVHVD(LN,K,ND)=0.
          UUUU(LE,K,ND) =0.
          VVVV(LN,K,ND) =0.
          DUU(LE,K,ND)  =0.
          DVV(LN,K,ND)  =0.
        ENDDO
      ENDDO
    ENDDO
    !$OMP END SINGLE
  ENDIF

  IF( ISTRAN(1) == 1 .AND. ISCDCA(1) < 4 )THEN
    !$OMP SINGLE PRIVATE(IT)
    !$  IT = OMP_GET_THREAD_NUM() + 1
    CALL CALTRAN(ISTL_,IS2TL_,1,1,SAL,SAL1,IT)
    !$OMP END SINGLE NOWAIT
  ENDIF

  IF( ISTRAN(2) == 1 .AND. ISCDCA(2) < 4 )THEN
    !$OMP SINGLE PRIVATE(IT)
    !$  IT = OMP_GET_THREAD_NUM() + 1
    CALL CALTRAN(ISTL_,IS2TL_,2,2,TEM,TEM1,IT)
    !$OMP END SINGLE NOWAIT

    IF( ISICE == 4 .AND. (IS2TL > 0 .OR. (IS2TL == 0 .AND. NCTBC /= NTSTBC)) )THEN
      IF( LFRAZIL )THEN
        ! *** FRAZIL ICE TRANSPORT
        !$OMP SINGLE PRIVATE(IT)
        !$  IT = OMP_GET_THREAD_NUM() + 1
        CALL CALTRANICE(ISTL_,FRAZILICE,FRAZILICE1,IT)
        NICE = NICE+1
        !$OMP END SINGLE NOWAIT

      ELSEIF( NICE > 0 )THEN
        ! *** ADVANCE THE FRAZIL ICE VARIABLE
        !$OMP SINGLE
        FRAZILICE1 = FRAZILICE
        NICE = 0
        !$OMP END SINGLE NOWAIT

      ENDIF
    ENDIF

  ENDIF

  IF( ISTRAN(3) == 1 .AND. ISCDCA(3) < 4 )THEN
    !$OMP SINGLE PRIVATE(IT)
    !$  IT = OMP_GET_THREAD_NUM() + 1
    CALL CALTRAN(ISTL_,IS2TL_,3,3,DYE,DYE1,IT)
    !$OMP END SINGLE NOWAIT
  ENDIF

  IF( ISTRAN(5) == 1 .AND. ISCDCA(5) < 4 )THEN
    !$OMP DO SCHEDULE(STATIC,1) PRIVATE(NT,M,IT)
    DO NT=1,NTOX
      M=MSVTOX(NT)
      !$  IT = OMP_GET_THREAD_NUM() + 1
      CALL CALTRAN(ISTL_,IS2TL_,5,M,TOX(1,1,NT),TOX1(1,1,NT),IT)
    ENDDO
    !$OMP END DO NOWAIT
  ENDIF

  IF( ISTRAN(6) == 1 .AND. ISCDCA(6) < 4 )THEN
    !$OMP DO SCHEDULE(STATIC,1) PRIVATE(NS,M,IT)
    DO NS=1,NSED
      M=MSVSED(NS)
      !$  IT = OMP_GET_THREAD_NUM() + 1
      CALL CALTRAN(ISTL_,IS2TL_,6,M,SED(1,1,NS),SED1(1,1,NS),IT)
    ENDDO
    !$OMP END DO NOWAIT
  ENDIF

  IF( ISTRAN(7) == 1 .AND. ISCDCA(7) < 4 )THEN
    !$OMP DO SCHEDULE(STATIC,1) PRIVATE(NS,M,IT)
    DO NS=1,NSND
      M=MSVSND(NS)
      !$  IT = OMP_GET_THREAD_NUM() + 1
      CALL CALTRAN(ISTL_,IS2TL_,7,M,SND(1,1,NS),SND1(1,1,NS),IT)
    ENDDO
    !$OMP END DO NOWAIT
  ENDIF
  !$OMP END PARALLEL
  TSADV=TSADV+(DSTIME(0)-TTDS)

  ! *** COSMIC TRANSPORT REMOVED 2014-07 PMC

  ! *** 1D ADVECTI0N TRANSPORT CALCULATION
  ! *** REMOVED 2004-09-19  PMC

  ! ******************************************************************************************
  ! *** VERTICAL DIFFUSION IMPLICIT HALF STEP CALCULATION
  IF( KC == 1 ) GOTO 1500
  TTDS=DSTIME(0)

  IF( LADRY > 0 )THEN
    DO LP=1,LADRY
      L=LDRY(LP)
      EEB(L)=1.
      CCLBTMP(L)=0.
    ENDDO

    DO K=1,KC
      DO LP=1,LADRY
        L=LDRY(LP)
        VCU(L,K)=0.
      ENDDO
    ENDDO
  ENDIF

  !$OMP PARALLEL DEFAULT(SHARED)

  !$OMP DO PRIVATE(ND,LF,LL,LP,L,K,CCUBTMP,CCMBTMP,RCDZKK,RCDZKMK,NT,NS) SCHEDULE(DYNAMIC,1)
  DO ND=1,NDM
    LF=(ND-1)*LDMWET+1
    LL=MIN(LF+LDMWET-1,LAWET)

    ! *** BOTTOM LAYER
    DO LP=1,LLWET(KS,ND)
      L=LKWET(LP,KS,ND)
      RCDZKK=-DELTD2*CDZKK(L,KSZ(L))
      CCUBTMP=RCDZKK*HPI(L)*AB(L,KSZ(L))
      CCMBTMP=1._8-CCUBTMP
      EEB(L)=1._8/CCMBTMP
      VCU(L,KSZ(L))=CCUBTMP*EEB(L)
    ENDDO
    IF( ISTRAN(1) >= 1 )THEN
      DO LP=1,LLWET(KS,ND)
        L=LKWET(LP,KS,ND)
        SAL(L,KSZ(L))=SAL(L,KSZ(L))*EEB(L)
      ENDDO
    ENDIF
    IF( ISTRAN(2) >= 1 )THEN
      DO LP=1,LLWET(KS,ND)
        L=LKWET(LP,KS,ND)
        TEM(L,KSZ(L))=TEM(L,KSZ(L))*EEB(L)
      ENDDO
    ENDIF
    IF( ISTRAN(3) >= 1 )THEN
      DO LP=1,LLWET(KS,ND)
        L=LKWET(LP,KS,ND)
        DYE(L,KSZ(L))=DYE(L,KSZ(L))*EEB(L)
      ENDDO
    ENDIF
    IF( ISTRAN(5) >= 1 )THEN
      DO NT=1,NTOX
        DO LP=1,LLWET(KS,ND)
          L=LKWET(LP,KS,ND)
          TOX(L,KSZ(L),NT)=TOX(L,KSZ(L),NT)*EEB(L)
        ENDDO
      ENDDO
    ENDIF
    IF( ISTRAN(6) >= 1 )THEN
      DO NS=1,NSED
        DO LP=1,LLWET(KS,ND)
          L=LKWET(LP,KS,ND)
          SED(L,KSZ(L),NS)=SED(L,KSZ(L),NS)*EEB(L)
        ENDDO
      ENDDO
    ENDIF
    IF( ISTRAN(7) >= 1 )THEN
      DO NS=1,NSND
        DO LP=1,LLWET(KS,ND)
          L=LKWET(LP,KS,ND)
          SND(L,KSZ(L),NS)=SND(L,KSZ(L),NS)*EEB(L)
        ENDDO
      ENDDO
    ENDIF

    ! *** MIDDLE LAYERS
    DO K=2,KS
      DO LP=1,LLWET(K-1,ND)
        L=LKWET(LP,K-1,ND)
        RCDZKMK=-DELTD2*CDZKMK(L,K)
        RCDZKK =-DELTD2*CDZKK(L,K)
        CCLBTMP(L)=RCDZKMK*HPI(L)*AB(L,K-1)
        CCUBTMP=RCDZKK*HPI(L)*AB(L,K)
        CCMBTMP=1._8-CCLBTMP(L)-CCUBTMP
        EEB(L)=1._8/(CCMBTMP-CCLBTMP(L)*VCU(L,K-1))
        VCU(L,K)=CCUBTMP*EEB(L)
      ENDDO
      IF( ISTRAN(1) >= 1 )THEN
        DO LP=1,LLWET(K-1,ND)
          L=LKWET(LP,K-1,ND)
          SAL(L,K)=(SAL(L,K)-CCLBTMP(L)*SAL(L,K-1))*EEB(L)
        ENDDO
      ENDIF
      IF( ISTRAN(2) >= 1 )THEN
        DO LP=1,LLWET(K-1,ND)
          L=LKWET(LP,K-1,ND)
          TEM(L,K)=(TEM(L,K)-CCLBTMP(L)*TEM(L,K-1))*EEB(L)
        ENDDO
      ENDIF
      IF( ISTRAN(3) >= 1 )THEN
        DO LP=1,LLWET(K-1,ND)
          L=LKWET(LP,K-1,ND)
          DYE(L,K)=(DYE(L,K)-CCLBTMP(L)*DYE(L,K-1))*EEB(L)
        ENDDO
      ENDIF
      IF( ISTRAN(5) >= 1 )THEN
        DO NT=1,NTOX
          DO LP=1,LLWET(K-1,ND)
            L=LKWET(LP,K-1,ND)
            TOX(L,K,NT)=(TOX(L,K,NT)-CCLBTMP(L)*TOX(L,K-1,NT))*EEB(L)
          ENDDO
        ENDDO
      ENDIF
      IF( ISTRAN(6) >= 1 )THEN
        DO NS=1,NSED
          DO LP=1,LLWET(K-1,ND)
            L=LKWET(LP,K-1,ND)
            SED(L,K,NS)=(SED(L,K,NS)-CCLBTMP(L)*SED(L,K-1,NS))*EEB(L)
          ENDDO
        ENDDO
      ENDIF
      IF( ISTRAN(7) >= 1 )THEN
        DO NS=1,NSND
          DO LP=1,LLWET(K-1,ND)
            L=LKWET(LP,K-1,ND)
            SND(L,K,NS)=(SND(L,K,NS)-CCLBTMP(L)*SND(L,K-1,NS))*EEB(L)
          ENDDO
        ENDDO
      ENDIF
    ENDDO

    ! *** TOP LAYER
    K=KC
    DO LP=1,LLWET(KS,ND)
      L=LKWET(LP,KS,ND)
      RCDZKMK=-DELTD2*CDZKMK(L,K)
      CCLBTMP(L)=RCDZKMK*HPI(L)*AB(L,K-1)
      CCMBTMP=1._8-CCLBTMP(L)
      EEB(L)=1._8/(CCMBTMP-CCLBTMP(L)*VCU(L,K-1))
    ENDDO
    IF( ISTRAN(1) >= 1 )THEN
      DO LP=1,LLWET(KS,ND)
        L=LKWET(LP,KS,ND)
        SAL(L,K)=(SAL(L,K)-CCLBTMP(L)*SAL(L,K-1))*EEB(L)
      ENDDO
    ENDIF
    IF( ISTRAN(2) >= 1 )THEN
      DO LP=1,LLWET(KS,ND)
        L=LKWET(LP,KS,ND)
        TEM(L,K)=(TEM(L,K)-CCLBTMP(L)*TEM(L,K-1))*EEB(L)
      ENDDO
    ENDIF
    IF( ISTRAN(3) >= 1 )THEN
      DO LP=1,LLWET(KS,ND)
        L=LKWET(LP,KS,ND)
        DYE(L,K)=(DYE(L,K)-CCLBTMP(L)*DYE(L,K-1))*EEB(L)
      ENDDO
    ENDIF
    IF( ISTRAN(5) >= 1 )THEN
      DO NT=1,NTOX
        DO LP=1,LLWET(KS,ND)
          L=LKWET(LP,KS,ND)
          TOX(L,K,NT)=(TOX(L,K,NT)-CCLBTMP(L)*TOX(L,K-1,NT))*EEB(L)
        ENDDO
      ENDDO
    ENDIF
    IF( ISTRAN(6) >= 1 )THEN
      DO NS=1,NSED
        DO LP=1,LLWET(KS,ND)
          L=LKWET(LP,KS,ND)
          SED(L,K,NS)=(SED(L,K,NS)-CCLBTMP(L)*SED(L,K-1,NS))*EEB(L)
        ENDDO
      ENDDO
    ENDIF
    IF( ISTRAN(7) >= 1 )THEN
      DO NS=1,NSND
        DO LP=1,LLWET(KS,ND)
          L=LKWET(LP,KS,ND)
          SND(L,K,NS)=(SND(L,K,NS)-CCLBTMP(L)*SND(L,K-1,NS))*EEB(L)
        ENDDO
      ENDDO
    ENDIF

    ! *** FINAL PASS
    DO K=KS,1,-1
      IF( ISTRAN(1) >= 1 )THEN
        DO LP=1,LLWET(K,ND)
          L=LKWET(LP,K,ND)
          SAL(L,K)=SAL(L,K)-VCU(L,K)*SAL(L,K+1)
        ENDDO
      ENDIF
      IF( ISTRAN(2) >= 1 )THEN
        DO LP=1,LLWET(K,ND)
          L=LKWET(LP,K,ND)
          TEM(L,K)=TEM(L,K)-VCU(L,K)*TEM(L,K+1)
        ENDDO
      ENDIF
      IF( ISTRAN(3) >= 1 )THEN
        DO LP=1,LLWET(K,ND)
          L=LKWET(LP,K,ND)
          DYE(L,K)=DYE(L,K)-VCU(L,K)*DYE(L,K+1)
        ENDDO
      ENDIF
      IF( ISTRAN(5) >= 1 )THEN
        DO NT=1,NTOX
          DO LP=1,LLWET(K,ND)
            L=LKWET(LP,K,ND)
            TOX(L,K,NT)=TOX(L,K,NT)-VCU(L,K)*TOX(L,K+1,NT)
          ENDDO
        ENDDO
      ENDIF
      IF( ISTRAN(6) >= 1 )THEN
        DO NS=1,NSED
          DO LP=1,LLWET(K,ND)
            L=LKWET(LP,K,ND)
            SED(L,K,NS)=SED(L,K,NS)-VCU(L,K)*SED(L,K+1,NS)
          ENDDO
        ENDDO
      ENDIF
      IF( ISTRAN(7) >= 1 )THEN
        DO NS=1,NSND
          DO LP=1,LLWET(K,ND)
            L=LKWET(LP,K,ND)
            SND(L,K,NS)=SND(L,K,NS)-VCU(L,K)*SND(L,K+1,NS)
          ENDDO
        ENDDO
      ENDIF
    ENDDO
  ENDDO     ! *** END OF DOMAIN
  !$OMP END DO
  !$OMP END PARALLEL

  TVDIF=TVDIF+(DSTIME(0)-TTDS)
  ! *** END OF VERTICAL DIFFUSION STEP

1500 CONTINUE

  ! ***  COMPUTE TOTALS FOR SEDIMENT TRANSPORT (AFTER ADVECTION/DIFFUSION)
  IF( ISTRAN(6) > 0 .OR. ISTRAN(7) > 0 )THEN
    !$OMP PARALLEL DEFAULT(SHARED)
    ! *** ZERO SED SEDIMENT ACCUMULATION ARRAYS
    IF( ISTRAN(6) > 0 )THEN
      !$OMP DO PRIVATE(ND,L,K,LF,LL,LP,NS) SCHEDULE(DYNAMIC)
      DO ND=1,NDM
        LF=(ND-1)*LDMWET+1
        LL=MIN(LF+LDMWET-1,LAWET)

        ! *** SEDIMENT BED  PMC - 2016-12-04 - NOT USED.  COMMENTED OUT
        !DO K=1,KB
        !  DO LP=LF,LL
        !    L=LWET(LP)
        !    SEDBT(L,K)=0.
        !  ENDDO
        !ENDDO

        ! *** WATER COLUMN
        DO K=1,KC
          DO LP=1,LLWET(K,ND)
            L=LKWET(LP,K,ND)
            SEDT(L,K) = 0.
          ENDDO
        ENDDO

        ! *** SEDIMENT BED  PMC - 2016-12-04 - NOT USED.  COMMENTED OUT
        !DO K=1,KB
        !  DO NS=1,NSED
        !    DO LP=LF,LL
        !      L=LWET(LP)
        !      SEDBT(L,K)=SEDBT(L,K)+SEDB(L,K,NS)
        !    ENDDO
        !  ENDDO
        !ENDDO
        DO NS=1,NSED
          DO K=1,KC
            DO LP=1,LLWET(K,ND)
              L=LKWET(LP,K,ND)
              SEDT(L,K) = SEDT(L,K) + SED(L,K,NS)
            ENDDO
          ENDDO
        ENDDO
      ENDDO
      !$OMP END DO
    ENDIF

    ! *** ZERO SND SEDIMENT ACCUMULATION ARRAYS
    IF( ISTRAN(7) > 0 )THEN
      !$OMP DO PRIVATE(ND,L,K,LF,LL,LP,NS) SCHEDULE(DYNAMIC)
      DO ND=1,NDM
        LF=(ND-1)*LDMWET+1
        LL=MIN(LF+LDMWET-1,LAWET)

        ! *** SEDIMENT BED  PMC - 2016-12-04 - NOT USED.  COMMENTED OUT
        !DO K=1,KB
        !  DO LP=LF,LL
        !    L=LWET(LP)
        !    SNDBT(L,K)=0.
        !  ENDDO
        !ENDDO

        ! *** WATER COLUMN
        DO K=1,KC
          DO LP=1,LLWET(K,ND)
            L=LKWET(LP,K,ND)
            SNDT(L,K) = 0.
          ENDDO
        ENDDO

        ! *** SEDIMENT BED  PMC - 2016-12-04 - NOT USED.  COMMENTED OUT
        !DO NS=1,NSND
        !  DO K=1,KB
        !    DO LP=LF,LL
        !      L=LWET(LP)
        !      SNDBT(L,K)=SNDBT(L,K)+SNDB(L,K,NS)
        !    ENDDO
        !  ENDDO
        !ENDDO

        DO NS=1,NSND
          DO K=1,KC
            DO LP=1,LLWET(K,ND)
              L=LKWET(LP,K,ND)
              SNDT(L,K) = SNDT(L,K) + SND(L,K,NS)
            ENDDO
          ENDDO
        ENDDO
      ENDDO
      !$OMP END DO
    ENDIF
    !$OMP END PARALLEL
  ENDIF

  ! ***************************************************************************************
  ! *** SURFACE AND INTERNAL HEAT SOURCE-SINK CALCULATION
  IF( ISTRAN(2) >= 1 )THEN
    TTDS=DSTIME(0)
    CALL CALHEAT(ISTL_)
    THEAT=THEAT+(DSTIME(0)-TTDS)
  ENDIF

  ! *** FULL IMPLICIT DYE AND TOXIC CONTAMINANT DECAY/GROWTH CALCULATION
  IF( ISTRAN(3) >= 1 )THEN
    IF( RKDYE == 1000.0 )THEN
      ! *** Age of Water
      DAGE=DELT/86400.
      !$OMP PARALLEL DO PRIVATE(ND,K,LP,L)
      DO ND=1,NDM

        DO K=1,KC
          DO LP=1,LLWET(K,ND)
            L=LKWET(LP,K,ND)
            DYE(L,K)=DYE(L,K)+DAGE
          ENDDO
        ENDDO
      ENDDO
      !$OMP END PARALLEL DO

    ELSEIF( RKDYE /= 0.0 )THEN
      IF( RKDYE < 0.0 )THEN
        CDYETMP=EXP(-RKDYE*DELT)
      ELSE
        CDYETMP=1./(1.+DELT*RKDYE)
      ENDIF

      !$OMP PARALLEL DO PRIVATE(ND,K,LP,L)
      DO ND=1,NDM

        DO K=1,KC
          DO LP=1,LLWET(K,ND)
            L=LKWET(LP,K,ND)
            DYE(L,K)=CDYETMP*DYE(L,K)
          ENDDO
        ENDDO
      ENDDO
      !$OMP END PARALLEL DO
    ENDIF
  ENDIF

  ! **************************************************************************************************
  ! *** BOTTOM AND INTERNAL SEDIMENT AND TOXIC CONTAMINATION SOURCE-SINK CALCULATION

  ! *** SEDIMENT AND TOXICS SETTLING,DEPOSITION,RESUSPENSION,ETC
  IF( ( ISTRAN(6) >= 1 .OR. ISTRAN(7) >= 1 ) .AND. TIMEDAY >= SEDSTART )THEN
    IF( IS2TIM >= 1 )THEN
      ! *** FOR TWO TIME LEVEL SIMULATION
      SEDTIME = SEDTIME + DELT
      IF( SEDTIME >= SEDSTEP )THEN
        DTSED = SEDTIME
        CALL SSEDTOX
        SEDTIME = 0.0
      ENDIF
    ELSE
      ! *** FOR THREE TIME LEVEL SIMULATION
      IF( NCTBC == 1 )THEN
        DTSED = FLOAT(NTSTBC)*DT
        CALL SSEDTOX
        SEDTIME = 0.0
      ENDIF
    ENDIF
  ENDIF

  ! *** OPTIONAL MASS BALANCE CALCULATION
  IF( IS2TIM == 0 )THEN
    IF( ISTL_/=2 .AND. ISBAL >= 1 )THEN
      CALL CALBAL2
      CALL CALBAL3
      NTMP=MOD(N,2)
      IF( NTMP == 0 )THEN
        CALL CBALEV2
        CALL CBALEV3
      ELSE
        CALL CBALOD2
        CALL CBALOD3
      ENDIF
    ENDIF
  ENDIF

  ! *** CALLS TO TWO-TIME LEVEL BALANCES
  !
  IF( IS2TIM >= 1 )THEN
    IF( ISBAL >= 1 )THEN
      CALL BAL2T2
      CALL BAL2T3B(1)
    ENDIF
  ENDIF

  ! *** SEDIMENT BUDGET CALCULATION    (DLK 10/15)
  !
  IF( IS2TIM == 0 )THEN
    IF( ISTL_/=2 .AND. ISSBAL >= 1 )THEN
      CALL BUDGET2
      CALL BUDGET3
    ENDIF
  ENDIF

  ! *** 2014-09 - REMOVED DATA ASSIMILATION CODE (PMC)

  RETURN

END

max_stars_repo_name: dsi-llc/EFDCPlus

max_forks_repo_name: dsi-llc/EFDCPlus


text:
      SUBROUTINE DGTTS2( ITRANS, N, NRHS, DL, D, DU, DU2, IPIV, B, LDB )
*
*  -- LAPACK auxiliary routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            ITRANS, LDB, N, NRHS
*     ..
*     .. Array Arguments ..
      INTEGER            IPIV( * )
      DOUBLE PRECISION   B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * )
*     ..
*
*  Purpose
*  =======
*
*  DGTTS2 solves one of the systems of equations
*     A*X = B  or  A'*X = B,
*  with a tridiagonal matrix A using the LU factorization computed
*  by DGTTRF.
*
*  Arguments
*  =========
*
*  ITRANS  (input) INTEGER
*          Specifies the form of the system of equations.
*          = 0:  A * X = B  (No transpose)
*          = 1:  A'* X = B  (Transpose)
*          = 2:  A'* X = B  (Conjugate transpose = Transpose)
*
*  N       (input) INTEGER
*          The order of the matrix A.
*
*  NRHS    (input) INTEGER
*          The number of right hand sides, i.e., the number of columns
*          of the matrix B.  NRHS >= 0.
*
*  DL      (input) DOUBLE PRECISION array, dimension (N-1)
*          The (n-1) multipliers that define the matrix L from the
*          LU factorization of A.
*
*  D       (input) DOUBLE PRECISION array, dimension (N)
*          The n diagonal elements of the upper triangular matrix U from
*          the LU factorization of A.
*
*  DU      (input) DOUBLE PRECISION array, dimension (N-1)
*          The (n-1) elements of the first super-diagonal of U.
*
*  DU2     (input) DOUBLE PRECISION array, dimension (N-2)
*          The (n-2) elements of the second super-diagonal of U.
*
*  IPIV    (input) INTEGER array, dimension (N)
*          The pivot indices; for 1 <= i <= n, row i of the matrix was
*          interchanged with row IPIV(i).  IPIV(i) will always be either
*          i or i+1; IPIV(i) = i indicates a row interchange was not
*          required.
*
*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
*          On entry, the matrix of right hand side vectors B.
*          On exit, B is overwritten by the solution vectors X.
*
*  LDB     (input) INTEGER
*          The leading dimension of the array B.  LDB >= max(1,N).
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            I, IP, J
      DOUBLE PRECISION   TEMP
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( N.EQ.0 .OR. NRHS.EQ.0 )
     $   RETURN
*
      IF( ITRANS.EQ.0 ) THEN
*
*        Solve A*X = B using the LU factorization of A,
*        overwriting each right hand side vector with its solution.
*
         IF( NRHS.LE.1 ) THEN
            J = 1
   10       CONTINUE
*
*           Solve L*x = b.
*
            DO 20 I = 1, N - 1
               IP = IPIV( I )
               TEMP = B( I+1-IP+I, J ) - DL( I )*B( IP, J )
               B( I, J ) = B( IP, J )
               B( I+1, J ) = TEMP
   20       CONTINUE
*
*           Solve U*x = b.
*
            B( N, J ) = B( N, J ) / D( N )
            IF( N.GT.1 )
     $         B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
     $                       D( N-1 )
            DO 30 I = N - 2, 1, -1
               B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
     $                     B( I+2, J ) ) / D( I )
   30       CONTINUE
            IF( J.LT.NRHS ) THEN
               J = J + 1
               GO TO 10
            END IF
         ELSE
            DO 60 J = 1, NRHS
*
*              Solve L*x = b.
*
               DO 40 I = 1, N - 1
                  IF( IPIV( I ).EQ.I ) THEN
                     B( I+1, J ) = B( I+1, J ) - DL( I )*B( I, J )
                  ELSE
                     TEMP = B( I, J )
                     B( I, J ) = B( I+1, J )
                     B( I+1, J ) = TEMP - DL( I )*B( I, J )
                  END IF
   40          CONTINUE
*
*              Solve U*x = b.
*
               B( N, J ) = B( N, J ) / D( N )
               IF( N.GT.1 )
     $            B( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) /
     $                          D( N-1 )
               DO 50 I = N - 2, 1, -1
                  B( I, J ) = ( B( I, J )-DU( I )*B( I+1, J )-DU2( I )*
     $                        B( I+2, J ) ) / D( I )
   50          CONTINUE
   60       CONTINUE
         END IF
      ELSE
*
*        Solve A' * X = B.
*
         IF( NRHS.LE.1 ) THEN
*
*           Solve U'*x = b.
*
            J = 1
   70       CONTINUE
            B( 1, J ) = B( 1, J ) / D( 1 )
            IF( N.GT.1 )
     $         B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
            DO 80 I = 3, N
               B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-DU2( I-2 )*
     $                     B( I-2, J ) ) / D( I )
   80       CONTINUE
*
*           Solve L'*x = b.
*
            DO 90 I = N - 1, 1, -1
               IP = IPIV( I )
               TEMP = B( I, J ) - DL( I )*B( I+1, J )
               B( I, J ) = B( IP, J )
               B( IP, J ) = TEMP
   90       CONTINUE
            IF( J.LT.NRHS ) THEN
               J = J + 1
               GO TO 70
            END IF
*
         ELSE
            DO 120 J = 1, NRHS
*
*              Solve U'*x = b.
*
               B( 1, J ) = B( 1, J ) / D( 1 )
               IF( N.GT.1 )
     $            B( 2, J ) = ( B( 2, J )-DU( 1 )*B( 1, J ) ) / D( 2 )
               DO 100 I = 3, N
                  B( I, J ) = ( B( I, J )-DU( I-1 )*B( I-1, J )-
     $                        DU2( I-2 )*B( I-2, J ) ) / D( I )
  100          CONTINUE
               DO 110 I = N - 1, 1, -1
                  IF( IPIV( I ).EQ.I ) THEN
                     B( I, J ) = B( I, J ) - DL( I )*B( I+1, J )
                  ELSE
                     TEMP = B( I+1, J )
                     B( I+1, J ) = B( I, J ) - DL( I )*TEMP
                     B( I, J ) = TEMP
                  END IF
  110          CONTINUE
  120       CONTINUE
         END IF
      END IF
*
*     End of DGTTS2
*
      END

ext: f

max_stars_repo_name: pkeuchel/oomph-lib

max_issues_repo_name: pkeuchel/oomph-lib

max_forks_repo_name: pkeuchel/oomph-lib


SUBROUTINE solve_band(pb,work,loads)
!
! This subroutine performs Gaussian forward and back-substitution
! on the reduced unsymmetric band matrix pb.
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::pb(:,:),work(:,:)
 REAL(iwp),INTENT(OUT)::loads(0:)
 INTEGER::iwp1,n,n1,i,iv,l,iq,iv1
 REAL(iwp)::s,pt5=0.5_iwp
 iwp1=(UBOUND(pb,2)-1)/2+1
 n=UBOUND(pb,1)
 iq=2*iwp1-1
 n1=n-1
 DO iv=1,n1
   i=INT(work(iwp1,iv)+pt5)
   IF(i/=iv)THEN
     s=loads(iv)
     loads(iv)=loads(i)
     loads(i)=s
   END IF
   l=iv+iwp1-1
   IF(l>n)l=n
   iv1=iv+1
   DO i=iv1,l
     loads(i)=loads(i)-work(i-iv,iv)*loads(iv)
   END DO
 END DO
 loads(n)=loads(n)/pb(n,1)
 iv=n-1
 DO WHILE(iv/=0)
   s=loads(iv)
   l=iq
   IF(iv+l-1>n)l=n-iv+1
   DO i=2,l
     s=s-pb(iv,i)*loads(iv+i-1)
     loads(iv)=s/pb(iv,1)
   END DO
 iv=iv-1
 END DO
RETURN
END SUBROUTINE solve_band


text:
So, youre in IRC!  Here are a few tips:

  Why isnt anybody answering me?   IRC is real time chat, but often people have the window open on their computer as they work through the day. Sometimes you get an instant reply, and sometimes it can be an hour or longer before somebody replies.  In other words, its vaguely like sending a text message to somebody.
  You can register your user name (your nick) with freenode.net so that other people cant use it.  First, make sure you have your nick set to the name you want to register.  You can type /nick xxxx to change it.  In your freenode.net tab/window, type /msg NickServ register <password> <email>.  For example, /msg NickServ register thisismypassword blahblah@email.whatever will register your current Nick with the password thisismypassword and send a confirmation email to the email you entered.  In the future, when you log in using that user name, youll be asked to authenticate.  As youll be prompted by NickServ, you just type /msg NickServ identify thisismypassword, and youre logged in and good to go!

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c=======================================================================
c
c     Linear Algebra - Level 2 BLAS Routines Interfaces
c
c-----------------------------------------------------------------------
c
c     FUNCTION LIST
c
c     (P)
c       PEVEV   : vector product,        z(n) <- x(n).*y(n)
c       PEVJEV  : vector product,        z(n) <- x(n)./y(n)
c       PMLV    : matrix vector product, y(n) <- A(n,n)*x(n) with A low tringular
c       PMTV    : matrix vector product, y(m) <- A'(n,m)*x(n)
c       PMV     : matrix vector product, y(n) <- A(n,m)*x(m)
c       PMVV    : matrix vector product, z(n) <- A(n,m)*x(m) + y(n)
c       PMVX    : matrix vector product, y(n)  <- alpha*A(n,m)*x(m)
c
c       PV      : product vector/vector, A(n,m) <- x(n)*y(m)'
c       PV1     : product vector/vector, A(n,n) <- x(n)*x(n)'
c
c=======================================================================
c
c     subroutine PEVEV
c
c     vector product, z(n) <- x(n).*y(n)
c
c-----------------------------------------------------------------------
      SUBROUTINE PEVEV (n, x, y, z)
c-----------------------------------------------------------------------
c
c     INPUT
c       n : vector size                                          integer
c       x : vector (n)                                            double
c       y : vector (n)                                            double
c
c     OUTPUT
c       z : vector (n)                                            double
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
      INTEGER n, i
      DOUBLE PRECISION x(*), y(*), z(*)
c
      DO i = 1,n
        z(i) = x(i)*y(i)
      ENDDO
      RETURN
      END
c
c=======================================================================
c
c     subroutine PEVJEV
c
c     vector product, z(n) <- x(n)./y(n)
c
c-----------------------------------------------------------------------
      SUBROUTINE PEVJEV (n, x, y, z, info)
c-----------------------------------------------------------------------
c
c     INPUT
c       n    : vector size                                       integer
c       x    : vector (nvect)                                     double
c       y    : vector (nvect)                                     double
c
c     OUTPUT
c       z    : vector (n)                                         double
c       info : diagnostic argument                               integer
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
      INTEGER n, info, i
      DOUBLE PRECISION EPS, x(*), y(*), z(*)
      PARAMETER ( EPS = 1.0E-30 )
c
      info = 0
      DO i = 1,n
         IF ( ABS(y(i)) .GT. EPS ) THEN
            z(i) = x(i) / y(i)
         ELSE
            info = -1
         ENDIF
      ENDDO
      RETURN
      END
c
c=======================================================================
c
c     subroutine PMLV
c
c     matrix vector product, y(n) <- A(n,n)*x(n) with A low tringular
c
c-----------------------------------------------------------------------
      SUBROUTINE PMLV ( n, A, x, y )
c-----------------------------------------------------------------------
c
c     INPUT
c       n : matrix size                                          integer
c       A : matrix (n*n)                                          double
c       x : vector (n)                                            double
c
c     OUTPUT
c       y : vector (n)                                            double
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
      INTEGER n, incz
      DOUBLE PRECISION A(n,*), x(*), y(*)
      PARAMETER (incz=1)
c
c     copy y <- x
      CALL YV ( n, x, y )
c
c     y <- A*y, A low triangular
      CALL DTRMV( 'L', 'N', 'N', n, A, n, y, incz)
      RETURN
      END
c
c=======================================================================
c
c     subroutine PMTV
c
c     matrix vector product, y(m) <- A'(n,m)*x(n)
c
c-----------------------------------------------------------------------
      SUBROUTINE PMTV ( n, m, A, x, y )
c-----------------------------------------------------------------------
c
c     INPUT
c       n  : number of rows of the matrix A                      integer
c       m  : number of columns of the matrix A                   integer
c       A  : matrix (n*m)                                         double
c       x  : vector (n)                                           double
c
c     OUTPUT
c       y  : vector(m) = A*x                                      double
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
      INTEGER n, m
      DOUBLE PRECISION A(n,*), x(*), y(*)
      INTEGER incx, incy
      DOUBLE PRECISION alpha, beta
      PARAMETER( incx=1, incy=1 )
      PARAMETER( alpha=1.0d0, beta=0.0d0 )
c
      CALL DGEMV ( 'T', n, m, alpha, A, n, x, incx, beta, y, incy)
      RETURN
      END
c
c=======================================================================
c
c     subroutine PMV
c
c     matrix vector product, y(n) <- A(n,m)*x(m)
c
c-----------------------------------------------------------------------
      SUBROUTINE PMV ( n, m, A, x, y )
c-----------------------------------------------------------------------
c
c     INPUT
c       n  : number of rows of the matrix A                      integer
c       m  : number of columns of the matrix A                   integer
c       A  : matrix (n*m)                                         double
c       x  : vector (m)                                           double
c
c     OUTPUT
c       y  : vector(n) = A*x (n)                                  double
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
      INTEGER n, m
      DOUBLE PRECISION A(n,*), x(*), y(*)
      INTEGER incx, incy
      DOUBLE PRECISION alpha, beta
      PARAMETER( incx=1, incy=1 )
      PARAMETER( alpha=1.0d0, beta=0.0d0 )
c
      CALL DGEMV ( 'N', n, m, alpha, A, n, x, incx, beta, y, incy)
      RETURN
      END
c
c=======================================================================
c
c     subroutine PMVV
c
c     matrix vector multiply: z(n) <- A(n,m)*x(m) + y(n)
c
c-----------------------------------------------------------------------
      SUBROUTINE PMVV ( n, m, A, x, y, z )
c-----------------------------------------------------------------------
c
c     INPUT
c       n : number of rows of the matrix A                       integer
c       m : number of columns of the matrix A                    integer
c       A : matrix (n*m)                                          double
c       x : vector (m)                                            double
c       y : vector (n)                                            double
c
c     OUTPUT
c       z : vector (n)                                            double
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
      INTEGER n, m
      DOUBLE PRECISION A(n,*), x(*), y(*), z(*)
      INTEGER incx, incz
      DOUBLE PRECISION alpha, beta
      PARAMETER( incx=1, incz=1 )
      PARAMETER( alpha=1.0d0, beta=1.0d0 )
c
c     copies z <- y
      CALL YV ( n, y, z )
c
c     z = A*x + z
      CALL DGEMV ( 'N', n, m, alpha, A, n, x, incx, beta, z, incz)
      RETURN
      END
c=======================================================================
c
c     subroutine PMVX
c
c     matrix vector product, y(n) <- alpha*A(n,m)*x(m)
c
c-----------------------------------------------------------------------
      SUBROUTINE PMVX(n, m, A, x, alpha, y, dwork, info)
c-----------------------------------------------------------------------
c
c     INPUT
c       n     : number of rows                                   integer
c       m     : number of columns                                integer
c       A     : matrix (n*m)                                      double
c       x     : vector (m)                                        double
c       alpha : scalar                                            double
c
c     WORKSPACE
c       dwork : m
c
c     OUTPUT
c       y     : vector (n)                                        double
c       info  : diagnostic argument                              integer
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
c
c     arguments i/o
      INTEGER n, m, info
      DOUBLE PRECISION alpha, A(n,*), x(*), y(*)
c
c     workspaces
      DOUBLE PRECISION dwork(*)
c
c     dwork(m) <- alpha*x(m)
      CALL PVX(m, x, alpha, dwork(1))
c
c     y(n) <- A(n,m)*dwork(m)
      CALL PMV(n, m, A, dwork(1), y)
      info = 0
      RETURN
      END
c
c=======================================================================
c
c     subroutine PV
c
c     product vector/vector A(n,m) <- x(n)*y(m)'
c
c-----------------------------------------------------------------------
      subroutine PV ( n, m, x, y, A )
c-----------------------------------------------------------------------
c
c     INPUT
c       n  : size of vector x                                    integer
c       m  : size of vector y                                    integer
c       x  : column vector (n)                                    double
c       y  : raw vector (m)                                       double
c
c     OUTPUT
c       A  : matrix (n,m)                                         double
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
      INTEGER n, m, incx, incy
      DOUBLE PRECISION alpha, ZERO, x(*), y(*), A(n,*)
      PARAMETER (alpha=1.0, incx=1, incy=1, ZERO = 0.0 )
c
c     A <- 0
      CALL PMX2 ( n, m, A, ZERO)
c
c     A <- alpha*x*y' + A
      CALL DGER(n, m, alpha, x, incx, y, incy, A, n)
      RETURN
      END
c
c=======================================================================
c
c     subroutine PV1
c
c     product vector/vector, A(n,n) <- x(n)*x(n)'
c
c-----------------------------------------------------------------------
      SUBROUTINE PV1 ( n, x, A )
c-----------------------------------------------------------------------
c
c     INPUT
c       n : vector size                                          integer
c       x : column vector (n)                                     double
c
c     OUTPUT
c       A : matrix (n,n)                                          double
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
      INTEGER n
      DOUBLE PRECISION x(*), A(n,*)
c
      CALL PV ( n, n, x, x, A )
      RETURN
      END

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


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

!> Compute with time quantities
!>
!> The advance_time program computes the resulting time when either
!> adding or subtracting time intervals.  The increments can be
!> expressed in days, hours, minutes or seconds.  The output can be
!> formatted as native WRF, CESM, Julian or Gregorian format.
!>
!> Reads input from standard input to be more portable, since older
!> versions of iargc() weren't standardized.
!>
!> Based on the WRF da_advance_cymdh utility.
!>
!> All time computations call DART time manager.
!>
!>   - has accuracy down to second,
!>   - can use day/hour/minute/second (with/without +/- sign) to advance time,
!>   - can digest various input date format if it still has the right order (ie. cc yy mm dd hh nn ss)
!>   - can digest flexible time increment
!>   - can output in wrf date format (ccyy-mm-dd_hh:nn:ss)
!>   - can specify output date format
!>   - can output Julian day
!>   - can output Gregorian days and seconds (since year 1601)
!>   - can output in CESM time format (ccyy-mm-dd-fffff where fffff is seconds of day)
!>
!> Examples
!>
!> - advance 12 h
!>
!>    echo 20070730      12         | advance_time
!>
!> - go back 1 day 2 hours 30 minutes and 30 seconds
!>
!>    echo 2007073012   -1d2h30m30s | advance_time
!>
!> - go back 3 hours 30 minutes less 1 second
!>
!>    echo 2007073012    1s-3h30m   | advance_time
!>
!> - advance 2 days and 1 second, output in wrf date format (three ways)
!>
!>    echo 200707301200        2d1s -w | advance_time
!>
!>    echo 2007-07-30_12:00:00 2d1s -w | advance_time
!>
!>    echo 200707301200        2d1s -f ccyy-mm-dd_hh:nn:ss | advance_time
!>
!> - advance 120 h, and print year and Julian day
!>
!>    echo 2007073006    120 -j     | advance_time
!>
!> - advance 120 h, print year, Julian day, hour, minute and second
!>
!>    echo 2007073006    120 -J     | advance_time
!>
!> - print Gregorian day and second (since year 1601)
!>
!>    echo 2007073006    0 -g       | advance_time
!>
!> - print CESM format time (ccyy-mm-dd-fffff where fffff is sec of day)
!>
!>    echo 2007073006    0 -c       | advance_time
!>
!> @todo if run with no inputs ... it just hangs. Can we make it fail straight away?

program advance_time

use time_manager_mod, only : time_type, set_calendar_type, GREGORIAN, &
                             increment_time, decrement_time, get_time, &
                             set_date, get_date, julian_day
use    utilities_mod, only : initialize_utilities, E_ERR, error_handler, &
                             finalize_utilities
use   parse_args_mod, only : get_args_from_string

implicit none

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

integer :: ccyy, mm, dd, hh, nn, ss, dday, dh, dn, ds, gday, gsec
integer :: nargum, i
character(len=80), dimension(10) :: argum
character(len=14) :: ccyymmddhhnnss
character(len=80) :: out_date_format, dtime
character(len=256) :: in_string
integer :: datelen
integer, parameter :: stdout=6
type(time_type) :: base_time

character(len=512) :: string1, string2

! Initialize modules used that require it, and be silent about it
call initialize_utilities('advance_time', output_flag = .false.)

call set_calendar_type(GREGORIAN)

! this routine reads a line from standard input and parses it up
! into blank-separated words.
read(*, '(A)') in_string
call get_args_from_string(in_string, nargum, argum)

if ( nargum < 2 ) then
   write(unit=stdout, fmt='(a)') &
   'Usage:   echo ccyymmddhh[nnss] [+|-]dt[d|h|m|s] [-w|-W|-wrf|-WRF] [-f|-F date_format] [-j|-J] [-g|-G] [-c|-C] | advance_time'
   write(unit=stdout, fmt='(a)') &
      'Option:  -w|-W|-wrf|-WRF  output in wrf date format as ccyy-mm-dd_hh:nn:ss'
   write(unit=stdout, fmt='(a)') &
      '         -f|-F  specify output date format, such as ccyy-mm-dd_hh:nn:ss, or ''ccyy/mm/dd  hh:nn:ss'''
   write(unit=stdout, fmt='(a)') &
      '         -j|-J  print Julian day'
   write(unit=stdout, fmt='(a)') &
      '         -g|-G  print Gregorian days and seconds (since year 1601)'
   write(unit=stdout, fmt='(a)') &
      'Example: echo 20070730      12         | advance_time   # advance 12 h'
   write(unit=stdout, fmt='(a)') &
      '         echo 2007073012   -1d2h30m30s | advance_time   # back 1 day 2 hours 30 min and 30 sec'
   write(unit=stdout, fmt='(a)') &
      '         echo 2007073012    1s-3h30m   | advance_time   # back 3 hours 30 minutes less 1 second'
   write(unit=stdout, fmt='(a)') &
      '         echo 200707301200  1d1s -w    | advance_time   # advance 1 day 1 sec, output in wrf date format'
   write(unit=stdout, fmt='(a)') &
      '         echo 2007-07-30_12:00:00 2d1s -w | advance_time              # same as previous example'
   write(unit=stdout, fmt='(a)') &
      '         echo 200707301200  2d1s -f ccyy-mm-dd_hh:nn:ss | advance_time # same as previous'
   write(unit=stdout, fmt='(a)') &
      '         echo 2007073006    120 -j     | advance_time    # advance 120 h, and print year and Julian day'
   write(unit=stdout, fmt='(a)') &
      '         echo 2007073006    120 -J     | advance_time    # advance 120 h, print year, Julian day, hour, minute and second'
   write(unit=stdout, fmt='(a)') &
      '         echo 2007073006    0 -g       | advance_time    # print Gregorian day and second (since year 1601)'
   write(unit=stdout, fmt='(a)') &
      '         echo 2007073006    0 -c       | advance_time    # print CESM format time (ccyy-mm-dd-fffff where fffff is sec of day)'
   write(unit=stdout, fmt='(a)') ''

   call error_handler(E_ERR,'advance_time','Invalid Usage', source, revision, revdate)
end if

ccyymmddhhnnss = parsedate(argum(1))
datelen = len_trim(ccyymmddhhnnss)

if (datelen == 8) then
   read(ccyymmddhhnnss(1:10), fmt='(i4, 2i2)')  ccyy, mm, dd
   hh = 0
   nn = 0
   ss = 0
else if (datelen == 10) then
   read(ccyymmddhhnnss(1:10), fmt='(i4, 3i2)')  ccyy, mm, dd, hh
   nn = 0
   ss = 0
else if (datelen == 12) then
   read(ccyymmddhhnnss(1:12), fmt='(i4, 4i2)')  ccyy, mm, dd, hh, nn
   ss = 0
else if (datelen == 14) then
   read(ccyymmddhhnnss(1:14), fmt='(i4, 5i2)')  ccyy, mm, dd, hh, nn, ss
elseif (datelen == 13) then
   read(ccyymmddhhnnss(1:13), fmt='(i4, 2i2, i5)')  ccyy, mm, dd, ss
   if (ss >= 86400) then
      write(string1,*)'seconds-of-day is ',ss,' as parsed from ',trim(ccyymmddhhnnss)
      write(string2,*)'seconds-of-day  must be less than 86400'
      call error_handler(E_ERR,'advance_time',string1, source, revision, revdate, text2=string2)
   endif
   hh = ss / 3600
   ss = ss - hh * 3600
   nn = ss / 60
   ss = ss - nn * 60
else
   write(string1,*)'unsupported format for ',trim(argum(1))
   call error_handler(E_ERR,'advance_time',string1, source, revision, revdate)
endif

base_time = set_date(ccyy, mm, dd, hh, nn, ss)

dtime = trim(argum(2))
call parsedt(dtime,dday,dh,dn,ds)

!print*, 'delta t: ', dday, dh, dn, ds

! each part can be positive or negative, or 0.
if (dday > 0) then
   base_time = increment_time(base_time, 0, dday)
else if (dday < 0) then
   base_time = decrement_time(base_time, 0, -dday)
endif

if (dh > 0) then
   base_time = increment_time(base_time, dh*3600)
else if (dh < 0) then
   base_time = decrement_time(base_time, -dh*3600)
endif

if (dn > 0) then
   base_time = increment_time(base_time, dn*60)
else if (dn < 0) then
   base_time = decrement_time(base_time, -dn*60)
endif

if (ds > 0) then
   base_time = increment_time(base_time, ds)
else if (ds < 0) then
   base_time = decrement_time(base_time, -ds)
endif

call get_date(base_time, ccyy, mm, dd, hh, nn, ss)

write(ccyymmddhhnnss(1:14), fmt='(i4, 5i2.2)')  ccyy, mm, dd, hh, nn, ss
if ( nargum == 2 ) then
   if (datelen == 13) datelen=10
   if (datelen < 14) then
      if(nn /= 0) datelen=12
      if(ss /= 0) datelen=14
   endif
   write(unit=stdout, fmt='(a)') ccyymmddhhnnss(1:datelen)
else if ( nargum > 2 ) then
   i = 3
   do while (i <= nargum)
     select case ( trim(argum(i)) )
        case ('-w', '-W', '-wrf','-WRF')
           out_date_format = 'ccyy-mm-dd_hh:nn:ss'
           write(unit=stdout, fmt='(a)') trim(formatdate(ccyymmddhhnnss, out_date_format))
           i = i+1
        case ('-f', '-F')
           out_date_format = trim(argum(i+1))
           write(unit=stdout, fmt='(a)') trim(formatdate(ccyymmddhhnnss, out_date_format))
           i = i+2
        case ('-j')
           write(unit=stdout, fmt='(I4,I4)') ccyy, julian_day(ccyy,mm,dd)
           i = i+1
        case ('-J')
           write(unit=stdout, fmt='(I4,I4,I3,I3,I3)') ccyy, julian_day(ccyy,mm,dd),hh,nn,ss
           i = i+1
        case ('-g','-G')
           call get_time(base_time, gsec, gday)
           write(unit=stdout, fmt='(I8,I8)') gday, gsec
           i = i+1
        case ('-c','-C')
           out_date_format = 'ccyy-mm-dd-fffff'
           write(unit=stdout, fmt='(a)') trim(formatCESMdate(ccyy,mm,dd,hh,nn,ss))
           i = i+1
        case default
           i = i+1
     end select
   end do
end if

call finalize_utilities()

contains

!-----------------------------------------------------------------------
!> removes non-numeric characters from a date string
!>
!> @param[in] datein character string containing the date. May include
!>               dashes, colons, etc.
!> @return character string with only numeric characters

function parsedate(datein)
character(len=*), intent(in) :: datein

character(len=14) :: parsedate
character(len=1 ) :: ch
integer :: n, i

parsedate = '00000000000000'

i=0
do n = 1, len_trim(datein)
   ch = datein(n:n)
   if (ch >= '0' .and. ch <= '9') then
      i=i+1
      parsedate(i:i)=ch
   end if
end do

if (i == 13) then
   parsedate(14:14) = ''
   return  ! CESM format
else if (parsedate(11:14) == '0000') then
   parsedate(11:14) = ''
else if(parsedate(13:14) == '00') then
   parsedate(13:14) = ''
end if

end function parsedate

!-----------------------------------------------------------------------
!> extracts the day,hour,minutes and seconds from the second input argument
!>
!> @param[in] dt character string with the temporal offset (the second input argument).
!> @param[out] dday the day
!> @param[out] dh the hour
!> @param[out] dn the minute
!> @param[out] ds the second

subroutine parsedt(dt,dday,dh,dn,ds)
character(len=*), intent(in)  :: dt
integer,          intent(out) :: dday
integer,          intent(out) :: dh
integer,          intent(out) :: dn
integer,          intent(out) :: ds

character(len=1 ) :: ch
integer :: n,i,d,s,nounit

! initialize time and sign
nounit=1
dday=0
dh=0
dn=0
ds=0
d=0
s=1

do n = 1, len_trim(dt)
   ch = dt(n:n)
   select case (ch)
      case ('0':'9')
        read(ch,fmt='(i1)') i
        d=d*10+i
      case ('-')
        s=-1
      case ('+')
        s=1
      case ('d')
        nounit=0
        dday=dday+d*s
        d=0
      case ('h')
        nounit=0
        dh=dh+d*s
        d=0
      case ('n','m')
        nounit=0
        dn=dn+d*s
        d=0
      case ('s')
        nounit=0
        ds=ds+d*s
        d=0
      case default
        continue
   end select
end do

if (nounit==1) dh=d*s

end subroutine parsedt

!-----------------------------------------------------------------------
!>
!> @param[in] datein date in the known character string of length 14
!> @param[in] dateform character string containing desired date format
!> @return character string containing desired date format

function formatdate(datein,dateform)
character(len=*), intent(in) :: datein
character(len=*), intent(in) :: dateform

character(len=80) :: formatdate
integer :: ic,iy,im,id,ih,in,is

ic=index(dateform,'cc')
iy=index(dateform,'yy')
im=index(dateform,'mm')
id=index(dateform,'dd')
ih=index(dateform,'hh')
in=index(dateform,'nn')
is=index(dateform,'ss')

formatdate=trim(dateform)

if (ic /= 0) formatdate(ic:ic+1) = datein(1:2)
if (iy /= 0) formatdate(iy:iy+1) = datein(3:4)
if (im /= 0) formatdate(im:im+1) = datein(5:6)
if (id /= 0) formatdate(id:id+1) = datein(7:8)
if (ih /= 0) formatdate(ih:ih+1) = datein(9:10)
if (in /= 0) formatdate(in:in+1) = datein(11:12)
if (is /= 0) formatdate(is:is+1) = datein(13:14)

end function formatdate

!-----------------------------------------------------------------------
!>
!> @param[in] ccyy the century and year
!> @param[in] mm the month
!> @param[in] dd the day
!> @param[in] hh the hour
!> @param[in] nn the minute
!> @param[in] ss the second
!> @return character string of date in CESM format YYYY-MM-DD-SSSSS

function formatCESMdate(ccyy,mm,dd,hh,nn,ss)
integer, intent(in) :: ccyy
integer, intent(in) :: mm
integer, intent(in) :: dd
integer, intent(in) :: hh
integer, intent(in) :: nn
integer, intent(in) :: ss
character(len=80) :: formatCESMdate

integer :: fffff

fffff = hh*3600 + nn*60 + ss
write(formatCESMdate, '(i4.4,1a,2(i2.2,1a),i5.5)') ccyy, '-', mm, '-', dd, '-', fffff

end function formatCESMdate

end program advance_time

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


PROGRAM test_random0
!
!  Purpose:
!    Subroutine test the random number generator random0.
!
!  Record of revisions:
!      Date       Programmer          Description of change
!      ====       ==========          =====================
!    11/22/06    S. J. Chapman        Original code
!
IMPLICIT NONE

! Data dictionary: declare variable types & definitions
REAL :: ave          ! Average of random numbers
INTEGER :: i         ! DO loop index
INTEGER :: iseed     ! Seed for random number sequence
INTEGER :: iseq      ! DO loop index
REAL :: ran          ! A random number
REAL :: sum          ! Sum of random numbers

! Get seed.
WRITE (*,*) 'Enter seed: '
READ (*,*) iseed

! Set seed.
CALL SEED ( iseed )

! Print out 10 random numbers.
WRITE (*,*) '10 random numbers: '
DO i = 1, 10
   CALL random0 ( ran )
   WRITE (*,'(3X,F16.6)') ran
END DO

! Average 5 consecutive 1000-value sequences.
WRITE (*,*) 'Averages of 5 consecutive 1000-sample sequences:'
DO iseq = 1, 5
   sum = 0.
   DO i = 1, 1000
      CALL random0 ( ran )
      sum = sum + ran
   END DO
   ave = sum / 1000.
   WRITE (*,'(3X,F16.6)') ave
END DO

END PROGRAM test_random0


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 02.01.02.02/01                        *
C  *    TEST TITLE : Deletion of a structure network using *
C  *                 a non-existent structure identifier   *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

C PHIGS parameter type for reference handling flag
      INTEGER    PDELE,     PKEEP
      PARAMETER (PDELE = 0, PKEEP = 1)

C PHIGS paramter types

      INTEGER NOSTID, REFFG

      LOGICAL CSSEQ

      CALL INITGL ('02.01.02.02/01')

C open PHIGS

      CALL XPOPPH (ERRFIL, MEMUN)

C creation of a structure network

      CALL POPST (100)
      CALL PEXST (103)
      CALL PEXST (104)
      CALL PCLST
      CALL POPST (103)
      CALL PEXST (105)
      CALL PEXST (105)
      CALL PCLST
      CALL POPST (104)
      CALL PEXST (106)
      CALL PCLST
      CALL POPST (102)
      CALL PEXST (106)
      CALL PEXST (109)
      CALL PCLST
      CALL POPST (106)
      CALL PEXST (107)
      CALL PEXST (108)
      CALL PCLST
      CALL POPST (109)
      CALL PEXST (108)
      CALL PEXST (110)
      CALL PCLST

      CALL SETMSG ('9', '<Delete structure network> with a '          //
     1             'non-existent structure identifier and the '       //
     2             'reference handling flag set to delete should '    //
     3             'have no effect on the contents of the CSS.')

C nostid = non-existent structure identifier

      NOSTID = 25
      REFFG = PDELE
      CALL PDSN (NOSTID, REFFG)
      CALL IFPF (CSSEQ ('100, 1, 103, 1, 105, 0/'          //
     1                  '100, 1, 103, 2, 105, 0/'          //
     2                  '100, 2, 104, 1, 106, 1, 107, 0/'  //
     3                  '100, 2, 104, 1, 106, 2, 108, 0/'  //
     4                  '102, 1, 106, 1, 107, 0/'          //
     5                  '102, 1, 106, 2, 108, 0/'          //
     6                  '102, 2, 109, 1, 108, 0/'          //
     7                  '102, 2, 109, 2, 110, 0/'))

      CALL SETMSG ('9', '<Delete structure network> with a '          //
     1             'non-existent structure identifier and the '       //
     2             'reference handling flag set to keep should have ' //
     3             'no effect on the contents of the CSS.')

      NOSTID = 40
      REFFG = PKEEP
      CALL PDSN (NOSTID, REFFG)
      CALL IFPF (CSSEQ ('100, 1, 103, 1, 105, 0/'          //
     1                  '100, 1, 103, 2, 105, 0/'          //
     2                  '100, 2, 104, 1, 106, 1, 107, 0/'  //
     3                  '100, 2, 104, 1, 106, 2, 108, 0/'  //
     4                  '102, 1, 106, 1, 107, 0/'          //
     5                  '102, 1, 106, 2, 108, 0/'          //
     6                  '102, 2, 109, 1, 108, 0/'          //
     7                  '102, 2, 109, 2, 110, 0/'))

      CALL ENDIT
      END


text:
C    @(#)busqua.f	20.5 7/29/96
      subroutine busqua (nb, pk, qk, array)
      dimension array(*)
C
C     This subroutine determines the bus quantities corresponding to
C     load, generation, shunt, unscheduled, and error according to Pk
C     and Qk.
C
C     Input paramters:
C
C     NB     - external bus number.
C     PK     - computed active generation.
C     QK     - computed reactive generation (includes ficticious
C              generation from adjustable shunts).
C
C     Output paramters:
C
C     ARRAY(1): PGEN   - active generation in MW.
C          (2): QGEN   - reactive generation in MVAR.
C          (3): PLOAD  - active load in MW.
C          (4): QLOAD  - reactive load in MVAR.
C          (5): PGNMAX - maximum active generation in MW.
C          (6): QGNMAX - maximum reactive generation in MVAR.
C          (7): PGNMIN - minimum active generation in MW.
C          (8): QGNMIN - minimum reactive generation in MVAR.
C          (9): SKCOND - total active shunt in MW.
C         (10): TOTCAP - total capacitive shunt in MVAR.
C         (11): (not used -- negative conductance has no significance.)
C         (12): TOTREK - total reactive shunt in MVAR.
C         (13): SKCOND - used active shunt in MW.
C         (14): USECAP - used capacitive shunt in MVAR.
C         (15): (not used -- negative conductance has no significance.)
C         (16): USEREK - used capacitive shunt in MVAR.
C         (17): PVLV - d-c converter losses in MWs.
C         (18): QVLV - d-c converter losses in MVARs.
C         (19): GEQUIV - equivalent shunt impedance in MWs.
C         (20): BEQUIV - equivalent shunt impedance in MVARs.
C         (21): (not used -- there is no unscheduled conductance.)
C         (22): UNSKED - reactive shunt unscheduled in MVAR.
C         (23): PERR   - active shunt mismatched in MW.
C         (24): QERR   - reactive shunt mismatched in MVAR.
C         (25): PLOAD  - Constant MVA load in MWs.
C         (26): QLOAD  - Constant MVA load in MVARs.
C         (27): PLOAD  - Constant current load in MWs.
C         (28): QLOAD  - Constant current load in MVARs.
C         (29): PLOAD  - Constant impedance load in MWs.
C         (30): QLOAD  - Constant impedance load in MVARs.
C

      include 'ipfinc/parametr.inc'

      include 'ipfinc/alpha2.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/branch.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/cbus.inc'
      include 'ipfinc/dc2t.inc'
      include 'ipfinc/dcmt.inc'
      include 'ipfinc/lodtyp.inc'
      include 'ipfinc/ordsta.inc'
      include 'ipfinc/tbx.inc'
      include 'ipfinc/xdata.inc'
C
C     Local Variables
C
      integer j1, j2, pointr
      logical capsw
      character kode2 * 1, kodeyr * 2, kowner * 3

      if (tbx_loaded .ne. ordtbx) then
         do i = 1, ntot
            ptrtbx(i)  = 0
         enddo
         do jt = 1, ntotb
            nb = tbx(2,jt)
            if (ordtbx .eq. 2) nb = opt2inp(nb)
            ptrtbx(nb) = jt
         enddo
         tbx_loaded = ordtbx
      endif
C
C     Initialize composite load quantities
C
      do i = 1, 6
         tottyp(i) = 0.0
      enddo
      gequiv = 0.0
      bequiv = 0.0
      bfixed = 0.0

      ktyp = kbsdta(1,nb)
      if (ktyp .eq. 5) then
C
C        Compute bus quantities for type BD bus.
C
         do 90 i = 1, 26
   90    array(i) = 0.0
         do 92 idc = 1,kdtot
            if (dc2t(1,idc) .eq. nb) then
               pdc=real(dc2t(42,idc))
               qdc=real(dc2t(44,idc))
               vdc=real(dc2t(40,idc))
               pvlv=abs(pdc-0.001*real(dc2t(39,idc))*real(dc2t(40,idc)))
               qvlv=qdc
               adc=real(dc2t(39,idc))
               array(17) = pvlv
               array(18) = qvlv
               go to 94
            else if (dc2t(3,idc) .eq. nb) then
               pdc=real(dc2t(43,idc))
               qdc=real(dc2t(45,idc))
               vdc=real(dc2t(41,idc))
               pvlv=abs(pdc+0.001*real(dc2t(39,idc))*real(dc2t(41,idc)))
               qvlv=qdc
               adc=-real(dc2t(39,idc))
               array(17) = pvlv
               array(18) = qvlv
               go to 94
            else
            endif
   92    continue
         call erexit
   94    continue
         pgen = 0.0
         qgen = 0.0
         pload = 0.0
         qload = 0.0
         pgnmax = 0.0
         qgnmax = 0.0
         pgnmin = 0.0
         qgnmin = 0.0
         skcond = 0.0
         sksusp = 0.0
         totcap = 0.0
         totrek = 0.0

      else if (ktyp .eq. 12) then
C
C        Compute bus quantities for type BM bus.
C
         do 96 idc = 1,mtdcbs
            if (dcmtbs(1,idc) .eq. nb) then
               pdc = real(dcmtbs(25,idc))
               qdc = real(dcmtbs(26,idc))
               vdc = real(dcmtbs(20,idc))
               pvlv = abs(pdc - real(dcmtbs(19,idc)))
               qvlv = qdc
               if (dcmtbs(20,idc) .eq. 0.0) then
                  adc = 0.0
               else
                  adc = 1000.0 * real(dcmtbs(19,idc)) /
     &                  real(dcmtbs(20,idc))
               endif
               array(17) = pvlv
               array(18) = qvlv
               go to 98
            endif
   96    continue
         call erexit
   98    continue
         pgen = 0.0
         qgen = 0.0
         pload = 0.0
         qload = 0.0
         pgnmax = 0.0
         qgnmax = 0.0
         pgnmin = 0.0
         qgnmin = 0.0
         skcond = 0.0
         sksusp = 0.0
         totcap = 0.0
         totrek = 0.0

      else
         array(17) = 0.0
         array(18) = 0.0
         pgnmax = busdta(7,nb)
         pgnmin = 0.0
         if (ktyp .eq. 1 .or. ktyp .eq. 4 .or. ktyp .eq. 10) then
            if (busdta(10,nb) .lt. 0.0) then
               qgnmax = 0.0
               qgnmin = 0.0
            else
               qgnmax = busdta(9,nb)
               qgnmin = busdta(9,nb)
            endif
         else
            qgnmax = busdta(9,nb)
            qgnmin = busdta(10,nb)
         endif
C
C        Compute bus quantities for all other type buses.
C
         kt = inp2opt(nb)
         pgen = pk * bmva
         qgen = qk * bmva
         vk = real( dsqrt(e(kt) ** 2 + f(kt) ** 2) )
         pload = busdta(3,nb)
         qload = busdta(4,nb)
         jtbx = ptrtbx(nb)
C
C        Obtain total susceptance.  Note: BX bus type store initial
C        value on B record.  The total is in array XDATA.
C
         if (ktyp .eq. 11) then
            if (jtbx .gt. 0) then
               nt = tbx(5,jtbx)
               totrek = xdata(3,nt) * vk ** 2
               totcap = xdata(4,nt) * vk ** 2
            else
               call erexit
            endif
         else
            skcond = busdta(5,nb) * vk ** 2
            sksusp = busdta(6,nb) * vk ** 2
            totcap = amax1 (sksusp, 0.0)
            totrek = amin1 (sksusp, 0.0)
C
         endif

      endif

      tottyp(1) = pload
      tottyp(2) = qload

      ncb = kbsdta(15,nb)
      do while (ncb .gt. 0)
         pload2 = bctbl(2,ncb)
         qload2 = bctbl(3,ncb)
         skcon2 = bctbl(4,ncb) * vk ** 2
         sksus2 = bctbl(5,ncb) * vk ** 2
         if (ktyp .eq. 1 .or. ktyp .eq. 4 .or. ktyp .eq. 10) then
            if (bctbl(12,ncb) .lt. 0.0) then
            else
               qgnmax = qgnmax + bctbl(11,ncb)
               qgnmin = qgnmin + bctbl(11,ncb)
            endif
         else
            qgnmax = qgnmax + bctbl(11,ncb)
            qgnmin = qgnmin + bctbl(12,ncb)
         endif
         call getchr(1,kode2,kbctbl(8,ncb))
         call getchr(2,kodeyr,kbctbl(9,ncb))
         call getchr(3,kowner,kbctbl(10,ncb))
C
C        Convert constant current and constant admittance
C        loads back to constant MVA.
C
         do 220 i = 1, 6
  220    lodtyp(i) = 0.0
         if (kode2 .eq. 'A') then
            if (kodeyr .eq. '01' .or. kodeyr .eq. '*I') then
               pload2 = pload2 * vk
               qload2 = qload2 * vk
               lodtyp(3) = pload2
               lodtyp(4) = qload2
               if (kowner .eq. '***') then
                  gequiv = gequiv + skcon2
                  bequiv = bequiv + sksus2
                  skcon2 = 0.0
                  sksus2 = 0.0
               else
                  lodtyp(5) = skcon2
                  lodtyp(6) = -sksus2
                  pload2 = pload2 + skcon2
                  qload2 = qload2 - sksus2
                  skcon2 = 0.0
                  sksus2 = 0.0
               endif
            else if (kodeyr .eq. '02' .or. kodeyr .eq. '*P') then
               lodtyp(1) = pload2
               lodtyp(2) = qload2
               if (kowner .eq. '***') then
                  gequiv = gequiv + skcon2
                  bequiv = bequiv + sksus2
                  skcon2 = 0.0
                  sksus2 = 0.0
               else
                  lodtyp(5) = skcon2
                  lodtyp(6) = -sksus2
                  pload2 = pload2 + skcon2
                  qload2 = qload2 - sksus2
                  skcon2 = 0.0
                  sksus2 = 0.0
               endif
            else
               lodtyp(1) = pload2
               lodtyp(2) = qload2
               bfixed = bfixed + sksus2
               sksus2 = 0.0
            endif
         else if (kodeyr .eq. '*I') then
            pload2 = pload2*vk
            qload2 = qload2*vk
            lodtyp(3) = pload2
            lodtyp(4) = qload2
            if (kowner .eq. '***') then
               gequiv = gequiv + skcon2
               bequiv = bequiv + sksus2
               skcon2 = 0.0
               sksus2 = 0.0
            else
               lodtyp(5) = skcon2
               lodtyp(6) = -sksus2
               pload2 = pload2 + skcon2
               qload2 = qload2 - sksus2
               skcon2 = 0.0
               sksus2 = 0.0
            endif
         else if (kodeyr .eq. '*P') then
            lodtyp(1) = pload2
            lodtyp(2) = qload2
            if (kowner .eq. '***') then
               gequiv = gequiv + skcon2
               bequiv = bequiv + sksus2
               skcon2 = 0.0
               sksus2 = 0.0
            else
               lodtyp(5) = skcon2
               lodtyp(6) = -sksus2
               pload2 = pload2 + skcon2
               qload2 = qload2 - sksus2
               skcon2 = 0.0
               sksus2 = 0.0
            endif
         else
            lodtyp(1) = pload2
            lodtyp(2) = qload2
         endif
         pload = pload + pload2
         qload = qload + qload2
         skcond = skcond + skcon2
         sksusp = sksusp + sksus2
         totcap = totcap + amax1 (sksus2, 0.0)
         totrek = totrek + amin1 (sksus2, 0.0)
         do 99 i = 1, 6
   99    tottyp(i) = tottyp(i) + lodtyp(i)
         ncb = bctbl_nxt(ncb)
      enddo
      totcap = totcap + amax1( bfixed, 0.0 )
      totrek = totrek + amin1( bfixed, 0.0 )
C
C     Remove switched shunt from BQ busses.
C
C     "CAPSW" = .TRUE./.FALSE.
C                  (shunt susceptance is fixed/adjustable)
C
      usecap = totcap
      userek = totrek
      capsw = .true.
      if (ktyp .eq. 1 .or. ktyp .eq. 4 .or. ktyp .eq. 10) then
      else if (usecap .gt. 0.0 .or. userek .lt. 0.0) then

         if (ktyp .eq. 7) then

            if (jtbx .gt. 0) then
               ltyp = tbx(1,jtbx)
               ityp = tbx(7,jtbx)
               if (ityp .eq. 3) then
C
C                 Q_min limit hit. Remove all capacitors.
C
                  usecap = amax1( 0.0, bfixed )
C
               else if (ityp .eq. 4) then
C
C                 Q_max limit hit. Remove all reactors.
C
                  userek = amin1( 0.0, bfixed )
C
               else
                  capsw = .false.
               endif
            else
               call erexit
            endif
C
C           "BG" shunt is held constant. CAPSW does not need to be
C           set to .TRUE., which implies fixed Qshunt, since .TRUE.
C           is its default value.
C
         else if (ktyp .eq. 11) then
C
C           "BX" shunt is explicitly defined.
C
            if (jtbx .gt. 0) then
               nt = tbx(5,jtbx)
               userek = xdata(5,nt) * vk ** 2 + amin1 (0.0, bfixed)
               usecap = xdata(6,nt) * vk ** 2 + amax1 (0.0, bfixed)
            else
               call erexit
            endif

         else if (ktyp .eq. 2 .or. ktyp .eq. 3) then
C
C           "BE" and "BS" have allocatable shunt also.
C
            capsw = .false.
         endif
      endif
C
C     Remove "ficticious generation" from BQ, BE, and BS busses.
C     This ficticious generation compensates any unused shunt.
C
      if (.not. capsw) then
         if (qgen * usecap .lt. 0.0) then
            if (abs(qgen) .le. abs(usecap) + 0.5) then
               usecap = usecap + qgen
               qgen = 0.0
            else
               qgen = qgen + usecap
               usecap = 0.0
            endif
         else if (qgen * userek .lt. 0.0) then
            if ( abs( qgen ) .le. abs( userek ) + 0.5 ) then
               userek = userek + qgen
               qgen = 0.0
            else
               qgen = qgen + userek
               userek = 0.0
            endif
         endif
      endif

 1280 unsked = dim (qgen, qgnmax) - dim(qgnmin, qgen)
      qgen = qgen - unsked
      if (abs(unsked) .lt. 0.5) unsked = 0.0

      perr = 0.0
      qerr = 0.0

      array(1) = pgen
      array(2) = qgen
      array(3) = pload
      array(4) = qload
      array(5) = pgnmax
      array(6) = qgnmax
      array(7) = pgnmin
      array(8) = qgnmin
      array(9) = skcond
      array(10) = totcap
      array(11) = 0.0
      array(12) = totrek
      array(13) = skcond
      array(14) = usecap
      array(15) = 0.0
      array(16) = userek
      array(19) = gequiv
      array(20) = bequiv
      array(21) = 0.0
      array(22) = unsked
      array(23) = perr
      array(24) = qerr
      array(25) = tottyp(1)
      array(26) = tottyp(2)
      array(27) = tottyp(3)
      array(28) = tottyp(4)
      array(29) = tottyp(5)
      array(30) = tottyp(6)
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2020 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine distattach_3d(xig,etg,zeg,pneigh,pnode,a,p,ratio,
     &          nterms)
!
!     calculates the distance between the node with coordinates
!     in "pnode" and the node with local coordinates xig and etg
!     in a face described by "nterms" nodes with coordinates
!     in pneigh
!
      implicit none
!
      integer nterms,i,j,n
!
      real*8 ratio(20),pneigh(3,*),pnode(3),a,xi,et,xig,etg,p(3),
     &  dummy,ze,zeg,omg,omh,omr,opg,oph,opr,dx(3),al,
     &  tpgphpr,tmgphpr,tmgmhpr,tpgmhpr,tpgphmr,tmgphmr,tmgmhmr,tpgmhmr,
     &  omgopg,omhoph,omropr,omgmopg,omhmoph,omrmopr
!
!
!
      n=3
!
      if(nterms.eq.4) then
         xi=(xig+1.d0)/2.d0
         et=(etg+1.d0)/2.d0
         ze=(zeg+1.d0)/2.d0
         dx(1)=xi
         dx(2)=et
         dx(3)=ze
         call insertsortd(dx,n)
c         call dsort(dx,iy,n,kflag)
         if(dx(3).gt.1.d-30) then
            al=dx(3)/(xi+et+ze)
            xi=al*xi
            et=al*et
            ze=al*ze
         endif
!
!        shape functions
!
         ratio(1)=1.d0-xi-et-ze
         ratio(2)=xi
         ratio(3)=et
         ratio(4)=ze
      elseif(nterms.eq.6) then
         xi=(xig+1.d0)/2.d0
         et=(etg+1.d0)/2.d0
         if(xi+et.gt.1.d0) then
            dummy=xi
            xi=1.d0-et
            et=1.d0-dummy
         endif
!
         ze=zeg

         a=1.d0-xi-et
!
!     shape functions
!
         ratio(1)=0.5d0*a *(1.d0-ze)
         ratio(2)=0.5d0*xi*(1.d0-ze)
         ratio(3)=0.5d0*et*(1.d0-ze)
         ratio(4)=0.5d0*a *(1.d0+ze)
         ratio(5)=0.5d0*xi*(1.d0+ze)
         ratio(6)=0.5d0*et*(1.d0+ze)
      elseif(nterms.eq.8) then
         xi=xig
         et=etg
         ze=zeg
!
         omg=1.d0-xi
         omh=1.d0-et
         omr=1.d0-ze
         opg=1.d0+xi
         oph=1.d0+et
         opr=1.d0+ze
!
!        shape functions
!
         ratio(1)=omg*omh*omr/8.d0
         ratio(2)=opg*omh*omr/8.d0
         ratio(3)=opg*oph*omr/8.d0
         ratio(4)=omg*oph*omr/8.d0
         ratio(5)=omg*omh*opr/8.d0
         ratio(6)=opg*omh*opr/8.d0
         ratio(7)=opg*oph*opr/8.d0
         ratio(8)=omg*oph*opr/8.d0
      elseif(nterms.eq.10) then
         xi=(xig+1.d0)/2.d0
         et=(etg+1.d0)/2.d0
         ze=(zeg+1.d0)/2.d0
         dx(1)=xi
         dx(2)=et
         dx(3)=ze
         call insertsortd(dx,n)
c         call dsort(dx,iy,n,kflag)
         if(dx(3).gt.1.d-30) then
            al=dx(3)/(xi+et+ze)
            xi=al*xi
            et=al*et
            ze=al*ze
         endif
c         if(xi+et+ze.gt.1.d0) then
c            dummy=2.d0*(1.d0-xi-et-ze)/3.d0
c            xi=dummy+xi
c            et=dummy+et
c            ze=dummy+ze
c         endif
!
!        shape functions
!
         a=1.d0-xi-et-ze
         ratio( 1)=(2.d0*a-1.d0)*a
         ratio( 2)=xi*(2.d0*xi-1.d0)
         ratio( 3)=et*(2.d0*et-1.d0)
         ratio( 4)=ze*(2.d0*ze-1.d0)
         ratio( 5)=4.d0*xi*a
         ratio( 6)=4.d0*xi*et
         ratio( 7)=4.d0*et*a
         ratio( 8)=4.d0*ze*a
         ratio( 9)=4.d0*xi*ze
         ratio(10)=4.d0*et*ze
      elseif(nterms.eq.15) then
         xi=(xig+1.d0)/2.d0
         et=(etg+1.d0)/2.d0
         if(xi+et.gt.1.d0) then
            dummy=xi
            xi=1.d0-et
            et=1.d0-dummy
         endif
!
         ze=zeg
         a=1.d0-xi-et
!
!     shape functions
!
         ratio(1)=-0.5d0*a*(1.d0-ze)*(2.d0*xi+2.d0*et+ze)
         ratio(2)=0.5d0*xi*(1.d0-ze)*(2.d0*xi-2.d0-ze)
         ratio(3)=0.5d0*et*(1.d0-ze)*(2.d0*et-2.d0-ze)
         ratio(4)=-0.5d0*a*(1.d0+ze)*(2.d0*xi+2.d0*et-ze)
         ratio(5)=0.5d0*xi*(1.d0+ze)*(2.d0*xi-2.d0+ze)
         ratio(6)=0.5d0*et*(1.d0+ze)*(2.d0*et-2.d0+ze)
         ratio(7)=2.d0*xi*a*(1.d0-ze)
         ratio(8)=2.d0*xi*et*(1.d0-ze)
         ratio(9)=2.d0*et*a*(1.d0-ze)
         ratio(10)=2.d0*xi*a*(1.d0+ze)
         ratio(11)=2.d0*xi*et*(1.d0+ze)
         ratio(12)=2.d0*et*a*(1.d0+ze)
         ratio(13)= a*(1.d0-ze*ze)
         ratio(14)=xi*(1.d0-ze*ze)
         ratio(15)=et*(1.d0-ze*ze)
      elseif(nterms.eq.20) then
         xi=xig
         et=etg
         ze=zeg
!
         omg=1.d0-xi
         omh=1.d0-et
         omr=1.d0-ze
         opg=1.d0+xi
         oph=1.d0+et
         opr=1.d0+ze
         tpgphpr=opg+oph+ze
         tmgphpr=omg+oph+ze
         tmgmhpr=omg+omh+ze
         tpgmhpr=opg+omh+ze
         tpgphmr=opg+oph-ze
         tmgphmr=omg+oph-ze
         tmgmhmr=omg+omh-ze
         tpgmhmr=opg+omh-ze
         omgopg=omg*opg/4.d0
         omhoph=omh*oph/4.d0
         omropr=omr*opr/4.d0
         omgmopg=(omg-opg)/4.d0
         omhmoph=(omh-oph)/4.d0
         omrmopr=(omr-opr)/4.d0
!
!     shape functions
!
         ratio( 1)=-omg*omh*omr*tpgphpr/8.d0
         ratio( 2)=-opg*omh*omr*tmgphpr/8.d0
         ratio( 3)=-opg*oph*omr*tmgmhpr/8.d0
         ratio( 4)=-omg*oph*omr*tpgmhpr/8.d0
         ratio( 5)=-omg*omh*opr*tpgphmr/8.d0
         ratio( 6)=-opg*omh*opr*tmgphmr/8.d0
         ratio( 7)=-opg*oph*opr*tmgmhmr/8.d0
         ratio( 8)=-omg*oph*opr*tpgmhmr/8.d0
         ratio( 9)=omgopg*omh*omr
         ratio(10)=omhoph*opg*omr
         ratio(11)=omgopg*oph*omr
         ratio(12)=omhoph*omg*omr
         ratio(13)=omgopg*omh*opr
         ratio(14)=omhoph*opg*opr
         ratio(15)=omgopg*oph*opr
         ratio(16)=omhoph*omg*opr
         ratio(17)=omropr*omg*omh
         ratio(18)=omropr*opg*omh
         ratio(19)=omropr*opg*oph
         ratio(20)=omropr*omg*oph
      else
         write(*,*) '*ERROR in distattach: case with ',nterms
         write(*,*) '       terms is not covered'
         call exit(201)
      endif
!
!     calculating the position in the face
!
      do i=1,3
         p(i)=0.d0
         do j=1,nterms
            p(i)=p(i)+ratio(j)*pneigh(i,j)
         enddo
      enddo
!
!     calculating the distance
!
      a=(pnode(1)-p(1))**2+(pnode(2)-p(2))**2+(pnode(3)-p(3))**2
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
Davis Motorsports is Davis only motorcycle shop, offering repairs, parts, and new/used scooter sales.

They are now under new management (former shop owner Jeff retired.)

  To learn more about motorcycles and scooters in Davis, visit our Motorcycles page.

20051202 08:06:38 nbsp Ive had discussions and some minor service done through Davis Motorsports and have nothing but good things to say thus far! Seems to be owned by the same man who does everything else in the place and hes very easy to speak with. The inventory may not look very big, but my guess is that its not worth maintain a large inventory for a bike shop, especially when he stands ready and willing to place orders for customers. Definitely would recommend others give this shop a try before heading out to Woodland, Fairfield, or Sacramento. Users/AlphaDog

20051210 14:10:45 nbsp Jeff, the head mechanic is awesome. Like AlphaDog said, they dont have very much in stock, but they can order anything you need. I went in to have my bike looked at, small oil leak, he had it for a day or two, and took a look and siad it wasnt a major issue and didnt charge me for the service. Users/JeffNakamura

20060406 13:48:27 nbsp Great shop, he does a good job on repairs and at least for me did the repairs in a cost effective way that avoided me having to pay for lots of expensive new parts. Users/AaronBurda

20070227 11:25:15 nbsp honest and fairly priced.best bike service i have ever had in 20 years........thanks jeff Users/RichardRusakowicz

20070711 16:38:33 nbsp They dont like gopeds, neither does JBs Users/StevenDaubert

20071005 09:57:08 nbsp Ive used Jeff as my motorcycle mechanic since 1990. (OMG 17 years!) He does it right the first time at a reasonable price. Hes by far the best mechanic I have worked with. Also, hell discuss your options and figure out the lowest cost way to get you back on the road. Finally, he sells one of the best scooter brands available! Users/TimDunn

20071017 21:12:47 nbsp Jeff worked on my bike today and I was impressed with his work, attitude, and price.  Top notch mechanic! Users/bordiga

20080308 21:32:13 nbsp Great service!!! honest and fair. All service departments ( motorcycle ,bike,car) should take a lesson from Jeff.. Users/zonepunk

20080619 17:41:10 nbsp Ive been there today for my first time... Im very impressed, very good work, nice people there, they gave me some nice suggestion. Very happy for their service! Users/Fede

20080820 14:36:16 nbsp I have been bringing my Vespa there for three years and have nothing but good things to say. Their work is always solid and their service is very personal. They know their stuff. Users/SteveKent

20100514 11:45:55 nbsp Jeff is a very good mechanic, he taught me several things about my bike and helped maintain it in good running condition. Users/LinaPlanuty

20110221 02:17:01 nbsp Ive gone to Davis Motorsports many times and Ive always been very impressed.The owner really listens to you and is always up front and honest as to what he think is wrong and what its going to cost. Hes fixed repair jobs that other places, namely Berkeley Performance Motorcycles, messed up. This place does an honest job at an honest price, Ive never been unsatisfied with the labor, estimates and turn around time. If you ride in Davis, or even Sacramento or Woodland this should be your place of business. After having brought my bike to this shop a handful of times the owner remembers me, and always treats me with respect, which as a young guy I dont often find from motorcycle shop mechanics. A ++ Users/ReneMacleay

20110420 09:30:55 nbsp Like everyone here pretty much says Jeff is a great guy that will be straight up with you and do an honest and fantastic job. If Im doing work on my bike and run into problems, I can call him up and he is willing to help out. Hes very knowledgeable when it comes to bikes, and what I assume quads and scooters as well. I have never felt that this place was unreasonable and the service is fast! I definitely recommend that you come here first when running into trouble. There really is no other competition when it comes to Davis Motorsports.

Keep the rubber side down. Users/DanielLebe

20110802 10:41:42 nbsp If the owner is reading this, they should set up and list a company email. I am abroad and cannot answer the phone but he just called me (he has my bike at the moment)
 Users/MattCorey

20111121 10:22:27 nbsp Jeff is a great guy who does honest work! He has tried to help me get the biggest bang for my buck with my scooter by keeping all repairs and maintenance directed at those most important while letting me know about the repairs I will need to address at a later time. We are fortunate to have such a trustworthy tradesman and shop owner in Davis!  Users/THOMASSLABAUGHII

20120419 10:50:35 nbsp This happened last year, but i was riding through downtown and my r6 would sporadically lose power and die. I kept starting it back up and was able to make it to the shop. I explained what happened to Jeff, and he simply tightened up my battery plugs and it was ready to go. No Charge or anything! Ive brought in my bike for service before too, and Jeff provides excellent customer service! I highly recommend his shop. Users/FiddyFetzer

20120619 15:26:16 nbsp Jeff has been gracious in providing after hours consultation with a prospective buyer of my sons motorcycle the problems with which Jeff diagnosed a year ago.  I was very pleased with his manner and his helpfulness.  Very easy going and capable.  I understand why he is so highly regardeda definite plus for Davis.  Thanks Jeff for your good heart.  Mikos Fabersunne Users/MikosFabersune

20130608 13:05:57 nbsp Jeff retired and sold the business to experienced mechanics. They opened on Memorial Day weekend, 2013. I am giving them a try today and hope they can help me as much as Jeff has done so in the past! Users/THOMASSLABAUGHII

20130610 11:51:55 nbsp Is Davis Motorsports still open? I keep calling but nobody answers and the place was closed when came by. Are there any places in Davis that can change a motorcycle tire?  Users/TomasSudnius

20130730 16:30:51 nbsp Nice people, but their incompetence was shown. Took my bike here for an oil leak problem 2 weeks ago. They keep making excuses as to why I havent gotten my bike back. Its just a startup shop good for quick and easy maintenance. In hindsight, I shouldve just shelled out the money, took it to the dealer, and payed the higher fee to get it back within a day. Ive had valve clearance inspections done in less time than this. Not happy, wont be returning.  Users/StevenTaste

20140811 10:10:53 nbsp Went there to buy and have installed new tires for my Yamaha 1300 Tourer. Talked to the new owner, ordered tires and at the same time, feeling very uneasy. After going home and really thinking about the conversation with this guy I decided not to have him do this work, phone called him and told him so.. He kept the $100 deposit which he had to do because of the expense of ordering the tires. Yes, it had been less than 24 hours. I wont be going back which is too bad as Jeff really took care of me when he had the place.. After seeing and hearing of others experiences, I am glad I got out of there...who knows, that 100 dollars may have saved my life,,, Users/frankdupree


text:
      subroutine gensdp
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
C***********************************************************************
CD1
CD1  PURPOSE
CD1
CD1  This subroutine solves the isothermal air-water equations with
CD1  full jacobian (unsymmetric, 2n by 2n).
CD1
C***********************************************************************
CD2
CD2  REVISION HISTORY
CD2
CD2 $Log:   /pvcs.config/fehm90/src/gensdp.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:43:06   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:06:32   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:09:24   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:24:08   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 12:02:30   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:41:42 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2    Rev 1.16   Mon Jul 01 14:17:22 1996   hend
CD2 Changed to use normal_dof
CD2
CD2    Rev 1.15   Fri May 31 15:10:44 1996   gaz
CD2 cleaned up some comments
CD2
CD2    Rev 1.14   Tue May 21 13:50:52 1996   gaz
CD2 added call to md_nodes
CD2
CD2    Rev 1.13   Thu Feb 15 10:48:24 1996   zvd
CD2 Added requirement.
CD2
CD2    Rev 1.12   Wed Jan 17 08:56:42 1996   hend
CD2 Removed nonconstant parameter statement for IBM
CD2
CD2    Rev 1.11   Tue Jan 16 14:31:04 1996   hend
CD2 Added capability for 5,6, and n degrees of freedom
CD2
CD2    Rev 1.10   Fri Jan 12 17:48:40 1996   llt
CD2 changed mmgetblk agruments
CD2
CD2    Rev 1.9   Wed Jan 10 12:27:52 1996   hend
CD2 Corrected ECD Number
CD2
CD2    Rev 1.8   Wed Jan 10 11:47:44 1996   hend
CD2 Corrected ECD Number
CD2
CD2    Rev 1.7   Wed Jan 10 10:28:18 1996   hend
CD2 Added Prolog
CD2
CD2    Rev 1.6   11/15/95 10:29:58   gaz
CD2 changes to get rdof working properly
CD2
CD2    Rev 1.5   03/23/95 23:55:40   gaz
CD2 gaz put abs(bp instead of bp ...
CD2
CD2    Rev 1.4   05/11/94 16:07:10   llt
CD2 bug fixes - gaz
CD2
CD2    Rev 1.3   03/24/94 12:53:12   llt
CD2 Changed mmgetblk for array b at bottom to mmrelblk. Moved
CD2 deallocation of array sto1 to after continue, so can't be missed.
CD2
CD2    Rev 1.2   03/23/94 14:41:06   robinson
CD2 Additional cleanup of memory management
CD2
CD2    Rev 1.1   03/18/94 15:47:30   gaz
CD2 Added solve_new and cleaned up memory management.
CD2
CD2    Rev 1.0   01/20/94 10:24:18   pvcs
CD2 original version in process of being certified
CD2
CD2    c 17-mar-94 gaz
CD2    c getting and releasing b matrix
CD2
C***********************************************************************
CD3
CD3  REQUIREMENTS TRACEABILITY
CD3
CD3  2.3.2 Heat- and mass-transfer equations
CD3  2.3.3 Noncondensible gas flow equations
CD3  2.5.2 Solve nonlinear equation set at each time step
CD3
C***********************************************************************
CD4
CD4  SPECIAL COMMENTS AND REFERENCES
CD4
CD4 Requirements from SDN: 10086-RD-2.20-00
CD4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
CD4   FEHM Application Version 2.20
CD4
C***********************************************************************

      use davidi
      use comhi
      use comgi
      use comfi
      use comei
      use comdi
      use comci
      use combi
      use comdti
      use comai
      implicit none

      integer nc(16)
      integer nmatd,ndex(6)
      integer i, id, idl, neqp1, nrhs1, nrhs2, nsizea
      real*8, allocatable :: dumn(:)
      real*8, allocatable :: sto1(:)
      real*8, allocatable :: sto5(:)
      real*8, allocatable :: dum(:)
      real*8  facr, fdm, fdum2, tollr, tolls
      parameter(fdm=20.0)

c call thermo for next level (i+neq)
c
      call dpdp(2)
c
c     adjust storage
c
      neqp1=neq+1
      fdum2=0.
      nmatd=nelm(neqp1)-neqp1
      nsizea=16*nmatd
      do i=1,nsizea
         a(i)=0.0
      enddo
      do i=1,neq
         bp(i+nrhs(1))=0.0
         bp(i+nrhs(2))=0.0
         bp(i+nrhs(3))=0.0
         bp(i+nrhs(4))=0.0
      enddo
      nrhs1=nrhs(1)
      nrhs2=nrhs(2)

      do id=1,neq
         nmat(1)=0
         nmat(2)=nmatd
         nmat(3)=nmatd*(5-1)
         nmat(4)=nmatd*(6-1)
         nrhs(1)=nrhs1
         nrhs(2)=nrhs2
c
c     decide on equation type
c
         if(ico2.lt.0) then
            if(ieos(id).ne.4) then
               call geneq2( id)
               call add_accumulation(id)
            else
               call geneq3( id)
               a(nelmdg(id)+nmat(1))=sx1(id)
               bp(id+nrhs(1))=0.0
            endif
         else if (ico2.eq.0) then
            if(ieos(id).ne.4) then
               call geneq1( id)
            else
               call geneq3( id)
               a(nelmdg(id)+nmat(1))=sx1(id)
               bp(id+nrhs(1))=0.0
            endif
         endif
c
c     decide on equation type(for node idl=id+neq)
c
         idl=id+neq
         nmat(1)=nmatd*(11-1)
         nmat(2)=nmatd*(12-1)
         nmat(3)=nmatd*(15-1)
         nmat(4)=nmatd*(16-1)
         nrhs(1)=nrhs(3)
         nrhs(2)=nrhs(4)
         if(ico2.lt.0) then
            if(ieos(idl).ne.4) then
               call geneq2(idl)
               call add_accumulation(idl)
            else
               call geneq3(idl)
               a(nelmdg(id)+nmat(1))=sx1(id)
               bp(id+nrhs(1))=0.0
            endif
         else if (ico2.eq.0) then
            if(ieos(idl).ne.4) then
               call geneq1(idl)
            else
               call geneq3(idl)
               a(nelmdg(id)+nmat(1))=sx1(id)
               bp(id+nrhs(1))=0.0
            endif
         endif
      enddo
      do i=1,16
         nmat(i)=nmatd*(i-1)
      enddo
      nrhs(1)=nrhs1
      nrhs(2)=nrhs2
c
c     call md_modes to complete equations for multiply defined nodes
c
c      if(imdnode.ne.0) call md_nodes(3,0,0)
c
c     now get coupling for influence terms
c
      if(ico2.eq.0) then
         call dpdpfh
      else if(ico2.lt.0) then
         call dpdpfa
      endif

c     normalize the equations
c
c     some changes to call sequence here
c
      if(islord.ne.0) then
         call switch(nmat,nmatb,islord,4,1,nmatd)
         call switchb(bp,nrhs,nrhsb,islord,4,1,neq)
      end if

      allocate(dumn(36))
c     call nrmlz4(neq,a,nmat,bp,nrhs,nelm,fdum2,sto1)
      call normal_dof(neq,a,bp,nelm,nmat,nrhs,nelmdg
     &    ,ndex,4,dumn(1),dumn(17),dumn(33),0,fdum2)

      deallocate(dumn)

c
c     return and stop if singular
c
      if(ndex(1).lt.0 ) then
         if (iout .ne. 0) then
            write(iout,*) '   '
            write(iout, 100)
            write(iout,*) '   '
            write(iout,'(i8,3g12.3)')
     &           abs(nmatb(1)),(cord(abs(nmatb(1)),i),i=1,3)
         end if
         if(iptty.gt.0) then
            write(iptty,*) '   '
            write(iptty,100)
            write(iptty,*) '   '
            write(iptty,'(i8,3g12.3)')
     &           abs(nmatb(1)),(cord(abs(nmatb(1)),i),i=1,3)
         endif
         iad=maxit
         return
      endif
 100  format ('* singular matrix found during normalization *')

      fdum=sqrt(fdum2)
      if(fdum.eq.0.0) go to 999
      mink=n
      if(iad.eq.0) then
         f0=max(fdum*epe,tmch)
      endif
      if(fdum1.lt.0.0.and.iad.ne.0) then
         do i=1,neq
            if(abs(bp(i+nrhs(1))).gt.tmch) go to 99
            if(abs(bp(i+nrhs(2))).gt.tmch) go to 99
            if(abs(bp(i+nrhs(3))).gt.tmch) go to 99
            if(abs(bp(i+nrhs(4))).gt.tmch) go to 99
         enddo
         fdum=-1.0
         go to 999
 99      continue
         f0=-1.0
      endif
      if(fdum.le.f0.and.iad.ne.0) goto 999
      facr=1.0
      tolls=max(facr*f0,min(g1*fdum,g2*fdum2))
      tollr=tolls*g3
      if(g3*tmch.gt.tollr) tollr=g3*tmch

c     set maximum iterations in solve_new
c
      iter=maxsolve
      call storage_derivatives(1,1)

      allocate(dum(neq*4*4))
      if(irdof.eq.0) then
c
c     full solution
c
         if (igauss .gt. 1) then
            call solve_new(neq,a,b,bp,nmat,nb,nrhs,nelm,nop,north
     *           ,tollr,irb,iirb,npvt,gmres,dum,piv
     *           ,h,c,ss,g,y,iter,iback,4,iptty,maxor,accm)
         else
            call solve_new(neq,a,b,bp,nmat,nmat,nrhs,nelm,nelm,north
     *           ,tollr,irb,iirb,npvt,gmres,dum,piv
     *           ,h,c,ss,g,y,iter,iback,4,iptty,maxor,accm)
         end if
      else
         allocate(sto5(neq*4))
         if (igauss .gt. 1) then
            call  rdof_new(neq,a,b,bp,nmat,nb,nrhs,nelm,nop,north
     *           ,tollr,irb,iirb,npvt,gmres,dum,piv
     *           ,h,c,ss,g,y,iter,iback,4,iptty,maxor
     *           ,overf,irdof,icoupl,0,sto5,accm)
         else
            call  rdof_new(neq,a,b,bp,nmat,nmat,nrhs,nelm,nelm,north
     *           ,tollr,irb,iirb,npvt,gmres,dum,piv
     *           ,h,c,ss,g,y,iter,iback,4,iptty,maxor
     *           ,overf,irdof,icoupl,0,sto5,accm)
         end if
         deallocate(sto5)
      endif

      deallocate(dum)
      call storage_derivatives(0,1)
      itert=itert+iter
      itotals=itotals+iter
      minkt=minkt+mink
 999  continue
      if(islord.ne.0) then
         call switch(nmat,nmatb,islord,4,2,nmatd)
         call switchb(bp,nrhs,nrhsb,islord,4,2,neq)
      end if

c     extract dpdp solution
c
      call dpdp(4)
c
c     times=times-second(0.0)
c     write(59,*) 'time to solve jacobian(full)', times

      return
      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
Davis Municipal Code Davis Municipal Code/1 Chapter 1
 GENERAL PROVISIONS
 1.01.150 Notice to appear in court for violations.

If any person is arrested for the violation of any provision of this Code or any ordinance of the city and such person is not immediately taken before a magistrate as prescribed in the Penal Code of the state, the arresting officer shall prepare in duplicate a written notice to appear in court, containing the name and address of such person, the offense charged, and the time and place where and when such person shall appear in court. (Code 1964, § 1 5.01.)