Cocal: /home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine_mpatch/iter_irrot_BNS_CF

00001 subroutine iter_irrot_BNS_CF_mpt(iter_count, iseq, iter_eps)
00002   use phys_constant, only :  long, nmpt
00003   use grid_parameter
00004   use grid_parameter_interpo, only : nrg_itp, ntg_itp, npg_itp
00005   use def_metric, only : psi, alph, alps, bvxd, bvyd, bvzd
00006   use coordinate_grav_r
00007   use grid_points_binary_excision, only : irg_exin, irg_exout
00008 !
00009   use def_matter
00010   use def_velocity_potential
00011   use def_matter_parameter
00012   use interface_interpo_fl2gr
00013   use interface_sourceterm_HaC_peos_irrot
00014   use interface_sourceterm_HaC_CF
00015   use interface_sourceterm_trG_peos_irrot
00016   use interface_sourceterm_trG_CF
00017   use interface_sourceterm_MoC_CF_with_divshift_peos_irrot
00018   use interface_sourceterm_MoC_CF_with_divshift_r3rd
00019 !  use interface_violation_midpoint_MoC_CF_peos_irrot
00020 !  use interface_violation_gridpoint_MoC_CF_peos_irrot
00021 !  use interface_violation_midpoint_HaC_CF_peos_irrot
00022 !  use interface_violation_gridpoint_HaC_CF_peos_irrot
00023   use interface_hydrostatic_eq_CF_peos_irrot
00024   use interface_calc_gradvep
00025   use interface_hydro_irbns_vep_CF_peos
00026   use interface_calc_surface
00027   use interface_update_matter
00028   use interface_update_parameter_BNS_irrot
00029   use interface_update_surface
00030   use interface_error_matter
00031   use interface_error_matter_type2
00032   use interface_error_velocity_potential
00033 !  use interface_IO_output_3D_general
00034 !  use interface_IO_output_1D_general
00035 !
00036   use def_binary_parameter, only : dis
00037   use copy_array_4dto3d_mpt
00038   use make_array_2d
00039   use make_array_3d
00040   use make_array_4d
00041   use interface_compute_alps2alph
00042   use interface_sourceterm_surface_int
00043   use interface_sourceterm_exsurf_binary_COCP
00044   use interface_sourceterm_outsurf_COCP_from_ARCP
00045   use interface_sourceterm_insurf_ARCP_from_COCP
00046   use interface_outer_boundary_d_BBH_CF
00047   use interface_outer_boundary_d_psi
00048   use interface_outer_boundary_d_alps
00049   use interface_outer_boundary_d_bvxd
00050   use interface_outer_boundary_d_bvyd
00051   use interface_outer_boundary_d_bvzd
00052 !  use interface_poisson_solver_binary_bhex_homosol
00053 !  use interface_poisson_solver_binary_star_homosol_plmex
00054   use interface_poisson_solver_binary_star_homosol
00055   use interface_poisson_solver_asymptotic_patch_homosol
00056   use interface_interpolation_fillup_binary_COCP
00057   use interface_update_grfield
00058   use interface_error_metric
00059   use interface_error_metric_type0
00060   use interface_error_metric_type2
00061   use interface_error_metric_type2_mpt
00062   implicit none
00063   real(long), pointer :: sou(:,:,:), pot(:,:,:), HaC_vio(:,:,:)
00064   real(long), pointer :: potx(:,:,:), poty(:,:,:), potz(:,:,:) 
00065   real(long), pointer :: pot_bak(:,:,:), work(:,:,:)
00066   real(long), pointer :: potx_bak(:,:,:), poty_bak(:,:,:), potz_bak(:,:,:) 
00067   real(long), pointer :: souvec(:,:,:,:), MoC_vio(:,:,:,:)
00068   real(long), pointer :: sou_exsurf(:,:), dsou_exsurf(:,:)
00069 !  real(long), pointer :: sou_bhsurf(:,:), dsou_bhsurf(:,:)
00070   real(long), pointer :: sou_insurf(:,:), dsou_insurf(:,:)
00071   real(long), pointer :: sou_outsurf(:,:), dsou_outsurf(:,:)
00072   real(long), pointer :: fnc_itp(:,:,:)
00073   real(long), pointer :: potf(:,:,:), emd_bak(:,:,:), vep_bak(:,:,:)
00074   real(long), pointer :: potrs(:,:)
00075   real(long), pointer :: potut(:,:,:), potux(:,:,:), potuy(:,:,:), potuz(:,:,:)
00076 
00077   real(long) :: count, iter_eps
00078   real(long) :: error_psi  = 0.0d0, error_alph = 0.0d0, error_tmp  = 0.0d0, 
00079                error_bvxd = 0.0d0, error_bvyd = 0.0d0, error_bvzd = 0.0d0, &
00080                error_ome  = 2.0d0, error_emd  = 0.0d0, error_metr = 0.0d0, &
00081                error_vep  = 0.0d0, error_emd_last, error_psi_last, error_alph_last, &
00082                error_vep_last
00083   real(long) :: pari, ome_bh_new, work_shift, error_metric_mpt_mid(10,10)
00084   integer    :: iter_count, iter_extra = 0, istep_niq = -1,     
00085                flag_tmp  = 0, flag_param = 99, flag_emd  = 0,  &
00086                flag_psi  = 0, flag_alph  = 0,  flag_vep  = 0,  &
00087                flag_bvxd = 0, flag_bvyd  = 0,  flag_bvzd = 0
00088   integer    :: irg, itg, ipg, ii, iseq, jj, sumjj, je
00089   integer    :: impt, impt_ex, impt_bin
00090   character(len=2) :: chgreen, chpa, chpB, cocp
00091   character(len=4) :: chbe
00092   integer    :: flag = 0, hydro_iter = 4, ishift, shift_iter=1, 
00093                flag_mid,   flag_mpt_mid(10)
00094   integer    :: irf, itf, ipf, ihy, npg_l, npg_r, icount
00095   character(30) :: char1, char2, char3, char4, char5
00096   real(long) :: start, finish
00097   character(len=1) :: char_1
00098 !
00099   call set_allocate_size_mpt
00100   call alloc_array4d(souvec ,0,nrg,0,ntg,0,npg,1,3)
00101   call alloc_array3d(sou,0,nrg,0,ntg,0,npg)
00102   call alloc_array3d(pot,0,nrg,0,ntg,0,npg)
00103   call alloc_array3d(potx,0,nrg,0,ntg,0,npg)
00104   call alloc_array3d(poty,0,nrg,0,ntg,0,npg)
00105   call alloc_array3d(potz,0,nrg,0,ntg,0,npg)
00106   call alloc_array3d(pot_bak,0,nrg,0,ntg,0,npg)
00107   call alloc_array3d(potx_bak,0,nrg,0,ntg,0,npg)
00108   call alloc_array3d(poty_bak,0,nrg,0,ntg,0,npg)
00109   call alloc_array3d(potz_bak,0,nrg,0,ntg,0,npg)
00110   call alloc_array3d(work,0,nrg,0,ntg,0,npg)
00111   call alloc_array3d(fnc_itp,0,nrg,0,ntg,0,npg)
00112   call alloc_array2d(sou_exsurf,0,ntg,0,npg)
00113   call alloc_array2d(dsou_exsurf,0,ntg,0,npg)
00114   call alloc_array2d(sou_insurf,0,ntg,0,npg)
00115   call alloc_array2d(dsou_insurf,0,ntg,0,npg)
00116   call alloc_array2d(sou_outsurf,0,ntg,0,npg)
00117   call alloc_array2d(dsou_outsurf,0,ntg,0,npg)
00118   call alloc_array3d( potf,0,nrf,0,ntf,0,npf)
00119   call alloc_array3d(potut,0,nrf,0,ntf,0,npf)
00120   call alloc_array3d(potux,0,nrf,0,ntf,0,npf)
00121   call alloc_array3d(potuy,0,nrf,0,ntf,0,npf)
00122   call alloc_array3d(potuz,0,nrf,0,ntf,0,npf)
00123   call alloc_array3d(emd_bak,0,nrf,0,ntf,0,npf)
00124   call alloc_array3d(vep_bak,0,nrf,0,ntf,0,npf)
00125   call alloc_array2d(potrs,0,ntf,0,npf)
00126 !
00127   iter_count = 0
00128   error_emd_last =1.0d0
00129   error_psi_last =1.0d0
00130   error_alph_last=1.0d0
00131   flag_mid=0
00132 !
00133   do
00134     iter_count = iter_count + 1
00135 !    iter_extra = iter_extra + 1
00136     count = dble(iter_count)
00137     flag_mpt_mid(1:10) =0    !  one flag for every patch
00138     write(6,'(a10,i4,a108)') "Iteration:",iter_count, &
00139     &  " ---psi------------alph------------bvxd------------bvyd------------bvzd------------emd------------vep-------"
00140 !
00141     do impt = 1, nmpt
00142       write(6,'(a8,i2,a1)', ADVANCE = "NO")  " Patch #", impt, ""
00143       write(6,*)  " "
00144       call copy_grid_parameter_from_mpt(impt)
00145       conv_gra = dmin1(conv0_gra,conv_ini*count)
00146       conv_den = dmin1(conv0_den,conv_ini*count)
00147       conv_vep = dmin1(conv0_vep,conv_ini*count)
00148       call copy_grid_parameter_to_mpt(impt)
00149 ! ---
00150       call copy_from_mpatch_all_BNS_CF(impt)
00151       call copy_def_metric_and_matter_from_mpt(impt)
00152       call copy_def_matter_irrot_from_mpt(impt)
00153 ! ---
00154 !      write(6,'(a20,i5,i5,1p,3e20.12)') "================", iter_count, impt, ome, ber, radi
00155 
00156 !!!!      call excurve
00157 !      emdg = 0.0d0
00158       if (impt==1 .or. impt==2) then
00159         cocp = 'ns'
00160         call interpo_fl2gr(emd, emdg)
00161         call interpo_fl2gr(vepxf, vepxg)
00162         call interpo_fl2gr(vepyf, vepyg)
00163         call interpo_fl2gr(vepzf, vepzg)
00164 !       ********************************   CHECK   ***************************************************
00165         if (mod(iter_count,5)==0.and.impt==10) then
00166           write(char1, '(i5)') iter_count
00167           char2 = adjustl(char1)
00168           write(char4, '(i5)') impt
00169           char5 = adjustl(char4)
00170 
00171           char3 = 'vep_iteration' // trim(char2) // '_mpt' // trim(char5) // '.txt'
00172           open(12,file=char3,status='unknown')
00173           itg=12;  ipg=6
00174           write(12,'(a1,a2,4a16,a30)') '#', 'ir', 'rg*rs', 'vepxf', 'vepyf', 'vepzf', 'plot using every :::0::0'
00175           do irg=0,nrf
00176             write(12,'(i3,1p,4e16.6)')  irg, rg(irg)*rs(itg,ipg), vepxf(irg,itg,ipg), vepyf(irg,itg,ipg), vepzf(irg,itg,ipg)
00177           end do
00178           write(12,*) '#------------------------------------------------------------------------------------------------------'
00179           write(12,*) ""
00180           write(12,'(a1,a2,4a16,a30)') '#', 'ir', 'rg', 'vepxg', 'vepyg', 'vepzg',  'plot using every :::1::1'
00181           do irg=0,nrg
00182             write(12,'(i3,1p,4e16.6)')  irg, rg(irg), vepxg(irg,itg,ipg), vepyg(irg,itg,ipg), vepzg(irg,itg,ipg)
00183           end do
00184           write(12,*) ""
00185           close(12)
00186         endif
00187 !       *********************************************************************************************
00188 
00189         call calc_vector_x_matter(2)
00190         call calc_vector_phi_matter(2)
00191         call excurve_CF(cocp)             !   3rd order midpoint from ir0=-2,...
00192         call excurve_CF_gridpoint         !   4th order from ir0=-2,...
00193         call copy_def_metric_and_matter_to_mpt(impt)
00194         call copy_def_matter_irrot_to_mpt(impt)
00195         if (iter_count==1 .and. impt==1) then
00196           open(2,file='iter_data.txt',status='unknown')
00197         end if
00198         write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi 
00199       else
00200         cocp = 'bh'
00201         call excurve_CF(cocp)             !   3rd order midpoint from ir0=0,...
00202         call excurve_CF_gridpoint_bhex    !   4th order from ir0=0,...
00203         call copy_def_metric_and_matter_to_mpt(impt)
00204         call copy_def_matter_irrot_to_mpt(impt)
00205         write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi 
00206         write(2,*)  '-----------------------------------------------------------------------'
00207       end if
00208 !
00209       if (impt.eq.nmpt) call calc_mass_asympto('ns')
00210       call copy_to_mpatch_all_BNS_CF(impt)
00211 !
00212 ! --- psi
00213       if(error_psi_last.le.1.0d-13) then
00214         write(6,'(a16)', ADVANCE = "NO")    ' ## SKIP psi  ##'
00215         write(6,*) " "
00216       else
00217       sou(0:nrg,0:ntg,0:npg)  = 0.0d0
00218       sou_exsurf(0:ntg,0:npg) = 0.0d0 ; dsou_exsurf(0:ntg,0:npg) = 0.0d0
00219       sou_outsurf(0:ntg,0:npg)= 0.0d0 ; dsou_outsurf(0:ntg,0:npg)= 0.0d0
00220       if (impt==1 .or. impt==2) then
00221         call sourceterm_HaC_peos_irrot(sou)
00222       else
00223         call sourceterm_HaC_CF(sou)
00224       end if
00225 
00226       if (impt==1 .or. impt==2) then
00227         call sourceterm_exsurf_binary_COCP(impt,'psi ', 'ev', sou_exsurf, dsou_exsurf)
00228         call sourceterm_outsurf_COCP_from_ARCP(impt,'psi ', 'ev', sou_outsurf, dsou_outsurf)
00229 
00230         call cpu_time(start)
00231         call poisson_solver_binary_star_homosol('sd', sou, sou_exsurf, dsou_exsurf, &
00232         &                                                  sou_outsurf,dsou_outsurf,pot)
00233 
00234 !        call poisson_solver_binary_star_homosol_plmex(impt,'sd', sou, sou_exsurf, dsou_exsurf, &
00235 !        &                                                             sou_outsurf,dsou_outsurf,pot)
00236         call cpu_time(finish)
00237         print '("Time = ",f6.3," seconds.")',finish-start
00238 
00239       else if (impt.eq.nmpt) then
00240         call sourceterm_insurf_ARCP_from_COCP(impt,'psi ', 'ev', sou_insurf, dsou_insurf)
00241         call outer_boundary_d_psi(sou_outsurf)
00242 !        call outer_boundary_d_BBH_CF(sou_outsurf,'psi ')
00243         call poisson_solver_asymptotic_patch_homosol('dd',sou, sou_insurf, dsou_insurf, &
00244         &                                                      sou_outsurf,dsou_outsurf,pot)
00245       end if
00246 
00247       pot_bak(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)
00248       call update_grfield(pot,psi,conv_gra)
00249 
00250       if ((impt==1 .or. impt==2) .and. (iter_count>5 .or. iseq>1 .or. indata_type.eq.'3D')) &
00251         &          call update_parameter_BNS_irrot(conv_gra)
00252 
00253       if (impt.eq.nmpt)       psi(nrg,0:ntg,0:npg) = 1.0d0 
00254 !      call reset_outer_boundary_BBH_CF('psi ')
00255       if (impt.ne.nmpt) call interpolation_fillup_binary_COCP(impt,'psi ','ev',psi)
00256 !!!      call error_metric_type0(psi,pot_bak,error_tmp,flag_tmp,cocp)
00257       call error_metric_type2_mpt(psi,pot_bak,error_tmp,flag_tmp,cocp,impt)
00258       flag_psi  =  max ( flag_psi, flag_tmp)
00259       error_psi = dmax1(error_psi,error_tmp)
00260 !!!      write(6,'(1p,e14.6,a2)', ADVANCE = "NO")  error_psi, "  "
00261       error_psi_last = error_psi
00262 !      call printout_error_metric_v1('psi ', iter_count,error_psi)
00263       call copy_to_mpatch_all_BNS_CF(impt)
00264       call copy_def_metric_and_matter_to_mpt(impt)
00265       call copy_def_matter_irrot_to_mpt(impt)
00266       end if
00267 !
00268 ! --- alps
00269       if(error_alph_last.le.1.0d-13) then
00270         write(6,'(a16)', ADVANCE = "NO")    ' ## SKIP alph ##'
00271         write(6,*) " "
00272       else
00273       sou(0:nrg,0:ntg,0:npg)  = 0.0d0
00274       sou_exsurf(0:ntg,0:npg) = 0.0d0 ; dsou_exsurf(0:ntg,0:npg) = 0.0d0
00275       sou_outsurf(0:ntg,0:npg)= 0.0d0 ; dsou_outsurf(0:ntg,0:npg)= 0.0d0
00276       alps(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)*psi(0:nrg,0:ntg,0:npg)
00277       if (impt==1 .or. impt==2) then
00278         call sourceterm_trG_peos_irrot(sou)
00279       else
00280         call sourceterm_trG_CF(sou)
00281       end if
00282       if (impt==1 .or. impt==2) then
00283         call sourceterm_exsurf_binary_COCP(impt,'alps','ev',sou_exsurf, dsou_exsurf)
00284         call sourceterm_outsurf_COCP_from_ARCP(impt,'alps','ev',sou_outsurf, dsou_outsurf)
00285 
00286         call poisson_solver_binary_star_homosol('sd',sou, sou_exsurf, dsou_exsurf, &
00287         &                                                 sou_outsurf,dsou_outsurf,pot)
00288 
00289 !        call poisson_solver_binary_star_homosol_plmex(impt,'sd',sou, sou_exsurf, dsou_exsurf, &
00290 !        &                                                            sou_outsurf,dsou_outsurf,pot)
00291       else if (impt.eq.nmpt) then
00292         call sourceterm_insurf_ARCP_from_COCP(impt,'alps','ev',sou_insurf, dsou_insurf)
00293         call outer_boundary_d_alps(sou_outsurf)
00294 !        call outer_boundary_d_BBH_CF(sou_outsurf,'alps')
00295         call poisson_solver_asymptotic_patch_homosol('dd',sou, sou_insurf,dsou_insurf, & 
00296         &                                                      sou_outsurf,dsou_outsurf,pot)
00297       end if
00298 
00299       call compute_alps2alph(pot,psi)
00300       pot_bak(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)
00301       call update_grfield(pot,alph,conv_gra)
00302 
00303       if ((impt==1 .or. impt==2) .and. (iter_count>5 .or. iseq>1 .or. indata_type.eq.'3D')) &
00304         &          call update_parameter_BNS_irrot(conv_gra)
00305 
00306       if (impt.eq.nmpt)       alph(nrg,0:ntg,0:npg) = 1.0d0
00307 !      call reset_outer_boundary_BBH_CF('alph')
00308       if (impt.ne.nmpt) call interpolation_fillup_binary_COCP(impt,'alph','ev',alph)
00309 !!!      call error_metric_type0(alph,pot_bak,error_tmp,flag_tmp,cocp)
00310       call error_metric_type2_mpt(alph,pot_bak,error_tmp,flag_tmp,cocp,impt)
00311       flag_alph  =  max ( flag_alph, flag_tmp)
00312       error_alph = dmax1(error_alph,error_tmp)
00313 !!!      write(6,'(1p,e14.6,a2)', ADVANCE = "NO")  error_alph, "  "
00314       error_alph_last = error_alph
00315 !      call printout_error_metric_v1('alph', iter_count,error_alph)
00316       alps(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)*psi(0:nrg,0:ntg,0:npg)
00317       call copy_to_mpatch_all_BNS_CF(impt)
00318       call copy_def_metric_and_matter_to_mpt(impt)
00319       call copy_def_matter_irrot_to_mpt(impt)
00320       end if
00321 ! 
00322 ! --- shift
00323       do ishift=1,shift_iter
00324       souvec(0:nrg,0:ntg,0:npg,1:3)  = 0.0d0
00325       if (impt==1 .or. impt==2) then
00326         call sourceterm_MoC_CF_with_divshift_peos_irrot(souvec)
00327       else
00328         call sourceterm_MoC_CF_with_divshift_r3rd(souvec)
00329       end if
00330 
00331       do ii = 1, 3
00332         sou_exsurf(0:ntg,0:npg) = 0.0d0 ; dsou_exsurf(0:ntg,0:npg) = 0.0d0
00333         sou_outsurf(0:ntg,0:npg)= 0.0d0 ; dsou_outsurf(0:ntg,0:npg)= 0.0d0
00334         sou(0:nrg,0:ntg,0:npg) = souvec(0:nrg,0:ntg,0:npg,ii)
00335         if (ii.eq.1) work(0:nrg,0:ntg,0:npg) = bvxd(0:nrg,0:ntg,0:npg)
00336         if (ii.eq.2) work(0:nrg,0:ntg,0:npg) = bvyd(0:nrg,0:ntg,0:npg)
00337         if (ii.eq.3) work(0:nrg,0:ntg,0:npg) = bvzd(0:nrg,0:ntg,0:npg)
00338         if (ii.eq.1) then; chbe = 'bvxd'; chpa = 'od'; end if
00339         if (ii.eq.2) then; chbe = 'bvyd'; chpa = 'od'; end if
00340         if (ii.eq.3) then; chbe = 'bvzd'; chpa = 'ev'; end if
00341         if (impt==1 .or. impt==2) then
00342           chpB = 'ev'
00343           if (impt.eq.2.and.ii.eq.1) chpB = 'od'
00344           if (impt.eq.2.and.ii.eq.2) chpB = 'od'
00345           call sourceterm_exsurf_binary_COCP(impt,chbe,chpa,sou_exsurf, dsou_exsurf)
00346           call sourceterm_outsurf_COCP_from_ARCP(impt,chbe,chpB,sou_outsurf, dsou_outsurf)
00347 
00348           call poisson_solver_binary_star_homosol('sd',sou, sou_exsurf,dsou_exsurf, &
00349           &                                                 sou_outsurf,dsou_outsurf,pot)
00350 
00351 !          call poisson_solver_binary_star_homosol_plmex(impt,'sd',sou,sou_exsurf,dsou_exsurf, &
00352 !          &                                                           sou_outsurf,dsou_outsurf,pot)
00353         else if (impt.eq.nmpt) then
00354           call sourceterm_insurf_ARCP_from_COCP(impt,chbe,chpa,sou_insurf, dsou_insurf)
00355           if (ii.eq.1) call outer_boundary_d_bvxd(sou_outsurf)
00356           if (ii.eq.2) call outer_boundary_d_bvyd(sou_outsurf)
00357           if (ii.eq.3) call outer_boundary_d_bvzd(sou_outsurf)
00358 !          if (ii.eq.1) call outer_boundary_d_BBH_CF(sou_outsurf,'bvxd')
00359 !          if (ii.eq.2) call outer_boundary_d_BBH_CF(sou_outsurf,'bvyd')
00360 !          if (ii.eq.3) call outer_boundary_d_BBH_CF(sou_outsurf,'bvzd')
00361           call poisson_solver_asymptotic_patch_homosol('dd',sou, sou_insurf,dsou_insurf, &
00362           &                                                      sou_outsurf,dsou_outsurf,pot)
00363         end if
00364         if (ii.eq.1) potx(0:nrg,0:ntg,0:npg) = pot(0:nrg,0:ntg,0:npg)
00365         if (ii.eq.2) poty(0:nrg,0:ntg,0:npg) = pot(0:nrg,0:ntg,0:npg)
00366         if (ii.eq.3) potz(0:nrg,0:ntg,0:npg) = pot(0:nrg,0:ntg,0:npg)
00367       end do
00368 
00369       potx_bak(0:nrg,0:ntg,0:npg) = bvxd(0:nrg,0:ntg,0:npg)
00370       poty_bak(0:nrg,0:ntg,0:npg) = bvyd(0:nrg,0:ntg,0:npg)
00371       potz_bak(0:nrg,0:ntg,0:npg) = bvzd(0:nrg,0:ntg,0:npg)
00372       call update_grfield(potx,bvxd,conv_gra)
00373       call update_grfield(poty,bvyd,conv_gra)
00374       call update_grfield(potz,bvzd,conv_gra)
00375 
00376       if ((impt==1 .or. impt==2) .and. (iter_count>5 .or. iseq>1 .or. indata_type.eq.'3D')) &
00377         &          call update_parameter_BNS_irrot(conv_gra)
00378 
00379       if (impt.ne.nmpt) then 
00380         call interpolation_fillup_binary_COCP(impt,'bvxd','od',bvxd)
00381         call interpolation_fillup_binary_COCP(impt,'bvyd','od',bvyd)
00382         call interpolation_fillup_binary_COCP(impt,'bvzd','ev',bvzd)
00383       end if
00384 !      call reset_metric_CF
00385       call error_metric_type2_mpt(bvxd,potx_bak,error_tmp,flag_tmp,cocp,impt)
00386       flag_bvxd  =  max ( flag_bvxd, flag_tmp)
00387       error_bvxd = dmax1(error_bvxd,error_tmp)
00388       call error_metric_type2_mpt(bvyd,poty_bak,error_tmp,flag_tmp,cocp,impt)
00389       flag_bvyd  =  max ( flag_bvyd, flag_tmp)
00390       error_bvyd = dmax1(error_bvyd,error_tmp)
00391       call error_metric_type2_mpt(bvzd,potz_bak,error_tmp,flag_tmp,cocp,impt)
00392       flag_bvzd  =  max ( flag_bvzd, flag_tmp)
00393       error_bvzd = dmax1(error_bvzd,error_tmp)
00394 !!!      if(ishift==shift_iter)  write(6,'(1p,e14.6,a2)', ADVANCE = "NO")  error_bvxd, "  "
00395 !!!      if(ishift==shift_iter)  write(6,'(1p,e14.6,a2)', ADVANCE = "NO")  error_bvyd, "  "
00396 !!!      if(ishift==shift_iter)  write(6,'(1p,e14.6,a2)', advance = "no")  error_bvzd, "  "
00397       end do   ! ==========>> ishift 
00398 
00399 !      call printout_error_metric_v1('bvxd', iter_count,error_bvxd)
00400 !      call printout_error_metric_v1('bvyd', iter_count,error_bvyd)
00401 !      call printout_error_metric_v1('bvzd', iter_count,error_bvzd)
00402       call copy_to_mpatch_all_BNS_CF(impt)
00403       call copy_def_metric_and_matter_to_mpt(impt)
00404       call copy_def_matter_irrot_to_mpt(impt)
00405 !
00406       error_metr = dmax1(error_psi,error_alph,error_bvxd,error_bvyd,error_bvzd)
00407 
00408       if(iter_count<5 .and. iseq<=1 .and. indata_type.eq.'1D') then
00409         write(6,'(a14)') '#SKIP HYDRO a#' 
00410         error_emd_last = 1.0d0 
00411         cycle
00412       end if
00413 ! ---
00414 ! -- Hydro equations.
00415 !      if (error_emd_last.le.eps*0.0001d0 .and. error_metr.gt.eps*2.0d0) then
00416       if (error_emd_last.le.eps*0.0001d0 .and. error_vep_last.le.eps*0.0001d0 ) then
00417         write(6,'(a14)')  '#SKIP HYDRO b#'
00418       else
00419         if (impt==1 .or. impt==2) then
00420           do ihy = 1, hydro_iter
00421             call interpo_gr2fl_metric_CF
00422             call hydrostatic_eq_CF_peos_irrot(potf,potut,potux,potuy,potuz)
00423             call calc_surface(potrs, potf)
00424             emd_bak = emd 
00425             call update_matter( potf,emd,conv_den)
00426 !            call update_matter(potut,utf,conv_den)
00427 !            call update_matter(potux,vepxf,conv_den)
00428 !            call update_matter(potuy,vepyf,conv_den)
00429 !            call update_matter(potuz,vepzf,conv_den)
00430             call update_surface(potrs,rs,conv_den)
00431 
00432             call reset_fluid_BNS
00433             call update_parameter_BNS_irrot(conv_den)
00434             call calc_vector_x_matter(2)
00435             call calc_vector_phi_matter(2)
00436 !
00437             call hydro_irbns_vep_CF_peos(potf,iter_count,impt,10)
00438             vep_bak = vep
00439             call update_matter(potf,vep,conv_vep)
00440             call reset_fluid_vep
00441             call calc_gradvep(vep,vepxf,vepyf,vepzf)
00442             call reset_fluid_gradvep
00443             call update_parameter_BNS_irrot(conv_vep)
00444 
00445             if (ihy.eq.hydro_iter) then
00446 !            if (ihy.eq.1) then
00447               call error_matter_type2(emd,emd_bak,error_tmp,flag_tmp)
00448               flag_emd  =  max ( flag_emd, flag_tmp)
00449               error_emd = dmax1(error_emd,error_tmp)
00450               error_emd_last = error_emd
00451 !!!              write(6,'(1p,e14.6,a2)', ADVANCE = "NO")  error_emd, "  "
00452 !              call printout_error_matter("emd ", iter_count,error_emd)
00453               call error_velocity_potential(vep,vep_bak,error_tmp,flag_tmp)
00454               flag_vep  =  max ( flag_vep, flag_tmp)
00455               error_vep = dmax1(error_vep,error_tmp)
00456               error_vep_last = error_vep
00457 !!!              write(6,'(1p,e14.6,a2)')  error_vep, "  "        
00458             end if
00459           end do   ! end hydro
00460           call copy_to_mpatch_all_BNS_CF(impt)
00461           call copy_def_metric_and_matter_to_mpt(impt)
00462           call copy_def_matter_irrot_to_mpt(impt)
00463         else
00464           write(6,*)   ' ##   NO HYDRO ##'
00465         end if
00466       end if
00467 !      error_emd_last = error_emd
00468 !
00469       if(impt>5) then
00470         if (mod(iter_count,1)==0.or.iter_count==1) then
00471           write(char1, '(i5)') iter_count
00472           char2 = adjustl(char1)
00473           write(char4, '(i5)') impt
00474           char5 = adjustl(char4)
00475           char3 = 'iteration' // trim(char2) // '_mpt' // trim(char5) // '.txt'
00476           call printout_axis_irrot_BNS_mpt(impt,char3)
00477         endif
00478       end if  !  ====> iteration.txt
00479 
00480 !     print once the solution when the error in shift for all patches is <0.1 
00481       if (flag_mid==0) then
00482         if(error_metr < iter_eps ) then
00483           flag_mpt_mid(impt) = 1
00484           error_metric_mpt_mid(impt,1) = error_psi
00485           error_metric_mpt_mid(impt,2) = error_alph
00486           error_metric_mpt_mid(impt,3) = error_bvxd
00487           error_metric_mpt_mid(impt,4) = error_bvyd
00488           error_metric_mpt_mid(impt,5) = error_bvzd
00489         end if
00490         sumjj=0
00491         do jj=1,nmpt 
00492           sumjj = sumjj + flag_mpt_mid(jj)
00493         end do
00494         if (sumjj==nmpt) then
00495           open(3,file='mid_error_data.txt',status='unknown')
00496           write(3,'(a11,i5)') 'Iteration =', iter_count
00497           write(3,'(a8,1p,5e17.6)') 'COCP #=1', (error_metric_mpt_mid(1,je), je=1,5)
00498           write(3,'(a8,1p,5e17.6)') 'COCP #=2', (error_metric_mpt_mid(2,je), je=1,5)
00499           write(3,'(a8,1p,5e17.6)') 'ARCP #  ', (error_metric_mpt_mid(3,je), je=1,5)
00500           close(3)
00501 !          do jj=1,nmpt
00502 !            call copy_from_mpatch_all_BNS_CF(jj)
00503 !            call copy_def_metric_and_matter_from_mpt(jj)
00504 !            call copy_def_matter_irrot_from_mpt(jj)
00505 !            call IO_output_solution_3D_CF_irrot_NS_mpt(jj,'.mid')
00506 !            if (jj==1 .or. jj==2)  call printout_NS_shape_mpt(jj)
00507 !          end do
00508           flag_mid = 1
00509         end if 
00510       end if
00511     end do    !  ====================>> end all patches
00512 !
00513 ! -- convergence of iteration is checked.
00514 !
00515 !    if (iter_count==6) then
00516 !      do impt=1,nmpt
00517 !        call copy_from_mpatch_all_BNS_CF(impt)
00518 !        write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi
00519 !      end do
00520 !      write(2,*) '-----------------------------------------------------------------------'
00521 !      close(2)
00522 !      exit
00523 !    end if
00524 
00525     if (flag_mid==1) then
00526       do impt=1,nmpt
00527         call copy_from_mpatch_all_BNS_CF(impt)
00528 !        call copy_def_metric_and_matter_from_mpt(impt)
00529 !        call copy_def_matter_irrot_from_mpt(impt)
00530         write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi
00531 
00532 !        MoC_vio(0:nrg,0:ntg,0:npg,1:3) = 0.0d0
00533 !        call violation_gridpoint_MoC_CF_peos_irrot(MoC_vio,'ns')
00534 !        pot = 0.0d0
00535 !        pot(0:nrg,0:ntg,0:npg) = MoC_vio(0:nrg,0:ntg,0:npg,2)   
00536 !        char1 = 'MoC_vio_by_3D_mpt1.txt'
00537 !        call IO_output_3D_general(char1,'g',pot)
00538 !        char1 = 'MoC_vio_by_1D_mpt1.txt'
00539 !        call IO_output_1D_general(char1,'g',pot,-1,ntg/2,0)
00540 
00541 !        HaC_vio = 0.0d0
00542 !        call violation_gridpoint_HaC_CF_peos_irrot(HaC_vio,'ns')
00543 !        char1 = 'HaC_vio_3D_mpt1.txt'
00544 !        call IO_output_3D_general(char1,'g',HaC_vio)
00545 !        char1 = 'HaC_vio_1D_mpt1.txt'
00546 !        call IO_output_1D_general(char1,'g',HaC_vio,-1,ntg/2,0)
00547       end do
00548       write(2,*) '-----------------------------------------------------------------------'
00549       close(2)
00550       exit
00551     end if
00552 
00553     if ( (flag_vep==0).and.(flag_emd==0).and.(flag_psi==0).and.(flag_alph==0).and. &
00554     &   (flag_bvxd==0).and.(flag_bvyd==0).and.(flag_bvzd==0))      then
00555       do impt=1, nmpt
00556         call copy_from_mpatch_all_BNS_CF(impt)
00557         write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi
00558       end do
00559       write(2,*) '-----------------------------------------------------------------------'
00560       close(2)
00561       exit
00562     endif
00563 
00564     if (iter_count >= iter_max) then
00565       write(6,*) "*************** iter_count==iter_max ....exiting *******************"
00566       stop
00567 !      exit
00568     end if
00569 !    error_emd_last = error_emd
00570 !
00571     flag_psi  = 0
00572     flag_alph = 0
00573     flag_bvxd = 0
00574     flag_bvyd = 0
00575     flag_bvzd = 0
00576     flag_emd  = 0
00577     flag_vep  = 0
00578     error_psi  = 0.0d0
00579     error_alph = 0.0d0
00580     error_bvxd = 0.0d0
00581     error_bvyd = 0.0d0
00582     error_bvzd = 0.0d0
00583     error_emd  = 0.0d0
00584     error_vep  = 0.0d0
00585   end do
00586 !
00587   deallocate(souvec)
00588   deallocate(sou)
00589   deallocate(pot)
00590   deallocate(potx)
00591   deallocate(poty)
00592   deallocate(potz)
00593   deallocate(pot_bak)
00594   deallocate(potx_bak)
00595   deallocate(poty_bak)
00596   deallocate(potz_bak)
00597   deallocate(work)
00598   deallocate(fnc_itp)
00599   deallocate(sou_exsurf)
00600   deallocate(dsou_exsurf)
00601   deallocate(sou_insurf)
00602   deallocate(dsou_insurf)
00603   deallocate(sou_outsurf)
00604   deallocate(dsou_outsurf)
00605   deallocate(potf)
00606   deallocate(potut)
00607   deallocate(potux)
00608   deallocate(potuy)
00609   deallocate(potuz)
00610   deallocate(emd_bak)
00611   deallocate(vep_bak)
00612   deallocate(potrs)
00613   
00614 end subroutine iter_irrot_BNS_CF_mpt
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine_mpatch/iter_irrot_BNS_CF_mpt.f90
