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

Add adaptive timestepping for lateral melt and growth #298

Merged
merged 8 commits into from
Feb 14, 2020
Merged

Add adaptive timestepping for lateral melt and growth #298

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
e6af1a2
Checking latm timestep
lettie-roach Jan 29, 2020
1c86722
Debugging
lettie-roach Jan 31, 2020
5800f93
Added in adaptive timestepping for lateral melt and growth
lettie-roach Feb 6, 2020
648a3d2
Clean up
lettie-roach Feb 10, 2020
bc01a07
Checked diffs to master
lettie-roach Feb 10, 2020
a081aad
Move get_subdt to icepack_fsd
lettie-roach Feb 10, 2020
0030e1f
Add note to docs on FSD adaptive timestepping
lettie-roach Feb 12, 2020
7f5630f
Add reference
lettie-roach Feb 12, 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
108 changes: 88 additions & 20 deletions columnphysics/icepack_fsd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ module icepack_fsd

private
public :: icepack_init_fsd_bounds, icepack_init_fsd, icepack_cleanup_fsd, &
fsd_lateral_growth, fsd_add_new_ice, fsd_weld_thermo
fsd_lateral_growth, fsd_add_new_ice, fsd_weld_thermo, get_subdt_fsd

real(kind=dbl_kind), dimension(:), allocatable :: &
floe_rad_h, & ! fsd size higher bound in m (radius)
Expand Down Expand Up @@ -565,6 +565,13 @@ subroutine fsd_lateral_growth (ncat, nfsd, &
+ c2*aicen(n)*afsdn(k,n)*G_radial*dt/floe_rad_c(k)
end do
end do ! n

! cannot expand ice laterally beyond lead region
if (SUM(d_an_latg(:)).ge.lead_area) then
d_an_latg(:) = d_an_latg(:)/SUM(d_an_latg(:))
d_an_latg(:) = d_an_latg(:)*lead_area
end if

endif ! vi0new_lat > 0

! Use remaining ice volume as in standard model,
Expand Down Expand Up @@ -649,44 +656,67 @@ subroutine fsd_add_new_ice (ncat, n, nfsd, &
d_afsd_newi ! new ice formation

integer (kind=int_kind) :: &
k ! floe size category index
k, & ! floe size category index
new_size, & ! index for floe size of new ice
nsubt ! number of subtimesteps

real (kind=dbl_kind) :: &
elapsed_t, subdt ! elapsed time, subtimestep (s)

real (kind=dbl_kind), dimension (nfsd,ncat) :: &
afsdn_latg ! fsd after lateral growth

real (kind=dbl_kind), dimension (nfsd) :: &
df_flx , & ! finite differences for G_r*tilda(L)
dafsd_tmp, & ! tmp FSD
df_flx , & ! finite differences for fsd
afsd_ni ! fsd after new ice added

real (kind=dbl_kind), dimension(nfsd+1) :: &
f_flx ! finite differences in floe size

integer (kind=int_kind) :: &
new_size ! index for floe size of new ice

character(len=*),parameter :: subname='(fsd_add_new_ice)'

afsdn_latg(:,n) = afsdn(:,n) ! default

if (d_an_latg(n) > puny) then ! lateral growth

df_flx(:) = c0 ! NB could stay zero if all in largest FS cat
f_flx (:) = c0
do k = 2, nfsd
f_flx(k) = G_radial * afsdn(k-1,n) / floe_binwidth(k-1)
end do
do k = 1, nfsd
df_flx(k) = f_flx(k+1) - f_flx(k)
end do
! adaptive timestep
elapsed_t = c0
nsubt = 0

DO WHILE (elapsed_t.lt.dt)

nsubt = nsubt + 1
if (nsubt.gt.100) print *, 'latg not converging'

! finite differences
df_flx(:) = c0 ! NB could stay zero if all in largest FS cat
f_flx (:) = c0
do k = 2, nfsd
f_flx(k) = G_radial * afsdn_latg(k-1,n) / floe_binwidth(k-1)
end do
do k = 1, nfsd
df_flx(k) = f_flx(k+1) - f_flx(k)
end do

! if (abs(sum(df_flx)) > puny) print*,'fsd_add_new ERROR df_flx /= 0'

afsdn_latg(:,n) = c0
do k = 1, nfsd
afsdn_latg(k,n) = afsdn(k,n) &
+ dt * (-df_flx(k) + c2 * G_radial * afsdn(k,n) &
* (c1/floe_rad_c(k) - SUM(afsdn(:,n)/floe_rad_c(:))) )
end do
dafsd_tmp(:) = c0
do k = 1, nfsd
dafsd_tmp(k) = (-df_flx(k) + c2 * G_radial * afsdn_latg(k,n) &
* (c1/floe_rad_c(k) - SUM(afsdn_latg(:,n)/floe_rad_c(:))) )

end do

! timestep required for this
subdt = get_subdt_fsd(nfsd, afsdn_latg(:,n), dafsd_tmp(:))
subdt = MIN(subdt, dt)

! update fsd and elapsed time
afsdn_latg(:,n) = afsdn_latg(:,n) + subdt*dafsd_tmp(:)
elapsed_t = elapsed_t + subdt

END DO

call icepack_cleanup_fsdn (nfsd, afsdn_latg(:,n))
trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsdn_latg(:,n)
Expand Down Expand Up @@ -981,6 +1011,44 @@ subroutine fsd_weld_thermo (ncat, nfsd, &

end subroutine fsd_weld_thermo

!=======================================================================
!
! Adaptive timestepping (process agnostic)
! See reference: Horvat & Tziperman (2017) JGR, Appendix A
!
! authors: 2018 Lettie Roach, NIWA/VUW
!
!
function get_subdt_fsd(nfsd, afsd_init, d_afsd) &
result(subdt)

integer (kind=int_kind), intent(in) :: &
nfsd ! number of floe size categories

real (kind=dbl_kind), dimension (nfsd), intent(in) :: &
afsd_init, d_afsd ! floe size distribution tracer

! output
real (kind=dbl_kind) :: &
subdt ! subcycle timestep (s)

! local variables
real (kind=dbl_kind), dimension (nfsd) :: &
check_dt ! to compute subcycle timestep (s)

integer (kind=int_kind) :: k

check_dt(:) = bignum
do k = 1, nfsd
if (d_afsd(k) > puny) check_dt(k) = (1-afsd_init(k))/d_afsd(k)
if (d_afsd(k) < -puny) check_dt(k) = afsd_init(k)/ABS(d_afsd(k))
end do

subdt = MINVAL(check_dt)

end function get_subdt_fsd


!=======================================================================

end module icepack_fsd
Expand Down
110 changes: 59 additions & 51 deletions columnphysics/icepack_therm_itd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ module icepack_therm_itd
use icepack_warnings, only: warnstr, icepack_warnings_add
use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted

use icepack_fsd, only: fsd_weld_thermo, icepack_cleanup_fsd
use icepack_fsd, only: fsd_weld_thermo, icepack_cleanup_fsd, get_subdt_fsd
use icepack_itd, only: reduce_area, cleanup_itd
use icepack_itd, only: aggregate_area, shift_ice
use icepack_itd, only: column_sum, column_conservation_check
Expand All @@ -50,7 +50,7 @@ module icepack_therm_itd
use icepack_therm_shared, only: hi_min
use icepack_zbgc, only: add_new_ice_bgc
use icepack_zbgc, only: lateral_melt_bgc

implicit none

private
Expand Down Expand Up @@ -915,7 +915,8 @@ subroutine lateral_melt (dt, ncat, &

integer (kind=int_kind) :: &
n , & ! thickness category index
k ! layer index
k , & ! layer index
nsubt ! sub timesteps for FSD tendency

real (kind=dbl_kind) :: &
dfhocn , & ! change in fhocn
Expand All @@ -942,7 +943,8 @@ subroutine lateral_melt (dt, ncat, &
afsdn_init ! initial value

real (kind=dbl_kind), dimension (nfsd) :: &
df_flx ! finite difference for G_r * areal mFSTD tilda
df_flx, & ! finite difference for FSD
afsd_tmp, d_afsd_tmp

real (kind=dbl_kind), dimension(nfsd+1) :: &
f_flx !
Expand All @@ -951,7 +953,9 @@ subroutine lateral_melt (dt, ncat, &
real (kind=dbl_kind), intent(in) :: &
sss
real (kind=dbl_kind) :: &
Ti, Si0, qi0
Ti, Si0, qi0, &
elapsed_t, & ! FSD subcycling
subdt ! FSD timestep (s)

character(len=*), parameter :: subname='(lateral_melt)'

Expand Down Expand Up @@ -1009,8 +1013,6 @@ subroutine lateral_melt (dt, ncat, &

if (qin(n) < -puny) G_radialn(n) = -fside/qin(n) ! negative

! if (G_radialn(n) > puny) stop 'G_radial positive'

if (G_radialn(n) < -puny) then


Expand All @@ -1032,7 +1034,7 @@ subroutine lateral_melt (dt, ncat, &
if (delta_an(n) > c0) print*,'ERROR delta_an > 0', delta_an(n)

! following original code, not necessary for fsd
if (aicen(n) > c0) rsiden(n) = -delta_an(n)/aicen(n)
if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n)/aicen(n),c1)

if (rsiden(n) < c0) print*,'ERROR rsiden < 0', rsiden(n)

Expand All @@ -1048,18 +1050,8 @@ subroutine lateral_melt (dt, ncat, &

if (flag) then ! grid cells with lateral melting.

! LR is this necessary?
!tmp = SUM(rsiden(:))
do n = 1, ncat

!if (tr_fsd) then
! if (tmp > c0) then
! rsiden(n) = rsiden(n)/tmp
! else
! rsiden(n) = c0
! end if
!end if

!-----------------------------------------------------------------
! Melt the ice and increment fluxes.
!-----------------------------------------------------------------
Expand Down Expand Up @@ -1090,40 +1082,56 @@ subroutine lateral_melt (dt, ncat, &
if (tr_fsd) then
if (rsiden(n) > puny) then
if (aicen(n) > puny) then
df_flx(:) = c0
f_flx (:) = c0
do k = 2, nfsd
f_flx(k) = G_radialn(n) * afsdn_init(k,n) / floe_binwidth(k)
end do

do k = 1, nfsd
df_flx(k) = f_flx(k+1) - f_flx(k)
end do

if (abs(sum(df_flx(:))) > puny) &
print*,'sum(df_flx)/=0'

tmp = SUM(afsdn_init(:,n)/floe_rad_c(:))
do k = 1, nfsd
afsdn (k,n) = afsdn_init(k,n) &
+ dt * (-df_flx(k) + c2 * G_radialn(n) * afsdn_init(k,n) &
* (c1/floe_rad_c(k) - tmp))
end do

if (abs(sum(afsdn(:,n))-c1) > puny) &
print*,'lateral_melt E afsdn not normed',sum(df_flx), sum(afsdn(:,n))-c1
if (any(afsdn < -puny)) &
print*,'lateral_melt: afsdn < 0'
if (any(afsdn > c1+puny)) &
print*,'lateral_melt: afsdn > 1'


!trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsdn(:,n)

!else ! aicen = 0

! trcrn(nt_fsd:nt_fsd+nfsd-1,n) = c0
! adaptive subtimestep
elapsed_t = c0
afsd_tmp(:) = afsdn_init(:,n)
d_afsd_tmp(:) = c0
nsubt = 0

DO WHILE (elapsed_t.lt.dt)

nsubt = nsubt + 1
if (nsubt.gt.100) &
print *, 'latm not converging'

! finite differences
df_flx(:) = c0
f_flx (:) = c0
do k = 2, nfsd
f_flx(k) = G_radialn(n) * afsd_tmp(k) / floe_binwidth(k)
end do

do k = 1, nfsd
df_flx(k) = f_flx(k+1) - f_flx(k)
end do

if (abs(sum(df_flx(:))) > puny) &
print*,'sum(df_flx)/=0'

! this term ensures area conservation
tmp = SUM(afsd_tmp(:)/floe_rad_c(:))

! fsd tendency
do k = 1, nfsd
d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsd_tmp(k) &
* (c1/floe_rad_c(k) - tmp)
end do

! timestep required for this
subdt = get_subdt_fsd(nfsd, afsd_tmp(:), d_afsd_tmp(:))
subdt = MIN(subdt, dt)

! update fsd and elapsed time
afsd_tmp(:) = afsd_tmp(:) + subdt*d_afsd_tmp(:)
elapsed_t = elapsed_t + subdt


END DO

afsdn(:,n) = afsd_tmp(:)


end if ! aicen
end if ! rside > 0, otherwise do nothing

Expand Down Expand Up @@ -1185,7 +1193,7 @@ subroutine lateral_melt (dt, ncat, &
endif ! flag

if (tr_fsd) then
!call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )

! diagnostics
do k = 1, nfsd
Expand Down
Loading