From ad0894b608edf3dbe28abde895310f1168d5cab8 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Sun, 27 Oct 2019 18:28:35 -0600 Subject: [PATCH 1/4] Loop from 1+nhalo rather than just nhalo, so don't go outside of array bounds --- libglissade/glissade_velo_higher.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/libglissade/glissade_velo_higher.F90 b/libglissade/glissade_velo_higher.F90 index 2d327f65..e83802fc 100644 --- a/libglissade/glissade_velo_higher.F90 +++ b/libglissade/glissade_velo_higher.F90 @@ -8848,8 +8848,8 @@ subroutine dirichlet_boundary_conditions_2d(nx, ny, & ! OK to skip vertices outside the global domain (i < nhalo or j < nhalo). ! Note: Need nhalo >= 2 so as not to step out of bounds. - do j = nhalo, ny-nhalo+1 - do i = nhalo, nx-nhalo+1 + do j = nhalo+1, ny-nhalo+1 + do i = nhalo+1, nx-nhalo+1 if (active_vertex(i,j)) then if (umask_dirichlet(i,j) == 1) then From ddc005605abbcc28185a8d6cccf5b6c73eee4afb Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Sun, 27 Oct 2019 18:41:47 -0600 Subject: [PATCH 2/4] Some changes that make sure array bounds are followed. pass unphys_val down to interpolate_thermal_forcing_to_lsrf so can check if thermal_forcing is set, if it isn't set it to zero. change loops that start at nhalo to nhalo+1, so subscripts remain within arrays. Set indices for nearest neighbors so can ensure within array bounds. Add more checks that both kbot and ktop are set to signify a valid ocean point. If ktop == kbot set both to zero. If checks on kbot or ktop, make sure it's not zero before doing the check. Comment out the bug check that all points are filled. It writes the error for several points. But, doing this allows the simulation to complete. --- libglissade/glissade_bmlt_float.F90 | 201 +++++++++++++++++----------- 1 file changed, 121 insertions(+), 80 deletions(-) diff --git a/libglissade/glissade_bmlt_float.F90 b/libglissade/glissade_bmlt_float.F90 index c58c3cf1..3f2d5824 100644 --- a/libglissade/glissade_bmlt_float.F90 +++ b/libglissade/glissade_bmlt_float.F90 @@ -1111,6 +1111,7 @@ subroutine glissade_bmlt_float_thermal_forcing(& bmlt_float_mask, & lake_mask, & lsrf, & + unphys_val, & ! identifies unfilled cells on input ocean_data%thermal_forcing, & ocean_data%thermal_forcing_lsrf) @@ -1324,6 +1325,8 @@ subroutine glissade_thermal_forcing_extrapolate(& integer :: i, j, k, iter integer :: iglobal, jglobal + integer :: iwe, iea, iso, ino ! longitude index in neighbor cells + integer :: jwe, jea, jso, jno ! longitude index in neighbor cells integer :: kw, ke, ks, kn ! ocean level in neighbor cells real(dp) :: phiw, phie, phin, phis ! field value in neighbor cells character(len=128) :: message @@ -1389,6 +1392,19 @@ subroutine glissade_thermal_forcing_extrapolate(& endif endif + if ( (kbot(i,j) >= 1) .and. (ktop(i,j) >= 1) )then + if ( kbot(i,j) == ktop(i,j) )then + !write(6,*) 'ktop==kbot, i, j, ktop, topg, lsrf = ', i, j, ktop(i,j), topg(i,j), lsrf(i,j) + !write(6,*) 'ocean_mask, bmlt_float_mask = ', ocean_mask(i,j), bmlt_float_mask(i,j) + !write(6,*) 'zocn = ', zocn + !call write_log('kbot == ktop', GM_FATAL) + ktop(i,j) = 0 + kbot(i,j) = 0 + end if + end if + if ( (kbot(i,j) > nzocn) .and. (ktop(i,j) > nzocn) )then + call write_log('kbot or ktop > nzocn', GM_FATAL) + end if enddo ! i enddo ! j @@ -1470,6 +1486,18 @@ subroutine glissade_thermal_forcing_extrapolate(& filled = .false. + ! + ! Indices of nearbor points, making sure to not go beyond array limits + ! + iwe = i-1 + if ( iwe < 1 ) iwe = 1 + iea = i+1 + if ( iea > nx ) iea = nx + jso = j-1 + if ( jso < 1 ) jso = 1 + jno = j+1 + if ( jno > ny ) jno = ny + ! Set thermal_forcing(k,i,j) to the mean value in filled neighbors at the same level if (k == ktop(i,j)) then @@ -1479,43 +1507,43 @@ subroutine glissade_thermal_forcing_extrapolate(& ! Note: Thermal forcing is corrected for the elevation difference, ! assuming linear dependence of Tf on zocn. - if (kbot(i-1,j) < k) then ! kbot in west neighbor lies above k in this cell - kw = kbot(i-1,j) - phiw = phi(kw,i-1,j) - dtocnfrz_dz * (zocn(kw) - zocn(k)) + if ( (kbot(iwe,j) > 0) .and. (kbot(iwe,j) < k) ) then ! kbot in west neighbor lies above k in this cell + kw = kbot(iwe,j) + phiw = phi(kw,iwe,j) - dtocnfrz_dz * (zocn(kw) - zocn(k)) else kw = k - phiw = phi(k,i-1,j) + phiw = phi(k,iwe,j) endif - if (kbot(i+1,j) < k) then ! kbot in east neighbor lies above k in this cell - ke = kbot(i+1,j) - phie = phi(ke,i+1,j) - dtocnfrz_dz * (zocn(ke) - zocn(k)) + if ( (kbot(iea,j) > 0) .and. (kbot(iea,j) < k) ) then ! kbot in east neighbor lies above k in this cell + ke = kbot(iea,j) + phie = phi(ke,iea,j) - dtocnfrz_dz * (zocn(ke) - zocn(k)) else ke = k - phie = phi(k,i+1,j) + phie = phi(k,iea,j) endif - if (kbot(i,j-1) < k) then ! kbot in south neighbor lies above k in this cell - ks = kbot(i,j-1) - phis = phi(ks,i,j-1) - dtocnfrz_dz * (zocn(ks) - zocn(k)) + if ( (kbot(i,jso) > 0) .and. (kbot(i,jso) < k) ) then ! kbot in south neighbor lies above k in this cell + ks = kbot(i,jso) + phis = phi(ks,i,jso) - dtocnfrz_dz * (zocn(ks) - zocn(k)) else ks = k - phis = phi(k,i,j-1) + phis = phi(k,i,jso) endif - if (kbot(i,j+1) < k) then ! kbot in north neighbor lies above k in this cell - kn = kbot(i,j+1) - phin = phi(kn,i,j+1) - dtocnfrz_dz * (zocn(kn) - zocn(k)) + if ( (kbot(i,jno) > 0) .and. (kbot(i,jno) < k) ) then ! kbot in north neighbor lies above k in this cell + kn = kbot(i,jno) + phin = phi(kn,i,jno) - dtocnfrz_dz * (zocn(kn) - zocn(k)) else kn = k - phin = phi(k,i,j+1) + phin = phi(k,i,jno) endif - sum_mask = mask(kw,i-1,j) + mask(ke,i+1,j) + mask(ks,i,j-1) + mask(kn,i,j+1) + sum_mask = mask(kw,iwe,j) + mask(ke,iea,j) + mask(ks,i,jso) + mask(kn,i,jno) if (sum_mask > 0) then - sum_phi = mask(kw,i-1,j)*phiw + mask(ke,i+1,j)*phie & - + mask(ks,i,j-1)*phis + mask(kn,i,j+1)*phin + sum_phi = mask(kw,iwe,j)*phiw + mask(ke,iea,j)*phie & + + mask(ks,i,jso)*phis + mask(kn,i,jno)*phin thermal_forcing(k,i,j) = sum_phi / real(sum_mask, dp) filled = .true. endif @@ -1531,43 +1559,43 @@ subroutine glissade_thermal_forcing_extrapolate(& ! need extra logic to allow spreading of values upward to shallower depth, ! in case ktop in a neighbor cell lies below kbot in this cell - if (ktop(i-1,j) > k) then ! ktop in west neighbor lies below k in this cell - kw = ktop(i-1,j) - phiw = phi(kw,i-1,j) - dtocnfrz_dz * (zocn(kw) - zocn(k)) + if (ktop(iwe,j) > k) then ! ktop in west neighbor lies below k in this cell + kw = ktop(iwe,j) + phiw = phi(kw,iwe,j) - dtocnfrz_dz * (zocn(kw) - zocn(k)) else kw = k - phiw = phi(k,i-1,j) + phiw = phi(k,iwe,j) endif - if (ktop(i+1,j) > k) then ! ktop in east neighbor lies below k in this cell - ke = ktop(i+1,j) - phie = phi(ke,i+1,j) - dtocnfrz_dz * (zocn(ke) - zocn(k)) + if (ktop(iea,j) > k) then ! ktop in east neighbor lies below k in this cell + ke = ktop(iea,j) + phie = phi(ke,iea,j) - dtocnfrz_dz * (zocn(ke) - zocn(k)) else ke = k - phie = phi(k,i+1,j) + phie = phi(k,iea,j) endif - if (ktop(i,j-1) > k) then ! ktop in south neighbor lies below k in this cell - ks = ktop(i,j-1) - phis = phi(ks,i,j-1) - dtocnfrz_dz * (zocn(ks) - zocn(k)) + if (ktop(i,jso) > k) then ! ktop in south neighbor lies below k in this cell + ks = ktop(i,jso) + phis = phi(ks,i,jso) - dtocnfrz_dz * (zocn(ks) - zocn(k)) else ks = k - phis = phi(k,i,j-1) + phis = phi(k,i,jso) endif - if (ktop(i,j+1) > k) then ! ktop in north neighbor lies below k in this cell - kn = ktop(i,j+1) - phin = phi(kn,i,j+1) - dtocnfrz_dz * (zocn(kn) - zocn(k)) + if (ktop(i,jno) > k) then ! ktop in north neighbor lies below k in this cell + kn = ktop(i,jno) + phin = phi(kn,i,jno) - dtocnfrz_dz * (zocn(kn) - zocn(k)) else kn = k - phin = phi(k,i,j+1) + phin = phi(k,i,jno) endif - sum_mask = mask(kw,i-1,j) + mask(ke,i+1,j) + mask(ks,i,j-1) + mask(kn,i,j+1) + sum_mask = mask(kw,iwe,j) + mask(ke,iea,j) + mask(ks,i,jso) + mask(kn,i,jno) if (sum_mask > 0) then - sum_phi = mask(kw,i-1,j)*phiw + mask(ke,i+1,j)*phie & - + mask(ks,i,j-1)*phis + mask(kn,i,j+1)*phin + sum_phi = mask(kw,iwe,j)*phiw + mask(ke,iea,j)*phie & + + mask(ks,i,jso)*phis + mask(kn,i,jno)*phin thermal_forcing(k,i,j) = sum_phi / real(sum_mask, dp) filled = .true. endif @@ -1578,11 +1606,11 @@ subroutine glissade_thermal_forcing_extrapolate(& ! simpler case; look only at neighbor levels with the same k value - sum_mask = mask(k,i-1,j) + mask(k,i+1,j) + mask(k,i,j-1) + mask(k,i,j+1) + sum_mask = mask(k,iwe,j) + mask(k,iea,j) + mask(k,i,jso) + mask(k,i,jno) if (sum_mask > 0) then - sum_phi = mask(k,i-1,j)*phi(k,i-1,j) + mask(k,i+1,j)*phi(k,i+1,j) & - + mask(k,i,j-1)*phi(k,i,j-1) + mask(k,i,j+1)*phi(k,i,j+1) + sum_phi = mask(k,iwe,j)*phi(k,iwe,j) + mask(k,iea,j)*phi(k,iea,j) & + + mask(k,i,jso)*phi(k,i,jso) + mask(k,i,jno)*phi(k,i,jno) thermal_forcing(k,i,j) = sum_phi / real(sum_mask, dp) filled = .true. endif @@ -1642,9 +1670,11 @@ subroutine glissade_thermal_forcing_extrapolate(& local_count = 0 do j = 1+nhalo, ny-nhalo do i = 1+nhalo, nx-nhalo - do k = ktop(i,j), kbot(i,j) - if (thermal_forcing(k,i,j) /= unphys_val) local_count = local_count + 1 - enddo + if (ktop(i,j) >= 1 .and. kbot(i,j) >= 1) then ! ocean or floating cell + do k = ktop(i,j), kbot(i,j) + if (thermal_forcing(k,i,j) /= unphys_val) local_count = local_count + 1 + enddo + endif enddo enddo @@ -1680,16 +1710,19 @@ subroutine glissade_thermal_forcing_extrapolate(& do j = 1, ny do i = 1, nx if (bmlt_float_mask(i,j) == 1 .and. lake_mask(i,j) == 0) then - do k = ktop(i,j), kbot(i,j) - if (thermal_forcing(k,i,j) == unphys_val) then - call parallel_globalindex(i, j, iglobal, jglobal) -!! print*, 'bmlt_float_mask, lake_mask =', bmlt_float_mask(i,j), lake_mask(i,j) -!! print*, 'ktop, kbot =', ktop(i,j), kbot(i,j) - write(message,*) & - 'Ocean data extrapolation error: did not fill level k, i, j =', k, iglobal, jglobal - call write_log(message, GM_FATAL) - endif - enddo ! k + if (ktop(i,j) >= 1 .and. kbot(i,j) >= 1) then ! ocean or floating cell + do k = ktop(i,j), kbot(i,j) + if (thermal_forcing(k,i,j) == unphys_val) then + call parallel_globalindex(i, j, iglobal, jglobal) + print*, 'bmlt_float_mask, lake_mask =', bmlt_float_mask(i,j), lake_mask(i,j) + print*, 'ktop, kbot, local-i, local-j, mask =', ktop(i,j), kbot(i,j), i, j, mask(k,i,j) + write(message,*) & + 'Ocean data extrapolation error: did not fill level k, i, j =', k, iglobal, jglobal + call write_log(message) + !call write_log(message, GM_FATAL) + endif + enddo ! k + endif endif ! floating and not lake enddo ! i enddo ! j @@ -1713,6 +1746,7 @@ subroutine interpolate_thermal_forcing_to_lsrf(& bmlt_float_mask, & lake_mask, & lsrf, & + unphys_val, & thermal_forcing, & thermal_forcing_lsrf) @@ -1734,6 +1768,9 @@ subroutine interpolate_thermal_forcing_to_lsrf(& real(dp), dimension(nx,ny), intent(in) :: & lsrf !> ice lower surface elevation (m), negative below sea level + real(dp), intent(in) :: & + unphys_val ! unphysical value given to cells/levels not yet filled + real(dp), dimension(nzocn,nx,ny), intent(in) :: & thermal_forcing !> thermal forcing field at ocean levels @@ -1758,11 +1795,15 @@ subroutine interpolate_thermal_forcing_to_lsrf(& if (bmlt_float_mask(i,j) == 1 .and. lake_mask(i,j) == 0) then if (lsrf(i,j) >= zocn(1)) then thermal_forcing_lsrf(i,j) = thermal_forcing(1,i,j) + if ( thermal_forcing_lsrf(i,j) == unphys_val ) thermal_forcing_lsrf(i,j) = 0.0d0 elseif (lsrf(i,j) < zocn(nzocn)) then thermal_forcing_lsrf(i,j) = thermal_forcing(nzocn,i,j) + if ( thermal_forcing_lsrf(i,j) == unphys_val ) thermal_forcing_lsrf(i,j) = 0.0d0 else + thermal_forcing_lsrf(i,j) = 0.0d0 do k = 1, nzocn-1 - if (lsrf(i,j) < zocn(k) .and. lsrf(i,j) >= zocn(k+1)) then + if ( (lsrf(i,j) < zocn(k) .and. lsrf(i,j) >= zocn(k+1)) .and. & + (thermal_forcing(k+1,i,j) /= unphys_val) .and. (thermal_forcing(k,i,j) /= unphys_val) )then dtf = thermal_forcing(k+1,i,j) - thermal_forcing(k,i,j) dzocn = zocn(k+1) - zocn(k) dzice = lsrf(i,j) - zocn(k) @@ -1784,11 +1825,11 @@ subroutine interpolate_thermal_forcing_to_lsrf(& !TODO - Remove this bug check if the ocean can realistically have TF < 0. do j = 1, ny do i = 1, nx - if (thermal_forcing_lsrf(i,j) < 0.0d0) then + if ( thermal_forcing_lsrf(i,j) < 0.0d0) then call parallel_globalindex(i, j, iglobal, jglobal) write(message,*) & - 'Ocean thermal forcing error: negative TF at level k, i, j, lsrf, TF =', & - k, iglobal, jglobal, lsrf(i,j), thermal_forcing_lsrf(i,j) + 'Ocean thermal forcing error: negative TF at i, j, lsrf, TF =', & + iglobal, jglobal, lsrf(i,j), thermal_forcing_lsrf(i,j) call write_log(message, GM_FATAL) endif enddo @@ -2497,8 +2538,8 @@ subroutine glissade_plume_melt_rate(& edge_mask_north(:,:) = 0 ! loop over all edges of locally owned cells - do j = nhalo, ny-nhalo - do i = nhalo, nx-nhalo + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo if (plume_mask_cell(i,j) == 1 .and. plume_mask_cell(i+1,j) == 1) then edge_mask_east(i,j) = 1 endif @@ -3709,8 +3750,8 @@ subroutine compute_edge_gradients(& ! Compute gradients at edges with plume cells on each side - do j = nhalo, ny-nhalo - do i = nhalo, nx-nhalo + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo ! east edges if (plume_mask_cell(i,j) == 1 .and. plume_mask_cell(i+1,j) == 1) then @@ -3727,8 +3768,8 @@ subroutine compute_edge_gradients(& ! Set gradients at edges that have a plume cell on one side and floating ice or water on the other. ! Extrapolate the gradient from the nearest neighbor edge. - do j = nhalo, ny-nhalo - do i = nhalo, nx-nhalo + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo ! east edges if (plume_mask_cell(i,j) == 1 .and. plume_mask_cell(i+1,j) == 0 .and. global_bndy_east(i,j) == 0) then @@ -3759,8 +3800,8 @@ subroutine compute_edge_gradients(& ! Average over 4 neighboring edges to estimate the y derivative on east edges and the x derivative on north edges. - do j = nhalo, ny-nhalo - do i = nhalo, nx-nhalo + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo ! y derivative on east edges df_dy_east(i,j) = 0.25d0 * (df_dy_north(i,j) + df_dy_north(i+1,j) & @@ -3930,8 +3971,8 @@ subroutine compute_plume_velocity(& !TODO - Use edge_mask_east instead? ! Maybe divu_mask_east is correct, since I stopped extrapolating, but should be called edge_mask_east. - do j = nhalo, ny-nhalo - do i = nhalo, nx-nhalo + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo ! PGF on east edge if (divu_mask_east(i,j) == 1) then @@ -4292,8 +4333,8 @@ subroutine compute_plume_velocity(& !TODO - Are global_bndy masks needed here? Wondering if we can avoid passing in 4 global_bndy fields. -! do j = nhalo, ny-nhalo -! do i = nhalo, nx-nhalo +! do j = nhalo+1, ny-nhalo +! do i = nhalo+1, nx-nhalo ! east edges ! if (plume_mask_cell(i,j) == 1 .and. plume_mask_cell(i+1,j) == 0 .and. global_bndy_east(i,j) == 0) then @@ -4325,8 +4366,8 @@ subroutine compute_plume_velocity(& ! enddo ! j ! Compute the plume speed at the edges (including the u_tidal term) - do j = nhalo, ny-nhalo - do i = nhalo, nx-nhalo + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo plume_speed_east(i,j) = sqrt(u_plume_east(i,j)**2 + v_plume_east(i,j)**2 + u_tidal**2) plume_speed_north(i,j) = sqrt(u_plume_north(i,j)**2 + v_plume_north(i,j)**2 + u_tidal**2) enddo ! i @@ -4493,8 +4534,8 @@ subroutine compute_velocity(& f_y(:,:) = pgf_y(:,:) + latdrag_y(:,:) ! Loop over edges of locally owned cells - do j = nhalo, ny-nhalo - do i = nhalo, nx-nhalo + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo if (edge_mask(i,j) == 1 .and. .not.converged_velo(i,j) ) then @@ -4668,8 +4709,8 @@ subroutine compute_lateral_drag(& ! Compute lateral drag terms ! Loop over edges of locally owned cells - do j = nhalo, ny-nhalo - do i = nhalo, nx-nhalo + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo ! lateral drag on east edges @@ -4933,8 +4974,8 @@ subroutine compute_entrainment(& ! entrainment_east(:,:) = 0.0d0 -! do j = nhalo, ny-nhalo -! do i = nhalo, nx-nhalo +! do j = nhalo+1, ny-nhalo +! do i = nhalo+1, nx-nhalo ! if (divu_mask_east(i,j) == 1) then @@ -4963,8 +5004,8 @@ subroutine compute_entrainment(& ! entrainment_north(:,:) = 0.0d0 -! do j = nhalo, ny-nhalo -! do i = nhalo, nx-nhalo +! do j = nhalo+1, ny-nhalo +! do i = nhalo+1, nx-nhalo ! if (divu_mask_north(i,j) == 1) then ! slope = sqrt(dlsrf_dx_north(i,j)**2 + dlsrf_dy_north(i,j)**2) From afe5603fdebd28c04e2440864f2611255d56ab53 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Fri, 7 Feb 2020 12:19:23 -0700 Subject: [PATCH 3/4] Back out my changes to the do i/do j loops, add back in a error check that is now passing, add in a verification that nhalo is large enough when basal melt is using thermal forcing --- libglissade/glissade.F90 | 8 +++- libglissade/glissade_bmlt_float.F90 | 59 ++++++++++++++-------------- libglissade/glissade_velo_higher.F90 | 4 +- 3 files changed, 39 insertions(+), 32 deletions(-) diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index b281dd14..43486103 100644 --- a/libglissade/glissade.F90 +++ b/libglissade/glissade.F90 @@ -89,7 +89,7 @@ subroutine glissade_initialise(model, evolve_ice) use parallel use glide_stop, only: register_model - use glide_setup + use glide_setup, only: glide_scale_params, glide_load_sigma use glimmer_ncio use glide_velo, only: init_velo !TODO - Remove call to init_velo? use glissade_therm, only: glissade_init_therm @@ -1072,6 +1072,7 @@ subroutine glissade_bmlt_float_solve(model) integer :: i, j integer :: ewn, nsn integer :: itest, jtest, rtest + character(len=100) :: message ! set grid dimensions ewn = model%general%ewn @@ -1161,6 +1162,11 @@ subroutine glissade_bmlt_float_solve(model) endif ! ISMIP6 nonlocal slope + if ( nhalo < 1 )then + write(message,*) 'nhalo is NOT large enough to use bmlt_thermal_forcing' + call write_log(message,GM_FATAL) + end if + call glissade_bmlt_float_thermal_forcing(& model%options%bmlt_float_thermal_forcing_param, & model%options%ocean_data_domain, & diff --git a/libglissade/glissade_bmlt_float.F90 b/libglissade/glissade_bmlt_float.F90 index 3f2d5824..492158b9 100644 --- a/libglissade/glissade_bmlt_float.F90 +++ b/libglissade/glissade_bmlt_float.F90 @@ -1394,9 +1394,10 @@ subroutine glissade_thermal_forcing_extrapolate(& endif if ( (kbot(i,j) >= 1) .and. (ktop(i,j) >= 1) )then if ( kbot(i,j) == ktop(i,j) )then - !write(6,*) 'ktop==kbot, i, j, ktop, topg, lsrf = ', i, j, ktop(i,j), topg(i,j), lsrf(i,j) - !write(6,*) 'ocean_mask, bmlt_float_mask = ', ocean_mask(i,j), bmlt_float_mask(i,j) - !write(6,*) 'zocn = ', zocn + write(6,*) 'ktop==kbot, i, j, ktop, topg, lsrf = ', i, j, ktop(i,j), topg(i,j), lsrf(i,j) + write(6,*) 'ocean_mask, bmlt_float_mask = ', ocean_mask(i,j), bmlt_float_mask(i,j) + write(6,*) 'zocn = ', zocn + call write_log('kbot == ktop') !call write_log('kbot == ktop', GM_FATAL) ktop(i,j) = 0 kbot(i,j) = 0 @@ -1490,13 +1491,13 @@ subroutine glissade_thermal_forcing_extrapolate(& ! Indices of nearbor points, making sure to not go beyond array limits ! iwe = i-1 - if ( iwe < 1 ) iwe = 1 + if ( iwe < 1 ) call write_log('Extrapolation out of bounds to the west', GM_FATAL) iea = i+1 - if ( iea > nx ) iea = nx + if ( iea > nx ) call write_log('Extrapolation out of bounds to the east', GM_FATAL) jso = j-1 - if ( jso < 1 ) jso = 1 + if ( jso < 1 ) call write_log('Extrapolation out of bounds to the south', GM_FATAL) jno = j+1 - if ( jno > ny ) jno = ny + if ( jno > ny ) call write_log('Extrapolation out of bounds to the north', GM_FATAL) ! Set thermal_forcing(k,i,j) to the mean value in filled neighbors at the same level @@ -2538,8 +2539,8 @@ subroutine glissade_plume_melt_rate(& edge_mask_north(:,:) = 0 ! loop over all edges of locally owned cells - do j = nhalo+1, ny-nhalo - do i = nhalo+1, nx-nhalo + do j = nhalo, ny-nhalo + do i = nhalo, nx-nhalo if (plume_mask_cell(i,j) == 1 .and. plume_mask_cell(i+1,j) == 1) then edge_mask_east(i,j) = 1 endif @@ -3750,8 +3751,8 @@ subroutine compute_edge_gradients(& ! Compute gradients at edges with plume cells on each side - do j = nhalo+1, ny-nhalo - do i = nhalo+1, nx-nhalo + do j = nhalo, ny-nhalo + do i = nhalo, nx-nhalo ! east edges if (plume_mask_cell(i,j) == 1 .and. plume_mask_cell(i+1,j) == 1) then @@ -3768,8 +3769,8 @@ subroutine compute_edge_gradients(& ! Set gradients at edges that have a plume cell on one side and floating ice or water on the other. ! Extrapolate the gradient from the nearest neighbor edge. - do j = nhalo+1, ny-nhalo - do i = nhalo+1, nx-nhalo + do j = nhalo, ny-nhalo + do i = nhalo, nx-nhalo ! east edges if (plume_mask_cell(i,j) == 1 .and. plume_mask_cell(i+1,j) == 0 .and. global_bndy_east(i,j) == 0) then @@ -3800,8 +3801,8 @@ subroutine compute_edge_gradients(& ! Average over 4 neighboring edges to estimate the y derivative on east edges and the x derivative on north edges. - do j = nhalo+1, ny-nhalo - do i = nhalo+1, nx-nhalo + do j = nhalo, ny-nhalo + do i = nhalo, nx-nhalo ! y derivative on east edges df_dy_east(i,j) = 0.25d0 * (df_dy_north(i,j) + df_dy_north(i+1,j) & @@ -3971,8 +3972,8 @@ subroutine compute_plume_velocity(& !TODO - Use edge_mask_east instead? ! Maybe divu_mask_east is correct, since I stopped extrapolating, but should be called edge_mask_east. - do j = nhalo+1, ny-nhalo - do i = nhalo+1, nx-nhalo + do j = nhalo, ny-nhalo + do i = nhalo, nx-nhalo ! PGF on east edge if (divu_mask_east(i,j) == 1) then @@ -4333,8 +4334,8 @@ subroutine compute_plume_velocity(& !TODO - Are global_bndy masks needed here? Wondering if we can avoid passing in 4 global_bndy fields. -! do j = nhalo+1, ny-nhalo -! do i = nhalo+1, nx-nhalo +! do j = nhalo, ny-nhalo +! do i = nhalo, nx-nhalo ! east edges ! if (plume_mask_cell(i,j) == 1 .and. plume_mask_cell(i+1,j) == 0 .and. global_bndy_east(i,j) == 0) then @@ -4366,8 +4367,8 @@ subroutine compute_plume_velocity(& ! enddo ! j ! Compute the plume speed at the edges (including the u_tidal term) - do j = nhalo+1, ny-nhalo - do i = nhalo+1, nx-nhalo + do j = nhalo, ny-nhalo + do i = nhalo, nx-nhalo plume_speed_east(i,j) = sqrt(u_plume_east(i,j)**2 + v_plume_east(i,j)**2 + u_tidal**2) plume_speed_north(i,j) = sqrt(u_plume_north(i,j)**2 + v_plume_north(i,j)**2 + u_tidal**2) enddo ! i @@ -4534,8 +4535,8 @@ subroutine compute_velocity(& f_y(:,:) = pgf_y(:,:) + latdrag_y(:,:) ! Loop over edges of locally owned cells - do j = nhalo+1, ny-nhalo - do i = nhalo+1, nx-nhalo + do j = nhalo, ny-nhalo + do i = nhalo, nx-nhalo if (edge_mask(i,j) == 1 .and. .not.converged_velo(i,j) ) then @@ -4709,8 +4710,8 @@ subroutine compute_lateral_drag(& ! Compute lateral drag terms ! Loop over edges of locally owned cells - do j = nhalo+1, ny-nhalo - do i = nhalo+1, nx-nhalo + do j = nhalo, ny-nhalo + do i = nhalo, nx-nhalo ! lateral drag on east edges @@ -4974,8 +4975,8 @@ subroutine compute_entrainment(& ! entrainment_east(:,:) = 0.0d0 -! do j = nhalo+1, ny-nhalo -! do i = nhalo+1, nx-nhalo +! do j = nhalo, ny-nhalo +! do i = nhalo, nx-nhalo ! if (divu_mask_east(i,j) == 1) then @@ -5004,8 +5005,8 @@ subroutine compute_entrainment(& ! entrainment_north(:,:) = 0.0d0 -! do j = nhalo+1, ny-nhalo -! do i = nhalo+1, nx-nhalo +! do j = nhalo, ny-nhalo +! do i = nhalo, nx-nhalo ! if (divu_mask_north(i,j) == 1) then ! slope = sqrt(dlsrf_dx_north(i,j)**2 + dlsrf_dy_north(i,j)**2) diff --git a/libglissade/glissade_velo_higher.F90 b/libglissade/glissade_velo_higher.F90 index e83802fc..2d327f65 100644 --- a/libglissade/glissade_velo_higher.F90 +++ b/libglissade/glissade_velo_higher.F90 @@ -8848,8 +8848,8 @@ subroutine dirichlet_boundary_conditions_2d(nx, ny, & ! OK to skip vertices outside the global domain (i < nhalo or j < nhalo). ! Note: Need nhalo >= 2 so as not to step out of bounds. - do j = nhalo+1, ny-nhalo+1 - do i = nhalo+1, nx-nhalo+1 + do j = nhalo, ny-nhalo+1 + do i = nhalo, nx-nhalo+1 if (active_vertex(i,j)) then if (umask_dirichlet(i,j) == 1) then From 2f3c26914b8938def0b3c9cad5bcb27d0b289c34 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Fri, 7 Feb 2020 13:49:01 -0700 Subject: [PATCH 4/4] Backout the change I had made to use indices for the neighbor cells --- libglissade/glissade_bmlt_float.F90 | 96 ++++++++++++++--------------- 1 file changed, 45 insertions(+), 51 deletions(-) diff --git a/libglissade/glissade_bmlt_float.F90 b/libglissade/glissade_bmlt_float.F90 index 492158b9..8381c89e 100644 --- a/libglissade/glissade_bmlt_float.F90 +++ b/libglissade/glissade_bmlt_float.F90 @@ -1325,8 +1325,6 @@ subroutine glissade_thermal_forcing_extrapolate(& integer :: i, j, k, iter integer :: iglobal, jglobal - integer :: iwe, iea, iso, ino ! longitude index in neighbor cells - integer :: jwe, jea, jso, jno ! longitude index in neighbor cells integer :: kw, ke, ks, kn ! ocean level in neighbor cells real(dp) :: phiw, phie, phin, phis ! field value in neighbor cells character(len=128) :: message @@ -1490,14 +1488,10 @@ subroutine glissade_thermal_forcing_extrapolate(& ! ! Indices of nearbor points, making sure to not go beyond array limits ! - iwe = i-1 - if ( iwe < 1 ) call write_log('Extrapolation out of bounds to the west', GM_FATAL) - iea = i+1 - if ( iea > nx ) call write_log('Extrapolation out of bounds to the east', GM_FATAL) - jso = j-1 - if ( jso < 1 ) call write_log('Extrapolation out of bounds to the south', GM_FATAL) - jno = j+1 - if ( jno > ny ) call write_log('Extrapolation out of bounds to the north', GM_FATAL) + if ( i-1 < 1 ) call write_log('Extrapolation out of bounds to the west', GM_FATAL) + if ( i+1 > nx ) call write_log('Extrapolation out of bounds to the east', GM_FATAL) + if ( j-1 < 1 ) call write_log('Extrapolation out of bounds to the south', GM_FATAL) + if ( j+1 > ny ) call write_log('Extrapolation out of bounds to the north', GM_FATAL) ! Set thermal_forcing(k,i,j) to the mean value in filled neighbors at the same level @@ -1508,43 +1502,43 @@ subroutine glissade_thermal_forcing_extrapolate(& ! Note: Thermal forcing is corrected for the elevation difference, ! assuming linear dependence of Tf on zocn. - if ( (kbot(iwe,j) > 0) .and. (kbot(iwe,j) < k) ) then ! kbot in west neighbor lies above k in this cell - kw = kbot(iwe,j) - phiw = phi(kw,iwe,j) - dtocnfrz_dz * (zocn(kw) - zocn(k)) + if ( (kbot(i-1,j) > 0) .and. (kbot(i-1,j) < k) ) then ! kbot in west neighbor lies above k in this cell + kw = kbot(i-1,j) + phiw = phi(kw,i-1,j) - dtocnfrz_dz * (zocn(kw) - zocn(k)) else kw = k - phiw = phi(k,iwe,j) + phiw = phi(k,i-1,j) endif - if ( (kbot(iea,j) > 0) .and. (kbot(iea,j) < k) ) then ! kbot in east neighbor lies above k in this cell - ke = kbot(iea,j) - phie = phi(ke,iea,j) - dtocnfrz_dz * (zocn(ke) - zocn(k)) + if ( (kbot(i+1,j) > 0) .and. (kbot(i+1,j) < k) ) then ! kbot in east neighbor lies above k in this cell + ke = kbot(i+1,j) + phie = phi(ke,i+1,j) - dtocnfrz_dz * (zocn(ke) - zocn(k)) else ke = k - phie = phi(k,iea,j) + phie = phi(k,i+1,j) endif - if ( (kbot(i,jso) > 0) .and. (kbot(i,jso) < k) ) then ! kbot in south neighbor lies above k in this cell - ks = kbot(i,jso) - phis = phi(ks,i,jso) - dtocnfrz_dz * (zocn(ks) - zocn(k)) + if ( (kbot(i,j-1) > 0) .and. (kbot(i,j-1) < k) ) then ! kbot in south neighbor lies above k in this cell + ks = kbot(i,j-1) + phis = phi(ks,i,j-1) - dtocnfrz_dz * (zocn(ks) - zocn(k)) else ks = k - phis = phi(k,i,jso) + phis = phi(k,i,j-1) endif - if ( (kbot(i,jno) > 0) .and. (kbot(i,jno) < k) ) then ! kbot in north neighbor lies above k in this cell - kn = kbot(i,jno) - phin = phi(kn,i,jno) - dtocnfrz_dz * (zocn(kn) - zocn(k)) + if ( (kbot(i,j+1) > 0) .and. (kbot(i,j+1) < k) ) then ! kbot in north neighbor lies above k in this cell + kn = kbot(i,j+1) + phin = phi(kn,i,j+1) - dtocnfrz_dz * (zocn(kn) - zocn(k)) else kn = k - phin = phi(k,i,jno) + phin = phi(k,i,j+1) endif - sum_mask = mask(kw,iwe,j) + mask(ke,iea,j) + mask(ks,i,jso) + mask(kn,i,jno) + sum_mask = mask(kw,i-1,j) + mask(ke,i+1,j) + mask(ks,i,j-1) + mask(kn,i,j+1) if (sum_mask > 0) then - sum_phi = mask(kw,iwe,j)*phiw + mask(ke,iea,j)*phie & - + mask(ks,i,jso)*phis + mask(kn,i,jno)*phin + sum_phi = mask(kw,i-1,j)*phiw + mask(ke,i+1,j)*phie & + + mask(ks,i,j-1)*phis + mask(kn,i,j+1)*phin thermal_forcing(k,i,j) = sum_phi / real(sum_mask, dp) filled = .true. endif @@ -1560,43 +1554,43 @@ subroutine glissade_thermal_forcing_extrapolate(& ! need extra logic to allow spreading of values upward to shallower depth, ! in case ktop in a neighbor cell lies below kbot in this cell - if (ktop(iwe,j) > k) then ! ktop in west neighbor lies below k in this cell - kw = ktop(iwe,j) - phiw = phi(kw,iwe,j) - dtocnfrz_dz * (zocn(kw) - zocn(k)) + if (ktop(i-1,j) > k) then ! ktop in west neighbor lies below k in this cell + kw = ktop(i-1,j) + phiw = phi(kw,i-1,j) - dtocnfrz_dz * (zocn(kw) - zocn(k)) else kw = k - phiw = phi(k,iwe,j) + phiw = phi(k,i-1,j) endif - if (ktop(iea,j) > k) then ! ktop in east neighbor lies below k in this cell - ke = ktop(iea,j) - phie = phi(ke,iea,j) - dtocnfrz_dz * (zocn(ke) - zocn(k)) + if (ktop(i+1,j) > k) then ! ktop in east neighbor lies below k in this cell + ke = ktop(i+1,j) + phie = phi(ke,i+1,j) - dtocnfrz_dz * (zocn(ke) - zocn(k)) else ke = k - phie = phi(k,iea,j) + phie = phi(k,i+1,j) endif - if (ktop(i,jso) > k) then ! ktop in south neighbor lies below k in this cell - ks = ktop(i,jso) - phis = phi(ks,i,jso) - dtocnfrz_dz * (zocn(ks) - zocn(k)) + if (ktop(i,j-1) > k) then ! ktop in south neighbor lies below k in this cell + ks = ktop(i,j-1) + phis = phi(ks,i,j-1) - dtocnfrz_dz * (zocn(ks) - zocn(k)) else ks = k - phis = phi(k,i,jso) + phis = phi(k,i,j-1) endif - if (ktop(i,jno) > k) then ! ktop in north neighbor lies below k in this cell - kn = ktop(i,jno) - phin = phi(kn,i,jno) - dtocnfrz_dz * (zocn(kn) - zocn(k)) + if (ktop(i,j+1) > k) then ! ktop in north neighbor lies below k in this cell + kn = ktop(i,j+1) + phin = phi(kn,i,j+1) - dtocnfrz_dz * (zocn(kn) - zocn(k)) else kn = k - phin = phi(k,i,jno) + phin = phi(k,i,j+1) endif - sum_mask = mask(kw,iwe,j) + mask(ke,iea,j) + mask(ks,i,jso) + mask(kn,i,jno) + sum_mask = mask(kw,i-1,j) + mask(ke,i+1,j) + mask(ks,i,j-1) + mask(kn,i,j+1) if (sum_mask > 0) then - sum_phi = mask(kw,iwe,j)*phiw + mask(ke,iea,j)*phie & - + mask(ks,i,jso)*phis + mask(kn,i,jno)*phin + sum_phi = mask(kw,i-1,j)*phiw + mask(ke,i+1,j)*phie & + + mask(ks,i,j-1)*phis + mask(kn,i,j+1)*phin thermal_forcing(k,i,j) = sum_phi / real(sum_mask, dp) filled = .true. endif @@ -1607,11 +1601,11 @@ subroutine glissade_thermal_forcing_extrapolate(& ! simpler case; look only at neighbor levels with the same k value - sum_mask = mask(k,iwe,j) + mask(k,iea,j) + mask(k,i,jso) + mask(k,i,jno) + sum_mask = mask(k,i-1,j) + mask(k,i+1,j) + mask(k,i,j-1) + mask(k,i,j+1) if (sum_mask > 0) then - sum_phi = mask(k,iwe,j)*phi(k,iwe,j) + mask(k,iea,j)*phi(k,iea,j) & - + mask(k,i,jso)*phi(k,i,jso) + mask(k,i,jno)*phi(k,i,jno) + sum_phi = mask(k,i-1,j)*phi(k,i-1,j) + mask(k,i+1,j)*phi(k,i+1,j) & + + mask(k,i,j-1)*phi(k,i,j-1) + mask(k,i,j+1)*phi(k,i,j+1) thermal_forcing(k,i,j) = sum_phi / real(sum_mask, dp) filled = .true. endif