sfc_climo_gen 1.14.0
Loading...
Searching...
No Matches
interp.F90
Go to the documentation of this file.
1
5
13 subroutine interp(localpet, method, input_file)
14
15 use esmf
16 use netcdf
17 use model_grid
18 use program_setup, only : fract_vegsoil_type
19 use source_grid
20 use utils
21 use mpi
22
23 implicit none
24
25 character(len=*), intent(in) :: input_file
26
27 integer :: rc, localpet
28 integer :: i, j, ij, tile, n, ncid, status
29 integer :: l(1), u(1), t
30 integer :: clb_mdl(2), cub_mdl(2)
31 integer :: varid, record
32 integer :: tile_num, pt_loc_this_tile
33 integer :: isrctermprocessing
34
35 integer(esmf_kind_i4), allocatable :: mask_mdl_one_tile(:,:)
36 integer(esmf_kind_i4), pointer :: unmapped_ptr(:)
37
38 real(esmf_kind_r4), pointer :: data_mdl_ptr(:,:)
39 real(esmf_kind_r4), pointer :: data_src_ptr(:,:)
40 real(esmf_kind_r4), allocatable :: data_src_global(:,:)
41 real(esmf_kind_r4), allocatable :: data_mdl_one_tile(:,:)
42 real(esmf_kind_r4), allocatable :: vegt_mdl_one_tile(:,:)
43 real(esmf_kind_r4), allocatable :: lat_mdl_one_tile(:,:)
44 real(esmf_kind_r4), allocatable :: lon_mdl_one_tile(:,:)
45
46 type(esmf_regridmethod_flag),intent(in) :: method
47 type(esmf_field) :: data_field_src
48 type(esmf_routehandle) :: regrid_data
49 type(esmf_polemethod_flag) :: pole
50
51 print*,"- CALL FieldCreate FOR SOURCE GRID DATA."
52 data_field_src = esmf_fieldcreate(grid_src, &
53 typekind=esmf_typekind_r4, &
54 indexflag=esmf_index_global, &
55 staggerloc=esmf_staggerloc_center, &
56 name="source data", &
57 rc=rc)
58 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
59 call error_handler("IN FieldCreate", rc)
60
61 print*,"- CALL FieldGet FOR SOURCE GRID DATA."
62 nullify(data_src_ptr)
63 call esmf_fieldget(data_field_src, &
64 farrayptr=data_src_ptr, &
65 rc=rc)
66 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
67 call error_handler("IN FieldGet", rc)
68
69 print*,"- CALL FieldGet FOR MODEL GRID DATA."
70 nullify(data_mdl_ptr)
71 call esmf_fieldget(data_field_mdl, &
72 farrayptr=data_mdl_ptr, &
73 computationallbound=clb_mdl, &
74 computationalubound=cub_mdl, &
75 rc=rc)
76 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
77 call error_handler("IN FieldGet", rc)
78
79 if (localpet == 0) then
80 allocate(data_src_global(i_src,j_src))
81 else
82 allocate(data_src_global(0,0))
83 endif
84
85 print*,'- OPEN SOURCE FILE ', trim(input_file)
86 status = nf90_open(input_file, nf90_nowrite, ncid)
87 call netcdf_err(status, "IN ROUTINE INTERP OPENING SOURCE FILE")
88
89 if (localpet == 0) then
90 allocate(data_mdl_one_tile(i_mdl,j_mdl))
91 allocate(mask_mdl_one_tile(i_mdl,j_mdl))
92 allocate(lat_mdl_one_tile(i_mdl,j_mdl))
93 allocate(lon_mdl_one_tile(i_mdl,j_mdl))
94 else
95 allocate(data_mdl_one_tile(0,0))
96 allocate(mask_mdl_one_tile(0,0))
97 allocate(lat_mdl_one_tile(0,0))
98 allocate(lon_mdl_one_tile(0,0))
99 endif
100
101 record = 0
102
103 time_loop : do t = 1, num_time_recs ! loop over each time period
104
105 field_loop : do n = 1, num_fields ! loop over each surface field.
106
107 record = record + 1
108
109 if (localpet == 0) then
110 status = nf90_inq_varid(ncid, field_names(n), varid)
111 call netcdf_err(status, "IN ROUTINE INTERP READING FIELD ID")
112 status = nf90_get_var(ncid, varid, data_src_global, start=(/1,1,t/), count=(/i_src,j_src,1/))
113 call netcdf_err(status, "IN ROUTINE INTERP READING FIELD")
114 endif
115
116 print*,"- CALL FieldScatter FOR SOURCE GRID DATA."
117 call esmf_fieldscatter(data_field_src, data_src_global, rootpet=0, rc=rc)
118 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
119 call error_handler("IN FieldScatter.", rc)
120
121 isrctermprocessing = 1
122
123 if (record == 1) then
124
125 if (method == esmf_regridmethod_bilinear) then
126 pole = esmf_polemethod_allavg
127 else
128 pole = esmf_polemethod_none
129 endif
130
131 print*,"- CALL FieldRegridStore."
132 nullify(unmapped_ptr)
133 call esmf_fieldregridstore(data_field_src, &
134 data_field_mdl, &
135 srcmaskvalues=(/0/), &
136 dstmaskvalues=(/0/), &
137 polemethod=pole, &
138 unmappedaction=esmf_unmappedaction_ignore, &
139 normtype=esmf_normtype_fracarea, &
140 srctermprocessing=isrctermprocessing, &
141 routehandle=regrid_data, &
142 regridmethod=method, &
143 unmappeddstlist=unmapped_ptr, &
144 rc=rc)
145 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
146 call error_handler("IN FieldRegridStore.", rc)
147
148 endif
149
150 print*,"- CALL Field_Regrid."
151 call esmf_fieldregrid(data_field_src, &
152 data_field_mdl, &
153 routehandle=regrid_data, &
154 termorderflag=esmf_termorder_srcseq, &
155 rc=rc)
156 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
157 call error_handler("IN FieldRegrid.", rc)
158
159!-----------------------------------------------------------------------
160! Unmapped points are stored in "unmapped_ptr". The pointer contains
161! "ij" global indices as follows:
162!
163! tile 1: 1 thru (itile*jtile)
164! tile n: (n-1)*(itile*jtile) thru n*(itile*jtile)
165!
166! This "ij" index is converted to the tile number and i/j index of the
167! field object. This logic assumes the model grid object was
168! created using "GLOBAL" indices.
169!
170! Unmapped data points are given the flag value of -9999.9, which
171! will be replaced in routine "search".
172!-----------------------------------------------------------------------
173
174 l = lbound(unmapped_ptr)
175 u = ubound(unmapped_ptr)
176 do ij = l(1), u(1)
177
178 tile_num = ((unmapped_ptr(ij)-1) / (i_mdl*j_mdl)) ! tile number minus 1
179 pt_loc_this_tile = unmapped_ptr(ij) - (tile_num * i_mdl * j_mdl)
180 ! "ij" location of point within tile.
181
182 j = (pt_loc_this_tile - 1) / i_mdl + 1
183 i = mod(pt_loc_this_tile, i_mdl)
184 if (i==0) i = i_mdl
185 data_mdl_ptr(i,j) = -9999.9
186
187 enddo
188
189! Adjust some fields at permanent land ice points. These points are identified by the
190! 'permanent ice' vegetation type category.
191!
192! When outputting the fraction of each vegetation type, land ice points are
193! not defined. So don't do this adjustment.
194
195 if (.not. fract_vegsoil_type) then
196 select case (trim(field_names(n)))
197 case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type','soil_color')
198 if (localpet == 0) then
199 allocate(vegt_mdl_one_tile(i_mdl,j_mdl))
200 else
201 allocate(vegt_mdl_one_tile(0,0))
202 endif
203 end select
204 endif
205
206 output_loop : do tile = 1, num_tiles
207
208 print*,"- CALL FieldGather FOR MODEL LATITUDE."
209 call esmf_fieldgather(latitude_field_mdl, lat_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
210 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
211 call error_handler("IN FieldGather.", rc)
212
213 print*,"- CALL FieldGather FOR MODEL LONGITUDE."
214 call esmf_fieldgather(longitude_field_mdl, lon_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
215 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
216 call error_handler("IN FieldGather.", rc)
217
218 print*,"- CALL FieldGather FOR MODEL GRID DATA."
219 call esmf_fieldgather(data_field_mdl, data_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
220 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
221 call error_handler("IN FieldGather.", rc)
222
223 print*,"- CALL FieldGather FOR MODEL GRID LAND MASK."
224 call esmf_fieldgather(mask_field_mdl, mask_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
225 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
226 call error_handler("IN FieldGather.", rc)
227
228 if (.not. fract_vegsoil_type) then
229 select case (trim(field_names(n)))
230 case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type','soil_color')
231 print*,"- CALL FieldGather FOR MODEL GRID VEG TYPE."
232 call esmf_fieldgather(vegt_field_mdl, vegt_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
233 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
234 call error_handler("IN FieldGather.", rc)
235 end select
236 endif
237
238 if (localpet == 0) then
239 print*,'- CALL SEARCH FOR TILE ',tile
240 call search (data_mdl_one_tile, mask_mdl_one_tile, i_mdl, j_mdl, tile, field_names(n))
241 if (.not. fract_vegsoil_type) then
242 select case (field_names(n))
243 case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type','soil_color')
244 call adjust_for_landice (data_mdl_one_tile, vegt_mdl_one_tile, i_mdl, j_mdl, field_names(n))
245 end select
246 endif
247 where(mask_mdl_one_tile == 0) data_mdl_one_tile = missing
248 call output (data_mdl_one_tile, lat_mdl_one_tile, lon_mdl_one_tile, i_mdl, j_mdl, tile, record, t, n)
249 endif
250
251 if (.not. fract_vegsoil_type) then
252 if (field_names(n) == 'vegetation_type') then
253 print*,"- CALL FieldScatter FOR MODEL GRID VEGETATION TYPE."
254 call esmf_fieldscatter(vegt_field_mdl, data_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
255 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
256 call error_handler("IN FieldScatter.", rc)
257 endif
258 endif
259
260 enddo output_loop
261
262 if (allocated(vegt_mdl_one_tile)) deallocate(vegt_mdl_one_tile)
263
264 enddo field_loop
265 enddo time_loop
266
267 status=nf90_close(ncid)
268
269 deallocate(data_mdl_one_tile, mask_mdl_one_tile)
270 deallocate(data_src_global, lat_mdl_one_tile, lon_mdl_one_tile)
271
272 print*,"- CALL FieldRegridRelease."
273 call esmf_fieldregridrelease(routehandle=regrid_data, rc=rc)
274
275 print*,"- CALL FieldDestroy."
276 call esmf_fielddestroy(data_field_src, rc=rc)
277
278 end subroutine interp
279
290 subroutine adjust_for_landice(field, vegt, idim, jdim, field_ch)
291
292 use esmf
293 use mpi
294
295 implicit none
296
297 character(len=*), intent(in) :: field_ch
298
299 integer, intent(in) :: idim, jdim
300
301 real(esmf_kind_i4), intent(in) :: vegt(idim,jdim)
302 real(esmf_kind_r4), intent(inout) :: field(idim,jdim)
303
304 integer, parameter :: landice=15
305
306 integer :: i, j, ierr
307
308 real :: landice_value
309
310 select case (field_ch)
311 case ('substrate_temperature') ! soil substrate temp
312 landice_value = 273.15
313 do j = 1, jdim
314 do i = 1, idim
315 if (nint(vegt(i,j)) == landice) then
316 field(i,j) = min(field(i,j), landice_value)
317 endif
318 enddo
319 enddo
320 case ('vegetation_greenness') ! vegetation greenness
321 landice_value = 0.01 ! 1.0% is bare ground
322 do j = 1, jdim
323 do i = 1, idim
324 if (nint(vegt(i,j)) == landice) then
325 field(i,j) = landice_value
326 endif
327 enddo
328 enddo
329 case ('leaf_area_index') ! leaf area index
330 landice_value = 0.0 ! bare ground
331 do j = 1, jdim
332 do i = 1, idim
333 if (nint(vegt(i,j)) == landice) then
334 field(i,j) = landice_value
335 endif
336 enddo
337 enddo
338 case ('slope_type') ! slope type
339 landice_value = 9.0
340 do j = 1, jdim
341 do i = 1, idim
342 if (nint(vegt(i,j)) == landice) then
343 field(i,j) = landice_value
344 else
345 if (nint(field(i,j)) == nint(landice_value)) field(i,j) = 2.0
346 endif
347 enddo
348 enddo
349 case ('soil_type') ! soil type
350 landice_value = 16.0
351 do j = 1, jdim
352 do i = 1, idim
353 if (nint(vegt(i,j)) == landice) then
354 field(i,j) = landice_value
355 else
356 if (nint(field(i,j)) == nint(landice_value)) field(i,j) = 6.0
357 endif
358 enddo
359 enddo
360 case ('soil_color') ! soil color
361 landice_value = 10.0
362 do j = 1, jdim
363 do i = 1, idim
364 if (nint(vegt(i,j)) == landice) then
365 field(i,j) = landice_value
366 endif
367 enddo
368 enddo
369 case default
370 print*,'- FATAL ERROR IN ROUTINE ADJUST_FOR_LANDICE. UNIDENTIFIED FIELD : ', field_ch
371 call mpi_abort(mpi_comm_world, 57, ierr)
372 end select
373
374 end subroutine adjust_for_landice
adjust_for_landice
subroutine adjust_for_landice(field, vegt, idim, jdim, field_ch)
Ensure consistent fields at land ice points.
Definition interp.F90:291
interp
subroutine interp(localpet, method, input_file)
Read the input source data and interpolate it to the model grid.
Definition interp.F90:14
model_grid
This module defines the model grid.
Definition model_grid.F90:10
model_grid::data_field_mdl
type(esmf_field), public data_field_mdl
ESMF field object that holds the data interpolated to model grid.
Definition model_grid.F90:29
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::missing
real(kind=4), public missing
Value assigned to undefined points (i.e., ocean points for a land field).
Definition model_grid.F90:25
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::longitude_field_mdl
type(esmf_field), public longitude_field_mdl
ESMF field object that holds the model grid longitude.
Definition model_grid.F90:43
model_grid::vegt_field_mdl
type(esmf_field), public vegt_field_mdl
ESMF field object that holds the vegetation type on the model grid.
Definition model_grid.F90:45
program_setup
Set up program execution.
Definition program_setup.f90:17
program_setup::fract_vegsoil_type
logical, public fract_vegsoil_type
When true, output the percentage of each soil and vegetation type category, and the dominant category...
Definition program_setup.f90:51
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::num_time_recs
integer, public num_time_recs
Number of time records.
Definition source_grid.F90:27
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::num_fields
integer, public num_fields
Number of fields in the file.
Definition source_grid.F90:28
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
output
subroutine output(data_one_tile, lat_one_tile, lon_one_tile, i_mdl, j_mdl, tile, record, time, field_idx)
Output model data for a single tile and a single record in netcdf format.
Definition output.f90:20
search
subroutine search(field, mask, idim, jdim, tile, field_name)
Replace undefined values on the model grid with a valid value at a nearby neighbor.
Definition search.f90:24