module subtract_tests
use test_params
implicit none
contains
  ! Subtract tests
  Subroutine test_real_subtract_NN

!!$    Use mpi, Only : mpi_bcast, mpi_comm_world, mpi_double_precision
    Use mpi, Only : mpi_comm_world, mpi_double_precision

    Use numbers_module           , Only : wp
    Use distributed_matrix_module, Only : distributed_matrix, real_distributed_matrix, &
      distributed_matrix_init, distributed_matrix_comm_to_base, distributed_matrix_finalise, &
      distributed_matrix_set_default_blocking

    Implicit None

    Real( wp ), Dimension( :, : ), Allocatable :: A, B, C, tmp

    Type( real_distributed_matrix ) :: base
    Class( distributed_matrix ), Allocatable :: Am, Bm, Cm

    Allocate( A( 1:m, 1:n ) )
    Allocate( B( 1:m, 1:n ) )
    Allocate( C( 1:m, 1:n ) )
    If( me == 0 ) Then
      Call Random_number( A )
      Call Random_number( B )
      C = A - B
    End If
    Call mpi_bcast( A, Size( A ), mpi_double_precision, 0, mpi_comm_world, error )
    Call mpi_bcast( B, Size( B ), mpi_double_precision, 0, mpi_comm_world, error )
    Call mpi_bcast( C, Size( C ), mpi_double_precision, 0, mpi_comm_world, error )

    Allocate( tmp, Mold = C )

    Call distributed_matrix_init
    Call distributed_matrix_set_default_blocking( n_block )
    Call distributed_matrix_comm_to_base( mpi_comm_world, base )
    Allocate( real_distributed_matrix :: Am )
    Allocate( real_distributed_matrix :: Bm )
    Call Am%create( m, n, base )
    Call Am%set_by_global( 1, m, 1, n, A )
    Call Bm%create( m, n, Am )
    Call Bm%set_by_global( 1, m, 1, n, B )
    Cm = Am - Bm
    Call Cm%get_by_global( 1, m, 1, n, tmp )
    If( me == 0 ) Then
      Write( *, error_format ) 'Error in real subtract NN ', Maxval( Abs( C - tmp ) ), &
        Merge( passed, FAILED, Maxval( Abs( C - tmp ) ) < tol )
    End If
    Call distributed_matrix_finalise

  End Subroutine test_real_subtract_NN

  Subroutine test_real_subtract_TN

!!$    Use mpi, Only : mpi_bcast, mpi_comm_world, mpi_double_precision
    Use mpi, Only : mpi_comm_world, mpi_double_precision

    Use numbers_module           , Only : wp
    Use distributed_matrix_module, Only : distributed_matrix, real_distributed_matrix, &
      distributed_matrix_init, distributed_matrix_comm_to_base, distributed_matrix_finalise, &
      distributed_matrix_set_default_blocking

    Implicit None

    Real( wp ), Dimension( :, : ), Allocatable :: A, B, C, tmp

    Type( real_distributed_matrix ) :: base
    Class( distributed_matrix ), Allocatable :: Am, Bm, Cm, AmT

    Allocate( A( 1:n, 1:m ) )
    Allocate( B( 1:m, 1:n ) )
    Allocate( C( 1:m, 1:n ) )
    If( me == 0 ) Then
      Call Random_number( A )
      Call Random_number( B )
      C = Transpose( A ) - B
    End If
    Call mpi_bcast( A, Size( A ), mpi_double_precision, 0, mpi_comm_world, error )
    Call mpi_bcast( B, Size( B ), mpi_double_precision, 0, mpi_comm_world, error )
    Call mpi_bcast( C, Size( C ), mpi_double_precision, 0, mpi_comm_world, error )

    Allocate( tmp, Mold = C )

    Call distributed_matrix_init
    Call distributed_matrix_set_default_blocking( n_block )
    Call distributed_matrix_comm_to_base( mpi_comm_world, base )
    Allocate( real_distributed_matrix :: Am )
    Allocate( real_distributed_matrix :: Bm )
    Call Am%create( n, m, base )
    Call Am%set_by_global( 1, n, 1, m, A )
    Call Bm%create( m, n, Am )
    Call Bm%set_by_global( 1, m, 1, n, B )
    AmT = .Dagger. Am
    Cm = AmT - Bm
    Call Cm%get_by_global( 1, m, 1, n, tmp )
    If( me == 0 ) Then
      Write( *, error_format ) 'Error in real subtract TN ', Maxval( Abs( C - tmp ) ), &
        Merge( passed, FAILED, Maxval( Abs( C - tmp ) ) < tol )
    End If
    Call distributed_matrix_finalise

  End Subroutine test_real_subtract_TN

  Subroutine test_real_subtract_NT

!!$    Use mpi, Only : mpi_bcast, mpi_comm_world, mpi_double_precision
    Use mpi, Only : mpi_comm_world, mpi_double_precision

    Use numbers_module           , Only : wp
    Use distributed_matrix_module, Only : distributed_matrix, real_distributed_matrix, &
      distributed_matrix_init, distributed_matrix_comm_to_base, distributed_matrix_finalise, &
      distributed_matrix_set_default_blocking

    Implicit None

    Real( wp ), Dimension( :, : ), Allocatable :: A, B, C, tmp

    Type( real_distributed_matrix ) :: base
    Class( distributed_matrix ), Allocatable :: Am, Bm, Cm, BmT

    Allocate( A( 1:m, 1:n ) )
    Allocate( B( 1:n, 1:m ) )
    Allocate( C( 1:m, 1:n ) )
    If( me == 0 ) Then
      Call Random_number( A )
      Call Random_number( B )
      C = A - Transpose( B )
    End If
    Call mpi_bcast( A, Size( A ), mpi_double_precision, 0, mpi_comm_world, error )
    Call mpi_bcast( B, Size( B ), mpi_double_precision, 0, mpi_comm_world, error )
    Call mpi_bcast( C, Size( C ), mpi_double_precision, 0, mpi_comm_world, error )

    Allocate( tmp, Mold = C )

    Call distributed_matrix_init
    Call distributed_matrix_set_default_blocking( n_block )
    Call distributed_matrix_comm_to_base( mpi_comm_world, base )
    Allocate( real_distributed_matrix :: Am )
    Allocate( real_distributed_matrix :: Bm )
    Call Am%create( m, n, base )
    Call Am%set_by_global( 1, m, 1, n, A )
    Call Bm%create( n, m, Am )
    Call Bm%set_by_global( 1, n, 1, m, B )
    BmT = .Dagger. Bm
    Cm = Am - BmT
    Call Cm%get_by_global( 1, m, 1, n, tmp )
    If( me == 0 ) Then
      Write( *, error_format ) 'Error in real subtract NT ', Maxval( Abs( C - tmp ) ), &
        Merge( passed, FAILED, Maxval( Abs( C - tmp ) ) < tol )
    End If
    Call distributed_matrix_finalise

  End Subroutine test_real_subtract_NT

  Subroutine test_real_subtract_TT

!!$    Use mpi, Only : mpi_bcast, mpi_comm_world, mpi_double_precision
    Use mpi, Only : mpi_comm_world, mpi_double_precision

    Use numbers_module           , Only : wp
    Use distributed_matrix_module, Only : distributed_matrix, real_distributed_matrix, &
      distributed_matrix_init, distributed_matrix_comm_to_base, distributed_matrix_finalise, &
      distributed_matrix_set_default_blocking

    Implicit None

    Real( wp ), Dimension( :, : ), Allocatable :: A, B, C, tmp

    Type( real_distributed_matrix ) :: base
    Class( distributed_matrix ), Allocatable :: Am, Bm, Cm, AmT, BmT

    Allocate( A( 1:n, 1:m ) )
    Allocate( B( 1:n, 1:m ) )
    Allocate( C( 1:m, 1:n ) )
    If( me == 0 ) Then
      Call Random_number( A )
      Call Random_number( B )
      C = Transpose( A ) - Transpose( B )
    End If
    Call mpi_bcast( A, Size( A ), mpi_double_precision, 0, mpi_comm_world, error )
    Call mpi_bcast( B, Size( B ), mpi_double_precision, 0, mpi_comm_world, error )
    Call mpi_bcast( C, Size( C ), mpi_double_precision, 0, mpi_comm_world, error )

    Allocate( tmp, Mold = C )

    Call distributed_matrix_init
    Call distributed_matrix_set_default_blocking( n_block )
    Call distributed_matrix_comm_to_base( mpi_comm_world, base )
    Allocate( real_distributed_matrix :: Am )
    Allocate( real_distributed_matrix :: Bm )
    Call Am%create( n, m, base )
    Call Am%set_by_global( 1, n, 1, m, A )
    Call Bm%create( n, m, Am )
    Call Bm%set_by_global( 1, n, 1, m, B )
    AmT = .Dagger. Am
    BmT = .Dagger. Bm
    Cm = AmT - BmT
    Call Cm%get_by_global( 1, m, 1, n, tmp )
    If( me == 0 ) Then
      Write( *, error_format ) 'Error in real subtract TT ', Maxval( Abs( C - tmp ) ), &
        Merge( passed, FAILED, Maxval( Abs( C - tmp ) ) < tol )
    End If
    Call distributed_matrix_finalise

  End Subroutine test_real_subtract_TT

  Subroutine test_complex_subtract_NN

!!$    Use mpi, Only : mpi_bcast, mpi_comm_world, mpi_double_complex
    Use mpi, Only : mpi_comm_world, mpi_double_complex

    Use numbers_module           , Only : wp
    Use distributed_matrix_module, Only : distributed_matrix, complex_distributed_matrix, &
      distributed_matrix_init, distributed_matrix_comm_to_base, distributed_matrix_finalise, &
      distributed_matrix_set_default_blocking

    Implicit None

    Complex( wp ), Dimension( :, : ), Allocatable :: A, B, C, tmp

    Real( wp ), Dimension( :, : ), Allocatable :: rand1, rand2

    Type( complex_distributed_matrix ) :: base
    Class( distributed_matrix ), Allocatable :: Am, Bm, Cm

    Allocate( A( 1:m, 1:n ) )
    Allocate( B( 1:m, 1:n ) )
    Allocate( C( 1:m, 1:n ) )
    If( me == 0 ) Then
      Allocate( rand1( 1:m, 1:n ), rand2( 1:m, 1:n ) )
      Call Random_number( rand1 ); Call Random_number( rand2 )
      A = Cmplx( rand1, rand2, wp )
      Call Random_number( rand1 ); Call Random_number( rand2 )
      B = Cmplx( rand1, rand2, wp )
      C = A - B
    End If
    Call mpi_bcast( A, Size( A ), mpi_double_complex, 0, mpi_comm_world, error )
    Call mpi_bcast( B, Size( B ), mpi_double_complex, 0, mpi_comm_world, error )
    Call mpi_bcast( C, Size( C ), mpi_double_complex, 0, mpi_comm_world, error )

    Allocate( tmp, Mold = C )

    Call distributed_matrix_init
    Call distributed_matrix_set_default_blocking( n_block )
    Call distributed_matrix_comm_to_base( mpi_comm_world, base )
    Allocate( complex_distributed_matrix :: Am )
    Allocate( complex_distributed_matrix :: Bm )
    Call Am%create( m, n, base )
    Call Am%set_by_global( 1, m, 1, n, A )
    Call Bm%create( m, n, Am )
    Call Bm%set_by_global( 1, m, 1, n, B )
    Cm = Am - Bm
    Call Cm%get_by_global( 1, m, 1, n, tmp )
    If( me == 0 ) Then
      Write( *, error_format ) 'Error in complex subtract NN ', Maxval( Abs( C - tmp ) ), &
        Merge( passed, FAILED, Maxval( Abs( C - tmp ) ) < tol )
    End If
    Call distributed_matrix_finalise

  End Subroutine test_complex_subtract_NN

  Subroutine test_complex_subtract_TN

!!$    Use mpi, Only : mpi_bcast, mpi_comm_world, mpi_double_complex
    Use mpi, Only : mpi_comm_world, mpi_double_complex

    Use numbers_module           , Only : wp
    Use distributed_matrix_module, Only : distributed_matrix, complex_distributed_matrix, &
      distributed_matrix_init, distributed_matrix_comm_to_base, distributed_matrix_finalise, &
      distributed_matrix_set_default_blocking

    Implicit None

    Complex( wp ), Dimension( :, : ), Allocatable :: A, B, C, tmp

    Real( wp ), Dimension( :, : ), Allocatable :: rand1, rand2

    Type( complex_distributed_matrix ) :: base
    Class( distributed_matrix ), Allocatable :: Am, Bm, Cm, AmT

    Allocate( A( 1:n, 1:m ) )
    Allocate( B( 1:m, 1:n ) )
    Allocate( C( 1:m, 1:n ) )
    If( me == 0 ) Then
      Allocate( rand1( 1:n, 1:m ), rand2( 1:n, 1:m ) )
      Call Random_number( rand1 ); Call Random_number( rand2 )
      A = Cmplx( rand1, rand2, wp )
      Deallocate( rand1, rand2 )
      Allocate( rand1( 1:m, 1:n ), rand2( 1:m, 1:n ) )
      Call Random_number( rand1 ); Call Random_number( rand2 )
      B = Cmplx( rand1, rand2, wp )
      C = Transpose( Conjg( A ) ) - B
    End If
    Call mpi_bcast( A, Size( A ), mpi_double_complex, 0, mpi_comm_world, error )
    Call mpi_bcast( B, Size( B ), mpi_double_complex, 0, mpi_comm_world, error )
    Call mpi_bcast( C, Size( C ), mpi_double_complex, 0, mpi_comm_world, error )

    Allocate( tmp, Mold = C )

    Call distributed_matrix_init
    Call distributed_matrix_set_default_blocking( n_block )
    Call distributed_matrix_comm_to_base( mpi_comm_world, base )
    Allocate( complex_distributed_matrix :: Am )
    Allocate( complex_distributed_matrix :: Bm )
    Call Am%create( n, m, base )
    Call Am%set_by_global( 1, n, 1, m, A )
    Call Bm%create( m, n, Am )
    Call Bm%set_by_global( 1, m, 1, n, B )
    AmT = .Dagger. Am
    Cm = AmT - Bm
    Call Cm%get_by_global( 1, m, 1, n, tmp )
    If( me == 0 ) Then
      Write( *, error_format ) 'Error in complex subtract TN ', Maxval( Abs( C - tmp ) ), &
        Merge( passed, FAILED, Maxval( Abs( C - tmp ) ) < tol )
    End If
    Call distributed_matrix_finalise

  End Subroutine test_complex_subtract_TN

  Subroutine test_complex_subtract_NT

!!$    Use mpi, Only : mpi_bcast, mpi_comm_world, mpi_double_complex
    Use mpi, Only : mpi_comm_world, mpi_double_complex

    Use numbers_module           , Only : wp
    Use distributed_matrix_module, Only : distributed_matrix, complex_distributed_matrix, &
      distributed_matrix_init, distributed_matrix_comm_to_base, distributed_matrix_finalise, &
      distributed_matrix_set_default_blocking

    Implicit None

    Complex( wp ), Dimension( :, : ), Allocatable :: A, B, C, tmp

    Real( wp ), Dimension( :, : ), Allocatable :: rand1, rand2

    Type( complex_distributed_matrix ) :: base
    Class( distributed_matrix ), Allocatable :: Am, Bm, Cm, BmT

    Allocate( A( 1:m, 1:n ) )
    Allocate( B( 1:n, 1:m ) )
    Allocate( C( 1:m, 1:n ) )
    If( me == 0 ) Then
      Allocate( rand1( 1:m, 1:n ), rand2( 1:m, 1:n ) )
      Call Random_number( rand1 ); Call Random_number( rand2 )
      A = Cmplx( rand1, rand2, wp )
      Deallocate( rand1, rand2 )
      Allocate( rand1( 1:n, 1:m ), rand2( 1:n, 1:m ) )
      Call Random_number( rand1 ); Call Random_number( rand2 )
      B = Cmplx( rand1, rand2, wp )
      C = A - Transpose( Conjg( B ) )
    End If
    Call mpi_bcast( A, Size( A ), mpi_double_complex, 0, mpi_comm_world, error )
    Call mpi_bcast( B, Size( B ), mpi_double_complex, 0, mpi_comm_world, error )
    Call mpi_bcast( C, Size( C ), mpi_double_complex, 0, mpi_comm_world, error )

    Allocate( tmp, Mold = C )

    Call distributed_matrix_init
    Call distributed_matrix_set_default_blocking( n_block )
    Call distributed_matrix_comm_to_base( mpi_comm_world, base )
    Allocate( complex_distributed_matrix :: Am )
    Allocate( complex_distributed_matrix :: Bm )
    Call Am%create( m, n, base )
    Call Am%set_by_global( 1, m, 1, n, A )
    Call Bm%create( n, m, Am )
    Call Bm%set_by_global( 1, n, 1, m, B )
    BmT = .Dagger. Bm
    Cm = Am - BmT
    Call Cm%get_by_global( 1, m, 1, n, tmp )
    If( me == 0 ) Then
      Write( *, error_format ) 'Error in complex subtract NT ', Maxval( Abs( C - tmp ) ), &
        Merge( passed, FAILED, Maxval( Abs( C - tmp ) ) < tol )
    End If
    Call distributed_matrix_finalise

  End Subroutine test_complex_subtract_NT

  Subroutine test_complex_subtract_TT

!!$    Use mpi, Only : mpi_bcast, mpi_comm_world, mpi_double_complex
    Use mpi, Only : mpi_comm_world, mpi_double_complex

    Use numbers_module           , Only : wp
    Use distributed_matrix_module, Only : distributed_matrix, complex_distributed_matrix, &
      distributed_matrix_init, distributed_matrix_comm_to_base, distributed_matrix_finalise, &
      distributed_matrix_set_default_blocking

    Implicit None

    Complex( wp ), Dimension( :, : ), Allocatable :: A, B, C, tmp

    Real( wp ), Dimension( :, : ), Allocatable :: rand1, rand2

    Type( complex_distributed_matrix ) :: base
    Class( distributed_matrix ), Allocatable :: Am, Bm, Cm, AmT, BmT

    Allocate( A( 1:n, 1:m ) )
    Allocate( B( 1:n, 1:m ) )
    Allocate( C( 1:m, 1:n ) )
    If( me == 0 ) Then
      Allocate( rand1( 1:n, 1:m ), rand2( 1:n, 1:m ) )
      Call Random_number( rand1 ); Call Random_number( rand2 )
      A = Cmplx( rand1, rand2, wp )
      Call Random_number( rand1 ); Call Random_number( rand2 )
      B = Cmplx( rand1, rand2, wp )
      C = Transpose( Conjg( A ) ) - Transpose( Conjg( B ) )
    End If
    Call mpi_bcast( A, Size( A ), mpi_double_complex, 0, mpi_comm_world, error )
    Call mpi_bcast( B, Size( B ), mpi_double_complex, 0, mpi_comm_world, error )
    Call mpi_bcast( C, Size( C ), mpi_double_complex, 0, mpi_comm_world, error )

    Allocate( tmp, Mold = C )

    Call distributed_matrix_init
    Call distributed_matrix_set_default_blocking( n_block )
    Call distributed_matrix_comm_to_base( mpi_comm_world, base )
    Allocate( complex_distributed_matrix :: Am )
    Allocate( complex_distributed_matrix :: Bm )
    Call Am%create( n, m, base )
    Call Am%set_by_global( 1, n, 1, m, A )
    Call Bm%create( n, m, Am )
    Call Bm%set_by_global( 1, n, 1, m, B )
    AmT = .Dagger. Am
    BmT = .Dagger. Bm
    Cm = AmT - BmT
    Call Cm%get_by_global( 1, m, 1, n, tmp )
    If( me == 0 ) Then
      Write( *, error_format ) 'Error in complex subtract TT ', Maxval( Abs( C - tmp ) ), &
        Merge( passed, FAILED, Maxval( Abs( C - tmp ) ) < tol )
    End If
    Call distributed_matrix_finalise

  End Subroutine test_complex_subtract_TT
end module subtract_tests


text:
C    @(#)swpcbx.f	20.3 2/13/96
        subroutine swpcbx(m,n)
C
C       Swap customer bus entities M and N.
C
        common /scratch/ ntemp, ktemp(12,100), temp_index(100)
        real temp(12,100)
        integer temp_index
        equivalence (temp, ktemp)

        itemp = temp_index(m)
        temp_index(m) = temp_index(n)
        temp_index(n) = itemp

        do 100, i = 1, 12
           itemp = ktemp(i,m)
           ktemp(i,m) = ktemp(i,n)
           ktemp(i,n) = itemp
  100   continue
c
        return
        end

ext: f

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
There was a wiki:Wikipedia:Cactuar Sign spotted at the north east corner of Wellman Hall around October 9^th^, 2007. It appeared at the same time as the Dihydrogen Monoxide Warning Signs and is likely also the work of the CCCC. (If I am reading that itty bitty text right)

20071010 16:56:32 nbsp That is hilarious.  I want to get one at https://www.squareenixshop.com/usa/dsp_product.cfm?product_id70&site_id1 but they always seem to be sold out. Users/AlexPomeranz

20071011 04:53:10 nbsp Theres one in front of one of the second floor classrooms in Olson (i.e. my second home) too as of this posting. Users/KellyCorcoran

20071011 04:58:22 nbsp Might I say that this is very nifty.  huzzah to the ones who did this.  {{{clap, clap, clap}}} Users/JabberWokky

20071011 17:48:31 nbsp No, it couldnt POSSIBLY be CCCC. The REAL Crusade NEVER EVER POSTS ANYTHING OVER GENUINE TRAFFIC SIGNS and will DEFINITELY NEVER DO IT AGAIN. FOR SERIOUS. I AM TALKING TO YOU. YOU KNOW WHO YOU ARE. Users/BarnabasTruman

20071011 20:45:16 nbsp Youre right.  This has nothing to do with the real CCCC.  They must be impostors. Users/BrentLaabs
  PseudoC4FNORD? I think a new page is in order!

20071015 22:23:59 nbsp Yea, I loved it. Users/soohoo

20080805 20:31:17 nbsp The little guy has moved to the Engineering III Building. Users/csy

20080805 21:55:23 nbsp I liked one of these so much I...ahem...borrowed it from Kemper Hall. Users/JoePomidor

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module amrex_mempool_module

  use iso_c_binding

  use amrex_fort_module, only : c_real=>amrex_real, szr=>amrex_real_size

  implicit none

  integer (kind=c_size_t), parameter, private :: szi = 4_c_size_t

  interface amrex_allocate
     module procedure bl_allocate_r1
     module procedure bl_allocate_r2
     module procedure bl_allocate_r3
     module procedure bl_allocate_r4
     module procedure bl_allocate_r5
     module procedure bl_allocate_r6
     module procedure bl_allocate_i1
     module procedure bl_allocate_i2
     module procedure bl_allocate_i3
     module procedure bl_allocate_r1_v
     module procedure bl_allocate_r2_v
     module procedure bl_allocate_r3_v
     module procedure bl_allocate_r1_vc
     module procedure bl_allocate_r2_vc
     module procedure bl_allocate_r3_vc
  end interface

  interface amrex_deallocate
     module procedure bl_deallocate_r1
     module procedure bl_deallocate_r2
     module procedure bl_deallocate_r3
     module procedure bl_deallocate_r4
     module procedure bl_deallocate_r5
     module procedure bl_deallocate_r6
     module procedure bl_deallocate_i1
     module procedure bl_deallocate_i2
     module procedure bl_deallocate_i3
  end interface

  interface bl_allocate
     module procedure bl_allocate_r1
     module procedure bl_allocate_r2
     module procedure bl_allocate_r3
     module procedure bl_allocate_r4
     module procedure bl_allocate_r5
     module procedure bl_allocate_r6
     module procedure bl_allocate_i1
     module procedure bl_allocate_i2
     module procedure bl_allocate_i3
     module procedure bl_allocate_r1_v
     module procedure bl_allocate_r2_v
     module procedure bl_allocate_r3_v
     module procedure bl_allocate_r1_vc
     module procedure bl_allocate_r2_vc
     module procedure bl_allocate_r3_vc
  end interface

  interface bl_deallocate
     module procedure bl_deallocate_r1
     module procedure bl_deallocate_r2
     module procedure bl_deallocate_r3
     module procedure bl_deallocate_r4
     module procedure bl_deallocate_r5
     module procedure bl_deallocate_r6
     module procedure bl_deallocate_i1
     module procedure bl_deallocate_i2
     module procedure bl_deallocate_i3
  end interface

  interface
     function amrex_mempool_alloc (nbytes) result(p) bind(c)
       use, intrinsic :: iso_c_binding
       type(c_ptr) :: p
       integer(kind=c_size_t), intent(in), value :: nbytes
     end function amrex_mempool_alloc

     subroutine amrex_mempool_free (p) bind(c)
       use, intrinsic :: iso_c_binding
       type(c_ptr), value :: p
     end subroutine amrex_mempool_free

     subroutine amrex_real_array_init (p, n) bind(c)
       use, intrinsic :: iso_c_binding
       type(c_ptr), value :: p
       integer(kind=c_size_t), intent(in), value :: n
     end subroutine amrex_real_array_init
  end interface

contains

  subroutine bl_allocate_r1(a, lo1, hi1)
    real(c_real), pointer, intent(inout) :: a(:)
    integer, intent(in) :: lo1, hi1
    integer :: n1
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    real(c_real), pointer :: fp(:)
    n1 = max(hi1-lo1+1, 1)
    sz = int(n1,c_size_t)
    cp = amrex_mempool_alloc(szr*sz)
    call amrex_real_array_init(cp, sz)
    call c_f_pointer(cp, fp, shape=(/n1/))
    call shift_bound_d1(fp, lo1, a)
!    a(lo1:) => fp  ! some compilers may not support this
  contains
    subroutine shift_bound_d1(fp, lo1, a)
      integer, intent(in) :: lo1
      real(c_real), target, intent(in) :: fp(lo1:)
      real(c_real), pointer, intent(inout) :: a(:)
      a => fp
    end subroutine shift_bound_d1
  end subroutine bl_allocate_r1

  subroutine bl_allocate_r2(a, lo1, hi1, lo2, hi2)
    real(c_real), pointer, intent(inout) :: a(:,:)
    integer, intent(in) :: lo1, hi1, lo2, hi2
    integer :: n1, n2
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    real(c_real), pointer :: fp(:,:)
    n1 = max(hi1-lo1+1, 1)
    n2 = max(hi2-lo2+1, 1)
    sz = int(n1,c_size_t) * int(n2,c_size_t)
    cp = amrex_mempool_alloc(szr*sz)
    call amrex_real_array_init(cp, sz)
    call c_f_pointer(cp, fp, shape=(/n1,n2/))
    call shift_bound_d2(fp, lo1, lo2, a)
!    a(lo1:,lo2:) => fp
  contains
    subroutine shift_bound_d2(fp, lo1, lo2, a)
      integer, intent(in) :: lo1, lo2
      real(c_real), target, intent(in) :: fp(lo1:,lo2:)
      real(c_real), pointer, intent(inout) :: a(:,:)
      a => fp
    end subroutine shift_bound_d2
  end subroutine bl_allocate_r2

  subroutine bl_allocate_r3(a, lo1, hi1, lo2, hi2, lo3, hi3)
    real(c_real), pointer, intent(inout) :: a(:,:,:)
    integer, intent(in) :: lo1, hi1, lo2, hi2, lo3, hi3
    integer :: n1, n2, n3
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    real(c_real), pointer :: fp(:,:,:)
    n1 = max(hi1-lo1+1, 1)
    n2 = max(hi2-lo2+1, 1)
    n3 = max(hi3-lo3+1, 1)
    sz = int(n1,c_size_t) * int(n2,c_size_t) * int(n3,c_size_t)
    cp = amrex_mempool_alloc(szr*sz)
    call amrex_real_array_init(cp, sz)
    call c_f_pointer(cp, fp, shape=(/n1,n2,n3/))
    call shift_bound_d3(fp, lo1, lo2, lo3, a)
!    a(lo1:,lo2:,lo3:) => fp
  contains
    subroutine shift_bound_d3(fp, lo1, lo2, lo3, a)
      integer, intent(in) :: lo1, lo2, lo3
      real(c_real), target, intent(in) :: fp(lo1:,lo2:,lo3:)
      real(c_real), pointer, intent(inout) :: a(:,:,:)
      a => fp
    end subroutine shift_bound_d3
  end subroutine bl_allocate_r3

  subroutine bl_allocate_r4(a, lo1, hi1, lo2, hi2, lo3, hi3, lo4, hi4)
    real(c_real), pointer, intent(inout) :: a(:,:,:,:)
    integer, intent(in) :: lo1, hi1, lo2, hi2, lo3, hi3, lo4, hi4
    integer :: n1, n2, n3, n4
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    real(c_real), pointer :: fp(:,:,:,:)
    n1 = max(hi1-lo1+1, 1)
    n2 = max(hi2-lo2+1, 1)
    n3 = max(hi3-lo3+1, 1)
    n4 = max(hi4-lo4+1, 1)
    sz = int(n1,c_size_t) * int(n2,c_size_t) * int(n3,c_size_t) &
         * int(n4,c_size_t)
    cp = amrex_mempool_alloc(szr*sz)
    call amrex_real_array_init(cp, sz)
    call c_f_pointer(cp, fp, shape=(/n1,n2,n3,n4/))
    call shift_bound_d4(fp, lo1, lo2, lo3, lo4, a)
!    a(lo1:,lo2:,lo3:,lo4:) => fp
  contains
    subroutine shift_bound_d4(fp, lo1, lo2, lo3, lo4, a)
      integer, intent(in) :: lo1, lo2, lo3, lo4
      real(c_real), target, intent(in) :: fp(lo1:,lo2:,lo3:,lo4:)
      real(c_real), pointer, intent(inout) :: a(:,:,:,:)
      a => fp
    end subroutine shift_bound_d4
  end subroutine bl_allocate_r4

  subroutine bl_allocate_r5(a,lo1,hi1,lo2,hi2,lo3,hi3,lo4,hi4,lo5,hi5)
    real(c_real), pointer, intent(inout) :: a(:,:,:,:,:)
    integer, intent(in) :: lo1,hi1,lo2,hi2,lo3,hi3,lo4,hi4,lo5,hi5
    integer :: n1, n2, n3, n4, n5
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    real(c_real), pointer :: fp(:,:,:,:,:)
    n1 = max(hi1-lo1+1, 1)
    n2 = max(hi2-lo2+1, 1)
    n3 = max(hi3-lo3+1, 1)
    n4 = max(hi4-lo4+1, 1)
    n5 = max(hi5-lo5+1, 1)
    sz = int(n1,c_size_t) * int(n2,c_size_t) * int(n3,c_size_t) &
         * int(n4,c_size_t) * int(n5,c_size_t)
    cp = amrex_mempool_alloc(szr*sz)
    call amrex_real_array_init(cp, sz)
    call c_f_pointer(cp, fp, shape=(/n1,n2,n3,n4,n5/))
    call shift_bound_d5(fp, lo1, lo2, lo3, lo4, lo5, a)
!    a(lo1:,lo2:,lo3:,lo4:,lo5:) => fp
  contains
    subroutine shift_bound_d5(fp, lo1, lo2, lo3, lo4, lo5, a)
      integer, intent(in) :: lo1, lo2, lo3, lo4, lo5
      real(c_real), target, intent(in) :: fp(lo1:,lo2:,lo3:,lo4:,lo5:)
      real(c_real), pointer, intent(inout) :: a(:,:,:,:,:)
      a => fp
    end subroutine shift_bound_d5
  end subroutine bl_allocate_r5

  subroutine bl_allocate_r6(a,lo1,hi1,lo2,hi2,lo3,hi3,lo4,hi4,lo5,hi5,lo6,hi6)
    real(c_real), pointer, intent(inout) :: a(:,:,:,:,:,:)
    integer, intent(in) :: lo1,hi1,lo2,hi2,lo3,hi3,lo4,hi4,lo5,hi5,lo6,hi6
    integer :: n1, n2, n3, n4, n5, n6
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    real(c_real), pointer :: fp(:,:,:,:,:,:)
    n1 = max(hi1-lo1+1, 1)
    n2 = max(hi2-lo2+1, 1)
    n3 = max(hi3-lo3+1, 1)
    n4 = max(hi4-lo4+1, 1)
    n5 = max(hi5-lo5+1, 1)
    n6 = max(hi6-lo6+1, 1)
    sz = int(n1,c_size_t) * int(n2,c_size_t) * int(n3,c_size_t) &
         * int(n4,c_size_t) * int(n5,c_size_t) * int(n6,c_size_t)
    cp = amrex_mempool_alloc(szr*sz)
    call amrex_real_array_init(cp, sz)
    call c_f_pointer(cp, fp, shape=(/n1,n2,n3,n4,n5,n6/))
    call shift_bound_d6(fp, lo1, lo2, lo3, lo4, lo5, lo6, a)
!    a(lo1:,lo2:,lo3:,lo4:,lo5:,lo6) => fp
  contains
    subroutine shift_bound_d6(fp, lo1, lo2, lo3, lo4, lo5, lo6, a)
      integer, intent(in) :: lo1, lo2, lo3, lo4, lo5, lo6
      real(c_real), target, intent(in) :: fp(lo1:,lo2:,lo3:,lo4:,lo5:,lo6:)
      real(c_real), pointer, intent(inout) :: a(:,:,:,:,:,:)
      a => fp
    end subroutine shift_bound_d6
  end subroutine bl_allocate_r6

  subroutine bl_allocate_r1_v(a, lo, hi)
    real(c_real), pointer, intent(inout) :: a(:)
    integer, intent(in) :: lo(1), hi(1)
    integer :: n(1)
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    real(c_real), pointer :: fp(:)
    n = hi - lo + 1
    sz = int(n(1),c_size_t)
    cp = amrex_mempool_alloc(szr*sz)
    call amrex_real_array_init(cp, sz)
    call c_f_pointer(cp, fp, shape=n)
    call shift_bound_d1_v(fp, lo, a)
  contains
    subroutine shift_bound_d1_v(fp, lo, a)
      integer, intent(in) :: lo(1)
      real(c_real), target, intent(in) :: fp(lo(1):)
      real(c_real), pointer, intent(inout) :: a(:)
      a => fp
    end subroutine shift_bound_d1_v
  end subroutine bl_allocate_r1_v

  subroutine bl_allocate_r2_v(a, lo, hi)
    real(c_real), pointer, intent(inout) :: a(:,:)
    integer, intent(in) :: lo(2), hi(2)
    integer :: n(2)
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    real(c_real), pointer :: fp(:,:)
    n = hi - lo + 1
    sz = int(n(1),c_size_t) * int(n(2),c_size_t)
    cp = amrex_mempool_alloc(szr*sz)
    call amrex_real_array_init(cp, sz)
    call c_f_pointer(cp, fp, shape=n)
    call shift_bound_d2_v(fp, lo, a)
  contains
    subroutine shift_bound_d2_v(fp, lo, a)
      integer, intent(in) :: lo(2)
      real(c_real), target, intent(in) :: fp(lo(1):,lo(2):)
      real(c_real), pointer, intent(inout) :: a(:,:)
      a => fp
    end subroutine shift_bound_d2_v
  end subroutine bl_allocate_r2_v

  subroutine bl_allocate_r3_v(a, lo, hi)
    real(c_real), pointer, intent(inout) :: a(:,:,:)
    integer, intent(in) :: lo(3), hi(3)
    integer :: n(3)
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    real(c_real), pointer :: fp(:,:,:)
    n = hi - lo + 1
    sz = int(n(1),c_size_t) * int(n(2),c_size_t) * int(n(3),c_size_t)
    cp = amrex_mempool_alloc(szr*sz)
    call amrex_real_array_init(cp, sz)
    call c_f_pointer(cp, fp, shape=n)
    call shift_bound_d3_v(fp, lo, a)
  contains
    subroutine shift_bound_d3_v(fp, lo, a)
      integer, intent(in) :: lo(3)
      real(c_real), target, intent(in) :: fp(lo(1):,lo(2):,lo(3):)
      real(c_real), pointer, intent(inout) :: a(:,:,:)
      a => fp
    end subroutine shift_bound_d3_v
  end subroutine bl_allocate_r3_v

  subroutine bl_allocate_r1_vc(a, lo, hi, ncomp)
    real(c_real), pointer, intent(inout) :: a(:,:)
    integer, intent(in) :: lo(1), hi(1), ncomp
    integer :: n(2)
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    real(c_real), pointer :: fp(:,:)
    n(1:1) = hi - lo + 1
    n(2) = ncomp
    sz = int(n(1),c_size_t) * int(n(2),c_size_t)
    cp = amrex_mempool_alloc(szr*sz)
    call amrex_real_array_init(cp, sz)
    call c_f_pointer(cp, fp, shape=n)
    call shift_bound_d1_vc(fp, lo, a)
  contains
    subroutine shift_bound_d1_vc(fp, lo, a)
      integer, intent(in) :: lo(1)
      real(c_real), target, intent(in) :: fp(lo(1):,1:)
      real(c_real), pointer, intent(inout) :: a(:,:)
      a => fp
    end subroutine shift_bound_d1_vc
  end subroutine bl_allocate_r1_vc

  subroutine bl_allocate_r2_vc(a, lo, hi, ncomp)
    real(c_real), pointer, intent(inout) :: a(:,:,:)
    integer, intent(in) :: lo(2), hi(2), ncomp
    integer :: n(3)
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    real(c_real), pointer :: fp(:,:,:)
    n(1:2) = hi - lo + 1
    n(3) = ncomp
    sz = int(n(1),c_size_t) * int(n(2),c_size_t) * int(n(3),c_size_t)
    cp = amrex_mempool_alloc(szr*sz)
    call amrex_real_array_init(cp, sz)
    call c_f_pointer(cp, fp, shape=n)
    call shift_bound_d2_vc(fp, lo, a)
  contains
    subroutine shift_bound_d2_vc(fp, lo, a)
      integer, intent(in) :: lo(2)
      real(c_real), target, intent(in) :: fp(lo(1):,lo(2):,1:)
      real(c_real), pointer, intent(inout) :: a(:,:,:)
      a => fp
    end subroutine shift_bound_d2_vc
  end subroutine bl_allocate_r2_vc

  subroutine bl_allocate_r3_vc(a, lo, hi, ncomp)
    real(c_real), pointer, intent(inout) :: a(:,:,:,:)
    integer, intent(in) :: lo(3), hi(3), ncomp
    integer :: n(4)
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    real(c_real), pointer :: fp(:,:,:,:)
    n(1:3) = hi - lo + 1
    n(4) = ncomp
    sz = int(n(1),c_size_t) * int(n(2),c_size_t) * int(n(3),c_size_t) &
         * int(n(4),c_size_t)
    cp = amrex_mempool_alloc(szr*sz)
    call amrex_real_array_init(cp, sz)
    call c_f_pointer(cp, fp, shape=n)
    call shift_bound_d3_vc(fp, lo, a)
  contains
    subroutine shift_bound_d3_vc(fp, lo, a)
      integer, intent(in) :: lo(3)
      real(c_real), target, intent(in) :: fp(lo(1):,lo(2):,lo(3):,1:)
      real(c_real), pointer, intent(inout) :: a(:,:,:,:)
      a => fp
    end subroutine shift_bound_d3_vc
  end subroutine bl_allocate_r3_vc

  subroutine bl_deallocate_r1(a)
    real(c_real), pointer, intent(inout) :: a(:)
    integer :: lo(1)
    type(c_ptr) :: cp
    lo = lbound(a)
    cp = c_loc(a(lo(1)))
    call amrex_mempool_free(cp)
    a => Null()
  end subroutine bl_deallocate_r1

  subroutine bl_deallocate_r2(a)
    real(c_real), pointer, intent(inout) :: a(:,:)
    integer :: lo(2)
    type(c_ptr) :: cp
    lo = lbound(a)
    cp = c_loc(a(lo(1),lo(2)))
    call amrex_mempool_free(cp)
    a => Null()
  end subroutine bl_deallocate_r2

  subroutine bl_deallocate_r3(a)
    real(c_real), pointer, intent(inout) :: a(:,:,:)
    integer :: lo(3)
    type(c_ptr) :: cp
    lo = lbound(a)
    cp = c_loc(a(lo(1),lo(2),lo(3)))
    call amrex_mempool_free(cp)
    a => Null()
  end subroutine bl_deallocate_r3

  subroutine bl_deallocate_r4(a)
    real(c_real), pointer, intent(inout) :: a(:,:,:,:)
    integer :: lo(4)
    type(c_ptr) :: cp
    lo = lbound(a)
    cp = c_loc(a(lo(1),lo(2),lo(3),lo(4)))
    call amrex_mempool_free(cp)
    a => Null()
  end subroutine bl_deallocate_r4

  subroutine bl_deallocate_r5(a)
    real(c_real), pointer, intent(inout) :: a(:,:,:,:,:)
    integer :: lo(5)
    type(c_ptr) :: cp
    lo = lbound(a)
    cp = c_loc(a(lo(1),lo(2),lo(3),lo(4),lo(5)))
    call amrex_mempool_free(cp)
    a => Null()
  end subroutine bl_deallocate_r5

  subroutine bl_deallocate_r6(a)
    real(c_real), pointer, intent(inout) :: a(:,:,:,:,:,:)
    integer :: lo(6)
    type(c_ptr) :: cp
    lo = lbound(a)
    cp = c_loc(a(lo(1),lo(2),lo(3),lo(4),lo(5),lo(6)))
    call amrex_mempool_free(cp)
    a => Null()
  end subroutine bl_deallocate_r6

  subroutine bl_allocate_i1(a, lo1, hi1)
    integer, pointer, intent(inout) :: a(:)
    integer, intent(in) :: lo1, hi1
    integer :: n1
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    integer, pointer :: fp(:)
    n1 = max(hi1-lo1+1, 1)
    sz = szi * int(n1,c_size_t)
    cp = amrex_mempool_alloc(sz)
    call c_f_pointer(cp, fp, shape=(/n1/))
    call shift_bound_i1(fp, lo1, a)
!    a(lo1:) => fp  ! some compilers may not support this
  contains
    subroutine shift_bound_i1(fp, lo1, a)
      integer, intent(in) :: lo1
      integer, target, intent(in) :: fp(lo1:)
      integer, pointer, intent(inout) :: a(:)
      a => fp
    end subroutine shift_bound_i1
  end subroutine bl_allocate_i1

  subroutine bl_allocate_i2(a, lo1, hi1, lo2, hi2)
    integer, pointer, intent(inout) :: a(:,:)
    integer, intent(in) :: lo1, hi1, lo2, hi2
    integer :: n1, n2
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    integer, pointer :: fp(:,:)
    n1 = max(hi1-lo1+1, 1)
    n2 = max(hi2-lo2+1, 1)
    sz = szi * int(n1,c_size_t) * int(n2,c_size_t)
    cp = amrex_mempool_alloc(sz)
    call c_f_pointer(cp, fp, shape=(/n1,n2/))
    call shift_bound_i2(fp, lo1, lo2, a)
!    a(lo1:,lo2:) => fp
  contains
    subroutine shift_bound_i2(fp, lo1, lo2, a)
      integer, intent(in) :: lo1, lo2
      integer, target, intent(in) :: fp(lo1:,lo2:)
      integer, pointer, intent(inout) :: a(:,:)
      a => fp
    end subroutine shift_bound_i2
  end subroutine bl_allocate_i2

  subroutine bl_allocate_i3(a, lo1, hi1, lo2, hi2, lo3, hi3)
    integer, pointer, intent(inout) :: a(:,:,:)
    integer, intent(in) :: lo1, hi1, lo2, hi2, lo3, hi3
    integer :: n1, n2, n3
    integer (kind=c_size_t) :: sz
    type(c_ptr) :: cp
    integer, pointer :: fp(:,:,:)
    n1 = max(hi1-lo1+1, 1)
    n2 = max(hi2-lo2+1, 1)
    n3 = max(hi3-lo3+1, 1)
    sz = szi * int(n1,c_size_t) * int(n2,c_size_t) * int(n3,c_size_t)
    cp = amrex_mempool_alloc(sz)
    call c_f_pointer(cp, fp, shape=(/n1,n2,n3/))
    call shift_bound_i3(fp, lo1, lo2, lo3, a)
!    a(lo1:,lo2:,lo3:) => fp
  contains
    subroutine shift_bound_i3(fp, lo1, lo2, lo3, a)
      integer, intent(in) :: lo1, lo2, lo3
      integer, target, intent(in) :: fp(lo1:,lo2:,lo3:)
      integer, pointer, intent(inout) :: a(:,:,:)
      a => fp
    end subroutine shift_bound_i3
  end subroutine bl_allocate_i3

  subroutine bl_deallocate_i1(a)
    integer, pointer, intent(inout) :: a(:)
    integer :: lo(1)
    type(c_ptr) :: cp
    lo = lbound(a)
    cp = c_loc(a(lo(1)))
    call amrex_mempool_free(cp)
    a => Null()
  end subroutine bl_deallocate_i1

  subroutine bl_deallocate_i2(a)
    integer, pointer, intent(inout) :: a(:,:)
    integer :: lo(2)
    type(c_ptr) :: cp
    lo = lbound(a)
    cp = c_loc(a(lo(1),lo(2)))
    call amrex_mempool_free(cp)
    a => Null()
  end subroutine bl_deallocate_i2

  subroutine bl_deallocate_i3(a)
    integer, pointer, intent(inout) :: a(:,:,:)
    integer :: lo(3)
    type(c_ptr) :: cp
    lo = lbound(a)
    cp = c_loc(a(lo(1),lo(2),lo(3)))
    call amrex_mempool_free(cp)
    a => Null()
  end subroutine bl_deallocate_i3

end module amrex_mempool_module

max_stars_repo_name: ylunalin/amrex

max_issues_repo_name: ylunalin/amrex

max_forks_repo_name: ylunalin/amrex


text:
!==============================================================================!
  subroutine User_Mod_End_Of_Compute_Scalar(flow, turb, mult, dt, ini)
!------------------------------------------------------------------------------!
!   This function is called at the end of Compute_Scalar function.             !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Field_Type),      target :: flow
  type(Turb_Type),       target :: turb
  type(Multiphase_Type), target :: mult
  real                          :: dt    ! time step
  integer                       :: ini   ! inner iteration
!==============================================================================!

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
module eval_conductivity
    use conductivity_def

contains

    subroutine conductivity(rq,rho,T,chi,Gamma,eta,mu_e,ionic,Tns,kappa)
        use exceptions_lib
        use constants_def
        use nucchem_def, only: composition_info_type
        use neutron_conductivity
        use electron_conductivity
        use photon_conductivity

        type(conductivity_general_info), pointer :: rq
        real(dp), intent(in) :: rho, T, chi, Gamma, eta, mu_e   ! mu_e is in MeV
        type(composition_info_type), intent(in) :: ionic
        real(dp), intent(in) :: Tns ! superfluid critical temperature, neutrons
        type(conductivity_components), intent(out) :: kappa
        real(dp) :: nn, nion, ne, kF, xF, eF, Gamma_e, kappa_e_pre
        real(dp) :: kappa_n_pre, kfn, mnstar, tau, v, R
        real(dp) :: nu, nu_c, K_opacity, lgrho, alpha
        ! consider putting these in the general info structure
        real(dp), parameter :: eQ_threshold = 1.0e-8_dp
        real(dp), parameter :: sf_reduction_threshold = 1.0e-10_dp
        integer :: ierr
        type(alert) :: scattering_options=alert(scope='conductivity')
        type(failure) :: nimp_error=failure(scope='conductivity')
        type(failure) :: nph_error=failure(scope='conductivity')

        call clear_kappa
        lgrho = log10(rho)

        ! electrons; always included
        ne = rho/amu*ionic%Ye
        kF = (threepisquare*ne)**onethird
        xF = hbar*kF/Melectron/clight
        eF = sqrt(1.0_dp+xF**2)
        Gamma_e = Gamma/ionic%Z53
        kappa_e_pre = onethird*pi**2*boltzmann**2*T*ne/Melectron/eF
        nu = 0.0_dp
        nu_c = 0.0_dp

        ! electron-electron scattering
        select case (rq% ee_scattering_fmla)
        case(icond_sy06)
            nu_c = ee_SY06(ne,T)
        case(icond_pcy)
            nu_c = ee_PCY(ne,T)
        case default
            call scattering_options% report( &
            &   'unknown option for ee scattering: using SY06')
            nu_c = ee_SY06(ne,T)
        end select
        nu = nu + nu_c
        if (nu_c > tiny(1.0_dp)) kappa% ee = kappa_e_pre/nu_c

        ! electron-ion scattering
        nu_c = eion( &
        &   kF,Gamma_e,eta,ionic%Ye,ionic%Z,ionic%Z2,ionic%Z53,ionic%A)
        nu = nu + nu_c
        if (nu_c > tiny(1.0_dp)) kappa% ei = kappa_e_pre/nu_c

        ! electron-impurity scattering
        select case(rq% eQ_scattering_fmla)
        case(icond_eQ_potekhin)
            nu_c = eQ(kF,T,ionic%Ye,ionic%Z,ionic%Z2,ionic%A,ionic%Q)
        case(icond_eQ_page)
            nu_c = eQ_page( &
            &    kF,T,ionic%Ye,ionic%Z,ionic%Z2,ionic%A,ionic%Q)
        case default
            call scattering_options% report( &
            &   'unknown option for eQ scattering: using Potekhin')
            nu_c = eQ(kF,T,ionic%Ye,ionic%Z,ionic%Z2,ionic%A,ionic%Q)
        end select
        nu = nu + nu_c
        if (ionic%Q > eQ_threshold) kappa% eQ = kappa_e_pre/nu_c

        kappa% electron_total = kappa_e_pre/nu

        ! neutrons
        if (ionic% Yn > 0.0_dp) then
            if (rq% include_neutrons) then
                nu = 0.0_dp
                nu_c = 0.0_dp
                nn = rho/amu*ionic% Yn/(1.0-chi)
                kFn = (threepisquare*nn)**onethird
                nion = rho/amu*(1.0-ionic%Yn)/ionic%A
                mnstar = Mneutron
                kappa_n_pre = onethird*pi**2*boltzmann**2*T*nn/mnstar

                ! superfluid reduction factor
                R = 1.0
                if (T < Tns) then
                    tau = T/Tns
                    v = sqrt(1.0-tau)*(1.456-0.157/sqrt(tau) + 1.764/tau)
                    R = (0.4186+sqrt(1.007**2+(0.5010*v)**2))**2.5 * &
                    &   exp(1.456-sqrt(1.456**2+v**2))
                    if (R < sf_reduction_threshold) R = 0.0
                end if

                nu_c = n_imp(nn,nion,T,ionic,ierr)
                if (nimp_error% raised(ierr)) then
                    nu_c = 0.0_dp
                else
                    kappa% nQ = kappa_n_pre/nu_c*R
                    nu = nu + nu_c
                end if

                nu_c = n_phonon(nn,nion,T,ionic, ierr)
                if (nph_error% raised(ierr)) then
                    nu_c = 0.0_dp
                else
                    kappa% np = kappa_n_pre/nu_c*R
                    nu = nu + nu_c
                end if

                if (nu > 0.0_dp) kappa% neutron_total = kappa_n_pre/nu*R
            end if

            ! superfluid phonons
            if (rq% include_superfluid_phonons .and. T < Tns) then
                kappa% sf =  sPh(nn/density_n,nion/density_n,T,ionic)
            end if
        end if

        ! photons
        K_opacity = 0.0_dp
        kappa% kap = 0.0_dp
        if (lgrho < rq% rad_full_off_lgrho) then
            kappa% kap = Rosseland_kappa(rho,T,mu_e,ionic)
            if (lgrho > rq% rad_full_on_lgrho) then
                alpha = (lgrho-rq% rad_full_on_lgrho)/(rq% rad_full_off_lgrho-rq% rad_full_on_lgrho)
                kappa% kap = kappa% kap * (1.0_dp-alpha)
            end if
            K_opacity = 4.0_dp*onethird*arad*clight*T**3/rho/kappa% kap
        end if

        kappa% total = kappa% electron_total  &
        &   + kappa% neutron_total + kappa% sf &
        &   + K_opacity

        contains
        subroutine clear_kappa()
            kappa% total = 0.0_dp
            kappa% ee  = 0.0_dp
            kappa% ei = 0.0_dp
            kappa% eQ = 0.0_dp
            kappa% sf = 0.0_dp
            kappa% kap = 0.0_dp
            kappa% nQ = 0.0_dp
            kappa% np = 0.0_dp
            kappa% electron_total = 0.0_dp
            kappa% neutron_total = 0.0_dp
        end subroutine clear_kappa
    end subroutine conductivity

end module eval_conductivity

max_stars_repo_name: sudhanvalalit/dStar

max_issues_repo_name: sudhanvalalit/dStar

max_forks_repo_name: sudhanvalalit/dStar


! { dg-do run }
! { dg-options "-fcray-pointer" }
!
use iso_c_binding
real target(10)
real pointee(10)
pointer (ipt, pointee)
integer(c_intptr_t) :: int_cptr
real :: x
if (c_sizeof(ipt) /= c_sizeof(int_cptr)) call abort()
if (c_sizeof(pointee) /= c_sizeof(x)*10) call abort()
end


    Program zhesv_example

!     ZHESV Example Program Text

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     .. Use Statements ..
      Use lapack_example_aux, Only: nagf_file_print_matrix_complex_gen_comp
      Use lapack_interfaces, Only: zhesv
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nb = 64, nin = 5, nout = 6
!     .. Local Scalars ..
      Integer :: i, ifail, info, lda, lwork, n
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), b(:), work(:)
      Integer, Allocatable :: ipiv(:)
      Character (1) :: clabs(1), rlabs(1)
!     .. Executable Statements ..
      Write (nout, *) 'ZHESV Example Program Results'
      Write (nout, *)
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n
      lda = n
      lwork = nb*n
      Allocate (a(lda,n), b(n), work(lwork), ipiv(n))

!     Read the upper triangular part of the matrix A from data file

      Read (nin, *)(a(i,i:n), i=1, n)

!     Read b from data file

      Read (nin, *) b(1:n)

!     Solve the equations Ax = b for x
      Call zhesv('Upper', n, 1, a, lda, ipiv, b, n, work, lwork, info)

      If (info==0) Then

!       Print solution

        Write (nout, *) 'Solution'
        Write (nout, 100) b(1:n)

!       Print details of factorization

        Write (nout, *)
        Flush (nout)

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call nagf_file_print_matrix_complex_gen_comp('Upper', &
          'Non-unit diagonal', n, n, a, lda, 'Bracketed', 'F7.4', &
          'Details of the factorization', 'Integer', rlabs, 'Integer', clabs, &
          80, 0, ifail)

!       Print pivot indices

        Write (nout, *)
        Write (nout, *) 'Pivot indices'
        Write (nout, 110) ipiv(1:n)

      Else
        Write (nout, 120) 'The diagonal block ', info, ' of D is zero'
      End If

100   Format ((3X,4(' (',F7.4,',',F7.4,')',:)))
110   Format (1X, 7I11)
120   Format (1X, A, I3, A)
    End Program


text:
use SoybeanClass
implicit none

type(Soybean_) :: soy
real(real64) :: biomass,totalweight

call soy%init(config="Tutorial/obj/realSoybeanConfig.json")

biomass = soy%getBiomass(stem=true)
print *, "Total Dry Weight (stem) : ",biomass*1000.0d0," g"
biomass = soy%getBiomass(leaf=true)
print *, "Total Dry Weight (leaf) : ",biomass*1000.0d0," g"
biomass = soy%getBiomass(leaf=true,stem=true)
print *, "Total Dry Weight (leaf+stem) :/0.20 m^2 ",biomass*1000.0d0*9.50d0," g/m^2"

totalweight = soy%gettotalweight(stem=true)
print *, "Total Dry Weight (stem) : ",totalweight*1000.0d0," g"
totalweight = soy%gettotalweight(leaf=true)
print *, "Total Dry Weight (leaf) : ",totalweight*1000.0d0," g"
totalweight = soy%gettotalweight(leaf=true,stem=true)
print *, "Total Dry Weight (leaf+stem) :/0.20 m^2 ",totalweight*1000.0d0*9.50d0," g/m^2"

print *, "Stem Volume (m^3):",soy%getVolume(stem=true)
print *, "Stem Volume (cm^3):",soy%getVolume(stem=true)*100.0d0*100.0d0*100.0d0
print *, "Leaf Volume (m^3):",soy%getVolume(leaf=true)
print *, "Leaf Volume (cm^3):",soy%getVolume(leaf=true)*100.0d0*100.0d0*100.0d0
!call soy%msh(name="soy")

end

max_stars_repo_name: kazulagi/plantFEM


text:
      SUBROUTINE STEMPERATURE(SOIL_TEMPERATURE,LIQUID_FRACTION,
     &                        TIME_STEP,SN,GC)
      IMPLICIT NONE
      INCLUDE 'pwbm.soil_temperature.h'

      REAL*8 SOIL_TEMPERATURE(SX:LX+1)        !Soil temperature       (C)   [In/Out]
      REAL*8 LIQUID_FRACTION(SX:LX+1)         !Liquid water fraction  (%)   [In/Out]
      REAL*8 TIME                             !Start of calculations  (days)[In]
      REAL*8 TIME_END                         !End of calculations    (days)[In]
      REAL*8 TIME_STEP                        !Current time step             (days)[In/Out]
      INTEGER SN                              !Number of the first node in the grid
      INTEGER GC                              !Counter of converged iterations

      REAL*8 OLD_TIME_STEP,DT_TEMP

      time = 0.0d0
      time_end = 1.0d0
      CALL STEMPERATURE_ATOMIC(SOIL_TEMPERATURE,LIQUID_FRACTION,
     &                         TIME_STEP,TIME,TIME_END,SN,GC)
      OLD_TIME_STEP=TIME_STEP
      DT_TEMP=TIME_END-TIME
      DO WHILE(1.0D-5.LT.DT_TEMP)
        CALL STEMPERATURE_ATOMIC(SOIL_TEMPERATURE,LIQUID_FRACTION,
     &                           DT_TEMP,TIME,TIME_END,SN,GC)
        DT_TEMP=TIME_END-TIME
      ENDDO
      TIME_STEP=OLD_TIME_STEP

      END

      SUBROUTINE STEMPERATURE_ATOMIC(SOIL_TEMPERATURE,
     &                               LIQUID_FRACTION,
     &                               DT,TIME,TIME_END,SN,GC)
      IMPLICIT NONE
      INCLUDE 'pwbm.soil_temperature.h'

      INTEGER SN,GC
      REAL*8 SOIL_TEMPERATURE(SX:LX+1)
      REAL*8 LIQUID_FRACTION(SX:LX+1)

      INTEGER IT,I
      INTEGER FLAG

      REAL*8 TIME,TIME_END,DT
      REAL*8 T(SX:LX+2)
      REAL*8 Wc, Wp, dTemperature, Capp_water, Capp

      REAL*8 DELTA
      REAL*8 C1,B1,D1,R1,L1
      REAL*8 TL(SX:LX+1),   TR(SX:LX+1)
      REAL*8 ATR(SX:LX+1),  BTR(SX:LX+1),  CTR(SX:LX+1)
      REAL*8 ATRNR(SX:LX+1),BTRNR(SX:LX+1),CTRNR(SX:LX+1)
      REAL*8 DTR(SX:LX+1),  DTRNR(SX:LX+1)

      REAL*8 C(SX:LX+1), DC_DT(SX:LX+1)                       !Heat capacity
      REAL*8 H(SX:LX+1), DH_DT(SX:LX+1), D2H_D2T(SX:LX+1)     !Enthalpy due to water
      REAL*8 K(SX:LX+1), DK_DT(SX:LX+1)                       !Thermal conductivity
      REAL*8 Y(SX:LX+1), DY_DT(SX:LX+1), D2Y_D2T(SX:LX+1)     !Liquid water fraction

      REAL*8 SURFACE_TEMPERTURE

      DO WHILE(TIME+DT.LE.TIME_END)
        DO I=SN,LX+1
          T(I)=SOIL_TEMPERATURE(I)
        ENDDO

        FLAG=1
        IT=0
        DO WHILE (FLAG.EQ.1)
          DO I=SN,0
            CALL SNOW_TH_PROPERTIES(T(I),Y(I),DY_DT(I),D2Y_D2T(I),
     &                              H(I),DH_DT(I),D2H_D2T(I),
     &                              C(I),DC_DT(I),
     &                              K(I),DK_DT(I),
     &               LAYER_W1(I),LAYER_CF(I),LAYER_CT(I),LAYER_LF(I),
     &               LAYER_LT(I),LAYER_B(I),LAYER_Tp(I),LAYER_Tp_1(I))
          ENDDO
          DO I=2,LX+1
            CALL SOIL_TH_PROPERTIES(T(I),Y(I),DY_DT(I),D2Y_D2T(I),
     &                              H(I),DH_DT(I),D2H_D2T(I),
     &                              C(I),DC_DT(I),
     &                              K(I),DK_DT(I),
     &               LAYER_W1(I),LAYER_CF(I),LAYER_CT(I),LAYER_LF(I),
     &               LAYER_LT(I),LAYER_B(I),LAYER_Tp(I),LAYER_Tp_1(I))

c            print*, I, LAYER_LF(I), LAYER_LT(I)
c            if(I .eq. 3) then
c            if(T(I) .lt. 0.0d0) then
c             write(103,*) T(I),Y(I),C(I),K(I),LAYER_CF(I),LAYER_LF(I)
c            else
c             write(104,*) T(I),Y(I),C(I),K(I),LAYER_CT(I),LAYER_LT(I)
c            endif
c            endif

c            if(I .eq. 10) then
c            if(T(I) .lt. 0.0d0) then
c             write(105,*) T(I),Y(I),C(I),K(I),LAYER_CF(I),LAYER_LF(I)
c            else
c             write(106,*) T(I),Y(I),C(I),K(I),LAYER_CT(I),LAYER_LT(I)
c            endif
c            endif
          ENDDO

          !Snow model
          DO I=SN,-1
            TL(I)=(K(I)+K(I+1))*DT*DX_1(I)*0.5D0

            BTR(I)  = TL(I)
            CTR(I)  =-TL(I)
            ATR(I+1)=-TL(I)

            BTRNR(I)  = TL(I)
            CTRNR(I)  =-TL(I)
            ATRNR(I+1)=-TL(I)
          ENDDO
          BTR(0)=0.0D0
          BTRNR(0)=0.0D0

          !Flux
          CTR(0)= 1.0D0
          CTRNR(0)= 1.0D0
          !Mortar
          BTR(1)= 0.0D0
          BTRNR(1)= 0.0D0
          ATR(1)= 1.0D0
          ATRNR(1)= 1.0D0
          CTR(1)=-1.0D0
          CTRNR(1)=-1.0D0
          DTR(1)= 0.0D0
          !Flux
          ATR(2)=-1.0D0; ATRNR(2)=-1.0D0

          ATR(SN)=0.0D0
          ATRNR(SN)=0.0D0
          !Soil
          DO I=2,LX
            TL(I) = (K(I)+K(I+1))*DT*DX_1(I-1)*0.5D0
            TR(I) = (T(I)-T(I+1))*DT*DX_1(I-1)*0.5D0

            BTR(I)  = TL(I)
            CTR(I)  =-TL(I)
            ATR(I+1)=-TL(I)

            BTRNR(I)  = TL(I)+TR(I)*DK_DT(I)
            CTRNR(I)  =-TL(I)+TR(I)*DK_DT(I+1)
            ATRNR(I+1)=-TL(I)-TR(I)*DK_DT(I)
          ENDDO
          BTR(LX+1)=0.0D0
          BTRNR(LX+1)=0.0D0
          CTR(LX+1)=0.0D0
          CTRNR(LX+1)=0.0D0

          !Snow
          DO I=SN,-1
            BTR(I+1)  =BTR(I+1)  +TL(I)
            BTRNR(I+1)=BTRNR(I+1)+TL(I)
          ENDDO
          !Soil
          DO I=2,LX
            BTR(I+1)  =BTR(I+1)  +TL(I)
            BTRNR(I+1)=BTRNR(I+1)+TL(I)-TR(I)*DK_DT(I+1)
          ENDDO

          !Snow
          DO I=SN,0
            BTR(I)=BTR(I)+C(I)*DIAG(I)
            BTRNR(I)=BTRNR(I)+C(I)*DIAG(I)
            DTR(I)=C(I)*DIAG(I)*SOIL_TEMPERATURE(I)
          ENDDO
          !Soil
          DO I=2,LX+1
            dTemperature=T(I)-SOIL_TEMPERATURE(I)
            IF(ABS(dTemperature)<dlt) THEN                        !Second term in the apparent heat capacity due to freezing water
              Capp_water=DH_DT(I)                                 !if temperature in the soil layer is the same we approximate it analytically
              BTRNR(I)=BTRNR(I)+D2H_D2T(I)*DIAG(I-1)*dTemperature
            ELSE
              Wp=LAYER_W0(I)*LIQUID_FRACTION(I)*La*TSCALE          !Enthalpy due to unfrozen water at the previous step
              Wc=H(I)                                             !Enthalpy due to unfrozen water at the current step
              Capp_water=(Wc-Wp)/dTemperature                     !Otherwise we use an approximation
              BTRNR(I)=BTRNR(I)- (Capp_water-DH_DT(I))*DIAG(I-1)
            ENDIF
            Capp=C(I)+Capp_water                                  !Apparent heat capacity

            BTR(I)=BTR(I)+Capp*DIAG(I-1)
            BTRNR(I)=BTRNR(I)+Capp*DIAG(I-1)
     &              +DC_DT(I)*DIAG(I-1)*(T(I)-SOIL_TEMPERATURE(I))

            DTR(I)=Capp*DIAG(I-1)*SOIL_TEMPERATURE(I)
          ENDDO

          DO I=2,LX+1
              DTR(I)=DTR(I)+(LAYER_W1(I)-LAYER_W0(I))*DT*DIAG(I-1)
     &                                             *La*TSCALE
          ENDDO

          R1=SURFACE_TEMPERTURE(TIME+DT)

          IF(IT<3) THEN  ! yyh: what does this mean?
            C1=CTR(SN)
            B1=BTR(SN)
            D1=DTR(SN)
            CTR(SN)=0.0D0
            BTR(SN)=1.0D0
            DTR(SN)=R1

            CALL TRIDIAG(SN,SX,LX+1,ATR,BTR,CTR,DTR)
            L1=D1-(B1*DTR(SN)+C1*DTR(SN+1))

            DELTA=1.0D0
!            DELTA=0.0D0
!            DO I=SN,LX+1
!              delta=DMAX1(DABS(DTR(I)-T(I)),DELTA)
!            ENDDO

            DO I=SN,LX+1
              T(I)=DTR(I)
            ENDDO
            T(LX+2)=L1
          ELSE
            DTRNR(SN)=BTR(SN)*T(SN)+CTR(SN)*T(SN+1)+T(LX+2)-DTR(SN)
            DO I=SN+1,LX
              DTRNR(I)=ATR(I)*T(I-1)+BTR(I)*T(I)
     &                +CTR(I)*T(I+1)-DTR(I)
            ENDDO
            DTRNR(LX+1)=ATR(LX+1)*T(LX)+BTR(LX+1)*T(LX+1)-DTR(LX+1)

            C1=CTRNR(SN)
            B1=BTRNR(SN)
            D1=DTRNR(SN)
            CTRNR(SN)=0.0D0
            BTRNR(SN) =1.0D0
            DTRNR(SN) =T(SN)-R1
            CALL TRIDIAG(SN,SX,LX+1,ATRNR,BTRNR,CTRNR,DTRNR)
            L1=D1-(B1*DTRNR(SN)+C1*DTRNR(SN+1))

            DELTA=0.0D0
            DO I=SN,LX
              T(I)=T(I)-DTRNR(I)
              delta=DMAX1(DABS(DTRNR(I)),DELTA)
            ENDDO
              T(LX+2)=T(LX+2)-L1
          ENDIF

C          WRITE(*,'(F10.6,F16.13 ,I3,F20.16,I4)') TIME+DT,DELTA,IT,DT,GC
          IF(DELTA.LT.MAX_DELTA) THEN
            !PRINT *,time,IT,DT
            TIME=TIME+DT
            GC=GC+1
            IF(GC.EQ.5) THEN
              DT=DT*2.0D0
              GC=0
            ENDIF
            DT=DMIN1(DT,MAX_DT)
            IT=0
            FLAG=0
          ELSE
            IT=IT+1
          ENDIF
          IF(IT.EQ.NUM_MAX_ITER) THEN
            DT=DT/2.0D0
            DO I=SN,LX+1
              T(I)=SOIL_TEMPERATURE(I)
            ENDDO
            GC=0
            IT=0
          ENDIF
        ENDDO

        DO I=SN,0
          SOIL_TEMPERATURE(I)=T(I)
          CALL SNOW_TH_PROPERTIES(SOIL_TEMPERATURE(I),
     &                            LIQUID_FRACTION(I),
     &                            DY_DT(I),D2Y_D2T(I),
     &                            H(I),DH_DT(I),D2H_D2T(I),
     &                            C(I),DC_DT(I),
     &                            K(I),DK_DT(I),
     &             LAYER_W1(I), LAYER_CF(I),LAYER_CT(I),LAYER_LF(I),
     &             LAYER_LT(I),LAYER_B(I), LAYER_Tp(I),LAYER_Tp_1(I))
          LAYER_L(I)=LAYER_W1(I)*LIQUID_FRACTION(I)
          LAYER_I(I)=LAYER_W1(I)*(1.0D0-LIQUID_FRACTION(I))
          LAYER_K(I)=K(I)
          LAYER_C(I)=C(I)/TSCALE
        ENDDO
        SOIL_TEMPERATURE(1)=T(1)
          DO I=2,LX+1
          SOIL_TEMPERATURE(I)=T(I)
          CALL SOIL_TH_PROPERTIES(SOIL_TEMPERATURE(I),
     &                            LIQUID_FRACTION(I),
     &                            DY_DT(I),D2Y_D2T(I),
     &                            H(I),DH_DT(I),D2H_D2T(I),
     &                            C(I),DC_DT(I),
     &                            K(I),DK_DT(I),
     &             LAYER_W1(I),LAYER_CF(I),LAYER_CT(I),LAYER_LF(I),
     &             LAYER_LT(I),LAYER_B(I),LAYER_Tp(I),LAYER_Tp_1(I))
          LAYER_L(I)=LAYER_W1(I)*LIQUID_FRACTION(I)
          LAYER_I(I)=LAYER_W1(I)*(1.0D0-LIQUID_FRACTION(I))
          LAYER_K(I)=K(I)
          LAYER_C(I)=C(I)/TSCALE
        ENDDO
        DO I=SX,SN-1
          SOIL_TEMPERATURE(I)=SOIL_TEMPERATURE(SN)
          LIQUID_FRACTION(I)=LIQUID_FRACTION(SN)
        ENDDO
      ENDDO
      END SUBROUTINE STEMPERATURE_ATOMIC

      SUBROUTINE SNOW_TH_PROPERTIES(T,Y,DY_DT,D2Y_D2T,
     &                              H,DH_DT,D2H_D2T,
     &                              C,DC_DT,
     &                              K,DK_DT,
     &                              W,CF,CT,LF,LT,B,Tp,Tp_1)
      IMPLICIT NONE
      INCLUDE 'pwbm.soil_temperature.h'

      REAL*8 T,W,LT,LF,CF,CT,B,Tp,Tp_1
      REAL*8      Y,DY_DT,D2Y_D2T
      REAL*8      H,DH_DT,D2H_D2T
      REAL*8      C,DC_DT
      REAL*8      K,DK_DT

      C       =CF      *TSCALE
      K       =LF

      DK_DT   =0.0D0
      DC_DT   =0.0D0

      Y       =1.0D0
      DY_DT   =0.0D0
      D2Y_D2T =0.0D0
      H       =0.0D0
      DH_DT   =0.0D0
      D2H_D2T =0.0D0

      END SUBROUTINE SNOW_TH_PROPERTIES

      SUBROUTINE SOIL_TH_PROPERTIES(T,Y,DY_DT,D2Y_D2T,
     &                              H,DH_DT,D2H_D2T,
     &                              C,DC_DT,
     &                              K,DK_DT,
     &                              W,CF,CT,LF,LT,B,Tp,Tp_1)
      IMPLICIT NONE
      INCLUDE 'pwbm.soil_temperature.h'

      REAL*8 T,W,LT,LF,CF,CT,B,Tp,Tp_1

      REAL*8      Y,DY_DT,D2Y_D2T
      REAL*8      H,DH_DT,D2H_D2T
      REAL*8      C,DC_DT
      REAL*8      K,DK_DT

      REAL*8 F1, DF1, DDF1, KWKI
      REAL*8 E2, Q1,I1,I2,Z1,TS_2,Tp_2,TS

      TS=T-XD-Tp
      IF(TS.LT.0.0D0) THEN
        TS_2=1.0D0/(TS*TS)
        Tp_2=Tp_1*Tp_1

        Q1=TS_2*XS_2
        E2=DEXP(-Q1)
        F1=TS*E2+Tp
        DF1=E2*(1.0D0+2.0D0*Q1)
        DDF1=DF1*(2.0D0*Q1/TS)+E2*(-4.0D0*Q1/TS)

        Z1=DLOG(F1*Tp_1)
        Y         =DEXP(B*Z1) !(F1/Tp)**B
        I1=Y*Tp/F1
        I2=I1*Tp/F1

        DY_DT     =B*Tp_1*dF1*I1
        D2Y_D2T   =B*Tp_1*ddF1*I1+B*(B-1.0D0)*Tp_2*dF1*dF1*I2

        C        =(CF+W*Y*(CW-CI))      *TSCALE
        DC_DT    =W*DY_DT*(CW-CI)       *TSCALE
        H        =W*Y*La                *TSCALE
        DH_DT    =W*DY_DT*La            *TSCALE
        D2H_D2T  =W*D2Y_D2T*La          *TSCALE

        KWKI     =DEXP(Y*W*LKWKI)
        K        =LF*KWKI
c        K = LF*KWKI*2.0d0
c        write(103,*) T, Y, KWKI, K, C
        DK_DT    =LF*KWKI*W*LKWKI*DY_DT
c        DK_DT = LF*KWKI*W*LKWKI*DY_DT*2.0d0
      ELSE
        Y        =1.0D0
        DY_DT    =0.0D0
        D2Y_D2T  =0.0D0

        C        =CT                    *TSCALE
        DC_DT    =0.0D0
        H        =W*La                  *TSCALE
        DH_DT    =0.0D0
        D2H_D2T  =0.0D0

        K         = LT
        DK_DT     = 0.0D0
c        write(104,*) T, K, C
      ENDIF
      END SUBROUTINE SOIL_TH_PROPERTIES

      SUBROUTINE SOILSAT(T,Y,B,Tp)
      IMPLICIT NONE
      INCLUDE 'pwbm.soil_temperature.h'

      REAL*8 Y,T,B,Tp,Tp_1
      REAL*8 F1,E2,Q1,Z1,TS_2,TS

      Tp_1=1.0D0/Tp
      TS=T-XD-Tp
      IF(TS.LT.0.0D0) THEN
        TS_2=1.0D0/(TS*TS)

        Q1=TS_2*XS_2
        E2=DEXP(-Q1)
        F1=TS*E2+Tp
        Z1=DLOG(F1*Tp_1)
        Y =DEXP(B*Z1) !(F1/Tp)**B
      ELSE
        Y=1.0D0
      ENDIF
      END SUBROUTINE SOILSAT

      !Solution of a tridiagonal matrix
      SUBROUTINE TRIDIAG(M,N0,N1,A,B,C,D)
      IMPLICIT NONE
      INTEGER M,N0,N1,I,K
      REAL*8 A(N0:N1),B(N0:N1),C(N0:N1),D(N0:N1), XM

      DO K=M+1,N1
        XM=A(K)/B(K-1)
        B(K)=B(K)-XM*C(K-1)
        D(K)=D(K)-XM*D(K-1)
      ENDDO

      D(N1)=D(N1)/B(N1)
      DO I=M+1,N1
        K=N1+M-I
        D(K)=(D(K)-C(K)*D(K+1))/B(K)
      ENDDO
      END SUBROUTINE TRIDIAG

      !Subroutine to compute an upper boundary condition
      REAL*8 FUNCTION SURFACE_TEMPERTURE(TIME)
      IMPLICIT NONE
      INCLUDE 'pwbm.soil_temperature.h'
      REAL*8 TIME

      SURFACE_TEMPERTURE=SURFACE_TEMPERATURE_BEGIN+
     &   (SURFACE_TEMPERATURE_END-SURFACE_TEMPERATURE_BEGIN)*TIME
