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_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 4291

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

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

Download in other formats: