/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine/calc_omega_drot_secant.f90
Go to the documentation of this file.00001 subroutine calc_omega_drot_secant(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_DRq, index_DR, & 00005 & B2DR, DRAT_B2DR, index_DRp 00006 implicit none 00007 real(long) :: vphiy, alphw, psiw, bvydw, hyydw, omega 00008 ! -- Jome: functional of Omega (specific angular momentum in Newtonian) 00009 real(long) :: Jome, omegatmp, Jometmp, pmfac 00010 real(long) :: x_st(0:2), f_st(0:1) 00011 real(long) :: delta, dfdx 00012 real(long) :: error_st, diff_st, epsilon = 1.0d-9, small = 1.0d-30 00013 integer :: ii, ic, nic, flag_st 00014 real(long) :: facp(5) = (/ 1.0d-1, 2.0d-1, 4.0d-1, 7.0d-1, 1.0d-0 /) 00015 00016 !------------------------------------------------------------- 00017 ! Assume axisymmetry. 00018 ! 00019 if (ome .eq.0.0d0) ome = 0.01d0 ! For 1D initial 00020 if (omega.eq.0.0d0) omega = 0.01d0 ! For 1D initial 00021 ! 00022 if (vphiy.lt.0.001d0*rg(1)) then 00023 omega = ome 00024 return 00025 end if 00026 ! 00027 nic = 40 00028 ! 00029 do ic = 0, nic 00030 ! 00031 !rotlaw_type0 if(omega.gt.ome) omega = ome 00032 !rotlaw_type2 if(omega.lt.ome) omega = ome 00033 ! 00034 call calc_J_utuphi(vphiy,alphw,psiw,bvydw,hyydw,omega,Jome) 00035 !rotlaw_type0 Jometmp = A2DR*omega*((ome/(omega+small))**index_DR - 1.0d0) 00036 !rotlaw_type0 diff_st = Jome - Jometmp 00037 !rotlaw_type0 error_st= dabs(diff_st/(Jome+small)) 00038 !rotlaw_type1 pmfac = -1.0d0 00039 !rotlaw_type2 pmfac = 1.0d0 00040 !!rotlaw_type12 Jometmp = A2DR*omega*(pmfac* & 00041 !!rotlaw_type12 & (omega/ome - 1.0d0)+small)**index_DR 00042 !rotlaw_type12 Jometmp = A2DR*omega*(pmfac* & 00043 !rotlaw_type12 & (omega/ome - 1.0d0)) 00044 !rotlaw_type12 diff_st = Jome**(1.0/index_DR) - Jometmp 00045 !!rotlaw_type12 error_st= dabs(diff_st/(Jometmp+small)) 00046 !rotlaw_type12 error_st= dabs((x_st(1)-x_st(0))/(x_st(1)+small)) 00047 !rotlaw_typeOJ omegatmp = ome*(1.0d0+(Jome/(B2DR*ome))**index_DRp) & 00048 !rotlaw_typeOJ & /(1.0d0+(Jome/(A2DR*ome))**(index_DRq+index_DRp)) 00049 !rotlaw_typeOJ diff_st = omega - omegatmp 00050 !rotlaw_typeOJ error_st= dabs(diff_st/omega) 00051 !rotlaw_typeJC omegatmp = ome*(1.0d0+(Jome/(B2DR*ome))**index_DRp) & 00052 !rotlaw_typeJC & *(1.0d0-(Jome/(A2DR*ome))) 00053 !rotlaw_typeJC diff_st = omega - omegatmp 00054 !rotlaw_typeJC error_st= dabs(diff_st/omega) 00055 !! write(6,'(1a6,1p,3e12.4)')'omeome', omega, omegatmp 00056 !!! write(6,'(1a6,1p,3e12.4)')'omeome', Jome, Jometmp 00057 !!! write(6,'(1a6,1p,3e12.4)')'omeome', diff_st, error_st 00058 !!! write(6,*) A2DR,omega,pmfac 00059 !!! write(6,*) ome,small,index_DR 00060 ! 00061 if (ic.ne.0) then 00062 x_st(0) = x_st(1) 00063 f_st(0) = f_st(1) 00064 end if 00065 ! 00066 if (ic.eq.0.or.error_st.gt.epsilon) then 00067 x_st(1) = omega 00068 f_st(1) = diff_st 00069 else 00070 exit 00071 end if 00072 ! 00073 if (ic.eq.0) then 00074 delta = omega*0.1d0 00075 if (delta.eq.0.0d0) delta = 0.01d0 00076 !rotlaw_type01 delta = - omega*0.1d0 00077 !rotlaw_type01 if (delta.eq.0.0d0) delta = - 0.01d0 00078 else 00079 dfdx = (f_st(1) - f_st(0))/(x_st(1) - x_st(0)) 00080 delta = - f_st(1)/dfdx 00081 end if 00082 ! 00083 ii = min0(5,ic+1) 00084 x_st(2) = x_st(1) + delta*facp(ii) 00085 omega = x_st(2) 00086 ! 00087 !!! write(6,'(1a6,i10)')'omseca', ic 00088 !!! write(6,'(1a6,1p,3e12.4)')'omseca', x_st(0),x_st(1),x_st(2) 00089 !!! write(6,'(1a6,1p,3e12.4)')'omseca', f_st(0),f_st(1) 00090 !!! write(6,'(1a6,1p,3e12.4)')'omseca', delta, dfdx, omega 00091 !!! write(6,*)'omseca', x_st(0),x_st(1) 00092 !!! write(6,*)'omseca', f_st(0),f_st(1) 00093 !!! write(6,*)'omseca', dfdx 00094 end do 00095 ! 00096 flag_st = 0 00097 if (ic.ge.nic-2) flag_st = 1 00098 !rotlaw_type01 if (omega.lt.0.0.or.omega.gt.ome) flag_st = 1 00099 !rotlaw_type2 if (omega.lt.ome) flag_st = 1 00100 ! 00101 !!write(6,*) ome,A2DR,B2DR 00102 !!write(6,*) index_DRq,index_DRp 00103 !!write(6,*) ic, nic, error_st 00104 ! 00105 !! if (ic.eq.nic) then 00106 !! write(6,*) ic, omega 00107 !! write(6,*) 'calc_omega_drot_secant.f90_typeOJ not converged' 00108 !!! stop 'calc_omega_drot_secant.f90_typeOJ' 00109 !! end if 00110 ! 00111 end subroutine calc_omega_drot_secant
