Merge pull request #461 from danieljprice/fm17-master

same as #460

AUTHORS

@@ -17,18 +17,18 @@ Mark Hutchison <markahutch@gmail.com>
 Fitz Hu <fhuu0005@student.monash.edu>
 Megha Sharma <megha.sharma@monash.edu>
 Rebecca Nealon <rebecca.nealon@warwick.ac.uk>
+Elisabeth Borchert <elisabeth.borchert@monash.edu>
 Ward Homan <ward.homan@kuleuven.be>
 Christophe Pinte <christophe.pinte@univ-grenoble-alpes.fr>
-Elisabeth Borchert <elisabeth.borchert@monash.edu>
 Fangyi (Fitz) Hu <fhuu0005@student.monash.edu>
 Megha Sharma <msha0023@student.monash.edu>
 Terrence Tricco <ttricco@cita.utoronto.ca>
 Mats Esseldeurs <matsesseldeurs@yahoo.com>
 Stephane Michoulier <stephane.michoulier@gmail.com>
 Simone Ceppi <simone.ceppi@unimi.it>
 MatsEsseldeurs <matsesseldeurs@yahoo.com>
-Enrico Ragusa <enr.ragusa@gmail.com>
 Caitlyn Hardiman <caitlyn.hardiman1@monash.edu>
+Enrico Ragusa <enr.ragusa@gmail.com>
 fhu <fhuu0005@student.monash.edu>
 Sergei Biriukov <svbiriukov@gmail.com>
 Cristiano Longarini <cristiano.longarini@unimi.it>
@@ -37,12 +37,12 @@ Roberto Iaconi <robertoiaconi1@gmail.com>
 Hauke Worpel <hworpel@aip.de>
 Alison Young <alison.young@ed.ac.uk>
 Simone Ceppi <simo.ceppi@gmail.com>
-Stephen Neilson <36410751+s-neilson@users.noreply.github.com>
 Amena Faruqi <afaruqi.97@gmail.com>
+Stephen Neilson <36410751+s-neilson@users.noreply.github.com>
 Martina Toscani <mtoscani94@gmail.com>
 Benedetta Veronesi <benedetta.veronesi@unimi.it>
-Simon Glover <glover@uni-heidelberg.de>
 Sahl Rowther <sahl.rowther@leicester.ac.uk>
+Simon Glover <glover@uni-heidelberg.de>
 Thomas Reichardt <thomas.reichardt@students.mq.edu.au>
 Jean-François Gonzalez <Jean-Francois.Gonzalez@ens-lyon.fr>
 Christopher Russell <crussell@udel.edu>
@@ -52,24 +52,25 @@ Jolien Malfait <jolien.malfait@kuleuven.be>
 Phantom benchmark bot <ubuntu@phantom-benchmarks.novalocal>
 Kieran Hirsh <kieran.hirsh1@monash.edu>
 Nicole Rodrigues <nicole.rodrigues@monash.edu>
-David Trevascus <dtre10@student.monash.edu>
 Amena Faruqi <Amena.Faruqi@warwick.ac.uk>
+David Trevascus <dtre10@student.monash.edu>
+Chris Nixon <cjn@leicester.ac.uk>
 Megha Sharma <megha@meghas-air.home>
 Nicolas Cuello <cuellonicolas@gmail.com>
-Chris Nixon <cjn@leicester.ac.uk>
-Orsola De Marco <orsola.demarco@mq.edu.au>
-Joe Fisher <jwfis1@student.monash.edu>
+Benoit Commercon <benoit.commercon@gmail.com>
+Farzana Meru <f.meru@warwick.ac.uk>
 Giulia Ballabio <giulia.ballabio2@studenti.unimi.it>
-s-neilson <36410751+s-neilson@users.noreply.github.com>
-Megha Sharma <megha@dyn-49-127-43-125.its.monash.edu.au>
+Joe Fisher <jwfis1@student.monash.edu>
 Maxime Lombart <maxime.lombart@ens-lyon.fr>
+Megha Sharma <megha@dyn-49-127-43-125.its.monash.edu.au>
+Orsola De Marco <orsola.demarco@mq.edu.au>
 Terrence Tricco <tstricco@mun.ca>
-Benoit Commercon <benoit.commercon@gmail.com>
 Zachary Pellow <zpel1@student.monash.edu>
+s-neilson <36410751+s-neilson@users.noreply.github.com>
+Alison Young <ayoung@astro.ex.ac.uk>
 Cox, Samuel <sc676@leicester.ac.uk>
-mats esseldeurs <matse@naos.ster.kuleuven.be>
-Nicolás Cuello <cuello.nicolas@gmail.com>
 Jorge Cuadra <jcuadra@astro.puc.cl>
+Nicolás Cuello <cuello.nicolas@gmail.com>
 Steven Rieder <steven@rieder.nl>
 Stéven Toupin <steven.toupin@gmail.com>
-Alison Young <ayoung@astro.ex.ac.uk>
+mats esseldeurs <matse@naos.ster.kuleuven.be>

src/main/checksetup.F90

@@ -422,7 +422,7 @@ subroutine check_setup(nerror,nwarn,restart)
 !
  call get_centreofmass(xcom,vcom,npart,xyzh,vxyzu,nptmass,xyzmh_ptmass,vxyz_ptmass)
 
- if (.not.h2chemistry .and. maxvxyzu >= 4 .and. icooling == 3 .and. iexternalforce/=iext_corotate) then
+ if (.not.h2chemistry .and. maxvxyzu >= 4 .and. icooling == 3 .and. iexternalforce/=iext_corotate .and. nptmass==0) then
     if (dot_product(xcom,xcom) >  1.e-2) then
        print*,'ERROR: Gammie (2001) cooling (icooling=3) assumes Omega = 1./r^1.5'
        print*,'                but the centre of mass is not at the origin!'

src/main/cooling_gammie.f90

@@ -17,7 +17,7 @@ module cooling_gammie
 ! :Runtime parameters:
 !   - beta_cool : *beta factor in Gammie (2001) cooling*
 !
-! :Dependencies: infile_utils, io
+! :Dependencies: infile_utils, io, part
 !
  implicit none
  real, private :: beta_cool  = 3.
@@ -29,13 +29,18 @@ module cooling_gammie
 !+
 !-----------------------------------------------------------------------
 subroutine cooling_Gammie_explicit(xi,yi,zi,ui,dudti)
-
+ use part, only:xyzmh_ptmass, nptmass
  real, intent(in)    :: ui,xi,yi,zi
  real, intent(inout) :: dudti
 
  real :: omegai,r2,tcool1
 
- r2     = xi*xi + yi*yi + zi*zi
+ if (nptmass > 0) then
+    r2     = (xi-xyzmh_ptmass(1,1))**2 + (yi-xyzmh_ptmass(2,1))**2 + (zi-xyzmh_ptmass(3,1))**2
+ else
+    r2     = xi*xi + yi*yi + zi*zi
+ endif
+
  Omegai = r2**(-0.75)
  tcool1 = Omegai/beta_cool
  dudti  = dudti - ui*tcool1

src/main/radiation_implicit.f90

@@ -18,8 +18,9 @@ module radiation_implicit
 !
 ! :Runtime parameters: None
 !
-! :Dependencies: boundary, dim, eos, io, kdtree, kernel, linklist, options,
-!   part, physcon, quartic, radiation_utils, units
+! :Dependencies: boundary, derivutils, dim, eos, implicit, io, kdtree,
+!   kernel, linklist, options, part, physcon, quartic, radiation_utils,
+!   timing, units
 !
  use part,            only:ikappa,ilambda,iedd,idkappa,iradxi,icv,ifluxx,ifluxy,ifluxz,igas,rhoh,massoftype,imu
  use eos,             only:iopacity_type

src/main/utils_implicit.f90

@@ -7,14 +7,14 @@
 module implicit
 !
 ! Utility routines for implicit radiative diffusion
-!   
+!
 ! :References:
 !
 ! :Owner: Mike Lau
 !
 ! :Runtime parameters: None
 !
-! :Dependencies:
+! :Dependencies: io, physcon
 !
  implicit none
  integer              :: ncompact,ncompactlocal,icompactmax,nneigh_average

src/utils/combinedustdumps.f90

@@ -28,7 +28,7 @@ program combinedustdumps
  use part,            only:xyzh,vxyzu,npart,hfact,iphase,npartoftype,massoftype,&
                            igas,idust,ndusttypes,ndustsmall,ndustlarge,set_particle_type,&
                            grainsize,graindens,iamtype,isdead_or_accreted
- use readwrite_dumps, only:read_dump,write_fulldump
+ use readwrite_dumps, only:read_dump,write_fulldump,init_readwrite_dumps
  use units,           only:set_units,select_unit,umass,udist,utime
  use memory,          only:allocate_memory
  use checksetup,      only:check_setup
@@ -38,10 +38,12 @@ program combinedustdumps
  real, allocatable :: xyzh_tmp(:,:,:),vxyzu_tmp(:,:,:),massofdust_tmp(:)
  real, allocatable :: grainsize_tmp(:),graindens_tmp(:)
  integer, allocatable :: npartofdust_tmp(:)
- integer :: i,j,counter,ipart,itype,ierr,nargs,idust_tmp,ninpdumps
- integer :: nwarn,nerror
+ integer :: i,j,ipart,itype,ierr,nargs,idust_tmp,ninpdumps
+ integer(kind=8) :: counter
+ integer :: nwarn,nerror,ndust
  real    :: time
  real(kind=8) :: utime_tmp,udist_tmp,umass_tmp
+ logical :: first_gas_only
 
  call set_io_unit_numbers
  iprint = 6
@@ -70,6 +72,8 @@ program combinedustdumps
  ! read first dumpfile HEADER ONLY: check idust, check MAXP
  ! we assume all dumps are from the same phantom version
  !
+
+ call init_readwrite_dumps()
  counter = 0
  idust_tmp = idust ! new dumps, location of first dust particle type
  do i=1,ninpdumps
@@ -81,6 +85,11 @@ program combinedustdumps
     endif
     counter = counter + npartoftype(idust_tmp)
  enddo
+ !
+ ! save the number of dust particles for later
+ !
+ ndust = npartoftype(idust_tmp)
+
  !
  !--sanity check array sizes
  !
@@ -91,17 +100,21 @@ program combinedustdumps
  !
  ! allocate memory
  !
- call allocate_memory(int(counter,kind=8))
+ call allocate_memory(counter)
  !
  ! read gas particles from first file
  !
+
  call read_dump(trim(indumpfiles(1)),time,hfact,idisk1,iprint,0,1,ierr)
+
  !
  ! allocate temporary arrays
  !
- allocate (xyzh_tmp(ninpdumps,4,npartoftype(idust_tmp)),stat=ierr)
+
+ allocate (xyzh_tmp(ninpdumps,4,ndust),stat=ierr)
+ print*,shape(xyzh_tmp),npartoftype(idust_tmp),idust_tmp
  if (ierr /= 0) stop 'error allocating memory to store positions'
- allocate (vxyzu_tmp(ninpdumps,maxvxyzu,npartoftype(idust_tmp)),stat=ierr)
+ allocate (vxyzu_tmp(ninpdumps,maxvxyzu,ndust),stat=ierr)
  if (ierr /= 0) stop 'error allocating memory to store velocities'
  allocate (npartofdust_tmp(ninpdumps),stat=ierr)
  if (ierr /= 0) stop 'error allocating memory to store number of dust particles'
@@ -115,6 +128,7 @@ program combinedustdumps
  !
  !--read dumps and get dust particle information
  !
+ first_gas_only = .false.
  do i=1,ninpdumps
     call read_dump(trim(indumpfiles(i)),time,hfact,idisk1,iprint,0,1,ierr)
     if (ierr /= 0) stop 'error reading dumpfile'
@@ -123,6 +137,7 @@ program combinedustdumps
     grainsize_tmp(i) = grainsize(1)
     graindens_tmp(i) = graindens(1)
     counter = 0
+    if (i == 1 .and. npartoftype(idust_tmp) == 0) first_gas_only = .true.
     do j=1,maxp
        if (iphase(j)==idust_tmp) then
           counter = counter + 1
@@ -156,6 +171,7 @@ program combinedustdumps
  !
  do i=2,ninpdumps
     itype = idust + i - 1
+    if (first_gas_only) itype = idust + i - 2
     npartoftype(itype) = npartofdust_tmp(i)
     massoftype(itype) = massofdust_tmp(i)
     grainsize(i) = grainsize_tmp(i)
@@ -176,6 +192,7 @@ program combinedustdumps
  !--dust properties
  !
  ndusttypes = ninpdumps
+ if (first_gas_only) ndusttypes = ninpdumps - 1
  ndustlarge = ndusttypes
  ndustsmall = 0