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

00001 subroutine iter_lecc_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_lecc_HaC_peos_irrot
00014   use interface_sourceterm_HaC_CF
00015   use interface_sourceterm_lecc_trG_peos_irrot
00016   use interface_sourceterm_trG_CF
00017   use interface_sourceterm_lecc_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_lecc_irrot
00024   use interface_calc_gradvep
00025   use interface_hydro_irbns_vep_CF_peos_lecc
00026   use interface_calc_surface
00027   use interface_update_matter
00028   use interface_update_parameter_BNS
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 
00165         call interpo_gr2fl_metric_CF
00166         call rotation_law  ! calculate omef and omeg
00167 
00168 !       ********************************   CHECK   ***************************************************
00169         if (mod(iter_count,5)==0.and.impt==10) then
00170           write(char1, '(i5)') iter_count
00171           char2 = adjustl(char1)
00172           write(char4, '(i5)') impt
00173           char5 = adjustl(char4)
00174 
00175           char3 = 'vep_iteration' // trim(char2) // '_mpt' // trim(char5) // '.txt'
00176           open(12,file=char3,status='unknown')
00177           itg=12;  ipg=6
00178           write(12,'(a1,a2,4a16,a30)') '#', 'ir', 'rg*rs', 'vepxf', 'vepyf', 'vepzf', 'plot using every :::0::0'
00179           do irg=0,nrf
00180             write(12,'(i3,1p,4e16.6)')  irg, rg(irg)*rs(itg,ipg), vepxf(irg,itg,ipg), vepyf(irg,itg,ipg), vepzf(irg,itg,ipg)
00181           end do
00182           write(12,*) '#------------------------------------------------------------------------------------------------------'
00183           write(12,*) ""
00184           write(12,'(a1,a2,4a16,a30)') '#', 'ir', 'rg', 'vepxg', 'vepyg', 'vepzg',  'plot using every :::1::1'
00185           do irg=0,nrg
00186             write(12,'(i3,1p,4e16.6)')  irg, rg(irg), vepxg(irg,itg,ipg), vepyg(irg,itg,ipg), vepzg(irg,itg,ipg)
00187           end do
00188           write(12,*) ""
00189           close(12)
00190         endif
00191 !       *********************************************************************************************
00192 
00193         call calc_vector_x_matter(2)
00194         call calc_vector_phi_matter(2)
00195         call excurve_CF(cocp)             !   3rd order midpoint from ir0=-2,...
00196         call excurve_CF_gridpoint         !   4th order from ir0=-2,...
00197         call copy_def_metric_and_matter_to_mpt(impt)
00198         call copy_def_matter_irrot_to_mpt(impt)
00199         if (iter_count==1 .and. impt==1) then
00200           open(2,file='iter_data.txt',status='unknown')
00201         end if
00202         write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi 
00203         write(6,'(a16,1p,2e20.12)') "(delome,delvel)=", delome, delvel
00204         write(6,'(a16,1p,2e20.12)') "(   ome,  radi)=", ome, radi
00205       else
00206         cocp = 'bh'
00207         call excurve_CF(cocp)             !   3rd order midpoint from ir0=0,...
00208         call excurve_CF_gridpoint_bhex    !   4th order from ir0=0,...
00209         call copy_def_metric_and_matter_to_mpt(impt)
00210         call copy_def_matter_irrot_to_mpt(impt)
00211         write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi 
00212         write(2,*)  '-----------------------------------------------------------------------'
00213       end if
00214 !
00215       if (impt.eq.nmpt) call calc_mass_asympto('ns')
00216       call copy_to_mpatch_all_BNS_CF(impt)
00217 !
00218 ! --- psi
00219       if(error_psi_last.le.1.0d-13) then
00220         write(6,'(a16)', ADVANCE = "NO")    ' ## SKIP psi  ##'
00221         write(6,*) " "
00222       else 
00223       write(6,'(a21,1p,4e20.12)') "(velx,ome,radi,emdc)=", velx, ome, radi, emdc
00224 
00225       sou(0:nrg,0:ntg,0:npg)  = 0.0d0
00226       sou_exsurf(0:ntg,0:npg) = 0.0d0 ; dsou_exsurf(0:ntg,0:npg) = 0.0d0
00227       sou_outsurf(0:ntg,0:npg)= 0.0d0 ; dsou_outsurf(0:ntg,0:npg)= 0.0d0
00228       if (impt==1 .or. impt==2) then
00229         call sourceterm_lecc_HaC_peos_irrot(sou)
00230       else
00231         call sourceterm_HaC_CF(sou)
00232       end if
00233 
00234       if (impt==1 .or. impt==2) then
00235         call sourceterm_exsurf_binary_COCP(impt,'psi ', 'ev', sou_exsurf, dsou_exsurf)
00236         call sourceterm_outsurf_COCP_from_ARCP(impt,'psi ', 'ev', sou_outsurf, dsou_outsurf)
00237 !        call poisson_solver_binary_star_homosol('sd', sou, sou_exsurf, dsou_exsurf, &
00238 !        &                                                  sou_outsurf,dsou_outsurf,pot)
00239 
00240         call cpu_time(start)
00241         call poisson_solver_binary_star_homosol_plmex(impt,'sd', sou, sou_exsurf, dsou_exsurf, &
00242         &                                                  sou_outsurf,dsou_outsurf,pot)
00243         call cpu_time(finish)
00244         print '("Time = ",f6.3," seconds.")',finish-start
00245       else if (impt.eq.nmpt) then
00246         call sourceterm_insurf_ARCP_from_COCP(impt,'psi ', 'ev', sou_insurf, dsou_insurf)
00247         call outer_boundary_d_psi(sou_outsurf)
00248 !        call outer_boundary_d_BBH_CF(sou_outsurf,'psi ')
00249         call poisson_solver_asymptotic_patch_homosol('dd',sou, sou_insurf, dsou_insurf, &
00250         &                                                      sou_outsurf,dsou_outsurf,pot)
00251       end if
00252 
00253       pot_bak(0:nrg,0:ntg,0:npg) = psi(0:nrg,0:ntg,0:npg)
00254       call update_grfield(pot,psi,conv_gra)
00255 
00256       if ((impt==1 .or. impt==2) .and. (iter_count>5 .or. iseq>1 .or. indata_type.eq.'3D')) &
00257         &          call update_parameter_BNS(conv_gra)
00258 
00259       if (impt.eq.nmpt)       psi(nrg,0:ntg,0:npg) = 1.0d0 
00260 !      call reset_outer_boundary_BBH_CF('psi ')
00261       if (impt.ne.nmpt) call interpolation_fillup_binary_COCP(impt,'psi ','ev',psi)
00262 !!!      call error_metric_type0(psi,pot_bak,error_tmp,flag_tmp,cocp)
00263       call error_metric_type2_mpt(psi,pot_bak,error_tmp,flag_tmp,cocp,impt)
00264       flag_psi  =  max ( flag_psi, flag_tmp)
00265       error_psi = dmax1(error_psi,error_tmp)
00266 !!!      write(6,'(1p,e14.6,a2)', ADVANCE = "NO")  error_psi, "  "
00267       error_psi_last = error_psi
00268 !      call printout_error_metric_v1('psi ', iter_count,error_psi)
00269       call copy_to_mpatch_all_BNS_CF(impt)
00270       call copy_def_metric_and_matter_to_mpt(impt)
00271       call copy_def_matter_irrot_to_mpt(impt)
00272       end if
00273 !
00274 ! --- alps
00275       if(error_alph_last.le.1.0d-13) then
00276         write(6,'(a16)', ADVANCE = "NO")    ' ## SKIP alph ##'
00277         write(6,*) " "
00278       else
00279       sou(0:nrg,0:ntg,0:npg)  = 0.0d0
00280       sou_exsurf(0:ntg,0:npg) = 0.0d0 ; dsou_exsurf(0:ntg,0:npg) = 0.0d0
00281       sou_outsurf(0:ntg,0:npg)= 0.0d0 ; dsou_outsurf(0:ntg,0:npg)= 0.0d0
00282       alps(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)*psi(0:nrg,0:ntg,0:npg)
00283       if (impt==1 .or. impt==2) then
00284         call sourceterm_lecc_trG_peos_irrot(sou)
00285       else
00286         call sourceterm_trG_CF(sou)
00287       end if
00288       if (impt==1 .or. impt==2) then
00289         call sourceterm_exsurf_binary_COCP(impt,'alps','ev',sou_exsurf, dsou_exsurf)
00290         call sourceterm_outsurf_COCP_from_ARCP(impt,'alps','ev',sou_outsurf, dsou_outsurf)
00291 !        call poisson_solver_binary_star_homosol('sd',sou, sou_exsurf, dsou_exsurf, &
00292 !        &                                                 sou_outsurf,dsou_outsurf,pot)
00293 
00294         call poisson_solver_binary_star_homosol_plmex(impt,'sd',sou, sou_exsurf, dsou_exsurf, &
00295         &                                                 sou_outsurf,dsou_outsurf,pot)
00296       else if (impt.eq.nmpt) then
00297         call sourceterm_insurf_ARCP_from_COCP(impt,'alps','ev',sou_insurf, dsou_insurf)
00298         call outer_boundary_d_alps(sou_outsurf)
00299 !        call outer_boundary_d_BBH_CF(sou_outsurf,'alps')
00300         call poisson_solver_asymptotic_patch_homosol('dd',sou, sou_insurf,dsou_insurf, & 
00301         &                                                      sou_outsurf,dsou_outsurf,pot)
00302       end if
00303 
00304       call compute_alps2alph(pot,psi)
00305       pot_bak(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)
00306       call update_grfield(pot,alph,conv_gra)
00307 
00308       if ((impt==1 .or. impt==2) .and. (iter_count>5 .or. iseq>1 .or. indata_type.eq.'3D')) &
00309         &          call update_parameter_BNS(conv_gra)
00310 
00311       if (impt.eq.nmpt)       alph(nrg,0:ntg,0:npg) = 1.0d0
00312 !      call reset_outer_boundary_BBH_CF('alph')
00313       if (impt.ne.nmpt) call interpolation_fillup_binary_COCP(impt,'alph','ev',alph)
00314 !!!      call error_metric_type0(alph,pot_bak,error_tmp,flag_tmp,cocp)
00315       call error_metric_type2_mpt(alph,pot_bak,error_tmp,flag_tmp,cocp,impt)
00316       flag_alph  =  max ( flag_alph, flag_tmp)
00317       error_alph = dmax1(error_alph,error_tmp)
00318 !!!      write(6,'(1p,e14.6,a2)', ADVANCE = "NO")  error_alph, "  "
00319       error_alph_last = error_alph
00320 !      call printout_error_metric_v1('alph', iter_count,error_alph)
00321       alps(0:nrg,0:ntg,0:npg) = alph(0:nrg,0:ntg,0:npg)*psi(0:nrg,0:ntg,0:npg)
00322       call copy_to_mpatch_all_BNS_CF(impt)
00323       call copy_def_metric_and_matter_to_mpt(impt)
00324       call copy_def_matter_irrot_to_mpt(impt)
00325       end if
00326 ! 
00327 ! --- shift
00328       do ishift=1,shift_iter
00329       souvec(0:nrg,0:ntg,0:npg,1:3)  = 0.0d0
00330       if (impt==1 .or. impt==2) then
00331         call sourceterm_lecc_MoC_CF_with_divshift_peos_irrot(souvec)
00332       else
00333         call sourceterm_MoC_CF_with_divshift_r3rd(souvec)
00334       end if
00335 
00336       do ii = 1, 3
00337         sou_exsurf(0:ntg,0:npg) = 0.0d0 ; dsou_exsurf(0:ntg,0:npg) = 0.0d0
00338         sou_outsurf(0:ntg,0:npg)= 0.0d0 ; dsou_outsurf(0:ntg,0:npg)= 0.0d0
00339         sou(0:nrg,0:ntg,0:npg) = souvec(0:nrg,0:ntg,0:npg,ii)
00340         if (ii.eq.1) work(0:nrg,0:ntg,0:npg) = bvxd(0:nrg,0:ntg,0:npg)
00341         if (ii.eq.2) work(0:nrg,0:ntg,0:npg) = bvyd(0:nrg,0:ntg,0:npg)
00342         if (ii.eq.3) work(0:nrg,0:ntg,0:npg) = bvzd(0:nrg,0:ntg,0:npg)
00343         if (ii.eq.1) then; chbe = 'bvxd'; chpa = 'od'; end if
00344         if (ii.eq.2) then; chbe = 'bvyd'; chpa = 'od'; end if
00345         if (ii.eq.3) then; chbe = 'bvzd'; chpa = 'ev'; end if
00346         if (impt==1 .or. impt==2) then
00347           chpB = 'ev'
00348           if (impt.eq.2.and.ii.eq.1) chpB = 'od'
00349           if (impt.eq.2.and.ii.eq.2) chpB = 'od'
00350           call sourceterm_exsurf_binary_COCP(impt,chbe,chpa,sou_exsurf, dsou_exsurf)
00351           call sourceterm_outsurf_COCP_from_ARCP(impt,chbe,chpB,sou_outsurf, dsou_outsurf)
00352 !          call poisson_solver_binary_star_homosol('sd',sou, sou_exsurf,dsou_exsurf, &
00353 !          &                                                 sou_outsurf,dsou_outsurf,pot)
00354 
00355           call poisson_solver_binary_star_homosol_plmex(impt,'sd',sou, sou_exsurf,dsou_exsurf, &
00356           &                                                 sou_outsurf,dsou_outsurf,pot)
00357         else if (impt.eq.nmpt) then
00358           call sourceterm_insurf_ARCP_from_COCP(impt,chbe,chpa,sou_insurf, dsou_insurf)
00359           if (ii.eq.1) call outer_boundary_d_bvxd(sou_outsurf)
00360           if (ii.eq.2) call outer_boundary_d_bvyd(sou_outsurf)
00361           if (ii.eq.3) call outer_boundary_d_bvzd(sou_outsurf)
00362 !          if (ii.eq.1) call outer_boundary_d_BBH_CF(sou_outsurf,'bvxd')
00363 !          if (ii.eq.2) call outer_boundary_d_BBH_CF(sou_outsurf,'bvyd')
00364 !          if (ii.eq.3) call outer_boundary_d_BBH_CF(sou_outsurf,'bvzd')
00365           call poisson_solver_asymptotic_patch_homosol('dd',sou, sou_insurf,dsou_insurf, &
00366           &                                                      sou_outsurf,dsou_outsurf,pot)
00367         end if
00368         if (ii.eq.1) potx(0:nrg,0:ntg,0:npg) = pot(0:nrg,0:ntg,0:npg)
00369         if (ii.eq.2) poty(0:nrg,0:ntg,0:npg) = pot(0:nrg,0:ntg,0:npg)
00370         if (ii.eq.3) potz(0:nrg,0:ntg,0:npg) = pot(0:nrg,0:ntg,0:npg)
00371       end do
00372 
00373       potx_bak(0:nrg,0:ntg,0:npg) = bvxd(0:nrg,0:ntg,0:npg)
00374       poty_bak(0:nrg,0:ntg,0:npg) = bvyd(0:nrg,0:ntg,0:npg)
00375       potz_bak(0:nrg,0:ntg,0:npg) = bvzd(0:nrg,0:ntg,0:npg)
00376       call update_grfield(potx,bvxd,conv_gra)
00377       call update_grfield(poty,bvyd,conv_gra)
00378       call update_grfield(potz,bvzd,conv_gra)
00379 
00380       if ((impt==1 .or. impt==2) .and. (iter_count>5 .or. iseq>1 .or. indata_type.eq.'3D')) &
00381         &          call update_parameter_BNS(conv_gra)
00382 
00383       if (impt.ne.nmpt) then 
00384         call interpolation_fillup_binary_COCP(impt,'bvxd','od',bvxd)
00385         call interpolation_fillup_binary_COCP(impt,'bvyd','od',bvyd)
00386         call interpolation_fillup_binary_COCP(impt,'bvzd','ev',bvzd)
00387       end if
00388 !      call reset_metric_CF
00389       call error_metric_type2_mpt(bvxd,potx_bak,error_tmp,flag_tmp,cocp,impt)
00390       flag_bvxd  =  max ( flag_bvxd, flag_tmp)
00391       error_bvxd = dmax1(error_bvxd,error_tmp)
00392       call error_metric_type2_mpt(bvyd,poty_bak,error_tmp,flag_tmp,cocp,impt)
00393       flag_bvyd  =  max ( flag_bvyd, flag_tmp)
00394       error_bvyd = dmax1(error_bvyd,error_tmp)
00395       call error_metric_type2_mpt(bvzd,potz_bak,error_tmp,flag_tmp,cocp,impt)
00396       flag_bvzd  =  max ( flag_bvzd, flag_tmp)
00397       error_bvzd = dmax1(error_bvzd,error_tmp)
00398 !!!      if(ishift==shift_iter)  write(6,'(1p,e14.6,a2)', ADVANCE = "NO")  error_bvxd, "  "
00399 !!!      if(ishift==shift_iter)  write(6,'(1p,e14.6,a2)', ADVANCE = "NO")  error_bvyd, "  "
00400 !!!      if(ishift==shift_iter)  write(6,'(1p,e14.6,a2)', advance = "no")  error_bvzd, "  "
00401       end do   ! ==========>> ishift 
00402 
00403 !      call printout_error_metric_v1('bvxd', iter_count,error_bvxd)
00404 !      call printout_error_metric_v1('bvyd', iter_count,error_bvyd)
00405 !      call printout_error_metric_v1('bvzd', iter_count,error_bvzd)
00406       call copy_to_mpatch_all_BNS_CF(impt)
00407       call copy_def_metric_and_matter_to_mpt(impt)
00408       call copy_def_matter_irrot_to_mpt(impt)
00409 !
00410       error_metr = dmax1(error_psi,error_alph,error_bvxd,error_bvyd,error_bvzd)
00411 
00412       if(iter_count<5 .and. iseq<=1 .and. indata_type.eq.'1D') then
00413         write(6,'(a14)') '#SKIP HYDRO a#' 
00414         error_emd_last = 1.0d0 
00415         cycle
00416       end if
00417 ! ---
00418 ! -- Hydro equations.
00419 !      if (error_emd_last.le.eps*0.0001d0 .and. error_metr.gt.eps*2.0d0) then
00420       if (error_emd_last.le.eps*0.0001d0 .and. error_vep_last.le.eps*0.0001d0 ) then
00421         write(6,'(a14)')  '#SKIP HYDRO b#'
00422       else
00423         if (impt==1 .or. impt==2) then
00424           do ihy = 1, hydro_iter
00425             call interpo_gr2fl_metric_CF
00426             call rotation_law
00427             call hydrostatic_eq_CF_peos_lecc_irrot(potf,potut,potux,potuy,potuz)
00428             call calc_surface(potrs, potf)
00429             emd_bak = emd 
00430             call update_matter( potf,emd,conv_den)
00431 !            call update_matter(potut,utf,conv_den)
00432 !            call update_matter(potux,vepxf,conv_den)
00433 !            call update_matter(potuy,vepyf,conv_den)
00434 !            call update_matter(potuz,vepzf,conv_den)
00435             call update_surface(potrs,rs,conv_den)
00436 
00437             call reset_fluid_BNS
00438             call update_parameter_BNS(conv_den)
00439             call calc_vector_x_matter(2)
00440             call calc_vector_phi_matter(2)
00441 !
00442             call hydro_irbns_vep_CF_peos_lecc(potf,iter_count,impt,10)
00443             vep_bak = vep
00444             call update_matter(potf,vep,conv_vep)
00445             call reset_fluid_vep
00446             call calc_gradvep(vep,vepxf,vepyf,vepzf)
00447             call reset_fluid_gradvep
00448             call update_parameter_BNS(conv_vep)
00449 
00450             if (ihy.eq.hydro_iter) then
00451 !            if (ihy.eq.1) then
00452               call error_matter_type2(emd,emd_bak,error_tmp,flag_tmp)
00453               flag_emd  =  max ( flag_emd, flag_tmp)
00454               error_emd = dmax1(error_emd,error_tmp)
00455               error_emd_last = error_emd
00456 !!!              write(6,'(1p,e14.6,a2)', ADVANCE = "NO")  error_emd, "  "
00457 !              call printout_error_matter("emd ", iter_count,error_emd)
00458               call error_velocity_potential(vep,vep_bak,error_tmp,flag_tmp)
00459               flag_vep  =  max ( flag_vep, flag_tmp)
00460               error_vep = dmax1(error_vep,error_tmp)
00461               error_vep_last = error_vep
00462 !!!              write(6,'(1p,e14.6,a2)')  error_vep, "  "        
00463             end if
00464           end do   ! end hydro
00465           call copy_to_mpatch_all_BNS_CF(impt)
00466           call copy_def_metric_and_matter_to_mpt(impt)
00467           call copy_def_matter_irrot_to_mpt(impt)
00468         else
00469           write(6,*)   ' ##   NO HYDRO ##'
00470         end if
00471       end if
00472 !      error_emd_last = error_emd
00473 !
00474       if(impt>5) then
00475         if (mod(iter_count,1)==0.or.iter_count==1) then
00476           write(char1, '(i5)') iter_count
00477           char2 = adjustl(char1)
00478           write(char4, '(i5)') impt
00479           char5 = adjustl(char4)
00480           char3 = 'iteration' // trim(char2) // '_mpt' // trim(char5) // '.txt'
00481           call printout_axis_irrot_BNS_mpt(impt,char3)
00482         endif
00483       end if  !  ====> iteration.txt
00484 
00485 !     print once the solution when the error in shift for all patches is <0.1 
00486       if (flag_mid==0) then
00487         if(error_metr < iter_eps ) then
00488           flag_mpt_mid(impt) = 1
00489           error_metric_mpt_mid(impt,1) = error_psi
00490           error_metric_mpt_mid(impt,2) = error_alph
00491           error_metric_mpt_mid(impt,3) = error_bvxd
00492           error_metric_mpt_mid(impt,4) = error_bvyd
00493           error_metric_mpt_mid(impt,5) = error_bvzd
00494         end if
00495         sumjj=0
00496         do jj=1,nmpt 
00497           sumjj = sumjj + flag_mpt_mid(jj)
00498         end do
00499         if (sumjj==nmpt) then
00500           open(3,file='mid_error_data.txt',status='unknown')
00501           write(3,'(a11,i5)') 'Iteration =', iter_count
00502           write(3,'(a8,1p,5e17.6)') 'COCP #=1', (error_metric_mpt_mid(1,je), je=1,5)
00503           write(3,'(a8,1p,5e17.6)') 'COCP #=2', (error_metric_mpt_mid(2,je), je=1,5)
00504           write(3,'(a8,1p,5e17.6)') 'ARCP #  ', (error_metric_mpt_mid(3,je), je=1,5)
00505           close(3)
00506 !          do jj=1,nmpt
00507 !            call copy_from_mpatch_all_BNS_CF(jj)
00508 !            call copy_def_metric_and_matter_from_mpt(jj)
00509 !            call copy_def_matter_irrot_from_mpt(jj)
00510 !            call IO_output_solution_3D_CF_irrot_NS_mpt(jj,'.mid')
00511 !            if (jj==1 .or. jj==2)  call printout_NS_shape_mpt(jj)
00512 !          end do
00513           flag_mid = 1
00514         end if 
00515       end if
00516     end do    !  ====================>> end all patches
00517 !
00518 ! -- convergence of iteration is checked.
00519 !
00520 !    if (iter_count==6) then
00521 !      do impt=1,nmpt
00522 !        call copy_from_mpatch_all_BNS_CF(impt)
00523 !        write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi
00524 !      end do
00525 !      write(2,*) '-----------------------------------------------------------------------'
00526 !      close(2)
00527 !      exit
00528 !    end if
00529 
00530     if (flag_mid==1) then
00531       do impt=1,nmpt
00532         call copy_from_mpatch_all_BNS_CF(impt)
00533 !        call copy_def_metric_and_matter_from_mpt(impt)
00534 !        call copy_def_matter_irrot_from_mpt(impt)
00535         write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi
00536 
00537 !        MoC_vio(0:nrg,0:ntg,0:npg,1:3) = 0.0d0
00538 !        call violation_gridpoint_MoC_CF_peos_irrot(MoC_vio,'ns')
00539 !        pot = 0.0d0
00540 !        pot(0:nrg,0:ntg,0:npg) = MoC_vio(0:nrg,0:ntg,0:npg,2)   
00541 !        char1 = 'MoC_vio_by_3D_mpt1.txt'
00542 !        call IO_output_3D_general(char1,'g',pot)
00543 !        char1 = 'MoC_vio_by_1D_mpt1.txt'
00544 !        call IO_output_1D_general(char1,'g',pot,-1,ntg/2,0)
00545 
00546 !        HaC_vio = 0.0d0
00547 !        call violation_gridpoint_HaC_CF_peos_irrot(HaC_vio,'ns')
00548 !        char1 = 'HaC_vio_3D_mpt1.txt'
00549 !        call IO_output_3D_general(char1,'g',HaC_vio)
00550 !        char1 = 'HaC_vio_1D_mpt1.txt'
00551 !        call IO_output_1D_general(char1,'g',HaC_vio,-1,ntg/2,0)
00552       end do
00553       write(2,*) '-----------------------------------------------------------------------'
00554       close(2)
00555       exit
00556     end if
00557 
00558     if ( (flag_vep==0).and.(flag_emd==0).and.(flag_psi==0).and.(flag_alph==0).and. &
00559     &   (flag_bvxd==0).and.(flag_bvyd==0).and.(flag_bvzd==0))      then
00560       do impt=1, nmpt
00561         call copy_from_mpatch_all_BNS_CF(impt)
00562         write(2,'(i5,i5,1p,3e20.12)') iter_count, impt, ome, ber, radi
00563       end do
00564       write(2,*) '-----------------------------------------------------------------------'
00565       close(2)
00566       exit
00567     endif
00568 
00569     if (iter_count >= iter_max) then
00570       write(6,*) "*************** iter_count==iter_max ....exiting *******************"
00571       stop
00572 !      exit
00573     end if
00574 !    error_emd_last = error_emd
00575 !
00576     flag_psi  = 0
00577     flag_alph = 0
00578     flag_bvxd = 0
00579     flag_bvyd = 0
00580     flag_bvzd = 0
00581     flag_emd  = 0
00582     flag_vep  = 0
00583     error_psi  = 0.0d0
00584     error_alph = 0.0d0
00585     error_bvxd = 0.0d0
00586     error_bvyd = 0.0d0
00587     error_bvzd = 0.0d0
00588     error_emd  = 0.0d0
00589     error_vep  = 0.0d0
00590   end do
00591 !
00592   deallocate(souvec)
00593   deallocate(sou)
00594   deallocate(pot)
00595   deallocate(potx)
00596   deallocate(poty)
00597   deallocate(potz)
00598   deallocate(pot_bak)
00599   deallocate(potx_bak)
00600   deallocate(poty_bak)
00601   deallocate(potz_bak)
00602   deallocate(work)
00603   deallocate(fnc_itp)
00604   deallocate(sou_exsurf)
00605   deallocate(dsou_exsurf)
00606   deallocate(sou_insurf)
00607   deallocate(dsou_insurf)
00608   deallocate(sou_outsurf)
00609   deallocate(dsou_outsurf)
00610   deallocate(potf)
00611   deallocate(potut)
00612   deallocate(potux)
00613   deallocate(potuy)
00614   deallocate(potuz)
00615   deallocate(emd_bak)
00616   deallocate(vep_bak)
00617   deallocate(potrs)
00618   
00619 end subroutine iter_lecc_irrot_BNS_CF_mpt
/home/uryu/Cocal_svn_checkinout/Cocal/code/Subroutine_mpatch/iter_lecc_irrot_BNS_CF_mpt.f90
