From b366d528a095acdac9f9c456f7f01170bf3a3ca1 Mon Sep 17 00:00:00 2001 From: Steven Boeing Date: Thu, 19 Jun 2025 10:04:49 +0100 Subject: [PATCH 1/7] First attempt. Runs but in trouble when splitting --- configure.ac | 1 + examples/robert2d.nml | 4 ++-- src/2d/Makefile.am | 7 +++++-- src/2d/epic2d.f90 | 4 ++-- src/2d/parcels/parcel_bc.f90 | 2 +- src/2d/parcels/parcel_correction.f90 | 2 +- src/2d/parcels/parcel_diagnostics.f90 | 10 +++++----- src/2d/parcels/parcel_diagnostics_netcdf.f90 | 2 +- src/2d/parcels/parcel_init.f90 | 7 +++---- src/2d/parcels/parcel_interpl.f90 | 4 ++-- src/2d/parcels/parcel_merge.f90 | 21 ++++++++++---------- src/2d/parcels/parcel_netcdf.f90 | 6 +++--- src/2d/parcels/parcel_split.f90 | 7 ++++--- src/2d/stepper/ls_rk4.f90 | 4 ++-- src/2d/utils/utils.f90 | 8 ++++++-- src/Makefile.am | 2 +- 16 files changed, 49 insertions(+), 42 deletions(-) diff --git a/configure.ac b/configure.ac index 6cef59532..1aad4e60a 100644 --- a/configure.ac +++ b/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 diff --git a/examples/robert2d.nml b/examples/robert2d.nml index b05ee1e84..537279a3b 100644 --- a/examples/robert2d.nml +++ b/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 diff --git a/src/2d/Makefile.am b/src/2d/Makefile.am index acaff8dfa..0e4033474 100644 --- a/src/2d/Makefile.am +++ b/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 diff --git a/src/2d/epic2d.f90 b/src/2d/epic2d.f90 index 3ac58a087..9dc3121d4 100644 --- a/src/2d/epic2d.f90 +++ b/src/2d/epic2d.f90 @@ -4,7 +4,7 @@ program epic2d use constants, only : zero use timer - use parcel_container + use dynamic_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) diff --git a/src/2d/parcels/parcel_bc.f90 b/src/2d/parcels/parcel_bc.f90 index a305c1aa9..da4bf82b6 100644 --- a/src/2d/parcels/parcel_bc.f90 +++ b/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 diff --git a/src/2d/parcels/parcel_correction.f90 b/src/2d/parcels/parcel_correction.f90 index 8ac68d18b..b99e19625 100644 --- a/src/2d/parcels/parcel_correction.f90 +++ b/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 diff --git a/src/2d/parcels/parcel_diagnostics.f90 b/src/2d/parcels/parcel_diagnostics.f90 index 385fbc702..6f7ac5394 100644 --- a/src/2d/parcels/parcel_diagnostics.f90 +++ b/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) diff --git a/src/2d/parcels/parcel_diagnostics_netcdf.f90 b/src/2d/parcels/parcel_diagnostics_netcdf.f90 index 9bdf50bf4..0c5829377 100644 --- a/src/2d/parcels/parcel_diagnostics_netcdf.f90 +++ b/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 diff --git a/src/2d/parcels/parcel_init.f90 b/src/2d/parcels/parcel_init.f90 index a626411b8..52cd591b1 100644 --- a/src/2d/parcels/parcel_init.f90 +++ b/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,7 +91,7 @@ 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 @@ -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 diff --git a/src/2d/parcels/parcel_interpl.f90 b/src/2d/parcels/parcel_interpl.f90 index 0aa42264f..04472cd9e 100644 --- a/src/2d/parcels/parcel_interpl.f90 +++ b/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 @@ -205,7 +205,7 @@ 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)) & diff --git a/src/2d/parcels/parcel_merge.f90 b/src/2d/parcels/parcel_merge.f90 index d4bae92dd..df27ddc44 100644 --- a/src/2d/parcels/parcel_merge.f90 +++ b/src/2d/parcels/parcel_merge.f90 @@ -5,10 +5,9 @@ module parcel_merge use parcel_nearest use constants, only : pi, zero, one, two, four - use parcel_container, only : parcel_container_type & - , n_parcels & - , parcel_replace & - , get_delx + use parcel_container, only : get_delx + use dynamic_parcels, only : n_parcels, parcels + use parcel_types, only : idealised_parcel_type use parcel_ellipse, only : get_B22, get_ab use options, only : parcel, verbose use parcel_bc @@ -27,7 +26,7 @@ module parcel_merge ! parcels which are close by. ! @param[inout] parcels is the parcel container subroutine merge_ellipses(parcels) - type(parcel_container_type), intent(inout) :: parcels + type(idealised_parcel_type), intent(inout) :: parcels integer, allocatable, dimension(:) :: isma integer, allocatable, dimension(:) :: iclo integer :: n_merge ! number of merges @@ -73,7 +72,7 @@ end subroutine merge_ellipses ! @param[out] B22m are the B22 matrix entries of the mergers ! @param[out] vm are the volumes of the mergers subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm) - type(parcel_container_type), intent(inout) :: parcels + type(idealised_parcel_type), intent(inout) :: parcels integer, intent(in) :: isma(0:) integer, intent(in) :: iclo(:) integer, intent(in) :: n_merge @@ -116,7 +115,7 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm) #ifndef ENABLE_DRY_MODE hum(l) = parcels%volume(ic) * parcels%humidity(ic) #endif - vortm(l) = parcels%volume(ic) * parcels%vorticity(ic) + vortm(l) = parcels%volume(ic) * parcels%vorticity(1, ic) B11m(l) = zero B12m(l) = zero @@ -143,7 +142,7 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm) #ifndef ENABLE_DRY_MODE hum(n) = hum(n) + parcels%volume(is) * parcels%humidity(is) #endif - vortm(n) = vortm(n) + parcels%volume(is) * parcels%vorticity(is) + vortm(n) = vortm(n) + parcels%volume(is) * parcels%vorticity(1, is) enddo ! Obtain the merged parcel centres @@ -201,7 +200,7 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm) #ifndef ENABLE_DRY_MODE parcels%humidity(ic) = hum(l) #endif - parcels%vorticity(ic) = vortm(l) + parcels%vorticity(1, ic) = vortm(l) endif @@ -232,7 +231,7 @@ end subroutine do_group_merge ! @param[in] iclo are the indices of the close parcels ! @param[in] n_merge is the array size of isma and iclo subroutine geometric_merge(parcels, isma, iclo, n_merge) - type(parcel_container_type), intent(inout) :: parcels + type(idealised_parcel_type), intent(inout) :: parcels integer, intent(in) :: isma(0:) integer, intent(in) :: iclo(:) integer, intent(in) :: n_merge @@ -307,7 +306,7 @@ subroutine pack_parcels(isma, n_merge) do while (m <= k) ! invalid parcel; overwrite *isma(m)* with last valid parcel *l* - call parcel_replace(isma(m), l) + call parcels%replace(isma(m), l) l = l - 1 diff --git a/src/2d/parcels/parcel_netcdf.f90 b/src/2d/parcels/parcel_netcdf.f90 index 9b86c1df9..d6cf55b7e 100644 --- a/src/2d/parcels/parcel_netcdf.f90 +++ b/src/2d/parcels/parcel_netcdf.f90 @@ -3,7 +3,7 @@ module parcel_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 parameters, only : nx, nz, extent, lower, max_num_parcels use config, only : package_version, cf_version use timer, only : start_timer, stop_timer @@ -208,7 +208,7 @@ subroutine write_netcdf_parcels(t) call write_netcdf_dataset(ncid, vol_id, parcels%volume(1:n_parcels), start, cnt) - call write_netcdf_dataset(ncid, vor_id, parcels%vorticity(1:n_parcels), start, cnt) + call write_netcdf_dataset(ncid, vor_id, parcels%vorticity(1, 1:n_parcels), start, cnt) call write_netcdf_dataset(ncid, buo_id, parcels%buoyancy(1:n_parcels), start, cnt) @@ -290,7 +290,7 @@ subroutine read_netcdf_parcels(fname) if (has_dataset(ncid, 'vorticity')) then l_valid = .true. call read_netcdf_dataset(ncid, 'vorticity', & - parcels%vorticity(1:n_parcels), start, cnt) + parcels%vorticity(1, 1:n_parcels), start, cnt) endif if (has_dataset(ncid, 'buoyancy')) then diff --git a/src/2d/parcels/parcel_split.f90 b/src/2d/parcels/parcel_split.f90 index cbe3c8b11..3891664e5 100644 --- a/src/2d/parcels/parcel_split.f90 +++ b/src/2d/parcels/parcel_split.f90 @@ -5,7 +5,8 @@ module parcel_split use options, only : verbose use constants, only : pi, three, f12, f14, f34 use parameters, only : vmax - use parcel_container, only : parcel_container_type, n_parcels + use parcel_types, only : idealised_parcel_type + use dynamic_parcels, only : n_parcels use parcel_bc, only : apply_reflective_bc use parcel_ellipse, only : get_eigenvalue & , get_eigenvector & @@ -24,7 +25,7 @@ module parcel_split ! @param[inout] parcels ! @param[in] threshold is the largest allowed aspect ratio subroutine split_ellipses(parcels, threshold) - type(parcel_container_type), intent(inout) :: parcels + type(idealised_parcel_type), intent(inout) :: parcels double precision, intent(in) :: threshold double precision :: B11 double precision :: B12 @@ -78,7 +79,7 @@ subroutine split_ellipses(parcels, threshold) parcels%B(:, n_thread_loc) = parcels%B(:, n) - parcels%vorticity(n_thread_loc) = parcels%vorticity(n) + parcels%vorticity(:, n_thread_loc) = parcels%vorticity(:, n) parcels%volume(n_thread_loc) = parcels%volume(n) parcels%buoyancy(n_thread_loc) = parcels%buoyancy(n) #ifndef ENABLE_DRY_MODE diff --git a/src/2d/stepper/ls_rk4.f90 b/src/2d/stepper/ls_rk4.f90 index f7b99c8d5..46909af7d 100644 --- a/src/2d/stepper/ls_rk4.f90 +++ b/src/2d/stepper/ls_rk4.f90 @@ -4,7 +4,7 @@ ! ============================================================================= module ls_rk4 use options, only : parcel - use parcel_container + use dynamic_parcels, only : parcels, n_parcels use parcel_bc use rk4_utils, only: get_B, get_time_step use utils, only : write_step @@ -164,7 +164,7 @@ subroutine ls_rk4_substep(dt, step) parcels%position(:, n) = parcels%position(:, n) & + cb * dt * delta_pos(:, n) - parcels%vorticity(n) = parcels%vorticity(n) + cb * dt * delta_vor(n) + parcels%vorticity(1, n) = parcels%vorticity(1, n) + cb * dt * delta_vor(n) parcels%B(:, n) = parcels%B(:, n) + cb * dt * delta_b(:, n) enddo !$omp end parallel do diff --git a/src/2d/utils/utils.f90 b/src/2d/utils/utils.f90 index ddc891627..0a046371c 100644 --- a/src/2d/utils/utils.f90 +++ b/src/2d/utils/utils.f90 @@ -17,7 +17,8 @@ module utils use parcel_diagnostics, only : calculate_parcel_diagnostics, calculate_peref use field_diagnostics, only : calculate_field_diagnostics use parcel_init, only : init_parcels - use parcel_container, only : n_parcels, parcel_alloc + use dynamic_parcels, only : parcels, n_parcels + use parcel_types, only : idealised_parcel_alloc use tri_inversion, only : vor2vel, vorticity_tendency use parcel_interpl, only : par2grid, grid2par use netcdf_reader, only : get_file_type, get_num_steps, get_time, get_netcdf_box @@ -213,7 +214,10 @@ end subroutine setup_domain_and_parameters subroutine setup_parcels character(len=16) :: file_type - call parcel_alloc(max_num_parcels) + call parcels%set_dimension(2) + parcels%is_moist=.true. + parcels%shape_type="ellipsoid2" + call parcels%alloc(max_num_parcels) if (l_restart) then call setup_restart(trim(restart_file), time%initial, file_type) diff --git a/src/Makefile.am b/src/Makefile.am index 1d9b6c307..d60cb7065 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,4 +1,4 @@ -SUBDIRS = mpi netcdf utils fft 2d +SUBDIRS = mpi netcdf parcel_types utils fft 2d if ENABLE_3D SUBDIRS += 3d From 3e282500c4a71f497136aa177e59956b3baad64b Mon Sep 17 00:00:00 2001 From: Steven Boeing Date: Thu, 19 Jun 2025 10:28:04 +0100 Subject: [PATCH 2/7] Nearest bugfix --- src/2d/parcels/parcel_nearest.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/2d/parcels/parcel_nearest.f90 b/src/2d/parcels/parcel_nearest.f90 index 0bd89fbe3..04affa83c 100644 --- a/src/2d/parcels/parcel_nearest.f90 +++ b/src/2d/parcels/parcel_nearest.f90 @@ -3,7 +3,8 @@ !============================================================================== module parcel_nearest use constants, only : pi, f12 - use parcel_container, only : parcels, n_parcels, get_delx + use dynamic_parcels, only : parcels, n_parcels + use parcel_container, only : get_delx use parameters, only : dx, dxi, vcell, hli, lower, extent, ncell, nx, nz, vmin, max_num_parcels use options, only : parcel use timer, only : start_timer, stop_timer From 651883f3e75c60f4ed85cadc3353803a3e230cd3 Mon Sep 17 00:00:00 2001 From: Steven Boeing Date: Thu, 19 Jun 2025 10:28:16 +0100 Subject: [PATCH 3/7] Adding the files for parcel types --- src/2d/parcels/dynamic_parcels.f90 | 9 + src/parcel_types/Makefile.am | 10 + src/parcel_types/parcel_container.f90 | 696 ++++++++++++++++++++++++++ src/parcel_types/parcel_ellipsoid.f90 | 163 ++++++ src/parcel_types/parcel_types.f90 | 207 ++++++++ 5 files changed, 1085 insertions(+) create mode 100644 src/2d/parcels/dynamic_parcels.f90 create mode 100644 src/parcel_types/Makefile.am create mode 100644 src/parcel_types/parcel_container.f90 create mode 100644 src/parcel_types/parcel_ellipsoid.f90 create mode 100644 src/parcel_types/parcel_types.f90 diff --git a/src/2d/parcels/dynamic_parcels.f90 b/src/2d/parcels/dynamic_parcels.f90 new file mode 100644 index 000000000..1208feb40 --- /dev/null +++ b/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 + diff --git a/src/parcel_types/Makefile.am b/src/parcel_types/Makefile.am new file mode 100644 index 000000000..fc2833a82 --- /dev/null +++ b/src/parcel_types/Makefile.am @@ -0,0 +1,10 @@ +AM_FCFLAGS = -I $(top_builddir)/src/mpi \ + -I $(top_builddir)/src/utils + +lib_LTLIBRARIES = libepic_parcel_types.la +libepic_parcel_types_la_SOURCES = \ + ../utils/datatypes.f90 \ + ../utils/armanip.f90 \ + parcel_container.f90 \ + parcel_ellipsoid.f90 \ + parcel_types.f90 diff --git a/src/parcel_types/parcel_container.f90 b/src/parcel_types/parcel_container.f90 new file mode 100644 index 000000000..70bc72828 --- /dev/null +++ b/src/parcel_types/parcel_container.f90 @@ -0,0 +1,696 @@ +module parcel_container + use datatypes, only : int64 + use armanip, only : resize_array + implicit none + + ! TODO for EPIC implementation + ! Replace commands by MPI versions + + integer :: resize_timer + + type attr_ptr + double precision, pointer :: aptr(:) + character(len=32) :: name + character(len=128) :: long_name + character(len=128) :: std_name + character(len=32) :: dtype ! will need to be MPI_Datatype + character(len=32) :: unit + logical :: write_to_netcdf = .false. + end type attr_ptr + + type int_attr_ptr + integer(kind=8), pointer :: aptr(:) + character(len=32) :: name + character(len=128) :: long_name + character(len=128) :: std_name + character(len=32) :: dtype ! will need to be MPI_Datatype + character(len=32) :: unit + logical :: write_to_netcdf = .false. + end type int_attr_ptr + + type, abstract :: base_parcel_type + double precision, allocatable, dimension(:,:) :: position + double precision, allocatable, dimension(:,:) :: delta_pos ! For time-integration of the position + ! Sedimentation may need to be done separately + type(attr_ptr), allocatable, dimension(:) :: attrib + type(int_attr_ptr), allocatable, dimension(:) :: int_attrib + integer :: attr_num ! number of parcel attributes for serialisation + integer :: int_attr_num ! number of integer(kind=8) parcel attributes for serialisation + integer :: local_num ! local number of parcels + integer(kind=int64) :: total_num ! global number of parcels (over all MPI ranks) + integer :: max_num ! capacity per attribute, i.e. maximum number of parcels + integer, private :: n_pos = -1 ! number of spatial dimensions + + contains + procedure :: base_alloc => base_parcel_alloc + procedure :: base_dealloc => base_parcel_dealloc + procedure :: base_resize => base_parcel_resize + procedure :: serialize => base_parcel_serialize + procedure :: deserialize => base_parcel_deserialize + procedure :: replace => base_parcel_replace + procedure :: pack => parcel_pack + procedure :: unpack => parcel_unpack + procedure :: delete => parcel_delete + procedure(parcel_alloc), deferred :: alloc + procedure(parcel_dealloc), deferred :: dealloc + procedure(parcel_resize), deferred :: resize + procedure :: print_me + procedure :: set_dimension + procedure :: register_attribute + procedure :: register_int_attribute + procedure :: reset_attribute + procedure :: reset_int_attribute + + end type + + ! This type is for parcels associated with dynamics (which will always have volume and voriticity in EPIC) + type, extends(base_parcel_type) :: dynamic_parcel_type ! add procedures + double precision, allocatable, dimension(:) :: volume + double precision, allocatable, dimension(:,:) :: vorticity + double precision, allocatable, dimension(:,:) :: delta_vor + double precision, allocatable, dimension(:) :: dilution + integer(kind=8), allocatable, dimension(:) :: label + logical :: has_labels = .false. + + contains + procedure :: alloc => dynamic_parcel_alloc + procedure :: dealloc => dynamic_parcel_dealloc + procedure :: resize => dynamic_parcel_resize + + end type + + interface + subroutine parcel_alloc(this, num) + import base_parcel_type + class(base_parcel_type), intent(inout) :: this + integer, intent(in) :: num + end subroutine parcel_alloc + + subroutine parcel_dealloc(this) + import base_parcel_type + class(base_parcel_type), intent(inout) :: this + end subroutine parcel_dealloc + + subroutine parcel_resize(this, new_size) + import base_parcel_type + class(base_parcel_type), intent(inout) :: this + integer, intent(in) :: new_size + end subroutine parcel_resize + + logical pure function parcel_is_small(this, n) + import :: base_parcel_type + class(base_parcel_type), intent(in) :: this + integer, intent(in) :: n + end function parcel_is_small + + end interface + + interface try_deallocate + module procedure :: try_deallocate_1d + module procedure :: try_deallocate_1d_integer + module procedure :: try_deallocate_2d + end interface try_deallocate + + contains + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine try_deallocate_1d(in_array) + double precision, allocatable, dimension(:) :: in_array + + if (allocated(in_array)) then + deallocate(in_array) + endif + + end subroutine try_deallocate_1d + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine try_deallocate_1d_integer(in_array) + integer(kind=8), allocatable, dimension(:) :: in_array + + if (allocated(in_array)) then + deallocate(in_array) + endif + + end subroutine try_deallocate_1d_integer + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine try_deallocate_2d(in_array) + double precision, allocatable, dimension(:,:) :: in_array + + if (allocated(in_array)) then + deallocate(in_array) + endif + + end subroutine try_deallocate_2d + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + ! Allocate parcel memory + ! ATTENTION: Extended types must allocate additional parcel attributes + ! in their own routine. + ! @param[in] num number of parcels + subroutine base_parcel_alloc(this, num) + class(base_parcel_type), intent(inout), target :: this + integer, intent(in) :: num + character(len=1) :: dir(3) + integer :: n + + this%max_num = num + this%local_num = num ! Usually this is done in parcel_init + this%attr_num = 0 + this%int_attr_num = 0 + + if (this%n_pos > 3) then + print *, "Only 3 dimensions allowed." + stop + endif + + dir = (/'x', 'y', 'z'/) + + allocate(this%position(this%n_pos, num)) + allocate(this%delta_pos(this%n_pos, num)) + + do n = 1, this%n_pos + call this%register_attribute(this%position(n, :), dir(n) // "_position", "m") + call this%register_attribute(this%delta_pos(n, :), dir(n) // "_position_rk_tendency", "m/s") + enddo + + end subroutine base_parcel_alloc + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + ! Deallocate parcel memory + ! ATTENTION: Extended types must deallocate additional parcel attributes + ! in their own routine. + subroutine base_parcel_dealloc(this) + class(base_parcel_type), intent(inout) :: this + + this%local_num = 0 + this%total_num = 0 + this%max_num = 0 + this%attr_num = 0 + this%int_attr_num = 0 + + call try_deallocate(this%position) + call try_deallocate(this%delta_pos) + + if (allocated(this%attrib)) then + deallocate(this%attrib) + endif + + if (allocated(this%int_attrib)) then + deallocate(this%int_attrib) + endif + + end subroutine base_parcel_dealloc + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine base_parcel_resize(this, new_size) + class(base_parcel_type), intent(inout), target :: this + integer, intent(in) :: new_size + character(len=1) :: dir(3) + integer :: n + + if (new_size < this%local_num) then + print *, "in parcel_container::base_parcel_resize: losing parcels when resizing." + stop + endif + + this%max_num = new_size + + if (this%n_pos > 3) then + print *, "Only 3 dimensions allowed." + stop + endif + + dir = (/'x', 'y', 'z'/) + + call resize_array(this%position, new_size, this%local_num) + call resize_array(this%delta_pos, new_size, this%local_num) + + do n = 1, this%n_pos + call this%reset_attribute(this%position(n, :), dir(n) // "_position") + call this%reset_attribute(this%delta_pos(n, :), dir(n) // "_position_rk_tendency") + enddo + + end subroutine base_parcel_resize + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine set_dimension(this, n_dim) + class(base_parcel_type), intent(inout) :: this + integer, intent(in) :: n_dim + + this%n_pos = n_dim + + if (this%n_pos /= -1) then + print *, "WARNING: Dimension already set." + return + endif + + if ((n_dim < 1) .or. (n_dim > 3)) then + print *, "Only 1-, 2- or 3-dimensional." + stop + endif + + end subroutine set_dimension + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine register_attribute(this, attr, name, unit) + class(base_parcel_type), intent(inout) :: this + double precision, dimension(:), target, intent(inout) :: attr + character(len=*) :: name + character(len=*) :: unit + type(attr_ptr), allocatable, dimension(:) :: tmp + integer :: n_attr, n, j + + attr=0.0 ! set to zero when registring, for safety + + ! check if name is unique + do j = 1, this%attr_num + if(trim(this%attrib(j)%name) == trim(name)) then + print *, "Attribute name not unique" + stop + end if + end do + + do j = 1, this%int_attr_num + if(trim(this%int_attrib(j)%name) == trim(name)) then + print *, "Attribute name not unique" + stop + end if + end do + + this%attr_num = this%attr_num + 1 + + if (.not. allocated(this%attrib)) then + allocate(this%attrib(1)) + this%attrib(1)%aptr => attr + this%attrib(1)%name = name + this%attrib(1)%unit = unit + else + n_attr = size(this%attrib) + + allocate(tmp(n_attr+1)) + + do n = 1, n_attr + tmp(n)%aptr => this%attrib(n)%aptr + tmp(n)%name = this%attrib(n)%name + tmp(n)%unit = this%attrib(n)%unit + enddo + + tmp(n_attr+1)%aptr => attr + tmp(n_attr+1)%name = name + tmp(n_attr+1)%unit = unit + + deallocate(this%attrib) + + call move_alloc(from=tmp, to=this%attrib) + + endif + + end subroutine register_attribute + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine reset_attribute(this, attr, name) + class(base_parcel_type), intent(inout) :: this + double precision, dimension(:), target, intent(inout) :: attr + character(len=*) :: name + integer :: j + logical :: l_success + + l_success= .false. + + ! check if name is unique + do j = 1, this%attr_num + if(trim(this%attrib(j)%name) == trim(name)) then + this%attrib(j)%aptr => attr + l_success=.true. + exit + end if + end do + + if(.not. l_success) then + print *, "Attribute not found." + stop + end if + + end subroutine reset_attribute + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine register_int_attribute(this, int_attr, name, unit) + class(base_parcel_type), intent(inout) :: this + integer(kind=8), dimension(:), target, intent(inout) :: int_attr + character(len=*) :: name + character(len=*) :: unit + type(int_attr_ptr), allocatable, dimension(:) :: tmp + integer :: n_int_attr, n, j + + int_attr=0 ! set to zero when registring, for safety + + ! check if name is unique + do j = 1, this%attr_num + if(trim(this%attrib(j)%name) == trim(name)) then + print *, "Attribute name not unique" + stop + end if + end do + + do j = 1, this%int_attr_num + if(trim(this%int_attrib(j)%name) == trim(name)) then + print *, "Attribute name not unique" + stop + end if + end do + + this%int_attr_num = this%int_attr_num + 1 + + n_int_attr = size(this%int_attrib) + + if (.not. allocated(this%int_attrib)) then + allocate(this%int_attrib(1)) + this%int_attrib(1)%aptr => int_attr + this%int_attrib(1)%name = name + this%int_attrib(1)%unit = unit + else + allocate(tmp(n_int_attr+1)) + + do n = 1, n_int_attr + tmp(n)%aptr => this%int_attrib(n)%aptr + tmp(n)%name = this%int_attrib(n)%name + tmp(n)%unit = this%int_attrib(n)%unit + enddo + + tmp(n_int_attr+1)%aptr => int_attr + tmp(n_int_attr+1)%name = name + tmp(n_int_attr+1)%unit = unit + + deallocate(this%int_attrib) + + call move_alloc(from=tmp, to=this%int_attrib) + + endif + + end subroutine register_int_attribute + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine reset_int_attribute(this, int_attr, name) + class(base_parcel_type), intent(inout) :: this + integer(kind=8), dimension(:), target, intent(inout) :: int_attr + character(len=*) :: name + integer :: j + logical :: l_success + + l_success= .false. + + ! check if name is unique + do j = 1, this%int_attr_num + if(trim(this%int_attrib(j)%name) == trim(name)) then + this%int_attrib(j)%aptr => int_attr + l_success=.true. + exit + end if + end do + + if(.not. l_success) then + print *, "Attribute not found." + stop + end if + + end subroutine reset_int_attribute + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine dynamic_parcel_alloc(this, num) + class(dynamic_parcel_type), intent(inout) :: this + integer, intent(in) :: num + + call this%base_alloc(num) + + allocate(this%volume(num)) + call this%register_attribute(this%volume, "volume", "m^3") + + if ((this%n_pos < 2) .or. (this%n_pos > 3)) then + print *, "Only 2- or 3-dimensional." + stop + endif + + if (this%n_pos == 2) then + allocate(this%vorticity(1, num)) + call this%register_attribute(this%vorticity(1, :), "z_vorticity", "1/s") + allocate(this%delta_vor(1, num)) + call this%register_attribute(this%delta_vor(1, :), "z_vorticity_rk_tendency", "1/s") + elseif (this%n_pos == 3) then + allocate(this%vorticity(3, num)) + call this%register_attribute(this%vorticity(1, :), "x_vorticity", "1/s") + call this%register_attribute(this%vorticity(2, :), "y_vorticity", "1/s") + call this%register_attribute(this%vorticity(3, :), "z_vorticity", "1/s") + allocate(this%delta_vor(3, num)) + call this%register_attribute(this%delta_vor(1, :), "x_vorticity_rk_tendency", "1/s") + call this%register_attribute(this%delta_vor(2, :), "y_vorticity_rk_tendency", "1/s") + call this%register_attribute(this%delta_vor(3, :), "z_vorticity_rk_tendency", "1/s") + endif + + if (this%has_labels) then + allocate(this%label(num)) + allocate(this%dilution(num)) + call this%register_int_attribute(this%label, "label", "-") + call this%register_attribute(this%dilution, "dilution", "-") + endif + + end subroutine dynamic_parcel_alloc + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine dynamic_parcel_dealloc(this) + class(dynamic_parcel_type), intent(inout) :: this + + call try_deallocate(this%volume) + call try_deallocate(this%vorticity) + call try_deallocate(this%delta_vor) + call try_deallocate(this%dilution) + call try_deallocate(this%label) + call this%base_dealloc + + end subroutine dynamic_parcel_dealloc + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine dynamic_parcel_resize(this, new_size) + class(dynamic_parcel_type), intent(inout) :: this + integer, intent(in) :: new_size + + call this%base_resize(new_size) + + call resize_array(this%volume, new_size, this%local_num) + call resize_array(this%vorticity, new_size, this%local_num) + call resize_array(this%delta_vor, new_size, this%local_num) + + call this%reset_attribute(this%volume, "volume") + if (this%n_pos == 2) then + call this%reset_attribute(this%vorticity(1, :), "z_vorticity") + call this%reset_attribute(this%delta_vor(1, :), "z_vorticity_rk_tendency") + elseif (this%n_pos == 3) then + call this%reset_attribute(this%vorticity(1, :), "x_vorticity") + call this%reset_attribute(this%vorticity(2, :), "y_vorticity") + call this%reset_attribute(this%vorticity(3, :), "z_vorticity") + call this%reset_attribute(this%delta_vor(1, :), "x_vorticity_rk_tendency") + call this%reset_attribute(this%delta_vor(2, :), "y_vorticity_rk_tendency") + call this%reset_attribute(this%delta_vor(3, :), "z_vorticity_rk_tendency") + endif + + if (this%has_labels) then + call resize_array(this%label, new_size, this%local_num) + call resize_array(this%dilution, new_size, this%local_num) + call this%reset_int_attribute(this%label, "label") + call this%reset_attribute(this%dilution, "dilution") + endif + + end subroutine dynamic_parcel_resize + + + ! Serialize all parcel attributes into a single buffer + subroutine base_parcel_serialize(this, n, buffer) + class(base_parcel_type), intent(in) :: this + integer, intent(in) :: n + integer :: j + double precision, intent(out) :: buffer(this%attr_num+this%int_attr_num) + + do j = 1, this%attr_num + buffer(j)=this%attrib(j)%aptr(n) + end do + + do j = 1, this%int_attr_num + buffer(j+this%attr_num)=this%int_attrib(j)%aptr(n) + end do + + end subroutine base_parcel_serialize + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + ! Deserialize parcel attributes into a single buffer + subroutine base_parcel_deserialize(this, n, buffer) + class(base_parcel_type), intent(inout) :: this + integer, intent(in) :: n + integer :: j + double precision, intent(in) :: buffer(this%attr_num+this%int_attr_num) + + do j = 1, this%attr_num + this%attrib(j)%aptr(n)=buffer(j) + end do + + do j = 1, this%int_attr_num + this%int_attrib(j+this%attr_num)%aptr(n)=nint(buffer(j)) + end do + + end subroutine base_parcel_deserialize + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine parcel_pack(this, pid, num, buffer) + class(base_parcel_type), intent(in) :: this + integer, intent(in) :: pid(:) + integer, intent(in) :: num + double precision, intent(out) :: buffer(:) + integer :: n, i, j + + do n = 1, num + i = 1 + (n-1) * (this%attr_num+this%int_attr_num) + j = n * (this%attr_num+this%int_attr_num) + call this%serialize(pid(n), buffer(i:j)) + enddo + end subroutine parcel_pack + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine parcel_unpack(this, num, buffer) + class(base_parcel_type), intent(inout) :: this + integer, intent(in) :: num + double precision, intent(in) :: buffer(:) + integer :: n, i, j + + do n = 1, num + i = 1 + (n-1) * (this%attr_num+this%int_attr_num) + j = n * (this%attr_num+this%int_attr_num) + call this%deserialize(this%local_num + n, buffer(i:j)) + enddo + + this%local_num = this%local_num + num + + end subroutine parcel_unpack + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + ! This algorithm replaces invalid parcels with valid parcels + ! from the end of the container + ! @param[in] pid are the parcel indices of the parcels to be deleted + ! @param[in] n_del is the array size of pid + ! @pre + ! - pid must be sorted in ascending order + ! - pid must be contiguously filled + ! The above preconditions must be fulfilled so that the + ! parcel pack algorithm works correctly. + subroutine parcel_delete(this, pid, n_del) + class(base_parcel_type), intent(inout) :: this + integer, intent(in) :: pid(0:) + integer, intent(in) :: n_del + integer :: k, l, m + + ! l points always to the last valid parcel + l = this%local_num + + ! k points always to last invalid parcel in pid + k = n_del + + ! find last parcel which is not invalid + do while ((k > 0) .and. (l == pid(k))) + l = l - 1 + k = k - 1 + enddo + + if (l == -1) then + print *, "in parcel_container::parcel_delete: more than all parcels are invalid." + stop + endif + + ! replace invalid parcels with the last valid parcel + m = 1 + + do while (m <= k) + ! invalid parcel; overwrite *pid(m)* with last valid parcel *l* + call this%replace(pid(m), l) + + l = l - 1 + + ! find next valid last parcel + do while ((k > 0) .and. (l == pid(k))) + l = l - 1 + k = k - 1 + enddo + + ! next invalid + m = m + 1 + enddo + + ! update number of valid parcels + this%local_num = this%local_num - n_del + + end subroutine parcel_delete + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + ! Overwrite parcel n with parcel m. This subroutine only replaces the + ! common types. + ! ATTENTION: Extended types must replace additional parcel attributes + ! in their own routine. + ! @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 base_parcel_replace(this, n, m) + class(base_parcel_type), intent(inout) :: this + integer, intent(in) :: n, m + integer :: j + + + do j = 1, this%attr_num + this%attrib(j)%aptr(n) = this%attrib(j)%aptr(m) + end do + + do j = 1, this%int_attr_num + this%int_attrib(j)%aptr(n) = this%int_attrib(j)%aptr(m) + end do + + end subroutine base_parcel_replace + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + ! Print parcel attributes + subroutine print_me(this) + class(base_parcel_type), intent(inout) :: this + integer :: j + + do j = 1, this%attr_num + print *, this%attrib(j)%name + print *, this%attrib(j)%unit + print *, this%attrib(j)%aptr + end do + + do j = 1, this%int_attr_num + print *, this%int_attrib(j)%name + print *, this%int_attrib(j)%unit + print *, this%int_attrib(j)%aptr + end do + + end subroutine print_me + +end module diff --git a/src/parcel_types/parcel_ellipsoid.f90 b/src/parcel_types/parcel_ellipsoid.f90 new file mode 100644 index 000000000..806ab0a58 --- /dev/null +++ b/src/parcel_types/parcel_ellipsoid.f90 @@ -0,0 +1,163 @@ + module parcel_ellipsoid + + use parcel_container + implicit none + + ! Adding the ellipsoid geomerty to the dynamics + type, extends(dynamic_parcel_type) :: ellipsoid_parcel_type ! add procedures for ellipsoid below + double precision, allocatable, dimension(:,:) :: B + double precision, allocatable, dimension(:,:) :: delta_B + double precision, allocatable, dimension(:,:) :: strain + double precision, allocatable, dimension(:,:) :: Vetas + double precision, allocatable, dimension(:,:) :: Vtaus + double precision, allocatable, dimension(:,:) :: evec ! for now add evec for 2D. Need to think about consistency with 3D + character(len=16) :: shape_type ! (e.g. "ellipsoid5", for B with 5 elements, strain with 8) + logical :: par2grid_1p = .false. ! use one point for par2grid + logical :: grid2par_1p = .false. ! use one point for grid2par + + contains + procedure :: alloc => ellipsoid_parcel_alloc + procedure :: dealloc => ellipsoid_parcel_dealloc + procedure :: resize => ellipsoid_parcel_resize + + ! Other ellipsoid procedures to go here + end type + + contains + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine ellipsoid_parcel_alloc(this, num) + class(ellipsoid_parcel_type), intent(inout) :: this + integer, intent(in) :: num + + call this%dynamic_parcel_type%alloc(num) + + + if (this%shape_type == "ellipsoid2") then + allocate(this%B(2, num)) + allocate(this%delta_B(2, num)) + allocate(this%strain(4, num)) + allocate(this%evec(2, num)) + call this%register_attribute(this%B(1, :), "B11", "m^2") + call this%register_attribute(this%B(2, :), "B12", "m^2") + call this%register_attribute(this%delta_B(1, :), "B11_rk_tendency", "m^2/s") + call this%register_attribute(this%delta_B(2, :), "B12_rk_tendency", "m^2/s") + call this%register_attribute(this%strain(1, :), "DUDX", "1/s") + call this%register_attribute(this%strain(2, :), "DUDY", "1/s") + call this%register_attribute(this%strain(3, :), "DVDX", "1/s") + call this%register_attribute(this%strain(4, :), "DVDY", "1/s") + call this%register_attribute(this%evec(1, :), "evec 1", "m") + call this%register_attribute(this%evec(2, :), "evec 2", "m") + elseif (this%shape_type == "ellipsoid5") then + allocate(this%B(5, num)) + allocate(this%delta_B(5, num)) + allocate(this%strain(8, num)) + allocate(this%Vetas(3, num)) + allocate(this%Vtaus(3, num)) + call this%register_attribute(this%B(1, :), "B11", "m^2") + call this%register_attribute(this%B(2, :), "B12", "m^2") + call this%register_attribute(this%B(3, :), "B13", "m^2") + call this%register_attribute(this%B(4, :), "B22", "m^2") + call this%register_attribute(this%B(5, :), "B23", "m^2") + call this%register_attribute(this%delta_B(1, :), "B11_rk_tendency", "m^2/s") + call this%register_attribute(this%delta_B(2, :), "B12_rk_tendency", "m^2/s") + call this%register_attribute(this%delta_B(3, :), "B13_rk_tendency", "m^2/s") + call this%register_attribute(this%delta_B(4, :), "B22_rk_tendency", "m^2/s") + call this%register_attribute(this%delta_B(5, :), "B23_rk_tendency", "m^2/s") + call this%register_attribute(this%strain(1, :), "DUDX", "1/s") + call this%register_attribute(this%strain(2, :), "DUDY", "1/s") + call this%register_attribute(this%strain(3, :), "DUDZ", "1/s") + call this%register_attribute(this%strain(4, :), "DVDX", "1/s") + call this%register_attribute(this%strain(5, :), "DVDY", "1/s") + call this%register_attribute(this%strain(6, :), "DVDZ", "1/s") + call this%register_attribute(this%strain(7, :), "DWDX", "1/s") + call this%register_attribute(this%strain(8, :), "DWDY", "1/s") + call this%register_attribute(this%Vetas(1, :), "Veta1", "m") + call this%register_attribute(this%Vetas(2, :), "Veta2", "m") + call this%register_attribute(this%Vetas(3, :), "Veta3", "m") + call this%register_attribute(this%Vtaus(1, :), "Vtau1", "m") + call this%register_attribute(this%Vtaus(2, :), "Vtau2", "m") + call this%register_attribute(this%Vtaus(3, :), "Vtau3", "m") + else + print *, "Ellipsoid type not known." + stop + endif + + end subroutine ellipsoid_parcel_alloc + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine ellipsoid_parcel_dealloc(this) + class(ellipsoid_parcel_type), intent(inout) :: this + + call try_deallocate(this%B) + call try_deallocate(this%delta_B) + call try_deallocate(this%strain) + call try_deallocate(this%Vetas) + call try_deallocate(this%Vtaus) + call try_deallocate(this%evec) + + call this%dynamic_parcel_type%dealloc + + end subroutine ellipsoid_parcel_dealloc + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine ellipsoid_parcel_resize(this, new_size) + class(ellipsoid_parcel_type), intent(inout) :: this + integer, intent(in) :: new_size + + call this%dynamic_parcel_type%resize(new_size) + + call resize_array(this%B, new_size, this%local_num) + call resize_array(this%delta_B, new_size, this%local_num) + call resize_array(this%strain, new_size, this%local_num) + + if (this%shape_type == "ellipsoid2") then + call resize_array(this%evec, new_size, this%local_num) + call this%reset_attribute(this%B(1, :), "B11") + call this%reset_attribute(this%B(2, :), "B12") + call this%reset_attribute(this%delta_B(1, :), "B11_rk_tendency") + call this%reset_attribute(this%delta_B(2, :), "B12_rk_tendency") + call this%reset_attribute(this%strain(1, :), "DUDX") + call this%reset_attribute(this%strain(2, :), "DUDY") + call this%reset_attribute(this%strain(3, :), "DVDX") + call this%reset_attribute(this%strain(4, :), "DVDY") + call this%reset_attribute(this%evec(1, :), "evec 1") + call this%reset_attribute(this%evec(2, :), "evec 2") + elseif (this%shape_type == "ellipsoid5") then + call resize_array(this%Vetas, new_size, this%local_num) + call resize_array(this%Vtaus, new_size, this%local_num) + call this%reset_attribute(this%B(1, :), "B11") + call this%reset_attribute(this%B(2, :), "B12") + call this%reset_attribute(this%B(3, :), "B13") + call this%reset_attribute(this%B(4, :), "B22") + call this%reset_attribute(this%B(5, :), "B23") + call this%reset_attribute(this%delta_B(1, :), "B11_rk_tendency") + call this%reset_attribute(this%delta_B(2, :), "B12_rk_tendency") + call this%reset_attribute(this%delta_B(3, :), "B13_rk_tendency") + call this%reset_attribute(this%delta_B(4, :), "B22_rk_tendency") + call this%reset_attribute(this%delta_B(5, :), "B23_rk_tendency") + call this%reset_attribute(this%strain(1, :), "DUDX") + call this%reset_attribute(this%strain(2, :), "DUDY") + call this%reset_attribute(this%strain(3, :), "DUDZ") + call this%reset_attribute(this%strain(4, :), "DVDX") + call this%reset_attribute(this%strain(5, :), "DVDY") + call this%reset_attribute(this%strain(6, :), "DVDZ") + call this%reset_attribute(this%strain(7, :), "DWDX") + call this%reset_attribute(this%strain(8, :), "DWDY") + call this%reset_attribute(this%Vetas(1, :), "Veta1") + call this%reset_attribute(this%Vetas(2, :), "Veta2") + call this%reset_attribute(this%Vetas(3, :), "Veta3") + call this%reset_attribute(this%Vtaus(1, :), "Vtau1") + call this%reset_attribute(this%Vtaus(2, :), "Vtau2") + call this%reset_attribute(this%Vtaus(3, :), "Vtau3") + else + print *, "Ellipsoid type not known." + stop + endif + + end subroutine ellipsoid_parcel_resize + +end module diff --git a/src/parcel_types/parcel_types.f90 b/src/parcel_types/parcel_types.f90 new file mode 100644 index 000000000..76b236c87 --- /dev/null +++ b/src/parcel_types/parcel_types.f90 @@ -0,0 +1,207 @@ + module parcel_types + + use parcel_ellipsoid + implicit none + + type, extends(ellipsoid_parcel_type) :: idealised_parcel_type ! add procedures + double precision, allocatable, dimension(:) :: humidity + double precision, allocatable, dimension(:) :: buoyancy + logical :: is_moist = .false. + + contains + procedure :: alloc => idealised_parcel_alloc + procedure :: dealloc=> idealised_parcel_dealloc + procedure :: resize => idealised_parcel_resize + + ! get_buoyancy added here + end type + + type, extends(ellipsoid_parcel_type) :: realistic_parcel_type ! add procedures + double precision, allocatable, dimension(:) :: qv + double precision, allocatable, dimension(:) :: ql + double precision, allocatable, dimension(:) :: theta + double precision, allocatable, dimension(:) :: Nl ! optional droplet number + logical :: is_moist = .false. + logical :: has_droplets = .false. + + contains + procedure :: alloc => realistic_parcel_alloc + procedure :: dealloc=> realistic_parcel_dealloc + procedure :: resize => realistic_parcel_resize + + ! get_buoyancy added here + end type + + type, extends(base_parcel_type) :: prec_parcel_type ! add procedures + double precision, allocatable, dimension(:) :: volume + double precision, allocatable, dimension(:) :: qr + double precision, allocatable, dimension(:) :: Nr ! droplet number + + contains + procedure :: alloc => prec_parcel_alloc + procedure :: dealloc=> prec_parcel_dealloc + procedure :: resize => prec_parcel_resize + + ! get_buoyancy added here + end type + + contains + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine idealised_parcel_alloc(this, num) + class(idealised_parcel_type), intent(inout) :: this + integer, intent(in) :: num + + call this%ellipsoid_parcel_type%alloc(num) + + allocate(this%humidity(num)) + allocate(this%buoyancy(num)) + + call this%register_attribute(this%humidity, "humidity", "kg/kg") + call this%register_attribute(this%buoyancy, "buoyancy", "m/s^2") + + end subroutine idealised_parcel_alloc + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine idealised_parcel_dealloc(this) + class(idealised_parcel_type), intent(inout) :: this + + call try_deallocate(this%humidity) + call try_deallocate(this%buoyancy) + call this%ellipsoid_parcel_type%dealloc + + end subroutine idealised_parcel_dealloc + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine idealised_parcel_resize(this, new_size) + class(idealised_parcel_type), intent(inout) :: this + integer, intent(in) :: new_size + + call this%ellipsoid_parcel_type%resize(new_size) + + call resize_array(this%buoyancy, new_size, this%local_num) + call this%reset_attribute(this%buoyancy, "buoyancy") + + if(this%is_moist) then + call resize_array(this%humidity, new_size, this%local_num) + call this%reset_attribute(this%humidity, "humidity") + endif + + end subroutine idealised_parcel_resize + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine realistic_parcel_alloc(this, num) + class(realistic_parcel_type), intent(inout) :: this + integer, intent(in) :: num + + call this%ellipsoid_parcel_type%alloc(num) + + allocate(this%theta(num)) + + call this%register_attribute(this%theta, "theta", "K") + + if(this%is_moist) then + allocate(this%qv(num)) + allocate(this%ql(num)) + call this%register_attribute(this%qv, "qv", "kg/kg") + call this%register_attribute(this%ql, "ql", "kg/kg") + endif + + if(this%has_droplets) then + allocate(this%Nl(num)) + call this%register_attribute(this%Nl, "Nl", "/m^3") + endif + + end subroutine realistic_parcel_alloc + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine realistic_parcel_dealloc(this) + class(realistic_parcel_type), intent(inout) :: this + + call try_deallocate(this%theta) + call try_deallocate(this%qv) + call try_deallocate(this%ql) + call try_deallocate(this%Nl) + + call this%ellipsoid_parcel_type%dealloc + + end subroutine realistic_parcel_dealloc + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine realistic_parcel_resize(this, new_size) + class(realistic_parcel_type), intent(inout) :: this + integer, intent(in) :: new_size + + call this%ellipsoid_parcel_type%resize(new_size) + + call resize_array(this%theta, new_size, this%local_num) + call resize_array(this%qv, new_size, this%local_num) + call resize_array(this%ql, new_size, this%local_num) + + call this%reset_attribute(this%theta, "theta") + call this%reset_attribute(this%qv, "qv") + call this%reset_attribute(this%ql, "ql") + + if(this%has_droplets) then + call resize_array(this%Nl, new_size, this%local_num) + call this%reset_attribute(this%Nl, "Nl") + endif + + end subroutine realistic_parcel_resize + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine prec_parcel_alloc(this, num) + class(prec_parcel_type), intent(inout) :: this + integer, intent(in) :: num + + call this%base_alloc(num) + + allocate(this%volume(num)) + allocate(this%qr(num)) + allocate(this%Nr(num)) + + call this%register_attribute(this%volume, "volume", "m^3") + call this%register_attribute(this%qr, "qr", "kg/kg") + call this%register_attribute(this%Nr, "Nr", "/m^3") + + end subroutine prec_parcel_alloc + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine prec_parcel_dealloc(this) + class(prec_parcel_type), intent(inout) :: this + + call try_deallocate(this%volume) + call try_deallocate(this%qr) + call try_deallocate(this%Nr) + + call this%base_dealloc + + end subroutine prec_parcel_dealloc + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + subroutine prec_parcel_resize(this, new_size) + class(prec_parcel_type), intent(inout) :: this + integer, intent(in) :: new_size + + call this%base_resize(new_size) + + call resize_array(this%volume, new_size, this%local_num) + call resize_array(this%qr, new_size, this%local_num) + call resize_array(this%Nr, new_size, this%local_num) + + call this%reset_attribute(this%volume, "volume") + call this%reset_attribute(this%qr, "qr") + call this%reset_attribute(this%Nr, "Nr") + + end subroutine prec_parcel_resize + +end module From 425c353ffc1c65971eccf9a8fac4762a1f931407 Mon Sep 17 00:00:00 2001 From: Steven Boeing Date: Tue, 24 Jun 2025 21:21:41 +0100 Subject: [PATCH 4/7] Work to generalise dimensions --- src/2d/epic2d.f90 | 2 +- src/2d/utils/utils.f90 | 4 +- src/parcel_types/parcel_container.f90 | 155 +++++++++++++------- src/parcel_types/parcel_ellipsoid.f90 | 197 ++++++++++++++++---------- src/parcel_types/parcel_types.f90 | 52 +++++-- 5 files changed, 273 insertions(+), 137 deletions(-) diff --git a/src/2d/epic2d.f90 b/src/2d/epic2d.f90 index 9dc3121d4..7a940250f 100644 --- a/src/2d/epic2d.f90 +++ b/src/2d/epic2d.f90 @@ -4,7 +4,7 @@ program epic2d use constants, only : zero use timer - use dynamic_parcels + use dynamic_parcels, only : parcels use parcel_bc use parcel_split, only : split_ellipses, split_timer use parcel_merge, only : merge_ellipses, merge_timer diff --git a/src/2d/utils/utils.f90 b/src/2d/utils/utils.f90 index 0a046371c..b3ee1847b 100644 --- a/src/2d/utils/utils.f90 +++ b/src/2d/utils/utils.f90 @@ -214,9 +214,9 @@ end subroutine setup_domain_and_parameters subroutine setup_parcels character(len=16) :: file_type - call parcels%set_dimension(2) + parcels%dim_string='xz' + call parcels%ellipsoid_dimensions() parcels%is_moist=.true. - parcels%shape_type="ellipsoid2" call parcels%alloc(max_num_parcels) if (l_restart) then diff --git a/src/parcel_types/parcel_container.f90 b/src/parcel_types/parcel_container.f90 index 70bc72828..8d76ec1e0 100644 --- a/src/parcel_types/parcel_container.f90 +++ b/src/parcel_types/parcel_container.f90 @@ -39,7 +39,9 @@ module parcel_container integer :: local_num ! local number of parcels integer(kind=int64) :: total_num ! global number of parcels (over all MPI ranks) integer :: max_num ! capacity per attribute, i.e. maximum number of parcels - integer, private :: n_pos = -1 ! number of spatial dimensions + integer :: n_pos = -1 ! number of spatial dimensions + character(len=1), allocatable, dimension(:) :: pos_names ! Names of directions for positions + character(len=8) :: dim_string ! Names of directions for positions contains procedure :: base_alloc => base_parcel_alloc @@ -55,27 +57,30 @@ module parcel_container procedure(parcel_dealloc), deferred :: dealloc procedure(parcel_resize), deferred :: resize procedure :: print_me - procedure :: set_dimension + procedure :: set_dimensions procedure :: register_attribute procedure :: register_int_attribute procedure :: reset_attribute procedure :: reset_int_attribute - end type ! This type is for parcels associated with dynamics (which will always have volume and voriticity in EPIC) - type, extends(base_parcel_type) :: dynamic_parcel_type ! add procedures + type, abstract, extends(base_parcel_type) :: dynamic_parcel_type ! add procedures double precision, allocatable, dimension(:) :: volume double precision, allocatable, dimension(:,:) :: vorticity double precision, allocatable, dimension(:,:) :: delta_vor double precision, allocatable, dimension(:) :: dilution integer(kind=8), allocatable, dimension(:) :: label + integer, private :: n_vor = -1 ! number of voriticity components + character(len=1), allocatable, dimension(:) :: vor_names ! Names of vorticity components logical :: has_labels = .false. contains - procedure :: alloc => dynamic_parcel_alloc - procedure :: dealloc => dynamic_parcel_dealloc - procedure :: resize => dynamic_parcel_resize + procedure :: dynamic_alloc => dynamic_parcel_alloc + procedure :: dynamic_dealloc => dynamic_parcel_dealloc + procedure :: dynamic_resize => dynamic_parcel_resize + procedure :: set_vorticity_dimensions + procedure(dynamic_parcel_get_buoyancy), deferred :: get_buoyancy end type @@ -109,8 +114,18 @@ end function parcel_is_small module procedure :: try_deallocate_1d module procedure :: try_deallocate_1d_integer module procedure :: try_deallocate_2d + module procedure :: try_deallocate_string end interface try_deallocate + interface + subroutine dynamic_parcel_get_buoyancy(this, num, buoyancy) + import dynamic_parcel_type + class(dynamic_parcel_type), intent(inout) :: this + integer, intent(in) :: num + double precision, intent(out) :: buoyancy + end subroutine dynamic_parcel_get_buoyancy + end interface + contains !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -148,6 +163,17 @@ end subroutine try_deallocate_2d !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + subroutine try_deallocate_string(in_array) + character(len=*), allocatable, dimension(:) :: in_array + + if (allocated(in_array)) then + deallocate(in_array) + endif + + end subroutine try_deallocate_string + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ! Allocate parcel memory ! ATTENTION: Extended types must allocate additional parcel attributes ! in their own routine. @@ -155,7 +181,6 @@ end subroutine try_deallocate_2d subroutine base_parcel_alloc(this, num) class(base_parcel_type), intent(inout), target :: this integer, intent(in) :: num - character(len=1) :: dir(3) integer :: n this%max_num = num @@ -168,14 +193,12 @@ subroutine base_parcel_alloc(this, num) stop endif - dir = (/'x', 'y', 'z'/) - allocate(this%position(this%n_pos, num)) allocate(this%delta_pos(this%n_pos, num)) do n = 1, this%n_pos - call this%register_attribute(this%position(n, :), dir(n) // "_position", "m") - call this%register_attribute(this%delta_pos(n, :), dir(n) // "_position_rk_tendency", "m/s") + call this%register_attribute(this%position(n, :), this%pos_names(n) // "_position", "m") + call this%register_attribute(this%delta_pos(n, :), this%pos_names(n) // "_position_rk_tendency", "m/s") enddo end subroutine base_parcel_alloc @@ -196,6 +219,7 @@ subroutine base_parcel_dealloc(this) call try_deallocate(this%position) call try_deallocate(this%delta_pos) + call try_deallocate(this%pos_names) if (allocated(this%attrib)) then deallocate(this%attrib) @@ -212,7 +236,6 @@ end subroutine base_parcel_dealloc subroutine base_parcel_resize(this, new_size) class(base_parcel_type), intent(inout), target :: this integer, intent(in) :: new_size - character(len=1) :: dir(3) integer :: n if (new_size < this%local_num) then @@ -222,30 +245,28 @@ subroutine base_parcel_resize(this, new_size) this%max_num = new_size - if (this%n_pos > 3) then - print *, "Only 3 dimensions allowed." + if ((this%n_pos < 1) .or. (this%n_pos > 3)) then + print *, "ERROR: base_parcel_resize:: only 1-, 2- or 3-dimensional parcels allowed." stop endif - dir = (/'x', 'y', 'z'/) - call resize_array(this%position, new_size, this%local_num) call resize_array(this%delta_pos, new_size, this%local_num) do n = 1, this%n_pos - call this%reset_attribute(this%position(n, :), dir(n) // "_position") - call this%reset_attribute(this%delta_pos(n, :), dir(n) // "_position_rk_tendency") + call this%reset_attribute(this%position(n, :), this%pos_names(n) // "_position") + call this%reset_attribute(this%delta_pos(n, :), this%pos_names(n) // "_position_rk_tendency") enddo end subroutine base_parcel_resize !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - subroutine set_dimension(this, n_dim) + subroutine set_dimensions(this) class(base_parcel_type), intent(inout) :: this - integer, intent(in) :: n_dim + integer :: n_dim, i_dim - this%n_pos = n_dim + n_dim = len_trim(this%dim_string) if (this%n_pos /= -1) then print *, "WARNING: Dimension already set." @@ -253,11 +274,50 @@ subroutine set_dimension(this, n_dim) endif if ((n_dim < 1) .or. (n_dim > 3)) then - print *, "Only 1-, 2- or 3-dimensional." + print *, "ERROR: set_dimension:: only 1-, 2- or 3-dimensional parcels allowed." stop endif - end subroutine set_dimension + this%n_pos=n_dim + + allocate(this%pos_names(this%n_pos)) + + do i_dim = 1,this%n_pos + this%pos_names(i_dim)=this%dim_string(i_dim:i_dim) + enddo + + end subroutine set_dimensions + + subroutine set_vorticity_dimensions(this) + class(dynamic_parcel_type), intent(inout) :: this + integer :: i_dim + + call set_dimensions(this) + + if (trim(this%dim_string) == trim('xy')) then + this%n_vor=1 + allocate(this%vor_names(this%n_vor)) + this%vor_names(1)='z' + elseif (trim(this%dim_string) == trim('xz')) then + this%n_vor=1 + allocate(this%vor_names(this%n_vor)) + this%vor_names(1)='y' + elseif(trim(this%dim_string) == trim('yz')) then + this%n_vor=1 + allocate(this%vor_names(this%n_vor)) + this%vor_names(1)='x' + elseif(trim(this%dim_string) == trim('xyz')) then + this%n_vor=3 + allocate(this%vor_names(this%n_vor)) + do i_dim = 1,this%n_vor + this%vor_names(i_dim)=this%dim_string(i_dim:i_dim) + enddo + else + print *, "ERROR: set_vorticiy_dimension:: only xy, xz, yz or xyz allowed." + stop + endif + + end subroutine set_vorticity_dimensions !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -275,6 +335,7 @@ subroutine register_attribute(this, attr, name, unit) do j = 1, this%attr_num if(trim(this%attrib(j)%name) == trim(name)) then print *, "Attribute name not unique" + print *, name stop end if end do @@ -282,6 +343,7 @@ subroutine register_attribute(this, attr, name, unit) do j = 1, this%int_attr_num if(trim(this%int_attrib(j)%name) == trim(name)) then print *, "Attribute name not unique" + print *, name stop end if end do @@ -432,32 +494,24 @@ end subroutine reset_int_attribute subroutine dynamic_parcel_alloc(this, num) class(dynamic_parcel_type), intent(inout) :: this integer, intent(in) :: num + integer :: i_dim call this%base_alloc(num) allocate(this%volume(num)) call this%register_attribute(this%volume, "volume", "m^3") - if ((this%n_pos < 2) .or. (this%n_pos > 3)) then - print *, "Only 2- or 3-dimensional." + if(.not. ((this%n_vor == 1) .or. (this%n_vor == 3))) then + print *, "Vorticity needs 1 or 3 components." stop endif - if (this%n_pos == 2) then - allocate(this%vorticity(1, num)) - call this%register_attribute(this%vorticity(1, :), "z_vorticity", "1/s") - allocate(this%delta_vor(1, num)) - call this%register_attribute(this%delta_vor(1, :), "z_vorticity_rk_tendency", "1/s") - elseif (this%n_pos == 3) then - allocate(this%vorticity(3, num)) - call this%register_attribute(this%vorticity(1, :), "x_vorticity", "1/s") - call this%register_attribute(this%vorticity(2, :), "y_vorticity", "1/s") - call this%register_attribute(this%vorticity(3, :), "z_vorticity", "1/s") - allocate(this%delta_vor(3, num)) - call this%register_attribute(this%delta_vor(1, :), "x_vorticity_rk_tendency", "1/s") - call this%register_attribute(this%delta_vor(2, :), "y_vorticity_rk_tendency", "1/s") - call this%register_attribute(this%delta_vor(3, :), "z_vorticity_rk_tendency", "1/s") - endif + allocate(this%vorticity(this%n_vor, num)) + allocate(this%delta_vor(this%n_vor, num)) + do i_dim=1,this%n_vor + call this%register_attribute(this%vorticity(i_dim, :), this%vor_names(i_dim)//"_vorticity", "1/s") + call this%register_attribute(this%delta_vor(i_dim, :), this%vor_names(i_dim)//"_vorticity_rk_tendency", "1/s") + enddo if (this%has_labels) then allocate(this%label(num)) @@ -478,6 +532,7 @@ subroutine dynamic_parcel_dealloc(this) call try_deallocate(this%delta_vor) call try_deallocate(this%dilution) call try_deallocate(this%label) + call try_deallocate(this%vor_names) call this%base_dealloc end subroutine dynamic_parcel_dealloc @@ -487,6 +542,8 @@ end subroutine dynamic_parcel_dealloc subroutine dynamic_parcel_resize(this, new_size) class(dynamic_parcel_type), intent(inout) :: this integer, intent(in) :: new_size + integer :: i_dim + call this%base_resize(new_size) @@ -495,17 +552,11 @@ subroutine dynamic_parcel_resize(this, new_size) call resize_array(this%delta_vor, new_size, this%local_num) call this%reset_attribute(this%volume, "volume") - if (this%n_pos == 2) then - call this%reset_attribute(this%vorticity(1, :), "z_vorticity") - call this%reset_attribute(this%delta_vor(1, :), "z_vorticity_rk_tendency") - elseif (this%n_pos == 3) then - call this%reset_attribute(this%vorticity(1, :), "x_vorticity") - call this%reset_attribute(this%vorticity(2, :), "y_vorticity") - call this%reset_attribute(this%vorticity(3, :), "z_vorticity") - call this%reset_attribute(this%delta_vor(1, :), "x_vorticity_rk_tendency") - call this%reset_attribute(this%delta_vor(2, :), "y_vorticity_rk_tendency") - call this%reset_attribute(this%delta_vor(3, :), "z_vorticity_rk_tendency") - endif + + do i_dim=1,this%n_vor + call this%reset_attribute(this%vorticity(i_dim, :), this%vor_names(i_dim)//"_vorticity") + call this%reset_attribute(this%delta_vor(i_dim, :), this%vor_names(i_dim)//"_vorticity_rk_tendency") + end do if (this%has_labels) then call resize_array(this%label, new_size, this%local_num) diff --git a/src/parcel_types/parcel_ellipsoid.f90 b/src/parcel_types/parcel_ellipsoid.f90 index 806ab0a58..cfb2e13cb 100644 --- a/src/parcel_types/parcel_ellipsoid.f90 +++ b/src/parcel_types/parcel_ellipsoid.f90 @@ -4,75 +4,133 @@ module parcel_ellipsoid implicit none ! Adding the ellipsoid geomerty to the dynamics - type, extends(dynamic_parcel_type) :: ellipsoid_parcel_type ! add procedures for ellipsoid below + type, abstract, extends(dynamic_parcel_type) :: ellipsoid_parcel_type ! add procedures for ellipsoid below double precision, allocatable, dimension(:,:) :: B double precision, allocatable, dimension(:,:) :: delta_B double precision, allocatable, dimension(:,:) :: strain double precision, allocatable, dimension(:,:) :: Vetas double precision, allocatable, dimension(:,:) :: Vtaus double precision, allocatable, dimension(:,:) :: evec ! for now add evec for 2D. Need to think about consistency with 3D - character(len=16) :: shape_type ! (e.g. "ellipsoid5", for B with 5 elements, strain with 8) logical :: par2grid_1p = .false. ! use one point for par2grid logical :: grid2par_1p = .false. ! use one point for grid2par - + integer, private :: n_strain = -1 ! number of strain components + character(len=4), allocatable, dimension(:) :: strain_names ! Names of directions for strain + integer, private :: n_shape = -1 ! number of shape components + character(len=3), allocatable, dimension(:) :: shape_names ! Names of shape components contains - procedure :: alloc => ellipsoid_parcel_alloc - procedure :: dealloc => ellipsoid_parcel_dealloc - procedure :: resize => ellipsoid_parcel_resize + procedure :: ellipsoid_alloc => ellipsoid_parcel_alloc + procedure :: ellipsoid_dealloc => ellipsoid_parcel_dealloc + procedure :: ellipsoid_resize => ellipsoid_parcel_resize + procedure :: ellipsoid_split => ellipsoid_parcel_split + procedure :: ellipsoid_dimensions => set_ellipsoid_dimensions + procedure(parcel_split), deferred :: split ! Other ellipsoid procedures to go here end type + interface + subroutine parcel_split(this, n, n_thread_loc) + import ellipsoid_parcel_type + class(ellipsoid_parcel_type), intent(inout) :: this + integer, intent(in) :: n + integer, intent(in) :: n_thread_loc + end subroutine parcel_split + end interface + contains + subroutine set_ellipsoid_dimensions(this) + class(ellipsoid_parcel_type), intent(inout) :: this + integer :: i_dim + + call this%set_vorticity_dimensions + + if (trim(this%dim_string) == 'xy') then + this%n_strain=4 + allocate(this%strain_names(this%n_strain)) + this%strain_names(1)='DUDX' + this%strain_names(2)='DUDY' + this%strain_names(3)='DVDX' + this%strain_names(4)='DVDY' + elseif (trim(this%dim_string) == 'xz') then + this%n_strain=4 + allocate(this%strain_names(this%n_strain)) + this%strain_names(1)='DUDX' + this%strain_names(2)='DUDZ' + this%strain_names(3)='DWDX' + this%strain_names(4)='DWDZ' + elseif(trim(this%dim_string) == 'yz') then + this%n_strain=4 + allocate(this%strain_names(this%n_strain)) + this%strain_names(1)='DVDY' + this%strain_names(2)='DVDZ' + this%strain_names(3)='DWDY' + this%strain_names(4)='DWDZ' + elseif(trim(this%dim_string) == 'xyz') then + this%n_strain=9 + allocate(this%strain_names(this%n_strain)) + this%strain_names(1)='DUDX' + this%strain_names(2)='DUDY' + this%strain_names(3)='DUDZ' + this%strain_names(4)='DVDX' + this%strain_names(5)='DVDY' + this%strain_names(6)='DVDZ' + this%strain_names(7)='DWDX' + this%strain_names(8)='DWDY' + this%strain_names(9)='DWDZ' + else + print *, "ERROR: set_ellipsoid_dimension:: only xy, xz, yz or xyz allowed." + stop + endif + + if (len_trim(this%dim_string) == 2) then + this%n_shape=2 + allocate(this%shape_names(this%n_shape)) + this%shape_names(1)='B11' + this%shape_names(2)='B12' + elseif (len_trim(this%dim_string) == 3) then + this%n_shape=5 + allocate(this%shape_names(this%n_shape)) + this%shape_names(1)='B11' + this%shape_names(2)='B12' + this%shape_names(3)='B13' + this%shape_names(4)='B22' + this%shape_names(5)='B23' + else + print *, "ERROR: set_ellipsoid_dimension:: wrong number of dimensions" + stop + endif + + end subroutine set_ellipsoid_dimensions !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: subroutine ellipsoid_parcel_alloc(this, num) class(ellipsoid_parcel_type), intent(inout) :: this integer, intent(in) :: num + integer :: i_dim + + call this%dynamic_alloc(num) + + allocate(this%B(this%n_shape, num)) + allocate(this%delta_B(this%n_shape, num)) + allocate(this%strain(this%n_strain, num)) - call this%dynamic_parcel_type%alloc(num) + do i_dim=1, this%n_shape + call this%register_attribute(this%B(i_dim, :), this%shape_names(i_dim), "1/s") + call this%register_attribute(this%delta_B(i_dim, :), this%shape_names(i_dim)//"_rk_tendency", "1/s") + enddo + do i_dim=1, this%n_strain + call this%register_attribute(this%strain(i_dim, :), this%strain_names(i_dim), "1/s") + enddo - if (this%shape_type == "ellipsoid2") then - allocate(this%B(2, num)) - allocate(this%delta_B(2, num)) - allocate(this%strain(4, num)) + if (this%n_pos == 2) then allocate(this%evec(2, num)) - call this%register_attribute(this%B(1, :), "B11", "m^2") - call this%register_attribute(this%B(2, :), "B12", "m^2") - call this%register_attribute(this%delta_B(1, :), "B11_rk_tendency", "m^2/s") - call this%register_attribute(this%delta_B(2, :), "B12_rk_tendency", "m^2/s") - call this%register_attribute(this%strain(1, :), "DUDX", "1/s") - call this%register_attribute(this%strain(2, :), "DUDY", "1/s") - call this%register_attribute(this%strain(3, :), "DVDX", "1/s") - call this%register_attribute(this%strain(4, :), "DVDY", "1/s") call this%register_attribute(this%evec(1, :), "evec 1", "m") call this%register_attribute(this%evec(2, :), "evec 2", "m") - elseif (this%shape_type == "ellipsoid5") then - allocate(this%B(5, num)) - allocate(this%delta_B(5, num)) - allocate(this%strain(8, num)) + elseif (this%n_pos == 3) then allocate(this%Vetas(3, num)) allocate(this%Vtaus(3, num)) - call this%register_attribute(this%B(1, :), "B11", "m^2") - call this%register_attribute(this%B(2, :), "B12", "m^2") - call this%register_attribute(this%B(3, :), "B13", "m^2") - call this%register_attribute(this%B(4, :), "B22", "m^2") - call this%register_attribute(this%B(5, :), "B23", "m^2") - call this%register_attribute(this%delta_B(1, :), "B11_rk_tendency", "m^2/s") - call this%register_attribute(this%delta_B(2, :), "B12_rk_tendency", "m^2/s") - call this%register_attribute(this%delta_B(3, :), "B13_rk_tendency", "m^2/s") - call this%register_attribute(this%delta_B(4, :), "B22_rk_tendency", "m^2/s") - call this%register_attribute(this%delta_B(5, :), "B23_rk_tendency", "m^2/s") - call this%register_attribute(this%strain(1, :), "DUDX", "1/s") - call this%register_attribute(this%strain(2, :), "DUDY", "1/s") - call this%register_attribute(this%strain(3, :), "DUDZ", "1/s") - call this%register_attribute(this%strain(4, :), "DVDX", "1/s") - call this%register_attribute(this%strain(5, :), "DVDY", "1/s") - call this%register_attribute(this%strain(6, :), "DVDZ", "1/s") - call this%register_attribute(this%strain(7, :), "DWDX", "1/s") - call this%register_attribute(this%strain(8, :), "DWDY", "1/s") call this%register_attribute(this%Vetas(1, :), "Veta1", "m") call this%register_attribute(this%Vetas(2, :), "Veta2", "m") call this%register_attribute(this%Vetas(3, :), "Veta3", "m") @@ -80,7 +138,7 @@ subroutine ellipsoid_parcel_alloc(this, num) call this%register_attribute(this%Vtaus(2, :), "Vtau2", "m") call this%register_attribute(this%Vtaus(3, :), "Vtau3", "m") else - print *, "Ellipsoid type not known." + print *, "ERROR: ellipsoid_parcel_alloc:: n_pos must be 2 or 3." stop endif @@ -97,8 +155,10 @@ subroutine ellipsoid_parcel_dealloc(this) call try_deallocate(this%Vetas) call try_deallocate(this%Vtaus) call try_deallocate(this%evec) + call try_deallocate(this%shape_names) + call try_deallocate(this%strain_names) - call this%dynamic_parcel_type%dealloc + call this%dynamic_dealloc end subroutine ellipsoid_parcel_dealloc @@ -107,46 +167,31 @@ end subroutine ellipsoid_parcel_dealloc subroutine ellipsoid_parcel_resize(this, new_size) class(ellipsoid_parcel_type), intent(inout) :: this integer, intent(in) :: new_size + integer :: i_dim - call this%dynamic_parcel_type%resize(new_size) + call this%dynamic_resize(new_size) call resize_array(this%B, new_size, this%local_num) call resize_array(this%delta_B, new_size, this%local_num) call resize_array(this%strain, new_size, this%local_num) - if (this%shape_type == "ellipsoid2") then + do i_dim=1,this%n_shape + call this%reset_attribute(this%B(i_dim, :), this%shape_names(i_dim)) + call this%reset_attribute(this%delta_B(i_dim, :), this%shape_names(i_dim)//"_rk_tendency") + end do + + do i_dim=1,this%n_strain + call this%reset_attribute(this%strain(i_dim, :), this%strain_names(i_dim)) + end do + + ! Distinguish between 2D and 3D ellipsoids + if (this%n_pos == 2) then call resize_array(this%evec, new_size, this%local_num) - call this%reset_attribute(this%B(1, :), "B11") - call this%reset_attribute(this%B(2, :), "B12") - call this%reset_attribute(this%delta_B(1, :), "B11_rk_tendency") - call this%reset_attribute(this%delta_B(2, :), "B12_rk_tendency") - call this%reset_attribute(this%strain(1, :), "DUDX") - call this%reset_attribute(this%strain(2, :), "DUDY") - call this%reset_attribute(this%strain(3, :), "DVDX") - call this%reset_attribute(this%strain(4, :), "DVDY") call this%reset_attribute(this%evec(1, :), "evec 1") call this%reset_attribute(this%evec(2, :), "evec 2") - elseif (this%shape_type == "ellipsoid5") then + elseif (this%n_pos == 3) then call resize_array(this%Vetas, new_size, this%local_num) call resize_array(this%Vtaus, new_size, this%local_num) - call this%reset_attribute(this%B(1, :), "B11") - call this%reset_attribute(this%B(2, :), "B12") - call this%reset_attribute(this%B(3, :), "B13") - call this%reset_attribute(this%B(4, :), "B22") - call this%reset_attribute(this%B(5, :), "B23") - call this%reset_attribute(this%delta_B(1, :), "B11_rk_tendency") - call this%reset_attribute(this%delta_B(2, :), "B12_rk_tendency") - call this%reset_attribute(this%delta_B(3, :), "B13_rk_tendency") - call this%reset_attribute(this%delta_B(4, :), "B22_rk_tendency") - call this%reset_attribute(this%delta_B(5, :), "B23_rk_tendency") - call this%reset_attribute(this%strain(1, :), "DUDX") - call this%reset_attribute(this%strain(2, :), "DUDY") - call this%reset_attribute(this%strain(3, :), "DUDZ") - call this%reset_attribute(this%strain(4, :), "DVDX") - call this%reset_attribute(this%strain(5, :), "DVDY") - call this%reset_attribute(this%strain(6, :), "DVDZ") - call this%reset_attribute(this%strain(7, :), "DWDX") - call this%reset_attribute(this%strain(8, :), "DWDY") call this%reset_attribute(this%Vetas(1, :), "Veta1") call this%reset_attribute(this%Vetas(2, :), "Veta2") call this%reset_attribute(this%Vetas(3, :), "Veta3") @@ -154,10 +199,16 @@ subroutine ellipsoid_parcel_resize(this, new_size) call this%reset_attribute(this%Vtaus(2, :), "Vtau2") call this%reset_attribute(this%Vtaus(3, :), "Vtau3") else - print *, "Ellipsoid type not known." + print *, "ERROR: ellipsoid_parcel_resize:: n_pos must be 2 or 3." stop endif end subroutine ellipsoid_parcel_resize + subroutine ellipsoid_parcel_split(this, n, n_thread_loc) + class(ellipsoid_parcel_type), intent(inout) :: this + integer, intent(in) :: n + integer, intent(in) :: n_thread_loc + end subroutine ellipsoid_parcel_split + end module diff --git a/src/parcel_types/parcel_types.f90 b/src/parcel_types/parcel_types.f90 index 76b236c87..da04f3a90 100644 --- a/src/parcel_types/parcel_types.f90 +++ b/src/parcel_types/parcel_types.f90 @@ -12,8 +12,9 @@ module parcel_types procedure :: alloc => idealised_parcel_alloc procedure :: dealloc=> idealised_parcel_dealloc procedure :: resize => idealised_parcel_resize + procedure :: split => idealised_parcel_split + procedure :: get_buoyancy => idealised_parcel_get_buoyancy - ! get_buoyancy added here end type type, extends(ellipsoid_parcel_type) :: realistic_parcel_type ! add procedures @@ -28,6 +29,8 @@ module parcel_types procedure :: alloc => realistic_parcel_alloc procedure :: dealloc=> realistic_parcel_dealloc procedure :: resize => realistic_parcel_resize + procedure :: split => realistic_parcel_split + procedure :: get_buoyancy => realistic_parcel_get_buoyancy ! get_buoyancy added here end type @@ -53,14 +56,17 @@ subroutine idealised_parcel_alloc(this, num) class(idealised_parcel_type), intent(inout) :: this integer, intent(in) :: num - call this%ellipsoid_parcel_type%alloc(num) + call this%ellipsoid_alloc(num) - allocate(this%humidity(num)) allocate(this%buoyancy(num)) - call this%register_attribute(this%humidity, "humidity", "kg/kg") call this%register_attribute(this%buoyancy, "buoyancy", "m/s^2") + if(this%is_moist) then + allocate(this%humidity(num)) + call this%register_attribute(this%humidity, "humidity", "kg/kg") + endif + end subroutine idealised_parcel_alloc !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: @@ -70,7 +76,7 @@ subroutine idealised_parcel_dealloc(this) call try_deallocate(this%humidity) call try_deallocate(this%buoyancy) - call this%ellipsoid_parcel_type%dealloc + call this%ellipsoid_dealloc end subroutine idealised_parcel_dealloc @@ -80,7 +86,7 @@ subroutine idealised_parcel_resize(this, new_size) class(idealised_parcel_type), intent(inout) :: this integer, intent(in) :: new_size - call this%ellipsoid_parcel_type%resize(new_size) + call this%ellipsoid_resize(new_size) call resize_array(this%buoyancy, new_size, this%local_num) call this%reset_attribute(this%buoyancy, "buoyancy") @@ -98,7 +104,7 @@ subroutine realistic_parcel_alloc(this, num) class(realistic_parcel_type), intent(inout) :: this integer, intent(in) :: num - call this%ellipsoid_parcel_type%alloc(num) + call this%ellipsoid_alloc(num) allocate(this%theta(num)) @@ -128,7 +134,7 @@ subroutine realistic_parcel_dealloc(this) call try_deallocate(this%ql) call try_deallocate(this%Nl) - call this%ellipsoid_parcel_type%dealloc + call this%ellipsoid_dealloc end subroutine realistic_parcel_dealloc @@ -138,7 +144,7 @@ subroutine realistic_parcel_resize(this, new_size) class(realistic_parcel_type), intent(inout) :: this integer, intent(in) :: new_size - call this%ellipsoid_parcel_type%resize(new_size) + call this%ellipsoid_resize(new_size) call resize_array(this%theta, new_size, this%local_num) call resize_array(this%qv, new_size, this%local_num) @@ -204,4 +210,32 @@ subroutine prec_parcel_resize(this, new_size) end subroutine prec_parcel_resize + subroutine realistic_parcel_split(this, n, n_thread_loc) + class(realistic_parcel_type), intent(inout) :: this + integer, intent(in) :: n + integer, intent(in) :: n_thread_loc + call this%ellipsoid_split(n, n_thread_loc) + end subroutine realistic_parcel_split + + subroutine idealised_parcel_split(this, n, n_thread_loc) + class(idealised_parcel_type), intent(inout) :: this + integer, intent(in) :: n + integer, intent(in) :: n_thread_loc + call this%ellipsoid_split(n, n_thread_loc) + end subroutine idealised_parcel_split + + pure subroutine realistic_parcel_get_buoyancy(this, num, buoyancy) + class(realistic_parcel_type), intent(inout) :: this + integer, intent(in) :: num + double precision, intent(out) :: buoyancy + buoyancy = 0.0 + end subroutine realistic_parcel_get_buoyancy + + pure subroutine idealised_parcel_get_buoyancy(this, num, buoyancy) + class(idealised_parcel_type), intent(inout) :: this + integer, intent(in) :: num + double precision, intent(out) :: buoyancy + buoyancy = 0.0 + end subroutine idealised_parcel_get_buoyancy + end module From 6b316c17715009c34143dc7d0de3d5849a9a38b3 Mon Sep 17 00:00:00 2001 From: Steven Boeing Date: Wed, 25 Jun 2025 08:45:30 +0100 Subject: [PATCH 5/7] Implement generalised buoyancy --- src/2d/parcels/parcel_interpl.f90 | 11 +---------- src/Makefile.am | 2 +- src/parcel_types/parcel_container.f90 | 13 +++++++++++++ src/parcel_types/parcel_types.f90 | 23 ++++++++++++++++++++--- src/utils/physics.f90 | 24 ++++++++++++++---------- 5 files changed, 49 insertions(+), 24 deletions(-) diff --git a/src/2d/parcels/parcel_interpl.f90 b/src/2d/parcels/parcel_interpl.f90 index 04472cd9e..22201017c 100644 --- a/src/2d/parcels/parcel_interpl.f90 +++ b/src/2d/parcels/parcel_interpl.f90 @@ -168,17 +168,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)) diff --git a/src/Makefile.am b/src/Makefile.am index d60cb7065..7c80c1a73 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,4 +1,4 @@ -SUBDIRS = mpi netcdf parcel_types utils fft 2d +SUBDIRS = mpi netcdf utils parcel_types fft 2d if ENABLE_3D SUBDIRS += 3d diff --git a/src/parcel_types/parcel_container.f90 b/src/parcel_types/parcel_container.f90 index 8d76ec1e0..16a136ca4 100644 --- a/src/parcel_types/parcel_container.f90 +++ b/src/parcel_types/parcel_container.f90 @@ -42,6 +42,9 @@ module parcel_container integer :: n_pos = -1 ! number of spatial dimensions character(len=1), allocatable, dimension(:) :: pos_names ! Names of directions for positions character(len=8) :: dim_string ! Names of directions for positions + integer :: x_dim = -1 ! For reverse lookup + integer :: y_dim = -1 ! For reverse lookup + integer :: z_dim = -1 ! For reverse lookup contains procedure :: base_alloc => base_parcel_alloc @@ -201,6 +204,16 @@ subroutine base_parcel_alloc(this, num) call this%register_attribute(this%delta_pos(n, :), this%pos_names(n) // "_position_rk_tendency", "m/s") enddo + do n = 1, this%n_pos + if(this%pos_names(n)=='x') then + this%x_dim=n + elseif(this%pos_names(n)=='y') then + this%y_dim=n + elseif(this%pos_names(n)=='z') then + this%z_dim=n + endif + enddo + end subroutine base_parcel_alloc !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: diff --git a/src/parcel_types/parcel_types.f90 b/src/parcel_types/parcel_types.f90 index da04f3a90..445c2eddb 100644 --- a/src/parcel_types/parcel_types.f90 +++ b/src/parcel_types/parcel_types.f90 @@ -1,5 +1,6 @@ module parcel_types - + use physics, only : glat, lambda_c, q_0, qv_dens_coeff, theta_0, gravity + use constants, only : zero, one use parcel_ellipsoid implicit none @@ -228,14 +229,30 @@ pure subroutine realistic_parcel_get_buoyancy(this, num, buoyancy) class(realistic_parcel_type), intent(inout) :: this integer, intent(in) :: num double precision, intent(out) :: buoyancy - buoyancy = 0.0 + + if(this%is_moist) then + ! total buoyancy (including effects of latent heating) + buoyancy = gravity*(this%theta(num)*(one+qv_dens_coeff*this%qv(num)-this%ql(num))/theta_0-one) + else + buoyancy = gravity*(this%theta(num)/theta_0-one) + endif end subroutine realistic_parcel_get_buoyancy pure subroutine idealised_parcel_get_buoyancy(this, num, buoyancy) class(idealised_parcel_type), intent(inout) :: this integer, intent(in) :: num double precision, intent(out) :: buoyancy - buoyancy = 0.0 + double precision :: q_c + + !! NOTE: need to check height offset (see previous code) + if(this%is_moist) then + q_c = this%humidity(num) & + - q_0 * exp(- lambda_c * this%position(this%z_dim, num)) + q_c = max(zero, q_c) + buoyancy = this%buoyancy(num) + glat * q_c + else + buoyancy = this%buoyancy(num) + endif end subroutine idealised_parcel_get_buoyancy end module diff --git a/src/utils/physics.f90 b/src/utils/physics.f90 index 00da4564b..33e320543 100644 --- a/src/utils/physics.f90 +++ b/src/utils/physics.f90 @@ -70,7 +70,12 @@ module physics ![m] MPIC specific, scale-height, H double precision, protected :: height_c = 1000.0d0 - ! + ![-] Molecular weight of dry air/ molecular weight of water - 1 + double precision, protected :: qv_dens_coeff = 0.608 + + ! gas constant of gry air + double precision, protected :: r_d = 287.05 + ! The following quantities are calculated: ! @@ -131,14 +136,14 @@ subroutine read_physical_quantities_from_namelist(fname) logical :: exists = .false. ! namelist definitions - namelist /PHYSICS/ gravity, & - L_v, & - c_p, & - theta_0, & - ang_vel, & + namelist /PHYSICS/ gravity, & + L_v, & + c_p, & + theta_0, & + ang_vel, & l_planetary_vorticity, & - lat_degrees, & - q_0, & + lat_degrees, & + q_0, & height_c ! check whether file exists @@ -171,7 +176,6 @@ subroutine update_physical_quantities glat = gravity * L_v / (c_p * theta_0) glati = one / glat - lambda_c = one / height_c if (l_planetary_vorticity) then @@ -288,7 +292,7 @@ subroutine print_physical_quantities call print_physical_quantity('temperature at sea level', theta_0, 'K') call print_physical_quantity('planetary angular velocity', ang_vel, 'rad/s') call print_physical_quantity('saturation specific humidity at ground level', q_0) - call print_physical_quantity('planetary vorticity', l_planetary_vorticity) + call print_physical_quantity('l_planetary vorticity', l_planetary_vorticity) call print_physical_quantity('vertical planetary vorticity', f_cor(3), '1/s') call print_physical_quantity('horizontal planetary vorticity', f_cor(2), '1/s') call print_physical_quantity('latitude degrees', lat_degrees, 'deg') From 170f4606c50f9edb04d80a3e8ceb838110395fec Mon Sep 17 00:00:00 2001 From: Steven Boeing Date: Wed, 25 Jun 2025 17:01:17 +0100 Subject: [PATCH 6/7] Generalised splitting --- src/2d/parcels/parcel_container.f90 | 27 +-------------------------- src/2d/parcels/parcel_split.f90 | 13 +------------ src/parcel_types/parcel_ellipsoid.f90 | 19 +++++++++++++++++-- src/parcel_types/parcel_types.f90 | 23 +++++++++++++++++++---- 4 files changed, 38 insertions(+), 44 deletions(-) diff --git a/src/2d/parcels/parcel_container.f90 b/src/2d/parcels/parcel_container.f90 index d2a17a3d5..7d8a16fac 100644 --- a/src/2d/parcels/parcel_container.f90 +++ b/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 diff --git a/src/2d/parcels/parcel_split.f90 b/src/2d/parcels/parcel_split.f90 index 3891664e5..ab44618a7 100644 --- a/src/2d/parcels/parcel_split.f90 +++ b/src/2d/parcels/parcel_split.f90 @@ -67,7 +67,6 @@ subroutine split_ellipses(parcels, threshold) parcels%B(2, n) = B12 - f34 * a2 * (evec(1) * evec(2)) h = f14 * sqrt(three * a2) - parcels%volume(n) = f12 * V !$omp critical n_thread_loc = n_parcels + 1 @@ -76,17 +75,7 @@ subroutine split_ellipses(parcels, threshold) n_parcels = n_parcels + 1 !$omp end critical - - parcels%B(:, n_thread_loc) = parcels%B(:, n) - - parcels%vorticity(:, n_thread_loc) = parcels%vorticity(:, n) - parcels%volume(n_thread_loc) = parcels%volume(n) - parcels%buoyancy(n_thread_loc) = parcels%buoyancy(n) -#ifndef ENABLE_DRY_MODE - parcels%humidity(n_thread_loc) = parcels%humidity(n) -#endif - parcels%position(:, n_thread_loc) = parcels%position(:, n) - h * evec - parcels%position(:, n) = parcels%position(:, n) + h * evec + call parcels%split(n, n_thread_loc, h*evec) ! child parcels need to be reflected into domain, if their center ! is inside the halo region diff --git a/src/parcel_types/parcel_ellipsoid.f90 b/src/parcel_types/parcel_ellipsoid.f90 index cfb2e13cb..59e6da021 100644 --- a/src/parcel_types/parcel_ellipsoid.f90 +++ b/src/parcel_types/parcel_ellipsoid.f90 @@ -1,6 +1,7 @@ module parcel_ellipsoid use parcel_container + use constants, only : f12 implicit none ! Adding the ellipsoid geomerty to the dynamics @@ -29,11 +30,12 @@ module parcel_ellipsoid end type interface - subroutine parcel_split(this, n, n_thread_loc) + subroutine parcel_split(this, n, n_thread_loc, d_pos_split) import ellipsoid_parcel_type class(ellipsoid_parcel_type), intent(inout) :: this integer, intent(in) :: n integer, intent(in) :: n_thread_loc + double precision, intent(in) :: d_pos_split(:) end subroutine parcel_split end interface @@ -205,10 +207,23 @@ subroutine ellipsoid_parcel_resize(this, new_size) end subroutine ellipsoid_parcel_resize - subroutine ellipsoid_parcel_split(this, n, n_thread_loc) + subroutine ellipsoid_parcel_split(this, n, n_thread_loc, d_pos_split) class(ellipsoid_parcel_type), intent(inout) :: this integer, intent(in) :: n integer, intent(in) :: n_thread_loc + double precision, intent(in) :: d_pos_split(:) + + this%volume(n) = f12 * this%volume(n) + this%B(:, n_thread_loc) = this%B(:, n) + this%vorticity(:, n_thread_loc) = this%vorticity(:, n) + this%volume(n_thread_loc) = this%volume(n) + this%position(:, n_thread_loc) = this%position(:, n) - d_pos_split + this%position(:, n) = this%position(:, n) + d_pos_split + if(this%has_labels) then + this%label(n_thread_loc) = this%label(n) + this%dilution(n_thread_loc) = this%dilution(n) + endif + end subroutine ellipsoid_parcel_split end module diff --git a/src/parcel_types/parcel_types.f90 b/src/parcel_types/parcel_types.f90 index 445c2eddb..039bdba39 100644 --- a/src/parcel_types/parcel_types.f90 +++ b/src/parcel_types/parcel_types.f90 @@ -211,18 +211,33 @@ subroutine prec_parcel_resize(this, new_size) end subroutine prec_parcel_resize - subroutine realistic_parcel_split(this, n, n_thread_loc) + subroutine realistic_parcel_split(this, n, n_thread_loc, d_pos_split) class(realistic_parcel_type), intent(inout) :: this integer, intent(in) :: n integer, intent(in) :: n_thread_loc - call this%ellipsoid_split(n, n_thread_loc) + double precision, intent(in) :: d_pos_split(:) + + call this%ellipsoid_split(n, n_thread_loc, d_pos_split) + + this%theta(n_thread_loc) = this%theta(n) + if(this%is_moist) then + this%qv(n_thread_loc) = this%qv(n) + this%ql(n_thread_loc) = this%ql(n) + endif end subroutine realistic_parcel_split - subroutine idealised_parcel_split(this, n, n_thread_loc) + subroutine idealised_parcel_split(this, n, n_thread_loc, d_pos_split) class(idealised_parcel_type), intent(inout) :: this integer, intent(in) :: n integer, intent(in) :: n_thread_loc - call this%ellipsoid_split(n, n_thread_loc) + double precision, intent(in) :: d_pos_split(:) + + call this%ellipsoid_split(n, n_thread_loc, d_pos_split) + + this%buoyancy(n_thread_loc) = this%buoyancy(n) + if(this%is_moist) then + this%humidity(n_thread_loc) = this%humidity(n) + endif end subroutine idealised_parcel_split pure subroutine realistic_parcel_get_buoyancy(this, num, buoyancy) From c8e457cc45148cc988c0ed2af98d8ebe59c8a9b0 Mon Sep 17 00:00:00 2001 From: Steven Boeing Date: Tue, 1 Jul 2025 12:00:16 +0100 Subject: [PATCH 7/7] Work towards generalised parcel types --- src/2d/parcels/parcel_init.f90 | 6 +- src/2d/parcels/parcel_interpl.f90 | 7 +- src/2d/parcels/parcel_merge.f90 | 119 ++++++-------- src/2d/utils/utils.f90 | 2 +- src/parcel_types/parcel_container.f90 | 8 +- src/parcel_types/parcel_ellipsoid.f90 | 144 ++++++++++++++++- src/parcel_types/parcel_types.f90 | 216 +++++++++++++++++++++++++- 7 files changed, 409 insertions(+), 93 deletions(-) diff --git a/src/2d/parcels/parcel_init.f90 b/src/2d/parcels/parcel_init.f90 index 52cd591b1..6fa97e92e 100644 --- a/src/2d/parcels/parcel_init.f90 +++ b/src/2d/parcels/parcel_init.f90 @@ -93,9 +93,9 @@ subroutine init_parcels(fname, tol) do n = 1, n_parcels 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 diff --git a/src/2d/parcels/parcel_interpl.f90 b/src/2d/parcels/parcel_interpl.f90 index 22201017c..d5c206d00 100644 --- a/src/2d/parcels/parcel_interpl.f90 +++ b/src/2d/parcels/parcel_interpl.f90 @@ -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) @@ -199,9 +196,9 @@ subroutine par2grid + 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)) & diff --git a/src/2d/parcels/parcel_merge.f90 b/src/2d/parcels/parcel_merge.f90 index df27ddc44..a6b9c581c 100644 --- a/src/2d/parcels/parcel_merge.f90 +++ b/src/2d/parcels/parcel_merge.f90 @@ -8,6 +8,7 @@ module parcel_merge use parcel_container, only : get_delx use dynamic_parcels, only : n_parcels, parcels use parcel_types, only : idealised_parcel_type + use parcel_ellipsoid, only : ellipsoid_parcel_type use parcel_ellipse, only : get_B22, get_ab use options, only : parcel, verbose use parcel_bc @@ -67,11 +68,7 @@ end subroutine merge_ellipses ! @param[in] isma are the indices of the small parcels ! @param[in] iclo are the indices of the close parcels ! @param[in] n_merge is the array size of isma and iclo - ! @param[out] B11m are the B11 matrix entries of the mergers - ! @param[out] B12m are the B12 matrix entries of the mergers - ! @param[out] B22m are the B22 matrix entries of the mergers - ! @param[out] vm are the volumes of the mergers - subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm) + subroutine do_group_merge(parcels, isma, iclo, n_merge) type(idealised_parcel_type), intent(inout) :: parcels integer, intent(in) :: isma(0:) integer, intent(in) :: iclo(:) @@ -79,13 +76,7 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm) integer :: m, ic, is, l, n integer :: loca(n_parcels) double precision :: x0(n_merge) - double precision :: posm(2, n_merge), delx, vmerge, dely, B22, mu - double precision :: buoym(n_merge), vortm(n_merge) -#ifndef ENABLE_DRY_MODE - double precision :: hum(n_merge) -#endif - double precision, intent(out) :: B11m(n_merge), B12m(n_merge), B22m(n_merge), & - vm(n_merge) + double precision :: delx, vmerge, dely, B22, mu loca = zero @@ -99,75 +90,61 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm) loca(ic) = l ! vm will contain the total volume of the merged parcel - vm(l) = parcels%volume(ic) + parcels%merge_volume(l) = parcels%volume(ic) !x0 stores the x centre of the other parcel x0(l) = parcels%position(1, ic) ! posm(1, l) will sum v(is)*(x(is)-x(ic)) modulo periodicity - posm(1, l) = zero + parcels%merge_position(1, l) = zero ! posm(2, l) will contain v(ic)*z(ic)+sum{v(is)*z(is)} - posm(2, l) = parcels%volume(ic) * parcels%position(2, ic) + parcels%merge_position(2, l) = parcels%volume(ic) * parcels%position(2, ic) - ! buoyancy and humidity - buoym(l) = parcels%volume(ic) * parcels%buoyancy(ic) -#ifndef ENABLE_DRY_MODE - hum(l) = parcels%volume(ic) * parcels%humidity(ic) -#endif - vortm(l) = parcels%volume(ic) * parcels%vorticity(1, ic) + call parcels%assign_p2m(ic, l, parcels%volume(ic)) - B11m(l) = zero - B12m(l) = zero - B22m(l) = zero + parcels%merge_B(1, l) = zero + parcels%merge_B(2, l) = zero + parcels%merge_B(3, l) = zero !B22 endif ! Sum up all the small parcels merging with a common other one: ! "is" refers to the small parcel index is = isma(m) !Small parcel n = loca(ic) !Index of merged parcel - vm(n) = vm(n) + parcels%volume(is) !Accumulate volume of merged parcel + parcels%merge_volume(n) = parcels%merge_volume(n) + parcels%volume(is) !Accumulate volume of merged parcel ! works across periodic edge delx = get_delx(parcels%position(1, is), x0(n)) ! Accumulate sum of v(is)*(x(is)-x(ic)) - posm(1, n) = posm(1, n) + parcels%volume(is) * delx + parcels%merge_position(1, n) = parcels%merge_position(1, n) + parcels%volume(is) * delx ! Accumulate v(ic)*z(ic)+sum{v(is)*z(is)} - posm(2, n) = posm(2, n) + parcels%volume(is) * parcels%position(2, is) + parcels%merge_position(2, n) = parcels%merge_position(2, n) + parcels%volume(is) * parcels%position(2, is) - ! Accumulate buoyancy and humidity - buoym(n) = buoym(n) + parcels%volume(is) * parcels%buoyancy(is) -#ifndef ENABLE_DRY_MODE - hum(n) = hum(n) + parcels%volume(is) * parcels%humidity(is) -#endif - vortm(n) = vortm(n) + parcels%volume(is) * parcels%vorticity(1, is) + call parcels%add_p2m(is, n, parcels%volume(is)) enddo ! Obtain the merged parcel centres ! (l = total number of merged parcels) do m = 1, l ! temporary scalar containing 1 / vm(m) - vmerge = one / vm(m) + vmerge = one / parcels%merge_volume(m) ! x centre of merged parcel, modulo periodicity - posm(1, m) = vmerge * posm(1, m) + parcels%merge_position(1, m) = vmerge * parcels%merge_position(1, m) - posm(1, m) = x0(m) + posm(1, m) + parcels%merge_position(1, m) = x0(m) + parcels%merge_position(1, m) ! z centre of merged parcel - posm(2, m) = vmerge * posm(2, m) + parcels%merge_position(2, m) = vmerge * parcels%merge_position(2, m) ! need to correct position - call apply_periodic_bc(posm(:, m)) + call apply_periodic_bc(parcels%merge_position(:, m)) - ! buoyancy and humidity - buoym(m) = vmerge * buoym(m) -#ifndef ENABLE_DRY_MODE - hum(m) = vmerge * hum(m) -#endif - vortm(m) = vmerge * vortm(m) + ! normalise other properties by vmerge + call parcels%assign_m2m(m, vmerge) enddo loca = zero @@ -180,46 +157,41 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, B11m, B12m, B22m, vm) l = l + 1 loca(ic) = l - vmerge = one / vm(l) + vmerge = one / parcels%merge_volume(l) B22 = get_B22(parcels%B(1, ic), parcels%B(2, ic), parcels%volume(ic)) - delx = get_delx(parcels%position(1, ic), posm(1, l)) - dely = parcels%position(2, ic) - posm(2, l) + delx = get_delx(parcels%position(1, ic), parcels%merge_position(1, l)) + dely = parcels%position(2, ic) - parcels%merge_position(2, l) mu = parcels%volume(ic) * vmerge - B11m(l) = mu * (four * delx ** 2 + parcels%B(1, ic)) - B12m(l) = mu * (four * delx * dely + parcels%B(2, ic)) - B22m(l) = mu * (four * dely ** 2 + B22) + parcels%merge_B(1, l) = mu * (four * delx ** 2 + parcels%B(1, ic)) + parcels%merge_B(2, l) = mu * (four * delx * dely + parcels%B(2, ic)) + parcels%merge_B(3, l) = mu * (four * dely ** 2 + B22) - parcels%volume(ic) = vm(l) - parcels%position(1, ic) = posm(1, l) - parcels%position(2, ic) = posm(2, l) - - parcels%buoyancy(ic) = buoym(l) -#ifndef ENABLE_DRY_MODE - parcels%humidity(ic) = hum(l) -#endif - parcels%vorticity(1, ic) = vortm(l) + parcels%volume(ic) = parcels%merge_volume(l) + parcels%position(1, ic) = parcels%merge_position(1, l) + parcels%position(2, ic) = parcels%merge_position(2, l) + call parcels%assign_m2p(l, ic) endif is = isma(m) n = loca(ic) - vmerge = one / vm(n) + vmerge = one / parcels%merge_volume(n) - delx = get_delx(parcels%position(1, is), posm(1, n)) - dely = parcels%position(2, is) - posm(2, n) + delx = get_delx(parcels%position(1, is), parcels%merge_position(1, n)) + dely = parcels%position(2, is) - parcels%merge_position(2, n) B22 = get_B22(parcels%B(1, is), parcels%B(2, is), parcels%volume(is)) ! volume fraction A_{is} / A mu = vmerge * parcels%volume(is) - B11m(n) = B11m(n) + mu * (four * delx ** 2 + parcels%B(1, is)) - B12m(n) = B12m(n) + mu * (four * delx * dely + parcels%B(2, is)) - B22m(n) = B22m(n) + mu * (four * dely ** 2 + B22) + parcels%merge_B(1, n) = parcels%merge_B(1, n) + mu * (four * delx ** 2 + parcels%B(1, is)) + parcels%merge_B(2, n) = parcels%merge_B(2, n) + mu * (four * delx * dely + parcels%B(2, is)) + parcels%merge_B(3, n) = parcels%merge_B(3, n) + mu * (four * dely ** 2 + B22) enddo end subroutine do_group_merge @@ -238,12 +210,10 @@ subroutine geometric_merge(parcels, isma, iclo, n_merge) integer :: m, ic, l integer :: loca(n_parcels) double precision :: factor - double precision :: B11(n_merge), & - B12(n_merge), & - B22(n_merge), & - V(n_merge) - call do_group_merge(parcels, isma, iclo, n_merge, B11, B12, B22, V) + call parcels%merge_alloc(n_merge) + + call do_group_merge(parcels, isma, iclo, n_merge) loca = zero @@ -258,15 +228,18 @@ subroutine geometric_merge(parcels, isma, iclo, n_merge) ! normalize such that determinant of the merger is (ab)**2 ! ab / sqrt(det(B)) - factor = get_ab(V(l)) / sqrt(B11(l) * B22(l) - B12(l) ** 2) + factor = get_ab(parcels%merge_volume(l)) / & + sqrt(parcels%merge_B(1, l) * parcels%merge_B(3, l) - parcels%merge_B(2, l) ** 2) - parcels%B(1, ic) = B11(l) * factor - parcels%B(2, ic) = B12(l) * factor + parcels%B(1, ic) = parcels%merge_B(1, l) * factor + parcels%B(2, ic) = parcels%merge_B(2, l) * factor call apply_periodic_bc(parcels%position(:, ic)) endif enddo + call parcels%merge_dealloc + end subroutine geometric_merge diff --git a/src/2d/utils/utils.f90 b/src/2d/utils/utils.f90 index b3ee1847b..08cc316c4 100644 --- a/src/2d/utils/utils.f90 +++ b/src/2d/utils/utils.f90 @@ -216,7 +216,7 @@ subroutine setup_parcels parcels%dim_string='xz' call parcels%ellipsoid_dimensions() - parcels%is_moist=.true. + parcels%is_moist=.false. call parcels%alloc(max_num_parcels) if (l_restart) then diff --git a/src/parcel_types/parcel_container.f90 b/src/parcel_types/parcel_container.f90 index 16a136ca4..ef72ca5e0 100644 --- a/src/parcel_types/parcel_container.f90 +++ b/src/parcel_types/parcel_container.f90 @@ -74,7 +74,7 @@ module parcel_container double precision, allocatable, dimension(:,:) :: delta_vor double precision, allocatable, dimension(:) :: dilution integer(kind=8), allocatable, dimension(:) :: label - integer, private :: n_vor = -1 ! number of voriticity components + integer :: n_vor = -1 ! number of voriticity components character(len=1), allocatable, dimension(:) :: vor_names ! Names of vorticity components logical :: has_labels = .false. @@ -543,8 +543,10 @@ subroutine dynamic_parcel_dealloc(this) call try_deallocate(this%volume) call try_deallocate(this%vorticity) call try_deallocate(this%delta_vor) - call try_deallocate(this%dilution) - call try_deallocate(this%label) + if(this%has_labels) then + call try_deallocate(this%dilution) + call try_deallocate(this%label) + endif call try_deallocate(this%vor_names) call this%base_dealloc diff --git a/src/parcel_types/parcel_ellipsoid.f90 b/src/parcel_types/parcel_ellipsoid.f90 index 59e6da021..47c848fcb 100644 --- a/src/parcel_types/parcel_ellipsoid.f90 +++ b/src/parcel_types/parcel_ellipsoid.f90 @@ -12,6 +12,15 @@ module parcel_ellipsoid double precision, allocatable, dimension(:,:) :: Vetas double precision, allocatable, dimension(:,:) :: Vtaus double precision, allocatable, dimension(:,:) :: evec ! for now add evec for 2D. Need to think about consistency with 3D + ! Variables for merging. Note these are not attributes as the merges take place on the local domain + ! Also reallocation is done every time step with length n_merge + double precision, allocatable, dimension(:,:) :: merge_position + double precision, allocatable, dimension(:,:) :: merge_B + double precision, allocatable, dimension(:,:) :: merge_vorticity + double precision, allocatable, dimension(:) :: merge_volume + double precision, allocatable, dimension(:) :: merge_dilution + integer(kind=8), allocatable, dimension(:) :: merge_label + logical :: par2grid_1p = .false. ! use one point for par2grid logical :: grid2par_1p = .false. ! use one point for grid2par integer, private :: n_strain = -1 ! number of strain components @@ -24,12 +33,31 @@ module parcel_ellipsoid procedure :: ellipsoid_resize => ellipsoid_parcel_resize procedure :: ellipsoid_split => ellipsoid_parcel_split procedure :: ellipsoid_dimensions => set_ellipsoid_dimensions + procedure :: ellipsoid_merge_alloc => ellipsoid_parcel_merge_alloc + procedure :: ellipsoid_merge_dealloc => ellipsoid_parcel_merge_dealloc procedure(parcel_split), deferred :: split + procedure(parcel_merge_alloc), deferred :: merge_alloc + procedure(parcel_merge_dealloc), deferred :: merge_dealloc + procedure(parcel_assign_p2m), deferred :: assign_p2m + procedure(parcel_add_p2m), deferred :: add_p2m + procedure(parcel_assign_m2m), deferred :: assign_m2m + procedure(parcel_assign_m2p), deferred :: assign_m2p ! Other ellipsoid procedures to go here end type interface + subroutine parcel_merge_alloc(this, n_merge) + import ellipsoid_parcel_type + class(ellipsoid_parcel_type), intent(inout) :: this + integer, intent(in) :: n_merge + end subroutine parcel_merge_alloc + + subroutine parcel_merge_dealloc(this) + import ellipsoid_parcel_type + class(ellipsoid_parcel_type), intent(inout) :: this + end subroutine parcel_merge_dealloc + subroutine parcel_split(this, n, n_thread_loc, d_pos_split) import ellipsoid_parcel_type class(ellipsoid_parcel_type), intent(inout) :: this @@ -37,13 +65,43 @@ subroutine parcel_split(this, n, n_thread_loc, d_pos_split) integer, intent(in) :: n_thread_loc double precision, intent(in) :: d_pos_split(:) end subroutine parcel_split + + subroutine parcel_assign_p2m(this, n_in, n_out, temp_volume) + import ellipsoid_parcel_type + class(ellipsoid_parcel_type), intent(inout) :: this + integer, intent(in) :: n_in + integer, intent(in) :: n_out + double precision, intent(in) :: temp_volume + end subroutine + + subroutine parcel_add_p2m(this, n_in, n_out, temp_volume) + import ellipsoid_parcel_type + class(ellipsoid_parcel_type), intent(inout) :: this + integer, intent(in) :: n_in + integer, intent(in) :: n_out + double precision, intent(in) :: temp_volume + end subroutine + + subroutine parcel_assign_m2m(this, n_inout, temp_volume) + import ellipsoid_parcel_type + class(ellipsoid_parcel_type), intent(inout) :: this + integer, intent(in) :: n_inout + double precision, intent(in) :: temp_volume + end subroutine + + subroutine parcel_assign_m2p(this,n_in,n_out) + import ellipsoid_parcel_type + class(ellipsoid_parcel_type), intent(inout) :: this + integer, intent(in) :: n_in + integer, intent(in) :: n_out + end subroutine + end interface contains subroutine set_ellipsoid_dimensions(this) class(ellipsoid_parcel_type), intent(inout) :: this - integer :: i_dim call this%set_vorticity_dimensions @@ -226,4 +284,88 @@ subroutine ellipsoid_parcel_split(this, n, n_thread_loc, d_pos_split) end subroutine ellipsoid_parcel_split + subroutine ellipsoid_parcel_merge_alloc(this, n_merge) + class(ellipsoid_parcel_type), intent(inout) :: this + integer, intent(in) :: n_merge + + allocate(this%merge_position(this%n_pos, n_merge)) + allocate(this%merge_B(this%n_shape+1, n_merge)) ! Extra element for B22/B33 + allocate(this%merge_vorticity(this%n_vor, n_merge)) + allocate(this%merge_volume(n_merge)) + if(this%has_labels) then + allocate(this%merge_dilution(n_merge)) + allocate(this%merge_label(n_merge)) + endif + + end subroutine ellipsoid_parcel_merge_alloc + + subroutine ellipsoid_parcel_merge_dealloc(this) + class(ellipsoid_parcel_type), intent(inout) :: this + + call try_deallocate(this%merge_position) + call try_deallocate(this%merge_B) + call try_deallocate(this%merge_vorticity) + call try_deallocate(this%merge_volume) + if(this%has_labels) then + call try_deallocate(this%merge_dilution) + call try_deallocate(this%merge_label) + endif + + end subroutine ellipsoid_parcel_merge_dealloc + + subroutine ellipsoid_assign_p2m(this, n_in, n_out, temp_volume) + class(ellipsoid_parcel_type), intent(inout) :: this + integer, intent(in) :: n_in + integer, intent(in) :: n_out + double precision, intent(in) :: temp_volume + + this%merge_vorticity(:, n_out) = temp_volume * this%vorticity(:, n_in) + if(this%has_labels) then + this%merge_dilution(n_out) = this%dilution(n_in) + this%merge_label(n_out) = this%label(n_in) + endif + end subroutine ellipsoid_assign_p2m + + subroutine ellipsoid_add_p2m(this, n_in, n_out, temp_volume) + class(ellipsoid_parcel_type), intent(inout) :: this + integer, intent(in) :: n_in + integer, intent(in) :: n_out + double precision, intent(in) :: temp_volume + double precision :: rn + + this%merge_vorticity(:, n_out) = this%merge_vorticity(:, n_out) + temp_volume * this%vorticity(:, n_in) + if(this%has_labels) then + ! Dilute the parcel when volume is added for now + ! This could be optimised moving it to later in code + call random_number(rn) + if(rn*this%merge_volume(n_out) idealised_parcel_alloc - procedure :: dealloc=> idealised_parcel_dealloc + procedure :: dealloc => idealised_parcel_dealloc procedure :: resize => idealised_parcel_resize procedure :: split => idealised_parcel_split procedure :: get_buoyancy => idealised_parcel_get_buoyancy - + procedure :: merge_alloc => idealised_merge_alloc + procedure :: merge_dealloc => idealised_merge_dealloc + procedure :: assign_p2m => idealised_assign_p2m + procedure :: add_p2m => idealised_add_p2m + procedure :: assign_m2m => idealised_assign_m2m + procedure :: assign_m2p => idealised_assign_m2p end type type, extends(ellipsoid_parcel_type) :: realistic_parcel_type ! add procedures @@ -23,6 +30,10 @@ module parcel_types double precision, allocatable, dimension(:) :: ql double precision, allocatable, dimension(:) :: theta double precision, allocatable, dimension(:) :: Nl ! optional droplet number + double precision, allocatable, dimension(:) :: merge_qv + double precision, allocatable, dimension(:) :: merge_ql + double precision, allocatable, dimension(:) :: merge_theta + double precision, allocatable, dimension(:) :: merge_Nl logical :: is_moist = .false. logical :: has_droplets = .false. @@ -32,8 +43,12 @@ module parcel_types procedure :: resize => realistic_parcel_resize procedure :: split => realistic_parcel_split procedure :: get_buoyancy => realistic_parcel_get_buoyancy - - ! get_buoyancy added here + procedure :: merge_alloc => realistic_merge_alloc + procedure :: merge_dealloc => realistic_merge_dealloc + procedure :: assign_p2m => realistic_assign_p2m + procedure :: add_p2m => realistic_add_p2m + procedure :: assign_m2m => realistic_assign_m2m + procedure :: assign_m2p => realistic_assign_m2p end type type, extends(base_parcel_type) :: prec_parcel_type ! add procedures @@ -131,9 +146,15 @@ subroutine realistic_parcel_dealloc(this) class(realistic_parcel_type), intent(inout) :: this call try_deallocate(this%theta) - call try_deallocate(this%qv) - call try_deallocate(this%ql) - call try_deallocate(this%Nl) + + if(this%is_moist) then + call try_deallocate(this%qv) + call try_deallocate(this%ql) + endif + + if(this%has_droplets) then + call try_deallocate(this%Nl) + endif call this%ellipsoid_dealloc @@ -224,6 +245,9 @@ subroutine realistic_parcel_split(this, n, n_thread_loc, d_pos_split) this%qv(n_thread_loc) = this%qv(n) this%ql(n_thread_loc) = this%ql(n) endif + if(this%has_droplets) then + this%Nl(n_thread_loc) = this%Nl(n) + endif end subroutine realistic_parcel_split subroutine idealised_parcel_split(this, n, n_thread_loc, d_pos_split) @@ -270,4 +294,182 @@ pure subroutine idealised_parcel_get_buoyancy(this, num, buoyancy) endif end subroutine idealised_parcel_get_buoyancy + subroutine realistic_merge_alloc(this, n_merge) + class(realistic_parcel_type), intent(inout) :: this + integer, intent(in) :: n_merge + + call this%ellipsoid_merge_alloc(n_merge) + allocate(this%merge_theta(n_merge)) + if(this%is_moist) then + allocate(this%merge_qv(n_merge)) + allocate(this%merge_ql(n_merge)) + endif + if(this%has_droplets) then + allocate(this%merge_Nl(n_merge)) + endif + end subroutine realistic_merge_alloc + + subroutine idealised_merge_alloc(this, n_merge) + class(idealised_parcel_type), intent(inout) :: this + integer, intent(in) :: n_merge + + call this%ellipsoid_merge_alloc(n_merge) + allocate(this%merge_buoyancy(n_merge)) + if(this%is_moist) then + allocate(this%merge_humidity(n_merge)) + endif + + end subroutine idealised_merge_alloc + + subroutine realistic_merge_dealloc(this) + class(realistic_parcel_type), intent(inout) :: this + + call this%ellipsoid_merge_dealloc + call try_deallocate(this%merge_theta) + if(this%is_moist) then + call try_deallocate(this%merge_qv) + call try_deallocate(this%merge_ql) + endif + if(this%has_droplets) then + call try_deallocate(this%merge_Nl) + endif + + end subroutine realistic_merge_dealloc + + subroutine idealised_merge_dealloc(this) + class(idealised_parcel_type), intent(inout) :: this + + call this%ellipsoid_merge_dealloc + call try_deallocate(this%merge_buoyancy) + if(this%is_moist) then + call try_deallocate(this%merge_humidity) + endif + + end subroutine idealised_merge_dealloc + + subroutine idealised_assign_p2m(this, n_in, n_out, temp_volume) + class(idealised_parcel_type), intent(inout) :: this + integer, intent(in) :: n_in + integer, intent(in) :: n_out + double precision, intent(in) :: temp_volume + + call ellipsoid_assign_p2m(this, n_in, n_out, temp_volume) + + this%merge_buoyancy(n_out) = temp_volume * this%buoyancy(n_in) + if(this%is_moist) then + this%merge_humidity(n_out) = temp_volume * this%humidity(n_in) + endif + end subroutine idealised_assign_p2m + + subroutine idealised_add_p2m(this, n_in, n_out, temp_volume) + class(idealised_parcel_type), intent(inout) :: this + integer, intent(in) :: n_in + integer, intent(in) :: n_out + double precision, intent(in) :: temp_volume + + call ellipsoid_add_p2m(this, n_in, n_out, temp_volume) + + this%merge_buoyancy(n_out) = this%merge_buoyancy(n_out) + temp_volume * this%buoyancy(n_in) + if(this%is_moist) then + this%merge_humidity(n_out) = this%merge_humidity(n_out) + temp_volume * this%humidity(n_in) + endif + end subroutine idealised_add_p2m + + subroutine idealised_assign_m2m(this, n_inout, temp_volume) + class(idealised_parcel_type), intent(inout) :: this + integer, intent(in) :: n_inout + double precision, intent(in) :: temp_volume + + call ellipsoid_assign_m2m(this, n_inout, temp_volume) + + this%merge_buoyancy(n_inout) = temp_volume * this%merge_buoyancy(n_inout) + if(this%is_moist) then + this%merge_humidity(n_inout) = temp_volume * this%merge_humidity(n_inout) + endif + end subroutine idealised_assign_m2m + + subroutine idealised_assign_m2p(this, n_in, n_out) + class(idealised_parcel_type), intent(inout) :: this + integer, intent(in) :: n_in + integer, intent(in) :: n_out + + call ellipsoid_assign_m2p(this, n_in, n_out) + + this%buoyancy(n_out) = this%merge_buoyancy(n_in) + if(this%is_moist) then + this%humidity(n_out) = this%merge_humidity(n_in) + endif + end subroutine idealised_assign_m2p + + subroutine realistic_assign_p2m(this, n_in, n_out, temp_volume) + class(realistic_parcel_type), intent(inout) :: this + integer, intent(in) :: n_in + integer, intent(in) :: n_out + double precision, intent(in) :: temp_volume + + call ellipsoid_assign_p2m(this, n_in, n_out, temp_volume) + + this%merge_theta(n_out) = temp_volume * this%theta(n_in) + if(this%is_moist) then + this%merge_qv(n_out) = temp_volume * this%qv(n_in) + this%merge_ql(n_out) = temp_volume * this%ql(n_in) + endif + if(this%has_droplets) then + this%merge_Nl(n_out) = temp_volume * this%Nl(n_in) + endif + end subroutine realistic_assign_p2m + + subroutine realistic_add_p2m(this, n_in, n_out, temp_volume) + class(realistic_parcel_type), intent(inout) :: this + integer, intent(in) :: n_in + integer, intent(in) :: n_out + double precision, intent(in) :: temp_volume + + call ellipsoid_add_p2m(this, n_in, n_out, temp_volume) + + this%merge_theta(n_out) = this%merge_theta(n_out) + temp_volume * this%theta(n_in) + if(this%is_moist) then + this%merge_qv(n_out) = this%merge_qv(n_out) + temp_volume * this%qv(n_in) + this%merge_ql(n_out) = this%merge_ql(n_out) + temp_volume * this%ql(n_in) + endif + if(this%has_droplets) then + this%merge_Nl(n_out) = this%merge_Nl(n_out) + temp_volume * this%Nl(n_in) + endif + end subroutine realistic_add_p2m + + subroutine realistic_assign_m2m(this, n_inout, temp_volume) + class(realistic_parcel_type), intent(inout) :: this + integer, intent(in) :: n_inout + double precision, intent(in) :: temp_volume + + call ellipsoid_assign_m2m(this, n_inout, temp_volume) + + this%merge_theta(n_inout) = temp_volume * this%merge_theta(n_inout) + if(this%is_moist) then + this%merge_qv(n_inout) = temp_volume * this%merge_qv(n_inout) + this%merge_ql(n_inout) = temp_volume * this%merge_ql(n_inout) + endif + if(this%has_droplets) then + this%merge_Nl(n_inout) = temp_volume * this%merge_Nl(n_inout) + endif + end subroutine realistic_assign_m2m + + subroutine realistic_assign_m2p(this, n_in, n_out) + class(realistic_parcel_type), intent(inout) :: this + integer, intent(in) :: n_in + integer, intent(in) :: n_out + + call ellipsoid_assign_m2p(this, n_in, n_out) + + this%theta(n_out) = this%merge_theta(n_in) + if(this%is_moist) then + this%qv(n_out) = this%merge_qv(n_in) + this%ql(n_out) = this%merge_ql(n_in) + endif + if(this%has_droplets) then + this%Nl(n_out) = this%merge_Nl(n_in) + endif + end subroutine realistic_assign_m2p + + end module