diff --git a/esmf_utils/InfoUtilities.F90 b/esmf_utils/InfoUtilities.F90 index d6f758a71894..b9b91a4680b2 100644 --- a/esmf_utils/InfoUtilities.F90 +++ b/esmf_utils/InfoUtilities.F90 @@ -49,6 +49,7 @@ module mapl3g_InfoUtilities interface MAPL_InfoCreateFromInternal procedure :: info_field_create_from_internal + procedure :: info_bundle_create_from_internal end interface MAPL_InfoCreateFromInternal interface MAPL_InfoCreateFromShared @@ -113,7 +114,9 @@ module mapl3g_InfoUtilities interface MAPL_InfoGetInternal procedure :: info_field_get_internal_string procedure :: info_field_get_internal_i4 - procedure :: info_get_bundle_internal_r4_1d + procedure :: info_bundle_get_internal_string + procedure :: info_bundle_get_internal_i4 + procedure :: info_bundle_get_internal_r4_1d procedure :: info_stateitem_get_internal_string procedure :: info_stateitem_get_internal_logical procedure :: info_stateitem_get_internal_i4 @@ -123,8 +126,13 @@ module mapl3g_InfoUtilities end interface MAPL_InfoGetInternal interface MAPL_InfoSetInternal + procedure :: info_field_set_internal_info procedure :: info_field_set_internal_string procedure :: info_field_set_internal_i4 + procedure :: info_bundle_set_internal_info + procedure :: info_bundle_set_internal_string + procedure :: info_bundle_set_internal_i4 + procedure :: info_bundle_set_internal_r4_1d procedure :: info_stateitem_set_internal_string procedure :: info_stateitem_set_internal_logical procedure :: info_stateitem_set_internal_i4 @@ -254,28 +262,57 @@ end subroutine info_get_r4_1d ! MAPL_InfoCreateFromInternal - function info_field_create_from_internal(field, rc) result(info) + function info_field_create_from_internal(field, key, rc) result(info) type(ESMF_Info) :: info type(ESMF_Field), intent(in) :: field + character(*), optional, intent(in) :: key integer, optional, intent(out) :: rc type(ESMF_Info) :: host_info integer :: status + character(:), allocatable :: key_ call ESMF_InfoGetFromHost(field, host_info, _RC) - info = ESMF_InfoCreate(host_info, key=INFO_INTERNAL_NAMESPACE, _RC) + + key_ = INFO_INTERNAL_NAMESPACE + if (present(key)) then + key_ = concat(key_, key) + end if + + info = ESMF_InfoCreate(host_info, key=key_, _RC) _RETURN(_SUCCESS) end function info_field_create_from_internal - function info_field_create_from_shared(field, rc) result(info) + function info_bundle_create_from_internal(bundle, key, rc) result(info) type(ESMF_Info) :: info - type(ESMF_Field), intent(in) :: field + type(ESMF_FieldBundle), intent(in) :: bundle + character(*), optional, intent(in) :: key integer, optional, intent(out) :: rc type(ESMF_Info) :: host_info + character(:), allocatable :: key_ integer :: status + + key_ = INFO_INTERNAL_NAMESPACE + if (present(key)) then + key_ = concat(key_, key) + end if + call ESMF_InfoGetFromHost(bundle, host_info, _RC) + info = ESMF_InfoCreate(host_info, key=key_, _RC) + + _RETURN(_SUCCESS) + end function info_bundle_create_from_internal + + function info_field_create_from_shared(field, rc) result(info) + type(ESMF_Info) :: info + type(ESMF_Field), intent(in) :: field + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_Info) :: host_info + call ESMF_InfoGetFromHost(field, host_info, _RC) info = ESMF_InfoCreate(host_info, key=INFO_SHARED_NAMESPACE, _RC) @@ -828,7 +865,37 @@ subroutine info_field_get_internal_i4(field, key, value, rc) _RETURN(_SUCCESS) end subroutine info_field_get_internal_i4 - subroutine info_get_bundle_internal_r4_1d(bundle, key, values, rc) + subroutine info_bundle_get_internal_string(bundle, key, value, rc) + type(ESMF_FieldBundle), intent(in) :: bundle + character(*), intent(in) :: key + character(:), allocatable, intent(out) :: value + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_Info) :: info + + call ESMF_InfoGetFromHost(bundle, info, _RC) + call MAPL_InfoGet(info, key=concat(INFO_INTERNAL_NAMESPACE,key), value=value, _RC) + + _RETURN(_SUCCESS) + end subroutine info_bundle_get_internal_string + + subroutine info_bundle_get_internal_i4(bundle, key, value, rc) + type(ESMF_FieldBundle), intent(in) :: bundle + character(*), intent(in) :: key + integer(kind=ESMF_KIND_I4), intent(out) :: value + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_Info) :: info + + call ESMF_InfoGetFromHost(bundle, info, _RC) + call MAPL_InfoGet(info, key=concat(INFO_INTERNAL_NAMESPACE,key), value=value, _RC) + + _RETURN(_SUCCESS) + end subroutine info_bundle_get_internal_i4 + + subroutine info_bundle_get_internal_r4_1d(bundle, key, values, rc) type(ESMF_FieldBundle), intent(in) :: bundle character(*), intent(in) :: key real(kind=ESMF_KIND_R4), allocatable, intent(out) :: values(:) @@ -841,7 +908,7 @@ subroutine info_get_bundle_internal_r4_1d(bundle, key, values, rc) call MAPL_InfoGet(info, key=concat(INFO_INTERNAL_NAMESPACE,key), values=values, _RC) _RETURN(_SUCCESS) - end subroutine info_get_bundle_internal_r4_1d + end subroutine info_bundle_get_internal_r4_1d subroutine info_stateitem_get_internal_string(state, short_name, key, value, rc) type(ESMF_State), intent(in) :: state @@ -941,6 +1008,21 @@ end subroutine info_stateitem_get_internal_r4_1d ! MAPL_InfoSetInternal + subroutine info_field_set_internal_info(field, key, value, rc) + type(ESMF_Field), intent(in) :: field + character(*), intent(in) :: key + type(ESMF_Info), intent(in) :: value + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_Info) :: field_info + + call ESMF_InfoGetFromHost(field, field_info, _RC) + call MAPL_InfoSet(field_info, key=concat(INFO_INTERNAL_NAMESPACE,key), value=value, _RC) + + _RETURN(_SUCCESS) + end subroutine info_field_set_internal_info + subroutine info_field_set_internal_string(field, key, value, rc) type(ESMF_Field), intent(in) :: field character(*), intent(in) :: key @@ -971,6 +1053,66 @@ subroutine info_field_set_internal_i4(field, key, value, rc) _RETURN(_SUCCESS) end subroutine info_field_set_internal_i4 + subroutine info_bundle_set_internal_info(bundle, key, value, rc) + type(ESMF_FieldBundle), intent(inout) :: bundle + character(*), intent(in) :: key + type(ESMF_Info), intent(in) :: value + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_Info) :: bundle_info + + call ESMF_InfoGetFromHost(bundle, bundle_info, _RC) + call MAPL_InfoSet(bundle_info, key=concat(INFO_INTERNAL_NAMESPACE,key), value=value, _RC) + + _RETURN(_SUCCESS) + end subroutine info_bundle_set_internal_info + + subroutine info_bundle_set_internal_string(bundle, key, value, rc) + type(ESMF_FieldBundle), intent(inout) :: bundle + character(*), intent(in) :: key + character(*), intent(in) :: value + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_Info) :: info + + call ESMF_InfoGetFromHost(bundle, info, _RC) + call MAPL_InfoSet(info, key=concat(INFO_INTERNAL_NAMESPACE,key), value=value, _RC) + + _RETURN(_SUCCESS) + end subroutine info_bundle_set_internal_string + + subroutine info_bundle_set_internal_i4(bundle, key, value, rc) + type(ESMF_FieldBundle), intent(inout) :: bundle + character(*), intent(in) :: key + integer, intent(in) :: value + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_Info) :: info + + call ESMF_InfoGetFromHost(bundle, info, _RC) + call MAPL_InfoSet(info, key=concat(INFO_INTERNAL_NAMESPACE,key), value=value, _RC) + + _RETURN(_SUCCESS) + end subroutine info_bundle_set_internal_i4 + + subroutine info_bundle_set_internal_r4_1d(bundle, key, values, rc) + type(ESMF_FieldBundle), intent(inout) :: bundle + character(*), intent(in) :: key + real(kind=ESMF_KIND_R4), dimension(:), intent(in) :: values + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_Info) :: info + + call ESMF_InfoGetFromHost(bundle, info, _RC) + call MAPL_InfoSet(info, key=concat(INFO_INTERNAL_NAMESPACE,key), values=values, _RC) + + _RETURN(_SUCCESS) + end subroutine info_bundle_set_internal_r4_1d + subroutine info_stateitem_set_internal_string(state, short_name, key, value, rc) type(ESMF_State), intent(in) :: state character(*), intent(in) :: short_name diff --git a/esmf_utils/UngriddedDims.F90 b/esmf_utils/UngriddedDims.F90 index 100e4203e90e..1441d9675eb5 100644 --- a/esmf_utils/UngriddedDims.F90 +++ b/esmf_utils/UngriddedDims.F90 @@ -225,6 +225,9 @@ function make_ungriddedDims(info, key, rc) result(ungridded_dims) do i = 1, num_ungridded_dims dim_key = make_dim_key(i, _RC) + if (present(key)) then + dim_key = key // dim_key + end if dim_info = ESMF_InfoCreate(info, key=dim_key, _RC) dim_specs(i) = make_ungriddedDim(dim_info, _RC) call ESMF_InfoDestroy(dim_info, _RC) diff --git a/esmf_utils/tests/Test_FieldDimensionInfo.pf b/esmf_utils/tests/Test_FieldDimensionInfo.pf index cdbee53eb7c7..c3388f6af2f7 100644 --- a/esmf_utils/tests/Test_FieldDimensionInfo.pf +++ b/esmf_utils/tests/Test_FieldDimensionInfo.pf @@ -212,7 +212,7 @@ contains call ESMF_InfoSet(info, KEY_UNGRIDDED_DIMS // KEY_NUM_UNGRIDDED_DIMS, num_ungridded, _RC) do i=1, num_ungridded - key = make_dim_key(i, _RC) + key = KEY_UNGRIDDED_DIMS // make_dim_key(i, _RC) name = names_(i) units = units_(i) coord = coordinates_(i, :) diff --git a/field_utils/CMakeLists.txt b/field_utils/CMakeLists.txt index fec2a17ccc3e..645099bb52da 100644 --- a/field_utils/CMakeLists.txt +++ b/field_utils/CMakeLists.txt @@ -4,12 +4,15 @@ set(srcs FieldUtils.F90 FieldBLAS.F90 FieldPointerUtilities.F90 + FieldDelta.F90 FieldUtilities.F90 FieldUnaryFunctions.F90 FieldBinaryOperations.F90 FieldUnits.F90 FieldCondensedArray.F90 FieldCondensedArray_private.F90 + FieldDelta.F90 + FieldBundleDelta.F90 ) # To use extended udunits2 procedures, udunits2.c must be built and linked. diff --git a/field_utils/FieldBundleDelta.F90 b/field_utils/FieldBundleDelta.F90 new file mode 100644 index 000000000000..1b19c638edfb --- /dev/null +++ b/field_utils/FieldBundleDelta.F90 @@ -0,0 +1,288 @@ +! This class is to support propagation of time-dependent Field +! attributes across couplers as well as to provide guidance to the +! containt Action objects on when to recompute internal items. + +#include "MAPL_Exceptions.h" +module mapl3g_FieldBundleDelta + use mapl3g_LU_Bound + use mapl3g_FieldDelta + use mapl3g_InfoUtilities + use mapl_FieldUtilities + use mapl_FieldPointerUtilities + use mapl3g_esmf_info_keys + use mapl_ErrorHandling + use mapl_KeywordEnforcer + use esmf + implicit none (type, external) + private + + public :: FieldBundleDelta + + type :: FieldBundleDelta + private + type(FieldDelta) :: field_delta ! constant across bundle + real(ESMF_KIND_R4), allocatable :: interpolation_weights(:) + contains + procedure :: initialize_bundle_delta + generic :: initialize => initialize_bundle_delta + procedure :: update_bundle + procedure :: reallocate_bundle + end type FieldBundleDelta + + + interface FieldBundleDelta + procedure new_FieldBundleDelta + procedure new_FieldBundleDelta_field_delta + end interface FieldBundleDelta + +contains + + function new_FieldBundleDelta(fieldCount, geom, typekind, num_levels, units, interpolation_weights) result(bundle_delta) + type(FieldBundleDelta) :: bundle_delta + integer, optional, intent(in) :: fieldCount + type(ESMF_Geom), optional, intent(in) :: geom + type(ESMF_TypeKind_Flag), optional, intent(in) :: typekind + integer, optional, intent(in) :: num_levels + character(*), optional, intent(in) :: units + real(ESMF_KIND_R4), intent(in), optional :: interpolation_weights(:) + + associate (field_delta => FieldDelta(geom=geom, typekind=typekind, num_levels=num_levels, units=units)) + bundle_delta = FieldBundleDelta(field_delta, fieldCount, interpolation_weights) + end associate + + end function new_FieldBundleDelta + + function new_FieldBundleDelta_field_delta(field_delta, fieldCount, interpolation_weights) result(bundle_delta) + type(FieldBundleDelta) :: bundle_delta + type(FieldDelta), intent(in) :: field_delta + integer, optional, intent(in) :: fieldCount + real(ESMF_KIND_R4), optional, intent(in) :: interpolation_weights(:) + + bundle_delta%field_delta = field_delta + + if (present(interpolation_weights)) then + bundle_delta%interpolation_weights = interpolation_weights + end if + + end function new_FieldBundleDelta_field_delta + + + ! delta = bundle_b - bundle_a + subroutine initialize_bundle_delta(this, bundle_a, bundle_b, rc) + class(FieldBundleDelta), intent(out) :: this + type(ESMF_FieldBundle), intent(in) :: bundle_a + type(ESMF_FieldBundle), intent(in) :: bundle_b + integer, optional, intent(out) :: rc + + integer :: status + + call compute_interpolation_weights_delta(this%interpolation_weights, bundle_a, bundle_b, _RC) + call compute_field_delta(this%field_delta, bundle_a, bundle_b, _RC) + + _RETURN(_SUCCESS) + + + contains + + subroutine compute_interpolation_weights_delta(interpolation_weights, bundle_a, bundle_b, rc) + real(ESMF_KIND_R4), allocatable, intent(out) :: interpolation_weights(:) + type(ESMF_FieldBundle), intent(in) :: bundle_a + type(ESMF_FieldBundle), intent(in) :: bundle_b + integer, optional, intent(out) :: rc + + integer :: status + real(ESMF_KIND_R4), allocatable :: weights_a(:), weights_b(:) + + call MAPL_InfoGetInternal(bundle_a, key=KEY_INTERPOLATION_WEIGHTS, values=weights_a, _RC) + call MAPL_InfoGetInternal(bundle_b, key=KEY_INTERPOLATION_WEIGHTS, values=weights_b, _RC) + + if (any(weights_a /= weights_b)) then + interpolation_weights = weights_b + end if + + _RETURN(_SUCCESS) + + end subroutine compute_interpolation_weights_delta + + subroutine compute_field_delta(field_delta, bundle_a, bundle_b, rc) + type(FieldDelta), intent(out) :: field_delta + type(ESMF_FieldBundle), intent(in) :: bundle_a + type(ESMF_FieldBundle), intent(in) :: bundle_b + integer, optional, intent(out) :: rc + + integer :: status + integer :: fieldCount_a, fieldCount_b + type(ESMF_Field), allocatable :: fieldList_a(:), fieldList_b(:) + + call ESMF_FieldBundleGet(bundle_a, fieldCount=fieldCount_a, _RC) + call ESMF_FieldBundleGet(bundle_b, fieldCount=fieldCount_b, _RC) + allocate(fieldList_a(fieldCount_a), fieldList_b(fieldCount_b)) + + if ((fieldCount_a > 0) .and. (fieldCount_b > 0)) then + call ESMF_FieldBundleGet(bundle_a, fieldList=fieldList_a, _RC) + call ESMF_FieldBundleGet(bundle_b, fieldList=fieldList_b, _RC) + call field_delta%initialize(fieldList_a(1), fieldList_b(1), _RC) + _RETURN(_SUCCESS) + end if + + if (fieldCount_b > 0) then + call ESMF_FieldBundleGet(bundle_b, fieldList=fieldList_b, _RC) + ! full FieldDelta + call field_delta%initialize(fieldList_b(1), _RC) + _RETURN(_SUCCESS) + end if + + ! Otherwise nothing to do. Fields are either going away + ! (n_fields_b = 0) or there are no fields on either side + ! (n_fields_a = 0 and n_fields_b = 0). + + _RETURN(_SUCCESS) + end subroutine compute_field_delta + + + end subroutine initialize_bundle_delta + + subroutine update_bundle(this, bundle, ignore, rc) + class(FieldBundleDelta), intent(in) :: this + type(ESMF_FieldBundle), intent(inout) :: bundle + character(*), intent(in), optional :: ignore + integer, optional, intent(out) :: rc + + integer :: status + character(:), allocatable :: ignore_ + type(ESMF_Field), allocatable :: fieldList(:) + + ignore_ = '' + if (present(ignore)) ignore_ = ignore + + call this%reallocate_bundle(bundle, ignore=ignore_, _RC) + call MAPL_FieldBundleGet(bundle, fieldList=fieldList, _RC) + call this%field_delta%update_fields(fieldList, ignore=ignore_, _RC) + + ! unique attribute in bundle + call update_interpolation_weights(this%interpolation_weights, bundle, ignore=ignore_, _RC) + + _RETURN(_SUCCESS) + contains + + + subroutine update_interpolation_weights(interpolation_weights, bundle, ignore, rc) + real(ESMF_KIND_R4), optional, intent(in) :: interpolation_weights(:) + type(ESMF_FieldBundle), intent(inout) :: bundle + character(*), intent(in) :: ignore + integer, optional, intent(out) :: rc + + integer :: status + + _RETURN_UNLESS(present(interpolation_weights)) + _RETURN_IF(ignore == 'interpolation_weights') + + call MAPL_InfoSetInternal(bundle, KEY_INTERPOLATION_WEIGHTS, values=interpolation_weights, _RC) + + _RETURN(_SUCCESS) + end subroutine update_interpolation_weights + + end subroutine update_bundle + + + ! If the size of the bundle is not changing, then any reallocation is + ! relegated to fields through the FieldDelta component. + ! Otherwise we need to create or destroy fields in the bundle. + + subroutine reallocate_bundle(this, bundle, ignore, unusable, rc) + class(FieldBundleDelta), intent(in) :: this + type(ESMF_FieldBundle), intent(inout) :: bundle + character(*), intent(in) :: ignore + class(KeywordEnforcer), optional, intent(in) :: unusable + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_Field), allocatable :: fieldList(:) + type(ESMF_Geom) :: bundle_geom + integer :: i + type(LU_Bound), allocatable :: bounds(:) + type(LU_Bound) :: vertical_bounds + type(ESMF_TypeKind_Flag) :: typekind + integer, allocatable :: ungriddedLbound(:), ungriddedUbound(:) + type(ESMF_Info) :: ungridded_info + type(ESMF_Info) :: vertical_info + integer :: old_field_count, new_field_count + integer :: num_levels + character(:), allocatable :: units, vloc + character(ESMF_MAXSTR), allocatable :: fieldNameList(:) + + ! Easy case 1: field count unchanged + call MAPL_FieldBundleGet(bundle, fieldList=fieldList, _RC) + _RETURN_UNLESS(allocated(this%interpolation_weights)) + ! The number of weights is always one larger than the number of fields to support a constant + ! offset. ("Weights" is a funny term in that case.) + new_field_count = size(this%interpolation_weights) - 1 + old_field_count = size(fieldList) + _RETURN_IF(new_field_count == old_field_count) + + ! Easy case 2: field count changing to zero + if (new_field_count == 0) then! "/dev/null" case + call destroy_fields(fieldList, _RC) + _RETURN(_SUCCESS) + end if + + ! Hard case: need to create new fields? + _ASSERT(size(fieldList) == 0, 'fieldCount should only change to or from zero. ExtData use case.') + deallocate(fieldList) + allocate(fieldList(new_field_count)) + + ! Need geom, typekind, and bounds to allocate fields before + call MAPL_FieldBundleGet(bundle, geom=bundle_geom, _RC) + call MAPL_FieldBundleGet(bundle, typekind=typekind, ungriddedUBound=ungriddedUbound, _RC) + ungriddedLBound = [(1, i = 1, size(ungriddedUBound))] + + ungridded_info = MAPL_InfoCreateFromInternal(bundle, key=KEY_UNGRIDDED_DIMS, _RC) + call MAPL_InfoGetInternal(bundle, KEY_UNITS, value=units, _RC) + + call MAPL_InfoGetInternal(bundle, KEY_VLOC, value=vloc, _RC) + if (vloc /= "VERTICAL_DIM_NONE") then + call MAPL_InfoGetInternal(bundle, KEY_NUM_LEVELS, value=num_levels, _RC) + end if + + do i = 1, new_field_count + fieldList(i) = ESMF_FieldEmptyCreate(_RC) + call ESMF_FieldEmptySet(fieldList(i), geom=bundle_geom, _RC) + call ESMF_FieldEmptyComplete(fieldList(i), typekind=typekind, & + ungriddedLbound=ungriddedLBound, ungriddedUbound=ungriddedUBound, _RC) + call MAPL_InfoSetInternal(fieldList(i), KEY_UNGRIDDED_DIMS, value=ungridded_info, _RC) + call MAPL_InfoSetInternal(fieldList(i), KEY_VLOC, value=vloc, _RC) + if (vloc /= "VERTICAL_DIM_NONE") then + call MAPL_InfoSetInternal(fieldList(i), KEY_NUM_LEVELS, value=num_levels, _RC) + end if + call MAPL_InfoSetInternal(fieldList(i), KEY_UNITS, value=units, _RC) + end do + + call ESMF_InfoDestroy(ungridded_info, _RC) + + allocate(fieldNameList(old_field_count)) + call ESMF_FieldBundleGet(bundle, fieldNameList=fieldNameList, _RC) + call ESMF_FieldBundleRemove(bundle, fieldNameList, multiflag=.true., _RC) + + call ESMF_FieldBundleAdd(bundle, fieldList, multiFlag=.true., relaxedFlag=.true., _RC) + + _RETURN(_SUCCESS) + + contains + + subroutine destroy_fields(fieldList, rc) + type(ESMF_Field), intent(inout) :: fieldList(:) + integer, optional, intent(out) :: rc + + integer :: status + integer :: i + + do i = 1, size(fieldList) + call ESMF_FieldDestroy(fieldList(i), _RC) + end do + + _RETURN(_SUCCESS) + end subroutine destroy_fields + + end subroutine reallocate_bundle + +end module mapl3g_FieldBundleDelta diff --git a/field_utils/FieldDelta.F90 b/field_utils/FieldDelta.F90 new file mode 100644 index 000000000000..3cfabb903226 --- /dev/null +++ b/field_utils/FieldDelta.F90 @@ -0,0 +1,497 @@ +! This class is to support propagation of time-dependent Field +! attributes across couplers as well as to provide guidance to the +! containt Action objects on when to recompute internal items. + +#include "MAPL_Exceptions.h" +module mapl3g_FieldDelta + use mapl3g_InfoUtilities + use mapl_FieldPointerUtilities + use mapl3g_esmf_info_keys + use mapl_ErrorHandling + use mapl_KeywordEnforcer + use esmf + implicit none + private + + public :: FieldDelta + public :: operator(==), operator(/=) + + ! Allocatable components are used to indicate that the delta involves a + ! change in the relevant quantity. Unallocated means unchanged. + type :: FieldDelta + private + ! intrinsic + type(ESMF_Geom), allocatable :: geom + type(ESMF_TypeKind_Flag), allocatable :: typekind + ! info attributes + integer, allocatable :: num_levels + character(:), allocatable :: units + +!# logical :: geom_coords_changed = .false. +!# logical :: vgrid_coords_changed = .false. + contains + procedure :: initialize_field_delta + procedure :: initialize_field_delta_degenerate + generic :: initialize => initialize_field_delta + generic :: initialize => initialize_field_delta_degenerate + procedure :: update_field + procedure :: update_fields + procedure :: reallocate_field + procedure :: reallocate_fields + end type FieldDelta + + + interface FieldDelta + procedure new_FieldDelta + end interface FieldDelta + + + ! Will be in next release of ESMF (8.8?) + interface operator(==) + procedure :: ESMF_GeomEqual + end interface operator(==) + + interface operator(/=) + procedure :: ESMF_GeomNotEqual + end interface operator(/=) + +contains + + function new_FieldDelta(geom, typekind, num_levels, units) result(field_delta) + type(FieldDelta) :: field_delta + + type(ESMF_Geom), intent(in), optional :: geom + type(ESMF_TypeKind_Flag), intent(in), optional :: typekind + integer, intent(in), optional :: num_levels + character(*), intent(in), optional :: units + + if (present(geom)) then + field_delta%geom = geom + end if + + if (present(typekind)) then + field_delta%typekind = typekind + end if + + if (present(num_levels)) then + field_delta%num_levels = num_levels + end if + + if (present(units)) then + field_delta%units = units + end if + + end function new_FieldDelta + + + ! delta = f_b - f_a + subroutine initialize_field_delta(this, f_a, f_b, rc) + class(FieldDelta), intent(out) :: this + type(ESMF_Field), intent(in) :: f_a + type(ESMF_Field), intent(in) :: f_b + integer, optional, intent(out) :: rc + + integer :: status + + call compute_geom_delta(this%geom, f_a, f_b, _RC) + call compute_typekind_delta(this%typekind, f_a, f_b, _RC) + call compute_num_levels_delta(this%num_levels, f_a, f_b, _RC) + call compute_units_delta(this%units, f_a, f_b, _RC) + + _RETURN(_SUCCESS) + + + contains + + subroutine compute_geom_delta(geom, f_a, f_b, rc) + type(ESMF_Geom), allocatable, intent(out) :: geom + type(ESMF_Field), intent(in) :: f_a + type(ESMF_Field), intent(in) :: f_b + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_Geom):: geom_a, geom_b + + call ESMF_FieldGet(f_a, geom=geom_a, _RC) + call ESMF_FieldGet(f_b, geom=geom_b, _RC) + + if (geom_a /= geom_b) then + geom = geom_b + end if + + _RETURN(_SUCCESS) + + end subroutine compute_geom_delta + + subroutine compute_typekind_delta(typekind, f_a, f_b, rc) + type(ESMF_TypeKind_Flag), allocatable, intent(out) :: typekind + type(ESMF_Field), intent(in) :: f_a + type(ESMF_Field), intent(in) :: f_b + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_TypeKind_Flag) :: typekind_a, typekind_b + + call ESMF_FieldGet(f_a, typekind=typekind_a, _RC) + call ESMF_FieldGet(f_b, typekind=typekind_b, _RC) + + if (typekind_a /= typekind_b) then + typekind = typekind_b + end if + + _RETURN(_SUCCESS) + + end subroutine compute_typekind_delta + + subroutine compute_num_levels_delta(num_levels, f_a, f_b, rc) + integer, allocatable, intent(out) :: num_levels + type(ESMF_Field), intent(in) :: f_a + type(ESMF_Field), intent(in) :: f_b + integer, optional, intent(out) :: rc + + integer :: status + integer :: num_levels_a, num_levels_b + + call MAPL_InfoGetInternal(f_a, key=KEY_NUM_LEVELS, value=num_levels_a, _RC) + call MAPL_InfoGetInternal(f_b, key=KEY_NUM_LEVELS, value=num_levels_b, _RC) + + if (num_levels_a /= num_levels_b) then + num_levels = num_levels_b + end if + + _RETURN(_SUCCESS) + + end subroutine compute_num_levels_delta + + subroutine compute_units_delta(units, f_a, f_b, rc) + character(:), allocatable, intent(out) :: units + type(ESMF_Field), intent(in) :: f_a + type(ESMF_Field), intent(in) :: f_b + integer, optional, intent(out) :: rc + + integer :: status + character(len=:), allocatable :: units_a, units_b + + call MAPL_InfoGetInternal(f_a, KEY_UNITS, value=units_a, _RC) + call MAPL_InfoGetInternal(f_b, KEY_UNITS, value=units_b, _RC) + + if (units_a /= units_b) then + allocate(character(len_trim(units_b)) :: units) + units = units_b + end if + + _RETURN(_SUCCESS) + + end subroutine compute_units_delta + + end subroutine initialize_field_delta + + ! delta = f + subroutine initialize_field_delta_degenerate(this, f, rc) + class(FieldDelta), intent(out) :: this + type(ESMF_Field), intent(in) :: f + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_TypeKind_Flag) :: typekind + + allocate(this%geom) + allocate(this%typekind) + call ESMF_FieldGet(f, geom=this%geom, typekind=typekind, _RC) + + allocate(this%num_levels) + call MAPL_InfoGetInternal(f, KEY_NUM_LEVELS, value=this%num_levels, _RC) + call MAPL_InfoGetInternal(f, KEY_UNITS, value=this%units, _RC) + + _RETURN(_SUCCESS) + end subroutine initialize_field_delta_degenerate + + + + + subroutine update_field(this, field, ignore, rc) + class(FieldDelta), intent(in) :: this + type(ESMF_Field), intent(inout) :: field + character(*), intent(in), optional :: ignore + integer, optional, intent(out) :: rc + + integer :: status + character(:), allocatable :: ignore_ + + ignore_ = '' + if (present(ignore)) ignore_ = ignore + + call this%reallocate_field(field, ignore=ignore_, _RC) + + call update_num_levels(this%num_levels, field, ignore=ignore_, _RC) + call update_units(this%units, field, ignore=ignore, _RC) + + _RETURN(_SUCCESS) + contains + + subroutine update_num_levels(num_levels, field, ignore, rc) + integer, optional, intent(in) :: num_levels + type(ESMF_Field), intent(inout) :: field + character(*), intent(in) :: ignore + integer, optional, intent(out) :: rc + + integer :: status + + _RETURN_UNLESS(present(num_levels)) + _RETURN_IF(ignore == 'num_levels') + + call MAPL_InfoSetInternal(field, key=KEY_NUM_LEVELS, value=num_levels, _RC) + + _RETURN(_SUCCESS) + end subroutine update_num_levels + + subroutine update_units(units, field, ignore, rc) + character(*), optional, intent(in) :: units + type(ESMF_Field), intent(inout) :: field + character(*), intent(in) :: ignore + integer, optional, intent(out) :: rc + + integer :: status + + _RETURN_UNLESS(present(units)) + _RETURN_IF(ignore == 'units') + + call MAPL_InfoSetInternal(field, key=KEY_UNITS, value=units, _RC) + + _RETURN(_SUCCESS) + end subroutine update_units + + end subroutine update_field + + subroutine update_fields(this, fieldList, ignore, rc) + class(FieldDelta), intent(in) :: this + type(ESMF_Field), intent(inout) :: fieldList(:) + character(*), intent(in) :: ignore + integer, optional, intent(out) :: rc + + integer :: status + integer :: i + + do i = 1, size(fieldList) + call this%update_field(fieldList(i), ignore, _RC) + end do + + _RETURN(_SUCCESS) + end subroutine update_fields + + subroutine reallocate_field(this, field, ignore, unusable, rc) + class(FieldDelta), intent(in) :: this + type(ESMF_Field), intent(inout) :: field + character(*), optional, intent(in) :: ignore + class(KeywordEnforcer), optional, intent(in) :: unusable + integer, optional, intent(out) :: rc + + integer :: status + + type(ESMF_Geom) :: current_geom, geom + type(ESMF_TypeKind_Flag) :: current_typekind, typekind + + integer :: i, rank + integer, allocatable :: ungriddedLBound(:), ungriddedUBound(:) + integer, allocatable :: localElementCount(:), current_ungriddedUBound(:) + character(:), allocatable :: ignore_ + logical :: new_array + type(ESMF_FieldStatus_Flag) :: field_status + + new_array = .false. + ignore_ = '' + if (present(ignore)) ignore_ = ignore + + + call ESMF_FieldGet(field, status=field_status, _RC) + _ASSERT(field_status == ESMF_FIELDSTATUS_COMPLETE, 'field must at least have a geom.') + call ESMF_FieldGet(field, geom=current_geom, _RC) + + call ESMF_FieldGet(field, typekind=current_typekind, _RC) + localElementCount = FieldGetLocalElementCount(field, _RC) + + call select_geom(geom, current_geom, this%geom, ignore_, new_array) + call select_typekind(typekind, current_typekind, this%typekind, ignore_, new_array) + call select_ungriddedUbound(ungriddedUbound, field, this%num_levels, ignore_, new_array, _RC) + ungriddedLBound = [(1, i=1, size(ungriddedUBound))] + + _RETURN_UNLESS(new_array) + + call MAPL_EmptyField(field, _RC) + call ESMF_FieldEmptySet(field, geom, _RC) + + call ESMF_FieldEmptyComplete(field, & + typekind=typekind, & + ungriddedLBound=ungriddedLBound, ungriddedUbound=ungriddedUBound, & + _RC) + + _RETURN(_SUCCESS) + + contains + + subroutine select_geom(geom, current_geom, new_geom, ignore, new_array) + type(ESMF_Geom), intent(out) :: geom + type(ESMF_Geom), intent(in) :: current_geom + type(ESMF_Geom), optional, intent(in) :: new_geom + character(*), intent(in) :: ignore + logical, intent(inout) :: new_array + + geom = current_geom + + if (ignore == 'geom') return + if (.not. present(new_geom)) return + + new_array = new_array .or. (new_geom /= current_geom) + geom = new_geom + + end subroutine select_geom + + subroutine select_typekind(typekind, current_typekind, new_typekind, ignore, new_array) + type(ESMF_TypeKind_Flag), intent(out) :: typekind + type(ESMF_TypeKind_Flag), intent(in) :: current_typekind + type(ESMF_TypeKind_Flag), optional, intent(in) :: new_typekind + character(*), intent(in) :: ignore + logical, intent(inout) :: new_array + + typekind = current_typekind + + if (ignore == 'typekind') return + if (.not. present(new_typekind)) return + + new_array = new_array .or. (new_typekind /= current_typekind) + typekind = new_typekind + + end subroutine select_typekind + + subroutine select_ungriddedUbound(ungriddedUbound, field, new_num_levels, ignore, new_array, rc) + integer, allocatable, intent(out) :: ungriddedUbound(:) + type(ESMF_Field), intent(inout) :: field + integer, optional, intent(in) :: new_num_levels + character(*), intent(in) :: ignore + logical, intent(inout) :: new_array + integer, optional, intent(inout) :: rc + + integer :: status + character(:), allocatable :: vloc + integer :: ungriddedDimCount + integer :: rank + integer :: current_num_levels + integer, allocatable :: localElementCount(:) + integer, allocatable :: current_ungriddedUBound(:) + + call ESMF_FieldGet(field, & + ungriddedDimCount=ungriddedDimCount, & + rank=rank, _RC) + localElementCount = FieldGetLocalElementCount(field, _RC) + current_ungriddedUBound = localElementCount(rank-ungriddedDimCount+1:) + ungriddedUbound = current_ungriddedUBound + + if (ignore == 'num_levels') return + if (.not. present(new_num_levels)) return + + call MAPL_InfoGetInternal(field, KEY_NUM_LEVELS, value=current_num_levels, _RC) + call MAPL_InfoGetInternal(field, KEY_VLOC, value=vloc, _RC) + + ! Surface fields are not impacted by change in vertical grid + _RETURN_IF(vloc == 'VERTICAL_DIM_NONE') + + new_array = new_array .or. (this%num_levels /= current_num_levels) + + select case (vloc) + case ('VERTICAL_DIM_CENTER') + ungriddedUBound(1) = this%num_levels + case ('VERTICAL_DIM_EDGE') + ungriddedUBound(1) = this%num_levels + 1 + case default + _FAIL('unsupported vertical location: '//vloc) + end select + + _RETURN(_SUCCESS) + end subroutine select_ungriddedUbound + + end subroutine reallocate_field + + + subroutine reallocate_fields(this, fieldList, ignore, rc) + class(FieldDelta), intent(in) :: this + type(ESMF_Field), intent(inout) :: fieldList(:) + character(*), intent(in) :: ignore + integer, optional, intent(out) :: rc + + integer :: status + integer :: i + + do i = 1, size(fieldList) + call this%reallocate_field(fieldList(i), ignore, _RC) + end do + + _RETURN(_SUCCESS) + end subroutine reallocate_fields + + ! TODO - delete when next ESMF release provides support. + + impure elemental logical function ESMF_GeomEqual(geom1, geom2) + type(ESMF_Geom), intent(in) :: geom1, geom2 + + type(ESMF_GeomType_Flag) :: geomtype1, geomtype2 + type(ESMF_Grid) :: grid1, grid2 + type(ESMF_LocStream) :: locstream1, locstream2 + type(ESMF_Mesh) :: mesh1, mesh2 + type(ESMF_XGrid) :: xgrid1, xgrid2 + + ESMF_GeomEqual = .false. + + call ESMF_GeomGet(geom1, geomtype=geomtype1) + call ESMF_GeomGet(geom2, geomtype=geomtype2) + + if (geomtype1 /= geomtype2) return + + if (geomtype1 == ESMF_GEOMTYPE_GRID) then + call ESMF_GeomGet(geom1, grid=grid1) + call ESMF_GeomGet(geom2, grid=grid2) + ESMF_GeomEqual = (grid1 == grid2) + return + end if + + if (geomtype1 == ESMF_GEOMTYPE_LOCSTREAM) then + call ESMF_GeomGet(geom1, locstream=locstream1) + call ESMF_GeomGet(geom2, locstream=locstream2) + ESMF_GeomEqual = (locstream1 == locstream2) + return + end if + + if (geomtype1 == ESMF_GEOMTYPE_MESH) then + call ESMF_GeomGet(geom1, mesh=mesh1) + call ESMF_GeomGet(geom2, mesh=mesh2) + ESMF_GeomEqual = (mesh1 == mesh2) + return + end if + + if (geomtype1 == ESMF_GEOMTYPE_XGRID) then + call ESMF_GeomGet(geom1, xgrid=xgrid1) + call ESMF_GeomGet(geom2, xgrid=xgrid2) + ESMF_GeomEqual = (xgrid1 == xgrid2) + return + end if + + end function ESMF_GeomEqual + + + impure elemental logical function ESMF_GeomNotEqual(geom1, geom2) + type(ESMF_Geom), intent(in) :: geom1, geom2 + ESMF_GeomNotEqual = .not. (geom1 == geom2) + end function ESMF_GeomNotEqual + + subroutine MAPL_EmptyField(field, rc) + type(ESMF_Field), intent(inout) :: field + integer, optional, intent(out) :: rc + + integer :: status + + field%ftypep%status = ESMF_FIELDSTATUS_GRIDSET + call ESMF_ArrayDestroy(field%ftypep%array, _RC) + + _RETURN(_SUCCESS) + end subroutine MAPL_EmptyField + +end module mapl3g_FieldDelta diff --git a/field_utils/FieldUtilities.F90 b/field_utils/FieldUtilities.F90 index ddd630f95dad..d66a96209f3c 100644 --- a/field_utils/FieldUtilities.F90 +++ b/field_utils/FieldUtilities.F90 @@ -6,31 +6,21 @@ module MAPL_FieldUtilities use MAPL_FieldPointerUtilities use mapl3g_esmf_info_keys use mapl3g_InfoUtilities + use mapl3g_UngriddedDims + use mapl3g_LU_Bound use mapl_KeywordEnforcer use esmf implicit none (type, external) private - public :: FieldUpdate - public :: FieldReallocate public :: FieldIsConstant public :: FieldSet public :: FieldNegate public :: FieldPow - ! TODO delete these operators once ESMF supports == for geom - ! objects. - public :: operator(==) - public :: operator(/=) - interface FieldUpdate - procedure FieldUpdate_from_attributes - procedure FieldUpdate_from_field - end interface FieldUpdate - - interface FieldReallocate - procedure field_reallocate - end interface FieldReallocate + public :: MAPL_FieldBundleGet + public :: MAPL_FieldBundleSet interface FieldIsConstant procedure FieldIsConstantR4 @@ -41,97 +31,8 @@ module MAPL_FieldUtilities procedure FieldSet_R8 end interface FieldSet - ! Should be in ESMF someday ... - interface operator(==) - procedure :: ESMF_GeomEqual - end interface operator(==) - - interface operator(/=) - procedure :: ESMF_GeomNotEqual - end interface operator(/=) - contains - - subroutine field_reallocate(field, unusable, geom, typekind, num_levels, rc) - type(ESMF_Field), intent(inout) :: field - class(KeywordEnforcer), optional, intent(in) :: unusable - type(ESMF_Geom), optional, intent(in) :: geom - type(ESMF_TypeKind_Flag), optional, intent(in) :: typekind - integer, optional, intent(in) :: num_levels - integer, optional, intent(out) :: rc - - integer :: status - type(ESMF_Geom) :: old_geom, geom_ - type(ESMF_TypeKind_Flag) :: old_typekind, typekind_ - integer :: old_num_levels, num_levels_ - - logical :: skip_reallocate - integer :: ungriddedDimCount, rank - integer, allocatable :: localElementCount(:) - integer, allocatable :: old_ungriddedUBound(:) - integer, allocatable :: ungriddedUBound_(:), ungriddedLBound_(:) - integer :: i - - skip_reallocate = .true. - - call ESMF_FieldGet(field, typekind=old_typekind, geom=old_geom, ungriddedDimCount=ungriddedDimCount, rank=rank, _RC) - localElementCount = FieldGetLocalElementCount(field, _RC) - - typekind_ = old_typekind - if (present(typekind)) typekind_ = typekind - - geom_ = old_geom - if (present(geom)) geom_ = geom - - old_ungriddedUBound = localElementCount(rank-ungriddedDimCount+1:) - ungriddedUBound_ = old_ungriddedUBound - - old_num_levels = get_num_levels(field, _RC) - num_levels_ = old_num_levels - if (present(num_levels)) then - _ASSERT(num_levels_ > 0, 'Cannot add vertical dimension to field after initialization.') - _ASSERT(num_levels > 0, 'Cannot remove vertical dimension to field after initialization.') - num_levels_ = num_levels - - ungriddedUBound_ = old_ungriddedUBound - ungriddedUBound_(1) = num_levels_ ! Vertical dimension is always 1st ungridded dimension - end if - - if (typekind_ /= old_typekind) skip_reallocate = .false. - if (geom_ /= old_geom) skip_reallocate = .false. - if (num_levels_ /= old_num_levels) skip_reallocate = .false. - _RETURN_IF(skip_reallocate) - - call MAPL_EmptyField(field, _RC) - - call ESMF_FieldEmptySet(field, geom=geom_, _RC) - ungriddedLBound_ = [(1, i=1, size(ungriddedUBound_))] - call ESMF_FieldEmptyComplete(field, typekind=typekind_, ungriddedLBound=ungriddedLBound_, ungriddedUbound=ungriddedUBound_, _RC) - - ! Update info - if (num_levels_ /= old_num_levels) then - call MAPL_InfoSetInternal(field, key=KEY_NUM_LEVELS, value=num_levels_, _RC) - end if - - _RETURN(_SUCCESS) - end subroutine field_reallocate - - subroutine MAPL_EmptyField(field, rc) - type(ESMF_Field), intent(inout) :: field - integer, optional, intent(out) :: rc - - integer :: status - - field%ftypep%status = ESMF_FIELDSTATUS_GRIDSET - call ESMF_ArrayDestroy(field%ftypep%array, _RC) - - _RETURN(_SUCCESS) - end subroutine MAPL_EmptyField - - - - function FieldIsConstantR4(field,constant_val,rc) result(field_is_constant) logical :: field_is_constant type(ESMF_Field), intent(inout) :: field @@ -303,116 +204,127 @@ subroutine FieldPow(field_out,field_in,expo,rc) _RETURN(ESMF_SUCCESS) end subroutine FieldPow - impure elemental logical function ESMF_GeomEqual(geom1, geom2) - type(ESMF_Geom), intent(in) :: geom1, geom2 - - type(ESMF_GeomType_Flag) :: geomtype1, geomtype2 - type(ESMF_Grid) :: grid1, grid2 - type(ESMF_LocStream) :: locstream1, locstream2 - type(ESMF_Mesh) :: mesh1, mesh2 - type(ESMF_XGrid) :: xgrid1, xgrid2 - - ESMF_GeomEqual = .false. - - call ESMF_GeomGet(geom1, geomtype=geomtype1) - call ESMF_GeomGet(geom2, geomtype=geomtype2) - - if (geomtype1 /= geomtype2) return - - if (geomtype1 == ESMF_GEOMTYPE_GRID) then - call ESMF_GeomGet(geom1, grid=grid1) - call ESMF_GeomGet(geom2, grid=grid2) - ESMF_GeomEqual = (grid1 == grid2) - return - end if - if (geomtype1 == ESMF_GEOMTYPE_LOCSTREAM) then - call ESMF_GeomGet(geom1, locstream=locstream1) - call ESMF_GeomGet(geom2, locstream=locstream2) - ESMF_GeomEqual = (locstream1 == locstream2) - return - end if + ! Supplement ESMF + subroutine MAPL_FieldBundleGet(fieldBundle, unusable, fieldList, geom, typekind, ungriddedUbound, rc) + type(ESMF_FieldBundle), intent(in) :: fieldBundle + class(KeywordEnforcer), optional, intent(in) :: unusable + type(ESMF_Field), optional, allocatable, intent(out) :: fieldList(:) + type(ESMF_Geom), optional, intent(out) :: geom + type(ESMF_TypeKind_Flag), optional, intent(out) :: typekind + integer, allocatable, optional, intent(out) :: ungriddedUbound(:) + integer, optional, intent(out) :: rc - if (geomtype1 == ESMF_GEOMTYPE_MESH) then - call ESMF_GeomGet(geom1, mesh=mesh1) - call ESMF_GeomGet(geom2, mesh=mesh2) - ESMF_GeomEqual = (mesh1 == mesh2) - return + integer :: status + integer :: fieldCount + type(ESMF_GeomType_Flag) :: geomtype + character(:), allocatable :: typekind_str + type(ESMF_Info) :: ungridded_info + type(UngriddedDims) :: ungridded_dims + type(LU_Bound), allocatable :: bounds(:) + integer :: num_levels + character(:), allocatable :: vloc + + if (present(fieldList)) then + call ESMF_FieldBundleGet(fieldBundle, fieldCount=fieldCount, _RC) + allocate(fieldList(fieldCount)) + call ESMF_FieldBundleGet(fieldBundle, fieldList=fieldList, itemOrderflag=ESMF_ITEMORDER_ADDORDER, _RC) end if - if (geomtype1 == ESMF_GEOMTYPE_XGRID) then - call ESMF_GeomGet(geom1, xgrid=xgrid1) - call ESMF_GeomGet(geom2, xgrid=xgrid2) - ESMF_GeomEqual = (xgrid1 == xgrid2) - return + if (present(geom)) then + call get_geom(fieldBundle, geom, rc) end if - - end function ESMF_GeomEqual + if (present(typekind)) then + call MAPL_InfoGetInternal(fieldBundle, key=KEY_TYPEKIND, value=typekind_str, _RC) + select case (typekind_str) + case ('R4') + typekind = ESMF_TYPEKIND_R4 + case ('R8') + typekind = ESMF_TYPEKIND_R8 + case ('I4') + typekind = ESMF_TYPEKIND_I4 + case ('I8') + typekind = ESMF_TYPEKIND_I8 + case ('LOGICAL') + typekind = ESMF_TYPEKIND_LOGICAL + case default + _FAIL('unsupported typekind') + end select + end if - impure elemental logical function ESMF_GeomNotEqual(geom1, geom2) - type(ESMF_Geom), intent(in) :: geom1, geom2 - ESMF_GeomNotEqual = .not. (geom1 == geom2) - end function ESMF_GeomNotEqual - + if (present(ungriddedUbound)) then + ungridded_info = MAPL_InfoCreateFromInternal(fieldBundle, _RC) + ungridded_dims = make_ungriddedDims(ungridded_info, KEY_UNGRIDDED_DIMS, _RC) + bounds = ungridded_dims%get_bounds() + + call MAPL_InfoGetInternal(fieldBundle, key=KEY_VLOC, value=vloc, _RC) + if (vloc /= 'VERTICAL_DIM_NONE') then + call MAPL_InfoGetInternal(fieldBundle, key=KEY_NUM_LEVELS, value=num_levels, _RC) + select case (vloc) + case ('VERTICAL_DIM_CENTER') + bounds = [LU_Bound(1, num_levels), bounds] + case ('VERTICAL_DIM_EDGE') + bounds = [LU_Bound(1, num_levels+1), bounds] + case default + _FAIL('unsupported vertical location') + end select + end if - subroutine FieldUpdate_from_attributes(field, unusable, geom, num_levels, typekind, units, rc) - type(ESMF_Field), intent(inout) :: field - class(KeywordEnforcer), optional, intent(in) :: unusable - type(ESMF_Geom), optional, intent(in) :: geom - integer, optional, intent(in) :: num_levels - type(ESMF_TypeKind_Flag), optional, intent(in) :: typekind - character(len=*), optional, intent(in) :: units - integer, optional, intent(out) :: rc + ungriddedUbound = bounds%upper + end if - integer :: status + _RETURN(_SUCCESS) - call FieldReallocate(field, geom=geom, typekind=typekind, num_levels=num_levels, rc=rc) + contains - if (present(units)) then - call MAPL_InfoSetInternal(field, key=KEY_UNITS, value=units, _RC) - end if + subroutine get_geom(fieldBundle, geom, rc) + type(ESMF_FieldBundle), intent(in) :: fieldBundle + type(ESMF_Geom), intent(inout) :: geom + integer, optional, intent(out) :: rc - _RETURN(_SUCCESS) - - end subroutine FieldUpdate_from_attributes + integer :: status + type(ESMF_GeomType_Flag) :: geomtype + type(ESMF_Grid) :: grid + call ESMF_FieldBundleGet(fieldBundle, geomtype=geomtype, _RC) + if (geomtype == ESMF_GEOMTYPE_GRID) then + call ESMF_FieldBundleGet(fieldBundle, grid=grid, _RC) + ! memory leak + geom = ESMF_GeomCreate(grid=grid, _RC) + _RETURN(_SUCCESS) + end if - subroutine FieldUpdate_from_field(field, reference_field, ignore, rc) - type(ESMF_Field), intent(inout) :: field - type(ESMF_Field), intent(in) :: reference_field - character(*), optional, intent(in) :: ignore - integer, intent(out), optional :: rc + _FAIL('unsupported geomtype; needs simple extension') - integer :: status - integer, allocatable :: num_levels - type(ESMF_Geom), allocatable :: geom - type(ESMF_TypeKind_Flag), allocatable :: typekind - character(:), allocatable :: units - - if (ignore /= 'geom') then - allocate(geom) - call ESMF_FieldGet(reference_field, geom=geom,_RC) - end if + _RETURN(_SUCCESS) + end subroutine get_geom - if (ignore /= 'typekind') then - allocate(typekind) - call ESMF_FieldGet(reference_field, typekind=typekind, _RC) - end if + end subroutine MAPL_FieldBundleGet - if (ignore /= 'units') then - call MAPL_InfoGetInternal(reference_field, key=KEY_UNITS, value=units, _RC) - end if + subroutine MAPL_FieldBundleSet(fieldBundle, unusable, geom, rc) + type(ESMF_FieldBundle), intent(inout) :: fieldBundle + class(KeywordEnforcer), optional, intent(in) :: unusable + type(ESMF_Geom), optional, intent(in) :: geom + integer, optional, intent(out) :: rc - if (ignore /= 'num_levels') then - num_levels = get_num_levels(reference_field, _RC) + integer :: status + type(ESMF_GeomType_Flag) :: geomtype + type(ESMF_Grid) :: grid + + if (present(geom)) then + call ESMF_GeomGet(geom, geomtype=geomtype, _RC) + if (geomtype == ESMF_GEOMTYPE_GRID) then + call ESMF_GeomGet(geom, grid=grid, _RC) + call ESMF_FieldBundleSet(fieldBundle, grid=grid, _RC) + _RETURN(_SUCCESS) + end if + _FAIL('unsupported geomtype') end if - call FieldUpdate(field, geom=geom, typekind=typekind, num_levels=num_levels, units=units, _RC) - _RETURN(_SUCCESS) - - end subroutine FieldUpdate_from_field + end subroutine MAPL_FieldBundleSet + end module MAPL_FieldUtilities diff --git a/field_utils/tests/CMakeLists.txt b/field_utils/tests/CMakeLists.txt index 5c982070df8f..acf2e9837803 100644 --- a/field_utils/tests/CMakeLists.txt +++ b/field_utils/tests/CMakeLists.txt @@ -5,7 +5,8 @@ set (test_srcs Test_FieldBLAS.pf Test_FieldArithmetic.pf Test_FieldCondensedArray_private.pf - Test_FieldUtilities.pf + Test_FieldDelta.pf + Test_FieldBundleDelta.pf ) diff --git a/field_utils/tests/Test_FieldBundleDelta.pf b/field_utils/tests/Test_FieldBundleDelta.pf new file mode 100644 index 000000000000..ef9e974d4b0d --- /dev/null +++ b/field_utils/tests/Test_FieldBundleDelta.pf @@ -0,0 +1,530 @@ +#include "MAPL_TestErr.h" +#include "unused_dummy.H" +module Test_FieldBundleDelta + use mapl3g_FieldBundleDelta + use mapl3g_FieldDelta + use mapl3g_ESMF_Info_Keys + use mapl3g_InfoUtilities + use mapl_FieldUtilities + use mapl3g_UngriddedDims + use mapl3g_UngriddedDim + use mapl3g_LU_Bound + use esmf + use ESMF_TestMethod_mod + use funit + implicit none (type, external) + + real, parameter :: FILL_VALUE = 99. + real, parameter :: DEFAULT_WEIGHTS(*) = [0.0, 0.5, 0.5] + integer, parameter :: FIELD_COUNT = 2 + integer, parameter :: NUM_LEVELS = 3 + integer, parameter :: NUM_RADII = 5 + +contains + + subroutine setup_geom(geom, im) + type(ESMF_Geom), intent(out) :: geom + integer, intent(in) :: im + + type(ESMF_Grid) :: grid + + grid = ESMF_GridCreateNoPeriDim(maxIndex=[IM,IM]) + geom = ESMF_GeomCreate(grid) + + end subroutine setup_geom + + subroutine teardown_geom(geom) + type(ESMF_Geom), intent(inout) :: geom + + type(ESMF_Grid) :: grid + + call ESMF_GeomGet(geom, grid=grid) + call ESMF_GridDestroy(grid) + call ESMF_GeomDestroy(geom) + + end subroutine teardown_geom + + subroutine setup_field(field, geom, typekind, units, with_ungridded) + type(ESMF_Field), intent(out) :: field + type(ESMF_Geom), intent(in) :: geom + type(ESMF_TypeKind_Flag), intent(in) :: typekind + character(len=*), intent(in) :: units + logical, optional, intent(in) :: with_ungridded + + type(UngriddedDims) :: ungridded_dims + type(ESMF_Info) :: ungridded_info + type(LU_Bound), allocatable :: bounds(:) + + field = ESMF_FieldEmptyCreate() + call ESMF_FieldEmptySet(field, geom=geom) + + call MAPL_InfoSetInternal(field, KEY_UNITS, units) + + call MAPL_InfoSetInternal(field, KEY_VLOC, "VERTICAL_DIM_NONE") + + ungridded_dims = UngriddedDims() + bounds = ungridded_dims%get_bounds() + if (present(with_ungridded)) then + if (with_ungridded) then + call MAPL_InfoSetInternal(field, KEY_VLOC, "VERTICAL_DIM_CENTER") + call MAPL_InfoSetInternal(field, KEY_NUM_LEVELS, NUM_LEVELS) + call ungridded_dims%add_dim(UngriddedDim(NUM_RADII, "radius", 'nm')) + bounds = [LU_Bound(1, NUM_LEVELS), ungridded_dims%get_bounds()] + end if + end if + + ungridded_info = ungridded_dims%make_info() + call MAPL_InfoSetInternal(field, KEY_UNGRIDDED_DIMS, value=ungridded_info) + + call ESMF_FieldEmptyComplete(field, typekind=typekind, ungriddedLBound=bounds%lower, ungriddedUBound=bounds%upper) + call FieldSet(field, FILL_VALUE) + + end subroutine setup_field + + subroutine teardown_field(field) + type(ESMF_Field), intent(inout) :: field + + call ESMF_FieldDestroy(field) + + end subroutine teardown_field + + subroutine setup_bundle(bundle, weights, geom, typekind, units, with_ungridded) + type(ESMF_FieldBundle), intent(out) :: bundle + real(kind=ESMF_KIND_R4), intent(in) :: weights(:) + type(ESMF_Geom), intent(in) :: geom + type(ESMF_TypeKind_Flag), intent(in) :: typekind + character(len=*), intent(in) :: units + logical, optional, intent(in) :: with_ungridded + + integer :: i + type(ESMF_Field) :: f + integer :: fieldCount + type(UngriddedDims) :: ungridded_dims + type(ESMF_Info) :: ungridded_info + + bundle = ESMF_FieldBundleCreate() + call MAPL_FieldBundleSet(bundle, geom=geom) + fieldCount = size(weights) - 1 + do i = 1, fieldCount + call setup_field(f, geom, typekind, units, with_ungridded=with_ungridded) + call ESMF_FieldBundleAdd(bundle, [f], multiflag=.true.) + end do + + ungridded_dims = UngriddedDims() + ungridded_info = ungridded_dims%make_info() + call MAPL_InfoSetInternal(bundle, KEY_UNGRIDDED_DIMS, value=ungridded_info) + + call MAPL_InfoSetInternal(bundle, KEY_INTERPOLATION_WEIGHTS, weights) + if (typekind == ESMF_TYPEKIND_R4) then + call MAPL_InfoSetInternal(bundle, KEY_TYPEKIND, "R4") + else + call MAPL_InfoSetInternal(bundle, KEY_TYPEKIND, "R8") + end if + call MAPL_InfoSetInternal(bundle, KEY_UNITS, units) + + call MAPL_InfoSetInternal(bundle, KEY_VLOC, "VERTICAL_DIM_NONE") + ungridded_dims = UngriddedDims() + + if (present(with_ungridded)) then + if (with_ungridded) then + call MAPL_InfoSetInternal(bundle, KEY_VLOC, "VERTICAL_DIM_CENTER") + call MAPL_InfoSetInternal(bundle, KEY_NUM_LEVELS, NUM_LEVELS) + call ungridded_dims%add_dim(UngriddedDim(NUM_RADII, "radius", 'nm')) + end if + end if + + ungridded_info = ungridded_dims%make_info() + call MAPL_InfoSetInternal(bundle, KEY_UNGRIDDED_DIMS, value=ungridded_info) + + end subroutine setup_bundle + + subroutine teardown_bundle(bundle) + type(ESMF_FieldBundle), intent(inout) :: bundle + type(ESMF_Field), allocatable :: fieldList(:) + + integer :: i + + call MAPL_FieldBundleGet(bundle, fieldList=fieldList) + do i = 1, size(fieldList) + call ESMF_FieldDestroy(fieldList(i)) + end do + call ESMF_FieldBundleDestroy(bundle) + + end subroutine teardown_bundle + + @test(type=ESMF_TestMethod, npes=[1]) + subroutine test_change_typekind(this) + class(ESMF_TestMethod), intent(inout) :: this + + integer :: status + type(Fieldbundledelta) :: delta + type(ESMF_Geom) :: geom + type(ESMF_FieldBundle) :: bundle + type(ESMF_Field), allocatable :: fieldList(:) + real(kind=ESMF_KIND_R8), pointer :: x_r8(:,:) + real(kind=ESMF_KIND_R4), allocatable :: weights(:) + integer :: i + + call setup_geom(geom, 4) + call setup_bundle(bundle, weights=DEFAULT_WEIGHTS, geom=geom, typekind=ESMF_TYPEKIND_R4, units='m') + + delta = FieldBundleDelta(FieldDelta(typekind=ESMF_TYPEKIND_R8)) + call delta%update_bundle(bundle, _RC) + + call MAPL_FieldBundleGet(bundle, fieldList=fieldList, _RC) + @assert_that(size(fieldList), is(FIELD_COUNT)) + + do i = 1, FIELD_COUNT + call ESMF_FieldGet(fieldList(i), fArrayPtr=x_r8, _RC) + @assert_that(shape(x_r8), is(equal_to([4,4]))) + end do + + call MAPL_InfoGetInternal(bundle, KEY_INTERPOLATION_WEIGHTS, values=weights, _RC) + @assert_that(weights, is(equal_to(DEFAULT_WEIGHTS))) + + call teardown_bundle(bundle) + call teardown_geom(geom) + + _UNUSED_DUMMY(this) + end subroutine test_change_typekind + + @test(type=ESMF_TestMethod, npes=[1]) + subroutine test_change_units(this) + class(ESMF_TestMethod), intent(inout) :: this + + integer :: status + type(Fieldbundledelta) :: delta + type(ESMF_Geom) :: geom + type(ESMF_FieldBundle) :: bundle + integer :: i + type(ESMF_Field), allocatable :: fieldList(:) + real(kind=ESMF_KIND_R4), pointer :: x_r4(:,:) + character(:), allocatable :: new_units + + call setup_geom(geom, 4) + call setup_bundle(bundle, weights=DEFAULT_WEIGHTS, geom=geom, typekind=ESMF_TYPEKIND_R4, units='km') + + delta = FieldBundleDelta(FieldDelta(units='m')) + call delta%update_bundle(bundle, _RC) ! must reallocate fields + + call MAPL_FieldBundleGet(bundle, fieldList=fieldList, _RC) + @assert_that(size(fieldList), is(FIELD_COUNT)) + + do i = 1, FIELD_COUNT + call ESMF_FieldGet(fieldList(i), fArrayPtr=x_r4, _RC) + @assert_that(shape(x_r4), is(equal_to([4,4]))) + @assert_that(x_r4, every_item(is(FILL_VALUE))) + + call MAPL_infoGetInternal(fieldList(i), KEY_UNITS, value=new_units, _RC) + @assertEqual('m', new_units) + end do + + call teardown_bundle(bundle) + call teardown_geom(geom) + + _UNUSED_DUMMY(this) + end subroutine test_change_units + + @test(type=ESMF_TestMethod, npes=[1]) + subroutine test_change_geom(this) + class(ESMF_TestMethod), intent(inout) :: this + + integer :: status + type(Fieldbundledelta) :: delta + type(ESMF_Geom) :: geom, new_geom, tmp_geom + type(ESMF_FieldBundle) :: bundle + integer :: i + type(ESMF_Field), allocatable :: fieldList(:) + real(kind=ESMF_KIND_R4), pointer :: x_r4(:,:) + character(:), allocatable :: new_units + + call setup_geom(geom, 4) + call setup_bundle(bundle, weights=DEFAULT_WEIGHTS, geom=geom, typekind=ESMF_TYPEKIND_R4, units='km') + + call setup_geom(new_geom, 6) + delta = FieldBundleDelta(FieldDelta(new_geom)) ! same geom + call delta%update_bundle(bundle, _RC) ! should reallocate fields + + call MAPL_FieldBundleGet(bundle, fieldList=fieldList, _RC) + @assert_that(size(fieldList), is(FIELD_COUNT)) + + do i = 1, FIELD_COUNT + call ESMF_FieldGet(fieldList(i), fArrayPtr=x_r4, _RC) + @assert_that(shape(x_r4), is(equal_to([6,6]))) + + call MAPL_InfoGetInternal(fieldList(i), KEY_UNITS, value=new_units, _RC) + @assertEqual('km', new_units) + + call ESMF_FieldGet(fieldList(i), geom=tmp_geom, _RC) + @assert_that(tmp_geom == new_geom, is(true())) + end do + + call teardown_bundle(bundle) + call teardown_geom(geom) + + _UNUSED_DUMMY(this) + end subroutine test_change_geom + + + @test(type=ESMF_TestMethod, npes=[1]) + subroutine test_same_geom(this) + class(ESMF_TestMethod), intent(inout) :: this + + integer :: status + type(Fieldbundledelta) :: delta + type(ESMF_Geom) :: geom, tmp_geom + type(ESMF_FieldBundle) :: bundle + integer :: i + type(ESMF_Field), allocatable :: fieldList(:) + real(kind=ESMF_KIND_R4), pointer :: x_r4(:,:) + character(:), allocatable :: new_units + + call setup_geom(geom, 4) + call setup_bundle(bundle, weights=DEFAULT_WEIGHTS, geom=geom, typekind=ESMF_TYPEKIND_R4, units='km') + + delta = FieldBundleDelta(FieldDelta(geom)) ! same geom + call delta%update_bundle(bundle, _RC) ! should not reallocate fields + + call MAPL_FieldBundleGet(bundle, fieldList=fieldList, _RC) + @assert_that(size(fieldList), is(FIELD_COUNT)) + + do i = 1, FIELD_COUNT + call ESMF_FieldGet(fieldList(i), fArrayPtr=x_r4, _RC) + @assert_that(shape(x_r4), is(equal_to([4,4]))) + @assert_that(x_r4, every_item(is(FILL_VALUE))) + + call MAPL_InfoGetInternal(fieldList(i), KEY_UNITS, value=new_units, _RC) + @assertEqual('km', new_units) + + call ESMF_FieldGet(fieldList(i), geom=tmp_geom, _RC) + @assert_that(tmp_geom == geom, is(true())) + end do + + call teardown_bundle(bundle) + call teardown_geom(geom) + + _UNUSED_DUMMY(this) + end subroutine test_same_geom + + @test(type=ESMF_TestMethod, npes=[1]) + subroutine test_change_weights(this) + class(ESMF_TestMethod), intent(inout) :: this + + integer :: status + type(Fieldbundledelta) :: delta + type(ESMF_Geom) :: geom, tmp_geom + type(ESMF_FieldBundle) :: bundle + integer :: i + type(ESMF_Field), allocatable :: fieldList(:) + real(kind=ESMF_KIND_R4), pointer :: x_r4(:,:) + character(:), allocatable :: new_units + real(kind=ESMF_KIND_R4), allocatable :: weights(:) + real(kind=ESMF_KIND_R4), parameter :: new_weights(*) = [0.,0.25,0.75] + + call setup_geom(geom, 4) + call setup_bundle(bundle, weights=DEFAULT_WEIGHTS, geom=geom, typekind=ESMF_TYPEKIND_R4, units='km') + + + delta = FieldBundleDelta(interpolation_weights=new_weights) + call delta%update_bundle(bundle, _RC) ! should not reallocate fields + + call MAPL_FieldBundleGet(bundle, fieldList=fieldList, _RC) + @assert_that(size(fieldList), is(FIELD_COUNT)) + + do i = 1, FIELD_COUNT + call ESMF_FieldGet(fieldList(i), fArrayPtr=x_r4, _RC) + @assert_that(shape(x_r4), is(equal_to([4,4]))) + @assert_that(x_r4, every_item(is(FILL_VALUE))) + + call MAPL_InfoGetInternal(fieldList(i), KEY_UNITS, value=new_units, _RC) + @assertEqual('km', new_units) + + call ESMF_FieldGet(fieldList(i), geom=tmp_geom, _RC) + @assert_that(tmp_geom == geom, is(true())) + end do + + call MAPL_InfoGetInternal(bundle, KEY_INTERPOLATION_WEIGHTS, values=weights, _RC) + @assert_that(weights, is(equal_to(new_weights))) + + call teardown_bundle(bundle) + call teardown_geom(geom) + + _UNUSED_DUMMY(this) + end subroutine test_change_weights + + @test(type=ESMF_TestMethod, npes=[1]) + subroutine test_change_weights_with_ungridded(this) + class(ESMF_TestMethod), intent(inout) :: this + + integer :: status + type(Fieldbundledelta) :: delta + type(ESMF_Geom) :: geom, tmp_geom + type(ESMF_FieldBundle) :: bundle + integer :: i + type(ESMF_Field), allocatable :: fieldList(:) + real(kind=ESMF_KIND_R4), pointer :: x_r4(:,:,:,:) + character(:), allocatable :: new_units + real(kind=ESMF_KIND_R4), allocatable :: weights(:) + real(kind=ESMF_KIND_R4), parameter :: new_weights(*) = [0.,0.25,0.75] + integer :: ndims, nlevels, rank + + call setup_geom(geom, 4) + call setup_bundle(bundle, weights=DEFAULT_WEIGHTS, geom=geom, typekind=ESMF_TYPEKIND_R4, units='km', with_ungridded=.true.) + + + delta = FieldBundleDelta(interpolation_weights=new_weights) + call delta%update_bundle(bundle, _RC) ! should not reallocate fields + + call MAPL_FieldBundleGet(bundle, fieldList=fieldList, _RC) + @assert_that(size(fieldList), is(FIELD_COUNT)) + + do i = 1, FIELD_COUNT + call ESMF_FieldGet(fieldList(i), rank=rank, _RC) + + call ESMF_FieldGet(fieldList(i), fArrayPtr=x_r4, _RC) + @assert_that(shape(x_r4), is(equal_to([4,4,NUM_LEVELS,NUM_RADII]))) + @assert_that(all(x_r4 == FILL_VALUE), is(true())) + + call MAPL_InfoGetInternal(fieldList(i), KEY_UNITS, value=new_units, _RC) + @assertEqual('km', new_units) + + call ESMF_FieldGet(fieldList(i), geom=tmp_geom, _RC) + @assert_that(tmp_geom == geom, is(true())) + + call MAPL_InfoGetInternal(fieldList(i), KEY_UNGRIDDED_DIMS//KEY_NUM_UNGRIDDED_DIMS, value=ndims, _RC) + @assert_that(ndims, is(1)) + + call MAPL_InfoGetInternal(fieldList(i), KEY_NUM_LEVELS, value=nlevels, _RC) + @assert_that(nlevels, is(NUM_LEVELS)) + + end do + + call MAPL_InfoGetInternal(bundle, KEY_INTERPOLATION_WEIGHTS, values=weights, _RC) + @assert_that(weights, is(equal_to(new_weights))) + + call MAPL_InfoGetInternal(bundle, KEY_UNGRIDDED_DIMS//KEY_NUM_UNGRIDDED_DIMS, value=ndims, _RC) + @assert_that(ndims, is(1)) + + call MAPL_InfoGetInternal(bundle, KEY_NUM_LEVELS, value=nlevels, _RC) + @assert_that(nlevels, is(NUM_LEVELS)) + + call teardown_bundle(bundle) + call teardown_geom(geom) + + _UNUSED_DUMMY(this) + end subroutine test_change_weights_with_ungridded + + @test(type=ESMF_TestMethod, npes=[1]) + ! This is the hard use case. Typically it arises when ExtData + ! starts with a rule which is a constant expression, but then later + ! becomes an ordinary interpolation rule. The bundle then goes + ! from 0 fields to 2 fields. The hard part is finding all the information that + ! is needed to create properly initialized fields. E.g., geom, units, ... + subroutine test_create_fields(this) + class(ESMF_TestMethod), intent(inout) :: this + + integer :: status + type(Fieldbundledelta) :: delta + type(ESMF_Geom) :: geom, tmp_geom + type(ESMF_FieldBundle) :: bundle + integer :: i + type(ESMF_Field), allocatable :: fieldList(:) + real(kind=ESMF_KIND_R4), pointer :: x_r4(:,:) + character(:), allocatable :: new_units + real(kind=ESMF_KIND_R4), allocatable :: weights(:) + real(kind=ESMF_KIND_R4), parameter :: new_weights(*) = [0.,0.25,0.75] + + call setup_geom(geom, 4) + call setup_bundle(bundle, weights=[5.], geom=geom, typekind=ESMF_TYPEKIND_R4, units='km') + + delta = FieldBundleDelta(interpolation_weights=new_weights) + call delta%update_bundle(bundle, _RC) ! should allocate fields + + call MAPL_FieldBundleGet(bundle, fieldList=fieldList, _RC) + @assert_that(size(fieldList), is(FIELD_COUNT)) + + do i = 1, FIELD_COUNT + call ESMF_FieldGet(fieldList(i), fArrayPtr=x_r4, _RC) + @assert_that(shape(x_r4), is(equal_to([4,4]))) + + call MAPL_InfoGetInternal(fieldList(i), KEY_UNITS, value=new_units, _RC) + @assertEqual('km', new_units) + + call ESMF_FieldGet(fieldList(i), geom=tmp_geom, _RC) + @assert_that(tmp_geom == geom, is(true())) + end do + + call MAPL_InfoGetInternal(bundle, KEY_INTERPOLATION_WEIGHTS, values=weights, _RC) + @assert_that(weights, is(equal_to(new_weights))) + + call teardown_bundle(bundle) + call teardown_geom(geom) + + _UNUSED_DUMMY(this) + end subroutine test_create_fields + + @test(type=ESMF_TestMethod, npes=[1]) + ! This is the hard use case. Typically it arises when ExtData + ! starts with a rule which is a constant expression, but then later + ! becomes an ordinary interpolation rule. The bundle then goes + ! from 0 fields to 2 fields. The hard part is finding all the information that + ! is needed to create properly initialized fields. E.g., geom, units, ... + subroutine test_create_fields_with_ungridded(this) + class(ESMF_TestMethod), intent(inout) :: this + + integer :: status + type(Fieldbundledelta) :: delta + type(ESMF_Geom) :: geom, tmp_geom + type(ESMF_FieldBundle) :: bundle + integer :: i + type(ESMF_Field), allocatable :: fieldList(:) + real(kind=ESMF_KIND_R4), pointer :: x_r4(:,:,:,:) + character(:), allocatable :: new_units + real(kind=ESMF_KIND_R4), allocatable :: weights(:) + real(kind=ESMF_KIND_R4), parameter :: new_weights(*) = [0.,0.25,0.75] + integer :: ndims, nlevels + + call setup_geom(geom, 4) + call setup_bundle(bundle, weights=[5.], geom=geom, typekind=ESMF_TYPEKIND_R4, units='km', & + with_ungridded=.true.) + + delta = FieldBundleDelta(interpolation_weights=new_weights) + call delta%update_bundle(bundle, _RC) ! should allocate fields + + call MAPL_FieldBundleGet(bundle, fieldList=fieldList, _RC) + @assert_that(size(fieldList), is(FIELD_COUNT)) + + do i = 1, FIELD_COUNT + call ESMF_FieldGet(fieldList(i), fArrayPtr=x_r4, _RC) + @assert_that(shape(x_r4), is(equal_to([4,4,NUM_LEVELS,NUM_RADII]))) + + call MAPL_InfoGetInternal(fieldList(i), KEY_UNITS, value=new_units, _RC) + @assertEqual('km', new_units) + + call ESMF_FieldGet(fieldList(i), geom=tmp_geom, _RC) + @assert_that(tmp_geom == geom, is(true())) + + call MAPL_InfoGetInternal(fieldList(i), KEY_UNGRIDDED_DIMS//KEY_NUM_UNGRIDDED_DIMS, value=ndims, _RC) + @assert_that(ndims, is(1)) + + call MAPL_InfoGetInternal(fieldList(i), KEY_NUM_LEVELS, value=nlevels, _RC) + @assert_that(nlevels, is(NUM_LEVELS)) + end do + + call MAPL_InfoGetInternal(bundle, KEY_INTERPOLATION_WEIGHTS, values=weights, _RC) + @assert_that(weights, is(equal_to(new_weights))) + + call MAPL_InfoGetInternal(bundle, KEY_UNGRIDDED_DIMS//KEY_NUM_UNGRIDDED_DIMS, value=ndims, _RC) + @assert_that(ndims, is(1)) + + call MAPL_InfoGetInternal(bundle, KEY_NUM_LEVELS, value=nlevels, _RC) + @assert_that(nlevels, is(NUM_LEVELS)) + + call teardown_bundle(bundle) + call teardown_geom(geom) + + _UNUSED_DUMMY(this) + end subroutine test_create_fields_with_ungridded + +end module Test_FieldBundleDelta + diff --git a/field_utils/tests/Test_FieldUtilities.pf b/field_utils/tests/Test_FieldDelta.pf similarity index 79% rename from field_utils/tests/Test_FieldUtilities.pf rename to field_utils/tests/Test_FieldDelta.pf index 4fcdfabcd329..9a58684634a6 100644 --- a/field_utils/tests/Test_FieldUtilities.pf +++ b/field_utils/tests/Test_FieldDelta.pf @@ -1,13 +1,13 @@ #include "MAPL_TestErr.h" #include "unused_dummy.H" -module Test_FieldUtilities - use mapl_FieldUtilities +module Test_FieldDelta + use mapl3g_FieldDelta use mapl3g_ESMF_Info_Keys use mapl3g_InfoUtilities use esmf use ESMF_TestMethod_mod use funit - implicit none + implicit none (type, external) integer, parameter :: ORIGINAL_NUM_LEVELS = 5 real, parameter :: FILL_VALUE = 99. @@ -26,6 +26,7 @@ contains integer :: status type(ESMF_FieldStatus_Flag) :: field_status type(ESMF_TypeKind_Flag) :: typekind + type(FieldDelta) :: delta grid = ESMF_GridCreateNoPeriDim(maxIndex=[4,4], name='I_AM_GROOT', _RC) geom = ESMF_GeomCreate(grid, _RC) @@ -34,7 +35,8 @@ contains call MAPL_InfoSetInternal(f, key=KEY_NUM_LEVELS, value=ORIGINAL_NUM_LEVELS, _RC) call MAPL_InfoSetInternal(f, key=KEY_VLOC, value='VERTICAL_DIM_EDGE', _RC) - call FieldReallocate(f, typekind=ESMF_TYPEKIND_R8, _RC) + delta = FieldDelta(typekind=ESMF_TYPEKIND_R8) + call delta%reallocate_field(f, _RC) call ESMF_FieldGet(f, status=field_status, typekind=typekind, _RC) @assert_that(field_status == ESMF_FIELDSTATUS_COMPLETE, is(true())) @@ -58,6 +60,7 @@ contains type(ESMF_FieldStatus_Flag) :: field_status type(ESMF_TypeKind_Flag) :: typekind real(kind=ESMF_KIND_R4), pointer :: x(:,:) + type(FieldDelta) :: delta grid = ESMF_GridCreateNoPeriDim(maxIndex=[4,4], name='I_AM_GROOT', _RC) geom = ESMF_GeomCreate(grid, _RC) @@ -68,7 +71,8 @@ contains call ESMF_FieldGet(f, fArrayPtr=x, _RC) x = FILL_VALUE - call FieldReallocate(f, typekind=ESMF_TYPEKIND_R4, _RC) + delta = FieldDelta(typekind=ESMF_TYPEKIND_R4) + call delta%reallocate_field(f, _RC) call ESMF_FieldGet(f, status=field_status, typekind=typekind, geom=other_geom, _RC) @assert_that(field_status == ESMF_FIELDSTATUS_COMPLETE, is(true())) @@ -96,6 +100,7 @@ contains type(ESMF_FieldStatus_Flag) :: field_status type(ESMF_TypeKind_Flag) :: typekind real(kind=ESMF_KIND_R4), pointer :: x(:,:) + type(FieldDelta) :: delta grid1 = ESMF_GridCreateNoPeriDim(maxIndex=[4,4], name='I_AM_GROOT', _RC) geom1 = ESMF_GeomCreate(grid1, _RC) @@ -105,7 +110,8 @@ contains grid2 = ESMF_GridCreateNoPeriDim(maxIndex=[3,5], name='I_AM_GROOT', _RC) geom2 = ESMF_GeomCreate(grid2, _RC) - call FieldReallocate(f, geom=geom2, _RC) ! same geom + delta = FieldDelta(geom=geom2) + call delta%reallocate_field(f, _RC) call ESMF_FieldGet(f, status=field_status, typekind=typekind, _RC) @assert_that(field_status == ESMF_FIELDSTATUS_COMPLETE, is(true())) @@ -134,7 +140,8 @@ contains type(ESMF_FieldStatus_Flag) :: field_status type(ESMF_TypeKind_Flag) :: typekind real(kind=ESMF_KIND_R4), pointer :: x(:,:) - + type(FieldDelta) :: delta + grid1 = ESMF_GridCreateNoPeriDim(maxIndex=[4,4], name='I_AM_GROOT', _RC) geom1 = ESMF_GeomCreate(grid1, _RC) f = ESMF_FieldCreate(geom1, typekind=ESMF_TYPEKIND_R4, name='in', _RC) @@ -151,8 +158,9 @@ contains x = FILL_VALUE geom2 = geom1 - call FieldReallocate(f, geom=geom2, _RC) ! same geom - + delta = FieldDelta(geom=geom2) + call delta%reallocate_field(f, _RC) + call ESMF_FieldGet(f, status=field_status, typekind=typekind, _RC) @assert_that(field_status == ESMF_FIELDSTATUS_COMPLETE, is(true())) @assert_that(typekind == ESMF_TYPEKIND_R4, is(true())) @@ -181,31 +189,77 @@ contains type(ESMF_FieldStatus_Flag) :: field_status type(ESMF_TypeKind_Flag) :: typekind real(ESMF_KIND_R4), pointer :: x(:,:,:,:) + type(FieldDelta) :: delta grid = ESMF_GridCreateNoPeriDim(maxIndex=[4,4], name='I_AM_GROOT', _RC) geom = ESMF_GeomCreate(grid, _RC) f = ESMF_FieldCreate(geom, typekind=ESMF_TYPEKIND_R4, name='in', & - ungriddedLbound=[1,1], ungriddedUbound=[ORIGINAL_NUM_LEVELS,3], _RC) + ungriddedLbound=[1,1], ungriddedUbound=[ORIGINAL_NUM_LEVELS+1,3], _RC) call MAPL_InfoSetInternal(f, key=KEY_NUM_LEVELS, value=ORIGINAL_NUM_LEVELS, _RC) call MAPL_InfoSetInternal(f, key=KEY_VLOC, value='VERTICAL_DIM_EDGE', _RC) + delta = FieldDelta(num_levels=4) + call delta%reallocate_field(f, _RC) + + call ESMF_FieldGet(f, status=field_status, typekind=typekind, _RC) + @assert_that(field_status == ESMF_FIELDSTATUS_COMPLETE, is(true())) + @assert_that(typekind == ESMF_TYPEKIND_R4, is(true())) - call FieldReallocate(f, num_levels=4, _RC) + call ESMF_FieldGet(f, fArrayPtr=x, _RC) + @assert_that(shape(x), is(equal_to([4,4,4+1,3]))) + + call ESMF_FieldDestroy(f, _RC) + call ESMF_GridDestroy(grid, _RC) + call ESMF_GeomDestroy(geom, _RC) + + _UNUSED_DUMMY(this) + end subroutine test_change_n_levels + + @test(type=ESMF_TestMethod, npes=[1]) + ! Surface fields should be unaffected when changing num_levels of + ! vertical grid. + subroutine test_change_n_levels_surface(this) + class(ESMF_TestMethod), intent(inout) :: this + type(ESMF_Field) :: f + type(ESMF_Grid) :: grid + type(ESMF_Geom) :: geom + + integer :: status + type(ESMF_FieldStatus_Flag) :: field_status + type(ESMF_TypeKind_Flag) :: typekind + real(ESMF_KIND_R4), pointer :: x(:,:,:,:) + type(FieldDelta) :: delta + + grid = ESMF_GridCreateNoPeriDim(maxIndex=[4,4], name='I_AM_GROOT', _RC) + geom = ESMF_GeomCreate(grid, _RC) + + ! Surface field + f = ESMF_FieldCreate(geom, typekind=ESMF_TYPEKIND_R4, name='in', & + ungriddedLbound=[1,1], ungriddedUbound=[2,3], _RC) + call MAPL_InfoSetInternal(f, key=KEY_NUM_LEVELS, value=ORIGINAL_NUM_LEVELS, _RC) + call MAPL_InfoSetInternal(f, key=KEY_VLOC, value='VERTICAL_DIM_NONE', _RC) + call ESMF_FieldGet(f, fArrayPtr=x, _RC) + x = FILL_VALUE + + delta = FieldDelta(num_levels=4) + call delta%reallocate_field(f, _RC) call ESMF_FieldGet(f, status=field_status, typekind=typekind, _RC) @assert_that(field_status == ESMF_FIELDSTATUS_COMPLETE, is(true())) @assert_that(typekind == ESMF_TYPEKIND_R4, is(true())) call ESMF_FieldGet(f, fArrayPtr=x, _RC) - @assert_that(shape(x), is(equal_to([4,4,4,3]))) + @assert_that(shape(x), is(equal_to([4,4,2,3]))) + @assert_that(all(x == FILL_VALUE), is(true())) call ESMF_FieldDestroy(f, _RC) call ESMF_GridDestroy(grid, _RC) call ESMF_GeomDestroy(geom, _RC) _UNUSED_DUMMY(this) - end subroutine test_change_n_levels + end subroutine test_change_n_levels_surface + @test(type=ESMF_TestMethod, npes=[1]) @@ -219,19 +273,21 @@ contains type(ESMF_FieldStatus_Flag) :: field_status real(ESMF_KIND_R4), pointer :: x(:,:,:,:) type(ESMF_TypeKind_Flag) :: typekind + type(FieldDelta) :: delta grid = ESMF_GridCreateNoPeriDim(maxIndex=[4,4], name='I_AM_GROOT', _RC) geom = ESMF_GeomCreate(grid, _RC) f = ESMF_FieldCreate(geom, typekind=ESMF_TYPEKIND_R4, name='in', & - ungriddedLbound=[1,1], ungriddedUbound=[ORIGINAL_NUM_LEVELS,3], _RC) + ungriddedLbound=[1,1], ungriddedUbound=[ORIGINAL_NUM_LEVELS+1,3], _RC) call MAPL_InfoSetInternal(f, key=KEY_NUM_LEVELS, value=ORIGINAL_NUM_LEVELS, _RC) call MAPL_InfoSetInternal(f, key=KEY_VLOC, value='VERTICAL_DIM_EDGE', _RC) call ESMF_FieldGet(f, fArrayPtr=x, _RC) x = FILL_VALUE - call FieldReallocate(f, num_levels=ORIGINAL_NUM_LEVELS, _RC) + delta = FieldDelta(num_levels=ORIGINAL_NUM_LEVELS) + call delta%reallocate_field(f, _RC) call ESMF_FieldGet(f, status=field_status, typekind=typekind, geom=other_geom, _RC) @assert_that(field_status == ESMF_FIELDSTATUS_COMPLETE, is(true())) @@ -262,6 +318,7 @@ contains type(ESMF_FieldStatus_Flag) :: field_status real(ESMF_KIND_R8), pointer :: x8(:,:,:,:) type(ESMF_TypeKind_Flag) :: typekind + type(FieldDelta) :: delta grid = ESMF_GridCreateNoPeriDim(maxIndex=[4,4], name='I_AM_GROOT', _RC) geom = ESMF_GeomCreate(grid, _RC) @@ -270,20 +327,21 @@ contains geom_ref = ESMF_GeomCreate(grid_ref, _RC) f = ESMF_FieldCreate(geom, typekind=ESMF_TYPEKIND_R4, name='in', & - ungriddedLbound=[1,1], ungriddedUbound=[ORIGINAL_NUM_LEVELS+1,3], _RC) + ungriddedLbound=[1,1], ungriddedUbound=[ORIGINAL_NUM_LEVELS,3], _RC) f_ref = ESMF_FieldCreate(geom_ref, typekind=ESMF_TYPEKIND_R8, name='in', & - ungriddedLbound=[1,1], ungriddedUbound=[ORIGINAL_NUM_LEVELS,3], _RC) + ungriddedLbound=[1,1], ungriddedUbound=[ORIGINAL_NUM_LEVELS-1,3], _RC) - call MAPL_InfoSetInternal(f, key=KEY_NUM_LEVELS, value=ORIGINAL_NUM_LEVELS+1, _RC) - call MAPL_InfoSetInternal(f, key=KEY_VLOC, value='VERTICAL_DIM_EDGE', _RC) + call MAPL_InfoSetInternal(f, key=KEY_NUM_LEVELS, value=ORIGINAL_NUM_LEVELS, _RC) + call MAPL_InfoSetInternal(f, key=KEY_VLOC, value='VERTICAL_DIM_CENTER', _RC) call MAPL_InfoSetInternal(f, key=KEY_UNITS, value=ORIGINAL_UNITS) - call MAPL_InfoSetInternal(f_ref, key=KEY_NUM_LEVELS, value=ORIGINAL_NUM_LEVELS, _RC) - call MAPL_InfoSetInternal(f_ref, key=KEY_VLOC, value='VERTICAL_DIM_EDGE', _RC) + call MAPL_InfoSetInternal(f_ref, key=KEY_NUM_LEVELS, value=ORIGINAL_NUM_LEVELS-1, _RC) + call MAPL_InfoSetInternal(f_ref, key=KEY_VLOC, value='VERTICAL_DIM_CENTER', _RC) call MAPL_InfoSetInternal(f_ref, key=KEY_UNITS, value=REFERENCE_UNITS) - call FieldUpdate(f, f_ref, ignore='geom', _RC) + call delta%initialize(f, f_ref, _RC) + call delta%update_field(f, ignore='geom', _RC) call ESMF_FieldGet(f, status=field_status, typekind=typekind, geom=new_geom, _RC) @assert_that(field_status == ESMF_FIELDSTATUS_COMPLETE, is(true())) @@ -295,7 +353,7 @@ contains ! check that field shape is changed due to new num levels call ESMF_FieldGet(f, fArrayPtr=x8, _RC) - @assert_that(shape(x8),is(equal_to([4,4,ORIGINAL_NUM_LEVELS,3]))) + @assert_that(shape(x8),is(equal_to([4,4,ORIGINAL_NUM_LEVELS-1,3]))) call ESMF_FieldDestroy(f, _RC) call ESMF_GridDestroy(grid, _RC) @@ -304,5 +362,4 @@ contains _UNUSED_DUMMY(this) end subroutine test_field_update_from_field_ignore_geom - -end module Test_FieldUtilities +end module Test_FieldDelta diff --git a/generic3g/actions/ExtensionAction.F90 b/generic3g/actions/ExtensionAction.F90 index 0ee10ddcce73..3c52fa3d4d76 100644 --- a/generic3g/actions/ExtensionAction.F90 +++ b/generic3g/actions/ExtensionAction.F90 @@ -1,6 +1,7 @@ #include "MAPL_Generic.h" module mapl3g_ExtensionAction use mapl_ErrorHandling + use ESMF implicit none private diff --git a/generic3g/couplers/BidirectionalObserver.F90 b/generic3g/couplers/BidirectionalObserver.F90 deleted file mode 100644 index d982438d701a..000000000000 --- a/generic3g/couplers/BidirectionalObserver.F90 +++ /dev/null @@ -1,107 +0,0 @@ -#include "MAPL_Generic.h" - -module mapl3g_BidirectionalObserver - use mapl3g_Observer - use mapl_ErrorHandlingMod - implicit none - private - - ! Class - public :: BidirectionalObserver - - - ! Ideally this will not be abstract, but for now it is - type, extends(Observer), abstract :: BidirectionalObserver - private - type(ObserverPtrVector) :: import_observers ! think couplers - type(ObserverPtrVector) :: export_observers ! think couplers - contains - procedure :: update - procedure :: invalidate - procedure :: update_imports - procedure :: invalidate_exports - end type BidirectionalObserver - - abstract interface - subroutine I_Notify(this, rc) - import :: BidirectionalObserver - class(Obserer), intent(inout) :: this - integer, optional, intent(out) :: rc - end subroutine I_Notify - end interface - -contains - - recursive function update(this, rc) - class(Observable), intent(inout) :: this - integer, optional, intent(out) :: rc - - integer :: status - logical :: is_up_to_date - - is_up_to_date = this%is_up_to_date() - _RETURN_IF(is_up_to_date) - - call this%update_imports(_RC) - call this%update_self(_RC) - - _RETURN(_SUCCESS) - end function update - - recursive function invalidate(this, rc) - class(Observable), intent(inout) :: this - integer, optional, intent(out) :: rc - - integer :: status - logical :: is_stale - - is_stale = this%is_up_to_date() - _RETURN_IF(is_up_to_date) - - call this%invalidate_self(_RC) - call this%invalidate_exports(_RC) - - _RETURN(_SUCCESS) - end function invalidate - - - recursive subroutine update_imports(this, rc) - class(BidirectionalObserver), intent(inout) :: this - integer, optional, intent(out) :: rc - - integer :: status - type(ObserverPtrVectorIterator) :: iter - class(ObserverPtr), pointer :: obsrvr - - associate(e => this%import_observers%ftn_end()) - iter = observers%ftn_begin() - do while (iter /= e) - call iter%next() - obsrvr => iter%of() - call obsrvr%ptr%update(_RC) - end do - end associate - - _RETURN(_SUCCESS) - end subroutine update_imports - - subroutine invalidate_exports(observers, rc) - class(BidirectionalObserver), intent(inout) :: observers - integer, optional, intent(out) :: rc - - integer :: status - - associate(e => this%export_observers%ftn_end()) - iter = observers%ftn_begin() - do while (iter /= e) - call iter%next() - obsrvr => iter%of() - call obsrvr%ptr%invalidate(_RC) - end do - end associate - - - _RETURN(_SUCCESS) - end subroutine invalidate_exports - -end module mapl3g_BidirectionalObserver diff --git a/generic3g/couplers/CouplerMetaComponent.F90 b/generic3g/couplers/CouplerMetaComponent.F90 index f26f57c49d23..a8b086c30b69 100644 --- a/generic3g/couplers/CouplerMetaComponent.F90 +++ b/generic3g/couplers/CouplerMetaComponent.F90 @@ -91,6 +91,27 @@ recursive subroutine initialize(this, importState, exportState, clock, rc) _RETURN(_SUCCESS) end subroutine initialize + ! Check if export item has been updated and update import item + ! accordingly. + recursive subroutine update_time_varying(this, importState, exportState, clock, rc) + class(CouplerMetaComponent), intent(inout) :: this + type(ESMF_State), intent(inout) :: importState + type(ESMF_State), intent(inout) :: exportState + type(ESMF_Clock), intent(inout) :: clock + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_Field) :: f_in, f_out + +!# _RETURN_UNLESS(this%export_is_time_varying()) + call ESMF_StateGet(importState, itemName='import[1]', field=f_in, _RC) + call ESMF_StateGet(exportState, itemName='export[1]', field=f_out, _RC) + +!# call FieldUpdate(f_in, from=f_out, ignore=this%action%get_ignore(), _RC) + + _RETURN(_SUCCESS) + end subroutine update_time_varying + recursive subroutine update(this, importState, exportState, clock, rc) class(CouplerMetaComponent), intent(inout) :: this type(ESMF_State), intent(inout) :: importState @@ -126,6 +147,27 @@ recursive subroutine update_sources(this, rc) _RETURN(_SUCCESS) end subroutine update_sources + ! Check if export item has been updated and update import item + ! accordingly. + recursive subroutine invalidate_time_varying(this, importState, exportState, clock, rc) + class(CouplerMetaComponent), intent(inout) :: this + type(ESMF_State), intent(inout) :: importState + type(ESMF_State), intent(inout) :: exportState + type(ESMF_Clock), intent(inout) :: clock + integer, optional, intent(out) :: rc + + integer :: status + type(ESMF_Field) :: f_in, f_out + +!# _RETURN_UNLESS(this%import_is_time_varying()) + call ESMF_StateGet(importState, itemName='import[1]', field=f_in, _RC) + call ESMF_StateGet(exportState, itemName='export[1]', field=f_out, _RC) + +!# call FieldUpdate(f_out, from=f_in, ignore=this%action%get_ignore(), _RC) + + _RETURN(_SUCCESS) + end subroutine invalidate_time_varying + recursive subroutine invalidate(this, importState, exportState, clock, rc) class(CouplerMetaComponent) :: this type(ESMF_State) :: importState diff --git a/generic3g/couplers/GenericCoupler.F90 b/generic3g/couplers/GenericCoupler.F90 index a7c9f3017db1..7859c29169a4 100644 --- a/generic3g/couplers/GenericCoupler.F90 +++ b/generic3g/couplers/GenericCoupler.F90 @@ -98,6 +98,7 @@ recursive subroutine update(gridcomp, importState, exportState, clock, rc) type(CouplerMetaComponent), pointer :: meta meta => get_coupler_meta(gridcomp, _RC) +!# call meta%update_time_varying(importState, exportState, _RC) call meta%update(importState, exportState, clock, _RC) _RETURN(_SUCCESS) @@ -115,6 +116,7 @@ recursive subroutine invalidate(gridcomp, importState, exportState, clock, rc) type(CouplerMetaComponent), pointer :: meta meta => get_coupler_meta(gridcomp, _RC) +!# call meta%invalidate_time_varying(importState, exportState, _RC) call meta%invalidate(importstate, exportState, clock, _RC) _RETURN(_SUCCESS) diff --git a/generic3g/couplers/HandlerMap.F90 b/generic3g/couplers/HandlerMap.F90 deleted file mode 100644 index 1c53a53c7fba..000000000000 --- a/generic3g/couplers/HandlerMap.F90 +++ /dev/null @@ -1,20 +0,0 @@ -module mapl3g_ComponentHandlerMap - use mapl3g_AbstractComponentHandler - ! Maybe should be VirtualConnectionPt instead? -#define Key __CHARACTER_DEFERRED -#define T AbstractComponentHandler -#define T_polymorphic -#define Map ComponentHandlerMap -#define MapIterator ComponentHandlerMapIterator -#define Pair ComponentHandlerPair - -#include "map/template.inc" - -#undef Pair -#undef MapIterator -#undef Map -#undef T_polymorphic -#undef T -#undef Key - -end module mapl3g_CouplerComponentVector diff --git a/generic3g/couplers/HandlerVector.F90 b/generic3g/couplers/HandlerVector.F90 deleted file mode 100644 index 5f73b6f48f9d..000000000000 --- a/generic3g/couplers/HandlerVector.F90 +++ /dev/null @@ -1,16 +0,0 @@ -module mapl3g_ComponentHandlerVector - use mapl3g_AbstractComponentHandler - -#define T AbstractComponentHandler -#define T_polymorphic -#define Vector ComponentHandlerVector -#define VectorIterator ComponentHandlerVectorIterator - -#include "vector/template.inc" - -#undef VectorIterator -#undef Vector -#undef T_polymorphic -#undef T - -end module mapl3g_ComponentHandlerVector diff --git a/generic3g/couplers/ImportCoupler.F90 b/generic3g/couplers/ImportCoupler.F90 deleted file mode 100644 index 66f230d910b9..000000000000 --- a/generic3g/couplers/ImportCoupler.F90 +++ /dev/null @@ -1,25 +0,0 @@ -module mapl3g_ImportCoupler - use mapl3g_GenericCoupler - implicit none - private - - public :: ImportCoupler - - type, extends :: GenericCoupler - contains - procedure :: update - end type GenericCoupler - -contains - - subroutine update(this) - class(ImportCoupler), intent(in) :: this - - alarm = ESMF_ClockGetAlarm(..., _RC) - is_ringing = ESMF_AlarmIsRinging(alarm, _RC) - _RETURN_UNLESS(is_ringing) - - call this%update_dependecies() - - -end module mapl3g_ImportCoupler diff --git a/generic3g/couplers/Observable.F90 b/generic3g/couplers/Observable.F90 deleted file mode 100644 index 5f844d568006..000000000000 --- a/generic3g/couplers/Observable.F90 +++ /dev/null @@ -1,84 +0,0 @@ -#include "MAPL_Generic.h" - -module mapl3g_Observable - use mapl_ErrorHandlingMod - implicit none - private - - ! Class - public :: Observable - ! procedures - public :: update_observable - public :: invalidate_observable - - - type, abstract :: Observable - private - logical :: stale = .true. - contains - procedure(I_Notify), deferred :: should_update ! ??? needed? - procedure(I_Notify), deferred :: update_self - procedure(I_Notify), deferred :: invalidate_self - - ! Accessors - procedure, non_overridable :: is_up_to_date - procedure, non_overridable :: is_stale - procedure, non_overridable :: set_up_to_date - procedure, non_overridable :: set_stale - end type Observable - - abstract interface - subroutine I_Notify(this, rc) - import :: Observable - class(Obserer), intent(inout) :: this - integer, optional, intent(out) :: rc - end subroutine I_Notify - end interface - -contains - - subroutine update_observable(this, rc) - class(Observable), intent(inout) :: this - integer, optional, intent(in) :: rc - - _RETURN_IF(this%is_up_to_date()) - - call this%update_self(_RC) - call this%set_up_to_date() - - _RETURN(_SUCCESS) - end subroutine update - - subroutine invalidate(this, rc) - class(Observable), intent(inout) :: this - integer, optional, intent(in) :: rc - - _RETURN_IF(this%is_stale()) - - call this%invalidate_self(_RC) - call this%set_stale() - - _RETURN(_SUCCESS) - end subroutine invalidate - - pure subroutine set_up_to_date(this) - class(Observable), intent(inout) :: this - this%up_to_date = .true - end subroutine set_up_to_date - - pure subroutine set_stale(this) - class(Observable), intent(inout) :: this - this%up_to_date = .false - end subroutine set_stale - - pure logical function is_up_to_date(this) - class(Observable), intent(in) :: this - is_up_to_date = this%up_to_date - end function is_up_to_date - - pure logical function is_stale(this) - class(Observable), intent(in) :: this - is_stale = .not. this%up_to_date - end function is_up_to_date - -end module mapl3g_Observable diff --git a/generic3g/couplers/ObservablePtrVector.F90 b/generic3g/couplers/ObservablePtrVector.F90 deleted file mode 100644 index af47dab70854..000000000000 --- a/generic3g/couplers/ObservablePtrVector.F90 +++ /dev/null @@ -1,14 +0,0 @@ -module mapl3g_ObservablePtrVector - use mapl3g_Observable - -#define T ObservablePtr -#define Vector ObservablePtrVector -#define VectorIterator ObservablePtrVectorIterator - -#include "vector/template.inc" - -#undef T -#undef Vector -#undef VectorIterator - -end module mapl3g_ObservablePtrVector diff --git a/generic3g/couplers/Observed.F90 b/generic3g/couplers/Observed.F90 deleted file mode 100644 index 62e23ebf3f3d..000000000000 --- a/generic3g/couplers/Observed.F90 +++ /dev/null @@ -1,35 +0,0 @@ -#include "MAPL_Generic.h" - -module mapl3g_Observable - use mapl3g_Observer - implicit none - private - - public :: Observable - - type :: Observable - type(ObserverPtrVector) :: observers - contains - procedure :: update_observers - end type Observable - -contains - - subroutine update_observers(this, rc) - class(Observable), intent(inout) :: this - integer, optional, intent(out) :: rc - - integer :: status - - associate (e => this%observers%end()) - iter = this%observers%begin() - do while (iter /= e) - call iter%next() - obsrvr => iter%of() - call obsrvr%update(_RC) - end do - end associate - _RETURN(_SUCCESS) - end subroutine update_observers - -end module mapl3g_Observable diff --git a/generic3g/couplers/Observer.F90 b/generic3g/couplers/Observer.F90 deleted file mode 100644 index 4e69ae57b927..000000000000 --- a/generic3g/couplers/Observer.F90 +++ /dev/null @@ -1,94 +0,0 @@ -#include "MAPL_Generic.h" - -module mapl3g_Observer - use mapl_ErrorHandlingMod - implicit none - private - - ! Class - public :: Observer - public :: ObserverPtr - - ! procedures - public :: update - public :: invalidate - - - type, abstract :: Observer - private - logical :: stale = .true. - contains - procedure(I_Notify), deferred :: should_update ! ??? needed? - procedure(I_Notify), deferred :: update_self - procedure(I_Notify), deferred :: invalidate_self - - ! Accessors - procedure, non_overridable :: is_up_to_date - procedure, non_overridable :: is_stale - procedure, non_overridable :: set_up_to_date - procedure, non_overridable :: set_stale - end type Observer - - type :: ObserverPtr - class(Observer), pointer :: ptr => null() - end type ObserverPtr - - abstract interface - subroutine I_Notify(this, rc) - import :: Observer - class(Observer), intent(inout) :: this - integer, optional, intent(out) :: rc - end subroutine I_Notify - end interface - -contains - - subroutine update(this, rc) - class(Observer), intent(inout) :: this - integer, optional, intent(out) :: rc - - integer :: status - - _RETURN_IF(this%is_up_to_date()) - - call this%update_self(_RC) - call this%set_up_to_date() - - _RETURN(_SUCCESS) - end subroutine update - - subroutine invalidate(this, rc) - class(Observer), intent(inout) :: this - integer, optional, intent(out) :: rc - - integer :: status - - _RETURN_IF(this%is_stale()) - - call this%invalidate_self(_RC) - call this%set_stale() - - _RETURN(_SUCCESS) - end subroutine invalidate - - pure subroutine set_up_to_date(this) - class(Observer), intent(inout) :: this - this%stale = .false. - end subroutine set_up_to_date - - pure subroutine set_stale(this) - class(Observer), intent(inout) :: this - this%stale = .true. - end subroutine set_stale - - pure logical function is_up_to_date(this) - class(Observer), intent(in) :: this - is_up_to_date = .not. this%stale - end function is_up_to_date - - pure logical function is_stale(this) - class(Observer), intent(in) :: this - is_stale = this%stale - end function is_stale - -end module mapl3g_Observer diff --git a/generic3g/couplers/ObserverPtrVector.F90 b/generic3g/couplers/ObserverPtrVector.F90 deleted file mode 100644 index 027cf5640a4e..000000000000 --- a/generic3g/couplers/ObserverPtrVector.F90 +++ /dev/null @@ -1,14 +0,0 @@ -module mapl3g_ObserverPtrVector - use mapl3g_Observer - -#define T ObserverPtr -#define Vector ObserverPtrVector -#define VectorIterator ObserverPtrVectorIterator - -#include "vector/template.inc" - -#undef T -#undef Vector -#undef VectorIterator - -end module mapl3g_ObserverPtrVector diff --git a/generic3g/couplers/outer.F90 b/generic3g/couplers/outer.F90 deleted file mode 100644 index 848f348e81bf..000000000000 --- a/generic3g/couplers/outer.F90 +++ /dev/null @@ -1,96 +0,0 @@ - - - type(ObserverPtrVector) :: export_couplers - type(ObserverPtrVector) :: import_couplers - - ! Connect E --> I - - sequence = cplr(E, I) - call src_comp%add_export_coupler(sequence%first()) - call dst_comp%add_import_coupler(sequence%last()) - - - ! (1) Trivial case: - ! No need to add coupler - ! I and E share field - - ! (2) Regrid - - cplr = Regrid(E, I) - call src_comp%add_export_coupler(cplr) - call dst_comp%add_import_coupler(cplr) - - - ! (3) Change units and then regrid - - cplr1 = ChangeUnits(E, E1) - cplr2 = Regrid(E1, I) - call cplr2%add_import(cplr1) - call cplr1%add_export(cplr2) - - call src_comp%add_export_coupler(cplr1) - call dst_comp%add_import_coupler(cplr2) - - ! dst comp runs - call update_all(dst_comp%import_couplers) - ! triggers - call update(cplr1) ! change units - call update(cplr2) ! regrid - - - ! parent is "this" - coupler = this%registry%connect(C1:E, C2:I) - - export_cplrs = this%get_export_couplers(c1) - import_cplrs => this%get_import_couplers(c2) - - export_cplr => export_cplrs(E) - import_cplr => import_cplrs(I) - - call import_cplr%add_import(export_cplr) ! does not work for complex sequence - call export_cplr%add_import(import_cplr) - - - ! coupler includes import dependencies - - ! always a new cplr for given import - it can only connect once. - ! (except wildcards) - import_cplrs = this%get_import_couplers(C2) ! imports of child C2 - call import_cplrs%push_back(coupler) ! careful not to break internal pointers! - - call i - cplr => this%export_couplers%at(E, _RC) ! extends mapping - if (cplr%size() == 0) then - cplr% - call cplr%add_export(new_couplers%first()) - - ! Child C1 gets the extensions - - - - - couplers is - - - - - subroutine connect(C_e, e, C_i, i) - - coupler_0 => C_e%export_couplers(e) ! possibly null() - - e_0 = e - do while (e_0 /= i) - e_1 => connect_one_step(e_0, i) - coupler_1 => NewCoupler(e_0, e_1) - call coupler_1%add_import(coupler_0) - call coupler_0%add_export(coupler_1) - - e_0 => e_1 - coupler_0 => coupler_1 ! memory leak - end do - - if (.associated(coupler_c)) then - call C_i%import_couplers%push_back(Ptr(last_coupler) - end if - - diff --git a/shared/MAPL_ESMF_InfoKeys.F90 b/shared/MAPL_ESMF_InfoKeys.F90 index c77c2d29a87f..b27657914fd9 100644 --- a/shared/MAPL_ESMF_InfoKeys.F90 +++ b/shared/MAPL_ESMF_InfoKeys.F90 @@ -11,9 +11,11 @@ module mapl3g_esmf_info_keys public :: KEY_UNGRIDDED_DIMS public :: KEY_VERT_DIM public :: KEY_VERT_GRID + public :: KEY_INTERPOLATION_WEIGHTS public :: KEY_UNITS public :: KEY_LONG_NAME public :: KEY_STANDARD_NAME + public :: KEY_TYPEKIND public :: KEY_NUM_LEVELS public :: KEY_VLOC public :: KEY_NUM_UNGRIDDED_DIMS @@ -35,8 +37,10 @@ module mapl3g_esmf_info_keys character(len=*), parameter :: KEY_VERT_DIM = '/vertical_dim' character(len=*), parameter :: KEY_VERT_GRID = '/vertical_grid' character(len=*), parameter :: KEY_UNITS = '/units' + character(len=*), parameter :: KEY_TYPEKIND = '/typekind' character(len=*), parameter :: KEY_LONG_NAME = '/long_name' character(len=*), parameter :: KEY_STANDARD_NAME = '/standard_name' + character(len=*), parameter :: KEY_INTERPOLATION_WEIGHTS = '/interpolation_weights' ! VerticalGeom info keys character(len=*), parameter :: KEY_NUM_LEVELS = KEY_VERT_GRID // '/num_levels'