00001 subroutine source_rest_mass(souf)
00002 use phys_constant, only : long, pi
00003 use grid_parameter, only : nrg, ntg, npg, nrf, ntf, npf
00004 use def_matter, only : emdg, emd
00005 use def_matter_parameter, only : radi, pinx, ber
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, rhow, hhw, utw, rhoHw, esseS
00014 real(long) :: zfac, small = 1.0d-15
00015 real(long) :: otermx, otermy, otermz
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 rhow = emdw**pinx
00032 hhw = 1.0d0 + (pinx+1.0d0)*emdw
00033 utw = hhw/ber
00034
00035 souf(ir,it,ip) = rhow*alphw*utw*psiw**6
00036 end do
00037 end do
00038 end do
00039
00040 deallocate(alphf)
00041 deallocate(psif)
00042 end subroutine source_rest_mass
