!!****************************************************************************** !! !! 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-2019 Grzegorz Kowal !! !! 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 . !! !!****************************************************************************** !! !! module: COORDINATES !! !! This module provides variables and subroutines handling the coordinates !! for all refinement levels. !! !!****************************************************************************** ! module coordinates ! module variables are not implicit by default ! implicit none ! MODULE PARAMETERS: ! ================= ! ! the number of interior cells along each block dimension (and half of it) ! integer , save :: ncells = 8, ncells_half = 4 ! the number of ghost zones at each side of the block (and half and double of it) ! integer , save :: nghosts = 2, nghosts_half = 1, nghosts_double = 4 ! the domain block dimensions including the ghost zones ! integer , save :: im = 12, jm = 12, km = 1 ! the domain block dimensions including the ghost zones ! integer , save :: it = 11, jt = 11, kt = 1 ! the domain division - the number of blocks in each direction at the base level ! integer, dimension(3), save :: domain_base_dims = 1 ! the limits of refinement level ! integer , save :: minlev = 1, maxlev = 1, toplev = 1 ! block indices ! integer , save :: ih = 6, jh = 6, kh = 1 integer , save :: ib = 3, jb = 3, kb = 1 integer , save :: ie = 10, je = 10, ke = 1 integer , save :: ibl = 2, jbl = 2, kbl = 1 integer , save :: ibu = 4, jbu = 4, kbu = 1 integer , save :: iel = 9, jel = 9, kel = 1 integer , save :: ieu = 11, jeu = 11, keu = 1 ! the domain bounds ! real(kind=8), save :: xmin = 0.0d+00 real(kind=8), save :: xmax = 1.0d+00 real(kind=8), save :: xlen = 1.0d+00 real(kind=8), save :: ymin = 0.0d+00 real(kind=8), save :: ymax = 1.0d+00 real(kind=8), save :: ylen = 1.0d+00 real(kind=8), save :: zmin = 0.0d+00 real(kind=8), save :: zmax = 1.0d+00 real(kind=8), save :: zlen = 1.0d+00 ! the domain volume and its inversion ! real(kind=8), save :: vol = 1.0d+00 real(kind=8), save :: voli = 1.0d+00 ! the domain boundary areas ! real(kind=8), save :: xarea = 1.0d+00 real(kind=8), save :: yarea = 1.0d+00 real(kind=8), save :: zarea = 1.0d+00 ! flags to indicate periodicity of boundaries ! logical, dimension(3), save :: periodic = .true. ! the block coordinates for all levels of refinement ! real(kind=8), dimension(:,:), allocatable, save :: ax , ay , az real(kind=8), dimension(: ), allocatable, save :: adx , ady , adz, adr real(kind=8), dimension(: ), allocatable, save :: adxi, adyi, adzi real(kind=8), dimension(: ), allocatable, save :: advol ! define type for rectangular subarray description ! type rectangular integer, dimension(NDIMS) :: l ! indices of the lower corner integer, dimension(NDIMS) :: u ! indices of the upper corner end type rectangular ! the subarray indices to ghost and domain areas used for boundary exchange ! ('c' for copy, 'p' for prolongation, 'r' for restriction) ! #if NDIMS == 2 type(rectangular), dimension(2,2,NDIMS) , save :: edges_dc , edges_gc type(rectangular), dimension(2,2,NDIMS) , save :: edges_dp , edges_gp type(rectangular), dimension(2,2,NDIMS) , save :: edges_dr , edges_gr type(rectangular), dimension(2,2) , save :: corners_dc, corners_gc type(rectangular), dimension(2,2) , save :: corners_dp, corners_gp type(rectangular), dimension(2,2) , save :: corners_dr, corners_gr #endif /* NDIMS == 2 */ #if NDIMS == 3 type(rectangular), dimension(2,2,2,NDIMS), save :: faces_dc , faces_gc type(rectangular), dimension(2,2,2,NDIMS), save :: faces_dp , faces_gp type(rectangular), dimension(2,2,2,NDIMS), save :: faces_dr , faces_gr type(rectangular), dimension(2,2,2,NDIMS), save :: edges_dc , edges_gc type(rectangular), dimension(2,2,2,NDIMS), save :: edges_dp , edges_gp type(rectangular), dimension(2,2,2,NDIMS), save :: edges_dr , edges_gr type(rectangular), dimension(2,2,2) , save :: corners_dc, corners_gc type(rectangular), dimension(2,2,2) , save :: corners_dp, corners_gp type(rectangular), dimension(2,2,2) , save :: corners_dr, corners_gr #endif /* NDIMS == 3 */ ! by default everything is private ! public !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! !=============================================================================== ! ! subroutine INITIALIZE_COORDINATES: ! --------------------------------- ! ! Subroutine initializes mesh coordinates and other coordinate parameters. ! ! Arguments: ! ! ni - the number of interior cells per dimension in a block; ! ng - the number of ghosts cells at each side in a block; ! tlevel - the maximum refinement level for coordinate generation; ! bdims - the base level dimension in blocks; ! xmn, xmx, ymn, ymx, zmn, zmx - ! the domain extrema in each direction; ! verbose - flag determining if the subroutine should be verbose; ! iret - return flag of the procedure execution status; ! !=============================================================================== ! subroutine initialize_coordinates(ni, ng, tlevel, bdims, & xmn, xmx, ymn, ymx, zmn, zmx, verbose, iret) ! include external procedures ! use parameters, only : get_parameter ! local variables are not implicit by default ! implicit none ! subroutine arguments ! logical , intent(in) :: verbose integer , intent(in) :: ni, ng, tlevel integer, dimension(3), intent(in) :: bdims real(kind=8) , intent(in) :: xmn, xmx, ymn, ymx, zmn, zmx integer , intent(inout) :: iret ! local variables ! character(len=32) :: lbndry, ubndry integer :: i, j, k, l, p, q, r integer :: fi, fj, fk integer :: nh, nd, nn, nm, np, nr, nq, ns, nt, nu logical :: info ! !------------------------------------------------------------------------------- ! ! set the number of cells along each block dimension and ! the number of ghost zones ! ncells = ni ncells_half = ni / 2 nghosts = ng nghosts_half = ng / 2 nghosts_double = ng * 2 ! local variables ! nh = nghosts_half nd = nghosts_double ! calculate the block dimensions including ghost cells ! im = ni + 2 * ng jm = ni + 2 * ng #if NDIMS == 3 km = ni + 2 * ng #endif /* NDIMS == 3 */ ! prepare indices ! it = im - nh + 1 jt = jm - nh + 1 #if NDIMS == 3 kt = km - nh + 1 #endif /* NDIMS == 3 */ ! calculate block indices ! ih = im / 2 ib = ng + 1 ie = ng + ni ibl = ib - 1 ibu = ib + ng - 1 iel = ie - ng + 1 ieu = ie + 1 jh = jm / 2 jb = ng + 1 je = ng + ni jbl = jb - 1 jbu = jb + ng - 1 jel = je - ng + 1 jeu = je + 1 #if NDIMS == 3 kh = km / 2 kb = ng + 1 ke = ng + ni kbl = kb - 1 kbu = kb + ng - 1 kel = ke - ng + 1 keu = ke + 1 #endif /* NDIMS == 3 */ ! determine domain periodicity ! lbndry = "periodic" ubndry = "periodic" call get_parameter("xlbndry", lbndry) call get_parameter("xubndry", ubndry) periodic(1) = lbndry == "periodic" .and. ubndry == "periodic" lbndry = "periodic" ubndry = "periodic" call get_parameter("ylbndry", lbndry) call get_parameter("yubndry", ubndry) periodic(2) = lbndry == "periodic" .and. ubndry == "periodic" #if NDIMS == 3 lbndry = "periodic" ubndry = "periodic" call get_parameter("ylbndry", lbndry) call get_parameter("yubndry", ubndry) periodic(3) = lbndry == "periodic" .and. ubndry == "periodic" #endif /* NDIMS == 3 */ ! obtain the refinement level bounds ! call get_parameter("minlev", minlev) call get_parameter("maxlev", maxlev) ! set the top level ! toplev = max(tlevel, maxlev) ! get the domain base dimensions ! domain_base_dims(1:NDIMS) = bdims(1:NDIMS) ! set the domain bounds ! xmin = xmn xmax = xmx ymin = ymn ymax = ymx #if NDIMS == 3 zmin = zmn zmax = zmx #endif /* NDIMS == 3 */ ! calculate the domain sizes ! xlen = xmax - xmin ylen = ymax - ymin #if NDIMS == 3 zlen = zmax - zmin #endif /* NDIMS == 3 */ ! calculate the domain volume ! vol = xlen * ylen * zlen voli = 1.0d+00 / vol ! calculate the boundary areas ! xarea = ylen * zlen yarea = xlen * zlen zarea = xlen * ylen ! allocate space for coordinate variables ! allocate(ax (toplev, im)) allocate(ay (toplev, jm)) allocate(az (toplev, km)) allocate(adx (toplev)) allocate(ady (toplev)) allocate(adz (toplev)) allocate(adr (toplev)) allocate(adxi (toplev)) allocate(adyi (toplev)) allocate(adzi (toplev)) allocate(advol(toplev)) ! reset all coordinate variables to initial values ! ax (:,:) = 0.0d+00 ay (:,:) = 0.0d+00 az (:,:) = 0.0d+00 adx (:) = 1.0d+00 ady (:) = 1.0d+00 adz (:) = 1.0d+00 adr (:) = 1.0d+00 adxi (:) = 1.0d+00 adyi (:) = 1.0d+00 adzi (:) = 1.0d+00 advol(:) = 1.0d+00 ! generate the coordinate variables for each level ! do l = 1, toplev ! calculate the block resolution at each level ! nn = ni * 2**(l - 1) ! calculate the cell sizes for each level ! adx (l) = xlen / (domain_base_dims(1) * nn) ady (l) = ylen / (domain_base_dims(2) * nn) #if NDIMS == 3 adz (l) = zlen / (domain_base_dims(3) * nn) #endif /* NDIMS == 3 */ #if NDIMS == 2 adr (l) = sqrt(adx(l)**2 + ady(l)**2) #endif /* NDIMS == 2 */ #if NDIMS == 3 adr (l) = sqrt(adx(l)**2 + ady(l)**2 + adz(l)**2) #endif /* NDIMS == 3 */ ! calculate the inverse of cell size ! adxi(l) = 1.0d+00 / adx(l) adyi(l) = 1.0d+00 / ady(l) #if NDIMS == 3 adzi(l) = 1.0d+00 / adz(l) #endif /* NDIMS == 3 */ ! calculate the block coordinates for each level ! ax(l,:) = ((/(i, i = 1, im)/) - ng - 5.0d-01) * adx(l) ay(l,:) = ((/(j, j = 1, jm)/) - ng - 5.0d-01) * ady(l) #if NDIMS == 3 az(l,:) = ((/(k, k = 1, km)/) - ng - 5.0d-01) * adz(l) #endif /* NDIMS == 3 */ ! calculate the cell volume at each level ! advol(l) = adx(l) * ady(l) * adz(l) end do ! l = 1, toplev ! initialize ghost subarray indices ! np = ni + ng nm = ni - ng nr = ni - nd nq = ni - nh ns = ncells_half nt = np / 2 nu = nm / 2 #if NDIMS == 2 do j = 1, 2 fj = j - 1 q = 3 - j do i = 1, 2 fi = i - 1 p = 3 - i ! for edge copy ! edges_gc(i,j,1)%l(1) = ib + fi * ns edges_gc(i,j,1)%l(2) = 1 + fj * np edges_gc(i,j,2)%l(1) = 1 + fi * np edges_gc(i,j,2)%l(2) = jb + fj * ns edges_dc(i,q,1)%l(1) = ib + fi * ns edges_dc(i,q,1)%l(2) = jb + fj * nm edges_dc(p,j,2)%l(1) = ib + fi * nm edges_dc(p,j,2)%l(2) = jb + fj * ns edges_gc(i,j,1)%u(:) = edges_gc(i,j,1)%l(:) + (/ ns, ng /) - 1 edges_gc(i,j,2)%u(:) = edges_gc(i,j,2)%l(:) + (/ ng, ns /) - 1 edges_dc(i,q,1)%u(:) = edges_dc(i,q,1)%l(:) + (/ ns, ng /) - 1 edges_dc(p,j,2)%u(:) = edges_dc(p,j,2)%l(:) + (/ ng, ns /) - 1 ! for edge restriction ! edges_gr(i,j,1)%l(1) = ib + fi * ns edges_gr(i,j,1)%l(2) = 1 + fj * np edges_gr(i,j,2)%l(1) = 1 + fi * np edges_gr(i,j,2)%l(2) = jb + fj * ns edges_dr(i,q,1)%l(1) = ib edges_dr(i,q,1)%l(2) = jb + fj * nr edges_dr(p,j,2)%l(1) = ib + fi * nr edges_dr(p,j,2)%l(2) = jb edges_gr(i,j,1)%u(:) = edges_gr(i,j,1)%l(:) + (/ ns, ng /) - 1 edges_gr(i,j,2)%u(:) = edges_gr(i,j,2)%l(:) + (/ ng, ns /) - 1 edges_dr(i,q,1)%u(:) = edges_dr(i,q,1)%l(:) + (/ ni, nd /) - 1 edges_dr(p,j,2)%u(:) = edges_dr(p,j,2)%l(:) + (/ nd, ni /) - 1 ! for edge prolongation ! edges_gp(i,j,1)%l(1) = ib - fi * ng edges_gp(i,j,1)%l(2) = 1 + fj * np edges_gp(i,j,2)%l(1) = 1 + fi * np edges_gp(i,j,2)%l(2) = jb - fj * ng edges_dp(i,q,1)%l(1) = ib + fi * nu edges_dp(i,q,1)%l(2) = jb + fj * nq edges_dp(p,j,2)%l(1) = ib + fi * nq edges_dp(p,j,2)%l(2) = jb + fj * nu edges_gp(i,j,1)%u(:) = edges_gp(i,j,1)%l(:) + (/ np, ng /) - 1 edges_gp(i,j,2)%u(:) = edges_gp(i,j,2)%l(:) + (/ ng, np /) - 1 edges_dp(i,q,1)%u(:) = edges_dp(i,q,1)%l(:) + (/ nt, nh /) - 1 edges_dp(p,j,2)%u(:) = edges_dp(p,j,2)%l(:) + (/ nh, nt /) - 1 ! for corner copy ! corners_gc(i,j)%l(1) = 1 + fi * np corners_gc(i,j)%l(2) = 1 + fj * np corners_dc(p,q)%l(1) = ib + fi * nm corners_dc(p,q)%l(2) = jb + fj * nm corners_gc(i,j)%u(:) = corners_gc(i,j)%l(:) + ng - 1 corners_dc(p,q)%u(:) = corners_dc(p,q)%l(:) + ng - 1 ! for corner restriction ! corners_gr(i,j)%l(1) = 1 + fi * np corners_gr(i,j)%l(2) = 1 + fj * np corners_dr(p,q)%l(1) = ib + fi * nr corners_dr(p,q)%l(2) = jb + fj * nr corners_gr(i,j)%u(:) = corners_gr(i,j)%l(:) + ng - 1 corners_dr(p,q)%u(:) = corners_dr(p,q)%l(:) + nd - 1 ! for corner prolongation ! corners_gp(i,j)%l(1) = 1 + fi * np corners_gp(i,j)%l(2) = 1 + fj * np corners_dp(p,q)%l(1) = ib + fi * nq corners_dp(p,q)%l(2) = jb + fj * nq corners_gp(i,j)%u(:) = corners_gp(i,j)%l(:) + ng - 1 corners_dp(p,q)%u(:) = corners_dp(p,q)%l(:) + nh - 1 end do ! i = 1, 2 end do ! j = 1, 2 #endif /* NDIMS == 2 */ #if NDIMS == 3 do k = 1, 2 fk = k - 1 r = 3 - k do j = 1, 2 fj = j - 1 q = 3 - j do i = 1, 2 fi = i - 1 p = 3 - i ! for face copy ! faces_gc(i,j,k,1)%l(1) = 1 + fi * np faces_gc(i,j,k,1)%l(2) = jb + fj * ns faces_gc(i,j,k,1)%l(3) = kb + fk * ns faces_gc(i,j,k,2)%l(1) = ib + fi * ns faces_gc(i,j,k,2)%l(2) = 1 + fj * np faces_gc(i,j,k,2)%l(3) = kb + fk * ns faces_gc(i,j,k,3)%l(1) = ib + fi * ns faces_gc(i,j,k,3)%l(2) = jb + fj * ns faces_gc(i,j,k,3)%l(3) = 1 + fk * np faces_dc(p,j,k,1)%l(1) = ib + fi * nm faces_dc(p,j,k,1)%l(2) = jb + fj * ns faces_dc(p,j,k,1)%l(3) = kb + fk * ns faces_dc(i,q,k,2)%l(1) = ib + fi * ns faces_dc(i,q,k,2)%l(2) = jb + fj * nm faces_dc(i,q,k,2)%l(3) = kb + fk * ns faces_dc(i,j,r,3)%l(1) = ib + fi * ns faces_dc(i,j,r,3)%l(2) = jb + fj * ns faces_dc(i,j,r,3)%l(3) = kb + fk * nm faces_gc(i,j,k,1)%u(:) = faces_gc(i,j,k,1)%l(:) + (/ ng, ns, ns /) - 1 faces_gc(i,j,k,2)%u(:) = faces_gc(i,j,k,2)%l(:) + (/ ns, ng, ns /) - 1 faces_gc(i,j,k,3)%u(:) = faces_gc(i,j,k,3)%l(:) + (/ ns, ns, ng /) - 1 faces_dc(p,j,k,1)%u(:) = faces_dc(p,j,k,1)%l(:) + (/ ng, ns, ns /) - 1 faces_dc(i,q,k,2)%u(:) = faces_dc(i,q,k,2)%l(:) + (/ ns, ng, ns /) - 1 faces_dc(i,j,r,3)%u(:) = faces_dc(i,j,r,3)%l(:) + (/ ns, ns, ng /) - 1 ! for face restriction ! faces_gr(i,j,k,1)%l(1) = 1 + fi * np faces_gr(i,j,k,1)%l(2) = jb + fj * ns faces_gr(i,j,k,1)%l(3) = kb + fk * ns faces_gr(i,j,k,2)%l(1) = ib + fi * ns faces_gr(i,j,k,2)%l(2) = 1 + fj * np faces_gr(i,j,k,2)%l(3) = kb + fk * ns faces_gr(i,j,k,3)%l(1) = ib + fi * ns faces_gr(i,j,k,3)%l(2) = jb + fj * ns faces_gr(i,j,k,3)%l(3) = 1 + fk * np faces_dr(p,j,k,1)%l(1) = ib + fi * nr faces_dr(p,j,k,1)%l(2) = jb faces_dr(p,j,k,1)%l(3) = kb faces_dr(i,q,k,2)%l(1) = ib faces_dr(i,q,k,2)%l(2) = jb + fj * nr faces_dr(i,q,k,2)%l(3) = kb faces_dr(i,j,r,3)%l(1) = ib faces_dr(i,j,r,3)%l(2) = jb faces_dr(i,j,r,3)%l(3) = kb + fk * nr faces_gr(i,j,k,1)%u(:) = faces_gr(i,j,k,1)%l(:) + (/ ng, ns, ns /) - 1 faces_gr(i,j,k,2)%u(:) = faces_gr(i,j,k,2)%l(:) + (/ ns, ng, ns /) - 1 faces_gr(i,j,k,3)%u(:) = faces_gr(i,j,k,3)%l(:) + (/ ns, ns, ng /) - 1 faces_dr(p,j,k,1)%u(:) = faces_dr(p,j,k,1)%l(:) + (/ nd, ni, ni /) - 1 faces_dr(i,q,k,2)%u(:) = faces_dr(i,q,k,2)%l(:) + (/ ni, nd, ni /) - 1 faces_dr(i,j,r,3)%u(:) = faces_dr(i,j,r,3)%l(:) + (/ ni, ni, nd /) - 1 ! for face prolongation ! faces_gp(i,j,k,1)%l(1) = 1 + fi * np faces_gp(i,j,k,1)%l(2) = jb - fj * ng faces_gp(i,j,k,1)%l(3) = kb - fk * ng faces_gp(i,j,k,2)%l(1) = ib - fi * ng faces_gp(i,j,k,2)%l(2) = 1 + fj * np faces_gp(i,j,k,2)%l(3) = kb - fk * ng faces_gp(i,j,k,3)%l(1) = ib - fi * ng faces_gp(i,j,k,3)%l(2) = jb - fj * ng faces_gp(i,j,k,3)%l(3) = 1 + fk * np faces_dp(p,j,k,1)%l(1) = ib + fi * nq faces_dp(p,j,k,1)%l(2) = jb + fj * nu faces_dp(p,j,k,1)%l(3) = kb + fk * nu faces_dp(i,q,k,2)%l(1) = ib + fi * nu faces_dp(i,q,k,2)%l(2) = jb + fj * nq faces_dp(i,q,k,2)%l(3) = kb + fk * nu faces_dp(i,j,r,3)%l(1) = ib + fi * nu faces_dp(i,j,r,3)%l(2) = jb + fj * nu faces_dp(i,j,r,3)%l(3) = kb + fk * nq faces_gp(i,j,k,1)%u(:) = faces_gp(i,j,k,1)%l(:) + (/ ng, np, np /) - 1 faces_gp(i,j,k,2)%u(:) = faces_gp(i,j,k,2)%l(:) + (/ np, ng, np /) - 1 faces_gp(i,j,k,3)%u(:) = faces_gp(i,j,k,3)%l(:) + (/ np, np, ng /) - 1 faces_dp(p,j,k,1)%u(:) = faces_dp(p,j,k,1)%l(:) + (/ nh, nt, nt /) - 1 faces_dp(i,q,k,2)%u(:) = faces_dp(i,q,k,2)%l(:) + (/ nt, nh, nt /) - 1 faces_dp(i,j,r,3)%u(:) = faces_dp(i,j,r,3)%l(:) + (/ nt, nt, nh /) - 1 ! for edge copy ! edges_gc(i,j,k,1)%l(1) = ib + fi * ns edges_gc(i,j,k,1)%l(2) = 1 + fj * np edges_gc(i,j,k,1)%l(3) = 1 + fk * np edges_gc(i,j,k,2)%l(1) = 1 + fi * np edges_gc(i,j,k,2)%l(2) = jb + fj * ns edges_gc(i,j,k,2)%l(3) = 1 + fk * np edges_gc(i,j,k,3)%l(1) = 1 + fi * np edges_gc(i,j,k,3)%l(2) = 1 + fj * np edges_gc(i,j,k,3)%l(3) = kb + fk * ns edges_dc(i,q,r,1)%l(1) = ib + fi * ns edges_dc(i,q,r,1)%l(2) = jb + fj * nm edges_dc(i,q,r,1)%l(3) = kb + fk * nm edges_dc(p,j,r,2)%l(1) = jb + fi * nm edges_dc(p,j,r,2)%l(2) = jb + fj * ns edges_dc(p,j,r,2)%l(3) = kb + fk * nm edges_dc(p,q,k,3)%l(1) = ib + fi * nm edges_dc(p,q,k,3)%l(2) = jb + fj * nm edges_dc(p,q,k,3)%l(3) = kb + fk * ns edges_gc(i,j,k,1)%u(:) = edges_gc(i,j,k,1)%l(:) + (/ ns, ng, ng /) - 1 edges_gc(i,j,k,2)%u(:) = edges_gc(i,j,k,2)%l(:) + (/ ng, ns, ng /) - 1 edges_gc(i,j,k,3)%u(:) = edges_gc(i,j,k,3)%l(:) + (/ ng, ng, ns /) - 1 edges_dc(i,q,r,1)%u(:) = edges_dc(i,q,r,1)%l(:) + (/ ns, ng, ng /) - 1 edges_dc(p,j,r,2)%u(:) = edges_dc(p,j,r,2)%l(:) + (/ ng, ns, ng /) - 1 edges_dc(p,q,k,3)%u(:) = edges_dc(p,q,k,3)%l(:) + (/ ng, ng, ns /) - 1 ! for edge restriction ! edges_gr(i,j,k,1)%l(1) = ib + fi * ns edges_gr(i,j,k,1)%l(2) = 1 + fj * np edges_gr(i,j,k,1)%l(3) = 1 + fk * np edges_gr(i,j,k,2)%l(1) = 1 + fi * np edges_gr(i,j,k,2)%l(2) = jb + fj * ns edges_gr(i,j,k,2)%l(3) = 1 + fk * np edges_gr(i,j,k,3)%l(1) = 1 + fi * np edges_gr(i,j,k,3)%l(2) = 1 + fj * np edges_gr(i,j,k,3)%l(3) = kb + fk * ns edges_dr(i,q,r,1)%l(1) = ib edges_dr(i,q,r,1)%l(2) = jb + fj * nr edges_dr(i,q,r,1)%l(3) = kb + fk * nr edges_dr(p,j,r,2)%l(1) = ib + fi * nr edges_dr(p,j,r,2)%l(2) = jb edges_dr(p,j,r,2)%l(3) = kb + fk * nr edges_dr(p,q,k,3)%l(1) = ib + fi * nr edges_dr(p,q,k,3)%l(2) = jb + fj * nr edges_dr(p,q,k,3)%l(3) = kb edges_gr(i,j,k,1)%u(:) = edges_gr(i,j,k,1)%l(:) + (/ ns, ng, ng /) - 1 edges_gr(i,j,k,2)%u(:) = edges_gr(i,j,k,2)%l(:) + (/ ng, ns, ng /) - 1 edges_gr(i,j,k,3)%u(:) = edges_gr(i,j,k,3)%l(:) + (/ ng, ng, ns /) - 1 edges_dr(i,q,r,1)%u(:) = edges_dr(i,q,r,1)%l(:) + (/ ni, nd, nd /) - 1 edges_dr(p,j,r,2)%u(:) = edges_dr(p,j,r,2)%l(:) + (/ nd, ni, nd /) - 1 edges_dr(p,q,k,3)%u(:) = edges_dr(p,q,k,3)%l(:) + (/ nd, nd, ni /) - 1 ! for edge prolongation ! edges_gp(i,j,k,1)%l(1) = ib - fi * ng edges_gp(i,j,k,1)%l(2) = 1 + fj * np edges_gp(i,j,k,1)%l(3) = 1 + fk * np edges_gp(i,j,k,2)%l(1) = 1 + fi * np edges_gp(i,j,k,2)%l(2) = jb - fj * ng edges_gp(i,j,k,2)%l(3) = 1 + fk * np edges_gp(i,j,k,3)%l(1) = 1 + fi * np edges_gp(i,j,k,3)%l(2) = 1 + fj * np edges_gp(i,j,k,3)%l(3) = kb - fk * ng edges_dp(i,q,r,1)%l(1) = ib + fi * nu edges_dp(i,q,r,1)%l(2) = jb + fj * nq edges_dp(i,q,r,1)%l(3) = kb + fk * nq edges_dp(p,j,r,2)%l(1) = ib + fi * nq edges_dp(p,j,r,2)%l(2) = jb + fj * nu edges_dp(p,j,r,2)%l(3) = kb + fk * nq edges_dp(p,q,k,3)%l(1) = ib + fi * nq edges_dp(p,q,k,3)%l(2) = jb + fj * nq edges_dp(p,q,k,3)%l(3) = kb + fk * nu edges_gp(i,j,k,1)%u(:) = edges_gp(i,j,k,1)%l(:) + (/ np, ng, ng /) - 1 edges_gp(i,j,k,2)%u(:) = edges_gp(i,j,k,2)%l(:) + (/ ng, np, ng /) - 1 edges_gp(i,j,k,3)%u(:) = edges_gp(i,j,k,3)%l(:) + (/ ng, ng, np /) - 1 edges_dp(i,q,r,1)%u(:) = edges_dp(i,q,r,1)%l(:) + (/ nt, nh, nh /) - 1 edges_dp(p,j,r,2)%u(:) = edges_dp(p,j,r,2)%l(:) + (/ nh, nt, nh /) - 1 edges_dp(p,q,k,3)%u(:) = edges_dp(p,q,k,3)%l(:) + (/ nh, nh, nt /) - 1 ! for corner copy ! corners_gc(i,j,k)%l(1) = 1 + fi * np corners_gc(i,j,k)%l(2) = 1 + fj * np corners_gc(i,j,k)%l(3) = 1 + fk * np corners_dc(p,q,r)%l(1) = ib + fi * nm corners_dc(p,q,r)%l(2) = jb + fj * nm corners_dc(p,q,r)%l(3) = kb + fk * nm corners_gc(i,j,k)%u(:) = corners_gc(i,j,k)%l(:) + ng - 1 corners_dc(p,q,r)%u(:) = corners_dc(p,q,r)%l(:) + ng - 1 ! for corner restriction ! corners_gr(i,j,k)%l(1) = 1 + fi * np corners_gr(i,j,k)%l(2) = 1 + fj * np corners_gr(i,j,k)%l(3) = 1 + fk * np corners_dr(p,q,r)%l(1) = ib + fi * nr corners_dr(p,q,r)%l(2) = jb + fj * nr corners_dr(p,q,r)%l(3) = kb + fk * nr corners_gr(i,j,k)%u(:) = corners_gr(i,j,k)%l(:) + ng - 1 corners_dr(p,q,r)%u(:) = corners_dr(p,q,r)%l(:) + nd - 1 ! for corner prolongation ! corners_gp(i,j,k)%l(1) = 1 + fi * np corners_gp(i,j,k)%l(2) = 1 + fj * np corners_gp(i,j,k)%l(3) = 1 + fk * np corners_dp(p,q,r)%l(1) = ib + fi * nq corners_dp(p,q,r)%l(2) = jb + fj * nq corners_dp(p,q,r)%l(3) = kb + fk * nq corners_gp(i,j,k)%u(:) = corners_gp(i,j,k)%l(:) + ng - 1 corners_dp(p,q,r)%u(:) = corners_dp(p,q,r)%l(:) + nh - 1 end do ! i = 1, 2 end do ! j = 1, 2 end do ! k = 1, 2 #endif /* NDIMS == 3 */ !------------------------------------------------------------------------------- ! end subroutine initialize_coordinates ! !=============================================================================== ! ! subroutine FINALIZE_COORDINATES: ! ------------------------------- ! ! Subroutine deallocates mesh coordinates. ! ! Arguments: ! ! iret - return flag of the procedure execution status; ! !=============================================================================== ! subroutine finalize_coordinates(iret) ! local variables are not implicit by default ! implicit none ! subroutine arguments ! integer, intent(inout) :: iret ! !------------------------------------------------------------------------------- ! ! deallocating coordinate variables ! if (allocated(ax) ) deallocate(ax) if (allocated(ay) ) deallocate(ay) if (allocated(az) ) deallocate(az) if (allocated(adx) ) deallocate(adx) if (allocated(ady) ) deallocate(ady) if (allocated(adz) ) deallocate(adz) if (allocated(adr) ) deallocate(adr) if (allocated(adxi) ) deallocate(adxi) if (allocated(adyi) ) deallocate(adyi) if (allocated(adzi) ) deallocate(adzi) if (allocated(advol)) deallocate(advol) !------------------------------------------------------------------------------- ! end subroutine finalize_coordinates ! !=============================================================================== ! ! subroutine PRINT_COORDINATES: ! ---------------------------- ! ! Subroutine print module parameters. ! ! Arguments: ! ! verbose - flag determining if the subroutine should be verbose; ! !=============================================================================== ! subroutine print_coordinates(verbose) ! include external procedures ! use helpers , only : print_section, print_parameter ! local variables are not implicit by default ! implicit none ! subroutine arguments ! logical, intent(in) :: verbose ! local variables ! character(len= 80) :: sfmt character(len=255) :: msg integer :: p, q ! local arrays ! integer(kind=4), dimension(3) :: bm, rm, cm, fm ! !------------------------------------------------------------------------------- ! if (verbose) then ! the base and top level block dimensions, and the base and effective resolution ! bm(:) = domain_base_dims(:) rm(:) = bm(:) * 2**(maxlev - 1) cm(:) = bm(:) * ncells fm(:) = rm(:) * ncells call print_section(verbose, "Geometry") call print_parameter(verbose, "refinement to level" , maxlev ) call print_parameter(verbose, "number of block cells", ncells ) call print_parameter(verbose, "number of ghost zones", nghosts) write(msg,"(i0)") maxval(fm(1:NDIMS)) p = len(trim(adjustl(msg))) write(msg,"(i0)") maxval(rm(1:NDIMS)) q = len(trim(adjustl(msg))) #if NDIMS == 3 write(sfmt,"(6(a,i0),a)") "('[',i", p, ",' x ',i", p, ",' x ',i", p, & ",'] ([',i", q, ",' x ',i", q, ",' x ',i", q, ",'] blocks)')" #else /* NDIMS == 3 */ write(sfmt,"(4(a,i0),a)") "('[',i", p, ",' x ',i", p, & ",'] ([',i", q, ",' x ',i", q, ",'] blocks)')" #endif /* NDIMS == 3 */ write(msg,sfmt) cm(1:NDIMS), bm(1:NDIMS) call print_parameter(verbose, "base resolution" , msg ) write(msg,sfmt) fm(1:NDIMS), rm(1:NDIMS) call print_parameter(verbose, "effective resolution" , msg ) sfmt = "(1es12.5,1x,'...',1x,1es12.5)" write(msg,sfmt) xmin, xmax call print_parameter(verbose, "X-bounds" , msg ) write(msg,sfmt) ymin, ymax call print_parameter(verbose, "Y-bounds" , msg ) #if NDIMS == 3 write(msg,sfmt) zmin, zmax call print_parameter(verbose, "Z-bounds" , msg ) #endif /* NDIMS */ end if !------------------------------------------------------------------------------- ! end subroutine print_coordinates !=============================================================================== ! end module