diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90 new file mode 100644 index 000000000..057cca65a --- /dev/null +++ b/columnphysics/icepack_fsd.F90 @@ -0,0 +1,1025 @@ +!========================================================================= +! +! This module contains the subroutines required to define +! a floe size distribution tracer for sea ice +! +! Theory based on: +! +! Horvat, C., & Tziperman, E. (2015). A prognostic model of the sea-ice +! floe size and thickness distribution. The Cryosphere, 9(6), 2119–2134. +! doi:10.5194/tc-9-2119-2015 +! +! and implementation described in: +! +! Roach, L. A., Horvat, C., Dean, S. M., & Bitz, C. M. (2018). An emergent +! sea ice floe size distribution in a global coupled ocean--sea ice model. +! Journal of Geophysical Research: Oceans, 123(6), 4322–4337. +! doi:10.1029/2017JC013692 +! +! with some modifications. +! +! For floe welding parameter and tensile mode parameter values, see +! +! Roach, L. A., Smith, M. M., & Dean, S. M. (2018). Quantifying +! growth of pancake sea ice floes using images from drifting buoys. +! Journal of Geophysical Research: Oceans, 123(4), 2851–2866. +! doi: 10.1002/2017JC013693 +! +! +! Variable naming convention +! for k = 1, nfsd and n = 1, ncat +! afsdn(k,n) = trcrn(:,:,nt_nfsd+k-1,n,:) +! afsd (k) is per thickness category or averaged over n +! afs is the associated scalar value for (k,n) +! +! authors: Lettie Roach, VUW/NIWA +! C. M. Bitz, UW +! +! 2016: CMB started +! 2016-8: LR worked on most of it +! 2019: ECH ported to Icepack + +!----------------------------------------------------------------- + + module icepack_fsd + + use icepack_kinds + use icepack_parameters, only: c0, c1, c2, c3, c4, p01, p1, p5, puny + use icepack_parameters, only: pi, floeshape, wave_spec, bignum, gravit, rhoi + 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 + + implicit none + + private + public :: icepack_init_fsd_bounds, icepack_init_fsd, icepack_cleanup_fsd, & + fsd_lateral_growth, fsd_add_new_ice, fsd_weld_thermo + + real(kind=dbl_kind), dimension(:), allocatable :: & + floe_rad_h, & ! fsd size higher bound 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) + + integer(kind=int_kind), dimension(:,:), allocatable, public :: & + iweld ! floe size categories that can combine + ! during welding (dimensionless) + +!======================================================================= + + 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 +! + 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 + + integer (kind=int_kind), intent(in) :: & + nfsd ! number of floe size categories + + 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) + + character (len=35), intent(out) :: & + c_fsd_range(nfsd) ! string for history output + + ! local variables + + integer (kind=int_kind) :: n, m, k + integer (kind=int_kind) :: ierr + + real (kind=dbl_kind) :: test + + real (kind=dbl_kind), dimension (nfsd+1) :: & + area_lims, area_lims_scaled + + real (kind=dbl_kind), dimension (0:nfsd) :: & + floe_rad + + real (kind=dbl_kind), dimension(:), allocatable :: & + lims + + character(len=8) :: c_fsd1,c_fsd2 + character(len=2) :: c_nf + character(len=*), parameter :: subname='(init_fsd_bounds)' + + if (nfsd.eq.24) then + + allocate(lims(24+1)) + + lims = (/ 6.65000000e-02, 5.31030847e+00, 1.42865861e+01, 2.90576686e+01, & + 5.24122136e+01, 8.78691405e+01, 1.39518470e+02, 2.11635752e+02, & + 3.08037274e+02, 4.31203059e+02, 5.81277225e+02, 7.55141047e+02, & + 9.45812834e+02, 1.34354446e+03, 1.82265364e+03, 2.47261361e+03, & + 3.35434988e+03, 4.55051413e+03, 6.17323164e+03, 8.37461170e+03, & + 1.13610059e+04, 1.54123510e+04, 2.09084095e+04, 2.83643675e+04, & + 3.84791270e+04 /) + + elseif (nfsd.eq.16) then + + allocate(lims(16+1)) + + lims = (/ 6.65000000e-02, 5.31030847e+00, 1.42865861e+01, 2.90576686e+01, & + 5.24122136e+01, 8.78691405e+01, 1.39518470e+02, 2.11635752e+02, & + 3.08037274e+02, 4.31203059e+02, 5.81277225e+02, 7.55141047e+02, & + 9.45812834e+02, 1.34354446e+03, 1.82265364e+03, 2.47261361e+03, & + 3.35434988e+03 /) + + else if (nfsd.eq.12) then + + allocate(lims(12+1)) + + lims = (/ 6.65000000e-02, 5.31030847e+00, 1.42865861e+01, 2.90576686e+01, & + 5.24122136e+01, 8.78691405e+01, 1.39518470e+02, 2.11635752e+02, & + 3.08037274e+02, 4.31203059e+02, 5.81277225e+02, 7.55141047e+02, & + 9.45812834e+02 /) + + else if (nfsd.eq.1) then ! default case + + allocate(lims(1+1)) + + lims = (/ 6.65000000e-02, 3.0e+02 /) + + else + + call icepack_warnings_add(subname//& + ' floe size categories not defined for nfsd') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + + end if + + allocate( & + floe_rad_h (nfsd), & ! fsd size higher bound 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) + iweld (nfsd, nfsd), & ! fsd categories that can weld + stat=ierr) + if (ierr/=0) then + call icepack_warnings_add(subname//' Out of Memory fsd') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + + floe_rad_l = lims(1:nfsd ) + floe_rad_h = lims(2:nfsd+1) + floe_rad_c = (floe_rad_h+floe_rad_l)/c2 + + floe_area_l = c4*floeshape*floe_rad_l**c2 + floe_area_c = c4*floeshape*floe_rad_c**c2 + floe_area_h = c4*floeshape*floe_rad_h**c2 + + floe_binwidth = floe_rad_h - floe_rad_l + + floe_area_binwidth = floe_area_h - floe_area_l + + ! floe size categories that can combine during welding + iweld(:,:) = -999 + do n = 1, nfsd + do m = 1, nfsd + test = floe_area_c(n) + floe_area_c(m) + + do k = 1, nfsd-1 + if ((test >= floe_area_l(k)) .and. (test < floe_area_h(k))) & + iweld(n,m) = k + end do + if (test >= floe_area_l(nfsd)) iweld(n,m) = nfsd + + + end do + end do + + if (allocated(lims)) deallocate(lims) + + ! write fsd bounds + floe_rad(0) = floe_rad_l(1) + do n = 1, nfsd + floe_rad(n) = floe_rad_h(n) + 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) + ! 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) + + 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) +! +! Otherwise initalize with a power law, following Perovich +! & Jones (2014). The basin-wide applicability of such a +! prescribed power law has not yet been tested. +! +! Perovich, D. K., & Jones, K. F. (2014). The seasonal evolution of +! sea ice floe size distribution. Journal of Geophysical Research: Oceans, +! 119(12), 8767–8777. doi:10.1002/2014JC010136 +! +! authors: Lettie Roach, NIWA/VUW + + subroutine icepack_init_fsd(nfsd, 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 + + integer(kind=int_kind), intent(in) :: & + nfsd + + 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 + + ! local variables + + real (kind=dbl_kind) :: alpha, totfrac + + integer (kind=int_kind) :: k + + real (kind=dbl_kind), dimension (nfsd) :: & + num_fsd ! number distribution of floes + + if (trim(ice_ic) == 'none') then + + afsd(:) = c0 + + else ! Perovich (2014) + + ! fraction of ice in each floe size and thickness category + ! same for ALL cells (even where no ice) initially + alpha = 2.1_dbl_kind + totfrac = c0 ! total fraction of floes + do k = 1, nfsd + num_fsd(k) = (2*floe_rad_c(k))**(-alpha-c1) ! number distribution of floes + afsd (k) = num_fsd(k)*floe_area_c(k)*floe_binwidth(k) ! fraction distribution of floes + totfrac = totfrac + afsd(k) + enddo + afsd = afsd/totfrac ! normalize + + endif ! ice_ic + + end subroutine icepack_init_fsd + +!======================================================================= +! +! Clean up small values and renormalize +! +! authors: Elizabeth Hunke, LANL +! + subroutine icepack_cleanup_fsd (ncat, nfsd, afsdn) + + integer (kind=int_kind), intent(in) :: & + ncat , & ! number of thickness categories + nfsd ! number of floe size categories + + real (kind=dbl_kind), dimension(:,:), intent(inout) :: & + afsdn ! floe size distribution tracer + + ! local variables + + integer (kind=int_kind) :: & + n ! thickness category index + + if (tr_fsd) then + + do n = 1, ncat + call icepack_cleanup_fsdn(nfsd, afsdn(:,n)) + enddo + + endif ! tr_fsd + + end subroutine icepack_cleanup_fsd + +!======================================================================= +! +! Clean up small values and renormalize -- per category +! +! authors: Elizabeth Hunke, LANL +! + + subroutine icepack_cleanup_fsdn (nfsd, afsd) + + integer (kind=int_kind), intent(in) :: & + nfsd ! number of floe size categories + + real (kind=dbl_kind), dimension(:), intent(inout) :: & + afsd ! floe size distribution tracer + + ! local variables + + integer (kind=int_kind) :: & + k ! floe size category index + + real (kind=dbl_kind) :: & + tot + + do k = 1, nfsd + if (afsd(k) < puny) afsd(k) = c0 + enddo + + tot = sum(afsd(:)) + if (tot > puny) then + do k = 1, nfsd + afsd(k) = afsd(k) / tot ! normalize + enddo + else ! represents ice-free cell, so set to zero + afsd(:) = c0 + endif + + end subroutine icepack_cleanup_fsdn + +!======================================================================= +! +! Given the joint ice thickness and floe size distribution, calculate +! the lead region and the total lateral surface area following Horvat +! and Tziperman (2015). +! +! authors: Lettie Roach, NIWA/VUW + + subroutine partition_area (ncat, nfsd, & + floe_rad_c, aice, & + aicen, vicen, & + afsdn, lead_area, & + latsurf_area) + + integer (kind=int_kind), intent(in) :: & + ncat , & ! number of thickness categories + nfsd ! number of floe size categories + + 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 + + real (kind=dbl_kind), dimension(:), intent(in) :: & + aicen , & ! fractional area of ice + vicen ! volume per unit area of ice (m) + + real (kind=dbl_kind), dimension(:,:), intent(in) :: & + afsdn ! floe size distribution tracer + + real (kind=dbl_kind), intent(out) :: & + lead_area , & ! the fractional area of the lead region + latsurf_area ! the area covered by lateral surface of floes + + ! local variables + + integer (kind=int_kind) :: & + n , & ! thickness category index + k ! floe size index + + real (kind=dbl_kind) :: & + width_leadreg, & ! width of lead region (m) + thickness ! actual thickness of ice in thickness cat (m) + + !----------------------------------------------------------------- + ! Initialize + !----------------------------------------------------------------- + lead_area = c0 + latsurf_area = c0 + + ! Set the width of the lead region to be the smallest + ! floe size category, as per Horvat & Tziperman (2015) + width_leadreg = floe_rad_c(1) + + ! Only calculate these areas where there is some ice + if (aice > puny) then + + ! lead area = sum of areas of annuli around all floes + do n = 1, ncat + do k = 1, nfsd + lead_area = lead_area + aicen(n) * afsdn(k,n) & + * (c2*width_leadreg /floe_rad_c(k) & + + width_leadreg**c2/floe_rad_c(k)**2) + enddo ! k + enddo ! n + + ! cannot be greater than the open water fraction + lead_area=MIN(lead_area, c1-aice) + if (lead_area < c0) then + if (lead_area < -puny) then +! stop 'lead_area lt0 in partition_area' + else + lead_area=MAX(lead_area,c0) + end if + end if + + ! Total fractional lateral surface area in each grid (per unit ocean area) + do n = 1, ncat + thickness = c0 + if (aicen(n) > c0) thickness = vicen(n)/aicen(n) + + do k = 1, nfsd + latsurf_area = latsurf_area & + + afsdn(k,n) * aicen(n) * c2 * thickness/floe_rad_c(k) + end do + end do + + ! check +! if (latsurf_area < c0) stop 'negative latsurf_area' +! if (latsurf_area /= latsurf_area) stop 'latsurf_area NaN' + + end if ! aice + + end subroutine partition_area + +!======================================================================= +! +! Lateral growth at the edges of floes +! +! Compute the portion of new ice growth that occurs at the edges of +! floes. The remainder will grow as new ice frazil ice in open water +! (away from floes). +! +! See Horvat & Tziperman (2015) and Roach, Horvat et al. (2018). +! +! authors: Lettie Roach, NIWA/VUW +! + subroutine fsd_lateral_growth (ncat, nfsd, & + dt, aice, & + aicen, vicen, & + vi0new, & + frazil, floe_rad_c, & + afsdn, & + lead_area, latsurf_area, & + G_radial, d_an_latg, & + tot_latg) + + integer (kind=int_kind), intent(in) :: & + ncat , & ! number of thickness categories + nfsd ! number of floe size categories + + real (kind=dbl_kind), intent(in) :: & + dt , & ! time step (s) + aice ! total concentration of ice + + real (kind=dbl_kind), dimension (:), intent(in) :: & + aicen , & ! concentration of ice + vicen ! volume per unit area of ice (m) + + real (kind=dbl_kind), dimension(:,:), intent(in) :: & + afsdn ! floe size distribution tracer + + real (kind=dbl_kind), intent(inout) :: & + 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 + + real (kind=dbl_kind), intent(out) :: & + G_radial , & ! lateral melt rate (m/s) + tot_latg ! total change in aice due to + ! lateral growth at the edges of floes + + + ! local variables + + integer (kind=int_kind) :: & + n, k ! ice category indices + + real (kind=dbl_kind) :: & + vi0new_lat ! volume of new ice added laterally to FSD (m) + + real (kind=dbl_kind), intent(out) :: & + lead_area , & ! the fractional area of the lead region + latsurf_area ! the area covered by lateral surface of floes + + character(len=*),parameter :: subname='(fsd_lateral_growth)' + + lead_area = c0 + latsurf_area = c0 + G_radial = c0 + tot_latg = c0 + d_an_latg = c0 + + ! partition volume into lateral growth and frazil + call partition_area (ncat, nfsd, & + floe_rad_c, aice, & + aicen, vicen, & + afsdn, lead_area, & + latsurf_area) + + vi0new_lat = c0 + if (latsurf_area > puny) then + 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 + do n = 1, ncat + + do k = 1, nfsd + d_an_latg(n) = d_an_latg(n) & + + c2*aicen(n)*afsdn(k,n)*G_radial*dt/floe_rad_c(k) + end do + end do ! n + 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 + +!======================================================================= +! +! Evolve the FSD subject to lateral growth and the growth of new ice +! in thickness category 1. +! +! If ocean surface wave forcing is provided, the floe size of new ice +! (grown away from floe edges) can computed from the wave field +! based on the tensile (stretching) stress limitation following +! Shen et al. (2001). Otherwise, new floes all grow in the smallest +! floe size category, representing pancake ice formation. +! +! Shen, H., Ackley, S., & Hopkins, M. (2001). A conceptual model +! for pancake-ice formation in a wave field. +! Annals of Glaciology, 33, 361-367. doi:10.3189/172756401781818239 +! +! authors: Lettie Roach, NIWA/VUW +! + + subroutine fsd_add_new_ice (ncat, n, nfsd, & + dt, ai0new, & + d_an_latg, d_an_newi, & + floe_rad_c, floe_binwidth, & + G_radial, area2, & + wave_sig_ht, & + wave_spectrum, & + wavefreq, & + dwavefreq, & + d_afsd_latg, & + d_afsd_newi, & + afsdn, aicen_init, & + aicen, trcrn) + + integer (kind=int_kind), intent(in) :: & + n , & ! thickness category number + ncat , & ! number of thickness categories + nfsd ! number of floe size categories + + real (kind=dbl_kind), intent(in) :: & + dt , & ! time step (s) + ai0new , & ! area of new ice added to cat 1 + G_radial , & ! lateral melt rate (m/s) + wave_sig_ht ! wave significant height (everywhere) (m) + + real (kind=dbl_kind), dimension(:), intent(in) :: & + wave_spectrum ! ocean surface wave spectrum as a function of frequency + ! power spectral density of surface elevation, E(f) (units m^2 s) + + real(kind=dbl_kind), dimension(:), intent(in) :: & + wavefreq, & ! wave frequencies (s^-1) + dwavefreq ! wave frequency bin widths (s^-1) + + real (kind=dbl_kind), dimension(:), intent(in) :: & + d_an_latg, d_an_newi ! change in aicen due to + ! lateral growth and frazil + ! ice formation + + 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) + + real (kind=dbl_kind), dimension (:,:), intent(in) :: & + afsdn ! floe size distribution tracer + + real (kind=dbl_kind), dimension (:), intent(in) :: & + area2 ! area after lateral growth and before new ice formation + + real (kind=dbl_kind), dimension (:,:), intent(inout) :: & + trcrn ! ice tracers + + real (kind=dbl_kind), dimension(:), intent(inout) :: & + ! change in floe size distribution (area) + d_afsd_latg , & ! due to fsd lateral growth + d_afsd_newi ! new ice formation + + integer (kind=int_kind) :: & + k ! floe size category index + + 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) + 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 + +! 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 + + call icepack_cleanup_fsdn (nfsd, afsdn_latg(:,n)) + trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsdn_latg(:,n) + + end if ! lat growth + + new_size = nfsd + if (n == 1) then + ! add new frazil ice to smallest thickness + if (d_an_newi(n) > puny) then + + afsd_ni(:) = c0 + + if (SUM(afsdn_latg(:,n)) > puny) then ! fsd exists + + if (wave_spec) then + if (wave_sig_ht > puny) & + call wave_dep_growth (nfsd, wave_spectrum, & + wavefreq, dwavefreq, & + new_size) + + ! grow in new_size category + afsd_ni(new_size) = (afsdn_latg(new_size,n)*area2(n) + ai0new) & + / (area2(n) + ai0new) + do k = 1, new_size-1 ! diminish other floe cats accordingly + afsd_ni(k) = afsdn_latg(k,n)*area2(n) / (area2(n) + ai0new) + end do + do k = new_size+1, nfsd ! diminish other floe cats accordingly + afsd_ni(k) = afsdn_latg(k,n)*area2(n) / (area2(n) + ai0new) + end do + + else ! grow in smallest floe size category + afsd_ni(1) = (afsdn_latg(1,n)*area2(n) + ai0new) & + / (area2(n) + ai0new) + do k = 2, nfsd ! diminish other floe cats accordingly + afsd_ni(k) = afsdn_latg(k,n)*area2(n) / (area2(n)+ai0new) + enddo + end if ! wave spec + + else ! no fsd, so entirely new ice + + if (wave_spec) then + if (wave_sig_ht > puny) & + call wave_dep_growth (nfsd, wave_spectrum, & + wavefreq, dwavefreq, & + new_size) + afsd_ni(new_size) = c1 + else + afsd_ni(1) = c1 + endif ! wave forcing + + endif ! entirely new ice or not + + trcrn(nt_fsd:nt_fsd+nfsd-1,n) = afsd_ni(:) + call icepack_cleanup_fsdn (nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,n)) + + endif ! d_an_newi > puny + endif ! n = 1 + + + ! history/diagnostics + do k = 1, nfsd + ! sum over n + d_afsd_latg(k) = d_afsd_latg(k) & + + area2(n)*afsdn_latg(k,n) & ! after latg + - aicen_init(n)*afsdn(k,n) ! at start + d_afsd_newi(k) = d_afsd_newi(k) & + + aicen(n)*trcrn(nt_fsd+k-1,n) & ! after latg and newi + - area2(n)*afsdn_latg(k,n) ! after latg + enddo ! k + + end subroutine fsd_add_new_ice + +!======================================================================= +! +! Given a wave spectrum, calculate size of new floes based on +! tensile failure, following Shen et al. (2001) +! +! The tensile mode parameter is based on in-situ measurements +! by Roach, Smith & Dean (2018). +! +! authors: Lettie Roach, NIWA/VUW +! + + subroutine wave_dep_growth (nfsd, local_wave_spec, & + wavefreq, dwavefreq, & + new_size) + + integer (kind=int_kind), intent(in) :: & + nfsd ! number of floe size categories + + real (kind=dbl_kind), dimension(:), intent(in) :: & + local_wave_spec ! ocean surface wave spectrum as a function of frequency + ! power spectral density of surface elevation, E(f) (units m^2 s) + ! dimension set in ice_forcing + + + real(kind=dbl_kind), dimension(:), intent(in) :: & + wavefreq, & ! wave frequencies (s^-1) + dwavefreq ! wave frequency bin widths (s^-1) + + integer (kind=int_kind), intent(out) :: & + new_size ! index of floe size category in which new floes will grow + + ! local variables + real (kind=dbl_kind), parameter :: & + tensile_param = 0.167_dbl_kind ! tensile mode parameter + ! value from Roach, Smith & Dean (2018) + ! units kg m^-1 s^-2 + + real (kind=dbl_kind) :: & + mom0, & ! zeroth moment of the spectrum (m) + h_sig, & ! significant wave height (m) + w_amp, & ! wave amplitude (m) + f_peak, & ! peak frequency (s^-1) + w_peak, & ! wavelength from peak freqency (m) + r_max ! floe radius (m) + + integer (kind=int_kind) :: k + + + mom0 = SUM(local_wave_spec*dwavefreq) ! zeroth moment + h_sig = c4*SQRT(mom0) ! sig wave height + w_amp = h_sig/c2 ! sig wave amplitude + f_peak = wavefreq(MAXLOC(local_wave_spec, DIM=1)) ! peak frequency + if (f_peak > puny) w_peak = gravit / (c2*pi*f_peak**c2) ! wavelength from peak freq + + ! tensile failure + if (w_amp > puny) then + r_max = SQRT(c2*tensile_param*w_peak**c2/(pi**c3*w_amp*gravit*rhoi))/c2 + else + r_max = bignum + end if + + new_size = nfsd + do k = nfsd-1, 1, -1 + if (r_max <= floe_rad_h(k)) new_size = k + end do + + end subroutine wave_dep_growth + +!======================================================================= +! +! Floes are perimitted to weld together in freezing conditions, according +! to their geometric probability of overlap if placed randomly on the +! domain. The rate per unit area c_weld is the total number +! of floes that weld with another, per square meter, per unit time, in the +! case of a fully covered ice surface (aice=1), equal to twice the reduction +! in total floe number. See Roach, Smith & Dean (2018). +! +! +! authors: Lettie Roach, NIWA/VUW +! + + subroutine fsd_weld_thermo (ncat, nfsd, & + dt, frzmlt, & + aicen, trcrn, & + d_afsd_weld) + + integer (kind=int_kind), intent(in) :: & + ncat , & ! number of thickness categories + nfsd ! number of floe size categories + + real (kind=dbl_kind), intent(in) :: & + dt ! time step (s) + + real (kind=dbl_kind), dimension (:), intent(in) :: & + aicen ! ice concentration + + real (kind=dbl_kind), intent(in) :: & + frzmlt ! freezing/melting potential (W/m^2) + + real (kind=dbl_kind), dimension (:,:), intent(inout) :: & + trcrn ! ice tracers + + real (kind=dbl_kind), dimension (:), intent(inout) :: & + d_afsd_weld ! change in fsd due to welding + + ! local variables + + real (kind=dbl_kind), parameter :: & + aminweld = p1 ! minimum ice concentration likely to weld + + real (kind=dbl_kind), parameter :: & + c_weld = 1.0e-8_dbl_kind + ! constant of proportionality for welding + ! total number of floes that weld with another, per square meter, + ! per unit time, in the case of a fully covered ice surface + ! units m^-2 s^-1, see documentation for details + + + + integer (kind=int_kind) :: & + nt , & ! time step index + n , & ! thickness category index + k, kx, ky, i, j ! floe size category indices + + real (kind=dbl_kind), dimension(nfsd,ncat) :: & + afsdn ! floe size distribution tracer + + real (kind=dbl_kind), dimension(nfsd,ncat) :: & + d_afsdn_weld ! change in afsdn due to welding + + real (kind=dbl_kind), dimension(nfsd) :: & + stability , & ! check for stability + nfsd_tmp , & ! number fsd + afsd_init , & ! initial values + afsd_tmp , & ! work array + gain, loss ! welding tendencies + + real(kind=dbl_kind) :: & + prefac , & ! multiplies kernel + kern , & ! kernel + subdt , & ! subcycling time step for stability (s) + elapsed_t ! elapsed subcycling time + + + afsdn (:,:) = c0 + afsd_init(:) = c0 + stability = c0 + prefac = p5 + + do n = 1, ncat + + d_afsd_weld (:) = c0 + d_afsdn_weld(:,n) = c0 + afsdn(:,n) = trcrn(nt_fsd:nt_fsd+nfsd-1,n) + call icepack_cleanup_fsdn (nfsd, afsdn(:,n)) + + ! If there is some ice in the lower (nfsd-1) categories + ! and there is freezing potential + if ((frzmlt > puny) .and. & ! freezing potential + (aicen(n) > aminweld) .and. & ! low concentrations area unlikely to weld + (SUM(afsdn(1:nfsd-1,n)) > puny)) then ! some ice in nfsd-1 categories + + afsd_init(:) = afsdn(:,n) ! save initial values + afsd_tmp (:) = afsd_init(:) ! work array + + ! in case of minor numerical errors + WHERE(afsd_tmp < puny) afsd_tmp = c0 + afsd_tmp = afsd_tmp/SUM(afsd_tmp) + + ! adaptive sub-timestep + elapsed_t = c0 + DO WHILE (elapsed_t < dt) + + ! calculate sub timestep + nfsd_tmp = afsd_tmp/floe_area_c + WHERE (afsd_tmp > puny) & + stability = nfsd_tmp/(c_weld*afsd_tmp*aicen(n)) + WHERE (stability < puny) stability = bignum + subdt = MINVAL(stability) + subdt = MIN(subdt,dt) + + loss(:) = c0 + gain(:) = c0 + + do i = 1, nfsd ! consider loss from this category + do j = 1, nfsd ! consider all interaction partners + + k = iweld(i,j) ! product of i+j + + if (k > i) then + + kern = c_weld * floe_area_c(i) * aicen(n) + + loss(i) = loss(i) + kern*afsd_tmp(i)*afsd_tmp(j) + + if (i.eq.j) prefac = c1 ! otherwise 0.5 + + gain(k) = gain(k) + prefac*kern*afsd_tmp(i)*afsd_tmp(j) + + end if + + end do + end do + + ! does largest category lose? +! if (loss(nfsd) > puny) stop 'weld, largest cat losing' +! if (gain(1) > puny) stop 'weld, smallest cat gaining' + + ! update afsd + afsd_tmp(:) = afsd_tmp(:) + subdt*(gain(:) - loss(:)) + + ! in case of minor numerical errors + WHERE(afsd_tmp < puny) afsd_tmp = c0 + afsd_tmp = afsd_tmp/SUM(afsd_tmp) + + ! update time + elapsed_t = elapsed_t + subdt + + ! stop if all in largest floe size cat + if (afsd_tmp(nfsd) > (c1-puny)) exit + + END DO ! time + + call icepack_cleanup_fsdn (nfsd, afsdn(:,n)) + + do k = 1, nfsd + afsdn(k,n) = afsd_tmp(k) + trcrn(nt_fsd+k-1,n) = afsdn(k,n) + ! history/diagnostics + d_afsdn_weld(k,n) = afsdn(k,n) - afsd_init(k) + enddo + + + endif ! try to weld + enddo ! n + + ! history/diagnostics + do k = 1, nfsd + d_afsd_weld(k) = c0 + do n = 1, ncat + d_afsd_weld(k) = d_afsd_weld(k) + aicen(n)*d_afsdn_weld(k,n) + end do ! n + end do ! k + + end subroutine fsd_weld_thermo + +!======================================================================= + + end module icepack_fsd + +!======================================================================= + diff --git a/columnphysics/icepack_intfc.F90 b/columnphysics/icepack_intfc.F90 index 1ed3f8e0c..d93edcdee 100644 --- a/columnphysics/icepack_intfc.F90 +++ b/columnphysics/icepack_intfc.F90 @@ -53,9 +53,16 @@ module icepack_intfc use icepack_itd, only: icepack_init_itd_hist use icepack_itd, only: icepack_aggregate + use icepack_fsd, only: icepack_init_fsd_bounds + use icepack_fsd, only: icepack_init_fsd + use icepack_fsd, only: icepack_cleanup_fsd + use icepack_mechred, only: icepack_step_ridge use icepack_mechred, only: icepack_ice_strength + use icepack_wavefracspec, only: icepack_init_wave + use icepack_wavefracspec, only: icepack_step_wavefracture + use icepack_shortwave, only: icepack_prep_radiation use icepack_shortwave, only: icepack_step_radiation diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90 index abf3838d7..db2d0affc 100644 --- a/columnphysics/icepack_parameters.F90 +++ b/columnphysics/icepack_parameters.F90 @@ -266,6 +266,22 @@ module icepack_parameters ! 2 = WMO standard ! 3 = asymptotic formula +!----------------------------------------------------------------------- +! Parameters for the floe size distribution +!----------------------------------------------------------------------- + + integer (kind=int_kind), public :: & + nfreq = 25 ! number of frequencies + + real (kind=dbl_kind), public :: & + floeshape = 0.666_dbl_kind ! constant from Steele (unitless) + + logical (kind=log_kind), public :: & + wave_spec = .false. ! if true, use wave forcing + + character (len=char_len), public :: & + wave_spec_type = 'constant' ! 'none', 'constant', or 'random' + !----------------------------------------------------------------------- ! Parameters for melt ponds !----------------------------------------------------------------------- @@ -369,6 +385,7 @@ subroutine icepack_init_parameters( & kalg_in, kstrength_in, krdg_partic_in, krdg_redist_in, mu_rdg_in, & atmbndy_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, & tfrz_option_in, kitd_in, kcatbound_in, hs0_in, frzpnd_in, & + floeshape_in, wave_spec_in, wave_spec_type_in, nfreq_in, & dpscale_in, rfracmin_in, rfracmax_in, pndaspect_in, hs1_in, hp1_in, & bgc_flux_type_in, z_tracers_in, scale_bgc_in, solve_zbgc_in, & modal_aero_in, skl_bgc_in, solve_zsal_in, grid_o_in, l_sk_in, & @@ -567,6 +584,22 @@ subroutine icepack_init_parameters( & ! 2 = WMO standard ! 3 = asymptotic formula +!----------------------------------------------------------------------- +! Parameters for the floe size distribution +!----------------------------------------------------------------------- + + integer (kind=int_kind), intent(in), optional :: & + nfreq_in ! number of frequencies + + real (kind=dbl_kind), intent(in), optional :: & + floeshape_in ! constant from Steele (unitless) + + logical (kind=log_kind), intent(in), optional :: & + wave_spec_in ! if true, use wave forcing + + character (len=char_len), intent(in), optional :: & + wave_spec_type_in ! type of wave spectrum forcing + !----------------------------------------------------------------------- ! Parameters for biogeochemistry !----------------------------------------------------------------------- @@ -745,6 +778,10 @@ subroutine icepack_init_parameters( & if (present(tfrz_option_in) ) tfrz_option = tfrz_option_in if (present(kitd_in) ) kitd = kitd_in if (present(kcatbound_in) ) kcatbound = kcatbound_in + if (present(floeshape_in) ) floeshape = floeshape_in + if (present(wave_spec_in) ) wave_spec = wave_spec_in + if (present(wave_spec_type_in) ) wave_spec_type = wave_spec_type_in + if (present(nfreq_in) ) nfreq = nfreq_in if (present(hs0_in) ) hs0 = hs0_in if (present(frzpnd_in) ) frzpnd = frzpnd_in if (present(dpscale_in) ) dpscale = dpscale_in @@ -828,6 +865,7 @@ subroutine icepack_query_parameters( & kalg_out, kstrength_out, krdg_partic_out, krdg_redist_out, mu_rdg_out, & atmbndy_out, calc_strair_out, formdrag_out, highfreq_out, natmiter_out, & tfrz_option_out, kitd_out, kcatbound_out, hs0_out, frzpnd_out, & + floeshape_out, wave_spec_out, wave_spec_type_out, nfreq_out, & dpscale_out, rfracmin_out, rfracmax_out, pndaspect_out, hs1_out, hp1_out, & bgc_flux_type_out, z_tracers_out, scale_bgc_out, solve_zbgc_out, & modal_aero_out, skl_bgc_out, solve_zsal_out, grid_o_out, l_sk_out, & @@ -1035,6 +1073,22 @@ subroutine icepack_query_parameters( & ! 2 = WMO standard ! 3 = asymptotic formula +!----------------------------------------------------------------------- +! Parameters for the floe size distribution +!----------------------------------------------------------------------- + + integer (kind=int_kind), intent(out), optional :: & + nfreq_out ! number of frequencies + + real (kind=dbl_kind), intent(out), optional :: & + floeshape_out ! constant from Steele (unitless) + + logical (kind=log_kind), intent(out), optional :: & + wave_spec_out ! if true, use wave forcing + + character (len=char_len), intent(out), optional :: & + wave_spec_type_out ! type of wave spectrum forcing + !----------------------------------------------------------------------- ! Parameters for biogeochemistry !----------------------------------------------------------------------- @@ -1254,6 +1308,10 @@ subroutine icepack_query_parameters( & if (present(tfrz_option_out) ) tfrz_option_out = tfrz_option if (present(kitd_out) ) kitd_out = kitd if (present(kcatbound_out) ) kcatbound_out = kcatbound + if (present(floeshape_out) ) floeshape_out = floeshape + if (present(wave_spec_out) ) wave_spec_out = wave_spec + if (present(wave_spec_type_out) ) wave_spec_type_out = wave_spec_type + if (present(nfreq_out) ) nfreq_out = nfreq if (present(hs0_out) ) hs0_out = hs0 if (present(frzpnd_out) ) frzpnd_out = frzpnd if (present(dpscale_out) ) dpscale_out = dpscale @@ -1423,6 +1481,10 @@ subroutine icepack_write_parameters(iounit) write(iounit,*) " tfrz_option = ", tfrz_option write(iounit,*) " kitd = ", kitd write(iounit,*) " kcatbound = ", kcatbound + write(iounit,*) " floeshape = ", floeshape + write(iounit,*) " wave_spec = ", wave_spec + write(iounit,*) " wave_spec_type= ", wave_spec_type + write(iounit,*) " nfreq = ", nfreq write(iounit,*) " hs0 = ", hs0 write(iounit,*) " frzpnd = ", frzpnd write(iounit,*) " dpscale = ", dpscale diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index b97f210c9..0f3390fe5 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -31,15 +31,16 @@ module icepack_therm_itd use icepack_tracers, only: ntrcr, nbtrcr use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero - use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY + use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY, nt_fsd use icepack_tracers, only: nt_alvl, nt_vlvl use icepack_tracers, only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo - use icepack_tracers, only: tr_iage, tr_FY, tr_lvl, tr_aero, tr_brine + use icepack_tracers, only: tr_iage, tr_FY, tr_lvl, tr_aero, tr_brine, tr_fsd use icepack_tracers, only: bio_index 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_itd, only: reduce_area, cleanup_itd use icepack_itd, only: aggregate_area, shift_ice use icepack_itd, only: column_sum, column_conservation_check @@ -57,7 +58,7 @@ module icepack_therm_itd lateral_melt, & icepack_step_therm2 - logical (kind=log_kind), parameter :: & + logical (kind=log_kind), parameter, public :: & l_conservation_check = .false. ! if true, check conservation ! (useful for debugging) @@ -847,6 +848,7 @@ end subroutine update_vertical_tracers ! ! author: C. M. Bitz, UW ! 2003: Modified by William H. Lipscomb and Elizabeth C. Hunke, LANL +! 2016 Lettie Roach, NIWA/VUW, added floe size dependence ! subroutine lateral_melt (dt, ncat, & nilyr, nslyr, & @@ -854,16 +856,20 @@ subroutine lateral_melt (dt, ncat, & fresh, fsalt, & fhocn, faero_ocn, & rside, meltl, & + fside, sss, & aicen, vicen, & vsnon, trcrn, & fzsal, flux_bio, & - nbtrcr, nblyr) + nbtrcr, nblyr, & + nfsd, d_afsd_latm,& + floe_rad_c, floe_binwidth) real (kind=dbl_kind), intent(in) :: & dt ! time step (s) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories + nfsd , & ! number of floe size categories nilyr , & ! number of ice layers nblyr , & ! number of bio layers nslyr , & ! number of snow layers @@ -875,11 +881,12 @@ subroutine lateral_melt (dt, ncat, & vicen , & ! volume per unit area of ice (m) vsnon ! volume per unit area of snow (m) - real (kind=dbl_kind), dimension (:,:), intent(in) :: & - trcrn ! tracer array + real (kind=dbl_kind), dimension (:,:), intent(inout) :: & + trcrn ! tracer array real (kind=dbl_kind), intent(in) :: & - rside ! fraction of ice that melts laterally + rside , & ! fraction of ice that melts laterally + fside ! lateral heat flux (W/m^2) real (kind=dbl_kind), intent(inout) :: & fpond , & ! fresh water flux to ponds (kg/m^2/s) @@ -889,6 +896,13 @@ subroutine lateral_melt (dt, ncat, & meltl , & ! lateral ice melt (m/step-->cm/day) fzsal ! salt flux from zsalinity (kg/m2/s) + real (kind=dbl_kind), dimension (:), intent(in) :: & + 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) :: & + d_afsd_latm ! change in fsd due to lateral melt (m) + real (kind=dbl_kind), dimension(nbtrcr), & intent(inout) :: & flux_bio ! biology tracer flux from layer bgc (mmol/m^2/s) @@ -908,16 +922,142 @@ subroutine lateral_melt (dt, ncat, & dfresh , & ! change in fresh dfsalt , & ! change in fsalt dvssl , & ! snow surface layer volume - dvint ! snow interior layer + dvint , & ! snow interior layer + cat1_arealoss, tmp ! + + logical (kind=log_kind) :: & + flag ! .true. if there could be lateral melting real (kind=dbl_kind), dimension (ncat) :: & - vicen_init ! volume per unit area of ice (m) + aicen_init, & ! initial area fraction + vicen_init, & ! volume per unit area of ice (m) + G_radialn , & ! rate of lateral melt (m/s) + delta_an , & ! change in the ITD + qin , & ! enthalpy + rsiden ! delta_an/aicen + + real (kind=dbl_kind), dimension (nfsd,ncat) :: & + afsdn , & ! floe size distribution tracer + afsdn_init ! initial value - character(len=*),parameter :: subname='(lateral_melt)' + real (kind=dbl_kind), dimension (nfsd) :: & + df_flx ! finite difference for G_r * areal mFSTD tilda - if (rside > c0) then ! grid cells with lateral melting. + real (kind=dbl_kind), dimension(nfsd+1) :: & + f_flx ! + +!echmod - for average qin + real (kind=dbl_kind), intent(in) :: & + sss + real (kind=dbl_kind) :: & + Ti, Si0, qi0 + + character(len=*), parameter :: subname='(lateral_melt)' + + flag = .false. + dfhocn = c0 + dfpond = c0 + dfresh = c0 + dfsalt = c0 + dvssl = c0 + dvint = c0 + cat1_arealoss = c0 + tmp = c0 + vicen_init = c0 + G_radialn = c0 + delta_an = c0 + qin = c0 + rsiden = c0 + afsdn = c1 + afsdn_init = c0 + df_flx = c0 + f_flx = c0 + + if (tr_fsd) then + call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:)) + afsdn = trcrn(nt_fsd:nt_fsd+nfsd-1,:) + aicen_init = aicen + afsdn_init = afsdn ! for diagnostics + d_afsd_latm(:) = c0 + end if + + if (tr_fsd .and. fside < c0) then + flag = .true. + + +! echmod - using category values would be preferable to the average value + ! Compute average enthalpy of ice (taken from add_new_ice) + if (sss > c2 * dSin0_frazil) then + Si0 = sss - dSin0_frazil + else + Si0 = sss**2 / (c4*dSin0_frazil) + endif + Ti = min(liquidus_temperature_mush(Si0/phi_init), -p1) + qi0 = enthalpy_mush(Ti, Si0) do n = 1, ncat + if (ktherm == 2) then ! mushy + do k = 1, nilyr + !qin(n) = qin(n) & + ! + trcrn(nt_qice+k-1,n)*vicen(n)/real(nilyr,kind=dbl_kind) + qin(n) = qi0 + enddo + else + qin(n) = -rhoi*Lfresh + endif + + 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 + + + if (any(afsdn(:,n) < c0)) print*,& + 'lateral_melt B afsd < 0',n + + cat1_arealoss = -trcrn(nt_fsd+1-1,n) * aicen(n) * dt & + * G_radialn(n) / floe_binwidth(1) + + delta_an(n) = c0 + do k = 1, nfsd + delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k))*aicen(n) & + * trcrn(nt_fsd+k-1,n)*G_radialn(n)*dt) ! delta_an < 0 + end do + + ! add negative area loss from fsd + delta_an(n) = delta_an(n) - cat1_arealoss + + 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 (rsiden(n) < c0) print*,'ERROR rsiden < 0', rsiden(n) + + end if ! G_radialn + enddo ! ncat + + else if (rside > c0) then ! original, non-fsd implementation + + flag = .true. + rsiden(:) = rside ! initialize + + endif + + 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. @@ -926,39 +1066,80 @@ subroutine lateral_melt (dt, ncat, & ! fluxes to coupler ! dfresh > 0, dfsalt > 0, dfpond > 0 - dfresh = (rhos*vsnon(n) + rhoi*vicen(n)) * rside / dt - dfsalt = rhoi*vicen(n)*ice_ref_salinity*p001 * rside / dt + dfresh = (rhoi*vicen(n) + rhos*vsnon(n)) * rsiden(n) / dt + dfsalt = rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt fresh = fresh + dfresh fsalt = fsalt + dfsalt if (tr_pond_topo) then - dfpond = aicen(n) & - * trcrn(nt_apnd,n) & - * trcrn(nt_hpnd,n) & - * rside - fpond = fpond - dfpond + dfpond = aicen(n)*trcrn(nt_apnd,n)*trcrn(nt_hpnd,n)*rsiden(n) + fpond = fpond - dfpond endif ! history diagnostics - meltl = meltl + vicen(n)*rside + meltl = meltl + vicen(n)*rsiden(n) ! state variables vicen_init(n) = vicen(n) - aicen(n) = aicen(n) * (c1 - rside) - vicen(n) = vicen(n) * (c1 - rside) - vsnon(n) = vsnon(n) * (c1 - rside) + aicen(n) = aicen(n) * (c1 - rsiden(n)) + vicen(n) = vicen(n) * (c1 - rsiden(n)) + vsnon(n) = vsnon(n) * (c1 - rsiden(n)) + + ! floe size distribution + 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 + + end if ! aicen + end if ! rside > 0, otherwise do nothing + + end if ! tr_fsd + + ! fluxes do k = 1, nilyr ! enthalpy tracers do not change (e/v constant) ! heat flux to coupler for ice melt (dfhocn < 0) - dfhocn = trcrn(nt_qice+k-1,n)*rside / dt & + dfhocn = trcrn(nt_qice+k-1,n)*rsiden(n) / dt & * vicen(n)/real(nilyr,kind=dbl_kind) fhocn = fhocn + dfhocn enddo ! nilyr do k = 1, nslyr ! heat flux to coupler for snow melt (dfhocn < 0) - dfhocn = trcrn(nt_qsno+k-1,n)*rside / dt & + dfhocn = trcrn(nt_qsno+k-1,n)*rsiden(n) / dt & * vsnon(n)/real(nslyr,kind=dbl_kind) fhocn = fhocn + dfhocn enddo ! nslyr @@ -971,7 +1152,7 @@ subroutine lateral_melt (dt, ncat, & + vicen(n) & *(trcrn(nt_aero+2+4*(k-1),n) & + trcrn(nt_aero+3+4*(k-1),n))) & - * rside / dt + * rsiden(n) / dt enddo ! k endif ! tr_aero @@ -986,7 +1167,7 @@ subroutine lateral_melt (dt, ncat, & flux_bio(k) = flux_bio(k) & + (trcrn(bio_index(k)+nblyr+1,n)*dvssl & + trcrn(bio_index(k)+nblyr+2,n)*dvint) & - * rside / dt + * rsiden(n) / dt enddo endif @@ -995,12 +1176,25 @@ subroutine lateral_melt (dt, ncat, & if (solve_zsal .or. z_tracers) & call lateral_melt_bgc(dt, & ncat, nblyr, & - rside, vicen_init, & + rside, vicen_init, & !echmod: use rsiden trcrn, fzsal, & flux_bio, nbtrcr) if (icepack_warnings_aborted(subname)) return - endif ! rside + endif ! flag + + if (tr_fsd) then + !call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) ) + + ! diagnostics + do k = 1, nfsd + d_afsd_latm(k) = c0 + do n = 1, ncat + d_afsd_latm(k) = d_afsd_latm(k) & + + afsdn(k,n)*aicen(n) - afsdn_init(k,n)*aicen_init(n) + end do + end do + end if end subroutine lateral_melt @@ -1020,13 +1214,15 @@ end subroutine lateral_melt ! added to only category 1, all formulations combine the new ice and ! existing ice tracers as bulk quantities. ! -! authors William H. Lipscomb, LANL -! Elizabeth C. Hunke, LANL -! Adrian Turner, LANL +! authors: William H. Lipscomb, LANL +! Elizabeth C. Hunke, LANL +! Adrian Turner, LANL +! Lettie Roach, NIWA/VUW ! - subroutine add_new_ice (ncat, nilyr, nblyr, & + subroutine add_new_ice (ncat, nilyr, & + nfsd, nblyr, & n_aero, dt, & - ntrcr, nltrcr, & + ntrcr, nltrcr, & hin_max, ktherm, & aicen, trcrn, & vicen, vsnon1, & @@ -1041,10 +1237,20 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & bgrid, cgrid, igrid, & nbtrcr, flux_bio, & ocean_bio, fzsal, & - frazil_diag ) + frazil_diag, & + wave_sig_ht, & + wave_spectrum, & + wavefreq, & + dwavefreq, & + d_afsd_latg, & + d_afsd_newi, & + floe_rad_c, floe_binwidth) + + use icepack_fsd, only: fsd_lateral_growth, fsd_add_new_ice integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories + nfsd , & ! number of floe size categories nilyr , & ! number of ice layers nblyr , & ! number of bio layers ntrcr , & ! number of tracers @@ -1114,12 +1320,38 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & ocean_bio ! ocean concentration of biological tracer ! zsalinity - real (kind=dbl_kind), intent(inout) :: & + real (kind=dbl_kind), intent(inout) :: & fzsal ! salt flux to ocean from zsalinity (kg/m^2/s) + ! floe size distribution + real (kind=dbl_kind), intent(in) :: & + wave_sig_ht ! significant height of waves globally (m) + + real (kind=dbl_kind), dimension(:), intent(in) :: & + wave_spectrum ! ocean surface wave spectrum, E(f) (m^2 s) + + real(kind=dbl_kind), dimension(:), intent(in) :: & + wavefreq, & ! wave frequencies (s^-1) + dwavefreq ! wave frequency bin widths (s^-1) + + real (kind=dbl_kind), dimension (:), intent(in) :: & + floe_rad_c , & ! fsd size bin centre in m (radius) + floe_binwidth ! fsd size bin width in m (radius) + + real (kind=dbl_kind), dimension(ncat) :: & ! for now + ! change in thickness distribution (area) + d_an_latg , & ! due to fsd lateral growth + d_an_newi ! new ice formation + + real (kind=dbl_kind), dimension(:), intent(out) :: & + ! change in thickness distribution (area) + d_afsd_latg , & ! due to fsd lateral growth + d_afsd_newi ! new ice formation + ! local variables integer (kind=int_kind) :: & + ncats , & ! max categories to which new ice is added, initially n , & ! ice category index k , & ! ice layer index it ! aerosol tracer index @@ -1127,6 +1359,7 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & real (kind=dbl_kind) :: & ai0new , & ! area of new ice added to cat 1 vi0new , & ! volume of new ice added to cat 1 + vi0new_lat , & ! volume of new ice added laterally to fsd hsurp , & ! thickness of new ice added to each cat fnew , & ! heat flx to open water for new ice (W/m^2) hi0new , & ! thickness of new ice @@ -1157,23 +1390,76 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & aicen_init, & ! fractional area of ice vicen_init ! volume per unit area of ice (m) + ! floe size distribution + real (kind=dbl_kind), dimension (nfsd,ncat) :: & + afsdn ! floe size distribution tracer (originally areal_mfstd_init) + +! real (kind=dbl_kind), dimension (nfsd) :: & +! afsd , & ! fsd tracer for each thickness category + + real (kind=dbl_kind), dimension (ncat) :: & + d_an_tot, & ! change in the ITD due to lateral growth and new ice + area2 ! area after lateral growth and before new ice formation + + real (kind=dbl_kind), dimension (ncat) :: & + vin0new ! volume of new ice added to any thickness cat + + real (kind=dbl_kind), dimension (nfsd) :: & + afsd_ni ! areal mFSTD after new ice added + + real (kind=dbl_kind) :: & + tmp, & + latsurf_area, & ! fractional area of ice on sides of floes + lead_area , & ! fractional area of ice in lead region + G_radial , & ! lateral melt rate (m/s) + tot_latg , & ! total fsd lateral growth in open water + ai0mod ! ai0new - tot_latg + character(len=*),parameter :: subname='(add_new_ice)' !----------------------------------------------------------------- ! initialize !----------------------------------------------------------------- + hsurp = c0 + hi0new = c0 + ai0new = c0 + afsdn(:,:) = c0 + d_an_latg(:) = c0 + d_an_tot(:) = c0 + d_an_newi(:) = c0 + + if (tr_fsd) then + d_afsd_latg(:) = c0 ! diagnostics + d_afsd_newi(:) = c0 + end if + + area2(:) = aicen(:) + lead_area = c0 + latsurf_area = c0 + G_radial = c0 + tot_latg = c0 if (ncat > 1) then hi0max = hin_max(1)*0.9_dbl_kind ! not too close to boundary else hi0max = bignum ! big number endif + if (tr_fsd) then + call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:)) + endif + do n = 1, ncat aicen_init(n) = aicen(n) vicen_init(n) = vicen(n) eicen(n) = c0 + if (tr_fsd) then + do k = 1, nfsd + afsdn(k,n) = trcrn(nt_fsd+k-1,n) + enddo + endif enddo +! if (any(afsdn < c0)) print*,'add_new B afsdn < 0' if (l_conservation_check) then @@ -1232,9 +1518,7 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & ! history diagnostics frazil = vi0new - if (solve_zsal) & - fzsal = fzsal & - - rhosi*vi0new/dt*p001*sss*salt_loss + if (solve_zsal) fzsal = fzsal - rhosi*vi0new/dt*p001*sss*salt_loss if (present(frz_onset) .and. present(yday)) then if (frazil > puny .and. frz_onset < puny) frz_onset = yday @@ -1251,8 +1535,8 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & if (ktherm <= 1) then dfresh = -rhoi*vi0new/dt dfsalt = ice_ref_salinity*p001*dfresh - fresh = fresh + dfresh - fsalt = fsalt + dfsalt + fresh = fresh + dfresh + fsalt = fsalt + dfsalt ! elseif (ktherm == 2) the fluxes are added elsewhere endif else ! update_ocn_f = false @@ -1271,20 +1555,33 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & ! Decide how to distribute the new ice. !----------------------------------------------------------------- - hsurp = c0 - ai0new = c0 - if (vi0new > c0) then + if (tr_fsd) & ! lateral growth of existing ice + ! calculate change in conc due to lateral growth + ! update vi0new, without change to afsdn or aicen + call fsd_lateral_growth (ncat, nfsd, & + dt, aice, & + aicen, vicen, & + vi0new, frazil, & + floe_rad_c, afsdn, & + lead_area, latsurf_area, & + G_radial, d_an_latg, & + tot_latg) + + ai0mod = aice0 + ! separate frazil ice growth from lateral ice growth + if (tr_fsd) ai0mod = aice0-tot_latg + ! new ice area and thickness ! hin_max(0) < new ice thickness < hin_max(1) - if (aice0 > puny) then - hi0new = max(vi0new/aice0, hfrazilmin) - if (hi0new > hi0max .and. aice0+puny < c1) then + if (ai0mod > puny) then + hi0new = max(vi0new/ai0mod, hfrazilmin) + if (hi0new > hi0max .and. ai0mod+puny < c1) then ! distribute excess volume over all categories (below) hi0new = hi0max - ai0new = aice0 - vsurp = vi0new - ai0new*hi0new + ai0new = ai0mod + vsurp = vi0new - ai0new*hi0new hsurp = vsurp / aice vi0new = ai0new*hi0new else @@ -1295,6 +1592,23 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & hsurp = vi0new/aice ! new thickness in each cat vi0new = c0 endif ! aice0 > puny + + ! volume added to each category from lateral growth + vin0new(:) = c0 + do n = 1, ncat + if (aicen(n) > c0) vin0new(n) = d_an_latg(n) * vicen(n)/aicen(n) + end do + + ! combine areal change from new ice growth and lateral growth + d_an_newi(1) = ai0new + d_an_tot(2:ncat) = d_an_latg(2:ncat) + d_an_tot(1) = d_an_latg(1) + d_an_newi(1) + if (tr_fsd) then + vin0new(1) = vin0new(1) + ai0new*hi0new ! not BFB + else + vin0new(1) = vi0new + endif + endif ! vi0new > puny !----------------------------------------------------------------- @@ -1375,79 +1689,103 @@ subroutine add_new_ice (ncat, nilyr, nblyr, & endif ! hsurp > 0 !----------------------------------------------------------------- - ! Combine new ice grown in open water with category 1 ice. + ! Combine new ice grown in open water with ice categories. + ! Using the floe size distribution, ice is added laterally to all + ! categories; otherwise it is added to category 1. ! Assume that vsnon and esnon are unchanged. ! The mushy formulation assumes salt from frazil is added uniformly ! to category 1, while the others use a salinity profile. !----------------------------------------------------------------- - if (vi0new > c0) then ! add ice to category 1 + ncats = 1 ! add new ice to category 1 by default + if (tr_fsd) ncats = ncat ! add new ice laterally to all categories + - area1 = aicen(1) ! save - vice1 = vicen(1) ! save - aicen(1) = aicen(1) + ai0new - aice0 = aice0 - ai0new - vicen(1) = vicen(1) + vi0new + do n = 1, ncats - trcrn(nt_Tsfc,1) = & - (trcrn(nt_Tsfc,1)*area1 + Tf*ai0new)/aicen(1) - trcrn(nt_Tsfc,1) = min (trcrn(nt_Tsfc,1), c0) + if (d_an_tot(n) > c0 .and. vin0new(n) > c0) then ! add ice to category n + + area1 = aicen(n) ! save + vice1 = vicen(n) ! save + area2(n) = aicen_init(n) + d_an_latg(n) ! save area after latg, before newi + aicen(n) = aicen(n) + d_an_tot(n) ! after lateral growth and new ice growth + + aice0 = aice0 - d_an_tot(n) + vicen(n) = vicen(n) + vin0new(n) + + trcrn(nt_Tsfc,n) = (trcrn(nt_Tsfc,n)*area1 + Tf*d_an_tot(n))/aicen(n) + trcrn(nt_Tsfc,n) = min (trcrn(nt_Tsfc,n), c0) if (tr_FY) then - trcrn(nt_FY,1) = & - (trcrn(nt_FY,1)*area1 + ai0new)/aicen(1) - trcrn(nt_FY,1) = min(trcrn(nt_FY,1), c1) + trcrn(nt_FY,n) = (trcrn(nt_FY,n)*area1 + d_an_tot(n))/aicen(n) + trcrn(nt_FY,n) = min(trcrn(nt_FY,n), c1) endif - if (vicen(1) > puny) then + if (tr_fsd) & ! evolve the floe size distribution + ! both new frazil ice and lateral growth + call fsd_add_new_ice (ncat, n, nfsd, & + dt, ai0new, & + d_an_latg, d_an_newi, & + floe_rad_c, floe_binwidth, & + G_radial, area2, & + wave_sig_ht, & + wave_spectrum, & + wavefreq, & + dwavefreq, & + d_afsd_latg, & + d_afsd_newi, & + afsdn, aicen_init, & + aicen, trcrn) + + if (vicen(n) > puny) then if (tr_iage) & - trcrn(nt_iage,1) = & - (trcrn(nt_iage,1)*vice1 + dt*vi0new)/vicen(1) + trcrn(nt_iage,n) = (trcrn(nt_iage,n)*vice1 + dt*vin0new(n))/vicen(n) if (tr_aero) then do it = 1, n_aero - trcrn(nt_aero+2+4*(it-1),1) = & - trcrn(nt_aero+2+4*(it-1),1)*vice1/vicen(1) - trcrn(nt_aero+3+4*(it-1),1) = & - trcrn(nt_aero+3+4*(it-1),1)*vice1/vicen(1) + trcrn(nt_aero+2+4*(it-1),n) = & + trcrn(nt_aero+2+4*(it-1),n)*vice1/vicen(n) + trcrn(nt_aero+3+4*(it-1),n) = & + trcrn(nt_aero+3+4*(it-1),n)*vice1/vicen(n) enddo endif if (tr_lvl) then - alvl = trcrn(nt_alvl,1) - trcrn(nt_alvl,1) = & - (trcrn(nt_alvl,1)*area1 + ai0new)/aicen(1) - trcrn(nt_vlvl,1) = & - (trcrn(nt_vlvl,1)*vice1 + vi0new)/vicen(1) + alvl = trcrn(nt_alvl,n) + trcrn(nt_alvl,n) = & + (trcrn(nt_alvl,n)*area1 + d_an_tot(n))/aicen(n) + trcrn(nt_vlvl,n) = & + (trcrn(nt_vlvl,n)*vice1 + vin0new(n))/vicen(n) endif if (tr_pond_cesm .or. tr_pond_topo) then - trcrn(nt_apnd,1) = & - trcrn(nt_apnd,1)*area1/aicen(1) + trcrn(nt_apnd,n) = & + trcrn(nt_apnd,n)*area1/aicen(n) elseif (tr_pond_lvl) then - if (trcrn(nt_alvl,1) > puny) then - trcrn(nt_apnd,1) = & - trcrn(nt_apnd,1) * alvl*area1 / (trcrn(nt_alvl,1)*aicen(1)) + if (trcrn(nt_alvl,n) > puny) then + trcrn(nt_apnd,n) = & + trcrn(nt_apnd,n) * alvl*area1 / (trcrn(nt_alvl,n)*aicen(n)) endif endif endif do k = 1, nilyr - - if (vicen(1) > c0) then + if (vicen(n) > c0) then ! factor of nilyr cancels out ! enthalpy - trcrn(nt_qice+k-1,1) = & - (trcrn(nt_qice+k-1,1)*vice1 + qi0new*vi0new)/vicen(1) + trcrn(nt_qice+k-1,n) = & + (trcrn(nt_qice+k-1,n)*vice1 + qi0new*vin0new(n))/vicen(n) ! salinity if (.NOT. solve_zsal)& - trcrn(nt_sice+k-1,1) = & - (trcrn(nt_sice+k-1,1)*vice1 + Sprofile(k)*vi0new)/vicen(1) + trcrn(nt_sice+k-1,n) = & + (trcrn(nt_sice+k-1,n)*vice1 + Sprofile(k)*vin0new(n))/vicen(n) endif enddo endif ! vi0new > 0 + enddo ! ncats + if (l_conservation_check) then do n = 1, ncat @@ -1514,6 +1852,7 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, nltrcr, & Tf, sss, & salinz, & rside, meltl, & + fside, & frzmlt, frazil, & frain, fpond, & fresh, fsalt, & @@ -1523,10 +1862,18 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, nltrcr, & first_ice, fzsal, & flux_bio, ocean_bio, & frazil_diag, & - frz_onset, yday) + frz_onset, yday, & + nfsd, wave_sig_ht, & + wave_spectrum, & + wavefreq, & + dwavefreq, & + d_afsd_latg, d_afsd_newi, & + d_afsd_latm, d_afsd_weld, & + floe_rad_c, floe_binwidth) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories + nfsd , & ! number of floe size categories nltrcr , & ! number of zbgc tracers nblyr , & ! number of bio layers nilyr , & ! number of ice layers @@ -1544,7 +1891,20 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, nltrcr, & Tf , & ! freezing temperature (C) sss , & ! sea surface salinity (ppt) rside , & ! fraction of ice that melts laterally - frzmlt ! freezing/melting potential (W/m^2) + fside , & ! lateral heat flux (W/m^2) + frzmlt , & ! freezing/melting potential (W/m^2) + wave_sig_ht ! significant height of waves in ice (m) + + real (kind=dbl_kind), dimension(:), intent(in) :: & + wave_spectrum ! ocean surface wave spectrum E(f) (m^2 s) + + real(kind=dbl_kind), dimension(:), intent(in) :: & + wavefreq, & ! wave frequencies (s^-1) + dwavefreq ! wave frequency bin widths (s^-1) + + real (kind=dbl_kind), dimension (:), intent(in) :: & + floe_rad_c , & ! fsd size bin centre in m (radius) + floe_binwidth ! fsd size bin width in m (radius) integer (kind=int_kind), dimension (:), intent(in) :: & trcr_depend, & ! = 0 for aicen tracers, 1 for vicen, 2 for vsnon @@ -1598,6 +1958,13 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, nltrcr, & logical (kind=log_kind), dimension(:), intent(inout) :: & first_ice ! true until ice forms + real (kind=dbl_kind), dimension(:), intent(out) :: & + ! change in floe size distribution (area) + d_afsd_latg , & ! due to fsd lateral growth + d_afsd_newi , & ! new ice formation + d_afsd_latm , & ! lateral melt + d_afsd_weld ! welding + real (kind=dbl_kind), intent(inout), optional :: & frz_onset ! day of year that freezing begins (congel or frazil) @@ -1608,6 +1975,7 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, nltrcr, & character(len=*),parameter :: subname='(icepack_step_therm2)' + !----------------------------------------------------------------- ! Let rain drain through to the ocean. !----------------------------------------------------------------- @@ -1666,7 +2034,7 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, nltrcr, & ! identify ice-ocean cells call add_new_ice (ncat, nilyr, & - nblyr, & + nfsd, nblyr, & n_aero, dt, & ntrcr, nltrcr, & hin_max, ktherm, & @@ -1683,7 +2051,12 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, nltrcr, & cgrid, igrid, & nbtrcr, flux_bio, & ocean_bio, fzsal, & - frazil_diag ) + frazil_diag, & + wave_sig_ht, & + wave_spectrum, & + wavefreq, dwavefreq, & + d_afsd_latg, d_afsd_newi, & + floe_rad_c, floe_binwidth) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- @@ -1696,10 +2069,21 @@ subroutine icepack_step_therm2 (dt, ncat, n_aero, nltrcr, & fresh, fsalt, & fhocn, faero_ocn, & rside, meltl, & + fside, sss, & aicen, vicen, & vsnon, trcrn, & fzsal, flux_bio, & - nbtrcr, nblyr) + nbtrcr, nblyr, & + nfsd, d_afsd_latm, & + floe_rad_c,floe_binwidth) + if (icepack_warnings_aborted(subname)) return + + ! Floe welding during freezing conditions + if (tr_fsd) & + call fsd_weld_thermo (ncat, nfsd, & + dt, frzmlt, & + aicen, trcrn, & + d_afsd_weld) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index f976c7c17..9645f436c 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -28,7 +28,7 @@ module icepack_therm_vertical use icepack_parameters, only: rfracmin, rfracmax, pndaspect, dpscale, frzpnd use icepack_parameters, only: phi_i_mushy - use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond + use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd use icepack_tracers, only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo use icepack_therm_shared, only: ferrmax, l_brine @@ -492,7 +492,8 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & fbot_xfer_type, & strocnxT, strocnyT, & Tbot, fbot, & - rside, Cdn_ocn) + rside, Cdn_ocn, & + fside) integer (kind=int_kind), intent(in) :: & ncat , & ! number of thickness categories @@ -526,7 +527,8 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & real (kind=dbl_kind), intent(out) :: & Tbot , & ! ice bottom surface temperature (deg C) fbot , & ! heat flux to ice bottom (W/m^2) - rside ! fraction of ice that melts laterally + rside , & ! fraction of ice that melts laterally + fside ! lateral heat flux (W/m^2) ! local variables @@ -536,7 +538,7 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & real (kind=dbl_kind) :: & etot , & ! total energy in column - fside ! lateral heat flux (W/m^2) + qavg ! average enthalpy in column (approximate) real (kind=dbl_kind) :: & deltaT , & ! SST - Tbot >= 0 @@ -566,13 +568,13 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & !----------------------------------------------------------------- rside = c0 + fside = c0 Tbot = Tf fbot = c0 - + wlat = c0 + if (aice > puny .and. frzmlt < c0) then ! ice can melt - fside = c0 - !----------------------------------------------------------------- ! Use boundary layer theory for fbot. ! See Maykut and McPhee (1995): JGR, 100, 24,691-24,703. @@ -616,19 +618,26 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & do n = 1, ncat etot = c0 + qavg = c0 ! melting energy/unit area in each column, etot < 0 do k = 1, nslyr etot = etot + qsnon(k,n) * vsnon(n)/real(nslyr,kind=dbl_kind) + qavg = qavg + qsnon(k,n) enddo do k = 1, nilyr etot = etot + qicen(k,n) * vicen(n)/real(nilyr,kind=dbl_kind) + qavg = qavg + qicen(k,n) enddo ! nilyr - ! lateral heat flux - fside = fside + rside*etot/dt ! fside < 0 + ! lateral heat flux, fside < 0 + if (tr_fsd) then ! floe size distribution + fside = fside + wlat*qavg + else ! default floe size + fside = fside + rside*etot/dt + endif enddo ! n @@ -640,6 +649,7 @@ subroutine frzmlt_bottom_lateral (dt, ncat, & xtmp = min(xtmp, c1) fbot = fbot * xtmp rside = rside * xtmp + fside = fside * xtmp endif @@ -2050,8 +2060,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, & sss , Tf , & strocnxT , strocnyT , & fbot , & - Tbot , Tsnice , & + Tbot , Tsnice , & frzmlt , rside , & + fside , & fsnow , frain , & fpond , & fsurf , fsurfn , & @@ -2163,10 +2174,11 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, & fbot , & ! ice-ocean heat flux at bottom surface (W/m^2) frzmlt , & ! freezing/melting potential (W/m^2) rside , & ! fraction of ice that melts laterally + fside , & ! lateral heat flux (W/m^2) sst , & ! sea surface temperature (C) Tf , & ! freezing temperature (C) Tbot , & ! ice bottom surface temperature (deg C) - Tsnice , & ! snow ice interface temperature (deg C) + Tsnice , & ! snow ice interface temperature (deg C) sss , & ! sea surface salinity (ppt) meltt , & ! top ice melt (m/step-->cm/day) melts , & ! snow melt (m/step-->cm/day) @@ -2274,7 +2286,9 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr, n_aero, & fbot_xfer_type, & strocnxT, strocnyT, & Tbot, fbot, & - rside, Cdn_ocn) + rside, Cdn_ocn, & + fside) + if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- diff --git a/columnphysics/icepack_tracers.F90 b/columnphysics/icepack_tracers.F90 index 3b2ed529e..71e9b3a6b 100644 --- a/columnphysics/icepack_tracers.F90 +++ b/columnphysics/icepack_tracers.F90 @@ -73,6 +73,7 @@ module icepack_tracers nt_apnd , & ! melt pond area fraction nt_hpnd , & ! melt pond depth nt_ipnd , & ! melt pond refrozen lid thickness + nt_fsd , & ! floe size distribution nt_aero , & ! starting index for aerosols in ice nt_bgc_Nit, & ! nutrients nt_bgc_Am, & ! @@ -94,7 +95,8 @@ module icepack_tracers tr_pond_lvl , & ! if .true., use level-ice pond tracer tr_pond_topo, & ! if .true., use explicit topography-based ponds tr_aero , & ! if .true., use aerosol tracers - tr_brine ! if .true., brine height differs from ice thickness + tr_brine , & ! if .true., brine height differs from ice thickness + tr_fsd ! if .true., use floe size distribution !----------------------------------------------------------------- ! biogeochemistry @@ -249,7 +251,7 @@ end subroutine icepack_write_tracer_sizes subroutine icepack_init_tracer_flags(& tr_iage_in, tr_FY_in, tr_lvl_in, & tr_pond_in, tr_pond_cesm_in, tr_pond_lvl_in, tr_pond_topo_in, & - tr_aero_in, tr_brine_in, tr_zaero_in, & + tr_fsd_in, tr_aero_in, tr_brine_in, tr_zaero_in, & tr_bgc_Nit_in, tr_bgc_N_in, tr_bgc_DON_in, tr_bgc_C_in, tr_bgc_chl_in, & tr_bgc_Am_in, tr_bgc_Sil_in, tr_bgc_DMS_in, tr_bgc_Fe_in, tr_bgc_hum_in, & tr_bgc_PON_in) @@ -262,6 +264,7 @@ subroutine icepack_init_tracer_flags(& tr_pond_cesm_in , & ! if .true., use cesm pond tracer tr_pond_lvl_in , & ! if .true., use level-ice pond tracer tr_pond_topo_in , & ! if .true., use explicit topography-based ponds + tr_fsd_in , & ! if .true., use floe size distribution tracers tr_aero_in , & ! if .true., use aerosol tracers tr_brine_in , & ! if .true., brine height differs from ice thickness tr_zaero_in , & ! if .true., black carbon is tracers (n_zaero) @@ -286,6 +289,7 @@ subroutine icepack_init_tracer_flags(& if (present(tr_pond_cesm_in)) tr_pond_cesm = tr_pond_cesm_in if (present(tr_pond_lvl_in) ) tr_pond_lvl = tr_pond_lvl_in if (present(tr_pond_topo_in)) tr_pond_topo = tr_pond_topo_in + if (present(tr_fsd_in) ) tr_fsd = tr_fsd_in if (present(tr_aero_in) ) tr_aero = tr_aero_in if (present(tr_brine_in) ) tr_brine = tr_brine_in if (present(tr_zaero_in) ) tr_zaero = tr_zaero_in @@ -310,7 +314,7 @@ end subroutine icepack_init_tracer_flags subroutine icepack_query_tracer_flags(& tr_iage_out, tr_FY_out, tr_lvl_out, & tr_pond_out, tr_pond_cesm_out, tr_pond_lvl_out, tr_pond_topo_out, & - tr_aero_out, tr_brine_out, tr_zaero_out, & + tr_fsd_out, tr_aero_out, tr_brine_out, tr_zaero_out, & tr_bgc_Nit_out, tr_bgc_N_out, tr_bgc_DON_out, tr_bgc_C_out, tr_bgc_chl_out, & tr_bgc_Am_out, tr_bgc_Sil_out, tr_bgc_DMS_out, tr_bgc_Fe_out, tr_bgc_hum_out, & tr_bgc_PON_out) @@ -323,6 +327,7 @@ subroutine icepack_query_tracer_flags(& tr_pond_cesm_out , & ! if .true., use cesm pond tracer tr_pond_lvl_out , & ! if .true., use level-ice pond tracer tr_pond_topo_out , & ! if .true., use explicit topography-based ponds + tr_fsd_out , & ! if .true., use floe size distribution tr_aero_out , & ! if .true., use aerosol tracers tr_brine_out , & ! if .true., brine height differs from ice thickness tr_zaero_out , & ! if .true., black carbon is tracers (n_zaero) @@ -347,6 +352,7 @@ subroutine icepack_query_tracer_flags(& if (present(tr_pond_cesm_out)) tr_pond_cesm_out = tr_pond_cesm if (present(tr_pond_lvl_out) ) tr_pond_lvl_out = tr_pond_lvl if (present(tr_pond_topo_out)) tr_pond_topo_out = tr_pond_topo + if (present(tr_fsd_out) ) tr_fsd_out = tr_fsd if (present(tr_aero_out) ) tr_aero_out = tr_aero if (present(tr_brine_out) ) tr_brine_out = tr_brine if (present(tr_zaero_out) ) tr_zaero_out = tr_zaero @@ -382,7 +388,8 @@ subroutine icepack_write_tracer_flags(iounit) write(iounit,*) " tr_pond_cesm = ",tr_pond_cesm write(iounit,*) " tr_pond_lvl = ",tr_pond_lvl write(iounit,*) " tr_pond_topo = ",tr_pond_topo - write(iounit,*) " tr_aero = ",tr_aero + write(iounit,*) " tr_fsd = ",tr_fsd + write(iounit,*) " tr_aero = ",tr_aero write(iounit,*) " tr_brine = ",tr_brine write(iounit,*) " tr_zaero = ",tr_zaero write(iounit,*) " tr_bgc_Nit = ",tr_bgc_Nit @@ -407,7 +414,7 @@ subroutine icepack_init_tracer_indices(& nt_Tsfc_in, nt_qice_in, nt_qsno_in, nt_sice_in, & nt_fbri_in, nt_iage_in, nt_FY_in, & nt_alvl_in, nt_vlvl_in, nt_apnd_in, nt_hpnd_in, nt_ipnd_in, & - nt_aero_in, nt_zaero_in, & + nt_fsd_in, nt_aero_in, nt_zaero_in, & nt_bgc_N_in, nt_bgc_chl_in, nt_bgc_DOC_in, nt_bgc_DON_in, & nt_bgc_DIC_in, nt_bgc_Fed_in, nt_bgc_Fep_in, nt_bgc_Nit_in, nt_bgc_Am_in, & nt_bgc_Sil_in, nt_bgc_DMSPp_in, nt_bgc_DMSPd_in, nt_bgc_DMS_in, nt_bgc_hum_in, & @@ -433,6 +440,7 @@ subroutine icepack_init_tracer_indices(& nt_apnd_in, & ! melt pond area fraction nt_hpnd_in, & ! melt pond depth nt_ipnd_in, & ! melt pond refrozen lid thickness + nt_fsd_in, & ! floe size distribution nt_aero_in, & ! starting index for aerosols in ice nt_bgc_Nit_in, & ! nutrients nt_bgc_Am_in, & ! @@ -516,6 +524,7 @@ subroutine icepack_init_tracer_indices(& if (present(nt_apnd_in)) nt_apnd = nt_apnd_in if (present(nt_hpnd_in)) nt_hpnd = nt_hpnd_in if (present(nt_ipnd_in)) nt_ipnd = nt_ipnd_in + if (present(nt_fsd_in) ) nt_fsd = nt_fsd_in if (present(nt_aero_in)) nt_aero = nt_aero_in if (present(nt_bgc_Nit_in) ) nt_bgc_Nit = nt_bgc_Nit_in if (present(nt_bgc_Am_in) ) nt_bgc_Am = nt_bgc_Am_in @@ -623,7 +632,7 @@ subroutine icepack_query_tracer_indices(& nt_Tsfc_out, nt_qice_out, nt_qsno_out, nt_sice_out, & nt_fbri_out, nt_iage_out, nt_FY_out, & nt_alvl_out, nt_vlvl_out, nt_apnd_out, nt_hpnd_out, nt_ipnd_out, & - nt_aero_out, nt_zaero_out, & + nt_fsd_out, nt_aero_out, nt_zaero_out, & nt_bgc_N_out, nt_bgc_C_out, nt_bgc_chl_out, nt_bgc_DOC_out, nt_bgc_DON_out, & nt_bgc_DIC_out, nt_bgc_Fed_out, nt_bgc_Fep_out, nt_bgc_Nit_out, nt_bgc_Am_out, & nt_bgc_Sil_out, nt_bgc_DMSPp_out, nt_bgc_DMSPd_out, nt_bgc_DMS_out, nt_bgc_hum_out, & @@ -648,6 +657,7 @@ subroutine icepack_query_tracer_indices(& nt_apnd_out, & ! melt pond area fraction nt_hpnd_out, & ! melt pond depth nt_ipnd_out, & ! melt pond refrozen lid thickness + nt_fsd_out, & ! floe size distribution nt_aero_out, & ! starting index for aerosols in ice nt_bgc_Nit_out, & ! nutrients nt_bgc_Am_out, & ! @@ -718,6 +728,7 @@ subroutine icepack_query_tracer_indices(& if (present(nt_apnd_out)) nt_apnd_out = nt_apnd if (present(nt_hpnd_out)) nt_hpnd_out = nt_hpnd if (present(nt_ipnd_out)) nt_ipnd_out = nt_ipnd + if (present(nt_fsd_out) ) nt_fsd_out = nt_fsd if (present(nt_aero_out)) nt_aero_out = nt_aero if (present(nt_bgc_Nit_out) ) nt_bgc_Nit_out = nt_bgc_Nit if (present(nt_bgc_Am_out) ) nt_bgc_Am_out = nt_bgc_Am @@ -789,6 +800,7 @@ subroutine icepack_write_tracer_indices(iounit) write(iounit,*) " nt_apnd = ",nt_apnd write(iounit,*) " nt_hpnd = ",nt_hpnd write(iounit,*) " nt_ipnd = ",nt_ipnd + write(iounit,*) " nt_fsd = ",nt_fsd write(iounit,*) " nt_aero = ",nt_aero write(iounit,*) " nt_bgc_Nit = ",nt_bgc_Nit write(iounit,*) " nt_bgc_Am = ",nt_bgc_Am diff --git a/columnphysics/icepack_wavefracspec.F90 b/columnphysics/icepack_wavefracspec.F90 new file mode 100644 index 000000000..0b16f87e1 --- /dev/null +++ b/columnphysics/icepack_wavefracspec.F90 @@ -0,0 +1,723 @@ +! This module contains the subroutines required to fracture sea ice +! by ocean surface waves +! +! Theory based on: +! +! Horvat, C., & Tziperman, E. (2015). A prognostic model of the sea-ice +! floe size and thickness distribution. The Cryosphere, 9(6), 2119–2134. +! doi:10.5194/tc-9-2119-2015 +! +! and implementation described in: +! +! Roach, L. A., Horvat, C., Dean, S. M., & Bitz, C. M. (2018). An emergent +! sea ice floe size distribution in a global coupled ocean--sea ice model. +! Journal of Geophysical Research: Oceans, 123(6), 4322–4337. +! doi:10.1029/2017JC013692 +! +! now with some modifications to allow direct input of ocean surface wave spectrum. +! +! We calculate the fractures that would occur if waves enter a fully ice-covered +! region defined in one dimension in the direction of propagation, and then apply +! the outcome proportionally to the ice-covered fraction in each grid cell. Assuming +! that sea ice flexes with the sea surface height field, strains are computed on this +! sub-grid-scale 1D domain. If the strain between successive extrema exceeds a critical +! value new floes are formed with diameters equal to the distance between the extrema. +! +! authors: 2016-8 Lettie Roach, NIWA/VUW +! +! + module icepack_wavefracspec + + use icepack_kinds + use icepack_parameters, only: p01, p5, c0, c1, c2, c3, c4, c10 + use icepack_parameters, only: bignum, puny, gravit, pi + use icepack_tracers, only: nt_fsd + use icepack_warnings, only: warnstr, icepack_warnings_add + + implicit none + private + public :: icepack_init_wave, icepack_step_wavefracture + + real (kind=dbl_kind), parameter :: & + swh_minval = 0.01_dbl_kind, & ! minimum value of wave height (m) + straincrit = 3.e-5_dbl_kind, & ! critical strain + D = 1.e4_dbl_kind, & ! domain size + dx = c1, & ! domain spacing + threshold = c10 ! peak-finding threshold - + ! points are defined to be extrema if they + ! are a local max or min over a distance + ! of 10m on both sides, based on the + ! observations of Toyota et al. (2011) who + ! find this to be the order of the smallest + ! floe size affected by wave fracture + integer (kind=int_kind) :: & + nx = 10000 ! number of points in domain + +!======================================================================= + + contains + +!======================================================================= +! +! Initialize the wave spectrum and frequencies for the FSD +! +! authors: 2018 Lettie Roach, NIWA/VUW + + subroutine icepack_init_wave(nfreq, & + wave_spectrum_profile, & + wavefreq, dwavefreq) + + integer(kind=int_kind), intent(in) :: & + nfreq ! number of wave frequencies + + real(kind=dbl_kind), dimension(:), intent(out) :: & + wave_spectrum_profile, & ! ocean surface wave spectrum as a function of frequency + ! power spectral density of surface elevation, E(f) (units m^2 s) + wavefreq, & ! wave frequencies (s^-1) + dwavefreq ! wave frequency bin widths (s^-1) + + ! local variables + integer (kind=int_kind) :: k + + real(kind=dbl_kind), dimension(100) :: & + wave_spectrum_data ! default values for nfreq profile + + ! set for 25 frequencies + + wave_spectrum_data = c0 + + ! FOR TESTING ONLY - do not use for actual runs!! + wave_spectrum_data(1) = 0.00015429197810590267 + wave_spectrum_data(2) = 0.002913531381636858 + wave_spectrum_data(3) = 0.02312942035496235 + wave_spectrum_data(4) = 0.07201970368623734 + wave_spectrum_data(5) = 0.06766948103904724 + wave_spectrum_data(6) = 0.005527883302420378 + wave_spectrum_data(7) = 3.326293881400488e-05 + wave_spectrum_data(8) = 6.815936703929992e-10 + wave_spectrum_data(9) = 2.419401186610744e-20 + + do k = 1, nfreq + wave_spectrum_profile(k) = wave_spectrum_data(k) + enddo + + ! hardwired for wave coupling with NIWA version of Wavewatch + ! From Wavewatch, f(n+1) = C*f(n) where C is a constant set by the user + ! These freq are for C = 1.1 + wavefreq = (/0.04118, 0.045298, 0.0498278, 0.05481058, 0.06029164, & + 0.06632081, 0.07295289, 0.08024818, 0.08827299, 0.09710029, & + 0.10681032, 0.11749136, 0.1292405, 0.14216454, 0.15638101, & + 0.17201911, 0.18922101, 0.20814312, 0.22895744, 0.25185317, & + 0.27703848, 0.30474234, 0.33521661, 0.36873826, 0.40561208/) + + ! boundaries of bin n are at f(n)*sqrt(1/C) and f(n)*sqrt(C) + dwavefreq(:) = wavefreq(:)*(SQRT(1.1_dbl_kind) - SQRT(c1/1.1_dbl_kind)) + + end subroutine icepack_init_wave + +!======================================================================= +! +! Calculate the change in the FSD arising from wave fracture +! +! authors: 2017 Lettie Roach, NIWA/VUW +! + function get_dafsd_wave(nfsd, afsd_init, fracture_hist, frac) & + result(d_afsd) + + integer (kind=int_kind), intent(in) :: & + nfsd ! number of floe size categories + + real (kind=dbl_kind), dimension (:), intent(in) :: & + afsd_init, fracture_hist + + real (kind=dbl_kind), dimension (:,:), intent(in) :: & + frac + + ! output + real (kind=dbl_kind), dimension (nfsd) :: & + d_afsd + + ! local variables + real (kind=dbl_kind), dimension (nfsd) :: & + loss, gain, omega + + integer (kind=int_kind) :: k + + character(len=*),parameter :: subname='(get_dafsd_wave)' + + do k = 1, nfsd + ! fracture_hist is already normalized + omega(k) = afsd_init(k)*SUM(fracture_hist(1:k-1)) + end do + + loss = omega + + do k =1,nfsd + gain(k) = SUM(omega*frac(:,k)) + end do + + d_afsd(:) = gain(:) - loss(:) + + if (SUM(d_afsd(:)) > puny) then + write(warnstr,*) subname, 'area not conserved, waves' + call icepack_warnings_add(warnstr) + endif + + WHERE (ABS(d_afsd).lt.puny) d_afsd = c0 + + end function get_dafsd_wave + +!======================================================================= +! +! Adaptive timestepping for wave fracture +! See reference: Horvat & Tziperman (2017) JGR, Appendix A +! +! authors: 2018 Lettie Roach, NIWA/VUW +! +! + function get_subdt_wave(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_wave + +!======================================================================= +! +! Given fracture histogram computed from local wave spectrum, evolve +! the floe size distribution +! +! authors: 2018 Lettie Roach, NIWA/VUW +! + subroutine icepack_step_wavefracture(wave_spec_type, & + dt, ncat, nfsd, & + nfreq, & + aice, vice, aicen, & + floe_rad_l, floe_rad_c, & + 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 + + integer (kind=int_kind), intent(in) :: & + nfreq, & ! number of wave frequency categories + ncat, & ! number of thickness categories + nfsd ! number of floe size categories + + real (kind=dbl_kind), intent(in) :: & + dt, & ! time step + aice, & ! ice area fraction + vice ! ice volume per unit area + + real (kind=dbl_kind), dimension(ncat), intent(in) :: & + aicen ! ice area fraction (categories) + + real(kind=dbl_kind), dimension(:), intent(in) :: & + floe_rad_l, & ! fsd size lower bound in m (radius) + floe_rad_c ! fsd size bin centre in m (radius) + + real (kind=dbl_kind), dimension (:), intent(in) :: & + wavefreq, & ! wave frequencies (s^-1) + dwavefreq ! wave frequency bin widths (s^-1) + + real (kind=dbl_kind), dimension(:), intent(in) :: & + wave_spectrum ! ocean surface wave spectrum as a function of frequency + ! power spectral density of surface elevation, E(f) (units m^2 s) + + real (kind=dbl_kind), dimension(:,:), intent(inout) :: & + trcrn ! tracer array + + real (kind=dbl_kind), dimension(:), intent(out) :: & + d_afsd_wave ! change in fsd due to waves + + real (kind=dbl_kind), dimension(nfsd,ncat) :: & + d_afsdn_wave ! change in fsd due to waves, per category + + ! local variables + integer (kind=int_kind) :: & + n, k, t, & + nsubt ! number of subcycles + + real (kind=dbl_kind), dimension(nfsd,ncat) :: & + afsdn ! floe size and thickness distribution + + real (kind=dbl_kind), dimension (nfsd, nfsd) :: & + frac + + real (kind=dbl_kind) :: & + hbar , & ! mean ice thickness + elapsed_t , & ! elapsed subcycling time + subdt , & ! subcycling time step + cons_error ! area conservation error + + real (kind=dbl_kind), dimension (nfsd) :: & + fracture_hist, & ! fracture histogram + afsd_init , & ! tracer array + afsd_tmp , & ! tracer array + d_afsd_tmp ! change + + character(len=*),parameter :: subname='(icepack_step_wavefracture)' + + !------------------------------------ + + ! initialize + d_afsd_wave (:) = c0 + d_afsdn_wave (:,:) = c0 + fracture_hist (:) = c0 + + ! if all ice is not in first floe size category + if (.NOT. ALL(trcrn(nt_fsd,:).ge.c1-puny)) then + + ! do not try to fracture for minimal ice concentration or zero wave spectrum + if ((aice > p01).and.(MAXVAL(wave_spectrum(:)) > puny)) then + + hbar = vice / aice + + ! calculate fracture histogram + call wave_frac(nfsd, nfreq, wave_spec_type, & + floe_rad_l, floe_rad_c, & + wavefreq, dwavefreq, & + hbar, wave_spectrum, fracture_hist) + + ! 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,:) ) + + do n = 1, ncat + + afsd_init(:) = trcrn(nt_fsd:nt_fsd+nfsd-1,n) + + ! if there is ice, and a FSD, and not all ice is the smallest floe size + if ((aicen(n) > puny) .and. (SUM(afsd_init(:)) > puny) & + .and. (afsd_init(1) < c1)) then + + afsd_tmp = afsd_init + + ! frac does not vary within subcycle + frac(:,:) = c0 + do k = 2, nfsd + frac(k,1:k-1) = fracture_hist(1:k-1) + end do + do k = 1, nfsd + if (SUM(frac(k,:)) > c0) frac(k,:) = frac(k,:)/SUM(frac(k,:)) + end do + + ! adaptive sub-timestep + elapsed_t = c0 + cons_error = c0 + nsubt = 0 + DO WHILE (elapsed_t < dt) + nsubt = nsubt + 1 + + ! 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 = MIN(subdt, dt) + + ! update afsd + afsd_tmp = afsd_tmp + subdt * d_afsd_tmp(:) + + ! check conservation and negatives + if (MINVAL(afsd_tmp) < -puny) then + write(warnstr,*) subname, 'wb, <0 loop' + call icepack_warnings_add(warnstr) + endif + if (MAXVAL(afsd_tmp) > c1+puny) then + write(warnstr,*) subname, 'wb, >1 loop' + call icepack_warnings_add(warnstr) + endif + + ! update time + elapsed_t = elapsed_t + subdt + + END DO ! elapsed_t < dt + + ! In some cases---particularly for strong fracturing---the equation + ! for wave fracture does not quite conserve area. + ! With the dummy wave forcing, this happens < 2% of the time (in + ! 1997) and is always less than 10^-7. + ! Simply renormalizing may cause the first floe size + ! category to reduce, which is not physically allowed + ! to happen. So we adjust here + cons_error = SUM(afsd_tmp) - c1 + + ! area loss: add to first category + if (cons_error.lt.c0) then + afsd_tmp(1) = afsd_tmp(1) - cons_error + else + ! area gain: take it from the largest possible category + do k = nfsd, 1, -1 + if (afsd_tmp(k).gt.cons_error) then + afsd_tmp(k) = afsd_tmp(k) - cons_error + EXIT + end if + end do + end if + + ! 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,:) ) + + + ! for diagnostics + d_afsdn_wave(:,n) = afsd_tmp(:) - afsd_init(:) + d_afsd_wave (:) = d_afsd_wave(:) + aicen(n)*d_afsdn_wave(:,n) + + endif ! aicen > puny + enddo ! n + endif ! fracture occurs + + endif ! aice > p01 + end if ! all small floes + + end subroutine icepack_step_wavefracture + +!======================================================================= +! +! Calculates functions to describe the change in the FSD when waves +! fracture ice, given a wave spectrum (1D frequency, nfreq (default 25) +! frequency bins) +! +! We calculate extrema and if these are successive maximum, +! minimum, maximum or vice versa, and have strain greater than a +! critical strain, break ice and create new floes with lengths equal +! to these distances. Based on MatLab code written by Chris Horvat, +! from Horvat & Tziperman (2015). +! +! Note that a realization of sea surface height requires a random phase. +! +! authors: 2018 Lettie Roach, NIWA/VUW + + subroutine wave_frac(nfsd, nfreq, wave_spec_type, & + floe_rad_l, floe_rad_c, & + wavefreq, dwavefreq, & + hbar, spec_efreq, frac_local) + + integer (kind=int_kind), intent(in) :: & + nfsd, & ! number of floe size categories + nfreq ! number of wave frequency categories + + character (len=char_len), intent(in) :: & + wave_spec_type ! type of wave spectrum forcing + + real (kind=dbl_kind), intent(in) :: & + hbar ! mean ice thickness (m) + + real(kind=dbl_kind), dimension(:), intent(in) :: & + floe_rad_l, & ! fsd size lower bound in m (radius) + floe_rad_c ! fsd size bin centre in m (radius) + + real (kind=dbl_kind), dimension (:), intent(in) :: & + wavefreq, & ! wave frequencies (s^-1) + dwavefreq, & ! wave frequency bin widths (s^-1) + spec_efreq ! wave spectrum (m^2 s) + + real (kind=dbl_kind), dimension (nfsd), intent(out) :: & + frac_local ! fracturing histogram + + ! local variables + + integer (kind=int_kind) :: i, j, k + + integer, parameter :: & + loopcts = 1 ! number of SSH realizations + + real (kind=dbl_kind), dimension(nfreq) :: & + spec_elambda, & ! spectrum as a function of wavelength (m^-1 s^-1) + reverse_spec_elambda, & ! reversed + reverse_lambda, lambda, & ! wavelengths (m) + reverse_dlambda, dlambda, & ! wavelength bin spacing (m) + spec_coeff, & + phi, rand_array, summand + + real (kind=dbl_kind), dimension(2*nx) :: & + fraclengths + + real (kind=dbl_kind), dimension(nx) :: & + X, & ! spatial domain (m) + eta ! sea surface height field (m) + + real (kind=dbl_kind), dimension(nfsd) :: & + frachistogram + + logical (kind=log_kind) :: & + e_stop ! if true, stop and return zero omega and fsdformed + + e_stop = .false. ! if true, aborts fracture calc + + ! spatial domain + do j = 1, nx + X(j)= j*dx + end do + + ! dispersion relation + reverse_lambda (:) = gravit/(c2*pi*wavefreq (:)**c2) + reverse_dlambda(:) = gravit/(c2*pi*dwavefreq(:)**c2) + ! convert to lambda spectrum + reverse_spec_elambda(:) = spec_efreq(:) & + *(p5 * (gravit/(c2*pi*reverse_lambda(:)**c3) )**p5) + ! reverse lambda + lambda (:) = reverse_lambda (nfreq:1:-1) + dlambda(:) = reverse_dlambda(nfreq:1:-1) + spec_elambda(:) = reverse_spec_elambda(nfreq:1:-1) + + ! spectral coefficients + spec_coeff = sqrt(c2*spec_elambda*dlambda) + + ! initialize fracture lengths + fraclengths(:) = c0 + + ! loop over n. realizations of SSH + do i = 1, loopcts + + ! 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)=='random') then + call RANDOM_NUMBER(rand_array) + else + rand_array(:) = p5 + endif + phi = c2*pi*rand_array + + do j = 1, nx + ! SSH field in space (sum over wavelengths, no attenuation) + summand = spec_coeff*COS(2*pi*X(j)/lambda+phi) + eta(j) = SUM(summand) + end do + + if ((SUM(ABS(eta)) > puny).and.(hbar > puny)) then + call get_fraclengths(X, eta, fraclengths, hbar, e_stop) + end if + end do + + frachistogram(:) = c0 + + if (.not. e_stop) then + + ! convert from diameter to radii + fraclengths(:) = fraclengths(:)/c2 + + ! bin into FS cats + ! highest cat cannot be fractured into + do j = 1, size(fraclengths) + 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 + + end do + end if + + do k = 1, nfsd + frac_local(k) = floe_rad_c(k)*frachistogram(k) + end do + + ! normalize + if (SUM(frac_local) /= c0) frac_local(:) = frac_local(:) / SUM(frac_local(:)) + + end subroutine wave_frac + +!=========================================================================== +! +! Given the (attenuated) sea surface height, find the strain across triplets +! of max, min, max or min, max, min (local extrema within 10m). +! If this strain is greater than the critical strain, ice can fracture +! and new floes are formed with sizes equal to the distances between +! extrema. Based on MatLab code written by Chris Horvat, +! from Horvat & Tziperman (2015). +! +! authors: 2016 Lettie Roach, NIWA/VUW +! + subroutine get_fraclengths(X, eta, fraclengths, hbar, e_stop) + + real (kind=dbl_kind) :: & + hbar ! mean thickness (m) + + real (kind=dbl_kind), intent(in), dimension (nx) :: & + X, & ! spatial domain (m) + eta ! sea surface height field (m) + + real (kind=dbl_kind), intent(inout), dimension (2*nx) :: & + fraclengths ! the biggest number of fraclengths we could have is + ! two floe pieces created at each subgridpoint ie. 2*nx + ! This will never actually happen - most of the array + ! will be zeros + + logical (kind=log_kind), intent(inout) :: & + e_stop ! if true, stop and return zero omega and fsdformed + + ! local variables + integer (kind=int_kind) :: & + spcing, & ! distance over which to search for extrema on each side of point + j, k, & ! indices to iterate over domain + first, last, & ! indices over which to search for extrema + j_neg, & ! nearest extrema backwards + j_pos, & ! nearest extrema forwards + n_above ! number of points where strain is above critical + + real (kind=dbl_kind), dimension(nx) :: & + strain, & ! the strain between triplets of extrema + frac_size_one, & ! + frac_size_two + + logical (kind=log_kind), dimension(nx) :: & + is_max, is_min,& ! arrays to hold whether each point is a local max or min + is_extremum, & ! or extremum + is_triplet ! or triplet of extrema + + real (kind=dbl_kind) :: & + delta, & ! difference in x between current and prev extrema + delta_pos ! difference in x between next and current extrema + + integer (kind=int_kind), dimension(1) :: & + maxj, minj ! indices of local max and min + + ! ------- equivalent of peakfinder2 + ! given eta and spcing, compute extremelocs in ascending order + spcing = nint(threshold/dx) + + is_max = .false. + is_min = .false. + is_extremum = .false. + is_triplet = .false. + strain = c0 + frac_size_one = c0 + frac_size_two = c0 + j_neg = 0 + j_pos = 0 + fraclengths(:) = c0 + + ! search for local max and min within spacing + ! on either side of each point + + do j = 1, nx + + first = MAX(1,j-spcing) + last = MIN(nx,j+spcing) + + maxj = MAXLOC(eta(first:last)) + minj = MINLOC(eta(first:last)) + +! if (COUNT(eta(first:last) == MAXVAL(eta(first:last))) > 1) & +! stop 'more than one max' +! if (COUNT(eta(first:last) == MINVAL(eta(first:last))) > 1) & +! stop 'more than one min' + + if (maxj(1)+first-1 == j) is_max(j) = .true. + if (minj(1)+first-1 == j) is_min(j) = .true. + +! if (is_min(j).and.is_max(j)) then +! print *, 'X ',X +! print *, 'eta ',eta +! print *, 'frst last' ,first, last +! print *, 'maxj, minj ',maxj,minj +! stop 'error in extrema' +! end if + if (is_min(j).or.is_max(j)) is_extremum(j) = .true. + end do + + do j = 2, nx-1 + if (is_extremum(j)) then + if (j == 2) then + if (is_extremum(1)) j_neg = 1 + else + do k = j-1, 1, -1 + if (is_extremum(k)) then + j_neg = k + EXIT + end if + end do + end if + + do k = j+1, nx + if (is_extremum(k)) then + j_pos = k + EXIT + end if + end do + + if ((j_neg > 0).and.(j_pos > 0)) then + if (is_max(j_neg).and.is_min(j).and.is_max(j_pos)) & + is_triplet(j) = .true. + if (is_min(j_neg).and.is_max(j).and.is_min(j_pos)) & + is_triplet(j) = .true. + end if + + ! calculate strain + if (is_triplet(j)) then + delta_pos = X(j_pos) - X(j ) + delta = X(j ) - X(j_neg) + + strain(j) = p5*hbar*(eta(j_neg) - eta(j)) & + / (delta*(delta+delta_pos)) + + if (strain(j) > straincrit) then + frac_size_one(j) = X(j_pos) - X(j ) + frac_size_two(j) = X(j ) - X(j_neg) + end if + end if + end if + + end do + + n_above = COUNT(strain > straincrit) + if (n_above > 0) then + fraclengths(1:nx) = frac_size_one(:) + fraclengths(nx+1:2*nx) = frac_size_two(:) + e_stop = .false. + else + e_stop = .true. + end if + + end subroutine get_fraclengths + +!======================================================================= + + end module icepack_wavefracspec + +!======================================================================= + + diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90 index 0c22d5bce..c5cfd35d5 100644 --- a/configuration/driver/icedrv_InitMod.F90 +++ b/configuration/driver/icedrv_InitMod.F90 @@ -29,15 +29,18 @@ module icedrv_InitMod subroutine icedrv_initialize use icedrv_arrays_column, only: hin_max, c_hi_range + use icedrv_arrays_column, only: floe_rad_l, floe_rad_c, & + floe_binwidth, c_fsd_range use icedrv_calendar, only: dt, time, istep, istep1, & init_calendar, calendar use icepack_intfc, only: icepack_init_itd, icepack_init_itd_hist + use icepack_intfc, only: icepack_init_fsd_bounds use icepack_intfc, only: icepack_warnings_flush - use icedrv_domain_size, only: ncat + use icedrv_domain_size, only: ncat, nfsd ! use icedrv_diagnostics, only: icedrv_diagnostics_debug use icedrv_flux, only: init_coupler_flux, init_history_therm, & init_flux_atm_ocn - use icedrv_forcing, only: init_forcing, get_forcing + use icedrv_forcing, only: init_forcing, get_forcing, get_wave_spec use icedrv_forcing_bgc, only: get_forcing_bgc, faero_default, init_forcing_bgc use icedrv_restart_shared, only: restart use icedrv_init, only: input_data, init_state, init_grid2 @@ -48,7 +51,8 @@ subroutine icedrv_initialize skl_bgc, & ! from icepack z_tracers, & ! from icepack tr_aero, & ! from icepack - tr_zaero ! from icepack + tr_zaero, & ! from icepack + tr_fsd, wave_spec character(len=*), parameter :: subname='(icedrv_initialize)' @@ -72,6 +76,19 @@ subroutine icedrv_initialize call icepack_init_itd_hist(ncat=ncat, c_hi_range=c_hi_range, hin_max=hin_max) ! output + call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted(subname)) then + 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 icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted(subname)) then call icedrv_system_abort(file=__FILE__,line=__LINE__) @@ -96,6 +113,7 @@ subroutine icedrv_initialize !-------------------------------------------------------------------- call icepack_query_parameters(skl_bgc_out=skl_bgc) call icepack_query_parameters(z_tracers_out=z_tracers) + call icepack_query_parameters(wave_spec_out=wave_spec) call icepack_query_tracer_flags(tr_aero_out=tr_aero) call icepack_query_tracer_flags(tr_zaero_out=tr_zaero) call icepack_warnings_flush(nu_diag) @@ -104,6 +122,7 @@ subroutine icedrv_initialize call init_forcing ! initialize forcing (standalone) if (skl_bgc .or. z_tracers) call init_forcing_bgc !cn + if (tr_fsd .and. wave_spec) call get_wave_spec ! wave spectrum in ice call get_forcing(istep1) ! get forcing from data arrays ! aerosols @@ -141,7 +160,8 @@ subroutine init_restart skl_bgc, & ! from icepack z_tracers, & ! from icepack solve_zsal, & ! from icepack - tr_brine ! from icepack + tr_brine, & ! from icepack + tr_fsd ! from icepack character(len=*), parameter :: subname='(init_restart)' @@ -152,7 +172,7 @@ subroutine init_restart call icepack_query_parameters(skl_bgc_out=skl_bgc) call icepack_query_parameters(z_tracers_out=z_tracers) call icepack_query_parameters(solve_zsal_out=solve_zsal) - call icepack_query_tracer_flags(tr_brine_out=tr_brine) + call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_fsd_out=tr_fsd) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) @@ -169,6 +189,10 @@ subroutine init_restart endif if (solve_zsal .or. skl_bgc .or. z_tracers) then + if (tr_fsd) then + write (nu_diag,*) 'FSD implementation incomplete for use with BGC' + call icedrv_system_abort(string=subname,file=__FILE__,line=__LINE__) + endif call init_bgc if (restart) call read_restart_bgc ! complete BGC initialization endif diff --git a/configuration/driver/icedrv_RunMod.F90 b/configuration/driver/icedrv_RunMod.F90 index e89a5eaba..3f86cacd5 100644 --- a/configuration/driver/icedrv_RunMod.F90 +++ b/configuration/driver/icedrv_RunMod.F90 @@ -32,11 +32,12 @@ module icedrv_RunMod subroutine icedrv_run use icedrv_calendar, only: istep, istep1, time, dt, stop_now, calendar - use icedrv_forcing, only: get_forcing + use icedrv_forcing, only: get_forcing, get_wave_spec use icedrv_forcing_bgc, only: faero_default, get_forcing_bgc use icedrv_flux, only: init_flux_atm_ocn - logical (kind=log_kind) :: skl_bgc, z_tracers, tr_aero, tr_zaero + logical (kind=log_kind) :: skl_bgc, z_tracers, tr_aero, tr_zaero, & + wave_spec, tr_fsd character(len=*), parameter :: subname='(icedrv_run)' @@ -44,7 +45,8 @@ subroutine icedrv_run ! timestep loop !-------------------------------------------------------------------- - call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero) + call icepack_query_tracer_flags(tr_aero_out=tr_aero, tr_zaero_out=tr_zaero, & + tr_fsd_out=tr_fsd) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) @@ -61,16 +63,18 @@ subroutine icedrv_run if (stop_now >= 1) exit timeLoop + call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers,& + wave_spec_out=wave_spec) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & + file=__FILE__,line= __LINE__) + + if (tr_fsd .and. wave_spec) call get_wave_spec ! wave spectrum in ice call get_forcing(istep1) ! get forcing from data arrays ! aerosols if (tr_aero .or. tr_zaero) call faero_default ! default values - call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers) - call icepack_warnings_flush(nu_diag) - if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & - file=__FILE__,line= __LINE__) - if (skl_bgc .or. z_tracers) call get_forcing_bgc ! biogeochemistry call init_flux_atm_ocn ! initialize atmosphere, ocean fluxes @@ -97,13 +101,14 @@ subroutine ice_step use icedrv_restart_bgc, only: write_restart_bgc use icedrv_step, only: prep_radiation, step_therm1, step_therm2, & update_state, step_dyn_ridge, step_radiation, & - biogeochemistry + biogeochemistry, step_dyn_wave integer (kind=int_kind) :: & k ! dynamics supercycling index logical (kind=log_kind) :: & - calc_Tsfc, skl_bgc, solve_zsal, z_tracers, tr_brine ! from icepack + calc_Tsfc, skl_bgc, solve_zsal, z_tracers, tr_brine, & ! from icepack + tr_fsd, wave_spec real (kind=dbl_kind) :: & offset ! d(age)/dt time offset @@ -117,8 +122,10 @@ subroutine ice_step !----------------------------------------------------------------- call icepack_query_parameters(skl_bgc_out=skl_bgc, z_tracers_out=z_tracers) - call icepack_query_parameters(solve_zsal_out=solve_zsal, calc_Tsfc_out=calc_Tsfc) - call icepack_query_tracer_flags(tr_brine_out=tr_brine) + call icepack_query_parameters(solve_zsal_out=solve_zsal, & + calc_Tsfc_out=calc_Tsfc, & + wave_spec_out=wave_spec) + call icepack_query_tracer_flags(tr_brine_out=tr_brine,tr_fsd_out=tr_fsd) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) @@ -159,6 +166,10 @@ subroutine ice_step call init_history_dyn + ! wave fracture of the floe size distribution + ! note this is called outside of the dynamics subcycling loop + if (tr_fsd .and. wave_spec) call step_dyn_wave(dt) + do k = 1, ndtd ! ridging diff --git a/configuration/driver/icedrv_arrays_column.F90 b/configuration/driver/icedrv_arrays_column.F90 index fef4d7bdf..fb372cc37 100644 --- a/configuration/driver/icedrv_arrays_column.F90 +++ b/configuration/driver/icedrv_arrays_column.F90 @@ -8,7 +8,7 @@ module icedrv_arrays_column use icedrv_kinds use icedrv_constants, only: nu_diag - use icedrv_domain_size, only: nx, ncat, nilyr, nslyr + use icedrv_domain_size, only: nx, ncat, nilyr, nslyr, nfsd, nfreq use icedrv_domain_size, only: nblyr, max_nsw , max_ntrcr use icepack_intfc, only: icepack_max_nbtrcr, icepack_max_algae, icepack_max_aero use icepack_intfc, only: icepack_nmodal1, icepack_nmodal2 @@ -254,6 +254,31 @@ module icedrv_arrays_column ice_bio_net, & ! depth integrated tracer (mmol/m^2) snow_bio_net ! depth integrated snow tracer (mmol/m^2) + ! floe size distribution + real(kind=dbl_kind), dimension(nfsd), public :: & + 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 (nx), public :: & + wave_sig_ht ! significant height of waves (m) + + real (kind=dbl_kind), dimension (nfreq), public :: & + wavefreq, & ! wave frequencies + dwavefreq ! wave frequency bin widths + + real (kind=dbl_kind), dimension (nx,nfreq), public :: & + wave_spectrum ! wave spectrum + + real (kind=dbl_kind), dimension (nx,nfsd), public :: & + ! change in floe size distribution due to processes + d_afsd_newi, d_afsd_latg, d_afsd_latm, d_afsd_wave, d_afsd_weld + + character (len=35), public, dimension(nfsd) :: & + c_fsd_range ! fsd floe_rad bounds (m) + + + !======================================================================= end module icedrv_arrays_column diff --git a/configuration/driver/icedrv_diagnostics.F90 b/configuration/driver/icedrv_diagnostics.F90 index 5962e6812..123ac8c56 100644 --- a/configuration/driver/icedrv_diagnostics.F90 +++ b/configuration/driver/icedrv_diagnostics.F90 @@ -55,13 +55,15 @@ module icedrv_diagnostics subroutine runtime_diags (dt) + use icedrv_arrays_column, only: floe_rad_c + use icedrv_domain_size, only: ncat, nfsd use icedrv_flux, only: evap, fsnow, frazil use icedrv_flux, only: fswabs, flw, flwout, fsens, fsurf, flat use icedrv_flux, only: frain use icedrv_flux, only: Tair, Qa, fsw, fcondtop use icedrv_flux, only: meltt, meltb, meltl, snoice use icedrv_flux, only: dsnow, congel, sst, sss, Tf, fhocn - use icedrv_state, only: aice, vice, vsno, trcr + use icedrv_state, only: aice, vice, vsno, trcr, trcrn, aicen real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -69,16 +71,16 @@ subroutine runtime_diags (dt) ! local variables integer (kind=int_kind) :: & - n + n, nc, k logical (kind=log_kind) :: & - calc_Tsfc + calc_Tsfc, tr_fsd ! fields at diagnostic points real (kind=dbl_kind) :: & pTair, pfsnow, pfrain, & paice, hiavg, hsavg, hbravg, psalt, pTsfc, & - pevap, pfhocn + pevap, pfhocn, fsdavg real (kind=dbl_kind), dimension (nx) :: & work1, work2 @@ -87,7 +89,7 @@ subroutine runtime_diags (dt) Tffresh, rhos, rhow, rhoi logical (kind=log_kind) :: tr_brine - integer (kind=int_kind) :: nt_fbri, nt_Tsfc + integer (kind=int_kind) :: nt_fbri, nt_Tsfc, nt_fsd character(len=*), parameter :: subname='(runtime_diags)' @@ -96,8 +98,9 @@ subroutine runtime_diags (dt) !----------------------------------------------------------------- call icepack_query_parameters(calc_Tsfc_out=calc_Tsfc) - call icepack_query_tracer_flags(tr_brine_out=tr_brine) - call icepack_query_tracer_indices(nt_fbri_out=nt_fbri, nt_Tsfc_out=nt_Tsfc) + call icepack_query_tracer_flags(tr_brine_out=tr_brine,tr_fsd_out=tr_fsd) + call icepack_query_tracer_indices(nt_fbri_out=nt_fbri, nt_Tsfc_out=nt_Tsfc,& + nt_fsd_out=nt_fsd) call icepack_query_parameters(Tffresh_out=Tffresh, rhos_out=rhos, & rhow_out=rhow, rhoi_out=rhoi) call icepack_warnings_flush(nu_diag) @@ -118,6 +121,7 @@ subroutine runtime_diags (dt) paice = aice(n) ! ice area hiavg = c0 ! avg snow/ice thickness + fsdavg = c0 ! FSD rep radius hsavg = c0 hbravg = c0 ! avg brine thickness psalt = c0 @@ -125,6 +129,16 @@ subroutine runtime_diags (dt) hiavg = vice(n)/paice hsavg = vsno(n)/paice if (tr_brine) hbravg = trcr(n,nt_fbri)* hiavg + if (tr_fsd) then + do nc = 1, ncat + do k = 1, nfsd + fsdavg = fsdavg & + + trcrn(n,nt_fsd+k-1,nc) * floe_rad_c(k) & + * aicen(n,nc) / paice + end do + end do + end if + endif if (vice(n) /= c0) psalt = work2(n)/vice(n) pTsfc = trcr(n,nt_Tsfc) ! ice/snow sfc temperature @@ -160,7 +174,10 @@ subroutine runtime_diags (dt) write(nu_diag_out+n-1,900) 'avg snow depth (m) = ',hsavg write(nu_diag_out+n-1,900) 'avg salinity (ppt) = ',psalt write(nu_diag_out+n-1,900) 'avg brine thickness (m)= ',hbravg - + if (tr_fsd) & + write(nu_diag_out+n-1,900) 'avg fsd rep radius (m) = ',fsdavg + + if (calc_Tsfc) then write(nu_diag_out+n-1,900) 'surface temperature(C) = ',pTsfc ! ice/snow write(nu_diag_out+n-1,900) 'absorbed shortwave flx = ',fswabs(n) @@ -346,7 +363,8 @@ subroutine icedrv_diagnostics_debug (plabeld) character (*), intent(in) :: plabeld - character(len=*), parameter :: subname='(icedrv_diagnostics_debug)' + character(len=*), parameter :: & + subname='(icedrv_diagnostics_debug)' ! printing info for routine print_state @@ -372,7 +390,7 @@ end subroutine icedrv_diagnostics_debug subroutine print_state(plabel,i) use icedrv_calendar, only: istep1, time - use icedrv_domain_size, only: ncat, nilyr, nslyr + use icedrv_domain_size, only: ncat, nilyr, nslyr, nfsd use icedrv_state, only: aice0, aicen, vicen, vsnon, uvel, vvel, trcrn use icedrv_flux, only: uatm, vatm, potT, Tair, Qa, flw, frain, fsnow use icedrv_flux, only: fsens, flat, evap, flwout @@ -397,7 +415,9 @@ subroutine print_state(plabel,i) integer (kind=int_kind) :: n, k - integer (kind=int_kind) :: nt_Tsfc, nt_qice, nt_qsno + integer (kind=int_kind) :: nt_Tsfc, nt_qice, nt_qsno, nt_fsd + + logical (kind=log_kind) :: tr_fsd character(len=*), parameter :: subname='(print_state)' @@ -405,8 +425,9 @@ subroutine print_state(plabel,i) ! query Icepack values !----------------------------------------------------------------- + call icepack_query_tracer_flags(tr_fsd_out=tr_fsd) call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_qice_out=nt_qice, & - nt_qsno_out=nt_qsno) + nt_qsno_out=nt_qsno,nt_fsd_out=nt_fsd) call icepack_query_parameters(puny_out=puny, Lfresh_out=Lfresh, cp_ice_out=cp_ice, & rhoi_out=rhoi, rhos_out=rhos) call icepack_warnings_flush(nu_diag) @@ -433,6 +454,7 @@ subroutine print_state(plabel,i) write(nu_diag,*) 'hsn', vsnon(i,n)/aicen(i,n) endif write(nu_diag,*) 'Tsfcn',trcrn(i,nt_Tsfc,n) + if (tr_fsd) write(nu_diag,*) 'afsdn',trcrn(i,nt_fsd,n) ! fsd cat 1 write(nu_diag,*) ' ' enddo ! n diff --git a/configuration/driver/icedrv_domain_size.F90 b/configuration/driver/icedrv_domain_size.F90 index f13db9997..37b825e27 100644 --- a/configuration/driver/icedrv_domain_size.F90 +++ b/configuration/driver/icedrv_domain_size.F90 @@ -15,6 +15,7 @@ module icedrv_domain_size integer (kind=int_kind), parameter, public :: & nx = NXGLOB , & ! vector length + nfsd = NFSDCAT , & ! number of floe size categories ncat = NICECAT , & ! number of categories nilyr = NICELYR , & ! number of ice layers per category nslyr = NSNWLYR , & ! number of snow layers per category @@ -26,6 +27,7 @@ module icedrv_domain_size n_don = TRDON , & ! number of DON pools in use n_fed = TRFED , & ! number of Fe pools in use dissolved Fe n_fep = TRFEP , & ! number of Fe pools in use particulate Fe + nfreq = 25 , & ! number of wave frequencies ! HARDWIRED FOR NOW nblyr = NBGCLYR , & ! number of bio/brine layers per category ! maximum number of biology tracers + aerosols ! *** add to kscavz in icepack_zbgc_shared.F90 @@ -40,6 +42,7 @@ module icedrv_domain_size + nilyr & ! ice enthalpy + nslyr & ! snow enthalpy !!!!! optional tracers: + + nfsd & ! number of floe size categories + TRAGE & ! age + TRFY & ! first-year area + TRLVL*2 & ! level/deformed ice diff --git a/configuration/driver/icedrv_flux.F90 b/configuration/driver/icedrv_flux.F90 index 4cbd7388e..0abecd448 100644 --- a/configuration/driver/icedrv_flux.F90 +++ b/configuration/driver/icedrv_flux.F90 @@ -266,6 +266,7 @@ module icedrv_flux real (kind=dbl_kind), dimension (nx), public :: & rside , & ! fraction of ice that melts laterally + fside , & ! lateral heat flux (W/m^2) fsw , & ! incoming shortwave radiation (W/m^2) coszen , & ! cosine solar zenith angle, < 0 for sun below horizon rdg_conv, & ! convergence term for ridging (1/s) diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90 index 126a78235..c0d83fb9f 100644 --- a/configuration/driver/icedrv_forcing.F90 +++ b/configuration/driver/icedrv_forcing.F90 @@ -14,7 +14,8 @@ module icedrv_forcing use icedrv_constants, only: c0, c1, c2, c10, c100, p5, c4, c24 use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted use icepack_intfc, only: icepack_query_parameters - use icepack_intfc, only: icepack_sea_freezing_temperature + use icepack_intfc, only: icepack_sea_freezing_temperature + use icepack_intfc, only: icepack_init_wave use icedrv_system, only: icedrv_system_abort use icedrv_flux, only: zlvl, Tair, potT, rhoa, uatm, vatm, wind, & strax, stray, fsw, swvdr, swvdf, swidr, swidf, Qa, flw, frain, & @@ -22,7 +23,8 @@ module icedrv_forcing implicit none private - public :: init_forcing, get_forcing, interp_coeff, interp_coeff_monthly + public :: init_forcing, get_forcing, interp_coeff, & + interp_coeff_monthly, get_wave_spec integer (kind=int_kind), parameter :: & ntime = 8760 ! number of data points in time @@ -1122,6 +1124,40 @@ subroutine ice_open_clos end subroutine ice_open_clos +!======================================================================= + + subroutine get_wave_spec + + use icedrv_arrays_column, only: wave_spectrum, wave_sig_ht, & + dwavefreq, wavefreq + use icedrv_domain_size, only: nfreq + +#ifdef ncdf + use netcdf +#endif + + ! local variables + integer (kind=int_kind) :: & + k + + real(kind=dbl_kind), dimension(nfreq) :: & + wave_spectrum_profile ! wave spectrum + + wave_spectrum(:,:) = c0 + + ! wave spectrum and frequencies + ! get hardwired frequency bin info and a dummy wave spectrum profile + call icepack_init_wave(nfreq=nfreq, & + wave_spectrum_profile=wave_spectrum_profile, & + wavefreq=wavefreq, dwavefreq=dwavefreq) + + do k = 1, nfreq + wave_spectrum(:,k) = wave_spectrum_profile(k) + enddo + + end subroutine get_wave_spec + + !======================================================================= end module icedrv_forcing diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90 index b7988ea83..67160469c 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/configuration/driver/icedrv_init.F90 @@ -11,6 +11,7 @@ module icedrv_init use icedrv_constants, only: c0, c1, c2, c3, p2, p5 use icedrv_domain_size, only: nx use icepack_intfc, only: icepack_init_parameters + use icepack_intfc, only: icepack_init_fsd use icepack_intfc, only: icepack_init_tracer_flags use icepack_intfc, only: icepack_init_tracer_numbers use icepack_intfc, only: icepack_init_tracer_indices @@ -56,7 +57,7 @@ module icedrv_init subroutine input_data use icedrv_diagnostics, only: diag_file, nx_names - use icedrv_domain_size, only: nilyr, nslyr, max_ntrcr, ncat, n_aero + use icedrv_domain_size, only: nilyr, nslyr, max_ntrcr, ncat, n_aero, nfsd use icedrv_calendar, only: year_init, istep0 use icedrv_calendar, only: dumpfreq, diagfreq, dump_last use icedrv_calendar, only: npt, dt, ndtd, days_per_year, use_leap_years @@ -92,15 +93,16 @@ subroutine input_data natmiter, kitd, kcatbound character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, & - tfrz_option, frzpnd, atmbndy + tfrz_option, frzpnd, atmbndy, wave_spec_type logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair integer (kind=int_kind) :: ntrcr - logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_pond, tr_aero - logical (kind=log_kind) :: tr_pond_cesm, tr_pond_lvl, tr_pond_topo + logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_pond, tr_aero, tr_fsd + logical (kind=log_kind) :: tr_pond_cesm, tr_pond_lvl, tr_pond_topo, wave_spec integer (kind=int_kind) :: nt_Tsfc, nt_sice, nt_qice, nt_qsno, nt_iage, nt_FY - integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, nt_aero + integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, & + nt_aero, nt_fsd real (kind=real_kind) :: rpcesm, rplvl, rptopo real (kind=dbl_kind) :: Cf, puny @@ -144,7 +146,7 @@ subroutine input_data update_ocn_f, l_mpond_fresh, ustar_min, & fbot_xfer_type, oceanmixed_ice, emissivity, & formdrag, highfreq, natmiter, & - tfrz_option, default_season, & + tfrz_option, default_season, wave_spec_type, & precip_units, fyear_init, ycycle, & atm_data_type, ocn_data_type, bgc_data_type, & atm_data_file, ocn_data_file, bgc_data_file, & @@ -159,7 +161,8 @@ subroutine input_data tr_pond_cesm, & tr_pond_lvl, & tr_pond_topo, & - tr_aero + tr_aero, & + tr_fsd !----------------------------------------------------------------- ! query Icepack values @@ -189,7 +192,8 @@ subroutine input_data phi_c_slow_mode_out=phi_c_slow_mode, & phi_i_mushy_out=phi_i_mushy, & tfrz_option_out=tfrz_option, kalg_out=kalg, & - fbot_xfer_type_out=fbot_xfer_type, puny_out=puny) + fbot_xfer_type_out=fbot_xfer_type, puny_out=puny, & + wave_spec_type_out=wave_spec_type) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__, line=__LINE__) @@ -227,6 +231,7 @@ subroutine input_data precip_units = 'mks' ! 'mm_per_month' or ! 'mm_per_sec' = 'mks' = kg/m^2 s oceanmixed_ice = .false. ! if true, use internal ocean mixed layer + wave_spec_type = 'none' ! type of wave spectrum forcing ocn_data_format = 'bin' ! file format ('bin'=binary or 'nc'=netcdf) ocn_data_type = 'default' ! source of ocean forcing data ocn_data_file = ' ' ! ocean forcing data file @@ -246,6 +251,7 @@ subroutine input_data tr_pond_lvl = .false. ! level-ice melt ponds tr_pond_topo = .false. ! explicit melt ponds (topographic) tr_aero = .false. ! aerosols + tr_fsd = .false. ! floe size distribution !----------------------------------------------------------------- ! read from input file @@ -444,6 +450,9 @@ subroutine input_data fbot_xfer_type = 'constant' endif + wave_spec = .false. + if (tr_fsd .and. (trim(wave_spec_type) /= 'none')) wave_spec = .true. + !----------------------------------------------------------------- ! spew !----------------------------------------------------------------- @@ -559,9 +568,12 @@ subroutine input_data write(nu_diag,*) ' default_season = ', trim(default_season) write(nu_diag,1010) ' update_ocn_f = ', update_ocn_f + write(nu_diag,1010) ' wave_spec = ', wave_spec + if (wave_spec) & + write(nu_diag,*) ' wave_spec_type = ', wave_spec_type write(nu_diag,1010) ' l_mpond_fresh = ', l_mpond_fresh write(nu_diag,1005) ' ustar_min = ', ustar_min - write(nu_diag, *) ' fbot_xfer_type = ', & + write(nu_diag,*) ' fbot_xfer_type = ', & trim(fbot_xfer_type) write(nu_diag,1010) ' oceanmixed_ice = ', oceanmixed_ice write(nu_diag,*) ' tfrz_option = ', & @@ -578,6 +590,7 @@ subroutine input_data write(nu_diag,1010) ' tr_pond_lvl = ', tr_pond_lvl write(nu_diag,1010) ' tr_pond_topo = ', tr_pond_topo write(nu_diag,1010) ' tr_aero = ', tr_aero + write(nu_diag,1010) ' tr_fsd = ', tr_fsd nt_Tsfc = 1 ! index tracers, starting with Tsfc = 1 ntrcr = 1 ! count tracers, starting with Tsfc = 1 @@ -630,6 +643,12 @@ subroutine input_data endif endif + nt_fsd = max_ntrcr + if (tr_fsd) then + nt_fsd = ntrcr + 1 ! floe size distribution + ntrcr = ntrcr + nfsd + end if + nt_aero = max_ntrcr - 4*n_aero if (tr_aero) then nt_aero = ntrcr + 1 @@ -653,6 +672,7 @@ subroutine input_data write(nu_diag,1020) 'ncat', ncat write(nu_diag,1020) 'nilyr', nilyr write(nu_diag,1020) 'nslyr', nslyr + write(nu_diag,1020) 'nfsd', nfsd write(nu_diag,*)' ' write(nu_diag,1020) 'nx', nx write(nu_diag,*)' ' @@ -711,19 +731,20 @@ subroutine input_data phi_c_slow_mode_in=phi_c_slow_mode, & phi_i_mushy_in=phi_i_mushy, & tfrz_option_in=tfrz_option, kalg_in=kalg, & - fbot_xfer_type_in=fbot_xfer_type) + fbot_xfer_type_in=fbot_xfer_type, & + wave_spec_type_in=wave_spec_type, wave_spec_in=wave_spec) call icepack_init_tracer_numbers(ntrcr_in=ntrcr) call icepack_init_tracer_flags(tr_iage_in=tr_iage, & tr_FY_in=tr_FY, tr_lvl_in=tr_lvl, tr_aero_in=tr_aero, & tr_pond_in=tr_pond, tr_pond_cesm_in=tr_pond_cesm, & tr_pond_lvl_in=tr_pond_lvl, & - tr_pond_topo_in=tr_pond_topo) + tr_pond_topo_in=tr_pond_topo, tr_fsd_in=tr_fsd) call icepack_init_tracer_indices(nt_Tsfc_in=nt_Tsfc, & nt_sice_in=nt_sice, nt_qice_in=nt_qice, & nt_qsno_in=nt_qsno, nt_iage_in=nt_iage, & nt_fy_in=nt_fy, nt_alvl_in=nt_alvl, nt_vlvl_in=nt_vlvl, & nt_apnd_in=nt_apnd, nt_hpnd_in=nt_hpnd, nt_ipnd_in=nt_ipnd, & - nt_aero_in=nt_aero) + nt_aero_in=nt_aero, nt_fsd_in=nt_fsd) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & @@ -790,7 +811,7 @@ end subroutine init_grid2 subroutine init_state use icepack_intfc, only: icepack_aggregate - use icedrv_domain_size, only: ncat, nilyr, nslyr, max_ntrcr, n_aero + use icedrv_domain_size, only: ncat, nilyr, nslyr, max_ntrcr, n_aero, nfsd use icedrv_flux, only: sst, Tf, Tair, salinz, Tmltz use icedrv_state, only: trcr_depend, aicen, trcrn, vicen, vsnon use icedrv_state, only: aice0, aice, vice, vsno, trcr, aice_init @@ -805,10 +826,11 @@ subroutine init_state heat_capacity ! from icepack integer (kind=int_kind) :: ntrcr - logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_aero + logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_aero, tr_fsd logical (kind=log_kind) :: tr_pond_cesm, tr_pond_lvl, tr_pond_topo integer (kind=int_kind) :: nt_Tsfc, nt_sice, nt_qice, nt_qsno, nt_iage, nt_fy - integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, nt_aero + integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, & + nt_ipnd, nt_aero, nt_fsd character(len=*), parameter :: subname='(init_state)' @@ -821,13 +843,13 @@ subroutine init_state call icepack_query_tracer_flags(tr_iage_out=tr_iage, & tr_FY_out=tr_FY, tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, & tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, & - tr_pond_topo_out=tr_pond_topo) + tr_pond_topo_out=tr_pond_topo, tr_fsd_out=tr_fsd) call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, & nt_sice_out=nt_sice, nt_qice_out=nt_qice, & nt_qsno_out=nt_qsno, nt_iage_out=nt_iage, nt_fy_out=nt_fy, & nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, & nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, & - nt_ipnd_out=nt_ipnd, nt_aero_out=nt_aero) + nt_ipnd_out=nt_ipnd, nt_aero_out=nt_aero, nt_fsd_out=nt_fsd) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) @@ -897,6 +919,11 @@ subroutine init_state trcr_depend(nt_hpnd) = 2+nt_apnd ! melt pond depth trcr_depend(nt_ipnd) = 2+nt_apnd ! refrozen pond lid endif + if (tr_fsd) then + do it = 1, nfsd + trcr_depend(nt_fsd + it - 1) = 0 ! area-weighted floe size distribution + enddo + endif if (tr_aero) then ! volume-weighted aerosols do it = 1, n_aero trcr_depend(nt_aero+(it-1)*4 ) = 2 ! snow @@ -1014,7 +1041,9 @@ subroutine set_state_var (nx, & vicen, vsnon) use icedrv_arrays_column, only: hin_max - use icedrv_domain_size, only: nilyr, nslyr, max_ntrcr, ncat + use icedrv_domain_size, only: nilyr, nslyr, max_ntrcr, ncat, nfsd + use icedrv_arrays_column, only: floe_rad_c, floe_binwidth + integer (kind=int_kind), intent(in) :: & nx ! number of grid cells @@ -1067,8 +1096,8 @@ subroutine set_state_var (nx, & real (kind=dbl_kind), parameter :: & hsno_init = 0.25_dbl_kind ! initial snow thickness (m) - logical (kind=log_kind) :: tr_brine, tr_lvl - integer (kind=int_kind) :: nt_Tsfc, nt_qice, nt_qsno, nt_sice + logical (kind=log_kind) :: tr_brine, tr_lvl, tr_fsd + integer (kind=int_kind) :: nt_Tsfc, nt_qice, nt_qsno, nt_sice, nt_fsd integer (kind=int_kind) :: nt_fbri, nt_alvl, nt_vlvl character(len=*), parameter :: subname='(set_state_var)' @@ -1077,9 +1106,10 @@ subroutine set_state_var (nx, & ! query Icepack values !----------------------------------------------------------------- - call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_lvl_out=tr_lvl) + call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_lvl_out=tr_lvl, & + tr_fsd_out=tr_fsd) call icepack_query_tracer_indices( nt_Tsfc_out=nt_Tsfc, nt_qice_out=nt_qice, & - nt_qsno_out=nt_qsno, nt_sice_out=nt_sice, & + nt_qsno_out=nt_qsno, nt_sice_out=nt_sice, nt_fsd_out=nt_fsd, & nt_fbri_out=nt_fbri, nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl) call icepack_query_parameters(rhos_out=rhos, Lfresh_out=Lfresh, puny_out=puny) call icepack_warnings_flush(nu_diag) @@ -1151,6 +1181,11 @@ subroutine set_state_var (nx, & nilyr=nilyr, nslyr=nslyr, & qin=qin(:), qsn=qsn(:)) + ! floe size distribution + if (tr_fsd) call icepack_init_fsd(nfsd=nfsd, ice_ic=ice_ic, & + floe_rad_c=floe_rad_c, & + floe_binwidth=floe_binwidth, & + afsd=trcrn(i,nt_fsd:nt_fsd+nfsd-1,n)) ! surface temperature trcrn(i,nt_Tsfc,n) = Tsfc ! deg C ! ice enthalpy, salinity @@ -1206,6 +1241,11 @@ subroutine set_state_var (nx, & Tsfc = Tsfc, & nilyr=nilyr, nslyr=nslyr, & qin=qin(:), qsn=qsn(:)) + ! floe size distribution + if (tr_fsd) call icepack_init_fsd(nfsd=nfsd, ice_ic=ice_ic, & + floe_rad_c=floe_rad_c, & + floe_binwidth=floe_binwidth, & + afsd=trcrn(i,nt_fsd:nt_fsd+nfsd-1,n)) ! surface temperature trcrn(i,nt_Tsfc,n) = Tsfc ! deg C diff --git a/configuration/driver/icedrv_restart.F90 b/configuration/driver/icedrv_restart.F90 index 994fd1ff8..0110df0f9 100644 --- a/configuration/driver/icedrv_restart.F90 +++ b/configuration/driver/icedrv_restart.F90 @@ -22,6 +22,7 @@ module icedrv_restart write_restart_lvl, read_restart_lvl, & write_restart_pond_cesm, read_restart_pond_cesm, & write_restart_pond_lvl, read_restart_pond_lvl, & + write_restart_fsd, read_restart_fsd, & write_restart_aero, read_restart_aero public :: dumpfile, restartfile, final_restart, & @@ -61,7 +62,7 @@ subroutine dumpfile logical (kind=log_kind) :: & tr_iage, tr_FY, tr_lvl, tr_aero, tr_brine, & - tr_pond_topo, tr_pond_cesm, tr_pond_lvl + tr_pond_topo, tr_pond_cesm, tr_pond_lvl, tr_fsd ! solve_zsal, skl_bgc, z_tracers character(len=char_len_long) :: filename @@ -80,7 +81,7 @@ subroutine dumpfile call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, tr_brine_out=tr_brine, & tr_pond_topo_out=tr_pond_topo, tr_pond_cesm_out=tr_pond_cesm, & - tr_pond_lvl_out=tr_pond_lvl) + tr_pond_lvl_out=tr_pond_lvl,tr_fsd_out=tr_fsd) ! call icepack_query_parameters(solve_zsal_out=solve_zsal, & ! skl_bgc_out=skl_bgc, z_tracers_out=z_tracers) call icepack_warnings_flush(nu_diag) @@ -137,6 +138,7 @@ subroutine dumpfile if (tr_pond_topo) call write_restart_pond_topo() ! topographic melt ponds if (tr_aero) call write_restart_aero() ! ice aerosol if (tr_brine) call write_restart_hbrine() ! brine height + if (tr_fsd) call write_restart_fsd() ! floe size distribution ! called in icedrv_RunMod.F90 to prevent circular dependencies ! if (solve_zsal .or. skl_bgc .or. z_tracers) & ! call write_restart_bgc ! biogeochemistry @@ -174,7 +176,7 @@ subroutine restartfile (ice_ic) logical (kind=log_kind) :: & tr_iage, tr_FY, tr_lvl, tr_aero, tr_brine, & - tr_pond_topo, tr_pond_cesm, tr_pond_lvl + tr_pond_topo, tr_pond_cesm, tr_pond_lvl, tr_fsd character(len=char_len_long) :: filename character(len=*), parameter :: subname='(restartfile)' @@ -196,7 +198,7 @@ subroutine restartfile (ice_ic) call icepack_query_tracer_flags(tr_iage_out=tr_iage, tr_FY_out=tr_FY, & tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, tr_brine_out=tr_brine, & tr_pond_topo_out=tr_pond_topo, tr_pond_cesm_out=tr_pond_cesm, & - tr_pond_lvl_out=tr_pond_lvl) + tr_pond_lvl_out=tr_pond_lvl,tr_fsd_out=tr_fsd) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & file=__FILE__,line= __LINE__) @@ -261,7 +263,7 @@ subroutine restartfile (ice_ic) if (tr_pond_topo) call read_restart_pond_topo() ! topographic melt ponds if (tr_aero) call read_restart_aero() ! ice aerosol if (tr_brine) call read_restart_hbrine ! brine height - + if (tr_fsd) call read_restart_fsd() ! floe size distribution !----------------------------------------------------------------- ! Ensure ice is binned in correct categories ! (should not be necessary unless restarting from a run with @@ -491,6 +493,54 @@ end subroutine read_restart_age !======================================================================= +! Dumps all values needed for restarting +! author Elizabeth C. Hunke, LANL + + subroutine write_restart_fsd() + + use icedrv_state, only: trcrn + use icedrv_domain_size, only: ncat, nfsd + integer (kind=int_kind) :: nt_fsd, k + character(len=*), parameter :: subname='(write_restart_fsd)' + + call icepack_query_tracer_indices(nt_fsd_out=nt_fsd) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & + file=__FILE__,line= __LINE__) + + do k = 1, nfsd + call write_restart_field(nu_dump,trcrn(:,nt_fsd+k-1,:),ncat) + end do + + end subroutine write_restart_fsd + +!======================================================================= + +! Reads all values needed for a FSD restart +! author Elizabeth C. Hunke, LANL + + subroutine read_restart_fsd() + + use icedrv_state, only: trcrn + use icedrv_domain_size, only: ncat, nfsd + integer (kind=int_kind) :: nt_fsd, k + character(len=*), parameter :: subname='(read_restart_fsd)' + + write(nu_diag,*) 'min/max fsd' + + call icepack_query_tracer_indices(nt_fsd_out=nt_fsd) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & + file=__FILE__,line= __LINE__) + + do k =1, nfsd + call read_restart_field(nu_restart,trcrn(:,nt_fsd+k-1,:),ncat) + end do + + end subroutine read_restart_fsd + +!======================================================================= + ! Dumps all values needed for restarting ! author Elizabeth C. Hunke, LANL diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90 index 8b03b61c1..de4d1ce00 100644 --- a/configuration/driver/icedrv_step.F90 +++ b/configuration/driver/icedrv_step.F90 @@ -6,7 +6,7 @@ module icedrv_step - use icedrv_constants, only: c0, nu_diag + use icedrv_constants, only: c0, nu_diag, c4 use icedrv_kinds ! use icedrv_calendar, only: istep1 use icedrv_forcing, only: ocn_data_type @@ -24,7 +24,7 @@ module icedrv_step public :: step_therm1, step_therm2, step_dyn_ridge, & prep_radiation, step_radiation, ocean_mixed_layer, & - update_state, biogeochemistry + update_state, biogeochemistry, step_dyn_wave !======================================================================= @@ -104,7 +104,8 @@ subroutine step_therm1 (dt) use icedrv_arrays_column, only: fswsfcn, fswintn, fswthrun, Sswabsn, Iswabsn use icedrv_calendar, only: yday use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nx - use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fbot, Tbot, Tsnice + use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fside, & + fbot, Tbot, Tsnice use icedrv_flux, only: meltsn, melttn, meltbn, congeln, snoicen, uatm, vatm use icedrv_flux, only: wind, rhoa, potT, Qa, zlvl, strax, stray, flatn use icedrv_flux, only: fsensn, fsurfn, fcondtopn, fcondbotn @@ -271,7 +272,7 @@ subroutine step_therm1 (dt) strocnxT = strocnxT(i), strocnyT = strocnyT(i), & fbot = fbot(i), frzmlt = frzmlt(i), & Tbot = Tbot(i), Tsnice = Tsnice(i), & - rside = rside(i), & + rside = rside(i), fside = fside(i), & fsnow = fsnow(i), frain = frain(i), & fpond = fpond(i), & fsurf = fsurf(i), fsurfn = fsurfn(i,:), & @@ -333,12 +334,17 @@ end subroutine step_therm1 subroutine step_therm2 (dt) - use icedrv_arrays_column, only: hin_max, fzsal, ocean_bio + use icedrv_arrays_column, only: hin_max, fzsal, ocean_bio, & + wave_sig_ht, wave_spectrum, & + wavefreq, dwavefreq, & + floe_rad_c, floe_binwidth, & + d_afsd_latg, d_afsd_newi, d_afsd_latm, d_afsd_weld use icedrv_arrays_column, only: first_ice, bgrid, cgrid, igrid use icedrv_calendar, only: yday - use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, nltrcr, nx + use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, & + nltrcr, nx, nfsd use icedrv_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset - use icedrv_flux, only: update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside + use icedrv_flux, only: update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside use icedrv_flux, only: meltl, frazil_diag, flux_bio, faero_ocn use icedrv_init, only: tmask use icedrv_state, only: aice, aicen, aice0, trcr_depend @@ -375,6 +381,8 @@ subroutine step_therm2 (dt) do i = 1, nx if (tmask(i)) then + ! wave_sig_ht - compute here to pass to add new ice + wave_sig_ht(i) = c4*SQRT(SUM(wave_spectrum(i,:)*dwavefreq(:))) call icepack_step_therm2(dt=dt, ncat=ncat, n_aero=n_aero, & nltrcr=nltrcr, nilyr=nilyr, nslyr=nslyr, & @@ -392,7 +400,7 @@ subroutine step_therm2 (dt) n_trcr_strata=n_trcr_strata(1:ntrcr), & nt_strata=nt_strata(1:ntrcr,:), & Tf=Tf(i), sss=sss(i), & - salinz=salinz(i,:), & + salinz=salinz(i,:), fside=fside(i), & rside=rside(i), meltl=meltl(i), & frzmlt=frzmlt(i), frazil=frazil(i), & frain=frain(i), fpond=fpond(i), & @@ -406,7 +414,17 @@ subroutine step_therm2 (dt) ocean_bio=ocean_bio(i,1:nbtrcr), & frazil_diag=frazil_diag(i), & frz_onset=frz_onset(i), & - yday=yday) + yday=yday, & + nfsd=nfsd, wave_sig_ht=wave_sig_ht(i), & + wave_spectrum=wave_spectrum(i,:), & + wavefreq=wavefreq(:), & + dwavefreq=dwavefreq(:), & + d_afsd_latg=d_afsd_latg(i,:), & + d_afsd_newi=d_afsd_newi(i,:), & + d_afsd_latm=d_afsd_latm(i,:), & + d_afsd_weld=d_afsd_weld(i,:), & + floe_rad_c=floe_rad_c(:), & + floe_binwidth=floe_binwidth(:)) endif ! tmask @@ -518,6 +536,62 @@ subroutine update_state (dt, daidt, dvidt, dagedt, offset) end subroutine update_state +!======================================================================= +! +! Run one time step of wave-fracturing the floe size distribution +! +! authors: Lettie Roach, NIWA +! Elizabeth C. Hunke, LANL + + subroutine step_dyn_wave (dt) + + use icedrv_arrays_column, only: wave_spectrum, wave_sig_ht, & + d_afsd_wave, floe_rad_l, floe_rad_c, wavefreq, dwavefreq + use icedrv_domain_size, only: ncat, nfsd, nfreq, nx + use icedrv_state, only: trcrn, aicen, aice, vice + use icepack_intfc, only: icepack_step_wavefracture + + real (kind=dbl_kind), intent(in) :: & + dt ! time step + + ! local variables + + integer (kind=int_kind) :: & + i, j, & ! horizontal indices + ntrcr, & ! + nbtrcr ! + + character (len=char_len) :: wave_spec_type + + character(len=*), parameter :: subname = '(step_dyn_wave)' + + call icepack_query_parameters(wave_spec_type_out=wave_spec_type) + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & + file=__FILE__,line= __LINE__) + + do i = 1, nx + d_afsd_wave(i,:) = c0 + call icepack_step_wavefracture (wave_spec_type=wave_spec_type, & + dt=dt, ncat=ncat, nfsd=nfsd, nfreq=nfreq, & + aice = aice (i), & + vice = vice (i), & + aicen = aicen (i,:), & + floe_rad_l = floe_rad_l (:), & + floe_rad_c = floe_rad_c (:), & + wave_spectrum = wave_spectrum(i,:), & + wavefreq = wavefreq (:), & + dwavefreq = dwavefreq (:), & + trcrn = trcrn (i,:,:), & + d_afsd_wave = d_afsd_wave (i,:)) + end do ! i + + call icepack_warnings_flush(nu_diag) + if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, & + file=__FILE__,line= __LINE__) + + end subroutine step_dyn_wave + !======================================================================= ! ! Run one time step of ridging. diff --git a/configuration/scripts/icepack.build b/configuration/scripts/icepack.build index cb92e941c..715b4a30a 100755 --- a/configuration/scripts/icepack.build +++ b/configuration/scripts/icepack.build @@ -24,8 +24,8 @@ endif if !(-d ${ICE_OBJDIR}) mkdir -p ${ICE_OBJDIR} cd ${ICE_OBJDIR} -#setenv ICE_CPPDEFS "-DNXGLOB=${ICE_NXGLOB} -DNYGLOB=${ICE_NYGLOB} -DBLCKX=${ICE_BLCKX} -DBLCKY=${ICE_BLCKY} -DMXBLCKS=${ICE_MXBLCKS} -DNICELYR=${NICELYR} -DNSNWLYR=${NSNWLYR} -DNICECAT=${NICECAT} -DTRAGE=${TRAGE} -DTRFY=${TRFY} -DTRLVL=${TRLVL} -DTRPND=${TRPND} -DTRBRI=${TRBRI} -DNTRAERO=${NTRAERO} -DTRZS=${TRZS} -DNBGCLYR=${NBGCLYR} -DTRALG=${TRALG} -DTRBGCZ=${TRBGCZ} -DTRDOC=${TRDOC} -DTRDOC=${TRDOC} -DTRDIC=${TRDIC} -DTRDON=${TRDON} -DTRFED=${TRFED} -DTRFEP=${TRFEP} -DTRZAERO=${TRZAERO} -DTRBGCS=${TRBGCS} -DNUMIN=${NUMIN} -DNUMAX=${NUMAX}" -setenv ICE_CPPDEFS "-DNXGLOB=${ICE_NXGLOB} -DNICELYR=${NICELYR} -DNSNWLYR=${NSNWLYR} -DNICECAT=${NICECAT} -DTRAGE=${TRAGE} -DTRFY=${TRFY} -DTRLVL=${TRLVL} -DTRPND=${TRPND} -DTRBRI=${TRBRI} -DNTRAERO=${NTRAERO} -DTRZS=${TRZS} -DNBGCLYR=${NBGCLYR} -DTRALG=${TRALG} -DTRBGCZ=${TRBGCZ} -DTRDOC=${TRDOC} -DTRDOC=${TRDOC} -DTRDIC=${TRDIC} -DTRDON=${TRDON} -DTRFED=${TRFED} -DTRFEP=${TRFEP} -DTRZAERO=${TRZAERO} -DTRBGCS=${TRBGCS} " +#setenv ICE_CPPDEFS "-DNXGLOB=${ICE_NXGLOB} -DNYGLOB=${ICE_NYGLOB} -DBLCKX=${ICE_BLCKX} -DBLCKY=${ICE_BLCKY} -DMXBLCKS=${ICE_MXBLCKS} -DNICELYR=${NICELYR} -DNSNWLYR=${NSNWLYR} -DNICECAT=${NICECAT} -DNFSDCAT=${NFSDCAT} -DTRAGE=${TRAGE} -DTRFY=${TRFY} -DTRLVL=${TRLVL} -DTRPND=${TRPND} -DTRBRI=${TRBRI} -DNTRAERO=${NTRAERO} -DTRZS=${TRZS} -DNBGCLYR=${NBGCLYR} -DTRALG=${TRALG} -DTRBGCZ=${TRBGCZ} -DTRDOC=${TRDOC} -DTRDOC=${TRDOC} -DTRDIC=${TRDIC} -DTRDON=${TRDON} -DTRFED=${TRFED} -DTRFEP=${TRFEP} -DTRZAERO=${TRZAERO} -DTRBGCS=${TRBGCS} -DNUMIN=${NUMIN} -DNUMAX=${NUMAX}" +setenv ICE_CPPDEFS "-DNXGLOB=${ICE_NXGLOB} -DNICELYR=${NICELYR} -DNSNWLYR=${NSNWLYR} -DNICECAT=${NICECAT} -DNFSDCAT=${NFSDCAT} -DTRAGE=${TRAGE} -DTRFY=${TRFY} -DTRLVL=${TRLVL} -DTRPND=${TRPND} -DTRBRI=${TRBRI} -DNTRAERO=${NTRAERO} -DTRZS=${TRZS} -DNBGCLYR=${NBGCLYR} -DTRALG=${TRALG} -DTRBGCZ=${TRBGCZ} -DTRDOC=${TRDOC} -DTRDOC=${TRDOC} -DTRDIC=${TRDIC} -DTRDON=${TRDON} -DTRFED=${TRFED} -DTRFEP=${TRFEP} -DTRZAERO=${TRZAERO} -DTRBGCS=${TRBGCS} " if (${ICE_NTASKS} == 1) then setenv ICE_COMMDIR serial diff --git a/configuration/scripts/icepack.settings b/configuration/scripts/icepack.settings index b1de66970..e08ffdd55 100755 --- a/configuration/scripts/icepack.settings +++ b/configuration/scripts/icepack.settings @@ -37,10 +37,11 @@ setenv ICE_QUEUE undefined setenv ICE_THREADED false if (${ICE_NTHRDS} > 1) setenv ICE_THREADED true -### Layers +### Layers and Categories setenv NICELYR 7 # number of vertical layers in the ice setenv NSNWLYR 1 # number of vertical layers in the snow setenv NICECAT 5 # number of ice thickness categories +setenv NFSDCAT 1 # number of floe size categories ### Tracers # match icepack_in tracer_nml to conserve memory setenv TRAGE 1 # set to 1 for ice age tracer diff --git a/configuration/scripts/icepack_in b/configuration/scripts/icepack_in index 81072bf9e..324267a6b 100644 --- a/configuration/scripts/icepack_in +++ b/configuration/scripts/icepack_in @@ -6,7 +6,7 @@ dt = 3600.0 npt = 8760 ndtd = 1 - ice_ic = 'none' + ice_ic = 'default' restart = .false. restart_dir = './restart/' dumpfreq = 'y' @@ -28,6 +28,7 @@ tr_pond_topo = .false. tr_pond_lvl = .true. tr_aero = .false. + tr_fsd = .false. / &thermo_nml @@ -83,6 +84,7 @@ l_mpond_fresh = .false. tfrz_option = 'mushy' oceanmixed_ice = .true. + wave_spec_type = 'none' restore_ocn = .false. trestore = 90 precip_units = 'mks' diff --git a/configuration/scripts/options/set_env.fsd1 b/configuration/scripts/options/set_env.fsd1 new file mode 100644 index 000000000..8a15c8ecd --- /dev/null +++ b/configuration/scripts/options/set_env.fsd1 @@ -0,0 +1,3 @@ +setenv NICECAT 5 # number of ice thickness categories +setenv NFSDCAT 1 # number of floe size categories + diff --git a/configuration/scripts/options/set_env.fsd12 b/configuration/scripts/options/set_env.fsd12 new file mode 100644 index 000000000..5e8edb40f --- /dev/null +++ b/configuration/scripts/options/set_env.fsd12 @@ -0,0 +1,3 @@ +setenv NICECAT 5 # number of ice thickness categories +setenv NFSDCAT 12 # number of floe size categories + diff --git a/configuration/scripts/options/set_nml.fsd1 b/configuration/scripts/options/set_nml.fsd1 new file mode 100644 index 000000000..f79927e5c --- /dev/null +++ b/configuration/scripts/options/set_nml.fsd1 @@ -0,0 +1,2 @@ +tr_fsd = .true. +wave_spec_type = 'none' diff --git a/configuration/scripts/options/set_nml.fsd12 b/configuration/scripts/options/set_nml.fsd12 new file mode 100644 index 000000000..f3534b8ce --- /dev/null +++ b/configuration/scripts/options/set_nml.fsd12 @@ -0,0 +1,2 @@ +tr_fsd = .true. +wave_spec_type = 'constant' diff --git a/configuration/scripts/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts index 990fdbe5c..14cfd3019 100644 --- a/configuration/scripts/tests/base_suite.ts +++ b/configuration/scripts/tests/base_suite.ts @@ -12,6 +12,8 @@ smoke col 1x1 debug,run1year,alt03 smoke col 1x1 debug,run1year,alt04 smoke col 1x1 debug,run1year,leap,dt30min smoke col 1x1 debug,run1year,dyn +smoke col 1x1 debug,run1year,fsd12 +smoke col 1x1 debug,run1year,fsd1 restart col 1x1 debug restart col 1x1 diag1 restart col 1x1 pondcesm @@ -25,3 +27,5 @@ restart col 1x1 alt02 restart col 1x1 alt03 restart col 1x1 alt04 restart col 1x1 dyn +restart col 1x1 fsd12 + diff --git a/configuration/scripts/tests/report_results.csh b/configuration/scripts/tests/report_results.csh index 7cf65228a..68ea6179d 100755 --- a/configuration/scripts/tests/report_results.csh +++ b/configuration/scripts/tests/report_results.csh @@ -169,7 +169,6 @@ if ( $fbuild != "" || $frun != "" || $ftest != "" ) then if (${tchkpass} == 1) then @ tpass = $tpass + 1 - else else if (${tchkpass} == 2) then @ tunkn = $tunkn + 1 diff --git a/doc/source/master_list.bib b/doc/source/master_list.bib index b22a03aff..534956e3f 100644 --- a/doc/source/master_list.bib +++ b/doc/source/master_list.bib @@ -597,6 +597,76 @@ @Article{Duarte17 pages = {1632-1654}, url = {http://dx.doi.org/10.1002/2016JG003660} } + +@article{Horvat15, + author = "C. Horvat and E. Tziperman", + journal = {The Cryosphere}, + number = {6}, + pages = {2119-2134}, + title = "{A prognostic model of the sea-ice floe size and thickness distribution}", + url = {http://dx.doi.org/10.5194/tc-9-2119-2015}, + volume = {9}, + year = {2015} +} + +@article{Roach18, + author = "L. A. Roach and C. Horvat and S. M. Dean and C. M. Bitz", + url = {http://dx.doi.org/10.1029/2017JC013692}, + journal = JGRO, + number = {6}, + pages = {4322--4337}, + title = "{An emergent sea ice floe size distribution in a global coupled ocean-sea ice model}", + volume = {123}, + year = {2018} +} + +@article{Rothrock84, + author = "D. A. Rothrock and A. S. Thorndike", + journal = JGR, + number = {C4}, + pages = {6477}, + title = {{Measuring the sea ice floe size distribution}}, + url = {http://doi.wiley.com/10.1029/JC089iC04p06477}, + volume = {89}, + year = {1984} +} + +@article{Perovich14, +author = "D. K. Perovich and K. F. Jones", +journal = JGRO, +number = {12}, +pages = {8767--8777}, +title = {{The seasonal evolution of sea ice floe size distribution}}, +url = {http://doi.wiley.com/10.1002/2014JC010136}, +volume = {119}, +year = {2014} +} + +@article{Roach18b, +author = "L. A. Roach and M. M. Smith and S. M. Dean", +url = {http://doi.wiley.com/10.1002/2017JC013693}, +journal = JGRO, +number = {4}, +pages = {2851--2866}, +title = {{Quantifying growth of pancake sea ice floes using images from drifting buoys}}, +volume = {123}, +year = {2018} +} + +@article{Shen01, +author = "H. H. Shen and S. F. Ackley and M. A. Hopkins", +url = {http://dx.doi.org/10.3189/172756401781818239}, +journal = {Annals of Glaciology}, +number = {2}, +pages = {361--367}, +title = {{A conceptual model for pancake-ice formation in a wave field}}, +volume = {33}, +year = {2001} +} + + + + % ********************************************** % For new entries, see example entry in BIB_TEMPLATE.txt % ********************************************** diff --git a/doc/source/science_guide/sg_itd.rst b/doc/source/science_guide/sg_itd.rst index 4b22aaf04..a382204d7 100755 --- a/doc/source/science_guide/sg_itd.rst +++ b/doc/source/science_guide/sg_itd.rst @@ -90,4 +90,89 @@ In the WMO case, the distribution used depends on the number of categories used. +----------------+------------+---------+--------+--------+--------+ | 7 | | | | | 2.00 | +----------------+------------+---------+--------+--------+--------+ + +Joint floe size and thickness distribution +============================ + +Sizes of individual sea ice floes vary over an extremely broad range, from centimeters +to hundreds of kilometers. The floe size distribution (FSD) is a probability function that +characterizes this variability :cite:`Rothrock84`. An option to include a +prognostic sea ice floe size distribution is available and used if ``tr_fsd`` is set to true. +The scheme is based on the theoretical framework described in :cite:`Horvat15` for a +*joint* floe size and thickness distribution (FSTD), and was implemented by :cite:`Roach18`. + +In this theory, individual floes are identified with a size :math:`r` and area :math:`x(r)`, where +:math:`x(r)=4\alpha r^2` for :math:`\alpha=0.66 < \pi/4` (:cite:`Rothrock84`). The probability +distribution :math:`f(r,h) dr dh` is the fraction of grid surface area +covered by ice with thickness between :math:`h` and :math:`h + dh` and lateral floe +size between :math:`r` and :math:`r + dr`. The FSTD integrates over all floe sizes and +ice thicknesses to unity; over all floe sizes to the ITD; and over all thicknesses to the FSD. + +For implementation in CICE, the continuous function :math:`f(r,h)` is replaced +with a product of two discrete variables: :math:`a_{in}` as defined above and :math:`F_{in,k}`. +:math:`F_{in,k}` is the fraction of ice belonging to thickness category :math:`n` with lateral +floe size belonging to floe size class :math:`k`, giving +:math:`\sum_{n=0}^{N_C}\sum_{k=0}^{N_f} a_{in} F_{in,k} = 1` and :math:`\sum_{k=0}^{N_f} F_{in,k} = 1`. +:math:`F_{in,k}` is carried as an area-weighted tracer. + +The FSD may be ignored when considering processes that only modify ice thickness +(eg. vertical thermodynamics), and the ITD can be ignored when considering processes that only modify floe sizes (eg. wave fracture). For processes that affect both the ITD and the FSD, (eg. lateral melt), +both :math:`a_{in}` and :math:`F_{in,k}` are evolved. + +The FSTD evolves subject to lateral growth, lateral melt, new ice growth, floe welding and +wave fracture, as described in :cite:`Roach18` and with some modifications described in +:cite:`Roach19`. The equation for time evolution of the FSTD is (:cite:`Horvat15`), + +:math:`\frac{\partial f(r,h)}{\partial t} = - \nabla \cdot (f(r,h)\mathbf{v}) + \mathcal{L}_T + \mathcal{L}_M + \mathcal{L}_W`, + +where the terms on the right hand side represent the effects of advection, thermodynamics, mechanical +redistribution and wave fracture respectively. Floe sizes do not explicitly appear in the equations of sea ice motion and therefore the FSTD is advected as an area tracer. We also assume that mechanical redistribution of sea ice through ridging does not impact floe sizes. Thus it remains only to compute the thermodynamic and wave fracture tendencies. + +Thermodynamic changes to the FSTD are given by + +:math:`\mathcal{L}_T(r,h)=-\nabla_{(r,h)} \cdot (f(r,h) \mathbf{G}) +\frac{2}{r}f(r,h)G_r + +\delta(r-r_{\text{min}})\delta(h-h_{\text{min}})\dot{A}_p + \beta_{\text{weld}}.` + +The first two terms on the right-hand side represent growth and melt of existing floes +in thickness and lateral size, at a rate :math:`\mathbf{G} = (G_r,G_h)`. The third +term represents growth of new ice: new floes are created at a rate :math:`\dot{A}_p` +in the smallest thickness category and a given lateral size category. If wave forcing +is provided, the size of newly formed floes is determined via a tensile stress limitation +arising from the wave field (:cite:`Shen2001`,:cite:`Roach2019`); otherwise, all floes +are presumed to grow as pancakes in the smallest floe size category resolved. +To allow for the joining of individual floes to one another, we represent +the welding together of floes in freezing conditions via the fourth term, +:math:`\beta_{\text{weld}}`, using a coagulation equation. + +To compute the impact of wave fracture of the FSD, given a local ocean surface wave +spectrum is provided, we generate a realization of the sea surface height field, which +is uniquely determined by the spectrum up to a phase. In :cite:`Horvat2015` this phase is +randomly chosen, and multiple realizations of the resulting surface height field are used to +obtain convergent statistics. However this stochastic component would lead to a model that is +not bit-for-bit reproducible. Users can choose in the namelist (via ``wave_spec_type``) +to run the model with the phase set to be constant to obtain bit-for-bit reproducibility or +to include the random phase. + +We calculate the number and length of fractures that would occur if waves enter a fully ice-covered +region defined in one dimension in the direction of propagation, and then apply +the outcome proportionally to the ice-covered fraction in each grid cell. +Assuming that sea ice flexes with the sea surface height field, strains are computed +on this sub-grid-scale 1D domain. If the strain between successive extrema exceeds +a critical value new floes are formed with diameters equal to the distance between +the extrema. + +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`, the number of floe size categories. It is assumed that +the floe size lies at the midpoint of each floe size category. Note that there is a sensitivity to the choice +of floe size categories inherent in the current floe welding scheme. + +If simulations begin without ice (``ice_init='none'``), the FSD can emerge without initialization. This +is the recommended initialization for studies on the FSD itself. If simulations begin with ice cover, +some initial FSD must be prescribed in ``init_fsd``. The default (used for ``ice_init='default'``) +is a simple relationship determined from point observations by :cite:`Perovich14`, but its basin-wide +applicability has not been tested. In Icepack, ``ice_init='default'`` is selected for the slab +and the full ITD cells. + + + diff --git a/doc/source/science_guide/sg_thermo.rst b/doc/source/science_guide/sg_thermo.rst index 74a74399d..c927e1512 100755 --- a/doc/source/science_guide/sg_thermo.rst +++ b/doc/source/science_guide/sg_thermo.rst @@ -1716,6 +1716,34 @@ to a specified minimum thickness; if the open water area is nearly zero or if there is more new ice than will fit into the thinnest ice category, then the new ice is spread over the entire cell. +If ``tr_fsd=true``, a floe size must be assigned to the new frazil ice. +If spectral ocean surface wave forcing is provided (and set using the +namelist option ``wave_spec_type``), this will be used +to calculate a tensile stress on new floes that determines their maximum +possible size :cite:`Shen2001,Roach2019`. If no ocean surface wave forcing +is provided, all floes are assumed to grow as pancakes, at the smallest +possible floe size. + +If ``tr_fsd=true``, lateral growth at the edges of exisiting floes may +also occur, calculated using the prognostic floe size distribution as +described in :cite:`Horvat15` and :cite:`Roach18`. The lateral growth +that occurs is a portion of the total new ice growth, depending on the +area of open water close to floe edges. Lateral growth +modifies the ITD and the FSD. + +If ``tr_fsd=true``, floes may weld together thermodynamically during +freezing conditions according to the probability that they overlap, +assuming they are replaced randomly on the domain. Evolution of the +FSD is described using a coagulation equation. +The total number of floes that weld with another, per square meter, +per unit time, in the case of a fully covered ice surface was estimated +from observations in :cite:`Roach18b`. In its original model implementation, +with 12 floe size categories, the tendency term for floe welding was divided by a +constant equal to the area of the largest floe, (approx 2 km^2), with this choice made +as the product of sensitivity studies to balance the climatological tendencies of +wave fracture and welding. So that results do not vary as the number or range of +floe size categories varies, we fix this scaling coefficient, c_weld. + If the latent heat flux is negative (i.e., latent heat is transferred from the ice to the atmosphere), snow or snow-free ice sublimates at the top surface. If the latent heat flux is positive, vapor from the @@ -1749,11 +1777,16 @@ old and new layers, respectively. The enthalpies of the new layers are .. math:: q_k = \frac{1}{\Delta h_i} \sum_{m=1}^{N_i} \eta_{km} q_m. -Lateral melting is accomplished by multiplying the state variables by +If ``tr_fsd=false``, lateral melting is accomplished by multiplying the state variables by :math:`1-r_{side}`, where :math:`r_{side}` is the fraction of ice melted laterally :cite:`Maykut87,Steele92`, and adjusting the ice energy and fluxes as appropriate. We assume a floe diameter of 300 m. +If ``tr_fsd=true``, lateral melting is accomplished using the :cite:`Maykut87` +lateral heat flux, but applied to the ice using the prognostic floe size distribution +as described in :cite:`Horvat15` and :cite:`Roach18`. Lateral melt modifies +the ITD and the FSD. + Snow-ice formation ------------------ @@ -1805,3 +1838,5 @@ into equal thicknesses while conserving energy and salt. calculations, unlike Bitz99, which does include it. This extra heat is returned to the ocean (or the atmosphere, in the case of evaporation) with the melt water. + + diff --git a/doc/source/science_guide/sg_tracers.rst b/doc/source/science_guide/sg_tracers.rst index c7af8f9ac..c884f5b78 100755 --- a/doc/source/science_guide/sg_tracers.rst +++ b/doc/source/science_guide/sg_tracers.rst @@ -10,7 +10,7 @@ required (surface temperature and thickness, salinity and enthalpy of ice and sn and many others are options. For instance, there are tracers to track the age of the ice; the area of first-year ice, fractions of ice area and volume that are level, from which the amount of deformed ice can be calculated; pond area, pond volume and volume of ice covering ponds; -aerosols and numerous other biogeochemical tracers. +a prognostic floe size distribution; aerosols and numerous other biogeochemical tracers. Most of these tracers are presented in later sections. Here we describe the ice age tracers and how tracers may depend on other tracers, using the pond tracers as an example. diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst index 790745ac0..e6aacf506 100755 --- a/doc/source/user_guide/ug_case_settings.rst +++ b/doc/source/user_guide/ug_case_settings.rst @@ -56,6 +56,7 @@ to support the CICE model. "NICELYR", " ", "number of vertical layers in the ice", "7" "NSNWLYR", " ", "number of vertical layers in the snow", "1" "NICECAT", " ", "number of ice thickness categories", "5" + "NFSDCAT", " ", "number of floe size categories", "1" "TRAGE", "0,1", "ice age tracer", "1" "TRFY", "0,1", "first-year ice area tracer", "1" "TRLVL", "0,1", "deformed ice tracer", "1" @@ -227,6 +228,9 @@ column physics. "``l_mpond_fresh``", "true", "retain (topo) pond water until ponds drain", "" "", "false", "release (topo) pond water immediately to ocean", "" "``oceanmixed_ice``", "true/false", "active ocean mixed layer calculation", "``.true.`` (if uncoupled)" + "``wave_spec_type``", "``none``", "no ocean wave spectrum data - no wave-ice interations","" + "", "``constant``", "ocean wave spectrum data present, sea surface height field generated using constant phase, for testing FSD", "" + "", "``random``", "ocean wave spectrum data present, sea surface height field generated using random phase", "" "``ocn_data_type``", "``default``", "constant values defined in the code", "" "", "``ISPOL``", "ISPOL experiment data (see :ref:`force`)", "" "", "``NICE``", "N-ICE experiment data (see :ref:`force`)", "" diff --git a/icepack.setup b/icepack.setup index 2d3aed74a..5cc0d7da9 100755 --- a/icepack.setup +++ b/icepack.setup @@ -741,7 +741,7 @@ EOF2 cat ${ICE_SCRIPTS}/options/set_env.${name} >> ${fsmods} -cat >> ${fimods} << EOF2 +cat >> ${fsmods} << EOF2 EOF2 echo "adding env mods set_env.${name}"