/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/gauge_hiju.f90

Go to the documentation of this file.

subroutine gauge_hiju
  use grid_parameter, only : nrg, ntg, npg
  use def_metric_hij, only : hxxd, hxyd, hxzd, hyyd, hyzd, hzzd, &
  &                          hxxu, hxyu, hxzu, hyyu, hyzu, hzzu
  use interface_poisson_solver
  use make_array_3d
  use interface_grgrad_4th_gridpoint
  use interface_grgrad_midpoint_type0
  implicit none
!
  real(8), pointer :: gx(:,:,:), gy(:,:,:), gz(:,:,:)
  real(8), pointer :: divg(:,:,:), 
                     sougx(:,:,:), sougy(:,:,:), sougz(:,:,:), &
                     bfn(:,:,:)
  real(8) :: dgxdx, dgxdy, dgxdz, dgydx, dgydy, dgydz, 
            dgzdx, dgzdy, dgzdz, &
            dhxxux, dhxxuy, dhxxuz, dhxyux, dhxyuy, dhxyuz, &
            dhxzux, dhxzuy, dhxzuz, &
            dhyxux, dhyxuy, dhyxuz, dhyyux, dhyyuy, dhyyuz, &
            dhyzux, dhyzuy, dhyzuz, &
            dhzxux, dhzxuy, dhzxuz, dhzyux, dhzyuy, dhzyuz, &
            dhzzux, dhzzuy, dhzzuz, &
            divgg, dlgxdx, dlgxdy, dlgxdz, dlgydy, dlgydz, dlgzdz, &
            fac23, hdhx, hdhy, hdhz, hxidiv, hxxuc, hxyuc, hxzuc, &
            hyidiv, hyxuc, hyyuc, hyzuc, hzidiv, hzxuc, hzyuc, hzzuc
  real(8) :: dbdx, dbdy, dbdz
  integer :: ipg, irg, itg
!
  call alloc_array3d(gx,0,nrg,0,ntg,0,npg)
  call alloc_array3d(gy,0,nrg,0,ntg,0,npg)
  call alloc_array3d(gz,0,nrg,0,ntg,0,npg)
  call alloc_array3d(divg,0,nrg,0,ntg,0,npg)
  call alloc_array3d(sougx,0,nrg,0,ntg,0,npg)
  call alloc_array3d(sougy,0,nrg,0,ntg,0,npg)
  call alloc_array3d(sougz,0,nrg,0,ntg,0,npg)
  call alloc_array3d(bfn,0,nrg,0,ntg,0,npg)

! --- Impose extended Dirac gauge condition.
!
  do ipg = 1, npg
    do itg = 1, ntg
      do irg = 1, nrg
!
        call grgrad_midpoint_type0(hxxu,dhxxux,dhxxuy,dhxxuz,irg,itg,ipg)
        call grgrad_midpoint_type0(hxyu,dhxyux,dhxyuy,dhxyuz,irg,itg,ipg)
        call grgrad_midpoint_type0(hxzu,dhxzux,dhxzuy,dhxzuz,irg,itg,ipg)
        call grgrad_midpoint_type0(hyyu,dhyyux,dhyyuy,dhyyuz,irg,itg,ipg)
        call grgrad_midpoint_type0(hyzu,dhyzux,dhyzuy,dhyzuz,irg,itg,ipg)
        call grgrad_midpoint_type0(hzzu,dhzzux,dhzzuy,dhzzuz,irg,itg,ipg)
        hxidiv = dhxxux + dhxyuy + dhxzuz
        hyidiv = dhxyux + dhyyuy + dhyzuz
        hzidiv = dhxzux + dhyzuy + dhzzuz
!
! -- original
        sougx(irg,itg,ipg) = hxidiv
        sougy(irg,itg,ipg) = hyidiv
        sougz(irg,itg,ipg) = hzidiv
!
      end do
    end do
  end do
!
  call poisson_solver(sougx,gx)
  call poisson_solver(sougy,gy)
  call poisson_solver(sougz,gz)
!
  do ipg = 1, npg
    do itg = 1, ntg
      do irg = 1, nrg
        call grgrad_midpoint_type0(gx,dgxdx,dgxdy,dgxdz,irg,itg,ipg)
        call grgrad_midpoint_type0(gy,dgydx,dgydy,dgydz,irg,itg,ipg)
        call grgrad_midpoint_type0(gz,dgzdx,dgzdy,dgzdz,irg,itg,ipg)
        divg(irg,itg,ipg) = dgxdx + dgydy + dgzdz
!shi      xxx = rg(irg)*sinthg(itg)*cosphig(ipg)
!shi      yyy = rg(irg)*sinthg(itg)*sinphig(ipg)
!shi      zzz = rg(irg)*costhg(itg)
!shi      divg(irg,itg,ipg) = xxx*gx(irg,itg,ipg)
!shi     &                  + yyy*gy(irg,itg,ipg)
!shi     &                  + zzz*gz(irg,itg,ipg)
      end do
    end do
  end do
!
  call poisson_solver(divg,bfn)
!
!shi      do 50 ipg = 0, npg
!shi      do 50 itg = 0, ntg
!shi      do 50 irg = 0, nrg
!shi      xxx = rg(irg)*sinthg(itg)*cosphig(ipg)
!shi      yyy = rg(irg)*sinthg(itg)*sinphig(ipg)
!shi      zzz = rg(irg)*costhg(itg)
!shi      bfn(irg,itg,ipg) = bfn(irg,itg,ipg) 
!shi     &                 - (xxx*gx(irg,itg,ipg)
!shi     &                  + yyy*gy(irg,itg,ipg)
!shi     &                  + zzz*gz(irg,itg,ipg))
!shi   50 continue
!
  do ipg = 0, npg
    do itg = 0, ntg
      do irg = 0, nrg
        call grgrad_4th_gridpoint(bfn,dbdx,dbdy,dbdz,irg,itg,ipg)
        gx(irg,itg,ipg) = gx(irg,itg,ipg) -0.25d0*dbdx
        gy(irg,itg,ipg) = gy(irg,itg,ipg) -0.25d0*dbdy
        gz(irg,itg,ipg) = gz(irg,itg,ipg) -0.25d0*dbdz
!shi      gx(irg,itg,ipg) = gx(irg,itg,ipg) + 0.125d0*dbdx
!shi      gy(irg,itg,ipg) = gy(irg,itg,ipg) + 0.125d0*dbdy
!shi      gz(irg,itg,ipg) = gz(irg,itg,ipg) + 0.125d0*dbdz
      end do
    end do
  end do
!
  fac23 = 2.0d0/3.0d0
  do ipg = 0, npg
    do itg = 0, ntg
      do irg = 0, nrg
!
        call grgrad_4th_gridpoint(gx,dgxdx,dgxdy,dgxdz,irg,itg,ipg)
        call grgrad_4th_gridpoint(gy,dgydx,dgydy,dgydz,irg,itg,ipg)
        call grgrad_4th_gridpoint(gz,dgzdx,dgzdy,dgzdz,irg,itg,ipg)
!
        divgg  = dgxdx + dgydy + dgzdz
        dlgxdx = dgxdx + dgxdx - fac23*divgg
        dlgxdy = dgxdy + dgydx
        dlgxdz = dgxdz + dgzdx
        dlgydy = dgydy + dgydy - fac23*divgg
        dlgydz = dgydz + dgzdy
        dlgzdz = dgzdz + dgzdz - fac23*divgg
!
        hxxu(irg,itg,ipg) = hxxu(irg,itg,ipg) - dlgxdx
        hxyu(irg,itg,ipg) = hxyu(irg,itg,ipg) - dlgxdy
        hxzu(irg,itg,ipg) = hxzu(irg,itg,ipg) - dlgxdz
        hyyu(irg,itg,ipg) = hyyu(irg,itg,ipg) - dlgydy
        hyzu(irg,itg,ipg) = hyzu(irg,itg,ipg) - dlgydz
        hzzu(irg,itg,ipg) = hzzu(irg,itg,ipg) - dlgzdz
!
      end do
    end do
  end do
!
  deallocate(gx)
  deallocate(gy)
  deallocate(gz)
  deallocate(divg)
  deallocate(sougx)
  deallocate(sougy)
  deallocate(sougz)
  deallocate(bfn)
end subroutine gauge_hiju

---