Changeset 3553 for trunk/LMDZ.COMMON/libf/evolution/math_mod.F90

Timestamp:

Dec 16, 2024, 5:30:48 PM (7 months ago)

Author:

jbclement

Message:

PEM:
Addition of the shifting of the soil temperature profile to follow the surface evolution due to ice condensation/sublimation + small cleanings/improvements.
JBC

File:

• trunk/LMDZ.COMMON/libf/evolution/math_mod.F90 (modified) (6 diffs)

trunk/LMDZ.COMMON/libf/evolution/math_mod.F90

@@ -28 +28 @@
 implicit none
 
-integer, intent(in) :: nz      ! number of layer
-real,    intent(in) :: z(nz)   ! depth layer
-real,    intent(in) :: y(nz)   ! function which needs to be derived
-real,    intent(in) :: y0,ybot ! boundary conditions
-real, intent(out) :: dzY(nz) ! derivative of y w.r.t depth
-! local
+! Inputs
+!-------
+integer,             intent(in) :: nz       ! number of layer
+real, dimension(nz), intent(in) :: z        ! depth layer
+real, dimension(nz), intent(in) :: y        ! function which needs to be derived
+real,                intent(in) :: y0, ybot ! boundary conditions
+! Outputs
+!--------
+real, dimension(nz), intent(out) :: dzY ! derivative of y w.r.t depth
+! Local variables
+!----------------
 integer :: j
 real    :: hm, hp, c1, c2, c3
@@ -69 +74 @@
 implicit none
 
-integer, intent(in) :: nz
-real,    intent(in) :: z(nz), y(nz), y0, yNp1
-real, intent(out) :: yp2(nz)
+! Inputs
+!-------
+integer,             intent(in) :: nz
+real, dimension(nz), intent(in) :: z, y
+real,                intent(in) :: y0, yNp1
+! Outputs
+!--------
+real, dimension(nz), intent(out) :: yp2
+! Local variables
+!----------------
 integer :: j
 real    :: hm, hp, c1, c2, c3
@@ -102 +114 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-  implicit none
-
-integer, intent(in) :: nz, j
-real,    intent(in) :: z(nz), y(nz)
+implicit none
+
+! Inputs
+!-------
+integer,             intent(in) :: nz, j
+real, dimension(nz), intent(in) :: z, y
+! Outputs
+!--------
 real, intent(out) :: dy_zj
+! Local viariables
+!-----------------
 real :: h1, h2, c1, c2, c3
 
@@ -135 +153 @@
 implicit none
 
-integer, intent(in) :: nz, i1, i2
-real,    intent(in) :: y(nz), z(nz)
+! Inputs
+!-------
+integer,             intent(in) :: nz, i1, i2
+real, dimension(nz), intent(in) :: y, z
+! Outputs
+!--------
 real, intent(out) :: integral
-integer i
-real dz(nz)
+! Local viariables
+!-----------------
+integer             :: i
+real, dimension(nz) :: dz
 
 dz(1) = (z(2) - 0.)/2
@@ -163 +187 @@
 implicit none
 
+! Inputs
+!-------
 real, intent(in) :: y1, y2 ! difference between surface water density and at depth  [kg/m^3]
 real, intent(in) :: z1, z2 ! depth [m]
-real, intent(out) :: zr    ! depth at which we have zero
+! Outputs
+!--------
+real, intent(out) :: zr ! depth at which we have zero
 
 zr = (y1*z2 - y2*z1)/(y1 - y2)
@@ -171 +199 @@
 END SUBROUTINE findroot
 
+!=======================================================================
+
+SUBROUTINE solve_tridiag(a,b,c,d,n,x,error)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Solve a tridiagonal system Ax = d using the Thomas' algorithm
+!!!          a: sub-diagonal
+!!!          b: main diagonal
+!!!          c: super-diagonal
+!!!          d: right-hand side
+!!!          x: solution
+!!!
+!!! Author: JBC
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+implicit none
+
+! Inputs
+!-------
+integer,                intent(in) :: n
+real, dimension(n),     intent(in) :: b, d
+real, dimension(n - 1), intent(in) :: a, c
+! Outputs
+!--------
+real, dimension(n), intent(out) :: x
+integer,            intent(out) :: error
+! Local viariables
+!-----------------
+integer            :: i
+real               :: m
+real, dimension(n) :: cp, dp
+
+! Check stability: diagonally dominant condition
+error = 0
+if (abs(b(1)) < abs(c(1))) then
+    error = 1
+    return
+endif
+do i = 2,n - 1
+    if (abs(b(i)) < abs(a(i - 1)) + abs(c(i))) then
+        error = 1
+        return
+    endif
+enddo
+if (abs(b(n)) < abs(a(n - 1))) then
+    error = 1
+    return
+endif
+
+! Initialization
+cp(1) = c(1)/b(1)
+dp(1) = d(1)/b(1)
+
+! Forward sweep
+do i = 2,n - 1
+    m = b(i) - a(i - 1)*cp(i - 1)
+    cp(i) = c(i)/m
+    dp(i) = (d(i) - a(i - 1)*dp(i - 1))/m
+enddo
+m = b(n) - a(n - 1)*cp(n - 1)
+dp(n) = (d(n) - a(n - 1)*dp(n - 1))/m
+
+! Backward substitution
+x(n) = dp(n)
+do i = n - 1,1,-1
+    x(i) = dp(i) - cp(i)*x(i + 1)
+enddo
+
+END SUBROUTINE solve_tridiag
+
+!=======================================================================
+
+SUBROUTINE solve_steady_heat(n,z,mz,kappa,mkappa,T_left,q_right,T)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!!!
+!!! Purpose: Solve 1D steady-state heat equation with space-dependent diffusivity
+!!!          Left boudary condition : prescribed temperature T_left
+!!!          Right boudary condition: prescribed thermal flux q_right
+!!!
+!!!          z     : grid points
+!!!          mz    : mid-grid points
+!!!          kappa : thermal diffusivity at grid points
+!!!          mkappa: thermal diffusivity at mid-grid points
+!!!
+!!! Author: JBC
+!!!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+use abort_pem_mod, only: abort_pem
+
+implicit none
+
+! Inputs
+!-------
+integer,                intent(in) :: n
+real, dimension(n),     intent(in) :: z, kappa
+real, dimension(n - 1), intent(in) :: mz, mkappa
+real,                   intent(in) :: T_left, q_right
+! Outputs
+!--------
+real, dimension(n), intent(out) :: T
+! Local viariables
+!-----------------
+integer                :: i, error
+real, dimension(n)     :: b, d
+real, dimension(n - 1) :: a, c
+
+! Initialization
+a = 0.; b = 0.; c = 0.; d = 0.
+
+! Left boundary condition (Dirichlet: prescribed temperature)
+b(1) = 1.
+d(1) = T_left
+
+! Internal points
+do i = 2,n - 1
+    a(i - 1) = -mkappa(i - 1)/((mz(i) - mz(i - 1))*(z(i) - z(i - 1)))
+    c(i) = -mkappa(i)/((mz(i) - mz(i - 1))*(z(i + 1) - z(i)))
+    b(i) = -(a(i - 1) + c(i))
+enddo
+
+! Right boundary condition (Neumann: prescribed temperature)
+a(n - 1) = kappa(n - 1)/(z(n) - z(n - 1))
+b(n) = -kappa(n)/(z(n) - z(n - 1))
+d(n) = q_right
+
+! Solve the tridiagonal linear system with the Thomas' algorithm
+call solve_tridiag(a,b,c,d,n,T,error)
+if (error /= 0) call abort_pem("solve_steady_heat","Unstable solving!",1)
+
+END SUBROUTINE solve_steady_heat
+
 end module math_mod