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 @ 4635

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

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

Download in other formats: