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

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

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

Download in other formats: