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)
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 1519223 : subroutine temperature_changes (dt, &
57 : rhoa, flw, &
58 : potT, Qa, &
59 : shcoef, lhcoef, &
60 : fswsfc, fswint, &
61 1519223 : Sswabs, Iswabs, &
62 : hilyr, hslyr, &
63 1519223 : zqin, zTin, &
64 1519223 : zqsn, zTsn, &
65 1519223 : zSin, &
66 : Tsf, Tbot, &
67 : fsensn, flatn, &
68 : flwoutn, fsurfn, &
69 : fcondtopn,fcondbot, &
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)
77 : flw , & ! incoming longwave radiation (W/m^2)
78 : potT , & ! air potential temperature (K)
79 : Qa , & ! specific humidity (kg/kg)
80 : shcoef , & ! transfer coefficient for sensible heat
81 : lhcoef , & ! transfer coefficient for latent heat
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)
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)
90 : hslyr , & ! snow layer thickness (m)
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
101 : fcondtopn , & ! downward cond flux at top surface (W m-2)
102 : fsensn , & ! surface downward sensible heat (W m-2)
103 : flatn , & ! surface downward latent heat (W m-2)
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)
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)
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
134 : niter , & ! iteration counter in temperature solver
135 : nmat ! matrix dimension
136 :
137 : logical (kind=log_kind) :: &
138 : l_snow , & ! true if snow temperatures are computed
139 : l_cold ! true if surface temperature is computed
140 :
141 : real (kind=dbl_kind) :: &
142 : Tsf_start , & ! Tsf at start of iteration
143 : dTsf , & ! Tsf - Tsf_start
144 : dTi1 , & ! Ti1(1) - Tin_start(1)
145 : dfsurf_dT , & ! derivative of fsurf wrt Tsf
146 : avg_Tsi , & ! = 1. if new snow/ice temps averaged w/starting temps
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
151 : dTi1_prev , & ! dTi1 from previous iteration
152 : dfsens_dT , & ! deriv of fsens wrt Tsf (W m-2 deg-1)
153 : dflat_dT , & ! deriv of flat wrt Tsf (W m-2 deg-1)
154 : dflwout_dT , & ! deriv of flwout wrt Tsf (W m-2 deg-1)
155 : dt_rhoi_hlyr, & ! dt/(rhoi*hilyr)
156 : einex , & ! excess energy from dqmat to ocean
157 : ferr ! energy conservation error (W m-2)
158 :
159 : real (kind=dbl_kind), dimension (nilyr) :: &
160 3038446 : Tin_init , & ! zTin at beginning of time step
161 3038446 : Tin_start , & ! zTin at start of iteration
162 3038446 : dTmat , & ! zTin - matrix solution before limiting
163 3038446 : dqmat , & ! associated enthalpy difference
164 3038446 : Tmlts ! melting temp, -depressT * salinity
165 :
166 : real (kind=dbl_kind), dimension (nslyr) :: &
167 3038446 : Tsn_init , & ! zTsn at beginning of time step
168 3038446 : Tsn_start , & ! zTsn at start of iteration
169 3038446 : etas ! dt / (rho * cp * h) for snow layers
170 :
171 : real (kind=dbl_kind), dimension (nilyr+nslyr+1) :: &
172 3038446 : sbdiag , & ! sub-diagonal matrix elements
173 1519223 : diag , & ! diagonal matrix elements
174 3038446 : spdiag , & ! super-diagonal matrix elements
175 3038446 : rhs , & ! rhs of tri-diagonal matrix equation
176 3038446 : Tmat ! matrix output temperatures
177 :
178 : real (kind=dbl_kind), dimension (nilyr) :: &
179 3038446 : etai ! dt / (rho * cp * h) for ice layers
180 :
181 : real (kind=dbl_kind), dimension(nilyr+nslyr+1):: &
182 3038446 : 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)
186 : avg_Tsf , & ! = 1. if Tsf averaged w/Tsf_start, else = 0.
187 : Iswabs_tmp , & ! energy to melt through fraction frac of layer
188 : Sswabs_tmp , & ! same for snow
189 : dswabs , & ! difference in swabs and swabs_tmp
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 1519223 : reduce_kh ! reduce conductivity when T exceeds Tmlt
197 :
198 : character(len=*),parameter :: subname='(temperature_changes)'
199 :
200 : !-----------------------------------------------------------------
201 : ! Initialize
202 : !-----------------------------------------------------------------
203 :
204 1519223 : converged = .false.
205 1519223 : l_snow = .false.
206 1519223 : l_cold = .true.
207 1519223 : fcondbot = c0
208 1519223 : dTsf_prev = c0
209 1519223 : dTi1_prev = c0
210 1519223 : dfsens_dT = c0
211 1519223 : dflat_dT = c0
212 1519223 : dflwout_dT = c0
213 1519223 : einex = c0
214 1519223 : dt_rhoi_hlyr = dt / (rhoi*hilyr) ! hilyr > 0
215 1519223 : if (hslyr > hs_min/real(nslyr,kind=dbl_kind)) &
216 885675 : l_snow = .true.
217 :
218 3944942 : do k = 1, nslyr
219 2425719 : Tsn_init (k) = zTsn(k) ! beginning of time step
220 2425719 : Tsn_start(k) = zTsn(k) ! beginning of iteration
221 3944942 : if (l_snow) then
222 1033011 : etas(k) = dt/(rhos*cp_ice*hslyr)
223 : else
224 1392708 : etas(k) = c0
225 : endif
226 : enddo ! k
227 :
228 10419676 : do k = 1, nilyr
229 8900453 : Tin_init (k) = zTin(k) ! beginning of time step
230 8900453 : Tin_start(k) = zTin(k) ! beginning of iteration
231 10419676 : 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, &
243 1519223 : zTin, kh, zSin)
244 1519223 : 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 1519223 : if (sw_redist) then
260 :
261 7026656 : do k = 1, nilyr
262 :
263 6148324 : Iswabs_tmp = c0 ! all Iswabs is moved into fswsfc
264 6148324 : if (Tin_init(k) <= Tmlts(k) - sw_dtemp) then
265 5972274 : if (l_brine) then
266 5972274 : ci = cp_ice - Lfresh * Tmlts(k) / (Tin_init(k)**2)
267 : Iswabs_tmp = min(Iswabs(k), &
268 5972274 : 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 6148324 : if (Iswabs_tmp < puny) Iswabs_tmp = c0
276 :
277 6148324 : dswabs = min(Iswabs(k) - Iswabs_tmp, fswint)
278 :
279 6148324 : fswsfc = fswsfc + dswabs
280 6148324 : fswint = fswint - dswabs
281 7026656 : Iswabs(k) = Iswabs_tmp
282 :
283 : enddo
284 :
285 2663160 : do k = 1, nslyr
286 2663160 : if (l_snow) then
287 :
288 653927 : Sswabs_tmp = c0
289 653927 : if (Tsn_init(k) <= -sw_dtemp) then
290 : Sswabs_tmp = min(Sswabs(k), &
291 601873 : -sw_frac*Tsn_init(k)/etas(k))
292 : endif
293 653927 : if (Sswabs_tmp < puny) Sswabs_tmp = c0
294 :
295 653927 : dswabs = min(Sswabs(k) - Sswabs_tmp, fswint)
296 :
297 653927 : fswsfc = fswsfc + dswabs
298 653927 : fswint = fswint - dswabs
299 653927 : 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 1519223 : converged = .false.
311 :
312 153441523 : do niter = 1, nitermax
313 :
314 : !-----------------------------------------------------------------
315 : ! Identify cells, if any, where calculation has not converged.
316 : !-----------------------------------------------------------------
317 :
318 153441523 : if (.not.converged) then
319 :
320 : !-----------------------------------------------------------------
321 : ! Allocate and initialize
322 : !-----------------------------------------------------------------
323 :
324 3143376 : converged = .true.
325 3143376 : dfsurf_dT = c0
326 3143376 : avg_Tsi = c0
327 3143376 : enew = c0
328 3143376 : 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 22793892 : do k = 1, nilyr
338 :
339 19650516 : if (l_brine) then
340 : ci = cp_ice - Lfresh*Tmlts(k) / &
341 19650516 : (zTin(k)*Tin_init(k))
342 : else
343 0 : ci = cp_ice
344 : endif
345 22793892 : etai(k) = dt_rhoi_hlyr / ci
346 :
347 : enddo
348 :
349 3143376 : 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 , &
359 : potT , Qa , &
360 : shcoef , lhcoef, &
361 : flwoutn, fsensn, &
362 3143376 : flatn , fsurfn)
363 3143376 : 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 , &
368 : dfsurf_dT, dflwout_dT, &
369 3143376 : dfsens_dT, dflat_dT )
370 3143376 : 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 3143376 : if (l_snow) then
379 1831119 : fcondtopn = kh(1) * (Tsf - zTsn(1))
380 : else
381 1312257 : fcondtopn = kh(1+nslyr) * (Tsf - zTin(1))
382 : endif
383 :
384 3143376 : if (Tsf >= c0 .and. fsurfn < fcondtopn) &
385 9851 : Tsf = -puny
386 :
387 : !-----------------------------------------------------------------
388 : ! Save surface temperature at start of iteration
389 : !-----------------------------------------------------------------
390 :
391 3143376 : Tsf_start = Tsf
392 :
393 3143376 : if (Tsf < c0) then
394 2571801 : l_cold = .true.
395 : else
396 571575 : 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, &
405 : Tsf, Tbot, &
406 : fsurfn, dfsurf_dT, &
407 : Tin_init, Tsn_init, &
408 : kh, Sswabs, &
409 : Iswabs, &
410 : etai, etas, &
411 : sbdiag, diag, &
412 3143376 : spdiag, rhs)
413 3143376 : if (icepack_warnings_aborted(subname)) return
414 :
415 : else
416 :
417 : call get_matrix_elements_know_Tsfc ( &
418 : l_snow, Tbot, &
419 : Tin_init, Tsn_init, &
420 : kh, Sswabs, &
421 : Iswabs, &
422 : etai, etas, &
423 : sbdiag, diag, &
424 : spdiag, rhs, &
425 0 : fcondtopn)
426 0 : 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 3143376 : nmat = nslyr + nilyr + 1 ! matrix dimension
435 :
436 : call tridiag_solver (nmat, sbdiag(:), &
437 : diag(:), spdiag(:), &
438 3143376 : rhs(:), Tmat(:))
439 3143376 : 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 3143376 : if (calc_Tsfc) then
469 :
470 : !-----------------------------------------------------------------
471 : ! Reload Tsf from matrix solution
472 : !-----------------------------------------------------------------
473 :
474 3143376 : if (l_cold) then
475 2571801 : if (l_snow) then
476 1787478 : Tsf = Tmat(1)
477 : else
478 784323 : Tsf = Tmat(1+nslyr)
479 : endif
480 : else ! melting surface
481 571575 : 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 3143376 : dTsf = Tsf - Tsf_start
492 3143376 : 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 3143376 : if (Tsf > puny) then
501 147215 : Tsf = c0
502 147215 : dTsf = -Tsf_start
503 147215 : if (l_brine) avg_Tsi = c1 ! avg with starting temp
504 147215 : 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 &
513 : .and. abs(dTsf) > puny &
514 : .and. abs(dTsf_prev) > puny &
515 2996161 : .and. -dTsf/(dTsf_prev+puny*puny) > p5) then
516 :
517 1372 : if (l_brine) then ! average with starting temp
518 1372 : avg_Tsf = c1
519 1372 : avg_Tsi = c1
520 : endif
521 1372 : dTsf = p5 * dTsf
522 1372 : 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 3143376 : + avg_Tsf * p5 * (Tsf_start - Tsf)
533 :
534 : endif ! calc_Tsfc
535 :
536 8289692 : do k = 1, nslyr
537 :
538 : !-----------------------------------------------------------------
539 : ! Reload zTsn from matrix solution
540 : !-----------------------------------------------------------------
541 :
542 5146316 : if (l_snow) then
543 2129463 : zTsn(k) = Tmat(k+1)
544 : else
545 3016853 : zTsn(k) = c0
546 : endif
547 5146316 : 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 5146316 : + avg_Tsi*p5*(Tsn_start(k)-zTsn(k))
555 :
556 : !-----------------------------------------------------------------
557 : ! Compute zqsn and increment new energy.
558 : !-----------------------------------------------------------------
559 5146316 : zqsn(k) = -rhos * (Lfresh - cp_ice*zTsn(k))
560 5146316 : enew = enew + hslyr * zqsn(k)
561 :
562 8289692 : Tsn_start(k) = zTsn(k) ! for next iteration
563 :
564 : enddo ! nslyr
565 :
566 22793892 : dTmat(:) = c0
567 22793892 : dqmat(:) = c0
568 22793892 : reduce_kh(:) = .false.
569 22793892 : do k = 1, nilyr
570 :
571 : !-----------------------------------------------------------------
572 : ! Reload zTin from matrix solution
573 : !-----------------------------------------------------------------
574 :
575 19650516 : zTin(k) = Tmat(k+1+nslyr)
576 :
577 19650516 : if (l_brine .and. zTin(k) > Tmlts(k) - puny) then
578 3 : dTmat(k) = zTin(k) - Tmlts(k)
579 : dqmat(k) = rhoi * dTmat(k) &
580 3 : * (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 3 : zTin(k) = Tmlts(k)
585 3 : 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 19650516 : if (k==1 .and. .not.calc_Tsfc) then
594 0 : dTi1 = zTin(k) - Tin_start(k)
595 :
596 : if (niter > 1 & ! condition 2b
597 : .and. abs(dTi1) > puny &
598 : .and. abs(dTi1_prev) > puny &
599 0 : .and. -dTi1/(dTi1_prev+puny*puny) > p5) then
600 :
601 0 : if (l_brine) avg_Tsi = c1
602 0 : dTi1 = p5 * dTi1
603 0 : converged = .false.
604 : endif
605 0 : 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 19650516 : + avg_Tsi*p5*(Tin_start(k)-zTin(k))
614 :
615 : !-----------------------------------------------------------------
616 : ! Compute zqin and increment new energy.
617 : !-----------------------------------------------------------------
618 19650516 : if (l_brine) then
619 : zqin(k) = -rhoi * (cp_ice*(Tmlts(k)-zTin(k)) &
620 : + Lfresh*(c1-Tmlts(k)/zTin(k)) &
621 19650516 : - cp_ocn*Tmlts(k))
622 : else
623 0 : zqin(k) = -rhoi * (-cp_ice*zTin(k) + Lfresh)
624 : endif
625 19650516 : enew = enew + hilyr * zqin(k)
626 19650516 : einex = einex + hilyr * dqmat(k)
627 :
628 22793892 : Tin_start(k) = zTin(k) ! for next iteration
629 :
630 : enddo ! nilyr
631 :
632 3143376 : if (calc_Tsfc) then
633 :
634 : !-----------------------------------------------------------------
635 : ! Condition 3: check for large change in Tsf
636 : !-----------------------------------------------------------------
637 :
638 3143376 : if (abs(dTsf) > Tsf_errmax) then
639 1106677 : converged = .false.
640 : endif
641 :
642 : !-----------------------------------------------------------------
643 : ! Condition 4: check for fsurfn < fcondtopn with Tsf >= 0
644 : !-----------------------------------------------------------------
645 :
646 3143376 : fsurfn = fsurfn + dTsf*dfsurf_dT
647 3143376 : if (l_snow) then
648 1831119 : fcondtopn = kh(1) * (Tsf-zTsn(1))
649 : else
650 1312257 : fcondtopn = kh(1+nslyr) * (Tsf-zTin(1))
651 : endif
652 :
653 3143376 : if (Tsf >= c0 .and. fsurfn < fcondtopn) then
654 4549 : converged = .false.
655 : endif
656 :
657 3143376 : 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 3143376 : (zTin(nilyr) - Tbot)
668 :
669 : ! Flux extra energy out of the ice
670 3143376 : fcondbot = fcondbot + einex/dt
671 :
672 : ferr = abs( (enew-einit)/dt &
673 3143376 : - (fcondtopn - fcondbot + fswint) )
674 :
675 : ! factor of 0.9 allows for roundoff errors later
676 3143376 : if (ferr > 0.9_dbl_kind*ferrmax) then ! condition (5)
677 :
678 1251978 : converged = .false.
679 :
680 : ! reduce conductivity for next iteration
681 9820830 : do k = 1, nilyr
682 9820830 : if (reduce_kh(k) .and. dqmat(k) > c0) then
683 3 : frac = max(0.5*(c1-ferr/abs(fcondtopn-fcondbot)),p1)
684 : ! frac = p1
685 3 : kh(k+nslyr+1) = kh(k+nslyr+1) * frac
686 3 : 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 1519223 : 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 1519223 : if (calc_Tsfc) then
778 :
779 : ! update fluxes that depend on Tsf
780 1519223 : flwoutn = flwoutn + dTsf_prev * dflwout_dT
781 1519223 : fsensn = fsensn + dTsf_prev * dfsens_dT
782 1519223 : 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 1519223 : subroutine conductivity (l_snow, &
799 : hilyr, hslyr, &
800 1519223 : 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)
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
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 3038446 : kilyr ! thermal cond at ice layer midpoints (W m-1 deg-1)
823 :
824 : real (kind=dbl_kind), dimension (nslyr) :: &
825 3038446 : 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 3944942 : do k = 1, nslyr
831 3944942 : kslyr(k) = ksno
832 : enddo ! nslyr
833 :
834 : ! interior ice layers
835 1519223 : if (conduct == 'MU71') then
836 : ! Maykut and Untersteiner 1971 form (with Wettlaufer 1991 constants)
837 2668480 : do k = 1, nilyr
838 2334920 : kilyr(k) = kice + betak*zSin(k)/min(-puny,zTin(k))
839 2668480 : kilyr(k) = max (kilyr(k), kimin)
840 : enddo ! nilyr
841 : else
842 : ! Pringle et al JGR 2007 'bubbly brine'
843 7751196 : 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))) &
846 6565533 : * rhoi / 917._dbl_kind
847 7751196 : 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 1519223 : if (l_snow) then
854 885675 : kh(1) = c2 * kslyr(1) / hslyr
855 : kh(1+nslyr) = c2 * kslyr(nslyr) * kilyr(1) / &
856 : ( kslyr(nslyr)*hilyr + &
857 885675 : kilyr(1 )*hslyr )
858 : else
859 633548 : kh(1) = c0
860 633548 : kh(1+nslyr) = c2 * kilyr(1) / hilyr
861 : endif
862 :
863 : ! bottom surface of bottom ice layer
864 1519223 : kh(1+nslyr+nilyr) = c2 * kilyr(nilyr) / hilyr
865 :
866 : ! interior snow interfaces
867 :
868 1519223 : if (nslyr > 1) then
869 1133120 : do k = 2, nslyr
870 1133120 : if (l_snow) then
871 : kh(k) = c2 * kslyr(k-1) * kslyr(k) / &
872 147336 : ((kslyr(k-1) + kslyr(k))*hslyr)
873 : else
874 759160 : kh(k) = c0
875 : endif
876 : enddo ! nilyr
877 : endif ! nslyr > 1
878 :
879 : ! interior ice interfaces
880 8900453 : do k = 2, nilyr
881 : kh(k+nslyr) = c2 * kilyr(k-1) * kilyr(k) / &
882 8611435 : ((kilyr(k-1) + kilyr(k))*hilyr)
883 : enddo ! nilyr
884 :
885 1519223 : 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, &
897 : potT, Qa, &
898 : shcoef, lhcoef, &
899 : flwoutn, fsensn, &
900 : flatn, fsurfn, &
901 : dflwout_dT, dfsens_dT, &
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)
909 : rhoa , & ! air density (kg/m^3)
910 : flw , & ! incoming longwave radiation (W/m^2)
911 : potT , & ! air potential temperature (K)
912 : Qa , & ! specific humidity (kg/kg)
913 : shcoef , & ! transfer coefficient for sensible heat
914 : lhcoef ! transfer coefficient for latent heat
915 :
916 : real (kind=dbl_kind), intent(inout) :: &
917 : fsensn , & ! surface downward sensible heat (W m-2)
918 : flatn , & ! surface downward latent heat (W m-2)
919 : flwoutn , & ! upward LW at surface (W m-2)
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)
924 : dflat_dT , & ! deriv of flat wrt Tsf (W m-2 deg-1)
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, &
935 : potT, Qa, &
936 : shcoef, lhcoef, &
937 : flwoutn, fsensn, &
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, &
944 : dfsurf_dT, dflwout_dT, &
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 3143376 : subroutine get_matrix_elements_calc_Tsfc ( &
963 : l_snow, l_cold, &
964 : Tsf, Tbot, &
965 : fsurfn, dfsurf_dT, &
966 0 : Tin_init, Tsn_init, &
967 3143376 : kh, Sswabs, &
968 3143376 : Iswabs, &
969 6286752 : etai, etas, &
970 3143376 : sbdiag, diag, &
971 3143376 : spdiag, rhs)
972 :
973 : logical (kind=log_kind), intent(in) :: &
974 : l_snow , & ! true if snow temperatures are computed
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)
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
989 : Tin_init , & ! ice temp at beginning of time step
990 : Sswabs , & ! SW radiation absorbed in snow layers (W m-2)
991 : Iswabs , & ! absorbed SW flux in ice layers
992 : etas , & ! dt / (rho*cp*h) for snow layers
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
1001 : diag , & ! diagonal matrix elements
1002 : spdiag , & ! super-diagonal matrix elements
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 11433068 : do k = 1, nslyr+1
1020 8289692 : sbdiag(k) = c0
1021 8289692 : diag (k) = c1
1022 8289692 : spdiag(k) = c0
1023 11433068 : 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 3143376 : if (l_cold) then
1041 2571801 : if (l_snow) then
1042 1787478 : k = 1
1043 : else ! no snow
1044 784323 : k = 1 + nslyr
1045 : endif
1046 2571801 : kr = k
1047 :
1048 2571801 : sbdiag(kr) = c0
1049 2571801 : diag (kr) = dfsurf_dT - kh(k)
1050 2571801 : spdiag(kr) = kh(k)
1051 2571801 : 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 3143376 : if (l_snow) then
1061 1831119 : if (l_cold) then
1062 1787478 : sbdiag(2) = -etas(1) * kh(1)
1063 1787478 : spdiag(2) = -etas(1) * kh(2)
1064 : diag (2) = c1 &
1065 1787478 : + etas(1) * (kh(1) + kh(2))
1066 : rhs (2) = Tsn_init(1) &
1067 1787478 : + etas(1) * Sswabs(1)
1068 : else ! melting surface
1069 43641 : sbdiag(2) = c0
1070 43641 : spdiag(2) = -etas(1) * kh(2)
1071 : diag (2) = c1 &
1072 43641 : + etas(1) * (kh(1) + kh(2))
1073 : rhs (2) = Tsn_init(1) &
1074 : + etas(1)*kh(1)*Tsf &
1075 43641 : + etas(1) * Sswabs(1)
1076 : endif ! l_cold
1077 : endif ! l_snow
1078 :
1079 : !-----------------------------------------------------------------
1080 : ! remaining snow layers
1081 : !-----------------------------------------------------------------
1082 :
1083 3143376 : if (nslyr > 1) then
1084 :
1085 2503675 : do k = 2, nslyr
1086 2002940 : kr = k + 1
1087 :
1088 2503675 : if (l_snow) then
1089 298344 : sbdiag(kr) = -etas(k) * kh(k)
1090 298344 : spdiag(kr) = -etas(k) * kh(k+1)
1091 : diag (kr) = c1 &
1092 298344 : + etas(k) * (kh(k) + kh(k+1))
1093 : rhs (kr) = Tsn_init(k) &
1094 298344 : + etas(k) * Sswabs(k)
1095 : endif
1096 : enddo ! nslyr
1097 :
1098 : endif ! nslyr > 1
1099 :
1100 3143376 : if (nilyr > 1) then
1101 :
1102 : !-----------------------------------------------------------------
1103 : ! top ice layer
1104 : !-----------------------------------------------------------------
1105 :
1106 2751190 : ki = 1
1107 2751190 : k = ki + nslyr
1108 2751190 : kr = k + 1
1109 :
1110 2751190 : if (l_snow .or. l_cold) then
1111 2223256 : sbdiag(kr) = -etai(ki) * kh(k)
1112 2223256 : spdiag(kr) = -etai(ki) * kh(k+1)
1113 : diag (kr) = c1 &
1114 2223256 : + etai(ki) * (kh(k) + kh(k+1))
1115 : rhs (kr) = Tin_init(ki) &
1116 2223256 : + etai(ki)*Iswabs(ki)
1117 : else ! no snow, warm surface
1118 527934 : sbdiag(kr) = c0
1119 527934 : spdiag(kr) = -etai(ki) * kh(k+1)
1120 : diag (kr) = c1 &
1121 527934 : + etai(ki) * (kh(k) + kh(k+1))
1122 : rhs (kr) = Tin_init(ki) &
1123 : + etai(ki)*Iswabs(ki) &
1124 527934 : + etai(ki)*kh(k)*Tsf
1125 : endif
1126 :
1127 : !-----------------------------------------------------------------
1128 : ! bottom ice layer
1129 : !-----------------------------------------------------------------
1130 :
1131 2751190 : ki = nilyr
1132 2751190 : k = ki + nslyr
1133 2751190 : kr = k + 1
1134 :
1135 2751190 : sbdiag(kr) = -etai(ki) * kh(k)
1136 2751190 : spdiag(kr) = c0
1137 : diag (kr) = c1 &
1138 2751190 : + etai(ki) * (kh(k) + kh(k+1))
1139 : rhs (kr) = Tin_init(ki) &
1140 : + etai(ki)*Iswabs(ki) &
1141 2751190 : + etai(ki)*kh(k+1)*Tbot
1142 :
1143 : else ! nilyr = 1
1144 :
1145 : !-----------------------------------------------------------------
1146 : ! single ice layer
1147 : !-----------------------------------------------------------------
1148 :
1149 392186 : ki = 1
1150 392186 : k = ki + nslyr
1151 392186 : kr = k + 1
1152 :
1153 392186 : if (l_snow .or. l_cold) then
1154 392186 : sbdiag(kr) = -etai(ki) * kh(k)
1155 392186 : spdiag(kr) = c0
1156 : diag (kr) = c1 &
1157 392186 : + etai(ki) * (kh(k) + kh(k+1))
1158 : rhs (kr) = Tin_init(ki) &
1159 : + etai(ki) * Iswabs(ki) &
1160 392186 : + 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) &
1168 : + etai(ki) * kh(k)*Tsf &
1169 0 : + etai(ki) * kh(k+1)*Tbot
1170 : endif
1171 :
1172 : endif ! nilyr > 1
1173 :
1174 : !-----------------------------------------------------------------
1175 : ! interior ice layers
1176 : !-----------------------------------------------------------------
1177 :
1178 16899326 : do ki = 2, nilyr-1
1179 :
1180 13755950 : k = ki + nslyr
1181 13755950 : kr = k + 1
1182 :
1183 13755950 : sbdiag(kr) = -etai(ki) * kh(k)
1184 13755950 : spdiag(kr) = -etai(ki) * kh(k+1)
1185 : diag (kr) = c1 &
1186 13755950 : + etai(ki) * (kh(k) + kh(k+1))
1187 : rhs (kr) = Tin_init(ki) &
1188 16507140 : + etai(ki)*Iswabs(ki)
1189 : enddo ! nilyr
1190 :
1191 3143376 : 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 0 : subroutine get_matrix_elements_know_Tsfc ( &
1206 : l_snow, Tbot, &
1207 0 : Tin_init, Tsn_init, &
1208 0 : kh, Sswabs, &
1209 0 : Iswabs, &
1210 0 : etai, etas, &
1211 0 : sbdiag, diag, &
1212 0 : spdiag, rhs, &
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
1223 : Tin_init , & ! ice temp at beginning of time step
1224 : Sswabs , & ! SW radiation absorbed in snow layers (W m-2)
1225 : Iswabs , & ! absorbed SW flux in ice layers
1226 : etas , & ! dt / (rho*cp*h) for snow layers
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
1235 : diag , & ! diagonal matrix elements
1236 : spdiag , & ! super-diagonal matrix elements
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 0 : do k = 1, nslyr+1
1257 0 : sbdiag(k) = c0
1258 0 : diag (k) = c1
1259 0 : spdiag(k) = c0
1260 0 : 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 0 : if (l_snow) then
1280 0 : sbdiag(2) = c0
1281 0 : spdiag(2) = -etas(1) * kh(2)
1282 : diag (2) = c1 &
1283 0 : + etas(1) * kh(2)
1284 : rhs (2) = Tsn_init(1) &
1285 : + etas(1) * Sswabs(1) &
1286 0 : + etas(1) * fcondtopn
1287 : endif ! l_snow
1288 :
1289 : !-----------------------------------------------------------------
1290 : ! remaining snow layers
1291 : !-----------------------------------------------------------------
1292 :
1293 0 : 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 0 : if (nilyr > 1) then
1312 :
1313 : !-----------------------------------------------------------------
1314 : ! top ice layer
1315 : !-----------------------------------------------------------------
1316 :
1317 0 : ki = 1
1318 0 : k = ki + nslyr
1319 0 : kr = k + 1
1320 :
1321 0 : if (l_snow) then
1322 :
1323 0 : sbdiag(kr) = -etai(ki) * kh(k)
1324 0 : spdiag(kr) = -etai(ki) * kh(k+1)
1325 : diag (kr) = c1 &
1326 0 : + etai(ki) * (kh(k) + kh(k+1))
1327 : rhs (kr) = Tin_init(ki) &
1328 0 : + etai(ki) * Iswabs(ki)
1329 : else
1330 0 : sbdiag(kr) = c0
1331 0 : spdiag(kr) = -etai(ki) * kh(k+1)
1332 : diag (kr) = c1 &
1333 0 : + etai(ki) * kh(k+1)
1334 : rhs (kr) = Tin_init(ki) &
1335 : + etai(ki) * Iswabs(ki) &
1336 0 : + etai(ki) * fcondtopn
1337 : endif ! l_snow
1338 :
1339 : !-----------------------------------------------------------------
1340 : ! bottom ice layer
1341 : !-----------------------------------------------------------------
1342 :
1343 0 : ki = nilyr
1344 0 : k = ki + nslyr
1345 0 : kr = k + 1
1346 :
1347 0 : sbdiag(kr) = -etai(ki) * kh(k)
1348 0 : spdiag(kr) = c0
1349 : diag (kr) = c1 &
1350 0 : + etai(ki) * (kh(k) + kh(k+1))
1351 : rhs (kr) = Tin_init(ki) &
1352 : + etai(ki)*Iswabs(ki) &
1353 0 : + 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) &
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) &
1380 : + etai(ki) * fcondtopn &
1381 0 : + etai(ki) * kh(k+1)*Tbot
1382 : endif
1383 :
1384 : endif ! nilyr > 1
1385 :
1386 : !-----------------------------------------------------------------
1387 : ! interior ice layers
1388 : !-----------------------------------------------------------------
1389 :
1390 0 : do ki = 2, nilyr-1
1391 :
1392 0 : k = ki + nslyr
1393 0 : kr = k + 1
1394 :
1395 0 : sbdiag(kr) = -etai(ki) * kh(k)
1396 0 : spdiag(kr) = -etai(ki) * kh(k+1)
1397 : diag (kr) = c1 &
1398 0 : + etai(ki) * (kh(k) + kh(k+1))
1399 : rhs (kr) = Tin_init(ki) &
1400 0 : + etai(ki)*Iswabs(ki)
1401 :
1402 : enddo ! nilyr
1403 :
1404 0 : 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 6286752 : diag, spdiag, &
1416 3143376 : 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
1423 : diag , & ! diagonal matrix elements
1424 : spdiag , & ! super-diagonal matrix elements
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 6286752 : wgamma ! temporary matrix variable
1440 :
1441 : character(len=*),parameter :: subname='(tridiag_solver)'
1442 :
1443 3143376 : wbeta = diag(1)
1444 3143376 : xout(1) = rhs(1) / wbeta
1445 :
1446 27940208 : do k = 2, nmat
1447 24796832 : wgamma(k) = spdiag(k-1) / wbeta
1448 24796832 : wbeta = diag(k) - sbdiag(k)*wgamma(k)
1449 : xout(k) = (rhs(k) - sbdiag(k)*xout(k-1)) &
1450 27940208 : / wbeta
1451 : enddo ! k
1452 :
1453 27940208 : do k = nmat-1, 1, -1
1454 27940208 : xout(k) = xout(k) - wgamma(k+1)*xout(k+1)
1455 : enddo ! k
1456 :
1457 3143376 : end subroutine tridiag_solver
1458 :
1459 : !=======================================================================
1460 :
1461 : end module icepack_therm_bl99
1462 :
1463 : !=======================================================================
|
