Bug found in the MQUASILAG=1 option, xmass array accessed with out-of-bounds...

Bug found in the MQUASILAG=1 option, xmass array accessed with out-of-bounds index. Quick fix implemented (see changes in timemanager and advance).

src/advance.f90

@@ -492,14 +492,15 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
   !*************************************************************************
 
       if (mdomainfill.eq.0) then
-        do nsp=1,nspec
-          if (xmass(nrelpoint,nsp).gt.eps2) goto 887
-        end do
-887     nsp=min(nsp,nspec)
-!!$        if (density(nsp).gt.0.) &
-!!$             call get_settling(itime,xts,yts,zt,nsp,settling)    !old
-        if (density(nsp).gt.0.) &
-             call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
+! ESO 05.2015  Changed this to fix MQUASILAG option, where nrelpoint is
+!              particle number and thus xmass array goes out of bounds
+!        do nsp=1,nspec
+!           if (xmass(nrelpoint,nsp).gt.eps2) goto 887
+!         end do
+! 887     nsp=min(nsp,nspec)
+        if (nspec.eq.1.and.density(1).gt.0.) then
+          call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
+        end if
         w=w+settling
       endif
 
@@ -653,14 +654,17 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
 
 
     if (mdomainfill.eq.0) then
-      do nsp=1,nspec
-        if (xmass(nrelpoint,nsp).gt.eps2) goto 888
-      end do
-888   nsp=min(nsp,nspec)
-!!$      if (density(nsp).gt.0.) &
-!!$           call get_settling(itime,xts,yts,zt,nsp,settling)    !old
-      if (density(nsp).gt.0.) &
-           call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
+! ESO 05.2015  Changed this to fix MQUASILAG option, where nrelpoint is
+!              particle number and thus xmass array goes out of bounds
+
+!      do nsp=1,nspec
+!         if (xmass(nrelpoint,nsp).gt.eps2) goto 888
+!       end do
+! 888   nsp=min(nsp,nspec)
+!        if (density(nsp).gt.0.) then
+      if (nspec.eq.1.and.density(1).gt.0.) then
+        call get_settling(itime,real(xt),real(yt),zt,nsp,settling)  !bugfix
+      end if
       w=w+settling
     endif
 
@@ -857,14 +861,16 @@ subroutine advance(itime,nrelpoint,ldt,up,vp,wp, &
   endif
 
   if (mdomainfill.eq.0) then
-    do nsp=1,nspec
-      if (xmass(nrelpoint,nsp).gt.eps2) goto 889
-    end do
-889   nsp=min(nsp,nspec)
-!!$    if (density(nsp).gt.0.) &
-!!$         call get_settling(itime+ldt,xts,yts,zt,nsp,settling)    !old
-    if (density(nsp).gt.0.) &
-         call get_settling(itime+ldt,real(xt),real(yt),zt,nsp,settling)  !bugfix
+! ESO 05.2015  Changed this to fix MQUASILAG option, where nrelpoint is
+!              particle number and thus xmass array goes out of bounds
+!    do nsp=1,nspec
+!       if (xmass(nrelpoint,nsp).gt.eps2) goto 889
+!     end do
+! 889   nsp=min(nsp,nspec)
+!      if (density(nsp).gt.0.) then
+    if (nspec.eq.1.and.density(1).gt.0.) then
+      call get_settling(itime+ldt,real(xt),real(yt),zt,nsp,settling)  !bugfix
+    end if
     w=w+settling
   endif
 

src/par_mod.f90

@@ -185,7 +185,7 @@ module par_mod
   ! Maximum number of particles, species, and similar
   !**************************************************
 
-  integer,parameter :: maxpart=20000000
+  integer,parameter :: maxpart=1000000
   integer,parameter :: maxspec=6
   real,parameter :: minmass=0.0001
 

src/partoutput_short.f90

@@ -145,7 +145,7 @@ subroutine partoutput_short(itime)
        (i4dump(i),(idump(j,i),j=1,3),i=1,numshortout)
 
 
-  write(*,*) numshortout,numshortall
+  !write(*,*) numshortout,numshortall
 
   close(unitshortpart)
 

src/timemanager.f90

@@ -588,8 +588,8 @@ subroutine timemanager
               xmass1(j,ks)=xmass1(j,ks)*decfact
             endif
 
-
-            if (mdomainfill.eq.0) then
+! Skip check on mass fraction when npoint represents particle number
+            if (mdomainfill.eq.0.and.mquasilag.eq.0) then
               if (xmass(npoint(j),ks).gt.0.) &
                    xmassfract=max(xmassfract,real(npart(npoint(j)))* &
                    xmass1(j,ks)/xmass(npoint(j),ks))
@@ -604,8 +604,8 @@ subroutine timemanager
                    !CGZ-lifetime: Check mass fraction left/save lifetime
 !ZHG 2015
             else
-              xmassfract=1.
-            endif
+              xmassfract=1.0
+            end if
           end do
 
           if (xmassfract.lt.minmass) then   ! terminate all particles carrying less mass

src/timemanager_mpi.f90

@@ -545,9 +545,14 @@ subroutine timemanager
 
 ! Decide whether to write an estimate of the number of particles released, 
 ! or exact number (require MPI reduce operation)
-        numpart_tot_mpi = numpart*mp_partgroup_np
+        if (mp_dev_mode) then
+          numpart_tot_mpi = numpart
+        else
+          numpart_tot_mpi = numpart*mp_partgroup_np
+        end if
 
-        if (mp_exact_numpart.and..not.(lmpreader.and.lmp_use_reader)) then
+        if (mp_exact_numpart.and..not.(lmpreader.and.lmp_use_reader).and.&
+             &.not.mp_dev_mode) then
           call MPI_Reduce(numpart, numpart_tot_mpi, 1, MPI_INTEGER, MPI_SUM, id_root, &
                & mp_comm_used, mp_ierr)
         endif
@@ -737,8 +742,8 @@ subroutine timemanager
               xmass1(j,ks)=xmass1(j,ks)*decfact
             endif
 
-
-            if (mdomainfill.eq.0) then
+! Skip check on mass fraction when npoint represents particle number
+            if (mdomainfill.eq.0.and.mquasilag.eq.0) then
               if (xmass(npoint(j),ks).gt.0.)then 
                    xmassfract=max(xmassfract,real(npart(npoint(j)))* &
                    xmass1(j,ks)/xmass(npoint(j),ks))