!========================================================================= ! ! Update ice and snow internal temperatures ! using Bitz and Lipscomb 1999 thermodynamics ! ! authors: William H. Lipscomb, LANL ! C. M. Bitz, UW ! Elizabeth C. Hunke, LANL ! ! 2012: Split from icepack_therm_vertical.F90 module icepack_therm_bl99 use icepack_kinds use icepack_parameters, only: c0, c1, c2, p1, p5, puny #ifdef CESMCOUPLED use icepack_parameters, only: p01 #endif use icepack_parameters, only: rhoi, rhos, hs_min, cp_ice, cp_ocn, depressT, Lfresh, ksno, kice use icepack_parameters, only: conduct, calc_Tsfc, solve_zsal use icepack_parameters, only: sw_redist, sw_frac, sw_dtemp use icepack_warnings, only: warnstr, icepack_warnings_add use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted use icepack_therm_shared, only: ferrmax, l_brine use icepack_therm_shared, only: surface_heat_flux, dsurface_heat_flux_dTsf implicit none private public :: surface_fluxes, temperature_changes real (kind=dbl_kind), parameter :: & betak = 0.13_dbl_kind, & ! constant in formula for k (W m-1 ppt-1) kimin = 0.10_dbl_kind ! min conductivity of saline ice (W m-1 deg-1) !======================================================================= contains !======================================================================= ! ! Compute new surface temperature and internal ice and snow ! temperatures. Include effects of salinity on sea ice heat ! capacity in a way that conserves energy (Bitz and Lipscomb, 1999). ! ! New temperatures are computed iteratively by solving a tridiagonal ! system of equations; heat capacity is updated with each iteration. ! Finite differencing is backward implicit. ! ! See Bitz, C.M., and W.H. Lipscomb, 1999: ! An energy-conserving thermodynamic model of sea ice, ! J. Geophys. Res., 104, 15,669-15,677. ! ! authors William H. Lipscomb, LANL ! C. M. Bitz, UW subroutine temperature_changes (dt, & nilyr, nslyr, & rhoa, flw, & potT, Qa, & shcoef, lhcoef, & fswsfc, fswint, & Sswabs, Iswabs, & hilyr, hslyr, & zqin, zTin, & zqsn, zTsn, & zSin, & Tsf, Tbot, & fsensn, flatn, & flwoutn, fsurfn, & fcondtopn,fcondbot, & einit ) integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers nslyr ! number of snow layers real (kind=dbl_kind), intent(in) :: & dt ! time step real (kind=dbl_kind), & intent(in) :: & rhoa , & ! air density (kg/m^3) flw , & ! incoming longwave radiation (W/m^2) potT , & ! air potential temperature (K) Qa , & ! specific humidity (kg/kg) shcoef , & ! transfer coefficient for sensible heat lhcoef , & ! transfer coefficient for latent heat Tbot ! ice bottom surface temperature (deg C) real (kind=dbl_kind), & intent(inout) :: & fswsfc , & ! SW absorbed at ice/snow surface (W m-2) fswint ! SW absorbed in ice interior below surface (W m-2) real (kind=dbl_kind), intent(in) :: & hilyr , & ! ice layer thickness (m) hslyr , & ! snow layer thickness (m) einit ! initial energy of melting (J m-2) real (kind=dbl_kind), dimension (nslyr), & intent(inout) :: & Sswabs ! SW radiation absorbed in snow layers (W m-2) real (kind=dbl_kind), dimension (nilyr), & intent(inout) :: & Iswabs ! SW radiation absorbed in ice layers (W m-2) real (kind=dbl_kind), intent(inout):: & fsurfn , & ! net flux to top surface, excluding fcondtopn fcondtopn , & ! downward cond flux at top surface (W m-2) fsensn , & ! surface downward sensible heat (W m-2) flatn , & ! surface downward latent heat (W m-2) flwoutn ! upward LW at surface (W m-2) real (kind=dbl_kind), intent(out):: & fcondbot ! downward cond flux at bottom surface (W m-2) real (kind=dbl_kind), & intent(inout):: & Tsf ! ice/snow surface temperature, Tsfcn real (kind=dbl_kind), dimension (nilyr), & intent(inout) :: & zqin , & ! ice layer enthalpy (J m-3) zTin ! internal ice layer temperatures real (kind=dbl_kind), dimension (nilyr), & intent(in) :: & zSin ! internal ice layer salinities real (kind=dbl_kind), dimension (nslyr), & intent(inout) :: & zqsn , & ! snow layer enthalpy (J m-3) zTsn ! internal snow layer temperatures ! local variables integer (kind=int_kind), parameter :: & nitermax = 200 ! max number of iterations in temperature solver real (kind=dbl_kind), parameter :: & Tsf_errmax = 5.e-4_dbl_kind ! max allowed error in Tsf ! recommend Tsf_errmax < 0.01 K integer (kind=int_kind) :: & k , & ! ice layer index niter , & ! iteration counter in temperature solver nmat ! matrix dimension logical (kind=log_kind) :: & l_snow , & ! true if snow temperatures are computed l_cold ! true if surface temperature is computed real (kind=dbl_kind) :: & Tsf_start , & ! Tsf at start of iteration dTsf , & ! Tsf - Tsf_start dTi1 , & ! Ti1(1) - Tin_start(1) dfsurf_dT , & ! derivative of fsurf wrt Tsf avg_Tsi , & ! = 1. if new snow/ice temps averaged w/starting temps enew ! new energy of melting after temp change (J m-2) real (kind=dbl_kind) :: & dTsf_prev , & ! dTsf from previous iteration dTi1_prev , & ! dTi1 from previous iteration dfsens_dT , & ! deriv of fsens wrt Tsf (W m-2 deg-1) dflat_dT , & ! deriv of flat wrt Tsf (W m-2 deg-1) dflwout_dT , & ! deriv of flwout wrt Tsf (W m-2 deg-1) dt_rhoi_hlyr, & ! dt/(rhoi*hilyr) einex , & ! excess energy from dqmat to ocean ferr ! energy conservation error (W m-2) real (kind=dbl_kind), dimension (nilyr) :: & Tin_init , & ! zTin at beginning of time step Tin_start , & ! zTin at start of iteration dTmat , & ! zTin - matrix solution before limiting dqmat , & ! associated enthalpy difference Tmlts ! melting temp, -depressT * salinity real (kind=dbl_kind), dimension (nslyr) :: & Tsn_init , & ! zTsn at beginning of time step Tsn_start , & ! zTsn at start of iteration etas ! dt / (rho * cp * h) for snow layers real (kind=dbl_kind), dimension (nilyr+nslyr+1) :: & sbdiag , & ! sub-diagonal matrix elements diag , & ! diagonal matrix elements spdiag , & ! super-diagonal matrix elements rhs , & ! rhs of tri-diagonal matrix equation Tmat ! matrix output temperatures real (kind=dbl_kind), dimension (nilyr) :: & etai ! dt / (rho * cp * h) for ice layers real (kind=dbl_kind), dimension(nilyr+nslyr+1):: & kh ! effective conductivity at interfaces (W m-2 deg-1) real (kind=dbl_kind) :: & ci , & ! specific heat of sea ice (J kg-1 deg-1) avg_Tsf , & ! = 1. if Tsf averaged w/Tsf_start, else = 0. Iswabs_tmp , & ! energy to melt through fraction frac of layer Sswabs_tmp , & ! same for snow dswabs , & ! difference in swabs and swabs_tmp frac logical (kind=log_kind) :: & converged ! = true when local solution has converged logical (kind=log_kind) , dimension (nilyr) :: & reduce_kh ! reduce conductivity when T exceeds Tmlt character(len=*),parameter :: subname='(temperature_changes)' !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- converged = .false. l_snow = .false. l_cold = .true. fcondbot = c0 dTsf_prev = c0 dTi1_prev = c0 dfsens_dT = c0 dflat_dT = c0 dflwout_dT = c0 einex = c0 dt_rhoi_hlyr = dt / (rhoi*hilyr) ! hilyr > 0 if (hslyr > hs_min/real(nslyr,kind=dbl_kind)) & l_snow = .true. do k = 1, nslyr Tsn_init (k) = zTsn(k) ! beginning of time step Tsn_start(k) = zTsn(k) ! beginning of iteration if (l_snow) then etas(k) = dt/(rhos*cp_ice*hslyr) else etas(k) = c0 endif enddo ! k do k = 1, nilyr Tin_init (k) = zTin(k) ! beginning of time step Tin_start(k) = zTin(k) ! beginning of iteration Tmlts (k) = -zSin(k) * depressT enddo !----------------------------------------------------------------- ! Compute thermal conductivity at interfaces (held fixed during ! subsequent iterations). ! Ice and snow interfaces are combined into one array (kh) to ! simplify the logic. !----------------------------------------------------------------- call conductivity (l_snow, & nilyr, nslyr, & hilyr, hslyr, & zTin, kh, zSin) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! Check for excessive absorbed solar radiation that may result in ! temperature overshoots. Convergence is particularly difficult ! if the starting temperature is already very close to the melting ! temperature and extra energy is added. In that case, or if the ! amount of energy absorbed is greater than the amount needed to ! melt through a given fraction of a layer, we put the extra ! energy into the surface. ! NOTE: This option is not available if the atmosphere model ! has already computed fsurf. (Unless we adjust fsurf here) !----------------------------------------------------------------- !mclaren: Should there be an if calc_Tsfc statement here then?? if (sw_redist) then if (solve_zsal) sw_dtemp = p1 ! lower tolerance with dynamic salinity do k = 1, nilyr Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc if (Tin_init(k) <= Tmlts(k) - sw_dtemp) then if (l_brine) then ci = cp_ice - Lfresh * Tmlts(k) / (Tin_init(k)**2) Iswabs_tmp = min(Iswabs(k), & sw_frac*(Tmlts(k)-Tin_init(k))*ci/dt_rhoi_hlyr) else ci = cp_ice Iswabs_tmp = min(Iswabs(k), & sw_frac*(-Tin_init(k))*ci/dt_rhoi_hlyr) endif endif if (Iswabs_tmp < puny) Iswabs_tmp = c0 dswabs = min(Iswabs(k) - Iswabs_tmp, fswint) fswsfc = fswsfc + dswabs fswint = fswint - dswabs Iswabs(k) = Iswabs_tmp enddo do k = 1, nslyr if (l_snow) then Sswabs_tmp = c0 if (Tsn_init(k) <= -sw_dtemp) then Sswabs_tmp = min(Sswabs(k), & -sw_frac*Tsn_init(k)/etas(k)) endif if (Sswabs_tmp < puny) Sswabs_tmp = c0 dswabs = min(Sswabs(k) - Sswabs_tmp, fswint) fswsfc = fswsfc + dswabs fswint = fswint - dswabs Sswabs(k) = Sswabs_tmp endif enddo endif !----------------------------------------------------------------- ! Solve for new temperatures. ! Iterate until temperatures converge with minimal energy error. !----------------------------------------------------------------- converged = .false. do niter = 1, nitermax !----------------------------------------------------------------- ! Identify cells, if any, where calculation has not converged. !----------------------------------------------------------------- if (.not.converged) then !----------------------------------------------------------------- ! Allocate and initialize !----------------------------------------------------------------- converged = .true. dfsurf_dT = c0 avg_Tsi = c0 enew = c0 einex = c0 !----------------------------------------------------------------- ! Update specific heat of ice layers. ! To ensure energy conservation, the specific heat is a function of ! both the starting temperature and the (latest guess for) the ! final temperature. !----------------------------------------------------------------- do k = 1, nilyr if (l_brine) then ci = cp_ice - Lfresh*Tmlts(k) / & (zTin(k)*Tin_init(k)) else ci = cp_ice endif etai(k) = dt_rhoi_hlyr / ci enddo if (calc_Tsfc) then !----------------------------------------------------------------- ! Update radiative and turbulent fluxes and their derivatives ! with respect to Tsf. !----------------------------------------------------------------- ! surface heat flux call surface_heat_flux(Tsf , fswsfc, & rhoa , flw , & potT , Qa , & shcoef , lhcoef, & flwoutn, fsensn, & flatn , fsurfn) if (icepack_warnings_aborted(subname)) return ! derivative of heat flux with respect to surface temperature call dsurface_heat_flux_dTsf(Tsf , rhoa , & shcoef , lhcoef , & dfsurf_dT, dflwout_dT, & dfsens_dT, dflat_dT ) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! Compute conductive flux at top surface, fcondtopn. ! If fsurfn < fcondtopn and Tsf = 0, then reset Tsf to slightly less ! than zero (but not less than -puny). !----------------------------------------------------------------- if (l_snow) then fcondtopn = kh(1) * (Tsf - zTsn(1)) else fcondtopn = kh(1+nslyr) * (Tsf - zTin(1)) endif if (Tsf >= c0 .and. fsurfn < fcondtopn) & Tsf = -puny !----------------------------------------------------------------- ! Save surface temperature at start of iteration !----------------------------------------------------------------- Tsf_start = Tsf if (Tsf < c0) then l_cold = .true. else l_cold = .false. endif !----------------------------------------------------------------- ! Compute elements of tridiagonal matrix. !----------------------------------------------------------------- call get_matrix_elements_calc_Tsfc (nilyr, nslyr, & l_snow, l_cold, & Tsf, Tbot, & fsurfn, dfsurf_dT, & Tin_init, Tsn_init, & kh, Sswabs, & Iswabs, & etai, etas, & sbdiag, diag, & spdiag, rhs) if (icepack_warnings_aborted(subname)) return else call get_matrix_elements_know_Tsfc (nilyr, nslyr, & l_snow, Tbot, & Tin_init, Tsn_init, & kh, Sswabs, & Iswabs, & etai, etas, & sbdiag, diag, & spdiag, rhs, & fcondtopn) if (icepack_warnings_aborted(subname)) return endif ! calc_Tsfc !----------------------------------------------------------------- ! Solve tridiagonal matrix to obtain the new temperatures. !----------------------------------------------------------------- nmat = nslyr + nilyr + 1 ! matrix dimension call tridiag_solver (nmat, sbdiag(:), & diag(:), spdiag(:), & rhs(:), Tmat(:)) if (icepack_warnings_aborted(subname)) return !----------------------------------------------------------------- ! Determine whether the computation has converged to an acceptable ! solution. Five conditions must be satisfied: ! ! (1) Tsf <= 0 C. ! (2) Tsf is not oscillating; i.e., if both dTsf(niter) and ! dTsf(niter-1) have magnitudes greater than puny, then ! dTsf(niter)/dTsf(niter-1) cannot be a negative number ! with magnitude greater than 0.5. ! (3) abs(dTsf) < Tsf_errmax ! (4) If Tsf = 0 C, then the downward turbulent/radiative ! flux, fsurfn, must be greater than or equal to the downward ! conductive flux, fcondtopn. ! (5) The net energy added to the ice per unit time must equal ! the net change in internal ice energy per unit time, ! withinic the prescribed error ferrmax. ! ! For briny ice (the standard case), zTsn and zTin are limited ! to prevent them from exceeding their melting temperatures. ! (Note that the specific heat formula for briny ice assumes ! that T < Tmlt.) ! For fresh ice there is no limiting, since there are cases ! when the only convergent solution has zTsn > 0 and/or zTin > 0. ! Above-zero temperatures are then reset to zero (with melting ! to conserve energy) in the thickness_changes subroutine. !----------------------------------------------------------------- if (calc_Tsfc) then !----------------------------------------------------------------- ! Reload Tsf from matrix solution !----------------------------------------------------------------- if (l_cold) then if (l_snow) then Tsf = Tmat(1) else Tsf = Tmat(1+nslyr) endif else ! melting surface Tsf = c0 endif !----------------------------------------------------------------- ! Initialize convergence flag (true until proven false), dTsf, ! and temperature-averaging coefficients. ! Average only if test 1 or 2 fails. ! Initialize energy. !----------------------------------------------------------------- dTsf = Tsf - Tsf_start avg_Tsf = c0 !----------------------------------------------------------------- ! Condition 1: check for Tsf > 0 ! If Tsf > 0, set Tsf = 0, then average zTsn and zTin to force ! internal temps below their melting temps. !----------------------------------------------------------------- if (Tsf > puny) then Tsf = c0 dTsf = -Tsf_start if (l_brine) avg_Tsi = c1 ! avg with starting temp converged = .false. !----------------------------------------------------------------- ! Condition 2: check for oscillating Tsf ! If oscillating, average all temps to increase rate of convergence. !----------------------------------------------------------------- elseif (niter > 1 & ! condition (2) .and. Tsf_start <= -puny & .and. abs(dTsf) > puny & .and. abs(dTsf_prev) > puny & .and. -dTsf/(dTsf_prev+puny*puny) > p5) then if (l_brine) then ! average with starting temp avg_Tsf = c1 avg_Tsi = c1 endif dTsf = p5 * dTsf converged = .false. endif !!! dTsf_prev = dTsf !----------------------------------------------------------------- ! If condition 2 failed, average new surface temperature with ! starting value. !----------------------------------------------------------------- Tsf = Tsf & + avg_Tsf * p5 * (Tsf_start - Tsf) endif ! calc_Tsfc do k = 1, nslyr !----------------------------------------------------------------- ! Reload zTsn from matrix solution !----------------------------------------------------------------- if (l_snow) then zTsn(k) = Tmat(k+1) else zTsn(k) = c0 endif if (l_brine) zTsn(k) = min(zTsn(k), c0) !----------------------------------------------------------------- ! If condition 1 or 2 failed, average new snow layer ! temperatures with their starting values. !----------------------------------------------------------------- zTsn(k) = zTsn(k) & + avg_Tsi*p5*(Tsn_start(k)-zTsn(k)) !----------------------------------------------------------------- ! Compute zqsn and increment new energy. !----------------------------------------------------------------- zqsn(k) = -rhos * (Lfresh - cp_ice*zTsn(k)) enew = enew + hslyr * zqsn(k) Tsn_start(k) = zTsn(k) ! for next iteration enddo ! nslyr dTmat(:) = c0 dqmat(:) = c0 reduce_kh(:) = .false. do k = 1, nilyr !----------------------------------------------------------------- ! Reload zTin from matrix solution !----------------------------------------------------------------- zTin(k) = Tmat(k+1+nslyr) if (l_brine .and. zTin(k) > Tmlts(k) - puny) then dTmat(k) = zTin(k) - Tmlts(k) dqmat(k) = rhoi * dTmat(k) & * (cp_ice - Lfresh * Tmlts(k)/zTin(k)**2) ! use this for the case that Tmlt changes by an amount dTmlt=Tmltnew-Tmlt(k) ! + rhoi * dTmlt & ! * (cp_ocn - cp_ice + Lfresh/zTin(k)) zTin(k) = Tmlts(k) reduce_kh(k) = .true. endif !----------------------------------------------------------------- ! Condition 2b: check for oscillating zTin(1) ! If oscillating, average all ice temps to increase rate of convergence. !----------------------------------------------------------------- if (k==1 .and. .not.calc_Tsfc) then dTi1 = zTin(k) - Tin_start(k) if (niter > 1 & ! condition 2b .and. abs(dTi1) > puny & .and. abs(dTi1_prev) > puny & .and. -dTi1/(dTi1_prev+puny*puny) > p5) then if (l_brine) avg_Tsi = c1 dTi1 = p5 * dTi1 converged = .false. endif dTi1_prev = dTi1 endif ! k = 1 .and. calc_Tsfc = F !----------------------------------------------------------------- ! If condition 1 or 2 failed, average new ice layer ! temperatures with their starting values. !----------------------------------------------------------------- zTin(k) = zTin(k) & + avg_Tsi*p5*(Tin_start(k)-zTin(k)) !----------------------------------------------------------------- ! Compute zqin and increment new energy. !----------------------------------------------------------------- if (l_brine) then zqin(k) = -rhoi * (cp_ice*(Tmlts(k)-zTin(k)) & + Lfresh*(c1-Tmlts(k)/zTin(k)) & - cp_ocn*Tmlts(k)) else zqin(k) = -rhoi * (-cp_ice*zTin(k) + Lfresh) endif enew = enew + hilyr * zqin(k) einex = einex + hilyr * dqmat(k) Tin_start(k) = zTin(k) ! for next iteration enddo ! nilyr if (calc_Tsfc) then !----------------------------------------------------------------- ! Condition 3: check for large change in Tsf !----------------------------------------------------------------- if (abs(dTsf) > Tsf_errmax) then converged = .false. endif !----------------------------------------------------------------- ! Condition 4: check for fsurfn < fcondtopn with Tsf >= 0 !----------------------------------------------------------------- fsurfn = fsurfn + dTsf*dfsurf_dT if (l_snow) then fcondtopn = kh(1) * (Tsf-zTsn(1)) else fcondtopn = kh(1+nslyr) * (Tsf-zTin(1)) endif if (Tsf >= c0 .and. fsurfn < fcondtopn) then converged = .false. endif dTsf_prev = dTsf endif ! calc_Tsfc !----------------------------------------------------------------- ! Condition 5: check for energy conservation error ! Change in internal ice energy should equal net energy input. !----------------------------------------------------------------- fcondbot = kh(1+nslyr+nilyr) * & (zTin(nilyr) - Tbot) ! Flux extra energy out of the ice fcondbot = fcondbot + einex/dt ferr = abs( (enew-einit)/dt & - (fcondtopn - fcondbot + fswint) ) ! factor of 0.9 allows for roundoff errors later if (ferr > 0.9_dbl_kind*ferrmax) then ! condition (5) converged = .false. ! reduce conductivity for next iteration do k = 1, nilyr if (reduce_kh(k) .and. dqmat(k) > c0) then frac = max(0.5*(c1-ferr/abs(fcondtopn-fcondbot)),p1) ! frac = p1 kh(k+nslyr+1) = kh(k+nslyr+1) * frac kh(k+nslyr) = kh(k+nslyr+1) endif enddo endif ! ferr endif ! convergence enddo ! temperature iteration niter !----------------------------------------------------------------- ! Check for convergence failures. !----------------------------------------------------------------- if (.not.converged) then write(warnstr,*) subname, 'Thermo iteration does not converge,' call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Ice thickness:', hilyr*nilyr call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Snow thickness:', hslyr*nslyr call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'dTsf, Tsf_errmax:',dTsf_prev, & Tsf_errmax call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Tsf:', Tsf call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'fsurf:', fsurfn call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'fcondtop, fcondbot, fswint', & fcondtopn, fcondbot, fswint call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'fswsfc', fswsfc call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Iswabs',(Iswabs(k),k=1,nilyr) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Flux conservation error =', ferr call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Initial snow temperatures:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (Tsn_init(k),k=1,nslyr) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Initial ice temperatures:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (Tin_init(k),k=1,nilyr) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Matrix ice temperature diff:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (dTmat(k),k=1,nilyr) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'dqmat*hilyr/dt:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (hilyr*dqmat(k)/dt,k=1,nilyr) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Final snow temperatures:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (zTsn(k),k=1,nslyr) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Matrix ice temperature diff:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (dTmat(k),k=1,nilyr) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'dqmat*hilyr/dt:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (hilyr*dqmat(k)/dt,k=1,nilyr) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Final ice temperatures:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (zTin(k),k=1,nilyr) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Ice melting temperatures:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (Tmlts(k),k=1,nilyr) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'Ice bottom temperature:', Tbot call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'dT initial:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (Tmlts(k)-Tin_init(k),k=1,nilyr) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'dT final:' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (Tmlts(k)-zTin(k),k=1,nilyr) call icepack_warnings_add(warnstr) write(warnstr,*) subname, 'zSin' call icepack_warnings_add(warnstr) write(warnstr,*) subname, (zSin(k),k=1,nilyr) call icepack_warnings_add(warnstr) call icepack_warnings_setabort(.true.,__FILE__,__LINE__) call icepack_warnings_add(subname//" temperature_changes: Thermo iteration does not converge" ) return endif if (calc_Tsfc) then ! update fluxes that depend on Tsf flwoutn = flwoutn + dTsf_prev * dflwout_dT fsensn = fsensn + dTsf_prev * dfsens_dT flatn = flatn + dTsf_prev * dflat_dT endif ! calc_Tsfc end subroutine temperature_changes !======================================================================= ! ! Compute thermal conductivity at interfaces (held fixed during ! the subsequent iteration). ! ! NOTE: Ice conductivity must be >= kimin ! ! authors William H. Lipscomb, LANL ! C. M. Bitz, UW subroutine conductivity (l_snow, & nilyr, nslyr, & hilyr, hslyr, & zTin, kh, zSin) logical (kind=log_kind), intent(in) :: & l_snow ! true if snow temperatures are computed integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers nslyr ! number of snow layers real (kind=dbl_kind), intent(in) :: & hilyr , & ! ice layer thickness (same for all ice layers) hslyr ! snow layer thickness (same for all snow layers) real (kind=dbl_kind), dimension (:), intent(in) :: & zTin , & ! internal ice layer temperatures zSin ! internal ice layer salinities real (kind=dbl_kind), dimension (nilyr+nslyr+1), & intent(out) :: & kh ! effective conductivity at interfaces (W m-2 deg-1) ! local variables integer (kind=int_kind) :: & k ! vertical index real (kind=dbl_kind), dimension (nilyr) :: & kilyr ! thermal cond at ice layer midpoints (W m-1 deg-1) real (kind=dbl_kind), dimension (nslyr) :: & kslyr ! thermal cond at snow layer midpoints (W m-1 deg-1) character(len=*),parameter :: subname='(conductivity)' ! interior snow layers (simple for now, but may be fancier later) do k = 1, nslyr kslyr(k) = ksno enddo ! nslyr ! interior ice layers if (conduct == 'MU71') then ! Maykut and Untersteiner 1971 form (with Wettlaufer 1991 constants) do k = 1, nilyr kilyr(k) = kice + betak*zSin(k)/min(-puny,zTin(k)) kilyr(k) = max (kilyr(k), kimin) enddo ! nilyr else ! Pringle et al JGR 2007 'bubbly brine' do k = 1, nilyr kilyr(k) = (2.11_dbl_kind - 0.011_dbl_kind*zTin(k) & + 0.09_dbl_kind*zSin(k)/min(-puny,zTin(k))) & * rhoi / 917._dbl_kind kilyr(k) = max (kilyr(k), kimin) enddo ! nilyr endif ! conductivity ! top snow interface, top and bottom ice interfaces ! top of snow layer; top surface of top ice layer if (l_snow) then kh(1) = c2 * kslyr(1) / hslyr kh(1+nslyr) = c2 * kslyr(nslyr) * kilyr(1) / & ( kslyr(nslyr)*hilyr + & kilyr(1 )*hslyr ) else kh(1) = c0 kh(1+nslyr) = c2 * kilyr(1) / hilyr endif ! bottom surface of bottom ice layer kh(1+nslyr+nilyr) = c2 * kilyr(nilyr) / hilyr ! interior snow interfaces if (nslyr > 1) then do k = 2, nslyr if (l_snow) then kh(k) = c2 * kslyr(k-1) * kslyr(k) / & ((kslyr(k-1) + kslyr(k))*hslyr) else kh(k) = c0 endif enddo ! nilyr endif ! nslyr > 1 ! interior ice interfaces do k = 2, nilyr kh(k+nslyr) = c2 * kilyr(k-1) * kilyr(k) / & ((kilyr(k-1) + kilyr(k))*hilyr) enddo ! nilyr end subroutine conductivity !======================================================================= ! ! Compute radiative and turbulent fluxes and their derivatives ! with respect to Tsf. ! ! authors William H. Lipscomb, LANL ! C. M. Bitz, UW subroutine surface_fluxes (Tsf, fswsfc, & rhoa, flw, & potT, Qa, & shcoef, lhcoef, & flwoutn, fsensn, & flatn, fsurfn, & dflwout_dT, dfsens_dT, & dflat_dT, dfsurf_dT) real (kind=dbl_kind), intent(in) :: & Tsf ! ice/snow surface temperature, Tsfcn real (kind=dbl_kind), intent(in) :: & fswsfc , & ! SW absorbed at ice/snow surface (W m-2) rhoa , & ! air density (kg/m^3) flw , & ! incoming longwave radiation (W/m^2) potT , & ! air potential temperature (K) Qa , & ! specific humidity (kg/kg) shcoef , & ! transfer coefficient for sensible heat lhcoef ! transfer coefficient for latent heat real (kind=dbl_kind), & intent(inout) :: & fsensn , & ! surface downward sensible heat (W m-2) flatn , & ! surface downward latent heat (W m-2) flwoutn , & ! upward LW at surface (W m-2) fsurfn ! net flux to top surface, excluding fcondtopn real (kind=dbl_kind), & intent(inout) :: & dfsens_dT , & ! deriv of fsens wrt Tsf (W m-2 deg-1) dflat_dT , & ! deriv of flat wrt Tsf (W m-2 deg-1) dflwout_dT ! deriv of flwout wrt Tsf (W m-2 deg-1) real (kind=dbl_kind), & intent(inout) :: & dfsurf_dT ! derivative of fsurfn wrt Tsf character(len=*),parameter :: subname='(surface_fluxes)' ! surface heat flux call surface_heat_flux(Tsf, fswsfc, & rhoa, flw, & potT, Qa, & shcoef, lhcoef, & flwoutn, fsensn, & flatn, fsurfn) if (icepack_warnings_aborted(subname)) return ! derivative of heat flux with respect to surface temperature call dsurface_heat_flux_dTsf(Tsf, rhoa, & shcoef, lhcoef, & dfsurf_dT, dflwout_dT, & dfsens_dT, dflat_dT) if (icepack_warnings_aborted(subname)) return end subroutine surface_fluxes !======================================================================= ! ! Compute terms in tridiagonal matrix that will be solved to find ! the new vertical temperature profile ! This routine is for the case in which Tsfc is being computed. ! ! authors William H. Lipscomb, LANL ! C. M. Bitz, UW ! ! March 2004 by William H. Lipscomb for multiple snow layers ! April 2008 by E. C. Hunke, divided into two routines based on calc_Tsfc subroutine get_matrix_elements_calc_Tsfc (nilyr, nslyr, & l_snow, l_cold, & Tsf, Tbot, & fsurfn, dfsurf_dT, & Tin_init, Tsn_init, & kh, Sswabs, & Iswabs, & etai, etas, & sbdiag, diag, & spdiag, rhs) integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers nslyr ! number of snow layers logical (kind=log_kind), & intent(in) :: & l_snow , & ! true if snow temperatures are computed l_cold ! true if surface temperature is computed real (kind=dbl_kind), intent(in) :: & Tsf ! ice/snow top surface temp (deg C) real (kind=dbl_kind), intent(in) :: & fsurfn , & ! net flux to top surface, excluding fcondtopn (W/m^2) Tbot ! ice bottom surface temperature (deg C) real (kind=dbl_kind), intent(in) :: & dfsurf_dT ! derivative of fsurf wrt Tsf real (kind=dbl_kind), dimension (:), intent(in) :: & etai , & ! dt / (rho*cp*h) for ice layers Tin_init , & ! ice temp at beginning of time step Sswabs , & ! SW radiation absorbed in snow layers (W m-2) Iswabs , & ! absorbed SW flux in ice layers etas , & ! dt / (rho*cp*h) for snow layers Tsn_init ! snow temp at beginning of time step ! Note: no absorbed SW in snow layers real (kind=dbl_kind), dimension (nslyr+nilyr+1), & intent(in) :: & kh ! effective conductivity at layer interfaces real (kind=dbl_kind), dimension (nslyr+nilyr+1), & intent(inout) :: & sbdiag , & ! sub-diagonal matrix elements diag , & ! diagonal matrix elements spdiag , & ! super-diagonal matrix elements rhs ! rhs of tri-diagonal matrix eqn. ! local variables integer (kind=int_kind) :: & k, ki, kr ! vertical indices and row counters character(len=*),parameter :: subname='(get_matrix_elements_calc_Tsrf)' !----------------------------------------------------------------- ! Initialize matrix elements. ! Note: When we do not need to solve for the surface or snow ! temperature, we solve dummy equations with solution T = 0. ! Ice layers are fully initialized below. !----------------------------------------------------------------- do k = 1, nslyr+1 sbdiag(k) = c0 diag (k) = c1 spdiag(k) = c0 rhs (k) = c0 enddo !----------------------------------------------------------------- ! Compute matrix elements ! ! Four possible cases to solve: ! (1) Cold surface (Tsf < 0), snow present ! (2) Melting surface (Tsf = 0), snow present ! (3) Cold surface (Tsf < 0), no snow ! (4) Melting surface (Tsf = 0), no snow !----------------------------------------------------------------- !----------------------------------------------------------------- ! Tsf equation for case of cold surface (with or without snow) !----------------------------------------------------------------- if (l_cold) then if (l_snow) then k = 1 else ! no snow k = 1 + nslyr endif kr = k sbdiag(kr) = c0 diag (kr) = dfsurf_dT - kh(k) spdiag(kr) = kh(k) rhs (kr) = dfsurf_dT*Tsf - fsurfn endif ! l_cold !----------------------------------------------------------------- ! top snow layer !----------------------------------------------------------------- ! k = 1 ! kr = 2 if (l_snow) then if (l_cold) then sbdiag(2) = -etas(1) * kh(1) spdiag(2) = -etas(1) * kh(2) diag (2) = c1 & + etas(1) * (kh(1) + kh(2)) rhs (2) = Tsn_init(1) & + etas(1) * Sswabs(1) else ! melting surface sbdiag(2) = c0 spdiag(2) = -etas(1) * kh(2) diag (2) = c1 & + etas(1) * (kh(1) + kh(2)) rhs (2) = Tsn_init(1) & + etas(1)*kh(1)*Tsf & + etas(1) * Sswabs(1) endif ! l_cold endif ! l_snow !----------------------------------------------------------------- ! remaining snow layers !----------------------------------------------------------------- if (nslyr > 1) then do k = 2, nslyr kr = k + 1 if (l_snow) then sbdiag(kr) = -etas(k) * kh(k) spdiag(kr) = -etas(k) * kh(k+1) diag (kr) = c1 & + etas(k) * (kh(k) + kh(k+1)) rhs (kr) = Tsn_init(k) & + etas(k) * Sswabs(k) endif enddo ! nslyr endif ! nslyr > 1 if (nilyr > 1) then !----------------------------------------------------------------- ! top ice layer !----------------------------------------------------------------- ki = 1 k = ki + nslyr kr = k + 1 if (l_snow .or. l_cold) then sbdiag(kr) = -etai(ki) * kh(k) spdiag(kr) = -etai(ki) * kh(k+1) diag (kr) = c1 & + etai(ki) * (kh(k) + kh(k+1)) rhs (kr) = Tin_init(ki) & + etai(ki)*Iswabs(ki) else ! no snow, warm surface sbdiag(kr) = c0 spdiag(kr) = -etai(ki) * kh(k+1) diag (kr) = c1 & + etai(ki) * (kh(k) + kh(k+1)) rhs (kr) = Tin_init(ki) & + etai(ki)*Iswabs(ki) & + etai(ki)*kh(k)*Tsf endif !----------------------------------------------------------------- ! bottom ice layer !----------------------------------------------------------------- ki = nilyr k = ki + nslyr kr = k + 1 sbdiag(kr) = -etai(ki) * kh(k) spdiag(kr) = c0 diag (kr) = c1 & + etai(ki) * (kh(k) + kh(k+1)) rhs (kr) = Tin_init(ki) & + etai(ki)*Iswabs(ki) & + etai(ki)*kh(k+1)*Tbot else ! nilyr = 1 !----------------------------------------------------------------- ! single ice layer !----------------------------------------------------------------- ki = 1 k = ki + nslyr kr = k + 1 if (l_snow .or. l_cold) then sbdiag(kr) = -etai(ki) * kh(k) spdiag(kr) = c0 diag (kr) = c1 & + etai(ki) * (kh(k) + kh(k+1)) rhs (kr) = Tin_init(ki) & + etai(ki) * Iswabs(ki) & + etai(ki) * kh(k+1)*Tbot else ! no snow, warm surface sbdiag(kr) = c0 spdiag(kr) = c0 diag (kr) = c1 & + etai(ki) * (kh(k) + kh(k+1)) rhs (kr) = Tin_init(ki) & + etai(ki) * Iswabs(ki) & + etai(ki) * kh(k)*Tsf & + etai(ki) * kh(k+1)*Tbot endif endif ! nilyr > 1 !----------------------------------------------------------------- ! interior ice layers !----------------------------------------------------------------- do ki = 2, nilyr-1 k = ki + nslyr kr = k + 1 sbdiag(kr) = -etai(ki) * kh(k) spdiag(kr) = -etai(ki) * kh(k+1) diag (kr) = c1 & + etai(ki) * (kh(k) + kh(k+1)) rhs (kr) = Tin_init(ki) & + etai(ki)*Iswabs(ki) enddo ! nilyr end subroutine get_matrix_elements_calc_Tsfc !======================================================================= ! ! Compute terms in tridiagonal matrix that will be solved to find ! the new vertical temperature profile ! This routine is for the case in which Tsfc is already known. ! ! authors William H. Lipscomb, LANL ! C. M. Bitz, UW ! ! March 2004 by William H. Lipscomb for multiple snow layers ! April 2008 by E. C. Hunke, divided into two routines based on calc_Tsfc subroutine get_matrix_elements_know_Tsfc (nilyr, nslyr, & l_snow, Tbot, & Tin_init, Tsn_init, & kh, Sswabs, & Iswabs, & etai, etas, & sbdiag, diag, & spdiag, rhs, & fcondtopn) integer (kind=int_kind), intent(in) :: & nilyr , & ! number of ice layers nslyr ! number of snow layers logical (kind=log_kind), & intent(in) :: & l_snow ! true if snow temperatures are computed real (kind=dbl_kind), intent(in) :: & Tbot ! ice bottom surface temperature (deg C) real (kind=dbl_kind), dimension (:), intent(in) :: & etai , & ! dt / (rho*cp*h) for ice layers Tin_init , & ! ice temp at beginning of time step Sswabs , & ! SW radiation absorbed in snow layers (W m-2) Iswabs , & ! absorbed SW flux in ice layers etas , & ! dt / (rho*cp*h) for snow layers Tsn_init ! snow temp at beginning of time step ! Note: no absorbed SW in snow layers real (kind=dbl_kind), dimension (nslyr+nilyr+1), & intent(in) :: & kh ! effective conductivity at layer interfaces real (kind=dbl_kind), dimension (nslyr+nilyr+1), & intent(inout) :: & sbdiag , & ! sub-diagonal matrix elements diag , & ! diagonal matrix elements spdiag , & ! super-diagonal matrix elements rhs ! rhs of tri-diagonal matrix eqn. real (kind=dbl_kind), intent(in), & optional :: & fcondtopn ! conductive flux at top sfc, positive down (W/m^2) ! local variables integer (kind=int_kind) :: & k, ki, kr ! vertical indices and row counters character(len=*),parameter :: subname='(get_matrix_elements_know_Tsrf)' !----------------------------------------------------------------- ! Initialize matrix elements. ! Note: When we do not need to solve for the surface or snow ! temperature, we solve dummy equations with solution T = 0. ! Ice layers are fully initialized below. !----------------------------------------------------------------- do k = 1, nslyr+1 sbdiag(k) = c0 diag (k) = c1 spdiag(k) = c0 rhs (k) = c0 enddo !----------------------------------------------------------------- ! Compute matrix elements ! ! Four possible cases to solve: ! (1) Cold surface (Tsf < 0), snow present ! (2) Melting surface (Tsf = 0), snow present ! (3) Cold surface (Tsf < 0), no snow ! (4) Melting surface (Tsf = 0), no snow !----------------------------------------------------------------- !----------------------------------------------------------------- ! top snow layer !----------------------------------------------------------------- ! k = 1 ! kr = 2 if (l_snow) then sbdiag(2) = c0 spdiag(2) = -etas(1) * kh(2) diag (2) = c1 & + etas(1) * kh(2) rhs (2) = Tsn_init(1) & + etas(1) * Sswabs(1) & + etas(1) * fcondtopn endif ! l_snow !----------------------------------------------------------------- ! remaining snow layers !----------------------------------------------------------------- if (nslyr > 1) then do k = 2, nslyr kr = k + 1 if (l_snow) then sbdiag(kr) = -etas(k) * kh(k) spdiag(kr) = -etas(k) * kh(k+1) diag (kr) = c1 & + etas(k) * (kh(k) + kh(k+1)) rhs (kr) = Tsn_init(k) & + etas(k) * Sswabs(k) endif enddo ! nslyr endif ! nslyr > 1 if (nilyr > 1) then !----------------------------------------------------------------- ! top ice layer !----------------------------------------------------------------- ki = 1 k = ki + nslyr kr = k + 1 if (l_snow) then sbdiag(kr) = -etai(ki) * kh(k) spdiag(kr) = -etai(ki) * kh(k+1) diag (kr) = c1 & + etai(ki) * (kh(k) + kh(k+1)) rhs (kr) = Tin_init(ki) & + etai(ki) * Iswabs(ki) else sbdiag(kr) = c0 spdiag(kr) = -etai(ki) * kh(k+1) diag (kr) = c1 & + etai(ki) * kh(k+1) rhs (kr) = Tin_init(ki) & + etai(ki) * Iswabs(ki) & + etai(ki) * fcondtopn endif ! l_snow !----------------------------------------------------------------- ! bottom ice layer !----------------------------------------------------------------- ki = nilyr k = ki + nslyr kr = k + 1 sbdiag(kr) = -etai(ki) * kh(k) spdiag(kr) = c0 diag (kr) = c1 & + etai(ki) * (kh(k) + kh(k+1)) rhs (kr) = Tin_init(ki) & + etai(ki)*Iswabs(ki) & + etai(ki)*kh(k+1)*Tbot else ! nilyr = 1 !----------------------------------------------------------------- ! single ice layer !----------------------------------------------------------------- ki = 1 k = ki + nslyr kr = k + 1 if (l_snow) then sbdiag(kr) = -etai(ki) * kh(k) spdiag(kr) = c0 diag (kr) = c1 & + etai(ki) * (kh(k) + kh(k+1)) rhs (kr) = Tin_init(ki) & + etai(ki) * Iswabs(ki) & + etai(ki) * kh(k+1)*Tbot else sbdiag(kr) = c0 spdiag(kr) = c0 diag (kr) = c1 & + etai(ki) * kh(k+1) rhs (kr) = Tin_init(ki) & + etai(ki) * Iswabs(ki) & + etai(ki) * fcondtopn & + etai(ki) * kh(k+1)*Tbot endif endif ! nilyr > 1 !----------------------------------------------------------------- ! interior ice layers !----------------------------------------------------------------- do ki = 2, nilyr-1 k = ki + nslyr kr = k + 1 sbdiag(kr) = -etai(ki) * kh(k) spdiag(kr) = -etai(ki) * kh(k+1) diag (kr) = c1 & + etai(ki) * (kh(k) + kh(k+1)) rhs (kr) = Tin_init(ki) & + etai(ki)*Iswabs(ki) enddo ! nilyr end subroutine get_matrix_elements_know_Tsfc !======================================================================= ! ! Tridiagonal matrix solver--used to solve the implicit vertical heat ! equation in ice and snow ! ! authors William H. Lipscomb, LANL ! C. M. Bitz, UW subroutine tridiag_solver (nmat, sbdiag, & diag, spdiag, & rhs, xout) integer (kind=int_kind), intent(in) :: & nmat ! matrix dimension real (kind=dbl_kind), dimension (:), intent(in) :: & sbdiag , & ! sub-diagonal matrix elements diag , & ! diagonal matrix elements spdiag , & ! super-diagonal matrix elements rhs ! rhs of tri-diagonal matrix eqn. real (kind=dbl_kind), dimension (:), intent(inout) :: & xout ! solution vector ! local variables integer (kind=int_kind) :: & k ! row counter real (kind=dbl_kind) :: & wbeta ! temporary matrix variable real (kind=dbl_kind), dimension(nmat) :: & wgamma ! temporary matrix variable character(len=*),parameter :: subname='(tridiag_solver)' wbeta = diag(1) xout(1) = rhs(1) / wbeta do k = 2, nmat wgamma(k) = spdiag(k-1) / wbeta wbeta = diag(k) - sbdiag(k)*wgamma(k) xout(k) = (rhs(k) - sbdiag(k)*xout(k-1)) & / wbeta enddo ! k do k = nmat-1, 1, -1 xout(k) = xout(k) - wgamma(k+1)*xout(k+1) enddo ! k end subroutine tridiag_solver !======================================================================= end module icepack_therm_bl99 !=======================================================================