Merge pull request #3154 from GEOS-ESM/feature/pchakrab/vertical-regr…

…idding

Feature/pchakrab/vertical regridding

generic3g/actions/VerticalRegridAction.F90

@@ -23,7 +23,7 @@ module mapl3g_VerticalRegridAction
 
    type, extends(ExtensionAction) :: VerticalRegridAction
       type(ESMF_Field) :: v_in_coord, v_out_coord
-      type(SparseMatrix_sp), allocatable :: matrix(:, :)
+      type(SparseMatrix_sp), allocatable :: matrix(:)
       type(GriddedComponentDriver), pointer :: v_in_coupler => null()
       type(GriddedComponentDriver), pointer :: v_out_coupler => null()
       type(VerticalRegridMethod) :: method = VERTICAL_REGRID_UNKNOWN
@@ -67,7 +67,7 @@ subroutine initialize(this, importState, exportState, clock, rc)
 
       real(ESMF_KIND_R4), pointer :: v_in(:, :, :), v_out(:, :, :)
       integer :: shape_in(3), shape_out(3), n_horz, n_ungridded
-      integer :: horz1, horz2, ungrd, status
+      integer :: horz, ungrd, status
 
       _ASSERT(this%method == VERTICAL_REGRID_LINEAR, "regrid method can only be linear")
 
@@ -89,16 +89,14 @@ subroutine initialize(this, importState, exportState, clock, rc)
       _ASSERT((shape_in(1) == shape_out(1)), "horz dims are expected to be equal")
       _ASSERT((shape_in(3) == shape_out(3)), "ungridded dims are expected to be equal")
 
-      allocate(this%matrix(n_horz, n_horz))
+      allocate(this%matrix(n_horz))
 
       ! TODO: Convert to a `do concurrent` loop
-      do horz1 = 1, n_horz
-         do horz2 = 1, n_horz
-            do ungrd = 1, n_ungridded
-               associate(src => v_in(horz1, :, ungrd), dst => v_out(horz2, :, ungrd))
-                 call compute_linear_map(src, dst, this%matrix(horz1, horz2), _RC)
-               end associate
-            end do
+      do horz = 1, n_horz
+         do ungrd = 1, n_ungridded
+            associate(src => v_in(horz, :, ungrd), dst => v_out(horz, :, ungrd))
+              call compute_linear_map(src, dst, this%matrix(horz), _RC)
+            end associate
          end do
       end do
 
@@ -117,7 +115,7 @@ subroutine update(this, importState, exportState, clock, rc)
       type(ESMF_Field) :: f_in, f_out
       real(ESMF_KIND_R4), pointer :: x_in(:,:,:), x_out(:,:,:)
       integer :: shape_in(3), shape_out(3), n_horz, n_ungridded
-      integer :: horz1, horz2, ungrd
+      integer :: horz, ungrd
 
       ! if (associated(this%v_in_coupler)) then
       !    call this%v_in_coupler%run(phase_idx=GENERIC_COUPLER_UPDATE, _RC)
@@ -140,8 +138,8 @@ subroutine update(this, importState, exportState, clock, rc)
       _ASSERT((shape_in(1) == shape_out(1)), "horz dims are expected to be equal")
       _ASSERT((shape_in(3) == shape_out(3)), "ungridded dims are expected to be equal")
 
-      do concurrent (horz1=1:n_horz, horz2=1:n_horz, ungrd=1:n_ungridded)
-         x_out(horz2, :, ungrd) = matmul(this%matrix(horz1, horz2), x_in(horz1, :, ungrd))
+      do concurrent (horz=1:n_horz, ungrd=1:n_ungridded)
+         x_out(horz, :, ungrd) = matmul(this%matrix(horz), x_in(horz, :, ungrd))
       end do
 
       _RETURN(_SUCCESS)

generic3g/tests/scenarios/vertical_regridding/A.yaml

@@ -16,8 +16,8 @@ mapl:
   states:
     import: {}
     export:
-      E:
-        standard_name: 'E'
-        units: 'm'
-        default_value: 1.
-        vertical_dim_spec: center # or edge
+      E_A:
+        standard_name: E_A
+        units: m
+        default_value: 15.
+        vertical_dim_spec: center

generic3g/tests/scenarios/vertical_regridding/B.yaml

@@ -15,9 +15,8 @@ mapl:
 
   states:
     import:
-      I:
-        standard_name: 'I'
-        units: 'm'
-        default_value: 1.
+      I_B:
+        standard_name: I_B
+        units: m
         vertical_dim_spec: center
     export: {}

generic3g/tests/scenarios/vertical_regridding/expectations.yaml

@@ -5,8 +5,8 @@
 
 - component: A
   export:
-    E: {status: complete}
+    E_A: {status: complete}
 
 - component: B
   import:
-    I: {status: complete}
+    I_B: {status: complete}

generic3g/tests/scenarios/vertical_regridding/parent.yaml

@@ -6,13 +6,14 @@ mapl:
       setServices: setservices_
       config_file: scenarios/vertical_regridding/A.yaml
     B:
-      dso: libsimple_leaf_gridcomp
+      sharedObj: libsimple_leaf_gridcomp
+      setServices: setservices_
       config_file: scenarios/vertical_regridding/B.yaml
 
   states: {}
 
   connections:
-    - src_name: E
-      dst_name: I
+    - src_name: E_A
+      dst_name: I_B
       src_comp: A
       dst_comp: B

generic3g/tests/scenarios/vertical_regridding_2/A.yaml

@@ -9,15 +9,15 @@ mapl:
       dateline: DC
     vertical_grid:
       class: model
-      short_name: "PLE"
-      units: "hPa"
+      short_name: PLE
+      units: hPa
       num_levels: 4
 
   states:
     import: {}
     export:
       PLE:
-        standard_name: "E"
-        units: "hPa"
+        standard_name: air_pressure
+        units: hPa
         default_value: 17.
-        vertical_dim_spec: "edge"
+        vertical_dim_spec: center

generic3g/tests/scenarios/vertical_regridding_2/B.yaml

@@ -15,9 +15,8 @@ mapl:
 
   states:
     import:
-      I:
-        standard_name: "I"
-        units: "hPa"
-        default_value: 1.
-        vertical_dim_spec: edge
+      I_B:
+        standard_name: I_B
+        units: hPa
+        vertical_dim_spec: center
     export: {}

generic3g/tests/scenarios/vertical_regridding_2/C.yaml

@@ -3,8 +3,8 @@ mapl:
   geometry:
     esmf_geom:
       class: latlon
-      im_world: 2
-      jm_world: 3
+      im_world: 12
+      jm_world: 13
       pole: PC
       dateline: DC
     vertical_grid:
@@ -17,7 +17,7 @@ mapl:
     import: {}
     export:
       ZLE:
-        standard_name: E
+        standard_name: height
         units: m
-        default_value: 17.
-        vertical_dim_spec: edge
+        default_value: 23.
+        vertical_dim_spec: center

generic3g/tests/scenarios/vertical_regridding_2/D.yaml

@@ -3,21 +3,20 @@ mapl:
   geometry:
     esmf_geom:
       class: latlon
-      im_world: 2
-      jm_world: 3
+      im_world: 12
+      jm_world: 13
       pole: PC
       dateline: DC
     vertical_grid:
       class: fixed_levels
       standard_name: height
       units: m
-      levels: [17.]
+      levels: [23.]
 
   states:
     import:
-      I:
-        standard_name: I
+      I_D:
+        standard_name: I_D
         units: m
-        default_value: 1.
-        vertical_dim_spec: edge
+        vertical_dim_spec: center
     export: {}

generic3g/tests/scenarios/vertical_regridding_2/expectations.yaml

@@ -9,12 +9,12 @@
 
 - component: B
   import:
-    I: {status: complete}
+    I_B: {status: complete}
 
 - component: C
   export:
     ZLE: {status: complete}
 
 - component: D
   import:
-    I: {status: complete}
+    I_D: {status: complete}

generic3g/tests/scenarios/vertical_regridding_2/parent.yaml

@@ -22,10 +22,10 @@ mapl:
 
   connections:
     - src_name: PLE
-      dst_name: I
+      dst_name: I_B
       src_comp: A
       dst_comp: B
     - src_name: ZLE
-      dst_name: I
+      dst_name: I_D
       src_comp: C
       dst_comp: D

generic3g/vertical/FixedLevelsVerticalGrid.F90

@@ -8,8 +8,6 @@ module mapl3g_FixedLevelsVerticalGrid
    use mapl3g_FieldCreate
    use mapl3g_GriddedComponentDriver
    use mapl3g_VerticalDimSpec
-   use mapl3g_InfoUtilities, only: MAPL_InfoSetInternal
-   use mapl3g_esmf_info_keys, only: KEY_VLOC, KEY_NUM_LEVELS
    use esmf
 
    use, intrinsic :: iso_fortran_env, only: REAL32
@@ -75,30 +73,30 @@ subroutine get_coordinate_field(this, field, coupler, standard_name, geom, typek
       type(VerticalDimSpec), intent(in) :: vertical_dim_spec
       integer, optional, intent(out) :: rc
 
-      real(kind=REAL32), allocatable :: adjusted_levels(:)
-      character(:), allocatable :: vloc
-      integer :: status
-
-      ! KLUDGE - for VERTICAL_DIM_EDGE, we simply extend the the size
-      ! [40, 30, 20, 10] -> [40, 30, 20, 10, 10]
-      ! Also, vloc assignment gets simpler once we have co-located description in VerticalDimSpec
-      if (vertical_dim_spec == VERTICAL_DIM_CENTER) then
-         adjusted_levels = this%levels
-         vloc = "VERTICAL_DIM_CENTER"
-      else if (vertical_dim_spec == VERTICAL_DIM_EDGE) then
-         adjusted_levels = [this%levels, this%levels(size(this%levels))]
-         vloc = "VERTICAL_DIM_EDGE"
-      else
-         _FAIL("invalid vertical_dim_spec")
-      end if
+      real(kind=REAL32), pointer :: farray3d(:, :, :)
+      integer, allocatable :: local_cell_count(:)
+      integer :: i, j, IM, JM, status
 
-      field = esmf_field_create_(geom, adjusted_levels, _RC)
+      field = MAPL_FieldCreate( &
+           geom=geom, &
+           typekind=ESMF_TYPEKIND_R4, &
+           num_levels=size(this%levels), &
+           vert_staggerloc=VERTICAL_STAGGER_CENTER, &
+           _RC)
+      ! Copy the 1D array, levels(:), to each point on the horz grid
+      call ESMF_FieldGet(field, fArrayPtr=farray3d, _RC)
+      call MAPL_GeomGet_(geom, localCellCount=local_cell_count, _RC)
+      IM = local_cell_count(1); JM = local_cell_count(2)
+      do concurrent (i=1:IM, j=1:JM)
+         farray3d(i, j, :) = this%levels(:)
+      end do
 
       _RETURN(_SUCCESS)
       _UNUSED_DUMMY(coupler)
       _UNUSED_DUMMY(standard_name)
       _UNUSED_DUMMY(typekind)
       _UNUSED_DUMMY(units)
+      _UNUSED_DUMMY(vertical_dim_spec)
    end subroutine get_coordinate_field
 
    logical function can_connect_to(this, src, rc)
@@ -149,55 +147,6 @@ impure elemental logical function not_equal_FixedLevelsVerticalGrid(a, b) result
       not_equal = .not. (a==b)
    end function not_equal_FixedLevelsVerticalGrid
 
-   ! Create an ESMF_Field containing a 3D array that is replicated from
-   ! a 1D array at each point of the horizontal grid
-   function esmf_field_create_(geom, farray1d, rc) result(field)
-      type(ESMF_Field) :: field ! result
-      type(ESMF_Geom), intent(in) :: geom
-      real(kind=REAL32), intent(in) :: farray1d(:)
-!#      character(len=*), intent(in) :: vloc
-      integer, optional, intent(out) :: rc
-
-      integer, allocatable :: local_cell_count(:)
-      real(kind=REAL32), pointer :: farray3d(:, :, :)
-      integer :: i, j, IM, JM, status
-
-!#      ! First, copy the 1D array, farray1d, to each point on the horz grid
-!#      allocate(farray3d(IM, JM, size(farray1d)))
-!#      do concurrent (i=1:IM, j=1:JM)
-!#         farray3d(i, j, :) = farray1d(:)
-!#      end do
-
-      ! Create an ESMF_Field containing farray3d
-      field = MAPL_FieldCreate( &
-           geom=geom, typekind=ESMF_TYPEKIND_R4, &
-           num_levels=size(farray1d), &
-           vert_staggerloc=VERTICAL_STAGGER_CENTER, &
-           _RC)
-
-!#      ! First, copy the 1D array, farray1d, to each point on the horz grid
-      call ESMF_FieldGet(field, fArrayPtr=farray3d, _RC)
-      call MAPL_GeomGet_(geom, localCellCount=local_cell_count, _RC)
-      IM = local_cell_count(1); JM = local_cell_count(2)
-      do concurrent (i=1:IM, j=1:JM)
-         farray3d(i, j, :) = farray1d(:)
-      end do
-
-!#      field = ESMF_FieldCreate( &
-!#           geom=geom, &
-!#           farray=farray3d, &
-!#           indexflag=ESMF_INDEX_DELOCAL, &
-!#           datacopyFlag=ESMF_DATACOPY_VALUE, &
-!#           ungriddedLBound=[1], &
-!#           ungriddedUBound=[size(farray1d)], &
-!#           _RC)
-!#      
-!#      call MAPL_InfoSetInternal(field, key=KEY_NUM_LEVELS, value=size(farray1d), _RC)
-!#      call MAPL_InfoSetInternal(field, key=KEY_VEVLOC, value=vloc, _RC)
-
-      _RETURN(_SUCCESS)
-   end function esmf_field_create_
-
    ! Temporary version here while the detailed MAPL_GeomGet utility gets developed
    subroutine MAPL_GeomGet_(geom, localCellCount, rc)
       use MAPLBase_Mod