Skip to content

Commit

Permalink
FSD fixes for conservation (#495)
Browse files Browse the repository at this point in the history 
This is a bug fix plus some rearranging for conservation in CICE. This impacts FSD cases as well as non-FSD cases.

There are two main aspects to the bug fix:

    Initial ice and snow volume values at the beginning of the lateral melt routine are now used for updating the snow and ice enthalpy as well other tracers for lateral melting. This is answer changing in all cases!

    The other fix is to move the computation of vi0new_lat (lateral growth from FSD) outside of the subroutine fsd_lateral_growth and a check to make sure the category vi0new (vin0new) does not exceed the total vi0new. This is a warning and it recomputes to ensure conservation. This is only answer changing in FSD cases.

    There is also a change in definition of floe_area_c, so that it is precisely half way between floe_area_l and floe_area_h. The original computation of floe_area_c based on floe_rad_c biases the area towards the lower value.

Update the fsd subroutine arguments (floe_rad_c, floe_rad_l, floe_binwidth, c_fsd_range) to store static data inside Icepack when computed in Icepack rather than pass them back and forth. Provide an ability to pass out the values from Icepack as optional arguments.

Move the computation of rsiden and get rid of fside which is no longer needed.

Clean up some comments and indentation

Update swccsm3 test, set ktherm=1

Update interface documentation

There is an associated CICE tag to go with this.
  • Loading branch information
@dabail10
dabail10 authored Oct 31, 2024
1 parent 6da5668 commit 286630f
Show file tree
Hide file tree
Showing 12 changed files with 408 additions and 447 deletions.
140 changes: 82 additions & 58 deletions columnphysics/icepack_fsd.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,13 +55,19 @@ module icepack_fsd
public :: icepack_init_fsd_bounds, icepack_init_fsd, icepack_cleanup_fsd, &
fsd_lateral_growth, fsd_add_new_ice, fsd_weld_thermo, get_subdt_fsd

real(kind=dbl_kind), dimension(:), allocatable :: &
real(kind=dbl_kind), dimension(:), allocatable, public :: &
floe_rad_h, & ! fsd size higher bound in m (radius)
floe_rad_c, & ! fsd size center in m (radius)
floe_rad_l, & ! fsd size lower bound in m (radius)
floe_binwidth, & ! fsd size binwidth in m (radius)
floe_area_l, & ! fsd area at lower bound (m^2)
floe_area_h, & ! fsd area at higher bound (m^2)
floe_area_c, & ! fsd area at bin centre (m^2)
floe_area_binwidth ! floe area bin width (m^2)

character (len=35), dimension(:), allocatable, public :: &
c_fsd_range ! string for history output

integer(kind=int_kind), dimension(:,:), allocatable, public :: &
iweld ! floe size categories that can combine
! during welding (dimensionless)
Expand All @@ -85,22 +91,22 @@ module icepack_fsd
! authors: Lettie Roach, NIWA/VUW and C. M. Bitz, UW

subroutine icepack_init_fsd_bounds( &
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
write_diags ) ! flag for writing diagnostics
floe_rad_l_out, & ! fsd size lower bound in m (radius)
floe_rad_c_out, & ! fsd size bin centre in m (radius)
floe_binwidth_out, & ! fsd size bin width in m (radius)
c_fsd_range_out, & ! string for history output
write_diags) ! flag for writing diagnostics

real(kind=dbl_kind), dimension(:), intent(inout) :: &
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)
real(kind=dbl_kind), dimension(:), intent(out), optional :: &
floe_rad_l_out, & ! fsd size lower bound in m (radius)
floe_rad_c_out, & ! fsd size bin centre in m (radius)
floe_binwidth_out ! fsd size bin width in m (radius)

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

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

!autodocument_end

