From eba6fe1b7c56f8303f1a59186300acfa5c72e73f Mon Sep 17 00:00:00 2001
From: Jesse Nusbaumer <nusbaume@ucar.edu>
Date: Thu, 7 Jan 2021 15:55:57 -0700
Subject: [PATCH 1/8] Replace 'get_dyn_grid_info' pointer method with
 simplified array-passing interface for Null dycore initialization.

---
 src/control/cam_comp.F90                  |  10 +-
 src/dynamics/none/dyn_grid.F90            | 169 +++++++++++-----------
 src/physics/utils/physics_column_type.F90 |  47 +++++-
 src/physics/utils/physics_grid.F90        |  81 ++++++-----
 4 files changed, 180 insertions(+), 127 deletions(-)

diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90
index 636f6b42..9123a2fd 100644
--- a/src/control/cam_comp.F90
+++ b/src/control/cam_comp.F90
@@ -67,8 +67,7 @@ subroutine cam_init(caseid, ctitle, model_doi_url,                         &
 
 !      use history_defaults,     only: bldfld
       use cam_initfiles,        only: cam_initfiles_open
-      use dyn_grid,             only: dyn_grid_init
-      use physics_grid,         only: phys_grid_init
+      use dyn_grid,             only: model_grid_init
       use phys_comp,            only: phys_init
       use dyn_comp,             only: dyn_init
 !      use cam_restart,          only: cam_read_restart
@@ -149,11 +148,8 @@ subroutine cam_init(caseid, ctitle, model_doi_url,                         &
       ! Open initial or restart file, and topo file if specified.
       call cam_initfiles_open()
 
-      ! Initialize grids and dynamics grid decomposition
-      call dyn_grid_init()
-
-      ! Initialize physics grid decomposition
-      call phys_grid_init()
+      ! Initialize model grids and decompositions
+      call model_grid_init()
 
       ! Initialize ghg surface values before default initial distributions
       ! are set in dyn_init
diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90
index 4467b714..fa9c6d5d 100644
--- a/src/dynamics/none/dyn_grid.F90
+++ b/src/dynamics/none/dyn_grid.F90
@@ -10,12 +10,12 @@ module dyn_grid
    private
    save
 
-   public dyn_grid_init
-   public get_dyn_grid_info
-   public physgrid_copy_attributes_d
+   public model_grid_init
 
-   type(physics_column_t), public, protected, allocatable :: local_columns(:)
    ! Private module variables
+
+   type(physics_column_t), allocatable :: dyn_columns(:)
+
    integer               :: num_levels         = -1
    integer               :: num_global_columns = -1
    integer               :: num_lats           = -1 ! Global
@@ -52,50 +52,57 @@ module dyn_grid
 CONTAINS
 !==============================================================================
 
-   subroutine dyn_grid_init()
-      use pio,            only: file_desc_t, var_desc_t, io_desc_t
-      use pio,            only: iMap=>PIO_OFFSET_KIND, PIO_DOUBLE
-      use pio,            only: PIO_BCAST_ERROR, pio_seterrorhandling
-      use pio,            only: pio_get_var, pio_freedecomp
-      use pio,            only: pio_read_darray
-      use spmd_utils,     only: npes, iam, masterprocid, mpicom
-      use cam_pio_utils,  only: cam_pio_handle_error, cam_pio_find_var
-      use cam_pio_utils,  only: cam_pio_var_info, pio_subsystem
-      use cam_pio_utils,  only: cam_pio_newdecomp
-      use cam_abortutils, only: endrun
-      use cam_logfile,    only: cam_log_multiwrite
-      use cam_initfiles,  only: initial_file_get_id
+   subroutine model_grid_init()
+      use shr_kind_mod,     only: SHR_KIND_CL
+      use pio,              only: file_desc_t, var_desc_t, io_desc_t
+      use pio,              only: iMap=>PIO_OFFSET_KIND, PIO_DOUBLE
+      use pio,              only: PIO_BCAST_ERROR, pio_seterrorhandling
+      use pio,              only: pio_get_var, pio_freedecomp
+      use pio,              only: pio_read_darray
+      use spmd_utils,       only: npes, iam, masterprocid, mpicom
+      use cam_pio_utils,    only: cam_pio_handle_error, cam_pio_find_var
+      use cam_pio_utils,    only: cam_pio_var_info, pio_subsystem
+      use cam_pio_utils,    only: cam_pio_newdecomp
+      use cam_abortutils,   only: endrun
+      use cam_logfile,      only: cam_log_multiwrite
+      use cam_initfiles,    only: initial_file_get_id
+      use physics_grid,     only: phys_grid_init
+      use cam_grid_support, only: hclen => max_hcoordname_len
 
-      ! Initialize a dynamics decomposition based on an input data file
+      ! Initializes a dynamics decomposition based on an input data file,
+      ! and then initializes the physics decomposition based on the dynamics
+      ! grid.
 
       ! Local variables
-      type(file_desc_t), pointer     :: fh_ini
-      type(var_desc_t)               :: lat_vardesc
-      type(var_desc_t)               :: lon_vardesc
-      type(var_desc_t)               :: vardesc
-      type(io_desc_t),   pointer     :: iodesc
-      integer                        :: err_handling
-      logical                        :: var_found
-      logical                        :: is_degrees
-      logical                        :: is_lat
-      integer                        :: num_var_dims
-      integer                        :: time_id
-      integer                        :: lindex
-      integer                        :: dimids(MAX_DIMS)
-      integer                        :: dimlens(MAX_DIMS)
-      integer                        :: col_mod ! Temp for calculating decomp
-      integer                        :: col_start, col_end
-      integer                        :: start(1), kount(1)
-      integer                        :: iret
-      integer(iMap),     allocatable :: ldof(:) ! For reading coordinates
-      real(r8),          allocatable :: temp_arr(:)
-      character(len=128)             :: var_name
-      character(len=256)             :: dimnames(MAX_DIMS)
-      character(len=8)               :: lat_dim_name
-      character(len=8)               :: lon_dim_name
-      character(len=128)             :: errormsg
-
-      character(len=*),  parameter   :: subname = 'dyn_grid_init'
+      type(file_desc_t), pointer        :: fh_ini
+      type(var_desc_t)                  :: lat_vardesc
+      type(var_desc_t)                  :: lon_vardesc
+      type(var_desc_t)                  :: vardesc
+      type(io_desc_t),   pointer        :: iodesc
+      integer                           :: err_handling
+      logical                           :: var_found
+      logical                           :: is_degrees
+      logical                           :: is_lat
+      integer                           :: num_var_dims
+      integer                           :: time_id
+      integer                           :: lindex
+      integer                           :: dimids(MAX_DIMS)
+      integer                           :: dimlens(MAX_DIMS)
+      integer                           :: col_mod ! Temp for calculating decomp
+      integer                           :: col_start, col_end
+      integer                           :: start(1), kount(1)
+      integer                           :: iret
+      integer(iMap),     allocatable    :: ldof(:) ! For reading coordinates
+      real(r8),          allocatable    :: temp_arr(:)
+      character(len=128)                :: var_name
+      character(len=hclen), allocatable :: grid_attribute_names(:)
+      character(len=SHR_KIND_CL)        :: dimnames(MAX_DIMS)
+      character(len=8)                  :: lat_dim_name
+      character(len=8)                  :: lon_dim_name
+      character(len=128)                :: errormsg
+
+      character(len=*),  parameter :: subname = 'dyn_grid_init'
+
 
       nullify(iodesc)
 
@@ -341,42 +348,53 @@ subroutine dyn_grid_init()
       ! Back to old error handling
       call pio_seterrorhandling(fh_ini, err_handling)
 
-   end subroutine dyn_grid_init
+      ! Set dyn_columns values:
+      call set_dyn_col_values()
+
+      ! Set dynamics grid attributes
+      allocate(grid_attribute_names(0))
+
+      ! Initialize physics grid decomposition:
+      call phys_grid_init(num_lons, num_lats, num_levels, 'NULL', &
+                          1, num_levels, dyn_columns, gridname, &
+                          grid_attribute_names)
+
+      ! Deallocate grid_attirbute_names, as it is no longer needed:
+      deallocate(grid_attribute_names)
+
+      ! Deallocate dyn_columns, as it is now stored in the
+      ! global phys_columns structure:
+      deallocate(dyn_columns)
+
+   end subroutine model_grid_init
 
    !===========================================================================
 
-   subroutine get_dyn_grid_info(hdim1_d, hdim2_d, num_lev,                    &
-        dycore_name, index_model_top_layer, index_surface_layer, dyn_columns)
+   subroutine set_dyn_col_values()
+
+      ! Sets the values stored in the "dyn_columns" structure,
+      ! which are the physics columns as they exist on the
+      ! dynamics decomposition.
+
       use shr_const_mod,  only: SHR_CONST_PI
       use cam_abortutils, only: endrun
       use spmd_utils,     only: iam
+
       ! Dummy arguments
-      integer,          intent(out)   :: hdim1_d ! # longitudes or grid size
-      integer,          intent(out)   :: hdim2_d ! # latitudes or 1
-      integer,          intent(out)   :: num_lev ! # levels
-      character(len=*), intent(out)   :: dycore_name
-      integer,          intent(out)   :: index_model_top_layer
-      integer,          intent(out)   :: index_surface_layer
-      type(physics_column_t), pointer :: dyn_columns(:) ! Phys col in Dyn decomp
-      ! Local variables
       integer                         :: lindex
       integer                         :: gindex
       integer                         :: lat_index, lat1
       integer                         :: lon_index
       real(r8),         parameter     :: radtodeg = 180.0_r8 / SHR_CONST_PI
       real(r8),         parameter     :: degtorad = SHR_CONST_PI / 180.0_r8
-      character(len=*), parameter     :: subname = 'get_dyn_grid_info'
+      character(len=*), parameter     :: subname = 'set_dyn_col_values'
 
-      if (associated(dyn_columns)) then
-         call endrun(subname//': dyn_columns must be unassociated pointer')
+      ! Allocate dyn_columns structure if not already allocated:
+      if (.not.allocated(dyn_columns)) then
+         allocate(dyn_columns(num_local_columns))
       end if
-      allocate(dyn_columns(num_local_columns))
-      hdim1_d = num_lons
-      hdim2_d = num_lats
-      num_lev = num_levels
-      dycore_name = 'NULL'
-      index_model_top_layer = 1
-      index_surface_layer = num_levels
+
+      ! Calculate dyn_columns variable values:
       lat1 = global_col_offset / num_lons
       do lindex = 1, num_local_columns
          if (grid_is_latlon) then
@@ -433,24 +451,7 @@ subroutine get_dyn_grid_info(hdim1_d, hdim2_d, num_lev,                    &
          dyn_columns(lindex)%dyn_block_index(1) = lindex
       end do
 
-   end subroutine get_dyn_grid_info
-
-   !===========================================================================
-
-   subroutine physgrid_copy_attributes_d(gridname_out, grid_attribute_names)
-      ! create list of attributes for the physics grid that should be copied
-      ! from the corresponding grid object on the dynamics decomposition
-
-      use cam_grid_support, only: hclen => max_hcoordname_len
-
-      ! Dummy arguments
-      character(len=hclen),          intent(out) :: gridname_out
-      character(len=hclen), pointer, intent(out) :: grid_attribute_names(:)
-
-      gridname_out = gridname
-      allocate(grid_attribute_names(0))
-
-   end subroutine physgrid_copy_attributes_d
+   end subroutine set_dyn_col_values
 
    !===========================================================================
 
diff --git a/src/physics/utils/physics_column_type.F90 b/src/physics/utils/physics_column_type.F90
index ba04655c..e665f6ce 100644
--- a/src/physics/utils/physics_column_type.F90
+++ b/src/physics/utils/physics_column_type.F90
@@ -1,6 +1,6 @@
 module physics_column_type
 
-   use shr_kind_mod, only: r8 => shr_kind_r8
+   use shr_kind_mod, only: r8 => shr_kind_r8, shr_kind_cl
 !   use ISO_FORTRAN_ENV, only: kind_phys
    use ccpp_kinds,   only: kind_phys
  
@@ -9,6 +9,12 @@ module physics_column_type
    private
    save
 
+   public :: assignment ( = )
+
+   interface assignment ( = )
+      module procedure copy_phys_col
+   end interface
+
 !> \section arg_table_physics_column_t  Argument Table
 !! \htmlinclude physics_column_t.html
 type, public :: physics_column_t
@@ -36,9 +42,46 @@ module physics_column_type
    integer              :: phys_chunk_index = -1   ! Index into physics chunk
 end type physics_column_t
 
-
 !==============================================================================
 CONTAINS
 !==============================================================================
 
+
+   subroutine copy_phys_col(phys_col_out, phys_col_in)
+
+      use cam_abortutils, only: endrun
+
+      ! Copy all of the values from one physics_column type
+      ! structure to another.
+
+      ! Dummy (input) variables:
+      type(physics_column_t), intent(out) :: phys_col_out
+      type(physics_column_t), intent(in)  :: phys_col_in
+
+      ! Copy values from input array to output array:
+
+      ! Column information
+      phys_col_out%lat_rad = phys_col_in%lat_rad
+      phys_col_out%lon_rad = phys_col_in%lon_rad
+      phys_col_out%lat_deg = phys_col_in%lat_deg
+      phys_col_out%lon_deg = phys_col_in%lon_deg
+      phys_col_out%area    = phys_col_in%area
+      phys_col_out%weight  = phys_col_in%weight
+
+      ! File decomposition
+      phys_col_out%global_col_num   = phys_col_in%global_col_num
+      phys_col_out%coord_indices(:) = phys_col_in%coord_indices(:)
+
+      ! Dynamics decomposition
+      phys_col_out%dyn_task         = phys_col_in%dyn_task
+      phys_col_out%local_dyn_block  = phys_col_in%local_dyn_block
+      phys_col_out%global_dyn_block = phys_col_in%global_dyn_block
+
+      ! Physics decomposition
+      phys_col_out%phys_task        = phys_col_in%phys_task
+      phys_col_out%local_phys_chunk = phys_col_in%local_phys_chunk
+      phys_col_out%phys_chunk_index = phys_col_in%phys_chunk_index
+
+   end subroutine copy_phys_col
+
 end module physics_column_type
diff --git a/src/physics/utils/physics_grid.F90 b/src/physics/utils/physics_grid.F90
index 5ccc694e..edb3c845 100644
--- a/src/physics/utils/physics_grid.F90
+++ b/src/physics/utils/physics_grid.F90
@@ -1,7 +1,7 @@
 module physics_grid
 
    use shr_kind_mod,        only: r8 => shr_kind_r8
-   use physics_column_type, only: physics_column_t
+   use physics_column_type, only: physics_column_t, assignment(=)
    use perf_mod,            only: t_adj_detailf, t_startf, t_stopf
 
    implicit none
@@ -35,7 +35,7 @@ module physics_grid
    character(len=8), protected, public :: dycore_name = ''
 
    ! Physics decomposition information
-   type(physics_column_t), pointer     :: phys_columns(:) => NULL()
+   type(physics_column_t), protected, public, allocatable :: phys_columns(:)
 
    ! These variables are last to provide a limited table to search
 
@@ -56,22 +56,36 @@ module physics_grid
 CONTAINS
 !==============================================================================
 
-   subroutine phys_grid_init()
+   subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
+                             index_top_layer_in, index_bottom_layer_in, &
+                             dyn_columns, dyn_gridname, dyn_attributes)
+
 !      use mpi,              only: MPI_reduce ! XXgoldyXX: Should this work?
       use mpi,              only: MPI_INTEGER, MPI_MIN
       use cam_abortutils,   only: endrun
-      use spmd_utils,       only: npes, mpicom
-      use dyn_grid,         only: get_dyn_grid_info, physgrid_copy_attributes_d
+      use spmd_utils,       only: npes, mpicom, iam
       use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register
       use cam_grid_support, only: iMap, hclen => max_hcoordname_len
       use cam_grid_support, only: horiz_coord_t, horiz_coord_create
       use cam_grid_support, only: cam_grid_attribute_copy, cam_grid_attr_exists
 
+      ! Dummy (input) variables:
+      integer, intent(in) :: hdim1_d_in            ! First dyn grid horizontal dimension
+      integer, intent(in) :: hdim2_d_in            ! Second dyn grid horizontal dimension
+      integer, intent(in) :: pver_in               ! Dyn grid vertical dimension
+      integer, intent(in) :: index_top_layer_in    ! Vertical index that represents model top
+      integer, intent(in) :: index_bottom_layer_in ! Vertical index that represents model surface
+
+      character(len=*), intent(in) :: dycore_name_in    ! Name of dycore
+      character(len=*), intent(in) :: dyn_gridname      ! Name of dynamics grid
+      character(len=*), intent(in) :: dyn_attributes(:) ! Dyanmics grid attributes
+
+      type(physics_column_t), intent(in) :: dyn_columns(:) ! physics columns structure on dynamics grid
+
       ! Local variables
       integer                             :: index
       integer                             :: col_index
       integer                             :: first_dyn_column, last_dyn_column
-      type(physics_column_t), pointer     :: dyn_columns(:) ! Dyn decomp
       ! Maps and values for physics grid
       real(r8),               pointer     :: lonvals(:)
       real(r8),               pointer     :: latvals(:)
@@ -84,29 +98,29 @@ subroutine phys_grid_init()
       logical                             :: unstructured
       real(r8)                            :: temp ! For MPI
       integer                             :: ierr ! For MPI
-      character(len=hclen),   pointer     :: copy_attributes(:)
-      character(len=hclen)                :: copy_gridname
 
-      nullify(dyn_columns)
       nullify(lonvals)
       nullify(latvals)
       nullify(grid_map)
       nullify(lat_coord)
       nullify(lon_coord)
       nullify(area_d)
-      nullify(copy_attributes)
 
       call t_adj_detailf(-2)
       call t_startf("phys_grid_init")
 
-      ! Gather info from the dycore
-      call get_dyn_grid_info(hdim1_d, hdim2_d, pver, dycore_name,             &
-           index_top_layer, index_bottom_layer, dyn_columns)
-      num_global_phys_cols = hdim1_d * hdim2_d
-      pverp = pver + 1
+      ! Set public variables:
+      hdim1_d            = hdim1_d_in
+      hdim2_d            = hdim2_d_in
+      pver               = pver_in
+      index_top_layer    = index_top_layer_in
+      index_bottom_layer = index_bottom_layer_in
+      dycore_name        = dycore_name_in
+
+      pverp            = pver + 1
       first_dyn_column = LBOUND(dyn_columns, 1)
-      last_dyn_column = UBOUND(dyn_columns, 1)
-      unstructured = hdim2_d <= 1
+      last_dyn_column  = UBOUND(dyn_columns, 1)
+      unstructured     = hdim2_d <= 1
       !!XXgoldyXX: Can we enforce interface numbering separate from dycore?
       !!XXgoldyXX: This will work for both CAM and WRF/MPAS physics
       !!XXgoldyXX: This only has a 50% chance of working on a single level model
@@ -118,15 +132,20 @@ subroutine phys_grid_init()
          index_top_interface = index_top_layer + 1
       end if
 
-      ! Set up the physics decomposition
-      columns_on_task = size(dyn_columns)
-      phys_columns => dyn_columns
+      !Calculate total number of physics columns:
+      num_global_phys_cols = hdim1_d * hdim2_d
 
-      ! Now that we are done settine up the physics decomposition, clean up
-      if (.not. associated(phys_columns, target=dyn_columns)) then
-         deallocate(dyn_columns)
+      ! Set columns_on_task and allocate phys_columns if
+      ! not already allocated:
+      if (.not. allocated(phys_columns)) then
+          columns_on_task = size(dyn_columns)
+          allocate(phys_columns(columns_on_task))
       end if
-      nullify(dyn_columns)
+
+      ! Set up the physics decomposition
+      do index = 1, columns_on_task
+         phys_columns(index) = dyn_columns(index)
+      end do
 
       ! Add physics-package grid to set of CAM grids
       ! physgrid always uses 'lat' and 'lon' as coordinate names; If dynamics
@@ -214,12 +233,11 @@ subroutine phys_grid_init()
       call cam_grid_register('physgrid', phys_decomp,                         &
            lat_coord, lon_coord, grid_map, src_in=(/ 1, 0 /),                 &
            unstruct=unstructured, block_indexed=.false.)
+
       ! Copy required attributes from the dynamics array
-      nullify(copy_attributes)
-      call physgrid_copy_attributes_d(copy_gridname, copy_attributes)
-      do index = 1, size(copy_attributes)
-         call cam_grid_attribute_copy(copy_gridname, 'physgrid',              &
-              copy_attributes(index))
+      do index = 1, size(dyn_attributes)
+         call cam_grid_attribute_copy(dyn_gridname, 'physgrid',              &
+              dyn_attributes(index))
       end do
 
       if ((.not. cam_grid_attr_exists('physgrid', 'area')) .and.              &
@@ -242,11 +260,6 @@ subroutine phys_grid_init()
       nullify(latvals)
       deallocate(lonvals)
       nullify(lonvals)
-      ! Cleanup, we are responsible for copy attributes
-      if (associated(copy_attributes)) then
-         deallocate(copy_attributes)
-         nullify(copy_attributes)
-      end if
 
       ! Set flag indicating physics grid is now set
       phys_grid_initialized = .true.

From ec31ef22e313e8143cab4fe41964b371bb374910 Mon Sep 17 00:00:00 2001
From: Jesse Nusbaumer <nusbaume@ucar.edu>
Date: Wed, 13 Jan 2021 20:45:26 -0700
Subject: [PATCH 2/8] Add missing variables found during code review, and
 improve error-checking.

---
 src/dynamics/none/dyn_grid.F90            | 92 +++++++++++++++++++----
 src/physics/utils/physics_column_type.F90 | 28 ++++++-
 src/physics/utils/physics_grid.F90        | 63 +++++++++++++---
 3 files changed, 155 insertions(+), 28 deletions(-)

diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90
index fa9c6d5d..4e621f63 100644
--- a/src/dynamics/none/dyn_grid.F90
+++ b/src/dynamics/none/dyn_grid.F90
@@ -101,7 +101,7 @@ subroutine model_grid_init()
       character(len=8)                  :: lon_dim_name
       character(len=128)                :: errormsg
 
-      character(len=*),  parameter :: subname = 'dyn_grid_init'
+      character(len=*),  parameter :: subname = 'model_grid_init'
 
 
       nullify(iodesc)
@@ -264,9 +264,24 @@ subroutine model_grid_init()
                  '(a,i4,i9,2i7)', (/ kount(1), start(1), start(1)+kount(1)-1/))
          end if
          if (is_degrees) then
-            allocate(local_lons_deg(num_lons))
-            allocate(local_lats_deg(kount(1)))
-            allocate(temp_arr(num_lats))
+            allocate(local_lons_deg(num_lons), stat=iret)
+            if (iret /= 0) then
+               write(errormsg, *) &
+                  subname//': allocate local_lons_deg(num_lons) failed with stat: ',iret
+               call endrun(errormsg)
+            end if
+            allocate(local_lats_deg(kount(1)), stat=iret)
+            if (iret /= 0) then
+               write(errormsg, *) &
+                  subname//': allocate local_lats_deg(kount) failed with stat: ',iret
+               call endrun(errormsg)
+            end if
+            allocate(temp_arr(num_lats), stat=iret)
+            if (iret /= 0) then
+               write(errormsg, *) &
+                  subname//': allocate temp_arr failed with stat: ',iret
+               call endrun(errormsg)
+            end if
             iret = pio_get_var(fh_ini, lat_vardesc, (/ 1 /), (/ num_lats /),  &
                  temp_arr)
             call cam_pio_handle_error(iret,                                   &
@@ -287,14 +302,34 @@ subroutine model_grid_init()
       else
          ! Do parallel read of lat and lon
          if (is_degrees) then
-            allocate(local_lats_deg(num_local_columns))
-            allocate(local_lons_deg(num_local_columns))
-            allocate(ldof(num_local_columns))
+            allocate(local_lats_deg(num_local_columns), stat=iret)
+            if (iret /= 0) then
+               write(errormsg, *) &
+                  subname//': allocate local_lats_deg(columns) failed with stat: ',iret
+               call endrun(errormsg)
+            end if
+            allocate(local_lons_deg(num_local_columns), stat=iret)
+            if (iret /= 0) then
+               write(errormsg, *) &
+                  subname//': allocate local_lats_deg(columns) failed with stat: ',iret
+               call endrun(errormsg)
+            end if
+            allocate(ldof(num_local_columns), stat=iret)
+            if (iret /= 0) then
+               write(errormsg, *) &
+                  subname//': allocate ldof failed with stat: ',iret
+               call endrun(errormsg)
+            end if
             ldof = 0_iMap
             do lindex = 1, num_local_columns
                ldof(lindex) = col_start + lindex - 1
             end do
-            allocate(iodesc)
+            allocate(iodesc, stat=iret)
+            if (iret /= 0) then
+               write(errormsg, *) &
+                  subname//': allocate iodesc failed with stat: ',iret
+               call endrun(errormsg)
+            end if
             call cam_pio_newdecomp(iodesc, (/ num_global_columns /), ldof,    &
                  PIO_DOUBLE)
             call pio_read_darray(fh_ini, lat_vardesc, iodesc, local_lats_deg, &
@@ -321,14 +356,24 @@ subroutine model_grid_init()
             call endrun(errormsg)
          end if
          if ((num_lats > 1) .and. (dimlens(1) == num_lats)) then
-            allocate(local_areas(num_lats))
+            allocate(local_areas(num_lats), stat=iret)
+            if (iret /= 0) then
+               write(errormsg, *) &
+                  subname//': allocate local_areas(num_lats) failed with stat: ',iret
+               call endrun(errormsg)
+            end if
             start(1) = 1
             kount(1) = num_lats
             iret = pio_get_var(fh_ini, vardesc, start, kount, local_areas)
             call cam_pio_handle_error(iret,                                   &
                  subname//': Unable to read '//trim(var_name))
          else if (dimlens(1) == num_global_columns) then
-            allocate(local_areas(num_local_columns))
+            allocate(local_areas(num_local_columns), stat=iret)
+            if (iret /= 0) then
+               write(errormsg, *) &
+                  subname//': allocate local_areas(columns) failed with stat: ',iret
+               call endrun(errormsg)
+            end if
             call pio_read_darray(fh_ini, vardesc, iodesc, local_areas, iret)
             call cam_pio_handle_error(iret, subname//': Unable to read areas')
          else
@@ -352,7 +397,12 @@ subroutine model_grid_init()
       call set_dyn_col_values()
 
       ! Set dynamics grid attributes
-      allocate(grid_attribute_names(0))
+      allocate(grid_attribute_names(0), stat=iret)
+      if (iret /= 0) then
+         write(errormsg, *) &
+            subname//': allocate grid_attribute_names failed with stat: ',iret
+         call endrun(errormsg)
+      end if
 
       ! Initialize physics grid decomposition:
       call phys_grid_init(num_lons, num_lats, num_levels, 'NULL', &
@@ -380,18 +430,26 @@ subroutine set_dyn_col_values()
       use cam_abortutils, only: endrun
       use spmd_utils,     only: iam
 
-      ! Dummy arguments
+      ! Local variables:
       integer                         :: lindex
       integer                         :: gindex
       integer                         :: lat_index, lat1
       integer                         :: lon_index
+      integer                         :: ierr
+      character(len=128)              :: emsg
+
       real(r8),         parameter     :: radtodeg = 180.0_r8 / SHR_CONST_PI
       real(r8),         parameter     :: degtorad = SHR_CONST_PI / 180.0_r8
       character(len=*), parameter     :: subname = 'set_dyn_col_values'
 
       ! Allocate dyn_columns structure if not already allocated:
       if (.not.allocated(dyn_columns)) then
-         allocate(dyn_columns(num_local_columns))
+         allocate(dyn_columns(num_local_columns), stat=ierr)
+         if (ierr /= 0) then
+            write(emsg, *) &
+               subname//': allocate dyn_columns failed with stat: ',ierr
+            call endrun(emsg)
+         end if
       end if
 
       ! Calculate dyn_columns variable values:
@@ -447,7 +505,13 @@ subroutine set_dyn_col_values()
          dyn_columns(lindex)%global_dyn_block = iam + 1
          !  If there is more than one block lindex, they are in the same order
          !    as in the dynamics block structure
-         allocate(dyn_columns(lindex)%dyn_block_index(1))
+         allocate(dyn_columns(lindex)%dyn_block_index(1), stat=ierr)
+         if (ierr /= 0) then
+            write(emsg, *) &
+               subname//': allocate dyn_block_index failed with stat: ',ierr
+            call endrun(emsg)
+         end if
+
          dyn_columns(lindex)%dyn_block_index(1) = lindex
       end do
 
diff --git a/src/physics/utils/physics_column_type.F90 b/src/physics/utils/physics_column_type.F90
index e665f6ce..c4f5af93 100644
--- a/src/physics/utils/physics_column_type.F90
+++ b/src/physics/utils/physics_column_type.F90
@@ -1,9 +1,8 @@
 module physics_column_type
 
    use shr_kind_mod, only: r8 => shr_kind_r8, shr_kind_cl
-!   use ISO_FORTRAN_ENV, only: kind_phys
    use ccpp_kinds,   only: kind_phys
- 
+
 
    implicit none
    private
@@ -58,6 +57,12 @@ subroutine copy_phys_col(phys_col_out, phys_col_in)
       type(physics_column_t), intent(out) :: phys_col_out
       type(physics_column_t), intent(in)  :: phys_col_in
 
+      ! Local variables:
+      integer :: istat
+
+      character(len=*),  parameter :: subname = 'copy_phys_col'
+      character(len=128)           :: emsg
+
       ! Copy values from input array to output array:
 
       ! Column information
@@ -77,6 +82,25 @@ subroutine copy_phys_col(phys_col_out, phys_col_in)
       phys_col_out%local_dyn_block  = phys_col_in%local_dyn_block
       phys_col_out%global_dyn_block = phys_col_in%global_dyn_block
 
+      ! Dynamics blocks
+      if (allocated(phys_col_in%dyn_block_index)) then
+         ! De-allocate output block indices if already allocated:
+         if (allocated(phys_col_out%dyn_block_index)) then
+            deallocate(phys_col_out%dyn_block_index)
+         end if
+
+         ! Allocate output to match size of input:
+         allocate(phys_col_out%dyn_block_index(size(phys_col_in%dyn_block_index)), &
+                  stat=istat)
+         if (istat /= 0) then
+           write(emsg, *) &
+              subname//': allocate dyn_block_index failed with stat: ',istat
+           call endrun(emsg)
+         end if
+
+         phys_col_out%dyn_block_index(:) = phys_col_in%dyn_block_index(:)
+      end if
+
       ! Physics decomposition
       phys_col_out%phys_task        = phys_col_in%phys_task
       phys_col_out%local_phys_chunk = phys_col_in%local_phys_chunk
diff --git a/src/physics/utils/physics_grid.F90 b/src/physics/utils/physics_grid.F90
index edb3c845..264a662a 100644
--- a/src/physics/utils/physics_grid.F90
+++ b/src/physics/utils/physics_grid.F90
@@ -63,7 +63,7 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
 !      use mpi,              only: MPI_reduce ! XXgoldyXX: Should this work?
       use mpi,              only: MPI_INTEGER, MPI_MIN
       use cam_abortutils,   only: endrun
-      use spmd_utils,       only: npes, mpicom, iam
+      use spmd_utils,       only: npes, mpicom
       use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register
       use cam_grid_support, only: iMap, hclen => max_hcoordname_len
       use cam_grid_support, only: horiz_coord_t, horiz_coord_create
@@ -98,6 +98,9 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
       logical                             :: unstructured
       real(r8)                            :: temp ! For MPI
       integer                             :: ierr ! For MPI
+      character(len=128)                  :: emsg
+
+      character(len=*),  parameter :: subname = 'phys_grid_init'
 
       nullify(lonvals)
       nullify(latvals)
@@ -118,8 +121,6 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
       dycore_name        = dycore_name_in
 
       pverp            = pver + 1
-      first_dyn_column = LBOUND(dyn_columns, 1)
-      last_dyn_column  = UBOUND(dyn_columns, 1)
       unstructured     = hdim2_d <= 1
       !!XXgoldyXX: Can we enforce interface numbering separate from dycore?
       !!XXgoldyXX: This will work for both CAM and WRF/MPAS physics
@@ -138,12 +139,21 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
       ! Set columns_on_task and allocate phys_columns if
       ! not already allocated:
       if (.not. allocated(phys_columns)) then
-          columns_on_task = size(dyn_columns)
-          allocate(phys_columns(columns_on_task))
+         columns_on_task = size(dyn_columns)
+         allocate(phys_columns(columns_on_task), stat=ierr)
+         if (ierr /= 0) then
+            write(emsg, *) &
+               subname//': allocate phys_columns failed with stat: ',ierr
+            call endrun(emsg)
+         end if
       end if
 
+      ! Set column index bounds:
+      first_dyn_column = 1
+      last_dyn_column  = columns_on_task
+
       ! Set up the physics decomposition
-      do index = 1, columns_on_task
+      do index = first_dyn_column, last_dyn_column
          phys_columns(index) = dyn_columns(index)
       end do
 
@@ -155,13 +165,30 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
       ! It's structure will depend on whether or not the physics grid is
       ! unstructured
       if (unstructured) then
-         allocate(grid_map(3, columns_on_task))
+         allocate(grid_map(3, columns_on_task), stat=ierr)
       else
-         allocate(grid_map(4, columns_on_task))
+         allocate(grid_map(4, columns_on_task), stat=ierr)
+      end if
+      if (ierr /= 0) then
+         write(emsg, *) &
+            subname//': allocate grid_map failed with stat: ',ierr
+         call endrun(emsg)
       end if
       grid_map = 0
-      allocate(latvals(size(grid_map, 2)))
-      allocate(lonvals(size(grid_map, 2)))
+
+      allocate(latvals(size(grid_map, 2)), stat=ierr)
+      if (ierr /= 0) then
+         write(emsg, *) &
+            subname//': allocate latvals failed with stat: ',ierr
+         call endrun(emsg)
+      end if
+
+      allocate(lonvals(size(grid_map, 2)), stat=ierr)
+      if (ierr /= 0) then
+         write(emsg, *) &
+            subname//': allocate lonvals failed with stat: ',ierr
+         call endrun(emsg)
+      end if
 
       lonmin = 1000.0_r8 ! Out of longitude range
       latmin = 1000.0_r8 ! Out of latitude range
@@ -199,7 +226,13 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
               'latitude', 'degrees_north', 1, size(latvals), latvals,         &
               map=grid_map(3,:))
       else
-         allocate(coord_map(size(grid_map, 2)))
+         allocate(coord_map(size(grid_map, 2)), stat=ierr)
+         if (ierr /= 0) then
+            write(emsg, *) &
+               subname//': allocate coord_map failed with stat: ',ierr
+            call endrun(emsg)
+         end if
+
          ! We need a global minimum longitude and latitude
          if (npes > 1) then
             temp = lonmin
@@ -246,7 +279,13 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
          !   from the dycore (i.e., physics and dynamics are on different
          !   grids), create that attribute here (Note, a separate physics
          !   grid is only supported for unstructured grids).
-         allocate(area_d(size(grid_map, 2)))
+         allocate(area_d(size(grid_map, 2)), stat=ierr)
+         if (ierr /= 0) then
+            write(emsg, *) &
+               subname//': allocate area_d failed with stat: ',ierr
+            call endrun(emsg)
+         end if
+
          do col_index = 1, columns_on_task
             area_d(col_index) = phys_columns(col_index)%area
          end do

From d3d599fd00f92742847d708e60b28f1d4d4bfdd8 Mon Sep 17 00:00:00 2001
From: Jesse Nusbaumer <nusbaume@ucar.edu>
Date: Thu, 14 Jan 2021 11:32:44 -0700
Subject: [PATCH 3/8] Fix NAG warnings.

---
 src/dynamics/none/dyn_grid.F90            | 2 +-
 src/physics/utils/physics_column_type.F90 | 2 +-
 src/physics/utils/physics_grid.F90        | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90
index 4e621f63..ba9028bb 100644
--- a/src/dynamics/none/dyn_grid.F90
+++ b/src/dynamics/none/dyn_grid.F90
@@ -59,7 +59,7 @@ subroutine model_grid_init()
       use pio,              only: PIO_BCAST_ERROR, pio_seterrorhandling
       use pio,              only: pio_get_var, pio_freedecomp
       use pio,              only: pio_read_darray
-      use spmd_utils,       only: npes, iam, masterprocid, mpicom
+      use spmd_utils,       only: npes, iam
       use cam_pio_utils,    only: cam_pio_handle_error, cam_pio_find_var
       use cam_pio_utils,    only: cam_pio_var_info, pio_subsystem
       use cam_pio_utils,    only: cam_pio_newdecomp
diff --git a/src/physics/utils/physics_column_type.F90 b/src/physics/utils/physics_column_type.F90
index c4f5af93..3d058f69 100644
--- a/src/physics/utils/physics_column_type.F90
+++ b/src/physics/utils/physics_column_type.F90
@@ -1,6 +1,6 @@
 module physics_column_type
 
-   use shr_kind_mod, only: r8 => shr_kind_r8, shr_kind_cl
+   use shr_kind_mod, only: r8 => shr_kind_r8
    use ccpp_kinds,   only: kind_phys
 
 
diff --git a/src/physics/utils/physics_grid.F90 b/src/physics/utils/physics_grid.F90
index 264a662a..d477a432 100644
--- a/src/physics/utils/physics_grid.F90
+++ b/src/physics/utils/physics_grid.F90
@@ -65,7 +65,7 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
       use cam_abortutils,   only: endrun
       use spmd_utils,       only: npes, mpicom
       use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register
-      use cam_grid_support, only: iMap, hclen => max_hcoordname_len
+      use cam_grid_support, only: iMap
       use cam_grid_support, only: horiz_coord_t, horiz_coord_create
       use cam_grid_support, only: cam_grid_attribute_copy, cam_grid_attr_exists
 

From 14f783c41cbe0c4b0ac16bd6229f3f49db59a7bf Mon Sep 17 00:00:00 2001
From: Jesse Nusbaumer <nusbaume@ucar.edu>
Date: Fri, 15 Jan 2021 14:00:56 -0700
Subject: [PATCH 4/8] Address reviewer comments, and replace local string utils
 with CIME versions.

---
 src/dynamics/none/dyn_grid.F90            | 60 ++++++---------
 src/physics/utils/physics_column_type.F90 | 26 ++++---
 src/physics/utils/physics_grid.F90        | 40 +++++-----
 src/utils/string_utils.F90                | 89 ++++++-----------------
 4 files changed, 82 insertions(+), 133 deletions(-)

diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90
index ba9028bb..e769a79b 100644
--- a/src/dynamics/none/dyn_grid.F90
+++ b/src/dynamics/none/dyn_grid.F90
@@ -5,6 +5,7 @@ module dyn_grid
    use cam_logfile,         only: iulog, debug_output
    use spmd_utils,          only: masterproc
    use physics_column_type, only: physics_column_t
+   use string_utils,        only: to_str
 
    implicit none
    private
@@ -266,21 +267,18 @@ subroutine model_grid_init()
          if (is_degrees) then
             allocate(local_lons_deg(num_lons), stat=iret)
             if (iret /= 0) then
-               write(errormsg, *) &
-                  subname//': allocate local_lons_deg(num_lons) failed with stat: ',iret
-               call endrun(errormsg)
+               call endrun(subname//': allocate local_lons_deg(num_lons) failed with stat: '//&
+                           to_str(iret))
             end if
             allocate(local_lats_deg(kount(1)), stat=iret)
             if (iret /= 0) then
-               write(errormsg, *) &
-                  subname//': allocate local_lats_deg(kount) failed with stat: ',iret
-               call endrun(errormsg)
+               call endrun(subname//': allocate local_lats_deg(kount) failed with stat: '//&
+                           to_str(iret))
             end if
             allocate(temp_arr(num_lats), stat=iret)
             if (iret /= 0) then
-               write(errormsg, *) &
-                  subname//': allocate temp_arr failed with stat: ',iret
-               call endrun(errormsg)
+               call endrun(subname//': allocate temp_arr failed with stat: '//&
+                           to_str(iret))
             end if
             iret = pio_get_var(fh_ini, lat_vardesc, (/ 1 /), (/ num_lats /),  &
                  temp_arr)
@@ -304,21 +302,17 @@ subroutine model_grid_init()
          if (is_degrees) then
             allocate(local_lats_deg(num_local_columns), stat=iret)
             if (iret /= 0) then
-               write(errormsg, *) &
-                  subname//': allocate local_lats_deg(columns) failed with stat: ',iret
-               call endrun(errormsg)
+               call endrun(subname//': allocate local_lats_deg(columns) failed with stat: '//&
+                           to_str(iret))
             end if
             allocate(local_lons_deg(num_local_columns), stat=iret)
             if (iret /= 0) then
-               write(errormsg, *) &
-                  subname//': allocate local_lats_deg(columns) failed with stat: ',iret
-               call endrun(errormsg)
+               call endrun(subname//': allocate local_lons_deg(columns) failed with stat: '//&
+                           to_str(iret))
             end if
             allocate(ldof(num_local_columns), stat=iret)
             if (iret /= 0) then
-               write(errormsg, *) &
-                  subname//': allocate ldof failed with stat: ',iret
-               call endrun(errormsg)
+               call endrun(subname//': allocate ldof failed with stat: '//to_str(iret))
             end if
             ldof = 0_iMap
             do lindex = 1, num_local_columns
@@ -326,9 +320,8 @@ subroutine model_grid_init()
             end do
             allocate(iodesc, stat=iret)
             if (iret /= 0) then
-               write(errormsg, *) &
-                  subname//': allocate iodesc failed with stat: ',iret
-               call endrun(errormsg)
+               call endrun(subname//': allocate iodesc failed with stat: '//&
+                           to_str(iret))
             end if
             call cam_pio_newdecomp(iodesc, (/ num_global_columns /), ldof,    &
                  PIO_DOUBLE)
@@ -358,9 +351,8 @@ subroutine model_grid_init()
          if ((num_lats > 1) .and. (dimlens(1) == num_lats)) then
             allocate(local_areas(num_lats), stat=iret)
             if (iret /= 0) then
-               write(errormsg, *) &
-                  subname//': allocate local_areas(num_lats) failed with stat: ',iret
-               call endrun(errormsg)
+               call endrun(subname//': allocate local_areas(num_lats) failed with stat: '//&
+                           to_str(iret))
             end if
             start(1) = 1
             kount(1) = num_lats
@@ -370,9 +362,8 @@ subroutine model_grid_init()
          else if (dimlens(1) == num_global_columns) then
             allocate(local_areas(num_local_columns), stat=iret)
             if (iret /= 0) then
-               write(errormsg, *) &
-                  subname//': allocate local_areas(columns) failed with stat: ',iret
-               call endrun(errormsg)
+               call endrun(subname//': allocate local_areas(columns) failed with stat: '//&
+                           to_str(iret))
             end if
             call pio_read_darray(fh_ini, vardesc, iodesc, local_areas, iret)
             call cam_pio_handle_error(iret, subname//': Unable to read areas')
@@ -399,9 +390,8 @@ subroutine model_grid_init()
       ! Set dynamics grid attributes
       allocate(grid_attribute_names(0), stat=iret)
       if (iret /= 0) then
-         write(errormsg, *) &
-            subname//': allocate grid_attribute_names failed with stat: ',iret
-         call endrun(errormsg)
+         call endrun(subname//': allocate grid_attribute_names failed with stat: '//&
+                     to_str(iret))
       end if
 
       ! Initialize physics grid decomposition:
@@ -446,9 +436,8 @@ subroutine set_dyn_col_values()
       if (.not.allocated(dyn_columns)) then
          allocate(dyn_columns(num_local_columns), stat=ierr)
          if (ierr /= 0) then
-            write(emsg, *) &
-               subname//': allocate dyn_columns failed with stat: ',ierr
-            call endrun(emsg)
+            call endrun(subname//': allocate dyn_columns failed with stat: '//&
+                        to_str(ierr))
          end if
       end if
 
@@ -507,9 +496,8 @@ subroutine set_dyn_col_values()
          !    as in the dynamics block structure
          allocate(dyn_columns(lindex)%dyn_block_index(1), stat=ierr)
          if (ierr /= 0) then
-            write(emsg, *) &
-               subname//': allocate dyn_block_index failed with stat: ',ierr
-            call endrun(emsg)
+            call endrun(subname//': allocate dyn_block_index failed with stat: '//&
+                        to_str(ierr))
          end if
 
          dyn_columns(lindex)%dyn_block_index(1) = lindex
diff --git a/src/physics/utils/physics_column_type.F90 b/src/physics/utils/physics_column_type.F90
index 3d058f69..ddf67f26 100644
--- a/src/physics/utils/physics_column_type.F90
+++ b/src/physics/utils/physics_column_type.F90
@@ -49,6 +49,7 @@ module physics_column_type
    subroutine copy_phys_col(phys_col_out, phys_col_in)
 
       use cam_abortutils, only: endrun
+      use string_utils,   only: to_str
 
       ! Copy all of the values from one physics_column type
       ! structure to another.
@@ -84,21 +85,28 @@ subroutine copy_phys_col(phys_col_out, phys_col_in)
 
       ! Dynamics blocks
       if (allocated(phys_col_in%dyn_block_index)) then
-         ! De-allocate output block indices if already allocated:
+         ! De-allocate output block indices allocated to incorrect size:
          if (allocated(phys_col_out%dyn_block_index)) then
-            deallocate(phys_col_out%dyn_block_index)
+            if (size(phys_col_out%dyn_block_index) /= &
+               size(phys_col_in%dyn_block_index)) then
+               deallocate(phys_col_out%dyn_block_index)
+            end if
          end if
 
-         ! Allocate output to match size of input:
-         allocate(phys_col_out%dyn_block_index(size(phys_col_in%dyn_block_index)), &
-                  stat=istat)
-         if (istat /= 0) then
-           write(emsg, *) &
-              subname//': allocate dyn_block_index failed with stat: ',istat
-           call endrun(emsg)
+         ! If necessary, allocate output to match size of input:
+         if (.not. allocated(phys_col_out%dyn_block_index)) then
+            allocate(phys_col_out%dyn_block_index(size(phys_col_in%dyn_block_index)), &
+                     stat=istat)
+            if (istat /= 0) then
+               call endrun(subname//': allocate dyn_block_index failed with stat: '//&
+                           to_str(istat))
+            end if
          end if
 
          phys_col_out%dyn_block_index(:) = phys_col_in%dyn_block_index(:)
+      else if (allocated(phys_col_out%dyn_block_index)) then
+         ! De-allocate output array if input array has not been allocated:
+         deallocate(phys_col_out%dyn_block_index)
       end if
 
       ! Physics decomposition
diff --git a/src/physics/utils/physics_grid.F90 b/src/physics/utils/physics_grid.F90
index d477a432..4cd28f94 100644
--- a/src/physics/utils/physics_grid.F90
+++ b/src/physics/utils/physics_grid.F90
@@ -64,6 +64,7 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
       use mpi,              only: MPI_INTEGER, MPI_MIN
       use cam_abortutils,   only: endrun
       use spmd_utils,       only: npes, mpicom
+      use string_utils,     only: to_str
       use cam_grid_support, only: cam_grid_register, cam_grid_attribute_register
       use cam_grid_support, only: iMap
       use cam_grid_support, only: horiz_coord_t, horiz_coord_create
@@ -133,18 +134,18 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
          index_top_interface = index_top_layer + 1
       end if
 
-      !Calculate total number of physics columns:
+      ! Calculate total number of physics columns:
       num_global_phys_cols = hdim1_d * hdim2_d
 
-      ! Set columns_on_task and allocate phys_columns if
-      ! not already allocated:
+      ! Calculate number of columns on tasks:
+      columns_on_task = size(dyn_columns)
+
+      ! Set allocate phys_columns if not already allocated:
       if (.not. allocated(phys_columns)) then
-         columns_on_task = size(dyn_columns)
          allocate(phys_columns(columns_on_task), stat=ierr)
          if (ierr /= 0) then
-            write(emsg, *) &
-               subname//': allocate phys_columns failed with stat: ',ierr
-            call endrun(emsg)
+            call endrun(subname//': allocate phys_columns failed with stat: '//&
+                        to_str(ierr))
          end if
       end if
 
@@ -170,24 +171,21 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
          allocate(grid_map(4, columns_on_task), stat=ierr)
       end if
       if (ierr /= 0) then
-         write(emsg, *) &
-            subname//': allocate grid_map failed with stat: ',ierr
-         call endrun(emsg)
+         call endrun(subname//': allocate grid_map failed with stat: '//&
+                     to_str(ierr))
       end if
       grid_map = 0
 
       allocate(latvals(size(grid_map, 2)), stat=ierr)
       if (ierr /= 0) then
-         write(emsg, *) &
-            subname//': allocate latvals failed with stat: ',ierr
-         call endrun(emsg)
+         call endrun(subname//': allocate latvals failed with stat: '//&
+                     to_str(ierr))
       end if
 
       allocate(lonvals(size(grid_map, 2)), stat=ierr)
       if (ierr /= 0) then
-         write(emsg, *) &
-            subname//': allocate lonvals failed with stat: ',ierr
-         call endrun(emsg)
+         call endrun(subname//': allocate lonvals failed with stat: '//&
+                     to_str(ierr))
       end if
 
       lonmin = 1000.0_r8 ! Out of longitude range
@@ -228,9 +226,8 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
       else
          allocate(coord_map(size(grid_map, 2)), stat=ierr)
          if (ierr /= 0) then
-            write(emsg, *) &
-               subname//': allocate coord_map failed with stat: ',ierr
-            call endrun(emsg)
+            call endrun(subname//': allocate coord_map failed with stat: '//&
+                        to_str(ierr))
          end if
 
          ! We need a global minimum longitude and latitude
@@ -281,9 +278,8 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
          !   grid is only supported for unstructured grids).
          allocate(area_d(size(grid_map, 2)), stat=ierr)
          if (ierr /= 0) then
-            write(emsg, *) &
-               subname//': allocate area_d failed with stat: ',ierr
-            call endrun(emsg)
+            call endrun(subname//': allocate area_d failed with stat: '//&
+                        to_str(ierr))
          end if
 
          do col_index = 1, columns_on_task
diff --git a/src/utils/string_utils.F90 b/src/utils/string_utils.F90
index b57c708d..7751da69 100644
--- a/src/utils/string_utils.F90
+++ b/src/utils/string_utils.F90
@@ -1,5 +1,7 @@
 module string_utils
 
+   use shr_string_mod, only: to_upper => shr_string_toUpper
+   use shr_string_mod, only: to_lower => shr_string_toLower
 
    implicit none
    private
@@ -10,6 +12,7 @@ module string_utils
    public :: to_lower         ! Convert character string to lower case
    public :: increment_string ! increments a string
    public :: last_sig_char    ! Position of last significant character in string
+   public :: to_str           ! convert integer to left justified string
 
    ! Private module variables
    integer, parameter :: lower_to_upper = iachar("A") - iachar("a")
@@ -17,72 +20,6 @@ module string_utils
 
 CONTAINS
 
-   function to_upper(str)
-
-      !-----------------------------------------------------------------------
-      ! Purpose:
-      ! Convert character string to upper case.
-      !
-      ! Method:
-      ! Use achar and iachar intrinsics to ensure use of ascii collating seq.
-      !
-      ! Author:  B. Eaton, July 2001
-      !
-      !-----------------------------------------------------------------------
-      character(len=*), intent(in) :: str ! String to convert to upper case
-      character(len=len(str))      :: to_upper
-
-      ! Local variables
-
-      integer          :: ind            ! Index
-      integer          :: aseq           ! ascii collating sequence
-      character(len=1) :: ctmp           ! Character temporary
-      !-----------------------------------------------------------------------
-      to_upper = ''
-      do ind = 1, len_trim(str)
-         ctmp = str(ind:ind)
-         aseq = iachar(ctmp)
-         if ((aseq >= iachar("a")) .and. (aseq <= iachar("z"))) then
-            ctmp = achar(aseq + lower_to_upper)
-         end if
-         to_upper(ind:ind) = ctmp
-      end do
-
-   end function to_upper
-
-   function to_lower(str)
-
-      !-----------------------------------------------------------------------
-      ! Purpose:
-      ! Convert character string to lower case.
-      !
-      ! Method:
-      ! Use achar and iachar intrinsics to ensure use of ascii collating seq.
-      !
-      ! Author:  B. Eaton, July 2001
-      !
-      !-----------------------------------------------------------------------
-      character(len=*), intent(in) :: str ! String to convert to lower case
-      character(len=len(str))      :: to_lower
-
-      ! Local variables
-
-      integer          :: ind            ! Index
-      integer          :: aseq           ! ascii collating sequence
-      character(len=1) :: ctmp           ! Character temporary
-      !-----------------------------------------------------------------------
-
-      do ind = 1, len(str)
-         ctmp = str(ind:ind)
-         aseq = iachar(ctmp)
-         if ((aseq >= iachar("A")) .and. (aseq <= iachar("Z"))) then
-            ctmp = achar(aseq + upper_to_lower)
-         end if
-         to_lower(ind:ind) = ctmp
-      end do
-
-   end function to_lower
-
    integer function increment_string(str, increment)
       !-----------------------------------------------------------------------
       ! 	... Increment a string whose ending characters are digits.
@@ -155,6 +92,8 @@ integer function increment_string(str, increment)
 
    end function increment_string
 
+!=========================================================================================
+
    integer function last_index(cstr)
       !-----------------------------------------------------------------------
       ! 	... Position of last non-digit in the first input token.
@@ -192,6 +131,8 @@ integer function last_index(cstr)
 
    end function last_index
 
+!=========================================================================================
+
    integer function last_sig_char(cstr)
       !-----------------------------------------------------------------------
       ! 	... Position of last significant character in string.
@@ -227,4 +168,20 @@ integer function last_sig_char(cstr)
 
    end function last_sig_char
 
+!=========================================================================================
+
+character(len=10) function to_str(n)
+
+   ! return default integer as a left justified string
+
+   ! arguments
+   integer, intent(in) :: n
+   !----------------------------------------------------------------------------
+
+   write(to_str,'(i0)') n
+
+end function to_str
+
+!=========================================================================================
+
 end module string_utils

From d1915b2bd280302ae449b3fb24db2472a4271e6b Mon Sep 17 00:00:00 2001
From: Jesse Nusbaumer <nusbaume@ucar.edu>
Date: Fri, 15 Jan 2021 15:06:26 -0700
Subject: [PATCH 5/8] Fix additional NAG warnings.

---
 src/dynamics/none/dyn_grid.F90            | 1 -
 src/physics/utils/physics_column_type.F90 | 1 -
 src/physics/utils/physics_grid.F90        | 1 -
 3 files changed, 3 deletions(-)

diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90
index e769a79b..bd5f221a 100644
--- a/src/dynamics/none/dyn_grid.F90
+++ b/src/dynamics/none/dyn_grid.F90
@@ -426,7 +426,6 @@ subroutine set_dyn_col_values()
       integer                         :: lat_index, lat1
       integer                         :: lon_index
       integer                         :: ierr
-      character(len=128)              :: emsg
 
       real(r8),         parameter     :: radtodeg = 180.0_r8 / SHR_CONST_PI
       real(r8),         parameter     :: degtorad = SHR_CONST_PI / 180.0_r8
diff --git a/src/physics/utils/physics_column_type.F90 b/src/physics/utils/physics_column_type.F90
index ddf67f26..7ed6bd3a 100644
--- a/src/physics/utils/physics_column_type.F90
+++ b/src/physics/utils/physics_column_type.F90
@@ -62,7 +62,6 @@ subroutine copy_phys_col(phys_col_out, phys_col_in)
       integer :: istat
 
       character(len=*),  parameter :: subname = 'copy_phys_col'
-      character(len=128)           :: emsg
 
       ! Copy values from input array to output array:
 
diff --git a/src/physics/utils/physics_grid.F90 b/src/physics/utils/physics_grid.F90
index 4cd28f94..f2f8fad7 100644
--- a/src/physics/utils/physics_grid.F90
+++ b/src/physics/utils/physics_grid.F90
@@ -99,7 +99,6 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
       logical                             :: unstructured
       real(r8)                            :: temp ! For MPI
       integer                             :: ierr ! For MPI
-      character(len=128)                  :: emsg
 
       character(len=*),  parameter :: subname = 'phys_grid_init'
 

From 4291fcb0bf66f83d1a340a18b9f556a7feba23df Mon Sep 17 00:00:00 2001
From: Jesse Nusbaumer <nusbaume@ucar.edu>
Date: Tue, 19 Jan 2021 11:08:44 -0700
Subject: [PATCH 6/8] Address additional code review comments.

---
 src/physics/utils/physics_column_type.F90 |  2 +-
 src/physics/utils/physics_grid.F90        | 17 ++++++++++-------
 2 files changed, 11 insertions(+), 8 deletions(-)

diff --git a/src/physics/utils/physics_column_type.F90 b/src/physics/utils/physics_column_type.F90
index 7ed6bd3a..92060450 100644
--- a/src/physics/utils/physics_column_type.F90
+++ b/src/physics/utils/physics_column_type.F90
@@ -84,7 +84,7 @@ subroutine copy_phys_col(phys_col_out, phys_col_in)
 
       ! Dynamics blocks
       if (allocated(phys_col_in%dyn_block_index)) then
-         ! De-allocate output block indices allocated to incorrect size:
+         ! De-allocate output block indices if allocated to incorrect size:
          if (allocated(phys_col_out%dyn_block_index)) then
             if (size(phys_col_out%dyn_block_index) /= &
                size(phys_col_in%dyn_block_index)) then
diff --git a/src/physics/utils/physics_grid.F90 b/src/physics/utils/physics_grid.F90
index f2f8fad7..048c5706 100644
--- a/src/physics/utils/physics_grid.F90
+++ b/src/physics/utils/physics_grid.F90
@@ -109,6 +109,11 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
       nullify(lon_coord)
       nullify(area_d)
 
+      ! Check that the physics grid is not already initialized:
+      if (phys_grid_initialized) then
+         call endrun(subname//": Physics grid is already initialized.")
+      end if
+
       call t_adj_detailf(-2)
       call t_startf("phys_grid_init")
 
@@ -139,13 +144,11 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
       ! Calculate number of columns on tasks:
       columns_on_task = size(dyn_columns)
 
-      ! Set allocate phys_columns if not already allocated:
-      if (.not. allocated(phys_columns)) then
-         allocate(phys_columns(columns_on_task), stat=ierr)
-         if (ierr /= 0) then
-            call endrun(subname//': allocate phys_columns failed with stat: '//&
-                        to_str(ierr))
-         end if
+      ! Allocate phys_columns:
+      allocate(phys_columns(columns_on_task), stat=ierr)
+      if (ierr /= 0) then
+         call endrun(subname//': allocate phys_columns failed with stat: '//&
+                     to_str(ierr))
       end if
 
       ! Set column index bounds:

From 75ba8f9ae895bc3cbdb50bae2a1d4e33f8f80796 Mon Sep 17 00:00:00 2001
From: Jesse Nusbaumer <nusbaume@ucar.edu>
Date: Thu, 21 Jan 2021 13:40:41 -0700
Subject: [PATCH 7/8] Modify allocation code as recommended by code reviewer.

---
 src/dynamics/none/dyn_grid.F90            | 42 ++++++++++++-----------
 src/physics/utils/physics_column_type.F90 |  4 +--
 src/physics/utils/physics_grid.F90        | 38 +++++++++++---------
 3 files changed, 45 insertions(+), 39 deletions(-)

diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90
index bd5f221a..cf877c19 100644
--- a/src/dynamics/none/dyn_grid.F90
+++ b/src/dynamics/none/dyn_grid.F90
@@ -277,7 +277,7 @@ subroutine model_grid_init()
             end if
             allocate(temp_arr(num_lats), stat=iret)
             if (iret /= 0) then
-               call endrun(subname//': allocate temp_arr failed with stat: '//&
+               call endrun(subname//': allocate temp_arr(num_lats) failed with stat: '//&
                            to_str(iret))
             end if
             iret = pio_get_var(fh_ini, lat_vardesc, (/ 1 /), (/ num_lats /),  &
@@ -302,17 +302,18 @@ subroutine model_grid_init()
          if (is_degrees) then
             allocate(local_lats_deg(num_local_columns), stat=iret)
             if (iret /= 0) then
-               call endrun(subname//': allocate local_lats_deg(columns) failed with stat: '//&
-                           to_str(iret))
+               call endrun(subname//': allocate local_lats_deg(num_local_columns) '//&
+                           'failed with stat: '//to_str(iret))
             end if
             allocate(local_lons_deg(num_local_columns), stat=iret)
             if (iret /= 0) then
-               call endrun(subname//': allocate local_lons_deg(columns) failed with stat: '//&
-                           to_str(iret))
+               call endrun(subname//': allocate local_lons_deg(num_local_columns) '//&
+                           'failed with stat: '//to_str(iret))
             end if
             allocate(ldof(num_local_columns), stat=iret)
             if (iret /= 0) then
-               call endrun(subname//': allocate ldof failed with stat: '//to_str(iret))
+               call endrun(subname//': allocate ldof(num_local_columns) '//&
+                           'failed with stat: '//to_str(iret))
             end if
             ldof = 0_iMap
             do lindex = 1, num_local_columns
@@ -362,8 +363,8 @@ subroutine model_grid_init()
          else if (dimlens(1) == num_global_columns) then
             allocate(local_areas(num_local_columns), stat=iret)
             if (iret /= 0) then
-               call endrun(subname//': allocate local_areas(columns) failed with stat: '//&
-                           to_str(iret))
+               call endrun(subname//': allocate local_areas(num_local_columns) '//&
+                           'failed with stat: '//to_str(iret))
             end if
             call pio_read_darray(fh_ini, vardesc, iodesc, local_areas, iret)
             call cam_pio_handle_error(iret, subname//': Unable to read areas')
@@ -384,13 +385,22 @@ subroutine model_grid_init()
       ! Back to old error handling
       call pio_seterrorhandling(fh_ini, err_handling)
 
+      ! Allocate dyn_columns structure if not already allocated:
+      if (.not.allocated(dyn_columns)) then
+         allocate(dyn_columns(num_local_columns), stat=iret)
+         if (iret /= 0) then
+            call endrun(subname//': allocate dyn_columns(num_local_columns) '//&
+                        'failed with stat: '//to_str(iret))
+         end if
+      end if
+
       ! Set dyn_columns values:
       call set_dyn_col_values()
 
-      ! Set dynamics grid attributes
+      ! The null dycore has no grid attributes, so allocate to size zero.
       allocate(grid_attribute_names(0), stat=iret)
       if (iret /= 0) then
-         call endrun(subname//': allocate grid_attribute_names failed with stat: '//&
+         call endrun(subname//': allocate grid_attribute_names(0) failed with stat: '//&
                      to_str(iret))
       end if
 
@@ -431,15 +441,6 @@ subroutine set_dyn_col_values()
       real(r8),         parameter     :: degtorad = SHR_CONST_PI / 180.0_r8
       character(len=*), parameter     :: subname = 'set_dyn_col_values'
 
-      ! Allocate dyn_columns structure if not already allocated:
-      if (.not.allocated(dyn_columns)) then
-         allocate(dyn_columns(num_local_columns), stat=ierr)
-         if (ierr /= 0) then
-            call endrun(subname//': allocate dyn_columns failed with stat: '//&
-                        to_str(ierr))
-         end if
-      end if
-
       ! Calculate dyn_columns variable values:
       lat1 = global_col_offset / num_lons
       do lindex = 1, num_local_columns
@@ -495,7 +496,8 @@ subroutine set_dyn_col_values()
          !    as in the dynamics block structure
          allocate(dyn_columns(lindex)%dyn_block_index(1), stat=ierr)
          if (ierr /= 0) then
-            call endrun(subname//': allocate dyn_block_index failed with stat: '//&
+            call endrun(subname//': allocate dyn_columns(lindex)%'//&
+                        'dyn_block_index(1) failed with stat: '//&
                         to_str(ierr))
          end if
 
diff --git a/src/physics/utils/physics_column_type.F90 b/src/physics/utils/physics_column_type.F90
index 92060450..697e8248 100644
--- a/src/physics/utils/physics_column_type.F90
+++ b/src/physics/utils/physics_column_type.F90
@@ -97,8 +97,8 @@ subroutine copy_phys_col(phys_col_out, phys_col_in)
             allocate(phys_col_out%dyn_block_index(size(phys_col_in%dyn_block_index)), &
                      stat=istat)
             if (istat /= 0) then
-               call endrun(subname//': allocate dyn_block_index failed with stat: '//&
-                           to_str(istat))
+               call endrun(subname//': allocate phys_col_out%dyn_block_index failed '//&
+                           'with stat: '//to_str(istat))
             end if
          end if
 
diff --git a/src/physics/utils/physics_grid.F90 b/src/physics/utils/physics_grid.F90
index 048c5706..ade984b9 100644
--- a/src/physics/utils/physics_grid.F90
+++ b/src/physics/utils/physics_grid.F90
@@ -147,8 +147,8 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
       ! Allocate phys_columns:
       allocate(phys_columns(columns_on_task), stat=ierr)
       if (ierr /= 0) then
-         call endrun(subname//': allocate phys_columns failed with stat: '//&
-                     to_str(ierr))
+         call endrun(subname//': allocate phys_columns(columns_on_task) '//&
+                     ' failed with stat: '//to_str(ierr))
       end if
 
       ! Set column index bounds:
@@ -169,25 +169,29 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
       ! unstructured
       if (unstructured) then
          allocate(grid_map(3, columns_on_task), stat=ierr)
+         if (ierr /= 0) then
+            call endrun(subname//': allocate grid_map(3, columns_on_task) '//&
+                        'failed with stat: '//to_str(ierr))
+         end if
       else
          allocate(grid_map(4, columns_on_task), stat=ierr)
-      end if
-      if (ierr /= 0) then
-         call endrun(subname//': allocate grid_map failed with stat: '//&
-                     to_str(ierr))
+         if (ierr /= 0) then
+            call endrun(subname//': allocate grid_map(4, columns_on_task) '//&
+                        'failed with stat: '//to_str(ierr))
+         end if
       end if
       grid_map = 0
 
-      allocate(latvals(size(grid_map, 2)), stat=ierr)
+      allocate(latvals(columns_on_task), stat=ierr)
       if (ierr /= 0) then
-         call endrun(subname//': allocate latvals failed with stat: '//&
-                     to_str(ierr))
+         call endrun(subname//': allocate latvals(columns_on_task) '//&
+                     'failed with stat: '//to_str(ierr))
       end if
 
-      allocate(lonvals(size(grid_map, 2)), stat=ierr)
+      allocate(lonvals(columns_on_task), stat=ierr)
       if (ierr /= 0) then
-         call endrun(subname//': allocate lonvals failed with stat: '//&
-                     to_str(ierr))
+         call endrun(subname//': allocate lonvals(columns_on_task) '//&
+                     'failed with stat: '//to_str(ierr))
       end if
 
       lonmin = 1000.0_r8 ! Out of longitude range
@@ -226,10 +230,10 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
               'latitude', 'degrees_north', 1, size(latvals), latvals,         &
               map=grid_map(3,:))
       else
-         allocate(coord_map(size(grid_map, 2)), stat=ierr)
+         allocate(coord_map(columns_on_task), stat=ierr)
          if (ierr /= 0) then
-            call endrun(subname//': allocate coord_map failed with stat: '//&
-                        to_str(ierr))
+            call endrun(subname//': allocate coord_map(columns_on_task) '//&
+                        'failed with stat: '//to_str(ierr))
          end if
 
          ! We need a global minimum longitude and latitude
@@ -278,9 +282,9 @@ subroutine phys_grid_init(hdim1_d_in, hdim2_d_in, pver_in, dycore_name_in, &
          !   from the dycore (i.e., physics and dynamics are on different
          !   grids), create that attribute here (Note, a separate physics
          !   grid is only supported for unstructured grids).
-         allocate(area_d(size(grid_map, 2)), stat=ierr)
+         allocate(area_d(columns_on_task), stat=ierr)
          if (ierr /= 0) then
-            call endrun(subname//': allocate area_d failed with stat: '//&
+            call endrun(subname//': allocate area_d(columns_on_task) failed with stat: '//&
                         to_str(ierr))
          end if
 

From c3d87a95f09dfe8744f44b5473190eb39ba047f8 Mon Sep 17 00:00:00 2001
From: Jesse Nusbaumer <nusbaume@ucar.edu>
Date: Thu, 21 Jan 2021 14:18:41 -0700
Subject: [PATCH 8/8] Print actual index number in allocation error message.

---
 src/dynamics/none/dyn_grid.F90 | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/dynamics/none/dyn_grid.F90 b/src/dynamics/none/dyn_grid.F90
index cf877c19..4f6247c9 100644
--- a/src/dynamics/none/dyn_grid.F90
+++ b/src/dynamics/none/dyn_grid.F90
@@ -496,9 +496,9 @@ subroutine set_dyn_col_values()
          !    as in the dynamics block structure
          allocate(dyn_columns(lindex)%dyn_block_index(1), stat=ierr)
          if (ierr /= 0) then
-            call endrun(subname//': allocate dyn_columns(lindex)%'//&
-                        'dyn_block_index(1) failed with stat: '//&
-                        to_str(ierr))
+            call endrun(subname//': allocate dyn_columns('//&
+                        to_str(lindex)//')%dyn_block_index(1)'//&
+                        ' failed with stat: '//to_str(ierr))
          end if
 
          dyn_columns(lindex)%dyn_block_index(1) = lindex