From 286630fdecc8efc4dcf4bd7ecf2315101c69c391 Mon Sep 17 00:00:00 2001 From: "David A. Bailey" Date: Wed, 30 Oct 2024 18:29:04 -0600 Subject: [PATCH] FSD fixes for conservation (#495) This is a bug fix plus some rearranging for conservation in CICE. This impacts FSD cases as well as non-FSD cases. There are two main aspects to the bug fix: Initial ice and snow volume values at the beginning of the lateral melt routine are now used for updating the snow and ice enthalpy as well other tracers for lateral melting. This is answer changing in all cases! The other fix is to move the computation of vi0new_lat (lateral growth from FSD) outside of the subroutine fsd_lateral_growth and a check to make sure the category vi0new (vin0new) does not exceed the total vi0new. This is a warning and it recomputes to ensure conservation. This is only answer changing in FSD cases. There is also a change in definition of floe_area_c, so that it is precisely half way between floe_area_l and floe_area_h. The original computation of floe_area_c based on floe_rad_c biases the area towards the lower value. Update the fsd subroutine arguments (floe_rad_c, floe_rad_l, floe_binwidth, c_fsd_range) to store static data inside Icepack when computed in Icepack rather than pass them back and forth. Provide an ability to pass out the values from Icepack as optional arguments. Move the computation of rsiden and get rid of fside which is no longer needed. Clean up some comments and indentation Update swccsm3 test, set ktherm=1 Update interface documentation There is an associated CICE tag to go with this. --- columnphysics/icepack_fsd.F90 | 140 +++--- columnphysics/icepack_therm_itd.F90 | 434 +++++++----------- columnphysics/icepack_therm_vertical.F90 | 122 ++++- columnphysics/icepack_wavefracspec.F90 | 19 +- configuration/driver/icedrv_InitMod.F90 | 13 +- configuration/driver/icedrv_arrays_column.F90 | 17 +- configuration/driver/icedrv_flux.F90 | 3 +- configuration/driver/icedrv_forcing.F90 | 1 + configuration/driver/icedrv_init.F90 | 5 - configuration/driver/icedrv_step.F90 | 34 +- configuration/scripts/options/set_nml.swccsm3 | 1 + doc/source/user_guide/interfaces.include | 66 +-- 12 files changed, 408 insertions(+), 447 deletions(-) diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90 index 28f12d922..2e5175afa 100644 --- a/columnphysics/icepack_fsd.F90 +++ b/columnphysics/icepack_fsd.F90 @@ -55,13 +55,19 @@ module icepack_fsd public :: icepack_init_fsd_bounds, icepack_init_fsd, icepack_cleanup_fsd, & fsd_lateral_growth, fsd_add_new_ice, fsd_weld_thermo, get_subdt_fsd - real(kind=dbl_kind), dimension(:), allocatable :: & + real(kind=dbl_kind), dimension(:), allocatable, public :: & floe_rad_h, & ! fsd size higher bound in m (radius) + floe_rad_c, & ! fsd size center in m (radius) + floe_rad_l, & ! fsd size lower bound in m (radius) + floe_binwidth, & ! fsd size binwidth in m (radius) floe_area_l, & ! fsd area at lower bound (m^2) floe_area_h, & ! fsd area at higher bound (m^2) floe_area_c, & ! fsd area at bin centre (m^2) floe_area_binwidth ! floe area bin width (m^2) + character (len=35), dimension(:), allocatable, public :: & + c_fsd_range ! string for history output + integer(kind=int_kind), dimension(:,:), allocatable, public :: & iweld ! floe size categories that can combine ! during welding (dimensionless) @@ -85,22 +91,22 @@ module icepack_fsd ! authors: Lettie Roach, NIWA/VUW and C. M. Bitz, UW subroutine icepack_init_fsd_bounds( & - floe_rad_l, & ! fsd size lower bound in m (radius) - floe_rad_c, & ! fsd size bin centre in m (radius) - floe_binwidth, & ! fsd size bin width in m (radius) - c_fsd_range, & ! string for history output - write_diags ) ! flag for writing diagnostics + floe_rad_l_out, & ! fsd size lower bound in m (radius) + floe_rad_c_out, & ! fsd size bin centre in m (radius) + floe_binwidth_out, & ! fsd size bin width in m (radius) + c_fsd_range_out, & ! string for history output + write_diags) ! flag for writing diagnostics - real(kind=dbl_kind), dimension(:), intent(inout) :: & - floe_rad_l, & ! fsd size lower bound in m (radius) - floe_rad_c, & ! fsd size bin centre in m (radius) - floe_binwidth ! fsd size bin width in m (radius) + real(kind=dbl_kind), dimension(:), intent(out), optional :: & + floe_rad_l_out, & ! fsd size lower bound in m (radius) + floe_rad_c_out, & ! fsd size bin centre in m (radius) + floe_binwidth_out ! fsd size bin width in m (radius) - character (len=35), intent(out) :: & - c_fsd_range(nfsd) ! string for history output + character (len=35), dimension(:), intent(out), optional :: & + c_fsd_range_out ! string for history output logical (kind=log_kind), intent(in), optional :: & - write_diags ! write diags flag + write_diags ! write diags flag !autodocument_end @@ -111,11 +117,8 @@ subroutine icepack_init_fsd_bounds( & real (kind=dbl_kind) :: test - real (kind=dbl_kind), dimension (0:nfsd) :: & - floe_rad - real (kind=dbl_kind), dimension(:), allocatable :: & - lims + lims, floe_rad character(len=8) :: c_fsd1,c_fsd2 character(len=2) :: c_nf @@ -169,10 +172,15 @@ subroutine icepack_init_fsd_bounds( & allocate( & floe_rad_h (nfsd), & ! fsd size higher bound in m (radius) + floe_rad_l (nfsd), & ! fsd size lower bound in m (radius) + floe_rad_c (nfsd), & ! fsd size center in m (radius) + floe_rad (0:nfsd), & ! fsd bounds in m (radius) floe_area_l (nfsd), & ! fsd area at lower bound (m^2) floe_area_h (nfsd), & ! fsd area at higher bound (m^2) floe_area_c (nfsd), & ! fsd area at bin centre (m^2) floe_area_binwidth (nfsd), & ! floe area bin width (m^2) + floe_binwidth (nfsd), & ! floe bin width (m) + c_fsd_range (nfsd), & ! iweld (nfsd, nfsd), & ! fsd categories that can weld stat=ierr) if (ierr/=0) then @@ -186,8 +194,11 @@ subroutine icepack_init_fsd_bounds( & floe_rad_c = (floe_rad_h+floe_rad_l)/c2 floe_area_l = c4*floeshape*floe_rad_l**2 - floe_area_c = c4*floeshape*floe_rad_c**2 floe_area_h = c4*floeshape*floe_rad_h**2 +! floe_area_c = c4*floeshape*floe_rad_c**2 +! This is exactly in the middle of floe_area_h and floe_area_l +! Whereas the above calculation is closer to floe_area_l. + floe_area_c = (floe_area_h+floe_area_l)/c2 floe_binwidth = floe_rad_h - floe_rad_l @@ -220,20 +231,56 @@ subroutine icepack_init_fsd_bounds( & enddo if (present(write_diags)) then - if (write_diags) then - write(warnstr,*) ' ' - call icepack_warnings_add(warnstr) - write(warnstr,*) subname - call icepack_warnings_add(warnstr) - write(warnstr,*) 'floe_rad(n-1) < fsd Cat n < floe_rad(n)' - call icepack_warnings_add(warnstr) - do n = 1, nfsd - write(warnstr,*) floe_rad(n-1),' < fsd Cat ',n, ' < ',floe_rad(n) + if (write_diags) then + write(warnstr,*) ' ' call icepack_warnings_add(warnstr) - enddo - write(warnstr,*) ' ' - call icepack_warnings_add(warnstr) + write(warnstr,*) subname + call icepack_warnings_add(warnstr) + write(warnstr,*) 'floe_rad(n-1) < fsd Cat n < floe_rad(n)' + call icepack_warnings_add(warnstr) + do n = 1, nfsd + write(warnstr,*) floe_rad(n-1),' < fsd Cat ',n, ' < ',floe_rad(n) + call icepack_warnings_add(warnstr) + enddo + write(warnstr,*) ' ' + call icepack_warnings_add(warnstr) + endif + endif + + if (present(floe_rad_l_out)) then + if (size(floe_rad_l_out) /= size(floe_rad_l)) then + call icepack_warnings_add(subname//' floe_rad_l_out incorrect size') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + floe_rad_l_out(:) = floe_rad_l(:) + endif + + if (present(floe_rad_c_out)) then + if (size(floe_rad_c_out) /= size(floe_rad_c)) then + call icepack_warnings_add(subname//' floe_rad_c_out incorrect size') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + floe_rad_c_out(:) = floe_rad_c(:) endif + + if (present(floe_binwidth_out)) then + if (size(floe_binwidth_out) /= size(floe_binwidth)) then + call icepack_warnings_add(subname//' floe_binwidth_out incorrect size') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + floe_binwidth_out(:) = floe_binwidth(:) + endif + + if (present(c_fsd_range_out)) then + if (size(c_fsd_range_out) /= size(c_fsd_range)) then + call icepack_warnings_add(subname//' c_fsd_range_out incorrect size') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + c_fsd_range_out(:) = c_fsd_range(:) endif end subroutine icepack_init_fsd_bounds @@ -256,18 +303,11 @@ end subroutine icepack_init_fsd_bounds ! ! authors: Lettie Roach, NIWA/VUW - subroutine icepack_init_fsd(ice_ic, & - floe_rad_c, & ! fsd size bin centre in m (radius) - floe_binwidth, & ! fsd size bin width in m (radius) - afsd) ! floe size distribution tracer + subroutine icepack_init_fsd(ice_ic, afsd) ! floe size distribution tracer character(len=char_len_long), intent(in) :: & ice_ic ! method of ice cover initialization - real(kind=dbl_kind), dimension(:), intent(inout) :: & - floe_rad_c, & ! fsd size bin centre in m (radius) - floe_binwidth ! fsd size bin width in m (radius) - real (kind=dbl_kind), dimension (:), intent(inout) :: & afsd ! floe size tracer: fraction distribution of floes @@ -323,7 +363,6 @@ subroutine icepack_cleanup_fsd (afsdn) character(len=*), parameter :: subname='(icepack_cleanup_fsd)' - if (tr_fsd) then do n = 1, ncat @@ -378,14 +417,11 @@ end subroutine icepack_cleanup_fsdn ! ! authors: Lettie Roach, NIWA/VUW - subroutine partition_area (floe_rad_c, aice, & + subroutine partition_area (aice, & aicen, vicen, & afsdn, lead_area, & latsurf_area) - real (kind=dbl_kind), dimension(:), intent(in) :: & - floe_rad_c ! fsd size bin centre in m (radius) - real (kind=dbl_kind), intent(in) :: & aice ! ice concentration @@ -476,7 +512,7 @@ end subroutine partition_area subroutine fsd_lateral_growth (dt, aice, & aicen, vicen, & vi0new, & - frazil, floe_rad_c, & + frazil, & afsdn, & lead_area, latsurf_area, & G_radial, d_an_latg, & @@ -497,10 +533,6 @@ subroutine fsd_lateral_growth (dt, aice, & vi0new , & ! volume of new ice added to cat 1 (m) frazil ! frazil ice growth (m/step-->cm/day) - ! floe size distribution - real (kind=dbl_kind), dimension (:), intent(in) :: & - floe_rad_c ! fsd size bin centre in m (radius) - real (kind=dbl_kind), dimension(ncat), intent(out) :: & d_an_latg ! change in aicen occuring due to lateral growth @@ -529,7 +561,7 @@ subroutine fsd_lateral_growth (dt, aice, & d_an_latg = c0 ! partition volume into lateral growth and frazil - call partition_area (floe_rad_c, aice, & + call partition_area (aice, & aicen, vicen, & afsdn, lead_area, & latsurf_area) @@ -540,9 +572,6 @@ subroutine fsd_lateral_growth (dt, aice, & vi0new_lat = vi0new * lead_area / (c1 + aice/latsurf_area) end if - ! for history/diagnostics - frazil = vi0new - vi0new_lat - ! lateral growth increment if (vi0new_lat > puny) then G_radial = vi0new_lat/dt @@ -563,8 +592,6 @@ subroutine fsd_lateral_growth (dt, aice, & endif ! vi0new_lat > 0 ! Use remaining ice volume as in standard model, - ! but ice cannot grow into the area that has grown laterally - vi0new = vi0new - vi0new_lat tot_latg = SUM(d_an_latg(:)) end subroutine fsd_lateral_growth @@ -589,7 +616,6 @@ end subroutine fsd_lateral_growth subroutine fsd_add_new_ice (n, & dt, ai0new, & d_an_latg, d_an_newi, & - floe_rad_c, floe_binwidth, & G_radial, area2, & wave_sig_ht, & wave_spectrum, & @@ -623,9 +649,7 @@ subroutine fsd_add_new_ice (n, & real (kind=dbl_kind), dimension (:), intent(in) :: & aicen_init , & ! fractional area of ice - aicen , & ! after update - floe_rad_c , & ! fsd size bin centre in m (radius) - floe_binwidth ! fsd size bin width in m (radius) + aicen ! after update real (kind=dbl_kind), dimension (:,:), intent(in) :: & afsdn ! floe size distribution tracer diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90 index 02a9ab4b8..0fcc96738 100644 --- a/columnphysics/icepack_therm_itd.F90 +++ b/columnphysics/icepack_therm_itd.F90 @@ -20,6 +20,9 @@ module icepack_therm_itd use icepack_kinds + + use icepack_fsd, only: floe_rad_c, floe_binwidth + use icepack_parameters, only: c0, c1, c2, c3, c4, c6, c10 use icepack_parameters, only: p001, p1, p333, p5, p666, puny, bignum use icepack_parameters, only: rhos, rhoi, Lfresh, ice_ref_salinity @@ -870,12 +873,11 @@ subroutine lateral_melt (dt, fpond, & fresh, fsalt, & fhocn, faero_ocn, & fiso_ocn, & - rside, meltl, & - fside, wlat, & + rsiden, meltl, & + wlat, & aicen, vicen, & vsnon, trcrn, & - flux_bio, d_afsd_latm,& - floe_rad_c, floe_binwidth) + flux_bio, d_afsd_latm) real (kind=dbl_kind), intent(in) :: & dt ! time step (s) @@ -888,15 +890,12 @@ subroutine lateral_melt (dt, fpond, & real (kind=dbl_kind), dimension (:,:), intent(inout) :: & trcrn ! tracer array - real (kind=dbl_kind), intent(in) :: & - rside ! fraction of ice that melts laterally + real (kind=dbl_kind), dimension (:), intent(in) :: & + rsiden ! fraction of ice that melts laterally real (kind=dbl_kind), intent(in), optional :: & wlat ! lateral melt rate (m/s) - real (kind=dbl_kind), intent(inout) :: & - fside ! lateral heat flux (W/m^2) - real (kind=dbl_kind), intent(inout) :: & fpond , & ! fresh water flux to ponds (kg/m^2/s) fresh , & ! fresh water flux to ocean (kg/m^2/s) @@ -913,10 +912,6 @@ subroutine lateral_melt (dt, fpond, & real (kind=dbl_kind), dimension(:), intent(inout), optional :: & fiso_ocn ! isotope flux to ocean (kg/m^2/s) - real (kind=dbl_kind), dimension (:), intent(in), optional :: & - floe_rad_c , & ! fsd size bin centre in m (radius) - floe_binwidth ! fsd size bin width in m (radius) - real (kind=dbl_kind), dimension (:), intent(out), optional :: & d_afsd_latm ! change in fsd due to lateral melt (m) @@ -934,19 +929,13 @@ subroutine lateral_melt (dt, fpond, & dfsalt , & ! change in fsalt dvssl , & ! snow surface layer volume dvint , & ! snow interior layer - bin1_arealoss, tmp ! - - logical (kind=log_kind) :: & - fsd_wlat, & ! .true. if wlat present and wlat > puny - flag ! .true. if there could be lateral melting + tmp real (kind=dbl_kind), dimension (ncat) :: & aicen_init, & ! initial area fraction vicen_init, & ! volume per unit area of ice (m) - vsnon_init, & ! initial volume of snow (m) - G_radialn , & ! rate of lateral melt (m/s) - delta_an , & ! change in the ITD - rsiden ! delta_an/aicen + vsnon_init, & ! volume per unit area of snow (m) + G_radialn ! rate of lateral melt (m/s) real (kind=dbl_kind), dimension (:,:), allocatable :: & afsdn , & ! floe size distribution tracer @@ -966,19 +955,18 @@ subroutine lateral_melt (dt, fpond, & character(len=*), parameter :: subname='(lateral_melt)' - flag = .false. + if (tr_fsd) d_afsd_latm = c0 + + if (.not. any(rsiden(:) > c0)) return ! no lateral melt, get out now + dfhocn = c0 dfpond = c0 dfresh = c0 dfsalt = c0 dvssl = c0 dvint = c0 - bin1_arealoss = c0 - tmp = c0 - vicen_init = vicen(:) - G_radialn = c0 - delta_an = c0 - rsiden = c0 + vicen_init(:) = vicen(:) + vsnon_init(:) = vsnon(:) if (tr_fsd) then call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:)) @@ -995,216 +983,141 @@ subroutine lateral_melt (dt, fpond, & afsdn = trcrn(nt_fsd:nt_fsd+nfsd-1,:) afsdn_init = afsdn ! for diagnostics df_flx = c0 - d_afsd_latm = c0 f_flx = c0 end if - ! fsd_wlat == if (tr_fsd .and. wlat > puny) - ! need fsd_wlat because wlat is optional - fsd_wlat = .false. - if (tr_fsd .and. present(wlat)) then - if (wlat > puny) fsd_wlat = .true. - endif - - if (fsd_wlat) then - flag = .true. - - ! for FSD rside and fside not yet computed correctly, redo here - fside = c0 - do n = 1, ncat - - G_radialn(n) = -wlat ! negative - - if (any(afsdn(:,n) < c0)) then - write(warnstr,*) subname, 'lateral_melt B afsd < 0 ',n - call icepack_warnings_add(warnstr) - endif - - bin1_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) - bin1_arealoss - - if (delta_an(n) > c0) then - write(warnstr,*) subname, 'ERROR delta_an > 0 ',delta_an(n) - call icepack_warnings_add(warnstr) - endif - - ! following original code, not necessary for fsd - if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n)/aicen(n),c1) - - if (rsiden(n) < c0) then - write(warnstr,*) subname, 'ERROR rsiden < 0 ',rsiden(n) - call icepack_warnings_add(warnstr) - endif - - ! melting energy/unit area in each column, etot < 0 - etot = c0 - do k = 1, nslyr - etot = etot + trcrn(nt_qsno+k-1,n) * vsnon(n)/real(nslyr,kind=dbl_kind) - enddo - - do k = 1, nilyr - etot = etot + trcrn(nt_qice+k-1,n) * vicen(n)/real(nilyr,kind=dbl_kind) - enddo ! nilyr - - ! lateral heat flux, fside < 0 - fside = fside + rsiden(n)*etot/dt - - 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. - - do n = 1, ncat + do n = 1, ncat !----------------------------------------------------------------- ! Melt the ice and increment fluxes. !----------------------------------------------------------------- - ! fluxes to coupler - ! dfresh > 0, dfsalt > 0, dfpond > 0 + ! fluxes to coupler + ! dfresh > 0, dfsalt > 0, dfpond > 0 - dfresh = (rhoi*vicen(n) + rhos*vsnon(n)) * rsiden(n) / dt - if (saltflux_option == 'prognostic') then - sicen = c0 - do k=1,nilyr - sicen = sicen + trcrn(nt_sice+k-1,n) / real(nilyr,kind=dbl_kind) - enddo - dfsalt = rhoi*vicen(n)*sicen*p001 * rsiden(n) / dt - else - dfsalt = rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt - endif - fresh = fresh + dfresh - fsalt = fsalt + dfsalt - - if (tr_pond_topo) then - dfpond = aicen(n)*trcrn(nt_apnd,n)*trcrn(nt_hpnd,n)*rsiden(n) - fpond = fpond - dfpond - endif - - ! history diagnostics - meltl = meltl + vicen(n)*rsiden(n) + dfresh = (rhoi*vicen(n) + rhos*vsnon(n)) * rsiden(n) / dt + if (saltflux_option == 'prognostic') then + sicen = c0 + do k=1,nilyr + sicen = sicen + trcrn(nt_sice+k-1,n) / real(nilyr,kind=dbl_kind) + enddo + dfsalt = rhoi*vicen(n)*sicen*p001 * rsiden(n) / dt + else + dfsalt = rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt + endif + fresh = fresh + dfresh + fsalt = fsalt + dfsalt - ! state variables - vicen_init(n) = vicen(n) - vsnon_init(n) = vsnon(n) - aicen(n) = aicen(n) * (c1 - rsiden(n)) - vicen(n) = vicen(n) * (c1 - rsiden(n)) - vsnon(n) = vsnon(n) * (c1 - rsiden(n)) + if (tr_pond_topo) then + dfpond = aicen(n)*trcrn(nt_apnd,n)*trcrn(nt_hpnd,n)*rsiden(n) + fpond = fpond - dfpond + endif - ! floe size distribution - if (tr_fsd) then - if (rsiden(n) > puny) then - if (aicen(n) > puny) then + ! history diagnostics + meltl = meltl + vicen_init(n)*rsiden(n) - ! adaptive subtimestep - elapsed_t = c0 - afsd_tmp(:) = afsdn_init(:,n) - d_afsd_tmp(:) = c0 - nsubt = 0 + ! state variables + aicen(n) = aicen(n) * (c1 - rsiden(n)) + vicen(n) = vicen(n) * (c1 - rsiden(n)) + vsnon(n) = vsnon(n) * (c1 - rsiden(n)) - DO WHILE (elapsed_t.lt.dt) + ! floe size distribution + if (tr_fsd) then + if (rsiden(n) > puny) then + if (aicen(n) > puny) then - nsubt = nsubt + 1 - if (nsubt.gt.100) then - write(warnstr,*) subname, 'latm not converging' - call icepack_warnings_add(warnstr) - endif + ! adaptive subtimestep + elapsed_t = c0 + afsd_tmp(:) = afsdn_init(:,n) + d_afsd_tmp(:) = c0 + nsubt = 0 - ! finite differences - df_flx(:) = c0 - f_flx (:) = c0 - do k = 2, nfsd - f_flx(k) = G_radialn(n) * afsd_tmp(k) / floe_binwidth(k) - end do + DO WHILE (elapsed_t.lt.dt) - do k = 1, nfsd - df_flx(k) = f_flx(k+1) - f_flx(k) - end do + nsubt = nsubt + 1 + if (nsubt.gt.100) then + write(warnstr,*) subname, 'latm not converging' + call icepack_warnings_add(warnstr) + endif - if (abs(sum(df_flx(:))) > puny) then - write(warnstr,*) subname, 'sum(df_flx) /= 0' - call icepack_warnings_add(warnstr) - endif + ! finite differences + df_flx(:) = c0 + f_flx (:) = c0 + G_radialn(n) = -wlat + do k = 2, nfsd + f_flx(k) = G_radialn(n) * afsd_tmp(k) / floe_binwidth(k) + end do - ! this term ensures area conservation - tmp = SUM(afsd_tmp(:)/floe_rad_c(:)) + do k = 1, nfsd + df_flx(k) = f_flx(k+1) - f_flx(k) + end do - ! fsd tendency - do k = 1, nfsd - d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsd_tmp(k) & - * (c1/floe_rad_c(k) - tmp) - end do + if (abs(sum(df_flx(:))) > puny) then + write(warnstr,*) subname, 'sum(df_flx) /= 0' + call icepack_warnings_add(warnstr) + endif - ! timestep required for this - subdt = get_subdt_fsd(afsd_tmp(:), d_afsd_tmp(:)) - subdt = MIN(subdt, dt) + ! this term ensures area conservation + tmp = SUM(afsd_tmp(:)/floe_rad_c(:)) - ! update fsd and elapsed time - afsd_tmp(:) = afsd_tmp(:) + subdt*d_afsd_tmp(:) - elapsed_t = elapsed_t + subdt + ! fsd tendency + do k = 1, nfsd + d_afsd_tmp(k) = -df_flx(k) + c2 * G_radialn(n) * afsd_tmp(k) & + * (c1/floe_rad_c(k) - tmp) + end do + ! timestep required for this + subdt = get_subdt_fsd(afsd_tmp(:), d_afsd_tmp(:)) + subdt = MIN(subdt, dt) - END DO + ! update fsd and elapsed time + afsd_tmp(:) = afsd_tmp(:) + subdt*d_afsd_tmp(:) + elapsed_t = elapsed_t + subdt - afsdn(:,n) = afsd_tmp(:) + END DO + afsdn(:,n) = afsd_tmp(:) - end if ! aicen - end if ! rside > 0, otherwise do nothing + end if ! aicen + end if ! rside > 0, otherwise do nothing - end if ! tr_fsd + 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)*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)*rsiden(n) / dt & - * vsnon(n)/real(nslyr,kind=dbl_kind) - fhocn = fhocn + dfhocn - enddo ! nslyr + ! 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)*rsiden(n) / dt & + * vicen_init(n)/real(nilyr,kind=dbl_kind) + fhocn = fhocn + dfhocn + enddo ! nilyr - if (tr_aero) then - do k = 1, n_aero - faero_ocn(k) = faero_ocn(k) + (vsnon(n) & - *(trcrn(nt_aero +4*(k-1),n) & - + trcrn(nt_aero+1+4*(k-1),n)) & - + vicen(n) & - *(trcrn(nt_aero+2+4*(k-1),n) & - + trcrn(nt_aero+3+4*(k-1),n))) & - * rsiden(n) / dt - enddo ! k - endif ! tr_aero - - if (tr_iso) then - do k = 1, n_iso - fiso_ocn(k) = fiso_ocn(k) & - + (vsnon(n)*trcrn(nt_isosno+k-1,n) & - + vicen(n)*trcrn(nt_isoice+k-1,n)) & - * rside / dt - enddo ! k - endif ! tr_iso + do k = 1, nslyr + ! heat flux to coupler for snow melt (dfhocn < 0) + dfhocn = trcrn(nt_qsno+k-1,n)*rsiden(n) / dt & + * vsnon_init(n)/real(nslyr,kind=dbl_kind) + fhocn = fhocn + dfhocn + enddo ! nslyr + + if (tr_aero) then + do k = 1, n_aero + faero_ocn(k) = faero_ocn(k) & + + (vsnon_init(n) * (trcrn(nt_aero +4*(k-1),n) & + + trcrn(nt_aero+1+4*(k-1),n)) & + + vicen_init(n) * (trcrn(nt_aero+2+4*(k-1),n) & + + trcrn(nt_aero+3+4*(k-1),n))) & + * rsiden(n) / dt + enddo ! k + endif ! tr_aero + + if (tr_iso) then + do k = 1, n_iso + fiso_ocn(k) = fiso_ocn(k) & + + (vsnon_init(n)*trcrn(nt_isosno+k-1,n) & + + vicen_init(n)*trcrn(nt_isoice+k-1,n)) & + * rsiden(n) / dt + enddo ! k + endif ! tr_iso !----------------------------------------------------------------- ! Biogeochemistry @@ -1229,32 +1142,31 @@ subroutine lateral_melt (dt, fpond, & trcrn, flux_bio) if (icepack_warnings_aborted(subname)) return - endif ! flag - if (tr_fsd) then + if (tr_fsd) then - trcrn(nt_fsd:nt_fsd+nfsd-1,:) = afsdn + trcrn(nt_fsd:nt_fsd+nfsd-1,:) = afsdn - call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:) ) - if (icepack_warnings_aborted(subname)) return + call icepack_cleanup_fsd (trcrn(nt_fsd:nt_fsd+nfsd-1,:) ) + if (icepack_warnings_aborted(subname)) return - ! 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 + ! 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 - deallocate(afsdn) - deallocate(afsdn_init) - deallocate(df_flx) - deallocate(afsd_tmp) - deallocate(d_afsd_tmp) - deallocate(f_flx) + deallocate(afsdn) + deallocate(afsdn_init) + deallocate(df_flx) + deallocate(afsd_tmp) + deallocate(d_afsd_tmp) + deallocate(f_flx) - end if + end if end subroutine lateral_melt @@ -1301,8 +1213,7 @@ subroutine add_new_ice (dt, & wavefreq, & dwavefreq, & d_afsd_latg, & - d_afsd_newi, & - floe_rad_c, floe_binwidth) + d_afsd_newi) use icepack_fsd, only: fsd_lateral_growth, fsd_add_new_ice @@ -1376,10 +1287,6 @@ subroutine add_new_ice (dt, & wavefreq, & ! wave frequencies (s^-1) dwavefreq ! wave frequency bin widths (s^-1) - real (kind=dbl_kind), dimension (:), intent(in), optional :: & - floe_rad_c , & ! fsd size bin centre in m (radius) - floe_binwidth ! fsd size bin width in m (radius) - real (kind=dbl_kind), dimension(:), intent(out), optional :: & ! change in thickness distribution (area) d_afsd_latg , & ! due to fsd lateral growth @@ -1608,11 +1515,30 @@ subroutine add_new_ice (dt, & call fsd_lateral_growth(dt, aice, & aicen, vicen, & vi0new, frazil, & - floe_rad_c, afsdn, & + afsdn, & lead_area, latsurf_area, & G_radial, d_an_latg, & tot_latg) if (icepack_warnings_aborted(subname)) return + + ! volume added to each category from lateral growth + do n = 1, ncat + if (aicen(n) > puny) then + vin0new(n) = d_an_latg(n) * vicen(n)/aicen(n) + endif + end do + + ! check lateral growth doesn't exceed total growth + ! if it does, adjust it + if (SUM(vin0new)>vi0new) then + write(warnstr,*) subname,'lateral growth warning ',vi0new,SUM(vin0new) + call icepack_warnings_add(warnstr) + vin0new(:) = vin0new(:)*vi0new/SUM(vin0new) + end if + + vi0new = vi0new - SUM(vin0new) + frazil = frazil - SUM(vin0new) + endif ai0mod = aice0 @@ -1639,11 +1565,6 @@ subroutine add_new_ice (dt, & vi0new = c0 endif ! aice0 > puny - ! volume added to each category from lateral growth - 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) @@ -1790,7 +1711,6 @@ subroutine add_new_ice (dt, & call fsd_add_new_ice (n, & dt, ai0new, & d_an_latg, d_an_newi, & - floe_rad_c, floe_binwidth, & G_radial, area2, & wave_sig_ht, & wave_spectrum, & @@ -1933,8 +1853,8 @@ subroutine icepack_step_therm2(dt, hin_max, & nt_strata, & Tf, sss, & salinz, & - rside, meltl, & - fside, wlat, & + rsiden, meltl, & + wlat, & frzmlt, frazil, & frain, fpond, & fresh, fsalt, & @@ -1951,8 +1871,7 @@ subroutine icepack_step_therm2(dt, hin_max, & wavefreq, & dwavefreq, & d_afsd_latg, d_afsd_newi, & - d_afsd_latm, d_afsd_weld, & - floe_rad_c, floe_binwidth) + d_afsd_latm, d_afsd_weld) use icepack_parameters, only: icepack_init_parameters @@ -1966,7 +1885,6 @@ subroutine icepack_step_therm2(dt, hin_max, & dt , & ! time step Tf , & ! freezing temperature (C) sss , & ! sea surface salinity (ppt) - rside , & ! fraction of ice that melts laterally frzmlt ! freezing/melting potential (W/m^2) integer (kind=int_kind), dimension (:), intent(in) :: & @@ -1981,13 +1899,13 @@ subroutine icepack_step_therm2(dt, hin_max, & nt_strata ! indices of underlying tracer layers real (kind=dbl_kind), dimension(:), intent(in) :: & + rsiden , & ! fraction of ice that melts laterally salinz , & ! initial salinity profile ocean_bio ! ocean concentration of biological tracer real (kind=dbl_kind), intent(inout) :: & aice , & ! sea ice concentration aice0 , & ! concentration of open water - fside , & ! lateral heat flux (W/m^2) frain , & ! rainfall rate (kg/m^2 s) fpond , & ! fresh water flux to ponds (kg/m^2/s) fresh , & ! fresh water flux to ocean (kg/m^2/s) @@ -2050,10 +1968,6 @@ subroutine icepack_step_therm2(dt, hin_max, & d_afsd_latm, & ! lateral melt d_afsd_weld ! welding - real (kind=dbl_kind), dimension (:), intent(in), optional :: & - floe_rad_c, & ! fsd size bin centre in m (radius) - floe_binwidth ! fsd size bin width in m (radius) - !autodocument_end ! local variables @@ -2090,9 +2004,7 @@ subroutine icepack_step_therm2(dt, hin_max, & present(d_afsd_latg) .and. & present(d_afsd_newi) .and. & present(d_afsd_latm) .and. & - present(d_afsd_weld) .and. & - present(floe_rad_c) .and. & - present(floe_binwidth))) then + present(d_afsd_weld))) then call icepack_warnings_add(subname//' error in FSD arguments, tr_fsd=T') call icepack_warnings_setabort(.true.,__FILE__,__LINE__) return @@ -2177,8 +2089,7 @@ subroutine icepack_step_therm2(dt, hin_max, & wave_sig_ht, & wave_spectrum, & wavefreq, dwavefreq, & - d_afsd_latg, d_afsd_newi, & - floe_rad_c, floe_binwidth) + d_afsd_latg, d_afsd_newi) if (icepack_warnings_aborted(subname)) return @@ -2190,13 +2101,12 @@ subroutine icepack_step_therm2(dt, hin_max, & fresh, fsalt, & fhocn, faero_ocn, & fiso_ocn, & - rside, meltl, & - fside, wlat, & + rsiden, meltl, & + wlat, & aicen, vicen, & vsnon, trcrn, & flux_bio, & - d_afsd_latm, & - floe_rad_c,floe_binwidth) + d_afsd_latm) if (icepack_warnings_aborted(subname)) return ! Floe welding during freezing conditions diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90 index 6dda193af..bf96096eb 100644 --- a/columnphysics/icepack_therm_vertical.F90 +++ b/columnphysics/icepack_therm_vertical.F90 @@ -20,7 +20,10 @@ module icepack_therm_vertical use icepack_kinds - use icepack_parameters, only: c0, c1, p001, p5, puny + + use icepack_fsd, only: floe_rad_c, floe_binwidth + + use icepack_parameters, only: c0, c1, c2, p001, p5, puny use icepack_parameters, only: pi, depressT, Lvap, hs_min, cp_ice, min_salin use icepack_parameters, only: cp_ocn, rhow, rhoi, rhos, Lfresh, rhofresh, ice_ref_salinity use icepack_parameters, only: ktherm, calc_Tsfc, rsnw_fall, rsnw_tmax @@ -30,7 +33,7 @@ module icepack_therm_vertical use icepack_parameters, only: saltflux_option, congel_freeze use icepack_parameters, only: icepack_chkoptargflag - use icepack_tracers, only: ncat, nilyr, nslyr + use icepack_tracers, only: ncat, nilyr, nslyr, nfsd use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd, tr_iso use icepack_tracers, only: tr_pond_lvl, tr_pond_topo use icepack_tracers, only: n_aero, n_iso @@ -484,8 +487,9 @@ subroutine frzmlt_bottom_lateral (dt, & fbot_xfer_type, & strocnxT, strocnyT, & Tbot, fbot, & - rside, Cdn_ocn, & - fside, wlat) + rsiden, Cdn_ocn, & + wlat, aicen, & + afsdn) real (kind=dbl_kind), intent(in) :: & dt ! time step @@ -513,14 +517,28 @@ subroutine frzmlt_bottom_lateral (dt, & 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 - fside ! lateral heat flux (W/m^2) + fbot ! heat flux to ice bottom (W/m^2) + + real (kind=dbl_kind), dimension(:), intent(out) :: & + rsiden ! fraction of ice that melts laterally real (kind=dbl_kind), intent(out), optional :: & wlat ! lateral melt rate (m/s) + real (kind=dbl_kind), dimension(:), intent(in) :: & + aicen ! ice concentration + + real (kind=dbl_kind), dimension (:,:), intent(in), optional :: & + afsdn ! area floe size distribution + ! local variables + real (kind=dbl_kind), dimension (ncat) :: & + delta_an , & ! change in the ITD + G_radialn ! lateral melt rate in FSD (m/s) + + real (kind=dbl_kind) :: & + rside , & ! + fside ! lateral heat flux (W/m^2) integer (kind=int_kind) :: & n , & ! thickness category index @@ -534,6 +552,7 @@ subroutine frzmlt_bottom_lateral (dt, & real (kind=dbl_kind) :: & deltaT , & ! SST - Tbot >= 0 ustar , & ! skin friction velocity for fbot (m/s) + bin1_arealoss, & xtmp ! temporary variable ! Parameters for bottom melting @@ -555,11 +574,11 @@ subroutine frzmlt_bottom_lateral (dt, & ! Identify grid cells where ice can melt. !----------------------------------------------------------------- - rside = c0 - fside = c0 + rsiden(:) = c0 Tbot = Tf fbot = c0 wlat_loc = c0 + if (present(wlat)) wlat=c0 if (aice > puny .and. frzmlt < c0) then ! ice can melt @@ -599,10 +618,52 @@ subroutine frzmlt_bottom_lateral (dt, & rside = wlat_loc*dt*pi/(floeshape*floediam) ! Steele rside = max(c0,min(rside,c1)) + if (rside == c0) return ! nothing more to do so get out + + rsiden(:) = rside + + if (tr_fsd) then ! alter rsiden now since floes are not of size floediam + + do n = 1, ncat + G_radialn(n) = -wlat_loc ! negative + + ! afsdn present check up the calling tree + if (any(afsdn(:,n) < c0)) then + write(warnstr,*) subname, 'lateral_melt B afsd < 0 ',n + call icepack_warnings_add(warnstr) + endif + + bin1_arealoss = -afsdn(1,n) / floe_binwidth(1) ! when scaled by *G_radialn(n)*dt*aicen(n) + + delta_an(n) = c0 + do k = 1, nfsd + ! this is delta_an(n) when scaled by *G_radialn(n)*dt*aicen(n) + delta_an(n) = delta_an(n) + ((c2/floe_rad_c(k)) * afsdn(k,n)) ! delta_an < 0 + end do + + ! add negative area loss from fsd + delta_an(n) = (delta_an(n) - bin1_arealoss)*G_radialn(n)*dt + + if (delta_an(n) > c0) then + write(warnstr,*) subname, 'ERROR delta_an > 0 ',delta_an(n) + call icepack_warnings_add(warnstr) + endif + + ! following original code, not necessary for fsd + if (aicen(n) > c0) rsiden(n) = MIN(-delta_an(n),c1) + + if (rsiden(n) < c0) then + write(warnstr,*) subname, 'ERROR rsiden < 0 ',rsiden(n) + call icepack_warnings_add(warnstr) + endif + enddo ! ncat + + endif ! if tr_fsd + !----------------------------------------------------------------- ! Compute heat flux associated with this value of rside. !----------------------------------------------------------------- - + fside = c0 do n = 1, ncat etot = c0 @@ -617,20 +678,23 @@ subroutine frzmlt_bottom_lateral (dt, & enddo ! nilyr ! lateral heat flux, fside < 0 - fside = fside + rside*etot/dt + fside = fside + rsiden(n)*etot/dt enddo ! n !----------------------------------------------------------------- ! Limit bottom and lateral heat fluxes if necessary. - ! FYI: fside is not yet correct for fsd, may need to adjust fbot further + ! Limit rside so we don't melt laterally more ice than frzmlt permits !----------------------------------------------------------------- xtmp = frzmlt/(fbot + fside - puny) xtmp = min(xtmp, c1) + xtmp = max(xtmp, c0) fbot = fbot * xtmp - rside = rside * xtmp - fside = fside * xtmp + + do n = 1, ncat + rsiden(n) = rsiden(n) * xtmp ! xtmp is almost always 1 so usually nothing happens here + enddo ! ncat endif @@ -2089,8 +2153,8 @@ subroutine icepack_step_therm1(dt, & strocnxT , strocnyT , & fbot , & Tbot , Tsnice , & - frzmlt , rside , & - fside , wlat , & + frzmlt , rsiden , & + wlat , & fsnow , frain , & fpond , fsloss , & fsurf , fsurfn , & @@ -2137,7 +2201,7 @@ subroutine icepack_step_therm1(dt, & lmask_n , lmask_s , & mlt_onset , frz_onset , & yday , prescribed_ice, & - zlvs) + zlvs , afsdn) real (kind=dbl_kind), intent(in) :: & dt , & ! time step @@ -2213,8 +2277,6 @@ subroutine icepack_step_therm1(dt, & strocnyT , & ! ice-ocean stress, y-direction 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) @@ -2227,7 +2289,7 @@ subroutine icepack_step_therm1(dt, & frz_onset ! day of year that freezing begins (congel or frazil) real (kind=dbl_kind), intent(out), optional :: & - wlat ! lateral melt rate (m/s) + wlat ! lateral melt rate (m/s) real (kind=dbl_kind), intent(inout), optional :: & fswthru_vdr , & ! vis dir shortwave penetrating to ocean (W/m^2) @@ -2261,6 +2323,9 @@ subroutine icepack_step_therm1(dt, & H2_18O_ocn , & ! ocean concentration of H2_18O (kg/kg) zlvs ! atm level height for scalars (if different than zlvl) (m) + real (kind=dbl_kind), dimension(:,:), intent(in), optional :: & + afsdn ! afsd tracer + real (kind=dbl_kind), dimension(:), intent(inout) :: & aicen_init , & ! fractional area of ice vicen_init , & ! volume per unit area of ice (m) @@ -2276,6 +2341,7 @@ subroutine icepack_step_therm1(dt, & ipnd , & ! melt pond refrozen lid thickness (m) iage , & ! volume-weighted ice age FY , & ! area-weighted first-year ice area + rsiden , & ! fraction of ice that melts laterally fsurfn , & ! net flux to top surface, excluding fcondtop fcondtopn , & ! downward cond flux at top surface (W m-2) fcondbotn , & ! downward cond flux at bottom surface (W m-2) @@ -2420,8 +2486,13 @@ subroutine icepack_step_therm1(dt, & return endif if (tr_fsd) then - if (.not.(present(wlat))) then - call icepack_warnings_add(subname//' error in FSD arguments, tr_fsd=T') + if (.not.present(afsdn)) then + call icepack_warnings_add(subname//' error missing afsdn argument, tr_fsd=T') + call icepack_warnings_setabort(.true.,__FILE__,__LINE__) + return + endif + if (size(afsdn,dim=1) /= nfsd .or. size(afsdn,dim=2) /= ncat) then + call icepack_warnings_add(subname//' error size of afsdn argument, tr_fsd=T') call icepack_warnings_setabort(.true.,__FILE__,__LINE__) return endif @@ -2492,7 +2563,7 @@ subroutine icepack_step_therm1(dt, & endif !----------------------------------------------------------------- - ! Adjust frzmlt to account for ice-ocean heat fluxes since last + ! Use frzmlt to account for ice-ocean heat fluxes since last ! call to coupler. ! Compute lateral and bottom heat fluxes. !----------------------------------------------------------------- @@ -2506,8 +2577,9 @@ subroutine icepack_step_therm1(dt, & fbot_xfer_type, & strocnxT, strocnyT, & Tbot, fbot, & - rside, Cdn_ocn, & - fside, wlat) + rsiden, Cdn_ocn, & + wlat, aicen, & + afsdn) if (icepack_warnings_aborted(subname)) return diff --git a/columnphysics/icepack_wavefracspec.F90 b/columnphysics/icepack_wavefracspec.F90 index 1beb7a234..34c4e448f 100644 --- a/columnphysics/icepack_wavefracspec.F90 +++ b/columnphysics/icepack_wavefracspec.F90 @@ -182,7 +182,6 @@ end function get_dafsd_wave subroutine icepack_step_wavefracture(wave_spec_type, & dt, nfreq, & aice, vice, aicen, & - floe_rad_l, floe_rad_c, & wave_spectrum, wavefreq, dwavefreq, & trcrn, d_afsd_wave) @@ -201,10 +200,6 @@ subroutine icepack_step_wavefracture(wave_spec_type, & 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) @@ -243,6 +238,9 @@ subroutine icepack_step_wavefracture(wave_spec_type, & afsd_tmp , & ! tracer array d_afsd_tmp ! change + real (kind=dbl_kind) :: & + local_sig_ht + character(len=*),parameter :: & subname='(icepack_step_wavefracture)' @@ -256,15 +254,15 @@ subroutine icepack_step_wavefracture(wave_spec_type, & ! if all ice is not in first floe size category if (.NOT. ALL(trcrn(nt_fsd,:).ge.c1-puny)) then - + local_sig_ht = c4*SQRT(SUM(wave_spectrum(:)*dwavefreq(:))) ! do not try to fracture for minimal ice concentration or zero wave spectrum - if ((aice > p01).and.(MAXVAL(wave_spectrum(:)) > puny)) then +! if ((aice > p01).and.(MAXVAL(wave_spectrum(:)) > puny)) then + if ((aice > p01).and.(local_sig_ht>0.1_dbl_kind)) then hbar = vice / aice ! calculate fracture histogram call wave_frac(nfreq, wave_spec_type, & - floe_rad_l, floe_rad_c, & wavefreq, dwavefreq, & hbar, wave_spectrum, fracture_hist) @@ -393,7 +391,6 @@ end subroutine icepack_step_wavefracture ! authors: 2018 Lettie Roach, NIWA/VUW subroutine wave_frac(nfreq, wave_spec_type, & - floe_rad_l, floe_rad_c, & wavefreq, dwavefreq, & hbar, spec_efreq, frac_local) @@ -406,10 +403,6 @@ subroutine wave_frac(nfreq, wave_spec_type, & 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) diff --git a/configuration/driver/icedrv_InitMod.F90 b/configuration/driver/icedrv_InitMod.F90 index c45d28022..ba1b1de6f 100644 --- a/configuration/driver/icedrv_InitMod.F90 +++ b/configuration/driver/icedrv_InitMod.F90 @@ -31,16 +31,14 @@ 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_arrays_column, only: hin_max, c_hi_range, floe_rad_c 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_init_snow use icepack_intfc, only: icepack_warnings_flush - use icedrv_domain_size, only: ncat, nfsd + use icedrv_domain_size, only: ncat ! use icedrv_diagnostics, only: icedrv_diagnostics_debug use icedrv_flux, only: init_coupler_flux, init_history_therm, & init_flux_atm_ocn @@ -98,12 +96,7 @@ subroutine icedrv_initialize endif if (tr_fsd) then - call icepack_init_fsd_bounds( & - 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 - write_diags=.true.) + call icepack_init_fsd_bounds(floe_rad_c_out=floe_rad_c, write_diags=.true. ) call icepack_warnings_flush(nu_diag) if (icepack_warnings_aborted(subname)) then call icedrv_system_abort(file=__FILE__,line=__LINE__) diff --git a/configuration/driver/icedrv_arrays_column.F90 b/configuration/driver/icedrv_arrays_column.F90 index b55f129f4..dbe1ab56a 100644 --- a/configuration/driver/icedrv_arrays_column.F90 +++ b/configuration/driver/icedrv_arrays_column.F90 @@ -221,29 +221,22 @@ module icedrv_arrays_column ! 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) + floe_rad_c ! fsd size bin centre in m (radius) real (kind=dbl_kind), dimension (nx), public :: & - wave_sig_ht ! significant height of waves (m) + wave_sig_ht ! significant height of waves (m) real (kind=dbl_kind), dimension (nfreq), public :: & - wavefreq, & ! wave frequencies - dwavefreq ! wave frequency bin widths + wavefreq, & ! wave frequencies + dwavefreq ! wave frequency bin widths real (kind=dbl_kind), dimension (nx,nfreq), public :: & - wave_spectrum ! wave spectrum + 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_flux.F90 b/configuration/driver/icedrv_flux.F90 index 39c090818..321a84912 100644 --- a/configuration/driver/icedrv_flux.F90 +++ b/configuration/driver/icedrv_flux.F90 @@ -248,6 +248,7 @@ module icedrv_flux real (kind=dbl_kind), & dimension (nx,ncat), public :: & + rsiden, & ! fraction of ice that melts laterally fsurfn, & ! category fsurf fcondtopn,& ! category fcondtop fcondbotn,& ! category fcondbot @@ -272,8 +273,6 @@ module icedrv_flux !----------------------------------------------------------------- real (kind=dbl_kind), dimension (nx), public :: & - rside , & ! fraction of ice that melts laterally - fside , & ! lateral heat flux (W/m^2) wlat , & ! lateral melt rate (m/s) fsw , & ! incoming shortwave radiation (W/m^2) coszen , & ! cosine solar zenith angle, < 0 for sun below horizon diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90 index b3f360f7a..98e390e22 100644 --- a/configuration/driver/icedrv_forcing.F90 +++ b/configuration/driver/icedrv_forcing.F90 @@ -1150,6 +1150,7 @@ subroutine get_wave_spec ! 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) diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90 index ddcd8b2f6..f342d74fe 100644 --- a/configuration/driver/icedrv_init.F90 +++ b/configuration/driver/icedrv_init.F90 @@ -1296,7 +1296,6 @@ subroutine set_state_var (nx, & use icedrv_arrays_column, only: hin_max 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 @@ -1437,8 +1436,6 @@ subroutine set_state_var (nx, & ! floe size distribution if (tr_fsd) call icepack_init_fsd(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 @@ -1507,8 +1504,6 @@ subroutine set_state_var (nx, & qin=qin(:), qsn=qsn(:)) ! floe size distribution if (tr_fsd) call icepack_init_fsd(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 diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90 index f99d47e4c..13f6f2f83 100644 --- a/configuration/driver/icedrv_step.F90 +++ b/configuration/driver/icedrv_step.F90 @@ -112,8 +112,8 @@ subroutine step_therm1 (dt) use icedrv_arrays_column, only: fswthrun, fswthrun_vdr, fswthrun_vdf, fswthrun_idr, fswthrun_idf use icedrv_arrays_column, only: meltsliqn, meltsliq use icedrv_calendar, only: yday - use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, n_iso, nx - use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rside, fside, wlat, & + use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, n_iso, nfsd, nx + use icedrv_flux, only: frzmlt, sst, Tf, strocnxT, strocnyT, rsiden, wlat, & fbot, Tbot, Tsnice use icedrv_flux, only: meltsn, melttn, meltbn, congeln, snoicen, uatm, vatm use icedrv_flux, only: wind, rhoa, potT, Qa, Qa_iso, zlvl, strax, stray, flatn @@ -151,7 +151,7 @@ subroutine step_therm1 (dt) integer (kind=int_kind) :: & ntrcr, nt_apnd, nt_hpnd, nt_ipnd, nt_alvl, nt_vlvl, nt_Tsfc, & - nt_iage, nt_FY, nt_qice, nt_sice, nt_qsno, & + nt_iage, nt_FY, nt_qice, nt_sice, nt_qsno, nt_fsd, & nt_aero, nt_isosno, nt_isoice, nt_rsnw, nt_smice, nt_smliq logical (kind=log_kind) :: & @@ -202,7 +202,7 @@ subroutine step_therm1 (dt) nt_qice_out=nt_qice, nt_sice_out=nt_sice, & nt_aero_out=nt_aero, nt_qsno_out=nt_qsno, & nt_rsnw_out=nt_rsnw, nt_smice_out=nt_smice, nt_smliq_out=nt_smliq, & - nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice) + nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice, 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__) @@ -325,7 +325,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), fside = fside(i), & + rsiden = rsiden(i,:), & wlat = wlat(i), & fsnow = fsnow(i), frain = frain(i), & fpond = fpond(i), fsloss = fsloss(i), & @@ -370,6 +370,7 @@ subroutine step_therm1 (dt) snoice = snoice(i), snoicen = snoicen(i,:), & dsnow = dsnow(i), dsnown = dsnown(i,:), & meltsliqn= meltsliqn(i,:), & + afsdn = trcrn (i,nt_fsd:nt_fsd+nfsd-1,:), & lmask_n = lmask_n(i), lmask_s = lmask_s(i), & mlt_onset=mlt_onset(i), frz_onset = frz_onset(i), & yday = yday, prescribed_ice = prescribed_ice) @@ -429,14 +430,13 @@ subroutine step_therm2 (dt) use icedrv_arrays_column, only: hin_max, 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 use icedrv_calendar, only: yday use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, & - nx, nfsd + nx use icedrv_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset - use icedrv_flux, only: fsalt, Tf, sss, salinz, fhocn, rside, fside, wlat + use icedrv_flux, only: fsalt, Tf, sss, salinz, fhocn, rsiden, wlat use icedrv_flux, only: meltl, frazil_diag, flux_bio, faero_ocn, fiso_ocn use icedrv_flux, only: HDO_ocn, H2_16O_ocn, H2_18O_ocn use icedrv_init, only: tmask @@ -497,9 +497,9 @@ 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,:), fside=fside(i), & + salinz=salinz(i,:), & wlat=wlat(i), & - rside=rside(i), meltl=meltl(i), & + rsiden=rsiden(i,:), meltl=meltl(i), & frzmlt=frzmlt(i), frazil=frazil(i), & frain=frain(i), fpond=fpond(i), & fresh=fresh(i), fsalt=fsalt(i), & @@ -522,9 +522,7 @@ subroutine step_therm2 (dt) 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(:)) + d_afsd_weld=d_afsd_weld(i,:)) endif ! tmask @@ -654,8 +652,8 @@ end subroutine update_state 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 + d_afsd_wave, wavefreq, dwavefreq + use icedrv_domain_size, only: ncat, nfreq, nx use icedrv_state, only: trcrn, aicen, aice, vice use icepack_intfc, only: icepack_step_wavefracture @@ -685,11 +683,9 @@ subroutine step_dyn_wave (dt) 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 (:), & + wavefreq = wavefreq (:), & + dwavefreq = dwavefreq (:), & trcrn = trcrn (i,:,:), & d_afsd_wave = d_afsd_wave (i,:)) end do ! i diff --git a/configuration/scripts/options/set_nml.swccsm3 b/configuration/scripts/options/set_nml.swccsm3 index ae4d64d96..4ec03fddb 100644 --- a/configuration/scripts/options/set_nml.swccsm3 +++ b/configuration/scripts/options/set_nml.swccsm3 @@ -1,3 +1,4 @@ +ktherm = 1 shortwave = 'ccsm3' albedo_type = 'ccsm3' calc_tsfc = .true. diff --git a/doc/source/user_guide/interfaces.include b/doc/source/user_guide/interfaces.include index ae640d1bf..b41f6055d 100644 --- a/doc/source/user_guide/interfaces.include +++ b/doc/source/user_guide/interfaces.include @@ -142,22 +142,22 @@ icepack_init_fsd_bounds ! authors: Lettie Roach, NIWA/VUW and C. M. Bitz, UW subroutine icepack_init_fsd_bounds( & - floe_rad_l, & ! fsd size lower bound in m (radius) - floe_rad_c, & ! fsd size bin centre in m (radius) - floe_binwidth, & ! fsd size bin width in m (radius) - c_fsd_range, & ! string for history output - write_diags ) ! flag for writing diagnostics + floe_rad_l_out, & ! fsd size lower bound in m (radius) + floe_rad_c_out, & ! fsd size bin centre in m (radius) + floe_binwidth_out, & ! fsd size bin width in m (radius) + c_fsd_range_out, & ! string for history output + write_diags) ! flag for writing diagnostics - real(kind=dbl_kind), dimension(:), intent(inout) :: & - floe_rad_l, & ! fsd size lower bound in m (radius) - floe_rad_c, & ! fsd size bin centre in m (radius) - floe_binwidth ! fsd size bin width in m (radius) + real(kind=dbl_kind), dimension(:), intent(out), optional :: & + floe_rad_l_out, & ! fsd size lower bound in m (radius) + floe_rad_c_out, & ! fsd size bin centre in m (radius) + floe_binwidth_out ! fsd size bin width in m (radius) - character (len=35), intent(out) :: & - c_fsd_range(nfsd) ! string for history output + character (len=35), dimension(:), intent(out), optional :: & + c_fsd_range_out ! string for history output logical (kind=log_kind), intent(in), optional :: & - write_diags ! write diags flag + write_diags ! write diags flag @@ -172,18 +172,11 @@ icepack_init_fsd ! ! authors: Lettie Roach, NIWA/VUW - subroutine icepack_init_fsd(ice_ic, & - floe_rad_c, & ! fsd size bin centre in m (radius) - floe_binwidth, & ! fsd size bin width in m (radius) - afsd) ! floe size distribution tracer + subroutine icepack_init_fsd(ice_ic, afsd) ! floe size distribution tracer character(len=char_len_long), intent(in) :: & ice_ic ! method of ice cover initialization - real(kind=dbl_kind), dimension(:), intent(inout) :: & - floe_rad_c, & ! fsd size bin centre in m (radius) - floe_binwidth ! fsd size bin width in m (radius) - real (kind=dbl_kind), dimension (:), intent(inout) :: & afsd ! floe size tracer: fraction distribution of floes @@ -2253,8 +2246,8 @@ icepack_step_therm2 nt_strata, & Tf, sss, & salinz, & - rside, meltl, & - fside, wlat, & + rsiden, meltl, & + wlat, & frzmlt, frazil, & frain, fpond, & fresh, fsalt, & @@ -2271,8 +2264,7 @@ icepack_step_therm2 wavefreq, & dwavefreq, & d_afsd_latg, d_afsd_newi, & - d_afsd_latm, d_afsd_weld, & - floe_rad_c, floe_binwidth) + d_afsd_latm, d_afsd_weld) use icepack_parameters, only: icepack_init_parameters @@ -2286,7 +2278,6 @@ icepack_step_therm2 dt , & ! time step Tf , & ! freezing temperature (C) sss , & ! sea surface salinity (ppt) - rside , & ! fraction of ice that melts laterally frzmlt ! freezing/melting potential (W/m^2) integer (kind=int_kind), dimension (:), intent(in) :: & @@ -2301,13 +2292,13 @@ icepack_step_therm2 nt_strata ! indices of underlying tracer layers real (kind=dbl_kind), dimension(:), intent(in) :: & + rsiden , & ! fraction of ice that melts laterally salinz , & ! initial salinity profile ocean_bio ! ocean concentration of biological tracer real (kind=dbl_kind), intent(inout) :: & aice , & ! sea ice concentration aice0 , & ! concentration of open water - fside , & ! lateral heat flux (W/m^2) frain , & ! rainfall rate (kg/m^2 s) fpond , & ! fresh water flux to ponds (kg/m^2/s) fresh , & ! fresh water flux to ocean (kg/m^2/s) @@ -2370,10 +2361,6 @@ icepack_step_therm2 d_afsd_latm, & ! lateral melt d_afsd_weld ! welding - real (kind=dbl_kind), dimension (:), intent(in), optional :: & - floe_rad_c, & ! fsd size bin centre in m (radius) - floe_binwidth ! fsd size bin width in m (radius) - icepack_therm_shared.F90 @@ -2563,8 +2550,8 @@ icepack_step_therm1 strocnxT , strocnyT , & fbot , & Tbot , Tsnice , & - frzmlt , rside , & - fside , wlat , & + frzmlt , rsiden , & + wlat , & fsnow , frain , & fpond , fsloss , & fsurf , fsurfn , & @@ -2611,7 +2598,7 @@ icepack_step_therm1 lmask_n , lmask_s , & mlt_onset , frz_onset , & yday , prescribed_ice, & - zlvs) + zlvs , afsdn) real (kind=dbl_kind), intent(in) :: & dt , & ! time step @@ -2687,8 +2674,6 @@ icepack_step_therm1 strocnyT , & ! ice-ocean stress, y-direction 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) @@ -2701,7 +2686,7 @@ icepack_step_therm1 frz_onset ! day of year that freezing begins (congel or frazil) real (kind=dbl_kind), intent(out), optional :: & - wlat ! lateral melt rate (m/s) + wlat ! lateral melt rate (m/s) real (kind=dbl_kind), intent(inout), optional :: & fswthru_vdr , & ! vis dir shortwave penetrating to ocean (W/m^2) @@ -2735,6 +2720,9 @@ icepack_step_therm1 H2_18O_ocn , & ! ocean concentration of H2_18O (kg/kg) zlvs ! atm level height for scalars (if different than zlvl) (m) + real (kind=dbl_kind), dimension(:,:), intent(in), optional :: & + afsdn ! afsd tracer + real (kind=dbl_kind), dimension(:), intent(inout) :: & aicen_init , & ! fractional area of ice vicen_init , & ! volume per unit area of ice (m) @@ -2750,6 +2738,7 @@ icepack_step_therm1 ipnd , & ! melt pond refrozen lid thickness (m) iage , & ! volume-weighted ice age FY , & ! area-weighted first-year ice area + rsiden , & ! fraction of ice that melts laterally fsurfn , & ! net flux to top surface, excluding fcondtop fcondtopn , & ! downward cond flux at top surface (W m-2) fcondbotn , & ! downward cond flux at bottom surface (W m-2) @@ -3380,7 +3369,6 @@ icepack_step_wavefracture subroutine icepack_step_wavefracture(wave_spec_type, & dt, nfreq, & aice, vice, aicen, & - floe_rad_l, floe_rad_c, & wave_spectrum, wavefreq, dwavefreq, & trcrn, d_afsd_wave) @@ -3399,10 +3387,6 @@ icepack_step_wavefracture 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)