diff --git a/sources/sources.F90 b/sources/sources.F90 index 6d48a69..041e9dd 100644 --- a/sources/sources.F90 +++ b/sources/sources.F90 @@ -268,7 +268,7 @@ module sources ! include external variables ! use blocks , only : block_data - use coordinates , only : im, jm, km + use coordinates , only : nb => bcells 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 @@ -283,13 +283,13 @@ module sources ! 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 + type(block_data), pointer , intent(inout) :: pdata + real(kind=8) , intent(in) :: t, dt + real(kind=8), dimension(:,:,:,:), intent(inout) :: du ! local variables ! - integer :: i , j , k + integer :: i, j, k = 1 real(kind=8) :: fc, gc real(kind=8) :: gx, gy, gz real(kind=8) :: dbx, dby, dbz @@ -300,11 +300,17 @@ module sources ! 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 + real(kind=8), dimension(nb) :: x + real(kind=8), dimension(nb) :: y +#if NDIMS == 3 + real(kind=8), dimension(nb) :: z + real(kind=8), dimension(nb,nb,nb) :: db + real(kind=8), dimension(3,3,nb,nb,nb) :: tmp +#else /* NDIMS == 3 */ + real(kind=8), dimension( 1) :: z + real(kind=8), dimension(nb,nb, 1) :: db + real(kind=8), dimension(3,3,nb,nb, 1) :: tmp +#endif /* NDIMS == 3 */ ! !------------------------------------------------------------------------------- ! @@ -320,17 +326,19 @@ module sources ! 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) + x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) + y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) #if NDIMS == 3 - z(1:km) = pdata%meta%zmin + az(pdata%meta%level,1:km) + z(:) = pdata%meta%zmin + az(pdata%meta%level,:) #endif /* NDIMS == 3 */ ! iterate over all positions in the YZ plane ! - do k = 1, km - do j = 1, jm - do i = 1, im +#if NDIMS == 3 + do k = 1, nb +#endif /* NDIMS == 3 */ + do j = 1, nb + do i = 1, nb ! get gravitational acceleration components ! @@ -368,9 +376,11 @@ module sources end if - end do ! i = 1, im - end do ! j = 1, jm - end do ! k = 1, km + end do ! i = 1, nb + end do ! j = 1, nb +#if NDIMS == 3 + end do ! k = 1, nb +#endif /* NDIMS == 3 */ end if ! gravity enabled @@ -386,18 +396,17 @@ module sources ! 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)) + call gradient(dh(:), pdata%q(ivx,:,:,:), tmp(inx,inx:inz,:,:,:)) + call gradient(dh(:), pdata%q(ivy,:,:,:), tmp(iny,inx:inz,:,:,:)) + call gradient(dh(:), pdata%q(ivz,:,:,:), tmp(inz,inx:inz,:,:,:)) ! iterate over all cells ! - do k = 1, km - do j = 1, jm - do i = 1, im +#if NDIMS == 3 + do k = 1, nb +#endif /* NDIMS == 3 */ + do j = 1, nb + do i = 1, nb ! prepare the νρ factor ! @@ -429,36 +438,35 @@ module sources 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 + end do ! i = 1, nb + end do ! j = 1, nb +#if NDIMS == 3 + end do ! k = 1, nb +#endif /* NDIMS == 3 */ ! 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)) + call divergence(dh(:), tmp(inx,inx:inz,:,:,:), db(:,:,:)) ! 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) + du(imx,:,:,:) = du(imx,:,:,:) + db(:,:,:) ! 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)) + call divergence(dh(:), tmp(iny,inx:inz,:,:,:), db(:,:,:)) ! 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) + du(imy,:,:,:) = du(imy,:,:,:) + db(:,:,:) ! 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)) + call divergence(dh(:), tmp(inz,inx:inz,:,:,:), db(:,:,:)) ! 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) + du(imz,:,:,:) = du(imz,:,:,:) + db(:,:,:) ! add viscous source term to total energy equation ! @@ -466,9 +474,11 @@ module sources ! iterate over all cells ! - do k = 1, km - do j = 1, jm - do i = 1, im +#if NDIMS == 3 + do k = 1, nb +#endif /* NDIMS == 3 */ + do j = 1, nb + do i = 1, nb ! calculate scalar product of v and viscous stress tensor τ ! @@ -488,18 +498,19 @@ module sources 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 + end do ! i = 1, nb + end do ! j = 1, nb +#if NDIMS == 3 + end do ! k = 1, nb +#endif /* NDIMS == 3 */ ! 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)) + call divergence(dh(:), tmp(inx,inx:inz,:,:,:), db(:,:,:)) ! 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) + du(ien,:,:,:) = du(ien,:,:,:) + db(:,:,:) end if ! ien > 0 @@ -583,26 +594,19 @@ module sources ! 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)) + call laplace(dh(:), pdata%q(ibx,:,:,:), tmp(inx,inx,:,:,:)) + call laplace(dh(:), pdata%q(iby,:,:,:), tmp(inx,iny,:,:,:)) + call laplace(dh(:), pdata%q(ibz,:,:,:), tmp(inx,inz,:,:,:)) ! 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) + tmp(iny,inx:inz,:,:,:) = resistivity * tmp(inx,inx:inz,:,:,:) ! 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) + du(ibx,:,:,:) = du(ibx,:,:,:) + tmp(iny,inx,:,:,:) + du(iby,:,:,:) = du(iby,:,:,:) + tmp(iny,iny,:,:,:) + du(ibz,:,:,:) = du(ibz,:,:,:) + tmp(iny,inz,:,:,:) ! update energy equation ! @@ -611,27 +615,24 @@ module sources ! 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)) + du(ien,:,:,:) = du(ien,:,:,:) & + + (pdata%q(ibx,:,:,:) * tmp(iny,inx,:,:,:) & + + pdata%q(iby,:,:,:) * tmp(iny,iny,:,:,:) & + + pdata%q(ibz,:,:,:) * tmp(iny,inz,:,:,:)) ! 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)) + call curl(dh(:), pdata%q(ibx:ibz,:,:,:), tmp(inz,inx:inz,:,:,:)) ! 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 + db(:,:,:) = tmp(inz,inx,:,:,:)**2 + tmp(inz,iny,:,:,:)**2 & + + tmp(inz,inz,:,:,:)**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) + du(ien,:,:,:) = du(ien,:,:,:) + resistivity * db(:,:,:) end if ! energy equation present @@ -641,7 +642,7 @@ module sources ! add user defined source terms ! - call update_sources_user(pdata, t, dt, du(1:nv,1:im,1:jm,1:km)) + call update_sources_user(pdata, t, dt, du(:,:,:,:)) #ifdef PROFILE ! stop accounting time for source terms