Cocal: /home/uryu/Cocal_svn_checkinout/Cocal/code/Main_code/Main_irrot_BNS_CF

00001 !______________________________________________
00002 include '../Include_file/include_modulefiles_BNS_CF_3mpt.f90'
00003 include '../Include_file/include_PEOS_modulefile.f90'
00004 include '../Include_file/include_interface_modulefiles_BNS_CF_3mpt.f90'
00005 include '../Include_file/include_subroutines_BNS_CF_3mpt.f90'
00006 include '../Include_file/include_PEOS_subroutines.f90'
00007 include '../Include_file/include_functions.f90'
00008 !______________________________________________
00009 !
00010 !              Main Program
00011 !______________________________________________
00012 PROGRAM Main_irrot_BNS_CF_mpt
00013 !
00014   use phys_constant, only : long, nmpt
00015   use grid_parameter
00016   use def_matter
00017   use def_matter_parameter, only : emdc, radi, ome
00018   use def_matter_parameter_mpt
00019   use def_quantities_mpt
00020   use def_peos_parameter, only : rhoini_cgs, emdini_gcm1
00021 !  use def_bh_parameter, only : bh_soltype
00022   use grid_parameter_binary_excision
00023   use grid_points_binary_excision
00024   use grid_points_asymptotic_patch
00025   use grid_points_binary_in_asympto
00026   use weight_midpoint_binary_excision
00027   use radial_green_fn_grav
00028   use radial_green_fn_grav_bhex_nb
00029   use radial_green_fn_grav_bhex_dd
00030 !  use radial_green_fn_grav_bhex_di
00031   use radial_green_fn_grav_bhex_nd
00032   use radial_green_fn_grav_bhex_dh
00033   use radial_green_fn_grav_bhex_nh
00034   use radial_green_fn_grav_bhex_sd
00035   use interface_copy_to_hgfn_and_gfnsf
00036   use interface_calc_gradvep
00037   use interface_calc_gradvep_from_corot_id
00038   implicit none
00039 
00040 !------------------------------------------------------------------------------------------------------------!
00041 !                                                                                                            !
00042 ! numseq=1 for a single solution with given TOV initial data or given 3D initial data                        !
00043 ! numseq>1 produces a sequence of solutions with increasing emdc using routine next_solution_BNS_mpt.f90     !
00044 ! In this case the iteration starts from TOV initial data (or 3D) computes the first irrotating solution     !
00045 ! and then increases emdc to move to the next irrotating solution.                                           !
00046 ! When computing parallel irrotating solutions starting from different TOV use numseq=1 and the script       !
00047 ! debug_BNS_CF_peos_seq_from_TOV.sh                                                                          !
00048 !                                                                                                            !
00049 ! sw_mass_iter = 'n' implements the above scenarios                                                          !
00050 ! sw_mass_iter = 'y' finds a solution for a given rest mass.                                                 !
00051 !                                                                                                            !
00052 ! print_sol_seq=0 for NOT printing solution data files                                                       !
00053 !                                                                                                            !
00054 !------------------------------------------------------------------------------------------------------------!
00055 
00056   integer :: impt, itg, print_sol_seq=1, numseq=1
00057   integer :: iseq, iter_count, total_iteration
00058   character(30) :: char1, char2, char3, char4, char5
00059   character(100) :: dircommand
00060   real(long) :: iter_eps
00061 !
00062 ! -- Read parameters
00063   call allocate_grid_parameter_mpt
00064   call allocate_grid_parameter_binary_excision_mpt
00065   call allocate_def_matter_parameter_mpt
00066   do impt = 1, nmpt
00067     call read_parameter_mpt(impt)
00068     call read_surf_parameter_mpt(impt)
00069     call copy_grid_parameter_to_mpt(impt)
00070     call read_parameter_binary_excision_mpt(impt)
00071     call copy_grid_parameter_binary_excision_to_mpt(impt)
00072     if (impt==1 .or. impt==2) then
00073       call peos_initialize_mpt(impt)
00074       call copy_def_peos_parameter_to_mpt(impt)
00075     end if
00076     call copy_grid_parameter_to_mpt(impt)
00077   end do
00078 
00079 !stop
00080 
00081 !
00082 ! -- Allocate arrays
00083   call set_allocate_size_mpt
00084 !
00085   call allocate_coordinate_patch_kit_grav_mpt
00086   call allocate_grid_points_binary_excision
00087   call allocate_grid_points_asymptotic_patch
00088   call allocate_grid_points_binary_in_asympto
00089   call allocate_weight_midpoint_binary_excision
00090   call allocate_hgfn_bhex
00091   call allocate_hgfn_bhex_nb
00092   call allocate_hgfn_bhex_dd
00093   call allocate_hgfn_bhex_nd
00094 !  call allocate_hgfn_bhex_di
00095   call allocate_hgfn_bhex_sd
00096   call allocate_BNS_CF
00097   call allocate_metric_on_SFC_CF
00098   call allocate_BNS_CF_irrot
00099 !
00100   call allocate_mpatch_all_BNS_CF
00101   call allocate_grid_points_asymptotic_patch_mpt
00102   call allocate_grid_points_binary_in_asympto_mpt
00103   call allocate_BNS_CF_mpt
00104   call allocate_BNS_CF_irrot_mpt
00105 !
00106   do impt = 1, nmpt
00107     call copy_grid_parameter_from_mpt(impt)
00108     call copy_grid_parameter_binary_excision_from_mpt(impt)
00109     call copy_def_peos_parameter_from_mpt(impt)
00110     write(6,*) "*****************************************************************************"
00111     write(6,*) "* Change integer input (1,2,3,4) ALSO in following subroutines:             *"
00112     write(6,*) "* code/Main_utility/interpolation_contour_potential_irrot_BNS_CF_3mpt.f90   *"
00113     write(6,*) "* spherical_data/initial_isotropic_coord/code/Import_isotr_peos_surf.f90    *"
00114     write(6,*) "*****************************************************************************"
00115     call coordinate_patch_kit_grav_grid_mpt(3)       ! 3:r_surf is used,  4:constant dr until r=7
00116 
00117     call calc_parameter_binary_excision
00118     call IO_printout_grid_data_mpt(impt)
00119     if (impt.ne.nmpt) then
00120       call calc_grid_points_binary_excision
00121       call calc_grid_points_binary_in_asympto(impt,nmpt)
00122       call copy_grid_points_binary_in_asympto_to_mpt(impt)
00123     end if
00124     call calc_weight_midpoint_binary_excision
00125 !      call calc_weight_midpoint_binary_excision_hybrid
00126     if (impt==1 .or. impt==2) then
00127       call calc_vector_x_grav(2)
00128       call calc_vector_phi_grav(2)
00129     else
00130       call calc_vector_x_grav(1)
00131       call calc_vector_phi_grav(1)
00132     end if
00133 !    call calc_vector_x_matter(2)
00134 !    call calc_vector_phi_matter(2)
00135     call copy_to_mpatch_all_BNS_CF(impt)
00136     call copy_def_metric_and_matter_to_mpt(impt)
00137     call copy_def_matter_irrot_to_mpt(impt)
00138   end do
00139   call copy_from_mpatch_all_BNS_CF(nmpt)
00140 ! -- coordinates of asymptotic patch in central patches
00141   do impt = 1, 2
00142     call calc_grid_points_asymptotic_patch(impt,nmpt)
00143     call copy_grid_points_asymptotic_patch_to_mpt(impt)
00144   end do
00145 !
00146 ! Initialization of variables
00147   do impt = 1, nmpt
00148     call copy_from_mpatch_all_BNS_CF(impt)
00149     call copy_def_metric_and_matter_from_mpt(impt)
00150     call copy_def_matter_irrot_from_mpt(impt)
00151     if (indata_type.eq.'3D') then
00152       if(NS_shape.eq.'IR')  then
00153         write(6,*) "########################### Reading 3D Irrotational Initial Data ####################################"
00154         call IO_input_initial_3D_CF_irrot_NS_mpt(impt)
00155 !       if (impt==1 .or. impt==2)     emdc = emdini_gcm1    ! from peos_parameter_mpt file
00156         if (impt==1 .or. impt==2)    then
00157           emdc = emd(0,0,0)
00158           call calc_gradvep(vep,vepxf,vepyf,vepzf)
00159 !          ======= check initial xaxis
00160 !          write(char1, '(i5)') 0
00161 !          char2 = adjustl(char1)
00162 !          write(char4, '(i5)') impt
00163 !          char5 = adjustl(char4)
00164 !          char3 = 'iteration' // trim(char2) // '_mpt' // trim(char5) // '.txt'
00165 !          call printout_axis_irrot_BNS_mpt(impt,char3)
00166         else
00167           emdc = 0.0d0
00168         end if 
00169       else if(NS_shape.eq.'CO')  then
00170         write(6,*) "########################### Reading 3D Corotating Initial Data ####################################"
00171         call IO_input_initial_3D_CF_NS_mpt(impt)
00172         if (impt==1 .or. impt==2)   then
00173           emdc = emd(0,0,0)
00174           call calc_gradvep_from_corot_id(vep,vepxf,vepyf,vepzf)
00175         else
00176           emdc = 0.0d0
00177         end if
00178       else
00179         write(6,*) "**************** indata_type=3D and NS_shape is neither IR nor CO....exiting"
00180         stop
00181       end if
00182     else
00183       if (impt==1 .or. impt==2)  then
00184         call  IO_input_initial_1D_CF_NS_mpt(impt) 
00185         emdc = emd(0,0,0)
00186         call initial_velocity_potential_NS_mpt(impt)
00187       else
00188         emdc = 0.0d0
00189       end if
00190     end if
00191     write (6,'(a8,i2,a16,1p,e20.12)')  "--Patch=",impt,"           emdc=",emdc
00192     if (impt==1 .or. impt==2)  then
00193 !      call initial_velocity_potential_NS_mpt(impt)
00194 !      call initial_vrot_NS_mpt(impt)
00195     end if
00196     call copy_def_metric_and_matter_to_mpt(impt)
00197     call copy_def_matter_irrot_to_mpt(impt)
00198     call copy_to_mpatch_all_BNS_CF(impt)
00199   end do  ! =========>>> end initialization of all patches
00200 !
00201 
00202   call copy_grid_parameter_from_mpt(1)  
00203   if ( sw_mass_iter=='y' ) then
00204     call calc_irrot_BNS_CF_mpt(total_iteration)
00205     call calc_physical_quantities_irrot_BNS_CF_mpt
00206     char3 = 'main_bnsphys_all_mpt.txt'
00207     call printout_physq_BNS_all_mpt(char3)
00208     call write_last_physq_BNS_mpt
00209     do impt = 1, nmpt
00210       call copy_from_mpatch_all_BNS_CF(impt)
00211       call copy_def_metric_and_matter_from_mpt(impt)
00212       call copy_def_matter_irrot_from_mpt(impt)
00213       if (outdata_type.eq.'3D')  call IO_output_solution_3D_CF_irrot_NS_mpt(impt,'.las')
00214       if (impt==1 .or. impt==2)  call printout_NS_shape_mpt(impt)
00215     end do
00216   else
00217     do iseq = 1, numseq   !num_sol_seq
00218       write(char1, '(i5)') iseq
00219       char2 = adjustl(char1)
00220       if (iseq < 10)  then
00221         char3 = 'iseq0' // trim(char2)
00222       else
00223         char3 = 'iseq'  // trim(char2)
00224       end if
00225       dircommand = 'mkdir ' // char3
00226       call system(dircommand)
00227 !      dircommand = 'cp -f peos_parameter_mpt* ./' // char3
00228 !      call system(dircommand)
00229       call chdir(char3)
00230 !
00231       write(6,'(a50,i2,a32)') '============================== Solution sequence #',iseq,' ============================== '
00232       write(6,'(a14,1p,2e20.12)')  "emdc COCP1,2: ", def_matter_param_real_(2,1), def_matter_param_real_(2,2)
00233 !     Change the argument iseq to a number bigger than 1 (ex 2) if you want to  compute hydro from first iteration. 
00234 !     Usually this is the case when the initial data are 3D.
00235       iter_eps = 1.0d-02
00236       call iter_irrot_BNS_CF_mpt(iter_count, iseq, iter_eps)
00237       call calc_physical_quantities_irrot_BNS_CF_mpt
00238       char3 = 'main_bnsphys_all_mpt.txt'
00239       call printout_physq_BNS_all_mpt(char3)
00240       call write_last_physq_BNS_mpt
00241       if (print_sol_seq==1) then
00242         do impt = 1, nmpt
00243           call copy_from_mpatch_all_BNS_CF(impt)
00244           call copy_def_metric_and_matter_from_mpt(impt)
00245           call copy_def_matter_irrot_from_mpt(impt)
00246           if (outdata_type.eq.'3D') call IO_output_solution_3D_CF_irrot_NS_mpt(impt,'.las')
00247           if (impt==1 .or. impt==2)  call printout_NS_shape_mpt(impt)
00248         end do
00249       end if
00250       call chdir('../')   ! now we are in work_area_BNS
00251       if ( iseq < numseq )   call next_solution_BNS_mpt(iseq)
00252     end do
00253   end if 
00254 !
00255 END PROGRAM Main_irrot_BNS_CF_mpt
/home/uryu/Cocal_svn_checkinout/Cocal/code/Main_code/Main_irrot_BNS_CF_mpt.f90
