/WPS/geogrid/src/output_module.F

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MODULE OUTPUT_MODULE
!
! This module handles the output of the fields that are generated by the main
!   geogrid routines. This output may include output to a console and output to 
!   the WRF I/O API.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module output_module

   use parallel_module
   use gridinfo_module
   use misc_definitions_module
   use module_debug
#ifdef IO_BINARY
   use module_internal_header_util
#endif
 
   integer, parameter :: MAX_DIMENSIONS = 3
 
#ifdef _GEOGRID
   ! Information about fields that will be written
   integer :: NUM_AUTOMATIC_FIELDS    ! Set later, but very near to a parameter
#endif
 
   integer :: NUM_FIELDS
 
   type field_info
      integer :: ndims, istagger
      integer, dimension(MAX_DIMENSIONS) :: dom_start, mem_start, patch_start
      integer, dimension(MAX_DIMENSIONS) :: dom_end, mem_end, patch_end
      logical :: is_subgrid
      real, pointer, dimension(:,:,:) :: rdata_arr
  
      character (len=128), dimension(MAX_DIMENSIONS) :: dimnames
      character (len=128) :: fieldname, mem_order, stagger, units, descr
   end type field_info
 
   type (field_info), pointer, dimension(:) :: fields 
 
   ! WRF I/O API related variables
   integer :: handle 
 
   contains
 
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
   ! Name: output_init
   ! 
   ! Purpose: To initialize the output module. Such initialization may include 
   !   opening an X window, and making initialization calls to an I/O API.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
   subroutine output_init(nest_number, title, datestr, grid_type, dynopt, &
                          corner_lats, corner_lons, &
                          start_dom_1, end_dom_1, start_dom_2, end_dom_2, &
                          start_patch_1, end_patch_1, start_patch_2, end_patch_2, &
                          start_mem_1, end_mem_1, start_mem_2, end_mem_2, &
                          extra_col, extra_row)
 
#ifdef _GEOGRID
      use llxy_module
      use source_data_module
#endif
  
      implicit none
  
      ! Arguments
      integer, intent(in) :: nest_number, dynopt, &
                             start_dom_1, end_dom_1, start_dom_2, end_dom_2, &
                             start_patch_1, end_patch_1, start_patch_2, end_patch_2, &
                             start_mem_1, end_mem_1, start_mem_2, end_mem_2
      real, dimension(16), intent(in) :: corner_lats, corner_lons
      logical, intent(in) :: extra_col, extra_row
      character (len=1), intent(in) :: grid_type
      character (len=19), intent(in) :: datestr
      character (len=*), intent(in) :: title
  
#include "wrf_io_flags.h"
#include "wrf_status_codes.h"
  
      ! Local variables
      integer :: i, istatus, save_domain, comm_1, comm_2
      integer :: sp1, ep1, sp2, ep2, ep1_stag, ep2_stag
      integer :: ngeo_flags
      integer :: num_land_cat, min_land_cat, max_land_cat
      real :: dx, dy, cen_lat, cen_lon, moad_cen_lat
      character (len=128) :: coption, temp_fldname
      character (len=128), dimension(1) :: geo_flags
      character (len=MAX_FILENAME_LEN) :: output_fname
      logical :: supports_training, supports_3d_fields
  
      call init_output_fields(nest_number, grid_type, &
                              start_dom_1, end_dom_1, start_dom_2, end_dom_2, &
                              start_patch_1, end_patch_1, start_patch_2, end_patch_2, &
                              start_mem_1, end_mem_1, start_mem_2, end_mem_2, &
                              extra_col, extra_row)
  
      if (my_proc_id == IO_NODE .or. do_tiled_output) then
      istatus = 0
#ifdef IO_BINARY
      if (io_form_output == BINARY) call ext_int_ioinit('sysdep info', istatus)
#endif
#ifdef IO_NETCDF
      if (io_form_output == NETCDF) call ext_ncd_ioinit('sysdep info', istatus)
#endif
#ifdef IO_GRIB1
      if (io_form_output == GRIB1) call ext_gr1_ioinit('sysdep info', istatus)
#endif
      call mprintf((istatus /= 0), ERROR, 'Error in ext_pkg_ioinit')
      
      ! Find out what this implementation of WRF I/O API supports
      istatus = 0
#ifdef IO_BINARY
      if (io_form_output == BINARY) coption = 'REQUIRE'
#endif
#ifdef IO_NETCDF
      if (io_form_output == NETCDF) call ext_ncd_inquiry('OPEN_COMMIT_WRITE',coption,istatus)
#endif
#ifdef IO_GRIB1
      if (io_form_output == GRIB1) call ext_gr1_inquiry('OPEN_COMMIT_WRITE',coption,istatus)
#endif
      call mprintf((istatus /= 0), ERROR, 'Error in ext_pkg_inquiry')
  
      if (index(coption,'ALLOW') /= 0) then
         supports_training = .false.
      else if (index(coption,'REQUIRE') /= 0) then
         supports_training = .true.
      else if (index(coption,'NO') /= 0) then
         supports_training = .false.
      end if 
  
      istatus = 0
#ifdef IO_BINARY
      if (io_form_output == BINARY) coption = 'YES'
#endif
#ifdef IO_NETCDF
      if (io_form_output == NETCDF) call ext_ncd_inquiry('SUPPORT_3D_FIELDS',coption,istatus)
#endif
#ifdef IO_GRIB1
      if (io_form_output == GRIB1) call ext_gr1_inquiry('SUPPORT_3D_FIELDS',coption,istatus)
#endif
      call mprintf((istatus /= 0), ERROR, 'Error in ext_pkg_inquiry')
  
      if (index(coption,'YES') /= 0) then
         supports_3d_fields = .true.
      else if (index(coption,'NO') /= 0) then
         supports_3d_fields = .false.
! BUG: What if we have no plans to write 3-d fields? We should take this into account.
         call mprintf(.true.,ERROR,'WRF I/O API implementation does NOT support 3-d fields.')
      end if
  
      comm_1 = 1
      comm_2 = 1
  
#ifdef _GEOGRID
      output_fname = ' '
      if (grid_type == 'C') then
#ifdef IO_BINARY
         if (io_form_output == BINARY) output_fname = trim(opt_output_from_geogrid_path)//'geo_em.d  .int'
#endif
#ifdef IO_NETCDF
         if (io_form_output == NETCDF) output_fname = trim(opt_output_from_geogrid_path)//'geo_em.d  .nc'
#endif
#ifdef IO_GRIB1
         if (io_form_output == GRIB1) output_fname = trim(opt_output_from_geogrid_path)//'geo_em.d  .grib'
#endif
         i = len_trim(opt_output_from_geogrid_path)
         write(output_fname(i+9:i+10),'(i2.2)') nest_number
      else if (grid_type == 'E') then
         if (nest_number == 1) then
#ifdef IO_BINARY
            if (io_form_output == BINARY) output_fname = trim(opt_output_from_geogrid_path)//'geo_nmm.d  .int'
#endif
#ifdef IO_NETCDF
            if (io_form_output == NETCDF) output_fname = trim(opt_output_from_geogrid_path)//'geo_nmm.d  .nc'
#endif
#ifdef IO_GRIB1
            if (io_form_output == GRIB1) output_fname = trim(opt_output_from_geogrid_path)//'geo_nmm.d  .grib'
#endif
            i = len_trim(opt_output_from_geogrid_path)
            write(output_fname(i+10:i+11),'(i2.2)') nest_number
         else
#ifdef IO_BINARY
            if (io_form_output == BINARY) output_fname = trim(opt_output_from_geogrid_path)//'geo_nmm_nest.l  .int'
#endif
#ifdef IO_NETCDF
            if (io_form_output == NETCDF) output_fname = trim(opt_output_from_geogrid_path)//'geo_nmm_nest.l  .nc'
#endif
#ifdef IO_GRIB1
            if (io_form_output == GRIB1) output_fname = trim(opt_output_from_geogrid_path)//'geo_nmm_nest.l  .grib'
#endif
            i = len_trim(opt_output_from_geogrid_path)
            write(output_fname(i+15:i+16),'(i2.2)') nest_number-1
         end if
      end if

      if (nprocs > 1 .and. do_tiled_output) then
         write(output_fname(len_trim(output_fname)+1:len_trim(output_fname)+5), '(a1,i4.4)') &
              '_', my_proc_id 
      end if
#endif
  
#ifdef _METGRID
      output_fname = ' '
      if (grid_type == 'C') then
#ifdef IO_BINARY
         if (io_form_output == BINARY) then
            output_fname = trim(opt_output_from_metgrid_path)//'met_em.d  .'//trim(datestr)//'.int'
         end if
#endif
#ifdef IO_NETCDF
         if (io_form_output == NETCDF) then
            output_fname = trim(opt_output_from_metgrid_path)//'met_em.d  .'//trim(datestr)//'.nc'
         end if
#endif
#ifdef IO_GRIB1
         if (io_form_output == GRIB1) then
            output_fname = trim(opt_output_from_metgrid_path)//'met_em.d  .'//trim(datestr)//'.grib'
         end if
#endif
         i = len_trim(opt_output_from_metgrid_path)
         write(output_fname(i+9:i+10),'(i2.2)') nest_number
      else if (grid_type == 'E') then
#ifdef IO_BINARY
         if (io_form_output == BINARY) then
            output_fname = trim(opt_output_from_metgrid_path)//'met_nmm.d  .'//trim(datestr)//'.int'
         end if
#endif
#ifdef IO_NETCDF
         if (io_form_output == NETCDF) then
            output_fname = trim(opt_output_from_metgrid_path)//'met_nmm.d  .'//trim(datestr)//'.nc'
         end if
#endif
#ifdef IO_GRIB1
         if (io_form_output == GRIB1) then
            output_fname = trim(opt_output_from_metgrid_path)//'met_nmm.d  .'//trim(datestr)//'.grib'
         end if
#endif
         i = len_trim(opt_output_from_metgrid_path)
         write(output_fname(i+10:i+11),'(i2.2)') nest_number
      end if

      if (nprocs > 1 .and. do_tiled_output) then
         write(output_fname(len_trim(output_fname)+1:len_trim(output_fname)+5), '(a1,i4.4)') &
              '_', my_proc_id 
      end if
#endif
      end if
  
      call parallel_bcast_logical(supports_training)
  
      ! If the implementation requires or supports open_for_write begin/commit semantics
      if (supports_training) then
  
         if (my_proc_id == IO_NODE .or. do_tiled_output) then
            istatus = 0
#ifdef IO_BINARY
            if (io_form_output == BINARY) then
               call ext_int_open_for_write_begin(trim(output_fname), comm_1, comm_2, 'sysdep info', handle, istatus)
            end if
#endif
#ifdef IO_NETCDF
            if (io_form_output == NETCDF) then
               call ext_ncd_open_for_write_begin(trim(output_fname), comm_1, comm_2, 'sysdep info', handle, istatus)
            end if
#endif
#ifdef IO_GRIB1
            if (io_form_output == GRIB1) then
               call ext_gr1_open_for_write_begin(trim(output_fname), comm_1, comm_2, 'sysdep info', handle, istatus)
            end if
#endif
            call mprintf((istatus /= 0), ERROR, 'Error in ext_pkg_open_for_write_begin.')
         end if
   
         do i=1,NUM_FIELDS
   
            allocate(fields(i)%rdata_arr(fields(i)%mem_start(1):fields(i)%mem_end(1), &
                                         fields(i)%mem_start(2):fields(i)%mem_end(2), &
                                         fields(i)%mem_start(3):fields(i)%mem_end(3)))
     
            call write_field(fields(i)%mem_start(1), fields(i)%mem_end(1), fields(i)%mem_start(2), &
                             fields(i)%mem_end(2), fields(i)%mem_start(3), fields(i)%mem_end(3), &
                             trim(fields(i)%fieldname), datestr, fields(i)%rdata_arr, is_training=.true.)
     
            deallocate(fields(i)%rdata_arr)
    
         end do
   
         if (my_proc_id == IO_NODE .or. do_tiled_output) then
            istatus = 0
#ifdef IO_BINARY
            if (io_form_output == BINARY) call ext_int_open_for_write_commit(handle, istatus)
#endif
#ifdef IO_NETCDF
            if (io_form_output == NETCDF) call ext_ncd_open_for_write_commit(handle, istatus)
#endif
#ifdef IO_GRIB1
            if (io_form_output == GRIB1) call ext_gr1_open_for_write_commit(handle, istatus)
#endif
            call mprintf((istatus /= 0), ERROR, 'Error in ext_pkg_write_commit')
         end if
  
      else ! No training required
  
         if (my_proc_id == IO_NODE .or. do_tiled_output) then
            istatus = 0
#ifdef IO_BINARY
            if (io_form_output == BINARY) then
               call ext_int_open_for_write(trim(output_fname), comm_1, comm_2, 'sysdep info', handle, istatus)
            end if
#endif
#ifdef IO_NETCDF
            if (io_form_output == NETCDF) then
               call ext_ncd_open_for_write(trim(output_fname), comm_1, comm_2, 'sysdep info', handle, istatus)
            end if
#endif
#ifdef IO_GRIB1
            if (io_form_output == GRIB1) then
               call mprintf(.true.,ERROR,'In output_init(), GRIB1 requires begin/commit open sequence.')
            end if
#endif
            call mprintf((istatus /= 0),ERROR,'Error in ext_pkg_open_for_write_begin')
         end if
  
      end if 
  
#ifdef _GEOGRID
      sp1 = start_patch_1
      ep1 = end_patch_1
      sp2 = start_patch_2
      ep2 = end_patch_2
  
      if (grid_type == 'C') then
         if (extra_col .or. (my_proc_id == IO_NODE .and. .not. do_tiled_output)) then
            ep1_stag = ep1 + 1
         else
            ep1_stag = ep1
         end if
         if (extra_row .or. (my_proc_id == IO_NODE .and. .not. do_tiled_output)) then
            ep2_stag = ep2 + 1
         else
            ep2_stag = ep2
         end if
      else if (grid_type == 'E') then
         ep1 = ep1
         ep2 = ep2
         ep1_stag = ep1
         ep2_stag = ep2
      end if

      if (grid_type == 'C') then
         dx = proj_stack(nest_number)%dx
         dy = proj_stack(nest_number)%dy

         save_domain = iget_selected_domain()

         ! Note: In the following, we use ixdim/2 rather than (ixdim+1)/2 because
         !       the i/j coordinates given to xytoll must be with respect to the
         !       mass grid, and ixdim and jydim are the full grid dimensions.

         ! Get MOAD_CEN_LAT
         call select_domain(1)
         call xytoll(real(ixdim(1))/2.,real(jydim(1))/2., moad_cen_lat, cen_lon, M)

         ! Get CEN_LAT and CEN_LON for this nest
         call select_domain(nest_number)
         call xytoll(real(ixdim(nest_number))/2.,real(jydim(nest_number))/2., cen_lat, cen_lon, M)

         call select_domain(save_domain)

         ngeo_flags = 1
         geo_flags(1) = 'FLAG_MF_XY' 
      else if (grid_type == 'E') then
         dx = dxkm / 3**(nest_number-1)   ! For NMM, nest_number is really nesting level
         dy = dykm / 3**(nest_number-1)
         moad_cen_lat = 0.
         cen_lat=known_lat
         cen_lon=known_lon

         ngeo_flags = 0
      end if

      write(temp_fldname,'(a)') 'LANDUSEF'
      call get_max_categories(temp_fldname, min_land_cat, max_land_cat, istatus)
      num_land_cat = max_land_cat - min_land_cat + 1
  
      ! We may now write global attributes to the file
      call write_global_attrs(title, datestr, grid_type, dynopt, ixdim(nest_number), jydim(nest_number), &
                              0, sp1, ep1, sp1, ep1_stag, sp2, ep2, sp2, ep2_stag,                       &
                              iproj_type, source_mminlu, num_land_cat, source_iswater, source_islake,    &
                              source_isice, source_isurban, source_isoilwater, nest_number,              &
                              parent_id(nest_number),                                                    &
                              nint(parent_ll_x(nest_number)), nint(parent_ll_y(nest_number)),            &
                              nint(parent_ur_x(nest_number)), nint(parent_ur_y(nest_number)),            &
                              dx, dy, cen_lat, moad_cen_lat,                                             &
                              cen_lon, stand_lon, truelat1, truelat2, pole_lat, pole_lon,                &
                              parent_grid_ratio(nest_number),                                            &
                              subgrid_ratio_x(nest_number), subgrid_ratio_y(nest_number),                &
                              corner_lats, corner_lons, flags=geo_flags, nflags=ngeo_flags)
#endif
 
   end subroutine output_init
 
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
   ! Name: init_output_fields
   ! 
   ! Purpose: To fill in structures describing each of the fields that will be 
   !   written to the I/O API 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
   subroutine init_output_fields(nest_num, grid_type, &
                                 start_dom_1, end_dom_1, start_dom_2, end_dom_2, &
                                 start_patch_1, end_patch_1, start_patch_2, end_patch_2, &
                                 start_mem_1, end_mem_1, start_mem_2, end_mem_2, &
                                 extra_col, extra_row)
 
 
      ! Modules
#ifdef _GEOGRID
      use source_data_module
#endif
#ifdef _METGRID
      use storage_module
#endif
      use parallel_module
      use gridinfo_module, only : subgrid_ratio_x, subgrid_ratio_y
  
      implicit none
  
      ! Arguments
      integer, intent(in) :: nest_num
      integer, intent(in) :: start_dom_1, end_dom_1, start_dom_2, end_dom_2, &
                             start_patch_1, end_patch_1, start_patch_2, end_patch_2, &
                             start_mem_1, end_mem_1, start_mem_2, end_mem_2
      logical, intent(in) :: extra_col, extra_row
      character (len=1), intent(in) :: grid_type
  
#include "wrf_io_flags.h"
#include "wrf_status_codes.h"
  
      ! Local variables
      integer :: i, istagger, ifieldstatus, &
                 nfields, min_category, max_category
      integer :: lhalo_width, rhalo_width, bhalo_width, thalo_width
      integer :: ndims
      character (len=128) :: fieldname
      character (len=128) :: memorder, units, description
      character (len=128), dimension(3) :: dimnames 
      integer :: sr_x, sr_y, istatus
      logical :: subgrid_var
 
      !
      ! Get local nest subgrid refinement ratio
      !
      sr_x=subgrid_ratio_x(nest_num)
      sr_y=subgrid_ratio_y(nest_num)

      !
      ! First find out how many fields there are
      !
      call reset_next_field()
  
      ifieldstatus = 0
      nfields = 0
      do while (ifieldstatus == 0)
         call get_next_output_fieldname(nest_num, fieldname, ndims, &
                                        min_category, max_category, &
                                        istagger, memorder, dimnames, &
                                        units, description, subgrid_var, ifieldstatus)
   
         if (ifieldstatus == 0) nfields = nfields + 1
      end do
  
#ifdef _METGRID
      NUM_FIELDS = nfields
#endif
  
#ifdef _GEOGRID
      if (grid_type == 'C') NUM_AUTOMATIC_FIELDS = 24
      if (grid_type == 'E') NUM_AUTOMATIC_FIELDS = 7
  
      NUM_FIELDS = nfields+NUM_AUTOMATIC_FIELDS
      allocate(fields(NUM_FIELDS))
  
      !
      ! There are some fields that will always be computed
      !   Initialize those fields first, followed by all user-specified fields
      !
      if (grid_type == 'C') then
         fields(1)%fieldname = 'XLAT_M'
         fields(1)%units = 'degrees latitude'
         fields(1)%descr = 'Latitude on mass grid'
   
         fields(2)%fieldname = 'XLONG_M'
         fields(2)%units = 'degrees longitude'
         fields(2)%descr = 'Longitude on mass grid'
   
         fields(3)%fieldname = 'XLAT_V'
         fields(3)%units = 'degrees latitude'
         fields(3)%descr = 'Latitude on V grid'
   
         fields(4)%fieldname = 'XLONG_V'
         fields(4)%units = 'degrees longitude'
         fields(4)%descr = 'Longitude on V grid'
   
         fields(5)%fieldname = 'XLAT_U'
         fields(5)%units = 'degrees latitude'
         fields(5)%descr = 'Latitude on U grid'
   
         fields(6)%fieldname = 'XLONG_U'
         fields(6)%units = 'degrees longitude'
         fields(6)%descr = 'Longitude on U grid'
  
         fields(7)%fieldname = 'CLAT'
         fields(7)%units = 'degrees latitude'
         fields(7)%descr = 'Computational latitude on mass grid'
   
         fields(8)%fieldname = 'CLONG'
         fields(8)%units = 'degrees longitude'
         fields(8)%descr = 'Computational longitude on mass grid'

         fields(9)%fieldname = 'MAPFAC_M'
         fields(9)%units = 'none'
         fields(9)%descr = 'Mapfactor on mass grid'

         fields(10)%fieldname = 'MAPFAC_V'
         fields(10)%units = 'none'
         fields(10)%descr = 'Mapfactor on V grid'

         fields(11)%fieldname = 'MAPFAC_U'
         fields(11)%units = 'none'
         fields(11)%descr = 'Mapfactor on U grid'
   
         fields(12)%fieldname = 'MAPFAC_MX'
         fields(12)%units = 'none'
         fields(12)%descr = 'Mapfactor (x-dir) on mass grid'
   
         fields(13)%fieldname = 'MAPFAC_VX'
         fields(13)%units = 'none'
         fields(13)%descr = 'Mapfactor (x-dir) on V grid'
   
         fields(14)%fieldname = 'MAPFAC_UX'
         fields(14)%units = 'none'
         fields(14)%descr = 'Mapfactor (x-dir) on U grid'

         fields(15)%fieldname = 'MAPFAC_MY'
         fields(15)%units = 'none'
         fields(15)%descr = 'Mapfactor (y-dir) on mass grid'
   
         fields(16)%fieldname = 'MAPFAC_VY'
         fields(16)%units = 'none'
         fields(16)%descr = 'Mapfactor (y-dir) on V grid'
   
         fields(17)%fieldname = 'MAPFAC_UY'
         fields(17)%units = 'none'
         fields(17)%descr = 'Mapfactor (y-dir) on U grid'
   
         fields(18)%fieldname = 'E'
         fields(18)%units = '-'
         fields(18)%descr = 'Coriolis E parameter'
   
         fields(19)%fieldname = 'F'
         fields(19)%units = '-'
         fields(19)%descr = 'Coriolis F parameter'
   
         fields(20)%fieldname = 'SINALPHA'
         fields(20)%units = 'none'
         fields(20)%descr = 'Sine of rotation angle'
   
         fields(21)%fieldname = 'COSALPHA'
         fields(21)%units = 'none'
         fields(21)%descr = 'Cosine of rotation angle'
   
         fields(22)%fieldname = 'LANDMASK'
         fields(22)%units = 'none'
         fields(22)%descr = 'Landmask : 1=land, 0=water'

         fields(23)%fieldname = 'FXLONG'
         fields(23)%units = 'degrees longitude'
         fields(23)%descr = 'Longitude on refined grid'

         fields(24)%fieldname = 'FXLAT'
         fields(24)%units = 'degrees latitude'
         fields(24)%descr = 'Latitude on refined grid'

      else if (grid_type == 'E') then
         fields(1)%fieldname = 'XLAT_M'
         fields(1)%units = 'degrees latitude'
         fields(1)%descr = 'Latitude on mass grid'
   
         fields(2)%fieldname = 'XLONG_M'
         fields(2)%units = 'degrees longitude'
         fields(2)%descr = 'Longitude on mass grid'
   
         fields(3)%fieldname = 'XLAT_V'
         fields(3)%units = 'degrees latitude'
         fields(3)%descr = 'Latitude on velocity grid'
   
         fields(4)%fieldname = 'XLONG_V'
         fields(4)%units = 'degrees longitude'
         fields(4)%descr = 'Longitude on velocity grid'
   
         fields(5)%fieldname = 'E'
         fields(5)%units = '-'
         fields(5)%descr = 'Coriolis E parameter'
   
         fields(6)%fieldname = 'F'
         fields(6)%units = '-'
         fields(6)%descr = 'Coriolis F parameter'
   
         fields(7)%fieldname = 'LANDMASK'
         fields(7)%units = 'none'
         fields(7)%descr = 'Landmask : 1=land, 0=water'
  
      end if
  
      !
      ! General defaults for "always computed" fields 
      !
      do i=1,NUM_AUTOMATIC_FIELDS
         fields(i)%ndims = 2
         fields(i)%dom_start(1) = start_dom_1 
         fields(i)%dom_start(2) = start_dom_2
         fields(i)%dom_start(3) = 1
         fields(i)%mem_start(1) = start_mem_1 
         fields(i)%mem_start(2) = start_mem_2
         fields(i)%mem_start(3) = 1
         fields(i)%patch_start(1) = start_patch_1
         fields(i)%patch_start(2) = start_patch_2
         fields(i)%patch_start(3) = 1
         fields(i)%dom_end(1) = end_dom_1
         fields(i)%dom_end(2) = end_dom_2
         fields(i)%dom_end(3) = 1
         fields(i)%mem_end(1) = end_mem_1
         fields(i)%mem_end(2) = end_mem_2
         fields(i)%mem_end(3) = 1
         fields(i)%patch_end(1) = end_patch_1
         fields(i)%patch_end(2) = end_patch_2
         fields(i)%patch_end(3) = 1
         fields(i)%dimnames(3) = ' '
         fields(i)%mem_order = 'XY'
         fields(i)%stagger = 'M'
         fields(i)%is_subgrid = .false.
         if (grid_type == 'C') then
            fields(i)%istagger = M
         else if (grid_type == 'E') then
            fields(i)%istagger = HH
         end if
         fields(i)%dimnames(1) = 'west_east'
         fields(i)%dimnames(2) = 'south_north'
      end do
  
      !
      ! Perform adjustments to metadata for non-mass-staggered "always computed" fields
      !
      if (grid_type == 'C') then
         ! Lat V
         if (extra_row) then
            fields(3)%dom_end(2) = fields(3)%dom_end(2) + 1
            fields(3)%mem_end(2) = fields(3)%mem_end(2) + 1
            fields(3)%patch_end(2) = fields(3)%patch_end(2) + 1
         else if (my_proc_id == IO_NODE .and. .not. do_tiled_output) then
            fields(3)%dom_end(2) = fields(3)%dom_end(2) + 1
         end if
         fields(3)%dimnames(2) = 'south_north_stag'
         fields(3)%stagger = 'V'
         fields(3)%istagger = V
   
         ! Lon V
         if (extra_row) then
            fields(4)%dom_end(2) = fields(4)%dom_end(2) + 1
            fields(4)%mem_end(2) = fields(4)%mem_end(2) + 1
            fields(4)%patch_end(2) = fields(4)%patch_end(2) + 1
         else if (my_proc_id == IO_NODE .and. .not. do_tiled_output) then
            fields(4)%dom_end(2) = fields(4)%dom_end(2) + 1
         end if
         fields(4)%dimnames(2) = 'south_north_stag'
         fields(4)%stagger = 'V'
         fields(4)%istagger = V
   
         ! Lat U
         if (extra_col) then
            fields(5)%dom_end(1) = fields(5)%dom_end(1) + 1
            fields(5)%mem_end(1) = fields(5)%mem_end(1) + 1
            fields(5)%patch_end(1) = fields(5)%patch_end(1) + 1
         else if (my_proc_id == IO_NODE .and. .not. do_tiled_output) then
            fields(5)%dom_end(1) = fields(5)%dom_end(1) + 1
         end if
         fields(5)%dimnames(1) = 'west_east_stag'
         fields(5)%stagger = 'U'
         fields(5)%istagger = U
   
         ! Lon U
         if (extra_col) then
            fields(6)%dom_end(1) = fields(6)%dom_end(1) + 1
            fields(6)%mem_end(1) = fields(6)%mem_end(1) + 1
            fields(6)%patch_end(1) = fields(6)%patch_end(1) + 1
         else if (my_proc_id == IO_NODE .and. .not. do_tiled_output) then
            fields(6)%dom_end(1) = fields(6)%dom_end(1) + 1
         end if
         fields(6)%dimnames(1) = 'west_east_stag'
         fields(6)%stagger = 'U'
         fields(6)%istagger = U

         ! Mapfac V
         if (extra_row) then
            fields(10)%dom_end(2) = fields(10)%dom_end(2) + 1
            fields(10)%mem_end(2) = fields(10)%mem_end(2) + 1
            fields(10)%patch_end(2) = fields(10)%patch_end(2) + 1
         else if (my_proc_id == IO_NODE .and. .not. do_tiled_output) then
            fields(10)%dom_end(2) = fields(10)%dom_end(2) + 1
         end if
         fields(10)%dimnames(2) = 'south_north_stag'
         fields(10)%stagger = 'V'
         fields(10)%istagger = V

         ! Mapfac U
         if (extra_col) then
            fields(11)%dom_end(1) = fields(11)%dom_end(1) + 1
            fields(11)%mem_end(1) = fields(11)%mem_end(1) + 1
            fields(11)%patch_end(1) = fields(11)%patch_end(1) + 1
         else if (my_proc_id == IO_NODE .and. .not. do_tiled_output) then
            fields(11)%dom_end(1) = fields(11)%dom_end(1) + 1
         end if
         fields(11)%dimnames(1) = 'west_east_stag'
         fields(11)%stagger = 'U' 
         fields(11)%istagger = U
   
         ! Mapfac V-X
         if (extra_row) then
            fields(13)%dom_end(2) = fields(13)%dom_end(2) + 1
            fields(13)%mem_end(2) = fields(13)%mem_end(2) + 1
            fields(13)%patch_end(2) = fields(13)%patch_end(2) + 1
         else if (my_proc_id == IO_NODE .and. .not. do_tiled_output) then
            fields(13)%dom_end(2) = fields(13)%dom_end(2) + 1
         end if
         fields(13)%dimnames(2) = 'south_north_stag'
         fields(13)%stagger = 'V'
         fields(13)%istagger = V
   
         ! Mapfac U-X
         if (extra_col) then
            fields(14)%dom_end(1) = fields(14)%dom_end(1) + 1
            fields(14)%mem_end(1) = fields(14)%mem_end(1) + 1
            fields(14)%patch_end(1) = fields(14)%patch_end(1) + 1
         else if (my_proc_id == IO_NODE .and. .not. do_tiled_output) then
            fields(14)%dom_end(1) = fields(14)%dom_end(1) + 1
         end if
         fields(14)%dimnames(1) = 'west_east_stag'
         fields(14)%stagger = 'U'
         fields(14)%istagger = U

         ! Mapfac V-Y
         if (extra_row) then
            fields(16)%dom_end(2) = fields(16)%dom_end(2) + 1
            fields(16)%mem_end(2) = fields(16)%mem_end(2) + 1
            fields(16)%patch_end(2) = fields(16)%patch_end(2) + 1
         else if (my_proc_id == IO_NODE .and. .not. do_tiled_output) then
            fields(16)%dom_end(2) = fields(16)%dom_end(2) + 1
         end if
         fields(16)%dimnames(2) = 'south_north_stag'
         fields(16)%stagger = 'V'
         fields(16)%istagger = V

         ! Mapfac U-Y
         if (extra_col) then
            fields(17)%dom_end(1) = fields(17)%dom_end(1) + 1
            fields(17)%mem_end(1) = fields(17)%mem_end(1) + 1
            fields(17)%patch_end(1) = fields(17)%patch_end(1) + 1
         else if (my_proc_id == IO_NODE .and. .not. do_tiled_output) then
            fields(17)%dom_end(1) = fields(17)%dom_end(1) + 1
         end if
         fields(17)%dimnames(1) = 'west_east_stag'
         fields(17)%stagger = 'U'
         fields(17)%istagger = U
  
         !
         ! Perform adjustments for subgrid lat/lon fields
         !
         do i=23,24
            fields(i)%is_subgrid=.true.
            call get_subgrid_dim_name(fields(i)%dimnames(1:2))
            fields(i)%mem_start(1)   = (fields(i)%mem_start(1) - 1)*sr_x + 1
            fields(i)%patch_start(1) = (fields(i)%patch_start(1) - 1)*sr_x + 1
            fields(i)%dom_start(1)   = (fields(i)%dom_start(1) - 1)*sr_x + 1

            fields(i)%mem_end(1)     = (fields(i)%mem_end(1) + 1)*sr_x
            fields(i)%patch_end(1)   = (fields(i)%patch_end(1) + 1)*sr_x
            fields(i)%dom_end(1)     = (fields(i)%dom_end(1) + 1)*sr_x

            fields(i)%mem_start(2)   = (fields(i)%mem_start(2) - 1)*sr_y + 1
            fields(i)%patch_start(2) = (fields(i)%patch_start(2) - 1)*sr_y + 1
            fields(i)%dom_start(2)   = (fields(i)%dom_start(2) - 1)*sr_y + 1

            fields(i)%mem_end(2)     = (fields(i)%mem_end(2) + 1)*sr_y
            fields(i)%patch_end(2)   = (fields(i)%patch_end(2) + 1)*sr_y
            fields(i)%dom_end(2)     = (fields(i)%dom_end(2) + 1)*sr_y
         enddo
 
      else if (grid_type == 'E') then
         ! Lat V
         fields(3)%stagger = 'V'
         fields(3)%istagger = VV
   
         ! Lon V
         fields(4)%stagger = 'V'
         fields(4)%istagger = VV
   
      end if
#endif
  
      !
      ! Now set up the field_info structure for each user-specified field
      !
      call reset_next_field()
  
      ifieldstatus = 0
#ifdef _GEOGRID
      nfields = NUM_AUTOMATIC_FIELDS+1
#endif
#ifdef _METGRID
      allocate(fields(NUM_FIELDS))
      nfields = 1
#endif
  
      do while (ifieldstatus == 0)  !{
         call get_next_output_fieldname(nest_num, fieldname, ndims, &
                                      min_category, max_category, &
                                      istagger, memorder, dimnames, &
                                      units, description, subgrid_var, ifieldstatus)
   
         if (ifieldstatus == 0) then !{
     
            fields(nfields)%ndims = ndims
            fields(nfields)%fieldname = fieldname
            fields(nfields)%istagger = istagger
            if (istagger == M) then
               fields(nfields)%stagger = 'M'
            else if (istagger == U) then
               fields(nfields)%stagger = 'U'
            else if (istagger == V) then
               fields(nfields)%stagger = 'V'
            else if (istagger == HH) then
               fields(nfields)%stagger = 'M'
            else if (istagger == VV) then
               fields(nfields)%stagger = 'V'
            end if
            fields(nfields)%mem_order = memorder
            fields(nfields)%dimnames(1) = dimnames(1)
            fields(nfields)%dimnames(2) = dimnames(2)
            fields(nfields)%dimnames(3) = dimnames(3)
            fields(nfields)%units = units
            fields(nfields)%descr = description
    
            fields(nfields)%dom_start(1)   = start_dom_1
            fields(nfields)%dom_start(2)   = start_dom_2
            fields(nfields)%dom_start(3)   = min_category
            fields(nfields)%mem_start(1)   = start_mem_1
            fields(nfields)%mem_start(2)   = start_mem_2
            fields(nfields)%mem_start(3)   = min_category
            fields(nfields)%patch_start(1) = start_patch_1
            fields(nfields)%patch_start(2) = start_patch_2
            fields(nfields)%patch_start(3) = min_category
    
            fields(nfields)%dom_end(1)   = end_dom_1
            fields(nfields)%dom_end(2)   = end_dom_2
            fields(nfields)%dom_end(3)   = max_category
            fields(nfields)%mem_end(1)   = end_mem_1
            fields(nfields)%mem_end(2)   = end_mem_2
            fields(nfields)%mem_end(3)   = max_category
            fields(nfields)%patch_end(1) = end_patch_1
            fields(nfields)%patch_end(2) = end_patch_2
            fields(nfields)%patch_end(3) = max_category

            fields(nfields)%is_subgrid=subgrid_var
    
            if (extra_col .and. (istagger == U .or. subgrid_var)) then !{
               fields(nfields)%dom_end(1)   = fields(nfields)%dom_end(1) + 1
               fields(nfields)%mem_end(1)   = fields(nfields)%mem_end(1) + 1
               fields(nfields)%patch_end(1) = fields(nfields)%patch_end(1) + 1
            else if ((istagger == U .or. subgrid_var) .and. my_proc_id == IO_NODE .and. .not. do_tiled_output) then
               fields(nfields)%dom_end(1)=fields(nfields)%dom_end(1) + 1
            end if !}
    
            if (extra_row .and. (istagger == V .or. subgrid_var)) then !{
               fields(nfields)%dom_end(2)   = fields(nfields)%dom_end(2) + 1
               fields(nfields)%mem_end(2)   = fields(nfields)%mem_end(2) + 1
               fields(nfields)%patch_end(2) = fields(nfields)%patch_end(2) + 1
            else if ((istagger == V .or. subgrid_var) .and. my_proc_id == IO_NODE .and. .not. do_tiled_output) then
               fields(nfields)%dom_end(2)=fields(nfields)%dom_end(2) + 1
            end if !}

#ifdef _METGRID
            lhalo_width = start_patch_1 - start_mem_1      ! Halo width on left of patch
            rhalo_width = end_mem_1     - end_patch_1      ! Halo width on right of patch
            bhalo_width = start_patch_2 - start_mem_2      ! Halo width on bottom of patch
            thalo_width = end_mem_2     - end_patch_2      ! Halo width on top of patch
#else
            lhalo_width = 0
            rhalo_width = 0
            bhalo_width = 0
            thalo_width = 0
#endif

            if (subgrid_var) then
               fields(nfields)%mem_start(1)   = (fields(nfields)%mem_start(1) + lhalo_width - 1)*sr_x + 1 - lhalo_width
               fields(nfields)%patch_start(1) = (fields(nfields)%patch_start(1)             - 1)*sr_x + 1
               fields(nfields)%dom_start(1)   = (fields(nfields)%dom_start(1)               - 1)*sr_x + 1

               fields(nfields)%mem_end(1)     = (fields(nfields)%mem_end(1) - rhalo_width)*sr_x + rhalo_width
               fields(nfields)%patch_end(1)   = (fields(nfields)%patch_end(1)            )*sr_x
               fields(nfields)%dom_end(1)     = (fields(nfields)%dom_end(1)              )*sr_x

               fields(nfields)%mem_start(2)   = (fields(nfields)%mem_start(2) + bhalo_width - 1)*sr_y + 1 - bhalo_width
               fields(nfields)%patch_start(2) = (fields(nfields)%patch_start(2)             - 1)*sr_y + 1
               fields(nfields)%dom_start(2)   = (fields(nfields)%dom_start(2)               - 1)*sr_y + 1

               fields(nfields)%mem_end(2)     = (fields(nfields)%mem_end(2) - thalo_width)*sr_y + thalo_width
               fields(nfields)%patch_end(2)   = (fields(nfields)%patch_end(2)            )*sr_y
               fields(nfields)%dom_end(2)     = (fields(nfields)%dom_end(2)              )*sr_y

           endif

            nfields = nfields + 1
   
         end if  ! the next field given by get_next_fieldname() is valid }
    
      end do  ! for each user-specified field }
  
   end subroutine init_output_fields
 
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
   ! Name: write_field
   !
   ! Purpose: This routine writes the provided field to any output devices or APIs
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
   subroutine write_field(start_mem_i, end_mem_i, &
                          start_mem_j, end_mem_j, &
                          start_mem_k, end_mem_k, &
                          cname, datestr, real_array, is_training)
 
      implicit none
  
      ! Arguments
      integer, intent(in) :: start_mem_i, end_mem_i, start_mem_j, end_mem_j, start_mem_k, end_mem_k
      real, target, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j, start_mem_k:end_mem_k), &
                              intent(in) :: real_array
      logical, intent(in), optional :: is_training
      character (len=19), intent(in) :: datestr
      character (len=*), intent(in) :: cname
  
#include "wrf_io_flags.h"
#include "wrf_status_codes.h"
  
      ! Local variables
      integer :: i
      integer :: istatus, comm_1, comm_2, domain_desc
      integer, dimension(3) :: sd, ed, sp, ep, sm, em
      real, pointer, dimension(:,:,:) :: real_dom_array
      logical :: allocated_real_locally
      integer :: is_subgrid
  
      allocated_real_locally = .false.
  
      ! If we are running distributed memory and need to gather all tiles onto a single processor for output
      if (nprocs > 1 .and. .not. do_tiled_output) then
         do i=1,NUM_FIELDS
            if ( (index(cname, trim(fields(i)%fieldname) ) /= 0) .and. &
                  (len_trim(cname) == len_trim(fields(i)%fieldname)) ) then
               istatus = 0
     
               ! For the gather routines below, the IO_NODE should give the full domain dimensions, but the
               !   memory and patch dimensions should indicate what the processor already has in its patch_array.
               ! This is because an array with dimensions of the full domain will be allocated, and the patch_array
               !   will be copied from local memory into the full domain array in the area specified by the patch
               !   dimensions.
               sd = fields(i)%dom_start
               ed = fields(i)%dom_end
               sp = fields(i)%patch_start
               ep = fields(i)%patch_end
               sm = fields(i)%mem_start
               em = fields(i)%mem_end
     
               allocate(real_dom_array(sd(1):ed(1),sd(2):ed(2),sd(3):ed(3)))
               allocated_real_locally = .true.
               call gather_whole_field_r(real_array, &
                                         sm(1), em(1), sm(2), em(2), sm(3), em(3), &
                                         sp(1), ep(1), sp(2), ep(2), sp(3), ep(3), &
                                         real_dom_array, &
                                         sd(1), ed(1), sd(2), ed(2), sd(3), ed(3))
               exit
            end if 
         end do
      else
         do i=1,NUM_FIELDS
            if ( (index(cname, trim(fields(i)%fieldname) ) /= 0) .and. &
                 (len_trim(cname) == len_trim(fields(i)%fieldname)) ) then
               istatus = 0
               real_dom_array => real_array
               exit
            end if 
         end do
      end if
  
      ! Now a write call is only done if each processor writes its own file, or if we are the IO_NODE
      if (my_proc_id == IO_NODE .or. do_tiled_output) then
         comm_1 = 1
         comm_2 = 1
         domain_desc = 0
   
         do i=1,NUM_FIELDS
            if ( (index(cname, trim(fields(i)%fieldname) ) /= 0) .and. &
                 (len_trim(cname) == len_trim(fields(i)%fieldname)) ) then
   
               if (fields(i)%is_subgrid) then
                 is_subgrid=1
               else
                 is_subgrid=0
               end if

               ! Here, the output array has dimensions of the full grid if it was gathered together
               !   from all processors
               if (my_proc_id == IO_NODE .and. nprocs > 1 .and. .not. do_tiled_output) then
                  sd = fields(i)%dom_start
                  ed = fields(i)%dom_end
                  sm = sd
                  em = ed
                  sp = sd  
                  ep = ed
               ! If we are writing one file per processor, then each processor only writes out the 
               !   part of the domain that it has in memory
               else
! BUG: Shouldn't we set sd/ed to be domain_start/domain_end?
!      Maybe not, since patch is already adjusted for staggering; but maybe so, and also adjust
!      for staggering if it is alright to pass true domain dimensions to write_field.
                  sd = fields(i)%patch_start
                  ed = fields(i)%patch_end
                  sp = fields(i)%patch_start
                  ep = fields(i)%patch_end
                  sm = fields(i)%mem_start
                  em = fields(i)%mem_end
               end if
     
               istatus = 0
#ifdef IO_BINARY
               if (io_form_output == BINARY) then
                  call ext_int_write_field(handle, datestr, trim(fields(i)%fieldname), &
                       real_dom_array, WRF_REAL, comm_1, comm_2, domain_desc, trim(fields(i)%mem_order), &
                       trim(fields(i)%stagger), fields(i)%dimnames, sd, ed, sm, em, sp, ep, istatus)
               end if
#endif
#ifdef IO_NETCDF
               if (io_form_output == NETCDF) then
                  call ext_ncd_write_field(handle, datestr, trim(fields(i)%fieldname), &
                       real_dom_array, WRF_REAL, comm_1, comm_2, domain_desc, trim(fields(i)%mem_order), &
                       trim(fields(i)%stagger), fields(i)%dimnames, sd, ed, sm, em, sp, ep, istatus)
               end if
#endif
#ifdef IO_GRIB1
               if (io_form_output == GRIB1) then
                  call ext_gr1_write_field(handle, datestr, trim(fields(i)%fieldname), &
                       real_dom_array, WRF_REAL, comm_1, comm_2, domain_desc, trim(fields(i)%mem_order), &
                       trim(fields(i)%stagger), fields(i)%dimnames, sd, ed, sm, em, sp, ep, istatus)
               end if
#endif
               call mprintf((istatus /= 0),ERROR,'Error in ext_pkg_write_field')

               if (present(is_training)) then
                  if (is_training) then
#ifdef IO_BINARY
                     if (io_form_output == BINARY) then
                        call ext_int_put_var_ti_char(handle, 'units', &
                                                trim(fields(i)%fieldname), trim(fields(i)%units), istatus)
                        call ext_int_put_var_ti_char(handle, 'description', &
                                                trim(fields(i)%fieldname), trim(fields(i)%descr), istatus)
                        call ext_int_put_var_ti_char(handle, 'stagger', &
                                                trim(fields(i)%fieldname), trim(fields(i)%stagger), istatus)
                        call ext_int_put_var_ti_integer(handle,'subgrid', &
                                                trim(fields(i)%fieldname),(/is_subgrid/),1, istatus)
                     end if
#endif
#ifdef IO_NETCDF
                     if (io_form_output == NETCDF) then
                        call ext_ncd_put_var_ti_char(handle, 'units', &
                                                trim(fields(i)%fieldname), trim(fields(i)%units), istatus)
                        call ext_ncd_put_var_ti_char(handle, 'description', &
                                                trim(fields(i)%fieldname), trim(fields(i)%descr), istatus)
                        call ext_ncd_put_var_ti_char(handle, 'stagger', &
                                                trim(fields(i)%fieldname), trim(fields(i)%stagger), istatus)
                        call ext_ncd_put_var_ti_integer(handle,'subgrid', &
                                                trim(fields(i)%fieldname),(/is_subgrid/),1, istatus)
                    end if
#endif
#ifdef IO_GRIB1
                     if (io_form_output == GRIB1) then
                        call ext_gr1_put_var_ti_char(handle, 'units', &
                                                trim(fields(i)%fieldname), trim(fields(i)%units), istatus)
                        call ext_gr1_put_var_ti_char(handle, 'description', &
                                                trim(fields(i)%fieldname), trim(fields(i)%descr), istatus)
                        call ext_gr1_put_var_ti_char(handle, 'stagger', &
                                                trim(fields(i)%fieldname), trim(fields(i)%stagger), istatus)
                        call ext_gr1_put_var_ti_integer(handle,'subgrid', &
                                                trim(fields(i)%fieldname),(/is_subgrid/),1, istatus)
                    end if
#endif
                  end if
               end if
               exit
            end if
         end do
   
      end if
  
      if (allocated_real_locally) deallocate(real_dom_array)
  
   end subroutine write_field
 
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
   ! Name: write_global_attrs
   !
   ! Purpose:
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
   subroutine write_global_attrs(title, start_date, grid_type, dyn_opt,                             &
                                west_east_dim, south_north_dim, bottom_top_dim,                     &
                                we_patch_s, we_patch_e, we_patch_s_stag, we_patch_e_stag,           &
                                sn_patch_s, sn_patch_e, sn_patch_s_stag, sn_patch_e_stag,           &
                                map_proj, cmminlu, num_land_cat, is_water, is_lake, is_ice,         &
                                is_urban, i_soilwater, grid_id, parent_id,                          &
                                i_parent_start, j_parent_start, i_parent_end, j_parent_end,         &
                                dx, dy, cen_lat, moad_cen_lat, cen_lon,                             &
                                stand_lon, truelat1, truelat2, pole_lat, pole_lon,                  &
                                parent_grid_ratio, sr_x, sr_y, corner_lats, corner_lons,            &
                                num_metgrid_soil_levs,                                              &
                                flags, nflags, flag_excluded_middle)
 
      implicit none
  
      ! Arguments
      integer, intent(in) :: dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
                 we_patch_s, we_patch_e, we_patch_s_stag, we_patch_e_stag,            &
                 sn_patch_s, sn_patch_e, sn_patch_s_stag, sn_patch_e_stag,            &
                 map_proj, is_water, is_lake, is_ice, is_urban, i_soilwater,          &
                 grid_id, parent_id, i_parent_start, j_parent_start,                  &
                 i_parent_end, j_parent_end, parent_grid_ratio, sr_x, sr_y, num_land_cat
      integer, intent(in), optional :: num_metgrid_soil_levs
      integer, intent(in), optional :: nflags
      integer, intent(in), optional :: flag_excluded_middle
      real, intent(in) :: dx, dy, cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
                 pole_lat, pole_lon
      real, dimension(16), intent(in) :: corner_lats, corner_lons
      character (len=*), intent(in) :: title, start_date, grid_type
      character (len=128), intent(in) :: cmminlu
      character (len=128), dimension(:), intent(in), optional :: flags
  
      ! Local variables
      integer :: local_we_patch_s, local_we_patch_s_stag, &
                 local_we_patch_e, local_we_patch_e_stag, &
                 local_sn_patch_s, local_sn_patch_s_stag, &
                 local_sn_patch_e, local_sn_patch_e_stag
      integer :: i
      real, dimension(16) :: local_corner_lats, local_corner_lons

      local_we_patch_s      = we_patch_s
      local_we_patch_s_stag = we_patch_s_stag 
      local_we_patch_e      = we_patch_e
      local_we_patch_e_stag = we_patch_e_stag 
      local_sn_patch_s      = sn_patch_s
      local_sn_patch_s_stag = sn_patch_s_stag 
      local_sn_patch_e      = sn_patch_e
      local_sn_patch_e_stag = sn_patch_e_stag 
      local_corner_lats = corner_lats
      local_corner_lons = corner_lons

      if (nprocs > 1) then

         if (.not. do_tiled_output) then
            call parallel_bcast_int(local_we_patch_s,      processors(0, 0))
            call parallel_bcast_int(local_we_patch_s_stag, processors(0, 0))
            call parallel_bcast_int(local_sn_patch_s,      processors(0, 0))
            call parallel_bcast_int(local_sn_patch_s_stag, processors(0, 0))

            call parallel_bcast_int(local_we_patch_e,      processors(nproc_x-1, nproc_y-1))
            call parallel_bcast_int(local_we_patch_e_stag, processors(nproc_x-1, nproc_y-1))
            call parallel_bcast_int(local_sn_patch_e,      processors(nproc_x-1, nproc_y-1))
            call parallel_bcast_int(local_sn_patch_e_stag, processors(nproc_x-1, nproc_y-1))
         end if

         call parallel_bcast_real(local_corner_lats(1),  processors(0,         0))
         call parallel_bcast_real(local_corner_lats(2),  processors(0,         nproc_y-1))
         call parallel_bcast_real(local_corner_lats(3),  processors(nproc_x-…