Draft PR for 2D parcel types, for reference

configure.ac

@@ -24,6 +24,7 @@ AC_CONFIG_FILES([
     src/mpi/Makefile
     src/netcdf/Makefile
     src/fft/Makefile
+    src/parcel_types/Makefile
     src/2d/Makefile
     src/3d/Makefile
     models/Makefile

examples/robert2d.nml

@@ -1,9 +1,9 @@
 &MODELS
 
  model      = 'Robert'
- ncfname    = 'robert2d_1024x1536.nc'
+ ncfname    = 'robert2d.nc'
 
- box%ncells = 1024, 1536
+ box%ncells = 64, 96
  box%extent = 1000, 1500
  box%origin = -500.0, 0.0
 

src/2d/Makefile.am

@@ -1,7 +1,8 @@
 AM_FCFLAGS =				\
 	-I $(top_builddir)/src/mpi      \
 	-I $(top_builddir)/src/utils	\
-	-I $(top_builddir)/src/netcdf
+	-I $(top_builddir)/src/netcdf   \
+	-I $(top_builddir)/src/parcel_types
 
 bin_PROGRAMS = epic2d
 epic2d_SOURCES = 				\
@@ -14,6 +15,7 @@ epic2d_SOURCES = 				\
 	fields/field_diagnostics.f90		\
 	fields/field_netcdf.f90			\
 	fields/field_diagnostics_netcdf.f90	\
+	parcels/dynamic_parcels.f90		\
 	parcels/parcel_ellipse.f90		\
 	parcels/parcel_container.f90 		\
 	parcels/parcel_bc.f90			\
@@ -35,7 +37,8 @@ epic2d_SOURCES = 				\
 epic2d_LDADD = 						\
 	$(top_builddir)/src/utils/libepic_utils.la 	\
 	$(top_builddir)/src/mpi/libepic_mpi.la          \
-	$(top_builddir)/src/netcdf/libepic_netcdf.la
+	$(top_builddir)/src/netcdf/libepic_netcdf.la    \
+	$(top_builddir)/src/parcel_types/libepic_parcel_types.la
 
 clean-local:
 	rm -f *.mod

src/2d/epic2d.f90

@@ -4,7 +4,7 @@
 program epic2d
     use constants, only : zero
     use timer
-    use parcel_container
+    use dynamic_parcels, only : parcels
     use parcel_bc
     use parcel_split, only : split_ellipses, split_timer
     use parcel_merge, only : merge_ellipses, merge_timer
@@ -138,7 +138,7 @@ end subroutine run
 
         subroutine post_run
             use options, only : output
-            call parcel_dealloc
+            call parcels%dealloc
             call ls_rk4_dealloc
 
             call stop_timer(epic_timer)

src/2d/parcels/dynamic_parcels.f90

@@ -0,0 +1,9 @@
+module  dynamic_parcels
+use parcel_types, only : idealised_parcel_type
+implicit none
+
+    type(idealised_parcel_type) :: parcels
+    integer :: n_parcels
+
+end module dynamic_parcels
+

src/2d/parcels/parcel_bc.f90

@@ -6,7 +6,7 @@
 module parcel_bc
     use constants, only : zero, two
     use parameters, only : lower, upper, extent, hli, center
-    use parcel_container, only : n_parcels
+    use dynamic_parcels, only : n_parcels
     use omp_lib
     implicit none
 

src/2d/parcels/parcel_container.f90

@@ -54,32 +54,7 @@ elemental function get_delx(x1, x2) result (delx)
             delx = delx - extent(1) * dble(nint(delx * extenti(1)))
         end function get_delx
 
-
-        ! Overwrite parcel n with parcel m
-        ! @param[in] n index of parcel to be replaced
-        ! @param[in] m index of parcel used to replace parcel at index n
-        ! @pre n and m must be valid parcel indices
-        subroutine parcel_replace(n, m)
-            integer, intent(in) :: n, m
-
-#ifdef ENABLE_VERBOSE
-            if (verbose) then
-                print '(a19, i0, a6, i0)', '    replace parcel ', n, ' with ', m
-            endif
-#endif
-
-            parcels%position(:, n) = parcels%position(:, m)
-
-            parcels%vorticity(n) = parcels%vorticity(m)
-
-            parcels%volume(n)  = parcels%volume(m)
-            parcels%buoyancy(n) = parcels%buoyancy(m)
-#ifndef ENABLE_DRY_MODE
-            parcels%humidity(n) = parcels%humidity(m)
-#endif
-            parcels%B(:, n) = parcels%B(:, m)
-
-        end subroutine parcel_replace
+        ! KEEP AROUND FOR NOW FOR UNIT TESTS
 
         ! Allocate parcel memory
         ! @param[in] num number of parcels

src/2d/parcels/parcel_correction.f90

@@ -20,7 +20,7 @@ module parcel_correction
     use constants
     use parameters, only : vcelli, nx, nz, dx, dxi
 
-    use parcel_container
+    use dynamic_parcels, only : parcels, n_parcels
 
     use timer, only : start_timer, stop_timer
 

src/2d/parcels/parcel_diagnostics.f90

@@ -5,7 +5,7 @@ module parcel_diagnostics
     use constants, only : zero, one, f12
     use merge_sort
     use parameters, only : extent, lower, vcell, vmin, nx, nz, vdomaini
-    use parcel_container, only : parcels, n_parcels
+    use dynamic_parcels, only : parcels, n_parcels
     use parcel_ellipse
     use omp_lib
     use physics, only : peref, ape_calculation
@@ -101,8 +101,8 @@ subroutine calculate_parcel_diagnostics(velocity)
             ! this way the result is reproducible
             bmin = minval(parcels%buoyancy(1:n_parcels))
             bmax = maxval(parcels%buoyancy(1:n_parcels))
-            vormin = minval(parcels%vorticity(1:n_parcels))
-            vormax = maxval(parcels%vorticity(1:n_parcels))
+            vormin = minval(parcels%vorticity(1, 1:n_parcels))
+            vormax = maxval(parcels%vorticity(1, 1:n_parcels))
 
             lsum = zero
             l2sum = zero
@@ -155,7 +155,7 @@ subroutine calculate_parcel_diagnostics(velocity)
                     n_small = n_small + 1
                 endif
 
-                rms_zeta = rms_zeta + vol * parcels%vorticity(n) ** 2
+                rms_zeta = rms_zeta + vol * parcels%vorticity(1, n) ** 2
 
             enddo
             !$omp end do
@@ -242,7 +242,7 @@ subroutine straka_diagnostics
                     xzbv = xzbv + bv * parcels%position(1, n) * parcels%position(2, n)
 
 
-                    vv = parcels%vorticity(n) * parcels%volume(n)
+                    vv = parcels%vorticity(1, n) * parcels%volume(n)
                     vvsum = vvsum + vv
                     xvv = xvv + vv * parcels%position(1, n)
                     zvv = zvv + vv * parcels%position(2, n)

src/2d/parcels/parcel_diagnostics_netcdf.f90

@@ -7,7 +7,7 @@ module parcel_diagnostics_netcdf
     use netcdf_utils
     use netcdf_writer
     use netcdf_reader
-    use parcel_container, only : parcels, n_parcels
+    use dynamic_parcels, only : parcels, n_parcels
     use parcel_diagnostics
     use parameters, only : lower, extent, nx, nz
     use config, only : package_version, cf_version

src/2d/parcels/parcel_init.f90

@@ -4,7 +4,7 @@
 module parcel_init
     use options, only : parcel, output, verbose, field_tol
     use constants, only : zero, two, one, f12
-    use parcel_container, only : parcels, n_parcels
+    use dynamic_parcels, only : parcels, n_parcels
     use parcel_ellipse, only : get_ab, get_B22, get_eigenvalue
     use parcel_split, only : split_ellipses
     use parcel_interpl, only : bilinear, ngp
@@ -58,7 +58,6 @@ subroutine init_parcels(fname, tol)
                 stop
             endif
 
-
             call init_regular_positions
 
             ! initialize the volume of each parcel
@@ -92,11 +91,11 @@ subroutine init_parcels(fname, tol)
             !$omp parallel default(shared)
             !$omp do private(n)
             do n = 1, n_parcels
-                parcels%vorticity(n) = zero
+                parcels%vorticity(1, n) = zero
                 parcels%buoyancy(n) = zero
-#ifndef ENABLE_DRY_MODE
-                parcels%humidity(n) = zero
-#endif
+                if(parcels%is_moist) then
+                    parcels%humidity(n) = zero
+                endif
             enddo
             !$omp end do
             !$omp end parallel
@@ -239,7 +238,7 @@ subroutine init_from_grids(ncfname, tol)
             if (has_dataset(ncid, 'vorticity')) then
                 buffer = zero
                 call read_netcdf_dataset(ncid, 'vorticity', buffer(0:nz, :), start=start, cnt=cnt)
-                call gen_parcel_scalar_attr(buffer, tol, parcels%vorticity)
+                call gen_parcel_scalar_attr(buffer, tol, parcels%vorticity(1, :))
             endif
 
 

src/2d/parcels/parcel_interpl.f90

@@ -7,7 +7,7 @@ module parcel_interpl
     use timer, only : start_timer, stop_timer
     use parameters, only : nx, nz, vmin
     use options, only : parcel
-    use parcel_container, only : parcels, n_parcels
+    use dynamic_parcels, only : parcels, n_parcels
     use parcel_bc, only : apply_periodic_bc
     use parcel_ellipse
     use fields
@@ -142,9 +142,6 @@ subroutine par2grid
             double precision :: points(2, 2)
             integer          :: n, p, l, i, j
             double precision :: pvol, weight, btot
-#ifndef ENABLE_DRY_MODE
-            double precision :: q_c
-#endif
 
             call start_timer(par2grid_timer)
 
@@ -168,17 +165,8 @@ subroutine par2grid
             do n = 1, n_parcels
                 pvol = parcels%volume(n)
 
-#ifndef ENABLE_DRY_MODE
-                ! liquid water content
-                q_c = parcels%humidity(n) &
-                    - q_0 * exp(lambda_c * (lower(2) - parcels%position(2, n)))
-                q_c = max(zero, q_c)
+                call parcels%get_buoyancy(n, btot)
 
-                ! total buoyancy (including effects of latent heating)
-                btot = parcels%buoyancy(n) + glat * q_c
-#else
-                btot = parcels%buoyancy(n)
-#endif
                 points = get_ellipse_points(parcels%position(:, n), &
                                             pvol, parcels%B(:, n))
 
@@ -205,12 +193,12 @@ subroutine par2grid
                         weight = f12 * weights(l) * pvol
 
                         vortg(js(l), is(l)) = vortg(js(l), is(l)) &
-                                            + weight * parcels%vorticity(n)
+                                            + weight * parcels%vorticity(1, n)
 
 #ifndef ENABLE_DRY_MODE
-                        dbuoyg(js(l), is(l)) = dbuoyg(js(l), is(l)) &
+                            dbuoyg(js(l), is(l)) = dbuoyg(js(l), is(l)) &
                                              + weight * parcels%buoyancy(n)
-                        humg(js(l), is(l)) = humg(js(l), is(l)) &
+                            humg(js(l), is(l)) = humg(js(l), is(l)) &
                                              + weight * parcels%humidity(n)
 #endif
                         tbuoyg(js(l), is(l)) = tbuoyg(js(l), is(l)) &