Fix dumpfields=true option by using ESMF_FieldBundleWrite

atmos_model.F90

@@ -815,7 +815,7 @@ subroutine atmos_model_init (Atmos, Time_init, Time, Time_step)
    if (mpp_pe() == mpp_root_pe()) print *,'in atmos_model, fhzero=',GFS_Control%fhzero, 'fhour=',sec/3600.,sec_lastfhzerofh/3600
 
    if (mod((sec-sec_lastfhzerofh),int(GFS_Control%fhzero*3600.)) /= 0) then
-     diag_time = Time - real_to_time_type(mod(int((GFS_Control%kdt - 1)*dt_phys-sec_lastfhzerofh),int(GFS_Control%fhzero))*3600.0)
+     diag_time = Time - real_to_time_type(real(mod(int((GFS_Control%kdt - 1)*dt_phys-sec_lastfhzerofh),int(GFS_Control%fhzero*3600.0))))
      if (mpp_pe() == mpp_root_pe()) print *,'Warning: in atmos_init,start at non multiple of fhzero'
    endif
    if (Atmos%iau_offset > zero) then

ccpp/data/GFS_typedefs.meta

@@ -9927,7 +9927,7 @@
   standard_name = auxiliary_2d_arrays
   long_name = auxiliary 2d arrays to output (for debugging)
   units = none
-  dimensions = (horizontal_loop_extent,number_of_xyz_dimensioned_auxiliary_arrays)
+  dimensions = (horizontal_loop_extent,number_of_xy_dimensioned_auxiliary_arrays)
   type = real
   kind = kind_phys
   active = (number_of_xy_dimensioned_auxiliary_arrays > 0)

cpl/module_cap_cpl.F90

@@ -17,27 +17,29 @@ module module_cap_cpl
   !-----------------------------------------------------------------------------
 
     subroutine diagnose_cplFields(gcomp, clock_fv3, fcstpe, &
-                                  statewrite_flag, stdiagnose_flag, state_tag)
+                                  statewrite_flag, stdiagnose_flag, state_tag, rc)
 
       type(ESMF_GridComp), intent(in)       :: gcomp
       type(ESMF_Clock),intent(in)           :: clock_fv3
       logical, intent(in)                   :: fcstpe
       logical, intent(in)                   :: statewrite_flag
       integer, intent(in)                   :: stdiagnose_flag
-      character(len=*),         intent(in)  :: state_tag                        !< Import or export.
+      character(len=*), intent(in)          :: state_tag                        !< "import" or "export".
+      integer, intent(out)                  :: rc
 
       character(len=*),parameter :: subname='(module_cap_cpl:diagnose_cplFields)'
       type(ESMF_Time) :: currTime
       type(ESMF_State) :: state
-      character(len=240) :: timestr
-      integer :: timeslice = 1
+      type(ESMF_TimeInterval) :: timeStep
+      character(len=240) :: import_timestr, export_timestr
       character(len=160) :: nuopcMsg
       character(len=160) :: filename
-      integer :: rc
 !
-      call ESMF_ClockGet(clock_fv3, currTime=currTime, rc=rc)
+      call ESMF_ClockGet(clock_fv3, currTime=currTime, timeStep=timestep, rc=rc)
+      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+      call ESMF_TimeGet(currTime, timestring=import_timestr, rc=rc)
       if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
-      call ESMF_TimeGet(currTime, timestring=timestr, rc=rc)
+      call ESMF_TimeGet(currTime+timestep, timestring=export_timestr, rc=rc)
       if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
 
       call ESMF_ClockPrint(clock_fv3, options="currTime", preString="current time: ", unit=nuopcMsg)
@@ -53,8 +55,8 @@ subroutine diagnose_cplFields(gcomp, clock_fv3, fcstpe, &
 
         ! Dump Fields out
         if (statewrite_flag) then
-          write(filename,'(A)') 'fv3_cap_import_'//trim(timestr)//'_'
-          call State_RWFields_tiles(state,trim(filename), timeslice, rc=rc)
+          write(filename,'(A)') 'fv3_cap_import_'//trim(import_timestr)//'.tile*.nc'
+          call State_RWFields_tiles(state,trim(filename), rc=rc)
           if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
         end if
       end if
@@ -69,8 +71,8 @@ subroutine diagnose_cplFields(gcomp, clock_fv3, fcstpe, &
 
         ! Dump Fields out
         if (statewrite_flag) then
-          write(filename,'(A)') 'fv3_cap_export_'//trim(timestr)//'_'
-          call State_RWFields_tiles(state,trim(filename), timeslice, rc=rc)
+          write(filename,'(A)') 'fv3_cap_export_'//trim(export_timestr)//'.tile*.nc'
+          call State_RWFields_tiles(state,trim(filename), rc=rc)
           if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
         end if
       end if
@@ -80,28 +82,36 @@ end subroutine diagnose_cplFields
   !-----------------------------------------------------------------------------
 
   ! This subroutine requires ESMFv8 - for coupled FV3
-    subroutine State_RWFields_tiles(state,filename,timeslice,rc)
+    subroutine State_RWFields_tiles(state,filename,rc)
 
       type(ESMF_State), intent(in)          :: state
       character(len=*), intent(in)          :: fileName
-      integer, intent(in)                   :: timeslice
       integer, intent(out)                  :: rc
 
-      ! local
-      type(ESMF_Field)                       :: firstESMFFLD
-      type(ESMF_Field),allocatable           :: flds(:)
-      type(ESMF_GridComp) :: IOComp
-      type(ESMF_Grid) :: gridFv3
-
-      character(len=256) :: msgString
-      integer                                :: i, icount, ifld
+      ! local variables
+      type(ESMF_Array)                       :: array
+      type(ESMF_Grid)                        :: grid
+      type(ESMF_FieldBundle)                 :: fieldbundle
+      type(ESMF_Field), allocatable          :: flds(:)
+      type(ESMF_DistGrid)                    :: distgrid
+      integer                                :: i, icount, ifld, id
       integer                                :: fieldcount, firstfld
+      integer                                :: fieldDimCount, gridDimCount, dimCount, tileCount, ungriddedDimCount
       character(64), allocatable             :: itemNameList(:), fldNameList(:)
       type(ESMF_StateItem_Flag), allocatable :: typeList(:)
+      integer, allocatable                   :: minIndexPTile(:,:), maxIndexPTile(:,:)
+      integer, allocatable                   :: ungriddedLBound(:), ungriddedUBound(:)
+      integer, allocatable                   :: fieldDimLen(:)
+      character(len=32), allocatable         :: gridded_dim_labels(:), ungridded_dim_labels(:)
 
-      character(len=*),parameter :: subname='(module_cap_cpl:State_RWFields_tiles)'
+      character(16), parameter :: convention = 'NetCDF'
+      character(16), parameter :: purpose = 'FV3'
 
-      ! local variables
+      integer, parameter :: max_n_axes = 4
+      integer, parameter :: max_n_dim = 16
+      integer, dimension(max_n_axes, max_n_dim) :: axes_dimcount = 0
+
+      character(len=*),parameter :: subname='(module_cap_cpl:State_RWFields_tiles)'
 
       rc = ESMF_SUCCESS
       !call ESMF_LogWrite(trim(subname)//trim(filename)//": called", ESMF_LOGMSG_INFO, rc=rc)
@@ -118,9 +128,6 @@ subroutine State_RWFields_tiles(state,filename,timeslice,rc)
         if(typeList(i) == ESMF_STATEITEM_FIELD) firstfld = i
         if(typeList(i) == ESMF_STATEITEM_FIELD) fieldcount = fieldcount + 1
       enddo
-      !write(msgString,*) trim(subname)//' icount = ',icount," fieldcount =
-      !",fieldcount," firstfld = ",firstfld
-      !call ESMF_LogWrite(trim(msgString), ESMF_LOGMSG_INFO, rc=rc)
 
       allocate(flds(fieldCount),fldNameList(fieldCount))
       ifld = 1
@@ -131,37 +138,152 @@ subroutine State_RWFields_tiles(state,filename,timeslice,rc)
         endif
       enddo
 
-      call ESMF_LogWrite(trim(subname)//": write "//trim(filename)//"tile1-tile6", ESMF_LOGMSG_INFO, rc=rc)
-      ! get first field
-      call ESMF_StateGet(state, itemName=itemNameList(firstfld), field=firstESMFFLD, rc=rc)
-      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
-          line=__LINE__, file=__FILE__)) return  ! bail out
-
-      call ESMF_FieldGet(firstESMFFLD, grid=gridFv3, rc=rc)
-      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
-          line=__LINE__, file=__FILE__)) return  ! bail out
+      fieldbundle = ESMF_FieldBundleCreate(rc=rc)
+      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
 
-      IOComp = ESMFIO_Create(gridFv3, rc=rc)
-      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
-          line=__LINE__, file=__FILE__)) return  ! bail out
+      call ESMF_LogWrite(trim(subname)//": write "//trim(filename), ESMF_LOGMSG_INFO, rc=rc)
 
       do ifld=1, fieldCount
         call ESMF_StateGet(state, itemName=fldNameList(ifld), field=flds(ifld), rc=rc)
+        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+        call ESMF_FieldGet(flds(ifld), grid=grid, dimCount=fieldDimCount, array=array, rc=rc)
+        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+        if (fieldDimCount > 4) then
+          call ESMF_LogWrite(trim(subname)//": fieldDimCount > 4 unsupported", ESMF_LOGMSG_ERROR, rc=rc)
+        end if
+
+        call ESMF_GridGet(grid, dimCount=gridDimCount, rc=rc)
+        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+        if (gridDimCount > 2) then
+          call ESMF_LogWrite(trim(subname)//": gridDimCount > 2 unsupported", ESMF_LOGMSG_ERROR, rc=rc)
+        end if
+
+        call ESMF_ArrayGet(array, distgrid=distgrid, dimCount=dimCount, tileCount=tileCount, rc=rc)
+        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+        ! skip 'cpl_scalars' field because it has tileCount == 1, while all other fields have 6.
+        ! This causes the following error:
+        ! 20240705 134459.788 ERROR            PET000 ESMCI_IO.C:1614 ESMCI::IO::checkNtiles() Wrong data value  - New number of tiles (6) does not match previously-set number of tiles (1) for this IO object. All arrays handled by a given IO object must have the same number of tiles.
+        if (trim(fldNameList(ifld)) == 'cpl_scalars') then
+          cycle
+        endif
+
+        allocate(fieldDimLen(fieldDimCount))
+
+        allocate(minIndexPTile(dimCount, tileCount))
+        allocate(maxIndexPTile(dimCount, tileCount))
+        call ESMF_DistGridGet(distgrid, minIndexPTile=minIndexPTile, maxIndexPTile=maxIndexPTile, rc=rc)
+        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+        allocate(gridded_dim_labels(gridDimCount))
+        do i = 1, gridDimCount
+          fieldDimLen(i) = maxIndexPTile(i,1) - minIndexPTile(i,1) + 1
+          id = find_axis_id_for_axis_count(i,fieldDimLen(i))
+          if (id < 1) then
+            call ESMF_LogWrite(trim(subname)//": id < 1", ESMF_LOGMSG_ERROR, rc=rc)
+          endif
+          if (i == 1) write(gridded_dim_labels(i),'(A,I0)') 'xaxis_',id
+          if (i == 2) write(gridded_dim_labels(i),'(A,I0)') 'yaxis_',id
+        end do
+
+        deallocate(minIndexPTile)
+        deallocate(maxIndexPTile)
+
+        ungriddedDimCount = fieldDimCount - gridDimCount
+        allocate(ungridded_dim_labels(ungriddedDimCount))
+        if (fieldDimCount > gridDimCount) then
+          allocate(ungriddedLBound(ungriddedDimCount))
+          allocate(ungriddedUBound(ungriddedDimCount))
+          call ESMF_FieldGet(flds(ifld), ungriddedLBound=ungriddedLBound, ungriddedUBound=ungriddedUBound, rc=rc)
+          if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+          do i=1,ungriddedDimCount
+            fieldDimLen(i+gridDimCount) = ungriddedUBound(i) - ungriddedLBound(i) + 1
+            id = find_axis_id_for_axis_count(i+gridDimCount, fieldDimLen(i+gridDimCount))
+            if (id < 1) then
+              write(0,*)'stop error', id, i, fieldDimLen(i+gridDimCount)
+            endif
+            if (i==1) write(ungridded_dim_labels(i),'(A,I0)') 'zaxis_',id
+            if (i==2) write(ungridded_dim_labels(i),'(A,I0)') 'taxis_',id
+          end do
+          deallocate(ungriddedLBound)
+          deallocate(ungriddedUBound)
+        end if
+
+        call ESMF_AttributeAdd(grid, convention=convention, purpose=purpose, attrList=(/ ESMF_ATT_GRIDDED_DIM_LABELS /), rc=rc)
+        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+        call ESMF_AttributeSet(grid, convention=convention, purpose=purpose, &
+                               name=ESMF_ATT_GRIDDED_DIM_LABELS, valueList=gridded_dim_labels, rc=rc)
+        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+        if (ungriddedDimCount > 0) then
+          call ESMF_AttributeAdd(flds(ifld), convention=convention, purpose=purpose, &
+                                 attrList=(/ ESMF_ATT_UNGRIDDED_DIM_LABELS /), rc=rc)
+          if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+          call ESMF_AttributeSet(flds(ifld), convention=convention, purpose=purpose, &
+                                 name=ESMF_ATT_UNGRIDDED_DIM_LABELS, valueList=ungridded_dim_labels, rc=rc)
+          if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+        end if
+
+        deallocate(fieldDimLen)
+        deallocate(gridded_dim_labels)
+        deallocate(ungridded_dim_labels)
+
+        call ESMF_FieldBundleAdd(fieldbundle, (/flds(ifld)/), rc=rc)
+        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
       enddo
 
-      call ESMFIO_Write(IOComp, filename, flds, filePath='./', rc=rc)
-      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
-          line=__LINE__, file=__FILE__)) return  ! bail out
+      call ESMF_FieldBundleWrite(fieldbundle, fileName=trim(filename), convention=convention, purpose=purpose, &
+           timeslice=1, overwrite=.true., rc=rc)
+      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
 
-! -- Finalize ESMFIO
+! -- Finalize
       deallocate(flds)
       deallocate(fldNameList)
-      call ESMFIO_Destroy(IOComp, rc=rc)
-      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
-      line=__LINE__, file=__FILE__)) call ESMF_Finalize()
+
+      call ESMF_FieldBundleDestroy(fieldbundle, rc=rc)
+      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
 
       !call ESMF_LogWrite(trim(subname)//trim(filename)//": finished", ESMF_LOGMSG_INFO, rc=rc)
 
+    contains
+
+      function find_axis_id_for_axis_count(axis, count) result(id)
+        integer, intent(in) :: axis, count
+
+        integer :: id
+        integer :: i
+
+        id = -1 ! not found
+
+        if (axis > max_n_axes) then
+          call ESMF_LogWrite('axis > max_n_axes. Increase max_n_axes in '//trim(subname), ESMF_LOGMSG_ERROR)
+          return
+        end if
+
+        do i =1, max_n_dim
+           if (axes_dimcount(axis, i) == 0) then
+             axes_dimcount(axis, i) = count
+             id = i
+             return
+           else
+             if (axes_dimcount(axis, i) == count) then
+                 id = i
+                 return
+             end if
+           end if
+        end do
+
+        call ESMF_LogWrite('Increase max_n_dim in '//trim(subname), ESMF_LOGMSG_ERROR)
+
+      end function find_axis_id_for_axis_count
+
     end subroutine State_RWFields_tiles
 
   !-----------------------------------------------------------------------------

cpl/module_cplscalars.F90

@@ -116,7 +116,7 @@ subroutine State_SetScalar(scalar_value, scalar_id, State, flds_scalar_name, fld
     if (mytask == 0) then
        call ESMF_FieldGet(lfield, farrayPtr = farrayptr, rc=rc)
       if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
-       if (scalar_id < 0 .or. scalar_id > flds_scalar_num) then
+       if (scalar_id <= 0 .or. scalar_id > flds_scalar_num) then
           call ESMF_LogWrite(trim(subname)//": ERROR in scalar_id", ESMF_LOGMSG_INFO)
           rc = ESMF_FAILURE
           return

fv3_cap.F90

@@ -288,7 +288,6 @@ subroutine InitializeAdvertise(gcomp, rc)
        call ESMF_LogWrite(trim(subname)//' flds_scalar_name = '//trim(flds_scalar_name), ESMF_LOGMSG_INFO)
        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
     endif
-
     call NUOPC_CompAttributeGet(gcomp, name="ScalarFieldCount", value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc)
     if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
     if (isPresent .and. isSet) then
@@ -313,14 +312,23 @@ subroutine InitializeAdvertise(gcomp, rc)
        call ESMF_LogWrite(trim(subname)//' : flds_scalar_index_ny = '//trim(msgString), ESMF_LOGMSG_INFO)
        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
     endif
-    call NUOPC_CompAttributeGet(gcomp, name="ScalarFieldIdxGridNTile", value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc)
-    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
-    if (isPresent .and. isSet) then
-       read(cvalue,*) flds_scalar_index_ntile
-       write(msgString,*) flds_scalar_index_ntile
-       call ESMF_LogWrite(trim(subname)//' : flds_scalar_index_ntile = '//trim(msgString), ESMF_LOGMSG_INFO)
+    ! tile index must be present if indices for nx and ny are non-zero
+    if (flds_scalar_index_nx /= 0 .and. flds_scalar_index_ny /=0 ) then
+       call NUOPC_CompAttributeGet(gcomp, name="ScalarFieldIdxGridNTile", isPresent=isPresent, isSet=isSet, rc=rc)
        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
-    endif
+       if (.not. isPresent .and. .not. isSet) then
+          if (mype == 0)write(*,*)'ERROR : ScalarFieldIdxGridNTile must be set'
+          call ESMF_LogWrite('ERROR : ScalarFieldIdxGridNTile must be set', ESMF_LOGMSG_ERROR)
+          rc = ESMF_FAILURE
+          if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+       else
+          call NUOPC_CompAttributeGet(gcomp, name="ScalarFieldIdxGridNTile", value=cvalue, rc=rc)
+          if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+          read(cvalue,*) flds_scalar_index_ntile
+          write(msgString,*) flds_scalar_index_ntile
+          call ESMF_LogWrite(trim(subname)//' : flds_scalar_index_ntile = '//trim(msgString), ESMF_LOGMSG_INFO)
+       endif
+    end if
 
 !------------------------------------------------------------------------
 ! get config variables
@@ -1081,7 +1089,8 @@ subroutine ModelAdvance_phase1(gcomp, rc)
     if( dbug > 0 .or. cplprint_flag ) then
          fcstpe = .false.
          if( mype < num_pes_fcst ) fcstpe = .true.
-         call diagnose_cplFields(gcomp, clock, fcstpe, cplprint_flag, dbug, 'import')
+         call diagnose_cplFields(gcomp, clock, fcstpe, cplprint_flag, dbug, 'import', rc=rc)
+         if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
     endif
 
     timep1re = MPI_Wtime()
@@ -1235,7 +1244,8 @@ subroutine ModelAdvance_phase2(gcomp, rc)
     if( dbug > 0 .or. cplprint_flag ) then
       fcstpe = .false.
       if( mype < num_pes_fcst ) fcstpe = .true.
-      call diagnose_cplFields(gcomp, clock_out, fcstpe, cplprint_flag, dbug, 'export')
+      call diagnose_cplFields(gcomp, clock_out, fcstpe, cplprint_flag, dbug, 'export', rc=rc)
+      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
     end if
 
     timep2re = MPI_Wtime()