      subroutine basin_output

      use time_module
      use hydrograph_module
      use calibration_data_module
      use output_landscape_module
      use basin_module

      implicit none

      integer :: ihru     !none      |counter
      integer :: iihru    !          |
      real :: const       !          |constant used for rate, days, etc

        ! summing subbasin output for the basin
        do ihru = 1, sp_ob%hru
          iihru = lsu_elem(ihru)%obtypno
          if (lsu_elem(iihru)%bsn_frac > 1.e-12) then
            const = lsu_elem(iihru)%bsn_frac
            if (lsu_elem(iihru)%obtyp == "hru") then
              bwb_d = bwb_d + hwb_d(iihru) * const
              bwb_d%sw_init = bwb_d%sw_init + hwb_d(iihru)%sw_init * const
              bwb_d%sw_final = bwb_d%sw_final + hwb_d(iihru)%sw_final * const
              bnb_d = bnb_d + hnb_d(iihru) * const
              bls_d = bls_d + hls_d(iihru) * const
              bpw_d = bpw_d + hpw_d(iihru) * const
            end if
          end if
        end do

            ! or if it is not routed and not in a subbasin
        do ihru = 1, sp_ob%hru_lte
          iihru = lsu_elem(ihru)%obtypno
          if (lsu_elem(iihru)%bsn_frac > 1.e-12) then
            const = lsu_elem(iihru)%bsn_frac
            if (lsu_elem(iihru)%obtyp == "hlt") then
              !const = lsu_elem(iihru)%bsn_frac
              bwb_d = bwb_d + hltwb_d(iihru) * const
              bwb_d%sw_init = bwb_d%sw_init + hltwb_d(iihru)%sw_init * const
              bwb_d%sw_final = bwb_d%sw_final + hltwb_d(iihru)%sw_final * const
              bnb_d = bnb_d + hltnb_d(iihru) * const
              bls_d = bls_d + hltls_d(iihru) * const
              bpw_d = bpw_d + hltpw_d(iihru) * const
            end if
          end if
        end do

        !! sum monthly variables
        bwb_m = bwb_m + bwb_d
        bnb_m = bnb_m + bnb_d
        bls_m = bls_m + bls_d
        bpw_m = bpw_m + bpw_d

