00001 subroutine source_proper_mass_peos(souf)
00002 use phys_constant, only : long, pi
00003 use grid_parameter, only : nrg, ntg, npg, nrf, ntf, npf
00004 use def_matter, only : emd, utf
00005 use def_matter_parameter, only : radi, ber, pinx
00006 use def_metric, only : tfkijkij, psi, alph
00007 use make_array_3d
00008 use interface_interpo_gr2fl
00009 implicit none
00010 real(long),pointer :: souf(:,:,:)
00011 real(long), pointer :: alphf(:,:,:), psif(:,:,:)
00012 real(long) :: psiwm6
00013 real(long) :: emdw, alphw, psiw, epsilonw, hhw, utw, rhoHw, esseS
00014 real(long) :: zfac, small = 1.0d-15
00015 real(long) :: otermx, otermy, otermz, prew, rhow
00016 integer :: ir,it,ip
00017
00018 call alloc_array3d(psif, 0, nrf, 0, ntf, 0, npf)
00019 call alloc_array3d(alphf, 0, nrf, 0, ntf, 0, npf)
00020 call interpo_gr2fl(alph, alphf)
00021 call interpo_gr2fl(psi, psif)
00022
00023
00024 do ip = 0, npf
00025 do it = 0, ntf
00026 do ir = 0, nrf
00027 emdw = emd(ir,it,ip)
00028 if (emdw <= small) emdw = small
00029 psiw = psif(ir,it,ip)
00030 alphw = alphf(ir,it,ip)
00031 call peos_q2hprho(emdw, hhw, prew, rhow, epsilonw)
00032 utw = utf(ir,it,ip)
00033
00034 souf(ir,it,ip) = epsilonw*alphw*utw*psiw**6
00035 end do
00036 end do
00037 end do
00038
00039 deallocate(alphf)
00040 deallocate(psif)
00041 end subroutine source_proper_mass_peos