Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update fms_mixedmode branch #76

Merged
5 changes: 3 additions & 2 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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}
Expand Down
129 changes: 106 additions & 23 deletions model/fv_regional_bc.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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 &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do we need to change the bounds?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@junwang-noaa For information, please see NOAA-GFDL#219

regional_bc_bounds%je_south_uvs=jed+1
!
regional_bc_bounds%is_south_uvw=isd
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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.
!
Expand All @@ -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
Expand All @@ -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.
Expand All @@ -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.
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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?
Expand Down Expand Up @@ -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?
Expand Down Expand Up @@ -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

Expand Down
2 changes: 1 addition & 1 deletion tools/external_ic.F90
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down