Skip to content

Commit cc929a8

Browse files
GeorgeGayno-NOAAGeorgeGayno-NOAA
committed
Update logic to work with fractional land/non-land grids.
Fixes ufs-community#709.
1 parent bb91f34 commit cc929a8

File tree

2 files changed

+81
-10
lines changed

2 files changed

+81
-10
lines changed

‎sorc/sfc_climo_gen.fd/interp2.F90

+39-4
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,8 @@ subroutine interp2(localpet, input_file)
1515
use netcdf
1616
use model_grid, only : grid_mdl, i_mdl, j_mdl, &
1717
num_tiles, latitude_field_mdl, &
18-
longitude_field_mdl, mask_field_mdl
18+
longitude_field_mdl, mask_field_mdl, &
19+
land_frac_field_mdl
1920
use source_grid
2021
use utils
2122
use mpi
@@ -43,6 +44,7 @@ subroutine interp2(localpet, input_file)
4344
real(esmf_kind_r4), allocatable :: lat_mdl_one_tile(:,:)
4445
real(esmf_kind_r4), allocatable :: sum_mdl_one_tile(:,:)
4546
real(esmf_kind_r4), allocatable :: lon_mdl_one_tile(:,:)
47+
real(esmf_kind_r4), allocatable :: land_frac_mdl_one_tile(:,:)
4648

4749
type(esmf_regridmethod_flag) :: method
4850
type(esmf_field) :: data_field_src
@@ -99,6 +101,7 @@ subroutine interp2(localpet, input_file)
99101
allocate(data_mdl_one_tile(i_mdl,j_mdl,num_categories))
100102
allocate(dom_cat_mdl_one_tile(i_mdl,j_mdl))
101103
allocate(mask_mdl_one_tile(i_mdl,j_mdl))
104+
allocate(land_frac_mdl_one_tile(i_mdl,j_mdl))
102105
allocate(lat_mdl_one_tile(i_mdl,j_mdl))
103106
allocate(sum_mdl_one_tile(i_mdl,j_mdl))
104107
allocate(lon_mdl_one_tile(i_mdl,j_mdl))
@@ -107,6 +110,7 @@ subroutine interp2(localpet, input_file)
107110
allocate(data_mdl_one_tile(0,0,0))
108111
allocate(dom_cat_mdl_one_tile(0,0))
109112
allocate(mask_mdl_one_tile(0,0))
113+
allocate(land_frac_mdl_one_tile(0,0))
110114
allocate(lat_mdl_one_tile(0,0))
111115
allocate(sum_mdl_one_tile(0,0))
112116
allocate(lon_mdl_one_tile(0,0))
@@ -183,25 +187,56 @@ subroutine interp2(localpet, input_file)
183187
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
184188
call error_handler("IN FieldGather.", rc)
185189

190+
print*,"- CALL FieldGather FOR MODEL GRID LAND FRACTION."
191+
call ESMF_FieldGather(land_frac_field_mdl, land_frac_mdl_one_tile, rootPet=0, tile=tile, rc=rc)
192+
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
193+
call error_handler("IN FieldGather.", rc)
194+
186195
if (localpet == 0) then
187196
print*,'- CALL SEARCH FOR TILE ',tile
188-
sum_mdl_one_tile = sum(data_mdl_one_tile, dim=3) ! use unused variable to now.
197+
198+
! Where sum is zero, the regridding did not find any input data for the model point
199+
! (ex. and isolated island). Call the search routine at these points.
200+
sum_mdl_one_tile = sum(data_mdl_one_tile, dim=3)
189201
do j = 1, j_mdl
190202
do i = 1, i_mdl
191203
if (mask_mdl_one_tile(i,j) == 1 .and. sum_mdl_one_tile(i,j) == 0.0) then
192-
data_mdl_one_tile(i,j,:) = -9999.9
204+
data_mdl_one_tile(i,j,:) = -9999.9 ! flag to tell search routine to search.
193205
endif
194206
enddo
195207
enddo
196208
call search2 (data_mdl_one_tile, mask_mdl_one_tile, i_mdl, j_mdl, num_categories, tile, field_names(1))
197209
print*,'after regrid ',data_mdl_one_tile(i_mdl/2,j_mdl/2,:)
210+
211+
! These points are all non-land. Set to 100% of the water category.
198212
do j = 1, j_mdl
199213
do i = 1, i_mdl
200214
if (mask_mdl_one_tile(i,j) == 0) then
201215
data_mdl_one_tile(i,j,water_category) = 1.0
202216
endif
203217
enddo
204218
enddo
219+
220+
! For fractional grids, need to rescale the category percentages by the
221+
! fraction of land in the model grid cell.
222+
223+
! When running with fractional grids, the land_frac_mdl_one_tile array will
224+
! contain a fraction between 0 and 1. When not running with fractional
225+
! grids, this array will contain negative fill values.
226+
227+
if (maxval(land_frac_mdl_one_tile) > 0.0) then
228+
print*,'before rescale ',land_frac_mdl_one_tile(658,95),data_mdl_one_tile(658,95,:)
229+
do j = 1, j_mdl
230+
do i = 1, i_mdl
231+
if (mask_mdl_one_tile(i,j) == 1) then
232+
data_mdl_one_tile(i,j,:) = data_mdl_one_tile(i,j,:) * land_frac_mdl_one_tile(i,j)
233+
data_mdl_one_tile(i,j,water_category) = 1.0 - land_frac_mdl_one_tile(i,j)
234+
endif
235+
enddo
236+
enddo
237+
print*,'after rescale ',land_frac_mdl_one_tile(658,95),data_mdl_one_tile(658,95,:)
238+
endif
239+
! under fractional grids, how do we define dominate category?
205240
dom_cat_mdl_one_tile = 0.0
206241
dom_cat_mdl_one_tile = maxloc(data_mdl_one_tile,dim=3)
207242
call output2 (data_mdl_one_tile, dom_cat_mdl_one_tile, lat_mdl_one_tile, lon_mdl_one_tile, i_mdl, j_mdl, num_categories, tile)
@@ -212,7 +247,7 @@ subroutine interp2(localpet, input_file)
212247
status=nf90_close(ncid)
213248

214249
deallocate(data_mdl_one_tile, dom_cat_mdl_one_tile, mask_mdl_one_tile, data_src_global2)
215-
deallocate(lat_mdl_one_tile, lon_mdl_one_tile, sum_mdl_one_tile)
250+
deallocate(lat_mdl_one_tile, lon_mdl_one_tile, sum_mdl_one_tile, land_frac_mdl_one_tile)
216251

217252
print*,"- CALL FieldRegridRelease."
218253
call ESMF_FieldRegridRelease(routehandle=regrid_data, rc=rc)

‎sorc/sfc_climo_gen.fd/model_grid.F90

+42-6
Original file line numberDiff line numberDiff line change
@@ -28,8 +28,16 @@ module model_grid
2828
type(esmf_grid), public :: grid_mdl !< ESMF grid object for the model grid.
2929
type(esmf_field), public :: data_field_mdl !< ESMF field object that holds the
3030
!! data interpolated to model grid.
31+
type(esmf_field), public :: land_frac_field_mdl !< ESMF field object that holds the
32+
!! model land fraction. When running
33+
!! with fractional grids, will be
34+
!! between zero and one. For non-
35+
!! fractional grids, will contain a
36+
!! fill value.
3137
type(esmf_field), public :: mask_field_mdl !< ESMF field object that holds the
32-
!! model land mask.
38+
!! model land mask. Equal to '1' if
39+
!! point is partial or all land. Equal
40+
!! to zero is point is all non-land.
3341
type(esmf_field), public :: latitude_field_mdl !< ESMF field object that holds the
3442
!! model grid latitude.
3543
type(esmf_field), public :: longitude_field_mdl !< ESMF field object that holds the
@@ -75,6 +83,7 @@ subroutine define_model_grid(localpet, npets)
7583

7684
real(esmf_kind_r4), allocatable :: latitude_one_tile(:,:)
7785
real(esmf_kind_r4), allocatable :: longitude_one_tile(:,:)
86+
real(esmf_kind_r4), allocatable :: land_frac_one_tile(:,:)
7887

7988
!-----------------------------------------------------------------------
8089
! Get the number of tiles from the mosaic file.
@@ -214,20 +223,31 @@ subroutine define_model_grid(localpet, npets)
214223
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
215224
call error_handler("IN FieldGet", rc)
216225

226+
print*,"- CALL FieldCreate FOR MODEL GRID LAND FRACTION."
227+
land_frac_field_mdl = ESMF_FieldCreate(grid_mdl, &
228+
typekind=ESMF_TYPEKIND_R4, &
229+
staggerloc=ESMF_STAGGERLOC_CENTER, &
230+
name="model grid land fraction", &
231+
rc=rc)
232+
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
233+
call error_handler("IN FieldCreate", rc)
234+
217235
if (localpet == 0) then
218236
allocate(mask_mdl_one_tile(i_mdl,j_mdl))
237+
allocate(land_frac_one_tile(i_mdl,j_mdl))
219238
allocate(latitude_one_tile(i_mdl,j_mdl))
220239
allocate(longitude_one_tile(i_mdl,j_mdl))
221240
else
222241
allocate(mask_mdl_one_tile(0,0))
242+
allocate(land_frac_one_tile(0,0))
223243
allocate(latitude_one_tile(0,0))
224244
allocate(longitude_one_tile(0,0))
225245
endif
226246

227247
do tile = 1, num_tiles
228248
if (localpet == 0) then
229249
the_file = trim(orog_dir_mdl) // trim(orog_files_mdl(tile))
230-
call get_model_info(trim(the_file), mask_mdl_one_tile, &
250+
call get_model_info(trim(the_file), mask_mdl_one_tile, land_frac_one_tile, &
231251
latitude_one_tile, longitude_one_tile, i_mdl, j_mdl)
232252
endif
233253

@@ -236,6 +256,11 @@ subroutine define_model_grid(localpet, npets)
236256
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
237257
call error_handler("IN FieldScatter", rc)
238258

259+
print*,"- CALL FieldScatter FOR MODEL GRID LAND FRACTION. TILE IS: ", tile
260+
call ESMF_FieldScatter(land_frac_field_mdl, land_frac_one_tile, rootpet=0, tile=tile, rc=rc)
261+
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
262+
call error_handler("IN FieldScatter", rc)
263+
239264
print*,"- CALL FieldScatter FOR MODEL LATITUDE. TILE IS: ", tile
240265
call ESMF_FieldScatter(latitude_field_mdl, latitude_one_tile, rootpet=0, tile=tile, rc=rc)
241266
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
@@ -282,7 +307,7 @@ end subroutine define_model_grid
282307
!! @param[in] idim i dimension of the model tile
283308
!! @param[in] jdim j dimension of the model tile
284309
!! @author George Gayno @date 2018
285-
subroutine get_model_info(orog_file, mask, lat2d, lon2d, idim, jdim)
310+
subroutine get_model_info(orog_file, mask, land_frac, lat2d, lon2d, idim, jdim)
286311

287312
use esmf
288313
use netcdf
@@ -297,6 +322,7 @@ subroutine get_model_info(orog_file, mask, lat2d, lon2d, idim, jdim)
297322

298323
real(esmf_kind_r4), intent(out) :: lat2d(idim,jdim)
299324
real(esmf_kind_r4), intent(out) :: lon2d(idim,jdim)
325+
real(esmf_kind_r4), intent(out) :: land_frac(idim,jdim)
300326

301327
integer :: error, lat, lon, i, j
302328
integer :: ncid, id_dim, id_var
@@ -330,11 +356,17 @@ subroutine get_model_info(orog_file, mask, lat2d, lon2d, idim, jdim)
330356
allocate(dummy(idim,jdim))
331357

332358
!-----------------------------------------------------------------------
333-
! If the lake maker was used, there will be a 'lake_frac' record.
359+
! If the lake maker was used, we are running with a fractional
360+
! land/non-land grid and there will be a 'lake_frac' record.
334361
! In that case, land/non-land is determined by 'land_frac'.
335362
!
336363
! If the lake maker was not used, use 'slmsk', which is defined
337364
! as the nint(land_frac).
365+
!
366+
! In summary, if 'mask' is one, the point is all land or
367+
! partial land and surface data will be mapped to it. Otherwise,
368+
! when 'mask' is zero, then the point is all non-land and
369+
! surface data will not be mapped to it.
338370
!-----------------------------------------------------------------------
339371

340372
error=nf90_inq_varid(ncid, 'lake_frac', id_var)
@@ -345,16 +377,17 @@ subroutine get_model_info(orog_file, mask, lat2d, lon2d, idim, jdim)
345377
error=nf90_get_var(ncid, id_var, dummy)
346378
call netcdf_err(error, "READING SLMSK")
347379
mask = nint(dummy)
380+
land_frac = -999.
348381
else
349382
print*,"- READ LAND FRACTION"
350383
error=nf90_inq_varid(ncid, 'land_frac', id_var)
351384
call netcdf_err(error, "READING LAND_FRAC ID")
352-
error=nf90_get_var(ncid, id_var, dummy)
385+
error=nf90_get_var(ncid, id_var, land_frac)
353386
call netcdf_err(error, "READING LAND_FRAC")
354387
mask = 0
355388
do j = 1, lat
356389
do i = 1, lon
357-
if (dummy(i,j) > 0.0) then
390+
if (land_frac(i,j) > 0.0) then
358391
mask(i,j) = 1
359392
endif
360393
enddo
@@ -398,6 +431,9 @@ subroutine model_grid_cleanup
398431
print*,"- CALL FieldDestroy FOR MODEL GRID LAND MASK."
399432
call ESMF_FieldDestroy(mask_field_mdl,rc=rc)
400433

434+
print*,"- CALL FieldDestroy FOR MODEL GRID LAND MASK."
435+
call ESMF_FieldDestroy(land_frac_field_mdl,rc=rc)
436+
401437
print*,"- CALL FieldDestroy FOR MODEL GRID DATA FIELD."
402438
call ESMF_FieldDestroy(data_field_mdl,rc=rc)
403439

0 commit comments

Comments
 (0)