From 91c08f48f2da4e2c26ad2c7d9914c4b8e1505bbc Mon Sep 17 00:00:00 2001 From: George Gayno Date: Thu, 21 Nov 2024 20:17:02 +0000 Subject: [PATCH 01/67] Update cpld_gridgen program to create C12, C18 and C24. Fixes #1000. --- sorc/cpld_gridgen.fd/grdvars.F90 | 5 +++-- ush/cpld_gridgen.sh | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/sorc/cpld_gridgen.fd/grdvars.F90 b/sorc/cpld_gridgen.fd/grdvars.F90 index a3abad291..431b76a4f 100644 --- a/sorc/cpld_gridgen.fd/grdvars.F90 +++ b/sorc/cpld_gridgen.fd/grdvars.F90 @@ -163,8 +163,9 @@ module grdvars !! rounded to minimum_depth real(kind=real_kind), parameter :: maximum_lat = 88.0 !< The maximum latitude for water points for WW3 - integer, parameter :: nar = 6 !< the number of possible ATM resolutions - integer, parameter, dimension(nar) :: catm = (/48, 96, 192, 384, 768, 1152/) !< the ATM resolutions for mapped ocean masks + integer, parameter :: nar = 3 !< the number of possible ATM resolutions + integer, parameter, dimension(nar) :: catm = (/12, 18, 24/) !< the ATM resolutions for mapped ocean masks + contains !> Allocate grid variables diff --git a/ush/cpld_gridgen.sh b/ush/cpld_gridgen.sh index 9bb08109a..8155a3657 100755 --- a/ush/cpld_gridgen.sh +++ b/ush/cpld_gridgen.sh @@ -20,7 +20,7 @@ export RESNAME=${RESNAME:-$1} export DEBUG=.false. export MASKEDIT=.false. export DO_POSTWGTS=.true. -export MOSAICDIR_PATH=${MOSAICDIR_PATH:-$PATHTR/fix/orog} +export MOSAICDIR_PATH=${MOSAICDIR_PATH:-$PATHTR/fix/orog.lowres} export FIXDIR_PATH=${MOM6_FIXDIR}/${RESNAME} APRUN=${APRUN:-"srun"} From 20d3b1ebad49e27066d729aeca193396f00692b9 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Thu, 21 Nov 2024 20:23:04 +0000 Subject: [PATCH 02/67] Remove some unused logic from the orog program. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 index d8e55c96a..1dee1dabd 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 @@ -419,8 +419,6 @@ SUBROUTINE MAKE_MASK(zslm,slm,land_frac, & real, intent(out) :: slm(im,jm) real, intent(out) :: land_frac(im,jm) - integer, parameter :: MAXSUM=20000000 - real, parameter :: D2R = 3.14159265358979/180. integer jst, jen @@ -486,11 +484,6 @@ SUBROUTINE MAKE_MASK(zslm,slm,land_frac, & XWATR_ALL = XWATR_ALL + FLOAT(1-ZSLM(ii,jj)) XNSUM_ALL = XNSUM_ALL + 1. nsum_all = nsum_all+1 - if(nsum_all > MAXSUM) then - print*, "FATAL ERROR: nsum_all is greater than MAXSUM," - print*, "increase MAXSUM." - call ABORT() - endif if(inside_a_polygon(LONI*D2R,LATI*D2R,4, & LONO_RAD,LATO_RAD))then @@ -499,11 +492,6 @@ SUBROUTINE MAKE_MASK(zslm,slm,land_frac, & XWATR = XWATR + FLOAT(1-ZSLM(ii,jj)) XNSUM = XNSUM + 1. nsum = nsum+1 - if(nsum > MAXSUM) then - print*, "FATAL ERROR: nsum is greater than MAXSUM," - print*, "increase MAXSUM." - call ABORT() - endif endif enddo ; enddo From 4612227d78b52810088cf4c15a2f763bb6b036dd Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 22 Nov 2024 16:15:26 +0000 Subject: [PATCH 03/67] Increase MAXSUM parameter so orog code will work with C12. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 index 1dee1dabd..d4ccfb09f 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 @@ -547,7 +547,7 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, & real, intent(out) :: oro(im,jm) real, intent(out) :: var(im,jm),var4(im,jm) - integer, parameter :: MAXSUM=20000000 + integer, parameter :: MAXSUM=65000000 real, parameter :: D2R = 3.14159265358979/180. real, dimension(:), allocatable :: hgt_1d, hgt_1d_all @@ -637,7 +637,7 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, & nsum_all = nsum_all+1 if(nsum_all > MAXSUM) then print*, "FATAL ERROR: nsum_all is greater than MAXSUM," - print*, "increase MAXSUM." + print*, "increase MAXSUM.", jst,jen call ABORT() endif hgt_1d_all(nsum_all) = HEIGHT_ALL @@ -656,7 +656,7 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, & nsum = nsum+1 if(nsum > MAXSUM) then print*, "FATAL ERROR: nsum is greater than MAXSUM," - print*, "increase MAXSUM." + print*, "increase MAXSUM.", jst,jen call ABORT() endif hgt_1d(nsum) = HEIGHT From 40a1c5077c4065b2ad43200a9ffd8bb07baa54df Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 22 Nov 2024 16:18:14 +0000 Subject: [PATCH 04/67] Update ocean merge script to point to location of low res 'grid' files. Fixes #1000. --- ush/fv3gfs_ocean_merge.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ush/fv3gfs_ocean_merge.sh b/ush/fv3gfs_ocean_merge.sh index 519a2e9da..ac8798f74 100755 --- a/ush/fv3gfs_ocean_merge.sh +++ b/ush/fv3gfs_ocean_merge.sh @@ -20,7 +20,7 @@ cat << EOF > input.nml &mask_nml - ocean_mask_dir="${home_dir}/fix/orog/C${res}/ocean_mask/${ocn}/" + ocean_mask_dir="${home_dir}/fix/orog.lowres/C${res}/ocean_mask/${ocn}/" ocnres="mx${ocn}" lake_mask_dir="${TEMP_DIR}/C${res}/orog/" atmres="C${res}" From 53a97b85b9d101cd80400da1e8c9147ecefa30fb Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 22 Nov 2024 18:37:04 +0000 Subject: [PATCH 05/67] Update gdas_init utility for low res cases. Fixes #1000. --- util/gdas_init/config | 14 +++++++------- util/gdas_init/run_v16.chgres.sh | 2 +- util/gdas_init/set_fixed_files.sh | 6 ++++++ 3 files changed, 14 insertions(+), 8 deletions(-) diff --git a/util/gdas_init/config b/util/gdas_init/config index 01a98bbf7..c7bc10e15 100644 --- a/util/gdas_init/config +++ b/util/gdas_init/config @@ -51,14 +51,14 @@ # #----------------------------------------------------------- -EXTRACT_DIR=/lfs/h2/emc/stmp/$USER/gdas.init/input -EXTRACT_DATA=yes +EXTRACT_DIR=/scratch2/NCEPDEV/stmp1/George.Gayno/gdas.init/input +EXTRACT_DATA=no RUN_CHGRES=yes -yy=2020 -mm=01 -dd=01 +yy=2024 +mm=11 +dd=16 hh=06 use_v16retro=no @@ -67,12 +67,12 @@ LEVS=65 CDUMP=gfs -CRES_HIRES=C192 +CRES_HIRES=C12 CRES_ENKF=C96 UFS_DIR=$PWD/../.. -OUTDIR=/lfs/h2/emc/stmp/$USER/gdas.init/output +OUTDIR=/scratch2/NCEPDEV/stmp1/George.Gayno/gdas.init/output #--------------------------------------------------------- # Dont touch anything below here. diff --git a/util/gdas_init/run_v16.chgres.sh b/util/gdas_init/run_v16.chgres.sh index 9b95843cb..e6a96efb3 100755 --- a/util/gdas_init/run_v16.chgres.sh +++ b/util/gdas_init/run_v16.chgres.sh @@ -13,7 +13,7 @@ set -x MEMBER=$1 FIX_FV3=$UFS_DIR/fix -FIX_ORO=${FIX_FV3}/orog +FIX_ORO=${FIX_FV3}/orog.lowres FIX_AM=${FIX_FV3}/am WORKDIR=${WORKDIR:-$OUTDIR/work.${MEMBER}} diff --git a/util/gdas_init/set_fixed_files.sh b/util/gdas_init/set_fixed_files.sh index 0c3e1b4fe..832614a59 100755 --- a/util/gdas_init/set_fixed_files.sh +++ b/util/gdas_init/set_fixed_files.sh @@ -17,6 +17,12 @@ elif [ ${CTAR} == 'C768' ]; then OCNRES='025' elif [ ${CTAR} == 'C1152' ]; then OCNRES='025' +elif [ ${CTAR} == 'C12' ]; then + OCNRES='100' +elif [ ${CTAR} == 'C18' ]; then + OCNRES='100' +elif [ ${CTAR} == 'C24' ]; then + OCNRES='100' else OCNRES='025' fi From fe71fa1b5498acdc84b459e0032ad6bdcf566080 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 22 Nov 2024 20:22:42 +0000 Subject: [PATCH 06/67] Dynamically compute MAXSUM used in the computation of GWD fields. Fixes #1000. --- .../orog.fd/mtnlm7_oclsm.F90 | 24 +++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 index d4ccfb09f..0eb10c09b 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 @@ -547,7 +547,8 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, & real, intent(out) :: oro(im,jm) real, intent(out) :: var(im,jm),var4(im,jm) - integer, parameter :: MAXSUM=65000000 +! integer, parameter :: MAXSUM=65000000 + integer :: MAXSUM real, parameter :: D2R = 3.14159265358979/180. real, dimension(:), allocatable :: hgt_1d, hgt_1d_all @@ -564,8 +565,6 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, & real XL1_ALL,XS1_ALL,XW1_ALL,XW2_ALL,XW4_ALL print*,'- CREATE OROGRAPHY AND CONVEXITY.' - allocate(hgt_1d(MAXSUM)) - allocate(hgt_1d_all(MAXSUM)) !---- GLOBAL XLAT AND XLON ( DEGREE ) ! JM1 = JM - 1 @@ -578,7 +577,24 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, & GLON(I) = 0. + (I-1) * DELXN + DELXN * 0.5 ENDDO -! land_frac(:,:) = 0.0 + MAXSUM=0 + DO J=1,JM + DO I=1,IM + LONO(1) = lon_c(i,j) + LONO(2) = lon_c(i+1,j) + LONO(3) = lon_c(i+1,j+1) + LONO(4) = lon_c(i,j+1) + LATO(1) = lat_c(i,j) + LATO(2) = lat_c(i+1,j) + LATO(3) = lat_c(i+1,j+1) + LATO(4) = lat_c(i,j+1) + call get_index(IMN,JMN,4,LONO,LATO,DELXN,jst,jen,ilist,numx) + MAXSUM=MAX(MAXSUM,(JEN-JST+1)*IMN) + print*,'test point ',i,j,jst,jen,imn,((JEN-JST+1)*IMN),maxsum + ENDDO + ENDDO + allocate(hgt_1d(MAXSUM)) + allocate(hgt_1d_all(MAXSUM)) ! !---- FIND THE AVERAGE OF THE MODES IN A GRID BOX ! From 73f9e7de8b0deb94bf837964c3dd0a9ad96a4770 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 22 Nov 2024 21:33:48 +0000 Subject: [PATCH 07/67] Add strawman unit test for orog routine get_index. Fixes #1000. --- tests/orog/CMakeLists.txt | 4 ++++ tests/orog/ftst_get_index.F90 | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+) create mode 100644 tests/orog/ftst_get_index.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index 863a4d47d..97dbd4ec2 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -21,3 +21,7 @@ target_link_libraries(ftst_minmax orog_lib) add_executable(ftst_get_ll_angle ftst_get_ll_angle.F90) add_test(NAME orog-ftst_get_ll_angle COMMAND ftst_get_ll_angle) target_link_libraries(ftst_get_ll_angle orog_lib) + +add_executable(ftst_get_index ftst_get_index.F90) +add_test(NAME orog-ftst_get_index COMMAND ftst_get_index) +target_link_libraries(ftst_get_index orog_lib) diff --git a/tests/orog/ftst_get_index.F90 b/tests/orog/ftst_get_index.F90 new file mode 100644 index 000000000..d321192ff --- /dev/null +++ b/tests/orog/ftst_get_index.F90 @@ -0,0 +1,32 @@ + program test_get_index + + use orog_utils, only : get_index + + implicit none + + integer, parameter :: imn=360*120 + integer, parameter :: jmn=180*120 + integer, parameter :: npts=4 + + integer :: jst, jen, ilist(imn), numx + + real :: lono(npts), lato(npts) + real, parameter :: delxn=360.0/imn + + print*,'hello world' + + lato(1) = 0.0; lono(1) = 0.0 + lato(2) = 1.0; lono(2) = 0.5 + lato(3) = 0.0; lono(3) = 1.0 + lato(4) = -1.0; lono(4) = 0.5 + + ilist = -999 + + call get_index(imn,jmn,npts,lonO,latO,delxn,jst,jen,ilist,numx) + + print*,jst,jen,numx + + print*,ilist(1:numx) + + + end program test_get_index From 546153f19b4a89086d4755beef10931e81456a48 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 25 Nov 2024 15:22:56 +0000 Subject: [PATCH 08/67] Add more test points to ftst_get_index.F90 unit test. Fixes #1000. --- .../orog_mask_tools.fd/orog.fd/orog_utils.F90 | 3 +- tests/orog/ftst_get_index.F90 | 69 +++++++++++++++++-- 2 files changed, 65 insertions(+), 7 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 index cae1f2bec..363dfe95f 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 @@ -699,7 +699,8 @@ end subroutine remove_isolated_pts !! @param[in] jmn 'j' dimension of the high-resolution orography !! data set. !! @param[in] npts Number of vertices to describe the cubed-sphere point. -!! @param[in] lonO The longitudes of the cubed-sphere vertices. +!! @param[in] lonO The longitudes of the cubed-sphere vertices. Must +!! range from 0 - 360. !! @param[in] latO The latitudes of the cubed-sphere vertices. !! @param[in] delxn Resolution of the high-resolution orography !! data set. diff --git a/tests/orog/ftst_get_index.F90 b/tests/orog/ftst_get_index.F90 index d321192ff..70bfd6bbf 100644 --- a/tests/orog/ftst_get_index.F90 +++ b/tests/orog/ftst_get_index.F90 @@ -1,5 +1,11 @@ program test_get_index +! Unit test for routine get_index, which finds the +! i/j location of the model point on the high-resolution +! mask and orography grids. +! +! Author George Gayno NCEP/EMC + use orog_utils, only : get_index implicit none @@ -8,12 +14,15 @@ program test_get_index integer, parameter :: jmn=180*120 integer, parameter :: npts=4 - integer :: jst, jen, ilist(imn), numx + integer :: i, ii, jst, jen, ilist(imn), numx real :: lono(npts), lato(npts) - real, parameter :: delxn=360.0/imn + real, parameter :: delxn=360.0/float(imn) + + print*,"Starting test of get_index." - print*,'hello world' +! Point 1 - At equator. Western edge of grid cell at +! Greenwich. lato(1) = 0.0; lono(1) = 0.0 lato(2) = 1.0; lono(2) = 0.5 @@ -22,11 +31,59 @@ program test_get_index ilist = -999 - call get_index(imn,jmn,npts,lonO,latO,delxn,jst,jen,ilist,numx) + call get_index(imn,jmn,npts,lono,lato,delxn,jst,jen,ilist,numx) + + if (jst /= 10676) stop 2 + if (jen /= 10925) stop 4 + if (numx /= 121) stop 6 + do i = 1, numx + if (ilist(i) /= i) stop 8 + enddo + +! Point 2 - At equator. Grid cell centered at +! Greenwich. + + lato(1) = 0.0; lono(1) = -0.5 + lato(2) = 1.0; lono(2) = 0.0 + lato(3) = 0.0; lono(3) = 0.5 + lato(4) = -1.0; lono(4) = 0.0 + + ilist = -999 + + call get_index(imn,jmn,npts,lono,lato,delxn,jst,jen,ilist,numx) + + if (jst /= 10676) stop 12 + if (jen /= 10925) stop 14 + if (numx /= 121) stop 16 + ii = 1 + do i = -59, 61, 1 + if (ilist(ii) /= i) stop 18 + ii = ii + 1 + enddo + +! Point 3 - At equator. Grid cell centered at +! the dateline. + + lato(1) = -1.0; lono(1) = 179.0 + lato(2) = 1.0; lono(2) = 179.0 + lato(3) = 1.0; lono(3) = 181.0 + lato(4) = -1.0; lono(4) = 181.0 + + ilist = -999 + + call get_index(imn,jmn,npts,lono,lato,delxn,jst,jen,ilist,numx) - print*,jst,jen,numx + if (jst /= 10676) stop 22 + if (jen /= 10925) stop 24 + if (numx /= 241) stop 26 + ii = 1 + do i = 21481, 21721 + if (ilist(ii) /= i) stop 28 + ii = ii + 1 + enddo - print*,ilist(1:numx) + print*,"OK" + print*,"SUCCESS" end program test_get_index From 6936faa92a7893f13ad3cbc242dafd7c2fd07961 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 25 Nov 2024 18:18:35 +0000 Subject: [PATCH 09/67] Add more test points to ftst_get_index.F90. Fixes #1000 --- tests/orog/ftst_get_index.F90 | 46 ++++++++++++++++++++++++++++++----- 1 file changed, 40 insertions(+), 6 deletions(-) diff --git a/tests/orog/ftst_get_index.F90 b/tests/orog/ftst_get_index.F90 index 70bfd6bbf..9786085a0 100644 --- a/tests/orog/ftst_get_index.F90 +++ b/tests/orog/ftst_get_index.F90 @@ -21,8 +21,13 @@ program test_get_index print*,"Starting test of get_index." +! The get_index routine contains three 'if' blocks. The test +! points are designed to check each block. + ! Point 1 - At equator. Western edge of grid cell at -! Greenwich. +! Greenwich. This tests the third 'if' block of the routine. + + print*,"Test point 1." lato(1) = 0.0; lono(1) = 0.0 lato(2) = 1.0; lono(2) = 0.5 @@ -40,10 +45,12 @@ program test_get_index if (ilist(i) /= i) stop 8 enddo -! Point 2 - At equator. Grid cell centered at -! Greenwich. +! Point 2 - At equator. Grid cell crosses Greenwich. +! This tests the second 'if' block of the routine. + + print*,"Test point 2." - lato(1) = 0.0; lono(1) = -0.5 + lato(1) = 0.0; lono(1) = 359.5 lato(2) = 1.0; lono(2) = 0.0 lato(3) = 0.0; lono(3) = 0.5 lato(4) = -1.0; lono(4) = 0.0 @@ -56,13 +63,19 @@ program test_get_index if (jen /= 10925) stop 14 if (numx /= 121) stop 16 ii = 1 - do i = -59, 61, 1 + do i = 43141, 43200 ! East of greenwich + if (ilist(ii) /= i) stop 18 + ii = ii + 1 + enddo + do i = 1, 61 ! West of greenwich if (ilist(ii) /= i) stop 18 ii = ii + 1 enddo ! Point 3 - At equator. Grid cell centered at -! the dateline. +! the dateline. This tests the third 'if' block. + + print*,"Test point 3." lato(1) = -1.0; lono(1) = 179.0 lato(2) = 1.0; lono(2) = 179.0 @@ -82,6 +95,27 @@ program test_get_index ii = ii + 1 enddo +! Point 4 - At North Pole. Grid cell centered at +! Greenwich. This tests the first 'if' block. + + print*,"Test point 4." + + lato(1) = 89.0; lono(1) = 359.5 + lato(2) = 90.0; lono(2) = 359.5 + lato(3) = 90.0; lono(3) = 0.5 + lato(4) = 89.0; lono(4) = 0.5 + + ilist = -999 + + call get_index(imn,jmn,npts,lono,lato,delxn,jst,jen,ilist,numx) + + if (jst /= 21476) stop 32 + if (jen /= 21600) stop 34 + if (numx /= 43200) stop 36 + do i = 1, numx + if (ilist(i) /= i) stop 38 + enddo + print*,"OK" print*,"SUCCESS" From 8ca9bdf57e03654a9c94b0ec9809a240c58e2560 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 25 Nov 2024 18:32:33 +0000 Subject: [PATCH 10/67] Add some comments to unit test. Fixes #1000. --- tests/orog/ftst_get_index.F90 | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/tests/orog/ftst_get_index.F90 b/tests/orog/ftst_get_index.F90 index 9786085a0..8831f739b 100644 --- a/tests/orog/ftst_get_index.F90 +++ b/tests/orog/ftst_get_index.F90 @@ -38,11 +38,11 @@ program test_get_index call get_index(imn,jmn,npts,lono,lato,delxn,jst,jen,ilist,numx) - if (jst /= 10676) stop 2 - if (jen /= 10925) stop 4 - if (numx /= 121) stop 6 + if (jst /= 10676) stop 2 ! Starting 'j' index on hi-res grid. + if (jen /= 10925) stop 4 ! Ending 'j' index on hi-res grid. + if (numx /= 121) stop 6 ! Number of 'i' indices on the hi-res grid. do i = 1, numx - if (ilist(i) /= i) stop 8 + if (ilist(i) /= i) stop 8 ! List of 'i' indices on the hi-res grid. enddo ! Point 2 - At equator. Grid cell crosses Greenwich. @@ -59,12 +59,12 @@ program test_get_index call get_index(imn,jmn,npts,lono,lato,delxn,jst,jen,ilist,numx) - if (jst /= 10676) stop 12 - if (jen /= 10925) stop 14 - if (numx /= 121) stop 16 + if (jst /= 10676) stop 12 ! Starting 'j' index on hi-res grid. + if (jen /= 10925) stop 14 ! Ending 'j' index on hi-res grid. + if (numx /= 121) stop 16 ! Number of 'i' indices on the hi-res grid. ii = 1 do i = 43141, 43200 ! East of greenwich - if (ilist(ii) /= i) stop 18 + if (ilist(ii) /= i) stop 18 ! List of 'i' indices on the hi-res grid. ii = ii + 1 enddo do i = 1, 61 ! West of greenwich @@ -86,12 +86,12 @@ program test_get_index call get_index(imn,jmn,npts,lono,lato,delxn,jst,jen,ilist,numx) - if (jst /= 10676) stop 22 - if (jen /= 10925) stop 24 - if (numx /= 241) stop 26 + if (jst /= 10676) stop 22 ! Starting 'j' index on hi-res grid. + if (jen /= 10925) stop 24 ! Ending 'j' index on hi-res grid. + if (numx /= 241) stop 26 ! Number of 'i' indices on the hi-res grid. ii = 1 do i = 21481, 21721 - if (ilist(ii) /= i) stop 28 + if (ilist(ii) /= i) stop 28 ! List of 'i' indices on the hi-res grid. ii = ii + 1 enddo @@ -109,11 +109,11 @@ program test_get_index call get_index(imn,jmn,npts,lono,lato,delxn,jst,jen,ilist,numx) - if (jst /= 21476) stop 32 - if (jen /= 21600) stop 34 - if (numx /= 43200) stop 36 + if (jst /= 21476) stop 32 ! Starting 'j' index on hi-res grid. + if (jen /= 21600) stop 34 ! Ending 'j' index on hi-res grid. + if (numx /= 43200) stop 36 ! Number of 'i' indices on the hi-res grid. do i = 1, numx - if (ilist(i) /= i) stop 38 + if (ilist(i) /= i) stop 38 ! List of 'i' indices on the hi-res grid. enddo print*,"OK" From c15de742744bda4ec26f20a37299a933ef6b2b62 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 26 Nov 2024 19:26:50 +0000 Subject: [PATCH 11/67] Baseline rough unit test for function inside_a_polygon. Fixes #1000. --- tests/orog/CMakeLists.txt | 4 +++ tests/orog/ftst_inside_polygon.F90 | 58 ++++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+) create mode 100644 tests/orog/ftst_inside_polygon.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index 97dbd4ec2..c55cad00b 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -25,3 +25,7 @@ target_link_libraries(ftst_get_ll_angle orog_lib) add_executable(ftst_get_index ftst_get_index.F90) add_test(NAME orog-ftst_get_index COMMAND ftst_get_index) target_link_libraries(ftst_get_index orog_lib) + +add_executable(ftst_inside_polygon ftst_inside_polygon.F90) +add_test(NAME orog-ftst_inside_polygon COMMAND ftst_inside_polygon) +target_link_libraries(ftst_inside_polygon orog_lib) diff --git a/tests/orog/ftst_inside_polygon.F90 b/tests/orog/ftst_inside_polygon.F90 new file mode 100644 index 000000000..df10766f4 --- /dev/null +++ b/tests/orog/ftst_inside_polygon.F90 @@ -0,0 +1,58 @@ + program inside_polygon + + use orog_utils, only : inside_a_polygon + + implicit none + + integer, parameter :: npts=4 + + real, parameter :: D2R = 3.14159265358979/180. + logical :: inside + + real :: lon1, lat1 + real :: lon2(npts), lat2(npts) + +! Test to trip the first 'if' range check + + print*, "Test point 1" + + lon1 = 90.0 * D2R + lat1 = 0.0 * D2R + + lon2(1) = 94.0 * D2R + lat2(1) = -1.0 * D2R + lon2(2) = 94.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 95.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 95.0 * D2R + lat2(4) = -1.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 2 + +! Test to trip the second 'if' range check + + print*, "Test point 2" + + lon1 = 90.0 * D2R + lat1 = 0.0 * D2R + + lon2(1) = 84.0 * D2R + lat2(1) = -1.0 * D2R + lon2(2) = 84.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 85.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 85.0 * D2R + lat2(4) = -1.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 4 + + print*,"OK" + print*,"SUCCSSS" + + end program inside_polygon From 35c06bd6c42fb27fa335d972a686fd6932396d23 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 26 Nov 2024 21:22:09 +0000 Subject: [PATCH 12/67] Add more test points. Fixes #1000. --- tests/orog/ftst_inside_polygon.F90 | 119 +++++++++++++++++++++++++---- 1 file changed, 106 insertions(+), 13 deletions(-) diff --git a/tests/orog/ftst_inside_polygon.F90 b/tests/orog/ftst_inside_polygon.F90 index df10766f4..11c86c4b4 100644 --- a/tests/orog/ftst_inside_polygon.F90 +++ b/tests/orog/ftst_inside_polygon.F90 @@ -4,22 +4,30 @@ program inside_polygon implicit none - integer, parameter :: npts=4 + integer, parameter :: npts=4 - real, parameter :: D2R = 3.14159265358979/180. - logical :: inside + real, parameter :: D2R = 3.14159265358979/180. - real :: lon1, lat1 - real :: lon2(npts), lat2(npts) + logical :: inside -! Test to trip the first 'if' range check + real :: lon1, lat1 + real :: lon2(npts), lat2(npts) + +! The first part of the function checks if the test point is outside +! the neighborhood of the polygon - i.e., a gross check. There +! are six separate checks. The first six tests are designed +! so that the polygon is far enough from the test point that +! each check is tripped. + +! Test to trip the first 'if' gross check. Is point outside the +! neighborhood in the plus 'x' direction? print*, "Test point 1" - lon1 = 90.0 * D2R + lon1 = 90.0 * D2R ! Test point. lat1 = 0.0 * D2R - lon2(1) = 94.0 * D2R + lon2(1) = 94.0 * D2R ! Polygon. lat2(1) = -1.0 * D2R lon2(2) = 94.0 * D2R lat2(2) = 1.0 * D2R @@ -30,16 +38,17 @@ program inside_polygon inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) - if (inside) stop 2 + if (inside) stop 2 ! Test point should be outside polygon. -! Test to trip the second 'if' range check +! Test to trip the second 'if' gross check. Is point outside the +! neighborhood in the minus 'x' direction? print*, "Test point 2" - lon1 = 90.0 * D2R + lon1 = 90.0 * D2R ! Test point. lat1 = 0.0 * D2R - lon2(1) = 84.0 * D2R + lon2(1) = 84.0 * D2R ! Polygon. lat2(1) = -1.0 * D2R lon2(2) = 84.0 * D2R lat2(2) = 1.0 * D2R @@ -50,7 +59,91 @@ program inside_polygon inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) - if (inside) stop 4 + if (inside) stop 4 ! Test point should be outside polygon. + +! Test to trip the third 'if' gross check. Is point outside the +! neighborhood in the plus 'y' direction? + + print*, "Test point 3" + + lon1 = 0.0 * D2R ! Test point. + lat1 = 0.0 * D2R + + lon2(1) = 354.0 * D2R ! Polygon. + lat2(1) = -1.0 * D2R + lon2(2) = 354.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 355.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 355.0 * D2R + lat2(4) = -1.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 6 ! Test point should be outside polygon. + +! Test to trip the fourth 'if' gross check. Is point outside the +! neighborhood in the minus 'y' direction? + + print*, "Test point 4" + + lon1 = 0.0 * D2R ! Test point. + lat1 = 0.0 * D2R + + lon2(1) = 4.0 * D2R ! Polygon. + lat2(1) = -1.0 * D2R + lon2(2) = 4.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 5.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 5.0 * D2R + lat2(4) = -1.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 8 ! Test point should be outside polygon. + +! Test to trip the fifth 'if' gross check. Is point outside the +! neighborhood in the plus 'z' direction? + + print*, "Test point 5" + + lon1 = 0.0 * D2R ! Test point. + lat1 = 0.0 * D2R + + lon2(1) = 359.0 * D2R ! Polygon. + lat2(1) = -6.0 * D2R + lon2(2) = 359.0 * D2R + lat2(2) = -5.0 * D2R + lon2(3) = 1.0 * D2R + lat2(3) = -5.0 * D2R + lon2(4) = 1.0 * D2R + lat2(4) = -6.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 10 ! Test point should be outside polygon. + +! Test to trip the sixth 'if' gross check. Is point outside the +! neighborhood in the minus 'z' direction? + + print*, "Test point 6" + + lon1 = 0.0 * D2R ! Test point. + lat1 = 0.0 * D2R + + lon2(1) = 359.0 * D2R ! Polygon. + lat2(1) = 5.0 * D2R + lon2(2) = 359.0 * D2R + lat2(2) = 6.0 * D2R + lon2(3) = 1.0 * D2R + lat2(3) = 6.0 * D2R + lon2(4) = 1.0 * D2R + lat2(4) = 5.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 12 ! Test point should be outside polygon. print*,"OK" print*,"SUCCSSS" From cc477f40247e1e7028df6e93cf6a4a142c568031 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 2 Dec 2024 19:55:38 +0000 Subject: [PATCH 13/67] Add point to test the case when the test point coincides with one of the polygon's corner points. Fixes #1000. --- tests/orog/ftst_inside_polygon.F90 | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/orog/ftst_inside_polygon.F90 b/tests/orog/ftst_inside_polygon.F90 index 11c86c4b4..c93cf8956 100644 --- a/tests/orog/ftst_inside_polygon.F90 +++ b/tests/orog/ftst_inside_polygon.F90 @@ -145,6 +145,27 @@ program inside_polygon if (inside) stop 12 ! Test point should be outside polygon. +! Test the case when the test point coincides with a corner point +! of the polygon. Result should be 'inside'. + + print*, "Test point 7" + + lon1 = 90.0 * D2R ! Test point. + lat1 = 0.0 * D2R + + lon2(1) = 90.0 * D2R ! Polygon. + lat2(1) = 0.0 * D2R + lon2(2) = 90.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 91.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 91.0 * D2R + lat2(4) = 0.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (.not.inside) stop 14 ! Test point should be inside polygon. + print*,"OK" print*,"SUCCSSS" From cfe59e188ee5d0f72b52dd2d8bf5f59ac7604587 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 2 Dec 2024 20:29:55 +0000 Subject: [PATCH 14/67] Add a test point located in the middle of the polygon. Fixes #1000. --- tests/orog/ftst_inside_polygon.F90 | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/tests/orog/ftst_inside_polygon.F90 b/tests/orog/ftst_inside_polygon.F90 index c93cf8956..94c2d850e 100644 --- a/tests/orog/ftst_inside_polygon.F90 +++ b/tests/orog/ftst_inside_polygon.F90 @@ -164,7 +164,30 @@ program inside_polygon inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) - if (.not.inside) stop 14 ! Test point should be inside polygon. + if (.not.inside) stop 14 ! Test point should be inside polygon. + +! Test the case when the test point is inside the polygon. +! Here, the test point is exactly in the middle of a square +! polygon. That means the computed angle to each corner +! is 90 degrees. + + print*, "Test point 8" + + lon1 = 90.5 * D2R ! Test point. + lat1 = 0.5 * D2R + + lon2(1) = 90.0 * D2R ! Polygon. + lat2(1) = 0.0 * D2R + lon2(2) = 90.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 91.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 91.0 * D2R + lat2(4) = 0.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (.not.inside) stop 16 ! Test point should be inside polygon. print*,"OK" print*,"SUCCSSS" From e8637dbee32253c05e99b4298f068604a5f1a618 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 2 Dec 2024 20:56:40 +0000 Subject: [PATCH 15/67] Finish unit test for inside_a_polygon. Fixes #1000. --- tests/orog/ftst_inside_polygon.F90 | 42 +++++++++++++++++++++++++----- 1 file changed, 35 insertions(+), 7 deletions(-) diff --git a/tests/orog/ftst_inside_polygon.F90 b/tests/orog/ftst_inside_polygon.F90 index 94c2d850e..b71e57c38 100644 --- a/tests/orog/ftst_inside_polygon.F90 +++ b/tests/orog/ftst_inside_polygon.F90 @@ -1,5 +1,11 @@ program inside_polygon +! Unit test for function inside_a_polygon, which determines +! whether a test point is located within an area specified +! by a polygon. +! +! Author George Gayno NCEP/EMC + use orog_utils, only : inside_a_polygon implicit none @@ -13,13 +19,15 @@ program inside_polygon real :: lon1, lat1 real :: lon2(npts), lat2(npts) -! The first part of the function checks if the test point is outside + print*,"Starting test of inside_a_polygon." + +! The first part of the function checks if the test point is well outside ! the neighborhood of the polygon - i.e., a gross check. There ! are six separate checks. The first six tests are designed ! so that the polygon is far enough from the test point that ! each check is tripped. -! Test to trip the first 'if' gross check. Is point outside the +! Test to trip the first 'if' gross check. Is point well outside the ! neighborhood in the plus 'x' direction? print*, "Test point 1" @@ -40,7 +48,7 @@ program inside_polygon if (inside) stop 2 ! Test point should be outside polygon. -! Test to trip the second 'if' gross check. Is point outside the +! Test to trip the second 'if' gross check. Is point well outside the ! neighborhood in the minus 'x' direction? print*, "Test point 2" @@ -61,7 +69,7 @@ program inside_polygon if (inside) stop 4 ! Test point should be outside polygon. -! Test to trip the third 'if' gross check. Is point outside the +! Test to trip the third 'if' gross check. Is point well outside the ! neighborhood in the plus 'y' direction? print*, "Test point 3" @@ -82,7 +90,7 @@ program inside_polygon if (inside) stop 6 ! Test point should be outside polygon. -! Test to trip the fourth 'if' gross check. Is point outside the +! Test to trip the fourth 'if' gross check. Is point well outside the ! neighborhood in the minus 'y' direction? print*, "Test point 4" @@ -103,7 +111,7 @@ program inside_polygon if (inside) stop 8 ! Test point should be outside polygon. -! Test to trip the fifth 'if' gross check. Is point outside the +! Test to trip the fifth 'if' gross check. Is point well outside the ! neighborhood in the plus 'z' direction? print*, "Test point 5" @@ -124,7 +132,7 @@ program inside_polygon if (inside) stop 10 ! Test point should be outside polygon. -! Test to trip the sixth 'if' gross check. Is point outside the +! Test to trip the sixth 'if' gross check. Is point well outside the ! neighborhood in the minus 'z' direction? print*, "Test point 6" @@ -189,6 +197,26 @@ program inside_polygon if (.not.inside) stop 16 ! Test point should be inside polygon. +! Test the case when the test point just outside the polygon. + + print*, "Test point 9" + + lon1 = 89.5 * D2R ! Test point. + lat1 = 0.5 * D2R + + lon2(1) = 90.0 * D2R ! Polygon. + lat2(1) = 0.0 * D2R + lon2(2) = 90.0 * D2R + lat2(2) = 1.0 * D2R + lon2(3) = 91.0 * D2R + lat2(3) = 1.0 * D2R + lon2(4) = 91.0 * D2R + lat2(4) = 0.0 * D2R + + inside=inside_a_polygon(lon1, lat1, npts, lon2, lat2) + + if (inside) stop 18 ! Test point should be outside polygon. + print*,"OK" print*,"SUCCSSS" From dc516b6a2d979f417dedc48adef205414ed85fb6 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 3 Dec 2024 09:10:30 -0600 Subject: [PATCH 16/67] Baseline a preliminary unit test for the two 'transpose' routines. Fixes #1000. --- tests/orog/CMakeLists.txt | 4 ++ tests/orog/ftst_transpose.F90 | 69 +++++++++++++++++++++++++++++++++++ 2 files changed, 73 insertions(+) create mode 100644 tests/orog/ftst_transpose.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index c55cad00b..b5af3c09f 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -29,3 +29,7 @@ target_link_libraries(ftst_get_index orog_lib) add_executable(ftst_inside_polygon ftst_inside_polygon.F90) add_test(NAME orog-ftst_inside_polygon COMMAND ftst_inside_polygon) target_link_libraries(ftst_inside_polygon orog_lib) + +add_executable(ftst_transpose ftst_transpose.F90) +add_test(NAME orog-ftst_transpose COMMAND ftst_transpose) +target_link_libraries(ftst_transpose orog_lib) diff --git a/tests/orog/ftst_transpose.F90 b/tests/orog/ftst_transpose.F90 new file mode 100644 index 000000000..6704de16e --- /dev/null +++ b/tests/orog/ftst_transpose.F90 @@ -0,0 +1,69 @@ + program transpose + + use orog_utils, only : transpose_mask + + implicit none + + integer, parameter :: imn = 360 + integer, parameter :: jmn = 181 + + integer(1) :: mask(imn,jmn) + integer :: i, ii, j, jj + + print*,"Starting test of transpose routines." + +! Transpose is from S to N to the NCEP standard N to S, +! and from the dateline to the NCEP standard Greenwich. + +! Set up a one-degree global mask. Although mask is a +! yes/no flag, for this test set each row to the +! latitude to simplify checking the answer. + + jj=0 + do j = -90, 90 ! row 1 is South Pole. + jj = jj + 1 + mask(:,jj) = j + enddo + + call transpose_mask(imn, jmn, mask) + + jj=0 + do j = 90, -90, -1 ! row 1 is North Pole. + jj = jj + 1 + do i = 1, imn + if (mask(i,jj) /= j) stop 2 + enddo + enddo + +! Now test the transpose in the E/W direction. +! Here, the East half of the domain is a flag value +! of minus 1 and the West half is plus 1. + + do i = 1, 180 + mask(i,:) = -1 + enddo + do i = 181, 360 + mask(i,:) = +1 + enddo + + call transpose_mask(imn, jmn, mask) + +! After the transpose, the East half should be plus 1 +! and the West half should be minus 1. + + do i = 1, 180 + do j = 1, jmn + if (mask(i,j) /= 1) stop 4 + enddo + enddo + do i = 181, 360 + do j = 1, jmn + if (mask(i,j) /= -1) stop 6 + enddo + enddo + + print*,"OK" + + print*,"SUCCESS" + + end program transpose From 1c21bcff53d7a85292c9b3f4d0e0c6da9098e3b6 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 3 Dec 2024 09:41:49 -0600 Subject: [PATCH 17/67] Update unit test to also call the transpose_orog routine. Fixes #1000. --- tests/orog/ftst_transpose.F90 | 31 ++++++++++++++++++++++++++----- 1 file changed, 26 insertions(+), 5 deletions(-) diff --git a/tests/orog/ftst_transpose.F90 b/tests/orog/ftst_transpose.F90 index 6704de16e..5d62a7a2e 100644 --- a/tests/orog/ftst_transpose.F90 +++ b/tests/orog/ftst_transpose.F90 @@ -1,6 +1,13 @@ program transpose - use orog_utils, only : transpose_mask +! Unit test for routines transpose_mask and transpose_orog. +! They are essentially the same routine - one takes +! one byte integer data, and one takes two byte integer +! data. +! +! Author George Gayno NCEP/EMC + + use orog_utils, only : transpose_mask, transpose_orog implicit none @@ -8,30 +15,37 @@ program transpose integer, parameter :: jmn = 181 integer(1) :: mask(imn,jmn) + integer(2) :: mask2(imn,jmn) integer :: i, ii, j, jj print*,"Starting test of transpose routines." ! Transpose is from S to N to the NCEP standard N to S, -! and from the dateline to the NCEP standard Greenwich. +! and, in the E/W direction, from the dateline to the +! NCEP standard Greenwich. -! Set up a one-degree global mask. Although mask is a -! yes/no flag, for this test set each row to the -! latitude to simplify checking the answer. +! Test N/S transpose. Set up a one-degree global mask. +! Although mask is a yes/no flag, for this test set each +! row to the latitude to simplify checking the answer. jj=0 do j = -90, 90 ! row 1 is South Pole. jj = jj + 1 mask(:,jj) = j + mask2(:,jj) = j enddo + print*,"Call transpose routines to test N/S flip." + call transpose_mask(imn, jmn, mask) + call transpose_orog(imn, jmn, mask2) jj=0 do j = 90, -90, -1 ! row 1 is North Pole. jj = jj + 1 do i = 1, imn if (mask(i,jj) /= j) stop 2 + if (mask2(i,jj) /= j) stop 3 enddo enddo @@ -41,12 +55,17 @@ program transpose do i = 1, 180 mask(i,:) = -1 + mask2(i,:) = -1 enddo do i = 181, 360 mask(i,:) = +1 + mask2(i,:) = +1 enddo + print*,"Call transpose routines to test E/W flip." + call transpose_mask(imn, jmn, mask) + call transpose_orog(imn, jmn, mask2) ! After the transpose, the East half should be plus 1 ! and the West half should be minus 1. @@ -54,11 +73,13 @@ program transpose do i = 1, 180 do j = 1, jmn if (mask(i,j) /= 1) stop 4 + if (mask2(i,j) /= 1) stop 5 enddo enddo do i = 181, 360 do j = 1, jmn if (mask(i,j) /= -1) stop 6 + if (mask2(i,j) /= -1) stop 7 enddo enddo From 92bac6847f460eeac9e25e1890d8334ffbb3e0f2 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 3 Dec 2024 12:45:10 -0600 Subject: [PATCH 18/67] Add unit test for routine find_poles. Fixes #1000. --- tests/orog/CMakeLists.txt | 4 ++ tests/orog/ftst_find_poles.F90 | 69 ++++++++++++++++++++++++++++++++++ 2 files changed, 73 insertions(+) create mode 100644 tests/orog/ftst_find_poles.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index b5af3c09f..8248d91bf 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -33,3 +33,7 @@ target_link_libraries(ftst_inside_polygon orog_lib) add_executable(ftst_transpose ftst_transpose.F90) add_test(NAME orog-ftst_transpose COMMAND ftst_transpose) target_link_libraries(ftst_transpose orog_lib) + +add_executable(ftst_find_poles ftst_find_poles.F90) +add_test(NAME orog-ftst_find_poles COMMAND ftst_find_poles) +target_link_libraries(ftst_find_poles orog_lib) diff --git a/tests/orog/ftst_find_poles.F90 b/tests/orog/ftst_find_poles.F90 new file mode 100644 index 000000000..3d98cb972 --- /dev/null +++ b/tests/orog/ftst_find_poles.F90 @@ -0,0 +1,69 @@ + program ftst_find_poles + +! Unit test for routine find_poles. +! +! Author George Gayno NCEP/EMC + + use orog_utils, only : find_poles + + implicit none + + integer, parameter :: nx=9 + integer, parameter :: ny=9 + + integer :: i_north_pole, j_north_pole + integer :: i_south_pole, j_south_pole + + real :: geolat(nx+1,ny+1) + + print*,"Starting test of find_poles." + +! Point 1 - North Pole at (4,5) + + print*,"North pole test." + + geolat = 89.0 + geolat(4,5) = 89.95 ! North pole + + call find_poles(geolat, nx, ny, i_north_pole, j_north_pole, & + i_south_pole, j_south_pole) + + if (i_north_pole /= 4) stop 2 + if (j_north_pole /= 5) stop 4 + if (i_south_pole /= 0) stop 6 + if (j_south_pole /= 0) stop 8 + +! Point 2 - South Pole at (2,8) + + print*,"South pole test." + + geolat = -89.0 + geolat(2,8) = -89.95 ! South pole + + call find_poles(geolat, nx, ny, i_north_pole, j_north_pole, & + i_south_pole, j_south_pole) + + if (i_north_pole /= 0) stop 12 + if (j_north_pole /= 0) stop 14 + if (i_south_pole /= 2) stop 16 + if (j_south_pole /= 8) stop 18 + +! Point 3 - No pole points. + + print*,"No pole test." + + geolat = -89.3 + + call find_poles(geolat, nx, ny, i_north_pole, j_north_pole, & + i_south_pole, j_south_pole) + + if (i_north_pole /= 0) stop 22 + if (j_north_pole /= 0) stop 24 + if (i_south_pole /= 0) stop 26 + if (j_south_pole /= 0) stop 28 + + print*,"OK" + + print*,"SUCCESS" + + end program ftst_find_poles From 58cd3bc6bcbc32f247da56ea9bfb62cbd4d3f209 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 3 Dec 2024 14:54:39 -0600 Subject: [PATCH 19/67] Add partial unit test for routine find_nearest_pole_points. Fixes #1000. --- tests/orog/CMakeLists.txt | 4 + tests/orog/ftst_find_nearest_pole_pts.F90 | 108 ++++++++++++++++++++++ 2 files changed, 112 insertions(+) create mode 100644 tests/orog/ftst_find_nearest_pole_pts.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index 8248d91bf..2a46e4f3e 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -37,3 +37,7 @@ target_link_libraries(ftst_transpose orog_lib) add_executable(ftst_find_poles ftst_find_poles.F90) add_test(NAME orog-ftst_find_poles COMMAND ftst_find_poles) target_link_libraries(ftst_find_poles orog_lib) + +add_executable(ftst_find_nearest_pole_pts ftst_find_nearest_pole_pts.F90) +add_test(NAME orog-ftst_find_nearest_pole_pts COMMAND ftst_find_nearest_pole_pts) +target_link_libraries(ftst_find_nearest_pole_pts orog_lib) diff --git a/tests/orog/ftst_find_nearest_pole_pts.F90 b/tests/orog/ftst_find_nearest_pole_pts.F90 new file mode 100644 index 000000000..b0da658da --- /dev/null +++ b/tests/orog/ftst_find_nearest_pole_pts.F90 @@ -0,0 +1,108 @@ + program find_nearest_pole_pts + + use orog_utils, only : find_nearest_pole_points + + implicit none + + integer, parameter :: im=48 + integer, parameter :: jm=48 + + integer :: i, j + integer :: i_north_pole, j_north_pole + integer :: i_south_pole, j_south_pole + + logical :: is_north_pole(im,jm) + logical :: is_south_pole(im,jm) + +! Test 1 - C48 uniform tile containing north pole. + + i_north_pole = 49 ! supergrid index + j_north_pole = 49 ! uniform grid + i_south_pole = 0 + j_south_pole = 0 + + call find_nearest_pole_points(i_north_pole, j_north_pole, & + i_south_pole, j_south_pole, im, jm, is_north_pole, & + is_south_pole) + + do j = 1, im + do i = 1, jm + if((i == 24 .and. j == 24) .or. & + (i == 24 .and. j == 25) .or. & + (i == 25 .and. j == 24) .or. & + (i == 25 .and. j == 25)) then + if (.not.is_north_pole(i,j)) stop 2 + else + if (is_north_pole(i,j)) stop 4 + endif + if (is_south_pole(i,j)) stop 8 + enddo + enddo + +! Test 2 - C48 uniform tile containing south pole. + + i_north_pole = 0 + j_north_pole = 0 + i_south_pole = 49 + j_south_pole = 49 + + call find_nearest_pole_points(i_north_pole, j_north_pole, & + i_south_pole, j_south_pole, im, jm, is_north_pole, & + is_south_pole) + + do j = 1, im + do i = 1, jm + if((i == 24 .and. j == 24) .or. & + (i == 24 .and. j == 25) .or. & + (i == 25 .and. j == 24) .or. & + (i == 25 .and. j == 25)) then + if (.not.is_south_pole(i,j)) stop 12 + else + if (is_south_pole(i,j)) stop 14 + endif + if (is_north_pole(i,j)) stop 18 + enddo + enddo + +! Test 3 - C48 uniform tile containing no pole. + + i_north_pole = 0 + j_north_pole = 0 + i_south_pole = 0 + j_south_pole = 0 + + call find_nearest_pole_points(i_north_pole, j_north_pole, & + i_south_pole, j_south_pole, im, jm, is_north_pole, & + is_south_pole) + + do j = 1, im + do i = 1, jm + if (is_south_pole(i,j)) stop 24 + if (is_north_pole(i,j)) stop 26 + enddo + enddo + +! Test 4 - C48 stretched grid tile containing south pole. + + i_north_pole = 0 + j_north_pole = 0 + i_south_pole = 10 + j_south_pole = 49 + + call find_nearest_pole_points(i_north_pole, j_north_pole, & + i_south_pole, j_south_pole, im, jm, is_north_pole, & + is_south_pole) + + do j = 1, im + do i = 1, jm + if((i == 5 .and. j == 24) .or. & + (i == 5 .and. j == 25)) then + if (.not.is_south_pole(i,j)) stop 32 + else + if (is_south_pole(i,j)) stop 34 + endif + if (is_north_pole(i,j)) stop 38 + enddo + enddo + + end program find_nearest_pole_pts From c8dad116266ce555e788b454527f7e163929fa44 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Thu, 5 Dec 2024 08:23:28 -0600 Subject: [PATCH 20/67] Complete unit test for "file_nearest_pole_points". Fixes #1000. --- tests/orog/ftst_find_nearest_pole_pts.F90 | 58 +++++++++++++++++++++-- 1 file changed, 53 insertions(+), 5 deletions(-) diff --git a/tests/orog/ftst_find_nearest_pole_pts.F90 b/tests/orog/ftst_find_nearest_pole_pts.F90 index b0da658da..4c6844790 100644 --- a/tests/orog/ftst_find_nearest_pole_pts.F90 +++ b/tests/orog/ftst_find_nearest_pole_pts.F90 @@ -1,5 +1,14 @@ program find_nearest_pole_pts +! Unit test for routine find_nearest_pole_points. +! +! The routine is passed the supergrid indices of +! the pole (from the 'grid' files) and returns the +! nearest pole points (using a logical array) on +! the model tile. +! +! Author George Gayno NCEP/EMC + use orog_utils, only : find_nearest_pole_points implicit none @@ -14,10 +23,14 @@ program find_nearest_pole_pts logical :: is_north_pole(im,jm) logical :: is_south_pole(im,jm) + print*,"Starting test of find_nearest_pole_points routine." + ! Test 1 - C48 uniform tile containing north pole. - i_north_pole = 49 ! supergrid index - j_north_pole = 49 ! uniform grid + print*,"Test number 1." + + i_north_pole = 49 ! Supergrid index + j_north_pole = 49 i_south_pole = 0 j_south_pole = 0 @@ -41,9 +54,11 @@ program find_nearest_pole_pts ! Test 2 - C48 uniform tile containing south pole. + print*,"Test number 2." + i_north_pole = 0 j_north_pole = 0 - i_south_pole = 49 + i_south_pole = 49 ! Supergrid index. j_south_pole = 49 call find_nearest_pole_points(i_north_pole, j_north_pole, & @@ -66,7 +81,9 @@ program find_nearest_pole_pts ! Test 3 - C48 uniform tile containing no pole. - i_north_pole = 0 + print*,"Test number 3." + + i_north_pole = 0 ! Zero indicates no pole. j_north_pole = 0 i_south_pole = 0 j_south_pole = 0 @@ -84,9 +101,11 @@ program find_nearest_pole_pts ! Test 4 - C48 stretched grid tile containing south pole. + print*,"Test number 4." + i_north_pole = 0 j_north_pole = 0 - i_south_pole = 10 + i_south_pole = 10 ! Supergrid index. j_south_pole = 49 call find_nearest_pole_points(i_north_pole, j_north_pole, & @@ -105,4 +124,33 @@ program find_nearest_pole_pts enddo enddo +! Test 5 - C48 stretched grid tile containing north pole. + + print*,"Test number 5." + + i_north_pole = 62 ! Supergrid index. + j_north_pole = 49 + i_south_pole = 0 + j_south_pole = 0 + + call find_nearest_pole_points(i_north_pole, j_north_pole, & + i_south_pole, j_south_pole, im, jm, is_north_pole, & + is_south_pole) + + do j = 1, im + do i = 1, jm + if((i == 31 .and. j == 24) .or. & + (i == 31 .and. j == 25)) then + if (.not.is_north_pole(i,j)) stop 32 + else + if (is_north_pole(i,j)) stop 34 + endif + if (is_south_pole(i,j)) stop 38 + enddo + enddo + + print*,"OK" + + print*,"SUCCESS" + end program find_nearest_pole_pts From 4df1f51316ccc60479dd01197bfb7621e5381294 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 13 Dec 2024 09:34:59 -0600 Subject: [PATCH 21/67] Baseline a test for function get_xnsum. Fixes #1000. --- tests/orog/CMakeLists.txt | 4 ++ tests/orog/ftst_get_xnsum.F90 | 88 +++++++++++++++++++++++++++++++++++ 2 files changed, 92 insertions(+) create mode 100644 tests/orog/ftst_get_xnsum.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index 2a46e4f3e..5cf771cd8 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -41,3 +41,7 @@ target_link_libraries(ftst_find_poles orog_lib) add_executable(ftst_find_nearest_pole_pts ftst_find_nearest_pole_pts.F90) add_test(NAME orog-ftst_find_nearest_pole_pts COMMAND ftst_find_nearest_pole_pts) target_link_libraries(ftst_find_nearest_pole_pts orog_lib) + +add_executable(ftst_get_xnsum ftst_get_xnsum.F90) +add_test(NAME orog-ftst_get_xnsum COMMAND ftst_get_xnsum) +target_link_libraries(ftst_get_xnsum orog_lib) diff --git a/tests/orog/ftst_get_xnsum.F90 b/tests/orog/ftst_get_xnsum.F90 new file mode 100644 index 000000000..399b3a8b7 --- /dev/null +++ b/tests/orog/ftst_get_xnsum.F90 @@ -0,0 +1,88 @@ + program check_get_xnsum + + use orog_utils, only : get_xnsum + + implicit none + + integer, parameter :: imn=360 + integer, parameter :: jmn=181 + + integer :: j + integer :: zavg(imn,jmn) + integer :: zslm(imn,jmn) + + real :: delxn=360.0/float(imn) + real :: glat(jmn) + real :: lon1,lat1,lon2,lat2 + real :: xnsum + + print*,"Begin test." + +! Set up a 'high-res' 1-degree grid. + + do j = 1, jmn + glat(j) = -90.0 + float(j-1) * delxn + print*,'j lat ',j,glat(j) + enddo + +! First test point. The high-res grid is all ocean. Since all points in +! the model grid box will be sea level, there will be no points +! higher than the average of 0 meters. + + print*,"Test point 1." + + zslm = 0 ! all water + zavg = -999 ! all sea level + +! Bounds of model grid box - straddles greenwich. + + lon1 = -2.5 + lon2 = 2.5 + lat1 = -1.5 + lat2 = 1.5 + + xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & + glat, zavg, zslm, delxn) + + if (nint(xnsum) /= 0) stop 2 + + print*,"Test point 2." + + zslm = 1 ! all land + zavg = 50 ! constant elevation of 50 meters. + +! Bounds of model grid box - straddles greenwich. + + lon1 = -2.5 + lon2 = 2.5 + lat1 = -1.5 + lat2 = 1.5 + + xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & + glat, zavg, zslm, delxn) + + if (nint(xnsum) /= 0) stop 4 + + print*,"Test point 3." + + zslm = 1 ! all land + zavg = 50 ! constant elevation of 50 meters. + zavg(359,91) = 100 + +! Bounds of model grid box - straddles greenwich. + + lon1 = -2.5 + lon2 = 2.5 + lat1 = -1.5 + lat2 = 1.5 + + xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & + glat, zavg, zslm, delxn) + + if (nint(xnsum) /= 1) stop 6 + + print*,"OK" + + print*,"SUCCESS" + + end program check_get_xnsum From d39963a757bd604068ec05a4b783eb6ed8798b12 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 13 Dec 2024 10:29:13 -0600 Subject: [PATCH 22/67] Add more test points to get_xnsum test. Fixes #1000. --- tests/orog/ftst_get_xnsum.F90 | 41 +++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/tests/orog/ftst_get_xnsum.F90 b/tests/orog/ftst_get_xnsum.F90 index 399b3a8b7..0781d214c 100644 --- a/tests/orog/ftst_get_xnsum.F90 +++ b/tests/orog/ftst_get_xnsum.F90 @@ -81,6 +81,47 @@ program check_get_xnsum if (nint(xnsum) /= 1) stop 6 + print*,"Test point 4." + + zslm = 0 ! all water + zavg = -999 ! set to sea level. + zavg(2,93) = 20 ! this represents an inland lake above sea level. + +! Bounds of model grid box - straddles greenwich. + + lon1 = -2.5 + lon2 = 2.5 + lat1 = -1.5 + lat2 = 1.5 + + xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & + glat, zavg, zslm, delxn) + + if (nint(xnsum) /= 1) stop 8 + + print*,"Test point 5." + + zslm = 0 ! all water + zavg = -999 ! set to sea level. + zslm(1:2,90:93) = 1 ! half points in grid box land + zavg(1,90:93) = 25 ! land points above sea level + zavg(2,90) = 100 ! land points above sea level + zavg(2,91) = 110 ! land points above sea level + zavg(2,92) = 107 ! land points above sea level + zavg(2,93) = 207 ! land points above sea level + +! Bounds of model grid box - straddles greenwich. + + lon1 = -2.5 + lon2 = 2.5 + lat1 = -1.5 + lat2 = 1.5 + + xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & + glat, zavg, zslm, delxn) + + if (nint(xnsum) /= 4) stop 10 + print*,"OK" print*,"SUCCESS" From 4f4de62ea068aca3ed6228de0384da4952800e69 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 13 Dec 2024 13:46:09 -0600 Subject: [PATCH 23/67] Finish unit test for function get_xnsum. Fixes #1000. --- tests/orog/ftst_get_xnsum.F90 | 116 ++++++++++++++++------------------ 1 file changed, 55 insertions(+), 61 deletions(-) diff --git a/tests/orog/ftst_get_xnsum.F90 b/tests/orog/ftst_get_xnsum.F90 index 0781d214c..0de97a43b 100644 --- a/tests/orog/ftst_get_xnsum.F90 +++ b/tests/orog/ftst_get_xnsum.F90 @@ -1,121 +1,115 @@ program check_get_xnsum +! Unit test for function get_xnsum, which computes the +! the number of high-resolution orography data points +! within a model grid box that are above the average +! terrain height. + use orog_utils, only : get_xnsum implicit none - integer, parameter :: imn=360 - integer, parameter :: jmn=181 + integer, parameter :: imn=360 ! i-dimension of high-res grid + integer, parameter :: jmn=181 ! j-dimension of high-res grid integer :: j - integer :: zavg(imn,jmn) - integer :: zslm(imn,jmn) + integer :: zavg(imn,jmn) ! High-res grid terrain + integer :: zslm(imn,jmn) ! High-res grid land mask - real :: delxn=360.0/float(imn) - real :: glat(jmn) + real :: delxn=360.0/float(imn) ! Resolution of high-res grid in degrees. + real :: glat(jmn) ! Latitude of each high-res grid row in degrees. real :: lon1,lat1,lon2,lat2 real :: xnsum - print*,"Begin test." + print*,"Begin test of function get_xnsum" -! Set up a 'high-res' 1-degree grid. +! The high-res grid is a global one-degree lat/lon grid. Point (1,1) +! is the south pole/greenwich. do j = 1, jmn glat(j) = -90.0 + float(j-1) * delxn - print*,'j lat ',j,glat(j) enddo +! Bounds of model grid box - straddles equator/greenwich. +! There are approximately 16 high-res points located within +! the model box. + + lon1 = -2.5 + lon2 = 2.5 + lat1 = -1.5 + lat2 = 1.5 + ! First test point. The high-res grid is all ocean. Since all points in ! the model grid box will be sea level, there will be no points ! higher than the average of 0 meters. print*,"Test point 1." - zslm = 0 ! all water - zavg = -999 ! all sea level - -! Bounds of model grid box - straddles greenwich. - - lon1 = -2.5 - lon2 = 2.5 - lat1 = -1.5 - lat2 = 1.5 + zslm = 0 ! high-res grid all water + zavg = -999 ! high-res grid all sea level xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & glat, zavg, zslm, delxn) if (nint(xnsum) /= 0) stop 2 +! Second test point. All high-res points are land with an elevation +! of 50 meters. There will be no points higher than the average. + print*,"Test point 2." - zslm = 1 ! all land + zslm = 1 ! high-res grid all land zavg = 50 ! constant elevation of 50 meters. -! Bounds of model grid box - straddles greenwich. - - lon1 = -2.5 - lon2 = 2.5 - lat1 = -1.5 - lat2 = 1.5 - xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & glat, zavg, zslm, delxn) if (nint(xnsum) /= 0) stop 4 - print*,"Test point 3." +! Third test point. All high-res points are land. One point is +! 100 meters, the rest are 50 meters. Therefore, there should +! be one point higher than the average. - zslm = 1 ! all land - zavg = 50 ! constant elevation of 50 meters. - zavg(359,91) = 100 - -! Bounds of model grid box - straddles greenwich. + print*,"Test point 3." - lon1 = -2.5 - lon2 = 2.5 - lat1 = -1.5 - lat2 = 1.5 + zslm = 1 ! all land + zavg = 50 ! initialize elevation to 50 meters. + zavg(359,91) = 100 ! one point within the model box is 100 meters. xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & glat, zavg, zslm, delxn) if (nint(xnsum) /= 1) stop 6 - print*,"Test point 4." - - zslm = 0 ! all water - zavg = -999 ! set to sea level. - zavg(2,93) = 20 ! this represents an inland lake above sea level. +! Fourth test point. All high-res points are water. One point is +! 20 meters, which represents an inland lake. The other points +! are ocean. Therefore, there should be one point higher than the average. -! Bounds of model grid box - straddles greenwich. + print*,"Test point 4." - lon1 = -2.5 - lon2 = 2.5 - lat1 = -1.5 - lat2 = 1.5 + zslm = 0 ! all water + zavg = -999 ! initialize to sea level. + zavg(2,93) = 20 ! this represents an inland lake above sea level. xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & glat, zavg, zslm, delxn) if (nint(xnsum) /= 1) stop 8 - print*,"Test point 5." - - zslm = 0 ! all water - zavg = -999 ! set to sea level. - zslm(1:2,90:93) = 1 ! half points in grid box land - zavg(1,90:93) = 25 ! land points above sea level - zavg(2,90) = 100 ! land points above sea level - zavg(2,91) = 110 ! land points above sea level - zavg(2,92) = 107 ! land points above sea level - zavg(2,93) = 207 ! land points above sea level +! Fifth test point. Half of the high-res points within the grid +! box are land and half are water. Four high-res points are +! above the average. -! Bounds of model grid box - straddles greenwich. + print*,"Test point 5." - lon1 = -2.5 - lon2 = 2.5 - lat1 = -1.5 - lat2 = 1.5 + zslm = 0 ! Initialize to all water + zavg = -999 ! Initialize to sea level. + zslm(1:2,90:93) = 1 ! Half the points in the grid box are land. + zavg(1,90:93) = 25 ! Set the elevation at the + zavg(2,90) = 100 ! land points so that + zavg(2,91) = 110 ! four points are above the average. + zavg(2,92) = 107 + zavg(2,93) = 207 xnsum = get_xnsum(lon1,lat1,lon2,lat2,imn,jmn, & glat, zavg, zslm, delxn) From 302a0822dfae10e0bdad50eb42358c10eff93430 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 17 Dec 2024 13:52:12 -0600 Subject: [PATCH 24/67] Add unit test for routine get_xnsum2. Fixes #1000. --- tests/orog/CMakeLists.txt | 4 ++ tests/orog/ftst_get_xnsum2.F90 | 72 ++++++++++++++++++++++++++++++++++ 2 files changed, 76 insertions(+) create mode 100644 tests/orog/ftst_get_xnsum2.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index 5cf771cd8..7fb423e13 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -45,3 +45,7 @@ target_link_libraries(ftst_find_nearest_pole_pts orog_lib) add_executable(ftst_get_xnsum ftst_get_xnsum.F90) add_test(NAME orog-ftst_get_xnsum COMMAND ftst_get_xnsum) target_link_libraries(ftst_get_xnsum orog_lib) + +add_executable(ftst_get_xnsum2 ftst_get_xnsum2.F90) +add_test(NAME orog-ftst_get_xnsum2 COMMAND ftst_get_xnsum2) +target_link_libraries(ftst_get_xnsum2 orog_lib) diff --git a/tests/orog/ftst_get_xnsum2.F90 b/tests/orog/ftst_get_xnsum2.F90 new file mode 100644 index 000000000..169db1d36 --- /dev/null +++ b/tests/orog/ftst_get_xnsum2.F90 @@ -0,0 +1,72 @@ + program check_get_xnsum2 + +! Unit test for routine get_xnsum2, which counts the +! number of high-resolution orography points within a +! model grid box, and counts the number of high-resolution +! points higher than a critical height. The critical +! height is a function of the standard deviation of height. + + use orog_utils, only : get_xnsum2 + + implicit none + + integer, parameter :: imn=360 ! i-dimension of high-res grid + integer, parameter :: jmn=181 ! j-dimension of high-res grid + + integer :: j + integer :: zavg(imn,jmn) ! High-res grid terrain + + real :: delxn=360.0/float(imn) ! Resolution of high-res grid in degrees. + real :: glat(jmn) ! Latitude of each high-res grid row in degrees. + real :: lon1,lat1,lon2,lat2,hc + real :: xnsum1,xnsum2 + +! Variables holding the expected test results. + + real, parameter :: EPSILON=0.001 + real :: expected_hc + integer :: expected_xnsum1, expected_xnsum2 + + data expected_hc /677.2/ ! Expected critical height in meters. + data expected_xnsum1 /8/ ! Expected number of high-res pts in model + ! grid box that are above the critical height. + data expected_xnsum2 /16/ ! Expected total number of high-res pts in model grid box. + + print*,"Begin test of routine get_xnsum2." + +! The high-res grid is a global one-degree lat/lon grid. Point (1,1) +! is the south pole/greenwich. + + do j = 1, jmn + glat(j) = -90.0 + float(j-1) * delxn + enddo + +! Bounds of model grid box - straddles equator/greenwich. +! There are 16 high-res points located within the model +! box. The i/j range of these 16 points is i=359,360,1,2 +! and j=0,91,92,93. + + lon1 = -2.5 ! in degrees. + lon2 = 2.5 + lat1 = -1.5 + lat2 = 1.5 + +! Initialize high-res orography. Half of the points +! within the model grid box are at sea level and the +! other half at 1000 meters. + + zavg = -999 ! Flag value for sea level. + zavg(359:360,90:93) = 1000 + + call get_xnsum2(lon1,lat1,lon2,lat2,imn,jmn, & + glat,zavg,delxn,xnsum1,xnsum2,hc) + + if (nint(xnsum1) /= expected_xnsum1) stop 2 + if (nint(xnsum2) /= expected_xnsum2) stop 4 + if (abs(hc-expected_hc) > EPSILON) stop 6 + + print*,"OK" + + print*,"SUCCESS" + + end program check_get_xnsum2 From 79c5d00875c32611479fc1c7fccf9375d614aaec Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 17 Dec 2024 14:50:50 -0600 Subject: [PATCH 25/67] Remove unused variable from ./tests/orog/ftst_transpose.F90 Fixes #1000. --- tests/orog/ftst_transpose.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/orog/ftst_transpose.F90 b/tests/orog/ftst_transpose.F90 index 5d62a7a2e..31f8a8d4d 100644 --- a/tests/orog/ftst_transpose.F90 +++ b/tests/orog/ftst_transpose.F90 @@ -16,7 +16,7 @@ program transpose integer(1) :: mask(imn,jmn) integer(2) :: mask2(imn,jmn) - integer :: i, ii, j, jj + integer :: i, j, jj print*,"Starting test of transpose routines." From d4ba27429cdd636f616fdab78456ebd20be4fde2 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 17 Dec 2024 15:04:48 -0600 Subject: [PATCH 26/67] Baseline unit test for routine get_xnsum3. Fixes #1000. --- tests/orog/CMakeLists.txt | 6 ++- tests/orog/ftst_get_xnsum3.F90 | 74 ++++++++++++++++++++++++++++++++++ 2 files changed, 79 insertions(+), 1 deletion(-) create mode 100644 tests/orog/ftst_get_xnsum3.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index 7fb423e13..1c5eeb150 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -3,7 +3,7 @@ # George Gayno, Ed Hartnett if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel|IntelLLVM)$") - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8 -warn unused") elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffree-line-length-0 -fdefault-real-8") endif() @@ -49,3 +49,7 @@ target_link_libraries(ftst_get_xnsum orog_lib) add_executable(ftst_get_xnsum2 ftst_get_xnsum2.F90) add_test(NAME orog-ftst_get_xnsum2 COMMAND ftst_get_xnsum2) target_link_libraries(ftst_get_xnsum2 orog_lib) + +add_executable(ftst_get_xnsum3 ftst_get_xnsum3.F90) +add_test(NAME orog-ftst_get_xnsum3 COMMAND ftst_get_xnsum3) +target_link_libraries(ftst_get_xnsum3 orog_lib) diff --git a/tests/orog/ftst_get_xnsum3.F90 b/tests/orog/ftst_get_xnsum3.F90 new file mode 100644 index 000000000..45fd76866 --- /dev/null +++ b/tests/orog/ftst_get_xnsum3.F90 @@ -0,0 +1,74 @@ + program check_get_xnsum3 + +! Unit test for routine get_xnsum3, which counts the +! number of high-resolution orography points within a +! model grid box, and counts the number of high-resolution +! points higher than a critical height. The critical +! height is passed into routine get_xnsum3, whereas in +! get_xnsum2 it is computed. + + use orog_utils, only : get_xnsum3 + + implicit none + + integer, parameter :: imn=360 ! i-dimension of high-res grid + integer, parameter :: jmn=181 ! j-dimension of high-res grid + + integer :: j + integer :: zavg(imn,jmn) ! High-res grid terrain + + real :: delxn=360.0/float(imn) ! Resolution of high-res grid in degrees. + real :: glat(jmn) ! Latitude of each high-res grid row in degrees. + real :: lon1,lat1,lon2,lat2,hc + real :: xnsum1,xnsum2 + +! Variables holding the expected test results. + + integer :: expected_xnsum1, expected_xnsum2 + + data expected_xnsum1 /4/ ! Expected number of high-res pts in model + ! grid box that are above the critical height. + data expected_xnsum2 /16/ ! Expected total number of high-res pts in model grid box. + + print*,"Begin test of routine get_xnsum2." + +! The high-res grid is a global one-degree lat/lon grid. Point (1,1) +! is the south pole/greenwich. + + do j = 1, jmn + glat(j) = -90.0 + float(j-1) * delxn + enddo + +! Bounds of model grid box - straddles equator/greenwich. +! There are 16 high-res points located within the model +! box. The i/j range of these 16 points is i=359,360,1,2 +! and j=0,91,92,93. + + lon1 = -2.5 ! in degrees. + lon2 = 2.5 + lat1 = -1.5 + lat2 = 1.5 + +! Initialize high-res orography. Half of the points +! within the model grid box are at sea level, one quarter +! are at 500 meters and one quarter are at 1000 meters. +! The critical height was chosen as 700 meters. So, +! four points should be above the critical value. + + zavg = -999 ! Flag value for sea level. + zavg(359,90:93) = 1000 + zavg(360,90:93) = 500 + + hc = 700.0 ! Critical height in meters. + + call get_xnsum3(lon1,lat1,lon2,lat2,imn,jmn, & + glat,zavg,delxn,xnsum1,xnsum2,hc) + + if (nint(xnsum1) /= expected_xnsum1) stop 2 + if (nint(xnsum2) /= expected_xnsum2) stop 4 + + print*,"OK" + + print*,"SUCCESS" + + end program check_get_xnsum3 From ed80ff3b6b0b736a0f734879c2d61a36cd87022d Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 30 Dec 2024 15:33:26 -0600 Subject: [PATCH 27/67] Create initial unit test for routine remove_isolated_pts. Fixes #1000. --- tests/orog/CMakeLists.txt | 4 ++ tests/orog/ftst_rm_isolated_pts.F90 | 69 +++++++++++++++++++++++++++++ 2 files changed, 73 insertions(+) create mode 100644 tests/orog/ftst_rm_isolated_pts.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index 1c5eeb150..b42f794af 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -53,3 +53,7 @@ target_link_libraries(ftst_get_xnsum2 orog_lib) add_executable(ftst_get_xnsum3 ftst_get_xnsum3.F90) add_test(NAME orog-ftst_get_xnsum3 COMMAND ftst_get_xnsum3) target_link_libraries(ftst_get_xnsum3 orog_lib) + +add_executable(ftst_rm_isolated_pts ftst_rm_isolated_pts.F90) +add_test(NAME orog-ftst_rm_isolated_pts COMMAND ftst_rm_isolated_pts) +target_link_libraries(ftst_rm_isolated_pts orog_lib) diff --git a/tests/orog/ftst_rm_isolated_pts.F90 b/tests/orog/ftst_rm_isolated_pts.F90 new file mode 100644 index 000000000..902af6e3a --- /dev/null +++ b/tests/orog/ftst_rm_isolated_pts.F90 @@ -0,0 +1,69 @@ + program rm_isolated_pts + +! Unit test for subroutine remove_isolated_pts. +! +! Author George Gayno NCEP/EMC + + use orog_utils, only : remove_isolated_pts + + implicit none + + integer :: im, jm, i, j, k + + real, allocatable :: slm(:,:), oro(:,:), var(:,:), & + var4(:,:), oa(:,:,:), ol(:,:,:) + + print*,"Starting test of remove_isolated_pts." + + im = 3 + jm = 3 + + allocate (slm(im,jm)) + allocate (oro(im,jm)) + allocate (var(im,jm)) + allocate (var4(im,jm)) + allocate (oa(im,jm,4)) + allocate (ol(im,jm,4)) + +! Initialize grid to all water. + + slm = 0.0 + oro = 0.0 + var = 0.0 + var4 = 0.0 + oa = 0.0 + ol = 0.0 + +! This is an isolated island. The island should be +! removed (slm set to 0.0) and all other fields +! should be the average of the surrounding points (which +! for this test case is zero. + + slm(2,2) = 1. + oro(2,2) = 50. + var(2,2) = 10. + var4(2,2) = 5. + oa(2,2,:) = -1. + ol(2,2,:) = 1. + + call remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) + + do j = 1, jm + do i = 1, im + if (slm(i,j) /= 0.0) stop 2 + if (oro(i,j) /= 0.0) stop 4 + if (var(i,j) /= 0.0) stop 6 + if (var4(i,j) /= 0.0) stop 8 + do k = 1, 4 + if (oa(i,j,k) /= 0.0) stop 10 + if (ol(i,j,k) /= 0.0) stop 12 + enddo + enddo + enddo + + deallocate (slm, oro, var, var4, oa, ol) + + print*,"OK" + print*,"SUCCSSS" + + end program rm_isolated_pts From 6c03013df4151ce4d69115b44653f72f0a68b9f3 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 31 Dec 2024 07:28:29 -0600 Subject: [PATCH 28/67] Add more realistic values for test point. Fixes #1000. --- tests/orog/ftst_rm_isolated_pts.F90 | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/tests/orog/ftst_rm_isolated_pts.F90 b/tests/orog/ftst_rm_isolated_pts.F90 index 902af6e3a..7018cc76f 100644 --- a/tests/orog/ftst_rm_isolated_pts.F90 +++ b/tests/orog/ftst_rm_isolated_pts.F90 @@ -25,7 +25,7 @@ program rm_isolated_pts allocate (oa(im,jm,4)) allocate (ol(im,jm,4)) -! Initialize grid to all water. +! Initialize grid to all ocean. slm = 0.0 oro = 0.0 @@ -39,12 +39,20 @@ program rm_isolated_pts ! should be the average of the surrounding points (which ! for this test case is zero. - slm(2,2) = 1. - oro(2,2) = 50. - var(2,2) = 10. - var4(2,2) = 5. - oa(2,2,:) = -1. - ol(2,2,:) = 1. + slm(2,2) = 1.0 + oro(2,2) = 50.0 + var(2,2) = 10.0 + var4(2,2) = 5.0 + + oa(2,2,1) = -1.0 + oa(2,2,2) = -0.5 + oa(2,2,3) = 0.5 + oa(2,2,4) = 1.0 + + ol(2,2,1) = 0.1 + ol(2,2,2) = 0.25 + ol(2,2,3) = 0.5 + ol(2,2,4) = 1.0 call remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) From 276a2ba3091e0362923d2f14125338a7c0cb0e67 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 31 Dec 2024 13:13:32 -0600 Subject: [PATCH 29/67] Add another test point to the rm_isolated_pts unit test. Fixes #1000. --- tests/orog/ftst_rm_isolated_pts.F90 | 134 +++++++++++++++++++++++++--- 1 file changed, 124 insertions(+), 10 deletions(-) diff --git a/tests/orog/ftst_rm_isolated_pts.F90 b/tests/orog/ftst_rm_isolated_pts.F90 index 7018cc76f..c19af0f90 100644 --- a/tests/orog/ftst_rm_isolated_pts.F90 +++ b/tests/orog/ftst_rm_isolated_pts.F90 @@ -10,9 +10,14 @@ program rm_isolated_pts integer :: im, jm, i, j, k + real, parameter :: EPSILON=0.001 + real, allocatable :: slm(:,:), oro(:,:), var(:,:), & var4(:,:), oa(:,:,:), ol(:,:,:) + real, allocatable :: slm_expected(:,:), oro_expected(:,:), var_expected(:,:), & + var4_expected(:,:), oa_expected(:,:,:), ol_expected(:,:,:) + print*,"Starting test of remove_isolated_pts." im = 3 @@ -25,8 +30,18 @@ program rm_isolated_pts allocate (oa(im,jm,4)) allocate (ol(im,jm,4)) +! Test Point 1. + +! This is an isolated island. The island should be +! removed (slm set to 0.0) and all other fields +! should be the average of the surrounding points (which +! for this test case is zero). All surrounding points +! should be unchanged. + ! Initialize grid to all ocean. + print*,'-Test point 1.' + slm = 0.0 oro = 0.0 var = 0.0 @@ -34,10 +49,7 @@ program rm_isolated_pts oa = 0.0 ol = 0.0 -! This is an isolated island. The island should be -! removed (slm set to 0.0) and all other fields -! should be the average of the surrounding points (which -! for this test case is zero. +! Initialize an island in the middle of the grid. slm(2,2) = 1.0 oro(2,2) = 50.0 @@ -54,22 +66,124 @@ program rm_isolated_pts ol(2,2,3) = 0.5 ol(2,2,4) = 1.0 + allocate (slm_expected(im,jm)) + allocate (oro_expected(im,jm)) + allocate (var_expected(im,jm)) + allocate (var4_expected(im,jm)) + allocate (oa_expected(im,jm,4)) + allocate (ol_expected(im,jm,4)) + +! All points should be ocean (all fields zero). + + slm_expected = 0.0 + oro_expected = 0.0 + var_expected = 0.0 + var4_expected = 0.0 + oa_expected = 0.0 + ol_expected = 0.0 + + call remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) + + do j = 1, jm + do i = 1, im + if (abs(slm(i,j)-slm_expected(i,j)) > EPSILON) stop 2 + if (abs(oro(i,j)-oro_expected(i,j)) > EPSILON) stop 4 + if (abs(var(i,j)-var_expected(i,j)) > EPSILON) stop 6 + if (abs(var4(i,j)-var4_expected(i,j)) > EPSILON) stop 8 + do k = 1, 4 + if (abs(oa(i,j,k)-oa_expected(i,j,k)) > EPSILON) stop 10 + if (abs(ol(i,j,k)-ol_expected(i,j,k)) > EPSILON) stop 12 + enddo + enddo + enddo + +! Test Point 2 + +! Now remove an isolated water point. + + print*,'-Test point 2.' + + slm = 1.0 + slm(2,2) = 0.0 ! water point + + oro(1,1) = 5.0 + oro(2,1) = 10.0 + oro(3,1) = 17.0 + oro(1,2) = 22.0 + oro(2,2) = 0.0 ! water point. + oro(3,2) = 100.0 + oro(1,3) = 34.0 + oro(2,3) = 55.0 + oro(3,3) = 68.0 + + var = 0.5 * oro + var4 = 0.25 * oro + + ol(1,1,1) = 0.1 + ol(2,1,1) = 0.2 + ol(3,1,1) = 0.3 + ol(1,2,1) = 0.4 + ol(2,2,1) = 0.0 ! water point + ol(3,2,1) = 0.6 + ol(1,3,1) = 0.7 + ol(2,3,1) = 0.8 + ol(3,3,1) = 0.9 + + do j = 1, jm + do i = 1, im + ol(i,j,2) = ol(i,j,1) * .75 + ol(i,j,3) = ol(i,j,1) * .5 + ol(i,j,4) = ol(i,j,1) * .25 + enddo + enddo + + oa = -(ol) + +! All points should be land. The former isolated ocean +! point should have values that are the average of the +! surrounding land points. All other points should remain +! unchanged. + + slm_expected = slm + oro_expected = oro + var_expected = var + var4_expected = var4 + ol_expected = ol + oa_expected = oa + +! Former ocean point should be average of the surrounding +! ocean points. + + slm_expected(2,2) = 1.0 ! land point + oro_expected(2,2) = 38.875 + var_expected(2,2) = 0.5 * oro_expected(2,2) + var4_expected(2,2) = 0.25 * oro_expected(2,2) + ol_expected(2,2,1) = 0.5 + ol_expected(2,2,2) = 0.375 + ol_expected(2,2,3) = 0.25 + ol_expected(2,2,4) = 0.125 + oa_expected(2,2,1) = -0.5 + oa_expected(2,2,2) = -0.375 + oa_expected(2,2,3) = -0.25 + oa_expected(2,2,4) = -0.125 + call remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) do j = 1, jm do i = 1, im - if (slm(i,j) /= 0.0) stop 2 - if (oro(i,j) /= 0.0) stop 4 - if (var(i,j) /= 0.0) stop 6 - if (var4(i,j) /= 0.0) stop 8 + if (abs(slm(i,j)-slm_expected(i,j)) > EPSILON) stop 22 + if (abs(oro(i,j)-oro_expected(i,j)) > EPSILON) stop 24 + if (abs(var(i,j)-var_expected(i,j)) > EPSILON) stop 26 + if (abs(var4(i,j)-var4_expected(i,j)) > EPSILON) stop 28 do k = 1, 4 - if (oa(i,j,k) /= 0.0) stop 10 - if (ol(i,j,k) /= 0.0) stop 12 + if (abs(oa(i,j,k)-oa_expected(i,j,k)) > EPSILON) stop 30 + if (abs(ol(i,j,k)-ol_expected(i,j,k)) > EPSILON) stop 32 enddo enddo enddo deallocate (slm, oro, var, var4, oa, ol) + deallocate (slm_expected, oro_expected, var_expected, var4_expected, oa_expected, ol_expected) print*,"OK" print*,"SUCCSSS" From 25c6b54109474087769ec4b9b0e686314ee9c382 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 31 Dec 2024 14:24:47 -0600 Subject: [PATCH 30/67] Add some diagnostic print to remove_isolated_pts. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 index 363dfe95f..7dfb45423 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 @@ -587,9 +587,12 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) JN=J-1 JS=J+1 i_loop : DO I=1,IM + print*,'check point ',i,j IW=MOD(I+IM-2,IM)+1 IE=MOD(I,IM)+1 SLMA=SLM(IW,J)+SLM(IE,J) + print*,'check points to the w ',iw,j + print*,'check points to the e ',ie,j OROA=ORO(IW,J)+ORO(IE,J) VARA=VAR(IW,J)+VAR(IE,J) VAR4A=VAR4(IW,J)+VAR4(IE,J) @@ -605,6 +608,9 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) INW=MOD(IN+IM-2,IM)+1 INE=MOD(IN,IM)+1 SLMA=SLMA+SLM(INW,JN)+SLM(IN,JN)+SLM(INE,JN) + print*,'check points to the nw ',inw,jn + print*,'check points to the n ',in,jn + print*,'check points to the ne ',ine,jn OROA=OROA+ORO(INW,JN)+ORO(IN,JN)+ORO(INE,JN) VARA=VARA+VAR(INW,JN)+VAR(IN,JN)+VAR(INE,JN) VAR4A=VAR4A+VAR4(INW,JN)+VAR4(IN,JN)+VAR4(INE,JN) @@ -614,6 +620,7 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) ENDDO WGTA=WGTA+3 ELSE + print*,'got here 2' INW=INT(XN) INE=MOD(INW,IM)+1 SLMA=SLMA+SLM(INW,JN)+SLM(INE,JN) @@ -628,9 +635,13 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) ENDIF XS=(I-1)+1 IF(ABS(XS-NINT(XS)).LT.1.E-2) THEN +! print*,'got here 3' IS=MOD(NINT(XS)-1,IM)+1 ISW=MOD(IS+IM-2,IM)+1 ISE=MOD(IS,IM)+1 + print*,'check points to the sw ',isw,js + print*,'check points to the s ',is,js + print*,'check points to the se ',ise,js SLMA=SLMA+SLM(ISW,JS)+SLM(IS,JS)+SLM(ISE,JS) OROA=OROA+ORO(ISW,JS)+ORO(IS,JS)+ORO(ISE,JS) VARA=VARA+VAR(ISW,JS)+VAR(IS,JS)+VAR(ISE,JS) @@ -640,7 +651,9 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) OLA(K)=OLA(K)+OL(ISW,JS,K)+OL(IS,JS,K)+OL(ISE,JS,K) ENDDO WGTA=WGTA+3 +! print*,'slm/wgta ',SLM(ISW,JS),SLM(IS,JS),SLM(ISE,JS),slma,wgta ELSE + print*,'got here 4' ISW=INT(XS) ISE=MOD(ISW,IM)+1 SLMA=SLMA+SLM(ISW,JS)+SLM(ISE,JS) @@ -660,6 +673,7 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) OAA(K)=OAA(K)/WGTA OLA(K)=OLA(K)/WGTA ENDDO +! print*,'before if - i/j/slm/slma/wgta ',i,j,SLM(I,J),slma,wgta IF(SLM(I,J).EQ.0..AND.SLMA.EQ.WGTA) THEN PRINT '(" - SEA ",2F8.0," MODIFIED TO LAND",2F8.0, & " AT ",2I8)',ORO(I,J),VAR(I,J),OROA,VARA,I,J From a7545ea68da795ea6887c9a9bfac9c038085a9aa Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 31 Dec 2024 15:16:29 -0600 Subject: [PATCH 31/67] Update routine remove_isolated_pts to not process the e/w edges of the tile. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 index 7dfb45423..88b6733a0 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 @@ -586,7 +586,7 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) iso_loop : DO J=2,JM-1 JN=J-1 JS=J+1 - i_loop : DO I=1,IM + i_loop : DO I=2,IM-1 print*,'check point ',i,j IW=MOD(I+IM-2,IM)+1 IE=MOD(I,IM)+1 From 11dda0c40064190e67b6f4ffc46feb6c20e77609 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 31 Dec 2024 15:27:50 -0600 Subject: [PATCH 32/67] Update first unit test point to be a 4x4 grid as to test both odd and even dimensioned grids. Fixes #1000. --- tests/orog/ftst_rm_isolated_pts.F90 | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/tests/orog/ftst_rm_isolated_pts.F90 b/tests/orog/ftst_rm_isolated_pts.F90 index c19af0f90..b82265ccd 100644 --- a/tests/orog/ftst_rm_isolated_pts.F90 +++ b/tests/orog/ftst_rm_isolated_pts.F90 @@ -20,8 +20,8 @@ program rm_isolated_pts print*,"Starting test of remove_isolated_pts." - im = 3 - jm = 3 + im = 4 + jm = 4 allocate (slm(im,jm)) allocate (oro(im,jm)) @@ -97,12 +97,25 @@ program rm_isolated_pts enddo enddo + deallocate (slm, oro, var, var4, oa, ol) + deallocate (slm_expected, oro_expected, var_expected, var4_expected, oa_expected, ol_expected) + ! Test Point 2 ! Now remove an isolated water point. print*,'-Test point 2.' + im = 3 + jm = 3 + + allocate (slm(im,jm)) + allocate (oro(im,jm)) + allocate (var(im,jm)) + allocate (var4(im,jm)) + allocate (oa(im,jm,4)) + allocate (ol(im,jm,4)) + slm = 1.0 slm(2,2) = 0.0 ! water point @@ -144,6 +157,13 @@ program rm_isolated_pts ! surrounding land points. All other points should remain ! unchanged. + allocate (slm_expected(im,jm)) + allocate (oro_expected(im,jm)) + allocate (var_expected(im,jm)) + allocate (var4_expected(im,jm)) + allocate (oa_expected(im,jm,4)) + allocate (ol_expected(im,jm,4)) + slm_expected = slm oro_expected = oro var_expected = var From 380a57cc438f10769b7aee15f3cb7fe7dfe7ba86 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Thu, 2 Jan 2025 12:16:18 -0600 Subject: [PATCH 33/67] Remove unused 'if' branch from 'remove_isolated_pts'. I suspect this branch was used for the old GFS reduced grid. Fixe #1000. --- .../orog_mask_tools.fd/orog.fd/orog_utils.F90 | 104 ++++++------------ 1 file changed, 33 insertions(+), 71 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 index 88b6733a0..b0e14e9d4 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/orog_utils.F90 @@ -574,7 +574,7 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) integer :: iw, ie, wgta, is, ise integer :: in, ine, inw, isw - real :: slma, oroa, vara, var4a, xn, xs + real :: slma, oroa, vara, var4a real, allocatable :: oaa(:), ola(:) ! REMOVE ISOLATED POINTS @@ -587,12 +587,10 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) JN=J-1 JS=J+1 i_loop : DO I=2,IM-1 - print*,'check point ',i,j - IW=MOD(I+IM-2,IM)+1 - IE=MOD(I,IM)+1 +! Check the points to the 'west' and 'east'. + IW=I-1 + IE=I+1 SLMA=SLM(IW,J)+SLM(IE,J) - print*,'check points to the w ',iw,j - print*,'check points to the e ',ie,j OROA=ORO(IW,J)+ORO(IE,J) VARA=VAR(IW,J)+VAR(IE,J) VAR4A=VAR4(IW,J)+VAR4(IE,J) @@ -602,70 +600,33 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) OLA(K)=OL(IW,J,K)+OL(IE,J,K) ENDDO WGTA=2 - XN=(I-1)+1 - IF(ABS(XN-NINT(XN)).LT.1.E-2) THEN - IN=MOD(NINT(XN)-1,IM)+1 - INW=MOD(IN+IM-2,IM)+1 - INE=MOD(IN,IM)+1 - SLMA=SLMA+SLM(INW,JN)+SLM(IN,JN)+SLM(INE,JN) - print*,'check points to the nw ',inw,jn - print*,'check points to the n ',in,jn - print*,'check points to the ne ',ine,jn - OROA=OROA+ORO(INW,JN)+ORO(IN,JN)+ORO(INE,JN) - VARA=VARA+VAR(INW,JN)+VAR(IN,JN)+VAR(INE,JN) - VAR4A=VAR4A+VAR4(INW,JN)+VAR4(IN,JN)+VAR4(INE,JN) - DO K=1,4 - OAA(K)=OAA(K)+OA(INW,JN,K)+OA(IN,JN,K)+OA(INE,JN,K) - OLA(K)=OLA(K)+OL(INW,JN,K)+OL(IN,JN,K)+OL(INE,JN,K) - ENDDO - WGTA=WGTA+3 - ELSE - print*,'got here 2' - INW=INT(XN) - INE=MOD(INW,IM)+1 - SLMA=SLMA+SLM(INW,JN)+SLM(INE,JN) - OROA=OROA+ORO(INW,JN)+ORO(INE,JN) - VARA=VARA+VAR(INW,JN)+VAR(INE,JN) - VAR4A=VAR4A+VAR4(INW,JN)+VAR4(INE,JN) - DO K=1,4 - OAA(K)=OAA(K)+OA(INW,JN,K)+OA(INE,JN,K) - OLA(K)=OLA(K)+OL(INW,JN,K)+OL(INE,JN,K) - ENDDO - WGTA=WGTA+2 - ENDIF - XS=(I-1)+1 - IF(ABS(XS-NINT(XS)).LT.1.E-2) THEN -! print*,'got here 3' - IS=MOD(NINT(XS)-1,IM)+1 - ISW=MOD(IS+IM-2,IM)+1 - ISE=MOD(IS,IM)+1 - print*,'check points to the sw ',isw,js - print*,'check points to the s ',is,js - print*,'check points to the se ',ise,js - SLMA=SLMA+SLM(ISW,JS)+SLM(IS,JS)+SLM(ISE,JS) - OROA=OROA+ORO(ISW,JS)+ORO(IS,JS)+ORO(ISE,JS) - VARA=VARA+VAR(ISW,JS)+VAR(IS,JS)+VAR(ISE,JS) - VAR4A=VAR4A+VAR4(ISW,JS)+VAR4(IS,JS)+VAR4(ISE,JS) - DO K=1,4 - OAA(K)=OAA(K)+OA(ISW,JS,K)+OA(IS,JS,K)+OA(ISE,JS,K) - OLA(K)=OLA(K)+OL(ISW,JS,K)+OL(IS,JS,K)+OL(ISE,JS,K) - ENDDO - WGTA=WGTA+3 -! print*,'slm/wgta ',SLM(ISW,JS),SLM(IS,JS),SLM(ISE,JS),slma,wgta - ELSE - print*,'got here 4' - ISW=INT(XS) - ISE=MOD(ISW,IM)+1 - SLMA=SLMA+SLM(ISW,JS)+SLM(ISE,JS) - OROA=OROA+ORO(ISW,JS)+ORO(ISE,JS) - VARA=VARA+VAR(ISW,JS)+VAR(ISE,JS) - VAR4A=VAR4A+VAR4(ISW,JS)+VAR4(ISE,JS) - DO K=1,4 - OAA(K)=OAA(K)+OA(ISW,JS,K)+OA(ISE,JS,K) - OLA(K)=OLA(K)+OL(ISW,JS,K)+OL(ISE,JS,K) - ENDDO - WGTA=WGTA+2 - ENDIF +! Check the points to the 'northwest', 'north' and 'northeast' + IN=I + INW=I-1 + INE=I+1 + SLMA=SLMA+SLM(INW,JN)+SLM(IN,JN)+SLM(INE,JN) + OROA=OROA+ORO(INW,JN)+ORO(IN,JN)+ORO(INE,JN) + VARA=VARA+VAR(INW,JN)+VAR(IN,JN)+VAR(INE,JN) + VAR4A=VAR4A+VAR4(INW,JN)+VAR4(IN,JN)+VAR4(INE,JN) + DO K=1,4 + OAA(K)=OAA(K)+OA(INW,JN,K)+OA(IN,JN,K)+OA(INE,JN,K) + OLA(K)=OLA(K)+OL(INW,JN,K)+OL(IN,JN,K)+OL(INE,JN,K) + ENDDO + WGTA=WGTA+3 +! Check the points to the 'southwest', 'south' and 'southeast' + IS=I + ISW=I-1 + ISE=I+1 + SLMA=SLMA+SLM(ISW,JS)+SLM(IS,JS)+SLM(ISE,JS) + OROA=OROA+ORO(ISW,JS)+ORO(IS,JS)+ORO(ISE,JS) + VARA=VARA+VAR(ISW,JS)+VAR(IS,JS)+VAR(ISE,JS) + VAR4A=VAR4A+VAR4(ISW,JS)+VAR4(IS,JS)+VAR4(ISE,JS) + DO K=1,4 + OAA(K)=OAA(K)+OA(ISW,JS,K)+OA(IS,JS,K)+OA(ISE,JS,K) + OLA(K)=OLA(K)+OL(ISW,JS,K)+OL(IS,JS,K)+OL(ISE,JS,K) + ENDDO + WGTA=WGTA+3 +! Take the average of the surrounding the points. OROA=OROA/WGTA VARA=VARA/WGTA VAR4A=VAR4A/WGTA @@ -673,7 +634,7 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) OAA(K)=OAA(K)/WGTA OLA(K)=OLA(K)/WGTA ENDDO -! print*,'before if - i/j/slm/slma/wgta ',i,j,SLM(I,J),slma,wgta +! Isolated water point. IF(SLM(I,J).EQ.0..AND.SLMA.EQ.WGTA) THEN PRINT '(" - SEA ",2F8.0," MODIFIED TO LAND",2F8.0, & " AT ",2I8)',ORO(I,J),VAR(I,J),OROA,VARA,I,J @@ -685,6 +646,7 @@ subroutine remove_isolated_pts(im,jm,slm,oro,var,var4,oa,ol) OA(I,J,K)=OAA(K) OL(I,J,K)=OLA(K) ENDDO +! Isolated land point. ELSEIF(SLM(I,J).EQ.1..AND.SLMA.EQ.0.) THEN PRINT '(" - LAND",2F8.0," MODIFIED TO SEA ",2F8.0, & " AT ",2I8)',ORO(I,J),VAR(I,J),OROA,VARA,I,J From 0c9f20b7ca7c6b36408ef8c60f1c995d4d25fcfb Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 3 Jan 2025 07:56:06 -0600 Subject: [PATCH 34/67] Some final updates to the ftst_rm_isolated_pts unit test. Fixes #1000. --- tests/orog/ftst_rm_isolated_pts.F90 | 34 +++++++++++++++++------------ 1 file changed, 20 insertions(+), 14 deletions(-) diff --git a/tests/orog/ftst_rm_isolated_pts.F90 b/tests/orog/ftst_rm_isolated_pts.F90 index b82265ccd..21cbdf42e 100644 --- a/tests/orog/ftst_rm_isolated_pts.F90 +++ b/tests/orog/ftst_rm_isolated_pts.F90 @@ -18,9 +18,11 @@ program rm_isolated_pts real, allocatable :: slm_expected(:,:), oro_expected(:,:), var_expected(:,:), & var4_expected(:,:), oa_expected(:,:,:), ol_expected(:,:,:) - print*,"Starting test of remove_isolated_pts." + print*,"- Starting test of remove_isolated_pts." - im = 4 +! A 5x4 grid. + + im = 5 jm = 4 allocate (slm(im,jm)) @@ -38,16 +40,16 @@ program rm_isolated_pts ! for this test case is zero). All surrounding points ! should be unchanged. -! Initialize grid to all ocean. + print*,'- Test point 1.' - print*,'-Test point 1.' +! Initialize grid to all ocean. - slm = 0.0 - oro = 0.0 - var = 0.0 - var4 = 0.0 - oa = 0.0 - ol = 0.0 + slm = 0.0 ! land/sea mask + oro = 0.0 ! orography + var = 0.0 ! standard deviation of orography + var4 = 0.0 ! convexity + oa = 0.0 ! orographic asymetry + ol = 0.0 ! orographic length scale ! Initialize an island in the middle of the grid. @@ -102,9 +104,13 @@ program rm_isolated_pts ! Test Point 2 -! Now remove an isolated water point. +! Now remove an isolated water point. The point should be changed +! to land (slm=1.0) and all other fields should be the average +! of the eight surrounding points. + + print*,'- Test point 2.' - print*,'-Test point 2.' +! A 3x3 grid. im = 3 jm = 3 @@ -116,8 +122,8 @@ program rm_isolated_pts allocate (oa(im,jm,4)) allocate (ol(im,jm,4)) - slm = 1.0 - slm(2,2) = 0.0 ! water point + slm = 1.0 ! Initialize grid to all land. + slm(2,2) = 0.0 ! Set interior point to water. oro(1,1) = 5.0 oro(2,1) = 10.0 From d63ecfd800534002bf186e1df5e389cfda7f680c Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 3 Jan 2025 09:13:34 -0600 Subject: [PATCH 35/67] Remove unused line from ./tests/orog/CMakeLists.txt Fixes #1000. --- tests/orog/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index b42f794af..6b709dce8 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -8,7 +8,7 @@ elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffree-line-length-0 -fdefault-real-8") endif() -include_directories(${PROJECT_SOURCE_DIR}) +#include_directories(${PROJECT_SOURCE_DIR}) add_executable(ftst_ll2xyz ftst_ll2xyz.F90) add_test(NAME orog-ftst_ll2xyz COMMAND ftst_ll2xyz) From be93997e03ebfb9c9e90d01b29465da32f362d0c Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 3 Jan 2025 12:13:06 -0600 Subject: [PATCH 36/67] Baseline a unit test for routine "read_global_mask" and a low-resolution version of the umd land use file. Fixes #1000. --- tests/orog/CMakeLists.txt | 11 +++++++- tests/orog/data/landcover.umd.lowres.nc | Bin 0 -> 13774 bytes tests/orog/ftst_read_global_mask.F90 | 36 ++++++++++++++++++++++++ 3 files changed, 46 insertions(+), 1 deletion(-) create mode 100644 tests/orog/data/landcover.umd.lowres.nc create mode 100644 tests/orog/ftst_read_global_mask.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index 6b709dce8..4fe257e9e 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -8,7 +8,12 @@ elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffree-line-length-0 -fdefault-real-8") endif() -#include_directories(${PROJECT_SOURCE_DIR}) +# Copy necessary test files from the source data directory to the +# build data directory. + +# Note, the "read_global_mask" routine expects the filename "landcover.umd.30s.nc". +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/landcover.umd.lowres.nc ${CMAKE_CURRENT_BINARY_DIR}/landcover.umd.30s.nc) add_executable(ftst_ll2xyz ftst_ll2xyz.F90) add_test(NAME orog-ftst_ll2xyz COMMAND ftst_ll2xyz) @@ -57,3 +62,7 @@ target_link_libraries(ftst_get_xnsum3 orog_lib) add_executable(ftst_rm_isolated_pts ftst_rm_isolated_pts.F90) add_test(NAME orog-ftst_rm_isolated_pts COMMAND ftst_rm_isolated_pts) target_link_libraries(ftst_rm_isolated_pts orog_lib) + +add_executable(ftst_read_global_mask ftst_read_global_mask.F90) +add_test(NAME orog-ftst_read_global_mask COMMAND ftst_read_global_mask) +target_link_libraries(ftst_read_global_mask orog_lib) diff --git a/tests/orog/data/landcover.umd.lowres.nc b/tests/orog/data/landcover.umd.lowres.nc new file mode 100644 index 0000000000000000000000000000000000000000..5652f438fd7188a41212d79d5ce55ec0b6fdf549 GIT binary patch literal 13774 zcmeGiZEPGzb=G@VJI*Fe(=={_qioU&Il(#KozEXhaO>ngV#gQT^*IPi%Jp)$K3{X* zwRZR77*R?Al>`VV4T?~!f-I4?sah$ZN>ypZ$B!To2_oX-2c$wRKOzy7!XF@1YI$#F z-udct!r6)~B=*^NZ+G_1o0<1H_jZo=#S*I;H#F8YG%$h>B>bvFv;^gYjsAZ+Lw&K4 zx>YZzbz69&xTZEwl(?~uDcG`R4Srwu;WuM@1!CG(yil12fAA=6X4}(_ zYZtOETLbeE)H|saWp3@)C6(-V?q1MhI zej>VSIBrM4%<_fvY4X_@tts1ZpuMq*WpKwlEG55o&M$)eh;fY!a8BW1Fxc8242F2; zcH)BUYUEO85#%hw!~%KMz@i|96YdUmceTRz4zVC&gZ$=U(>+HZ$Od~$9^pB%Z4aVg z`2+zoZJIGZsMRObIv5?|Dt$~EfoZtWSVA9~H0%uM+tI<29IzCMecNgfW@P z8+wl}Ez!(c0cR`#9b5INm8*D?FSz5SLxHR{(FR2L0f)K+29of?)VWuL$1+7(hmCcL zqNws|8Lva{x@+^Mvcoksu%;)0%Fn%5Z4}-+(;_-PBm3Xe%sZ%lMiFB!);@-$Dolx6 zcr>{?mEyv{Xn^(|>)#(AO!W^9rc=?~$+!#?9e(V_*tkFp*`V95mSJMhaI?9wNUShX z$W;Xj-eu)-1vD#m>ip~%+1n@Y;d|BNmb4mM<%T+xS0u5McH?4u(YyO23lGy*z7J6rpi4KG)#j&w^{ ztxOERk{^$NRzTDkgO;%kU&fu+}4Cp*>p zm4<4q$*d>c3ml?GtqpAx>zb=I-%U+-kq-AgNVW)z(i&PjxMmgnY4gZ8su~y5u=I+H zu=r>9{k&Rn^>Vs|YgJsAj2%m^xJMhs3>T$19{=Z2+&Aw#e_}Cus_L`!itA#>bASI) zS#b~G6zU+yQEDs4NjbOjy@n|d_n0E8H!y!aUaH&h{wa|G?$t!EZ^_P-Ic|~f;z+8p zb&#ka-KJ%tidvSPE0`JBLVY-0K{Dpv1DyyX1Du0)zsN5#b(YM6>+l4Fu8PE2V3)`J zqD+?1)aS1cbD8fDedRV6wwn30nVAq<xqOnf^b>scu>*1^ z6Z-Ev8T}zh$_TwIvXVtU0Bi8C7U~u4`0P7-8`dILm4N)6BNxg@&$kG(yLr!r=y^VI z2m#+%h0GW;1!7uK2*kL{40V& zDPfr#7pTCoIW>l5Q6TzV%9g@Sk>!*>Y9OqBsm@V))A|#khOr@S-&$ z)kzgAs1(h@xlfP2m>B6B=8bhch{}ibvlV_-aG;+VJc?lpdVB>?(W|0q-b--ca!li|A-%%vnr013D8FptLe=P_9E6ono|oD$ zM^a%eW!Az9%I+$Rt-M(i_B8d(wJu`r;boBkrYhPdH_KGbcg--mb`WEU3K6N+eJWWB zVT?btPsQuXd!+f-nb)`OjSeMwoZ>GKC}xShFOaZzq`eVOvYF?lrI9&OJJFL}QhQ3; z9vLXL6CbivYS$<{<0JqclUqS?=Y4w;@Fd_#z>|O{0Z#&+1Uw0N67VG8Nx+kUCjn0a z3z0x&1%m3d^1UoC1*`1BR>AM9#P$<(ar^Q8AHC+xHzRs=h-GM%?UJ?e<(Sd$fBT!~ z_zQ*4Q8gxrZw*J{iTFr-Fd9$uw}zci9>D?%9$<6~clPyR`BV~Xz~`bWqVr>7R8b*( zL7p`rslGyLkPWA=KCx>)*~r2T(~9(7^LEXU*M(qo!udIe8P5Ikp-*zzBx&SYD^urLD@3vWFTA?EQm2`;kT|8(>;CcwyHe_>Wt9$I?Q*!iZgN9_((xM zJS|F9L&ZnlQ!fwdQk_xpQDw4To>HZXqT-_=R^Y*k^hCxvOunc1s6MH}Yk2ZN#Yff2 zMrcio^Bm+1ijS)5DxB~j^HF?M$ZUjG*=4R=am^lh67VG8Nx+kUCxQPz35c_G{2R~N z9X|%J?%`f82*5|D0PcRc_tHfGxQ`(KKK~{Fw6`Mwe(7ZZdqA{11HDLo(H-c$@)kg8 vp!cuW0gem^zDEao^(_EEAHn`W&-nuV=m7wr&y4{9Js*E+pckig4Ep^CnyjEv literal 0 HcmV?d00001 diff --git a/tests/orog/ftst_read_global_mask.F90 b/tests/orog/ftst_read_global_mask.F90 new file mode 100644 index 000000000..9a4a67229 --- /dev/null +++ b/tests/orog/ftst_read_global_mask.F90 @@ -0,0 +1,36 @@ + program read_gbl_mask + +! Test routine "read_global_mask" using a reduced-size +! version (6 x 3 vs 43200 x 21600) of the umd land mask file. + + use io_utils, only : read_global_mask + + implicit none + + integer, parameter :: im=6 + integer, parameter :: jm=3 + + integer :: i, j + + integer(kind=1) :: mask(im,jm) + integer(kind=1) :: mask_expected(im,jm) + +! Note the routine flips the i and j directions. + data mask_expected /1,1,1,0,0,1, & + 1,1,1,0,0,0, & + 1,1,1,0,0,0/ + + print*,"- Begin test of read_global_mask" + + call read_global_mask(im, jm, mask) + + do j = 1, jm + do i = 1, im + if (mask(i,j) /= mask_expected(i,j)) stop 3 + enddo + enddo + + print*,"OK" + print*,"SUCCSSS" + + end program read_gbl_mask From 31d54863839e8ee4ef6f7e85d0a9de281093aa3e Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 3 Jan 2025 12:48:37 -0600 Subject: [PATCH 37/67] Convert 'read_global_mask' to use f90 version of netcdf. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 index 51a646779..ac15fd631 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 @@ -589,28 +589,25 @@ end subroutine read_global_orog subroutine read_global_mask(imn, jmn, mask) use orog_utils, only : transpose_mask + use netcdf implicit none - include 'netcdf.inc' - integer, intent(in) :: imn, jmn integer(1), intent(out) :: mask(imn,jmn) - integer :: ncid, fsize, id_var, error - - fsize = 65536 + integer :: ncid, id_var, error print*,"- OPEN AND READ ./landcover.umd.30s.nc" - error=NF__OPEN("./landcover.umd.30s.nc",NF_NOWRITE,fsize,ncid) + error=nf90_open("./landcover.umd.30s.nc",nf90_nowrite,ncid) call netcdf_err(error, 'Open file landcover.umd.30s.nc' ) - error=nf_inq_varid(ncid, 'land_mask', id_var) + error=nf90_inq_varid(ncid, 'land_mask', id_var) call netcdf_err(error, 'Inquire varid of land_mask') - error=nf_get_var_int1(ncid, id_var, mask) + error=nf90_get_var(ncid, id_var, mask) call netcdf_err(error, 'Inquire data of land_mask') - error = nf_close(ncid) + error = nf90_close(ncid) call transpose_mask(imn,jmn,mask) From ddfbc613980397229c449d4ad780c056f9c6c68c Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 3 Jan 2025 13:15:23 -0600 Subject: [PATCH 38/67] Add check to 'read_global_mask' to check the i/j dimensions of the umd land mask file. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 | 27 ++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 index ac15fd631..9f3b743fa 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 @@ -597,16 +597,39 @@ subroutine read_global_mask(imn, jmn, mask) integer(1), intent(out) :: mask(imn,jmn) - integer :: ncid, id_var, error + integer :: ncid, id_var, id_dim, error, idim, jdim print*,"- OPEN AND READ ./landcover.umd.30s.nc" error=nf90_open("./landcover.umd.30s.nc",nf90_nowrite,ncid) - call netcdf_err(error, 'Open file landcover.umd.30s.nc' ) + call netcdf_err(error, 'Opening file landcover.umd.30s.nc' ) + + error=nf90_inq_dimid(ncid, 'idim', id_dim) + call netcdf_err(error, 'Inquire dimid of idim' ) + + error=nf90_inquire_dimension(ncid,id_dim,len=idim) + call netcdf_err(error, 'Reading idim' ) + + if (imn /= idim) then + print*,"FATAL ERROR: i-dimensions do not match." + endif + + error=nf90_inq_dimid(ncid, 'jdim', id_dim) + call netcdf_err(error, 'Inquire dimid of jdim' ) + + error=nf90_inquire_dimension(ncid,id_dim,len=jdim) + call netcdf_err(error, 'Reading jdim' ) + + if (jmn /= jdim) then + print*,"FATAL ERROR: j-dimensions do not match." + endif + error=nf90_inq_varid(ncid, 'land_mask', id_var) call netcdf_err(error, 'Inquire varid of land_mask') + error=nf90_get_var(ncid, id_var, mask) call netcdf_err(error, 'Inquire data of land_mask') + error = nf90_close(ncid) call transpose_mask(imn,jmn,mask) From 53250b2baa4a7066fdcc5b906b3d42a979e4cfdc Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 3 Jan 2025 14:25:30 -0600 Subject: [PATCH 39/67] Add unit test for 'read_global_orog'. Fixes #1000. --- tests/orog/CMakeLists.txt | 8 ++++ .../orog/data/topography.gmted2010.lowres.nc | Bin 0 -> 13565 bytes tests/orog/ftst_read_global_orog.F90 | 36 ++++++++++++++++++ 3 files changed, 44 insertions(+) create mode 100644 tests/orog/data/topography.gmted2010.lowres.nc create mode 100644 tests/orog/ftst_read_global_orog.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index 4fe257e9e..6eb655aa1 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -15,6 +15,10 @@ endif() execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/data/landcover.umd.lowres.nc ${CMAKE_CURRENT_BINARY_DIR}/landcover.umd.30s.nc) +# Note, the "read_global_mask" routine expects the filename "topography.gmted2010.30s.nc". +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/topography.gmted2010.lowres.nc ${CMAKE_CURRENT_BINARY_DIR}/topography.gmted2010.30s.nc) + add_executable(ftst_ll2xyz ftst_ll2xyz.F90) add_test(NAME orog-ftst_ll2xyz COMMAND ftst_ll2xyz) target_link_libraries(ftst_ll2xyz orog_lib) @@ -66,3 +70,7 @@ target_link_libraries(ftst_rm_isolated_pts orog_lib) add_executable(ftst_read_global_mask ftst_read_global_mask.F90) add_test(NAME orog-ftst_read_global_mask COMMAND ftst_read_global_mask) target_link_libraries(ftst_read_global_mask orog_lib) + +add_executable(ftst_read_global_orog ftst_read_global_orog.F90) +add_test(NAME orog-ftst_read_global_orog COMMAND ftst_read_global_orog) +target_link_libraries(ftst_read_global_orog orog_lib) diff --git a/tests/orog/data/topography.gmted2010.lowres.nc b/tests/orog/data/topography.gmted2010.lowres.nc new file mode 100644 index 0000000000000000000000000000000000000000..80c71b1a5da21da3c9c19feebc3efb7a9f05d7f4 GIT binary patch literal 13565 zcmeHNeQaA-6+bU_KI*z@y0s0BHuhFkS~iU9*iQ2yRXFuW+{R7PIGNU|#dGZ!$HTFM z{n9oFt!NV3N-3xmRV-thn)0!ATD1wLfm&@7(k9r(A43xof(#)LAJSBYU}#XmIOpDb zoHkB2CnOOP?<9Wr`Mz_{J@5SPx##-adomh|Y;d$VY8o1V=>tXo^2-{JdhpE?ry2*M zp`n@$_u6WgG?mF*-fPl>Y;(2l?s!^lmXP!$MHj4Bk|t2-n)M|9P)bFJo@m! z;6U$C&tUYR9_ktHaqc888?2Y@CyIspjg*-!2K0~1CuiRhQ~qeiiW z<3?|Hx4C!etHGw-+E687ieRiK9v8+S52n!`+Pgp8AKyFBpNI#0VqxXle^RbJ zMy-)Z4Eb84a_xS@Or=Ntq3n1zUE$oKMWeXha2eTCL>-eA472r6OjJyE)$fT*~3+CLTqsd}6 zc_e4JbtWI&^y!w1mCz}}D4@`^RYEu2D&Jy0boND!uuSM?)aY9lbi%?5K{x#1tknO% zyskcQeh7bV-;XPSW2fY@THstFssYRgZpD4>9*mJh%L0xDJbGsZ8Sgg4BJp?_4yh2X zIW_nju|0F);2$d37pYi%)ONvH*Fz1JsMRZtxU3en7Z*Wo#ocX_Lu$Az)WrP10BZKh z%P%bg#nP%*AGi)Unf}X*61cXH4srUYaJD82fPX#sKR^rAHd64j-KB;!g!ldkt$dy6w}#!5|?T>#w(1*GOWl&sp?0C!B^MI z?|D`2jFD&?Y=kBfJ{LPWmIeey=P3tKmm%~fmG{hn|Yfh#O=XM#cRn(u%XU&p`jpRnu zGHqmvhEaOA{m{QFeQuEuw?ya;DCl(Tj%ZDwBa_O-!$TWHufWjB+CZzU#ifI!TE2mk zXHQR)bdjwMJ!6;ge*F4AOIzR6K}pyw;TWq#r!e={QDC>9INy_4sl$n zNppz+epUgzs*K%hnq2oBmK;F0_a$FxwzrAZ8xX(3`avyOMDc_)2)}(B6wZ7m-kQI z@qNKnM_e@JAlDWtV~%U&{ddM9ag+u87V2eG&)g@a9R}ezYC;Ai+NxGHG;a`bIey^W zTYH9M5vXmyGkP1)kE>?*#sjfdxQMq~;mcPOoj^mp1v(zyF$EF9@Y82!Ps1_6;Q!|@ z&%(=idn=q7FHHlzcMF`Zzx53u^_$_POzurMEf`*j&%S~;7X5EoK;X~Q0=hVA*E~0% zB%tLXE(ZWlE^27dR47P}iSlAb!d@0FhZKhtHEPtr1U2?v+ji|7OS8%~u|RKdASPxB zHCytW!pNLcX!980`U47Q9t>RXPzdv=;QE)<4UY$|XZbuZxSr*(5$;nfG5&6qmIW*e zSQfA>U|GPjfMo&80+t0V3s@GgEMQr{vOx6~pj0()bn>|!r)o}FU>p2Q?LvKth2G6K zD__oJXGcZyJhqa1oRqsbd{6Ybh2E8u#LgV|=%HQJp>QNT6z&g(6Jl3&7m6UxbkawLfMT1cZ?l*X%X zqi5%3c?w@S{QUJA{X*?2?O=@#F>2NBg>xQ@OUifgUO45ZuK!LssZ+_=71F@@d3l^~ zBdpZ;U=6sOa;TS5MNw^FajpL9>!r2PBK9i3pG0{Q?CxfazawGvq9~2ZjM1_3Bh8XCyT{ zcVWZxxr~Bl=R~#xYZEE4J;yGr**T4A!$T2lF3rwKOb6DgCOo1wJ7+fSv}wgArPs8` zvVdg)%L0}KEDQWUTj0-EkzN}NWJZy&e0vt@CxZbW(d&bOr!OL5{dOWOryoJW=TB07 zu`h7tGE%88@CK3D7x)#mpX>{~LvoMy1%z#n_XTunkNkd;i~J!X_%7bfz Date: Fri, 3 Jan 2025 15:08:50 -0600 Subject: [PATCH 40/67] Convert routine 'read_global_orog' to use f90 netcdf library. Add check to ensure grid dimensions of the orog file are correct. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 | 44 ++++++++++++++------ 1 file changed, 32 insertions(+), 12 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 index 9f3b743fa..5de578465 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 @@ -550,28 +550,48 @@ end subroutine read_mdl_grid_file subroutine read_global_orog(imn,jmn,glob) use orog_utils, only : transpose_orog + use netcdf implicit none - include 'netcdf.inc' - integer, intent(in) :: imn, jmn integer*2, intent(out) :: glob(imn,jmn) - integer :: ncid, error, id_var, fsize - - fsize=65536 + integer :: ncid, error, id_dim, id_var, idim, jdim print*,"- OPEN AND READ ./topography.gmted2010.30s.nc" - error=NF__OPEN("./topography.gmted2010.30s.nc", & - NF_NOWRITE,fsize,ncid) - call netcdf_err(error, 'Open file topography.gmted2010.30s.nc' ) - error=nf_inq_varid(ncid, 'topo', id_var) + error=nf90_open("./topography.gmted2010.30s.nc", & + nf90_nowrite, ncid) + call netcdf_err(error, 'Opening file topography.gmted2010.30s.nc' ) + + error=nf90_inq_dimid(ncid, 'idim', id_dim) + call netcdf_err(error, 'Inquire dimid of idim' ) + + error=nf90_inquire_dimension(ncid,id_dim,len=idim) + call netcdf_err(error, 'Reading idim' ) + + if (imn /= idim) then + print*,"FATAL ERROR: i-dimensions do not match." + endif + + error=nf90_inq_dimid(ncid, 'jdim', id_dim) + call netcdf_err(error, 'Inquire dimid of jdim' ) + + error=nf90_inquire_dimension(ncid,id_dim,len=jdim) + call netcdf_err(error, 'Reading jdim' ) + + if (jmn /= jdim) then + print*,"FATAL ERROR: j-dimensions do not match." + endif + + error=nf90_inq_varid(ncid, 'topo', id_var) call netcdf_err(error, 'Inquire varid of topo') - error=nf_get_var_int2(ncid, id_var, glob) - call netcdf_err(error, 'Read topo') - error = nf_close(ncid) + + error=nf90_get_var(ncid, id_var, glob) + call netcdf_err(error, 'Reading topo') + + error = nf90_close(ncid) print*,"- MAX/MIN OF OROGRAPHY DATA ",maxval(glob),minval(glob) From e23b2fd030bac332dc0c9827bb92e49586ad1d4e Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 7 Jan 2025 15:42:41 -0600 Subject: [PATCH 41/67] Add unit test for 'read_mdl_dims'. Fixes #1000. --- tests/orog/CMakeLists.txt | 8 ++++++++ tests/orog/data/C12_grid.tile1.nc | Bin 0 -> 51102 bytes tests/orog/ftst_read_mdl_dims.F90 | 23 +++++++++++++++++++++++ 3 files changed, 31 insertions(+) create mode 100644 tests/orog/data/C12_grid.tile1.nc create mode 100644 tests/orog/ftst_read_mdl_dims.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index 6eb655aa1..c0ec2ea2a 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -19,6 +19,10 @@ execute_process( COMMAND ${CMAKE_COMMAND} -E copy execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/data/topography.gmted2010.lowres.nc ${CMAKE_CURRENT_BINARY_DIR}/topography.gmted2010.30s.nc) +# Note, the "read_mdl_dims" test expects this file +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/C12_grid.tile1.nc ${CMAKE_CURRENT_BINARY_DIR}/C12_grid.tile1.nc) + add_executable(ftst_ll2xyz ftst_ll2xyz.F90) add_test(NAME orog-ftst_ll2xyz COMMAND ftst_ll2xyz) target_link_libraries(ftst_ll2xyz orog_lib) @@ -74,3 +78,7 @@ target_link_libraries(ftst_read_global_mask orog_lib) add_executable(ftst_read_global_orog ftst_read_global_orog.F90) add_test(NAME orog-ftst_read_global_orog COMMAND ftst_read_global_orog) target_link_libraries(ftst_read_global_orog orog_lib) + +add_executable(ftst_read_mdl_dims ftst_read_mdl_dims.F90) +add_test(NAME orog-ftst_read_mdl_dims COMMAND ftst_read_mdl_dims) +target_link_libraries(ftst_read_mdl_dims orog_lib) diff --git a/tests/orog/data/C12_grid.tile1.nc b/tests/orog/data/C12_grid.tile1.nc new file mode 100644 index 0000000000000000000000000000000000000000..5bcfb11514828e109fee9a770c4112b06ef9356c GIT binary patch literal 51102 zcmeHw2S8NUw)Rjg2qG#-vw&hRAfRGF_8FQ5ODteV1XL730To;93ijU77<)JN&bG#0 zV((pJ@4fx!%wFHj@g%+^zMI^4@9(^KU9->Gd#!Knvf3VFdX1ovoH<-_m}k$f!{-xy zX8tUmc+rah{h1Z@Y6LYf&pFo2ys+V+g)P0Sl3RQ^Nc1Vu9OgPR<&B-4^5so-iden{T6kEluZBlI6gy@vS z#Q5a0ZITibQp=Y2a`Uj6rtq|or5WV4k4sKTOzL3})rZD=c27*|=RM9<`?Zn2(?LmEalZd5a@p=V4|oD4HEIXWq>b4s#jxpJPdsj<95NwqoB`#OYTF( z;^@rG?3IV+W@cue=q;#)y)u!x-3uA#)-k)zv9H`Ki(#QGa)^bBDk}z#qwx;;4FfE) z=&(7o)TXr3)k213n5)1rsbNFIOeqVmEJ{$#I>B`t)~r`IvSFYU7EC*$u(D?Ir@bIL zg^jXqwzMPa#->EawDk>&>m;kI-4li}*nz}AA>9^R)vXs;qxK)87TPYWBnCGG4+5o| zdV5$6g^^7&vGKA9I!A6w2%e7q}E(A_OjeDH4<&xzoG{@KN=?JWlW?TdFt9MG>@JcC*I zoy7|rUaTuA{=Xo9_o~T9Qa7XExU>JdRlE34t=cfM8eetfl-({L&fk1mA2aX&n}>ht z!M|_CoCzoBpFL;UJjFomIm`a9Hj@72o4;zs47TidRxGcs#mc|5K)>z$_b*umT}j&` z;lDj(SSrvzyJTI<{XatfR$}!8HM|o&rreg}# z=@8J_6c_Q5mB89AHnCG|3dN|(7F1v+cWxgm$J9~rHYJqVvCvtI&IyT0DeWUWC&mk{ z45;MpUgf;pD$9S^$T6U$sdwk3#1661DTuXm67%S;ZDJ!4vf5OZSJE&Ah9EE|P7Zlv zQ{s9?8HYL0u=vDw3ZnQPkx4SssDyT1I(aUAt#!ViWDA+aR%}0I4#a0>sE7@p zQ6xNoo?{nUe%Z(eD_#H6E>u=Y`G6`0K(X5bf+d-PE=m{WL$*aUbP=W8)9LHi5rY|* zLtrJf`r@8fQVIp0yu;chMRjg3=O=Pz(k?EgYfP+7G3A{x_V1bymy&E)Gvx#;gJNRa zCB?=@#zrNlj5zY(&_A;DDioBK;6TZA~$Wd_g(Egeu=PBD?<@1j(oMLmuk?WAp zfIt^7i1x^JcC}v_(x66t+1@LYtJrn=AT3@&Z6|}>mK(|vC?8PK`n^z^-h57Ei(nin z+=iQCvHZ|We|%eaeWMydviX(!Dd#SxL>Yn0;rz0YeyudGThWpm2C48;N^s7J(+ zwBhCVAq{0-I%pvz3a12#3boYL*QM&*$B6tQy!v|J?iWIr-O8s@-SE3VO#T7xBfUQwwDcP-^tBypMUub7}noZuo5dMwv8-Hte0_% ziAd8hh|LA^HjIv=r*7}1uRE{+b zb7DxE&Tg^}Nj6d*a$&|9aw|AT*wm}mF@O>HNU0>va5 zuJZX4xtXN&`%781a8pR1qfbdQS7X%`w;FU$adk_3Q3 zrQr}gP*zv@;Q2iIf}jV_pHiHIa0k!PA8`+!qdy`ZJV$?oK6pNifFS6>b4KA7q8~g5 ze-)7+341Ak`Z!_?12ENU}w;A|01K(!g+YEf0f&cvssDHeH%fa`QV=H%! z7P_YL4X*MTCWuZKjxYauy4UaF^Un}5AtQ~_ZkrYgqQCMR7NyLod_a|Q27be$sHtu% zwT!|+`&WO4fF>nIey4oIoIImB*(4tRt#fjwXMvQm<1ZQtEU?#?LCIM&P(wAozRsg;H@-BTy=B&1QKY z5;A}Lry3b7U2EzlPW5>K{bp`y?b#jcC>$F&oC8@`e3X+pDaO z@=;URUZwtde5u*QkcJI|b@&yEnS2@kozi=r`H5Z1NYh6_X46~z;xCuZOZf$3gW!$`Yg^Q-HF}ol|M#3V#i%Kgv3iO-el+-XykRbu;&-M0 zs#U=+CqOq-u^+m3rP?^!ST(fU(yb-h2kG}dXiunKy4pXpicN5RAr?pHo1eZA*<#g{buft#7 zQ)ErrXI|Y*9AUYYRr%X^3GHxsVy65c#_%yuNX1Q#nI%rQc|s3d`|07xI-;*zh%|~GiqjuJ0A0xR7|UxB`&>(fzS5<{!rfTWi4m?tW4f`!Ebq~lugzIwyO%eXPIhMT$Ma>5 z{mFx+=%RKMZ(oLG%$xLF_o6&2xcjR0h^ZA=tCmqY6WuB^KR=s(cInla*)Ct!>3D$P z-%k=$)}!y*VT13g^jKiMJkoi&Nsm4E)1Dqw>0#I3X8tgh9($~;8x>UPF>}Y8<`Y$V zwDHe6w}VQLnfZIYJ^9$2H6lH7TuQ9IWVr`RCOu~ObUE$ytUSw4ddv(il%+#KlOEX? zHd+vIEPxrpN|heF+5A~nfu{8sxo}F=N-8}zCR(T8R_Rf8=G>oWsr2xFGVNVyl^(e_ zMCO~O(&J(14|if!dbEui>TQ14v>tI);;pkSHLXV{%bJNppP1I8T=%XaBl4T{SoLgr zmbOQJyB-$$*tW%0dZZrRp1-k5kE$Oo`&hV}_SHzgMm^W5^yr`1@?8~`9$rt9BW9`e zIBGj$VQZBhB?4-EH{hB{k5&&;*Lp4TU|nfF>Xp&$IQYo49=%^*^I7g<(&OyqB~#N6 z1u#SXrmtyVty&_D3smW`ukm`1E-F2C*YbIpU!}(~-mT{Wl^%_P3mmMW(xZ}JztVG6 zdhGkT?}_FrJ(`8Ld@}XCNss1ti~lrczG*$yH#{)I>!E2q;^z!o^T5fZ$C;p;=l2}= z?Rwm6awxKeN{`Uk&4-7n^k`q|R@(|HJxaKgG`phGqx_i7lfzYd#Ck2WUZm1v`jH0D z!&Q28zFn#CniD2H9&NgC?$aF8dh~3XF?P}Y@~i-zj}c31+BbBp%qr12Z_=fdMV0-( z)uTzd0>x8QdL$**tTsob$EQy#285~f=n=4>&NG!BF$Ja!h*s&5dgS`lWhy<=yz;kg zpwffad9u3tVUr%`?VM^)nr&JSIRiXd>7Hplu4Ip^deFhN9;q&6M{V8vTRrwgS^qpl zrN_YUyS-nh(&I|aTstCFdMwOSYrn;Glil^~JCi4<^mw}d-n5@odgQ8MA6r+YM@)OZ za@Jmxy?U{D>$)McOnSVFvc7C_*R&po9jEw>u{W(p^7Wdd$N&6WJu00Za(;?R5APMH zbM9B^ad~&c>Pad+u3VouGOtRHy}O*c_fY8(W&JpAtxAts&vQ+xrP6~>&qxm1X<85S zL4kEfWtjA^_;IJpgYC;mEwE^%yx~Txims-|BIC?%L?NDm~7|{@Ul1N{{G- zhoiq&>9Nt$Q!1#^giMnjAKMo?S80Yx zkNwWhW{xU7y8PI~eM%nFdW>B)Z`a1%ztyAbJ+@(mN{?~PdM~-A(&NbTwiAb`^q83E z)Z!8dn!FncABgkuhQe<=rY;7RCYFdw_h0|Xixn)|9B_Ce~*H!7^l0Pd; z`yD<0>GjC}2kQ~~mFw}}az6gA*{cJ8V6R4h!Cw8h_>X_LKY#iM{`}CN^ymLA@zp;a z?SF~B^jD}iM6I)lA)=WWcK!KY_Syuj@-edaaM^E+ul%YrNedB@-ZnH%Q>@(~fYuUk|M;*lTQ`MPfk;z#z^ z_P6&8<}=Ez8hdwqFn^qNtoN?l!F)rXS+VC!hVWH0S0|5X8p6YSPFT{STL{-FpMIt` z-lIRpPa}Td4|$LeJizBlc!7UVRIMs2OENzF!pP0X&M@A)PvH;c`$#-ui@8U*qn;l= zK562Lxq5yo=eizoMFV;L+{@>MP7dUGRvbwx{Wg$WnYI2{CoG8PYwVXNVPp`mdwKJl z23vyol+4@pj-3hOxos`ng07j`c<)60bI|zQ#1H%-5AuNr_<$Gq>-^N@QtoJf-hOt_ zt6M)+=PoyN)2iPO;G=WasdjA`<978!LYwE6_=4uo?PrXU_^}8(+g_g}-Y7G-Z-dr) zUgTZ5{ij#z`TlwRJ70gS=l2Tj>g!x6kgqI%`(v@Hf!s}ruY&^lPt?ZyebgW0!4Lc) z5AuNr_<$Gq$IZ%LV`owo-s9uB;wu(c<)dz4IcjKsU1YuSEl zUx_#E`>ybUG>NBB8}HE{<3ABU@E=3+ARl;u4|sw9%*ZQ0r3|RVS6fX9Yd^=2X9>7` zbjhm9{NtGDp@-L3;WsnAq?QY+^6vXi9G^C%8qfN$q2wFp&u8Bl`SI;rf8IB9$D7(C zs`J*H*Js?<1@JLT-Z|xN9Kb`|)+fZI1@Iz<`cnX}L2bNmMg1`z{J_5l$%B000Y2ab z{`|e?x)vMf%ZGTkoZn|+1%7b+InR4(75O5MDou~4RN`GX9lTMzjvt?r?Y?s!r^gy6AGG*XgQ!l_9S) zweem@{V{$Y@taTlArJDG5T1R6Z!qBn{x$}CRi0=1T&|ox)SHj~_2j5hkv@FVi8|>! z%KLK5we8cMKK126d+VCDT~vVwP2W-Xdaa5)eXV1Mrspg24uyu+Oleq&H&{RS;=DDL z_@%IeO?SSn#N#ivir84xj|cSgx|-M5&)CL$^vC#N#1H%fNgm__5AXpm@QZU&mR}fK z(6#xm<#>molLlQa<;9;mGn>+vy!drLTbBs~%JY{)?#`KB(3^kA$n|9ULT_GgL+ZWx z9zMMG)@gly9PYzQ4)W~~dD(}b**mKJE5#GFaR_?c#t&j; zA>jpnp^qbf+|y~vuzgOvZT+cUQ{FrC$o3vPel6?Di{(1#{4O{@?;Y)#mb+mAUf7M> zRh0_zBd?tY_I57BKU$SO7<8x*-~O!l75Tr!&4;;lE+oae@l$uU=Xf>Pjjx~kw8Ql= zrZ(QAKgJ&;e&7#zkRL^OfDd?q|EytOW#e`wvz;HgFFWaB#oN#5)AxB=><&@jd!uJowRi80A4e@Bkn10zd69M*EHQaAtHqNRRxC_ABX8kkR>I zw7(hccf#q$=zNeKrZ$}qV}FdN^I_x=h$>#mHVUvRB6b7*F=k$e-+$Q9kehAMgUd(5E;fd(6l_ zGqTr=>^CF-!N|TdviFSaKO=v@$UiXhAB_Bmv5oiWkMZCK{^UQ5@_`5VfEW15AF^}8 zpEL58jQl4ff6mCiGV-^K{5d0k%*a18^4Dyd@aK&DIV1mWe2@MZUq|?JBmZZ@pEL63 zMm)er{+t0n#V3s76-Mz3qj-iOM~vbnM)4D)_=-_{Wo+X;`eQuB zS4RF6Um4{C5AXpm@Ke0VDBfih?=p&a8O587;!j5LE~EIAQM}41-enZevMVCqWfbok z+jx)u7{6b{yNu#pHd(~iM)?%)GKzPNc!8hdcSi9%qxhatyw529XOw?1$`2Uj3yksy zM)?P${DM*b!6^S=lz$lCQ~qHb4}RcJ`G-+H@Bkn10zc(bELG%#jPfl;`4^*nj8Q(w zC|_fgzcI?^80B}2@;yfRAftScQ9fvVKThP2#_`|>{*(_IM)e6s^$JGy3r6(}M)eX#^%6$)4@UJ6M)eX#^%7$n@6n&?B}RVW4|$Le zJirIMz)$rbM)fO3^&v*}B1ZKiM)f2{^(992D@OGvM)fF0^(jX6D`p7Y`h^(Vc#r-V zPxUJ!f5?M;-~m3WUoqfcJoe@8yHD~;Gl)(kI-Y0;qHBrHB084n1fqwCCK1gb+JxwO zqSJ`3COVxc0rKdr5`l6fdrSqXpCLg_QFV%~=aQvlAdwKTieB+J1H2ZXZ+gCCzkY#t? zw%Ht{eve#xcT`b4??C+)(YR*B=P>a-O>(-BT(|ZDM-bs!L^z##3fw;@1vO6GVlUYx z3!2_!XX!(Q9Hfy%+g0g#>0oyU=_pa}gP*$8%yf{N5PfmAw(IUE4pP=@g3dZo;m~dC ze9}px=|t}lElB-7(YQe3Q-kPKl9NGnAmJEC^e*9jKs1kQ;TqE|9i%Bg_pY6q#X<7! z=Uw8#GzY0($lLx4r#ncc<_~XdRV1IZ)G5MxhD;|9m}$0qd_F1Xkl9gXCghW>D_rU~ z?M6Q7=QlfiPRq2jd+7n!Z5^fQcUpIPCDS@Yx6$(v)K8*ud5F&h;yZxkgp%BfgrhOx znnO50x(M9*y;IT`ggQuji3Zsg+%nN7pHz$J#?lE-mhgPiMWVgCjP%GG;V6wHdSUs= zzLR!3N^^;ZJS#UP?Y)efXc*B@qMNDT4x)pI&k&+JNX}8BA%vqQ(SC$;7*U?yBFDn6 z4${MQ9UdR*<{}F;?wzi|xtk8wvb~&7+ zT@xl=bF+4m28LcJmOjZz+E(#c^N%tueX3#bwkuAOe`xr3w^S6*Pf@?gG;R*@=}&yW zBRP<}hj6qfTyF{IYQkN6!AkRKpsT*i2${XxQF3h{=<%^O zEsrj8l3oz)Sn=ePrQbPAONst`%&Gb1#m>?d4aIZxA5HzQ5}(^dFOs~6M6C!%UZTl_ zvpZ3J%^|h#Zgh}N#uq;rw%I|_SzoQu_iaAu(wu~yH{a!x=0DHuc)qlg^j-fqx{_s` zq?b`H2V;GmrJR*Ux~Eofmb{mpDEj)fv-F|<%2CVSI7=ID9DY}^jf>RAH)~oq6~*(7 z)GwLFeX1z-J!Ky=2o z5BYKpaFU!e13$KN)=K>}3aOP?rZY=Y~5<7UL{x*J|P@NiNfi zRoLkybx-Ya{>UyTDQd@(T@Nf=r0Ku*XjmtUi?propzQuHT%<6|)6>elbdib`+!i&U zyQ@@rMZcw|WqN+Opo8dnR_d3fo*1{0_~?nRAITX=a`%lDI93y`vV;@3ZzW@b&)0#op5tPoywuE(kr4X$8CvU zbIeup8ZYQkYTHr!C+fGK=yKw-f~XtG$wPE1;Q%d6I6a9@*qmCU+XDwFsmsj%E{_}} zkBg~sHwHLL9n1oTEgj@2rFk|Q=$gw}@_U^-&ow(|Dg6BEZi5nBq>7iK`sD8HBE{c% zP;f=6tMsf%gNC~CD!){!pqp=7neHzjXtT@W`BLgPfyQ|gp9tb>O>%xB zx#@%>8{uk7IF}Lb@D-2uXMgMRb$SD^XFx@MvIU;iiP|37{I;~W>$*dNP(ll`eq_NVt3?2jAaTuFBGKi&Q_@?)pI z;D6_lA8hnD`QJC>e{+)m%}H{uliY>me|^8;fB)0tpT8P^Mo^r%_bcPilN5hWq4;w) z_z>R)6n{>jxK|s0zM=SYEybVz>G6Lq%D<-l$^2_7kIkUe|r9_tv@(Yo*zVc{~y&KzNb3E%)eNFSn`GX!+EMd#8dsjm2lMjLjB=C zz5X*vQwPF+LR|vq8+Hoz7I_=|AN)At5yVx93(>}V^v8JkA@GMh$j3PaKJ1$*g#SmP zsAHiXgt`jq6Ug(CFC*_n{)RjV`2^x}=z%ugBW}TX_)+kOJjjQg0Y2cJL-?^CsKcS2 zg}M;xE2vW--$&kz{1bT`@*(6Eh|kf+d-TV6@B@FygM8otKAcMW>NAn!;1j64+i9P%RM7ii->;zo=IKk$b<$Oj(agWW-07W)cz2=zqNHx^6k(VOBLmThWAL9{6fM*Ehpe}&? z8g0Bse~bq|@JCz;`G|9X4|s7tpa<>&P)|o)8ueY&Nl`yU-469H)R9n+L0tp&0krWR z{V^W=z#sA;A9#Qdal;1MS5t}N9s~CQsH>wsjXE#trKmfieup|3>Pe`}puT}N-lIRp zgCF=q9^?ZL@FDJ5Liphiao>Xb4BQK#zK(h{>bj_pqRxmq9_nGJE1@2PHr}H@#)BXD zLmuP<5AXpm?yYc6aF2xh7~{PL-3Oq~j(RofzNnv~4vBgm>SCxbp^f+GkMZCK{*Z_J z3E%-f;O$B4q16NTNw}B6{Ra3#4&(v{>ZqtkqOOPaKpXGTALGFf{2>qWfd}}Y$M7%g zE8K_SUJ3UzSZAy|_5t<<>cFU{qArQ`KpXGTALGFf{2>qWfd}}ouLh7FIOn*B#C;y_ z#c;obdm7ji+#8_&jyg8#!KkaEK8ZHoqd&%jANV6LhkWFTzz4iI=dioDN5y?4?)7j# zhI=O5*Wlg*_XnuMqn?erFzTzQQ{p}PV?6kQKk5pQk31LnkaxlELJ!;n<31JllDOZ) zJsIwsaBqYA58NZ59*?><>cgnB;ywCfJnB&3kGcrt0}t>4FYeb7C*mF(_rbVV#r-7i z`EXx`dnerA;2s3`38>4XzKys9@6jLQQAY#+fmGLleBc2-;ML|IxQEAmHtvOSzlwWG z-1p($4EImC$H9FF?iEm4|NCpjstN!eg}cy6X0GR z_uHt`;l363mbm}JJsR$zaIb@V5w!6h{V^VOSn!8D$VZ(D_)s^QM|Kx_;CCGOJp_JN zfctsebK|}k_pZ3V#62Kz;$8}RppEzFkMXD@gFoa!KI&}1hq@2)TI?(Q4h6sG!0#gP z`vTn4>avgzJirIMh)1v<_#F*?kAmNI z;P(;uodNFaaqq3YKgN9~?geqbhBn^g-UQ>p5ByPAhJ4@wKH!}Rf8M-H{uga4>Tmb> zH8OXVO8VNtOLJa!_0yj^zhS5UQ9ph4R-IlK8(mqSKd{lw>26i@2?MWuj9ynoZ!`Sq zf&llb`pOM^S@%h+s(;&YGe5Jhsy_5ZuGeWFtLk@bv$&9EU(HBu?u+N>m(oIv`=!6& z1HO<0xw&r&9E0l$T);V1gS+>qg)jwHXbV9eM)%24$9e6%}mA`&i;I8Rm=c?=d{91N(jt|i1e>ihi9yg|U@tPL5 zq7Ktr4KG@$$Y7?A&K+B5=>(?VZ6he2qaVhBPrWLF?=L%r9LNO@;PSj9<`?UM^?{zo z{*lA$R_a(yKgBtZ!{Kbz^%-Rg3~ao;y8ijlPrbc!GX1mUL(7{rWcugx3s~np%Jg;? zo@*|@WBRS#p@z zC(s+?r%xODQ?}ep?=iQ0@a%M^U$x$xt1NGJ)e>u9lxQXI9`hA=w#=Q*^e86|;EFtIPL?IVA%4`+;062jg=LP!? z>xX#&?wW^oEO_>j>En9)ZoH$@>qn#?9un0;uP=Ij+p#wJWt>js?6-#p>OC_a=tgc0 z)JIuNx%1<(K>ep|g5o*)wZ0+7Jsc$XWa}^ZLJsr*4(wmx1a6#roFD9OtS9CLxEGa5 zUv$7fP(N4Ktn|^9f%<;+9PT>I4$^mTz4%3z`$75#2PO^oNeI^OscY3XB+W#J9Tm^f zuRtR)?vj<@1HO<0xxfKj!0FRd*g@DY*nwn8oXat(Vt-?wV4kEbLSCbYYegHm2k8S= z%&5E3BUryKcb+SEItA-(ENX4L`8HVJ?Zm?GBJ4wqG~r`S!*lezxLS-`-$3vIU&w)6 zoDbjvPT+?9hCSQ#RM-=oFYHt76U;Yo=W4zu{Zjo9{l)`@Zk}8aqW_S-uSj%QsJ`zC z?Mj9hb$v!nq3nQP>CIg1v_SfInHZMc8lHH`o)LH|$%i zKjsm5>J@r5ZPbSl{Y>Xadn^`*>YvV8^W1uGsJ`(mLGgT*i|7||ROkyn_OW9BLC(|J z!X7NoChQt;!GC17`P^SRt`YtPc6w!RVFzJ1aQ?85u^%wsRSOE-QT_844T-6t@35wK z!ODYc=nt(BbXRZj9Q~$^66*mz;CpbkIRB6fy9-?KkMQSdYlWYOpM>98U@7ePBJ<-0 zI|zFc(@>mC>~HK7%r|h4X@7o(+58%MQ0Rg6z&P;1d50Xx1rGQ@_;ti5e-%IEKrV0q zKX8`q^!xbH597cG^Pzy+MZ4f_eZ6!GeF z`w2Vtwe}PH1ojhl2lf+w0CpL2V7Fn%VLxH#VL#yq;1}R0zSe)io?u^}JtF)Vo@3pB z8-DPw`mf2UpXUwm%ka~PGvUYK*Wu^C6hB}eh(rTFdkvCrc-_<5{1_ASmI>c;zq=8zy+L$-w=l)-&wa< z#Hol|5yyfb<^}Ny&IkPT?){&~&#>coj(*tRh@Tx+i?|$dC*pSC04{C(j64AOBJzX* zU&vn&4A|C`!tpUoeU=OTYZ9*n#g`6cwkdi_!Uj`$CL5PlN=9?xOFVFzLNI;;`r z68jtb890DTdMWA*z^$!AKrgHh&Iis5&OPcUu)`I)ejX3Pzk|ZhqaXYZ>^JyAjzwcp zx554fE-mf|+sUu4RrcfOvsc|z*sBuXcE9D>>-j72`lX9>t#9VT-`iwcJN}9nuYJ8k zpPaqQ@*%4%R^_cwl5Z=PKIXtv55C~iRga`~Zahz7lbMIjocY2TsWXPovEy@U<~Z*? zQ^)b2fS@)jf1qX3wcc%6xfPSVTJ>tn4v6uL__5A{KMND`*aIP-5gvA4;A4cB*$e!v zWPG*brGhGPZ!hb;CF)h=&ztuw?HuLHf9{w4W7Rls-qmk|2@@Ub!iFMHe3fS;Ay-hA2a z+7$8mYQ=J@a%4bNqPlpnG zvq!E&afLnjtCMTvv$(tQ9mP5qI#A7tckWrQO{aAVk@11>fHF<&_zoUl4;%CJ(rYR-0dFI+ChL!Ihn zODkQ3?^$!W{f=D)d9MaFg0GiwiZk{a0Xa3Tn%#mG1xZSc} z;ug(w;|WEI_I!HSg?pEHSI3jtbCImA)gfvHQ-@{kEIH{%u?WIGotQfxn9!u&N=+3t#EqE}YbV2^Z zy!@Z4G^!eHO-O_eX%=WxC>Ot**n{zD|t&;WVeDrmdDY71Y z#dt>i*l@w0ofPs|TOpqj9#%u(V}zGA6Zl!rB);YQv8CkyZe91Yy;qjxBVK2`k33bJ zkB@uOsqvkn{NeLu-`zf1gqK^|kF{Rl&U5$;_KaLwke`iRRbcr8XYRbSv2N=_d-*@h zUhtvi`dmE8=iI)L)w1xFGq>+P*z|=_kG!)C?SMKK&jwk3-vSh5>m(_qKyD{Tt!lI`XmS~%XBpM6rKY^5a*eE){$;q?dQ=KEgQ zwzQp;gKx=fUf%Vn8DF~W)QPSSo*4C@c0jH{c?X8edSrBclssJ44_uLtcX1ehM0sW)4IlA#H+fR77KP$j_*0LoAX1MaOvO`Yj ztett`&9+%)w9d!Nzu7x&Q}4X|-d>$^(Lr|n^p^ISVe517YzH#?bPToR;kt>tt>tsE zrq!KsMO)m<%vPaw>oq%_{2#V7v;(Sjznr{T*27`*g#B(Z-6qB};>XGe{>)v-W6453 z`&r;&TLeC~PvB)<0zZ3l@VhP}CgtM|>lM3M>$(FUwfxDACvENd#Wp$hR%i2YtKAm1 z^9$zYH}+&X>%730cenJ&dwWw3-fY6`w8=kZD$Uc9o@gJrA+6uX#{{rEgNgwFvNBTm4C$S!k)+4}|=8M)NfYW@^dIV&M^$5WHKkhZA z+N)i!w_tr}JvglgUqkCj>ycTH)|=KN^SxLPwN7_$xK8*W~&WCY7()nO? zJ{X-3<9?>|!RUMh(D?|!dB9ZTywLgJbUrwp594{G^TFwSa5^8x^GfH#c%JEe@cMM# z$zHJ?WDm$*v1Vj1$X>DPWKUpUU~kA?F|tKRU&&x z_KJ5Sdr9_+&m((E_KJs+y(Rx)w8yZ|u-D{282JxI{=;bR$$v2NAB_A5gTKIfz@L!+ z;N(9z`46K%BLBh3e{k|2Mt=tX27gEXgSQ}mNdBBnCx1!)oc%)nl>9lXPW~4Dm;5;+ zf6kcj=ZyS$z^Z#qt4D5k@+W_uN&cLZKj)Px9zc8$Nbv&sbDl}@1o?B`m*Nfb=bZdG zr})YkuTXqtjAtmmViaF7im#0E5XD!F;wwh+6+=9Q^+3Eu@fD}|ic@@LjMpf>;uK$T zim#0E9>rJ2c#z^NPVp|Ic-I(DQoPG3-enZ;vN(!IDc)r#C|;#_m)%)oh<6#{Tj+uK zm*QPc@h)FX@iN7`d=SOc6z}qi6mL_!%hM?yr+Am=p?ID04`Vz}`G+yyr~HFa{=q2! zFy;%Ce=y2F808-f`3B`5{>Vot|KOB=aLPZ7`3&VBobnG&`G+wdqWr^{FH!!%{V1P$ zMfuealy6Z!$S5CVln=7Tl&?`f$i`7VNBJP5e2^jkqkM2H@#F8b58lWG2f^BoKt?zDL*&X3n)J~))Oc{=Tt9YR4-xqsUAV~5*AMN3aXbd zs+TaTm#|q>@1S}KvlR6bhWZH32kIwOFX2=#;b~NFp?V4Tp?VC}OL!F3Yp7nra|wK$ z>Lr}&C7kM4#(EIduZ;B~s$VgxUoom*8S70{zhYFsVpP9ks8>M`)U&95#i@S9seWaw zhf)2CQ~iol{mNKRqxzMx-bVE+ZjjS2oDIA?D%84LI755s+>I+1s@nVdZTw+hchmbL zQClAl=%H#a$a(Z=FH<}8`P0fHzc;meTIbx8)Zf&u-sGvL=b&)<_YlLG=S+R>K10Hp zg#tjvQ@U?RWW}^_X80Gc<#wg!A@w8DP2)#7blLD+ZogLk0l&Vc=imK&DKtujd(d8& z)fOt;=j(6i<(_Klf9qYHkV7gxuJ8Z(rFn1D`=))Rw=-3GEbO-}OO#3v`wD(m9|oHG z7qeWo(OcF7^H_RC=%7nOP3`Y1-3)s#>xK7Q{M9K1-!oGnjJ(@){KafM#!|;!h%j*R_`sp4wPE_H3>vL$rc~!ey z@0stntMtfSeN5%bDn0UykA4-P(qro4hr?H?^r+c5cS;?V9#`XMWalb9?ru-qvq#ng zcyVxXVKqx%;Ycb~m5AcbKZZb**hiL)AVil+!|* zuhJtt=Yv`UReE&y9A9&vN{>!Q5=g?&)VK9J)S&$zpbfC55HPZ+ANj#8P;di zhxV1a%6j4a$b0((zE2f1l^$Dnj~ecyIxoFDI_4V}uCVO>Jm(9DK%6=_6xBgI~O0ORIhi}oT^k}&zyF+jDzV?IG1q##UA1E?fO=e%1U(U)66)RO#_1=u%M^l|9V9+5Oy2)p?v% zaBJ`wSr5pc^~?4w*H!D&YkR3_yJWrazU}EJMJucH8a>>7V1%lTe!qIIjbE+8nQQ*l z8v9l48dp1QQ`@7mnLRh0SLso!(YTlQReEGeOuiMY(j(v45zd=bdaPg4$ttf(kDuA# zP!E+JYYL8&6nplb_WzTn1=m+l>1A`Y*7z`08+tTfQlZI86|TfRxpr(;wGY;rSmdLs z9az5T-k~boea79|)Izn7(q7Km7_ZV}oK@?P@hUwAcXAG%quOWl4xg)(s?y76Rfh{^ zD!pb@e(WElYOiawe?dJ}8+r^`U0OP*!gcz6o^B^p?J=DiI-99*c6Tg&T#b9nv!8wL zs@^|}YMC-qrAMBbKMs#l=~30j@1TQ9kMX~rD5KV^-;M)fymRPB$y9x0birRS%+ou=nkwXaxhELBL=h93WJ zsvHA#oP>4wtGHfh^fSNrx&9V+KA)d+n)y4cnP2VvYu87+et9(OYvK6$`kOm^zCYYF z`v>~=)a)-lGwoNjhiUe6kY>MY&x7{7?56W1{{8Q8mY_Lrq~F2U0$O_Ek;3uvn8%q-*kx37UK) zT9dEXY4Vw$Dc^y8N~C-U_6zwE&JS?)*W_2(H2K$YO@8)BlfOOGl`7-ip^q;E9x4UTa@gz;Y-j?!t*sHvh z?_>T;YU%-jntH(kswZgcCCfGShu)g{L@`bM!c|k>D5$A_XzL@|`pJDwedU~{{&G@N zpBbR3-vn#wJDATtR1ZQskEUMaNA)DUFHiL*w7sbwh4w~Gy=sW2p0$JOU9eBbsUC*? mf0pWHnD5=1dRh|I+fXk_?i3Xt9~m8&6dfOH!0`XqAO8;keCB5W literal 0 HcmV?d00001 diff --git a/tests/orog/ftst_read_mdl_dims.F90 b/tests/orog/ftst_read_mdl_dims.F90 new file mode 100644 index 000000000..3f1aee111 --- /dev/null +++ b/tests/orog/ftst_read_mdl_dims.F90 @@ -0,0 +1,23 @@ + program read_model_dims + + use io_utils, only : read_mdl_dims + + implicit none + + character(len=19) :: mdl_grid_file + + integer :: im, jm + + print*,"- Begin test of routine read_mdl_dims." + + mdl_grid_file="./C12_grid.tile1.nc" + + call read_mdl_dims(mdl_grid_file,im,jm) + + print*,'im/jm ',im,jm + + print*,"OK" + + print*,"SUCCESS" + + end program read_model_dims From 06a6c30f0b54c0b106af7a5b406895bb0cf304f1 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Wed, 8 Jan 2025 07:17:19 -0600 Subject: [PATCH 42/67] Finish unit test for 'read_mdl_dims'. Update prolog for several tests. Fixes #1000. --- tests/orog/ftst_get_xnsum.F90 | 2 ++ tests/orog/ftst_get_xnsum2.F90 | 2 ++ tests/orog/ftst_get_xnsum3.F90 | 2 ++ tests/orog/ftst_read_global_mask.F90 | 2 ++ tests/orog/ftst_read_global_orog.F90 | 2 ++ tests/orog/ftst_read_mdl_dims.F90 | 11 ++++++++++- 6 files changed, 20 insertions(+), 1 deletion(-) diff --git a/tests/orog/ftst_get_xnsum.F90 b/tests/orog/ftst_get_xnsum.F90 index 0de97a43b..8ad59a27e 100644 --- a/tests/orog/ftst_get_xnsum.F90 +++ b/tests/orog/ftst_get_xnsum.F90 @@ -4,6 +4,8 @@ program check_get_xnsum ! the number of high-resolution orography data points ! within a model grid box that are above the average ! terrain height. +! +! Author George Gayno NCEP/EMC use orog_utils, only : get_xnsum diff --git a/tests/orog/ftst_get_xnsum2.F90 b/tests/orog/ftst_get_xnsum2.F90 index 169db1d36..3ad78d395 100644 --- a/tests/orog/ftst_get_xnsum2.F90 +++ b/tests/orog/ftst_get_xnsum2.F90 @@ -5,6 +5,8 @@ program check_get_xnsum2 ! model grid box, and counts the number of high-resolution ! points higher than a critical height. The critical ! height is a function of the standard deviation of height. +! +! Author George Gayno NCEP/EMC use orog_utils, only : get_xnsum2 diff --git a/tests/orog/ftst_get_xnsum3.F90 b/tests/orog/ftst_get_xnsum3.F90 index 45fd76866..b0f46aa02 100644 --- a/tests/orog/ftst_get_xnsum3.F90 +++ b/tests/orog/ftst_get_xnsum3.F90 @@ -6,6 +6,8 @@ program check_get_xnsum3 ! points higher than a critical height. The critical ! height is passed into routine get_xnsum3, whereas in ! get_xnsum2 it is computed. +! +! Author George Gayno NCEP/EMC use orog_utils, only : get_xnsum3 diff --git a/tests/orog/ftst_read_global_mask.F90 b/tests/orog/ftst_read_global_mask.F90 index 9a4a67229..53af1fba3 100644 --- a/tests/orog/ftst_read_global_mask.F90 +++ b/tests/orog/ftst_read_global_mask.F90 @@ -2,6 +2,8 @@ program read_gbl_mask ! Test routine "read_global_mask" using a reduced-size ! version (6 x 3 vs 43200 x 21600) of the umd land mask file. +! +! Author George Gayno NCEP/EMC use io_utils, only : read_global_mask diff --git a/tests/orog/ftst_read_global_orog.F90 b/tests/orog/ftst_read_global_orog.F90 index b382a9740..a5a854c0a 100644 --- a/tests/orog/ftst_read_global_orog.F90 +++ b/tests/orog/ftst_read_global_orog.F90 @@ -2,6 +2,8 @@ program read_gbl_orog ! Test routine "read_global_orog" using a reduced-size ! version (6 x 3 vs 43200 x 21600) of the gmted 2010 orog file. +! +! Author George Gayno NCEP/EMC use io_utils, only : read_global_orog diff --git a/tests/orog/ftst_read_mdl_dims.F90 b/tests/orog/ftst_read_mdl_dims.F90 index 3f1aee111..a072e9652 100644 --- a/tests/orog/ftst_read_mdl_dims.F90 +++ b/tests/orog/ftst_read_mdl_dims.F90 @@ -1,5 +1,13 @@ program read_model_dims +! Test routine "read_mdl_dims", which gets the +! i/j dimensions of the model tile from the +! 'grid' file. This test uses a global uniform +! C12 'grid' file. The dimensions returned from +! the routine should be i=12, j=12. +! +! Author George Gayno NCEP/EMC + use io_utils, only : read_mdl_dims implicit none @@ -14,7 +22,8 @@ program read_model_dims call read_mdl_dims(mdl_grid_file,im,jm) - print*,'im/jm ',im,jm + if (im /= 12) stop 3 + if (jm /= 12) stop 6 print*,"OK" From 5d560797324b825f9c657b8c985c7b072989ba77 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Wed, 8 Jan 2025 10:04:24 -0600 Subject: [PATCH 43/67] Convert routine 'read_mdl_dims' to use the f90 versions of netcdf. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 index 5de578465..14d4fc02d 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 @@ -398,33 +398,32 @@ end subroutine read_mask !! @author George Gayno NOAA/EMC subroutine read_mdl_dims(mdl_grid_file, im, jm) + use netcdf + implicit none - include "netcdf.inc" character(len=*), intent(in) :: mdl_grid_file integer, intent(out) :: im, jm - integer ncid, error, fsize, id_dim, nx, ny - - fsize = 66536 + integer ncid, error, id_dim, nx, ny print*, "- READ MDL GRID DIMENSIONS FROM= ", trim(mdl_grid_file) - error=NF__OPEN(mdl_grid_file,NF_NOWRITE,fsize,ncid) + error=nf90_open(mdl_grid_file, nf90_nowrite, ncid) call netcdf_err(error, 'Opening file '//trim(mdl_grid_file) ) - error=nf_inq_dimid(ncid, 'nx', id_dim) + error=nf90_inq_dimid(ncid, 'nx', id_dim) call netcdf_err(error, 'inquire dimension nx from file '// trim(mdl_grid_file) ) - error=nf_inq_dimlen(ncid,id_dim,nx) + error=nf90_inquire_dimension(ncid, id_dim, len=nx) call netcdf_err(error, 'inquire nx from file '//trim(mdl_grid_file) ) - error=nf_inq_dimid(ncid, 'ny', id_dim) + error=nf90_inq_dimid(ncid, 'ny', id_dim) call netcdf_err(error, 'inquire dimension ny from file '// trim(mdl_grid_file) ) - error=nf_inq_dimlen(ncid,id_dim,ny) + error=nf90_inquire_dimension(ncid, id_dim, len=ny) call netcdf_err(error, 'inquire ny from file '//trim(mdl_grid_file) ) - error=nf_close(ncid) + error=nf90_close(ncid) IM = nx/2 JM = ny/2 From e7d2ae6d94000f993901491bce0c65ba2d9bf9f3 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Wed, 8 Jan 2025 10:07:50 -0600 Subject: [PATCH 44/67] Baseline test for routine 'read_mdl_grid_file'. Fixes #1000. --- tests/orog/CMakeLists.txt | 6 +- tests/orog/ftst_read_mdl_grid_file.F90 | 77 ++++++++++++++++++++++++++ 2 files changed, 82 insertions(+), 1 deletion(-) create mode 100644 tests/orog/ftst_read_mdl_grid_file.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index c0ec2ea2a..7986143ea 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -19,7 +19,7 @@ execute_process( COMMAND ${CMAKE_COMMAND} -E copy execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/data/topography.gmted2010.lowres.nc ${CMAKE_CURRENT_BINARY_DIR}/topography.gmted2010.30s.nc) -# Note, the "read_mdl_dims" test expects this file +# Note, the "read_mdl_dims" and "read_mdl_grid_file" tests expects this file execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/data/C12_grid.tile1.nc ${CMAKE_CURRENT_BINARY_DIR}/C12_grid.tile1.nc) @@ -82,3 +82,7 @@ target_link_libraries(ftst_read_global_orog orog_lib) add_executable(ftst_read_mdl_dims ftst_read_mdl_dims.F90) add_test(NAME orog-ftst_read_mdl_dims COMMAND ftst_read_mdl_dims) target_link_libraries(ftst_read_mdl_dims orog_lib) + +add_executable(ftst_read_mdl_grid_file ftst_read_mdl_grid_file.F90) +add_test(NAME orog-ftst_read_mdl_grid_file COMMAND ftst_read_mdl_grid_file) +target_link_libraries(ftst_read_mdl_grid_file orog_lib) diff --git a/tests/orog/ftst_read_mdl_grid_file.F90 b/tests/orog/ftst_read_mdl_grid_file.F90 new file mode 100644 index 000000000..952212be5 --- /dev/null +++ b/tests/orog/ftst_read_mdl_grid_file.F90 @@ -0,0 +1,77 @@ + program read_model_grid_file + +! Test routine 'read_mdl_grid_file' by reading a +! C12 'grid' file from a global uniform grid. + +! Author George Gayno NCEP/EMC + + use io_utils, only : read_mdl_grid_file + + implicit none + + character(len=19) :: mdl_grid_file + + integer, parameter :: im = 12 + integer, parameter :: jm = 12 + integer :: i, j + + logical :: is_north_pole(im,jm) + logical :: is_south_pole(im,jm) + + real, parameter :: EPSILON=0.01 + real :: geolat(im,jm) + real :: geolat_c(im+1,jm+1) + real :: geolon(im,jm) + real :: geolon_c(im+1,jm+1) + real :: dx(im,jm), dy(im,jm) + + print*,"- Begin test of routine read_mdl_grid_file." + + mdl_grid_file="./C12_grid.tile1.nc" + +! Initialize all output variables to flag values. + + is_north_pole=.true. + is_south_pole=.true. + + geolat = -999.9 + geolat_c = -999.9 + geolon = -999.9 + geolon_c = -999.9 + dx = -999.9 + dy = -999.9 + + call read_mdl_grid_file(mdl_grid_file, im, jm, & + geolon, geolon_c, geolat, geolat_c, dx, dy, & + is_north_pole, is_south_pole) + +! Check values at selected points. + + if (abs(geolon_c(1,1)-305.0) > EPSILON) stop 2 + if (abs(geolon_c(13,13)-35.0) > EPSILON) stop 4 + if (abs(geolat_c(13,1)-(-35.2644)) > EPSILON) stop 6 + if (abs(geolat_c(1,13)-35.2644) > EPSILON) stop 8 + if (abs(geolon(5,5)-337.656) > EPSILON) stop 10 + if (abs(geolon(10,10)-17.9245) > EPSILON) stop 12 + if (abs(geolat(7,3)-(-27.84)) > EPSILON) stop 14 + if (abs(geolat(2,9)-16.8678) > EPSILON) stop 16 + if (abs(dx(1,1)-631596.076) > EPSILON) stop 18 + if (abs(dx(12,12)-631596.076) > EPSILON) stop 20 + +! There is no pole on this tile, so the pole variables +! should be .false. The routine sets dx equal to dy, +! so they should match. + + do j = 1, jm + do i = 1, im + if (abs(dx(i,j)-dy(i,j)) > EPSILON) stop 22 + if (is_north_pole(i,j)) stop 24 + if (is_south_pole(i,j)) stop 26 + enddo + enddo + + print*,"OK" + + print*,"SUCCESS" + + end program read_model_grid_file From f272d8c8c833dfb8d9ab0e81358f7dedb2368dd1 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Wed, 8 Jan 2025 13:14:04 -0600 Subject: [PATCH 45/67] Update routine 'read_mdl_grid_file' to use f90 version of netcdf. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 | 22 ++++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 index 14d4fc02d..a89d6e238 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 @@ -450,10 +450,11 @@ subroutine read_mdl_grid_file(mdl_grid_file, im, jm, & geolon, geolon_c, geolat, geolat_c, dx, dy, & is_north_pole, is_south_pole) + use netcdf + use orog_utils, only : find_poles, find_nearest_pole_points implicit none - include "netcdf.inc" character(len=*), intent(in) :: mdl_grid_file @@ -469,12 +470,11 @@ subroutine read_mdl_grid_file(mdl_grid_file, im, jm, & real, intent(out) :: dx(im,jm), dy(im,jm) integer :: i, j - integer :: ncid, error, fsize, id_var, nx, ny + integer :: ncid, error, id_var, nx, ny integer :: i_south_pole,j_south_pole integer :: i_north_pole,j_north_pole real, allocatable :: tmpvar(:,:) - fsize = 66536 nx = 2*im ny = 2*jm @@ -483,12 +483,12 @@ subroutine read_mdl_grid_file(mdl_grid_file, im, jm, & print*, "- OPEN AND READ= ", trim(mdl_grid_file) - error=NF__OPEN(mdl_grid_file,NF_NOWRITE,fsize,ncid) + error=nf90_open(mdl_grid_file, nf90_nowrite, ncid) call netcdf_err(error, 'Opening file '//trim(mdl_grid_file) ) - error=nf_inq_varid(ncid, 'x', id_var) + error=nf90_inq_varid(ncid, 'x', id_var) call netcdf_err(error, 'inquire varid of x from file ' // trim(mdl_grid_file)) - error=nf_get_var_double(ncid, id_var, tmpvar) + error=nf90_get_var(ncid, id_var, tmpvar) call netcdf_err(error, 'inquire data of x from file ' // trim(mdl_grid_file)) ! Adjust lontitude to be between 0 and 360. @@ -502,9 +502,9 @@ subroutine read_mdl_grid_file(mdl_grid_file, im, jm, & geolon(1:IM,1:JM) = tmpvar(2:nx:2,2:ny:2) geolon_c(1:IM+1,1:JM+1) = tmpvar(1:nx+1:2,1:ny+1:2) - error=nf_inq_varid(ncid, 'y', id_var) + error=nf90_inq_varid(ncid, 'y', id_var) call netcdf_err(error, 'inquire varid of y from file ' // trim(mdl_grid_file)) - error=nf_get_var_double(ncid, id_var, tmpvar) + error=nf90_get_var(ncid, id_var, tmpvar) call netcdf_err(error, 'inquire data of y from file ' // trim(mdl_grid_file)) geolat(1:IM,1:JM) = tmpvar(2:nx:2,2:ny:2) @@ -521,12 +521,12 @@ subroutine read_mdl_grid_file(mdl_grid_file, im, jm, & allocate(tmpvar(nx,ny)) - error=nf_inq_varid(ncid, 'area', id_var) + error=nf90_inq_varid(ncid, 'area', id_var) call netcdf_err(error, 'inquire varid of area from file ' // trim(mdl_grid_file)) - error=nf_get_var_double(ncid, id_var, tmpvar) + error=nf90_get_var(ncid, id_var, tmpvar) call netcdf_err(error, 'inquire data of area from file ' // trim(mdl_grid_file)) - error = nf_close(ncid) + error = nf90_close(ncid) do j = 1, jm do i = 1, im From 1bade53d19bb7c531e5bf16273684bf749d10b17 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Thu, 9 Jan 2025 10:28:06 -0600 Subject: [PATCH 46/67] Baseline unit test for routine 'read_mask'. Fix doxygen in that routine. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 | 6 +- tests/orog/CMakeLists.txt | 8 ++ tests/orog/data/C48.mx500.tile1.nc | Bin 0 -> 55384 bytes tests/orog/ftst_read_mask.F90 | 82 +++++++++++++++++++ 4 files changed, 93 insertions(+), 3 deletions(-) create mode 100644 tests/orog/data/C48.mx500.tile1.nc create mode 100644 tests/orog/ftst_read_mask.F90 diff --git a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 index a89d6e238..d6dc508cc 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 @@ -343,9 +343,9 @@ end subroutine write_mask_netcdf !> Read the land mask file !! !! @param[in] merge_file path -!! @param[in] slm Land-sea mask. -!! @param[in] land_frac Land fraction. -!! @param[in] lake_frac Lake fraction +!! @param[out] slm Land-sea mask. +!! @param[out] land_frac Land fraction. +!! @param[out] lake_frac Lake fraction !! @param[in] im 'i' dimension of a model grid tile. !! @param[in] jm 'j' dimension of a model grid tile. !! @author George Gayno NOAA/EMC diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index 7986143ea..ca02d8215 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -23,6 +23,10 @@ execute_process( COMMAND ${CMAKE_COMMAND} -E copy execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/data/C12_grid.tile1.nc ${CMAKE_CURRENT_BINARY_DIR}/C12_grid.tile1.nc) +# Note, the "read_mask" test expects this file. +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/C48.mx500.tile1.nc ${CMAKE_CURRENT_BINARY_DIR}/C48.mx500.tile1.nc) + add_executable(ftst_ll2xyz ftst_ll2xyz.F90) add_test(NAME orog-ftst_ll2xyz COMMAND ftst_ll2xyz) target_link_libraries(ftst_ll2xyz orog_lib) @@ -86,3 +90,7 @@ target_link_libraries(ftst_read_mdl_dims orog_lib) add_executable(ftst_read_mdl_grid_file ftst_read_mdl_grid_file.F90) add_test(NAME orog-ftst_read_mdl_grid_file COMMAND ftst_read_mdl_grid_file) target_link_libraries(ftst_read_mdl_grid_file orog_lib) + +add_executable(ftst_read_mask ftst_read_mask.F90) +add_test(NAME orog-ftst_read_mask COMMAND ftst_read_mask) +target_link_libraries(ftst_read_mask orog_lib) diff --git a/tests/orog/data/C48.mx500.tile1.nc b/tests/orog/data/C48.mx500.tile1.nc new file mode 100644 index 0000000000000000000000000000000000000000..51b5fc0a40eaee21145fa6b33a3dc215d82f0648 GIT binary patch literal 55384 zcmeI13vg7`8OQHtlMq6H@J7nxij0K;HX#HEA$u=jvmszeBtbxE37hPOEbJpD3mD`P zM@vPhP_T%V5zz`Z0y^a(uNIh5d{Y4jkSQ%Q7DfjZP#9}%rRUsxKlZ^s*xhq)vYqe3 z?sw1Y`~K&59(!}{s*SxNa%e%yV*1wo*`k->PZlB}ElU9gpae21=Rw zRHY2s1lsf@PDiU`X43C=diu?A#j@czXim}@U(ifFmqz|2b6lOh(c^At;<}N{56LHq zd}i3J{A7!Z=kxdxBlAa&%FoT?bFyYX ztap3YTz>auxEMa6u1y{*AQFf{UuYz>Wew~HTT$Tp1T z2bZ{=wtCAv8*k+6Jf8aetSncJ+t65RAPtA{Cb!k6WS-qqx3Jn!<95rR4lNjjYl1{$osI~}|(_YR-~LN1=`mB#g2LGDtdVb!P0!`l02QZct)He8Rs zQUiTEE^sCE4SXn6npwom7bf_ZK`4#nbX=q(hd@2`(i;-=D=s{_c1)rpSW?cT2m*_Tyy*nsTc@c7ysLNVKiFP7CJoGBd@0q5_ z9OVT$U@Onj@nMJ_+A{8SN;I`6m8wgYwsdl*ceq`9l}VY}TP{gBTWB2YrzQ{%Iebcng z>BhaICz{GjC5_Y>xuknym7(Q+R5HLP62Y z>DKLhs4sPbd}VirJG4ju2_OL^fCP{L5qBG-_8>ZdcSAjDl*_6+FBjxPgxqI2S8HX*Dmust3IRfNBDWpLuq(Ld92d5ohJgvfP))CsU zw}9maL4ShVnI-=D*-8>D(Dbi`DI|adkN^@u0!RP}Ab~_B5c6}R!@tMyya*Pb8~yoU z3JD+qB!C2v01`j~iB2Hq=f-vZJ%{H-u=w2Q&j(XT00|%gB!C2v01`-a0x>@~?(y$A zJTHR9=SF`%m_hk1tCRw1}lnjAn5~@lvsA^kA|5u?a?XVzBW# z;oQZsT}yHTSkpaJ{&YEetxf9x`I=3}f0v#Y#-5%k_P=BaIQR7Ge{6f((Z8eJo5mMs zye<57(-?8s*Gs&NuXy!e-BL=%!mEE(adN=euo?GgyllRIhj3|YsknR0^4pqyHU6jB zCOX}_g)!M?anF}4{ObI3zEhmrua{WyrBe`FCXS0HE~~nGyLR`}ohj{pS#dH`YUdAW5R##+pVMj+D$h4(|`G(SaPC+IC$w(ny%?nU$M~-Yvt4CYBHYc zULz`gDtO}PZv^L!OTvcZQ<=;k(0_5M`ueBZa{cykh|R(dY*Jr;%|q(x-?C7B{guV) z>7U)8zJ6n#diqDFyTtboY!wc?l_$89N2v0fHg>3b`hC~XZ}-!^js0aAppE|W-RkK7 zsJk}$5A9M%|Jn}P=)cgSj{ZOEwb6fQyE^)7|2vZT?^pFJeFe*oE`o98WMM|rGNDuM zMxl4n1XZ4ep8pqn{f}HJ2v`fJ@$5ZX z(m4J0Z^+h3IZU55tCy>6L9q1|ZQg;ePBPyQL-!+Eb)oAoKk=Gz zh5d%{qT|Pe)bFmRYtw(ds-F(cGujffJF)eD>r8C@iyeurzuT$@#D|`$6R&>a&~P0@ zJIf>2zj;0}_fP2AjI)WYzw>Jo1n2q7!Vm5X8vG6bHMOkTzXxjR2EU1|f6X6tLVjj% z;mVZ3D(y+YNnQJR|6X^c|JM3K@q^FIBDdW8&My2jh`PGLd)JS9-j)9Qzd0-YjeYTj z9OKP*HVWO1y~Uy<)-Y`$g1WZW`R+J=-_yq2vt5Of@y$Z&l_SE2(@El=RvN_Yl-Z(b zwfw&Ex%0AqUp3b=*Nvm^6zRJ~`i{{%9=^3&_x{Vk7CbH+yFb;hE0V$F$l!sD}E4*GYM{Qf!`kFnRCKcttq>EI&q*K;i5 z-};-xA3F15?$Yms`7IlSozJ`<^xRxYKle?B@tM=a<+K)m;* zgTk5Pe~|i9+al|K#ZBkDeCfBZHFnzYhEO=y6?9CX9&5XMczFHZHQ6jaZ&)T?t}=wC};uhm*N#EVxqvR=oUh1AVyh7@W-Mi~sffK*EE>bLiCsx- zEyHOERz{-)%ZI$!RK(s2;@+(dra$)Rp<4#i1GR&d#ghkqS=wNFz<;o^c=Et63mcD) zXxg&S8q7yD8se!NEH;b0cr-=RmW5XE5lviHbs;tj%Sux;epu)RAJN2RRTpBju&gvi z=11Euu)db*K EPSILON) stop 12 + if (abs(lake_frac(i,j)-lake_frac_chk(ij)) > EPSILON) stop 15 + if (abs(slm(i,j)-slm_chk(ij)) > EPSILON) stop 18 + enddo + + print*,"OK" + + print*,"SUCCESS" + + end program read_mask_test From f140e7a7b846f60963ae9eb3459abceee6c5ee04 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Thu, 9 Jan 2025 11:48:44 -0600 Subject: [PATCH 47/67] Convert 'read_mdl' routine to use f90 versions of netcdf library. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 | 23 ++++++++++---------- 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 index d6dc508cc..51f08e0db 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 @@ -352,8 +352,9 @@ end subroutine write_mask_netcdf subroutine read_mask(merge_file,slm,land_frac,lake_frac,im,jm) + use netcdf + implicit none - include "netcdf.inc" character(len=*), intent(in) :: merge_file @@ -363,30 +364,28 @@ subroutine read_mask(merge_file,slm,land_frac,lake_frac,im,jm) real, intent(out) :: lake_frac(im,jm) real, intent(out) :: slm(im,jm) - integer ncid, error, fsize, id_var - - fsize = 66536 + integer ncid, error, id_var print*,'- READ IN EXTERNAL LANDMASK FILE: ',trim(merge_file) - error=NF__OPEN(merge_file,NF_NOWRITE,fsize,ncid) + error=nf90_open(merge_file,nf90_nowrite,ncid) call netcdf_err(error, 'Open file '//trim(merge_file) ) - error=nf_inq_varid(ncid, 'land_frac', id_var) + error=nf90_inq_varid(ncid, 'land_frac', id_var) call netcdf_err(error, 'inquire varid of land_frac') - error=nf_get_var_double(ncid, id_var, land_frac) + error=nf90_get_var(ncid, id_var, land_frac) call netcdf_err(error, 'inquire data of land_frac') - error=nf_inq_varid(ncid, 'slmsk', id_var) + error=nf90_inq_varid(ncid, 'slmsk', id_var) call netcdf_err(error, 'inquire varid of slmsk') - error=nf_get_var_double(ncid, id_var, slm) + error=nf90_get_var(ncid, id_var, slm) call netcdf_err(error, 'inquire data of slmsk') - error=nf_inq_varid(ncid, 'lake_frac', id_var) + error=nf90_inq_varid(ncid, 'lake_frac', id_var) call netcdf_err(error, 'inquire varid of lake_frac') - error=nf_get_var_double(ncid, id_var, lake_frac) + error=nf90_get_var(ncid, id_var, lake_frac) call netcdf_err(error, 'inquire data of lake_frac') - error = nf_close(ncid) + error = nf90_close(ncid) end subroutine read_mask From 8a358381614a2f019ecedceeaf724dd7f9aa9d03 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Thu, 9 Jan 2025 15:17:56 -0600 Subject: [PATCH 48/67] Update 'qc_orog_by_ramp' to read the 'j' dimension of the ramp data from the file. This will make unit testing easier. Also, update the routine to use the f90 netcdf routines. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 | 33 ++++++++++++-------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 index 51f08e0db..576e628b9 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 @@ -664,42 +664,49 @@ end subroutine read_global_mask !! @author G. Gayno subroutine qc_orog_by_ramp(imn, jmn, zavg, zslm) - implicit none + use netcdf - include 'netcdf.inc' + implicit none integer, intent(in) :: imn, jmn integer, intent(inout) :: zavg(imn,jmn) integer, intent(inout) :: zslm(imn,jmn) - integer :: i, j, error, ncid, id_var, fsize + integer :: i, j, error, ncid, id_var, id_dim, jramp real(4), allocatable :: gice(:,:) - fsize = 65536 - - allocate (GICE(IMN+1,3601)) - ! Read 30-sec Antarctica RAMP data. Points scan from South ! to North, and from Greenwich to Greenwich. print*,"- OPEN/READ RAMP DATA ./topography.antarctica.ramp.30s.nc" - error=NF__OPEN("./topography.antarctica.ramp.30s.nc", & - NF_NOWRITE,fsize,ncid) + error=nf90_open("./topography.antarctica.ramp.30s.nc", & + nf90_nowrite, ncid) call netcdf_err(error, 'Opening RAMP topo file' ) - error=nf_inq_varid(ncid, 'topo', id_var) + + error=nf90_inq_dimid(ncid, 'jdim', id_dim) + call netcdf_err(error, 'Inquire dimid of jdim' ) + + error=nf90_inquire_dimension(ncid, id_dim, len=jramp) + call netcdf_err(error, 'Reading jdim' ) + + allocate (GICE(IMN+1,jramp)) + + error=nf90_inq_varid(ncid, 'topo', id_var) call netcdf_err(error, 'Inquire varid of RAMP topo') - error=nf_get_var_real(ncid, id_var, GICE) + + error=nf90_get_var(ncid, id_var, GICE) call netcdf_err(error, 'Inquire data of RAMP topo') - error = nf_close(ncid) + + error = nf90_close(ncid) print*,"- QC GLOBAL OROGRAPHY DATA WITH RAMP." ! If RAMP values are valid, replace the global value with the RAMP ! value. Invalid values are less than or equal to 0 (0, -1, or -99). - do j = 1, 3601 + do j = 1, jramp do i = 1, IMN if( GICE(i,j) .ne. -99. .and. GICE(i,j) .ne. -1.0 ) then if ( GICE(i,j) .gt. 0.) then From 2dc24d60a9fffe168a732b98612ac8af0cc2ae31 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 10 Jan 2025 16:12:55 +0000 Subject: [PATCH 49/67] Baseline test for routine "qc_orog_by_ramp". Fixes #1000. --- tests/orog/CMakeLists.txt | 8 ++ .../data/topography.antarctica.ramp.lowres.nc | Bin 0 -> 13633 bytes tests/orog/ftst_qc_orog_by_ramp.F90 | 97 ++++++++++++++++++ 3 files changed, 105 insertions(+) create mode 100644 tests/orog/data/topography.antarctica.ramp.lowres.nc create mode 100644 tests/orog/ftst_qc_orog_by_ramp.F90 diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index ca02d8215..584c65320 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -27,6 +27,10 @@ execute_process( COMMAND ${CMAKE_COMMAND} -E copy execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/data/C48.mx500.tile1.nc ${CMAKE_CURRENT_BINARY_DIR}/C48.mx500.tile1.nc) +# Note, the "qc_orog_by_ramp" test expects this file. +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/topography.antarctica.ramp.lowres.nc ${CMAKE_CURRENT_BINARY_DIR}/topography.antarctica.ramp.30s.nc) + add_executable(ftst_ll2xyz ftst_ll2xyz.F90) add_test(NAME orog-ftst_ll2xyz COMMAND ftst_ll2xyz) target_link_libraries(ftst_ll2xyz orog_lib) @@ -94,3 +98,7 @@ target_link_libraries(ftst_read_mdl_grid_file orog_lib) add_executable(ftst_read_mask ftst_read_mask.F90) add_test(NAME orog-ftst_read_mask COMMAND ftst_read_mask) target_link_libraries(ftst_read_mask orog_lib) + +add_executable(ftst_qc_orog_by_ramp ftst_qc_orog_by_ramp.F90) +add_test(NAME orog-ftst_qc_orog_by_ramp COMMAND ftst_qc_orog_by_ramp) +target_link_libraries(ftst_qc_orog_by_ramp orog_lib) diff --git a/tests/orog/data/topography.antarctica.ramp.lowres.nc b/tests/orog/data/topography.antarctica.ramp.lowres.nc new file mode 100644 index 0000000000000000000000000000000000000000..1023adf6e562af9431c42ebcfedb8a023ead7a22 GIT binary patch literal 13633 zcmeHNdu)@}6+ez05)u;^vY1FAwxp=@+dHmxmmmFRO~bXukkt(sa@?5}B6I<=ct@y|pz)UvgMopbL!Nr;m* z$wD@1-%0$P``!0FzjN=o-?{f-n?JCmY-L$dX(=#WQ{-)wG4(3p)5i|)>~8aiik5ha zifToqd9g~1SINLfimzHG-^+@?M4`$`dLA&-6QETxux_A$8Laqc0qD)doV|HHmyGF9 ziT%MSBn^u2Nwrva5;sE$lu}wgJCuth;3iDjhCj*`^7_4@us5Q4yCU9@FVgPQI=#VQ zdsnL#40Yex;)`f?A#Z1JgBIxycDIJS!M3|Kzc=EwZl(&%P%f(;%w_LS#Pnn~BdTAj zs?R0*hf>j;mWt||Q`wC5Gvoq@ zYp8E%+^%I}59GDm;#xAEOxqkb%%EO+2l373aM1hxl-S^My6DlwF_eQNw`qEIFx#Jt z4h{_4qZvJ#!*0i-_FObQXm{4-g_=^Rc?feA z3%xRzNxeF@P8NJZmA`qFC}1uI*5%3Ex{0!yu}eY}z7B6VEQCRNN}A+v?`-J`w|95- zhJD_S7Nx*%D+L~)(nwSdsZE7a;I4!oi}$(w$#gPapu(a=g?L`~aOQ(6x_xb*=Uy+G zw%2mB>qOrH!d#1(VbFf^!@p5A)_deEDOMS_ir*# zcQFtSw?MNpR3@tlE>pJkzr6BDLH21PZa8X=dN@^GQeB8zx#|(~nWA=e8r0@}x_w+K zdpIZ5#40-lYE`#C(J&1Zvx~m*z&YUT`}_S8ILBw3IL(c#HIWBg<>0!85-6@9<7eKe z1c_2x>+Ka!$Z(%khCk;q(A;@>$nTBT*6(Ni>~60P+C)+N+n&1;A`BGY0@a z`u^F0LIBysF`pZN7ZiX&C4vP3kVGjFSg+Q_tZH7yE>kPP6afBiE~l z74obP9MAQ$VsJdy!*bZA%sIas+>n4F0Yd_Y1Plon5-=oSNWhSQApt`Ih6D@=7!sJd z1jtp*2c7ICRjxjk2~!xY93hVE*~H&z~(5&hw7>JmVzah3&ULyFArcSExZ) zQ`}-bovMaf0xhAIE?-NpI91((Rgi-u9@V4ygg$<}I(D50N-14>J+vuLU03pPifZx7 zraycz|6?V8OkYv@Pg9MxPPI%9bu)Zt0dYP0$M(^wdPM$3q8l&$b#0 zTZ=y$p1L);qrryv>+m+<4wrqKy}>pR?`z!Qusd9KXZ&-w9}cwnn}Ki5Qp>MgZi<@N zMF|$V`N@t*lZ87|LeY5rlUFL=tZ}v*Z;D+PtdDw(Z8y% zD}KZ+>p!91PpWo)tm03p_MTViFQ|53Qu$s{?f=yMvi6Lpzvq-^|HIFDp!tU$`0F=4 zaQu)5_6~VK>+rytT2E^5n)}(|tL|t2Ug@cPsn)Zqvf2~wEAqVE`JQ{%lRt1b1)Xls N>(zJ2G5mL){{mapfdK#j literal 0 HcmV?d00001 diff --git a/tests/orog/ftst_qc_orog_by_ramp.F90 b/tests/orog/ftst_qc_orog_by_ramp.F90 new file mode 100644 index 000000000..b15145120 --- /dev/null +++ b/tests/orog/ftst_qc_orog_by_ramp.F90 @@ -0,0 +1,97 @@ + program qc_orog_ramp + +! Test routine "qc_orog_by_ramp", which adjusts +! the global terrain in the vicinity of Antarctica +! using 'RAMP' data. +! +! In OPS, the global data is 30-sec with dimensions +! 43200 x 21600. The RAMP data is 30-sec with dimension +! 43201 x 3601. For this test, reduced versions of +! both grids are used: global - 9x7; RAMP 10x5. +! +! Author George Gayno NCEP/EMC + + use io_utils, only : qc_orog_by_ramp + + implicit none + +! Dimensions of the global grid. + + integer, parameter :: imn = 9 + integer, parameter :: jmn = 7 + + integer :: i, j + +! The terrain height (zavg) and land mask (zslm) +! of the global grid. + + integer :: zavg(imn,jmn) + integer :: zslm(imn,jmn) + +! The expected values for a successful test. + + integer :: zavg_expected(imn,jmn) + integer :: zslm_expected(imn,jmn) + + print*,'- Starting test of qc_orog_by_ramp.' + +! Initialize the global grid to all ocean. + + zslm = 0 ! water mask + zavg = 0 ! sea level + +! These global grid points are outside the 'ramp' grid, +! so they should not change. + + zslm_expected(:,6:7) = 0 + zavg_expected(:,6:7) = 0 + +! These global grid points are located within the 'ramp' +! grid. For this test, the first two rows of the 'ramp' +! data have non-zero terrain. So, rows 1 and 2 of the global +! grid will become land, located above sea level. Rows +! 3,4 and 5 of the RAMP data are ocean. So, rows 3,4 and +! 5 of the global grid will remain ocean. + + zslm_expected(:,3:5) = 0 ! ocean mask + zavg_expected(:,3:5) = 0 ! terrain height + + zslm_expected(:,1:2) = 1 ! becomes land + + zavg_expected(1,1) = 5 ! acquires non-zero terrain. + zavg_expected(2,1) = 5 + zavg_expected(3,1) = 5 + zavg_expected(4,1) = 5 + zavg_expected(5,1) = 5 + zavg_expected(6,1) = 4 + zavg_expected(7,1) = 4 + zavg_expected(8,1) = 3 + zavg_expected(9,1) = 3 + + zavg_expected(1,2) = 2 + zavg_expected(2,2) = 2 + zavg_expected(3,2) = 3 + zavg_expected(4,2) = 2 + zavg_expected(5,2) = 2 + zavg_expected(6,2) = 2 + zavg_expected(7,2) = 1 + zavg_expected(8,2) = 1 + zavg_expected(9,2) = 0 ! Note: this 'ramp' point has non-zero terrain of + ! 0.14, which rounds down to zero. + +! Note: The location of the RAMP data is set in the routine. + + call qc_orog_by_ramp(imn, jmn, zavg, zslm) + + do i = 1, imn + do j = 1, jmn + if (zavg(i,j) /= zavg_expected(i,j)) stop 4 + if (zslm(i,j) /= zslm_expected(i,j)) stop 8 + enddo + enddo + + print*,"OK" + + print*,"SUCCESS" + + end program qc_orog_ramp From 92b232a2a187b3ddb8736c55493757cbd93d75ba Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 10 Jan 2025 16:46:09 +0000 Subject: [PATCH 50/67] Remove unneeded 'if' statement from 'qc_orog_by_ramp'. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 index 576e628b9..abb026c4a 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 @@ -708,11 +708,9 @@ subroutine qc_orog_by_ramp(imn, jmn, zavg, zslm) do j = 1, jramp do i = 1, IMN - if( GICE(i,j) .ne. -99. .and. GICE(i,j) .ne. -1.0 ) then - if ( GICE(i,j) .gt. 0.) then - ZAVG(i,j) = int( GICE(i,j) + 0.5 ) - ZSLM(i,j) = 1 - endif + if ( GICE(i,j) .gt. 0.) then + ZAVG(i,j) = int( GICE(i,j) + 0.5 ) + ZSLM(i,j) = 1 endif enddo enddo From 8a15fc3650f61de49dc6f189672e2ea4acc49545 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 10 Jan 2025 21:09:53 +0000 Subject: [PATCH 51/67] Update routine "write_mask_netcdf" to use f90 netcdf routines. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 | 41 ++++++++++---------- 1 file changed, 21 insertions(+), 20 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 index abb026c4a..6e3fb904a 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 @@ -262,6 +262,7 @@ end subroutine netcdf_err !! @author George Gayno NOAA/EMC subroutine write_mask_netcdf(im, jm, slm, land_frac, ntiles, tile, geolon, geolat) + use netcdf implicit none integer, intent(in):: im, jm, ntiles, tile real, intent(in), dimension(im,jm) :: slm, geolon, geolat, land_frac @@ -273,7 +274,6 @@ subroutine write_mask_netcdf(im, jm, slm, land_frac, ntiles, tile, geolon, geola integer :: dim_lon, dim_lat integer :: id_geolon,id_geolat integer :: id_slmsk,id_land_frac - include "netcdf.inc" if(ntiles > 1) then write(outfile, '(a,i4.4,a)') 'out.oro.tile', tile, '.nc' @@ -285,57 +285,58 @@ subroutine write_mask_netcdf(im, jm, slm, land_frac, ntiles, tile, geolon, geola dim2=jm !--- open the file - error = NF__CREATE(outfile, IOR(NF_NETCDF4,NF_CLASSIC_MODEL), inital, fsize, ncid) + error = nf90_create(outfile, IOR(NF90_NETCDF4,NF90_CLASSIC_MODEL), ncid, & + initialsize=inital, chunksize=fsize) call netcdf_err(error, 'Creating file '//trim(outfile) ) !--- define dimension - error = nf_def_dim(ncid, 'lon', dim1, dim_lon) + error = nf90_def_dim(ncid, 'lon', dim1, dim_lon) call netcdf_err(error, 'define dimension lon for file='//trim(outfile) ) - error = nf_def_dim(ncid, 'lat', dim2, dim_lat) + error = nf90_def_dim(ncid, 'lat', dim2, dim_lat) call netcdf_err(error, 'define dimension lat for file='//trim(outfile) ) !--- define field !---geolon - error = nf_def_var(ncid, 'geolon', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_geolon) + error = nf90_def_var(ncid, 'geolon', NF90_FLOAT, (/dim_lon,dim_lat/), id_geolon) call netcdf_err(error, 'define var geolon for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolon, "long_name", 9, "Longitude") + error = nf90_put_att(ncid, id_geolon, "long_name", "Longitude") call netcdf_err(error, 'define geolon name for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolon, "units", 12, "degrees_east") + error = nf90_put_att(ncid, id_geolon, "units", "degrees_east") call netcdf_err(error, 'define geolon units for file='//trim(outfile) ) !---geolat - error = nf_def_var(ncid, 'geolat', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_geolat) + error = nf90_def_var(ncid, 'geolat', NF90_FLOAT, (/dim_lon,dim_lat/), id_geolat) call netcdf_err(error, 'define var geolat for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolat, "long_name", 8, "Latitude") + error = nf90_put_att(ncid, id_geolat, "long_name", "Latitude") call netcdf_err(error, 'define geolat name for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolat, "units", 13, "degrees_north") + error = nf90_put_att(ncid, id_geolat, "units", "degrees_north") call netcdf_err(error, 'define geolat units for file='//trim(outfile) ) !---slmsk - error = nf_def_var(ncid, 'slmsk', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_slmsk) + error = nf90_def_var(ncid, 'slmsk', NF90_FLOAT, (/dim_lon,dim_lat/), id_slmsk) call netcdf_err(error, 'define var slmsk for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_slmsk, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_slmsk, "coordinates", "geolon geolat") call netcdf_err(error, 'define slmsk coordinates for file='//trim(outfile) ) !--- land_frac - error = nf_def_var(ncid, 'land_frac', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_land_frac) + error = nf90_def_var(ncid, 'land_frac', NF90_FLOAT, (/dim_lon,dim_lat/), id_land_frac) call netcdf_err(error, 'define var land_frac for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_land_frac, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_land_frac, "coordinates", "geolon geolat") call netcdf_err(error, 'define land_frac coordinates for file='//trim(outfile) ) - error = nf__enddef(ncid, header_buffer_val,4,0,4) + error = nf90_enddef(ncid, header_buffer_val,4,0,4) call netcdf_err(error, 'end meta define for file='//trim(outfile) ) !--- write out data - error = nf_put_var_double( ncid, id_geolon, geolon(:dim1,:dim2)) + error = nf90_put_var( ncid, id_geolon, geolon(:dim1,:dim2)) call netcdf_err(error, 'write var geolon for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_geolat, geolat(:dim1,:dim2)) + error = nf90_put_var( ncid, id_geolat, geolat(:dim1,:dim2)) call netcdf_err(error, 'write var geolat for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_slmsk, slm(:dim1,:dim2)) + error = nf90_put_var( ncid, id_slmsk, slm(:dim1,:dim2)) call netcdf_err(error, 'write var slmsk for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_land_frac, land_frac(:dim1,:dim2)) + error = nf90_put_var( ncid, id_land_frac, land_frac(:dim1,:dim2)) call netcdf_err(error, 'write var land_frac for file='//trim(outfile) ) - error = nf_close(ncid) + error = nf90_close(ncid) call netcdf_err(error, 'close file='//trim(outfile) ) end subroutine write_mask_netcdf From 99a2e4f74115e0fc7a381ec931b1b2b743574854 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 10 Jan 2025 21:52:12 +0000 Subject: [PATCH 52/67] Convert routine 'write_netcdf' to use f90 netcdf routines. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 | 137 ++++++++++--------- 1 file changed, 69 insertions(+), 68 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 index 6e3fb904a..f5ef6ff04 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 @@ -39,6 +39,7 @@ module io_utils !! @param[in] lat Latitude of the first column of the model grid tile. !! @author Jordan Alpert NOAA/EMC GFDL Programmer subroutine write_netcdf(im, jm, slm, land_frac, oro, hprime, ntiles, tile, geolon, geolat, lon, lat) + use netcdf implicit none integer, intent(in):: im, jm, ntiles, tile real, intent(in) :: lon(im), lat(jm) @@ -56,7 +57,6 @@ subroutine write_netcdf(im, jm, slm, land_frac, oro, hprime, ntiles, tile, geolo integer :: id_oa1,id_oa2,id_oa3,id_oa4 integer :: id_ol1,id_ol2,id_ol3,id_ol4 integer :: id_theta,id_gamma,id_sigma,id_elvmax - include "netcdf.inc" if(ntiles > 1) then write(outfile, '(a,i4.4,a)') 'out.oro.tile', tile, '.nc' @@ -68,164 +68,165 @@ subroutine write_netcdf(im, jm, slm, land_frac, oro, hprime, ntiles, tile, geolo dim2=size(lat,1) !--- open the file - error = NF__CREATE(outfile, IOR(NF_NETCDF4,NF_CLASSIC_MODEL), inital, fsize, ncid) + error = nf90_create(outfile, IOR(NF90_NETCDF4,NF90_CLASSIC_MODEL), ncid, & + initialsize=inital, chunksize=fsize) call netcdf_err(error, 'Creating file '//trim(outfile) ) !--- define dimension - error = nf_def_dim(ncid, 'lon', dim1, dim_lon) + error = nf90_def_dim(ncid, 'lon', dim1, dim_lon) call netcdf_err(error, 'define dimension lon for file='//trim(outfile) ) - error = nf_def_dim(ncid, 'lat', dim2, dim_lat) + error = nf90_def_dim(ncid, 'lat', dim2, dim_lat) call netcdf_err(error, 'define dimension lat for file='//trim(outfile) ) !--- define field !---geolon - error = nf_def_var(ncid, 'geolon', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_geolon) + error = nf90_def_var(ncid, 'geolon', NF90_FLOAT, (/dim_lon,dim_lat/), id_geolon) call netcdf_err(error, 'define var geolon for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolon, "long_name", 9, "Longitude") + error = nf90_put_att(ncid, id_geolon, "long_name", "Longitude") call netcdf_err(error, 'define geolon name for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolon, "units", 12, "degrees_east") + error = nf90_put_att(ncid, id_geolon, "units", "degrees_east") call netcdf_err(error, 'define geolon units for file='//trim(outfile) ) !---geolat - error = nf_def_var(ncid, 'geolat', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_geolat) + error = nf90_def_var(ncid, 'geolat', NF90_FLOAT, (/dim_lon,dim_lat/), id_geolat) call netcdf_err(error, 'define var geolat for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolat, "long_name", 8, "Latitude") + error = nf90_put_att(ncid, id_geolat, "long_name", "Latitude") call netcdf_err(error, 'define geolat name for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_geolat, "units", 13, "degrees_north") + error = nf90_put_att(ncid, id_geolat, "units", "degrees_north") call netcdf_err(error, 'define geolat units for file='//trim(outfile) ) !---slmsk - error = nf_def_var(ncid, 'slmsk', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_slmsk) + error = nf90_def_var(ncid, 'slmsk', NF90_FLOAT,(/dim_lon,dim_lat/), id_slmsk) call netcdf_err(error, 'define var slmsk for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_slmsk, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_slmsk, "coordinates", "geolon geolat") call netcdf_err(error, 'define slmsk coordinates for file='//trim(outfile) ) !--- land_frac - error = nf_def_var(ncid, 'land_frac', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_land_frac) + error = nf90_def_var(ncid, 'land_frac', NF90_FLOAT, (/dim_lon,dim_lat/), id_land_frac) call netcdf_err(error, 'define var land_frac for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_land_frac, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_land_frac, "coordinates", "geolon geolat") call netcdf_err(error, 'define land_frac coordinates for file='//trim(outfile) ) !---orography - raw - error = nf_def_var(ncid, 'orog_raw', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_orog_raw) + error = nf90_def_var(ncid, 'orog_raw', NF90_FLOAT, (/dim_lon,dim_lat/), id_orog_raw) call netcdf_err(error, 'define var orog_raw for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_orog_raw, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_orog_raw, "coordinates", "geolon geolat") call netcdf_err(error, 'define orog_raw coordinates for file='//trim(outfile) ) !---orography - filtered - error = nf_def_var(ncid, 'orog_filt', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_orog_filt) + error = nf90_def_var(ncid, 'orog_filt', NF90_FLOAT, (/dim_lon,dim_lat/), id_orog_filt) call netcdf_err(error, 'define var orog_filt for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_orog_filt, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_orog_filt, "coordinates", "geolon geolat") call netcdf_err(error, 'define orog_filt coordinates for file='//trim(outfile) ) !---stddev - error = nf_def_var(ncid, 'stddev', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_stddev) + error = nf90_def_var(ncid, 'stddev', NF90_FLOAT, (/dim_lon,dim_lat/), id_stddev) call netcdf_err(error, 'define var stddev for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_stddev, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_stddev, "coordinates", "geolon geolat") call netcdf_err(error, 'define stddev coordinates for file='//trim(outfile) ) !---convexity - error = nf_def_var(ncid, 'convexity', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_convex) + error = nf90_def_var(ncid, 'convexity', NF90_FLOAT, (/dim_lon,dim_lat/), id_convex) call netcdf_err(error, 'define var convexity for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_convex, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_convex, "coordinates", "geolon geolat") call netcdf_err(error, 'define convexity coordinates for file='//trim(outfile) ) !---oa1 -> oa4 - error = nf_def_var(ncid, 'oa1', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_oa1) + error = nf90_def_var(ncid, 'oa1', NF90_FLOAT, (/dim_lon,dim_lat/), id_oa1) call netcdf_err(error, 'define var oa1 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_oa1, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_oa1, "coordinates", "geolon geolat") call netcdf_err(error, 'define oa1 coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'oa2', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_oa2) + error = nf90_def_var(ncid, 'oa2', NF90_FLOAT, (/dim_lon,dim_lat/), id_oa2) call netcdf_err(error, 'define var oa2 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_oa2, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_oa2, "coordinates", "geolon geolat") call netcdf_err(error, 'define oa2 coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'oa3', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_oa3) + error = nf90_def_var(ncid, 'oa3', NF90_FLOAT, (/dim_lon,dim_lat/), id_oa3) call netcdf_err(error, 'define var oa3 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_oa3, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_oa3, "coordinates", "geolon geolat") call netcdf_err(error, 'define oa3 coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'oa4', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_oa4) + error = nf90_def_var(ncid, 'oa4', NF90_FLOAT, (/dim_lon,dim_lat/), id_oa4) call netcdf_err(error, 'define var oa4 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_oa4, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_oa4, "coordinates", "geolon geolat") call netcdf_err(error, 'define oa4 coordinates for file='//trim(outfile) ) !---ol1 -> ol4 - error = nf_def_var(ncid, 'ol1', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_ol1) + error = nf90_def_var(ncid, 'ol1', NF90_FLOAT, (/dim_lon,dim_lat/), id_ol1) call netcdf_err(error, 'define var ol1 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_ol1, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_ol1, "coordinates", "geolon geolat") call netcdf_err(error, 'define ol1 coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'ol2', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_ol2) + error = nf90_def_var(ncid, 'ol2', NF90_FLOAT, (/dim_lon,dim_lat/), id_ol2) call netcdf_err(error, 'define var ol2 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_ol2, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_ol2, "coordinates", "geolon geolat") call netcdf_err(error, 'define ol2 coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'ol3', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_ol3) + error = nf90_def_var(ncid, 'ol3', NF90_FLOAT, (/dim_lon,dim_lat/), id_ol3) call netcdf_err(error, 'define var ol3 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_ol3, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_ol3, "coordinates", "geolon geolat") call netcdf_err(error, 'define ol3 coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'ol4', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_ol4) + error = nf90_def_var(ncid, 'ol4', NF90_FLOAT, (/dim_lon,dim_lat/), id_ol4) call netcdf_err(error, 'define var ol4 for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_ol4, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_ol4, "coordinates", "geolon geolat") call netcdf_err(error, 'define ol4 coordinates for file='//trim(outfile) ) !---theta gamma sigma elvmax - error = nf_def_var(ncid, 'theta', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_theta) + error = nf90_def_var(ncid, 'theta', NF90_FLOAT, (/dim_lon,dim_lat/), id_theta) call netcdf_err(error, 'define var theta for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_theta, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_theta, "coordinates", "geolon geolat") call netcdf_err(error, 'define theta coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'gamma', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_gamma) + error = nf90_def_var(ncid, 'gamma', NF90_FLOAT, (/dim_lon,dim_lat/), id_gamma) call netcdf_err(error, 'define var gamma for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_gamma, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_gamma, "coordinates", "geolon geolat") call netcdf_err(error, 'define gamma coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'sigma', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_sigma) + error = nf90_def_var(ncid, 'sigma', NF90_FLOAT, (/dim_lon,dim_lat/), id_sigma) call netcdf_err(error, 'define var sigma for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_sigma, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_sigma, "coordinates", "geolon geolat") call netcdf_err(error, 'define sigma coordinates for file='//trim(outfile) ) - error = nf_def_var(ncid, 'elvmax', NF_FLOAT, 2, (/dim_lon,dim_lat/), id_elvmax) + error = nf90_def_var(ncid, 'elvmax', NF90_FLOAT, (/dim_lon,dim_lat/), id_elvmax) call netcdf_err(error, 'define var elvmax for file='//trim(outfile) ) - error = nf_put_att_text(ncid, id_elvmax, "coordinates", 13, "geolon geolat") + error = nf90_put_att(ncid, id_elvmax, "coordinates", "geolon geolat") call netcdf_err(error, 'define elvmax coordinates for file='//trim(outfile) ) - error = nf__enddef(ncid, header_buffer_val,4,0,4) + error = nf90_enddef(ncid, header_buffer_val,4,0,4) call netcdf_err(error, 'end meta define for file='//trim(outfile) ) !--- write out data - error = nf_put_var_double( ncid, id_geolon, geolon(:dim1,:dim2)) + error = nf90_put_var( ncid, id_geolon, geolon(:dim1,:dim2)) call netcdf_err(error, 'write var geolon for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_geolat, geolat(:dim1,:dim2)) + error = nf90_put_var( ncid, id_geolat, geolat(:dim1,:dim2)) call netcdf_err(error, 'write var geolat for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_slmsk, slm(:dim1,:dim2)) + error = nf90_put_var( ncid, id_slmsk, slm(:dim1,:dim2)) call netcdf_err(error, 'write var slmsk for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_land_frac, land_frac(:dim1,:dim2)) + error = nf90_put_var( ncid, id_land_frac, land_frac(:dim1,:dim2)) call netcdf_err(error, 'write var land_frac for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_orog_raw, oro(:dim1,:dim2)) + error = nf90_put_var( ncid, id_orog_raw, oro(:dim1,:dim2)) call netcdf_err(error, 'write var orog_raw for file='//trim(outfile) ) ! We no longer filter the orog, so the raw and filtered records are the same. - error = nf_put_var_double( ncid, id_orog_filt, oro(:dim1,:dim2)) + error = nf90_put_var( ncid, id_orog_filt, oro(:dim1,:dim2)) call netcdf_err(error, 'write var orog_filt for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_stddev, hprime(:dim1,:dim2,1)) + error = nf90_put_var( ncid, id_stddev, hprime(:dim1,:dim2,1)) call netcdf_err(error, 'write var stddev for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_convex, hprime(:dim1,:dim2,2)) + error = nf90_put_var( ncid, id_convex, hprime(:dim1,:dim2,2)) call netcdf_err(error, 'write var convex for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_oa1, hprime(:dim1,:dim2,3)) + error = nf90_put_var( ncid, id_oa1, hprime(:dim1,:dim2,3)) call netcdf_err(error, 'write var oa1 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_oa2, hprime(:dim1,:dim2,4)) + error = nf90_put_var( ncid, id_oa2, hprime(:dim1,:dim2,4)) call netcdf_err(error, 'write var oa2 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_oa3, hprime(:dim1,:dim2,5)) + error = nf90_put_var( ncid, id_oa3, hprime(:dim1,:dim2,5)) call netcdf_err(error, 'write var oa3 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_oa4, hprime(:dim1,:dim2,6)) + error = nf90_put_var( ncid, id_oa4, hprime(:dim1,:dim2,6)) call netcdf_err(error, 'write var oa4 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_ol1, hprime(:dim1,:dim2,7)) + error = nf90_put_var( ncid, id_ol1, hprime(:dim1,:dim2,7)) call netcdf_err(error, 'write var ol1 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_ol2, hprime(:dim1,:dim2,8)) + error = nf90_put_var( ncid, id_ol2, hprime(:dim1,:dim2,8)) call netcdf_err(error, 'write var ol2 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_ol3, hprime(:dim1,:dim2,9)) + error = nf90_put_var( ncid, id_ol3, hprime(:dim1,:dim2,9)) call netcdf_err(error, 'write var ol3 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_ol4, hprime(:dim1,:dim2,10)) + error = nf90_put_var( ncid, id_ol4, hprime(:dim1,:dim2,10)) call netcdf_err(error, 'write var ol4 for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_theta, hprime(:dim1,:dim2,11)) + error = nf90_put_var( ncid, id_theta, hprime(:dim1,:dim2,11)) call netcdf_err(error, 'write var theta for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_gamma, hprime(:dim1,:dim2,12)) + error = nf90_put_var( ncid, id_gamma, hprime(:dim1,:dim2,12)) call netcdf_err(error, 'write var gamma for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_sigma, hprime(:dim1,:dim2,13)) + error = nf90_put_var( ncid, id_sigma, hprime(:dim1,:dim2,13)) call netcdf_err(error, 'write var sigma for file='//trim(outfile) ) - error = nf_put_var_double( ncid, id_elvmax, hprime(:dim1,:dim2,14)) + error = nf90_put_var( ncid, id_elvmax, hprime(:dim1,:dim2,14)) call netcdf_err(error, 'write var elvmax for file='//trim(outfile) ) - error = nf_close(ncid) + error = nf90_close(ncid) call netcdf_err(error, 'close file='//trim(outfile) ) end subroutine write_netcdf From 1e57a0f5d570a22fb1fff29f007796681fa3ea4b Mon Sep 17 00:00:00 2001 From: George Gayno Date: Fri, 10 Jan 2025 22:04:36 +0000 Subject: [PATCH 53/67] Convert netcdf error routine to use f90 versions of netcdf. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 index f5ef6ff04..7ccb86765 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/io_utils.F90 @@ -237,18 +237,19 @@ end subroutine write_netcdf !! @param[in] string The NetCDF error message !! @author Jordan Alpert NOAA/EMC subroutine netcdf_err( err, string ) + use netcdf + implicit none integer, intent(in) :: err character(len=*), intent(in) :: string character(len=256) :: errmsg - include "netcdf.inc" - if( err.EQ.NF_NOERR )return - errmsg = NF_STRERROR(err) + if( err.EQ.NF90_NOERR )return + errmsg = NF90_STRERROR(err) print*, 'FATAL ERROR: ', trim(string), ': ', trim(errmsg) call abort return - end subroutine netcdf_err + end subroutine netcdf_err !> Write the land mask file !! From 6a68569815c011acc9e9fe8df3d1303833ec688f Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 13 Jan 2025 17:07:26 +0000 Subject: [PATCH 54/67] Remove some diagnostic print from the orog code. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 index 0eb10c09b..474e07972 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 @@ -547,14 +547,12 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, & real, intent(out) :: oro(im,jm) real, intent(out) :: var(im,jm),var4(im,jm) -! integer, parameter :: MAXSUM=65000000 - integer :: MAXSUM real, parameter :: D2R = 3.14159265358979/180. real, dimension(:), allocatable :: hgt_1d, hgt_1d_all real GLAT(JMN), GLON(IMN) - integer JST, JEN + integer JST, JEN, maxsum real LONO(4),LATO(4),LONI,LATI real LONO_RAD(4), LATO_RAD(4) real HEIGHT @@ -590,9 +588,9 @@ SUBROUTINE MAKEMT2(ZAVG,ZSLM,ORO,SLM,VAR,VAR4, & LATO(4) = lat_c(i,j+1) call get_index(IMN,JMN,4,LONO,LATO,DELXN,jst,jen,ilist,numx) MAXSUM=MAX(MAXSUM,(JEN-JST+1)*IMN) - print*,'test point ',i,j,jst,jen,imn,((JEN-JST+1)*IMN),maxsum ENDDO ENDDO + print*,"- MAXSUM IS ", maxsum allocate(hgt_1d(MAXSUM)) allocate(hgt_1d_all(MAXSUM)) ! From 8acb8a80bb99636a34b6b5951656634787f7e457 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 13 Jan 2025 18:58:04 +0000 Subject: [PATCH 55/67] Revert cpld_gridgen changes. These will be handled by Denise. Fixes #1000. --- sorc/cpld_gridgen.fd/grdvars.F90 | 5 ++--- ush/cpld_gridgen.sh | 2 +- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/sorc/cpld_gridgen.fd/grdvars.F90 b/sorc/cpld_gridgen.fd/grdvars.F90 index 431b76a4f..a3abad291 100644 --- a/sorc/cpld_gridgen.fd/grdvars.F90 +++ b/sorc/cpld_gridgen.fd/grdvars.F90 @@ -163,9 +163,8 @@ module grdvars !! rounded to minimum_depth real(kind=real_kind), parameter :: maximum_lat = 88.0 !< The maximum latitude for water points for WW3 - integer, parameter :: nar = 3 !< the number of possible ATM resolutions - integer, parameter, dimension(nar) :: catm = (/12, 18, 24/) !< the ATM resolutions for mapped ocean masks - + integer, parameter :: nar = 6 !< the number of possible ATM resolutions + integer, parameter, dimension(nar) :: catm = (/48, 96, 192, 384, 768, 1152/) !< the ATM resolutions for mapped ocean masks contains !> Allocate grid variables diff --git a/ush/cpld_gridgen.sh b/ush/cpld_gridgen.sh index 8155a3657..9bb08109a 100755 --- a/ush/cpld_gridgen.sh +++ b/ush/cpld_gridgen.sh @@ -20,7 +20,7 @@ export RESNAME=${RESNAME:-$1} export DEBUG=.false. export MASKEDIT=.false. export DO_POSTWGTS=.true. -export MOSAICDIR_PATH=${MOSAICDIR_PATH:-$PATHTR/fix/orog.lowres} +export MOSAICDIR_PATH=${MOSAICDIR_PATH:-$PATHTR/fix/orog} export FIXDIR_PATH=${MOM6_FIXDIR}/${RESNAME} APRUN=${APRUN:-"srun"} From a70379df2514f647f244c36c677ad64e61b0becf Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 13 Jan 2025 20:28:59 +0000 Subject: [PATCH 56/67] Minor cleanup to ./tests/orog/CMakeLists.txt. Fixes #1000. --- tests/orog/CMakeLists.txt | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/tests/orog/CMakeLists.txt b/tests/orog/CMakeLists.txt index 584c65320..7c2d894e1 100644 --- a/tests/orog/CMakeLists.txt +++ b/tests/orog/CMakeLists.txt @@ -3,31 +3,34 @@ # George Gayno, Ed Hartnett if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel|IntelLLVM)$") - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8 -warn unused") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8") elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffree-line-length-0 -fdefault-real-8") endif() # Copy necessary test files from the source data directory to the -# build data directory. +# build test directory. -# Note, the "read_global_mask" routine expects the filename "landcover.umd.30s.nc". +# The "read_global_mask" test needs this file. Note that the copy command changes the +# file name to "landcover.umd.30s.nc", which is the name expected by the "read_global_mask" routine. execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/data/landcover.umd.lowres.nc ${CMAKE_CURRENT_BINARY_DIR}/landcover.umd.30s.nc) -# Note, the "read_global_mask" routine expects the filename "topography.gmted2010.30s.nc". +# The "read_global_orog" test needs this file. Note that the copy command changes the +# file name to "topography.gmted2010.30s.nc", which is the name expected by the "read_global_orog" routine. execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/data/topography.gmted2010.lowres.nc ${CMAKE_CURRENT_BINARY_DIR}/topography.gmted2010.30s.nc) -# Note, the "read_mdl_dims" and "read_mdl_grid_file" tests expects this file +# The "read_mdl_dims" and "read_mdl_grid_file" tests expects this file. execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/data/C12_grid.tile1.nc ${CMAKE_CURRENT_BINARY_DIR}/C12_grid.tile1.nc) -# Note, the "read_mask" test expects this file. +# The "read_mask" test expects this file. execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/data/C48.mx500.tile1.nc ${CMAKE_CURRENT_BINARY_DIR}/C48.mx500.tile1.nc) -# Note, the "qc_orog_by_ramp" test expects this file. +# The "qc_orog_by_ramp" test needs this file. Note that the copy command changes the +# file name to "topography.antarctica.ramp.30s.nc", which is the name expected by the "qc_orog_by_ramp" routine. execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/data/topography.antarctica.ramp.lowres.nc ${CMAKE_CURRENT_BINARY_DIR}/topography.antarctica.ramp.30s.nc) From bc2dac8e9ec173060388ba2e48dcd5cc5934c413 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 13 Jan 2025 20:58:25 +0000 Subject: [PATCH 57/67] Add some comments to ./tests/orog/ftst_find_poles.F90. Fixes #1000. --- tests/orog/ftst_find_poles.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/orog/ftst_find_poles.F90 b/tests/orog/ftst_find_poles.F90 index 3d98cb972..b0b9e953f 100644 --- a/tests/orog/ftst_find_poles.F90 +++ b/tests/orog/ftst_find_poles.F90 @@ -2,6 +2,9 @@ program ftst_find_poles ! Unit test for routine find_poles. ! +! Checks if a model tile contains a 'pole' +! point based on the latitude of each point. +! ! Author George Gayno NCEP/EMC use orog_utils, only : find_poles From 9147bf385ef0a002ebc9305c9f2bce377ec0cd33 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 13 Jan 2025 21:23:11 +0000 Subject: [PATCH 58/67] Update comments in ./tests/orog/ftst_get_index.F90. Fixes #1000. --- tests/orog/ftst_get_index.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/orog/ftst_get_index.F90 b/tests/orog/ftst_get_index.F90 index 8831f739b..3ebf86d42 100644 --- a/tests/orog/ftst_get_index.F90 +++ b/tests/orog/ftst_get_index.F90 @@ -1,6 +1,6 @@ program test_get_index -! Unit test for routine get_index, which finds the +! Unit test for routine "get_index", which finds the ! i/j location of the model point on the high-resolution ! mask and orography grids. ! @@ -21,11 +21,11 @@ program test_get_index print*,"Starting test of get_index." -! The get_index routine contains three 'if' blocks. The test -! points are designed to check each block. +! The get_index routine contains a main 'if/elseif/else' statement. The test +! points are designed to check each branch. ! Point 1 - At equator. Western edge of grid cell at -! Greenwich. This tests the third 'if' block of the routine. +! Greenwich. This tests the 'else' branch. print*,"Test point 1." @@ -46,7 +46,7 @@ program test_get_index enddo ! Point 2 - At equator. Grid cell crosses Greenwich. -! This tests the second 'if' block of the routine. +! This tests the 'elseif' branch. print*,"Test point 2." @@ -73,7 +73,7 @@ program test_get_index enddo ! Point 3 - At equator. Grid cell centered at -! the dateline. This tests the third 'if' block. +! the dateline. This tests the 'else' branch. print*,"Test point 3." @@ -96,7 +96,7 @@ program test_get_index enddo ! Point 4 - At North Pole. Grid cell centered at -! Greenwich. This tests the first 'if' block. +! Greenwich. This tests the 'if' branch. print*,"Test point 4." From 5f831896112739b4b8fdf2b2cc33dd3307d254db Mon Sep 17 00:00:00 2001 From: George Gayno Date: Mon, 13 Jan 2025 21:46:09 +0000 Subject: [PATCH 59/67] Fix print statement in ./tests/orog/ftst_get_xnsum3.F90 Fixes #1000. --- tests/orog/ftst_get_xnsum3.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/orog/ftst_get_xnsum3.F90 b/tests/orog/ftst_get_xnsum3.F90 index b0f46aa02..75a0a8fe8 100644 --- a/tests/orog/ftst_get_xnsum3.F90 +++ b/tests/orog/ftst_get_xnsum3.F90 @@ -32,7 +32,7 @@ program check_get_xnsum3 ! grid box that are above the critical height. data expected_xnsum2 /16/ ! Expected total number of high-res pts in model grid box. - print*,"Begin test of routine get_xnsum2." + print*,"Begin test of routine get_xnsum3." ! The high-res grid is a global one-degree lat/lon grid. Point (1,1) ! is the south pole/greenwich. From 2197dff6f1edb419adbafdb5333105e2184520f2 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 14 Jan 2025 08:49:07 -0600 Subject: [PATCH 60/67] Update comments in ./tests/orog/ftst_qc_orog_by_ramp.F90 Fixes #1000. --- tests/orog/ftst_qc_orog_by_ramp.F90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/orog/ftst_qc_orog_by_ramp.F90 b/tests/orog/ftst_qc_orog_by_ramp.F90 index b15145120..06c4ad4cf 100644 --- a/tests/orog/ftst_qc_orog_by_ramp.F90 +++ b/tests/orog/ftst_qc_orog_by_ramp.F90 @@ -6,8 +6,9 @@ program qc_orog_ramp ! ! In OPS, the global data is 30-sec with dimensions ! 43200 x 21600. The RAMP data is 30-sec with dimension -! 43201 x 3601. For this test, reduced versions of -! both grids are used: global - 9x7; RAMP 10x5. +! 43201 x 3601 and only covers Antarctica. For this test, +! reduced versions of both grids are used: +! global - 9x7; RAMP 10x5. ! ! Author George Gayno NCEP/EMC @@ -79,7 +80,7 @@ program qc_orog_ramp zavg_expected(9,2) = 0 ! Note: this 'ramp' point has non-zero terrain of ! 0.14, which rounds down to zero. -! Note: The location of the RAMP data is set in the routine. +! Note: The path/name of the RAMP data is set in the routine. call qc_orog_by_ramp(imn, jmn, zavg, zslm) From 97314f466117dd7a2c20c4131003a1bd7694dccd Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 14 Jan 2025 09:13:33 -0600 Subject: [PATCH 61/67] Fix typo in four tests. Fixes #1000. --- tests/orog/ftst_inside_polygon.F90 | 2 +- tests/orog/ftst_read_global_mask.F90 | 2 +- tests/orog/ftst_read_global_orog.F90 | 2 +- tests/orog/ftst_rm_isolated_pts.F90 | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/orog/ftst_inside_polygon.F90 b/tests/orog/ftst_inside_polygon.F90 index b71e57c38..c9f12e99e 100644 --- a/tests/orog/ftst_inside_polygon.F90 +++ b/tests/orog/ftst_inside_polygon.F90 @@ -218,6 +218,6 @@ program inside_polygon if (inside) stop 18 ! Test point should be outside polygon. print*,"OK" - print*,"SUCCSSS" + print*,"SUCCESS" end program inside_polygon diff --git a/tests/orog/ftst_read_global_mask.F90 b/tests/orog/ftst_read_global_mask.F90 index 53af1fba3..9672a9391 100644 --- a/tests/orog/ftst_read_global_mask.F90 +++ b/tests/orog/ftst_read_global_mask.F90 @@ -33,6 +33,6 @@ program read_gbl_mask enddo print*,"OK" - print*,"SUCCSSS" + print*,"SUCCESS" end program read_gbl_mask diff --git a/tests/orog/ftst_read_global_orog.F90 b/tests/orog/ftst_read_global_orog.F90 index a5a854c0a..0f196243a 100644 --- a/tests/orog/ftst_read_global_orog.F90 +++ b/tests/orog/ftst_read_global_orog.F90 @@ -33,6 +33,6 @@ program read_gbl_orog enddo print*,"OK" - print*,"SUCCSSS" + print*,"SUCCESS" end program read_gbl_orog diff --git a/tests/orog/ftst_rm_isolated_pts.F90 b/tests/orog/ftst_rm_isolated_pts.F90 index 21cbdf42e..04375fd6b 100644 --- a/tests/orog/ftst_rm_isolated_pts.F90 +++ b/tests/orog/ftst_rm_isolated_pts.F90 @@ -212,6 +212,6 @@ program rm_isolated_pts deallocate (slm_expected, oro_expected, var_expected, var4_expected, oa_expected, ol_expected) print*,"OK" - print*,"SUCCSSS" + print*,"SUCCESS" end program rm_isolated_pts From e90b3f3efcec7fae73fe69b0bb1c35da22ebc461 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 14 Jan 2025 10:45:03 -0600 Subject: [PATCH 62/67] Remove some unused variables from routine make_mask. Fixes #1000. --- sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 index 474e07972..57cf5f2b7 100644 --- a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 +++ b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 @@ -425,7 +425,7 @@ SUBROUTINE MAKE_MASK(zslm,slm,land_frac, & real GLAT(JMN), GLON(IMN) real LONO(4),LATO(4),LONI,LATI real LONO_RAD(4), LATO_RAD(4) - integer JM1,i,j,nsum,nsum_all,ii,jj,numx,i2 + integer JM1,i,j,ii,jj,numx,i2 integer ilist(IMN) real DELXN,XNSUM,XLAND,XWATR,XL1,XS1,XW1 real XNSUM_ALL,XLAND_ALL,XWATR_ALL @@ -449,20 +449,18 @@ SUBROUTINE MAKE_MASK(zslm,slm,land_frac, & ! ! (*j*) for hard wired zero offset (lambda s =0) for terr05 !$omp parallel do & -!$omp private (j,i,xnsum,xland,xwatr,nsum,xl1,xs1,xw1,lono, & +!$omp private (j,i,xnsum,xland,xwatr,xl1,xs1,xw1,lono, & !$omp lato,lono_rad,lato_rad,jst,jen,ilist,numx,jj,i2,ii,loni,lati, & -!$omp xnsum_all,xland_all,xwatr_all,nsum_all) +!$omp xnsum_all,xland_all,xwatr_all) ! DO J=1,JM DO I=1,IM XNSUM = 0.0 XLAND = 0.0 XWATR = 0.0 - nsum = 0 XNSUM_ALL = 0.0 XLAND_ALL = 0.0 XWATR_ALL = 0.0 - nsum_all = 0 LONO(1) = lon_c(i,j) LONO(2) = lon_c(i+1,j) @@ -483,7 +481,6 @@ SUBROUTINE MAKE_MASK(zslm,slm,land_frac, & XLAND_ALL = XLAND_ALL + FLOAT(ZSLM(ii,jj)) XWATR_ALL = XWATR_ALL + FLOAT(1-ZSLM(ii,jj)) XNSUM_ALL = XNSUM_ALL + 1. - nsum_all = nsum_all+1 if(inside_a_polygon(LONI*D2R,LATI*D2R,4, & LONO_RAD,LATO_RAD))then @@ -491,7 +488,6 @@ SUBROUTINE MAKE_MASK(zslm,slm,land_frac, & XLAND = XLAND + FLOAT(ZSLM(ii,jj)) XWATR = XWATR + FLOAT(1-ZSLM(ii,jj)) XNSUM = XNSUM + 1. - nsum = nsum+1 endif enddo ; enddo From 55d45faa9931edfb7dc71bb3a848d10845ae892a Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 14 Jan 2025 12:41:16 -0600 Subject: [PATCH 63/67] Remove temporary logic from ./ush/fv3gfs_ocean_merge.sh. Fixes #1000. --- ush/fv3gfs_ocean_merge.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ush/fv3gfs_ocean_merge.sh b/ush/fv3gfs_ocean_merge.sh index ac8798f74..519a2e9da 100755 --- a/ush/fv3gfs_ocean_merge.sh +++ b/ush/fv3gfs_ocean_merge.sh @@ -20,7 +20,7 @@ cat << EOF > input.nml &mask_nml - ocean_mask_dir="${home_dir}/fix/orog.lowres/C${res}/ocean_mask/${ocn}/" + ocean_mask_dir="${home_dir}/fix/orog/C${res}/ocean_mask/${ocn}/" ocnres="mx${ocn}" lake_mask_dir="${TEMP_DIR}/C${res}/orog/" atmres="C${res}" From 28e0e0d69898631d9c33b436434228dc8927097b Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 14 Jan 2025 12:46:17 -0600 Subject: [PATCH 64/67] Remove temporary logic from ./util/gdas_init/run_v16.chgres.sh. Fixes #1000. --- util/gdas_init/run_v16.chgres.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/util/gdas_init/run_v16.chgres.sh b/util/gdas_init/run_v16.chgres.sh index e6a96efb3..d7d427c8a 100755 --- a/util/gdas_init/run_v16.chgres.sh +++ b/util/gdas_init/run_v16.chgres.sh @@ -1,4 +1,4 @@ -#!/bin/bash +!/bin/bash #--------------------------------------------------------------------------- # Run chgres using v16 netcdf history data as input. These history @@ -13,7 +13,7 @@ set -x MEMBER=$1 FIX_FV3=$UFS_DIR/fix -FIX_ORO=${FIX_FV3}/orog.lowres +FIX_ORO=${FIX_FV3}/orog FIX_AM=${FIX_FV3}/am WORKDIR=${WORKDIR:-$OUTDIR/work.${MEMBER}} From ac55299c082fba06533ddff7e14cbc3efb0caab9 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 14 Jan 2025 12:49:02 -0600 Subject: [PATCH 65/67] Fix typo in ./util/gdas_init/run_v16.chgres.sh. Fixes #1000. --- util/gdas_init/run_v16.chgres.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/util/gdas_init/run_v16.chgres.sh b/util/gdas_init/run_v16.chgres.sh index d7d427c8a..9b95843cb 100755 --- a/util/gdas_init/run_v16.chgres.sh +++ b/util/gdas_init/run_v16.chgres.sh @@ -1,4 +1,4 @@ -!/bin/bash +#!/bin/bash #--------------------------------------------------------------------------- # Run chgres using v16 netcdf history data as input. These history From a074ec42d046cc452bc8b808acc45664c1820d58 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 14 Jan 2025 12:51:29 -0600 Subject: [PATCH 66/67] For gdas_init utility, set the default ocean resolution to 9-degrees for the coarsest atmospheric resolutions. Fixes #1000. --- util/gdas_init/set_fixed_files.sh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/util/gdas_init/set_fixed_files.sh b/util/gdas_init/set_fixed_files.sh index 832614a59..451dd0bed 100755 --- a/util/gdas_init/set_fixed_files.sh +++ b/util/gdas_init/set_fixed_files.sh @@ -18,11 +18,11 @@ elif [ ${CTAR} == 'C768' ]; then elif [ ${CTAR} == 'C1152' ]; then OCNRES='025' elif [ ${CTAR} == 'C12' ]; then - OCNRES='100' + OCNRES='900' elif [ ${CTAR} == 'C18' ]; then - OCNRES='100' + OCNRES='900' elif [ ${CTAR} == 'C24' ]; then - OCNRES='100' + OCNRES='900' else OCNRES='025' fi From 64d26a45eae96b6d5808a9d4ffe012f1a4e10f79 Mon Sep 17 00:00:00 2001 From: George Gayno Date: Tue, 14 Jan 2025 12:58:42 -0600 Subject: [PATCH 67/67] Revert ./util/gdas_init/config to original settings. Fixes #1000. --- util/gdas_init/config | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/util/gdas_init/config b/util/gdas_init/config index c7bc10e15..01a98bbf7 100644 --- a/util/gdas_init/config +++ b/util/gdas_init/config @@ -51,14 +51,14 @@ # #----------------------------------------------------------- -EXTRACT_DIR=/scratch2/NCEPDEV/stmp1/George.Gayno/gdas.init/input -EXTRACT_DATA=no +EXTRACT_DIR=/lfs/h2/emc/stmp/$USER/gdas.init/input +EXTRACT_DATA=yes RUN_CHGRES=yes -yy=2024 -mm=11 -dd=16 +yy=2020 +mm=01 +dd=01 hh=06 use_v16retro=no @@ -67,12 +67,12 @@ LEVS=65 CDUMP=gfs -CRES_HIRES=C12 +CRES_HIRES=C192 CRES_ENKF=C96 UFS_DIR=$PWD/../.. -OUTDIR=/scratch2/NCEPDEV/stmp1/George.Gayno/gdas.init/output +OUTDIR=/lfs/h2/emc/stmp/$USER/gdas.init/output #--------------------------------------------------------- # Dont touch anything below here.