C      SURFACE_TEMPERTURE=1.0d0
      END FUNCTION SURFACE_TEMPERTURE

      SUBROUTINE UPDATE_PROPERTIES(X, SN,
     &                  SNOW_HEIGHT,SNOW_THCND,SNOW_HEATC,
     &                  PEAT_FRAC,LT_PEAT_DRY,LT_PEAT_SAT,
     &                  LF_PEAT_DRY,LF_PEAT_SAT,
     &                  SOIL_POROSITY,THCND_SOIL_DRY,
     &                  THCND_SOLIDS,HEATC_SOIL_DRY,
     &                  LIQUID_FRACTION,SATURATION0,SATURATION1,
     &                  SOIL_FREEZING_PD,SOIL_CLEANESS,IWRITE_FLAG)
      IMPLICIT NONE
      INCLUDE 'pwbm.soil_temperature.h'

      REAL*8 SNOW_HEIGHT                               !Snow height
      REAL*8 SNOW_THCND, SNOW_HEATC                    !Snow thermal conductivity and heat capacity
      REAL*8 SOIL_POROSITY(2:LX+1)                      !Soil porosity
      REAL*8 THCND_SOIL_SAT,KE
      REAL*8 PEAT_FRAC(2:LX+1)
      REAL*8 LT_PEAT_DRY,LT_PEAT_SAT       ! thermal conductivity of unfrozen peat (dry & saturation)
      REAL*8 LF_PEAT_DRY,LF_PEAT_SAT       ! thermal conductivity of frozen peat (dry & saturation)

      REAL*8 THCND_SOLIDS(2:LX+1)     ! Thermal conductivity of soil solids
      REAL*8 HEATC_SOIL_DRY(2:LX+1),THCND_SOIL_DRY(2:LX+1)   !Thermal conductivity of dry peat and soil
      REAL*8 X(SX:LX)
      INTEGER SN
      INTEGER I,J
      REAL*8 DX
      INTEGER SOIL_RLAYER                              !Number of layer in a rootzone
      INTEGER IWRITE_FLAG

      REAL*8 LIQUID_FRACTION(SX:LX+1)
      REAL*8 SATURATION0(SX:LX+1),SATURATION1(SX:LX+1)

      REAL*8 SOIL_CLEANESS(2:LX+1), SOIL_FREEZING_PD(2:LX+1)

      !Provided the snow height,
      !this code analyzes the grid, and inserts a new node representing the snow surface
      IF(1.0D-5.LT.SNOW_HEIGHT) THEN
        IF(DABS(XREF(SX))<SNOW_HEIGHT) THEN
          X(SX)=-SNOW_HEIGHT
          SN=SX
        ELSE
          I=SX
          DO WHILE (I<0)
            IF((DABS(XREF(I+1))<SNOW_HEIGHT).AND.
     &               (SNOW_HEIGHT<=DABS(XREF(I))+1.0d-10)) THEN
              IF((DABS(DABS(XREF(I+1))-SNOW_HEIGHT)<
     &                 DABS(XREF(I+1)-XREF(I))*0.5D0).AND.(I<-1)) THEN
                X(I+1)=-SNOW_HEIGHT
                SN=I+1
              ELSE
                X(I)=-SNOW_HEIGHT
                SN=I
              ENDIF
              I=1
            ELSE
              I=I+1
            ENDIF
          ENDDO
        ENDIF
      ELSE
        SN=2
      ENDIF

      DO J=SX,SN-1
        X(J)=XREF(J)
      ENDDO
      DO J=SN+1,2
        X(J)=XREF(J)
      ENDDO
      DO J=2,LX
        X(J)=XREF(J)
      ENDDO
      !Computes distances between each node in the grid in order to speed up consequent computations
      DO I=SN,1
        IF (I.EQ.SN) THEN
          DX=(X(SN+1)-X(SN))/2.0D0
        ELSEIF (I.EQ.LX) THEN
          DX=(X(LX)-X(LX-1))/2.0D0
        ELSEIF((I.LT.LX).AND.(SX.LT.I)) THEN
          DX=(X(I)-X(I-1))/2.0D0+(X(I+1)-X(I))/2.0D0
        ENDIF
        DIAG(I)=DX
        IF(I<LX) DX_1(I)=1.0D0/(X(I+1)-X(I))
      ENDDO

      !Update the snow thermal properties
      DO I=SN,0
        LAYER_CF(I)=SNOW_HEATC                                !Heat capacity of snow/ice
        LAYER_LF(I)=SNOW_THCND                                !thermal conductivity of snow
        LAYER_S(I)=0.0D0
        LAYER_A(I)=0.0D0
        LAYER_W0(I)=0.0D0                                      !Water content=0,No phase change in snow
        LAYER_W1(I)=0.0D0                                      !Water content=0,No phase change in snow
        LAYER_CT(I)=LAYER_CF(I)                               !Heat capacity is the same
        LAYER_B(I)=-2.0D0                                     !Parameterization of the unfrozen water content
        LAYER_Tp(I)=-0.03D0                                   !Freezing point depression
        LAYER_Tp_1(I)=1.0D0/LAYER_Tp(I)                       !Need it to speed up computations
        LAYER_LT(I)=LAYER_LF(I)                               !Thermal conductivity of frozen soil
      ENDDO

      !Update thermal properties at each soil layer
       DO I=1,LX
        LAYER_S(I+1)=1.0D0-SOIL_POROSITY(I+1)                    !Volumetric content of soil particles
        LAYER_A(I+1)=SOIL_POROSITY(I+1)*(1.0D0-SATURATION1(I+1)) !Volumetric air content

        LAYER_W0(I+1)=SOIL_POROSITY(I+1)*SATURATION0(I+1)        !Volumetric water content at the previous time step
        LAYER_W1(I+1)=SOIL_POROSITY(I+1)*SATURATION1(I+1)        !Volumetric water content at the previous time step

       ! treat organic and mineral soil separately
        if(PEAT_FRAC(I+1) .ge. 0.5d0) then
           ! thawed peat soil
           LAYER_LT(I+1) = LT_PEAT_DRY+(LT_PEAT_SAT-LT_PEAT_DRY)*
     &       SATURATION1(I+1)**2.0d0
!            LAYER_LT(I+1) = LT_PEAT_DRY+(LT_PEAT_SAT-LT_PEAT_DRY)*
!    &       SATURATION1(I+1)
           if(LAYER_LT(I+1) .lt. LT_PEAT_DRY .or. LAYER_LT(I+1) .gt.
     &       LT_PEAT_SAT) then
             print*, LAYER_LT(I+1)
             stop
           endif
           ! frozen peat soil
           LAYER_LF(I+1) = LF_PEAT_DRY*(LF_PEAT_SAT/LF_PEAT_DRY)
     &       **SATURATION1(I+1)
           if(LAYER_LF(I+1) .lt. LF_PEAT_DRY .or. LAYER_LF(I+1) .gt.
     &       LF_PEAT_SAT) then
             print*, LAYER_LF(I+1)
             stop
           endif
        else
           !thawed mineral soil
           THCND_SOIL_SAT=THCND_SOLIDS(I+1)**(1.0-SOIL_POROSITY(I+1))*
     &      (KW**SOIL_POROSITY(I+1))

          if(SATURATION1(I+1) .lt. 0.1) then
             KE = 0.0d0
          else
             KE = log10(SATURATION1(I+1))+1.0
          endif

           if(KE .lt. 0.0 .or. KE .gt. 1.0d0) then
              print*, KE
              stop
           endif

           LAYER_LT(I+1)=THCND_SOIL_DRY(I+1)*(1.0-KE)+
     &       KE*THCND_SOIL_SAT  ! thermal conductivity for thawed soil

           ! frozen mineral soil
           THCND_SOIL_SAT=THCND_SOLIDS(I+1)**(1.0-SOIL_POROSITY(I+1))*
     &       (KI**SOIL_POROSITY(I+1))
           LAYER_LF(I+1)=THCND_SOIL_DRY(I+1)*(1.0-SATURATION1(I+1))+
     &       THCND_SOIL_SAT*SATURATION1(I+1)
        endif

        ! heat capacity
        LAYER_CT(I+1)=HEATC_SOIL_DRY(I+1)+CW*LAYER_W1(I+1)!heat capacity of thawed peat
        LAYER_CF(I+1)=HEATC_SOIL_DRY(I+1)+CI*LAYER_W1(I+1)  ! changed by YYH

        !for the base rock (>3.25 m)
        if(I .gt. 15) then
           LAYER_LT(I+1) = 1.5d0
           LAYER_LF(I+1) = 2.0d0
           LAYER_CT(I+1) = 2.5d0
           LAYER_CF(I+1) = 2.2d0
        endif

c        if(iwrite_flag .eq. 1) then
c           print*, I,SOIL_POROSITY(I+1),SATURATION1(I+1),LAYER_CT(I+1),
c     &       LAYER_CF(I+1),LAYER_LT(I+1),LAYER_LF(I+1)
c        endif

        LAYER_B(I+1)=SOIL_CLEANESS(I+1)                        !Parameterization of the unfrozen water content
        LAYER_Tp(I+1)=SOIL_FREEZING_PD(I+1)                    !Freezing point depression

        !For testing agains the so-called classicall Stefan solution uncomment the following five lines
        !The soil is considered to be homogeneous
        LAYER_Tp_1(I+1)=1.0D0/LAYER_Tp(I+1)               !Need it to speed up computations
c        print*, I, LAYER_LT(I), LAYER_LF(I)
      ENDDO

      !Just for the sake of pretty looking graphics
      DO I=SX,SN-1
        LAYER_CF(I)=-1.0D0                                    !Heat capacity of ice
        LAYER_LF(I)=-1.0D0                                    !Thermal conductivity of thawed soil
        LAYER_K(I)=-1.0D0
        LAYER_C(I)=-1.0D0
                                                              !!! All other parameters are fictitious and introduced to keep computing algorithm efficient
        LAYER_S(I)=0.0D0
        LAYER_A(I)=0.0D0
        LAYER_W0(I)=0.0D0                                     !Porosity=0,No phase change
        LAYER_W1(I)=0.0D0                                     !Porosity=0,No phase change
        LAYER_CT(I)=LAYER_CF(I)                               !Heat capacity of ice
        LAYER_B(I)=-0.7D0                                     !Parameterization of the unfrozen water content
        LAYER_Tp(I)=-0.03D0                                   !Freezing point depression
        LAYER_Tp_1(I)=1.0D0/LAYER_Tp(I)                       !Need it to speed up computations
        LAYER_LT(I)=LAYER_LF(I)                               !Thermal conductivity of frozen soil
      ENDDO

      END SUBROUTINE UPDATE_PROPERTIES

      SUBROUTINE CHECK_SOILFROZEN(LIQ_FRACTION,FLG)
      Implicit none
      INCLUDE 'pwbm.soil_temperature.h'
      REAL*8 LIQ_FRACTION(SX:LX+1)
      INTEGER I, FLG

      FLG=1

      I=2
      DO WHILE(I.LE.LX)
        IF(LIQ_FRACTION(I)>CR_LIQ_FRACTION)THEN
          FLG=0
          I=LX+1
        ELSE
          I=I+1
        ENDIF
      ENDDO

      END SUBROUTINE CHECK_SOILFROZEN

      SUBROUTINE INITIALIZE_SOIL(SOIL_TEMPERATURE,LIQUID_FRACTION,
     &                           SOIL_DT,GC,SN,
     &                           SOIL_TIME,SOIL_FRONT,SOIL_NODES)
      IMPLICIT NONE
      INCLUDE 'pwbm.soil_temperature.h'
      REAL*8  SOIL_TEMPERATURE(SX:LX+1)                       !Temperature
      REAL*8  LIQUID_FRACTION(SX:LX+1)                        !Liquid water fraction
      REAL*8  SOIL_DT                                         !Automatically updated timestep
      REAL*8  SNOW_AGE                                        !Snow age
      INTEGER GC                                              !Automatically updated variable
      INTEGER SN                                              !Number of snow layers + 1
      REAL*8  SOIL_TIME                                       !Time in soil, for testing
      REAL*8  SOIL_FRONT                                      !Location of the freezing front closest to the surface

      INTEGER I
      REAL*8 X(SX:LX),DX,SOIL_NODES(SX:LX)                                     !Local variables

      DO I=SX,-1
        XREF(I)=-0.2D0*(DEXP(0.3D0*DBLE(ABS(I)))-1.0D0)
      ENDDO
      XREF(0)=0.0D0
      XREF(1)  =  0.00D0
      XREF(2:lx) = soil_nodes(2:lx)

      DO I=SX,LX
        X(I)=XREF(I)
      ENDDO

      DO I=SX,LX
        IF (I.EQ.SX) THEN
          DX=(X(SX+1)-X(SX))/2.0D0
        ELSEIF (I.EQ.LX) THEN
          DX=(X(LX)-X(LX-1))/2.0D0
        ELSEIF((I.LT.LX).AND.(SX.LT.I)) THEN
          DX=(X(I)-X(I-1))/2.0D0+(X(I+1)-X(I))/2.0D0
          ENDIF
        DIAG(I)=DX
        IF(I<LX) DX_1(I)=1.0D0/(X(I+1)-X(I))
      ENDDO

      DO I=SX,LX+1
        SOIL_TEMPERATURE(I)=1.0D0
        LIQUID_FRACTION(I) =1.0D0
      ENDDO

      SOIL_TIME =0.0D0
      SOIL_FRONT=0.0D0

      GC=0
      SOIL_DT=MAX_DT
      SN=0
      END SUBROUTINE INITIALIZE_SOIL


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 04.02.01.01/02                        *
C  *    TEST TITLE : Appearance of individual attributes   *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      CALL INITGL ('04.02.01.01/02')
      CALL PLATTR ('i')
      CALL ENDIT

      END


text: I like to play with my kitty and other dogs.  I love chewin stuff. I love destroying stuffed animals and other toys, hence my nickname. I was adopted from the Yolo County SPCA SPCA when I was a wee pup and now my person and I foster other dogs for the SPCA. Life is good!