/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/calc_omega_drot.f90

Go to the documentation of this file.

subroutine calc_omega_drot(vphiy,alphw,psiw,bvxdw,bvydw,omega,Jome,Jome_int)
  use phys_constant, only  : long
  use def_matter_parameter
!  use def_vector_phi, only : vec_phi
  implicit none
  real(long), intent(in):: vphiy, alphw, psiw, bvxdw, bvydw
  real(long), intent(inout):: omega
  real(long), intent(out):: Jome, Jome_int
! -- Jome: functional of Omega (specific angular momentum in Newtonian)
! -- Jome_int: integration of Jome w.r.t. omega
!
  real(long) :: sgg, gg
  real(long) :: p4a2, ovyu, ov2, ovphi, term1, term2
  real(long) :: dJomedo, dterm1do, dterm2do
  real(long) :: omegaold, ddomega 
  real(long) :: facfac, numero, error
!-------------------------------------------------------------
! Assyme axisymmetry.  
! Compute only on xz plane then copy to the other phi=const planes.
!
  sgg = 0.0d0
  gg  = 0.0d0
!
  p4a2 = psiw**4/alphw**2
  numero = 1
  do 
    Jome = A2DR*omega*((ome/omega)**index_DR - 1.0d0)
    ovyu = bvydw + omega*vphiy
! This is in fact only the \omega^y component.
!      ov2 = ovyu**2+bvxu(nrf,ntgxy,0)**2+bvzu(nrf,ntgxy,0)**2
    ov2   = ovyu**2
    ovphi = ovyu*vphiy
    term1 = 1.0d0 - p4a2*ov2
    term2 = p4a2*ovphi
!
    dJomedo  = - A2DR + A2DR*(ome/omega)**index_DR*dble(-index_DR+1)
    dterm1do = - p4a2*2.0d0*ovyu*vphiy
    dterm2do =   p4a2*vphiy*vphiy
!
    sgg = -(Jome*term1  - term2)
    gg  = dJomedo*term1 + Jome*dterm1do - dterm2do
!
    ddomega = sgg/gg
    facfac = dmin1(dble(numero)/5.0d0,1.0d0)
    omega = omega + ddomega*facfac
    error = dabs(ddomega/ome)
    numero = numero + 1
    if (numero>1000) then
      write(6,*)' numero = ', numero, '   error =',error
    end if
!      if (iter<10.and.numero>10) exit
    if (numero>1010) exit
    if (dabs(error)<1.0d-09) exit
!!write(6,*) Jome, ovyu, omega
!!write(6,*) Jome1, dJome1do
!!write(6,*) Jome2, dJome2do
!!write(6,*) omega,ddomega,sgg
!!pause
  end do
!!  end if
!--- computing JOmg and JOmg_intg
!
!  If (index_DR==2.00d0) then ! "v-const."-law
!     JOmg_intg = A2DR * ome**index_DR * log(abs(Omgfluid))&
!          &   -5.00d-1*A2DR*Omgfluid**2&
!          & - A2DR*ome**index_DR * log(abs(ome))&
!          &   +5.00d-1*A2DR*ome**2 
!  Else
    Jome = A2DR*omega*((ome/omega)**index_DR - 1.00d0)
    Jome_int = A2DR*ome**index_DR * omega**(2.00d0-index_DR)&
          &                  /(2.00d0-index_DR)&
          &              -5.00d-1*A2DR*omega**2&
          &              -index_DR*A2DR*ome**2/(4.00d0-2.00d0*index_DR)
!  EndIf
!
end subroutine calc_omega_drot

---