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

00001 subroutine update_parameter_spherical_peos(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
00029   real(long)     ::    radiold, berold, omeold
00030   real(long)     ::    ddradi, ddber, ddome
00031   real(long)     ::    facfac, numero, error, hc, pre, rho, ene
00032   integer        ::    nin
00033 !##############################################
00034   real(long)     ::    ut, hh, emdtest
00035 !##############################################
00036   nin = nrf_deform
00037 !
00038   sgg = 0.0d0
00039   gg  = 0.0d0
00040 !
00041   call peos_q2hprho(emdc, hc, pre, rho, ene)
00042   loghc = dlog(hc)
00043   vphiyin  = vec_phig(nin,0,0,2)
00044   vphiymid = vec_phig(0,0,0,2)
00045 !out  vphiyout = vec_phig(nrf,ntgxy,0,2)
00046   omeold = ome
00047   berold = ber
00048   radiold = radi
00049   numero = 1
00050 ! ################################################################ 
00051   ain  = dlog(alph(nin,0,0)**(1.d0/radi**2))
00052   pain = (psi(nin,0,0)**2/alph(nin,0,0))**(1.d0/radi**2)
00053   bin  = bvyd(nin,0,0)
00054 !
00055   amid = dlog(alph(0,0,0)**(1.d0/radi**2))
00056   pamid= (psi(0,0,0)**2/alph(0,0,0))**(1.d0/radi**2)
00057   bmid = bvyd(0,0,0)
00058 !
00059 !out  aout = dlog(alph(nrf,ntgxy,0)**(1.d0/radi**2))
00060 !out  paout= (psi(nrf,ntgxy,0)**2/alph(nrf,ntgxy,0))**(1.d0/radi**2)
00061 !out  bout = bvyd(nrf,ntgxy,0)
00062 ! write (6,*) "test values"
00063 ! write (6,'(1p,3e12.4)') ome, ber, radi
00064 ! write (6,'(1p,3e12.4)') psi(nin,0,0),psi(0,0,0)
00065 ! write (6,'(1p,3e12.4)') alph(nin,0,0),alph(0,0,0)
00066 ! write (6,'(1p,3e12.4)') error
00067 ! write (6,'(1p,3e12.4)') pain, pamid
00068 !! ################################################################ 
00069  ome = 0.0d0
00070   do
00071     ovin  = bin  + ome*vphiyin
00072     ovmid = bmid + ome*vphiymid
00073 !out    ovout = bout + ome*vphiyout
00074 ! This is in fact only the \omega^y component. 
00075 ! Since only \phi^y is not 0. (We are at x axis.) 
00076 ! This is also the reason why I have those bvxd down there.
00077 !      ov2in  = ovin**2 +bvxd(nin,0,0)**2+bvzd(nin,0,0)**2
00078 !      ov2mid = ovmid**2+bvxd(0,0,0)**2  +bvzd(0,0,0)**2
00079 !      ov2out = ovout**2+bvxd(nrf,ntgxy,0)**2+bvzd(nrf,ntgxy,0)**2
00080     ov2in  = ovin**2 
00081     ov2mid = ovmid**2
00082 !out    ov2out = ovout**2
00083     termin  = 1.0d0 - pain**( 2.d0*radi**2)*ov2in
00084     termmid = 1.0d0 - pamid**(2.d0*radi**2)*ov2mid
00085 !out write(6,*)termin, termmid
00086 !outstop
00087 !out    termout = 1.0d0 - paout**(2.d0*radi**2)*ov2out
00088 ! termin is 1-(\phi^4/\alpha^2)^((R_0/R_i)^2) \omega^2            
00089 !out    dtermdoin = - pain**(2.d0*radi**2)*2.0d0*ovin*vphiyin
00090 !out    dtermdomid= - pamid**(2.d0*radi**2)*2.0d0*ovmid*vphiymid
00091 !out    dtermdoout= - paout**(2.d0*radi**2)*2.0d0*ovout*vphiyout
00092     dtermdrin =-dlog(pain)*4.d0*radi*pain**(2.d0*radi**2)*ov2in
00093     dtermdrmid=-dlog(pamid)*4.d0*radi*pamid**(2.d0*radi**2)*ov2mid
00094 !outwrite(6,*)    dtermdrin ,     dtermdrmid
00095 !outstop
00096 !out    dtermdrout=-dlog(paout)*4.d0*radi*paout**(2.d0*radi**2)*ov2out
00097 !
00098 !out    sgg(1) = -(radi**2*ain  + 0.5d0*dlog(termin)  - dlog(ber))
00099 !out    sgg(2) = -(radi**2*aout + 0.5d0*dlog(termout) - dlog(ber))
00100 !out    sgg(3) = -(radi**2*amid + 0.5d0*dlog(termmid) - dlog(ber) &
00101 !out             + loghc)
00102 !    sgg(1) = -(radi**2*ain  + 0.5d0*dlog(termin)  - dlog(ber))
00103 !    sgg(2) = -(radi**2*amid + 0.5d0*dlog(termmid) - dlog(ber) &
00104 !    &        + loghc)
00105     sgg(1) = -(radi**2*ain  - dlog(ber))
00106     sgg(2) = -(radi**2*amid - 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     gg(1,1) = - 1.0d0/ber
00115     gg(1,2) = 2.0d0*radi*ain
00116 !
00117 !out    gg(2,1) =   0.5d0*dtermdoout/termout
00118 !out    gg(2,2) = - 1.0d0/ber
00119 !out    gg(2,3) = 2.0d0*radi*aout + 0.5d0*dtermdrout/termout
00120 !
00121 !out    gg(3,1) =   0.5d0*dtermdomid/termmid
00122 !out    gg(3,2) = - 1.0d0/ber
00123 !out    gg(3,3) = 2.0d0*radi*amid + 0.5d0*dtermdrmid/termmid
00124 !    gg(2,1) = - 1.0d0/ber
00125 !    gg(2,2) = 2.0d0*radi*amid + 0.5d0*dtermdrmid/termmid
00126     gg(2,1) = - 1.0d0/ber
00127     gg(2,2) = 2.0d0*radi*amid
00128 !
00129 !    test
00130 !if (numero.le.10) then
00131 ! write (6,'(1p,3e12.4)') sgg(1), sgg(2)
00132 ! write (6,'(1p,3e12.4)') gg(1,1), gg(1,2)
00133 ! write (6,'(1p,3e12.4)') gg(2,1), gg(2,2)
00134 !end if
00135 !if (numero.eq.10) stop
00136 !!stop
00137 !!pause
00138 !    test
00139 !out    call minv(gg,sgg,3,4)
00140 !out    ddome = sgg(1)
00141 !out    ddber = sgg(2)
00142 !out    ddradi = sgg(3)
00143     call minv(gg,sgg,2,4)
00144     ddber = sgg(1)
00145     ddradi = sgg(2)
00146     facfac = dmin1(dble(numero)/5.0d0,1.0d0)
00147 !out    ome = ome + ddome*facfac
00148     ome = 0.0d0
00149     ber = ber + ddber*facfac
00150     radi = radi + ddradi*facfac
00151 !out    error = dmax1(dabs(ddome/ome),dabs(ddber/ber),dabs(ddradi/radi))
00152     error = dmax1(dabs(ddber/ber),dabs(ddradi/radi))
00153     numero = numero + 1
00154     if (numero>1000) then
00155       write(6,*)' numero = ', numero, '   error =',error
00156     end if
00157 !      if (iter<10.and.numero>10) exit
00158     if (numero>1010) exit
00159     if (dabs(error)<1.0d-08) exit
00160   end do
00161 !    test
00162 ! write (6,*) "test values"
00163 ! write (6,*) sgg(1), sgg(2), sgg(3)
00164 ! write (6,*) ome, ber, radi
00165 !    test
00166 !      write(6,*)' numero = ', numero, '  error = ',error
00167 !
00168   if (radi < 0.d0) write(6,*) ' ### radi minus ###'
00169   if (radi < 0.d0) radi = - radi
00170   if (ome < 0.d0) write(6,*) ' ### ome minus ###'
00171   if (ome < 0.d0) ome = - ome
00172   if (ber < 0.d0) write(6,*) ' ### ber minus ###'
00173   if (ber < 0.d0) ber = - ber
00174 !out  ome = convf*ome + (1.d0-convf)*omeold
00175   ber = convf*ber + (1.d0-convf)*berold
00176   radi = convf*radi + (1.d0-convf)*radiold
00177 !
00178 !
00179 ! --  Improving alpha and psi.  
00180 !
00181   alph = alph**((radi/radiold)**2)
00182   psi = psi**((radi/radiold)**2)
00183 !##############################################
00184   ut = 1/sqrt(alph(nrf,ntgxy,0)**2 & 
00185      -psi(nrf,ntgxy,0)**4*ov2out)
00186   hh = ber*ut
00187   emdtest = 1.0d0/(pinx+1.0d0)*(hh-1.0d0)
00188 !write (6,*) "h  = ", hh
00189 !
00190 end subroutine update_parameter_spherical_peos
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/update_parameter_spherical_peos.f90
