! $Id$
!
! This module solves the radiative diffusion implicitly thanks
! to an Alternate Direction Implicit Scheme (ADI) in a D'Yakonov
! form
! lambda_x T(n+1/2) = lambda_x + lambda_z
! lambda_z T(n+1) = T(n+1/2)
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lADI = .true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 1
! COMMUNICATED AUXILIARIES 1
!
!***************************************************************
module ImplicitPhysics
!
use Cparam
use Cdata
use General, only: keep_compiler_quiet
use Messages, only: svn_id, fatal_error, warning
use General, only: tridag, cyclic
!
implicit none
!
include 'implicit_physics.h'
!
interface heatcond_TT ! Overload subroutine `hcond_TT' function
module procedure heatcond_TT_0d ! get one value (hcond, dhcond)
module procedure heatcond_TT_1d ! get 1d-arrays (hcond, dhcond)
module procedure heatcond_TT_2d ! get 2d-arrays (hcond, dhcond)
end interface
!
real, pointer :: hcond0, Fbot, hcond1, hcond2, widthlnTT
logical, pointer :: lADI_mixed, lmultilayer
real :: Tbump, Kmax, Kmin, hole_slope, hole_width, hole_alpha
real :: dx_2, dy_2, dz_2, cp1
logical :: lyakonov=.true.
!
real, dimension(mz) :: hcondz, dhcondz
!
contains
!***********************************************************************
subroutine register_implicit_physics()
!
! Initialise variables which should know that we solve the
! compressible hydro equations: ilnrho; increase nvar accordingly.
!
! 03-mar-2010/dintrans: coded
!
use FArrayManager, only: farray_register_auxiliary
!
call farray_register_auxiliary('TTold',iTTold,communicated=.true.)
print*, 'iTTold=', iTTold
!
! Identify version number (generated automatically by SVN).
!
if (lroot) call svn_id( &
"$Id$")
!
endsubroutine register_implicit_physics
!***********************************************************************
subroutine initialize_implicit_physics(f)
!
use SharedVariables, only: get_shared_variable
use MpiComm, only: stop_it
use EquationOfState, only: get_cp1
use Gravity, only: z1, z2
use Sub, only: step,der_step,write_zprof
!
implicit none
!
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(:), pointer :: hole_params
real, dimension(mz) :: profz
!
call get_shared_variable('hcond0', hcond0, caller='initialize_implicit_physics')
call get_shared_variable('hcond1', hcond1)
call get_shared_variable('hcond2', hcond2)
call get_shared_variable('widthlnTT', widthlnTT)
call get_shared_variable('lmultilayer', lmultilayer)
print*,'***********************************'
print*, hcond0, hcond1, hcond2, widthlnTT, lmultilayer
print*,'***********************************'
call get_shared_variable('Fbot', Fbot)
call get_shared_variable('lADI_mixed', lADI_mixed)
call get_shared_variable('hole_params', hole_params)
Tbump=hole_params(1)
Kmin=hole_params(2)
Kmax=hole_params(3)
hole_slope=hole_params(4)
hole_width=hole_params(5)
hole_alpha=(Kmax-Kmin)/(pi/2.+atan(hole_slope*hole_width**2))
if (lroot .and. ldebug) then
print*, '************ hole parameters ************'
print*,'Tbump, Kmax, Kmin, hole_slope, hole_width, hole_alpha=', &
Tbump, Kmax, Kmin, hole_slope, hole_width, hole_alpha
print*, '*****************************************'
endif
!
if (lrun) then
! hcondADI is dynamically shared with boundcond() for the 'c3' BC
call heatcond_TT(f(:,4,n1,ilnTT), hcondADI)
else
hcondADI=spread(Kmax, 1, mx)
endif
!
! variables that are needed everywhere in this module
!
call get_cp1(cp1)
if (dx>0.) then
dx_2 = 1.0 / dx**2
else
dx_2 = 0.0
endif
if (dy>0.) then
dy_2 = 1.0 / dy**2
else
dy_2 = 0.0
endif
if (dz>0.) then
dz_2 = 1.0 / dz**2
else
dz_2 = 0.0
endif
!
if (lrun) then
if (lmultilayer) then
profz = 1. + (hcond1-1.)*step(z,z1,-widthlnTT) &
+ (hcond2-1.)*step(z,z2,widthlnTT)
hcondz = hcond0*profz
dhcondz = (hcond1-1.)*der_step(z,z1,-widthlnTT) &
+ (hcond2-1.)*der_step(z,z2,widthlnTT)
dhcondz = hcond0*dhcondz
else
hcondz=hcond0
dhcondz=0.0
endif
call write_zprof('hcond',hcondz)
call write_zprof('dhcond',dhcondz)
endif
!
endsubroutine initialize_implicit_physics
!***********************************************************************
subroutine calc_heatcond_ADI(f)
!
! 10-sep-07/gastine+dintrans: wrapper to the two possible ADI subroutines
! ADI_Kconst: constant radiative conductivity
! ADI_Kprof: radiative conductivity depends on T, i.e. hcond(T)
! 02/05/14-dintrans: added the polytropic superposed layers (MPI or not)
!
implicit none
!
real, dimension(mx,my,mz,mfarray) :: f
!
if (hcond0 /= impossible) then
!
! polytropic setup (single layer or superposed layers)
!
if (nx == 1) then
if (lmultilayer) then
call crank_Kprof(f)
else
call crank_Kconst(f)
endif
else
if (nprocz>1) then
if (lmultilayer) then
call ADI_poly_MPI(f)
else
call ADI_Kconst_MPI(f)
endif
else
if (lyakonov) then
if (lmultilayer) then
call ADI_poly(f)
else
call ADI_Kconst_yakonov(f)
endif
else
call ADI_Kconst(f)
endif
endif
endif
else
!
! kappa-mechanism with a conductivity hollow
!
if (nx == 1) then
if (lADI_mixed) then
call ADI_Kprof_1d_mixed(f)
else
call ADI_Kprof_1d(f)
endif
else
if (nprocz>1) then
call ADI_Kprof_MPI(f)
else
if (lADI_mixed) then
call ADI_Kprof_mixed(f)
else
call ADI_Kprof(f)
endif
endif
endif
endif
!
endsubroutine calc_heatcond_ADI
!***********************************************************************
subroutine ADI_Kconst(f)
!
! 08-Sep-07/gastine+dintrans: coded
! 2-D ADI scheme for the radiative diffusion term (see
! Peaceman & Rachford 1955). Each direction is solved implicitly:
!
! (1-dt/2*Lambda_x)*T^(n+1/2) = (1+dt/2*Lambda_y)*T^n + source/2
! (1-dt/2*Lambda_y)*T^(n+1) = (1+dt/2*Lambda_x)*T^(n+1/2) + source/2
!
! where Lambda_x and Lambda_y denote diffusion operators and the source
! term comes from the explicit advance.
!
use EquationOfState, only: gamma, gamma_m1, cs2bot, cs2top
use Boundcond, only: update_ghosts
!
implicit none
!
integer :: i,j
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mx,mz) :: finter, source, rho, TT
real, dimension(nx) :: ax, bx, cx, wx, rhsx, workx
real, dimension(nz) :: az, bz, cz, wz, rhsz, workz
real :: aalpha, bbeta
logical :: err
character(len=255) :: msg
!
! first update all the ghost zones in the f-array
!
call update_ghosts(f)
!
TT=f(:,4,:,iTTold)
source=(f(:,4,:,ilnTT)-TT)/dt
if (ldensity) then
rho=exp(f(:,4,:,ilnrho))
else
rho=1.
endif
!
! row dealt implicitly
!
do j=n1,n2
wx=dt*gamma*hcond0*cp1/rho(l1:l2,j)
ax=-wx*dx_2/2.
bx=1.+wx*dx_2
cx=ax
rhsx=TT(l1:l2,j)+wx*dz_2/2.* &
(TT(l1:l2,j+1)-2.*TT(l1:l2,j)+TT(l1:l2,j-1)) &
+dt/2.*source(l1:l2,j)
!
! x boundary conditions: periodic
!
aalpha=cx(nx) ; bbeta=ax(1)
call cyclic(ax, bx, cx, aalpha, bbeta, rhsx, workx, nx)
finter(l1:l2,j)=workx
enddo
!
! finter must be periodic in the x-direction
!
finter(1:l1-1,:)=finter(l2i:l2,:)
finter(l2+1:mx,:)=finter(l1:l1i,:)
!
! columns dealt implicitly
!
do i=l1,l2
wz=dt*gamma*hcond0*cp1/rho(i,n1:n2)
az=-wz*dz_2/2.
bz=1.+wz*dz_2
cz=az
rhsz=finter(i,n1:n2)+wz*dx_2/2.* &
(finter(i+1,n1:n2)-2.*finter(i,n1:n2)+finter(i-1,n1:n2)) &
+dt/2.*source(i,n1:n2)
!
! z boundary conditions
! Always constant temperature at the top
!
bz(nz)=1. ; az(nz)=0.
rhsz(nz)=cs2top/gamma_m1
select case (bcz12(ilnTT,1))
! Constant temperature at the bottom
case ('cT')
bz(1)=1. ; cz(1)=0.
rhsz(1)=cs2bot/gamma_m1
! Constant flux at the bottom
case ('c1')
bz(1)=1. ; cz(1)=-1
rhsz(1)=dz*Fbot/hcond0
! we can use here the second-order relation for the first derivative:
! (T_{j+1}-T_{j_1})/2dz = dT/dz --> T_{j-1} = T_{j+1} - 2*dz*dT/dz
! and insert this expression in the difference relation to eliminate T_{j-1}:
! a_{j-1}*T_{j-1} + b_j T_j + c_{j+1}*T_{j+1} = RHS
! cz(1)=cz(1)+az(1)
! rhsz(1)=rhsz(1)-2.*az(1)*dz*Fbot/hcond0
case default
call fatal_error('ADI_Kconst','bcz on TT must be cT or c1')
endselect
!
call tridag(az, bz, cz, rhsz, workz, err, msg)
if (err) call warning('ADI_Kconst', trim(msg))
f(i,4,n1:n2,ilnTT)=workz
enddo
!
endsubroutine ADI_Kconst
!***********************************************************************
subroutine ADI_Kprof(f)
!
! 10-Sep-07/gastine+dintrans: coded
! 2-D ADI scheme for the radiative diffusion term where the radiative
! conductivity depends on T (uses heatcond_TT to compute hcond _and_
! dhcond). The ADI scheme is of Yakonov's form:
!
! (1-dt/2*J_x)*lambda = f_x(T^n) + f_y(T^n) + source
! (1-dt/2*J_y)*beta = lambda
! T^(n+1) = T^n + dt*beta
!
! where J_x and J_y denote Jacobian matrices df/dT.
!
use EquationOfState, only: gamma
!
implicit none
!
integer :: i,j
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mx,mz) :: source, hcond, dhcond, finter, val, TT, rho
real, dimension(nx) :: ax, bx, cx, wx, rhsx, workx
real, dimension(nz) :: az, bz, cz, wz, rhsz, workz
real :: aalpha, bbeta
logical :: err
character(len=255) :: msg
!
source=(f(:,4,:,ilnTT)-f(:,4,:,iTTold))/dt
! BC important not for the x-direction (always periodic) but for
! the z-direction as we must impose the 'c3' BC at the 2nd-order
! before going in the implicit stuff
call heatcond_TT(f(:,4,:,iTTold), hcond, dhcond)
call boundary_ADI(f(:,4,:,iTTold), hcond(:,n1))
TT=f(:,4,:,iTTold)
if (ldensity) then
rho=exp(f(:,4,:,ilnrho))
else
rho=1.
endif
!
! rows dealt implicitly
!
do j=n1,n2
wx=cp1*gamma/rho(l1:l2,j)
! ax=-dt/2*J_x for i=i-1 (lower diagonal)
ax=-dt*wx*dx_2/4.*(dhcond(l1-1:l2-1,j) &
*(TT(l1-1:l2-1,j)-TT(l1:l2,j)) &
+hcond(l1-1:l2-1,j)+hcond(l1:l2,j))
! bx=1-dt/2*J_x for i=i (main diagonal)
bx=1.+dt*wx*dx_2/4.*(dhcond(l1:l2,j) &
*(2.*TT(l1:l2,j)-TT(l1-1:l2-1,j) &
-TT(l1+1:l2+1,j))+2.*hcond(l1:l2,j) &
+hcond(l1+1:l2+1,j)+hcond(l1-1:l2-1,j))
! cx=-dt/2*J_x for i=i+1 (upper diagonal)
cx=-dt*wx*dx_2/4.*(dhcond(l1+1:l2+1,j) &
*(TT(l1+1:l2+1,j)-TT(l1:l2,j)) &
+hcond(l1:l2,j)+hcond(l1+1:l2+1,j))
! rhsx=f_y(T^n) + f_x(T^n) (Eq. 3.6)
! do first f_y(T^n)
rhsx=wx*dz_2/2.*((hcond(l1:l2,j+1) &
+hcond(l1:l2,j))*(TT(l1:l2,j+1) &
-TT(l1:l2,j))-(hcond(l1:l2,j) &
+hcond(l1:l2,j-1)) &
*(TT(l1:l2,j)-TT(l1:l2,j-1)))
! then add f_x(T^n)
rhsx=rhsx+wx*dx_2/2.*((hcond(l1+1:l2+1,j) &
+hcond(l1:l2,j))*(TT(l1+1:l2+1,j)-TT(l1:l2,j)) &
-(hcond(l1:l2,j)+hcond(l1-1:l2-1,j)) &
*(TT(l1:l2,j)-TT(l1-1:l2-1,j)))+source(l1:l2,j)
!
! x boundary conditions: periodic
aalpha=cx(nx) ; bbeta=ax(1)
call cyclic(ax,bx,cx,aalpha,bbeta,rhsx,workx,nx)
finter(l1:l2,j)=workx(1:nx)
enddo
!
! columns dealt implicitly
!
do i=l1,l2
wz=dt*cp1*gamma*dz_2/rho(i,n1:n2)
az=-wz/4.*(dhcond(i,n1-1:n2-1) &
*(TT(i,n1-1:n2-1)-TT(i,n1:n2)) &
+hcond(i,n1-1:n2-1)+hcond(i,n1:n2))
!
bz=1.+wz/4.*(dhcond(i,n1:n2)* &
(2.*TT(i,n1:n2)-TT(i,n1-1:n2-1) &
-TT(i,n1+1:n2+1))+2.*hcond(i,n1:n2) &
+hcond(i,n1+1:n2+1)+hcond(i,n1-1:n2-1))
!
cz=-wz/4.*(dhcond(i,n1+1:n2+1) &
*(TT(i,n1+1:n2+1)-TT(i,n1:n2)) &
+hcond(i,n1:n2)+hcond(i,n1+1:n2+1))
!
rhsz=finter(i,n1:n2)
!
! z boundary conditions
! Constant temperature at the top: T^(n+1)-T^n=0
bz(nz)=1. ; az(nz)=0.
rhsz(nz)=0.
! bottom
select case (bcz12(ilnTT,1))
! Constant temperature at the bottom: T^(n+1)-T^n=0
case ('cT')
bz(1)=1. ; cz(1)=0.
rhsz(1)=0.
! Constant flux at the bottom
case ('c3')
bz(1)=1. ; cz(1)=-1.
rhsz(1)=0.
case default
call fatal_error('ADI_Kprof','bcz on TT must be cT or c3')
endselect
!
call tridag(az,bz,cz,rhsz,workz,err,msg)
if (err) call warning('ADI_Kprof', trim(msg))
val(i,n1:n2)=workz(1:nz)
enddo
!
f(:,4,:,ilnTT)=f(:,4,:,iTTold)+dt*val
!
! update hcond used for the 'c3' condition in boundcond.f90
!
call heatcond_TT(f(:,4,n1,ilnTT), hcondADI)
!
endsubroutine ADI_Kprof
!***********************************************************************
subroutine ADI_Kprof_MPI(f)
!
! 15-jan-10/gastine: coded
! 2-D ADI scheme for the radiative diffusion term where the radiative
! conductivity depends on T (uses heatcond_TT to compute hcond _and_
! dhcond). The ADI scheme is of Yakonov's form:
!
! (1-dt/2*J_x)*lambda = f_x(T^n) + f_z(T^n) + source
! (1-dt/2*J_z)*beta = lambda
! T^(n+1) = T^n + dt*beta
!
! where J_x and J_z denote Jacobian matrices df/dT.
! 08-mar-2010/dintrans: added the case of a non-square domain (ibox-loop)
! 21-aug-2010/dintrans: simplified version that uses Anders' original
! transp_xz and transp_zx subroutines
!
use EquationOfState, only: gamma
use Mpicomm, only: transp_xz, transp_zx
use Boundcond, only: update_ghosts
!
implicit none
!
integer, parameter :: mzt=nzgrid+2*nghost
integer, parameter :: n1t=nghost+1, n2t=n1t+nzgrid-1
integer, parameter :: nxt=nx/nprocz
integer :: i ,j
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mx,mz) :: source, hcond, dhcond, finter, TT, rho
real, dimension(mzt,nxt) :: hcondt, dhcondt, fintert, TTt, rhot, valt
real, dimension(nx,nz) :: val
real, dimension(nx) :: ax, bx, cx, wx, rhsx, workx
real, dimension(nzgrid) :: az, bz, cz, wz, rhsz, workz
real :: aalpha, bbeta
logical :: err
character(len=255) :: msg
!
! It is necessary to communicate ghost-zones points between
! processors to ensure a correct transposition of these ghost
! zones. It is needed by rho,rhot and source,sourcet.
!
call update_ghosts(f)
source=(f(:,4,:,ilnTT)-f(:,4,:,iTTold))/dt
!
! BC important not for the x-direction (always periodic) but for
! the z-direction as we must impose the 'c3' BC at the 2nd-order
! before going in the implicit stuff
!
TT=f(:,4,:,iTTold)
call heatcond_TT(TT, hcond, dhcond)
call boundary_ADI(TT, hcond(:,n1))
if (ldensity) then
rho=exp(f(:,4,:,ilnrho))
else
rho=1.
endif
!
! rows dealt implicitly
!
do j=n1,n2
wx=cp1*gamma/rho(l1:l2,j)
! ax=-dt/2*J_x for i=i-1 (lower diagonal)
ax=-dt*wx*dx_2/4.*(dhcond(l1-1:l2-1,j) &
*(TT(l1-1:l2-1,j)-TT(l1:l2,j)) &
+hcond(l1-1:l2-1,j)+hcond(l1:l2,j))
! bx=1-dt/2*J_x for i=i (main diagonal)
bx=1.+dt*wx*dx_2/4.*(dhcond(l1:l2,j) &
*(2.*TT(l1:l2,j)-TT(l1-1:l2-1,j) &
-TT(l1+1:l2+1,j))+2.*hcond(l1:l2,j) &
+hcond(l1+1:l2+1,j)+hcond(l1-1:l2-1,j))
! cx=-dt/2*J_x for i=i+1 (upper diagonal)
cx=-dt*wx*dx_2/4.*(dhcond(l1+1:l2+1,j) &
*(TT(l1+1:l2+1,j)-TT(l1:l2,j)) &
+hcond(l1:l2,j)+hcond(l1+1:l2+1,j))
! rhsx=f_z(T^n) + f_x(T^n) (Eq. 3.6)
! do first f_z(T^n)
rhsx=wx*dz_2/2.*((hcond(l1:l2,j+1) &
+hcond(l1:l2,j))*(TT(l1:l2,j+1) &
-TT(l1:l2,j))-(hcond(l1:l2,j) &
+hcond(l1:l2,j-1)) &
*(TT(l1:l2,j)-TT(l1:l2,j-1)))
! then add f_x(T^n)
rhsx=rhsx+wx*dx_2/2.*((hcond(l1+1:l2+1,j) &
+hcond(l1:l2,j))*(TT(l1+1:l2+1,j)-TT(l1:l2,j)) &
-(hcond(l1:l2,j)+hcond(l1-1:l2-1,j)) &
*(TT(l1:l2,j)-TT(l1-1:l2-1,j)))+source(l1:l2,j)
!
! periodic boundary conditions in the x-direction
!
aalpha=cx(nx) ; bbeta=ax(1)
call cyclic(ax, bx, cx, aalpha, bbeta, rhsx, workx, nx)
finter(l1:l2,j)=workx(1:nx)
enddo
!
! do the transpositions x <--> z
!
call transp_xz(finter(l1:l2,n1:n2), fintert(n1t:n2t,:))
call transp_xz(rho(l1:l2,n1:n2), rhot(n1t:n2t,:))
call transp_xz(TT(l1:l2,n1:n2), TTt(n1t:n2t,:))
call heatcond_TT(TTt, hcondt, dhcondt)
!
do i=1,nxt
wz=dt*cp1*gamma*dz_2/rhot(n1t:n2t,i)
az=-wz/4.*(dhcondt(n1t-1:n2t-1,i) &
*(TTt(n1t-1:n2t-1,i)-TTt(n1t:n2t,i)) &
+hcondt(n1t-1:n2t-1,i)+hcondt(n1t:n2t,i))
!
bz=1.+wz/4.*(dhcondt(n1t:n2t,i)* &
(2.*TTt(n1t:n2t,i)-TTt(n1t-1:n2t-1,i) &
-TTt(n1t+1:n2t+1,i))+2.*hcondt(n1t:n2t,i) &
+hcondt(n1t+1:n2t+1,i)+hcondt(n1t-1:n2t-1,i))
!
cz=-wz/4.*(dhcondt(n1t+1:n2t+1,i) &
*(TTt(n1t+1:n2t+1,i)-TTt(n1t:n2t,i)) &
+hcondt(n1t:n2t,i)+hcondt(n1t+1:n2t+1,i))
!
rhsz=fintert(n1t:n2t,i)
!
! z boundary conditions
! Constant temperature at the top: T^(n+1)-T^n=0
!
bz(nzgrid)=1. ; az(nzgrid)=0.
rhsz(nzgrid)=0.
! bottom
select case (bcz12(ilnTT,1))
! Constant temperature at the bottom: T^(n+1)-T^n=0
case ('cT')
bz(1)=1. ; cz(1)=0.
rhsz(1)=0.
! Constant flux at the bottom
case ('c3')
bz(1)=1. ; cz(1)=-1.
rhsz(1)=0.
case default
call fatal_error('ADI_Kprof','bcz on TT must be cT or c3')
endselect
call tridag(az, bz, cz, rhsz, workz, err, msg)
if (err) call warning('ADI_Kprof', trim(msg))
valt(n1t:n2t,i)=workz(1:nzgrid)
enddo ! i
!
! come back on the grid (x,z)
!
call transp_zx(valt(n1t:n2t,:), val)
f(l1:l2,4,n1:n2,ilnTT)=f(l1:l2,4,n1:n2,iTTold)+dt*val
!
! update hcond used for the 'c3' condition in boundcond.f90
!
if (iproc==0) call heatcond_TT(f(:,4,n1,ilnTT), hcondADI)
!
endsubroutine ADI_Kprof_MPI
!***********************************************************************
subroutine boundary_ADI(f_2d, hcond)
!
! 13-Sep-07/gastine: computed two different types of boundary
! conditions for the implicit solver:
! - Always periodic in x-direction
! - Possibility to choose between 'cT' and 'c3' in z direction
! Note: 'c3' means that the flux is constant at the *bottom only*
!
implicit none
!
real, dimension(mx,mz) :: f_2d
real, dimension(mx), optional :: hcond
!
! x-direction: periodic
!
f_2d(1:l1-1,:)=f_2d(l2i:l2,:)
f_2d(l2+1:mx,:)=f_2d(l1:l1i,:)
!
! top boundary condition z=z(n2): always constant temperature
!
if (llast_proc_z) then
f_2d(:,n2+1)=2.*f_2d(:,n2)-f_2d(:,n2-1)
endif
!
! bottom bondary condition z=z(n1): constant T or imposed flux dT/dz
!
if (iproc==0) then
select case (bcz12(ilnTT,1))
case ('cT') ! constant temperature
f_2d(:,n1-1)=2.*f_2d(:,n1)-f_2d(:,n1+1)
case ('c3') ! constant flux
if (.not. present(hcond)) then
f_2d(:,n1-1)=f_2d(:,n1+1)+2.*dz*Fbot/hcond0
else
f_2d(:,n1-1)=f_2d(:,n1+1)+2.*dz*Fbot/hcond(:)
endif
endselect
endif
!
endsubroutine boundary_ADI
!***********************************************************************
subroutine crank_Kconst(f)
!
! 18-sep-07/dintrans: coded
! Implicit Crank Nicolson scheme in 1-D for a constant K.
!
use EquationOfState, only: gamma, gamma_m1, cs2bot, cs2top
use Boundcond, only: update_ghosts
!
implicit none
!
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mz) :: TT
real, dimension(nz) :: az, bz, cz, rhsz, wz
logical :: err
character(len=255) :: msg
!
call update_ghosts(f,iTT)
TT=f(4,4,:,iTT)
!
wz(:)=dt*dz_2*gamma*cp1*hcond0*exp(-f(4,4,n1:n2,ilnrho))
az(:)=-0.5*wz
bz(:)=1.+wz
cz(:)=az
do n=n1,n2
rhsz(n-nghost)=TT(n)+0.5*wz(n-nghost)*(TT(n+1)-2.*TT(n)+TT(n-1))
enddo
bz(nz)=1. ; az(nz)=0. ; rhsz(nz)=cs2top/gamma_m1 ! T = Ttop
if (bcz12(iTT,1)=='cT') then
bz(1)=1. ; cz(1)=0. ; rhsz(1)=cs2bot/gamma_m1 ! T = Tbot
else
! cz(1)=2.*cz(1) ; rhsz(1)=rhsz(1)+wz(1)*dz*Fbot/hcond0 ! T' = -Fbot/K
bz(1)=1. ; cz(1)=-1. ; rhsz(1)=dz*Fbot/hcond0 ! T' = -Fbot/K
endif
call tridag(az, bz, cz, rhsz, f(4,4,n1:n2,iTT), err, msg)
if (err) call warning('crank_Kconst', trim(msg))
!
endsubroutine crank_Kconst
!***********************************************************************
subroutine ADI_Kprof_1d(f)
!
! 18-sep-07/dintrans: coded
! Implicit 1-D case for a temperature-dependent conductivity K(T).
! Not really an ADI but keep the generic name for commodity.
!
use EquationOfState, only: gamma
!
implicit none
!
integer :: j, jj
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mz) :: source, rho, TT, hcond, dhcond
real, dimension(nz) :: a, b, c, rhs, work
real :: wz, hcondp, hcondm
logical :: err
character(len=255) :: msg
!
source=(f(4,4,:,ilnTT)-f(4,4,:,iTTold))/dt
rho=exp(f(4,4,:,ilnrho))
!
! need to set up the 'c3' BC at the 2nd-order before the implicit stuff
!
call heatcond_TT(f(4,4,:,iTTold), hcond, dhcond)
hcondADI=spread(hcond(1), 1, mx)
call boundary_ADI(f(:,4,:,iTTold), hcondADI)
TT=f(4,4,:,iTTold)
!
do j=n1,n2
jj=j-nghost
wz=dt*dz_2*gamma*cp1/rho(j)
hcondp=hcond(j+1)+hcond(j)
hcondm=hcond(j)+hcond(j-1)
!
a(jj)=-wz/4.*(hcondm-dhcond(j-1)*(TT(j)-TT(j-1)))
b(jj)=1.-wz/4.*(-hcondp-hcondm+dhcond(j)*(TT(j+1)-2.*TT(j)+TT(j-1)))
c(jj)=-wz/4.*(hcondp+dhcond(j+1)*(TT(j+1)-TT(j)))
rhs(jj)=wz/2.*(hcondp*(TT(j+1)-TT(j))-hcondm*(TT(j)-TT(j-1))) &
+dt*source(j)
!
! Always constant temperature at the top: T^(n+1)-T^n = 0
!
b(nz)=1. ; a(nz)=0.
rhs(nz)=0.
if (bcz12(ilnTT,1)=='cT') then
! Constant temperature at the bottom
b(1)=1. ; c(1)=0.
rhs(1)=0.
else
! Constant flux at the bottom: d/dz [T^(n+1)-T^n] = 0
b(1)=1. ; c(1)=-1.
rhs(1)=0.
endif
enddo
call tridag(a, b, c, rhs, work, err, msg)
if (err) call warning('ADI_Kprof_1d', trim(msg))
f(4,4,n1:n2,ilnTT)=f(4,4,n1:n2,iTTold)+work
!
! Update the bottom value of hcond used for the 'c3' BC in boundcond
!
call heatcond_TT(f(:,4,n1,ilnTT), hcondADI)
!
endsubroutine ADI_Kprof_1d
!***********************************************************************
subroutine ADI_Kconst_MPI(f)
!
! 04-sep-2009/dintrans: coded
! Parallel version of the ADI scheme for the K=cte case.
! Note: this is the parallelisation of the Yakonov form *only*.
!
use EquationOfState, only: gamma, gamma_m1, cs2bot, cs2top
use Mpicomm, only: transp_xz, transp_zx
use Boundcond, only: update_ghosts
!
implicit none
!
integer, parameter :: nxt=nx/nprocz
integer :: i,j
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mx,mz) :: finter, rho1, TT
real, dimension(nzgrid,nxt) :: fintert, rho1t, wtmp
real, dimension(nx) :: ax, bx, cx, wx, rhsx
real, dimension(nzgrid) :: az, bz, cz, wz, rhsz
real :: aalpha, bbeta
logical :: err
character(len=255) :: msg
!
call update_ghosts(f,iTT)
TT=f(:,4,:,iTT)
if (ldensity) then
rho1=exp(-f(:,4,:,ilnrho))
else
rho1=1.
endif
!
! Rows dealt implicitly
!
do j=n1,n2
wx=0.5*dt*cp1*gamma*hcond0*rho1(l1:l2,j)
ax=-wx*dx_2
bx=1.+2.*wx*dx_2
cx=ax
rhsx=TT(l1:l2,j) &
+wx*dx_2*(TT(l1+1:l2+1,j)-2.*TT(l1:l2,j)+TT(l1-1:l2-1,j)) &
+wx*dz_2*(TT(l1:l2,j+1)-2.*TT(l1:l2,j)+TT(l1:l2,j-1))
!
! x boundary conditions: periodic
!
aalpha=cx(nx) ; bbeta=ax(1)
call cyclic(ax, bx, cx, aalpha, bbeta, rhsx, finter(l1:l2,j), nx)
enddo
!
! Do the transpositions x <--> z
!
call transp_xz(finter(l1:l2,n1:n2), fintert)
call transp_xz(rho1(l1:l2,n1:n2), rho1t)
!
! Columns dealt implicitly
!
do i=1,nxt
wz=0.5*dz_2*dt*cp1*gamma*hcond0*rho1t(:,i)
az=-wz
bz=1.+2.*wz
cz=az
rhsz=fintert(:,i)
!
! z boundary conditions
! Always constant temperature at the top
!
bz(nzgrid)=1. ; az(nzgrid)=0. ; rhsz(nzgrid)=cs2top/gamma_m1
select case (bcz12(iTT,1))
case ('cT') ! Constant temperature at the bottom
bz(1)=1. ; cz(1)=0. ; rhsz(1)=cs2bot/gamma_m1
case ('c1') ! Constant flux at the bottom
bz(1)=1. ; cz(1)=-1. ; rhsz(1)=dz*Fbot/hcond0
case default
call fatal_error('ADI_Kconst_MPI','bcz on TT must be cT or c1')
endselect
!
call tridag(az, bz, cz, rhsz, wtmp(:,i), err, msg)
if (err) call warning('ADI_Kconst_MPI', trim(msg))
enddo
call transp_zx(wtmp, f(l1:l2,4,n1:n2,iTT))
!
endsubroutine ADI_Kconst_MPI
!***********************************************************************
subroutine heatcond_TT_2d(TT, hcond, dhcond)
!
! 07-Sep-07/gastine: computed 2-D radiative conductivity hcond(T) with
! its derivative dhcond=dhcond(T)/dT.
!
implicit none
!
real, dimension(:,:), intent(in) :: TT
real, dimension(:,:), intent(out) :: hcond
real, dimension(:,:), optional :: dhcond
!
hcond=hole_slope*(TT-Tbump-hole_width)*(TT-Tbump+hole_width)
if (present(dhcond)) &
dhcond=2.*hole_alpha/(1.+hcond**2)*hole_slope*(TT-Tbump)
hcond=Kmax+hole_alpha*(-pi/2.+atan(hcond))
!
endsubroutine heatcond_TT_2d
!***********************************************************************
subroutine heatcond_TT_1d(TT, hcond, dhcond)
!
! 18-Sep-07/dintrans: computed 1-D radiative conductivity
! hcond(T) with its derivative dhcond=dhcond(T)/dT.
!
implicit none
!
real, dimension(:), intent(in) :: TT
real, dimension(:), intent(out) :: hcond
real, dimension(:), optional :: dhcond
!
hcond=hole_slope*(TT-Tbump-hole_width)*(TT-Tbump+hole_width)
if (present(dhcond)) &
dhcond=2.*hole_alpha/(1.+hcond**2)*hole_slope*(TT-Tbump)
hcond=Kmax+hole_alpha*(-pi/2.+atan(hcond))
!
endsubroutine heatcond_TT_1d
!***********************************************************************
subroutine heatcond_TT_0d(TT, hcond, dhcond)
!
! 07-Sep-07/gastine: computed the radiative conductivity hcond(T)
! with its derivative dhcond=dhcond(T)/dT at a given temperature.
!
implicit none
!
real, intent(in) :: TT
real, intent(out) :: hcond
real, optional :: dhcond
!
hcond=hole_slope*(TT-Tbump-hole_width)*(TT-Tbump+hole_width)
if (present(dhcond)) &
dhcond=2.*hole_alpha/(1.+hcond**2)*hole_slope*(TT-Tbump)
hcond=Kmax+hole_alpha*(-pi/2.+atan(hcond))
!
endsubroutine heatcond_TT_0d
!***********************************************************************
subroutine ADI_Kprof_1d_mixed(f)
!
! 28-feb-10/dintrans: coded
! Simpler version where a part of the radiative diffusion term is
! computed during the explicit advance.
!
use EquationOfState, only: gamma
!
implicit none
!
integer :: j, jj
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mz) :: source, TT, hcond, dhcond, dLnhcond, chi
real, dimension(nz) :: a, b, c, rhs, work
real :: wz
logical :: err
character(len=255) :: msg
!
source=(f(4,4,:,ilnTT)-f(4,4,:,iTTold))/dt
call heatcond_TT(f(4,4,:,iTTold), hcond, dhcond)
!
! need to set up the 'c3' BC at the 2nd-order before the implicit stuff
!
hcondADI=spread(hcond(1), 1, mx)
call boundary_ADI(f(:,4,:,iTTold), hcondADI)
TT=f(4,4,:,iTTold)
if (ldensity) then
chi=cp1*hcond/exp(f(4,4,:,ilnrho))
else
chi=cp1*hcond
endif
dLnhcond=dhcond/hcond
!
do j=n1,n2
jj=j-nghost
wz=dt*dz_2*gamma*chi(j)
!
a(jj)=-wz/2.
b(jj)=1.-wz/2.*(-2.+dLnhcond(j)*(TT(j+1)-2.*TT(j)+TT(j-1)))
c(jj)=-wz/2.
rhs(jj)=wz*(TT(j+1)-2.*TT(j)+TT(j-1))+dt*source(j)
!
! Always constant temperature at the top: T^(n+1)-T^n = 0
!
b(nz)=1. ; a(nz)=0.
rhs(nz)=0.
if (bcz12(ilnTT,1)=='cT') then
! Constant temperature at the bottom
b(1)=1. ; c(1)=0.
rhs(1)=0.
else
! Constant flux at the bottom: d/dz [T^(n+1)-T^n] = 0
b(1)=1. ; c(1)=-1.
rhs(1)=0.
endif
enddo
!
call tridag(a, b, c, rhs, work, err, msg)
if (err) call warning('ADI_Kprof_1d_mixed', trim(msg))
f(4,4,n1:n2,ilnTT)=f(4,4,n1:n2,iTTold)+work
!
! Update the bottom value of hcond used for the 'c3' BC in boundcond
!
call heatcond_TT(f(:,4,n1,ilnTT), hcondADI)
!
endsubroutine ADI_Kprof_1d_mixed
!***********************************************************************
subroutine ADI_Kprof_mixed(f)
!
! 28-fev-2010/dintrans: coded
! simpler version where one part of the radiative diffusion term is
! computed during the explicit step. The implicit part remains
! of Yakonov's form:
!
! (1-dt/2*J_x)*lambda = f_x(T^n) + f_z(T^n) + source
! (1-dt/2*J_y)*beta = lambda
! T^(n+1) = T^n + dt*beta
!
! where J_x and J_y denote Jacobian matrices df/dT.
!
use EquationOfState, only: gamma
use Boundcond, only: update_ghosts
!
implicit none
!
integer :: i,j
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mx,mz) :: source, hcond, dhcond, finter, val, TT, &
chi, dLnhcond
real, dimension(nx) :: ax, bx, cx, wx, rhsx, workx
real, dimension(nz) :: az, bz, cz, wz, rhsz, workz
real :: aalpha, bbeta
logical :: err
character(len=255) :: msg
!
call update_ghosts(f)
!
source=(f(:,4,:,ilnTT)-f(:,4,:,iTTold))/dt
! BC important not for the x-direction (always periodic) but for
! the z-direction as we must impose the 'c3' BC at the 2nd-order
! before going in the implicit stuff
call heatcond_TT(f(:,4,:,iTTold), hcond, dhcond)
call boundary_ADI(f(:,4,:,iTTold), hcond(:,n1))
TT=f(:,4,:,iTTold)
if (ldensity) then
chi=cp1*hcond/exp(f(:,4,:,ilnrho))
! chi=cp1*hcond0/exp(f(:,4,:,ilnrho))
else
chi=cp1*hcond
endif
dLnhcond=dhcond/hcond
! dLnhcond=0.
!
! rows in the x-direction dealt implicitly
!
do j=n1,n2
wx=gamma*chi(l1:l2,j)
ax=-dt/2.*wx*dx_2
bx=1.-dt/2.*wx*dx_2*(-2.+dLnhcond(l1:l2,j)* &
(TT(l1+1:l2+1,j)-2.*TT(l1:l2,j)+TT(l1-1:l2-1,j)))
cx=-dt/2.*wx*dx_2
! rhsx=f_x(T^n) + f_z(T^n) + source
! do first f_z(T^n)
rhsx=wx*dz_2*(TT(l1:l2,j+1)-2.*TT(l1:l2,j)+TT(l1:l2,j-1))
! then add f_x(T^n) + source
rhsx=rhsx+wx*dx_2*(TT(l1+1:l2+1,j)-2.*TT(l1:l2,j)+TT(l1-1:l2-1,j)) &
+source(l1:l2,j)
!
! periodic boundary conditions in x --> cyclic matrix
!
aalpha=cx(nx) ; bbeta=ax(1)
call cyclic(ax,bx,cx,aalpha,bbeta,rhsx,workx,nx)
finter(l1:l2,j)=workx
enddo
!
! columns in the z-direction dealt implicitly
!
do i=l1,l2
wz=dt*gamma*dz_2*chi(i,n1:n2)
az=-wz/2.
bz=1.-wz/2.*(-2.+dLnhcond(i,n1:n2)* &
(TT(i,n1+1:n2+1)-2.*TT(i,n1:n2)+TT(i,n1-1:n2-1)))
cz=-wz/2.
rhsz=finter(i,n1:n2)
!
! z boundary conditions
! Constant temperature at the top: T^(n+1)-T^n=0
bz(nz)=1. ; az(nz)=0.
rhsz(nz)=0.
! bottom
select case (bcz12(ilnTT,1))
! Constant temperature at the bottom: T^(n+1)-T^n=0
case ('cT')
bz(1)=1. ; cz(1)=0.
rhsz(1)=0.
! Constant flux at the bottom
case ('c3')
bz(1)=1. ; cz(1)=-1.
rhsz(1)=0.
case default
call fatal_error('ADI_Kprof_mixed','bcz on TT must be cT or c3')
endselect
!
call tridag(az,bz,cz,rhsz,workz,err,msg)
if (err) call warning('ADI_Kprof_mixed', trim(msg))
val(i,n1:n2)=workz(1:nz)
enddo
!
f(:,4,:,ilnTT)=f(:,4,:,iTTold)+dt*val
!
! update hcond used for the 'c3' condition in boundcond.f90
!
call heatcond_TT(f(:,4,n1,ilnTT), hcondADI)
!
endsubroutine ADI_Kprof_mixed
!***********************************************************************
subroutine ADI_Kconst_yakonov(f)
!
! 26-Jan-2011/dintrans: coded
! 2-D ADI scheme for the radiative diffusion term for a constant
! radiative conductivity K. The ADI scheme is of Yakonov's form:
!
! (1-dt/2*Lamba_x)*T^(n+1/2) = Lambda_x(T^n) + Lambda_z(T^n)
! (1-dt/2*Lamba_z)*T^(n+1) = T^(n+1/2)
!
! where Lambda_x and Lambda_z denote diffusion operators.
! Note: this form is more adapted for a parallelisation compared the
! Peaceman & Rachford one.
!
use EquationOfState, only: gamma, gamma_m1, cs2bot, cs2top
use Boundcond, only: update_ghosts
!
implicit none
!
integer :: i,j
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mx,mz) :: finter, TT, rho1
real, dimension(nx) :: ax, bx, cx, wx, rhsx, workx
real, dimension(nz) :: az, bz, cz, wz, rhsz, workz
real :: aalpha, bbeta
logical :: err
character(len=255) :: msg
!
call update_ghosts(f)
TT=f(:,4,:,iTT)
if (ldensity) then
rho1=exp(-f(:,4,:,ilnrho))
else
rho1=1.
endif
!
! rows dealt implicitly
!
do j=n1,n2
wx=dt*cp1*gamma*hcond0*rho1(l1:l2,j)
ax=-wx*dx_2/2.
bx=1.+wx*dx_2
cx=ax
rhsx=TT(l1:l2,j)+ &
wx*dz_2/2.*(TT(l1:l2,j+1)-2.*TT(l1:l2,j)+TT(l1:l2,j-1))
rhsx=rhsx+wx*dx_2/2.* &
(TT(l1+1:l2+1,j)-2.*TT(l1:l2,j)+TT(l1-1:l2-1,j))
!
! x boundary conditions: periodic
!
aalpha=cx(nx) ; bbeta=ax(1)
call cyclic(ax,bx,cx,aalpha,bbeta,rhsx,workx,nx)
finter(l1:l2,j)=workx
enddo
!
! columns dealt implicitly
!
do i=l1,l2
wz=dt*cp1*gamma*hcond0*dz_2*rho1(i,n1:n2)
az=-wz/2.
bz=1.+wz
cz=az
rhsz=finter(i,n1:n2)
!
! z boundary conditions
!
! Constant temperature at the top
bz(nz)=1. ; az(nz)=0.
rhsz(nz)=cs2top/gamma_m1
! bottom
select case (bcz12(iTT,1))
! Constant temperature at the bottom
case ('cT')
bz(1)=1. ; cz(1)=0.
rhsz(1)=cs2bot/gamma_m1
! Constant flux at the bottom: c1 condition
case ('c1')
bz(1)=1. ; cz(1)=-1.
rhsz(1)=dz*Fbot/hcond0
case default
call fatal_error('ADI_Kconst_yakonov','bcz on TT must be cT or c1')
endselect
!
call tridag(az,bz,cz,rhsz,workz,err,msg)
if (err) call warning('ADI_Kconst_yakonov', trim(msg))
f(i,4,n1:n2,iTT)=workz
enddo
!
endsubroutine ADI_Kconst_yakonov
!***********************************************************************
subroutine crank_Kprof(f)
!
! 01-may-14/dintrans: coded
!
use EquationOfState, only: gamma, gamma_m1, cs2bot, cs2top
use Boundcond, only: update_ghosts
!
implicit none
!
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mz) :: TT, rho1
real, dimension(nz) :: az, bz, cz, rhsz, chi, dchi
logical :: err
character(len=255) :: msg
!
call update_ghosts(f,iTT)
TT=f(4,4,:,iTT)
rho1=exp(-f(4,4,:,ilnrho))
!
chi=0.5*dt*dz_2*gamma*cp1*hcondz(n1:n2)*rho1(n1:n2)
dchi=0.25*dt/dz*gamma*cp1*dhcondz(n1:n2)*rho1(n1:n2)
az=-chi+dchi
bz=1.+2.*chi
cz=-chi-dchi
do n=n1,n2
rhsz(n-nghost)=TT(n)+chi(n-nghost)*(TT(n+1)-2.*TT(n)+TT(n-1)) &
+dchi(n-nghost)*(TT(n+1)-TT(n-1))
enddo
bz(nz)=1. ; az(nz)=0. ; rhsz(nz)=cs2top/gamma_m1 ! T = Ttop
if (bcz12(iTT,1)=='cT') then
bz(1)=1. ; cz(1)=0. ; rhsz(1)=cs2bot/gamma_m1 ! T = Tbot
else
! cz(1)=2.*cz(1) ; rhsz(1)=rhsz(1)+wz(n1)*dz*Fbot/hcondz(n1) ! T' = -Fbot/K
bz(1)=1. ; cz(1)=-1. ; rhsz(1)=dz*Fbot/hcondz(n1) ! T' = -Fbot/K
endif
call tridag(az, bz, cz, rhsz, f(4,4,n1:n2,iTT), err, msg)
if (err) call warning('crank_Kprof', trim(msg))
!
endsubroutine crank_Kprof
!***********************************************************************
subroutine ADI_poly(f)
!
! 01-may-14/dintrans: coded
! 2-D ADI scheme for the radiative diffusion term for a variable
! radiative conductivity K(z). The ADI scheme is of Yakonov's form:
!
! (1-dt/2*Lamba_x)*T^(n+1/2) = T^n + Lambda_x(T^n) + Lambda_z(T^n)
! (1-dt/2*Lamba_z)*T^(n+1) = T^(n+1/2)
!
! where Lambda_x and Lambda_z denote diffusion operators.
!
use EquationOfState, only: gamma, gamma_m1, cs2bot, cs2top
use Boundcond, only: update_ghosts
!
implicit none
!
integer :: i,j
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mx,mz) :: finter, rho1, TT, chi, dchi
real, dimension(nx) :: ax, bx, cx, wx, rhsx, wx1
real, dimension(nz) :: az, bz, cz, wz, rhsz, wz1
real :: aalpha, bbeta
logical :: err
character(len=255) :: msg
!
call update_ghosts(f)
TT=f(:,4,:,iTT)
rho1=exp(-f(:,4,:,ilnrho))
!
! x-direction
!
do j=n1,n2
chi(l1:l2,j)=dt*cp1*gamma*hcondz(j)*rho1(l1:l2,j)
dchi(l1:l2,j)=dt*cp1*gamma*dhcondz(j)*rho1(l1:l2,j)
wx=0.5*chi(l1:l2,j)
wx1=0.25*dchi(l1:l2,j)
ax=-wx*dx_2
bx=1.+2.*wx*dx_2
cx=ax
rhsx=TT(l1:l2,j) &
+wx*dx_2*(TT(l1+1:l2+1,j)-2.*TT(l1:l2,j)+TT(l1-1:l2-1,j)) &
+wx*dz_2*(TT(l1:l2,j+1)-2.*TT(l1:l2,j)+TT(l1:l2,j-1)) &
+wx1/dz*(TT(l1:l2,j+1)-TT(l1:l2,j-1))
!
! x boundary conditions: periodic
!
aalpha=cx(nx) ; bbeta=ax(1)
call cyclic(ax,bx,cx,aalpha,bbeta,rhsx,finter(l1:l2,j),nx)
enddo
!
! z-direction
!
do i=l1,l2
wz=0.5*dz_2*chi(i,n1:n2)
wz1=0.25/dz*dchi(i,n1:n2)
az=-wz+wz1
bz=1.+2.*wz
cz=-wz-wz1
rhsz=finter(i,n1:n2)
!
! z boundary conditions
!
! Constant temperature at the top
bz(nz)=1. ; az(nz)=0. ; rhsz(nz)=cs2top/gamma_m1
! bottom
select case (bcz12(iTT,1))
! Constant temperature at the bottom
case ('cT')
bz(1)=1. ; cz(1)=0. ; rhsz(1)=cs2bot/gamma_m1
! Constant flux at the bottom: c1 condition
case ('c1')
bz(1)=1. ; cz(1)=-1 ; rhsz(1)=dz*Fbot/hcondz(n1)
case default
call fatal_error('ADI_poly','bcz on TT must be cT or c1')
endselect
!
call tridag(az,bz,cz,rhsz,f(i,4,n1:n2,iTT),err,msg)
if (err) call warning('ADI_poly', trim(msg))
enddo
!
endsubroutine ADI_poly
!***********************************************************************
subroutine ADI_poly_MPI(f)
!
! 01-may-14/dintrans: coded
! Parallel version in the z-direction of ADI_poly.
!
use EquationOfState, only: gamma, gamma_m1, cs2bot, cs2top
use Boundcond, only: update_ghosts
use Mpicomm, only: transp_xz, transp_zx
!
implicit none
!
integer :: i,j
integer, parameter :: nxt=nx/nprocz
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mx,mz) :: finter, rho1, TT, chi, dchi
real, dimension(nx) :: ax, bx, cx, wx, rhsx, wx1
real, dimension(nzgrid) :: az, bz, cz, wz, rhsz, wz1
real, dimension(nzgrid,nxt) :: fintert, chit, dchit, wtmp
real :: aalpha, bbeta
logical :: err
character(len=255) :: msg
!
call update_ghosts(f,iTT)
TT=f(:,4,:,iTT)
rho1=exp(-f(:,4,:,ilnrho))
!
! x-direction
!
do j=n1,n2
chi(l1:l2,j)=dt*cp1*gamma*hcondz(j)*rho1(l1:l2,j)
dchi(l1:l2,j)=dt*cp1*gamma*dhcondz(j)*rho1(l1:l2,j)
wx=0.5*chi(l1:l2,j)
wx1=0.25*dchi(l1:l2,j)
ax=-wx*dx_2
bx=1.+2.*wx*dx_2
cx=ax
rhsx=TT(l1:l2,j) &
+wx*dx_2*(TT(l1+1:l2+1,j)-2.*TT(l1:l2,j)+TT(l1-1:l2-1,j)) &
+wx*dz_2*(TT(l1:l2,j+1)-2.*TT(l1:l2,j)+TT(l1:l2,j-1)) &
+wx1/dz*(TT(l1:l2,j+1)-TT(l1:l2,j-1))
!
! x boundary conditions: periodic
!
aalpha=cx(nx) ; bbeta=ax(1)
call cyclic(ax,bx,cx,aalpha,bbeta,rhsx,finter(l1:l2,j),nx)
enddo
!
! do the transpositions x <--> z
!
call transp_xz(finter(l1:l2,n1:n2), fintert)
call transp_xz(chi(l1:l2,n1:n2), chit)
call transp_xz(dchi(l1:l2,n1:n2), dchit)
!
! z-direction
!
do i=1,nxt
wz=0.5*dz_2*chit(:,i)
wz1=0.25/dz*dchit(:,i)
az=-wz+wz1
bz=1.+2.*wz
cz=-wz-wz1
rhsz=fintert(:,i)
!
! z boundary conditions
!
! Constant temperature at the top
bz(nzgrid)=1. ; az(nzgrid)=0. ; rhsz(nzgrid)=cs2top/gamma_m1
! bottom
select case (bcz12(iTT,1))
! Constant temperature at the bottom
case ('cT')
bz(1)=1. ; cz(1)=0. ; rhsz(1)=cs2bot/gamma_m1
! Constant flux at the bottom: c1 condition
case ('c1')
! bz(1)=1. ; cz(1)=-1 ; rhsz(1)=dz*Fbot/hcondz(n1)
bz(1)=1. ; cz(1)=-1. ; rhsz(1)=dz*Fbot/hcond0
case default
call fatal_error('ADI_poly_MPI','bcz on TT must be cT or c1')
endselect
!
call tridag(az,bz,cz,rhsz,wtmp(:,i),err,msg)
if (err) call warning('ADI_poly_MPI', trim(msg))
enddo
call transp_zx(wtmp, f(l1:l2,4,n1:n2,iTT))
!
endsubroutine ADI_poly_MPI
!***********************************************************************
subroutine ADI3D(f)
!
! 02-may-14/dintrans: coded
! 3-D ADI scheme for the radiative diffusion term for a variable
! radiative conductivity K(z). The ADI scheme is of Yakonov's form:
!
! (1-dt/2*L_x)*T^(n+1/3) = T^n + L_x(T^n) + + L_y(T^n) + L_z(T^n) + source
! (1-dt/2*L_y)*T^(n+2/3) = T^(n+1/3)
! (1-dt/2*L_z)*T^(n+1) = T^(n+2/3)
!
! where L_x, L_y and L_z denote diffusion operators and the source
! term comes from the explicit advance.
!
use EquationOfState, only: gamma, gamma_m1, cs2bot, cs2top
use Boundcond, only: update_ghosts
!
implicit none
!
integer :: l,m,n
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mx,my,mz) :: source, finterx, fintery, TT, chi, dchi
real, dimension(nx) :: ax, bx, cx, rhsx, wx, wx1
real, dimension(ny) :: ay, by, cy, rhsy, wy, wy1
real, dimension(nz) :: az, bz, cz, rhsz, wz, wz1
real :: aalpha, bbeta
logical :: err
character(len=255) :: msg
!
call update_ghosts(f)
!
source=(f(:,:,:,ilnTT)-f(:,:,:,iTTold))/dt
TT=f(:,:,:,iTTold)
!
! x-direction
!
do n=n1,n2
do m=m1,m2
chi(l1:l2,m,n)=dt*cp1*gamma*hcondz(n)/exp(f(l1:l2,m,n,ilnrho))
dchi(l1:l2,m,n)=dt*cp1*gamma*dhcondz(n)/exp(f(l1:l2,m,n,ilnrho))
wx=0.5*chi(l1:l2,m,n)
wx1=0.25*dchi(l1:l2,m,n)
ax=-wx*dx_2
bx=1.+2.*wx*dx_2
cx=ax
rhsx=TT(l1:l2,m,n)+ &
wx*dz_2*(TT(l1:l2,m,n+1)-2.*TT(l1:l2,m,n)+TT(l1:l2,m,n-1)) &
+wx1/dz*(TT(l1:l2,m,n+1)-TT(l1:l2,m,n-1)) &
+wx*dy_2*(TT(l1:l2,m+1,n)-2.*TT(l1:l2,m,n)+TT(l1:l2,m-1,n))
rhsx=rhsx+wx*dx_2* &
(TT(l1+1:l2+1,m,n)-2.*TT(l1:l2,m,n)+TT(l1-1:l2-1,m,n)) &
+dt*source(l1:l2,m,n)
!
! x boundary conditions: periodic
!
aalpha=cx(nx) ; bbeta=ax(1)
call cyclic(ax,bx,cx,aalpha,bbeta,rhsx,finterx(l1:l2,m,n),nx)
enddo
enddo
!
! y-direction
!
do n=n1,n2
do l=l1,l2
wy=0.5*chi(l,m1:m2,n)
wy1=0.25*dchi(l,m1:m2,n)
ay=-wy*dy_2
by=1.+2.*wy*dy_2
cy=ay
rhsy=finterx(l,m1:m2,n)
!
! y boundary conditions: periodic
!
aalpha=cy(ny) ; bbeta=ay(1)
call cyclic(ay,by,cy,aalpha,bbeta,rhsy,fintery(l,m1:m2,n),ny)
enddo
enddo
!
! z-direction
!
do m=m1,m2
do l=l1,l2
wz=0.5*dz_2*chi(l,m,n1:n2)
wz1=0.25/dz*dchi(l,m,n1:n2)
az=-wz+wz1
bz=1.+2.*wz
cz=-wz-wz1
rhsz=fintery(l,m,n1:n2)
!
! z boundary conditions
!
! Constant temperature at the top
bz(nz)=1. ; az(nz)=0. ; rhsz(nz)=cs2top/gamma_m1
! bottom
select case (bcz12(ilnTT,1))
! Constant temperature at the bottom
case ('cT')
bz(1)=1. ; cz(1)=0. ; rhsz(1)=cs2bot/gamma_m1
! Constant flux at the bottom: c1 condition
case ('c1')
bz(1)=1. ; cz(1)=-1 ; rhsz(1)=dz*Fbot/hcondz(n1)
case default
call fatal_error('ADI_poly','bcz on TT must be cT or c1')
endselect
!
call tridag(az,bz,cz,rhsz,f(l,m,n1:n2,ilnTT),err,msg)
if (err) call warning('ADI_poly', trim(msg))
enddo
enddo
!
endsubroutine ADI3D
!***********************************************************************
subroutine ADI3D_MPI(f)
!
! 02-may-14/dintrans: coded
! 3-D ADI scheme for the radiative diffusion term for a variable
! radiative conductivity K(z). The ADI scheme is of Yakonov's form:
!
! (1-dt/2*L_x)*T^(n+1/3) = T^n + L_x(T^n) + + L_y(T^n) + L_z(T^n) + source
! (1-dt/2*L_y)*T^(n+2/3) = T^(n+1/3)
! (1-dt/2*L_z)*T^(n+1) = T^(n+2/3)
!
! where L_x, L_y and L_z denote diffusion operators and the source
! term comes from the explicit advance.
!
use EquationOfState, only: gamma, gamma_m1, cs2bot, cs2top
use Boundcond, only: update_ghosts
use Mpicomm, only: transp_xz, transp_zx
!
implicit none
!
integer, parameter :: nxt=nx/nprocz
integer :: l,m,n
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(mx,my,mz) :: source, finterx, fintery, TT, chi, dchi
real, dimension(nx) :: ax, bx, cx, rhsx, wx, wx1
real, dimension(ny) :: ay, by, cy, rhsy, wy, wy1
real, dimension(nzgrid) :: az, bz, cz, rhsz, wz, wz1
real, dimension(nzgrid,nxt) :: finteryt, chit, dchit, wtmp
real :: aalpha, bbeta
logical :: err
character(len=255) :: msg
!
call update_ghosts(f)
!
source=(f(:,:,:,ilnTT)-f(:,:,:,iTTold))/dt
TT=f(:,:,:,iTTold)
!
! x-direction
!
do n=n1,n2
do m=m1,m2
chi(l1:l2,m,n)=dt*cp1*gamma*hcondz(n)/exp(f(l1:l2,m,n,ilnrho))
dchi(l1:l2,m,n)=dt*cp1*gamma*dhcondz(n)/exp(f(l1:l2,m,n,ilnrho))
wx=0.5*chi(l1:l2,m,n)
wx1=0.25*dchi(l1:l2,m,n)
ax=-wx*dx_2
bx=1.+2.*wx*dx_2
cx=ax
rhsx=TT(l1:l2,m,n)+ &
wx*dz_2*(TT(l1:l2,m,n+1)-2.*TT(l1:l2,m,n)+TT(l1:l2,m,n-1)) &
+wx1/dz*(TT(l1:l2,m,n+1)-TT(l1:l2,m,n-1)) &
+wx*dy_2*(TT(l1:l2,m+1,n)-2.*TT(l1:l2,m,n)+TT(l1:l2,m-1,n))
rhsx=rhsx+wx*dx_2* &
(TT(l1+1:l2+1,m,n)-2.*TT(l1:l2,m,n)+TT(l1-1:l2-1,m,n)) &
+dt*source(l1:l2,m,n)
!
! x boundary conditions: periodic
!
aalpha=cx(nx) ; bbeta=ax(1)
call cyclic(ax,bx,cx,aalpha,bbeta,rhsx,finterx(l1:l2,m,n),nx)
enddo
enddo
!
! y-direction
!
do n=n1,n2
do l=l1,l2
wy=0.5*chi(l,m1:m2,n)
wy1=0.25*dchi(l,m1:m2,n)
ay=-wy*dy_2
by=1.+2.*wy*dy_2
cy=ay
rhsy=finterx(l,m1:m2,n)
!
! y boundary conditions: periodic
!
aalpha=cy(ny) ; bbeta=ay(1)
call cyclic(ay,by,cy,aalpha,bbeta,rhsy,fintery(l,m1:m2,n),ny)
enddo
enddo
!
! z-direction
!
do m=m1,m2
call transp_xz(fintery(l1:l2,m,n1:n2), finteryt)
call transp_xz(chi(l1:l2,m,n1:n2), chit)
call transp_xz(dchi(l1:l2,m,n1:n2), dchit)
do l=1,nxt
wz=0.5*dz_2*chit(:,l)
wz1=0.25/dz*dchit(:,l)
az=-wz+wz1
bz=1.+2.*wz
cz=-wz-wz1
rhsz=finteryt(:,l)
!
! z boundary conditions
!
! Constant temperature at the top
bz(nzgrid)=1. ; az(nzgrid)=0. ; rhsz(nzgrid)=cs2top/gamma_m1
! bottom
select case (bcz12(ilnTT,1))
! Constant temperature at the bottom
case ('cT')
bz(1)=1. ; cz(1)=0. ; rhsz(1)=cs2bot/gamma_m1
! Constant flux at the bottom: c1 condition
case ('c1')
! bz(1)=1. ; cz(1)=-1 ; rhsz(1)=dz*Fbot/hcondz(n1)
bz(1)=1. ; cz(1)=-1 ; rhsz(1)=dz*Fbot/hcond0
case default
call fatal_error('ADI_poly','bcz on TT must be cT or c1')
endselect
!
call tridag(az,bz,cz,rhsz,wtmp(:,l),err,msg)
if (err) call warning('ADI_poly', trim(msg))
enddo
call transp_zx(wtmp, f(l1:l2,m,n1:n2,ilnTT))
enddo
!
endsubroutine ADI3D_MPI
!***********************************************************************
endmodule ImplicitPhysics