Merge pull request #634 from billsacks/handle_rank3_fields

For water tracers: Add handling of rank-3 fields and refactor med_phases_prep_glc_map_lnd2glc

### Description of changes

Two sets of changes that will be needed for water tracers - and in particular, for the tracer version of Flgl_qice_elev, which will be the first field in CMEPS that has two ungridded dimensions:

(1) Extend various methods to support rank-3 fields

(2) Refactor med_phases_prep_glc_map_lnd2glc to pre-compute vertical interpolation weights; this will be helpful when we introduce Flgl_qice_elev_wtracers to avoid duplicating this calculation of vertical interpolation weights between the non-tracer (rank-2 on the lnd grid and rank-1 on the glc grid) and tracer (rank-3 on the lnd grid and rank-2 on the glc grid) fields.

### Specific notes

Contributors other than yourself, if any: Claude Code did a lot of the writing of this code, but I have reviewed it all carefully

CMEPS Issues Fixed (include github issue #):

Are changes expected to change answers? (specify if bfb, different at roundoff, more substantial) bfb when using #633 as a baseline

Any User Interface Changes (namelist or namelist defaults changes)? No

### Testing performed

In the context of cesm3_0_alpha08e, ran:

(1) aux_glc on derecho plus a few extra tests, with comparisons against baselines generated with #633 .

The extra tests beyond aux_glc were:

```
SMS_Ly2.f09_g17_gris20.T1850Gg.derecho_intel
SMS_Ld5.f10_f10_ais8gris4_mg37.I1850Clm50SpGag.derecho_intel.cism-test_coupling
SMS_Lm13.f10_f10_mg37.I1850Clm50SpG.derecho_intel
ERS_Ld5.ne30pg3_t232.B1850C_LTso.derecho_intel.allactive-decstart
ERS_Ld5.ne30pg3_t232.BHISTC_LTso.derecho_intel.allactive-decstart
```

All tests passed and were bit-for-bit, except for this failure that also failed in the baseline: `FAIL NCK_Ly3.f09_g17_gris20.T1850Gg.derecho_gnu COMPARE_base_multiinst`

(2) Full prealpha testing on derecho and izumi with comparison against cesm3_0_alpha08e

Tests passed except for tests that also failed in the baselines or seemed to have passed due to tweaks made in the baseline (e.g., increasing wallclock time for `SUB_D_Ln9.ne3pg3_ne3pg3_mt232.FHIST.izumi_nag.cam-outfrq9s`).

Baseline failures were all as expected (except that a few B compset tests failed baseline comparison due to a difference in which cpl.hx.ww3 files were present in the run vs. baseline, which presumably is due to something outside of this PR):

```
ERS_Ld5.ne30pg3_t232.B1850C_LTso.derecho_intel.allactive-decstart
ERS_Ld5.ne30pg3_t232.BHISTC_LTso.derecho_intel.allactive-decstart
ERS_Ly7.f09_g17_gris4.T1850Gg.derecho_intel
MULTINOAIS_Ly2.f10_f10_ais8gris4_mg37.I1850Clm50SpRsGag.derecho_intel.cism-change_params
SMS_Lm13.f10_f10_mg37.I1850Clm50SpG.derecho_intel
SMS_D_Ly1.f09_g17_ais8.T1850Ga.derecho_gnu
SMS_D_Ld5_P24x1.f10_f10_ais8gris4_mg37.I1850Clm50SpGag.izumi_nag.cism-test_coupling
```

Note that I showed that these differences go away when comparing against baselines generated with #633 , for similar or identical tests.

So I think it's safe to conclude that this is bit-for-bit with #633 .

Author: billsacks

mediator/med_io_mod.F90

@@ -729,7 +729,7 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
     type(ESMF_CoordSys_Flag)      :: coordsys
     integer                       :: rcode
     integer                       :: nf,ns,ng
-    integer                       :: k,n
+    integer                       :: k,n,n2
     integer                       :: ndims, nelements
     integer    ,target            :: dimid2(2)
     integer    ,target            :: dimid3(3)
@@ -754,15 +754,22 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
     integer, pointer              :: Dof(:)
     real(r8), pointer             :: fldptr1(:)
     real(r8), pointer             :: fldptr2(:,:)
+    real(r8), pointer             :: fldptr3(:,:,:)
     real(r8), allocatable         :: ownedElemCoords(:), ownedElemCoords_x(:), ownedElemCoords_y(:)
     character(CS)                 :: cnumber
+    character(CS)                 :: cnumber2
     character(CL)                 :: tmpstr
     type(ESMF_Field)              :: lfield
     integer                       :: rank
-    integer                       :: ungriddedUBound(1) ! currently the size must equal 1 for rank 2 fields
-    integer                       :: gridToFieldMap(1)  ! currently the size must equal 1 for rank 2 fields
     logical                       :: tiles
     character(CL), allocatable    :: fieldNameList(:)
+
+    ! For a single ungridded dimension, there will be 1 element in ungriddedUBound and 1
+    ! element in gridToFieldMap; for two ungridded dimensions, there will be 2 elements in
+    ! ungriddedUBound but still 1 element in gridToFieldMap.
+    integer                       :: ungriddedUBound(2)
+    integer                       :: gridToFieldMap(1)
+
     character(*),parameter :: subName = '(med_io_write_FB) '
     !-------------------------------------------------------------------------------
 
@@ -935,12 +942,47 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
 
           ! TODO (mvertens, 2019-03-13): this is a temporary mod to NOT write hgt
           if (trim(itemc) /= "hgt") then
-             if (rank == 2) then
+             if (rank == 3) then
                 call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound, rc=rc)
                 if (chkerr(rc,__LINE__,u_FILE_u)) return
                 write(cnumber,'(i0)') ungriddedUbound(1)
+                write(cnumber2,'(i0)') ungriddedUbound(2)
                 call ESMF_LogWrite(trim(subname)//':'//'field '//trim(itemc)// &
-                     ' has an griddedUBound of  '//trim(cnumber), ESMF_LOGMSG_INFO)
+                     ' has an ungriddedUBound of  '//trim(cnumber)//' x '//trim(cnumber2), ESMF_LOGMSG_INFO)
+
+                ! Create a new output variable for each element of the 2 ungridded dimensions
+                do n2 = 1,ungriddedUBound(2)
+                   do n = 1,ungriddedUBound(1)
+                      write(cnumber,'(i0)') n
+                      write(cnumber2,'(i0)') n2
+                      name1 = trim(lpre)//'_'//trim(itemc)//trim(cnumber)//'_'//trim(cnumber2)
+                      call ESMF_LogWrite(trim(subname)//': defining '//trim(name1), ESMF_LOGMSG_INFO)
+                      if (luse_float) then
+                         rcode = pio_def_var(io_file, trim(name1), PIO_REAL, dimid, varid)
+                         rcode = pio_put_att(io_file, varid,"_FillValue",real(lfillvalue,r4))
+                      else
+                         rcode = pio_def_var(io_file, trim(name1), PIO_DOUBLE, dimid, varid)
+                         rcode = pio_put_att(io_file,varid,"_FillValue",lfillvalue)
+                      end if
+                      if (NUOPC_FieldDictionaryHasEntry(trim(itemc))) then
+                         call NUOPC_FieldDictionaryGetEntry(itemc, canonicalUnits=cunit, rc=rc)
+                         if (chkerr(rc,__LINE__,u_FILE_u)) return
+                         rcode = pio_put_att(io_file, varid, "units", trim(cunit))
+                      end if
+                      rcode = pio_put_att(io_file, varid, "standard_name", trim(name1))
+                      if (present(tavg)) then
+                         if (tavg) then
+                            rcode = pio_put_att(io_file, varid, "cell_methods", "time: mean")
+                         endif
+                      endif
+                   end do
+                end do
+             else if (rank == 2) then
+                call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound, rc=rc)
+                if (chkerr(rc,__LINE__,u_FILE_u)) return
+                write(cnumber,'(i0)') ungriddedUbound(1)
+                call ESMF_LogWrite(trim(subname)//':'//'field '//trim(itemc)// &
+                     ' has an ungriddedUBound of  '//trim(cnumber), ESMF_LOGMSG_INFO)
 
                 ! Create a new output variable for each element of the undistributed dimension
                 do n = 1,ungriddedUBound(1)
@@ -958,7 +1000,7 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
                       if (NUOPC_FieldDictionaryHasEntry(trim(itemc))) then
                          call NUOPC_FieldDictionaryGetEntry(itemc, canonicalUnits=cunit, rc=rc)
                          if (chkerr(rc,__LINE__,u_FILE_u)) return
-                         rcode = pio_put_att(io_file, varid, "units"        , trim(cunit))
+                         rcode = pio_put_att(io_file, varid, "units", trim(cunit))
                       end if
                       rcode = pio_put_att(io_file, varid, "standard_name", trim(name1))
                       if (present(tavg)) then
@@ -968,15 +1010,15 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
                       endif
                    end if
                 end do
-             else
+             else if (rank == 1) then
                 name1 = trim(lpre)//'_'//trim(itemc)
-                call ESMF_LogWrite(trim(subname)//':'//trim(itemc)//':'//trim(name1),ESMF_LOGMSG_INFO)
+                call ESMF_LogWrite(trim(subname)//': defining '//trim(name1), ESMF_LOGMSG_INFO)
                 if (luse_float) then
                    rcode = pio_def_var(io_file, trim(name1), PIO_REAL, dimid, varid)
-                   rcode = pio_put_att(io_file, varid, "_FillValue", real(lfillvalue, r4))
+                   rcode = pio_put_att(io_file, varid,"_FillValue",real(lfillvalue,r4))
                 else
                    rcode = pio_def_var(io_file, trim(name1), PIO_DOUBLE, dimid, varid)
-                   rcode = pio_put_att(io_file, varid, "_FillValue", lfillvalue)
+                   rcode = pio_put_att(io_file,varid,"_FillValue",lfillvalue)
                 end if
                 if (NUOPC_FieldDictionaryHasEntry(trim(itemc))) then
                    call NUOPC_FieldDictionaryGetEntry(itemc, canonicalUnits=cunit, rc=rc)
@@ -988,7 +1030,10 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
                    if (tavg) then
                       rcode = pio_put_att(io_file, varid, "cell_methods", "time: mean")
                    endif
-                end if
+                endif
+             else
+                call shr_log_error(subname//' ERROR: unhandled rank', line=__LINE__, file=u_FILE_u, rc=rc)
+                return
              end if
           end if
        end do
@@ -1039,12 +1084,43 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
           end if
 
           call FB_getFldPtr(FB, itemc, &
-               fldptr1=fldptr1, fldptr2=fldptr2, rank=rank, rc=rc)
+               fldptr1=fldptr1, fldptr2=fldptr2, fldptr3=fldptr3, rank=rank, rc=rc)
           if (chkerr(rc,__LINE__,u_FILE_u)) return
 
           ! TODO (mvertens, 2019-03-13): this is a temporary mod to NOT write hgt
           if (trim(itemc) /= "hgt") then
-             if (rank == 2) then
+             if (rank == 3) then
+
+                ! Determine the size of the ungridded dimensions and the index where the distributed dimension is located
+                call ESMF_FieldBundleGet(FB, itemc,  field=lfield, rc=rc)
+                if (chkerr(rc,__LINE__,u_FILE_u)) return
+                call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound, gridToFieldMap=gridToFieldMap, rc=rc)
+                if (chkerr(rc,__LINE__,u_FILE_u)) return
+
+                if (gridToFieldMap(1) /= 3) then
+                   call shr_log_error( &
+                        subname//' ERROR: For rank-3 fields, currently only gridToFieldMap(1)==3 is supported', &
+                        line=__LINE__, file=u_FILE_u, rc=rc)
+                   return
+                end if
+
+                ! Output for each combination of ungriddedUbound indices
+                do n2 = 1,ungriddedUBound(2)
+                   do n = 1,ungriddedUBound(1)
+                      write(cnumber,'(i0)') n
+                      write(cnumber2,'(i0)') n2
+                      name1 = trim(lpre)//'_'//trim(itemc)//trim(cnumber)//'_'//trim(cnumber2)
+                      rcode = pio_inq_varid(io_file, trim(name1), varid)
+                      call pio_setframe(io_file,varid,frame)
+
+                      if (luse_float) then
+                         call pio_write_darray(io_file, varid, iodesc, real(fldptr3(n,n2,:),r4), rcode, fillval=real(lfillvalue,r4))
+                      else
+                         call pio_write_darray(io_file, varid, iodesc, fldptr3(n,n2,:), rcode, fillval=lfillvalue)
+                      end if
+                   end do
+                end do
+             else if (rank == 2) then
 
                 ! Determine the size of the ungridded dimension and the index where the undistributed dimension is located
                 call ESMF_FieldBundleGet(FB, itemc,  field=lfield, rc=rc)
@@ -1489,7 +1565,7 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc)
     type(ESMF_Field)              :: lfield
     integer                       :: rcode
     integer                       :: nf
-    integer                       :: k,n,l
+    integer                       :: k,n,n2,l
     type(file_desc_t)             :: pioid
     type(var_desc_t)              :: varid
     type(io_desc_t)               :: iodesc
@@ -1500,11 +1576,18 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc)
     integer                       :: rank, lsize
     real(r8), pointer             :: fldptr1(:), fldptr1_tmp(:)
     real(r8), pointer             :: fldptr2(:,:)
+    real(r8), pointer             :: fldptr3(:,:,:)
     character(CL)                 :: tmpstr
     character(len=16)             :: cnumber
+    character(len=16)             :: cnumber2
     integer(kind=Pio_Offset_Kind) :: lframe
-    integer                       :: ungriddedUBound(1) ! currently the size must equal 1 for rank 2 fieldds
-    integer                       :: gridToFieldMap(1)  ! currently the size must equal 1 for rank 2 fieldds
+
+    ! For a single ungridded dimension, there will be 1 element in ungriddedUBound and 1
+    ! element in gridToFieldMap; for two ungridded dimensions, there will be 2 elements in
+    ! ungriddedUBound but still 1 element in gridToFieldMap.
+    integer                       :: ungriddedUBound(2)
+    integer                       :: gridToFieldMap(1)
+
     character(*),parameter :: subName = '(med_io_read_FB) '
     !-------------------------------------------------------------------------------
     rc = ESMF_Success
@@ -1569,7 +1652,9 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc)
           if (chkerr(rc,__LINE__,u_FILE_u)) return
           call ESMF_FieldGet(lfield, rank=rank, rc=rc)
           if (chkerr(rc,__LINE__,u_FILE_u)) return
-          if (rank == 2) then
+          if (rank == 3) then
+             name1 = trim(lpre)//'_'//trim(itemc)//'1_1'
+          else if (rank == 2) then
              name1 = trim(lpre)//'_'//trim(itemc)//'1'
           else if (rank == 1) then
              name1 = trim(lpre)//'_'//trim(itemc)
@@ -1582,12 +1667,61 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc)
        if (chkerr(rc,__LINE__,u_FILE_u)) return
 
        ! Get pointer to field bundle field
-       ! Field bundle might be 2d or 1d - but field on mediator history or restart file will always be 1d
+       ! Field bundle might be 3d, 2d or 1d - but field on mediator history or restart file will always be 1d
        call FB_getFldPtr(FB, itemc, &
-            fldptr1=fldptr1, fldptr2=fldptr2, rank=rank, rc=rc)
+            fldptr1=fldptr1, fldptr2=fldptr2, fldptr3=fldptr3, rank=rank, rc=rc)
        if (chkerr(rc,__LINE__,u_FILE_u)) return
 
-       if (rank == 2) then
+       if (rank == 3) then
+
+          ! Determine the size of the ungridded dimensions and the
+          ! index where the distributed dimension is located
+          call ESMF_FieldBundleGet(FB, itemc,  field=lfield, rc=rc)
+          if (chkerr(rc,__LINE__,u_FILE_u)) return
+          call ESMF_FieldGet(lfield, ungriddedUBound=ungriddedUBound, gridToFieldMap=gridToFieldMap, rc=rc)
+          if (chkerr(rc,__LINE__,u_FILE_u)) return
+
+          if (gridToFieldMap(1) /= 3) then
+             call shr_log_error( &
+                  subname//' ERROR: For rank-3 fields, currently only gridToFieldMap(1)==3 is supported', &
+                  line=__LINE__, file=u_FILE_u, rc=rc)
+             return
+          end if
+
+          lsize = size(fldptr3, dim=3)
+          allocate(fldptr1_tmp(lsize))
+
+          do n2 = 1,ungriddedUBound(2)
+             do n = 1,ungriddedUBound(1)
+                ! Create a name for the 1d field on the mediator history or restart file based on the
+                ! ungridded dimension indices of the field bundle 3d field
+                write(cnumber,'(i0)') n
+                write(cnumber2,'(i0)') n2
+                name1 = trim(lpre)//'_'//trim(itemc)//trim(cnumber)//'_'//trim(cnumber2)
+
+                rcode = pio_inq_varid(pioid, trim(name1), varid)
+                if (rcode == pio_noerr) then
+                   call ESMF_LogWrite(trim(subname)//' read field '//trim(name1), ESMF_LOGMSG_INFO)
+                   if (chkerr(rc,__LINE__,u_FILE_u)) return
+                   call pio_setframe(pioid, varid, lframe)
+                   call pio_read_darray(pioid, varid, iodesc, fldptr1_tmp, rcode)
+                   rcode = pio_get_att(pioid, varid, "_FillValue", lfillvalue)
+                   if (rcode /= pio_noerr) then
+                      lfillvalue = fillvalue
+                   endif
+                   do l = 1,size(fldptr1_tmp)
+                      if (fldptr1_tmp(l) == lfillvalue) fldptr1_tmp(l) = 0.0_r8
+                   enddo
+                else
+                   fldptr1_tmp = 0.0_r8
+                endif
+                fldptr3(n,n2,:) = fldptr1_tmp(:)
+             end do
+          end do
+
+          deallocate(fldptr1_tmp)
+
+       else if (rank == 2) then
 
           ! Determine the size of the ungridded dimension and the
           ! index where the undistributed dimension is located