merge pull request #21

build/Makefile

@@ -864,10 +864,11 @@ endif
 
 ifeq ($(SETUP), jet)
 #   Jet simulation from Price, Tricco & Bate (2012)
-    FPPFLAGS= -DPERIODIC -DGRAVITY
     SETUPFILE= setup_sphereinbox.f90
     #ANALYSIS= analysis_jet.f90
     ANALYSIS= ${SRCNIMHD} analysis_protostar_environ.F90
+    PERIODIC=yes
+    GRAVITY=yes
     ISOTHERMAL=yes
     MHD=yes
     IND_TIMESTEPS=yes
@@ -963,7 +964,7 @@ endif
 
 ifeq ($(SETUP), testgrav)
 #   self-gravity unit tests
-    FPPFLAGS= -DGRAVITY
+    GRAVITY=yes
     KNOWN_SETUP=yes
     CONST_ARTRES=yes
     CURLV=yes

src/main/energies.F90

@@ -38,7 +38,7 @@ module energies
  real,            public    :: xmom,ymom,zmom
  real,            public    :: totlum
  integer,         public    :: iquantities
- integer(kind=8), public    :: ndead,np_cs_eq_0,np_e_eq_0
+ integer(kind=8), public    :: ndead,npartall,np_cs_eq_0,np_e_eq_0
  integer,         public    :: iev_time,iev_ekin,iev_etherm,iev_emag,iev_epot,iev_etot,iev_totmom,iev_com(3),&
                                iev_angmom,iev_rho,iev_dt,iev_dtx,iev_entrop,iev_rmsmach,iev_vrms,iev_rhop(6),&
                                iev_alpha,iev_B,iev_divB,iev_hdivB,iev_beta,iev_temp,iev_etaar,iev_etao(2),iev_etah(4),&
@@ -605,9 +605,10 @@ subroutine compute_energies(t)
 !$omp end parallel
 
  !--Determing the number of active gas particles
- nptot = reduce_fn('+',np)
- npgas = reduce_fn('+',npgas)
- ndead = npart - nptot
+ nptot    = reduce_fn('+',np)
+ npgas    = reduce_fn('+',npgas)
+ npartall = reduce_fn('+',npart)
+ ndead    = npartall - nptot
  if (nptot > 0) then
     dnptot = 1./real(nptot)
  else

src/main/evwrite.F90

@@ -47,12 +47,12 @@
 !+
 !--------------------------------------------------------------------------
 module evwrite
- use io,             only: fatal
+ use io,             only: fatal,iverbose
  use options,        only: iexternalforce
  use timestep,       only: dtmax_dratio
  use externalforces, only: iext_binary,was_accreted
  use energies,       only: inumev,iquantities,ev_data
- use energies,       only: ndead
+ use energies,       only: ndead,npartall
  use energies,       only: gas_only,track_mass,track_lum
  use energies,       only: iev_sum,iev_max,iev_min,iev_ave
  use energies,       only: iev_time,iev_ekin,iev_etherm,iev_emag,iev_epot,iev_etot,iev_totmom,iev_com,&
@@ -427,15 +427,12 @@ subroutine write_evlog(iprint)
  character(len=120)  :: string,Mdust_label(maxdusttypes)
  integer             :: i
 
-!***Uncomment this once debugging is complete
-! if (ndead > 0 .or. nptmass > 0 .or. icreate_sinks > 0 .or. particles_are_injected) then
-!    write(iprint,"(1x,4(a,I10))") 'npart=',npart,', n_alive=',npart-ndead,', n_dead_or_accreted=',ndead,', nptmass=',nptmass
-! endif
-!***Remove the following once debugging is complete
- write(iprint,"(1x,4(a,I10))") 'npart=',npart,', n_alive=',npart-ndead,', n_dead_or_accreted=',ndead,', nptmass=',nptmass
+ if (ndead > 0 .or. nptmass > 0 .or. icreate_sinks > 0 .or. particles_are_injected .or. iverbose > 0) then
+    write(iprint,"(1x,4(a,I10))") 'npart=',npartall,', n_alive=',npartall-ndead, &
+                                  ', n_dead_or_accreted=',ndead,', nptmass=',nptmass
+ endif
 
- write(iprint,"(1x,3('E',a,'=',es10.3,', '),('E',a,'=',es10.3))") &
-      'tot',etot,'kin',ekin,'therm',etherm,'pot',epot
+ write(iprint,"(1x,3('E',a,'=',es10.3,', '),('E',a,'=',es10.3))") 'tot',etot,'kin',ekin,'therm',etherm,'pot',epot
 
  if (mhd)          write(iprint,"(1x,('E',a,'=',es10.3))") 'mag',emag
  if (do_radiation) write(iprint,"(1x,('E',a,'=',es10.3))") 'rad',erad

src/main/phantom.F90

@@ -34,13 +34,10 @@ program phantom
  use io,              only:id,master,nprocs,set_io_unit_numbers,die
  use evolve,          only:evol
  implicit none
- integer :: nargs,i,ntests,npass,nfail
+ integer            :: nargs
  character(len=120) :: infile,logfile,evfile,dumpfile
 
  id = 0
- ntests = 0
- npass  = 0
- nfail  = 0
 
  call init_mpi(id,nprocs)
  call set_io_unit_numbers

src/main/utils_indtimesteps.F90

@@ -385,17 +385,17 @@ subroutine write_binsummary(npart,nbinmax,dtmax,timeperbin,iphase,ibin,xyzh)
     endif
  enddo over_part
 
- !ninbin(:) = int(reduce_mpi('+',ninbin(:)))
+ ninbin(:) = reduce_mpi('+',ninbin(:))
  ntypesprint = 0
  itypelist(:) = 0
  do itype=1,maxtypes
-    !noftypeinbin(:,itype) = int(reduce_mpi('+',noftypeinbin(:,itype)))
+    noftypeinbin(:,itype) = reduce_mpi('+',noftypeinbin(:,itype))
     if (any(noftypeinbin(:,itype) > 0)) then
        ntypesprint = ntypesprint + 1
        itypelist(ntypesprint) = itype
     endif
  enddo
- ntot = np !reduce_mpi('+',np)
+ ntot = reduce_mpi('+',np)
  timeperbintot(:) = reduce_mpi('+',timeperbin(:))
 
  if (id==master) then

src/tests/test_nonidealmhd.F90

@@ -402,7 +402,7 @@ subroutine test_standingshock(ntests,npass)
  gamma    =  1.0
  dt       =  2.0d-3
  nsteps   =  500
- t        =  0
+ t        =  0.
  iverbose =  0
  tmax     = nsteps*dt
  eta_const_type           = icnstsemi
@@ -610,10 +610,10 @@ subroutine test_etaval(ntests,npass)
  ion_rays     = .true.
  ion_thermal  = .true.
  use_sts      = .false.
- itmp  = 1 ! avoids compiler warning
+ itmp         = 1 ! avoids compiler warning
 
  ! initialise eos, & the Nicil library
- call init_eos(8,ierr)
+ call init_eos(ieos,ierr)
  call nicil_initialise(real(utime),real(udist),real(umass),real(unit_Bfield),ierr)
  !
  !--Loop over both sets of calculations
@@ -649,8 +649,8 @@ subroutine test_etaval(ntests,npass)
     tempi = get_temperature(ieos,xyzh(1:3,itmp),rhoi,vxyzu(:,itmp))
 
     print*, ' '
-    write(*,'(1x,a,3Es18.11)') 'Used   rho,B_z,temp (cgs): ',rhoi*unit_density,Bi*unit_Bfield,tempi
-    write(*,'(1x,a,3Es18.11)') 'eta_ohm, eta_hall, eta_ambi (cgs): ', eta_nimhd(1:3,itmp)*unit_eta
+    write(*,'(1x,a,3Es18.10)') 'Used   rho,B_z,temp (cgs): ',rhoi*unit_density,Bi*unit_Bfield,tempi
+    write(*,'(1x,a,3Es18.10)') 'eta_ohm, eta_hall, eta_ambi (cgs): ', eta_nimhd(1:3,itmp)*unit_eta
     call checkval(eta_nimhd(iohm, itmp)*unit_eta,eta_act(1,k),tol,nerr(3*(k-1)+1),'calculated non-constant eta_ohm')
     call checkval(eta_nimhd(ihall,itmp)*unit_eta,eta_act(2,k),tol,nerr(3*(k-1)+2),'calculated non-constant eta_hall')
     call checkval(eta_nimhd(iambi,itmp)*unit_eta,eta_act(3,k),tol,nerr(3*(k-1)+3),'calculated non-constant eta_ambi')

src/utils/analysis_MWpdf.f90

@@ -45,8 +45,9 @@ subroutine do_analysis(dumpfile,num,xyzh,vxyzu,particlemass,npart,time,iunit)
  real :: vx2,vy2,vz2,rhoi,rmsv
  real :: rhomean, rho(maxp), hi, pmassi
 
- rhomin = huge(rhomin)
- rhomax = 0.
+ rho     = 0.
+ rhomin  = huge(rhomin)
+ rhomax  = 0.
  rhomean = 0.
  totmass = 0.
  print*,'hfact = ',hfact

src/utils/analysis_clumpfind.F90

@@ -411,7 +411,7 @@ subroutine read_analysis_options(dumpfile)
     write(10,*) trim(dumpfile), "    Previous SPH dump analysed"
     close(10)
     previousdumpfile = ""
-    runningclumpmax  = ""
+    runningclumpmax  = 0
  endif
 
  if (checkbound) then