search

esmf-org / esmf

The Earth System Modeling Framework (ESMF) is a suite of software tools for developing high-performance, multi-component Earth science modeling applications.
https://earthsystemmodeling.org/
Other
156 stars 75 forks source link

Need improved interface(s) for contructing field "slices" #189

Open tclune opened 11 months ago

tclune commented 11 months ago

Because dynamic masking does not work with fields that have ungridded dimensions after the geom dimensions, I need to loop over the ungridded dimensions and create temporary fields that are associated with slices of the original data.

This presents several challenges with the existing interfaces:

  1. Needs to be typekind and rank agnostic - results in unfortunate nested block of code and declaration of large numbers of variant local array pointers.
  2. Cannot directly copy the geom.
    • must first query the field for the type of geom
    • must then get the geom subtype
    • must then wrap the subtype in a new geom object
  3. There is apparently no mechanism to detect the staggering of the original field, so the staggering cannot be maintained in this process.

Below is a bit of client code that addresses (2), but skips (1). I suspect that a bit more infrastructure on the ESMF side will make this much cleaner.

     function get_slice(f, k, rc) result(f_slice)
         type(ESMF_Field) :: f_slice
         type(ESMF_Field), intent(inout) :: f
         integer, intent(in) :: k
         integer, optional, intent(out) :: rc

         integer :: status
         real(kind=ESMF_KIND_R4), pointer :: x(:,:,:)
         real(kind=ESMF_KIND_R4), pointer :: x_slice(:,:)
         type(ESMF_Geom) :: geom
         type(ESMF_GeomType_Flag) :: geomtype
         type(ESMF_Grid) :: grid
         type(ESMF_Mesh) :: mesh
         type(ESMF_XGrid) :: xgrid
         type(ESMF_LocStream) :: locstream

         call ESMF_FieldGet(f, farrayptr=x, _RC)
         call ESMF_FieldGet(f, geomtype=geomtype, _RC)

         if (geomtype == ESMF_GEOMTYPE_GRID) then
            call ESMF_FieldGet(f, grid=grid, _RC)
            geom = ESMF_GeomCreate(grid, _RC)
         elseif (geomtype == ESMF_GEOMTYPE_MESH) then
            call ESMF_FieldGet(f, mesh=mesh, _RC)
            geom = ESMF_GeomCreate(mesh, _RC)
         elseif (geomtype == ESMF_GEOMTYPE_XGRID) then
            call ESMF_FieldGet(f, xgrid=xgrid, _RC)
            geom = ESMF_GeomCreate(xgrid, _RC)
         elseif (geomtype == ESMF_GEOMTYPE_LOCSTREAM) then
            call ESMF_FieldGet(f, locstream=locstream, _RC)
            geom = ESMF_GeomCreate(locstream, _RC)
         else
            _FAIL('Invalid geometry type.')
         end if

         x_slice => x(:,:,k)
         f_slice = ESMF_FieldCreate(geom, &
              datacopyflag=ESMF_DATACOPY_REFERENCE, &
              farrayptr=x_slice, _RC)

         call ESMF_GeomDestroy(geom, _RC)

         _RETURN(_SUCCESS)
      end function get_slice