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