New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

limthd_dif.F90 in branches/dev_002_LIM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/dev_002_LIM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 830

Last change on this file since 830 was 825, checked in by ctlod, 16 years ago
dev_002_LIM : add the LIM 3.0 component, see ticketr: #71
File size: 39.0 KB

Line
1	MODULE limthd_dif
2	#if defined key_lim3
3	!!======================================================================
4	!! * MODULE limthd_dif *
5	!! heat diffusion in sea ice
6	!! computation of surface and inner T
7	!!======================================================================
8
9	!!----------------------------------------------------------------------
10	!! * Modules used
11	USE par_oce ! ocean parameters
12	USE phycst ! physical constants (ocean directory)
13	USE ice_oce ! ice variables
14	USE thd_ice
15	USE iceini
16	USE limicepoints
17	USE limistate
18	USE in_out_manager
19	USE ice
20	USE par_ice
21
22	IMPLICIT NONE
23	PRIVATE
24
25	!! * Routine accessibility
26	PUBLIC lim_thd_dif ! called by lim_thd
27
28	!! * Module variables
29	REAL(wp) :: & ! constant values
30	epsi20 = 1e-20 , &
31	epsi13 = 1e-13 , &
32	zzero = 0.e0 , &
33	zone = 1.e0
34
35	!!----------------------------------------------------------------------
36	!! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2005)
37	!! (c) UCL-ASTR and Martin Vancoppenolle
38	!!----------------------------------------------------------------------
39
40	CONTAINS
41
42	SUBROUTINE lim_thd_dif( kideb , kiut , jl )
43	!!------------------------------------------------------------------
44	!! * ROUTINE lim_thd_dif *
45	!! ** Purpose :
46	!! This routine determines the time evolution of snow and sea-ice
47	!! temperature profiles.
48	!! ** Method :
49	!! This is done by solving the heat equation diffusion with
50	!! a Neumann boundary condition at the surface and a Dirichlet one
51	!! at the bottom. Solar radiation is partially absorbed into the ice.
52	!! The specific heat and thermal conductivities depend on ice salinity
53	!! and temperature to take into account brine pocket melting. The
54	!! numerical
55	!! scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid
56	!! in the ice and snow system.
57	!!
58	!! The successive steps of this routine are
59	!! 1. Thermal conductivity at the interfaces of the ice layers
60	!! 2. Internal absorbed radiation
61	!! 3. Scale factors due to non-uniform grid
62	!! 4. Kappa factors
63	!! Then iterative procedure begins
64	!! 5. specific heat in the ice
65	!! 6. eta factors
66	!! 7. surface flux computation
67	!! 8. tridiagonal system terms
68	!! 9. solving the tridiagonal system with Gauss elimination
69	!! Iterative procedure ends according to a criterion on evolution
70	!! of temperature
71	!!
72	!! ** Arguments :
73	!! kideb , kiut : Starting and ending points on which the
74	!! the computation is applied
75	!!
76	!! ** Inputs / Ouputs : (global commons)
77	!! surface temperature : t_su_b
78	!! ice/snow temperatures : t_i_b, t_s_b
79	!! ice salinities : s_i_b
80	!! number of layers in the ice/snow: nlay_i, nlay_s
81	!! profile of the ice/snow layers : z_i, z_s
82	!! total ice/snow thickness : ht_i_b, ht_s_b
83	!!
84	!! ** External :
85	!!
86	!! ** References :
87	!!
88	!! ** History :
89	!! (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium
90	!! (06-2005) Martin Vancoppenolle still!!! for the 3d version...
91	!! (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR)
92	!! (04-2007) Energy conservation tested by M. Vancoppenolle
93	!!
94	!!------------------------------------------------------------------
95	!! * Arguments
96
97	INTEGER , INTENT (in) :: &
98	kideb , & ! Start point on which the the computation is applied
99	kiut , & ! End point on which the the computation is applied
100	jl ! Category number
101
102	!! * Local variables
103	INTEGER :: ji, & ! spatial loop index
104	zji, zjj, & ! temporary dummy loop index
105	numeq, & ! current reference number of equation
106	layer, & ! vertical dummy loop index
107	nconv, & ! number of iterations in iterative procedure
108	! nconvmax = 50, &! maximum number of iterations for temperature computations (50)
109	index, & !
110	minnumeqmin, & !
111	maxnumeqmax
112
113	INTEGER , DIMENSION(jpij) :: &
114	numeqmin, & ! reference number of top equation
115	numeqmax, & ! reference number of bottom equation
116	isnow ! switch for presence (1) or absence (0) of snow
117
118	!! * New local variables
119	REAL(wp) , DIMENSION(jpij,0:nlay_i) :: &
120	ztcond_i, & !Ice thermal conductivity
121	zradtr_i, & !Radiation transmitted through the ice
122	zradab_i, & !Radiation absorbed in the ice
123	zkappa_i !Kappa factor in the ice
124
125	REAL(wp) , DIMENSION(jpij,0:nlay_s) :: &
126	zradtr_s, & !Radiation transmited through the snow
127	zradab_s, & !Radiation absorbed in the snow
128	zkappa_s !Kappa factor in the snow
129
130	REAL(wp) , DIMENSION(jpij,0:nlay_i) :: &
131	ztiold, & !Old temperature in the ice
132	zeta_i, & !Eta factor in the ice
133	ztitemp, & !Temporary temperature in the ice to check the convergence
134	zspeche_i, & !Ice specific heat
135	z_i !Vertical cotes of the layers in the ice
136
137	REAL(wp) , DIMENSION(jpij,0:nlay_s) :: &
138	zeta_s, & !Eta factor in the snow
139	ztstemp, & !Temporary temperature in the snow to check the convergence
140	ztsold, & !Temporary temperature in the snow
141	z_s !Vertical cotes of the layers in the snow
142
143	REAL(wp) , DIMENSION(jpij,jkmax+2) :: &
144	zindterm, & ! Independent term
145	zindtbis, & ! temporary independent term
146	zdiagbis
147
148	REAL(wp) , DIMENSION(jpij,jkmax+2,3) :: &
149	ztrid ! tridiagonal system terms
150
151	REAL(wp), DIMENSION(jpij) :: &
152	ztfs , & ! ice melting point
153	ztsuold , & ! old surface temperature (before the iterative
154	! procedure )
155	ztsuoldit, & ! surface temperature at previous iteration
156	zh_i , & !ice layer thickness
157	zh_s , & !snow layer thickness
158	zfsw , & !solar radiation absorbed at the surface
159	zf , & ! surface flux function
160	dzf , & ! derivative of the surface flux function
161	zksn !effective snow conductivity
162
163	REAL(wp) :: & ! constant values
164	zeps = 1.0e-10, & !
165	zg1s = 2.0, & !: for the tridiagonal system
166	zg1 = 2.0, &
167	zgamma = 18009.0, & !: for specific heat
168	zbeta = 0.117, & !: for thermal conductivity (could be 0.13)
169	zraext_s = 1.0e08, & !: extinction coefficient of radiation in the snow
170	! zraext_i = 1.5, & !extinction coefficient of radiation in the ice
171	! zraext_i = 1.0, & !extinction coefficient of radiation in the ice
172	! MY 0.87, FY 1.2 (Grenfell, Perovich and Light
173	! 07)
174	! zerrmax = 1.0e-4 , & !: error criterion for convergence of the iterative process
175	zfc_int , & !: conductive heat flux at the interface
176	zkimin = 0.10 , & !: minimum ice thermal conductivity
177	zht_smin = 1.0e-4 !: minimum snow depth
178
179	REAL(wp) :: & ! local variables
180	zheshth, & ! = zhe(ji) / thth
181	ztmelt_i, & ! ice melting temperature
182	zerritmax, & ! current maximal error on temperature
183	zexp
184
185	REAL(wp), DIMENSION(jpij) :: &
186	zerrit, & ! current error on temperature
187	zdifcase, & ! case of the equation resolution (1->4)
188	zghe, & ! correction factor of the thermal conductivity
189	zhe, & ! effective thickness for compu. of equ. thermal conductivity
190	zftrice, & ! solar radiation transmitted through the ice
191	zihic, zihe, zhsu
192
193	CHARACTER (len=50) :: charout
194	!!-- End of declarations
195	!!----------------------------------------------------------------------------------------------
196
197	! new namelist parameters
198	! maxer_i_thd = maxer_i_thd
199	! thcon_i_swi = thcon_i_swi
200	! kappa_i = kappa_i
201	! nconv_i_thd = nconv_i_thd
202
203	IF(lwp) WRITE(numout,*)'lim_thd_dif : Heat diffusion in sea ice for cat :', jl
204
205	thcon_i_swi = 1 ! Pringle et al. (2007)
206	!
207	!------------------------------------------------------------------------------!
208	! 1) Initialization !
209	!------------------------------------------------------------------------------!
210	!
211	DO ji = kideb , kiut
212	! is there snow or not
213	isnow(ji)= INT ( 1.0 - MAX( 0.0 , SIGN (1.0, - ht_s_b(ji) ) ) )
214	! surface temperature of fusion
215	ztfs(ji) = isnow(ji) * rtt + (1.0-isnow(ji)) * rtt
216	! layer thickness
217	zh_i(ji) = ht_i_b(ji) / nlay_i
218	zh_s(ji) = ht_s_b(ji) / nlay_s
219	END DO
220
221	!--------------------
222	! Ice / snow layers
223	!--------------------
224
225	z_s(:,0) = 0.0 ! vert. coord. of the up. lim. of the 1st snow layer
226	z_i(:,0) = 0.0 ! vert. coord. of the up. lim. of the 1st ice layer
227
228	DO layer = 1, nlay_s
229	DO ji = kideb , kiut
230	! vert. coord of the up. lim. of the layer-th snow layer
231	z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / nlay_s
232	END DO
233	END DO
234
235	DO layer = 1, nlay_i
236	DO ji = kideb , kiut
237	! vert. coord of the up. lim. of the layer-th ice layer
238	z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / nlay_i
239	END DO
240	END DO
241	!
242	!------------------------------------------------------------------------------\|
243	! 2) Radiations \|
244	!------------------------------------------------------------------------------\|
245	!
246	!-------------------
247	! Computation of i0
248	!-------------------
249	! i0 describes the fraction of solar radiation which does not contribute
250	! to the surface energy budget but rather penetrates inside the ice.
251	! We assume that no radiation is transmitted through the snow
252	! If there is no no snow
253	! zfsw = (1-i0).qsr_ice is absorbed at the surface
254	! zftrice = io.qsr_ice is below the surface
255	! fstbif = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
256
257	DO ji = kideb , kiut
258	! switches
259	isnow(ji) = INT ( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) )
260	! hs > 0, isnow = 1
261	zhsu(ji) = hnzst !threshold for the computation of i0
262	zihic(ji) = MAX( zzero , 1.0 - ( ht_i_b(ji) / zhsu(ji) ) )
263
264	i0(ji) = ( 1.0 - isnow(ji) ) * &
265	( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
266	!fr1_i0_1d = i0 for a thin ice surface
267	!fr1_i0_2d = i0 for a thick ice surface
268	! a function of the cloud cover
269	!
270	!i0(ji) = (1.0-FLOAT(isnow(ji)))3.0/(100ht_s_b(ji)+10.0)
271	!formula used in Cice
272	END DO
273
274	!-------------------------------------------------------
275	! Solar radiation absorbed / transmitted at the surface
276	! Derivative of the non solar flux
277	!-------------------------------------------------------
278	DO ji = kideb , kiut
279
280	! Shortwave radiation absorbed at surface
281	zfsw(ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) )
282
283	! Solar radiation transmitted below the surface layer
284	zftrice(ji)= qsr_ice_1d(ji) * i0(ji)
285
286	! derivative of incoming nonsolar flux
287	dzf(ji) = dqns_ice_1d(ji)
288
289	END DO
290
291	!---------------------------------------------------------
292	! Transmission - absorption of solar radiation in the ice
293	!---------------------------------------------------------
294
295	DO ji = kideb , kiut
296	! Initialization
297	zradtr_s(ji,0) = zftrice(ji) ! radiation penetrating through snow
298	END DO
299
300	! Radiation through snow
301	DO layer = 1, nlay_s
302	DO ji = kideb , kiut
303	! radiation transmitted below the layer-th snow layer
304	zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP ( - zraext_s * ( MAX ( 0.0 , &
305	z_s(ji,layer) ) ) )
306	! radiation absorbed by the layer-th snow layer
307	zradab_s(ji,layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)
308	END DO
309	END DO
310
311	! Radiation through ice
312	DO ji = kideb , kiut
313	zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + &
314	zftrice(ji) * ( 1 - isnow(ji) )
315	END DO
316
317	DO layer = 1, nlay_i
318	DO ji = kideb , kiut
319	! radiation transmitted below the layer-th ice layer
320	zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP ( - kappa_i * ( MAX ( 0.0 , &
321	z_i(ji,layer) ) ) )
322	! radiation absorbed by the layer-th ice layer
323	zradab_i(ji,layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)
324	END DO
325	END DO
326
327	! Radiation transmitted below the ice
328	DO ji = kideb , kiut
329	fstbif_1d(ji) = fstbif_1d(ji) + &
330	zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji)
331	END DO
332
333	! +++++
334	! just to check energy conservation
335	DO ji = kideb , kiut
336	zji = MOD( npb(ji) - 1, jpi ) + 1
337	zjj = ( npb(ji) - 1 ) / jpi + 1
338	fstroc(zji,zjj,jl) = &
339	zradtr_i(ji,nlay_i)
340	END DO
341	DO layer = 1, nlay_i
342	DO ji = kideb , kiut
343	radab(ji,layer) = zradab_i(ji,layer)
344	END DO
345	END DO
346
347	! +++++
348
349	!
350	!------------------------------------------------------------------------------\|
351	! 3) Iterative procedure begins \|
352	!------------------------------------------------------------------------------\|
353	!
354	! Old surface temperature
355	DO ji = kideb, kiut
356	ztsuold(ji) = t_su_b(ji) ! temperature at the beg of iter pr.
357	ztsuoldit(ji) = t_su_b(ji) ! temperature at the previous iter
358	t_su_b(ji) = MIN(t_su_b(ji),ztfs(ji)-0.00001) !necessary
359	zerrit(ji) = 1000.0 ! initial value of error
360	END DO
361
362	! Old snow temperature
363	DO layer = 1, nlay_s
364	DO ji = kideb , kiut
365	ztsold(ji,layer) = t_s_b(ji,layer)
366	END DO
367	END DO
368
369	! Old ice temperature
370	DO layer = 1, nlay_i
371	DO ji = kideb , kiut
372	ztiold(ji,layer) = t_i_b(ji,layer)
373	END DO
374	END DO
375
376	nconv = 0 ! number of iterations
377	zerritmax = 1000.0 ! maximal value of error on all points
378
379	DO WHILE ((zerritmax > maxer_i_thd).AND.(nconv < nconv_i_thd))
380
381	! ! +++++
382	! DO ji = kideb , kiut
383	! zji = MOD( npb(ji) - 1, jpi ) + 1
384	! zjj = ( npb(ji) - 1 ) / jpi + 1
385	! IF ( (zji.EQ.jiindex) .AND. (zjj.EQ.jjindex) ) THEN
386	! WRITE(numout,*) ' Iteration : ', nconv
387	! WRITE(numout,*) ' t_su : ', t_su_b(ji)
388	! WRITE(numout,*) ' t_s : ', t_s_b(ji,1)
389	! WRITE(numout,*) ' t_i : ', t_i_b(ji,:)
390	! WRITE(numout,*) ' difcase : ', zdifcase(ji)
391	! ENDIF
392	! END DO
393	! ! +++++
394
395	nconv = nconv+1
396
397	!
398	!------------------------------------------------------------------------------\|
399	! 4) Sea ice thermal conductivity \|
400	!------------------------------------------------------------------------------\|
401	!
402	IF ( thcon_i_swi .EQ. 0 ) THEN
403	! Untersteiner (1964) formula
404	DO ji = kideb , kiut
405	ztcond_i(ji,0) = rcdic + zbeta*s_i_b(ji,1) / &
406	MIN(-zeps,t_i_b(ji,1)-rtt)
407	ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin)
408	END DO
409	ENDIF
410
411	IF ( thcon_i_swi .EQ. 1 ) THEN
412	! Pringle et al formula included, should be tested
413	! 2.11 + 0.09 S/T - 0.011.T
414	DO ji = kideb , kiut
415	ztcond_i(ji,0) = rcdic + 0.09*s_i_b(ji,1) / &
416	MIN(-zeps,t_i_b(ji,1)-rtt) - &
417	0.011* ( t_i_b(ji,1) - rtt )
418	ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin)
419	END DO
420	ENDIF
421
422	IF ( thcon_i_swi .EQ. 0 ) THEN ! Untersteiner
423	DO layer = 1, nlay_i-1
424	DO ji = kideb , kiut
425	ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) &
426	+ s_i_b(ji,layer+1) ) / MIN(-zeps, &
427	t_i_b(ji,layer)+t_i_b(ji,layer+1)-2.0*rtt)
428	ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin)
429	END DO
430	END DO
431	ENDIF
432
433	IF ( thcon_i_swi .EQ. 1 ) THEN ! Pringle
434	DO layer = 1, nlay_i-1
435	DO ji = kideb , kiut
436	ztcond_i(ji,layer) = rcdic + 0.09*( s_i_b(ji,layer) &
437	+ s_i_b(ji,layer+1) ) / MIN(-zeps, &
438	t_i_b(ji,layer)+t_i_b(ji,layer+1)-2.0*rtt) - &
439	0.011* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )
440	ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin)
441	END DO
442	END DO
443	ENDIF
444
445	IF ( thcon_i_swi .EQ. 0 ) THEN ! Untersteiner
446	DO ji = kideb , kiut
447	ztcond_i(ji,nlay_i) = rcdic + zbeta*s_i_b(ji,nlay_i) / &
448	MIN(-zeps,t_bo_b(ji)-rtt)
449	ztcond_i(ji,nlay_i) = MAX(ztcond_i(ji,nlay_i),zkimin)
450	END DO
451	ENDIF
452
453	IF ( thcon_i_swi .EQ. 1 ) THEN ! Pringle
454	DO ji = kideb , kiut
455	ztcond_i(ji,nlay_i) = rcdic + 0.09*s_i_b(ji,nlay_i) / &
456	MIN(-zeps,t_bo_b(ji)-rtt) - &
457	0.011* ( t_bo_b(ji) - rtt )
458	ztcond_i(ji,nlay_i) = MAX(ztcond_i(ji,nlay_i),zkimin)
459	END DO
460	ENDIF
461	!
462	!------------------------------------------------------------------------------\|
463	! 5) kappa factors \|
464	!------------------------------------------------------------------------------\|
465	!
466	DO ji = kideb, kiut
467
468	!-- Snow kappa factors
469	zkappa_s(ji,0) = rcdsn / MAX(zeps,zh_s(ji))
470	zkappa_s(ji,nlay_s) = rcdsn / MAX(zeps,zh_s(ji))
471	END DO
472
473	DO layer = 1, nlay_s-1
474	DO ji = kideb , kiut
475	zkappa_s(ji,layer) = 2.0 * rcdsn / &
476	MAX(zeps,2.0*zh_s(ji))
477	END DO
478	END DO
479
480	DO layer = 1, nlay_i-1
481	DO ji = kideb , kiut
482	!-- Ice kappa factors
483	zkappa_i(ji,layer) = 2.0*ztcond_i(ji,layer)/ &
484	MAX(zeps,2.0*zh_i(ji))
485	END DO
486	END DO
487
488	DO ji = kideb , kiut
489	zkappa_i(ji,0) = ztcond_i(ji,0)/MAX(zeps,zh_i(ji))
490	zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(zeps,zh_i(ji))
491
492	!-- Interface
493	zkappa_s(ji,nlay_s) = 2.0rcdsnztcond_i(ji,0)/MAX(zeps, &
494	(ztcond_i(ji,0)zh_s(ji) + rcdsnzh_i(ji)))
495	zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*isnow(ji) &
496	+ zkappa_i(ji,0)*(1.0-isnow(ji))
497	END DO
498	!
499	!------------------------------------------------------------------------------\|
500	! 6) Sea ice specific heat, eta factors \|
501	!------------------------------------------------------------------------------\|
502	!
503	DO layer = 1, nlay_i
504	DO ji = kideb , kiut
505	ztitemp(ji,layer) = t_i_b(ji,layer)
506	zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ &
507	MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),zeps)
508	zeta_i(ji,layer) = rdt_ice / MAX(rhoiczspeche_i(ji,layer)zh_i(ji), &
509	zeps)
510	END DO
511	END DO
512
513	DO layer = 1, nlay_s
514	DO ji = kideb , kiut
515	ztstemp(ji,layer) = t_s_b(ji,layer)
516	zeta_s(ji,layer) = rdt_ice / MAX(rhosncpiczh_s(ji),zeps)
517	END DO
518	END DO
519	!
520	!------------------------------------------------------------------------------\|
521	! 7) surface flux computation \|
522	!------------------------------------------------------------------------------\|
523	!
524	DO ji = kideb , kiut
525
526	! update of the non solar flux according to the update in T_su
527	qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * &
528	( t_su_b(ji) - ztsuoldit(ji) )
529
530	! update incoming flux
531	zf(ji) = zfsw(ji) & ! net absorbed solar radiation
532	+ qnsr_ice_1d(ji) ! non solar total flux
533	! (LWup, LWdw, SH, LH)
534
535	END DO
536
537	!
538	!------------------------------------------------------------------------------\|
539	! 8) tridiagonal system terms \|
540	!------------------------------------------------------------------------------\|
541	!
542	!!layer denotes the number of the layer in the snow or in the ice
543	!!numeq denotes the reference number of the equation in the tridiagonal
544	!!system, terms of tridiagonal system are indexed as following :
545	!!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
546
547	!!ice interior terms (top equation has the same form as the others)
548
549	DO numeq=1,jkmax+2
550	DO ji = kideb , kiut
551	ztrid(ji,numeq,1) = 0.
552	ztrid(ji,numeq,2) = 0.
553	ztrid(ji,numeq,3) = 0.
554	zindterm(ji,numeq)= 0.
555	zindtbis(ji,numeq)= 0.
556	zdiagbis(ji,numeq)= 0.
557	ENDDO
558	ENDDO
559
560	DO numeq = nlay_s + 2, nlay_s + nlay_i
561	DO ji = kideb , kiut
562	layer = numeq - nlay_s - 1
563	ztrid(ji,numeq,1) = - zeta_i(ji,layer)*zkappa_i(ji,layer-1)
564	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,layer)*(zkappa_i(ji,layer-1) + &
565	zkappa_i(ji,layer))
566	ztrid(ji,numeq,3) = - zeta_i(ji,layer)*zkappa_i(ji,layer)
567	zindterm(ji,numeq) = ztiold(ji,layer) + zeta_i(ji,layer)* &
568	zradab_i(ji,layer)
569	END DO
570	ENDDO
571
572	numeq = nlay_s + nlay_i + 1
573	DO ji = kideb , kiut
574	!!ice bottom term
575	ztrid(ji,numeq,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)
576	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i)( zkappa_i(ji,nlay_i)zg1 &
577	+ zkappa_i(ji,nlay_i-1) )
578	ztrid(ji,numeq,3) = 0.0
579	zindterm(ji,numeq) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* &
580	( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 &
581	* t_bo_b(ji) )
582	ENDDO
583
584
585	DO ji = kideb , kiut
586	IF ( ht_s_b(ji).gt.0.0 ) THEN
587	!
588	!------------------------------------------------------------------------------\|
589	! snow-covered cells \|
590	!------------------------------------------------------------------------------\|
591	!
592	!!snow interior terms (bottom equation has the same form as the others)
593	DO numeq = 3, nlay_s + 1
594	layer = numeq - 1
595	ztrid(ji,numeq,1) = - zeta_s(ji,layer)*zkappa_s(ji,layer-1)
596	ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,layer)*( zkappa_s(ji,layer-1) + &
597	zkappa_s(ji,layer) )
598	ztrid(ji,numeq,3) = - zeta_s(ji,layer)*zkappa_s(ji,layer)
599	zindterm(ji,numeq) = ztsold(ji,layer) + zeta_s(ji,layer)* &
600	zradab_s(ji,layer)
601	END DO
602
603	!!case of only one layer in the ice (ice equation is altered)
604	IF ( nlay_i.eq.1 ) THEN
605	ztrid(ji,nlay_s+2,3) = 0.0
606	zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* &
607	t_bo_b(ji)
608	ENDIF
609
610	IF ( t_su_b(ji) .LT. rtt ) THEN
611
612	!------------------------------------------------------------------------------\|
613	! case 1 : no surface melting - snow present \|
614	!------------------------------------------------------------------------------\|
615	zdifcase(ji) = 1.0
616	numeqmin(ji) = 1
617	numeqmax(ji) = nlay_i + nlay_s + 1
618	!!surface equation
619	ztrid(ji,1,1) = 0.0
620	ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0)
621	ztrid(ji,1,3) = zg1s*zkappa_s(ji,0)
622	! bug again ???
623	! zindterm(ji,1) = dzf(ji)*ztsuold(ji) - zf(ji)
624	zindterm(ji,1) = dzf(ji)*t_su_b(ji) - zf(ji)
625
626	!!first layer of snow equation
627	ztrid(ji,2,1) = - zkappa_s(ji,0)zg1szeta_s(ji,1)
628	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)(zkappa_s(ji,1) + zkappa_s(ji,0)zg1s)
629	ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1)
630	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)
631
632	ELSE
633	!
634	!------------------------------------------------------------------------------\|
635	! case 2 : surface is melting - snow present \|
636	!------------------------------------------------------------------------------\|
637	!
638	zdifcase(ji) = 2.0
639	numeqmin(ji) = 2
640	numeqmax(ji) = nlay_i + nlay_s + 1
641
642	!!first layer of snow equation
643	ztrid(ji,2,1) = 0.0
644	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + &
645	zkappa_s(ji,0) * zg1s )
646	ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1)
647	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * &
648	( zradab_s(ji,1) + &
649	zkappa_s(ji,0) * zg1s * t_su_b(ji) )
650	ENDIF
651	ELSE
652	!
653	!------------------------------------------------------------------------------\|
654	! cells without snow \|
655	!------------------------------------------------------------------------------\|
656	!
657	IF (t_su_b(ji) .LT. rtt) THEN
658	!
659	!------------------------------------------------------------------------------\|
660	! case 3 : no surface melting - no snow \|
661	!------------------------------------------------------------------------------\|
662	!
663	zdifcase(ji) = 3.0
664	numeqmin(ji) = nlay_s + 1
665	numeqmax(ji) = nlay_i + nlay_s + 1
666
667	!!surface equation
668	ztrid(ji,numeqmin(ji),1) = 0.0
669	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1
670	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1
671	zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_b(ji) - zf(ji)
672
673	!!first layer of ice equation
674	ztrid(ji,numeqmin(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
675	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) &
676	+ zkappa_i(ji,0) * zg1 )
677	ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1)
678	zindterm(ji,numeqmin(ji)+1)= ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)
679
680	!!case of only one layer in the ice (surface & ice equations are altered)
681
682	IF (nlay_i.eq.1) THEN
683	ztrid(ji,numeqmin(ji),1) = 0.0
684	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*2.0
685	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*2.0
686	ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0)2.0zeta_i(ji,1)
687	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,0)2.0 + &
688	zkappa_i(ji,1))
689	ztrid(ji,numeqmin(ji)+1,3) = 0.0
690
691	zindterm(ji,numeqmin(ji)+1) = ztiold(ji,1) + zeta_i(ji,1)* &
692	( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji) )
693	ENDIF
694
695	ELSE
696
697	!
698	!------------------------------------------------------------------------------\|
699	! case 4 : surface is melting - no snow \|
700	!------------------------------------------------------------------------------\|
701	!
702	zdifcase(ji) = 4.0
703	numeqmin(ji) = nlay_s + 2
704	numeqmax(ji) = nlay_i + nlay_s + 1
705
706	!!first layer of ice equation
707	ztrid(ji,numeqmin(ji),1) = 0.0
708	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,1) + zkappa_i(ji,0) &
709	zg1)
710	ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)
711	zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &
712	zkappa_i(ji,0) * zg1 * t_su_b(ji) )
713
714	!!case of only one layer in the ice (surface & ice equations are altered)
715	IF (nlay_i.eq.1) THEN
716	ztrid(ji,numeqmin(ji),1) = 0.0
717	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,0)2.0 + &
718	zkappa_i(ji,1))
719	ztrid(ji,numeqmin(ji),3) = 0.0
720	zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)* &
721	(zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji)) &
722	+ t_su_b(ji)zeta_i(ji,1)zkappa_i(ji,0)*2.0
723	ENDIF
724
725	ENDIF
726	ENDIF
727
728	END DO
729
730	!
731	!------------------------------------------------------------------------------\|
732	! 9) tridiagonal system solving \|
733	!------------------------------------------------------------------------------\|
734	!
735
736	! Solve the tridiagonal system with Gauss elimination method.
737	! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
738	! McGraw-Hill 1984.
739
740	maxnumeqmax = 0
741	minnumeqmin = jkmax+4
742
743	DO ji = kideb , kiut
744	zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji))
745	zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2)
746	minnumeqmin = MIN(numeqmin(ji),minnumeqmin)
747	maxnumeqmax = MAX(numeqmax(ji),maxnumeqmax)
748	END DO
749
750	DO layer = minnumeqmin+1, maxnumeqmax
751	DO ji = kideb , kiut
752	numeq = min(max(numeqmin(ji)+1,layer),numeqmax(ji))
753	zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* &
754	ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1)
755	zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* &
756	zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1)
757	END DO
758	END DO
759
760	DO ji = kideb , kiut
761	! ice temperatures
762	t_i_b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))
763	END DO
764
765	DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1
766	DO ji = kideb , kiut
767	layer = numeq - nlay_s - 1
768	t_i_b(ji,layer) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &
769	t_i_b(ji,layer+1))/zdiagbis(ji,numeq)
770	END DO
771	END DO
772
773	DO ji = kideb , kiut
774	! snow temperatures
775	IF (ht_s_b(ji).GT.0) &
776	t_s_b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
777	* t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) &
778	* MAX(0.0,SIGN(1.0,ht_s_b(ji)-zeps))
779
780	! surface temperature
781	isnow(ji) = INT(1.0-max(0.0,sign(1.0,-ht_s_b(ji))))
782	ztsuoldit(ji) = t_su_b(ji)
783	IF (t_su_b(ji) .LT. ztfs(ji)) &
784	t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* &
785	( isnow(ji)*t_s_b(ji,1) + &
786	(1.0-isnow(ji))*t_i_b(ji,1) ) ) / &
787	zdiagbis(ji,numeqmin(ji))
788	END DO
789	!
790	!--------------------------------------------------------------------------
791	! 10) Has the scheme converged ?, end of the iterative procedure \|
792	!--------------------------------------------------------------------------
793	!
794	! check that nowhere it has started to melt
795	! zerrit(ji) is a measure of error, it has to be under maxer_i_thd
796	DO ji = kideb , kiut
797	t_su_b(ji) = MAX(MIN(t_su_b(ji),ztfs(ji)),190.0)
798	zerrit(ji) = ABS(t_su_b(ji)-ztsuoldit(ji))
799	END DO
800
801	DO layer = 1, nlay_s
802	DO ji = kideb , kiut
803	zji = MOD( npb(ji) - 1, jpi ) + 1
804	zjj = ( npb(ji) - 1 ) / jpi + 1
805	t_s_b(ji,layer) = MAX(MIN(t_s_b(ji,layer),rtt),190.0)
806	zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji,layer) &
807	- ztstemp(ji,layer)))
808	END DO
809	END DO
810
811	DO layer = 1, nlay_i
812	DO ji = kideb , kiut
813	ztmelt_i = -tmut*s_i_b(ji,layer) +rtt
814	t_i_b(ji,layer) = MAX(MIN(t_i_b(ji,layer),ztmelt_i),190.0)
815	zerrit(ji) = MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer)))
816	END DO
817	END DO
818
819	! Compute spatial maximum over all errors
820	! note that this could be optimized substantially by iterating only
821	! the non-converging points
822	zerritmax = 0.0
823	DO ji = kideb , kiut
824	zerritmax = MAX(zerritmax,zerrit(ji))
825	END DO
826
827	END DO ! End of the do while iterative procedure
828
829	WRITE(numout,*) ' zerritmax : ', zerritmax
830	WRITE(numout,*) ' nconv : ', nconv
831
832	!
833	!--------------------------------------------------------------------------
834	! 11) Fluxes at the interfaces \|
835	!--------------------------------------------------------------------------
836	!
837	DO ji = kideb, kiut
838	! update of latent heat fluxes
839	qla_ice_1d (ji) = qla_ice_1d (ji) + &
840	dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) )
841
842	! surface ice conduction flux
843	isnow(ji) = int(1.0-max(0.0,sign(1.0,-ht_s_b(ji))))
844	fc_su(ji) = - isnow(ji)zkappa_s(ji,0)zg1s*(t_s_b(ji,1) - &
845	t_su_b(ji)) &
846	- (1.0-isnow(ji))zkappa_i(ji,0)zg1* &
847	(t_i_b(ji,1) - t_su_b(ji))
848
849	! bottom ice conduction flux
850	fc_bo_i(ji) = - zkappa_i(ji,nlay_i)* &
851	( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
852
853	END DO
854
855	IF ( con_i ) THEN
856	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
857	! ++++ heat conservation ++++ '
858	! Internal fluxes ... goes to trash!!!
859
860	DO ji = kideb, kiut
861	! Upper snow value
862	fc_s(ji,0) = - isnow(ji)* &
863	zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - &
864	t_su_b(ji) )
865	! Bott. snow value
866	fc_s(ji,1) = - isnow(ji)* &
867	zkappa_s(ji,1) * ( t_i_b(ji,1) - &
868	t_s_b(ji,1) )
869	END DO
870
871	! Upper ice layer
872	DO ji = kideb, kiut
873	fc_i(ji,0) = - isnow(ji) * & ! interface flux if there is snow
874	( zkappa_i(ji,0) * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) &
875	- ( 1.0 - isnow(ji) ) * ( zkappa_i(ji,0) * &
876	zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not
877	END DO
878
879	! Internal ice layers
880	DO layer = 1, nlay_i - 1
881	DO ji = kideb, kiut
882	fc_i(ji,layer) = - zkappa_i(ji,layer) * ( t_i_b(ji,layer+1) - &
883	t_i_b(ji,layer) )
884	zji = MOD( npb(ji) - 1, jpi ) + 1
885	zjj = ( npb(ji) - 1 ) / jpi + 1
886	! IF ( ( zji.EQ.jiindex ) .AND. ( zjj.EQ.jjindex) ) THEN
887	! WRITE(numout,*) ' layer : ', layer
888	! WRITE(numout,*) ' fc_i : ', fc_i(ji,layer)
889	! WRITE(numout,*) ' zkappa_i : ', zkappa_i(ji,layer)
890	! WRITE(numout,*) ' ztcond_i : ', ztcond_i(ji,layer)
891	! WRITE(numout,*) ' zspeche_i: ', zspeche_i(ji,layer)
892	! WRITE(numout,*) ' ti2 : ', t_i_b(ji,layer+1)
893	! WRITE(numout,*) ' ti1 : ', t_i_b(ji,layer)
894	! ENDIF
895	END DO
896	END DO
897
898	! Bottom ice layers
899	DO ji = kideb, kiut
900	fc_i(ji,nlay_i) = - zkappa_i(ji,nlay_i)* &
901	( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
902	END DO
903
904	ENDIF
905
906	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
907
908	END SUBROUTINE lim_thd_dif
909
910	#else
911	!!======================================================================
912	!! * MODULE limthd_dif *
913	!! no sea ice model
914	!!======================================================================
915	CONTAINS
916	SUBROUTINE lim_thd_dif ! Empty routine
917	END SUBROUTINE lim_thd_dif
918	#endif
919	END MODULE limthd_dif

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

Download in other formats: