1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 ! Module: source_data_module
5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
6 module source_data_module
11 use misc_definitions_module
14 integer, parameter :: RETURN_LANDMASK = 0, &
15 RETURN_DOMCAT_LM = 1, &
18 RETURN_FIELDNAME = 4, &
24 integer :: num_entries
25 integer :: next_field = 1
26 integer :: output_field_state = RETURN_LANDMASK
27 integer, pointer, dimension(:) :: source_proj, source_wordsize, source_endian, source_fieldtype, &
28 source_dest_fieldtype, source_priority, source_tile_x, source_tile_y, &
29 source_tile_z, source_tile_z_start, source_tile_z_end, source_tile_bdr, &
30 source_category_min, source_category_max, source_landmask_water, &
31 source_landmask_land, source_smooth_option, &
32 source_smooth_passes, source_output_stagger, source_row_order
33 integer :: source_iswater, source_isice, source_isurban, source_isoilwater
34 real, pointer, dimension(:) :: source_dx, source_dy, source_known_x, source_known_y, &
35 source_known_lat, source_known_lon, source_masked, source_truelat1, source_truelat2, &
36 source_stdlon, source_scale_factor, source_missing_value, source_fill_missing
37 character (len=128), pointer, dimension(:) :: source_fieldname, source_path, source_interp_string, &
38 source_dominant_category, source_dominant_only, source_dfdx, source_dfdy, &
39 source_z_dim_name, source_units, source_descr
40 character (len=128) :: source_mminlu
41 logical, pointer, dimension(:) :: is_proj, is_wordsize, is_endian, is_fieldtype, &
42 is_dest_fieldtype, is_priority, is_tile_x, is_tile_y, is_tile_z, &
43 is_tile_z_start, is_tile_z_end, is_tile_bdr, is_category_min, &
44 is_category_max, is_landmask_water, is_landmask_land, is_masked, &
45 is_dx, is_dy, is_known_x, is_known_y, &
46 is_known_lat, is_known_lon, is_truelat1, is_truelat2, is_stdlon, &
47 is_scale_factor, is_fieldname, is_path, is_dominant_category, &
48 is_dominant_only, is_dfdx, is_dfdy, is_z_dim_name, &
49 is_smooth_option, is_smooth_passes, is_signed, is_missing_value, &
50 is_fill_missing, is_halt_missing, is_output_stagger, is_row_order, &
51 is_units, is_descr, use_subgrid
53 type (list), pointer, dimension(:) :: source_res_path, source_interp_option
54 type (hashtable) :: bad_files ! Track which files produce errors when we try to open them
55 type (hashtable) :: duplicate_fnames ! Remember which output fields we have returned
61 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
65 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
66 subroutine get_datalist()
73 integer, parameter :: BUFSIZE = 256
76 integer :: nparams, idx, eos, ispace, i, funit
77 logical :: have_specification, is_used
78 character (len=128) :: res_string, path_string, interp_string
79 character (len=BUFSIZE) :: buffer
85 inquire(unit=funit, opened=is_used)
86 if (.not. is_used) exit
88 open(funit,file=trim(opt_geogrid_tbl_path)//'GEOGRID.TBL',form='formatted',status='old',err=1000)
91 ! We will first go through the file to determine how many field
92 ! specifications there are
94 10 read(funit,'(a)',end=20,err=1001) buffer
97 ! Is this line a comment?
98 if (buffer(1:1) == '#') then
100 ! Are we beginning a new field specification?
101 else if (index(buffer,'=====') /= 0) then
102 if (nparams > 0) num_entries = num_entries + 1
106 eos = index(buffer,'#')
107 if (eos /= 0) buffer(eos:BUFSIZE) = ' '
109 ! Does this line contain at least one parameter specification?
110 if (index(buffer,'=') /= 0) then
111 nparams = nparams + 1
118 ! In case the last entry ended without a ======== line
119 if (nparams > 0) num_entries = num_entries + 1
121 call mprintf(.true.,STDOUT,'Parsed %i entries in GEOGRID.TBL', i1=num_entries)
124 ! Now that we know how many fields the user has specified, allocate
125 ! the properly sized arrays
127 allocate(source_wordsize(num_entries))
128 allocate(source_endian(num_entries))
129 allocate(source_fieldtype(num_entries))
130 allocate(source_dest_fieldtype(num_entries))
131 allocate(source_proj(num_entries))
132 allocate(source_priority(num_entries))
133 allocate(source_dx(num_entries))
134 allocate(source_dy(num_entries))
135 allocate(source_known_x(num_entries))
136 allocate(source_known_y(num_entries))
137 allocate(source_known_lat(num_entries))
138 allocate(source_known_lon(num_entries))
139 allocate(source_truelat1(num_entries))
140 allocate(source_truelat2(num_entries))
141 allocate(source_stdlon(num_entries))
142 allocate(source_fieldname(num_entries))
143 allocate(source_path(num_entries))
144 allocate(source_interp_string(num_entries))
145 allocate(source_tile_x(num_entries))
146 allocate(source_tile_y(num_entries))
147 allocate(source_tile_z(num_entries))
148 allocate(source_tile_z_start(num_entries))
149 allocate(source_tile_z_end(num_entries))
150 allocate(source_category_min(num_entries))
151 allocate(source_category_max(num_entries))
152 allocate(source_tile_bdr(num_entries))
153 allocate(source_landmask_water(num_entries))
154 allocate(source_landmask_land(num_entries))
155 allocate(source_masked(num_entries))
156 allocate(source_output_stagger(num_entries))
157 allocate(source_row_order(num_entries))
158 allocate(source_dominant_category(num_entries))
159 allocate(source_dominant_only(num_entries))
160 allocate(source_dfdx(num_entries))
161 allocate(source_dfdy(num_entries))
162 allocate(source_scale_factor(num_entries))
163 allocate(source_z_dim_name(num_entries))
164 allocate(source_units(num_entries))
165 allocate(source_descr(num_entries))
166 allocate(source_smooth_option(num_entries))
167 allocate(source_smooth_passes(num_entries))
168 allocate(source_missing_value(num_entries))
169 allocate(source_fill_missing(num_entries))
170 allocate(source_res_path(num_entries))
171 allocate(source_interp_option(num_entries))
173 call list_init(source_res_path(i))
174 call list_init(source_interp_option(i))
177 allocate(is_wordsize(num_entries))
178 allocate(is_endian(num_entries))
179 allocate(is_fieldtype(num_entries))
180 allocate(is_dest_fieldtype(num_entries))
181 allocate(is_proj(num_entries))
182 allocate(is_priority(num_entries))
183 allocate(is_dx(num_entries))
184 allocate(is_dy(num_entries))
185 allocate(is_known_x(num_entries))
186 allocate(is_known_y(num_entries))
187 allocate(is_known_lat(num_entries))
188 allocate(is_known_lon(num_entries))
189 allocate(is_truelat1(num_entries))
190 allocate(is_truelat2(num_entries))
191 allocate(is_stdlon(num_entries))
192 allocate(is_fieldname(num_entries))
193 allocate(is_path(num_entries))
194 allocate(is_tile_x(num_entries))
195 allocate(is_tile_y(num_entries))
196 allocate(is_tile_z(num_entries))
197 allocate(is_tile_z_start(num_entries))
198 allocate(is_tile_z_end(num_entries))
199 allocate(is_category_min(num_entries))
200 allocate(is_category_max(num_entries))
201 allocate(is_tile_bdr(num_entries))
202 allocate(is_landmask_water(num_entries))
203 allocate(is_landmask_land(num_entries))
204 allocate(is_masked(num_entries))
205 allocate(is_halt_missing(num_entries))
206 allocate(is_output_stagger(num_entries))
207 allocate(is_row_order(num_entries))
208 allocate(is_dominant_category(num_entries))
209 allocate(is_dominant_only(num_entries))
210 allocate(is_dfdx(num_entries))
211 allocate(is_dfdy(num_entries))
212 allocate(is_scale_factor(num_entries))
213 allocate(is_z_dim_name(num_entries))
214 allocate(is_units(num_entries))
215 allocate(is_descr(num_entries))
216 allocate(is_smooth_option(num_entries))
217 allocate(is_smooth_passes(num_entries))
218 allocate(is_signed(num_entries))
219 allocate(is_missing_value(num_entries))
220 allocate(is_fill_missing(num_entries))
221 allocate(use_subgrid(num_entries))
223 write(source_mminlu,'(a4)') 'USGS'
227 source_isoilwater = 14
230 ! Actually read and save the specifications
235 read(funit,'(a)',end=40,err=1001) buffer
238 ! Is this line a comment?
239 if (buffer(1:1) == '#') then
242 ! Are we beginning a new field specification?
243 else if (index(buffer,'=====') /= 0) then !{
244 if (nparams > 0) i = i + 1
246 if (i <= num_entries) then
248 is_landmask_water(i) = .false.
249 is_landmask_land(i) = .false.
250 is_masked(i) = .false.
251 is_halt_missing(i) = .false.
252 is_output_stagger(i) = .false.
253 is_dominant_category(i) = .false.
254 is_dominant_only(i) = .false.
257 is_dest_fieldtype(i) = .false.
258 is_priority(i) = .false.
259 is_z_dim_name(i) = .false.
260 is_smooth_option(i) = .false.
261 is_smooth_passes(i) = .false.
262 is_fill_missing(i) = .false.
263 use_subgrid(i) = .false.
267 ! Check whether the current line is a comment
268 if (buffer(1:1) /= '#') then
269 have_specification = .true.
271 have_specification = .false.
274 ! If only part of the line is a comment, just turn the comment into spaces
275 eos = index(buffer,'#')
276 if (eos /= 0) buffer(eos:BUFSIZE) = ' '
278 do while (have_specification) !{
280 ! If this line has no semicolon, it may contain a single specification,
281 ! so we set have_specification = .false. to prevent the line from being
282 ! processed again and "pretend" that the last character was a semicolon
283 eos = index(buffer,';')
285 have_specification = .false.
289 idx = index(buffer(1:eos-1),'=')
291 if (idx /= 0) then !{
292 nparams = nparams + 1
294 if (index('name',trim(buffer(1:idx-1))) /= 0) then
296 do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
299 is_fieldname(i) = .true.
300 source_fieldname(i) = ' '
301 source_fieldname(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
302 ! print*,i,' ',buffer(idx+1:ispace-1)
304 else if (index('priority',trim(buffer(1:idx-1))) /= 0) then
305 is_priority(i) = .true.
306 read(buffer(idx+1:eos-1),'(i10)') source_priority(i)
308 else if (index('dest_type',trim(buffer(1:idx-1))) /= 0) then
309 if (index('continuous',trim(buffer(idx+1:eos-1))) /= 0) then
310 is_dest_fieldtype(i) = .true.
311 source_dest_fieldtype(i) = CONTINUOUS
312 else if (index('categorical',trim(buffer(idx+1:eos-1))) /= 0) then
313 is_dest_fieldtype(i) = .true.
314 source_dest_fieldtype(i) = CATEGORICAL
317 else if (index('interp_option',trim(buffer(1:idx-1))) /= 0) then
319 do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
323 interp_string(1:ispace-idx-1) = buffer(idx+1:ispace-1)
324 ispace = index(interp_string,':')
325 if (ispace /= 0) then
326 write(res_string,'(a)') interp_string(1:ispace-1)
328 res_string = 'default'
330 write(interp_string,'(a)') trim(interp_string(ispace+1:128))
331 if (list_search(source_interp_option(i), ckey=res_string, cvalue=interp_string)) then
332 call mprintf(.true., WARN, &
333 'In GEOGRID.TBL entry %i, multiple interpolation methods are '// &
334 'given for resolution %s. %s will be used.', &
335 i1=i, s1=trim(res_string), s2=trim(interp_string))
337 call list_insert(source_interp_option(i), ckey=res_string, cvalue=interp_string)
340 else if (index('smooth_option',trim(buffer(1:idx-1))) /= 0) then
341 if ((index('1-2-1',trim(buffer(idx+1:eos-1))) /= 0) .and. &
342 (len_trim(buffer(idx+1:eos-1)) == len('1-2-1'))) then
343 is_smooth_option(i) = .true.
344 source_smooth_option(i) = ONETWOONE
345 else if ((index('smth-desmth',trim(buffer(idx+1:eos-1))) /= 0) .and. &
346 (len_trim(buffer(idx+1:eos-1)) == len('smth-desmth'))) then
347 is_smooth_option(i) = .true.
348 source_smooth_option(i) = SMTHDESMTH
349 else if ((index('smth-desmth_special',trim(buffer(idx+1:eos-1))) /= 0) .and. &
350 (len_trim(buffer(idx+1:eos-1)) == len('smth-desmth_special'))) then
351 is_smooth_option(i) = .true.
352 source_smooth_option(i) = SMTHDESMTH_SPECIAL
355 else if (index('smooth_passes',trim(buffer(1:idx-1))) /= 0) then
356 is_smooth_passes(i) = .true.
357 read(buffer(idx+1:eos-1),'(i10)') source_smooth_passes(i)
359 else if (index('rel_path',trim(buffer(1:idx-1))) /= 0) then
361 do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
365 path_string(1:ispace-idx-1) = buffer(idx+1:ispace-1)
366 if (path_string(ispace-idx-1:ispace-idx-1) /= '/') &
367 path_string(ispace-idx:ispace-idx) = '/'
368 if (path_string(1:1) == '/') then
369 path_string(1:127) = path_string(2:128)
370 path_string(128:128) = ' '
372 ispace = index(path_string,':')
373 if (ispace /= 0) then
374 write(res_string,'(a)') path_string(1:ispace-1)
376 res_string = 'default'
378 write(path_string,'(a)') trim(geog_data_path)//trim(path_string(ispace+1:128))
379 if (list_search(source_res_path(i), ckey=res_string, cvalue=path_string)) then
380 call mprintf(.true., WARN, &
381 'In GEOGRID.TBL entry %i, multiple paths are given for '// &
382 'resolution %s. %s will be used.', &
383 i1=i, s1=trim(res_string), s2=trim(path_string))
385 call list_insert(source_res_path(i), ckey=res_string, cvalue=path_string)
388 else if (index('abs_path',trim(buffer(1:idx-1))) /= 0) then
390 do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
394 path_string(1:ispace-idx-1) = buffer(idx+1:ispace-1)
395 if (path_string(ispace-idx-1:ispace-idx-1) /= '/') &
396 path_string(ispace-idx:ispace-idx) = '/'
397 ispace = index(path_string,':')
398 if (ispace /= 0) then
399 write(res_string,'(a)') path_string(1:ispace-1)
401 res_string = 'default'
403 write(path_string,'(a)') trim(path_string(ispace+1:128))
404 if (list_search(source_res_path(i), ckey=res_string, cvalue=path_string)) then
405 call mprintf(.true., WARN, &
406 'In GEOGRID.TBL entry %i, multiple paths are given for '// &
407 'resolution %s. %s will be used.', &
408 i1=i, s1=trim(res_string), s2=trim(path_string))
410 call list_insert(source_res_path(i), ckey=res_string, cvalue=path_string)
413 else if (index('output_stagger',trim(buffer(1:idx-1))) /= 0) then
414 if (index('M',trim(buffer(idx+1:eos-1))) /= 0) then
415 is_output_stagger(i) = .true.
416 source_output_stagger(i) = M
417 else if (index('U',trim(buffer(idx+1:eos-1))) /= 0) then
418 is_output_stagger(i) = .true.
419 source_output_stagger(i) = U
420 else if (index('V',trim(buffer(idx+1:eos-1))) /= 0) then
421 is_output_stagger(i) = .true.
422 source_output_stagger(i) = V
423 else if (index('HH',trim(buffer(idx+1:eos-1))) /= 0) then
424 is_output_stagger(i) = .true.
425 source_output_stagger(i) = HH
426 else if (index('VV',trim(buffer(idx+1:eos-1))) /= 0) then
427 is_output_stagger(i) = .true.
428 source_output_stagger(i) = VV
431 else if ((index('landmask_water',trim(buffer(1:idx-1))) /= 0) .and. &
432 (len_trim(buffer(1:idx-1)) == 14)) then
433 is_landmask_water(i) = .true.
434 read(buffer(idx+1:eos-1),'(i10)') source_landmask_water(i)
436 else if ((index('landmask_land',trim(buffer(1:idx-1))) /= 0) .and. &
437 (len_trim(buffer(1:idx-1)) == 13)) then
438 is_landmask_land(i) = .true.
439 read(buffer(idx+1:eos-1),'(i10)') source_landmask_land(i)
441 else if ((index('masked',trim(buffer(1:idx-1))) /= 0) .and. &
442 (len_trim(buffer(1:idx-1)) == 6)) then
443 if (index('water',trim(buffer(idx+1:eos-1))) /= 0) then
444 is_masked(i) = .true.
445 source_masked(i) = 0.
446 else if (index('land',trim(buffer(idx+1:eos-1))) /= 0) then
447 is_masked(i) = .true.
448 source_masked(i) = 1.
451 else if (index('fill_missing',trim(buffer(1:idx-1))) /= 0) then
452 is_fill_missing(i) = .true.
453 read(buffer(idx+1:eos-1),*) source_fill_missing(i)
455 else if (index('halt_on_missing',trim(buffer(1:idx-1))) /= 0) then
456 if (index('yes',trim(buffer(idx+1:eos-1))) /= 0) then
457 is_halt_missing(i) = .true.
458 else if (index('no',trim(buffer(idx+1:eos-1))) /= 0) then
459 is_halt_missing(i) = .false.
462 else if (index('dominant_category',trim(buffer(1:idx-1))) /= 0) then
464 do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
467 is_dominant_category(i) = .true.
468 source_dominant_category(i) = ' '
469 source_dominant_category(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
471 else if (index('dominant_only',trim(buffer(1:idx-1))) /= 0) then
473 do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
476 is_dominant_only(i) = .true.
477 source_dominant_only(i) = ' '
478 source_dominant_only(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
480 else if (index('df_dx',trim(buffer(1:idx-1))) /= 0) then
482 do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
487 source_dfdx(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
489 else if (index('df_dy',trim(buffer(1:idx-1))) /= 0) then
491 do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
496 source_dfdy(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
498 else if (index('z_dim_name',trim(buffer(1:idx-1))) /= 0) then
500 do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
503 is_z_dim_name(i) = .true.
504 source_z_dim_name(i) = ' '
505 source_z_dim_name(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
507 else if(index('subgrid',trim(buffer(1:idx-1))) /= 0) then
509 do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
512 use_subgrid(i)=.true.
515 call mprintf(.true., WARN, 'In GEOGRID.TBL, unrecognized option %s in '// &
516 'entry %i.',i1=idx, s1=buffer(i:idx-1))
520 end if !} index(buffer(1:eos-1),'=') /= 0
522 buffer = buffer(eos+1:BUFSIZE)
523 end do ! while eos /= 0 }
525 end if !} index(buffer, '=====') /= 0
530 ! Check the user specifications for gross errors
531 if ( .not. check_data_specification() ) then
532 call datalist_destroy()
533 call mprintf(.true.,ERROR,'Errors were found in either index files or GEOGRID.TBL.')
536 call hash_init(bad_files)
540 1000 call mprintf(.true.,ERROR,'Could not open GEOGRID.TBL')
542 1001 call mprintf(.true.,ERROR,'Encountered error while reading GEOGRID.TBL')
544 end subroutine get_datalist
547 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
548 ! Name: get_source_params
551 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
552 subroutine get_source_params(res_string)
558 integer, parameter :: BUFSIZE = 256
561 character (len=128), intent(in) :: res_string
564 integer :: idx, i, ispace, eos, iquoted, funit
565 character (len=128) :: temp_str
566 character (len=256) :: test_fname
567 character (len=BUFSIZE) :: buffer
568 logical :: have_specification, is_used
571 is_wordsize(idx) = .false.
572 is_endian(idx) = .false.
573 is_row_order(idx) = .false.
574 is_fieldtype(idx) = .false.
575 is_proj(idx) = .false.
578 is_known_x(idx) = .false.
579 is_known_y(idx) = .false.
580 is_known_lat(idx) = .false.
581 is_known_lon(idx) = .false.
582 is_truelat1(idx) = .false.
583 is_truelat2(idx) = .false.
584 is_stdlon(idx) = .false.
585 is_tile_x(idx) = .false.
586 is_tile_y(idx) = .false.
587 is_tile_z(idx) = .false.
588 is_tile_z_start(idx) = .false.
589 is_tile_z_end(idx) = .false.
590 is_category_min(idx) = .false.
591 is_category_max(idx) = .false.
592 is_tile_bdr(idx) = .false.
593 is_fieldname(idx) = .false.
594 is_scale_factor(idx) = .false.
595 is_units(idx) = .false.
596 is_descr(idx) = .false.
597 is_signed(idx) = .false.
598 is_missing_value(idx) = .false.
600 if (.not. list_search(source_interp_option(idx), ckey=res_string, cvalue=source_interp_string(idx))) then
601 call mprintf(.true., INFORM, 'For %s, couldn''t find interpolator sequence for '// &
602 'resolution %s. Trying default.', &
603 s1=trim(source_fieldname(idx)), s2=trim(res_string))
605 if (list_search(source_interp_option(idx), ckey=temp_str, cvalue=source_interp_string(idx))) then
606 call mprintf(.true., INFORM, 'Using default interpolator sequence for %s.', &
607 s1=trim(source_fieldname(idx)))
609 call mprintf(.true., ERROR, 'Could not find default or %s interpolator sequence.'// &
610 ' The source specified in namelist.wps is not listed in GEOGRID.TBL.', &
615 if (.not. list_search(source_res_path(idx), ckey=res_string, cvalue=source_path(idx))) then
616 call mprintf(.true., INFORM, 'For %s, couldn''t find %s source. Trying '// &
617 'default.', s1=trim(source_fieldname(idx)), s2=trim(res_string))
619 if (list_search(source_res_path(idx), ckey=temp_str, cvalue=source_path(idx))) then
620 call mprintf(.true., INFORM, 'Using default data source for %s.', &
621 s1=trim(source_fieldname(idx)))
623 call mprintf(.true., ERROR, 'Could not find default or %s data source. '// &
624 'The source specified in namelist.wps is not listed in GEOGRID.TBL.', &
628 write(test_fname, '(a)') trim(source_path(idx))//'index'
631 inquire(unit=funit, opened=is_used)
632 if (.not. is_used) exit
634 open(funit,file=trim(test_fname),form='formatted',status='old',err=1000)
637 read(funit,'(a)',end=40, err=1001) buffer
640 ! Is this line a comment?
641 if (buffer(1:1) == '#') then
645 have_specification = .true.
647 ! If only part of the line is a comment, just turn the comment into spaces
648 eos = index(buffer,'#')
649 if (eos /= 0) buffer(eos:BUFSIZE) = ' '
651 do while (have_specification) !{
653 ! If this line has no semicolon, it may contain a single specification,
654 ! so we set have_specification = .false. to prevent the line from being
655 ! processed again and pretend that the last character was a semicolon
656 eos = index(buffer,';')
658 have_specification = .false.
662 i = index(buffer(1:eos-1),'=')
666 if (index('projection',trim(buffer(1:i-1))) /= 0) then
667 if (index('lambert',trim(buffer(i+1:eos-1))) /= 0) then
668 is_proj(idx) = .true.
669 source_proj(idx) = PROJ_LC
670 else if (index('polar_wgs84',trim(buffer(i+1:eos-1))) /= 0 .and. &
671 len_trim('polar_wgs84') == len_trim(buffer(i+1:eos-1))) then
672 is_proj(idx) = .true.
673 source_proj(idx) = PROJ_PS_WGS84
674 else if (index('albers_nad83',trim(buffer(i+1:eos-1))) /= 0 .and. &
675 len_trim('albers_nad83') == len_trim(buffer(i+1:eos-1))) then
676 is_proj(idx) = .true.
677 source_proj(idx) = PROJ_ALBERS_NAD83
678 else if (index('polar',trim(buffer(i+1:eos-1))) /= 0 .and. &
679 len_trim('polar') == len_trim(buffer(i+1:eos-1))) then
680 is_proj(idx) = .true.
681 source_proj(idx) = PROJ_PS
682 else if (index('mercator',trim(buffer(i+1:eos-1))) /= 0) then
683 is_proj(idx) = .true.
684 source_proj(idx) = PROJ_MERC
685 else if (index('regular_ll',trim(buffer(i+1:eos-1))) /= 0) then
686 is_proj(idx) = .true.
687 source_proj(idx) = PROJ_LATLON
690 else if (index('type',trim(buffer(1:i-1))) /= 0) then
691 if (index('continuous',trim(buffer(i+1:eos-1))) /= 0) then
692 is_fieldtype(idx) = .true.
693 source_fieldtype(idx) = CONTINUOUS
694 else if (index('categorical',trim(buffer(i+1:eos-1))) /= 0) then
695 is_fieldtype(idx) = .true.
696 source_fieldtype(idx) = CATEGORICAL
699 else if (index('signed',trim(buffer(1:i-1))) /= 0) then
700 if (index('yes',trim(buffer(i+1:eos-1))) /= 0) then
701 is_signed(idx) = .true.
702 else if (index('no',trim(buffer(i+1:eos-1))) /= 0) then
703 is_signed(idx) = .false.
706 else if (index('units',trim(buffer(1:i-1))) /= 0) then
709 do while (((ispace < eos) .and. (buffer(ispace:ispace) /= ' ')) .or. (iquoted == 1))
710 if (buffer(ispace:ispace) == '"' .or. buffer(ispace:ispace) == '''') iquoted = mod(iquoted+1,2)
713 is_units(idx) = .true.
714 source_units(idx) = ' '
715 if (buffer(i+1:i+1) == '"' .or. buffer(i+1:i+1) == '''') i = i + 1
716 if (buffer(ispace-1:ispace-1) == '"' .or. buffer(ispace-1:ispace-1) == '''') ispace = ispace - 1
717 source_units(idx)(1:ispace-i) = buffer(i+1:ispace-1)
719 else if (index('description',trim(buffer(1:i-1))) /= 0) then
722 do while (((ispace < eos) .and. (buffer(ispace:ispace) /= ' ')) .or. (iquoted == 1))
723 if (buffer(ispace:ispace) == '"' .or. buffer(ispace:ispace) == '''') iquoted = mod(iquoted+1,2)
726 is_descr(idx) = .true.
727 source_descr(idx) = ' '
728 if (buffer(i+1:i+1) == '"' .or. buffer(i+1:i+1) == '''') i = i + 1
729 if (buffer(ispace-1:ispace-1) == '"' .or. buffer(ispace-1:ispace-1) == '''') ispace = ispace - 1
730 source_descr(idx)(1:ispace-i) = buffer(i+1:ispace-1)
732 else if (index('mminlu',trim(buffer(1:i-1))) /= 0) then
735 do while (((ispace < eos) .and. (buffer(ispace:ispace) /= ' ')) .or. (iquoted == 1))
736 if (buffer(ispace:ispace) == '"' .or. buffer(ispace:ispace) == '''') iquoted = mod(iquoted+1,2)
739 if (buffer(i+1:i+1) == '"' .or. buffer(i+1:i+1) == '''') i = i + 1
740 if (buffer(ispace-1:ispace-1) == '"' .or. buffer(ispace-1:ispace-1) == '''') ispace = ispace - 1
741 source_mminlu(1:ispace-i) = buffer(i+1:ispace-1)
743 else if (index('iswater',trim(buffer(1:i-1))) /= 0) then
744 read(buffer(i+1:eos-1),*) source_iswater
746 else if (index('isice',trim(buffer(1:i-1))) /= 0) then
747 read(buffer(i+1:eos-1),*) source_isice
749 else if (index('isurban',trim(buffer(1:i-1))) /= 0) then
750 read(buffer(i+1:eos-1),*) source_isurban
752 else if (index('isoilwater',trim(buffer(1:i-1))) /= 0) then
753 read(buffer(i+1:eos-1),*) source_isoilwater
755 else if (index('dx',trim(buffer(1:i-1))) /= 0) then
757 read(buffer(i+1:eos-1),*) source_dx(idx)
759 else if (index('dy',trim(buffer(1:i-1))) /= 0) then
761 read(buffer(i+1:eos-1),*) source_dy(idx)
763 else if (index('known_x',trim(buffer(1:i-1))) /= 0) then
764 is_known_x(idx) = .true.
765 read(buffer(i+1:eos-1),*) source_known_x(idx)
767 else if (index('known_y',trim(buffer(1:i-1))) /= 0) then
768 is_known_y(idx) = .true.
769 read(buffer(i+1:eos-1),*) source_known_y(idx)
771 else if (index('known_lat',trim(buffer(1:i-1))) /= 0) then
772 is_known_lat(idx) = .true.
773 read(buffer(i+1:eos-1),*) source_known_lat(idx)
775 else if (index('known_lon',trim(buffer(1:i-1))) /= 0) then
776 is_known_lon(idx) = .true.
777 read(buffer(i+1:eos-1),*) source_known_lon(idx)
779 else if (index('stdlon',trim(buffer(1:i-1))) /= 0) then
780 is_stdlon(idx) = .true.
781 read(buffer(i+1:eos-1),*) source_stdlon(idx)
783 else if (index('truelat1',trim(buffer(1:i-1))) /= 0) then
784 is_truelat1(idx) = .true.
785 read(buffer(i+1:eos-1),*) source_truelat1(idx)
787 else if (index('truelat2',trim(buffer(1:i-1))) /= 0) then
788 is_truelat2(idx) = .true.
789 read(buffer(i+1:eos-1),*) source_truelat2(idx)
791 else if (index('wordsize',trim(buffer(1:i-1))) /= 0) then
792 is_wordsize(idx) = .true.
793 read(buffer(i+1:eos-1),'(i10)') source_wordsize(idx)
795 else if (index('endian',trim(buffer(1:i-1))) /= 0) then
796 if (index('big',trim(buffer(i+1:eos-1))) /= 0) then
797 is_endian(idx) = .true.
798 source_endian(idx) = BIG_ENDIAN
799 else if (index('little',trim(buffer(i+1:eos-1))) /= 0) then
800 is_endian(idx) = .true.
801 source_endian(idx) = LITTLE_ENDIAN
803 call mprintf(.true.,WARN,'Invalid value for keyword ''endian'' '// &
804 'specified in index file. BIG_ENDIAN will be used.')
807 else if (index('row_order',trim(buffer(1:i-1))) /= 0) then
808 if (index('bottom_top',trim(buffer(i+1:eos-1))) /= 0) then
809 is_row_order(idx) = .true.
810 source_row_order(idx) = BOTTOM_TOP
811 else if (index('top_bottom',trim(buffer(i+1:eos-1))) /= 0) then
812 is_row_order(idx) = .true.
813 source_row_order(idx) = TOP_BOTTOM
816 else if (index('tile_x',trim(buffer(1:i-1))) /= 0) then
817 is_tile_x(idx) = .true.
818 read(buffer(i+1:eos-1),'(i10)') source_tile_x(idx)
820 else if (index('tile_y',trim(buffer(1:i-1))) /= 0) then
821 is_tile_y(idx) = .true.
822 read(buffer(i+1:eos-1),'(i10)') source_tile_y(idx)
824 else if (index('tile_z',trim(buffer(1:i-1))) /= 0) then
825 is_tile_z(idx) = .true.
826 read(buffer(i+1:eos-1),'(i10)') source_tile_z(idx)
828 else if (index('tile_z_start',trim(buffer(1:i-1))) /= 0) then
829 is_tile_z_start(idx) = .true.
830 read(buffer(i+1:eos-1),'(i10)') source_tile_z_start(idx)
832 else if (index('tile_z_end',trim(buffer(1:i-1))) /= 0) then
833 is_tile_z_end(idx) = .true.
834 read(buffer(i+1:eos-1),'(i10)') source_tile_z_end(idx)
836 else if (index('category_min',trim(buffer(1:i-1))) /= 0) then
837 is_category_min(idx) = .true.
838 read(buffer(i+1:eos-1),'(i10)') source_category_min(idx)
840 else if (index('category_max',trim(buffer(1:i-1))) /= 0) then
841 is_category_max(idx) = .true.
842 read(buffer(i+1:eos-1),'(i10)') source_category_max(idx)
844 else if (index('tile_bdr',trim(buffer(1:i-1))) /= 0) then
845 is_tile_bdr(idx) = .true.
846 read(buffer(i+1:eos-1),'(i10)') source_tile_bdr(idx)
848 else if (index('missing_value',trim(buffer(1:i-1))) /= 0) then
849 is_missing_value(idx) = .true.
850 read(buffer(i+1:eos-1),*) source_missing_value(idx)
852 else if (index('scale_factor',trim(buffer(1:i-1))) /= 0) then
853 is_scale_factor(idx) = .true.
854 read(buffer(i+1:eos-1),*) source_scale_factor(idx)
857 call mprintf(.true., WARN, 'In %s, unrecognized option %s in entry %i.', &
858 s1=trim(test_fname), s2=buffer(i:i-1), i1=i)
862 end if !} index(buffer(1:eos-1),'=') /= 0
864 buffer = buffer(eos+1:BUFSIZE)
865 end do ! while eos /= 0 }
867 end if !} index(buffer, '=====') /= 0
877 1000 call mprintf(.true.,ERROR,'Could not open %s', s1=trim(test_fname))
878 1001 call mprintf(.true.,ERROR,'Encountered error while reading %s', s1=trim(test_fname))
880 end subroutine get_source_params
883 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
884 ! Name: datalist_destroy()
886 ! Purpose: This routine deallocates any memory that was allocated by the
887 ! get_datalist() subroutine.
888 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
889 subroutine datalist_destroy()
896 if (associated(source_wordsize)) then
897 deallocate(source_wordsize)
898 deallocate(source_endian)
899 deallocate(source_fieldtype)
900 deallocate(source_dest_fieldtype)
901 deallocate(source_proj)
902 deallocate(source_priority)
903 deallocate(source_dx)
904 deallocate(source_dy)
905 deallocate(source_known_x)
906 deallocate(source_known_y)
907 deallocate(source_known_lat)
908 deallocate(source_known_lon)
909 deallocate(source_truelat1)
910 deallocate(source_truelat2)
911 deallocate(source_stdlon)
912 deallocate(source_fieldname)
913 deallocate(source_path)
914 deallocate(source_interp_string)
915 deallocate(source_tile_x)
916 deallocate(source_tile_y)
917 deallocate(source_tile_z)
918 deallocate(source_tile_z_start)
919 deallocate(source_tile_z_end)
920 deallocate(source_tile_bdr)
921 deallocate(source_category_min)
922 deallocate(source_category_max)
923 deallocate(source_landmask_water)
924 deallocate(source_landmask_land)
925 deallocate(source_masked)
926 deallocate(source_output_stagger)
927 deallocate(source_row_order)
928 deallocate(source_dominant_category)
929 deallocate(source_dominant_only)
930 deallocate(source_dfdx)
931 deallocate(source_dfdy)
932 deallocate(source_scale_factor)
933 deallocate(source_z_dim_name)
934 deallocate(source_smooth_option)
935 deallocate(source_smooth_passes)
936 deallocate(source_units)
937 deallocate(source_descr)
938 deallocate(source_missing_value)
939 deallocate(source_fill_missing)
941 call list_destroy(source_res_path(i))
942 call list_destroy(source_interp_option(i))
944 deallocate(source_res_path)
945 deallocate(source_interp_option)
947 deallocate(is_wordsize)
948 deallocate(is_endian)
949 deallocate(is_fieldtype)
950 deallocate(is_dest_fieldtype)
952 deallocate(is_priority)
955 deallocate(is_known_x)
956 deallocate(is_known_y)
957 deallocate(is_known_lat)
958 deallocate(is_known_lon)
959 deallocate(is_truelat1)
960 deallocate(is_truelat2)
961 deallocate(is_stdlon)
962 deallocate(is_fieldname)
964 deallocate(is_tile_x)
965 deallocate(is_tile_y)
966 deallocate(is_tile_z)
967 deallocate(is_tile_z_start)
968 deallocate(is_tile_z_end)
969 deallocate(is_tile_bdr)
970 deallocate(is_category_min)
971 deallocate(is_category_max)
972 deallocate(is_landmask_water)
973 deallocate(is_landmask_land)
974 deallocate(is_masked)
975 deallocate(is_halt_missing)
976 deallocate(is_output_stagger)
977 deallocate(is_row_order)
978 deallocate(is_dominant_category)
979 deallocate(is_dominant_only)
982 deallocate(is_scale_factor)
983 deallocate(is_z_dim_name)
984 deallocate(is_smooth_option)
985 deallocate(is_smooth_passes)
986 deallocate(is_signed)
989 deallocate(is_missing_value)
990 deallocate(is_fill_missing)
993 call hash_destroy(bad_files)
995 end subroutine datalist_destroy
998 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
999 ! Name: reset_next_field
1001 ! Purpose: To reset the pointer to the next field in the list of fields
1002 ! specified by the user.
1003 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1004 subroutine reset_next_field()
1010 end subroutine reset_next_field
1013 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1014 ! Name: get_next_fieldname
1016 ! Purpose: Calling this routine results in field_name being set to the name of
1017 ! the field currently pointed to. If istatus /= 0 upon return, an error
1018 ! occurred, and the value of field_name is undefined.
1019 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1020 subroutine get_next_fieldname(field_name, istatus)
1025 integer, intent(out) :: istatus
1026 character (len=128), intent(out) :: field_name
1030 if (next_field <= num_entries) then
1032 field_name = source_fieldname(next_field)
1033 next_field = next_field + 1
1038 end subroutine get_next_fieldname
1041 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1042 ! Name: get_next_output_fieldname
1045 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1046 recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, &
1048 istagger, memorder, dimnames, units, &
1049 description, sr_x, sr_y, istatus)
1053 #include "wrf_io_flags.h"
1056 integer, intent(in) :: nest_num
1057 integer, intent(out) :: istatus, ndims, istagger, min_cat, max_cat
1058 integer, intent(out) :: sr_x, sr_y
1059 character (len=128), intent(out) :: memorder, field_name, units, description
1060 character (len=128), dimension(3), intent(out) :: dimnames
1063 integer :: landmask, temp_fieldtype
1064 logical :: is_water_mask, is_dom_only
1065 character (len=128) :: domcat_name, dfdx_name, dfdy_name
1066 character (len=256) :: temphash
1070 if (output_field_state == RETURN_LANDMASK) then
1071 call hash_init(duplicate_fnames)
1072 call get_landmask_field(field_name, is_water_mask, landmask, istatus)
1073 if (istatus == 0) then
1074 temphash(129:256) = ' '
1075 temphash(1:128) = field_name(1:128)
1076 call hash_insert(duplicate_fnames, temphash)
1077 call get_domcategory_name(field_name, domcat_name, is_dom_only, istatus)
1078 ! We will only save the dominant category
1079 if (is_dom_only .and. (istatus == 0)) then
1080 output_field_state = RETURN_DOMCAT_LM
1081 call get_next_output_fieldname(nest_num, field_name, ndims, &
1082 min_cat, max_cat, istagger, &
1083 memorder, dimnames, units, description, &
1084 sr_x, sr_y, istatus)
1090 temp_fieldtype = iget_fieldtype(field_name, istatus)
1091 if (istatus == 0) then
1092 if (temp_fieldtype == CONTINUOUS) then
1093 call get_max_levels(field_name, min_cat, max_cat, istatus)
1094 else if (temp_fieldtype == CATEGORICAL) then
1095 call get_max_categories(field_name, min_cat, max_cat, istatus)
1097 if (max_cat - min_cat > 0) ndims = 3
1099 call get_output_stagger(field_name, istagger, istatus)
1100 if (istagger == M) then
1101 dimnames(1) = 'west_east'
1102 dimnames(2) = 'south_north'
1103 else if (istagger == U) then
1104 dimnames(1) = 'west_east_stag'
1105 dimnames(2) = 'south_north'
1106 else if (istagger == V) then
1107 dimnames(1) = 'west_east'
1108 dimnames(2) = 'south_north_stag'
1109 else if (istagger == HH) then
1110 dimnames(1) = 'west_east'
1111 dimnames(2) = 'south_north'
1112 else if (istagger == VV) then
1113 dimnames(1) = 'west_east'
1114 dimnames(2) = 'south_north'
1116 if (ndims == 2) then
1119 else if (ndims == 3) then
1121 call get_z_dim_name(field_name, dimnames(3), istatus)
1127 call get_subgrid_dim_name(nest_num, field_name,dimnames(1:2), sr_x,sr_y,istatus)
1128 call get_source_units(field_name, 1, units, istatus)
1129 if (istatus /= 0) units = '-'
1130 call get_source_descr(field_name, 1, description, istatus)
1131 if (istatus /= 0) description = '-'
1133 output_field_state = RETURN_DOMCAT_LM
1136 output_field_state = RETURN_FIELDNAME
1137 call get_next_output_fieldname(nest_num, field_name, ndims, &
1138 min_cat, max_cat, istagger, &
1139 memorder, dimnames, units, description, &
1140 sr_x, sr_y, istatus)
1144 else if (output_field_state == RETURN_FIELDNAME) then
1145 call get_next_fieldname(field_name, istatus)
1146 temphash(129:256) = ' '
1147 temphash(1:128) = field_name(1:128)
1148 if (istatus == 0 .and. (.not. hash_search(duplicate_fnames, temphash))) then
1149 call hash_insert(duplicate_fnames, temphash)
1150 call get_domcategory_name(field_name, domcat_name, is_dom_only, istatus)
1151 ! We will only save the dominant category
1152 if (is_dom_only .and. (istatus == 0)) then
1153 output_field_state = RETURN_DOMCAT
1154 call get_next_output_fieldname(nest_num, field_name, ndims, &
1155 min_cat, max_cat, istagger, &
1156 memorder, dimnames, units, description, &
1157 sr_x, sr_y, istatus)
1160 ! Return the fractional field
1165 temp_fieldtype = iget_fieldtype(field_name, istatus)
1166 if (istatus == 0) then
1167 if (temp_fieldtype == CONTINUOUS) then
1168 call get_max_levels(field_name, min_cat, max_cat, istatus)
1169 else if (temp_fieldtype == CATEGORICAL) then
1170 call get_max_categories(field_name, min_cat, max_cat, istatus)
1172 if (max_cat - min_cat > 0) ndims = 3
1174 call get_output_stagger(field_name, istagger, istatus)
1175 if (istagger == M) then
1176 dimnames(1) = 'west_east'
1177 dimnames(2) = 'south_north'
1178 else if (istagger == U) then
1179 dimnames(1) = 'west_east_stag'
1180 dimnames(2) = 'south_north'
1181 else if (istagger == V) then
1182 dimnames(1) = 'west_east'
1183 dimnames(2) = 'south_north_stag'
1184 else if (istagger == HH) then
1185 dimnames(1) = 'west_east'
1186 dimnames(2) = 'south_north'
1187 else if (istagger == VV) then
1188 dimnames(1) = 'west_east'
1189 dimnames(2) = 'south_north'
1191 if (ndims == 2) then
1194 else if (ndims == 3) then
1196 call get_z_dim_name(field_name, dimnames(3), istatus)
1202 call get_subgrid_dim_name(nest_num,field_name,dimnames(1:2),sr_x,sr_y, istatus)
1203 call get_source_units(field_name, 1, units, istatus)
1204 if (istatus /= 0) units = '-'
1205 call get_source_descr(field_name, 1, description, istatus)
1206 if (istatus /= 0) description = '-'
1208 output_field_state = RETURN_DOMCAT
1210 else if (istatus /= 0) then
1211 output_field_state = RETURN_LANDMASK
1212 call hash_destroy(duplicate_fnames)
1214 else if (hash_search(duplicate_fnames, temphash)) then
1215 call get_next_output_fieldname(nest_num,field_name, ndims, &
1216 min_cat, max_cat, istagger, &
1217 memorder, dimnames, units, description, &
1218 sr_x, sr_y, istatus)
1222 else if (output_field_state == RETURN_DOMCAT .or. &
1223 output_field_state == RETURN_DOMCAT_LM ) then
1224 if (output_field_state == RETURN_DOMCAT) then
1225 next_field = next_field - 1
1226 call get_next_fieldname(field_name, istatus)
1228 call get_landmask_field(field_name, is_water_mask, landmask, istatus)
1230 if (istatus == 0) then
1231 call get_domcategory_name(field_name, domcat_name, is_dom_only, istatus)
1232 if (istatus == 0) then
1236 call get_output_stagger(field_name, istagger, istatus)
1237 if (istagger == M) then
1238 dimnames(1) = 'west_east'
1239 dimnames(2) = 'south_north'
1240 else if (istagger == U) then
1241 dimnames(1) = 'west_east_stag'
1242 dimnames(2) = 'south_north'
1243 else if (istagger == V) then
1244 dimnames(1) = 'west_east'
1245 dimnames(2) = 'south_north_stag'
1246 else if (istagger == HH) then
1247 dimnames(1) = 'west_east'
1248 dimnames(2) = 'south_north'
1249 else if (istagger == VV) then
1250 dimnames(1) = 'west_east'
1251 dimnames(2) = 'south_north'
1256 call get_subgrid_dim_name(nest_num,field_name,dimnames(1:2),sr_x,sr_y, istatus)
1257 field_name = domcat_name
1259 description = 'Dominant category'
1260 if (output_field_state == RETURN_DOMCAT) then
1261 output_field_state = RETURN_DFDX
1263 output_field_state = RETURN_DFDX_LM
1266 if (output_field_state == RETURN_DOMCAT) then
1267 output_field_state = RETURN_DFDX
1269 output_field_state = RETURN_DFDX_LM
1271 call get_next_output_fieldname(nest_num,field_name, ndims, &
1272 min_cat, max_cat, istagger, &
1273 memorder, dimnames, units, description, &
1274 sr_x, sr_y, istatus)
1278 call mprintf(.true., ERROR, 'get_next_output_fieldname(): In state DOMCAT, '// &
1279 'but no field name is found.')
1282 else if (output_field_state == RETURN_DFDX .or. &
1283 output_field_state == RETURN_DFDX_LM) then
1284 if (output_field_state == RETURN_DFDX) then
1285 next_field = next_field - 1
1286 call get_next_fieldname(field_name, istatus)
1288 call get_landmask_field(field_name, is_water_mask, landmask, istatus)
1290 if (istatus == 0) then
1291 call get_dfdx_name(field_name, dfdx_name, istatus)
1292 if (istatus == 0) then
1296 temp_fieldtype = iget_fieldtype(field_name, istatus)
1297 if (istatus == 0) then
1298 if (temp_fieldtype == CONTINUOUS) then
1299 call get_max_levels(field_name, min_cat, max_cat, istatus)
1300 else if (temp_fieldtype == CATEGORICAL) then
1301 call get_max_categories(field_name, min_cat, max_cat, istatus)
1303 if (max_cat - min_cat > 0) ndims = 3
1305 call get_output_stagger(field_name, istagger, istatus)
1306 if (istagger == M) then
1307 dimnames(1) = 'west_east'
1308 dimnames(2) = 'south_north'
1309 else if (istagger == U) then
1310 dimnames(1) = 'west_east_stag'
1311 dimnames(2) = 'south_north'
1312 else if (istagger == V) then
1313 dimnames(1) = 'west_east'
1314 dimnames(2) = 'south_north_stag'
1315 else if (istagger == HH) then
1316 dimnames(1) = 'west_east'
1317 dimnames(2) = 'south_north'
1318 else if (istagger == VV) then
1319 dimnames(1) = 'west_east'
1320 dimnames(2) = 'south_north'
1322 if (ndims == 2) then
1325 else if (ndims == 3) then
1327 call get_z_dim_name(field_name, dimnames(3), istatus)
1333 field_name = dfdx_name
1336 call get_subgrid_dim_name(nest_num,field_name,dimnames(1:2),sr_x,sr_y, istatus)
1337 description = 'df/dx'
1338 if (output_field_state == RETURN_DFDX) then
1339 output_field_state = RETURN_DFDY
1341 output_field_state = RETURN_DFDY_LM
1344 if (output_field_state == RETURN_DFDX) then
1345 output_field_state = RETURN_DFDY
1347 output_field_state = RETURN_DFDY_LM
1349 call get_next_output_fieldname(nest_num,field_name, ndims, &
1350 min_cat, max_cat, istagger, &
1351 memorder, dimnames, units, description, &
1352 sr_x, sr_y, istatus)
1356 call mprintf(.true., ERROR, 'get_next_output_fieldname(): In state DFDX, '// &
1357 'but no field name is found.')
1360 else if (output_field_state == RETURN_DFDY .or. &
1361 output_field_state == RETURN_DFDY_LM) then
1362 if (output_field_state == RETURN_DFDY) then
1363 next_field = next_field - 1
1364 call get_next_fieldname(field_name, istatus)
1366 call get_landmask_field(field_name, is_water_mask, landmask, istatus)
1368 if (istatus == 0) then
1369 call get_dfdy_name(field_name, dfdy_name, istatus)
1370 if (istatus == 0) then
1374 temp_fieldtype = iget_fieldtype(field_name, istatus)
1375 if (istatus == 0) then
1376 if (temp_fieldtype == CONTINUOUS) then
1377 call get_max_levels(field_name, min_cat, max_cat, istatus)
1378 else if (temp_fieldtype == CATEGORICAL) then
1379 call get_max_categories(field_name, min_cat, max_cat, istatus)
1381 if (max_cat - min_cat > 0) ndims = 3
1383 call get_output_stagger(field_name, istagger, istatus)
1384 if (istagger == M) then
1385 dimnames(1) = 'west_east'
1386 dimnames(2) = 'south_north'
1387 else if (istagger == U) then
1388 dimnames(1) = 'west_east_stag'
1389 dimnames(2) = 'south_north'
1390 else if (istagger == V) then
1391 dimnames(1) = 'west_east'
1392 dimnames(2) = 'south_north_stag'
1393 else if (istagger == HH) then
1394 dimnames(1) = 'west_east'
1395 dimnames(2) = 'south_north'
1396 else if (istagger == VV) then
1397 dimnames(1) = 'west_east'
1398 dimnames(2) = 'south_north'
1400 if (ndims == 2) then
1403 else if (ndims == 3) then
1405 call get_z_dim_name(field_name, dimnames(3), istatus)
1412 call get_subgrid_dim_name(nest_num,field_name,dimnames(1:2),sr_x,sr_y, istatus)
1413 field_name = dfdy_name
1415 description = 'df/dy'
1416 output_field_state = RETURN_FIELDNAME
1418 output_field_state = RETURN_FIELDNAME
1419 call get_next_output_fieldname(nest_num,field_name, ndims, &
1420 min_cat, max_cat, istagger, &
1421 memorder, dimnames, units, description, &
1422 sr_x, sr_y, istatus)
1426 call mprintf(.true., ERROR, 'get_next_output_fieldname(): In state DFDY, but no '// &
1427 'field name is found.')
1432 end subroutine get_next_output_fieldname
1435 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1436 ! Name: get_landmask_field
1438 ! Purpose: To return the name of the field from which the landmask is to be
1439 ! computed. If no error occurs, is_water_mask is .true. if the landmask
1440 ! value specifies the value of water, and .false. if the landmask value
1441 ! specifies the value of land.
1442 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1443 subroutine get_landmask_field(landmask_name, is_water_mask, landmask, istatus)
1448 integer, intent(out) :: landmask, istatus
1449 logical, intent(out) :: is_water_mask
1450 character (len=128), intent(out) :: landmask_name
1457 do idx=1,num_entries
1459 if (is_landmask_water(idx)) then
1460 is_water_mask = .true.
1461 landmask = source_landmask_water(idx)
1462 landmask_name = source_fieldname(idx)
1466 else if (is_landmask_land(idx)) then
1467 is_water_mask = .false.
1468 landmask = source_landmask_land(idx)
1469 landmask_name = source_fieldname(idx)
1477 end subroutine get_landmask_field
1480 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1481 ! Name: get_missing_value
1483 ! Pupose: Return the value used in the source data to indicate missing data.
1484 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1485 subroutine get_missing_value(fieldnm, ilevel, rmissing, istatus)
1490 integer, intent(in) :: ilevel
1491 integer, intent(out) :: istatus
1492 real, intent(out) :: rmissing
1493 character (len=128), intent(in) :: fieldnm
1500 do idx=1,num_entries
1501 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
1502 (len_trim(source_fieldname(idx)) == len_trim(fieldnm)) .and. &
1503 (source_priority(idx) == ilevel)) then
1505 if (is_missing_value(idx)) then
1506 rmissing = source_missing_value(idx)
1514 end subroutine get_missing_value
1517 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1518 ! Name: get_source_units
1520 ! Pupose: Return a string giving the units of the specified source data.
1521 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1522 subroutine get_source_units(fieldnm, ilevel, cunits, istatus)
1527 integer, intent(in) :: ilevel
1528 integer, intent(out) :: istatus
1529 character (len=128), intent(in) :: fieldnm
1530 character (len=128), intent(out) :: cunits
1537 do idx=1,num_entries
1538 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
1539 (len_trim(source_fieldname(idx)) == len_trim(fieldnm)) .and. &
1540 (source_priority(idx) == ilevel)) then
1542 if (is_units(idx)) then
1543 cunits = source_units(idx)
1551 end subroutine get_source_units
1554 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1555 ! Name: get_source_descr
1557 ! Pupose: Return a string giving a description of the specified source data.
1558 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1559 subroutine get_source_descr(fieldnm, ilevel, descr, istatus)
1564 integer, intent(in) :: ilevel
1565 integer, intent(out) :: istatus
1566 character (len=128), intent(in) :: fieldnm
1567 character (len=128), intent(out) :: descr
1574 do idx=1,num_entries
1575 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
1576 (len_trim(source_fieldname(idx)) == len_trim(fieldnm)) .and. &
1577 (source_priority(idx) == ilevel)) then
1579 if (is_units(idx)) then
1580 descr = source_descr(idx)
1588 end subroutine get_source_descr
1591 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1592 ! Name: get_missing_fill_value
1594 ! Pupose: Return the value to fill missing points with.
1595 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1596 subroutine get_missing_fill_value(fieldnm, rmissing, istatus)
1601 integer, intent(out) :: istatus
1602 real, intent(out) :: rmissing
1603 character (len=128), intent(in) :: fieldnm
1610 do idx=1,num_entries
1611 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
1612 (len_trim(source_fieldname(idx)) == len_trim(fieldnm)) ) then
1614 if (is_fill_missing(idx)) then
1615 rmissing = source_fill_missing(idx)
1623 end subroutine get_missing_fill_value
1626 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1627 ! Name: get_halt_on_missing
1629 ! Pupose: Returns 1 if the program should halt upon encountering a missing
1630 ! value in the final output field, and 0 otherwise.
1631 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1632 subroutine get_halt_on_missing(fieldnm, halt, istatus)
1637 integer, intent(out) :: istatus
1638 logical, intent(out) :: halt
1639 character (len=128), intent(in) :: fieldnm
1647 do idx=1,num_entries
1648 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
1649 (len_trim(source_fieldname(idx)) == len_trim(fieldnm)) ) then
1651 if (is_halt_missing(idx)) halt = .true.
1656 end subroutine get_halt_on_missing
1659 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1660 ! Name: get_masked_value
1662 ! Pupose: If the field is to be masked by the landmask, returns 0 if the field
1663 ! is masked over water and 1 if the field is masked over land. If no mask is
1664 ! to be applied, -1 is returned.
1665 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1666 subroutine get_masked_value(fieldnm, ilevel, masked, istatus)
1671 integer, intent(in) :: ilevel
1672 integer, intent(out) :: istatus
1673 real, intent(out) :: masked
1674 character (len=128), intent(in) :: fieldnm
1682 do idx=1,num_entries
1683 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
1684 (len_trim(source_fieldname(idx)) == len_trim(fieldnm)) .and. &
1685 (source_priority(idx) == ilevel)) then
1687 if (is_masked(idx)) then
1688 masked = source_masked(idx)
1695 end subroutine get_masked_value
1698 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1699 ! Name: get_max_levels
1701 ! Purpose: Returns the number of levels for the field given by fieldnm.
1702 ! The number of levels will generally be specified for continuous fields,
1703 ! whereas min/max category will generally be specified for categorical ones.
1704 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1705 subroutine get_max_levels(fieldnm, min_level, max_level, istatus)
1710 integer, intent(out) :: min_level, max_level, istatus
1711 character (len=128), intent(in) :: fieldnm
1715 logical :: have_found_field
1717 have_found_field = .false.
1720 do idx=1,num_entries
1721 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
1722 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
1724 if (is_dest_fieldtype(idx) .and. (source_dest_fieldtype(idx) /= CONTINUOUS)) then
1725 call mprintf(.true., WARN, 'In GEOGRID.TBL, destination field type for %s is '// &
1726 'not continuous and min/max levels specified.', s1=trim(fieldnm))
1728 if (.not. have_found_field) then
1729 if (is_tile_z_start(idx) .and. is_tile_z_end(idx)) then
1730 have_found_field = .true.
1732 min_level = source_tile_z_start(idx)
1733 max_level = source_tile_z_end(idx)
1734 else if (is_tile_z(idx)) then
1735 have_found_field = .true.
1738 max_level = source_tile_z(idx)
1741 if (.not. (is_tile_z_start(idx) .and. is_tile_z_end(idx))) then
1742 if (is_tile_z_start(idx) .or. is_tile_z_end(idx)) then
1743 call mprintf(.true., ERROR, 'In GEOGRID.TBL, only one of tile_z_start '// &
1744 'and tile_z_end specified for entry %i.',i1=idx)
1748 if (is_tile_z_start(idx) .and. is_tile_z_end(idx)) then
1749 if (source_tile_z_start(idx) < min_level) min_level = source_tile_z_start(idx)
1750 if (source_tile_z_end(idx) > max_level) max_level = source_tile_z_end(idx)
1751 else if (is_tile_z(idx)) then
1752 if (min_level > 1) min_level = 1
1753 if (source_tile_z(idx) > max_level) max_level = source_tile_z(idx)
1756 if (.not. (is_tile_z_start(idx) .and. is_tile_z_end(idx))) then
1757 if (is_tile_z_start(idx) .or. is_tile_z_end(idx)) then
1758 call mprintf(.true., ERROR, 'In GEOGRID.TBL, only one of tile_z_start '// &
1759 'and tile_z_end specified for entry %i.',i1=idx)
1767 end subroutine get_max_levels
1770 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1771 ! Name: get_source_levels
1773 ! Purpose: Return the min and max z-index for the source data for fieldname
1774 ! at a specified priority level (compared with the min/max level over
1775 ! all priority levels, as given by get_max_levels).
1776 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1777 subroutine get_source_levels(fieldnm, ilevel, min_level, max_level, istatus)
1782 integer, intent(in) :: ilevel
1783 integer, intent(out) :: min_level, max_level, istatus
1784 character (len=128), intent(in) :: fieldnm
1791 do idx=1,num_entries
1792 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
1793 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
1794 if (ilevel == source_priority(idx)) then
1796 if (is_tile_z_start(idx) .and. is_tile_z_end(idx)) then
1798 min_level = source_tile_z_start(idx)
1799 max_level = source_tile_z_end(idx)
1800 else if (is_tile_z(idx)) then
1803 max_level = source_tile_z(idx)
1806 if (.not. (is_tile_z_start(idx) .and. is_tile_z_end(idx))) then
1807 if (is_tile_z_start(idx) .or. is_tile_z_end(idx)) then
1808 call mprintf(.true., ERROR, 'In GEOGRID.TBL, only one of tile_z_start '// &
1809 'and tile_z_end specified for entry %i.',i1=idx)
1817 end subroutine get_source_levels
1820 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1821 ! Name: get_max_categories
1823 ! Purpose: Returns the minimum category and the maximum category for the field
1825 ! Min/max category will generally be specified for categorical fields,
1826 ! whereas the number of levels will generally be specified for continuous
1828 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1829 subroutine get_max_categories(fieldnm, min_category, max_category, istatus)
1834 integer, intent(out) :: min_category, max_category, istatus
1835 character (len=128), intent(in) :: fieldnm
1839 logical :: have_found_field
1841 have_found_field = .false.
1844 do idx=1,num_entries
1845 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
1846 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
1848 if (is_dest_fieldtype(idx) .and. (source_dest_fieldtype(idx) /= CATEGORICAL)) then
1849 call mprintf(.true., WARN, &
1850 'In GEOGRID.TBL, cannot get min/max categories for continuous '// &
1851 'field %s at entry %i. Perhaps the user has requested to '// &
1852 'perform a strange operation on the field.', s1=trim(fieldnm), i1=idx)
1854 if (.not. have_found_field) then
1855 if (is_category_min(idx) .and. is_category_max(idx)) then
1856 have_found_field = .true.
1858 min_category = source_category_min(idx)
1859 max_category = source_category_max(idx)
1860 else if (is_category_min(idx) .or. is_category_max(idx)) then
1861 call mprintf(.true., ERROR, 'In GEOGRID.TBL, only one of min_category and '// &
1862 'max_category specified for entry %i.',i1=idx)
1865 if (is_category_min(idx) .and. is_category_max(idx)) then
1866 if (source_category_min(idx) < min_category) min_category = source_category_min(idx)
1867 if (source_category_max(idx) > max_category) max_category = source_category_max(idx)
1868 else if (is_category_min(idx) .or. is_category_max(idx)) then
1869 call mprintf(.true., ERROR, 'In GEOGRID.TBL, only one of min_category and '// &
1870 'max_category specified for entry %i.',i1=idx)
1877 end subroutine get_max_categories
1880 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1881 ! Name: get_source_categories
1883 ! Purpose: Return the min and max category for the source data for fieldname
1884 ! at a specified priority level (compared with the min/max category over
1885 ! all priority levels, as given by get_max_categories).
1886 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1887 subroutine get_source_categories(fieldnm, ilevel, min_category, max_category, istatus)
1892 integer, intent(in) :: ilevel
1893 integer, intent(out) :: min_category, max_category, istatus
1894 character (len=128), intent(in) :: fieldnm
1901 do idx=1,num_entries
1902 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
1903 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
1904 if (ilevel == source_priority(idx)) then
1906 if (is_category_min(idx) .and. is_category_max(idx)) then
1908 min_category = source_category_min(idx)
1909 max_category = source_category_max(idx)
1910 else if (is_category_min(idx) .or. is_category_max(idx)) then
1911 call mprintf(.true., ERROR, 'In GEOGRID.TBL, only one of min_category '// &
1912 'and max_category specified for entry %i.',i1=idx)
1919 end subroutine get_source_categories
1922 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1923 ! Name: get_domcategory_name
1926 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1927 subroutine get_domcategory_name(fieldnm, domcat_name, ldominant_only, istatus)
1932 integer, intent(out) :: istatus
1933 logical, intent(out) :: ldominant_only
1934 character (len=128), intent(in) :: fieldnm
1935 character (len=128), intent(out) :: domcat_name
1941 ldominant_only = .false.
1943 do idx=1,num_entries
1944 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
1945 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
1947 if (is_dominant_category(idx)) then
1948 domcat_name = source_dominant_category(idx)
1950 if (is_dominant_only(idx)) then
1951 call mprintf(.true., WARN, 'In GEOGRID.TBL, both dominant_category and '// &
1952 'dominant_only are specified in entry %i. Using specification '// &
1953 'for dominant_category.',i1=idx)
1954 is_dominant_only(idx) = .false.
1958 else if (is_dominant_only(idx)) then
1959 domcat_name = source_dominant_only(idx)
1960 ldominant_only = .true.
1968 end subroutine get_domcategory_name
1971 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1972 ! Name: get_dfdx_name
1975 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1976 subroutine get_dfdx_name(fieldnm, dfdx_name, istatus)
1981 integer, intent(out) :: istatus
1982 character (len=128), intent(in) :: fieldnm
1983 character (len=128), intent(out) :: dfdx_name
1990 do idx=1,num_entries
1991 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
1992 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
1994 if (is_dfdx(idx)) then
1995 dfdx_name = source_dfdx(idx)
2003 end subroutine get_dfdx_name
2006 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2007 ! Name: get_dfdy_name
2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2011 subroutine get_dfdy_name(fieldnm, dfdy_name, istatus)
2016 integer, intent(out) :: istatus
2017 character (len=128), intent(in) :: fieldnm
2018 character (len=128), intent(out) :: dfdy_name
2025 do idx=1,num_entries
2026 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
2027 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
2029 if (is_dfdy(idx)) then
2030 dfdy_name = source_dfdy(idx)
2038 end subroutine get_dfdy_name
2041 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2042 ! Name: get_z_dim_name
2045 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2046 subroutine get_z_dim_name(fieldnm, z_dim, istatus)
2051 integer, intent(out) :: istatus
2052 character (len=128), intent(in) :: fieldnm
2053 character (len=128), intent(out) :: z_dim
2060 do idx=1,num_entries
2061 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
2062 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
2063 if (is_z_dim_name(idx)) then
2064 z_dim = source_z_dim_name(idx)
2071 end subroutine get_z_dim_name
2074 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2075 ! Name: get_field_scale_factor
2078 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2079 subroutine get_field_scale_factor(fieldnm, ilevel, scale_factor, istatus)
2084 integer, intent(in) :: ilevel
2085 integer, intent(out) :: istatus
2086 real, intent(out) :: scale_factor
2087 character (len=128), intent(in) :: fieldnm
2094 do idx=1,num_entries
2095 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
2096 (len_trim(source_fieldname(idx)) == len_trim(fieldnm)) .and. &
2097 (ilevel == source_priority(idx))) then
2099 if (is_scale_factor(idx)) then
2100 scale_factor = source_scale_factor(idx)
2107 end subroutine get_field_scale_factor
2110 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2111 ! Name: get_output_stagger
2114 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2115 subroutine get_output_stagger(fieldnm, istagger, istatus)
2122 integer, intent(out) :: istatus, istagger
2123 character (len=128), intent(in) :: fieldnm
2129 do idx=1,num_entries
2130 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
2131 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
2133 if (is_output_stagger(idx)) then
2135 istagger = source_output_stagger(idx)
2138 if (gridtype == 'C') then
2142 else if (gridtype == 'E') then
2152 end subroutine get_output_stagger
2155 subroutine get_subgrid_dim_name(nest_num,field_name,dimnames, &
2156 sub_x,sub_y,istatus)
2160 integer,intent(in)::nest_num
2161 integer,intent(out)::sub_x,sub_y,istatus
2162 character(len=128),intent(in)::field_name
2163 character(len=128),dimension(2),intent(inout)::dimnames
2170 do idx=1,num_entries
2171 if ((index(source_fieldname(idx),trim(field_name)) /= 0) .and. &
2172 (len_trim(source_fieldname(idx)) == len_trim(field_name))) then
2173 if (use_subgrid(idx)) then
2175 if(is_output_stagger(idx))then
2176 call mprintf(.true.,ERROR,'Cannot use subgrids on variables with staggered grids')
2178 dimnames(1)=trim(dimnames(1))//"_subgrid"
2179 dimnames(2)=trim(dimnames(2))//"_subgrid"
2180 sub_x=sr_x(nest_num)
2181 sub_y=sr_y(nest_num)
2185 end subroutine get_subgrid_dim_name
2188 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2189 ! Name: get_interp_option
2192 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2193 subroutine get_interp_option(fieldnm, ilevel, interp_opt, istatus)
2198 integer, intent(in) :: ilevel
2199 integer, intent(out) :: istatus
2200 character (len=128), intent(in) :: fieldnm
2201 character (len=128), intent(out) :: interp_opt
2207 do idx=1,num_entries
2208 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
2209 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
2210 if (ilevel == source_priority(idx)) then
2212 interp_opt = source_interp_string(idx)
2220 end subroutine get_interp_option
2223 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2224 ! Name: get_gcel_threshold
2227 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2228 subroutine get_gcell_threshold(interp_opt, threshold, istatus)
2233 integer, intent(out) :: istatus
2234 real, intent(out) :: threshold
2235 character (len=128), intent(in) :: interp_opt
2238 integer :: i, p1, p2
2243 i = index(interp_opt,'average_gcell')
2245 ! Move the "average_gcell" option to the beginning
2250 ! Check for a threshold
2251 p1 = index(interp_opt(i:128),'(')
2252 p2 = index(interp_opt(i:128),')')
2253 if (p1 /= 0 .and. p2 /= 0) then
2254 read(interp_opt(p1+1:p2-1),*,err=1000) threshold
2256 call mprintf(.true., WARN, 'Problem with specified threshold '// &
2257 'for average_gcell interp option. Setting threshold to 0.0.')
2265 1000 call mprintf(.true., ERROR, 'Threshold option to average_gcell interpolator '// &
2266 'must be a real number, enclosed in parentheses immediately '// &
2267 'after keyword "average_gcell"')
2269 end subroutine get_gcell_threshold
2272 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2273 ! Name: get_smooth_option
2276 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2277 subroutine get_smooth_option(fieldnm, smooth_opt, smooth_passes, istatus)
2282 integer, intent(out) :: istatus, smooth_opt, smooth_passes
2283 character (len=128), intent(in) :: fieldnm
2290 do idx=1,num_entries
2291 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
2292 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
2294 if (is_smooth_option(idx)) then
2296 smooth_opt = source_smooth_option(idx)
2298 if (is_smooth_passes(idx)) then
2299 smooth_passes = source_smooth_passes(idx)
2310 end subroutine get_smooth_option
2313 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2314 ! Name: iget_fieldtype
2317 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2318 function iget_fieldtype(fieldnm, istatus)
2323 integer, intent(out) :: istatus
2324 character (len=128), intent(in) :: fieldnm
2330 integer :: iget_fieldtype
2334 do idx=1,num_entries
2335 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
2336 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
2338 if (is_dest_fieldtype(idx)) then
2339 iget_fieldtype = source_dest_fieldtype(idx)
2347 end function iget_fieldtype
2350 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2351 ! Name: iget_source_fieldtype
2353 ! Purpose: Given a resolution, in degrees, and the name of a field, this
2354 ! function returns the type (categorical, continuous, etc.) of the source
2355 ! data that will be used. This may, in general, depend on the field name
2356 ! and the resolution; for example, near 30 second resolution, land use data
2357 ! may come from a categorical field, whereas for lower resolutions, it may
2358 ! come from arrays of land use fractions for each category.
2359 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2360 function iget_source_fieldtype(fieldnm, ilevel, istatus)
2365 integer, intent(in) :: ilevel
2366 integer, intent(out) :: istatus
2367 character (len=128), intent(in) :: fieldnm
2370 integer :: iget_source_fieldtype
2375 ! Find information about the source tiles for the specified fieldname
2377 do idx=1,num_entries
2378 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
2379 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
2380 if (ilevel == source_priority(idx)) then
2382 iget_source_fieldtype = source_fieldtype(idx)
2388 end function iget_source_fieldtype
2391 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2392 ! Name: get_data_tile
2395 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2396 subroutine get_data_tile(xlat, xlon, ilevel, field_name, &
2397 file_name, array, start_x_dim, end_x_dim, start_y_dim, &
2398 end_y_dim, start_z_dim, end_z_dim, npts_bdr, &
2404 integer, intent(in) :: ilevel
2405 integer, intent(out) :: istatus
2406 integer, intent(out) :: start_x_dim, end_x_dim, &
2407 start_y_dim, end_y_dim, &
2408 start_z_dim, end_z_dim, &
2410 real, intent(in) :: xlat, xlon ! Location that tile should contain
2411 real, pointer, dimension(:,:,:) :: array ! The array to be allocated by this routine
2412 character (len=128), intent(in) :: field_name
2413 character (len=256), intent(out) :: file_name
2417 integer :: local_wordsize, local_endian, sign_convention, irow_order, strlen
2418 integer :: xdim,ydim,zdim
2420 real, allocatable, dimension(:) :: temprow
2422 call get_tile_fname(file_name, xlat, xlon, ilevel, field_name, istatus)
2424 if (index(file_name, 'OUTSIDE') /= 0) then
2427 else if (hash_search(bad_files, file_name)) then
2432 call get_tile_dimensions(xlat, xlon, start_x_dim, end_x_dim, start_y_dim, end_y_dim, &
2433 start_z_dim, end_z_dim, npts_bdr, local_wordsize, local_endian, &
2434 sign_convention, ilevel, field_name, istatus)
2436 xdim = (end_x_dim-start_x_dim+1)
2437 ydim = (end_y_dim-start_y_dim+1)
2438 zdim = (end_z_dim-start_z_dim+1)
2440 if (associated(array)) deallocate(array)
2441 allocate(array(xdim,ydim,zdim))
2443 call get_row_order(field_name, ilevel, irow_order, istatus)
2444 if (istatus /= 0) irow_order = BOTTOM_TOP
2446 call s_len(file_name,strlen)
2450 call read_geogrid(file_name, strlen, array, xdim, ydim, zdim, &
2451 sign_convention, local_endian, scalefac, local_wordsize, istatus)
2453 if (irow_order == TOP_BOTTOM) then
2454 allocate(temprow(xdim))
2457 if (ydim-j+1 <= j) exit
2458 temprow(1:xdim) = array(1:xdim,j,k)
2459 array(1:xdim,j,k) = array(1:xdim,ydim-j+1,k)
2460 array(1:xdim,ydim-j+1,k) = temprow(1:xdim)
2466 if (istatus /= 0) then
2467 start_x_dim = INVALID
2468 start_y_dim = INVALID
2472 call hash_insert(bad_files, file_name)
2475 end subroutine get_data_tile
2478 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2479 ! Name: get_row_order
2480 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2481 subroutine get_row_order(fieldnm, ilevel, irow_order, istatus)
2486 integer, intent(in) :: ilevel
2487 character (len=128), intent(in) :: fieldnm
2488 integer, intent(out) :: irow_order, istatus
2494 do idx=1,num_entries
2495 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
2496 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
2497 if (ilevel == source_priority(idx)) then
2498 if (is_row_order(idx)) then
2499 irow_order = source_row_order(idx)
2507 end subroutine get_row_order
2510 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2511 ! Name: get_tile_dimensions
2512 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2513 subroutine get_tile_dimensions(xlat, xlon, start_x_dim, end_x_dim, start_y_dim, end_y_dim, &
2514 start_z_dim, end_z_dim, npts_bdr, bytes_per_datum, endianness, &
2515 sign_convention, ilevel, fieldnm, istatus)
2522 integer, intent(in) :: ilevel
2523 integer, intent(out) :: start_x_dim, end_x_dim, &
2524 start_y_dim, end_y_dim, &
2525 start_z_dim, end_z_dim, &
2527 bytes_per_datum, endianness, &
2528 sign_convention, istatus
2529 real, intent(in) :: xlat, xlon
2530 character (len=128), intent(in) :: fieldnm
2533 integer :: idx, current_domain
2537 do idx=1,num_entries
2538 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
2539 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
2540 if (ilevel == source_priority(idx)) then
2547 if (istatus /= 0) then
2558 current_domain = iget_selected_domain()
2559 call select_domain(SOURCE_PROJ)
2560 call lltoxy(xlat, xlon, rx, ry, M)
2561 call select_domain(current_domain)
2563 start_x_dim = source_tile_x(idx) * nint(real(floor((rx-0.5) / real(source_tile_x(idx))))) + 1
2564 end_x_dim = start_x_dim + source_tile_x(idx) - 1
2566 start_y_dim = source_tile_y(idx) * nint(real(floor((ry-0.5) / real(source_tile_y(idx))))) + 1
2567 end_y_dim = start_y_dim + source_tile_y(idx) - 1
2569 if (is_tile_z_start(idx) .and. is_tile_z_end(idx)) then
2570 start_z_dim = source_tile_z_start(idx)
2571 end_z_dim = source_tile_z_end(idx)
2572 else if (is_tile_z(idx)) then
2574 end_z_dim = source_tile_z(idx)
2577 if (.not. (is_tile_z_start(idx) .and. is_tile_z_end(idx))) then
2578 if (is_tile_z_start(idx) .or. is_tile_z_end(idx)) then
2579 call mprintf(.true., ERROR, 'In GEOGRID.TBL, only one of tile_z_start and '// &
2580 'tile_z_end specified for entry %i.',i1=idx)
2584 if (is_tile_bdr(idx)) then
2585 npts_bdr = source_tile_bdr(idx)
2590 start_x_dim = start_x_dim - npts_bdr
2591 end_x_dim = end_x_dim + npts_bdr
2592 start_y_dim = start_y_dim - npts_bdr
2593 end_y_dim = end_y_dim + npts_bdr
2595 if (is_wordsize(idx)) then
2596 bytes_per_datum = source_wordsize(idx)
2599 call mprintf(.true.,ERROR,'In GEOGRID.TBL, no wordsize specified for data in entry %i.',i1=idx)
2602 if (is_endian(idx)) then
2603 endianness = source_endian(idx)
2605 endianness = BIG_ENDIAN
2608 if (is_signed(idx)) then
2614 end subroutine get_tile_dimensions
2617 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2618 ! Name: get_tile_fname
2621 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2622 subroutine get_tile_fname(test_fname, xlat, xlon, ilevel, fieldname, istatus)
2630 integer, intent(in) :: ilevel
2631 integer, intent(out) :: istatus
2632 real, intent(in) :: xlat, xlon
2633 character (len=*), intent(in) :: fieldname
2634 character (len=256), intent(out) :: test_fname
2637 integer :: current_domain, idx
2641 write(test_fname, '(a)') 'TILE.OUTSIDE.DOMAIN'
2643 do idx=1,num_entries
2644 if ((index(source_fieldname(idx),trim(fieldname)) /= 0) .and. &
2645 (len_trim(source_fieldname(idx)) == len_trim(fieldname))) then
2646 if (ilevel == source_priority(idx)) then
2653 if (istatus /= 0) return
2655 current_domain = iget_selected_domain()
2656 call select_domain(SOURCE_PROJ)
2657 call lltoxy(xlat, xlon, rx, ry, M)
2658 call select_domain(current_domain)
2660 ! rx = real(source_tile_x(idx)) * real(floor((rx-0.5*source_dx(idx))/ real(source_tile_x(idx)))) + 1.0
2661 ! ry = real(source_tile_y(idx)) * real(floor((ry-0.5*source_dy(idx))/ real(source_tile_y(idx)))) + 1.0
2662 rx = real(source_tile_x(idx)) * real(floor((rx-0.5) / real(source_tile_x(idx)))) + 1.0
2663 ry = real(source_tile_y(idx)) * real(floor((ry-0.5) / real(source_tile_y(idx)))) + 1.0
2665 if (rx > 0. .and. ry > 0.) then
2666 write(test_fname, '(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(source_path(idx)), &
2667 nint(rx),'-',nint(rx)+source_tile_x(idx)-1,'.',nint(ry),'-',nint(ry)+source_tile_y(idx)-1
2670 end subroutine get_tile_fname
2673 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2674 ! Name: get_source_resolution
2677 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2678 subroutine get_source_resolution(fieldnm, ilevel, src_dx, src_dy, istatus)
2683 integer, intent(in) :: ilevel
2684 integer, intent(out) :: istatus
2685 real, intent(out) :: src_dx, src_dy
2686 character (len=*), intent(in) :: fieldnm
2692 do idx=1,num_entries
2693 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
2694 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
2695 if (ilevel == source_priority(idx)) then
2696 if (is_dx(idx) .and. is_dy(idx)) then
2697 src_dx = source_dx(idx)
2698 src_dy = source_dy(idx)
2699 if (source_proj(idx) /= PROJ_LATLON) then
2700 src_dx = src_dx / 111000.
2701 src_dy = src_dy / 111000.
2710 end subroutine get_source_resolution
2713 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2714 ! Name: get_data_projection
2716 ! Purpose: To acquire the parameters necessary in defining the grid on which
2717 ! the user-specified data for field 'fieldnm' are given.
2719 ! NOTES: If the routine successfully acquires values for all necessary
2720 ! parameters, istatus is set to 0. In case of an error,
2721 ! OR IF THE USER HAS NOT SPECIFIED A TILE OF DATA FOR FIELDNM,
2722 ! istatus is set to 1.
2723 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2724 subroutine get_data_projection(fieldnm, iproj, stand_lon, truelat1, truelat2, &
2725 dxkm, dykm, known_x, known_y, known_lat, known_lon, ilevel, istatus)
2730 integer, intent(in) :: ilevel
2731 integer, intent(out) :: iproj, istatus
2732 real, intent(out) :: stand_lon, truelat1, truelat2, dxkm, dykm, &
2733 known_x, known_y, known_lat, known_lon
2734 character (len=*), intent(in) :: fieldnm
2740 do idx=1,num_entries
2741 if ((index(source_fieldname(idx),trim(fieldnm)) /= 0) .and. &
2742 (len_trim(source_fieldname(idx)) == len_trim(fieldnm))) then
2743 if (ilevel == source_priority(idx)) then
2745 if (is_proj(idx)) then
2746 iproj = source_proj(idx)
2749 call mprintf(.true., ERROR, &
2750 'In GEOGRID.TBL, no specification for projection in entry %i.',i1=idx)
2752 if (is_known_x(idx)) then
2753 known_x = source_known_x(idx)
2756 call mprintf(.true., ERROR, &
2757 'In GEOGRID.TBL, no specification for known_x in entry %i.',i1=idx)
2759 if (is_known_y(idx)) then
2760 known_y = source_known_y(idx)
2763 call mprintf(.true., ERROR, &
2764 'In GEOGRID.TBL, no specification for known_y in entry %i.',i1=idx)
2766 if (is_known_lat(idx)) then
2767 known_lat = source_known_lat(idx)
2770 call mprintf(.true., ERROR, &
2771 'In GEOGRID.TBL, no specification for known_lat in entry %i.',i1=idx)
2773 if (is_known_lon(idx)) then
2774 known_lon = source_known_lon(idx)
2777 call mprintf(.true., ERROR, &
2778 'In GEOGRID.TBL, no specification for known_lon in entry %i.',i1=idx)
2780 if (is_truelat1(idx)) then
2781 truelat1 = source_truelat1(idx)
2782 else if (is_proj(idx) .and. source_proj(idx) /= PROJ_LATLON) then
2784 call mprintf(.true., WARN, &
2785 'In GEOGRID.TBL, no specification for truelat1 in entry %i.',i1=idx)
2787 if (is_truelat2(idx)) then
2788 truelat2 = source_truelat2(idx)
2789 else if (is_proj(idx) .and. source_proj(idx) /= PROJ_LATLON) then
2791 call mprintf(.true., WARN, &
2792 'In GEOGRID.TBL, no specification for truelat2 in entry %i.',i1=idx)
2794 if (is_stdlon(idx)) then
2795 stand_lon = source_stdlon(idx)
2796 else if (is_proj(idx) .and. source_proj(idx) /= PROJ_LATLON) then
2798 call mprintf(.true., WARN, &
2799 'In GEOGRID.TBL, no specification for stdlon in entry %i.',i1=idx)
2801 if (is_dx(idx)) then
2802 dxkm = source_dx(idx)
2805 call mprintf(.true., ERROR, &
2806 'In GEOGRID.TBL, no specification for dx in entry %i.',i1=idx)
2808 if (is_dy(idx)) then
2809 dykm = source_dy(idx)
2812 call mprintf(.true., ERROR, &
2813 'In GEOGRID.TBL, no specification for dy in entry %i.',i1=idx)
2820 end subroutine get_data_projection
2823 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2824 ! Name: check_data_specification
2826 ! Purpose: To check for obvious errors in the user source data specifications.
2827 ! Returns .true. if specification passes all checks, and .false. otherwise.
2828 ! For failing checks, diagnostic messages are printed.
2829 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2830 function check_data_specification( )
2835 logical :: check_data_specification
2838 integer :: i, j, istatus
2839 integer, pointer, dimension(:) :: priorities
2841 logical :: begin_priority, halt
2842 character (len=128) :: cur_name
2844 check_data_specification = .false.
2846 ! Check that each specification has a name, priority level, and path
2848 if (.not. is_fieldname(i)) then
2849 call mprintf(.true., ERROR, &
2850 'In GEOGRID.TBL, specification %i does not have a name.',i1=i)
2852 if (.not. is_priority(i)) then
2853 call mprintf(.true., ERROR, &
2854 'In GEOGRID.TBL, specification %i does not have a priority.',i1=i)
2856 if (list_length(source_res_path(i)) == 0) then
2857 call mprintf(.true., ERROR, &
2858 'In GEOGRID.TBL, no path (relative or absolute) is specified '// &
2859 'for entry %i.',i1=i)
2863 ! The fill_missing and halt_on_missing options are mutually exclusive
2865 call get_halt_on_missing(source_fieldname(i), halt, istatus)
2866 call get_missing_fill_value(source_fieldname(i), rmissing, istatus)
2867 if (halt .and. (istatus == 0)) then
2868 call mprintf(.true., ERROR, 'In GEOGRID.TBL, the halt_on_missing and fill_missing '// &
2869 'options are mutually exclusive, but both are given for field %s', &
2870 s1=trim(source_fieldname(i)))
2874 ! Check that the field from which landmask is calculated is not output on a staggering
2876 if (is_landmask_land(i) .or. is_landmask_water(i)) then
2877 if (is_output_stagger(i)) then
2878 if (source_output_stagger(i) /= M) then
2879 call mprintf(.true., ERROR, 'In GEOGRID.TBL, landmask cannot be derived from '// &
2880 'a field that is computed on a staggered grid at entry %i.',i1=i)
2886 ! Also check that any field that is to be masked by the landmask is not output on a staggering
2888 if (is_masked(i) .and. is_output_stagger(i)) then
2889 if (source_output_stagger(i) /= M) then
2890 call mprintf(.true., ERROR, 'In GEOGRID.TBL, landmask cannot be used with '// &
2891 'a field that is computed on a staggered grid at entry %i.',i1=i)
2896 allocate(priorities(num_entries))
2898 ! Now check that priorities for each source are unique and in the interval [1,n], n <= num_entries
2901 cur_name = source_fieldname(i)
2903 if (source_fieldname(j) == cur_name) then
2905 if (source_priority(j) > num_entries .or. source_priority(j) < 1) then
2906 call mprintf(.true., ERROR, 'In GEOGRID.TBL, priorities for %s do not '// &
2907 'form a sequence 1,2,...,n.', s1=trim(cur_name))
2910 if (priorities(source_priority(j)) == 1) then
2911 call mprintf(.true., ERROR, 'In GEOGRID.TBL, more than one entry for %s '// &
2912 'has priority %i.', s1=trim(cur_name), i1=source_priority(j))
2915 priorities(source_priority(j)) = 1
2922 begin_priority = .false.
2923 do j=num_entries,1,-1
2924 if (.not.begin_priority .and. priorities(j) == 1) then
2925 begin_priority = .true.
2926 else if (begin_priority .and. priorities(j) == 0) then
2927 call mprintf(.true., ERROR, 'In GEOGRID.TBL, no entry for %s has '// &
2928 'priority %i, but an entry has priority %i.', &
2929 s1=trim(cur_name), i1=j, i2=j+1)
2934 deallocate(priorities)
2936 ! Units must match for all priority levels of a field
2938 if (source_priority(i) == 1) then
2940 if ((source_fieldname(i) == source_fieldname(j)) .and. &
2941 (source_units(i) /= source_units(j))) then
2942 call mprintf(.true., ERROR, 'In GEOGRID.TBL, units for %s at entry %i '// &
2943 'do not match units at entry %i (%s)', &
2944 s1=trim(source_fieldname(i)), i1=j, i2=i, s2=trim(source_units(i)))
2950 ! Make sure that user has not asked to calculate landmask from a continuous field
2952 if (is_dest_fieldtype(i)) then
2953 if (source_dest_fieldtype(i) == CONTINUOUS) then
2954 if (is_landmask_water(i) .or. is_landmask_land(i)) then
2955 call mprintf(.true., ERROR, 'In GEOGRID.TBL, landmask cannot be '// &
2956 'calculated from a continuous destination field at entry %i.',i1=i)
2962 ! If either min_category or max_category is specified, then both must be specified
2964 if (is_category_min(i) .or. is_category_max(i)) then
2965 if (.not. is_category_min(i)) then
2966 call mprintf(.true., ERROR, 'In GEOGRID.TBL, for index file of data at '// &
2967 'entry %i, category_max is specified, but category_min is '// &
2968 'not. Both must be specified.',i1=i)
2969 else if (.not. is_category_max(i)) then
2970 call mprintf(.true., ERROR, 'In GEOGRID.TBL, for index file of data at '// &
2971 'entry %i, category_min is specified, but category_max is '// &
2972 'not. Both must be specified.',i1=i)
2977 ! For continuous data, (category_max - category_min + 1) should equal tile_z
2979 if (is_fieldtype(i)) then
2980 if (source_fieldtype(i) == CONTINUOUS) then
2981 if (is_category_max(i) .and. is_category_min(i) .and. is_tile_z(i)) then
2982 if (source_tile_z(i) /= (source_category_max(i) - source_category_min(i) + 1)) then
2983 call mprintf(.true., ERROR, 'In GEOGRID.TBL, tile_z must equal '// &
2984 '(category_max - category_min + 1) at entry %i.',i1=i)
2986 else if (is_category_max(i) .and. is_category_min(i) .and. &
2987 is_tile_z_start(i) .and. is_tile_z_end(i)) then
2988 if (source_tile_z_end(i) /= source_category_max(i) .or. &
2989 source_tile_z_start(i) /= source_category_min(i)) then
2990 call mprintf(.true., ERROR, 'In GEOGRID.TBL, tile_z_end must equal '// &
2991 'category_max, and tile_z_start must equal category_min '// &
2992 'at entry %i.',i1=i)
2999 ! Make sure that user has not named a dominant category or computed slope field
3000 ! the same as a fractional field
3002 if (source_dominant_category(i) == source_fieldname(i)) then
3003 call mprintf(.true., ERROR, 'In GEOGRID.TBL, dominant category cannot have '// &
3004 'the same name as the field at entry %i.',i1=i)
3008 if (.not. is_dominant_only(i)) then
3009 if (is_dfdx(j)) then
3010 if (source_dfdx(j) == source_fieldname(i)) then
3011 call mprintf(.true., ERROR, 'In GEOGRID.TBL, field name at entry %i '// &
3012 'cannot have the same name as the slope field df_dx at entry %i.', &
3016 if (is_dfdy(j)) then
3017 if (source_dfdy(j) == source_fieldname(i)) then
3018 call mprintf(.true., ERROR, 'In GEOGRID.TBL, field name at entry %i '// &
3019 'cannot have the same name as the slope field df_dy at entry %i.', &
3023 if (is_dfdx(j) .and. is_dominant_category(i)) then
3024 if (source_dfdx(j) == source_dominant_category(i)) then
3025 call mprintf(.true., ERROR, 'In GEOGRID.TBL, dominant field name at '// &
3026 'entry %i cannot have the same name as the slope field df_dx '// &
3027 'at entry %i.',i1=i, i2=j)
3030 if (is_dfdy(j) .and. is_dominant_category(i)) then
3031 if (source_dfdy(j) == source_dominant_category(i)) then
3032 call mprintf(.true., ERROR, 'In GEOGRID.TBL, dominant field name at '// &
3033 'entry %i cannot have the same name as the slope field '// &
3034 'df_dy at entry %i.',i1=i, i2=j)
3038 if (is_dfdx(j)) then
3039 if (source_dfdx(j) == source_dominant_only(i)) then
3040 call mprintf(.true., ERROR, 'In GEOGRID.TBL, dominant field name at '// &
3041 'entry %i cannot have the same name as the slope field '// &
3042 'df_dx at entry %i.',i1=i, i2=j)
3045 if (is_dfdy(j)) then
3046 if (source_dfdy(j) == source_dominant_only(i)) then
3047 call mprintf(.true., ERROR, 'In GEOGRID.TBL, dominant field name at '// &
3048 'entry %i cannot have the same name as the slope field '// &
3049 'df_dy at entry %i.',i1=i, i2=j)
3054 if (is_dfdx(i)) then
3055 if (is_dfdx(j)) then
3056 if (source_dfdx(j) == source_dfdx(i)) then
3057 call mprintf(.true., ERROR, 'In GEOGRID.TBL, slope field df_dx at '// &
3058 'entry %i cannot have the same name as the slope '// &
3059 'field df_dx at entry %i.',i1=i, i2=j)
3062 if (is_dfdy(j)) then
3063 if (source_dfdy(j) == source_dfdx(i)) then
3064 call mprintf(.true., ERROR, 'In GEOGRID.TBL, slope field df_dx at '// &
3065 'entry %i cannot have the same name as the slope field '// &
3066 'df_dy at entry %i.',i1=i, i2=j)
3070 if (is_dfdy(i)) then
3071 if (is_dfdx(j)) then
3072 if (source_dfdx(j) == source_dfdy(i)) then
3073 call mprintf(.true., ERROR, 'In GEOGRID.TBL, slope field df_dy at '// &
3074 'entry %i cannot have the same name as the slope field '// &
3075 'df_dx at entry %i.',i1=i, i2=j)
3078 if (is_dfdy(j)) then
3079 if (source_dfdy(j) == source_dfdy(i)) then
3080 call mprintf(.true., ERROR, 'In GEOGRID.TBL, slope field df_dy at '// &
3081 'entry %i cannot have the same name as the slope field '// &
3082 'df_dy at entry %i.',i1=i, i2=j)
3086 if (is_dominant_category(i)) then
3087 if (source_dominant_category(i) == source_fieldname(j)) then ! Possible exception
3088 if (.not. (is_dominant_only(j) .and. (source_dominant_category(i) /= source_dominant_only(j)))) then
3089 call mprintf(.true., ERROR, 'In GEOGRID.TBL, dominant category at '// &
3090 'entry %i cannot have the same name as the field at '// &
3091 'entry %i.',i1=i, i2=j)
3093 else if (is_dominant_category(j) .and. &
3094 (source_dominant_category(i) == source_dominant_category(j))) then
3095 call mprintf(.true., ERROR, 'In GEOGRID.TBL, dominant category at entry '// &
3096 '%i cannot have the same name as dominant category at '// &
3097 'entry %i.',i1=i, i2=j)
3098 else if (is_dominant_only(j) .and. &
3099 (source_dominant_category(i) == source_dominant_only(j))) then
3100 call mprintf(.true., ERROR, 'In GEOGRID.TBL, dominant category at '// &
3101 'entry %i cannot have the same name as dominant_only '// &
3102 'category at entry %i.',i1=i, i2=j)
3104 else if (is_dominant_only(i)) then
3105 if (source_dominant_only(i) == source_fieldname(j)) then ! Possible exception
3106 if (.not. (is_dominant_only(j) .and. (source_dominant_only(i) /= source_dominant_only(j)))) then
3107 call mprintf(.true., ERROR, 'In GEOGRID.TBL, dominant_only category '// &
3108 'at entry %i cannot have the same name as the field at '// &
3109 'entry %i.',i1=i, i2=j)
3111 else if (is_dominant_category(j) .and. &
3112 (source_dominant_only(i) == source_dominant_category(j))) then
3113 call mprintf(.true., ERROR, 'In GEOGRID.TBL, dominant_only category '// &
3114 'at entry %i cannot have the same name as dominant '// &
3115 'category at entry %i.',i1=i, i2=j)
3116 else if (is_dominant_only(j) .and. &
3117 (source_dominant_only(i) == source_dominant_only(j))) then
3118 call mprintf(.true., ERROR, 'In GEOGRID.TBL, dominant_only category '// &
3119 'at entry %i cannot have the same name as dominant_only '// &
3120 'category at entry %i.',i1=i, i2=j)
3127 check_data_specification = .true.
3128 end function check_data_specification
3131 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3134 ! Purpose: This routine receives a fortran string, and returns the number of
3135 ! characters in the string before the first "space" is encountered. It
3136 ! considers ascii characters 33 to 126 to be valid characters, and ascii
3137 ! 0 to 32, and 127 to be "space" characters.
3138 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3139 subroutine s_len(string, s_length)
3144 character (len=*), intent(in) :: string
3145 integer, intent(out) :: s_length
3148 integer :: i, len_str, aval
3153 len_str = len(string)
3155 do while ((i .le. len_str) .and. (.not. space))
3156 aval = ichar(string(i:i))
3157 if ((aval .lt. 33) .or. (aval .gt. 126)) then
3164 end subroutine s_len
3166 end module source_data_module