!####################################################################### ! ### ! Main Program ### ! ### ! mod 2011.08.21 ### !####################################################################### program main use mod_parameter use Mod_Variable implicit none !---- Variables -------------------------------------------------------- real(8) :: time ! Time integer(4) :: i ! Particle ID integer(4) :: k ! x,y,z components real(8) :: dt ! Time Step integer(4) :: n_solve ! Iteretion Count integer(4) :: n_step ! Number of Time Steps !----------------------------------------------------------------------- real(8) :: tot_eng0=0.0d0 ! Initial Total Energy real(8) :: kin_eng0 ! Initial Total Kinetic Energy real(8) :: int_eng0 ! Initial Total Internal Energy real(8) :: pot_eng0 ! Initial Total Potential Energy real(8) :: tot_eng ! Total Energy real(8) :: kin_eng ! Total Kinetic Energy real(8) :: int_eng ! Total Internal Energy real(8) :: pot_eng ! Total Potential Energy real(8) :: err_tot_eng ! Error of Total Energy real(8) :: xc(3) ! Position of Center of Mass (x,y,z) real(8) :: vc(3) ! Velocity of Center of Mass (x,y,z) real(8) :: amom(3)=0.0d0 ! Total Angular Momentum (x,y,z) real(8) :: amom0(3)=0.0d0 ! Initial amom(1-3) real(8) :: amom0_max=0.0d0 ! Max of amom0(1-3) real(8) :: err_amom(3) ! Error of Total Angular Momentum (x,y,z) real(8) :: time_out integer(4) :: flag_out integer(4) :: n_out real(8) :: m_tot,m_ave !---- Variables for Calculation Statistics ----------------------------- real(8) :: h_ite_ave1=0.0d0 ! Average Iteration Number for Calc h integer(4) :: h_ite_max1=0 ! Max Iteration Number for Calc h real(8) :: np_ave_g_tree real(8) :: np_ave_g_tree1=0.0d0 !---- Variables for calculating CPU TIME ------------------------------- real(4) :: etime,te(2) real(8) :: cpu_tot,cpu_tot1,cpu_tot2 real(8) :: cpu_tree = 0.0d0 real(8) :: cpu_h = 0.0d0 real(8) :: cpu_d = 0.0d0 real(8) :: cpu_p = 0.0d0 real(8) :: cpu_a = 0.0d0 real(8) :: cpu_g = 0.0d0 real(8) :: cpu_t = 0.0d0 real(8) :: cpu_t1,cpu_t2 real(8) :: cpu_h1,cpu_h2 real(8) :: cpu_d1,cpu_d2 real(8) :: cpu_p1,cpu_p2 real(8) :: cpu_a1,cpu_a2 real(8) :: cpu_g1,cpu_g2 real(8) :: cpu_tree1,cpu_tree2 !---- Variables for Clock Time ----------------------------------------- integer(4) :: t1,t2,t_rate,t_max,t_diff !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% real(8) :: tmp real(8) :: dr_tb integer(4) :: i_tb,i_tb_max real(8) :: r_tb(0:1000),e_tb(0:1000) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% character(50) :: com_mkdir character(9) :: DIR_OUT = 'DATA_PTCL' com_mkdir='mkdir '//DIR_OUT call system(com_mkdir) !======================================================================= !==== Main Program ===================================================== !======================================================================= call system_clock(t1) cpu_tot1=etime(te) !==== Open Files ======================================================= open(90,file='STEP_INFO.dat') open(91,file='NUM_ERROR.dat') open(92,file='ALERT.dat') !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if(FLAG_SAUMON_EOS == 1) then call sub_eos_saumon_read_table endif if(FLAG_5PHASE_EOS == 1) then call sub_eos_5phase_read_table endif !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if(FLAG_WADA_INI == 1) then open(40,file=FILE_TABLE_WADA) read(40,'(e15.7)') dr_tb i_tb=0 400 continue read(40,'(3e15.7)') r_tb(i_tb),tmp,e_tb(i_tb) if (tmp .gt. 1.0d-5) then i_tb=i_tb+1 go to 400 endif i_tb_max=i_tb close(40) endif !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !==== Read Initial Particle Data ======================================= call sub_read_initial_data !==== Create rho-Ec Table ============================================== call sub_create_rho_ec_table !==== Create t-E Table ================================================= call sub_create_t_e_table !==== Initialization =================================================== time=time_sta n_step=0 n_solve=1 time=time_sta time_out=time_sta+(TIME_END-time_sta)/dble(OUTP_NUM) flag_out=1 n_out=0 nb=NB_STD !---- Setting of h_max (max smoothing length) -------------------------- h_max=0.0d0 if(FLAG_H_MAX == 1) then !---- Total Mass ---------- m_tot=sum(m) !-------------------------- m_ave=m_tot/dble(np) h_max=0.5d0*(3.0d0/4.0d0/PAI/DENS_2H_LIM*m_ave)**(1.0d0/3.0d0) endif !======================================================================= !==== Start of Time Evolution ========================================== !======================================================================= 555 continue !==== Making Tree List ================================================= cpu_tree1=etime(te) call sub_create_tree_list cpu_tree2=etime(te) cpu_tree = cpu_tree & +cpu_tree2-cpu_tree1 !==== Calculating Smoothing Length ===================================== cpu_h1=etime(te) if(FLAG_hg_CONST == 1) then call sub_calc_smoothing_length endif if(FLAG_SOLVE .eq. 1 .and. n_solve .eq. 1) then h_ite_ave1=h_ite_ave h_ite_max1=h_ite_max endif cpu_h2=etime(te) cpu_h =cpu_h+cpu_h2-cpu_h1 !==== Reconstructing Tree List ========================================= cpu_tree1=etime(te) call sub_reconst_tree_list cpu_tree2=etime(te) cpu_tree = cpu_tree & +cpu_tree2-cpu_tree1 !==== Calculating Density ============================================== cpu_d1=etime(te) call sub_calc_density cpu_d2=etime(te) cpu_d =cpu_d+cpu_d2-cpu_d1 !==== Calculating Pressure & Sound Velocity ============================ cpu_p1=etime(te) call sub_calc_pressure_sound_velocity(flag_out) cpu_p2=etime(te) cpu_p =cpu_p+cpu_p2-cpu_p1 !==== Calculating Temperature ========================================== if(n_solve .eq. 1) call sub_calc_temperature_tillotson !==== Output Data ====================================================== if (flag_out .eq. 1 .and. n_solve .eq. 1) then !==== Calculating Kinetic and Internal Energies ==== call sub_calc_energies(tot_eng,kin_eng,int_eng) !==== Calculating Potential Energy ================= pot_eng=0.0d0 !---- Case for Direct Calculation ----------- if (FLAG_GRAV .eq. 1) call sub_calc_pot_energy_direct(pot_eng) !---- Case for Tree Method ------------------ if (FLAG_GRAV .eq. 2) call sub_calc_pot_energy_tree(pot_eng) tot_eng=tot_eng+pot_eng !==== Calculating Position and Velocity of COM ===== call sub_calc_COM(xc,vc) !==== Calculating Angular Momentum ================= call sub_calc_angular_momentum(xc,vc,amom) if (n_out .eq. 0) then !---- Initial Energies ----------------- tot_eng0=tot_eng kin_eng0=kin_eng int_eng0=int_eng pot_eng0=pot_eng amom0(1)=amom(1) amom0(2)=amom(2) amom0(3)=amom(3) amom0_max=dmax1(dabs(amom(1)), & & dabs(amom(2)), & & dabs(amom(3))) if(amom0_max .eq. 0.0d0) amom0_max=1.0d100 !--------------------------------------- !---- Output Header -------------------- write(91,'(a7,$)') '# OUTP' write(91,'(a7,$)') ' STEP' write(91,'(a13,$)') ' TIME[s]' write(91,'(a13,$)') ' ENERGY ERROR' write(91,'(a13,$)') ' TOT ENG [J]' write(91,'(a13,$)') ' KINT ENG [J]' write(91,'(a13,$)') ' INTR ENG [J]' write(91,'(a13,$)') ' POT ENG [J]' write(91,'(3a13,$)') ' COM X [m]', & & ' COM Y [m]', & & ' COM Z [m]' write(91,'(3a13,$)') ' COM VZ [m/s]', & & ' COM VY [m/s]', & & ' COM VZ [m/s]' write(91,'(3a13,$)') ' AMOM ERR X', & & ' AMOM ERR Y', & & ' AMOM ERR Z' write(91,'(3a13)') ' AMOM X [SI]', & & ' AMOM Y [SI]', & & ' AMOM Z [SI]' !--------------------------------------- endif !---- Calculation of Energy Error ---------------- err_tot_eng=dabs((tot_eng0-tot_eng)/tot_eng0) !---- Calculation of Angular Momentum Error ------ err_amom(1)=dabs((amom0(1)-amom(1))/amom0_max) err_amom(2)=dabs((amom0(2)-amom(2))/amom0_max) err_amom(3)=dabs((amom0(3)-amom(3))/amom0_max) !---- Output Data -------------------------------- write(91,'(i7,$)') n_out write(91,'(i7,$)') n_step write(91,'(e13.5,$)') time write(91,'(e13.5,$)') err_tot_eng write(91,'(4e13.5,$)') tot_eng,kin_eng,int_eng,pot_eng write(91,'(3e13.5,$)') xc(1),xc(2),xc(3) write(91,'(3e13.5,$)') vc(1),vc(2),vc(3) write(91,'(3e13.5,$)') err_amom(1),err_amom(2),err_amom(3) write(91,'(3e13.5)') amom(1),amom(2),amom(3) call flush(91) !------------------------------------------------- call sub_output(n_out,time) !---- End of Output ----- flag_out=0 n_out=n_out+1 !------------------------ endif !==== Calculating Acceleration ========================================= cpu_a1=etime(te) call sub_calc_accel cpu_a2=etime(te) cpu_a =cpu_a+cpu_a2-cpu_a1 !==== Calculating Gravitational Acceleration =========================== cpu_g1=etime(te) !---- Case for Direct Calculation -------------- if (FLAG_GRAV .eq. 1) call sub_calc_gravity_direct !---- Case for Tree Method --------------------- if (FLAG_GRAV .eq. 2) then call sub_calc_gravity_tree(np_ave_g_tree) if(FLAG_SOLVE .eq. 1 .and. n_solve .eq. 1) then np_ave_g_tree1=np_ave_g_tree endif endif !---- Add Gravitational Acceleration ----------- if (FLAG_GRAV .ne. 0) then do i=1,np do k=1,3 a(k,i)=a(k,i)+GC*a_grav(k,i) enddo enddo endif cpu_g2=etime(te) cpu_g =cpu_g+cpu_g2-cpu_g1 !==== Calculating Time Step ============================================ cpu_t1=etime(te) if (n_solve .eq. 1) then call sub_calc_time_step(dt) !---- Ajustment of Time Step ------------------------ if (time+dt .ge. time_out) then flag_out=1 dt=time_out-time time_out=time_sta & & +(TIME_END-time_sta)/dble(OUTP_NUM)*dble(n_out+1) endif !---------------------------------------------------- endif cpu_t2=etime(te) cpu_t =cpu_t+cpu_t2-cpu_t1 !==== Update x, v, e =================================================== if(FLAG_WADA_INI == 1) then call sub_next_step_e_table(n_solve,dt, & & i_tb_max,dr_tb,r_tb,e_tb) else call sub_next_step(n_step,n_solve,dt) endif !==== Update the Time ================================================== if (FLAG_SOLVE .eq. 0) time=time+dt if (FLAG_SOLVE .eq. 1) then if (n_solve .eq. 2) time=time+dt endif !==== Determination of Iteration ================== if (FLAG_SOLVE .eq. 1 .and. n_solve .eq. 1) then n_solve=2 go to 555 endif !================================================== n_step=n_step+1 !==== Output Time Step Data ============================================ !---- Output to Screen ------------------------------------------------- if (FLAG_DISP_TIME .eq. 1) then !---- Output Header --------------------- if (mod(n_step,10) .eq. 1) then write(*,'(a7,$)') '# STEP' write(*,'(a13,$)') ' TIME' write(*,'(a13,$)') ' DT' write(*,'(a13,$)') ' ENERGY ERROR' write(*,*) ' ' endif !---- Output Data ----------------------- write(*,'(i7,$)') n_step write(*,'(e13.5,$)') time write(*,'(e13.5,$)') dt write(*,'(e13.5,$)') err_tot_eng write(*,*) ' ' endif !---- Output to File --------------------------------------------------- if(FLAG_GRAV .eq. 0) np_ave_g_tree=0.0d0 if(FLAG_GRAV .eq. 1) np_ave_g_tree=dble(np) if(FLAG_SOLVE .eq. 1) then h_ite_ave=0.5d0*(h_ite_ave1+h_ite_ave) h_ite_max=max0(h_ite_max1,h_ite_max) if(FLAG_GRAV .eq. 2) & & np_ave_g_tree=0.5d0*(np_ave_g_tree1+np_ave_g_tree) endif !---- Output Header ------------------------ if (n_step .eq. 1) then write(90,'(a7,$)') '# STEP' write(90,'(a13,$)') ' TIME' write(90,'(a13,$)') ' DT' write(90,'(a13,$)') ' ENERGY ERROR' write(90,'(a13,$)') ' AVE h ITE' write(90,'(a7,$)') ' MAX' write(90,'(a10,$)') ' AVE GNODE' write(90,'(a8,$)') ' EFF[%]' write(90,*) ' ' endif !---- Output Data -------------------------- write(90,'(i7,$)') n_step write(90,'(e13.5,$)') time write(90,'(e13.5,$)') dt write(90,'(e13.5,$)') err_tot_eng write(90,'(e13.5,$)') h_ite_ave write(90,'(i7,$)') h_ite_max write(90,'(f10.1,$)') np_ave_g_tree write(90,'(f8.3,$)') np_ave_g_tree/dble(np)*100.0d0 write(90,*) ' ' call flush(90) !======================================================================= if (FLAG_SOLVE .eq. 1 .and. n_solve .eq. 2) n_solve=1 !==== Determination of End ============================================= if (time .le. TIME_END) go to 555 call system_clock(t2,t_rate,t_max) if(t2 < t1) then t_diff=t_max-t1+t2 else t_diff=t2-t1 endif !======================================================================= !==== End of Time Evolution ============================================ !======================================================================= cpu_tot2=etime(te) cpu_tot =cpu_tot2-cpu_tot1 !==== Analyze and Output of CPU TIME =================================== call sub_output_CPU(t_diff/dble(t_rate),cpu_tot,cpu_tree,cpu_h, & & cpu_d,cpu_p,cpu_a,cpu_g,cpu_t) close(90) close(91) close(92) stop end program main