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.F90 in trunk/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMO/LIM_SRC_3/limthd.F90 @ 838

Last change on this file since 838 was 834, checked in by ctlod, 16 years ago
Clean comments and useless lines, see ticket:#72
File size: 44.3 KB

Line
1	MODULE limthd
2	!!======================================================================
3	!! * MODULE limthd *
4	!! LIM thermo ice model : ice thermodynamic
5	!!======================================================================
6	#if defined key_lim3
7	!!----------------------------------------------------------------------
8	!! 'key_lim3' LIM3 sea-ice model
9	!!----------------------------------------------------------------------
10	!! lim_thd : thermodynamic of sea ice
11	!! lim_thd_init : initialisation of sea-ice thermodynamic
12	!!----------------------------------------------------------------------
13	!! * Modules used
14	USE phycst ! physical constants
15	USE dom_oce ! ocean space and time domain variables
16	USE lbclnk
17	USE in_out_manager ! I/O manager
18	USE ice ! LIM sea-ice variables
19	USE ice_oce ! sea-ice/ocean variables
20	USE flx_oce ! sea-ice/ocean forcings variables
21	USE thd_ice ! LIM thermodynamic sea-ice variables
22	USE dom_ice ! LIM sea-ice domain
23	USE iceini
24	USE limthd_dif
25	USE limthd_dh
26	USE limthd_sal
27	USE limthd_ent
28	USE limtab
29	USE par_ice
30	USE limvar
31
32	IMPLICIT NONE
33	PRIVATE
34
35	!! * Routine accessibility
36	PUBLIC lim_thd ! called by lim_step
37
38	!! * Module variables
39	REAL(wp) :: & ! constant values
40	epsi20 = 1e-20 , &
41	epsi16 = 1e-16 , &
42	epsi06 = 1e-06 , &
43	epsi04 = 1e-04 , &
44	zzero = 0.e0 , &
45	zone = 1.e0
46
47	!! * Substitutions
48	# include "domzgr_substitute.h90"
49	# include "vectopt_loop_substitute.h90"
50	!!----------------------------------------------------------------------
51	!! LIM 2.0, UCL-LOCEAN-IPSL (2005)
52	!! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limthd.F90,v 1.6 2005/03/27 18:34:42 opalod Exp $
53	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
54	!!----------------------------------------------------------------------
55
56	CONTAINS
57
58	SUBROUTINE lim_thd
59	!!-------------------------------------------------------------------
60	!! * ROUTINE lim_thd *
61	!!
62	!! ** Purpose : This routine manages the ice thermodynamic.
63	!!
64	!! ** Action : - Initialisation of some variables
65	!! - Some preliminary computation (oceanic heat flux
66	!! at the ice base, snow acc.,heat budget of the leads)
67	!! - selection of the icy points and put them in an array
68	!! - call lim_vert_ther for vert ice thermodynamic
69	!! - back to the geographic grid
70	!! - selection of points for lateral accretion
71	!! - call lim_lat_acc for the ice accretion
72	!! - back to the geographic grid
73	!!
74	!! ** References :
75	!! H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
76	!!
77	!! History :
78	!! 1.0 ! 00-01 (M.A. Morales Maqueda, H. Goosse, T. Fichefet)
79	!! 2.0 ! 02-07 (C. Ethe, G. Madec) F90
80	!! 3.0 ! 05-11 (M. Vancoppenolle ) Multi-layer thermodynamics,
81	!! salinity variations
82	!!---------------------------------------------------------------------
83	!! * Local variables
84	INTEGER :: ji, jj, jk, jl, &
85	zji , zjj, & ! dummy loop indices
86	nbpb , & ! nb of icy pts for thermo. cal.
87	index
88
89	REAL(wp) :: &
90	zfric_umin = 5e-03 , & ! lower bound for the friction velocity
91	zfric_umax = 2e-02 ! upper bound for the friction velocity
92
93	REAL(wp) :: &
94	zinda , & ! switch for test. the val. of concen.
95	zindb, & ! switches for test. the val of arg
96	zthsnice , &
97	zfric_u , & ! friction velocity
98	zfnsol , & ! total non solar heat
99	zfontn , & ! heat flux from snow thickness
100	zfntlat, zpareff , & ! test. the val. of lead heat budget
101	zeps
102
103	REAL(wp) :: &
104	zareamin
105
106	REAL(wp), DIMENSION(jpi,jpj) :: &
107	zhicifp , & ! ice thickness for outputs
108	zqlbsbq ! link with lead energy budget qldif
109
110	!!-------------------------------------------------------------------
111
112	IF( numit == nstart ) CALL lim_thd_init ! Initialization (first time-step only)
113
114	WRITE(numout,*) 'limthd : Ice Thermodynamics'
115	WRITE(numout,*) '~~~~~~'
116
117	!------------------------------------------------------------------------------!
118	! 1) Initialization of diagnostic variables !
119	!------------------------------------------------------------------------------!
120	zeps = 1.0e-10
121	tatm_ice(:,:) = tatm_ice(:,:) + 273.15 ! convert C to K
122
123	!--------------------
124	! 1.2) Heat content
125	!--------------------
126	! Change the units of heat content; from global units to
127	! J.m3
128
129	DO jl = 1, jpl
130	DO jk = 1, nlay_i
131	DO jj = 1, jpj
132	DO ji = 1, jpi
133	!Energy of melting q(S,T) [J.m-3]
134	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / area(ji,jj) / &
135	MAX( v_i(ji,jj,jl) , epsi06 ) * nlay_i
136	!0 if no ice and 1 if yes
137	zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i(ji,jj,jl) ) )
138	!convert units ! very important that this line is here
139	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb
140	END DO
141	END DO
142	END DO
143	END DO
144
145	DO jl = 1, jpl
146	DO jk = 1, nlay_s
147	DO jj = 1, jpj
148	DO ji = 1, jpi
149	!Energy of melting q(S,T) [J.m-3]
150	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / area(ji,jj) / &
151	MAX( v_s(ji,jj,jl) , epsi06 ) * nlay_s
152	!0 if no ice and 1 if yes
153	zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s(ji,jj,jl) ) )
154	!convert units ! very important that this line is here
155	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb
156	END DO
157	END DO
158	END DO
159	END DO
160
161	!-----------------------------
162	! 1.3) Set some dummies to 0
163	!-----------------------------
164	rdvosif(:,:) = 0.e0 ! variation of ice volume at surface
165	rdvobif(:,:) = 0.e0 ! variation of ice volume at bottom
166	fdvolif(:,:) = 0.e0 ! total variation of ice volume
167	rdvonif(:,:) = 0.e0 ! lateral variation of ice volume
168	fstric (:,:) = 0.e0 ! part of solar radiation transmitted through the ice
169	ffltbif(:,:) = 0.e0 ! linked with fstric
170	qfvbq (:,:) = 0.e0 ! linked with fstric
171	rdmsnif(:,:) = 0.e0 ! variation of snow mass per unit area
172	rdmicif(:,:) = 0.e0 ! variation of ice mass per unit area
173	hicifp (:,:) = 0.e0 ! daily thermodynamic ice production.
174	fsbri (:,:) = 0.e0 ! brine flux contribution to salt flux to the ocean
175	fhbri (:,:) = 0.e0 ! brine flux contribution to heat flux to the ocean
176	fseqv (:,:) = 0.e0 ! equivalent salt flux to the ocean due to ice/growth decay
177
178	!-----------------------------------
179	! 1.4) Compute global heat content
180	!-----------------------------------
181	qt_i_in(:,:) = 0.0
182	qt_s_in(:,:) = 0.0
183	qt_i_fin(:,:) = 0.0
184	qt_s_fin(:,:) = 0.0
185	sum_fluxq(:,:) = 0.0
186	fatm(:,:) = 0.0
187
188	! 2) Partial computation of forcing for the thermodynamic sea ice model. !
189	!-----------------------------------------------------------------------------!
190
191	! !CDIR NOVERRCHK
192	DO jj = 1, jpj
193	! !CDIR NOVERRCHK
194	DO ji = 1, jpi
195	zthsnice = SUM( ht_s(ji,jj,1:jpl) ) + SUM( ht_i(ji,jj,1:jpl) )
196	zindb = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) )
197	phicif(ji,jj) = vt_i(ji,jj)
198	pfrld(ji,jj) = 1.0 - at_i(ji,jj)
199	zinda = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) )
200
201	! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
202	! ! practically no "direct lateral ablation"
203	!
204	! ! net downward heat flux from the ice to the ocean, expressed as a function of ocean
205	! ! temperature and turbulent mixing (McPhee, 1992)
206	! friction velocity
207	zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )
208
209	! here the drag will depend on ice thickness and type (0.006)
210	fdtcn(ji,jj) = zindb * rau0 * rcp * 0.006 * zfric_u * ( sst_io(ji,jj) - t_bo(ji,jj) )
211	! also category dependent
212	! !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead
213	qdtcn(ji,jj) = zindb * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice
214	!
215
216	! still need to be updated : fdtcn !!!!
217	! !-- Lead heat budget (part 1, next one is in limthd_dh
218	! !-- qldif -- (or qldif_1d in 1d routines)
219	zfontn = sprecip(ji,jj) * lfus ! energy of melting
220	zfnsol = qnsr_oce(ji,jj) ! total non solar flux
221	qldif(ji,jj) = tms(ji,jj) * ( qsr_oce(ji,jj) &
222	& + zfnsol + fdtcn(ji,jj) - zfontn &
223	& + ( 1.0 - zindb ) * fsbbq(ji,jj) ) &
224	& * ( 1.0 - at_i(ji,jj) ) * rdt_ice
225
226	! Positive heat budget is used for bottom ablation
227	zfntlat = 1.0 - MAX( zzero , SIGN( zone , - qldif(ji,jj) ) )
228	!= 1 if positive heat budget
229	zpareff = 1.0 - zinda * zfntlat
230	!= 0 if ice and positive heat budget and 1 if one of those two is
231	!false
232	zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / &
233	MAX( at_i(ji,jj) * rdt_ice , epsi16 )
234
235	! Heat budget of the lead, energy transferred from ice to ocean
236	qldif (ji,jj) = zpareff * qldif(ji,jj)
237	qdtcn (ji,jj) = zpareff * qdtcn(ji,jj)
238
239	! Energy needed to bring ocean surface layer until its freezing
240	! qcmif, limflx
241	qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - sst_io(ji,jj) ) * ( 1. - zinda )
242
243	! calculate oceanic heat flux (limthd_dh)
244	fbif (ji,jj) = zindb * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) )
245
246	! computation of the daily thermodynamic ice production (only needed for output)
247	zhicifp(ji,jj) = ht_i(ji,jj,1) * at_i(ji,jj)
248	zhicifp(ji,jj) = ht_i(ji,jj,1) * at_i(ji,jj)
249	END DO
250	END DO
251
252	!------------------------------------------------------------------------------!
253	! 3) Select icy points and fulfill arrays for the vectorial grid.
254	!------------------------------------------------------------------------------!
255
256	DO jl = 1, jpl !loop over ice categories
257
258	WRITE(numout,*) ' lim_thd : transfer to 1D vectors. Category no : ', jl
259	WRITE(numout,*) ' ~~~~~~~~'
260
261	zareamin = 1.0e-10
262	nbpb = 0
263	DO jj = 1, jpj
264	DO ji = 1, jpi
265	IF ( a_i(ji,jj,jl) .gt. zareamin ) THEN
266	nbpb = nbpb + 1
267	npb(nbpb) = (jj - 1) * jpi + ji
268	ENDIF
269	! debug point to follow
270	IF ( (ji.eq.jiindex).AND.(jj.eq.jjindex) ) THEN
271	jiindex_1d = nbpb
272	ENDIF
273	END DO
274	END DO
275
276	!------------------------------------------------------------------------------!
277	! 4) Thermodynamic computation
278	!------------------------------------------------------------------------------!
279
280	IF (nbpb > 0) THEN ! If there is no ice, do nothing.
281
282	!-------------------------
283	! 4.1 Move to 1D arrays
284	!-------------------------
285
286	CALL tab_2d_1d( nbpb, at_i_b (1:nbpb) , at_i , jpi, jpj, npb(1:nbpb) )
287	CALL tab_2d_1d( nbpb, a_i_b (1:nbpb) , a_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
288	CALL tab_2d_1d( nbpb, ht_i_b (1:nbpb) , ht_i(:,:,jl) , jpi, jpj, npb(1:nbpb) )
289	CALL tab_2d_1d( nbpb, ht_s_b (1:nbpb) , ht_s(:,:,jl) , jpi, jpj, npb(1:nbpb) )
290
291	CALL tab_2d_1d( nbpb, t_su_b (1:nbpb) , t_su(:,:,jl), jpi, jpj, npb(1:nbpb) )
292	CALL tab_2d_1d( nbpb, sm_i_b (1:nbpb) , sm_i(:,:,jl), jpi, jpj, npb(1:nbpb) )
293	DO jk = 1, nlay_s
294	CALL tab_2d_1d( nbpb, t_s_b(1:nbpb,jk) , t_s(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
295	CALL tab_2d_1d( nbpb, q_s_b(1:nbpb,jk) , e_s(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
296	END DO
297	DO jk = 1, nlay_i
298	CALL tab_2d_1d( nbpb, t_i_b(1:nbpb,jk) , t_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
299	CALL tab_2d_1d( nbpb, q_i_b(1:nbpb,jk) , e_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
300	CALL tab_2d_1d( nbpb, s_i_b(1:nbpb,jk) , s_i(:,:,jk,jl), jpi, jpj, npb(1:nbpb) )
301	END DO
302
303	CALL tab_2d_1d( nbpb, tatm_ice_1d(1:nbpb) , tatm_ice(:,:) , jpi, jpj, npb(1:nbpb) )
304	CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb) , qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
305	CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb) , fr1_i0 , jpi, jpj, npb(1:nbpb) )
306	CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb) , fr2_i0 , jpi, jpj, npb(1:nbpb) )
307	CALL tab_2d_1d( nbpb, qnsr_ice_1d(1:nbpb) , qnsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
308
309	#if ! defined key_coupled
310	CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb) , qla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
311	CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb) , dqla_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
312	#endif
313
314	CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb) , dqns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
315	CALL tab_2d_1d( nbpb, t_bo_b (1:nbpb) , t_bo , jpi, jpj, npb(1:nbpb) )
316	CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb) , sprecip , jpi, jpj, npb(1:nbpb) )
317	CALL tab_2d_1d( nbpb, fbif_1d (1:nbpb) , fbif , jpi, jpj, npb(1:nbpb) )
318	CALL tab_2d_1d( nbpb, qldif_1d (1:nbpb) , qldif , jpi, jpj, npb(1:nbpb) )
319	CALL tab_2d_1d( nbpb, rdmicif_1d (1:nbpb) , rdmicif , jpi, jpj, npb(1:nbpb) )
320	CALL tab_2d_1d( nbpb, dmgwi_1d (1:nbpb) , dmgwi , jpi, jpj, npb(1:nbpb) )
321	CALL tab_2d_1d( nbpb, qlbbq_1d (1:nbpb) , zqlbsbq , jpi, jpj, npb(1:nbpb) )
322
323	CALL tab_2d_1d( nbpb, fseqv_1d (1:nbpb) , fseqv , jpi, jpj, npb(1:nbpb) )
324	CALL tab_2d_1d( nbpb, fsbri_1d (1:nbpb) , fsbri , jpi, jpj, npb(1:nbpb) )
325	CALL tab_2d_1d( nbpb, fhbri_1d (1:nbpb) , fhbri , jpi, jpj, npb(1:nbpb) )
326	CALL tab_2d_1d( nbpb, fstbif_1d (1:nbpb) , fstric , jpi, jpj, npb(1:nbpb) )
327	CALL tab_2d_1d( nbpb, qfvbq_1d (1:nbpb) , qfvbq , jpi, jpj, npb(1:nbpb) )
328
329	!--------------------------------
330	! 4.3) Thermodynamic processes
331	!--------------------------------
332
333	IF ( con_i ) CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting
334	IF ( con_i ) CALL lim_thd_glohec( qt_i_in , qt_s_in , &
335	q_i_layer_in , 1 , nbpb , jl )
336
337	!---------------------------------!
338	CALL lim_thd_dif(1,nbpb,jl) ! Ice/Snow Temperature profile !
339	!---------------------------------!
340
341	CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting
342	! compulsory for limthd_dh
343
344	IF ( con_i ) CALL lim_thd_glohec( qt_i_fin , qt_s_fin , &
345	q_i_layer_fin , 1 , nbpb , jl )
346	IF ( con_i ) CALL lim_thd_con_dif( 1 , nbpb , jl )
347
348	!---------------------------------!
349	CALL lim_thd_dh(1,nbpb,jl) ! Ice/Snow thickness !
350	!---------------------------------!
351
352	!---------------------------------!
353	CALL lim_thd_ent(1,nbpb,jl) ! Ice/Snow enthalpy remapping !
354	!---------------------------------!
355
356	!---------------------------------!
357	CALL lim_thd_sal(1,nbpb,jl) ! Ice salinity computation !
358	!---------------------------------!
359
360	! CALL lim_thd_enmelt(1,nbpb) ! computes sea ice energy of melting
361	IF ( con_i ) CALL lim_thd_glohec( qt_i_fin, qt_s_fin, &
362	q_i_layer_fin , 1 , nbpb , jl )
363	IF ( con_i ) CALL lim_thd_con_dh ( 1 , nbpb , jl )
364
365	!--------------------------------
366	! 4.4) Move 1D to 2D vectors
367	!--------------------------------
368
369	CALL tab_1d_2d( nbpb, at_i , npb, at_i_b (1:nbpb), jpi, jpj )
370	CALL tab_1d_2d( nbpb, ht_i(:,:,jl), npb, ht_i_b(1:nbpb), jpi, jpj )
371	CALL tab_1d_2d( nbpb, ht_s(:,:,jl), npb, ht_s_b(1:nbpb), jpi, jpj )
372	CALL tab_1d_2d( nbpb, a_i (:,:,jl), npb, a_i_b(1:nbpb) , jpi, jpj )
373	CALL tab_1d_2d( nbpb, t_su(:,:,jl), npb, t_su_b(1:nbpb), jpi, jpj )
374	CALL tab_1d_2d( nbpb, sm_i(:,:,jl), npb, sm_i_b(1:nbpb), jpi, jpj )
375
376	DO jk = 1, nlay_s
377	CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_b(1:nbpb,jk), jpi, jpj)
378	CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_b(1:nbpb,jk), jpi, jpj)
379	END DO
380
381	DO jk = 1, nlay_i
382	CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_b(1:nbpb,jk), jpi, jpj)
383	CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_b(1:nbpb,jk), jpi, jpj)
384	CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_b(1:nbpb,jk), jpi, jpj)
385	END DO
386
387	CALL tab_1d_2d( nbpb, fstric , npb, fstbif_1d (1:nbpb) , jpi, jpj )
388	CALL tab_1d_2d( nbpb, qldif , npb, qldif_1d (1:nbpb) , jpi, jpj )
389	CALL tab_1d_2d( nbpb, qfvbq , npb, qfvbq_1d (1:nbpb) , jpi, jpj )
390	CALL tab_1d_2d( nbpb, rdmicif , npb, rdmicif_1d(1:nbpb) , jpi, jpj )
391	CALL tab_1d_2d( nbpb, dmgwi , npb, dmgwi_1d (1:nbpb) , jpi, jpj )
392	CALL tab_1d_2d( nbpb, rdmsnif , npb, rdmsnif_1d(1:nbpb) , jpi, jpj )
393	CALL tab_1d_2d( nbpb, rdvosif , npb, dvsbq_1d (1:nbpb) , jpi, jpj )
394	CALL tab_1d_2d( nbpb, rdvobif , npb, dvbbq_1d (1:nbpb) , jpi, jpj )
395	CALL tab_1d_2d( nbpb, fdvolif , npb, dvlbq_1d (1:nbpb) , jpi, jpj )
396	CALL tab_1d_2d( nbpb, rdvonif , npb, dvnbq_1d (1:nbpb) , jpi, jpj )
397	CALL tab_1d_2d( nbpb, fseqv , npb, fseqv_1d (1:nbpb) , jpi, jpj )
398
399	IF ( num_sal .EQ. 2 ) THEN
400	CALL tab_1d_2d( nbpb, fsbri , npb, fsbri_1d (1:nbpb) , jpi, jpj )
401	CALL tab_1d_2d( nbpb, fhbri , npb, fhbri_1d (1:nbpb) , jpi, jpj )
402	ENDIF
403
404	!+++++
405	!temporary stuff for a dummyversion
406	CALL tab_1d_2d( nbpb, dh_i_surf2D, npb, dh_i_surf(1:nbpb) , jpi, jpj )
407	CALL tab_1d_2d( nbpb, dh_i_bott2D, npb, dh_i_bott(1:nbpb) , jpi, jpj )
408	CALL tab_1d_2d( nbpb, fsup2D , npb, fsup (1:nbpb) , jpi, jpj )
409	CALL tab_1d_2d( nbpb, focea2D , npb, focea (1:nbpb) , jpi, jpj )
410	CALL tab_1d_2d( nbpb, s_i_newice , npb, s_i_new (1:nbpb) , jpi, jpj )
411	CALL tab_1d_2d( nbpb, izero(:,:,jl) , npb, i0 (1:nbpb) , jpi, jpj )
412	CALL tab_1d_2d( nbpb, qnsr_ice(:,:,jl), npb, qnsr_ice_1d(1:nbpb), jpi, jpj)
413	!+++++
414
415	ENDIF ! nbpb
416
417	END DO ! jl
418
419	!------------------------------------------------------------------------------!
420	! 5) Global variables, diagnostics
421	!------------------------------------------------------------------------------!
422
423	!------------------------
424	! 5.1) Ice heat content
425	!------------------------
426
427	! Enthalpies are global variables we have to readjust the units
428	DO jl = 1, jpl
429	DO jk = 1, nlay_i
430	DO jj = 1, jpj
431	DO ji = 1, jpi
432	! Change dimensions
433	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac
434
435	! Multiply by volume, divide by nlayers so that heat content in 10^9 Joules
436	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &
437	area(ji,jj) * a_i(ji,jj,jl) * &
438	ht_i(ji,jj,jl) / nlay_i
439	END DO !ji
440	END DO !jj
441	END DO !jk
442	END DO !jl
443
444	!------------------------
445	! 5.2) Snow heat content
446	!------------------------
447
448	! Enthalpies are global variables we have to readjust the units
449	DO jl = 1, jpl
450	DO jk = 1, nlay_s
451	DO jj = 1, jpj
452	DO ji = 1, jpi
453	! Change dimensions
454	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
455	! Multiply by volume, so that heat content in 10^9 Joules
456	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * &
457	a_i(ji,jj,jl) * ht_s(ji,jj,jl) / nlay_s
458	END DO !ji
459	END DO !jj
460	END DO !jk
461	END DO !jl
462
463	!----------------------------------
464	! 5.3) Change thickness to volume
465	!----------------------------------
466	CALL lim_var_eqv2glo
467
468	!--------------------------------------------
469	! 5.4) Diagnostic thermodynamic growth rates
470	!--------------------------------------------
471	d_v_i_thd (:,:,:) = v_i(:,:,:) - old_v_i(:,:,:) ! ice volumes
472	dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) / rdt_ice * 86400.0
473
474	IF ( con_i ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)
475
476	END SUBROUTINE lim_thd
477
478	!===============================================================================
479
480	SUBROUTINE lim_thd_glohec(eti,ets,etilayer,kideb,kiut,jl)
481	!!-----------------------------------------------------------------------
482	!! * ROUTINE lim_thd_glohec *
483	!!
484	!! ** Purpose : Compute total heat content for each category
485	!! Works with 1d vectors only
486	!!
487	!! history :
488	!! 9.9 ! 07-04 (M.Vancoppenolle) original code
489	!!-----------------------------------------------------------------------
490	!! * Local variables
491	INTEGER, INTENT(in) :: &
492	kideb, kiut, & ! bounds for the spatial loop
493	jl ! category number
494
495	REAL(wp), DIMENSION (jpij,jpl), INTENT(out) :: &
496	eti, ets ! vertically-summed heat content for ice /snow
497
498	REAL(wp), DIMENSION (jpij,jkmax), INTENT(out) :: &
499	etilayer ! heat content for ice layers
500
501	REAL(wp) :: &
502	zdes, & ! snow heat content increment (dummy)
503	zeps ! very small value (1.e-10)
504
505	INTEGER :: &
506	ji,jj,jk ! loop indices
507
508	!!-----------------------------------------------------------------------
509	eti(:,:) = 0.0
510	ets(:,:) = 0.0
511	zeps = 1.0e-10
512
513	! total q over all layers, ice [J.m-2]
514	DO jk = 1, nlay_i
515	DO ji = kideb, kiut
516	etilayer(ji,jk) = q_i_b(ji,jk) &
517	* ht_i_b(ji) / nlay_i
518	eti(ji,jl) = eti(ji,jl) + etilayer(ji,jk)
519	END DO
520	END DO
521
522	! total q over all layers, snow [J.m-2]
523	DO ji = kideb, kiut
524	zdes = q_s_b(ji,1) * ht_s_b(ji) / nlay_s
525	ets(ji,jl) = ets(ji,jl) + zdes
526	END DO
527
528	WRITE(numout,*) ' lim_thd_glohec '
529	WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) / rdt_ice
530	WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) / rdt_ice
531	WRITE(numout,*) ' qt_in : ', ( eti(jiindex_1d,jl) + &
532	ets(jiindex_1d,jl) ) / rdt_ice
533
534	END SUBROUTINE lim_thd_glohec
535
536	!===============================================================================
537
538	SUBROUTINE lim_thd_con_dif(kideb,kiut,jl)
539	!!-----------------------------------------------------------------------
540	!! * ROUTINE lim_thd_con_dif *
541	!!
542	!! ** Purpose : Test energy conservation after heat diffusion
543	!!
544	!! history :
545	!! 9.9 ! 07-04 (M.Vancoppenolle) original code
546	!!-------------------------------------------------------------------
547	!! * Local variables
548	INTEGER, INTENT(in) :: &
549	kideb, kiut, & !: bounds for the spatial loop
550	jl !: category number
551
552	REAL(wp) :: & !: ! goes to trash
553	meance, & !: mean conservation error
554	max_cons_err, & !: maximum tolerated conservation error
555	max_surf_err !: maximum tolerated surface error
556
557	INTEGER :: &
558	numce !: number of points for which conservation
559	! is violated
560	INTEGER :: &
561	ji,jj,jk, & !: loop indices
562	zji, zjj
563	!!---------------------------------------------------------------------
564
565	max_cons_err = 1.0
566	max_surf_err = 0.001
567
568	!--------------------------
569	! Increment of energy
570	!--------------------------
571	! global
572	DO ji = kideb, kiut
573	dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl) &
574	+ qt_s_fin(ji,jl) - qt_s_in(ji,jl)
575	END DO
576	! layer by layer
577	dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:)
578
579	!----------------------------------------
580	! Atmospheric heat flux, ice heat budget
581	!----------------------------------------
582
583	DO ji = kideb, kiut
584	zji = MOD( npb(ji) - 1, jpi ) + 1
585	zjj = ( npb(ji) - 1 ) / jpi + 1
586
587	fatm(ji,jl) = &
588	qnsr_ice_1d(ji) + & ! atm non solar
589	(1.0-i0(ji))*qsr_ice_1d(ji) ! atm solar
590
591	sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji)*i0(ji) &
592	- fstroc(zji,zjj,jl)
593	END DO
594
595	!--------------------
596	! Conservation error
597	!--------------------
598
599	DO ji = kideb, kiut
600	cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) )
601	END DO
602
603	numce = 0
604	meance = 0.0
605	DO ji = kideb, kiut
606	IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
607	numce = numce + 1
608	meance = meance + cons_error(ji,jl)
609	ENDIF
610	ENDDO
611	IF (numce .GT. 0 ) meance = meance / numce
612
613	WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err
614	WRITE(numout,*) ' After lim_thd_dif, category : ', jl
615	WRITE(numout,*) ' Mean conservation error on big error points ', meance, &
616	numit
617	WRITE(numout,*) ' Number of points where there is a cons err gt than c.e. : ', numce, numit
618
619	!-------------------------------------------------------
620	! Surface error due to imbalance between Fatm and Fcsu
621	!-------------------------------------------------------
622	numce = 0.0
623	meance = 0.0
624
625	DO ji = kideb, kiut
626	surf_error(ji,jl) = ABS ( fatm(ji,jl) - fc_su(ji) )
627	IF ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. &
628	max_surf_err ) ) THEN
629	numce = numce + 1
630	meance = meance + surf_error(ji,jl)
631	ENDIF
632	ENDDO
633	IF (numce .GT. 0 ) meance = meance / numce
634
635	WRITE(numout,*) ' Maximum tolerated surface error : ', max_surf_err
636	WRITE(numout,*) ' After lim_thd_dif, category : ', jl
637	WRITE(numout,*) ' Mean surface error on big error points ', meance, numit
638	WRITE(numout,*) ' Number of points where there is a surf err gt than surf_err : ', numce, numit
639
640	IF (jiindex_1D.GT.0) WRITE(numout,*) ' fc_su : ', fc_su(jiindex_1d)
641	IF (jiindex_1D.GT.0) WRITE(numout,*) ' fatm : ', fatm(jiindex_1d,jl)
642	IF (jiindex_1D.GT.0) WRITE(numout,*) ' t_su : ', t_su_b(jiindex_1d)
643
644	!---------------------------------------
645	! Write ice state in case of big errors
646	!---------------------------------------
647
648	DO ji = kideb, kiut
649	IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. &
650	( cons_error(ji,jl) .GT. max_cons_err ) ) THEN
651	zji = MOD( npb(ji) - 1, jpi ) + 1
652	zjj = ( npb(ji) - 1 ) / jpi + 1
653
654	WRITE(numout,*) ' alerte 1 '
655	WRITE(numout,*) ' Untolerated conservation / surface error after '
656	WRITE(numout,*) ' heat diffusion in the ice '
657	WRITE(numout,*) ' Category : ', jl
658	WRITE(numout,*) ' zji , zjj : ', zji, zjj
659	WRITE(numout,*) ' lat, lon : ', gphit(zji,zjj), glamt(zji,zjj)
660	WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)
661	WRITE(numout,*) ' surf_error : ', surf_error(ji,jl)
662	WRITE(numout,*) ' dq_i : ', - dq_i(ji,jl) / rdt_ice
663	WRITE(numout,*) ' Fdt : ', sum_fluxq(ji,jl)
664	WRITE(numout,*)
665	! WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl)
666	! WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl)
667	! WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl)
668	! WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl)
669	! WRITE(numout,*) ' qt : ', qt_i_fin(ji,jl) + &
670	! qt_s_fin(ji,jl)
671	WRITE(numout,*) ' ht_i : ', ht_i_b(ji)
672	WRITE(numout,*) ' ht_s : ', ht_s_b(ji)
673	WRITE(numout,*) ' t_su : ', t_su_b(ji)
674	WRITE(numout,*) ' t_s : ', t_s_b(ji,1)
675	WRITE(numout,*) ' t_i : ', t_i_b(ji,1:nlay_i)
676	WRITE(numout,*) ' t_bo : ', t_bo_b(ji)
677	WRITE(numout,*) ' q_i : ', q_i_b(ji,1:nlay_i)
678	WRITE(numout,*) ' s_i : ', s_i_b(ji,1:nlay_i)
679	WRITE(numout,) ' tmelts : ', rtt - tmuts_i_b(ji,1:nlay_i)
680	WRITE(numout,*)
681	WRITE(numout,*) ' Fluxes '
682	WRITE(numout,*) ' ~~~~~~ '
683	WRITE(numout,*) ' fatm : ', fatm(ji,jl)
684	WRITE(numout,*) ' fc_su : ', fc_su (ji)
685	WRITE(numout,) ' fstr_inice : ', qsr_ice_1d(ji)i0(ji)
686	WRITE(numout,*) ' fc_bo : ', - fc_bo_i (ji)
687	WRITE(numout,*) ' foc : ', fbif_1d(ji)
688	WRITE(numout,*) ' fstroc : ', fstroc (zji,zjj,jl)
689	WRITE(numout,*) ' i0 : ', i0(ji)
690	WRITE(numout,) ' fsolar : ', (1.0-i0(ji))qsr_ice_1d(ji)
691	WRITE(numout,*) ' fnsolar : ', qnsr_ice_1d(ji)
692	WRITE(numout,*) ' Conduction fluxes : '
693	WRITE(numout,*) ' fc_s : ', fc_s(ji,0:nlay_s)
694	WRITE(numout,*) ' fc_i : ', fc_i(ji,0:nlay_i)
695	WRITE(numout,*)
696	WRITE(numout,*) ' Layer by layer ... '
697	WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - &
698	qt_s_in(ji,jl) ) &
699	/ rdt_ice
700	WRITE(numout,*) ' dfc_snow : ', fc_s(ji,1) - &
701	fc_s(ji,0)
702	DO jk = 1, nlay_i
703	WRITE(numout,*) ' layer : ', jk
704	WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) / rdt_ice
705	WRITE(numout,*) ' radab : ', radab(ji,jk)
706	WRITE(numout,*) ' dfc_i : ', fc_i(ji,jk) - &
707	fc_i(ji,jk-1)
708	WRITE(numout,*) ' tot f : ', fc_i(ji,jk) - &
709	fc_i(ji,jk-1) - radab(ji,jk)
710	END DO
711
712	ENDIF
713
714	END DO
715
716	END SUBROUTINE lim_thd_con_dif
717
718	!==============================================================================
719
720	SUBROUTINE lim_thd_con_dh(kideb,kiut,jl)
721	!!-----------------------------------------------------------------------
722	!! * ROUTINE lim_thd_con_dh *
723	!!
724	!! ** Purpose : Test energy conservation after enthalpy redistr.
725	!!
726	!! history :
727	!! 9.9 ! 07-04 (M.Vancoppenolle) original code
728	!!-----------------------------------------------------------------------
729	!! * Local variables
730	INTEGER, INTENT(in) :: &
731	kideb, kiut, & !: bounds for the spatial loop
732	jl !: category number
733
734	REAL(wp) :: & !: ! goes to trash
735	meance, & !: mean conservation error
736	max_cons_err !: maximum tolerated conservation error
737
738	INTEGER :: &
739	numce !: number of points for which conservation
740	! is violated
741	INTEGER :: &
742	ji,jj,jk, & !: loop indices
743	zji, zjj
744
745	!!---------------------------------------------------------------------
746
747	max_cons_err = 1.0
748
749	!--------------------------
750	! Increment of energy
751	!--------------------------
752	! global
753	DO ji = kideb, kiut
754	dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl) &
755	+ qt_s_fin(ji,jl) - qt_s_in(ji,jl)
756	END DO
757	! layer by layer
758	dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:)
759
760	!----------------------------------------
761	! Atmospheric heat flux, ice heat budget
762	!----------------------------------------
763
764	DO ji = kideb, kiut
765	zji = MOD( npb(ji) - 1, jpi ) + 1
766	zjj = ( npb(ji) - 1 ) / jpi + 1
767
768	fatm(ji,jl) = &
769	qnsr_ice_1d(ji) + & ! atm non solar
770	! (1.0-i0(ji))*qsr_ice_1d(ji) ! atm solar
771	qsr_ice_1d(ji) ! atm solar
772
773	sum_fluxq(ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) &
774	- fstroc(zji,zjj,jl)
775	cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) )
776	END DO
777
778	!--------------------
779	! Conservation error
780	!--------------------
781
782	DO ji = kideb, kiut
783	cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) )
784	END DO
785
786	numce = 0
787	meance = 0.0
788	DO ji = kideb, kiut
789	IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
790	numce = numce + 1
791	meance = meance + cons_error(ji,jl)
792	ENDIF
793	ENDDO
794	IF (numce .GT. 0 ) meance = meance / numce
795
796	WRITE(numout,*) ' Error report - Category : ', jl
797	WRITE(numout,*) ' ~~~~~~~~~~~~ '
798	WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err
799	WRITE(numout,*) ' After lim_thd_ent, category : ', jl
800	WRITE(numout,*) ' Mean conservation error on big error points ', meance, &
801	numit
802	WRITE(numout,*) ' Number of points where there is a cons err gt than 0.1 W/m2 : ', numce, numit
803
804	!---------------------------------------
805	! Write ice state in case of big errors
806	!---------------------------------------
807
808	DO ji = kideb, kiut
809	IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
810	zji = MOD( npb(ji) - 1, jpi ) + 1
811	zjj = ( npb(ji) - 1 ) / jpi + 1
812
813	WRITE(numout,*) ' alerte 1 - category : ', jl
814	WRITE(numout,*) ' Untolerated conservation error after limthd_ent '
815	WRITE(numout,*) ' zji , zjj : ', zji, zjj
816	WRITE(numout,*) ' lat, lon : ', gphit(zji,zjj), glamt(zji,zjj)
817	WRITE(numout,) ' '
818	WRITE(numout,*) ' Ftotal : ', sum_fluxq(ji,jl)
819	WRITE(numout,*) ' dq_t : ', - dq_i(ji,jl) / rdt_ice
820	WRITE(numout,*) ' dq_i : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) / rdt_ice
821	WRITE(numout,*) ' dq_s : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice
822	WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)
823	WRITE(numout,) ' '
824	WRITE(numout,*) ' Fluxes --- : '
825	WRITE(numout,*) ' fatm : ', fatm(ji,jl)
826	WRITE(numout,*) ' foce : ', fbif_1d(ji)
827	WRITE(numout,*) ' fres : ', ftotal_fin(ji)
828	WRITE(numout,*) ' fhbri : ', fhbricat(zji,zjj,jl)
829	WRITE(numout,) ' '
830	WRITE(numout,*) ' Heat contents --- : '
831	WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl) / rdt_ice
832	WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl) / rdt_ice
833	WRITE(numout,*) ' qt_in : ', ( qt_i_in(ji,jl) + &
834	qt_s_in(ji,jl) ) / rdt_ice
835	WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl) / rdt_ice
836	WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl) / rdt_ice
837	WRITE(numout,*) ' qt_fin : ', ( qt_i_fin(ji,jl) + &
838	qt_s_fin(ji,jl) ) / rdt_ice
839	WRITE(numout,) ' '
840	WRITE(numout,*) ' Ice variables --- : '
841	WRITE(numout,*) ' ht_i : ', ht_i_b(ji)
842	WRITE(numout,*) ' ht_s : ', ht_s_b(ji)
843	WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji)
844	WRITE(numout,*) ' dh_snowice: ', dh_snowice(ji)
845	WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji)
846	WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
847
848	ENDIF
849
850	END DO
851
852	END SUBROUTINE lim_thd_con_dh
853	!==============================================================================
854
855	SUBROUTINE lim_thd_enmelt(kideb,kiut)
856	!!-----------------------------------------------------------------------
857	!! * ROUTINE lim_thd_enmelt *
858	!!
859	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3)
860	!!
861	!! ** Method : Formula (Bitz and Lipscomb, 1999)
862	!!
863	!! history : Martin Vancoppenolle, May 2007
864	!!-------------------------------------------------------------------
865	INTEGER, INTENT(in) :: &
866	kideb, kiut !: bounds for the spatial loop
867
868	REAL(wp) :: & !: goes to trash
869	ztmelts , & !: sea ice freezing point in K
870	zeps
871
872	INTEGER :: &
873	ji, & !: spatial loop index
874	jk !: vertical index
875
876	!!-------------------------------------------------------------------
877	zeps = 1.0e-10
878
879	! Sea ice energy of melting
880	DO jk = 1, nlay_i
881	DO ji = kideb, kiut
882	ztmelts = - tmut * s_i_b(ji,jk) + rtt
883	q_i_b(ji,jk) = rhoic( cpic ( ztmelts - t_i_b(ji,jk) ) &
884	+ lfus * ( 1.0 - (ztmelts-rtt)/MIN((t_i_b(ji,jk)-rtt),-zeps) ) &
885	- rcp * ( ztmelts-rtt ) )
886	END DO !ji
887	END DO !jk
888
889	! Snow energy of melting
890	DO jk = 1, nlay_s
891	DO ji = kideb,kiut
892	q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
893	END DO !ji
894	END DO !jk
895
896	END SUBROUTINE lim_thd_enmelt
897
898	!==============================================================================
899
900	SUBROUTINE lim_thd_init
901
902	!!-----------------------------------------------------------------------
903	!! * ROUTINE lim_thd_init *
904	!!
905	!! ** Purpose : Physical constants and parameters linked to the ice
906	!! thermodynamics
907	!!
908	!! ** Method : Read the namicethd namelist and check the ice-thermo
909	!! parameter values called at the first timestep (nit000)
910	!!
911	!! ** input : Namelist namicether
912	!!
913	!! history :
914	!! 8.5 ! 03-08 (C. Ethe) original code
915	!!-------------------------------------------------------------------
916	NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, &
917	& hicmin, hiclim, amax , &
918	& sbeta , parlat, hakspl, hibspl, exld, &
919	& hakdif, hnzst , thth , parsub, alphs, betas, &
920	& kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi
921	!!-------------------------------------------------------------------
922
923	! Define the initial parameters
924	! -------------------------
925	REWIND( numnam_ice )
926	READ ( numnam_ice , namicethd )
927	IF (lwp) THEN
928	WRITE(numout,*)
929	WRITE(numout,*)'lim_thd_init : ice parameters for ice thermodynamic computation '
930	WRITE(numout,*)'~~~~~~~~~~~~'
931	WRITE(numout,*)' maximum melting at the bottom hmelt = ', hmelt
932	WRITE(numout,*)' ice thick. for lateral accretion in NH (SH) hiccrit(1/2) = ', hiccrit
933	WRITE(numout,*)' Frazil ice thickness as a function of wind or not fraz_swi = ', fraz_swi
934	WRITE(numout,*)' Maximum proportion of frazil ice collecting at bottom maxfrazb = ', maxfrazb
935	WRITE(numout,*)' Thresold relative drift speed for collection of frazil vfrazb = ', vfrazb
936	WRITE(numout,*)' Squeezing coefficient for collection of frazil Cfrazb = ', Cfrazb
937	WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin
938	WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim
939	WRITE(numout,*)' maximum lead fraction amax = ', amax
940	WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice '
941	WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta
942	WRITE(numout,*)' percentage of energy used for lateral ablation parlat = ', parlat
943	WRITE(numout,*)' slope of distr. for Hakkinen-Mellor lateral melting hakspl = ', hakspl
944	WRITE(numout,*)' slope of distribution for Hibler lateral melting hibspl = ', hibspl
945	WRITE(numout,*)' exponent for leads-closure rate exld = ', exld
946	WRITE(numout,*)' coefficient for diffusions of ice and snow hakdif = ', hakdif
947	WRITE(numout,*)' threshold thick. for comp. of eq. thermal conductivity zhth = ', thth
948	WRITE(numout,*)' thickness of the surf. layer in temp. computation hnzst = ', hnzst
949	WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub
950	WRITE(numout,*)' coefficient for snow density when snow ice formation alphs = ', alphs
951	WRITE(numout,*)' coefficient for ice-lead partition of snowfall betas = ', betas
952	WRITE(numout,*)' extinction radiation parameter in sea ice (1.0) kappa_i = ', kappa_i
953	WRITE(numout,*)' maximal n. of iter. for heat diffusion computation nconv_i_thd = ', nconv_i_thd
954	WRITE(numout,*)' maximal err. on T for heat diffusion computation maxer_i_thd = ', maxer_i_thd
955	WRITE(numout,*)' switch for comp. of thermal conductivity in the ice thcon_i_swi = ', thcon_i_swi
956	WRITE(numout,*)
957	ENDIF
958
959	rcdsn = hakdif * rcdsn
960	rcdic = hakdif * rcdic
961
962
963	END SUBROUTINE lim_thd_init
964
965	#else
966	!!======================================================================
967	!! * MODULE limthd *
968	!! No sea ice model
969	!!======================================================================
970	CONTAINS
971	SUBROUTINE lim_thd ! Empty routine
972	END SUBROUTINE lim_thd
973	SUBROUTINE lim_thd_con_dif
974	END SUBROUTINE lim_thd_con_dif
975	#endif
976
977	!!======================================================================
978	END MODULE limthd

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

Download in other formats: