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

Some changes needed for the inversion_ismip6 branch for ice shelves #30

Open
wants to merge 5 commits into
base: lipscomb/inversion_ismip6
Choose a base branch
from
Open

Some changes needed for the inversion_ismip6 branch for ice shelves #30

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
ad0894b
Loop from 1+nhalo rather than just nhalo, so don't go outside of arra…
ekluzek Oct 28, 2019
ddc0056
Some changes that make sure array bounds are followed.
ekluzek Oct 28, 2019
dc11b1c
Merge branch 'lipscomb/inversion_ismip6' of https://github.com/ESCOMP…
ekluzek Oct 28, 2019
afe5603
Back out my changes to the do i/do j loops, add back in a error check…
ekluzek Feb 7, 2020
2f3c269
Backout the change I had made to use indices for the neighbor cells
ekluzek Feb 7, 2020
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
8 changes: 7 additions & 1 deletion libglissade/glissade.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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, &
Expand Down
78 changes: 57 additions & 21 deletions libglissade/glissade_bmlt_float.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -1389,6 +1390,20 @@ 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')
!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
Expand Down Expand Up @@ -1470,6 +1485,14 @@ subroutine glissade_thermal_forcing_extrapolate(&

filled = .false.

!
! Indices of nearbor points, making sure to not go beyond array limits
!
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

if (k == ktop(i,j)) then
Expand All @@ -1479,31 +1502,31 @@ 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
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,i-1,j)
endif

if (kbot(i+1,j) < k) then ! kbot in east neighbor lies above k in this cell
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,i+1,j)
endif

if (kbot(i,j-1) < k) then ! kbot in south neighbor lies above k in this cell
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,j-1)
endif

if (kbot(i,j+1) < k) then ! kbot in north neighbor lies above k in this cell
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
Expand Down Expand Up @@ -1642,9 +1665,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

Expand Down Expand Up @@ -1680,16 +1705,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
Expand All @@ -1713,6 +1741,7 @@ subroutine interpolate_thermal_forcing_to_lsrf(&
bmlt_float_mask, &
lake_mask, &
lsrf, &
unphys_val, &
thermal_forcing, &
thermal_forcing_lsrf)

Expand All @@ -1734,6 +1763,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

Expand All @@ -1758,11 +1790,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)
Expand All @@ -1784,11 +1820,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
Expand Down