Skip to content

Commit

Permalink
Update wavespec convergence algorithm (#305)
Browse files Browse the repository at this point in the history 
* update wavespec convergence algorithm, get rid of large allfraclengths array, minor cleanup of a few other things

* add fsd interface documentation and update documentation in general

* update documentation

* update documentation

* update documentation
  • Loading branch information
@apcraig
apcraig authored Mar 12, 2020
1 parent 11e4ec0 commit edb8c34
Show file tree
Hide file tree
Showing 8 changed files with 477 additions and 125 deletions.
60 changes: 38 additions & 22 deletions columnphysics/icepack_fsd.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,24 +71,25 @@ module icepack_fsd
contains

!=======================================================================
!
! Initialize ice fsd bounds (call whether or not restarting)
! Define the bounds, midpoints and widths of floe size
! categories in area and radius
!
! Note that radius widths cannot be larger than twice previous
! category width or floe welding will not have an effect
!
! Note also that the bound of the lowest floe size category is used
! to define the lead region width and the domain spacing for wave fracture
!
! authors: Lettie Roach, NIWA/VUW and C. M. Bitz, UW
!autodocument_start icepack_init_fsd_bounds
! Initialize ice fsd bounds (call whether or not restarting)
! Define the bounds, midpoints and widths of floe size
! categories in area and radius
!
! authors: Lettie Roach, NIWA/VUW and C. M. Bitz, UW

subroutine icepack_init_fsd_bounds(nfsd, &
floe_rad_l, & ! fsd size lower bound in m (radius)
floe_rad_c, & ! fsd size bin centre in m (radius)
floe_binwidth, & ! fsd size bin width in m (radius)
c_fsd_range) ! string for history output
c_fsd_range, & ! string for history output
write_diags ) ! flag for writing diagnostics

