sfc_climo_gen  1.13.0
 All Data Structures Files Functions Variables Pages
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
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:290
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
program_setup
Set up program execution.
Definition: program_setup.f90:17
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