Merge pull request #36 from jameswurster/master

formatting changes

src/main/force.F90

@@ -163,7 +163,7 @@ subroutine force(icall,npart,xyzh,vxyzu,fxyzu,divcurlv,divcurlB,Bevol,dBevol,&
  use io,           only:iprint,fatal,iverbose,id,master,real4,warning,error,nprocs
  use linklist,     only:ncells,get_neighbour_list,get_hmaxcell,get_cell_location
  use options,      only:iresistive_heating
- use part,         only:rhoh,dhdrho,rhoanddhdrho,alphaind,nabundances,iactive,gradh,&
+ use part,         only:rhoh,dhdrho,rhoanddhdrho,alphaind,iactive,gradh,&
                         hrho,iphase,igas,maxgradh,dvdx,eta_nimhd,deltav,poten,iamtype
  use timestep,     only:dtcourant,dtforce,dtrad,bignumber,dtdiff
  use io_summary,   only:summary_variable, &
@@ -2516,7 +2516,7 @@ subroutine finish_cell_and_store_results(icall,cell,fxyzu,xyzh,vxyzu,poten,dt,dv
  use dim,            only:mhd,mhd_nonideal,lightcurve,use_dust,maxdvdx,maxBevol,use_dustgrowth,gr,use_krome
  use eos,            only:use_entropy,gamma,ieos
  use options,        only:alpha,icooling,ipdv_heating,ishock_heating,psidecayfac,overcleanfac,hdivbbmax_max,use_dustfrac,damp
- use part,           only:h2chemistry,rhoanddhdrho,abundance,iboundary,igas,maxphase,maxvxyzu,nabundances, &
+ use part,           only:h2chemistry,rhoanddhdrho,iboundary,igas,maxphase,maxvxyzu, &
                           massoftype,get_partinfo,tstop,strain_from_dvdx,ithick
 #ifdef IND_TIMESTEPS
  use part,           only:ibin

src/main/initial.F90

@@ -651,16 +651,18 @@ subroutine startrun(infile,logfile,evfile,dumpfile)
 !
 !--Print box sizes and masses
 !
- if (get_conserv > 0.0) then
-    write(iprint,'(1x,a)') 'Initial mass and box size (in code units):'
- else
-    write(iprint,'(1x,a)') 'Mass and box size (in code units) of this restart:'
+ if (id==master) then
+    if (get_conserv > 0.0) then
+       write(iprint,'(1x,a)') 'Initial mass and box size (in code units):'
+    else
+       write(iprint,'(1x,a)') 'Mass and box size (in code units) of this restart:'
+    endif
+    write(iprint,'(2x,a,es18.6)') 'Total mass: ', mtot
+    write(iprint,'(2x,a,es18.6)') 'x-Box size: ', dx
+    write(iprint,'(2x,a,es18.6)') 'y-Box size: ', dy
+    write(iprint,'(2x,a,es18.6)') 'z-Box size: ', dz
+    write(iprint,'(a)') ' '
  endif
- write(iprint,'(2x,a,es18.6)') 'Total mass: ', mtot
- write(iprint,'(2x,a,es18.6)') 'x-Box size: ', dx
- write(iprint,'(2x,a,es18.6)') 'y-Box size: ', dy
- write(iprint,'(2x,a,es18.6)') 'z-Box size: ', dz
- write(iprint,'(a)') ' '
 !
 !--Set initial values for continual verification of conservation laws
 !  get_conserve=0.5: update centre of mass only; get_conserve=1: update all; get_conserve=-1: update none
@@ -670,23 +672,28 @@ subroutine startrun(infile,logfile,evfile,dumpfile)
     angtot_in = angtot
     totmom_in = totmom
     mdust_in  = mdust
-    write(iprint,'(1x,a)') 'Setting initial values (in code units) to verify conservation laws:'
+    if (id==master) write(iprint,'(1x,a)') 'Setting initial values (in code units) to verify conservation laws:'
  else
-    write(iprint,'(1x,a)') 'Reading initial values (in code units) to verify conservation laws from previous run:'
+    if (id==master) then
+       write(iprint,'(1x,a)') 'Reading initial values (in code units) to verify conservation laws from previous run:'
+    endif
  endif
- write(iprint,'(2x,a,es18.6)') 'Initial total energy:     ', etot_in
- write(iprint,'(2x,a,es18.6)') 'Initial angular momentum: ', angtot_in
- write(iprint,'(2x,a,es18.6)') 'Initial linear momentum:  ', totmom_in
- if (use_dust) then
-    dust_label = 'dust'
-    call make_tags_unique(ndusttypes,dust_label)
-    do i=1,ndusttypes
-       if (mdust_in(i) > 0.) write(iprint,'(2x,a,es18.6)') 'Initial '//trim(dust_label(i))//' mass:     ',mdust_in(i)
-    enddo
-    write(iprint,'(2x,a,es18.6)') 'Initial total dust mass:  ', sum(mdust_in(:))
+ if (id==master) then
+    write(iprint,'(2x,a,es18.6)') 'Initial total energy:     ', etot_in
+    write(iprint,'(2x,a,es18.6)') 'Initial angular momentum: ', angtot_in
+    write(iprint,'(2x,a,es18.6)') 'Initial linear momentum:  ', totmom_in
+    if (use_dust) then
+       dust_label = 'dust'
+       call make_tags_unique(ndusttypes,dust_label)
+       do i=1,ndusttypes
+          if (mdust_in(i) > 0.) write(iprint,'(2x,a,es18.6)') 'Initial '//trim(dust_label(i))//' mass:     ',mdust_in(i)
+       enddo
+       write(iprint,'(2x,a,es18.6)') 'Initial total dust mass:  ', sum(mdust_in(:))
+    endif
  endif
 !
 !--Print warnings of units if values are not reasonable
+!
  tolu = 1.0d2
  toll = 1.0d-2
  if (get_conserv > 0.0) then

src/main/readwrite_dumps.F90

@@ -879,7 +879,7 @@ subroutine read_dump(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ierr,heade
 #ifdef INJECT_PARTICLES
        call allocate_memory(maxp_hard)
 #else
-       call allocate_memory(int(nparttot / nprocs) + 1)
+       call allocate_memory(int( min(nprocs,2)*nparttot / nprocs))
 #endif
     endif
 !
@@ -1071,7 +1071,7 @@ subroutine read_smalldump(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ierr,
 #ifdef INJECT_PARTICLES
  call allocate_memory(maxp_hard)
 #else
- call allocate_memory(int(nparttot / nprocs) + 1)
+ call allocate_memory(int( min(nprocs,2)*nparttot / nprocs))
 #endif
 
 !
@@ -1991,16 +1991,16 @@ subroutine unfill_rheader(hdr,phantomdump,ntypesinfile,nptmass,&
  call extract('rhozero',rhozero,hdr,ierr)
  call extract('RK2',rk2,hdr,ierr)
  polyk = 2./3.*rk2
- if (maxvxyzu >= 4) then
-    if (use_krome) then
-       if (id==master) write(iprint,*) 'KROME eos: initial gamma = 1.666667'
+ if (id==master) then
+    if (maxvxyzu >= 4) then
+       if (use_krome) then
+          write(iprint,*) 'KROME eos: initial gamma = 1.666667'
+       else
+          write(iprint,*) 'adiabatic eos: gamma = ',gamma
+       endif
     else
-       if (id==master) write(iprint,*) 'adiabatic eos: gamma = ',gamma
-    endif
- else
-    write(iprint,*) 'setting isothermal sound speed^2 (polyk) = ',polyk,' gamma = ',gamma
-    if (polyk <= tiny(polyk)) then
-       write(iprint,*) 'WARNING! sound speed zero in dump!, polyk = ',polyk
+       write(iprint,*) 'setting isothermal sound speed^2 (polyk) = ',polyk,' gamma = ',gamma
+       if (polyk <= tiny(polyk)) write(iprint,*) 'WARNING! sound speed zero in dump!, polyk = ',polyk
     endif
  endif
 
@@ -2085,10 +2085,12 @@ subroutine unfill_rheader(hdr,phantomdump,ntypesinfile,nptmass,&
     call extract('Bextx',Bextx,hdr,ierrs(1))
     call extract('Bexty',Bexty,hdr,ierrs(2))
     call extract('Bextz',Bextz,hdr,ierrs(3))
-    if (any(ierrs(1:3) /= 0)) then
-       write(*,*) 'ERROR reading external field (setting to zero)'
-    else
-       write(*,*) 'External field found, Bext = ',Bextx,Bexty,Bextz
+    if (id==master) then
+       if (any(ierrs(1:3) /= 0)) then
+          write(*,*) 'ERROR reading external field (setting to zero)'
+       else
+          write(*,*) 'External field found, Bext = ',Bextx,Bexty,Bextz
+       endif
     endif
  endif
 

src/main/readwrite_dumps_hdf5.F90

@@ -595,7 +595,7 @@ subroutine read_dump(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ierr,heade
 #ifdef INJECT_PARTICLES
  call allocate_memory(maxp_hard)
 #else
- call allocate_memory(int(npart / nprocs) + 1)
+ call allocate_memory(int( min(nprocs,2)*nparttot / nprocs))
 #endif
 
  if (periodic) then
@@ -796,7 +796,7 @@ subroutine read_smalldump(dumpfile,tfile,hfactfile,idisk1,iprint,id,nprocs,ierr,
 #ifdef INJECT_PARTICLES
  call allocate_memory(maxp_hard)
 #else
- call allocate_memory(int(npart / nprocs) + 1)
+ call allocate_memory(int( min(nprocs,2)*nparttot / nprocs))
 #endif
 
  if (periodic) then

src/setup/setup_dustsettle.f90

@@ -51,7 +51,7 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  use externalforces, only:Rdisc,iext_discgravity
  use options,        only:iexternalforce,use_dustfrac
  use timestep,       only:dtmax,tmax
- use units,          only:set_units,udist,umass
+ use units,          only:set_units,udist,umass,utime
  use dust,           only:init_drag,idrag,get_ts
  use set_dust,       only:set_dustfrac,set_dustbinfrac
  use table_utils,    only:logspace
@@ -183,8 +183,8 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  if (maxvxyzu >= 4) then
     gamma = 5./3.
  else
-    gamma  = 1.
-    polyk  = cs**2
+    gamma = 1.
+    polyk = cs**2
  endif
 !
 !--get stopping time information
@@ -246,11 +246,21 @@ subroutine setpart(id,npart,npartoftype,xyzh,massoftype,vxyzu,polyk,gamma,hfact,
  massoftype(itype)  = totmass/npartoftype(itype)*(1. + dtg)
 
  if (id==master) then
-    print*,' npart                 = ',npart,npart_total
+    print*,' npart,npart_total     = ',npart,npart_total
     print*,' particle mass         = ',massoftype(itype),'code units'
     print*,' particle mass         = ',massoftype(itype)*umass,'g'
     print*,' total mass            = ',npart*massoftype(itype)*umass,'g'
     print*,' mid-plane gas density = ',rhozero*umass/udist**3,'g/cm^3'
+    print*,' Rdisc                 = ',Rdisc*udist/au,'au'
+    print*,' H0                    = ',H0*udist/au,'au'
+    print*,' H/R                   = ',HonR
+    print*,' cs                    = ',cs*udist/utime,'cm/s'
+ endif
+ if (HonR > 0.1) then
+    print*, ' '
+    print*, 'WARNING! This disc is hot, and *may* blow apart rather than settle.  A smaller ratio of H/R may be required.'
+    print*, '         The default value of 0.05 yields stable results'
+    print*, ' '
  endif
 
 end subroutine setpart