text:
*DECK CG
      SUBROUTINE CG (NM, N, AR, AI, WR, WI, MATZ, ZR, ZI, FV1, FV2, FV3,
     +   IERR)
C***BEGIN PROLOGUE  CG
C***PURPOSE  Compute the eigenvalues and, optionally, the eigenvectors
C            of a complex general matrix.
C***LIBRARY   SLATEC (EISPACK)
C***CATEGORY  D4A4
C***TYPE      COMPLEX (RG-S, CG-C)
C***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
C***AUTHOR  Smith, B. T., et al.
C***DESCRIPTION
C
C     This subroutine calls the recommended sequence of
C     subroutines from the eigensystem subroutine package (EISPACK)
C     to find the eigenvalues and eigenvectors (if desired)
C     of a COMPLEX GENERAL matrix.
C
C     On INPUT
C
C        NM must be set to the row dimension of the two-dimensional
C          array parameters, AR, AI, ZR and ZI, as declared in the
C          calling program dimension statement.  NM is an INTEGER
C          variable.
C
C        N is the order of the matrix A=(AR,AI).  N is an INTEGER
C          variable.  N must be less than or equal to NM.
C
C        AR and AI contain the real and imaginary parts, respectively,
C          of the complex general matrix.  AR and AI are two-dimensional
C          REAL arrays, dimensioned AR(NM,N) and AI(NM,N).
C
C        MATZ is an INTEGER variable set equal to zero if only
C          eigenvalues are desired.  Otherwise, it is set to any
C          non-zero integer for both eigenvalues and eigenvectors.
C
C     On OUTPUT
C
C        WR and WI contain the real and imaginary parts, respectively,
C          of the eigenvalues.  WR and WI are one-dimensional REAL
C          arrays, dimensioned WR(N) and WI(N).
C
C        ZR and ZI contain the real and imaginary parts, respectively,
C          of the eigenvectors if MATZ is not zero.  ZR and ZI are
C          two-dimensional REAL arrays, dimensioned ZR(NM,N) and
C          ZI(NM,N).
C
C        IERR is an INTEGER flag set to
C          Zero       for normal return,
C          10*N       if N is greater than NM,
C          J          if the J-th eigenvalue has not been
C                     determined after a total of 30 iterations.
C                     The eigenvalues should be correct for indices
C                     IERR+1, IERR+2, ..., N, but no eigenvectors are
C                     computed.
C
C        FV1, FV2, and FV3 are one-dimensional REAL arrays used for
C          temporary storage, dimensioned FV1(N), FV2(N), and FV3(N).
C
C     Questions and comments should be directed to B. S. Garbow,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C     ------------------------------------------------------------------
C
C***REFERENCES  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
C                 Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
C                 system Routines - EISPACK Guide, Springer-Verlag,
C                 1976.
C***ROUTINES CALLED  CBABK2, CBAL, COMQR, COMQR2, CORTH
C***REVISION HISTORY  (YYMMDD)
C   760101  DATE WRITTEN
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  CG
C
      INTEGER N,NM,IS1,IS2,IERR,MATZ
      REAL AR(NM,*),AI(NM,*),WR(*),WI(*),ZR(NM,*),ZI(NM,*)
      REAL FV1(*),FV2(*),FV3(*)
C
C***FIRST EXECUTABLE STATEMENT  CG
      IF (N .LE. NM) GO TO 10
      IERR = 10 * N
      GO TO 50
C
   10 CALL  CBAL(NM,N,AR,AI,IS1,IS2,FV1)
      CALL  CORTH(NM,N,IS1,IS2,AR,AI,FV2,FV3)
      IF (MATZ .NE. 0) GO TO 20
C     .......... FIND EIGENVALUES ONLY ..........
      CALL  COMQR(NM,N,IS1,IS2,AR,AI,WR,WI,IERR)
      GO TO 50
C     .......... FIND BOTH EIGENVALUES AND EIGENVECTORS ..........
   20 CALL  COMQR2(NM,N,IS1,IS2,FV2,FV3,AR,AI,WR,WI,ZR,ZI,IERR)
      IF (IERR .NE. 0) GO TO 50
      CALL  CBABK2(NM,N,IS1,IS2,FV1,N,ZR,ZI)
   50 RETURN
      END


program integer_kinds
    implicit none
    integer :: range
    integer :: kind

    do range = 0, 45
        kind = selected_int_kind(range)
        print *, range, kind
    end do
end program integer_kinds


text:
program main
    use std ! standard package of plantFEM library
    implicit none

    ! This example utilizes Array Class, where
    ! You can handle external files.

    ! How to use:

    ! First, create the instance
    real(real64),allocatable :: a(:,:)
    real(real64),allocatable :: b(:,:)
    real(real64),allocatable :: ab(:,:)
    real(real64),allocatable :: copyobj(:,:)
    real(real64),allocatable :: ret(:,:)
    integer(int32) :: i

    type(IO_) :: f

    ! #1 create array
    ! ----> allocate array
    allocate(ret(3,3) )
    ! ----> set identity matrix
    ret = identity_matrix(3)
    ! show array on terminal.
    print *, "ret = "
    call showArray(ret)

    ! #2 Array Input-Output
    ! ----> save matrix as "./test.txt" file
    call savetxt(ret,"./","test",".txt")
    ! ----> load a .txt or .csv file.
    a = loadtxt("./","test",".txt")
    b = loadtxt("./","test",".txt")

    ! -----> also, you can do it for multiple arrays as
    call f%open("./"//"arrays"//".txt")
    ! save "a" array
    call saveArray(f%fh,a)
    ! save "b" array
    call saveArray(f%fh,b)
    ! save "ret" array
    call saveArray(f%fh,ret)
    call f%close()

    ! -----> and,
    call f%open("./"//"arrays"//".txt")
    ! load "a" array
    call loadArray(f%fh,a)
    ! load "b" array
    call loadArray(f%fh,b)
    ! load "ret" array
    call loadArray(f%fh,ret)
    call f%close()

    ! -----> It shows the array.
    print *, "a = "
    call showArray(a)

    ! -----> copy array
    a(1,2)=20.0d0
    call copyarray(a,copyobj)
    ! -----> It shows the array.
    print *, "a (original) = "
    call showArray(a)
    print *, "copy = "
    call showArray(copyobj)

    ! Importance Index 7 / 10 : [*******   ]

end program

max_stars_repo_name: kazulagi/plantFEM


text:
	SUBROUTINE GQMFIL  ( mapnam, iret )
C************************************************************************
C* GQMFIL								*
C*									*
C* This subroutine returns the current map file name to be used by	*
C* GDRMAP to draw a map.						*
C*									*
C* GQMFIL  ( MAPNAM, IRET )						*
C*									*
C* Output parameters:							*
C*	MAPNAM 		CHAR*   	Map file name 			*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* I. Graffman/RDS	 3/85	GEMPLT Version 3.0			*
C* M. Goodman/RDS	 7/85	GEMPLT Version 3.1			*
C* M. desJardins/GSFC	 6/88	Cleaned up				*
C* M. desJardins/GSFC	 3/89	Eliminate the attempts to build the	*
C*				map file name since this was never	*
C*				really implemented			*
C* M. desJardins/GSFC	 4/90	Eliminate check for map mode		*
C************************************************************************
	INCLUDE 	'ERROR.PRM'
	INCLUDE 	'XYDEF.CMN'
C*
	CHARACTER*(*) 	mapnam
C------------------------------------------------------------------------
	iret   =  NORMAL
C
C*	Return the map file name from the common area.
C
	mapnam = mpfil
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
! Copyright (c) 2005-2010, 2012-2013, Andrew Hang Chen and contributors,
! All rights reserved.
! Licensed under the 3-clause BSD license.

!------------------------
!
! Test for FORTRAN unit test codes
!
! Author: Andrew H. Chen chena@westinghouse.com
!
!------------------------
!
module fruit_test
  use FRUIT
  implicit none

contains
  ! Run all the test sub routines
  ! -----------------------------
  subroutine fruit_test_package
    implicit none

    ! test printing and formats
    !--------------------------
    call test_should_add_success_or_failed_actions

    call test_show_output
    call fruit_summary_test
    call showOutputForReport

    ! the following should give 100% pass
    !------------------------------------
    call fruit_summary

    call test_assert_true

    call test_add_success
    call test_add_fail
    call test_is_all_successful
    call test_obsolete_message
    call fruit_summary

    call demonstrate_case_summary
    call fruit_summary

  end subroutine fruit_test_package

  subroutine test_should_add_success_or_failed_actions
    write (*,*) "Should see . here:"
    call add_success
    write (*,*) "Should see .. here:"
    call add_success
    call add_success

    write (*,*)

    write (*,*) "Should see F here:"
    call failed_assert_action('','')

    write (*,*) "Should see FF here:"
    call failed_assert_action('','')
    call failed_assert_action('','')

    write (*,*)
  end subroutine test_should_add_success_or_failed_actions

  subroutine test_assert_true
    implicit none

    call init_fruit
    write (*,*) 'Should see 1 successful case'
    call assert_true (.true.)

    write (*,*) 'Should see 1 failed case'
    call assert_true (.FALSE.)

    call fruit_summary
  end subroutine test_assert_true

  ! Test assert_true with message
  ! ----------------------------
  subroutine assert_true_message_test
    implicit none

    call init_fruit
    call assert_true (.true., 'message in assert_true_message_test true test')

    call assert_true (.FALSE., 'message in assert_true_message_test false test')

    call fruit_summary
  end subroutine assert_true_message_test

  ! Test add_success and is_all_successful
  ! -----------------------------------
  subroutine test_add_success
    implicit none

    logical :: result = .FALSE.

    call init_fruit

    call add_success
    call is_all_successful (result)

    IF (result .neqv. .true.) then
      write (*, *)
      write (*, *) 'FAILED add_success!!!'
      write (*, *)
    end IF

    call fruit_summary

  end subroutine test_add_success

  ! Test add_success and is_all_successful
  ! -----------------------------------
  subroutine test_add_fail
    implicit none

    logical :: result = .FALSE.

    call init_fruit
    call add_fail
    call is_all_successful (result)

    IF (result .neqv. .FALSE.) then
      write (*, *)
      write (*, *) 'FAILED add_fail !!!'
      write (*, *)
    end IF

    call add_fail('test_add_fail', 'testing add_fail')
    call fruit_summary

  end subroutine test_add_fail

  ! Test add_success and is_all_successful
  ! -----------------------------------
  subroutine test_add_fail_message
    implicit none

    logical :: result = .FALSE.

    call init_fruit
    call add_fail ('Add a failed case')
    call add_fail ('test_add_fail_message', 'Add a failed case')
    call is_all_successful (result)

    IF (result .neqv. .FALSE.) then
      write (*, *)
      write (*, *) 'FAILED add_fail !!!'
      write (*, *)
    end IF

  end subroutine test_add_fail_message

  ! Test successful and failed summary
  ! ----------------------------------
  subroutine fruit_summary_test
    implicit none

    write (*,*) 'Summary for successful cases without message: '
    write (*,*) 'Should see 1 successful case'
    call init_fruit
    call add_success
    call fruit_summary

    write (*,*) 'Summary for successful cases and 2 messages: '
    write (*,*) 'Should see 3 successful cases'
    call init_fruit
    call add_success
    call add_success
    call add_success
    call fruit_summary

    write (*,*) 'Summary for failed cases: '
    write (*,*) 'Should see 2 failed case amd 1 message'
    call init_fruit
    call add_fail
    call add_fail('Fail message from test case.')
    call fruit_summary

  end subroutine fruit_summary_test

  ! Add one failed case and assert false
  ! Add all success case, and assert true
  ! -------------------------------------
  subroutine test_is_all_successful
    implicit none

    logical :: result = .true.

    call init_fruit
    call add_success
    call add_fail
    call add_success
    call is_all_successful (result)

    IF (result .neqv. .FALSE.) then
      write (*, *)
      write (*, *) 'FAILED test_is_all_successful!!! (add_fail)'
      write (*, *)
    end IF

    call fruit_summary

    call init_fruit
    call add_success
    call add_success
    call add_success
    call is_all_successful (result)

    IF (result .neqv. .true.) then
      write (*, *)
      write (*, *) 'FAILED test_is_all_successful!!! (add_success)'
      write (*, *)
    end IF
    call fruit_summary

    call init_fruit
    IF (result .neqv. .true.) then
      write (*, *)
      write (*, *) 'FAILED test_is_all_successful!!! (add_success)'
      write (*, *)
    end IF
    call fruit_summary

  end subroutine test_is_all_successful

  subroutine test_show_output
    implicit none
    integer :: i

    call init_fruit
    DO i=1,5
       call assert_true (.true.)
    END DO
    call assert_equals (trim(get_last_message()), '')

    DO i=1,5
       call assert_true (.false.)
    END DO
    call assert_equals (trim(get_last_message()), 'Expected T got F')
    call fruit_summary

  end subroutine test_show_output

  subroutine showOutputForReport
    implicit none
    logical :: trueValue = .TRUE.
    logical :: falseValue = .FALSE.
    integer :: i
    integer :: count
    !character (100), DIMENSION (3)  :: msgs;

    call init_fruit
    DO i=1,2
      call assert_true (trueValue, 'msg')
    END DO

    DO i=1,2
      call assert_true (falseValue, 'msg')
    END DO

    call get_total_count (count)
    write (*, *) 'Total count is: ' , count;

    call get_failed_count (count)
    write (*, *) 'Failed count is: ' , count;

    ! to be implemented
    !call getMessages (msgs)
    !write (*, *) 'Messages are: ';
    !write (*, *) msgs;

    call fruit_summary
  end subroutine showOutputForReport

  subroutine test_obsolete_message
    print *, "Warning on obsolete subroutine name should appear below:"
    call assertTrue(.true.)
  end subroutine test_obsolete_message

  subroutine test_multiple_cases_1
     call assert_true (.true.)
     call assert_true (.false.)
  end subroutine test_multiple_cases_1

  subroutine test_multiple_cases_2
     call assert_true (.true.)
  end subroutine test_multiple_cases_2

  subroutine demonstrate_case_summary
     call run_test_case(test_multiple_cases_1)
     call run_test_case(test_multiple_cases_2)
  end subroutine demonstrate_case_summary

end module fruit_test

max_stars_repo_name: thomasms/fortran-examples

max_issues_repo_name: thomasms/fortran-examples

max_forks_repo_name: thomasms/fortran-examples


! Extracting substrings from a substrings

program subString
    implicit none

    CHARACTER(LEN=11) :: hello

    hello = "Hello World!"

   print*, hello(7:11)
end program subString


!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006, 2010, 2015, 2017
!        Salvatore Filippone    Cranfield University
!        Alfredo Buttari        CNRS-IRIT, Toulouse
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
!
! File: psi_crea_bnd_elem.f90
!
! Subroutine: psi_crea_bnd_elem
!    Extracts a list of boundary indices. If no boundary is present in
!    the distribution the output vector is put in the unallocated state,
!    otherwise its size is equal to the number of boundary indices on the
!    current (calling) process.
!
! Arguments:
!    bndel(:) - integer(psb_ipk_), allocatable      Array containing the output list
!    desc_a   - type(psb_desc_type).    The communication descriptor.
!    info     - integer.                  return code.
!
subroutine psi_crea_bnd_elem(bndel,desc_a,info)
  use psi_mod, psb_protect_name => psi_crea_bnd_elem
  use psb_realloc_mod
  use psb_desc_mod
  use psb_error_mod
  use psb_serial_mod
  implicit none

  integer(psb_ipk_), allocatable :: bndel(:)
  type(psb_desc_type), intent(in)  :: desc_a
  integer(psb_ipk_), intent(out) :: info

  integer(psb_ipk_), allocatable :: work(:)
  integer(psb_ipk_) :: i, j, nr, ns, k, err_act
  character(len=20)    :: name

  info = psb_success_
  name='psi_crea_bnd_elem'
  call psb_erractionsave(err_act)

  allocate(work(size(desc_a%halo_index)),stat=info)
  if (info /= psb_success_ ) then
    info = psb_err_alloc_dealloc_
    call psb_errpush(info,name)
    goto 9999
  end if

  i=0
  j=1
  do while(desc_a%halo_index(j) /= -1)

    nr = desc_a%halo_index(j+1)
    ns = desc_a%halo_index(j+1+nr+1)
    do k=1, ns
      i = i + 1
      work(i) = desc_a%halo_index(j+1+nr+1+k)
    enddo
    j  = j + 1 + ns + 1 + nr + 1
  enddo

  call psb_msort_unique(work(1:i),j)

  if (.true.) then
    if (j>=0) then
      call psb_realloc(j,bndel,info)
      if (info /= psb_success_) then
        call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
        goto 9999
      end if
      bndel(1:j) = work(1:j)
    else
      if (allocated(bndel)) then
        deallocate(bndel)
      end if
    end if
  else
    call psb_realloc(j+1,bndel,info)
    if (info /= psb_success_) then
      call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
      goto 9999
    end if
    bndel(1:j) = work(1:j)
    bndel(j+1) = -1
  endif

  deallocate(work)
  call psb_erractionrestore(err_act)
  return

9999 call psb_error_handler(err_act)

  return

end subroutine psi_crea_bnd_elem


text:
/*
vim:syntax=sml
*/

let min (a:Int) (b:Int) =
  if[c:Int.(a>=c) && (b>=c)] a >= b then
    b
  else
    a;;

let max (a:Int) (b:Int) =
  if[c:Int.(c>=a) && (c>=b)] a >= b then
    a
  else
    b;;

let abs (x:Int) =
  if[r:Int.r>=0] x <= 0 then
    x
  else
    (0 - x);;

let Pos = (x:Int. x > 0);;

let rec rem (n:Int) (d:Pos) : (r:Int.(r <= n) && (r < d)) =
  if[r:Int.(r <= n) && (r < d)] n < d then
    n
  else
    rem (n - d) d;;

let pred (x:Int) : (r:Int.r < x) = x - 1;;

let succ (x:Int) : (r:Int.r > x) = x + 1;;

abs (0 - 5);;

min 7 22;;

max 12 41;;

rem 6 4;;

pred 18;;

succ 9;;

max_stars_repo_name: ucsc-proglang/sage

max_issues_repo_name: ucsc-proglang/sage

max_forks_repo_name: ucsc-proglang/sage


text:
      module fmpi
        include "mpif.h"
      end module
!****************************************************************
!  This is a simple send/receive program in MPI
!  Processor 0 sends an integer to processor 1,
!  while processor 1 receives the integer from proc. 0
!****************************************************************
      program hello
      use fmpi
!     include "mpif.h"
      integer myid, ierr,numprocs
      integer tag,source,destination,count
      integer buffer
      integer status(MPI_STATUS_SIZE)
      call MPI_INIT( ierr )
      call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
      call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr )
      tag=1234
      source=0
      destination=1
      count=1
      do i =1,3
      if(myid .eq. source)then
         buffer=5678-i
         Call MPI_Send(buffer, count, MPI_INTEGER,destination,&
          tag, MPI_COMM_WORLD, ierr)
         write(*,*)"Fortran processor ",myid," sent ",buffer
      endif
      if(myid .eq. destination)then
         Call MPI_Recv(buffer, count, MPI_INTEGER,source,&
          tag, MPI_COMM_WORLD, status,ierr)
         write(*,*)"Fortran processor ",myid," got ",buffer
      endif
      enddo
      call MPI_FINALIZE(ierr)
      end

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
! Test for Ticket 561: the value of I is destroyed by the call to INC
! although the formal parameter is an array

      PROGRAM MAIN
      INTEGER I
      I = 2
      CALL INC(I)
      PRINT *, I
      END

      SUBROUTINE INC(A)
      INTEGER A(1)

      A(1) = A(1) + 1

      RETURN
      END


text:
program foo
  integer  :: a, b
  a = 10
  b = 20

  !$omp parallel
  print *, "From OpenMP ", a, b
  !$omp end parallel

  print *, "From outsitde ", a, b
end program foo

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


subroutine multigrid_defect(u_local, res1, res2, aP_u, aW_u, aE_u, aS_u, aN_u, f_u, dx_mg, dy_mg, nx, ny, alpha)
use kind_parameters
implicit none

integer:: i, j, nx, ny
real(kind=dp), intent(in), dimension(nx,ny):: u_local, aP_u, f_u, aW_u, aE_u, aS_u, aN_u
real(kind=dp), intent(out), dimension(nx,ny):: res1, res2
real(kind=dp), dimension(nx):: dx_mg
real(kind=dp), dimension(ny):: dy_mg
real(kind=dp):: alpha

! Calculate the residual r = f - Au, where u_local is the guessed value
! Input
! u_local      Matrix for variables
! f      Right Hand Side
! hx     grid spacing in x direction
! hy     grid spacing in y direction
!
! Output
! r      Residual

! Author: Pratanu Roy
! History:
! First Written: July 20, 2012

alpha = 1.0

do j = 2,ny-1

    do i = 2,nx-1

        res1(i,j) = alpha/(dx_mg(i)*dy_mg(j))*(aW_u(i,j)*u_local(i-1,j) + aE_u(i,j)*u_local(i+1,j) + aS_u(i,j)*u_local(i,j-1) + aN_u(i,j)*u_local(i,j+1) - aP_u(i,j)*u_local(i,j))
        res2(i,j) = f_u(i,j) + res1(i,j)

    end do

end do

end subroutine


    !-------------------------------------------------------------------------------------------------------------
    !
    !> \file    CompExcessGibbsEnergySUBL.f90
    !> \brief   Compute the partial molar excess Gibbs energy of mixing of solution phase constituents in a SUBL
    !!           or SUBLM solution phase.
    !> \author  M.H.A. Piro
    !> \date    January 17, 2013
    !> \sa      CompExcessGibbsEnergy.f90
    !> \sa      CompExcessGibbsEnergyRKMP.f90
    !> \sa      CompExcessGibbsEnergyQKTO.f90
    !
    !
    ! Revisions:
    ! ==========
    !
    !   Date            Programmer          Description of change
    !   ----            ----------          ---------------------
    !   01/17/2013      M.H.A. Piro         Original code.
    !   02/11/2013      M.H.A. Piro         Added capability of handling non-ideal terms.
    !   02/14/2013      M.H.A. Piro         Corrected handling of non-ideal terms with higher order
    !                                        mixing parameters (happy valentine's day).
    !   02/15/2013      M.H.A. Piro         Fix bug in handling higher order terms from yesterday.
    !   02/27/2013      M.H.A. Piro         The sum of site fractions may not necessarily equal unity. Apply
    !                                        an appropriate correction.
    !
    !
    ! Purpose:
    ! ========
    !
    !> \details The purpose of this subroutine is to compute the partial molar excess Gibbs energy of mixing
    !! (dPartialExcessGibbs) of all constituents in a non-ideal solution phase designated as 'SUBL'
    !! (Compound Energy Formalism) or 'SUBLM' (SUBL-Magnetic) phases.  For more information on the
    !! SUBL model and the derivation of equations, the reader is referred to the following literature:
    !!
    !!      B. Sundman and J. Agren, "A Regular Solution Model for Phases with Several Components
    !!      and Sublattices, Suitable for Computer Applications," Journal of Physical Chemistry of
    !!      Solids, 42 (1981) 297-301.
    !!
    !!      M. Hillert, "The Compound Energy Formalism," Journal of Alloys and Compounds, 320 (2001)
    !!      161-176.
    !!
    !! An important distinction between a phase represented by the SUBL model and any other model
    !! is that the phase components and the constituents are not the same entity (refer to above
    !! literature for a theoretical discussion).  As a result of this, the contribution to the chemical
    !! potential term from the reference molar Gibbs energy is NOT \f$ g_i^{\circ} \f$, but rather a
    !! more complicated expression involved all components in this phase.
    !!
    !! The chemical potential of a component in a SUBL phase is defined by the following equation:
    !!
    !! 	\f$ \mu_{i(\lambda)} = \sum_{j=1}^{N_{\lambda}} x_{j(\lambda)} g_{j(\lambda)}^{\circ}
    !!      \left(1 - N_s + \sum_{s=1}^{N_s} \frac{\delta_{i,j}}{y_{j(s)}} \right)
    !!      + RT \sum_{s=1}^{N_s} a_s \mathrm{ln} (y_{i(s)}) + g_{i(\lambda)}^{ex} \f$
    !!
    !! The partial molar excess Gibbs energy of mixing of a component of a SUBL phase is:
    !!
    !! 	\f$ g_{i(\lambda)}^{ex} = \sum_{p=1}^{N_p} \left( \prod_{m=1} y_{m(s)}  \right)
    !!      \sum_{z=0} {^zL_{j,k}} \left(  1 - (N_s + z) + \sum_{s=1}^{N_s}
    !!      \frac{\delta _{i,p}}{y_{i(s)}}   \right) \f$
    !!
    !! The mole fraction of any component in a SUBL phase is related to the site fractions of the
    !! constituents through the following multiplicative relationship:
    !!
    !! \f$ x_{i(\lambda)} = \prod_{s=1}^{N_s} y_{i(s)} \f$
    !!
    !! Similarly, the site fraction of any constituent is related to the mole fraction through the following
    !! summation:
    !!
    !! \f$ y_{c(s)} = \sum_{i=1}^{N_{\lambda}} x_{i(\lambda)} \delta_{i,c(s)} \f$
    !!
    !! Since the molar Gibbs energy terms are defined by the phase components and NOT the constituents, the
    !! numerics necessarily works with mole fractions (corresponding to components) rather than site fractions
    !! (corresponding to constituents).  Therefore, the mole fractions are modified at any iteration and the
    !! site fractions are computed from the mole fractions.  Note that the computation of the site fractions
    !! is a summation operation and it is thus insensitive to components with relatively small mole fractions.
    !! However, the calculation of a mole fraction from the site fractions is sensitive because it is a
    !! multiplicative function. A correction is applied to the mole fractions of all components after the site
    !! fractions are computed by defining the mole fraction as a function of the site fractions.
    !!
    !
    ! Pertinent variables:
    ! ====================
    !
    !> \param[in] iSolnIndex    Absolute index of a solution phase
    !
    ! nSpeciesPhase             An integec vector representing the last index of a species in a particular
    !                            solution phase
    ! nParamPhase               An integer vector representing the number of parameters for a particular phase.
    ! nSublattice               An integer scalar representing the number of sublattices associated with this
    !                            particular phase (used for convenience).
    ! nSublatticePhase          An integer vector representing the number of sublattices for each charged phase.
    ! iChargedPhaseID           An integer scalar representing the relative index of the charged phase.
    ! iParam                    An integer scalar representing the index number of a parameter.
    ! iPhaseSublattice          An integer vector representing the relative index of a charged phase corresponding
    !                            to each solution phase in the system.  Thus, it is equal to zero for a phase
    !                            that does not contain any sublattices.
    ! iRegularParam             An integer matrix representing the relative indices of species mixing for a
    !                            particular parmater.
    ! iFirst             An integer scalar representing the absolute index of the first species in a phase.
    ! iLast              An integer scalar representing the absolute index of the last species in a phase.
    ! dChemicalPotential        A double real vector representing the chemical potential for every species in the
    !                            system.
    ! dMolFraction              A double real vector representing hte mole fraction of all species in the system.
    ! dSiteFraction             A double real array representing the site fraction of each constituent on each
    !                            sublattice.  The first dimension corresponds to the charged phase index, the
    !                            second dimension corresponds to the sublattice index and the third dimension
    !                            corresponds to the constituent index.
    ! dExcessGibbsParam         A double real vector representing the molar excess Gibbs energy of mixing for
    !                            each subsystem.
    ! dPartialExcessGibbs       Partial molar excess Gibbs energy of mixing of species.
    ! dMolFraction              Current estimated mole fraction.
    ! cSolnPhaseType            A character vector representing the solution phase type.
    ! iConstituentSublattice    An integer array representing the constituent index for a particular sublattice
    !                            phase.  The first dimension refers to the relative phase index for charged phases,
    !                            the second dimension represents the sublattice index and the third dimension
    !                            corresponds to the relative component index (not absolute).
    ! iFirstParam               An integer scalar representing the constituent index of the first constituent
    !                            that is being mixed for a mixing parameter.
    ! iSeconParam               An integer scalar representing the constituent index of the second constituent
    !                            that is being mixed for a mixing parameter.
    ! iSubParam                 An integer scalar representing the sublattice index corresponding to the mixing
    !                            parameter.
    !
    !-------------------------------------------------------------------------------------------------------------

subroutine CompExcessGibbsEnergySUBL(iSolnIndex)

    USE ModuleThermo
    USE ModuleThermoIO
    USE ModuleGEMSolver

    implicit none

    integer :: i, j, k, l, m, n, s, c, d
    integer :: iSolnIndex, nSublattice, iChargedPhaseID, iMixType
    integer :: iFirst, iLast, iFirstParam, iSecondParam, iSubParam, iFirstParam2, iSecondParam2, iSubParam2
    integer :: iTempParam1, iTempParam2, iExponent, iTempSub, iThirdParam, iTernaryCon
    real(8) :: dTemp, dPreFactor, dFirstParam, dSecondParam, dFirstParam2, dSecondParam2, dThirdParam
    real(8) :: dTempParam1, dTempParam2, KD, dSum
    real(8), dimension(nMaxSublatticeSys):: dTempVec

    ! Only proceed if the correct phase type is selected:
    IF_SUBL: if ((cSolnPhaseType(iSolnIndex) == 'SUBL').OR.(cSolnPhaseType(iSolnIndex) == 'SUBLM')) then

        ! Define temporary variables for sake of convenience:
        iChargedPhaseID = iPhaseSublattice(iSolnIndex)
        nSublattice     = nSublatticePhase(iChargedPhaseID)
        iFirst          = nSpeciesPhase(iSolnIndex-1) + 1
        iLast           = nSpeciesPhase(iSolnIndex)

        ! Initialize variables:
        dSiteFraction(iChargedPhaseID,1:nSublattice,1:nMaxConstituentSys) = 0D0
        dChemicalPotential(iFirst:iLast)                                  = 0D0
        dPartialExcessGibbs(iFirst:iLast)                                 = 0D0
        dTempVec                                                          = 0D0

        ! Compute site fractions on each sublattice:
        LOOP_SiteFraction: do i = nSpeciesPhase(iSolnIndex-1) + 1, nSpeciesPhase(iSolnIndex)

            ! Relative component index:
            m = i - iFirst + 1

            ! Loop through sublattices:
            do s = 1, nSublattice
                ! Store constituent index on sublattice s:
                c = iConstituentSublattice(iChargedPhaseID,s,m)
                dSiteFraction(iChargedPhaseID,s,c) = dSiteFraction(iChargedPhaseID,s,c) + dMolFraction(i)
            end do
        end do LOOP_SiteFraction

        ! Compute sum of site fractions on each sublattice:
        do s = 1, nSublattice
            do c = 1, nMaxConstituentSys
                dTempVec(s) = dTempVec(s) + dSiteFraction(iChargedPhaseID,s,c)
            end do
            dTempVec(s) = 1D0 / dTempVec(s)
        end do

        ! Correct the mole fractions of phase components by the site fractions of the constituents
        ! (see top for details):
        LOOP_CorrectX: do i = iFirst, iLast
            dTemp = 1D0
            m = i - iFirst + 1

            ! Loop through sublattices:
            do s = 1, nSublattice
                c     = iConstituentSublattice(iChargedPhaseID,s,m)
                dTemp = dTemp * dSiteFraction(iChargedPhaseID,s,c) * dTempVec(s)
            end do

            dMolFraction(i)  = dTemp
        end do LOOP_CorrectX

        ! REFERENCE GIBBS ENERGY AND IDEAL MIXING
        ! ---------------------------------------

        ! Compute the chemical potential for each phase component assuming ideal mixing:
        LOOP_Ideal: do i = iFirst, iLast
            ! Relative species index:
            m = i - iFirst + 1

            ! Loop through components:
            LOOP_Ideal_Components: do j = iFirst, iLast
                ! Relative species index:
                n = j - iFirst + 1

                ! Compute pre-factor term:
                dTemp = 1D0 - DFLOAT(nSublattice)

                ! Loop through sublattices:
                do s = 1, nSublattice
                    ! Store constituent indices:
                    k = iConstituentSublattice(iChargedPhaseID,s,m)
                    l = iConstituentSublattice(iChargedPhaseID,s,n)

                    ! Effectively apply Kronecker-Delta term to pre-factor:
                    if (k == l)  dTemp = dTemp + 1D0 / dSiteFraction(iChargedPhaseID,s,k)
                end do

                ! Update the reference molar Gibbs energy:
                dChemicalPotential(i) = dChemicalPotential(i) + dTemp * dMolFraction(j) * dStdGibbsEnergy(j)
            end do LOOP_Ideal_Components

            ! Add ideal mixing contribution:
            do s = 1, nSublattice
                c = iConstituentSublattice(iChargedPhaseID,s,m)
                dChemicalPotential(i) = dChemicalPotential(i) + dStoichSublattice(iChargedPhaseID,s) &
                    * DLOG(dSiteFraction(iChargedPhaseID,s,c))
            end do

            ! Sum stoichiometry and add large penalty if this species is vacancies only and a large mole fraction
            dSum = 0D0
            do j = 1, nElements
                dSum = dSum + ABS(dStoichSpecies(i,j))
            end do
            if ((dSum == 0D0) .AND. (dMolFraction(i) > 0.9D0)) dChemicalPotential(i) = dChemicalPotential(i) + 1D3
        end do LOOP_Ideal

        ! MAGNETIC TERMS
        ! --------------

        ! Compute magnetic ordering terms if this phase is magnetic:
        if (cSolnPhaseType(iSolnIndex) == 'SUBLM') call CompGibbsMagneticSoln(iSolnIndex)

        ! EXCESS TERMS
        ! ------------

        ! Return control to the parent subroutine if there aren't any interaction parameters for this phase:
        if (nParamPhase(iSolnIndex) - nParamPhase(iSolnIndex-1) == 0) return

        ! Loop through parameters:
        LOOP_Param: do l = nParamPhase(iSolnIndex-1) + 1, nParamPhase(iSolnIndex)

            ! Reinitialize temporary variable:
            dPreFactor = 1D0
            iFirstParam = 0
            iSecondParam = 0
            iThirdParam = 0
            dFirstParam = 0D0
            dSecondParam = 0D0
            dThirdParam = 0D0
            iFirstParam2 = 0
            iSecondParam2 = 0

            iMixType = 0

            ! Store the number of constituents involved in this parameter:
            n = iRegularParam(l,1)
            ! Determine the mixing parameter type
            if ((iSUBLParamData(l,1) == 1) .AND. (iSUBLParamData(l,3) == 2)) then
                iMixType = 2

                ! Loop through constituents associated with this parameter:
                do k = 2, n + 1
                    ! Determine constituent and sublattice indices:
                    c = MOD(iRegularParam(l,k), 10000)
                    s = iRegularParam(l,k) - c
                    s = s / 10000

                    ! Compute prefactor term:
                    dPreFactor = dPreFactor * dSiteFraction(iChargedPhaseID,s,c)

                    ! Store the first and second site fractions:
                    if (k == iSUBLParamData(l,2)) then
                        dFirstParam = dSiteFraction(iChargedPhaseID,s,c)
                        iFirstParam = c
                        iSubParam   = s
                    else if (k == iSUBLParamData(l,2) + 1) then
                        dSecondParam = dSiteFraction(iChargedPhaseID,s,c)
                        iSecondParam = c
                    end if
                end do

                ! Multiply prefactor term by excess Gibbs energy parameter:
                iExponent = iRegularParam(l,n+2)
                dPreFactor = dPreFactor * dExcessGibbsParam(l) * (dFirstParam - dSecondParam)**(iExponent)
            else if ((iSUBLParamData(l,1) == 1) .AND. (iSUBLParamData(l,3) == 3)) then
                iMixType = 3
                iTernaryCon = iRegularParam(l,n+2)

                ! Loop through constituents associated with this parameter:
                do k = 2, n + 1
                    ! Determine constituent and sublattice indices:
                    c = MOD(iRegularParam(l,k), 10000)
                    s = iRegularParam(l,k) - c
                    s = s / 10000

                    ! Compute prefactor term:
                    dPreFactor = dPreFactor * dSiteFraction(iChargedPhaseID,s,c)

                    ! Store the first and second site fractions:
                    if (k == iSUBLParamData(l,2) + iTernaryCon) then
                        dFirstParam = dSiteFraction(iChargedPhaseID,s,c)
                        iFirstParam = c
                        iSubParam   = s
                    else if (k == iSUBLParamData(l,2) + MOD(iTernaryCon + 1, 3)) then
                        dSecondParam = dSiteFraction(iChargedPhaseID,s,c)
                        iSecondParam = c
                    else if (k == iSUBLParamData(l,2) + MOD(iTernaryCon + 2, 3)) then
                        dThirdParam = dSiteFraction(iChargedPhaseID,s,c)
                        iThirdParam = c
                    end if
                end do

                ! Multiply prefactor term by excess Gibbs energy parameter:
                iExponent = 0
                dPreFactor = dPreFactor * dExcessGibbsParam(l)
                dPreFactor = dPreFactor * (dFirstParam + (1D0 - dFirstParam - dSecondParam - dThirdParam) / 3D0)
            else if ((iSUBLParamData(l,1) == 2) .AND. (iSUBLParamData(l,3) == 2) .AND. (iSUBLParamData(l,5) == 2)) then
                iMixType = 4
                ! Loop through constituents associated with this parameter:
                do k = 2, n + 1

                    ! Determine constituent and sublattice indices:
                    c = MOD(iRegularParam(l,k), 10000)
                    s = iRegularParam(l,k) - c
                    s = s / 10000

                    ! Compute prefactor term:
                    dPreFactor = dPreFactor * dSiteFraction(iChargedPhaseID,s,c)

                    ! Store the site fractions:
                    if (k == iSUBLParamData(l,2)) then
                        dFirstParam = dSiteFraction(iChargedPhaseID,s,c)
                        iFirstParam = c
                        iSubParam   = s
                    else if (k == iSUBLParamData(l,2) + 1) then
                        dSecondParam = dSiteFraction(iChargedPhaseID,s,c)
                        iSecondParam = c
                    else if (k == iSUBLParamData(l,4)) then
                        dFirstParam2 = dSiteFraction(iChargedPhaseID,s,c)
                        iFirstParam2 = c
                        iSubParam2   = s
                    else if (k == iSUBLParamData(l,4) + 1) then
                        dSecondParam2 = dSiteFraction(iChargedPhaseID,s,c)
                        iSecondParam2 = c
                    end if
                end do

                ! Multiply prefactor term by excess Gibbs energy parameter:
                if (MOD(iRegularParam(l,n+2),2) == 0) then
                    iExponent = iRegularParam(l,n+2)/2
                else if (MOD(iRegularParam(l,n+2),2) == 1) then
                    iExponent = (iRegularParam(l,n+2)-1)/2
                    iTempParam1 = iFirstParam
                    iTempParam2 = iSecondParam
                    dTempParam1 = dFirstParam
                    dTempParam2 = dSecondParam

                    iFirstParam = iFirstParam2
                    iSecondParam = iSecondParam2
                    dFirstParam = dFirstParam2
                    dSecondParam = dSecondParam2

                    iFirstParam2 = iTempParam1
                    iSecondParam2 = iTempParam2
                    dFirstParam2 = dTempParam1
                    dSecondParam2 = dTempParam2

                    iTempSub = iSubParam
                    iSubParam = iSubParam2
                    iSubParam2 = iTempSub
                end if

                dPreFactor = dPreFactor * dExcessGibbsParam(l) * (dFirstParam - dSecondParam)**(iExponent)

                ! print *, dExcessGibbsParam(l), dFirstParam, dSecondParam, dFirstParam2, dSecondParam2
            else
                print *, 'Unrecognized excess mixing term in SUBL phase ', cSolnPhaseName(iSolnIndex)
                INFOThermo = 36
                return
            end if

            ! Loop through species in phase:
            LOOP_Param_Species: do i = iFirst, iLast

                ! Reinitialize variables:
                KD    = 0D0
                m     = i - iFirst + 1
                if ((iMixType == 2) .OR. (iMixType == 4)) then
                    dTemp = -DFLOAT(nSublattice + iExponent)
                else if (iMixType == 3) then
                    dTemp = -DFLOAT(nSublattice + 1)
                end if

                ! Loop through sublattices associated with this phase:
                LOOP_Param_Sub: do s = 1, nSublattice
                    ! Store constituent index corresponding to component i on sublattice s:
                    c = iConstituentSublattice(iChargedPhaseID,s,m)

                    ! Assign Kronecker-Delta term if this component contains the constituent corresponding
                    ! to the mixing parameter:
                    if (s == iSubParam) then
                        if (c == iFirstParam) then
                            ! This is the first mixing constituent:
                            if ((iMixType == 2) .OR. (iMixType == 4)) then
                                KD = 1D0
                            else if (iMixType == 3) then
                                KD = 2D0/3D0
                            end if
                        elseif (c == iSecondParam) then
                            ! This is the second mixing constituent:
                            if ((iMixType == 2) .OR. (iMixType == 4)) then
                                KD = -1D0
                            else if (iMixType == 3) then
                                KD = -1D0/3D0
                            end if
                        elseif (c == iThirdParam) then
                            ! This is the third mixing constituent:
                            if (iMixType == 3) then
                                KD = -1D0/3D0
                            end if
                        else
                            ! The constituents don't match:
                            KD = 0D0
                        end if
                        if (iMixType == 3) KD = KD - (dFirstParam + (- dFirstParam - dSecondParam - dThirdParam) / 3D0)
                    end if

                    ! Loop through constituents involved in this mixing parameter:
                    LOOP_Param_Const: do j = 2, n + 1
                        k = MOD(iRegularParam(l,j), 10000)  ! Constituent index corresponding to parameter.
                        d = iRegularParam(l,j) - k
                        d = d / 10000                       ! Sublattice index corresponding to parameter.

                        ! Cycle if they are on different sublattices:
                        if (d /= s) cycle LOOP_Param_Const

                        ! Cycle if they are different constituents:
                        if (k /= c) cycle LOOP_Param_Const

                        ! Include contribution from this site fraction:
                        dTemp = dTemp + 1D0 / dSiteFraction(iChargedPhaseID,s,c)

                    end do LOOP_Param_Const ! j
                end do LOOP_Param_Sub       ! s

                ! Apply higher order terms
                if (iMixType == 2) then
                    ! only if dFirstParam and dSecondParam are not the same:
                    if (dFirstParam /= dSecondParam) then
                        dTemp = dTemp + KD * DFLOAT(iExponent) / (dFirstParam - dSecondParam)
                    end if
                else if (iMixType == 3) then
                    dTemp = dTemp + KD / (dFirstParam + (1D0 - dFirstParam - dSecondParam - dThirdParam) / 3D0)
                else if (iMixType == 4) then
                    ! only if dFirstParam and dSecondParam are not the same:
                    if (dFirstParam /= dSecondParam) then
                        dTemp = dTemp + KD * DFLOAT(iExponent) / (dFirstParam - dSecondParam)
                    end if
                    dTemp = dTemp - 1D0
                end if

                ! Apply partial molar excess Gibbs energy of mixing:
                dPartialExcessGibbs(i) = dPartialExcessGibbs(i) + dPreFactor * dTemp

            end do LOOP_Param_Species       ! i

        end do LOOP_Param                   ! l

    end if IF_SUBL

    return

end subroutine CompExcessGibbsEnergySUBL


text:
enum E { A, B, C }
constant c = E

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text:
      program TRADM
      integer B, M
      if (.TRUE.) then
         call TRAD(B)
      endif
      end

      subroutine TRAD(B)
      integer B
      read *,B
      end


text:
      integer pid_pdg(mxpart)
      common/pid_pdg/pid_pdg

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
!***********************************************************************
! Computes streamflow and other variables in subbasins
!
!     version: 3.1 (rsregan, April 2008 ) changed code to assume
!                  nhru=nssr=ngw
!     version: 3.0 (rsregan, April 2008 ) removed all but subbasin code
!     version: 2.1 (rsregan, June 2007)
!     version: 2 (lhay, May 2007)
!     version: 1.1 (rsregan, March 2007)
!     version: 1 (lhay, March 2007)
!***********************************************************************

      MODULE PRMS_SUBBASIN
      IMPLICIT NONE
!   Local Variables
      CHARACTER(LEN=8), SAVE :: MODNAME
      DOUBLE PRECISION, SAVE, ALLOCATABLE :: Qsub(:), Sub_area(:), Laststor(:)
      INTEGER, SAVE, ALLOCATABLE :: Tree(:, :)
!   Declared Variables
      DOUBLE PRECISION, SAVE, ALLOCATABLE :: Sub_inq(:), Sub_cfs(:), Sub_cms(:)
      DOUBLE PRECISION, SAVE, ALLOCATABLE :: Sub_interflow(:), Sub_gwflow(:), Sub_sroff(:)
      DOUBLE PRECISION, SAVE, ALLOCATABLE :: Subinc_interflow(:), Subinc_gwflow(:), Subinc_sroff(:)
      DOUBLE PRECISION, SAVE, ALLOCATABLE :: Subinc_precip(:), Subinc_snowmelt(:), Subinc_pkweqv(:)
      DOUBLE PRECISION, SAVE, ALLOCATABLE :: Subinc_actet(:), Subinc_potet(:)
      DOUBLE PRECISION, SAVE, ALLOCATABLE :: Subinc_swrad(:), Subinc_snowcov(:)
      DOUBLE PRECISION, SAVE, ALLOCATABLE :: Subinc_tmaxc(:), Subinc_tminc(:), Subinc_tavgc(:)
      DOUBLE PRECISION, SAVE, ALLOCATABLE :: Subinc_deltastor(:), Subinc_wb(:)
      DOUBLE PRECISION, SAVE, ALLOCATABLE :: Subinc_rain(:), Subinc_snow(:), Subinc_stor(:)
      DOUBLE PRECISION, SAVE, ALLOCATABLE :: Subinc_recharge(:), Subinc_szstor_frac(:), Subinc_capstor_frac(:)
!   Declared Parameters
      INTEGER, SAVE, ALLOCATABLE :: Subbasin_down(:), Hru_subbasin(:)
      END MODULE PRMS_SUBBASIN

!***********************************************************************
!     Main daily stream flow routine
!***********************************************************************
      INTEGER FUNCTION subbasin()
      USE PRMS_MODULE, ONLY: Process
      IMPLICIT NONE
! Functions
      INTEGER, EXTERNAL :: subdecl, subinit, subrun
!***********************************************************************
      subbasin = 0

      IF ( Process(:3)=='run' ) THEN
        subbasin = subrun()
      ELSEIF ( Process(:4)=='decl' ) THEN
        subbasin = subdecl()
      ELSEIF ( Process(:4)=='init' ) THEN
        subbasin = subinit()
      ENDIF

      END FUNCTION subbasin

!***********************************************************************
!     subdecl - set up parameters for streamflow, lake computations,
!               and subbasin flow
!   Declared Parameters
!     hru_area, subbasin_down, hru_subbasin
!***********************************************************************
      INTEGER FUNCTION subdecl()
      USE PRMS_SUBBASIN
      USE PRMS_MODULE, ONLY: Model, Nsub, Nhru
      IMPLICIT NONE
! Functions
      INTEGER, EXTERNAL :: declparam, declvar
      EXTERNAL read_error, print_module
! Local Variables
      CHARACTER(LEN=80), SAVE :: Version_subbasin
!***********************************************************************
      subdecl = 0

      Version_subbasin = 'subbasin.f90 2017-11-15 15:11:00Z'
      CALL print_module(Version_subbasin, 'Output Summary              ', 90)
      MODNAME = 'subbasin'

      IF ( Nsub==0 ) THEN
        IF ( Model/=99 ) STOP 'ERROR, nsub=0 when subbasin module called'
        Nsub = 1
      ENDIF

! Declared Variables
      ALLOCATE ( Sub_interflow(Nsub) )
      IF ( declvar(MODNAME, 'sub_interflow', 'nsub', Nsub, 'double', &
     &     'Area-weighted average interflow to each subbasin from'// &
     &     ' associated HRUs and from upstream subbasins', &
     &     'cfs', Sub_interflow)/=0 ) CALL read_error(3, 'sub_interflow')

      IF ( Model/=0 ) THEN
        ALLOCATE ( Sub_gwflow(Nsub) )
        IF ( declvar(MODNAME, 'sub_gwflow', 'nsub', Nsub, 'double', &
     &       'Area-weighted average groundwater discharge from'// &
     &       ' associated GWRs to each subbasin and from upstream subbasins', &
     &       'cfs', Sub_gwflow)/=0 ) CALL read_error(3, 'sub_gwflow')
        ALLOCATE ( Subinc_gwflow(Nsub) )
        IF ( declvar(MODNAME, 'subinc_gwflow', 'nsub', Nsub, 'double', &
     &       'Area-weighted average groundwater discharge from associated  GWRs to each subbasin', &
     &       'cfs', Subinc_gwflow)/=0 ) CALL read_error(3, 'subinc_gwflow')
      ENDIF

      ALLOCATE ( Sub_sroff(Nsub) )
      IF ( declvar(MODNAME, 'sub_sroff', 'nsub', Nsub, 'double', &
     &     'Area-weighted average surface runoff from associated HRUs'// &
     &     ' to each subbasin and from upstream subbasins', &
     &     'cfs', Sub_sroff)/=0 ) CALL read_error(3, 'sub_sroff')

      ALLOCATE ( Subinc_snowcov(Nsub) )
      IF ( declvar(MODNAME, 'subinc_snowcov', 'nsub', Nsub, 'double', &
     &     'Area-weighted average snow-covered area for associated HRUs for each subbasin', &
     &     'decimal fraction', Subinc_snowcov)/=0 ) CALL read_error(3, 'subinc_snowcov')

      ALLOCATE ( Subinc_interflow(Nsub) )
      IF ( declvar(MODNAME, 'subinc_interflow', 'nsub', Nsub, 'double', &
     &     'Area-weighted average interflow from associated HRUs to each subbasin', &
     &     'cfs', Subinc_interflow)/=0 ) CALL read_error(3, 'subinc_interflow')

      ALLOCATE ( Subinc_sroff(Nsub) )
      IF ( declvar(MODNAME, 'subinc_sroff', 'nsub', Nsub, 'double', &
     &    'Area-weighted average surface runoff from associated HRUs to each subbasin', &
     &    'cfs', Subinc_sroff)/=0 ) CALL read_error(3, 'subinc_sroff')

      ALLOCATE ( Subinc_precip(Nsub) )
      IF ( declvar(MODNAME, 'subinc_precip', 'nsub', Nsub, 'double', &
     &     'Area-weighted average precipitation from associated HRUs to each subbasin', &
     &     'inches', Subinc_precip)/=0 ) CALL read_error(3, 'subinc_precip')

      ALLOCATE ( Subinc_rain(Nsub) )
      IF ( declvar(MODNAME, 'subinc_rain', 'nsub', Nsub, 'double', &
     &     'Area-weighted average rain on associated HRUs to each subbasin', &
     &     'inches', Subinc_rain)/=0 ) CALL read_error(3, 'subinc_rain')

      ALLOCATE ( Subinc_snow(Nsub) )
      IF ( declvar(MODNAME, 'subinc_snow', 'nsub', Nsub, 'double', &
     &     'Area-weighted average snow on associated HRUs to each subbasin', &
     &     'inches', Subinc_snow)/=0 ) CALL read_error(3, 'subinc_snow')

      ALLOCATE ( Subinc_actet(Nsub) )
      IF ( declvar(MODNAME, 'subinc_actet', 'nsub', Nsub, 'double', &
     &     'Area-weighted average actual ET from associated HRUs to each subbasin', &
     &     'inches', Subinc_actet)/=0 ) CALL read_error(3, 'subinc_actet')

      ALLOCATE ( Subinc_potet(Nsub) )
      IF ( declvar(MODNAME, 'subinc_potet', 'nsub', Nsub, 'double', &
     &     'Area-weighted average potential ET from associated HRUs to each subbasin', &
     &     'inches', Subinc_potet)/=0 ) CALL read_error(3, 'subinc_potet')

      ALLOCATE ( Subinc_swrad(Nsub) )
      IF ( declvar(MODNAME, 'subinc_swrad', 'nsub', Nsub, 'double', &
     &     'Area-weighted average shortwave radiation distributed'// &
     &     ' to associated HRUs of each subbasin', &
     &     'Langleys', Subinc_swrad)/=0 ) CALL read_error(3, 'subinc_swrad')

      ALLOCATE ( Subinc_tminc(Nsub) )
      IF ( declvar(MODNAME, 'subinc_tminc', 'nsub', Nsub, 'double', &
     &     'Area-weighted average minimum air temperature for'// &
     &     ' associated HRUs to each subbasin', &
     &     'degrees Celsius', Subinc_tminc)/=0 ) CALL read_error(3, 'subinc_tminc')

      ALLOCATE ( Subinc_tmaxc(Nsub) )
      IF ( declvar(MODNAME, 'subinc_tmaxc', 'nsub', Nsub, 'double', &
     &     'Area-weighted average maximum air temperature for'// &
     &     ' associated HRUs to each subbasin', &
     &     'degrees Celsius', Subinc_tmaxc)/=0 ) CALL read_error(3, 'subinc_tmaxc')

      ALLOCATE ( Subinc_tavgc(Nsub) )
      IF ( declvar(MODNAME, 'subinc_tavgc', 'nsub', Nsub, 'double', &
     &     'Area-weighted average air temperature for associated HRUs to each subbasin', &
     &     'degrees Celsius', Subinc_tavgc)/=0 ) CALL read_error(3, 'subinc_tavgc')

      ALLOCATE ( Subinc_wb(Nsub) )
      IF ( declvar(MODNAME, 'subinc_wb', 'nsub', Nsub, 'double', &
     &     'Water balance for each subbasin', &
     &     'inches', Subinc_wb)/=0 ) CALL read_error(3, 'subinc_wb')

      ALLOCATE ( Subinc_deltastor(Nsub) )
      IF ( declvar(MODNAME, 'subinc_deltastor', 'nsub', Nsub, 'double', &
     &     'Change in storage for each subbasin', &
     &     'inches', Subinc_deltastor)/=0 ) CALL read_error(3, 'subinc_deltastor')

      ALLOCATE ( Subinc_snowmelt(Nsub) )
      IF ( declvar(MODNAME, 'subinc_snowmelt', 'nsub', Nsub, 'double', &
     &     'Area-weighted average snowmelt from associated HRUs of each subbasin', &
     &     'inches', Subinc_snowmelt)/=0 ) CALL read_error(3, 'subinc_snowmelt')

      ALLOCATE ( Subinc_pkweqv(Nsub) )
      IF ( declvar(MODNAME, 'subinc_pkweqv', 'nsub', Nsub, 'double', &
     &     'Area-weighted average snowpack water equivalent from'// &
     &     ' associated HRUs of each subbasin', &
     &     'inches', Subinc_pkweqv)/=0 ) CALL read_error(3, 'subinc_pkweqv')

      ALLOCATE ( Subinc_recharge(Nsub) )
      IF ( declvar(MODNAME, 'subinc_recharge', 'nsub', Nsub, 'double', &
     &     'Area-weighted average recharge from associated HRUs of each subbasin', &
     &     'inches', Subinc_recharge)/=0 ) CALL read_error(3, 'subinc_recharge')

      ALLOCATE ( Subinc_szstor_frac(Nsub) )
      IF ( declvar(MODNAME, 'subinc_szstor_frac', 'nsub', Nsub, 'double', &
     &     'Area-weighted average fraction of soil-zone water content storage for associated HRUs of each subbasin', &
     &     'decimal fraction', Subinc_szstor_frac)/=0 ) CALL read_error(3, 'subinc_szstor_frac')

      ALLOCATE ( Subinc_capstor_frac(Nsub) )
      IF ( declvar(MODNAME, 'subinc_capstor_frac', 'nsub', Nsub, 'double', &
     &     'Area-weighted average fraction of capillary reservoir water content storage for associated HRUs of each subbasin', &
     &     'decimal fraction', Subinc_capstor_frac)/=0 ) CALL read_error(3, 'subinc_capstor_frac')

      ALLOCATE ( Subinc_stor(Nsub) )
      IF ( declvar(MODNAME, 'subinc_stor', 'nsub', Nsub, 'double', &
     &     'Area-weighted average total water content in storage reservoirs associated HRUs of each subbasin', &
     &     'inches', subinc_stor)/=0 ) CALL read_error(3, 'subinc_stor')

      ALLOCATE ( Qsub(Nsub), Tree(Nsub, Nsub), Sub_inq(Nsub) )
      IF ( declvar(MODNAME, 'sub_inq', 'nsub', Nsub, 'double', &
     &     'Sum of streamflow from upstream subbasins to each subbasin', &
     &     'cfs', Sub_inq)/=0 ) CALL read_error(3, 'sub_inq')

      ALLOCATE ( Sub_cfs(Nsub) )
      IF ( declvar(MODNAME, 'sub_cfs', 'nsub', Nsub, 'double', &
     &     'Total streamflow leaving each subbasin', &
     &     'cfs', Sub_cfs)/=0 ) CALL read_error(3, 'sub_cfs')

      ALLOCATE ( Sub_cms(Nsub) )
      IF ( declvar(MODNAME, 'sub_cms', 'nsub', Nsub, 'double', &
     &     'Total streamflow leaving each subbasin', &
     &     'cms', Sub_cms)/=0 ) CALL read_error(3, 'sub_cms')

      ALLOCATE ( Subbasin_down(Nsub) )
      IF ( declparam(MODNAME, 'subbasin_down', 'nsub', 'integer', &
     &     '0', 'bounded', 'nsub', &
     &     'Downstream subbasin for each subbasin', &
     &     'Index number for the downstream subbasin whose inflow is outflow from this subbasin', &
     &     'none')/=0 ) CALL read_error(1, 'subbasin_down')

      ALLOCATE ( Hru_subbasin(Nhru) )
      IF ( declparam(MODNAME, 'hru_subbasin', 'nhru', 'integer', &
     &     '0', 'bounded', 'nsub', &
     &     'Index of subbasin assigned to each HRU', &
     &     'Index of subbasin assigned to each HRU', &
     &     'none')/=0 ) CALL read_error(1, 'hru_subbasin')

! Allocate arrays for variables
      ALLOCATE ( Sub_area(Nsub), Laststor(Nsub) )

      END FUNCTION subdecl

!***********************************************************************
!     subinit - Initialize subbasin module - get parameter values,
!               compute initial values
!***********************************************************************
      INTEGER FUNCTION subinit()
      USE PRMS_SUBBASIN
      USE PRMS_MODULE, ONLY: Model, Nsub, Nhru, Print_debug, Inputerror_flag, Dprst_flag, Cascade_flag !, Lake_route_flag
      USE PRMS_BASIN, ONLY: Hru_area_dble, Active_hrus, Hru_route_order, &
     &    Hru_type, Hru_frac_perv, DNEARZERO, Cfs2cms_conv !, Lake_hru_id
      USE PRMS_FLOWVARS, ONLY: Ssres_stor, Soil_moist, Pkwater_equiv, Gwres_stor, Sroff, Ssres_flow
      USE PRMS_SET_TIME, ONLY: Cfs_conv !, Cfs2inches
      USE PRMS_INTCP, ONLY: Hru_intcpstor
      USE PRMS_SRUNOFF, ONLY: Hru_impervstor, Hortonian_lakes, Dprst_stor_hru
      USE PRMS_SOILZONE, ONLY: Lakein_sz
      USE PRMS_GWFLOW, ONLY: Gwres_flow
!      USE PRMS_LAKE_ROUTE, ONLY: Lake_vol, Lake_outcfs
      IMPLICIT NONE
! Functions
      INTRINSIC :: DBLE
      INTEGER, EXTERNAL :: getparam
      EXTERNAL read_error, PRMS_open_module_file
! Local Variables
      INTEGER :: i, j, k, kk, TREEUNIT
      DOUBLE PRECISION :: harea, gwstor, soilstor, snowstor, landstor, srq, ssq, gwq
!***********************************************************************
      subinit = 0

      IF ( getparam(MODNAME, 'hru_subbasin', Nhru, 'integer', Hru_subbasin)/=0 ) CALL read_error(2, 'hru_subbasin')
      IF ( getparam(MODNAME, 'subbasin_down', Nsub, 'integer', Subbasin_down)/=0 ) CALL read_error(2, 'subbasin_down')

! Determine the tree structure for the internal nodes
      Tree = 0

      DO j = 1, Nsub
        kk = Subbasin_down(j)
        IF ( kk/=0 ) Tree(kk, j) = 1
      ENDDO

      Sub_cfs = 0.0D0
      Sub_cms = 0.0D0
      Sub_inq = 0.0D0
      Subinc_interflow = 0.0D0
      IF ( Model/=0 ) THEN
        Subinc_gwflow = 0.0D0
        Sub_gwflow = 0.0D0
      ENDIF
      Subinc_sroff = 0.0D0
      Subinc_precip = 0.0D0
      Subinc_rain = 0.0D0
      Subinc_snow = 0.0D0
      Subinc_snowmelt = 0.0D0
      Subinc_pkweqv = 0.0D0
      Subinc_actet = 0.0D0
      Subinc_snowcov = 0.0D0
      Subinc_swrad = 0.0D0
      Subinc_tminc = 0.0D0
      Subinc_tmaxc = 0.0D0
      Subinc_tavgc = 0.0D0
      Subinc_potet = 0.0D0
      Subinc_wb = 0.0D0
      Subinc_deltastor = 0.0D0
      Subinc_recharge = 0.0D0
      Subinc_szstor_frac = 0.0D0
      Subinc_capstor_frac = 0.0D0
      Sub_interflow = 0.0D0
      Sub_sroff = 0.0D0

      IF ( Print_debug==14 ) THEN
        CALL PRMS_open_module_file(TREEUNIT, 'tree_structure')
        WRITE ( TREEUNIT, 9001 ) (j, j=1, Nsub)
        DO j = 1, Nsub
          WRITE ( TREEUNIT, 9002 ) j, (Tree(j,k), k=1, Nsub)
        ENDDO
      ENDIF

      DO j = 1, Nsub
        DO k = 1, Nsub
          IF ( Tree(j,k)>0 ) THEN
            DO i = 1, Nsub
              Tree(j, i) = Tree(j, i) + Tree(k, i)
              IF ( Tree(j,i)>1 ) THEN
                Tree(j, i) = 1
              ENDIF
            ENDDO
          ENDIF
        ENDDO
      ENDDO

      DO j = Nsub, 1, -1
        DO k = 1, Nsub
          IF ( Tree(j,k)>0 ) THEN
            DO i = 1, Nsub
              Tree(j, i) = Tree(j, i) + Tree(k, i)
              IF ( Tree(j,i)>1 ) Tree(j,i)=1
            ENDDO
          ENDIF
        ENDDO
      ENDDO

      IF ( Print_debug==14 ) THEN
        WRITE ( TREEUNIT, 9003 ) (j, j=1, Nsub)
        DO j = 1, Nsub
          WRITE ( TREEUNIT, 9002 ) j, (Tree(j,k), k=1, Nsub)
        ENDDO
        CLOSE ( TREEUNIT )
      ENDIF

! added some code to allow for restart, but not climate states and fluxes and subinc_deltastor
      Subinc_stor = 0.0D0
      Sub_area = 0.0D0
      gwstor = 0.0D0
      gwq = 0.0D0
      Qsub = 0.0D0
      DO i = 1, Active_hrus
        j = Hru_route_order(i)
        ! k indicates which HRU is in which subbasin
        k = Hru_subbasin(j)
        IF ( k>0 ) THEN
          harea = Hru_area_dble(j)
          IF ( Model/=0 ) THEN
            gwstor = Gwres_stor(j)*harea
            gwq = DBLE(Gwres_flow(j))*harea
          ENDIF
          IF ( Hru_type(j)/=2 ) THEN
            srq = DBLE(Sroff(j))*harea
            ssq = DBLE(Ssres_flow(j))*harea
            soilstor = DBLE(Soil_moist(j)*Hru_frac_perv(j) + Ssres_stor(j))*harea
            snowstor = Pkwater_equiv(j)*harea
            landstor = DBLE(Hru_intcpstor(j)+Hru_impervstor(j))*harea
            IF ( Dprst_flag==1 ) landstor = landstor + Dprst_stor_hru(j)*harea
          ELSE
            soilstor = 0.0D0
            snowstor = 0.0D0
            landstor = 0.0D0
            ! wrong if multiple HRUs for any lake
!            IF ( Lake_route_flag==1 ) THEN
!              landstor = Lake_vol(Lake_hru_id(j))*12.0D0
!              srq = Lake_outcfs(Lake_hru_id(j))*Cfs2inches
!              ssq = 0.0D0
!            ELSEIF ( Cascade_flag==1 ) THEN
            IF ( Cascade_flag==1 ) THEN
              srq = Hortonian_lakes(j)*harea
              ssq = Lakein_sz(j)*harea
            ELSE
              srq = 0.0D0
              ssq = 0.0D0
            ENDIF
          ENDIF
          Qsub(k) = Qsub(k) + srq + ssq + gwq
          Subinc_interflow(k) = Subinc_interflow(k) + ssq
          Subinc_sroff(k) = Subinc_sroff(k) + srq
          IF ( Model/=0 ) Subinc_gwflow(k) = Subinc_gwflow(k) + gwq
          Subinc_stor(k) = Subinc_stor(k) + soilstor + gwstor + snowstor + landstor
          Sub_area(k) = Sub_area(k) + harea
        ENDIF
      ENDDO

      DO i = 1, Nsub
        IF ( Sub_area(i)<DNEARZERO ) THEN
          PRINT *, 'ERROR, subbasin:', i, ' does not include any HRUs'
          Inputerror_flag = 1
        ELSE
          Subinc_stor(i) = Subinc_stor(i)/Sub_area(i)
        ENDIF
        Sub_inq(i) = Qsub(i)*Cfs_conv
        Subinc_interflow(i) = Subinc_interflow(i)*Cfs_conv
        IF ( Model/=0 ) Subinc_gwflow(i) = Subinc_gwflow(i)*Cfs_conv
        Subinc_sroff(i) = Subinc_sroff(i)*Cfs_conv
        ! water balance off if lake or muskingum routing
      ENDDO

! allow for possible restart
      !get cummulative subbasin flows
      DO j = 1, Nsub
        IF ( Model/=0 ) Sub_gwflow(j) = Subinc_gwflow(j)
        Sub_sroff(j) = Subinc_sroff(j)
        Sub_interflow(j) = Subinc_interflow(j)
        Sub_cfs(j) = Sub_inq(j)
        DO k = 1, Nsub
          IF ( Tree(j,k)/=0 ) THEN
            IF ( Model/=0 ) Sub_gwflow(j) = Sub_gwflow(j) + Subinc_gwflow(k)
            Sub_sroff(j) = Sub_sroff(j) + Subinc_sroff(k)
            Sub_interflow(j) = Sub_interflow(j) + Subinc_interflow(k)
            Sub_cfs(j) = Sub_cfs(j) + Sub_inq(k)
          ENDIF
        ENDDO
        Sub_cms(j) = Sub_cfs(j)*CFS2CMS_CONV
      ENDDO

 9001 FORMAT ('Initial Tree Structure for Internal Subbasins', /, &
     &        ' 1 indicates a node that flows directly into subbasin', /, 5X, 500I3)
 9002 FORMAT (I3, 2X, 500I3)
 9003 FORMAT (/, 'Tree Structure for Internal Subbasins', /, &
     &        ' 1 indicates a node that eventually flows into subbasin', /, 5X, 500I3)

      END FUNCTION subinit

!***********************************************************************
!     subrun - Computes basin streamflow and on-channel reservoir
!                  storage and outflows
!***********************************************************************
      INTEGER FUNCTION subrun()
      USE PRMS_SUBBASIN
      USE PRMS_MODULE, ONLY: Model, Nsub, Cascade_flag, Dprst_flag !, Lake_route_flag
      USE PRMS_BASIN, ONLY: Hru_area_dble, Active_hrus, Hru_route_order, &
     &    Hru_type, CFS2CMS_CONV, Hru_frac_perv !, Lake_hru_id
      USE PRMS_SET_TIME, ONLY: Cfs_conv !, Cfs2inches
      USE PRMS_SNOW, ONLY: Snowcov_area, Snowmelt
      USE PRMS_CLIMATEVARS, ONLY: Hru_ppt, Swrad, Potet, Tminc, Tmaxc, Tavgc, Hru_rain, Hru_snow
      USE PRMS_FLOWVARS, ONLY: Hru_actet, Ssres_flow, Sroff, &
     &    Ssres_stor, Soil_moist, Pkwater_equiv, Gwres_stor, Soil_moist, Soil_moist_max
      USE PRMS_INTCP, ONLY: Hru_intcpstor
      USE PRMS_SRUNOFF, ONLY: Hru_impervstor, Hortonian_lakes, Dprst_stor_hru
!      USE PRMS_SOILZONE, ONLY: Lakein_sz, Soil_moist_frac, Cpr_stor_frac, Recharge
      USE PRMS_SOILZONE, ONLY: Lakein_sz, Recharge, Soil_moist_tot, Soil_zone_max
      USE PRMS_GWFLOW, ONLY: Gwres_flow
!      USE PRMS_LAKE_ROUTE, ONLY: Lake_vol, Lake_outcfs
      IMPLICIT NONE
! Functions
      INTRINSIC :: DBLE
! Local Variables
      INTEGER :: j, jj, k
      DOUBLE PRECISION :: harea, srq, ssq, gwq, dmy, dmy1, subarea
      DOUBLE PRECISION :: soilstor, snowstor, landstor, dmy2
!     DOUBLE PRECISION :: conv
!***********************************************************************
      subrun = 0

!   Daily time step.
!   Compute reservoir routing and basin outflow

      !rsr, not getting groundwater storage and flow for GSFLOW mode

      DO j = 1, Nsub
        Qsub(j) = 0.0D0
        Sub_inq(j) = 0.0D0
        Subinc_interflow(j) = 0.0D0
        Subinc_sroff(j) = 0.0D0
        Subinc_precip(j) = 0.0D0
        Subinc_rain(j) = 0.0D0
        Subinc_snow(j) = 0.0D0
        Subinc_snowmelt(j) = 0.0D0
        Subinc_pkweqv(j) = 0.0D0
        Subinc_actet(j) = 0.0D0
        Subinc_snowcov(j) = 0.0D0
        Subinc_swrad(j) = 0.0D0
        Subinc_tminc(j) = 0.0D0
        Subinc_tmaxc(j) = 0.0D0
        Subinc_tavgc(j) = 0.0D0
        Subinc_potet(j) = 0.0D0
        Subinc_wb(j) = 0.0D0
        Subinc_recharge(j) = 0.0D0
        Subinc_szstor_frac(j) = 0.0D0
        Subinc_capstor_frac(j) = 0.0D0
      ENDDO
      IF ( Model/=0 ) Subinc_gwflow = 0.0D0

      Laststor = Subinc_stor
      Subinc_stor = 0.0D0

      DO jj = 1, Active_hrus
        j = Hru_route_order(jj)
        ! k indicates which HRU is in which subbasin
        k = Hru_subbasin(j)
        IF ( k>0 ) THEN
          harea = Hru_area_dble(j)
          IF ( Hru_type(j)/=2 ) THEN
            srq = DBLE(Sroff(j))*harea
            ssq = DBLE(Ssres_flow(j))*harea
            soilstor = DBLE(Soil_moist(j)*Hru_frac_perv(j) + Ssres_stor(j))*harea
            snowstor = Pkwater_equiv(j)*harea
            landstor = DBLE(Hru_intcpstor(j)+Hru_impervstor(j))*harea
            IF ( Dprst_flag==1 ) landstor = landstor + Dprst_stor_hru(j)*harea
          ELSE
            soilstor = 0.0D0
            snowstor = 0.0D0
            landstor = 0.0D0
            ! wrong if multiple HRUs for any lake
!            IF ( Lake_route_flag==1 ) THEN
!              landstor = Lake_vol(Lake_hru_id(j))*12.0D0
!              srq = Lake_outcfs(Lake_hru_id(j))*Cfs2inches
!              ssq = 0.0D0
!            ELSEIF ( Cascade_flag==1 ) THEN
            IF ( Cascade_flag==1 ) THEN
              srq = Hortonian_lakes(j)*harea
              ssq = Lakein_sz(j)*harea
            ELSE
              srq = 0.0D0
              ssq = 0.0D0
            ENDIF
          ENDIF
          Qsub(k) = Qsub(k) + srq + ssq
          Subinc_interflow(k) = Subinc_interflow(k) + ssq
          Subinc_sroff(k) = Subinc_sroff(k) + srq
          Subinc_precip(k) = Subinc_precip(k) + DBLE(Hru_ppt(j))*harea
          Subinc_rain(k) = Subinc_rain(k) + DBLE(Hru_rain(j))*harea
          Subinc_snow(k) = Subinc_snow(k) + DBLE(Hru_snow(j))*harea
          Subinc_actet(k) = Subinc_actet(k) + DBLE(Hru_actet(j))*harea
          Subinc_snowmelt(k) = Subinc_snowmelt(k) + DBLE(Snowmelt(j))*harea
          Subinc_pkweqv(k) = Subinc_pkweqv(k) + Pkwater_equiv(j)*harea
          Subinc_snowcov(k) = Subinc_snowcov(k) + DBLE(Snowcov_area(j))*harea
          Subinc_potet(k) = Subinc_potet(k) + DBLE(Potet(j))*harea
          Subinc_swrad(k) = Subinc_swrad(k) + DBLE(Swrad(j))*harea
          Subinc_tminc(k) = Subinc_tminc(k) + DBLE(Tminc(j))*harea
          Subinc_tmaxc(k) = Subinc_tmaxc(k) + DBLE(Tmaxc(j))*harea
          Subinc_tavgc(k) = Subinc_tavgc(k) + DBLE(Tavgc(j))*harea
          Subinc_recharge(k) = Subinc_recharge(k) + Recharge(j)*harea
!          Subinc_szstor_frac(k) = Subinc_szstor_frac(k) + Soil_moist_frac(j)*harea
!          Subinc_capstor_frac(k) = Subinc_capstor_frac(k) + Cpr_stor_frac(k)*harea
          Subinc_szstor_frac(k) = Subinc_szstor_frac(k) + Soil_moist_tot(j)/Soil_zone_max(j)*harea
          Subinc_capstor_frac(k) = Subinc_capstor_frac(k) + Soil_moist(j)/Soil_moist_max(j)*harea
          Subinc_stor(k) = Subinc_stor(k) + soilstor + snowstor + landstor
          IF ( Model/=0 ) THEN
            gwq = DBLE(Gwres_flow(j))*harea
            Qsub(k) = Qsub(k) + gwq
            Subinc_gwflow(k) = Subinc_gwflow(k) + gwq
            Subinc_stor(k) = Subinc_stor(k) + Gwres_stor(j)*harea
          ENDIF
        ENDIF
      ENDDO

      !convert first as subbasins don't have to be in order
      dmy1 = 0.0D0
      DO j = 1, Nsub
        subarea = Sub_area(j)
        Sub_inq(j) = Qsub(j)*Cfs_conv
        dmy = Subinc_interflow(j)/subarea
        dmy2 = Subinc_sroff(j)/subarea
        Subinc_interflow(j) = Subinc_interflow(j)*Cfs_conv
        Subinc_sroff(j) = Subinc_sroff(j)*Cfs_conv
        Subinc_precip(j) = Subinc_precip(j)/subarea
        Subinc_rain(j) = Subinc_rain(j)/subarea
        Subinc_snow(j) = Subinc_snow(j)/subarea
        Subinc_actet(j) = Subinc_actet(j)/subarea
        Subinc_snowmelt(j) = Subinc_snowmelt(j)/subarea
        Subinc_pkweqv(j) = Subinc_pkweqv(j)/subarea
        Subinc_snowcov(j) = Subinc_snowcov(j)/subarea
        Subinc_potet(j) = Subinc_potet(j)/subarea
        Subinc_swrad(j) = Subinc_swrad(j)/subarea
        Subinc_tminc(j) = Subinc_tminc(j)/subarea
        Subinc_tmaxc(j) = Subinc_tmaxc(j)/subarea
        Subinc_tavgc(j) = Subinc_tavgc(j)/subarea
        Subinc_stor(j) = Subinc_stor(j)/subarea
        Subinc_recharge(j) = Subinc_recharge(j)/subarea
        Subinc_szstor_frac(j) = Subinc_szstor_frac(j)/subarea
        Subinc_capstor_frac(j) = Subinc_capstor_frac(j)/subarea
        Subinc_deltastor(j) = Laststor(j) - Subinc_stor(j)
        IF ( Model/=0 ) THEN
          dmy1 = Subinc_gwflow(j)/subarea
          Subinc_gwflow(j) = Subinc_gwflow(j)*Cfs_conv
        ENDIF
        ! water balance off if lake or muskingum routing
        Subinc_wb(j) = Subinc_precip(j) - Subinc_actet(j) - &
     &                 dmy - dmy1 - dmy2 + Subinc_deltastor(j)
      ENDDO

      !get cummulative subbasin flows
      IF ( Model/=0 ) Sub_gwflow = Subinc_gwflow
      DO j = 1, Nsub
        Sub_sroff(j) = Subinc_sroff(j)
        Sub_interflow(j) = Subinc_interflow(j)
        DO k = 1, Nsub
          IF ( Tree(j,k)/=0 ) THEN
            IF ( Model/=0 ) Sub_gwflow(j) = Sub_gwflow(j) + Subinc_gwflow(k)
            Sub_sroff(j) = Sub_sroff(j) + Subinc_sroff(k)
            Sub_interflow(j) = Sub_interflow(j) + Subinc_interflow(k)
          ENDIF
        ENDDO
      ENDDO

      DO j = 1, Nsub
        Sub_cfs(j) = Sub_inq(j)
        DO k = 1, Nsub
          IF ( Tree(j,k)/=0 ) Sub_cfs(j) = Sub_cfs(j) + Sub_inq(k)
        ENDDO
        Sub_cms(j) = Sub_cfs(j)*CFS2CMS_CONV
      ENDDO

      END FUNCTION subrun

max_stars_repo_name: PJ-Thompson/GSFLOWurban

max_issues_repo_name: PJ-Thompson/GSFLOWurban

max_forks_repo_name: PJ-Thompson/GSFLOWurban