Update orography code for low-res grids (#1012)

Fix logic that caused crash for very low-resolution grids. Update I/O routines to use F90 NetCDF calls. Add unit tests for 20 subroutines/functions. For complete details see #1000 and #1012. Fixes #1000.
ufs-community · Jan 17, 2025 · d4f0db2 · d4f0db2
1 parent 21c51d7
commit d4f0db2
Show file tree

Hide file tree

Showing 25 changed files with 1,886 additions and 254 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
diff --git a/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90 b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.F90
@@ -419,15 +419,13 @@ 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
       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
@@ -451,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) 
@@ -485,25 +481,13 @@ 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(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
 
                XLAND = XLAND + FLOAT(ZSLM(ii,jj))
                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
 
@@ -559,13 +543,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=20000000
       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
@@ -576,8 +559,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
@@ -590,7 +571,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)
+      ENDDO
+      ENDDO
+      print*,"- MAXSUM IS ", maxsum
+      allocate(hgt_1d(MAXSUM))
+      allocate(hgt_1d_all(MAXSUM))
 !
 !---- FIND THE AVERAGE OF THE MODES IN A GRID BOX
 !
@@ -649,7 +647,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
@@ -668,7 +666,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

diff --git 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
@@ -586,9 +586,10 @@ 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
-     IW=MOD(I+IM-2,IM)+1
-     IE=MOD(I,IM)+1
+   i_loop : DO I=2,IM-1
+! Check the points to the 'west' and 'east'.
+     IW=I-1
+     IE=I+1
      SLMA=SLM(IW,J)+SLM(IE,J)
      OROA=ORO(IW,J)+ORO(IE,J)
      VARA=VAR(IW,J)+VAR(IE,J)
@@ -599,67 +600,41 @@ 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)
-       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
-       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
-       IS=MOD(NINT(XS)-1,IM)+1
-       ISW=MOD(IS+IM-2,IM)+1
-       ISE=MOD(IS,IM)+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
-     ELSE
-       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
      DO K=1,4
        OAA(K)=OAA(K)/WGTA
        OLA(K)=OLA(K)/WGTA
      ENDDO
+! 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
@@ -671,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
@@ -699,7 +675,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/CMakeLists.txt b/tests/orog/CMakeLists.txt
@@ -8,7 +8,31 @@ 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 test directory.
+
+# 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)
+
+# 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)
+
+# 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)
+
+# 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)
+
+# 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)
 
 add_executable(ftst_ll2xyz ftst_ll2xyz.F90)
 add_test(NAME orog-ftst_ll2xyz COMMAND ftst_ll2xyz)
@@ -21,3 +45,63 @@ 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)
+
+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)
+
+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)
+
+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)
+
+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)
+
+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)
+
+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)
+
+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/C12_grid.tile1.nc b/tests/orog/data/C12_grid.tile1.nc
diff --git a/tests/orog/data/C48.mx500.tile1.nc b/tests/orog/data/C48.mx500.tile1.nc
diff --git a/tests/orog/data/landcover.umd.lowres.nc b/tests/orog/data/landcover.umd.lowres.nc
diff --git a/tests/orog/data/topography.antarctica.ramp.lowres.nc b/tests/orog/data/topography.antarctica.ramp.lowres.nc
diff --git a/tests/orog/data/topography.gmted2010.lowres.nc b/tests/orog/data/topography.gmted2010.lowres.nc