Line data Source code
1 : !=========================================================================
2 : !
3 : ! Update ice and snow internal temperatures
4 : ! using Bitz and Lipscomb 1999 thermodynamics
5 : !
6 : ! authors: William H. Lipscomb, LANL
7 : ! C. M. Bitz, UW
8 : ! Elizabeth C. Hunke, LANL
9 : !
10 : ! 2012: Split from icepack_therm_vertical.F90
11 :
12 : module icepack_therm_bl99
13 :
14 : use icepack_kinds
15 : use icepack_parameters, only: c0, c1, c2, p1, p5, puny
16 : #ifdef CESMCOUPLED
17 : use icepack_parameters, only: p01
18 : #endif
19 : use icepack_parameters, only: rhoi, rhos, hs_min, cp_ice, cp_ocn, depressT, Lfresh, ksno, kice
20 : use icepack_parameters, only: conduct, calc_Tsfc, solve_zsal
21 : use icepack_parameters, only: sw_redist, sw_frac, sw_dtemp
22 : use icepack_warnings, only: warnstr, icepack_warnings_add
23 : use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
24 :
25 : use icepack_therm_shared, only: ferrmax, l_brine
26 : use icepack_therm_shared, only: surface_heat_flux, dsurface_heat_flux_dTsf
27 :
28 : implicit none
29 :
30 : private
31 : public :: surface_fluxes, temperature_changes
32 :
33 : real (kind=dbl_kind), parameter :: &
34 : betak = 0.13_dbl_kind, & ! constant in formula for k (W m-1 ppt-1) ! LCOV_EXCL_LINE
35 : kimin = 0.10_dbl_kind ! min conductivity of saline ice (W m-1 deg-1)
36 :
37 : !=======================================================================
38 :
39 : contains
40 :
41 : !=======================================================================
42 : !
43 : ! Compute new surface temperature and internal ice and snow
44 : ! temperatures. Include effects of salinity on sea ice heat
45 : ! capacity in a way that conserves energy (Bitz and Lipscomb, 1999).
46 : !
47 : ! New temperatures are computed iteratively by solving a tridiagonal
48 : ! system of equations; heat capacity is updated with each iteration.
49 : ! Finite differencing is backward implicit.
50 : !
51 : ! See Bitz, C.M., and W.H. Lipscomb, 1999:
52 : ! An energy-conserving thermodynamic model of sea ice,
53 : ! J. Geophys. Res., 104, 15,669-15,677.
54 : !
55 : ! authors William H. Lipscomb, LANL
56 : ! C. M. Bitz, UW
57 :
58 59683339 : subroutine temperature_changes (dt, &
59 : nilyr, nslyr, & ! LCOV_EXCL_LINE
60 : rhoa, flw, & ! LCOV_EXCL_LINE
61 : potT, Qa, & ! LCOV_EXCL_LINE
62 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
63 : fswsfc, fswint, & ! LCOV_EXCL_LINE
64 : Sswabs, Iswabs, & ! LCOV_EXCL_LINE
65 : hilyr, hslyr, & ! LCOV_EXCL_LINE
66 : zqin, zTin, & ! LCOV_EXCL_LINE
67 : zqsn, zTsn, & ! LCOV_EXCL_LINE
68 : zSin, & ! LCOV_EXCL_LINE
69 : Tsf, Tbot, & ! LCOV_EXCL_LINE
70 : fsensn, flatn, & ! LCOV_EXCL_LINE
71 : flwoutn, fsurfn, & ! LCOV_EXCL_LINE
72 : fcondtopn,fcondbot, & ! LCOV_EXCL_LINE
73 : einit )
74 :
75 : integer (kind=int_kind), intent(in) :: &
76 : nilyr , & ! number of ice layers ! LCOV_EXCL_LINE
77 : nslyr ! number of snow layers
78 :
79 : real (kind=dbl_kind), intent(in) :: &
80 : dt ! time step
81 :
82 : real (kind=dbl_kind), &
83 : intent(in) :: & ! LCOV_EXCL_LINE
84 : rhoa , & ! air density (kg/m^3) ! LCOV_EXCL_LINE
85 : flw , & ! incoming longwave radiation (W/m^2) ! LCOV_EXCL_LINE
86 : potT , & ! air potential temperature (K) ! LCOV_EXCL_LINE
87 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
88 : shcoef , & ! transfer coefficient for sensible heat ! LCOV_EXCL_LINE
89 : lhcoef , & ! transfer coefficient for latent heat ! LCOV_EXCL_LINE
90 : Tbot ! ice bottom surface temperature (deg C)
91 :
92 : real (kind=dbl_kind), &
93 : intent(inout) :: & ! LCOV_EXCL_LINE
94 : fswsfc , & ! SW absorbed at ice/snow surface (W m-2) ! LCOV_EXCL_LINE
95 : fswint ! SW absorbed in ice interior below surface (W m-2)
96 :
97 : real (kind=dbl_kind), intent(in) :: &
98 : hilyr , & ! ice layer thickness (m) ! LCOV_EXCL_LINE
99 : hslyr , & ! snow layer thickness (m) ! LCOV_EXCL_LINE
100 : einit ! initial energy of melting (J m-2)
101 :
102 : real (kind=dbl_kind), dimension (nslyr), &
103 : intent(inout) :: & ! LCOV_EXCL_LINE
104 : Sswabs ! SW radiation absorbed in snow layers (W m-2)
105 :
106 : real (kind=dbl_kind), dimension (nilyr), &
107 : intent(inout) :: & ! LCOV_EXCL_LINE
108 : Iswabs ! SW radiation absorbed in ice layers (W m-2)
109 :
110 : real (kind=dbl_kind), intent(inout):: &
111 : fsurfn , & ! net flux to top surface, excluding fcondtopn ! LCOV_EXCL_LINE
112 : fcondtopn , & ! downward cond flux at top surface (W m-2) ! LCOV_EXCL_LINE
113 : fsensn , & ! surface downward sensible heat (W m-2) ! LCOV_EXCL_LINE
114 : flatn , & ! surface downward latent heat (W m-2) ! LCOV_EXCL_LINE
115 : flwoutn ! upward LW at surface (W m-2)
116 :
117 : real (kind=dbl_kind), intent(out):: &
118 : fcondbot ! downward cond flux at bottom surface (W m-2)
119 :
120 : real (kind=dbl_kind), &
121 : intent(inout):: & ! LCOV_EXCL_LINE
122 : Tsf ! ice/snow surface temperature, Tsfcn
123 :
124 : real (kind=dbl_kind), dimension (nilyr), &
125 : intent(inout) :: & ! LCOV_EXCL_LINE
126 : zqin , & ! ice layer enthalpy (J m-3) ! LCOV_EXCL_LINE
127 : zTin ! internal ice layer temperatures
128 :
129 : real (kind=dbl_kind), dimension (nilyr), &
130 : intent(in) :: & ! LCOV_EXCL_LINE
131 : zSin ! internal ice layer salinities
132 :
133 : real (kind=dbl_kind), dimension (nslyr), &
134 : intent(inout) :: & ! LCOV_EXCL_LINE
135 : zqsn , & ! snow layer enthalpy (J m-3) ! LCOV_EXCL_LINE
136 : zTsn ! internal snow layer temperatures
137 :
138 : ! local variables
139 :
140 : integer (kind=int_kind), parameter :: &
141 : nitermax = 100 ! max number of iterations in temperature solver
142 :
143 : real (kind=dbl_kind), parameter :: &
144 : Tsf_errmax = 5.e-4_dbl_kind ! max allowed error in Tsf
145 : ! recommend Tsf_errmax < 0.01 K
146 :
147 : integer (kind=int_kind) :: &
148 : k , & ! ice layer index ! LCOV_EXCL_LINE
149 : niter , & ! iteration counter in temperature solver ! LCOV_EXCL_LINE
150 : nmat ! matrix dimension
151 :
152 : logical (kind=log_kind) :: &
153 : l_snow , & ! true if snow temperatures are computed ! LCOV_EXCL_LINE
154 : l_cold ! true if surface temperature is computed
155 :
156 : real (kind=dbl_kind) :: &
157 : Tsf_start , & ! Tsf at start of iteration ! LCOV_EXCL_LINE
158 : dTsf , & ! Tsf - Tsf_start ! LCOV_EXCL_LINE
159 : dTi1 , & ! Ti1(1) - Tin_start(1) ! LCOV_EXCL_LINE
160 : dfsurf_dT , & ! derivative of fsurf wrt Tsf ! LCOV_EXCL_LINE
161 : avg_Tsi , & ! = 1. if new snow/ice temps averaged w/starting temps ! LCOV_EXCL_LINE
162 5427462 : enew ! new energy of melting after temp change (J m-2)
163 :
164 : real (kind=dbl_kind) :: &
165 : dTsf_prev , & ! dTsf from previous iteration ! LCOV_EXCL_LINE
166 : dTi1_prev , & ! dTi1 from previous iteration ! LCOV_EXCL_LINE
167 : dfsens_dT , & ! deriv of fsens wrt Tsf (W m-2 deg-1) ! LCOV_EXCL_LINE
168 : dflat_dT , & ! deriv of flat wrt Tsf (W m-2 deg-1) ! LCOV_EXCL_LINE
169 : dflwout_dT , & ! deriv of flwout wrt Tsf (W m-2 deg-1) ! LCOV_EXCL_LINE
170 : dt_rhoi_hlyr, & ! dt/(rhoi*hilyr) ! LCOV_EXCL_LINE
171 : einex , & ! excess energy from dqmat to ocean ! LCOV_EXCL_LINE
172 5427462 : ferr ! energy conservation error (W m-2)
173 :
174 : real (kind=dbl_kind), dimension (nilyr) :: &
175 : Tin_init , & ! zTin at beginning of time step ! LCOV_EXCL_LINE
176 : Tin_start , & ! zTin at start of iteration ! LCOV_EXCL_LINE
177 : dTmat , & ! zTin - matrix solution before limiting ! LCOV_EXCL_LINE
178 : dqmat , & ! associated enthalpy difference ! LCOV_EXCL_LINE
179 151931450 : Tmlts ! melting temp, -depressT * salinity
180 :
181 : real (kind=dbl_kind), dimension (nslyr) :: &
182 : Tsn_init , & ! zTsn at beginning of time step ! LCOV_EXCL_LINE
183 : Tsn_start , & ! zTsn at start of iteration ! LCOV_EXCL_LINE
184 119366678 : etas ! dt / (rho * cp * h) for snow layers
185 :
186 : real (kind=dbl_kind), dimension (nilyr+nslyr+1) :: &
187 : sbdiag , & ! sub-diagonal matrix elements ! LCOV_EXCL_LINE
188 : diag , & ! diagonal matrix elements ! LCOV_EXCL_LINE
189 : spdiag , & ! super-diagonal matrix elements ! LCOV_EXCL_LINE
190 : rhs , & ! rhs of tri-diagonal matrix equation ! LCOV_EXCL_LINE
191 162786374 : Tmat ! matrix output temperatures
192 :
193 : real (kind=dbl_kind), dimension (nilyr) :: &
194 151931450 : etai ! dt / (rho * cp * h) for ice layers
195 :
196 : real (kind=dbl_kind), dimension(nilyr+nslyr+1):: &
197 162786374 : kh ! effective conductivity at interfaces (W m-2 deg-1)
198 :
199 : real (kind=dbl_kind) :: &
200 : ci , & ! specific heat of sea ice (J kg-1 deg-1) ! LCOV_EXCL_LINE
201 : avg_Tsf , & ! = 1. if Tsf averaged w/Tsf_start, else = 0. ! LCOV_EXCL_LINE
202 : Iswabs_tmp , & ! energy to melt through fraction frac of layer ! LCOV_EXCL_LINE
203 : Sswabs_tmp , & ! same for snow ! LCOV_EXCL_LINE
204 : dswabs , & ! difference in swabs and swabs_tmp ! LCOV_EXCL_LINE
205 5427462 : frac
206 :
207 : logical (kind=log_kind) :: &
208 : converged ! = true when local solution has converged
209 :
210 : logical (kind=log_kind) , dimension (nilyr) :: &
211 65110801 : reduce_kh ! reduce conductivity when T exceeds Tmlt
212 :
213 : character(len=*),parameter :: subname='(temperature_changes)'
214 :
215 : !-----------------------------------------------------------------
216 : ! Initialize
217 : !-----------------------------------------------------------------
218 :
219 59683339 : converged = .false.
220 59683339 : l_snow = .false.
221 59683339 : l_cold = .true.
222 59683339 : fcondbot = c0
223 59683339 : dTsf_prev = c0
224 59683339 : dTi1_prev = c0
225 59683339 : dfsens_dT = c0
226 59683339 : dflat_dT = c0
227 59683339 : dflwout_dT = c0
228 59683339 : einex = c0
229 59683339 : dt_rhoi_hlyr = dt / (rhoi*hilyr) ! hilyr > 0
230 59683339 : if (hslyr > hs_min/real(nslyr,kind=dbl_kind)) &
231 56027048 : l_snow = .true.
232 :
233 119366678 : do k = 1, nslyr
234 59683339 : Tsn_init (k) = zTsn(k) ! beginning of time step
235 59683339 : Tsn_start(k) = zTsn(k) ! beginning of iteration
236 119366678 : if (l_snow) then
237 56027048 : etas(k) = dt/(rhos*cp_ice*hslyr)
238 : else
239 3656291 : etas(k) = c0
240 : endif
241 : enddo ! k
242 :
243 477466712 : do k = 1, nilyr
244 417783373 : Tin_init (k) = zTin(k) ! beginning of time step
245 417783373 : Tin_start(k) = zTin(k) ! beginning of iteration
246 477466712 : Tmlts (k) = -zSin(k) * depressT
247 : enddo
248 :
249 : !-----------------------------------------------------------------
250 : ! Compute thermal conductivity at interfaces (held fixed during
251 : ! subsequent iterations).
252 : ! Ice and snow interfaces are combined into one array (kh) to
253 : ! simplify the logic.
254 : !-----------------------------------------------------------------
255 :
256 : call conductivity (l_snow, &
257 : nilyr, nslyr, & ! LCOV_EXCL_LINE
258 : hilyr, hslyr, & ! LCOV_EXCL_LINE
259 59683339 : zTin, kh, zSin)
260 59683339 : if (icepack_warnings_aborted(subname)) return
261 :
262 : !-----------------------------------------------------------------
263 : ! Check for excessive absorbed solar radiation that may result in
264 : ! temperature overshoots. Convergence is particularly difficult
265 : ! if the starting temperature is already very close to the melting
266 : ! temperature and extra energy is added. In that case, or if the
267 : ! amount of energy absorbed is greater than the amount needed to
268 : ! melt through a given fraction of a layer, we put the extra
269 : ! energy into the surface.
270 : ! NOTE: This option is not available if the atmosphere model
271 : ! has already computed fsurf. (Unless we adjust fsurf here)
272 : !-----------------------------------------------------------------
273 : !mclaren: Should there be an if calc_Tsfc statement here then??
274 :
275 59683339 : if (sw_redist) then
276 :
277 59683339 : if (solve_zsal) sw_dtemp = p1 ! lower tolerance with dynamic salinity
278 :
279 477466712 : do k = 1, nilyr
280 :
281 417783373 : Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc
282 417783373 : if (Tin_init(k) <= Tmlts(k) - sw_dtemp) then
283 417677887 : if (l_brine) then
284 417677887 : ci = cp_ice - Lfresh * Tmlts(k) / (Tin_init(k)**2)
285 37990022 : Iswabs_tmp = min(Iswabs(k), &
286 417677887 : sw_frac*(Tmlts(k)-Tin_init(k))*ci/dt_rhoi_hlyr)
287 : else
288 0 : ci = cp_ice
289 0 : Iswabs_tmp = min(Iswabs(k), &
290 0 : sw_frac*(-Tin_init(k))*ci/dt_rhoi_hlyr)
291 : endif
292 : endif
293 417783373 : if (Iswabs_tmp < puny) Iswabs_tmp = c0
294 :
295 417783373 : dswabs = min(Iswabs(k) - Iswabs_tmp, fswint)
296 :
297 417783373 : fswsfc = fswsfc + dswabs
298 417783373 : fswint = fswint - dswabs
299 477466712 : Iswabs(k) = Iswabs_tmp
300 :
301 : enddo
302 :
303 119366678 : do k = 1, nslyr
304 119366678 : if (l_snow) then
305 :
306 56027048 : Sswabs_tmp = c0
307 56027048 : if (Tsn_init(k) <= -sw_dtemp) then
308 5261928 : Sswabs_tmp = min(Sswabs(k), &
309 55940960 : -sw_frac*Tsn_init(k)/etas(k))
310 : endif
311 56027048 : if (Sswabs_tmp < puny) Sswabs_tmp = c0
312 :
313 56027048 : dswabs = min(Sswabs(k) - Sswabs_tmp, fswint)
314 :
315 56027048 : fswsfc = fswsfc + dswabs
316 56027048 : fswint = fswint - dswabs
317 56027048 : Sswabs(k) = Sswabs_tmp
318 :
319 : endif
320 : enddo
321 :
322 : endif
323 :
324 : !-----------------------------------------------------------------
325 : ! Solve for new temperatures.
326 : ! Iterate until temperatures converge with minimal energy error.
327 : !-----------------------------------------------------------------
328 59683339 : converged = .false.
329 :
330 6028017239 : do niter = 1, nitermax
331 :
332 : !-----------------------------------------------------------------
333 : ! Identify cells, if any, where calculation has not converged.
334 : !-----------------------------------------------------------------
335 :
336 6028017239 : if (.not.converged) then
337 :
338 : !-----------------------------------------------------------------
339 : ! Allocate and initialize
340 : !-----------------------------------------------------------------
341 :
342 134079170 : converged = .true.
343 134079170 : dfsurf_dT = c0
344 134079170 : avg_Tsi = c0
345 134079170 : enew = c0
346 134079170 : einex = c0
347 :
348 : !-----------------------------------------------------------------
349 : ! Update specific heat of ice layers.
350 : ! To ensure energy conservation, the specific heat is a function of
351 : ! both the starting temperature and the (latest guess for) the
352 : ! final temperature.
353 : !-----------------------------------------------------------------
354 :
355 1072633360 : do k = 1, nilyr
356 :
357 938554190 : if (l_brine) then
358 78299613 : ci = cp_ice - Lfresh*Tmlts(k) / &
359 938554190 : (zTin(k)*Tin_init(k))
360 : else
361 0 : ci = cp_ice
362 : endif
363 1072633360 : etai(k) = dt_rhoi_hlyr / ci
364 :
365 : enddo
366 :
367 134079170 : if (calc_Tsfc) then
368 :
369 : !-----------------------------------------------------------------
370 : ! Update radiative and turbulent fluxes and their derivatives
371 : ! with respect to Tsf.
372 : !-----------------------------------------------------------------
373 :
374 : ! surface heat flux
375 : call surface_heat_flux(Tsf , fswsfc, &
376 : rhoa , flw , & ! LCOV_EXCL_LINE
377 : potT , Qa , & ! LCOV_EXCL_LINE
378 : shcoef , lhcoef, & ! LCOV_EXCL_LINE
379 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
380 74706776 : flatn , fsurfn)
381 74706776 : if (icepack_warnings_aborted(subname)) return
382 :
383 : ! derivative of heat flux with respect to surface temperature
384 : call dsurface_heat_flux_dTsf(Tsf , rhoa , &
385 : shcoef , lhcoef , & ! LCOV_EXCL_LINE
386 : dfsurf_dT, dflwout_dT, & ! LCOV_EXCL_LINE
387 74706776 : dfsens_dT, dflat_dT )
388 74706776 : if (icepack_warnings_aborted(subname)) return
389 :
390 : !-----------------------------------------------------------------
391 : ! Compute conductive flux at top surface, fcondtopn.
392 : ! If fsurfn < fcondtopn and Tsf = 0, then reset Tsf to slightly less
393 : ! than zero (but not less than -puny).
394 : !-----------------------------------------------------------------
395 :
396 74706776 : if (l_snow) then
397 64507405 : fcondtopn = kh(1) * (Tsf - zTsn(1))
398 : else
399 10199371 : fcondtopn = kh(1+nslyr) * (Tsf - zTin(1))
400 : endif
401 :
402 74706776 : if (Tsf >= c0 .and. fsurfn < fcondtopn) &
403 1123256 : Tsf = -puny
404 :
405 : !-----------------------------------------------------------------
406 : ! Save surface temperature at start of iteration
407 : !-----------------------------------------------------------------
408 :
409 74706776 : Tsf_start = Tsf
410 :
411 74706776 : if (Tsf < c0) then
412 54336220 : l_cold = .true.
413 : else
414 20370556 : l_cold = .false.
415 : endif
416 :
417 : !-----------------------------------------------------------------
418 : ! Compute elements of tridiagonal matrix.
419 : !-----------------------------------------------------------------
420 :
421 : call get_matrix_elements_calc_Tsfc (nilyr, nslyr, &
422 : l_snow, l_cold, & ! LCOV_EXCL_LINE
423 : Tsf, Tbot, & ! LCOV_EXCL_LINE
424 : fsurfn, dfsurf_dT, & ! LCOV_EXCL_LINE
425 : Tin_init, Tsn_init, & ! LCOV_EXCL_LINE
426 : kh, Sswabs, & ! LCOV_EXCL_LINE
427 : Iswabs, & ! LCOV_EXCL_LINE
428 : etai, etas, & ! LCOV_EXCL_LINE
429 : sbdiag, diag, & ! LCOV_EXCL_LINE
430 74706776 : spdiag, rhs)
431 74706776 : if (icepack_warnings_aborted(subname)) return
432 :
433 : else
434 :
435 : call get_matrix_elements_know_Tsfc (nilyr, nslyr, &
436 : l_snow, Tbot, & ! LCOV_EXCL_LINE
437 : Tin_init, Tsn_init, & ! LCOV_EXCL_LINE
438 : kh, Sswabs, & ! LCOV_EXCL_LINE
439 : Iswabs, & ! LCOV_EXCL_LINE
440 : etai, etas, & ! LCOV_EXCL_LINE
441 : sbdiag, diag, & ! LCOV_EXCL_LINE
442 : spdiag, rhs, & ! LCOV_EXCL_LINE
443 59372394 : fcondtopn)
444 59372394 : if (icepack_warnings_aborted(subname)) return
445 :
446 : endif ! calc_Tsfc
447 :
448 : !-----------------------------------------------------------------
449 : ! Solve tridiagonal matrix to obtain the new temperatures.
450 : !-----------------------------------------------------------------
451 :
452 134079170 : nmat = nslyr + nilyr + 1 ! matrix dimension
453 :
454 0 : call tridiag_solver (nmat, sbdiag(:), &
455 : diag(:), spdiag(:), & ! LCOV_EXCL_LINE
456 134079170 : rhs(:), Tmat(:))
457 134079170 : if (icepack_warnings_aborted(subname)) return
458 :
459 : !-----------------------------------------------------------------
460 : ! Determine whether the computation has converged to an acceptable
461 : ! solution. Five conditions must be satisfied:
462 : !
463 : ! (1) Tsf <= 0 C.
464 : ! (2) Tsf is not oscillating; i.e., if both dTsf(niter) and
465 : ! dTsf(niter-1) have magnitudes greater than puny, then
466 : ! dTsf(niter)/dTsf(niter-1) cannot be a negative number
467 : ! with magnitude greater than 0.5.
468 : ! (3) abs(dTsf) < Tsf_errmax
469 : ! (4) If Tsf = 0 C, then the downward turbulent/radiative
470 : ! flux, fsurfn, must be greater than or equal to the downward
471 : ! conductive flux, fcondtopn.
472 : ! (5) The net energy added to the ice per unit time must equal
473 : ! the net change in internal ice energy per unit time,
474 : ! withinic the prescribed error ferrmax.
475 : !
476 : ! For briny ice (the standard case), zTsn and zTin are limited
477 : ! to prevent them from exceeding their melting temperatures.
478 : ! (Note that the specific heat formula for briny ice assumes
479 : ! that T < Tmlt.)
480 : ! For fresh ice there is no limiting, since there are cases
481 : ! when the only convergent solution has zTsn > 0 and/or zTin > 0.
482 : ! Above-zero temperatures are then reset to zero (with melting
483 : ! to conserve energy) in the thickness_changes subroutine.
484 : !-----------------------------------------------------------------
485 :
486 134079170 : if (calc_Tsfc) then
487 :
488 : !-----------------------------------------------------------------
489 : ! Reload Tsf from matrix solution
490 : !-----------------------------------------------------------------
491 :
492 74706776 : if (l_cold) then
493 54336220 : if (l_snow) then
494 51373951 : Tsf = Tmat(1)
495 : else
496 2962269 : Tsf = Tmat(1+nslyr)
497 : endif
498 : else ! melting surface
499 20370556 : Tsf = c0
500 : endif
501 :
502 : !-----------------------------------------------------------------
503 : ! Initialize convergence flag (true until proven false), dTsf,
504 : ! and temperature-averaging coefficients.
505 : ! Average only if test 1 or 2 fails.
506 : ! Initialize energy.
507 : !-----------------------------------------------------------------
508 :
509 74706776 : dTsf = Tsf - Tsf_start
510 74706776 : avg_Tsf = c0
511 :
512 : !-----------------------------------------------------------------
513 : ! Condition 1: check for Tsf > 0
514 : ! If Tsf > 0, set Tsf = 0, then average zTsn and zTin to force
515 : ! internal temps below their melting temps.
516 : !-----------------------------------------------------------------
517 :
518 74706776 : if (Tsf > puny) then
519 3687187 : Tsf = c0
520 3687187 : dTsf = -Tsf_start
521 3687187 : if (l_brine) avg_Tsi = c1 ! avg with starting temp
522 3687187 : converged = .false.
523 :
524 : !-----------------------------------------------------------------
525 : ! Condition 2: check for oscillating Tsf
526 : ! If oscillating, average all temps to increase rate of convergence.
527 : !-----------------------------------------------------------------
528 :
529 : elseif (niter > 1 & ! condition (2)
530 : .and. Tsf_start <= -puny & ! LCOV_EXCL_LINE
531 : .and. abs(dTsf) > puny & ! LCOV_EXCL_LINE
532 : .and. abs(dTsf_prev) > puny & ! LCOV_EXCL_LINE
533 71019589 : .and. -dTsf/(dTsf_prev+puny*puny) > p5) then
534 :
535 20129 : if (l_brine) then ! average with starting temp
536 20129 : avg_Tsf = c1
537 20129 : avg_Tsi = c1
538 : endif
539 20129 : dTsf = p5 * dTsf
540 20129 : converged = .false.
541 : endif
542 :
543 : !!! dTsf_prev = dTsf
544 :
545 : !-----------------------------------------------------------------
546 : ! If condition 2 failed, average new surface temperature with
547 : ! starting value.
548 : !-----------------------------------------------------------------
549 : Tsf = Tsf &
550 74706776 : + avg_Tsf * p5 * (Tsf_start - Tsf)
551 :
552 : endif ! calc_Tsfc
553 :
554 268158340 : do k = 1, nslyr
555 :
556 : !-----------------------------------------------------------------
557 : ! Reload zTsn from matrix solution
558 : !-----------------------------------------------------------------
559 :
560 134079170 : if (l_snow) then
561 122651473 : zTsn(k) = Tmat(k+1)
562 : else
563 11427697 : zTsn(k) = c0
564 : endif
565 134079170 : if (l_brine) zTsn(k) = min(zTsn(k), c0)
566 :
567 : !-----------------------------------------------------------------
568 : ! If condition 1 or 2 failed, average new snow layer
569 : ! temperatures with their starting values.
570 : !-----------------------------------------------------------------
571 11185659 : zTsn(k) = zTsn(k) &
572 134079170 : + avg_Tsi*p5*(Tsn_start(k)-zTsn(k))
573 :
574 : !-----------------------------------------------------------------
575 : ! Compute zqsn and increment new energy.
576 : !-----------------------------------------------------------------
577 134079170 : zqsn(k) = -rhos * (Lfresh - cp_ice*zTsn(k))
578 134079170 : enew = enew + hslyr * zqsn(k)
579 :
580 268158340 : Tsn_start(k) = zTsn(k) ! for next iteration
581 :
582 : enddo ! nslyr
583 :
584 1072633360 : dTmat(:) = c0
585 1072633360 : dqmat(:) = c0
586 1072633360 : reduce_kh(:) = .false.
587 1072633360 : do k = 1, nilyr
588 :
589 : !-----------------------------------------------------------------
590 : ! Reload zTin from matrix solution
591 : !-----------------------------------------------------------------
592 :
593 938554190 : zTin(k) = Tmat(k+1+nslyr)
594 :
595 938554190 : if (l_brine .and. zTin(k) > Tmlts(k) - puny) then
596 22475 : dTmat(k) = zTin(k) - Tmlts(k)
597 2876 : dqmat(k) = rhoi * dTmat(k) &
598 22475 : * (cp_ice - Lfresh * Tmlts(k)/zTin(k)**2)
599 : ! use this for the case that Tmlt changes by an amount dTmlt=Tmltnew-Tmlt(k)
600 : ! + rhoi * dTmlt &
601 : ! * (cp_ocn - cp_ice + Lfresh/zTin(k))
602 22475 : zTin(k) = Tmlts(k)
603 22475 : reduce_kh(k) = .true.
604 : endif
605 :
606 : !-----------------------------------------------------------------
607 : ! Condition 2b: check for oscillating zTin(1)
608 : ! If oscillating, average all ice temps to increase rate of convergence.
609 : !-----------------------------------------------------------------
610 :
611 938554190 : if (k==1 .and. .not.calc_Tsfc) then
612 59372394 : dTi1 = zTin(k) - Tin_start(k)
613 :
614 : if (niter > 1 & ! condition 2b
615 : .and. abs(dTi1) > puny & ! LCOV_EXCL_LINE
616 : .and. abs(dTi1_prev) > puny & ! LCOV_EXCL_LINE
617 59372394 : .and. -dTi1/(dTi1_prev+puny*puny) > p5) then
618 :
619 2883 : if (l_brine) avg_Tsi = c1
620 2883 : dTi1 = p5 * dTi1
621 2883 : converged = .false.
622 : endif
623 59372394 : dTi1_prev = dTi1
624 : endif ! k = 1 .and. calc_Tsfc = F
625 :
626 : !-----------------------------------------------------------------
627 : ! If condition 1 or 2 failed, average new ice layer
628 : ! temperatures with their starting values.
629 : !-----------------------------------------------------------------
630 78299613 : zTin(k) = zTin(k) &
631 938554190 : + avg_Tsi*p5*(Tin_start(k)-zTin(k))
632 :
633 : !-----------------------------------------------------------------
634 : ! Compute zqin and increment new energy.
635 : !-----------------------------------------------------------------
636 938554190 : if (l_brine) then
637 78299613 : zqin(k) = -rhoi * (cp_ice*(Tmlts(k)-zTin(k)) &
638 : + Lfresh*(c1-Tmlts(k)/zTin(k)) & ! LCOV_EXCL_LINE
639 938554190 : - cp_ocn*Tmlts(k))
640 : else
641 0 : zqin(k) = -rhoi * (-cp_ice*zTin(k) + Lfresh)
642 : endif
643 938554190 : enew = enew + hilyr * zqin(k)
644 938554190 : einex = einex + hilyr * dqmat(k)
645 :
646 1072633360 : Tin_start(k) = zTin(k) ! for next iteration
647 :
648 : enddo ! nilyr
649 :
650 134079170 : if (calc_Tsfc) then
651 :
652 : !-----------------------------------------------------------------
653 : ! Condition 3: check for large change in Tsf
654 : !-----------------------------------------------------------------
655 :
656 74706776 : if (abs(dTsf) > Tsf_errmax) then
657 28493554 : converged = .false.
658 : endif
659 :
660 : !-----------------------------------------------------------------
661 : ! Condition 4: check for fsurfn < fcondtopn with Tsf >= 0
662 : !-----------------------------------------------------------------
663 :
664 74706776 : fsurfn = fsurfn + dTsf*dfsurf_dT
665 74706776 : if (l_snow) then
666 64507405 : fcondtopn = kh(1) * (Tsf-zTsn(1))
667 : else
668 10199371 : fcondtopn = kh(1+nslyr) * (Tsf-zTin(1))
669 : endif
670 :
671 74706776 : if (Tsf >= c0 .and. fsurfn < fcondtopn) then
672 829138 : converged = .false.
673 : endif
674 :
675 74706776 : dTsf_prev = dTsf
676 :
677 : endif ! calc_Tsfc
678 :
679 : !-----------------------------------------------------------------
680 : ! Condition 5: check for energy conservation error
681 : ! Change in internal ice energy should equal net energy input.
682 : !-----------------------------------------------------------------
683 :
684 11185659 : fcondbot = kh(1+nslyr+nilyr) * &
685 145264829 : (zTin(nilyr) - Tbot)
686 :
687 : ! Flux extra energy out of the ice
688 134079170 : fcondbot = fcondbot + einex/dt
689 :
690 : ferr = abs( (enew-einit)/dt &
691 134079170 : - (fcondtopn - fcondbot + fswint) )
692 :
693 : ! factor of 0.9 allows for roundoff errors later
694 134079170 : if (ferr > 0.9_dbl_kind*ferrmax) then ! condition (5)
695 :
696 65015894 : converged = .false.
697 :
698 : ! reduce conductivity for next iteration
699 520127152 : do k = 1, nilyr
700 520127152 : if (reduce_kh(k) .and. dqmat(k) > c0) then
701 22475 : frac = max(0.5*(c1-ferr/abs(fcondtopn-fcondbot)),p1)
702 : ! frac = p1
703 22475 : kh(k+nslyr+1) = kh(k+nslyr+1) * frac
704 22475 : kh(k+nslyr) = kh(k+nslyr+1)
705 : endif
706 : enddo
707 :
708 : endif ! ferr
709 :
710 : endif ! convergence
711 :
712 : enddo ! temperature iteration niter
713 :
714 : !-----------------------------------------------------------------
715 : ! Check for convergence failures.
716 : !-----------------------------------------------------------------
717 59683339 : if (.not.converged) then
718 0 : write(warnstr,*) subname, 'Thermo iteration does not converge,'
719 0 : call icepack_warnings_add(warnstr)
720 0 : write(warnstr,*) subname, 'Ice thickness:', hilyr*nilyr
721 0 : call icepack_warnings_add(warnstr)
722 0 : write(warnstr,*) subname, 'Snow thickness:', hslyr*nslyr
723 0 : call icepack_warnings_add(warnstr)
724 0 : write(warnstr,*) subname, 'dTsf, Tsf_errmax:',dTsf_prev, &
725 0 : Tsf_errmax
726 0 : call icepack_warnings_add(warnstr)
727 0 : write(warnstr,*) subname, 'Tsf:', Tsf
728 0 : call icepack_warnings_add(warnstr)
729 0 : write(warnstr,*) subname, 'fsurf:', fsurfn
730 0 : call icepack_warnings_add(warnstr)
731 0 : write(warnstr,*) subname, 'fcondtop, fcondbot, fswint', &
732 0 : fcondtopn, fcondbot, fswint
733 0 : call icepack_warnings_add(warnstr)
734 0 : write(warnstr,*) subname, 'fswsfc', fswsfc
735 0 : call icepack_warnings_add(warnstr)
736 0 : write(warnstr,*) subname, 'Iswabs',(Iswabs(k),k=1,nilyr)
737 0 : call icepack_warnings_add(warnstr)
738 0 : write(warnstr,*) subname, 'Flux conservation error =', ferr
739 0 : call icepack_warnings_add(warnstr)
740 0 : write(warnstr,*) subname, 'Initial snow temperatures:'
741 0 : call icepack_warnings_add(warnstr)
742 0 : write(warnstr,*) subname, (Tsn_init(k),k=1,nslyr)
743 0 : call icepack_warnings_add(warnstr)
744 0 : write(warnstr,*) subname, 'Initial ice temperatures:'
745 0 : call icepack_warnings_add(warnstr)
746 0 : write(warnstr,*) subname, (Tin_init(k),k=1,nilyr)
747 0 : call icepack_warnings_add(warnstr)
748 0 : write(warnstr,*) subname, 'Matrix ice temperature diff:'
749 0 : call icepack_warnings_add(warnstr)
750 0 : write(warnstr,*) subname, (dTmat(k),k=1,nilyr)
751 0 : call icepack_warnings_add(warnstr)
752 0 : write(warnstr,*) subname, 'dqmat*hilyr/dt:'
753 0 : call icepack_warnings_add(warnstr)
754 0 : write(warnstr,*) subname, (hilyr*dqmat(k)/dt,k=1,nilyr)
755 0 : call icepack_warnings_add(warnstr)
756 0 : write(warnstr,*) subname, 'Final snow temperatures:'
757 0 : call icepack_warnings_add(warnstr)
758 0 : write(warnstr,*) subname, (zTsn(k),k=1,nslyr)
759 0 : call icepack_warnings_add(warnstr)
760 0 : write(warnstr,*) subname, 'Matrix ice temperature diff:'
761 0 : call icepack_warnings_add(warnstr)
762 0 : write(warnstr,*) subname, (dTmat(k),k=1,nilyr)
763 0 : call icepack_warnings_add(warnstr)
764 0 : write(warnstr,*) subname, 'dqmat*hilyr/dt:'
765 0 : call icepack_warnings_add(warnstr)
766 0 : write(warnstr,*) subname, (hilyr*dqmat(k)/dt,k=1,nilyr)
767 0 : call icepack_warnings_add(warnstr)
768 0 : write(warnstr,*) subname, 'Final ice temperatures:'
769 0 : call icepack_warnings_add(warnstr)
770 0 : write(warnstr,*) subname, (zTin(k),k=1,nilyr)
771 0 : call icepack_warnings_add(warnstr)
772 0 : write(warnstr,*) subname, 'Ice melting temperatures:'
773 0 : call icepack_warnings_add(warnstr)
774 0 : write(warnstr,*) subname, (Tmlts(k),k=1,nilyr)
775 0 : call icepack_warnings_add(warnstr)
776 0 : write(warnstr,*) subname, 'Ice bottom temperature:', Tbot
777 0 : call icepack_warnings_add(warnstr)
778 0 : write(warnstr,*) subname, 'dT initial:'
779 0 : call icepack_warnings_add(warnstr)
780 0 : write(warnstr,*) subname, (Tmlts(k)-Tin_init(k),k=1,nilyr)
781 0 : call icepack_warnings_add(warnstr)
782 0 : write(warnstr,*) subname, 'dT final:'
783 0 : call icepack_warnings_add(warnstr)
784 0 : write(warnstr,*) subname, (Tmlts(k)-zTin(k),k=1,nilyr)
785 0 : call icepack_warnings_add(warnstr)
786 0 : write(warnstr,*) subname, 'zSin'
787 0 : call icepack_warnings_add(warnstr)
788 0 : write(warnstr,*) subname, (zSin(k),k=1,nilyr)
789 0 : call icepack_warnings_add(warnstr)
790 0 : call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
791 0 : call icepack_warnings_add(subname//" temperature_changes: Thermo iteration does not converge" )
792 0 : return
793 : endif
794 :
795 59683339 : if (calc_Tsfc) then
796 :
797 : ! update fluxes that depend on Tsf
798 27633970 : flwoutn = flwoutn + dTsf_prev * dflwout_dT
799 27633970 : fsensn = fsensn + dTsf_prev * dfsens_dT
800 27633970 : flatn = flatn + dTsf_prev * dflat_dT
801 :
802 : endif ! calc_Tsfc
803 :
804 : end subroutine temperature_changes
805 :
806 : !=======================================================================
807 : !
808 : ! Compute thermal conductivity at interfaces (held fixed during
809 : ! the subsequent iteration).
810 : !
811 : ! NOTE: Ice conductivity must be >= kimin
812 : !
813 : ! authors William H. Lipscomb, LANL
814 : ! C. M. Bitz, UW
815 :
816 59683339 : subroutine conductivity (l_snow, &
817 : nilyr, nslyr, & ! LCOV_EXCL_LINE
818 : hilyr, hslyr, & ! LCOV_EXCL_LINE
819 59683339 : zTin, kh, zSin)
820 :
821 : logical (kind=log_kind), intent(in) :: &
822 : l_snow ! true if snow temperatures are computed
823 :
824 : integer (kind=int_kind), intent(in) :: &
825 : nilyr , & ! number of ice layers ! LCOV_EXCL_LINE
826 : nslyr ! number of snow layers
827 :
828 : real (kind=dbl_kind), intent(in) :: &
829 : hilyr , & ! ice layer thickness (same for all ice layers) ! LCOV_EXCL_LINE
830 : hslyr ! snow layer thickness (same for all snow layers)
831 :
832 : real (kind=dbl_kind), dimension (:), intent(in) :: &
833 : zTin , & ! internal ice layer temperatures ! LCOV_EXCL_LINE
834 : zSin ! internal ice layer salinities
835 :
836 : real (kind=dbl_kind), dimension (nilyr+nslyr+1), &
837 : intent(out) :: & ! LCOV_EXCL_LINE
838 : kh ! effective conductivity at interfaces (W m-2 deg-1)
839 :
840 : ! local variables
841 :
842 : integer (kind=int_kind) :: &
843 : k ! vertical index
844 :
845 : real (kind=dbl_kind), dimension (nilyr) :: &
846 151931450 : kilyr ! thermal cond at ice layer midpoints (W m-1 deg-1)
847 :
848 : real (kind=dbl_kind), dimension (nslyr) :: &
849 119366678 : kslyr ! thermal cond at snow layer midpoints (W m-1 deg-1)
850 :
851 : character(len=*),parameter :: subname='(conductivity)'
852 :
853 : ! interior snow layers (simple for now, but may be fancier later)
854 119366678 : do k = 1, nslyr
855 119366678 : kslyr(k) = ksno
856 : enddo ! nslyr
857 :
858 : ! interior ice layers
859 59683339 : if (conduct == 'MU71') then
860 : ! Maykut and Untersteiner 1971 form (with Wettlaufer 1991 constants)
861 186577328 : do k = 1, nilyr
862 163255162 : kilyr(k) = kice + betak*zSin(k)/min(-puny,zTin(k))
863 186577328 : kilyr(k) = max (kilyr(k), kimin)
864 : enddo ! nilyr
865 : else
866 : ! Pringle et al JGR 2007 'bubbly brine'
867 290889384 : do k = 1, nilyr
868 36713005 : kilyr(k) = (2.11_dbl_kind - 0.011_dbl_kind*zTin(k) &
869 : + 0.09_dbl_kind*zSin(k)/min(-puny,zTin(k))) & ! LCOV_EXCL_LINE
870 254528211 : * rhoi / 917._dbl_kind
871 290889384 : kilyr(k) = max (kilyr(k), kimin)
872 : enddo ! nilyr
873 : endif ! conductivity
874 :
875 : ! top snow interface, top and bottom ice interfaces
876 : ! top of snow layer; top surface of top ice layer
877 59683339 : if (l_snow) then
878 56027048 : kh(1) = c2 * kslyr(1) / hslyr
879 10532730 : kh(1+nslyr) = c2 * kslyr(nslyr) * kilyr(1) / &
880 : ( kslyr(nslyr)*hilyr + & ! LCOV_EXCL_LINE
881 71826143 : kilyr(1 )*hslyr )
882 : else
883 3656291 : kh(1) = c0
884 3656291 : kh(1+nslyr) = c2 * kilyr(1) / hilyr
885 : endif
886 :
887 : ! bottom surface of bottom ice layer
888 59683339 : kh(1+nslyr+nilyr) = c2 * kilyr(nilyr) / hilyr
889 :
890 : ! interior snow interfaces
891 :
892 59683339 : if (nslyr > 1) then
893 0 : do k = 2, nslyr
894 0 : if (l_snow) then
895 0 : kh(k) = c2 * kslyr(k-1) * kslyr(k) / &
896 0 : ((kslyr(k-1) + kslyr(k))*hslyr)
897 : else
898 0 : kh(k) = c0
899 : endif
900 : enddo ! nilyr
901 : endif ! nslyr > 1
902 :
903 : ! interior ice interfaces
904 417783373 : do k = 2, nilyr
905 65129544 : kh(k+nslyr) = c2 * kilyr(k-1) * kilyr(k) / &
906 482912917 : ((kilyr(k-1) + kilyr(k))*hilyr)
907 : enddo ! nilyr
908 :
909 59683339 : end subroutine conductivity
910 :
911 : !=======================================================================
912 : !
913 : ! Compute radiative and turbulent fluxes and their derivatives
914 : ! with respect to Tsf.
915 : !
916 : ! authors William H. Lipscomb, LANL
917 : ! C. M. Bitz, UW
918 :
919 345756478 : subroutine surface_fluxes (Tsf, fswsfc, &
920 : rhoa, flw, & ! LCOV_EXCL_LINE
921 : potT, Qa, & ! LCOV_EXCL_LINE
922 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
923 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
924 : flatn, fsurfn, & ! LCOV_EXCL_LINE
925 : dflwout_dT, dfsens_dT, & ! LCOV_EXCL_LINE
926 : dflat_dT, dfsurf_dT)
927 :
928 : real (kind=dbl_kind), intent(in) :: &
929 : Tsf ! ice/snow surface temperature, Tsfcn
930 :
931 : real (kind=dbl_kind), intent(in) :: &
932 : fswsfc , & ! SW absorbed at ice/snow surface (W m-2) ! LCOV_EXCL_LINE
933 : rhoa , & ! air density (kg/m^3) ! LCOV_EXCL_LINE
934 : flw , & ! incoming longwave radiation (W/m^2) ! LCOV_EXCL_LINE
935 : potT , & ! air potential temperature (K) ! LCOV_EXCL_LINE
936 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
937 : shcoef , & ! transfer coefficient for sensible heat ! LCOV_EXCL_LINE
938 : lhcoef ! transfer coefficient for latent heat
939 :
940 : real (kind=dbl_kind), &
941 : intent(inout) :: & ! LCOV_EXCL_LINE
942 : fsensn , & ! surface downward sensible heat (W m-2) ! LCOV_EXCL_LINE
943 : flatn , & ! surface downward latent heat (W m-2) ! LCOV_EXCL_LINE
944 : flwoutn , & ! upward LW at surface (W m-2) ! LCOV_EXCL_LINE
945 : fsurfn ! net flux to top surface, excluding fcondtopn
946 :
947 : real (kind=dbl_kind), &
948 : intent(inout) :: & ! LCOV_EXCL_LINE
949 : dfsens_dT , & ! deriv of fsens wrt Tsf (W m-2 deg-1) ! LCOV_EXCL_LINE
950 : dflat_dT , & ! deriv of flat wrt Tsf (W m-2 deg-1) ! LCOV_EXCL_LINE
951 : dflwout_dT ! deriv of flwout wrt Tsf (W m-2 deg-1)
952 :
953 : real (kind=dbl_kind), &
954 : intent(inout) :: & ! LCOV_EXCL_LINE
955 : dfsurf_dT ! derivative of fsurfn wrt Tsf
956 :
957 : character(len=*),parameter :: subname='(surface_fluxes)'
958 :
959 : ! surface heat flux
960 : call surface_heat_flux(Tsf, fswsfc, &
961 : rhoa, flw, & ! LCOV_EXCL_LINE
962 : potT, Qa, & ! LCOV_EXCL_LINE
963 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
964 : flwoutn, fsensn, & ! LCOV_EXCL_LINE
965 345756478 : flatn, fsurfn)
966 345756478 : if (icepack_warnings_aborted(subname)) return
967 :
968 : ! derivative of heat flux with respect to surface temperature
969 : call dsurface_heat_flux_dTsf(Tsf, rhoa, &
970 : shcoef, lhcoef, & ! LCOV_EXCL_LINE
971 : dfsurf_dT, dflwout_dT, & ! LCOV_EXCL_LINE
972 345756478 : dfsens_dT, dflat_dT)
973 345756478 : if (icepack_warnings_aborted(subname)) return
974 :
975 : end subroutine surface_fluxes
976 :
977 : !=======================================================================
978 : !
979 : ! Compute terms in tridiagonal matrix that will be solved to find
980 : ! the new vertical temperature profile
981 : ! This routine is for the case in which Tsfc is being computed.
982 : !
983 : ! authors William H. Lipscomb, LANL
984 : ! C. M. Bitz, UW
985 : !
986 : ! March 2004 by William H. Lipscomb for multiple snow layers
987 : ! April 2008 by E. C. Hunke, divided into two routines based on calc_Tsfc
988 :
989 74706776 : subroutine get_matrix_elements_calc_Tsfc (nilyr, nslyr, &
990 : l_snow, l_cold, & ! LCOV_EXCL_LINE
991 : Tsf, Tbot, & ! LCOV_EXCL_LINE
992 : fsurfn, dfsurf_dT, & ! LCOV_EXCL_LINE
993 : Tin_init, Tsn_init, & ! LCOV_EXCL_LINE
994 : kh, Sswabs, & ! LCOV_EXCL_LINE
995 : Iswabs, & ! LCOV_EXCL_LINE
996 : etai, etas, & ! LCOV_EXCL_LINE
997 : sbdiag, diag, & ! LCOV_EXCL_LINE
998 74706776 : spdiag, rhs)
999 :
1000 : integer (kind=int_kind), intent(in) :: &
1001 : nilyr , & ! number of ice layers ! LCOV_EXCL_LINE
1002 : nslyr ! number of snow layers
1003 :
1004 : logical (kind=log_kind), &
1005 : intent(in) :: & ! LCOV_EXCL_LINE
1006 : l_snow , & ! true if snow temperatures are computed ! LCOV_EXCL_LINE
1007 : l_cold ! true if surface temperature is computed
1008 :
1009 : real (kind=dbl_kind), intent(in) :: &
1010 : Tsf ! ice/snow top surface temp (deg C)
1011 :
1012 : real (kind=dbl_kind), intent(in) :: &
1013 : fsurfn , & ! net flux to top surface, excluding fcondtopn (W/m^2) ! LCOV_EXCL_LINE
1014 : Tbot ! ice bottom surface temperature (deg C)
1015 :
1016 : real (kind=dbl_kind), intent(in) :: &
1017 : dfsurf_dT ! derivative of fsurf wrt Tsf
1018 :
1019 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1020 : etai , & ! dt / (rho*cp*h) for ice layers ! LCOV_EXCL_LINE
1021 : Tin_init , & ! ice temp at beginning of time step ! LCOV_EXCL_LINE
1022 : Sswabs , & ! SW radiation absorbed in snow layers (W m-2) ! LCOV_EXCL_LINE
1023 : Iswabs , & ! absorbed SW flux in ice layers ! LCOV_EXCL_LINE
1024 : etas , & ! dt / (rho*cp*h) for snow layers ! LCOV_EXCL_LINE
1025 : Tsn_init ! snow temp at beginning of time step
1026 : ! Note: no absorbed SW in snow layers
1027 :
1028 : real (kind=dbl_kind), dimension (nslyr+nilyr+1), &
1029 : intent(in) :: & ! LCOV_EXCL_LINE
1030 : kh ! effective conductivity at layer interfaces
1031 :
1032 : real (kind=dbl_kind), dimension (nslyr+nilyr+1), &
1033 : intent(inout) :: & ! LCOV_EXCL_LINE
1034 : sbdiag , & ! sub-diagonal matrix elements ! LCOV_EXCL_LINE
1035 : diag , & ! diagonal matrix elements ! LCOV_EXCL_LINE
1036 : spdiag , & ! super-diagonal matrix elements ! LCOV_EXCL_LINE
1037 : rhs ! rhs of tri-diagonal matrix eqn.
1038 :
1039 : ! local variables
1040 :
1041 : integer (kind=int_kind) :: &
1042 : k, ki, kr ! vertical indices and row counters
1043 :
1044 : character(len=*),parameter :: subname='(get_matrix_elements_calc_Tsrf)'
1045 :
1046 : !-----------------------------------------------------------------
1047 : ! Initialize matrix elements.
1048 : ! Note: When we do not need to solve for the surface or snow
1049 : ! temperature, we solve dummy equations with solution T = 0.
1050 : ! Ice layers are fully initialized below.
1051 : !-----------------------------------------------------------------
1052 :
1053 224120328 : do k = 1, nslyr+1
1054 149413552 : sbdiag(k) = c0
1055 149413552 : diag (k) = c1
1056 149413552 : spdiag(k) = c0
1057 224120328 : rhs (k) = c0
1058 : enddo
1059 :
1060 : !-----------------------------------------------------------------
1061 : ! Compute matrix elements
1062 : !
1063 : ! Four possible cases to solve:
1064 : ! (1) Cold surface (Tsf < 0), snow present
1065 : ! (2) Melting surface (Tsf = 0), snow present
1066 : ! (3) Cold surface (Tsf < 0), no snow
1067 : ! (4) Melting surface (Tsf = 0), no snow
1068 : !-----------------------------------------------------------------
1069 :
1070 : !-----------------------------------------------------------------
1071 : ! Tsf equation for case of cold surface (with or without snow)
1072 : !-----------------------------------------------------------------
1073 :
1074 74706776 : if (l_cold) then
1075 54336220 : if (l_snow) then
1076 51373951 : k = 1
1077 : else ! no snow
1078 2962269 : k = 1 + nslyr
1079 : endif
1080 54336220 : kr = k
1081 :
1082 54336220 : sbdiag(kr) = c0
1083 54336220 : diag (kr) = dfsurf_dT - kh(k)
1084 54336220 : spdiag(kr) = kh(k)
1085 54336220 : rhs (kr) = dfsurf_dT*Tsf - fsurfn
1086 : endif ! l_cold
1087 :
1088 : !-----------------------------------------------------------------
1089 : ! top snow layer
1090 : !-----------------------------------------------------------------
1091 : ! k = 1
1092 : ! kr = 2
1093 :
1094 74706776 : if (l_snow) then
1095 64507405 : if (l_cold) then
1096 51373951 : sbdiag(2) = -etas(1) * kh(1)
1097 51373951 : spdiag(2) = -etas(1) * kh(2)
1098 2507186 : diag (2) = c1 &
1099 51373951 : + etas(1) * (kh(1) + kh(2))
1100 2507186 : rhs (2) = Tsn_init(1) &
1101 51373951 : + etas(1) * Sswabs(1)
1102 : else ! melting surface
1103 13133454 : sbdiag(2) = c0
1104 13133454 : spdiag(2) = -etas(1) * kh(2)
1105 572686 : diag (2) = c1 &
1106 13133454 : + etas(1) * (kh(1) + kh(2))
1107 572686 : rhs (2) = Tsn_init(1) &
1108 : + etas(1)*kh(1)*Tsf & ! LCOV_EXCL_LINE
1109 13133454 : + etas(1) * Sswabs(1)
1110 : endif ! l_cold
1111 : endif ! l_snow
1112 :
1113 : !-----------------------------------------------------------------
1114 : ! remaining snow layers
1115 : !-----------------------------------------------------------------
1116 :
1117 74706776 : if (nslyr > 1) then
1118 :
1119 0 : do k = 2, nslyr
1120 0 : kr = k + 1
1121 :
1122 0 : if (l_snow) then
1123 0 : sbdiag(kr) = -etas(k) * kh(k)
1124 0 : spdiag(kr) = -etas(k) * kh(k+1)
1125 0 : diag (kr) = c1 &
1126 0 : + etas(k) * (kh(k) + kh(k+1))
1127 0 : rhs (kr) = Tsn_init(k) &
1128 0 : + etas(k) * Sswabs(k)
1129 : endif
1130 : enddo ! nslyr
1131 :
1132 : endif ! nslyr > 1
1133 :
1134 74706776 : if (nilyr > 1) then
1135 :
1136 : !-----------------------------------------------------------------
1137 : ! top ice layer
1138 : !-----------------------------------------------------------------
1139 :
1140 74706776 : ki = 1
1141 74706776 : k = ki + nslyr
1142 74706776 : kr = k + 1
1143 :
1144 74706776 : if (l_snow .or. l_cold) then
1145 67469674 : sbdiag(kr) = -etai(ki) * kh(k)
1146 67469674 : spdiag(kr) = -etai(ki) * kh(k+1)
1147 3270253 : diag (kr) = c1 &
1148 67469674 : + etai(ki) * (kh(k) + kh(k+1))
1149 3270253 : rhs (kr) = Tin_init(ki) &
1150 67469674 : + etai(ki)*Iswabs(ki)
1151 : else ! no snow, warm surface
1152 7237102 : sbdiag(kr) = c0
1153 7237102 : spdiag(kr) = -etai(ki) * kh(k+1)
1154 80260 : diag (kr) = c1 &
1155 7237102 : + etai(ki) * (kh(k) + kh(k+1))
1156 80260 : rhs (kr) = Tin_init(ki) &
1157 : + etai(ki)*Iswabs(ki) & ! LCOV_EXCL_LINE
1158 7237102 : + etai(ki)*kh(k)*Tsf
1159 : endif
1160 :
1161 : !-----------------------------------------------------------------
1162 : ! bottom ice layer
1163 : !-----------------------------------------------------------------
1164 :
1165 74706776 : ki = nilyr
1166 74706776 : k = ki + nslyr
1167 74706776 : kr = k + 1
1168 :
1169 74706776 : sbdiag(kr) = -etai(ki) * kh(k)
1170 74706776 : spdiag(kr) = c0
1171 3350513 : diag (kr) = c1 &
1172 74706776 : + etai(ki) * (kh(k) + kh(k+1))
1173 3350513 : rhs (kr) = Tin_init(ki) &
1174 : + etai(ki)*Iswabs(ki) & ! LCOV_EXCL_LINE
1175 74706776 : + etai(ki)*kh(k+1)*Tbot
1176 :
1177 : else ! nilyr = 1
1178 :
1179 : !-----------------------------------------------------------------
1180 : ! single ice layer
1181 : !-----------------------------------------------------------------
1182 :
1183 0 : ki = 1
1184 0 : k = ki + nslyr
1185 0 : kr = k + 1
1186 :
1187 0 : if (l_snow .or. l_cold) then
1188 0 : sbdiag(kr) = -etai(ki) * kh(k)
1189 0 : spdiag(kr) = c0
1190 0 : diag (kr) = c1 &
1191 0 : + etai(ki) * (kh(k) + kh(k+1))
1192 0 : rhs (kr) = Tin_init(ki) &
1193 : + etai(ki) * Iswabs(ki) & ! LCOV_EXCL_LINE
1194 0 : + etai(ki) * kh(k+1)*Tbot
1195 : else ! no snow, warm surface
1196 0 : sbdiag(kr) = c0
1197 0 : spdiag(kr) = c0
1198 0 : diag (kr) = c1 &
1199 0 : + etai(ki) * (kh(k) + kh(k+1))
1200 0 : rhs (kr) = Tin_init(ki) &
1201 : + etai(ki) * Iswabs(ki) & ! LCOV_EXCL_LINE
1202 : + etai(ki) * kh(k)*Tsf & ! LCOV_EXCL_LINE
1203 0 : + etai(ki) * kh(k+1)*Tbot
1204 : endif
1205 :
1206 : endif ! nilyr > 1
1207 :
1208 : !-----------------------------------------------------------------
1209 : ! interior ice layers
1210 : !-----------------------------------------------------------------
1211 :
1212 448240656 : do ki = 2, nilyr-1
1213 :
1214 373533880 : k = ki + nslyr
1215 373533880 : kr = k + 1
1216 :
1217 373533880 : sbdiag(kr) = -etai(ki) * kh(k)
1218 373533880 : spdiag(kr) = -etai(ki) * kh(k+1)
1219 16752565 : diag (kr) = c1 &
1220 373533880 : + etai(ki) * (kh(k) + kh(k+1))
1221 16752565 : rhs (kr) = Tin_init(ki) &
1222 448240656 : + etai(ki)*Iswabs(ki)
1223 : enddo ! nilyr
1224 :
1225 74706776 : end subroutine get_matrix_elements_calc_Tsfc
1226 :
1227 : !=======================================================================
1228 : !
1229 : ! Compute terms in tridiagonal matrix that will be solved to find
1230 : ! the new vertical temperature profile
1231 : ! This routine is for the case in which Tsfc is already known.
1232 : !
1233 : ! authors William H. Lipscomb, LANL
1234 : ! C. M. Bitz, UW
1235 : !
1236 : ! March 2004 by William H. Lipscomb for multiple snow layers
1237 : ! April 2008 by E. C. Hunke, divided into two routines based on calc_Tsfc
1238 :
1239 59372394 : subroutine get_matrix_elements_know_Tsfc (nilyr, nslyr, &
1240 : l_snow, Tbot, & ! LCOV_EXCL_LINE
1241 : Tin_init, Tsn_init, & ! LCOV_EXCL_LINE
1242 : kh, Sswabs, & ! LCOV_EXCL_LINE
1243 : Iswabs, & ! LCOV_EXCL_LINE
1244 : etai, etas, & ! LCOV_EXCL_LINE
1245 : sbdiag, diag, & ! LCOV_EXCL_LINE
1246 : spdiag, rhs, & ! LCOV_EXCL_LINE
1247 : fcondtopn)
1248 :
1249 : integer (kind=int_kind), intent(in) :: &
1250 : nilyr , & ! number of ice layers ! LCOV_EXCL_LINE
1251 : nslyr ! number of snow layers
1252 :
1253 : logical (kind=log_kind), &
1254 : intent(in) :: & ! LCOV_EXCL_LINE
1255 : l_snow ! true if snow temperatures are computed
1256 :
1257 : real (kind=dbl_kind), intent(in) :: &
1258 : Tbot ! ice bottom surface temperature (deg C)
1259 :
1260 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1261 : etai , & ! dt / (rho*cp*h) for ice layers ! LCOV_EXCL_LINE
1262 : Tin_init , & ! ice temp at beginning of time step ! LCOV_EXCL_LINE
1263 : Sswabs , & ! SW radiation absorbed in snow layers (W m-2) ! LCOV_EXCL_LINE
1264 : Iswabs , & ! absorbed SW flux in ice layers ! LCOV_EXCL_LINE
1265 : etas , & ! dt / (rho*cp*h) for snow layers ! LCOV_EXCL_LINE
1266 : Tsn_init ! snow temp at beginning of time step
1267 : ! Note: no absorbed SW in snow layers
1268 :
1269 : real (kind=dbl_kind), dimension (nslyr+nilyr+1), &
1270 : intent(in) :: & ! LCOV_EXCL_LINE
1271 : kh ! effective conductivity at layer interfaces
1272 :
1273 : real (kind=dbl_kind), dimension (nslyr+nilyr+1), &
1274 : intent(inout) :: & ! LCOV_EXCL_LINE
1275 : sbdiag , & ! sub-diagonal matrix elements ! LCOV_EXCL_LINE
1276 : diag , & ! diagonal matrix elements ! LCOV_EXCL_LINE
1277 : spdiag , & ! super-diagonal matrix elements ! LCOV_EXCL_LINE
1278 : rhs ! rhs of tri-diagonal matrix eqn.
1279 :
1280 : real (kind=dbl_kind), intent(in), &
1281 : optional :: & ! LCOV_EXCL_LINE
1282 : fcondtopn ! conductive flux at top sfc, positive down (W/m^2)
1283 :
1284 : ! local variables
1285 :
1286 : integer (kind=int_kind) :: &
1287 : k, ki, kr ! vertical indices and row counters
1288 :
1289 : character(len=*),parameter :: subname='(get_matrix_elements_know_Tsrf)'
1290 :
1291 : !-----------------------------------------------------------------
1292 : ! Initialize matrix elements.
1293 : ! Note: When we do not need to solve for the surface or snow
1294 : ! temperature, we solve dummy equations with solution T = 0.
1295 : ! Ice layers are fully initialized below.
1296 : !-----------------------------------------------------------------
1297 :
1298 178117182 : do k = 1, nslyr+1
1299 118744788 : sbdiag(k) = c0
1300 118744788 : diag (k) = c1
1301 118744788 : spdiag(k) = c0
1302 178117182 : rhs (k) = c0
1303 : enddo
1304 :
1305 : !-----------------------------------------------------------------
1306 : ! Compute matrix elements
1307 : !
1308 : ! Four possible cases to solve:
1309 : ! (1) Cold surface (Tsf < 0), snow present
1310 : ! (2) Melting surface (Tsf = 0), snow present
1311 : ! (3) Cold surface (Tsf < 0), no snow
1312 : ! (4) Melting surface (Tsf = 0), no snow
1313 : !-----------------------------------------------------------------
1314 :
1315 : !-----------------------------------------------------------------
1316 : ! top snow layer
1317 : !-----------------------------------------------------------------
1318 : ! k = 1
1319 : ! kr = 2
1320 :
1321 59372394 : if (l_snow) then
1322 58144068 : sbdiag(2) = c0
1323 58144068 : spdiag(2) = -etas(1) * kh(2)
1324 7658066 : diag (2) = c1 &
1325 58144068 : + etas(1) * kh(2)
1326 7658066 : rhs (2) = Tsn_init(1) &
1327 : + etas(1) * Sswabs(1) & ! LCOV_EXCL_LINE
1328 58144068 : + etas(1) * fcondtopn
1329 : endif ! l_snow
1330 :
1331 : !-----------------------------------------------------------------
1332 : ! remaining snow layers
1333 : !-----------------------------------------------------------------
1334 :
1335 59372394 : if (nslyr > 1) then
1336 :
1337 0 : do k = 2, nslyr
1338 0 : kr = k + 1
1339 :
1340 0 : if (l_snow) then
1341 0 : sbdiag(kr) = -etas(k) * kh(k)
1342 0 : spdiag(kr) = -etas(k) * kh(k+1)
1343 0 : diag (kr) = c1 &
1344 0 : + etas(k) * (kh(k) + kh(k+1))
1345 0 : rhs (kr) = Tsn_init(k) &
1346 0 : + etas(k) * Sswabs(k)
1347 : endif
1348 :
1349 : enddo ! nslyr
1350 :
1351 : endif ! nslyr > 1
1352 :
1353 59372394 : if (nilyr > 1) then
1354 :
1355 : !-----------------------------------------------------------------
1356 : ! top ice layer
1357 : !-----------------------------------------------------------------
1358 :
1359 59372394 : ki = 1
1360 59372394 : k = ki + nslyr
1361 59372394 : kr = k + 1
1362 :
1363 59372394 : if (l_snow) then
1364 :
1365 58144068 : sbdiag(kr) = -etai(ki) * kh(k)
1366 58144068 : spdiag(kr) = -etai(ki) * kh(k+1)
1367 7658066 : diag (kr) = c1 &
1368 58144068 : + etai(ki) * (kh(k) + kh(k+1))
1369 7658066 : rhs (kr) = Tin_init(ki) &
1370 58144068 : + etai(ki) * Iswabs(ki)
1371 : else
1372 1228326 : sbdiag(kr) = c0
1373 1228326 : spdiag(kr) = -etai(ki) * kh(k+1)
1374 177080 : diag (kr) = c1 &
1375 1228326 : + etai(ki) * kh(k+1)
1376 177080 : rhs (kr) = Tin_init(ki) &
1377 : + etai(ki) * Iswabs(ki) & ! LCOV_EXCL_LINE
1378 1228326 : + etai(ki) * fcondtopn
1379 : endif ! l_snow
1380 :
1381 : !-----------------------------------------------------------------
1382 : ! bottom ice layer
1383 : !-----------------------------------------------------------------
1384 :
1385 59372394 : ki = nilyr
1386 59372394 : k = ki + nslyr
1387 59372394 : kr = k + 1
1388 :
1389 59372394 : sbdiag(kr) = -etai(ki) * kh(k)
1390 59372394 : spdiag(kr) = c0
1391 7835146 : diag (kr) = c1 &
1392 59372394 : + etai(ki) * (kh(k) + kh(k+1))
1393 7835146 : rhs (kr) = Tin_init(ki) &
1394 : + etai(ki)*Iswabs(ki) & ! LCOV_EXCL_LINE
1395 59372394 : + etai(ki)*kh(k+1)*Tbot
1396 :
1397 : else ! nilyr = 1
1398 :
1399 : !-----------------------------------------------------------------
1400 : ! single ice layer
1401 : !-----------------------------------------------------------------
1402 :
1403 0 : ki = 1
1404 0 : k = ki + nslyr
1405 0 : kr = k + 1
1406 :
1407 0 : if (l_snow) then
1408 0 : sbdiag(kr) = -etai(ki) * kh(k)
1409 0 : spdiag(kr) = c0
1410 0 : diag (kr) = c1 &
1411 0 : + etai(ki) * (kh(k) + kh(k+1))
1412 0 : rhs (kr) = Tin_init(ki) &
1413 : + etai(ki) * Iswabs(ki) & ! LCOV_EXCL_LINE
1414 0 : + etai(ki) * kh(k+1)*Tbot
1415 : else
1416 0 : sbdiag(kr) = c0
1417 0 : spdiag(kr) = c0
1418 0 : diag (kr) = c1 &
1419 0 : + etai(ki) * kh(k+1)
1420 0 : rhs (kr) = Tin_init(ki) &
1421 : + etai(ki) * Iswabs(ki) & ! LCOV_EXCL_LINE
1422 : + etai(ki) * fcondtopn & ! LCOV_EXCL_LINE
1423 0 : + etai(ki) * kh(k+1)*Tbot
1424 : endif
1425 :
1426 : endif ! nilyr > 1
1427 :
1428 : !-----------------------------------------------------------------
1429 : ! interior ice layers
1430 : !-----------------------------------------------------------------
1431 :
1432 356234364 : do ki = 2, nilyr-1
1433 :
1434 296861970 : k = ki + nslyr
1435 296861970 : kr = k + 1
1436 :
1437 296861970 : sbdiag(kr) = -etai(ki) * kh(k)
1438 296861970 : spdiag(kr) = -etai(ki) * kh(k+1)
1439 39175730 : diag (kr) = c1 &
1440 296861970 : + etai(ki) * (kh(k) + kh(k+1))
1441 39175730 : rhs (kr) = Tin_init(ki) &
1442 356234364 : + etai(ki)*Iswabs(ki)
1443 :
1444 : enddo ! nilyr
1445 :
1446 59372394 : end subroutine get_matrix_elements_know_Tsfc
1447 :
1448 : !=======================================================================
1449 : !
1450 : ! Tridiagonal matrix solver--used to solve the implicit vertical heat
1451 : ! equation in ice and snow
1452 : !
1453 : ! authors William H. Lipscomb, LANL
1454 : ! C. M. Bitz, UW
1455 :
1456 134079170 : subroutine tridiag_solver (nmat, sbdiag, &
1457 : diag, spdiag, & ! LCOV_EXCL_LINE
1458 134079170 : rhs, xout)
1459 :
1460 : integer (kind=int_kind), intent(in) :: &
1461 : nmat ! matrix dimension
1462 :
1463 : real (kind=dbl_kind), dimension (:), intent(in) :: &
1464 : sbdiag , & ! sub-diagonal matrix elements ! LCOV_EXCL_LINE
1465 : diag , & ! diagonal matrix elements ! LCOV_EXCL_LINE
1466 : spdiag , & ! super-diagonal matrix elements ! LCOV_EXCL_LINE
1467 : rhs ! rhs of tri-diagonal matrix eqn.
1468 :
1469 : real (kind=dbl_kind), dimension (:), intent(inout) :: &
1470 : xout ! solution vector
1471 :
1472 : ! local variables
1473 :
1474 : integer (kind=int_kind) :: &
1475 : k ! row counter
1476 :
1477 : real (kind=dbl_kind) :: &
1478 11185659 : wbeta ! temporary matrix variable
1479 :
1480 : real (kind=dbl_kind), dimension(nmat) :: &
1481 357643612 : wgamma ! temporary matrix variable
1482 :
1483 : character(len=*),parameter :: subname='(tridiag_solver)'
1484 :
1485 134079170 : wbeta = diag(1)
1486 134079170 : xout(1) = rhs(1) / wbeta
1487 :
1488 1206712530 : do k = 2, nmat
1489 1072633360 : wgamma(k) = spdiag(k-1) / wbeta
1490 1072633360 : wbeta = diag(k) - sbdiag(k)*wgamma(k)
1491 178970544 : xout(k) = (rhs(k) - sbdiag(k)*xout(k-1)) &
1492 1296197802 : / wbeta
1493 : enddo ! k
1494 :
1495 1206712530 : do k = nmat-1, 1, -1
1496 1206712530 : xout(k) = xout(k) - wgamma(k+1)*xout(k+1)
1497 : enddo ! k
1498 :
1499 134079170 : end subroutine tridiag_solver
1500 :
1501 : !=======================================================================
1502 :
1503 : end module icepack_therm_bl99
1504 :
1505 : !=======================================================================
|
