sfc_climo_gen  1.10.0
model_grid.F90
Go to the documentation of this file.
1 
4 
10  module model_grid
11 
12  use esmf
13 
14  implicit none
15 
16  private
17 
18  character(len=5), allocatable, public :: grid_tiles(:)
19 
20  integer, public :: i_mdl
21  integer, public :: j_mdl
22  integer, public :: ij_mdl
23  integer, public :: num_tiles
24 
25  real(kind=4), public :: missing = -999.
27 
28  type(esmf_grid), public :: grid_mdl
29  type(esmf_field), public :: data_field_mdl
31  type(esmf_field), public :: land_frac_field_mdl
37  type(esmf_field), public :: mask_field_mdl
41  type(esmf_field), public :: latitude_field_mdl
43  type(esmf_field), public :: longitude_field_mdl
45  type(esmf_field), public :: vegt_field_mdl
47 
48  public :: define_model_grid
49  public :: model_grid_cleanup
50 
51  contains
52 
61  subroutine define_model_grid(localpet, npets)
62 
63  use esmf
64  use netcdf
65  use program_setup
66  use utils
67  use mpi
68 
69  implicit none
70 
71  integer, intent(in) :: localpet, npets
72 
73  character(len=500) :: the_file
74 
75  integer :: error, id_dim, id_tiles, ncid
76  integer :: id_grid_tiles, ierr
77  integer :: extra, rc, tile
78  integer, allocatable :: decomptile(:,:)
79 
80  integer(esmf_kind_i4), allocatable :: mask_mdl_one_tile(:,:)
81  integer(esmf_kind_i4), pointer :: mask_field_mdl_ptr(:,:)
82  integer(esmf_kind_i4), pointer :: mask_mdl_ptr(:,:)
83 
84  real(esmf_kind_r4), allocatable :: latitude_one_tile(:,:)
85  real(esmf_kind_r4), allocatable :: longitude_one_tile(:,:)
86  real(esmf_kind_r4), allocatable :: land_frac_one_tile(:,:)
87 
88 !-----------------------------------------------------------------------
89 ! Get the number of tiles from the mosaic file.
90 !-----------------------------------------------------------------------
91 
92  print*,'- OPEN MODEL GRID MOSAIC FILE: ',trim(mosaic_file_mdl)
93  error=nf90_open(trim(mosaic_file_mdl),nf90_nowrite,ncid)
94  call netcdf_err(error, "OPENING MODEL GRID MOSAIC FILE")
95 
96  print*,"- READ NUMBER OF TILES"
97  error=nf90_inq_dimid(ncid, 'ntiles', id_tiles)
98  call netcdf_err(error, "READING NTILES ID")
99  error=nf90_inquire_dimension(ncid,id_tiles,len=num_tiles)
100  call netcdf_err(error, "READING NTILES")
101  error=nf90_inq_varid(ncid, 'gridtiles', id_grid_tiles)
102  call netcdf_err(error, "READING GRIDTILES ID")
103  allocate(grid_tiles(num_tiles))
104  grid_tiles="NULL"
105  print*,"- READ TILE NAMES"
106  error=nf90_get_var(ncid, id_grid_tiles, grid_tiles)
107  call netcdf_err(error, "READING GRIDTILES")
108 
109  error = nf90_close(ncid)
110 
111  print*,'- NUMBER OF TILES, MODEL GRID IS ', num_tiles
112 
113  if (mod(npets,num_tiles) /= 0) then
114  print*,'- FATAL ERROR: MUST RUN THIS PROGRAM WITH A TASK COUNT THAT'
115  print*,'- IS A MULTIPLE OF THE NUMBER OF TILES.'
116  call mpi_abort(mpi_comm_world, 44, ierr)
117  endif
118 
119 !-----------------------------------------------------------------------
120 ! Get the model grid specs and land mask from the orography files.
121 !-----------------------------------------------------------------------
122 
123  orog_dir_mdl = trim(orog_dir_mdl) // '/'
124  the_file = trim(orog_dir_mdl) // trim(orog_files_mdl(1))
125 
126  print*,'- OPEN FIRST MODEL GRID OROGRAPHY FILE: ',trim(the_file)
127  error=nf90_open(trim(the_file),nf90_nowrite,ncid)
128  call netcdf_err(error, "OPENING MODEL GRID OROGRAPHY FILE")
129  print*,"- READ GRID DIMENSIONS"
130  error=nf90_inq_dimid(ncid, 'lon', id_dim)
131  call netcdf_err(error, "READING MODEL LON ID")
132  error=nf90_inquire_dimension(ncid,id_dim,len=i_mdl)
133  call netcdf_err(error, "READING MODEL LON")
134  error=nf90_inq_dimid(ncid, 'lat', id_dim)
135  call netcdf_err(error, "READING MODEL LAT ID")
136  error=nf90_inquire_dimension(ncid,id_dim,len=j_mdl)
137  call netcdf_err(error, "READING MODEL LAT")
138  error = nf90_close(ncid)
139 
140  print*,"- I/J DIMENSIONS OF THE MODEL GRID TILES ", i_mdl, j_mdl
141 
142  ij_mdl = i_mdl * j_mdl
143 
144 !-----------------------------------------------------------------------
145 ! Create ESMF grid object for the model grid.
146 !-----------------------------------------------------------------------
147 
148  extra = npets / num_tiles
149 
150  allocate(decomptile(2,num_tiles))
151 
152  do tile = 1, num_tiles
153  decomptile(:,tile)=(/1,extra/)
154  enddo
155 
156  print*,"- CALL GridCreateMosaic FOR MODEL GRID"
157  grid_mdl = esmf_gridcreatemosaic(filename=trim(mosaic_file_mdl), &
158  regdecompptile=decomptile, &
159  staggerloclist=(/esmf_staggerloc_center, &
160  esmf_staggerloc_corner/), &
161  indexflag=esmf_index_global, &
162  tilefilepath=trim(orog_dir_mdl), &
163  rc=rc)
164  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
165  call error_handler("IN GridCreateMosaic", rc)
166 
167  print*,"- CALL FieldCreate FOR DATA INTERPOLATED TO MODEL GRID."
168  data_field_mdl = esmf_fieldcreate(grid_mdl, &
169  typekind=esmf_typekind_r4, &
170  staggerloc=esmf_staggerloc_center, &
171  name="data on model grid", &
172  rc=rc)
173  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
174  call error_handler("IN FieldCreate", rc)
175 
176  if (.not. fract_vegsoil_type) then
177  print*,"- CALL FieldCreate FOR VEGETATION TYPE INTERPOLATED TO MODEL GRID."
178  vegt_field_mdl = esmf_fieldcreate(grid_mdl, &
179  typekind=esmf_typekind_r4, &
180  staggerloc=esmf_staggerloc_center, &
181  name="veg type on model grid", &
182  rc=rc)
183  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
184  call error_handler("IN FieldCreate", rc)
185  endif
186 
187  print*,"- CALL FieldCreate FOR MODEL GRID LATITUDE."
188  latitude_field_mdl = esmf_fieldcreate(grid_mdl, &
189  typekind=esmf_typekind_r4, &
190  staggerloc=esmf_staggerloc_center, &
191  name="latitude on model grid", &
192  rc=rc)
193  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
194  call error_handler("IN FieldCreate", rc)
195 
196  print*,"- CALL FieldCreate FOR MODEL GRID LONGITUDE."
197  longitude_field_mdl = esmf_fieldcreate(grid_mdl, &
198  typekind=esmf_typekind_r4, &
199  staggerloc=esmf_staggerloc_center, &
200  name="longitude on model grid", &
201  rc=rc)
202  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
203  call error_handler("IN FieldCreate", rc)
204 
205 !-----------------------------------------------------------------------
206 ! Set model land mask.
207 !-----------------------------------------------------------------------
208 
209  print*,"- CALL FieldCreate FOR MODEL GRID LANDMASK."
210  mask_field_mdl = esmf_fieldcreate(grid_mdl, &
211  typekind=esmf_typekind_i4, &
212  staggerloc=esmf_staggerloc_center, &
213  name="model grid land mask", &
214  rc=rc)
215  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
216  call error_handler("IN FieldCreate", rc)
217 
218  print*,"- CALL FieldGet FOR MODEL GRID LANDMASK."
219  nullify(mask_field_mdl_ptr)
220  call esmf_fieldget(mask_field_mdl, &
221  farrayptr=mask_field_mdl_ptr, &
222  rc=rc)
223  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
224  call error_handler("IN FieldGet", rc)
225 
226  print*,"- CALL FieldCreate FOR MODEL GRID LAND FRACTION."
227  land_frac_field_mdl = esmf_fieldcreate(grid_mdl, &
228  typekind=esmf_typekind_r4, &
229  staggerloc=esmf_staggerloc_center, &
230  name="model grid land fraction", &
231  rc=rc)
232  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
233  call error_handler("IN FieldCreate", rc)
234 
235  if (localpet == 0) then
236  allocate(mask_mdl_one_tile(i_mdl,j_mdl))
237  allocate(land_frac_one_tile(i_mdl,j_mdl))
238  allocate(latitude_one_tile(i_mdl,j_mdl))
239  allocate(longitude_one_tile(i_mdl,j_mdl))
240  else
241  allocate(mask_mdl_one_tile(0,0))
242  allocate(land_frac_one_tile(0,0))
243  allocate(latitude_one_tile(0,0))
244  allocate(longitude_one_tile(0,0))
245  endif
246 
247  do tile = 1, num_tiles
248  if (localpet == 0) then
249  the_file = trim(orog_dir_mdl) // trim(orog_files_mdl(tile))
250  call get_model_info(trim(the_file), mask_mdl_one_tile, land_frac_one_tile, &
251  latitude_one_tile, longitude_one_tile, i_mdl, j_mdl)
252  endif
253 
254  print*,"- CALL FieldScatter FOR MODEL GRID MASK. TILE IS: ", tile
255  call esmf_fieldscatter(mask_field_mdl, mask_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
256  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
257  call error_handler("IN FieldScatter", rc)
258 
259  print*,"- CALL FieldScatter FOR MODEL GRID LAND FRACTION. TILE IS: ", tile
260  call esmf_fieldscatter(land_frac_field_mdl, land_frac_one_tile, rootpet=0, tile=tile, rc=rc)
261  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
262  call error_handler("IN FieldScatter", rc)
263 
264  print*,"- CALL FieldScatter FOR MODEL LATITUDE. TILE IS: ", tile
265  call esmf_fieldscatter(latitude_field_mdl, latitude_one_tile, rootpet=0, tile=tile, rc=rc)
266  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
267  call error_handler("IN FieldScatter", rc)
268 
269  print*,"- CALL FieldScatter FOR MODEL LONGITUDE. TILE IS: ", tile
270  call esmf_fieldscatter(longitude_field_mdl, longitude_one_tile, rootpet=0, tile=tile, rc=rc)
271  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
272  call error_handler("IN FieldScatter", rc)
273 
274  enddo
275 
276  deallocate(mask_mdl_one_tile, latitude_one_tile, longitude_one_tile, land_frac_one_tile)
277 
278  print*,"- CALL GridAddItem FOR MODEL GRID."
279  call esmf_gridadditem(grid_mdl, &
280  itemflag=esmf_griditem_mask, &
281  staggerloc=esmf_staggerloc_center, &
282  rc=rc)
283  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
284  call error_handler("IN GridAddItem", rc)
285 
286  print*,"- CALL GridGetItem FOR MODEL GRID."
287  nullify(mask_mdl_ptr)
288  call esmf_gridgetitem(grid_mdl, &
289  itemflag=esmf_griditem_mask, &
290  farrayptr=mask_mdl_ptr, &
291  rc=rc)
292  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
293  call error_handler("IN GridGetItem", rc)
294 
295  mask_mdl_ptr = mask_field_mdl_ptr
296 
297  end subroutine define_model_grid
298 
311  subroutine get_model_info(orog_file, mask, land_frac, lat2d, lon2d, idim, jdim)
312 
313  use esmf
314  use netcdf
315  use utils
316 
317  implicit none
318 
319  character(len=*), intent(in) :: orog_file
320 
321  integer, intent(in) :: idim, jdim
322  integer(esmf_kind_i4), intent(out) :: mask(idim,jdim)
323 
324  real(esmf_kind_r4), intent(out) :: lat2d(idim,jdim)
325  real(esmf_kind_r4), intent(out) :: lon2d(idim,jdim)
326  real(esmf_kind_r4), intent(out) :: land_frac(idim,jdim)
327 
328  integer :: error, lat, lon, i, j
329  integer :: ncid, id_dim, id_var
330 
331  real(kind=4), allocatable :: dummy(:,:)
332 
333  print*,"- READ MODEL OROGRAPHY FILE"
334 
335  print*,'- OPEN FILE: ', orog_file
336  error=nf90_open(orog_file,nf90_nowrite,ncid)
337  call netcdf_err(error, "OPENING MODEL OROGRAPHY FILE")
338 
339  print*,"- READ I-DIMENSION"
340  error=nf90_inq_dimid(ncid, 'lon', id_dim)
341  call netcdf_err(error, "READING LON ID")
342  error=nf90_inquire_dimension(ncid,id_dim,len=lon)
343  call netcdf_err(error, "READING LON")
344 
345  print*,"- READ J-DIMENSION"
346  error=nf90_inq_dimid(ncid, 'lat', id_dim)
347  call netcdf_err(error, "READING LAT ID")
348  error=nf90_inquire_dimension(ncid,id_dim,len=lat)
349  call netcdf_err(error, "READING LAT")
350 
351  print*,"- I/J DIMENSIONS: ", lon, lat
352 
353  if ((lon /= idim) .or. (lat /= jdim)) then
354  call error_handler("MISMATCH IN DIMENSIONS.")
355  endif
356 
357  allocate(dummy(idim,jdim))
358 
359 !-----------------------------------------------------------------------
360 ! If the lake maker was used, we are running with a fractional
361 ! land/non-land grid and there will be a 'lake_frac' record.
362 ! In that case, land/non-land is determined by 'land_frac'.
363 !
364 ! If the lake maker was not used, use 'slmsk', which is defined
365 ! as the nint(land_frac).
366 !
367 ! In summary, if 'mask' is one, the point is all land or
368 ! partial land and surface data will be mapped to it. Otherwise,
369 ! when 'mask' is zero, then the point is all non-land and
370 ! surface data will not be mapped to it.
371 !-----------------------------------------------------------------------
372 
373  error=nf90_inq_varid(ncid, 'lake_frac', id_var)
374  if (error /= 0) then
375  print*,"- READ LAND MASK (SLMSK)"
376  error=nf90_inq_varid(ncid, 'slmsk', id_var)
377  call netcdf_err(error, "READING SLMSK ID")
378  error=nf90_get_var(ncid, id_var, dummy)
379  call netcdf_err(error, "READING SLMSK")
380  mask = nint(dummy)
381  land_frac = -999.
382  else
383  print*,"- READ LAND FRACTION"
384  error=nf90_inq_varid(ncid, 'land_frac', id_var)
385  call netcdf_err(error, "READING LAND_FRAC ID")
386  error=nf90_get_var(ncid, id_var, land_frac)
387  call netcdf_err(error, "READING LAND_FRAC")
388  mask = 0
389  do j = 1, lat
390  do i = 1, lon
391  if (land_frac(i,j) > 0.0) then
392  mask(i,j) = 1
393  endif
394  enddo
395  enddo
396  endif
397 
398  print*,"- READ LATITUDE"
399  error=nf90_inq_varid(ncid, 'geolat', id_var)
400  call netcdf_err(error, "READING GEOLAT ID")
401  error=nf90_get_var(ncid, id_var, dummy)
402  call netcdf_err(error, "READING GEOLAT")
403  lat2d=dummy
404 
405  print*,"- READ LONGITUDE"
406  error=nf90_inq_varid(ncid, 'geolon', id_var)
407  call netcdf_err(error, "READING GEOLON ID")
408  error=nf90_get_var(ncid, id_var, dummy)
409  call netcdf_err(error, "READING GEOLON")
410  lon2d=dummy
411 
412  error = nf90_close(ncid)
413 
414  deallocate (dummy)
415 
416  end subroutine get_model_info
417 
423  subroutine model_grid_cleanup
424 
425  implicit none
426 
427  integer :: rc
428 
429  print*,"- CALL GridDestroy FOR MODEL GRID."
430  call esmf_griddestroy(grid_mdl,rc=rc)
431 
432  print*,"- CALL FieldDestroy FOR MODEL GRID LAND MASK."
433  call esmf_fielddestroy(mask_field_mdl,rc=rc)
434 
435  print*,"- CALL FieldDestroy FOR MODEL GRID LAND FRACTION."
436  call esmf_fielddestroy(land_frac_field_mdl,rc=rc)
437 
438  print*,"- CALL FieldDestroy FOR MODEL GRID DATA FIELD."
439  call esmf_fielddestroy(data_field_mdl,rc=rc)
440 
441  if (esmf_fieldiscreated(vegt_field_mdl)) then
442  print*,"- CALL FieldDestroy FOR MODEL GRID VEGETATION TYPE."
443  call esmf_fielddestroy(vegt_field_mdl,rc=rc)
444  endif
445 
446  print*,"- CALL FieldDestroy FOR MODEL GRID LATITUDE."
447  call esmf_fielddestroy(latitude_field_mdl,rc=rc)
448 
449  print*,"- CALL FieldDestroy FOR MODEL GRID LONGITUDE."
450  call esmf_fielddestroy(longitude_field_mdl,rc=rc)
451 
452  end subroutine model_grid_cleanup
453 
454  end module model_grid
program_setup
Set up program execution.
Definition: program_setup.f90:17
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
model_grid::model_grid_cleanup
subroutine, public model_grid_cleanup
Model grid cleanup.
Definition: model_grid.F90:424
program_setup::orog_files_mdl
character(len=500), dimension(6), public orog_files_mdl
Model grid orography filenames.
Definition: program_setup.f90:37
model_grid::grid_tiles
character(len=5), dimension(:), allocatable, public grid_tiles
Array of model grid tile names.
Definition: model_grid.F90:18
model_grid::define_model_grid
subroutine, public define_model_grid(localpet, npets)
Define model grid.
Definition: model_grid.F90:62
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::longitude_field_mdl
type(esmf_field), public longitude_field_mdl
ESMF field object that holds the model grid longitude.
Definition: model_grid.F90:43
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
model_grid
This module defines the model grid.
Definition: model_grid.F90:10
model_grid::j_mdl
integer, public j_mdl
j dimension of model tile.
Definition: model_grid.F90:21
model_grid::get_model_info
subroutine get_model_info(orog_file, mask, land_frac, lat2d, lon2d, idim, jdim)
Get model information.
Definition: model_grid.F90:312
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::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::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::num_tiles
integer, public num_tiles
Total number of model grid tiles.
Definition: model_grid.F90:23
model_grid::grid_mdl
type(esmf_grid), public grid_mdl
ESMF grid object for the model grid.
Definition: model_grid.F90:28
utils::error_handler
subroutine, public error_handler(string, rc)
Handle errors.
Definition: utils.f90:49
model_grid::ij_mdl
integer, public ij_mdl
Total number of points on a model tile.
Definition: model_grid.F90:22
utils::netcdf_err
subroutine, public netcdf_err(err, string)
Handle netCDF error codes.
Definition: utils.f90:23
utils
Utilities.
Definition: utils.f90:8
program_setup::orog_dir_mdl
character(len=500), public orog_dir_mdl
Directory containing the model grid orography files.
Definition: program_setup.f90:36
program_setup::mosaic_file_mdl
character(len=500), public mosaic_file_mdl
Model grid mosaic file.
Definition: program_setup.f90:35
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