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

00001 subroutine calc_omega_drot_Newton(vphiy,alphw,psiw,bvydw,hyydw,omega,flag_st)
00002   use phys_constant, only  : long
00003   use coordinate_grav_r, only : rg
00004   use def_matter_parameter, only : ome, A2DR, DRAT_A2DR, index_DR
00005   implicit none
00006   real(long) :: vphiy, alphw, psiw, bvydw, hyydw, omega
00007 ! -- Jome: functional of Omega (specific angular momentum in Newtonian)
00008   real(long) :: Jome, pmfac
00009   real(long) :: sgg, gg
00010   real(long) :: p4a2, ovyu, ov2, ovphi, term1, term2, term3
00011   real(long) :: dJomedo, dterm1do, dterm2do
00012   real(long) :: ddomega
00013   real(long) :: facfac, error, small = 1.0d-30
00014   integer    :: ic, nic, flag_st
00015   real(long) :: facp(5) = (/ 1.0d-1, 2.0d-1, 4.0d-1, 7.0d-1, 1.0d-0 /)
00016 !
00017 !-------------------------------------------------------------
00018 ! Assume axisymmetry.  
00019 !
00020   if (ome  .eq.0.0d0) ome   = 0.01d0 ! For 1D initial
00021   if (omega.eq.0.0d0) omega = 0.01d0 ! For 1D initial
00022 !
00023   if (vphiy.lt.0.001d0*rg(1)) then
00024     omega = ome
00025     return
00026   end if
00027 !
00028 !rotlaw_type1  pmfac = -1.0d0
00029 !rotlaw_type2  pmfac =  1.0d0
00030 !
00031   nic = 30
00032 !
00033   do ic = 0, nic
00034 !
00035     Jome = A2DR*omega*((ome/(omega+small))**index_DR - 1.0d0)
00036     ovyu = bvydw + omega*vphiy
00037 ! This is in fact only the \omega^y component.
00038 !      ov2 = ovyu**2+ovxu**2+ovzu**2
00039     ov2   = ovyu**2*(1.0d0 + hyydw)
00040     ovphi = ovyu*vphiy*(1.0d0 + hyydw)
00041     term1 = 1.0d0 - p4a2*ov2
00042     term2 = p4a2*ovphi
00043 !
00044     dJomedo  = - A2DR + A2DR*(ome/(omega+small))**index_DR*(-index_DR+1.0d0)
00045     dterm1do = - p4a2*2.0d0*ovyu*vphiy
00046     dterm2do =   p4a2*vphiy*vphiy
00047 !
00048     sgg = -(Jome*term1  - term2)
00049     gg  = dJomedo*term1 + Jome*dterm1do - dterm2do
00050     ddomega = sgg/gg
00051     facfac = facp(min(ic+1,5))
00052     omega = omega + ddomega*facfac
00053 !
00054     if (omega.gt.ome) omega = ome
00055     if (omega.le.0.0d0) omega = small
00056     error = dabs(ddomega/omega)
00057     if (dabs(error)<1.0d-14) exit
00058 !write(6,*) Jome, ovyu, omega
00059 !write(6,*) dJomedo
00060 !!write(6,*) Jome1, dJome1do
00061 !!write(6,*) Jome2, dJome2do
00062 !write(6,*) omega,ddomega,sgg
00063 !!pause
00064   end do
00065 !
00066   flag_st = 0
00067   if (ic.ge.nic-2) flag_st = 1
00068 !!write(6,*) ome,A2DR,B2DR
00069 !!write(6,*) index_DRq,index_DRp
00070 !!write(6,*) ic, nic, error_st
00071 !
00072 !!  if (ic.eq.nic) then 
00073 !!    write(6,*) ic, omega
00074 !!    write(6,*) 'calc_omega_drot_secant.f90_typeOJ not converged'
00075 !!!    stop 'calc_omega_drot_Newton.f90_type1type2'
00076 !!  end if
00077 !
00078 end subroutine calc_omega_drot_Newton
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/calc_omega_drot_Newton.f90
