Feature/gefs soilt fix

sorc/chgres_cube.fd/input_data.F90

@@ -48,6 +48,8 @@ module input_data
 
  private
 
+ public :: check_soilt
+
 ! Fields associated with the atmospheric model.
 
  type(esmf_field), public              :: dzdt_input_grid       !< vert velocity
@@ -69,7 +71,8 @@ module input_data
  integer, public                 :: veg_type_landice_input = 15 !< NOAH land ice option
                                                                 !< defined at this veg type.
                                                                 !< Default is igbp.
-
+ integer, parameter              :: icet_default = 265.0    !< Default value of soil and skin
+                                                            !< temperature (K) over ice.
  type(esmf_field), public        :: canopy_mc_input_grid    !< canopy moist content
  type(esmf_field), public        :: f10m_input_grid         !< log((z0+10)*1/z0)
  type(esmf_field), public        :: ffmm_input_grid         !< log((z0+z1)*1/z0)
@@ -4597,15 +4600,15 @@ subroutine read_input_sfc_grib2_file(localpet)
 
    integer                               :: rc, varnum, iret, i, j,k
    integer                               :: ncid2d, varid, varsize
-   integer, parameter                    :: icet_default = 265.0
+   !integer, parameter                    :: icet_default = 265.0
 
    logical                               :: exist, rap_latlon
 
    real(esmf_kind_r4)                    :: value
 
-   real(esmf_kind_r4), allocatable       :: dummy2d(:,:),tsk_save(:,:),icec_save(:,:)
+   real(esmf_kind_r4), allocatable       :: dummy2d(:,:),icec_save(:,:)
    real(esmf_kind_r4), allocatable       :: dummy1d(:)
-   real(esmf_kind_r8), allocatable       :: dummy2d_8(:,:),dummy2d_82(:,:)
+   real(esmf_kind_r8), allocatable       :: dummy2d_8(:,:),dummy2d_82(:,:),tsk_save(:,:)
    real(esmf_kind_r8), allocatable       :: dummy3d(:,:,:), dummy3d_stype(:,:,:)
    integer(esmf_kind_i4), allocatable    :: slmsk_save(:,:)
    integer(esmf_kind_i8), allocatable    :: dummy2d_i(:,:)
@@ -4834,7 +4837,7 @@ subroutine read_input_sfc_grib2_file(localpet)
    print*,"- READ SKIN TEMPERATURE."
    rc = grb2_inq(the_file, inv_file, ':TMP:',':surface:', data2=dummy2d)
    if (rc <= 0 ) call error_handler("READING SKIN TEMPERATURE.", rc)
-   tsk_save(:,:) = dummy2d
+   tsk_save(:,:) = real(dummy2d,esmf_kind_r8)
    dummy2d_8 = real(dummy2d,esmf_kind_r8)
    do j = 1, j_input
      do i = 1, i_input
@@ -5390,14 +5393,7 @@ subroutine read_input_sfc_grib2_file(localpet)
    vname = "soilt"
    vname_file = ":TSOIL:"
    call read_grib_soil(the_file,inv_file,vname,vname_file,dummy3d,rc)
-   do k=1,lsoil_input
-     do j = 1, j_input
-       do i = 1, i_input
-         if (slmsk_save(i,j) == 0_esmf_kind_i4 ) dummy3d(i,j,k) = tsk_save(i,j)
-         if (slmsk_save(i,j) == 2_esmf_kind_i4 ) dummy3d(i,j,k) = icet_default
-       enddo
-     enddo
-   enddo
+   call check_soilt(dummy3d,slmsk_save,tsk_save)
    print*,'soilt ',maxval(dummy3d),minval(dummy3d)
 
    deallocate(tsk_save, slmsk_save)
@@ -6607,4 +6603,33 @@ recursive subroutine quicksort(a, first, last)
   if (j+1 < last)  call quicksort(a, j+1, last)
 end subroutine quicksort
 
+!> Check for and replace erroneous values in soil temperature.
+!!
+!! @param soilt    [inout] 3-dimensional soil temperature arrray
+!! @param landmask [in]    landmask of the input grid
+!! @param skint    [in]    2-dimensional skin temperature array
+!! @author Larissa Reames CIMMS/NSSL
+
+subroutine check_soilt(soilt, landmask, skint)
+  implicit none
+  real(esmf_kind_r8), intent(inout) ::  soilt(i_input,j_input,lsoil_input)
+  real(esmf_kind_r8), intent(in)    ::  skint(i_input,j_input)
+  integer(esmf_kind_i4), intent(in)    ::  landmask(i_input,j_input)
+
+  integer                           :: i, j, k
+
+  do k=1,lsoil_input
+    do j = 1, j_input
+      do i = 1, i_input
+        if (landmask(i,j) == 0_esmf_kind_i4 ) then 
+          soilt(i,j,k) = skint(i,j)
+        else if (landmask(i,j) == 1_esmf_kind_i4 .and. soilt(i,j,k) > 350.0_esmf_kind_r8) then 
+          soilt(i,j,k) = skint(i,j)
+        else if (landmask(i,j) == 2_esmf_kind_i4 ) then 
+          soilt(i,j,k) = icet_default
+        endif
+      enddo
+    enddo
+  enddo
+end subroutine check_soilt
  end module input_data

