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 trunk/NEMO/LIM_SRC_2 – NEMO

Context Navigation

source: trunk/NEMO/LIM_SRC_2/limthd_2.F90 @ 1649

Last change on this file since 1649 was 1565, checked in by rblod, 15 years ago
Use appropiate e3 for divergence computation and for sea-ice, see ticket #507
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 28.1 KB

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

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

Download in other formats: