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

Go to the documentation of this file.

subroutine poisson_solver_homogeneous_sol(surp,vpot_b)
  use grid_parameter, only : nrf, ntf, npf
  use def_metric_on_SFC, only : hxxdf, hxydf, hxzdf, hyydf, hyzdf, hzzdf
  use trigonometry_grav_theta, only : hcosecg
  use trigonometry_grav_phi, only : sinmpg, hsinmpg, hcosmpg
  use def_matter, only : rs
  use def_velocity_potential, only : alm
  use def_vector_phi, only : hvec_phif_surface
  use interface_interpo_linear_surface_type0
  use interface_flgrad_midpoint_surface_type0
  use interface_flgrad_midpoint_type0_2Dsurf_sph
  use interface_interpo_linear_type0_2Dsurf
  use make_array1d
  use make_array2d
  implicit none
  real(long), pointer :: surp(:,:), vpot_b(:,:,:)
  real(long), pointer :: sgg(:), gg(:,:)
  real(long) :: hhxxu, hhxyu, hhxzu, hhyxu, hhyyu, hhyzu, 
               hhzxu, hhzyu, hhzzu
  real(long) :: gamxxu, gamxyu, gamxzu, gamyxu, gamyyu, gamyzu, 
               gamzxu, gamzyu, gamzzu
  real(long) :: dxvep, dyvep, dzvep, dtrs, dprs, hrs
  real(long) :: rnr, rnth, rnphi, rnx, rny, rnz
  real(long) :: erx, ery, erz, ethx, ethy, ethz, ephix, ephiy, ephiz
  integer    :: ir, it, ip
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     Coefficients of harmonic solution is computed.
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
  nbou = (nlg/2+1)*nlg/2
  call allocate_array1d(sgg,1,nbou)
  call allocate_array2d(gg,1,nbou,1,nbou)
!
  sgg(1:nbou) = 0.0d0
   gg(1:nbou,1:nbou) = 0.0d0
  alm(1:nbou,1:nbou) = 0.0d0
!
! --- Boundary equations are solved.
!
  ir = nrf
  ieqs = 0
  do il = 1, nlg
    do im = 2-mod(il,2), il, 2
      ieqs = ieqs + 1
!
      do ip = 1, npf
        do it = 1, ntf
!
          call interpo_linear_surface_type0(emdfc,emd,ir,it,ip)
          call interpo_linear_surface_type0(utfc,utf,ir,it,ip)
          call peos_q2hprho(emdfc, hhfc, prefc, rhofc, enefc)
          hhut = hhfc*utfc
          call interpo_linear_surface_type0(psifc,psif,ir,it,ip)
          psif4 = psifc**4
          call flgrad_midpoint_type0_2Dsurf_sph(rs,dtrs,dprs,it,ip)
          call interpo_linear_type0_2Dsurf(hrs,rs,it,ip)
          call calc_surface_normal_midpoint(rs,rnx,rny,rnz,it,ip)
          call interpo_linear_surface_type0(bvxufc,bvxuf,ir,it,ip)
          call interpo_linear_surface_type0(bvyufc,bvyuf,ir,it,ip)
          call interpo_linear_surface_type0(bvzufc,bvzuf,ir,it,ip)
          rotshx = bvxufc + ome*hvec_phif_surface(it,ip,1)
          rotshy = bvyufc + ome*hvec_phif_surface(it,ip,2)
          rotshz = bvzufc + ome*hvec_phif_surface(it,ip,3)
!
          flm = dble(il)*hrs**(il-1)*hylmg(il,im,it)*hsinmpg(im,ip) - &
          &     hrs**(il-2)*dtrs*hdtyplmg(il,im,it)*hsinmpg(im,ip) -  &
          &     hrs**(il-2)*hcosecg(it)**2*dprs*hyplmg(il,im,it)*    &
          &     dble(im)*hcosmpg(im,ip)
!
          sgg(ieqs) = sgg(ieqs) + ( - surp(it,ip) + &
          &   hhut*psif4*(rotshx*rnx + rotshy*rny + rotshz*rnz) )*flm
!
          ivec = 0
          do ill = 1, nlg
            do imm = 2-mod(ill,2), ill, 2
              ivec = ivec + 1
!
              fllmm =
              & dble(ill)*hrs**(ill-1)*hyplmg(ill,imm,it)*hsinmpg(imm,ip) -
              & hrs**(ill-2)*dtrs*hdtyplmg(ill,imm,it)*hsinmpg(imm,ip) - 
              & hrs**(ill-2)*hcosecg(it)**2*dprs*hyplmg(ill,imm,it)*
              & dble(imm)*hcosmpg(imm,ip)
!
              gg(ieqs,ivec) = gg(ieqs,ivec) + fllmm*flm
!
            end do
          end do
!
          if(ivec.ne.nbou) write(6,*) ' harmonic ivec wrong ', ivec, nbou
!
        end do
      end do
!
    end do
  end do
!
  if(ieqs.ne.nbou) write(6,*) ' harmonic ieqs wrong ', ieqs, nbou
!
  call minv(gg,sgg,nbou,nbou)
!
  ieqs = 0
  do il = 1, nlg
    do im = 2-mod(il,2), il, 2
      ieqs = ieqs + 1
      alm(il,im) = sgg(ieqs)
    end do
  end do
  if(ieqs.ne.nbou) write(6,*) ' harmonic alm wrong ', ieqs, nbou
!
! --  calculate homogeneous part.
!
  do ip = 0, npf
    do it = 0, ntf
      do ir = 0, nrf
!
        surf = 0.0d0
        do il = 1, nlg
          do im = 2-mod(il,2), il, 2
            surf = surf + alm(il,im)*(r(ir)*rs(it,ip))**il* &
            &             yplmg(il,im,it)*sinmp(im,ip) 
         end do
        end do
!
        vpot_b(ir,it,ip) = surf
!
      end do
    end do
  end do
!
!  if (iter.eq.itmx) then
!    do 141 il = 1, nlg
!      do 141 im = 2-mod(il,2), il, 2
!      write(24,*)il,im,alm(il,im)
! 141  continue
!      end if
!
  deallocate(sgg)
  deallocate(gg)
!
end subroutine poisson_solver_homogeneous_sol

---