Add QC checks for ICE source files

sorc/ocnice_prep.fd/docs/ocnice_prep.md

@@ -117,6 +117,12 @@ from their orientation along model dimensions (IJ) to eastward-northward (EN) di
 4. Vector fields are re-rotated back to model index direction (EN->IJ) and remapped back to their native stagger locations.
 5. Fields are written into a warm-start file.
 
+For CICE6, an additional 2-step QC is performed on the source restart files. This procedure accounts for possible inconsistent fields
+that may have been generated as part of the marine DA process. The QC process affects 4 fields only: ``aicen``, ``vicen``,
+``vsnon`` and ``Tsfcn``. The first step of the QC is to remove any values of the 4 fields which are located where the CICE grid has
+land. The second step is to remove any ice or snow volume where there is no associated ice coverage. The QC process does not impact
+source files which do not contain these physical inconsistencies.
+
 ## Note on rotation angles
 
 The rotation angle retrieved from the master grid file is ``anglet``, which is defined in the same sense as used MOM6. This

sorc/ocnice_prep.fd/ocniceprep.F90

@@ -26,6 +26,7 @@ program ocniceprep
   use arrays_mod ,     only : nbilin2d, nbilin3d, nconsd2d, bilin2d, bilin3d, consd2d
   use arrays_mod ,     only : mask3d, hmin, maskspval, eta
   use utils_mod  ,     only : getfield, packarrays, remap, dumpnc, nf90_err
+  use utils_mod  ,     only : zero_out_land_ice, zero_out_phantom_ice
   use utils_esmf_mod , only : createRH, remapRH, ChkErr, rotremap
   use restarts_mod ,   only : setup_icerestart, setup_ocnrestart
   use ocncalc_mod ,    only : calc_eta, vfill
@@ -47,13 +48,16 @@ program ocniceprep
   real(kind=8), allocatable, dimension(:,:)   :: out2d !< 2D destination grid output array
   real(kind=8), allocatable, dimension(:,:,:) :: out3d !< 3D destination grid output array
 
+  ! ice mask array for QC of ice files
+  integer(kind=4), allocatable, dimension(:,:) :: a2d  !< 2D land mask for ice
+  integer(kind=4), allocatable, dimension(:)   :: kmt  !< 1D land mask for ice
+  integer :: idx1, idx2, idx3, idx4
+
   character(len=120) :: errmsg
   character(len=120) :: meshfsrc, meshfdst
-  integer            :: nvalid
+  integer            :: nvalid, icnt
   integer            :: k,n,nn,rc,ncid,varid
   character(len=20)  :: vname
-  ! debug
-  integer :: i,j
 
   character(len=*), parameter :: u_FILE_u = __FILE__
 
@@ -121,6 +125,26 @@ program ocniceprep
   call getfield(trim(gridfile),  'depth', dims=(/nxr,nyr/), field=bathydst)
   call nf90_err(nf90_close(ncid), 'close: '//trim(gridfile))
 
+  if (.not. do_ocnprep) then
+  ! -----------------------------------------------------------------------------
+  ! obtain the land mask for the ice model on the source grid to QC the ice fields
+  ! -----------------------------------------------------------------------------
+
+     allocate(a2d(nxt,nyt)); a2d = 0
+     allocate(kmt(nxt*nyt)); kmt = 0
+
+     ! obtain the land mask for the source grid
+     vname = 'wet'
+     gridfile = trim(griddir)//fsrc(3:5)//'/'//'tripole.'//trim(fsrc)//'.nc'
+     call nf90_err(nf90_open(trim(gridfile), nf90_nowrite, ncid), &
+          'open: '//trim(gridfile))
+     call nf90_err(nf90_inq_varid(ncid, trim(vname), varid), 'get variable ID: '//trim(vname))
+     call nf90_err(nf90_get_var(ncid, varid, a2d), 'get variable: '//trim(vname))
+     call nf90_err(nf90_close(ncid), 'close: '//trim(gridfile))
+     kmt(:) = reshape(a2d, (/nxt*nyt/))
+     deallocate(a2d)
+  end if
+
   ! -----------------------------------------------------------------------------
   ! get the 3rd (vertical or ncat) dimension and variable attributes for the
   ! ocean file
@@ -248,6 +272,41 @@ program ocniceprep
   if (allocated(bilin3d))then
      call packarrays(trim(input_file), trim(wgtsdir)//fsrc(3:5)//'/',           &
           cos(angsrc), sin(angsrc), b3d, dims=(/nxt,nyt,nlevs/), nflds=nbilin3d, fields=bilin3d)
+
+     if (.not. do_ocnprep) then
+        ! -----------------------------------------------------------------------------
+        ! QC the source ice files after packing. Not every possible inconsistency is
+        ! checked but these are those known to create issues in the coupled model
+        ! -----------------------------------------------------------------------------
+
+        do n = 1,nbilin3d
+           if (trim(b3d(n)%var_name) == 'aicen')idx1 = n
+           if (trim(b3d(n)%var_name) == 'vicen')idx2 = n
+           if (trim(b3d(n)%var_name) == 'vsnon')idx3 = n
+           if (trim(b3d(n)%var_name) == 'Tsfcn')idx4 = n
+        end do
+
+        ! remove land values, if they exist
+        call zero_out_land_ice(kmt, bilin3d(idx1,:,:), icnt)
+        write(logunit, '(a,i8,a)')'removed ',icnt,' locations of '//trim(b3d(idx1)%var_name)//' from land'
+        call zero_out_land_ice(kmt, bilin3d(idx2,:,:), icnt)
+        write(logunit, '(a,i8,a)')'removed ',icnt,' locations of '//trim(b3d(idx2)%var_name)//' from land'
+        call zero_out_land_ice(kmt, bilin3d(idx3,:,:), icnt)
+        write(logunit, '(a,i8,a)')'removed ',icnt,' locations of '//trim(b3d(idx3)%var_name)//' from land'
+        call zero_out_land_ice(kmt, bilin3d(idx4,:,:), icnt)
+        write(logunit, '(a,i8,a)')'removed ',icnt,' locations of '//trim(b3d(idx4)%var_name)//' from land'
+
+        ! remove phantom-ice (aice=0, vice or vsno /= 0)
+        call zero_out_phantom_ice(bilin3d(idx1,:,:), bilin3d(idx2,:,:), icnt)
+        write(logunit, '(a,i8,a)')'removed ',icnt,' locations of phantom '//trim(b3d(idx2)%var_name)
+        call zero_out_phantom_ice(bilin3d(idx1,:,:), bilin3d(idx3,:,:), icnt)
+        write(logunit, '(a,i8,a)')'removed ',icnt,' locations of phantom '//trim(b3d(idx3)%var_name)
+        ! remove phantom-ice (vicen=0, aicen /=0);  do not QC vsnon=0, aicen /=0)
+        call zero_out_phantom_ice(bilin3d(idx2,:,:), bilin3d(idx1,:,:), icnt)
+        write(logunit, '(a,i8,a)')'removed ',icnt,' locations of phantom '//trim(b3d(idx1)%var_name)
+        deallocate(kmt)
+     end if
+
      rgb3d = 0.0
      do k = 1,nlevs
         if (do_ocnprep) then

sorc/ocnice_prep.fd/restarts_mod.F90

@@ -66,11 +66,11 @@ subroutine setup_icerestart(fin, fout)
        enddo
     end if
     if (allocated(c2d)) then
-      do n = 1,nconsd2d
-         vname = trim(c2d(n)%var_name)
-         call nf90_err(nf90_def_var(ncid, vname, nf90_double, (/idimid,jdimid/), varid),         &
-              'define variable: '// vname)
-      enddo
+       do n = 1,nconsd2d
+          vname = trim(c2d(n)%var_name)
+          call nf90_err(nf90_def_var(ncid, vname, nf90_double, (/idimid,jdimid/), varid),         &
+               'define variable: '// vname)
+       enddo
     end if
     if (allocated(b3d)) then
        do n = 1,nbilin3d
@@ -174,7 +174,7 @@ subroutine setup_ocnrestart(fin, fout, bathy)
                'put variable attribute: long_name')
        enddo
     end if
-     if (allocated(c2d)) then
+    if (allocated(c2d)) then
        do n = 1,nconsd2d
           vname = trim(c2d(n)%var_name)
           vunit = trim(c2d(n)%units)

sorc/ocnice_prep.fd/utils_mod.F90

@@ -47,6 +47,8 @@ module utils_mod
   public :: remap
   public :: dumpnc
   public :: nf90_err
+  public :: zero_out_land_ice
+  public :: zero_out_phantom_ice
 
 contains
   !> Pack 2D fields into arrays by mapping type
@@ -699,6 +701,59 @@ subroutine dumpnc1d(fname, vname, dims, field)
 
   end subroutine dumpnc1d
 
+  !> Reset field values to zero on land
+  !!
+  !! @param[in]     mask   the land mask
+  !! @param[inout]  fout   the field value
+  !! @param[out]    icnt   the number spatial points reset
+  !!
+  !! @author [email protected]
+  subroutine zero_out_land_ice(mask, fout, icnt)
+
+    integer(kind=4), intent(in)    :: mask(:)
+    real(kind=8),    intent(inout) :: fout(:,:)
+    integer,         intent(out)   :: icnt
+
+    !local variables
+    integer :: ij
+    !----------------------------------------------------------------------------
+
+    icnt = 0
+    do ij = 1,size(fout,2)
+       if ( mask(ij) .eq. 0 .and. sum(fout(:,ij)) .ne. 0.0)  then
+          icnt = icnt + 1
+          fout(:,ij) = 0.0
+       end if
+    end do
+  end subroutine zero_out_land_ice
+
+  !> Ensure that when fin contains zeros for all ncat, fout will also contain zeros
+  !! for all ncat
+  !!
+  !! @param[in]     fin    the field to test against
+  !! @param[inout]  fout   the field tested
+  !! @param[out]    icnt   the number spatial points reset
+  !!
+  !! @author [email protected]
+  subroutine zero_out_phantom_ice(fin, fout, icnt)
+
+    real(kind=8), intent(in)    :: fin(:,:)
+    real(kind=8), intent(inout) :: fout(:,:)
+    integer,      intent(out)   :: icnt
+
+    !local variables
+    integer :: ij
+    !----------------------------------------------------------------------------
+
+    icnt = 0
+    do ij = 1,size(fout,2)
+       if (sum(fin(:,ij)) .eq. 0.0 .and. sum(fout(:,ij)) .ne. 0.0 ) then
+          icnt = icnt + 1
+          fout(:,ij) = 0.0
+       end if
+    end do
+  end subroutine zero_out_phantom_ice
+
   !> Handle netcdf errors
   !!
   !! @param[in]  ierr        the error code
@@ -712,12 +767,13 @@ subroutine nf90_err(ierr, string)
     !----------------------------------------------------------------------------
 
     if (ierr /= nf90_noerr) then
-      write(0, '(a)') 'FATAL ERROR: ' // trim(string)// ' : ' // trim(nf90_strerror(ierr))
-      ! This fails on WCOSS2 with Intel 19 compiler. See
-      ! https://community.intel.com/t5/Intel-Fortran-Compiler/STOP-and-ERROR-STOP-with-variable-stop-codes/m-p/1182521#M149254
-      ! When WCOSS2 moves to Intel 2020+, uncomment the next line and remove stop 99
-      !stop ierr
-      stop 99
+       write(0, '(a)') 'FATAL ERROR: ' // trim(string)// ' : ' // trim(nf90_strerror(ierr))
+       ! This fails on WCOSS2 with Intel 19 compiler. See
+       ! https://community.intel.com/t5/Intel-Fortran-Compiler/STOP-and-ERROR-STOP-with-variable-stop-codes/m-p/1182521#M149254
+       ! When WCOSS2 moves to Intel 2020+, uncomment the next line and remove stop 99
+       !stop ierr
+       stop 99
     end if
   end subroutine nf90_err
+
 end module utils_mod

tests/ocnice_prep/CMakeLists.txt

@@ -36,5 +36,9 @@ execute_process( COMMAND ${CMAKE_COMMAND} -E copy
   ${CMAKE_CURRENT_SOURCE_DIR}/data/ocean.badvecgrid.csv ${CMAKE_CURRENT_BINARY_DIR}/data/ocean.badvecgrid.csv)
 
 add_executable(ocnice_ftst_program_setup ftst_program_setup.F90)
-add_test(NAME ocnice-ftst_program_setup COMMAND ocnice_ftst_program_setup)
+add_test(NAME ocnice_ftst_program_setup COMMAND ocnice_ftst_program_setup)
 target_link_libraries(ocnice_ftst_program_setup ocnice_prep_lib)
+
+add_executable(ocnice_ftst_qcice ftst_qcice.F90)
+add_test(NAME ocnice_ftst_qcice COMMAND ocnice_ftst_qcice)
+target_link_libraries(ocnice_ftst_qcice ocnice_prep_lib)