tests/chgres_cube/CMakeLists.txt

@@ -1,7 +1,7 @@
 # This is the cmake build file for the tests directory of the
 # UFS_UTILS project.
 #
-# George Gayno, Lin Gan, Ed Hartnett
+# George Gayno, Lin Gan, Ed Hartnett, Larissa Reames
 
 # Include cmake to allow parallel I/O tests.
 include (LibMPI)
@@ -58,3 +58,13 @@ add_mpi_test(ftst_search_util
   EXECUTABLE ${CMAKE_CURRENT_BINARY_DIR}/ftst_search_util
   NUMPROCS 1
   TIMEOUT 60)
+
+add_executable(ftst_sfc_input_data ftst_sfc_input_data.F90)
+target_link_libraries(ftst_sfc_input_data
+  chgres_cube_lib)
+
+# Cause test to be run with MPI.
+add_mpi_test(ftst_sfc_input_data
+  EXECUTABLE ${CMAKE_CURRENT_BINARY_DIR}/ftst_sfc_input_data
+  NUMPROCS 4
+  TIMEOUT 60)

tests/chgres_cube/ftst_sfc_input_data.F90

@@ -0,0 +1,94 @@
+ program test_sfc_input_data
+
+! Test sfc_input_data. For now only one test is performed
+! on the subroutine soilt_check using a sample 3x3 matrix.
+!
+! author: Larissa Reames (larissa.reames@noaa.gov)
+
+ use mpi
+ use esmf
+ use input_data
+ use model_grid
+
+ implicit none
+
+ integer :: ierr
+
+ integer(esmf_kind_i4), allocatable :: mask(:,:)
+ real(esmf_kind_r8), allocatable    :: soilt_bad(:,:,:), &
+                                       soilt_updated(:,:,:), &
+                                       soilt_correct(:,:,:)
+ real(esmf_kind_r8), allocatable    :: skint(:,:)
+
+ call mpi_init(ierr)
+
+ lsoil_input = 2
+ i_input = 3
+ j_input = 3
+
+ allocate(mask(i_input,j_input))
+ allocate(skint(i_input,j_input))
+ allocate(soilt_bad(i_input,j_input,lsoil_input))
+ allocate(soilt_updated(i_input,j_input,lsoil_input))
+ allocate(soilt_correct(i_input,j_input,lsoil_input))
+
+!--------------------------------------------------------
+! These variables are used to test all three functions of
+! the check_soilt routine (i.e., replace water point soil 
+! temperature and excessive land point soil temperature 
+! with skin temperature and store a default ice column
+! temperature in the soil temperature array since this
+! array isn't available in grib2 files)
+!--------------------------------------------------------
+
+! Definition of the mask. The '1' indicates an isolated
+! island or lake depending on the field type and '2' 
+! indicates sea ice.
+
+ mask = reshape((/0, 1, 0,  0, 0, 1, 0, 0, 2/),(/3,3/))
+
+! Soil temperature array with some values that are all incorrect
+! based on landmask type. This will be the input soil array
+
+ soilt_bad(1:i_input,1:j_input,1) = reshape((/0., 9999999.9, 0., &
+                                              0., 0., 295.0, &
+                                              0., 25.0, 25.0 /),(/3,3/))
+ soilt_bad(1:i_input,1:j_input,2) = reshape((/0., 9999999.9, 0., &
+                                              0., 0., 295.0, &
+                                              0., 25.0, 25.0 /),(/3,3/))
+
+! Subjective, reasonable skin temperature array.
+
+ skint = reshape((/285.0, 295.0, 280.0, 281.0, 282.0, 283.0, 284.0, 285.0, 260.0/), &
+                 (/3,3/))
+
+!--------------------------------------------------------
+! This is the corrected soil temperature array that 
+! should be passed back from the check_soilt routine
+!--------------------------------------------------------
+
+ soilt_correct(1:i_input,1:j_input,1) = reshape( (/285.0, 295.0, 280.0, &
+                                                   281.0, 282.0, 295.0, &
+                                                   284.0, 285.0, 265.0/), (/3,3/))
+ soilt_correct(1:i_input,1:j_input,2) = reshape( (/285.0, 295.0, 280.0, &
+                                                   281.0, 282.0, 295.0, &
+                                                   284.0, 285.0, 265.0/), (/3,3/))
+
+
+ print*,"Starting test of check_soilt subroutine."
+
+ soilt_updated = soilt_bad
+ call check_soilt(soilt_updated,mask,skint)
+
+ if (any(soilt_updated /= soilt_correct)) then
+   print*,'TEST FAILED '
+   print*,'ARRAY SHOULD BE:', soilt_correct
+   print*,'ARRAY FROM TEST:', soilt_updated
+   stop 2
+ endif
+
+ call mpi_finalize(ierr)
+
+ print*,"SUCCESS!"
+
+ end program test_sfc_input_data