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_lac.F90 in branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 4332

Last change on this file since 4332 was 4332, checked in by clem, 10 years ago
update LIM3 to fix remaining bugs. Now working in global and regional config.
Property svn:keywords set to `Id`
File size: 33.2 KB

Line
1	MODULE limthd_lac
2	!!======================================================================
3	!! * MODULE limthd_lac *
4	!! lateral thermodynamic growth of the ice
5	!!======================================================================
6	!! History : LIM ! 2005-12 (M. Vancoppenolle) Original code
7	!! - ! 2006-01 (M. Vancoppenolle) add ITD
8	!! 3.0 ! 2007-07 (M. Vancoppenolle) Mass and energy conservation tested
9	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
10	!!----------------------------------------------------------------------
11	#if defined key_lim3
12	!!----------------------------------------------------------------------
13	!! 'key_lim3' LIM3 sea-ice model
14	!!----------------------------------------------------------------------
15	!! lim_lat_acr : lateral accretion of ice
16	!!----------------------------------------------------------------------
17	USE par_oce ! ocean parameters
18	USE dom_oce ! domain variables
19	USE phycst ! physical constants
20	USE sbc_oce ! Surface boundary condition: ocean fields
21	USE sbc_ice ! Surface boundary condition: ice fields
22	USE thd_ice ! LIM thermodynamics
23	USE dom_ice ! LIM domain
24	USE par_ice ! LIM parameters
25	USE ice ! LIM variables
26	USE limtab ! LIM 2D <==> 1D
27	USE limcons ! LIM conservation
28	USE in_out_manager ! I/O manager
29	USE lib_mpp ! MPP library
30	USE wrk_nemo ! work arrays
31	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
32
33	IMPLICIT NONE
34	PRIVATE
35
36	PUBLIC lim_thd_lac ! called by lim_thd
37
38	REAL(wp) :: epsi10 = 1.e-10_wp !
39	REAL(wp) :: zzero = 0._wp !
40	REAL(wp) :: zone = 1._wp !
41
42	!!----------------------------------------------------------------------
43	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
44	!! $Id$
45	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
46	!!----------------------------------------------------------------------
47	CONTAINS
48
49	SUBROUTINE lim_thd_lac
50	!!-------------------------------------------------------------------
51	!! * ROUTINE lim_thd_lac *
52	!!
53	!! ** Purpose : Computation of the evolution of the ice thickness and
54	!! concentration as a function of the heat balance in the leads.
55	!! It is only used for lateral accretion
56	!!
57	!! ** Method : Ice is formed in the open water when ocean lose heat
58	!! (heat budget of open water Bl is negative) .
59	!! Computation of the increase of 1-A (ice concentration) fol-
60	!! lowing the law :
61	!! (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
62	!! where - h0 is the thickness of ice created in the lead
63	!! - a is a minimum fraction for leads
64	!! - F is a monotonic non-increasing function defined as:
65	!! F(X)=( 1 - Xexld )(1.0/exld)
66	!! - exld is the exponent closure rate (=2 default val.)
67	!!
68	!! ** Action : - Adjustment of snow and ice thicknesses and heat
69	!! content in brine pockets
70	!! - Updating ice internal temperature
71	!! - Computation of variation of ice volume and mass
72	!! - Computation of frldb after lateral accretion and
73	!! update ht_s_b, ht_i_b and tbif_1d(:,:)
74	!!------------------------------------------------------------------------
75	INTEGER :: ji,jj,jk,jl,jm ! dummy loop indices
76	INTEGER :: layer, nbpac ! local integers
77	INTEGER :: ii, ij, iter ! - -
78	REAL(wp) :: ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zinda, zde ! local scalars
79	REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - -
80	REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - -
81	LOGICAL :: iterate_frazil ! iterate frazil ice collection thickness
82	CHARACTER (len = 15) :: fieldid
83	!
84	INTEGER , POINTER, DIMENSION(:) :: zcatac ! indexes of categories where new ice grows
85	REAL(wp), POINTER, DIMENSION(:) :: zswinew ! switch for new ice or not
86
87	REAL(wp), POINTER, DIMENSION(:) :: zv_newice ! volume of accreted ice
88	REAL(wp), POINTER, DIMENSION(:) :: za_newice ! fractional area of accreted ice
89	REAL(wp), POINTER, DIMENSION(:) :: zh_newice ! thickness of accreted ice
90	REAL(wp), POINTER, DIMENSION(:) :: ze_newice ! heat content of accreted ice
91	REAL(wp), POINTER, DIMENSION(:) :: zs_newice ! salinity of accreted ice
92	REAL(wp), POINTER, DIMENSION(:) :: zo_newice ! age of accreted ice
93	REAL(wp), POINTER, DIMENSION(:) :: zdv_res ! residual volume in case of excessive heat budget
94	REAL(wp), POINTER, DIMENSION(:) :: zda_res ! residual area in case of excessive heat budget
95	REAL(wp), POINTER, DIMENSION(:) :: zat_i_ac ! total ice fraction
96	REAL(wp), POINTER, DIMENSION(:) :: zat_i_lev ! total ice fraction for level ice only (type 1)
97	REAL(wp), POINTER, DIMENSION(:) :: zdh_frazb ! accretion of frazil ice at the ice bottom
98	REAL(wp), POINTER, DIMENSION(:) :: zvrel_ac ! relative ice / frazil velocity (1D vector)
99
100	REAL(wp), POINTER, DIMENSION(:,:) :: zhice_old ! previous ice thickness
101	REAL(wp), POINTER, DIMENSION(:,:) :: zdummy ! dummy thickness of new ice
102	REAL(wp), POINTER, DIMENSION(:,:) :: zdhicbot ! thickness of new ice which is accreted vertically
103	REAL(wp), POINTER, DIMENSION(:,:) :: zv_old ! old volume of ice in category jl
104	REAL(wp), POINTER, DIMENSION(:,:) :: za_old ! old area of ice in category jl
105	REAL(wp), POINTER, DIMENSION(:,:) :: za_i_ac ! 1-D version of a_i
106	REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_ac ! 1-D version of v_i
107	REAL(wp), POINTER, DIMENSION(:,:) :: zoa_i_ac ! 1-D version of oa_i
108	REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_ac ! 1-D version of smv_i
109
110	REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_ac !: 1-D version of e_i
111
112	REAL(wp), POINTER, DIMENSION(:) :: zqbgow ! heat budget of the open water (negative)
113	REAL(wp), POINTER, DIMENSION(:) :: zdhex ! excessively thick accreted sea ice (hlead-hice)
114
115	REAL(wp), POINTER, DIMENSION(:,:,:) :: zqm0 ! old layer-system heat content
116	REAL(wp), POINTER, DIMENSION(:,:,:) :: zthick0 ! old ice thickness
117
118	REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories
119	REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories
120	REAL(wp), POINTER, DIMENSION(:,:) :: et_i_init, et_i_final ! ice energy summed over categories
121	REAL(wp), POINTER, DIMENSION(:,:) :: et_s_init ! snow energy summed over categories
122	REAL(wp), POINTER, DIMENSION(:,:) :: zvrel ! relative ice / frazil velocity
123	!!-----------------------------------------------------------------------!
124
125	CALL wrk_alloc( jpij, zcatac ) ! integer
126	CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
127	CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
128	CALL wrk_alloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )
129	CALL wrk_alloc( jpij,jkmax,jpl, ze_i_ac )
130	CALL wrk_alloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
131	CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel )
132
133	et_i_init(:,:) = 0._wp
134	et_s_init(:,:) = 0._wp
135	vt_i_init(:,:) = 0._wp
136	vt_s_init(:,:) = 0._wp
137
138	!------------------------------------------------------------------------------!
139	! 1) Conservation check and changes in each ice category
140	!------------------------------------------------------------------------------!
141	IF( con_i ) THEN
142	CALL lim_column_sum ( jpl, v_i , vt_i_init)
143	CALL lim_column_sum ( jpl, v_s , vt_s_init)
144	CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init)
145	CALL lim_column_sum ( jpl, e_s(:,:,1,:) , et_s_init)
146	ENDIF
147
148	!------------------------------------------------------------------------------\|
149	! 2) Convert units for ice internal energy
150	!------------------------------------------------------------------------------\|
151	DO jl = 1, jpl
152	DO jk = 1, nlay_i
153	DO jj = 1, jpj
154	DO ji = 1, jpi
155	!Energy of melting q(S,T) [J.m-3]
156	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * REAL( nlay_i )
157	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes
158	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb
159	END DO
160	END DO
161	END DO
162	END DO
163
164	!------------------------------------------------------------------------------!
165	! 3) Collection thickness of ice formed in leads and polynyas
166	!------------------------------------------------------------------------------!
167	! hicol is the thickness of new ice formed in open water
168	! hicol can be either prescribed (frazswi = 0)
169	! or computed (frazswi = 1)
170	! Frazil ice forms in open water, is transported by wind
171	! accumulates at the edge of the consolidated ice edge
172	! where it forms aggregates of a specific thickness called
173	! collection thickness.
174
175	! Note : the following algorithm currently breaks vectorization
176	!
177
178	zvrel(:,:) = 0._wp
179
180	! Default new ice thickness
181	hicol(:,:) = hiccrit(1)
182
183	IF( fraz_swi == 1._wp ) THEN
184
185	!--------------------
186	! Physical constants
187	!--------------------
188	hicol(:,:) = 0._wp
189
190	zhicrit = 0.04 ! frazil ice thickness
191	ztwogp = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav
192	zsqcd = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)
193	zgamafr = 0.03
194
195	DO jj = 1, jpj
196	DO ji = 1, jpi
197
198	IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
199	!-------------
200	! Wind stress
201	!-------------
202	! C-grid wind stress components
203	ztaux = ( utau_ice(ji-1,jj ) * tmu(ji-1,jj ) &
204	& + utau_ice(ji ,jj ) * tmu(ji ,jj ) ) * 0.5_wp
205	ztauy = ( vtau_ice(ji ,jj-1) * tmv(ji ,jj-1) &
206	& + vtau_ice(ji ,jj ) * tmv(ji ,jj ) ) * 0.5_wp
207	! Square root of wind stress
208	ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
209
210	!---------------------
211	! Frazil ice velocity
212	!---------------------
213	zvfrx = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10)
214	zvfry = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10)
215
216	!-------------------
217	! Pack ice velocity
218	!-------------------
219	! C-grid ice velocity
220	zindb = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) ) )
221	zvgx = zindb * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) &
222	& + u_ice(ji,jj ) * tmu(ji ,jj ) ) * 0.5_wp
223	zvgy = zindb * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) &
224	& + v_ice(ji,jj ) * tmv(ji ,jj ) ) * 0.5_wp
225
226	!-----------------------------------
227	! Relative frazil/pack ice velocity
228	!-----------------------------------
229	! absolute relative velocity
230	zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) &
231	& + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 )
232	zvrel(ji,jj) = SQRT( zvrel2 )
233
234	!---------------------
235	! Iterative procedure
236	!---------------------
237	hicol(ji,jj) = zhicrit + 0.1
238	hicol(ji,jj) = zhicrit + hicol(ji,jj) &
239	& / ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) * ztwogp * zvrel2
240
241	!!gm better coding: above: hicol(ji,jj) * hicol(ji,jj) = (zhicrit + 0.1)*(zhicrit + 0.1)
242	!!gm = zhicrit*2 + 0.2zhicrit +0.01
243	!!gm therefore the 2 lines with hicol can be replaced by 1 line:
244	!!gm hicol(ji,jj) = zhicrit + (zhicrit + 0.1) / ( 0.2 * zhicrit + 0.01 ) * ztwogp * zvrel2
245	!!gm further more (zhicrit + 0.1)/(0.2 * zhicrit + 0.01 )*ztwogp can be computed one for all outside the DO loop
246
247	iter = 1
248	iterate_frazil = .true.
249
250	DO WHILE ( iter .LT. 100 .AND. iterate_frazil )
251	zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)2 - zhicrit2 ) &
252	- hicol(ji,jj) * zhicrit * ztwogp * zvrel2
253	zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0*hicol(ji,jj) + zhicrit ) &
254	- zhicrit * ztwogp * zvrel2
255	zhicol_new = hicol(ji,jj) - zf/zfp
256	hicol(ji,jj) = zhicol_new
257
258	iter = iter + 1
259
260	END DO ! do while
261
262	ENDIF ! end of selection of pixels where ice forms
263
264	END DO ! loop on ji ends
265	END DO ! loop on jj ends
266
267	ENDIF ! End of computation of frazil ice collection thickness
268
269	!------------------------------------------------------------------------------!
270	! 4) Identify grid points where new ice forms
271	!------------------------------------------------------------------------------!
272
273	!-------------------------------------
274	! Select points for new ice formation
275	!-------------------------------------
276	! This occurs if open water energy budget is negative
277	nbpac = 0
278	DO jj = 1, jpj
279	DO ji = 1, jpi
280	IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0._wp ) THEN
281	nbpac = nbpac + 1
282	npac( nbpac ) = (jj - 1) * jpi + ji
283	ENDIF
284	END DO
285	END DO
286
287	! debug point to follow
288	jiindex_1d = 0
289	IF( ln_nicep ) THEN
290	DO ji = mi0(jiindx), mi1(jiindx)
291	DO jj = mj0(jjindx), mj1(jjindx)
292	IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0._wp ) THEN
293	jiindex_1d = (jj - 1) * jpi + ji
294	ENDIF
295	END DO
296	END DO
297	ENDIF
298
299	IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
300
301	!------------------------------
302	! Move from 2-D to 1-D vectors
303	!------------------------------
304	! If ocean gains heat do nothing
305	! 0therwise compute new ice formation
306
307	IF ( nbpac > 0 ) THEN
308
309	CALL tab_2d_1d( nbpac, zat_i_ac (1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) )
310	DO jl = 1, jpl
311	CALL tab_2d_1d( nbpac, za_i_ac (1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
312	CALL tab_2d_1d( nbpac, zv_i_ac (1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
313	CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
314	CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
315	DO jk = 1, nlay_i
316	CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
317	END DO ! jk
318	END DO ! jl
319
320	CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) )
321	CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) )
322	CALL tab_2d_1d( nbpac, t_bo_b (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) )
323	CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac) , sfx_thd, jpi, jpj, npac(1:nbpac) )
324	CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac) , rdm_ice, jpi, jpj, npac(1:nbpac) )
325	CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) )
326	CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) )
327
328	!------------------------------------------------------------------------------!
329	! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
330	!------------------------------------------------------------------------------!
331
332	!----------------------
333	! Thickness of new ice
334	!----------------------
335	DO ji = 1, nbpac
336	zh_newice(ji) = hiccrit(1)
337	END DO
338	IF( fraz_swi == 1.0 ) zh_newice(:) = hicol_b(:)
339
340	!----------------------
341	! Salinity of new ice
342	!----------------------
343
344	SELECT CASE ( num_sal )
345	CASE ( 1 ) ! Sice = constant
346	zs_newice(:) = bulk_sal
347	CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)]
348	DO ji = 1, nbpac
349	ii = MOD( npac(ji) - 1 , jpi ) + 1
350	ij = ( npac(ji) - 1 ) / jpi + 1
351	zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij) )
352	END DO
353	CASE ( 3 ) ! Sice = F(z) [multiyear ice]
354	zs_newice(:) = 2.3
355	END SELECT
356
357
358	!-------------------------
359	! Heat content of new ice
360	!-------------------------
361	! We assume that new ice is formed at the seawater freezing point
362	DO ji = 1, nbpac
363	ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K)
364	ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) &
365	& + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) ) &
366	& - rcp * ( ztmelts - rtt ) )
367	ze_newice(ji) = MAX( ze_newice(ji) , 0._wp ) &
368	& + MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) ) * rhoic * lfus
369	END DO ! ji
370	!----------------
371	! Age of new ice
372	!----------------
373	DO ji = 1, nbpac
374	zo_newice(ji) = 0._wp
375	END DO ! ji
376
377	!--------------------------
378	! Open water energy budget
379	!--------------------------
380	DO ji = 1, nbpac
381	zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) !<0
382	END DO ! ji
383
384	!-------------------
385	! Volume of new ice
386	!-------------------
387	DO ji = 1, nbpac
388	zv_newice(ji) = - zqbgow(ji) / ze_newice(ji)
389
390	! A fraction zfrazb of frazil ice is accreted at the ice bottom
391	zfrazb = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
392	zdh_frazb(ji) = zfrazb * zv_newice(ji)
393	zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
394	END DO
395
396	!------------------------------------
397	! Diags for energy conservation test
398	!------------------------------------
399	DO ji = 1, nbpac
400	ii = MOD( npac(ji) - 1 , jpi ) + 1
401	ij = ( npac(ji) - 1 ) / jpi + 1
402	!
403	zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji)
404	!
405	vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji) ! volume
406	et_i_init(ii,ij) = et_i_init(ii,ij) + zde ! Energy
407
408	END DO
409
410	! keep new ice volume in memory
411	CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj )
412
413	!-----------------
414	! Area of new ice
415	!-----------------
416	DO ji = 1, nbpac
417	ii = MOD( npac(ji) - 1 , jpi ) + 1
418	ij = ( npac(ji) - 1 ) / jpi + 1
419	za_newice(ji) = zv_newice(ji) / zh_newice(ji)
420	diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem
421	END DO !ji
422
423	!------------------------------------------------------------------------------!
424	! 6) Redistribute new ice area and volume into ice categories !
425	!------------------------------------------------------------------------------!
426
427	!-----------------------------------------
428	! Keep old ice areas and volume in memory
429	!-----------------------------------------
430	zv_old(:,:) = zv_i_ac(:,:)
431	za_old(:,:) = za_i_ac(:,:)
432
433	!-------------------------------------------
434	! Compute excessive new ice area and volume
435	!-------------------------------------------
436	! If lateral ice growth gives an ice concentration gt 1, then
437	! we keep the excessive volume in memory and attribute it later to bottom accretion
438	DO ji = 1, nbpac
439	IF ( za_newice(ji) > ( amax - zat_i_ac(ji) ) ) THEN
440	zda_res(ji) = za_newice(ji) - ( amax - zat_i_ac(ji) )
441	zdv_res(ji) = zda_res (ji) * zh_newice(ji)
442	za_newice(ji) = za_newice(ji) - zda_res (ji)
443	zv_newice(ji) = zv_newice(ji) - zdv_res (ji)
444	ELSE
445	zda_res(ji) = 0._wp
446	zdv_res(ji) = 0._wp
447	ENDIF
448	END DO ! ji
449
450	!------------------------------------------------
451	! Laterally redistribute new ice volume and area
452	!------------------------------------------------
453	zat_i_ac(:) = 0._wp
454	DO jl = 1, jpl
455	DO ji = 1, nbpac
456	IF( hi_max (jl-1) < zh_newice(ji) .AND. &
457	& zh_newice(ji) <= hi_max (jl) ) THEN
458	za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji)
459	zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji)
460	zat_i_ac(ji) = zat_i_ac(ji) + za_i_ac (ji,jl)
461	zcatac (ji) = jl
462	ENDIF
463	END DO
464	END DO
465
466	!----------------------------------
467	! Heat content - lateral accretion
468	!----------------------------------
469	DO ji = 1, nbpac
470	jl = zcatac(ji) ! categroy in which new ice is put
471	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) ) ! zindb=1 if ice =0 otherwise
472	zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb ! old ice thickness
473	zdhex (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) ) ! difference in thickness
474	zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) ) ! ice totally new in jl category
475	END DO
476
477	DO jk = 1, nlay_i
478	DO ji = 1, nbpac
479	jl = zcatac(ji)
480	zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ]
481	zalphai = MIN( zhice_old(ji,jl) * REAL( jk ) / REAL( nlay_i ), zh_newice(ji) ) &
482	& - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) )
483	ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji) &
484	+ ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / REAL( nlay_i ) &
485	+ za_newice(ji) * ze_newice(ji) * zalphai &
486	+ za_newice(ji) * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) )
487	END DO
488	END DO
489
490	!-----------------------------------------------
491	! Add excessive volume of new ice at the bottom
492	!-----------------------------------------------
493	! If the ice concentration exceeds 1, the remaining volume of new ice
494	! is equally redistributed among all ice categories in which there is
495	! ice
496
497	! Fraction of level ice
498	jm = 1
499	zat_i_lev(:) = 0._wp
500
501	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
502	DO ji = 1, nbpac
503	zat_i_lev(ji) = zat_i_lev(ji) + za_i_ac(ji,jl)
504	END DO
505	END DO
506
507	IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
508	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
509	DO ji = 1, nbpac
510	zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) )
511	zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi10 ) ) ! clem
512	zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi10 )
513	END DO
514	END DO
515	IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
516
517	!---------------------------------
518	! Heat content - bottom accretion
519	!---------------------------------
520	jm = 1
521	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
522	DO ji = 1, nbpac
523	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) ) ! zindb=1 if ice =0 otherwise
524	zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
525	zdhicbot (ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb &
526	& + zindb * zdh_frazb(ji) ! frazil ice may coalesce
527	zdummy(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb ! thickness of residual ice
528	END DO
529	END DO
530
531	! old layers thicknesses and enthalpies
532	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
533	DO jk = 1, nlay_i
534	DO ji = 1, nbpac
535	zthick0(ji,jk,jl) = zhice_old(ji,jl) / REAL( nlay_i )
536	zqm0 (ji,jk,jl) = ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)
537	END DO
538	END DO
539	END DO
540	!!gm ??? why the previous do loop if ocerwriten by the following one ?
541	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
542	DO ji = 1, nbpac
543	zthick0(ji,nlay_i+1,jl) = zdhicbot(ji,jl)
544	zqm0 (ji,nlay_i+1,jl) = ze_newice(ji) * zdhicbot(ji,jl)
545	END DO ! ji
546	END DO ! jl
547
548	! Redistributing energy on the new grid
549	ze_i_ac(:,:,:) = 0._wp
550	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
551	DO jk = 1, nlay_i
552	DO layer = 1, nlay_i + 1
553	DO ji = 1, nbpac
554	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) )
555	! Redistributing energy on the new grid
556	zweight = MAX ( MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) ) &
557	& - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp ) &
558	& /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb
559	ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl)
560	END DO ! ji
561	END DO ! layer
562	END DO ! jk
563	END DO ! jl
564
565	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
566	DO jk = 1, nlay_i
567	DO ji = 1, nbpac
568	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )
569	ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) &
570	& / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb
571	END DO
572	END DO
573	END DO
574
575	!------------
576	! Update age
577	!------------
578	DO jl = 1, jpl
579	DO ji = 1, nbpac
580	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes
581	zoa_i_ac(ji,jl) = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
582	END DO
583	END DO
584
585	!-----------------
586	! Update salinity
587	!-----------------
588	!clem IF( num_sal == 2 ) THEN
589	DO jl = 1, jpl
590	DO ji = 1, nbpac
591	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes
592	zdv = zv_i_ac(ji,jl) - zv_old(ji,jl)
593	zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif
594	END DO
595	END DO
596	!clem ENDIF
597
598	!--------------------------------
599	! Update mass/salt fluxes (clem)
600	!--------------------------------
601	DO jl = 1, jpl
602	DO ji = 1, nbpac
603	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes
604	zdv = zv_i_ac(ji,jl) - zv_old(ji,jl)
605	rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb
606	sfx_thd_1d(ji) = sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb
607	END DO
608	END DO
609
610	!------------------------------------------------------------------------------!
611	! 8) Change 2D vectors to 1D vectors
612	!------------------------------------------------------------------------------!
613	DO jl = 1, jpl
614	CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ac (1:nbpac,jl), jpi, jpj )
615	CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj )
616	CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj )
617	!clem IF ( num_sal == 2 ) &
618	CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj )
619	DO jk = 1, nlay_i
620	CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj )
621	END DO
622	END DO
623	CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj )
624	CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj )
625	!
626	ENDIF ! nbpac > 0
627
628	!------------------------------------------------------------------------------!
629	! 9) Change units for e_i
630	!------------------------------------------------------------------------------!
631	DO jl = 1, jpl
632	DO jk = 1, nlay_i ! heat content in 10^9 Joules
633	e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i ) / unit_fac
634	END DO
635	END DO
636
637	!------------------------------------------------------------------------------\|
638	! 10) Conservation check and changes in each ice category
639	!------------------------------------------------------------------------------\|
640	IF( con_i ) THEN
641	CALL lim_column_sum (jpl, v_i, vt_i_final)
642	fieldid = 'v_i, limthd_lac'
643	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
644	!
645	CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)
646	fieldid = 'e_i, limthd_lac'
647	CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid)
648	!
649	CALL lim_column_sum (jpl, v_s, vt_s_final)
650	fieldid = 'v_s, limthd_lac'
651	CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)
652	!
653	! CALL lim_column_sum (jpl, e_s(:,:,1,:) , et_s_init)
654	! fieldid = 'e_s, limthd_lac'
655	! CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)
656	IF( ln_nicep ) THEN
657	DO ji = mi0(jiindx), mi1(jiindx)
658	DO jj = mj0(jjindx), mj1(jjindx)
659	WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj)
660	WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj)
661	WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj)
662	WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj)
663	END DO
664	END DO
665	ENDIF
666	!
667	ENDIF
668	!
669	CALL wrk_dealloc( jpij, zcatac ) ! integer
670	CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
671	CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
672	CALL wrk_dealloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )
673	CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac )
674	CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
675	CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel )
676	!
677	END SUBROUTINE lim_thd_lac
678
679	#else
680	!!----------------------------------------------------------------------
681	!! Default option NO LIM3 sea-ice model
682	!!----------------------------------------------------------------------
683	CONTAINS
684	SUBROUTINE lim_thd_lac ! Empty routine
685	END SUBROUTINE lim_thd_lac
686	#endif
687
688	!!======================================================================
689	END MODULE limthd_lac

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

Download in other formats: