text:
!
!  Copyright (c) 2019 National Technology & Engineering Solutions of
!  Sandia, LLC (NTESS). Under the terms of Contract DE-NA0003525 with
!  NTESS, the U.S. Government retains certain rights in this software.
!
   SUBROUTINE Z36_Z40
     USE NWC_DATABASE
     USE NWC_2000
     IMPLICIT NONE

     !This subroutine contains elements 36-40 (krypton, rubidium, strontium,
     !yttrium, and zirconium) and their naturally occuring isotopes.
     !Created with the Current Nuclear Wallet Cards (2004)
     !on May 26, 2004 by Karen Kajder.
     !
     !Modified to latest NWC electronic version values and ENDF/B-VII cross
     !sections by K. Russell DePriest on December 14, 2011.
     !
     ! - ENDF/B-VII Release 1 cross sections (that is, .80c) made the default
     !   on 05/09/2014 by K. Russell DePriest.
     !
     ! - ENDF/B-VIII.0 cross sections at 293.6 K (that is, .00c) made the default
     !   on 07/22/2020 by K. Russell DePriest.
     !
     ! - Nuclear Wallet Card values verified on 07/22/2020 using
     !   "Nuclear Wallet Cards database version of 7/10/2019" on nndc.bnl.gov website.

     !Local variables
     REAL::holder, mass_of_element
     INTEGER:: i

     !Number of Krypton Isotopes
     kr_size = 6

     krypton_array(1)%z = 36
     krypton_array(1)%a = 78
     krypton_array(1)%symbol ="Kr-78 "
     krypton_array(1)%mass_defect = -74.1782
     krypton_array(1)%atom_percent = 0.355
     krypton_array(1)%zaid = "36078.00c"

     krypton_array(2)%z = 36
     krypton_array(2)%a = 80
     krypton_array(2)%symbol ="Kr-80 "
     krypton_array(2)%mass_defect = -77.8933
     krypton_array(2)%atom_percent = 2.286
     krypton_array(2)%zaid = "36080.00c"

     krypton_array(3)%z = 36
     krypton_array(3)%a = 82
     krypton_array(3)%symbol ="Kr-82 "
     krypton_array(3)%mass_defect = -80.5917
     krypton_array(3)%atom_percent = 11.593
     krypton_array(3)%zaid = "36082.00c"

     krypton_array(4)%z = 36
     krypton_array(4)%a = 83
     krypton_array(4)%symbol ="Kr-83 "
     krypton_array(4)%mass_defect = -79.9906
     krypton_array(4)%atom_percent = 11.500
     krypton_array(4)%zaid = "36083.00c"

     krypton_array(5)%z = 36
     krypton_array(5)%a = 84
     krypton_array(5)%symbol ="Kr-84 "
     krypton_array(5)%mass_defect = -82.4393
     krypton_array(5)%atom_percent = 56.987
     krypton_array(5)%zaid = "36084.00c"

     krypton_array(6)%z = 36
     krypton_array(6)%a = 86
     krypton_array(6)%symbol ="Kr-86 "
     krypton_array(6)%mass_defect = -83.2656
     krypton_array(6)%atom_percent = 17.279
     krypton_array(6)%zaid = "36086.00c"

     !Calculate isotopic mass for Krypton
     DO i=1,kr_size
          holder = krypton_array(i)%a + (krypton_array(i)%mass_defect/MeV_amu)

          krypton_array(i)%isotopic_mass = holder
     END DO

     !Set mass and holder to zero.
     kr_mass = 0.0
     holder = 0.0

     !Calculate mass of Krypton
     DO i=1, kr_size
          kr_mass = holder + (krypton_array(i)%isotopic_mass* &
          &  krypton_array(i)%atom_percent)/100

          holder = kr_mass
     END DO

     !Calculate weight percent of Krypton
     DO i=1, kr_size
          krypton_array(i)%weight_percent = (krypton_array(i)%isotopic_mass &
          &   * krypton_array(i)%atom_percent)/kr_mass
     END DO

     !Number of Rubidium Isotopes
     rb_size = 2

     rubidium_array(1)%z = 37
     rubidium_array(1)%a = 85
     rubidium_array(1)%symbol ="Rb-85 "
     rubidium_array(1)%mass_defect = -82.1673
     rubidium_array(1)%atom_percent = 72.17
     rubidium_array(1)%zaid = "37085.00c"

     rubidium_array(2)%z = 37
     rubidium_array(2)%a = 87
     rubidium_array(2)%symbol ="Rb-87 "
     rubidium_array(2)%mass_defect = -84.5977
     rubidium_array(2)%atom_percent = 27.83
     rubidium_array(2)%zaid = "37087.00c"

     !Calculate isotopic mass for Rubidium
     DO i=1,rb_size
          holder = rubidium_array(i)%a + (rubidium_array(i)%mass_defect/MeV_amu)

          rubidium_array(i)%isotopic_mass = holder
     END DO

     !Set mass and holder to zero.
     rb_mass = 0.0
     holder = 0.0

     !Calculate mass of Rubidium
     DO i=1, rb_size
          rb_mass = holder + (rubidium_array(i)%isotopic_mass* &
          &  rubidium_array(i)%atom_percent)/100

          holder = rb_mass
     END DO

     !Calculate weight percent of Rubidium
     DO i=1, rb_size
          rubidium_array(i)%weight_percent = (rubidium_array(i)%isotopic_mass &
          &   * rubidium_array(i)%atom_percent)/rb_mass
     END DO

     !Number of Strontium Isotopes
     sr_size = 4

     strontium_array(1)%z = 38
     strontium_array(1)%a = 84
     strontium_array(1)%symbol ="Sr-84 "
     strontium_array(1)%mass_defect = -80.6495
     strontium_array(1)%atom_percent = 0.56
     strontium_array(1)%zaid = "38084.00c"

     strontium_array(2)%z = 38
     strontium_array(2)%a = 86
     strontium_array(2)%symbol ="Sr-86 "
     strontium_array(2)%mass_defect = -84.5230
     strontium_array(2)%atom_percent = 9.86
     strontium_array(2)%zaid = "38086.00c"

     strontium_array(3)%z = 38
     strontium_array(3)%a = 87
     strontium_array(3)%symbol ="Sr-87 "
     strontium_array(3)%mass_defect = -84.8800
     strontium_array(3)%atom_percent = 7.00
     strontium_array(3)%zaid = "38087.00c"

     strontium_array(4)%z = 38
     strontium_array(4)%a = 88
     strontium_array(4)%symbol ="Sr-88 "
     strontium_array(4)%mass_defect = -87.9216
     strontium_array(4)%atom_percent = 82.58
     strontium_array(4)%zaid = "38088.00c"

     !Calculate isotopic mass for Strontium
     DO i=1,sr_size
          holder = strontium_array(i)%a + (strontium_array(i)%mass_defect/MeV_amu)

          strontium_array(i)%isotopic_mass = holder
     END DO

     !Set mass and holder to zero.
     sr_mass = 0.0
     holder = 0.0

     !Calculate mass of Strontium
     DO i=1, sr_size
          sr_mass = holder + (strontium_array(i)%isotopic_mass* &
          & strontium_array(i)%atom_percent)/100

          holder = sr_mass
     END DO

     !Calculate weight percent of Strontium
     DO i=1, sr_size
          strontium_array(i)%weight_percent = (strontium_array(i)%isotopic_mass &
          &   * strontium_array(i)%atom_percent)/sr_mass
     END DO

     !Number of Yttrium Isotopes
     y_size = 1

     yttrium_array(1)%z = 39
     yttrium_array(1)%a = 89
     yttrium_array(1)%symbol =" Y-89 "
     yttrium_array(1)%mass_defect = -87.7083
     yttrium_array(1)%atom_percent = 100.0
     yttrium_array(1)%zaid = "39089.00c"

     !Calculate isotopic mass for Yttrium
     DO i=1,y_size
          holder = yttrium_array(i)%a + (yttrium_array(i)%mass_defect/MeV_amu)

          yttrium_array(i)%isotopic_mass = holder
     END DO

     !Set mass and holder to zero.
     y_mass = 0.0
     holder = 0.0

     !Calculate mass of Yttrium
     DO i=1, y_size
          y_mass = holder + (yttrium_array(i)%isotopic_mass* &
          & yttrium_array(i)%atom_percent)/100

          holder = y_mass
     END DO

     !Calculate weight percent of Yttrium
     DO i=1, y_size
          yttrium_array(i)%weight_percent = (yttrium_array(i)%isotopic_mass &
          &   * yttrium_array(i)%atom_percent)/y_mass
     END DO

     !Number of Zirconium Isotopes
     zr_size = 5

     zirconium_array(1)%z = 40
     zirconium_array(1)%a = 90
     zirconium_array(1)%symbol ="Zr-90 "
     zirconium_array(1)%mass_defect = -88.7725
     zirconium_array(1)%atom_percent = 51.45
     zirconium_array(1)%zaid = "40090.00c"

     zirconium_array(2)%z = 40
     zirconium_array(2)%a = 91
     zirconium_array(2)%symbol ="Zr-91 "
     zirconium_array(2)%mass_defect = -87.8955
     zirconium_array(2)%atom_percent = 11.22
     zirconium_array(2)%zaid = "40091.00c"

     zirconium_array(3)%z = 40
     zirconium_array(3)%a = 92
     zirconium_array(3)%symbol ="Zr-92 "
     zirconium_array(3)%mass_defect = -88.4490
     zirconium_array(3)%atom_percent = 17.15
     zirconium_array(3)%zaid = "40092.00c"

     zirconium_array(4)%z = 40
     zirconium_array(4)%a = 94
     zirconium_array(4)%symbol ="Zr-94 "
     zirconium_array(4)%mass_defect = -87.2693
     zirconium_array(4)%atom_percent = 17.38
     zirconium_array(4)%zaid = "40094.00c"

     zirconium_array(5)%z = 40
     zirconium_array(5)%a = 96
     zirconium_array(5)%symbol ="Zr-96 "
     zirconium_array(5)%mass_defect = -85.4388
     zirconium_array(5)%atom_percent = 2.80
     zirconium_array(5)%zaid = "40096.00c"

     !Calculate isotopic mass for Zirconium
     DO i=1,zr_size
          holder = zirconium_array(i)%a + (zirconium_array(i)%mass_defect/MeV_amu)

          zirconium_array(i)%isotopic_mass = holder
     END DO

     !Set mass and holder to zero.
     zr_mass = 0.0
     holder = 0.0

     !Calculate mass of Zirconium
     DO i=1, zr_size
          zr_mass = holder + (zirconium_array(i)%isotopic_mass* &
          & zirconium_array(i)%atom_percent)/100

          holder = zr_mass
     END DO

     !Calculate weight percent of Zirconium
     DO i=1, zr_size
          zirconium_array(i)%weight_percent = (zirconium_array(i)%isotopic_mass &
          &   * zirconium_array(i)%atom_percent)/zr_mass
     END DO

   END SUBROUTINE Z36_Z40

max_stars_repo_name: krdepri/MatMCNP

max_issues_repo_name: sandialabs/MatMCNP

max_forks_repo_name: sandialabs/MatMCNP


text:
! Author  : wansooha@gmail.com
! https://github.com/pkgpl/cfgio
module string_conv_mod
  use, intrinsic :: iso_fortran_env

  private
  integer(int32), parameter :: MXNSTR = 256
  integer(int32), parameter :: sp = kind(0.0)
  integer(int32), parameter :: dp = kind(0.d0)

  ! subroutines
  interface from_string
    module procedure from_string_i
    module procedure from_string_f
    module procedure from_string_d
    module procedure from_string_c
    module procedure from_string_z
    module procedure from_string_b
  end interface from_string

  interface to_string
    module procedure to_string_i
    module procedure to_string_f
    module procedure to_string_d
    module procedure to_string_c
    module procedure to_string_z
    module procedure to_string_b
  end interface to_string

  ! functions
  interface tostr
    module procedure i2s
    module procedure f2s
    module procedure d2s
    module procedure c2s
    module procedure z2s
    module procedure b2s
    module procedure s2s
  end interface tostr

  interface tolist
    module procedure to_list_i
    module procedure to_list_f
    module procedure to_list_d
    module procedure to_list_c
    module procedure to_list_z
    module procedure to_list_b
    module procedure to_list_s
  end interface tolist

  public :: from_string, to_string
  public :: tostr, tolist
  public :: list_size, list_size_cmplx, list_size_str
  public :: quote, unquote
  public :: s2i, s2f, s2d, s2c, s2z, s2b

contains

  pure function quote(str) result(v)
    character(len=*), intent(in) :: str
    character(len=:), allocatable :: v
    if(len_trim(str) == 0) then
      v = '""'
      return
    endif
    if(str(1:1) == '"' .or. str(1:1) == "'") then
      v = str
    else
      v = '"'//trim(str)//'"'
    endif
  end function quote

  pure function unquote(str) result(v)
    character(len=*), intent(in) :: str
    character(len=:), allocatable :: v
    integer(int32) :: l
    if(str(1:1) == '"' .or. str(1:1) == "'") then
      l = len_trim(str)
      v = str(2:l - 1)
    else
      v = str
    endif
  end function unquote

  ! string to other types
  pure function s2i(str) result(v)
    character(len=*), intent(in) :: str
    integer(int32) :: v
    read(str, *) v
  end function s2i

  pure function s2f(str) result(v)
    character(len=*), intent(in) :: str
    real(real32) :: v
    read(str, *) v
  end function s2f

  pure function s2d(str) result(v)
    character(len=*), intent(in) :: str
    real(real64) :: v
    read(str, *) v
  end function s2d

  pure function s2c(str) result(v)
    character(len=*), intent(in) :: str
    complex(real32) :: v
    read(str, *) v
  end function s2c

  pure function s2z(str) result(v)
    character(len=*), intent(in) :: str
    complex(real64) :: v
    read(str, *) v
  end function s2z

  pure function s2b(str) result(v)
    character(len=*), intent(in) :: str
    logical :: v
    select case(str)
    case("TRUE", "True", "true", "T", "t")
      v = .true.
    case("YES", "Yes", "yes", "Y", "y")
      v = .true.
    case("ON", "On", "on")
      v = .true.
    case(".TRUE.", ".true.")
      v = .true.
    case default
      v = .false.
    end select
  end function s2b

  subroutine from_string_i(str, v)
    integer(int32), intent(out) :: v
    character(len=*), intent(in) :: str
    read(str, *) v
  end subroutine from_string_i

  subroutine from_string_f(str, v)
    real(real32), intent(out) :: v
    character(len=*), intent(in) :: str
    read(str, *) v
  end subroutine from_string_f

  subroutine from_string_d(str, v)
    real(real64), intent(out) :: v
    character(len=*), intent(in) :: str
    read(str, *) v
  end subroutine from_string_d

  subroutine from_string_c(str, v)
    complex(real32), intent(out) :: v
    character(len=*), intent(in) :: str
    read(str, *) v
  end subroutine from_string_c

  subroutine from_string_z(str, v)
    complex(real64), intent(out) :: v
    character(len=*), intent(in) :: str
    read(str, *) v
  end subroutine from_string_z

  subroutine from_string_b(str, v)
    logical, intent(out) :: v
    character(len=*), intent(in) :: str
    v = s2b(str)
  end subroutine from_string_b

  ! to string
  pure function i2s(v) result(str)
    integer(int32), intent(in) :: v
    character(len=MXNSTR) :: str
    write(str, *) v
  end function i2s

  pure function f2s(v) result(str)
    real(real32), intent(in) :: v
    character(len=MXNSTR) :: str
    write(str, *) v
  end function f2s

  pure function d2s(v) result(str)
    real(real64), intent(in) :: v
    character(len=MXNSTR) :: str
    write(str, *) v
  end function d2s

  pure function c2s(v) result(str)
    complex(real32), intent(in) :: v
    character(len=MXNSTR) :: str
    write(str, *) v
  end function c2s

  pure function z2s(v) result(str)
    complex(real64), intent(in) :: v
    character(len=MXNSTR) :: str
    write(str, *) v
  end function z2s

  pure function b2s(v) result(str)
    logical, intent(in) :: v
    character(len=MXNSTR) :: str
    write(str, *) v
  end function b2s

  pure function s2s(v) result(str)
    character(len=*), intent(in) :: v
    character(len=MXNSTR) :: str
    write(str, *) v
  end function s2s

  subroutine to_string_i(v, str)
    integer(int32), intent(in) :: v
    character(len=*), intent(out) :: str
    write(str, *) v
  end subroutine to_string_i

  subroutine to_string_f(v, str)
    real(real32), intent(in) :: v
    character(len=*), intent(out) :: str
    write(str, *) v
  end subroutine to_string_f

  subroutine to_string_d(v, str)
    real(real64), intent(in) :: v
    character(len=*), intent(out) :: str
    write(str, *) v
  end subroutine to_string_d

  subroutine to_string_c(v, str)
    complex(real32), intent(in) :: v
    character(len=*), intent(out) :: str
    write(str, *) v
  end subroutine to_string_c

  subroutine to_string_z(v, str)
    complex(real64), intent(in) :: v
    character(len=*), intent(out) :: str
    write(str, *) v
  end subroutine to_string_z

  subroutine to_string_b(v, str)
    logical, intent(in) :: v
    character(len=*), intent(out) :: str
    write(str, *) v
  end subroutine to_string_b

  ! array to string
  function to_list_i(arr) result(str)
    integer(int32), intent(in) :: arr(:)
    character(len=MXNSTR) :: str
    character :: sep
    integer(int32) i
    str = ''
    do i = 1, size(arr)
      if(i == 1) then
        sep = ''
      else
        sep = ','
      endif
      str = trim(adjustl(str))//sep//trim(adjustl(tostr(arr(i))))
    enddo
  end function to_list_i

  function to_list_f(arr) result(str)
    real(real32), intent(in) :: arr(:)
    character(len=MXNSTR) :: str
    character :: sep
    integer(int32) i
    str = ''
    do i = 1, size(arr)
      if(i == 1) then
        sep = ''
      else
        sep = ','
      endif
      str = trim(adjustl(str))//sep//trim(adjustl(tostr(arr(i))))
    enddo
  end function to_list_f

  function to_list_d(arr) result(str)
    real(real64), intent(in) :: arr(:)
    character(len=MXNSTR) :: str
    character :: sep
    integer(int32) i
    str = ''
    do i = 1, size(arr)
      if(i == 1) then
        sep = ''
      else
        sep = ','
      endif
      str = trim(adjustl(str))//sep//trim(adjustl(tostr(arr(i))))
    enddo
  end function to_list_d

  function to_list_c(arr) result(str)
    complex(real32), intent(in) :: arr(:)
    character(len=MXNSTR) :: str
    character :: sep
    integer(int32) i
    str = ''
    do i = 1, size(arr)
      if(i == 1) then
        sep = ''
      else
        sep = ','
      endif
      str = trim(adjustl(str))//sep//trim(adjustl(tostr(arr(i))))
    enddo
  end function to_list_c

  function to_list_z(arr) result(str)
    complex(real64), intent(in) :: arr(:)
    character(len=MXNSTR) :: str
    character :: sep
    integer(int32) i
    str = ''
    do i = 1, size(arr)
      if(i == 1) then
        sep = ''
      else
        sep = ','
      endif
      str = trim(adjustl(str))//sep//trim(adjustl(tostr(arr(i))))
    enddo
  end function to_list_z

  function to_list_b(arr) result(str)
    logical, intent(in) :: arr(:)
    character(len=MXNSTR) :: str
    character :: sep
    integer(int32) i
    str = ''
    do i = 1, size(arr)
      if(i == 1) then
        sep = ''
      else
        sep = ','
      endif
      str = trim(adjustl(str))//sep//trim(adjustl(tostr(arr(i))))
    enddo
  end function to_list_b

  function to_list_s(arr) result(str)
    character(len=*), intent(in) :: arr(:)
    character(len=MXNSTR) :: str
    character :: sep
    integer(int32) i
    str = ''
    do i = 1, size(arr)
      if(i == 1) then
        sep = ''
      else
        sep = ','
      endif
      str = trim(adjustl(str))//sep//trim(adjustl(quote(tostr(arr(i))))) !! quote string
    enddo
  end function to_list_s

  integer(int32) function list_size(text) result(npar)
    character(len=*), intent(in) :: text
    character :: sep = ','
    integer(int32) i
    npar = 1
    do i = 1, len(text)
      if(text(i:i) == sep) then
        npar = npar + 1
      endif
    enddo
  end function list_size

  integer(int32) function list_size_str(text) result(npar)
    character(len=*), intent(in) :: text
    character :: sep = ','
    logical :: in_sq = .false., in_bq = .false.
    integer(int32) i
    npar = 1
    do i = 1, len(text)
      if(text(i:i) == "'") in_sq = .not. in_sq
      if(text(i:i) == '"') in_bq = .not. in_bq
      if(.not. in_sq .and. .not. in_bq .and. text(i:i) == sep) then
        npar = npar + 1
      endif
    enddo
  end function list_size_str

  integer(int32) function list_size_cmplx(text) result(npar)
    character(len=*), intent(in) :: text
    character :: sep = ','
    integer(int32) i
    logical :: inside = .false.
    npar = 1
    do i = 1, len(text)
      if(text(i:i) == '(') inside = .true.
      if(text(i:i) == ')') inside = .false.
      if(.not. inside .and. text(i:i) == sep) then
        npar = npar + 1
      endif
    enddo
  end function list_size_cmplx

end module string_conv_mod

max_stars_repo_name: smillerc/lbm

max_issues_repo_name: smillerc/lbm

max_forks_repo_name: smillerc/lbm


text:
Daniel Dan Wolk is a member of the Davis City Council. Dan has served on the City Council since 2011, and voters returned him to the Council in 2012. He was originally appointed by the City Council to fill the unfinished term of Former Mayor Don Saylor who was elected to the Yolo County Board of Supervisors.

Immediately prior to serving on the Davis City Council, Dan was the Chair of the Davis Social Services Commission.

He is a 1995 graduate of Davis Senior High School. He holds a B.A. in economics from Stanford University and a law degree from UC Berkeley Boalt Hall School of Law. He works as Deputy County Counsel for Solano County.

Dan is the son of California State Senator and former Davis Mayor Lois Wolk and former Dean of the UC Davis School of Law Bruce Wolk.

On June 15, 2013, Dan announced his candidacy for the California State Assembly. He would replace Mariko Yamada at the end of her term.

 About Dan
Dan grew up in Davis.  The only city council member in decades to have done so, Dan attended Davis schools (including taking Mr. Willett himself in sixth grade at what is now Willett Elementary School) and graduated from Davis High School in 1995. Dan also participated in many of the activities children in Davis still participate in, including Rainbow Summer, AYSO, Camp Putah, Cub Scouts, Davis High School Madrigal Singers Madrigals, Jazz Choir, and the Nutcracker Davis Children’s Nutcracker.

Dan holds a B.A. in economics with honors and Phi Beta Kappa from Stanford University and a law degree from U.C. Berkeley (Boalt Hall).

After college, Dan worked at the Federal Reserve and then the Public Policy Institute of California as a research assistant. After law school he worked with Davis City Attorney Harriet Steiner in the public law division of McDonough, Holland & Allen.

In 2009, he was appointed by the Davis City Council to the Social Services Commission. He was later elected Chair by his commission colleagues in November, 2010.

Dan Wolk has served on the Davis City Council since his appointment in 2011. Voters returned him to the Council in 2012 with more votes than any other candidate in Davis’ history.

As a City Councilmember, he has led efforts to adopt the surface water project, pushed for greater investment in the city’s roads and bike paths and explored sustainable energy initiatives.

He works as Deputy County Counsel for Solano County, handling public finance, public contracting, and water issues.

Dan is also the founder of the Legal Clinic of Yolo County, a legal services provider for lowincome families.

He lives in Davis with his wife, Jamima, a former professional triathlete and small business owner, and their two young daughters, Avery and Layla.

In his spare time, Dan enjoys cycling, running, camping, hiking, watching the Oakland A’s, and watching ducks at the UC Davis Arboretum with his daughters.

 Davis City Council Campaign

Dan was a candidate in the 2012 City Council Election. He took first pace in the election with 10,212 votes. According to Yolo County Elections, this is the most votes any candidate has ever received in a Davis City Council election. Dan also took first place in every precinct in the city.

 Platform
Davis is my home town; I grew up here, and I want my daughters to enjoy growing up here as much as I did. If I am elected to continue my service on the City Council, I will do everything I can to make a stronger, better Davis for the next generation.

While I have accomplished much during my short tenure, we still face critical issues and decisions today that will define what kind of Davis we will leave to the future. If elected, I will:

Address Longterm Structural Issues in the Budget. We face massive, unfunded liabilities and structural, longterm budgetary challenges. I will work aggressively to ensure our City’s longterm fiscal sustainability in addressing the City’s future obligations, such as pensions and retiree health care. I will also seek ways to restore sustainable revenues by encouraging new and innovative forms of economic development. When cuts are necessary, I will prioritize our future by preserving city funding for children’s programs wherever possible. This is a new era of innovation, compromise, and change. We must be ready for new ideas and shared burdens to make our City sustainable.
Safeguard our Water and Infrastructure. I will explore ways of funding new and existing city infrastructure improvements, and I will work to preserve and maintain our roads, bike lanes, and city green belts. Securing a sustainable water supply is an immediate issue for our city today, but tomorrow it will be a crucial issue for all of California’s communities. I will pursue surface water for our city in a way that is smart, sustainable, and affordable.

Prioritize Economic Development and the Downtown. Davis boasts a worldclass university, an excellent K12 school system, and is centrally located in one of the most fertile agricultural regions in the world. I see in our city the ingredients for a new Silicon Valley of Sustainability, an innovation hub for green technologies. I will make sure that the city stands behind, and not in front of economic development by cutting red tape, streamlining approval processes, furthering technology transfers from UC Davis and keeping those businesses in our community, and appointing a stafflevel champion for business development. Moreover, I will work to make our downtown even better, by creating an Arts and Entertainment District, encouraging new entrepreneurial activity, making streetscape and building improvements, and bettering pedestrian, bicycle, and automobile circulation and parking.

Focus on Environmental Sustainability. I will develop smart, sustainable practices today to ensure a healthy, vibrant Davis tomorrow. I will explore innovative ways of expanding our community’s energy independence and a more efficient use of energy, such as employing a Community Choice Aggregation model, reducing our peak load through weatherization and other upgrades, expanding residential solar power, and pursuing our own public energy utility. I will manage our waste sustainably by partnering with DWR and others to pursue a communitywide waste reduction and composting program. I will follow a policy of smart, limited growth that helps us achieve economic and environmental independence, selfsufficiency and sustainability.

Pursue Affordable Housing for Seniors and Young Families, and Combat Homelessness. The demographic data are clear: our city is aging, and there are fewer young households. While there is no shortage of parents who want to move to Davis because of our safe and outstanding parks and schools, our existing housing stock is unaffordable and often not well suited to new families. No affordable housing for families means less funding for our schools and even school closures. At the same time, many aging householders who wish to retire and remain in Davis want smaller, more accessible housing. If elected, I will spearhead these issues and make attracting and serving young families and developing a “senior housing strategy” a top priority. I will also work to end homelessness in our community by 2020.
Our community faces budget challenges and large unfunded liabilities in infrastructure, pensions, and health care; critical questions about environmental sustainability, such as where we will reliably source our water; and a changing economic environment, where we must find a way to develop a sustainable local economy while providing adequate and affordable housing for young families and seniors.

Only a short while ago, I asked to serve on the Council to make sure that we address these issues now instead of leaving them to our children. While I’ve done a lot in my short tenure to address these issues, there is much left to do. I hope to continue serving, and make a better Davis for all of us. danwolk.org (archived)

 Endorsements

In addition to individual endorsements, Dan Wolk was endorsed by the The Davis Enterprise and Sacramento Bee newspapers.

 The Davis Song

As part of his 2012 City Council Campaign, Dan starred in a music video of The Davis Song. The video can be seen here: http://youtu.be/If1LgiYOmx0

 Campaign Literature

 California State Assembly Campaign
On June 15, 2013, Dan Wolk announced his decision to run for California State Assembly at Sycamore Park, next to Willett Elementary School, which Dan attended as a child. Dan was joined by his family and nearly 100 supporters at the event. He was introduced by his wife, Jamima, and talked about the top priorities of his campaign.

In his announcement speech, Dan said, “I am running for Assembly to restore the California dream for our children and grandchildren. Our state has been neglecting their future and that is unacceptable. In the Assembly, I will work to ensure we fully invest in public education, our infrastructure, our environment and our middle class.”

Assembly District 4 consists of Yolo, Napa and Lake Counties, as well as portions of Sonoma, Solano and Colusa Counties. It is currently represented by Assemblywoman Mariko Yamada.

 Endorsements
Dan has been endorsed by more than 90 community, education, public safety and other leaders from all over the district. The list includes all five Davis Board of Education School Board members and both members of the Yolo County Board of Supervisors who represent Davis.

State Elected Officials
Jim Beall, State Senator
Kevin de León, State Senator
Mark DeSaulnier, State Senator
Loni Hancock, State Senator
Lois Wolk, State Senator
Jim Frazier, Assemblymember
Delaine Eastin, Former State Superintendent of Public Instruction
Fabian Núñez, Former State Assembly Speaker
Tom Hannigan, Former Assemblymember and Director of the California Department of Water Resources
Helen Thomson, Former State Assemblymember and Yolo County Supervisor

County Elected Officials
Tom Ferrara, Solano County Sheriff
Chuck Lomeli, Solano County TreasurerTax CollectorCounty Clerk
Marc Tonnesen, Solano County AssessorRecorder

David Sundstrom, Sonoma County AuditorController, Treasurer Tax Collector

Ed Prieto, Yolo County Sheriff
Jeff Reisig, Yolo County District Attorney
Cass Sylvia, Yolo County Public Guardian
Howard Newens, Yolo County AuditorController, Treasurer and Tax Collector

Linda Seifert, Solano County Board of Supervisors
Skip Thomson, Solano County Board of Supervisors
Mike McGowan, Vice Chair, Yolo County Board of Supervisors
Jim Provenza, Yolo County Board of Supervisors
Don Saylor, Yolo County Board of Supervisors
Duane Kromm, Former Solano Supervisor

City Elected Officials
Jack Batchelor, Jr., Mayor, City of Dixon
Harry Price, Mayor, City of Fairfield
Steve Hardy, Mayor, City of Vacaville
Dilenna Harris,Vice Mayor, City of Vacaville
Pete Mott, Vice Mayor, City of Napa
Jerry Adler, Former Davis Mayor
Ann M. Evans, Former Davis Mayor
Debbie Nichols Poulos, Former Davis Mayor
Ted Puntillo, Former Davis Mayor
Maynard Skinner, Former Davis Mayor
Gary Sandy, Former Woodland Mayor, Yuba Community College District Trustee

Denise Loustalot, Clearlake City Council Member
Joey Luiz, Clearlake City Council Member
Jerry Castanon, Dixon City Council Member
Mike Healy, Petaluma City Council Member
Mitch Mashburn, Vacaville City Council Member
Sean Denny, Woodland City Council Member
Jim Hilliard, Woodland City Council Member

School Boards and Education Leaders
Wally Holbrook, Lake County Superintendent of Schools
Jay Speck, Solano County Superintendent of Schools
Dr. Jorge Ayala, Yolo County Superintendent of Schools

Shelton Yip, Yolo County Board of Education Member
Dana Dean, Solano County Board of Education Member
Rosemary Thurston, Solano Community College District Governing Board

Sheila Allen, Davis School Board
Gina Daldeiden, Davis School Board
Susan Lovenburg, Davis School Board
Nancy Peterson, Davis School Board
Tim Taylor, Davis School Board
Herb Cross, Dixon School Board
Caitlin O’Halloran, Dixon School Board MemberElect
Judi Honeychurch, FairfieldSuisun School Board
Michael Kitzes, Vacaville School Board
Marty West, Former Davis School Board

Julie Gallelo, Executive Director, First 5 Yolo Children & Families Commissioner

Davis Campbell, Former Executive Director, California School Boards Association

Community Leaders
Mark Braly, Davis Planning Commissioner
Herman Boschken, Davis Planning Commissioner
Ananya Choudhuri, Davis Planning Commissioner
George Hague, Davis Planning Commissioner
Dawn LaBar, Vacaville Planning Commissioner

Yvonne Clinton, Davis Human Relations Commissioner
Tilahun Yilma, Davis Human Relations Commissioner

Judy Wolf, Chair, Davis Social Services Commissioner
Matthew Wise, Davis Social Services Commissioner

Laura Westrup, Chair, Davis Tree Commissioner

Gene Wilson, Chair, Davis Natural Resources Commissioner
Ben Bourne, Davis Natural Resources Commissioner

Steve Boschken, Chair, Davis Business and Economic Development Commissioner
Brian Horsfield, Davis Business and Economic Development Commissioner
Jag Nagendra, Davis Business and Economic Development Commissioner
David Robert, Davis Business and Economic Development Commissioner

Michael Bartolic, Davis Water Advisory Committee Member
Petrea Marchand, Davis Water Advisory Committee Member
Matt Williams, Jr., Davis Water Advisory Committee Member

Beth Bourne, Davis Bicycle Advisory Commissioner

Bill Bernheim, Chair, Davis Safety and Parking Advisory Commissioner
Jim Watson, Davis Safety and Parking Advisory Commissioner

Dr. Geni Bennetts, M.D.
Steve Perry, Napa County Educator
Cindy Water, Napa County Educator
Sandee Dickson, Retired Napa County Educator
Tobin Booth, Founder and CEO, Blue Oak Energy
Dick Bourne, Davis Energy Group
Russ Lester, Dixon Ridge Farms
Dorothy Peterson, Founder, Davis Farm to School

Organizations
Davis Police Officers Association

 Campaign Literature

20110301 20:43:44 nbsp He was officially sworn in tonight and it was actually a nice, albeit short, ceremony. Congratulations to newly sworn Councilman Wolk. Users/WesP

20131027 21:18:03 nbsp Thumbs up for the Davis YouTube clip http://youtu.be/If1LgiYOmx0 The Davis Song (with city council candidate Dan Wolk)! Do you have this in HD too? IMO You could get more out of it if you do. So many Davis images we all recognize and then this great end with a variation on http://youtu.be/4kjmIXr6MoI Empire State of Mind! Users/ConstantiaOomen

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!----------------------Function----------------------------------------!
!
      FUNCTION FNHGBL( NBS,NBE,M )
!
!-------------------------Disclaimer-----------------------------------!
!
!     This material was prepared as an account of work sponsored by
!     an agency of the United States Government. Neither the
!     United States Government nor the United States Department of
!     Energy, nor Battelle, nor any of their employees, makes any
!     warranty, express or implied, or assumes any legal liability or
!     responsibility for the accuracy, completeness, or usefulness
!     of any information, apparatus, product, software or process
!     disclosed, or represents that its use would not infringe
!     privately owned rights.
!
!----------------------Acknowledgement---------------------------------!
!
!     This software and its documentation were produced with Government
!     support under Contract Number DE-AC06-76RLO-1830 awarded by the
!     United Department of Energy. The Government retains a paid-up
!     non-exclusive, irrevocable worldwide license to reproduce,
!     prepare derivative works, perform publicly and display publicly
!     by or for the Government, including the right to distribute to
!     other Government contractors.
!
!---------------------Copyright Notices--------------------------------!
!
!            Copyright Battelle Memorial Institute, 1996
!                    All Rights Reserved.
!
!----------------------Description-------------------------------------!
!
!     Aqueous hydraulic gradient boundary condition function.
!
!----------------------Authors-----------------------------------------!
!
!     Written by MD White, Battelle, PNL, June 1994.
!     Last Modified by MD White, Battelle, PNL, September 10, 1997.
!     Last Modified by MD White, Battelle, PNNL, July 12, 2000.
!     $Id: fnhgbl.F 1080 2017-03-14 16:22:02Z d3c002 $
!

!----------------------Fortran 90 Modules------------------------------!
!
      USE GLB_PAR
      USE SOLTN
      USE GRID
      USE CONST
      USE BCVP
      USE BCV
!

!----------------------Implicit Double Precision-----------------------!
!
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
!
!----------------------Parameter Statements----------------------------!
!

!
!----------------------Common Blocks-----------------------------------!
!

!
!----------------------Executable Lines--------------------------------!
!
      ISUB_LOG = ISUB_LOG+1
      SUB_LOG(ISUB_LOG) = '/FNHGBL'
      IF( INDEX(SVN_ID(57)(1:1),'$').EQ.0 ) SVN_ID(57) =
     & '$Id: fnhgbl.F 1080 2017-03-14 16:22:02Z d3c002 $'
!
!---  Set starting-point pressures and temperatures  ---
!
      PLX = PLB(M,NBS)
      PGX = PGB(M,NBS)
      PNX = PLB(M,NBS)
      IF( IEQO.GT.0 ) PNX = PNB(M,NBS)
      TX = TB(M,NBS)
      PX = MAX( PLX,PGX,PNX ) + PATM
      CALL WATLQD( TX,PX,RHOLX )
      NS = IBCN(NBS)
      NE = IBCN(NBE)
      ISX = ID(NS)
      IEX = ID(NE)
      JSX = JD(NS)
      JEX = JD(NE)
      KSX = KD(NS)
      KEX = KD(NE)
      II = SIGN(1,IEX-ISX)
      JI = SIGN(1,JEX-JSX)
      KI = SIGN(1,KEX-KSX)
!
!---  Convert boundary pressure to node-centroid pressure  ---
!
      IF( IBCD(NBS).EQ.-3 ) THEN
        NPZ = NSZ(NS)
        PLX = PLX - 5.D-1*DZGF(NS)*RHOLX*GRVZ(NPZ)
      ELSEIF( IBCD(NBS).EQ.-2 ) THEN
        NPY = NSY(NS)
        PLX = PLX - 5.D-1*DYGF(NS)*RHOLX*GRVY(NPY)*RP(ID(NS))
      ELSEIF( IBCD(NBS).EQ.-1 ) THEN
        NPX = NSX(NS)
        PLX = PLX - 5.D-1*DXGF(NS)*RHOLX*GRVX(NPX)
      ELSEIF( IBCD(NBS).EQ.1 ) THEN
        NQX = NSX(NS)+1
        PLX = PLX + 5.D-1*DXGF(NS)*RHOLX*GRVX(NQX)
      ELSEIF( IBCD(NBS).EQ.2 ) THEN
        NQY = NSY(NS)+IFLD
        PLX = PLX + 5.D-1*DYGF(NS)*RHOLX*GRVY(NQY)*RP(ID(NS))
      ELSEIF( IBCD(NBS).EQ.3 ) THEN
        NQZ = NSZ(NS)+IJFLD
        PLX = PLX + 5.D-1*DZGF(NS)*RHOLX*GRVZ(NQZ)
      ENDIF
!
!---  Loop over nodes in the x direction
!
      I = ISX
      J = JSX
      K = KSX
      NX = ND(I,J,K)
      DO 100 I = ISX+II,IEX,II
        N = ND(I,J,K)
        NPX = NSX(N)
        GB = (XP(N)-XP(NX))*GRVX(NPX)
        PX = MAX( PLX,PGX,PNX ) + PATM
        CALL WATLQD( TX,PX,RHOLX )
        PLX = PLX - RHOLX*GB
        NX = N
  100 CONTINUE
!
!---  Loop over nodes in the y direction
!
      I = IEX
      J = JSX
      K = KSX
      NX = ND(I,J,K)
      DO 110 J = JSX+JI,JEX,JI
        N = ND(I,J,K)
        NPY = NSY(N)
        GB = (YP(N)*RP(I)-YP(NX)*RP(I))*GRVY(NPY)
        PX = MAX( PLX,PGX,PNX ) + PATM
        CALL WATLQD( TX,PX,RHOLX )
        PLX = PLX - RHOLX*GB
        NX = N
  110 CONTINUE
!
!---  Loop over nodes in the z direction
!
      I = IEX
      J = JEX
      K = KSX
      NX = ND(I,J,K)
      DO 120 K = KSX+KI,KEX,KI
        N = ND(I,J,K)
        NPZ = NSZ(N)
        GB = (ZP(N)-ZP(NX))*GRVZ(NPZ)
        PX = MAX( PLX,PGX,PNX ) + PATM
        CALL WATLQD( TX,PX,RHOLX )
        PLX = PLX - RHOLX*GB
        NX = N
  120 CONTINUE
!
!---  Convert node-centroid pressure to boundary pressure  ---
!
      NE = IBCN(NBE)
      IF( IBCD(NBE).EQ.-3 ) THEN
        NPZ = NSZ(NE)
        PLX = PLX + 5.D-1*DZGF(NE)*RHOLX*GRVZ(NPZ)
      ELSEIF( IBCD(NBE).EQ.-2 ) THEN
        NPY = NSY(NE)
        PLX = PLX + 5.D-1*DYGF(NE)*RHOLX*GRVY(NPY)*RP(ID(NE))
      ELSEIF( IBCD(NBE).EQ.-1 ) THEN
        NPX = NSX(NE)
        PLX = PLX + 5.D-1*DXGF(NE)*RHOLX*GRVX(NPX)
      ELSEIF( IBCD(NBE).EQ.1 ) THEN
        NQX = NSX(NE)+1
        PLX = PLX - 5.D-1*DXGF(NE)*RHOLX*GRVX(NQX)
      ELSEIF( IBCD(NBE).EQ.2 ) THEN
        NQY = NSY(NE)+IFLD
        PLX = PLX - 5.D-1*DYGF(NE)*RHOLX*GRVY(NQY)*RP(ID(NE))
      ELSEIF( IBCD(NBE).EQ.3 ) THEN
        NQZ = NSZ(NE)+IJFLD
        PLX = PLX - 5.D-1*DZGF(NE)*RHOLX*GRVZ(NQZ)
      ENDIF
      FNHGBL = PLX
      ISUB_LOG = ISUB_LOG-1
!
!---  End of FNHGBL group
!
      RETURN
      END

max_stars_repo_name: pnnl/stomp-w

max_issues_repo_name: pnnl/stomp-w

max_forks_repo_name: pnnl/stomp-w


text:
      subroutine wconstruct(etvec,plept,lepton,etmissing)
      implicit none
      include 'constants.f'
      include 'masses.f'
      double precision etvec(2),plept(4),pw(4),etmissing(4),
     . pte,ptn,ptwsq,a,b,discr,plnp,plnm,Ee,ple,mtsq,Enp,Enm
      double precision yrap1,yrap2
      integer j
      character*2 lepton

      do j=1,2
        pw(j)=plept(j)+etvec(j)
        etmissing(j)=etvec(j)
      enddo

      pte=dsqrt(plept(1)**2+plept(2)**2)
      ptn=dsqrt(etvec(1)**2+etvec(2)**2)
      ptwsq=(pw(1)**2+pw(2)**2)
      Ee=plept(4)
      ple=plept(3)
      mtsq=(pte+ptn)**2-ptwsq

c  + (Mw^2 - MT^2) * Ee^2 * ( 4*pte*ptn + Mw^2 - Mt^2 )
c
c - pln
c  * ple*( 2*pte*ptn + Mw^2 - MT^2 )
c
c + pln^2
c  * ( pte^2 ) + 0.

c      write(6,*) 'wmass**2,mtsq',wmass**2,mtsq
      if (wmass**2 .gt. mtsq) then
        discr=Ee*dsqrt((wmass**2-mtsq)*(wmass**2-mtsq+4d0*pte*ptn))
        a=pte**2
        b=-(2d0*pte*ptn+wmass**2-mtsq)*ple
c      c=-0.25d0*(wmass**2-mtsq+2*pte*ptn)**2
c     # +ptn**2*Ee**2
c      write(6,*) 'discr',sqrt(b**2-4*a*c),discr
      else
        discr=0d0
        a=pte**2
        b=-2*pte*ptn*ple
      endif

c      wsq=pw(4)**2-pw(1)**2-pw(2)**2-pw(3)**2
c      write(6,*) 'cal:mw',wsq,plntrue
      plnp=(-b+discr)/(two*a)
      plnm=(-b-discr)/(two*a)

      Enm=dsqrt(plnm**2+ptn**2)
      Enp=dsqrt(plnp**2+ptn**2)
      yrap1=(Enm+plnm)/(Enm-plnm)
      yrap2=(Enp+plnp)/(Enp-plnp)
      if (yrap1 .lt. 1d-13) then
        yrap1=100d0
      else
        yrap1=0.5d0*dlog(yrap1)
      endif
      if (yrap2 .lt. 1d-13) then
        yrap2=100d0
      else
        yrap2=0.5d0*dlog(yrap2)
      endif

      if     ((lepton.eq.'ea') .or. (lepton.eq.'ma')) then
c--- choose the solution with the larger rapidity
        if (yrap1 .gt. yrap2) then
          etmissing(4)=Enm
          etmissing(3)=plnm
        else
          etmissing(4)=Enp
          etmissing(3)=plnp
        endif
      elseif ((lepton.eq.'el') .or. (lepton.eq.'ml')) then
c--- choose the solution with the smaller rapidity
        if (yrap1 .lt. yrap2) then
          etmissing(4)=Enm
          etmissing(3)=plnm
        else
          etmissing(4)=Enp
          etmissing(3)=plnp
        endif
      else
        write(6,*) 'Problem: not a lepton in wconstruct!'
        stop
      endif

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
c READIM -- Read through an image and count the lines (used for timing tests).
c Tests line sequential i/o.
c
c	usage:  readim image
c ----------------------------------------------------------------------------

	program readim

	character*80	image, errmsg
	integer		ncols, nlines, nbands, j, k
	integer		im, ier, axlen(7), naxis, dtype
	integer*2	pix(8192)

c --- Get image name.
	call clargc (1, image, ier)
	if (ier .ne. 0) then
	    write (*, '('' enter image name: '',$)')
	    read (*,*) image
	endif

c --- Open the image.
	call imopen (image, 1, im, ier)
	if (ier .ne. 0) goto 91
	call imgsiz (im, axlen, naxis, dtype, ier)
	if (ier .ne. 0) goto 91

	ncols  = axlen(1)
	nlines = axlen(2)
	nbands = axlen(3)

c --- Read through the image.
	do 20 k = 1, nbands
	    do 10 j = 1, nlines
		call imgl3s (im, pix, j, k, ier)
		if (ier .ne. 0) goto 91
 10	    continue
 20	continue

c --- Clean up.
	call imclos (im, ier)
	if (ier .ne. 0) goto 91

	print 81, nlines, image
 81	format (' read ', i4, ' lines from image ', a64)

	stop

c --- Error exit.
 91	call imemsg (ier, errmsg)
	write (*, '('' Error: '', a80)') errmsg

	stop
	end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
\ utils.f      - debug utility functions
\ ------------------------------------------------------------------------

  .( butils.f )

\ ------------------------------------------------------------------------

  0 var app-ip              \ address where debugger is debugging
  0 var app-rp              \ applications return stack pointer
  0 var app-sp              \ applications parameter stack pointer
  0 var memaddr             \ memory window address
  0 var bug-stacks          \ points to stacks buffer for debugger
  0 var app-rp0             \ apps rp address on entry into debugger
  0 var bug-base            \ system radix on entry into debugger
  0 var mem-addr            \ memory window address
  0 var bug-stacks          \ pointer to allocated stacks for debugger
  0 var app-rp0             \ applications rp on entry into debugger
  0 var halted              \ true if single stepping is halted
  0 var stepping            \ true if single stepping (auto)
  0 var updating            \ true if updating during auto step
  0 var stepto              \ address to auto step to
  0 var break0
$c8 var step-delay          \ auto step delay between steps
  0 var app-out             \ true if application output window is visible

\ ------------------------------------------------------------------------

  0 var []xu                \ execution unit array
  0 var #xu                 \ total number of items in array
  0 var csr-ix              \ current cursor index within xu array
  0 var csr-line            \ which line of decompilation is sursor on
  0 var ip-index            \ current ip index within []xu array
  0 var mid-point           \ mid line of code window
  0 var in-xu               \ spaces inside single xu to be highlighted

