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 @ 913

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

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

Download in other formats: