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

00001 subroutine update_parameter_spherical_qeos(convf)
00002   use phys_constant, only  : long
00003   use grid_parameter, only : nrg, ntg, npg, nrf_deform, nrf, ntgxy
00004   use def_metric, only : psi, alph, bvxd, bvyd, bvzd, tfkij, tfkijkij
00005   use def_matter_parameter
00006   use coordinate_grav_r, only : rg
00007   use trigonometry_grav_theta, only : sinthg, costhg
00008   use trigonometry_grav_phi, only   : sinphig, cosphig
00009   use def_vector_phi, only : vec_phig
00010   use interface_minv
00011   implicit none
00012   real(long), intent(in) :: convf
00013   real(long)     ::    sgg(1:4), gg(1:4,1:4)
00014   real(long)     ::    bin, bmid, bout
00015 ! "bin" is beta(in) and so on.
00016   real(long)     ::    ain, amid, aout
00017 ! ain is log(\alpha^((1/R_i)^2)) (in). R_in is in fact the previous R.
00018   real(long)     ::    pain, pamid, paout
00019 ! pain = (\phi^2/\alpha)^((1/R_i)^2) (in).
00020   real(long)     ::    ovin, ovmid, ovout
00021 ! ovin is \omega(in). \omega = \beta + \Omega \phi.
00022   real(long)     ::    ov2in, ov2mid, ov2out
00023 ! "ov2in" is f_{ab} \omega^a \omega^b.
00024   real(long)     ::    termin, termmid, termout
00025   real(long)     ::    dtermdoin, dtermdomid, dtermdoout
00026   real(long)     ::    dtermdrin, dtermdrmid, dtermdrout      
00027   real(long)     ::    vphiyin, vphiymid, vphiyout
00028   real(long)     ::    loghc, loghs
00029   real(long)     ::    radiold, berold, omeold
00030   real(long)     ::    ddradi, ddber, ddome
00031   real(long)     ::    facfac, numero, error, hc, hs, pre, rho, ene
00032   real(long)     ::    dummy
00033   integer        ::    nin
00034 !##############################################
00035   real(long)     ::    ut, hh, emdtest
00036 !##############################################
00037   nin = nrf_deform
00038 !
00039   sgg = 0.0d0
00040   gg  = 0.0d0
00041 !
00042   call quark_rho2phenedpdrho(rhoc_qs, pre, hc, ene, dummy)
00043   loghc = dlog(hc)
00044   call quark_rho2phenedpdrho(rhos_qs, pre, hs, ene, dummy)
00045   loghs = dlog(hs)
00046   vphiyin  = vec_phig(nin,0,0,2)
00047   vphiymid = vec_phig(0,0,0,2)
00048 !out  vphiyout = vec_phig(nrf,ntgxy,0,2)
00049   omeold = ome
00050   berold = ber
00051   radiold = radi
00052   numero = 1
00053 ! ################################################################ 
00054   ain  = dlog(alph(nin,0,0)**(1.d0/radi**2))
00055   pain = (psi(nin,0,0)**2/alph(nin,0,0))**(1.d0/radi**2)
00056   bin  = bvyd(nin,0,0)
00057 !
00058   amid = dlog(alph(0,0,0)**(1.d0/radi**2))
00059   pamid= (psi(0,0,0)**2/alph(0,0,0))**(1.d0/radi**2)
00060   bmid = bvyd(0,0,0)
00061 !
00062 !out  aout = dlog(alph(nrf,ntgxy,0)**(1.d0/radi**2))
00063 !out  paout= (psi(nrf,ntgxy,0)**2/alph(nrf,ntgxy,0))**(1.d0/radi**2)
00064 !out  bout = bvyd(nrf,ntgxy,0)
00065 ! write (6,*) "test values"
00066 ! write (6,'(1p,3e12.4)') ome, ber, radi
00067 ! write (6,'(1p,3e12.4)') psi(nin,0,0),psi(0,0,0)
00068 ! write (6,'(1p,3e12.4)') alph(nin,0,0),alph(0,0,0)
00069 ! write (6,'(1p,3e12.4)') error
00070 ! write (6,'(1p,3e12.4)') pain, pamid
00071 !! ################################################################ 
00072  ome = 0.0d0
00073   do
00074     ovin  = bin  + ome*vphiyin
00075     ovmid = bmid + ome*vphiymid
00076 !out    ovout = bout + ome*vphiyout
00077 ! This is in fact only the \omega^y component. 
00078 ! Since only \phi^y is not 0. (We are at x axis.) 
00079 ! This is also the reason why I have those bvxd down there.
00080 !      ov2in  = ovin**2 +bvxd(nin,0,0)**2+bvzd(nin,0,0)**2
00081 !      ov2mid = ovmid**2+bvxd(0,0,0)**2  +bvzd(0,0,0)**2
00082 !      ov2out = ovout**2+bvxd(nrf,ntgxy,0)**2+bvzd(nrf,ntgxy,0)**2
00083     ov2in  = ovin**2 
00084     ov2mid = ovmid**2
00085 !out    ov2out = ovout**2
00086     termin  = 1.0d0 - pain**( 2.d0*radi**2)*ov2in
00087     termmid = 1.0d0 - pamid**(2.d0*radi**2)*ov2mid
00088 !out write(6,*)termin, termmid
00089 !outstop
00090 !out    termout = 1.0d0 - paout**(2.d0*radi**2)*ov2out
00091 ! termin is 1-(\phi^4/\alpha^2)^((R_0/R_i)^2) \omega^2            
00092 !out    dtermdoin = - pain**(2.d0*radi**2)*2.0d0*ovin*vphiyin
00093 !out    dtermdomid= - pamid**(2.d0*radi**2)*2.0d0*ovmid*vphiymid
00094 !out    dtermdoout= - paout**(2.d0*radi**2)*2.0d0*ovout*vphiyout
00095     dtermdrin =-dlog(pain)*4.d0*radi*pain**(2.d0*radi**2)*ov2in
00096     dtermdrmid=-dlog(pamid)*4.d0*radi*pamid**(2.d0*radi**2)*ov2mid
00097 !outwrite(6,*)    dtermdrin ,     dtermdrmid
00098 !outstop
00099 !out    dtermdrout=-dlog(paout)*4.d0*radi*paout**(2.d0*radi**2)*ov2out
00100 !
00101 !out    sgg(1) = -(radi**2*ain  + 0.5d0*dlog(termin)  - dlog(ber))
00102 !out    sgg(2) = -(radi**2*aout + 0.5d0*dlog(termout) - dlog(ber))
00103 !out    sgg(3) = -(radi**2*amid + 0.5d0*dlog(termmid) - dlog(ber) &
00104 !out             + loghc)
00105     sgg(1) = -(radi**2*ain  + 0.5d0*dlog(termin)  - dlog(ber) + loghs)
00106     sgg(2) = -(radi**2*amid + 0.5d0*dlog(termmid) - dlog(ber) + loghc)
00107 !
00108 !  ain is ln(\alpha)
00109 !out    gg(1,1) =   0.5d0*dtermdoin/termin
00110 !out    gg(1,2) = - 1.0d0/ber
00111 !out    gg(1,3) = 2.0d0*radi*ain + 0.5d0*dtermdrin/termin
00112     gg(1,1) = - 1.0d0/ber
00113     gg(1,2) = 2.0d0*radi*ain + 0.5d0*dtermdrin/termin
00114 !
00115 !out    gg(2,1) =   0.5d0*dtermdoout/termout
00116 !out    gg(2,2) = - 1.0d0/ber
00117 !out    gg(2,3) = 2.0d0*radi*aout + 0.5d0*dtermdrout/termout
00118 !
00119 !out    gg(3,1) =   0.5d0*dtermdomid/termmid
00120 !out    gg(3,2) = - 1.0d0/ber
00121 !out    gg(3,3) = 2.0d0*radi*amid + 0.5d0*dtermdrmid/termmid
00122     gg(2,1) = - 1.0d0/ber
00123     gg(2,2) = 2.0d0*radi*amid + 0.5d0*dtermdrmid/termmid
00124 !
00125 !    test
00126 !if (numero.le.10) then
00127 ! write (6,'(1p,3e12.4)') sgg(1), sgg(2)
00128 ! write (6,'(1p,3e12.4)') gg(1,1), gg(1,2)
00129 ! write (6,'(1p,3e12.4)') gg(2,1), gg(2,2)
00130 !end if
00131 !if (numero.eq.10) stop
00132 !!stop
00133 !!pause
00134 !    test
00135 !out    call minv(gg,sgg,3,4)
00136 !out    ddome = sgg(1)
00137 !out    ddber = sgg(2)
00138 !out    ddradi = sgg(3)
00139     call minv(gg,sgg,2,4)
00140     ddber = sgg(1)
00141     ddradi = sgg(2)
00142     facfac = dmin1(dble(numero)/5.0d0,1.0d0)
00143 !out    ome = ome + ddome*facfac
00144     ome = 0.0d0
00145     ber = ber + ddber*facfac
00146     radi = radi + ddradi*facfac
00147 !out    error = dmax1(dabs(ddome/ome),dabs(ddber/ber),dabs(ddradi/radi))
00148     error = dmax1(dabs(ddber/ber),dabs(ddradi/radi))
00149     numero = numero + 1
00150     if (numero>1000) then
00151       write(6,*)' numero = ', numero, '   error =',error
00152     end if
00153 !      if (iter<10.and.numero>10) exit
00154     if (numero>1010) exit
00155     if (dabs(error)<1.0d-08) exit
00156   end do
00157 !    test
00158 ! write (6,*) "test values"
00159 ! write (6,*) sgg(1), sgg(2), sgg(3)
00160 ! write (6,*) ome, ber, radi
00161 !    test
00162 !      write(6,*)' numero = ', numero, '  error = ',error
00163 !
00164   if (radi < 0.d0) write(6,*) ' ### radi minus ###'
00165   if (radi < 0.d0) radi = - radi
00166   if (ome < 0.d0) write(6,*) ' ### ome minus ###'
00167   if (ome < 0.d0) ome = - ome
00168   if (ber < 0.d0) write(6,*) ' ### ber minus ###'
00169   if (ber < 0.d0) ber = - ber
00170 !out  ome = convf*ome + (1.d0-convf)*omeold
00171   ber = convf*ber + (1.d0-convf)*berold
00172   radi = convf*radi + (1.d0-convf)*radiold
00173 !
00174 !
00175 ! --  Improving alpha and psi.  
00176 !
00177   alph = alph**((radi/radiold)**2)
00178   psi = psi**((radi/radiold)**2)
00179 !##############################################
00180   ut = 1/sqrt(alph(nrf,ntgxy,0)**2 & 
00181      -psi(nrf,ntgxy,0)**4*ov2out)
00182   hh = ber*ut
00183   emdtest = 1.0d0/(pinx+1.0d0)*(hh-1.0d0)
00184 !write (6,*) "h  = ", hh
00185 !
00186 end subroutine update_parameter_spherical_qeos
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/update_parameter_spherical_qeos.f90
