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 branches/dev_002_LIM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

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

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

Download in other formats: