Sync w/ ESCOMP, add cpl_scalars for CSG and regional ATM domains

cesm/driver/esm.F90

@@ -1224,13 +1224,17 @@ subroutine esm_set_single_column_attributes(compname, gcomp, rc)
     real (r8), allocatable :: lats(:)        ! temporary
     real (r8), allocatable :: lons(:)        ! temporary
     real (r8), allocatable :: pos_lons(:)    ! temporary
+    real (r8), allocatable :: pos_lats(:)    ! temporary
+    real (r8), allocatable :: cols(:)        ! temporary
     real (r8), allocatable :: glob_grid(:,:) ! temporary
     real (r8)              :: pos_scol_lon   ! temporary
+    real (r8)              :: pos_scol_lat   ! temporary
     real (r8)              :: scol_data(1)
     integer                :: iscol_data(1)
     integer                :: petcount
     character(len=CL)      :: cvalue
     character(len=*), parameter :: subname= ' (esm_get_single_column_attributes) '
+    logical                :: unstructured = .false.
     !-------------------------------------------------------------------------------
 
     rc = ESMF_SUCCESS
@@ -1324,7 +1328,15 @@ subroutine esm_set_single_column_attributes(compname, gcomp, rc)
           if (status /= nf90_noerr) call shr_sys_abort (subname//' inq_varid frac')
 
           ! Read in domain file for single column
-          allocate(lats(nj))
+          ! Check for unstructured data ni>1 and nj==1
+          if (ni.gt.1 .and. nj == 1) unstructured=.true.
+
+          if (unstructured) then
+             allocate(lats(ni))
+             allocate(pos_lats(ni))
+          else
+             allocate(lats(nj))
+          end if
           allocate(lons(ni))
           allocate(pos_lons(ni))
           allocate(glob_grid(ni,nj))
@@ -1334,28 +1346,37 @@ subroutine esm_set_single_column_attributes(compname, gcomp, rc)
           count3=(/ni,nj,1/)
           status = nf90_get_var(ncid, varid_xc, glob_grid, start3, count3)
           if (status /= nf90_noerr) call shr_sys_abort (subname//' get_var xc')
-          do i = 1,ni
-             lons(i) = glob_grid(i,1)
-          end do
+          lons(1:ni) = glob_grid(1:ni,1)
           status = nf90_get_var(ncid, varid_yc, glob_grid, start3, count3)
           if (status /= nf90_noerr) call shr_sys_abort (subname//' get_var yc')
-          do j = 1,nj
-             lats(j) = glob_grid(1,j)
-          end do
-
+          if (unstructured) then
+             lats(1:ni) = glob_grid(1:ni,1)
+          else
+             lats(1:nj) = glob_grid(1,1:nj)
+          end if
           ! find nearest neighbor indices of scol_lon and scol_lat in single_column_lnd_domain file
           ! convert lons array and scol_lon to 0,360 and find index of value closest to 0
           ! and obtain single-column longitude/latitude indices to retrieve
-          pos_lons(:)  = mod(lons(:)  + 360._r8, 360._r8)
-          pos_scol_lon = mod(scol_lon + 360._r8, 360._r8)
-          start(1) = (MINLOC(abs(pos_lons - pos_scol_lon), dim=1))
-          start(2) = (MINLOC(abs(lats      -scol_lat    ), dim=1))
-
+          if (unstructured) then
+             allocate(cols(ni))
+             pos_lons(:)  = mod(lons(:)  + 360._r8, 360._r8)
+             pos_scol_lon = mod(scol_lon + 360._r8, 360._r8)
+             pos_lats(:)  = lats(:)  + 90._r8
+             pos_scol_lat = scol_lat + 90._r8
+             cols=abs(pos_lons - pos_scol_lon)+abs(pos_lats - pos_scol_lat)
+             start(1) = MINLOC(cols, dim=1)
+             start(2) = 1
+             deallocate(cols)
+          else
+             pos_lons(:)  = mod(lons(:)  + 360._r8, 360._r8)
+             pos_scol_lon = mod(scol_lon + 360._r8, 360._r8)
+             start(1) = (MINLOC(abs(pos_lons - pos_scol_lon), dim=1))
+             start(2) = (MINLOC(abs(lats      -scol_lat    ), dim=1))
+          end if
           deallocate(lats)
           deallocate(lons)
           deallocate(pos_lons)
           deallocate(glob_grid)
-
           ! read in value of nearest neighbor lon and RESET scol_lon and scol_lat
           ! also get area of gridcell, mask and frac
           status = nf90_get_var(ncid, varid_xc, scol_lon, start)

mediator/med.F90

@@ -510,7 +510,7 @@ subroutine SetServices(gcomp, rc)
 
 #ifdef CDEPS_INLINE
     !------------------
-    ! phase routine for cdeps inline capabilty 
+    ! phase routine for cdeps inline capability
     !------------------
 
     call NUOPC_CompSetEntryPoint(gcomp, ESMF_METHOD_RUN, &
@@ -832,10 +832,10 @@ subroutine AdvertiseFields(gcomp, importState, exportState, clock, rc)
     if (trim(coupling_mode) == 'cesm') then
        call esmFldsExchange_cesm(gcomp, phase='advertise', rc=rc)
        if (ChkErr(rc,__LINE__,u_FILE_u)) return
-    else if (trim(coupling_mode(1:3)) == 'ufs') then
+    else if (coupling_mode(1:3) == 'ufs') then
        call esmFldsExchange_ufs(gcomp, phase='advertise', rc=rc)
        if (ChkErr(rc,__LINE__,u_FILE_u)) return
-    else if (trim(coupling_mode(1:4)) == 'hafs') then
+    else if (coupling_mode(1:4) == 'hafs') then
        call esmFldsExchange_hafs(gcomp, phase='advertise', rc=rc)
        if (ChkErr(rc,__LINE__,u_FILE_u)) return
     else
@@ -867,13 +867,22 @@ subroutine AdvertiseFields(gcomp, importState, exportState, clock, rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
     read(cvalue,*) is_local%wrap%flds_scalar_index_ny
 
+    call NUOPC_CompAttributeGet(gcomp, name="ScalarFieldIdxGridNTile", value=cvalue, &
+         isPresent=isPresent, isSet=isSet,rc=rc)
+    if (ChkErr(rc,__LINE__,u_FILE_u)) return
+    if (isPresent .and. isSet) then
+       read(cvalue,*) is_local%wrap%flds_scalar_index_ntile
+    else
+       is_local%wrap%flds_scalar_index_ntile = 0
+    end if
+
     call NUOPC_CompAttributeGet(gcomp, name="ScalarFieldIdxNextSwCday", value=cvalue, &
          isPresent=isPresent, isSet=isSet, rc=rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
     if (isPresent .and. isSet) then
        read(cvalue,*) is_local%wrap%flds_scalar_index_nextsw_cday
     else
-       is_local%wrap%flds_scalar_index_nextsw_cday = spval
+       is_local%wrap%flds_scalar_index_nextsw_cday = 0
     end if
 
     call NUOPC_CompAttributeGet(gcomp, name="ScalarFieldIdxPrecipFactor", value=cvalue, isPresent=isPresent, isSet=isSet, rc=rc)
@@ -962,7 +971,7 @@ subroutine AdvertiseFields(gcomp, importState, exportState, clock, rc)
           endif
           if (maintask) then
              write(logunit,*) trim(compname(ncomp))//'_use_data_first_import is ', is_local%wrap%med_data_force_first(ncomp)
-          endif    
+          endif
        end if
     end do
 
@@ -1067,7 +1076,7 @@ subroutine ModifyDecompofMesh(gcomp, importState, exportState, clock, rc)
     if (ChkErr(rc,__LINE__,u_FILE_u)) return
 
     !------------------
-    ! Recieve Grids
+    ! Receive Grids
     !------------------
 
     do n1 = 1,ncomps
@@ -1644,7 +1653,7 @@ subroutine DataInitialize(gcomp, rc)
     logical                            :: read_restart
     logical                            :: allDone = .false.
     logical,save                       :: first_call = .true.
-    real(r8)                           :: real_nx, real_ny
+    real(r8)                           :: real_nx, real_ny, real_ntile
     character(len=CX)                  :: msgString
     character(len=*), parameter :: subname = '('//__FILE__//':DataInitialize)'
     !-----------------------------------------------------------
@@ -1832,7 +1841,7 @@ subroutine DataInitialize(gcomp, rc)
       if (trim(coupling_mode) == 'cesm') then
          call esmFldsExchange_cesm(gcomp, phase='initialize', rc=rc)
          if (ChkErr(rc,__LINE__,u_FILE_u)) return
-      else if (trim(coupling_mode(1:3)) == 'ufs') then
+      else if (coupling_mode(1:3) == 'ufs') then
          call esmFldsExchange_ufs(gcomp, phase='initialize', rc=rc)
          if (ChkErr(rc,__LINE__,u_FILE_u)) return
       else if (coupling_mode(1:4) == 'hafs') then
@@ -2128,11 +2137,22 @@ subroutine DataInitialize(gcomp, rc)
                   flds_scalar_name=is_local%wrap%flds_scalar_name, &
                   flds_scalar_num=is_local%wrap%flds_scalar_num, rc=rc)
              if (ChkErr(rc,__LINE__,u_FILE_u)) return
+             if (is_local%wrap%flds_scalar_index_ntile > 0) then
+                call State_GetScalar(scalar_value=real_ntile, &
+                     scalar_id=is_local%wrap%flds_scalar_index_ntile, &
+                     state=is_local%wrap%NstateImp(n1), &
+                     flds_scalar_name=is_local%wrap%flds_scalar_name, &
+                     flds_scalar_num=is_local%wrap%flds_scalar_num, rc=rc)
+                if (ChkErr(rc,__LINE__,u_FILE_u)) return
+                is_local%wrap%ntile(n1) = nint(real_ntile)
+             else
+                is_local%wrap%ntile(n1) = 0
+             end if
              is_local%wrap%nx(n1) = nint(real_nx)
              is_local%wrap%ny(n1) = nint(real_ny)
-             write(msgString,'(2i8,2l4)') is_local%wrap%nx(n1), is_local%wrap%ny(n1)
+             write(msgString,'(3i8)') is_local%wrap%nx(n1), is_local%wrap%ny(n1), is_local%wrap%ntile(n1)
              if (maintask) then
-                write(logunit,'(a)') 'global nx,ny sizes for '//trim(compname(n1))//":"//trim(msgString)
+                write(logunit,'(a)') 'global nx,ny,ntile sizes for '//trim(compname(n1))//":"//trim(msgString)
              end if
              call ESMF_LogWrite(trim(subname)//":"//trim(compname(n1))//":"//trim(msgString), ESMF_LOGMSG_INFO)
           end if

mediator/med_fraction_mod.F90

@@ -293,7 +293,7 @@ subroutine med_fraction_init(gcomp, rc)
           ! If ice and atm are on the same mesh - a redist route handle has already been created
           maptype = mapfcopy
        else
-          if (trim(coupling_mode(1:9)) == 'ufs.nfrac' ) then
+          if (coupling_mode(1:9) == 'ufs.nfrac' ) then
              maptype = mapnstod_consd
           else
              maptype = mapconsd
@@ -345,7 +345,7 @@ subroutine med_fraction_init(gcomp, rc)
           ! If ocn and atm are on the same mesh - a redist route handle has already been created
           maptype = mapfcopy
        else
-          if (trim(coupling_mode(1:9)) == 'ufs.nfrac' ) then
+          if (coupling_mode(1:9) == 'ufs.nfrac' ) then
              maptype = mapnstod_consd
           else
              maptype = mapconsd
@@ -756,7 +756,7 @@ subroutine med_fraction_set(gcomp, rc)
 
           call t_startf('MED:'//trim(subname)//' fbfrac(compatm)')
           ! Determine maptype
-          if (trim(coupling_mode(1:9)) == 'ufs.nfrac' ) then
+          if (coupling_mode(1:9) == 'ufs.nfrac' ) then
              maptype = mapnstod_consd
           else
              if (med_map_RH_is_created(is_local%wrap%RH(compice,compatm,:),mapfcopy, rc=rc)) then

mediator/med_internalstate_mod.F90

@@ -121,7 +121,7 @@ module med_internalstate_mod
     ! Present/allowed coupling/active coupling logical flags
     logical, pointer :: comp_present(:)               ! comp present flag
     logical, pointer :: med_coupling_active(:,:)      ! computes the active coupling
-    logical, pointer :: med_data_active(:,:)          ! uses stream data to provide background fill 
+    logical, pointer :: med_data_active(:,:)          ! uses stream data to provide background fill
     logical, pointer :: med_data_force_first(:)       ! force to use stream data for first coupling timestep
     integer          :: num_icesheets                 ! obtained from attribute
     logical          :: ocn2glc_coupling = .false.    ! obtained from attribute
@@ -133,12 +133,15 @@ module med_internalstate_mod
 
     ! Global nx,ny dimensions of input arrays (needed for mediator history output)
     integer, pointer   :: nx(:), ny(:)
+    ! Number of nx*ny domains (needed for cubed-sphere and regional domains)
+    integer, pointer   :: ntile(:)
 
     ! Import/Export Scalars
     character(len=CL) :: flds_scalar_name = ''
     integer           :: flds_scalar_num = 0
     integer           :: flds_scalar_index_nx = 0
     integer           :: flds_scalar_index_ny = 0
+    integer           :: flds_scalar_index_ntile = 0
     integer           :: flds_scalar_index_nextsw_cday = 0
     integer           :: flds_scalar_index_precip_factor = 0
     real(r8)          :: flds_scalar_precip_factor = 1._r8  ! actual value of precip factor from ocn
@@ -312,6 +315,7 @@ subroutine med_internalstate_init(gcomp, rc)
     allocate(is_local%wrap%med_data_force_first(ncomps))
     allocate(is_local%wrap%nx(ncomps))
     allocate(is_local%wrap%ny(ncomps))
+    allocate(is_local%wrap%ntile(ncomps))
     allocate(is_local%wrap%NStateImp(ncomps))
     allocate(is_local%wrap%NStateExp(ncomps))
     allocate(is_local%wrap%FBImp(ncomps,ncomps))
@@ -601,7 +605,7 @@ subroutine med_internalstate_defaultmasks(gcomp, rc)
     if (is_local%wrap%comp_present(compocn)) defaultMasks(compocn,:) = 0
     if (is_local%wrap%comp_present(compice)) defaultMasks(compice,:) = 0
     if (is_local%wrap%comp_present(compwav)) defaultMasks(compwav,:) = 0
-    if ( trim(coupling_mode(1:3)) == 'ufs') then
+    if ( coupling_mode(1:3) == 'ufs') then
        if (is_local%wrap%comp_present(compatm)) defaultMasks(compatm,:) = 1
     endif
     if ( trim(coupling_mode) == 'hafs') then

mediator/med_io_mod.F90

@@ -698,7 +698,7 @@ end function med_io_sec2hms
 
   !===============================================================================
   subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
-       fillval, pre, flds, tavg, use_float, tilesize, rc)
+       fillval, pre, flds, tavg, use_float, ntile, rc)
 
     !---------------
     ! Write FB to netcdf file
@@ -728,7 +728,7 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
     character(len=*), optional , intent(in) :: flds(:)   ! specific fields to write out
     logical,          optional , intent(in) :: tavg      ! is this a tavg
     logical,          optional , intent(in) :: use_float ! write output as float rather than double
-    integer,          optional , intent(in) :: tilesize  ! if non-zero, write atm component on tiles
+    integer,          optional , intent(in) :: ntile     ! number of nx * ny tiles
     integer                    , intent(out):: rc
 
     ! local variables
@@ -754,7 +754,7 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
     character(CS)                 :: coordvarnames(2)   ! coordinate variable names
     character(CS)                 :: coordnames(2)      ! coordinate long names
     character(CS)                 :: coordunits(2)      ! coordinate units
-    integer                       :: lnx,lny
+    integer                       :: lnx,lny,lntile
     logical                       :: luse_float
     real(r8)                      :: lfillvalue
     integer, pointer              :: minIndexPTile(:,:)
@@ -770,8 +770,7 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
     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                       :: atmtiles
-    integer                       :: ntiles = 1
+    logical                       :: tiles
     character(CL), allocatable    :: fieldNameList(:)
     character(*),parameter :: subName = '(med_io_write_FB) '
     !-------------------------------------------------------------------------------
@@ -785,9 +784,9 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
     luse_float = .false.
     if (present(use_float)) luse_float = use_float
 
-    atmtiles = .false.
-    if (present(tilesize)) then
-      if (tilesize > 0) atmtiles = .true.
+    tiles = .false.
+    if (present(ntile)) then
+      if (ntile > 0) tiles = .true.
     end if
 
     ! Error check
@@ -870,14 +869,14 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
     ! all the global grid values in the distgrid - e.g. CTSM
 
     ng = maxval(maxIndexPTile)
-    if (atmtiles) then
-      lnx = tilesize
-      lny = tilesize
-      ntiles = ng/(lnx*lny)
-      write(tmpstr,*) subname, 'ng,lnx,lny,ntiles = ',ng,lnx,lny,ntiles
+    if (tiles) then
+      lnx = nx
+      lny = ny
+      lntile = ng/(lnx*lny)
+      write(tmpstr,*) subname, 'ng,lnx,lny,lntile = ',ng,lnx,lny,lntile
       call ESMF_LogWrite(trim(tmpstr), ESMF_LOGMSG_INFO)
-      if (ntiles /= 6) then
-         call ESMF_LogWrite(trim(subname)//' ERROR: only cubed sphere atm tiles valid ', ESMF_LOGMSG_INFO)
+      if (lntile /= ntile) then
+         call ESMF_LogWrite(trim(subname)//' ERROR: grid2d size and ntile are not consistent ', ESMF_LOGMSG_INFO)
          call ESMF_Finalize(endflag=ESMF_END_ABORT)
       endif
     else
@@ -900,10 +899,10 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
 
     ! Write header
     if (whead) then
-      if (atmtiles) then
+      if (tiles) then
        rcode = pio_def_dim(io_file, trim(lpre)//'_nx', lnx, dimid3(1))
        rcode = pio_def_dim(io_file, trim(lpre)//'_ny', lny, dimid3(2))
-       rcode = pio_def_dim(io_file, trim(lpre)//'_ntiles', ntiles, dimid3(3))
+       rcode = pio_def_dim(io_file, trim(lpre)//'_ntile', ntile, dimid3(3))
        if (present(nt)) then
           dimid4(1:3) = dimid3
           rcode = pio_inq_dimid(io_file, 'time', dimid4(4))
@@ -1020,8 +1019,8 @@ subroutine med_io_write_FB(io_file, FB, whead, wdata, nx, ny, nt, &
        call ESMF_DistGridGet(distgrid, localDE=0, seqIndexList=dof, rc=rc)
        write(tmpstr,*) subname,' dof = ',ns,size(dof),dof(1),dof(ns)  !,minval(dof),maxval(dof)
        call ESMF_LogWrite(trim(tmpstr), ESMF_LOGMSG_INFO)
-       if (atmtiles) then
-          call pio_initdecomp(io_subsystem, pio_double, (/lnx,lny,ntiles/), dof, iodesc)
+       if (tiles) then
+          call pio_initdecomp(io_subsystem, pio_double, (/lnx,lny,ntile/), dof, iodesc)
        else
           call pio_initdecomp(io_subsystem, pio_double, (/lnx,lny/), dof, iodesc)
          !call pio_writedof(lpre, (/lnx,lny/), int(dof,kind=PIO_OFFSET_KIND), mpicom)
@@ -1579,8 +1578,8 @@ subroutine med_io_read_FB(filename, vm, FB, pre, frame, rc)
           allocate(fldptr1_tmp(lsize))
 
           do n = 1,ungriddedUBound(1)
-             ! Creat a name for the 1d field on the mediator history or restart file based on the
-             ! ungridded dimension index of the field bundle 2d fiedl
+             ! Create a name for the 1d field on the mediator history or restart file based on the
+             ! ungridded dimension index of the field bundle 2d field
              write(cnumber,'(i0)') n
              name1 = trim(lpre)//'_'//trim(itemc)//trim(cnumber)