Merge pull request ESCOMP#489 from peverwhee/f2000climo_ncold

cam6_3_043: Read ncdata file with either ncol or ncol_d dimension (SE dycore)

doc/ChangeLog

@@ -1,5 +1,61 @@
 ===============================================================
 
+Tag name: cam6_3_043
+Originator(s): peverwhee
+Date: 11 Jan 2022
+One-line Summary: Enable use of initial data file with either horizontal
+dimension name "ncol" or "ncol_d"
+
+Github PR URL: https://github.com/ESCOMP/CAM/pull/489
+
+Purpose of changes (include the issue number and title text for each relevant GitHub issue):
+
+Fixes issue #339 -- SE dycore not reading default ncdata for ne30np4 F2000
+
+Describe any changes made to build system:
+
+Describe any changes made to the namelist:
+
+List any changes to the defaults for the boundary datasets:
+
+Describe any substantial timing or memory changes:
+
+Code reviewed by: gold2718, cacraigucar, fvitt, nusbaume
+
+List all files eliminated:
+
+List all files added and what they do:
+
+List all existing files that have been modified, and describe the changes:
+
+M	src/dynamics/se/dyn_comp.F90
+ - use grid dimension name variable from dyn_grid
+
+M	src/dynamics/se/dyn_grid.F90
+ - make grid dimension name variable publicly available
+ - add logic to add '_d' to attributes lat, lon, and area if the dimension is 'ncol_d'
+ - make get_hdim_name dummy variables generic
+
+If there were any failures reported from running test_driver.sh on any test
+platform, and checkin with these failures has been OK'd by the gatekeeper,
+then copy the lines from the td.*.status files for the failed tests to the
+appropriate machine below.  All failed tests must be justified.
+
+cheyenne/intel/aux_cam: ALL PASS
+
+izumi/nag/aux_cam: 
+ DAE_Vmct.f45_f45_mg37.FHS94.izumi_nag.cam-dae (Overall:FAIL) details:
+  - pre-existing failure
+
+izumi/pgi/aux_cam: ALL PASS
+
+CAM tag used for the baseline comparison tests if different than previous
+tag:
+
+Summarize any changes to answers: None, all BFB
+
+===============================================================
+
 Tag name: cam6_3_042
 Originator(s): fvitt, tilmes, emmons
 Date: 10 Jan 2022

src/dynamics/se/dyn_comp.F90

@@ -11,7 +11,8 @@ module dyn_comp
 use cam_control_mod,        only: initial_run
 use cam_initfiles,          only: initial_file_get_id, topo_file_get_id, pertlim
 use phys_control,           only: use_gw_front, use_gw_front_igw, waccmx_is
-use dyn_grid,               only: ini_grid_name, timelevel, hvcoord, edgebuf
+use dyn_grid,               only: ini_grid_name, timelevel, hvcoord, edgebuf, &
+                                  ini_grid_hdim_name
 
 use cam_grid_support,       only: cam_grid_id, cam_grid_get_gcid, &
                                   cam_grid_dimensions, cam_grid_get_dim_names, &
@@ -1201,7 +1202,7 @@ subroutine read_inidat(dyn_in)
    integer                          :: kptr, m_cnst
    type(EdgeBuffer_t)               :: edge
 
-   character(len=max_fieldname_len) :: dimname, varname
+   character(len=max_fieldname_len) :: varname
    integer                          :: ierr
 
    integer                          :: rndm_seed_sz
@@ -1351,18 +1352,18 @@ subroutine read_inidat(dyn_in)
       allocate(dbuf3(npsq,nlev,nelemd))
 
       ! Check that columns in IC file match grid definition.
-      call check_file_layout(fh_ini, elem, dyn_cols, 'ncdata', .true., dimname)
+      call check_file_layout(fh_ini, elem, dyn_cols, 'ncdata', .true.)
 
       ! Read 2-D field
 
       fieldname  = 'PS'
       fieldname2 = 'PSDRY'
       if (dyn_field_exists(fh_ini, trim(fieldname), required=.false.)) then
          inic_wet = .true.
