/WPS/geogrid/src/process_tile_module.F
FORTRAN Legacy | 2084 lines | 1265 code | 310 blank | 509 comment | 366 complexity | 039d3d2e0118ec2527dcd4f13fac8fbe MD5 | raw file
Possible License(s): AGPL-1.0
Large files files are truncated, but you can click here to view the full file
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Module: process_tile
- !
- ! Description:
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- module process_tile_module
- use module_debug
- contains
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Name: process_tile
- !
- ! Purpose: To process a tile, whose lower-left corner is at
- ! (tile_i_min, tile_j_min) and whose upper-right corner is at
- ! (tile_i_max, tile_j_max), of the model grid given by which_domain
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- subroutine process_tile(which_domain, grid_type, dynopt, &
- dummy_start_dom_i, dummy_end_dom_i, &
- dummy_start_dom_j, dummy_end_dom_j, &
- dummy_start_patch_i, dummy_end_patch_i, &
- dummy_start_patch_j, dummy_end_patch_j, &
- extra_col, extra_row)
-
- use bitarray_module
- use hash_module
- use llxy_module
- use misc_definitions_module
- use output_module
- use smooth_module
- use source_data_module
- use gridinfo_module, only : subgrid_ratio_x, subgrid_ratio_y
-
- implicit none
- ! Arguments
- integer, intent(in) :: which_domain, dynopt, &
- dummy_start_dom_i, dummy_end_dom_i, dummy_start_dom_j, dummy_end_dom_j, &
- dummy_start_patch_i, dummy_end_patch_i, dummy_start_patch_j, dummy_end_patch_j
- logical, intent(in) :: extra_col, extra_row
- character (len=1), intent(in) :: grid_type
-
- ! Local variables
- integer :: i, j, k, kk, istatus, ifieldstatus, idomcatstatus, field_count
- integer :: min_category, max_category, min_level, max_level, &
- smth_opt, smth_passes, num_landmask_categories
- integer :: start_dom_i, end_dom_i, start_dom_j, end_dom_j, end_dom_stag_i, end_dom_stag_j
- integer :: start_patch_i, end_patch_i, start_patch_j, end_patch_j, end_patch_stag_i, end_patch_stag_j
- integer :: start_mem_i, end_mem_i, start_mem_j, end_mem_j, end_mem_stag_i, end_mem_stag_j
- integer :: start_mem_sub_i, end_mem_sub_i, start_mem_sub_j, end_mem_sub_j
- integer :: sm1, em1, sm2, em2
- integer :: istagger
- integer, dimension(MAX_LANDMASK_CATEGORIES) :: landmask_value
- real :: sum, dominant, msg_fill_val, topo_flag_val, mass_flag, land_total, water_total
- real, dimension(16) :: corner_lats, corner_lons
- real, pointer, dimension(:,:) :: xlat_array, xlon_array, &
- xlat_array_u, xlon_array_u, &
- xlat_array_v, xlon_array_v, &
- clat_array, clon_array, &
- xlat_array_subgrid, xlon_array_subgrid, &
- f_array, e_array, &
- mapfac_array_m_x, mapfac_array_u_x, mapfac_array_v_x, &
- mapfac_array_m_y, mapfac_array_u_y, mapfac_array_v_y, &
- mapfac_array_x_subgrid, mapfac_array_y_subgrid, &
- sina_array, cosa_array
- real, pointer, dimension(:,:) :: xlat_ptr, xlon_ptr, mapfac_ptr_x, mapfac_ptr_y, landmask, dominant_field
- real, pointer, dimension(:,:,:) :: field, slp_field
- logical :: is_water_mask, only_save_dominant, halt_on_missing
- character (len=19) :: datestr
- character (len=128) :: fieldname, gradname, domname, landmask_name
- character (len=256) :: temp_string
- type (bitarray) :: processed_domain
- type (hashtable) :: processed_fieldnames
- character (len=128), dimension(2) :: dimnames
- integer :: sr_x, sr_y
- logical :: subgrid_var
-
- sr_x=subgrid_ratio_x(which_domain)
- sr_y=subgrid_ratio_y(which_domain)
- datestr = '0000-00-00_00:00:00'
- field_count = 0
- mass_flag=1.0
-
- ! The following pertains primarily to the C grid
- ! Determine whether only (n-1)th rows/columns should be computed for variables
- ! on staggered grid. In a distributed memory situation, not every tile should
- ! have only (n-1)th rows/columns computed, or we end up with (n-k)
- ! rows/columns when there are k patches in the y/x direction
- if (extra_col) then
- start_patch_i = dummy_start_patch_i ! The seemingly pointless renaming of start
- end_patch_i = dummy_end_patch_i - 1 ! naming convention with modified end_patch variables,
- end_patch_stag_i = dummy_end_patch_i ! variables is so that we can maintain consistent
- ! which are marked as intent(in)
- start_mem_i = start_patch_i - HALO_WIDTH
- end_mem_i = end_patch_i + HALO_WIDTH
- end_mem_stag_i = end_patch_stag_i + HALO_WIDTH
- else
- start_patch_i = dummy_start_patch_i
- end_patch_i = dummy_end_patch_i
- end_patch_stag_i = dummy_end_patch_i
- start_mem_i = start_patch_i - HALO_WIDTH
- end_mem_i = end_patch_i + HALO_WIDTH
- end_mem_stag_i = end_patch_stag_i + HALO_WIDTH
- end if
-
- if (extra_row) then
- start_patch_j = dummy_start_patch_j
- end_patch_j = dummy_end_patch_j - 1
- end_patch_stag_j = dummy_end_patch_j
- start_mem_j = start_patch_j - HALO_WIDTH
- end_mem_j = end_patch_j + HALO_WIDTH
- end_mem_stag_j = end_patch_stag_j + HALO_WIDTH
- else
- start_patch_j = dummy_start_patch_j
- end_patch_j = dummy_end_patch_j
- end_patch_stag_j = dummy_end_patch_j
- start_mem_j = start_patch_j - HALO_WIDTH
- end_mem_j = end_patch_j + HALO_WIDTH
- end_mem_stag_j = end_patch_stag_j + HALO_WIDTH
- end if
- start_dom_i = dummy_start_dom_i
- if (grid_type == 'C') then
- end_dom_i = dummy_end_dom_i - 1
- end_dom_stag_i = dummy_end_dom_i
- else if (grid_type == 'E') then
- end_dom_i = dummy_end_dom_i
- end_dom_stag_i = dummy_end_dom_i
- end if
- start_dom_j = dummy_start_dom_j
- if (grid_type == 'C') then
- end_dom_j = dummy_end_dom_j - 1
- end_dom_stag_j = dummy_end_dom_j
- else if (grid_type == 'E') then
- end_dom_j = dummy_end_dom_j
- end_dom_stag_j = dummy_end_dom_j
- end if
- start_mem_sub_i = (start_mem_i - 1 ) * sr_x + 1
- if(.not.extra_col)then
- end_mem_sub_i = (end_mem_i) * sr_x
- else
- end_mem_sub_i = (end_mem_i + 1) * sr_x
- end if
- start_mem_sub_j = (start_mem_j -1 ) * sr_y + 1
- if(.not.extra_row)then
- end_mem_sub_j = (end_mem_j) * sr_y
- else
- end_mem_sub_j = (end_mem_j + 1) * sr_y
- end if
-
- ! Allocate arrays to hold all lat/lon fields; these will persist for the duration of
- ! the process_tile routine
- ! For C grid, we have M, U, and V points
- ! For E grid, we have only M and V points
- allocate(xlat_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
- allocate(xlon_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
- allocate(xlat_array_v(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
- allocate(xlon_array_v(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
- if (grid_type == 'C') then
- allocate(xlat_array_u(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
- allocate(xlon_array_u(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
- allocate(clat_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
- allocate(clon_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
- allocate(xlat_array_subgrid(start_mem_sub_i:end_mem_sub_i,start_mem_sub_j:end_mem_sub_j))
- allocate(xlon_array_subgrid(start_mem_sub_i:end_mem_sub_i,start_mem_sub_j:end_mem_sub_j))
- allocate(mapfac_array_x_subgrid(start_mem_sub_i:end_mem_sub_i,start_mem_sub_j:end_mem_sub_j))
- allocate(mapfac_array_y_subgrid(start_mem_sub_i:end_mem_sub_i,start_mem_sub_j:end_mem_sub_j))
- end if
-
- ! Initialize hash table to track which fields have been processed
- call hash_init(processed_fieldnames)
-
- !
- ! Calculate lat/lon for every point in the tile (XLAT and XLON)
- ! The xlat_array and xlon_array arrays will be used in processing other fields
- !
- call mprintf(.true.,STDOUT,' Processing XLAT and XLONG')
-
- if (grid_type == 'C') then
- call get_lat_lon_fields(xlat_array, xlon_array, start_mem_i, &
- start_mem_j, end_mem_i, end_mem_j, M)
- call get_lat_lon_fields(xlat_array_v, xlon_array_v, start_mem_i, &
- start_mem_j, end_mem_i, end_mem_stag_j, V)
- call get_lat_lon_fields(xlat_array_u, xlon_array_u, start_mem_i, &
- start_mem_j, end_mem_stag_i, end_mem_j, U)
- call get_lat_lon_fields(clat_array, clon_array, start_mem_i, &
- start_mem_j, end_mem_i, end_mem_j, M, comp_ll=.true.)
- call get_lat_lon_fields(xlat_array_subgrid, xlon_array_subgrid, &
- start_mem_sub_i, start_mem_sub_j, end_mem_sub_i, end_mem_sub_j, M, sub_x=sr_x, sub_y=sr_y)
- call get_map_factor(xlat_array_subgrid, xlon_array_subgrid, mapfac_array_x_subgrid, &
- mapfac_array_y_subgrid, start_mem_sub_i, start_mem_sub_j, end_mem_sub_i, end_mem_sub_j)
- corner_lats(1) = xlat_array(start_patch_i,start_patch_j)
- corner_lats(2) = xlat_array(start_patch_i,end_patch_j)
- corner_lats(3) = xlat_array(end_patch_i,end_patch_j)
- corner_lats(4) = xlat_array(end_patch_i,start_patch_j)
-
- corner_lats(5) = xlat_array_u(start_patch_i,start_patch_j)
- corner_lats(6) = xlat_array_u(start_patch_i,end_patch_j)
- corner_lats(7) = xlat_array_u(end_patch_stag_i,end_patch_j)
- corner_lats(8) = xlat_array_u(end_patch_stag_i,start_patch_j)
-
- corner_lats(9) = xlat_array_v(start_patch_i,start_patch_j)
- corner_lats(10) = xlat_array_v(start_patch_i,end_patch_stag_j)
- corner_lats(11) = xlat_array_v(end_patch_i,end_patch_stag_j)
- corner_lats(12) = xlat_array_v(end_patch_i,start_patch_j)
-
- call xytoll(real(start_patch_i)-0.5, real(start_patch_j)-0.5, corner_lats(13), corner_lons(13), M)
- call xytoll(real(start_patch_i)-0.5, real(end_patch_j)+0.5, corner_lats(14), corner_lons(14), M)
- call xytoll(real(end_patch_i)+0.5, real(end_patch_j)+0.5, corner_lats(15), corner_lons(15), M)
- call xytoll(real(end_patch_i)+0.5, real(start_patch_j)-0.5, corner_lats(16), corner_lons(16), M)
- corner_lons(1) = xlon_array(start_patch_i,start_patch_j)
- corner_lons(2) = xlon_array(start_patch_i,end_patch_j)
- corner_lons(3) = xlon_array(end_patch_i,end_patch_j)
- corner_lons(4) = xlon_array(end_patch_i,start_patch_j)
-
- corner_lons(5) = xlon_array_u(start_patch_i,start_patch_j)
- corner_lons(6) = xlon_array_u(start_patch_i,end_patch_j)
- corner_lons(7) = xlon_array_u(end_patch_stag_i,end_patch_j)
- corner_lons(8) = xlon_array_u(end_patch_stag_i,start_patch_j)
-
- corner_lons(9) = xlon_array_v(start_patch_i,start_patch_j)
- corner_lons(10) = xlon_array_v(start_patch_i,end_patch_stag_j)
- corner_lons(11) = xlon_array_v(end_patch_i,end_patch_stag_j)
- corner_lons(12) = xlon_array_v(end_patch_i,start_patch_j)
-
- else if (grid_type == 'E') then
- call get_lat_lon_fields(xlat_array, xlon_array, start_mem_i, &
- start_mem_j, end_mem_i, end_mem_j, HH)
- call get_lat_lon_fields(xlat_array_v, xlon_array_v, start_mem_i, &
- start_mem_j, end_mem_i, end_mem_stag_j, VV)
-
- corner_lats(1) = xlat_array(start_patch_i,start_patch_j)
- corner_lats(2) = xlat_array(start_patch_i,end_patch_j)
- corner_lats(3) = xlat_array(end_patch_i,end_patch_j)
- corner_lats(4) = xlat_array(end_patch_i,start_patch_j)
-
- corner_lats(5) = xlat_array_v(start_patch_i,start_patch_j)
- corner_lats(6) = xlat_array_v(start_patch_i,end_patch_stag_j)
- corner_lats(7) = xlat_array_v(end_patch_i,end_patch_stag_j)
- corner_lats(8) = xlat_array_v(end_patch_i,start_patch_j)
-
- corner_lats(9) = 0.0
- corner_lats(10) = 0.0
- corner_lats(11) = 0.0
- corner_lats(12) = 0.0
-
- corner_lats(13) = 0.0
- corner_lats(14) = 0.0
- corner_lats(15) = 0.0
- corner_lats(16) = 0.0
-
- corner_lons(1) = xlon_array(start_patch_i,start_patch_j)
- corner_lons(2) = xlon_array(start_patch_i,end_patch_j)
- corner_lons(3) = xlon_array(end_patch_i,end_patch_j)
- corner_lons(4) = xlon_array(end_patch_i,start_patch_j)
-
- corner_lons(5) = xlon_array_v(start_patch_i,start_patch_j)
- corner_lons(6) = xlon_array_v(start_patch_i,end_patch_stag_j)
- corner_lons(7) = xlon_array_v(end_patch_i,end_patch_stag_j)
- corner_lons(8) = xlon_array_v(end_patch_i,start_patch_j)
-
- corner_lons(9) = 0.0
- corner_lons(10) = 0.0
- corner_lons(11) = 0.0
- corner_lons(12) = 0.0
-
- corner_lons(13) = 0.0
- corner_lons(14) = 0.0
- corner_lons(15) = 0.0
- corner_lons(16) = 0.0
-
- end if
- ! Initialize the output module now that we have the corner point lats/lons
- call output_init(which_domain, 'OUTPUT FROM GEOGRID V3.3', '0000-00-00_00:00:00', grid_type, dynopt, &
- corner_lats, corner_lons, &
- start_dom_i, end_dom_i, start_dom_j, end_dom_j, &
- start_patch_i, end_patch_i, start_patch_j, end_patch_j, &
- start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
- extra_col, extra_row)
-
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
- 'XLAT_M', datestr, real_array = xlat_array)
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
- 'XLONG_M', datestr, real_array = xlon_array)
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, &
- 'XLAT_V', datestr, real_array = xlat_array_v)
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, &
- 'XLONG_V', datestr, real_array = xlon_array_v)
- if (grid_type == 'C') then
- call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, &
- 'XLAT_U', datestr, real_array = xlat_array_u)
- call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, &
- 'XLONG_U', datestr, real_array = xlon_array_u)
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
- 'CLAT', datestr, real_array = clat_array)
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
- 'CLONG', datestr, real_array = clon_array)
- call write_field(start_mem_sub_i, end_mem_sub_i, start_mem_sub_j, end_mem_sub_j, 1, 1, &
- 'FXLAT', datestr, real_array = xlat_array_subgrid)
- call write_field(start_mem_sub_i, end_mem_sub_i, start_mem_sub_j, end_mem_sub_j, 1, 1, &
- 'FXLONG', datestr, real_array = xlon_array_subgrid)
- if (associated(clat_array)) deallocate(clat_array)
- if (associated(clon_array)) deallocate(clon_array)
- end if
-
- !
- ! Calculate map factor for current domain
- !
- if (grid_type == 'C') then
- call mprintf(.true.,STDOUT,' Processing MAPFAC')
-
- allocate(mapfac_array_m_x(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
- allocate(mapfac_array_m_y(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
- call get_map_factor(xlat_array, xlon_array, mapfac_array_m_x, mapfac_array_m_y, start_mem_i, &
- start_mem_j, end_mem_i, end_mem_j)
- ! Global WRF uses map scale factors in X and Y directions, but "regular" WRF uses a single MSF
- ! on each staggering. In the case of regular WRF, we can assume that MAPFAC_MX = MAPFAC_MY = MAPFAC_M,
- ! and so we can simply write MAPFAC_MX as the MAPFAC_M field. Ultimately, when global WRF is
- ! merged into the WRF trunk, we will need only two map scale factor fields for each staggering,
- ! in the x and y directions, and these will be the same in the case of non-Cassini projections
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_M', &
- datestr, real_array = mapfac_array_m_x)
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_MX', &
- datestr, real_array = mapfac_array_m_x)
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_MY', &
- datestr, real_array = mapfac_array_m_y)
-
- allocate(mapfac_array_v_x(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
- allocate(mapfac_array_v_y(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
- call get_map_factor(xlat_array_v, xlon_array_v, mapfac_array_v_x, mapfac_array_v_y, start_mem_i, &
- start_mem_j, end_mem_i, end_mem_stag_j)
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_V', &
- datestr, real_array = mapfac_array_v_x)
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_VX', &
- datestr, real_array = mapfac_array_v_x)
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_VY', &
- datestr, real_array = mapfac_array_v_y)
-
- allocate(mapfac_array_u_x(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
- allocate(mapfac_array_u_y(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
- call get_map_factor(xlat_array_u, xlon_array_u, mapfac_array_u_x, mapfac_array_u_y, start_mem_i, &
- start_mem_j, end_mem_stag_i, end_mem_j)
- call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_U', &
- datestr, real_array = mapfac_array_u_x)
- call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_UX', &
- datestr, real_array = mapfac_array_u_x)
- call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_UY', &
- datestr, real_array = mapfac_array_u_y)
- end if
-
-
- !
- ! Coriolis parameters (E and F)
- !
- call mprintf(.true.,STDOUT,' Processing F and E')
-
- allocate(f_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
- allocate(e_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
-
- call get_coriolis_parameters(xlat_array, f_array, e_array, &
- start_mem_i, start_mem_j, end_mem_i, end_mem_j)
-
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'E', &
- datestr, real_array = e_array)
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'F', &
- datestr, real_array = f_array)
-
- if (associated(f_array)) deallocate(f_array)
- if (associated(e_array)) deallocate(e_array)
-
-
- !
- ! Rotation angle (SINALPHA and COSALPHA)
- !
- if (grid_type == 'C') then
- call mprintf(.true.,STDOUT,' Processing ROTANG')
-
- allocate(sina_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
- allocate(cosa_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
-
- call get_rotang(xlat_array, xlon_array, cosa_array, sina_array, &
- start_mem_i, start_mem_j, end_mem_i, end_mem_j)
-
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'SINALPHA', &
- datestr, real_array = sina_array)
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'COSALPHA', &
- datestr, real_array = cosa_array)
-
- if (associated(sina_array)) deallocate(sina_array)
- if (associated(cosa_array)) deallocate(cosa_array)
- end if
-
- ! Every field up until now should probably just be processed regardless of what the user
- ! has specified for fields to be processed.
- ! Hereafter, we process user-specified fields
-
- !
- ! First process the field that we will derive a landmask from
- !
- call get_landmask_field(geog_data_res(which_domain), landmask_name, is_water_mask, landmask_value, istatus)
- do kk=1,MAX_LANDMASK_CATEGORIES
- if (landmask_value(kk) == INVALID) then
- num_landmask_categories = kk-1
- exit
- end if
- end do
- if (kk > MAX_LANDMASK_CATEGORIES) num_landmask_categories = MAX_LANDMASK_CATEGORIES
- if (istatus /= 0) then
- call mprintf(.true.,WARN,'No field specified for landmask calculation. Will set landmask=1 at every grid point.')
-
- allocate(landmask(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
- landmask = 1.
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'LANDMASK', &
- datestr, landmask)
-
- else
-
- allocate(landmask(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
- landmask = 1.
-
- call mprintf(.true.,STDOUT,' Processing %s', s1=trim(landmask_name))
-
- call get_missing_fill_value(landmask_name, msg_fill_val, istatus)
- if (istatus /= 0) msg_fill_val = NAN
-
- call get_halt_on_missing(landmask_name, halt_on_missing, istatus)
- if (istatus /= 0) halt_on_missing = .false.
-
- ! Do we calculate a dominant category for this field?
- call get_domcategory_name(landmask_name, domname, only_save_dominant, idomcatstatus)
-
- temp_string = ' '
- temp_string(1:128) = landmask_name
- call hash_insert(processed_fieldnames, temp_string)
-
- call get_max_categories(landmask_name, min_category, max_category, istatus)
- allocate(field(start_mem_i:end_mem_i, start_mem_j:end_mem_j, min_category:max_category))
- if (.not. only_save_dominant) then
- field_count = field_count + 1
- call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
- i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=landmask_name)
- else
- field_count = field_count + 1
- call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
- i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
- end if
- if (grid_type == 'C') then
- call calc_field(landmask_name, field, xlat_array, xlon_array, M, &
- start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
- min_category, max_category, processed_domain, 1, landmask=landmask, sr_x=1, sr_y=1)
- else if (grid_type == 'E') then
- call calc_field(landmask_name, field, xlat_array, xlon_array, HH, &
- start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
- min_category, max_category, processed_domain, 1, landmask=landmask, sr_x=1, sr_y=1)
- end if
-
- ! Find fractions
- do i=start_mem_i, end_mem_i
- do j=start_mem_j, end_mem_j
- sum = 0.
- do k=min_category,max_category
- sum = sum + field(i,j,k)
- end do
- if (sum > 0.0) then
- do k=min_category,max_category
- field(i,j,k) = field(i,j,k) / sum
- end do
- else
- do k=min_category,max_category
- field(i,j,k) = msg_fill_val
- end do
- end if
- end do
- end do
-
- ! If user wants to halt when a missing value is found in output field, check now
- if (halt_on_missing) then
- do i=start_mem_i, end_mem_i
- do j=start_mem_j, end_mem_j
- ! Only need to examine k=1
- if (field(i,j,1) == msg_fill_val) then
- call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
- end if
- end do
- end do
- end if
-
- if (is_water_mask) then
- call mprintf(.true.,STDOUT,' Calculating landmask from %s ( WATER =', &
- newline=.false.,s1=trim(landmask_name))
- else
- call mprintf(.true.,STDOUT,' Calculating landmask from %s ( LAND =', &
- newline=.false.,s1=trim(landmask_name))
- end if
- do k = 1, num_landmask_categories
- call mprintf(.true.,STDOUT,' %i',newline=.false.,i1=landmask_value(k))
- if (k == num_landmask_categories) call mprintf(.true.,STDOUT,')')
- end do
-
- ! Calculate landmask
- if (is_water_mask) then
- do i=start_mem_i, end_mem_i
- do j=start_mem_j, end_mem_j
- water_total = -1.
- do k=1,num_landmask_categories
- if (landmask_value(k) >= min_category .and. landmask_value(k) <= max_category) then
- if (field(i,j,landmask_value(k)) /= msg_fill_val) then
- if (water_total < 0.) water_total = 0.
- water_total = water_total + field(i,j,landmask_value(k))
- else
- water_total = -1.
- exit
- end if
- end if
- end do
- if (water_total >= 0.0) then
- if (water_total < 0.50) then
- landmask(i,j) = 1.
- else
- landmask(i,j) = 0.
- end if
- else
- landmask(i,j) = -1.
- end if
- end do
- end do
- else
- do i=start_mem_i, end_mem_i
- do j=start_mem_j, end_mem_j
- land_total = -1.
- do k=1,num_landmask_categories
- if (landmask_value(k) >= min_category .and. landmask_value(k) <= max_category) then
- if (field(i,j,landmask_value(k)) /= msg_fill_val) then
- if (land_total < 0.) land_total = 0.
- land_total = land_total + field(i,j,landmask_value(k))
- else
- land_total = -1.
- exit
- end if
- end if
- end do
- if (land_total >= 0.0) then
- if (land_total > 0.50) then
- landmask(i,j) = 1.
- else
- landmask(i,j) = 0.
- end if
- else
- landmask(i,j) = -1.
- end if
- end do
- end do
- end if
-
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'LANDMASK', &
- datestr, landmask)
-
- ! If we should only save the dominant category, then no need to write out fractional field
- if (.not.only_save_dominant .or. (idomcatstatus /= 0)) then
-
- ! Finally, we may be asked to smooth the fractional field
- call get_smooth_option(landmask_name, smth_opt, smth_passes, istatus)
- if (istatus == 0) then
- if (grid_type == 'C') then
- if (smth_opt == ONETWOONE) then
- call one_two_one(field, &
- start_patch_i, end_patch_i, &
- start_patch_j, end_patch_j, &
- start_mem_i, end_mem_i, &
- start_mem_j, end_mem_j, &
- min_category, max_category, &
- smth_passes, msg_fill_val)
- else if (smth_opt == SMTHDESMTH) then
- call smth_desmth(field, &
- start_patch_i, end_patch_i, &
- start_patch_j, end_patch_j, &
- start_mem_i, end_mem_i, &
- start_mem_j, end_mem_j, &
- min_category, max_category, &
- smth_passes, msg_fill_val)
- else if (smth_opt == SMTHDESMTH_SPECIAL) then
- call smth_desmth_special(field, &
- start_patch_i, end_patch_i, &
- start_patch_j, end_patch_j, &
- start_mem_i, end_mem_i, &
- start_mem_j, end_mem_j, &
- min_category, max_category, &
- smth_passes, msg_fill_val)
- end if
- else if (grid_type == 'E') then
- if (smth_opt == ONETWOONE) then
- call one_two_one_egrid(field, &
- start_patch_i, end_patch_i, &
- start_patch_j, end_patch_j, &
- start_mem_i, end_mem_i, &
- start_mem_j, end_mem_j, &
- min_category, max_category, &
- smth_passes, msg_fill_val, 1.0)
- else if (smth_opt == SMTHDESMTH) then
- call smth_desmth_egrid(field, &
- start_patch_i, end_patch_i, &
- start_patch_j, end_patch_j, &
- start_mem_i, end_mem_i, &
- start_mem_j, end_mem_j, &
- min_category, max_category, &
- smth_passes, msg_fill_val, 1.0)
- else if (smth_opt == SMTHDESMTH_SPECIAL) then
- call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
- 'No smoothing will be done.')
- end if
- end if
- end if
-
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
- min_category, max_category, trim(landmask_name), &
- datestr, real_array=field)
- end if
-
- if (idomcatstatus == 0) then
- allocate(dominant_field(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
-
- if (.not. only_save_dominant) then
- field_count = field_count + 1
- call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
- i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
- end if
- do i=start_mem_i, end_mem_i
- do j=start_mem_j, end_mem_j
- if ((landmask(i,j) == 1. .and. is_water_mask) .or. &
- (landmask(i,j) == 0. .and. .not.is_water_mask)) then
- dominant = 0.
- dominant_field(i,j) = real(min_category-1)
- do k=min_category,max_category
- do kk=1,num_landmask_categories
- if (k == landmask_value(kk)) exit
- end do
- if (field(i,j,k) > dominant .and. kk > num_landmask_categories) then
- dominant_field(i,j) = real(k)
- dominant = field(i,j,k)
- end if
- end do
- else
- dominant = 0.
- dominant_field(i,j) = real(min_category-1)
- do k=min_category,max_category
- do kk=1,num_landmask_categories
- if (field(i,j,k) > dominant .and. k == landmask_value(kk)) then
- dominant_field(i,j) = real(k)
- dominant = field(i,j,k)
- end if
- end do
- end do
- end if
- end do
- end do
-
- call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, trim(domname), &
- datestr, dominant_field)
-
- deallocate(dominant_field)
- end if
-
- deallocate(field)
- end if
-
- !
- ! Now process all other fields specified by the user
- !
- call reset_next_field()
- ifieldstatus = 0
- do while (ifieldstatus == 0)
- call get_next_fieldname(fieldname, ifieldstatus)
-
- ! There is another field in the GEOGRID.TBL file
- if (ifieldstatus == 0) then
- temp_string(1:128) = fieldname
-
- ! If this field is still to be processed
- if (.not. hash_search(processed_fieldnames, temp_string)) then
-
- call hash_insert(processed_fieldnames, temp_string)
- call mprintf(.true.,STDOUT,' Processing %s', s1=trim(fieldname))
-
- call get_output_stagger(fieldname, istagger, istatus)
- dimnames(:) = 'null'
- subgrid_var=is_subgrid_var(fieldname)
- if( subgrid_var ) call get_subgrid_dim_name(dimnames)
-
- if (istagger == M .and. .not. subgrid_var ) then
- sm1 = start_mem_i
- em1 = end_mem_i
- sm2 = start_mem_j
- em2 = end_mem_j
- xlat_ptr => xlat_array
- xlon_ptr => xlon_array
- mapfac_ptr_x => mapfac_array_m_x
- mapfac_ptr_y => mapfac_array_m_y
- else if ( subgrid_var ) then
- sm1 = start_mem_sub_i
- em1 = end_mem_sub_i
- sm2 = start_mem_sub_j
- em2 = end_mem_sub_j
- xlat_ptr => xlat_array_subgrid
- xlon_ptr => xlon_array_subgrid
- mapfac_ptr_x => mapfac_array_x_subgrid
- mapfac_ptr_y => mapfac_array_y_subgrid
- else if (istagger == U) then ! In the case that extra_cols = .false.
- sm1 = start_mem_i ! we should have that end_mem_stag is
- em1 = end_mem_stag_i ! the same as end_mem, so we do not need
- sm2 = start_mem_j ! to check extra_cols or extra rows here
- em2 = end_mem_j
- xlat_ptr => xlat_array_u
- xlon_ptr => xlon_array_u
- mapfac_ptr_x => mapfac_array_u_x
- mapfac_ptr_y => mapfac_array_u_y
- else if (istagger == V) then
- sm1 = start_mem_i
- em1 = end_mem_i
- sm2 = start_mem_j
- em2 = end_mem_stag_j
- xlat_ptr => xlat_array_v
- xlon_ptr => xlon_array_v
- mapfac_ptr_x => mapfac_array_v_x
- mapfac_ptr_y => mapfac_array_v_y
- else if (istagger == HH) then ! E grid
- sm1 = start_mem_i
- em1 = end_mem_i
- sm2 = start_mem_j
- em2 = end_mem_j
- xlat_ptr => xlat_array
- xlon_ptr => xlon_array
- mapfac_ptr_x => mapfac_array_m_x
- mapfac_ptr_y => mapfac_array_m_y
- else if (istagger == VV) then ! E grid
- sm1 = start_mem_i
- em1 = end_mem_i
- sm2 = start_mem_j
- em2 = end_mem_stag_j
- xlat_ptr => xlat_array_v
- xlon_ptr => xlon_array_v
- mapfac_ptr_x => mapfac_array_v_x
- mapfac_ptr_y => mapfac_array_v_y
- end if
-
- call get_missing_fill_value(fieldname, msg_fill_val, istatus)
- if (istatus /= 0) msg_fill_val = NAN
-
- call get_halt_on_missing(fieldname, halt_on_missing, istatus)
- if (istatus /= 0) halt_on_missing = .false.
-
- ! Destination field type is CONTINUOUS
- if (iget_fieldtype(fieldname,istatus) == CONTINUOUS) then
- call get_max_levels(fieldname, min_level, max_level, istatus)
- allocate(field(sm1:em1, sm2:em2, min_level:max_level))
- field_count = field_count + 1
- call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
- i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=fieldname)
- if ( subgrid_var ) then
- call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
- sm1, em1, sm2, em2, min_level, max_level, &
- processed_domain, 1, sr_x=sr_x, sr_y=sr_y)
- else
- call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
- sm1, em1, sm2, em2, min_level, max_level, &
- processed_domain, 1, landmask=landmask)
- end if
-
- ! If user wants to halt when a missing value is found in output field, check now
- if (halt_on_missing) then
- do i=sm1, em1
- do j=sm2, em2
- ! Only need to examine k=1
- if (field(i,j,1) == msg_fill_val) then
- call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
- end if
- end do
- end do
- end if
-
- ! We may be asked to smooth the fractional field
- call get_smooth_option(fieldname, smth_opt, smth_passes, istatus)
- if (istatus == 0) then
- if (grid_type == 'C') then
- if (smth_opt == ONETWOONE) then
- call one_two_one(field, &
- start_patch_i, end_patch_i, &
- start_patch_j, end_patch_j, &
- sm1, em1, &
- sm2, em2, &
- min_level, max_level, &
- smth_passes, msg_fill_val)
- else if (smth_opt == SMTHDESMTH) then
- call smth_desmth(field, &
- start_patch_i, end_patch_i, &
- start_patch_j, end_patch_j, &
- sm1, em1, &
- sm2, em2, &
- min_level, max_level, &
- smth_passes, msg_fill_val)
- else if (smth_opt == SMTHDESMTH_SPECIAL) then
- call smth_desmth_special(field, &
- start_patch_i, end_patch_i, &
- start_patch_j, end_patch_j, &
- sm1, em1, &
- sm2, em2, &
- min_level, max_level, &
- smth_passes, msg_fill_val)
- end if
-
- else if (grid_type == 'E') then
-
- if (trim(fieldname) == 'HGT_M' ) then
- topo_flag_val=1.0
- mass_flag=1.0
- else if (trim(fieldname) == 'HGT_V') then
- topo_flag_val=1.0
- mass_flag=0.0
- else
- topo_flag_val=0.0
- end if
-
- if (smth_opt == ONETWOONE) then
- call one_two_one_egrid(field, &
- start_patch_i, end_patch_i, &
- start_patch_j, end_patch_j, &
- sm1, em1, &
- sm2, em2, &
- min_level, max_level, &
- smth_passes, topo_flag_val, mass_flag)
- else if (smth_opt == SMTHDESMTH) then
- call smth_desmth_egrid(field, &
- start_patch_i, end_patch_i, &
- start_patch_j, end_patch_j, &
- sm1, em1, &
- sm2, em2, &
- min_level, max_level, &
- smth_passes, topo_flag_val, mass_flag)
- else if (smth_opt == SMTHDESMTH_SPECIAL) then
- call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
- 'No smoothing will be done.')
- end if
-
- end if
- end if
- call write_field(sm1, em1, sm2, em2, &
- min_level, max_level, trim(fieldname), datestr, real_array=field)
-
- ! Do we calculate directional derivatives from this field?
- call get_dfdx_name(fieldname, gradname, istatus)
- if (istatus == 0) then
- allocate(slp_field(sm1:em1,sm2:em2,min_level:max_level))
- field_count = field_count + 1
- call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
- i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=gradname)
- if (grid_type == 'C') then
- call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, &
- which_domain, mapfac_ptr_x, subgrid_var)
- else if (grid_type == 'E') then
- call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, &
- which_domain, subgrid_var=subgrid_var)
- end if
- call write_field(sm1, em1, sm2, em2, &
- min_level, max_level, trim(gradname), datestr, real_array=slp_field)
- deallocate(slp_field)
- end if
- call get_dfdy_name(fieldname, gradname, istatus)
- if (istatus == 0) then
- allocate(slp_field(sm1:em1,sm2:em2,min_level:max_level))
- field_count = field_count + 1
- call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
- i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=gradname)
- if (grid_type == 'C') then
- call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, &
- which_domain, mapfac_ptr_y, subgrid_var)
- else if (grid_type == 'E') then
- call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, &
- which_domain, subgrid_var=subgrid_var)
- end if
- call write_field(sm1, em1, sm2, em2, &
- min_level, max_level, trim(gradname), datestr, real_array=slp_field)
- deallocate(slp_field)
- end if
-
- deallocate(field)
-
- ! Destination field type is CATEGORICAL
- else
- call get_max_categories(fieldname, min_category, max_category, istatus)
- allocate(field(sm1:em1, sm2:em2, min_category:max_category))
- ! Do we calculate a dominant category for this field?
- call get_domcategory_name(fieldname, domname, only_save_dominant, idomcatstatus)
-
- if (.not. only_save_dominant) then
- field_count = field_count + 1
- call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
- i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=fieldname)
- else
- field_count = field_count + 1
- call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
- i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
- end if
- if ( subgrid_var ) then
- call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
- sm1, em1, sm2, em2, min_category, max_category, &
- processed_domain, 1, sr_x=sr_x, sr_y=sr_y)
- else
- call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
- sm1, em1, sm2, em2, min_category, max_category, &
- processed_domain, 1, landmask=landmask)
- end if
-
- ! Find fractions
- do i=sm1, em1
- do j=sm2, em2
- sum = 0.
- do k=min_category,max_category
- sum = sum + field(i,j,k)
- end do
- if (sum > 0.0) then
- do k=min_category,max_category
- field(i,j,k) = field(i,j,k) / sum
- end do
- else
- do k=min_category,max_category
- field(i,j,k) = msg_fill_val
- end do
- end if
- end do
- end do
-
- ! If user wants to halt when a missing value is found in output field, check now
- if (halt_on_missing) then
- do i=sm1, em1
- do j=sm2, em2
- ! Only need to examine k=1
- if (field(i,j,1) == msg_fill_val) then
- call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
- end if
- end do
- end do
- end if
-
- ! If we should only save the dominant category, then no need to write out fractional field
- if (.not.only_save_dominant .or. (idomcatstatus /= 0)) then
-
- ! Finally, we may be asked to smooth the fractional field
- call get_smooth_option(fieldname, smth_opt, smth_passes, istatus)
- if (istatus == 0) then
- if (grid_type == 'C') then
- if (smth_opt == ONETWOONE) then
- call one_two_one(field, &
- start_patch_i, end_patch_i, &
- start_patch_j, end_patch_j, &
- sm1, em1, &
- sm2, em2, &
- min_category, max_category, &
- smth_passes, msg_fill_val)
- else if (smth_opt…
Large files files are truncated, but you can click here to view the full file