!!!!! daily print - BASIN
         if (pco%day_print == "y" .and. pco%int_day_cur == pco%int_day) then
          if (pco%wb_bsn%d == "y") then
            write (2050,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bwb_d  !! waterbal
            if (pco%csvout == "y") then
              write (2054,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bwb_d  !! waterbal
              bwb_d%sw_init = bwb_d%sw_final
            end if
          end if
          if (pco%nb_bsn%d == "y") then
            write (2060,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bnb_d  !! nutrient bal
            if (pco%csvout == "y") then
              write (2064,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bnb_d  !! nutrient bal
            end if
          end if
          if (pco%ls_bsn%d == "y") then
            write (2070,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bls_d  !! losses
            if (pco%csvout == "y") then
              write (2074,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bls_d  !! losses
            end if
          end if
          if (pco%pw_bsn%d == "y") then
            write (2080,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bpw_d  !! plant weather
            if (pco%csvout == "y") then
              write (2084,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bpw_d  !! plant weather
            end if
          end if
        end if

      !! zero daily basin outputs
      bwb_d = hwbz
      bnb_d = hnbz
      bls_d = hlsz
      bpw_d = hpwz

!!!!! monthly print - BASIN
        if (time%end_mo == 1) then
          bwb_y = bwb_y + bwb_m
          bnb_y = bnb_y + bnb_m
          bls_y = bls_y + bls_m
          bpw_y = bpw_y + bpw_m

          const = float (ndays(time%mo + 1) - ndays(time%mo))
          bwb_m = bwb_m // const
          bpw_m = bpw_m // const
          if (pco%wb_bsn%m == "y") then
            write (2051,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bwb_m
            if (pco%csvout == "y") then
              write (2055,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bwb_m
              bwb_m%sw_init = bwb_m%sw_final
            end if
          end if
          if (pco%nb_bsn%m == "y") then
            write (2061,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bnb_m
            if (pco%csvout == "y") then
              write (2065,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bnb_m
            end if
          end if
          if (pco%ls_bsn%m == "y") then
            write (2071,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bls_m
            if (pco%csvout == "y") then
              write (2075,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bls_m
            end if
          end if
          if (pco%pw_bsn%m == "y") then
            write (2081,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bpw_m
            if (pco%csvout == "y") then
              write (2085,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bpw_m
            end if
          end if

          bwb_m = hwbz
          bnb_m = hnbz
          bls_m = hlsz
          bpw_m = hpwz
        end if

!!!!! yearly print - BASIN
        if (time%end_yr == 1) then
           bwb_a = bwb_a + bwb_y
           bnb_a = bnb_a + bnb_y
           bls_a = bls_a + bls_y
           bpw_a = bpw_a + bpw_y

           const = time%day_end_yr
           bwb_y = bwb_y // const
           bpw_y = bpw_y // const
           if (pco%wb_bsn%y == "y") then
             write (2052,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bwb_y
             if (pco%csvout == "y") then
                write (2056,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bwb_y
                bwb_y%sw_init = bwb_y%sw_final
             end if
           end if
           if (pco%nb_bsn%y == "y") then
             write (2062,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bnb_y
             if (pco%csvout == "y") then
               write (2066,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bnb_y
             end if
           end if
           if (pco%ls_bsn%y == "y") then
             write (2072,100)time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bls_y
             if (pco%csvout == "y") then
                write (2076,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bls_y
             end if
           end if
           if (pco%pw_bsn%y == "y") then
             write (2082,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bpw_y
             if (pco%csvout == "y") then
               write (2086,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bpw_y
             end if
           end if

!!!!! zero yearly variables
          bwb_y = hwbz
          bnb_y = hnbz
          bls_y = hlsz
          bpw_y = hpwz
        end if

!!!!! average annual print - BASIN
      if (time%end_sim == 1 .and. pco%wb_bsn%a == "y") then
        bwb_a = bwb_a / time%yrs_prt
        bwb_a = bwb_a // time%days_prt
        write (2053,103) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bwb_a, cal_sim
        if (pco%csvout == "y") then
          write (2057,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bwb_a, cal_sim
        end if
        bwb_a = hwbz
      end if
      if (time%end_sim == 1 .and. pco%nb_bsn%a == "y") then
        bnb_a = bnb_a / time%days_prt
        write (2063,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bnb_a
        if (pco%csvout == "y") then
          write (2067,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bnb_a
        end if
        bnb_a = hnbz
      end if
      if (time%end_sim == 1 .and. pco%ls_bsn%a == "y") then
        bls_a = bls_a / time%yrs_prt
        write (2073,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bls_a
        if (pco%csvout == "y") then
          write (2077,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bls_a
        end if
        bls_a = hlsz
      end if
      if (time%end_sim == 1 .and. pco%pw_bsn%a == "y") then
        bpw_a = bpw_a / time%yrs_prt
        bpw_a = bpw_a // time%days_prt
        write (2083,100) time%day, time%mo, time%day_mo, time%yrc, "       1", "       1", bsn%name, bpw_a
        if (pco%csvout == "y") then
          write (2087,'(*(G0.3,:","))') time%day, time%mo, time%day_mo, time%yrc, "       1", "        1", bsn%name, bpw_a
        end if
        bpw_a = hpwz
      end if

      return

100   format (4i6,2a,2x,a16,30f12.3)
103   format (4i6,2x,2a,2x,a16,30f12.3,a)

      end subroutine basin_output


text:
  !=========================================================================================
  !Copyright (c) 2009-2015, The Regents of the University of Massachusetts, Amherst.
  !E. Polizzi research lab
  !All rights reserved.
  !
  !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. Neither the name of the University nor the names of its contributors may be used to endorse or promote
  !    products derived from this software without specific prior written permission.
  !
  !THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``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 AUTHOR 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.
  !==========================================================================================

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!! SPARSE PRIMITIVES !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!! SINGLE PRECISION VERSION

!!!!! List of routines

!saddcsr, caddcsr           :: add two square csr matrices and return a new one (C=alpha*A+beta*B)
!csaddcsr                   :: same as zaddcsr but A, B are real;  alpha, beta, C are complex

!scsr2csr_up,  ccsr2csr_up  :: extract upper triangular part of a full csr matrix
!scsr2csr_low, ccsr2csr_low :: extract lower triangular part of a full csr matrix

!scsrmm, ccsrmm             :: mat-vec multiplication =>  b=alpha*A*x+beta*b
!chcsrmm                    :: mat-vec multiplication =>  b=alpha*A*x+beta*b (A is complex Hermitian if UPLO=L,U)
!sccsrmm                    :: same as zcsrmm but A is real;  alpha, beta, b are complex

!scsr_transpose, ccsr_transpose :: transpose a csr matrix
!ccsr_htranspose                :: transpose conjugate a csr matrix

!scsr_uplo_to_csr, ccsr_uplo_to_csr  :: convert a upper or lower csr matrix into a full csr one
!chcsr_uplo_to_csr                   :: convert a upper or lower (Hermitian) csr matrix into a full csr one

subroutine saddcsr(N,opt,alpha,sa,isa,jsa,beta,sb,isb,jsb,sc,isc,jsc)
  !  Purpose
  !  =======
  !
  !  Addition of two (square) sparse CSR matrices A and B and return a new (square) CSR matrix C
  !  C=alpha*A+beta*B
  !
  !  A,B,C are real
  !
  !  Arguments
  !  =========
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,B,C.  N >= 0.
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs addition, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine computes only values of the array isc of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value ic(n+1) - 1 is the actual
  !         number of the elements in the arrays sc and jsc.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array jsc
  !
  !         if opt=3, the routine has been called previously with the parameter opt=2. On exit it will
  !         return the array sc
  !
  !
  !  alpha  (input) REAL
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa REAL)
  !
  !  beta  (input) REAL
  !
  !  sb,isb,jsb (input) CSR format for the matrix B (sb REAL)
  !
  !  sc,isc,jsc (input/output) CSR format for the matrix C (sc REAL).
  !             see definition of 'opt' to determine input or output.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  real,dimension(*) :: sa,sb,sc
  integer,dimension(*) :: isa,isb,isc,jsa,jsb,jsc
  real :: alpha,beta
  integer :: opt,N

  integer :: i,ka,kb,kc,ja,jb
  logical :: next,row,col,val

  row=.false.
  col=.false.
  val=.false.

  if (opt==0) then
     row=.true.
     col=.true.
     val=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     col=.true.
  elseif (opt==3) then
     val=.true.
  endif

  kc=1
  if (row) isc(1)=kc
  do i=1,N
     ka=isa(i)
     kb=isb(i)

     next=.false.
     do while (.not.next)

        if (ka<=isa(i+1)-1) then
           ja=jsa(ka)
        else
           ja=n+1
        endif

        if (kb<=isb(i+1)-1) then
           jb=jsb(kb)
        else
           jb=n+1
        endif

        if (ja==jb) then
           if (val) sc(kc)=alpha*sa(ka)+beta*sb(kb)
           if (col) jsc(kc)=ja
           ka=ka+1
           kb=kb+1
           kc=kc+1
        elseif (ja<jb) then
           if (val) sc(kc)=alpha*sa(ka)
           if (col) jsc(kc)=ja
           ka=ka+1
           kc=kc+1
        elseif (jb<ja) then
           if (val) sc(kc)=beta*sb(kb)
           if (col) jsc(kc)=jb
           kb=kb+1
           kc=kc+1
        end if

        if ((ka>isa(i+1)-1).and.(kb>isb(i+1)-1)) next=.true.

     end do

     if (row) isc(i+1)=kc

  end do

end subroutine saddcsr

subroutine caddcsr(N,opt,alpha,sa,isa,jsa,beta,sb,isb,jsb,sc,isc,jsc)
  !  Purpose
  !  =======
  !
  !  Addition of two (square) sparse CSR matrices A and B and return a new (square) CSR matrix C
  !  C=alpha*A+beta*B
  !
  !  A,B,C are complex
  !
  !  Arguments
  !  =========
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,B,C.  N >= 0.
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs addition, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine computes only values of the array isc of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value ic(n+1) - 1 is the actual
  !         number of the elements in the arrays sc and jsc.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array jsc
  !
  !         if opt=3, the routine has been called previously with the parameter opt=2. On exit it will
  !         return the array sc
  !
  !
  !  alpha  (input) COMPLEX
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa COMPLEX)
  !
  !  beta  (input) COMPLEX
  !
  !  sb,isb,jsb (input) CSR format for the matrix B (sb COMPLEX)
  !
  !  sc,isc,jsc (input/output) CSR format for the matrix C (sc COMPLEX).
  !             see definition of 'opt' to determine input or output.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  complex ,dimension(*) :: sa,sb,sc
  integer,dimension(*) :: isa,isb,isc,jsa,jsb,jsc
  complex  :: alpha,beta
  integer :: opt,N
  integer :: i,ka,kb,kc,ja,jb
  logical :: next,row,col,val

  row=.false.
  col=.false.
  val=.false.

  if (opt==0) then
     row=.true.
     col=.true.
     val=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     col=.true.
  elseif (opt==3) then
     val=.true.
  endif

  kc=1
  if (row) isc(1)=kc

  do i=1,N
     ka=isa(i)
     kb=isb(i)

     next=.false.
     do while (.not.next)

        if (ka<=isa(i+1)-1) then
           ja=jsa(ka)
        else
           ja=n+1
        endif

        if (kb<=isb(i+1)-1) then
           jb=jsb(kb)
        else
           jb=n+1
        endif

        if (ja==jb) then
           if (val) sc(kc)=alpha*sa(ka)+beta*sb(kb)
           if (col) jsc(kc)=ja
           ka=ka+1
           kb=kb+1
           kc=kc+1
        elseif (ja<jb) then
           if (val) sc(kc)=alpha*sa(ka)
           if (col) jsc(kc)=ja
           ka=ka+1
           kc=kc+1
        elseif (jb<ja) then
           if (val) sc(kc)=beta*sb(kb)
           if (col) jsc(kc)=jb
           kb=kb+1
           kc=kc+1
        end if

        if ((ka>isa(i+1)-1).and.(kb>isb(i+1)-1)) next=.true.

     end do

     if (row) isc(i+1)=kc

  end do

end subroutine caddcsr

subroutine csaddcsr(N,opt,alpha,sa,isa,jsa,beta,sb,isb,jsb,sc,isc,jsc)
  !  Purpose
  !  =======
  !
  !  Addition of two (square) sparse CSR matrices A and B and return a new (square) CSR matrix C
  !  C=alpha*A+beta*B
  !
  !  A is complex, B and C are real, (alpha,beta are COMPLEX)
  !
  !  Arguments
  !  =========
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,B,C.  N >= 0.
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs addition, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine computes only values of the array isc of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value ic(n+1) - 1 is the actual
  !         number of the elements in the arrays sc and jsc.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array jsc
  !
  !         if opt=3, the routine has been called previously with the parameter opt=2. On exit it will
  !         return the array sc
  !
  !
  !  alpha  (input) COMPLEX
  !
  !  sa,isa,jsa (input) CSR format for the matrix A (sa REAL)
  !
  !  beta  (input) COMPLEX
  !
  !  sb,isb,jsb (input) CSR format for the matrix B (sb REAL)
  !
  !  sc,isc,jsc (input/output) CSR format for the matrix C (sc COMPLEX).
  !             see definition of 'opt' to determine input or output.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  real,dimension(*) :: sa,sb
  complex ,dimension(*) :: sc
  integer,dimension(*) :: isa,isb,isc,jsa,jsb,jsc
  complex  :: alpha,beta
  integer :: opt,N
  integer :: i,ka,kb,kc,ja,jb
  logical :: next,row,col,val

  row=.false.
  col=.false.
  val=.false.

  if (opt==0) then
     row=.true.
     col=.true.
     val=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     col=.true.
  elseif (opt==3) then
     val=.true.
  endif

  kc=1
  if (row) isc(1)=kc

  do i=1,N
     ka=isa(i)
     kb=isb(i)

     next=.false.
     do while (.not.next)

        if (ka<=isa(i+1)-1) then
           ja=jsa(ka)
        else
           ja=n+1
        endif

        if (kb<=isb(i+1)-1) then
           jb=jsb(kb)
        else
           jb=n+1
        endif

        if (ja==jb) then
           if (val) sc(kc)=alpha*sa(ka)+beta*sb(kb)
           if (col) jsc(kc)=ja
           ka=ka+1
           kb=kb+1
           kc=kc+1
        elseif (ja<jb) then
           if (val) sc(kc)=alpha*sa(ka)
           if (col) jsc(kc)=ja
           ka=ka+1
           kc=kc+1
        elseif (jb<ja) then
           if (val) sc(kc)=beta*sb(kb)
           if (col) jsc(kc)=jb
           kb=kb+1
           kc=kc+1
        end if

        if ((ka>isa(i+1)-1).and.(kb>isb(i+1)-1)) next=.true.

     end do

     if (row) isc(i+1)=kc

  end do

end subroutine csaddcsr

subroutine scsr2csr_up(opt,N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine extract the upper triangular part (including diagonal) of a real CSR matrix A
  !      into a new CSR matrix sA
  !      the size of sa, and sja for csr upper format may be overestimated
  !
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs extraction, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine returns only values of the array sia of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value sia(n+1) - 1 is the actual
  !         number of the elements in the arrays sa and sja.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array sa and sja
  !
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a REAL)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa REAL)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  real,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  real,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k
  logical :: row,colval
!!!!
  row=.false.
  colval=.false.
  if (opt==0) then
     row=.true.
     colval=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     colval=.true.
  end if

  k=0
  if (row) sia(1)=1
  do i=1,n
     if (row) sia(i+1)=sia(i)
     do j=ia(i),ia(i+1)-1
        if (ja(j)>=i) then
           k=k+1
           if (row) sia(i+1)=sia(i+1)+1
           if (colval) then
              sja(k)=ja(j)
              sa(k)=a(j)
           end if
        endif
     end do
  end do
end subroutine scsr2csr_up

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine ccsr2csr_up(opt,N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine extract the upper triangular part (including diagonal) of a real CSR matrix A
  !      into a new CSR matrix sA
  !      the size of sa, and sja for csr upper format may be overestimated
  !
  !
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs extraction, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine returns only values of the array sia of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value sia(n+1) - 1 is the actual
  !         number of the elements in the arrays sa and sja.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array sa and sja
  !
  !
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  complex ,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex ,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k
  logical :: row,colval
!!!!
  row=.false.
  colval=.false.
  if (opt==0) then
     row=.true.
     colval=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     colval=.true.
  end if

  k=0
  if (row) sia(1)=1
  do i=1,n
     if (row) sia(i+1)=sia(i)
     do j=ia(i),ia(i+1)-1
        if (ja(j)>=i) then
           k=k+1
           if (row) sia(i+1)=sia(i+1)+1
           if (colval) then
              sja(k)=ja(j)
              sa(k)=a(j)
           end if
        endif
     end do
  end do
end subroutine ccsr2csr_up

subroutine scsr2csr_low(opt,N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine extract the lower triangular part (including diagonal) of a real CSR matrix A
  !      into a new CSR matrix sA
  !      the size of sa, and sja for csr lower format may be overestimated
  !
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs extraction, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine returns only values of the array sia of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value sia(n+1) - 1 is the actual
  !         number of the elements in the arrays sa and sja.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array sa and sja
  !
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a REAL)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa REAL)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  real,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  real,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k
  logical :: row,colval
!!!!
  row=.false.
  colval=.false.
  if (opt==0) then
     row=.true.
     colval=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     colval=.true.
  end if

  k=0
  if (row) sia(1)=1
  do i=1,n
     if (row) sia(i+1)=sia(i)
     do j=ia(i),ia(i+1)-1
        if (ja(j)<=i) then
           k=k+1
           if (row) sia(i+1)=sia(i+1)+1
           if (colval) then
              sja(k)=ja(j)
              sa(k)=a(j)
           end if
        endif
     end do
  end do
end subroutine scsr2csr_low
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine ccsr2csr_low(opt,N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine extract the lower triangular part (including diagonal) of a COMPLEX CSR matrix A
  !      into a new CSR matrix sA
  !      the size of sa, and sja for csr lower format may be overestimated
  !
  !
  !  opt    (input) INTEGER
  !         If opt=0, the routine performs extraction, the memory for the output arrays ic, jc, c
  !         must be allocated beforehand.
  !
  !         if opt=1 the routine returns only values of the array sia of length n + 1, the memory
  !         for this array must be allocated beforehand. On exit the value sia(n+1) - 1 is the actual
  !         number of the elements in the arrays sa and sja.
  !
  !         if opt=2, the routine has been called previously with the parameter opt=1. On exit it will
  !         return the array sa and sja
  !
  !
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: opt,N
  complex ,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex ,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k
  logical :: row,colval
!!!!
  row=.false.
  colval=.false.
  if (opt==0) then
     row=.true.
     colval=.true.
  elseif (opt==1) then
     row=.true.
  elseif (opt==2) then
     colval=.true.
  end if

  k=0
  if (row) sia(1)=1
  do i=1,n
     if (row) sia(i+1)=sia(i)
     do j=ia(i),ia(i+1)-1
        if (ja(j)<=i) then
           k=k+1
           if (row) sia(i+1)=sia(i+1)+1
           if (colval) then
              sja(k)=ja(j)
              sa(k)=a(j)
           end if
        endif
     end do
  end do
end subroutine ccsr2csr_low

SUBROUTINE scsrmm(UPLO,TRANS,N,M,rhs,alpha,a,ia,ja,X,beta,B)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine peforms
  !
  !      b=alpha*o(A)*x+beta*b
  !
  !      where the NxM matrix A is stored in CSR format (a,ia,ja),
  !      and X is Mxrhs and B is Nxrhs
  !
  !      REAL SINGLE PRECISION version
  !
  !      Arguments:
  !
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !       if UPLO='U' only the upper part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !
  !  TRANS (input)  CHARACTER*1
  !       if TRANS='N' o(A)=A -- normal mode
  !       if TRANS='T' o(A)=A^T -- transpose mode
  !
  !  N    (input) INTEGER
  !        The number of row of the matrix A and row of matrix B.  N >= 0.
  !  M    (input) INTEGER
  !        The number of column of matrix A and row of matrix X. M>=0; M=N (square matrix) for UPLO=L,U
  !  rhs  (input) INTEGER
  !        The number of column of matrices B and C
  !
  !  alpha (input) SINGLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a SINGLE PRECISION)
  !
  !  X     (input) DOUBLE PRECISION
  !        matrix of size (Mxrhs)
  !
  !  beta  (input) SINGLE PRECISION
  !
  !  B     (input/output) SINGLE PRECISION
  !        matrix of size (Nxrhs). On exit contains the solution.
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  character(len=1) :: UPLO,TRANS
  integer :: N,M,rhs
  real :: alpha,beta
  real,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  real,dimension(M,*):: x
  real,dimension(N,*) ::b
!!!!!!!!

  integer ::i,j,k
  real,parameter :: SZERO=0.0e0

!!!!!!!! Initialization

if (beta/=SZERO) then
call SSCAL(N*rhs,b,beta,1)
else
b(1:N,1:rhs)=SZERO
end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (UPLO=='F') then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     if (TRANS=='N') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
           end do
        end do

     elseif (TRANS=='T') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

     end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  elseif ((UPLO=='L').or.(UPLO=='U')) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
              if (j/=i) b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

  end if

end SUBROUTINE scsrmm

SUBROUTINE ccsrmm(UPLO,TRANS,N,M,rhs,alpha,a,ia,ja,X,beta,B)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine peforms
  !
  !      b=alpha*o(A)*x+beta*b
  !
  !      where the NxM matrix A is stored in CSR format (a,ia,ja),
  !      and X is Mxrhs and B is Nxrhs
  !
  !      COMPLEX SINGLE PRECISION version
  !
  !
  !      Arguments:
  !
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !       if UPLO='U' only the upper part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !
  !  TRANS (input)  CHARACTER*1
  !       if TRANS='N' o(A)=A   -- normal mode
  !       if TRANS='T' o(A)=A^T -- transpose mode
  !       if TRANS='C' o(A)=A^H -- transpose conjugate mode
  !
  !  N    (input) INTEGER
  !        The number of row of the matrix A and row of matrix B.  N >= 0.
  !  M    (input) INTEGER
  !        The number of column of matrix A and row of matrix X. M>=0; M=N (square matrix) for UPLO=L,U
  !  rhs  (input) INTEGER
  !        The number of column of matrices B and C
  !
  !  alpha (input) SINGLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a SINGLE PRECISION)
  !
  !  X     (input) SINGLE PRECISION
  !        matrix of size (Mxrhs)
  !
  !  beta  (input) SINGLE PRECISION
  !
  !  B     (input/output) SINGLE PRECISION
  !        matrix of size (Nxrhs). On exit contains the solution.
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  character(len=1) :: UPLO,TRANS
  integer :: N,M,rhs
  complex  :: alpha,beta
  complex ,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex ,dimension(M,*):: x
  complex ,dimension(N,*) ::b
!!!!!!!!

  integer ::i,j,k
  complex,parameter :: ZEROC=(0.0e0,0.0e0)

!!!!!!!! Initialization

if (beta/=ZEROC) then
call CSCAL(N*rhs,b,beta,1)
else
b(1:N,1:rhs)=ZEROC
end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (UPLO=='F') then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     if (TRANS=='N') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
           end do
        end do

     elseif (TRANS=='T') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

     elseif (TRANS=='C') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*conjg(a(k))*x(i,1:rhs)
           end do
        end do

     end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  elseif ((UPLO=='L').or.(UPLO=='U')) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

if ((TRANS=='N').or.(TRANS=='T')) then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
              if (j/=i) b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

 elseif (TRANS=='C') then

         do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*conjg(a(k))*x(i,1:rhs)
              if (j/=i) b(i,1:rhs)=b(i,1:rhs)+alpha*conjg(a(k))*x(j,1:rhs)
           end do
        end do

     end if

  end if

END SUBROUTINE ccsrmm

SUBROUTINE chcsrmm(UPLO,TRANS,N,M,rhs,alpha,a,ia,ja,X,beta,B)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine peforms
  !
  !      b=alpha*o(A)*x+beta*b
  !
  !      where the NxM matrix A is stored in CSR format (a,ia,ja),
  !      and X is Mxrhs and B is Nxrhs
  !
  !      Note: A is complex Hermitian if UPLO='L or UPLO='U'
  !
  !      Arguments:
  !
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided- symmetric/hermitian CSR (N must be equal to M)
  !       if UPLO='U' only the upper part of the matrix A is provided- symmetric/hermitian CSR (N must be equal to M)
  !
  !  TRANS (input)  CHARACTER*1
  !       if TRANS='N' o(A)=A   -- normal mode
  !       if TRANS='T' o(A)=A^T -- transpose mode
  !       if TRANS='C' o(A)=A^H -- transpose conjugate mode
  !
  !  N    (input) INTEGER
  !        The number of row of the matrix A and row of matrix B.  N >= 0.
  !  M    (input) INTEGER
  !        The number of column of matrix A and row of matrix X. M>=0; M=N (square matrix) for UPLO=L,U
  !  rhs  (input) INTEGER
  !        The number of column of matrices B and C
  !
  !  alpha (input) COMPLEX SINGLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX SINGLE PRECISION)
  !
  !  X     (input) COMPLEX SINGLE PRECISION
  !        matrix of size (Mxrhs)
  !
  !  beta  (input) COMPLEX SINGLE PRECISION
  !
  !  B     (input/output) COMPLEX SINGLE PRECISION
  !        matrix of size (Nxrhs). On exit contains the solution.
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  character(len=1) :: UPLO,TRANS
  integer :: N,M,rhs
  complex :: alpha,beta
  complex,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex,dimension(M,*) :: x
  complex,dimension(N,*) :: b
!!!!!!!!

  integer ::i,j,k
  real,parameter :: SZERO=0.0e0
  complex,parameter :: ZEROC=(SZERO,SZERO)

!!!!!!!! Initialization

if (beta/=ZEROC) then
call CSCAL(N*rhs,b,beta,1)
else
b(1:N,1:rhs)=ZEROC
end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (UPLO=='F') then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     if (TRANS=='N') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
           end do
        end do

     elseif (TRANS=='T') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

     elseif (TRANS=='C') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*conjg(a(k))*x(i,1:rhs)
           end do
        end do

     end if

  elseif ((UPLO=='L').or.(UPLO=='U')) then

if ((TRANS=='N').or.(TRANS=='C')) then

     do i=1,N
        do k=ia(i),ia(i+1)-1
           j=ja(k)
           b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
   if (j/=i) b(j,1:rhs)=b(j,1:rhs)+alpha*conjg(a(k))*x(i,1:rhs)
        end do
     end do

elseif (TRANS=='T') then

         do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*conjg(a(k))*x(i,1:rhs)
              if (j/=i) b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
           end do
        end do

     end if

  end if

end SUBROUTINE chcsrmm

SUBROUTINE sccsrmm(UPLO,TRANS,N,M,rhs,alpha,a,ia,ja,X,beta,B)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine peforms
  !
  !      b=alpha*o(A)*x+beta*b
  !
  !      where the NxM matrix A is stored in CSR format (a,ia,ja),
  !      and X is Mxrhs and B is Nxrhs
  !
  !      SINGLE PRECISION version
  !
  !      Note: A is real
  !
  !      Arguments:
  !
  !  UPLO (input) CHARACTER*1
  !       if UPLO='F' all the matrix A is provided in CSR format ('F'ull format)
  !       if UPLO='L' only the lower part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !       if UPLO='U' only the upper part of the matrix A is provided- symmetric CSR (N must be equal to M)
  !
  !  TRANS (input)  CHARACTER*1
  !       if TRANS='N' o(A)=A   -- normal mode
  !       if TRANS='T' o(A)=A^T -- transpose mode
  !       if TRANS='C' o(A)=A^H -- transpose conjugate mode
  !
  !  N    (input) INTEGER
  !        The number of row of the matrix A and row of matrix B.  N >= 0.
  !  M    (input) INTEGER
  !        The number of column of matrix A and row of matrix X. M>=0; M=N (square matrix) for UPLO=L,U
  !  rhs  (input) INTEGER
  !        The number of column of matrices B and C
  !
  !  alpha (input) SINGLE PRECISION
  !
  !  a,ia,ja (input) CSR format for the matrix A (a SINGLE PRECISION)
  !
  !  X     (input) SINGLE PRECISION
  !        matrix of size (Mxrhs)
  !
  !  beta  (input) SINGLE PRECISION
  !
  !  B     (input/output) SINGLE PRECISION
  !        matrix of size (Nxrhs). On exit contains the solution.
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  character(len=1) :: UPLO,TRANS
  integer :: N,M,rhs
  complex :: alpha,beta
  real,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex,dimension(M,*):: x
  complex,dimension(N,*) ::b
!!!!!!!!

  integer ::i,j,k
  complex,parameter :: ZEROC=(0.0e0,0.0e0)

!!!!!!!! Initialization

if (beta/=ZEROC) then
call CSCAL(N*rhs,b,beta,1)
else
b(1:N,1:rhs)=ZEROC
end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (UPLO=='F') then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     if (TRANS=='N') then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
           end do
        end do

     elseif ((TRANS=='T').or.(TRANS=='C')) then

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

     end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  elseif ((UPLO=='L').or.(UPLO=='U')) then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        do i=1,N
           do k=ia(i),ia(i+1)-1
              j=ja(k)
              b(i,1:rhs)=b(i,1:rhs)+alpha*a(k)*x(j,1:rhs)
              if (j/=i) b(j,1:rhs)=b(j,1:rhs)+alpha*a(k)*x(i,1:rhs)
           end do
        end do

  end if

END SUBROUTINE sccsrmm

subroutine scsr_transpose(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine transpose a  real CSR matrix A
  !      into a new CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a REAL)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa REAL)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: N
  real,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  real,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k,istart
  integer,dimension(1:N) :: tab
!!!!

  !! #of elements by row
  tab=0
  istart=N
  do i=1,N
     do j=ia(i),ia(i+1)-1
        tab(ja(j))=tab(ja(j))+1
        if (ja(j)<istart) istart=ja(j)
     end do
  end do

  !! find sia
  do i=1,istart-1
     sia(i)=0
  end do
  sia(istart)=1
  do i=istart+1,N+1
     sia(i)=sia(i-1)+tab(i-1)
  end do

  tab=sia(1:N)
  !! find jsa,sa
  do i=1,N
     do k=ia(i),ia(i+1)-1
        sja(tab(ja(k)))=i
        sa(tab(ja(k)))=a(k)
        tab(ja(k))=tab(ja(k))+1
     end do
  end do

end subroutine scsr_transpose

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine ccsr_transpose(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine performs the transpose conjugate of
  !      a double complex CSR matrix A into a new CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: N
  complex,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k,istart
  integer,dimension(1:N) :: tab
!!!!

  !! #of elements by row
  tab=0
  istart=N
  do i=1,N
     do j=ia(i),ia(i+1)-1
        tab(ja(j))=tab(ja(j))+1
        if (ja(j)<istart) istart=ja(j)
     end do
  end do

  !! find sia
  do i=1,istart-1
     sia(i)=0
  end do
  sia(istart)=1
  do i=istart+1,N+1
     sia(i)=sia(i-1)+tab(i-1)
  end do

  tab=sia(1:N)
  !! find jsa,sa
  do i=1,N
     do k=ia(i),ia(i+1)-1
        sja(tab(ja(k)))=i
  !      if (ja(k)/=i) then
           sa(tab(ja(k)))=a(k)
  !      else
  !         sa(tab(ja(k)))=a(k)
  !      end if
        tab(ja(k))=tab(ja(k))+1
     end do
  end do

end subroutine ccsr_transpose

subroutine ccsr_htranspose(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine performs the transpose conjugate of
  !      a  complex CSR matrix A into a new CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  integer :: N
  complex ,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex ,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,j,k,istart
  integer,dimension(1:N) :: tab
!!!!

  !! #of elements by row
  tab=0
  istart=N
  do i=1,N
     do j=ia(i),ia(i+1)-1
        tab(ja(j))=tab(ja(j))+1
        if (ja(j)<istart) istart=ja(j)
     end do
  end do

  !! find sia
  do i=1,istart-1
     sia(i)=0
  end do
  sia(istart)=1
  do i=istart+1,N+1
     sia(i)=sia(i-1)+tab(i-1)
  end do

  tab=sia(1:N)
  !! find jsa,sa
  do i=1,N
     do k=ia(i),ia(i+1)-1
        sja(tab(ja(k)))=i
        if (ja(k)/=i) then
           sa(tab(ja(k)))=conjg(a(k))
        else
           sa(tab(ja(k)))=a(k)
        end if
        tab(ja(k))=tab(ja(k))+1
     end do
  end do

end subroutine ccsr_htranspose

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine scsr_uplo_to_csr(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine convert a real upper or lower CSR matrix A
  !      into a new full CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a REAL)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa REAL)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  include 'f90_noruntime_interface.fi'
  integer :: N
  real,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  real,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,k,nnza,info,opt
  real,dimension(:),pointer ::sb
  integer,dimension(:),pointer :: sib
  integer,dimension(:),pointer :: sjb
  real,parameter :: SONE=1.0e0,SZERO=0.0e0
!!!!
  nnza=ia(n+1)-1

  call wallocate_1s(sb, nnza, info)
  call wallocate_1i(sib, N+1, info)
  call wallocate_1i(sjb, nnza, info)

!!! transpose
  call scsr_transpose(N,a,ia,ja,sb,sib,sjb)
!!! nullify diagonal
  do i=1,N
     do k=sib(i),sib(i+1)-1
        if (sjb(k)==i) sb(k)=SZERO
     end do
  enddo
  !! add the two matrices to form the full matrix
  opt=0
  call saddcsr(N,opt,SONE,a,ia,ja,SONE,sb,sib,sjb,sa,sia,sja)

  call wdeallocate_1s(sb)
  call wdeallocate_1i(sib)
  call wdeallocate_1i(sjb)

end subroutine scsr_uplo_to_csr

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine ccsr_uplo_to_csr(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine convert a complex symetric upper or lower CSR matrix A
  !      into a new full CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  include 'f90_noruntime_interface.fi'
  integer :: N
  complex,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,k,nnza,info,opt
  complex,dimension(:),pointer ::sb
  integer,dimension(:),pointer :: sib
  integer,dimension(:),pointer :: sjb
  real,parameter :: SONE=1.0e0,SZERO=0.0e0
  complex,parameter :: ONEC=(SONE,SZERO),ZEROC=(SZERO,SZERO)
!!!!
  nnza=ia(n+1)-1

  call wallocate_1c(sb, nnza, info)
  call wallocate_1i(sib, N+1, info)
  call wallocate_1i(sjb, nnza, info)

!!! transpose
  call ccsr_transpose(N,a,ia,ja,sb,sib,sjb)
!!! nullify diagonal
  do i=1,N
     do k=sib(i),sib(i+1)-1
        if (sjb(k)==i) sb(k)=ZEROC
     end do
  enddo
  !! add the two matrices to form the full matrix
  opt=0
  call caddcsr(N,opt,ONEC,a,ia,ja,ONEC,sb,sib,sjb,sa,sia,sja)

  call wdeallocate_1c(sb)
  call wdeallocate_1i(sib)
  call wdeallocate_1i(sjb)

end subroutine ccsr_uplo_to_csr

subroutine chcsr_uplo_to_csr(N,a,ia,ja,sa,sia,sja)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !      Purpose: This subroutine convert a complex Hermitian upper or lower CSR matrix A
  !      into a new full CSR matrix sA
  !
  !  N      (input) INTEGER
  !          The number of rows/columns of the matrices A,sA.  N >= 0.
  !
  !  a,ia,ja (input) CSR format for the matrix A (a COMPLEX)
  !
  !  sa,isa,jsa (output) CSR format for the matrix sA (sa COMPLEX)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  implicit none
  include 'f90_noruntime_interface.fi'
  integer :: N
  complex ,dimension(*) ::a
  integer,dimension(*) :: ia
  integer,dimension(*) :: ja
  complex ,dimension(*) ::sa
  integer,dimension(*) :: sia
  integer,dimension(*) :: sja
!!!!
  integer :: i,k,nnza,info,opt
  complex ,dimension(:),pointer ::sb
  integer,dimension(:),pointer :: sib
  integer,dimension(:),pointer :: sjb
  real,parameter :: SONE=1.0e0,SZERO=0.0e0
  complex ,parameter :: ONEC=(SONE,SZERO),ZEROC=(SZERO,SZERO)
!!!!
  nnza=ia(n+1)-1

  call wallocate_1c(sb, nnza, info)
  call wallocate_1i(sib, N+1, info)
  call wallocate_1i(sjb, nnza, info)

!!! transpose
  call ccsr_htranspose(N,a,ia,ja,sb,sib,sjb)
!!! nullify diagonal
  do i=1,N
     do k=sib(i),sib(i+1)-1
        if (sjb(k)==i) sb(k)=ZEROC
     end do
  enddo
  !! add the two matrices to form the full matrix
  opt=0
  call caddcsr(N,opt,ONEC,a,ia,ja,ONEC,sb,sib,sjb,sa,sia,sja)

  call wdeallocate_1c(sb)
  call wdeallocate_1i(sib)
  call wdeallocate_1i(sjb)

end subroutine chcsr_uplo_to_csr

max_stars_repo_name: lcirrott/Kratos

max_issues_repo_name: PFEM/feast

max_forks_repo_name: PFEM/feast


          program demo_gulp
          use M_io,      only : gulp
          implicit none
          character(len=4096)          :: FILENAME   ! file to read
          character(len=:),allocatable :: pageout(:) ! array to hold file in memory
          integer                      :: longest, lines, i
             ! get a filename
             call get_command_argument(1, FILENAME)
             ! allocate character array and copy file into it
             call gulp(FILENAME,pageout)
             if(.not.allocated(pageout))then
                write(*,*)'*demo_gulp* failed to load file '//FILENAME
             else
                ! write file from last line to first line
                longest=len(pageout)
                lines=size(pageout)
                write(*,*)'number of lines is ',lines
                write(*,*)'and length of lines is ',longest
                write(*,'(a)')repeat('%',longest+2)
                write(*,'("%",a,"%")')(trim(pageout(i)),i=lines,1,-1)
                write(*,'(a)')repeat('%',longest+2)
                deallocate(pageout)  ! release memory
             endif
          end program demo_gulp


text:
Describe Users/courtwilson here.

20120402 18:40:55 nbsp Welcome to the Wiki, courtwilson. Please dont remove information like you have been doing with Sutter Medical Group. Thank you. Users/PeterBoulay

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module artificial_viscosity_bc_operator
    use mod_kinds,              only: rk,ik
    use mod_constants,          only: HALF, ONE, TWO, THREE, FIVE
    use type_operator,          only: operator_t
    use type_chidg_worker,      only: chidg_worker_t
    use type_properties,        only: properties_t
    use DNAD_D
    implicit none

    private

    !> Implementation of the artificial viscosity boundary averate flux.
    !!
    !!  - At a boundary interface, the solution states Q- and Q+ exists on opposite
    !!    sides of the boundary. The average flux is computed as Favg = 1/2(F(Q-) + F(Q+))
    !!
    !!  @author Nathan A. Wukie
    !!  @date   01/31/2017
    !!
    !-------------------------------------------------------------------------------------------
    type, extends(operator_t), public :: artificial_viscosity_bc_operator_t

    contains

        procedure   :: init
        procedure   :: compute

    end type artificial_viscosity_bc_operator_t
    !********************************************************************************************

contains

    !>
    !!
    !!  @author Nathan A. Wukie
    !!  @date   01/31/2017
    !!
    !--------------------------------------------------------------------------------------------
    subroutine init(self)
        class(artificial_viscosity_bc_operator_t),   intent(inout) :: self

        !
        ! Set operator name
        !
        call self%set_name('Artificial Viscosity BC Operator')

        !
        ! Set operator type
        !
        call self%set_operator_type('BC Diffusive Operator')

        !
        ! Set operator equations
        !
        call self%add_primary_field('Artificial Viscosity')

    end subroutine init
    !*******************************************************************************************

    !>  Boundary Flux routine for Artificial Viscosity
    !!
    !!  @author Nathan A. Wukie
    !!  @date   01/31/2017
    !!
    !!------------------------------------------------------------------------------------------
    subroutine compute(self,worker,prop)
        class(artificial_viscosity_bc_operator_t),  intent(inout)   :: self
        type(chidg_worker_t),                           intent(inout)   :: worker
        class(properties_t),                            intent(inout)   :: prop

        ! Storage at quadrature nodes
        type(AD_D), allocatable,    dimension(:) ::         &
            lamda,  tau,                                    &
            diffusion_1,    diffusion_2,    diffusion_3,    &
            deps_dx,        deps_dy,        deps_dz,        &
            flux_1,         flux_2,         flux_3,         &
            integrand

        integer(ik) :: order
        real(rk)    :: C1, C2, h(3),  &
                       eta_1, eta_2, eta_3

        !
        ! Interpolate solution to quadrature nodes
        !
        deps_dx = worker%get_field('Artificial Viscosity', 'grad1', 'boundary')
        deps_dy = worker%get_field('Artificial Viscosity', 'grad2', 'boundary')
        deps_dz = worker%get_field('Artificial Viscosity', 'grad3', 'boundary')

        !
        ! Get model field for Maximum Wave Speed
        !
        lamda = worker%get_field('Maximum Wave Speed', 'value', 'boundary')

        !
        ! Get element resolution information
        !
        h     = worker%element_size('interior')
        order = worker%solution_order('interior')

        !
        ! Compute time constant and conductivity
        !
        C1 = THREE
        C2 = FIVE

        tau = minval(h)/(C1 * real(order+1,rk) * lamda)
        eta_1 = C2*(h(1)**TWO)
        eta_2 = C2*(h(2)**TWO)
        eta_3 = C2*(h(3)**TWO)

        !
        ! Compute diffusion coefficient tensor diagonal
        !
        diffusion_1 = eta_1/tau
        diffusion_2 = eta_2/tau
        diffusion_3 = eta_3/tau

        !================================
        !             FLUX
        !================================
        flux_1 = -(diffusion_1 * deps_dx)
        flux_2 = -(diffusion_2 * deps_dy)
        flux_3 = -(diffusion_3 * deps_dz)

        call worker%integrate_boundary_condition('Artificial Viscosity','Diffusion',flux_1,flux_2,flux_3)

    end subroutine compute
    !*********************************************************************************************

end module artificial_viscosity_bc_operator

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


text:
      real *8 function ftime (iparam1, jparam2, kparam3, fm, concv)
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
!***********************************************************************
!D1
!D1 PURPOSE
!D1
!D1 To calculate particle time delay using type curves.
!D1
!***********************************************************************
!D2
!D2 REVISION HISTORY
!D2
!D2 FEHM Version 2.20 [10086-STN-2.20-00]
!D2 Initial implementation: 28-Oct-02, Programmer: Z. Dash
!D2
!D2 $Log:   /pvcs.config/fehm90/src/ftime.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:43:02   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:05:28   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!***********************************************************************
!D3
!D3 REQUIREMENTS TRACEABILITY
!D3
!D3 2.3.5 Cell-based particle-tracking module
!D3 2.3.6 Streamline particle-tracking module
!D3
!***********************************************************************
!D4
!D4 SPECIAL COMMENTS AND REFERENCES
!D4
!D4 Requirements from SDN: 10086-RD-2.20-00
!D4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
!D4   FEHM Application Version 2.20
!D4
!***********************************************************************

      use compfrac
      implicit none
      real*8 concv, cadj, ctime, concv_sc
      integer fm, iparam1, jparam2, kparam3, low, high, midpoint, numv

      numv = nump (iparam1, jparam2, kparam3, fm)
      if (numparams .le. 2) then
! Don't scale for Sudicky-Frind curves
         concv_sc = concv
      else
         concv_sc = concv*conc(iparam1, jparam2, kparam3, fm, numv)
      endif
      if (concv_sc .le. conc(iparam1, jparam2, kparam3, fm, 1)) then
         low = 1
         high = 1
         ftime = dtime(iparam1, jparam2, kparam3, fm, 1)
         if(numparams.le.2) then
            ftime = dlog10(ftime)
         end if
         return
      else if (concv_sc .ge. conc(iparam1, jparam2, kparam3, fm, numv))
     .        then
         low = numv
         high = numv
         ftime = dtime(iparam1, jparam2, kparam3, fm, numv)
         if(numparams.le.2) then
            ftime = dlog10(ftime)
         end if
         return
      else
         low = 1
         high = numv
 20      midpoint = (low + high) / 2
         if (concv_sc .eq.
     .        conc(iparam1, jparam2, kparam3, fm, midpoint))
     .        then
            low = midpoint
            high = low
            ftime = dtime(iparam1, jparam2, kparam3, fm,
     .           midpoint)
            if(numparams.le.2) then
               ftime = dlog10(ftime)
            end if
            return
         else if (concv_sc .lt.
     .           conc(iparam1, jparam2, kparam3, fm, midpoint)) then
            high = midpoint
         else
            low = midpoint
         end if
         if (high - low .eq. 1) goto 30
         goto 20
      endif

 30   continue
      cadj = (concv_sc - conc(iparam1, jparam2, kparam3, fm, low)) /
     .     (conc(iparam1, jparam2, kparam3, fm, high) -
     .     conc(iparam1, jparam2, kparam3, fm, low))
      if(numparams.le.2) then
      ctime = dlog10 (dtime(iparam1, jparam2, kparam3, fm, low)) +
     .     cadj * (dlog10 (dtime(iparam1, jparam2, kparam3, fm, high)) -
     .     dlog10 (dtime(iparam1, jparam2, kparam3, fm, low)))
      ftime = ctime
      else
      ctime = dtime(iparam1, jparam2, kparam3, fm, low) +
     .     cadj * (dtime(iparam1, jparam2, kparam3, fm, high) -
     .     dtime(iparam1, jparam2, kparam3, fm, low))
      ftime = ctime
      end if

      return
      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
      program zerotrip

c     The loop body should not disturb the postcondition

      real t(10)

      if(n.ge.1) then

         k = 3
         do i = 1, n
            k = 2*k
            t(i) = 0
         enddo
         print *, k, i

      endif

      if(n.ge.0) then

         k = 3
         do i = 1, n
            k = 2*k
            t(i) = 0
         enddo
         print *, k, i

      endif

      if(n.lt.0) then

         k = 3
         do i = 1, n
            k = 2*k
            t(i) = 0
         enddo
         print *, k, i

      endif

      end


!
!
!                             AMG4PSBLAS version 1.0
!    Algebraic Multigrid Package
!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
!    (C) Copyright 2021
!
!        Salvatore Filippone
!        Pasqua D'Ambra
!        Fabio Durastante
!
!    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 AMG4PSBLAS 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 AMG4PSBLAS 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: amg_d_base_smoother_mod.f90
!
! Module: amg_d_base_smoother_mod
!
!  This module defines:
!  - the amg_d_base_smoother_type data structure containing the
!    smoother  and related data structures;
!
!  It contains routines for
!  - Building and applying;
!  - checking if the smoother is correctly defined;
!  - printing a	description of the preconditioner;
!  - deallocating the preconditioner data structure.
!
!  What is the difference between a smoother and a solver?
!  In the mathematics literature the two concepts are treated
!  essentially as synonymous, but here we are using them in a more
!  computer-science oriented fashion. In particular, a SMOOTHER object
!  contains a SOLVER object: the SOLVER operates locally within the
!  current process, whereas the SMOOTHER object accounts for (possible)
!  interactions between processes.
!  Some solvers (MUMPS and SuperLU_DIST) can also operate on the entire
!  distributed matrix, in which case the smoother object essentially
!  becomes transparent.
!
module amg_d_base_smoother_mod

  use amg_d_base_solver_mod
  use psb_base_mod, only : psb_desc_type, psb_dspmat_type, psb_epk_,&
       & psb_d_vect_type, psb_d_base_vect_type, psb_d_base_sparse_mat, &
       & psb_dpk_, psb_i_base_vect_type, psb_erractionsave, psb_error_handler

  !
  !
  !
  ! Type: amg_T_base_smoother_type.
  !
  !  It holds the smoother a single level. Its only mandatory component is a solver
  !  object which holds a local solver; this decoupling allows to have the same solver
  !  e.g ILU to work with Jacobi with multiple sweeps as well as with any AS variant.
  !
  !  type  amg_T_base_smoother_type
  !    class(amg_T_base_solver_type), allocatable :: sv
  !  end type amg_T_base_smoother_type
  !
  !   Methods:
  !
  !    build      -   Compute the actual contents of the smoother; includes
  !                   invocation of the build method on the solver component.
  !    free       -   Release memory
  !    apply      -   Apply the smoother to a vector (or to an array); includes
  !                   invocation of the apply method on the solver component.
  !    descr      -   Prints a description of the object.
  !    default    -   Set default values
  !    dump       -   Dump to file object contents
  !    set        -   Sets various parameters; when a request is unknown
  !                   it is passed to the solver object for further processing.
  !    check      -   Sanity checks.
  !    sizeof     -   Total memory occupation in bytes
  !    get_nzeros -   Number of nonzeros
  !    stringval  -   convert string to val for internal parms
  !    get_fmt    -   short string descriptor
  !    get_id     -   numeric id descriptro
  !    get_wrksz  -   How many workspace vector does apply_vect need
  !
  !
  !

  type  amg_d_base_smoother_type
    class(amg_d_base_solver_type), allocatable :: sv
  contains
    procedure, pass(sm) :: apply_v => amg_d_base_smoother_apply_vect
    procedure, pass(sm) :: apply_a => amg_d_base_smoother_apply
    generic, public     :: apply => apply_a, apply_v
    procedure, pass(sm) :: check => amg_d_base_smoother_check
    procedure, pass(sm) :: dump  => amg_d_base_smoother_dmp
    procedure, pass(sm) :: clone => amg_d_base_smoother_clone
    procedure, pass(sm) :: build => amg_d_base_smoother_bld
    procedure, pass(sm) :: cnv   => amg_d_base_smoother_cnv
    procedure, pass(sm) :: free  => amg_d_base_smoother_free
    procedure, pass(sm) :: clone_settings => amg_d_base_smoother_clone_settings
    procedure, pass(sm) :: clear_data     => amg_d_base_smoother_clear_data
    procedure, pass(sm) :: cseti => amg_d_base_smoother_cseti
    procedure, pass(sm) :: csetc => amg_d_base_smoother_csetc
    procedure, pass(sm) :: csetr => amg_d_base_smoother_csetr
    generic, public     :: set   => cseti, csetc, csetr
    procedure, pass(sm) :: default => d_base_smoother_default
    procedure, pass(sm) :: descr =>   amg_d_base_smoother_descr
    procedure, pass(sm) :: sizeof =>  d_base_smoother_sizeof
    procedure, pass(sm) :: get_nzeros => d_base_smoother_get_nzeros
    procedure, pass(sm) :: get_wrksz => d_base_smoother_get_wrksize
    procedure, nopass   :: stringval => amg_stringval
    procedure, nopass   :: get_fmt   => d_base_smoother_get_fmt
    procedure, nopass   :: get_id    => d_base_smoother_get_id
  end type amg_d_base_smoother_type

  private :: d_base_smoother_sizeof, d_base_smoother_get_fmt, &
       &  d_base_smoother_default, d_base_smoother_get_nzeros, &
       & d_base_smoother_get_id, d_base_smoother_get_wrksize

  interface
    subroutine amg_d_base_smoother_apply(alpha,sm,x,beta,y,desc_data,&
         & trans,sweeps,work,info,init,initu)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_
      type(psb_desc_type), intent(in)             :: desc_data
      class(amg_d_base_smoother_type), intent(inout) :: sm
      real(psb_dpk_),intent(inout)                :: x(:)
      real(psb_dpk_),intent(inout)                :: y(:)
      real(psb_dpk_),intent(in)                   :: alpha,beta
      character(len=1),intent(in)                  :: trans
      integer(psb_ipk_), intent(in)                :: sweeps
      real(psb_dpk_),target, intent(inout)        :: work(:)
      integer(psb_ipk_), intent(out)               :: info
      character, intent(in), optional       :: init
      real(psb_dpk_),intent(inout), optional :: initu(:)
    end subroutine amg_d_base_smoother_apply
  end interface

  interface
    subroutine amg_d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
         &  trans,sweeps,work,wv,info,init,initu)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_
      type(psb_desc_type), intent(in)                :: desc_data
      class(amg_d_base_smoother_type), intent(inout) :: sm
      type(psb_d_vect_type),intent(inout)            :: x
      type(psb_d_vect_type),intent(inout)            :: y
      real(psb_dpk_),intent(in)                       :: alpha,beta
      character(len=1),intent(in)                      :: trans
      integer(psb_ipk_), intent(in)                    :: sweeps
      real(psb_dpk_),target, intent(inout)            :: work(:)
      type(psb_d_vect_type),intent(inout)            :: wv(:)
      integer(psb_ipk_), intent(out)                   :: info
      character, intent(in), optional                :: init
      type(psb_d_vect_type),intent(inout), optional   :: initu
    end subroutine amg_d_base_smoother_apply_vect
  end interface

  interface
    subroutine amg_d_base_smoother_check(sm,info)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_
      ! Arguments
      class(amg_d_base_smoother_type), intent(inout) :: sm
      integer(psb_ipk_), intent(out)                   :: info
    end subroutine amg_d_base_smoother_check
  end interface

  interface
    subroutine amg_d_base_smoother_cseti(sm,what,val,info,idx)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_
      ! Arguments
      class(amg_d_base_smoother_type), intent(inout) :: sm
      character(len=*), intent(in)                     :: what
      integer(psb_ipk_), intent(in)                    :: val
      integer(psb_ipk_), intent(out)                   :: info
      integer(psb_ipk_), intent(in), optional          :: idx
    end subroutine amg_d_base_smoother_cseti
  end interface

  interface
    subroutine amg_d_base_smoother_csetc(sm,what,val,info,idx)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_
      class(amg_d_base_smoother_type), intent(inout) :: sm
      character(len=*), intent(in)                     :: what
      character(len=*), intent(in)                     :: val
      integer(psb_ipk_), intent(out)                   :: info
      integer(psb_ipk_), intent(in), optional          :: idx
    end subroutine amg_d_base_smoother_csetc
  end interface

  interface
    subroutine amg_d_base_smoother_csetr(sm,what,val,info,idx)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_
      ! Arguments
      class(amg_d_base_smoother_type), intent(inout) :: sm
      character(len=*), intent(in)                     :: what
      real(psb_dpk_), intent(in)                        :: val
      integer(psb_ipk_), intent(out)                   :: info
      integer(psb_ipk_), intent(in), optional          :: idx
    end subroutine amg_d_base_smoother_csetr
  end interface

  interface
    subroutine amg_d_base_smoother_bld(a,desc_a,sm,info,amold,vmold,imold)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_, psb_i_base_vect_type
      ! Arguments
      type(psb_dspmat_type), intent(in), target     :: a
      Type(psb_desc_type), Intent(inout)              :: desc_a
      class(amg_d_base_smoother_type), intent(inout) :: sm
      integer(psb_ipk_), intent(out)                  :: info
      class(psb_d_base_sparse_mat), intent(in), optional :: amold
      class(psb_d_base_vect_type), intent(in), optional  :: vmold
      class(psb_i_base_vect_type), intent(in), optional  :: imold
    end subroutine amg_d_base_smoother_bld
  end interface

  interface
    subroutine amg_d_base_smoother_cnv(sm,info,amold,vmold,imold)
      import :: psb_d_base_sparse_mat, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_, psb_i_base_vect_type
      ! Arguments
      class(amg_d_base_smoother_type), intent(inout) :: sm
      integer(psb_ipk_), intent(out)                  :: info
      class(psb_d_base_sparse_mat), intent(in), optional :: amold
      class(psb_d_base_vect_type), intent(in), optional  :: vmold
      class(psb_i_base_vect_type), intent(in), optional  :: imold
    end subroutine amg_d_base_smoother_cnv
  end interface

  interface
    subroutine amg_d_base_smoother_free(sm,info)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_
      ! Arguments
      class(amg_d_base_smoother_type), intent(inout) :: sm
      integer(psb_ipk_), intent(out)                   :: info
    end subroutine amg_d_base_smoother_free
  end interface

  interface
    subroutine amg_d_base_smoother_descr(sm,info,iout,coarse)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_
      ! Arguments
      class(amg_d_base_smoother_type), intent(in) :: sm
      integer(psb_ipk_), intent(out)                :: info
      integer(psb_ipk_), intent(in), optional       :: iout
      logical, intent(in), optional                 :: coarse
    end subroutine amg_d_base_smoother_descr
  end interface

  interface
    subroutine amg_d_base_smoother_dmp(sm,desc,level,info,prefix,head,smoother,solver,global_num)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_
      class(amg_d_base_smoother_type), intent(in) :: sm
      type(psb_desc_type), intent(in)               :: desc
      integer(psb_ipk_), intent(in)              :: level
      integer(psb_ipk_), intent(out)             :: info
      character(len=*), intent(in), optional :: prefix, head
      logical, optional, intent(in)    :: smoother, solver, global_num
    end subroutine amg_d_base_smoother_dmp
  end interface

  interface
    subroutine amg_d_base_smoother_clone(sm,smout,info)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_
      Implicit None

      ! Arguments
      class(amg_d_base_smoother_type), intent(inout)              :: sm
      class(amg_d_base_smoother_type), allocatable, intent(inout) :: smout
      integer(psb_ipk_), intent(out)                 :: info
    end subroutine amg_d_base_smoother_clone
  end interface

  interface
    subroutine amg_d_base_smoother_clone_settings(sm,smout,info)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_
      Implicit None

      ! Arguments
      class(amg_d_base_smoother_type), intent(inout) :: sm
      class(amg_d_base_smoother_type), intent(inout) :: smout
      integer(psb_ipk_), intent(out)                   :: info
    end subroutine amg_d_base_smoother_clone_settings
  end interface

  interface
    subroutine amg_d_base_smoother_clear_data(sm,info)
      import :: psb_desc_type, psb_dspmat_type,  psb_d_base_sparse_mat, &
           & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
           & amg_d_base_smoother_type, psb_ipk_
      Implicit None

      ! Arguments
      class(amg_d_base_smoother_type), intent(inout) :: sm
      integer(psb_ipk_), intent(out)                   :: info
    end subroutine amg_d_base_smoother_clear_data
  end interface

contains
  !
  ! Function returning the size of the amg_prec_type data structure
  ! in bytes or in number of nonzeros of the operator(s) involved.
  !

  function d_base_smoother_get_nzeros(sm) result(val)
    implicit none
    class(amg_d_base_smoother_type), intent(in) :: sm
    integer(psb_epk_) :: val
    integer(psb_ipk_)             :: i
    val = 0
    if (allocated(sm%sv)) &
         &  val =  sm%sv%get_nzeros()
  end function d_base_smoother_get_nzeros

  function d_base_smoother_sizeof(sm) result(val)
    implicit none
    ! Arguments
    class(amg_d_base_smoother_type), intent(in) :: sm
    integer(psb_epk_)                    :: val
    integer(psb_ipk_)             :: i

    val = 0
    if (allocated(sm%sv)) then
      val = sm%sv%sizeof()
    end if

    return
  end function d_base_smoother_sizeof

  !
  ! Set sensible defaults.
  ! To be called immediately after allocation
  !
  subroutine d_base_smoother_default(sm)
    implicit none
    ! Arguments
    class(amg_d_base_smoother_type), intent(inout) :: sm
    ! Do nothing for base version

    if (allocated(sm%sv)) call sm%sv%default()

    return
  end subroutine d_base_smoother_default

  function d_base_smoother_get_wrksize(sm) result(val)
    implicit none
    class(amg_d_base_smoother_type), intent(inout) :: sm
    integer(psb_ipk_)  :: val

    val = 0
    if (allocated(sm%sv)) val = val + sm%sv%get_wrksz()

  end function d_base_smoother_get_wrksize

  function d_base_smoother_get_fmt() result(val)
    implicit none
    character(len=32)  :: val

    val = "Base smoother"
  end function d_base_smoother_get_fmt

  function d_base_smoother_get_id() result(val)
    implicit none
    integer(psb_ipk_)  :: val

    val = amg_base_smooth_
  end function d_base_smoother_get_id

end module amg_d_base_smoother_mod


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20090705 19:20:19 nbsp Welcome to the Wiki. So is it Taqueria Guadalajara The Guad or Taqueria Guadalajara Guads? Is it possible that in your circle it is called Guads, but known as The Guad by others? Users/JasonAller

20090707 17:39:12 nbsp hey, on the guad, irregardless of the hullabaloo over pseudonyms, ppl are gonna get peved about reverts blowin out comments and uninvolved edits. Users/BlueDot

20090707 17:44:28 nbsp I am a current student. :) Users/WilliamLewis

20090707 17:48:28 nbsp Anyway, youre talking about changing something that goes back to the http://daviswiki.org/Taqueria_Guadalajara?actionrecall&version3 third revision of the page in 2004. Nobody has seemed to have a problem with it until early 2008. Yes, calling the place the guad is stupid, but my friends do it all the time. It pretty much is the default burrito joint that people in my social circle go to (although I much prefer Taqueria Davis). I personally prefer just guad, but whatever. Users/WilliamLewis

20090707 20:26:29 nbsp I am a current student, and Ive never heard it called Guads. Maybe we should only call it The Guad, then? Users/JoePomidor

20090707 23:56:06 nbsp The only true name is the one on their sign or their business license. All others are, as the page says, also known as and so there is no point in saying one is more correct, they are not the true name. Its like arguing if 2+2 equals 5 or 6. A different situation might be if The Guad or Guads used to be on the sign and they changed their name, but I havent seen any evidence of this. Users/NickSchmalenberger

20091004 13:10:48 nbsp If not that house, what is the oldest structure in Davis? Users/JabberWokky

20100308 13:59:56 nbsp Great post on the protest criticism page! Users/TomGarberson

20100308 14:54:10 nbsp Im on here all the time too.  Crossplatform communications huzzah! Users/OliviaY

20100309 23:51:42 nbsp where did you hear about Trader Joes?  Last time I checked it was a rumor. Users/OliviaY


text:
!
!     CalculiX - A 3-dimensional finite element program
!     Copyright (C) 1998-2021 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine restartshort(nset,nload,nbody,nforc,nboun,nk,ne,
     &     nmpc,nalset,nmat,ntmat_,npmat_,norien,nam,nprint,mi,
     &     ntrans,ncs_,namtot,ncmat_,memmpc_,ne1d,ne2d,nflow,
     &     set,meminset,rmeminset,jobnamec,irestartstep,icntrl,ithermal,
     &     nener,nstate_,ntie,nslavs,nkon,mcs,nprop,mortar,
     &     ifacecount,nintpoint,infree,nef,mpcend,nheading_,network,
     &     nfc,ndc,iprestr)
!
!     istartset := meminset
!     iendset := rmeminset
!
      implicit none
!
      character*80 version
      character*81 set(*)
      character*132 fnrstrt,jobnamec(*)
!
      integer istep,nset,nload,nforc,nboun,nk,ne,nmpc,nalset,
     &     nmat,ntmat_,npmat_,norien,nam,nprint,mi(*),ntrans,ncs_,
     &     nprop,namtot,ncmat_,memmpc_,ne1d,ne2d,nflow,infree(4),mortar,
     &     nmethod,iperturb(2),meminset(*),rmeminset(*),nintpoint,
     &     i,j,k,ipos,icntrl,nener,irestartstep,im0,im1,im2,mem,iact,
     &     istat,nkon,nlabel,iplas,ithermal(*),nstate_,iprestr,
     &     mcs,ntie,nbody,nslavs,ifacecount,iversion,nef,mpcend,
     &     maxlenmpc,nheading_,network,nfc,ndc
!
      if(icntrl.eq.0) then
!
!     this branch is called from readinput.c
!     its purpose is to read the value of nset
!
!     determining the name of the restart file
!
        ipos=index(jobnamec(1),char(0))
        fnrstrt(1:ipos-1)=jobnamec(1)(1:ipos-1)
        fnrstrt(ipos:ipos+3)=".rin"
        do i=ipos+4,132
          fnrstrt(i:i)=' '
        enddo
!
!     opening the restart file
!
        open(15,file=fnrstrt,ACCESS='SEQUENTIAL',FORM='UNFORMATTED',
     &       err=15)
!
        iversion=0
!
        do
!
          read(15,iostat=istat) version
          if(istat.lt.0) then
            if(irestartstep.eq.0) then
!
!     reading the last step
!
              irestartstep=istep
              close(15)
              open(15,file=fnrstrt,ACCESS='SEQUENTIAL',
     &             FORM='UNFORMATTED',err=15)
              read(15) version
            else
              write(*,*) '*ERROR in restartshort: requested step'
              write(*,*) '       is not in the restart file'
              call exit(201)
            endif
          endif
!
          if(iversion.eq.0) then
            write(*,*)
            write(*,*) '*INFO: restart file ',fnrstrt
            write(*,*) '       has been opened for reading.'
            write(*,*) '       it was created with CalculiX ',version
            iversion=1
          endif
!
          read(15)istep
!
!     reading the number of sets
!
          read(15)nset
!
          if(istep.eq.irestartstep) exit
!
          read(15)nalset
!
!     load size
!
          read(15)nload
          read(15)nbody
          read(15)nforc
          read(15)nboun
          read(15)nflow
!
!     mesh size
!
          read(15)nk
          read(15)ne
          read(15)nef
          read(15)nkon
          read(15)(mi(i),i=1,3)
!
!     constraint size
!
          read(15)nmpc
          read(15)mpcend
          read(15)maxlenmpc
          read(15)memmpc_
!
!     coupling, distributing
!
          read(15)nfc
          read(15)ndc
!
!     material size
!
          read(15)nmat
          read(15)ntmat_
          read(15)npmat_
          read(15)ncmat_
!
!     property info
!
          read(15)nprop
!
!     transformation size
!
          read(15)norien
          read(15)ntrans
!
!     amplitude size
!
          read(15)nam
          read(15)namtot
!
!     print size
!
          read(15)nprint
          read(15)nlabel
!
!     tie size
!
          read(15)ntie
!
!     cyclic symmetry size
!
          read(15)ncs_
          read(15)mcs
!
!     1d and 2d element size
!
          read(15)ne1d
          read(15)ne2d
          read(15)(infree(i),i=1,4)
!
!     procedure info
!
          read(15)nmethod
          read(15)(iperturb(i),i=1,2)
          read(15)nener
          read(15)iplas
          read(15)ithermal(1)
          read(15)nstate_
          read(15)nslavs
          read(15)iprestr
          read(15)mortar
          if(mortar.eq.1) then
            read(15)ifacecount
            read(15)nintpoint
          endif
          read(15)nheading_
          read(15)network
!
!     skipping the next entries
!
          call skip(nset,nalset,nload,nbody,
     &         nforc,nboun,nk,ne,nkon,
     &         mi,nmpc,mpcend,nmat,ntmat_,npmat_,ncmat_,norien,
     &         ntrans,nam,nprint,nlabel,ncs_,ne1d,ne2d,infree,
     &         nmethod,iperturb,nener,ithermal,nstate_,iprestr,
     &         mcs,ntie,nslavs,nprop,mortar,ifacecount,nintpoint,
     &         nef,nheading_,nfc,ndc)
!
        enddo
!
        close(15)
!
        return
      endif
!
!     determining the name of the restart file
!
      ipos=index(jobnamec(1),char(0))
      fnrstrt(1:ipos-1)=jobnamec(1)(1:ipos-1)
      fnrstrt(ipos:ipos+3)=".rin"
      do i=ipos+4,132
        fnrstrt(i:i)=' '
      enddo
!
!     opening the restart file
!
      open(15,file=fnrstrt,ACCESS='SEQUENTIAL',FORM='UNFORMATTED',
     &     err=15)
!
      do
!
        read(15,iostat=istat) version
        if(istat.lt.0) then
          if(irestartstep.eq.0) then
!
!     reading the last step
!
            irestartstep=istep
            close(15)
            open(15,file=fnrstrt,ACCESS='SEQUENTIAL',
     &           FORM='UNFORMATTED',err=15)
            read(15) version
          else
            write(*,*) '*ERROR in restartshort: requested step'
            write(*,*) '       is not in the restart file'
            call exit(201)
          endif
        endif
        read(15)istep
!
!     set size
!
        read(15)nset
        read(15)nalset
!
!     load size
!
        read(15)nload
        read(15)nbody
        read(15)nforc
        read(15)nboun
        read(15)nflow
!
!     mesh size
!
        read(15)nk
        read(15)ne
        read(15)nef
        read(15)nkon
        read(15)(mi(i),i=1,3)
!
!     constraint size
!
        read(15)nmpc
        read(15)mpcend
        read(15)maxlenmpc
        read(15)memmpc_
!
!     coupling, distributing
!
        read(15)nfc
        read(15)ndc
!
!     material size
!
        read(15)nmat
        read(15)ntmat_
        read(15)npmat_
        read(15)ncmat_
!
!     property info
!
        read(15)nprop
!
!     transformation size
!
        read(15)norien
        read(15)ntrans
!
!     amplitude size
!
        read(15)nam
        read(15)namtot
!
!     print size
!
        read(15)nprint
        read(15)nlabel
!
!     tie size
!
        read(15)ntie
!
!     cyclic symmetry size
!
        read(15)ncs_
        read(15)mcs
!
!     1d and 2d element size
!
        read(15)ne1d
        read(15)ne2d
        read(15)(infree(i),i=1,4)
!
!     procedure info
!
        read(15)nmethod
        read(15)(iperturb(i),i=1,2)
        read(15)nener
        read(15)iplas
        read(15)ithermal(1)
        read(15)nstate_
        read(15)nslavs
        read(15)iprestr
        read(15)mortar
        if(mortar.eq.1) then
          read(15)ifacecount
          read(15)nintpoint
        endif
        read(15)nheading_
        read(15)network
!
        if(istep.eq.irestartstep) exit
!
!     skipping the next entries
!
        call skip(nset,nalset,nload,nbody,nforc,nboun,nk,ne,
     &       nkon,mi,nmpc,mpcend,nmat,ntmat_,npmat_,ncmat_,norien,
     &       ntrans,nam,nprint,nlabel,ncs_,ne1d,ne2d,infree,nmethod,
     &       iperturb,nener,ithermal,nstate_,iprestr,mcs,ntie,
     &       nslavs,nprop,mortar,ifacecount,nintpoint,nef,nheading_,
     &       nfc,ndc)
!
      enddo
!
!     sets
!
      read(15)(set(i),i=1,nset)
!
!     the contents of istartset is temporarily stored in meminset
!
      read(15)(meminset(i),i=1,nset)
!
!     the contents of iendset is temporarily stored in rmeminset
!
      read(15)(rmeminset(i),i=1,nset)
!
!     reordering the information of istartset, iendset and ialset
!     into meminset and rmeminset
!
      iact=0
      do j=1,nalset
        if(iact.eq.0) then
          do k=1,nset
            if(meminset(k).eq.j) then
              meminset(k)=0
              mem=0
              iact=1
              exit
            endif
          enddo
          if(k.gt.nset) cycle
        endif
        mem=mem+1
        im2=im1
        im1=im0
        read(15) im0
        if(im0.gt.0) then
          meminset(k)=meminset(k)+1
        else
!
!     im0<0 and two elements are already stored
!
          meminset(k)=meminset(k)+(im2-im1)/im0-1
        endif
        if(rmeminset(k).eq.j) then
          iact=0
          rmeminset(k)=mem
!
!     make set k ineligible in further iterations
!
          meminset(k)=-meminset(k)
        endif
      enddo
!
!     restore the sign of meminset
!
      do k=1,nset
        meminset(k)=-meminset(k)
      enddo
!
      close(15)
!
      return
!
 15   write(*,*) '*ERROR in restartshort: could not open file ',fnrstrt
      call exit(201)
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
C
C************************** SURF_FLOWTRA  ******************************
C
C surface water routing and solute transport procedure:
C
C FOR FLOW: takes OVFLNOD as input from FLOW3D and returns `routed'
C           PONDNOD as BC for FLOW3D
C
C FOR TRANSPORT: takes TRAFLNOD as input from subsurface transport
C                and returns routed CONCNOD that is used to define the
C                Cauchy BC for subsurface transport
C***********************************************************************
C
      subroutine surf_flowtra(ncell,nnod,nrow,ncol,ntri,dostep,numres,
     1                       ncell_coarse,ncelnl,cell,
     2                       indcel,indcelwl,tipo_r,reservr,dem_map,
     3                       lakes_map,cellcol,cellrow,
     4                       celtype,cells_r,tp2d,triang,
     5                       time,deltat,arenod,OVFLNOD,PONDNOD,
     6                       ovflcel,pondcel,TRAFLNOD,TRAFLCEL,CONCNOD,
     7                       CONCCEL,N_HA,NSURF,NSURFT,NSURFT_TB,TRANSP,
     8                       SOURCE_MIXING,MIX_CORRECT,SURFACE_MIX_SN)
C
      implicit none
      include 'CATHY.H'
      INCLUDE 'SURFWATER.H'
      INCLUDE 'RIVERNETWORK.H'
      INCLUDE 'TRANSPSURF.H'
      LOGICAL  TRANSP
      integer  ncell,nnod,nrow,ncol,ntri,dostep,numres
      integer  ncell_coarse,ncelnl
      INTEGER  NSURF,NSURFT,NSURFT_TB
      integer  i,j,k
      integer  cell(5,maxcel),indcel(rowmax,colmax)
      integer  indcelwl(rowmax,colmax)
      integer  tipo_r(*),reservr(*)
      integer  lakes_map(rowmax,*)
      integer  cellcol(*),cellrow(*)
      integer  celtype(*),cells_r(*)
      integer  tp2d(*),triang(4,*)
      integer  n_hA(*)
      real*8   time,deltat
      real*8   dtsurf
      real*8   arenod(*)
      real*8   OVFLNOD(*),PONDNOD(*)
      real*8   TRAFLNOD(*),CONCNOD(*),SOURCE_MIXING(*)
      real*8   ovflcel(*),pondcel(*),dem_map(rowmax,*)
      real*8   traflcel(*),conccel(*),MIX_CORRECT(*)
      real*8   traflcel_mix(maxcel),surface_mix_sn(maxcel)
C
      call init0r(maxcel,traflcel_mix)
      call init0r(maxcel,surface_mix_sn)
c
c  from node to cell
c
cm    call nod_cell(ncell,nrow,ncol,dostep,ncell_coarse,
cm   1              nnod,cell,dem_map,indcelwl,cellcoarse,OVFLNOD,
cm   2              ovflcel,arenod,delta_x,delta_y,pondnod)
      call nod_cell(ncell,nrow,ncol,dostep,ncell_coarse,
     1              nnod,cell,dem_map,indcelwl,cellcoarse,OVFLNOD,
     2              ovflcel,arenod,delta_x,delta_y)
c
c   trasferimento informazioni al surf_route (cambia la numerazione
c   delle celle se ci sono laghi!)
c
      call transfer_f3d_surf(nrow,ncol,indcel,indcelwl,
     1                       surface_water_sn,ovflcel,
     2                       reservr,lakes_map)
c
      IF (TRANSP) THEN
         call nod_cell(ncell,nrow,ncol,dostep,ncell_coarse,
     1        nnod,cell,dem_map,indcelwl,cellcoarse,traflnod,
     2        traflcel,arenod,delta_x,delta_y)

         call nod_cell(ncell,nrow,ncol,dostep,ncell_coarse,
     1        nnod,cell,dem_map,indcelwl,cellcoarse,source_mixing,
     2        traflcel_mix,arenod,delta_x,delta_y)

         call transfer_f3d_surf(nrow,ncol,indcel,indcelwl,
     1        surface_conc_sn,traflcel,
     2        reservr,lakes_map)

         call transfer_f3d_surf(nrow,ncol,indcel,indcelwl,
     1        surface_mix_sn,traflcel_mix,
     2        reservr,lakes_map)
      END IF
c
c  surface routing
c
      cu_max=ak_max*deltat
      CUTRGT = 1.0d0
      if(cu_max .gt. CUTRGT) then
         dtsurf = CUTRGT/ak_max
         nsurf = IDINT((deltat/dtsurf))+1
         dtsurf = deltat/nsurf
      else
         dtsurf = deltat
         nsurf=1
      end if
      NSURFT=NSURFT + NSURF
      NSURFT_TB=NSURFT_TB + NSURF
C
      DO I=1,NSURF
C------------------------WATER ROUTING -----------------------------
         CALL ROUTE(NROW,NCELL,TIPO_R,RESERVR,
     1        CELLS_R,N_HA,TIME,DTSURF)
C     CALCULATION OF SURFACE WATER HEIGHTS
         CALL ALTEZZE(NCELNL,CELTYPE,DELTA_X,
     1        DELTA_Y,Q_IN_KK_SN,Q_IN_KKP1_SN,
     2        Q_OUT_KK_SN_1,Q_OUT_KK_SN_2,
     3        Q_OUT_KKP1_SN_1,Q_OUT_KKP1_SN_2,
     4        VOLUME_KK_SN,VOLUME_KKP1_SN,SURFACE_WATER_SN,
     5        H_WATER_KKP1_SN,DTSURF,
     6        H_POOL_KKP1_VEC,RESERVR,ELEVATION)

C
C---------------------------SURFACE TRANSPORT -------------------------
         IF (TRANSP) THEN
C     CONCENTRATION ROUTING
            CALL ROUTE_TRA(NROW,NCELL,TIPO_R,RESERVR,
     1           CELLS_R,N_HA,TIME,DTSURF)
c
C     CALCULATION OF CONCENTRATION FOR CELLS
            CALL MASBIL(NCELL,TIPO_R,DELTA_X,DELTA_Y,
     1           QMASS_IN_KK_SN,QMASS_IN_KKP1_SN,
     2           QMASS_OUT_KK_SN_1,QMASS_OUT_KK_SN_2,
     3           QMASS_OUT_KKP1_SN_1,QMASS_OUT_KKP1_SN_2,
     4           MASS_KK_SN,MASS_KKP1_SN,SURFACE_CONC_SN,
     5           CONC_KKP1_SN,VOLUME_KKP1_SN,DTSURF,MIX_CORRECT,
     6           surface_mix_sn)
c
         END IF
C-----------------------------------------------------------------------

C INTERNAL UPDATES OF SURFACE ROUTING VARIABLES (I.E., FLOW RATE, WATER HEIGHTS)
            IF(NSURF.GT.1 .AND. I.LT.NSURF) THEN
               CALL VCOPYR(MAXCEL,Q_IN_KK_SN,Q_IN_KKP1_SN)
               CALL INIT0R(MAXCEL,Q_IN_KKP1_SN)
               CALL VCOPYR(MAXCEL,Q_OUT_KK_SN_1,Q_OUT_KKP1_SN_1)
               CALL INIT0R(MAXCEL,Q_OUT_KKP1_SN_1)
               CALL VCOPYR(MAXCEL,Q_OUT_KK_SN_2,Q_OUT_KKP1_SN_2)
               CALL INIT0R(MAXCEL,Q_OUT_KKP1_SN_2)
               CALL VCOPYR(MAXCEL,VOLUME_KK_SN,VOLUME_KKP1_SN)
               CALL INIT0R(MAXCEL,VOLUME_KKP1_SN)
C
               IF (TRANSP) THEN
                  CALL VCOPYR(MAXCEL,QMASS_IN_KK_SN,QMASS_IN_KKP1_SN)
                  CALL INIT0R(MAXCEL,QMASS_IN_KKP1_SN)
                  CALL VCOPYR(MAXCEL,QMASS_OUT_KK_SN_1,
     1                 QMASS_OUT_KKP1_SN_1)
                  CALL INIT0R(MAXCEL,QMASS_OUT_KKP1_SN_1)
                  CALL VCOPYR(MAXCEL,QMASS_OUT_KK_SN_2,
     1                 QMASS_OUT_KKP1_SN_2)
                  CALL INIT0R(MAXCEL,QMASS_OUT_KKP1_SN_2)
                  CALL VCOPYR(MAXCEL,MASS_KK_SN,MASS_KKP1_SN)
                  CALL INIT0R(MAXCEL,MASS_KKP1_SN)
               END IF
            END IF
         END DO

c
c trasferimento informazioni per il flow3d
c
      call transfer_surf_f3d(nrow,ncol,indcel,indcelwl,
     1                       ncell,h_water_kkp1_sn,
     2                       pondcel,h_pool_kkp1_vec,
     3                       lakes_map,elevation_with_lakes)
c
c  calcola PONDNOD da PONDCEL passando per i triangoli
c
cm    call cell_nod(ncell,nrow,ncol,nnod,ntri,dostep,
cm   1                    cell,tp2d,triang,indcel,dem_map,
cm   2                    PONDCEL,cellcoarse,pondnod,
cm   3                    delta_x,delta_y,arenod)
      call cell_nod(ncell,nrow,ncol,nnod,ntri,dostep,
     1              tp2d,triang,indcelwl,dem_map,
     2              pondcel,cellcoarse,PONDNOD)
c
      IF (TRANSP) THEN
         call transfer_surf_f3d(nrow,ncol,indcel,indcelwl,
     1        ncell,CONC_KKP1_SN,
     2        CONCCEL,h_pool_kkp1_vec,
     3        lakes_map,elevation_with_lakes)
c
         call cell_nod_tra(ntri,nnod,tp2d,triang,CONCCEL,concnod,
     3                    delta_x,delta_y,arenod,pondcel,pondnod)
cm       call cell_nod(ncell,nrow,ncol,nnod,ntri,dostep,
cm   1        tp2d,triang,indcelwl,dem_map,
cm   2        CONCCEL,cellcoarse,CONCNOD)
c
      END IF
      return
      end