added README_changes.txt
[wrffire.git] / WPS / geogrid / src / source_data_module.F90
blobbcb7d20a65a1c8329065405c2a1c51ffc285fa12
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 ! Module: source_data_module
4 ! Description: 
5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
6 module source_data_module
8    use hash_module
9    use list_module
10    use module_debug
11    use misc_definitions_module
13    ! Parameters
14    integer, parameter :: RETURN_LANDMASK = 0, &
15                          RETURN_DOMCAT_LM = 1, &
16                          RETURN_DFDX_LM = 2, &
17                          RETURN_DFDY_LM = 3, &
18                          RETURN_FIELDNAME = 4, &
19                          RETURN_DOMCAT = 5, &
20                          RETURN_DFDX = 6, &
21                          RETURN_DFDY = 7
23    ! Module variables
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 
58    contains
61    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
62    ! Name: get_datalist
63    !
64    ! Purpose:
65    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
66    subroutine get_datalist()
68       use gridinfo_module
69   
70       implicit none
71   
72       ! Parameters
73       integer, parameter :: BUFSIZE = 256
74   
75       ! Local variables
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
80   
81       nparams = 0
82       num_entries = 0
83   
84       do funit=10,100
85          inquire(unit=funit, opened=is_used)
86          if (.not. is_used) exit
87       end do
88       open(funit,file=trim(opt_geogrid_tbl_path)//'GEOGRID.TBL',form='formatted',status='old',err=1000)
89   
90       !
91       ! We will first go through the file to determine how many field
92       !   specifications there are
93       !
94     10 read(funit,'(a)',end=20,err=1001) buffer
95       call despace(buffer)
96   
97       ! Is this line a comment?
98       if (buffer(1:1) == '#') then
99   
100       ! Are we beginning a new field specification?
101       else if (index(buffer,'=====') /= 0) then
102          if (nparams > 0) num_entries = num_entries + 1  
103          nparams = 0
104   
105       else
106          eos = index(buffer,'#')
107          if (eos /= 0) buffer(eos:BUFSIZE) = ' ' 
108    
109          ! Does this line contain at least one parameter specification?
110          if (index(buffer,'=') /= 0) then
111             nparams = nparams + 1
112          end if
113       end if
114       go to 10
115   
116     20 rewind(funit)
117   
118       ! In case the last entry ended without a ======== line
119       if (nparams > 0) num_entries = num_entries + 1  
120   
121       call mprintf(.true.,STDOUT,'Parsed %i entries in GEOGRID.TBL', i1=num_entries)
122   
123       !
124       ! Now that we know how many fields the user has specified, allocate
125       !   the properly sized arrays
126       !
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))
172       do i=1,num_entries
173          call list_init(source_res_path(i))
174          call list_init(source_interp_option(i))
175       end do
176   
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'
224       source_iswater = 16
225       source_isice = 24
226       source_isurban = 1
227       source_isoilwater = 14
228   
229       ! 
230       ! Actually read and save the specifications
231       !
232       nparams = 0
233       i = 1
234     30 buffer = ' '
235       read(funit,'(a)',end=40,err=1001) buffer
236       call despace(buffer)
237   
238       ! Is this line a comment?
239       if (buffer(1:1) == '#') then
240          ! Do nothing.
241   
242       ! Are we beginning a new field specification?
243       else if (index(buffer,'=====') /= 0) then   !{
244          if (nparams > 0) i = i + 1  
245          nparams = 0
246          if (i <= num_entries) then
247             is_path(i) = .false.
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.
255             is_dfdx(i) = .false.
256             is_dfdy(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.
264          end if
265   
266       else
267          ! Check whether the current line is a comment
268          if (buffer(1:1) /= '#') then
269             have_specification = .true. 
270          else
271             have_specification = .false.
272          end if
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) = ' ' 
277    
278          do while (have_specification)   !{
279    
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,';')
284             if (eos == 0) then
285                have_specification = .false.
286                eos = BUFSIZE
287             end if
288     
289             idx = index(buffer(1:eos-1),'=')
290     
291             if (idx /= 0) then   !{
292                nparams = nparams + 1
293      
294                if (index('name',trim(buffer(1:idx-1))) /= 0) then
295                   ispace = idx+1
296                   do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
297                      ispace = ispace + 1
298                   end do 
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)
303      
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)
307      
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
315                   end if
316      
317                else if (index('interp_option',trim(buffer(1:idx-1))) /= 0) then
318                   ispace = idx+1
319                   do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
320                      ispace = ispace + 1
321                   end do 
322                   interp_string = ' '
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)
327                   else
328                      res_string = 'default'
329                   end if
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))
336                   else
337                      call list_insert(source_interp_option(i), ckey=res_string, cvalue=interp_string)
338                   end if
339      
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
353                   end if
354      
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)
358        
359                else if (index('rel_path',trim(buffer(1:idx-1))) /= 0) then
360                   ispace = idx+1
361                   do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
362                      ispace = ispace + 1
363                   end do 
364                   path_string = ' '
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) = ' '
371                   end if
372                   ispace = index(path_string,':')
373                   if (ispace /= 0) then
374                      write(res_string,'(a)') path_string(1:ispace-1)
375                   else
376                      res_string = 'default'
377                   end if
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))
384                   else
385                      call list_insert(source_res_path(i), ckey=res_string, cvalue=path_string)
386                   end if
387      
388                else if (index('abs_path',trim(buffer(1:idx-1))) /= 0) then
389                   ispace = idx+1
390                   do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
391                      ispace = ispace + 1
392                   end do 
393                   path_string = ' '
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)
400                   else
401                      res_string = 'default'
402                   end if
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))
409                   else
410                      call list_insert(source_res_path(i), ckey=res_string, cvalue=path_string)
411                   end if
412     
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
429                   end if
430      
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)
435      
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)
440      
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.
449                   end if
450      
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)
454        
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.
460                   end if
461      
462                else if (index('dominant_category',trim(buffer(1:idx-1))) /= 0) then
463                   ispace = idx+1
464                   do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
465                      ispace = ispace + 1
466                   end do 
467                   is_dominant_category(i) = .true.
468                   source_dominant_category(i) = ' '
469                   source_dominant_category(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
470       
471                else if (index('dominant_only',trim(buffer(1:idx-1))) /= 0) then
472                   ispace = idx+1
473                   do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
474                      ispace = ispace + 1
475                   end do 
476                   is_dominant_only(i) = .true.
477                   source_dominant_only(i) = ' '
478                   source_dominant_only(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
479      
480                else if (index('df_dx',trim(buffer(1:idx-1))) /= 0) then
481                   ispace = idx+1
482                   do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
483                      ispace = ispace + 1
484                   end do 
485                   is_dfdx(i) = .true.
486                   source_dfdx(i) = ' '
487                   source_dfdx(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
488      
489                else if (index('df_dy',trim(buffer(1:idx-1))) /= 0) then
490                   ispace = idx+1
491                   do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
492                      ispace = ispace + 1
493                   end do 
494                   is_dfdy(i) = .true.
495                   source_dfdy(i) = ' '
496                   source_dfdy(i)(1:ispace-idx) = buffer(idx+1:ispace-1)
497      
498                else if (index('z_dim_name',trim(buffer(1:idx-1))) /= 0) then
499                   ispace = idx+1
500                   do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
501                      ispace = ispace + 1
502                   end do 
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
508                    ispace = idx+1
509                    do while ((ispace < eos) .and. (buffer(ispace:ispace) /= ' '))
510                      ispace = ispace +1
511                    end do
512                    use_subgrid(i)=.true.
513      
514                else
515                    call mprintf(.true., WARN, 'In GEOGRID.TBL, unrecognized option %s in '// &
516                                 'entry %i.',i1=idx, s1=buffer(i:idx-1))
517        
518                end if
519     
520             end if   !} index(buffer(1:eos-1),'=') /= 0
521     
522             buffer = buffer(eos+1:BUFSIZE)
523          end do   ! while eos /= 0 }
524   
525       end if   !} index(buffer, '=====') /= 0
526       go to 30
527   
528     40 close(funit)
529   
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.')
534       end if
535   
536       call hash_init(bad_files)
537   
538       return
539   
540     1000 call mprintf(.true.,ERROR,'Could not open GEOGRID.TBL')
541   
542     1001 call mprintf(.true.,ERROR,'Encountered error while reading GEOGRID.TBL')
544    end subroutine get_datalist
547    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
548    ! Name: get_source_params
549    !
550    ! Purpose: 
551    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
552    subroutine get_source_params(res_string)
555       implicit none
556   
557       ! Parameters
558       integer, parameter :: BUFSIZE = 256
559   
560       ! Arguments
561       character (len=128), intent(in) :: res_string
562   
563       ! Local variables
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
569   
570       do idx=1,num_entries
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.
576          is_dx(idx) = .false.
577          is_dy(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.
599    
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))
604             temp_str = 'default'
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)))
608             else
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.', &
611                             s1=trim(res_string))
612             end if
613          end if
614    
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))
618             temp_str = 'default'
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)))
622             else
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.', &
625                             s1=trim(res_string))
626             end if
627          end if
628          write(test_fname, '(a)') trim(source_path(idx))//'index'
629      
630          do funit=10,100
631             inquire(unit=funit, opened=is_used)
632             if (.not. is_used) exit
633          end do
634          open(funit,file=trim(test_fname),form='formatted',status='old',err=1000)
635    
636       30 buffer = ' '
637          read(funit,'(a)',end=40, err=1001) buffer
638          call despace(buffer)
639      
640          ! Is this line a comment?
641          if (buffer(1:1) == '#') then
642             ! Do nothing.
643      
644          else
645             have_specification = .true. 
646       
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) = ' ' 
650       
651             do while (have_specification)   !{
652       
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,';')
657                if (eos == 0) then
658                   have_specification = .false.
659                   eos = BUFSIZE
660                end if
661        
662                i = index(buffer(1:eos-1),'=')
663        
664                if (i /= 0) then   !{
665        
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
688                      end if
689         
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
697                      end if
698         
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.
704                      end if
705         
706                   else if (index('units',trim(buffer(1:i-1))) /= 0) then
707                      ispace = i+1
708                      iquoted = 0
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)
711                         ispace = ispace + 1
712                      end do 
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)
718         
719                   else if (index('description',trim(buffer(1:i-1))) /= 0) then
720                      ispace = i+1
721                      iquoted = 0
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)
724                         ispace = ispace + 1
725                      end do 
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)
731         
732                   else if (index('mminlu',trim(buffer(1:i-1))) /= 0) then
733                      ispace = i+1
734                      iquoted = 0
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)
737                         ispace = ispace + 1
738                      end do 
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)
742         
743                   else if (index('iswater',trim(buffer(1:i-1))) /= 0) then
744                      read(buffer(i+1:eos-1),*) source_iswater
745           
746                   else if (index('isice',trim(buffer(1:i-1))) /= 0) then
747                      read(buffer(i+1:eos-1),*) source_isice
748           
749                   else if (index('isurban',trim(buffer(1:i-1))) /= 0) then
750                      read(buffer(i+1:eos-1),*) source_isurban
751           
752                   else if (index('isoilwater',trim(buffer(1:i-1))) /= 0) then
753                      read(buffer(i+1:eos-1),*) source_isoilwater
754           
755                   else if (index('dx',trim(buffer(1:i-1))) /= 0) then
756                      is_dx(idx) = .true.
757                      read(buffer(i+1:eos-1),*) source_dx(idx)
758           
759                   else if (index('dy',trim(buffer(1:i-1))) /= 0) then
760                      is_dy(idx) = .true.
761                      read(buffer(i+1:eos-1),*) source_dy(idx)
762           
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)
766           
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)
770           
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)
774           
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)
778           
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)
782           
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)
786           
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)
790           
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)
794         
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
802                      else
803                         call mprintf(.true.,WARN,'Invalid value for keyword ''endian'' '// &
804                                      'specified in index file. BIG_ENDIAN will be used.')
805                      end if
806           
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
814                      end if
815           
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)
819         
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)
823         
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)
827         
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)
831         
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)
835         
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)
839         
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)
843         
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)
847         
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)
851         
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)
855           
856                   else
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)
859           
860                   end if
861        
862                end if   !} index(buffer(1:eos-1),'=') /= 0
863        
864                buffer = buffer(eos+1:BUFSIZE)
865             end do   ! while eos /= 0 }
866      
867          end if   !} index(buffer, '=====') /= 0
868      
869          go to 30
870     
871     40 close(funit)
872   
873       end do
874   
875       return
876   
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
882    
883    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
884    ! Name: datalist_destroy()
885    !
886    ! Purpose: This routine deallocates any memory that was allocated by the 
887    !   get_datalist() subroutine.
888    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
889    subroutine datalist_destroy()
891       implicit none
892   
893       ! Local variables
894       integer :: i
895   
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)
940          do i=1,num_entries
941             call list_destroy(source_res_path(i))
942             call list_destroy(source_interp_option(i))
943          end do
944          deallocate(source_res_path)
945          deallocate(source_interp_option)
946      
947          deallocate(is_wordsize)
948          deallocate(is_endian)
949          deallocate(is_fieldtype)
950          deallocate(is_dest_fieldtype)
951          deallocate(is_proj)
952          deallocate(is_priority)
953          deallocate(is_dx)
954          deallocate(is_dy)
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)
963          deallocate(is_path)
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)
980          deallocate(is_dfdx)
981          deallocate(is_dfdy)
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)
987          deallocate(is_units)
988          deallocate(is_descr)
989          deallocate(is_missing_value)
990          deallocate(is_fill_missing)
991       end if
992   
993       call hash_destroy(bad_files)
995    end subroutine datalist_destroy
998    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
999    ! Name: reset_next_field
1000    !
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()
1006       implicit none
1007   
1008       next_field = 1
1010    end subroutine reset_next_field
1013    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1014    ! Name: get_next_fieldname
1015    !
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)
1022       implicit none
1023   
1024       ! Arguments
1025       integer, intent(out) :: istatus
1026       character (len=128), intent(out) :: field_name
1027   
1028       istatus = 1
1029   
1030       if (next_field <= num_entries) then
1031   
1032          field_name = source_fieldname(next_field)
1033          next_field = next_field + 1
1034          istatus = 0
1035   
1036       end if
1037      
1038    end subroutine get_next_fieldname
1041    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1042    ! Name: get_next_output_fieldname
1043    !
1044    ! Purpose: 
1045    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1046    recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, &
1047                                             min_cat, max_cat, &
1048                                             istagger, memorder, dimnames, units, &
1049                                             description, sr_x, sr_y, istatus)
1051       implicit none
1052   
1053 #include "wrf_io_flags.h"
1054   
1055       ! Arguments
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
1061   
1062       ! Local variables
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
1067   
1068       istatus = 1
1069   
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)
1085                return
1086             else
1087                ndims = 2
1088                min_cat = 1
1089                max_cat = 1
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)
1096                   end if
1097                   if (max_cat - min_cat > 0) ndims = 3
1098                end if
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' 
1115                end if
1116                if (ndims == 2) then
1117                   memorder = 'XY ' 
1118                   dimnames(3) = ' '
1119                else if (ndims == 3) then
1120                   memorder = 'XYZ' 
1121                   call get_z_dim_name(field_name, dimnames(3), istatus)
1122                   istatus = 0 
1123                else
1124                   memorder = '   ' 
1125                   dimnames(3) = ' '
1126                end if
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 = '-'
1132                istatus = 0
1133                output_field_state = RETURN_DOMCAT_LM
1134             end if
1135          else
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)
1141             return
1142          end if
1143   
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)
1158                return
1159      
1160             ! Return the fractional field
1161             else
1162                ndims = 2
1163                min_cat = 1
1164                max_cat = 1
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)
1171                   end if
1172                   if (max_cat - min_cat > 0) ndims = 3
1173                end if
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' 
1190                end if
1191                if (ndims == 2) then
1192                   memorder = 'XY ' 
1193                   dimnames(3) = ' '
1194                else if (ndims == 3) then
1195                   memorder = 'XYZ' 
1196                   call get_z_dim_name(field_name, dimnames(3), istatus)
1197                   istatus = 0 
1198                else
1199                   memorder = '   ' 
1200                   dimnames(3) = ' '
1201                end if
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 = '-'
1207                istatus = 0
1208                output_field_state = RETURN_DOMCAT  
1209             end if 
1210          else if (istatus /= 0) then
1211             output_field_state = RETURN_LANDMASK
1212             call hash_destroy(duplicate_fnames)
1213             return 
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)
1219             return
1220          end if
1221   
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)
1227          else
1228             call get_landmask_field(field_name, is_water_mask, landmask, istatus)
1229          end if
1230          if (istatus == 0) then
1231             call get_domcategory_name(field_name, domcat_name, is_dom_only, istatus)
1232             if (istatus == 0) then
1233                ndims = 2
1234                min_cat = 1
1235                max_cat = 1
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' 
1252                end if
1253                dimnames(3) = ' ' 
1254                memorder = 'XY ' 
1256                call get_subgrid_dim_name(nest_num,field_name,dimnames(1:2),sr_x,sr_y, istatus)
1257                field_name = domcat_name
1258                units = 'category'
1259                description = 'Dominant category'
1260                if (output_field_state == RETURN_DOMCAT) then
1261                   output_field_state = RETURN_DFDX
1262                else
1263                   output_field_state = RETURN_DFDX_LM
1264                end if
1265             else
1266                if (output_field_state == RETURN_DOMCAT) then
1267                   output_field_state = RETURN_DFDX
1268                else
1269                   output_field_state = RETURN_DFDX_LM
1270                end if
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)
1275               return
1276             end if 
1277          else
1278             call mprintf(.true., ERROR, 'get_next_output_fieldname(): In state DOMCAT, '// &
1279                          'but no field name is found.')
1280          end if
1281   
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)
1287          else
1288             call get_landmask_field(field_name, is_water_mask, landmask, istatus)
1289          end if
1290          if (istatus == 0) then
1291             call get_dfdx_name(field_name, dfdx_name, istatus)
1292             if (istatus == 0) then
1293                ndims = 2
1294                min_cat = 1
1295                max_cat = 1
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)
1302                   end if
1303                   if (max_cat - min_cat > 0) ndims = 3
1304                end if
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' 
1321                end if
1322                if (ndims == 2) then
1323                   memorder = 'XY ' 
1324                   dimnames(3) = ' '
1325                else if (ndims == 3) then
1326                   memorder = 'XYZ' 
1327                   call get_z_dim_name(field_name, dimnames(3), istatus)
1328                   istatus = 0 
1329                else
1330                   memorder = '   ' 
1331                   dimnames(3) = ' '
1332                end if
1333                field_name = dfdx_name
1334                units = '-'
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
1340                else
1341                   output_field_state = RETURN_DFDY_LM
1342                end if
1343             else 
1344                if (output_field_state == RETURN_DFDX) then
1345                   output_field_state = RETURN_DFDY
1346                else
1347                   output_field_state = RETURN_DFDY_LM
1348                end if
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)
1353                return
1354             end if 
1355          else
1356             call mprintf(.true., ERROR, 'get_next_output_fieldname(): In state DFDX, '// &
1357                          'but no field name is found.')
1358          end if
1359   
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)
1365          else
1366             call get_landmask_field(field_name, is_water_mask, landmask, istatus)
1367          end if
1368          if (istatus == 0) then
1369             call get_dfdy_name(field_name, dfdy_name, istatus)
1370             if (istatus == 0) then
1371                ndims = 2
1372                min_cat = 1
1373                max_cat = 1
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)
1380                   end if
1381                   if (max_cat - min_cat > 0) ndims = 3
1382                end if
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' 
1399                end if
1400                if (ndims == 2) then
1401                   memorder = 'XY ' 
1402                   dimnames(3) = ' '
1403                else if (ndims == 3) then
1404                   memorder = 'XYZ' 
1405                   call get_z_dim_name(field_name, dimnames(3), istatus)
1406                   istatus = 0 
1407                else
1408                   memorder = '   ' 
1409                   dimnames(3) = ' '
1410                end if
1411                
1412                call get_subgrid_dim_name(nest_num,field_name,dimnames(1:2),sr_x,sr_y, istatus)
1413                field_name = dfdy_name
1414                units = '-'
1415                description = 'df/dy'
1416                output_field_state = RETURN_FIELDNAME
1417             else
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)
1423                return
1424             end if 
1425          else
1426             call mprintf(.true., ERROR, 'get_next_output_fieldname(): In state DFDY, but no '// &
1427                          'field name is found.')
1428          end if
1429   
1430       end if
1431   
1432    end subroutine get_next_output_fieldname
1435    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1436    ! Name: get_landmask_field
1437    !
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)
1445       implicit none
1446   
1447       ! Arguments
1448       integer, intent(out) :: landmask, istatus
1449       logical, intent(out) :: is_water_mask
1450       character (len=128), intent(out) :: landmask_name
1451   
1452       ! Local variables
1453       integer :: idx
1454   
1455       istatus = 1
1456   
1457       do idx=1,num_entries
1458   
1459          if (is_landmask_water(idx)) then
1460             is_water_mask = .true.
1461             landmask = source_landmask_water(idx) 
1462             landmask_name = source_fieldname(idx)
1463             istatus = 0
1464             exit
1465    
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)
1470             istatus = 0
1471             exit
1473          end if
1474   
1475       end do
1477    end subroutine get_landmask_field
1480    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1481    ! Name: get_missing_value
1482    !
1483    ! Pupose: Return the value used in the source data to indicate missing data.
1484    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1485    subroutine get_missing_value(fieldnm, ilevel, rmissing, istatus)
1487       implicit none
1488   
1489       ! Arguments
1490       integer, intent(in) :: ilevel
1491       integer, intent(out) :: istatus
1492       real, intent(out) :: rmissing
1493       character (len=128), intent(in) :: fieldnm
1494   
1495       ! Local variables
1496       integer :: idx
1497   
1498       istatus = 1
1499   
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
1504    
1505             if (is_missing_value(idx)) then
1506                rmissing = source_missing_value(idx)
1507                istatus = 0
1508                exit
1509             end if
1510    
1511          end if
1512       end do
1513   
1514    end subroutine get_missing_value
1517    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1518    ! Name: get_source_units
1519    !
1520    ! Pupose: Return a string giving the units of the specified source data.
1521    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1522    subroutine get_source_units(fieldnm, ilevel, cunits, istatus)
1524       implicit none
1525   
1526       ! Arguments
1527       integer, intent(in) :: ilevel
1528       integer, intent(out) :: istatus
1529       character (len=128), intent(in) :: fieldnm
1530       character (len=128), intent(out) :: cunits
1531   
1532       ! Local variables
1533       integer :: idx
1534   
1535       istatus = 1
1536   
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
1541    
1542             if (is_units(idx)) then
1543                cunits = source_units(idx)
1544                istatus = 0
1545                exit
1546             end if
1547    
1548          end if
1549       end do
1550   
1551    end subroutine get_source_units
1554    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1555    ! Name: get_source_descr
1556    !
1557    ! Pupose: Return a string giving a description of the specified source data.
1558    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1559    subroutine get_source_descr(fieldnm, ilevel, descr, istatus)
1561       implicit none
1562   
1563       ! Arguments
1564       integer, intent(in) :: ilevel
1565       integer, intent(out) :: istatus
1566       character (len=128), intent(in) :: fieldnm
1567       character (len=128), intent(out) :: descr
1568   
1569       ! Local variables
1570       integer :: idx
1571   
1572       istatus = 1
1573   
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
1578    
1579             if (is_units(idx)) then
1580                descr = source_descr(idx)
1581                istatus = 0
1582                exit
1583             end if
1584    
1585          end if
1586       end do
1588    end subroutine get_source_descr
1591    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1592    ! Name: get_missing_fill_value
1593    !
1594    ! Pupose: Return the value to fill missing points with.
1595    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1596    subroutine get_missing_fill_value(fieldnm, rmissing, istatus)
1598       implicit none
1599   
1600       ! Arguments
1601       integer, intent(out) :: istatus
1602       real, intent(out) :: rmissing
1603       character (len=128), intent(in) :: fieldnm
1604   
1605       ! Local variables
1606       integer :: idx
1607   
1608       istatus = 1
1609   
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 
1613    
1614             if (is_fill_missing(idx)) then
1615                rmissing = source_fill_missing(idx)
1616                istatus = 0
1617                exit
1618             end if
1619    
1620          end if
1621       end do
1623    end subroutine get_missing_fill_value
1626    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1627    ! Name: get_halt_on_missing
1628    !
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)
1634       implicit none
1635   
1636       ! Arguments
1637       integer, intent(out) :: istatus
1638       logical, intent(out) :: halt
1639       character (len=128), intent(in) :: fieldnm
1640   
1641       ! Local variables
1642       integer :: idx
1643   
1644       istatus = 0
1645       halt = .false.
1646   
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 
1650    
1651             if (is_halt_missing(idx)) halt = .true.
1652    
1653          end if
1654       end do
1656    end subroutine get_halt_on_missing
1659    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1660    ! Name: get_masked_value
1661    !
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)
1668       implicit none
1669   
1670       ! Arguments
1671       integer, intent(in) :: ilevel
1672       integer, intent(out) :: istatus
1673       real, intent(out) :: masked
1674       character (len=128), intent(in) :: fieldnm
1675   
1676       ! Local variables
1677       integer :: idx
1678   
1679       istatus = 0
1680       masked = -1.
1681   
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
1686    
1687             if (is_masked(idx)) then
1688                masked = source_masked(idx)
1689                exit
1690             end if
1691    
1692          end if
1693       end do
1695    end subroutine get_masked_value
1698    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1699    ! Name: get_max_levels
1700    !
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)
1707       implicit none
1708   
1709       ! Arguments
1710       integer, intent(out) :: min_level, max_level, istatus
1711       character (len=128), intent(in) :: fieldnm
1712   
1713       ! Local variables
1714       integer :: idx
1715       logical :: have_found_field
1716   
1717       have_found_field = .false.
1718       istatus = 1
1719   
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
1723    
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))
1727             end if
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.
1731                   istatus = 0
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.
1736                   istatus = 0
1737                   min_level = 1
1738                   max_level = source_tile_z(idx)
1739                end if
1740      
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)
1745                   end if
1746                end if
1747             else
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)
1754                end if
1755      
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)
1760                   end if
1761                end if
1762             end if
1763    
1764          end if
1765       end do
1767    end subroutine get_max_levels
1770    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1771    ! Name: get_source_levels
1772    !
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)
1778      
1779       implicit none
1780   
1781       ! Arguments
1782       integer, intent(in) :: ilevel
1783       integer, intent(out) :: min_level, max_level, istatus
1784       character (len=128), intent(in) :: fieldnm
1785   
1786       ! Local variables
1787       integer :: idx
1788   
1789       istatus = 1
1790   
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
1795     
1796                if (is_tile_z_start(idx) .and. is_tile_z_end(idx)) then
1797                   istatus = 0
1798                   min_level = source_tile_z_start(idx)
1799                   max_level = source_tile_z_end(idx)
1800                else if (is_tile_z(idx)) then
1801                   istatus = 0
1802                   min_level = 1
1803                   max_level = source_tile_z(idx)
1804                end if
1805      
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)
1810                   end if 
1811                end if
1812     
1813             end if
1814          end if
1815       end do
1817    end subroutine get_source_levels
1819    
1820    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1821    ! Name: get_max_categories
1822    !
1823    ! Purpose: Returns the minimum category and the maximum category for the field
1824    !   given by fieldnm.
1825    !   Min/max category will generally be specified for categorical fields, 
1826    !   whereas the number of levels will generally be specified for continuous 
1827    !   fields. 
1828    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1829    subroutine get_max_categories(fieldnm, min_category, max_category, istatus)
1831       implicit none
1832   
1833       ! Arguments
1834       integer, intent(out) :: min_category, max_category, istatus
1835       character (len=128), intent(in) :: fieldnm
1836   
1837       ! Local variables
1838       integer :: idx
1839       logical :: have_found_field
1840   
1841       have_found_field = .false.
1842       istatus = 1
1843   
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
1847    
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)
1853             end if
1854             if (.not. have_found_field) then
1855                if (is_category_min(idx) .and. is_category_max(idx)) then
1856                   have_found_field = .true.
1857                   istatus = 0
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)
1863                end if
1864             else
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)
1871                end if
1872             end if
1873    
1874          end if
1875       end do
1876   
1877    end subroutine get_max_categories
1880    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1881    ! Name: get_source_categories
1882    !
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)
1888      
1889       implicit none
1890   
1891       ! Arguments
1892       integer, intent(in) :: ilevel
1893       integer, intent(out) :: min_category, max_category, istatus
1894       character (len=128), intent(in) :: fieldnm
1895   
1896       ! Local variables
1897       integer :: idx
1898   
1899       istatus = 1
1900   
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
1905     
1906                if (is_category_min(idx) .and. is_category_max(idx)) then
1907                   istatus = 0
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)
1913                end if 
1914     
1915             end if
1916          end if
1917       end do
1919    end subroutine get_source_categories
1922    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1923    ! Name: get_domcategory_name
1924    !
1925    ! Purpose:
1926    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1927    subroutine get_domcategory_name(fieldnm, domcat_name, ldominant_only, istatus)
1929       implicit none
1930   
1931       ! Arguments
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
1936   
1937       ! Local variables
1938       integer :: idx
1939   
1940       istatus = 1
1941       ldominant_only = .false.
1942   
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
1946    
1947             if (is_dominant_category(idx)) then
1948                domcat_name = source_dominant_category(idx) 
1949                istatus = 0
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.
1955                end if
1956                exit
1957     
1958             else if (is_dominant_only(idx)) then
1959                domcat_name = source_dominant_only(idx) 
1960                ldominant_only = .true.
1961                istatus = 0
1962                exit
1963             end if
1964    
1965          end if
1966       end do
1968    end subroutine get_domcategory_name
1971    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1972    ! Name: get_dfdx_name
1973    !
1974    ! Purpose:
1975    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1976    subroutine get_dfdx_name(fieldnm, dfdx_name, istatus)
1978       implicit none
1979   
1980       ! Arguments
1981       integer, intent(out) :: istatus
1982       character (len=128), intent(in) :: fieldnm
1983       character (len=128), intent(out) :: dfdx_name
1984   
1985       ! Local variables
1986       integer :: idx
1987   
1988       istatus = 1
1989   
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
1993    
1994             if (is_dfdx(idx)) then
1995                dfdx_name = source_dfdx(idx) 
1996                istatus = 0
1997                exit
1998             end if
1999    
2000          end if
2001       end do
2003    end subroutine get_dfdx_name
2006    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2007    ! Name: get_dfdy_name
2008    !
2009    ! Purpose:
2010    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2011    subroutine get_dfdy_name(fieldnm, dfdy_name, istatus)
2013       implicit none
2014   
2015       ! Arguments
2016       integer, intent(out) :: istatus
2017       character (len=128), intent(in) :: fieldnm
2018       character (len=128), intent(out) :: dfdy_name
2019   
2020       ! Local variables
2021       integer :: idx
2022   
2023       istatus = 1
2024   
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
2028    
2029             if (is_dfdy(idx)) then
2030                dfdy_name = source_dfdy(idx) 
2031                istatus = 0
2032                exit
2033             end if
2034    
2035          end if
2036       end do
2038    end subroutine get_dfdy_name
2041    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2042    ! Name: get_z_dim_name
2043    !
2044    ! Purpose:
2045    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2046    subroutine get_z_dim_name(fieldnm, z_dim, istatus)
2047    
2048       implicit none
2049   
2050       ! Arguments
2051       integer, intent(out) :: istatus
2052       character (len=128), intent(in) :: fieldnm
2053       character (len=128), intent(out) :: z_dim
2054   
2055       ! Local variables
2056       integer :: idx
2057   
2058       istatus = 1
2059   
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)
2065                istatus = 0
2066                exit
2067             end if
2068          end if
2069       end do
2070      
2071    end subroutine get_z_dim_name
2074    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2075    ! Name: get_field_scale_factor
2076    !
2077    ! Purpose:
2078    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2079    subroutine get_field_scale_factor(fieldnm, ilevel, scale_factor, istatus)
2080    
2081       implicit none
2082   
2083       ! Arguments
2084       integer, intent(in) :: ilevel
2085       integer, intent(out) :: istatus
2086       real, intent(out) :: scale_factor
2087       character (len=128), intent(in) :: fieldnm
2088   
2089       ! Local variables
2090       integer :: idx
2091   
2092       istatus = 1
2093   
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
2098    
2099             if (is_scale_factor(idx)) then
2100                scale_factor = source_scale_factor(idx) 
2101                istatus = 0
2102             end if
2103    
2104          end if
2105       end do
2107    end subroutine get_field_scale_factor
2110    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2111    ! Name: get_output_stagger
2112    !
2113    ! Pupose:
2114    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2115    subroutine get_output_stagger(fieldnm, istagger, istatus)
2117       use gridinfo_module
2118     
2119       implicit none
2120   
2121       ! Arguments
2122       integer, intent(out) :: istatus, istagger
2123       character (len=128), intent(in) :: fieldnm
2124   
2125       ! Local variables
2126       integer :: idx
2127   
2128       istatus = 1
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
2132    
2133             if (is_output_stagger(idx)) then
2134                istatus = 0
2135                istagger = source_output_stagger(idx)
2136                exit
2137             else
2138                if (gridtype == 'C') then
2139                   istatus = 0
2140                   istagger = M
2141                   exit
2142                else if (gridtype == 'E') then
2143                   istatus = 0
2144                   istagger = HH
2145                   exit
2146                end if
2147             end if
2148    
2149          end if
2150       end do
2152    end subroutine get_output_stagger
2155    subroutine get_subgrid_dim_name(nest_num,field_name,dimnames, &
2156                                    sub_x,sub_y,istatus)
2157    use gridinfo_module
2158    use module_debug
2159    implicit none
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
2164    integer :: idx,nlen
2166    sub_x=1
2167    sub_y=1
2169    istatus = 0
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
2174              istatus = 0
2175              if(is_output_stagger(idx))then
2176                 call mprintf(.true.,ERROR,'Cannot use subgrids on variables with staggered grids')
2177              endif
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)
2182          endif
2183       endif
2184    enddo
2185    end subroutine get_subgrid_dim_name
2188    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2189    ! Name: get_interp_option
2190    !
2191    ! Pupose:
2192    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2193    subroutine get_interp_option(fieldnm, ilevel, interp_opt, istatus)
2194    
2195       implicit none
2196   
2197       ! Arguments
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
2202   
2203       ! Local variables
2204       integer :: idx
2205   
2206       istatus = 1
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
2211     
2212                interp_opt = source_interp_string(idx)
2213                istatus = 0
2214                exit
2215     
2216             end if
2217          end if
2218       end do
2219   
2220    end subroutine get_interp_option
2223    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2224    ! Name: get_gcel_threshold
2225    !
2226    ! Pupose:
2227    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2228    subroutine get_gcell_threshold(interp_opt, threshold, istatus)
2229    
2230       implicit none
2231   
2232       ! Arguments
2233       integer, intent(out) :: istatus
2234       real, intent(out) :: threshold
2235       character (len=128), intent(in) :: interp_opt
2237       ! Local variables
2238       integer :: i, p1, p2
2240       istatus = 1
2241       threshold = 1.0
2242     
2243       i = index(interp_opt,'average_gcell')
2244       if (i /= 0) then
2245          ! Move the "average_gcell" option to the beginning
2246 !         if (i /= 1) then
2247 !            p1 =  
2248 !         end if
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
2255          else
2256             call mprintf(.true., WARN, 'Problem with specified threshold '// &
2257                          'for average_gcell interp option. Setting threshold to 0.0.')
2258             threshold = 0.0
2259          end if
2260       end if
2261       istatus = 0
2262     
2263       return
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"')
2268   
2269    end subroutine get_gcell_threshold
2272    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2273    ! Name: get_smooth_option
2274    !
2275    ! Pupose:
2276    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2277    subroutine get_smooth_option(fieldnm, smooth_opt, smooth_passes, istatus)
2279       implicit none
2280   
2281       ! Arguments
2282       integer, intent(out) :: istatus, smooth_opt, smooth_passes
2283       character (len=128), intent(in) :: fieldnm
2284   
2285       ! Local variables
2286       integer :: idx
2287   
2288       istatus = 1
2289   
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
2293    
2294             if (is_smooth_option(idx)) then
2295                istatus = 0
2296                smooth_opt = source_smooth_option(idx)
2297      
2298                if (is_smooth_passes(idx)) then
2299                   smooth_passes = source_smooth_passes(idx)
2300                else
2301                   smooth_passes = 1
2302                end if
2303      
2304                exit
2305             end if
2306    
2307          end if
2308       end do
2310    end subroutine get_smooth_option
2313    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2314    ! Name: iget_fieldtype
2315    !
2316    ! Purpose:
2317    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2318    function iget_fieldtype(fieldnm, istatus)
2320       implicit none
2321     
2322       ! Arguments
2323       integer, intent(out) :: istatus
2324       character (len=128), intent(in) :: fieldnm
2325   
2326       ! Local variables
2327       integer :: idx
2328   
2329       ! Return value
2330       integer :: iget_fieldtype
2331      
2332       istatus = 1
2333   
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
2337    
2338             if (is_dest_fieldtype(idx)) then
2339                iget_fieldtype = source_dest_fieldtype(idx) 
2340                istatus = 0
2341                exit
2342             end if
2343    
2344          end if
2345       end do
2347    end function iget_fieldtype
2350    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2351    ! Name: iget_source_fieldtype
2352    !
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)
2362       implicit none
2363     
2364       ! Arguments
2365       integer, intent(in) :: ilevel
2366       integer, intent(out) :: istatus
2367       character (len=128), intent(in) :: fieldnm
2368   
2369       ! Return value
2370       integer :: iget_source_fieldtype
2371   
2372       ! Local variables
2373       integer :: idx
2374   
2375       ! Find information about the source tiles for the specified fieldname
2376       istatus = 1
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
2381                istatus = 0
2382                iget_source_fieldtype = source_fieldtype(idx)
2383                exit
2384             end if
2385          end if
2386       end do
2388    end function iget_source_fieldtype
2391    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2392    ! Name: get_data_tile
2393    !
2394    ! Purpose:
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, &
2399                             istatus)
2401       implicit none
2402   
2403       ! Arguments
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, &
2409                               npts_bdr
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
2414   
2415       ! Local variables
2416       integer :: j, k
2417       integer :: local_wordsize, local_endian, sign_convention, irow_order, strlen
2418       integer :: xdim,ydim,zdim
2419       real :: scalefac
2420       real, allocatable, dimension(:) :: temprow
2421   
2422       call get_tile_fname(file_name, xlat, xlon, ilevel, field_name, istatus)
2424       if (index(file_name, 'OUTSIDE') /= 0) then
2425          istatus = 1
2426          return
2427       else if (hash_search(bad_files, file_name)) then
2428          istatus = 1
2429          return
2430       end if
2431   
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)
2435   
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))
2442   
2443       call get_row_order(field_name, ilevel, irow_order, istatus)
2444       if (istatus /= 0) irow_order = BOTTOM_TOP
2445   
2446       call s_len(file_name,strlen)
2448       scalefac = 1.0
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))
2455          do k=1,zdim
2456             do j=1,ydim
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)
2461             end do
2462          end do
2463          deallocate(temprow)
2464       end if
2465   
2466       if (istatus /= 0) then
2467          start_x_dim = INVALID
2468          start_y_dim = INVALID
2469          end_x_dim   = INVALID
2470          end_y_dim   = INVALID
2472          call hash_insert(bad_files, file_name)
2473       end if
2475    end subroutine get_data_tile
2478    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2479    ! Name: get_row_order
2480    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2481    subroutine get_row_order(fieldnm, ilevel, irow_order, istatus)
2483       implicit none
2485       ! Arguments
2486       integer, intent(in) :: ilevel
2487       character (len=128), intent(in) :: fieldnm
2488       integer, intent(out) :: irow_order, istatus
2490       ! Local variables
2491       integer :: idx
2493       istatus = 1
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)
2500                   istatus = 0
2501                   exit
2502                end if
2503             end if
2504          end if
2505       end do
2507    end subroutine get_row_order
2508   
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)
2517       use llxy_module
2518   
2519       implicit none
2520   
2521       ! Arguments
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, &
2526                               npts_bdr, &
2527                               bytes_per_datum, endianness, &
2528                               sign_convention, istatus
2529       real, intent(in) :: xlat, xlon
2530       character (len=128), intent(in) :: fieldnm
2531   
2532       ! Local variables
2533       integer :: idx, current_domain
2534       real :: rx, ry
2535   
2536       istatus = 1
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
2541                istatus = 0
2542                exit
2543             end if
2544          end if
2545       end do
2546   
2547       if (istatus /= 0) then
2548          start_x_dim = 1
2549          start_y_dim = 1
2550          start_z_dim = 1
2551          end_x_dim = 1
2552          end_y_dim = 1
2553          end_z_dim = 1
2554          bytes_per_datum = 0
2555          return
2556       end if
2557   
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)
2562   
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
2565   
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
2568   
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
2573          start_z_dim = 1
2574          end_z_dim = source_tile_z(idx)
2575       end if
2576   
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)
2581          end if 
2582       end if
2583   
2584       if (is_tile_bdr(idx)) then
2585          npts_bdr = source_tile_bdr(idx)
2586       else
2587          npts_bdr = 0
2588       end if
2589   
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
2594   
2595       if (is_wordsize(idx)) then
2596          bytes_per_datum = source_wordsize(idx)
2597       else
2598          bytes_per_datum = 1
2599          call mprintf(.true.,ERROR,'In GEOGRID.TBL, no wordsize specified for data in entry %i.',i1=idx)
2600       end if
2602       if (is_endian(idx)) then
2603          endianness = source_endian(idx)
2604       else
2605          endianness = BIG_ENDIAN
2606       end if
2607   
2608       if (is_signed(idx)) then
2609          sign_convention = 1
2610       else
2611          sign_convention = 0
2612       end if
2614    end subroutine get_tile_dimensions
2617    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
2618    ! Name: get_tile_fname
2619    !
2620    ! Purpose: 
2621    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
2622    subroutine get_tile_fname(test_fname, xlat, xlon, ilevel, fieldname, istatus)
2624       use llxy_module
2625       use gridinfo_module
2626   
2627       implicit none
2628   
2629       ! Arguments
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
2635   
2636       ! Local variables
2637       integer :: current_domain, idx
2638       real :: rx, ry
2639   
2640       istatus = 1
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
2647                istatus = 0
2648                exit
2649             end if
2650          end if
2651       end do
2652   
2653       if (istatus /= 0) return
2654   
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
2668       end if
2670    end subroutine get_tile_fname
2673    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
2674    ! Name: get_source_resolution
2675    !
2676    ! Purpose:
2677    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
2678    subroutine get_source_resolution(fieldnm, ilevel, src_dx, src_dy, istatus)
2680       implicit none
2682       ! Arguments
2683       integer, intent(in) :: ilevel
2684       integer, intent(out) :: istatus
2685       real, intent(out) :: src_dx, src_dy
2686       character (len=*), intent(in) :: fieldnm
2687   
2688       ! Local variables
2689       integer :: idx
2690   
2691       istatus = 1
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.
2702                   end if
2703                   istatus = 0
2704                   exit
2705                end if 
2706             end if
2707          end if
2708       end do
2710    end subroutine get_source_resolution
2713    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
2714    ! Name: get_data_projection
2715    !
2716    ! Purpose: To acquire the parameters necessary in defining the grid on which
2717    !   the user-specified data for field 'fieldnm' are given.
2718    ! 
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)
2727       implicit none
2728   
2729       ! Arguments
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
2735   
2736       ! Local variables
2737       integer :: idx
2738   
2739       istatus = 1
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
2744                istatus = 0
2745                if (is_proj(idx)) then
2746                   iproj = source_proj(idx) 
2747                else
2748                   iproj = 1 
2749                   call mprintf(.true., ERROR, &
2750                                'In GEOGRID.TBL, no specification for projection in entry %i.',i1=idx)
2751                end if
2752                if (is_known_x(idx)) then
2753                   known_x = source_known_x(idx) 
2754                else
2755                   known_x = 1. 
2756                   call mprintf(.true., ERROR, &
2757                                'In GEOGRID.TBL, no specification for known_x in entry %i.',i1=idx)
2758                end if
2759                if (is_known_y(idx)) then
2760                   known_y =  source_known_y(idx)
2761                else
2762                   known_y = 1. 
2763                   call mprintf(.true., ERROR, &
2764                                'In GEOGRID.TBL, no specification for known_y in entry %i.',i1=idx)
2765                end if
2766                if (is_known_lat(idx)) then
2767                   known_lat = source_known_lat(idx)
2768                else
2769                   known_lat = 1.
2770                   call mprintf(.true., ERROR, &
2771                                'In GEOGRID.TBL, no specification for known_lat in entry %i.',i1=idx)
2772                end if
2773                if (is_known_lon(idx)) then
2774                   known_lon = source_known_lon(idx)
2775                else
2776                   known_lon = 1.
2777                   call mprintf(.true., ERROR, &
2778                                'In GEOGRID.TBL, no specification for known_lon in entry %i.',i1=idx)
2779                end if
2780                if (is_truelat1(idx)) then
2781                   truelat1 = source_truelat1(idx)
2782                else if (is_proj(idx) .and. source_proj(idx) /= PROJ_LATLON) then
2783                   truelat1 = 1.
2784                   call mprintf(.true., WARN, &
2785                                'In GEOGRID.TBL, no specification for truelat1 in entry %i.',i1=idx)
2786                end if
2787                if (is_truelat2(idx)) then
2788                   truelat2 = source_truelat2(idx)
2789                else if (is_proj(idx) .and. source_proj(idx) /= PROJ_LATLON) then
2790                   truelat2 = 1.
2791                   call mprintf(.true., WARN, &
2792                                'In GEOGRID.TBL, no specification for truelat2 in entry %i.',i1=idx)
2793                end if
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
2797                   stand_lon = 1.
2798                   call mprintf(.true., WARN, &
2799                                'In GEOGRID.TBL, no specification for stdlon in entry %i.',i1=idx)
2800                end if
2801                if (is_dx(idx)) then
2802                   dxkm = source_dx(idx)
2803                else
2804                   dxkm = 1. 
2805                   call mprintf(.true., ERROR, &
2806                                'In GEOGRID.TBL, no specification for dx in entry %i.',i1=idx)
2807                end if
2808                if (is_dy(idx)) then
2809                   dykm = source_dy(idx)
2810                else
2811                   dykm = 1. 
2812                   call mprintf(.true., ERROR, &
2813                                'In GEOGRID.TBL, no specification for dy in entry %i.',i1=idx)
2814                end if
2815                exit
2816             end if
2817          end if
2818       end do
2820    end subroutine get_data_projection
2823    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2824    ! Name: check_data_specification
2825    !
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( )
2831   
2832       implicit none
2833   
2834       ! Return value
2835       logical :: check_data_specification
2836   
2837       ! Local variables
2838       integer :: i, j, istatus
2839       integer, pointer, dimension(:) :: priorities
2840       real :: rmissing
2841       logical :: begin_priority, halt
2842       character (len=128) :: cur_name
2843   
2844       check_data_specification = .false.
2845   
2846       ! Check that each specification has a name, priority level, and path
2847       do i=1,num_entries
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)
2851          end if
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)
2855          end if
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)
2860          end if
2861       end do
2862   
2863       ! The fill_missing and halt_on_missing options are mutually exclusive
2864       do i=1,num_entries
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)))
2871          end if
2872       end do
2873   
2874       ! Check that the field from which landmask is calculated is not output on a staggering
2875       do i=1,num_entries
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)
2881                end if
2882             end if
2883          end if
2884       end do
2885   
2886       ! Also check that any field that is to be masked by the landmask is not output on a staggering
2887       do i=1,num_entries
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)
2892             end if
2893          end if
2894       end do
2895   
2896       allocate(priorities(num_entries))
2897   
2898       ! Now check that priorities for each source are unique and in the interval [1,n], n <= num_entries 
2899       do i=1,num_entries
2900          priorities = 0
2901          cur_name = source_fieldname(i)
2902          do j=1,num_entries
2903             if (source_fieldname(j) == cur_name) then
2904     
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))
2908      
2909                else 
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))
2913                 
2914                   else
2915                      priorities(source_priority(j)) = 1
2916                   end if
2917                end if
2918     
2919             end if
2920          end do
2921    
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)
2930             end if
2931          end do
2932       end do
2933   
2934       deallocate(priorities)
2935   
2936       ! Units must match for all priority levels of a field
2937       do i=1,num_entries
2938          if (source_priority(i) == 1) then
2939             do j=1,num_entries
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)))
2945                end if
2946             end do
2947          end if
2948       end do
2949   
2950       ! Make sure that user has not asked to calculate landmask from a continuous field
2951       do i=1,num_entries
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)
2957                end if
2958             end if
2959          end if
2960       end do
2961   
2962       ! If either min_category or max_category is specified, then both must be specified
2963       do i=1,num_entries
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)
2973             end if
2974          end if
2975       end do
2976   
2977       ! For continuous data, (category_max - category_min + 1) should equal tile_z
2978       do i=1,num_entries
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)
2985                   end if
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)
2993                   end if
2994                end if
2995             end if
2996          end if
2997       end do
2998   
2999       ! Make sure that user has not named a dominant category or computed slope field 
3000       !    the same as a fractional field
3001       do i=1,num_entries
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)
3005          end if
3006    
3007          do j=1,num_entries
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.', &
3013                                   i1=i, i2=j)
3014                   end if
3015                end if
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.', &
3020                                   i1=i, i2=j)
3021                   end if
3022                end if
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)
3028                   end if
3029                end if
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)
3035                   end if
3036                end if
3037             else
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)
3043                   end if
3044                end if
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)
3050                   end if
3051                end if
3052             end if
3053             if (i /= j) then
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)
3060                      end if
3061                   end if
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)
3067                      end if
3068                   end if
3069                end if
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)
3076                      end if
3077                   end if
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)
3083                      end if
3084                   end if
3085                end if
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)
3092                      end if
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)
3103                   end if
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)
3110                      end if
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)
3121                   end if
3122                end if
3123             end if
3124          end do
3125       end do
3126   
3127       check_data_specification = .true.
3128    end function check_data_specification
3131    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
3132    ! Name: s_len
3133    !
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)
3141       implicit none
3142   
3143       ! Arguments
3144       character (len=*), intent(in) :: string
3145       integer, intent(out) :: s_length
3146   
3147       ! Local variables
3148       integer :: i, len_str, aval
3149       logical :: space
3150   
3151       space = .false.
3152       i = 1
3153       len_str = len(string)
3154       s_length = len_str     
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
3158             s_length = i - 1
3159             space = .true.
3160          endif
3161          i = i + 1
3162       enddo
3164    end subroutine s_len
3166 end module source_data_module