1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 ! MODULE GRIDINFO_MODULE
4 ! This module handles (i.e., acquires, stores, and makes available) all data
5 ! describing the model domains to be processed.
6 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
10 use misc_definitions_module
14 integer, parameter :: MAX_DOMAINS = 21
17 integer :: iproj_type, n_domains, io_form_output, dyn_opt
18 integer, dimension(MAX_DOMAINS) :: parent_grid_ratio, parent_id, ixdim, jydim
19 integer, dimension(MAX_DOMAINS) :: sr_x, sr_y
20 real :: known_lat, known_lon, pole_lat, pole_lon, stand_lon, truelat1, truelat2, &
21 known_x, known_y, dxkm, dykm, phi, lambda, ref_lat, ref_lon, ref_x, ref_y, &
23 real, dimension(MAX_DOMAINS) :: parent_ll_x, parent_ll_y, parent_ur_x, parent_ur_y
24 character (len=MAX_FILENAME_LEN) :: geog_data_path, opt_output_from_geogrid_path, opt_geogrid_tbl_path
26 character (len=128), dimension(MAX_DOMAINS) :: geog_data_res
27 character (len=1) :: gridtype
28 logical :: do_tiled_output
29 integer :: debug_level
33 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
34 ! Name: get_grid_params
36 ! Purpose: This subroutine retrieves all parameters regarding the model domains
37 ! to be processed by geogrid.exe. This includes map parameters, domain
38 ! size and location, and nest information.
39 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
40 subroutine get_grid_params()
45 integer :: i, j, max_dom, funit, io_form_geogrid, interval_seconds
46 real :: dx, dy, rparent_gridpts
47 integer, dimension(MAX_DOMAINS) :: i_parent_start, j_parent_start, &
48 s_we, e_we, s_sn, e_sn, &
49 start_year, start_month, start_day, start_hour, start_minute, start_second, &
50 end_year, end_month, end_day, end_hour, end_minute, end_second
51 character (len=128) :: map_proj
52 character (len=128), dimension(MAX_DOMAINS) :: start_date, end_date
53 character (len=3) :: wrf_core
54 logical :: is_used, nest_outside
56 namelist /share/ wrf_core, max_dom, start_date, end_date, &
57 start_year, end_year, start_month, end_month, &
58 start_day, end_day, start_hour, end_hour, &
59 start_minute, end_minute, start_second, end_second, &
61 io_form_geogrid, opt_output_from_geogrid_path, debug_level, &
63 namelist /geogrid/ parent_id, parent_grid_ratio, &
64 i_parent_start, j_parent_start, s_we, e_we, s_sn, e_sn, &
65 map_proj, ref_x, ref_y, ref_lat, ref_lon, &
66 pole_lat, pole_lon, truelat1, truelat2, stand_lon, &
67 dx, dy, geog_data_res, geog_data_path, opt_geogrid_tbl_path
69 ! Set defaults for namelist variables
74 geog_data_path = 'NOT_SPECIFIED'
88 geog_data_res(i) = 'default'
90 parent_grid_ratio(i) = INVALID
107 start_date(i) = '0000-00-00_00:00:00'
108 end_date(i) = '0000-00-00_00:00:00'
112 opt_output_from_geogrid_path = './'
113 opt_geogrid_tbl_path = 'geogrid/'
114 interval_seconds = INVALID
116 ! Read parameters from Fortran namelist
118 inquire(unit=funit, opened=is_used)
119 if (.not. is_used) exit
121 open(funit,file='namelist.wps',status='old',form='formatted',err=1000)
126 ! BUG: More properly handle debug_level in module_debug
127 if (debug_level.gt.100) then
128 call set_debug_level(DEBUG)
130 call set_debug_level(WARN)
133 call mprintf(.true.,LOGFILE,'Using the following namelist variables:')
134 call mprintf(.true.,LOGFILE,'&SHARE')
135 call mprintf(.true.,LOGFILE,' WRF_CORE = %s',s1=wrf_core)
136 call mprintf(.true.,LOGFILE,' MAX_DOM = %i',i1=max_dom)
137 call mprintf(.true.,LOGFILE,' START_YEAR = %i',i1=start_year(1))
139 call mprintf(.true.,LOGFILE,' = %i',i1=start_year(i))
141 call mprintf(.true.,LOGFILE,' START_MONTH = %i',i1=start_month(1))
143 call mprintf(.true.,LOGFILE,' = %i',i1=start_month(i))
145 call mprintf(.true.,LOGFILE,' START_DAY = %i',i1=start_day(1))
147 call mprintf(.true.,LOGFILE,' = %i',i1=start_day(i))
149 call mprintf(.true.,LOGFILE,' START_HOUR = %i',i1=start_hour(1))
151 call mprintf(.true.,LOGFILE,' = %i',i1=start_hour(i))
153 call mprintf(.true.,LOGFILE,' START_MINUTE = %i',i1=start_minute(1))
155 call mprintf(.true.,LOGFILE,' = %i',i1=start_minute(i))
157 call mprintf(.true.,LOGFILE,' START_SECOND = %i',i1=start_second(1))
159 call mprintf(.true.,LOGFILE,' = %i',i1=start_second(i))
161 call mprintf(.true.,LOGFILE,' END_YEAR = %i',i1=end_year(1))
163 call mprintf(.true.,LOGFILE,' = %i',i1=end_year(i))
165 call mprintf(.true.,LOGFILE,' END_MONTH = %i',i1=end_month(1))
167 call mprintf(.true.,LOGFILE,' = %i',i1=end_month(i))
169 call mprintf(.true.,LOGFILE,' END_DAY = %i',i1=end_day(1))
171 call mprintf(.true.,LOGFILE,' = %i',i1=end_day(i))
173 call mprintf(.true.,LOGFILE,' END_HOUR = %i',i1=end_hour(1))
175 call mprintf(.true.,LOGFILE,' = %i',i1=end_hour(i))
177 call mprintf(.true.,LOGFILE,' END_MINUTE = %i',i1=end_minute(1))
179 call mprintf(.true.,LOGFILE,' = %i',i1=end_minute(i))
181 call mprintf(.true.,LOGFILE,' END_SECOND = %i',i1=end_second(1))
183 call mprintf(.true.,LOGFILE,' = %i',i1=end_second(i))
185 call mprintf(.true.,LOGFILE,' START_DATE = %s',s1=start_date(1))
187 call mprintf(.true.,LOGFILE,' = %s',s1=start_date(i))
189 call mprintf(.true.,LOGFILE,' END_DATE = %s',s1=end_date(1))
191 call mprintf(.true.,LOGFILE,' = %s',s1=end_date(i))
193 call mprintf(.true.,LOGFILE,' INTERVAL_SECONDS = %i',i1=interval_seconds)
194 call mprintf(.true.,LOGFILE,' IO_FORM_GEOGRID = %i',i1=io_form_geogrid)
195 call mprintf(.true.,LOGFILE,' OPT_OUTPUT_FROM_GEOGRID_PATH = %s',s1=opt_output_from_geogrid_path)
196 call mprintf(.true.,LOGFILE,' SR_X = %i',i1=sr_x(1))
198 call mprintf(.true.,LOGFILE,' = %i',i1=sr_x(i))
200 call mprintf(.true.,LOGFILE,' SR_Y = %i',i1=sr_y(1))
202 call mprintf(.true.,LOGFILE,' = %i',i1=sr_y(i))
205 call mprintf(.true.,LOGFILE,' DEBUG_LEVEL = %i',i1=debug_level)
206 call mprintf(.true.,LOGFILE,'/')
208 call mprintf(.true.,LOGFILE,'&GEOGRID')
209 call mprintf(.true.,LOGFILE,' PARENT_ID = %i',i1=parent_id(1))
211 call mprintf(.true.,LOGFILE,' = %i',i1=parent_id(i))
213 call mprintf(.true.,LOGFILE,' PARENT_GRID_RATIO = %i',i1=parent_grid_ratio(1))
215 call mprintf(.true.,LOGFILE,' = %i',i1=parent_grid_ratio(i))
217 call mprintf(.true.,LOGFILE,' I_PARENT_START = %i',i1=i_parent_start(1))
219 call mprintf(.true.,LOGFILE,' = %i',i1=i_parent_start(i))
221 call mprintf(.true.,LOGFILE,' J_PARENT_START = %i',i1=j_parent_start(1))
223 call mprintf(.true.,LOGFILE,' = %i',i1=j_parent_start(i))
225 call mprintf(.true.,LOGFILE,' S_WE = %i',i1=s_we(1))
227 call mprintf(.true.,LOGFILE,' = %i',i1=s_we(i))
229 call mprintf(.true.,LOGFILE,' E_WE = %i',i1=e_we(1))
231 call mprintf(.true.,LOGFILE,' = %i',i1=e_we(i))
233 call mprintf(.true.,LOGFILE,' S_SN = %i',i1=s_sn(1))
235 call mprintf(.true.,LOGFILE,' = %i',i1=s_sn(i))
237 call mprintf(.true.,LOGFILE,' E_SN = %i',i1=e_sn(1))
239 call mprintf(.true.,LOGFILE,' = %i',i1=e_sn(i))
241 call mprintf(.true.,LOGFILE,' GEOG_DATA_RES = %s',s1=geog_data_res(1))
243 call mprintf(.true.,LOGFILE,' = %s',s1=geog_data_res(i))
245 call mprintf(.true.,LOGFILE,' DX = %f',f1=dx)
246 call mprintf(.true.,LOGFILE,' DY = %f',f1=dy)
247 call mprintf(.true.,LOGFILE,' MAP_PROJ = %s',s1=map_proj)
248 call mprintf(.true.,LOGFILE,' POLE_LAT = %f',f1=pole_lat)
249 call mprintf(.true.,LOGFILE,' POLE_LON = %f',f1=pole_lon)
250 call mprintf(.true.,LOGFILE,' REF_LAT = %f',f1=ref_lat)
251 call mprintf(.true.,LOGFILE,' REF_LON = %f',f1=ref_lon)
252 call mprintf(.true.,LOGFILE,' REF_X = %f',f1=ref_x)
253 call mprintf(.true.,LOGFILE,' REF_Y = %f',f1=ref_y)
254 call mprintf(.true.,LOGFILE,' TRUELAT1 = %f',f1=truelat1)
255 call mprintf(.true.,LOGFILE,' TRUELAT2 = %f',f1=truelat2)
256 call mprintf(.true.,LOGFILE,' STAND_LON = %f',f1=stand_lon)
257 call mprintf(.true.,LOGFILE,' GEOG_DATA_PATH = %s',s1=geog_data_path)
258 call mprintf(.true.,LOGFILE,' OPT_GEOGRID_TBL_PATH = %s',s1=opt_geogrid_tbl_path)
259 call mprintf(.true.,LOGFILE,'/')
269 ! Convert wrf_core to uppercase letters
271 if (ichar(wrf_core(i:i)) >= 97) wrf_core(i:i) = char(ichar(wrf_core(i:i))-32)
274 ! Before doing anything else, we must have a valid grid type
276 if (wrf_core == 'ARW') then
279 else if (wrf_core == 'NMM') then
284 ! Next, if this is NMM, we need to subtract 1 from the specified E_WE and E_SN;
285 ! for some reason, these two variables need to be set to 1 larger than they
286 ! really ought to be in the WRF namelist, so, to be consistent, we will do
287 ! the same in the WPS namelist
288 if (gridtype == 'E') then
290 e_we(i) = e_we(i) - 1
291 e_sn(i) = e_sn(i) - 1
295 call mprintf(gridtype /= 'C' .and. gridtype /= 'E', ERROR, &
296 'A valid wrf_core must be specified in the namelist. '// &
297 'Currently, only "ARW" and "NMM" are supported.')
299 call mprintf(max_dom > MAX_DOMAINS, ERROR, &
300 'In namelist, max_dom must be <= %i. To run with more'// &
301 ' than %i domains, increase the MAX_DOMAINS parameter.', &
302 i1=MAX_DOMAINS, i2=MAX_DOMAINS)
304 ! Every domain must have a valid parent id
306 call mprintf(parent_id(i) <= 0 .or. parent_id(i) >= i, ERROR, &
307 'In namelist, the parent_id of domain %i must be in '// &
308 'the range 1 to %i.', i1=i, i2=i-1)
311 ! Check for valid geog_data_path
313 do i=1,len(geog_data_path)
314 geog_data_path(j:j) = geog_data_path(i:i)
315 if (geog_data_path(i:i) /= ' ') j = j + 1
317 if (geog_data_path(1:1) == ' ') then
318 call mprintf(.true.,ERROR,'In namelist, geog_data_path must be specified.')
320 j = len_trim(geog_data_path)
321 if (j >= MAX_FILENAME_LEN) then
322 call mprintf(.true.,ERROR, &
323 'In namelist, geog_data_path must be strictly less '// &
324 'than 128 characters in length.')
326 if (geog_data_path(j:j) /= '/') then
327 geog_data_path(j+1:j+1) = '/'
331 ! Paths need to end with a /
332 j = len_trim(opt_geogrid_tbl_path)
333 if (opt_geogrid_tbl_path(j:j) /= '/') then
334 opt_geogrid_tbl_path(j+1:j+1) = '/'
337 j = len_trim(opt_output_from_geogrid_path)
338 if (opt_output_from_geogrid_path(j:j) /= '/') then
339 opt_output_from_geogrid_path(j+1:j+1) = '/'
342 ! Handle IOFORM+100 to do tiled IO
343 if (io_form_geogrid > 100) then
344 do_tiled_output = .true.
345 io_form_geogrid = io_form_geogrid - 100
347 do_tiled_output = .false.
350 ! Check for valid io_form_geogrid
353 io_form_geogrid /= BINARY .and. &
356 io_form_geogrid /= NETCDF .and. &
359 io_form_geogrid /= GRIB1 .and. &
362 call mprintf(.true.,WARN,'Valid io_form_geogrid values are:')
364 call mprintf(.true.,WARN,' %i (=BINARY)',i1=BINARY)
367 call mprintf(.true.,WARN,' %i (=NETCDF)',i1=NETCDF)
370 call mprintf(.true.,WARN,' %i (=GRIB1)',i1=GRIB1)
372 call mprintf(.true.,ERROR,'No valid value for io_form_geogrid was specified in the namelist.')
374 io_form_output = io_form_geogrid
376 ! Convert map_proj to uppercase letters
378 if (ichar(map_proj(i:i)) >= 97) map_proj(i:i) = char(ichar(map_proj(i:i))-32)
381 ! Assign parameters to module variables
382 if ((index(map_proj, 'LAMBERT') /= 0) .and. &
383 (len_trim(map_proj) == len('LAMBERT'))) then
386 else if ((index(map_proj, 'MERCATOR') /= 0) .and. &
387 (len_trim(map_proj) == len('MERCATOR'))) then
388 iproj_type = PROJ_MERC
390 else if ((index(map_proj, 'POLAR') /= 0) .and. &
391 (len_trim(map_proj) == len('POLAR'))) then
394 else if ((index(map_proj, 'ROTATED_LL') /= 0) .and. &
395 (len_trim(map_proj) == len('ROTATED_LL'))) then
396 iproj_type = PROJ_ROTLL
398 else if ((index(map_proj, 'LAT-LON') /= 0) .and. &
399 (len_trim(map_proj) == len('LAT-LON'))) then
400 iproj_type = PROJ_CASSINI
403 call mprintf(.true.,ERROR,&
404 'In namelist, invalid map_proj specified. Valid '// &
405 'projections are "lambert", "mercator", "polar", '// &
406 '"lat-lon", and "rotated_ll".')
409 ! For Cassini / lat-lon projections
410 if (iproj_type == PROJ_CASSINI) then
412 ! If no dx,dy specified, assume global grid
413 if (dx == NAN .and. dy == NAN) then
414 dlondeg = 360. / (e_we(1)-s_we(1)) ! Here, we really do not want e_we-s_we+1
415 dlatdeg = 180. / (e_sn(1)-s_sn(1)) ! Here, we really do not want e_we-s_we+1
418 known_lon = stand_lon + dlondeg/2.
419 known_lat = -90. + dlatdeg/2.
420 dxkm = EARTH_RADIUS_M * PI * 2.0 / (e_we(1)-s_we(1))
421 dykm = EARTH_RADIUS_M * PI / (e_sn(1)-s_sn(1))
423 ! If dx,dy specified, however, assume regional grid
427 dxkm = dlondeg * EARTH_RADIUS_M * PI * 2.0 / 360.0
428 dykm = dlatdeg * EARTH_RADIUS_M * PI * 2.0 / 360.0
429 if (known_lat == NAN .or. known_lon == NAN) then
430 call mprintf(.true.,ERROR,'For lat-lon projection, if dx/dy are specified, '// &
431 'a regional domain is assumed, and a ref_lat,ref_lon must also be specified')
436 ! Manually set truelat2 = truelat1 if truelat2 not specified for Lambert
437 if (iproj_type == PROJ_LC .and. truelat2 == NAN) then
438 call mprintf ((truelat1 == NAN), ERROR, "No TRUELAT1 specified for Lambert conformal projection.")
445 ! For C grid, let ixdim and jydim be the number of velocity points in
446 ! each direction; for E grid, we will put the row and column back
447 ! later; maybe this should be changed to be more clear, though.
449 ixdim(i) = e_we(i) - s_we(i) + 1
450 jydim(i) = e_sn(i) - s_sn(i) + 1
453 if (gridtype == 'E') then
454 phi = dykm*real(jydim(1)-1)/2.
455 lambda = dxkm*real(ixdim(1)-1)
458 ! If the user hasn't supplied a known_x and known_y, assume the center of domain 1
459 if (gridtype == 'E' .and. (known_x /= NAN .or. known_y /= NAN)) then
460 call mprintf(.true.,WARN, &
461 'Namelist variables ref_x and ref_y cannot be used for NMM grids.'// &
462 ' (ref_lat, ref_lon) will refer to the center of the coarse grid.')
463 else if (gridtype == 'C') then
464 if (known_x == NAN .and. known_y == NAN) then
465 known_x = ixdim(1) / 2.
466 known_y = jydim(1) / 2.
467 else if (known_x == NAN .or. known_y == NAN) then
468 call mprintf(.true.,ERROR, &
469 'In namelist.wps, neither or both of ref_x, ref_y must be specified.')
473 ! Checks specific to E grid
474 if (gridtype == 'E') then
476 ! E grid supports only the rotated lat/lon projection
477 if (iproj_type /= PROJ_ROTLL) then
478 call mprintf(.true., WARN, &
479 'For the NMM core, projection type must be rotated '// &
480 'lat/lon (map_proj=rotated_ll)')
481 call mprintf(.true.,WARN,'Projection will be set to rotated_ll')
482 iproj_type = PROJ_ROTLL
485 ! In the following check, add back the 1 that we had to subtract above
486 ! for the sake of being consistent with WRF namelist
487 call mprintf(mod(e_sn(1)+1,2) /= 0, ERROR, &
488 'For the NMM core, E_SN must be an even number for grid %i.', i1=1)
491 call mprintf((parent_grid_ratio(i) /= 3), WARN, &
492 'For the NMM core, the parent_grid_ratio must be 3 for '// &
493 'domain %i. A ratio of 3 will be assumed.', i1=i)
494 parent_grid_ratio(i) = 3
497 ! Check that nests have an acceptable number of grid points in each dimension
499 ! rparent_gridpts = real(ixdim(i)+2)/real(parent_grid_ratio(i))
500 ! if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then
501 ! call mprintf(.true.,ERROR,'For nest %i, e_we must be 3n-2 '// &
502 ! 'for some integer n > 1.', &
505 ! rparent_gridpts = real(jydim(i)+2)/real(parent_grid_ratio(i))
506 ! if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then
507 ! call mprintf(.true.,ERROR,'For nest %i, e_sn must be 3n-2 '// &
508 ! 'for some odd integer n > 1.', &
518 ! Checks specific to C grid
519 else if (gridtype == 'C') then
521 ! C grid does not support the rotated lat/lon projection
522 call mprintf((iproj_type == PROJ_ROTLL), ERROR, &
523 'Rotated lat/lon projection is not supported for the ARW core. '// &
524 'Valid projecitons are "lambert", "mercator", "polar", and "lat-lon".')
526 ! Check that nests have an acceptable number of grid points in each dimension
528 rparent_gridpts = real(ixdim(i)-1)/real(parent_grid_ratio(i))
529 if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then
530 call mprintf(.true.,ERROR,'For nest %i, (e_we-s_we+1) must be one greater '// &
531 'than an interger multiple of the parent_grid_ratio of %i.', &
532 i1=i, i2=parent_grid_ratio(i))
534 rparent_gridpts = real(jydim(i)-1)/real(parent_grid_ratio(i))
535 if (floor(rparent_gridpts) /= ceiling(rparent_gridpts)) then
536 call mprintf(.true.,ERROR,'For nest %i, (e_sn-s_sn+1) must be one greater '// &
537 'than an interger multiple of the parent_grid_ratio of %i.', &
538 i1=i, i2=parent_grid_ratio(i))
542 ! Check that a nest does not extend outside of its parent grid
543 nest_outside = .false.
545 if (i_parent_start(i) >= ixdim(parent_id(i))) then
546 call mprintf(.true.,WARN,'Nest %i cannot have i_parent_start outside of parent domain.',i1=i)
547 nest_outside = .true.
548 else if (i_parent_start(i) + (ixdim(i) - 1)/parent_grid_ratio(i) > ixdim(parent_id(i))) then
549 call mprintf(.true.,WARN,'Nest %i extends beyond its parent grid in the west-east direction.',i1=i)
550 call mprintf(.true.,WARN,' Maximum allowable e_we for current i/j_parent_start is %i.', &
551 i1=(ixdim(parent_id(i))-i_parent_start(i))*parent_grid_ratio(i)+1 )
552 nest_outside = .true.
554 if (j_parent_start(i) >= jydim(parent_id(i))) then
555 call mprintf(.true.,WARN,'Nest %i cannot have j_parent_start outside of parent domain.',i1=i)
556 nest_outside = .true.
557 else if (j_parent_start(i) + (jydim(i) - 1)/parent_grid_ratio(i) > jydim(parent_id(i))) then
558 call mprintf(.true.,WARN,'Nest %i extends beyond its parent grid in the south-north direction.',i1=i)
559 call mprintf(.true.,WARN,' Maximum allowable e_sn for current i/j_parent_start is %i.', &
560 i1=(jydim(parent_id(i))-j_parent_start(i))*parent_grid_ratio(i)+1 )
561 nest_outside = .true.
564 if (nest_outside) then
565 call mprintf(.true.,ERROR,'One or more nested domains extend beyond their parent domains.')
569 parent_ll_x(i) = real(i_parent_start(i))
570 parent_ll_y(i) = real(j_parent_start(i))
571 parent_ur_x(i) = real(i_parent_start(i))+real(ixdim(i))/real(parent_grid_ratio(i))-1.
572 parent_ur_y(i) = real(j_parent_start(i))+real(jydim(i))/real(parent_grid_ratio(i))-1.
579 1000 call mprintf(.true.,ERROR,'Error opening file namelist.wps')
581 end subroutine get_grid_params
583 end module gridinfo_module