From e6af1a27f7b21b9d9b4d40638066bcd283b9352f Mon Sep 17 00:00:00 2001 From: lettie_roach Date: Tue, 28 Jan 2020 17:47:27 -0700 Subject: [PATCH 1/8] Checking latm timestep --- columnphysics/icepack_therm_itd.F90 | 27 +++++++++++++++----------- columnphysics/icepack_wavefracspec.F90 | 26 +++++++++++-------------- 2 files changed, 27 insertions(+), 26 deletions(-) diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 1f7e3de93..ce9d8ab9d 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -50,7 +50,8 @@ 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 - + use icepack_wavefracspec, only: get_subdt_fsd + implicit none private @@ -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 G_r * areal mFSTD tilda + afsd_tmp, d_afsd_tmp real (kind=dbl_kind), dimension(nfsd+1) :: & f_flx ! @@ -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)' @@ -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 @@ -1104,6 +1106,15 @@ subroutine lateral_melt (dt, ncat, & print*,'sum(df_flx)/=0' tmp = SUM(afsdn_init(:,n)/floe_rad_c(:)) + do k = 1, nfsd + d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsdn_init(k,n) & + * (c1/floe_rad_c(k) - tmp) + end do + + subdt = get_subdt_fsd(nfsd, afsdn_init(:,n), d_afsd_tmp(:)) + subdt = MIN(subdt, dt) + print *, 'subdt =',subdt + do k = 1, nfsd afsdn (k,n) = afsdn_init(k,n) & + dt * (-df_flx(k) + c2 * G_radialn(n) * afsdn_init(k,n) & @@ -1118,12 +1129,6 @@ subroutine lateral_melt (dt, ncat, & 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 - end if ! aicen end if ! rside > 0, otherwise do nothing diff --git a/columnphysics/icepack_wavefracspec.F90 b/columnphysics/icepack_wavefracspec.F90 index 12c01eadb..f659cfaba 100644 --- a/columnphysics/icepack_wavefracspec.F90 +++ b/columnphysics/icepack_wavefracspec.F90 @@ -36,7 +36,7 @@ module icepack_wavefracspec implicit none private - public :: icepack_init_wave, icepack_step_wavefracture + public :: icepack_init_wave, icepack_step_wavefracture, get_subdt_fsd real (kind=dbl_kind), parameter :: & swh_minval = 0.01_dbl_kind, & ! minimum value of wave height (m) @@ -169,13 +169,13 @@ end function get_dafsd_wave !======================================================================= ! -! Adaptive timestepping for wave fracture +! Adaptive timestepping (process agnostic) ! See reference: Horvat & Tziperman (2017) JGR, Appendix A ! ! authors: 2018 Lettie Roach, NIWA/VUW ! ! - function get_subdt_wave(nfsd, afsd_init, d_afsd) & + function get_subdt_fsd(nfsd, afsd_init, d_afsd) & result(subdt) integer (kind=int_kind), intent(in) :: & @@ -202,7 +202,7 @@ function get_subdt_wave(nfsd, afsd_init, d_afsd) & subdt = MINVAL(check_dt) - end function get_subdt_wave + end function get_subdt_fsd !======================================================================= ! @@ -281,7 +281,8 @@ subroutine icepack_step_wavefracture(wave_spec_type, & afsd_tmp , & ! tracer array d_afsd_tmp ! change - character(len=*),parameter :: subname='(icepack_step_wavefracture)' + character(len=*),parameter :: & + subname='(icepack_step_wavefracture)' !------------------------------------ @@ -335,23 +336,18 @@ subroutine icepack_step_wavefracture(wave_spec_type, & DO WHILE (elapsed_t < dt) nsubt = nsubt + 1 + ! check in case wave fracture struggles to converge + if (nsubt>100) print *, & + 'wave frac taking a while to converge....' + ! if all floes in smallest category already, exit if (afsd_tmp(1).ge.c1-puny) EXIT ! calculate d_afsd using current afstd d_afsd_tmp = get_dafsd_wave(nfsd, afsd_tmp, fracture_hist, frac) - ! check in case wave fracture struggles to converge - if (nsubt>100) then - print *, 'afsd_tmp ',afsd_tmp - print *, 'dafsd_tmp ',d_afsd_tmp - print *, 'subt ',nsubt - print *, & - 'wave frac taking a while to converge....' - end if - ! required timestep - subdt = get_subdt_wave(nfsd, afsd_tmp, d_afsd_tmp) + subdt = get_subdt_fsd(nfsd, afsd_tmp, d_afsd_tmp) subdt = MIN(subdt, dt) ! update afsd From 1c86722a545defec8d35feb9c9fb2741d4bf20a5 Mon Sep 17 00:00:00 2001 From: lettie_roach Date: Fri, 31 Jan 2020 11:57:23 -0700 Subject: [PATCH 2/8] Debugging --- columnphysics/icepack_fsd.F90 | 9 +++++++++ columnphysics/icepack_therm_itd.F90 | 2 +- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90 index 7500bcebd..4c10ae683 100644 --- a/columnphysics/icepack_fsd.F90 +++ b/columnphysics/icepack_fsd.F90 @@ -655,6 +655,7 @@ subroutine fsd_add_new_ice (ncat, n, nfsd, & afsdn_latg ! fsd after lateral growth real (kind=dbl_kind), dimension (nfsd) :: & + dafsd_tmp, & ! df_flx , & ! finite differences for G_r*tilda(L) afsd_ni ! fsd after new ice added @@ -682,6 +683,14 @@ subroutine fsd_add_new_ice (ncat, n, nfsd, & ! if (abs(sum(df_flx)) > puny) print*,'fsd_add_new ERROR df_flx /= 0' afsdn_latg(:,n) = c0 + do k = 1, nfsd + dafsd_tmp(k) = (-df_flx(k) + c2 * G_radial * afsdn(k,n) & + * (c1/floe_rad_c(k) - SUM(afsdn(:,n)/floe_rad_c(:))) ) + + end do + + if (.NOT. ALL(dafsd_tmp.eq.c0)) print *, 'latg ',dafsd_tmp(:) + do k = 1, nfsd afsdn_latg(k,n) = afsdn(k,n) & + dt * (-df_flx(k) + c2 * G_radial * afsdn(k,n) & diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index ce9d8ab9d..c64746871 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -1110,10 +1110,10 @@ subroutine lateral_melt (dt, ncat, & d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsdn_init(k,n) & * (c1/floe_rad_c(k) - tmp) end do + print *, 'latm ',d_afsd_tmp(:) subdt = get_subdt_fsd(nfsd, afsdn_init(:,n), d_afsd_tmp(:)) subdt = MIN(subdt, dt) - print *, 'subdt =',subdt do k = 1, nfsd afsdn (k,n) = afsdn_init(k,n) & From 5800f933fc95b1b0f2776e24099ad96adcc6d983 Mon Sep 17 00:00:00 2001 From: lettie_roach Date: Wed, 5 Feb 2020 18:53:48 -0700 Subject: [PATCH 3/8] Added in adaptive timestepping for lateral melt and growth --- columnphysics/icepack_fsd.F90 | 73 +++++++++++++-------- columnphysics/icepack_therm_itd.F90 | 88 +++++++++++++++----------- columnphysics/icepack_wavefracspec.F90 | 6 +- 3 files changed, 101 insertions(+), 66 deletions(-) diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90 index 4c10ae683..2133db2d9 100644 --- a/columnphysics/icepack_fsd.F90 +++ b/columnphysics/icepack_fsd.F90 @@ -48,6 +48,7 @@ module icepack_fsd use icepack_tracers, only: nt_fsd, tr_fsd use icepack_warnings, only: warnstr, icepack_warnings_add use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted + use icepack_wavefracspec, only: get_subdt_fsd implicit none @@ -565,6 +566,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, @@ -649,53 +657,68 @@ 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) :: & - dafsd_tmp, & ! - 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 - dafsd_tmp(k) = (-df_flx(k) + c2 * G_radial * afsdn(k,n) & - * (c1/floe_rad_c(k) - SUM(afsdn(:,n)/floe_rad_c(:))) ) + 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 + end do + WHERE (abs(dafsd_tmp).lt.puny) dafsd_tmp = c0 - if (.NOT. ALL(dafsd_tmp.eq.c0)) print *, 'latg ',dafsd_tmp(:) + ! 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 - 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 + END DO call icepack_cleanup_fsdn (nfsd, afsdn_latg(:,n)) trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsdn_latg(:,n) diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index c64746871..8ed2bb29c 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -916,7 +916,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 @@ -943,7 +944,7 @@ 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) :: & @@ -1034,7 +1035,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) @@ -1050,18 +1051,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. !----------------------------------------------------------------- @@ -1092,34 +1083,55 @@ 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 - d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsdn_init(k,n) & + + ! 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 - print *, 'latm ',d_afsd_tmp(:) + end do + WHERE (abs(d_afsd_tmp).lt.puny) d_afsd_tmp = c0 + + ! timestep required for this + subdt = get_subdt_fsd(nfsd, afsd_tmp(:), d_afsd_tmp(:)) + subdt = MIN(subdt, dt) - subdt = get_subdt_fsd(nfsd, afsdn_init(:,n), 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 - 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 + + END DO + + afsdn(:,n) = afsd_tmp(:) if (abs(sum(afsdn(:,n))-c1) > puny) & print*,'lateral_melt E afsdn not normed',sum(df_flx), sum(afsdn(:,n))-c1 diff --git a/columnphysics/icepack_wavefracspec.F90 b/columnphysics/icepack_wavefracspec.F90 index f659cfaba..d275c6c20 100644 --- a/columnphysics/icepack_wavefracspec.F90 +++ b/columnphysics/icepack_wavefracspec.F90 @@ -219,7 +219,7 @@ subroutine icepack_step_wavefracture(wave_spec_type, & wave_spectrum, wavefreq, dwavefreq, & trcrn, d_afsd_wave) - use icepack_fsd, only: icepack_cleanup_fsd + character (len=char_len), intent(in) :: & wave_spec_type ! type of wave spectrum forcing @@ -308,7 +308,7 @@ subroutine icepack_step_wavefracture(wave_spec_type, & ! if fracture occurs if (MAXVAL(fracture_hist) > puny) then ! protect against small numerical errors - 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,:) ) do n = 1, ncat @@ -392,7 +392,7 @@ subroutine icepack_step_wavefracture(wave_spec_type, & ! update trcrn trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsd_tmp/SUM(afsd_tmp) - 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,:) ) ! for diagnostics From 648a3d270caf5809e0ad4986b44d9ec7bfe3d9e4 Mon Sep 17 00:00:00 2001 From: lettie_roach Date: Mon, 10 Feb 2020 13:12:39 -0700 Subject: [PATCH 4/8] Clean up --- columnphysics/icepack_fsd.F90 | 1 - columnphysics/icepack_therm_itd.F90 | 12 ++---------- 2 files changed, 2 insertions(+), 11 deletions(-) diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90 index 2133db2d9..b61a780b8 100644 --- a/columnphysics/icepack_fsd.F90 +++ b/columnphysics/icepack_fsd.F90 @@ -708,7 +708,6 @@ subroutine fsd_add_new_ice (ncat, n, nfsd, & * (c1/floe_rad_c(k) - SUM(afsdn_latg(:,n)/floe_rad_c(:))) ) end do - WHERE (abs(dafsd_tmp).lt.puny) dafsd_tmp = c0 ! timestep required for this subdt = get_subdt_fsd(nfsd, afsdn_latg(:,n), dafsd_tmp(:)) diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 8ed2bb29c..0c5ee100a 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -1118,7 +1118,6 @@ subroutine lateral_melt (dt, ncat, & d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsd_tmp(k) & * (c1/floe_rad_c(k) - tmp) end do - WHERE (abs(d_afsd_tmp).lt.puny) d_afsd_tmp = c0 ! timestep required for this subdt = get_subdt_fsd(nfsd, afsd_tmp(:), d_afsd_tmp(:)) @@ -1133,14 +1132,7 @@ subroutine lateral_melt (dt, ncat, & afsdn(:,n) = afsd_tmp(:) - 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' - - + end if ! aicen end if ! rside > 0, otherwise do nothing @@ -1202,7 +1194,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 From bc01a079b0d354feaf54a4e4bf993c616547ea2b Mon Sep 17 00:00:00 2001 From: lettie_roach Date: Mon, 10 Feb 2020 13:15:21 -0700 Subject: [PATCH 5/8] Checked diffs to master --- columnphysics/icepack_wavefracspec.F90 | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/columnphysics/icepack_wavefracspec.F90 b/columnphysics/icepack_wavefracspec.F90 index d275c6c20..8b3849bdf 100644 --- a/columnphysics/icepack_wavefracspec.F90 +++ b/columnphysics/icepack_wavefracspec.F90 @@ -219,7 +219,7 @@ subroutine icepack_step_wavefracture(wave_spec_type, & wave_spectrum, wavefreq, dwavefreq, & trcrn, d_afsd_wave) - + use icepack_fsd, only: icepack_cleanup_fsd character (len=char_len), intent(in) :: & wave_spec_type ! type of wave spectrum forcing @@ -281,8 +281,7 @@ subroutine icepack_step_wavefracture(wave_spec_type, & afsd_tmp , & ! tracer array d_afsd_tmp ! change - character(len=*),parameter :: & - subname='(icepack_step_wavefracture)' + character(len=*),parameter :: subname='(icepack_step_wavefracture)' !------------------------------------ @@ -308,7 +307,7 @@ subroutine icepack_step_wavefracture(wave_spec_type, & ! if fracture occurs if (MAXVAL(fracture_hist) > puny) then ! protect against small numerical errors - !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,:) ) do n = 1, ncat @@ -336,16 +335,21 @@ subroutine icepack_step_wavefracture(wave_spec_type, & DO WHILE (elapsed_t < dt) nsubt = nsubt + 1 - ! check in case wave fracture struggles to converge - if (nsubt>100) print *, & - 'wave frac taking a while to converge....' - ! if all floes in smallest category already, exit if (afsd_tmp(1).ge.c1-puny) EXIT ! calculate d_afsd using current afstd d_afsd_tmp = get_dafsd_wave(nfsd, afsd_tmp, fracture_hist, frac) + ! check in case wave fracture struggles to converge + if (nsubt>100) then + print *, 'afsd_tmp ',afsd_tmp + print *, 'dafsd_tmp ',d_afsd_tmp + print *, 'subt ',nsubt + print *, & + 'wave frac taking a while to converge....' + end if + ! required timestep subdt = get_subdt_fsd(nfsd, afsd_tmp, d_afsd_tmp) subdt = MIN(subdt, dt) @@ -392,7 +396,7 @@ subroutine icepack_step_wavefracture(wave_spec_type, & ! update trcrn trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsd_tmp/SUM(afsd_tmp) - !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,:) ) ! for diagnostics From a081aadd4e4d749116e22ecbc6e9805abe73d432 Mon Sep 17 00:00:00 2001 From: lettie_roach Date: Mon, 10 Feb 2020 16:23:18 -0700 Subject: [PATCH 6/8] Move get_subdt to icepack_fsd --- columnphysics/icepack_fsd.F90 | 41 ++++++++++++++++++++++-- columnphysics/icepack_therm_itd.F90 | 3 +- columnphysics/icepack_wavefracspec.F90 | 43 ++------------------------ 3 files changed, 43 insertions(+), 44 deletions(-) diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90 index b61a780b8..e6d6ca897 100644 --- a/columnphysics/icepack_fsd.F90 +++ b/columnphysics/icepack_fsd.F90 @@ -48,13 +48,12 @@ module icepack_fsd use icepack_tracers, only: nt_fsd, tr_fsd use icepack_warnings, only: warnstr, icepack_warnings_add use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted - use icepack_wavefracspec, only: get_subdt_fsd implicit none 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) @@ -1012,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 diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 0c5ee100a..41510ae0e 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -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 @@ -50,7 +50,6 @@ 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 - use icepack_wavefracspec, only: get_subdt_fsd implicit none diff --git a/columnphysics/icepack_wavefracspec.F90 b/columnphysics/icepack_wavefracspec.F90 index 8b3849bdf..16c1fe695 100644 --- a/columnphysics/icepack_wavefracspec.F90 +++ b/columnphysics/icepack_wavefracspec.F90 @@ -33,10 +33,11 @@ module icepack_wavefracspec use icepack_parameters, only: bignum, puny, gravit, pi use icepack_tracers, only: nt_fsd use icepack_warnings, only: warnstr, icepack_warnings_add - + use icepack_fsd + implicit none private - public :: icepack_init_wave, icepack_step_wavefracture, get_subdt_fsd + public :: icepack_init_wave, icepack_step_wavefracture real (kind=dbl_kind), parameter :: & swh_minval = 0.01_dbl_kind, & ! minimum value of wave height (m) @@ -167,43 +168,6 @@ function get_dafsd_wave(nfsd, afsd_init, fracture_hist, frac) & end function get_dafsd_wave -!======================================================================= -! -! 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 - !======================================================================= ! ! Given fracture histogram computed from local wave spectrum, evolve @@ -219,7 +183,6 @@ subroutine icepack_step_wavefracture(wave_spec_type, & wave_spectrum, wavefreq, dwavefreq, & trcrn, d_afsd_wave) - use icepack_fsd, only: icepack_cleanup_fsd character (len=char_len), intent(in) :: & wave_spec_type ! type of wave spectrum forcing From 0030e1fcdb560d229c9202436482ee06e2ac7682 Mon Sep 17 00:00:00 2001 From: Lettie Roach Date: Wed, 12 Feb 2020 09:50:21 -0800 Subject: [PATCH 7/8] Add note to docs on FSD adaptive timestepping --- doc/source/science_guide/sg_itd.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/doc/source/science_guide/sg_itd.rst b/doc/source/science_guide/sg_itd.rst index df0848c47..74c883a83 100755 --- a/doc/source/science_guide/sg_itd.rst +++ b/doc/source/science_guide/sg_itd.rst @@ -161,6 +161,9 @@ on this sub-grid-scale 1D domain. If the strain between successive extrema excee a critical value new floes are formed with diameters equal to the distance between the extrema. +Note that tendencies in the FSD are computed using adaptive timestepping to ensure that +the FSD is bounded by zero and one (see :cite:`Horvat17`). + Floe size categories are set in *init\_fsd\_bounds* using an exponential spacing, beginning at 0.5 m with the largest size resolved set by choice of :math:`N_f` (``nfsd``), the number of floe size categories. Icepack currently supports ``nfsd = 1, 12, 16, 24``. Although ``nfsd = 1`` tracks the same ice floe diameter as From 7f5630fcde76a4b6308130b4e0aa6a3fc3bfdef4 Mon Sep 17 00:00:00 2001 From: Lettie Roach Date: Wed, 12 Feb 2020 09:54:18 -0800 Subject: [PATCH 8/8] Add reference --- doc/source/master_list.bib | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/doc/source/master_list.bib b/doc/source/master_list.bib index 1bd3e6bfb..5cf827ff4 100644 --- a/doc/source/master_list.bib +++ b/doc/source/master_list.bib @@ -669,10 +669,26 @@ @article{Roach19 title = {{Advances in modelling interactions between sea ice and ocean surface waves}}, journal = {Journal of Advances in Modeling Earth Systems}, url = {http://doi.wiley.com/10.1029/2019MS001836}, +volume = {11}, +number = {12}, +pages = {4167-4181}, year={2019} } +@article{Horvat17, +author = "C. Horvat and E. Tziperman", +title = {The evolution of scaling laws in the sea ice floe size distribution}, +journal = {Journal of Geophysical Research: Oceans}, +volume = {122}, +number = {9}, +pages = {7630-7650}, +url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2016JC012573}, +year = {2017} +} + + + % **********************************************