sfc_climo_gen 1.14.0
Loading...
Searching...
No Matches
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
model_grid
This module defines the model grid.
Definition model_grid.F90:10
model_grid::define_model_grid
subroutine, public define_model_grid(localpet, npets)
Define model grid.
Definition model_grid.F90:62
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::grid_mdl
type(esmf_grid), public grid_mdl
ESMF grid object for the model grid.
Definition model_grid.F90:28
model_grid::i_mdl
integer, public i_mdl
i dimension of model tile.
Definition model_grid.F90:20
model_grid::ij_mdl
integer, public ij_mdl
Total number of points on a model tile.
Definition model_grid.F90:22
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::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::model_grid_cleanup
subroutine, public model_grid_cleanup
Model grid cleanup.
Definition model_grid.F90:424
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::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
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::mosaic_file_mdl
character(len=500), public mosaic_file_mdl
Model grid mosaic file.
Definition program_setup.f90:35
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
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::orog_files_mdl
character(len=500), dimension(6), public orog_files_mdl
Model grid orography filenames.
Definition program_setup.f90:37
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