Add ability to read snapshot variable with extra dimension

src/data/write_init_files.py

@@ -789,20 +789,34 @@ def get_dimension_info(hvar):
         # based on constituent name.
         # This case will be handled separately.
         legal_dims = True
-    elif (ldims > 2) or ((ldims > 1) and (not levnm)):
+    elif (ldims > 3) or ((ldims > 1) and (not levnm)):
         # The regular case where the second dimension must be vertical,
-        # and higher dimensions are unsupported.
+        # and dimensions greater than three are unsupported.
         legal_dims = False
         unsupp = []
-        for dim in dims:
-            if ((not is_horizontal_dimension(dim)) and
-                (not is_vertical_dimension(dim))):
-                if dim[0:18] == "ccpp_constant_one:":
-                    rdim = dim[18:]
+        vert_idx = -1
+        for dim_idx, dim in enumerate(dims):
+            # Skip the first dimension if it is horizontal:
+            if (not is_horizontal_dimension(dim) or dim_idx > 0):
+                # Save (and skip) the vertical dimension index,
+                # which must be the second dimension present:
+                if (is_vertical_dimension(dim) and dim_idx == 1):
+                    vert_idx = dim_idx
                 else:
-                    rdim = dim
+                    # An extra, non-horizontal dimension is allowed
+                    # immediately after the vertical dimension,
+                    # so check if this is not the case:
+                    if (vert_idx < 0) or (dim_idx != (vert_idx + 1)):
+                        if dim[0:18] == "ccpp_constant_one:":
+                            rdim = dim[18:]
+                        else:
+                            rdim = dim
+                        # end if
+                        #Update dimension string with dimension index:
+                        rdim = rdim + f" (dimension {dim_idx+1})"
+                        unsupp.append(rdim)
+                    # end if
                 # end if
-                unsupp.append(rdim)
             # end if
         # end for
         if len(unsupp) > 1:
@@ -949,8 +963,8 @@ def write_phys_read_subroutine(outfile, host_dict, host_vars, host_imports,
                                    "no_exist_idx", "init_mark_idx",
                                    "prot_no_init_idx", "const_idx",
                                    "read_constituent_dimensioned_field"]],
-                 ["cam_ccpp_cap", ["ccpp_physics_suite_variables", 
-                                   "cam_constituents_array", 
+                 ["cam_ccpp_cap", ["ccpp_physics_suite_variables",
+                                   "cam_constituents_array",
                                    "cam_model_const_properties"]],
                  ["ccpp_kinds", ["kind_phys"]],
                  [phys_check_fname_str, ["phys_var_num", "phys_var_stdnames",
@@ -1149,7 +1163,7 @@ def write_phys_read_subroutine(outfile, host_dict, host_vars, host_imports,
     # End suite loop:
     outfile.write(" end do !CCPP suites", 2)
     outfile.blank_line()
-    
+
     # Read in constituent data
     outfile.comment("Read in constituent variables if not using init variables", 2)
     outfile.write("field_data_ptr => cam_constituents_array()", 2)

src/utils/cam_field_read.F90

@@ -107,7 +107,7 @@ subroutine print_input_field_info(dimlens, ndims, min_ndims, max_ndims,    &
               trim(varname), ', ', ndims
          call safe_endrun(subname//trim(errormsg))
       else if (num_bounds < 1) then
-         call safe_endrun(subname//': too few dimension boundss for '//trim(varname))
+         call safe_endrun(subname//': too few dimension bounds for '//trim(varname))
       else if (num_bounds > 3) then
          write(errormsg, '(3a,i0)') ': too many dimension bounds for, ',      &
               trim(varname), ', ', num_bounds
@@ -186,6 +186,7 @@ subroutine infld_real8_1d(varname, ncid, field, readvar, gridname,         &
       ! Local array, <field> is 1D
       !
 
+      use shr_kind_mod,     only: cl=>shr_kind_cl
       use pio,              only: pio_read_darray, PIO_NOERR
       use pio,              only: PIO_MAX_NAME, pio_inq_dimname
       use cam_grid_support, only: cam_grid_get_decomp, cam_grid_is_unstructured
@@ -234,7 +235,7 @@ subroutine infld_real8_1d(varname, ncid, field, readvar, gridname,         &
       integer                                   :: dim_bounds(2, 2)
       character(len=PIO_MAX_NAME)               :: dim1name, dim2name
       character(len=PIO_MAX_NAME)               :: tmpname
-      character(len=128)                        :: errormsg
+      character(len=cl)                         :: errormsg
       ! Logical to determine whether to print output messages
       logical                                   :: log_read_field
       logical                                   :: readvar_tmp
@@ -327,7 +328,9 @@ subroutine infld_real8_1d(varname, ncid, field, readvar, gridname,         &
          arraydimsize(1) = (dim_bounds(1,2) - dim_bounds(1,1) + 1)
          if (arraydimsize(1) /= size(field, 1)) then
             write(errormsg, '(4a,i0)') ': Mismatch between array bounds ',    &
-                 'and field size for ', trim(varname), ', dimension ', 1
+                 'and field size for ', trim(varname), ', dimension 1', &
+                 ', variable array bounds = ', arraydimsize(1), ', file field size = ', &
+                 size(field, 1)
             call safe_endrun(subname//errormsg)
          end if
 
@@ -409,6 +412,7 @@ subroutine infld_real8_2d(varname, ncid, field, readvar, gridname,         &
       ! Local array, <field> is 2D
       !
 
+      use shr_kind_mod,     only: cl=>shr_kind_cl
       use pio,              only: pio_read_darray, PIO_NOERR
       use pio,              only: PIO_MAX_NAME, pio_inq_dimname
       use cam_grid_support, only: cam_grid_get_decomp, cam_grid_is_unstructured
@@ -425,7 +429,7 @@ subroutine infld_real8_2d(varname, ncid, field, readvar, gridname,         &
       ! gridname: Name of variable's grid (default 'physgrid')
       character(len=*), optional, intent(in)  :: gridname
       integer,          optional, intent(in)  :: timelevel
-      ! dim3name: Name of vertical dimension, if field reprsents a 3D quantity
+      ! dim3name: Name of vertical dimension, if field represents a 3D quantity
       character(len=*), optional, intent(in)  :: dim3name
       ! dim3_bnds: Bounds of vertical dimension, if field is 3D
       integer,          optional, intent(in)  :: dim3_bnds(2)
@@ -463,7 +467,7 @@ subroutine infld_real8_2d(varname, ncid, field, readvar, gridname,         &
       character(len=PIO_MAX_NAME)             :: dim1name, dim2name
       character(len=PIO_MAX_NAME)             :: file_dnames(PIO_MAX_VAR_DIMS)
       character(len=PIO_MAX_NAME)             :: tmpname
-      character(len=128)                      :: errormsg
+      character(len=cl)                       :: errormsg
       ! Logical to determine whether to print output messages
       logical                                 :: log_read_field
       logical                                 :: readvar_tmp
@@ -575,7 +579,9 @@ subroutine infld_real8_2d(varname, ncid, field, readvar, gridname,         &
                write(errormsg, '(2(a,i0),3a,i0)')                             &
                     ': Mismatch between array size (', arraydimsize(jndex),   &
                     ') and field size (', size(field, jndex), ') for ',       &
-                    trim(varname), ', dimension = ', jndex
+                    trim(varname), ', dimension = ', jndex,                   &
+                    ', variable array bounds = ', arraydimsize(jndex),        &
+                    ', file field size = ', size(field, jndex)
                call safe_endrun(subname//errormsg)
             end if
          end do
@@ -664,22 +670,21 @@ subroutine infld_real8_2d(varname, ncid, field, readvar, gridname,         &
 
       readvar = readvar_tmp
 
-      return
-
    end subroutine infld_real8_2d
 
    !
    ! ROUTINE: infld_real8_3d
    !
    subroutine infld_real8_3d(varname, ncid, field, readvar, dim3name,         &
-        dim3_bnds, dim3_pos, gridname, timelevel, log_output, fillvalue)
+        dim3_bnds, dim3_pos, dim4_bnds, gridname, timelevel, log_output, fillvalue)
       !
       ! infld_real8_3d:
       ! Netcdf I/O of 8-byte real field from netCDF file
-      ! Field on file is 3D
+      ! Field on file is 3D or 4D
       ! Local array, <field> is 3D
       !
 
+      use shr_kind_mod,     only: cl=>shr_kind_cl
       use pio,              only: pio_read_darray, PIO_NOERR
       use pio,              only: PIO_MAX_NAME, pio_inq_dimname
       use cam_grid_support, only: cam_grid_get_decomp, cam_grid_is_unstructured
@@ -698,6 +703,8 @@ subroutine infld_real8_3d(varname, ncid, field, readvar, dim3name,         &
       ! dim3_bnds: Bounds of vertical dimension, if field is 3D
       integer,                    intent(in)  :: dim3_bnds(2)
       integer,          optional, intent(in)  :: dim3_pos
+      ! dim4_bnds: Bounds of 4th dimension, if field is 4D (e.g., size bins dimension for aerosol fields)
+      integer,          optional, intent(in)  :: dim4_bnds(2)
       ! gridname: Name of variable's grid (default 'physgrid')
       character(len=*), optional, intent(in)  :: gridname
       integer,          optional, intent(in)  :: timelevel
@@ -719,6 +726,8 @@ subroutine infld_real8_3d(varname, ncid, field, readvar, dim3name,         &
 
       ! ndims: The number of dimensions of the field in the file
       integer                                 :: ndims
+      ! pdims: The number of dimensions of the field in the file to print
+      integer                                 :: pdims
       ! target_ndims: The number of expected dimensions for field on file
       integer                                 :: target_ndims
       ! dimids: file variable dims
@@ -733,7 +742,7 @@ subroutine infld_real8_3d(varname, ncid, field, readvar, dim3name,         &
       character(len=PIO_MAX_NAME)             :: dim1name, dim2name
       character(len=PIO_MAX_NAME)             :: file_dnames(PIO_MAX_VAR_DIMS)
       character(len=PIO_MAX_NAME)             :: tmpname
-      character(len=128)                      :: errormsg
+      character(len=cl)                       :: errormsg
       ! Logical to determine whether to print output messages
       logical                                 :: log_read_field
       logical                                 :: readvar_tmp
@@ -769,7 +778,11 @@ subroutine infld_real8_3d(varname, ncid, field, readvar, dim3name,         &
       block_indexed = cam_grid_is_block_indexed(grid_id)
       ! Is this an unstructured grid (i.e., one column dimension on file)?
       unstruct = cam_grid_is_unstructured(grid_id)
-      target_ndims = num_target_dims(3, unstruct)
+      if (present(dim4_bnds)) then
+         target_ndims = num_target_dims(4, unstruct)
+      else
+         target_ndims = num_target_dims(3, unstruct)
+      end if   
       if ((debug_output > 0) .and. masterproc) then
          if (present(timelevel)) then
             write(errormsg, '(a,i0)') ', timelevel = ', timelevel
@@ -826,6 +839,21 @@ subroutine infld_real8_3d(varname, ncid, field, readvar, dim3name,         &
          dim_bounds(index,1) = dim3_bnds(1)
          dim_bounds(index,2) = dim3_bnds(2)
       end if
+
+      !If the variable has an extra (non-spatial) dimension as indicated
+      !by the presence of dim4_bnds, then set the dim bounds of the extra
+      !dimension as provided by the caller, with the assumption that it
+      !is always immediately after the vertical dimension:
+      if (present(dim4_bnds)) then
+         if (present(dim3_pos)) then
+            index = dim3_pos + 1
+         else
+            index = 3
+         end if
+         dim_bounds(index,1) = dim4_bnds(1)
+         dim_bounds(index,2) = dim4_bnds(2)
+      end if
+
       !
       ! Read netCDF file
       !
@@ -838,8 +866,7 @@ subroutine infld_real8_3d(varname, ncid, field, readvar, dim3name,         &
       ! If field is on file:
       !
       if (readvar_tmp) then
-         call print_input_field_info(dimlens, ndims, 2, 4, dim_bounds, 3,     &
-              varname, subname)
+         pdims = ndims
          ! Check to make sure that any 'extra' dimension is time
          if (ndims > target_ndims + 1) then
             call safe_endrun(subname//': too many dimensions for '//trim(varname))
@@ -855,9 +882,12 @@ subroutine infld_real8_3d(varname, ncid, field, readvar, dim3name,         &
                   call safe_endrun(subname//errormsg)
                end if
             end if
+            pdims = target_ndims
          else if (ndims < target_ndims) then
             call safe_endrun(subname//': too few dimensions for '//trim(varname))
          end if ! No else, things are okay
+         call print_input_field_info(dimlens, pdims, 2, 4, dim_bounds, 3,     &
+              varname, subname)
          !
          ! Get array dimension id's and sizes
          !
@@ -866,8 +896,10 @@ subroutine infld_real8_3d(varname, ncid, field, readvar, dim3name,         &
          arraydimsize(3) = (dim_bounds(3,2) - dim_bounds(3,1) + 1)
          do jndex = 1, 3
             if (arraydimsize(jndex) /= size(field, jndex)) then
-               write(errormsg, '(4a,i0)') ': Mismatch between array bounds ', &
-                    'and field size for ', trim(varname), ', dimension ', jndex
+               write(errormsg, '(4a,i0)') ': Mismatch between array bounds ',     &
+                    'and field size for ', trim(varname), ', dimension ', jndex,  &
+                    ', variable array bounds = ', arraydimsize(jndex),            &
+                    ', file field size = ', size(field, jndex)
                call safe_endrun(subname//errormsg)
             end if
          end do
@@ -947,8 +979,6 @@ subroutine infld_real8_3d(varname, ncid, field, readvar, dim3name,         &
 
       readvar = readvar_tmp
 
-      return
-
    end subroutine infld_real8_3d
 
 end module cam_field_read

test/unit/python/sample_files/write_init_files/phys_vars_init_check_4D.F90

@@ -33,19 +33,20 @@ module phys_vars_init_check_4D
    integer, public, parameter ::          PARAM = 2
    integer, public, parameter :: READ_FROM_FILE = 3
    ! Total number of physics-related variables:
-   integer, public, parameter :: phys_var_num = 2
+   integer, public, parameter :: phys_var_num = 3
    integer, public, parameter :: phys_const_num = 16
 
    !Max length of physics-related variable standard names:
    integer, public, parameter :: std_name_len = 25
 
    ! Max length of input (IC) file variable names:
-   integer, public, parameter :: ic_name_len = 5
+   integer, public, parameter :: ic_name_len = 8
 
    ! Physics-related input variable standard names:
    character(len=25), public, protected :: phys_var_stdnames(phys_var_num) = (/ &
       'potential_temperature    ', &
-      'air_pressure_at_sea_level' /)
+      'air_pressure_at_sea_level', &
+      'eddy_length_scale        ' /)
 
    character(len=36), public, protected :: phys_const_stdnames(phys_const_num) = (/ &
       "ccpp_constituent_minimum_values     ", &
@@ -65,17 +66,20 @@ module phys_vars_init_check_4D
       "suite_name                          ", &
       "suite_part                          " /)
    !Array storing all registered IC file input names for each variable:
-   character(len=5), public, protected :: input_var_names(1, phys_var_num) = reshape((/ &
-      'theta', &
-      'slp  ' /), (/1, phys_var_num/))
+   character(len=8), public, protected :: input_var_names(1, phys_var_num) = reshape((/ &
+      'theta   ', &
+      'slp     ', &
+      'eddy_len' /), (/1, phys_var_num/))
 
    ! Array indicating whether or not variable is protected:
    logical, public, protected :: protected_vars(phys_var_num)= (/ &
+      .false., &
       .false., &
       .false. /)
 
    ! Variable state (UNINITIALIZED, INTIIALIZED, PARAM or READ_FROM_FILE):
    integer, public, protected :: initialized_vars(phys_var_num)= (/ &
+      UNINITIALIZED, &
       UNINITIALIZED, &
       UNINITIALIZED /)
 

test/unit/python/sample_files/write_init_files/physics_inputs_4D.F90

@@ -40,7 +40,7 @@ subroutine physics_read_data(file, suite_names, timestep, read_initialized_varia
       use phys_vars_init_check_4D,   only: phys_var_num, phys_var_stdnames, input_var_names, std_name_len, is_initialized
       use ccpp_constituent_prop_mod, only: ccpp_constituent_prop_ptr_t
       use cam_logfile,               only: iulog
-      use physics_types_4D,          only: slp, theta
+      use physics_types_4D,          only: eddy_len, slp, theta
 
       ! Dummy arguments
       type(file_desc_t),          intent(inout) :: file
@@ -150,7 +150,10 @@ subroutine physics_read_data(file, suite_names, timestep, read_initialized_varia
                         call read_field(file, 'potential_temperature', input_var_names(:,name_idx), 'lev', timestep, theta)
 
                      case ('air_pressure_at_sea_level')
-                        call endrun('Cannot read slp from file'//', slp has unsupported dimension, timestep_for_physics.')
+                        call read_field(file, 'air_pressure_at_sea_level', input_var_names(:,name_idx), 'lev', timestep, slp)
+
+                     case ('eddy_length_scale')
+                        call endrun('Cannot read eddy_len from file'//', eddy_len has unsupported dimension, vertical_layer_dimension (dimension 4).')
 
                   end select !read variables
                end select !special indices

test/unit/python/sample_files/write_init_files/physics_inputs_bvd.F90

@@ -150,7 +150,7 @@ subroutine physics_read_data(file, suite_names, timestep, read_initialized_varia
                         call read_field(file, 'potential_temperature', input_var_names(:,name_idx), 'lev', timestep, theta)
 
                      case ('air_pressure_at_sea_level')
-                        call endrun('Cannot read slp from file'//', slp has unsupported dimension, band_number.')
+                        call endrun('Cannot read slp from file'//', slp has unsupported dimension, band_number (dimension 2).')
 
                   end select !read variables
                end select !special indices

test/unit/python/sample_files/write_init_files/temp_adjust_4D.F90

@@ -18,14 +18,15 @@ MODULE temp_adjust
 !! \htmlinclude arg_table_temp_adjust_run.html
 !!
   SUBROUTINE temp_adjust_run(nbox, lev, temp_layer,    &
-    slp, timestep, errmsg, errflg)
+    slp, eddy_len, timestep, errmsg, errflg)
 !----------------------------------------------------------------
    IMPLICIT NONE
 !----------------------------------------------------------------
 
    integer,            intent(in)    :: nbox, lev
    REAL(kind_phys),    intent(inout) :: temp_layer(:, :)
    real(kind_phys),    intent(in)    :: slp(:,:,:)
+   real(kind_phys),    intent(in)    :: eddy_len(:,:,:,:)
    real(kind_phys),    intent(in)    :: timestep
    character(len=512), intent(out)   :: errmsg
    integer,            intent(out)   :: errflg
@@ -45,6 +46,9 @@ SUBROUTINE temp_adjust_run(nbox, lev, temp_layer,    &
           !Add a made-up term which uses slp:
           temp_layer(box_index, lev_index) = temp_layer(box_index, lev_index) &
                + 0._kind_phys*slp(box_index, lev_index, 0)
+          !Add a made-up term which uses eddy_len:
+          temp_layer(box_index, lev_index) = temp_layer(box_index, lev_index) &
+               + 1._kind_phys*eddy_len(box_index, lev_index, 0, 0)
        end do
     end do
 

test/unit/python/sample_files/write_init_files/temp_adjust_4D.meta

@@ -30,6 +30,13 @@
   type = real
   kind = kind_phys
   intent = in
+[ eddy_len ]
+  standard_name = eddy_length_scale
+  units = m
+  dimensions = (horizontal_loop_extent, vertical_layer_dimension, timestep_for_physics, vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
 [ timestep ]
   standard_name = timestep_for_physics
   long_name = time step

test/unit/python/sample_files/write_init_files/var_4D_reg.xml

@@ -15,7 +15,7 @@
     </variable>
     <variable local_name="eddy_len" standard_name="eddy_length_scale"
               units="m" type="real" kind="kind_phys"  allocatable="pointer">
-      <dimensions>horizontal_dimension</dimensions>
+      <dimensions>horizontal_dimension vertical_layer_dimension timestep_for_physics vertical_layer_dimension</dimensions>
       <ic_file_input_names>eddy_len</ic_file_input_names>
     </variable>
   </file>