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