Cocal: /home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/update_parameter_BNS_peos

00001 subroutine update_parameter_BNS_peos_irrot(convf)
00002   use phys_constant, only  : long
00003   use grid_parameter, only : nrg, ntg, npg, nrf_deform, nrf, ntgxy, npgxzm
00004   use def_metric, only : psi, alph, bvxd, bvyd, bvzd, tfkij, tfkijkij
00005   use def_matter_parameter
00006 !  use def_velocity_rot
00007   use def_velocity_potential
00008   use def_matter
00009   use coordinate_grav_r, only : rg
00010   use trigonometry_grav_theta, only : sinthg, costhg
00011   use trigonometry_grav_phi, only   : sinphig, cosphig
00012   use def_vector_phi, only : vec_phig
00013   use interface_minv
00014   use interface_flgrad_4th_gridpoint
00015   implicit none
00016   real(long), intent(in) :: convf
00017   real(long) ::    sgg(1:4), gg(1:4,1:4)
00018   real(long) ::    ain, amid, aout, pin, pmid, pout
00019 ! ain is log(\alpha^((1/R_i)^2)) (in). R_in is in fact the previous R.
00020   real(long) :: a2ppin,  a2p4pin,  pm4in,  ppin,  p4pin,  bxin,  byin,  bzin
00021   real(long) :: a2ppmid, a2p4pmid, pm4mid, ppmid, p4pmid, bxmid, bymid, bzmid  
00022   real(long) :: a2ppout, a2p4pout, pm4out, ppout, p4pout, bxout, byout, bzout 
00023   real(long) :: dxvepin,  dyvepin,  dzvepin,  wxin,  wyin,  wzin
00024   real(long) :: dxvepmid, dyvepmid, dzvepmid, wxmid, wymid, wzmid
00025   real(long) :: dxvepout, dyvepout, dzvepout, wxout, wyout, wzout
00026   real(long) :: ovin, ovmid, ovout, lamin, lammid, lamout, lam2in, lam2mid, lam2out
00027   real(long) :: vphidvepin, vphidvepmid, vphidvepout
00028   real(long) :: dvep2in,  wdvepin,  w2in,  term1in,  term2in,  term3in,  term4in,  term5in,  term6in
00029   real(long) :: dvep2mid, wdvepmid, w2mid, term1mid, term2mid, term3mid, term4mid, term5mid, term6mid
00030   real(long) :: dvep2out, wdvepout, w2out, term1out, term2out, term3out, term4out, term5out, term6out
00031   real(long) :: vphiyin, vphiymid, vphiyout
00032   real(long) :: loghc
00033   real(long) :: radiold, berold, omeold
00034   real(long) :: ddradi, ddber, ddome
00035   real(long) :: facfac, numero, error, hc, pre, rho, ene
00036   integer    :: nin
00037   integer    :: irin, itin, ipin, irmid, itmid, ipmid, irout, itout, ipout
00038 !##############################################
00039   real(long) ::    ut, hh, emdtest
00040 !##############################################
00041   nin = nrf_deform
00042 !
00043   sgg = 0.0d0
00044   gg  = 0.0d0
00045 !
00046   irin = nrf;   itin = ntgxy;   ipin = 0;
00047   irmid= 0  ;   itmid= 0    ;   ipmid= 0;
00048   irout= nrf;   itout= ntgxy;   ipout= npgxzm;
00049 
00050   call peos_q2hprho(emdc, hc, pre, rho, ene)
00051   loghc = dlog(hc)
00052 
00053   vphiyin  = vec_phig(irin,  itin,  ipin,  2)
00054   vphiymid = vec_phig(irmid, itmid, ipmid, 2)
00055   vphiyout = vec_phig(irout, itout, ipout, 2)
00056 
00057   omeold = ome
00058   berold = ber
00059   radiold = radi
00060 
00061 !  write(6,'(1p,6e12.4)') ome, ber, radi, vphiyin, vphiymid, vphiyout
00062 
00063   numero = 1
00064 ! ################################################################
00065   ain     = dlog(alph(irin,itin,ipin)**(1.d0/radi**2))
00066   pin     = dlog( psi(irin,itin,ipin)**(1.d0/radi**2))
00067 !  a2ppin  = ( alph(irin,itin,ipin)**2 * psi(irin,itin,ipin)**(confpow)               )**(1.d0/radi**2)
00068 !  a2p4pin = ( alph(irin,itin,ipin)**2 * psi(irin,itin,ipin)**(4.0d0 + 2.0d0*confpow) )**(1.d0/radi**2)
00069   pm4in   = ( psi(irin,itin,ipin)**(-4)                                              )**(1.d0/radi**2)
00070 !  ppin    = ( psi(irin,itin,ipin)**(confpow)                                         )**(1.d0/radi**2)
00071 !  p4pin   = ( psi(irin,itin,ipin)**(4.0d0 + 2.0d0*confpow)                           )**(1.d0/radi**2)
00072   bxin    = bvxd(irin,itin,ipin)
00073   byin    = bvyd(irin,itin,ipin)
00074   bzin    = bvzd(irin,itin,ipin)
00075 !
00076 !  write(6,*)  ain, pin, confpow, a2ppin
00077 
00078   amid     = dlog(alph(irmid,itmid,ipmid)**(1.d0/radi**2))
00079   pmid     = dlog( psi(irmid,itmid,ipmid)**(1.d0/radi**2))
00080 !  a2ppmid  = ( alph(irmid,itmid,ipmid)**2 * psi(irmid,itmid,ipmid)**(confpow)               )**(1.d0/radi**2)
00081 !  a2p4pmid = ( alph(irmid,itmid,ipmid)**2 * psi(irmid,itmid,ipmid)**(4.0d0 + 2.0d0*confpow) )**(1.d0/radi**2)
00082   pm4mid   = ( psi(irmid,itmid,ipmid)**(-4)                                                 )**(1.d0/radi**2)
00083 !  ppmid    = ( psi(irmid,itmid,ipmid)**(confpow)                                            )**(1.d0/radi**2)
00084 !  p4pmid   = ( psi(irmid,itmid,ipmid)**(4.0d0 + 2.0d0*confpow)                              )**(1.d0/radi**2)
00085   bxmid    = bvxd(irmid,itmid,ipmid)
00086   bymid    = bvyd(irmid,itmid,ipmid)
00087   bzmid    = bvzd(irmid,itmid,ipmid)
00088 !
00089   aout     = dlog(alph(irout,itout,ipout)**(1.d0/radi**2))
00090   pout     = dlog( psi(irout,itout,ipout)**(1.d0/radi**2))
00091 !  a2ppout  = ( alph(irout,itout,ipout)**2 * psi(irout,itout,ipout)**(confpow)               )**(1.d0/radi**2)
00092 !  a2p4pout = ( alph(irout,itout,ipout)**2 * psi(irout,itout,ipout)**(4.0d0 + 2.0d0*confpow) )**(1.d0/radi**2)
00093   pm4out   = ( psi(irout,itout,ipout)**(-4)                                                 )**(1.d0/radi**2)
00094 !  ppout    = ( psi(irout,itout,ipout)**(confpow)                                            )**(1.d0/radi**2)
00095 !  p4pout   = ( psi(irout,itout,ipout)**(4.0d0 + 2.0d0*confpow)                              )**(1.d0/radi**2)
00096   bxout    = bvxd(irout,itout,ipout)
00097   byout    = bvyd(irout,itout,ipout)
00098   bzout    = bvzd(irout,itout,ipout)
00099 
00100 ! write (6,*) "test values"
00101 ! write (6,'(1p,3e12.4)') ome, ber, radi
00102 ! write (6,'(1p,3e12.4)') psi(irin,itin,ipin) ,psi(irmid,itmid,ipmid) ,psi(irout,itout,ipout)
00103 ! write (6,'(1p,3e12.4)') alph(irin,itin,ipin),alph(irmid,itmid,ipmid),alph(irout,itout,ipout)
00104 ! write (6,'(1p,3e12.4)') error
00105 ! write (6,'(1p,3e12.4)') pain, pamid, paout
00106 !  write(6,'(a14,f10.7,a19,f10.7,a12,f10.7,a7,f10.7,a6)') &
00107 !  &   "ff=-Log[ber]+(",ain,")*R*R+0.5*Log[1.0-(",pain,")^(2*R*R)*((",bin,")+ome*(",vphiyin,"))^2];"
00108 !  write(6,'(a3,f10.7,a11,f10.7,a19,f10.7,a12,f10.7,a7,f10.7,a6)') "gg=",loghc,    &
00109 !  &   "-Log[ber]+(",amid,")*R*R+0.5*Log[1.0-(",pamid,")^(2*R*R)*((",bmid,")+ome*(",vphiymid,"))^2];"
00110 !  write(6,'(a14,f10.7,a19,f10.7,a12,f10.7,a7,f10.7,a6)') &
00111 !  &   "hh=-Log[ber]+(",aout,")*R*R+0.5*Log[1.0-(",paout,")^(2*R*R)*((",bout,")+ome*(",vphiyout,"))^2];"
00112 !  write(6,'(a7,f10.7,a7,f10.7,a5,f10.7,a2)') ",{{ome,",ome,"},{ber,",ber,"},{R,",radi,"}}"
00113 ! ################################################################
00114   do
00115     dxvepin = vepxf(irin,itin,ipin)
00116     dyvepin = vepyf(irin,itin,ipin)
00117     dzvepin = vepzf(irin,itin,ipin)
00118    
00119     dxvepmid = vepxf(irmid,itmid,ipmid)
00120     dyvepmid = vepyf(irmid,itmid,ipmid)
00121     dzvepmid = vepzf(irmid,itmid,ipmid)
00122 
00123     dxvepout = vepxf(irout,itout,ipout)
00124     dyvepout = vepyf(irout,itout,ipout)
00125     dzvepout = vepzf(irout,itout,ipout)
00126 
00127     ovin  = byin  + ome*vphiyin
00128     ovmid = bymid + ome*vphiymid
00129     ovout = byout + ome*vphiyout
00130 
00131 !    lamin  = ber + (bxin +0)*dxvepin  +  ovin*dyvepin  + (bzin +0)*dzvepin
00132 !    lammid = ber + (bxmid+0)*dxvepmid + ovmid*dyvepmid + (bzmid+0)*dzvepmid
00133 !    lamout = ber + (bxout+0)*dxvepout + ovout*dyvepout + (bzout+0)*dzvepout
00134 
00135     lamin  = ber + ovin*dyvepin 
00136     lammid = ber + ovmid*dyvepmid 
00137     lamout = ber + ovout*dyvepout
00138 
00139     lam2in  = lamin*lamin
00140     lam2mid = lammid*lammid
00141     lam2out = lamout*lamout
00142  
00143     vphidvepin  =  vphiyin*dyvepin
00144     vphidvepmid = vphiymid*dyvepmid
00145     vphidvepout = vphiyout*dyvepout
00146 
00147 !    write(6,*)  vphiyin, dxvepin, dyvepin, dzvepin
00148 
00149     dvep2in = (dxvepin**2 + dyvepin**2 + dzvepin**2)*pm4in**(radi**2)
00150 
00151     term2in = dvep2in 
00152 
00153     term3in = 0.25d0 
00154 
00155     term4in = 1.0d0 + term2in/1.0d0
00156 
00157 
00158     dvep2mid = (dxvepmid**2 + dyvepmid**2 + dzvepmid**2)*pm4mid**(radi**2)
00159 
00160     term2mid = dvep2mid 
00161 
00162     term3mid = 0.25d0 
00163 
00164     term4mid = 1.0d0 + term2mid/hc/hc
00165 
00166 
00167     dvep2out = (dxvepout**2 + dyvepout**2 + dzvepout**2)*pm4out**(radi**2)
00168 
00169     term2out = dvep2out 
00170 
00171     term3out = 0.25d0 
00172 
00173     term4out = 1.0d0 + term2out/1.0d0
00174 !
00175     if(lamin<0 .or. lamout<0 .or. lammid<0) then
00176       write(6,*)  "Quantities lam inside log negative. Exiting..."
00177       stop
00178     end if
00179     if(term3in<0 .or. term3out<0 .or. term3mid<0) then
00180       write(6,*)  "Quantities inside sqrt negative. Exiting..."
00181       stop
00182     end if
00183     if(term4in<0 .or. term4out<0 .or. term4mid<0) then
00184       write(6,*)  "Quantities inside log negative. Exiting..."
00185       stop
00186     else
00187       sgg(1) = -( -dlog(lamin)  + radi**2*ain   + 0.5d0*dlog(term4in))
00188 
00189       sgg(2) = -( -dlog(lamout) + radi**2*aout  + 0.5d0*dlog(term4out))
00190 
00191       sgg(3) = -( -dlog(lammid) + radi**2*amid  + 0.5d0*dlog(term4mid) + loghc)
00192     end if
00193 !
00194     term6in  = dvep2in*2.0d0*radi*dlog(pm4in) 
00195 
00196     term6mid = dvep2mid*2.0d0*radi*dlog(pm4mid)
00197 
00198     term6out = dvep2out*2.0d0*radi*dlog(pm4out)
00199 
00200 !  ain is ln(\alpha)
00201     gg(1,1) = vphidvepin*(-1.0d0 )/lamin
00202 
00203     gg(1,2) = (-1.0d0 )/lamin
00204 
00205     gg(1,3) = 2.0d0*radi*ain + term6in/term4in/(2.0d0*1.0d0)
00206 
00207     gg(2,1) = vphidvepout*(-1.0d0 )/lamout
00208 
00209     gg(2,2) = (-1.0d0 )/lamout
00210 
00211     gg(2,3) = 2.0d0*radi*aout + term6out/term4out/(2.0d0*1.0d0)
00212 
00213     gg(3,1) = vphidvepmid*(-1.0d0 )/lammid
00214 
00215     gg(3,2) = (-1.0d0 )/lammid
00216 
00217     gg(3,3) = 2.0d0*radi*amid + term6mid/term4mid/(2.0d0*hc*hc)
00218 !
00219 !    test
00220 !if (numero.le.20) then
00221 ! write (6,'(1p,8e12.4)') sgg(1), sgg(2), sgg(3), pamid, radi, ov2mid, termmid
00222 ! write(6,*) "--------------------------------------------------------"
00223 ! write (6,'(1p,4f20.12)') gg(1,1), gg(1,2), gg(1,3), sgg(1)
00224 ! write (6,'(1p,4f20.12)') gg(2,1), gg(2,2), gg(2,3), sgg(2)
00225 ! write (6,'(1p,4f20.12)') gg(3,1), gg(3,2), gg(3,3), sgg(3)
00226 ! write(6,*) "--------------------------------------------------------"
00227 !end if
00228 !stop
00229 !!pause
00230 !    test
00231     call minv(gg,sgg,3,4)
00232 
00233 !    write (6,'(1p,4e12.4)') numero, sgg(1), sgg(2), sgg(3)
00234 !    write(6,*) "Old:", numero, ome, ber, radi
00235 
00236     ddome = sgg(1)
00237     ddber = sgg(2)
00238     ddradi = sgg(3)
00239     facfac = dmin1(dble(numero)/5.0d0,1.0d0)
00240     ome = ome + ddome*facfac
00241     ber = ber + ddber*facfac
00242     radi = radi + ddradi*facfac
00243 
00244 !    write(6,*) "New:", ome, ber, radi
00245     error = dmax1(dabs(ddome/ome),dabs(ddber/ber),dabs(ddradi/radi))
00246 !    numero = numero + 1
00247 !    if (numero<1000) then
00248 !      write(6,*)' numero = ', numero, '   error =',error
00249 !    end if
00250 !      if (iter<10.and.numero>10) exit
00251     if (numero>1010) exit
00252 !    if (dabs(error)<1.0d-08) exit
00253     if (dabs(error)<1.0d-12) then
00254 !      write(6,*) "New:", numero, ome, ber, radi
00255       exit
00256     end if
00257     numero = numero + 1
00258   end do
00259 !    test
00260 ! write (6,*) "test values"
00261 ! write (6,*) sgg(1), sgg(2), sgg(3)
00262 ! write (6,*) ome
00263 !    test
00264 !      write(6,*)' numero = ', numero, '  error = ',error
00265 !
00266   if (radi < 0.d0) write(6,*) ' ### radi minus ###'
00267   if (radi < 0.d0) radi = - radi
00268   if (ome < 0.d0) write(6,*) ' ### ome minus ###'
00269   if (ome < 0.d0) ome = - ome
00270   if (ber < 0.d0) write(6,*) ' ### ber minus ###'
00271   if (ber < 0.d0) ber = - ber
00272   ome = convf*ome + (1.d0-convf)*omeold
00273   ber = convf*ber + (1.d0-convf)*berold
00274   radi = convf*radi + (1.d0-convf)*radiold
00275 !
00276 !
00277 ! --  Improving alpha and psi.
00278 !
00279   alph = alph**((radi/radiold)**2)
00280   psi = psi**((radi/radiold)**2)
00281 !##############################################
00282 !
00283 !  stop
00284 end subroutine update_parameter_BNS_peos_irrot
00285
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/update_parameter_BNS_peos_irrot.f90
