00001 subroutine hydrostatic_eq_peos(emd)
00002 use phys_constant, only : long
00003 use grid_parameter
00004 use def_matter, only : rs, omef, jomef, jomef_int, &
00005 & omeg, jomeg, jomeg_int
00006 use def_matter_parameter, only : ome, ber, ROT_LAW
00007 use def_metric, only : alph, psi, bvxd, bvyd, bvzd
00008 use coordinate_grav_r, only : rg
00009 use trigonometry_grav_theta, only : sinthg, costhg
00010 use trigonometry_grav_phi, only : sinphig, cosphig
00011 use def_vector_phi, only : vec_phif
00012 use def_vector_x, only : vec_xf, vec_xg
00013 use make_array_3d
00014 use interface_interpo_gr2fl
00015 implicit none
00016 real(long), pointer :: emd(:,:,:)
00017 real(long), pointer :: alphf(:,:,:), psif(:,:,:)
00018 real(long), pointer :: bvxdf(:,:,:), bvydf(:,:,:), bvzdf(:,:,:)
00019 real(long) :: vphif(3)
00020 real(long) :: ovufc(3), ovufc2
00021 real(long) :: alp_tmp, psi_tmp, bx_tmp, by_tmp
00022 real(long) :: omefc, jomefc, jomef_intfc, omegc, jomegc, jomeg_intgc
00023 real(long) :: hh, ut, pre, rho, ene, qq
00024 real(long) :: xx, yy, zz, Rcyl
00025 integer :: irf, itf, ipf
00026 call alloc_array3d(psif, 0, nrf, 0, ntf, 0, npf)
00027 call alloc_array3d(alphf, 0, nrf, 0, ntf, 0, npf)
00028 call alloc_array3d(bvxdf, 0, nrf, 0, ntf, 0, npf)
00029 call alloc_array3d(bvydf, 0, nrf, 0, ntf, 0, npf)
00030 call alloc_array3d(bvzdf, 0, nrf, 0, ntf, 0, npf)
00031 call interpo_gr2fl(alph, alphf)
00032 call interpo_gr2fl(psi, psif)
00033 call interpo_gr2fl(bvxd, bvxdf)
00034 call interpo_gr2fl(bvyd, bvydf)
00035 call interpo_gr2fl(bvzd, bvzdf)
00036
00037 if (ROT_LAW.eq.'DR') then
00038 do itf = 0, ntf
00039 do irf = 0, nrf
00040 ipf = 0
00041 xx = vec_xf(irf,itf,ipf,1)
00042 yy = vec_xf(irf,itf,ipf,2)
00043 zz = vec_xf(irf,itf,ipf,3)
00044 Rcyl = sqrt(xx**2 + yy**2)
00045 omefc = omef(irf,itf,ipf)
00046
00047 alp_tmp = alphf(irf,itf,ipf)
00048 psi_tmp = psif(irf,itf,ipf)
00049 bx_tmp = bvxdf(irf,itf,ipf)
00050 by_tmp = bvydf(irf,itf,ipf)
00051
00052 call calc_omega_drot(Rcyl,alp_tmp,psi_tmp,bx_tmp,by_tmp&
00053 & ,omefc,jomefc,jomef_intfc)
00054 omef(irf,itf,0:npf) = omefc
00055 jomef(irf,itf,0:npf) = jomefc
00056 jomef_int(irf,itf,0:npf) = jomef_intfc
00057
00058 xx = vec_xg(irf,itf,ipf,1)
00059 yy = vec_xg(irf,itf,ipf,2)
00060 zz = vec_xg(irf,itf,ipf,3)
00061 Rcyl = sqrt(xx**2 + yy**2)
00062 omegc = omeg(irf,itf,ipf)
00063 alp_tmp = alph(irf,itf,ipf)
00064 psi_tmp = psi(irf,itf,ipf)
00065 bx_tmp = bvxd(irf,itf,ipf)
00066 by_tmp = bvyd(irf,itf,ipf)
00067 call calc_omega_drot(Rcyl,alp_tmp,psi_tmp,bx_tmp,by_tmp&
00068 & ,omegc,jomegc,jomeg_intgc)
00069 omeg(irf,itf,0:npf) = omegc
00070 jomeg(irf,itf,0:npf) = jomegc
00071 jomeg_int(irf,itf,0:npf) = jomeg_intgc
00072 end do
00073 end do
00074 else
00075 omef(0:nrf,0:ntf,0:npf) = ome
00076 jomef(0:nrf,0:ntf,0:npf) = 0.0d0
00077 jomef_int(0:nrf,0:ntf,0:npf) = 0.0d0
00078 omeg(0:nrf,0:ntf,0:npf) = ome
00079 jomeg(0:nrf,0:ntf,0:npf) = 0.0d0
00080 jomeg_int(0:nrf,0:ntf,0:npf) = 0.0d0
00081 end if
00082
00083 do ipf = 0, npf
00084 do itf = 0, ntf
00085 do irf = 0, nrf
00086 vphif(1) = vec_phif(irf,itf,ipf,1)
00087 vphif(2) = vec_phif(irf,itf,ipf,2)
00088 vphif(3) = vec_phif(irf,itf,ipf,3)
00089 omefc = omef(irf,itf,ipf)
00090 ovufc(1) = bvxdf(irf,itf,ipf) + omefc*vphif(1)
00091 ovufc(2) = bvydf(irf,itf,ipf) + omefc*vphif(2)
00092 ovufc(3) = bvzdf(irf,itf,ipf) + omefc*vphif(3)
00093 ovufc2 = ovufc(1)**2 + ovufc(2)**2 + ovufc(3)**2
00094 ut = 1.0d0/sqrt(alphf(irf,itf,ipf)**2 &
00095 & - psif(irf,itf,ipf)**4*ovufc2)
00096 jomef_intfc = jomef_int(irf,itf,ipf)
00097 hh = ber*ut*dexp(-jomef_intfc)
00098 call peos_h2qprho(hh, qq, pre, rho, ene)
00099 emd(irf,itf,ipf) = qq
00100 end do
00101 end do
00102 end do
00103
00104 deallocate(alphf)
00105 deallocate(psif)
00106 deallocate(bvxdf)
00107 deallocate(bvydf)
00108 deallocate(bvzdf)
00109 end subroutine hydrostatic_eq_peos
