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..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 parcel_container + use dynamic_parcels, only : parcels use parcel_bc use parcel_split, only : split_ellipses, split_timer use parcel_merge, only : merge_ellipses, merge_timer @@ -138,7 +138,7 @@ end subroutine run subroutine post_run use options, only : output - call parcel_dealloc + call parcels%dealloc call ls_rk4_dealloc call stop_timer(epic_timer) 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/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_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_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..6fa97e92e 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,11 +91,11 @@ subroutine init_parcels(fname, tol) !$omp parallel default(shared) !$omp do private(n) do n = 1, n_parcels - parcels%vorticity(n) = zero + parcels%vorticity(1, n) = zero parcels%buoyancy(n) = zero -#ifndef ENABLE_DRY_MODE - parcels%humidity(n) = zero -#endif + if(parcels%is_moist) then + parcels%humidity(n) = zero + endif enddo !$omp end do !$omp end parallel @@ -239,7 +238,7 @@ subroutine init_from_grids(ncfname, tol) if (has_dataset(ncid, 'vorticity')) then buffer = zero call read_netcdf_dataset(ncid, 'vorticity', buffer(0:nz, :), start=start, cnt=cnt) - call gen_parcel_scalar_attr(buffer, tol, parcels%vorticity) + call gen_parcel_scalar_attr(buffer, tol, parcels%vorticity(1, :)) endif diff --git a/src/2d/parcels/parcel_interpl.f90 b/src/2d/parcels/parcel_interpl.f90 index 0aa42264f..d5c206d00 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 @@ -142,9 +142,6 @@ subroutine par2grid double precision :: points(2, 2) integer :: n, p, l, i, j double precision :: pvol, weight, btot -#ifndef ENABLE_DRY_MODE - double precision :: q_c -#endif call start_timer(par2grid_timer) @@ -168,17 +165,8 @@ subroutine par2grid do n = 1, n_parcels pvol = parcels%volume(n) -#ifndef ENABLE_DRY_MODE - ! liquid water content - q_c = parcels%humidity(n) & - - q_0 * exp(lambda_c * (lower(2) - parcels%position(2, n))) - q_c = max(zero, q_c) + call parcels%get_buoyancy(n, btot) - ! total buoyancy (including effects of latent heating) - btot = parcels%buoyancy(n) + glat * q_c -#else - btot = parcels%buoyancy(n) -#endif points = get_ellipse_points(parcels%position(:, n), & pvol, parcels%B(:, n)) @@ -205,12 +193,12 @@ subroutine par2grid weight = f12 * weights(l) * pvol vortg(js(l), is(l)) = vortg(js(l), is(l)) & - + weight * parcels%vorticity(n) + + weight * parcels%vorticity(1, n) #ifndef ENABLE_DRY_MODE - dbuoyg(js(l), is(l)) = dbuoyg(js(l), is(l)) & + dbuoyg(js(l), is(l)) = dbuoyg(js(l), is(l)) & + weight * parcels%buoyancy(n) - humg(js(l), is(l)) = humg(js(l), is(l)) & + humg(js(l), is(l)) = humg(js(l), is(l)) & + weight * parcels%humidity(n) #endif tbuoyg(js(l), is(l)) = tbuoyg(js(l), is(l)) & diff --git a/src/2d/parcels/parcel_merge.f90 b/src/2d/parcels/parcel_merge.f90 index d4bae92dd..a6b9c581c 100644 --- a/src/2d/parcels/parcel_merge.f90 +++ b/src/2d/parcels/parcel_merge.f90 @@ -5,10 +5,10 @@ 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_ellipsoid, only : ellipsoid_parcel_type use parcel_ellipse, only : get_B22, get_ab use options, only : parcel, verbose use parcel_bc @@ -27,7 +27,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 @@ -68,25 +68,15 @@ 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) - type(parcel_container_type), intent(inout) :: parcels + 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(:) integer, intent(in) :: n_merge 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 @@ -100,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(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(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 @@ -181,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(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 @@ -232,19 +203,17 @@ 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 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 @@ -259,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 @@ -307,7 +279,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_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 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..ab44618a7 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 @@ -66,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 @@ -75,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/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..08cc316c4 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) + parcels%dim_string='xz' + call parcels%ellipsoid_dimensions() + parcels%is_moist=.false. + 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..7c80c1a73 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,4 +1,4 @@ -SUBDIRS = mpi netcdf utils fft 2d +SUBDIRS = mpi netcdf utils parcel_types fft 2d if ENABLE_3D SUBDIRS += 3d 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..ef72ca5e0 --- /dev/null +++ b/src/parcel_types/parcel_container.f90 @@ -0,0 +1,762 @@ +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 :: 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 + 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_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, 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 :: n_vor = -1 ! number of voriticity components + character(len=1), allocatable, dimension(:) :: vor_names ! Names of vorticity components + logical :: has_labels = .false. + + contains + 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 + + 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 + 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 + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + 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 + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + 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. + ! @param[in] num number of parcels + subroutine base_parcel_alloc(this, num) + class(base_parcel_type), intent(inout), target :: this + integer, intent(in) :: num + 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 + + 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, :), this%pos_names(n) // "_position", "m") + 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 + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + ! 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) + call try_deallocate(this%pos_names) + + 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 + 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 < 1) .or. (this%n_pos > 3)) then + print *, "ERROR: base_parcel_resize:: only 1-, 2- or 3-dimensional parcels allowed." + stop + endif + + 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, :), 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_dimensions(this) + class(base_parcel_type), intent(inout) :: this + integer :: n_dim, i_dim + + n_dim = len_trim(this%dim_string) + + if (this%n_pos /= -1) then + print *, "WARNING: Dimension already set." + return + endif + + if ((n_dim < 1) .or. (n_dim > 3)) then + print *, "ERROR: set_dimension:: only 1-, 2- or 3-dimensional parcels allowed." + stop + endif + + 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 + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + 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" + print *, name + 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" + print *, name + 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 + integer :: i_dim + + call this%base_alloc(num) + + allocate(this%volume(num)) + call this%register_attribute(this%volume, "volume", "m^3") + + if(.not. ((this%n_vor == 1) .or. (this%n_vor == 3))) then + print *, "Vorticity needs 1 or 3 components." + stop + 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)) + 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) + 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 + + 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) + + 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") + + 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) + 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..47c848fcb --- /dev/null +++ b/src/parcel_types/parcel_ellipsoid.f90 @@ -0,0 +1,371 @@ + module parcel_ellipsoid + + use parcel_container + use constants, only : f12 + implicit none + + ! Adding the ellipsoid geomerty to the dynamics + 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 + ! 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 + 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 :: 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 :: 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 + integer, intent(in) :: n + 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 + + 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)) + + 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%n_pos == 2) then + allocate(this%evec(2, num)) + call this%register_attribute(this%evec(1, :), "evec 1", "m") + call this%register_attribute(this%evec(2, :), "evec 2", "m") + elseif (this%n_pos == 3) then + allocate(this%Vetas(3, num)) + allocate(this%Vtaus(3, num)) + 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 *, "ERROR: ellipsoid_parcel_alloc:: n_pos must be 2 or 3." + 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 try_deallocate(this%shape_names) + call try_deallocate(this%strain_names) + + call this%dynamic_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 + integer :: i_dim + + 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) + + 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%evec(1, :), "evec 1") + call this%reset_attribute(this%evec(2, :), "evec 2") + 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%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 *, "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, 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 + + 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 :: 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 + double precision, allocatable, dimension(:) :: qv + 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. + + contains + 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 + 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 + 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_alloc(num) + + allocate(this%buoyancy(num)) + + 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 + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + 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_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_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_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) + + 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 + + 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_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 + + 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 + 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 + 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) + class(idealised_parcel_type), intent(inout) :: this + integer, intent(in) :: n + integer, intent(in) :: 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) + class(realistic_parcel_type), intent(inout) :: this + integer, intent(in) :: num + double precision, intent(out) :: buoyancy + + 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 + 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 + + 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 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')