1 | SUBROUTINE ice_th_diff(nlay_s,nlay_i,kideb,kiut) |
---|
2 | !!------------------------------------------------------------------ |
---|
3 | !! *** ROUTINE ice_th_diff *** |
---|
4 | !! ** Purpose : |
---|
5 | !! This routine determines the time evolution of snow and sea-ice |
---|
6 | !! temperature profiles. |
---|
7 | !! ** Method : |
---|
8 | !! This is done by solving the heat equation diffusion with |
---|
9 | !! a Neumann boundary condition at the surface and a Dirichlet one |
---|
10 | !! at the bottom. Solar radiation is partially absorbed into the ice. |
---|
11 | !! The specific heat and thermal conductivities depend on ice salinity |
---|
12 | !! and temperature to take into account brine pocket melting. The |
---|
13 | !! numerical |
---|
14 | !! scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid |
---|
15 | !! in the ice and snow system. |
---|
16 | !! The successive steps of this routine are |
---|
17 | !! Vertical grid |
---|
18 | !! 1. Thermal conductivity at the interfaces of the ice layers |
---|
19 | !! 2. Internal absorbed radiation |
---|
20 | !! 3. Scale factors due to non-uniform grid |
---|
21 | !! 4. Kappa factors |
---|
22 | !! Then iterative procedure begins |
---|
23 | !! 5. specific heat in the ice |
---|
24 | !! 6. eta factors |
---|
25 | !! 7. surface flux computation |
---|
26 | !! 8. tridiagonal system terms |
---|
27 | !! 9. solving the tridiagonal system with Gauss elimination |
---|
28 | !! Iterative procedure ends according to a criterion on evolution |
---|
29 | !! of temperature |
---|
30 | !! |
---|
31 | !! ** Arguments : |
---|
32 | !! kideb , kiut : Starting and ending points on which the |
---|
33 | !! the computation is applied |
---|
34 | !! |
---|
35 | !! ** Inputs / Ouputs : (global commons) |
---|
36 | !! surface temperature : t_su_b |
---|
37 | !! ice/snow temperatures : t_i_b, t_s_b |
---|
38 | !! ice salinities : s_i_b |
---|
39 | !! number of layers in the ice/snow: nlay_i, nlay_s |
---|
40 | !! total ice/snow thickness : ht_i_b, ht_s_b |
---|
41 | !! |
---|
42 | !! ** External : |
---|
43 | !! |
---|
44 | !! ** References : |
---|
45 | !! |
---|
46 | !! ** History : |
---|
47 | !! (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium |
---|
48 | !! |
---|
49 | |
---|
50 | USE lib_fortran |
---|
51 | |
---|
52 | INCLUDE 'type.com' |
---|
53 | INCLUDE 'para.com' |
---|
54 | INCLUDE 'const.com' |
---|
55 | INCLUDE 'ice.com' |
---|
56 | INCLUDE 'thermo.com' |
---|
57 | |
---|
58 | ! Local variables |
---|
59 | INTEGER numeqmin, numeqmax, numeq |
---|
60 | DIMENSION ztcond_i(0:nlay_i), |
---|
61 | & zkappa_s(0:nlay_s), |
---|
62 | & zkappa_i(0:nlay_i),ztstemp(0:nlay_s),ztitemp(0:nlay_i), |
---|
63 | & zspeche_i(0:nlay_i),ztsold(0:nlay_s),ztiold(0:nlay_i), |
---|
64 | & zeta_s(0:nlay_s),zeta_i(0:nlay_i),ztrid(2*maxnlay+1,3), |
---|
65 | & zindterm(2*maxnlay+1),zindtbis(2*maxnlay+1), |
---|
66 | & zdiagbis(2*maxnlay+1),ykn(nbpt),ykg(nbpt), |
---|
67 | & zrchu1(nbpt),zrchu2(nbpt),zqsat(nbpt) |
---|
68 | |
---|
69 | LOGICAL :: ln_write |
---|
70 | |
---|
71 | ! Local constants |
---|
72 | zeps = 1.0d-20 |
---|
73 | zg1s = 2.0 |
---|
74 | zg1 = 2.0 |
---|
75 | zbeta = 0.117 ! factor for thermal conductivity |
---|
76 | zerrmax = 1.0d-11 ! max error at the surface |
---|
77 | |
---|
78 | ! new lines |
---|
79 | zerrmax = 1.0e-4 |
---|
80 | nconv_max = 50 |
---|
81 | |
---|
82 | ln_write = .TRUE. |
---|
83 | |
---|
84 | DO 5 ji = kideb, kiut |
---|
85 | |
---|
86 | IF ( ln_write ) THEN |
---|
87 | WRITE(numout,*) ' ** ice_th_diff : ' |
---|
88 | WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ ' |
---|
89 | WRITE(numout,*) ' nlay_i: ', nlay_i |
---|
90 | WRITE(numout,*) ' nlay_s: ', nlay_s |
---|
91 | WRITE(numout,*) ' t_su_b: ', t_su_b(ji) |
---|
92 | WRITE(numout,*) ' t_s_b : ', ( t_s_b(ji,layer), |
---|
93 | & layer = 1, nlay_s ) |
---|
94 | WRITE(numout,*) ' t_i_b : ', ( t_i_b(ji,layer), |
---|
95 | & layer = 1, nlay_i ) |
---|
96 | WRITE(numout,*) ' t_bo_b : ', t_bo_b(ji) |
---|
97 | WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) |
---|
98 | WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) |
---|
99 | ENDIF |
---|
100 | |
---|
101 | ! Switches |
---|
102 | ! isnow equals 1 if snow is present and 0 if absent |
---|
103 | isnow = int(1.0-max(0.0,sign(1.0d0,-ht_s_b(ji)))) |
---|
104 | ! imelt equals 1 if surface is melting and 0 if not |
---|
105 | imelt = int(max(0.0,sign(1.0d0,t_su_b(ji)-tpw))) |
---|
106 | |
---|
107 | tfs(ji) = real(isnow)*tfsn+ (1.0-real(isnow))*tfsg |
---|
108 | |
---|
109 | ! Oceanic heat flux and precipitations |
---|
110 | fbbqb(ji) = oce_flx |
---|
111 | |
---|
112 | IF ( ln_write ) THEN |
---|
113 | WRITE(numout,*) ' isnow : ', isnow |
---|
114 | WRITE(numout,*) ' imelt : ', imelt |
---|
115 | WRITE(numout,*) ' tfs : ', tfs(ji) |
---|
116 | ENDIF |
---|
117 | |
---|
118 | 5 CONTINUE |
---|
119 | ! |
---|
120 | !------------------------------------------------------------------------------ |
---|
121 | ! 1) Thermal conductivity at the ice interfaces |
---|
122 | !------------------------------------------------------------------------------ |
---|
123 | ! |
---|
124 | ! Pringle et al., JGR 2007 formula |
---|
125 | ! 2.11 + 0.09 S/T - 0.011.T |
---|
126 | |
---|
127 | DO 10 ji = kideb, kiut |
---|
128 | |
---|
129 | ! thermal conductivity in the snow |
---|
130 | ykn(ji) = xkn |
---|
131 | zkimin = 0.1 |
---|
132 | |
---|
133 | ztcond_i(0) = xkg + betak1*s_i_b(ji,1) |
---|
134 | & / MIN( -zeps , t_i_b(ji,1) - tpw ) |
---|
135 | & - betak2* ( t_i_b(ji,1) - tpw ) |
---|
136 | ztcond_i(0) = MAX( ztcond_i(0) , zkimin ) |
---|
137 | |
---|
138 | DO layer = 1, nlay_i-1 |
---|
139 | ztcond_i(layer) = xkg + betak1*( s_i_b(ji,layer) |
---|
140 | & + s_i_b(ji,layer+1) ) / MIN(-zeps, |
---|
141 | & t_i_b(ji,layer)+t_i_b(ji,layer+1)-2.0*tpw) |
---|
142 | & - betak2 * 0.5 * ( t_i_b(ji,layer) + ! bugfix fred dupont add 0.5 |
---|
143 | & t_i_b(ji,layer+1) - 2.0*tpw ) |
---|
144 | ztcond_i(layer) = MAX( ztcond_i(layer) , zkimin ) |
---|
145 | END DO |
---|
146 | |
---|
147 | ztcond_i(nlay_i) = xkg + betak1*s_i_b(ji,nlay_i) / |
---|
148 | & MIN( -zeps , t_bo_b(ji) - tpw ) |
---|
149 | & - betak2 * ( t_bo_b(ji) - tpw ) |
---|
150 | ztcond_i(nlay_i) = MAX( ztcond_i(nlay_i) , zkimin ) |
---|
151 | |
---|
152 | IF ( ln_write ) THEN |
---|
153 | WRITE(numout,*) ' ztcond_i : ', ztcond_i(0:nlay_i) |
---|
154 | ENDIF |
---|
155 | |
---|
156 | 10 CONTINUE |
---|
157 | ! |
---|
158 | !------------------------------------------------------------------------------ |
---|
159 | ! 3) kappa factors |
---|
160 | !------------------------------------------------------------------------------ |
---|
161 | ! |
---|
162 | DO 30 ji = kideb, kiut |
---|
163 | ! snow |
---|
164 | zkappa_s(0) = ykn(ji)/max(zeps,deltaz_s_phy(1)) |
---|
165 | do layer = 1, nlay_s-1 |
---|
166 | zkappa_s(layer) = 2.0*ykn(ji)/ |
---|
167 | & max(zeps,deltaz_s_phy(layer)+deltaz_s_phy(layer+1)) |
---|
168 | end do |
---|
169 | zkappa_s(nlay_s) = ykn(ji)/max(zeps,deltaz_s_phy(nlay_s)) |
---|
170 | |
---|
171 | ! ice |
---|
172 | zkappa_i(0) = ztcond_i(0)/max(zeps,deltaz_i_phy(1)) |
---|
173 | do layer = 1, nlay_i-1 |
---|
174 | zkappa_i(layer) = 2.0*ztcond_i(layer)/ |
---|
175 | & max(zeps,deltaz_i_phy(layer)+deltaz_i_phy(layer+1)) |
---|
176 | end do |
---|
177 | zkappa_i(nlay_i) = ztcond_i(nlay_i) / |
---|
178 | & MAX(zeps,deltaz_i_phy(nlay_i)) |
---|
179 | |
---|
180 | ! interface |
---|
181 | zkappa_s(nlay_s) = 2.0*ykn(ji)*ztcond_i(0)/max(zeps, |
---|
182 | & (ztcond_i(0)*deltaz_s_phy(nlay_s) + |
---|
183 | & ykn(ji)*deltaz_i_phy(1))) |
---|
184 | zkappa_i(0) = zkappa_s(nlay_s)*real(isnow) |
---|
185 | & + zkappa_i(0)*(1.0-real(isnow)) |
---|
186 | |
---|
187 | IF ( ln_write ) THEN |
---|
188 | WRITE(numout,*) ' nlay_s : ', nlay_s |
---|
189 | WRITE(numout,*) ' zkappa_s : ', zkappa_s(0:nlay_s) |
---|
190 | WRITE(numout,*) ' zkappa_i : ', zkappa_i(0:nlay_i) |
---|
191 | ENDIF |
---|
192 | |
---|
193 | 30 CONTINUE |
---|
194 | ! |
---|
195 | !------------------------------------------------------------------------------| |
---|
196 | ! 4) iterative procedure begins | |
---|
197 | !------------------------------------------------------------------------------| |
---|
198 | ! |
---|
199 | DO 40 ji = kideb, kiut |
---|
200 | |
---|
201 | !------------------------------ |
---|
202 | ! keeping old values in memory |
---|
203 | !------------------------------ |
---|
204 | |
---|
205 | ztsuold = t_su_b(ji) |
---|
206 | t_su_b(ji) = min(t_su_b(ji),tfs(ji)-0.00001) |
---|
207 | DO layer = 1, nlay_s |
---|
208 | ztsold(layer) = t_s_b(ji,layer) |
---|
209 | END DO |
---|
210 | DO layer = 1, nlay_i |
---|
211 | ztiold(layer) = t_i_b(ji,layer) |
---|
212 | ti_old(layer) = ztiold(layer) |
---|
213 | END DO |
---|
214 | |
---|
215 | nconv = 0 |
---|
216 | |
---|
217 | 40 CONTINUE |
---|
218 | |
---|
219 | !------------------------------ |
---|
220 | ! Beginning of the loop |
---|
221 | !------------------------------ |
---|
222 | |
---|
223 | zerrit = 10000.0 |
---|
224 | |
---|
225 | ! DO WHILE ( zerrit .GT. zerrmax ) |
---|
226 | DO WHILE ( ( zerrit .GT. zerrmax ) .AND. ( nconv < nconv_max ) ) |
---|
227 | |
---|
228 | do 42 ji = kideb, kiut |
---|
229 | |
---|
230 | !45 CONTINUE |
---|
231 | |
---|
232 | nconv = nconv+1 |
---|
233 | |
---|
234 | ztsutemp = t_su_b(ji) |
---|
235 | DO layer = 1, nlay_s |
---|
236 | ztstemp(layer) = t_s_b(ji,layer) |
---|
237 | END DO |
---|
238 | DO layer = 1, nlay_i |
---|
239 | ztitemp(layer) = t_i_b(ji,layer) |
---|
240 | END DO |
---|
241 | |
---|
242 | IF ( ln_write ) THEN |
---|
243 | WRITE(numout,*) ' zerrit : ', zerrit |
---|
244 | WRITE(numout,*) ' nconv : ', nconv |
---|
245 | WRITE(numout,*) ' t_s_b : ', t_s_b(ji,1) |
---|
246 | WRITE(numout,*) ' t_i_b : ', ( t_i_b(ji,layer), layer = 1, |
---|
247 | & nlay_i ) |
---|
248 | ENDIF |
---|
249 | |
---|
250 | ! |
---|
251 | !------------------------------------------------------------------------------| |
---|
252 | ! 5) specific heat in the ice | |
---|
253 | !------------------------------------------------------------------------------| |
---|
254 | ! |
---|
255 | DO layer = 1, nlay_i |
---|
256 | zspeche_i(layer) = cpg + lfus*tmut*s_i_b(ji,layer)/ |
---|
257 | & MAX((t_i_b(ji,layer)-tpw)*(ztiold(layer)-tpw),zeps) |
---|
258 | END DO |
---|
259 | ! |
---|
260 | !------------------------------------------------------------------------------| |
---|
261 | ! 6) eta factors | |
---|
262 | !------------------------------------------------------------------------------| |
---|
263 | ! |
---|
264 | DO layer = 1, nlay_s |
---|
265 | zeta_s(layer) = ddtb / max(rhon*cpg*deltaz_s_phy(layer),zeps) |
---|
266 | END DO |
---|
267 | |
---|
268 | DO layer = 1, nlay_i |
---|
269 | zeta_i(layer) = ddtb / max(rhog*deltaz_i_phy(layer) * |
---|
270 | & zspeche_i(layer) |
---|
271 | & ,zeps) |
---|
272 | END DO |
---|
273 | |
---|
274 | ! |
---|
275 | !------------------------------------------------------------------------------| |
---|
276 | ! 7) surface flux computation | |
---|
277 | !------------------------------------------------------------------------------| |
---|
278 | ! |
---|
279 | t_su_b(ji) = min(t_su_b(ji),tfs(ji)) |
---|
280 | imelt = int(max(0.0,sign(1.0d0,t_su_b(ji)-tpw))) |
---|
281 | ! pressure of water vapor saturation (Pa) |
---|
282 | es = 611.0*10.0**(9.5*(t_su_b(ji)-273.16)/ |
---|
283 | & (t_su_b(ji)-7.66)) |
---|
284 | ! net longwave radiative flux |
---|
285 | ! MV BUG |
---|
286 | ! fratsb(ji) = emig*(ratbqb(ji)- |
---|
287 | ! & stefan*t_su_b(ji)*t_su_b(ji)*t_su_b(ji)*t_su_b(ji)) |
---|
288 | fratsb(ji) = ratbqb(ji) - emig * |
---|
289 | & stefan*t_su_b(ji)*t_su_b(ji)*t_su_b(ji)*t_su_b(ji) |
---|
290 | ! sensible and latent heat flux |
---|
291 | CALL flx(ht_i_b(ji),t_su_b(ji),tabqb(ji),qabqb(ji),fcsb(ji), |
---|
292 | & fleb(ji),qsfcb(ji),zrchu1(ji),zrchu2(ji),vabqb(ji),zref) |
---|
293 | fcsb(ji) = -fcsb(ji) |
---|
294 | fleb(ji) = MIN( -fleb(ji) , 0.0 ) ! always negative, as precip |
---|
295 | ! energy already added |
---|
296 | |
---|
297 | ! qsfcb(ji) = zqsat(ji) |
---|
298 | |
---|
299 | ! intermediate variable |
---|
300 | zssdqw = qsfcb(ji)*qsfcb(ji)*psbqb(ji)/ |
---|
301 | & (0.622*es)*alog(10.0)*9.5* |
---|
302 | & ((273.16-7.66)/(t_su_b(ji)-7.66)**2.0) |
---|
303 | ! derivative of the surface atmospheric net flux |
---|
304 | dzf = -4.0*emig*stefan*t_su_b(ji) |
---|
305 | & *t_su_b(ji)*t_su_b(ji) |
---|
306 | & -(zrchu1(ji)+zrchu2(ji)*zssdqw) |
---|
307 | ! surface atmospheric net flux |
---|
308 | zf = ab(ji)*fsolgb(ji)+fratsb(ji) |
---|
309 | & +fcsb(ji)+fleb(ji) |
---|
310 | |
---|
311 | ! |
---|
312 | !------------------------------------------------------------------------------| |
---|
313 | ! 8) tridiagonal system terms | |
---|
314 | !------------------------------------------------------------------------------| |
---|
315 | ! |
---|
316 | ! layer denotes the number of the layer in the snow or in the ice |
---|
317 | ! numeq denotes the reference number of the equation in the tridiagonal |
---|
318 | ! system, terms of tridiagonal system are indexed as following : |
---|
319 | ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one |
---|
320 | |
---|
321 | ! ice interior terms (top equation has the same form as the others) |
---|
322 | DO numeq = 1, maxnlay |
---|
323 | ztrid(numeq,1) = 0.0 |
---|
324 | ztrid(numeq,2) = 0.0 |
---|
325 | ztrid(numeq,3) = 0.0 |
---|
326 | zindterm(numeq) = 0.0 |
---|
327 | zindtbis(numeq) = 0.0 |
---|
328 | zdiagbis(numeq) = 0.0 |
---|
329 | END DO |
---|
330 | DO numeq = nlay_s + 2, nlay_s + nlay_i |
---|
331 | layer = numeq - nlay_s - 1 |
---|
332 | ztrid(numeq,1) = - zeta_i(layer)*zkappa_i(layer-1) |
---|
333 | ztrid(numeq,2) = 1.0 + zeta_i(layer)*(zkappa_i(layer-1) + |
---|
334 | & zkappa_i(layer)) |
---|
335 | ztrid(numeq,3) = - zeta_i(layer)*zkappa_i(layer) |
---|
336 | zindterm(numeq) = ztiold(layer) + zeta_i(layer)* |
---|
337 | ! & ( radab_phy_i(layer) + radab_alg_i(layer) ) |
---|
338 | & radab_i(layer) |
---|
339 | END DO |
---|
340 | ! ice bottom terms |
---|
341 | numeq = nlay_s + nlay_i + 1 |
---|
342 | ztrid(numeq,1) = - zeta_i(nlay_i)*zkappa_i(nlay_i-1) |
---|
343 | ztrid(numeq,2) = 1.0 + zeta_i(nlay_i)*( zkappa_i(nlay_i)*zg1 |
---|
344 | & + zkappa_i(nlay_i-1) ) |
---|
345 | ztrid(numeq,3) = 0.0 |
---|
346 | zindterm(numeq) = ztiold(nlay_i) + zeta_i(nlay_i)* |
---|
347 | ! & ( radab_phy_i(nlay_i) + radab_alg_i(nlay_i) |
---|
348 | & ( radab_i(nlay_i) |
---|
349 | & + zkappa_i(nlay_i)*zg1 |
---|
350 | & *t_bo_b(ji) ) |
---|
351 | |
---|
352 | IF (ht_s_b(ji).GT.0.0) THEN |
---|
353 | ! |
---|
354 | !------------------------------------------------------------------------------| |
---|
355 | ! snow-covered cells | |
---|
356 | !------------------------------------------------------------------------------| |
---|
357 | ! |
---|
358 | ! snow interior terms (bottom equation has the same form as the others) |
---|
359 | do numeq = 3, nlay_s + 1 |
---|
360 | layer = numeq - 1 |
---|
361 | ztrid(numeq,1) = - zeta_s(layer)*zkappa_s(layer-1) |
---|
362 | ztrid(numeq,2) = 1.0 + zeta_s(layer)*( zkappa_s(layer-1) + |
---|
363 | & zkappa_s(layer) ) |
---|
364 | ztrid(numeq,3) = - zeta_s(layer)*zkappa_s(layer) |
---|
365 | zindterm(numeq) = ztsold(layer) + zeta_s(layer)* |
---|
366 | & radab_s(layer) |
---|
367 | end do |
---|
368 | |
---|
369 | ! case of only one layer in the ice (ice equation is altered) |
---|
370 | if (nlay_i.eq.1) then |
---|
371 | ztrid(nlay_s+2,3) = 0.0 |
---|
372 | zindterm(nlay_s+2) = zindterm(nlay_s+2) + zkappa_i(1)* |
---|
373 | & t_bo_b(ji) |
---|
374 | endif |
---|
375 | |
---|
376 | IF (t_su_b(ji).LT.tpw) THEN |
---|
377 | ! |
---|
378 | !------------------------------------------------------------------------------| |
---|
379 | ! case 1 : no surface melting - snow present | |
---|
380 | !------------------------------------------------------------------------------| |
---|
381 | ! |
---|
382 | zdifcase = 1.0 |
---|
383 | numeqmin = 1 |
---|
384 | numeqmax = nlay_i + nlay_s + 1 |
---|
385 | |
---|
386 | ! surface equation |
---|
387 | ztrid(1,1) = 0.0 |
---|
388 | ztrid(1,2) = dzf - zg1s*zkappa_s(0) |
---|
389 | ztrid(1,3) = zg1s*zkappa_s(0) |
---|
390 | zindterm(1) = dzf*t_su_b(ji) - zf |
---|
391 | |
---|
392 | ! first layer of snow equation |
---|
393 | ztrid(2,1) = - zkappa_s(0)*zg1s*zeta_s(1) |
---|
394 | ztrid(2,2) = 1.0 + zeta_s(1)*(zkappa_s(1) + zkappa_s(0)*zg1s) |
---|
395 | ztrid(2,3) = - zeta_s(1)* zkappa_s(1) |
---|
396 | zindterm(2) = ztsold(1) + zeta_s(1)*radab_s(1) |
---|
397 | |
---|
398 | else |
---|
399 | ! |
---|
400 | !------------------------------------------------------------------------------| |
---|
401 | ! case 2 : surface is melting - snow present | |
---|
402 | !------------------------------------------------------------------------------| |
---|
403 | ! |
---|
404 | zdifcase = 2.0 |
---|
405 | numeqmin = 2 |
---|
406 | numeqmax = nlay_i + nlay_s + 1 |
---|
407 | |
---|
408 | ! first layer of snow equation |
---|
409 | ztrid(2,1) = 0.0 |
---|
410 | ztrid(2,2) = 1.0 + zeta_s(1)*(zkappa_s(1) + zkappa_s(0)*zg1s) |
---|
411 | ztrid(2,3) = - zeta_s(1)*zkappa_s(1) |
---|
412 | zindterm(2) = ztsold(1) + zeta_s(1)*(radab_s(1) + zkappa_s(0)* |
---|
413 | & zg1s*t_su_b(ji)) |
---|
414 | |
---|
415 | ENDIF |
---|
416 | ELSE |
---|
417 | ! |
---|
418 | !------------------------------------------------------------------------------| |
---|
419 | ! cells without snow | |
---|
420 | !------------------------------------------------------------------------------| |
---|
421 | ! |
---|
422 | IF (t_su_b(ji).lt.tpw) THEN |
---|
423 | ! |
---|
424 | !------------------------------------------------------------------------------| |
---|
425 | ! case 3 : no surface melting - no snow | |
---|
426 | !------------------------------------------------------------------------------| |
---|
427 | ! |
---|
428 | zdifcase = 3.0 |
---|
429 | numeqmin = nlay_s + 1 |
---|
430 | numeqmax = nlay_i + nlay_s + 1 |
---|
431 | |
---|
432 | ! surface equation |
---|
433 | ztrid(numeqmin,1) = 0.0 |
---|
434 | ztrid(numeqmin,2) = dzf - zkappa_i(0)*zg1 |
---|
435 | ztrid(numeqmin,3) = zkappa_i(0)*zg1 |
---|
436 | zindterm(numeqmin) = dzf*t_su_b(ji) - zf |
---|
437 | |
---|
438 | ! first layer of ice equation |
---|
439 | ztrid(numeqmin+1,1) = - zkappa_i(0)*zg1*zeta_i(1) |
---|
440 | ztrid(numeqmin+1,2) = 1.0 + zeta_i(1)*(zkappa_i(1) + zkappa_i(0)* |
---|
441 | & zg1) |
---|
442 | ztrid(numeqmin+1,3) = - zeta_i(1)*zkappa_i(1) |
---|
443 | zindterm(numeqmin+1)= ztiold(1) + zeta_i(1)*radab_i(1) ! ( radab_phy_i(1) + |
---|
444 | ! & radab_alg_i(1) ) |
---|
445 | |
---|
446 | ! case of only one layer in the ice (surface & ice equations are altered) |
---|
447 | if (nlay_i.eq.1) then |
---|
448 | ztrid(numeqmin,1) = 0.0 |
---|
449 | ztrid(numeqmin,2) = dzf - zkappa_i(0)*2.0 |
---|
450 | ztrid(numeqmin,3) = zkappa_i(0)*2.0 |
---|
451 | |
---|
452 | ztrid(numeqmin+1,1) = -zkappa_i(0)*2.0*zeta_i(1) |
---|
453 | ztrid(numeqmin+1,2) = 1.0 + zeta_i(1)*(zkappa_i(0)*2.0 + |
---|
454 | & zkappa_i(1)) |
---|
455 | ztrid(numeqmin+1,3) = 0.0 |
---|
456 | |
---|
457 | zindterm(numeqmin+1) = ztiold(1) + zeta_i(1)* |
---|
458 | ! & ( radab_phy_i(1) + radab_alg_i(1) + |
---|
459 | & ( radab_i(1) + |
---|
460 | & zkappa_i(1)*t_bo_b(ji) ) |
---|
461 | endif |
---|
462 | |
---|
463 | else |
---|
464 | ! |
---|
465 | !------------------------------------------------------------------------------| |
---|
466 | ! case 4 : surface is melting - no snow | |
---|
467 | !------------------------------------------------------------------------------| |
---|
468 | ! |
---|
469 | zdifcase = 4.0 |
---|
470 | numeqmin = nlay_s + 2 |
---|
471 | numeqmax = nlay_i + nlay_s + 1 |
---|
472 | |
---|
473 | ! first layer of ice equation |
---|
474 | ztrid(numeqmin,1) = 0.0 |
---|
475 | ztrid(numeqmin,2) = 1.0 + zeta_i(1)*(zkappa_i(1) + zkappa_i(0)* |
---|
476 | & zg1) |
---|
477 | ztrid(numeqmin,3) = - zeta_i(1)* zkappa_i(1) |
---|
478 | zindterm(numeqmin) = ztiold(1) + zeta_i(1)* ( radab_i(1) + !(radab_phy_i(1) + |
---|
479 | ! & radab_alg_i(1) + |
---|
480 | & zkappa_i(0)*zg1*t_su_b(ji) ) |
---|
481 | |
---|
482 | ! case of only one layer in the ice (surface & ice equations are altered) |
---|
483 | if (nlay_i.eq.1) then |
---|
484 | ztrid(numeqmin,1) = 0.0 |
---|
485 | ztrid(numeqmin,2) = 1.0 + zeta_i(1)*(zkappa_i(0)*2.0 + |
---|
486 | & zkappa_i(1)) |
---|
487 | ztrid(numeqmin,3) = 0.0 |
---|
488 | zindterm(numeqmin) = ztiold(1) + zeta_i(1)* |
---|
489 | ! & (radab_phy_i(1) + radab_alg_i(1) |
---|
490 | & ( radab_i(1) + |
---|
491 | & + zkappa_i(1)*t_bo_b(ji) ) |
---|
492 | & + t_su_b(ji)*zeta_i(1)*zkappa_i(0)*2.0 |
---|
493 | endif |
---|
494 | |
---|
495 | endif |
---|
496 | endif |
---|
497 | ! |
---|
498 | !------------------------------------------------------------------------------| |
---|
499 | ! 9) tridiagonal system solving | |
---|
500 | !------------------------------------------------------------------------------| |
---|
501 | ! |
---|
502 | ! solving the tridiagonal system with Gauss elimination |
---|
503 | ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, |
---|
504 | ! McGraw-Hill 1984. |
---|
505 | ! computing temporary results |
---|
506 | zindtbis(numeqmin) = zindterm(numeqmin) |
---|
507 | zdiagbis(numeqmin) = ztrid(numeqmin,2) |
---|
508 | |
---|
509 | do numeq = numeqmin+1, numeqmax |
---|
510 | zdiagbis(numeq) = ztrid(numeq,2) - ztrid(numeq,1)* |
---|
511 | & ztrid(numeq-1,3)/zdiagbis(numeq-1) |
---|
512 | zindtbis(numeq) = zindterm(numeq) - ztrid(numeq,1)* |
---|
513 | & zindtbis(numeq-1)/zdiagbis(numeq-1) |
---|
514 | end do |
---|
515 | |
---|
516 | ! ice temperatures |
---|
517 | t_i_b(ji,nlay_i) = zindtbis(numeqmax)/zdiagbis(numeqmax) |
---|
518 | ! do numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1 |
---|
519 | do numeq = nlay_i + nlay_s, nlay_s + 2, -1 |
---|
520 | layer = numeq - nlay_s - 1 |
---|
521 | t_i_b(ji,layer) = (zindtbis(numeq) - ztrid(numeq,3)* |
---|
522 | & t_i_b(ji,layer+1))/zdiagbis(numeq) |
---|
523 | end do |
---|
524 | |
---|
525 | ! snow temperatures |
---|
526 | if (ht_s_b(ji).gt.0.0) then |
---|
527 | t_s_b(ji,nlay_s) = (zindtbis(nlay_s+1) - ztrid(nlay_s+1,3)* |
---|
528 | & t_i_b(ji,1))/zdiagbis(nlay_s+1) |
---|
529 | if (nlay_s.gt.1) then |
---|
530 | do numeq = nlay_s, 2, -1 |
---|
531 | layer = numeq - 1 |
---|
532 | t_s_b(ji,layer) = (zindtbis(numeq) - ztrid(numeq,3)* |
---|
533 | & t_s_b(ji,layer+1))/zdiagbis(numeq) |
---|
534 | end do |
---|
535 | endif |
---|
536 | endif |
---|
537 | |
---|
538 | ! surface temperature |
---|
539 | if (t_su_b(ji).lt.tfs(ji)) then |
---|
540 | t_su_b(ji) = ( zindtbis(numeqmin) - ztrid(numeqmin,3)* |
---|
541 | & ( real(isnow)*t_s_b(ji,1) + |
---|
542 | & (1.0-real(isnow))*t_i_b(ji,1) ) ) / |
---|
543 | & zdiagbis(numeqmin) |
---|
544 | endif |
---|
545 | ! |
---|
546 | !-------------------------------------------------------------------------- |
---|
547 | ! 10) Has the scheme converged ?, end of the iterative procedure | |
---|
548 | !-------------------------------------------------------------------------- |
---|
549 | ! |
---|
550 | ! we verify that nothing has started to melt |
---|
551 | t_su_b(ji) = MIN( t_su_b(ji) , tfs(ji) ) |
---|
552 | DO layer = 1, nlay_i |
---|
553 | ztmelt_i = MIN( -tmut*s_i_b(ji,layer) + tpw , |
---|
554 | & 273.149999999d0) |
---|
555 | t_i_b(ji,layer) = MIN( t_i_b(ji,layer) ,ztmelt_i ) |
---|
556 | END DO |
---|
557 | |
---|
558 | ! zerrit is a residual which has to be under zerrmax |
---|
559 | zerrit = ABS(t_su_b(ji)-ztsutemp) |
---|
560 | DO layer = 1, nlay_s |
---|
561 | zerrit = MAX(zerrit,ABS(t_s_b(ji,layer) |
---|
562 | & - ztstemp(layer))) |
---|
563 | END DO |
---|
564 | DO layer = 1, nlay_i |
---|
565 | zerrit = max(zerrit,abs(t_i_b(ji,layer) - ztitemp(layer))) |
---|
566 | END DO |
---|
567 | |
---|
568 | 42 CONTINUE |
---|
569 | |
---|
570 | END DO |
---|
571 | |
---|
572 | DO 45 ji = kideb, kiut |
---|
573 | ! compute available energy for internal melt |
---|
574 | f_s_im(ji) = 0. |
---|
575 | DO layer = 1, nlay_s |
---|
576 | IF ( t_s_b(ji,layer) .GE. tfs(ji) ) THEN |
---|
577 | f_s_im(ji) = - rhon * ( cpg*( tfs(ji) - t_s_b(ji,layer))) |
---|
578 | & * ht_s_b(ji) / ddtb |
---|
579 | t_s_b(ji,layer) = MIN( t_s_b(ji,layer) ,tfs(ji) ) |
---|
580 | ENDIF |
---|
581 | END DO |
---|
582 | |
---|
583 | 45 CONTINUE |
---|
584 | ! |
---|
585 | !-------------------------------------------------------------------------- |
---|
586 | ! 11) Heat conduction fluxes | |
---|
587 | !-------------------------------------------------------------------------- |
---|
588 | ! |
---|
589 | |
---|
590 | DO 50 ji = kideb, kiut |
---|
591 | ! surface conduction flux |
---|
592 | fc_su(ji) = - real(isnow)*zkappa_s(0)*zg1s*(t_s_b(ji,1) - |
---|
593 | & t_su_b(ji)) |
---|
594 | & - (1.0-real(isnow))*zkappa_i(0)*zg1* |
---|
595 | & (t_i_b(ji,1) - t_su_b(ji)) |
---|
596 | |
---|
597 | ! bottom conduction flux |
---|
598 | fc_bo_i(ji) = - zkappa_i(nlay_i)* |
---|
599 | & ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) |
---|
600 | |
---|
601 | ! internal conduction fluxes : snow |
---|
602 | !--upper snow value |
---|
603 | fc_s(ji,0) = - isnow* |
---|
604 | & zkappa_s(0) * zg1s * ( t_s_b(ji,1) - |
---|
605 | & t_su_b(ji) ) |
---|
606 | !--basal snow value |
---|
607 | fc_s(ji,1) = - isnow* |
---|
608 | & zkappa_s(1) * ( t_i_b(ji,1) - |
---|
609 | & t_s_b(ji,1) ) |
---|
610 | |
---|
611 | ! internal conduction fluxes : ice |
---|
612 | !--upper layer |
---|
613 | fc_i(ji,0) = - isnow * ! interface flux if there is snow |
---|
614 | & ( zkappa_i(0) * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) |
---|
615 | & - ( 1.0 - isnow ) * ( zkappa_i(0) * |
---|
616 | & zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if no |
---|
617 | !--internal ice layers |
---|
618 | DO layer = 1, nlay_i - 1 |
---|
619 | fc_i(ji,layer) = - zkappa_i(layer) * ( t_i_b(ji,layer+1) - |
---|
620 | & t_i_b(ji,layer) ) |
---|
621 | END DO |
---|
622 | !--under the basal ice layer |
---|
623 | fc_i(ji,nlay_i) = fc_bo_i(ji) |
---|
624 | |
---|
625 | ! case of only one layer in the ice |
---|
626 | IF (nlay_i.EQ.1) THEN |
---|
627 | fc_su(ji) = -real(isnow)*(zkappa_s(0)*(zg1s*(t_s_b(ji,1)- |
---|
628 | & t_su_b(ji)))) |
---|
629 | & -(1.0-real(isnow))*zkappa_i(0)*zg1*(t_i_b(ji,1)- |
---|
630 | & t_su_b(ji)) |
---|
631 | fc_bo_i(ji) = -zkappa_i(nlay_i)*zg1* |
---|
632 | & (t_bo_b(ji)-t_i_b(ji,nlay_i)) |
---|
633 | ENDIF |
---|
634 | ! |
---|
635 | !-------------------------------------------------------------------------- |
---|
636 | ! 12) Update atmospheric heat fluxes and energy of melting | |
---|
637 | !-------------------------------------------------------------------------- |
---|
638 | ! |
---|
639 | ! pressure of water vapor saturation (Pa) |
---|
640 | es = 611.0*10.0**(9.5*(t_su_b(ji)-273.16)/ |
---|
641 | & (t_su_b(ji)-7.66)) |
---|
642 | ! net longwave radiative flux |
---|
643 | ! MV BUG |
---|
644 | ! fratsb(ji) = emig*(ratbqb(ji)- |
---|
645 | ! & stefan*t_su_b(ji)*t_su_b(ji)*t_su_b(ji)*t_su_b(ji)) |
---|
646 | fratsb(ji) = ratbqb(ji) - emig * |
---|
647 | & stefan*t_su_b(ji)*t_su_b(ji)*t_su_b(ji)*t_su_b(ji) |
---|
648 | ! sensible and latent heat flux |
---|
649 | CALL flx(ht_i_b(ji),t_su_b(ji),tabqb(ji),qabqb(ji),fcsb(ji), |
---|
650 | & fleb(ji),qsfcb(ji),zrchu1(ji),zrchu2(ji),vabqb(ji),zref) |
---|
651 | |
---|
652 | fcsb(ji) = -fcsb(ji) |
---|
653 | fleb(ji) = MIN( -fleb(ji) , 0.0 ) ! always negative, as precip |
---|
654 | ! energy already added |
---|
655 | |
---|
656 | ! ice energy of melting |
---|
657 | CALL ice_th_enmelt(kideb, kiut, nlay_s, nlay_i) |
---|
658 | |
---|
659 | IF ( ln_write ) THEN |
---|
660 | WRITE(numout,*) ' nconv : ', nconv |
---|
661 | WRITE(numout,*) ' zerrit : ', zerrit |
---|
662 | WRITE(numout,*) ' t_su_b: ', t_su_b(ji) |
---|
663 | WRITE(numout,*) ' t_s_b : ', ( t_s_b(ji,layer), |
---|
664 | & layer = 1, nlay_s ) |
---|
665 | WRITE(numout,*) ' t_i_b : ', ( t_i_b(ji,layer), |
---|
666 | & layer = 1, nlay_i ) |
---|
667 | WRITE(numout,*) ' t_bo_b : ', t_bo_b(ji) |
---|
668 | WRITE(numout,*) |
---|
669 | ENDIF |
---|
670 | |
---|
671 | 50 CONTINUE |
---|
672 | ! |
---|
673 | !------------------------------------------------------------------------------ |
---|
674 | ! End of ice_th_diff |
---|
675 | END SUBROUTINE |
---|