Merge pull request #535 from danieljprice/chinchen

Chinese coin problem / bug fixes and tests added for sink particles in external potentials

AUTHORS

@@ -24,48 +24,49 @@ Christophe Pinte <christophe.pinte@univ-grenoble-alpes.fr>
 Terrence Tricco <tstricco@mun.ca>
 Stephane Michoulier <stephane.michoulier@gmail.com>
 Simone Ceppi <simo.ceppi@gmail.com>
+Yrisch <yannbernard13@gmail.com>
 Spencer Magnall <spencermagnall@gmail.com>
-Caitlyn Hardiman <caitlyn.hardiman1@monash.edu>
 Enrico Ragusa <enr.ragusa@gmail.com>
+Caitlyn Hardiman <caitlyn.hardiman1@monash.edu>
 Sergei Biriukov <svbiriukov@gmail.com>
 Cristiano Longarini <cristiano.longarini@unimi.it>
 Giovanni Dipierro <giovanni.dipierro@leicester.ac.uk>
 Roberto Iaconi <robertoiaconi1@gmail.com>
-Hauke Worpel <hworpel@aip.de>
 Amena Faruqi <Amena.Faruqi@warwick.ac.uk>
+Hauke Worpel <hworpel@aip.de>
 Alison Young <alison.young@ed.ac.uk>
 Stephen Neilson <36410751+s-neilson@users.noreply.github.com>
 Martina Toscani <mtoscani94@gmail.com>
 Benedetta Veronesi <benedetta.veronesi@unimi.it>
-Sahl Rowther <sahl.rowther@leicester.ac.uk>
 Simon Glover <glover@uni-heidelberg.de>
+Sahl Rowther <sahl.rowther@leicester.ac.uk>
 Thomas Reichardt <thomas.reichardt@students.mq.edu.au>
 Jean-François Gonzalez <Jean-Francois.Gonzalez@ens-lyon.fr>
 Christopher Russell <crussell@udel.edu>
-Alessia Franchini <alessia.franchini@unlv.edu>
 Alex Pettitt <alex@astro1.sci.hokudai.ac.jp>
+Alessia Franchini <alessia.franchini@unlv.edu>
 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>
-Nicolás Cuello <cuello.nicolas@gmail.com>
+Kieran Hirsh <kieran.hirsh1@monash.edu>
 Farzana Meru <f.meru@warwick.ac.uk>
+Nicolás Cuello <cuello.nicolas@gmail.com>
+David Trevascus <dtre10@student.monash.edu>
 Mike Lau <mike.lau@h-its.org>
-Chris Nixon <cjn@leicester.ac.uk>
 Miguel Gonzalez-Bolivar <miguelgb@astro.unam.mx>
+Chris Nixon <cjn@leicester.ac.uk>
 Orsola De Marco <orsola.demarco@mq.edu.au>
 Maxime Lombart <maxime.lombart@ens-lyon.fr>
 Joe Fisher <jwfis1@student.monash.edu>
-Zachary Pellow <zpel1@student.monash.edu>
-Benoit Commercon <benoit.commercon@gmail.com>
 Giulia Ballabio <giulia.ballabio2@studenti.unimi.it>
+Benoit Commercon <benoit.commercon@gmail.com>
+Zachary Pellow <zpel1@student.monash.edu>
 s-neilson <36410751+s-neilson@users.noreply.github.com>
-MICHOULIER Stephane <smichoulier@sidus11>
 Steven Rieder <steven@rieder.nl>
-Jeremy Smallwood <drjeremysmallwood@gmail.com>
-Cox, Samuel <sc676@leicester.ac.uk>
-Jorge Cuadra <jcuadra@astro.puc.cl>
-Stéven Toupin <steven.toupin@gmail.com>
 Taj Jankovič <taj.jankovic@gmail.com>
 Chunliang Mu <Mclliang7326@outlook.com>
+MICHOULIER Stephane <smichoulier@sidus11>
+Jorge Cuadra <jcuadra@astro.puc.cl>
+Cox, Samuel <sc676@leicester.ac.uk>
+Jeremy Smallwood <drjeremysmallwood@gmail.com>
+Stéven Toupin <steven.toupin@gmail.com>

scripts/bots.sh

@@ -207,6 +207,8 @@ for edittype in $bots_to_run; do
                'shout' )
                  sed -e 's/SQRT(/sqrt(/g' \
                      -e 's/NINT(/nint(/g' \
+                     -e 's/REAL(/real(/g' \
+                     -e 's/DBLE(/dble(/g' \
                      -e 's/ STOP/ stop/g' \
                      -e 's/ATAN/atan/g' \
                      -e 's/ACOS(/acos(/g' \
@@ -216,6 +218,17 @@ for edittype in $bots_to_run; do
                      -e 's/EXP(/exp(/g' \
                      -e 's/LOG(/log(/g' \
                      -e 's/READ(/read(/g' \
+                     -e 's/OPEN(/open(/g' \
+                     -e 's/OPEN (/open(/g' \
+                     -e 's/CLOSE(/close(/g' \
+                     -e 's/CLOSE (/close(/g' \
+                     -e 's/INDEX(/index(/g' \
+                     -e 's/ANY(/any(/g' \
+                     -e 's/, STATUS=/,status=/g' \
+                     -e 's/,STATUS=/,status=/g' \
+                     -e 's/, FORM=/,form=/g' \
+                     -e 's/,FORM=/,form=/g' \
+                     -e 's/TRIM(/trim(/g' \
                      -e 's/IF (/if (/g' \
                      -e 's/) THEN/) then/g' \
                      -e 's/ENDDO/enddo/g' \
@@ -267,6 +280,22 @@ for edittype in $bots_to_run; do
                  sed -e 's/end if/endif/g' \
                      -e 's/end do/enddo/g' \
                      -e 's/else if/elseif/g' \
+                     -e 's/open (/open(/g' \
+                     -e 's/, file = /,file=/g' \
+                     -e 's/, file=/,file=/g' \
+                     -e 's/, status = /,status=/g' \
+                     -e 's/, status=/,status=/g' \
+                     -e 's/, iostat = /,iostat=/g' \
+                     -e 's/, iostat=/,iostat=/g' \
+                     -e 's/, access = /,access=/g' \
+                     -e 's/, access=/,access=/g' \
+                     -e 's/, form = /,form=/g' \
+                     -e 's/, form=/,form=/g' \
+                     -e 's/, action = /,action=/g' \
+                     -e 's/, action=/,action=/g' \
+                     -e 's/, iomsg = /,iomsg=/g' \
+                     -e 's/, iomsg=/,iomsg=/g' \
+                     -e 's/(unit =/(unit=/g' \
                      -e 's/if(/if (/g' \
                      -e 's/)then/) then/g' $file > $out;;
                'header' )

src/main/checksetup.f90

@@ -16,7 +16,7 @@ module checksetup
 !
 ! :Dependencies: boundary, boundary_dyn, centreofmass, dim, dust, eos,
 !   externalforces, io, metric_tools, nicil, options, part, physcon,
-!   ptmass_radiation, sortutils, timestep, units, utils_gr
+!   ptmass, ptmass_radiation, sortutils, timestep, units, utils_gr
 !
  implicit none
  public :: check_setup

src/main/cooling_molecular.f90

@@ -210,21 +210,21 @@ subroutine loadCoolingTable(data_array)
 
  iunit = 1
  filename = find_phantom_datafile('radcool_all.dat','cooling')
- OPEN(unit=iunit, file=trim(filename), STATUS="OLD", ACTION="read", &
+ open(unit=iunit,file=trim(filename),status="OLD", ACTION="read", &
             iostat=istat, IOMSG=imsg)
 
  ! Begin loading in data
  openif: if (istat == 0) then
     !!! Skip header
     rewind(unit=iunit)
     do o = 1, headerLines
-       read(iunit, *, iostat=istat, IOMSG = imsg)
+       read(iunit, *,iostat=istat, IOMSG = imsg)
     enddo
 
     ! Read data
     skipheaderif: if ((istat == 0)) then
        readdo: do
-          read(iunit, *, iostat=istat) i, j, k, T, n_H, N_coolant, lambda_CO, lambda_H2O, lambda_HCN
+          read(iunit, *,iostat=istat) i, j, k, T, n_H, N_coolant, lambda_CO, lambda_H2O, lambda_HCN
           if (istat /= 0) exit
           data_array(i, j, k, :) = [T, n_H, N_coolant, lambda_CO, lambda_H2O, lambda_HCN]
 
@@ -275,20 +275,20 @@ subroutine loadCDTable(data_array)
 
  iunit = 1
  filename = find_phantom_datafile('table_cd.dat','cooling')
- open(unit=iunit, file=filename, STATUS="OLD", iostat=istat, IOMSG=imsg)
+ open(unit=iunit,file=filename,status="OLD",iostat=istat, IOMSG=imsg)
 
  ! Begin loading in data
  openif: if (istat == 0) then
     !!! Skip header
     rewind(unit=iunit)
     do o = 1, headerLines
-       read(iunit, *, iostat=istat, IOMSG = imsg)
+       read(iunit, *,iostat=istat, IOMSG = imsg)
     enddo
 
     !!! Read data
     skipheaderif: if ((istat == 0)) then
        readdo: do
-          read(iunit, *, iostat=istat) i, j, k, l, r_part, widthLine, m_exp, r_sep, N_H
+          read(iunit, *,iostat=istat) i, j, k, l, r_part, widthLine, m_exp, r_sep, N_H
           if (istat /= 0) exit
           data_array(i, j, k, l, :) = [r_part, widthLine, m_exp, r_sep, N_H]
 

src/main/eos_mesa_microphysics.f90

@@ -69,7 +69,7 @@ subroutine get_opacity_constants_mesa
  opacs_file = find_phantom_datafile(filename,'eos/mesa')
 
  ! Read the constants from the header of the opacity file
- open(newunit=fnum, file=trim(opacs_file), status='old', action='read', form='unformatted')
+ open(newunit=fnum,file=trim(opacs_file),status='old',action='read',form='unformatted')
  read(fnum) mesa_opacs_nz,mesa_opacs_nx,mesa_opacs_nr,mesa_opacs_nt
  close(fnum)
 
@@ -102,7 +102,7 @@ subroutine read_opacity_mesa(x,z)
  filename = trim(mesa_opacs_dir)//'opacs'//trim(mesa_opacs_suffix)//'.bindata'
 ! filename = trim(mesa_opacs_dir)//'/'//'opacs'//trim(mesa_opacs_suffix)//'.bindata'
  opacs_file = find_phantom_datafile(filename,'eos/mesa')
- open(unit=fnum, file=trim(opacs_file), status='old', action='read', form='unformatted')
+ open(unit=fnum,file=trim(opacs_file),status='old',action='read',form='unformatted')
  read(fnum) mesa_opacs_nz,mesa_opacs_nx,mesa_opacs_nr,mesa_opacs_nt
 
  ! Read in the size of the table and the data
@@ -308,7 +308,7 @@ subroutine get_eos_constants_mesa(ierr)
  filename = find_phantom_datafile(filename,'eos/mesa')
 
  ! Read constants from the header of first EoS tables
- open(unit=fnum, file=trim(filename), status='old', action='read', form='unformatted',iostat=ierr)
+ open(unit=fnum,file=trim(filename),status='old',action='read',form='unformatted',iostat=ierr)
  if (ierr /= 0) return
  read(fnum) mesa_eos_ne, mesa_eos_nv, mesa_eos_nvar2
  close(fnum)
@@ -364,7 +364,7 @@ subroutine read_eos_mesa(x,z,ierr)
        ! Read in the size of the tables and the data
        ! i and j hold the Z and X values respectively
        ! k, l and m hold the values of V, Eint and the data respectively
-       open(unit=fnum, file=trim(filename), status='old', action='read', form='unformatted')
+       open(unit=fnum,file=trim(filename),status='old',action='read',form='unformatted')
        read(fnum) mesa_eos_ne, mesa_eos_nv, mesa_eos_nvar2
        read(fnum)(mesa_eos_logVs(k),k=1,mesa_eos_nv)
        read(fnum)(mesa_eos_logEs(l),l=1,mesa_eos_ne)

src/main/extern_binary.f90

@@ -19,6 +19,7 @@ module extern_binary
 !   - accradius2 : *accretion radius of secondary (if iexternalforce=binary)*
 !   - eps_soft1  : *Plummer softening of primary*
 !   - eps_soft2  : *Plummer softening of secondary*
+!   - mass1      : *m1 of central binary system (if iexternalforce=binary)*
 !   - mass2      : *m2 of central binary system (if iexternalforce=binary)*
 !   - ramp       : *ramp up mass of secondary over first 5 orbits?*
 !
@@ -234,6 +235,7 @@ subroutine write_options_externbinary(iunit)
  use infile_utils,   only:write_inopt
  integer, intent(in) :: iunit
 
+ call write_inopt(mass1,'mass1','m1 of central binary system (if iexternalforce=binary)',iunit)
  call write_inopt(mass2,'mass2','m2 of central binary system (if iexternalforce=binary)',iunit)
  call write_inopt(accradius1,'accradius1','accretion radius of primary',iunit)
  call write_inopt(accradius2,'accradius2','accretion radius of secondary (if iexternalforce=binary)',iunit)
@@ -259,6 +261,12 @@ subroutine read_options_externbinary(name,valstring,imatch,igotall,ierr)
  imatch  = .true.
  igotall = .false.
  select case(trim(name))
+ case('mass1')
+    read(valstring,*,iostat=ierr) mass1
+    ngot = ngot + 1
+    if (mass1 < 0.) then
+       call fatal(where,'invalid setting for m1 (<0)')
+    endif
  case('mass2')
     read(valstring,*,iostat=ierr) mass2
     ngot = ngot + 1

src/main/extern_densprofile.f90

@@ -137,13 +137,13 @@ subroutine read_rhotab(filename, rsize, rtab, rhotab, nread, polyk, gamma, rhoc,
  endif
 
  ! First line: # K gamma rhoc
- read(iunit, *, iostat=ierr) hash,polyk, gamma, rhoc
+ read(iunit, *,iostat=ierr) hash,polyk, gamma, rhoc
  if (ierr /= 0) then
     call error('extern_densityprofile','Error reading first line of header from '//trim(filename))
     return
  endif
  ! Second line: # nentries  (number of r density entries in file)
- read(iunit,*, iostat=ierr) hash,nread
+ read(iunit,*,iostat=ierr) hash,nread
  if (ierr /= 0) then
     call error('extern_densityprofile','Error reading second line of header from '//trim(filename))
     return
@@ -155,7 +155,7 @@ subroutine read_rhotab(filename, rsize, rtab, rhotab, nread, polyk, gamma, rhoc,
  endif
  ! Loop over 'n' lines: r and density separated by space
  do i = 1,nread
-    read(iunit,*, iostat=ierr) rtab(i), rhotab(i)
+    read(iunit,*,iostat=ierr) rtab(i), rhotab(i)
     if (ierr /= 0) then
        call error('extern_densityprofile','Error reading data from '//trim(filename))
        return

src/main/extern_spiral.f90

@@ -377,20 +377,20 @@ subroutine initialise_spiral(ierr)
        spiralsum(jj)=0.0d0
        !-Loop over spheroids
        do j=1,Nt
-          Rspheroids(jj,j) = Ri+(DBLE(j)-1.d0)*d_0
+          Rspheroids(jj,j) = Ri+(dble(j)-1.d0)*d_0
           shapefn(jj,j)    = (cotalpha/Nshape) * &
-    log(1.d0+(Rspheroids(jj,j)/Rsarms)**Nshape) + jj*2.0d0*pi/DBLE(NNi)
+    log(1.d0+(Rspheroids(jj,j)/Rsarms)**Nshape) + jj*2.0d0*pi/dble(NNi)
           !print*,jj,j,Rspheroids(jj,j),shapefn(jj,j)
           select case(iarms)
           case(2,4)
              !--For a linear density drop off from galactic centre:
              den0(jj,j)       = (Rf-Rspheroids(jj,j))*3.d0*Mspiral &
-      / (DBLE(NNi)*pi*a_0*a_0*c_0)
+      / (dble(NNi)*pi*a_0*a_0*c_0)
              spiralsum(jj) = spiralsum(jj) + (Rf-Rspheroids(jj,j))
           case(3)
              !--For a log density drop off from galactic centre:
              den0(jj,j)       = exp((Ri-Rspheroids(jj,j))/Rl)*3.d0*Mspiral &
-      / (DBLE(NNi)*pi*a_0*a_0*c_0)
+      / (dble(NNi)*pi*a_0*a_0*c_0)
              spiralsum(jj) = spiralsum(jj) + exp((Ri-Rspheroids(jj,j))/Rl)
           end select
        enddo
@@ -420,9 +420,9 @@ subroutine initialise_spiral(ierr)
  case(1)
     potfilename = 'pot3D.bin'
     if (id==master) print*,'Reading in potential from an external file (BINARY): ',potfilename
-    open (unit =1, file = TRIM(potfilename), status='old', form='UNFORMATTED', access='SEQUENTIAL', iostat=ios)
+    open(unit=1,file=trim(potfilename),status='old',form='UNFORMATTED',access='SEQUENTIAL',iostat=ios)
     if (ios /= 0 .and. id==master) then
-       print*, 'Error opening file:', TRIM(potfilename)
+       print*, 'Error opening file:', trim(potfilename)
     endif
     !Read in the grid lengths if they exist in the header.
     read(1) potlenz,potlenx,potleny
@@ -1281,7 +1281,7 @@ subroutine Wang_bar(ri,phii,thetai,pot)
     allocate(PlmA(l+1))
     call legendre_associated(l,m,cos(thetai),PlmA)
     Plm=PlmA(l+1)
-    thisphi = Anlm(i) * (s**REAL(l))/((1.+s)**(2.*REAL(l)+1.)) * Gnl * Plm * cos(REAL(m)*(phii))
+    thisphi = Anlm(i) * (s**real(l))/((1.+s)**(2.*real(l)+1.)) * Gnl * Plm * cos(real(m)*(phii))
     AlmnSum = AlmnSum + thisphi
 
     deallocate(GnlA,PlmA)

src/main/forcing.F90

@@ -1076,8 +1076,8 @@ subroutine read_stirring_data_from_file(infile, time, timeinfile)
 
   my_file = find_phantom_datafile(infile,'forcing')
 
-  open (unit=42, file=my_file, iostat=ierr, status='old', action='read', &
-        access='sequential', form='unformatted')
+  open(unit=42,file=my_file,iostat=ierr,status='old',action='read', &
+        access='sequential',form='unformatted')
   ! header contains number of times and number of modes, end time, autocorrelation time, ...
   if (ierr==0) then
      if (Debug) write (*,'(A)') 'reading header...'

src/main/kernel_WendlandC2.f90

@@ -17,9 +17,6 @@ module kernel
 !
 ! :Dependencies: physcon
 !
-! :Generated: 2024-04-08 15:21:28.635699
-!
-!--------------------------------------------------------------------------
  use physcon, only:pi
  implicit none
  character(len=17), public :: kernelname = 'Wendland 2/3D C^2'

src/main/kernel_WendlandC4.f90

@@ -17,9 +17,6 @@ module kernel
 !
 ! :Dependencies: physcon
 !
-! :Generated: 2024-04-08 15:21:39.886138
-!
-!--------------------------------------------------------------------------
  use physcon, only:pi
  implicit none
  character(len=17), public :: kernelname = 'Wendland 2/3D C^4'

src/main/kernel_WendlandC6.f90

@@ -17,9 +17,6 @@ module kernel
 !
 ! :Dependencies: physcon
 !
-! :Generated: 2024-04-08 15:21:50.637883
-!
-!--------------------------------------------------------------------------
  use physcon, only:pi
  implicit none
  character(len=17), public :: kernelname = 'Wendland 2/3D C^6'

src/main/kernel_quartic.f90

@@ -17,9 +17,6 @@ module kernel
 !
 ! :Dependencies: physcon
 !
-! :Generated: 2024-04-08 15:20:17.993158
-!
-!--------------------------------------------------------------------------
  use physcon, only:pi
  implicit none
  character(len=11), public :: kernelname = 'M_5 quartic'

src/main/kernel_quintic.f90

@@ -17,9 +17,6 @@ module kernel
 !
 ! :Dependencies: physcon
 !
-! :Generated: 2024-04-22 14:12:57.936556
-!
-!--------------------------------------------------------------------------
  use physcon, only:pi
  implicit none
  character(len=11), public :: kernelname = 'M_6 quintic'