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 !
00106 !  ain is ln(\alpha)
00107 !out    gg(1,1) =   0.5d0*dtermdoin/termin
00108 !out    gg(1,2) = - 1.0d0/ber
00109 !out    gg(1,3) = 2.0d0*radi*ain + 0.5d0*dtermdrin/termin
00110     gg(1,1) = - 1.0d0/ber
00111     gg(1,2) = 2.0d0*radi*ain + 0.5d0*dtermdrin/termin
00112 !
00113 !out    gg(2,1) =   0.5d0*dtermdoout/termout
00114 !out    gg(2,2) = - 1.0d0/ber
00115 !out    gg(2,3) = 2.0d0*radi*aout + 0.5d0*dtermdrout/termout
00116 !
00117 !out    gg(3,1) =   0.5d0*dtermdomid/termmid
00118 !out    gg(3,2) = - 1.0d0/ber
00119 !out    gg(3,3) = 2.0d0*radi*amid + 0.5d0*dtermdrmid/termmid
00120     gg(2,1) = - 1.0d0/ber
00121     gg(2,2) = 2.0d0*radi*amid + 0.5d0*dtermdrmid/termmid
00122 !
00123 !    test
00124 !if (numero.le.10) then
00125 ! write (6,'(1p,3e12.4)') sgg(1), sgg(2)
00126 ! write (6,'(1p,3e12.4)') gg(1,1), gg(1,2)
00127 ! write (6,'(1p,3e12.4)') gg(2,1), gg(2,2)
00128 !end if
00129 !if (numero.eq.10) stop
00130 !!stop
00131 !!pause
00132 !    test
00133 !out    call minv(gg,sgg,3,4)
00134 !out    ddome = sgg(1)
00135 !out    ddber = sgg(2)
00136 !out    ddradi = sgg(3)
00137     call minv(gg,sgg,2,4)
00138     ddber = sgg(1)
00139     ddradi = sgg(2)
00140     facfac = dmin1(dble(numero)/5.0d0,1.0d0)
00141 !out    ome = ome + ddome*facfac
00142     ome = 0.0d0
00143     ber = ber + ddber*facfac
00144     radi = radi + ddradi*facfac
00145 !out    error = dmax1(dabs(ddome/ome),dabs(ddber/ber),dabs(ddradi/radi))
00146     error = dmax1(dabs(ddber/ber),dabs(ddradi/radi))
00147     numero = numero + 1
00148     if (numero>1000) then
00149       write(6,*)' numero = ', numero, '   error =',error
00150     end if
00151 !      if (iter<10.and.numero>10) exit
00152     if (numero>1010) exit
00153     if (dabs(error)<1.0d-08) exit
00154   end do
00155 !    test
00156 ! write (6,*) "test values"
00157 ! write (6,*) sgg(1), sgg(2), sgg(3)
00158 ! write (6,*) ome, ber, radi
00159 !    test
00160 !      write(6,*)' numero = ', numero, '  error = ',error
00161 !
00162   if (radi < 0.d0) write(6,*) ' ### radi minus ###'
00163   if (radi < 0.d0) radi = - radi
00164   if (ome < 0.d0) write(6,*) ' ### ome minus ###'
00165   if (ome < 0.d0) ome = - ome
00166   if (ber < 0.d0) write(6,*) ' ### ber minus ###'
00167   if (ber < 0.d0) ber = - ber
00168 !out  ome = convf*ome + (1.d0-convf)*omeold
00169   ber = convf*ber + (1.d0-convf)*berold
00170   radi = convf*radi + (1.d0-convf)*radiold
00171 !
00172 !
00173 ! --  Improving alpha and psi.  
00174 !
00175   alph = alph**((radi/radiold)**2)
00176   psi = psi**((radi/radiold)**2)
00177 !##############################################
00178   ut = 1/sqrt(alph(nrf,ntgxy,0)**2 & 
00179      -psi(nrf,ntgxy,0)**4*ov2out)
00180   hh = ber*ut
00181   emdtest = 1.0d0/(pinx+1.0d0)*(hh-1.0d0)
00182 !write (6,*) "h  = ", hh
00183 !
00184 end subroutine update_parameter_spherical_peos
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/update_parameter_spherical_peos.f90
