text:
!==============================================================================!
  subroutine Source_Kin_K_Eps(grid)
!------------------------------------------------------------------------------!
!   Computes the source terms in kin transport equation for k-epsilon model    !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use all_mod
  use pro_mod
  use les_mod
  use rans_mod
  use Grid_Mod
  use Parameters_Mod
  use Work_Mod, only: kin_x => r_cell_01,  &
                      kin_y => r_cell_02,  &
                      kin_z => r_cell_03
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type) :: grid
!-----------------------------------[Locals]-----------------------------------!
  integer :: c, c1, c2, s
  real    :: UtotSq, Unor, UnorSq, Utan
!==============================================================================!

  !-----------------------------------------------!
  !  Compute the sources in the near wall cells   !
  !-----------------------------------------------!
  if(MODE == HIGH_RE) then
    do s = 1, grid % n_faces
      c1=grid % faces_c(1,s)
      c2=grid % faces_c(2,s)

      if(c2 < 0 .and. TypeBC(c2) /= BUFFER) then
        if(TypeBC(c2)==WALL .or. TypeBC(c2)==WALLFL) then

          ! Compute tangential velocity component
          UtotSq = U % n(c1) * U % n(c1) &
                 + V % n(c1) * V % n(c1) &
                 + W % n(c1) * W % n(c1)
          Unor = ( U % n(c1) * grid % sx(s)     &
                 + V % n(c1) * grid % sy(s)     &
                 + W % n(c1) * grid % sz(s) )   &
                 / sqrt(  grid % sx(s)*grid % sx(s)    &
                        + grid % sy(s)*grid % sy(s)    &
                        + grid % sz(s)*grid % sz(s))
          UnorSq = Unor*Unor

          if(UtotSq  >  UnorSq) then
            Utan = sqrt(UtotSq - UnorSq)
          else
            Utan = TINY
          end if

          ! Compute nondimensional wall distance and wall-shear stress
          if(ROUGH == NO) then
            Ynd(c1) = sqrt(kin % n(c1)) * Cmu25 * WallDs(c1)  &
                    / VISc
            TauWall(c1) = abs(DENc(material(c1))                   &
                        * kappa * sqrt(kin % n(c1)) * Cmu25 * Utan &
                        / (log(Elog*Ynd(c1))))

            ! Compute production in the first wall cell
            Pk(c1) = TauWall(c1) * Cmu25 * sqrt(kin % n(c1)) &
                   / (kappa*WallDs(c1))

          else if(ROUGH==YES) then
            Ynd(c1) = sqrt(kin % n(c1)) * Cmu25 * (WallDs(c1)+Zo)  &
                    / VISc
            TauWall(c1) = abs(DENc(material(c1))                     &
                        * kappa * sqrt(kin % n(c1)) * Cmu25 * Utan   &
                        / (log((WallDs(c1)+Zo)/Zo)))

            Pk(c1) = TauWall(c1) * Cmu25 * sqrt(kin % n(c1)) &
                   / (kappa*(WallDs(c1)+Zo))
            kin % n(c2) = TauWall(c1)/0.09**0.5
          end if

          ! Filling up the source term
          b(c1) = b(c1) + Pk(c1) * grid % vol(c1)
        end if  ! TypeBC(c2)==WALL or WALLFL
      end if    ! c2 < 0
    end do

    !-----------------------------------------!
    !   Compute the sources in the interior   !
    !-----------------------------------------!
    do c=1,grid % n_cells

      ! IsNearWall ensures not to put Pk twice into the near wall cells
      if(.NOT. IsNearWall(c)) then

        ! Production:
        Pk(c)= VISt(c) * Shear(c)*Shear(c)
        b(c) = b(c) + Pk(c) * grid % vol(c)
      end if

      ! Dissipation:
      A % val(A % dia(c)) = A % val(A % dia(c))                          &
                          + DENc(material(c)) * eps % n(c) / kin % n(c)  &
                          * grid % vol(c)
    end do
  end if    ! end if mode = wf

  !--------------------------------------------------------!
  !   Jones-Launder model and Launder-Sharma + Yap model   !
  !--------------------------------------------------------!
  if(MODE == LOW_RE) then
    do c = 1, grid % n_cells

      ! Production:
      Pk(c)= VISt(c) * Shear(c) * Shear(c)
      b(c) = b(c) + Pk(c) * grid % vol(c)

      ! Dissipation:
      A % val(A % dia(c)) = A % val(A % dia(c)) + &
      DENc(material(c))*eps%n(c)/(kin%n(c)+TINY)*grid % vol(c)

      ! Buoyancy contribution:
      if (BUOY == YES) &
      b(c) = b(c) + Pbuoy(c) * grid % vol(c)

      ! Preparation of kin for the boundary condition. kin variable is temporaraly borrowed.
      kin % n(c) = sqrt(kin % n(c))
    end do

    call GraPhi(kin % n, 1, kin_x, .TRUE.)  ! dK/dx
    call GraPhi(kin % n, 2, kin_y, .TRUE.)  ! dK/dy
    call GraPhi(kin % n, 3, kin_z, .TRUE.)  ! dK/dz

    do c = 1, grid % n_cells

      ! Turning kin back to its real value
      kin % n(c) = kin % n(c) * kin % n(c)
      A % val(A % dia(c)) = A % val(A % dia(c))           &
                          + 2.0 * VISc*(  kin_x(c)**2     &
                                        + kin_y(c)**2     &
                                        + kin_z(c)**2 )   &
                           * grid % vol(c) / (kin % n(c) + TINY)
    end do
  end if

  if(SIMULA == HYB_PITM) then
    do c = 1, grid % n_cells

      ! Production:
      Pk(c)= VISt(c) * Shear(c) * Shear(c)
      b(c) = b(c) + Pk(c) * grid % vol(c)

      ! Dissipation:
      A % val(A % dia(c)) = A % val(A % dia(c))      &
                          + DENc(material(c))        &
                          * eps % n(c) / kin % n(c)  &
                          * grid % vol(c)
    end do
  end if

  call Exchange(grid, kin % n)

  end subroutine

max_stars_repo_name: mhadziabdic/T-FlowS

max_issues_repo_name: mhadziabdic/T-FlowS

max_forks_repo_name: mhadziabdic/T-FlowS


text:
      SUBROUTINE fouri (grpmn, gsource, amatrix, amatsq, bvec,
     &                  bvecNS, ndim)
      USE vacmod
      USE parallel_include_module
      USE timer_sub
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: ndim
      REAL(dp), DIMENSION(mnpd2,nuv3), INTENT(IN) :: grpmn
      REAL(dp), DIMENSION(nuv), INTENT(IN) :: gsource
      REAL(dp), DIMENSION(mnpd,mnpd,ndim**2), INTENT(OUT) :: amatrix
      REAL(dp), DIMENSION(mnpd2,mnpd2), INTENT(OUT) :: amatsq
      REAL(dp), DIMENSION(mnpd,ndim), INTENT(INOUT) :: bvec, bvecNS
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
C     interior  (int_ext=-1), exterior  (int_ext=+1)  neumann problem
      REAL(dp), PARAMETER :: int_ext = 1
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: i, j, k, m, mn, mn0, n
      REAL(dp), ALLOCATABLE, DIMENSION(:) :: source
      REAL(dp) :: ton, toff, tfourion, tfourioff
!-----------------------------------------------
!
!     AMATRIX(,1) = A(Sin)(Sin');  AMATRIX(,2) = A(Sin)(Cos');
!     AMATRIX(,3) = A(Cos)(Sin');  AMATRIX(,4) = A(Cos)(Cos')
!
!     ARG OF TRIG FUNCTIONS: m'u' - n'v' CORRESPONDING TO THE PRIMED MESH IN
!     PKM EQ.(2.14), AND mu  - nv (UNPRIMED MESH) IN EQ. (2.16)
!
!     ON ENTRY, GRPMN(MN,...) HAS BEEN FOURIER-TRANSFORMED OVER THE UNPRIMED
!     COORDINATES (IN FOURP), SO THE FIRST INDEX OF GRPMN CORRESPONDS TO THE FIRST
!     INDEX OF AMATRIX. THUS, THE FOURIER TRANSFORMS OF GRPMN HERE ARE OVER THE PRIMED MESH.
!
!     IN CONTRAST, THE INTEGRAL OF THE SOURCE TERM OVER THE PRIMED MESH WAS ALREADY
!     DONE (IN SCALPOT), SO HERE THE FT ARE OVER THE UNPRIMED MESH FOR THE SOURCE.
!
      CALL second0(tfourion)

      ALLOCATE (source(nuv3), stat=i)
      IF (i .NE. 0) STOP 'Allocation error in fouri'
!
!     STELLARATOR-SYMMETRIZE SOURCE TERMS (with respect to u,v and -u,-v)
!     INDEX (1) IS ANTI-SYMMETRIC, INDEX (2) IS SYMMETRIC
!
!     GSOURCE = - (2pi)**2 B * dS (h - hsing) * WINT
!
!     WINT: needed to normalize integral over angles to unity
!
      IF (lasym) THEN
         source(nuv3min:nuv3max) = onp*gsource(nuv3min:nuv3max)
      ELSE
         k = 0
         DO i = 1, nu2
            DO j = 1, nv
               k = k + 1
               IF (nuv3min.LE.k .AND. k.LE.nuv3max) THEN
                  source(k) = p5*onp*(gsource(k) - gsource(imirr(k)))
               END IF
            END DO
         END DO
      END IF

!
!     INITIALIZE RUNNING-SUM ARRAYS
!
      bvecNS = 0
      amatrix = 0

!
!     PERFORM M,N TRANSFORMS
!
      mn = 0
      NLOOP2: DO n = -nf, nf
         MLOOP: DO m = 0, mf
            mn = mn + 1
            j = 0
            IF (m.EQ.0 .AND. n.LT.0) CYCLE
            DO i = nuv3min, nuv3max
               j = j + 1
               bvecNS(mn,1) = bvecNS(mn,1) + sinmni(j,mn)*source(i)

               amatrix(:,mn,1) = amatrix(:,mn,1)
     &                         + sinmni(j,mn)*grpmn(:mnpd,i)

               IF (.NOT.lasym) CYCLE

               bvecNS(mn,2) = bvecNS(mn,2) + cosmni(j,mn)*source(i)

               amatrix(:,mn,2) = amatrix(:,mn,2)
     &                         + cosmni(j,mn)*grpmn(:mnpd,i)
               amatrix(:,mn,3) = amatrix(:,mn,3)
     &                         + sinmni(j,mn)*grpmn(mnpd+1:,i)
               amatrix(:,mn,4) = amatrix(:,mn,4)
     &                         + cosmni(j,mn)*grpmn(mnpd+1:,i)
            END DO
         END DO MLOOP
      END DO NLOOP2

      IF (vlactive) THEN
         CALL second0(ton)
         CALL MPI_Allreduce(MPI_IN_PLACE, amatrix, SIZE(amatrix),
     &                      MPI_REAL8, MPI_SUM, VAC_COMM, MPI_ERR)
         CALL second0(toff)
         allreduce_time = allreduce_time + (toff - ton)
      END IF
!
!     ADD (still not reduced) ANALYTIC AND NON-SINGULAR PARTS
!
      bvec = bvec + bvecNS

      DEALLOCATE (source, stat=i)

      mn0 = 1 + mf1*nf                                            !Index of m,n=(0,0)

!SANITY CHECKS
!      IF (ANY(bvec(1:mn0-mf1:mf1,:) .NE. 0._dp)) STOP 'BVEC != 0'
!      IF (ANY(amatrix(:,1:mn0-mf1:mf1,:) .ne. 0._dp)) STOP 'AMAT1 != 0'
!      IF (ANY(amatrix(1:mn0-mf1:mf1,:,:) .ne. 0._dp)) STOP 'AMAT2 != 0'
!
!     ADD DIAGONAL TERMS TO AMATRIX [THE FIRST TERM IN EQ(3.2) OF PKM]
!
      DO mn = 1, mnpd
         amatrix(mn,mn,1) = amatrix(mn,mn,1) + pi3*int_ext
      END DO

      IF (lasym) THEN
         DO mn = 1, mnpd
            amatrix(mn,mn,4) = amatrix(mn,mn,4) + pi3*int_ext
         END DO
         amatrix(mn0,mn0,4) = amatrix(mn0,mn0,4) + pi3*int_ext     !!COS(0-0) mode *2
      END IF

