PageRenderTime 52ms CodeModel.GetById 17ms RepoModel.GetById 0ms app.codeStats 1ms

/WPS/metgrid/src/process_domain_module.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2839 lines | 1753 code | 456 blank | 630 comment | 397 complexity | 2d29f8bcd4bc0d9d1f57c097777cb2a6 MD5 | raw file
Possible License(s): AGPL-1.0

Large files files are truncated, but you can click here to view the full file

    

  1. module process_domain_module
    2. 

  2. 

  3.    contains
    4. 

  4. 

  5.    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    6. 

  6.    ! Name: process_domain
    7. 

  7.    !
    8. 

  8.    ! Purpose:
    9. 

  9.    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    10. 

  10.    subroutine process_domain(n, extra_row, extra_col)
    11. 

  11.    
    12. 

  12.       use date_pack
    13. 

  13.       use gridinfo_module
    14. 

  14.       use interp_option_module
    15. 

  15.       use misc_definitions_module
    16. 

  16.       use module_debug
    17. 

  17.       use storage_module
    18. 

  18.    
    19. 

  19.       implicit none
    20. 

  20.    
    21. 

  21.       ! Arguments
    22. 

  22.       integer, intent(in) :: n
    23. 

  23.       logical, intent(in) :: extra_row, extra_col
    24. 

  24.    
    25. 

  25.       ! Local variables
    26. 

  26.       integer :: i, t, dyn_opt, &
    27. 

  27.                  we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
    28. 

  28.                  we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
    29. 

  29.                  sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
    30. 

  30.                  we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
    31. 

  31.                  sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
    32. 

  32.                  idiff, n_times, &
    33. 

  33.                  west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
    34. 

  34.                  is_water, is_lake, is_ice, is_urban, i_soilwater, &
    35. 

  35.                  grid_id, parent_id, i_parent_start, j_parent_start, &
    36. 

  36.                  i_parent_end, j_parent_end, parent_grid_ratio, sub_x, sub_y, num_land_cat, process_bdy_width
    37. 

  37.       real :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
    38. 

  38.               dom_dx, dom_dy, pole_lat, pole_lon
    39. 

  39.       real, dimension(16) :: corner_lats, corner_lons
    40. 

  40.       real, pointer, dimension(:,:) :: landmask
    41. 

  41.       real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
    42. 

  42.       logical, allocatable, dimension(:) :: got_this_field, got_const_field
    43. 

  43.       character (len=19) :: valid_date, temp_date
    44. 

  44.       character (len=128) :: title, mminlu
    45. 

  45.       character (len=128), allocatable, dimension(:) :: output_flags, td_output_flags
    46. 

  46. 

  47.       ! Compute number of times that we will process
    48. 

  48.       call geth_idts(end_date(n), start_date(n), idiff)
    49. 

  49.       call mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist for domain %i.', i1=n)
    50. 

  50.    
    51. 

  51.       n_times = idiff / interval_seconds
    52. 

  52.    
    53. 

  53.       ! Check that the interval evenly divides the range of times to process
    54. 

  54.       call mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
    55. 

  55.                    'In namelist, interval_seconds does not evenly divide '// &
    56. 

  56.                    '(end_date - start_date) for domain %i. Only %i time periods '// &
    57. 

  57.                    'will be processed.', i1=n, i2=n_times)
    58. 

  58.    
    59. 

  59.       ! Initialize the storage module
    60. 

  60.       call mprintf(.true.,LOGFILE,'Initializing storage module')
    61. 

  61.       call storage_init()
    62. 

  62.    
    63. 

  63.       ! 
    64. 

  64.       ! Do time-independent processing
    65. 

  65.       ! 
    66. 

  66.       call get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
    67. 

  67.                     we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
    68. 

  68.                     we_patch_s,      we_patch_e, &
    69. 

  69.                     we_patch_stag_s, we_patch_stag_e, &
    70. 

  70.                     sn_patch_s,      sn_patch_e, &
    71. 

  71.                     sn_patch_stag_s, sn_patch_stag_e, &
    72. 

  72.                     we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
    73. 

  73.                     sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
    74. 

  74.                     mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
    75. 

  75.                     grid_id, parent_id, &
    76. 

  76.                     i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
    77. 

  77.                     parent_grid_ratio, sub_x, sub_y, &
    78. 

  78.                     cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
    79. 

  79.                     pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
    80. 

  80.                     xlat_v, xlon_v, corner_lats, corner_lons, title)
    81. 

  81. 

  82. 

  83.       allocate(output_flags(num_entries))
    84. 

  84.       allocate(got_const_field(num_entries))
    85. 

  85. 

  86.       do i=1,num_entries
    87. 

  87.          output_flags(i)    = ' '
    88. 

  88.          got_const_field(i) = .false.
    89. 

  89.       end do
    90. 

  90.    
    91. 

  91.       ! This call is to process the constant met fields (SST or SEAICE, for example)
    92. 

  92.       ! That we process constant fields is indicated by the first argument
    93. 

  93.       call process_single_met_time(.true., temp_date, n, extra_row, extra_col, xlat, xlon, &
    94. 

  94.                           xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
    95. 

  95.                           title, dyn_opt, &
    96. 

  96.                           west_east_dim, south_north_dim, &
    97. 

  97.                           we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
    98. 

  98.                           we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
    99. 

  99.                           sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
    100. 

  100.                           we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
    101. 

  101.                           sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
    102. 

  102.                           got_const_field, &
    103. 

  103.                           map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
    104. 

  104.                           grid_id, parent_id, i_parent_start, &
    105. 

  105.                           j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
    106. 

  106.                           cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
    107. 

  107.                           pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
    108. 

  108.                           corner_lats, corner_lons, output_flags, 0)
    109. 

  109. 

  110.       !
    111. 

  111.       ! Begin time-dependent processing
    112. 

  112.       !
    113. 

  113. 

  114.       allocate(td_output_flags(num_entries))
    115. 

  115.       allocate(got_this_field (num_entries))
    116. 

  116.    
    117. 

  117.       ! Loop over all times to be processed for this domain
    118. 

  118.       do t=0,n_times
    119. 

  119.    
    120. 

  120.          call geth_newdate(valid_date, trim(start_date(n)), t*interval_seconds)
    121. 

  121.          temp_date = ' '
    122. 

  122. 

  123.          if (mod(interval_seconds,3600) == 0) then
    124. 

  124.             write(temp_date,'(a13)') valid_date(1:10)//'_'//valid_date(12:13)
    125. 

  125.          else if (mod(interval_seconds,60) == 0) then
    126. 

  126.             write(temp_date,'(a16)') valid_date(1:10)//'_'//valid_date(12:16)
    127. 

  127.          else
    128. 

  128.             write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
    129. 

  129.          end if
    130. 

  130.    
    131. 

  131.          call mprintf(.true.,STDOUT, ' Processing %s', s1=trim(temp_date))
    132. 

  132.          call mprintf(.true.,LOGFILE, 'Preparing to process output time %s', s1=temp_date)
    133. 

  133.    
    134. 

  134.          do i=1,num_entries
    135. 

  135.             td_output_flags(i) = output_flags(i)
    136. 

  136.             got_this_field(i)  = got_const_field(i)
    137. 

  137.          end do
    138. 

  138. 

  139.          if (t > 0) then
    140. 

  140.             process_bdy_width = process_only_bdy
    141. 

  141.          else
    142. 

  142.             process_bdy_width = 0
    143. 

  143.          end if
    144. 

  144.    
    145. 

  145.          call process_single_met_time(.false., temp_date, n, extra_row, extra_col, xlat, xlon, &
    146. 

  146.                              xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
    147. 

  147.                              title, dyn_opt, &
    148. 

  148.                              west_east_dim, south_north_dim, &
    149. 

  149.                              we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
    150. 

  150.                              we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
    151. 

  151.                              sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
    152. 

  152.                              we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
    153. 

  153.                              sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
    154. 

  154.                              got_this_field, &
    155. 

  155.                              map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
    156. 

  156.                              grid_id, parent_id, i_parent_start, &
    157. 

  157.                              j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
    158. 

  158.                              cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
    159. 

  159.                              pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
    160. 

  160.                              corner_lats, corner_lons, td_output_flags, process_bdy_width)
    161. 

  161.    
    162. 

  162.       end do  ! Loop over n_times
    163. 

  163. 

  164. 

  165.       deallocate(td_output_flags)
    166. 

  166.       deallocate(got_this_field)
    167. 

  167. 

  168.       deallocate(output_flags)
    169. 

  169.       deallocate(got_const_field)
    170. 

  170.    
    171. 

  171.       call storage_delete_all()
    172. 

  172.    
    173. 

  173.    end subroutine process_domain
    174. 

  174. 

  175. 

  176.    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    177. 

  177.    ! Name: get_static_fields
    178. 

  178.    !
    179. 

  179.    ! Purpose:
    180. 

  180.    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    181. 

  181.    subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
    182. 

  182.                     map_proj,                                                               &
    183. 

  183.                     we_dom_s,   we_dom_e,   sn_dom_s,        sn_dom_e,                      &
    184. 

  184.                     we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e,               &
    185. 

  185.                     sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e,               &
    186. 

  186.                     we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e,                       &
    187. 

  187.                     sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e,                       &
    188. 

  188.                     mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
    189. 

  189.                     grid_id, parent_id,                                                     &
    190. 

  190.                     i_parent_start, j_parent_start, i_parent_end, j_parent_end,             &
    191. 

  191.                     parent_grid_ratio, sub_x, sub_y,                                        &
    192. 

  192.                     cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2,          &
    193. 

  193.                     pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
    194. 

  194.                     xlat_v, xlon_v, corner_lats, corner_lons, title)
    195. 

  195. 

  196.       use gridinfo_module
    197. 

  197.       use input_module
    198. 

  198.       use llxy_module
    199. 

  199.       use parallel_module
    200. 

  200.       use storage_module
    201. 

  201.       use module_debug
    202. 

  202. 

  203.       implicit none
    204. 

  204. 

  205.       ! Arguments
    206. 

  206.       integer, intent(in) :: n
    207. 

  207.       integer, intent(inout) :: dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
    208. 

  208.                                 map_proj, &
    209. 

  209.                                 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
    210. 

  210.                                 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
    211. 

  211.                                 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
    212. 

  212.                                 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
    213. 

  213.                                 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
    214. 

  214.                                 is_water, is_lake, is_ice, is_urban, i_soilwater, grid_id, parent_id, &
    215. 

  215.                                 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
    216. 

  216.                                 parent_grid_ratio, sub_x, sub_y, num_land_cat
    217. 

  217.       real, pointer, dimension(:,:) :: landmask
    218. 

  218.       real, intent(inout) :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
    219. 

  219.                              dom_dx, dom_dy, pole_lat, pole_lon
    220. 

  220.       real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
    221. 

  221.       real, dimension(16), intent(out) :: corner_lats, corner_lons
    222. 

  222.       character (len=128), intent(inout) :: title, mminlu
    223. 

  223.     
    224. 

  224.       ! Local variables
    225. 

  225.       integer :: istatus, i, j, k, sp1, ep1, sp2, ep2, sp3, ep3, &
    226. 

  226.                  lh_mult, rh_mult, bh_mult, th_mult
    227. 

  227.       integer :: we_mem_subgrid_s, we_mem_subgrid_e, &
    228. 

  228.                  sn_mem_subgrid_s, sn_mem_subgrid_e
    229. 

  229.       integer :: we_patch_subgrid_s, we_patch_subgrid_e, &
    230. 

  230.                  sn_patch_subgrid_s, sn_patch_subgrid_e
    231. 

  231.       real, pointer, dimension(:,:,:) :: real_array
    232. 

  232.       character (len=3) :: memorder
    233. 

  233.       character (len=128) :: grid_type, datestr, cname, stagger, cunits, cdesc
    234. 

  234.       character (len=128), dimension(3) :: dimnames
    235. 

  235.       type (fg_input) :: field
    236. 

  236.       logical :: is_subgrid
    237. 

  237. 

  238.       ! Initialize the input module to read static input data for this domain
    239. 

  239.       call mprintf(.true.,LOGFILE,'Opening static input file.')
    240. 

  240.       call input_init(n, istatus)
    241. 

  241.       call mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input for domain %i.', i1=n)
    242. 

  242.    
    243. 

  243.       ! Read global attributes from the static data input file 
    244. 

  244.       call mprintf(.true.,LOGFILE,'Reading static global attributes.')
    245. 

  245.       call read_global_attrs(title, datestr, grid_type, dyn_opt, west_east_dim,          &
    246. 

  246.                              south_north_dim, bottom_top_dim,                            &
    247. 

  247.                              we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e,   &
    248. 

  248.                              sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e,   &
    249. 

  249.                              map_proj, mminlu, num_land_cat,                             &
    250. 

  250.                              is_water, is_lake, is_ice, is_urban, i_soilwater,           &
    251. 

  251.                              grid_id, parent_id, i_parent_start,                         &
    252. 

  252.                              j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
    253. 

  253.                              cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1,        &
    254. 

  254.                              truelat2, pole_lat, pole_lon, parent_grid_ratio,            &
    255. 

  255.                              corner_lats, corner_lons, sub_x, sub_y)
    256. 

  256. 

  257.       we_dom_s = 1
    258. 

  258.       sn_dom_s = 1
    259. 

  259.       if (grid_type(1:1) == 'C') then
    260. 

  260.          we_dom_e = west_east_dim   - 1
    261. 

  261.          sn_dom_e = south_north_dim - 1
    262. 

  262.       else if (grid_type(1:1) == 'E') then
    263. 

  263.          we_dom_e = west_east_dim 
    264. 

  264.          sn_dom_e = south_north_dim
    265. 

  265.       end if
    266. 

  266.      
    267. 

  267.       ! Given the full dimensions of this domain, find out the range of indices 
    268. 

  268.       !   that will be worked on by this processor. This information is given by 
    269. 

  269.       !   my_minx, my_miny, my_maxx, my_maxy
    270. 

  270.       call parallel_get_tile_dims(west_east_dim, south_north_dim)
    271. 

  271. 

  272.       ! Must figure out patch dimensions from info in parallel module
    273. 

  273.       if (nprocs > 1 .and. .not. do_tiled_input) then
    274. 

  274. 

  275.          we_patch_s      = my_minx
    276. 

  276.          we_patch_stag_s = my_minx
    277. 

  277.          we_patch_e      = my_maxx - 1
    278. 

  278.          sn_patch_s      = my_miny
    279. 

  279.          sn_patch_stag_s = my_miny
    280. 

  280.          sn_patch_e      = my_maxy - 1
    281. 

  281. 

  282.          if (gridtype == 'C') then
    283. 

  283.             if (my_x /= nproc_x - 1) then
    284. 

  284.                we_patch_e = we_patch_e + 1
    285. 

  285.                we_patch_stag_e = we_patch_e
    286. 

  286.             else
    287. 

  287.                we_patch_stag_e = we_patch_e + 1
    288. 

  288.             end if
    289. 

  289.             if (my_y /= nproc_y - 1) then
    290. 

  290.                sn_patch_e = sn_patch_e + 1
    291. 

  291.                sn_patch_stag_e = sn_patch_e
    292. 

  292.             else
    293. 

  293.                sn_patch_stag_e = sn_patch_e + 1
    294. 

  294.             end if
    295. 

  295.          else if (gridtype == 'E') then
    296. 

  296.             we_patch_e = we_patch_e + 1
    297. 

  297.             sn_patch_e = sn_patch_e + 1
    298. 

  298.             we_patch_stag_e = we_patch_e
    299. 

  299.             sn_patch_stag_e = sn_patch_e
    300. 

  300.          end if
    301. 

  301. 

  302.       end if
    303. 

  303. 

  304.       ! Compute multipliers for halo width; these must be 0/1
    305. 

  305.       if (my_x /= 0) then
    306. 

  306.         lh_mult = 1
    307. 

  307.       else
    308. 

  308.         lh_mult = 0
    309. 

  309.       end if
    310. 

  310.       if (my_x /= (nproc_x-1)) then
    311. 

  311.         rh_mult = 1
    312. 

  312.       else
    313. 

  313.         rh_mult = 0
    314. 

  314.       end if
    315. 

  315.       if (my_y /= 0) then
    316. 

  316.         bh_mult = 1
    317. 

  317.       else
    318. 

  318.         bh_mult = 0
    319. 

  319.       end if
    320. 

  320.       if (my_y /= (nproc_y-1)) then
    321. 

  321.         th_mult = 1
    322. 

  322.       else
    323. 

  323.         th_mult = 0
    324. 

  324.       end if
    325. 

  325. 

  326.       we_mem_s = we_patch_s - HALO_WIDTH*lh_mult
    327. 

  327.       we_mem_e = we_patch_e + HALO_WIDTH*rh_mult
    328. 

  328.       sn_mem_s = sn_patch_s - HALO_WIDTH*bh_mult
    329. 

  329.       sn_mem_e = sn_patch_e + HALO_WIDTH*th_mult
    330. 

  330.       we_mem_stag_s = we_patch_stag_s - HALO_WIDTH*lh_mult
    331. 

  331.       we_mem_stag_e = we_patch_stag_e + HALO_WIDTH*rh_mult
    332. 

  332.       sn_mem_stag_s = sn_patch_stag_s - HALO_WIDTH*bh_mult
    333. 

  333.       sn_mem_stag_e = sn_patch_stag_e + HALO_WIDTH*th_mult
    334. 

  334. 

  335.       ! Initialize a proj_info type for the destination grid projection. This will
    336. 

  336.       !   primarily be used for rotating Earth-relative winds to grid-relative winds
    337. 

  337.       call set_domain_projection(map_proj, stand_lon, truelat1, truelat2, &
    338. 

  338.                                  dom_dx, dom_dy, dom_dx, dom_dy, west_east_dim, &
    339. 

  339.                                  south_north_dim, real(west_east_dim)/2., &
    340. 

  340.                                  real(south_north_dim)/2.,cen_lat, cen_lon, &
    341. 

  341.                                  cen_lat, cen_lon)
    342. 

  342.    
    343. 

  343.       ! Read static fields using the input module; we know that there are no more
    344. 

  344.       !   fields to be read when read_next_field() returns a non-zero status.
    345. 

  345.       istatus = 0
    346. 

  346.       do while (istatus == 0)  
    347. 

  347.         call read_next_field(sp1, ep1, sp2, ep2, sp3, ep3, cname, cunits, cdesc, &
    348. 

  348.                              memorder, stagger, dimnames, is_subgrid, &
    349. 

  349.                              real_array, istatus)
    350. 

  350.         if (istatus == 0) then
    351. 

  351. 

  352.           call mprintf(.true.,LOGFILE, 'Read in static field %s.',s1=cname)
    353. 

  353.    
    354. 

  354.           ! We will also keep copies in core of the lat/lon arrays, for use in 
    355. 

  355.           !    interpolation of the met fields to the model grid.
    356. 

  356.           ! For now, we assume that the lat/lon arrays will have known field names
    357. 

  357.           if (index(cname, 'XLAT_M') /= 0 .and. &
    358. 

  358.               len_trim(cname) == len_trim('XLAT_M')) then
    359. 

  359.              allocate(xlat(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
    360. 

  360.              xlat(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
    361. 

  361.              call exchange_halo_r(xlat, & 
    362. 

  362.                                   we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
    363. 

  363.                                   we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
    364. 

  364. 

  365.           else if (index(cname, 'XLONG_M') /= 0 .and. &
    366. 

  366.                    len_trim(cname) == len_trim('XLONG_M')) then
    367. 

  367.              allocate(xlon(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
    368. 

  368.              xlon(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
    369. 

  369.              call exchange_halo_r(xlon, & 
    370. 

  370.                                   we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
    371. 

  371.                                   we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
    372. 

  372. 

  373.           else if (index(cname, 'XLAT_U') /= 0 .and. &
    374. 

  374.                    len_trim(cname) == len_trim('XLAT_U')) then
    375. 

  375.              allocate(xlat_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
    376. 

  376.              xlat_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
    377. 

  377.              call exchange_halo_r(xlat_u, & 
    378. 

  378.                                   we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
    379. 

  379.                                   we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
    380. 

  380. 

  381.           else if (index(cname, 'XLONG_U') /= 0 .and. &
    382. 

  382.                    len_trim(cname) == len_trim('XLONG_U')) then
    383. 

  383.              allocate(xlon_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
    384. 

  384.              xlon_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
    385. 

  385.              call exchange_halo_r(xlon_u, & 
    386. 

  386.                                   we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
    387. 

  387.                                   we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
    388. 

  388. 

  389.           else if (index(cname, 'XLAT_V') /= 0 .and. &
    390. 

  390.                    len_trim(cname) == len_trim('XLAT_V')) then
    391. 

  391.              allocate(xlat_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
    392. 

  392.              xlat_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
    393. 

  393.              call exchange_halo_r(xlat_v, & 
    394. 

  394.                                   we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
    395. 

  395.                                   we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
    396. 

  396. 

  397.           else if (index(cname, 'XLONG_V') /= 0 .and. &
    398. 

  398.                    len_trim(cname) == len_trim('XLONG_V')) then
    399. 

  399.              allocate(xlon_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
    400. 

  400.              xlon_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
    401. 

  401.              call exchange_halo_r(xlon_v, & 
    402. 

  402.                                   we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
    403. 

  403.                                   we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
    404. 

  404. 

  405.           else if (index(cname, 'LANDMASK') /= 0 .and. &
    406. 

  406.                    len_trim(cname) == len_trim('LANDMASK')) then
    407. 

  407.              allocate(landmask(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
    408. 

  408.              landmask(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
    409. 

  409.              call exchange_halo_r(landmask, & 
    410. 

  410.                                   we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
    411. 

  411.                                   we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
    412. 

  412. 

  413.           end if
    414. 

  414. 

  415.           we_mem_subgrid_s   = (we_mem_s                 + HALO_WIDTH*lh_mult - 1)*sub_x - HALO_WIDTH*lh_mult + 1
    416. 

  416.           we_mem_subgrid_e   = (we_mem_e   + (1-rh_mult) - HALO_WIDTH*rh_mult    )*sub_x + HALO_WIDTH*rh_mult
    417. 

  417.           we_patch_subgrid_s = (we_patch_s                                    - 1)*sub_x                      + 1
    418. 

  418.           we_patch_subgrid_e = (we_patch_e + (1-rh_mult)                         )*sub_x
    419. 

  419.           sn_mem_subgrid_s   = (sn_mem_s                 + HALO_WIDTH*bh_mult - 1)*sub_y - HALO_WIDTH*bh_mult + 1
    420. 

  420.           sn_mem_subgrid_e   = (sn_mem_e   + (1-th_mult) - HALO_WIDTH*th_mult    )*sub_y + HALO_WIDTH*th_mult
    421. 

  421.           sn_patch_subgrid_s = (sn_patch_s                                    - 1)*sub_y                      + 1
    422. 

  422.           sn_patch_subgrid_e = (sn_patch_e + (1-th_mult)                         )*sub_y
    423. 

  423.     
    424. 

  424.           ! Having read in a field, we write each level individually to the
    425. 

  425.           !   storage module; levels will be reassembled later on when they
    426. 

  426.           !   are written.
    427. 

  427.           do k=sp3,ep3
    428. 

  428.              field%header%sr_x=sub_x
    429. 

  429.              field%header%sr_y=sub_y
    430. 

  430.              field%header%is_subgrid=is_subgrid
    431. 

  431.              field%header%version = 1
    432. 

  432.              field%header%date = start_date(n) 
    433. 

  433.              field%header%time_dependent = .false.
    434. 

  434.              field%header%mask_field = .false.
    435. 

  435.              field%header%forecast_hour = 0.0
    436. 

  436.              field%header%fg_source = 'geogrid_model'
    437. 

  437.              field%header%field = cname
    438. 

  438.              field%header%units = cunits
    439. 

  439.              field%header%description = cdesc
    440. 

  440.              field%header%vertical_coord = dimnames(3) 
    441. 

  441.              field%header%vertical_level = k
    442. 

  442.              field%header%array_order = memorder
    443. 

  443.              field%header%is_wind_grid_rel = .true.
    444. 

  444.              field%header%array_has_missing_values = .false.
    445. 

  445.              if (gridtype == 'C') then
    446. 

  446.                 if (is_subgrid) then
    447. 

  447.                    field%map%stagger = M
    448. 

  448.                    field%header%dim1(1) = we_mem_subgrid_s
    449. 

  449.                    field%header%dim1(2) = we_mem_subgrid_e
    450. 

  450.                    field%header%dim2(1) = sn_mem_subgrid_s
    451. 

  451.                    field%header%dim2(2) = sn_mem_subgrid_e
    452. 

  452.                 else if (trim(stagger) == 'M') then
    453. 

  453.                    field%map%stagger = M
    454. 

  454.                    field%header%dim1(1) = we_mem_s
    455. 

  455.                    field%header%dim1(2) = we_mem_e
    456. 

  456.                    field%header%dim2(1) = sn_mem_s
    457. 

  457.                    field%header%dim2(2) = sn_mem_e
    458. 

  458.                 else if (trim(stagger) == 'U') then
    459. 

  459.                    field%map%stagger = U
    460. 

  460.                    field%header%dim1(1) = we_mem_stag_s
    461. 

  461.                    field%header%dim1(2) = we_mem_stag_e
    462. 

  462.                    field%header%dim2(1) = sn_mem_s
    463. 

  463.                    field%header%dim2(2) = sn_mem_e
    464. 

  464.                 else if (trim(stagger) == 'V') then
    465. 

  465.                    field%map%stagger = V
    466. 

  466.                    field%header%dim1(1) = we_mem_s
    467. 

  467.                    field%header%dim1(2) = we_mem_e
    468. 

  468.                    field%header%dim2(1) = sn_mem_stag_s
    469. 

  469.                    field%header%dim2(2) = sn_mem_stag_e
    470. 

  470.                 end if
    471. 

  471.              else if (gridtype == 'E') then
    472. 

  472.                 if (trim(stagger) == 'M') then
    473. 

  473.                    field%map%stagger = HH
    474. 

  474.                 else if (trim(stagger) == 'V') then
    475. 

  475.                    field%map%stagger = VV
    476. 

  476.                 end if
    477. 

  477.                 field%header%dim1(1) = we_mem_s
    478. 

  478.                 field%header%dim1(2) = we_mem_e
    479. 

  479.                 field%header%dim2(1) = sn_mem_s
    480. 

  480.                 field%header%dim2(2) = sn_mem_e
    481. 

  481.              end if
    482. 

  482. 

  483.              allocate(field%valid_mask)
    484. 

  484. 

  485.              if (is_subgrid) then
    486. 

  486.                 allocate(field%r_arr(we_mem_subgrid_s:we_mem_subgrid_e,&
    487. 

  487.                                      sn_mem_subgrid_s:sn_mem_subgrid_e))
    488. 

  488.                 field%r_arr(we_patch_subgrid_s:we_patch_subgrid_e,sn_patch_subgrid_s:sn_patch_subgrid_e) = &
    489. 

  489.                            real_array(sp1:ep1,sp2:ep2,k)
    490. 

  490.                 call exchange_halo_r(field%r_arr, &
    491. 

  491.                            we_mem_subgrid_s, we_mem_subgrid_e, sn_mem_subgrid_s, sn_mem_subgrid_e, 1, 1, &
    492. 

  492.                            we_patch_subgrid_s, we_patch_subgrid_e, sn_patch_subgrid_s, sn_patch_subgrid_e, 1, 1)
    493. 

  493.                 call bitarray_create(field%valid_mask, &
    494. 

  494.                                      (we_mem_subgrid_e-we_mem_subgrid_s)+1, &
    495. 

  495.                                      (sn_mem_subgrid_e-sn_mem_subgrid_s)+1)
    496. 

  496.                 do j=1,(sn_mem_subgrid_e-sn_mem_subgrid_s)+1
    497. 

  497.                    do i=1,(we_mem_subgrid_e-we_mem_subgrid_s)+1
    498. 

  498.                       call bitarray_set(field%valid_mask, i, j)     
    499. 

  499.                    end do
    500. 

  500.                 end do
    501. 

  501. 

  502.              else if (field%map%stagger == M  .or. & 
    503. 

  503.                  field%map%stagger == HH .or. &
    504. 

  504.                  field%map%stagger == VV) then
    505. 

  505.                 allocate(field%r_arr(we_mem_s:we_mem_e,&
    506. 

  506.                                      sn_mem_s:sn_mem_e))
    507. 

  507.                 field%r_arr(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
    508. 

  508.                 call exchange_halo_r(field%r_arr, &
    509. 

  509.                            we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
    510. 

  510.                            we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
    511. 

  511.                 call bitarray_create(field%valid_mask, &
    512. 

  512.                                      (we_mem_e-we_mem_s)+1, &
    513. 

  513.                                      (sn_mem_e-sn_mem_s)+1)
    514. 

  514.                 do j=1,(sn_mem_e-sn_mem_s)+1
    515. 

  515.                    do i=1,(we_mem_e-we_mem_s)+1
    516. 

  516.                       call bitarray_set(field%valid_mask, i, j)     
    517. 

  517.                    end do
    518. 

  518.                 end do
    519. 

  519.              else if (field%map%stagger == U) then
    520. 

  520.                 allocate(field%r_arr(we_mem_stag_s:we_mem_stag_e,&
    521. 

  521.                                      sn_mem_s:sn_mem_e))
    522. 

  522.                 field%r_arr(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
    523. 

  523.                 call exchange_halo_r(field%r_arr, &
    524. 

  524.                            we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
    525. 

  525.                            we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
    526. 

  526.                 call bitarray_create(field%valid_mask, &
    527. 

  527.                                      (we_mem_stag_e-we_mem_stag_s)+1, &
    528. 

  528.                                      (sn_mem_e-sn_mem_s)+1)
    529. 

  529.                 do j=1,(sn_mem_e-sn_mem_s)+1
    530. 

  530.                    do i=1,(we_mem_stag_e-we_mem_stag_s)+1
    531. 

  531.                       call bitarray_set(field%valid_mask, i, j)     
    532. 

  532.                    end do
    533. 

  533.                 end do
    534. 

  534.              else if (field%map%stagger == V) then
    535. 

  535.                 allocate(field%r_arr(we_mem_s:we_mem_e,&
    536. 

  536.                                      sn_mem_stag_s:sn_mem_stag_e))
    537. 

  537.                 field%r_arr(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(sp1:ep1,sp2:ep2,k)
    538. 

  538.                 call exchange_halo_r(field%r_arr, &
    539. 

  539.                            we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
    540. 

  540.                            we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
    541. 

  541.                 call bitarray_create(field%valid_mask, &
    542. 

  542.                                      (we_mem_e-we_mem_s)+1, &
    543. 

  543.                                      (sn_mem_stag_e-sn_mem_stag_s)+1)
    544. 

  544.                 do j=1,(sn_mem_stag_e-sn_mem_stag_s)+1
    545. 

  545.                    do i=1,(we_mem_e-we_mem_s)+1
    546. 

  546.                       call bitarray_set(field%valid_mask, i, j)     
    547. 

  547.                    end do
    548. 

  548.                 end do
    549. 

  549.              end if
    550. 

  550. 

  551.              nullify(field%modified_mask)
    552. 

  552.      
    553. 

  553.              call storage_put_field(field)
    554. 

  554.     
    555. 

  555.           end do
    556. 

  556.     
    557. 

  557.         end if
    558. 

  558.       end do
    559. 

  559.     
    560. 

  560.       ! Done reading all static fields for this domain
    561. 

  561.       call input_close()
    562. 

  562. 

  563.    end subroutine get_static_fields
    564. 

  564. 

  565. 

  566.    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    567. 

  567.    ! Name: process_single_met_time
    568. 

  568.    !
    569. 

  569.    ! Purpose:
    570. 

  570.    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    571. 

  571.    subroutine process_single_met_time(do_const_processing, &
    572. 

  572.                              temp_date, n, extra_row, extra_col, xlat, xlon, &
    573. 

  573.                              xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
    574. 

  574.                              title, dyn_opt, &
    575. 

  575.                              west_east_dim, south_north_dim, &
    576. 

  576.                              we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
    577. 

  577.                              we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
    578. 

  578.                              sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
    579. 

  579.                              we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
    580. 

  580.                              sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
    581. 

  581.                              got_this_field, &
    582. 

  582.                              map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
    583. 

  583.                              grid_id, parent_id, i_parent_start, &
    584. 

  584.                              j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
    585. 

  585.                              cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
    586. 

  586.                              pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
    587. 

  587.                              corner_lats, corner_lons, output_flags, process_bdy_width)
    588. 

  588.    
    589. 

  589.       use bitarray_module
    590. 

  590.       use gridinfo_module
    591. 

  591.       use interp_module
    592. 

  592.       use interp_option_module
    593. 

  593.       use llxy_module
    594. 

  594.       use misc_definitions_module
    595. 

  595.       use module_debug
    596. 

  596.       use output_module
    597. 

  597.       use parallel_module
    598. 

  598.       use read_met_module
    599. 

  599.       use rotate_winds_module
    600. 

  600.       use storage_module
    601. 

  601.    
    602. 

  602.       implicit none
    603. 

  603.    
    604. 

  604.       ! Arguments
    605. 

  605.       logical, intent(in) :: do_const_processing
    606. 

  606.       integer, intent(in) :: n, dyn_opt, west_east_dim, south_north_dim, map_proj, &
    607. 

  607.                  we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
    608. 

  608.                  we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
    609. 

  609.                  sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
    610. 

  610.                  we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
    611. 

  611.                  sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
    612. 

  612.                  is_water, is_lake, is_ice, is_urban, i_soilwater, &
    613. 

  613.                  grid_id, parent_id, i_parent_start, j_parent_start, &
    614. 

  614.                  i_parent_end, j_parent_end, parent_grid_ratio, sub_x, sub_y, num_land_cat, &
    615. 

  615.                  process_bdy_width
    616. 

  616. ! BUG: Should we be passing these around as pointers, or just declare them as arrays?
    617. 

  617.       real, pointer, dimension(:,:) :: landmask
    618. 

  618.       real, intent(in) :: dom_dx, dom_dy, cen_lat, moad_cen_lat, cen_lon, stand_lon, &
    619. 

  619.                  truelat1, truelat2, pole_lat, pole_lon
    620. 

  620.       real, dimension(16), intent(in) :: corner_lats, corner_lons
    621. 

  621.       real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
    622. 

  622.       logical, intent(in) :: extra_row, extra_col
    623. 

  623.       logical, dimension(:), intent(inout) :: got_this_field
    624. 

  624.       character (len=19), intent(in) :: temp_date
    625. 

  625.       character (len=128), intent(in) :: mminlu
    626. 

  626.       character (len=128), dimension(:), intent(inout) :: output_flags
    627. 

  627. 

  628. ! BUG: Move this constant to misc_definitions_module?
    629. 

  629. integer, parameter :: BDR_WIDTH = 3
    630. 

  630.    
    631. 

  631.       ! Local variables
    632. 

  632.       integer :: istatus, iqstatus, fg_idx, idx, idxt, i, j, bottom_top_dim, &
    633. 

  633.                  sm1, em1, sm2, em2, sm3, em3, &
    634. 

  634.                  sp1, ep1, sp2, ep2, sp3, ep3, &
    635. 

  635.                  sd1, ed1, sd2, ed2, sd3, ed3, &
    636. 

  636.                  u_idx, bdr_wdth
    637. 

  637.       integer :: nmet_flags
    638. 

  638.       integer :: num_metgrid_soil_levs
    639. 

  639.       integer, pointer, dimension(:) :: soil_levels
    640. 

  640.       real :: rx, ry
    641. 

  641.       real :: threshold
    642. 

  642.       logical :: do_gcell_interp
    643. 

  643.       integer, pointer, dimension(:) :: u_levels, v_levels
    644. 

  644.       real, pointer, dimension(:,:) :: halo_slab
    645. 

  645.       real, pointer, dimension(:,:,:) :: real_array
    646. 

  646.       character (len=19) :: output_date
    647. 

  647.       character (len=128) :: cname, title
    648. 

  648.       character (len=MAX_FILENAME_LEN) :: input_name
    649. 

  649.       character (len=128), allocatable, dimension(:) :: met_flags
    650. 

  650.       type (fg_input) :: field, u_field, v_field
    651. 

  651.       type (met_data) :: fg_data
    652. 

  652. 

  653. 

  654.       ! For this time, we need to process all first-guess filename roots. When we 
    655. 

  655.       !   hit a root containing a '*', we assume we have hit the end of the list
    656. 

  656.       fg_idx = 1
    657. 

  657.       if (do_const_processing) then
    658. 

  658.          input_name = constants_name(fg_idx)
    659. 

  659.       else
    660. 

  660.          input_name = fg_name(fg_idx)
    661. 

  661.       end if
    662. 

  662.       do while (input_name /= '*')
    663. 

  663.    
    664. 

  664.          call mprintf(.true.,STDOUT, '    %s', s1=input_name)
    665. 

  665.          call mprintf(.true.,LOGFILE, 'Getting input fields from %s', s1=input_name)
    666. 

  666. 

  667.          ! Do a first pass through this fg source to get all mask fields used
    668. 

  668.          !   during interpolation
    669. 

  669.          call get_interp_masks(trim(input_name), do_const_processing, temp_date)
    670. 

  670.    
    671. 

  671.          istatus = 0
    672. 

  672. 

  673.          ! Initialize the module for reading in the met fields
    674. 

  674.          call read_met_init(trim(input_name), do_const_processing, temp_date, istatus)
    675. 

  675. 

  676.          if (istatus == 0) then
    677. 

  677.    
    678. 

  678.             ! Process all fields and levels from the current file; read_next_met_field()
    679. 

  679.             !   will return a non-zero status when there are no more fields to be read.
    680. 

  680.             do while (istatus == 0) 
    681. 

  681.       
    682. 

  682.                call read_next_met_field(fg_data, istatus)
    683. 

  683.       
    684. 

  684.                if (istatus == 0) then
    685. 

  685.       
    686. 

  686.                   ! Find index into fieldname, interp_method, masked, and fill_missing
    687. 

  687.                   !   of the current field
    688. 

  688.                   idxt = num_entries + 1
    689. 

  689.                   do idx=1,num_entries
    690. 

  690.                      if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
    691. 

  691.                          (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
    692. 

  692. 

  693.                         got_this_field(idx) = .true.
    694. 

  694. 

  695.                         if (index(input_name,trim(from_input(idx))) /= 0 .or. &
    696. 

  696.                            (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
    697. 

  697.                            idxt = idx
    698. 

  698.                         end if
    699. 

  699. 

  700.                      end if
    701. 

  701.                   end do
    702. 

  702.                   idx = idxt
    703. 

  703.                   if (idx > num_entries) idx = num_entries ! The last entry is a default
    704. 

  704. 

  705.                   ! Do we need to rename this field?
    706. 

  706.                   if (output_name(idx) /= ' ') then
    707. 

  707.                      fg_data%field = output_name(idx)(1:9)
    708. 

  708. 

  709.                      idxt = num_entries + 1
    710. 

  710.                      do idx=1,num_entries
    711. 

  711.                         if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
    712. 

  712.                             (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
    713. 

  713.    
    714. 

  714.                            got_this_field(idx) = .true.
    715. 

  715.    
    716. 

  716.                            if (index(input_name,trim(from_input(idx))) /= 0 .or. &
    717. 

  717.                               (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
    718. 

  718.                               idxt = idx
    719. 

  719.                            end if
    720. 

  720.    
    721. 

  721.                         end if
    722. 

  722.                      end do
    723. 

  723.                      idx = idxt
    724. 

  724.                      if (idx > num_entries) idx = num_entries ! The last entry is a default
    725. 

  725.                   end if
    726. 

  726. 

  727.                   ! Do a simple check to see whether this is a global dataset
    728. 

  728.                   ! Note that we do not currently support regional Gaussian grids
    729. 

  729.                   if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
    730. 

  730.                        .or. (fg_data%iproj == PROJ_GAUSS)) then
    731. 

  731.                      bdr_wdth = BDR_WIDTH
    732. 

  732.                      allocate(halo_slab(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
    733. 

  733. 

  734.                      halo_slab(1:fg_data%nx,                      1:fg_data%ny) = &
    735. 

  735.                                fg_data%slab(1:fg_data%nx,              1:fg_data%ny)
    736. 

  736. 

  737.                      halo_slab(1-BDR_WIDTH:0,                     1:fg_data%ny) = &
    738. 

  738.                                fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
    739. 

  739. 

  740.                      halo_slab(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
    741. 

  741.                                fg_data%slab(1:BDR_WIDTH,       1:fg_data%ny)
    742. 

  742. 

  743.                      deallocate(fg_data%slab)
    744. 

  744.                   else
    745. 

  745.                      bdr_wdth = 0
    746. 

  746.                      halo_slab => fg_data%slab
    747. 

  747.                      nullify(fg_data%slab)
    748. 

  748.                   end if
    749. 

  749. 

  750.                   call mprintf(.true.,LOGFILE,'Processing %s at level %f.',s1=fg_data%field,f1=fg_data%xlvl)               
    751. 

  751.    
    752. 

  752.                   call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
    753. 

  753.                                               fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
    754. 

  754.                                               fg_data%deltalon, fg_data%starti, fg_data%startj, &
    755. 

  755.                                               fg_data%startlat, fg_data%startlon, earth_radius=fg_data%earth_radius*1000.)
    756. 

  756.       
    757. 

  757.                   ! Initialize fg_input structure to store the field
    758. 

  758.                   field%header%version = 1
    759. 

  759.                   field%header%date = fg_data%hdate//'        '
    760. 

  760.                   if (do_const_processing) then
    761. 

  761.                      field%header%time_dependent = .false.
    762. 

  762.                   else
    763. 

  763.                      field%header%time_dependent = .true.
    764. 

  764.                   end if
    765. 

  765.                   field%header%forecast_hour = fg_data%xfcst 
    766. 

  766.                   field%header%fg_source = 'FG'
    767. 

  767.                   field%header%field = ' '
    768. 

  768.                   field%header%field(1:9) = fg_data%field
    769. 

  769.                   field%header%units = ' '
    770. 

  770.                   field%header%units(1:25) = fg_data%units
    771. 

  771.                   field%header%description = ' '
    772. 

  772.                   field%header%description(1:46) = fg_data%desc
    773. 

  773.                   call get_z_dim_name(fg_data%field,field%header%vertical_coord)
    774. 

  774.                   field%header%vertical_level = nint(fg_data%xlvl) 
    775. 

  775.                   field%header%sr_x = sub_x
    776. 

  776.                   field%header%sr_y = sub_y
    777. 

  777.                   field%header%is_subgrid = .false.
    778. 

  778.                   field%header%array_order = 'XY ' 
    779. 

  779.                   field%header%is_wind_grid_rel = fg_data%is_wind_grid_rel 
    780. 

  780.                   field%header%array_has_missing_values = .false.
    781. 

  781.                   nullify(field%r_arr)
    782. 

  782.                   nullify(field%valid_mask)
    783. 

  783.                   nullify(field%modified_mask)
    784. 

  784. 

  785.                   if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
    786. 

  786.                      output_flags(idx) = flag_in_output(idx)
    787. 

  787.                   end if
    788. 

  788. 

  789.                   ! If we should not output this field, just list it as a mask field
    790. 

  790.                   if (output_this_field(idx)) then
    791. 

  791.                      field%header%mask_field = .false.
    792. 

  792.                   else
    793. 

  793.                      field%header%mask_field = .true.
    794. 

  794.                   end if
    795. 

  795.       
    796. 

  796.                   !
    797. 

  797.                   ! Before actually doing any interpolation to the model grid, we must check
    798. 

  798.                   !    whether we will be using the average_gcell interpolator that averages all 
    799. 

  799.                   !    source points in each model grid cell
    800. 

  800.                   !
    801. 

  801.                   do_gcell_interp = .false.
    802. 

  802.                   if (index(interp_method(idx),'average_gcell') /= 0) then
    803. 

  803.       
    804. 

  804.                      call get_gcell_threshold(interp_method(idx), threshold, istatus)
    805. 

  805.                      if (istatus == 0) then
    806. 

  806.                         if (fg_data%dx == 0. .and. fg_data%dy == 0. .and. &
    807. 

  807.                             fg_data%deltalat /= 0. .and. fg_data%deltalon /= 0.) then
    808. 

  808.                            fg_data%dx = abs(fg_data%deltalon)
    809. 

  809.                            fg_data%dy = abs(fg_data%deltalat)
    810. 

  810.                         else
    811. 

  811. ! BUG: Need to more correctly handle dx/dy in meters.
    812. 

  812.                            fg_data%dx = fg_data%dx / 111000.  ! Convert meters to approximate degrees
    813. 

  813.                            fg_data%dy = fg_data%dy / 111000.
    814. 

  814.                         end if
    815. 

  815.                         if (gridtype == 'C') then
    816. 

  816.                            if (threshold*max(fg_data%dx,fg_data%dy)*111. <= max(dom_dx,dom_dy)/1000.) &
    817. 

  817.                               do_gcell_interp = .true. 
    818. 

  818.                         else if (gridtype == 'E') then
    819. 

  819.                            if (threshold*max(fg_data%dx,fg_data%dy) <= max(dom_dx,dom_dy)) &
    820. 

  820.                               do_gcell_interp = .true. 
    821. 

  821.                         end if
    822. 

  822.                      end if
    823. 

  823.                   end if
    824. 

  824.       
    825. 

  825.                   ! Interpolate to U staggering
    826. 

  826.                   if (output_stagger(idx) == U) then
    827. 

  827.    
    828. 

  828.                      call storage_query_field(field, iqstatus)
    829. 

  829.                      if (iqstatus == 0) then
    830. 

  830.                         call storage_get_field(field, iqstatus)
    831. 

  831.                         call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
    832. 

  832.                                      ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
    833. 

  833.                         if (associated(field%modified_mask)) then
    834. 

  834.                            call bitarray_destroy(field%modified_mask)
    835. 

  835.                            nullify(field%modified_mask)
    836. 

  836.                         end if
    837. 

  837.                      else
    838. 

  838.                         allocate(field%valid_mask)
    839. 

  839.                         call bitarray_create(field%valid_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
    840. 

  840.                      end if
    841. 

  841.       
    842. 

  842.                      ! Save a copy of the fg_input structure for the U field so that we can find it later
    843. 

  843.                      if (is_u_field(idx)) call dup(field, u_field)
    844. 

  844.    
    845. 

  845.                      allocate(field%modified_mask)
    846. 

  846.                      call bitarray_create(field%modified_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
    847. 

  847.    
    848. 

  848.                      if (do_const_processing .or. field%header%time_dependent) then
    849. 

  849.                         call interp_met_field(input_name, fg_data%field, U, M, &
    850. 

  850.                                      field, xlat_u, xlon_u, we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, &
    851. 

  851.                                      we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
    852. 

  852.                                      halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
    853. 

  853.                                      field%modified_mask, process_bdy_width)
    854. 

  854.                      else
    855. 

  855.                         call mprintf(.true.,INFORM,' - already processed this field from constant file.')
    856. 

  856.                      end if
    857. 

  857.    
    858. 

  858.                   ! Interpolate to V staggering
    859. 

  859.                   else if (output_stagger(idx) == V) then
    860. 

  860.    
    861. 

  861.                      call storage_query_field(field, iqstatus)
    862. 

  862.                      if (iqstatus == 0) then
    863. 

  863.                         call storage_get_field(field, iqstatus)
    864. 

  864.                         call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
    865. 

  865.                                      ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
    866. 

  866.                         if (associated(field%modified_mask)) then
    867. 

  867.                            call bitarray_destroy(field%modified_mask)
    868. 

  868.                            nullify(field%modified_mask)
    869. 

  869.                         end if
    870. 

  870.                      else
    871. 

  871.                         allocate(field%valid_mask)
    872. 

  872.                         call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
    873. 

  873.                      end if
    874. 

  874.       
    875. 

  875.                      ! Save a copy of the fg_input structure for the V field so that we can find it later
    876. 

  876.                      if (is_v_field(idx)) call dup(field, v_field)
    877. 

  877.    
    878. 

  878.                      allocate(field%modified_mask)
    879. 

  879.                      call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
    880. 

  880.    
    881. 

  881.                      if (do_const_processing .or. field%header%time_dependent) then
    882. 

  882.                         call interp_met_field(input_name, fg_data%field, V, M, &
    883. 

  883.                                      field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
    884. 

  884.                                      we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
    885. 

  885.                                      halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
    886. 

  886.                                      field%modified_mask, process_bdy_width)
    887. 

  887.                      else
    888. 

  888.                         call mprintf(.true.,INFORM,' - already processed this field from constant file.')
    889. 

  889.                      end if
    890. 

  890.              
    891. 

  891.                   ! Interpolate to VV staggering
    892. 

  892.                   else if (output_stagger(idx) == VV) then
    893. 

  893.    
    894. 

  894.                      call storage_query_field(field, iqstatus)
    895. 

  895.                      if (iqstatus == 0) then
    896. 

  896.                         call storage_get_field(field, iqstatus)
    897. 

  897.                         call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
    898. 

  898.                                      ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
    899. 

  899.                         if (associated(field%modified_mask)) then
    900. 

  900.                            call bitarray_destroy(field%modified_mask)
    901. 

  901.                            nullify(field%modified_mask)
    902. 

  902.                         end if
    903. 

  903.                      else
    904. 

  904.                         allocate(field%valid_mask)
    905. 

  905.                         call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
    906. 

  906.                      end if
    907. 

  907.       
    908. 

  908.                      ! Save a copy of the fg_input structure for the U field so that we can find it later
    909. 

  909.                      if (is_u_field(idx)) call dup(field, u_field)
    910. 

  910. 

  911.                      ! Save a copy of the fg_input structure for the V field so that we can find it later
    912. 

  912.                      if (is_v_field(idx)) call dup(field, v_field)
    913. 

  913.    
    914. 

  914.                      allocate(field%modified_mask)
    915. 

  915.                      call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
    916. 

  916.    
    917. 

  917.                      if (do_const_processing .or. field%header%time_dependent) then
    918. 

  918.                         call interp_met_field(input_name, fg_data%field, VV, M, &
    919. 

  919.                                      field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
    920. 

  920.                                      we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
    921. 

  921.                                      halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
    922. 

  922.                                      field%modified_mask, process_bdy_width)
    923. 

  923.                      else
    924. 

  924.                         call mprintf(.true.,INFORM,' - already processed this field from constant file.')
    925. 

  925.                      end if
    926. 

  926.    
    927. 

  927.                   ! All other fields interpolated to M staggering for C grid, H staggering for E grid
    928. 

  928.                   else
    929. 

  929.    
    930. 

  930.                      call storage_query_field(field, iqstatus)
    931. 

  931.                      if (iqstatus == 0) then
    932. 

  932.                         call storage_get_field(field, iqstatus)
    933. 

  933.                         call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
    934. 

  934.                                      ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
    935. 

  935.                         if (associated(field%modified_mask)) then
    936. 

  936.                            call bitarray_destroy(field%modified_mask)
    937. 

  937.                            nullify(field%modified_mask)
    938. 

  938.                         end if
    939. 

  939.                      else
    940. 

  940.                         allocate(field%valid_mask)
    941. 

  941.                         call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
    942. 

  942.                      end if
    943. 

  943.    
    944. 

  944.                      allocate(field%modified_mask)
    945. 

  945.                      call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
    946. 

  946.    
    947. 

  947.                      if (do_const_processing .or. field%header%time_dependent) then
    948. 

  948.                         if (gridtype == 'C') then
    949. 

  949.                            call interp_met_field(input_name, fg_data%field, M, M, &
    950. 

  950.                                         field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
    951. 

  951.                                         we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
    952. 

  952.                                         halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
    953. 

  953.                                         field%modified_mask, process_bdy_width, landmask)
    954. 

  954.       
    955. 

  955.                         else if (gridtype == 'E') then
    956. 

  956.                            call interp_met_field(input_name, fg_data%field, HH, M, &
    957. 

  957.                                         field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
    958. 

  958.                                         we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
    959. 

  959.                                         halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
    960. 

  960.                                         field%modified_mask, process_bdy_width, landmask)
    961. 

  961.                         end if
    962. 

  962.                      else
    963. 

  963.                         call mprintf(.true.,INFORM,' - already processed this field from constant file.')
    964. 

  964.                      end if
    965. 

  965.    
    966. 

  966.                   end if
    967. 

  967.    
    968. 

  968.                   call bitarray_merge(field%valid_mask, field%modified_mask)
    969. 

  969. 

  970.                   deallocate(halo_slab)
    971. 

  971.                                
    972. 

  972.                   ! Store the interpolated field
    973. 

  973.                   call storage_put_f…
    974.

Large files files are truncated, but you can click here to view the full file