Update version number to 4.0.3
[WPS-merge.git] / geogrid / src / process_tile_module.F
blob0007436e2eaf90ea635ce73bb22ef2cb73398aa7
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 ! Module: process_tile
4 ! Description:
5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
6 module process_tile_module
8    use module_debug
11    contains
14    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15    ! Name: process_tile
16    !
17    ! Purpose: To process a tile, whose lower-left corner is at 
18    !       (tile_i_min, tile_j_min) and whose upper-right corner is at
19    !       (tile_i_max, tile_j_max), of the model grid given by which_domain 
20    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
21    subroutine process_tile(which_domain, grid_type, dynopt,        &
22                            dummy_start_dom_i, dummy_end_dom_i,     &
23                            dummy_start_dom_j, dummy_end_dom_j,     &
24                            dummy_start_patch_i, dummy_end_patch_i, &
25                            dummy_start_patch_j, dummy_end_patch_j, &
26                            extra_col, extra_row)
27    
28       use bitarray_module
29       use hash_module
30       use llxy_module
31       use misc_definitions_module
32       use output_module
33       use smooth_module
34       use source_data_module
35     
36       implicit none
38       ! Arguments
39       integer, intent(in) :: which_domain, dynopt, &
40                              dummy_start_dom_i, dummy_end_dom_i, dummy_start_dom_j, dummy_end_dom_j, &
41                              dummy_start_patch_i, dummy_end_patch_i, dummy_start_patch_j, dummy_end_patch_j
42       logical, intent(in) :: extra_col, extra_row
43       character (len=1), intent(in) :: grid_type
44     
45       ! Local variables
46       integer :: i, j, k, kk, istatus, ifieldstatus, idomcatstatus, field_count
47       integer :: min_category, max_category, min_level, max_level, &
48                  smth_opt, smth_passes, num_landmask_categories
49       integer :: start_dom_i, end_dom_i, start_dom_j, end_dom_j, end_dom_stag_i, end_dom_stag_j
50       integer :: start_patch_i, end_patch_i, start_patch_j, end_patch_j, end_patch_stag_i, end_patch_stag_j
51       integer :: start_mem_i, end_mem_i, start_mem_j, end_mem_j, end_mem_stag_i, end_mem_stag_j
52       integer :: sm1, em1, sm2, em2
53       integer :: istagger
54       integer, dimension(MAX_LANDMASK_CATEGORIES) :: landmask_value
55       real :: sum, dominant, msg_fill_val, topo_flag_val, mass_flag, land_total, water_total
56       real, dimension(16) :: corner_lats, corner_lons
57       real, pointer, dimension(:,:) :: xlat_array,   xlon_array,   &
58                                        xlat_array_u, xlon_array_u, &
59                                        xlat_array_v, xlon_array_v, &
60                                        xlat_array_corner, xlon_array_corner, &
61                                        clat_array,   clon_array,   &
62                                        xlat_array_subgrid, xlon_array_subgrid, &
63                                        f_array, e_array, &
64                                        mapfac_array_m_x, mapfac_array_u_x, mapfac_array_v_x, &
65                                        mapfac_array_m_y, mapfac_array_u_y, mapfac_array_v_y, &
66                                        mapfac_array_x_subgrid, mapfac_array_y_subgrid,       &
67                                        sina_array, cosa_array
68       real, pointer, dimension(:,:) :: xlat_ptr, xlon_ptr, mapfac_ptr_x, mapfac_ptr_y, landmask, dominant_field
69       real, pointer, dimension(:,:,:) :: field, slp_field
70       logical :: is_water_mask, only_save_dominant, halt_on_missing
71       character (len=19) :: datestr
72       character (len=128) :: fieldname, gradname, domname, landmask_name
73       character (len=256) :: temp_string
74       type (bitarray) :: processed_domain 
75       type (hashtable) :: processed_fieldnames
76       character (len=128), dimension(2) :: dimnames
77       integer :: sub_x, sub_y
78       integer :: opt_status
80       ! Probably not all of these nullify statements are needed...
81       nullify(xlat_array)
82       nullify(xlon_array)
83       nullify(xlat_array_u)
84       nullify(xlon_array_u)
85       nullify(xlat_array_v)
86       nullify(xlon_array_v)
87       nullify(xlat_array_corner)
88       nullify(xlon_array_corner)
89       nullify(clat_array)
90       nullify(clon_array)
91       nullify(xlat_array_subgrid)
92       nullify(xlon_array_subgrid)
93       nullify(f_array)
94       nullify(e_array)
95       nullify(mapfac_array_m_x)
96       nullify(mapfac_array_u_x)
97       nullify(mapfac_array_v_x)
98       nullify(mapfac_array_m_y)
99       nullify(mapfac_array_u_y)
100       nullify(mapfac_array_v_y)
101       nullify(mapfac_array_x_subgrid)
102       nullify(mapfac_array_y_subgrid)
103       nullify(sina_array)
104       nullify(cosa_array)
105       nullify(xlat_ptr)
106       nullify(xlon_ptr)
107       nullify(mapfac_ptr_x)
108       nullify(mapfac_ptr_y)
109       nullify(landmask)
110       nullify(dominant_field)
111       nullify(field)
112       nullify(slp_field)
113    
114       datestr = '0000-00-00_00:00:00'
115       field_count = 0
116       mass_flag=1.0
117     
118       ! The following pertains primarily to the C grid
119       ! Determine whether only (n-1)th rows/columns should be computed for variables
120       !   on staggered grid. In a distributed memory situation, not every tile should
121       !   have only (n-1)th rows/columns computed, or we end up with (n-k) 
122       !   rows/columns when there are k patches in the y/x direction
123       if (extra_col) then
124          start_patch_i    = dummy_start_patch_i    ! The seemingly pointless renaming of start
125          end_patch_i      = dummy_end_patch_i - 1  !   naming convention with modified end_patch variables, 
126          end_patch_stag_i = dummy_end_patch_i      !   variables is so that we can maintain consistent
127                                                    !   which are marked as intent(in)
128          start_mem_i    = start_patch_i    - HALO_WIDTH
129          end_mem_i      = end_patch_i      + HALO_WIDTH
130          end_mem_stag_i = end_patch_stag_i + HALO_WIDTH
131       else                                     
132          start_patch_i    = dummy_start_patch_i
133          end_patch_i      = dummy_end_patch_i
134          end_patch_stag_i = dummy_end_patch_i
136          start_mem_i    = start_patch_i    - HALO_WIDTH
137          end_mem_i      = end_patch_i      + HALO_WIDTH
138          end_mem_stag_i = end_patch_stag_i + HALO_WIDTH
139       end if
140     
141       if (extra_row) then
142          start_patch_j    = dummy_start_patch_j
143          end_patch_j      = dummy_end_patch_j - 1
144          end_patch_stag_j = dummy_end_patch_j
146          start_mem_j    = start_patch_j    - HALO_WIDTH
147          end_mem_j      = end_patch_j      + HALO_WIDTH
148          end_mem_stag_j = end_patch_stag_j + HALO_WIDTH
149       else
150          start_patch_j    = dummy_start_patch_j
151          end_patch_j      = dummy_end_patch_j
152          end_patch_stag_j = dummy_end_patch_j
154          start_mem_j    = start_patch_j    - HALO_WIDTH
155          end_mem_j      = end_patch_j      + HALO_WIDTH
156          end_mem_stag_j = end_patch_stag_j + HALO_WIDTH
157       end if
159       start_dom_i = dummy_start_dom_i
160       if (grid_type == 'C') then
161          end_dom_i      = dummy_end_dom_i - 1
162          end_dom_stag_i = dummy_end_dom_i
163       else if (grid_type == 'E') then
164          end_dom_i      = dummy_end_dom_i
165          end_dom_stag_i = dummy_end_dom_i
166       end if
168       start_dom_j = dummy_start_dom_j
169       if (grid_type == 'C') then
170          end_dom_j      = dummy_end_dom_j - 1
171          end_dom_stag_j = dummy_end_dom_j
172       else if (grid_type == 'E') then
173          end_dom_j      = dummy_end_dom_j
174          end_dom_stag_j = dummy_end_dom_j
175       end if
176     
177       ! Allocate arrays to hold all lat/lon fields; these will persist for the duration of
178       !   the process_tile routine
179       ! For C grid, we have M, U, and V points
180       ! For E grid, we have only M and V points
181       allocate(xlat_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
182       allocate(xlon_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
183       allocate(xlat_array_v(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
184       allocate(xlon_array_v(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
185       if (grid_type == 'C') then
186          allocate(xlat_array_u(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
187          allocate(xlon_array_u(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
188          allocate(clat_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
189          allocate(clon_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
190          allocate(xlat_array_corner(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_stag_j))
191          allocate(xlon_array_corner(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_stag_j))
192       end if
193       nullify(xlat_array_subgrid)
194       nullify(xlon_array_subgrid)
195       nullify(mapfac_array_x_subgrid)
196       nullify(mapfac_array_y_subgrid)
197     
198       ! Initialize hash table to track which fields have been processed
199       call hash_init(processed_fieldnames)
200     
201       !
202       ! Calculate lat/lon for every point in the tile (XLAT and XLON)
203       ! The xlat_array and xlon_array arrays will be used in processing other fields
204       !
205       call mprintf(.true.,STDOUT,'  Processing XLAT and XLONG')
206     
207       if (grid_type == 'C') then
208          call get_lat_lon_fields(xlat_array, xlon_array, start_mem_i, &
209                                start_mem_j, end_mem_i, end_mem_j, M)
210          call get_lat_lon_fields(xlat_array_v, xlon_array_v, start_mem_i, &
211                                start_mem_j, end_mem_i, end_mem_stag_j, V)
212          call get_lat_lon_fields(xlat_array_u, xlon_array_u, start_mem_i, &
213                                start_mem_j, end_mem_stag_i, end_mem_j, U)
214          call get_lat_lon_fields(xlat_array_corner, xlon_array_corner, start_mem_i, &
215                                start_mem_j, end_mem_stag_i, end_mem_stag_j, CORNER)
216          call get_lat_lon_fields(clat_array, clon_array, start_mem_i, &
217                                start_mem_j, end_mem_i, end_mem_j, M, comp_ll=.true.)
219          corner_lats(1) = xlat_array(start_patch_i,start_patch_j)
220          corner_lats(2) = xlat_array(start_patch_i,end_patch_j)
221          corner_lats(3) = xlat_array(end_patch_i,end_patch_j)
222          corner_lats(4) = xlat_array(end_patch_i,start_patch_j)
223      
224          corner_lats(5) = xlat_array_u(start_patch_i,start_patch_j)
225          corner_lats(6) = xlat_array_u(start_patch_i,end_patch_j)
226          corner_lats(7) = xlat_array_u(end_patch_stag_i,end_patch_j)
227          corner_lats(8) = xlat_array_u(end_patch_stag_i,start_patch_j)
228      
229          corner_lats(9)  = xlat_array_v(start_patch_i,start_patch_j)
230          corner_lats(10) = xlat_array_v(start_patch_i,end_patch_stag_j)
231          corner_lats(11) = xlat_array_v(end_patch_i,end_patch_stag_j)
232          corner_lats(12) = xlat_array_v(end_patch_i,start_patch_j)
233      
234          call xytoll(real(start_patch_i)-0.5, real(start_patch_j)-0.5, corner_lats(13), corner_lons(13), M)
235          call xytoll(real(start_patch_i)-0.5, real(end_patch_j)+0.5, corner_lats(14), corner_lons(14), M)
236          call xytoll(real(end_patch_i)+0.5, real(end_patch_j)+0.5, corner_lats(15), corner_lons(15), M)
237          call xytoll(real(end_patch_i)+0.5, real(start_patch_j)-0.5, corner_lats(16), corner_lons(16), M)
239          corner_lons(1) = xlon_array(start_patch_i,start_patch_j)
240          corner_lons(2) = xlon_array(start_patch_i,end_patch_j)
241          corner_lons(3) = xlon_array(end_patch_i,end_patch_j)
242          corner_lons(4) = xlon_array(end_patch_i,start_patch_j)
243      
244          corner_lons(5) = xlon_array_u(start_patch_i,start_patch_j)
245          corner_lons(6) = xlon_array_u(start_patch_i,end_patch_j)
246          corner_lons(7) = xlon_array_u(end_patch_stag_i,end_patch_j)
247          corner_lons(8) = xlon_array_u(end_patch_stag_i,start_patch_j)
248      
249          corner_lons(9)  = xlon_array_v(start_patch_i,start_patch_j)
250          corner_lons(10) = xlon_array_v(start_patch_i,end_patch_stag_j)
251          corner_lons(11) = xlon_array_v(end_patch_i,end_patch_stag_j)
252          corner_lons(12) = xlon_array_v(end_patch_i,start_patch_j)
253      
254       else if (grid_type == 'E') then
255          call get_lat_lon_fields(xlat_array, xlon_array, start_mem_i, &
256                                start_mem_j, end_mem_i, end_mem_j, HH)
257          call get_lat_lon_fields(xlat_array_v, xlon_array_v, start_mem_i, &
258                                start_mem_j, end_mem_i, end_mem_stag_j, VV)
259    
260          corner_lats(1) = xlat_array(start_patch_i,start_patch_j)
261          corner_lats(2) = xlat_array(start_patch_i,end_patch_j)
262          corner_lats(3) = xlat_array(end_patch_i,end_patch_j)
263          corner_lats(4) = xlat_array(end_patch_i,start_patch_j)
264      
265          corner_lats(5) = xlat_array_v(start_patch_i,start_patch_j)
266          corner_lats(6) = xlat_array_v(start_patch_i,end_patch_stag_j)
267          corner_lats(7) = xlat_array_v(end_patch_i,end_patch_stag_j)
268          corner_lats(8) = xlat_array_v(end_patch_i,start_patch_j)
269      
270          corner_lats(9)  = 0.0
271          corner_lats(10) = 0.0
272          corner_lats(11) = 0.0
273          corner_lats(12) = 0.0
274      
275          corner_lats(13) = 0.0
276          corner_lats(14) = 0.0
277          corner_lats(15) = 0.0
278          corner_lats(16) = 0.0
279      
280          corner_lons(1) = xlon_array(start_patch_i,start_patch_j)
281          corner_lons(2) = xlon_array(start_patch_i,end_patch_j)
282          corner_lons(3) = xlon_array(end_patch_i,end_patch_j)
283          corner_lons(4) = xlon_array(end_patch_i,start_patch_j)
284      
285          corner_lons(5) = xlon_array_v(start_patch_i,start_patch_j)
286          corner_lons(6) = xlon_array_v(start_patch_i,end_patch_stag_j)
287          corner_lons(7) = xlon_array_v(end_patch_i,end_patch_stag_j)
288          corner_lons(8) = xlon_array_v(end_patch_i,start_patch_j)
289      
290          corner_lons(9)  = 0.0
291          corner_lons(10) = 0.0
292          corner_lons(11) = 0.0
293          corner_lons(12) = 0.0
294      
295          corner_lons(13) = 0.0
296          corner_lons(14) = 0.0
297          corner_lons(15) = 0.0
298          corner_lons(16) = 0.0
299     
300       end if
302       ! Initialize the output module now that we have the corner point lats/lons
303       call output_init(which_domain, 'OUTPUT FROM GEOGRID V4.0.3', '0000-00-00_00:00:00', grid_type, dynopt, &
304                        corner_lats, corner_lons, &
305                        start_dom_i,   end_dom_i,   start_dom_j,   end_dom_j, &
306                        start_patch_i, end_patch_i, start_patch_j, end_patch_j, &
307                        start_mem_i,   end_mem_i,   start_mem_j,   end_mem_j, &
308                        extra_col, extra_row)
309     
310       call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
311                        'XLAT_M', datestr, real_array = xlat_array)
312       call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
313                        'XLONG_M', datestr, real_array = xlon_array)
314       call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, &
315                        'XLAT_V', datestr, real_array = xlat_array_v)
316       call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, &
317                        'XLONG_V', datestr, real_array = xlon_array_v)
318       if (grid_type == 'C') then
319          call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, &
320                         'XLAT_U', datestr, real_array = xlat_array_u)
321          call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, &
322                         'XLONG_U', datestr, real_array = xlon_array_u)
323          call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_stag_j, 1, 1, &
324                         'XLAT_C', datestr, real_array = xlat_array_corner)
325          call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_stag_j, 1, 1, &
326                         'XLONG_C', datestr, real_array = xlon_array_corner)
327          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
328                         'CLAT', datestr, real_array = clat_array)
329          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
330                         'CLONG', datestr, real_array = clon_array)
332          if (associated(clat_array)) deallocate(clat_array)
333          if (associated(clon_array)) deallocate(clon_array)
335       end if
336    
338       !
339       ! Calculate map factor for current domain
340       !
341       if (grid_type == 'C') then
342          call mprintf(.true.,STDOUT,'  Processing MAPFAC')
343     
344          allocate(mapfac_array_m_x(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
345          allocate(mapfac_array_m_y(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
346          call get_map_factor(xlat_array, xlon_array, mapfac_array_m_x, mapfac_array_m_y, start_mem_i, &
347                              start_mem_j, end_mem_i, end_mem_j)
348 ! Global WRF uses map scale factors in X and Y directions, but "regular" WRF uses a single MSF
349 !    on each staggering. In the case of regular WRF, we can assume that MAPFAC_MX = MAPFAC_MY = MAPFAC_M,
350 !    and so we can simply write MAPFAC_MX as the MAPFAC_M field. Ultimately, when global WRF is
351 !    merged into the WRF trunk, we will need only two map scale factor fields for each staggering,
352 !    in the x and y directions, and these will be the same in the case of non-Cassini projections
353          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_M',  &
354                           datestr, real_array = mapfac_array_m_x)
355          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_MX', &
356                           datestr, real_array = mapfac_array_m_x)
357          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_MY', &
358                          datestr, real_array = mapfac_array_m_y)
359     
360          allocate(mapfac_array_v_x(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
361          allocate(mapfac_array_v_y(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
362          call get_map_factor(xlat_array_v, xlon_array_v, mapfac_array_v_x, mapfac_array_v_y, start_mem_i, &
363                              start_mem_j, end_mem_i, end_mem_stag_j)
364          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_V',  &
365                           datestr, real_array = mapfac_array_v_x)
366          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_VX', &
367                           datestr, real_array = mapfac_array_v_x)
368          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_VY', &
369                           datestr, real_array = mapfac_array_v_y)
370     
371          allocate(mapfac_array_u_x(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
372          allocate(mapfac_array_u_y(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
373          call get_map_factor(xlat_array_u, xlon_array_u, mapfac_array_u_x, mapfac_array_u_y, start_mem_i, &
374                              start_mem_j, end_mem_stag_i, end_mem_j)
375          call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_U',  &
376                           datestr, real_array = mapfac_array_u_x)
377          call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_UX', &
378                           datestr, real_array = mapfac_array_u_x)
379          call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_UY', &
380                           datestr, real_array = mapfac_array_u_y)
382       end if
383     
384     
385       !
386       ! Coriolis parameters (E and F)
387       !
388       call mprintf(.true.,STDOUT,'  Processing F and E')
389     
390       allocate(f_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
391       allocate(e_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
392     
393       call get_coriolis_parameters(xlat_array, f_array, e_array, &
394                                    start_mem_i, start_mem_j, end_mem_i, end_mem_j)
395     
396       call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'E', &
397                        datestr, real_array = e_array)
398       call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'F', &
399                        datestr, real_array = f_array)
400     
401       if (associated(f_array)) deallocate(f_array)
402       if (associated(e_array)) deallocate(e_array)
403     
404     
405       !
406       ! Rotation angle (SINALPHA and COSALPHA)
407       !
408       if (grid_type == 'C') then
409          call mprintf(.true.,STDOUT,'  Processing ROTANG')
410     
411          ! Mass-staggered points
412          allocate(sina_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
413          allocate(cosa_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
414     
415          call get_rotang(xlat_array, xlon_array, cosa_array, sina_array, &
416                          start_mem_i, start_mem_j, end_mem_i, end_mem_j)
417     
418          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'SINALPHA', &
419                           datestr, real_array = sina_array)
420          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'COSALPHA', &
421                           datestr, real_array = cosa_array)
422     
423          if (associated(sina_array)) deallocate(sina_array)
424          if (associated(cosa_array)) deallocate(cosa_array)
426          ! U-staggered points
427          allocate(sina_array(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
428          allocate(cosa_array(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
429     
430          call get_rotang(xlat_array_u, xlon_array_u, cosa_array, sina_array, &
431                          start_mem_i, start_mem_j, end_mem_stag_i, end_mem_j)
432     
433          call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'SINALPHA_U', &
434                           datestr, real_array = sina_array)
435          call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'COSALPHA_U', &
436                           datestr, real_array = cosa_array)
437     
438          if (associated(sina_array)) deallocate(sina_array)
439          if (associated(cosa_array)) deallocate(cosa_array)
440     
441          ! V-staggered points
442          allocate(sina_array(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
443          allocate(cosa_array(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
444     
445          call get_rotang(xlat_array_v, xlon_array_v, cosa_array, sina_array, &
446                          start_mem_i, start_mem_j, end_mem_i, end_mem_stag_j)
447     
448          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'SINALPHA_V', &
449                           datestr, real_array = sina_array)
450          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'COSALPHA_V', &
451                           datestr, real_array = cosa_array)
452     
453          if (associated(sina_array)) deallocate(sina_array)
454          if (associated(cosa_array)) deallocate(cosa_array)
455       end if
456     
457       ! Every field up until now should probably just be processed regardless of what the user 
458       !   has specified for fields to be processed.
459       ! Hereafter, we process user-specified fields
460     
461       !
462       ! First process the field that we will derive a landmask from
463       !
464       call get_landmask_field(geog_data_res(which_domain), landmask_name, is_water_mask, landmask_value, istatus)
466       do kk=1,MAX_LANDMASK_CATEGORIES
467          if (landmask_value(kk) == INVALID) then
468             num_landmask_categories = kk-1
469             exit
470          end if
471       end do
472       if (kk > MAX_LANDMASK_CATEGORIES) num_landmask_categories = MAX_LANDMASK_CATEGORIES
474       if (istatus /= 0) then
475          call mprintf(.true.,WARN,'No field specified for landmask calculation. Will set landmask=1 at every grid point.')
476      
477          allocate(landmask(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
478          landmask = 1.
479          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'LANDMASK', &
480                           datestr, landmask)
481     
482       else
483     
484          allocate(landmask(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
485          landmask = 1.
486      
487          call mprintf(.true.,STDOUT,'  Processing %s', s1=trim(landmask_name))
488      
489          call get_missing_fill_value(landmask_name, msg_fill_val, istatus)
490          if (istatus /= 0) msg_fill_val = NAN
491      
492          call get_halt_on_missing(landmask_name, halt_on_missing, istatus)
493          if (istatus /= 0) halt_on_missing = .false.
494      
495          ! Do we calculate a dominant category for this field?
496          call get_domcategory_name(landmask_name, domname, only_save_dominant, idomcatstatus)
497      
498          temp_string = ' '
499          temp_string(1:128) = landmask_name
500          call hash_insert(processed_fieldnames, temp_string)
501      
502          call get_max_categories(landmask_name, min_category, max_category, istatus)
503          allocate(field(start_mem_i:end_mem_i, start_mem_j:end_mem_j, min_category:max_category))
505          if (.not. only_save_dominant) then
506             field_count = field_count + 1
507             call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
508                          i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=landmask_name)
509          else
510             field_count = field_count + 1
511             call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
512                          i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
513          end if
515          if (grid_type == 'C') then
516             call calc_field(landmask_name, field, xlat_array, xlon_array, M, &
517                             start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
518                             min_category, max_category, processed_domain, 1, landmask=landmask, sr_x=1, sr_y=1)
519          else if (grid_type == 'E') then
520             call calc_field(landmask_name, field, xlat_array, xlon_array, HH, &
521                             start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
522                             min_category, max_category, processed_domain, 1, landmask=landmask, sr_x=1, sr_y=1)
523          end if
524      
525          ! If user wants to halt when a missing value is found in output field, check now
526          if (halt_on_missing) then
527             do i=start_mem_i, end_mem_i
528                do j=start_mem_j, end_mem_j
529                   ! Only need to examine k=1
530                   if (field(i,j,1) == msg_fill_val) then
531                      call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
532                   end if
533                end do
534             end do
535          end if
536      
537          ! Find fractions
538          do i=start_mem_i, end_mem_i
539             do j=start_mem_j, end_mem_j
540                sum = 0.
541                do k=min_category,max_category
542                   sum = sum + field(i,j,k)
543                end do
544                if (sum > 0.0) then
545                   do k=min_category,max_category
546                      field(i,j,k) = field(i,j,k) / sum
547                   end do
548                else
549                   do k=min_category,max_category
550                      field(i,j,k) = msg_fill_val
551                   end do
552                end if
553             end do
554          end do
555      
556          if (is_water_mask) then
557             call mprintf(.true.,STDOUT,'  Calculating landmask from %s ( WATER =', &
558                          newline=.false.,s1=trim(landmask_name))
559          else
560             call mprintf(.true.,STDOUT,'  Calculating landmask from %s ( LAND =', &
561                          newline=.false.,s1=trim(landmask_name))
562          end if
563          do k = 1, num_landmask_categories
564             call mprintf(.true.,STDOUT,' %i',newline=.false.,i1=landmask_value(k))
565             if (k == num_landmask_categories) call mprintf(.true.,STDOUT,')')
566          end do
567      
568          ! Calculate landmask
569          if (is_water_mask) then
570             do i=start_mem_i, end_mem_i
571                do j=start_mem_j, end_mem_j
572                   water_total = -1.
573                   do k=1,num_landmask_categories
574                      if (landmask_value(k) >= min_category .and. landmask_value(k) <= max_category) then 
575                         if (field(i,j,landmask_value(k)) /= msg_fill_val) then
576                            if (water_total < 0.) water_total = 0.
577                            water_total = water_total + field(i,j,landmask_value(k)) 
578                         else
579                            water_total = -1.
580                            exit
581                         end if
582                      end if
583                   end do
584                   if (water_total >= 0.0) then
585                      if (water_total < 0.50) then
586                         landmask(i,j) = 1.
587                      else
588                         landmask(i,j) = 0.
589                      end if
590                   else
591                      landmask(i,j) = -1.
592                   end if
593                end do
594             end do
595          else
596             do i=start_mem_i, end_mem_i
597                do j=start_mem_j, end_mem_j
598                   land_total = -1.
599                   do k=1,num_landmask_categories
600                      if (landmask_value(k) >= min_category .and. landmask_value(k) <= max_category) then 
601                         if (field(i,j,landmask_value(k)) /= msg_fill_val) then
602                            if (land_total < 0.) land_total = 0.
603                            land_total = land_total + field(i,j,landmask_value(k)) 
604                         else
605                            land_total = -1.
606                            exit
607                         end if
608                      end if
609                   end do
610                   if (land_total >= 0.0) then
611                      if (land_total > 0.50) then
612                         landmask(i,j) = 1.
613                      else
614                         landmask(i,j) = 0.
615                      end if
616                   else
617                      landmask(i,j) = -1.
618                   end if
619                end do
620             end do
621          end if
622     
623          call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'LANDMASK', &
624                           datestr, landmask)
625      
626          ! If we should only save the dominant category, then no need to write out fractional field
627          if (.not.only_save_dominant .or. (idomcatstatus /= 0)) then
628      
629             ! Finally, we may be asked to smooth the fractional field
630             call get_smooth_option(landmask_name, smth_opt, smth_passes, istatus)
631             if (istatus == 0) then
633                if (grid_type == 'C') then
634                   if (smth_opt == ONETWOONE) then
635                      call one_two_one(field,                      &
636                                       start_patch_i, end_patch_i, &
637                                       start_patch_j, end_patch_j, &
638                                       start_mem_i, end_mem_i,     &
639                                       start_mem_j, end_mem_j,     &
640                                       min_category, max_category, &
641                                       smth_passes, msg_fill_val)
642                   else if (smth_opt == SMTHDESMTH) then
643                      call smth_desmth(field,                      &
644                                       start_patch_i, end_patch_i, &
645                                       start_patch_j, end_patch_j, &
646                                       start_mem_i, end_mem_i,     &
647                                       start_mem_j, end_mem_j,     &
648                                       min_category, max_category, &
649                                       smth_passes, msg_fill_val)
650                   else if (smth_opt == SMTHDESMTH_SPECIAL) then
651                      call smth_desmth_special(field,              &
652                                       start_patch_i, end_patch_i, &
653                                       start_patch_j, end_patch_j, &
654                                       start_mem_i, end_mem_i,     &
655                                       start_mem_j, end_mem_j,     &
656                                       min_category, max_category, &
657                                       smth_passes, msg_fill_val)
658                   end if
659                else if (grid_type == 'E') then
660                   if (smth_opt == ONETWOONE) then
661                      call one_two_one_egrid(field,                &
662                                       start_patch_i, end_patch_i, &
663                                       start_patch_j, end_patch_j, &
664                                       start_mem_i, end_mem_i,     &
665                                       start_mem_j, end_mem_j,     &
666                                       min_category, max_category, &
667                                       smth_passes, msg_fill_val, 1.0)
668                   else if (smth_opt == SMTHDESMTH) then
669                      call smth_desmth_egrid(field,                &
670                                       start_patch_i, end_patch_i, &
671                                       start_patch_j, end_patch_j, &
672                                       start_mem_i, end_mem_i,     &
673                                       start_mem_j, end_mem_j,     &
674                                       min_category, max_category, &
675                                       smth_passes, msg_fill_val, 1.0)
676                   else if (smth_opt == SMTHDESMTH_SPECIAL) then
677                      call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
678                                               'No smoothing will be done.')
679                   end if
680                end if
682             end if
683       
684             call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
685                              min_category, max_category, trim(landmask_name), &
686                              datestr, real_array=field)
687          end if
688      
689          if (idomcatstatus == 0) then
690             allocate(dominant_field(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
691      
692             if (.not. only_save_dominant) then
693                field_count = field_count + 1
694                call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
695                             i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
696             end if
698             do i=start_mem_i, end_mem_i
699                do j=start_mem_j, end_mem_j
700                   if ((landmask(i,j) == 1. .and. is_water_mask) .or. &
701                       (landmask(i,j) == 0. .and. .not.is_water_mask)) then
702                      dominant = 0.
703                      dominant_field(i,j) = real(min_category-1)
704                      do k=min_category,max_category
705                         do kk=1,num_landmask_categories
706                            if (k == landmask_value(kk)) exit 
707                         end do
708                         if (field(i,j,k) > dominant .and. kk > num_landmask_categories) then
709                            dominant_field(i,j) = real(k)
710                            dominant = field(i,j,k)
711                         end if
712                      end do
713                   else
714                      dominant = 0.
715                      dominant_field(i,j) = real(min_category-1)
716                      do k=min_category,max_category
717                         do kk=1,num_landmask_categories
718                            if (field(i,j,k) > dominant .and. k == landmask_value(kk)) then
719                               dominant_field(i,j) = real(k)
720                               dominant = field(i,j,k)
721                            end if
722                         end do
723                      end do
724                   end if
725                end do
726             end do
727      
728             call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, trim(domname), &
729                              datestr, dominant_field)
730           
731             deallocate(dominant_field)
732          end if
733      
734          deallocate(field)
735       end if
736    
737       !
738       ! Now process all other fields specified by the user
739       !
740       call reset_next_field()
741       ifieldstatus = 0
742       do while (ifieldstatus == 0) 
743          call get_next_fieldname(fieldname, ifieldstatus)
744      
745          ! There is another field in the GEOGRID.TBL file 
746          if (ifieldstatus == 0) then
747             temp_string(1:128) = fieldname
749             call get_source_opt_status(fieldname, 0, opt_status)
750       
751             ! If this field is still to be processed
752             if (.not. hash_search(processed_fieldnames, temp_string) .and. opt_status == 0) then
753      
754                call hash_insert(processed_fieldnames, temp_string)
755                call mprintf(.true.,STDOUT,'  Processing %s', s1=trim(fieldname))
756        
757                call get_output_stagger(fieldname, istagger, istatus)
758                dimnames(:) = 'null'
759                call get_subgrid_dim_name(which_domain, fieldname, dimnames, &
760                                          sub_x, sub_y, istatus)
761        
762                if (istagger == M .or. (sub_x > 1) .or. (sub_y > 1)) then
763                   sm1 = start_mem_i
764                   em1 = end_mem_i
765                   sm2 = start_mem_j
766                   em2 = end_mem_j
767                   xlat_ptr => xlat_array 
768                   xlon_ptr => xlon_array 
769                   mapfac_ptr_x => mapfac_array_m_x
770                   mapfac_ptr_y => mapfac_array_m_y
771                else if (istagger == U) then ! In the case that extra_cols = .false.
772                   sm1 = start_mem_i          ! we should have that end_mem_stag is
773                   em1 = end_mem_stag_i       ! the same as end_mem, so we do not need
774                   sm2 = start_mem_j          ! to check extra_cols or extra rows here
775                   em2 = end_mem_j
776                   xlat_ptr => xlat_array_u 
777                   xlon_ptr => xlon_array_u 
778                   mapfac_ptr_x => mapfac_array_u_x
779                   mapfac_ptr_y => mapfac_array_u_y
780                else if (istagger == V) then
781                   sm1 = start_mem_i
782                   em1 = end_mem_i
783                   sm2 = start_mem_j
784                   em2 = end_mem_stag_j
785                   xlat_ptr => xlat_array_v 
786                   xlon_ptr => xlon_array_v 
787                   mapfac_ptr_x => mapfac_array_v_x
788                   mapfac_ptr_y => mapfac_array_v_y
789                else if (istagger == HH) then   ! E grid
790                   sm1 = start_mem_i
791                   em1 = end_mem_i
792                   sm2 = start_mem_j
793                   em2 = end_mem_j
794                   xlat_ptr => xlat_array 
795                   xlon_ptr => xlon_array 
796                   mapfac_ptr_x => mapfac_array_m_x
797                   mapfac_ptr_y => mapfac_array_m_y
798                else if (istagger == VV) then   ! E grid 
799                   sm1 = start_mem_i
800                   em1 = end_mem_i
801                   sm2 = start_mem_j
802                   em2 = end_mem_stag_j
803                   xlat_ptr => xlat_array_v 
804                   xlon_ptr => xlon_array_v 
805                   mapfac_ptr_x => mapfac_array_v_x
806                   mapfac_ptr_y => mapfac_array_v_y
807                end if
809                if (sub_x > 1) then
810                   sm1 = (start_mem_i - 1)*sub_x + 1
811                   if (extra_col) then
812                      em1 = (end_mem_i + 1)*sub_x
813                   else
814                      em1 = (end_mem_i    )*sub_x
815                   end if
816                end if
817                if (sub_y > 1)then
818                   sm2 = (start_mem_j - 1)*sub_y + 1
819                   if (extra_row) then
820                      em2 = (end_mem_j + 1)*sub_y
821                   else
822                      em2 = (end_mem_j    )*sub_y
823                   end if
824                end if
826 !BUG: This should probably be moved up to where other lat/lon fields are calculated, and we should
827 !     just determine whether we will have any subgrids or not at that point
828                if ((sub_x > 1) .or. (sub_y > 1)) then
829 !                  if (associated(xlat_array_subgrid))     deallocate(xlat_array_subgrid)
830 !                  if (associated(xlon_array_subgrid))     deallocate(xlon_array_subgrid)
831 !                  if (associated(mapfac_array_x_subgrid)) deallocate(mapfac_array_x_subgrid)
832 !                  if (associated(mapfac_array_y_subgrid)) deallocate(mapfac_array_y_subgrid)
833                   allocate(xlat_array_subgrid(sm1:em1,sm2:em2))
834                   allocate(xlon_array_subgrid(sm1:em1,sm2:em2))
835                   allocate(mapfac_array_x_subgrid(sm1:em1,sm2:em2))
836                   allocate(mapfac_array_y_subgrid(sm1:em1,sm2:em2))
837                   call get_lat_lon_fields(xlat_array_subgrid, xlon_array_subgrid, &
838                                    sm1, sm2, em1, em2, M, sub_x=sub_x, sub_y=sub_y)
839                   xlat_ptr => xlat_array_subgrid
840                   xlon_ptr => xlon_array_subgrid
841                   call get_map_factor(xlat_ptr, xlon_ptr, mapfac_array_x_subgrid, &
842                                       mapfac_array_y_subgrid, sm1, sm2, em1, em2)
843                   mapfac_ptr_x => mapfac_array_x_subgrid
844                   mapfac_ptr_y => mapfac_array_y_subgrid
845                end if
846        
847                call get_missing_fill_value(fieldname, msg_fill_val, istatus)
848                if (istatus /= 0) msg_fill_val = NAN 
849        
850                call get_halt_on_missing(fieldname, halt_on_missing, istatus)
851                if (istatus /= 0) halt_on_missing = .false.
852        
853                ! Destination field type is CONTINUOUS
854                if (iget_fieldtype(fieldname,istatus) == CONTINUOUS) then
855                   call get_max_levels(fieldname, min_level, max_level, istatus)
856                   allocate(field(sm1:em1, sm2:em2, min_level:max_level))
858                   field_count = field_count + 1
859                   call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
860                                i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=fieldname)
862                   if ((sub_x > 1) .or. (sub_y > 1)) then
863                      call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
864                                 sm1, em1, sm2, em2, min_level, max_level, &
865                                 processed_domain, 1, sr_x=sub_x, sr_y=sub_y)
866                   else
867                      call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
868                                 sm1, em1, sm2, em2, min_level, max_level, &
869                                 processed_domain, 1, landmask=landmask, sr_x=sub_x, sr_y=sub_y)
870                   end if
871         
872                   ! If user wants to halt when a missing value is found in output field, check now
873                   if (halt_on_missing) then
874                      do i=sm1, em1
875                         do j=sm2, em2
876                            ! Only need to examine k=1
877                            if (field(i,j,1) == msg_fill_val) then
878                               call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
879                            end if
880                         end do
881                      end do
882                   end if
883         
884                   ! We may be asked to smooth the fractional field
885                   call get_smooth_option(fieldname, smth_opt, smth_passes, istatus)
886                   if (istatus == 0) then
888                      if (grid_type == 'C') then
889                         if (smth_opt == ONETWOONE) then
890                            call one_two_one(field,                      &
891                                             start_patch_i, end_patch_i, &
892                                             start_patch_j, end_patch_j, &
893                                             sm1,         em1,           &
894                                             sm2,         em2,           &
895                                             min_level, max_level,       &
896                                             smth_passes, msg_fill_val)
897                         else if (smth_opt == SMTHDESMTH) then
898                            call smth_desmth(field,                      &
899                                             start_patch_i, end_patch_i, &
900                                             start_patch_j, end_patch_j, &
901                                             sm1,         em1,           &
902                                             sm2,         em2,           &
903                                             min_level, max_level,       &
904                                             smth_passes, msg_fill_val)
905                         else if (smth_opt == SMTHDESMTH_SPECIAL) then
906                            call smth_desmth_special(field,              &
907                                             start_patch_i, end_patch_i, &
908                                             start_patch_j, end_patch_j, &
909                                             sm1,         em1,           &
910                                             sm2,         em2,           &
911                                             min_level, max_level,       &
912                                             smth_passes, msg_fill_val)
913                        end if
914   
915                      else if (grid_type == 'E') then
917                         if (trim(fieldname) == 'HGT_M' ) then
918                            topo_flag_val=1.0
919                            mass_flag=1.0
920                         else if (trim(fieldname) == 'HGT_V') then
921                            topo_flag_val=1.0
922                            mass_flag=0.0
923                         else
924                            topo_flag_val=0.0
925                         end if
926   
927                         if (smth_opt == ONETWOONE) then
928                            call one_two_one_egrid(field,                &
929                                             start_patch_i, end_patch_i, &
930                                             start_patch_j, end_patch_j, &
931                                             sm1,         em1,           &
932                                             sm2,         em2,           &
933                                             min_level, max_level,       &
934                                             smth_passes, topo_flag_val, mass_flag)
935                         else if (smth_opt == SMTHDESMTH) then
936                            call smth_desmth_egrid(field,                &
937                                             start_patch_i, end_patch_i, &
938                                             start_patch_j, end_patch_j, &
939                                             sm1,         em1,           &
940                                             sm2,         em2,           &
941                                             min_level, max_level,       &
942                                             smth_passes, topo_flag_val, mass_flag)
943                         else if (smth_opt == SMTHDESMTH_SPECIAL) then
944                            call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
945                                                     'No smoothing will be done.')
946                         end if
947   
948                      end if
950                   end if
952                   call write_field(sm1, em1, sm2, em2, &
953                                    min_level, max_level, trim(fieldname), datestr, real_array=field)
954         
955                   ! Do we calculate directional derivatives from this field?
956                   call get_dfdx_name(fieldname, gradname, istatus)
957                   if (istatus == 0) then
958                      allocate(slp_field(sm1:em1,sm2:em2,min_level:max_level))
960                      field_count = field_count + 1
961                      call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
962                                   i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=gradname)
964                      if (grid_type == 'C') then
965                         call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, mapfac_ptr_x)
966                      else if (grid_type == 'E') then
967                         call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level)
968                      end if
969                      call write_field(sm1, em1, sm2, em2, &
970                                       min_level, max_level, trim(gradname), datestr, real_array=slp_field)
971                      deallocate(slp_field)
972                   end if
973                   call get_dfdy_name(fieldname, gradname, istatus)
974                   if (istatus == 0) then
975                      allocate(slp_field(sm1:em1,sm2:em2,min_level:max_level))
977                      field_count = field_count + 1
978                      call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
979                                   i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=gradname)
981                      if (grid_type == 'C') then
982                         call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, mapfac_ptr_y)
983                      else if (grid_type == 'E') then
984                         call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level)
985                      end if
986                      call write_field(sm1, em1, sm2, em2, &
987                                       min_level, max_level, trim(gradname), datestr, real_array=slp_field)
988                      deallocate(slp_field)
989                   end if
990         
991                   deallocate(field)
992       
993                ! Destination field type is CATEGORICAL
994                else
995                   call get_max_categories(fieldname, min_category, max_category, istatus)
996                   allocate(field(sm1:em1, sm2:em2, min_category:max_category))
998                   ! Do we calculate a dominant category for this field?
999                   call get_domcategory_name(fieldname, domname, only_save_dominant, idomcatstatus)
1000         
1001                   if (.not. only_save_dominant) then
1002                      field_count = field_count + 1
1003                      call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
1004                                   i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=fieldname)
1005                   else
1006                      field_count = field_count + 1
1007                      call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
1008                                   i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
1009                   end if
1011                   if ((sub_x > 1) .or. (sub_y > 1)) then
1012                      call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
1013                                      sm1, em1, sm2, em2, min_category, max_category, &
1014                                      processed_domain, 1, sr_x=sub_x, sr_y=sub_y)
1015                   else
1016                      call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
1017                                      sm1, em1, sm2, em2, min_category, max_category, &
1018                                      processed_domain, 1, landmask=landmask, sr_x=sub_x, sr_y=sub_y)
1019                   end if
1020         
1021                   ! If user wants to halt when a missing value is found in output field, check now
1022                   if (halt_on_missing) then
1023                      do i=sm1, em1
1024                         do j=sm2, em2
1025                            ! Only need to examine k=1
1026                            if (field(i,j,1) == msg_fill_val) then
1027                               call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
1028                            end if
1029                         end do
1030                      end do
1031                   end if
1032        
1033                   ! Find fractions
1034                   do i=sm1, em1
1035                      do j=sm2, em2
1036                         sum = 0.
1037                         do k=min_category,max_category
1038                            sum = sum + field(i,j,k)
1039                         end do
1040                         if (sum > 0.0) then
1041                            do k=min_category,max_category
1042                               field(i,j,k) = field(i,j,k) / sum
1043                            end do
1044                         else
1045                            do k=min_category,max_category
1046                               field(i,j,k) = msg_fill_val
1047                            end do
1048                         end if
1049                      end do
1050                   end do
1051         
1052                   ! If we should only save the dominant category, then no need to write out fractional field
1053                   if (.not.only_save_dominant .or. (idomcatstatus /= 0)) then
1054         
1055                      ! Finally, we may be asked to smooth the fractional field
1056                      call get_smooth_option(fieldname, smth_opt, smth_passes, istatus)
1057                      if (istatus == 0) then
1058                         if (grid_type == 'C') then
1059                            if (smth_opt == ONETWOONE) then
1060                               call one_two_one(field,                    &
1061                                              start_patch_i, end_patch_i, &
1062                                              start_patch_j, end_patch_j, &
1063                                              sm1,         em1,           &
1064                                              sm2,         em2,           &
1065                                              min_category, max_category, &
1066                                              smth_passes, msg_fill_val)
1067                            else if (smth_opt == SMTHDESMTH) then
1068                               call smth_desmth(field,                    &
1069                                              start_patch_i, end_patch_i, &
1070                                              start_patch_j, end_patch_j, &
1071                                              sm1,         em1,           &
1072                                              sm2,         em2,           &
1073                                              min_category, max_category, &
1074                                              smth_passes, msg_fill_val)
1075                            else if (smth_opt == SMTHDESMTH_SPECIAL) then
1076                               call smth_desmth_special(field,            &
1077                                              start_patch_i, end_patch_i, &
1078                                              start_patch_j, end_patch_j, &
1079                                              sm1,         em1,           &
1080                                              sm2,         em2,           &
1081                                              min_category, max_category, &
1082                                              smth_passes, msg_fill_val)
1083                            end if
1084                         else if (grid_type == 'E') then
1085                            if (smth_opt == ONETWOONE) then
1086                               call one_two_one_egrid(field,              &
1087                                              start_patch_i, end_patch_i, &
1088                                              start_patch_j, end_patch_j, &
1089                                              sm1,         em1,           &
1090                                              sm2,         em2,           &
1091                                              min_category, max_category, &
1092                                              smth_passes, msg_fill_val, 1.0)
1093                            else if (smth_opt == SMTHDESMTH) then
1094                               call smth_desmth_egrid(field,              &
1095                                              start_patch_i, end_patch_i, &
1096                                              start_patch_j, end_patch_j, &
1097                                              sm1,         em1,           &
1098                                              sm2,         em2,           &
1099                                              min_category, max_category, &
1100                                              smth_passes, msg_fill_val, 1.0)
1101                            else if (smth_opt == SMTHDESMTH_SPECIAL) then
1102                               call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
1103                                                        'No smoothing will be done.')
1104                            end if
1105                         end if
1106                      end if
1107          
1108                      call write_field(sm1, em1, sm2, em2, &
1109                                       min_category, max_category, trim(fieldname), datestr, real_array=field)
1110                   end if
1111         
1112                   if (idomcatstatus == 0) then
1113                      call mprintf(.true.,STDOUT,'  Processing %s', s1=trim(domname))
1114                      allocate(dominant_field(sm1:em1, sm2:em2))
1116                      if (.not. only_save_dominant) then
1117                         field_count = field_count + 1
1118                         call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
1119                                      i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
1120                      end if
1122                      do i=sm1, em1
1123                         do j=sm2, em2
1124                            dominant = 0.
1125                            dominant_field(i,j) = real(min_category-1)
1126                            do k=min_category,max_category
1127                               if (field(i,j,k) > dominant .and. field(i,j,k) /= msg_fill_val) then 
1128                                  dominant_field(i,j) = real(k)
1129                                  dominant = field(i,j,k)
1130              !                  else
1131              !                    dominant_field(i,j) = nint(msg_fill_val)
1132 ! Maybe we should put an else clause here to set the category equal to the missing fill value?
1133 ! BUG: The problem here seems to be that, when we set a fraction equal to the missing fill value
1134 !   above, if the last fractional index we process here has been filled, we think that the dominant
1135 !   category should be set to the missing fill value. Perhaps we could do some check to only
1136 !   assign the msg_fill_val if no other valid category has been assigned? But this may still not
1137 !   work if the missing fill value is something like 0.5. Somehow use bitarrays, perhaps, to remember
1138 !   which points are missing and which just happen to have the missing fill value?
1139                               end if
1140                            end do
1141                            if (dominant_field(i,j) == real(min_category-1)) dominant_field(i,j) = msg_fill_val
1142                         end do
1143                      end do
1144                      call write_field(sm1, em1, sm2, em2, 1, 1, &
1145                                       trim(domname), datestr, dominant_field)
1146                      deallocate(dominant_field)
1147                   end if
1148        
1149                   deallocate(field)
1151                   if ((sub_x > 1) .or. (sub_y > 1)) then
1152                      if (associated(xlat_array_subgrid))     deallocate(xlat_array_subgrid)
1153                      if (associated(xlon_array_subgrid))     deallocate(xlon_array_subgrid)
1154                      if (associated(mapfac_array_x_subgrid)) deallocate(mapfac_array_x_subgrid)
1155                      if (associated(mapfac_array_y_subgrid)) deallocate(mapfac_array_y_subgrid)
1156                   end if
1157    
1158                end if
1159       
1160             end if
1161          end if
1162     
1163       end do
1165       ! Close output
1166       call output_close()
1167     
1168       call hash_destroy(processed_fieldnames)
1169     
1170       ! Free up memory
1171       if (associated(xlat_array)) deallocate(xlat_array)
1172       if (associated(xlon_array)) deallocate(xlon_array)
1173       if (grid_type == 'C') then
1174          if (associated(xlat_array_u)) deallocate(xlat_array_u)
1175          if (associated(xlon_array_u)) deallocate(xlon_array_u)
1176          if (associated(xlat_array_corner)) deallocate(xlat_array_corner)
1177          if (associated(xlon_array_corner)) deallocate(xlon_array_corner)
1178          if (associated(mapfac_array_u_x)) deallocate(mapfac_array_u_x)
1179          if (associated(mapfac_array_u_y)) deallocate(mapfac_array_u_y)
1180       end if
1181       if (associated(xlat_array_v)) deallocate(xlat_array_v)
1182       if (associated(xlon_array_v)) deallocate(xlon_array_v)
1183       if (associated(mapfac_array_m_x)) deallocate(mapfac_array_m_x)
1184       if (associated(mapfac_array_m_y)) deallocate(mapfac_array_m_y)
1185       if (associated(mapfac_array_v_x)) deallocate(mapfac_array_v_x)
1186       if (associated(mapfac_array_v_y)) deallocate(mapfac_array_v_y)
1187       if (associated(landmask)) deallocate(landmask)
1188       if (associated(xlat_array_subgrid))     deallocate(xlat_array_subgrid)
1189       if (associated(xlon_array_subgrid))     deallocate(xlon_array_subgrid)
1190       if (associated(mapfac_array_x_subgrid)) deallocate(mapfac_array_x_subgrid)
1191       if (associated(mapfac_array_y_subgrid)) deallocate(mapfac_array_y_subgrid)
1193       nullify(xlat_ptr)
1194       nullify(xlon_ptr)
1195    
1196    end subroutine process_tile
1197    
1198    
1199    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1200    ! Name: calc_field
1201    !
1202    ! Purpose: This routine fills in the "field" array with interpolated source 
1203    !   data. When multiple resolutions of source data are available, an appropriate
1204    !   resolution is chosen automatically. The specified field may either be a 
1205    !   continuous field or a categorical field.
1206    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1207    recursive subroutine calc_field(fieldname, field, xlat_array, xlon_array, istagger, &
1208                                    start_i, end_i, start_j, end_j, start_k, end_k, &
1209                                    processed_domain, ilevel, landmask, sr_x, sr_y)
1210    
1211       use bitarray_module
1212       use interp_module
1213       use llxy_module
1214       use misc_definitions_module
1215       use proc_point_module
1216       use queue_module
1217       use source_data_module
1218     
1219       implicit none
1220     
1221       ! Arguments
1222       integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, ilevel, istagger
1223       real, dimension(start_i:end_i, start_j:end_j), intent(in) :: xlat_array, xlon_array
1224       real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: field
1225       real, dimension(start_i:end_i, start_j:end_j), intent(in), optional :: landmask
1226       integer, intent(in), optional :: sr_x, sr_y
1227       character (len=128), intent(in) :: fieldname
1228       type (bitarray), intent(inout) :: processed_domain
1229     
1230       ! Local variables
1231       integer :: start_src_k, end_src_k
1232       integer :: i, j, k, ix, iy, itype
1233       integer :: user_iproj, istatus
1234       integer :: opt_status
1235       real :: mask_val
1236       real :: temp
1237       real :: scale_factor
1238       real :: msg_val, msg_fill_val, threshold, src_dx, src_dy, dom_dx, dom_dy
1239       real :: user_stand_lon, user_truelat1, user_truelat2, user_dxkm, user_dykm, &
1240               user_known_x, user_known_y, user_known_lat, user_known_lon
1241       real, pointer, dimension(:,:,:) :: data_count
1242       integer, pointer, dimension(:) :: interp_type
1243       integer, pointer, dimension(:) :: interp_opts
1244       character (len=128) :: interp_string
1245       type (bitarray) :: bit_domain, level_domain
1246       type (queue)    :: point_queue, tile_queue
1247       type (q_data)   :: current_pt
1249       nullify(data_count)
1250       nullify(interp_type)
1251       nullify(interp_opts)
1253       ! If this is the first trip through this routine, we need to allocate the bit array that
1254       !  will persist through all recursive calls, tracking which grid points have been assigned
1255       !  a value. 
1256       if (ilevel == 1) call bitarray_create(processed_domain, end_i-start_i+1, end_j-start_j+1)
1258       ! Find out if this "priority level" (given by ilevel) exists
1259       call check_priority_level(fieldname, ilevel, istatus)
1260     
1261       ! A bad status indicates that that no data for priority level ilevel is available, and thus, that
1262       !   no further levels will be specified. We are done processing for this level.
1263       if (istatus /= 0) then
1264          if (ilevel == 1) call bitarray_destroy(processed_domain)
1265          return 
1266       end if
1267     
1268       ! Before proceeding with processing for this level, though, process for the next highest priority level
1269       !   of source data
1270       call calc_field(fieldname, field, xlat_array, xlon_array, istagger, start_i, end_i, &
1271                       start_j, end_j, start_k, end_k, processed_domain, ilevel+1, landmask, sr_x, sr_y)
1273       ! At this point, all levels of source data with higher priority have been processed, and we can assign
1274       !   values to all grid points that have not already been given values from higher-priority data
1276       call get_source_opt_status(fieldname, ilevel, opt_status)
1277       if (opt_status == 0) then
1279          ! Find out the projection of the data for this "priority level" (given by ilevel)
1280          call get_data_projection(fieldname, user_iproj, user_stand_lon, user_truelat1, user_truelat2, &
1281                                   user_dxkm, user_dykm, user_known_x, user_known_y, user_known_lat, &
1282                                   user_known_lon, ilevel, istatus)
1283        
1284          ! A good status indicates that there is data for this priority level, so we store the projection
1285          !   of that data on a stack. The projection will be on the top of the stack (and hence will be 
1286          !   the "active" projection) once all higher-priority levels have been processed
1287          call push_source_projection(user_iproj, user_stand_lon, user_truelat1, user_truelat2, &
1288                                   user_dxkm, user_dykm, user_dykm, user_dxkm, user_known_x, user_known_y, &
1289                                   user_known_lat, user_known_lon)
1290        
1291          ! Initialize point processing module
1292          call proc_point_init()
1293        
1294          ! Initialize queues
1295          call q_init(point_queue)
1296          call q_init(tile_queue)
1297        
1298          ! Determine whether we will be processing categorical data or continuous data
1299          itype = iget_source_fieldtype(fieldname, ilevel, istatus)
1300          call get_interp_option(fieldname, ilevel, interp_string, istatus)
1301          interp_type => interp_array_from_string(interp_string)
1302          interp_opts => interp_options_from_string(interp_string)
1303    
1304          ! Also, check whether we will be using the cell averaging interpolator for continuous fields
1305          if (index(interp_string,'average_gcell') /= 0 .and. itype == CONTINUOUS) then
1306             call get_gcell_threshold(interp_string, threshold, istatus)
1307             if (istatus == 0) then
1308                call get_source_resolution(fieldname, ilevel, src_dx, src_dy, istatus)
1309                if (istatus == 0) then
1310                   call get_domain_resolution(dom_dx, dom_dy)
1311                   if (gridtype == 'C') then
1312                      if (threshold*max(src_dx,src_dy)*111. <= max(dom_dx,dom_dy)/1000.) then
1313                         itype = SP_CONTINUOUS
1314                         allocate(data_count(start_i:end_i,start_j:end_j,start_k:end_k))
1315                         data_count = 0.
1316                      end if
1317                   else if (gridtype == 'E') then
1318                      if (max(src_dx,src_dy) >= threshold*max(dom_dx,dom_dy)) then
1319                         itype = SP_CONTINUOUS
1320                         allocate(data_count(start_i:end_i,start_j:end_j,start_k:end_k))
1321                         data_count = 0.
1322                      end if
1323                   end if
1324                end if
1325             end if
1326          end if
1327    
1328          call get_missing_value(fieldname, ilevel, msg_val, istatus)
1329          if (istatus /= 0) msg_val = NAN
1330          call get_missing_fill_value(fieldname, msg_fill_val, istatus)
1331          if (istatus /= 0) msg_fill_val = NAN
1332        
1333          call get_masked_value(fieldname, ilevel, mask_val, istatus)
1334          if (istatus /= 0) mask_val = -1.
1335        
1336          if (itype == CONTINUOUS .or. itype == SP_CONTINUOUS) then
1337             call get_source_levels(fieldname, ilevel, start_src_k, end_src_k, istatus)
1338             if (istatus /= 0) return
1339          end if
1340        
1341          ! Initialize bitarray used to track which points have been visited and assigned values while 
1342          !   processing *this* priority level of data
1343          call bitarray_create(bit_domain, end_i-start_i+1, end_j-start_j+1)
1344          call bitarray_create(level_domain, end_i-start_i+1, end_j-start_j+1)
1345        
1346          ! Begin by placing a point in the tile_queue
1347          current_pt%lat = xlat_array(start_i,start_j)
1348          current_pt%lon = xlon_array(start_i,start_j)
1349          current_pt%x = start_i
1350          current_pt%y = start_j
1351          call q_insert(tile_queue, current_pt)
1352        
1353          ! While there are still grid points in tiles that have not yet been processed
1354          do while (q_isdata(tile_queue))
1355        
1356             ! Take a point from the outer queue and place it in the point_queue for processing 
1357             current_pt = q_remove(tile_queue)
1358         
1359             ! If this level of data is categorical (i.e., is given as an array of category indices),
1360             !   then first try to process the entire tile in one call to accum_categorical. Any grid
1361             !   points that are not given values by accum_categorical and that lie within the current
1362             !   tile of source data are individually assigned values in the inner loop
1363             if (itype == CATEGORICAL) then
1364         
1365                ! Have we already visited this point? If so, this tile has already been processed by 
1366                !   accum_categorical. 
1367                if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
1368                   call q_insert(point_queue, current_pt) 
1369                   if (.not. have_processed_tile(current_pt%lat, current_pt%lon, fieldname, ilevel)) then
1370                      call accum_categorical(current_pt%lat, current_pt%lon, istagger, field, &
1371                                             start_i, end_i, start_j, end_j, start_k, end_k, &
1372                                             fieldname, processed_domain, level_domain, &
1373                                             ilevel, msg_val, mask_val, sr_x, sr_y)
1374 ! BUG: Where do we mask out those points that are on land/water when masked=land/water is set? 
1375                   end if
1376                   call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
1377                end if
1378         
1379             else if (itype == SP_CONTINUOUS) then
1380    
1381                ! Have we already visited this point? If so, this tile has already been processed by 
1382                !   accum_continuous. 
1383                if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
1384                   call q_insert(point_queue, current_pt) 
1385                   if (.not. have_processed_tile(current_pt%lat, current_pt%lon, fieldname, ilevel)) then
1386                      call accum_continuous(current_pt%lat, current_pt%lon, istagger, field, data_count, &
1387                                             start_i, end_i, start_j, end_j, start_k, end_k, &
1388                                             fieldname, processed_domain, level_domain, &
1389                                             ilevel, msg_val, mask_val, sr_x, sr_y)
1390 ! BUG: Where do we mask out those points that are on land/water when masked=land/water is set? 
1391                   end if
1392                   call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
1393                end if
1394    
1395             else if (itype == CONTINUOUS) then
1396         
1397                ! Have we already visited this point? If so, the tile containing this point has already been
1398                !   processed in the inner loop.
1399                if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
1400                   call q_insert(point_queue, current_pt) 
1401                   call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
1402                end if
1403         
1404             end if
1405         
1406             ! This inner loop, where all grid points contained in the current source tile are processed
1407             do while (q_isdata(point_queue))
1408                current_pt = q_remove(point_queue)
1409                ix = current_pt%x
1410                iy = current_pt%y
1411          
1412                ! Process the current point
1413                if (itype == CONTINUOUS .or. itype == SP_CONTINUOUS) then
1414          
1415                   ! Have we already assigned this point a value from this priority level?
1416                   if (.not. bitarray_test(level_domain, ix-start_i+1, iy-start_j+1)) then
1417            
1418                      ! If the point was already assigned a value from a higher-priority level, no 
1419                      !   need to assign a new value
1420                      if (bitarray_test(processed_domain, ix-start_i+1, iy-start_j+1)) then
1421                         call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1422             
1423                      ! Otherwise, need to assign values from this level of source data if we can
1424                      else
1425                         if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1426                            if (landmask(ix,iy) /= mask_val) then
1427                               do k=start_src_k,end_src_k
1428                                  temp = get_point(current_pt%lat, current_pt%lon, k, &
1429                                                   fieldname, ilevel, interp_type, interp_opts, msg_val)
1430                                  if (temp /= msg_val) then
1431                                     field(ix, iy, k) = temp
1432                                     call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1433                                     if (itype == SP_CONTINUOUS) data_count(ix, iy, k) = 1.0
1434                                  else
1435                                     field(ix, iy, k) = msg_fill_val
1436                                  end if
1437                               end do
1438                            else
1439                               do k=start_k,end_k
1440                                  field(ix,iy,k) = msg_fill_val
1441                               end do
1442                            end if
1443                         else
1444                            do k=start_src_k,end_src_k
1445                               temp = get_point(current_pt%lat, current_pt%lon, k, &
1446                                                fieldname, ilevel, interp_type, interp_opts, msg_val)
1447                               if (temp /= msg_val) then
1448                                  field(ix, iy, k) = temp
1449                                  call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1450                                  if (itype == SP_CONTINUOUS) data_count(ix, iy, k) = 1.0
1451                               else
1452                                  field(ix, iy, k) = msg_fill_val
1453                               end if
1454                            end do
1455                         end if
1456                      end if
1457                   end if
1458          
1459                else if (itype == CATEGORICAL) then
1460          
1461                   ! Have we already assigned this point a value from this priority level?
1462                   if (.not.bitarray_test(level_domain, ix-start_i+1, iy-start_j+1)) then
1463           
1464                      ! If the point was already assigned a value from a higher-priority level, no 
1465                      !   need to assign a new value
1466                      if (bitarray_test(processed_domain, ix-start_i+1, iy-start_j+1)) then
1467                         call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1468             
1469                      ! Otherwise, the point was apparently not given a value when accum_categorical
1470                      !   was called for the current tile, and we need to assign values from this 
1471                      !   level of source data if we can
1472                      else
1473                         if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1474                            if (landmask(ix,iy) /= mask_val) then
1475                               temp = get_point(current_pt%lat, current_pt%lon, 1, &
1476                                                fieldname, ilevel, interp_type, interp_opts, msg_val)
1477          
1478                               do k=start_k,end_k
1479                                  field(ix,iy,k) = 0.
1480                               end do
1481          
1482                               if (temp /= msg_val) then
1483                                  if (int(temp) >= start_k .and. int(temp) <= end_k) then
1484                                     field(ix, iy, int(temp)) = field(ix, iy, int(temp)) + 1.
1485                                  else
1486                                     call mprintf(.true.,WARN,' Attempted to assign an invalid category '// &
1487                                                  '%i to grid point (%i, %i)', i1=int(temp), i2=ix, i3=iy)
1488                                  end if
1489                                  call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1490                               end if
1491       
1492                            else
1493                               do k=start_k,end_k
1494                                  field(ix,iy,k) = 0.
1495                               end do
1496                            end if
1497                         else
1498                            temp = get_point(current_pt%lat, current_pt%lon, 1, &
1499                                             fieldname, ilevel, interp_type, interp_opts, msg_val)
1500          
1501                            do k=start_k,end_k
1502                               field(ix,iy,k) = 0.
1503                            end do
1504          
1505                            if (temp /= msg_val) then
1506                               if (int(temp) >= start_k .and. int(temp) <= end_k) then
1507                                  field(ix, iy, int(temp)) = field(ix, iy, int(temp)) + 1.
1508                               else
1509                                  call mprintf(.true.,WARN,' Attempted to assign an invalid category '// &
1510                                               '%i to grid point (%i, %i)', i1=int(temp), i2=ix, i3=iy)
1511                               end if
1512                               call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1513                            end if
1514                         end if
1515                      end if
1516                   end if
1517          
1518                end if
1519          
1520                ! Scan neighboring points, adding them to the appropriate queue based on whether they
1521                !   are in the current tile or not
1522                if (iy > start_j) then
1523                   if (ix > start_i) then
1524            
1525                      ! Neighbor with relative position (-1,-1)
1526                      call process_neighbor(ix-1, iy-1, bit_domain, point_queue, tile_queue, &
1527                                            xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1528                   end if
1529           
1530                   ! Neighbor with relative position (0,-1) 
1531                   call process_neighbor(ix, iy-1, bit_domain, point_queue, tile_queue, &
1532                                         xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1533           
1534                   if (ix < end_i) then
1535           
1536                      ! Neighbor with relative position (+1,-1)
1537                      call process_neighbor(ix+1, iy-1, bit_domain, point_queue, tile_queue, &
1538                                            xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1539                   end if
1540                end if
1541          
1542                if (ix > start_i) then
1543          
1544                   ! Neighbor with relative position (-1,0)
1545                   call process_neighbor(ix-1, iy, bit_domain, point_queue, tile_queue, &
1546                                         xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1547                end if
1548          
1549                if (ix < end_i) then
1550          
1551                   ! Neighbor with relative position (+1,0)
1552                   call process_neighbor(ix+1, iy, bit_domain, point_queue, tile_queue, &
1553                                         xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1554                end if
1555          
1556                if (iy < end_j) then
1557                   if (ix > start_i) then
1558           
1559                      ! Neighbor with relative position (-1,+1)
1560                      call process_neighbor(ix-1, iy+1, bit_domain, point_queue, tile_queue, &
1561                                            xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1562                   end if
1563           
1564                   ! Neighbor with relative position (0,+1)
1565                   call process_neighbor(ix, iy+1, bit_domain, point_queue, tile_queue, &
1566                                         xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1567                   if (ix < end_i) then
1568           
1569                      ! Neighbor with relative position (+1,+1)
1570                      call process_neighbor(ix+1, iy+1, bit_domain, point_queue, tile_queue, &
1571                                            xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1572                   end if
1573                end if
1574          
1575             end do
1576         
1577          end do
1578    
1579          if (itype == SP_CONTINUOUS) then
1580             itype = CONTINUOUS
1581             if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1582                do j=start_j,end_j
1583                   do i=start_i,end_i
1584                      if (landmask(i,j) /= mask_val) then
1585                         do k=start_k,end_k
1586                            if (data_count(i,j,k) > 0.) then
1587                               field(i,j,k) = field(i,j,k) / data_count(i,j,k)
1588                            else
1589                               if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
1590                                  field(i,j,k) = msg_fill_val
1591                               end if
1592                            end if
1593                         end do
1594                      else
1595                         if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
1596                            do k=start_k,end_k
1597                               field(i,j,k) = msg_fill_val
1598                            end do
1599                         end if
1600                      end if
1601                   end do
1602                end do
1603             else
1604                do k=start_k,end_k
1605                   do j=start_j,end_j
1606                      do i=start_i,end_i
1607                         if (data_count(i,j,k) > 0.) then
1608                            field(i,j,k) = field(i,j,k) / data_count(i,j,k)
1609                         else
1610                            if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
1611                               field(i,j,k) = msg_fill_val
1612                            end if
1613                         end if
1614                      end do
1615                   end do
1616                end do
1617             end if
1618             deallocate(data_count)
1619    
1620          else if (itype == CATEGORICAL) then
1621             if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1622                do j=start_j,end_j
1623                   do i=start_i,end_i
1624                      if (landmask(i,j) == mask_val) then
1625                         do k=start_k,end_k
1626                            field(i,j,k) = 0.
1627                         end do
1628                      end if
1629                   end do
1630                end do
1631             end if
1632          end if
1633    
1634          deallocate(interp_type)
1635          deallocate(interp_opts)
1636    
1637    
1638          ! We may need to scale this field by a constant
1639          call get_field_scale_factor(fieldname, ilevel, scale_factor, istatus)
1640          if (istatus == 0) then
1641             do i=start_i, end_i
1642                do j=start_j, end_j
1643                   if (bitarray_test(level_domain,i-start_i+1,j-start_j+1) .and. &
1644                       .not. bitarray_test(processed_domain,i-start_i+1,j-start_j+1)) then
1645                      do k=start_k,end_k
1646                         if (field(i,j,k) /= msg_fill_val) then
1647                            field(i,j,k) = field(i,j,k) * scale_factor
1648                         end if
1649                      end do
1650                   end if
1651                end do
1652             end do
1653          end if
1654    
1655        
1656          ! Now add the points that were assigned values at this priority level to the complete array
1657          !   of points that have been assigned values
1658          call bitarray_merge(processed_domain, level_domain)
1659        
1660          call bitarray_destroy(bit_domain)
1661          call bitarray_destroy(level_domain)
1662          call q_destroy(point_queue)
1663          call q_destroy(tile_queue)
1664          call proc_point_shutdown()
1665        
1666          ! Remove the projection of the current level of source data from the stack, thus "activating" 
1667          !   the projection of the next highest level
1668          call pop_source_projection()
1670       else
1671          call mprintf(.true.,STDOUT,'   Important note: could not open input dataset for priority level %i, '// &
1672                                     'but this level is optional.', i1=ilevel)
1673          call mprintf(.true.,LOGFILE,'   Important note: could not open input dataset for priority level %i, '// &
1674                                     'but this level is optional.', i1=ilevel)
1675       end if
1677       ! If this is the last level of the recursion, we can also deallocate processed_domain
1678       if (ilevel == 1) call bitarray_destroy(processed_domain)
1679    
1680    end subroutine calc_field
1681    
1682    
1683    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1684    ! Name: get_lat_lon_fields
1685    !
1686    ! Purpose: To calculate the latitude and longitude for every gridpoint in the
1687    !   tile of the model domain. The caller may specify that the grid for which 
1688    !   values are computed is staggered or unstaggered using the "stagger" 
1689    !   argument.
1690    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1691    subroutine get_lat_lon_fields(xlat_arr, xlon_arr, start_mem_i, &
1692                                  start_mem_j, end_mem_i, end_mem_j, stagger, comp_ll, &
1693                                  sub_x, sub_y)
1694    
1695       use llxy_module
1696       use misc_definitions_module
1697     
1698       implicit none
1699     
1700       ! Arguments
1701       integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, &
1702                              end_mem_j, stagger
1703       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: xlat_arr, xlon_arr
1704       logical, optional, intent(in) :: comp_ll
1705       integer, optional, intent(in) :: sub_x, sub_y
1707       ! Local variables
1708       integer :: i, j
1709       real :: rx, ry
1710     
1711       rx = 1.0
1712       ry = 1.0
1713       if (present(sub_x)) rx = real(sub_x)
1714       if (present(sub_y)) ry = real(sub_y)
1716       do i=start_mem_i, end_mem_i
1717          do j=start_mem_j, end_mem_j
1718             call xytoll(real(i-1)/rx+1.0, real(j-1)/ry+1.0, &
1719                         xlat_arr(i,j), xlon_arr(i,j), stagger, comp_ll=comp_ll)
1720          end do
1721       end do
1723    end subroutine get_lat_lon_fields
1724    
1726    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1727    ! Name: get_map_factor
1728    !
1729    ! Purpose: Given the latitude field, this routine calculates map factors for 
1730    !   the grid points of the specified domain. For different grids (e.g., C grid, 
1731    !   E grid), the latitude array should provide the latitudes of the points for
1732    !   which map factors are to be calculated. 
1733    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1734    subroutine get_map_factor(xlat_arr, xlon_arr, mapfac_arr_x, mapfac_arr_y, &
1735                              start_mem_i, start_mem_j, end_mem_i, end_mem_j)
1736    
1737       use constants_module
1738       use gridinfo_module
1739       use misc_definitions_module
1740       use map_utils
1741     
1742       implicit none
1743     
1744       ! Arguments
1745       integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
1746       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr, xlon_arr
1747       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: mapfac_arr_x
1748       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: mapfac_arr_y
1749     
1750       ! Local variables
1751       integer :: i, j
1752       real :: n, colat, colat0, colat1, colat2, comp_lat, comp_lon
1753     
1754       !
1755       ! Equations for map factor given in Principles of Meteorological Analysis,
1756       ! Walter J. Saucier, pp. 32-33 
1757       !
1758     
1759       ! Lambert conformal projection
1760       if (iproj_type == PROJ_LC) then
1761          if (truelat1 /= truelat2) then
1762             colat1 = rad_per_deg*(90.0 - truelat1)
1763             colat2 = rad_per_deg*(90.0 - truelat2)
1764             n = (log(sin(colat1)) - log(sin(colat2))) &
1765                 / (log(tan(colat1/2.0)) - log(tan(colat2/2.0)))
1766       
1767             do i=start_mem_i, end_mem_i
1768                do j=start_mem_j, end_mem_j
1769                   colat = rad_per_deg*(90.0 - xlat_arr(i,j))
1770                   mapfac_arr_x(i,j) = sin(colat2)/sin(colat)*(tan(colat/2.0)/tan(colat2/2.0))**n
1771                   mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1772                end do
1773             end do
1774      
1775          else
1776             colat0 = rad_per_deg*(90.0 - truelat1)
1777       
1778             do i=start_mem_i, end_mem_i
1779                do j=start_mem_j, end_mem_j
1780                   colat = rad_per_deg*(90.0 - xlat_arr(i,j))
1781                   mapfac_arr_x(i,j) = sin(colat0)/sin(colat)*(tan(colat/2.0)/tan(colat0/2.0))**cos(colat0)
1782                   mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1783                end do
1784             end do
1785     
1786          end if
1787     
1788       ! Polar stereographic projection
1789       else if (iproj_type == PROJ_PS) then
1790     
1791          do i=start_mem_i, end_mem_i
1792             do j=start_mem_j, end_mem_j
1793                mapfac_arr_x(i,j) = (1.0 + sin(rad_per_deg*abs(truelat1)))/(1.0 + sin(rad_per_deg*sign(1.,truelat1)*xlat_arr(i,j)))
1794                mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1795             end do
1796          end do
1797     
1798       ! Mercator projection 
1799       else if (iproj_type == PROJ_MERC) then
1800          colat0 = rad_per_deg*(90.0 - truelat1)
1801      
1802          do i=start_mem_i, end_mem_i
1803             do j=start_mem_j, end_mem_j
1804                colat = rad_per_deg*(90.0 - xlat_arr(i,j))
1805                mapfac_arr_x(i,j) = sin(colat0) / sin(colat) 
1806                mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1807             end do
1808          end do
1809     
1810       ! Global cylindrical projection
1811       else if (iproj_type == PROJ_CYL) then
1812      
1813          do i=start_mem_i, end_mem_i
1814             do j=start_mem_j, end_mem_j
1815                if (abs(xlat_arr(i,j)) == 90.0) then
1816                   mapfac_arr_x(i,j) = 0.    ! MSF actually becomes infinite at poles, but 
1817                                             !   the values should never be used there; by
1818                                             !   setting to 0, we hope to induce a "divide
1819                                             !   by zero" error if they are
1820                else
1821                   mapfac_arr_x(i,j) = 1.0 / cos(xlat_arr(i,j)*rad_per_deg) 
1822                end if
1823                mapfac_arr_y(i,j) = 1.0
1824             end do
1825          end do
1826     
1827       ! Rotated global cylindrical projection
1828       else if (iproj_type == PROJ_CASSINI) then
1829      
1830          if (abs(pole_lat) == 90.) then
1831             do i=start_mem_i, end_mem_i
1832                do j=start_mem_j, end_mem_j
1833                   if (abs(xlat_arr(i,j)) >= 90.0) then
1834                      mapfac_arr_x(i,j) = 0.    ! MSF actually becomes infinite at poles, but 
1835                                                !   the values should never be used there; by
1836                                                !   setting to 0, we hope to induce a "divide
1837                                                !   by zero" error if they are
1838                   else
1839                      mapfac_arr_x(i,j) = 1.0 / cos(xlat_arr(i,j)*rad_per_deg) 
1840                   end if
1841                   mapfac_arr_y(i,j) = 1.0
1842                end do
1843             end do
1844          else
1845             do i=start_mem_i, end_mem_i
1846                do j=start_mem_j, end_mem_j
1847                   call rotate_coords(xlat_arr(i,j),xlon_arr(i,j), &
1848                                      comp_lat, comp_lon, &
1849                                      pole_lat, pole_lon, stand_lon, &
1850                                      -1)
1851                   if (abs(comp_lat) >= 90.0) then
1852                      mapfac_arr_x(i,j) = 0.    ! MSF actually becomes infinite at poles, but 
1853                                                !   the values should never be used there; by
1854                                                !   setting to 0, we hope to induce a "divide
1855                                                !   by zero" error if they are
1856                   else
1857                      mapfac_arr_x(i,j) = 1.0 / cos(comp_lat*rad_per_deg) 
1858                   end if
1859                   mapfac_arr_y(i,j) = 1.0
1860                end do
1861             end do
1862          end if
1863     
1864       else if (iproj_type == PROJ_ROTLL) then
1865     
1866          do i=start_mem_i, end_mem_i
1867             do j=start_mem_j, end_mem_j
1868                mapfac_arr_x(i,j) = 1.0
1869                mapfac_arr_y(i,j) = 1.0
1870             end do
1871          end do
1872     
1873       end if
1874    
1875    end subroutine get_map_factor
1876    
1877    
1878    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1879    ! Name: get_coriolis_parameters
1880    !
1881    ! Purpose: To calculate the Coriolis parameters f and e for every gridpoint in
1882    !   the tile of the model domain 
1883    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1884    subroutine get_coriolis_parameters(xlat_arr, f, e, &
1885                                       start_mem_i, start_mem_j, end_mem_i, end_mem_j)
1886      
1887       use constants_module
1888     
1889       implicit none
1890     
1891       ! Arguments
1892       integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
1893       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr
1894       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: f, e
1895     
1896       ! Local variables
1897       integer :: i, j
1898     
1899       do i=start_mem_i, end_mem_i
1900          do j=start_mem_j, end_mem_j
1901      
1902             f(i,j) = 2.0*OMEGA_E*sin(rad_per_deg*xlat_arr(i,j))
1903             e(i,j) = 2.0*OMEGA_E*cos(rad_per_deg*xlat_arr(i,j))
1904      
1905          end do
1906       end do
1907       
1908    end subroutine get_coriolis_parameters
1909    
1910    
1911    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1912    ! Name: get_rotang
1913    !
1914    ! Purpose: To calculate the sine and cosine of rotation angle. 
1915    !
1916    ! NOTES: The formulas used in this routine come from those in the 
1917    !   vecrot_rotlat() routine of the original WRF SI.
1918    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1919    subroutine get_rotang(xlat_arr, xlon_arr, cosa, sina, &
1920                          start_mem_i, start_mem_j, end_mem_i, end_mem_j)
1921    
1922       use constants_module
1923       use gridinfo_module
1924     
1925       implicit none
1926     
1927       ! Arguments
1928       integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
1929       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr, xlon_arr
1930       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: cosa, sina
1931     
1932       ! Local variables
1933       integer :: i, j
1934       real :: alpha, d_lon
1936       do i=start_mem_i, end_mem_i
1937          do j=start_mem_j+1, end_mem_j-1
1938             d_lon = xlon_arr(i,j+1)-xlon_arr(i,j-1)
1939             if (d_lon > 180.) then
1940                d_lon = d_lon - 360.
1941             else if (d_lon < -180.) then
1942                d_lon = d_lon + 360.
1943             end if
1945             alpha = atan2(-cos(xlat_arr(i,j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
1946                             ((xlat_arr(i,j+1)-xlat_arr(i,j-1))*RAD_PER_DEG))
1947             sina(i,j) = sin(alpha)
1948             cosa(i,j) = cos(alpha)
1949          end do
1950       end do
1952       do i=start_mem_i, end_mem_i
1953          d_lon = xlon_arr(i,start_mem_j+1)-xlon_arr(i,start_mem_j)
1954          if (d_lon > 180.) then
1955             d_lon = d_lon - 360.
1956          else if (d_lon < -180.) then
1957             d_lon = d_lon + 360.
1958          end if
1960          alpha = atan2(-cos(xlat_arr(i,start_mem_j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
1961                        ((xlat_arr(i,start_mem_j+1)-xlat_arr(i,start_mem_j))*RAD_PER_DEG))
1962          sina(i,start_mem_j) = sin(alpha)
1963          cosa(i,start_mem_j) = cos(alpha)
1964       end do
1966       do i=start_mem_i, end_mem_i
1967          d_lon = xlon_arr(i,end_mem_j)-xlon_arr(i,end_mem_j-1)
1968          if (d_lon > 180.) then
1969             d_lon = d_lon - 360.
1970          else if (d_lon < -180.) then
1971             d_lon = d_lon + 360.
1972          end if
1974          alpha = atan2(-cos(xlat_arr(i,end_mem_j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
1975                        ((xlat_arr(i,end_mem_j)-xlat_arr(i,end_mem_j-1))*RAD_PER_DEG))
1976          sina(i,end_mem_j) = sin(alpha)
1977          cosa(i,end_mem_j) = cos(alpha)
1978       end do
1979     
1980    end subroutine get_rotang
1981    
1982    
1983    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1984    ! Name: process_neighbor
1985    !
1986    ! Purpose: This routine, give the x/y location of a point, determines whether
1987    !   the point has already been processed, and if not, which processing queue
1988    !   the point should be placed in.
1989    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1990    subroutine process_neighbor(ix, iy, bit_domain, point_queue, tile_queue, &
1991                                xlat_array, xlon_array, &
1992                                start_i, end_i, start_j, end_j, ilevel)
1993    
1994       use bitarray_module
1995       use misc_definitions_module
1996       use proc_point_module
1997       use queue_module
1998     
1999       implicit none
2000     
2001       ! Arguments
2002       integer, intent(in) :: ix, iy, start_i, end_i, start_j, end_j, ilevel
2003       real, dimension(start_i:end_i, start_j:end_j), intent(in) :: xlat_array, xlon_array
2004       type (bitarray), intent(inout) :: bit_domain
2005       type (queue), intent(inout) :: point_queue, tile_queue
2006     
2007       ! Local variables
2008       type (q_data) :: process_pt
2009       logical :: is_in_tile
2010     
2011       ! If the point has already been visited, no need to do anything more.
2012       if (.not. bitarray_test(bit_domain, ix-start_i+1, iy-start_j+1)) then 
2013     
2014          ! Create a queue item for the current point
2015          process_pt%lat = xlat_array(ix,iy)            
2016          process_pt%lon = xlon_array(ix,iy)            
2017          process_pt%x = ix
2018          process_pt%y = iy
2019      
2020          is_in_tile = is_point_in_tile(process_pt%lat, process_pt%lon, ilevel)
2021      
2022          ! If the point is in the current tile, add it to the list of points
2023          !   to be processed in the inner loop
2024          if (is_in_tile) then
2025             call q_insert(point_queue, process_pt)
2026             call bitarray_set(bit_domain, ix-start_i+1, iy-start_j+1) 
2027      
2028          ! Otherwise, we will process this point later. Add it to the list for 
2029          !   the outer loop.
2030          else
2031             call q_insert(tile_queue, process_pt)
2032          end if
2033     
2034       end if
2035    
2036    end subroutine process_neighbor
2037    
2038    
2039    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2040    ! Name: calc_dfdy
2041    !
2042    ! Purpose: This routine calculates df/dy for the field in src_arr, and places
2043    !   the result in dst_array.
2044    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2045    subroutine calc_dfdy(src_arr, dst_arr, start_mem_i, start_mem_j, start_mem_k, &
2046                         end_mem_i, end_mem_j, end_mem_k, mapfac)
2047    
2048       ! Modules
2049       use gridinfo_module
2050     
2051       implicit none
2052     
2053       ! Arguments
2054       integer, intent(in) :: start_mem_i, start_mem_j, start_mem_k, end_mem_i, end_mem_j, end_mem_k
2055       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j,start_mem_k:end_mem_k), intent(in) :: src_arr
2056       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j,start_mem_k:end_mem_k), intent(out) :: dst_arr
2057       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in), optional :: mapfac
2058     
2059       ! Local variables
2060       integer :: i, j, k
2061     
2062       if (present(mapfac)) then
2063          do k=start_mem_k,end_mem_k
2064             do i=start_mem_i, end_mem_i
2065                do j=start_mem_j+1, end_mem_j-1
2066                   dst_arr(i,j,k) = (src_arr(i,j+1,k) - src_arr(i,j-1,k))/(2.*dykm*mapfac(i,j))
2067                end do
2068             end do
2069      
2070             do i=start_mem_i, end_mem_i
2071                dst_arr(i,start_mem_j,k) = (src_arr(i,start_mem_j+1,k) - src_arr(i,start_mem_j,k))/(dykm*mapfac(i,j))
2072             end do
2073      
2074             do i=start_mem_i, end_mem_i
2075                dst_arr(i,end_mem_j,k) = (src_arr(i,end_mem_j,k) - src_arr(i,end_mem_j-1,k))/(dykm*mapfac(i,j))
2076             end do
2077          end do
2078       else
2079          do k=start_mem_k,end_mem_k
2080             do i=start_mem_i, end_mem_i
2081                do j=start_mem_j+1, end_mem_j-1
2082                   dst_arr(i,j,k) = (src_arr(i,j+1,k) - src_arr(i,j-1,k))/(2.*dykm)
2083                end do
2084             end do
2085      
2086             do i=start_mem_i, end_mem_i
2087                dst_arr(i,start_mem_j,k) = (src_arr(i,start_mem_j+1,k) - src_arr(i,start_mem_j,k))/(dykm)
2088             end do
2089      
2090             do i=start_mem_i, end_mem_i
2091                dst_arr(i,end_mem_j,k) = (src_arr(i,end_mem_j,k) - src_arr(i,end_mem_j-1,k))/(dykm)
2092             end do
2093          end do
2094       end if
2095    
2096    end subroutine calc_dfdy
2097    
2098    
2099    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2100    ! Name: calc_dfdx
2101    !
2102    ! Purpose: This routine calculates df/dx for the field in src_arr, and places
2103    !   the result in dst_array.
2104    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2105    subroutine calc_dfdx(src_arr, dst_arr, start_mem_i, start_mem_j, &
2106                         start_mem_k, end_mem_i, end_mem_j, end_mem_k, mapfac)
2107    
2108       ! Modules
2109       use gridinfo_module
2110     
2111       implicit none
2112     
2113       ! Arguments
2114       integer, intent(in) :: start_mem_i, start_mem_j, start_mem_k, end_mem_i, end_mem_j, end_mem_k
2115       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j, start_mem_k:end_mem_k), intent(in) :: src_arr
2116       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j, start_mem_k:end_mem_k), intent(out) :: dst_arr
2117       real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in), optional :: mapfac
2118     
2119       ! Local variables
2120       integer :: i, j, k
2121     
2122       if (present(mapfac)) then
2123          do k=start_mem_k, end_mem_k
2124             do i=start_mem_i+1, end_mem_i-1
2125                do j=start_mem_j, end_mem_j
2126                   dst_arr(i,j,k) = (src_arr(i+1,j,k) - src_arr(i-1,j,k))/(2.*dxkm*mapfac(i,j))
2127                end do
2128             end do
2129      
2130             do j=start_mem_j, end_mem_j
2131                dst_arr(start_mem_i,j,k) = (src_arr(start_mem_i+1,j,k) - src_arr(start_mem_i,j,k))/(dxkm*mapfac(i,j))
2132             end do
2133      
2134             do j=start_mem_j, end_mem_j
2135                dst_arr(end_mem_i,j,k) = (src_arr(end_mem_i,j,k) - src_arr(end_mem_i-1,j,k))/(dxkm*mapfac(i,j))
2136             end do
2137          end do
2138       else
2139          do k=start_mem_k, end_mem_k
2140             do i=start_mem_i+1, end_mem_i-1
2141                do j=start_mem_j, end_mem_j
2142                   dst_arr(i,j,k) = (src_arr(i+1,j,k) - src_arr(i-1,j,k))/(2.*dxkm)
2143                end do
2144             end do
2145      
2146             do j=start_mem_j, end_mem_j
2147                dst_arr(start_mem_i,j,k) = (src_arr(start_mem_i+1,j,k) - src_arr(start_mem_i,j,k))/(dxkm)
2148             end do
2149      
2150             do j=start_mem_j, end_mem_j
2151                dst_arr(end_mem_i,j,k) = (src_arr(end_mem_i,j,k) - src_arr(end_mem_i-1,j,k))/(dxkm)
2152             end do
2153          end do
2154       end if
2155    
2156    end subroutine calc_dfdx
2158 end module process_tile_module