Merge pull request #436 from danieljprice/growth

build/Makefile

@@ -96,11 +96,11 @@ ifndef SETUPFILE
 endif
 
 ifndef SRCNIMHD
-    SRCNIMHD = nicil.F90 nicil_supplement.F90
+    SRCNIMHD = nicil.F90 nicil_supplement.f90
 endif
 
 ifndef SRCDUST
-    SRCDUST = dust.F90 growth.F90
+    SRCDUST = dust.f90 growth.F90
 endif
 
 ifdef SMOL
@@ -471,13 +471,13 @@ SRCPOTS= extern_corotate.f90 \
          extern_binary.f90 \
          extern_spiral.f90 \
          extern_lensethirring.f90 \
-         extern_gnewton.F90 \
-         lumin_nsdisc.F90 extern_prdrag.F90 \
+         extern_gnewton.f90 \
+         lumin_nsdisc.F90 extern_prdrag.f90 \
          extern_Bfield.f90 \
          extern_densprofile.f90 \
          extern_staticsine.f90 \
          extern_gwinspiral.f90 \
-         externalforces.F90
+         externalforces.f90
 endif
 ifeq (X$(SRCPOT), X)
 SRCPOT=${SRCPOTS}
@@ -486,7 +486,7 @@ endif
 # metrics for GR
 #
 ifeq ($(GR),yes)
-    SRCPOT=extern_gr.F90 $(SRCPOTS:externalforces.F90=externalforces_gr.F90)
+    SRCPOT=extern_gr.F90 $(SRCPOTS:externalforces.f90=externalforces_gr.f90)
 endif
 ifdef METRIC
     SRCMETRIC= metric_${METRIC}.f90
@@ -505,12 +505,12 @@ SRCCHEM= fs_data.f90 mol_data.f90 utils_spline.f90 \
          cooling_molecular.f90 \
          cooling_functions.f90 \
          cooling_solver.f90 \
-         h2chem.f90 cooling.F90
+         h2chem.f90 cooling.f90
 #
 # equations of state
 #
 SRCMESA= eos_mesa_microphysics.f90 eos_mesa.f90
-SRCEOS = eos_barotropic.f90 eos_stratified.f90 eos_piecewise.f90 ${SRCMESA} eos_shen.f90 eos_helmholtz.f90 eos_idealplusrad.f90 ionization.F90 eos_gasradrec.f90 eos.F90
+SRCEOS = eos_barotropic.f90 eos_stratified.f90 eos_piecewise.f90 ${SRCMESA} eos_shen.f90 eos_helmholtz.f90 eos_idealplusrad.f90 ionization.F90 eos_gasradrec.f90 eos.f90
 
 ifeq ($(HDF5), yes)
     SRCREADWRITE_DUMPS= utils_hdf5.f90 utils_dumpfiles_hdf5.f90 readwrite_dumps_common.F90 readwrite_dumps_fortran.F90 readwrite_dumps_hdf5.F90 readwrite_dumps.F90
@@ -535,14 +535,14 @@ SOURCES= physcon.f90 ${CONFIG} ${SRCKERNEL} io.F90 units.f90 \
          checkoptions.F90 viscosity.f90 damping.f90 options.f90 cons2primsolver.f90 radiation_utils.f90 cons2prim.f90 \
          centreofmass.f90 ${SRCPOT} checkconserved.f90 \
          utils_filenames.f90 utils_summary.F90 ${SRCCHEM} ${SRCDUST} \
-         mpi_memory.F90 mpi_derivs.F90 mpi_tree.F90 kdtree.F90 linklist_kdtree.F90 utils_healpix.f90 utils_raytracer.f90 \
-         partinject.F90 utils_inject.f90 dust_formation.F90 ptmass_radiation.F90 ptmass_heating.F90 \
+         mpi_memory.f90 mpi_derivs.F90 mpi_tree.F90 kdtree.F90 linklist_kdtree.F90 utils_healpix.f90 utils_raytracer.f90 \
+         partinject.F90 utils_inject.f90 dust_formation.f90 ptmass_radiation.f90 ptmass_heating.f90 \
          radiation_implicit.f90 ${SRCTURB} \
          ${SRCPHOTO} ${SRCINJECT_DEPS} ${SRCINJECT} \
          ${SRCKROME} memory.F90 ${SRCREADWRITE_DUMPS}  \
          quitdump.f90 ptmass.F90 \
-         readwrite_infile.F90 dens.F90 force.F90 utils_deriv.f90 deriv.F90 energies.F90 sort_particles.F90 \
-         utils_shuffleparticles.F90 evwrite.F90 step_leapfrog.F90 writeheader.F90 ${SRCAN} step_supertimestep.F90 \
+         readwrite_infile.F90 dens.F90 force.F90 utils_deriv.f90 deriv.F90 energies.F90 sort_particles.f90 \
+         utils_shuffleparticles.F90 evwrite.f90 step_leapfrog.F90 writeheader.F90 ${SRCAN} step_supertimestep.F90 \
          mf_write.f90 evolve.F90 utils_orbits.f90 utils_linalg.f90 \
          checksetup.F90 initial.F90
 
@@ -612,7 +612,7 @@ edit: checksetup
 SRCDUMP= physcon.f90 ${CONFIG} ${SRCKERNEL} utils_omp.F90 io.F90 units.f90 boundary.f90 boundary_dynamic.f90 mpi_utils.F90 \
          utils_timing.f90 utils_infiles.f90 dtype_kdtree.f90 utils_allocate.f90 part.F90 ${DOMAIN} \
          mpi_dens.F90 mpi_force.F90 \
-         mpi_balance.F90 mpi_memory.F90 mpi_derivs.F90 mpi_tree.F90 kdtree.F90 linklist_kdtree.F90 \
+         mpi_balance.F90 mpi_memory.f90 mpi_derivs.F90 mpi_tree.F90 kdtree.F90 linklist_kdtree.F90 \
          utils_dumpfiles.f90 utils_vectors.f90 utils_mathfunc.f90 \
          utils_datafiles.f90 utils_filenames.f90 utils_system.f90 utils_tables.f90 datafiles.f90 gitinfo.f90 \
          centreofmass.f90 \
@@ -677,14 +677,14 @@ cleansetup:
 # phantom test suite
 #
 ifeq ($(MPI),yes)
-    SRCTESTMPI = test_mpi.F90
+    SRCTESTMPI = test_mpi.f90
 else
     SRCTESTMPI =
 endif
 
 SRCTESTS=utils_testsuite.f90 ${TEST_FASTMATH} test_kernel.f90 \
-         test_dust.F90 test_growth.F90 test_smol.F90 \
-         test_nonidealmhd.F90 directsum.f90 test_gravity.F90 \
+         test_dust.F90 test_growth.f90 test_smol.F90 \
+         test_nonidealmhd.F90 directsum.f90 test_gravity.f90 \
          test_derivs.F90 test_cooling.f90 test_eos_stratified.f90 \
          test_eos.f90 test_externf.f90 test_rwdump.f90 \
          test_step.F90 test_indtstep.F90 set_disc.F90 test_setdisc.F90 \

build/Makefile_fastmath

@@ -29,7 +29,7 @@ endif
 
 ifeq ($(FASTMATH), yes)
     SRCFASTMATH=fastmath.o
-    TEST_FASTMATH=test_fastmath.F90
+    TEST_FASTMATH=test_fastmath.f90
     FPPFLAGS+=-DFINVSQRT
 else
     SRCFASTMATH=
@@ -38,7 +38,7 @@ endif
 
 fastmath.o: fastmath.f90
 	$(FC) $(FFLAGS) -o $@ -c $< || ${MAKE} fastmathlinkerr
-test_fastmath.o: test_fastmath.F90
+test_fastmath.o: test_fastmath.f90
 	$(FC) $(FFLAGS) -o $@ -c $< || ${MAKE} fastmathlinkerr
 getmathflags.o: getmathflags.f90
 	$(FC) $(FFLAGS) -o $@ -c $< || ${MAKE} fastmathlinkerr

build/Makefile_setups

@@ -732,7 +732,7 @@ ifeq ($(SETUP), star)
 #   import stellar model from 1D stellar evolution code
     SETUPFILE= setup_star.f90
     MODFILE= utils_binary.f90 set_binary.f90 moddump_binary.f90
-    ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.F90
+    ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.f90
     KNOWN_SETUP=yes
     MAXP=10000000
     GRAVITY=yes
@@ -745,7 +745,7 @@ ifeq ($(SETUP), grstar)
     IND_TIMESTEPS=yes
     SETUPFILE= setup_star.f90
     MODFILE= moddump_tidal.f90
-    ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.F90
+    ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.f90
     KNOWN_SETUP=yes
     MAXP=100000000
     GRAVITY=yes
@@ -756,7 +756,7 @@ ifeq ($(SETUP), dustystar)
     FPPFLAGS= -DDUST_NUCLEATION -DSINK_RADIATION -DSTAR
     SETUPFILE= setup_star.f90
     MODFILE= utils_binary.f90 set_binary.f90 moddump_binary.f90
-    ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.F90 dust_formation.F90
+    ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.f90 dust_formation.F90
     KNOWN_SETUP=yes
     MAXP=10000000
     GRAVITY=yes
@@ -766,7 +766,7 @@ ifeq ($(SETUP), radstar)
 #   setup a star as in the star setup but with radiation
     SETUPFILE= setup_star.f90
     MODFILE= utils_binary.f90 set_binary.f90 moddump_binary.f90
-    ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.F90
+    ANALYSIS= ${SRCNIMHD} utils_summary.o utils_omp.o ptmass.o energies.o analysis_common_envelope.f90
     KNOWN_SETUP=yes
     MAXP=10000000
     GRAVITY=yes
@@ -795,13 +795,13 @@ endif
 ifeq ($(SETUP), radwind)
 #   wind setup with dust nucleation
     FPPFLAGS= -DWIND -DSINK_RADIATION -DDUST_NUCLEATION
-    ANALYSIS= utils_getneighbours.f90 utils_raytracer_all.F90 analysis_raytracer.f90
+    ANALYSIS= utils_getneighbours.f90 utils_raytracer_all.f90 analysis_raytracer.f90
     WIND=yes
 endif
 
 ifeq ($(WIND), yes)
 #   wind particle injection setup
-    SETUPFILE= setup_wind.F90
+    SETUPFILE= setup_wind.f90
     SRCINJECT= utils_binary.f90 set_binary.f90 wind_equations.F90 wind.F90 inject_wind.f90
     KNOWN_SETUP=yes
 endif

src/main/eos.F90 → src/main/eos.f90

@@ -1374,6 +1374,7 @@ end subroutine read_headeropts_eos
 !+
 !-----------------------------------------------------------------------
 subroutine write_options_eos(iunit)
+ use dim,            only:use_krome
  use infile_utils,   only:write_inopt
  use eos_helmholtz,  only:eos_helmholtz_write_inopt
  use eos_barotropic, only:write_options_eos_barotropic
@@ -1383,9 +1384,11 @@ subroutine write_options_eos(iunit)
 
  write(iunit,"(/,a)") '# options controlling equation of state'
  call write_inopt(ieos,'ieos','eqn of state (1=isoth;2=adiab;3=locally iso;8=barotropic)',iunit)
-#ifndef KROME
- if (.not. eos_outputs_mu(ieos)) call write_inopt(gmw,'mu','mean molecular weight',iunit)
-#endif
+
+ if (.not.use_krome .or. .not.eos_outputs_mu(ieos)) then
+    call write_inopt(gmw,'mu','mean molecular weight',iunit)
+ endif
+
  select case(ieos)
  case(8)
     call write_options_eos_barotropic(iunit)

src/main/evwrite.F90 → src/main/evwrite.f90

@@ -247,7 +247,7 @@ subroutine init_evfile(iunit,evfile,open_file)
  endif
 
 end subroutine init_evfile
-!
+
 !----------------------------------------------------------------
 !+
 !  creates up to three lables per input value, and fills the required
@@ -300,7 +300,7 @@ subroutine fill_ev_tag(ev_fmt,itag,label,cmd,i,j)
        endif
     enddo
  enddo
- !
+
 end subroutine fill_ev_tag
 !----------------------------------------------------------------
 !+
@@ -350,11 +350,9 @@ end subroutine fill_ev_header
 subroutine write_evfile(t,dt)
  use energies,      only:compute_energies,ev_data_update
  use io,            only:id,master,ievfile
-#ifndef GR
  use timestep,      only:dtmax_user
  use options,       only:iexternalforce
- use extern_binary, only:accretedmass1,accretedmass2
-#endif
+ use externalforces,only:accretedmass1,accretedmass2
  real, intent(in)  :: t,dt
  integer           :: i,j
  real              :: ev_data_out(ielements)
@@ -364,14 +362,12 @@ subroutine write_evfile(t,dt)
 
  if (id==master) then
     !--fill in additional details that are not calculated in energies.f
-#ifndef GR
     ev_data(iev_sum,iev_dt)  = dt
     ev_data(iev_sum,iev_dtx) = dtmax_user
     if (iexternalforce==iext_binary) then
        ev_data(iev_sum,iev_maccsink(1)) = accretedmass1
        ev_data(iev_sum,iev_maccsink(2)) = accretedmass2
     endif
-#endif
     ! Fill in the data_out array
     j = 1
     do i = 1,iquantities
@@ -405,7 +401,6 @@ subroutine write_evfile(t,dt)
     call flush(ievfile)
  endif
 
- return
 end subroutine write_evfile
 !----------------------------------------------------------------
 !+
@@ -499,7 +494,6 @@ subroutine write_evlog(iprint)
  endif
  write(iprint,"(/)")
 
- return
 end subroutine write_evlog
-!----------------------------------------------------------------
+
 end module evwrite

src/main/externalforces.F90 → src/main/externalforces.f90

@@ -23,14 +23,14 @@ module externalforces
 !   extern_lensethirring, extern_prdrag, extern_spiral, extern_staticsine,
 !   infile_utils, io, lumin_nsdisc, part, units
 !
- use extern_binary,   only:accradius1,mass1
+ use extern_binary,   only:accradius1,mass1,accretedmass1,accretedmass2
  use extern_corotate, only:omega_corotate  ! so public from this module
  implicit none
 
  private
  public :: externalforce,externalforce_vdependent
  public :: accrete_particles,was_accreted
- public :: accradius1,omega_corotate
+ public :: accradius1,omega_corotate,accretedmass1,accretedmass2
  public :: write_options_externalforces,read_options_externalforces
  public :: initialise_externalforces,is_velocity_dependent
  public :: update_vdependent_extforce_leapfrog

src/main/externalforces_gr.F90 → src/main/externalforces_gr.f90

@@ -39,6 +39,8 @@ module externalforces
  public :: mass1  ! exported from metric module
  real, public :: accradius1 = 0.
  real, public :: accradius1_hard = 0.
+ real, public :: accretedmass1 = 0.
+ real, public :: accretedmass2 = 0.
 
  logical, public :: extract_iextern_from_hdr = .false.
 

src/main/step_leapfrog.F90

@@ -100,7 +100,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew)
                           iphase,iamtype,massoftype,maxphase,igas,idust,mhd,&
                           iamboundary,get_ntypes,npartoftypetot,&
                           dustfrac,dustevol,ddustevol,eos_vars,alphaind,nptmass,&
-                          dustprop,ddustprop,dustproppred,ndustsmall,pxyzu,dens,metrics,ics
+                          dustprop,ddustprop,dustproppred,pxyzu,dens,metrics,ics
  use options,        only:avdecayconst,alpha,ieos,alphamax
  use deriv,          only:derivs
  use timestep,       only:dterr,bignumber,tolv
@@ -114,7 +114,7 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew)
  use part,           only:gamma_chem
 #endif
  use timestep,       only:dtmax,dtmax_ifactor,dtdiff
- use timestep_ind,   only:get_dt,nbinmax,decrease_dtmax,ibinnow,dt_too_small
+ use timestep_ind,   only:get_dt,nbinmax,decrease_dtmax,dt_too_small
  use timestep_sts,   only:sts_get_dtau_next,use_sts,ibin_sts,sts_it_n
  use part,           only:ibin,ibin_old,twas,iactive,ibin_wake
 #ifdef GR
@@ -138,8 +138,8 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew)
  real,    intent(out)   :: dtnew
  integer            :: i,its,np,ntypes,itype,nwake,nvfloorp,nvfloorps,nvfloorc,ialphaloc
  real               :: timei,erri,errmax,v2i,errmaxmean
- real               :: vxi,vyi,vzi,eni,vxoldi,vyoldi,vzoldi,hdtsph,pmassi
- real               :: alphaloci,divvdti,source,tdecay1,hi,rhoi,ddenom,spsoundi
+ real               :: vxi,vyi,vzi,eni,hdtsph,pmassi
+ real               :: alphaloci,source,tdecay1,hi,rhoi,ddenom,spsoundi
  real               :: v2mean,hdti
  real(kind=4)       :: t1,t2,tcpu1,tcpu2
  real               :: pxi,pyi,pzi,p2i,p2mean
@@ -265,12 +265,12 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew)
 !$omp shared(pxyzu,ppred) &
 !$omp shared(Bevol,dBevol,Bpred,dtsph,massoftype,iphase) &
 !$omp shared(dustevol,ddustprop,dustprop,dustproppred,dustfrac,ddustevol,dustpred,use_dustfrac) &
-!$omp shared(alphaind,ieos,alphamax,ndustsmall,ialphaloc) &
+!$omp shared(alphaind,ieos,alphamax,ialphaloc) &
 !$omp shared(eos_vars,ufloor) &
 !$omp shared(twas,timei) &
 !$omp shared(rad,drad,radpred)&
 !$omp private(hi,rhoi,tdecay1,source,ddenom,hdti) &
-!$omp private(i,spsoundi,alphaloci,divvdti) &
+!$omp private(i,spsoundi,alphaloci) &
 !$omp firstprivate(pmassi,itype,avdecayconst,alpha) &
 !$omp reduction(+:nvfloorps)
  predict_sph: do i=1,npart
