00001 subroutine adjust_omega_BBH(flag_omega,error_ome)
00002 use grid_parameter, only : mass_eps, eps
00003 use def_quantities, only : admmass_asymp, komarmass_asymp
00004 use def_bh_parameter, only : ome_bh, spin_bh
00005 implicit none
00006 integer :: flag_omega
00007 integer, save :: count_adj, ii
00008 real(8), save :: fnc_val, fnc_val_bak, dfdval,
00009 ome_bh_bak, delta_ome_bh, delta_ome_bh_bak
00010 real(8), save :: facp(5) = (/ 1.0d-1, 2.0d-1, 4.0d-1, 7.0d-1, 1.0d-0 /)
00011
00012 real(8) :: error_ome, fac0 = 1.0d0
00013
00014
00015
00016
00017
00018 if (flag_omega.eq.99) then
00019 if (admmass_asymp.le.komarmass_asymp) delta_ome_bh = -0.05d0*ome_bh
00020 if (admmass_asymp.ge.komarmass_asymp) delta_ome_bh = 0.05d0*ome_bh
00021 ome_bh = ome_bh + delta_ome_bh
00022 fnc_val = (admmass_asymp - komarmass_asymp)/admmass_asymp
00023 flag_omega = 1
00024 count_adj = 0
00025 return
00026 end if
00027 fnc_val_bak = fnc_val
00028 delta_ome_bh_bak = delta_ome_bh
00029 fnc_val = (admmass_asymp - komarmass_asymp)/admmass_asymp
00030 dfdval = (fnc_val - fnc_val_bak)/delta_ome_bh
00031 delta_ome_bh = - fnc_val/dfdval
00032
00033 error_ome = dabs(delta_ome_bh/ome_bh)
00034
00035 if (error_ome.le.eps) then
00036 fnc_val = fnc_val_bak
00037 delta_ome_bh = delta_ome_bh_bak
00038 flag_omega = 0
00039 return
00040 end if
00041
00042 count_adj = count_adj + 1
00043 ii = min0(5,count_adj)
00044 fac0 = facp(ii)
00045 delta_ome_bh = fac0*delta_ome_bh
00046 ome_bh_bak = ome_bh
00047 ome_bh = ome_bh + delta_ome_bh
00048
00049 write(6,'(1a24,2(7x,i5))') ' -- ome iteration # -- ',count_adj
00050 write(6,'(1a24,1p,2e12.4)')' -- ome OLD -> NEW -- ',ome_bh_bak, ome_bh
00051 write(6,*)'-- ADM and Komar --', admmass_asymp, komarmass_asymp
00052 write(6,*)'-- 1 - MK/MADM --', fnc_val_bak, fnc_val
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 return
00064 end subroutine adjust_omega_BBH