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

Context Navigation

source: trunk/NEMO/LIM_SRC_3/limthd_lac.F90 @ 865

Last change on this file since 865 was 865, checked in by rblod, 16 years ago
Correct a bug and clean comments in limthd_lac, see ticket #79
File size: 34.3 KB

Line
1	MODULE limthd_lac
2	#if defined key_lim3
3	!!----------------------------------------------------------------------
4	!! 'key_lim3' LIM3 sea-ice model
5	!!----------------------------------------------------------------------
6	!!======================================================================
7	!! * MODULE limthd_lac *
8	!! lateral thermodynamic growth of the ice
9	!!======================================================================
10
11	!!----------------------------------------------------------------------
12	!! lim_lat_acr : lateral accretion of ice
13	!! * Modules used
14	USE par_oce ! ocean parameters
15	USE dom_oce
16	USE in_out_manager
17	USE phycst
18	USE ice_oce ! ice variables
19	USE thd_ice
20	USE dom_ice
21	USE par_ice
22	USE ice
23	USE iceini
24	USE limtab
25	USE taumod
26	USE blk_oce
27	USE limcons
28
29	IMPLICIT NONE
30	PRIVATE
31
32	!! * Routine accessibility
33	PUBLIC lim_thd_lac ! called by lim_thd
34
35	!! * Module variables
36	REAL(wp) :: & ! constant values
37	epsi20 = 1.e-20 , &
38	epsi13 = 1.e-13 , &
39	epsi11 = 1.e-13 , &
40	epsi03 = 1.e-03 , &
41	epsi06 = 1.e-06 , &
42	zeps = 1.e-10 , &
43	zzero = 0.e0 , &
44	zone = 1.e0
45
46	!!----------------------------------------------------------------------
47	!! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008)
48	!! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limthd_lac.F90,v 1.5 2005/03/27 18:34:42 opalod Exp $
49	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
50	!!----------------------------------------------------------------------
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	!! ** References : Not available yet
81	!!
82	!! History :
83	!! 3.0 ! 12-05 (M. Vancoppenolle) Thorough rewrite of the routine
84	!! Salinity variations in sea ice,
85	!! Multi-layer code
86	!! 3.1 ! 01-06 (M. Vancoppenolle) ITD
87	!! 3.2 ! 04-07 (M. Vancoppenolle) Mass and energy conservation tested
88	!!------------------------------------------------------------------------
89	!! * Arguments
90	!! * Local variables
91	INTEGER :: &
92	ji,jj,jk,jl,jm , & !: dummy loop indices
93	layer , & !: layer index
94	nbpac !: nb of pts for lateral accretion
95
96	INTEGER :: &
97	zji , & !: ji of dummy test point
98	zjj , & !: jj of dummy test point
99	iter !: iteration for frazil ice computation
100
101	INTEGER, DIMENSION(jpij) :: &
102	zcatac , & !: indexes of categories where new ice grows
103	zswinew !: switch for new ice or not
104
105	REAL(wp), DIMENSION(jpij) :: &
106	zv_newice , & !: volume of accreted ice
107	za_newice , & !: fractional area of accreted ice
108	zh_newice , & !: thickness of accreted ice
109	ze_newice , & !: heat content of accreted ice
110	zs_newice , & !: salinity of accreted ice
111	zo_newice , & !: age of accreted ice
112	zdv_res , & !: residual volume in case of excessive heat budget
113	zda_res , & !: residual area in case of excessive heat budget
114	zat_i_ac , & !: total ice fraction
115	zat_i_lev , & !: total ice fraction for level ice only (type 1)
116	zdh_frazb , & !: accretion of frazil ice at the ice bottom
117	zvrel_ac !: relative ice / frazil velocity (1D vector)
118
119	REAL(wp), DIMENSION(jpij,jpl) :: &
120	zhice_old , & !: previous ice thickness
121	zdummy , & !: dummy thickness of new ice
122	zdhicbot , & !: thickness of new ice which is accreted vertically
123	zv_old , & !: old volume of ice in category jl
124	za_old , & !: old area of ice in category jl
125	za_i_ac , & !: 1-D version of a_i
126	zv_i_ac , & !: 1-D version of v_i
127	zoa_i_ac , & !: 1-D version of oa_i
128	zsmv_i_ac !: 1-D version of smv_i
129
130	REAL(wp), DIMENSION(jpij,jkmax,jpl) :: &
131	ze_i_ac !: 1-D version of e_i
132
133	REAL(wp), DIMENSION(jpij) :: &
134	zqbgow , & !: heat budget of the open water (negative)
135	zdhex !: excessively thick accreted sea ice (hlead-hice)
136
137	REAL(wp) :: &
138	ztmelts , & !: melting point of an ice layer
139	zdv , & !: increase in ice volume in each category
140	zfrazb !: fraction of frazil ice accreted at the ice bottom
141
142	! Redistribution of energy after bottom accretion
143	REAL(wp) :: & !: Energy redistribution
144	zqold , & !: old ice enthalpy
145	zweight , & !: weight of redistribution
146	zeps6 , & !: epsilon value
147	zalphai , & !: factor describing how old and new layers overlap each other [m]
148	zindb
149
150	REAL(wp), DIMENSION(jpij,jkmax+1,jpl) :: &
151	zqm0 , & !: old layer-system heat content
152	zthick0 !: old ice thickness
153
154	! Frazil ice collection thickness
155	LOGICAL :: & !: iterate frazil ice collection thickness
156	iterate_frazil
157
158	REAL(wp), DIMENSION(jpi,jpj) :: &
159	zvrel !: relative ice / frazil velocity
160
161	REAL(wp) :: &
162	zgamafr , & !: mult. coeff. between frazil vel. and wind speed
163	ztenagm , & !: square root of wind stress
164	zvfrx , & !: x-component of frazil velocity
165	zvfry , & !: y-component of frazil velocity
166	zvgx , & !: x-component of ice velocity
167	zvgy , & !: y-component of ice velocity
168	ztaux , & !: x-component of wind stress
169	ztauy , & !: y-component of wind stress
170	ztwogp , & !: dummy factor including reduced gravity
171	zvrel2 , & !: square of the relative ice-frazil velocity
172	zf , & !: F for Newton-Raphson procedure
173	zfp , & !: dF for Newton-Raphson procedure
174	zhicol_new , & !: updated collection thickness
175	zsqcd , & !: 1 / square root of ( airdensity * drag )
176	zhicrit !: minimum thickness of frazil ice
177
178	! Variables for energy conservation
179	REAL (wp), DIMENSION(jpi,jpj) :: & !
180	vt_i_init, vt_i_final, & ! ice volume summed over categories
181	vt_s_init, vt_s_final, & ! snow volume summed over categories
182	et_i_init, et_i_final, & ! ice energy summed over categories
183	et_s_init, et_s_final ! snow energy summed over categories
184
185	REAL(wp) :: &
186	zde ! :increment of energy in category jl
187
188	CHARACTER (len = 15) :: fieldid
189
190	!!-----------------------------------------------------------------------!
191
192	et_i_init(:,:) = 0.0
193	et_s_init(:,:) = 0.0
194	vt_i_init(:,:) = 0.0
195	vt_s_init(:,:) = 0.0
196	zeps6 = 1.0e-6
197
198	!------------------------------------------------------------------------------!
199	! 1) Conservation check and changes in each ice category
200	!------------------------------------------------------------------------------!
201	IF ( con_i ) THEN
202	CALL lim_column_sum (jpl, v_i, vt_i_init)
203	CALL lim_column_sum (jpl, v_s, vt_s_init)
204	CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init)
205	CALL lim_column_sum (jpl, e_s(:,:,1,:) , et_s_init)
206	ENDIF
207
208	!------------------------------------------------------------------------------\|
209	! 2) Convert units for ice internal energy
210	!------------------------------------------------------------------------------\|
211	DO jl = 1, jpl
212	DO jk = 1, nlay_i
213	DO jj = 1, jpj
214	DO ji = 1, jpi
215	!Energy of melting q(S,T) [J.m-3]
216	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / &
217	MAX( area(ji,jj) * v_i(ji,jj,jl) , zeps ) * &
218	nlay_i
219	zindb = 1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl))) !0 if no ice and 1 if yes
220	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl)unit_faczindb
221	END DO
222	END DO
223	END DO
224	END DO
225
226	!------------------------------------------------------------------------------!
227	! 3) Collection thickness of ice formed in leads and polynyas
228	!------------------------------------------------------------------------------!
229	! hicol is the thickness of new ice formed in open water
230	! hicol can be either prescribed (frazswi = 0)
231	! or computed (frazswi = 1)
232	! Frazil ice forms in open water, is transported by wind
233	! accumulates at the edge of the consolidated ice edge
234	! where it forms aggregates of a specific thickness called
235	! collection thickness.
236
237	! Note : the following algorithm currently breaks vectorization
238	!
239
240	zvrel(:,:) = 0.0
241
242	! Default new ice thickness
243	DO jj = 1, jpj
244	DO ji = 1, jpi
245	hicol(ji,jj) = hiccrit(1)
246	END DO
247	END DO
248
249	IF (fraz_swi.eq.1.0) THEN
250
251	!--------------------
252	! Physical constants
253	!--------------------
254	hicol(:,:) = 0.0
255
256	zhicrit = 0.04 ! frazil ice thickness
257	ztwogp = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav
258	zsqcd = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)
259	zgamafr = 0.03
260
261	DO jj = 1, jpj
262	DO ji = 1, jpi
263
264	IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
265	!-------------
266	! Wind stress
267	!-------------
268	! C-grid wind stress components
269	ztaux = ( gtaux(ji-1,jj ) * tmu(ji-1,jj ) &
270	+ gtaux(ji ,jj ) * tmu(ji ,jj ) ) / 2.0
271	ztauy = ( gtauy(ji ,jj-1) * tmv(ji ,jj-1) &
272	+ gtauy(ji ,jj ) * tmv(ji ,jj ) ) / 2.0
273	! Square root of wind stress
274	ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
275
276	!---------------------
277	! Frazil ice velocity
278	!---------------------
279	zvfrx = zgamafr * zsqcd * ztaux / MAX(ztenagm,zeps)
280	zvfry = zgamafr * zsqcd * ztauy / MAX(ztenagm,zeps)
281
282	!-------------------
283	! Pack ice velocity
284	!-------------------
285	! C-grid ice velocity
286	zindb = MAX(0.0, SIGN(1.0, at_i(ji,jj) ))
287	zvgx = zindb * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) &
288	+ u_ice(ji,jj ) * tmu(ji ,jj ) ) / 2.0
289	zvgy = zindb * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) &
290	+ v_ice(ji,jj ) * tmv(ji ,jj ) ) / 2.0
291
292	!-----------------------------------
293	! Relative frazil/pack ice velocity
294	!-----------------------------------
295	! absolute relative velocity
296	zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) + &
297	( zvfry - zvgy ) * ( zvfry - zvgy ) &
298	, 0.15 * 0.15 )
299	zvrel(ji,jj) = SQRT(zvrel2)
300
301	!---------------------
302	! Iterative procedure
303	!---------------------
304	hicol(ji,jj) = zhicrit + 0.1
305	hicol(ji,jj) = zhicrit + hicol(ji,jj) / &
306	( hicol(ji,jj) * hicol(ji,jj) - &
307	zhicrit * zhicrit ) * ztwogp * zvrel2
308
309	iter = 1
310	iterate_frazil = .true.
311
312	DO WHILE ( iter .LT. 100 .AND. iterate_frazil )
313	zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)2 - zhicrit2 ) &
314	- hicol(ji,jj) * zhicrit * ztwogp * zvrel2
315	zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0*hicol(ji,jj) + zhicrit ) &
316	- zhicrit * ztwogp * zvrel2
317	zhicol_new = hicol(ji,jj) - zf/zfp
318	hicol(ji,jj) = zhicol_new
319
320	iter = iter + 1
321
322	END DO ! do while
323
324	ENDIF ! end of selection of pixels where ice forms
325
326	END DO ! loop on ji ends
327	END DO ! loop on jj ends
328
329	ENDIF ! End of computation of frazil ice collection thickness
330
331	!------------------------------------------------------------------------------!
332	! 4) Identify grid points where new ice forms
333	!------------------------------------------------------------------------------!
334
335	!-------------------------------------
336	! Select points for new ice formation
337	!-------------------------------------
338	! This occurs if open water energy budget is negative
339	nbpac = 0
340	DO jj = 1, jpj
341	DO ji = 1, jpi
342	IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
343	nbpac = nbpac + 1
344	npac( nbpac ) = (jj - 1) * jpi + ji
345	IF ( (ji.eq.jiindex).AND.(jj.eq.jjindex) ) THEN
346	jiindex_1d = nbpac
347	ENDIF
348	ENDIF
349	END DO
350	END DO
351
352	IF(lwp) THEN
353	WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
354	ENDIF
355
356	!------------------------------
357	! Move from 2-D to 1-D vectors
358	!------------------------------
359	! If ocean gains heat do nothing
360	! 0therwise compute new ice formation
361
362	IF ( nbpac > 0 ) THEN
363
364	CALL tab_2d_1d( nbpac, zat_i_ac (1:nbpac) , at_i , &
365	jpi, jpj, npac(1:nbpac) )
366	DO jl = 1, jpl
367	CALL tab_2d_1d( nbpac, za_i_ac(1:nbpac,jl) , a_i(:,:,jl) , &
368	jpi, jpj, npac(1:nbpac) )
369	CALL tab_2d_1d( nbpac, zv_i_ac(1:nbpac,jl) , v_i(:,:,jl) , &
370	jpi, jpj, npac(1:nbpac) )
371	CALL tab_2d_1d( nbpac, zoa_i_ac(1:nbpac,jl) , oa_i(:,:,jl) , &
372	jpi, jpj, npac(1:nbpac) )
373	CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), &
374	jpi, jpj, npac(1:nbpac) )
375	DO jk = 1, nlay_i
376	CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , &
377	jpi, jpj, npac(1:nbpac) )
378	END DO ! jk
379	END DO ! jl
380
381	CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , &
382	jpi, jpj, npac(1:nbpac) )
383	CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , &
384	jpi, jpj, npac(1:nbpac) )
385	CALL tab_2d_1d( nbpac, t_bo_b (1:nbpac) , t_bo , &
386	jpi, jpj, npac(1:nbpac) )
387	CALL tab_2d_1d( nbpac, fseqv_1d (1:nbpac) , fseqv , &
388	jpi, jpj, npac(1:nbpac) )
389	CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , &
390	jpi, jpj, npac(1:nbpac) )
391	CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , &
392	jpi, jpj, npac(1:nbpac) )
393
394	!------------------------------------------------------------------------------!
395	! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
396	!------------------------------------------------------------------------------!
397
398	!----------------------
399	! Thickness of new ice
400	!----------------------
401	DO ji = 1, nbpac
402	zh_newice(ji) = hiccrit(1)
403	END DO
404	IF ( fraz_swi .EQ. 1.0 ) zh_newice(:) = hicol_b(:)
405
406	!----------------------
407	! Salinity of new ice
408	!----------------------
409
410	IF ( num_sal .EQ. 1 ) THEN
411	zs_newice(:) = bulk_sal
412	ENDIF ! num_sal
413
414	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
415
416	DO ji = 1, nbpac
417	zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max )
418	zji = MOD( npac(ji) - 1, jpi ) + 1
419	zjj = ( npac(ji) - 1 ) / jpi + 1
420	zs_newice(ji) = MIN( 0.5*sss_io(zji,zjj) , zs_newice(ji) )
421	END DO ! jl
422
423	ENDIF ! num_sal
424
425	IF ( num_sal .EQ. 3 ) THEN
426	zs_newice(:) = 2.3
427	ENDIF ! num_sal
428
429	!-------------------------
430	! Heat content of new ice
431	!-------------------------
432	! We assume that new ice is formed at the seawater freezing point
433	DO ji = 1, nbpac
434	ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K)
435	ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) &
436	+ lfus * ( 1.0 - ( ztmelts - rtt ) &
437	/ ( t_bo_b(ji) - rtt ) ) &
438	- rcp * ( ztmelts-rtt ) )
439	ze_newice(ji) = MAX( ze_newice(ji) , 0.0 ) + &
440	MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) ) &
441	* rhoic * lfus
442	END DO ! ji
443	!----------------
444	! Age of new ice
445	!----------------
446	DO ji = 1, nbpac
447	zo_newice(ji) = 0.0
448	END DO ! ji
449
450	!--------------------------
451	! Open water energy budget
452	!--------------------------
453	DO ji = 1, nbpac
454	zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) !<0
455	END DO ! ji
456
457	!-------------------
458	! Volume of new ice
459	!-------------------
460	DO ji = 1, nbpac
461	zv_newice(ji) = - zqbgow(ji) / ze_newice(ji)
462
463	! A fraction zfrazb of frazil ice is accreted at the ice bottom
464	zfrazb = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) &
465	+ 1.0 ) / 2.0 * maxfrazb
466	zdh_frazb(ji) = zfrazb*zv_newice(ji)
467	zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
468	END DO
469
470	!---------------------------------
471	! Salt flux due to new ice growth
472	!---------------------------------
473	IF ( ( num_sal .EQ. 4 ) ) THEN
474	DO ji = 1, nbpac
475	zji = MOD( npac(ji) - 1, jpi ) + 1
476	zjj = ( npac(ji) - 1 ) / jpi + 1
477	fseqv_1d(ji) = fseqv_1d(ji) + &
478	( sss_io(zji,zjj) - bulk_sal ) * rhoic * &
479	zv_newice(ji) / rdt_ice
480	END DO
481	ELSE
482	DO ji = 1, nbpac
483	zji = MOD( npac(ji) - 1, jpi ) + 1
484	zjj = ( npac(ji) - 1 ) / jpi + 1
485	fseqv_1d(ji) = fseqv_1d(ji) + &
486	( sss_io(zji,zjj) - zs_newice(ji) ) * rhoic * &
487	zv_newice(ji) / rdt_ice
488	END DO ! ji
489	ENDIF
490
491	!------------------------------------
492	! Diags for energy conservation test
493	!------------------------------------
494	DO ji = 1, nbpac
495	! Volume
496	zji = MOD( npac(ji) - 1, jpi ) + 1
497	zjj = ( npac(ji) - 1 ) / jpi + 1
498	vt_i_init(zji,zjj) = vt_i_init(zji,zjj) + zv_newice(ji)
499	! Energy
500	zde = ze_newice(ji) / unit_fac
501	zde = zde * area(zji,zjj) * zv_newice(ji)
502	et_i_init(zji,zjj) = et_i_init(zji,zjj) + zde
503	END DO
504
505	! keep new ice volume in memory
506	CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , &
507	jpi, jpj )
508
509	!-----------------
510	! Area of new ice
511	!-----------------
512	DO ji = 1, nbpac
513	za_newice(ji) = zv_newice(ji) / zh_newice(ji)
514	! diagnostic
515	zji = MOD( npac(ji) - 1, jpi ) + 1
516	zjj = ( npac(ji) - 1 ) / jpi + 1
517	diag_lat_gr(zji,zjj) = zv_newice(ji) / rdt_ice
518	END DO !ji
519
520	!------------------------------------------------------------------------------!
521	! 6) Redistribute new ice area and volume into ice categories !
522	!------------------------------------------------------------------------------!
523
524	!-----------------------------------------
525	! Keep old ice areas and volume in memory
526	!-----------------------------------------
527	zv_old(:,:) = zv_i_ac(:,:)
528	za_old(:,:) = za_i_ac(:,:)
529
530	!-------------------------------------------
531	! Compute excessive new ice area and volume
532	!-------------------------------------------
533	! If lateral ice growth gives an ice concentration gt 1, then
534	! we keep the excessive volume in memory and attribute it later
535	! to bottom accretion
536	DO ji = 1, nbpac
537	! vectorize
538	IF ( za_newice(ji) .GT. ( 1.0 - zat_i_ac(ji) ) ) THEN
539	zda_res(ji) = za_newice(ji) - (1.0 - zat_i_ac(ji) )
540	zdv_res(ji) = zda_res(ji) * zh_newice(ji)
541	za_newice(ji) = za_newice(ji) - zda_res(ji)
542	zv_newice(ji) = zv_newice(ji) - zdv_res(ji)
543	ELSE
544	zda_res(ji) = 0.0
545	zdv_res(ji) = 0.0
546	ENDIF
547	END DO ! ji
548
549	!------------------------------------------------
550	! Laterally redistribute new ice volume and area
551	!------------------------------------------------
552	zat_i_ac(:) = 0.0
553
554	DO jl = 1, jpl
555	DO ji = 1, nbpac
556	! vectorize
557	IF ( ( hi_max(jl-1) .LT. zh_newice(ji) ) &
558	.AND. ( zh_newice(ji) .LE. hi_max(jl) ) ) THEN
559	za_i_ac(ji,jl) = za_i_ac(ji,jl) + za_newice(ji)
560	zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zv_newice(ji)
561	zat_i_ac(ji) = zat_i_ac(ji) + za_i_ac(ji,jl)
562	zcatac(ji) = jl
563	ENDIF
564	END DO ! ji
565	END DO ! jl
566
567	!----------------------------------
568	! Heat content - lateral accretion
569	!----------------------------------
570	DO ji = 1, nbpac
571	jl = zcatac(ji) ! categroy in which new ice is put
572	! zindb = 0 if no ice and 1 if yes
573	zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , -za_old(ji,jl) ) )
574	! old ice thickness
575	zhice_old(ji,jl) = zv_old(ji,jl) &
576	/ MAX ( za_old(ji,jl) , zeps ) * zindb
577	! difference in thickness
578	zdhex(ji) = MAX( 0.0, zh_newice(ji) - zhice_old(ji,jl) )
579	! is ice totally new in category jl ?
580	zswinew(ji) = MAX( 0.0, SIGN( 1.0 , - za_old(ji,jl) + epsi11 ) )
581	END DO
582
583	DO jk = 1, nlay_i
584	DO ji = 1, nbpac
585	jl = zcatac(ji)
586	zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ]
587	zalphai = MIN( zhice_old(ji,jl) * jk / nlay_i , &
588	zh_newice(ji) ) &
589	- MIN( zhice_old(ji,jl) * ( jk - 1 ) &
590	/ nlay_i , zh_newice(ji) )
591	ze_i_ac(ji,jk,jl) = &
592	zswinew(ji) * ze_newice(ji) &
593	+ ( 1.0 - zswinew(ji) ) * &
594	( za_old(ji,jl) * zqold * zhice_old(ji,jl) / nlay_i &
595	+ za_newice(ji) * ze_newice(ji) * zalphai &
596	+ za_newice(ji) * ze_newice(ji) * zdhex(ji) / nlay_i ) / &
597	( ( zv_i_ac(ji,jl) ) / nlay_i )
598
599	END DO !ji
600	END DO !jl
601
602	!-----------------------------------------------
603	! Add excessive volume of new ice at the bottom
604	!-----------------------------------------------
605	! If the ice concentration exceeds 1, the remaining volume of new ice
606	! is equally redistributed among all ice categories in which there is
607	! ice
608
609	! Fraction of level ice
610	jm = 1
611	zat_i_lev(:) = 0.0
612
613	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
614	DO ji = 1, nbpac
615	zat_i_lev(ji) = zat_i_lev(ji) + za_i_ac(ji,jl)
616	END DO
617	END DO
618
619	WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex, 1:jpl)
620	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
621	DO ji = 1, nbpac
622	zindb = MAX( 0.0, SIGN( 1.0, zdv_res(ji) ) )
623	zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + &
624	zindb * zdv_res(ji) * za_i_ac(ji,jl) / &
625	MAX( zat_i_lev(ji) , epsi06 )
626	END DO ! ji
627	END DO ! jl
628	WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex, 1:jpl)
629
630	!---------------------------------
631	! Heat content - bottom accretion
632	!---------------------------------
633	jm = 1
634	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
635	DO ji = 1, nbpac
636	! zindb = 0 if no ice and 1 if yes
637	zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 &
638	, - za_i_ac(ji,jl ) ) )
639	zhice_old(ji,jl) = zv_i_ac(ji,jl) / &
640	MAX( za_i_ac(ji,jl) , zeps ) * zindb
641	zdhicbot(ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , zeps ) &
642	* zindb &
643	+ zindb * zdh_frazb(ji) ! frazil ice
644	! may coalesce
645	! thickness of residual ice
646	zdummy(ji,jl) = zv_i_ac(ji,jl)/MAX(za_i_ac(ji,jl),zeps)*zindb
647	END DO !ji
648	END DO !jl
649
650	! old layers thicknesses and enthalpies
651	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
652	DO jk = 1, nlay_i
653	DO ji = 1, nbpac
654	zthick0(ji,jk,jl)= zhice_old(ji,jl) / nlay_i
655	zqm0 (ji,jk,jl)= ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)
656	END DO !ji
657	END DO !jk
658	END DO !jl
659
660	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
661	DO ji = 1, nbpac
662	zthick0(ji,nlay_i+1,jl) = zdhicbot(ji,jl)
663	zqm0 (ji,nlay_i+1,jl) = ze_newice(ji)*zdhicbot(ji,jl)
664	END DO ! ji
665	END DO ! jl
666
667	! Redistributing energy on the new grid
668	ze_i_ac(:,:,:) = 0.0
669	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
670	DO jk = 1, nlay_i
671	DO layer = 1, nlay_i + 1
672	DO ji = 1, nbpac
673	zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , &
674	- za_i_ac(ji,jl ) ) )
675	! Redistributing energy on the new grid
676	zweight = MAX ( &
677	MIN( zhice_old(ji,jl) * layer , zdummy(ji,jl) * jk ) - &
678	MAX( zhice_old(ji,jl) * ( layer - 1 ) , zdummy(ji,jl) * &
679	( jk - 1 ) ) , 0.0 ) &
680	/ ( MAX(nlay_i * zthick0(ji,layer,jl),zeps) ) * zindb
681	ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + &
682	zweight * zqm0(ji,layer,jl)
683	END DO ! ji
684	END DO ! layer
685	END DO ! jk
686	END DO ! jl
687
688	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
689	DO jk = 1, nlay_i
690	DO ji = 1, nbpac
691	zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 &
692	, - zv_i_ac(ji,jl) ) ) !0 if no ice
693	ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) / &
694	MAX( zv_i_ac(ji,jl) , zeps) &
695	* za_i_ac(ji,jl) * nlay_i * zindb
696	END DO
697	END DO
698	END DO
699
700	!------------
701	! Update age
702	!------------
703	DO jl = 1, jpl
704	DO ji = 1, nbpac
705	!--ice age
706	zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - &
707	za_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes
708	zoa_i_ac(ji,jl) = za_old(ji,jl) * zoa_i_ac(ji,jl) / &
709	MAX( za_i_ac(ji,jl) , zeps ) * zindb
710	END DO ! ji
711	END DO ! jl
712
713	!-----------------
714	! Update salinity
715	!-----------------
716	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
717
718	DO jl = 1, jpl
719	DO ji = 1, nbpac
720	!zindb = 0 if no ice and 1 if yes
721	zindb = 1.0 - MAX( 0.0 , SIGN( 1.0 , - &
722	zv_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes
723	zdv = zv_i_ac(ji,jl) - zv_old(ji,jl)
724	zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * &
725	zindb
726	END DO ! ji
727	END DO ! jl
728
729	ENDIF ! num_sal
730
731
732	!------------------------------------------------------------------------------!
733	! 8) Change 2D vectors to 1D vectors
734	!------------------------------------------------------------------------------!
735
736	DO jl = 1, jpl
737	CALL tab_1d_2d( nbpac, a_i(:,:,jl) , npac(1:nbpac) , &
738	za_i_ac(1:nbpac,jl) , jpi, jpj )
739	CALL tab_1d_2d( nbpac, v_i(:,:,jl) , npac(1:nbpac) , &
740	zv_i_ac(1:nbpac,jl) , jpi, jpj )
741	CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac) , &
742	zoa_i_ac(1:nbpac,jl), jpi, jpj )
743	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
744	CALL tab_1d_2d( nbpac, smv_i(:,:,jl) , npac(1:nbpac) , &
745	zsmv_i_ac(1:nbpac,jl) , jpi, jpj )
746	DO jk = 1, nlay_i
747	CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl) , npac(1:nbpac), &
748	ze_i_ac(1:nbpac,jk,jl), jpi, jpj )
749	END DO ! jk
750	END DO !jl
751	CALL tab_1d_2d( nbpac, fseqv , npac(1:nbpac), fseqv_1d (1:nbpac) , &
752	jpi, jpj )
753
754	ENDIF ! nbpac > 0
755
756	!------------------------------------------------------------------------------!
757	! 9) Change units for e_i
758	!------------------------------------------------------------------------------!
759
760	DO jl = 1, jpl
761	DO jk = 1, nlay_i
762	DO jj = 1, jpj
763	DO ji = 1, jpi
764	! Correct dimensions to avoid big values
765	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac
766
767	! Mutliply by ice volume, and divide by number
768	! of layers to get heat content in 10^9 Joules
769	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &
770	area(ji,jj) * v_i(ji,jj,jl) / &
771	nlay_i
772	END DO
773	END DO
774	END DO
775	END DO
776
777	!------------------------------------------------------------------------------\|
778	! 10) Conservation check and changes in each ice category
779	!------------------------------------------------------------------------------\|
780
781	IF ( con_i ) THEN
782	CALL lim_column_sum (jpl, v_i, vt_i_final)
783	fieldid = 'v_i, limthd_lac'
784	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
785
786	CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)
787	fieldid = 'e_i, limthd_lac'
788	CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid)
789
790	CALL lim_column_sum (jpl, v_s, vt_s_final)
791	fieldid = 'v_s, limthd_lac'
792	CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)
793
794	! CALL lim_column_sum (jpl, e_s(:,:,1,:) , et_s_init)
795	! fieldid = 'e_s, limthd_lac'
796	! CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)
797
798	WRITE(numout,*) ' vt_i_init : ', vt_i_init(jiindex,jjindex)
799	WRITE(numout,*) ' vt_i_final: ', vt_i_final(jiindex,jjindex)
800	WRITE(numout,*) ' et_i_init : ', et_i_init(jiindex,jjindex)
801	WRITE(numout,*) ' et_i_final: ', et_i_final(jiindex,jjindex)
802
803	ENDIF
804
805	END SUBROUTINE lim_thd_lac
806
807	#else
808	!!======================================================================
809	!! * MODULE limthd_lac *
810	!! no sea ice model
811	!!======================================================================
812	CONTAINS
813	SUBROUTINE lim_thd_lac ! Empty routine
814	END SUBROUTINE lim_thd_lac
815	#endif
816	END MODULE limthd_lac

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

Download in other formats: