diff --git a/columnphysics/icepack_fsd.F90 b/columnphysics/icepack_fsd.F90 index 28f12d92..2e5175af 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 02a9ab4b..0fcc9673 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 6dda193a..bf96096e 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 1beb7a23..34c4e448 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 c45d2802..ba1b1de6 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 b55f129f..dbe1ab56 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 39c09081..321a8491 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 b3f360f7..98e390e2 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 ddcd8b2f..f342d74f 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 f99d47e4..13f6f2f8 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 ae4d64d9..4ec03fdd 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 ae640d1b..b41f6055 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)