From b6024119f278c45b11cff06ffd0db1b13481fd71 Mon Sep 17 00:00:00 2001 From: Ratko Vasic Date: Fri, 20 Aug 2021 15:43:21 +0000 Subject: [PATCH 01/10] Add extrapolation for geolat and geolon to halo points. This commit references issue #105 --- .../filter_topo.fd/filter_topo.F90 | 39 +++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 b/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 index b93c77edb..fd8d40b6f 100644 --- a/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 +++ b/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 @@ -1144,6 +1147,42 @@ subroutine fill_cubic_grid_halo(data, data2, halo, ioff, joff, sign1, sign2) end subroutine fill_cubic_grid_halo + !> 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] data ??? + !! @param[in] halo ??? + !! @author Ratko + 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=1 + i_ed=npx + j_st=1 + j_ed=npy + + 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 + !> ??? !! !! @param[in] is ??? From dc72965ab44c53217193f192b4574f2ead61ee1e Mon Sep 17 00:00:00 2001 From: Ratko Vasic Date: Mon, 23 Aug 2021 15:32:12 +0000 Subject: [PATCH 02/10] Add comments for doxygen. --- sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 b/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 index fd8d40b6f..a5a18f840 100644 --- a/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 +++ b/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 @@ -1151,9 +1151,9 @@ end subroutine fill_cubic_grid_halo !! regional standalone grid. Halo points are needed for dxc and dyc !! calculation. !! - !! @param[in] data ??? - !! @param[in] halo ??? - !! @author Ratko + !! @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 From 06cc53024462cc60a33af64110c8de33416b21a2 Mon Sep 17 00:00:00 2001 From: Ratko Vasic Date: Wed, 25 Aug 2021 02:22:49 +0000 Subject: [PATCH 03/10] Fix bugs in FV3_zs_filter and two_delta_filter. --- sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 b/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 index a5a18f840..9a7c6e39d 100644 --- a/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 +++ b/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 @@ -1223,7 +1223,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 @@ -1334,6 +1334,8 @@ subroutine two_delta_filter(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles integer:: i,j, nt, t integer:: is1, ie2, js1, je2 + a2 = 0. + if ( .not. nested .and. grid_type<3 ) then is1 = max(3,is-1); ie2 = min(npx-2,ie+2) js1 = max(3,js-1); je2 = min(npy-2,je+2) From 1ec6689a21867fc4966f5a3028013b03a670db13 Mon Sep 17 00:00:00 2001 From: Ratko Vasic Date: Wed, 25 Aug 2021 20:16:14 +0000 Subject: [PATCH 04/10] Fix oro and mask fields on halo regions in regional application. --- .../filter_topo.fd/filter_topo.F90 | 30 +++++++++++++++---- 1 file changed, 24 insertions(+), 6 deletions(-) diff --git a/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 b/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 index 9a7c6e39d..910ce183d 100644 --- a/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 +++ b/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 @@ -1044,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 @@ -1159,10 +1165,10 @@ subroutine fill_regional_halo(data, halo) real, dimension(1-halo:,1-halo:,:), intent(inout) :: data integer :: h, i_st, i_ed, j_st, j_ed - i_st=1 - i_ed=npx - j_st=1 - j_ed=npy + 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 @@ -1334,8 +1340,6 @@ subroutine two_delta_filter(is, ie, js, je, isd, ied, jsd, jed, npx, npy, ntiles integer:: i,j, nt, t integer:: is1, ie2, js1, je2 - a2 = 0. - if ( .not. nested .and. grid_type<3 ) then is1 = max(3,is-1); ie2 = min(npx-2,ie+2) js1 = max(3,js-1); je2 = min(npy-2,je+2) @@ -1487,6 +1491,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 From f87f1f57666eb293221c598a8453eeea5d69745d Mon Sep 17 00:00:00 2001 From: George Gayno Date: Thu, 26 Aug 2021 15:10:03 +0000 Subject: [PATCH 05/10] Move new routine 'ftst_fill_regional_halo' to utils.F90 to allow for using it in unit tests. Create simple unit test for that routine. Fixes #105 --- .../filter_topo.fd/filter_topo.F90 | 36 ------------------- sorc/grid_tools.fd/filter_topo.fd/utils.F90 | 36 +++++++++++++++++++ tests/filter_topo/CMakeLists.txt | 4 +++ tests/filter_topo/ftst_fill_regional_halo.F90 | 23 ++++++++++++ 4 files changed, 63 insertions(+), 36 deletions(-) create mode 100644 tests/filter_topo/ftst_fill_regional_halo.F90 diff --git a/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 b/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 index 910ce183d..86f6b95ae 100644 --- a/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 +++ b/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 @@ -1153,42 +1153,6 @@ subroutine fill_cubic_grid_halo(data, data2, halo, ioff, joff, sign1, sign2) end subroutine fill_cubic_grid_halo - !> 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 - !> ??? !! !! @param[in] is ??? diff --git a/sorc/grid_tools.fd/filter_topo.fd/utils.F90 b/sorc/grid_tools.fd/filter_topo.fd/utils.F90 index 7927e769b..1444495ae 100644 --- a/sorc/grid_tools.fd/filter_topo.fd/utils.F90 +++ b/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. diff --git a/tests/filter_topo/CMakeLists.txt b/tests/filter_topo/CMakeLists.txt index 9c4501438..0928b3e24 100644 --- a/tests/filter_topo/CMakeLists.txt +++ b/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) diff --git a/tests/filter_topo/ftst_fill_regional_halo.F90 b/tests/filter_topo/ftst_fill_regional_halo.F90 new file mode 100644 index 000000000..76b553425 --- /dev/null +++ b/tests/filter_topo/ftst_fill_regional_halo.F90 @@ -0,0 +1,23 @@ +! Unit test for filter_topo routine "fill_regional_halo". +! +! Author Ratko Vasic + + program fill_halo + + use utils + + implicit none + + integer :: halo + + real, allocatable :: testdata(:,:,:) + + print*, "Starting test of filter_topo routine fill_regional_halo" + + call fill_regional_halo(testdata, halo) + + + print*, "OK" + print*, "SUCCESS!" + + end program fill_halo From 88fa6d9adb3794e6fab018acdb54c451876749cc Mon Sep 17 00:00:00 2001 From: George Gayno Date: Thu, 26 Aug 2021 17:45:40 +0000 Subject: [PATCH 06/10] Comment out call to routine in new unit test. Fixes #105 --- tests/filter_topo/ftst_fill_regional_halo.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/filter_topo/ftst_fill_regional_halo.F90 b/tests/filter_topo/ftst_fill_regional_halo.F90 index 76b553425..003f465c3 100644 --- a/tests/filter_topo/ftst_fill_regional_halo.F90 +++ b/tests/filter_topo/ftst_fill_regional_halo.F90 @@ -14,7 +14,7 @@ program fill_halo print*, "Starting test of filter_topo routine fill_regional_halo" - call fill_regional_halo(testdata, halo) +!call fill_regional_halo(testdata, halo) print*, "OK" From 852405db77aa61358b5c8821d020f19f930c5f02 Mon Sep 17 00:00:00 2001 From: Ratko Vasic Date: Wed, 1 Sep 2021 04:25:04 +0000 Subject: [PATCH 07/10] Fixed coefficient a1 in subroutine two_delta_filter. Corners were used, but not initialized. --- sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 b/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 index 86f6b95ae..d4b5d672a 100644 --- a/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 +++ b/sorc/grid_tools.fd/filter_topo.fd/filter_topo.F90 @@ -1400,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 From 9b73442a459ccac981364f1f96180be93979e98e Mon Sep 17 00:00:00 2001 From: Ratko Vasic Date: Fri, 10 Sep 2021 17:15:09 +0000 Subject: [PATCH 08/10] Add test for routine fill_regional_halo. --- tests/filter_topo/ftst_fill_regional_halo.F90 | 43 ++++++++++++++++--- 1 file changed, 37 insertions(+), 6 deletions(-) diff --git a/tests/filter_topo/ftst_fill_regional_halo.F90 b/tests/filter_topo/ftst_fill_regional_halo.F90 index 003f465c3..5c291c712 100644 --- a/tests/filter_topo/ftst_fill_regional_halo.F90 +++ b/tests/filter_topo/ftst_fill_regional_halo.F90 @@ -8,16 +8,47 @@ program fill_halo implicit none - integer :: halo + integer :: halo,ix,jx,i,j real, allocatable :: testdata(:,:,:) print*, "Starting test of filter_topo routine fill_regional_halo" -!call fill_regional_halo(testdata, halo) - - - print*, "OK" - print*, "SUCCESS!" + 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 end program fill_halo From 831344201bbf20f0e1a276c104ba88192978a733 Mon Sep 17 00:00:00 2001 From: Ratko Vasic Date: Fri, 10 Sep 2021 17:28:00 +0000 Subject: [PATCH 09/10] Fix 'module use' in test file. --- tests/filter_topo/ftst_fill_regional_halo.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/filter_topo/ftst_fill_regional_halo.F90 b/tests/filter_topo/ftst_fill_regional_halo.F90 index 5c291c712..ffb83c07c 100644 --- a/tests/filter_topo/ftst_fill_regional_halo.F90 +++ b/tests/filter_topo/ftst_fill_regional_halo.F90 @@ -4,7 +4,7 @@ program fill_halo - use utils +!use utils implicit none From bab5f9be0f41477c4f6dcfe2cebed02a939d799d Mon Sep 17 00:00:00 2001 From: Ratko Vasic Date: Fri, 10 Sep 2021 17:41:06 +0000 Subject: [PATCH 10/10] Add deallocate. --- tests/filter_topo/ftst_fill_regional_halo.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/filter_topo/ftst_fill_regional_halo.F90 b/tests/filter_topo/ftst_fill_regional_halo.F90 index ffb83c07c..de57b6b00 100644 --- a/tests/filter_topo/ftst_fill_regional_halo.F90 +++ b/tests/filter_topo/ftst_fill_regional_halo.F90 @@ -4,7 +4,7 @@ program fill_halo -!use utils + use utils implicit none @@ -51,4 +51,6 @@ program fill_halo stop 22 endif + deallocate(testdata) + end program fill_halo