sfc_climo_gen  1.9.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 :: mask_field_mdl
33  type(esmf_field), public :: latitude_field_mdl
35  type(esmf_field), public :: longitude_field_mdl
37  type(esmf_field), public :: vegt_field_mdl
39 
40  public :: define_model_grid
41  public :: model_grid_cleanup
42 
43  contains
44 
53  subroutine define_model_grid(localpet, npets)
54 
55  use esmf
56  use netcdf
57  use program_setup
58  use utils
59  use mpi
60 
61  implicit none
62 
63  integer, intent(in) :: localpet, npets
64 
65  character(len=500) :: the_file
66 
67  integer :: error, id_dim, id_tiles, ncid
68  integer :: id_grid_tiles, ierr
69  integer :: extra, rc, tile
70  integer, allocatable :: decomptile(:,:)
71 
72  integer(esmf_kind_i4), allocatable :: mask_mdl_one_tile(:,:)
73  integer(esmf_kind_i4), pointer :: mask_field_mdl_ptr(:,:)
74  integer(esmf_kind_i4), pointer :: mask_mdl_ptr(:,:)
75 
76  real(esmf_kind_r4), allocatable :: latitude_one_tile(:,:)
77  real(esmf_kind_r4), allocatable :: longitude_one_tile(:,:)
78 
79 !-----------------------------------------------------------------------
80 ! Get the number of tiles from the mosaic file.
81 !-----------------------------------------------------------------------
82 
83  print*,'- OPEN MODEL GRID MOSAIC FILE: ',trim(mosaic_file_mdl)
84  error=nf90_open(trim(mosaic_file_mdl),nf90_nowrite,ncid)
85  call netcdf_err(error, "OPENING MODEL GRID MOSAIC FILE")
86 
87  print*,"- READ NUMBER OF TILES"
88  error=nf90_inq_dimid(ncid, 'ntiles', id_tiles)
89  call netcdf_err(error, "READING NTILES ID")
90  error=nf90_inquire_dimension(ncid,id_tiles,len=num_tiles)
91  call netcdf_err(error, "READING NTILES")
92  error=nf90_inq_varid(ncid, 'gridtiles', id_grid_tiles)
93  call netcdf_err(error, "READING GRIDTILES ID")
94  allocate(grid_tiles(num_tiles))
95  grid_tiles="NULL"
96  print*,"- READ TILE NAMES"
97  error=nf90_get_var(ncid, id_grid_tiles, grid_tiles)
98  call netcdf_err(error, "READING GRIDTILES")
99 
100  error = nf90_close(ncid)
101 
102  print*,'- NUMBER OF TILES, MODEL GRID IS ', num_tiles
103 
104  if (mod(npets,num_tiles) /= 0) then
105  print*,'- FATAL ERROR: MUST RUN THIS PROGRAM WITH A TASK COUNT THAT'
106  print*,'- IS A MULTIPLE OF THE NUMBER OF TILES.'
107  call mpi_abort(mpi_comm_world, 44, ierr)
108  endif
109 
110 !-----------------------------------------------------------------------
111 ! Get the model grid specs and land mask from the orography files.
112 !-----------------------------------------------------------------------
113 
114  orog_dir_mdl = trim(orog_dir_mdl) // '/'
115  the_file = trim(orog_dir_mdl) // trim(orog_files_mdl(1))
116 
117  print*,'- OPEN FIRST MODEL GRID OROGRAPHY FILE: ',trim(the_file)
118  error=nf90_open(trim(the_file),nf90_nowrite,ncid)
119  call netcdf_err(error, "OPENING MODEL GRID OROGRAPHY FILE")
120  print*,"- READ GRID DIMENSIONS"
121  error=nf90_inq_dimid(ncid, 'lon', id_dim)
122  call netcdf_err(error, "READING MODEL LON ID")
123  error=nf90_inquire_dimension(ncid,id_dim,len=i_mdl)
124  call netcdf_err(error, "READING MODEL LON")
125  error=nf90_inq_dimid(ncid, 'lat', id_dim)
126  call netcdf_err(error, "READING MODEL LAT ID")
127  error=nf90_inquire_dimension(ncid,id_dim,len=j_mdl)
128  call netcdf_err(error, "READING MODEL LAT")
129  error = nf90_close(ncid)
130 
131  print*,"- I/J DIMENSIONS OF THE MODEL GRID TILES ", i_mdl, j_mdl
132 
133  ij_mdl = i_mdl * j_mdl
134 
135 !-----------------------------------------------------------------------
136 ! Create ESMF grid object for the model grid.
137 !-----------------------------------------------------------------------
138 
139  extra = npets / num_tiles
140 
141  allocate(decomptile(2,num_tiles))
142 
143  do tile = 1, num_tiles
144  decomptile(:,tile)=(/1,extra/)
145  enddo
146 
147  print*,"- CALL GridCreateMosaic FOR MODEL GRID"
148  grid_mdl = esmf_gridcreatemosaic(filename=trim(mosaic_file_mdl), &
149  regdecompptile=decomptile, &
150  staggerloclist=(/esmf_staggerloc_center, &
151  esmf_staggerloc_corner/), &
152  indexflag=esmf_index_global, &
153  tilefilepath=trim(orog_dir_mdl), &
154  rc=rc)
155  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
156  call error_handler("IN GridCreateMosaic", rc)
157 
158  print*,"- CALL FieldCreate FOR DATA INTERPOLATED TO MODEL GRID."
159  data_field_mdl = esmf_fieldcreate(grid_mdl, &
160  typekind=esmf_typekind_r4, &
161  staggerloc=esmf_staggerloc_center, &
162  name="data on model grid", &
163  rc=rc)
164  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
165  call error_handler("IN FieldCreate", rc)
166 
167  print*,"- CALL FieldCreate FOR VEGETATION TYPE INTERPOLATED TO MODEL GRID."
168  vegt_field_mdl = esmf_fieldcreate(grid_mdl, &
169  typekind=esmf_typekind_r4, &
170  staggerloc=esmf_staggerloc_center, &
171  name="veg type 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  print*,"- CALL FieldCreate FOR MODEL GRID LATITUDE."
177  latitude_field_mdl = esmf_fieldcreate(grid_mdl, &
178  typekind=esmf_typekind_r4, &
179  staggerloc=esmf_staggerloc_center, &
180  name="latitude on model grid", &
181  rc=rc)
182  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
183  call error_handler("IN FieldCreate", rc)
184 
185  print*,"- CALL FieldCreate FOR MODEL GRID LONGITUDE."
186  longitude_field_mdl = esmf_fieldcreate(grid_mdl, &
187  typekind=esmf_typekind_r4, &
188  staggerloc=esmf_staggerloc_center, &
189  name="longitude on model grid", &
190  rc=rc)
191  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
192  call error_handler("IN FieldCreate", rc)
193 
194 !-----------------------------------------------------------------------
195 ! Set model land mask.
196 !-----------------------------------------------------------------------
197 
198  print*,"- CALL FieldCreate FOR MODEL GRID LANDMASK."
199  mask_field_mdl = esmf_fieldcreate(grid_mdl, &
200  typekind=esmf_typekind_i4, &
201  staggerloc=esmf_staggerloc_center, &
202  name="model grid land mask", &
203  rc=rc)
204  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
205  call error_handler("IN FieldCreate", rc)
206 
207  print*,"- CALL FieldGet FOR MODEL GRID LANDMASK."
208  nullify(mask_field_mdl_ptr)
209  call esmf_fieldget(mask_field_mdl, &
210  farrayptr=mask_field_mdl_ptr, &
211  rc=rc)
212  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
213  call error_handler("IN FieldGet", rc)
214 
215  if (localpet == 0) then
216  allocate(mask_mdl_one_tile(i_mdl,j_mdl))
217  allocate(latitude_one_tile(i_mdl,j_mdl))
218  allocate(longitude_one_tile(i_mdl,j_mdl))
219  else
220  allocate(mask_mdl_one_tile(0,0))
221  allocate(latitude_one_tile(0,0))
222  allocate(longitude_one_tile(0,0))
223  endif
224 
225  do tile = 1, num_tiles
226  if (localpet == 0) then
227  the_file = trim(orog_dir_mdl) // trim(orog_files_mdl(tile))
228  call get_model_info(trim(the_file), mask_mdl_one_tile, &
229  latitude_one_tile, longitude_one_tile, i_mdl, j_mdl)
230  endif
231 
232  print*,"- CALL FieldScatter FOR MODEL GRID MASK. TILE IS: ", tile
233  call esmf_fieldscatter(mask_field_mdl, mask_mdl_one_tile, rootpet=0, tile=tile, rc=rc)
234  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
235  call error_handler("IN FieldScatter", rc)
236 
237  print*,"- CALL FieldScatter FOR MODEL LATITUDE. TILE IS: ", tile
238  call esmf_fieldscatter(latitude_field_mdl, latitude_one_tile, rootpet=0, tile=tile, rc=rc)
239  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
240  call error_handler("IN FieldScatter", rc)
241 
242  print*,"- CALL FieldScatter FOR MODEL LONGITUDE. TILE IS: ", tile
243  call esmf_fieldscatter(longitude_field_mdl, longitude_one_tile, rootpet=0, tile=tile, rc=rc)
244  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
245  call error_handler("IN FieldScatter", rc)
246 
247  enddo
248 
249  deallocate(mask_mdl_one_tile, latitude_one_tile, longitude_one_tile)
250 
251  print*,"- CALL GridAddItem FOR MODEL GRID."
252  call esmf_gridadditem(grid_mdl, &
253  itemflag=esmf_griditem_mask, &
254  staggerloc=esmf_staggerloc_center, &
255  rc=rc)
256  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
257  call error_handler("IN GridAddItem", rc)
258 
259  print*,"- CALL GridGetItem FOR MODEL GRID."
260  nullify(mask_mdl_ptr)
261  call esmf_gridgetitem(grid_mdl, &
262  itemflag=esmf_griditem_mask, &
263  farrayptr=mask_mdl_ptr, &
264  rc=rc)
265  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
266  call error_handler("IN GridGetItem", rc)
267 
268  mask_mdl_ptr = mask_field_mdl_ptr
269 
270  end subroutine define_model_grid
271 
283  subroutine get_model_info(orog_file, mask, lat2d, lon2d, idim, jdim)
284 
285  use esmf
286  use netcdf
287  use utils
288 
289  implicit none
290 
291  character(len=*), intent(in) :: orog_file
292 
293  integer, intent(in) :: idim, jdim
294  integer(esmf_kind_i4), intent(out) :: mask(idim,jdim)
295 
296  real(esmf_kind_r4), intent(out) :: lat2d(idim,jdim)
297  real(esmf_kind_r4), intent(out) :: lon2d(idim,jdim)
298 
299  integer :: error, lat, lon, i, j
300  integer :: ncid, id_dim, id_var
301 
302  real(kind=4), allocatable :: dummy(:,:)
303 
304  print*,"- READ MODEL OROGRAPHY FILE"
305 
306  print*,'- OPEN FILE: ', orog_file
307  error=nf90_open(orog_file,nf90_nowrite,ncid)
308  call netcdf_err(error, "OPENING MODEL OROGRAPHY FILE")
309 
310  print*,"- READ I-DIMENSION"
311  error=nf90_inq_dimid(ncid, 'lon', id_dim)
312  call netcdf_err(error, "READING LON ID")
313  error=nf90_inquire_dimension(ncid,id_dim,len=lon)
314  call netcdf_err(error, "READING LON")
315 
316  print*,"- READ J-DIMENSION"
317  error=nf90_inq_dimid(ncid, 'lat', id_dim)
318  call netcdf_err(error, "READING LAT ID")
319  error=nf90_inquire_dimension(ncid,id_dim,len=lat)
320  call netcdf_err(error, "READING LAT")
321 
322  print*,"- I/J DIMENSIONS: ", lon, lat
323 
324  if ((lon /= idim) .or. (lat /= jdim)) then
325  call error_handler("MISMATCH IN DIMENSIONS.")
326  endif
327 
328  allocate(dummy(idim,jdim))
329 
330 !-----------------------------------------------------------------------
331 ! If the lake maker was used, there will be a 'lake_frac' record.
332 ! In that case, land/non-land is determined by 'land_frac'.
333 !
334 ! If the lake maker was not used, use 'slmsk', which is defined
335 ! as the nint(land_frac).
336 !-----------------------------------------------------------------------
337 
338  error=nf90_inq_varid(ncid, 'lake_frac', id_var)
339  if (error /= 0) then
340  print*,"- READ LAND MASK (SLMSK)"
341  error=nf90_inq_varid(ncid, 'slmsk', id_var)
342  call netcdf_err(error, "READING SLMSK ID")
343  error=nf90_get_var(ncid, id_var, dummy)
344  call netcdf_err(error, "READING SLMSK")
345  mask = nint(dummy)
346  else
347  print*,"- READ LAND FRACTION"
348  error=nf90_inq_varid(ncid, 'land_frac', id_var)
349  call netcdf_err(error, "READING LAND_FRAC ID")
350  error=nf90_get_var(ncid, id_var, dummy)
351  call netcdf_err(error, "READING LAND_FRAC")
352  mask = 0
353  do j = 1, lat
354  do i = 1, lon
355  if (dummy(i,j) > 0.0) then
356  mask(i,j) = 1
357  endif
358  enddo
359  enddo
360  endif
361 
362  print*,"- READ LATITUDE"
363  error=nf90_inq_varid(ncid, 'geolat', id_var)
364  call netcdf_err(error, "READING GEOLAT ID")
365  error=nf90_get_var(ncid, id_var, dummy)
366  call netcdf_err(error, "READING GEOLAT")
367  lat2d=dummy
368 
369  print*,"- READ LONGITUDE"
370  error=nf90_inq_varid(ncid, 'geolon', id_var)
371  call netcdf_err(error, "READING GEOLON ID")
372  error=nf90_get_var(ncid, id_var, dummy)
373  call netcdf_err(error, "READING GEOLON")
374  lon2d=dummy
375 
376  error = nf90_close(ncid)
377 
378  deallocate (dummy)
379 
380  end subroutine get_model_info
381 
387  subroutine model_grid_cleanup
388 
389  implicit none
390 
391  integer :: rc
392 
393  print*,"- CALL GridDestroy FOR MODEL GRID."
394  call esmf_griddestroy(grid_mdl,rc=rc)
395 
396  print*,"- CALL FieldDestroy FOR MODEL GRID LAND MASK."
397  call esmf_fielddestroy(mask_field_mdl,rc=rc)
398 
399  print*,"- CALL FieldDestroy FOR MODEL GRID DATA FIELD."
400  call esmf_fielddestroy(data_field_mdl,rc=rc)
401 
402  print*,"- CALL FieldDestroy FOR MODEL GRID VEGETATION TYPE."
403  call esmf_fielddestroy(vegt_field_mdl,rc=rc)
404 
405  print*,"- CALL FieldDestroy FOR MODEL GRID LATITUDE."
406  call esmf_fielddestroy(latitude_field_mdl,rc=rc)
407 
408  print*,"- CALL FieldDestroy FOR MODEL GRID LONGITUDE."
409  call esmf_fielddestroy(longitude_field_mdl,rc=rc)
410 
411  end subroutine model_grid_cleanup
412 
413  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:37
model_grid::model_grid_cleanup
subroutine, public model_grid_cleanup
Model grid cleanup.
Definition: model_grid.F90:388
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:54
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:35
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::latitude_field_mdl
type(esmf_field), public latitude_field_mdl
ESMF field object that holds the model grid latitude.
Definition: model_grid.F90:33
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:31
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
model_grid::get_model_info
subroutine get_model_info(orog_file, mask, lat2d, lon2d, idim, jdim)
Get model information.
Definition: model_grid.F90:284
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