@@ -428,10 +428,10 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew)
 !$omp shared(xyzmh_ptmass,vxyz_ptmass,fxyz_ptmass,nptmass,massoftype) &
 !$omp shared(dtsph,ieos,ufloor) &
 !$omp shared(ibin,ibin_old,ibin_sts,twas,timei,use_sts,dtsph_next,ibin_wake,sts_it_n) &
-!$omp shared(ibin_dts,nbinmax,ibinnow) &
+!$omp shared(ibin_dts,nbinmax) &
 !$omp private(dti,hdti) &
 !$omp shared(rad,radpred,drad)&
-!$omp private(i,vxi,vyi,vzi,vxoldi,vyoldi,vzoldi) &
+!$omp private(i,vxi,vyi,vzi) &
 !$omp private(pxi,pyi,pzi,p2i) &
 !$omp private(erri,v2i,eni) &
 !$omp reduction(max:errmax) &
@@ -546,7 +546,6 @@ subroutine step(npart,nactive,t,dtsph,dtextforce,dtnew)
                 if (maxvxyzu >= 4) eni = vxyzu(4,i) + hdtsph*fxyzu(4,i)
 
                 erri = (vxi - vpred(1,i))**2 + (vyi - vpred(2,i))**2 + (vzi - vpred(3,i))**2
-                !if (erri > errmax) print*,id,' errmax = ',erri,' part ',i,vxi,vxoldi,vyi,vyoldi,vzi,vzoldi
                 errmax = max(errmax,erri)
 
                 v2i    = vxi*vxi + vyi*vyi + vzi*vzi
@@ -1029,15 +1028,15 @@ end subroutine step_extern_gr
 !+
 !----------------------------------------------------------------
 subroutine step_extern_sph(dt,npart,xyzh,vxyzu)
- use part, only:isdead_or_accreted,iboundary,iphase,iamtype
+ use part, only:isdead_or_accreted
  real,    intent(in)    :: dt
  integer, intent(in)    :: npart
  real,    intent(inout) :: xyzh(:,:)
  real,    intent(in)    :: vxyzu(:,:)
  integer :: i
 
  !$omp parallel do default(none) &
- !$omp shared(npart,xyzh,vxyzu,dt,iphase) &
+ !$omp shared(npart,xyzh,vxyzu,dt) &
  !$omp private(i)
  do i=1,npart
     if (.not.isdead_or_accreted(xyzh(4,i))) then
@@ -1097,11 +1096,11 @@ subroutine step_extern(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,fext,fxyzu,time,
  real,            intent(inout) :: xyzmh_ptmass(:,:),vxyz_ptmass(:,:),fxyz_ptmass(:,:)
  integer(kind=1), intent(in)    :: nbinmax
  integer(kind=1), intent(inout) :: ibin_wake(:)
- integer         :: i,itype,nsubsteps,ichem,naccreted,nfail,nfaili,merge_n
+ integer         :: i,itype,nsubsteps,naccreted,nfail,nfaili,merge_n
  integer         :: merge_ij(nptmass)
  integer(kind=1) :: ibin_wakei
  real            :: timei,hdt,fextx,fexty,fextz,fextxi,fextyi,fextzi,phii,pmassi
- real            :: dtphi2,dtphi2i,vxhalfi,vyhalfi,vzhalfi,fxi,fyi,fzi,deni
+ real            :: dtphi2,dtphi2i,vxhalfi,vyhalfi,vzhalfi,fxi,fyi,fzi
  real            :: dudtcool,fextv(3),poti,ui,rhoi
  real            :: dt,dtextforcenew,dtsinkgas,fonrmax,fonrmaxi
  real            :: dtf,accretedmass,t_end_step,dtextforce_min
@@ -1204,8 +1203,8 @@ subroutine step_extern(npart,ntypes,dtsph,dtextforce,xyzh,vxyzu,fext,fxyzu,time,
 #endif
     !$omp private(dphot,abundi,gmwvar) &
     !$omp private(ui,rhoi) &
-    !$omp private(i,ichem,dudtcool,fxi,fyi,fzi,phii) &
-    !$omp private(fextx,fexty,fextz,fextxi,fextyi,fextzi,poti,deni,fextv,accreted) &
+    !$omp private(i,dudtcool,fxi,fyi,fzi,phii) &
+    !$omp private(fextx,fexty,fextz,fextxi,fextyi,fextzi,poti,fextv,accreted) &
     !$omp private(fonrmaxi,dtphi2i,dtf) &
     !$omp private(vxhalfi,vyhalfi,vzhalfi) &
     !$omp firstprivate(pmassi,itype) &