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_2.F90 in branches/dev_002_LIM/NEMO/LIM_SRC_2 – NEMO

Context Navigation

source: branches/dev_002_LIM/NEMO/LIM_SRC_2/limthd_2.F90 @ 823

Last change on this file since 823 was 823, checked in by rblod, 16 years ago
Final step to rename LIM_SRC in LIM_SRC_2
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 24.9 KB

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

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

Download in other formats: