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

Context Navigation

source: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 4506

Last change on this file since 4506 was 4506, checked in by vancop, 10 years ago
[Heat conservation in LIM3, part HC1 (LIM_SRC_3_HC17)]
Property svn:keywords set to `Id`
File size: 34.6 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	REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
85	REAL(wp) :: zEi ! sea ice specific enthalpy (J/kg)
86	REAL(wp) :: zEw ! seawater specific enthalpy (J/kg)
87	REAL(wp) :: zfmdt ! mass flux x time step (kg/m2, >0 towards ocean)
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) + epsi10 ) ) !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	ENDIF
289	END DO
290	END DO
291
292	! debug point to follow
293	jiindex_1d = 0
294	IF( ln_nicep ) THEN
295	DO ji = mi0(jiindx), mi1(jiindx)
296	DO jj = mj0(jjindx), mj1(jjindx)
297	IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0._wp ) THEN
298	jiindex_1d = (jj - 1) * jpi + ji
299	ENDIF
300	END DO
301	END DO
302	ENDIF
303
304	IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
305
306	!------------------------------
307	! Move from 2-D to 1-D vectors
308	!------------------------------
309	! If ocean gains heat do nothing
310	! 0therwise compute new ice formation
311
312	IF ( nbpac > 0 ) THEN
313
314	CALL tab_2d_1d( nbpac, zat_i_ac (1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) )
315	DO jl = 1, jpl
316	CALL tab_2d_1d( nbpac, za_i_ac (1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
317	CALL tab_2d_1d( nbpac, zv_i_ac (1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
318	CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
319	CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
320	DO jk = 1, nlay_i
321	CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
322	END DO ! jk
323	END DO ! jl
324
325	CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) )
326	CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) )
327	CALL tab_2d_1d( nbpac, t_bo_b (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) )
328	CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac) , sfx_thd, jpi, jpj, npac(1:nbpac) )
329	CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac) , rdm_ice, jpi, jpj, npac(1:nbpac) )
330	CALL tab_2d_1d( nbpac, rdq_ice_1d(1:nbpac) , rdq_ice, jpi, jpj, npac(1:nbpac) )
331	CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) )
332	CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) )
333
334	!------------------------------------------------------------------------------!
335	! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
336	!------------------------------------------------------------------------------!
337
338	!----------------------
339	! Thickness of new ice
340	!----------------------
341	DO ji = 1, nbpac
342	zh_newice(ji) = hiccrit(1)
343	END DO
344	IF( fraz_swi == 1.0 ) zh_newice(:) = hicol_b(:)
345
346	!----------------------
347	! Salinity of new ice
348	!----------------------
349
350	SELECT CASE ( num_sal )
351	CASE ( 1 ) ! Sice = constant
352	zs_newice(:) = bulk_sal
353	CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)]
354	DO ji = 1, nbpac
355	ii = MOD( npac(ji) - 1 , jpi ) + 1
356	ij = ( npac(ji) - 1 ) / jpi + 1
357	zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij) )
358	END DO
359	CASE ( 3 ) ! Sice = F(z) [multiyear ice]
360	zs_newice(:) = 2.3
361	END SELECT
362
363
364	!-------------------------
365	! Heat content of new ice
366	!-------------------------
367	! We assume that new ice is formed at the seawater freezing point
368	DO ji = 1, nbpac
369	ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K)
370	ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) &
371	& + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) ) &
372	& - rcp * ( ztmelts - rtt ) )
373	! MV HC 2014 comment I dont see why this line below is here... ?
374	! This implies that ze_newice gets to rhoic*Lfus if it was negative, but this should never happen
375	ze_newice(ji) = MAX( ze_newice(ji) , 0._wp ) &
376	& + MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) ) * rhoic * lfus
377	END DO ! ji
378	!----------------
379	! Age of new ice
380	!----------------
381	DO ji = 1, nbpac
382	zo_newice(ji) = 0._wp
383	END DO ! ji
384
385	!--------------------------
386	! Open water energy budget
387	!--------------------------
388	DO ji = 1, nbpac
389	zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) !<0
390	END DO ! ji
391
392	!-------------------
393	! Volume of new ice
394	!-------------------
395	DO ji = 1, nbpac
396
397	zEi = - ze_newice(ji) / rhoic ! specific enthalpy of forming ice [J/kg]
398
399	zEw = rcp * ( t_bo_b(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_b [J/kg]
400
401	zdE = zEi - zEw ! specific enthalpy difference [J/kg]
402
403	zfmdt = - zqbgow(ji) / zdE ! Fm.dt [kg/m2] (<0)
404
405	zv_newice(ji) = - zfmdt / rhoic
406
407	zQm = zfmdt * zEw ! heat to the ocean >0 associated with mass flux
408
409	! Contribution to energy flux to the ocean [J/m2], >0
410	rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm
411
412	! A fraction zfrazb of frazil ice is accreted at the ice bottom
413	zfrazb = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
414	zdh_frazb(ji) = zfrazb * zv_newice(ji)
415	zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
416	END DO
417
418	!------------------------------------
419	! Diags for energy conservation test
420	!------------------------------------
421	DO ji = 1, nbpac
422	ii = MOD( npac(ji) - 1 , jpi ) + 1
423	ij = ( npac(ji) - 1 ) / jpi + 1
424	!
425	zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji)
426	!
427	vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji) ! volume
428	et_i_init(ii,ij) = et_i_init(ii,ij) + zde ! Energy
429
430	END DO
431
432	! keep new ice volume in memory
433	CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj )
434
435	!-----------------
436	! Area of new ice
437	!-----------------
438	DO ji = 1, nbpac
439	ii = MOD( npac(ji) - 1 , jpi ) + 1
440	ij = ( npac(ji) - 1 ) / jpi + 1
441	za_newice(ji) = zv_newice(ji) / zh_newice(ji)
442	diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem
443	END DO !ji
444
445	!------------------------------------------------------------------------------!
446	! 6) Redistribute new ice area and volume into ice categories !
447	!------------------------------------------------------------------------------!
448
449	!-----------------------------------------
450	! Keep old ice areas and volume in memory
451	!-----------------------------------------
452	zv_old(:,:) = zv_i_ac(:,:)
453	za_old(:,:) = za_i_ac(:,:)
454
455	!-------------------------------------------
456	! Compute excessive new ice area and volume
457	!-------------------------------------------
458	! If lateral ice growth gives an ice concentration gt 1, then
459	! we keep the excessive volume in memory and attribute it later to bottom accretion
460	DO ji = 1, nbpac
461	IF ( za_newice(ji) > ( amax - zat_i_ac(ji) ) ) THEN
462	zda_res(ji) = za_newice(ji) - ( amax - zat_i_ac(ji) )
463	zdv_res(ji) = zda_res (ji) * zh_newice(ji)
464	za_newice(ji) = za_newice(ji) - zda_res (ji)
465	zv_newice(ji) = zv_newice(ji) - zdv_res (ji)
466	ELSE
467	zda_res(ji) = 0._wp
468	zdv_res(ji) = 0._wp
469	ENDIF
470	END DO ! ji
471
472	!------------------------------------------------
473	! Laterally redistribute new ice volume and area
474	!------------------------------------------------
475	zat_i_ac(:) = 0._wp
476	DO jl = 1, jpl
477	DO ji = 1, nbpac
478	IF( hi_max (jl-1) < zh_newice(ji) .AND. &
479	& zh_newice(ji) <= hi_max (jl) ) THEN
480	za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji)
481	zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji)
482	zat_i_ac(ji) = zat_i_ac(ji) + za_i_ac (ji,jl)
483	zcatac (ji) = jl
484	ENDIF
485	END DO
486	END DO
487
488	!----------------------------------
489	! Heat content - lateral accretion
490	!----------------------------------
491	DO ji = 1, nbpac
492	jl = zcatac(ji) ! categroy in which new ice is put
493	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) ) ! zindb=1 if ice =0 otherwise
494	zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb ! old ice thickness
495	zdhex (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) ) ! difference in thickness
496	zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) ) ! ice totally new in jl category
497	END DO
498
499	DO jk = 1, nlay_i
500	DO ji = 1, nbpac
501	jl = zcatac(ji)
502	zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ]
503	zalphai = MIN( zhice_old(ji,jl) * REAL( jk ) / REAL( nlay_i ), zh_newice(ji) ) &
504	& - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) )
505	ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji) &
506	+ ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / REAL( nlay_i ) &
507	+ za_newice(ji) * ze_newice(ji) * zalphai &
508	+ za_newice(ji) * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) )
509	END DO
510	END DO
511
512	!-----------------------------------------------
513	! Add excessive volume of new ice at the bottom
514	!-----------------------------------------------
515	! If the ice concentration exceeds 1, the remaining volume of new ice
516	! is equally redistributed among all ice categories in which there is
517	! ice
518
519	! Fraction of level ice
520	jm = 1
521	zat_i_lev(:) = 0._wp
522
523	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
524	DO ji = 1, nbpac
525	zat_i_lev(ji) = zat_i_lev(ji) + za_i_ac(ji,jl)
526	END DO
527	END DO
528
529	IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
530	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
531	DO ji = 1, nbpac
532	zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) )
533	zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi10 ) ) ! clem
534	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 )
535	END DO
536	END DO
537	IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
538
539	!---------------------------------
540	! Heat content - bottom accretion
541	!---------------------------------
542	jm = 1
543	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
544	DO ji = 1, nbpac
545	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) ) ! zindb=1 if ice =0 otherwise
546	zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
547	zdhicbot (ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb &
548	& + zindb * zdh_frazb(ji) ! frazil ice may coalesce
549	zdummy(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb ! thickness of residual ice
550	END DO
551	END DO
552
553	! old layers thicknesses and enthalpies
554	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
555	DO jk = 1, nlay_i
556	DO ji = 1, nbpac
557	zthick0(ji,jk,jl) = zhice_old(ji,jl) / REAL( nlay_i )
558	zqm0 (ji,jk,jl) = ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)
559	END DO
560	END DO
561	END DO
562	!!gm ??? why the previous do loop if ocerwriten by the following one ?
563	!! MV HC 2014 just because it'snot the same index and the expression is different
564	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
565	DO ji = 1, nbpac
566	zthick0(ji,nlay_i+1,jl) = zdhicbot(ji,jl)
567	zqm0 (ji,nlay_i+1,jl) = ze_newice(ji) * zdhicbot(ji,jl)
568	END DO ! ji
569	END DO ! jl
570
571	! Redistributing energy on the new grid
572	ze_i_ac(:,:,:) = 0._wp
573	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
574	DO jk = 1, nlay_i
575	DO layer = 1, nlay_i + 1
576	DO ji = 1, nbpac
577	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) )
578	! Redistributing energy on the new grid
579	zweight = MAX ( MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) ) &
580	& - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp ) &
581	& /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb
582	ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl)
583	END DO ! ji
584	END DO ! layer
585	END DO ! jk
586	END DO ! jl
587
588	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
589	DO jk = 1, nlay_i
590	DO ji = 1, nbpac
591	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )
592	ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) &
593	& / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb
594	END DO
595	END DO
596	END DO
597
598	!------------
599	! Update age
600	!------------
601	DO jl = 1, jpl
602	DO ji = 1, nbpac
603	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes
604	zoa_i_ac(ji,jl) = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
605	END DO
606	END DO
607
608	!-----------------
609	! Update salinity
610	!-----------------
611	!clem IF( num_sal == 2 ) THEN
612	DO jl = 1, jpl
613	DO ji = 1, nbpac
614	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes
615	zdv = zv_i_ac(ji,jl) - zv_old(ji,jl)
616	zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif
617	END DO
618	END DO
619	!clem ENDIF
620
621	!--------------------------------
622	! Update mass/salt fluxes (clem)
623	!--------------------------------
624	DO jl = 1, jpl
625	DO ji = 1, nbpac
626	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes
627	zdv = zv_i_ac(ji,jl) - zv_old(ji,jl)
628	rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb
629	sfx_thd_1d(ji) = sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb
630	END DO
631	END DO
632
633	!------------------------------------------------------------------------------!
634	! 8) Change 2D vectors to 1D vectors
635	!------------------------------------------------------------------------------!
636	DO jl = 1, jpl
637	CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ac (1:nbpac,jl), jpi, jpj )
638	CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj )
639	CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj )
640	!clem IF ( num_sal == 2 ) &
641	CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj )
642	DO jk = 1, nlay_i
643	CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj )
644	END DO
645	END DO
646	CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj )
647	CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj )
648	CALL tab_1d_2d( nbpac, rdq_ice, npac(1:nbpac), rdq_ice_1d(1:nbpac), jpi, jpj )
649	!
650	ENDIF ! nbpac > 0
651
652	!------------------------------------------------------------------------------!
653	! 9) Change units for e_i
654	!------------------------------------------------------------------------------!
655	DO jl = 1, jpl
656	DO jk = 1, nlay_i ! heat content in 10^9 Joules
657	e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i ) / unit_fac
658	END DO
659	END DO
660
661	!------------------------------------------------------------------------------\|
662	! 10) Conservation check and changes in each ice category
663	!------------------------------------------------------------------------------\|
664	IF( con_i ) THEN
665	CALL lim_column_sum (jpl, v_i, vt_i_final)
666	fieldid = 'v_i, limthd_lac'
667	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
668	!
669	CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)
670	fieldid = 'e_i, limthd_lac'
671	CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid)
672	!
673	CALL lim_column_sum (jpl, v_s, vt_s_final)
674	fieldid = 'v_s, limthd_lac'
675	CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)
676	!
677	! CALL lim_column_sum (jpl, e_s(:,:,1,:) , et_s_init)
678	! fieldid = 'e_s, limthd_lac'
679	! CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)
680	IF( ln_nicep ) THEN
681	DO ji = mi0(jiindx), mi1(jiindx)
682	DO jj = mj0(jjindx), mj1(jjindx)
683	WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj)
684	WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj)
685	WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj)
686	WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj)
687	END DO
688	END DO
689	ENDIF
690	!
691	ENDIF
692	!
693	CALL wrk_dealloc( jpij, zcatac ) ! integer
694	CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
695	CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
696	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 )
697	CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac )
698	CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
699	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 )
700	!
701	END SUBROUTINE lim_thd_lac
702
703	#else
704	!!----------------------------------------------------------------------
705	!! Default option NO LIM3 sea-ice model
706	!!----------------------------------------------------------------------
707	CONTAINS
708	SUBROUTINE lim_thd_lac ! Empty routine
709	END SUBROUTINE lim_thd_lac
710	#endif
711
712	!!======================================================================
713	END MODULE limthd_lac

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

Download in other formats: