chgres_cube 1.14.0
Loading...
Searching...
No Matches
model_grid.F90
Go to the documentation of this file.
1
4
9 module model_grid
10
11 use esmf
12 use esmf_logpublicmod
13
14 use utilities, only : error_handler, netcdf_err
15 implicit none
16
17 private
18
19 character(len=5), allocatable, public :: tiles_target_grid(:)
20
21 character(len=50), public :: input_grid_type = "latlon"
22
23
24 ! Made lsoil_target non-parameter to allow for RAP land surface initiation
25 integer, public :: lsoil_target = 4 ! # soil layers
26
27 integer, public :: i_input
28
30 integer, public :: j_input
31
33 integer, public :: ip1_input
34
35 integer, public :: jp1_input
36
37 integer, public :: i_target
38
40 integer, public :: j_target
41
43 integer, public :: ip1_target
44
45 integer, public :: jp1_target
46
47 integer, public :: num_tiles_input_grid
48
49 integer, public :: num_tiles_target_grid
50
51
52 type(esmf_grid), public :: input_grid
53
54 type(esmf_grid), public :: target_grid
55
56
57 type(esmf_field), public :: latitude_input_grid
58
59 type(esmf_field), public :: longitude_input_grid
60
61 type(esmf_field), public :: latitude_s_input_grid
62
64 type(esmf_field), public :: longitude_s_input_grid
65
67 type(esmf_field), public :: latitude_w_input_grid
68
70 type(esmf_field), public :: longitude_w_input_grid
71
73 type(esmf_field), public :: landmask_target_grid
74
76 type(esmf_field), public :: land_frac_target_grid
77
78 type(esmf_field), public :: latitude_target_grid
79
80 type(esmf_field), public :: latitude_s_target_grid
81
83 type(esmf_field), public :: latitude_w_target_grid
84
86 type(esmf_field), public :: longitude_target_grid
87
88 type(esmf_field), public :: longitude_s_target_grid
89
91 type(esmf_field), public :: longitude_w_target_grid
92
94 type(esmf_field), public :: seamask_target_grid
95
97 type(esmf_field), public :: terrain_target_grid
98
99
100 public :: define_target_grid
101 public :: define_input_grid
102 public :: cleanup_input_target_grid_data
103
104 contains
105
117 subroutine define_input_grid(localpet, npets)
118
119 use program_setup, only : input_type
120
121 implicit none
122
123 integer, intent(in) :: localpet, npets
124
125 if (trim(input_type) == "gaussian_nemsio" .or. &
126 trim(input_type) == "gfs_gaussian_nemsio" .or. &
127 trim(input_type) == "gfs_sigio" .or. &
128 trim(input_type) == "gaussian_netcdf") then
129 call define_input_grid_gaussian(localpet,npets)
130 elseif (trim(input_type) == "grib2") then
131 call define_input_grid_grib2(localpet, npets)
132 else
133 call define_input_grid_mosaic(localpet, npets)
134 endif
135
136 end subroutine define_input_grid
137
149 subroutine define_input_grid_gaussian(localpet, npets)
150
151#ifdef CHGRES_ALL
152 use nemsio_module
153#endif
154
155 use program_setup, only : data_dir_input_grid, &
156 atm_files_input_grid, &
157 sfc_files_input_grid, &
158 input_type, &
159 convert_atm, convert_sfc
160
161#ifdef CHGRES_ALL
162 use sfcio_module
163 use sigio_module
164#endif
165 use netcdf
166
167 implicit none
168
169 integer, intent(in) :: localpet, npets
170
171 character(len=250) :: the_file
172
173 integer :: i, j, rc, clb(2), cub(2), ncid, id_grid, idum(2)
174#ifdef CHGRES_ALL
175 integer(sfcio_intkind) :: rc2
176 integer(sigio_intkind) :: rc3
177#endif
178
179 real(esmf_kind_r8), allocatable :: latitude(:,:)
180 real(esmf_kind_r8), allocatable :: longitude(:,:)
181 real(esmf_kind_r8), pointer :: lat_src_ptr(:,:)
182 real(esmf_kind_r8), pointer :: lon_src_ptr(:,:)
183 real(esmf_kind_r8), pointer :: lat_corner_src_ptr(:,:)
184 real(esmf_kind_r8), pointer :: lon_corner_src_ptr(:,:)
185 real(esmf_kind_r8) :: deltalon
186 real(esmf_kind_r8), allocatable :: slat(:), wlat(:)
187
188#ifdef CHGRES_ALL
189 type(nemsio_gfile) :: gfile
190#endif
191 type(esmf_polekind_flag) :: polekindflag(2)
192#ifdef CHGRES_ALL
193 type(sfcio_head) :: sfchead
194 type(sigio_head) :: sighead
195#endif
196
197 type(esmf_vm) :: vm
198
199 print*,"- DEFINE INPUT GRID OBJECT FOR GAUSSIAN DATA."
200
201 num_tiles_input_grid = 1
202
203 if (convert_sfc) then
204 the_file=trim(data_dir_input_grid) // "/" // trim(sfc_files_input_grid(1))
205 elseif (convert_atm) then
206 the_file=trim(data_dir_input_grid) // "/" // trim(atm_files_input_grid(1))
207 endif
208
209 if (trim(input_type) == "gaussian_netcdf") then
210
211 call esmf_vmgetglobal(vm, rc=rc)
212 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
213 call error_handler("IN VMGetGlobal", rc)
214
215 if (localpet == 0) then
216 print*,'- OPEN AND READ: ',trim(the_file)
217 rc=nf90_open(trim(the_file),nf90_nowrite,ncid)
218 call netcdf_err(rc, 'opening file')
219
220 print*,"- READ grid_xt"
221 rc=nf90_inq_dimid(ncid, 'grid_xt', id_grid)
222 call netcdf_err(rc, 'reading grid_xt id')
223 rc=nf90_inquire_dimension(ncid,id_grid,len=idum(1))
224 call netcdf_err(rc, 'reading grid_xt')
225
226 print*,"- READ grid_yt"
227 rc=nf90_inq_dimid(ncid, 'grid_yt', id_grid)
228 call netcdf_err(rc, 'reading grid_yt id')
229 rc=nf90_inquire_dimension(ncid,id_grid,len=idum(2))
230 call netcdf_err(rc, 'reading grid_yt')
231
232 rc = nf90_close(ncid)
233 endif
234
235 call esmf_vmbroadcast(vm, idum, 2, 0, rc=rc)
236 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
237 call error_handler("IN VMBroadcast", rc)
238
239 i_input = idum(1)
240 j_input = idum(2)
241
242#ifdef CHGRES_ALL
243 elseif (trim(input_type) == "gfs_sigio") then ! sigio/sfcio format, used by
244
245 if (convert_sfc) then ! sfcio format
246 print*,"- OPEN AND READ ", trim(the_file)
247 call sfcio_sropen(21, trim(the_file), rc2)
248 if (rc2 /= 0) call error_handler("OPENING FILE", rc2)
249 call sfcio_srhead(21, sfchead, rc2)
250 if (rc2 /= 0) call error_handler("READING FILE", rc2)
251 call sfcio_sclose(21, rc2)
252 i_input = sfchead%lonb
253 j_input = sfchead%latb
254 elseif (convert_atm) then ! sigio format
255 print*,"- OPEN AND READ ", trim(the_file)
256 call sigio_sropen(21, trim(the_file), rc3)
257 if (rc3 /= 0) call error_handler("OPENING FILE", rc3)
258 call sigio_srhead(21, sighead, rc3)
259 if (rc3 /= 0) call error_handler("READING FILE", rc3)
260 call sigio_sclose(21, rc3)
261 i_input = sighead%lonb
262 j_input = sighead%latb
263 endif
264
265 else ! nemsio format
266
267 call nemsio_init(iret=rc)
268
269 print*,"- OPEN AND READ ", trim(the_file)
270 call nemsio_open(gfile, the_file, "read", iret=rc)
271 if (rc /= 0) call error_handler("OPENING FILE", rc)
272
273 call nemsio_getfilehead(gfile, iret=rc, dimx=i_input, dimy=j_input)
274 if (rc /= 0) call error_handler("READING FILE", rc)
275
276 call nemsio_close(gfile)
277
278#endif
279
280 endif
281
282 ip1_input = i_input + 1
283 jp1_input = j_input + 1
284
285 polekindflag(1:2) = esmf_polekind_monopole
286
287 print*,"- CALL GridCreate1PeriDim FOR INPUT GRID."
288 input_grid = esmf_gridcreate1peridim(minindex=(/1,1/), &
289 maxindex=(/i_input,j_input/), &
290 polekindflag=polekindflag, &
291 periodicdim=1, &
292 poledim=2, &
293 coordsys=esmf_coordsys_sph_deg, &
294 regdecomp=(/1,npets/), &
295 indexflag=esmf_index_global, rc=rc)
296 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
297 call error_handler("IN GridCreate1PeriDim", rc)
298
299 print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE."
300 latitude_input_grid = esmf_fieldcreate(input_grid, &
301 typekind=esmf_typekind_r8, &
302 staggerloc=esmf_staggerloc_center, &
303 name="input_grid_latitude", rc=rc)
304 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
305 call error_handler("IN FieldCreate", rc)
306
307 print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE."
308 longitude_input_grid = esmf_fieldcreate(input_grid, &
309 typekind=esmf_typekind_r8, &
310 staggerloc=esmf_staggerloc_center, &
311 name="input_grid_longitude", rc=rc)
312 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
313 call error_handler("IN FieldCreate", rc)
314
315 allocate(longitude(i_input,j_input))
316 allocate(latitude(i_input,j_input))
317
318 deltalon = 360.0_esmf_kind_r8 / real(i_input,kind=esmf_kind_r8)
319 do i = 1, i_input
320 longitude(i,:) = real((i-1),kind=esmf_kind_r8) * deltalon
321 enddo
322
323 allocate(slat(j_input))
324 allocate(wlat(j_input))
325 call splat(4, j_input, slat, wlat)
326
327 do i = 1, j_input
328 latitude(:,i) = 90.0_esmf_kind_r8 - (acos(slat(i))* 180.0_esmf_kind_r8 / &
329 (4.0_esmf_kind_r8*atan(1.0_esmf_kind_r8)))
330 enddo
331
332 deallocate(slat, wlat)
333
334 print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE."
335 call esmf_fieldscatter(longitude_input_grid, longitude, rootpet=0, rc=rc)
336 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
337 call error_handler("IN FieldScatter", rc)
338
339 print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE."
340 call esmf_fieldscatter(latitude_input_grid, latitude, rootpet=0, rc=rc)
341 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
342 call error_handler("IN FieldScatter", rc)
343
344 print*,"- CALL GridAddCoord FOR INPUT GRID."
345 call esmf_gridaddcoord(input_grid, &
346 staggerloc=esmf_staggerloc_center, rc=rc)
347 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
348 call error_handler("IN GridAddCoord", rc)
349
350 print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
351 nullify(lon_src_ptr)
352 call esmf_gridgetcoord(input_grid, &
353 staggerloc=esmf_staggerloc_center, &
354 coorddim=1, &
355 farrayptr=lon_src_ptr, rc=rc)
356 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
357 call error_handler("IN GridGetCoord", rc)
358
359 print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
360 nullify(lat_src_ptr)
361 call esmf_gridgetcoord(input_grid, &
362 staggerloc=esmf_staggerloc_center, &
363 coorddim=2, &
364 computationallbound=clb, &
365 computationalubound=cub, &
366 farrayptr=lat_src_ptr, rc=rc)
367 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
368 call error_handler("IN GridGetCoord", rc)
369
370 do j = clb(2), cub(2)
371 do i = clb(1), cub(1)
372 lon_src_ptr(i,j) = longitude(i,j)
373 if (lon_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_src_ptr(i,j) = lon_src_ptr(i,j) - 360.0_esmf_kind_r8
374 lat_src_ptr(i,j) = latitude(i,j)
375 enddo
376 enddo
377
378 print*,"- CALL GridAddCoord FOR INPUT GRID."
379 call esmf_gridaddcoord(input_grid, &
380 staggerloc=esmf_staggerloc_corner, rc=rc)
381 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
382 call error_handler("IN GridAddCoord", rc)
383
384 print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
385 nullify(lon_corner_src_ptr)
386 call esmf_gridgetcoord(input_grid, &
387 staggerloc=esmf_staggerloc_corner, &
388 coorddim=1, &
389 farrayptr=lon_corner_src_ptr, rc=rc)
390 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
391 call error_handler("IN GridGetCoord", rc)
392
393 print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
394 nullify(lat_corner_src_ptr)
395 call esmf_gridgetcoord(input_grid, &
396 staggerloc=esmf_staggerloc_corner, &
397 coorddim=2, &
398 computationallbound=clb, &
399 computationalubound=cub, &
400 farrayptr=lat_corner_src_ptr, rc=rc)
401 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
402 call error_handler("IN GridGetCoord", rc)
403
404 do j = clb(2), cub(2)
405 do i = clb(1), cub(1)
406 lon_corner_src_ptr(i,j) = longitude(i,1) - (0.5_esmf_kind_r8*deltalon)
407 if (lon_corner_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_corner_src_ptr(i,j) = lon_corner_src_ptr(i,j) - 360.0_esmf_kind_r8
408 if (j == 1) then
409 lat_corner_src_ptr(i,j) = 90.0_esmf_kind_r8
410 cycle
411 endif
412 if (j == jp1_input) then
413 lat_corner_src_ptr(i,j) = -90.0_esmf_kind_r8
414 cycle
415 endif
416 lat_corner_src_ptr(i,j) = 0.5_esmf_kind_r8 * (latitude(i,j-1)+ latitude(i,j))
417 enddo
418 enddo
419
420 deallocate(latitude,longitude)
421
422 end subroutine define_input_grid_gaussian
423
430 subroutine define_input_grid_mosaic(localpet, npets)
431
432 use netcdf
433 use program_setup, only : mosaic_file_input_grid, &
434 orog_dir_input_grid, &
435 orog_files_input_grid
436
437 implicit none
438
439 character(len=500) :: the_file
440
441 integer, intent(in) :: localpet, npets
442
443 integer :: id_tiles, id_dim, tile
444 integer :: extra, error, ncid, idum(1), idum2(2)
445 integer, allocatable :: decomptile(:,:)
446
447 real(esmf_kind_r8), allocatable :: latitude_one_tile(:,:)
448 real(esmf_kind_r8), allocatable :: latitude_s_one_tile(:,:)
449 real(esmf_kind_r8), allocatable :: latitude_w_one_tile(:,:)
450 real(esmf_kind_r8), allocatable :: longitude_one_tile(:,:)
451 real(esmf_kind_r8), allocatable :: longitude_s_one_tile(:,:)
452 real(esmf_kind_r8), allocatable :: longitude_w_one_tile(:,:)
453
454 type(esmf_vm) :: vm
455
456 call esmf_vmgetglobal(vm, rc=error)
457 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
458 call error_handler("IN VMGetGlobal", error)
459
460 if (localpet == 0) then
461 print*,'- OPEN INPUT GRID MOSAIC FILE: ',trim(mosaic_file_input_grid)
462 error=nf90_open(trim(mosaic_file_input_grid),nf90_nowrite,ncid)
463 call netcdf_err(error, 'opening grid mosaic file')
464 print*,"- READ NUMBER OF TILES"
465 error=nf90_inq_dimid(ncid, 'ntiles', id_tiles)
466 call netcdf_err(error, 'reading ntiles id')
467 error=nf90_inquire_dimension(ncid,id_tiles,len=idum(1))
468 call netcdf_err(error, 'reading ntiles')
469 error = nf90_close(ncid)
470 endif
471
472 call esmf_vmbroadcast(vm, idum, 1, 0, rc=error)
473 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
474 call error_handler("IN VMBroadcast", error)
475
476 num_tiles_input_grid = idum(1)
477
478 print*,'- NUMBER OF TILES, INPUT MODEL GRID IS ', num_tiles_input_grid
479
480 if (mod(npets,num_tiles_input_grid) /= 0) then
481 call error_handler("MUST RUN WITH A TASK COUNT THAT IS A MULTIPLE OF 6.", 1)
482 endif
483
484!-----------------------------------------------------------------------
485! Create ESMF grid object for the model grid.
486!-----------------------------------------------------------------------
487
488 extra = npets / num_tiles_input_grid
489
490 allocate(decomptile(2,num_tiles_input_grid))
491
492 do tile = 1, num_tiles_input_grid
493 decomptile(:,tile)=(/1,extra/)
494 enddo
495
496 print*,"- CALL GridCreateMosaic FOR INPUT MODEL GRID"
497 input_grid = esmf_gridcreatemosaic(filename=trim(mosaic_file_input_grid), &
498 regdecompptile=decomptile, &
499 staggerloclist=(/esmf_staggerloc_center, esmf_staggerloc_corner, &
500 esmf_staggerloc_edge1, esmf_staggerloc_edge2/), &
501 indexflag=esmf_index_global, &
502 tilefilepath=trim(orog_dir_input_grid), &
503 rc=error)
504 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
505 call error_handler("IN GridCreateMosaic", error)
506
507!-----------------------------------------------------------------------
508! Read the mask and lat/lons.
509!-----------------------------------------------------------------------
510
511 print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE."
512 latitude_input_grid = esmf_fieldcreate(input_grid, &
513 typekind=esmf_typekind_r8, &
514 staggerloc=esmf_staggerloc_center, &
515 name="input_grid_latitude", &
516 rc=error)
517 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
518 call error_handler("IN FieldCreate", error)
519
520 print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE."
521 longitude_input_grid = esmf_fieldcreate(input_grid, &
522 typekind=esmf_typekind_r8, &
523 staggerloc=esmf_staggerloc_center, &
524 name="input_grid_longitude", &
525 rc=error)
526 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
527 call error_handler("IN FieldCreate", error)
528
529 print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE_S."
530 latitude_s_input_grid = esmf_fieldcreate(input_grid, &
531 typekind=esmf_typekind_r8, &
532 staggerloc=esmf_staggerloc_edge2, &
533 name="input_grid_latitude_s", &
534 rc=error)
535 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
536 call error_handler("IN FieldCreate", error)
537
538 print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE_S."
539 longitude_s_input_grid = esmf_fieldcreate(input_grid, &
540 typekind=esmf_typekind_r8, &
541 staggerloc=esmf_staggerloc_edge2, &
542 name="input_grid_longitude_s", &
543 rc=error)
544 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
545 call error_handler("IN FieldCreate", error)
546
547 print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE_W."
548 latitude_w_input_grid = esmf_fieldcreate(input_grid, &
549 typekind=esmf_typekind_r8, &
550 staggerloc=esmf_staggerloc_edge1, &
551 name="input_grid_latitude_w", &
552 rc=error)
553 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
554 call error_handler("IN FieldCreate", error)
555
556 print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE_W."
557 longitude_w_input_grid = esmf_fieldcreate(input_grid, &
558 typekind=esmf_typekind_r8, &
559 staggerloc=esmf_staggerloc_edge1, &
560 name="input_grid_longitude_w", &
561 rc=error)
562 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
563 call error_handler("IN FieldCreate", error)
564
565 the_file = trim(orog_dir_input_grid) // trim(orog_files_input_grid(1))
566
567 if (localpet == 0) then
568 print*,'- OPEN FIRST INPUT GRID OROGRAPHY FILE: ',trim(the_file)
569 error=nf90_open(trim(the_file),nf90_nowrite,ncid)
570 call netcdf_err(error, 'opening ororgraphy file')
571 print*,"- READ GRID DIMENSIONS"
572 error=nf90_inq_dimid(ncid, 'lon', id_dim)
573 call netcdf_err(error, 'reading lon id')
574 error=nf90_inquire_dimension(ncid,id_dim,len=idum2(1))
575 call netcdf_err(error, 'reading lon')
576 error=nf90_inq_dimid(ncid, 'lat', id_dim)
577 call netcdf_err(error, 'reading lat id')
578 error=nf90_inquire_dimension(ncid,id_dim,len=idum2(2))
579 call netcdf_err(error, 'reading lat')
580 error = nf90_close(ncid)
581 endif
582
583 call esmf_vmbroadcast(vm, idum2, 2, 0, rc=error)
584 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
585 call error_handler("IN VMBroadcast", error)
586
587 i_input = idum2(1)
588 j_input = idum2(2)
589
590 print*,"- I/J DIMENSIONS OF THE INPUT GRID TILES ", i_input, j_input
591
592 ip1_input = i_input + 1
593 jp1_input = j_input + 1
594
595 if (localpet == 0) then
596 allocate(longitude_one_tile(i_input,j_input))
597 allocate(longitude_s_one_tile(i_input,jp1_input))
598 allocate(longitude_w_one_tile(ip1_input,j_input))
599 allocate(latitude_one_tile(i_input,j_input))
600 allocate(latitude_s_one_tile(i_input,jp1_input))
601 allocate(latitude_w_one_tile(ip1_input,j_input))
602 else
603 allocate(longitude_one_tile(0,0))
604 allocate(longitude_s_one_tile(0,0))
605 allocate(longitude_w_one_tile(0,0))
606 allocate(latitude_one_tile(0,0))
607 allocate(latitude_s_one_tile(0,0))
608 allocate(latitude_w_one_tile(0,0))
609 endif
610
611 do tile = 1, num_tiles_input_grid
612 if (localpet == 0) then
613 call get_model_latlons(mosaic_file_input_grid, orog_dir_input_grid, num_tiles_input_grid, tile, &
614 i_input, j_input, ip1_input, jp1_input, latitude_one_tile, &
615 latitude_s_one_tile, latitude_w_one_tile, longitude_one_tile, &
616 longitude_s_one_tile, longitude_w_one_tile)
617 endif
618 print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE. TILE IS: ", tile
619 call esmf_fieldscatter(latitude_input_grid, latitude_one_tile, rootpet=0, tile=tile, rc=error)
620 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
621 call error_handler("IN FieldScatter", error)
622 print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE. TILE IS: ", tile
623 call esmf_fieldscatter(longitude_input_grid, longitude_one_tile, rootpet=0, tile=tile, rc=error)
624 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
625 call error_handler("IN FieldScatter", error)
626 print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE_S. TILE IS: ", tile
627 call esmf_fieldscatter(latitude_s_input_grid, latitude_s_one_tile, rootpet=0, tile=tile, rc=error)
628 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
629 call error_handler("IN FieldScatter", error)
630 print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE_S. TILE IS: ", tile
631 call esmf_fieldscatter(longitude_s_input_grid, longitude_s_one_tile, rootpet=0, tile=tile, rc=error)
632 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
633 call error_handler("IN FieldScatter", error)
634 print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE_W. TILE IS: ", tile
635 call esmf_fieldscatter(latitude_w_input_grid, latitude_w_one_tile, rootpet=0, tile=tile, rc=error)
636 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
637 call error_handler("IN FieldScatter", error)
638 print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE_W. TILE IS: ", tile
639 call esmf_fieldscatter(longitude_w_input_grid, longitude_w_one_tile, rootpet=0, tile=tile, rc=error)
640 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
641 call error_handler("IN FieldScatter", error)
642 enddo
643
644 deallocate(longitude_one_tile)
645 deallocate(longitude_s_one_tile)
646 deallocate(longitude_w_one_tile)
647 deallocate(latitude_one_tile)
648 deallocate(latitude_s_one_tile)
649 deallocate(latitude_w_one_tile)
650
651 end subroutine define_input_grid_mosaic
652
660 subroutine define_input_grid_grib2(localpet, npets)
661
662 use grib_mod
663 use gdswzd_mod
664 use program_setup, only : grib2_file_input_grid, data_dir_input_grid
665
666 implicit none
667
668 integer, intent(in) :: localpet, npets
669
670 character(len=500) :: the_file
671
672 integer :: i, j, k, jdisc, jgdtn, jpdtn, lugb, lugi
673 integer :: jids(200), jgdt(200), jpdt(200), rc
674 integer :: kgds(200), nret, clb(2), cub(2), idum(3)
675
676 logical :: unpack
677
678 real :: res
679 real, allocatable :: rlon(:,:),rlat(:,:),xpts(:,:),ypts(:,:)
680 real, allocatable :: rlon_corner(:,:),rlat_corner(:,:)
681 real, allocatable :: xpts_corner(:,:),ypts_corner(:,:)
682 real(esmf_kind_r8), allocatable :: latitude(:,:)
683 real(esmf_kind_r8), allocatable :: longitude(:,:)
684 real(esmf_kind_r8), allocatable :: latitude_corner(:,:)
685 real(esmf_kind_r8), allocatable :: longitude_corner(:,:)
686 real(esmf_kind_r8), pointer :: lat_src_ptr(:,:), lat_ptr(:,:)
687 real(esmf_kind_r8), pointer :: lat_corner_src_ptr(:,:)
688 real(esmf_kind_r8), pointer :: lon_src_ptr(:,:), lon_ptr(:,:)
689 real(esmf_kind_r8), pointer :: lon_corner_src_ptr(:,:)
690
691 type(esmf_field) :: lat_corner, lon_corner
692
693 type(esmf_polekind_flag) :: polekindflag(2)
694
695 type(gribfield) :: gfld
696
697 type(esmf_vm) :: vm
698
699 call esmf_vmgetglobal(vm, rc=rc)
700 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
701 call error_handler("IN VMGetGlobal", rc)
702
703 the_file = trim(data_dir_input_grid) // "/" // grib2_file_input_grid
704
705 lugb=12
706
707 if (localpet == 0) then
708
709 print*,"- OPEN AND READ INPUT DATA GRIB2 FILE: ", trim(the_file)
710 call baopenr(lugb,the_file,rc)
711 if (rc /= 0) call error_handler("OPENING FILE", rc)
712
713! Read the first record and get the grid definition template.
714
715 j = 0 ! Search at beginning of file
716 lugi = 0 ! No grib index file
717 jdisc = -1 ! Search for any discipline
718 jpdtn = -1 ! Search for any product definition template number
719 jgdtn = -1 ! Search for any grid definition template number
720 jids = -9999 ! Array of values in identification section, set to wildcard.
721 jgdt = -9999 ! Array of values in grid definition template, set to wildcard.
722 jpdt = -9999 ! Array of values in product definition template, set to wildcard.
723 unpack = .false. ! unpack data
724
725 call getgb2(lugb, lugi, j, jdisc, jids, jpdtn, jpdt, jgdtn, jgdt, &
726 unpack, k, gfld, rc)
727 if (rc /= 0) call error_handler("DEGRIBBING INPUT FILE.", rc)
728
729 call baclose(lugb,rc)
730
731 kgds = 0
732 call gdt_to_gds(gfld%igdtnum, gfld%igdtmpl, gfld%igdtlen, kgds, idum(1), idum(2), res)
733
734 idum(3) = gfld%igdtnum
735
736 endif
737
738 call esmf_vmbroadcast(vm, idum, 3, 0, rc=rc)
739 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
740 call error_handler("IN VMBroadcast", rc)
741
742 i_input = idum(1)
743 j_input = idum(2)
744
745 ip1_input = i_input + 1
746 jp1_input = j_input + 1
747
748 if (idum(3) == 0) then
749 print*,"- INPUT DATA ON LAT/LON GRID."
750 input_grid_type = 'latlon'
751 elseif (idum(3) == 30) then
752 print*,"- INPUT DATA ON LAMBERT CONFORMAL GRID."
753 input_grid_type = 'lambert'
754 elseif (idum(3) == 32769 .or. idum(3) == 1) then
755 print*,"- INPUT DATA ON ROTATED LAT/LON GRID."
756 input_grid_type = 'rotated_latlon'
757 else
758 call error_handler("INPUT GRID TEMPLATE NOT SUPPORTED.", 2)
759 endif
760
761 if (localpet == 0) then
762
763 allocate(rlat(i_input,j_input))
764 allocate(rlon(i_input,j_input))
765 allocate(xpts(i_input,j_input))
766 allocate(ypts(i_input,j_input))
767 allocate(rlat_corner(ip1_input,jp1_input))
768 allocate(rlon_corner(ip1_input,jp1_input))
769 allocate(xpts_corner(ip1_input,jp1_input))
770 allocate(ypts_corner(ip1_input,jp1_input))
771
772 do j = 1, j_input
773 do i = 1, i_input
774 xpts(i,j) = float(i)
775 ypts(i,j) = float(j)
776 enddo
777 enddo
778
779 print*,"- COMPUTE GRID CELL CENTER COORDINATES."
780 call gdswzd(kgds,1,(i_input*j_input),-9999.,xpts,ypts,rlon,rlat,nret)
781
782 if (nret /= (i_input*j_input)) then
783 call error_handler("GDSWZD RETURNED WRONG NUMBER OF POINTS.", 2)
784 endif
785
786 deallocate(xpts, ypts)
787
788 do j = 1, jp1_input
789 do i = 1, ip1_input
790 xpts_corner(i,j) = float(i) - 0.5
791 ypts_corner(i,j) = float(j) - 0.5
792 enddo
793 enddo
794
795 print*,"- COMPUTE GRID CELL CORNER COORDINATES."
796 call gdswzd(kgds,1,(ip1_input*jp1_input),-9999.,xpts_corner,ypts_corner,rlon_corner,rlat_corner,nret)
797
798 if (nret /= (ip1_input*jp1_input)) then
799 call error_handler("GDSWZD RETURNED WRONG NUMBER OF POINTS.", 2)
800 endif
801
802 deallocate(xpts_corner, ypts_corner)
803
804 end if
805
806 if (trim(input_grid_type) == 'latlon') then ! gfs lat/lon data
807
808 print*,"- CALL GridCreate1PeriDim FOR INPUT GRID."
809
810 polekindflag(1:2) = esmf_polekind_monopole
811
812 input_grid = esmf_gridcreate1peridim(minindex=(/1,1/), &
813 maxindex=(/i_input,j_input/), &
814 polekindflag=polekindflag, &
815 periodicdim=1, &
816 poledim=2, &
817 coordsys=esmf_coordsys_sph_deg, &
818 regdecomp=(/1,npets/), &
819 indexflag=esmf_index_global, rc=rc)
820
821 else
822
823 print*,"- CALL GridCreateNoPeriDim FOR INPUT GRID."
824
825 input_grid = esmf_gridcreatenoperidim(maxindex=(/i_input,j_input/), &
826 indexflag=esmf_index_global, &
827 rc=rc)
828
829 endif
830
831 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
832 call error_handler("IN GridCreate1PeriDim", rc)
833
834 print*,"- CALL FieldCreate FOR INPUT GRID LATITUDE."
835 latitude_input_grid = esmf_fieldcreate(input_grid, &
836 typekind=esmf_typekind_r8, &
837 staggerloc=esmf_staggerloc_center, &
838 name="input_grid_latitude", rc=rc)
839 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
840 call error_handler("IN FieldCreate", rc)
841
842 print*,"- CALL FieldCreate FOR INPUT GRID LONGITUDE."
843 longitude_input_grid = esmf_fieldcreate(input_grid, &
844 typekind=esmf_typekind_r8, &
845 staggerloc=esmf_staggerloc_center, &
846 name="input_grid_longitude", rc=rc)
847 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
848 call error_handler("IN FieldCreate", rc)
849
850 if (localpet == 0) then
851 allocate(latitude(i_input,j_input))
852 allocate(longitude(i_input,j_input))
853 latitude = rlat
854 longitude = rlon
855 deallocate (rlat, rlon)
856 else
857 allocate(latitude(0,0))
858 allocate(longitude(0,0))
859 endif
860
861 print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE."
862 call esmf_fieldscatter(longitude_input_grid, longitude, rootpet=0, rc=rc)
863 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
864 call error_handler("IN FieldScatter", rc)
865
866 print*,"- CALL FieldScatter FOR INPUT GRID LATITUDE."
867 call esmf_fieldscatter(latitude_input_grid, latitude, rootpet=0, rc=rc)
868 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
869 call error_handler("IN FieldScatter", rc)
870
871 deallocate(latitude, longitude)
872
873 print*,"- CALL GridAddCoord FOR INPUT GRID."
874 call esmf_gridaddcoord(input_grid, &
875 staggerloc=esmf_staggerloc_center, rc=rc)
876 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
877 call error_handler("IN GridAddCoord", rc)
878
879 print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
880 nullify(lon_src_ptr)
881 call esmf_gridgetcoord(input_grid, &
882 staggerloc=esmf_staggerloc_center, &
883 coorddim=1, &
884 farrayptr=lon_src_ptr, rc=rc)
885 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
886 call error_handler("IN GridGetCoord", rc)
887
888 print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
889 nullify(lat_src_ptr)
890 call esmf_gridgetcoord(input_grid, &
891 staggerloc=esmf_staggerloc_center, &
892 coorddim=2, &
893 computationallbound=clb, &
894 computationalubound=cub, &
895 farrayptr=lat_src_ptr, rc=rc)
896 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
897 call error_handler("IN GridGetCoord", rc)
898
899 nullify(lat_ptr)
900 print*,"- CALL FieldGet FOR latitude."
901 call esmf_fieldget(latitude_input_grid, &
902 farrayptr=lat_ptr, rc=rc)
903 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
904 call error_handler("IN FieldGet", rc)
905
906 nullify(lon_ptr)
907 print*,"- CALL FieldGet FOR longitude."
908 call esmf_fieldget(longitude_input_grid, &
909 farrayptr=lon_ptr, rc=rc)
910 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
911 call error_handler("IN FieldGet", rc)
912
913 do j = clb(2), cub(2)
914 do i = clb(1), cub(1)
915 lon_src_ptr(i,j) = lon_ptr(i,j)
916 if (lon_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_src_ptr(i,j) = lon_src_ptr(i,j) - 360.0_esmf_kind_r8
917 lat_src_ptr(i,j) = lat_ptr(i,j)
918 enddo
919 enddo
920
921 if (localpet == 0) then
922 if (trim(input_grid_type) == 'latlon') then ! gfs lat/lon data. grid is periodic in 'i'
923 ! direction.
924 allocate(latitude_corner(i_input,jp1_input))
925 allocate(longitude_corner(i_input,jp1_input))
926 latitude_corner(1:i_input,:) = rlat_corner(1:i_input,:)
927 longitude_corner(1:i_input,:) = rlon_corner(1:i_input,:)
928 else
929 allocate(latitude_corner(ip1_input,jp1_input))
930 allocate(longitude_corner(ip1_input,jp1_input))
931 longitude_corner = rlon_corner
932 latitude_corner = rlat_corner
933 endif
934 deallocate (rlat_corner, rlon_corner)
935 else
936 allocate(latitude_corner(0,0))
937 allocate(longitude_corner(0,0))
938 endif
939
940 print*,"- CALL FieldCreate FOR INPUT GRID CORNER LATITUDE."
941 lat_corner = esmf_fieldcreate(input_grid, &
942 typekind=esmf_typekind_r8, &
943 staggerloc=esmf_staggerloc_corner, &
944 name="input_grid_corner_latitude", rc=rc)
945 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
946 call error_handler("IN FieldCreate", rc)
947
948 print*,"- CALL FieldCreate FOR INPUT GRID CORNER LONGITUDE."
949 lon_corner = esmf_fieldcreate(input_grid, &
950 typekind=esmf_typekind_r8, &
951 staggerloc=esmf_staggerloc_corner, &
952 name="input_grid_corner_longitude", rc=rc)
953 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
954 call error_handler("IN FieldCreate", rc)
955
956 print*,"- CALL FieldScatter FOR INPUT GRID CORNER LONGITUDE."
957 call esmf_fieldscatter(lon_corner, longitude_corner, rootpet=0, rc=rc)
958 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
959 call error_handler("IN FieldScatter", rc)
960
961 print*,"- CALL FieldScatter FOR INPUT GRID CORNER LATITUDE."
962 call esmf_fieldscatter(lat_corner, latitude_corner, rootpet=0, rc=rc)
963 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
964 call error_handler("IN FieldScatter", rc)
965
966 deallocate(latitude_corner, longitude_corner)
967
968 print*,"- CALL GridAddCoord FOR INPUT GRID."
969 call esmf_gridaddcoord(input_grid, &
970 staggerloc=esmf_staggerloc_corner, rc=rc)
971 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
972 call error_handler("IN GridAddCoord", rc)
973
974 print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
975 nullify(lon_corner_src_ptr)
976 call esmf_gridgetcoord(input_grid, &
977 staggerloc=esmf_staggerloc_corner, &
978 coorddim=1, &
979 farrayptr=lon_corner_src_ptr, rc=rc)
980 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
981 call error_handler("IN GridGetCoord", rc)
982
983 print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
984 nullify(lat_corner_src_ptr)
985 call esmf_gridgetcoord(input_grid, &
986 staggerloc=esmf_staggerloc_corner, &
987 coorddim=2, &
988 computationallbound=clb, &
989 computationalubound=cub, &
990 farrayptr=lat_corner_src_ptr, rc=rc)
991 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
992 call error_handler("IN GridGetCoord", rc)
993
994 nullify(lat_ptr)
995 print*,"- CALL FieldGet FOR corner latitude."
996 call esmf_fieldget(lat_corner, &
997 farrayptr=lat_ptr, rc=rc)
998 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
999 call error_handler("IN FieldGet", rc)
1000
1001 nullify(lon_ptr)
1002 print*,"- CALL FieldGet FOR corner longitude."
1003 call esmf_fieldget(lon_corner, &
1004 farrayptr=lon_ptr, rc=rc)
1005 if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1006 call error_handler("IN FieldGet", rc)
1007
1008 do j = clb(2), cub(2)
1009 do i = clb(1), cub(1)
1010 lon_corner_src_ptr(i,j) = lon_ptr(i,j)
1011 if (lon_corner_src_ptr(i,j) > 360.0_esmf_kind_r8) lon_corner_src_ptr(i,j) = lon_corner_src_ptr(i,j) - 360.0_esmf_kind_r8
1012 lat_corner_src_ptr(i,j) = lat_ptr(i,j)
1013 enddo
1014 enddo
1015
1016 call esmf_fielddestroy(lon_corner, rc=rc)
1017 call esmf_fielddestroy(lat_corner, rc=rc)
1018
1019 end subroutine define_input_grid_grib2
1020
1026 subroutine define_target_grid(localpet, npets)
1027
1028 use netcdf
1029 use program_setup, only : mosaic_file_target_grid, &
1030 orog_dir_target_grid, &
1031 orog_files_target_grid, &
1032 nsoill_out
1033
1034 implicit none
1035
1036 integer, intent(in) :: localpet, npets
1037
1038 character(len=500) :: the_file
1039
1040 integer :: error, ncid, extra
1041 integer :: id_tiles
1042 integer :: id_dim, id_grid_tiles
1043 integer :: tile, idum1(1), idum2(2)
1044 integer, allocatable :: decomptile(:,:)
1045 integer(esmf_kind_i8), allocatable :: landmask_one_tile(:,:)
1046 integer(esmf_kind_i8), allocatable :: seamask_one_tile(:,:)
1047
1048 real(esmf_kind_r8), allocatable :: land_frac_one_tile(:,:)
1049 real(esmf_kind_r8), allocatable :: latitude_one_tile(:,:)
1050 real(esmf_kind_r8), allocatable :: latitude_s_one_tile(:,:)
1051 real(esmf_kind_r8), allocatable :: latitude_w_one_tile(:,:)
1052 real(esmf_kind_r8), allocatable :: longitude_one_tile(:,:)
1053 real(esmf_kind_r8), allocatable :: longitude_s_one_tile(:,:)
1054 real(esmf_kind_r8), allocatable :: longitude_w_one_tile(:,:)
1055 real(esmf_kind_r8), allocatable :: terrain_one_tile(:,:)
1056
1057 type(esmf_vm) :: vm
1058
1059 lsoil_target = nsoill_out
1060
1061 call esmf_vmgetglobal(vm, rc=error)
1062 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1063 call error_handler("IN VMGetGlobal", error)
1064
1065 if (localpet == 0) then
1066 print*,'- OPEN TARGET GRID MOSAIC FILE: ',trim(mosaic_file_target_grid)
1067 error=nf90_open(trim(mosaic_file_target_grid),nf90_nowrite,ncid)
1068 call netcdf_err(error, 'opening grid mosaic file')
1069 print*,"- READ NUMBER OF TILES"
1070 error=nf90_inq_dimid(ncid, 'ntiles', id_tiles)
1071 call netcdf_err(error, 'reading ntile id')
1072 error=nf90_inquire_dimension(ncid,id_tiles,len=idum1(1))
1073 call netcdf_err(error, 'reading ntiles')
1074 endif
1075
1076 call esmf_vmbroadcast(vm, idum1, 1, 0, rc=error)
1077 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1078 call error_handler("IN VMBroadcast", error)
1079
1080 num_tiles_target_grid = idum1(1)
1081
1082 allocate(tiles_target_grid(num_tiles_target_grid))
1083 tiles_target_grid="NULL"
1084
1085 if (localpet == 0) then
1086 error=nf90_inq_varid(ncid, 'gridtiles', id_grid_tiles)
1087 call netcdf_err(error, 'reading gridtiles id')
1088 print*,"- READ TILE NAMES"
1089 error=nf90_get_var(ncid, id_grid_tiles, tiles_target_grid)
1090 call netcdf_err(error, 'reading gridtiles')
1091 error = nf90_close(ncid)
1092 endif
1093
1094 call esmf_vmbroadcast(vm, tiles_target_grid, num_tiles_target_grid, 0, rc=error)
1095 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1096 call error_handler("IN VMBroadcast", error)
1097
1098 print*,'- NUMBER OF TILES, TARGET MODEL GRID IS ', num_tiles_target_grid
1099
1100 if (mod(npets,num_tiles_target_grid) /= 0) then
1101 call error_handler("MUST RUN WITH TASK COUNT THAT IS A MULTIPLE OF # OF TILES.", 1)
1102 endif
1103
1104!-----------------------------------------------------------------------
1105! Get the model grid specs and land mask from the orography files.
1106!-----------------------------------------------------------------------
1107
1108 the_file = trim(orog_dir_target_grid) // trim(orog_files_target_grid(1))
1109
1110 if (localpet == 0) then
1111 print*,'- OPEN FIRST TARGET GRID OROGRAPHY FILE: ',trim(the_file)
1112 error=nf90_open(trim(the_file),nf90_nowrite,ncid)
1113 call netcdf_err(error, 'opening orography file')
1114 print*,"- READ GRID DIMENSIONS"
1115 error=nf90_inq_dimid(ncid, 'lon', id_dim)
1116 call netcdf_err(error, 'reading lon id')
1117 error=nf90_inquire_dimension(ncid,id_dim,len=idum2(1))
1118 call netcdf_err(error, 'reading lon')
1119 error=nf90_inq_dimid(ncid, 'lat', id_dim)
1120 call netcdf_err(error, 'reading lat id')
1121 error=nf90_inquire_dimension(ncid,id_dim,len=idum2(2))
1122 call netcdf_err(error, 'reading lat')
1123 error = nf90_close(ncid)
1124 endif
1125
1126 call esmf_vmbroadcast(vm, idum2, 2, 0, rc=error)
1127 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1128 call error_handler("IN VMBroadcast", error)
1129
1130 i_target = idum2(1)
1131 j_target = idum2(2)
1132
1133 print*,"- I/J DIMENSIONS OF THE TARGET GRID TILES ", i_target, j_target
1134
1135 ip1_target = i_target + 1
1136 jp1_target = j_target + 1
1137
1138!-----------------------------------------------------------------------
1139! Create ESMF grid object for the model grid.
1140!-----------------------------------------------------------------------
1141
1142 extra = npets / num_tiles_target_grid
1143
1144 allocate(decomptile(2,num_tiles_target_grid))
1145
1146 do tile = 1, num_tiles_target_grid
1147 decomptile(:,tile)=(/1,extra/)
1148 enddo
1149
1150 print*,"- CALL GridCreateMosaic FOR TARGET GRID"
1151 target_grid = esmf_gridcreatemosaic(filename=trim(mosaic_file_target_grid), &
1152 regdecompptile=decomptile, &
1153 staggerloclist=(/esmf_staggerloc_center, esmf_staggerloc_corner, &
1154 esmf_staggerloc_edge1, esmf_staggerloc_edge2/), &
1155 indexflag=esmf_index_global, &
1156 tilefilepath=trim(orog_dir_target_grid), rc=error)
1157 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1158 call error_handler("IN GridCreateMosaic", error)
1159
1160!-----------------------------------------------------------------------
1161! Set target model landmask (1 - land, 0 - not land) and
1162! seamask (1 - non-land, 0 -land). Read lat/lon on target grid.
1163!-----------------------------------------------------------------------
1164
1165 print*,"- CALL FieldCreate FOR TARGET GRID LAND FRACTION."
1166 land_frac_target_grid = esmf_fieldcreate(target_grid, &
1167 typekind=esmf_typekind_r8, &
1168 staggerloc=esmf_staggerloc_center, &
1169 name="target_grid_landmask", rc=error)
1170 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1171 call error_handler("IN FieldCreate", error)
1172
1173 print*,"- CALL FieldCreate FOR TARGET GRID LANDMASK."
1174 landmask_target_grid = esmf_fieldcreate(target_grid, &
1175 typekind=esmf_typekind_i8, &
1176 staggerloc=esmf_staggerloc_center, &
1177 name="target_grid_landmask", rc=error)
1178 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1179 call error_handler("IN FieldCreate", error)
1180
1181 print*,"- CALL FieldCreate FOR TARGET GRID SEAMASK."
1182 seamask_target_grid = esmf_fieldcreate(target_grid, &
1183 typekind=esmf_typekind_i8, &
1184 staggerloc=esmf_staggerloc_center, &
1185 name="target_grid_seamask", rc=error)
1186 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1187 call error_handler("IN FieldCreate", error)
1188
1189 print*,"- CALL FieldCreate FOR TARGET GRID LATITUDE."
1190 latitude_target_grid = esmf_fieldcreate(target_grid, &
1191 typekind=esmf_typekind_r8, &
1192 staggerloc=esmf_staggerloc_center, &
1193 name="target_grid_latitude", rc=error)
1194 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1195 call error_handler("IN FieldCreate", error)
1196
1197 print*,"- CALL FieldCreate FOR TARGET GRID LATITUDE_S."
1198 latitude_s_target_grid = esmf_fieldcreate(target_grid, &
1199 typekind=esmf_typekind_r8, &
1200 staggerloc=esmf_staggerloc_edge2, &
1201 name="target_grid_latitude_s", rc=error)
1202 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1203 call error_handler("IN FieldCreate", error)
1204
1205 print*,"- CALL FieldCreate FOR TARGET GRID LATITUDE_W."
1206 latitude_w_target_grid = esmf_fieldcreate(target_grid, &
1207 typekind=esmf_typekind_r8, &
1208 staggerloc=esmf_staggerloc_edge1, &
1209 name="target_grid_latitude_w", &
1210 rc=error)
1211 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1212 call error_handler("IN FieldCreate", error)
1213
1214 print*,"- CALL FieldCreate FOR TARGET GRID LONGITUDE."
1215 longitude_target_grid = esmf_fieldcreate(target_grid, &
1216 typekind=esmf_typekind_r8, &
1217 staggerloc=esmf_staggerloc_center, &
1218 name="target_grid_longitude", &
1219 rc=error)
1220 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1221 call error_handler("IN FieldCreate", error)
1222
1223 print*,"- CALL FieldCreate FOR TARGET GRID LONGITUDE_S."
1224 longitude_s_target_grid = esmf_fieldcreate(target_grid, &
1225 typekind=esmf_typekind_r8, &
1226 staggerloc=esmf_staggerloc_edge2, &
1227 name="target_grid_longitude_s", &
1228 rc=error)
1229 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1230 call error_handler("IN FieldCreate", error)
1231
1232 print*,"- CALL FieldCreate FOR TARGET GRID LONGITUDE_W."
1233 longitude_w_target_grid = esmf_fieldcreate(target_grid, &
1234 typekind=esmf_typekind_r8, &
1235 staggerloc=esmf_staggerloc_edge1, &
1236 name="target_grid_longitude_w", &
1237 rc=error)
1238 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1239 call error_handler("IN FieldCreate", error)
1240
1241 print*,"- CALL FieldCreate FOR TARGET GRID TERRAIN."
1242 terrain_target_grid = esmf_fieldcreate(target_grid, &
1243 typekind=esmf_typekind_r8, &
1244 staggerloc=esmf_staggerloc_center, &
1245 name="target_grid_terrain", &
1246 rc=error)
1247 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1248 call error_handler("IN FieldCreate", error)
1249
1250 if (localpet == 0) then
1251 allocate(land_frac_one_tile(i_target,j_target))
1252 allocate(landmask_one_tile(i_target,j_target))
1253 allocate(seamask_one_tile(i_target,j_target))
1254 allocate(latitude_one_tile(i_target,j_target))
1255 allocate(latitude_s_one_tile(i_target,jp1_target))
1256 allocate(latitude_w_one_tile(ip1_target,j_target))
1257 allocate(longitude_one_tile(i_target,j_target))
1258 allocate(longitude_s_one_tile(i_target,jp1_target))
1259 allocate(longitude_w_one_tile(ip1_target,j_target))
1260 allocate(terrain_one_tile(i_target,j_target))
1261 else
1262 allocate(land_frac_one_tile(0,0))
1263 allocate(landmask_one_tile(0,0))
1264 allocate(seamask_one_tile(0,0))
1265 allocate(longitude_one_tile(0,0))
1266 allocate(longitude_s_one_tile(0,0))
1267 allocate(longitude_w_one_tile(0,0))
1268 allocate(latitude_one_tile(0,0))
1269 allocate(latitude_s_one_tile(0,0))
1270 allocate(latitude_w_one_tile(0,0))
1271 allocate(terrain_one_tile(0,0))
1272 endif
1273
1274 do tile = 1, num_tiles_target_grid
1275 if (localpet == 0) then
1276 the_file = trim(orog_dir_target_grid) // trim(orog_files_target_grid(tile))
1277 call get_model_mask_terrain(trim(the_file), i_target, j_target, landmask_one_tile, &
1278 terrain_one_tile, land_frac_one_tile)
1279
1280 seamask_one_tile = 0 ! all land
1281 where(floor(land_frac_one_tile) == 0) seamask_one_tile = 1 ! at least some non-land.
1282 landmask_one_tile = 0 ! all non-land
1283 where(ceiling(land_frac_one_tile) == 1) landmask_one_tile = 1 ! at least some land
1284
1285 call get_model_latlons(mosaic_file_target_grid, orog_dir_target_grid, num_tiles_target_grid, tile, &
1286 i_target, j_target, ip1_target, jp1_target, latitude_one_tile, &
1287 latitude_s_one_tile, latitude_w_one_tile, longitude_one_tile, &
1288 longitude_s_one_tile, longitude_w_one_tile)
1289 endif
1290 print*,"- CALL FieldScatter FOR TARGET GRID LAND FRACTION. TILE IS: ", tile
1291 call esmf_fieldscatter(land_frac_target_grid, land_frac_one_tile, rootpet=0, tile=tile, rc=error)
1292 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1293 call error_handler("IN FieldScatter", error)
1294 print*,"- CALL FieldScatter FOR TARGET GRID LANDMASK. TILE IS: ", tile
1295 call esmf_fieldscatter(landmask_target_grid, landmask_one_tile, rootpet=0, tile=tile, rc=error)
1296 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1297 call error_handler("IN FieldScatter", error)
1298 print*,"- CALL FieldScatter FOR TARGET GRID SEAMASK. TILE IS: ", tile
1299 call esmf_fieldscatter(seamask_target_grid, seamask_one_tile, rootpet=0, tile=tile, rc=error)
1300 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1301 call error_handler("IN FieldScatter", error)
1302 print*,"- CALL FieldScatter FOR TARGET GRID LONGITUDE. TILE IS: ", tile
1303 call esmf_fieldscatter(longitude_target_grid, longitude_one_tile, rootpet=0, tile=tile, rc=error)
1304 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1305 call error_handler("IN FieldScatter", error)
1306 print*,"- CALL FieldScatter FOR TARGET GRID LONGITUDE_S. TILE IS: ", tile
1307 call esmf_fieldscatter(longitude_s_target_grid, longitude_s_one_tile, rootpet=0, tile=tile, rc=error)
1308 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1309 call error_handler("IN FieldScatter", error)
1310 print*,"- CALL FieldScatter FOR TARGET GRID LONGITUDE_W. TILE IS: ", tile
1311 call esmf_fieldscatter(longitude_w_target_grid, longitude_w_one_tile, rootpet=0, tile=tile, rc=error)
1312 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1313 call error_handler("IN FieldScatter", error)
1314 print*,"- CALL FieldScatter FOR TARGET GRID LATITUDE. TILE IS: ", tile
1315 call esmf_fieldscatter(latitude_target_grid, latitude_one_tile, rootpet=0, tile=tile, rc=error)
1316 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1317 call error_handler("IN FieldScatter", error)
1318 print*,"- CALL FieldScatter FOR TARGET GRID LATITUDE_S. TILE IS: ", tile
1319 call esmf_fieldscatter(latitude_s_target_grid, latitude_s_one_tile, rootpet=0, tile=tile, rc=error)
1320 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1321 call error_handler("IN FieldScatter", error)
1322 print*,"- CALL FieldScatter FOR TARGET GRID LATITUDE_W. TILE IS: ", tile
1323 call esmf_fieldscatter(latitude_w_target_grid, latitude_w_one_tile, rootpet=0, tile=tile, rc=error)
1324 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1325 call error_handler("IN FieldScatter", error)
1326 print*,"- CALL FieldScatter FOR TARGET GRID TERRAIN. TILE IS: ", tile
1327 call esmf_fieldscatter(terrain_target_grid, terrain_one_tile, rootpet=0, tile=tile, rc=error)
1328 if(esmf_logfounderror(rctocheck=error,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
1329 call error_handler("IN FieldScatter", error)
1330 enddo
1331
1332 deallocate(land_frac_one_tile)
1333 deallocate(landmask_one_tile)
1334 deallocate(seamask_one_tile)
1335 deallocate(longitude_one_tile)
1336 deallocate(longitude_s_one_tile)
1337 deallocate(longitude_w_one_tile)
1338 deallocate(latitude_one_tile)
1339 deallocate(latitude_s_one_tile)
1340 deallocate(latitude_w_one_tile)
1341 deallocate(terrain_one_tile)
1342
1343 end subroutine define_target_grid
1344
1363 subroutine get_model_latlons(mosaic_file, orog_dir, num_tiles, tile, &
1364 i_tile, j_tile, ip1_tile, jp1_tile, &
1365 latitude, latitude_s, latitude_w, &
1366 longitude, longitude_s, longitude_w)
1367
1368 use netcdf
1369
1370 implicit none
1371
1372 character(len=*), intent(in) :: mosaic_file, orog_dir
1373
1374 integer, intent(in) :: num_tiles, tile
1375 integer, intent(in) :: i_tile, j_tile
1376 integer, intent(in) :: ip1_tile, jp1_tile
1377
1378 real(esmf_kind_r8), intent(out) :: latitude(i_tile, j_tile)
1379 real(esmf_kind_r8), intent(out) :: latitude_s(i_tile, jp1_tile)
1380 real(esmf_kind_r8), intent(out) :: latitude_w(ip1_tile, j_tile)
1381 real(esmf_kind_r8), intent(out) :: longitude(i_tile, j_tile)
1382 real(esmf_kind_r8), intent(out) :: longitude_s(i_tile, jp1_tile)
1383 real(esmf_kind_r8), intent(out) :: longitude_w(ip1_tile, j_tile)
1384
1385 character(len=50) :: grid_files(num_tiles)
1386 character(len=255) :: grid_file
1387
1388 integer :: error, id_var, ncid
1389 integer :: id_dim, nxp, nyp, i, j, ii, jj
1390
1391 real(esmf_kind_r8), allocatable :: tmpvar(:,:)
1392
1393 print*,"- READ MODEL GRID FILE"
1394
1395 print*,'- OPEN MOSAIC FILE: ', trim(mosaic_file)
1396 error=nf90_open(trim(mosaic_file), nf90_nowrite, ncid)
1397 call netcdf_err(error, 'opening mosaic file')
1398
1399 print*,"- READ GRID FILE NAMES"
1400 error=nf90_inq_varid(ncid, 'gridfiles', id_var)
1401 call netcdf_err(error, 'reading gridfiles id')
1402 error=nf90_get_var(ncid, id_var, grid_files)
1403 call netcdf_err(error, 'reading gridfiles')
1404
1405 error = nf90_close(ncid)
1406
1407 grid_file = trim(orog_dir) // trim(grid_files(tile))
1408
1409 print*,'- OPEN GRID FILE: ', trim(grid_file)
1410 error=nf90_open(trim(grid_file), nf90_nowrite, ncid)
1411 call netcdf_err(error, 'opening grid file')
1412
1413 print*,'- READ NXP ID'
1414 error=nf90_inq_dimid(ncid, 'nxp', id_dim)
1415 call netcdf_err(error, 'reading nxp id')
1416
1417 print*,'- READ NXP'
1418 error=nf90_inquire_dimension(ncid,id_dim,len=nxp)
1419 call netcdf_err(error, 'reading nxp')
1420
1421 print*,'- READ NYP ID'
1422 error=nf90_inq_dimid(ncid, 'nyp', id_dim)
1423 call netcdf_err(error, 'reading nyp id')
1424
1425 print*,'- READ NYP'
1426 error=nf90_inquire_dimension(ncid,id_dim,len=nyp)
1427 call netcdf_err(error, 'reading nyp')
1428
1429 if ((nxp/2 /= i_tile) .or. (nyp/2 /= j_tile)) then
1430 call error_handler("DIMENSION MISMATCH IN GRID FILE.", 1)
1431 endif
1432
1433 allocate(tmpvar(nxp,nyp))
1434
1435 print*,'- READ LONGITUDE ID'
1436 error=nf90_inq_varid(ncid, 'x', id_var)
1437 call netcdf_err(error, 'reading longitude id')
1438
1439 print*,'- READ LONGITUDE'
1440 error=nf90_get_var(ncid, id_var, tmpvar)
1441 call netcdf_err(error, 'reading longitude')
1442
1443 do j = 1, j_tile
1444 do i = 1, i_tile
1445 ii = 2*i
1446 jj = 2*j
1447 longitude(i,j) = tmpvar(ii,jj)
1448 enddo
1449 enddo
1450
1451 do j = 1, jp1_tile
1452 do i = 1, i_tile
1453 ii = 2*i
1454 jj = (2*j) - 1
1455 longitude_s(i,j) = tmpvar(ii,jj)
1456 enddo
1457 enddo
1458
1459 do j = 1, j_tile
1460 do i = 1, ip1_tile
1461 ii = (2*i) - 1
1462 jj = 2*j
1463 longitude_w(i,j) = tmpvar(ii,jj)
1464 enddo
1465 enddo
1466
1467 print*,'- READ LATITUDE ID'
1468 error=nf90_inq_varid(ncid, 'y', id_var)
1469 call netcdf_err(error, 'reading latitude id')
1470
1471 print*,'- READ LATIITUDE'
1472 error=nf90_get_var(ncid, id_var, tmpvar)
1473 call netcdf_err(error, 'reading latitude')
1474
1475 do j = 1, j_tile
1476 do i = 1, i_tile
1477 ii = 2*i
1478 jj = 2*j
1479 latitude(i,j) = tmpvar(ii,jj)
1480 enddo
1481 enddo
1482
1483 do j = 1, jp1_tile
1484 do i = 1, i_tile
1485 ii = 2*i
1486 jj = (2*j) - 1
1487 latitude_s(i,j) = tmpvar(ii,jj)
1488 enddo
1489 enddo
1490
1491 do j = 1, j_tile
1492 do i = 1, ip1_tile
1493 ii = (2*i) - 1
1494 jj = 2*j
1495 latitude_w(i,j) = tmpvar(ii,jj)
1496 enddo
1497 enddo
1498
1499 deallocate(tmpvar)
1500
1501 error = nf90_close(ncid)
1502
1503 end subroutine get_model_latlons
1504
1515 subroutine get_model_mask_terrain(orog_file, idim, jdim, mask, terrain, land_frac)
1516
1517 use netcdf
1518
1519 implicit none
1520
1521 character(len=*), intent(in) :: orog_file
1522
1523 integer, intent(in) :: idim, jdim
1524 integer(esmf_kind_i8), intent(out) :: mask(idim,jdim)
1525
1526 real(esmf_kind_i8), intent(out) :: terrain(idim,jdim)
1527 real(esmf_kind_i8), intent(out) :: land_frac(idim,jdim)
1528
1529 integer :: error, lat, lon
1530 integer :: ncid, id_dim, id_var
1531
1532 real(kind=4), allocatable :: dummy(:,:)
1533
1534 print*,"- READ MODEL LAND MASK FILE"
1535
1536 print*,'- OPEN LAND MASK FILE: ', orog_file
1537 error=nf90_open(orog_file,nf90_nowrite,ncid)
1538 call netcdf_err(error, 'opening land mask file')
1539
1540 print*,"- READ I-DIMENSION"
1541 error=nf90_inq_dimid(ncid, 'lon', id_dim)
1542 call netcdf_err(error, 'reading idim id')
1543 error=nf90_inquire_dimension(ncid,id_dim,len=lon)
1544 call netcdf_err(error, 'reading idim')
1545
1546 print*,"- READ J-DIMENSION"
1547 error=nf90_inq_dimid(ncid, 'lat', id_dim)
1548 call netcdf_err(error, 'reading jdim id')
1549 error=nf90_inquire_dimension(ncid,id_dim,len=lat)
1550 call netcdf_err(error, 'reading jdim')
1551
1552 print*,"- I/J DIMENSIONS: ", lon, lat
1553
1554 if ((lon /= idim) .or. (lat /= jdim)) then
1555 call error_handler("MISMATCH IN DIMENSIONS.", 1)
1556 endif
1557
1558 allocate(dummy(idim,jdim))
1559
1560 print*,"- READ LAND MASK"
1561 error=nf90_inq_varid(ncid, 'slmsk', id_var)
1562 call netcdf_err(error, 'reading slmsk id')
1563 error=nf90_get_var(ncid, id_var, dummy)
1564 call netcdf_err(error, 'reading slmsk')
1565 mask = nint(dummy)
1566
1567 print*,"- READ RAW OROGRAPHY."
1568 error=nf90_inq_varid(ncid, 'orog_raw', id_var)
1569 call netcdf_err(error, 'reading orog_raw id')
1570 error=nf90_get_var(ncid, id_var, dummy)
1571 call netcdf_err(error, 'reading orog_raw')
1572 terrain = dummy
1573
1574 print*,"- READ LAND FRACTION."
1575 error=nf90_inq_varid(ncid, 'land_frac', id_var)
1576 call netcdf_err(error, 'reading land_frac id')
1577 error=nf90_get_var(ncid, id_var, dummy)
1578 call netcdf_err(error, 'reading orog_raw')
1579 land_frac = dummy
1580
1581 error = nf90_close(ncid)
1582
1583 deallocate (dummy)
1584
1585 end subroutine get_model_mask_terrain
1586
1590 subroutine cleanup_input_target_grid_data
1591
1592 implicit none
1593
1594 integer :: rc
1595
1596 print*,"- DESTROY MODEL DATA."
1597
1598 call esmf_fielddestroy(latitude_input_grid,rc=rc)
1599 call esmf_fielddestroy(longitude_input_grid,rc=rc)
1600 if (esmf_fieldiscreated(latitude_s_input_grid)) then
1601 call esmf_fielddestroy(latitude_s_input_grid, rc=rc)
1602 endif
1603 if (esmf_fieldiscreated(latitude_w_input_grid)) then
1604 call esmf_fielddestroy(latitude_w_input_grid, rc=rc)
1605 endif
1606 if (esmf_fieldiscreated(longitude_s_input_grid)) then
1607 call esmf_fielddestroy(longitude_s_input_grid, rc=rc)
1608 endif
1609 if (esmf_fieldiscreated(longitude_w_input_grid)) then
1610 call esmf_fielddestroy(longitude_w_input_grid, rc=rc)
1611 endif
1612 call esmf_fielddestroy(land_frac_target_grid, rc=rc)
1613 call esmf_fielddestroy(landmask_target_grid, rc=rc)
1614 call esmf_fielddestroy(latitude_target_grid, rc=rc)
1615 if (esmf_fieldiscreated(latitude_s_target_grid)) then
1616 call esmf_fielddestroy(latitude_s_target_grid, rc=rc)
1617 endif
1618 if (esmf_fieldiscreated(latitude_w_target_grid)) then
1619 call esmf_fielddestroy(latitude_w_target_grid, rc=rc)
1620 endif
1621 call esmf_fielddestroy(longitude_target_grid, rc=rc)
1622 if (esmf_fieldiscreated(longitude_s_target_grid)) then
1623 call esmf_fielddestroy(longitude_s_target_grid, rc=rc)
1624 endif
1625 if (esmf_fieldiscreated(longitude_w_target_grid)) then
1626 call esmf_fielddestroy(longitude_w_target_grid, rc=rc)
1627 endif
1628 call esmf_fielddestroy(seamask_target_grid, rc=rc)
1629 call esmf_fielddestroy(terrain_target_grid, rc=rc)
1630 call esmf_griddestroy(input_grid, rc=rc)
1631 call esmf_griddestroy(target_grid, rc=rc)
1632
1633 end subroutine cleanup_input_target_grid_data
1634
1646 subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)
1647
1648 implicit none
1649
1650 integer, intent(in ) :: igdtnum, igdtlen, igdstmpl(igdtlen)
1651 integer, intent( out) :: kgds(200), ni, nj
1652 integer :: iscale, i
1653
1654 real, intent( out) :: res
1655
1656 real :: clatr, slatr, clonr, dpr, slat
1657 real :: slat0, clat0, clat, clon, rlat
1658 real :: rlon0, rlon, hs
1659
1660 kgds=0
1661
1662 if (igdtnum.eq.32769) then ! rot lat/lon b grid
1663
1664 iscale=igdstmpl(10)*igdstmpl(11)
1665 if (iscale == 0) iscale = 1e6
1666 kgds(1)=205 ! oct 6, rotated lat/lon for Non-E
1667 ! Stagger grid
1668 kgds(2)=igdstmpl(8) ! octs 7-8, Ni
1669 ni = kgds(2)
1670 kgds(3)=igdstmpl(9) ! octs 9-10, Nj
1671 nj = kgds(3)
1672 kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of
1673 ! 1st grid point
1674 kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of
1675 ! 1st grid point
1676
1677 kgds(6)=0 ! oct 17, resolution and component flags
1678 if (igdstmpl(1)==2 ) kgds(6)=64
1679 if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
1680 if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
1681
1682 kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 18-20,
1683 ! Lat of cent of rotation
1684 kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 21-23,
1685 ! Lon of cent of rotation
1686 kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25,
1687 ! Di
1688 kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27,
1689 ! Dj
1690
1691 kgds(11) = 0 ! oct 28, scan mode
1692 if (btest(igdstmpl(19),7)) kgds(11) = 128
1693 if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64
1694 if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32
1695
1696 kgds(12)=nint(float(igdstmpl(20))/float(iscale)*1000.) ! octs 29-31, Lat of
1697 ! last grid point
1698 kgds(13)=nint(float(igdstmpl(21))/float(iscale)*1000.) ! octs 32-34, Lon of
1699 ! last grid point
1700
1701 kgds(19)=0 ! oct 4, # vert coordinate parameters
1702 kgds(20)=255 ! oct 5, used for thinned grids, set to 255
1703
1704 res = ((float(kgds(9)) / 1000.0) + (float(kgds(10)) / 1000.0)) &
1705 * 0.5 * 111.0
1706
1707 elseif (igdtnum.eq.1) then ! rot lat/lon b grid using standard wmo
1708 ! template.
1709
1710 iscale=igdstmpl(10)*igdstmpl(11)
1711 if (iscale == 0) iscale = 1e6
1712 kgds(1)=205 ! oct 6, rotated lat/lon for Non-E
1713 ! Stagger grid
1714 kgds(2)=igdstmpl(8) ! octs 7-8, Ni
1715 ni = kgds(2)
1716 kgds(3)=igdstmpl(9) ! octs 9-10, Nj
1717 nj = kgds(3)
1718
1719 kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of
1720 ! 1st grid point
1721 kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of
1722 ! 1st grid point
1723
1724 kgds(6)=0 ! oct 17, resolution and component flags
1725 if (igdstmpl(1)==2 ) kgds(6)=64
1726 if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
1727 if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
1728
1729 kgds(7)=nint(float(igdstmpl(20))/float(iscale)*1000.) ! octs 18-20,
1730 ! Lat of cent of rotation
1731 kgds(8)=nint(float(igdstmpl(21))/float(iscale)*1000.) ! octs 21-23,
1732 ! Lon of cent of rotation
1733 kgds(7) = kgds(7) + 90000.0
1734 print*, "INPUT LAT, LON CENTER ", kgds(7), kgds(8)
1735
1736 dpr = 180.0/3.1415926
1737 clatr=cos((float(kgds(4))/1000.0)/dpr)
1738 slatr=sin((float(kgds(4))/1000.0)/dpr)
1739 clonr=cos((float(kgds(5))/1000.0)/dpr)
1740 slat0=sin((float(kgds(7))/1000.0)/dpr)
1741 clat0=cos((float(kgds(7))/1000.0)/dpr)
1742
1743 slat=clat0*slatr+slat0*clatr*clonr
1744 clat=sqrt(1-slat**2)
1745 clon=(clat0*clatr*clonr-slat0*slatr)/clat
1746 clon=min(max(clon,-1.0),1.0)
1747
1748 rlat=dpr*asin(slat)
1749
1750 rlon0=float(kgds(8))/1000.0
1751
1752 if ((kgds(5)-kgds(8)) > 0) then
1753 hs = -1.0
1754 else
1755 hs = 1.0
1756 endif
1757
1758 rlon = mod(rlon0+hs*dpr*acos(clon)+3600,360.0)
1759
1760 kgds(4)=nint(rlat*1000.) ! octs 11-13, Lat of
1761 kgds(5)=nint(rlon*1000.) ! octs 14-16, Lon of
1762
1763 kgds(12)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 29-31, Lat of
1764 ! last grid point
1765 kgds(13)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 32-34, Lon of
1766 ! last grid point
1767
1768 clatr=cos((float(kgds(12))/1000.0)/dpr)
1769 slatr=sin((float(kgds(12))/1000.0)/dpr)
1770 clonr=cos((float(kgds(13))/1000.0)/dpr)
1771
1772 slat=clat0*slatr+slat0*clatr*clonr
1773 rlat=dpr*asin(slat)
1774
1775 clat=sqrt(1-slat**2)
1776 clon=(clat0*clatr*clonr-slat0*slatr)/clat
1777 clon=min(max(clon,-1.0),1.0)
1778
1779 if ((kgds(13)-kgds(8)) > 0) then
1780 hs = -1.0
1781 else
1782 hs = 1.0
1783 endif
1784
1785 rlon = mod(rlon0+hs*dpr*acos(clon)+3600,360.0)
1786
1787 print*,'got here last point ',kgds(12), kgds(13)
1788 print*,'got here last point rotated ', rlat, rlon
1789
1790 kgds(12)=nint(rlat*1000.) ! octs 11-13, Lat of
1791 kgds(13)=nint(rlon*1000.) ! octs 14-16, Lon of
1792
1793 kgds(9)=igdstmpl(17)
1794 kgds(10)=igdstmpl(18)
1795
1796 kgds(11) = 0 ! oct 28, scan mode
1797 if (btest(igdstmpl(19),7)) kgds(11) = 128
1798 if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64
1799 if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32
1800
1801 kgds(19)=0 ! oct 4, # vert coordinate parameters
1802 kgds(20)=255 ! oct 5, used for thinned grids, set to 255
1803
1804 res = ((float(kgds(9)) / 1.e6) + (float(kgds(10)) / 1.e6)) &
1805 * 0.5 * 111.0
1806
1807
1808 do i = 1, 25
1809 print*,'final kgds ',i,kgds(i)
1810 enddo
1811
1812
1813 elseif(igdtnum==30) then
1814
1815 kgds(1)=3 ! oct 6, lambert conformal
1816 kgds(2)=igdstmpl(8) ! octs 7-8, Ni
1817 ni = kgds(2)
1818 kgds(3)=igdstmpl(9) ! octs 9-10, Nj
1819 nj = kgds(3)
1820
1821 iscale = 1e6
1822 kgds(4) = nint(float(igdstmpl(10))/1000.0)
1823 kgds(5) = nint(float(igdstmpl(11))/1000.0)
1824
1825 kgds(6)=0 ! oct 17, resolution and component flags
1826 if (igdstmpl(1)==2 ) kgds(6)=64
1827 if ( btest(igdstmpl(12),4).OR.btest(igdstmpl(12),5) ) kgds(6)=kgds(6)+128
1828 if ( btest(igdstmpl(12),3) ) kgds(6)=kgds(6)+8
1829
1830 kgds(7) = nint(float(igdstmpl(14))/1000.0)
1831 kgds(8) = nint(float(igdstmpl(15))/1000.0)
1832 kgds(9) = nint(float(igdstmpl(16))/1000.0)
1833 kgds(10) = 0
1834
1835 kgds(11) = 0 ! oct 28, scan mode
1836 if (btest(igdstmpl(18),7)) kgds(11) = 128
1837 if (btest(igdstmpl(18),6)) kgds(11) = kgds(11) + 64
1838 if (btest(igdstmpl(18),5)) kgds(11) = kgds(11) + 32
1839
1840 kgds(12) = nint(float(igdstmpl(19))/1000.0)
1841 kgds(13) = nint(float(igdstmpl(20))/1000.0)
1842 kgds(14) = -90
1843 kgds(15) = 0
1844
1845 elseif(igdtnum==0) then ! lat/lon grid
1846
1847 iscale=igdstmpl(10)*igdstmpl(11)
1848 if (iscale == 0) iscale = 1e6
1849 kgds(1)=0 ! oct 6, data representation type.
1850 kgds(2)=igdstmpl(8) ! octs 7-8, Ni
1851 ni = kgds(2)
1852 kgds(3)=igdstmpl(9) ! octs 9-10, Nj
1853 nj = kgds(3)
1854 kgds(4)=nint(float(igdstmpl(12))/float(iscale)*1000.) ! octs 11-13, Lat of 1st grid point
1855 kgds(5)=nint(float(igdstmpl(13))/float(iscale)*1000.) ! octs 14-16, Lon of 1st grid point
1856
1857 kgds(6)=0 ! oct 17, resolution and component flags
1858 if (igdstmpl(1)==2 ) kgds(6)=64
1859 if ( btest(igdstmpl(14),4).OR.btest(igdstmpl(14),5) ) kgds(6)=kgds(6)+128
1860 if ( btest(igdstmpl(14),3) ) kgds(6)=kgds(6)+8
1861
1862 kgds(7)=nint(float(igdstmpl(15))/float(iscale)*1000.) ! octs 18-20, Lat of last grid point
1863 kgds(8)=nint(float(igdstmpl(16))/float(iscale)*1000.) ! octs 21-23, Lon of last grid point
1864 kgds(9)=nint(float(igdstmpl(17))/float(iscale)*1000.) ! octs 24-25, "i" resolution.
1865 kgds(10)=nint(float(igdstmpl(18))/float(iscale)*1000.) ! octs 26-27, "j" resolution.
1866
1867 kgds(11) = 0 ! oct 28, scan mode
1868 if (btest(igdstmpl(19),7)) kgds(11) = 128
1869 if (btest(igdstmpl(19),6)) kgds(11) = kgds(11) + 64
1870 if (btest(igdstmpl(19),5)) kgds(11) = kgds(11) + 32
1871
1872 kgds(12)=0 ! octs 29-32, reserved
1873 kgds(19)=0 ! oct 4, # vert coordinate parameters
1874 kgds(20)=255 ! oct 5, used for thinned grids, set to 255
1875
1876 else
1877
1878 call error_handler("UNRECOGNIZED INPUT GRID TYPE ", 1)
1879
1880 endif
1881
1882 end subroutine gdt_to_gds
1883
1884 end module model_grid
model_grid
Sets up the ESMF grid objects for the input data grid and target FV3 grid.
Definition model_grid.F90:9
model_grid::define_target_grid
subroutine, public define_target_grid(localpet, npets)
Setup the esmf grid object for the target grid.
Definition model_grid.F90:1027
model_grid::latitude_input_grid
type(esmf_field), public latitude_input_grid
latitude of grid center, input grid
Definition model_grid.F90:57
model_grid::ip1_input
integer, public ip1_input
i_input plus 1
Definition model_grid.F90:33
model_grid::cleanup_input_target_grid_data
subroutine, public cleanup_input_target_grid_data
Deallocate all esmf grid objects.
Definition model_grid.F90:1591
model_grid::latitude_target_grid
type(esmf_field), public latitude_target_grid
latitude of grid center, target grid
Definition model_grid.F90:78
model_grid::longitude_input_grid
type(esmf_field), public longitude_input_grid
longitude of grid center, input grid
Definition model_grid.F90:59
model_grid::i_target
integer, public i_target
i dimension of each global tile, or of a nest, target grid.
Definition model_grid.F90:37
model_grid::define_input_grid
subroutine, public define_input_grid(localpet, npets)
Driver routine to setup the esmf grid object for the input grid.
Definition model_grid.F90:118
model_grid::land_frac_target_grid
type(esmf_field), public land_frac_target_grid
land fraction, target grid
Definition model_grid.F90:76
model_grid::num_tiles_target_grid
integer, public num_tiles_target_grid
Number of tiles, target grid.
Definition model_grid.F90:49
model_grid::latitude_w_target_grid
type(esmf_field), public latitude_w_target_grid
latitude of 'west' edge of grid box, target grid
Definition model_grid.F90:83
model_grid::define_input_grid_grib2
subroutine define_input_grid_grib2(localpet, npets)
Define input grid object for grib2 input data.
Definition model_grid.F90:661
model_grid::landmask_target_grid
type(esmf_field), public landmask_target_grid
land mask target grid - '1' some or all land; '0' all non-land
Definition model_grid.F90:73
model_grid::input_grid
type(esmf_grid), public input_grid
input grid esmf grid object
Definition model_grid.F90:52
model_grid::longitude_w_target_grid
type(esmf_field), public longitude_w_target_grid
longitude of 'west' edge of grid box, target grid
Definition model_grid.F90:91
model_grid::define_input_grid_mosaic
subroutine define_input_grid_mosaic(localpet, npets)
Define input grid for tiled data using the 'mosaic', 'grid' and orography files.
Definition model_grid.F90:431
model_grid::jp1_input
integer, public jp1_input
j_input plus 1
Definition model_grid.F90:35
model_grid::terrain_target_grid
type(esmf_field), public terrain_target_grid
terrain height target grid
Definition model_grid.F90:97
model_grid::i_input
integer, public i_input
i-dimension of input grid (or of each global tile)
Definition model_grid.F90:27
model_grid::target_grid
type(esmf_grid), public target_grid
target grid esmf grid object.
Definition model_grid.F90:54
model_grid::longitude_s_input_grid
type(esmf_field), public longitude_s_input_grid
longitude of 'south' edge of grid box, input grid
Definition model_grid.F90:64
model_grid::get_model_mask_terrain
subroutine get_model_mask_terrain(orog_file, idim, jdim, mask, terrain, land_frac)
Read the model land mask and terrain for a single tile from the orography file.
Definition model_grid.F90:1516
model_grid::latitude_s_input_grid
type(esmf_field), public latitude_s_input_grid
latitude of 'south' edge of grid box, input grid
Definition model_grid.F90:61
model_grid::define_input_grid_gaussian
subroutine define_input_grid_gaussian(localpet, npets)
Define grid object for input data on global gaussian grids.
Definition model_grid.F90:150
model_grid::longitude_w_input_grid
type(esmf_field), public longitude_w_input_grid
longitude of 'west' edge of grid box, input grid
Definition model_grid.F90:70
model_grid::get_model_latlons
subroutine get_model_latlons(mosaic_file, orog_dir, num_tiles, tile, i_tile, j_tile, ip1_tile, jp1_tile, latitude, latitude_s, latitude_w, longitude, longitude_s, longitude_w)
Read model lat/lons for a single tile from the "grid" specificaton file.
Definition model_grid.F90:1367
model_grid::input_grid_type
character(len=50), public input_grid_type
map projection of input grid
Definition model_grid.F90:21
model_grid::longitude_target_grid
type(esmf_field), public longitude_target_grid
longitude of grid center, target grid
Definition model_grid.F90:86
model_grid::seamask_target_grid
type(esmf_field), public seamask_target_grid
sea mask target grid - '1' some or all non-land; '0' all land
Definition model_grid.F90:94
model_grid::gdt_to_gds
subroutine gdt_to_gds(igdtnum, igdstmpl, igdtlen, kgds, ni, nj, res)
Convert the GRIB2 grid description template to to the GRIB1 grid description section.
Definition model_grid.F90:1647
model_grid::tiles_target_grid
character(len=5), dimension(:), allocatable, public tiles_target_grid
Tile names of target grid.
Definition model_grid.F90:19
model_grid::num_tiles_input_grid
integer, public num_tiles_input_grid
Number of tiles, input grid.
Definition model_grid.F90:47
model_grid::latitude_w_input_grid
type(esmf_field), public latitude_w_input_grid
latitude of 'west' edge of grid box, input grid
Definition model_grid.F90:67
model_grid::j_target
integer, public j_target
j dimension of each global tile, or of a nest, target grid.
Definition model_grid.F90:40
model_grid::jp1_target
integer, public jp1_target
jp1_target plus 1
Definition model_grid.F90:45
model_grid::lsoil_target
integer, public lsoil_target
Number of soil layers, target grid.
Definition model_grid.F90:25
model_grid::j_input
integer, public j_input
j-dimension of input grid (or of each global tile)
Definition model_grid.F90:30
model_grid::longitude_s_target_grid
type(esmf_field), public longitude_s_target_grid
longitude of 'south' edge of grid box, target grid
Definition model_grid.F90:88
model_grid::ip1_target
integer, public ip1_target
ip1_target plus 1
Definition model_grid.F90:43
model_grid::latitude_s_target_grid
type(esmf_field), public latitude_s_target_grid
latitude of 'south' edge of grid box, target grid
Definition model_grid.F90:80
program_setup
This module contains code to read the setup namelist file, handle the varmap file for GRIB2 data,...
Definition program_setup.F90:9
program_setup::sfc_files_input_grid
character(len=500), dimension(6), public sfc_files_input_grid
File names containing input surface data.
Definition program_setup.F90:42
program_setup::orog_dir_input_grid
character(len=500), public orog_dir_input_grid
Directory containing the input grid orography files.
Definition program_setup.F90:37
program_setup::grib2_file_input_grid
character(len=500), public grib2_file_input_grid
REQUIRED.
Definition program_setup.F90:32
program_setup::orog_files_target_grid
character(len=500), dimension(6), public orog_files_target_grid
Target grid orography files.
Definition program_setup.F90:41
program_setup::orog_files_input_grid
character(len=500), dimension(6), public orog_files_input_grid
Input grid orography files.
Definition program_setup.F90:39
program_setup::orog_dir_target_grid
character(len=500), public orog_dir_target_grid
Directory containing the target grid orography files.
Definition program_setup.F90:40
program_setup::mosaic_file_target_grid
character(len=500), public mosaic_file_target_grid
Target grid mosaic file.
Definition program_setup.F90:30
program_setup::mosaic_file_input_grid
character(len=500), public mosaic_file_input_grid
Input grid mosaic file.
Definition program_setup.F90:29
program_setup::convert_sfc
logical, public convert_sfc
Convert sfc data when true.
Definition program_setup.F90:108
program_setup::nsoill_out
integer, public nsoill_out
Number of soil levels desired in the output data.
Definition program_setup.F90:103
program_setup::atm_files_input_grid
character(len=500), dimension(6), public atm_files_input_grid
File names of input atmospheric data.
Definition program_setup.F90:19
program_setup::convert_atm
logical, public convert_atm
Convert atmospheric data when true.
Definition program_setup.F90:106
program_setup::data_dir_input_grid
character(len=500), public data_dir_input_grid
Directory containing input atm or sfc files.
Definition program_setup.F90:27
program_setup::input_type
character(len=25), public input_type
Input data type:
Definition program_setup.F90:47