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

00001 subroutine update_parameter_axisym_peos_drot(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, bvyd, tfkij, tfkijkij
00005   use def_matter, only : omef
00006   use def_matter_parameter
00007   use coordinate_grav_r, only : rg
00008   use trigonometry_grav_theta, only : sinthg, costhg
00009   use trigonometry_grav_phi, only   : sinphig, cosphig
00010   use def_vector_phi, only : vec_phig
00011   use interface_minv
00012   implicit none
00013   real(long), intent(in) :: convf
00014   real(long)     ::    sgg(1:4), gg(1:4,1:4)
00015 ! "bin" is beta(in) and so on.
00016   real(long)     ::    bin, bmid, bout
00017 ! ain is log(\alpha^((1/R_i)^2)) (in). R_in is in fact the previous R.
00018   real(long)     ::    ain, amid, aout
00019 ! pain = (\phi^2/\alpha)^((1/R_i)^2) (in).
00020   real(long)     ::    pain, pamid, paout
00021 ! ovin is \omega(in). \omega = \beta + \Omega \phi.
00022   real(long)     ::    ovin, ovmid, ovout
00023 ! "ov2in" is f_{ab} \omega^a \omega^b.
00024   real(long)     ::    ov2in, ov2mid, ov2out
00025 ! "jin" is Integral of J(Omega).
00026   real(long)     ::    jin, jmid, jout, jint_in, jint_mid, jint_out
00027   real(long)     ::    alin, almid, alout
00028   real(long)     ::    psin, psmid, psout
00029   real(long)     ::    hyin, hymid, hyout
00030   real(long)     ::    bxin, bxmid, bxout
00031   real(long)     ::    omein, omemid, omeout
00032 !
00033   real(long)     ::    termin, termmid, termout
00034   real(long)     ::    dtermdoin, dtermdomid, dtermdoout
00035   real(long)     ::    dtermdrin, dtermdrmid, dtermdrout      
00036   real(long)     ::    djintdoin, djintdomid, djintdoout
00037   real(long)     ::    dodoc_in, dodoc_mid, dodoc_out
00038   real(long)     ::    vphiyin, vphiymid, vphiyout
00039   real(long)     ::    loghc
00040   real(long)     ::    radiold, berold, omeold
00041   real(long)     ::    ddradi, ddber, ddome
00042   real(long)     ::    facfac, numero, error, hc, pre, rho, ene
00043   real(long)     ::    ome_mx, ome_eq, Jome, Jome_int, dum, small = 1.0d-30
00044   integer        ::    nin, nmx
00045 !##############################################
00046   real(long)     ::    ut, hh, emdtest
00047 !##############################################
00048   nin = nrf_deform
00049 !
00050   sgg = 0.0d0
00051   gg  = 0.0d0
00052 !
00053   call peos_q2hprho(emdc, hc, pre, rho, ene)
00054   call search_emdmax_xaxis_grid(nmx)
00055 !  nmx = 0
00056 !testtest
00057   loghc = dlog(hc)
00058   vphiyin  = vec_phig(nin,0,0,2)
00059   vphiymid = vec_phig(nmx,ntgxy,0,2)
00060   vphiyout = vec_phig(nrf,ntgxy,0,2)
00061   omeold = ome
00062   berold = ber
00063   radiold = radi
00064   numero = 1
00065 ! ################################################################ 
00066   ain  = dlog(alph(nin,0,0)**(1.d0/radi**2))
00067   pain = (psi(nin,0,0)**2/alph(nin,0,0))**(1.d0/radi**2)
00068   bin  = bvyd(nin,0,0)
00069   hyin = 0.0d0
00070   psin = psi(nin,0,0)
00071   alin = alph(nin,0,0)
00072   omein = omef(nrf,0,0)
00073 !
00074   amid = dlog(alph(nmx,ntgxy,0)**(1.d0/radi**2))
00075   pamid= (psi(nmx,ntgxy,0)**2/alph(nmx,ntgxy,0))**(1.d0/radi**2)
00076   bmid = bvyd(nmx,ntgxy,0)
00077   hymid= 0.0d0
00078   psmid= psi(nmx,ntgxy,0)
00079   almid= alph(nmx,ntgxy,0)
00080   omemid=omef(nmx,ntgxy,0)
00081 !
00082   aout = dlog(alph(nrf,ntgxy,0)**(1.d0/radi**2))
00083   paout= (psi(nrf,ntgxy,0)**2/alph(nrf,ntgxy,0))**(1.d0/radi**2)
00084   bout = bvyd(nrf,ntgxy,0)
00085   hyout= 0.0d0
00086   psout= psi(nrf,ntgxy,0)
00087   alout= alph(nrf,ntgxy,0)
00088   omeout=omef(nrf,ntgxy,0)
00089 !
00090   ome_mx = omemid
00091   ome_eq = omeout
00092 !
00093 !  write (6,*) "test values"
00094 !  write (6,'(1p,3e12.4)') ome, ber, radi
00095 ! ! write (6,'(1p,3e12.4)') psi(nin,0,0),psi(nmx,ntgxy,0),psi(nrf,ntgxy,0)
00096 ! ! write (6,'(1p,3e12.4)') alph(nin,0,0),alph(nmx,ntgxy,0),alph(nrf,ntgxy,0)
00097 ! ! write (6,'(1p,3e12.4)') error
00098 !  write (6,'(1p,3e12.4)') vphiyin, vphiymid, vphiyout
00099 !  write (6,'(1p,3e12.4)') alin, almid, alout
00100 !  write (6,'(1p,3e12.4)') psin, psmid, psout
00101 !  write (6,'(1p,3e12.4)') bin, bmid, bout
00102 !  write (6,'(1p,3e12.4)') omein, omemid, omeout
00103 !  write (6,'(1p,3e12.4)') jin, jmid, jout
00104 !  write (6,'(1p,3e12.4)') jint_in, jint_mid, jint_out
00105 !  write (6,'(1p,3e12.4)') ome_mx, ome_eq
00106 !! stop
00107 ! ################################################################ 
00108   do
00109     call calc_omega_drot(vphiyin,alin,psin,bin,hyin,omein,Jome,Jome_int)
00110     jin  = Jome ; jint_in  = Jome_int
00111     call calc_omega_drot(vphiymid,almid,psmid,bmid,hymid,omemid,Jome,Jome_int)
00112     jmid = Jome ; jint_mid = Jome_int
00113     call calc_omega_drot(vphiyout,alout,psout,bout,hyout,omeout,Jome,Jome_int)
00114     jout = Jome ; jint_out = Jome_int
00115 !
00116     ovin  = bin  + omein*vphiyin
00117     ovmid = bmid + omemid*vphiymid
00118     ovout = bout + omeout*vphiyout
00119 ! This is in fact only the \omega^y component. 
00120 ! Since only \phi^y is not 0. (We are at x axis.) 
00121 ! This is also the reason why I have those bvxd down there.
00122 !      ov2in  = ovin**2 +bvxd(nin,0,0)**2+bvzd(nin,0,0)**2
00123 !      ov2mid = ovmid**2+bvxd(nmx,ntgxy,0)**2+bvzd(nmx,ntgxy,0)**2
00124 !      ov2out = ovout**2+bvxd(nrf,ntgxy,0)**2+bvzd(nrf,ntgxy,0)**2
00125     ov2in  = ovin**2 
00126     ov2mid = ovmid**2
00127     ov2out = ovout**2
00128     termin  = 1.0d0 - pain**( 2.d0*radi**2)*ov2in
00129     termmid = 1.0d0 - pamid**(2.d0*radi**2)*ov2mid
00130     termout = 1.0d0 - paout**(2.d0*radi**2)*ov2out
00131 !
00132 !    dodoc_in  = A2DR*index_DR*termin &
00133 !    &         /(pain**(2.d0*radi**2)*vphiyin**2 &
00134 !    &         + A2DR*index_DR*termin &
00135 !    &         + 2.0d0*pain**(2.d0*radi**2)*vphiyin*ome**2*ovin*jin &
00136 !    &         - ome*termin*jin)
00137 !    dodoc_mid = A2DR*index_DR*termmid*(ome/(omemid+small))**(index_DR-1.0d0) &
00138 !    &         /(pamid**(2.d0*radi**2)*vphiymid**2 &
00139 !    &         + A2DR*index_DR*termmid*(ome/(omemid+small))**index_DR &
00140 !    &         + 2.0d0*pamid**(2.d0*radi**2)*vphiymid*omemid**2*ovmid*jmid &
00141 !    &         - omemid*termmid*jmid)
00142 !    dodoc_out = A2DR*index_DR*termout*(ome/(omeout+small))**(index_DR-1.0d0) &
00143 !    &         /(paout**(2.d0*radi**2)*vphiyout**2 &
00144 !    &         + A2DR*index_DR*termout*(ome/(omeout+small))**index_DR &
00145 !    &         + 2.0d0*paout**(2.d0*radi**2)*vphiyout*omeout**2*ovout*jout &
00146 !    &         - omeout*termout*jout)
00147 !
00148     dodoc_in  = -((A2DR*index_DR*omein*termin) &
00149     &         /(ome*(-(A2DR*index_DR*termin) &
00150     &         +(jin*termin)/omein  &
00151     &         - 2.0d0*jin*ovin*pain*vphiyin - pain*vphiyin**2)))
00152     dodoc_mid = -((A2DR*index_DR*(ome/omemid)**index_DR*omemid*termmid) &
00153     &         /(ome*(-(A2DR*index_DR*(ome/omemid)**index_DR*termmid) &
00154     &         +(jmid*termmid)/omemid &
00155     &         - 2.0d0*jmid*ovmid*pamid*vphiymid - pamid*vphiymid**2)))
00156     dodoc_out = -((A2DR*index_DR*(ome/omeout)**index_DR*omeout*termout) &
00157     &         /(ome*(-(A2DR*index_DR*(ome/omeout)**index_DR*termout) &
00158     &         +(jout*termout)/omeout &
00159     &         - 2.0d0*jout*ovout*paout*vphiyout - paout*vphiyout**2)))
00160 !
00161 !!write(6,*) 'dodoc'
00162 !!write(6,*) ome,omemid,omeout
00163 !!write(6,*) dodoc_in, dodoc_mid, dodoc_out
00164 ! termin is 1-(\phi^4/\alpha^2)^((R_0/R_i)^2) \omega^2            
00165   if (index_DR.ne.2.0d0) then
00166     djintdoin  = 0.0d0  !+ jin *dodoc_in
00167     djintdomid = A2DR*index_DR*omemid*(ome/(omemid+small))**(index_DR-1.0d0) &
00168     &                             /(2.0d0-index_DR) &
00169     &          - A2DR*index_DR*ome/(2.0d0-index_DR) + jmid*dodoc_mid
00170     djintdoout = A2DR*index_DR*omeout*(ome/(omeout+small))**(index_DR-1.0d0) &
00171     &                             /(2.0d0-index_DR) &
00172     &          - A2DR*index_DR*ome/(2.0d0-index_DR) + jout*dodoc_out
00173   else 
00174     djintdoin  = 0.0d0 ! + jin*dodoc_in
00175     djintdomid = - 2.0d0*A2DR*ome*dlog(ome/(omemid+small)) + jmid*dodoc_mid
00176     djintdoout = - 2.0d0*A2DR*ome*dlog(ome/(omeout+small)) + jout*dodoc_out
00177   end if
00178 !
00179     dtermdoin = - pain**(2.d0*radi**2)*2.0d0*ovin*vphiyin*dodoc_in
00180     dtermdomid= - pamid**(2.d0*radi**2)*2.0d0*ovmid*vphiymid*dodoc_mid
00181     dtermdoout= - paout**(2.d0*radi**2)*2.0d0*ovout*vphiyout*dodoc_out
00182     dtermdrin =-dlog(pain)*4.d0*radi*pain**(2.d0*radi**2)*ov2in
00183     dtermdrmid=-dlog(pamid)*4.d0*radi*pamid**(2.d0*radi**2)*ov2mid
00184     dtermdrout=-dlog(paout)*4.d0*radi*paout**(2.d0*radi**2)*ov2out
00185 !
00186     sgg(1) = -(radi**2*ain  + 0.5d0*dlog(termin)  + jint_in  - dlog(ber))
00187     sgg(2) = -(radi**2*aout + 0.5d0*dlog(termout) + jint_out - dlog(ber))
00188     sgg(3) = -(radi**2*amid + 0.5d0*dlog(termmid) + jint_mid - dlog(ber) &
00189              + loghc)
00190 !
00191 !  ain is ln(\alpha)
00192     gg(1,1) =   0.5d0*dtermdoin/termin + djintdoin
00193     gg(1,2) = - 1.0d0/ber
00194     gg(1,3) = 2.0d0*radi*ain + 0.5d0*dtermdrin/termin
00195 !
00196     gg(2,1) =   0.5d0*dtermdoout/termout + djintdoout
00197     gg(2,2) = - 1.0d0/ber
00198     gg(2,3) = 2.0d0*radi*aout + 0.5d0*dtermdrout/termout
00199 !
00200     gg(3,1) =   0.5d0*dtermdomid/termmid + djintdomid
00201     gg(3,2) = - 1.0d0/ber
00202     gg(3,3) = 2.0d0*radi*amid + 0.5d0*dtermdrmid/termmid
00203 !
00204 !    test
00205 !if (numero.eq.1) then
00206 ! write (6,'(1p,3e12.4)') sgg(1), sgg(2), sgg(3)
00207 ! write (6,'(1p,3e12.4)') gg(1,1), gg(1,2), gg(1,3)
00208 ! write (6,'(1p,3e12.4)') gg(2,1), gg(2,2), gg(2,3)
00209 ! write (6,'(1p,3e12.4)') gg(3,1), gg(3,2), gg(3,3)
00210 !end if
00211 !stop
00212 !!pause
00213 !    test
00214     call minv(gg,sgg,3,4)
00215     ddome = sgg(1)
00216     ddber = sgg(2)
00217     ddradi = sgg(3)
00218     facfac = dmin1(dble(numero)/5.0d0,1.0d0)
00219     ome = ome + ddome*facfac
00220     ber = ber + ddber*facfac
00221     radi = radi + ddradi*facfac
00222     error = dmax1(dabs(ddome/ome),dabs(ddber/ber),dabs(ddradi/radi))
00223     numero = numero + 1
00224     if (numero>1000) then
00225       write(6,*)' numero = ', numero, '   error =',error
00226     end if
00227 !      if (iter<10.and.numero>10) exit
00228     if (numero>1010) exit
00229 !###    if (dabs(error)<1.0d-08) exit
00230     if (dabs(error)<1.0d-12) exit
00231   end do
00232 !    test
00233 ! write (6,*) "test values"
00234 ! write (6,*) sgg(1), sgg(2), sgg(3)
00235 ! write (6,*) 'update paramter', ome
00236 !    test
00237 !      write(6,*)' numero = ', numero, '  error = ',error
00238 !
00239   if (radi < 0.d0) write(6,*) ' ### radi minus ###'
00240   if (radi < 0.d0) radi = - radi
00241   if (ome < 0.d0) write(6,*) ' ### ome minus ###'
00242   if (ome < 0.d0) ome = - ome
00243   if (ber < 0.d0) write(6,*) ' ### ber minus ###'
00244   if (ber < 0.d0) ber = - ber
00245   ome = convf*ome + (1.d0-convf)*omeold
00246   ber = convf*ber + (1.d0-convf)*berold
00247   radi = convf*radi + (1.d0-convf)*radiold
00248 !
00249 !
00250 ! --  Improving alpha and psi.  
00251 !
00252   alph = alph**((radi/radiold)**2)
00253   psi = psi**((radi/radiold)**2)
00254 !##############################################
00255   ut = 1/sqrt(alph(nrf,ntgxy,0)**2 & 
00256      -psi(nrf,ntgxy,0)**4*ov2out)
00257   hh = ber*ut
00258   emdtest = 1.0d0/(pinx+1.0d0)*(hh-1.0d0)
00259 !!write (6,*) "h  = ", hh
00260 !
00261 end subroutine update_parameter_axisym_peos_drot
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/update_parameter_axisym_peos_drot.f90