-         call read_dyn_var(trim(fieldname), fh_ini, dimname, dbuf2)
+         call read_dyn_var(trim(fieldname), fh_ini, ini_grid_hdim_name, dbuf2)
       elseif (dyn_field_exists(fh_ini, trim(fieldname2), required=.false.)) then
          inic_wet = .false.
-         call read_dyn_var(trim(fieldname2), fh_ini, dimname, dbuf2)
+         call read_dyn_var(trim(fieldname2), fh_ini, ini_grid_hdim_name, dbuf2)
       else
          call endrun(trim(sub)//': PS or PSDRY must be on GLL grid')
       end if
@@ -1386,7 +1387,7 @@ subroutine read_inidat(dyn_in)
       ! Read in 3-D fields
 
       if (dyn_field_exists(fh_ini, 'U')) then
-         call read_dyn_var('U', fh_ini, dimname, dbuf3)
+         call read_dyn_var('U', fh_ini, ini_grid_hdim_name, dbuf3)
       else
          call endrun(trim(sub)//': U not found')
       end if
@@ -1402,7 +1403,7 @@ subroutine read_inidat(dyn_in)
       end do
 
       if (dyn_field_exists(fh_ini, 'V')) then
-         call read_dyn_var('V', fh_ini, dimname, dbuf3)
+         call read_dyn_var('V', fh_ini, ini_grid_hdim_name, dbuf3)
       else
          call endrun(trim(sub)//': V not found')
       end if
@@ -1417,7 +1418,7 @@ subroutine read_inidat(dyn_in)
       end do
 
       if (dyn_field_exists(fh_ini, 'T')) then
-         call read_dyn_var('T', fh_ini, dimname, dbuf3)
+         call read_dyn_var('T', fh_ini, ini_grid_hdim_name, dbuf3)
       else
          call endrun(trim(sub)//': T not found')
       end if
@@ -1494,7 +1495,7 @@ subroutine read_inidat(dyn_in)
    do m_cnst = 1, pcnst
       if (cnst_read_iv(m_cnst) .and. .not. cnst_is_a_water_species(cnst_name(m_cnst))) then
          if (dyn_field_exists(fh_ini, trim(cnst_name(m_cnst)), required=.false.)) then
-            call check_file_layout(fh_ini, elem, dyn_cols, 'ncdata', .true., dimname)
+            call check_file_layout(fh_ini, elem, dyn_cols, 'ncdata', .true.)
             exit
          end if
       end if
@@ -1512,7 +1513,7 @@ subroutine read_inidat(dyn_in)
       end if
 
       if (found) then
-         call read_dyn_var(trim(cnst_name(m_cnst)), fh_ini, dimname, dbuf3)
+         call read_dyn_var(trim(cnst_name(m_cnst)), fh_ini, ini_grid_hdim_name, dbuf3)
       else
          call cnst_init_default(m_cnst, latvals, lonvals, dbuf3, pmask)
       end if
@@ -2015,7 +2016,7 @@ end subroutine set_phis
 
 !========================================================================================
 
-subroutine check_file_layout(file, elem, dyn_cols, file_desc, dyn_ok, dimname)
+subroutine check_file_layout(file, elem, dyn_cols, file_desc, dyn_ok)
 
    ! This routine is only called when data will be read from the initial file.  It is not
    ! called when the initial file is only supplying vertical coordinate info.
@@ -2025,7 +2026,6 @@ subroutine check_file_layout(file, elem, dyn_cols, file_desc, dyn_ok, dimname)
    integer,           intent(in)    :: dyn_cols
    character(len=*),  intent(in)    :: file_desc
    logical,           intent(in)    :: dyn_ok ! .true. iff ncol_d is okay
-   character(len=*),  intent(out)   :: dimname
 
    integer                          :: ncol_did, ncol_size
    integer                          :: ierr
@@ -2034,18 +2034,20 @@ subroutine check_file_layout(file, elem, dyn_cols, file_desc, dyn_ok, dimname)
    integer                          :: indx
    real(r8)                         :: dbuf2(npsq, nelemd)
    logical                          :: found
-   character(len=max_fieldname_len) :: dimname2, coordname
+   character(len=max_fieldname_len) :: coordname
 
    character(len=*), parameter      :: sub = 'check_file_layout'
    !----------------------------------------------------------------------------
 
    ! Check that number of columns in IC file matches grid definition.
+   if (trim(ini_grid_hdim_name) == 'none') then
+      call endrun(sub//': ERROR: no horizontal dimension in initial data file. &
+         &Cannot read data from file')
+   end if
 
-   call cam_grid_get_dim_names(cam_grid_id(ini_grid_name), dimname, dimname2)
-
-   ierr = pio_inq_dimid(file, trim(dimname), ncol_did)
+   ierr = pio_inq_dimid(file, trim(ini_grid_hdim_name), ncol_did)
    if (ierr /= PIO_NOERR) then
-      call endrun(sub//': ERROR: either ncol or ncol_d dimension not found in ' &
+      call endrun(sub//': ERROR: '//trim(ini_grid_hdim_name)//' dimension not found in ' &
          //trim(file_desc)//' file')
    end if
 
@@ -2055,18 +2057,18 @@ subroutine check_file_layout(file, elem, dyn_cols, file_desc, dyn_ok, dimname)
          write(iulog, '(a,2(a,i0))') trim(sub), ': ncol_size=', ncol_size, &
              ' : dyn_cols=', dyn_cols
       end if
-      call endrun(sub//': ERROR: dimension ncol size not same as in ncdata file')
+      call endrun(sub//': ERROR: dimension '//trim(ini_grid_hdim_name)//' size not same as in ncdata file')
    end if
 
-   ! Set coordinate name associated with dimname.
-   if (dimname == 'ncol') then
+   ! Set coordinate name associated with ini_grid_hdim_name.
+   if (trim(ini_grid_hdim_name) == 'ncol') then
       coordname = 'lat'
    else
       coordname = 'lat_d'
    end if
 
    !! Check to make sure file is in correct order
-   call read_dyn_var(coordname, file, dimname, dbuf2)
+   call read_dyn_var(coordname, file, ini_grid_hdim_name, dbuf2)
    found = .true.
    do ie = 1, nelemd
       indx = 1
@@ -2092,13 +2094,13 @@ subroutine check_file_layout(file, elem, dyn_cols, file_desc, dyn_ok, dimname)
       call endrun("ncdata file latitudes not in correct column order")
    end if
 
-   if (dimname == 'ncol') then
+   if (trim(ini_grid_hdim_name) == 'ncol') then
       coordname = 'lon'
    else
       coordname = 'lon_d'
    end if
 
-   call read_dyn_var(coordname, file, dimname, dbuf2)
+   call read_dyn_var(coordname, file, ini_grid_hdim_name, dbuf2)
    do ie = 1, nelemd
       indx = 1
       do j = 1, np

src/dynamics/se/dyn_grid.F90

@@ -62,6 +62,9 @@ module dyn_grid
 
 character(len=3), protected :: ini_grid_name
 
+! Name of horizontal grid dimension in initial file.
+character(len=6), protected :: ini_grid_hdim_name = ''
+
 integer, parameter :: ptimelevels = 2
 
 type (TimeLevel_t)         :: TimeLevel     ! main time level struct (used by tracers)
@@ -71,6 +74,7 @@ module dyn_grid
 
 public :: dyn_decomp
 public :: ini_grid_name
+public :: ini_grid_hdim_name
 public :: ptimelevels
 public :: TimeLevel
 public :: hvcoord
@@ -99,9 +103,6 @@ module dyn_grid
 character(len=shr_kind_cl), public :: se_grid_filename     = ''
 logical,                    public :: se_write_gll_corners = .false.
 
-! Name of horizontal grid dimension in initial file.
-character(len=6) :: ini_grid_hdim_name = ' '
-
 type(physics_column_t), allocatable, target :: local_dyn_columns(:)
 
 ! number of global dynamics columns. Set by SE dycore init.
@@ -548,7 +549,11 @@ subroutine physgrid_copy_attributes_d(gridname, grid_attribute_names)
       allocate(grid_attribute_names(3))
       ! For standard CAM-SE, we need to copy the area attribute.
       ! For physgrid, the physics grid will create area (GLL has area_d)
-      grid_attribute_names(1) = 'area'
+      if (trim(ini_grid_hdim_name) == 'ncol_d') then
+         grid_attribute_names(1) = 'area_d'
+      else
+         grid_attribute_names(1) = 'area'
+      end if
       grid_attribute_names(2) = 'np'
       grid_attribute_names(3) = 'ne'
    end if
@@ -661,19 +666,19 @@ end subroutine dyn_grid_get_elem_coords
 ! Private routines.
 !=========================================================================================
 
-subroutine get_hdim_name(fh_ini, ini_grid_hdim_name)
+subroutine get_hdim_name(fh_ptr, grid_hdim_name)
    use pio, only: pio_inq_dimid, pio_seterrorhandling
    use pio, only: PIO_BCAST_ERROR, PIO_NOERR
 
-   ! Determine whether the initial file uses 'ncol' or 'ncol_d' horizontal
+   ! Determine whether the supplied file uses 'ncol' or 'ncol_d' horizontal
    ! dimension in the unstructured grid.  It is also possible when using
-   ! analytic initial conditions that the initial file only contains
+   ! analytic initial conditions that the file only contains
    ! vertical coordinates.
    ! Return 'none' if that is the case.
 
    ! Arguments
-   type(file_desc_t),   pointer  :: fh_ini
-   character(len=6), intent(out) :: ini_grid_hdim_name ! horizontal dimension name
+   type(file_desc_t),   pointer  :: fh_ptr
+   character(len=6), intent(out) :: grid_hdim_name ! horizontal dimension name
 
    ! local variables
    integer  :: ierr, pio_errtype
@@ -683,33 +688,33 @@ subroutine get_hdim_name(fh_ini, ini_grid_hdim_name)
    !----------------------------------------------------------------------------
 
    ! Set PIO to return error flags.
-   call pio_seterrorhandling(fh_ini, PIO_BCAST_ERROR, pio_errtype)
+   call pio_seterrorhandling(fh_ptr, PIO_BCAST_ERROR, pio_errtype)
 
-   ! Check for ncol_d first just in case the initial file also contains fields on
+   ! Check for ncol_d first just in case the file also contains fields on
    ! the physics grid.
-   ierr = pio_inq_dimid(fh_ini, 'ncol_d', ncol_did)
+   ierr = pio_inq_dimid(fh_ptr, 'ncol_d', ncol_did)
    if (ierr == PIO_NOERR) then
 
-      ini_grid_hdim_name = 'ncol_d'
+      grid_hdim_name = 'ncol_d'
 
    else
 
       ! if 'ncol_d' not in file, check for 'ncol'
-      ierr = pio_inq_dimid(fh_ini, 'ncol', ncol_did)
+      ierr = pio_inq_dimid(fh_ptr, 'ncol', ncol_did)
 
       if (ierr == PIO_NOERR) then
 
-         ini_grid_hdim_name = 'ncol'
+         grid_hdim_name = 'ncol'
 
       else
 
-         ini_grid_hdim_name = 'none'
+         grid_hdim_name = 'none'
 
       end if
    end if
 
    ! Return PIO to previous error handling.
-   call pio_seterrorhandling(fh_ini, pio_errtype)
+   call pio_seterrorhandling(fh_ptr, pio_errtype)
 
 end subroutine get_hdim_name
 
@@ -800,7 +805,7 @@ subroutine define_cam_grids()
    ! '_d' suffixes and the physics grid will use the unadorned names.
    ! This allows fields on both the GLL and physics grids to be written to history
    ! output files.
-   if (fv_nphys > 0) then
+   if (trim(ini_grid_hdim_name) == 'ncol_d') then
       latname  = 'lat_d'
       lonname  = 'lon_d'
       ncolname = 'ncol_d'
@@ -840,7 +845,7 @@ subroutine define_cam_grids()
    ! With CSLAM if the initial file uses the horizontal dimension 'ncol' rather than
    ! 'ncol_d' then we need a grid object with the names ncol,lat,lon to read it.
    ! Create that grid object here if it's needed.
-   if (fv_nphys > 0 .and. ini_grid_hdim_name == 'ncol') then
+   if (fv_nphys > 0 .and. trim(ini_grid_hdim_name) == 'ncol') then
 
       lat_coord => horiz_coord_create('lat', 'ncol', ngcols_d,  &
          'latitude', 'degrees_north', 1, size(pelat_deg), pelat_deg, map=pemap)