Test_FieldDelta.pf Source File

#include "MAPL_TestErr.h"
#include "unused_dummy.H"
module Test_FieldDelta
   use mapl3g_FieldDelta
   use mapl3g_FieldGet
   use mapl3g_FieldCreate
   use mapl3g_FieldInfo
   use mapl3g_UngriddedDims
   use mapl3g_UngriddedDim
   use mapl3g_InfoUtilities
   use mapl3g_VerticalStaggerLoc
   use mapl3g_FieldInfo
   use esmf
   use ESMF_TestMethod_mod
   use funit
   implicit none (type, external)

   integer, parameter :: ORIG_VGRID_LEVELS = 5
   real, parameter :: FILL_VALUE = 99.
   character(*), parameter :: ORIGINAL_UNITS = 'm'
   character(*), parameter :: REFERENCE_UNITS = 'km'
   
contains

   @test(type=ESMF_TestMethod, npes=[1])
   subroutine test_change_typekind(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
      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', _RC)
      call MAPL_FieldSet(f, num_levels=ORIG_VGRID_LEVELS+1, vert_staggerloc=VERTICAL_STAGGER_EDGE, _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()))
      @assert_that(typekind == ESMF_TYPEKIND_R8, is(true()))

      call ESMF_FieldDestroy(f, _RC)
      call ESMF_GridDestroy(grid, _RC)
      call ESMF_GeomDestroy(geom, _RC)

      _UNUSED_DUMMY(this)
   end subroutine test_change_typekind

   @test(type=ESMF_TestMethod, npes=[1])
   subroutine test_same_typekind_do_not_reallocate(this)
      class(ESMF_TestMethod), intent(inout) :: this
      type(ESMF_Field) :: f
      type(ESMF_Grid) :: grid
      type(ESMF_Geom) :: geom, other_geom

      integer :: status
      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)
      f = ESMF_FieldCreate(geom, typekind=ESMF_TYPEKIND_R4, name='in', _RC)
      call MAPL_FieldSet(f, num_levels=ORIG_VGRID_LEVELS+1, vert_staggerloc=VERTICAL_STAGGER_EDGE, _RC)

      call ESMF_FieldGet(f, fArrayPtr=x, _RC)
      x = FILL_VALUE

      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()))
      @assert_that(typekind == ESMF_TYPEKIND_R4, is(true()))
      @assert_that(other_geom == geom, is(true()))

      call ESMF_FieldGet(f, fArrayPtr=x, _RC)
      @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_same_typekind_do_not_reallocate

   @test(type=ESMF_TestMethod, npes=[1])
   subroutine test_change_geom(this)
      class(ESMF_TestMethod), intent(inout) :: this
      type(ESMF_Field) :: f
      type(ESMF_Grid), target :: grid1, grid2
      type(ESMF_Geom) :: geom1, geom2

      integer :: status
      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)
      call MAPL_FieldSet(f, num_levels=ORIG_VGRID_LEVELS+1, vert_staggerloc=VERTICAL_STAGGER_EDGE, _RC)

      grid2 = ESMF_GridCreateNoPeriDim(maxIndex=[3,5], name='I_AM_GROOT', _RC)
      geom2 = ESMF_GeomCreate(grid2, _RC)
      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()))

      call ESMF_FieldGet(f, fArrayPtr=x, _RC)
      @assert_that(shape(x),is(equal_to([3,5])))

      call ESMF_FieldDestroy(f, _RC)
      call ESMF_GridDestroy(grid1, _RC)
      call ESMF_GridDestroy(grid2, _RC)
      call ESMF_GeomDestroy(geom2, _RC)

       _UNUSED_DUMMY(this)
  end subroutine test_change_geom

   @test(type=ESMF_TestMethod, npes=[1])
   subroutine test_same_geom_do_not_reallocate(this)
      class(ESMF_TestMethod), intent(inout) :: this
      type(ESMF_Field) :: f
      type(ESMF_Grid), target :: grid1
      type(ESMF_Geom) :: geom1
      type(ESMF_Geom) :: geom2

      integer :: status
      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)
      call MAPL_FieldSet(f, num_levels=ORIG_VGRID_LEVELS+1, vert_staggerloc=VERTICAL_STAGGER_EDGE, _RC)
      call ESMF_FieldGet(f, fArrayPtr=x, _RC)
      x = FILL_VALUE

      geom2 = geom1
      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()))

      call ESMF_FieldGet(f, fArrayPtr=x, _RC)
      @assert_that(all(x == FILL_VALUE), is(true()))

      call ESMF_FieldDestroy(f, _RC)
      call ESMF_GridDestroy(grid1, _RC)
      call ESMF_GeomDestroy(geom2, _RC)

      _UNUSED_DUMMY(this)
   end subroutine test_same_geom_do_not_reallocate


    @test(type=ESMF_TestMethod, npes=[1])
   ! Probably exceedingly rare, but MAPL3 allows the vertical grid to change with time
   ! which could change the number of levels ...
   subroutine test_change_n_levels(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
      integer, parameter :: NEW_NUM_LEVELS = 7

      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=[ORIG_VGRID_LEVELS+1,3], _RC)
      call MAPL_FieldSet(f, num_levels=ORIG_VGRID_LEVELS+1, vert_staggerloc=VERTICAL_STAGGER_EDGE, _RC)

      delta = FieldDelta(num_levels=NEW_NUM_LEVELS+1) ! edge
      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,NEW_NUM_LEVELS+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_FieldSet(f, num_levels=0, vert_staggerloc=VERTICAL_STAGGER_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,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_surface
   


    @test(type=ESMF_TestMethod, npes=[1])
   subroutine test_same_n_levels_do_not_reallocate(this)
      class(ESMF_TestMethod), intent(inout) :: this
      type(ESMF_Field) :: f
      type(ESMF_Grid) :: grid
      type(ESMF_Geom) :: geom, other_geom

      integer :: status
      type(ESMF_FieldStatus_Flag) :: field_status
      real(ESMF_KIND_R4), pointer :: x(:,:,:,:)
      type(ESMF_TypeKind_Flag) :: typekind
      type(FieldDelta) :: delta
      type(UngriddedDims) :: ungridded_dims

      grid = ESMF_GridCreateNoPeriDim(maxIndex=[4,4], name='I_AM_GROOT', _RC)
      geom = ESMF_GeomCreate(grid, _RC)

      call ungridded_dims%add_dim(UngriddedDim(3))
      f = MAPL_FieldCreate(geom=geom, typekind=ESMF_TYPEKIND_R4, &
          num_levels=ORIG_VGRID_LEVELS,  ungridded_dims=ungridded_dims, &
          vert_staggerloc=VERTICAL_STAGGER_CENTER, _RC)
      call ESMF_FieldGet(f, fArrayPtr=x, _RC)
      x = FILL_VALUE

      delta = FieldDelta(num_levels=ORIG_VGRID_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()))
      @assert_that(typekind == ESMF_TYPEKIND_R4, is(true()))
      @assert_that(other_geom == geom, is(true()))

      ! Check that Field data is unchanged
      call ESMF_FieldGet(f, fArrayPtr=x, _RC)
      @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_same_n_levels_do_not_reallocate


    @test(type=ESMF_TestMethod, npes=[1])
   subroutine test_field_update_from_field_ignore_geom(this)
      class(ESMF_TestMethod), intent(inout) :: this
      type(ESMF_Field) :: f, f_ref
      type(ESMF_Grid) :: grid, grid_ref
      type(ESMF_Geom) :: geom, geom_ref, new_geom
      character(:), allocatable :: new_units

      integer :: status
      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)

      grid_ref = ESMF_GridCreateNoPeriDim(maxIndex=[7,7], name='I_AM_GROOT', _RC)
      geom_ref = ESMF_GeomCreate(grid_ref, _RC)

      f = ESMF_FieldCreate(geom, typekind=ESMF_TYPEKIND_R4, name='in', &
           ungriddedLbound=[1,1], ungriddedUbound=[ORIG_VGRID_LEVELS,3], _RC)

      f_ref = ESMF_FieldCreate(geom_ref, typekind=ESMF_TYPEKIND_R8, name='in', &
           ungriddedLbound=[1,1], ungriddedUbound=[ORIG_VGRID_LEVELS-1,3], _RC)

      call MAPL_FieldSet(f, num_levels=ORIG_VGRID_LEVELS, vert_staggerloc=VERTICAL_STAGGER_CENTER, &
           units=ORIGINAL_UNITS, _RC)
      call MAPL_FieldSet(f_ref, num_levels=ORIG_VGRID_LEVELS, vert_staggerloc=VERTICAL_STAGGER_CENTER, &
           units=REFERENCE_UNITS, _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()))
      @assert_that(typekind == ESMF_TYPEKIND_R8, is(true()))
      @assert_that(new_geom == geom, is(true()))

      call MAPL_FieldGet(f, units=new_units, _RC)
      @assertEqual(REFERENCE_UNITS, new_units)

      ! 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,ORIG_VGRID_LEVELS,3])))

      call ESMF_FieldDestroy(f, _RC)
      call ESMF_GridDestroy(grid, _RC)
      call ESMF_GeomDestroy(geom, _RC)

      _UNUSED_DUMMY(this)
   end subroutine test_field_update_from_field_ignore_geom

end module Test_FieldDelta