!        Generated by TAPENADE     (INRIA, Ecuador team)
!  tapenade 3.x
!
!  Differentiation of lz in forward (tangent) mode:
!   variations   of useful results: lz
!   with respect to varying inputs: *temp *propunit
!   Plus diff mem management of: temp:in propunit:in
! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!>    @brief calculates effective thermal conductivity of the two phase
!>    @param[in] i cell index, direction I0
!>    @param[in] j cell index, direction J0
!>    @param[in] k cell index, direction K0
!>    @param[in] ismpl local sample index
!>    @return  thermal conductivity                lz[W/(m*K)]
!>    @details
!>    calculates effective thermal conductivity of the two phase\n
!>    system matrix-porosity .\n
!>    input:\n
!>      porosity                            porlocal [-]\n
!>      pressure                            plocal [Mpa]\n
!>      temperature                         tlocal in [C]\n
DOUBLE PRECISION FUNCTION g_LZ(i, j, k, ismpl, lz)
  USE ARRAYS

  USE g_ARRAYS

  USE ICE
  USE MOD_TEMP
  IMPLICIT NONE
  INTEGER :: i, j, k, ui, ismpl
  DOUBLE PRECISION :: tlocal, lsolid, lfluid, porlocal, fm, fi, ff
  DOUBLE PRECISION :: g_tlocal, g_lsolid, g_lfluid, g_porlocal, &
& g_fm, g_fi, g_ff
  DOUBLE PRECISION :: t0, theta, dtheta, w0, lice
  DOUBLE PRECISION :: g_theta, g_dtheta, g_lice
  EXTERNAL LAMF, LAMM, LAMI, POR
  EXTERNAL g_LAMF, g_LAMM, g_LAMI, g_POR
  DOUBLE PRECISION :: LAMM, LAMF, LAMI, POR
  DOUBLE PRECISION :: g_LAMM, g_LAMF, g_LAMI, g_POR
  INTRINSIC ABS
  DOUBLE PRECISION :: abs0
  DOUBLE PRECISION :: pwr1
  DOUBLE PRECISION :: g_pwr1
  DOUBLE PRECISION :: pwr2
  DOUBLE PRECISION :: g_pwr2
  DOUBLE PRECISION :: pwr3
  DOUBLE PRECISION :: g_pwr3
  DOUBLE PRECISION :: tmpresult
  DOUBLE PRECISION :: temp0
  DOUBLE PRECISION :: lz
  ui = uindex(i, j, k)
  g_tlocal = g_temp(i, j, k, ismpl)
  tlocal = temp(i, j, k, ismpl)
  t0 = liq(i, j, k)
  IF (liq(i, j, k) - sol(i, j, k) .GE. 0.) THEN
    abs0 = liq(i, j, k) - sol(i, j, k)
  ELSE
    abs0 = -(liq(i, j, k)-sol(i, j, k))
  END IF
  w0 = abs0/2.d0
  CALL g_FTHETA(tlocal, g_tlocal, t0, w0, theta, g_theta, dtheta, &
&           g_dtheta, ismpl)
  g_lfluid = g_LAMF(i, j, k, ismpl, lfluid)
  g_lice = g_LAMI(i, j, k, ismpl, lice)
  g_lsolid = g_propunit(ui, idx_lz, ismpl)
  lsolid = propunit(ui, idx_lz, ismpl)
  g_lsolid = g_LAMM(lsolid, g_lsolid, tlocal, g_tlocal, tref, &
&   ismpl, tmpresult)
  lsolid = tmpresult
  g_porlocal = g_POR(i, j, k, ismpl, porlocal)
  g_fm = -g_porlocal
  fm = 1.d0 - porlocal
  g_ff = theta*g_porlocal + porlocal*g_theta
  ff = porlocal*theta
  g_fi = g_porlocal - g_ff
  fi = porlocal - ff
  temp0 = lsolid**fm
  IF (lsolid .LE. 0.0 .AND. (fm .EQ. 0.0 .OR. fm .NE. INT(fm))) THEN
    g_pwr1 = 0.D0
  ELSE IF (lsolid .LE. 0.0) THEN
    g_pwr1 = fm*lsolid**(fm-1)*g_lsolid
  ELSE
    g_pwr1 = fm*lsolid**(fm-1)*g_lsolid + temp0*LOG(lsolid)*g_fm
  END IF
  pwr1 = temp0
  temp0 = lfluid**ff
  IF (lfluid .LE. 0.0 .AND. (ff .EQ. 0.0 .OR. ff .NE. INT(ff))) THEN
    g_pwr2 = 0.D0
  ELSE IF (lfluid .LE. 0.0) THEN
    g_pwr2 = ff*lfluid**(ff-1)*g_lfluid
  ELSE
    g_pwr2 = ff*lfluid**(ff-1)*g_lfluid + temp0*LOG(lfluid)*g_ff
  END IF
  pwr2 = temp0
  temp0 = lice**fi
  IF (lice .LE. 0.0 .AND. (fi .EQ. 0.0 .OR. fi .NE. INT(fi))) THEN
    g_pwr3 = 0.D0
  ELSE IF (lice .LE. 0.0) THEN
    g_pwr3 = fi*lice**(fi-1)*g_lice
  ELSE
    g_pwr3 = fi*lice**(fi-1)*g_lice + temp0*LOG(lice)*g_fi
  END IF
  pwr3 = temp0
  g_lz = pwr3*(pwr2*g_pwr1+pwr1*g_pwr2) + pwr1*pwr2*g_pwr3
  lz = pwr1*pwr2*pwr3
  RETURN
END FUNCTION g_LZ

