Skip to content
New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Jump to bottom

Fixed an indexing bug in blending/remap_dwinds that was causing NaNs. #63

Merged
merged 3 commits into from
Feb 8, 2024
Merged

Fixed an indexing bug in blending/remap_dwinds that was causing NaNs. #63

Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0dcab6f
Fixed an indexing bug in blending/remap_dwinds that was causing NaNs.
Feb 6, 2024
1055ed9
additional index error
Feb 7, 2024
4581d09
cleaned up some debug statements.
Feb 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
170 changes: 129 additions & 41 deletions blending.fd/remap_dwinds/remap_dwinds.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,110 +1,198 @@
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
Copy link
Contributor

Choose a reason for hiding this comment

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

Do we need those checks?

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

!------
! 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
Copy link
Contributor

Choose a reason for hiding this comment

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

Do we need 1:3951 as comment here? this number is for NA 3km only.

Atm_v(i,j,k) = qn1(i,k)
enddo
enddo

endif


5000 continue

if(any(isnan(Atm_u))) then
Copy link
Contributor

Choose a reason for hiding this comment

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

Do we need those error check? Those are for debug, right?

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

!-------------------------------------------------------------------------------------------------
Expand Down
Loading