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