Expand All @@ -111,11 +117,8 @@ subroutine icepack_init_fsd_bounds( &

real (kind=dbl_kind) :: test

real (kind=dbl_kind), dimension (0:nfsd) :: &
floe_rad

real (kind=dbl_kind), dimension(:), allocatable :: &
lims
lims, floe_rad

character(len=8) :: c_fsd1,c_fsd2
character(len=2) :: c_nf
Expand Down Expand Up @@ -169,10 +172,15 @@ subroutine icepack_init_fsd_bounds( &

allocate( &
floe_rad_h (nfsd), & ! fsd size higher bound in m (radius)
floe_rad_l (nfsd), & ! fsd size lower bound in m (radius)
floe_rad_c (nfsd), & ! fsd size center in m (radius)
floe_rad (0:nfsd), & ! fsd bounds in m (radius)
floe_area_l (nfsd), & ! fsd area at lower bound (m^2)
floe_area_h (nfsd), & ! fsd area at higher bound (m^2)
floe_area_c (nfsd), & ! fsd area at bin centre (m^2)
floe_area_binwidth (nfsd), & ! floe area bin width (m^2)
floe_binwidth (nfsd), & ! floe bin width (m)
c_fsd_range (nfsd), & !
iweld (nfsd, nfsd), & ! fsd categories that can weld
stat=ierr)
if (ierr/=0) then
Expand All @@ -186,8 +194,11 @@ subroutine icepack_init_fsd_bounds( &
floe_rad_c = (floe_rad_h+floe_rad_l)/c2

floe_area_l = c4*floeshape*floe_rad_l**2
floe_area_c = c4*floeshape*floe_rad_c**2
floe_area_h = c4*floeshape*floe_rad_h**2
! floe_area_c = c4*floeshape*floe_rad_c**2
! This is exactly in the middle of floe_area_h and floe_area_l
! Whereas the above calculation is closer to floe_area_l.
floe_area_c = (floe_area_h+floe_area_l)/c2

floe_binwidth = floe_rad_h - floe_rad_l

Expand Down Expand Up @@ -220,20 +231,56 @@ subroutine icepack_init_fsd_bounds( &
enddo

if (present(write_diags)) then
if (write_diags) then
write(warnstr,*) ' '
call icepack_warnings_add(warnstr)
write(warnstr,*) subname
call icepack_warnings_add(warnstr)
write(warnstr,*) 'floe_rad(n-1) < fsd Cat n < floe_rad(n)'
call icepack_warnings_add(warnstr)
do n = 1, nfsd
write(warnstr,*) floe_rad(n-1),' < fsd Cat ',n, ' < ',floe_rad(n)
if (write_diags) then
write(warnstr,*) ' '
call icepack_warnings_add(warnstr)
enddo
write(warnstr,*) ' '
call icepack_warnings_add(warnstr)
write(warnstr,*) subname
call icepack_warnings_add(warnstr)
write(warnstr,*) 'floe_rad(n-1) < fsd Cat n < floe_rad(n)'
call icepack_warnings_add(warnstr)
do n = 1, nfsd
write(warnstr,*) floe_rad(n-1),' < fsd Cat ',n, ' < ',floe_rad(n)
call icepack_warnings_add(warnstr)
enddo
write(warnstr,*) ' '
call icepack_warnings_add(warnstr)
endif
endif

if (present(floe_rad_l_out)) then
if (size(floe_rad_l_out) /= size(floe_rad_l)) then
call icepack_warnings_add(subname//' floe_rad_l_out incorrect size')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
floe_rad_l_out(:) = floe_rad_l(:)
endif

if (present(floe_rad_c_out)) then
if (size(floe_rad_c_out) /= size(floe_rad_c)) then
call icepack_warnings_add(subname//' floe_rad_c_out incorrect size')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
floe_rad_c_out(:) = floe_rad_c(:)
endif

if (present(floe_binwidth_out)) then
if (size(floe_binwidth_out) /= size(floe_binwidth)) then
call icepack_warnings_add(subname//' floe_binwidth_out incorrect size')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
floe_binwidth_out(:) = floe_binwidth(:)
endif

if (present(c_fsd_range_out)) then
if (size(c_fsd_range_out) /= size(c_fsd_range)) then
call icepack_warnings_add(subname//' c_fsd_range_out incorrect size')
call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
return
endif
c_fsd_range_out(:) = c_fsd_range(:)
endif

end subroutine icepack_init_fsd_bounds
Expand All @@ -256,18 +303,11 @@ end subroutine icepack_init_fsd_bounds
!
! authors: Lettie Roach, NIWA/VUW

subroutine icepack_init_fsd(ice_ic, &
floe_rad_c, & ! fsd size bin centre in m (radius)
floe_binwidth, & ! fsd size bin width in m (radius)
afsd) ! floe size distribution tracer
subroutine icepack_init_fsd(ice_ic, afsd) ! floe size distribution tracer

character(len=char_len_long), intent(in) :: &
ice_ic ! method of ice cover initialization

real(kind=dbl_kind), dimension(:), intent(inout) :: &
floe_rad_c, & ! fsd size bin centre in m (radius)
floe_binwidth ! fsd size bin width in m (radius)

real (kind=dbl_kind), dimension (:), intent(inout) :: &
afsd ! floe size tracer: fraction distribution of floes

Expand Down Expand Up @@ -323,7 +363,6 @@ subroutine icepack_cleanup_fsd (afsdn)

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


if (tr_fsd) then

do n = 1, ncat
Expand Down Expand Up @@ -378,14 +417,11 @@ end subroutine icepack_cleanup_fsdn
!
! authors: Lettie Roach, NIWA/VUW

subroutine partition_area (floe_rad_c, aice, &
subroutine partition_area (aice, &
aicen, vicen, &
afsdn, lead_area, &
latsurf_area)

real (kind=dbl_kind), dimension(:), intent(in) :: &
floe_rad_c ! fsd size bin centre in m (radius)

real (kind=dbl_kind), intent(in) :: &
aice ! ice concentration

Expand Down Expand Up @@ -476,7 +512,7 @@ end subroutine partition_area
subroutine fsd_lateral_growth (dt, aice, &
aicen, vicen, &
vi0new, &
frazil, floe_rad_c, &
frazil, &
afsdn, &
lead_area, latsurf_area, &
G_radial, d_an_latg, &
Expand All @@ -497,10 +533,6 @@ subroutine fsd_lateral_growth (dt, aice, &
vi0new , & ! volume of new ice added to cat 1 (m)
frazil ! frazil ice growth (m/step-->cm/day)

! floe size distribution
real (kind=dbl_kind), dimension (:), intent(in) :: &
floe_rad_c ! fsd size bin centre in m (radius)

real (kind=dbl_kind), dimension(ncat), intent(out) :: &
d_an_latg ! change in aicen occuring due to lateral growth

Expand Down Expand Up @@ -529,7 +561,7 @@ subroutine fsd_lateral_growth (dt, aice, &
d_an_latg = c0

! partition volume into lateral growth and frazil
call partition_area (floe_rad_c, aice, &
call partition_area (aice, &
aicen, vicen, &
afsdn, lead_area, &
latsurf_area)
Expand All @@ -540,9 +572,6 @@ subroutine fsd_lateral_growth (dt, aice, &
vi0new_lat = vi0new * lead_area / (c1 + aice/latsurf_area)
end if

! for history/diagnostics
frazil = vi0new - vi0new_lat

! lateral growth increment
if (vi0new_lat > puny) then
G_radial = vi0new_lat/dt
Expand All @@ -563,8 +592,6 @@ subroutine fsd_lateral_growth (dt, aice, &
endif ! vi0new_lat > 0

! Use remaining ice volume as in standard model,
! but ice cannot grow into the area that has grown laterally
vi0new = vi0new - vi0new_lat
tot_latg = SUM(d_an_latg(:))

end subroutine fsd_lateral_growth
Expand All @@ -589,7 +616,6 @@ end subroutine fsd_lateral_growth
subroutine fsd_add_new_ice (n, &
dt, ai0new, &
d_an_latg, d_an_newi, &
floe_rad_c, floe_binwidth, &
G_radial, area2, &
wave_sig_ht, &
wave_spectrum, &
Expand Down Expand Up @@ -623,9 +649,7 @@ subroutine fsd_add_new_ice (n, &

real (kind=dbl_kind), dimension (:), intent(in) :: &
aicen_init , & ! fractional area of ice
aicen , & ! after update
floe_rad_c , & ! fsd size bin centre in m (radius)
floe_binwidth ! fsd size bin width in m (radius)
aicen ! after update

real (kind=dbl_kind), dimension (:,:), intent(in) :: &
afsdn ! floe size distribution tracer
Expand Down
Loading

0 comments on commit 286630f

Please sign in to comment.