integer (kind=int_kind), intent(in) :: &
nfsd ! number of floe size categories
Expand All @@ -99,7 +100,12 @@ subroutine icepack_init_fsd_bounds(nfsd, &
floe_binwidth ! fsd size bin width in m (radius)

character (len=35), intent(out) :: &
c_fsd_range(nfsd) ! string for history output
c_fsd_range(nfsd) ! string for history output

logical (kind=log_kind), intent(in), optional :: &
write_diags ! write diags flag

!autodocument_end

! local variables

Expand All @@ -117,9 +123,17 @@ subroutine icepack_init_fsd_bounds(nfsd, &
real (kind=dbl_kind), dimension(:), allocatable :: &
lims

logical (kind=log_kind) :: &
l_write_diags ! local write diags

character(len=8) :: c_fsd1,c_fsd2
character(len=2) :: c_nf
character(len=*), parameter :: subname='(init_fsd_bounds)'
character(len=*), parameter :: subname='(icepack_init_fsd_bounds)'

l_write_diags = .true.
if (present(write_diags)) then
l_write_diags = write_diags
endif

if (nfsd.eq.24) then

Expand Down Expand Up @@ -212,8 +226,14 @@ subroutine icepack_init_fsd_bounds(nfsd, &
floe_rad(0) = floe_rad_l(1)
do n = 1, nfsd
floe_rad(n) = floe_rad_h(n)
! Save character string to write to history file
write (c_nf, '(i2)') n
write (c_fsd1, '(f6.3)') floe_rad(n-1)
write (c_fsd2, '(f6.3)') floe_rad(n)
c_fsd_range(n)=c_fsd1//'m < fsd Cat '//c_nf//' < '//c_fsd2//'m'
enddo

if (l_write_diags) then
write(warnstr,*) ' '
call icepack_warnings_add(warnstr)
write(warnstr,*) subname
Expand All @@ -223,26 +243,14 @@ subroutine icepack_init_fsd_bounds(nfsd, &
do n = 1, nfsd
write(warnstr,*) floe_rad(n-1),' < fsd Cat ',n, ' < ',floe_rad(n)
call icepack_warnings_add(warnstr)
! Write integer n to character string
write (c_nf, '(i2)') n

! Write floe_rad to character string
write (c_fsd1, '(f6.3)') floe_rad(n-1)
write (c_fsd2, '(f6.3)') floe_rad(n)

! Save character string to write to history file
c_fsd_range(n)=c_fsd1//'m < fsd Cat '//c_nf//' < '//c_fsd2//'m'
enddo

write(warnstr,*) ' '
call icepack_warnings_add(warnstr)
endif

end subroutine icepack_init_fsd_bounds

!=======================================================================
!
! Initialize the FSD
!
! When growing from no-ice conditions, initialize to zero.
! This allows the FSD to emerge, as described in Roach, Horvat et al. (2018)
!
Expand All @@ -254,6 +262,10 @@ end subroutine icepack_init_fsd_bounds
! sea ice floe size distribution. Journal of Geophysical Research: Oceans,
! 119(12), 8767–8777. doi:10.1002/2014JC010136
!
!autodocument_start icepack_init_fsd
!
! Initialize the FSD
!
! authors: Lettie Roach, NIWA/VUW

subroutine icepack_init_fsd(nfsd, ice_ic, &
Expand All @@ -274,6 +286,8 @@ subroutine icepack_init_fsd(nfsd, ice_ic, &
real (kind=dbl_kind), dimension (:), intent(inout) :: &
afsd ! floe size tracer: fraction distribution of floes

!autodocument_end

! local variables

real (kind=dbl_kind) :: alpha, totfrac
Expand Down Expand Up @@ -305,6 +319,7 @@ subroutine icepack_init_fsd(nfsd, ice_ic, &
end subroutine icepack_init_fsd

!=======================================================================
!autodocument_start icepack_cleanup_fsd
!
! Clean up small values and renormalize
!
Expand All @@ -319,6 +334,7 @@ subroutine icepack_cleanup_fsd (ncat, nfsd, afsdn)
real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
afsdn ! floe size distribution tracer

!autodocument_end
! local variables

integer (kind=int_kind) :: &
Expand Down
64 changes: 29 additions & 35 deletions columnphysics/icepack_wavefracspec.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ module icepack_wavefracspec
contains

!=======================================================================
!
!autodocument_start icepack_init_wave
! Initialize the wave spectrum and frequencies for the FSD
!
! authors: 2018 Lettie Roach, NIWA/VUW
Expand All @@ -82,6 +82,7 @@ subroutine icepack_init_wave(nfreq, &
wavefreq, & ! wave frequencies (s^-1)
dwavefreq ! wave frequency bin widths (s^-1)

!autodocument_end
! local variables
integer (kind=int_kind) :: k

Expand Down Expand Up @@ -174,6 +175,7 @@ function get_dafsd_wave(nfsd, afsd_init, fracture_hist, frac) &
end function get_dafsd_wave

!=======================================================================
!autodocument_start icepack_step_wavefracture
!
! Given fracture histogram computed from local wave spectrum, evolve
! the floe size distribution
Expand Down Expand Up @@ -226,6 +228,7 @@ subroutine icepack_step_wavefracture(wave_spec_type, &
real (kind=dbl_kind), dimension(nfsd,ncat) :: &
d_afsdn_wave ! change in fsd due to waves, per category

!autodocument_end
! local variables
integer (kind=int_kind) :: &
n, k, t, &
Expand Down Expand Up @@ -319,9 +322,8 @@ subroutine icepack_step_wavefracture(wave_spec_type, &
write(warnstr,*) subname, &
'warning: step_wavefracture struggling to converge'
call icepack_warnings_add(warnstr)
endif
endif


! required timestep
subdt = get_subdt_fsd(nfsd, afsd_tmp, d_afsd_tmp)
subdt = MIN(subdt, dt)
Expand Down Expand Up @@ -444,9 +446,6 @@ subroutine wave_frac(nfsd, nfreq, wave_spec_type, &
real (kind=dbl_kind), dimension(nx) :: &
fraclengths

real (kind=dbl_kind), dimension(max_no_iter*nx) :: &
allfraclengths

real (kind=dbl_kind), dimension(nx) :: &
X, & ! spatial domain (m)
eta ! sea surface height field (m)
Expand Down Expand Up @@ -480,19 +479,22 @@ subroutine wave_frac(nfsd, nfreq, wave_spec_type, &
! initialize frac lengths
fraclengths(:) = c0
prev_frac_local(:) = c0
frachistogram(:) = c0
fracerror = bignum

! loop while fracerror greater than error tolerance
DO iter = 1, loop_max_iter
iter = 0
do while (iter < loop_max_iter .and. fracerror > errortol)
iter = iter + 1

! Phase for each Fourier component may be constant or
! a random phase that varies in each i loop
! See documentation for discussion
if (trim(wave_spec_type).eq.'random') then
call RANDOM_NUMBER(rand_array)
if (icepack_warnings_aborted(subname)) return
call RANDOM_NUMBER(rand_array)
if (icepack_warnings_aborted(subname)) return
else
rand_array(:) = p5
rand_array(:) = p5
endif
phi = c2*pi*rand_array

Expand All @@ -510,29 +512,26 @@ subroutine wave_frac(nfsd, nfreq, wave_spec_type, &

! convert from diameter to radii
fraclengths(:) = fraclengths(:)/c2


if (ALL(fraclengths.lt.floe_rad_l(1))) then
frac_local(:) = c0
frac_local(:) = c0
else
frachistogram(:) = c0
allfraclengths((iter-1)*nx+1:(iter)*nx) = fraclengths(1:nx)

! bin into FS cats
! accumulate the frac histogram each iteration
do j = 1, size(fraclengths)
if (allfraclengths(j).gt.floe_rad_l(1)) then
do k = 1, nfsd-1
if ((allfraclengths(j) >= floe_rad_l(k)) .and. &
(allfraclengths(j) < floe_rad_l(k+1))) then
frachistogram(k) = frachistogram(k) + 1
if (fraclengths(j).gt.floe_rad_l(1)) then
do k = 1, nfsd-1
if ((fraclengths(j) >= floe_rad_l(k)) .and. &
(fraclengths(j) < floe_rad_l(k+1))) then
frachistogram(k) = frachistogram(k) + 1
end if
end do
if (fraclengths(j)>floe_rad_l(nfsd)) frachistogram(nfsd) = frachistogram(nfsd) + 1
end if
end do
if (allfraclengths(j)>floe_rad_l(nfsd)) frachistogram(nfsd) = frachistogram(nfsd) + 1
end if
end do

do k = 1, nfsd
frac_local(k) = floe_rad_c(k)*frachistogram(k)
frac_local(k) = floe_rad_c(k)*frachistogram(k)
end do

! normalize
Expand All @@ -546,22 +545,18 @@ subroutine wave_frac(nfsd, nfreq, wave_spec_type, &
! check avg frac local against previous iteration
fracerror = SUM(ABS(frac_local - prev_frac_local))/nfsd

if (fracerror.lt.errortol) EXIT

if (iter.gt.100) then
write(warnstr,*) subname, &
'warning: wave_frac struggling to converge'
call icepack_warnings_add(warnstr)
endif

! save histogram for next iteration
prev_frac_local = frac_local


end if

END DO

if (iter >= max_no_iter) then
write(warnstr,*) subname,'warning: wave_frac struggling to converge'
call icepack_warnings_add(warnstr)
endif

end subroutine wave_frac

!===========================================================================
Expand All @@ -577,7 +572,7 @@ end subroutine wave_frac
!
subroutine get_fraclengths(X, eta, fraclengths, hbar)

real (kind=dbl_kind) :: &
real (kind=dbl_kind), intent(in) :: &
hbar ! mean thickness (m)

real (kind=dbl_kind), intent(in), dimension (nx) :: &
Expand Down Expand Up @@ -686,7 +681,6 @@ subroutine get_fraclengths(X, eta, fraclengths, hbar)
delta_pos = X(j_pos) - X(j )
delta = X(j ) - X(j_neg)


! This equation differs from HT2015 by a factor 2 in numerator
! and eta(j_pos). This is the correct form of the equation.

Expand Down
23 changes: 12 additions & 11 deletions configuration/driver/icedrv_InitMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,18 +93,19 @@ subroutine icedrv_initialize
call icedrv_system_abort(file=__FILE__,line=__LINE__)
endif

if (tr_fsd) call icepack_init_fsd_bounds( &
nfsd=nfsd, & ! floe size distribution
floe_rad_l=floe_rad_l, & ! fsd size lower bound in m (radius)
floe_rad_c=floe_rad_c, & ! fsd size bin centre in m (radius)
floe_binwidth=floe_binwidth, & ! fsd size bin width in m (radius)
c_fsd_range=c_fsd_range) ! string for history output
call init_fsd

call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted(subname)) then
call icedrv_system_abort(file=__FILE__,line=__LINE__)
if (tr_fsd) then
call icepack_init_fsd_bounds( &
nfsd=nfsd, & ! floe size distribution
floe_rad_l=floe_rad_l, & ! fsd size lower bound in m (radius)
floe_rad_c=floe_rad_c, & ! fsd size bin centre in m (radius)
floe_binwidth=floe_binwidth, & ! fsd size bin width in m (radius)
c_fsd_range=c_fsd_range) ! string for history output
call icepack_warnings_flush(nu_diag)
if (icepack_warnings_aborted(subname)) then
call icedrv_system_abort(file=__FILE__,line=__LINE__)
endif
endif
call init_fsd

call calendar(time) ! determine the initial date

Expand Down
9 changes: 8 additions & 1 deletion doc/generate_interfaces.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ mv ${rstfile} ${rstfile}.orig
for file in ${inpfiles}; do

filename=`basename $file`
firstfileintfc=0

while IFS= read -r line; do
if [[ $line =~ .*$endline.* ]]; then
Expand All @@ -68,11 +69,17 @@ filename=`basename $file`
exit -9
fi
echo "$filename $title"
if [ $firstfileintfc = 0 ]; then
firstfileintfc=1
echo "" >> $rstfile
echo "${filename}" >> $rstfile
echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" >> $rstfile
fi
echo "" >> $rstfile
echo ".. _${title}:" >> $rstfile
echo "" >> $rstfile
echo "${title}" >> $rstfile
echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" >> $rstfile
echo "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^" >> $rstfile
echo ".. code-block:: fortran" >> $rstfile
echo "" >> $rstfile
fi
Expand Down
Loading

0 comments on commit edb8c34

Please sign in to comment.