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