sfc_climo_gen  1.13.0
 All Data Structures Files Functions Variables Pages
interp_frac_cats.F90
Go to the documentation of this file.
1 
6 
14  subroutine interp_frac_cats(localpet, input_file)
15 
16  use esmf
17  use netcdf
18  use model_grid, only : grid_mdl, i_mdl, j_mdl, &
19  num_tiles, latitude_field_mdl, &
20  longitude_field_mdl, mask_field_mdl, &
21  land_frac_field_mdl
22  use source_grid
23  use output_frac_cats, only : output_driver
24  use utils
25  use mpi
26 
27  implicit none
28 
29  character(len=*), intent(in) :: input_file
30 
31  integer, intent(in) :: localpet
32 
33  integer :: i, j, tile, cat, ncid, status, rc
34  integer :: varid, water_category, max_cat
35  integer :: isrctermprocessing
36  integer :: category, num_categories
37 
38  integer(esmf_kind_i4), allocatable :: mask_mdl_one_tile(:,:)
39  integer(esmf_kind_i4), pointer :: unmapped_ptr(:)
40 
41  real(esmf_kind_r4), allocatable :: data_src_global(:,:)
42  real(esmf_kind_r4), allocatable :: data_src_global2(:,:,:)
43  real(esmf_kind_r4), allocatable :: data_mdl_one_tile(:,:,:)
44  real(esmf_kind_r4), allocatable :: dom_cat_mdl_one_tile(:,:)
45  real(esmf_kind_r4), allocatable :: lat_mdl_one_tile(:,:)
46  real(esmf_kind_r4), allocatable :: sum_mdl_one_tile(:,:)
47  real(esmf_kind_r4), allocatable :: lon_mdl_one_tile(:,:)
48  real(esmf_kind_r4), allocatable :: land_frac_mdl_one_tile(:,:)
49  real(esmf_kind_r4) :: max_frac
50 
51  type(esmf_regridmethod_flag) :: method
52  type(esmf_field) :: data_field_src
53  type(esmf_field) :: data_field_mdl
54  type(esmf_routehandle) :: regrid_data
55  type(esmf_polemethod_flag) :: pole
56 
57  if (localpet == 0) then
58  allocate(data_src_global(i_src,j_src))
59  else
60  allocate(data_src_global(0,0))
61  endif
62 
63  if (localpet == 0) then
64  print*,'- OPEN SOURCE FILE ', trim(input_file)
65  status = nf90_open(input_file, nf90_nowrite, ncid)
66  call netcdf_err(status, "IN ROUTINE INTERP OPENING SOURCE FILE")
67  status = nf90_inq_varid(ncid, field_names(1), varid)
68  call netcdf_err(status, "IN ROUTINE INTERP READING FIELD ID")
69  status = nf90_get_var(ncid, varid, data_src_global, start=(/1,1,1/), count=(/i_src,j_src,1/))
70  call netcdf_err(status, "IN ROUTINE INTERP READING FIELD")
71  print*,'number of cats ',maxval(data_src_global)
72  num_categories = nint(maxval(data_src_global))
73  status = nf90_get_att(ncid, varid, 'water_category', water_category)
74  call netcdf_err(status, "IN ROUTINE INTERP READING water_category")
75  print*,'water cat ',water_category
76  endif
77 
78  call mpi_bcast(num_categories,1,mpi_integer,0,mpi_comm_world,status)
79 
80  print*,"- CALL FieldCreate FOR SOURCE GRID DATA."
81  data_field_src = esmf_fieldcreate(grid_src, &
82  typekind=esmf_typekind_r4, &
83  indexflag=esmf_index_global, &
84  staggerloc=esmf_staggerloc_center, &
85  ungriddedlbound=(/1/), &
86  ungriddedubound=(/num_categories/), &
87  name="source data", &
88  rc=rc)
89  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
90  call error_handler("IN FieldCreate", rc)
91 
92  print*,"- CALL FieldCreate FOR model GRID veg DATA."
93  data_field_mdl = esmf_fieldcreate(grid_mdl, &
94  typekind=esmf_typekind_r4, &
95  indexflag=esmf_index_global, &
96  staggerloc=esmf_staggerloc_center, &
97  ungriddedlbound=(/1/), &
98  ungriddedubound=(/num_categories/), &
99  name="mdl data", &
100  rc=rc)
101  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
102  call error_handler("IN FieldCreate", rc)
103 
104  if (localpet == 0) then
105  allocate(data_src_global2(i_src,j_src,num_categories))
106  allocate(data_mdl_one_tile(i_mdl,j_mdl,num_categories))
107  allocate(dom_cat_mdl_one_tile(i_mdl,j_mdl))
108  allocate(mask_mdl_one_tile(i_mdl,j_mdl))
109  allocate(land_frac_mdl_one_tile(i_mdl,j_mdl))
110  allocate(lat_mdl_one_tile(i_mdl,j_mdl))
111  allocate(sum_mdl_one_tile(i_mdl,j_mdl))
112  allocate(lon_mdl_one_tile(i_mdl,j_mdl))
113  else
114  allocate(data_src_global2(0,0,0))
115  allocate(data_mdl_one_tile(0,0,0))
116  allocate(dom_cat_mdl_one_tile(0,0))
117  allocate(mask_mdl_one_tile(0,0))
118  allocate(land_frac_mdl_one_tile(0,0))
119  allocate(lat_mdl_one_tile(0,0))
120  allocate(sum_mdl_one_tile(0,0))
121  allocate(lon_mdl_one_tile(0,0))
122  endif
123 
124  if (localpet == 0) then
125  data_src_global2 = 0.0
126  do j = 1, j_src
127  do i = 1, i_src
128  category = nint(data_src_global(i,j))
129  if (category < 1) cycle
130  data_src_global2(i,j,category) = 1.0
131  enddo
132  enddo
133  endif
134 
135  deallocate(data_src_global)
136 
137  print*,"- CALL FieldScatter FOR SOURCE GRID DATA."
138  call esmf_fieldscatter(data_field_src, data_src_global2, rootpet=0, rc=rc)
139  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
140  call error_handler("IN FieldScatter.", rc)
141 
142  deallocate(data_src_global2)
143 
144  isrctermprocessing = 1
145 
146  method = esmf_regridmethod_conserve
147  pole = esmf_polemethod_none
148 
149  print*,"- CALL FieldRegridStore."
150  nullify(unmapped_ptr)
151  call esmf_fieldregridstore(data_field_src, &
152  data_field_mdl, &
153  srcmaskvalues=(/0/), &
154  dstmaskvalues=(/0/), &
155  polemethod=pole, &
156  unmappedaction=esmf_unmappedaction_ignore, &
157  normtype=esmf_normtype_fracarea, &
158  srctermprocessing=isrctermprocessing, &
159  routehandle=regrid_data, &
160  regridmethod=method, &
161  unmappeddstlist=unmapped_ptr, &
162  rc=rc)
163  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
164  call error_handler("IN FieldRegridStore.", rc)
165 
166  print*,"- CALL Field_Regrid."
167  call esmf_fieldregrid(data_field_src, &
168  data_field_mdl, &
169  routehandle=regrid_data, &
170  termorderflag=esmf_termorder_srcseq, &
171  rc=rc)
172  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
173  call error_handler("IN FieldRegrid.", rc)
174 
175  print*,"- CALL FieldRegridRelease."
176  call esmf_fieldregridrelease(routehandle=regrid_data, rc=rc)
177 
178  print*,"- CALL FieldDestroy."
179  call esmf_fielddestroy(data_field_src, rc=rc)
180 
181  output_loop : do tile = 1, num_tiles
182 
183  print*,"- CALL FieldGather FOR MODEL LATITUDE."
184  call esmf_fieldgather(latitude_field_mdl, lat_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
185  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
186  call error_handler("IN FieldGather.", rc)
187 
188  print*,"- CALL FieldGather FOR MODEL LONGITUDE."
189  call esmf_fieldgather(longitude_field_mdl, lon_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
190  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
191  call error_handler("IN FieldGather.", rc)
192 
193  print*,"- CALL FieldGather FOR MODEL GRID DATA."
194  call esmf_fieldgather(data_field_mdl, data_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
195  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
196  call error_handler("IN FieldGather.", rc)
197 
198  print*,"- CALL FieldGather FOR MODEL GRID LAND MASK."
199  call esmf_fieldgather(mask_field_mdl, mask_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
200  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
201  call error_handler("IN FieldGather.", rc)
202 
203  print*,"- CALL FieldGather FOR MODEL GRID LAND FRACTION."
204  call esmf_fieldgather(land_frac_field_mdl, land_frac_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
205  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
206  call error_handler("IN FieldGather.", rc)
207 
208  if (localpet == 0) then
209  print*,'- CALL SEARCH FOR TILE ',tile
210 
211 ! Where sum is zero, the regridding did not find any input data for the model point
212 ! (ex. and isolated island). Call the search routine at these points.
213  sum_mdl_one_tile = sum(data_mdl_one_tile, dim=3)
214  do j = 1, j_mdl
215  do i = 1, i_mdl
216  if (mask_mdl_one_tile(i,j) == 1 .and. sum_mdl_one_tile(i,j) == 0.0) then
217  data_mdl_one_tile(i,j,:) = -9999.9 ! flag to tell search routine to search.
218  endif
219  enddo
220  enddo
221  call search_frac_cats(data_mdl_one_tile, mask_mdl_one_tile, i_mdl, j_mdl, num_categories, tile, field_names(1))
222  print*,'after regrid ',data_mdl_one_tile(i_mdl/2,j_mdl/2,:)
223 
224 ! These points are all non-land. Set to 100% of the water category.
225  do j = 1, j_mdl
226  do i = 1, i_mdl
227  if (mask_mdl_one_tile(i,j) == 0) then
228  data_mdl_one_tile(i,j,water_category) = 1.0
229  endif
230  enddo
231  enddo
232 
233 ! For fractional grids, need to rescale the category percentages by the
234 ! fraction of land in the model grid cell.
235 
236 ! When running with fractional grids, the land_frac_mdl_one_tile array will
237 ! contain a fraction between 0 and 1. When not running with fractional
238 ! grids, this array will contain negative fill values.
239 
240  if (maxval(land_frac_mdl_one_tile) > 0.0) then
241  print*,'before rescale ',tile,land_frac_mdl_one_tile(42,82),data_mdl_one_tile(42,82,:)
242  do j = 1, j_mdl
243  do i = 1, i_mdl
244  if (mask_mdl_one_tile(i,j) == 1) then
245  data_mdl_one_tile(i,j,:) = data_mdl_one_tile(i,j,:) * land_frac_mdl_one_tile(i,j)
246  data_mdl_one_tile(i,j,water_category) = 1.0 - land_frac_mdl_one_tile(i,j)
247  max_frac = 0.0
248  max_cat = -9
249  do cat = 1, num_categories
250  if (cat == water_category) cycle
251  if(data_mdl_one_tile(i,j,cat) > max_frac) then
252  max_frac = data_mdl_one_tile(i,j,cat)
253  max_cat = cat
254  endif
255  enddo
256  dom_cat_mdl_one_tile(i,j) = max_cat
257  else
258  dom_cat_mdl_one_tile(i,j) = water_category
259  endif
260  enddo
261  enddo
262  else ! non-fractonal grids.
263  dom_cat_mdl_one_tile = 0.0
264  dom_cat_mdl_one_tile = maxloc(data_mdl_one_tile,dim=3)
265  endif
266  call output_driver(data_mdl_one_tile, dom_cat_mdl_one_tile, lat_mdl_one_tile, lon_mdl_one_tile, i_mdl, j_mdl, num_categories, tile)
267  endif
268 
269  enddo output_loop
270 
271  status=nf90_close(ncid)
272 
273  deallocate(data_mdl_one_tile, dom_cat_mdl_one_tile, mask_mdl_one_tile)
274  deallocate(lat_mdl_one_tile, lon_mdl_one_tile, sum_mdl_one_tile, land_frac_mdl_one_tile)
275 
276  print*,"- CALL FieldDestroy."
277  call esmf_fielddestroy(data_field_mdl, rc=rc)
278 
279  call mpi_barrier(mpi_comm_world, rc)
280 
281  end subroutine interp_frac_cats
search_frac_cats
subroutine search_frac_cats(field, mask, idim, jdim, num_categories, tile, field_name)
Replace undefined values on the model grid with valid values at a nearby neighbor.
Definition: search_frac_cats.f90:24
output_frac_cats::output_driver
subroutine, public output_driver(data_one_tile, dom_cat_one_tile, lat_one_tile, lon_one_tile, i_mdl, j_mdl, num_categories, tile)
Driver routine to output model categorical data.
Definition: output_frac_cats.F90:31
model_grid
This module defines the model grid.
Definition: model_grid.F90:10
interp_frac_cats
subroutine interp_frac_cats(localpet, input_file)
Read the input source data and interpolate it to the model grid.
Definition: interp_frac_cats.F90:14
utils::netcdf_err
subroutine, public netcdf_err(err, string)
Handle netCDF error codes.
Definition: utils.f90:22
utils::error_handler
subroutine, public error_handler(string, rc)
Handle errors.
Definition: utils.f90:48
source_grid
Read grid specs, date information and land/sea mask for the source data that will be interpolated to ...
Definition: source_grid.F90:12
output_frac_cats
Output categorical data such as vegetation type.
Definition: output_frac_cats.F90:10
utils
Utilities.
Definition: utils.f90:8