Various chgres_cube bug fixes and code improvements

Updated logic for reading soil type from the geogrid file.
If the user provides a geogrid file and any of the calls in this 
process fail, the code will now fail with an error. Also, added a 
check to ensure input and geogrid grid sizes match.

Add a sanity check to fail if the variable mapping file is empty 
or improperly formatted such that the code finds no variable entries to read.

Bug fix when using HRRR or RAP data with climo surface fields that 
caused erroneous soil types.

Clean up extraneous print statements.

Addresses #228

sorc/chgres_cube.fd/atmosphere.F90

@@ -1163,9 +1163,6 @@ subroutine read_vcoord_info
  endif
 
  print*
- do k = 1, levp1_target
-    print*,'VCOORD FOR LEV ', k, 'IS: ', vcoord_target(k,:)
- enddo
 
  close(14)
 

sorc/chgres_cube.fd/input_data.F90

@@ -4603,7 +4603,7 @@ subroutine read_input_sfc_grib2_file(localpet)
    character(len=50)                      :: method
 
    integer                               :: rc, varnum, iret, i, j,k
-   integer                               :: ncid2d, varid
+   integer                               :: ncid2d, varid, varsize
    integer, parameter                    :: icet_default = 265.0
 
    logical                               :: exist, rap_latlon
@@ -4875,40 +4875,41 @@ subroutine read_input_sfc_grib2_file(localpet)
    vname=":SOTYP:"                                     
    rc = grb2_inq(the_file, inv_file, vname,slev, data2=dummy2d)
    !failed => rc = 0
-   if (rc <= 0 .and. (trim(to_upper(external_model))=="HRRR" .or. rap_latlon)) then
+   if (rc <= 0 .and. (trim(to_upper(external_model))=="HRRR" .or. rap_latlon) .and. geo_file .ne. "NULL")  then
      ! Some HRRR and RAP files don't have dominant soil type in the output, but the geogrid files
      ! do, so this gives users the option to provide the geogrid file and use input soil
      ! type 
      print*, "OPEN GEOGRID FILE ", trim(geo_file)
      rc = nf90_open(geo_file,NF90_NOWRITE,ncid2d)
-     ! failed => rc < 0
-     if (rc == 0) then
-       print*, "INQUIRE ABOUT SOIL TYPE FROM GEOGRID FILE"
-       rc = nf90_inq_varid(ncid2d,"SCT_DOM",varid)
-       ! failed => rc < 0
-       if (rc<0) print*, "ERROR FINDING SCT_DOM IN GEOGRID FILE"
-       if (rc == 0) then
-         print*, "READ SOIL TYPE FROM GEOGRID FILE "
-         rc = nf90_get_var(ncid2d,varid,dummy2d)
-         ! failed => rc < 0
-         if (rc<0) print*, "ERROR READING SCT_DOM FROM FILE"
-         print*, "min max dummy2d = ", minval(dummy2d), maxval(dummy2d)
-       endif
-       print*, "INQUIRE ABOUT SOIL TYPE FRACTIONS FROM GEOGRID FILE"
-       rc = nf90_inq_varid(ncid2d,"SOILCTOP",varid)
-       ! failed => rc < 0
-       if (rc<0) print*, "ERROR FINDING SOILCTOP IN GEOGRID FILE"
-       if (rc == 0) then
-         print*, "READ SOIL TYPE FRACTIONS FROM GEOGRID FILE "
-         rc = nf90_get_var(ncid2d,varid,dummy3d_stype)
-         ! failed => rc < 0
-         if (rc<0) print*, "ERROR READING SCT_DOM FROM FILE"
-         print*, "min max dummy3d_stype = ", minval(dummy3d_stype), maxval(dummy3d_stype)
-       endif
+     call netcdf_err(rc,"READING GEOGRID FILE")
 
-       print*, "CLOSE GEOGRID FILE "
-       iret = nf90_close(ncid2d)
-     endif
+     print*, "INQURE ABOUT DIM IDS"
+     rc = nf90_inq_dimid(ncid2d,"west_east",varid)
+     call netcdf_err(rc,"READING west_east DIMENSION FROM GEOGRID FILE")
+
+     rc = nf90_inquire_dimension(ncid2d,varid,len=varsize)
+     call netcdf_err(rc,"READING west_east DIMENSION SIZE")
+     if (varsize .ne. i_input) call error_handler ("GEOGRID FILE GRID SIZE DIFFERS FROM INPUT DATA.", -1)
+
+     print*, "INQUIRE ABOUT SOIL TYPE FROM GEOGRID FILE"
+     rc = nf90_inq_varid(ncid2d,"SCT_DOM",varid)
+     call netcdf_err(rc,"FINDING SCT_DOM IN GEOGRID FILE")
+
+     print*, "READ SOIL TYPE FROM GEOGRID FILE "
+     rc = nf90_get_var(ncid2d,varid,dummy2d)
+     call netcdf_err(rc,"READING SCT_DOM FROM FILE")
+
+     print*, "INQUIRE ABOUT SOIL TYPE FRACTIONS FROM GEOGRID FILE"
+     rc = nf90_inq_varid(ncid2d,"SOILCTOP",varid)
+     call netcdf_err(rc,"FINDING SOILCTOP IN GEOGRID FILE")
+
+     print*, "READ SOIL TYPE FRACTIONS FROM GEOGRID FILE "
+     rc = nf90_get_var(ncid2d,varid,dummy3d_stype)
+     call netcdf_err(rc,"READING SCT_DOM FROM FILE")
+
+     print*, "CLOSE GEOGRID FILE "
+     iret = nf90_close(ncid2d)
+
 
      ! There's an issue with the geogrid file containing soil type water at land points. 
      ! This correction replaces the soil type at these points with the soil type with
@@ -4956,6 +4957,8 @@ subroutine read_input_sfc_grib2_file(localpet)
      where(slmsk_save == 1) dummy2d_i = 1
 
      call search(dummy2d_8,dummy2d_i,i_input,j_input,1,230)
+   else
+      dummy2d_8=real(dummy2d,esmf_kind_r8)
    endif
 
    print*,'sotype ',maxval(dummy2d_8),minval(dummy2d_8)

sorc/chgres_cube.fd/model_grid.F90

@@ -373,8 +373,6 @@ subroutine define_input_grid_gaussian(localpet, npets)
  if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
     call error_handler("IN GridGetCoord", rc)
 
- print*,'bounds for corners ',localpet,clb(1),cub(1),clb(2),cub(2)
-
  do j = clb(2), cub(2)
    do i = clb(1), cub(1)
      lon_corner_src_ptr(i,j) = longitude(i,1) - (0.5_esmf_kind_r8*deltalon)
@@ -699,7 +697,6 @@ subroutine define_input_grid_gfs_grib2(localpet, npets)
  deallocate(lat4, lon4)
 
  deltalon = abs(longitude(2,1)-longitude(1,1))
- if(localpet==0) print*, "deltalon = ", deltalon
 
  print*,"- CALL FieldScatter FOR INPUT GRID LONGITUDE."
  call ESMF_FieldScatter(longitude_input_grid, longitude, rootpet=0, rc=rc)
@@ -745,9 +742,6 @@ subroutine define_input_grid_gfs_grib2(localpet, npets)
    enddo
  enddo
 
- if(localpet==0) print*, "lon first = ", lon_src_ptr(1:10,1)
- if(localpet==0) print*, "lat first = ", lat_src_ptr(1,1:10)
-
  print*,"- CALL GridAddCoord FOR INPUT GRID."
  call ESMF_GridAddCoord(input_grid, &
                         staggerloc=ESMF_STAGGERLOC_CORNER, rc=rc)
@@ -846,7 +840,6 @@ subroutine define_input_grid_grib2(localpet, npets)
 
   ! Wgrib2 can't properly read the lat/lon arrays of data on NCEP rotated lat/lon
   ! grids, so read in lat/lon from fixed coordinate file
-  print*, 'temp num =', temp_num
   if (trim(temp_num)=="3.32769" .or. trim(temp_num)=="32769") then
 
      input_grid_type = "rotated_latlon"
@@ -1011,7 +1004,6 @@ subroutine define_input_grid_grib2(localpet, npets)
    if(ESMF_logFoundError(rcToCheck=error,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
       call error_handler("IN GridGetCoord", error)
 
-   print*,'bounds for corners ',localpet,clb(1),cub(1),clb(2),cub(2)
 
   ! If we have data on a lat/lon or lambert grid, create staggered coordinates
   if(trim(input_grid_type)=="latlon" .or. trim(input_grid_type) == "lambert") then
@@ -1048,7 +1040,6 @@ subroutine define_input_grid_grib2(localpet, npets)
          i = index(temp_msg, "Dx ") + len("Dx ")
          j = index(temp_msg," m Dy")
          read(temp_msg(i:j-1),"(F9.6)") dx
-         print*, "DX = ", dx
        endif
        call MPI_BARRIER(MPI_COMM_WORLD,error)
        call MPI_BCAST(dx,1,MPI_REAL8,0,MPI_COMM_WORLD,error)
@@ -1071,8 +1062,6 @@ subroutine define_input_grid_grib2(localpet, npets)
     call netcdf_err(error, 'reading field id' )
     error=nf90_get_var(ncid, id_var, lat_corners)
     call netcdf_err(error, 'reading field' )
-    print*, 'min, max lat corners =', minval(lat_corners), maxval(lat_corners)
-    print*, 'min, max lon corners =', minval(lon_corners), maxval(lon_corners)
 
     do j = clb(2),cub(2)
       do i = clb(1), cub(1)

sorc/chgres_cube.fd/program_setup.f90

@@ -444,7 +444,7 @@ subroutine read_varmap
      if ( trim(line) .eq. '' ) cycle
      nvars = nvars+1
    enddo
-
+   if ( nvars == 0) call error_handler("VARMAP FILE IS EMPTY.", -1)
 
    allocate(chgres_var_names(nvars))
    allocate(field_var_names(nvars))