\ an execution unit is either a single execution token or an execution
\ token and its operands (for example (.") and the following string are
\ what im calling a single execution unit.

\ ------------------------------------------------------------------------

  20 stack: see-stack
  20 stack: step-stack

\ ------------------------------------------------------------------------
\ flush all entries from execution unit array

: clear-[]xu    ( --- )
  []xu 4096 erase           \ clear xu array
  off> #xu ;

\ ------------------------------------------------------------------------
\ add current decompilation address to array of xu addresses

: +xu           ( a1 --- )
  []xu #xu []!
  incr> #xu ;

\ ------------------------------------------------------------------------

: app-ip@       ( --- xt ) app-ip @ ;
: app-ip++      ( --- )    cell +!> app-ip ;

\ ------------------------------------------------------------------------
\ allocate array of code execution unit addresses

: alloc-xu     ( --- )
  4096 allocate drop
  !> []xu ;                 \ 4k of xu's should be enough for 1 definition

\ ------------------------------------------------------------------------
\ deallocate xu array

: free-xu       ( --- )
  []xu free drop ;

\ ------------------------------------------------------------------------
\ does cfa reference one of the following?

: (?is-:)       ( cfa --- f1 )  ?cfa      ['] nest = ;
: (?is-defer)   ( cfa --- f1 )  ?cfa      ['] dodefer = ;
: (?is-does)    ( cfa --- f1 )  ?cfa ?cfa ['] dodoes = ;

\ -------------------------------------------------------------------------
\ what is the next xt to be stepped by application?

: ?is-:         ( --- f1 )  app-ip@ (?is-:) ;
: ?is-rep       ( --- f1 )  app-ip@ ['] dorep = ;

\ -------------------------------------------------------------------------

: isbreak       ( a1 --- f1 ) break0 = ;

\ -------------------------------------------------------------------------

: app-rp@       ( --- n1 ) app-rp @ ;
: >app-rp       ( n1 --- ) [ cell negate ]# +!> app-rp app-rp ! ;
: app-rp>       ( --- n1 ) app-rp@ cell +!> app-rp ;

: app-sp@       ( --- n1 ) app-sp @ ;
: >app-sp       ( n1 --- ) [ cell negate ]# +!> app-sp app-sp ! ;
: app-sp>       ( --- n1 ) app-sp@ cell +!> app-sp ;

\ =========================================================================


text:
subroutine FOURSAIL(rho, tau,             & !! Leaf spectra
        LIDFa,LIDFb,TypeLIDF, & !! LIDF parameters and function
        lai,q,                & !! LAI and hot spot
        tts,tto,psi,          & !! Instrument geometry
        rsoil,                & !! Soil reflectance

        !! Outputs
        rddt,                 & !! Bi-hemispherical reflectance
        rsdt,                 & !! Directional-hemispherical reflectance for solar incident flux
        rdot,                 & !! Hemispherical-directional reflectance in viewing direction
        rsot)                   !! Bi-directional reflectance factor

        use mod_dataSpec_wavelength
        implicit none
        real*8,intent(in),dimension(nw) :: rho,tau,rsoil
        REAL*8,INTENT(in) :: LIDFa,LIDFb,lai,q,tts,tto,psi
        INTEGER,INTENT(in) :: TypeLIDF

        common /angle/ pi,rd
        real*8 :: pi,rd

        real*8,intent(out),dimension(nw) :: rddt,rsdt,rdot,rsot

        INTEGER*4 :: na
        parameter(na = 13)
        real*8 :: cts,cto,ctscto,dso,lidf(na),litab(na)
        real*8 :: ks,ko,sob,sof,sdb,sdf,dob,dof,ddb,ddf,tss,too
        real*8 :: tsstoo,sumint
        real*8,dimension(nw) :: sigb,att,m,sb,sf,vb,vf,w
        real*8,dimension(nw) :: rdd,tdd,tsd,rsd,tdo,rdo,rsod

        data litab/5.,15.,25.,35.,45.,55.,65.,75.,81.,83.,85.,87.,89./

        pi = 4 * atan(1.)
        rd = pi / 180.

        !! FLAG 1 - Geometry: Compute geometric quantities
        CALL sail_sensgeom(tts,tto,psi, &
            cts,cto,ctscto,dso)

        !! FLAG 2 - Leaf angle distribution function
        CALL LIDF_fun(TypeLIDF, na, LIDFa, LIDFb, &
            lidf)

        !! FLAG 3 (if 1 or 2): Angular distance -- compensation of shadow length
        CALL SUITS(na,litab,lidf,tts,tto,cts,cto,psi,ctscto, &
            ks,ko,sob,sof,sdb,sdf,dob,dof,ddb,ddf)

        !! FLAG 4 (if 3 or leaf refl) -- Geometric leaf reflectance
        !!      Input rho and tau (from PROSPECT)
        !!      Correct using SUITS factors
        CALL RTgeom(rho,tau,ddb,ddf,sdb,sdf,dob,dof,sob,sof, &
            sigb,att,m,sb,sf,vb,vf,w)

        !! FLAG 5 (if 4 or LAI or soil) -- Canopy reflectance and transmittance
        !!      Input LAI and leaf rho and tau
        !!      Also give geometric factors
        CALL ReflTrans(rho,tau,lai,att,m,sigb,ks,ko,sf,sb,vf,vb, &
            rdd,tdd,tsd,rsd,tdo,rdo,tss,too,rsod)

        !! FLAG 6 (if 3 or LAI or hot) -- Hot spot effect
        CALL HotSpot(lai,q,tss,ks,ko,dso, &
            tsstoo,sumint)

        !! Bidirectional reflectance
        CALL sail_BDRF(w,lai,sumint,tsstoo,rsoil, &
                rdd,tdd,tsd,rsd,tdo,rdo,tss,too,rsod, &
                rddt,rsdt,rdot,rsot)

        return
    end subroutine FOURSAIL

max_stars_repo_name: ankurdesai/pecan

max_stars_repo_licenses: NCSA Unlicense

max_issues_repo_name: ankurdesai/pecan

max_issues_repo_licenses: NCSA Unlicense

max_forks_repo_name: ankurdesai/pecan

max_forks_repo_licenses: NCSA Unlicense


text:
mangerson can describe himself or herself here.

20100718 23:09:04 nbsp Hello and welcome to the Wiki!  Users/hankim


module WaterAbsorptionClass
    use fem
    use DiffusionEquationClass
    use FiniteDeformationClass
    implicit none

    type :: WaterAbsorption_
        type(FEMDomain_),pointer:: Water, Tissue
        type(MaterialProp_),pointer:: a_Psi
        type(MaterialProp_),pointer:: a_P
        type(MaterialProp_),pointer:: theta_eq
        type(MaterialProp_),pointer:: Psi_eq
        type(MaterialProp_),pointer:: a_E
        type(MaterialProp_),pointer:: a_v
        type(MaterialProp_),pointer:: E_eq
        type(MaterialProp_),pointer:: v_eq

        type(DiffusionEq_)::DiffusionEq
        type(FiniteDeform_)::FiniteDeform

        real(real64),allocatable :: WaterAbsorbingPower(:)
        real(real64),allocatable :: WaterPotential(:)
        real(real64),allocatable :: TurgorPressure(:)
        real(real64),allocatable :: WaterContent(:)
        real(real64),allocatable :: WaterMass(:)
        real(real64),allocatable :: volume(:)
        real(real64),allocatable :: Conductivity(:)
        real(real64),allocatable :: YoungsModulus(:)
        real(real64),allocatable :: PoissonsRatio(:)
        real(real64),allocatable :: Permiability(:)
        real(real64),allocatable :: PorePressure(:)
        ! Parameters(A,B)
        ! A : Element ID

        real(real64),allocatable :: a_Psi_val(:)
        real(real64),allocatable :: a_P_val(:)
        real(real64),allocatable :: theta_eq_val(:)
        real(real64),allocatable :: Psi_eq_val(:)
        real(real64),allocatable :: a_E_val(:)
        real(real64),allocatable :: a_v_val(:)
        real(real64),allocatable :: E_eq_val(:)
        real(real64),allocatable :: v_eq_val(:)

        real(real64),allocatable :: theta_ps_val(:) ! computed from other variables

        character(200) :: Name= " "
        integer(int32) :: timestep=0
        integer(int32) :: flushstep=0
        integer(int32) :: ID=0
        real(real64) :: dt
    contains
        procedure, public :: import=> importWaterAbsorption
        procedure, public :: assign=> importWaterAbsorption
        procedure, public :: init=> initWaterAbsorption
        procedure, public :: run => runWaterAbsorption
        procedure, public :: update=> updateWaterAbsorption
        procedure, public :: export=> exportWaterAbsorption
        procedure, public :: display => displayWaterAbsorption
        procedure, public :: bake => bakeWaterAbsorption
        procedure, public :: gnuplot => gnuplotWaterAbsorption
        procedure, public :: updatePermiability => updatePermiabilityWA
        procedure, public :: updateStiffness => updateStiffnessWA
        procedure, public :: UpdateVolume => UpdateVolumeWA

        ! <<<<<<<< IO >>>>>>>>>>
        procedure, public :: save =>saveWaterAbsorption
        procedure, public :: open => openWaterAbsorption
        procedure, public :: link => linkWaterAbsorption
        procedure, public :: result => resultWaterAbsorption
        procedure, public :: remove => removeWaterAbsorption
        ! <<<<<<<< IO >>>>>>>>>>

    end type
contains

!#####################################
subroutine removeWaterAbsorption(obj)
    class(WaterAbsorption_),intent(inout) :: obj

    if(associated(obj%Water)) nullify(obj%Water)
    if(associated(obj%Tissue)) nullify(obj%Tissue)

    if(associated(obj%a_Psi))  nullify(obj%a_Psi)
    if(associated(obj%a_P))  nullify(obj%a_P)
    if(associated(obj%theta_eq))  nullify(obj%theta_eq)
    if(associated(obj%Psi_eq))  nullify(obj%Psi_eq)
    if(associated(obj%a_E))  nullify(obj%a_E)
    if(associated(obj%a_v))  nullify(obj%a_v)
    if(associated(obj%E_eq))  nullify(obj%E_eq)
    if(associated(obj%v_eq))  nullify(obj%v_eq)

    call obj%DiffusionEq%remove()
    call obj%FiniteDeform%remove()

    if(allocated(obj%WaterAbsorbingPower)) deallocate(obj%WaterAbsorbingPower)
    if(allocated(obj%WaterPotential)) deallocate(obj%WaterPotential)
    if(allocated(obj%TurgorPressure)) deallocate(obj%TurgorPressure)
    if(allocated(obj%WaterContent)) deallocate(obj%WaterContent)
    if(allocated(obj%WaterMass)) deallocate(obj%WaterMass)
    if(allocated(obj%Conductivity)) deallocate(obj%Conductivity)
    if(allocated(obj%YoungsModulus)) deallocate(obj%YoungsModulus)
    if(allocated(obj%PoissonsRatio)) deallocate(obj%PoissonsRatio)
    if(allocated(obj%Permiability)) deallocate(obj%Permiability)
    if(allocated(obj%PorePressure)) deallocate(obj%PorePressure)
    ! Parameters(A,B)
    ! A : Element ID

    if(allocated(obj%a_Psi_val)) deallocate(obj%a_Psi_val)
    if(allocated(obj%a_P_val)) deallocate(obj%a_P_val)
    if(allocated(obj%theta_eq_val)) deallocate(obj%theta_eq_val)
    if(allocated(obj%Psi_eq_val)) deallocate(obj%Psi_eq_val)
    if(allocated(obj%a_E_val)) deallocate(obj%a_E_val)
    if(allocated(obj%a_v_val)) deallocate(obj%a_v_val)
    if(allocated(obj%E_eq_val)) deallocate(obj%E_eq_val)
    if(allocated(obj%v_eq_val)) deallocate(obj%v_eq_val)
    if(allocated(obj%theta_ps_val)) deallocate(obj%theta_ps_val)

    obj%Name=" "
    obj%timestep=0
    obj%flushstep=0
    obj%dt=1.0d0

end subroutine
!#####################################

!#####################################
subroutine resultWaterAbsorption(obj,path,name,step)
    class(WaterAbsorption_),intent(inout) :: obj
    character(*),intent(in) :: path, name
    integer(int32),intent(in) :: step

    call obj%DiffusionEq%display(OptionalStep=step,Name=trim(path)//"/output/"//trim(adjustl(name)))
    call obj%FiniteDeform%display(OptionalStep=step,Name=trim(path)//"/output/"//trim(adjustl(name)))

end subroutine
!#####################################

!#####################################
subroutine linkWaterAbsorption(obj,Water, Tissue,a_Psi,a_P,theta_eq,Psi_eq,a_E,a_v,E_eq,v_eq)
    class(WaterAbsorption_),intent(inout) :: obj
    type(FEMDomain_),optional,target :: Water, Tissue
    type(MaterialProp_),optional,target :: a_Psi,a_P,theta_eq,Psi_eq,a_E,a_v,E_eq,v_eq

    if(present(Water) )then
        if(associated(obj%Water) ) nullify(obj%Water)
        obj%Water => Water
    endif
    if(present(Tissue) )then
        if(associated(obj%Tissue) ) nullify(obj%Tissue)
        obj%Tissue => Tissue
    endif
    if(present(a_Psi) )then
        if(associated(obj%a_Psi) ) nullify(obj%a_Psi)
        obj%a_Psi => a_Psi
    endif
    if(present(a_P) )then
        if(associated(obj%a_P) ) nullify(obj%a_P)
        obj%a_P => a_P
    endif
    if(present(theta_eq) )then
        if(associated(obj%theta_eq) ) nullify(obj%theta_eq)
        obj%theta_eq => theta_eq
    endif
    if(present(Psi_eq) )then
        if(associated(obj%Psi_eq) ) nullify(obj%Psi_eq)
        obj%Psi_eq => Psi_eq
    endif
    if(present(a_E) )then
        if(associated(obj%a_E) ) nullify(obj%a_E)
        obj%a_E => a_E
    endif
    if(present(a_v) )then
        if(associated(obj%a_v) ) nullify(obj%a_v)
        obj%a_v => a_v
    endif
    if(present(E_eq) )then
        if(associated(obj%E_eq) ) nullify(obj%E_eq)
        obj%E_eq => E_eq
    endif
    if(present(v_eq) )then
        if(associated(obj%v_eq) ) nullify(obj%v_eq)
        obj%v_eq => v_eq
    endif

    if(present(Water) )then
        if(associated(obj%DiffusionEq%FEMDomain) ) nullify(obj%DiffusionEq%FEMDomain)
        obj%DiffusionEq%FEMDomain => Water
    endif

    if(present(Tissue) )then
        if(associated(obj%FiniteDeform%FEMDomain) ) nullify(obj%FiniteDeform%FEMDomain)
        obj%FiniteDeform%FEMDomain => Tissue
    endif
end subroutine
!#####################################

!#####################################
subroutine openWaterAbsorption(obj,path,name)
    class(WaterAbsorption_),intent(inout) :: obj
    character(*),optional,intent(in)::path
    character(*),optional,intent(in) :: Name
    character(200) ::fname
    integer(int32) :: fh
    logical :: restart=.true.
    type(IO_) :: f

    if(.not. present(path) )then
        print *, " exportWaterAbsorption ERROR >> .not. present(path)"
        stop
    endif

    if(present(name) )then
        call system("mkdir -p "//trim(path)//"/"//trim(adjustl(name)) )
	    fname=trim(path)//"/"//trim(adjustl(name))
        call f%open("./",trim(path)//"/"//trim(adjustl(name)) ,"/WaterAbsorption.prop")
        call obj%FiniteDeform%open(path=trim(path)//"/"//trim(adjustl(name)),name="FiniteDeform")
        call obj%DiffusionEq%open(path=trim(path)//"/"//trim(adjustl(name)),name="DiffusionEq")
    else
        call system("mkdir -p "//trim(path)//"/WaterAbsorption")
        call f%open("./",trim(path)//"/WaterAbsorption","/WaterAbsorption.prop")
        call obj%FiniteDeform%open(path=trim(path)//"/WaterAbsorption",name="FiniteDeform")
        call obj%DiffusionEq%open(path=trim(path)//"/WaterAbsorption",name="DiffusionEq")
        fname=trim(path)//"/WaterAbsorption"
    endif

    call openArray(f%fh, obj%WaterAbsorbingPower)
    call openArray(f%fh, obj%WaterPotential)
    call openArray(f%fh, obj%TurgorPressure)
    call openArray(f%fh, obj%WaterContent)
    call openArray(f%fh, obj%WaterMass)
    call openArray(f%fh, obj%Conductivity)
    call openArray(f%fh, obj%YoungsModulus)
    call openArray(f%fh, obj%PoissonsRatio)
    call openArray(f%fh, obj%Permiability)
    call openArray(f%fh, obj%PorePressure)
    call openArray(f%fh, obj%a_Psi_val)
    call openArray(f%fh, obj%a_P_val)
    call openArray(f%fh, obj%theta_eq_val)
    call openArray(f%fh, obj%Psi_eq_val)
    call openArray(f%fh, obj%a_E_val)
    call openArray(f%fh, obj%a_v_val)
    call openArray(f%fh, obj%E_eq_val)
    call openArray(f%fh, obj%v_eq_val)
    call openArray(f%fh, obj%theta_ps_val) ! computed from other )
    read(f%fh, '(A)' ) obj%Name
    read(f%fh,*) obj%timestep
    read(f%fh,*) obj%dt
    call f%close()

end subroutine
! ################################################

!#####################################
subroutine flushWaterAbsorption(obj,path,name)
    class(WaterAbsorption_),intent(inout) :: obj
    character(*),optional,intent(in)::path
    character(*),optional,intent(in) :: Name
    character(200) ::fname
    character(200) ::fs
    integer(int32) :: fh
    logical :: restart=.true.
    type(IO_) :: f

    fs=trim(str(obj%flushstep) )
    if(.not. present(path) )then
        print *, " exportWaterAbsorption ERROR >> .not. present(path)"
        stop
    endif

    if(present(name) )then
        fname=trim(path)//"/"//trim(adjustl(name))
        call system("mkdir -p "//trim(path)//"/"//trim(adjustl(name)) )

        if(associated(obj%water) )then
            call obj%Water%save(path = "./"//trim(fname)//"/water", name="water")
        endif

        if(associated(obj%tissue) )then
            call obj%Tissue%save(path = "./"//trim(fname)//"/", name="tissue")
        endif

        call obj%DiffusionEq%export(path="./"//trim(fname)//"/DiffusionEq",restart=restart)
        call obj%FiniteDeform%export(path="./"//trim(fname)//"/FiniteDeform",restart=restart)

        if(associated(obj%a_Psi) )then
            call system("mkdir -p ./"//trim(fname)//"/a_Psi")
            call obj%a_Psi%export(      path=trim(fname)//"/a_Psi"  ,restart=restart)
        endif

        if(associated(obj%a_P) )then
            call system("mkdir -p ./"//trim(fname)//"/a_P")
            call obj%a_P%export(        path=trim(fname)//"/a_P"    ,restart=restart)
        endif

        if(associated(obj%theta_eq) )then
            call system("mkdir -p ./"//trim(fname)//"/theta_eq")
            call obj%theta_eq%export(   path=trim(fname)//"/theta_eq"   ,restart=restart)
        endif

        if(associated(obj%Psi_eq) )then
            call system("mkdir -p ./"//trim(fname)//"/Psi_eq")
            call obj%Psi_eq%export(     path=trim(fname)//"/Psi_eq"     ,restart=restart)
        endif

        if(associated(obj%a_E) )then
            call system("mkdir -p ./"//trim(fname)//"/a_E")
            call obj%a_E%export(        path=trim(fname)//"/a_E"    ,restart=restart)
        endif

        if(associated(obj%a_v) )then
            call system("mkdir -p ./"//trim(fname)//"/a_v")
            call obj%a_v%export(        path=trim(fname)//"/a_v"    ,restart=restart)
        endif

        if(associated(obj%E_eq) )then
            call system("mkdir -p ./"//trim(fname)//"/E_eq")
            call obj%E_eq%export(       path=trim(fname)//"/E_eq"   ,restart=restart)
        endif

        if(associated(obj%v_eq) )then
            call system("mkdir -p ./"//trim(fname)//"/v_eq")
            call obj%v_eq%export(       path=trim(fname)//"/v_eq"   ,restart=restart)
        endif

        call f%open("./",trim(path)//"/"//trim(adjustl(name)) ,"/WaterAbsorption.prop")

        write(f%fh,*) obj%WaterAbsorbingPower(:)
        write(f%fh,*) obj%WaterPotential(:)
        write(f%fh,*) obj%TurgorPressure(:)
        write(f%fh,*) obj%WaterContent(:)
        write(f%fh,*) obj%WaterMass(:)
        write(f%fh,*) obj%Conductivity(:)
        write(f%fh,*) obj%YoungsModulus(:)
        write(f%fh,*) obj%PoissonsRatio(:)
        write(f%fh,*) obj%Permiability(:)
        write(f%fh,*) obj%PorePressure(:)
        write(f%fh,*) obj%a_Psi_val(:)
        write(f%fh,*) obj%a_P_val(:)
        write(f%fh,*) obj%theta_eq_val(:)
        write(f%fh,*) obj%Psi_eq_val(:)
        write(f%fh,*) obj%a_E_val(:)
        write(f%fh,*) obj%a_v_val(:)
        write(f%fh,*) obj%E_eq_val(:)
        write(f%fh,*) obj%v_eq_val(:)
        write(f%fh,*) obj%theta_ps_val(:) ! computed from other variables
        write(f%fh,*) obj%Name
        write(f%fh,*) obj%timestep
        write(f%fh,*) obj%dt
        call f%close()
    else
        call system("mkdir -p "//trim(path)//"/WaterAbsorption")
        call f%open("./",trim(path)//"/WaterAbsorption","/WaterAbsorption.prop")
        fname=trim(path)//"/WaterAbsorption"
        write(f%fh,*) obj%WaterAbsorbingPower(:)
        write(f%fh,*) obj%WaterPotential(:)
        write(f%fh,*) obj%TurgorPressure(:)
        write(f%fh,*) obj%WaterContent(:)
        write(f%fh,*) obj%WaterMass(:)
        write(f%fh,*) obj%Conductivity(:)
        write(f%fh,*) obj%YoungsModulus(:)
        write(f%fh,*) obj%PoissonsRatio(:)
        write(f%fh,*) obj%Permiability(:)
        write(f%fh,*) obj%PorePressure(:)
        write(f%fh,*) obj%a_Psi_val(:)
        write(f%fh,*) obj%a_P_val(:)
        write(f%fh,*) obj%theta_eq_val(:)
        write(f%fh,*) obj%Psi_eq_val(:)
        write(f%fh,*) obj%a_E_val(:)
        write(f%fh,*) obj%a_v_val(:)
        write(f%fh,*) obj%E_eq_val(:)
        write(f%fh,*) obj%v_eq_val(:)
        write(f%fh,*) obj%theta_ps_val(:) ! computed from other variables
        write(f%fh,*) obj%Name
        write(f%fh,*) obj%timestep
        write(f%fh,*) obj%dt
        call f%close()

        if(associated(obj%water) )then
            call obj%Water%save(path = "./"//trim(fname)//"/water", name="water")
        endif

        if(associated(obj%tissue) )then
            call obj%Tissue%save(path = "./"//trim(fname)//"/", name="tissue")
        endif

        call obj%DiffusionEq%export(path="./"//trim(fname)//"/DiffusionEq",restart=restart)
        call obj%FiniteDeform%export(path="./"//trim(fname)//"/FiniteDeform",restart=restart)

        if(associated(obj%a_Psi) )then
            call system("mkdir -p ./"//trim(fname)//"/a_Psi")
            call obj%a_Psi%export(      path=trim(fname)//"/a_Psi"  ,restart=restart)
        endif

        if(associated(obj%a_P) )then
            call system("mkdir -p ./"//trim(fname)//"/a_P")
            call obj%a_P%export(        path=trim(fname)//"/a_P"    ,restart=restart)
        endif

        if(associated(obj%theta_eq) )then
            call system("mkdir -p ./"//trim(fname)//"/theta_eq")
            call obj%theta_eq%export(   path=trim(fname)//"/theta_eq"   ,restart=restart)
        endif

        if(associated(obj%Psi_eq) )then
            call system("mkdir -p ./"//trim(fname)//"/Psi_eq")
            call obj%Psi_eq%export(     path=trim(fname)//"/Psi_eq"     ,restart=restart)
        endif

        if(associated(obj%a_E) )then
            call system("mkdir -p ./"//trim(fname)//"/a_E")
            call obj%a_E%export(        path=trim(fname)//"/a_E"    ,restart=restart)
        endif

        if(associated(obj%a_v) )then
            call system("mkdir -p ./"//trim(fname)//"/a_v")
            call obj%a_v%export(        path=trim(fname)//"/a_v"    ,restart=restart)
        endif

        if(associated(obj%E_eq) )then
            call system("mkdir -p ./"//trim(fname)//"/E_eq")
            call obj%E_eq%export(       path=trim(fname)//"/E_eq"   ,restart=restart)
        endif

        if(associated(obj%v_eq) )then
            call system("mkdir -p ./"//trim(fname)//"/v_eq")
            call obj%v_eq%export(       path=trim(fname)//"/v_eq"   ,restart=restart)
        endif

    endif

end subroutine
! ################################################

!#####################################
subroutine saveWaterAbsorption(obj,path,name)
    class(WaterAbsorption_),intent(inout) :: obj
    character(*),optional,intent(in)::path
    character(*),optional,intent(in) :: Name
    character(200) ::fname
    integer(int32) :: fh
    logical :: restart=.true.
    type(IO_) :: f

    if(.not. present(path) )then
        print *, " exportWaterAbsorption ERROR >> .not. present(path)"
        stop
    endif

    if(present(name) )then
        call system("mkdir -p "//trim(path)//"/"//trim(adjustl(name)) )
	    fname=trim(path)//"/"//trim(adjustl(name))
        call f%open("./",trim(path)//"/"//trim(adjustl(name)) ,"/WaterAbsorption.prop")
        call obj%FiniteDeform%save(path=trim(path)//"/"//trim(adjustl(name)),name="FiniteDeform")
        call obj%DiffusionEq%save(path=trim(path)//"/"//trim(adjustl(name)),name="DiffusionEq")

    else
        call system("mkdir -p "//trim(path)//"/WaterAbsorption")
        call f%open("./",trim(path)//"/WaterAbsorption","/WaterAbsorption.prop")
        call obj%FiniteDeform%save(path=trim(path)//"/WaterAbsorption",name="FiniteDeform")
        call obj%DiffusionEq%save(path=trim(path)//"/WaterAbsorption",name="DiffusionEq")
        fname=trim(path)//"/WaterAbsorption"
    endif

    call WriteArray(f%fh, obj%WaterAbsorbingPower)
    call WriteArray(f%fh, obj%WaterPotential)
    call WriteArray(f%fh, obj%TurgorPressure)
    call WriteArray(f%fh, obj%WaterContent)
    call WriteArray(f%fh, obj%WaterMass)
    call WriteArray(f%fh, obj%Conductivity)
    call WriteArray(f%fh, obj%YoungsModulus)
    call WriteArray(f%fh, obj%PoissonsRatio)
    call WriteArray(f%fh, obj%Permiability)
    call WriteArray(f%fh, obj%PorePressure)
    call WriteArray(f%fh, obj%a_Psi_val)
    call WriteArray(f%fh, obj%a_P_val)
    call WriteArray(f%fh, obj%theta_eq_val)
    call WriteArray(f%fh, obj%Psi_eq_val)
    call WriteArray(f%fh, obj%a_E_val)
    call WriteArray(f%fh, obj%a_v_val)
    call WriteArray(f%fh, obj%E_eq_val)
    call WriteArray(f%fh, obj%v_eq_val)
    call WriteArray(f%fh, obj%theta_ps_val) ! computed from other )
    write(f%fh,*) obj%Name
    write(f%fh,*) obj%timestep
    write(f%fh,*) obj%dt
    call f%close()

end subroutine
! ################################################

! ################################################
subroutine importWaterAbsorption(obj,Water,Tissue,a_Psi,a_P,theta_eq,theta_ps,&
    Psi_eq,a_E,a_v,E_eq,v_eq)
    class(WaterAbsorption_),intent(inout) :: obj
    type(FEMDomain_),target,optional,intent(inout) :: Water,Tissue

    type(MaterialProp_),target,optional,intent(in):: a_Psi
    type(MaterialProp_),target,optional,intent(in):: a_P
    type(MaterialProp_),target,optional,intent(in):: theta_eq
    type(MaterialProp_),target,optional,intent(in):: theta_ps
    type(MaterialProp_),target,optional,intent(in):: Psi_eq
    type(MaterialProp_),target,optional,intent(in):: a_E
    type(MaterialProp_),target,optional,intent(in):: a_v
    type(MaterialProp_),target,optional,intent(in):: E_eq
    type(MaterialProp_),target,optional,intent(in):: v_eq

    if(associated(obj%DiffusionEq%FEMDomain) )then
        nullify(obj%DiffusionEq%FEMDomain)
    endif

    if(associated(obj%FiniteDeform%FEMDomain) )then
        nullify(obj%FiniteDeform%FEMDomain)
    endif
    if(present(water) )then
        obj%DiffusionEq%FEMDomain => Water
        obj%water => water
        Water%SolverType  = "DiffusionEq_"
    endif
    if(present(tissue))then
        obj%FiniteDeform%FEMDomain => Tissue
        Tissue%SolverType = "FiniteDeform_"
        obj%Tissue => Tissue
    endif

    print *, "[ImportFile]>"

    if(present(a_Psi) )then
        if(associated(obj%a_Psi) )then
            nullify(obj%a_Psi)
        endif
        obj%a_Psi => a_Psi
    endif

    if(present(a_P) )then
        if(associated(obj%a_P) )then
            nullify(obj%a_P)
        endif
        obj%a_P => a_P
    endif

    if(present(theta_eq) )then
        if(associated(obj%theta_eq) )then
            nullify(obj%theta_eq)
        endif
        obj%theta_eq => theta_eq
    endif
    if(present(Psi_eq) )then
        if(associated(obj%Psi_eq) )then
            nullify(obj%Psi_eq)
        endif
        obj%Psi_eq => Psi_eq
    endif
    if(present(a_E) )then
        if(associated(obj%a_E) )then
            nullify(obj%a_E)
        endif
        obj%a_E => a_E
    endif
    if(present(a_v) )then
        if(associated(obj%a_v) )then
            nullify(obj%a_v)
        endif
        obj%a_v => a_v
    endif

    if(present(E_eq) )then
        if(associated(obj%E_eq) )then
            nullify(obj%E_eq)
        endif
        obj%E_eq => E_eq
    endif

    if(present(v_eq) )then
        if(associated(obj%v_eq) )then
            nullify(obj%v_eq)
        endif
        obj%v_eq => v_eq
    endif

end subroutine importWaterAbsorption
!#####################################

!#####################################
subroutine runWaterAbsorption(obj,timestep,dt,SolverType,onlyInit,Only1st,Display,nr_tol,ReducedIntegration,&
    infinitesimal,interval,Name,restart,ID)
    class(WaterAbsorption_),intent(inout) :: obj
    character(*),intent(in) :: SolverType
    character(*),optional,intent(in) :: Name
    integer(int32) :: i,n,interv,coun
    logical,optional,intent(in) :: onlyInit,Only1st,Display,&
    ReducedIntegration,infinitesimal,restart
    real(real64) ,optional,intent(in) :: nr_tol
    integer(int32),intent(in) :: timestep
    integer(int32),optional,intent(in) :: interval,ID
    real(real64),optional,intent(in) :: dt

    if( present(restart) )then
        if(restart .eqv. .true. )then
            do i=1,timestep
                call obj%update(SolverType,Display=.true.,step=i,nr_tol=nr_tol,restart=.true.)
                print *, "Timestep :: ",i,"/",timestep
            enddo
            return
        endif
    endif
    if(present(Name) )then
        obj%water%FileName=Name
        obj%water%Name=Name
        obj%tissue%FileName=Name
        obj%tissue%Name=Name
    endif
    interv=input(default=0,option=interval)
    if(present(ReducedIntegration) )then
        obj%FiniteDeform%ReducedIntegration = ReducedIntegration
    endif

    obj%dt=input(default=1.0d0,option=dt)
    !obj%timestep=timestep

    obj%DiffusionEq%dt  = obj%dt
    obj%FiniteDeform%dt = obj%dt

    print *, "[ImportFile]>>>[Initialize]>"

    call obj%init(SolverType,Display=Display,nr_tol=nr_tol,infinitesimal=infinitesimal)
    print *, "[ImportFile]>>>[Initialize]>>>[1st Step]>"
    if(present(onlyInit) )then
        return
    endif
    ! Repeat over time
    coun=0
    do i=2,timestep

        if(coun==interv)then
            coun=0
            call obj%update(SolverType,Display=Display,step=i,nr_tol=nr_tol)
        else
            coun=coun+1
            call obj%update(SolverType,Display=.false.,step=i,nr_tol=nr_tol)
        endif
        print *, "Timestep :: ",i,"/",timestep

        !n=1
        !call showArray(obj%FiniteDeform%DeformVecGloTot,&
        !Name="test"//trim(adjustl(fstring(n))//".txt")  )

        if(present(only1st) )then
            return
        endif
    enddo

end subroutine runWaterAbsorption
!#####################################

!#####################################
subroutine initWaterAbsorption(obj,SolverType,Display,nr_tol,infinitesimal)
    class(WaterAbsorption_),intent(inout) :: obj
    type(Term_)             :: term
    character(*),intent(in) :: SolverType
    logical,optional,intent(in) :: Display,infinitesimal
    real(real64) ,optional,intent(in) :: nr_tol
    integer(int32) :: n

    call term%init()
    obj%volume = obj%FiniteDeform%getVolume()

    !open(113,file=trim(obj%tissue%name)//"x_length.txt",status="replace")
    !open(114,file=trim(obj%tissue%name)//"y_length.txt",status="replace")
    !open(115,file=trim(obj%tissue%name)//"z_length.txt",status="replace")
    call obj%updatePermiability()
    ! ###### Diffusion Part ###################
    call obj%DiffusionEq%Setup()
    call obj%DiffusionEq%Solve(SolverType=SolverType)
    if(present(Display) )then
        if(Display .eqv. .true.)then
            call DisplayDiffusionEq(obj%DiffusionEq,&
            DisplayMode=term%Gmsh,OptionalStep=1)
        endif
    endif

    ! ###### Diffusion Part ###################
    print *, "[ImportFile]>>>[Initialize]>>>[1st Step]>> water >>"

    ! debug
    ! only diffusion

    ! ###### Finite deformation part #############################
    !call obj%updateStiffness()
    !call obj%updateStiffness()
    !call obj%FiniteDeform%import(YoungsModulus=obj%YoungsModulus)
    !call obj%FiniteDeform%import(PoissonsRatio=obj%PoissonsRatio)
    !call obj%FiniteDeform%import(PorePressure=obj%PorePressure)
    ! ###### Finite deformation part #############################
    if(present(infinitesimal) )then
        obj%FiniteDeform%infinitesimal = infinitesimal
    endif
    call obj%FiniteDeform%DivideBC()
    call obj%FiniteDeform%Solve(SolverType=SolverType,nr_tol=nr_tol)
    call DisplayReactionForce(obj%FiniteDeform,obj%ID)
    if(present(Display) )then
        if(Display .eqv. .true.)then
            call DisplayDeformStress(obj%FiniteDeform,&
                DisplayMode="gmsh",OptionalStep=1)
            n=obj%FiniteDeform%step
            !call showArray(obj%FiniteDeform%DeformVecGloTot,&
            !Name="test"//trim(adjustl(fstring(n))//".txt")  )

        endif
    endif
    call obj%UpdateVolume()

    ! update mesh of Diffusion analysis added at 20200827
    obj%DiffusionEq%FEMDomain%Mesh%NodCoord = obj%FiniteDeform%FEMDomain%Mesh%NodCoord
    ! ###### Finite deformation part #############################
    print *, "[ImportFile]>>>[Initialize]>>>[1st Step]>>>>["

end subroutine initWaterAbsorption
!#####################################

!#####################################
subroutine updateWaterAbsorption(obj,SolverType,Display,step,nr_tol,restart)
    class(WaterAbsorption_),intent(inout) :: obj
    character(*),intent(in) :: SolverType
    logical,optional,intent(in) :: Display,restart
    integer(int32),optional,intent(in) :: step
    type(Term_)             :: term
    real(real64),optional,intent(in) :: nr_tol
    integer(int32) :: i,n,m

    call term%init()
    obj%volume = obj%FiniteDeform%getVolume()

    obj%DiffusionEq%FEMDomain%name="DiffusionEq"
    obj%FiniteDeform%FEMDomain%name="FiniteDeform"
    obj%DiffusionEq%FEMDomain%filename="DiffusionEq"
    obj%FiniteDeform%FEMDomain%filename="FiniteDeform"
    ! ###### update DIffusion parameters ############
    call obj%updatePermiability()
    call obj%DiffusionEq%import(Permiability=obj%Permiability)
    !do i=1,size(obj%water%mesh%nodCoord,1)
    !    write(1234,*) obj%water%mesh%nodCoord(i,1),obj%DiffusionEq%dt*dble(step), obj%DiffusionEq%UnknownVec(i)
    !enddo
    !write(1234,*) "  "
    ! ###### Update Diffusion Field over timesteps ###################
    call obj%DiffusionEq%Update()
    call obj%DiffusionEq%Solve(SolverType=SolverType,restart=restart)
    obj%DiffusionEq%step=step
    if(present(Display) )then
        if(Display .eqv. .true.)then

            call DisplayDiffusionEq(obj%DiffusionEq,&
            DisplayMode=term%Gmsh,OptionalStep=step,name="water")
        endif
    endif
    ! ###### Update Diffusion Field over timesteps ###################

    ! only diffusion debug
    call obj%updateStiffness()
    ! ###### update DIffusion parameters ############
    call obj%FiniteDeform%import(YoungsModulus=obj%YoungsModulus)
    call obj%FiniteDeform%import(PoissonsRatio=obj%PoissonsRatio)
    call obj%FiniteDeform%import(PorePressure=obj%PorePressure)
    ! ###### Update Finite Deformation over timesteps ###################
    obj%FiniteDeform%step=step

    call obj%FiniteDeform%UpdateInitConfig()
    call obj%FiniteDeform%UpdateBC()
    call obj%FiniteDeform%Solve(SolverType=SolverType,nr_tol=nr_tol,restart=restart)

    ! update mesh of Diffusion analysis added at 20200827
    obj%DiffusionEq%FEMDomain%Mesh%NodCoord = obj%FiniteDeform%FEMDomain%Mesh%NodCoord

    call DisplayReactionForce(obj%FiniteDeform,obj%ID)
    print *, "Timestep ",step
    if(present(Display) )then
        if(Display.eqv. .true.)then
            call DisplayDeformStress(obj%FiniteDeform,&
                DisplayMode="gmsh",OptionalStep=step,name="tissue")

            call obj%export(displacement=.true.)
        endif
    endif

    ! ###### Update Finite Deformation over timesteps ###################

end subroutine updateWaterAbsorption
!#####################################

subroutine updatePermiabilityWA(obj)
    class(WaterAbsorption_),intent(inout) :: obj
    integer(int32) :: i,j,n,m
    real(real64) :: k_0, a_hat, val(2),x
    n=size(obj%a_Psi_val)

    if(.not. allocated(obj%Permiability) )then
        allocate(obj%Permiability(n) )
    endif

    if(.not. allocated(obj%WaterContent) )then
        allocate(obj%WaterContent(n) )
        allocate(obj%WaterMass(n) )
        obj%WaterContent(:) = 0.0d0
        obj%WaterMass(:) = 0.0d0
    endif

    if(allocated(obj%DiffusionEq%UnknownValue) )then
        if(size(obj%DiffusionEq%UnknownValue,1)==n)then
            m=size(obj%DiffusionEq%UnknownValue,2)
            do i=1,n
                obj%WaterContent(i) = 0.0d0
                do j=1,m
                    ! get avarage value of wate content theta
                    obj%WaterContent(i) =obj%WaterContent(i)+ obj%DiffusionEq%UnknownValue(i,j)/dble(m)
                enddo
            enddo
        endif
    endif

    ! determine permiability for each element
    do i=1,n
        k_0=obj%Water%MaterialProp%matPara(obj%Water%Mesh%ElemMat(i),1 )
        val(1)=obj%a_Psi_val(i)
        val(2)=obj%a_Psi_val(i) - abs(obj%a_P_val(i))
        ! Hebiside step function
        !if(obj%WaterContent(i) <=obj%theta_ps_val(i) )then
        !    a_hat = val(1)
        !else
        !    a_hat=val(2)
        !endif
        ! use normalization
        x = 2.0d0 * (obj%WaterContent(i) - 0.50d0 )*100.0d0
        a_hat = obj%a_Psi_val(i) - 1.0d0 / (1.0d0 + exp(- x ) ) * abs(obj%a_P_val(i))
        !write(220,*) obj%WaterContent(i), a_hat*k_0
        obj%Permiability(i) = k_0 * a_hat
    enddo

    !write(100,*) obj%Permiability(:)
end subroutine updatePermiabilityWA
!#####################################

!#####################################
subroutine updateStiffnessWA(obj)
    class(WaterAbsorption_),intent(inout) :: obj
    integer(int32) :: i,j,n
    real(real64) :: a_E, E, E_eq, theta,theta_eq,a_P,theta_ps
    real(real64) :: a_v, v, v_eq,p,val(2)
    n=size(obj%water%mesh%elemNod,1)
    if(.not. allocated(obj%YoungsModulus) )then
        allocate(obj%YoungsModulus(n) )
    endif

    if(.not. allocated(obj%PoissonsRatio) )then
        allocate(obj%PoissonsRatio(n) )
    endif

    if(.not. allocated(obj%PorePressure) )then
        allocate(obj%PorePressure(n) )
        obj%PorePressure=0.0d0
    endif

    ! update these parameters considering volumetric water content
    do i=1,n
        ! import data for each element
        a_E = obj%a_E_val(i)
        a_v = obj%a_v_val(i)
        E_eq = obj%E_eq_val(i)
        v_eq = obj%v_eq_val(i)
        theta = obj%WaterContent(i)
        theta_eq = obj%theta_eq_val(i)
        theta_ps = obj%theta_ps_val(i)
        a_P = obj%a_P_val(i)

        ! compute parameters
        if(a_E > 0)then
            a_E= - abs(a_E)
        endif

        if( theta > 1.0d0)then
            theta = 1.0d0
        endif
        E = a_E * (theta - theta_eq) + E_eq
        v = a_v * (theta - theta_eq) + v_eq
        if(theta > theta_ps)then
            p=a_P*(theta - theta_ps)
        else
            p=0.0d0
        endif

        ! export parameters
        obj%YoungsModulus(i)= E
        obj%PoissonsRatio(i)= v
        obj%PorePressure(i) = p !- obj%PorePressure(i)
    enddo
    print *,"maxval(obj%a_P_val(:))",maxval(obj%a_P_val(:)),minval(obj%a_P_val(:))
    print *,"maxval(obj%WaterContent(:))",maxval(obj%WaterContent(:)),minval(obj%WaterContent(:))
    print *,"maxval(obj%theta_ps_val(:))",maxval(obj%theta_ps_val(:)),minval(obj%theta_ps_val(:))
    print *, "maxval(obj%PorePressure(:))",maxval(obj%PorePressure(:)),minval(obj%PorePressure(:))
    call obj%FiniteDeform%import( obj%YoungsModulus, obj%PoissonsRatio, obj%PorePressure )

end subroutine updateStiffnessWA
!#####################################

!#####################################
subroutine exportWaterAbsorption(obj,OptionalFileFormat,OptionalProjectName,FileHandle,&
    SolverType,MeshDimension,FileName,Name,regacy,with, displacement,restart,path)
    class(WaterAbsorption_),intent(inout) :: obj
    class(FEMDomain_),optional,intent(inout)::with
    character(*),optional,intent(in)::OptionalFileFormat,path
    character(*),optional,intent(in)::OptionalProjectName,SolverType,FileName
	character*4::FileFormat
	character(*),optional,intent(in) :: Name
	logical,optional,intent(in) :: regacy,displacement,restart
    character*200::ProjectName
	character*200 ::iFileName,fname
    integer(int32),allocatable::IntMat(:,:)
    real(real64),allocatable::RealMat(:,:)
    integer(int32),optional,intent(in)::FileHandle,MeshDimension
    integer(int32) :: fh,i,j,k,NumOfDomain,n,m,DimNum,GpNum,nn,fh_
    character*70 Msg
    type(IO_) :: f

    if(present(restart) )then
        if(.not. present(path) )then
            print *, " exportWaterAbsorption ERROR >> .not. present(path)"
            stop
        endif

        if(present(name) )then
            call system("mkdir -p "//trim(path)//"/"//trim(adjustl(name)) )
		    call f%open("./",trim(path)//"/"//trim(adjustl(name)) ,"/WaterAbsorption.prop")
            fname=trim(path)//"/"//trim(adjustl(name))
            write(f%fh,*) obj%WaterAbsorbingPower(:)
            write(f%fh,*) obj%WaterPotential(:)
            write(f%fh,*) obj%TurgorPressure(:)
            write(f%fh,*) obj%WaterContent(:)
            write(f%fh,*) obj%Conductivity(:)
            write(f%fh,*) obj%YoungsModulus(:)
            write(f%fh,*) obj%PoissonsRatio(:)
            write(f%fh,*) obj%Permiability(:)
            write(f%fh,*) obj%PorePressure(:)
            write(f%fh,*) obj%a_Psi_val(:)
            write(f%fh,*) obj%a_P_val(:)
            write(f%fh,*) obj%theta_eq_val(:)
            write(f%fh,*) obj%Psi_eq_val(:)
            write(f%fh,*) obj%a_E_val(:)
            write(f%fh,*) obj%a_v_val(:)
            write(f%fh,*) obj%E_eq_val(:)
            write(f%fh,*) obj%v_eq_val(:)
            write(f%fh,*) obj%theta_ps_val(:) ! computed from other variables
            write(f%fh,*) obj%Name
            write(f%fh,*) obj%timestep
            write(f%fh,*) obj%dt
            call f%close()

            call system("mkdir -p ./"//trim(fname)//"/FEMDomain")
            call f%open("./",trim(fname),"/FEMDomain/Water.prop")
            call obj%Water%export(path = "./"//trim(fname)//"/FEMDomain", restart=restart)
            call f%close()

            call f%open("./",trim(fname),"/FEMDomain/Tissue.prop")
            call obj%Tissue%export(path = "./"//trim(fname)//"/FEMDomain", restart=restart)
            call f%close()

            call system("mkdir -p ./"//trim(fname)//"/DiffusionEq")
            call f%open("./",trim(fname),"/DiffusionEq/DiffusionEq.prop")
            call obj%DiffusionEq%export(path="./"//trim(fname)//"/DiffusionEq",restart=restart)
            call f%close()

            call system("mkdir -p ./"//trim(fname)//"/FiniteDeform")
            call f%open("./",trim(fname),"/FiniteDeform/FiniteDeform.prop")
            call obj%FiniteDeform%export(path="./"//trim(fname)//"/FiniteDeform",restart=restart)
            call f%close()

            call system("mkdir -p ./"//trim(fname)//"/a_Psi")
            call obj%a_Psi%export(      path=trim(fname)//"/a_Psi"  ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/a_P")
            call obj%a_P%export(        path=trim(fname)//"/a_P"    ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/theta_eq")
            call obj%theta_eq%export(   path=trim(fname)//"/theta_eq"   ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/Psi_eq")
            call obj%Psi_eq%export(     path=trim(fname)//"/Psi_eq"     ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/a_E")
            call obj%a_E%export(        path=trim(fname)//"/a_E"    ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/a_v")
            call obj%a_v%export(        path=trim(fname)//"/a_v"    ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/E_eq")
            call obj%E_eq%export(       path=trim(fname)//"/E_eq"   ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/v_eq")
            call obj%v_eq%export(       path=trim(fname)//"/v_eq"   ,restart=restart)

        else
            call system("mkdir -p "//trim(path)//"/WaterAbsorption")
            call f%open("./",trim(path)//"/WaterAbsorption","/WaterAbsorption.prop")
            fname=trim(path)//"/WaterAbsorption"
            write(f%fh,*) obj%WaterAbsorbingPower(:)
            write(f%fh,*) obj%WaterPotential(:)
            write(f%fh,*) obj%TurgorPressure(:)
            write(f%fh,*) obj%WaterContent(:)
            write(f%fh,*) obj%Conductivity(:)
            write(f%fh,*) obj%YoungsModulus(:)
            write(f%fh,*) obj%PoissonsRatio(:)
            write(f%fh,*) obj%Permiability(:)
            write(f%fh,*) obj%PorePressure(:)
            write(f%fh,*) obj%a_Psi_val(:)
            write(f%fh,*) obj%a_P_val(:)
            write(f%fh,*) obj%theta_eq_val(:)
            write(f%fh,*) obj%Psi_eq_val(:)
            write(f%fh,*) obj%a_E_val(:)
            write(f%fh,*) obj%a_v_val(:)
            write(f%fh,*) obj%E_eq_val(:)
            write(f%fh,*) obj%v_eq_val(:)
            write(f%fh,*) obj%theta_ps_val(:) ! computed from other variables
            write(f%fh,*) obj%Name
            write(f%fh,*) obj%timestep
            write(f%fh,*) obj%dt
            call f%close()

            call system("mkdir -p ./"//trim(fname)//"/FEMDomain")
            call f%open("./",trim(fname),"/FEMDomain/Water.prop")
            call obj%Water%export(path = "./"//trim(fname)//"/FEMDomain", restart=restart)
            call f%close()

            call f%open("./",trim(fname),"/FEMDomain/Tissue.prop")
            call obj%Tissue%export(path = "./"//trim(fname)//"/FEMDomain", restart=restart)
            call f%close()

            call system("mkdir -p ./"//trim(fname)//"/DiffusionEq")
            call f%open("./",trim(fname),"/DiffusionEq/DiffusionEq.prop")
            call obj%DiffusionEq%export(path="./"//trim(fname)//"/DiffusionEq",restart=restart)
            call f%close()

            call system("mkdir -p ./"//trim(fname)//"/FiniteDeform")
            call f%open("./",trim(fname),"/FiniteDeform/FiniteDeform.prop")
            call obj%FiniteDeform%export(path="./"//trim(fname)//"/FiniteDeform",restart=restart)
            call f%close()

            call system("mkdir -p ./"//trim(fname)//"/a_Psi")
            call obj%a_Psi%export(      path=trim(fname)//"/a_Psi"  ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/a_P")
            call obj%a_P%export(        path=trim(fname)//"/a_P"    ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/theta_eq")
            call obj%theta_eq%export(   path=trim(fname)//"/theta_eq"   ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/Psi_eq")
            call obj%Psi_eq%export(     path=trim(fname)//"/Psi_eq"     ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/a_E")
            call obj%a_E%export(        path=trim(fname)//"/a_E"    ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/a_v")
            call obj%a_v%export(        path=trim(fname)//"/a_v"    ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/E_eq")
            call obj%E_eq%export(       path=trim(fname)//"/E_eq"   ,restart=restart)

            call system("mkdir -p ./"//trim(fname)//"/v_eq")
            call obj%v_eq%export(       path=trim(fname)//"/v_eq"   ,restart=restart)

        endif

        return
    endif

    if(present(displacement) )then
        if(displacement .eqv. .true.)then
            if(present(Name) )then
                open(12,file="x"//trim(adjustl(name)))
                open(13,file="y"//trim(adjustl(name)))
                open(14,file="z"//trim(adjustl(name)))
            else
                open(12,file="x_Disp.txt")
                open(13,file="y_Disp.txt")
                open(14,file="z_Disp.txt")
            endif
            n=size(obj%water%mesh%nodCoord,1)
            do i=1,n
                write(12,*) obj%water%mesh%nodCoord(i,1), obj%tissue%mesh%nodCoord(i,1)-obj%water%mesh%nodCoord(i,1)
                write(13,*) obj%water%mesh%nodCoord(i,2), obj%tissue%mesh%nodCoord(i,2)-obj%water%mesh%nodCoord(i,2)
                write(14,*) obj%water%mesh%nodCoord(i,3), obj%tissue%mesh%nodCoord(i,3)-obj%water%mesh%nodCoord(i,3)
            enddo
            close(12)
            close(13)
            close(14)
        endif
    endif
    ! bug exsits
    !call obj%Tissue%export(SolverType="FiniteDeform",Name=Name)
    !call obj%Water%export(SolverType="DiffusionEq",Name=Name)

end subroutine exportWaterAbsorption
!#####################################

!#####################################
subroutine displayWaterAbsorption(obj)
    class(WaterAbsorption_),intent(inout) :: obj

    ! implement how to display the results.

end subroutine displayWaterAbsorption
!#####################################

!#####################################
subroutine bakeWaterAbsorption(obj)
    class(WaterAbsorption_),intent(inout) :: obj
    integer(int32) :: i,j,k,l,n,m,NumOfMaterial,layer,in_num,NumOfLayer
    real(real64),allocatable :: matPara(:,:),info(:,:)
    integer(int32),allocatable :: key(:)
    type(Rectangle_) :: rect,mrect
    logical :: in_case
    real(real64) :: matparaval,coord(3),x_max(3),x_min(3)

    !call obj%tissue%bake(template="FiniteDeform_")
    !call obj%water%bake(template="DiffusionEq_")

    n=size(obj%water%mesh%ElemNod,1)

    ! initialize all parameters
    if(allocated(obj%WaterAbsorbingPower) )then
        deallocate(obj%WaterAbsorbingPower)
    endif
    allocate(obj%WaterAbsorbingPower(n))
    if(allocated(obj%WaterPotential) )then
        deallocate(obj%WaterPotential)
    endif
    allocate(obj%WaterPotential(n))
    if(allocated(obj%TurgorPressure) )then
        deallocate(obj%TurgorPressure)
    endif
    allocate(obj%TurgorPressure(n))
    if(allocated(obj%WaterContent) )then
        deallocate(obj%WaterContent)
    endif
    allocate(obj%WaterContent(n))
    if(allocated(obj%Conductivity) )then
        deallocate(obj%Conductivity)
    endif
    allocate(obj%Conductivity(n))
    if(allocated(obj%YoungsModulus) )then
        deallocate(obj%YoungsModulus)
    endif
    allocate(obj%YoungsModulus(n))
    if(allocated(obj%PoissonsRatio) )then
        deallocate(obj%PoissonsRatio)
    endif
    allocate(obj%PoissonsRatio(n))

    ! initialize parameters
    obj%WaterAbsorbingPower(:) = 0.0d0 ! kPa
    obj%WaterPotential(:) = 0.0d0 ! water
    obj%TurgorPressure(:) = 0.0d0 ! changes with theta
    obj%WaterContent(:) = 1.0d0 ! theta
    obj%Conductivity(:) = 1.0d0 ! changes with theta
    obj%YoungsModulus(:) = 1.0d0 ! E
    obj%PoissonsRatio(:) = 0.30d0 ! v
    if(.not. associated(obj%a_Psi) )then
        print *, "a_Psi is not associated."
        stop
    endif
    if(.not. associated(obj%a_P) )then
        print *, "a_P is not associated."
        stop
    endif
    if(.not. associated(obj%theta_eq) )then
        print *, "theta_eq is not associated."
        stop
    endif
    if(.not. associated(obj%Psi_eq) )then
        print *, "Psi_eq is not associated."
        stop
    endif
    if(.not. associated(obj%a_E) )then
        print *, "a_E is not associated."
        stop
    endif
    if(.not. associated(obj%a_v) )then
        print *, "a_v is not associated."
        stop
    endif
    if(.not. associated(obj%E_eq) )then
        print *, "E_eq is not associated."
        stop
    endif
    if(.not. associated(obj%v_eq) )then
        print *, "v_eq is not associated."
        stop
    endif

    ! a_Psi
    call obj%a_Psi%getValues(mesh=obj%Tissue%Mesh,Values=obj%a_Psi_val)
    call obj%a_P%getValues(mesh=obj%Tissue%Mesh,Values=obj%a_P_val)
    call obj%theta_eq%getValues(mesh=obj%Tissue%Mesh,Values=obj%theta_eq_val)
    call obj%Psi_eq%getValues(mesh=obj%Tissue%Mesh,Values=obj%Psi_eq_val)
    call obj%a_E%getValues(mesh=obj%Tissue%Mesh,Values=obj%a_E_val)
    call obj%a_v%getValues(mesh=obj%Tissue%Mesh,Values=obj%a_v_val)
    call obj%E_eq%getValues(mesh=obj%Tissue%Mesh,Values=obj%E_eq_val)
    call obj%v_eq%getValues(mesh=obj%Tissue%Mesh,Values=obj%v_eq_val)

    ! a_Psi should be negative values
    n=size(obj%theta_eq_val)
    do i=1,n
        obj%a_Psi_val(i) = - abs(obj%a_Psi_val(i))
    enddo

    ! a_Psi should be negative values
    n=size(obj%theta_eq_val)
    do i=1,n
        obj%a_E_val(i) = - abs(obj%a_E_val(i))
    enddo

    ! a_Psi should be negative values
    n=size(obj%theta_eq_val)
    do i=1,n
        obj%a_v_val(i) = abs(obj%a_v_val(i))
    enddo

    n=size(obj%theta_eq_val)
    if(.not.allocated(obj%theta_ps_val) )then
        allocate(obj%theta_ps_val(n) )
    endif

    do i=1,n
        obj%theta_ps_val(i)  =  (-  obj%Psi_eq_val(i) + obj%a_P_val(i) *  obj%theta_eq_val(i) )/obj%a_P_val(i)
    enddo

    !call showarray(mat=obj%WaterAbsorbingPower, Name="test.txt")

end subroutine bakeWaterAbsorption
!#####################################

!#####################################
subroutine gnuplotWaterAbsorption(obj,mode,ElemID)
    class(WaterAbsorption_),intent(in) :: obj
    character(*),intent(in) :: mode
    integer(int32),optional,intent(in) :: ElemID
    integer(int32) :: i,j,n,m
    real(real64) :: theta,psi(2)

    n=input(default=1,option=ElemID)
    if(mode=="OWCC" .or. mode=="all")then
        open(20,file="OWCC.txt")
        do i=0,100
            theta=dble(i)*obj%theta_eq_val(n)/100.0d0
            write(20,*) theta, (theta-obj%theta_eq_val(n))*obj%a_Psi_val(n) +obj%Psi_eq_val(n)
        enddo
        write(20,*) " "
        do i=0,100
            theta=dble(i)*obj%theta_eq_val(n)/100.0d0
            psi(1)=0.0d0
            psi(2)=obj%a_P_val(n)*(theta-obj%theta_ps_val(n))
            write(20,*) theta, maxval(psi)
        enddo
        write(20,*) " "
        close(20)
    endif

    if(mode=="YoungsModulus" .or. mode=="all")then
        open(20,file="YoungsModulus.txt")
        do i=0,100
            theta=dble(i)*obj%theta_eq_val(n)/100.0d0
            write(20,*) theta, (theta-obj%theta_eq_val(n))*obj%a_E_val(n) +obj%E_eq_val(n)
        enddo
        close(20)
    endif

    if(mode=="PoissonsRatio" .or. mode=="all")then
        open(20,file="PoissonsRatio.txt")
        do i=0,100
            theta=dble(i)*obj%theta_eq_val(n)/100.0d0
            write(20,*) theta, (theta-obj%theta_eq_val(n))*obj%a_v_val(n) +obj%v_eq_val(n)
        enddo
        close(20)
    endif
end subroutine gnuplotWaterAbsorption
!#####################################

!#####################################
subroutine UpdateVolumeWA(obj)
    class(WaterAbsorption_),intent(inout) :: obj
    real(real64),allocatable :: volume(:)
    real(real64) :: theta_old, theta_new,mw,vol_old, vol_new
    integer(int32) :: i,j,n,numelem

    numelem = size(obj%FiniteDeform%FEMDomain%Mesh%ElemNod,1)
    volume = obj%FiniteDeform%getVolume()
    if(.not. allocated(obj%WaterMass) )then
        allocate(obj%WaterMass(numelem))
    endif

    do i=1,numelem
        theta_old = obj%WaterContent(i)
        vol_old = obj%volume(i)
        mw = theta_old*vol_old
        theta_new = mw/volume(i)
        obj%WaterMass(i) = mw
        obj%WaterContent(i) = theta_new
    enddo
    obj%volume = obj%FiniteDeform%getVolume()

end subroutine
!#####################################

end module