filter_topo - Add extrapolation of regional fields to halo points (uf…

…s-community#571)

Fixes a bug when running the filter_topo program in regional mode.

Several variables are allocated, but their values along the halo are 
never initialized. This lead to unstable results. Add a routine
to extrapolate fields from the grid interior to the halo.

Fixes ufs-community#105.

sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90

@@ -748,6 +748,9 @@ subroutine read_grid_file(regional)
       call fill_cubic_grid_halo(geolat_c, geolat_c, ng, 1, 1, 1, 1)    
       if(.not. nested) call fill_bgrid_scalar_corners(geolon_c, ng, npx, npy, isd, jsd, XDir)
       if(.not. nested) call fill_bgrid_scalar_corners(geolat_c, ng, npx, npy, isd, jsd, YDir)
+    else
+      call fill_regional_halo(geolon_c, ng)
+      call fill_regional_halo(geolat_c, ng)
     endif
 
     !--- compute grid cell center
@@ -1041,6 +1044,12 @@ subroutine read_topo_file(regional)
       call fill_cubic_grid_halo(mask, mask, ng, 0, 0, 1, 1)
     endif
 
+    if( regional ) then
+      call fill_regional_halo(oro, ng)
+      oro(:,:,:) = max(oro(:,:,:),0.)
+      call fill_regional_halo(mask, ng)
+      mask(:,:,:) = min(max(mask(:,:,:),0.),1.)
+    endif
 
 
   end subroutine read_topo_file
@@ -1184,7 +1193,7 @@ subroutine FV3_zs_filter (is, ie, js, je, isd, ied, jsd, jed, npx, npy, npx_glob
     real, intent(IN):: sin_sg(4,isd:ied,jsd:jed,ntiles)
     real, intent(IN):: stretch_fac
     logical, intent(IN) :: nested, regional
-    real, intent(inout):: phis(isd:ied,jsd,jed,ntiles)
+    real, intent(inout):: phis(isd:ied,jsd:jed,ntiles)
     real:: cd2
     integer mdim, n_del2, n_del4
 
@@ -1391,6 +1400,16 @@ subroutine two_delta_filter(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles
                 a1(npx+1) = c3*q(npx,j,t) + c2*q(npx+1,j,t) + c1*q(npx+2,j,t)
              endif
 
+             if ( regional .and. grid_type<3 ) then
+                a1(0) = c1*q(-1,j,t) + c2*q(-1,j,t) + c3*q(0,j,t)
+                a1(2) = c3*q(1,j,t) + c2*q(2,j,t) + c1*q(3,j,t)
+                a1(1) = 0.5*(a1(0) + a1(2))
+
+                a1(npx-1) = c1*q(npx-3,j,t) + c2*q(npx-2,j,t) + c3*q(npx-1,j,t)
+                a1(npx+1) = c3*q(npx,j,t) + c2*q(npx+1,j,t) + c1*q(npx+2,j,t)
+                a1(npx) = 0.5*(a1(npx-1)+a1(npx+1))
+             endif
+
              if ( filter_type == 0 ) then
                 do i=is-1, ie+1
                    if( abs(3.*(a1(i)+a1(i+1)-2.*q(i,j,t))) > abs(a1(i)-a1(i+1)) ) then
@@ -1446,6 +1465,20 @@ subroutine two_delta_filter(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles
              enddo
           endif
 
+          if ( regional .and. grid_type<3 ) then
+             do i=is,ie
+                a2(i,0) = c1*q(i,-2,t) + c2*q(i,-1,t) + c3*q(i,0,t)
+                a2(i,2) = c3*q(i,1,t) + c2*q(i,2,t) + c1*q(i,3,t)
+                a2(i,1) = 0.5*(a2(i,0) + a2(i,2))
+             enddo
+
+             do i=is,ie
+                a2(i,npy-1) = c1*q(i,npy-3,t) + c2*q(i,npy-2,t) + c3*q(i,npy-1,t)
+                a2(i,npy+1) = c3*q(i,npy,t) + c2*q(i,npy+1,t) + c1*q(i,npy+2,t)
+                a2(i,npy) = 0.5*(a2(i,npy-1)+a2(i,npy+1))
+             enddo
+          endif
+
           if ( filter_type == 0 ) then
              do j=js-1,je+1
                 do i=is,ie

sorc/grid_tools.fd/filter_topo.fd/utils.F90

@@ -60,6 +60,42 @@ subroutine read_namelist
 
  end subroutine read_namelist
 
+!> This routine extrapolate geolat_c and geolon_c halo points for the
+!! regional standalone grid. Halo points are needed for dxc and dyc
+!! calculation.
+!!
+!! @param[in,out] data - field to be extrapolated
+!! @param[in] halo - number of halo rows/columns
+!! @author Ratko Vasic (NCEP/EMC)
+  subroutine fill_regional_halo(data, halo)
+    integer, intent(in)                               :: halo
+    real, dimension(1-halo:,1-halo:,:), intent(inout) :: data
+    integer :: h, i_st, i_ed, j_st, j_ed
+
+    i_st=lbound(data,1)+halo
+    i_ed=ubound(data,1)-halo
+    j_st=lbound(data,2)+halo
+    j_ed=ubound(data,2)-halo
+
+    do h = 1, halo
+      data(i_st:i_ed, j_st-1   , :) = 2* data(i_st:i_ed, j_st     , :) - data(i_st:i_ed, j_st+1   , :)! north
+      data(i_st:i_ed, j_ed+1   , :) = 2* data(i_st:i_ed, j_ed     , :) - data(i_st:i_ed, j_ed-1   , :)! south
+      data(i_st-1   , j_st:j_ed, :) = 2* data(i_st     , j_st:j_ed, :) - data(i_st+1   , j_st:j_ed, :)! east
+      data(i_ed+1   , j_st:j_ed, :) = 2* data(i_ed     , j_st:j_ed, :) - data(i_ed-1   , j_st:j_ed, :)! west
+
+      data(i_st-1, j_st-1, :) = (data(i_st-1, j_st, :) + data(i_st, j_st-1, :))*0.5  !NW Corner
+      data(i_ed+1, j_st-1, :) = (data(i_ed+1, j_st, :) + data(i_ed, j_st-1, :))*0.5  !NE Corner
+      data(i_st-1, j_ed+1, :) = (data(i_st-1, j_ed, :) + data(i_st, j_ed+1, :))*0.5  !SW Corner
+      data(i_ed+1, j_ed+1, :) = (data(i_ed+1, j_ed, :) + data(i_ed, j_ed+1, :))*0.5  !SE Corner
+
+      i_st=i_st-1
+      i_ed=i_ed+1
+      j_st=j_st-1
+      j_ed=j_ed+1
+    enddo
+
+  end subroutine fill_regional_halo
+
 !> Prints an error message to standard output,
 !! then halts program execution with a
 !! bad status.

tests/filter_topo/CMakeLists.txt

@@ -21,3 +21,7 @@ execute_process( COMMAND ${CMAKE_COMMAND} -E copy
 add_executable(ftst_read_filter_topo_nml ftst_readnml.F90)
 add_test(NAME filter_topo-ftst_read_namelist COMMAND ftst_read_filter_topo_nml)
 target_link_libraries(ftst_read_filter_topo_nml filter_topo_lib)
+
+add_executable(ftst_fill_regional_halo ftst_fill_regional_halo.F90)
+add_test(NAME filter_topo-ftst_fill_regional_halo COMMAND ftst_fill_regional_halo)
+target_link_libraries(ftst_fill_regional_halo filter_topo_lib)

tests/filter_topo/ftst_fill_regional_halo.F90

@@ -0,0 +1,56 @@
+! Unit test for filter_topo routine "fill_regional_halo".
+!
+! Author Ratko Vasic
+
+ program fill_halo
+
+ use utils
+
+ implicit none
+
+ integer :: halo,ix,jx,i,j
+
+ real, allocatable :: testdata(:,:,:)
+
+ print*, "Starting test of filter_topo routine fill_regional_halo"
+
+ halo = 3
+ ix=6
+ jx=9
+
+ allocate(testdata(1-halo:ix+halo,1-halo:jx+halo,1))
+
+! Initialize whole domain (including halo points) to zero
+ testdata(:,:,:)=0.
+
+! Initialize inner domain
+ do j=1,jx
+ do i=1,ix
+   testdata(i,j,1)=10*i+j
+ enddo
+ enddo
+
+!!! do j=1-halo,jx+halo
+!!!  write(0,101) testdata(1-halo:ix+halo,j,1)
+!!! enddo
+!!!101 format(12(f9.5,x))
+
+! Fill halo points
+ call fill_regional_halo(testdata, halo)
+
+!!! print*
+!!! do j=1-halo,jx+halo
+!!!  write(0,101) testdata(1-halo:ix+halo,j,1)
+!!! enddo
+
+ if(testdata(-2,-2,1)==  5.5 .and. testdata(-2,12,1)== 14.5 .and. &
+    testdata( 9,-2,1)== 65.5 .and. testdata( 9,12,1)== 74.5) then
+  print*, "OK"
+  print*, "SUCCESS!"
+ else
+  stop 22
+ endif
+
+ deallocate(testdata)
+
+ end program fill_halo