mod_lfff.t

Go to the documentation of this file.

!> Program to extrapolate linear force-free fields in 3D Cartesian coordinates,
!> based on exact Green function method (Chiu & Hilton 1977 ApJ 212,873).
!>
!> Usage:
!> 1 In the subroutine usr_set_parameters of mod_usr.t:
!>  To extrapolate a linear force free field from a observed magnetogram 
!>  prepared in a data file, e.g., 'hmiM720sxxxx.dat' replace 
!>  call init_bc_fff_data('hmiM720sxxxx.dat',unit_length,unit_magneticfield)
!>  'hmiM720sxxxx.dat' must be a binary file containing nx1,nx2,xc1,xc2,dxm1,
!>  dxm2, Bz0(nx1,nx2). Integers nx1 and nx2 give the resolution of the 
!>  uniform-grid magentogram. Others are double-precision floats. xc1 and xc2
!>  are coordinates of the central point of the magnetogram. dxm1 and dxm2 
!>  are the cell sizes for each direction, Bz0 is the vertical conponent 
!>  of magetic field on the solar surface from observations.
!>2 In the subroutine usr_init_one_grid of mod_usr.t,
!>  add lines like:
!>
!>  double precision :: Bf(ixG^S,1:ndir), alpha, zshift
!>
!>  alpha=0.d0     ! potential field
!>  !alpha=0.08d0  ! non-potential linear force-free field
!>  zshift=0.05d0  ! lift your box zshift heigher to the bottom magnetogram
!>  call calc_lin_fff(ixG^L,ix^L,Bf,x,alpha,zshift) 
!>
!>3 Notice that the resolution of input magnetogram must be better than the best
!>  resolution of your AMR grid to have a good behavior close to the bottom layer
module mod_lfff
  implicit none
  
  double precision, save :: bzmax,darea
  double precision, allocatable, save :: bz0(:,:)
  double precision, allocatable, save :: xa1(:),xa2(:)
  integer, save :: nx1,nx2
  
contains
 
  subroutine init_b_fff_data(magnetogramname,qLunit,qBunit)
    use mod_global_parameters
    use mod_comm_lib, only: mpistop
 
    double precision, intent(in) :: qLunit,qBunit
    double precision :: xc1,xc2,dxm1,dxm2
    integer, dimension(MPI_STATUS_SIZE) :: statuss
    integer :: file_handle,i
    logical :: aexist
    character(len=*), intent(in) :: magnetogramname
    ! nx1,nx2 are numbers of cells for each direction
    ! xc1,xc2 are coordinates of the central point of the magnetogram
    ! dxm1,dxm2 are cell sizes for each direction
    ! Bz0 is the 2D Bz magnetogram
    inquire(file=magnetogramname,exist=aexist)
    if(.not. aexist) then
      if(mype==0) write(*,'(2a)') "can not find file:",magnetogramname
      call mpistop("no input magnetogram----init_b_fff_data")
    end if
    call mpi_file_open(icomm,magnetogramname,mpi_mode_rdonly,mpi_info_null,&
                       file_handle,ierrmpi)
    call mpi_file_read_all(file_handle,nx1,1,mpi_integer,statuss,ierrmpi)
    call mpi_file_read_all(file_handle,nx2,1,mpi_integer,statuss,ierrmpi)
    allocate(bz0(nx1,nx2))
    call mpi_file_read_all(file_handle,xc1,1,mpi_double_precision,statuss,ierrmpi)
    call mpi_file_read_all(file_handle,xc2,1,mpi_double_precision,statuss,ierrmpi)
    call mpi_file_read_all(file_handle,dxm1,1,mpi_double_precision,statuss,ierrmpi)
    call mpi_file_read_all(file_handle,dxm2,1,mpi_double_precision,statuss,ierrmpi)
    call mpi_file_read_all(file_handle,bz0,nx1*nx2,mpi_double_precision,&
                           statuss,ierrmpi)
    call mpi_file_close(file_handle,ierrmpi)
    allocate(xa1(nx1))
    allocate(xa2(nx2))
    do i=1,nx1 
      xa1(i) = xc1 + (dble(i) - dble(nx1)/2.d0 - 0.5d0)*dxm1
    enddo
    do i=1,nx2
      xa2(i) = xc2 + (dble(i) - dble(nx2)/2.d0 - 0.5d0)*dxm2
    enddo
    ! declare and define global variables Lunit and Bunit to be your length unit in
    ! cm and magnetic strength unit in Gauss first
    dxm1=dxm1/qlunit
    dxm2=dxm2/qlunit
    xa1=xa1/qlunit
    xa2=xa2/qlunit
    darea=dxm1*dxm2
    bz0=bz0/qbunit
    bzmax=maxval(dabs(bz0(:,:)))
    
    ! normalize b
    bz0=bz0/bzmax
    if(mype==0) then
      print*,'magnetogram xrange:',minval(xa1),maxval(xa1)
      print*,'magnetogram yrange:',minval(xa2),maxval(xa2)
    end if
    
    if(mype==0) then
      print*,'extrapolating 3D force-free field from an observed Bz '
      print*,'magnetogram of',nx1,'by',nx2,'pixels. Bzmax=',bzmax
    endif
  
  end subroutine init_b_fff_data
 
{^ifthreed
  subroutine calc_lin_fff(ixI^L,ixO^L,Bf,x,alpha,zshift,idir)
  ! PURPOSE: 
  ! Calculation to determine linear FFF from the field on 
  ! the lower boundary (Chiu and Hilton 1977 ApJ 212,873). 
  ! NOTE: Only works for Cartesian coordinates 
  ! INPUT: Bf,x
  ! OUTPUT: updated b in w 
    use mod_global_parameters
 
    integer, intent(in) :: ixI^L, ixO^L
    integer, optional, intent(in) :: idir
    double precision, intent(in) :: x(ixI^S,1:ndim),alpha,zshift
    double precision, intent(inout) :: Bf(ixI^S,1:ndir)
 
    double precision, dimension(ixO^S) :: cos_az,sin_az,zk,bigr,r,r2,r3,cos_ar,sin_ar,g,dgdz
    double precision :: gx(ixO^S,1:ndim),twopiinv
    integer :: idim,ixp1,ixp2
 
    bf=0.d0
    twopiinv = 0.5d0/dpi*bzmax*darea
    ! get cos and sin arrays
    zk(ixo^s)=x(ixo^s,3)-xprobmin3+zshift
    cos_az(ixo^s)=dcos(alpha*zk(ixo^s))
    sin_az(ixo^s)=dsin(alpha*zk(ixo^s))
    ! looping Bz0 pixels
    do ixp2=1,nx2
      do ixp1=1,nx1
        bigr(ixo^s)=dsqrt((x(ixo^s,1)-xa1(ixp1))**2+&
                          (x(ixo^s,2)-xa2(ixp2))**2)
        r2=bigr**2+zk**2
        r=dsqrt(r2)
        r3=r**3
        cos_ar=dcos(alpha*r)
        sin_ar=dsin(alpha*r)
        where(bigr/=0.d0)
          bigr=1.d0/bigr
        end where
        where(r/=0.d0)
          r=1.d0/r
        end where
        where(r2/=0.d0)
          r2=1.d0/r2
        end where
        where(r3/=0.d0)
          r3=1.d0/r3
        end where
        g=(zk*cos_ar*r-cos_az)*bigr
        dgdz=(cos_ar*(r-zk**2*r3)-alpha*zk**2*sin_ar*r2+alpha*sin_az)*bigr
        do idim=1,ndim
          if(present(idir)) then
            if(idim/=idir) cycle
          end if
          select case(idim)
          case(1)
            bf(ixo^s,1)=bf(ixo^s,1)+bz0(ixp1,ixp2)*((x(ixo^s,1)-xa1(ixp1))*dgdz(ixo^s)&
                     +alpha*g(ixo^s)*(x(ixo^s,2)-xa2(ixp2)))*bigr(ixo^s)
          case(2)
            bf(ixo^s,2)=bf(ixo^s,2)+bz0(ixp1,ixp2)*((x(ixo^s,2)-xa2(ixp2))*dgdz(ixo^s)&
                     -alpha*g(ixo^s)*(x(ixo^s,1)-xa1(ixp1)))*bigr(ixo^s)
          case(3)
            bf(ixo^s,3)=bf(ixo^s,3)+bz0(ixp1,ixp2)*(zk(ixo^s)*cos_ar(ixo^s)*r3(ixo^s)+alpha*&
                                        zk(ixo^s)*sin_ar(ixo^s)*r2(ixo^s))
          end select
        end do
      end do
    end do
    bf(ixo^s,:)=bf(ixo^s,:)*twopiinv
 
  end subroutine calc_lin_fff
 
  subroutine get_potential_field_potential(ixI^L,ixO^L,potential,x,zshift)
  ! PURPOSE: 
  ! Calculation scalar potential of potential field given
  ! Bz at photosphere (Schmidt 1964 NASSP). 
  ! NOTE: Only works for Cartesian coordinates 
  ! INPUT: x,zshift
  ! OUTPUT: potential
    use mod_global_parameters
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: x(ixI^S,1:ndim),zshift
    double precision, intent(inout) :: potential(ixI^S)
 
    double precision :: zk
    integer :: ixp1,ixp2,ix^D
 
    potential=0.d0
    ! looping Bz0 pixels see equation (2)
    !$OMP PARALLEL DO
    do ix3=ixomin3,ixomax3
      zk=x(ixomin1,ixomin2,ix3,3)-xprobmin3+zshift
      do ix2=ixomin2,ixomax2
        do ix1=ixomin1,ixomax1
          do ixp2=1,nx2
            do ixp1=1,nx1
              potential(ix^d)=potential(ix^d)+0.5d0*bz0(ixp1,ixp2)*darea/&
                 (dpi*dsqrt((x(ix^d,1)-xa1(ixp1))**2+(x(ix^d,2)-xa2(ixp2))**2+zk**2))
            end do
          end do
        end do
      end do
    end do
    !$OMP END PARALLEL DO
  end subroutine get_potential_field_potential
 
  subroutine get_potential_field_potential_sphere(ixI^L,x,potential,nth,nph,magnetogram,theta,phi,r_sphere)
  ! PURPOSE: 
  ! Calculation scalar potential of potential field given
  ! Bz at photosphere (Schmidt 1964 NASSP). 
  ! NOTE: Only works for spherical coordinates 
  ! OUTPUT: potential
    use mod_global_parameters
    integer, intent(in) :: ixI^L,nth,nph
    real*8, intent(in) :: x(ixi^s,1:ndim)
    ! magnetogram Br on photosphere
    real*8, intent(in) :: magnetogram(nth,nph)
    ! theta and phi grid of the photospheric magnetogram
    real*8, intent(in) :: theta(nth),phi(nph)
    ! radius of photosphere
    real*8, intent(in) :: r_sphere
    real*8, intent(out) :: potential(ixi^s)
 
    real*8 :: area(nth),dtheta_half,dphi,inv2pi
    integer :: ixp1,ixp2,ix^D
 
    potential=0.d0
    ! assume uniformly discretized theta and phi
    dtheta_half=0.5d0*(theta(2)-theta(1))
    dphi=phi(2)-phi(1)
    area(1:nth)=2.d0*r_sphere**2*sin(theta(1:nth))*sin(dtheta_half)*sin(dphi)
    inv2pi=1.d0/(2.d0*dpi)
 
    !$OMP PARALLEL DO
    do ix3=iximin3,iximax3
      do ix2=iximin2,iximax2
        do ix1=iximin1,iximax1
          do ixp2=1,nph
            do ixp1=1,nth
              potential(ix^d)=potential(ix^d)+inv2pi*magnetogram(ixp1,ixp2)*area(ixp1)/&
                dsqrt(x(ix^d,1)**2+r_sphere**2-2.d0*x(ix^d,1)*r_sphere*&
                (dsin(x(ix^d,2))*dsin(theta(ixp1))*dcos(phi(ixp2)-x(ix^d,3))+dcos(x(ix^d,2))*&
                 dcos(theta(ixp1))))
            end do
          end do
        end do
      end do
    end do
    !$OMP END PARALLEL DO
 
  end subroutine get_potential_field_potential_sphere
 
  !> get potential magnetic field energy given normal B on all boundaries
  subroutine potential_field_energy_mg(benergy)
  ! NOTE: Solve Poisson equation of scalar potential using multigrid solver
  ! Only works for 3D Cartesian coordinates 
    use mod_global_parameters
    use mod_multigrid_coupling
    use mod_forest
    use mod_geometry
    real*8, intent(out) :: benergy
    real*8  :: block_benergy(ixm^t), mype_benergy
    real*8, allocatable :: tmp(:,:,:)
    integer :: iigrid, igrid, idir
    integer :: id,nc, lvl, ixI^L
 
    call get_potential_field_potential_mg()
    iximin^d=ixmlo^d-1;
    iximax^d=ixmhi^d+1;
    allocate(tmp(ixi^s))
    mype_benergy=0.d0
    do iigrid=1,igridstail; igrid=igrids(iigrid);
      block=>ps(igrid)
      ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
      id = igrid_to_node(igrid, mype)%node%id
      lvl= mg%boxes(id)%lvl
      nc = mg%box_size_lvl(lvl)
      block_benergy=0.d0
      do idir=1,ndim
        call gradient(mg%boxes(id)%cc({0:nc+1},mg_iphi),ixi^l,ixm^ll,idir,tmp)
        block_benergy(ixm^t)=block_benergy(ixm^t)+tmp(ixm^t)**2
      end do
      mype_benergy=mype_benergy+sum(0.5d0*block_benergy(ixm^t)*block%dvolume(ixm^t))
    end do
 
    call mpi_allreduce(mype_benergy,benergy,1,mpi_double_precision,mpi_sum,icomm,ierrmpi)
 
  end subroutine potential_field_energy_mg
 
  !>  Solve Poisson equation of scalar potential using multigrid solver
  subroutine get_potential_field_potential_mg()
    use mod_global_parameters
    use mod_comm_lib, only: mpistop
    use mod_multigrid_coupling
    double precision :: max_residual, res
    integer :: i, max_its
 
    mg%operator_type = mg_laplacian
    max_its=50
    max_residual=1.d-8
 
    ! Set boundary conditions
    mg%bc(:, mg_iphi)%bc_type = mg_bc_neumann
 
    do i=1,2*ndim
       mg%bc(i, mg_iphi)%boundary_cond => multigrid_bc
    end do
 
    ! Solve laplacian(phi) = 0
    do i=1,max_its
       !call mg_fas_vcycle(mg, max_res=res)
       call mg_fas_fmg(mg,.true.,max_res=res)
       if(mype==0) write(*,*) 'mg iteration',i,'residual:', res
       if(res < max_residual) exit
    end do
    if(res > max_residual) call mpistop("get potential multigrid: no convergence")
 
  end subroutine get_potential_field_potential_mg
 
  !> To set boundary condition on physical boundaries for mg Poisson solver
  subroutine multigrid_bc(box, nc, iv, nb, bc_type, bc)
    use mod_global_parameters
    use mod_multigrid_coupling
    type(mg_box_t), intent(in)    :: box
    integer, intent(in)           :: nc
    integer, intent(in)           :: iv      !< Index of variable
    integer, intent(in)           :: nb      !< number of boundary from 1 to 6 for 3D
    integer, intent(out)          :: bc_type !< Type of b.c.
    ! mg boundary values
    double precision, intent(out) :: bc(nc, nc)
    ! mg coordinates of boundary layer
    double precision              :: rr(nc, nc, 3)
    double precision :: rmina,rminb,rmaxa,rmaxb,xmina,xminb,xmaxa,xmaxb
    ! store normal magnetic field on the boundary
    double precision :: wbn(ixG^T)
    double precision, allocatable :: xcoarse(:,:,:)
    integer :: iigrid,igrid,ix^D,idir,ixbca,ixbcb,ixbcn,dlvl,wnc
 
    bc_type = mg_bc_neumann
 
    call mg_get_face_coords(box, nb, nc, rr)
 
    idir=(nb+1)/2
    bc=0.d0
    ! send normal B on boundaries from AMRVAC to mg Poisson solver
    ! for Neumann boundary condition
    select case(idir)
    case(1)
      if(mod(nb,2)==0) then
        ixbcn=ixmhi1
      else
        ixbcn=ixmlo1-1
      end if
      rmina=rr(1,1,2)-0.5d0*box%dr(2)
      rmaxa=rr(nc,1,2)+0.5d0*box%dr(2)
      rminb=rr(1,1,3)-0.5d0*box%dr(3)
      rmaxb=rr(1,nc,3)+0.5d0*box%dr(3)
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        block=>ps(igrid)
        if(.not.block%is_physical_boundary(nb)) cycle
        if(stagger_grid) then
          !wbn(ixbcn^%1ixG^T)=block%ws(ixbcn^%1ixG^T,idir)
          wbn(ixbcn^%1ixg^t)=block%ws(ixbcn^%1ixg^t,idir)
        else
          wbn(ixbcn^%1ixg^t)=half*(block%w(ixbcn^%1ixg^t,iw_mag(idir))+block%w(ixbcn+1^%1ixg^t,iw_mag(idir)))
        end if
        xmina=block%x(1,1,1,2)-0.5d0*rnode(rpdx2_,igrid)
        xmaxa=block%x(1,ixghi2,1,2)+0.5d0*rnode(rpdx2_,igrid)
        xminb=block%x(1,1,1,3)-0.5d0*rnode(rpdx3_,igrid)
        xmaxb=block%x(1,1,ixghi3,3)+0.5d0*rnode(rpdx3_,igrid)
        if(xmina<rr(1,1,2) .and. xmaxa>rr(nc,1,2) .and. &
           xminb<rr(1,1,3) .and. xmaxb>rr(1,nc,3)) then
          do ix2=1,nc
            do ix1=1,nc
              ixbca=ceiling((rr(ix1,ix2,2)-xmina)/rnode(rpdx2_,igrid))
              ixbcb=ceiling((rr(ix1,ix2,3)-xminb)/rnode(rpdx3_,igrid))
              bc(ix1,ix2)=wbn(ixbcn,ixbca,ixbcb)
            end do
          end do
        else if(block%x(1,ixmlo2,1,2)>rmina .and. block%x(1,ixmhi2,1,2)<rmaxa .and. &
                block%x(1,1,ixmlo3,3)>rminb .and. block%x(1,1,ixmhi3,3)<rmaxb) then
          dlvl=node(plevel_,igrid)-box%lvl
          wnc=nc/2**dlvl
          allocate(xcoarse(wnc,wnc,2))
          do ix2=1,wnc
            do ix1=1,wnc
              xcoarse(ix1,ix2,1)=sum(block%x(1,(ix1-1)*2**dlvl+1+nghostcells:ix1*2**dlvl+nghostcells,1,2))/dble(2**dlvl)
              xcoarse(ix1,ix2,2)=sum(block%x(1,1,(ix2-1)*2**dlvl+1+nghostcells:ix2*2**dlvl+nghostcells,3))/dble(2**dlvl)
              ixbca=ceiling((xcoarse(ix1,ix2,1)-rmina)/box%dr(2))
              ixbcb=ceiling((xcoarse(ix1,ix2,2)-rminb)/box%dr(3))
              bc(ixbca,ixbcb)=sum(wbn(ixbcn,(ix1-1)*2**dlvl+1+nghostcells:ix1*2**dlvl+nghostcells,&
                                  (ix2-1)*2**dlvl+1+nghostcells:ix2*2**dlvl+nghostcells))/dble(2**(2*dlvl))
            end do
          end do
          deallocate(xcoarse)
        end if
      end do
    case(2)
      if(mod(nb,2)==0) then
        ixbcn=ixmhi2
      else
        ixbcn=ixmlo2-1
      end if
      rmina=rr(1,1,1)-0.5d0*box%dr(1)
      rmaxa=rr(nc,1,1)+0.5d0*box%dr(1)
      rminb=rr(1,1,3)-0.5d0*box%dr(3)
      rmaxb=rr(1,nc,3)+0.5d0*box%dr(3)
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        block=>ps(igrid)
        if(.not.block%is_physical_boundary(nb)) cycle
        if(stagger_grid) then
          wbn(ixbcn^%2ixg^t)=block%ws(ixbcn^%2ixg^t,idir)
        else
          wbn(ixbcn^%2ixg^t)=half*(block%w(ixbcn^%2ixg^t,iw_mag(idir))+block%w(ixbcn+1^%2ixg^t,iw_mag(idir)))
        end if
        xmina=block%x(1,1,1,1)-0.5d0*rnode(rpdx1_,igrid)
        xmaxa=block%x(ixghi1,1,1,1)+0.5d0*rnode(rpdx1_,igrid)
        xminb=block%x(1,1,1,3)-0.5d0*rnode(rpdx3_,igrid)
        xmaxb=block%x(1,1,ixghi3,3)+0.5d0*rnode(rpdx3_,igrid)
        if(xmina<rr(1,1,1) .and. xmaxa>rr(nc,1,1) .and. &
           xminb<rr(1,1,3) .and. xmaxb>rr(1,nc,3)) then
          do ix2=1,nc
            do ix1=1,nc
              ixbca=ceiling((rr(ix1,ix2,1)-xmina)/rnode(rpdx1_,igrid))
              ixbcb=ceiling((rr(ix1,ix2,3)-xminb)/rnode(rpdx3_,igrid))
              bc(ix1,ix2)=wbn(ixbca,ixbcn,ixbcb)
            end do
          end do
        else if(block%x(ixmlo1,1,1,1)>rmina .and. block%x(ixmhi1,1,1,1)<rmaxa .and. &
                block%x(1,1,ixmlo3,3)>rminb .and. block%x(1,1,ixmhi3,3)<rmaxb) then
          dlvl=node(plevel_,igrid)-box%lvl
          wnc=nc/2**dlvl
          allocate(xcoarse(wnc,wnc,2))
          do ix2=1,wnc
            do ix1=1,wnc
              xcoarse(ix1,ix2,1)=sum(block%x((ix1-1)*2**dlvl+1+nghostcells:ix1*2**dlvl+nghostcells,1,1,1))/dble(2**dlvl)
              xcoarse(ix1,ix2,2)=sum(block%x(1,1,(ix2-1)*2**dlvl+1+nghostcells:ix2*2**dlvl+nghostcells,3))/dble(2**dlvl)
              ixbca=ceiling((xcoarse(ix1,ix2,1)-rmina)/box%dr(1))
              ixbcb=ceiling((xcoarse(ix1,ix2,2)-rminb)/box%dr(3))
              bc(ixbca,ixbcb)=sum(wbn((ix1-1)*2**dlvl+1+nghostcells:ix1*2**dlvl+nghostcells,ixbcn,&
                                  (ix2-1)*2**dlvl+1+nghostcells:ix2*2**dlvl+nghostcells))/dble(2**(2*dlvl))
            end do
          end do
          deallocate(xcoarse)
        end if
      end do
    case(3)
      if(mod(nb,2)==0) then
        ixbcn=ixmhi3
      else
        ixbcn=ixmlo3-1
      end if
      rmina=rr(1,1,1)-0.5d0*box%dr(1)
      rmaxa=rr(nc,1,1)+0.5d0*box%dr(1)
      rminb=rr(1,1,2)-0.5d0*box%dr(2)
      rmaxb=rr(1,nc,2)+0.5d0*box%dr(2)
      do iigrid=1,igridstail; igrid=igrids(iigrid);
        block=>ps(igrid)
        if(.not.block%is_physical_boundary(nb)) cycle
        if(stagger_grid) then
          wbn(ixbcn^%3ixg^t)=block%ws(ixbcn^%3ixg^t,idir)
        else
          wbn(ixbcn^%3ixg^t)=half*(block%w(ixbcn^%3ixg^t,iw_mag(idir))+block%w(ixbcn+1^%3ixg^t,iw_mag(idir)))
        end if
        xmina=block%x(1,1,1,1)-0.5d0*rnode(rpdx1_,igrid)
        xmaxa=block%x(ixghi1,1,1,1)+0.5d0*rnode(rpdx1_,igrid)
        xminb=block%x(1,1,1,2)-0.5d0*rnode(rpdx2_,igrid)
        xmaxb=block%x(1,ixghi2,1,2)+0.5d0*rnode(rpdx2_,igrid)
        if(xmina<rr(1,1,1) .and. xmaxa>rr(nc,1,1) .and. &
           xminb<rr(1,1,2) .and. xmaxb>rr(1,nc,2)) then
          do ix2=1,nc
            do ix1=1,nc
              ixbca=ceiling((rr(ix1,ix2,1)-xmina)/rnode(rpdx1_,igrid))
              ixbcb=ceiling((rr(ix1,ix2,2)-xminb)/rnode(rpdx2_,igrid))
              bc(ix1,ix2)=wbn(ixbca,ixbcb,ixbcn)
            end do
          end do
        else if(block%x(ixmlo1,1,1,1)>rmina .and. block%x(ixmhi1,1,1,1)<rmaxa .and. &
                block%x(1,ixmlo2,1,2)>rminb .and. block%x(1,ixmhi2,1,2)<rmaxb) then
          dlvl=node(plevel_,igrid)-box%lvl
          wnc=nc/2**dlvl
          allocate(xcoarse(wnc,wnc,2))
          do ix2=1,wnc
            do ix1=1,wnc
              xcoarse(ix1,ix2,1)=sum(block%x((ix1-1)*2**dlvl+1+nghostcells:ix1*2**dlvl+nghostcells,1,1,1))/dble(2**dlvl)
              xcoarse(ix1,ix2,2)=sum(block%x(1,(ix2-1)*2**dlvl+1+nghostcells:ix2*2**dlvl+nghostcells,1,2))/dble(2**dlvl)
              ixbca=ceiling((xcoarse(ix1,ix2,1)-rmina)/box%dr(1))
              ixbcb=ceiling((xcoarse(ix1,ix2,2)-rminb)/box%dr(2))
              bc(ixbca,ixbcb)=sum(wbn((ix1-1)*2**dlvl+1+nghostcells:ix1*2**dlvl+nghostcells,&
                                  (ix2-1)*2**dlvl+1+nghostcells:ix2*2**dlvl+nghostcells,ixbcn))/dble(2**(2*dlvl))
            end do
          end do
          deallocate(xcoarse)
        end if
      end do
    end select
 
  end subroutine multigrid_bc
}
end module mod_lfff