Test_FieldBundleDelta.pf

Source Code

#include "MAPL_TestErr.h"
#include "unused_dummy.H"
module Test_FieldBundleDelta
   use mapl3g_FieldBundleDelta
   use mapl3g_FieldBundleGet
   use mapl3g_FieldDelta
   use mapl3g_FieldGet
   use mapl3g_FieldCreate
   use mapl3g_FieldInfo
   use mapl3g_esmf_info_keys
   use mapl3g_VerticalStaggerLoc
   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_VGRID = 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(LU_Bound), allocatable :: bounds(:)
      type(VerticalStaggerLoc) :: vert_staggerloc
      integer, allocatable :: num_levels

      ungridded_dims = UngriddedDims()
      bounds = ungridded_dims%get_bounds()

      vert_staggerloc = VERTICAL_STAGGER_NONE
      if (present(with_ungridded)) then
         if (with_ungridded) then
            vert_staggerloc = VERTICAL_STAGGER_CENTER
            num_levels = NUM_LEVELS_VGRID
            call ungridded_dims%add_dim(UngriddedDim(NUM_RADII, "radius", 'nm'))
         end if
      end if
      field = ESMF_FieldEmptyCreate()
      call ESMF_FieldEmptySet(field, geom=geom)
      call MAPL_FieldEmptyComplete(field, typekind=typekind, &
           num_levels=num_levels, vert_staggerloc=vert_staggerloc, &
           ungridded_dims=ungridded_dims, &
           units=units)
      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(VerticalStaggerLoc) :: vert_staggerloc

      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

      call MAPL_FieldBundleSet(bundle, interpolation_weights=weights, typekind=typekind, units=units)

      vert_staggerloc = VERTICAL_STAGGER_NONE
      ungridded_dims = UngriddedDims()
      if (present(with_ungridded)) then
         if (with_ungridded) then
            vert_staggerloc = VERTICAL_STAGGER_CENTER
            call MAPL_FieldBundleSet(bundle, num_levels=NUM_LEVELS_VGRID)
            call ungridded_dims%add_dim(UngriddedDim(NUM_RADII, "radius", 'nm'))
         end if
      end if
      call MAPL_FieldBundleSet(bundle, vert_staggerloc=vert_staggerloc)

      call MAPL_FieldBundleSet(bundle, ungridded_dims=ungridded_dims)

   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_FieldBundleGet(bundle, interpolation_weights=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_FieldGet(fieldList(i), units=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_FieldGet(fieldList(i), units=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_FieldGet(fieldList(i), units=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_FieldGet(fieldList(i), units=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_FieldBundleGet(bundle, interpolation_weights=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 :: nlevels, rank
      type(UngriddedDims) :: ungridded_dims

      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_VGRID,NUM_RADII])))
         @assert_that(all(x_r4 == FILL_VALUE), is(true()))

         call MAPL_FieldGet(fieldList(i), units=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_FieldGet(fieldList(i), ungridded_dims=ungridded_dims, _RC)
         @assert_that(ungridded_dims%get_num_ungridded(), is(1))

         call MAPL_FieldGet(fieldList(i), num_levels=nlevels, _RC)
         @assert_that(nlevels, is(NUM_LEVELS_VGRID))

      end do

      call MAPL_FieldBundleGet(bundle, interpolation_weights=weights, _RC)
      @assert_that(weights, is(equal_to(new_weights)))

      call MAPL_FieldBundleGet(bundle, ungridded_dims=ungridded_dims, _RC)
      @assert_that(ungridded_dims%get_num_ungridded(), is(1))

      call MAPL_FieldBundleGet(bundle, num_levels=nlevels, _RC)
      @assert_that(nlevels, is(NUM_LEVELS_VGRID))

      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_FieldGet(fieldList(i), units=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_FieldBundleGet(bundle, interpolation_weights=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 :: nlevels
      type(UngriddedDims) :: ungridded_dims

      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_VGRID,NUM_RADII])))

         call MAPL_FieldGet(fieldList(i), units=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_FieldGet(fieldList(i), ungridded_dims=ungridded_dims, _RC)
         @assert_that(ungridded_dims%get_num_ungridded(), is(1))

         call MAPL_FieldGet(fieldList(i), num_levels=nlevels, _RC)
         @assert_that(nlevels, is(NUM_LEVELS_VGRID))
      end do

      call MAPL_FieldBundleGet(bundle, interpolation_weights=weights, _RC)
      @assert_that(weights, is(equal_to(new_weights)))

      call MAPL_FieldBundleGet(bundle, ungridded_dims=ungridded_dims, _RC)
      @assert_that(ungridded_dims%get_num_ungridded(), is(1))

      call MAPL_FieldBundleGet(bundle, num_levels=nlevels, _RC)
      @assert_that(nlevels, is(NUM_LEVELS_VGRID))

      call teardown_bundle(bundle)
      call teardown_geom(geom)

      _UNUSED_DUMMY(this)
   end subroutine test_create_fields_with_ungridded

end module Test_FieldBundleDelta
Test_FieldBundleDelta.pf Source File

Contents

Modules

Source Code

This file depends on

Source Code