!
!     PUT ELEMENTS INTO SQUARE MATRIX
!
      amatsq(:mnpd,:mnpd) = amatrix(:,:,1)                      !Sin-Sin'

      IF (lasym) THEN
         amatsq(:mnpd,1+mnpd:mnpd2) = amatrix(:,:,2)               !Sin-Cos'
         amatsq(1+mnpd:mnpd2,:mnpd) = amatrix(:,:,3)               !Cos-Sin'
         amatsq(1+mnpd:mnpd2,1+mnpd:mnpd2) = amatrix(:,:,4)        !Cos-Cos'
      END IF

      CALL second0(tfourioff)
      timer_vac(tfouri) = timer_vac(tfouri) + (tfourioff-tfourion)
      fouri_time = timer_vac(tfouri)

      END SUBROUTINE fouri

ext: f

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT


text:
! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail,com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE GET_ELGP ( INT_ELEM_ID )

! Gets number of grid points for a given element based on the element's internal ID

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG
      USE IOUNT1, ONLY                :  ERR, F04, F06, WRT_ERR, WRT_LOG
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FATAL_ERR, MEFE, MELGP, METYPE
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  GET_ELGP_BEGEND
      USE MODEL_STUF, ONLY            :  EDAT, EID, ELGP, ELMTYP, etype, EMG_IFE, EPNT, ERR_SUB_NAM, NELGP, NUM_EMG_FATAL_ERRS, TYPE

      USE GET_ELGP_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'GET_ELGP'
      CHARACTER( 1*BYTE)              :: FOUND       = 'N' ! 'Y' if we find the requested element tyoe

      INTEGER(LONG), INTENT(IN)       :: INT_ELEM_ID       ! Internal element ID
      INTEGER(LONG)                   :: EPNTK             ! Value from array EPNT at the row for this internal elem ID. It is the
!                                                            row number in array EDAT where data begins for this element.
      INTEGER(LONG)                   :: I                 ! DO loop index
      INTEGER(LONG)                   :: NG                ! Number of GRID's for USERIN elem
      INTEGER(LONG)                   :: NS                ! Number of SPOINT's for USERIN elem
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = GET_ELGP_BEGEND

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
      EPNTK = EPNT(INT_ELEM_ID)
      FOUND = 'N'
      TYPE  = ETYPE(INT_ELEM_ID)
      EID   = EDAT(EPNTK)

      IF (TYPE == 'USERIN  ') THEN

         NG = EDAT(EPNTK+2)
         NS = EDAT(EPNTK+3)
         ELGP = NG + NS
         FOUND = 'Y'

      ELSE

         DO I=1,METYPE
            IF (TYPE == ELMTYP(I)) THEN
               ELGP  = NELGP(I)
               FOUND = 'Y'
            ENDIF
         ENDDO

      ENDIF

! If we didn't find a valid element type write error and quit

      IF (FOUND == 'N') THEN
         NUM_EMG_FATAL_ERRS = NUM_EMG_FATAL_ERRS + 1
         FATAL_ERR = FATAL_ERR + 1
         IF (WRT_ERR > 0) THEN
            WRITE(ERR,1940) SUBR_NAME, TYPE
            WRITE(F06,1940) SUBR_NAME, TYPE
         ELSE
            IF (NUM_EMG_FATAL_ERRS <= MEFE) THEN
               ERR_SUB_NAM(NUM_EMG_FATAL_ERRS) = SUBR_NAME
               EMG_IFE(NUM_EMG_FATAL_ERRS,1) = 1940
            ENDIF
         ENDIF
         CALL OUTA_HERE ( 'Y' )
      ENDIF

! Check ELGP against max allowable (for coding error)

      IF ((ELGP < 1) .OR. (ELGP > MELGP)) THEN
         NUM_EMG_FATAL_ERRS = NUM_EMG_FATAL_ERRS + 1
         FATAL_ERR = FATAL_ERR + 1
         IF (WRT_ERR > 0) THEN
            WRITE(ERR,1933) SUBR_NAME, TYPE, EID, ELGP, MELGP
            WRITE(F06,1933) SUBR_NAME, TYPE, EID, ELGP, MELGP
         ELSE
            IF (NUM_EMG_FATAL_ERRS <= MEFE) THEN
               ERR_SUB_NAM(NUM_EMG_FATAL_ERRS) = SUBR_NAME
               EMG_IFE(NUM_EMG_FATAL_ERRS,1) = 1933
               EMG_IFE(NUM_EMG_FATAL_ERRS,2) = ELGP
               EMG_IFE(NUM_EMG_FATAL_ERRS,3) = MELGP
            ENDIF
         ENDIF
         CALL OUTA_HERE ( 'Y' )
      ENDIF

