Merge pull request #298 from matt-frey/2d-merge-periodic-on-posm

Apply pull request #297 to master branch

src/2d/parcels/parcel_container.f90

@@ -33,11 +33,23 @@ module parcel_container
         ! @param[in] x1 first horizontal position
         ! @param[in] x2 second horizontal position
         ! @returns delx = x1 - x2
+        ! WARNING input needs to be between lower and upper (see debug statement)
+#ifndef NDEBUG
+        function get_delx(x1, x2) result (delx)
+#else
         elemental function get_delx(x1, x2) result (delx)
+#endif
             double precision, intent(in) :: x1, x2
             double precision             :: delx
 
             delx = x1 - x2
+#ifndef NDEBUG
+            if ((x1 < lower(1)) .or. (x2 < lower(1)) .or. (x1 > upper(1)) .or. (x2 > upper(1))) then
+                write(*,*) 'point outside domain was fed into get_delx'
+                write(*,*) 'x1, x2, lower(1), upper(1)'
+                write(*,*) x1, x2, lower(1), upper(1)
+            endif
+#endif
             ! works across periodic edge
             delx = delx - extent(1) * dble(int(delx * hli(1)))
         end function get_delx

src/2d/parcels/parcel_merge.f90

@@ -167,12 +167,20 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm)
                 ! temporary scalar containing 1 / vm(m)
                 vmerge = one / vm(m)
 
+                ! need to sanitise input and output, but first to determine input
+                posm(m, 1) = - vmerge * posm(m, 1)
+
+                call apply_periodic_bc(posm(m, :))
+
                 ! x centre of merged parcel, modulo periodicity
-                posm(m, 1) = get_delx(x0(m), - vmerge * posm(m, 1))
+                posm(m, 1) = get_delx(x0(m), posm(m, 1))
 
                 ! z centre of merged parcel
                 posm(m, 2) = vmerge * posm(m, 2)
 
+                ! need to correct position
+                call apply_periodic_bc(posm(m, :))
+
                 ! buoyancy and humidity
                 buoym(m) = vmerge * buoym(m)
 #ifndef ENABLE_DRY_MODE