Skip to content

Commit

Permalink
Sync w/ ESCOMP, add cpl_scalars for CSG and regional ATM domains (#115)
Browse files Browse the repository at this point in the history 
* add cpl_scalar for tiled grids, other minor fixes
* add new cpl_scalar for mediator history files for tiled gridded domains
* remove unnecessary trims, fix minor typos and indentation
* set ntile=0 when ntile scalar doesn't exist
* modify dstmask for lnd->atm in UFS


Co-authored-by: uturuncoglu <[email protected]>
  • Loading branch information
@DeniseWorthen @uturuncoglu
DeniseWorthen and uturuncoglu authored Apr 3, 2024
1 parent 758491e commit 4e19850
Show file tree
Hide file tree
Showing 8 changed files with 170 additions and 132 deletions.
49 changes: 35 additions & 14 deletions cesm/driver/esm.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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))
Expand All @@ -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)
Expand Down
40 changes: 30 additions & 10 deletions mediator/med.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)'
!-----------------------------------------------------------
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
6 changes: 3 additions & 3 deletions mediator/med_fraction_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
8 changes: 6 additions & 2 deletions mediator/med_internalstate_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand Down
41 changes: 20 additions & 21 deletions mediator/med_io_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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(:,:)
Expand All @@ -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) '
!-------------------------------------------------------------------------------
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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))
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)

Expand Down
Loading

0 comments on commit 4e19850

Please sign in to comment.