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