Grzegorz Kowal
3cd24e3b0f
In the sum of the J² we used Fortran subroutine sum() and the sum was done along the 2nd index. However, the interpretation of this index is uncertain when the first array rank is 1. Therefore, instead of using subroutine sum(), the calculation is done directly. After this change there are no strange effects appearing on the boundaries of the blocks. Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
620 lines
18 KiB
Fortran
620 lines
18 KiB
Fortran
!!******************************************************************************
|
||
!!
|
||
!! This file is part of the AMUN source code, a program to perform
|
||
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
||
!! adaptive mesh.
|
||
!!
|
||
!! Copyright (C) 2008-2017 Grzegorz Kowal <grzegorz@amuncode.org>
|
||
!!
|
||
!! This program is free software: you can redistribute it and/or modify
|
||
!! it under the terms of the GNU General Public License as published by
|
||
!! the Free Software Foundation, either version 3 of the License, or
|
||
!! (at your option) any later version.
|
||
!!
|
||
!! This program is distributed in the hope that it will be useful,
|
||
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||
!! GNU General Public License for more details.
|
||
!!
|
||
!! You should have received a copy of the GNU General Public License
|
||
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||
!!
|
||
!!******************************************************************************
|
||
!!
|
||
!! module: SOURCES
|
||
!!
|
||
!! This modules adds source terms.
|
||
!!
|
||
!!******************************************************************************
|
||
!
|
||
module sources
|
||
|
||
#ifdef PROFILE
|
||
! include external procedures
|
||
!
|
||
use timers, only : set_timer, start_timer, stop_timer
|
||
#endif /* PROFILE */
|
||
|
||
! module variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
#ifdef PROFILE
|
||
! timer indices
|
||
!
|
||
integer, save :: imi, imu
|
||
#endif /* PROFILE */
|
||
|
||
! GLM-MHD source terms type (1 - EGLM, 2 - HEGLM)
|
||
!
|
||
integer , save :: glm_type = 0
|
||
|
||
! viscosity coefficient
|
||
!
|
||
real(kind=8), save :: viscosity = 0.0d+00
|
||
|
||
! resistivity coefficient
|
||
!
|
||
real(kind=8), save :: resistivity = 0.0d+00
|
||
|
||
! by default everything is private
|
||
!
|
||
private
|
||
|
||
! declare public subroutines
|
||
!
|
||
public :: initialize_sources, finalize_sources
|
||
public :: update_sources
|
||
public :: viscosity, resistivity
|
||
|
||
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
||
!
|
||
contains
|
||
!
|
||
!===============================================================================
|
||
!!
|
||
!!*** PUBLIC SUBROUTINES *****************************************************
|
||
!!
|
||
!===============================================================================
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine INITIALIZE_SOURCES:
|
||
! -----------------------------
|
||
!
|
||
! Subroutine initializes module SOURCES.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! verbose - a logical flag turning the information printing;
|
||
! iret - an integer flag for error return value;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine initialize_sources(verbose, iret)
|
||
|
||
! include external procedures and variables
|
||
!
|
||
use parameters , only : get_parameter_string, get_parameter_real
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
logical, intent(in) :: verbose
|
||
integer, intent(inout) :: iret
|
||
|
||
! local variables
|
||
!
|
||
character(len=8) :: tglm = "none"
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
#ifdef PROFILE
|
||
! set timer descriptions
|
||
!
|
||
call set_timer('sources:: initialize', imi)
|
||
call set_timer('sources:: update' , imu)
|
||
|
||
! start accounting time for module initialization/finalization
|
||
!
|
||
call start_timer(imi)
|
||
#endif /* PROFILE */
|
||
|
||
! get the type of the GLM source terms
|
||
!
|
||
call get_parameter_string("glm_source_terms", tglm)
|
||
|
||
! set the glm_type variable to correct value
|
||
!
|
||
select case(trim(tglm))
|
||
case("eglm", "EGLM")
|
||
glm_type = 1
|
||
|
||
case("heglm", "HEGLM")
|
||
glm_type = 2
|
||
|
||
case default
|
||
glm_type = 0
|
||
end select
|
||
|
||
! get viscosity coefficient
|
||
!
|
||
call get_parameter_real("viscosity" , viscosity)
|
||
|
||
! get resistivity coefficient
|
||
!
|
||
call get_parameter_real("resistivity", resistivity)
|
||
|
||
! print information about the Riemann solver
|
||
!
|
||
if (verbose) then
|
||
|
||
write (*,"(4x,a,1x,a) ") "glm source terms =", trim(tglm)
|
||
write (*,"(4x,a,1x,1e9.2)") "viscosity =", viscosity
|
||
write (*,"(4x,a,1x,1e9.2)") "resistivity =", resistivity
|
||
|
||
end if
|
||
|
||
#ifdef PROFILE
|
||
! stop accounting time for module initialization/finalization
|
||
!
|
||
call stop_timer(imi)
|
||
#endif /* PROFILE */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine initialize_sources
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine FINALIZE_SOURCES:
|
||
! ---------------------------
|
||
!
|
||
! Subroutine releases memory used by the module.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! iret - an integer flag for error return value;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine finalize_sources(iret)
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
integer, intent(inout) :: iret
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
#ifdef PROFILE
|
||
! start accounting time for module initialization/finalization
|
||
!
|
||
call start_timer(imi)
|
||
#endif /* PROFILE */
|
||
|
||
#ifdef PROFILE
|
||
! stop accounting time for module initialization/finalization
|
||
!
|
||
call stop_timer(imi)
|
||
#endif /* PROFILE */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine finalize_sources
|
||
!
|
||
!===============================================================================
|
||
!
|
||
! subroutine UPDATE_SOURCES:
|
||
! -------------------------
|
||
!
|
||
! Subroutine add the source terms.
|
||
!
|
||
! Arguments:
|
||
!
|
||
! pdata - the pointer to a data block;
|
||
! t, dt - the time and time increment;
|
||
! du - the array of variable increment;
|
||
!
|
||
!===============================================================================
|
||
!
|
||
subroutine update_sources(pdata, t, dt, du)
|
||
|
||
! include external variables
|
||
!
|
||
use blocks , only : block_data
|
||
use coordinates , only : im, jm, km
|
||
use coordinates , only : ax, ay, az, adx, ady, adz
|
||
use equations , only : nv, inx, iny, inz
|
||
use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ien
|
||
use equations , only : ibx, iby, ibz, ibp
|
||
use gravity , only : gravity_enabled, gravitational_acceleration
|
||
use operators , only : divergence, gradient, laplace, curl
|
||
use user_problem , only : update_sources_user
|
||
|
||
! local variables are not implicit by default
|
||
!
|
||
implicit none
|
||
|
||
! subroutine arguments
|
||
!
|
||
type(block_data), pointer , intent(inout) :: pdata
|
||
real(kind=8) , intent(in) :: t, dt
|
||
real(kind=8), dimension(nv,im,jm,km), intent(inout) :: du
|
||
|
||
! local variables
|
||
!
|
||
integer :: i , j , k
|
||
real(kind=8) :: fc, gc
|
||
real(kind=8) :: gx, gy, gz
|
||
real(kind=8) :: dbx, dby, dbz
|
||
real(kind=8) :: dvxdx, dvxdy, dvxdz, divv
|
||
real(kind=8) :: dvydx, dvydy, dvydz
|
||
real(kind=8) :: dvzdx, dvzdy, dvzdz
|
||
|
||
! local arrays
|
||
!
|
||
real(kind=8), dimension(3) :: ga, dh
|
||
real(kind=8), dimension(im) :: x
|
||
real(kind=8), dimension(jm) :: y
|
||
real(kind=8), dimension(km) :: z
|
||
real(kind=8), dimension(im,jm,km) :: db
|
||
real(kind=8), dimension(3,3,im,jm,km) :: tmp
|
||
!
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
#ifdef PROFILE
|
||
! start accounting time for source terms
|
||
!
|
||
call start_timer(imu)
|
||
#endif /* PROFILE */
|
||
|
||
! proceed only if the gravitational term is enabled
|
||
!
|
||
if (gravity_enabled) then
|
||
|
||
! prepare block coordinates
|
||
!
|
||
x(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im)
|
||
y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm)
|
||
#if NDIMS == 3
|
||
z(1:km) = pdata%meta%zmin + az(pdata%meta%level,1:km)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! iterate over all positions in the YZ plane
|
||
!
|
||
do k = 1, km
|
||
do j = 1, jm
|
||
do i = 1, im
|
||
|
||
! get gravitational acceleration components
|
||
!
|
||
call gravitational_acceleration(t, dt, x(i), y(j), z(k), ga(1:3))
|
||
|
||
! calculate the gravitational source terms
|
||
!
|
||
gx = pdata%q(idn,i,j,k) * ga(1)
|
||
gy = pdata%q(idn,i,j,k) * ga(2)
|
||
#if NDIMS == 3
|
||
gz = pdata%q(idn,i,j,k) * ga(3)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! add source terms to momentum equations
|
||
!
|
||
du(imx,i,j,k) = du(imx,i,j,k) + gx
|
||
du(imy,i,j,k) = du(imy,i,j,k) + gy
|
||
#if NDIMS == 3
|
||
du(imz,i,j,k) = du(imz,i,j,k) + gz
|
||
#endif /* NDIMS == 3 */
|
||
|
||
! add source terms to total energy equation
|
||
!
|
||
if (ien > 0) then
|
||
|
||
#if NDIMS == 2
|
||
du(ien,i,j,k) = du(ien,i,j,k) + gx * pdata%q(ivx,i,j,k) &
|
||
+ gy * pdata%q(ivy,i,j,k)
|
||
#endif /* NDIMS == 2 */
|
||
#if NDIMS == 3
|
||
du(ien,i,j,k) = du(ien,i,j,k) + gx * pdata%q(ivx,i,j,k) &
|
||
+ gy * pdata%q(ivy,i,j,k) &
|
||
+ gz * pdata%q(ivz,i,j,k)
|
||
#endif /* NDIMS == 3 */
|
||
|
||
end if
|
||
|
||
end do ! i = 1, im
|
||
end do ! j = 1, jm
|
||
end do ! k = 1, km
|
||
|
||
end if ! gravity enabled
|
||
|
||
! proceed only if the viscosity coefficient is not zero
|
||
!
|
||
if (viscosity > 0.0d+00) then
|
||
|
||
! prepare coordinate increments
|
||
!
|
||
dh(1) = adx(pdata%meta%level)
|
||
dh(2) = ady(pdata%meta%level)
|
||
dh(3) = adz(pdata%meta%level)
|
||
|
||
! calculate the velocity Jacobian
|
||
!
|
||
call gradient(dh(:), pdata%q(ivx,1:im,1:jm,1:km) &
|
||
, tmp(inx,inx:inz,1:im,1:jm,1:km))
|
||
call gradient(dh(:), pdata%q(ivy,1:im,1:jm,1:km) &
|
||
, tmp(iny,inx:inz,1:im,1:jm,1:km))
|
||
call gradient(dh(:), pdata%q(ivz,1:im,1:jm,1:km) &
|
||
, tmp(inz,inx:inz,1:im,1:jm,1:km))
|
||
|
||
! iterate over all cells
|
||
!
|
||
do k = 1, km
|
||
do j = 1, jm
|
||
do i = 1, im
|
||
|
||
! prepare the νρ factor
|
||
!
|
||
gc = viscosity * pdata%q(idn,i,j,k)
|
||
fc = 2.0d+00 * gc
|
||
|
||
! get the velocity Jacobian elements
|
||
!
|
||
dvxdx = tmp(inx,inx,i,j,k)
|
||
dvxdy = tmp(inx,iny,i,j,k)
|
||
dvxdz = tmp(inx,inz,i,j,k)
|
||
dvydx = tmp(iny,inx,i,j,k)
|
||
dvydy = tmp(iny,iny,i,j,k)
|
||
dvydz = tmp(iny,inz,i,j,k)
|
||
dvzdx = tmp(inz,inx,i,j,k)
|
||
dvzdy = tmp(inz,iny,i,j,k)
|
||
dvzdz = tmp(inz,inz,i,j,k)
|
||
divv = (dvxdx + dvydy + dvzdz) / 3.0d+00
|
||
|
||
! calculate elements of the viscous stress tensor
|
||
!
|
||
tmp(inx,inx,i,j,k) = fc * (dvxdx - divv)
|
||
tmp(iny,iny,i,j,k) = fc * (dvydy - divv)
|
||
tmp(inz,inz,i,j,k) = fc * (dvzdz - divv)
|
||
tmp(inx,iny,i,j,k) = gc * (dvxdy + dvydx)
|
||
tmp(inx,inz,i,j,k) = gc * (dvxdz + dvzdx)
|
||
tmp(iny,inz,i,j,k) = gc * (dvydz + dvzdy)
|
||
tmp(iny,inx,i,j,k) = tmp(inx,iny,i,j,k)
|
||
tmp(inz,inx,i,j,k) = tmp(inx,inz,i,j,k)
|
||
tmp(inz,iny,i,j,k) = tmp(iny,inz,i,j,k)
|
||
|
||
end do ! i = 1, im
|
||
end do ! j = 1, jm
|
||
end do ! k = 1, km
|
||
|
||
! calculate the divergence of the first tensor row
|
||
!
|
||
call divergence(dh(:), tmp(inx,inx:inz,1:im,1:jm,1:km) &
|
||
, db(1:im,1:jm,1:km))
|
||
|
||
! add viscous source terms to the X momentum equation
|
||
!
|
||
du(imx,1:im,1:jm,1:km) = du(imx,1:im,1:jm,1:km) + db(1:im,1:jm,1:km)
|
||
|
||
! calculate the divergence of the second tensor row
|
||
!
|
||
call divergence(dh(:), tmp(iny,inx:inz,1:im,1:jm,1:km) &
|
||
, db(1:im,1:jm,1:km))
|
||
|
||
! add viscous source terms to the Y momentum equation
|
||
!
|
||
du(imy,1:im,1:jm,1:km) = du(imy,1:im,1:jm,1:km) + db(1:im,1:jm,1:km)
|
||
|
||
! calculate the divergence of the third tensor row
|
||
!
|
||
call divergence(dh(:), tmp(inz,inx:inz,1:im,1:jm,1:km) &
|
||
, db(1:im,1:jm,1:km))
|
||
|
||
! add viscous source terms to the Z momentum equation
|
||
!
|
||
du(imz,1:im,1:jm,1:km) = du(imz,1:im,1:jm,1:km) + db(1:im,1:jm,1:km)
|
||
|
||
! add viscous source term to total energy equation
|
||
!
|
||
if (ien > 0) then
|
||
|
||
! iterate over all cells
|
||
!
|
||
do k = 1, km
|
||
do j = 1, jm
|
||
do i = 1, im
|
||
|
||
! calculate scalar product of v and viscous stress tensor τ
|
||
!
|
||
gx = pdata%q(ivx,i,j,k) * tmp(inx,inx,i,j,k) &
|
||
+ pdata%q(ivy,i,j,k) * tmp(inx,iny,i,j,k) &
|
||
+ pdata%q(ivz,i,j,k) * tmp(inx,inz,i,j,k)
|
||
gy = pdata%q(ivx,i,j,k) * tmp(iny,inx,i,j,k) &
|
||
+ pdata%q(ivy,i,j,k) * tmp(iny,iny,i,j,k) &
|
||
+ pdata%q(ivz,i,j,k) * tmp(iny,inz,i,j,k)
|
||
gz = pdata%q(ivx,i,j,k) * tmp(inz,inx,i,j,k) &
|
||
+ pdata%q(ivy,i,j,k) * tmp(inz,iny,i,j,k) &
|
||
+ pdata%q(ivz,i,j,k) * tmp(inz,inz,i,j,k)
|
||
|
||
! update (v.τ), use the first row of the tensor tmp
|
||
!
|
||
tmp(inx,inx,i,j,k) = gx
|
||
tmp(inx,iny,i,j,k) = gy
|
||
tmp(inx,inz,i,j,k) = gz
|
||
|
||
end do ! i = 1, im
|
||
end do ! j = 1, jm
|
||
end do ! k = 1, km
|
||
|
||
! calculate the divergence of (v.τ)
|
||
!
|
||
call divergence(dh(:), tmp(inx,inx:inz,1:im,1:jm,1:km) &
|
||
, db(1:im,1:jm,1:km))
|
||
|
||
! update the energy increment
|
||
!
|
||
du(ien,1:im,1:jm,1:km) = du(ien,1:im,1:jm,1:km) + db(1:im,1:jm,1:km)
|
||
|
||
end if ! ien > 0
|
||
|
||
end if ! viscosity is not zero
|
||
|
||
!=== add magnetic field related source terms ===
|
||
!
|
||
if (ibx > 0) then
|
||
|
||
! prepare coordinate increments
|
||
!
|
||
dh(1) = adx(pdata%meta%level)
|
||
dh(2) = ady(pdata%meta%level)
|
||
dh(3) = adz(pdata%meta%level)
|
||
|
||
! add the EGLM-MHD source terms
|
||
!
|
||
if (glm_type > 0) then
|
||
|
||
! calculate the magnetic field divergence
|
||
!
|
||
call divergence(dh(:), pdata%q(ibx:ibz,:,:,:), db(:,:,:))
|
||
|
||
! update the momentum component increments, i.e.
|
||
! d/dt (ρv) + ∇.F = - (∇.B)B
|
||
!
|
||
du(imx,:,:,:) = du(imx,:,:,:) - db(:,:,:) * pdata%q(ibx,:,:,:)
|
||
du(imy,:,:,:) = du(imy,:,:,:) - db(:,:,:) * pdata%q(iby,:,:,:)
|
||
du(imz,:,:,:) = du(imz,:,:,:) - db(:,:,:) * pdata%q(ibz,:,:,:)
|
||
|
||
! update the energy equation
|
||
!
|
||
if (ien > 0 .and. ibp > 0) then
|
||
|
||
! calculate the gradient of divergence potential
|
||
!
|
||
call gradient(dh(:), pdata%q(ibp,:,:,:), tmp(inx:inz,inx,:,:,:))
|
||
|
||
! add the divergence potential source term to the energy equation, i.e.
|
||
! d/dt E + ∇.F = - B.(∇ψ)
|
||
!
|
||
du(ien,:,:,:) = du(ien,:,:,:) &
|
||
- sum(pdata%q(ibx:ibz,:,:,:) * tmp(inx:inz,inx,:,:,:), 1)
|
||
|
||
end if ! ien > 0
|
||
|
||
! add the HEGLM-MHD source terms
|
||
!
|
||
if (glm_type > 1) then
|
||
|
||
! update magnetic field component increments, i.e.
|
||
! d/dt B + ∇.F = - (∇.B)v
|
||
!
|
||
du(ibx,:,:,:) = du(ibx,:,:,:) - db(:,:,:) * pdata%q(ivx,:,:,:)
|
||
du(iby,:,:,:) = du(iby,:,:,:) - db(:,:,:) * pdata%q(ivy,:,:,:)
|
||
du(ibz,:,:,:) = du(ibz,:,:,:) - db(:,:,:) * pdata%q(ivz,:,:,:)
|
||
|
||
! update the energy equation
|
||
!
|
||
if (ien > 0) then
|
||
|
||
! calculate scalar product of velocity and magnetic field
|
||
!
|
||
tmp(inx,inx,:,:,:) = sum(pdata%q(ivx:ivz,:,:,:) &
|
||
* pdata%q(ibx:ibz,:,:,:), 1)
|
||
|
||
! add the divergence potential source term to the energy equation, i.e.
|
||
! d/dt E + ∇.F = - (∇.B) (v.B)
|
||
!
|
||
du(ien,:,:,:) = du(ien,:,:,:) - db(:,:,:) * tmp(inx,inx,:,:,:)
|
||
|
||
end if ! ien > 0
|
||
|
||
end if ! glm_type > 1
|
||
|
||
end if ! glmtype > 0
|
||
|
||
! proceed only if the resistivity coefficient is not zero
|
||
!
|
||
if (resistivity > 0.0d+00) then
|
||
|
||
! calculate the Laplace operator of B, i.e. Δ(B)
|
||
!
|
||
call laplace(dh(:), pdata%q(ibx,1:im,1:jm,1:km) &
|
||
, tmp(inx,inx,1:im,1:jm,1:km))
|
||
call laplace(dh(:), pdata%q(iby,1:im,1:jm,1:km) &
|
||
, tmp(inx,iny,1:im,1:jm,1:km))
|
||
call laplace(dh(:), pdata%q(ibz,1:im,1:jm,1:km) &
|
||
, tmp(inx,inz,1:im,1:jm,1:km))
|
||
|
||
! multiply by the resistivity coefficient
|
||
!
|
||
tmp(iny,inx:inz,1:im,1:jm,1:km) = &
|
||
resistivity * tmp(inx,inx:inz,1:im,1:jm,1:km)
|
||
|
||
! update magnetic field component increments
|
||
!
|
||
du(ibx,1:im,1:jm,1:km) = du(ibx,1:im,1:jm,1:km) &
|
||
+ tmp(iny,inx,1:im,1:jm,1:km)
|
||
du(iby,1:im,1:jm,1:km) = du(iby,1:im,1:jm,1:km) &
|
||
+ tmp(iny,iny,1:im,1:jm,1:km)
|
||
du(ibz,1:im,1:jm,1:km) = du(ibz,1:im,1:jm,1:km) &
|
||
+ tmp(iny,inz,1:im,1:jm,1:km)
|
||
|
||
! update energy equation
|
||
!
|
||
if (ien > 0) then
|
||
|
||
! add the first resistive source term to the energy equation, i.e.
|
||
! d/dt E + ∇.F = η B.[Δ(B)]
|
||
!
|
||
du(ien,1:im,1:jm,1:km) = du(ien,1:im,1:jm,1:km) &
|
||
+ (pdata%q(ibx,1:im,1:jm,1:km) * tmp(iny,inx,1:im,1:jm,1:km) &
|
||
+ pdata%q(iby,1:im,1:jm,1:km) * tmp(iny,iny,1:im,1:jm,1:km) &
|
||
+ pdata%q(ibz,1:im,1:jm,1:km) * tmp(iny,inz,1:im,1:jm,1:km))
|
||
|
||
! calculate current density J = ∇xB
|
||
!
|
||
call curl(dh(:), pdata%q(ibx:ibz,1:im,1:jm,1:km) &
|
||
, tmp(inz,inx:inz,1:im,1:jm,1:km))
|
||
|
||
! calculate J²
|
||
!
|
||
db(1:im,1:jm,1:km) = tmp(inz,inx,1:im,1:jm,1:km)**2 &
|
||
+ tmp(inz,iny,1:im,1:jm,1:km)**2 &
|
||
+ tmp(inz,inz,1:im,1:jm,1:km)**2
|
||
|
||
! add the second resistive source term to the energy equation, i.e.
|
||
! d/dt E + ∇.F = η J²
|
||
!
|
||
du(ien,1:im,1:jm,1:km) = du(ien,1:im,1:jm,1:km) &
|
||
+ resistivity * db(1:im,1:jm,1:km)
|
||
|
||
end if ! energy equation present
|
||
|
||
end if ! resistivity is not zero
|
||
|
||
end if ! ibx > 0
|
||
|
||
! add user defined source terms
|
||
!
|
||
call update_sources_user(pdata, t, dt, du(1:nv,1:im,1:jm,1:km))
|
||
|
||
#ifdef PROFILE
|
||
! stop accounting time for source terms
|
||
!
|
||
call stop_timer(imu)
|
||
#endif /* PROFILE */
|
||
|
||
!-------------------------------------------------------------------------------
|
||
!
|
||
end subroutine update_sources
|
||
|
||
!===============================================================================
|
||
!
|
||
end module sources
|