sfc_climo_gen 1.14.0
Loading...
Searching...
No Matches
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
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:15
model_grid
This module defines the model grid.
Definition model_grid.F90:10
model_grid::grid_mdl
type(esmf_grid), public grid_mdl
ESMF grid object for the model grid.
Definition model_grid.F90:28
model_grid::i_mdl
integer, public i_mdl
i dimension of model tile.
Definition model_grid.F90:20
model_grid::mask_field_mdl
type(esmf_field), public mask_field_mdl
ESMF field object that holds the model land mask.
Definition model_grid.F90:37
model_grid::num_tiles
integer, public num_tiles
Total number of model grid tiles.
Definition model_grid.F90:23
model_grid::latitude_field_mdl
type(esmf_field), public latitude_field_mdl
ESMF field object that holds the model grid latitude.
Definition model_grid.F90:41
model_grid::j_mdl
integer, public j_mdl
j dimension of model tile.
Definition model_grid.F90:21
model_grid::land_frac_field_mdl
type(esmf_field), public land_frac_field_mdl
ESMF field object that holds the model land fraction.
Definition model_grid.F90:31
model_grid::longitude_field_mdl
type(esmf_field), public longitude_field_mdl
ESMF field object that holds the model grid longitude.
Definition model_grid.F90:43
output_frac_cats
Output categorical data such as vegetation type.
Definition output_frac_cats.F90:10
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:33
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
source_grid::field_names
character(len=50), dimension(:), allocatable, public field_names
Names of fields to be processed.
Definition source_grid.F90:21
source_grid::grid_src
type(esmf_grid), public grid_src
ESMF grid object for the source grid.
Definition source_grid.F90:34
source_grid::j_src
integer, public j_src
j dimension of the source grid.
Definition source_grid.F90:25
source_grid::i_src
integer, public i_src
i dimension of the source grid.
Definition source_grid.F90:24
utils
Utilities.
Definition utils.f90:8
utils::netcdf_err
subroutine, public netcdf_err(err, string)
Handle netCDF error codes.
Definition utils.f90:23
utils::error_handler
subroutine, public error_handler(string, rc)
Handle errors.
Definition utils.f90:49
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:25