!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! ! ! THREE-LAYER VERTICAL THERMODYNAMICS ! ! ! ! Reference: M. Winton, 2000: "A reformulated three-layer sea ice model", ! ! Journal of Atmospheric and Oceanic Technology, 17, 525-531. ! ! ! ! ! ! -> +---------+ <- ts - diagnostic surface temperature ( <= 0C ) ! ! / | | ! ! hs | snow | <- 0-heat capacity snow layer ! ! \ | | ! ! => +---------+ ! ! / | | ! ! / | | <- t1 - upper 1/2 ice temperature; this layer has ! ! / | | a variable (T/S dependent) heat capacity ! ! hi |...ice...| ! ! \ | | ! ! \ | | <- t2 - lower 1/2 ice temp. (fixed heat capacity) ! ! \ | | ! ! -> +---------+ <- base of ice fixed at seawater freezing temp. ! ! ! ! Mike Winton (Michael.Winton)! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! TK mod: SLAB_ICE treatment modified to follow supersource ! (after Bryan 1969). The conductive heat ! flux from ice to atmosphere is computed based on ! an effective ice thickness which ensures a minimum ! thickness of 1.7cm for the calculation. module ice_thm_mod #if 0 use constants_mod, only : LI => hlf ! latent heat of fusion - 334e3 J/(kg-ice) #endif implicit none integer, parameter, private :: rk = kind(1.d0) #if 1 real(rk) :: LI=334.e3 #endif private public :: DS, DI, DW, MU_TS, TFI, ice_optics, ice3lay_temp, ice3lay_resize, & thm_pack, thm_unpack, ice_thm_param, e_to_melt ! ! properties of ice, snow, and seawater (NCAR CSM values) ! real(rk) :: KS = 0.31 ! conductivity of snow - 0.31 W/(mK) real(rk), parameter :: DS = 330.0 ! density of snow - 330 kg/(m^3) real(rk), parameter :: KI = 2.03 ! conductivity of ice - 2.03 W/(mK) real(rk), parameter :: DI = 905.0 ! density of ice - 905 kg/(m^3) real(rk), parameter :: CI = 21e2 ! heat cap. of fresh ice - 2100 J/(kg K) real(rk), parameter :: SI = 1.0 ! salinity of sea ice real(rk), parameter :: MU_TS = 0.054 ! relates freezing temp. to salinity real(rk), parameter :: TFI = -MU_TS*SI ! sea ice freezing temp. = -mu*salinity real(rk), parameter :: CW = 4.2e3 ! heat capacity of seawater real(rk), parameter :: DW = 1030.0 ! density of water for waterline - kg/(m^3) ! albedos are from CSIM4 assumming 0.53 visible and 0.47 near-ir insolation real(rk) :: ALB_SNO = 0.85 ! albedo of snow (not melting) real(rk) :: ALB_ICE = 0.5826 ! albedo of ice (not melting) real(rk) :: PEN_ICE = 0.3 ! ice surface penetrating solar fraction real(rk) :: OPT_DEP_ICE = 0.67 ! ice optical depth (m) real(rk) :: T_RANGE_MELT = 1.0 ! melt albedos scaled in below melting T real(rk) :: H_LO_LIM = 0.0 ! hi/hs lower limit for temp. calc. logical :: SLAB_ICE = .false. ! should we do old style GFDL slab ice? ! ! slab ice specific parameters ! real(rk), parameter :: CRIT_THICKNESS = 1.00 real(rk), parameter :: T_RANGE = 10.0 real(rk), parameter :: MIN_ICE_ALB = 0.55 ! coupled model uses 0.55 real(rk), parameter :: MAX_ICE_ALB = 0.80 real(rk), parameter :: ALB_OCEAN = 0.10 ! logical :: CM2_BUGS = .false. ! keep cm2 bugs for reproducibility contains !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! ice_thm_param - set ice thermodynamic parameters ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine ice_thm_param(alb_sno_in, alb_ice_in, pen_ice_in, opt_dep_ice_in, & slab_ice_in, t_range_melt_in, cm2_bugs_in, ks_in, & h_lo_lim_in ) real(rk), intent(in) :: alb_sno_in, alb_ice_in, pen_ice_in real(rk), intent(in) :: opt_dep_ice_in, t_range_melt_in logical, intent(in) :: slab_ice_in logical, intent(in) :: cm2_bugs_in real(rk), intent(in) :: ks_in real(rk), intent(in) :: h_lo_lim_in ALB_SNO = alb_sno_in ALB_ICE = alb_ice_in PEN_ICE = pen_ice_in OPT_DEP_ICE = opt_dep_ice_in SLAB_ICE = slab_ice_in T_RANGE_MELT = t_range_melt_in CM2_BUGS = cm2_bugs_in KS = ks_in H_LO_LIM = h_lo_lim_in end subroutine ice_thm_param !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! ice_optics - set albedo, penetrating solar, and ice/snow transmissivity ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine ice_optics(alb, pen, trn, hs, hi, ts, tfw) real(rk), intent( out) :: alb ! ice surface albedo (0-1) real(rk), intent( out) :: pen ! fraction of down solar penetrating the ice real(rk), intent( out) :: trn ! ratio of down solar at bottom to top of ice real(rk), intent(in ) :: hs ! snow thickness (m-snow) real(rk), intent(in ) :: hi ! ice thickness (m-ice) real(rk), intent(in ) :: ts ! surface temperature real(rk), intent(in ) :: tfw ! seawater freezing temperature real(rk) :: as, ai, cs real(rk) :: thick_ice_alb, tcrit, fh if (SLAB_ICE) then tcrit = tfw - T_RANGE if (ts <= tcrit) then thick_ice_alb = MAX_ICE_ALB else if (ts >= tfw) then thick_ice_alb = MIN_ICE_ALB else thick_ice_alb = MAX_ICE_ALB + (MIN_ICE_ALB-MAX_ICE_ALB)*(ts-tcrit)/T_RANGE endif if (hi >= crit_thickness) then alb = THICK_ICE_ALB else alb = ALB_OCEAN + (thick_ice_alb-ALB_OCEAN)*sqrt(hi/CRIT_THICKNESS) endif pen = 0.0 trn = 0.0 !! check for ice albdeos out of range (0 to 1) ! if (alb.lt.0.0 .or.alb.gt.1.0) then ! print *,'ice_optics: albedo out of range, alb_in=',alb_in, 'alb=',alb ! print *,'ts=',ts, 'tfw=',tfw, 'tcrit=',tcrit ! print *,'hi=',hi, 'thick_ice_alb=',thick_ice_alb ! print *,'ALB_OCEAN=',ALB_OCEAN ! print *,'MIN_ICE_ALB=',MIN_ICE_ALB, 'MAX_ICE_ALB=',MAX_ICE_ALB ! print *,'T_RANGE,=',T_RANGE, 'CRIT_THICKNESS=',CRIT_THICKNESS ! stop ! end if return endif !! 2007/04/11 Fix for thin ice negative ice albedos from Mike Winton !! Move ai calculation after if test as = ALB_SNO; ai = ALB_ICE cs = hs/(hs+0.02) ! thin snow partially covers ice fh = min(atan(5.0*hi)/atan(5.0*0.5),1.0) ! use this form from CSIM4 to ! reduce albedo for thin ice if (CM2_BUGS) then ai = fh*ai+(1-fh)*0.06 ! reduce albedo for thin ice if (ts+T_RANGE_MELT > TFI) then ! reduce albedo for melting as in ! CSIM4 assuming 0.53/0.47 vis/ir as = as-0.1235*min((ts+T_RANGE_MELT-TFI)/T_RANGE_MELT,1.0) ai = ai-0.075 *min((ts+T_RANGE_MELT-TFI)/T_RANGE_MELT,1.0) endif else if (ts+T_RANGE_MELT > TFI) then ! reduce albedo for melting as in ! CSIM4 assuming 0.53/0.47 vis/ir as = as-0.1235*min((ts+T_RANGE_MELT-TFI)/T_RANGE_MELT,1.0) ai = ai-0.075 *min((ts+T_RANGE_MELT-TFI)/T_RANGE_MELT,1.0) endif ai = fh*ai+(1-fh)*0.06 ! reduce albedo for thin ice end if alb = cs*as+(1-cs)*ai pen = (1-cs)*PEN_ICE trn = exp(-hi/OPT_DEP_ICE); !! check for ice albdeos out of range (0 to 1) ! if (alb.lt.0.0 .or. alb.gt.1.0) then ! print *,'ice_optics: albedo out of range, alb=',alb ! print *,'cs=',cs, 'as=',as, 'ai=',ai ! print *,'ts=',ts, 'fh=',fh, 'hs=',hs, 'hi=',hi, 'tfw=',tfw ! print *,'ALB_SNO=',ALB_SNO, 'ALB_ICE=',ALB_ICE, 'T_RANGE_MELT,=',T_RANGE_MELT, 'TFI=',TFI ! stop ! end if return end subroutine ice_optics !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! ice3lay_temp - ice & snow temp. change [Winton (2000) section 2.a] ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine ice3lay_temp(hs, hi, t1, t2, ts, A, B, I, tfw, fb, dt, tmelt, bmelt) implicit none real(rk), intent(in ) :: hs ! snow thickness (m) real(rk), intent(in ) :: hi ! ice thickness (m) real(rk), intent( out) :: ts ! surface temperature (deg-C) real(rk), intent(inout) :: t1 ! upper ice temperature (deg-C) real(rk), intent(inout) :: t2 ! lower ice temperature (deg-C) real(rk), intent(in ) :: A ! net surface heat flux (+ up) at ts=0 (W/m^2) real(rk), intent(in ) :: B ! d(sfc heat flux)/d(ts) [W/(m^2 deg-C)] real(rk), intent(in ) :: I ! solar absorbed by upper ice (W/m^2) real(rk), intent(in ) :: tfw ! seawater freezing temperature (deg-C) real(rk), intent(in ) :: fb ! heat flux from ocean to ice bottom (W/m^2) real(rk), intent(in ) :: dt ! timestep (sec) real(rk), intent(inout) :: tmelt ! accumulated top melting energy (J/m^2) real(rk), intent(inout) :: bmelt ! accumulated bottom melting energy (J/m^2) ! ! variables for temperature calculation [see Winton (1999) section II.A.] ! note: here equations are multiplied by hi to improve thin ice accuracy ! real(rk) :: tsf, k12, hi2, a10, b10, a1, b1, c1 ! TK Mods: real(rk) :: hi_effective, hie real(rk) :: KI_over_eps = 1.7065e-2 ! 5/2.93 from Bryan (1969); ! Value used in SS tsc.F (1.7065 cm) ! converted to meters... if (SLAB_ICE) then hi_effective = hi + KI_over_eps ! TK added ts = (KI*tfw-A*hi_effective)/(KI+B*hi_effective) ! TK mod if (ts > 0.0) then ! surface melting conditions ts = 0.0 if (hi>0.0) tmelt = tmelt + (KI*tfw/hi_effective-A)*dt ! TK mod endif if (hi>0.0) then bmelt = bmelt + (fb-KI*(tfw-ts)/hi_effective)*dt ! TK mod else bmelt = bmelt + (fb-A-B*tfw)*dt endif return endif if (hs > 0.0) then TSF = 0.0 else TSF = TFI endif hie = max(hi, H_LO_LIM); ! prevent thin ice inaccuracy (mw) ! ! Compute upper ice and surface temperatures ! K12 = 4*KI*KS/(KS+4*KI*hs/hie) hi2 = hie*hie A10 = DI*hi2*CI/(2*dt) + 2*KI*(4*dt*2*KI+DI*hi2*CI)/(6*dt*2*KI+DI*hi2*CI) B10 = -DI*hi2*(CI*t1+LI*TFI/t1)/(2*dt) - I*hie & -2*KI*(4*dt*2*KI*tfw+DI*hi2*CI*t2)/(6*dt*2*KI+DI*hi2*CI) A1 = A10+K12*B*hie/(K12+B*hie) B1 = B10+A*K12*hie/(K12+B*hie) C1 = DI*hi2*LI*TFI/(2*dt) t1 = -(sqrt(B1*B1-4*A1*C1)+B1)/(2*A1) ts = (K12*t1-A*hie)/(K12+B*hie) ! if (ts > tsf) then ! slightly different equation for melting conditions if (ts > -0.054) then ! slightly different equation for melting conditions A1 = A10+K12 B1 = B10-K12*tsf t1 = -(sqrt(B1*B1-4*A1*C1)+B1)/(2*A1) ts = tsf tmelt = tmelt + (K12*(t1-ts)/hie-(A+B*ts))*dt endif ! ! set lower ice temp. -- use tfw as reference for thin ice precision ! t1 = t1-tfw; t2 = t2-tfw; t2 = (2*dt*2*KI*t1+DI*hi2*CI*t2)/(6*dt*2*KI+DI*hi2*CI) t1 = t1+tfw; t2 = t2+tfw; bmelt = bmelt + (fb+4*KI*(t2-tfw)/hie)*dt if (t2 > TFI) then ! put excess lower ice energy into bmelt bmelt = bmelt + e_to_melt(h2=hie/2,t2=TFI) - e_to_melt(h2=hie/2,t2=t2) t2 = TFI endif if (t1 > TFI) then ! put excess upper ice energy into tmelt tmelt = tmelt + e_to_melt(h1=hie/2,t1=TFI) - e_to_melt(h1=hie/2,t1=t1) t1 = TFI endif ! if (tmelt<0) then ! print *,'neg. tmelt=',tmelt,ts,t1,hs,hi,K12*(t1-ts)/hi,-(A+B*ts) ! print *,'K12=',K12 ! print *,'A/B/I=',A,B,I ! print *,'A/B/C=',A1,B1,C1,A1*t1*t1+B1*t1+C1,A1*t1*t1,B1*t1,C1 ! stop ! end if !KB call thm_checkout(ts, hs, hi, t1, t2, bmelt, tmelt) return end subroutine ice3lay_temp subroutine thm_checkout(ts, hs, hi, t1, t2, bmelt, tmelt) real(rk), intent(in) :: ts, hs, hi, t1, t2, bmelt, tmelt if (ts>0 .or. t1>TFI .or. t2>0 .or. hs<0 .or. hs>100 .or. hi<0 .or. hi>100 & .or. abs(bmelt)>100*DI*LI .or. tmelt<0 .or. tmelt>100*DI*LI) then print *,'UNREASONABLE ICE: hs=',hs,'hi=',hi,'t1=',t1,'t2=',t2,'ts=', & ts,'tmelt=',tmelt,'bmelt=',bmelt end if end subroutine thm_checkout !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! e_to_melt - energy needed to melt a given snow/ice configuration ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! function e_to_melt(hs, h1, t1, h2, t2) real(rk), intent(in), optional :: hs, h1, t1, h2, t2 real(rk) :: e_to_melt e_to_melt = 0.0 if (present(hs)) e_to_melt = e_to_melt+DS*LI*hs if (present(h1).and.present(t1)) then e_to_melt = e_to_melt+DI*h1*(CI-LI/t1)*(TFI-t1) endif if (present(h2).and.present(t2)) then e_to_melt = e_to_melt+DI*h2*(LI+CI*(TFI-t2)) endif end function e_to_melt !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! add_to_top - add some ice to the top ice layer ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine add_to_top(h, t, h1, t1) real(rk), intent(in) :: h, t real(rk), intent(inout) :: h1, t1 real(rk) :: f1 f1 = h1/(h1+h) t1 = f1*(t1+LI*TFI/(CI*t1))+(1-f1)*t t1 = (t1-sqrt(t1*t1-4*TFI*LI/CI))/2 h1 = h1+h return end subroutine add_to_top !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! add_to_bot - add some ice to the bottom ice layer ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine add_to_bot(h, t, h2, t2) real(rk), intent(in) :: h, t real(rk), intent(inout) :: h2, t2 t2 = (h2*t2+h*t)/(h2+h) h2 = h2+h return end subroutine add_to_bot !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! even_up - transfer mass/energy between ice layers to maintain equal thickness! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine even_up(h1, t1, h2, t2) real(rk),intent(inout) :: h1, t1, h2, t2 real(rk) :: dh if (h1 > (h1+h2)/2) then call add_to_bot(h1-(h1+h2)/2, t1+LI*TFI/(CI*t1), h2, t2) h1 = h2 else if (h2 > (h1+h2)/2) then call add_to_top(h2-(h1+h2)/2, t2, h1, t1) h2 = h1 endif if (t2>TFI) then ! use extra energy to melt both layers evenly dh = h2*CI*(t2-TFI)*t1/(LI*t1+(CI*t1-LI)*(TFI-t1)) t2 = TFI h1 = h1-dh h2 = h2-dh endif return end subroutine even_up !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! ice3lay_resize - ice & snow thickness change [Winton (1998) section II.B.] ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine ice3lay_resize(hs, hi, t1, t2, snow, frazil, evap, tmelt, bmelt, & tfw, heat_to_ocn, h2o_to_ocn, h2o_from_ocn, & snow_to_ice, bablt ) real(rk), intent(inout) :: hs ! snow thickness (m-snow) real(rk), intent(inout) :: hi ! ice thickness (m-ice) real(rk), intent(inout) :: t1 ! temperature of upper ice (deg-C) real(rk), intent(inout) :: t2 ! temperature of lower ice (deg-C) real(rk), intent(in ) :: snow ! new snow (kg/m^2-snow) real(rk), intent(in ) :: frazil ! frazil in energy units real(rk), intent(in ) :: evap ! ice evaporation (kg/m^2) real(rk), intent(in ) :: tmelt ! top melting energy (J/m^2) real(rk), intent(in ) :: bmelt ! bottom melting energy (J/m^2) real(rk), intent(in ) :: tfw ! seawater freezing temperature (deg-C) real(rk), intent( out) :: heat_to_ocn ! energy left after ice all melted (J/m^2) real(rk), intent( out) :: h2o_to_ocn ! liquid water flux to ocean (kg/m^2) real(rk), intent( out) :: h2o_from_ocn! evaporation flux from ocean (kg/m^2) real(rk), intent( out) :: snow_to_ice ! snow below waterline becomes ice real(rk), intent( out), optional :: bablt ! bottom ablation (kg/m^2) real(rk) h1, h2, dh, hw real(rk) h1o, h2o heat_to_ocn = 0.0 h2o_to_ocn = DS*hs+DI*hi+snow-evap ! - from ice at end gives ocean h2o flux h2o_from_ocn = 0.0 snow_to_ice = 0.0 if (SLAB_ICE) then ! ! add snow and frazil ! hi = hi + snow/DI + frazil/(DI*LI) t1 = tfw; t2 = tfw; ! ! atmospheric evaporation ! if (evap <= hi*DI) then hi = hi - evap/DI else h2o_from_ocn = evap-hi*DI hi = 0.0 end if ! ! ... melting ! hi = hi - (tmelt+bmelt)/(DI*LI) if (hi<0.0) then heat_to_ocn = -hi*DI*LI hi = 0.0 else heat_to_ocn = 0.0 endif h2o_to_ocn = h2o_to_ocn+h2o_from_ocn ! reset mark for leftover evap thru ice h2o_to_ocn = h2o_to_ocn-DI*hi-DS*hs ! hs should be zero if (present(bablt)) bablt = bmelt/LI return endif h1 = hi/2 h2 = hi/2 ! ! add snow ... ! hs = hs + snow/DS ! ! ... and frazil ! call add_to_bot(frazil/e_to_melt(h2=1._rk,t2=tfw), tfw, h2, t2) ! ! atmospheric evaporation ! #if 0 if (evap <= hs*DS) then hs = hs - evap/DS else if (evap-hs*DS<=h1*DI) then hs = 0.0 h1 = h1-(evap-DS*hs)/DI else if (evap-hs*DS-h1*DI<=h2*DI) then hs = 0.0 h1 = 0.0 h2 = h2 - (evap-hs*DS-h1*DI)/DI else h2o_from_ocn = evap-hs*DS-(h1+h2)*DI hs = 0.0 h1 = 0.0 h2 = 0.0 end if #endif if (bmelt < 0.0) then h2o = h2 call add_to_bot(-bmelt/e_to_melt(h2=1._rk,t2=tfw), tfw, h2, t2) end if if (h1 == 0.0) t1 = tfw ! need this, below we divide by t1 even when h1 == 0 ! ! apply energy fluxes ... top ... ! h1o = h1 if (tmelt <= e_to_melt(hs)) then !KB hs = hs - tmelt/e_to_melt(hs=1._rk) else if (tmelt <= e_to_melt(hs,h1,t1)) then h1 = h1 - (tmelt-e_to_melt(hs))/e_to_melt(h1=1._rk,t1=t1) hs = 0.0 else if (tmelt <= e_to_melt(hs,h1,t1,h2,t2)) then h2 = h2 - (tmelt-e_to_melt(hs,h1,t1))/e_to_melt(h2=1._rk,t2=t2) hs = 0.0 h1 = 0.0 else heat_to_ocn = heat_to_ocn+tmelt-e_to_melt(hs,h1,t1,h2,t2) hs = 0.0 h1 = 0.0 h2 = 0.0 endif ! ! ... and bottom ! if (present(bablt)) bablt = DS*hs+DI*(h1+h2) if (bmelt > 0.0) then h2o = h2 if (bmelt < e_to_melt(h2=h2,t2=t2)) then h2 = h2 - bmelt/e_to_melt(h2=1._rk,t2=t2) else if (bmelt < e_to_melt(h1=h1,t1=t1,h2=h2,t2=t2)) then h1 = h1-(bmelt-e_to_melt(h2=h2,t2=t2))/e_to_melt(h1=1._rk,t1=t1) h2 = 0.0 else if (bmelt < e_to_melt(hs,h1,t1,h2,t2)) then !KB hs = hs - (bmelt-e_to_melt(h1=h1,t1=t1,h2=h2,t2=t2))/e_to_melt(hs=1._rk) h1 = 0.0 h2 = 0.0 else heat_to_ocn = heat_to_ocn+bmelt-e_to_melt(hs,h1,t1,h2,t2) hs = 0.0 h1 = 0.0 h2 = 0.0 endif endif if (present(bablt)) bablt = bablt-DS*hs-DI*(h1+h2) hi = h1 + h2 hw = (DI*hi+DS*hs)/DW if (hw>hi) then ! convert snow to ice to maintain ice at waterline snow_to_ice = (hw-hi)*DI !KB hs = hs - snow_to_ice/DS call add_to_top(hw-hi, TFI, h1, t1) endif call even_up(h1, t1, h2, t2) hi = h1+h2 if (hi==0.0) then t1 = 0.0 t2 = 0.0 endif h2o_to_ocn = h2o_to_ocn+h2o_from_ocn ! correct mark for leftover evap thru ice h2o_to_ocn = h2o_to_ocn-DS*hs-DI*hi return end subroutine ice3lay_resize !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! thm_pack - form conserved quantities for two-dimensional advection ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine thm_pack(cn, hs, hi, t1, t2) real(rk), dimension(:,:,:), intent(inout) :: cn, hs, hi, t1, t2 integer :: i, j, k real(rk) :: tmp tmp = MU_TS*SI*LI/CI do k = 1, size(hi,3) do j = 1, size(hi,2) do i = 1, size(hi,1) if (hi(i,j,k)>0.0) then hi(i,j,k) = cn(i,j,k)*hi(i,j,k) hs(i,j,k) = cn(i,j,k)*hs(i,j,k) t1(i,j,k) = (t1(i,j,k)-tmp/t1(i,j,k))*hi(i,j,k) t2(i,j,k) = t2(i,j,k)*hi(i,j,k) else cn(i,j,k) = 0.0 hi(i,j,k) = 0.0 hs(i,j,k) = 0.0 t1(i,j,k) = 0.0 t2(i,j,k) = 0.0 endif enddo enddo enddo return end subroutine thm_pack !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! thm_unpack - reform ice properties from advectively conserved quantities ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! subroutine thm_unpack(cn, hs, hi, t1, t2) real(rk), dimension(:,:,:), intent(inout) :: cn, hs, hi, t1, t2 integer :: i, j, k real(rk) :: tmp tmp = 4*MU_TS*SI*LI/CI do k = 1, size(hi,3) do j = 1, size(hi,2) do i = 1, size(hi,1) if (hi(i,j,k)>0.0) then t1(i,j,k) = t1(i,j,k)/hi(i,j,k) t1(i,j,k) = 0.5*(t1(i,j,k)-sqrt(t1(i,j,k)*t1(i,j,k)+tmp)) t2(i,j,k) = t2(i,j,k)/hi(i,j,k) hi(i,j,k) = hi(i,j,k)/cn(i,j,k) hs(i,j,k) = hs(i,j,k)/cn(i,j,k) else cn(i,j,k) = 0.0 hi(i,j,k) = 0.0 hs(i,j,k) = 0.0 t1(i,j,k) = 0.0 t2(i,j,k) = 0.0 endif enddo enddo enddo return end subroutine thm_unpack end module ice_thm_mod