Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

ice_th_diff.f @ 6

Last change on this file since 6 was 6, checked in by vancop, 8 years ago
initial import of v3.20 /Users/ioulianikolskaia/Boulot/CODES/LIM1D/ARCHIVE/TMP/LIM1D_v3.20/
File size: 25.8 KB

Line
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(2maxnlay+1),zindtbis(2maxnlay+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.0ykn(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 + lfustmuts_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(rhoncpgdeltaz_s_phy(layer),zeps)
266	END DO
267
268	DO layer = 1, nlay_i
269	zeta_i(layer) = ddtb / max(rhogdeltaz_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.010.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	! & stefant_su_b(ji)t_su_b(ji)t_su_b(ji)t_su_b(ji))
288	fratsb(ji) = ratbqb(ji) - emig *
289	& stefant_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.622es)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.0emigstefan*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)zg1szeta_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)zg1zeta_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.0zeta_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)zg1t_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.010.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	! & stefant_su_b(ji)t_su_b(ji)t_su_b(ji)t_su_b(ji))
646	fratsb(ji) = ratbqb(ji) - emig *
647	& stefant_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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: tags/LIM1D_v3.20/SOURCES/source_3.20/ice_th_diff.f @ 6

Download in other formats: