From 59635994682bb7797e218d860c9a476e3e0ed755 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Sun, 18 Sep 2022 22:10:42 +0000 Subject: [PATCH 01/10] fixes and workarounds for uninitialized memory in fv_regional_bc --- model/fv_regional_bc.F90 | 549 +++++++++++++++++++++------------------ 1 file changed, 303 insertions(+), 246 deletions(-) diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index f9d6e8c2f..1e0f91191 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -263,6 +263,7 @@ subroutine setup_regional_BC(Atm & !----------------------------------------------------------------------- !*** Regional boundary data is obtained from the external BC file. !----------------------------------------------------------------------- + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite use netcdf !----------------------------------------------------------------------- implicit none @@ -765,8 +766,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 @@ -1307,8 +1308,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 +1403,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) @@ -1552,6 +1553,7 @@ subroutine regional_bc_data(Atm,bc_hour & !*** Regional boundary data is obtained from the external BC file. !----------------------------------------------------------------------- !----------------------------------------------------------------------- + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite implicit none !----------------------------------------------------------------------- ! @@ -1908,7 +1910,45 @@ 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 +! +!----------------------------------------------------------------------- +!*** Whoever did this intentionally added a bug to the code and hid it. +!*** I'm retaining this code, commented out, as a reminder to others. +!----------------------------------------------------------------------- +! + !ps_reg=-9999999 ! for now don't set to snan until remap dwinds is changed + ps_reg=real_snan +! +!----------------------------------------------------------------------- +!*** Initialize the domain boundary surface pressure to the input +!*** surface pressure on points remap dwinds can't handle. +!*** When data is missing or invalid, use nearby points. If all nearby +!*** points are invalid, use standard atmospheric pressure. +!----------------------------------------------------------------------- +! + do j=js_input,je_input + do i=is_input,ie_input + if(ieee_is_finite(ps_input(i,j,1))) then + ps_reg(i,j) = ps_input(i,j,1) + + ! When there is invalid data at a point, search nearby points. + else if(ieee_is_finite(ps_input(max(is_input,i-1),j,1))) then + ps_reg(i,j) = ps_input(max(is_input,i-1),j,1) + else if(ieee_is_finite(ps_input(min(ie_input,i+1),j,1))) then + ps_reg(i,j) = ps_input(min(ie_input,i+1),j,1) + else if(ieee_is_finite(ps_input(i,max(js_input,j-1),1))) then + ps_reg(i,j) = ps_input(i,max(js_input,j-1),1) + else if(ieee_is_finite(ps_input(i,min(je_input,j+1),1))) then + ps_reg(i,j) = ps_input(i,min(je_input,j+1),1) + + else + ! Handle severe bugs in reader or input conditions with + ! standard atmospheric pressure. (Should never get here.) + ps_reg(i,j) = 1e5 + endif + enddo + enddo ! !----------------------------------------------------------------------- !*** We have the boundary variables from the BC file on the levels @@ -1942,16 +1982,7 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif ! - if(nside==2)then - if(south_bc)then - call_remap=.true. - side='south' - bc_side_t1=>BC_t1%south - bc_side_t0=>BC_t0%south - endif - endif -! - if(nside==3)then + if(nside==2)then if(east_bc)then call_remap=.true. side='east ' @@ -1960,7 +1991,7 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif ! - if(nside==4)then + if(nside==3)then if(west_bc)then call_remap=.true. side='west ' @@ -1968,6 +1999,15 @@ subroutine regional_bc_data(Atm,bc_hour & bc_side_t0=>BC_t0%west endif endif +! + if(nside==4)then + if(south_bc)then + call_remap=.true. + side='south' + bc_side_t1=>BC_t1%south + bc_side_t0=>BC_t0%south + endif + endif ! if(call_remap)then call remap_scalar_nggps_regional_bc(Atm & @@ -2045,57 +2085,9 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif - if(nside==2)then - if(south_bc)then - if (is == 1) then - istart = 1 - else - istart = isd - endif - if (ie == npx-1) then - iend = npx-1 - else - iend = ied - endif - - do k=1,npz - do j=npy,jed - do i=istart,iend - delz_regBC%north_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) - delz_regBC%north_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) - enddo - enddo - enddo - - ! North, south include all corners - if (is == 1) then - do k=1,npz - do j=npy,jed - do i=isd,0 - delz_regBC%west_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) - delz_regBC%west_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) - enddo - enddo - enddo - endif - - - if (ie == npx-1) then - do k=1,npz - do j=npy,jed - do i=npx,ied - delz_regBC%east_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) - delz_regBC%east_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) - enddo - enddo - enddo - endif - endif - endif - ! - if(nside==3)then + if(nside==2)then if(east_bc)then if (js == 1) then jstart = 1 @@ -2120,7 +2112,7 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif - if(nside==4)then + if(nside==3)then if(west_bc)then if (js == 1) then jstart = 1 @@ -2145,6 +2137,54 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif + if(nside==4)then + if(south_bc)then + if (is == 1) then + istart = 1 + else + istart = isd + endif + if (ie == npx-1) then + iend = npx-1 + else + iend = ied + endif + + do k=1,npz + do j=npy,jed + do i=istart,iend + delz_regBC%north_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) + delz_regBC%north_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) + enddo + enddo + enddo + + ! North, south include all corners + if (is == 1) then + do k=1,npz + do j=npy,jed + do i=isd,0 + delz_regBC%west_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) + delz_regBC%west_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) + enddo + enddo + enddo + endif + + + if (ie == npx-1) then + do k=1,npz + do j=npy,jed + do i=npx,ied + delz_regBC%east_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) + delz_regBC%east_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) + enddo + enddo + enddo + endif + endif + endif + endif ! !----------------------------------------------------------------------- @@ -2225,23 +2265,6 @@ subroutine regional_bc_data(Atm,bc_hour & endif ! if(nside==2)then - if(south_bc)then - call_remap=.true. - bc_side_t1=>BC_t1%south -! - is_u=Atm%regional_bc_bounds%is_south_uvs - ie_u=Atm%regional_bc_bounds%ie_south_uvs - js_u=Atm%regional_bc_bounds%js_south_uvs - je_u=Atm%regional_bc_bounds%je_south_uvs -! - is_v=Atm%regional_bc_bounds%is_south_uvw - ie_v=Atm%regional_bc_bounds%ie_south_uvw - js_v=Atm%regional_bc_bounds%js_south_uvw - je_v=Atm%regional_bc_bounds%je_south_uvw - endif - endif -! - if(nside==3)then if(east_bc)then call_remap=.true. bc_side_t1=>BC_t1%east @@ -2258,7 +2281,7 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif ! - if(nside==4)then + if(nside==3)then if(west_bc)then call_remap=.true. bc_side_t1=>BC_t1%west @@ -2274,6 +2297,23 @@ subroutine regional_bc_data(Atm,bc_hour & je_v=Atm%regional_bc_bounds%je_west_uvw endif endif +! + if(nside==4)then + if(south_bc)then + call_remap=.true. + bc_side_t1=>BC_t1%south +! + is_u=Atm%regional_bc_bounds%is_south_uvs + ie_u=Atm%regional_bc_bounds%ie_south_uvs + js_u=Atm%regional_bc_bounds%js_south_uvs + je_u=Atm%regional_bc_bounds%je_south_uvs +! + is_v=Atm%regional_bc_bounds%is_south_uvw + ie_v=Atm%regional_bc_bounds%ie_south_uvw + js_v=Atm%regional_bc_bounds%js_south_uvw + je_v=Atm%regional_bc_bounds%je_south_uvw + endif + endif ! if(call_remap)then ! @@ -2767,17 +2807,6 @@ subroutine fill_divgd_BC endif ! if(nside==2)then - if(south_bc)then - call_set=.true. - bc_side_t1=>BC_t1%south - is0=lbound(BC_t1%south%divgd_BC,1) - ie0=ubound(BC_t1%south%divgd_BC,1) - js0=lbound(BC_t1%south%divgd_BC,2) - je0=ubound(BC_t1%south%divgd_BC,2) - endif - endif -! - if(nside==3)then if(east_bc)then call_set=.true. bc_side_t1=>BC_t1%east @@ -2788,7 +2817,7 @@ subroutine fill_divgd_BC endif endif ! - if(nside==4)then + if(nside==3)then if(west_bc)then call_set=.true. bc_side_t1=>BC_t1%west @@ -2798,6 +2827,17 @@ subroutine fill_divgd_BC je0=ubound(BC_t1%west%divgd_BC,2) endif endif +! + if(nside==4)then + if(south_bc)then + call_set=.true. + bc_side_t1=>BC_t1%south + is0=lbound(BC_t1%south%divgd_BC,1) + ie0=ubound(BC_t1%south%divgd_BC,1) + js0=lbound(BC_t1%south%divgd_BC,2) + je0=ubound(BC_t1%south%divgd_BC,2) + endif + endif ! if(call_set)then do k=1,klev_out @@ -2860,17 +2900,6 @@ subroutine fill_q_con_BC endif ! if(nside==2)then - if(south_bc)then - call_set=.true. - bc_side_t1=>BC_t1%south - is0=lbound(BC_t1%south%q_con_BC,1) - ie0=ubound(BC_t1%south%q_con_BC,1) - js0=lbound(BC_t1%south%q_con_BC,2) - je0=ubound(BC_t1%south%q_con_BC,2) - endif - endif -! - if(nside==3)then if(east_bc)then call_set=.true. bc_side_t1=>BC_t1%east @@ -2881,7 +2910,7 @@ subroutine fill_q_con_BC endif endif ! - if(nside==4)then + if(nside==3)then if(west_bc)then call_set=.true. bc_side_t1=>BC_t1%west @@ -2891,6 +2920,17 @@ subroutine fill_q_con_BC je0=ubound(BC_t1%west%q_con_BC,2) endif endif +! + if(nside==4)then + if(south_bc)then + call_set=.true. + bc_side_t1=>BC_t1%south + is0=lbound(BC_t1%south%q_con_BC,1) + ie0=ubound(BC_t1%south%q_con_BC,1) + js0=lbound(BC_t1%south%q_con_BC,2) + je0=ubound(BC_t1%south%q_con_BC,2) + endif + endif ! if(call_set)then do k=1,klev_out @@ -2948,23 +2988,23 @@ subroutine fill_cappa_BC endif ! if(nside==2)then - if(south_bc)then + if(east_bc)then call_compute=.true. - bc_side_t1=>BC_t1%south + bc_side_t1=>BC_t1%east endif endif ! if(nside==3)then - if(east_bc)then + if(west_bc)then call_compute=.true. - bc_side_t1=>BC_t1%east + bc_side_t1=>BC_t1%west endif endif ! if(nside==4)then - if(west_bc)then + if(south_bc)then call_compute=.true. - bc_side_t1=>BC_t1%west + bc_side_t1=>BC_t1%south endif endif ! @@ -3179,34 +3219,11 @@ subroutine read_regional_bc_file(is_input,ie_input & endif endif ! -!----------- -!*** South -!----------- -! - if(nside==2)then - if(south_bc)then - call_get_var=.true. -! - var_name=trim(var_name_root)//"_top" -! - i_start_array=is_input - i_end_array =ie_input - j_start_array=je_input-nhalo_data-nrows_blend+1 - j_end_array =je_input -! - i_start_data=i_start_array+nhalo_data - i_count=i_end_array-i_start_array+1 - j_start_data=1 - j_count=j_end_array-j_start_array+1 -! - endif - endif -! !---------- !*** East !---------- ! - if(nside==3)then + if(nside==2)then if(east_bc)then call_get_var=.true. ! @@ -3250,7 +3267,7 @@ subroutine read_regional_bc_file(is_input,ie_input & !*** West !---------- ! - if(nside==4)then + if(nside==3)then if(west_bc)then call_get_var=.true. ! @@ -3286,6 +3303,29 @@ subroutine read_regional_bc_file(is_input,ie_input & endif endif ! +!----------- +!*** South +!----------- +! + if(nside==4)then + if(south_bc)then + call_get_var=.true. +! + var_name=trim(var_name_root)//"_top" +! + i_start_array=is_input + i_end_array =ie_input + j_start_array=je_input-nhalo_data-nrows_blend+1 + j_end_array =je_input +! + i_start_data=i_start_array+nhalo_data + i_count=i_end_array-i_start_array+1 + j_start_data=1 + j_count=j_end_array-j_start_array+1 +! + endif + endif +! !----------------------------------------------------------------------- !*** Fill this task's subset of boundary data for this 3-D !*** or 4-D variable. This includes the data in the domain's @@ -3545,6 +3585,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm & je=js_bc+nhalo_data+nrows_blend-1 endif ! + pn0 = 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 +3616,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 @@ -4431,61 +4473,11 @@ subroutine regional_boundary_update(array & endif endif ! -!----------- -!*** South -!----------- -! - if(nside==2)then - if(south_bc)then - call_interp=.true. - bc_side_t0=>bc_south_t0 - bc_side_t1=>bc_south_t1 -! - i1=isd - i2=ied - if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then - i2=ied+1 - endif -! - j1=je+1 - j2=jed - if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then - j1=je+2 - j2=jed+1 - endif -! - i1_blend=is - i2_blend=ie - IF ( east_bc ) THEN - i1_blend=is - ELSE - i1_blend=isd !is-nhalo_model - ENDIF - IF ( west_bc ) THEN - i2_blend=ie - ELSE - i2_blend=ied ! ie+nhalo_model - ENDIF - if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then -! i2_blend=ie+1 - i2_blend=i2_blend+1 - endif - j2_blend=je - if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then - j2_blend=je+1 - endif - j1_blend=j2_blend-nrows_blend_user+1 - i_bc=-9e9 - j_bc=j1 -! - endif - endif -! !---------- !*** East !---------- ! - if(nside==3)then + if(nside==2)then if(east_bc)then call_interp=.true. bc_side_t0=>bc_east_t0 @@ -4541,7 +4533,7 @@ subroutine regional_boundary_update(array & !*** West !---------- ! - if(nside==4)then + if(nside==3)then if(west_bc)then call_interp=.true. bc_side_t0=>bc_west_t0 @@ -4596,6 +4588,56 @@ subroutine regional_boundary_update(array & endif endif ! +!----------- +!*** South +!----------- +! + if(nside==4)then + if(south_bc)then + call_interp=.true. + bc_side_t0=>bc_south_t0 + bc_side_t1=>bc_south_t1 +! + i1=isd + i2=ied + if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then + i2=ied+1 + endif +! + j1=je+1 + j2=jed + if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then + j1=je+2 + j2=jed+1 + endif +! + i1_blend=is + i2_blend=ie + IF ( east_bc ) THEN + i1_blend=is + ELSE + i1_blend=isd !is-nhalo_model + ENDIF + IF ( west_bc ) THEN + i2_blend=ie + ELSE + i2_blend=ied ! ie+nhalo_model + ENDIF + if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then +! i2_blend=ie+1 + i2_blend=i2_blend+1 + endif + j2_blend=je + if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then + j2_blend=je+1 + endif + j1_blend=j2_blend-nrows_blend_user+1 + i_bc=-9e9 + j_bc=j1 +! + endif + endif +! !--------------------------------------------------------------------- !*** Get the pointers pointing at the boundary arrays holding the !*** two time levels of the given prognostic array's boundary region @@ -4758,6 +4800,7 @@ subroutine bc_time_interpolation(array & !*** forecast time that is within the interval bracketed by the !*** two current boundary region states. !--------------------------------------------------------------------- + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite implicit none !--------------------------------------------------------------------- ! @@ -4864,6 +4907,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(.not.ieee_is_finite(array(i,j,k))) then + cycle ! Bad data from caller + 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. ! @@ -4873,34 +4919,19 @@ subroutine bc_time_interpolation(array & enddo endif ! -!----------- -!*** South -!----------- -! - if(nside==2.and.south_bc)then - rdenom=1./real(Max(1,j_bc-j1_blend-1)) - do k=1,ubnd_z - 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 - 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 - enddo - enddo - enddo - endif -! !---------- !*** East !---------- ! - if(nside==3.and.east_bc)then + if(nside==2.and.east_bc)then rdenom=1./real(Max(1,i2_blend-i_bc-1)) do k=1,ubnd_z do j=j1_blend,j2_blend do i=i1_blend,i2_blend ! + if(.not.ieee_is_finite(array(i,j,k))) then + cycle ! Bad data from caller + 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. ! @@ -4916,12 +4947,15 @@ subroutine bc_time_interpolation(array & !*** West !---------- ! - if(nside==4.and.west_bc)then + if(nside==3.and.west_bc)then rdenom=1./real(Max(1, i_bc-i1_blend-1)) do k=1,ubnd_z do j=j1_blend,j2_blend do i=i1_blend,i2_blend ! + if(.not.ieee_is_finite(array(i,j,k))) then + cycle ! Bad data from caller + 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. ! @@ -4933,6 +4967,27 @@ subroutine bc_time_interpolation(array & enddo endif ! +!----------- +!*** South +!----------- +! + if(nside==4.and.south_bc)then + rdenom=1./real(Max(1,j_bc-j1_blend-1)) + do k=1,ubnd_z + 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(.not.ieee_is_finite(array(i,j,k))) then + cycle ! Bad data from caller + 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 + enddo + enddo + enddo + endif +! !--------------------------------------------------------------------- ! end subroutine bc_time_interpolation @@ -5011,29 +5066,6 @@ subroutine regional_bc_t1_to_t0(BC_t1,BC_t0 & endif ! if(nside==2)then - if(south_bc)then - move=.true. - bc_side_t0=>BC_t0%south - bc_side_t1=>BC_t1%south -! - is_c=bnds%is_south !<-- - ie_c=bnds%ie_south ! South BC index limits - js_c=bnds%js_south ! for centers of grid cells. - je_c=bnds%je_south !<-- -! - is_s=bnds%is_south_uvs !<-- - ie_s=bnds%ie_south_uvs ! South BC index limits - js_s=bnds%js_south_uvs ! for winds on N/S sides of grid cells. - je_s=bnds%je_south_uvs !<-- -! - is_w=bnds%is_south_uvw !<-- - ie_w=bnds%ie_south_uvw ! South BC index limits - js_w=bnds%js_south_uvw ! for winds on E/W sides of grid cells. - je_w=bnds%je_south_uvw !<-- - endif - endif -! - if(nside==3)then if(east_bc)then move=.true. bc_side_t0=>BC_t0%east @@ -5056,7 +5088,7 @@ subroutine regional_bc_t1_to_t0(BC_t1,BC_t0 & endif endif ! - if(nside==4)then + if(nside==3)then if(west_bc)then move=.true. bc_side_t0=>BC_t0%west @@ -5078,6 +5110,29 @@ subroutine regional_bc_t1_to_t0(BC_t1,BC_t0 & je_w=bnds%je_west_uvw !<-- endif endif +! + if(nside==4)then + if(south_bc)then + move=.true. + bc_side_t0=>BC_t0%south + bc_side_t1=>BC_t1%south +! + is_c=bnds%is_south !<-- + ie_c=bnds%ie_south ! South BC index limits + js_c=bnds%js_south ! for centers of grid cells. + je_c=bnds%je_south !<-- +! + is_s=bnds%is_south_uvs !<-- + ie_s=bnds%ie_south_uvs ! South BC index limits + js_s=bnds%js_south_uvs ! for winds on N/S sides of grid cells. + je_s=bnds%je_south_uvs !<-- +! + is_w=bnds%is_south_uvw !<-- + ie_w=bnds%ie_south_uvw ! South BC index limits + js_w=bnds%js_south_uvw ! for winds on E/W sides of grid cells. + je_w=bnds%je_south_uvw !<-- + endif + endif ! if(move)then do k=1,nlev @@ -5727,7 +5782,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 +5901,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 +6209,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 +6314,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 +6671,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 From 3766d58ea26491b09db56694224b29270abdc6b6 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Tue, 20 Sep 2022 12:23:00 +0000 Subject: [PATCH 02/10] remove workarounds and fix remaining known bugs in ps_reg --- model/fv_regional_bc.F90 | 50 ++++++++++++++-------------------------- 1 file changed, 17 insertions(+), 33 deletions(-) diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index 1e0f91191..da19be025 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -1913,42 +1913,21 @@ subroutine regional_bc_data(Atm,bc_hour & allocate(ps_reg(is_input:ie_input,js_input:je_input)) !<-- Sfc pressure in domain's boundary region derived from BC files ! !----------------------------------------------------------------------- -!*** Whoever did this intentionally added a bug to the code and hid it. -!*** I'm retaining this code, commented out, as a reminder to others. +!*** Whoever did this intentionally added a bug to the code and hid it, +!*** to avoid fixing other bugs elsewhere. +!*** The result is that RRFS debug mode was broken for a long time +!*** and some boundary data was gibberish. +!*** I'm retaining this code, commented out, as a reminder to others +!*** of what NOT to do. !----------------------------------------------------------------------- ! !ps_reg=-9999999 ! for now don't set to snan until remap dwinds is changed - ps_reg=real_snan ! !----------------------------------------------------------------------- -!*** Initialize the domain boundary surface pressure to the input -!*** surface pressure on points remap dwinds can't handle. -!*** When data is missing or invalid, use nearby points. If all nearby -!*** points are invalid, use standard atmospheric pressure. +!*** Initialize ps_reg to real_snan to detect errors later in code. !----------------------------------------------------------------------- ! - do j=js_input,je_input - do i=is_input,ie_input - if(ieee_is_finite(ps_input(i,j,1))) then - ps_reg(i,j) = ps_input(i,j,1) - - ! When there is invalid data at a point, search nearby points. - else if(ieee_is_finite(ps_input(max(is_input,i-1),j,1))) then - ps_reg(i,j) = ps_input(max(is_input,i-1),j,1) - else if(ieee_is_finite(ps_input(min(ie_input,i+1),j,1))) then - ps_reg(i,j) = ps_input(min(ie_input,i+1),j,1) - else if(ieee_is_finite(ps_input(i,max(js_input,j-1),1))) then - ps_reg(i,j) = ps_input(i,max(js_input,j-1),1) - else if(ieee_is_finite(ps_input(i,min(je_input,j+1),1))) then - ps_reg(i,j) = ps_input(i,min(je_input,j+1),1) - - else - ! Handle severe bugs in reader or input conditions with - ! standard atmospheric pressure. (Should never get here.) - ps_reg(i,j) = 1e5 - endif - enddo - enddo + ps_reg=real_snan ! !----------------------------------------------------------------------- !*** We have the boundary variables from the BC file on the levels @@ -3586,6 +3565,11 @@ subroutine remap_scalar_nggps_regional_bc(Atm & endif ! 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 @@ -3636,10 +3620,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(ie,max(is,lbound(Atm%ps,1))) + ie=min(ie,max(is,ubound(Atm%ps,1))) + js=min(je,max(js,lbound(Atm%ps,2))) + je=min(je,max(js,ubound(Atm%ps,2))) ! do j=js,je do i=is,ie From 12c9bd1eb6386137664c9db49a0ef4eb343fa5b8 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Tue, 20 Sep 2022 15:30:25 +0000 Subject: [PATCH 03/10] a few more surface pressure bug fixes; now the test case runs in debug mode --- model/fv_regional_bc.F90 | 124 +++++++++++++++++++++++++++++++++------ 1 file changed, 106 insertions(+), 18 deletions(-) diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index da19be025..e734ccd1f 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -3482,7 +3482,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm & ,npz & !<-- # of levels in final 3-D integration variables ,ncnst !<-- # of tracer variables real, intent(in):: ak0(km+1), bk0(km+1) - real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc):: psc + real, intent(inout), dimension(is_bc:ie_bc,js_bc:je_bc):: psc real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km):: t_in real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km):: omga real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km,ncnst):: qa @@ -3510,7 +3510,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm & real(kind=R_GRID), dimension(is_bc:ie_bc,km+1):: pn0 real(kind=R_GRID):: pst !!! High-precision - integer i,ie,is,j,je,js,k,l,m, k2,iq + integer i,ie,is,j,je,js,k,l,m, k2,iq, iss,ies,jss,jes, isv,iev,jsv,jev integer sphum, o3mr, liq_wat, ice_wat, rainwat, snowwat, graupel, hailwat, cld_amt ! !--------------------------------------------------------------------------------- @@ -3544,26 +3544,77 @@ subroutine remap_scalar_nggps_regional_bc(Atm & !*** This is needed to obtain pressures that will surround the wind points. !--------------------------------------------------------------------------------- ! - is=is_bc + isv=is_bc + iss=isv if(side=='west')then - is=ie_bc-nhalo_data-nrows_blend+1 + isv=ie_bc-nhalo_data-nrows_blend + iss=isv+1 endif ! - ie=ie_bc + iev=ie_bc + ies=iev if(side=='east')then - ie=is_bc+nhalo_data+nrows_blend-1 + iev=is_bc+nhalo_data+nrows_blend + ies=iev-1 endif ! - js=js_bc + jsv=js_bc + jss=jsv if(side=='south')then - js=je_bc-nhalo_data-nrows_blend+1 + jsv=je_bc-nhalo_data-nrows_blend + jss=jsv+1 endif ! - je=je_bc + jev=je_bc + jes=jev if(side=='north')then - je=js_bc+nhalo_data+nrows_blend-1 + jev=js_bc+nhalo_data+nrows_blend + jes=jev-1 endif ! +!--------------------------------------------------------------------------------- +!*** Workaround for incomplete 'ps' in boundary condition files. Copy from +!*** scalar halo regions (s) to velocity halo regions (v) beyond them. +!--------------------------------------------------------------------------------- +! + if(side=='west')then + do j=jss,jes + psc(isv,j) = psc(iss,j) + enddo + psc(iev,jsv) = psc(ies,jss) + psc(iev,jev) = psc(ies,jes) + endif +! + if(side=='east')then + do j=jss,jes + psc(iev,j) = psc(ies,j) + enddo + psc(iev,jev) = psc(ies,jes) + psc(iev,jev) = psc(ies,jes) + endif +! + if(side=='south')then + do i=iss,ies + psc(i,jsv) = psc(i,jss) + enddo + psc(isv,jsv) = psc(iss,jss) + psc(iev,jsv) = psc(ies,jss) + endif +! + if(side=='north')then + do i=iss,ies + psc(i,jev) = psc(i,jes) + enddo + psc(isv,jev) = psc(iss,jes) + psc(iev,jev) = psc(ies,jes) + endif +! + is = iss + js = jss + ie = ies + je = jes +! +! Ensure uninitialized memory isn't used pn0 = real_snan pn1 = real_snan gz_fv = real_snan @@ -3614,16 +3665,53 @@ subroutine remap_scalar_nggps_regional_bc(Atm & !--------------------------------------------------------------------------------- enddo jloop1 !--------------------------------------------------------------------------------- +! +!--------------------------------------------------------------------------------- +!*** Workaround for incomplete 'ps' in boundary condition files. Copy from +!*** scalar halo regions (s) to velocity halo regions (v) beyond them. +!--------------------------------------------------------------------------------- +! + if(side=='west')then + do j=jss,jes + ps(isv,j) = ps(iss,j) + enddo + ps(iev,jsv) = ps(ies,jss) + ps(iev,jev) = ps(ies,jes) + endif +! + if(side=='east')then + do j=jss,jes + ps(iev,j) = ps(ies,j) + enddo + ps(iev,jev) = ps(ies,jes) + ps(iev,jev) = ps(ies,jes) + endif +! + if(side=='south')then + do i=iss,ies + ps(i,jsv) = ps(i,jss) + enddo + ps(isv,jsv) = ps(iss,jss) + ps(iev,jsv) = ps(ies,jss) + endif +! + if(side=='north')then + do i=iss,ies + ps(i,jev) = ps(i,jes) + enddo + ps(isv,jev) = ps(iss,jes) + ps(iev,jev) = ps(ies,jes) + endif !--------------------------------------------------------------------------------- !*** Transfer values from the expanded boundary array for sfc pressure into !*** the Atm object. !--------------------------------------------------------------------------------- ! - is=min(ie,max(is,lbound(Atm%ps,1))) - ie=min(ie,max(is,ubound(Atm%ps,1))) - js=min(je,max(js,lbound(Atm%ps,2))) - je=min(je,max(js,ubound(Atm%ps,2))) + is=min(ubound(Atm%ps,1),max(lbound(Atm%ps,1),isv)) + ie=min(ubound(Atm%ps,1),max(lbound(Atm%ps,1),iev)) + js=min(ubound(Atm%ps,2),max(lbound(Atm%ps,2),jsv)) + je=min(ubound(Atm%ps,2),max(lbound(Atm%ps,2),jev)) ! do j=js,je do i=is,ie @@ -3960,12 +4048,12 @@ subroutine remap_dwinds_regional_bc(Atm & !------ do k=1,km+1 do i=is_u,ie_u - pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i,j-1)+psc(i,j)) + pe0(i,k) = ak0(k) + bk0(k)*psc(i,j) enddo enddo do k=1,npz+1 do i=is_u,ie_u - pe1(i,k) = Atm%ak(k) + Atm%bk(k)*0.5*(psc(i,j-1)+psc(i,j)) + pe1(i,k) = Atm%ak(k) + Atm%bk(k)*psc(i,j) enddo enddo call mappm(km, pe0(is_u:ie_u,1:km+1), ud(is_u:ie_u,j,1:km), npz, pe1(is_u:ie_u,1:npz+1), & @@ -4001,12 +4089,12 @@ subroutine remap_dwinds_regional_bc(Atm & do k=1,km+1 do i=is_v,ie_v - pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i-1,j)+psc(i,j)) + pe0(i,k) = ak0(k) + bk0(k)*psc(i,j) enddo enddo do k=1,npz+1 do i=is_v,ie_v - pe1(i,k) = Atm%ak(k) + Atm%bk(k)*0.5*(psc(i-1,j)+psc(i,j)) + pe1(i,k) = Atm%ak(k) + Atm%bk(k)*psc(i,j) enddo enddo call mappm(km, pe0(is_v:ie_v,1:km+1), vd(is_v:ie_v,j,1:km), npz, pe1(is_v:ie_v,1:npz+1), & From 73b7e5945eb78a44c8fe93e420cf52941e14fc50 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Tue, 20 Sep 2022 17:00:46 +0000 Subject: [PATCH 04/10] bug fixes to my bug fixes --- model/fv_regional_bc.F90 | 41 ++++++++++++---------------------------- 1 file changed, 12 insertions(+), 29 deletions(-) diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index e734ccd1f..7cb071350 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -1914,11 +1914,10 @@ subroutine regional_bc_data(Atm,bc_hour & ! !----------------------------------------------------------------------- !*** Whoever did this intentionally added a bug to the code and hid it, -!*** to avoid fixing other bugs elsewhere. -!*** The result is that RRFS debug mode was broken for a long time -!*** and some boundary data was gibberish. -!*** I'm retaining this code, commented out, as a reminder to others -!*** of what NOT to do. +!*** to conceal fixing other bugs elsewhere. This produced gibberish +!*** initial surface pressure, and made it impossible to test RRFS in +!*** debug mode for a long time. I'm retaining this code, commented out, +!*** as a reminder to others of what NOT to do. !----------------------------------------------------------------------- ! !ps_reg=-9999999 ! for now don't set to snan until remap dwinds is changed @@ -3578,35 +3577,27 @@ subroutine remap_scalar_nggps_regional_bc(Atm & !--------------------------------------------------------------------------------- ! if(side=='west')then - do j=jss,jes + do j=jsv,jev psc(isv,j) = psc(iss,j) enddo - psc(iev,jsv) = psc(ies,jss) - psc(iev,jev) = psc(ies,jes) endif ! if(side=='east')then - do j=jss,jes + do j=jsv,jev psc(iev,j) = psc(ies,j) enddo - psc(iev,jev) = psc(ies,jes) - psc(iev,jev) = psc(ies,jes) endif ! if(side=='south')then - do i=iss,ies + do i=isv,iev psc(i,jsv) = psc(i,jss) enddo - psc(isv,jsv) = psc(iss,jss) - psc(iev,jsv) = psc(ies,jss) endif ! if(side=='north')then - do i=iss,ies + do i=isv,iev psc(i,jev) = psc(i,jes) enddo - psc(isv,jev) = psc(iss,jes) - psc(iev,jev) = psc(ies,jes) endif ! is = iss @@ -3672,35 +3663,27 @@ subroutine remap_scalar_nggps_regional_bc(Atm & !--------------------------------------------------------------------------------- ! if(side=='west')then - do j=jss,jes + do j=jsv,jev ps(isv,j) = ps(iss,j) enddo - ps(iev,jsv) = ps(ies,jss) - ps(iev,jev) = ps(ies,jes) endif ! if(side=='east')then - do j=jss,jes + do j=jsv,jev ps(iev,j) = ps(ies,j) enddo - ps(iev,jev) = ps(ies,jes) - ps(iev,jev) = ps(ies,jes) endif ! if(side=='south')then - do i=iss,ies + do i=isv,iev ps(i,jsv) = ps(i,jss) enddo - ps(isv,jsv) = ps(iss,jss) - ps(iev,jsv) = ps(ies,jss) endif ! if(side=='north')then - do i=iss,ies + do i=isv,iev ps(i,jev) = ps(i,jes) enddo - ps(isv,jev) = ps(iss,jes) - ps(iev,jev) = ps(ies,jes) endif !--------------------------------------------------------------------------------- From 8c67895997b10a1a063f038f26c55d5304b81ec6 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Wed, 21 Sep 2022 14:40:48 +0000 Subject: [PATCH 05/10] workarounds and bug fixes from gnu compiler testing --- model/fv_regional_bc.F90 | 160 ++++++++++++++++++++++++++++++++++----- 1 file changed, 141 insertions(+), 19 deletions(-) diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index 7cb071350..b39651364 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -254,6 +254,75 @@ 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 does not work with signaling NaN in gfortran < version 12. +!*** It has the added advantage of being easily inlined and vectorized. +!*** It will detect all possible infinite and not-a-number (NaN) values +!*** in standard-conforming 32 bit and 64 bit floating-point. This is +!*** used to detect access beyond boundary regions in the blending code, +!*** which has no control over how the caller stores data in memory. +!*** +!*** The gfortran bug is documented here: +!*** +!*** https://gcc.gnu.org/bugzilla/show_bug.cgi?id=82207 +!*** +!*** A sufficiently intelligent compiler will optimize this function +!*** away, into three assembly instructions that can be vectorized. +!*** +!*** Technically, this routine needs Fortran 2008 due to shiftr() but +!*** that intrinsic is widely implemented in compilers that lack +!*** support for Fortran versions of the current millennium. +!----------------------------------------------------------------------- +! + use, intrinsic :: iso_c_binding, only: c_int32_t, c_int64_t + implicit none + intrinsic shiftr, transfer, iand ! <--- for the benefit of 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 standard-conforming floating point numbers, non-finite numbers +! 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. This will work with IEEE standard floating-point types +! besides 32 and 64, but we don't use those other types. +!----------------------------------------------------------------------- +! + is_not_finite = iand(shiftr(transfer(val,check),shift),check)==check +! + end function is_not_finite +! +!----------------------------------------------------------------------- + !----------------------------------------------------------------------- ! subroutine setup_regional_BC(Atm & @@ -263,7 +332,6 @@ subroutine setup_regional_BC(Atm & !----------------------------------------------------------------------- !*** Regional boundary data is obtained from the external BC file. !----------------------------------------------------------------------- - use, intrinsic :: ieee_arithmetic, only: ieee_is_finite use netcdf !----------------------------------------------------------------------- implicit none @@ -1553,7 +1621,6 @@ subroutine regional_bc_data(Atm,bc_hour & !*** Regional boundary data is obtained from the external BC file. !----------------------------------------------------------------------- !----------------------------------------------------------------------- - use, intrinsic :: ieee_arithmetic, only: ieee_is_finite implicit none !----------------------------------------------------------------------- ! @@ -3577,25 +3644,53 @@ subroutine remap_scalar_nggps_regional_bc(Atm & !--------------------------------------------------------------------------------- ! if(side=='west')then - do j=jsv,jev + do j=jss,jes psc(isv,j) = psc(iss,j) enddo endif ! if(side=='east')then - do j=jsv,jev + do j=jss,jes psc(iev,j) = psc(ies,j) enddo endif ! if(side=='south')then - do i=isv,iev + ! Special handling of corners + if(west_bc) then + is=iss + psc(isv,jsv) = psc(iss,jss) + else + is=isv + endif + if(east_bc) then + ie=ies + psc(iev,jsv) = psc(ies,jss) + else + ie=iev + endif + ! The rest of the south boundary + do i=is,ie psc(i,jsv) = psc(i,jss) enddo endif ! if(side=='north')then - do i=isv,iev + ! Special handling of corners + if(west_bc) then + is=iss + psc(isv,jev) = psc(iss,jes) + else + is=isv + endif + if(east_bc) then + ie=ies + psc(iev,jev) = psc(ies,jes) + else + ie=iev + endif + ! The rest of the north boundary + do i=is,ie psc(i,jev) = psc(i,jes) enddo endif @@ -3663,25 +3758,53 @@ subroutine remap_scalar_nggps_regional_bc(Atm & !--------------------------------------------------------------------------------- ! if(side=='west')then - do j=jsv,jev + do j=jss,jes ps(isv,j) = ps(iss,j) enddo endif ! if(side=='east')then - do j=jsv,jev + do j=jss,jes ps(iev,j) = ps(ies,j) enddo endif ! if(side=='south')then - do i=isv,iev + ! Special handling of corners + if(west_bc) then + is=iss + ps(isv,jsv) = ps(iss,jss) + else + is=isv + endif + if(east_bc) then + ie=ies + ps(iev,jsv) = ps(ies,jss) + else + ie=iev + endif + ! The rest of the south boundary + do i=is,ie ps(i,jsv) = ps(i,jss) enddo endif ! if(side=='north')then - do i=isv,iev + ! Special handling of corners + if(west_bc) then + is=iss + ps(isv,jev) = ps(iss,jes) + else + is=isv + endif + if(east_bc) then + ie=ies + ps(iev,jev) = ps(ies,jes) + else + ie=iev + endif + ! The rest of the north boundary + do i=is,ie ps(i,jev) = ps(i,jes) enddo endif @@ -4855,7 +4978,6 @@ subroutine bc_time_interpolation(array & !*** forecast time that is within the interval bracketed by the !*** two current boundary region states. !--------------------------------------------------------------------- - use, intrinsic :: ieee_arithmetic, only: ieee_is_finite implicit none !--------------------------------------------------------------------- ! @@ -4962,8 +5084,8 @@ 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(.not.ieee_is_finite(array(i,j,k))) then - cycle ! Bad data from caller + 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. @@ -4984,8 +5106,8 @@ subroutine bc_time_interpolation(array & do j=j1_blend,j2_blend do i=i1_blend,i2_blend ! - if(.not.ieee_is_finite(array(i,j,k))) then - cycle ! Bad data from caller + 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. @@ -5008,8 +5130,8 @@ subroutine bc_time_interpolation(array & do j=j1_blend,j2_blend do i=i1_blend,i2_blend ! - if(.not.ieee_is_finite(array(i,j,k))) then - cycle ! Bad data from caller + 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. @@ -5032,8 +5154,8 @@ 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(.not.ieee_is_finite(array(i,j,k))) then - cycle ! Bad data from caller + 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. From 6d04c313e37b71b898dbfc895ffce6f23c30e904 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Wed, 21 Sep 2022 15:24:59 +0000 Subject: [PATCH 06/10] remove -9999999 commented-out code --- model/fv_regional_bc.F90 | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index b39651364..4e3d9220d 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -1980,16 +1980,6 @@ subroutine regional_bc_data(Atm,bc_hour & allocate(ps_reg(is_input:ie_input,js_input:je_input)) !<-- Sfc pressure in domain's boundary region derived from BC files ! !----------------------------------------------------------------------- -!*** Whoever did this intentionally added a bug to the code and hid it, -!*** to conceal fixing other bugs elsewhere. This produced gibberish -!*** initial surface pressure, and made it impossible to test RRFS in -!*** debug mode for a long time. I'm retaining this code, commented out, -!*** as a reminder to others of what NOT to do. -!----------------------------------------------------------------------- -! - !ps_reg=-9999999 ! for now don't set to snan until remap dwinds is changed -! -!----------------------------------------------------------------------- !*** Initialize ps_reg to real_snan to detect errors later in code. !----------------------------------------------------------------------- ! From 229dbe5b6139733be924c42f9cee2c3520761268 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Wed, 21 Sep 2022 16:04:43 +0000 Subject: [PATCH 07/10] quiet the NaNs passed to Atmp%ps --- model/fv_regional_bc.F90 | 39 +++++++++++++++++++++++++++++---------- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index 4e3d9220d..e29ce7dd8 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -236,11 +236,14 @@ module fv_regional_mod ,oro_data ='oro_data.tile7.halo4.nc' #ifdef OVERLOAD_R4 - real, parameter:: real_snan=real(Z'FFBFFFFF') + real, parameter:: real_snan=real(Z'FFBFFFFF') ! <-- One of many ways to express 32-bit signalling NaN + real, parameter:: real_qnan=real(Z'FFFFFFFF') ! <-- One of many ways to express 32-bit quiet NaN #else - real, parameter:: real_snan=real(Z'FFF7FFFFFFFFFFFF') + real, parameter:: real_snan=real(Z'FFF7FFFFFFFFFFFF') ! <-- One of many ways to express 64-bit signalling NaN + real, parameter:: real_qnan=real(Z'FFFFFFFFFFFFFFFF') ! <-- One of many ways to express 64-bit quiet NaN #endif - real(kind=R_GRID), parameter:: dbl_snan=real(Z'FFF7FFFFFFFFFFFF',kind=R_GRID) + real(kind=R_GRID), parameter:: dbl_snan=real(Z'FFF7FFFFFFFFFFFF',kind=R_GRID) ! <-- One of many ways to express 64-bit signalling NaN + real(kind=R_GRID), parameter:: dbl_qnan=real(Z'FFFFFFFFFFFFFFFF',kind=R_GRID) ! <-- One of many ways to express 64-bit quiet NaN interface dump_field module procedure dump_field_3d @@ -260,12 +263,13 @@ logical function is_not_finite(val) ! !----------------------------------------------------------------------- !*** This routine is equivalent to ".not. ieee_is_finite(val)" -!*** which does not work with signaling NaN in gfortran < version 12. -!*** It has the added advantage of being easily inlined and vectorized. -!*** It will detect all possible infinite and not-a-number (NaN) values -!*** in standard-conforming 32 bit and 64 bit floating-point. This is -!*** used to detect access beyond boundary regions in the blending code, -!*** which has no control over how the caller stores data in memory. +!*** which does not work with signaling Not a Number (NaN) in gfortran +!*** before version 12. It has the added advantage of being easily +!*** inlined and vectorized. It will detect all possible infinite and +!*** not-a-number (NaN) values in standard-conforming 32 bit and 64 bit +!*** floating-point. This is used to detect access beyond boundary +!*** regions in the blending code, which has no control over how the +!*** caller stores data in memory. !*** !*** The gfortran bug is documented here: !*** @@ -310,7 +314,7 @@ logical function is_not_finite(val) #endif ! !----------------------------------------------------------------------- -! For standard-conforming floating point numbers, non-finite numbers +! For standard-conforming 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. This will work with IEEE standard floating-point types @@ -2225,6 +2229,21 @@ subroutine regional_bc_data(Atm,bc_hour & !----------------------------------------------------------------------- enddo sides_scalars !----------------------------------------------------------------------- + +! +!----------------------------------------------------------------------- +!*** The caller may eventually print or write Atm%ps to disk, so it +!*** cannot contain signalling NaN. Here, we replace all non-finite +!*** data with quiet NaN instead. +!----------------------------------------------------------------------- +! + do j=lbound(Atm%ps,2),ubound(Atm%ps,2) + do i=lbound(Atm%ps,1),ubound(Atm%ps,1) + if(is_not_finite(Atm%ps(i,j))) then + Atm%ps(i,j) = real_qnan + endif + enddo + enddo ! !----------------------------------------------------------------------- !*** Now that we have the pressure throughout the boundary region From b923e9f442d2fc6037e9ff781aead3855c878a52 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Wed, 21 Sep 2022 18:25:52 +0000 Subject: [PATCH 08/10] simplify comments and explain snan --- model/fv_regional_bc.F90 | 62 +++++++++++++++++++++++----------------- 1 file changed, 36 insertions(+), 26 deletions(-) diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index e29ce7dd8..5cbda84e5 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -235,6 +235,29 @@ module fv_regional_mod character(len=100) :: grid_data='grid.tile7.halo4.nc' & ,oro_data ='oro_data.tile7.halo4.nc' +! +!----------------------------------------------------------------------- +!*** All arrays and variables that contain boundary data should be +!*** initialized to real_snan or dbl_snan (as appropriate) because that +!*** is the only way to detect errors in the boundary generation. +!*** Usual approaches like detecting out-of-bounds accesses won't work +!*** because the arrays may be larger than the boundaries. +!*** +!*** Regions of the arrays that should have valid boundary data should +!*** never contain a NaN. If they do, then there is an error that has +!*** not yet been identified. In particular, the blending code has an +!*** unidentified problem. +!*** +!*** Presently, the blending code is asked to blend data in regions +!*** that contain NaN. This means either: +!*** +!*** 1) not all boundary areas are filled in, OR +!*** 2) the blending code is accessing beyond the boundaries +!*** +!*** The reasons for this are unknown, but the workaround does work: +!*** detect the NaNs and don't use those points in blending. +!----------------------------------------------------------------------- +! #ifdef OVERLOAD_R4 real, parameter:: real_snan=real(Z'FFBFFFFF') ! <-- One of many ways to express 32-bit signalling NaN real, parameter:: real_qnan=real(Z'FFFFFFFF') ! <-- One of many ways to express 32-bit quiet NaN @@ -263,28 +286,22 @@ logical function is_not_finite(val) ! !----------------------------------------------------------------------- !*** This routine is equivalent to ".not. ieee_is_finite(val)" -!*** which does not work with signaling Not a Number (NaN) in gfortran -!*** before version 12. It has the added advantage of being easily -!*** inlined and vectorized. It will detect all possible infinite and -!*** not-a-number (NaN) values in standard-conforming 32 bit and 64 bit -!*** floating-point. This is used to detect access beyond boundary -!*** regions in the blending code, which has no control over how the -!*** caller stores data in memory. -!*** -!*** The gfortran bug is documented here: +!*** 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 !*** -!*** A sufficiently intelligent compiler will optimize this function -!*** away, into three assembly instructions that can be vectorized. -!*** -!*** Technically, this routine needs Fortran 2008 due to shiftr() but -!*** that intrinsic is widely implemented in compilers that lack -!*** support for Fortran versions of the current millennium. +!*** 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 ! <--- for the benefit of older compilers ! !----------------------------------------------------------------------- @@ -314,11 +331,9 @@ logical function is_not_finite(val) #endif ! !----------------------------------------------------------------------- -! For standard-conforming 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. This will work with IEEE standard floating-point types -! besides 32 and 64, but we don't use those other types. +! 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 @@ -1982,12 +1997,7 @@ subroutine regional_bc_data(Atm,bc_hour & !----------------------------------------------------------------------- ! allocate(ps_reg(is_input:ie_input,js_input:je_input)) !<-- Sfc pressure in domain's boundary region derived from BC files -! -!----------------------------------------------------------------------- -!*** Initialize ps_reg to real_snan to detect errors later in code. -!----------------------------------------------------------------------- -! - ps_reg=real_snan + ps_reg=real_snan !<-- detect access of uninitialized pressures ! !----------------------------------------------------------------------- !*** We have the boundary variables from the BC file on the levels From 68a833fe0bf502737bda03cf32ba4ee41ba25464 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Tue, 27 Sep 2022 12:45:04 +0000 Subject: [PATCH 09/10] use i-1 & j-1 for two-point averages, when available --- model/fv_regional_bc.F90 | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index 5cbda84e5..4443d04ca 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -256,6 +256,10 @@ module fv_regional_mod !*** !*** The reasons for this are unknown, but the workaround does work: !*** detect the NaNs and don't use those points in blending. +!*** +!*** Also, the wind remapping code needs surface pressure beyond the +!*** wind boundary regions, but the surface pressure boundary data +!*** doesn't extend that far. !----------------------------------------------------------------------- ! #ifdef OVERLOAD_R4 @@ -4173,12 +4177,22 @@ subroutine remap_dwinds_regional_bc(Atm & !------ do k=1,km+1 do i=is_u,ie_u - pe0(i,k) = ak0(k) + bk0(k)*psc(i,j) + if(is_not_finite(psc(i,j-1))) then + ! Workaround for bug: PS not available in some velocity areas + pe0(i,k) = ak0(k) + bk0(k)*psc(i,j) + else + pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i,j-1)+psc(i,j)) + endif enddo enddo do k=1,npz+1 do i=is_u,ie_u - pe1(i,k) = Atm%ak(k) + Atm%bk(k)*psc(i,j) + if(is_not_finite(psc(i,j-1))) then + ! Workaround for bug: PS not available in some velocity areas + pe1(i,k) = Atm%ak(k) + Atm%bk(k)*psc(i,j) + else + pe1(i,k) = Atm%ak(k) + Atm%bk(k)*0.5*(psc(i,j-1)+psc(i,j)) + endif enddo enddo call mappm(km, pe0(is_u:ie_u,1:km+1), ud(is_u:ie_u,j,1:km), npz, pe1(is_u:ie_u,1:npz+1), & @@ -4214,12 +4228,22 @@ subroutine remap_dwinds_regional_bc(Atm & do k=1,km+1 do i=is_v,ie_v - pe0(i,k) = ak0(k) + bk0(k)*psc(i,j) + if(is_not_finite(psc(i-1,j))) then + ! Workaround for bug: PS not available in some velocity areas + pe0(i,k) = ak0(k) + bk0(k)*psc(i,j) + else + pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i-1,j)+psc(i,j)) + endif enddo enddo do k=1,npz+1 do i=is_v,ie_v - pe1(i,k) = Atm%ak(k) + Atm%bk(k)*psc(i,j) + if(is_not_finite(psc(i-1,j))) then + ! Workaround for bug: PS not available in some velocity areas + pe1(i,k) = Atm%ak(k) + Atm%bk(k)*psc(i,j) + else + pe1(i,k) = Atm%ak(k) + Atm%bk(k)*0.5*(psc(i-1,j)+psc(i,j)) + endif enddo enddo call mappm(km, pe0(is_v:ie_v,1:km+1), vd(is_v:ie_v,j,1:km), npz, pe1(is_v:ie_v,1:npz+1), & From 9824ef37504e23f61eda762f9cdce5d01cd1ef2b Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Mon, 3 Oct 2022 19:01:43 +0000 Subject: [PATCH 10/10] Replace many changes with PR #220 --- model/fv_regional_bc.F90 | 730 ++++++++++++++------------------------- 1 file changed, 268 insertions(+), 462 deletions(-) diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90 index 4443d04ca..b7c25e3a3 100644 --- a/model/fv_regional_bc.F90 +++ b/model/fv_regional_bc.F90 @@ -235,42 +235,12 @@ module fv_regional_mod character(len=100) :: grid_data='grid.tile7.halo4.nc' & ,oro_data ='oro_data.tile7.halo4.nc' -! -!----------------------------------------------------------------------- -!*** All arrays and variables that contain boundary data should be -!*** initialized to real_snan or dbl_snan (as appropriate) because that -!*** is the only way to detect errors in the boundary generation. -!*** Usual approaches like detecting out-of-bounds accesses won't work -!*** because the arrays may be larger than the boundaries. -!*** -!*** Regions of the arrays that should have valid boundary data should -!*** never contain a NaN. If they do, then there is an error that has -!*** not yet been identified. In particular, the blending code has an -!*** unidentified problem. -!*** -!*** Presently, the blending code is asked to blend data in regions -!*** that contain NaN. This means either: -!*** -!*** 1) not all boundary areas are filled in, OR -!*** 2) the blending code is accessing beyond the boundaries -!*** -!*** The reasons for this are unknown, but the workaround does work: -!*** detect the NaNs and don't use those points in blending. -!*** -!*** Also, the wind remapping code needs surface pressure beyond the -!*** wind boundary regions, but the surface pressure boundary data -!*** doesn't extend that far. -!----------------------------------------------------------------------- -! #ifdef OVERLOAD_R4 - real, parameter:: real_snan=real(Z'FFBFFFFF') ! <-- One of many ways to express 32-bit signalling NaN - real, parameter:: real_qnan=real(Z'FFFFFFFF') ! <-- One of many ways to express 32-bit quiet NaN + real, parameter:: real_snan=real(Z'FFBFFFFF') #else - real, parameter:: real_snan=real(Z'FFF7FFFFFFFFFFFF') ! <-- One of many ways to express 64-bit signalling NaN - real, parameter:: real_qnan=real(Z'FFFFFFFFFFFFFFFF') ! <-- One of many ways to express 64-bit quiet NaN + real, parameter:: real_snan=real(Z'FFF7FFFFFFFFFFFF') #endif - real(kind=R_GRID), parameter:: dbl_snan=real(Z'FFF7FFFFFFFFFFFF',kind=R_GRID) ! <-- One of many ways to express 64-bit signalling NaN - real(kind=R_GRID), parameter:: dbl_qnan=real(Z'FFFFFFFFFFFFFFFF',kind=R_GRID) ! <-- One of many ways to express 64-bit quiet NaN + real(kind=R_GRID), parameter:: dbl_snan=real(Z'FFF7FFFFFFFFFFFF',kind=R_GRID) interface dump_field module procedure dump_field_3d @@ -284,6 +254,7 @@ module fv_regional_mod logical :: data_source_fv3gfs contains + !----------------------------------------------------------------------- ! logical function is_not_finite(val) @@ -306,7 +277,7 @@ logical function is_not_finite(val) ! supported in older compilers. !----------------------------------------------------------------------- ! - intrinsic shiftr, transfer, iand ! <--- for the benefit of older compilers + intrinsic shiftr, transfer, iand ! <--- declare intrinsic to help older compilers ! !----------------------------------------------------------------------- ! Use value-based argument passing instead of reference-based to avoid @@ -344,8 +315,6 @@ logical function is_not_finite(val) ! end function is_not_finite ! -!----------------------------------------------------------------------- - !----------------------------------------------------------------------- ! subroutine setup_regional_BC(Atm & @@ -1043,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 @@ -1060,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 @@ -1120,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 @@ -2035,7 +2004,16 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif ! - if(nside==2)then + if(nside==2)then + if(south_bc)then + call_remap=.true. + side='south' + bc_side_t1=>BC_t1%south + bc_side_t0=>BC_t0%south + endif + endif +! + if(nside==3)then if(east_bc)then call_remap=.true. side='east ' @@ -2044,7 +2022,7 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif ! - if(nside==3)then + if(nside==4)then if(west_bc)then call_remap=.true. side='west ' @@ -2052,15 +2030,6 @@ subroutine regional_bc_data(Atm,bc_hour & bc_side_t0=>BC_t0%west endif endif -! - if(nside==4)then - if(south_bc)then - call_remap=.true. - side='south' - bc_side_t1=>BC_t1%south - bc_side_t0=>BC_t0%south - endif - endif ! if(call_remap)then call remap_scalar_nggps_regional_bc(Atm & @@ -2138,9 +2107,57 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif + if(nside==2)then + if(south_bc)then + if (is == 1) then + istart = 1 + else + istart = isd + endif + if (ie == npx-1) then + iend = npx-1 + else + iend = ied + endif + + do k=1,npz + do j=npy,jed + do i=istart,iend + delz_regBC%north_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) + delz_regBC%north_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) + enddo + enddo + enddo + + ! North, south include all corners + if (is == 1) then + do k=1,npz + do j=npy,jed + do i=isd,0 + delz_regBC%west_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) + delz_regBC%west_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) + enddo + enddo + enddo + endif + + + if (ie == npx-1) then + do k=1,npz + do j=npy,jed + do i=npx,ied + delz_regBC%east_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) + delz_regBC%east_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) + enddo + enddo + enddo + endif + endif + endif + ! - if(nside==2)then + if(nside==3)then if(east_bc)then if (js == 1) then jstart = 1 @@ -2165,7 +2182,7 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif - if(nside==3)then + if(nside==4)then if(west_bc)then if (js == 1) then jstart = 1 @@ -2190,74 +2207,11 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif - if(nside==4)then - if(south_bc)then - if (is == 1) then - istart = 1 - else - istart = isd - endif - if (ie == npx-1) then - iend = npx-1 - else - iend = ied - endif - - do k=1,npz - do j=npy,jed - do i=istart,iend - delz_regBC%north_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) - delz_regBC%north_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) - enddo - enddo - enddo - - ! North, south include all corners - if (is == 1) then - do k=1,npz - do j=npy,jed - do i=isd,0 - delz_regBC%west_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) - delz_regBC%west_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) - enddo - enddo - enddo - endif - - - if (ie == npx-1) then - do k=1,npz - do j=npy,jed - do i=npx,ied - delz_regBC%east_t1(i,j,k) = bc_side_t1%delz_BC(i,j,k) - delz_regBC%east_t0(i,j,k) = bc_side_t0%delz_BC(i,j,k) - enddo - enddo - enddo - endif - endif - endif - endif ! !----------------------------------------------------------------------- enddo sides_scalars !----------------------------------------------------------------------- - -! -!----------------------------------------------------------------------- -!*** The caller may eventually print or write Atm%ps to disk, so it -!*** cannot contain signalling NaN. Here, we replace all non-finite -!*** data with quiet NaN instead. -!----------------------------------------------------------------------- -! - do j=lbound(Atm%ps,2),ubound(Atm%ps,2) - do i=lbound(Atm%ps,1),ubound(Atm%ps,1) - if(is_not_finite(Atm%ps(i,j))) then - Atm%ps(i,j) = real_qnan - endif - enddo - enddo ! !----------------------------------------------------------------------- !*** Now that we have the pressure throughout the boundary region @@ -2333,6 +2287,23 @@ subroutine regional_bc_data(Atm,bc_hour & endif ! if(nside==2)then + if(south_bc)then + call_remap=.true. + bc_side_t1=>BC_t1%south +! + is_u=Atm%regional_bc_bounds%is_south_uvs + ie_u=Atm%regional_bc_bounds%ie_south_uvs + js_u=Atm%regional_bc_bounds%js_south_uvs + je_u=Atm%regional_bc_bounds%je_south_uvs +! + is_v=Atm%regional_bc_bounds%is_south_uvw + ie_v=Atm%regional_bc_bounds%ie_south_uvw + js_v=Atm%regional_bc_bounds%js_south_uvw + je_v=Atm%regional_bc_bounds%je_south_uvw + endif + endif +! + if(nside==3)then if(east_bc)then call_remap=.true. bc_side_t1=>BC_t1%east @@ -2349,7 +2320,7 @@ subroutine regional_bc_data(Atm,bc_hour & endif endif ! - if(nside==3)then + if(nside==4)then if(west_bc)then call_remap=.true. bc_side_t1=>BC_t1%west @@ -2365,23 +2336,6 @@ subroutine regional_bc_data(Atm,bc_hour & je_v=Atm%regional_bc_bounds%je_west_uvw endif endif -! - if(nside==4)then - if(south_bc)then - call_remap=.true. - bc_side_t1=>BC_t1%south -! - is_u=Atm%regional_bc_bounds%is_south_uvs - ie_u=Atm%regional_bc_bounds%ie_south_uvs - js_u=Atm%regional_bc_bounds%js_south_uvs - je_u=Atm%regional_bc_bounds%je_south_uvs -! - is_v=Atm%regional_bc_bounds%is_south_uvw - ie_v=Atm%regional_bc_bounds%ie_south_uvw - js_v=Atm%regional_bc_bounds%js_south_uvw - je_v=Atm%regional_bc_bounds%je_south_uvw - endif - endif ! if(call_remap)then ! @@ -2875,6 +2829,17 @@ subroutine fill_divgd_BC endif ! if(nside==2)then + if(south_bc)then + call_set=.true. + bc_side_t1=>BC_t1%south + is0=lbound(BC_t1%south%divgd_BC,1) + ie0=ubound(BC_t1%south%divgd_BC,1) + js0=lbound(BC_t1%south%divgd_BC,2) + je0=ubound(BC_t1%south%divgd_BC,2) + endif + endif +! + if(nside==3)then if(east_bc)then call_set=.true. bc_side_t1=>BC_t1%east @@ -2885,7 +2850,7 @@ subroutine fill_divgd_BC endif endif ! - if(nside==3)then + if(nside==4)then if(west_bc)then call_set=.true. bc_side_t1=>BC_t1%west @@ -2895,17 +2860,6 @@ subroutine fill_divgd_BC je0=ubound(BC_t1%west%divgd_BC,2) endif endif -! - if(nside==4)then - if(south_bc)then - call_set=.true. - bc_side_t1=>BC_t1%south - is0=lbound(BC_t1%south%divgd_BC,1) - ie0=ubound(BC_t1%south%divgd_BC,1) - js0=lbound(BC_t1%south%divgd_BC,2) - je0=ubound(BC_t1%south%divgd_BC,2) - endif - endif ! if(call_set)then do k=1,klev_out @@ -2968,6 +2922,17 @@ subroutine fill_q_con_BC endif ! if(nside==2)then + if(south_bc)then + call_set=.true. + bc_side_t1=>BC_t1%south + is0=lbound(BC_t1%south%q_con_BC,1) + ie0=ubound(BC_t1%south%q_con_BC,1) + js0=lbound(BC_t1%south%q_con_BC,2) + je0=ubound(BC_t1%south%q_con_BC,2) + endif + endif +! + if(nside==3)then if(east_bc)then call_set=.true. bc_side_t1=>BC_t1%east @@ -2978,7 +2943,7 @@ subroutine fill_q_con_BC endif endif ! - if(nside==3)then + if(nside==4)then if(west_bc)then call_set=.true. bc_side_t1=>BC_t1%west @@ -2988,17 +2953,6 @@ subroutine fill_q_con_BC je0=ubound(BC_t1%west%q_con_BC,2) endif endif -! - if(nside==4)then - if(south_bc)then - call_set=.true. - bc_side_t1=>BC_t1%south - is0=lbound(BC_t1%south%q_con_BC,1) - ie0=ubound(BC_t1%south%q_con_BC,1) - js0=lbound(BC_t1%south%q_con_BC,2) - je0=ubound(BC_t1%south%q_con_BC,2) - endif - endif ! if(call_set)then do k=1,klev_out @@ -3056,23 +3010,23 @@ subroutine fill_cappa_BC endif ! if(nside==2)then - if(east_bc)then + if(south_bc)then call_compute=.true. - bc_side_t1=>BC_t1%east + bc_side_t1=>BC_t1%south endif endif ! if(nside==3)then - if(west_bc)then + if(east_bc)then call_compute=.true. - bc_side_t1=>BC_t1%west + bc_side_t1=>BC_t1%east endif endif ! if(nside==4)then - if(south_bc)then + if(west_bc)then call_compute=.true. - bc_side_t1=>BC_t1%south + bc_side_t1=>BC_t1%west endif endif ! @@ -3287,23 +3241,46 @@ subroutine read_regional_bc_file(is_input,ie_input & endif endif ! -!---------- -!*** East -!---------- +!----------- +!*** South +!----------- ! if(nside==2)then - if(east_bc)then + if(south_bc)then call_get_var=.true. ! - var_name=trim(var_name_root)//"_left" + var_name=trim(var_name_root)//"_top" ! - j_start_array=js_input + i_start_array=is_input + i_end_array =ie_input + j_start_array=je_input-nhalo_data-nrows_blend+1 j_end_array =je_input ! - i_start_array=is_input + i_start_data=i_start_array+nhalo_data + i_count=i_end_array-i_start_array+1 + j_start_data=1 + j_count=j_end_array-j_start_array+1 ! - if(trim(var_name_root)=='u_w'.or.trim(var_name_root)=='v_w')then - i_end_array=is_input+nhalo_data+nrows_blend + endif + endif +! +!---------- +!*** East +!---------- +! + if(nside==3)then + if(east_bc)then + call_get_var=.true. +! + var_name=trim(var_name_root)//"_left" +! + j_start_array=js_input + j_end_array =je_input +! + i_start_array=is_input +! + if(trim(var_name_root)=='u_w'.or.trim(var_name_root)=='v_w')then + i_end_array=is_input+nhalo_data+nrows_blend else i_end_array=is_input+nhalo_data+nrows_blend-1 endif @@ -3335,7 +3312,7 @@ subroutine read_regional_bc_file(is_input,ie_input & !*** West !---------- ! - if(nside==3)then + if(nside==4)then if(west_bc)then call_get_var=.true. ! @@ -3371,29 +3348,6 @@ subroutine read_regional_bc_file(is_input,ie_input & endif endif ! -!----------- -!*** South -!----------- -! - if(nside==4)then - if(south_bc)then - call_get_var=.true. -! - var_name=trim(var_name_root)//"_top" -! - i_start_array=is_input - i_end_array =ie_input - j_start_array=je_input-nhalo_data-nrows_blend+1 - j_end_array =je_input -! - i_start_data=i_start_array+nhalo_data - i_count=i_end_array-i_start_array+1 - j_start_data=1 - j_count=j_end_array-j_start_array+1 -! - endif - endif -! !----------------------------------------------------------------------- !*** Fill this task's subset of boundary data for this 3-D !*** or 4-D variable. This includes the data in the domain's @@ -3571,7 +3525,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm & ,npz & !<-- # of levels in final 3-D integration variables ,ncnst !<-- # of tracer variables real, intent(in):: ak0(km+1), bk0(km+1) - real, intent(inout), dimension(is_bc:ie_bc,js_bc:je_bc):: psc + real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc):: psc real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km):: t_in real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km):: omga real, intent(in), dimension(is_bc:ie_bc,js_bc:je_bc,km,ncnst):: qa @@ -3599,7 +3553,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm & real(kind=R_GRID), dimension(is_bc:ie_bc,km+1):: pn0 real(kind=R_GRID):: pst !!! High-precision - integer i,ie,is,j,je,js,k,l,m, k2,iq, iss,ies,jss,jes, isv,iev,jsv,jev + integer i,ie,is,j,je,js,k,l,m, k2,iq integer sphum, o3mr, liq_wat, ice_wat, rainwat, snowwat, graupel, hailwat, cld_amt ! !--------------------------------------------------------------------------------- @@ -3633,95 +3587,25 @@ subroutine remap_scalar_nggps_regional_bc(Atm & !*** This is needed to obtain pressures that will surround the wind points. !--------------------------------------------------------------------------------- ! - isv=is_bc - iss=isv - if(side=='west')then - isv=ie_bc-nhalo_data-nrows_blend - iss=isv+1 - endif -! - iev=ie_bc - ies=iev - if(side=='east')then - iev=is_bc+nhalo_data+nrows_blend - ies=iev-1 - endif -! - jsv=js_bc - jss=jsv - if(side=='south')then - jsv=je_bc-nhalo_data-nrows_blend - jss=jsv+1 - endif -! - jev=je_bc - jes=jev - if(side=='north')then - jev=js_bc+nhalo_data+nrows_blend - jes=jev-1 - endif -! -!--------------------------------------------------------------------------------- -!*** Workaround for incomplete 'ps' in boundary condition files. Copy from -!*** scalar halo regions (s) to velocity halo regions (v) beyond them. -!--------------------------------------------------------------------------------- -! + is=is_bc if(side=='west')then - do j=jss,jes - psc(isv,j) = psc(iss,j) - enddo + is=ie_bc-nhalo_data-nrows_blend+1 endif ! + ie=ie_bc if(side=='east')then - do j=jss,jes - psc(iev,j) = psc(ies,j) - enddo + ie=is_bc+nhalo_data+nrows_blend-1 endif ! + js=js_bc if(side=='south')then - ! Special handling of corners - if(west_bc) then - is=iss - psc(isv,jsv) = psc(iss,jss) - else - is=isv - endif - if(east_bc) then - ie=ies - psc(iev,jsv) = psc(ies,jss) - else - ie=iev - endif - ! The rest of the south boundary - do i=is,ie - psc(i,jsv) = psc(i,jss) - enddo + js=je_bc-nhalo_data-nrows_blend+1 endif ! + je=je_bc if(side=='north')then - ! Special handling of corners - if(west_bc) then - is=iss - psc(isv,jev) = psc(iss,jes) - else - is=isv - endif - if(east_bc) then - ie=ies - psc(iev,jev) = psc(ies,jes) - else - ie=iev - endif - ! The rest of the north boundary - do i=is,ie - psc(i,jev) = psc(i,jes) - enddo + je=js_bc+nhalo_data+nrows_blend-1 endif -! - is = iss - js = jss - ie = ies - je = jes ! ! Ensure uninitialized memory isn't used pn0 = real_snan @@ -3729,7 +3613,6 @@ subroutine remap_scalar_nggps_regional_bc(Atm & 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 @@ -3774,73 +3657,16 @@ subroutine remap_scalar_nggps_regional_bc(Atm & !--------------------------------------------------------------------------------- enddo jloop1 !--------------------------------------------------------------------------------- -! -!--------------------------------------------------------------------------------- -!*** Workaround for incomplete 'ps' in boundary condition files. Copy from -!*** scalar halo regions (s) to velocity halo regions (v) beyond them. -!--------------------------------------------------------------------------------- -! - if(side=='west')then - do j=jss,jes - ps(isv,j) = ps(iss,j) - enddo - endif -! - if(side=='east')then - do j=jss,jes - ps(iev,j) = ps(ies,j) - enddo - endif -! - if(side=='south')then - ! Special handling of corners - if(west_bc) then - is=iss - ps(isv,jsv) = ps(iss,jss) - else - is=isv - endif - if(east_bc) then - ie=ies - ps(iev,jsv) = ps(ies,jss) - else - ie=iev - endif - ! The rest of the south boundary - do i=is,ie - ps(i,jsv) = ps(i,jss) - enddo - endif -! - if(side=='north')then - ! Special handling of corners - if(west_bc) then - is=iss - ps(isv,jev) = ps(iss,jes) - else - is=isv - endif - if(east_bc) then - ie=ies - ps(iev,jev) = ps(ies,jes) - else - ie=iev - endif - ! The rest of the north boundary - do i=is,ie - ps(i,jev) = ps(i,jes) - enddo - endif !--------------------------------------------------------------------------------- !*** Transfer values from the expanded boundary array for sfc pressure into !*** the Atm object. !--------------------------------------------------------------------------------- ! - is=min(ubound(Atm%ps,1),max(lbound(Atm%ps,1),isv)) - ie=min(ubound(Atm%ps,1),max(lbound(Atm%ps,1),iev)) - js=min(ubound(Atm%ps,2),max(lbound(Atm%ps,2),jsv)) - je=min(ubound(Atm%ps,2),max(lbound(Atm%ps,2),jev)) + 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 @@ -4177,22 +4003,12 @@ subroutine remap_dwinds_regional_bc(Atm & !------ do k=1,km+1 do i=is_u,ie_u - if(is_not_finite(psc(i,j-1))) then - ! Workaround for bug: PS not available in some velocity areas - pe0(i,k) = ak0(k) + bk0(k)*psc(i,j) - else - pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i,j-1)+psc(i,j)) - endif + pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i,j-1)+psc(i,j)) enddo enddo do k=1,npz+1 do i=is_u,ie_u - if(is_not_finite(psc(i,j-1))) then - ! Workaround for bug: PS not available in some velocity areas - pe1(i,k) = Atm%ak(k) + Atm%bk(k)*psc(i,j) - else - pe1(i,k) = Atm%ak(k) + Atm%bk(k)*0.5*(psc(i,j-1)+psc(i,j)) - endif + pe1(i,k) = Atm%ak(k) + Atm%bk(k)*0.5*(psc(i,j-1)+psc(i,j)) enddo enddo call mappm(km, pe0(is_u:ie_u,1:km+1), ud(is_u:ie_u,j,1:km), npz, pe1(is_u:ie_u,1:npz+1), & @@ -4228,22 +4044,12 @@ subroutine remap_dwinds_regional_bc(Atm & do k=1,km+1 do i=is_v,ie_v - if(is_not_finite(psc(i-1,j))) then - ! Workaround for bug: PS not available in some velocity areas - pe0(i,k) = ak0(k) + bk0(k)*psc(i,j) - else - pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i-1,j)+psc(i,j)) - endif + pe0(i,k) = ak0(k) + bk0(k)*0.5*(psc(i-1,j)+psc(i,j)) enddo enddo do k=1,npz+1 do i=is_v,ie_v - if(is_not_finite(psc(i-1,j))) then - ! Workaround for bug: PS not available in some velocity areas - pe1(i,k) = Atm%ak(k) + Atm%bk(k)*psc(i,j) - else - pe1(i,k) = Atm%ak(k) + Atm%bk(k)*0.5*(psc(i-1,j)+psc(i,j)) - endif + pe1(i,k) = Atm%ak(k) + Atm%bk(k)*0.5*(psc(i-1,j)+psc(i,j)) enddo enddo call mappm(km, pe0(is_v:ie_v,1:km+1), vd(is_v:ie_v,j,1:km), npz, pe1(is_v:ie_v,1:npz+1), & @@ -4694,11 +4500,61 @@ subroutine regional_boundary_update(array & endif endif ! +!----------- +!*** South +!----------- +! + if(nside==2)then + if(south_bc)then + call_interp=.true. + bc_side_t0=>bc_south_t0 + bc_side_t1=>bc_south_t1 +! + i1=isd + i2=ied + if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then + i2=ied+1 + endif +! + j1=je+1 + j2=jed + if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then + j1=je+2 + j2=jed+1 + endif +! + i1_blend=is + i2_blend=ie + IF ( east_bc ) THEN + i1_blend=is + ELSE + i1_blend=isd !is-nhalo_model + ENDIF + IF ( west_bc ) THEN + i2_blend=ie + ELSE + i2_blend=ied ! ie+nhalo_model + ENDIF + if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then +! i2_blend=ie+1 + i2_blend=i2_blend+1 + endif + j2_blend=je + if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then + j2_blend=je+1 + endif + j1_blend=j2_blend-nrows_blend_user+1 + i_bc=-9e9 + j_bc=j1 +! + endif + endif +! !---------- !*** East !---------- ! - if(nside==2)then + if(nside==3)then if(east_bc)then call_interp=.true. bc_side_t0=>bc_east_t0 @@ -4754,7 +4610,7 @@ subroutine regional_boundary_update(array & !*** West !---------- ! - if(nside==3)then + if(nside==4)then if(west_bc)then call_interp=.true. bc_side_t0=>bc_west_t0 @@ -4809,56 +4665,6 @@ subroutine regional_boundary_update(array & endif endif ! -!----------- -!*** South -!----------- -! - if(nside==4)then - if(south_bc)then - call_interp=.true. - bc_side_t0=>bc_south_t0 - bc_side_t1=>bc_south_t1 -! - i1=isd - i2=ied - if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then - i2=ied+1 - endif -! - j1=je+1 - j2=jed - if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then - j1=je+2 - j2=jed+1 - endif -! - i1_blend=is - i2_blend=ie - IF ( east_bc ) THEN - i1_blend=is - ELSE - i1_blend=isd !is-nhalo_model - ENDIF - IF ( west_bc ) THEN - i2_blend=ie - ELSE - i2_blend=ied ! ie+nhalo_model - ENDIF - if(trim(bc_vbl_name)=='uc'.or.trim(bc_vbl_name)=='v'.or.trim(bc_vbl_name)=='divgd')then -! i2_blend=ie+1 - i2_blend=i2_blend+1 - endif - j2_blend=je - if(trim(bc_vbl_name)=='u'.or.trim(bc_vbl_name)=='vc'.or.trim(bc_vbl_name)=='divgd')then - j2_blend=je+1 - endif - j1_blend=j2_blend-nrows_blend_user+1 - i_bc=-9e9 - j_bc=j1 -! - endif - endif -! !--------------------------------------------------------------------- !*** Get the pointers pointing at the boundary arrays holding the !*** two time levels of the given prognostic array's boundary region @@ -5139,24 +4945,21 @@ subroutine bc_time_interpolation(array & enddo endif ! -!---------- -!*** East -!---------- +!----------- +!*** South +!----------- ! - if(nside==2.and.east_bc)then - rdenom=1./real(Max(1,i2_blend-i_bc-1)) + if(nside==2.and.south_bc)then + rdenom=1./real(Max(1,j_bc-j1_blend-1)) do k=1,ubnd_z 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. -! - factor_dist=exp(-(blend_exp1+blend_exp2*(i-i_bc-1)*rdenom)) !<-- Exponential falloff of blending weights. -! + 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 enddo enddo @@ -5164,22 +4967,22 @@ subroutine bc_time_interpolation(array & endif ! !---------- -!*** West +!*** East !---------- ! - if(nside==3.and.west_bc)then - rdenom=1./real(Max(1, i_bc-i1_blend-1)) + if(nside==3.and.east_bc)then + rdenom=1./real(Max(1,i2_blend-i_bc-1)) 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. ! - factor_dist=exp(-(blend_exp1+blend_exp2*(i_bc-i-1)*rdenom)) !<-- Exponential falloff of blending weights. + factor_dist=exp(-(blend_exp1+blend_exp2*(i-i_bc-1)*rdenom)) !<-- Exponential falloff of blending weights. ! array(i,j,k)=(1.-factor_dist)*array(i,j,k)+factor_dist*blend_value enddo @@ -5187,21 +4990,24 @@ subroutine bc_time_interpolation(array & enddo endif ! -!----------- -!*** South -!----------- +!---------- +!*** West +!---------- ! - if(nside==4.and.south_bc)then - rdenom=1./real(Max(1,j_bc-j1_blend-1)) + if(nside==4.and.west_bc)then + rdenom=1./real(Max(1, i_bc-i1_blend-1)) do k=1,ubnd_z 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. +! + 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. +! + factor_dist=exp(-(blend_exp1+blend_exp2*(i_bc-i-1)*rdenom)) !<-- Exponential falloff of blending weights. +! array(i,j,k)=(1.-factor_dist)*array(i,j,k)+factor_dist*blend_value enddo enddo @@ -5286,6 +5092,29 @@ subroutine regional_bc_t1_to_t0(BC_t1,BC_t0 & endif ! if(nside==2)then + if(south_bc)then + move=.true. + bc_side_t0=>BC_t0%south + bc_side_t1=>BC_t1%south +! + is_c=bnds%is_south !<-- + ie_c=bnds%ie_south ! South BC index limits + js_c=bnds%js_south ! for centers of grid cells. + je_c=bnds%je_south !<-- +! + is_s=bnds%is_south_uvs !<-- + ie_s=bnds%ie_south_uvs ! South BC index limits + js_s=bnds%js_south_uvs ! for winds on N/S sides of grid cells. + je_s=bnds%je_south_uvs !<-- +! + is_w=bnds%is_south_uvw !<-- + ie_w=bnds%ie_south_uvw ! South BC index limits + js_w=bnds%js_south_uvw ! for winds on E/W sides of grid cells. + je_w=bnds%je_south_uvw !<-- + endif + endif +! + if(nside==3)then if(east_bc)then move=.true. bc_side_t0=>BC_t0%east @@ -5308,7 +5137,7 @@ subroutine regional_bc_t1_to_t0(BC_t1,BC_t0 & endif endif ! - if(nside==3)then + if(nside==4)then if(west_bc)then move=.true. bc_side_t0=>BC_t0%west @@ -5330,29 +5159,6 @@ subroutine regional_bc_t1_to_t0(BC_t1,BC_t0 & je_w=bnds%je_west_uvw !<-- endif endif -! - if(nside==4)then - if(south_bc)then - move=.true. - bc_side_t0=>BC_t0%south - bc_side_t1=>BC_t1%south -! - is_c=bnds%is_south !<-- - ie_c=bnds%ie_south ! South BC index limits - js_c=bnds%js_south ! for centers of grid cells. - je_c=bnds%je_south !<-- -! - is_s=bnds%is_south_uvs !<-- - ie_s=bnds%ie_south_uvs ! South BC index limits - js_s=bnds%js_south_uvs ! for winds on N/S sides of grid cells. - je_s=bnds%je_south_uvs !<-- -! - is_w=bnds%is_south_uvw !<-- - ie_w=bnds%ie_south_uvw ! South BC index limits - js_w=bnds%js_south_uvw ! for winds on E/W sides of grid cells. - je_w=bnds%je_south_uvw !<-- - endif - endif ! if(move)then do k=1,nlev