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 – NEMO

Context Navigation

source: trunk/NEMO/LIM_SRC/limthd.F90 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 23.7 KB

Line
1	MODULE limthd
2	!!======================================================================
3	!! * MODULE limthd *
4	!! LIM thermo ice model : ice thermodynamic
5	!!======================================================================
6	#if defined key_ice_lim
7	!!----------------------------------------------------------------------
8	!! 'key_ice_lim' : 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
19	USE ice ! LIM sea-ice variables
20	USE ice_oce ! sea-ice/ocean variables
21	USE flx_oce ! sea-ice/ocean forcings variables
22	USE thd_ice ! LIM thermodynamic sea-ice variables
23	USE dom_ice ! LIM sea-ice domain
24	USE iceini
25	USE limthd_zdf
26	USE limthd_lac
27	USE limtab
28
29	IMPLICIT NONE
30	PRIVATE
31
32	!! * Routine accessibility
33	PUBLIC lim_thd ! called by lim_step
34
35	!! * Module variables
36	REAL(wp) :: & ! constant values
37	epsi20 = 1e-20 , &
38	epsi16 = 1e-16 , &
39	epsi04 = 1e-04 , &
40	zzero = 0.0 , &
41	zone = 1.0
42	!! * Substitutions
43	# include "vectopt_loop_substitute.h90"
44	!!-------- -------------------------------------------------------------
45	!! LIM 2.0 , UCL-LODYC-IPSL (2003)
46	!!----------------------------------------------------------------------
47
48	CONTAINS
49
50	SUBROUTINE lim_thd
51	!!-------------------------------------------------------------------
52	!! * ROUTINE lim_thd *
53	!!
54	!! ** Purpose : This routine manages the ice thermodynamic.
55	!!
56	!! ** Action : - Initialisation of some variables
57	!! - Some preliminary computation (oceanic heat flux
58	!! at the ice base, snow acc.,heat budget of the leads)
59	!! - selection of the icy points and put them in an array
60	!! - call lim_vert_ther for vert ice thermodynamic
61	!! - back to the geographic grid
62	!! - selection of points for lateral accretion
63	!! - call lim_lat_acc for the ice accretion
64	!! - back to the geographic grid
65	!!
66	!!
67	!! ** References :
68	!! H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
69	!!
70	!! History :
71	!! 1.0 ! 00-01 (LIM)
72	!! 2.0 ! 02-07 (C. Ethe, G. Madec) F90
73	!!---------------------------------------------------------------------
74	!! * Local variables
75	INTEGER :: ji, jj, & ! dummy loop indices
76	nbpb , & ! nb of icy pts for thermo. cal.
77	nbpac ! nb of pts for lateral accretion
78
79	REAL(wp) :: &
80	zfric_umin = 5e-03 , & ! lower bound for the friction velocity
81	zfric_umax = 2e-02 ! upper bound for the friction velocity
82
83	REAL(wp) :: &
84	zinda , & ! switch for test. the val. of concen.
85	zindb, zindg , & ! switches for test. the val of arg
86	za , zh, zthsnice , &
87	zfric_u , & ! friction velocity
88	zfnsol , & ! total non solar heat
89	zfontn , & ! heat flux from snow thickness
90	zfntlat, zpareff ! test. the val. of lead heat budget
91
92	REAL(wp), DIMENSION(jpi,jpj) :: &
93	zhicifp , & ! ice thickness for outputs
94	zqlbsbq ! link with lead energy budget qldif
95	!!-------------------------------------------------------------------
96
97	IF( numit == nstart ) CALL lim_thd_init ! Initialization (first time-step only)
98
99	!-------------------------------------------!
100	! Initilization of diagnostic variables !
101	!-------------------------------------------!
102
103	!i est-ce utile? oui au moins en partie
104	rdvosif(:,:) = 0.0 ! variation of ice volume at surface
105	rdvobif(:,:) = 0.0 ! variation of ice volume at bottom
106	fdvolif(:,:) = 0.0 ! total variation of ice volume
107	rdvonif(:,:) = 0.0 ! lateral variation of ice volume
108	fstric (:,:) = 0.0 ! part of solar radiation absorbing inside the ice
109	fscmbq (:,:) = 0.0 ! linked with fstric
110	ffltbif(:,:) = 0.0 ! linked with fstric
111	qfvbq (:,:) = 0.0 ! linked with fstric
112	rdmsnif(:,:) = 0.0 ! variation of snow mass per unit area
113	rdmicif(:,:) = 0.0 ! variation of ice mass per unit area
114	hicifp (:,:) = 0.0 ! daily thermodynamic ice production.
115
116	DO jj = 1, jpj
117	DO ji = 1, jpi
118	hsnif(ji,jj) = hsnif(ji,jj) * MAX( zzero, SIGN( zone , hsnif(ji,jj) - epsi04 ) )
119	END DO
120	END DO
121	IF( l_ctl .AND. lwp ) WRITE(numout,*) 'lim_thd : ', SUM( hsnif(:,:) ) , ' hsnif'
122
123
124	!-----------------------------------!
125	! Treatment of particular cases !
126	!-----------------------------------!
127
128	DO jj = 1, jpj
129	DO ji = 1, jpi
130	! snow is transformed into ice if the original ice cover disappears.
131	zindg = tms(ji,jj) * MAX( zzero , SIGN( zone , -hicif(ji,jj) ) )
132	hicif(ji,jj) = hicif(ji,jj) + zindg * rhosn * hsnif(ji,jj) / rau0
133	hsnif(ji,jj) = ( zone - zindg ) * hsnif(ji,jj) + zindg * hicif(ji,jj) * ( rau0 - rhoic ) / rhosn
134	dmgwi(ji,jj) = zindg * (1.0 - frld(ji,jj)) * rhoic * hicif(ji,jj) ! snow/ice mass
135
136	! the lead fraction, frld, must be little than or equal to amax (ice ridging).
137	zthsnice = hsnif(ji,jj) + hicif(ji,jj)
138	zindb = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) )
139	za = zindb * MIN( zone, ( 1.0 - frld(ji,jj) ) * uscomi )
140	hsnif (ji,jj) = hsnif(ji,jj) * za
141	hicif (ji,jj) = hicif(ji,jj) * za
142	qstoif(ji,jj) = qstoif(ji,jj) * za
143	frld (ji,jj) = 1.0 - zindb * ( 1.0 - frld(ji,jj) ) / MAX( za , epsi20 )
144
145	! the in situ ice thickness, hicif, must be equal to or greater than hiclim.
146	zh = MAX( zone , zindb * hiclim / MAX( hicif(ji,jj) , epsi20 ) )
147	hsnif (ji,jj) = hsnif(ji,jj) * zh
148	hicif (ji,jj) = hicif(ji,jj) * zh
149	qstoif(ji,jj) = qstoif(ji,jj) * zh
150	frld (ji,jj) = ( frld(ji,jj) + ( zh - 1.0 ) ) / zh
151	END DO
152	END DO
153	IF( l_ctl .AND. lwp ) THEN
154	WRITE(numout,*) 'lim_thd: hicif : ', SUM( hicif ), ' hsnif ', SUM( hsnif )
155	WRITE(numout,*) 'lim_thd: dmgwi : ', SUM( dmgwi ), ' qstoif ', SUM( qstoif )
156	WRITE(numout,*) 'lim_thd: frld : ', SUM( frld )
157	ENDIF
158
159
160	!-------------------------------!
161	! Thermodynamics of sea ice !
162	!-------------------------------!
163
164	! Partial computation of forcing for the thermodynamic sea ice model.
165	!--------------------------------------------------------------------------
166
167	!CDIR NOVERRCHK
168	DO jj = 1, jpj
169	!CDIR NOVERRCHK
170	DO ji = 1, jpi
171	zthsnice = hsnif(ji,jj) + hicif(ji,jj)
172	zindb = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) )
173	pfrld(ji,jj) = frld(ji,jj)
174	zinda = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) )
175
176	! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
177	thcm(ji,jj) = 0.0
178
179	! net downward heat flux from the ice to the ocean, expressed as a function of ocean
180	! temperature and turbulent mixing (McPhee, 1992)
181	zfric_u = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) ! friction velocity
182	fdtcn(ji,jj) = zindb * rau0 * rcp * 0.006 * zfric_u * ( sst_io(ji,jj) - tfu(ji,jj) )
183	qdtcn(ji,jj) = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice
184
185	! partial computation of the lead energy budget (qldif)
186	zfontn = ( sprecip(ji,jj) / rhosn ) * xlsn ! energy for melting
187	zfnsol = qnsr_oce(ji,jj) ! total non solar flux
188	qldif(ji,jj) = tms(ji,jj) * ( qsr_oce(ji,jj) * ( 1.0 - thcm(ji,jj) ) &
189	& + zfnsol + fdtcn(ji,jj) - zfontn &
190	& + ( 1.0 - zindb ) * fsbbq(ji,jj) ) &
191	& * frld(ji,jj) * rdt_ice
192	! parlat : percentage of energy used for lateral ablation (0.0)
193	zfntlat = 1.0 - MAX( zzero , SIGN( zone , - qldif(ji,jj) ) )
194	zpareff = 1.0 + ( parlat - 1.0 ) * zinda * zfntlat
195	zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( (1.0 - frld(ji,jj)) * rdt_ice , epsi16 )
196	qldif (ji,jj) = zpareff * qldif(ji,jj)
197	qdtcn (ji,jj) = zpareff * qdtcn(ji,jj)
198
199	! energy needed to bring ocean surface layer until its freezing
200	qcmif (ji,jj) = rau0 * rcp * dz * ( tfu(ji,jj) - sst_io(ji,jj) ) * ( 1 - zinda )
201
202	! calculate oceanic heat flux.
203	fbif (ji,jj) = zindb * ( fsbbq(ji,jj) / MAX( (1.0 - frld(ji,jj)) , epsi20 ) + fdtcn(ji,jj) )
204
205	! computation of the daily thermodynamic ice production (only needed for output)
206	zhicifp(ji,jj) = hicif(ji,jj) * ( 1.0 - frld(ji,jj) )
207	END DO
208	END DO
209
210
211	! Select icy points and fulfill arrays for the vectorial grid.
212	!----------------------------------------------------------------------
213	nbpb = 0
214	DO jj = 1, jpj
215	DO ji = 1, jpi
216	IF ( frld(ji,jj) < 1.0 ) THEN
217	nbpb = nbpb + 1
218	npb(nbpb) = (jj - 1) * jpi + ji
219	ENDIF
220	END DO
221	END DO
222	IF( l_ctl .AND. lwp ) THEN
223	WRITE(numout,*) 'lim_thd: pfrld ', SUM( pfrld ), ' thcm ', SUM( thcm )
224	WRITE(numout,*) 'lim_thd: fdtcn ', SUM( fdtcn ), ' qdtcn ', SUM( qdtcn )
225	WRITE(numout,*) 'lim_thd: qldif ', SUM( qldif ), ' zqlbsbq ', SUM( zqlbsbq )
226	WRITE(numout,*) 'lim_thd: qcmif ', SUM( qcmif ), ' fbif ', SUM( fbif )
227	WRITE(numout,) 'lim_thd: qcmif ', SUM( qcmiftms )
228	WRITE(numout,*) 'lim_thd: zhicifp', SUM( zhicifp )
229	WRITE(numout,*) 'limthd : nbpb = ', nbpb
230	ENDIF
231
232
233	! If there is no ice, do nothing. Otherwise, compute Top and Bottom accretion/ablation
234	!------------------------------------------------------------------------------------
235
236	IF ( nbpb > 0) THEN
237
238	! put the variable in a 1-D array for thermodynamics process
239	CALL tab_2d_1d( nbpb, frld_1d (1:nbpb) , frld , jpi, jpj, npb(1:nbpb) )
240	CALL tab_2d_1d( nbpb, h_ice_1d (1:nbpb) , hicif , jpi, jpj, npb(1:nbpb) )
241	CALL tab_2d_1d( nbpb, h_snow_1d (1:nbpb) , hsnif , jpi, jpj, npb(1:nbpb) )
242	CALL tab_2d_1d( nbpb, sist_1d (1:nbpb) , sist , jpi, jpj, npb(1:nbpb) )
243	CALL tab_2d_1d( nbpb, tbif_1d (1:nbpb , 1 ), tbif(:,:,1), jpi, jpj, npb(1:nbpb) )
244	CALL tab_2d_1d( nbpb, tbif_1d (1:nbpb , 2 ), tbif(:,:,2), jpi, jpj, npb(1:nbpb) )
245	CALL tab_2d_1d( nbpb, tbif_1d (1:nbpb , 3 ), tbif(:,:,3), jpi, jpj, npb(1:nbpb) )
246	CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb) , qsr_ice , jpi, jpj, npb(1:nbpb) )
247	CALL tab_2d_1d( nbpb, fr1_i0_1d (1:nbpb) , fr1_i0 , jpi, jpj, npb(1:nbpb) )
248	CALL tab_2d_1d( nbpb, fr2_i0_1d (1:nbpb) , fr2_i0 , jpi, jpj, npb(1:nbpb) )
249	CALL tab_2d_1d( nbpb, qnsr_ice_1d(1:nbpb) , qnsr_ice , jpi, jpj, npb(1:nbpb) )
250	CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb) , qla_ice , jpi, jpj, npb(1:nbpb) )
251	CALL tab_2d_1d( nbpb, dqla_ice_1d(1:nbpb) , dqla_ice , jpi, jpj, npb(1:nbpb) )
252	CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb) , dqns_ice , jpi, jpj, npb(1:nbpb) )
253	CALL tab_2d_1d( nbpb, tfu_1d (1:nbpb) , tfu , jpi, jpj, npb(1:nbpb) )
254	CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb) , sprecip , jpi, jpj, npb(1:nbpb) )
255	CALL tab_2d_1d( nbpb, fbif_1d (1:nbpb) , fbif , jpi, jpj, npb(1:nbpb) )
256	CALL tab_2d_1d( nbpb, thcm_1d (1:nbpb) , thcm , jpi, jpj, npb(1:nbpb) )
257	CALL tab_2d_1d( nbpb, qldif_1d (1:nbpb) , qldif , jpi, jpj, npb(1:nbpb) )
258	CALL tab_2d_1d( nbpb, qstbif_1d (1:nbpb) , qstoif , jpi, jpj, npb(1:nbpb) )
259	CALL tab_2d_1d( nbpb, rdmicif_1d (1:nbpb) , rdmicif , jpi, jpj, npb(1:nbpb) )
260	CALL tab_2d_1d( nbpb, dmgwi_1d (1:nbpb) , dmgwi , jpi, jpj, npb(1:nbpb) )
261	CALL tab_2d_1d( nbpb, qlbbq_1d (1:nbpb) , zqlbsbq , jpi, jpj, npb(1:nbpb) )
262
263
264
265	! call the ice growth routine.
266	CALL lim_thd_zdf( 1, nbpb )
267
268	! back to the geographic grid.
269	CALL tab_1d_2d( nbpb, frld , npb, frld_1d (1:nbpb) , jpi, jpj )
270	CALL tab_1d_2d( nbpb, hicif , npb, h_ice_1d (1:nbpb) , jpi, jpj )
271	CALL tab_1d_2d( nbpb, hsnif , npb, h_snow_1d (1:nbpb) , jpi, jpj )
272	CALL tab_1d_2d( nbpb, sist , npb, sist_1d (1:nbpb) , jpi, jpj )
273	CALL tab_1d_2d( nbpb, tbif(:,:,1), npb, tbif_1d (1:nbpb , 1 ), jpi, jpj )
274	CALL tab_1d_2d( nbpb, tbif(:,:,2), npb, tbif_1d (1:nbpb , 2 ), jpi, jpj )
275	CALL tab_1d_2d( nbpb, tbif(:,:,3), npb, tbif_1d (1:nbpb , 3 ), jpi, jpj )
276	CALL tab_1d_2d( nbpb, fscmbq , npb, fscbq_1d (1:nbpb) , jpi, jpj )
277	CALL tab_1d_2d( nbpb, ffltbif , npb, fltbif_1d (1:nbpb) , jpi, jpj )
278	CALL tab_1d_2d( nbpb, fstric , npb, fstbif_1d (1:nbpb) , jpi, jpj )
279	CALL tab_1d_2d( nbpb, qldif , npb, qldif_1d (1:nbpb) , jpi, jpj )
280	CALL tab_1d_2d( nbpb, qfvbq , npb, qfvbq_1d (1:nbpb) , jpi, jpj )
281	CALL tab_1d_2d( nbpb, qstoif , npb, qstbif_1d (1:nbpb) , jpi, jpj )
282	CALL tab_1d_2d( nbpb, rdmicif , npb, rdmicif_1d(1:nbpb) , jpi, jpj )
283	CALL tab_1d_2d( nbpb, dmgwi , npb, dmgwi_1d (1:nbpb) , jpi, jpj )
284	CALL tab_1d_2d( nbpb, rdmsnif , npb, rdmsnif_1d(1:nbpb) , jpi, jpj )
285	CALL tab_1d_2d( nbpb, rdvosif , npb, dvsbq_1d (1:nbpb) , jpi, jpj )
286	CALL tab_1d_2d( nbpb, rdvobif , npb, dvbbq_1d (1:nbpb) , jpi, jpj )
287	CALL tab_1d_2d( nbpb, fdvolif , npb, dvlbq_1d (1:nbpb) , jpi, jpj )
288	CALL tab_1d_2d( nbpb, rdvonif , npb, dvnbq_1d (1:nbpb) , jpi, jpj )
289
290
291	ENDIF
292
293
294	! Up-date sea ice thickness.
295	!---------------------------------
296	DO jj = 1, jpj
297	DO ji = 1, jpi
298	phicif(ji,jj) = hicif(ji,jj)
299	hicif(ji,jj) = hicif(ji,jj) * ( zone - MAX( zzero, SIGN( zone, - ( 1.0 - frld(ji,jj) ) ) ) )
300	END DO
301	END DO
302
303
304	! Tricky trick : add 2 to frld in the Southern Hemisphere.
305	!----------------------------------------------------------
306	DO jj = 1, jeqm1 !ibug in mpp
307	DO ji = 1, jpi
308	frld(ji,jj) = frld(ji,jj) + 2.0
309	END DO
310	END DO
311
312
313	! Select points for lateral accretion (this occurs when heat exchange
314	! between ice and ocean is negative; ocean losing heat)
315	!-----------------------------------------------------------------
316	nbpac = 0
317	DO jj = 1, jpj
318	DO ji = 1, jpi
319	!i yes! IF ( ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.0 ) THEN
320	IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.0 ) THEN
321	nbpac = nbpac + 1
322	npac( nbpac ) = (jj - 1) * jpi + ji
323	ENDIF
324	END DO
325	END DO
326
327	IF( l_ctl .AND. lwp ) THEN
328	WRITE(numout,*) 'lim_thd : phicif ', SUM( phicif ), ' hicif ', SUM( hicif )
329	WRITE(numout,*) 'lim_thd : nbpac = ', nbpac
330	ENDIF
331
332
333	!
334	! If ocean gains heat do nothing ; otherwise, one performs lateral accretion
335	!--------------------------------------------------------------------------------
336
337	IF ( nbpac > 0 ) THEN
338
339	!...Put the variable in a 1-D array for lateral accretion
340	CALL tab_2d_1d( nbpac, frld_1d (1:nbpac) , frld , jpi, jpj, npac(1:nbpac) )
341	CALL tab_2d_1d( nbpac, h_snow_1d (1:nbpac) , hsnif , jpi, jpj, npac(1:nbpac) )
342	CALL tab_2d_1d( nbpac, h_ice_1d (1:nbpac) , hicif , jpi, jpj, npac(1:nbpac) )
343	CALL tab_2d_1d( nbpac, tbif_1d (1:nbpac , 1 ), tbif(:,:,1), jpi, jpj, npac(1:nbpac) )
344	CALL tab_2d_1d( nbpac, tbif_1d (1:nbpac , 2 ), tbif(:,:,2), jpi, jpj, npac(1:nbpac) )
345	CALL tab_2d_1d( nbpac, tbif_1d (1:nbpac , 3 ), tbif(:,:,3), jpi, jpj, npac(1:nbpac) )
346	CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) )
347	CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) )
348	CALL tab_2d_1d( nbpac, qstbif_1d (1:nbpac) , qstoif , jpi, jpj, npac(1:nbpac) )
349	CALL tab_2d_1d( nbpac, rdmicif_1d(1:nbpac) , rdmicif , jpi, jpj, npac(1:nbpac) )
350	CALL tab_2d_1d( nbpac, dvlbq_1d (1:nbpac) , fdvolif , jpi, jpj, npac(1:nbpac) )
351	CALL tab_2d_1d( nbpac, tfu_1d (1:nbpac) , tfu , jpi, jpj, npac(1:nbpac) )
352
353	! call lateral accretion routine.
354	CALL lim_thd_lac( 1 , nbpac )
355
356	! back to the geographic grid
357	CALL tab_1d_2d( nbpac, frld , npac(1:nbpac), frld_1d (1:nbpac) , jpi, jpj )
358	CALL tab_1d_2d( nbpac, hsnif , npac(1:nbpac), h_snow_1d (1:nbpac) , jpi, jpj )
359	CALL tab_1d_2d( nbpac, hicif , npac(1:nbpac), h_ice_1d (1:nbpac) , jpi, jpj )
360	CALL tab_1d_2d( nbpac, tbif(:,:,1), npac(1:nbpac), tbif_1d (1:nbpac , 1 ), jpi, jpj )
361	CALL tab_1d_2d( nbpac, tbif(:,:,2), npac(1:nbpac), tbif_1d (1:nbpac , 2 ), jpi, jpj )
362	CALL tab_1d_2d( nbpac, tbif(:,:,3), npac(1:nbpac), tbif_1d (1:nbpac , 3 ), jpi, jpj )
363	CALL tab_1d_2d( nbpac, qstoif , npac(1:nbpac), qstbif_1d (1:nbpac) , jpi, jpj )
364	CALL tab_1d_2d( nbpac, rdmicif , npac(1:nbpac), rdmicif_1d(1:nbpac) , jpi, jpj )
365	CALL tab_1d_2d( nbpac, fdvolif , npac(1:nbpac), dvlbq_1d (1:nbpac) , jpi, jpj )
366
367
368	ENDIF
369
370
371	! Recover frld values between 0 and 1 in the Southern Hemisphere (tricky trick)
372	! Update daily thermodynamic ice production.
373	!------------------------------------------------------------------------------
374
375	DO jj = 1, jpj
376	DO ji = 1, jpi
377	frld (ji,jj) = MIN( frld(ji,jj), ABS( frld(ji,jj) - 2.0 ) )
378	hicifp(ji,jj) = hicif(ji,jj) * ( 1.0 - frld(ji,jj) ) - zhicifp(ji,jj) + hicifp(ji,jj)
379	END DO
380	END DO
381
382	IF( l_ctl .AND. lwp ) THEN
383	WRITE(numout,*) ' lim_thd end '
384	WRITE(numout,*) ' hicif ', SUM( hicif ), ' hsnif ', SUM( hsnif )
385	WRITE(numout,*) ' frld ', SUM( frld ), ' hicifp', SUM( hicifp )
386	WRITE(numout,*) ' phicif', SUM( phicif ), ' pfrld ', SUM( pfrld )
387	WRITE(numout,*) ' sist ', SUM( sist ), ' tbif 1', SUM( tbif (:,:,1) )
388	WRITE(numout,*) ' tbif 2', SUM( tbif(:,:,2) ), ' tbif 3', SUM( tbif (:,:,3) )
389	WRITE(numout,*) ' fdtcn ', SUM( fdtcn ), ' qdtcn ', SUM( qdtcn )
390	WRITE(numout,*) ' qstoif', SUM( qstoif ), ' fsbbq ', SUM( fsbbq )
391	ENDIF
392
393	END SUBROUTINE lim_thd
394
395	SUBROUTINE lim_thd_init
396	!!-------------------------------------------------------------------
397	!! * ROUTINE lim_thd_init *
398	!!
399	!! ** Purpose : Physical constants and parameters linked to the ice
400	!! thermodynamics
401	!!
402	!! ** Method : Read the namicethd namelist and check the ice-thermo
403	!! parameter values called at the first timestep (nit000)
404	!!
405	!! ** input : Namelist namicether
406	!!
407	!! history :
408	!! 8.5 ! 03-08 (C. Ethe) original code
409	!!-------------------------------------------------------------------
410	NAMELIST/namicethd/ hmelt , hiccrit, hicmin, hiclim, amax , &
411	& swiqst, sbeta , parlat, hakspl, hibspl, exld, &
412	& hakdif, hnzst , thth , parsub, alphs
413	!!-------------------------------------------------------------------
414
415
416	! Define the initial parameters
417	! -------------------------
418	REWIND( numnam_ice )
419	READ ( numnam_ice , namicethd )
420	IF (lwp) THEN
421	WRITE(numout,*)
422	WRITE(numout,*)'lim_thd_init : ice parameters for ice thermodynamic computation '
423	WRITE(numout,*)'~~~~~~~~~~~~'
424	WRITE(numout,*)' maximum melting at the bottom hmelt = ', hmelt
425	WRITE(numout,*)' ice thick. for lateral accretion in NH (SH) hiccrit(1/2) = ', hiccrit
426	WRITE(numout,*)' ice thick. corr. to max. energy stored in brine pocket hicmin = ', hicmin
427	WRITE(numout,*)' minimum ice thickness hiclim = ', hiclim
428	WRITE(numout,*)' maximum lead fraction amax = ', amax
429	WRITE(numout,*)' energy stored in brine pocket (=1) or not (=0) swiqst = ', swiqst
430	WRITE(numout,*)' numerical carac. of the scheme for diffusion in ice '
431	WRITE(numout,*)' Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) sbeta = ', sbeta
432	WRITE(numout,*)' percentage of energy used for lateral ablation parlat = ', parlat
433	WRITE(numout,*)' slope of distr. for Hakkinen-Mellor lateral melting hakspl = ', hakspl
434	WRITE(numout,*)' slope of distribution for Hibler lateral melting hibspl = ', hibspl
435	WRITE(numout,*)' exponent for leads-closure rate exld = ', exld
436	WRITE(numout,*)' coefficient for diffusions of ice and snow hakdif = ', hakdif
437	WRITE(numout,*)' threshold thick. for comp. of eq. thermal conductivity zhth = ', thth
438	WRITE(numout,*)' thickness of the surf. layer in temp. computation hnzst = ', hnzst
439	WRITE(numout,*)' switch for snow sublimation (=1) or not (=0) parsub = ', parsub
440	WRITE(numout,*)' coefficient for snow density when snow ice formation alphs = ', alphs
441	ENDIF
442
443	uscomi = 1.0 / ( 1.0 - amax ) ! inverse of minimum lead fraction
444	rcdsn = hakdif * rcdsn
445	rcdic = hakdif * rcdic
446
447	IF ( ( hsndif > 100.0 ) .OR. ( hicdif > 100.0 ) ) THEN
448	cnscg = 0.e0
449	ELSE
450	cnscg = rcpsn / rcpic ! ratio rcpsn/rcpic
451	ENDIF
452
453	END SUBROUTINE lim_thd_init
454
455	#else
456	!!======================================================================
457	!! * MODULE limthd *
458	!! No sea ice model
459	!!======================================================================
460	CONTAINS
461	SUBROUTINE lim_thd ! Empty routine
462	END SUBROUTINE lim_thd
463	#endif
464
465	!!======================================================================
466	END MODULE limthd

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

Download in other formats: