Update fms_mixedmode branch (#76)

* fixes io performance issues by making everyone a reader when io_layout=1,1 adds capability to use FMS feature to ignore data integrity checksums in restarts * rename enforce_rst_cksum to ignore_rst_cksum and change the default value for compatibility * updated UFS/GFS atmosphere.F90 driver as per @BinLiu-NOAA and @junwang-noaa * Regional decomposition test fix (when nrows_blend > 0) (NOAA-GFDL#194) * Add missing instance for hailwat * Regional bc blend changes to extend into interior halos and overlap on corners. Still not working for u and v. * atmosphere.F90 : add hailwat to check dyn_core.F90 : Fix from Jun Wang to correct sync of u,v fv_regional_bc.F90 : add check for nrows_blend > tile size; fix error when nrows_blend=1 * Explanatory comment added * Removed commented code * Clean old code * In fv_fill.F90, use kind_phys for kind_phys instead of hard-coding 8 byte reals. (NOAA-GFDL#193) * Expose remap_scalar and remap_dwinds to fv3-jedi (NOAA-GFDL#199) * changed interface to public * added public * removed source * mods for jedi build * Transfer changes from PR NOAA-GFDL#202 to NOAA-GFDL#199 Made small changes from PR NOAA-GFDL#202 manually. * returned ignore checksum * fixed ignore checksum * Fix several bugs in fv_regional_bc.F90 relating to uninitialized or incorrectly initialized memory. (NOAA-GFDL#219) * fixes and workarounds for uninitialized memory in fv_regional_bc * remove workarounds and fix remaining known bugs in ps_reg * a few more surface pressure bug fixes; now the test case runs in debug mode * workarounds and bug fixes from gnu compiler testing * remove -9999999 commented-out code * quiet the NaNs passed to Atmp%ps * simplify comments and explain snan * use i-1 & j-1 for two-point averages, when available * Replace many changes with PR NOAA-GFDL#220 Co-authored-by: Rusty.Benson <[email protected]> Co-authored-by: Ted Mansell <[email protected]> Co-authored-by: Rusty Benson <[email protected]> Co-authored-by: Samuel Trahan (NOAA contractor) <[email protected]> Co-authored-by: Mark Potts <[email protected]>
NOAA-EMC · Oct 21, 2022 · 0ab12b6 · 0ab12b6
1 parent dcff74d
commit 0ab12b6
Show file tree

Hide file tree

Showing 3 changed files with 110 additions and 26 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -54,7 +54,8 @@ endif()
 
 if(NOT FMS_FOUND)
   find_package(FMS REQUIRED COMPONENTS ${kind})
-  add_library(fms ALIAS FMS::fms_${kind})
+  string(TOLOWER ${kind} kind_lower)
+  add_library(fms ALIAS FMS::fms_${kind_lower})
 endif()
 
 list(APPEND moving_srcs
@@ -189,7 +190,7 @@ set_property(SOURCE model/fv_mapz.F90  APPEND_STRING PROPERTY COMPILE_FLAGS "${F
 set_property(SOURCE tools/fv_diagnostics.F90  APPEND_STRING PROPERTY COMPILE_FLAGS "-O0")
 
 set_target_properties(fv3 PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/include/fv3)
-target_include_directories(fv3 INTERFACE $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}/include/fv3>
+target_include_directories(fv3 PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}/include/fv3>
                                          $<INSTALL_INTERFACE:include/fv3>)
 
 target_include_directories(fv3 PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}

diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90
@@ -254,6 +254,67 @@ module fv_regional_mod
       logical :: data_source_fv3gfs
 contains
 
+
+!-----------------------------------------------------------------------
+!
+  logical function is_not_finite(val)
+!
+!-----------------------------------------------------------------------
+!*** This routine is equivalent to ".not. ieee_is_finite(val)"
+!*** which returns .true. for infinite and Not a Number (NaN), or
+!*** .false. otherwise. It's here as a workaround for this gfortran bug:
+!***
+!***    https://gcc.gnu.org/bugzilla/show_bug.cgi?id=82207
+!***
+!*** The compiler must use IEEE-standard floating point for this to work
+!-----------------------------------------------------------------------
+!
+    use, intrinsic :: iso_c_binding, only: c_int32_t, c_int64_t
+    implicit none
+!
+!-----------------------------------------------------------------------
+! Portability note: shiftr() is part of Fortran 2008, but it is widely
+! supported in older compilers.
+!-----------------------------------------------------------------------
+!
+    intrinsic shiftr, transfer, iand ! <--- declare intrinsic to help older compilers
+!
+!-----------------------------------------------------------------------
+! Use value-based argument passing instead of reference-based to avoid
+! signaling a NaN on conversion to addressable storage.
+!-----------------------------------------------------------------------
+!
+    real, value :: val ! <-- bit pattern to test for infinity or NaN
+!
+!-----------------------------------------------------------------------
+! Bit manipulation constants for testing 32-bit floating-point
+! non-finite values.
+!-----------------------------------------------------------------------
+!
+#ifdef OVERLOAD_R4
+    integer(c_int32_t), parameter :: check = 255 ! <-- all bits on, size of exponent (8 bits)
+    integer, parameter :: shift = 23 ! <-- number of mantissa bits except sign
+!
+!-----------------------------------------------------------------------
+! Bit manipulation constants for testing 64-bit floating-point
+! non-finite values.
+!-----------------------------------------------------------------------
+!
+#else
+    integer(c_int64_t), parameter :: check = 2047 ! <-- all bits on, size of exponent (11 bits)
+    integer, parameter :: shift = 52 ! <-- number of mantissa bits except sign
+#endif
+!
+!-----------------------------------------------------------------------
+! For IEEE standard floating point numbers, non-finite values follow
+! a mandatory bit pattern. They have the mantissa sign bit on, and all
+! exponent bits on, except the exponent sign which can be on or off.
+!-----------------------------------------------------------------------
+!
+    is_not_finite = iand(shiftr(transfer(val,check),shift),check)==check
+!
+  end function is_not_finite
+!
 !-----------------------------------------------------------------------
 !
       subroutine setup_regional_BC(Atm                                  &
@@ -765,8 +826,8 @@ subroutine setup_regional_BC(Atm                                  &
 !***  reference pressure profile.  Compute it now.
 !-----------------------------------------------------------------------
 !
-      allocate(pref(npz+1))
-      allocate(dum1d(npz+1))
+      allocate(pref(npz+1)) ; pref=real_snan
+      allocate(dum1d(npz+1)) ; dum1d=real_snan
 !
       ps1=101325.
       pref(npz+1)=ps1
@@ -951,7 +1012,7 @@ subroutine compute_regional_bc_indices(regional_bc_bounds)
         regional_bc_bounds%ie_north_uvs=ied
 !
         regional_bc_bounds%js_north_uvs=jsd
-        regional_bc_bounds%je_north_uvs=nrows_blend+1
+        regional_bc_bounds%je_north_uvs=nrows_blend
 !
         regional_bc_bounds%is_north_uvw=isd
         regional_bc_bounds%ie_north_uvw=ied+1
@@ -968,7 +1029,7 @@ subroutine compute_regional_bc_indices(regional_bc_bounds)
         regional_bc_bounds%is_south_uvs=isd
         regional_bc_bounds%ie_south_uvs=ied
 !
-        regional_bc_bounds%js_south_uvs=jed-nhalo_model-nrows_blend+1
+        regional_bc_bounds%js_south_uvs=jed-nhalo_model-nrows_blend+2
         regional_bc_bounds%je_south_uvs=jed+1
 !
         regional_bc_bounds%is_south_uvw=isd
@@ -1028,7 +1089,7 @@ subroutine compute_regional_bc_indices(regional_bc_bounds)
           regional_bc_bounds%je_west_uvs=jed-nhalo_model+1
         endif
 !
-        regional_bc_bounds%is_west_uvw=ied-nhalo_model-nrows_blend+1
+        regional_bc_bounds%is_west_uvw=ied-nhalo_model-nrows_blend+2
         regional_bc_bounds%ie_west_uvw=ied+1
 !
         regional_bc_bounds%js_west_uvw=jsd
@@ -1307,8 +1368,8 @@ subroutine start_regional_cold_start(Atm, ak, bk, levp            &
       enddo
       enddo
 !
-      allocate (ak_in(1:levp+1))                                           !<-- Save the input vertical structure for
-      allocate (bk_in(1:levp+1))                                           !    remapping BC updates during the forecast.
+      allocate (ak_in(1:levp+1)) ; ak_in=real_snan                         !<-- Save the input vertical structure for
+      allocate (bk_in(1:levp+1)) ; bk_in=real_snan                         !    remapping BC updates during the forecast.
       do k=1,levp+1
         ak_in(k)=ak(k)
         bk_in(k)=bk(k)
@@ -1402,9 +1463,9 @@ subroutine start_regional_restart(Atm                             &
                             ,isd, ied, jsd, jed                         &
                             ,Atm%npx, Atm%npy )
 !
-      allocate (wk2(levp+1,2))
-      allocate (ak_in(levp+1))                                               !<-- Save the input vertical structure for
-      allocate (bk_in(levp+1))                                               !    remapping BC updates during the forecast.
+      allocate (wk2(levp+1,2)) ; wk2=real_snan
+      allocate (ak_in(levp+1)) ; ak_in=real_snan                             !<-- Save the input vertical structure for
+      allocate (bk_in(levp+1)) ; bk_in=real_snan                             !    remapping BC updates during the forecast.
       if (Atm%flagstruct%hrrrv3_ic) then
         if (open_file(Grid_input, 'INPUT/hrrr_ctrl.nc', "read", pelist=pes)) then
           call read_data(Grid_input,'vcoord',wk2)
@@ -1908,7 +1969,8 @@ subroutine regional_bc_data(Atm,bc_hour                           &
 !***  the integration levels.
 !-----------------------------------------------------------------------
 !
-        allocate(ps_reg(is_input:ie_input,js_input:je_input)) ; ps_reg=-9999999 ! for now don't set to snan until remap dwinds is changed
+        allocate(ps_reg(is_input:ie_input,js_input:je_input))              !<-- Sfc pressure in domain's boundary region derived from BC files
+        ps_reg=real_snan !<-- detect access of uninitialized pressures
 !
 !-----------------------------------------------------------------------
 !***  We have the boundary variables from the BC file on the levels
@@ -3545,6 +3607,12 @@ subroutine remap_scalar_nggps_regional_bc(Atm                         &
         je=js_bc+nhalo_data+nrows_blend-1
       endif
 !
+! Ensure uninitialized memory isn't used
+      pn0 = real_snan
+      pn1 = real_snan
+      gz_fv = real_snan
+      gz = real_snan
+      pn = real_snan
       allocate(pe0(is:ie,km+1)) ; pe0=real_snan
       allocate(qn1(is:ie,npz)) ; qn1=real_snan
       allocate(dp2(is:ie,npz)) ; dp2=real_snan
@@ -3575,13 +3643,14 @@ subroutine remap_scalar_nggps_regional_bc(Atm                         &
            pn(k) = 2.*pn(km+1) - pn(l)
         enddo
 
+        pst = real_snan
         do k=km+k2-1, 2, -1
           if( phis_reg(i,j).le.gz(k) .and. phis_reg(i,j).ge.gz(k+1) ) then
             pst = pn(k) + (pn(k+1)-pn(k))*(gz(k)-phis_reg(i,j))/(gz(k)-gz(k+1))
-            go to 123
+            exit
           endif
         enddo
-  123   ps(i,j) = exp(pst)
+        ps(i,j) = exp(pst)
 
      enddo   ! i-loop
 
@@ -3594,10 +3663,10 @@ subroutine remap_scalar_nggps_regional_bc(Atm                         &
 !***  the Atm object.
 !---------------------------------------------------------------------------------
 !
-      is=lbound(Atm%ps,1)
-      ie=ubound(Atm%ps,1)
-      js=lbound(Atm%ps,2)
-      je=ubound(Atm%ps,2)
+      is=min(ubound(Atm%ps,1),max(lbound(Atm%ps,1),is))
+      ie=min(ubound(Atm%ps,1),max(lbound(Atm%ps,1),ie))
+      js=min(ubound(Atm%ps,2),max(lbound(Atm%ps,2),js))
+      je=min(ubound(Atm%ps,2),max(lbound(Atm%ps,2),je))
 !
       do j=js,je
       do i=is,ie
@@ -4864,6 +4933,9 @@ subroutine bc_time_interpolation(array                               &
           do j=j1_blend,j2_blend
             factor_dist=exp(-(blend_exp1+blend_exp2*(j-j_bc-1)*rdenom)) !<-- Exponential falloff of blending weights.
             do i=i1_blend,i2_blend
+              if(is_not_finite(array(i,j,k))) then
+                cycle ! Outside boundary
+              endif
               blend_value=bc_t0(i,j,k)                                &  !<-- Blend data interpolated
                          +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval  !    between t0 and t1.
 !
@@ -4883,6 +4955,9 @@ subroutine bc_time_interpolation(array                               &
           do j=j1_blend,j2_blend
             factor_dist=exp(-(blend_exp1+blend_exp2*(j_bc-j-1)*rdenom)) !<-- Exponential falloff of blending weights.
             do i=i1_blend,i2_blend
+              if(is_not_finite(array(i,j,k))) then
+                cycle ! Outside boundary
+              endif
               blend_value=bc_t0(i,j,k)                                &  !<-- Blend data interpolated
                          +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval  !    between t0 and t1.
               array(i,j,k)=(1.-factor_dist)*array(i,j,k)+factor_dist*blend_value
@@ -4900,6 +4975,9 @@ subroutine bc_time_interpolation(array                               &
         do k=1,ubnd_z
           do j=j1_blend,j2_blend
             do i=i1_blend,i2_blend
+              if(is_not_finite(array(i,j,k))) then
+                cycle ! Outside boundary
+              endif
 !
               blend_value=bc_t0(i,j,k)                                  &  !<-- Blend data interpolated
                          +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval    !    between t0 and t1.
@@ -4921,6 +4999,9 @@ subroutine bc_time_interpolation(array                               &
         do k=1,ubnd_z
           do j=j1_blend,j2_blend
             do i=i1_blend,i2_blend
+              if(is_not_finite(array(i,j,k))) then
+                cycle ! Outside boundary
+              endif
 !
               blend_value=bc_t0(i,j,k)                                  &  !<-- Blend data interpolated
                          +(bc_t1(i,j,k)-bc_t0(i,j,k))*fraction_interval    !    between t0 and t1.
@@ -5727,7 +5808,7 @@ subroutine dump_field_3d (domain, name, field, isd, ied, jsd, jed, nlev, stag)
     nyg = npy + 2*halo + jext
     nz = size(field,dim=3)
 
-    allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext, 1:nz) )
+    allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext, 1:nz) ) ; glob_field=real_snan
 
     isection_s = is
     isection_e = ie
@@ -5846,7 +5927,7 @@ subroutine dump_field_2d (domain, name, field, isd, ied, jsd, jed, stag)
     nxg = npx + 2*halo + iext
     nyg = npy + 2*halo + jext
 
-    allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext) )
+    allocate( glob_field(isg-halo:ieg+halo+iext, jsg-halo:jeg+halo+jext) ) ; glob_field=real_snan
 
     isection_s = is
     isection_e = ie
@@ -6154,6 +6235,7 @@ subroutine write_full_fields(Atm)
         nz=size(fields_core(nv)%ptr,3)
 !
         allocate( global_field(istart_g:iend_g, jstart_g:jend_g, 1:nz) )
+        global_field=real_snan
 !
 !-----------------------------------------------------------------------
 !***  What is the local extent of the variable on the task subdomain?
@@ -6258,6 +6340,7 @@ subroutine write_full_fields(Atm)
       nz=size(fields_tracers(1)%ptr,3)
 !
       allocate( global_field(istart_g:iend_g, jstart_g:jend_g, 1:nz) )
+      global_field=real_snan
 !
 !-----------------------------------------------------------------------
 !***  What is the local extent of the variable on the task subdomain?
@@ -6614,10 +6697,10 @@ subroutine exch_uv(domain, bd, npz, u, v)
     ! buf1 and buf4 must be of the same size (sim. for buf2 and buf3).
     ! Changes to the code below should be tested with debug flags
     ! enabled (out-of-bounds reads/writes).
-    allocate(buf1(1:24*npz))
-    allocate(buf2(1:36*npz))
-    allocate(buf3(1:36*npz))
-    allocate(buf4(1:24*npz))
+    allocate(buf1(1:24*npz)) ; buf1=real_snan
+    allocate(buf2(1:36*npz)) ; buf2=real_snan
+    allocate(buf3(1:36*npz)) ; buf3=real_snan
+    allocate(buf4(1:24*npz)) ; buf4=real_snan
 
 ! FIXME: MPI_COMM_WORLD
 

diff --git a/tools/external_ic.F90 b/tools/external_ic.F90
@@ -206,7 +206,7 @@ module external_ic_mod
 ! Include variable "version" to be written to log file.
 #include<file_version.h>
 
-   public get_external_ic, get_cubed_sphere_terrain
+   public get_external_ic, get_cubed_sphere_terrain, remap_scalar, remap_dwinds
 
 contains