From 0dcab6f51c7e7dad5a73c330160464989559d287 Mon Sep 17 00:00:00 2001 From: donald lippi Date: Tue, 6 Feb 2024 19:03:06 +0000 Subject: [PATCH 1/3] Fixed an indexing bug in blending/remap_dwinds that was causing NaNs. --- blending.fd/remap_dwinds/remap_dwinds.f90 | 170 ++++++++++++++++------ 1 file changed, 129 insertions(+), 41 deletions(-) diff --git a/blending.fd/remap_dwinds/remap_dwinds.f90 b/blending.fd/remap_dwinds/remap_dwinds.f90 index c0a06ec..41be69c 100644 --- a/blending.fd/remap_dwinds/remap_dwinds.f90 +++ b/blending.fd/remap_dwinds/remap_dwinds.f90 @@ -1,70 +1,135 @@ subroutine main(km, npz, ak0, bk0, Atm_ak, Atm_bk, psc, ud, vd, is, ie, js, je, Atm_u, Atm_v, Atm_ps) use ISO_FORTRAN_ENV use omp_lib - !use, intrinsic :: ieee_arithmetic + use, intrinsic :: ieee_arithmetic implicit none integer, parameter :: r8_kind = selected_real_kind(15) ! 15 decimal digits - integer, intent(IN) :: is, ie, js, je - integer, intent(IN) :: km ! 128 - integer, intent(IN) :: npz ! 127 - real(kind=8), dimension(:), intent(IN) :: ak0 ! (129,) - real(kind=8), dimension(:), intent(IN) :: bk0 ! (129,) - real(kind=8), dimension(:,:), intent(IN) :: psc ! (768, 768) - real(kind=8), dimension(:,:), intent(IN) :: Atm_ps ! - real(kind=8), dimension(:,:,:), intent(IN) :: ud - real(kind=8), dimension(:,:,:), intent(IN) :: vd - real(kind=8), dimension(:), intent(IN) :: Atm_ak ! (128,) - real(kind=8), dimension(:), intent(IN) :: Atm_bk ! (128,) - real(kind=8), dimension(:,:,:), intent(INOUT) :: Atm_u ! - real(kind=8), dimension(:,:,:), intent(INOUT) :: Atm_v ! - real(kind=8) :: Atm_ptop - real(kind=8), dimension(is:ie,js:je) :: psd - real(kind=8), dimension(is:ie+1,1:npz) :: qn1 - real(kind=8), dimension(is:ie+1,1:km+1) :: pe0 - real(kind=8), dimension(is:ie+1,1:npz+1) :: pe1 + integer, intent(IN) :: is ! 1 + integer, intent(IN) :: ie ! 3950 + integer, intent(IN) :: js ! 1 + integer, intent(IN) :: je ! 2700 + integer, intent(IN) :: km ! 66 + integer, intent(IN) :: npz ! 65 + real(kind=8), dimension(:), intent(IN) :: ak0 !(67) + real(kind=8), dimension(:), intent(IN) :: bk0 !(67) + real(kind=8), dimension(:,:), intent(IN) :: psc !(3950, 2700) + real(kind=8), dimension(:,:), intent(IN) :: Atm_ps !(3950, 2700) + real(kind=8), dimension(:,:,:), intent(IN) :: ud !(3950, 2701, 66) + real(kind=8), dimension(:,:,:), intent(IN) :: vd !(3951, 2700, 66) + real(kind=8), dimension(:), intent(IN) :: Atm_ak !(66) + real(kind=8), dimension(:), intent(IN) :: Atm_bk !(66) + real(kind=8), dimension(:,:,:), intent(INOUT) :: Atm_u !(3950, 2701, 65) + real(kind=8), dimension(:,:,:), intent(INOUT) :: Atm_v !(3951, 2700, 65) + real(kind=8) :: Atm_ptop ! + real(kind=8), dimension(is:ie,js:je) :: psd !(3950, 2700) + real(kind=8), dimension(is:ie+1,1:npz) :: qn1 !(3951, 65) + real(kind=8), dimension(is:ie+1,1:km+1) :: pe0 !(3951, 67) + real(kind=8), dimension(is:ie+1,1:npz+1) :: pe1 !(3951, 66) integer :: i,j,k,itoa logical :: no_boundary + real (kind=8) :: NaN + NaN=IEEE_VALUE(NaN, IEEE_SIGNALING_NAN) + + write(6,'("Info: ",6I6)'),is,ie,js,je,npz,km + Atm_ptop = NaN + psd = NaN + qn1 = NaN + pe0 = NaN + pe1 = NaN + itoa = km - npz + 1 Atm_ptop = Atm_ak(1) + write(6,'("Info: psc ",2I6)'),size(psc,1), size(psc,2) psd = psc + + if(any(isnan(psd))) then + write(6,'("Warning: Found some NaNs in psc array before 5000 loop")') + do j=js,je+1 + do i=is,ie + if(isnan(psd(i,j))) write(6,'("Error: Found NaN in psd at ",2I6)'),i,j + enddo + enddo + endif + !psd = Atm_ps + if(any(isnan(ak0))) then + write(6,'("Warning: Found some NaNs in ak0 array before 5000 loop")') + endif + if(any(isnan(bk0))) then + write(6,'("Warning: Found some NaNs in bk0 array before 5000 loop")') + endif + if(any(isnan(psc))) then + write(6,'("Warning: Found some NaNs in psc array before 5000 loop")') + endif + if(any(isnan(Atm_ak))) then + write(6,'("Warning: Found some NaNs in Atm_ak array before 5000 loop")') + endif + if(any(isnan(Atm_bk))) then + write(6,'("Warning: Found some NaNs in Atm_bk array before 5000 loop")') + endif + + if(any(isnan(Atm_ps))) then + write(6,'("Warning: Found some NaNs in Atm_ps array before 5000 loop")') + endif + + if(any(isnan(ud))) then + write(6,'("Warning: Found some NaNs in ud array before 5000 loop")') + endif + if(any(isnan(vd))) then + write(6,'("Warning: Found some NaNs in vd array before 5000 loop")') + endif + + if(any(isnan(Atm_u))) then + write(6,'("Warning: Found some NaNs in Atm_u array before 5000 loop")') + endif + + if(any(isnan(Atm_v))) then + write(6,'("Warning: Found some NaNs in Atm_v array before 5000 loop")') + endif + + write(6,'("Info: ak ",2I6)'),size(ak0), size(Atm_ak) ! 67 66 + write(6,'("Info: bk ",2I6)'),size(bk0), size(Atm_bk) ! 67 66 + + !$OMP parallel do default(none) & !$OMP shared(is,ie,js,je,npz,km,ak0,bk0,psc,psd,ud,vd,Atm_ak, & !$OMP Atm_bk,Atm_u,Atm_v,Atm_ps,Atm_ptop) & !$OMP private(pe1,pe0,qn1) - do 5000 j=js,je+1 + do 5000 j=js,je+1 ! 1:2701 !------ ! map u !------ !pressure at layer edges (from model top to bottom surface) in the original vertical coordinate - do k=1,km+1 - do i=is,ie - if(j==1) then - pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i,j)+psd(i,j)) + do k=1,km+1 ! 1:67 + do i=is,ie ! 1:3950 + if(j==js) then + pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i,j )+psd(i,j )) + elseif(j==je+1) then + pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i,j-1)+psd(i,j-1)) else - pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i,j-1)+psd(i,j)) + pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i,j-1)+psd(i,j )) endif enddo enddo !pressure at layer edges (from model top to bottom surface) in the new vertical coordinate - do k=1,npz+1 - do i=is,ie + do k=1,npz+1 ! 1:66 + do i=is,ie ! 1:3950 if(j==1) then - pe1(i,k) = Atm_ak(k) + Atm_bk(k)*0.5*(Atm_ps(i,j)+Atm_ps(i,j)) + pe1(i,k) = Atm_ak(k) + Atm_bk(k)*0.5*(Atm_ps(i,j )+Atm_ps(i,j)) else - pe1(i,k) = Atm_ak(k) + Atm_bk(k)*0.5*(Atm_ps(i,j-1)+Atm_ps(i,j)) + pe1(i,k) = Atm_ak(k) + Atm_bk(k)*0.5*(Atm_ps(i,j-1)+Atm_ps(i,j)) endif enddo enddo call mappm(km, pe0(is:ie,1:km+1), ud(is:ie,j,1:km), npz, pe1(is:ie,1:npz+1), & qn1(is:ie,1:npz), is,ie, -1, 8, Atm_ptop) - do k=1,npz - do i=is,ie + do k=1,npz ! 1:65 + do i=is,ie ! 1:3950 Atm_u(i,j,k) = qn1(i,k) enddo enddo @@ -72,39 +137,62 @@ subroutine main(km, npz, ak0, bk0, Atm_ak, Atm_bk, psc, ud, vd, is, ie, js, je, !------ ! map v !------ - if ( j/=(je+1) ) then + if ( j/=(je+1) ) then ! 1:2700 do k=1,km+1 do i=is,ie+1 if(i==1) then - pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i,j)+psd(i,j)) + pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i ,j)+psd(i ,j)) + elseif(i==ie+1) then + pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i-1,j)+psd(i-1,j)) else - pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i-1,j)+psd(i,j)) + pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i-1,j)+psd(i ,j)) endif enddo enddo - do k=1,npz+1 - do i=is,ie+1 + do k=1,npz+1 ! 1:66 + do i=is,ie+1 ! 1:3951 if(i==1) then - pe1(i,k) = Atm_ak(k) + Atm_bk(k)*0.5*(Atm_ps(i,j)+Atm_ps(i,j)) + pe1(i,k) = Atm_ak(k) + Atm_bk(k)*0.5*(Atm_ps(i,j)+Atm_ps(i,j)) else - pe1(i,k) = Atm_ak(k) + Atm_bk(k)*0.5*(Atm_ps(i-1,j)+Atm_ps(i,j)) + pe1(i,k) = Atm_ak(k) + Atm_bk(k)*0.5*(Atm_ps(i-1,j)+Atm_ps(i,j)) endif enddo enddo call mappm(km, pe0(is:ie+1,1:km+1), vd(is:ie+1,j,1:km), npz, pe1(is:ie+1,1:npz+1), & qn1(is:ie+1,1:npz), is,ie+1, -1, 8, Atm_ptop) - do k=1,npz - do i=is,ie+1 + do k=1,npz ! 1:65 + do i=is,ie+1 ! 1:3951 Atm_v(i,j,k) = qn1(i,k) enddo enddo endif - 5000 continue + if(any(isnan(Atm_u))) then + write(6,'("Warning: Found some NaNs in Atm_u array after 5000 loop",3I6)'),size(Atm_u,1),size(Atm_u,2),size(Atm_u,3) + do k=1,npz !size(Atm_u,3) + do j=js,je+1 + do i=is,ie + if(isnan(Atm_u(i,j,k))) write(6,'("Error: Found NaN in Atm_u at ",3I6)'),i,j,k + enddo + enddo + enddo + endif + + if(any(isnan(Atm_v))) then + write(6,'("Warning: Found some NaNs in Atm_v array after 5000 loop",3I6)'),size(Atm_v,1),size(Atm_v,2),size(Atm_v,3) + do k=1,npz !size(Atm_v,3) + do j=js,je + do i=is,ie+1 + if(isnan(Atm_u(i,j,k))) write(6,'("Error: Found NaN in Atm_v at ",3I6)'),i,j,k + enddo + enddo + enddo + endif + end subroutine main !------------------------------------------------------------------------------------------------- From 1055ed99b49e3942cd30d1178326ed02382b0e23 Mon Sep 17 00:00:00 2001 From: donald lippi Date: Wed, 7 Feb 2024 02:14:29 +0000 Subject: [PATCH 2/3] additional index error --- blending.fd/remap_dwinds/remap_dwinds.f90 | 24 +++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/blending.fd/remap_dwinds/remap_dwinds.f90 b/blending.fd/remap_dwinds/remap_dwinds.f90 index 41be69c..dd31552 100644 --- a/blending.fd/remap_dwinds/remap_dwinds.f90 +++ b/blending.fd/remap_dwinds/remap_dwinds.f90 @@ -109,10 +109,10 @@ subroutine main(km, npz, ak0, bk0, Atm_ak, Atm_bk, psc, ud, vd, is, ie, js, je, do i=is,ie ! 1:3950 if(j==js) then pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i,j )+psd(i,j )) - elseif(j==je+1) then - pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i,j-1)+psd(i,j-1)) - else + elseif(j Date: Wed, 7 Feb 2024 02:57:48 +0000 Subject: [PATCH 3/3] cleaned up some debug statements. --- blending.fd/remap_dwinds/remap_dwinds.f90 | 108 +++------------------- 1 file changed, 13 insertions(+), 95 deletions(-) diff --git a/blending.fd/remap_dwinds/remap_dwinds.f90 b/blending.fd/remap_dwinds/remap_dwinds.f90 index dd31552..8012d18 100644 --- a/blending.fd/remap_dwinds/remap_dwinds.f90 +++ b/blending.fd/remap_dwinds/remap_dwinds.f90 @@ -1,7 +1,7 @@ subroutine main(km, npz, ak0, bk0, Atm_ak, Atm_bk, psc, ud, vd, is, ie, js, je, Atm_u, Atm_v, Atm_ps) use ISO_FORTRAN_ENV use omp_lib - use, intrinsic :: ieee_arithmetic + !use, intrinsic :: ieee_arithmetic implicit none integer, parameter :: r8_kind = selected_real_kind(15) ! 15 decimal digits integer, intent(IN) :: is ! 1 @@ -29,84 +29,24 @@ subroutine main(km, npz, ak0, bk0, Atm_ak, Atm_bk, psc, ud, vd, is, ie, js, je, integer :: i,j,k,itoa logical :: no_boundary - real (kind=8) :: NaN - NaN=IEEE_VALUE(NaN, IEEE_SIGNALING_NAN) - - write(6,'("Info: ",6I6)'),is,ie,js,je,npz,km - Atm_ptop = NaN - psd = NaN - qn1 = NaN - pe0 = NaN - pe1 = NaN - itoa = km - npz + 1 Atm_ptop = Atm_ak(1) - write(6,'("Info: psc ",2I6)'),size(psc,1), size(psc,2) psd = psc - - if(any(isnan(psd))) then - write(6,'("Warning: Found some NaNs in psc array before 5000 loop")') - do j=js,je+1 - do i=is,ie - if(isnan(psd(i,j))) write(6,'("Error: Found NaN in psd at ",2I6)'),i,j - enddo - enddo - endif - !psd = Atm_ps - if(any(isnan(ak0))) then - write(6,'("Warning: Found some NaNs in ak0 array before 5000 loop")') - endif - if(any(isnan(bk0))) then - write(6,'("Warning: Found some NaNs in bk0 array before 5000 loop")') - endif - if(any(isnan(psc))) then - write(6,'("Warning: Found some NaNs in psc array before 5000 loop")') - endif - if(any(isnan(Atm_ak))) then - write(6,'("Warning: Found some NaNs in Atm_ak array before 5000 loop")') - endif - if(any(isnan(Atm_bk))) then - write(6,'("Warning: Found some NaNs in Atm_bk array before 5000 loop")') - endif - - if(any(isnan(Atm_ps))) then - write(6,'("Warning: Found some NaNs in Atm_ps array before 5000 loop")') - endif - - if(any(isnan(ud))) then - write(6,'("Warning: Found some NaNs in ud array before 5000 loop")') - endif - if(any(isnan(vd))) then - write(6,'("Warning: Found some NaNs in vd array before 5000 loop")') - endif - - if(any(isnan(Atm_u))) then - write(6,'("Warning: Found some NaNs in Atm_u array before 5000 loop")') - endif - - if(any(isnan(Atm_v))) then - write(6,'("Warning: Found some NaNs in Atm_v array before 5000 loop")') - endif - - write(6,'("Info: ak ",2I6)'),size(ak0), size(Atm_ak) ! 67 66 - write(6,'("Info: bk ",2I6)'),size(bk0), size(Atm_bk) ! 67 66 - - !$OMP parallel do default(none) & !$OMP shared(is,ie,js,je,npz,km,ak0,bk0,psc,psd,ud,vd,Atm_ak, & !$OMP Atm_bk,Atm_u,Atm_v,Atm_ps,Atm_ptop) & !$OMP private(pe1,pe0,qn1) - do 5000 j=js,je+1 ! 1:2701 + do 5000 j=js,je+1 !------ ! map u !------ !pressure at layer edges (from model top to bottom surface) in the original vertical coordinate - do k=1,km+1 ! 1:67 - do i=is,ie ! 1:3950 + do k=1,km+1 + do i=is,ie if(j==js) then pe0(i,k) = ak0(k) + bk0(k)*0.5*(psd(i,j )+psd(i,j )) elseif(j