! **********************************************************************************************************************************
 1933 FORMAT(' *ERROR  1933: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' THE NUMBER OF GRID POINTS FOR ',A,' ELEMENT ',I8,', IS: ',I8                                          &
                    ,/,14X,' THE MINIMUM NUMBER IS 1 AND THE MAXIMUM NUMBER ALLOWED FOR ANY ELEMENT IS: ',I8)

 1940 FORMAT(' *ERROR  1940: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' ELEMENT TYPE "',A,'" NOT FOUND IN ARRAY ELMTYP')

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

      END SUBROUTINE GET_ELGP

max_stars_repo_name: luzpaz/MYSTRAN

max_issues_repo_name: luzpaz/MYSTRAN

max_forks_repo_name: luzpaz/MYSTRAN


text:
c
c----67--1---------2---------3---------4---------5---------6---------7-2
c Write the date and time to the given file unit number.
c Useful to mark the log files to check that the file is current.
c
      subroutine writeDate(io)
      implicit none
c
      character     day*12,tim*8, a*20,b*20,c*20
      character*4   month(12)
      integer       value(10), io,i,j,k
c
      call date_and_time(a,b,c,value)
c
c  Date
c
      data month /'Jan ','Feb ','Mar ','Apr ',
     &            'May ','Jun ','Jul ','Aug ',
     &            'Sep ','Oct ','Nov ','Dec '/
      i = value(2)
      j = value(3)
      k = value(1)
      if (j.lt.10) then
         write(day,'(a4,a1,i1,a1,i4,a1)') month(i),' ',j,' ',k,' '
      else
         write(day,'(a4,i2,a1,i4,a1)')    month(i),j,' ',k,' '
      endif
      write(io,'(t1,a,a)') 'current date = ',day
c
c  Time
c
      i = value(5)
      j = value(6)
      k = value(7)
      tim = '00:00:00'
      if (i.lt.10) then
         write( tim(2:2), '(i1)' ) i
      else
         write( tim(1:2), '(i2)' ) i
      endif
      if (j.lt.10) then
         write( tim(5:5), '(i1)' ) j
      else
         write( tim(4:5), '(i2)' ) j
      endif
      if (k.lt.10) then
         write( tim(8:8), '(i1)' ) k
      else
         write( tim(7:8), '(i2)' ) k
      endif
      write(io,'(t1,a,a)') 'current time = ',tim
c
      return
      end
c
c----67--1---------2---------3---------4---------5---------6---------7-2
c


text:
      SUBROUTINE zabort6 (IFLTAB, IERR, CSUB, JERR, IADD, CCOMET)

C
C     Records abort messages
C     On Harris records errors at the end of file SYST*DSS6-ERR
C
C     Written by Bill Charley at HEC, 1989
C
C     Error Codes:
C        10 - Pointer or address array incorrect
C        11 - Pointer or address array incorrect (Wrong Pathname)
C        20 - Illegal Unit number
C        30 - Error on Physical Read
C        40 - Error on Physical Write
C        41 - Disk Space Exceeded
C        50 - Corrupt IFLTAB array (Keys don't match)
C        60 - Incomplete Buffered write
C        70 - DSS File Not Opened
C       100 - Illegal KTABLE vaiable (IFLTAB Corrupt)
C       110 - Illegal Number of Characters per Machine word Set
C       120 - Time-series header length too small
C       130 - Excess write on read only file errors
C       200 - Unable to make shared assignment
C       210 - Unable to lock file
C       220 - Unalbe to unlock file
C       300 - Uncompatiable Versions
C       310 - Data Compression Used in old code
C       320 - Insufficient Header Space in zwrbuf6
C       330 - Insufficient Data Space in zwrbuf6
C
      INTEGER IFLTAB(*)
      CHARACTER CSUB*(*), CCOMET*(*)
C
      INCLUDE 'zdsskz.h'
C
      INCLUDE 'zdsscz.h'
C
      INCLUDE 'zdssmz.h'
C
      INCLUDE 'zdsslz.h'
C
C
C
C
      WRITE (MUNIT,20) IERR, CSUB
 20   FORMAT (' Error Code',I5,',  detected in function ',A)
      WRITE (MUNIT,22) IADD
 22   FORMAT (' Error detected at address ',I10)
      WRITE (MUNIT,23) CCOMET
 23   FORMAT (' Error: ',A)
      CALL CHRLNB (CKPATH, N)
      IF (N.EQ.0) N = 1
      WRITE (MUNIT, 30) CKPATH(1:N)
 30   FORMAT (' Last Pathname Accessed: ',A)
      WRITE (*,20) IERR, CSUB
      WRITE (*, 30) CKPATH(1:N)
C
C
C
      IF (LNOABT) THEN
         IF (IERRMS.EQ.0) THEN
            IERRMS = IERR
            WRITE(CERRMS, 140)IERR, CCOMET, CSUB
 140        FORMAT ("Error:",I4,"; ",A,", detected in ",A)
         ENDIF
         IF (IERR.NE.11) THEN
            WRITE (MUNIT, 120)
 120        FORMAT (/,' ----DSS:  File place in read access mode only.')
            IFLTAB(KREADO) = 1
            IFLTAB(KSTAT) = IERR
         ENDIF
         RETURN
      ENDIF
C
      CALL FLUSH(MUNIT)                                                 Mu
      CALL ABORT
C
      END

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


text:
Grace Garden raises vegetables and fruits to help feed the less fortunate and hungry of our community.  Ninety percent of our harvest is donated to Friday’s Harvest, a weekly giveaway at Davis Korean Church and to Davis Community Meals.  In 2013 we harvested over 1655 pounds of vegetables and fruits.  We start all our own seeds in the hoop house.

Each year we have two fundraisers, the Annual Plant Sale and a fall calendar sale.  We have Master Gardener classes in the winter, spring and fall.

Mission statement:  Grace Garden is a faith based ministry of Davis United Methodist Church created to feed the hungry and educate the community about gardening, the environment and sustainable practices through teaching and outreach.  We grow community which nourishes the body and spirit.

Workdays are Monday & Saturday 810 am, Thursday 6:308:30 pm.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine startup(tajj, tasii)
!***********************************************************************
!  Copyright, 1993, 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 Perform miscellaneous startup calculations.
CD1
C***********************************************************************
CD2
CD2 REVISION HISTORY
CD2
CD2 Revision                    ECD
CD2 Date         Programmer     Number  Comments
CD2
CD2 07-DEC-93    Z. Dash        22      Initial implementation.
CD2
CD2 $Log:   /pvcs.config/fehm90/src/startup.f_a  $
CD2
!D2
!D2    Rev 2.5   06 Jan 2004 10:44:00   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:17:22   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:14:44   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:27:00   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 12:10:32   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:45:58 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2    Rev 1.33   Fri May 31 15:28:24 1996   gaz
CD2 decrease allocation for b array in 6dof dpdp
CD2
CD2    Rev 1.32   Fri May 17 13:00:24 1996   hend
CD2 Added optional parameter to time macro for initial time
CD2
CD2    Rev 1.31   Fri Apr 26 16:14:36 1996   gaz
CD2 changes for mdnodes
CD2
CD2    Rev 1.30   Wed Apr 03 15:15:52 1996   hend
CD2 Removed New Allocation Unless Tracer Problem
CD2
CD2    Rev 1.29   Thu Mar 21 13:19:24 1996   hend
CD2 Fixed allocation of sehvariables for trac
CD2
CD2    Rev 1.28   Mon Mar 04 16:15:46 1996   hend
CD2 Removed uneccesary calculations from coneq1 and added trac input
CD2 option
CD2
CD2    Rev 1.27   Fri Feb 16 11:36:12 1996   zvd
CD2 Modified requirements.
CD2
CD2    Rev 1.26   Fri Feb 02 12:01:30 1996   hend
CD2 Updated Requirements Traceability
CD2
CD2    Rev 1.25   Fri Jan 12 17:58:50 1996   llt
CD2 changed mmgetblk arguments
CD2
CD2    Rev 1.24   Thu Jan 11 12:51:58 1996   gaz
CD2 fixed requirements for reduced degree of freedom
CD2
CD2    Rev 1.23   Tue Jan 09 14:14:06 1996   llt
CD2 gaz changes
CD2
CD2    Rev 1.22   12/13/95 10:31:18   robinson
CD2 Changed allocation of a matrix to handle multi-species problems
CD2
CD2    Rev 1.21   12/11/95 15:28:50   gaz
CD2 more rdof changes
CD2
CD2    Rev 1.20   12/11/95 14:17:48   gaz
CD2 changed mdof ,idof definition to accomdate reduced degree of
CD2 freedom dpdp
CD2
CD2    Rev 1.19   11/15/95 15:33:18   gaz
CD2 changes to handle irdof with 4dof
CD2
CD2    Rev 1.18   09/18/95 08:49:22   gaz
CD2 corrected (again) estimate for nnop
CD2
CD2    Rev 1.17   09/15/95 09:23:56   gaz
CD2 corrected error in nnop estimate
CD2
CD2    Rev 1.16   09/15/95 09:08:52   gaz
CD2 new estimates for nnop and nbd
CD2
CD2    Rev 1.15   08/16/95 11:00:08   zvd
CD2 Corrected write to iatty when unassigned.
CD2
CD2    Rev 1.14   08/08/95 08:57:26   awolf
CD2 rarng taken out by SEH and moved to infiles.f
CD2
CD2    Rev 1.13   08/03/95 17:08:10   gaz
CD2 added calls to md_nodes for multiply defined nodes
CD2
CD2    Rev 1.12   08/03/95 08:23:02   robinson
CD2 Fixed problem with the writing of the AVS geo file
CD2
CD2    Rev 1.11   08/02/95 16:47:50   llt
CD2 changed allocation of ib array from nnop to nbd
CD2
CD2    Rev 1.10   08/02/95 13:00:02   awolf
CD2 Allocate a_axy here instead of in allocmem now
CD2
CD2    Rev 1.9   06/23/95 14:44:06   gaz
CD2 made nbd and nemx larger
CD2
CD2    Rev 1.8   06/01/95 16:45:50   gaz
CD2 made change to allow idof=6 for h-m-a
CD2
CD2    Rev 1.7   03/29/95 12:35:44   llt
CD2 changed allocation of bp array
CD2
CD2    Rev 1.6   03/20/95 13:34:26   gaz
CD2 moved call to rarng before call to anonp
CD2
CD2    Rev 1.2   03/28/94 16:41:44   robinson
CD2 Removed unneeded array.
CD2
CD2    Rev 1.1   03/18/94 15:55:50   gaz
CD2 Added solve_new and cleaned up memory management.
CD2
CD2    Rev 1.0   01/20/94 10:28:14   pvcs
CD2 original version in process of being certified
CD2
c 17-mar-93
c because we are only allocating space for b at the
c call to the solvers,we must allocate space for array
c ib here
c also the realloc routines are commented out
c we need to address this
c 12/16/94 gaz alloc space for nop just before slvesu
c 12/22/94 gaz read directly in nar for fill-in level
c 12/23/94 gaz got rid of ldn1,ldn2
c          cal resize nop after slvesu
c 12/23/94 gaz printout resized bmatrix size
c 1/3/95   gaz called thikness(1) for 2-d problems l 410
c 1/3/95   gaz commented out x3>x2 l 310
C***********************************************************************
CD3
CD3 INTERFACES
CD3
CD3 Formal Calling Parameters
CD3
CD3   Identifier      Type     Use  Description
CD3
CD3   tajj            REAL*8   I    Elapsed cpu time (for reading input and
CD3                                   coefficient generation)
CD3   tasii           REAL*8   O    Cpu time at start of solution computations
CD3
CD3 Interface Tables
CD3
CD3   None
CD3
CD3 Files
CD3
CD3   Name                     Use  Description
CD3
CD3   iout                     O    File used for general code output
CD3
C***********************************************************************
CD4
CD4 GLOBAL OBJECTS
CD4
CD4 Global Constants
CD4
CD4   None
CD4
CD4 Global Types
CD4
CD4   None
CD4
CD4 Global Variables
CD4
CD4                            COMMON
CD4   Identifier      Type     Block  Description
CD4
CD4   cord            REAL*8   fbs    Contains the coordinates of each node
CD4   cpr             REAL*8   fdd    Rock specific heat at each node
CD4   denei           REAL*8   fcc    Energy accumulation term
CD4   deneh           REAL*8   fdd    Last time step energy accumulation term
CD4                                     at each node
CD4   denej           REAL*8   fdd    Last time step energy accumulation time
CD4                                     derivative at each node
CD4   denh            REAL*8   fdd    Last time step mass accumulation term at
CD4                                     each node
CD4   deni            REAL*8   fcc    Mass accumulation term
CD4   denj            REAL*8   fdd    Last time step mass accumulation time
CD4                                     derivative at each node
CD4   denr            REAL*8   fdd    Rock density at each node
CD4   dit             REAL*8   fdd1   Array containing time step changes
CD4   dstm            REAL*8   fcc    Steam mass
CD4   icnl            INT      faai   Problem dimension
CD4   ig              INT      fhh    Variable order lu decomposition
CD4                                     information
CD4   iieos           INT      fddi1  Thermodynamics set at each node
CD4   iirb            INT      fhh    Inverse of irb
CD4   iout            INT      faai   Unit number for output file
CD4   ipa             POINTER  fee    Pointer to variable array a
CD4   ipistrw         POINTER  fbb    Pointer to variable array istrw
CD4   ipnelm          POINTER  fbb    Pointer to variable array nelm
CD4   ipnop           POINTER  fbb    Pointer to variable array nop
CD4   iptty           INT      faai   Unit number for selected tty output
CD4   irb             INT      fhh    Array containing the reordered node
CD4                                     numbers
CD4   irdof           INT      david1 Reduced degree of freedom model used
CD4   istrw           INT      fbb    Starting positions in sx(nr,9) array of
CD4                                     finite element coefficients for each
CD4   mdof            INT      comrxni The maximum number of degrees
CD4                                    of freedom necessary to solve
CD4                                    the heat & mass and the
CD4                                    tracer solutions
CD4
CD4   nar             INT      fbb    Array containing gauss elimination order
CD4                                     for each node
CD4   neigh           INT      faai   Maximum number of neighbors occur in
CD4                                     tracer solution
CD4   nelm            INT      fbb    Initially information about nodes in each
CD4                                     element, later nodal connectivity
CD4   nemx            INT      faai   Number of unique (geometrically) elements
CD4   neq             INT      faai   Number of nodes, not including dual
CD4                                     porosity nodes
CD4   nop             INT      fbb    Matrix sparsity structure for lu
CD4                                     decomposition
CD4   nopt            INT      fhh    Array indicating active variables
CD4   npvt            INT      fhh    Pivot information for the lu
CD4                                     decomposition matrix
CD4   pflow           REAL*8   fdd    Flowing pressure at each source node
CD4   phi             REAL*8   fdd    Pressure at each node
CD4   phini           REAL*8   fdd2   Initial pressure at each node
CD4   pho             REAL*8   fdd    Last time step pressure at each node
CD4   pnx             REAL*8   fdd    Permeability in the x-direction, liquid
CD4                                     velocity in the x-direction, vapor
CD4                                     velocity in the x-direction
CD4   pny             REAL*8   fdd    Permeability in the y-direction liquid
CD4                                     velocity in the y direction, vapor
CD4                                     velocity in the y-direction
CD4   pnz             REAL*8   fdd    Permeability in the z-direction, liquid
CD4                                     velocity in the z-direction, vapor
CD4                                     velocity in the z-direction
CD4   ps              REAL*8   fdd    Porosity at each node
CD4   psini           REAL*8   fdd2   Initial porosity at each node
CD4   s               REAL*8   fdd    Liquid saturation at each node
CD4   so              REAL*8   fdd    Last time step saturation at each node
CD4   sx1             REAL*8   fbc    Contains volume associated with each node
CD4   t               REAL*8   fdd    Temperature at each node
CD4   thx             REAL*8   fdd    Thermal conductivity x-direction
CD4   thy             REAL*8   fdd    Thermal conductivity y-direction
CD4   thz             REAL*8   fdd    Thermal conductivity z-direction
CD4   tini            REAL*8   fdd2   Initial temperature at each node
CD4   to              REAL*8   fdd    Last time step temperature at each node
CD4   volume          REAL*8   fdd    Volume associated at each node
CD4
CD4 Global Subprograms
CD4
CD4   Identifier      Type     Description
CD4
CD4   airctr                   Manage the isothermal air-water calculations
CD4   anonp                    ?
CD4   bcon                     Adjust/manage boundary conditions

CD4   co2ctr                   Control isothermal air-water simulation
CD4   coeffc                   Change the coefficients of the equation of
CD4                              state fits
CD4   contr                    Write contour plot data
CD4   disk                     Read/write initial/final state data
CD4   dpdp                     Control dual porosity/dual permeability solution
CD4   dual                     Find dual porosity contributions to nodes
CD4   peint                    Initializes pressures and temperatures
CD4   plot                     Write history plot data
CD4   porosi                   Calculate pressure dependant porosity
CD4   radius                   Figure out radius in radial problem
CD4   rarng                    Rearrange coefficients in 2-d problem
CD4   setord                   Set up the order of solution for the equations
CD4                              at each node
CD4   sice                     Control ice simulation
CD4   slvesu                   Perform symbolic factorization and calculate
CD4                              computer storage necessary for the solvers
CD4   split                    Split rectangles(bricks) into 4 triangles(12
CD4                              tetrahedrals)
CD4   steady                   Calculate steady state solution
CD4   storsx                   Store and retrieve element coefficients
CD4   thickness                Adjust thicness in 2-d problems
CD4   tyming          REAL*8   Calculate elapsed cpu time
CD4   varchk                   Determine variable set and make n-r corrections
CD4
C***********************************************************************
CD5
CD5 LOCAL IDENTIFIERS
CD5
CD5 Local Constants
CD5
CD5   None
CD5
CD5 Local Types
CD5
CD5   None
CD5
CD5 Local variables
CD5
CD5   Identifier      Type     Description
CD5
CD5   caz             REAL*8   Dummy argument to function tyming
CD5   i               INT      Loop index
CD5   ib              INT      Temporary array for slvesu
CD5   idum1           INT      Temporary array for slvesu
CD5   idum2           INT      Temporary array for slvesu
CD5   icode            INT      variable used in memory management
CD5   imm             INT      Counter
CD5   nbytes          INT      Space to be allocated for variable arrays
CD5
CD5 Local Subprograms
CD5
CD5   None
CD5
C***********************************************************************
CD6
CD6 FUNCTIONAL DESCRIPTION
CD6
C***********************************************************************
CD7
CD7 ASSUMPTIONS AND LIMITATIONS
CD7
CD7 None
CD7
C***********************************************************************
CD8
CD8 SPECIAL COMMENTS
CD8
CD8 This routine initializes various startup conditions which are
CD8 necessary for the code.
CD8
CD8  Requirements from SDN: 10086-RD-2.20-00
CD8    SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
CD8    FEHM Application Version 2.20
CD8
C***********************************************************************
CD9
CD9 REQUIREMENTS TRACEABILITY
CD9
CD9 2.7   Provide Restart Capability
CD9 2.7.3 Resume the calculation
CD9 2.8   Provide Multiple Realization Option
CD9
C***********************************************************************
CDA
CDA REFERENCES
CDA
CDA None
CDA
C***********************************************************************
CPS
CPS PSEUDOCODE
CPS
CPS BEGIN startup
CPS   ???
CPS END startup
CPS
C***********************************************************************

      use comai
      use combi
      use comci
      use comco2
      use comcouple
      use comdi
      use comdti
      use comei
      use comfi
      use comflow
      use comgi
      use comii
      use commeth
      use compart, only : ptrak
      use comrlp, only : ishisrlp
      use comsi, only : idof_stress
      use comsptr, only : sptrak
      use avsio, only : iovapor
      use comwt
      use comzone
      use davidi
      implicit none

      integer icall, dummyint, j, iconv_tmp
      integer mdof,ldnmax,ldnmaxtmp,ncont,neqp1,igmax
      integer nstorederiv,nstoresolve,i,icode,irdum
      integer ncon_size,ncon_size_orig
      integer open_file
      integer idb_file, if_debug
      real*8 tajj, tasii, tyming
      real*4 caz(2)
      real*8 rnwn1, rnwn2, rnwn3, rnwn, rnwd1, rnwd2, rnwd3, rnwd
      real*8 xmin_bbox, ymin_bbox, zmin_bbox
      real*8 xmax_bbox, ymax_bbox, zmax_bbox, scale_factor
      integer, allocatable :: ib(:)
      integer, allocatable :: idum1(:)
      integer, allocatable :: idum2(:)
      integer, allocatable :: ncon_temp(:)
      integer, allocatable :: nop_temp(:)
      integer i1,i2,ipiv,icnt
      integer count, num_zone, ic_sym
      integer nsbb,idofdum,neidum,idof_save
      real*8 very_small, fac_nop, fac_mult
c symmetry check on = 1, off = 0
      parameter (ic_sym = 0)
      parameter (very_small= 1.d-30)
c
c
c     Calculate rho1grav
c      rho1grav = crl(1,1)*(9.81d-6)
c      rho1grav = 997.*9.81d-6
c**** Deallocate printout information array for variable porosity model
c**** if allocated and not a variable porosity problem
      if (iporos .ne. -4) then
         if (allocated(nskw3)) deallocate(nskw3)
      endif
c
c gaz 032916 turn off vapor flux output for richards equation
c same for avs
c restart uses 0.0 for vapor (in case restart is 2-phase)
c
       if(jswitch.ne.0) then
c gaz debug 041416 (affects rich flxz and his output)
        vflux_flag = .false.
c gaz debug 041416 (affects rich avs and other contour output!)
c gaz debug 041416 now left as is for correct verification output
c        iovapor = 0
       endif

c**** calculate number of neighbors and connected elements ****
      allocate(idum1(n0),idum2(n0))
      if (icnl .eq. 0)  then
c**** three dimensional problem ****
         nemx = 200
         neigh = 200
      else
c**** two   dimensional problem ****
         nemx = 200
         neigh = 200
      end if

      neigh = max0(neigh, nemx)

c**** set z coordinate = y coordinate for 2-d problems ****
c**** also insure gravity direction is correct ****
      if (icnl .ne. 0)  then
         if (igrav .eq. 3)  igrav = 2
      end if

c**** set idof = 3 if stress solution enabled and icnl = 0 (3-d) ****
      if (icnl .eq. 0 .and. idof_stress .lt. 0)  idof = 3

      if (n0 .lt. neq)  then
         write(ierr, 6000)  n0, neq
         if (iout .ne. 0) write(iout, 6000)  n0, neq
         if (iptty .gt. 0)  write(iptty, 6000)  n0, neq
 6000    format(/, 1x, '**** n0(', i6, ') .lt. neq(', i7, ') **** ',
     *        'check parameter statements ****')
         stop
      end if
c**** set idof = 3 if air-water-heat
      if (ico2 .gt. 0) idof = max(idof,3)
      if (idof_co2 .gt. 0) idof = max(idof,idof_co2)
      if (icnl.eq.0) i1 = 3
      if (icnl.ne.0) i1 = 2
      if (idof_stress.gt.3) idof = max(idof + i1,5)
      if (istrs.ne.0) idof = max(idof,i1)
c****

      if(idpdp.ne.0.and.idof_stress.ne.0) then
        if(iptty.ne.0) then
         write(iptty,*) '>>>>> dpdp not implemented with stress <<<<'
         write(iptty,*) 'stopping in startup'
        endif
        if(iout.ne.0) then
         write(iout,*) '>>>>> dpdp not implemented with stress <<<<'
         write(iout,*) 'stopping in startup'
        endif
        stop
      endif
c********* Make dpdp, two-phase work properly
c check for dpdp(ngas or 2dof)
      if(idpdp .ne. 0) then
c reorder id irdof gt 0
         if(irdof.gt.0.and.islord.eq.0) islord=1
         if(ico2 .gt. 0) then
            idof=6
         else
            idof=4
         endif
      endif

c calculate new idof based on irdof
      mdof=idof
      if(idof.eq.2) then
         if(irdof.gt.0) mdof=1
      else if(idof.eq.3) then
         if(irdof.eq.1) mdof=1
         if(irdof.eq.2) mdof=2
      else if(idof.eq.4) then
         if(irdof.gt.0) mdof=2
      else if(idof.eq.6) then
         if(irdof.eq.1) mdof=2
         if(irdof.eq.2) mdof=4
      endif
c
c Find the maximum degrees of freedom necessary by taking the
c max of the degrees of freedom needed for the heat and mass solution
c and the max degrees of freedom necessary for the tracer solution
      if (compute_flow) then
         mdof = max(mdof,mdof_sol)
      else
         if (iccen .eq. 1)  then
            mdof = mdof_sol
         else
            mdof = 1
         end if
      end if
c
c for stress simulation isotropic coefficients not allowed
c
      if (istrs.ne.0) then
        isoy = 2
        isoz = 3
      endif
c**** set eos cofficient set eq to 1 since fit was good ****
c****          the entire range only done if iieos(i) = 0 ****
      if (iwelb .ne. 1)  then
         do i = 1, n
            if (iieos(i) .eq. 0)  iieos(i) = 1
         end do
      end if

c**** compatability scaling (viscosity is in pa-sec, not in mpa-sec) ***
      do i = 1, n
         if (idoff .gt. 0) then
            pnx(i) = pnx(i) * 1.0d+06
            pny(i) = pny(i) * 1.0d+06
            pnz(i) = pnz(i) * 1.0d+06
         end if
         if(ico2.ge.0.or.ice.ne.0) then
           thx(i) = thx(i) * 1.0d-06
           thy(i) = thy(i) * 1.0d-06
           thz(i) = thz(i) * 1.0d-06
         endif
         if (cpr(i) .gt. 1.0)  cpr(i) = cpr(i) * 1.0d-06
      end do

c**** set permeabilities = 0 for porosity = 0.0 ****
      neq_active = 0
      do i = 1, n
         if (abs(ps(i)) .le. zero_t .and. idoff .gt. 0)  then
            pnx(i) = zero_t
            pny(i) = zero_t
            pnz(i) = zero_t
            neq_active = neq_active +1
         end if
      end do

      days = 0.0
      mink = neq - neq_active

c**** complete element information ****
!      if(contim.ge.0) then
         call contr (0)
!      else
!         call contr_days (0)
!      endif

c**** complete pest information(if interpolation is used) ****
c
      call pest(-1)
c
c gaz temp call to split 072613
c      ivf = -1
c      call split(0)
c      write(ierr,*) 'stopping startup line 579'
c      if(iout.ne.0) write(iout,*) 'stopping startup line 579'
c      stop
c
      if (mlz .ne. 0. and. irun .eq. 1 .and. lda .le. 0.and.ianpe.eq.0)
     &    call split(0)
      if (lda .le. 0.and.irun.eq.1) then
         if(ivf.ne.-1) then
            if(ianpe.eq.0) then
             call anonp
            else if(ianpe.ne.0) then
             call anonp_ani
             if (ico2.ge.0) call fd_calc_heat(1)
            endif
         else
            call structured(1)
            call structured(3)
         endif
      endif
c mlz is now the size of sx(mlz, j)
c      if (mlz .ne. 0 .and. irun .eq. 1 .and. lda .le. 0) call split(0)
      if (mlz .ne. 0 .and. irun .eq. 1 .and. lda .le. 0) call split(1)
      if (lda .ne. 0.and.irun.eq.1) call storsx
c gaz 11-09-2001
      if(lda.le.0.and.imdnode.ne.0.and.irun.eq.1) then
         ncon_size=nelm(neq_primary+1)
         ncon_size_orig = ncon_size
         nr=nr+1
         allocate(ncon_temp(ncon_size))
         do i=1,ncon_size
            ncon_temp(i)=nelm(i)
         enddo
         deallocate(nelm)
         call md_nodes(4,0,ncon_size)
c gaz 1-24-2003
         nelmd = ncon_size
         allocate(nelm(ncon_size))
         do i=1,ncon_size_orig
            nelm(i)=ncon_temp(i)
         enddo
         deallocate(ncon_temp)
         call md_nodes(5,0,ncon_size)
      endif

c gaz 11-11-2001 moved this before call to add_gdpm
      if (icnl .ne. 0)  then
         call radius
         call thickness(1)
      end if
c
c subroutine called to eliminate negative areas. This is
c only done after the coefficients are generated or read.
c the modified coeefficients are never saved to a "stor" file
c
c gaz 11-25-11
      if(icoef_neg.ne.0) then
       call coeff_management(1)
      endif
c
c
c     Call routine to adjust the connectivity array if needed
c     to add gdpm nodes
c     RJP 12/13/06 modified to add river_ctr for wellbore
c trying to maintain compatibility of add_gdpm and river_ctr
      if (irun.eq.1) then
c add gdpm connections
c calculate the dimensions on a control volume
c neq_primary still defined as gdpm primary grid
c gaz 122311 moved lower
c         call area_length_calc(3)
c gaz 051616 change calls for gdkm
c need drxg,dryg,drzg here
c gaz 091118 removed this line
c          i = sx(1,1)

          if (gdkm_flag .ne. 0.and.gdkm_flag .le. 3)
c calculate gridblock dimensions
     &       call gdkm_volume_fraction_interface(1)
         if (gdpm_flag .ne. 0) call add_gdpm
c allocate memory for gdkm volume factors
         if (gdkm_flag .ne. 0.and.gdkm_flag .le. 3)
     &       call gdkm_volume_fraction_interface(-1)
c         if (gdkm_flag .ne. 0.and.gdkm_flag .le. 3) call gdkm_connect(3)
c        gdkm_connect now called from add_gdpm1
c gaz 081617
         if (gdkm_flag .ne. 0.and.gdkm_flag .le. 3)
     &       call gdkm_volume_fraction_interface(3)
c
c calculate connectivity for generalized permeability mode
c
         if (gdkm_flag .ne. 0.and.gdkm_flag .le. 3) call gdkm_calc(1)
c
c neq_primary now defined as full grid
          if (gdkm_flag .ne. 0.and.gdkm_flag .le. 3)
     &       call gdkm_volume_fraction_interface(0)
         if (gdkm_flag .ne. 0.and.gdkm_flag .le. 3)
     &       call gdkm_volume_fraction_interface(4)
c
c gaz 081116 apply volume fraction to some flow and transport parameters (thx,thy,thz,diff)
c         call gdkm_volume_fraction_apply(1)
c gaz 11-09-2001 allocation of istrw_itfc and istrw_cold
c done here after call to anonp,storsx, or structured
c still could ne changed in add_gdpm
       ncon_size=nelm(neq+1)
       nelmd = ncon_size
      if(idpdp.eq.0) then
         if (.not. allocated (istrw_itfc)) then
          allocate(istrw_itfc(ncon_size))
          istrw_itfc = 0
         endif
         if (.not. allocated (istrw_cold)) then
            allocate(istrw_cold(ncon_size))
            istrw_cold = 0
         endif
      else
         if (.not. allocated (istrw_itfc)) then
             allocate(istrw_itfc(2*ncon_size))
             istrw_itfc = 0
         endif
         if (.not. allocated (istrw_cold)) then
             allocate(istrw_cold(2*ncon_size))
             istrw_cold = 0
         endif
      end if
c add river or well connections
         if (nriver .ne. 0) then
c river_ctr(1) call is made in incoord
c        call river_ctr(1)
            call river_ctr(2)
            call river_ctr(3)
c printout to file volumes,areas etc.
            call river_ctr(-3)
             call river_ctr(4)
c printout to file connectivities
            call river_ctr(-4)
         endif
      end if
c gaz moved 122311 from above
       call area_length_calc(3)
c
c  calculate gridblock height (if necessary) for subsidence calculations
c  (ico2 chance so embedded call to airctr in subsidence) goes through
      i = ico2
      ico2 = -1
      call subsidence(-1)
      ico2 = i
c
      if(interface_flag.ne.0) call setconnarray
c gaz 090618
c call sx_combine to break connections to fixed type BCs
c      if(ianpe.eq.0) then
c         if (irun.eq.1.and.inobr.eq.0) call sx_combine(1)
c      else
c         if (irun.eq.1.and.inobr.eq.0) call sx_combine_ani(1)
c      endif
c
c call fluxo now to calculate neighbors if necessary
c
      call flxo(1)

c allocate space for velocity and alpha connection terms (seh)
      if (iccen .ne. 0) then
         sehsize=0
         do i=1,neq
            i1=nelm(i)+1
            i2=nelm(i+1)
            ipiv=nelmdg(i)
            do icnt=ipiv+1,i2
               sehsize=sehsize+1
            enddo
         enddo
         if (idpdp.ne.0) sehsize=2*sehsize
         if(irun.eq.1) then
            allocate(alphaconl(sehsize),alphaconv(sehsize))
            allocate(sehvell(sehsize),sehvelv(sehsize))
         end if
         alphaconl=0
         alphaconv=0
         sehvell=0
         sehvelv=0
         sehdonevel=0
      endif

      iad = 0
c     Moved the head calls to after the initialization flag w/ iflg=-1
c**** print coordinates and elements to graphics files ****
c     Moved this call until after rarng is called - BAR 8-3-95
c
c     Moved earlier so that a_axy can be read in from disk
c
c       PHS  3/24/2004  adding a_wvxy to store water vapor mass flux
c
      neqp1=neq+1
      ncont=nelm(neqp1)
      ldna=ncont-neqp1
      if(irun.eq.1) then

         if(idof_co2.gt.0) then
            allocate(c_axy(ldna))
            allocate(c_vxy(ldna))
            c_axy = 0.d0
            c_vxy = 0.d0
         endif

         if (idoff .gt. 0) then
            if (idpdp .ne. 0) then
               allocate(a_axy(2*ldna+neq))
            else
               allocate(a_axy(ldna))
            end if

c s kelkar 3 July 2014, for calculating heat flow vectors
            if(flag_heat_out) then
               if (idpdp .ne. 0) then
                  flag_heat_out = .false.
       write(iptty,*),'heat out option currently only single porosity'
       write(ierr,*),'heat out option currently only single porosity'
               else
                  allocate(e_axy_adv(ldna))
                  allocate(e_axy_cond(ldna))
                  allocate(e_adv_nodal(neq,3))
                  allocate(e_cond_nodal(neq,3))
                  e_axy_adv = 0
                  e_axy_cond = 0
                  e_adv_nodal = 0
                  e_cond_nodal = 0
               end if
            endif

            if (irdof .ne. 13) then
               if (idpdp .ne. 0) then
                  allocate(a_vxy(2*ldna+neq))
                  if (iadif .eq. 1) then
                     allocate(a_wvxy(2*ldna+neq))
                  else
                     allocate(a_wvxy(1))
                  end if
               else
                  allocate(a_vxy(ldna))
                  if (iadif .eq. 1) then
                     allocate(a_wvxy(ldna))
                  else
                     allocate(a_wvxy(1))
                  end if
               end if
            else
               allocate(a_vxy(1), a_wvxy(1))
            end if
         else
            allocate(a_axy(1),a_vxy(1),a_wvxy(1))
            if (idoff .eq. -1 .and. flag_heat_out) then
	        allocate(e_axy_adv(l))
	        allocate(e_axy_cond(ldna))
	        allocate(e_adv_nodal(1,3))
	        allocate(e_cond_nodal(neq,3))
	        e_axy_adv = 0
	        e_axy_cond = 0
	        e_adv_nodal = 0
	        e_cond_nodal = 0
            end if
         end if
      end if
      a_axy=0
      a_vxy=0
      a_wvxy=0

c new placement of peint
      if (i_init.ne.0) call peint
      if (iread .gt. 0) call diskread
      if(isalt.ne.0) call saltctr(20,0,0.0d00,0.0d00)
      close (1010)

c added new feature -- set initial time if requested in time macro
       if (irsttime.ne.0) days=rsttime
c gaz 092111 have moved  above
c      if (iread .le. 0 .and. (tin0 .gt. 0.0 .or. tin1.gt.0.0))
c    &     call peint
c      if (iread.le.0 .and. igrad. ne.0) call gradctr(1)
c     gaz 9-26-04
      if (igrad. ne.0) call gradctr(1)
      if (iconv. ne.0) call convctr(1)
      if (isubwd. ne.0) call wellimped_ctr(1)
c
c check that source terms are ok on restart for wellbore model
c      call welbor(-1)

c      set clathrate properties equal to water properties for a test
c
      call icectr(-2,0)

      call icectrco2(-2,0)

c gaz 11-11-2001 moved before call to add_gdpm
c     if (icnl .ne. 0)  then
c        call radius
c        call thickness(1)
c     end if

      dtot = day * 86400.0

      if (sx1(1) .gt. 0.0)  then
         upwgt = 1.0 - upwgt
         upwgta = 1.0 - upwgta
      end if
      dnwgt = 1.0 - upwgt
      dnwgta = 1.0 - upwgta

c Compute space for gmres

c check if bcgs solution requested
      if(accm.eq.'bcgs') then
         kgmres = (3 + 2 * ( abs(north)+1 ) )*mdof*neq
         if (iout .ne. 0) write(iout,*) 'bcgstab method chosen '
         if (iatty.ne.0) write(iatty,*) 'bcgstab method chosen '
      elseif(accm.eq.'gmre') then
         kgmres = (north + 1) * mdof * neq
         if (iout .ne. 0) write(iout,*) 'gmres method chosen '
         if (iatty.ne.0) write(iatty,*) 'gmres method chosen '
      endif
c adjust space for a  matrix when it is called
      ldn=ldna*idof**2

c**** set up equation solvers ****
c  estimate parameters for LU arrays
      igmax=0
      do i=1,neq
         igmax=max(igmax,nar(i))
      enddo
c      nnop = neq*(ldna/neq+1)*igmax**2.1 +neqp1
c new estimate GAZ 080707
       if(igmax.le.1) then
	  nnop = ncont
	 else
	  fac_nop = 1.5
	  fac_mult = 2.2
        nnop = fac_nop*ldna*fac_mult**(igmax-1) + neqp1
	 endif
      if(ncont.gt.nnop) nnop=ncont
      nbd  = nnop*mdof**2
      if(irdof.gt.0) then
         irdum=mdof
      else
         irdum=0
      endif
      ldnmax = max(ldn,ldna*mdof_sol**2)
      ldnmaxtmp = ldnmax

c     ib needs only nnop storage (dof independent)
      if(irun.eq.1) then
         allocate(ib(nnop),nop(nnop))
c
c check for input inconsistencies
c
         call diagnostics(0)
c
c Compute approximate bounding box of each node
c before modification/simplification of connectivity.
c
         allocate(dx_bbox(neq))
         allocate(dy_bbox(neq))
         allocate(dz_bbox(neq))
C         if((icnl .eq. 4).or.(icnl .eq. 5).or.(icnl .eq. 6))then
         allocate(da_bbox(neq))
C         endif
         do i = 1, neq
            call node_midedge_box_size
     &           (dx_bbox(i),dy_bbox(i),dz_bbox(i),
     &           xmin_bbox, ymin_bbox, zmin_bbox,
     &           xmax_bbox, ymax_bbox, zmax_bbox, i)
c
c        Check if the prblem is radial r,z instead of xyz
c
            if((icnl .eq. 4) .or. (icnl .eq. 5))then
               da_bbox(i) = (xmax_bbox**2 - xmin_bbox**2)*acos(-1.0)
            elseif(icnl .eq. 6)then
               da_bbox(i) = (ymax_bbox**2 - ymin_bbox**2)*acos(-1.0)
            else
               da_bbox(i) = dx_bbox(i)*dy_bbox(i)
            endif
         enddo
C
         if_debug = 0
         if(if_debug .ne. 0)then
            idb_file = open_file('tmp_debug_dx_dy_dz.tbl','unknown')
            write(idb_file,*)'dx_bbox, real'
            write(idb_file,*)'dy_bbox, real'
            write(idb_file,*)'dz_bbox, real'
            write(idb_file,*)'da_bbox, real'
            write(idb_file,*)'dv_bbox, real'
            write(idb_file,*)'sx1    , real'
            write(idb_file,*)'vrat   , real'
            do i = 1, neq
               write(idb_file,'(i8,9e15.6)')
     1              i,dx_bbox(i),dy_bbox(i),dz_bbox(i),
     1              da_bbox(i), dx_bbox(i)*dy_bbox(i)*dz_bbox(i),
     2              sx1(i),(dx_bbox(i)*dy_bbox(i)*dz_bbox(i))/sx1(i)
            enddo
            close(idb_file)
         endif
         if_debug = 0
c
c s kelkar may 20 09 moved call to ptrac1 here to do the geometry
c before simplify_ncon for -ve porosities. But moved the call to
c load_omr_flux_array from ptrac1 to fehmn.f where ptrac1 used to be called
         if(sptrak) then
            call ptrac1
         endif

c simplify connectivity based on porosity=0.0
c only for isotropic problems (ianpe=0)
c
        i1 = 0
        do i = 1,n
         if(nelm(i)+1.eq.nelm(i+1).and.ps(i).gt.0.0d0) then
          i1 = i1 + 1
          write(ierr,*) 'WARNING isolated node (w por> 0) ',i
	   endif
       enddo
       if(i1.gt.0) then
        if(iptty.ne.0) write(iptty,*)
     &   'WARNING isolated nodes with por> 0 exist: see fehmn.err'
        if(iout.ne.0) write(iout,*)
     &   'WARNING isolated nodes with por> 0 exist: see fehmn.err'
       endif
         if(ianpe.eq.0) then
            call   simplify_ncon
     &           (0,nelm,nelmdg,nop,istrw,ib,neq,idof,ka,ps,i)
            if (connect_out) call connections_list
            if (igauss .le. 1) then
               deallocate (nop)
            else
               nop = 0
            end if
            neq_active = neq
            if(i.gt.0) then
               neq_active = neq -i
               if (iout .ne. 0) write(iout, 100) i
               if (ischk .ne. 0) write(ischk, 100) i
               if (iptty .ne. 0) write(iptty, 100) i
 100           format ('>>>> ', i8, ' nodes eliminated (porosity <0)')
               do i=1,neq
                  if(nelm(i)+1.eq.nelm(i+1)) then
                     ka(i) = 0.0
                     sk(i) = 0.0
                     qh(i) = 0.0
                     qc(i) = 0.0
                     if (idoff .gt. 0) then
                        pnx(i)=very_small
                        pny(i)=very_small
                        pnz(i)=very_small
                     end if
                     pho(i)=0.1
c gaz 042119 set to to 0.0 (will be changed later to tref
                      to(i)=0.0
                     if (irdof .ne. 13 .or. ifree .ne. 0) s(i)=1.0
                  endif
               enddo
            endif
         end if
      end if
c
c check connectivity
c
       if(ic_sym.eq.1)call connectivity_symmetry_test(1)
c
c call reordering algorithm here
c
      call renum(1)
c
      if(irun.eq.1) then
 11      if(nbd.lt.0) then
            nbd = 1.5*abs(nbd)
            if(iout.ne.0) write(iout,*) '>> nbd resized to ',nbd,' <<'
            if(iptty.ne.0) write(iptty,*) '>> nbd resized to ',nbd,' <<'
         endif
         if (igauss .gt. 1) then
            if (gdpm_flag.ne.0) then
               call slvesu (neq_primary, ib, nelm_primary, nop,
     *              nelmdg_primary, npvt, nar, irb, iirb,
     *              idum1, idum2, ldnmaxtmp, nbd, kgmres, nelmd, nnop,
     *              mdof, abs(north), ireord, iout, iatty,
     *              irdum, nbnd,accm)
               if(nbd.lt.0) go to 11
            else
               call slvesu (neq, ib, nelm, nop, nelmdg, npvt, nar, irb,
     *              iirb, idum1, idum2, ldnmaxtmp, nbd, kgmres,
     *              nelmd, nnop, mdof, abs(north), ireord,
     *              iout, iatty, irdum, nbnd,accm)
               if(nbd.lt.0) go to 11
            endif
            allocate (nop_temp(nnop))
            nop_temp(1:nnop) = nop(1:nnop)
            deallocate (nop)
            allocate (nop(nnop))
            nop = nop_temp
            deallocate (nop_temp)
         else
! Skip call to slvesu
            do j=1,neq
               iirb(j) =  j
               irb (j) =  j
            end do

            nbnd   =  0
            do j=2,neq
               nbnd =  max0( nelm(j+1)-nelm(j-1),nbnd )
            end do
            if(gdpm_flag.ne.0) then
               do j=1,neq_primary
                  npvt(j) = nelmdg_primary(j)
               enddo
            else
               do j=1,neq
	          npvt(j) = nelmdg(j)
               enddo
            endif
         endif
         deallocate(ib)
      end if
c     reset some array sizes
c GAZ 10/9/98 adjust space for 1dof problem
c GAZ 02/6/01 adjust space for gdpm problem
c
      if(irdof.eq.13) then
         if (iout .ne. 0) write(iout, 6009)
         if (iptty .ne. 0) write(iptty, 6009)
 6009    format('SZ only simulation, true size of a matrix :')
      endif
      if(gdpm_flag.ne.0) then
         if (iout .ne. 0) write(iout, 6012)
         if (iptty .ne. 0) write(iptty, 6012)
 6012    format('>>>> gdpm invoked, true size of matrices :')
      endif
      if(gdpm_flag.ne.0) then
         ldnmax = max(ldn,ldna*mdof_sol**2)
         if (iout .ne. 0) write(iout, 6011) ldnmax
         if (iptty .ne. 0) write(iptty, 6011) ldnmax
      endif
      if(irdof.eq.13) then
         ldn=ldna
         ldnmax = max(ldn,ldna*mdof_sol**2)
         if (iout .ne. 0) write(iout, 6011) ldnmax
         if (iptty .ne. 0) write(iptty, 6011) ldnmax
      endif
 6011 format(1x,'storage available for a matrix resized to ',
     &     i10, '<<<<<<')
c  gaz 032208 change to richards eq
      if(irun.eq.1) then
         if(jswitch.ne.0) then
            call airctr(-2,0)
         else if (compute_flow .or. iccen .eq. 1) then
            allocate(a(ldnmax))
         endif
      end if
c modify storage of LU factorization matrix for 6-2
c 3-1 full GMRES schemes
      if(compute_flow .or. iccen .eq. 1) then
         if (igauss .gt. 1) then
          if(gdpm_flag.ne.0) then
            nsbb = nop(neq_primary+1)-(neq_primary+1)
         else
            nsbb = nop(neq+1)-(neq+1)
          endif
         else
           if(gdpm_flag.ne.0) then
            nsbb = nelm(neq_primary+1)-(neq_primary+1)
           else
            nsbb = nelm(neq+1)-(neq+1)
           endif
         end if
         if (jswitch.eq.1 .and. idpdp .eq. 0) then
            nbd = nsbb
         else if (jswitch.eq.1 .and. idpdp .ne. 0) then
            nbd = 4*nsbb
         else if(irdof.ge.0) then
            nbd = nsbb*mdof**2
         else if(irdof.eq.-3) then
            nbd = nsbb*1**2
         else if(irdof.eq.-6) then
            nbd = nsbb*2**2
         else
            nbd = (nelm(neqp1)-neqp1)*mdof**2
         end if
      end if
      if(.not.compute_flow .and. allocated(nop)) then
         if (iccen .eq. 1 .and. igauss .gt. 1) then
! Do nothing, we need the nop array
         else
            deallocate(nop)
         end if
      end if
      if(irun.eq.1) then
         if (compute_flow .or. iccen .eq. 1) then
            if (ice .eq.0) then
               allocate(bp(max((mdof*n0),(idof*n0),2*n0)))
               bp = 0.0
            else
               allocate(bp(max((mdof*n0),((idof+1)*n0),2*n0)))
               bp = 0.0
            end if
c
c  make sure size for stress solution is at least 6*neq
c  for CO2 app (hopefully less for sequentially coupled)
c
c             if (istrs.ne.0) then
c             if(icnl.ne.0.and.idof_co2.ge.2) idofdum = 5
c             if(icnl.eq.0.and.idof_co2.ge.2) idofdum = 6
c             if(icnl.eq.0.and.ico2.eq.0) idofdum = 5
c             if(icnl.eq.0.and.ico2.ge.1) idofdum = 6
c             if(idof_stress.le.3)then
c              idofdum = max(idof_stress,idof_co2)
c             endif
c             if (max((mdof*n0),2*n0).lt.idofdum*n0) then
c	        deallocate(bp)
c	        allocate(bp(n0*idofdum))
c             endif
c             idof = idofdum
c	      endif
c
c            bp = 0
c

         end if
      end if
c     New allocation of bigblock
      if(compute_flow .or. iccen .eq. 1) then
         if(irdof.ne.13) then
            nstorederiv = 40*n0 + n7a
         else
c there are 15 arrays we can eliminate for saturated only
            nstorederiv = 25*n0 + 15*1 + n7a
         endif
         nstorepiv = n0*mdof**2
         nstoresolve = nbd + kgmres + nstorepiv
         icall = -1
      else
c    try only minimal memory
c         nstorederiv = 39*n0 + 3*n7 + n7a
         nstorederiv = 39*1 + 3*1 + 1
         nstorepiv = 1
         nstoresolve = 1
         icall = -2
      end if
      nbigblock = max(nstorederiv,nstoresolve)
      call storage_derivatives(icall,irun)
c     write out new size for b matrix
      if(compute_flow .or. iccen .eq. 1) then
         if (iout .ne. 0) write(iout, 6010) nbd
         if (iatty .ne. 0) write(iatty, 6010) nbd
      else
         if (iout .ne. 0) write(iout, 6010) 1
         if (iatty .ne. 0) write(iatty, 6010)
      end if
 6010 format(1x,'storage available for b matrix resized to ',
     &     i10, '<<<<<<')

c     Initialize eos values as necessary (for air)
c     first set phase state for wtsi
c      call wtsictr(-1)
      if (ico2.lt.0.and.ice.eq.0) then
         call airctr(-1,0)
!     else if (ico2.lt.0.and.ice.ne.0) then
!        call icectr(-1,0)
      endif
c     Calculate rho1grav
      if(ico2.lt.0) then
         rho1grav = crl(1,1)*(9.81d-6)
      else
         rho1grav = rol0*9.81d-6
      endif
c      rho1grav = 997.*9.81d-6

c     Moved calls to airctr for head option to below the iflg = -1 call
c if head input has been used,convert to pressures
c
c      if(ihead.ne.0.or.ichead.ne.0) then
c gaz 083119 this could possibly cause problems if ichead does more that I think  <<<<
       if(ihead.ne.0) then
         if(head0.gt.0) then
            phi_inc = head0*crl(1,1)*(-grav)
         end if

c convert from head to pressure if no initial value file
         if(iread.le.0 .or. .not. pres_read) then
c
c this is new(set IC heads equal to lowest fixed head
c
c            call wtsictr(12)
c
            call airctr(7,0)
         else
c now add pressure inc. equal to the head inc.
c in reading pressures from intial value field
            do i= 1,n
               pho(i) = pho(i) + phi_inc
            enddo
         endif
c convert from head to pressure for boundary nodes
         call airctr(9,0)
      endif

      if (ichead.ne.0 .and. head0 .gt. 0.) then
         phi_inc = head0*rol0*(-grav)
      end if

c**** set array ordering if necessary ****
      if(compute_flow .or. iccen .eq. 1) then
       if(icarb.eq.0) then
         call setord
       else
         idof_save = idof
         idof = 3
         call setord
         idof = idof_save
       endif
      end if
c**** initialize porosities in subsidence calcs ****
      call subsidence(-2)
c**** initialize coeffients adjust volumes in dual porosity calcs ****
      call dual (2)
c gaz 101518  debug only
c      to(1) = 85.
c      to(20) =5.
c**** initialize coefficients adjust volumes in dpdp calcs ****
      call dpdp (1)
      do i = 1, n
         volume(i) = sx1(i)
c gaz 01-18-2011 change code to allow for negative temperatures
c         if (to (i) .le. zero_t)   to (i) = tin0
         if(irdof.ne.13) then
            if (pho(i) .le. zero_t)   pho (i) = pein
         end if
         phi (i) = pho(i)
         t (i) = to(i)
         vtot = vtot + volume(i)
      end do

      if (mass_read) then
         iconv_tmp = iconv
         iconv = 1
         call convctr(3)
         deallocate (mass_var)
         iconv = iconv_tmp
      end if
c initialize pressures and temperatures
      if(ipini.eq.0.and.istrs.ne.0) then
         phini = pho
         psini = ps
         tini = to
      else if(iread.eq.0) then
         phini = pho
         psini = ps
         tini = to
      else
c gaz debug 110314 do nothing except for cases above
c if requested psini read from salt restart
c phini and tini have be read in from disk, just set porosity
c gaz debug 110314 psini remains equal to input file values
c except for salt restart file which reads in psini
         if(isalt.eq.0)  psini = ps
         phini = pho
         tini = to
      endif
c initialize porosity model -5 if necessary
      if(iporos.ne.0) call porosi(5)
c apply distributed forces
      call stressctr(2,0)
c this call handled in call to porosi above
c
c  find permeability nodes and startup operations(allocate memory)
c
      call stressctr(16,0)
c
c   combine multiple generalzed head BCs
c
      call inflogh(1)
c
c
c if boundary node fix initial value to fixed value
c
      if(ico2.lt.0.and.iread.le.0.and.ice.eq.0) then
         do i = 1, n
            if(ka(i).eq.-1.or.ka(i).eq.-2) then
               pho(i) = pflow(i)
               phi (i) = pho(i)
            endif
         end do
      endif
c  Set uz correction for water table
c
	call uz_wtctr(100)
c
c     These arrays are present only if the flow field is
c     being computed
      if(compute_flow) then
         do i = 1, n
            if (ps (i) .le. zero_t)   deni (i) = 0.0
            if(irdof.ne.13) then
               if (ps (i) .le. zero_t)   dstm (i) = 0.0
            endif
         end do
      end if

c modify volumes and heat capacities to affect boundary conditions
      call bcon(1)

c
c**** initialize ngas varibles ****
      call co2ctr (6)

c bookeeping for air-water or methane
      if(ico2.lt.0.and.ice.eq.0) then
c get cell lengths for wtsi if necessary
         call airctr(11, 0)
         call airctr(6, 0)
c check temperature used in isothermal calculations
         call airctr(13, 0)
      else if(ico2.lt.0.and.ice.ne.0) then
         call icectr(6, 0)
      endif
      if(icarb.eq.1) then
         call icectrco2(14,0)
         call icectrco2(6,0)
      endif

c change porosity and permeability if necessary to
c establish initial values
      if(iporos.ne.0) call porosi(3)

c**** determine initial coefficients for thermo fits ****
      call coeffc
c
c gaz debug 090113
c**** insure anl = anlo for cden to get density right for firsat timestep****
       call concen(6,0)
c**** determine initial variable state ****
c gaz 10-18-2001     call sice (1)
      if(compute_flow .or. iccen .eq. 1) then
         if(ice.eq.0) then
            if(icarb.eq.1) then
               call icectrco2(-1,0)
               call icectrco2(1,0)
               call icectrco2(-34,0)
               call icectrco2(3,0)
               call icectrco2(-3,0)
               call icectrco2(-33,0)
               if (ishisrlp .ne. 0) call check_rlp_carb
               call icectrco2(-35,0)
            else
               call varchk (0, 0)
            end if
         else
            call icectr (-1,0)
            call icectr (1,0)
c added check for hydrate line
            call icectr (-6,0)
c allocate space check properties, dellocate space
            if(idof_meth.ne.7) then
               call icectr (-34,0)
               call icectr (3,0)
               call icectr (-3,0)
               call icectr (-33,0)
               call icectr (-35,0)
	    else
	      call icectr (-34,0)
c     id mobile methane and water
              call hydrate_equil(2,0)
              call icectr (3,0)
              call icectr (-3,0)
              call icectr (-33,0)
              call  hydrate_equil(3,0)
              call icectr (-35,0)

           endif
        endif
      else
         if(ico2.lt.0.and.ice.eq.0) then
c     air water problem
            do i = 1, n0
               rolf(i)=crl(1,1)*
     2              (1.0+crl(3,1)*(phi(i)-crl(4,1)))
            end do
         else if(ico2.lt.0.and.ice.ne.0) then
c gaz 10-18-2001 following used ase dummy parameters
c rnwn1, rnwn2, rnwn3, rnwn, rnwd1, rnwd2, rnwd3, rnwd
            call methane_properties(2,2,phi(i),t(i),rnwn1, rnwn2,
     &           rolf(i),rnwn3, rnwn, rnwd1, rnwd2)
         else
c     nonisothermal problem
            do i = 1, n0
               call water_density(t(i),phi(i),rolf(i))
            end do
         end if
      end if
      call dual (1)
      if(compute_flow) then
         call dpdp ( 2 )
      end if
c
c calculate areas for use in flow through boundary conditions
c gaz 12222001
c      call area_length_calc(1)
c
      if(compute_flow ) then
c**** initialize some variables ****
         am0 = 0.0
         ame = 0.0
         astmo = 0.0
c
c change volumes so mass and energy comes out correct (GAZ 1-29-09)
c
      call vboun(1,0)
c
c**** calculate initial mass and energy ****
         do i = 1, n
            if(irdof.ne.13) then
               deneh(i) = denei(i) * dtot
               ame = ame + deneh(i) * volume(i)
               astmo = astmo + dstm (i)
               denei(i) = 0.0
            endif
	    if (ifree .ne. 0) then
               so(i) = rlxyf(i)
            else if (irdof.ne.13) then
               so (i) = s(i)
            end if
            denh (i) = deni (i) * dtot
            am0 = am0 + denh (i) * volume(i)
            deni (i) = 0.0
            nopt (i) = 0
         end do
         if(ice.ne.0) then
c**** calculate initial mass and energy for component system ****
            call icectr(7,0)
         elseif(icarb.eq.1) then
            call icectrco2(7,0)
            call icectrco2(9,0)
         endif
         amass = am0
         aener = ame
         asteam = astmo
      end if
c
c**** initialize noncondensible and ice arrays if applicable ****
      call co2ctr (1)
c added call 122314 split funcionality
      call co2ctr (-1)
      dtotdm = dtot
      if (icgts .gt. 0)  then
         do nicg = 1, icgts
            if (dit(nicg) .gt. days)  then
               ditnd = dit(nicg)
               go to 81
            end if
         end do
         nicg = 1
         icgts = 0
         ditnd = 1.0d+30
      else
         ditnd = 1.0d+30
      end if
 81   continue
c gaz 05-08-2009 moved here so stress initcalcs show up at time = 0
c calculate initial stress field and displacements
c
         call stress_perm(-3,0)
         call stress_uncoupled(1)
c
c reset boundary conditions for principal stresses (fraction of lithostatic)
c
         call stressctr(3,0)
c
c
c set neq_primary back to true primary for gdkm
c
      if(gdkm_flag.ne.0) then
       neq_primary = neq_gdkm
      endif
c
!      if(contim.ge.0) then
         call contr (-1)
         call contr ( 1)
c gaz debug 021114
c output map of active nodes
         call active_nodes_ctr(3)
!      else
!         call contr_days (-1)
!         call contr_days ( 1)
!      endif
c
c**** initalize concentration ****
      if (iccen .ne. 0)  call concen (-1,0)

      call pest(1)

! Moved history file setup below zone volume calculations
!      if (hist_flag) then
!         call plot_new (0, 0.0d00, 0.0d00, 0.0d00, 0.0d00)
!      else
!         call plot (0, 0.0d00, 0.0d00)
!      end if

      tajj = tyming(caz) - tajj

      if (iout .ne. 0) write(iout, 6020)  tajj
      if (iptty .gt. 0)  write(iptty, 6020)  tajj
 6020 format(//, 1x, 'time for reading input, forming coefficients ',
     *     g10.3, //)

      tasii = tyming(caz)

c**** initialize pflow if needed ****
      do i = 1, n
         if (abs(pflow(i)) .lt. zero_t.or.
     &        pflow(i).eq.-999.)  pflow(i) = pho(i)
         if (icarb.ne.0) then
            if(abs(pflowco2(i)) .lt. zero_t.or.
     &           pflowco2(i).eq.-999.)  pflowco2(i) = phoco2(i)
         endif
         if(ico2.ge.0.or.ice.ne.0) then
            if (qflux(i) .eq. -999 .and. qflxm(i).gt.0.0)
     &           qflux(i) = to(i)
         endif
      end do
c
c call calculation of areas,weights,etc for use in boundary
c conditions
c
      call area_vol_weightctr(1)
c
c apply areas, weights etc.
c
      call area_vol_weightctr(2)
c
c modify area/d and volumes from read-in data
c
      call coef_replace_ctr(1)
      call coef_replace_ctr(2)
c correct air pressures that are less than pref
c correct appropriate boundary pressures and change to seepage face
c conditions if necessary
c
      if(ihead.ne.0) call airctr(10,0)
c
      daynew = day
c release memory for original coordinates
c     BAR - leave around since multiple simulations is a possibility
c      call rarng(1)

      if(sptrak .and. .not. compute_flow) ihf = 1

! zvd - Compute zone volumes for average output
      if (out_zones) then
         allocate (zone_volume(node_azones))
         zone_volume = 0.
         node_ptr_num => node_head_num
         node_ptr => node_head
         num_zone = 0
         outer: do
            if (.not. associated(node_ptr_num)) exit
            num_zone = num_zone + 1
! Number of nodes in current zone
            j = node_ptr_num%node_number
            count = 0
            if (j .ne. 0) then
               sum_values: do
                  if (.not. associated(node_ptr)) exit
! Number of current node in the zone
                  i = node_ptr%node_number
                  zone_volume(num_zone) = sx1(i)*pnx(i) +
     &                 zone_volume(num_zone)
                  node_ptr =>node_ptr%nnp
                  count = count + 1
                  if (count .eq. j) exit
               end do sum_values
            else
               zone_volume(num_zone) = 0
            end if
            node_ptr_num =>node_ptr_num%nnp
         end do outer
      end if

c call gdpm_corr to add connections if necessary
c
       call gdpm_corr(-1)

! zvd - moved history file setup after zone volume calculations
      if (hist_flag) then
! Initial call for file setup
         call plot_new (0, 0.0d00, 0.0d00, 0.0d00, 0.0d00)
      else
         call plot (0, 0.0d00, 0.0d00)
      end if

! zvd - moved initial call here to accomodate new output options
! for trac macro
      if (.not. ptrak .and. .not. sptrak) call plotc1(0,0)

! zvd - Oct 7, 2004 Need to set sk values after we are sure all
! modifications have been made to nelmdg
      if (.not. compute_flow) then
         if (flux_flag(1:3) .eq. 'all' .or. flux_flag(1:3) .eq.
     &        'liq') then
            do i = 1, neq
               sk(i) = a_axy(nelmdg(i) - neq - 1)
            end do
         end if
         if (flux_flag(1:3) .eq. 'all' .or. flux_flag(1:3) .eq.
     &        'vap') then
            do i = 1, neq
               if (ico2 .gt. 0 ) then
                  qc(i) = a_vxy(nelmdg(i) - neq - 1)
               else if (ico2 .lt. 0) then
                  qh(i) = a_vxy(nelmdg(i) - neq - 1)
               end if
            end do
         end if
      end if
c
! zvd - 07-Oct-2010 Add optional printout of rlp table
      if (ishisrlp .ne. 0) then
         if (icarb .eq. 0) then
            call check_rlp
         end if
      end if

      deallocate(idum1,idum2)

      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


program test
use activations
implicit none

call test_cost()

contains

subroutine test_cost()
	double precision, dimension(5) :: out
	double precision, dimension(5) :: target
	double precision :: cost_output

	out = (/ 1.0, 2.0, 3.0, 4.0, 5.0/)
	target = (/2.0, 3.0, 4.0, 5.0, 6.0/)

	cost_output = cost(out, target)
	if(cost_output == 2.50) then
		write(*,*) "test_cost: success"
	else
		write(*,*) "test_cost: failure"
	end if

end subroutine test_cost

end program test


text:
module define_bc_module

  use bl_types
  use ml_layout_module
  use bc_module
  use inhomogeneous_bc_val_module

  implicit none

  type bc_level

     integer, pointer :: phys_bc_level_array(:,:,:) => Null()
     integer, pointer ::  adv_bc_level_array(:,:,:,:) => Null()
     integer, pointer ::  ell_bc_level_array(:,:,:,:) => Null()

  end type bc_level

  type bc_tower

     integer :: max_level_built = 0
     type(bc_level), pointer :: bc_tower_array(:) => Null()
     integer       , pointer :: domain_bc(:,:) => Null()

  end type bc_tower

  private

  public :: bc_level, bc_tower, &
            initialize_bc, bc_tower_init, bc_tower_level_build, bc_tower_destroy

contains

  subroutine initialize_bc(the_bc_tower,num_levs,dm,pmask,num_scal_bc_in,num_tran_bc_in)

     use bc_module

     use probin_common_module, only : bc_lo, bc_hi

     type(bc_tower), intent(  out) :: the_bc_tower
     integer       , intent(in   ) :: num_levs,dm,num_scal_bc_in,num_tran_bc_in
     logical       , intent(in   ) :: pmask(:)

     integer :: domain_phys_bc(dm,2)

     num_scal_bc  = num_scal_bc_in
     num_tran_bc  = num_tran_bc_in

     vel_bc_comp  = 1
     pres_bc_comp = vel_bc_comp+dm
     scal_bc_comp = pres_bc_comp+1
     tran_bc_comp = scal_bc_comp+num_scal_bc

     ! Define the physical boundary conditions on the domain
     ! Put the bc values from the inputs file into domain_phys_bc
     domain_phys_bc(1,1) = bc_lo(1)
     domain_phys_bc(1,2) = bc_hi(1)
     if (pmask(1)) then
        domain_phys_bc(1,:) = BC_PER
        if (bc_lo(1) .ne. -1 .or. bc_hi(1) .ne. -1) &
             call bl_error('MUST HAVE BCX = -1 if PMASK = T')
     end if
     if (dm > 1) then
        domain_phys_bc(2,1) = bc_lo(2)
        domain_phys_bc(2,2) = bc_hi(2)
        if (pmask(2)) then
           domain_phys_bc(2,:) = BC_PER
           if (bc_lo(2) .ne. -1 .or. bc_hi(2) .ne. -1) &
                call bl_error('MUST HAVE BCY = -1 if PMASK = T')
        end if
     end if
     if (dm > 2) then
        domain_phys_bc(3,1) = bc_lo(3)
        domain_phys_bc(3,2) = bc_hi(3)
        if (pmask(3)) then
           domain_phys_bc(3,:) = BC_PER
           if (bc_lo(3) .ne. -1 .or. bc_hi(3) .ne. -1) &
                call bl_error('MUST HAVE BCZ = -1 if PMASK = T')
        end if
     end if

     ! Initialize the_bc_tower object.
     call bc_tower_init(the_bc_tower,num_levs,dm,domain_phys_bc)

  end subroutine initialize_bc

  subroutine bc_tower_init(bct,num_levs,dm,phys_bc_in)

    type(bc_tower ), intent(  out) :: bct
    integer        , intent(in   ) :: num_levs
    integer        , intent(in   ) :: dm
    integer        , intent(in   ) :: phys_bc_in(:,:)

    allocate(bct%bc_tower_array(num_levs))
    allocate(bct%domain_bc(dm,2))

    bct%domain_bc(:,:) = phys_bc_in(:,:)

  end subroutine bc_tower_init

  subroutine bc_tower_level_build(bct,n,la)

    type(bc_tower ), intent(inout) :: bct
    integer        , intent(in   ) :: n
    type(layout)   , intent(in   ) :: la

    integer :: ngrids,dm
    integer :: default_value

    if (associated(bct%bc_tower_array(n)%phys_bc_level_array)) then
      deallocate(bct%bc_tower_array(n)%phys_bc_level_array)
      deallocate(bct%bc_tower_array(n)%adv_bc_level_array)
      deallocate(bct%bc_tower_array(n)%ell_bc_level_array)
      bct%bc_tower_array(n)%phys_bc_level_array => NULL()
      bct%bc_tower_array(n)%adv_bc_level_array => NULL()
      bct%bc_tower_array(n)%ell_bc_level_array => NULL()
    end if

    ngrids = layout_nlocal(la)
    dm = layout_dim(la)

    allocate(bct%bc_tower_array(n)%phys_bc_level_array(0:ngrids,dm,2))
    default_value = INTERIOR
    call phys_bc_level_build(bct%bc_tower_array(n)%phys_bc_level_array,la, &
                             bct%domain_bc,default_value)

    ! Here we allocate (dm+1)=pres_bc_comp components for x_u and x_p
    !                  num_scal_bc components for scalars
    !                  num_tran_bc components for transport coefficients
    allocate(bct%bc_tower_array(n)%adv_bc_level_array(0:ngrids,dm,2,pres_bc_comp+num_scal_bc+num_tran_bc))
    default_value = INTERIOR
    call adv_bc_level_build(bct%bc_tower_array(n)%adv_bc_level_array, &
                            bct%bc_tower_array(n)%phys_bc_level_array,default_value)

    ! Here we allocate (dm+1)=pres_bc_comp components for x_u and x_p
    !                  num_scal_bc components for scalars
    !                  num_tran_bc components for transport coefficients
    allocate(bct%bc_tower_array(n)%ell_bc_level_array(0:ngrids,dm,2,pres_bc_comp+num_scal_bc+num_tran_bc))
    default_value = BC_INT
    call ell_bc_level_build(bct%bc_tower_array(n)%ell_bc_level_array, &
                            bct%bc_tower_array(n)%phys_bc_level_array,default_value)

     bct%max_level_built = n

  end subroutine bc_tower_level_build

  subroutine bc_tower_destroy(bct)

    type(bc_tower), intent(inout) :: bct

    integer :: n

    do n = 1,bct%max_level_built
       deallocate(bct%bc_tower_array(n)%phys_bc_level_array)
       deallocate(bct%bc_tower_array(n)%adv_bc_level_array)
       deallocate(bct%bc_tower_array(n)%ell_bc_level_array)
       bct%bc_tower_array(n)%phys_bc_level_array => NULL()
       bct%bc_tower_array(n)%adv_bc_level_array => NULL()
       bct%bc_tower_array(n)%ell_bc_level_array => NULL()
    end do
    deallocate(bct%bc_tower_array)

    deallocate(bct%domain_bc)

  end subroutine bc_tower_destroy

  subroutine phys_bc_level_build(phys_bc_level,la_level,domain_bc,default_value)

    integer     , intent(inout) :: phys_bc_level(0:,:,:)
    integer     , intent(in   ) :: domain_bc(:,:)
    type(layout), intent(in   ) :: la_level
    integer     , intent(in   ) :: default_value
    type(box) :: bx,pd
    integer :: d,i

    pd = layout_get_pd(la_level)

    phys_bc_level = default_value

    i = 0
    do d = 1,layout_dim(la_level)
       phys_bc_level(i,d,1) = domain_bc(d,1)
       phys_bc_level(i,d,2) = domain_bc(d,2)
    end do

    do i = 1,layout_nlocal(la_level)
       bx = layout_get_box(la_level,global_index(la_level,i))
       do d = 1,layout_dim(la_level)
          if (lwb(bx,d) == lwb(pd,d)) phys_bc_level(i,d,1) = domain_bc(d,1)
          if (upb(bx,d) == upb(pd,d)) phys_bc_level(i,d,2) = domain_bc(d,2)
       end do
    end do

  end subroutine phys_bc_level_build

  subroutine adv_bc_level_build(adv_bc_level,phys_bc_level,default_value)

    integer  , intent(inout) ::  adv_bc_level(0:,:,:,:)
    integer  , intent(in   ) :: phys_bc_level(0:,:,:)
    integer  , intent(in   ) :: default_value

    integer :: dm,igrid,d,lohi

    ! these boundary conditions are written for a stokes system where we
    ! solve for the INCREMENT to velocity and pressure

    adv_bc_level = default_value

    dm = size(adv_bc_level,dim=2)

    do igrid = 0, size(adv_bc_level,dim=1)-1
    do d = 1, dm
    do lohi = 1, 2

       if (phys_bc_level(igrid,d,lohi) == PERIODIC .or. &
           phys_bc_level(igrid,d,lohi) == INTERIOR ) then

          ! retain the default value of INTERIOR

       else if ( (phys_bc_level(igrid,d,lohi) >= NO_SLIP_START) .and.  (phys_bc_level(igrid,d,lohi) <= NO_SLIP_END) ) then

          ! for normal velocity we impose a Dirichlet velocity condition
          ! for transverse velocity we imposie a Dirichlet velocity condition
          ! for pressure we use extrapolation
          adv_bc_level(igrid,d,lohi,1:dm)                    = DIR_VEL  ! normal and transverse velocity
          adv_bc_level(igrid,d,lohi,pres_bc_comp)                    = FOEXTRAP ! pressure

       else if ( (phys_bc_level(igrid,d,lohi) >= SLIP_START) .and.  (phys_bc_level(igrid,d,lohi) <= SLIP_END) ) then

          ! for normal velocity we impose a Dirichlet velocity condition
          ! for transverse velocity we imposie a Dirichlet stress condition
          ! for pressure we use extrapolation
          adv_bc_level(igrid,d,lohi,vel_bc_comp:vel_bc_comp+dm-1) = DIR_TRACT ! transverse velocity
          adv_bc_level(igrid,d,lohi,vel_bc_comp+d-1)              = DIR_VEL   ! normal velocity
          adv_bc_level(igrid,d,lohi,pres_bc_comp)                 = FOEXTRAP  ! pressure

       else

          print*,'adv_bc_level_build',igrid,d,lohi,phys_bc_level(igrid,d,lohi)
          call bl_error('BC TYPE NOT SUPPORTED 1')

       end if

       ! The scalars can be handled in many different ways, so defer to application-specific code:
       call scalar_bc(phys_bc_level(igrid,d,lohi),adv_bc_level(igrid,d,lohi, &
                      scal_bc_comp:scal_bc_comp+num_scal_bc-1))

       ! The transport coefficients can be handled in many different ways, so defer to application-specific code:
       call transport_bc(phys_bc_level(igrid,d,lohi),adv_bc_level(igrid,d,lohi, &
                         tran_bc_comp:tran_bc_comp+num_tran_bc-1))

       if (any(adv_bc_level(igrid,d,lohi,scal_bc_comp:scal_bc_comp+num_scal_bc-1) .eq. -999)) then

          print*,'adv_bc_level_build',igrid,d,lohi,phys_bc_level(igrid,d,lohi)
          call bl_error('BC TYPE NOT SUPPORTED 2')

       end if

    end do
    end do
    end do

  end subroutine adv_bc_level_build

  subroutine ell_bc_level_build(ell_bc_level,phys_bc_level,default_value)

    integer  , intent(inout) ::  ell_bc_level(0:,:,:,:)
    integer  , intent(in   ) :: phys_bc_level(0:,:,:)
    integer  , intent(in   ) :: default_value

    integer :: dm
    integer :: igrid,d,lohi

    ell_bc_level = default_value

    dm = size(ell_bc_level,dim=2)

    do igrid = 0, size(ell_bc_level,dim=1)-1
    do d = 1, dm
    do lohi = 1, 2

       if (phys_bc_level(igrid,d,lohi) == INTERIOR) then

          ! retain the default value of INTERIOR

       else if (phys_bc_level(igrid,d,lohi) == PERIODIC) then

          ! pressure is periodic
          ell_bc_level(igrid,d,lohi,pres_bc_comp) = BC_PER

       else

          ! pressure and temperature are homogeneous neumann
          ell_bc_level(igrid,d,lohi,pres_bc_comp) = BC_NEU

       end if

    end do
    end do
    end do

  end subroutine ell_bc_level_build

end module define_bc_module

max_stars_repo_name: BoxLib-Codes/LowMachMixtures

max_issues_repo_name: BoxLib-Codes/LowMachMixtures

max_forks_repo_name: BoxLib-Codes/LowMachMixtures


text:
! Initialize some arrays from vars module:

subroutine init()

  use vars
  use tracers

  implicit none

  fzero = 0.

  ttend = 0.
  qtend = 0.
  wsub = 0.
  unudge = 0.
  vnudge = 0.
  tnudge = 0.
  qnudge = 0.
  qlsvadv = 0.
  tlsvadv = 0.
  ulsvadv = 0.
  vlsvadv = 0.
  qstor = 0.
  tstor = 0.
  ustor = 0.
  vstor = 0.

  radlwup = 0.
  radlwdn = 0.
  radswup = 0.
  radswdn = 0.
  radqrlw = 0.
  radqrsw = 0.

  tlat = 0.
  tlatqi = 0.
  tadv = 0.
  tdiff = 0.
  qifall = 0.
  qpfall = 0.

  trwle = 0.
  trwsb = 0.
  tradv = 0.
  trdiff = 0.
  trphys = 0.

  gamt0 = 0.
  gamq0 = 0.

end subroutine init

max_stars_repo_name: kbren/uwnet

max_issues_repo_name: kbren/uwnet

max_forks_repo_name: kbren/uwnet


text:
subroutine print_cvector(A,name)
complex(kind=8),  intent(in) :: A(:)
character(len=*), intent(in) :: name
character(Len=1) :: tab = char(9)
integer :: i, n
n = size(A)
write(*,'(a,a)') name," = "
do i=1,n
    write(*,"("//FMT//",a,"//FMT//",a)") real(A(i),8), tab, aimag(A(i)), "i"
enddo
write(*,*) ""
end subroutine print_cvector

max_stars_repo_name: buvoli/epbm

max_issues_repo_name: buvoli/epbm

max_forks_repo_name: buvoli/epbm


text:
Brianne Crabtree is a 5th year WFCB major / Geology minor and a Damento Jugglers Damento at large.  She can be found juggling twice a week in Davis and Sacramento.  Bri also likes unicycling and other circuslike tricks.

20051025 00:59:45 nbsp I was wondering what Damento Jugglers that move was called :) Users/CraigBrozinsky

20060207 20:20:23 nbsp I went over to Everson tonight, and there was no one there!  Hopefully Ill see you guys next week. Users/AllisonEriksen

20060320 19:13:38 nbsp Gotem Users/BradRodriguez

20060326 20:42:33 nbsp When/where does juggling start again?  Ive been practicing like crazy, but seem to have hit the shallow part of the learning curve... Users/AllisonEriksen


text:
!> @ingroup quick_sort_recursive_dec
!>
!> @brief
!> Run a recursive three way quick sort on a vector of double precision values.
!>
!> @details
!> Run a decending on a vector of double precision values to sort them in the specified order.
!>
!> @param vec
!> The vector double precision values to be sorted.
!>
!> @author Arin R. Bratt
!> @date 11/07/2014
pure subroutine quick_sort_recursive_dec_dble(vec)
  ! Arguments
  double precision, intent(inout) :: vec(:)

  ! Variables
  integer :: low
  integer :: high
  integer :: low_eq
  integer :: high_eq
  double precision :: pivot
  double precision :: temp

  ! Go to interchange sort if list is small
  if (size(vec,1) <= max_interchange_sort_size) then
     call interchange_sort_dec_dble(vec)
  else
     ! Initilize values
     low = 1
     high = size(vec,1)
     low_eq = 1
     high_eq = size(vec,1)

     ! Build pivot
     pivot = (vec(low) + vec(high)) / 2
     if ( vec(low) < vec(high)) then
        call swap(vec(low),vec(high))
     end if

     ! Main loop with three way partition
     do
        ! Find pivot points
        do
           low = low + 1
           if (.not. pivot < vec(low)) exit
        end do
        do
           high = high - 1
           if (.not. vec(high) < pivot .or. high == 1) exit
        end do

        ! Exit if they are on the wrong side of the pivot
        if (low >= high) exit

        ! Run main quick sort swap
        call swap(vec(low),vec(high))

        ! Check for if values are equal to pivot, if so swap
        if (.not. vec(high) < pivot) then
           call swap(vec(high),vec(high_eq))
           high_eq = high_eq - 1
        end if

        if (.not. pivot < vec(low)) then
           call swap(vec(low),vec(low_eq))
           low_eq = low_eq + 1
        end if
     end do
     high = low - 1

     ! Move equals to the middle
     do while(1 < low_eq)
        low_eq = low_eq - 1
        call swap(vec(low_eq),vec(high))
        high = high - 1
     end do

     do while(high_eq < size(vec,1))
        high_eq = high_eq + 1
        call swap(vec(high_eq),vec(low))
        low = low + 1
     end do

     ! Recursive calls
     if (size(vec,1) - low > high - 1) then
        if (high > 1) call quick_sort_recursive_dec_dble(vec(:high))
        if (low < size(vec,1)) call quick_sort_recursive_dec_dble(vec(low:))
     else
        if (low < size(vec,1)) call quick_sort_recursive_dec_dble(vec(low:))
        if (high > 1) call quick_sort_recursive_dec_dble(vec(:high))
     end if
  end if
end subroutine quick_sort_recursive_dec_dble

max_stars_repo_name: arinrb/bsort

max_issues_repo_name: arinrb/bsort

max_forks_repo_name: arinrb/bsort


text:
PROGRAM bounds
IMPLICIT NONE
  INTEGER :: i
  REAL, DIMENSION(5) :: a = 0

  CALL sub1(a, 5, 6)

  DO i = 1, 6
    WRITE(*, 100) i, a(i)
    100 FORMAT ('A(',I1, ') = ', F6.2)
  END DO

END PROGRAM bounds

SUBROUTINE sub1(a, ndim, n)
IMPLICIT NONE
  INTEGER, INTENT(IN) :: ndim
  REAL, INTENT(OUT), DIMENSION(ndim) :: a
  INTEGER, INTENT(IN) :: n
  INTEGER :: i

  DO i = 1, n
    a(i) = i
  END DO

END SUBROUTINE sub1

max_stars_repo_name: evanmacbride/fortran-practice

max_issues_repo_name: evanmacbride/fortran-practice

max_forks_repo_name: evanmacbride/fortran-practice