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

Context Navigation

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 2833

Last change on this file since 2833 was 2777, checked in by smasson, 13 years ago
LIM3 compiling and (partly?) running in v3_3_1, see ticket#817
Property svn:keywords set to `Id`
File size: 32.0 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 wrk_nemo ! workspace manager
29	USE in_out_manager ! I/O manager
30	USE lib_mpp ! MPP library
31
32	IMPLICIT NONE
33	PRIVATE
34
35	PUBLIC lim_thd_lac ! called by lim_thd
36
37	REAL(wp) :: epsi20 = 1e-20_wp ! constant values
38	REAL(wp) :: epsi13 = 1e-13_wp !
39	REAL(wp) :: epsi11 = 1e-11_wp !
40	REAL(wp) :: epsi10 = 1e-10_wp !
41	REAL(wp) :: epsi06 = 1e-06_wp !
42	REAL(wp) :: epsi03 = 1e-03_wp !
43	REAL(wp) :: zzero = 0._wp !
44	REAL(wp) :: zone = 1._wp !
45
46	!!----------------------------------------------------------------------
47	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
48	!! $Id$
49	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
50	!!----------------------------------------------------------------------
51	CONTAINS
52
53	SUBROUTINE lim_thd_lac
54	!!-------------------------------------------------------------------
55	!! * ROUTINE lim_thd_lac *
56	!!
57	!! ** Purpose : Computation of the evolution of the ice thickness and
58	!! concentration as a function of the heat balance in the leads.
59	!! It is only used for lateral accretion
60	!!
61	!! ** Method : Ice is formed in the open water when ocean lose heat
62	!! (heat budget of open water Bl is negative) .
63	!! Computation of the increase of 1-A (ice concentration) fol-
64	!! lowing the law :
65	!! (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
66	!! where - h0 is the thickness of ice created in the lead
67	!! - a is a minimum fraction for leads
68	!! - F is a monotonic non-increasing function defined as:
69	!! F(X)=( 1 - Xexld )(1.0/exld)
70	!! - exld is the exponent closure rate (=2 default val.)
71	!!
72	!! ** Action : - Adjustment of snow and ice thicknesses and heat
73	!! content in brine pockets
74	!! - Updating ice internal temperature
75	!! - Computation of variation of ice volume and mass
76	!! - Computation of frldb after lateral accretion and
77	!! update ht_s_b, ht_i_b and tbif_1d(:,:)
78	!!------------------------------------------------------------------------
79	USE wrk_nemo, ONLY : vt_i_init => wrk_2d_1 , vt_i_final => wrk_2d_4 , et_i_init => wrk_2d_7
80	USE wrk_nemo, ONLY : vt_s_init => wrk_2d_2 , vt_s_final => wrk_2d_5 , et_s_init => wrk_2d_8
81	USE wrk_nemo, ONLY : zvrel => wrk_2d_3 , et_i_final => wrk_2d_6
82	!
83	INTEGER :: ji,jj,jk,jl,jm ! dummy loop indices
84	INTEGER :: layer, nbpac ! local integers
85	INTEGER :: zji, zjj, iter ! - -
86	REAL(wp) :: ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde ! local scalars
87	REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - -
88	REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - -
89	LOGICAL :: iterate_frazil ! iterate frazil ice collection thickness
90	CHARACTER (len = 15) :: fieldid
91	!
92	INTEGER, DIMENSION(jpij) :: zcatac ! indexes of categories where new ice grows
93
94	REAL(wp), DIMENSION(jpij,jpl) :: zhice_old ! previous ice thickness
95	REAL(wp), DIMENSION(jpij,jpl) :: zdummy ! dummy thickness of new ice
96	REAL(wp), DIMENSION(jpij,jpl) :: zdhicbot ! thickness of new ice which is accreted vertically
97	REAL(wp), DIMENSION(jpij,jpl) :: zv_old ! old volume of ice in category jl
98	REAL(wp), DIMENSION(jpij,jpl) :: za_old ! old area of ice in category jl
99	REAL(wp), DIMENSION(jpij,jpl) :: za_i_ac ! 1-D version of a_i
100	REAL(wp), DIMENSION(jpij,jpl) :: zv_i_ac ! 1-D version of v_i
101	REAL(wp), DIMENSION(jpij,jpl) :: zoa_i_ac ! 1-D version of oa_i
102	REAL(wp), DIMENSION(jpij,jpl) :: zsmv_i_ac ! 1-D version of smv_i
103
104	REAL(wp), DIMENSION(jpij,jkmax ,jpl) :: ze_i_ac !: 1-D version of e_i
105	REAL(wp), DIMENSION(jpij,jkmax+1,jpl) :: zqm0 ! old layer-system heat content
106	REAL(wp), DIMENSION(jpij,jkmax+1,jpl) :: zthick0 ! old ice thickness
107
108	REAL(wp), POINTER, DIMENSION(:) :: zv_newice, zh_newice, zs_newice, zdv_res, zat_i_ac , zdh_frazb, zqbgow
109	REAL(wp), POINTER, DIMENSION(:) :: za_newice, ze_newice, zo_newice, zda_res, zat_i_lev, zvrel_ac , zdhex
110	REAL(wp), POINTER, DIMENSION(:) :: zswinew
111	!!-----------------------------------------------------------------------!
112
113	IF( wrk_in_use(1, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15) .OR. &
114	& wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN
115	CALL ctl_stop('lim_thd_lac : requestead workspace arrays unavailable.') ; RETURN
116	END IF
117	! Set-up pointers to sub-arrays of workspace arrays
118	zv_newice => wrk_1d_1 (1:jpij) ! volume of accreted ice
119	za_newice => wrk_1d_2 (1:jpij) ! fractional area of accreted ice
120	zh_newice => wrk_1d_3 (1:jpij) ! thickness of accreted ice
121	ze_newice => wrk_1d_4 (1:jpij) ! heat content of accreted ice
122	zs_newice => wrk_1d_5 (1:jpij) ! salinity of accreted ice
123	zo_newice => wrk_1d_6 (1:jpij) ! age of accreted ice
124	zdv_res => wrk_1d_7 (1:jpij) ! residual volume in case of excessive heat budget
125	zda_res => wrk_1d_8 (1:jpij) ! residual area in case of excessive heat budget
126	zat_i_ac => wrk_1d_9 (1:jpij) ! total ice fraction
127	zat_i_lev => wrk_1d_10(1:jpij) ! total ice fraction for level ice only (type 1)
128	zdh_frazb => wrk_1d_11(1:jpij) ! accretion of frazil ice at the ice bottom
129	zvrel_ac => wrk_1d_12(1:jpij) ! relative ice / frazil velocity (1D vector)
130	zqbgow => wrk_1d_13(1:jpij) ! heat budget of the open water (negative)
131	zdhex => wrk_1d_14(1:jpij) ! excessively thick accreted sea ice (hlead-hice)
132	zswinew => wrk_1d_15(1:jpij) ! switch for new ice or not
133
134
135
136	et_i_init(:,:) = 0._wp
137	et_s_init(:,:) = 0._wp
138	vt_i_init(:,:) = 0._wp
139	vt_s_init(:,:) = 0._wp
140
141	!------------------------------------------------------------------------------!
142	! 1) Conservation check and changes in each ice category
143	!------------------------------------------------------------------------------!
144	IF ( con_i ) THEN
145	CALL lim_column_sum (jpl, v_i, vt_i_init)
146	CALL lim_column_sum (jpl, v_s, vt_s_init)
147	CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init)
148	CALL lim_column_sum (jpl, e_s(:,:,1,:) , et_s_init)
149	ENDIF
150
151	!------------------------------------------------------------------------------\|
152	! 2) Convert units for ice internal energy
153	!------------------------------------------------------------------------------\|
154	DO jl = 1, jpl
155	DO jk = 1, nlay_i
156	DO jj = 1, jpj
157	DO ji = 1, jpi
158	!Energy of melting q(S,T) [J.m-3]
159	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / &
160	MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * &
161	nlay_i
162	zindb = 1.0-MAX(0.0,SIGN(1.0,-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_faczindb
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.0
184
185	! Default new ice thickness
186	DO jj = 1, jpj
187	DO ji = 1, jpi
188	hicol(ji,jj) = hiccrit(1)
189	END DO
190	END DO
191
192	IF (fraz_swi.eq.1.0) THEN
193
194	!--------------------
195	! Physical constants
196	!--------------------
197	hicol(:,:) = 0.0
198
199	zhicrit = 0.04 ! frazil ice thickness
200	ztwogp = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav
201	zsqcd = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)
202	zgamafr = 0.03
203
204	DO jj = 1, jpj
205	DO ji = 1, jpi
206
207	IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
208	!-------------
209	! Wind stress
210	!-------------
211	! C-grid wind stress components
212	ztaux = ( utau_ice(ji-1,jj ) * tmu(ji-1,jj ) &
213	& + utau_ice(ji ,jj ) * tmu(ji ,jj ) ) / 2.0
214	ztauy = ( vtau_ice(ji ,jj-1) * tmv(ji ,jj-1) &
215	& + vtau_ice(ji ,jj ) * tmv(ji ,jj ) ) / 2.0
216	! Square root of wind stress
217	ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
218
219	!---------------------
220	! Frazil ice velocity
221	!---------------------
222	zvfrx = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10)
223	zvfry = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10)
224
225	!-------------------
226	! Pack ice velocity
227	!-------------------
228	! C-grid ice velocity
229	zindb = MAX(0.0, SIGN(1.0, at_i(ji,jj) ))
230	zvgx = zindb * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) &
231	+ u_ice(ji,jj ) * tmu(ji ,jj ) ) / 2.0
232	zvgy = zindb * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) &
233	+ v_ice(ji,jj ) * tmv(ji ,jj ) ) / 2.0
234
235	!-----------------------------------
236	! Relative frazil/pack ice velocity
237	!-----------------------------------
238	! absolute relative velocity
239	zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) + &
240	( zvfry - zvgy ) * ( zvfry - zvgy ) &
241	, 0.15 * 0.15 )
242	zvrel(ji,jj) = SQRT(zvrel2)
243
244	!---------------------
245	! Iterative procedure
246	!---------------------
247	hicol(ji,jj) = zhicrit + 0.1
248	hicol(ji,jj) = zhicrit + hicol(ji,jj) / &
249	( hicol(ji,jj) * hicol(ji,jj) - &
250	zhicrit * zhicrit ) * ztwogp * zvrel2
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.e0 ) THEN
286	nbpac = nbpac + 1
287	npac( nbpac ) = (jj - 1) * jpi + ji
288	IF ( (ji.eq.jiindx).AND.(jj.eq.jjindx) ) THEN
289	jiindex_1d = nbpac
290	ENDIF
291	ENDIF
292	END DO
293	END DO
294
295	IF( ln_nicep ) THEN
296	WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
297	ENDIF
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 , &
308	jpi, jpj, npac(1:nbpac) )
309	DO jl = 1, jpl
310	CALL tab_2d_1d( nbpac, za_i_ac(1:nbpac,jl) , a_i(:,:,jl) , &
311	jpi, jpj, npac(1:nbpac) )
312	CALL tab_2d_1d( nbpac, zv_i_ac(1:nbpac,jl) , v_i(:,:,jl) , &
313	jpi, jpj, npac(1:nbpac) )
314	CALL tab_2d_1d( nbpac, zoa_i_ac(1:nbpac,jl) , oa_i(:,:,jl) , &
315	jpi, jpj, npac(1:nbpac) )
316	CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), &
317	jpi, jpj, npac(1:nbpac) )
318	DO jk = 1, nlay_i
319	CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , &
320	jpi, jpj, npac(1:nbpac) )
321	END DO ! jk
322	END DO ! jl
323
324	CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , &
325	jpi, jpj, npac(1:nbpac) )
326	CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , &
327	jpi, jpj, npac(1:nbpac) )
328	CALL tab_2d_1d( nbpac, t_bo_b (1:nbpac) , t_bo , &
329	jpi, jpj, npac(1:nbpac) )
330	CALL tab_2d_1d( nbpac, fseqv_1d (1:nbpac) , fseqv , &
331	jpi, jpj, npac(1:nbpac) )
332	CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , &
333	jpi, jpj, npac(1:nbpac) )
334	CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , &
335	jpi, jpj, npac(1:nbpac) )
336
337	!------------------------------------------------------------------------------!
338	! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
339	!------------------------------------------------------------------------------!
340
341	!----------------------
342	! Thickness of new ice
343	!----------------------
344	DO ji = 1, nbpac
345	zh_newice(ji) = hiccrit(1)
346	END DO
347	IF ( fraz_swi .EQ. 1.0 ) zh_newice(:) = hicol_b(:)
348
349	!----------------------
350	! Salinity of new ice
351	!----------------------
352
353	IF ( num_sal .EQ. 1 ) THEN
354	zs_newice(:) = bulk_sal
355	ENDIF ! num_sal
356
357	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
358
359	DO ji = 1, nbpac
360	zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max )
361	zji = MOD( npac(ji) - 1, jpi ) + 1
362	zjj = ( npac(ji) - 1 ) / jpi + 1
363	zs_newice(ji) = MIN( 0.5*sss_m(zji,zjj) , zs_newice(ji) )
364	END DO ! jl
365
366	ENDIF ! num_sal
367
368	IF ( num_sal .EQ. 3 ) THEN
369	zs_newice(:) = 2.3
370	ENDIF ! num_sal
371
372	!-------------------------
373	! Heat content of new ice
374	!-------------------------
375	! We assume that new ice is formed at the seawater freezing point
376	DO ji = 1, nbpac
377	ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K)
378	ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) &
379	+ lfus * ( 1.0 - ( ztmelts - rtt ) &
380	/ ( t_bo_b(ji) - rtt ) ) &
381	- rcp * ( ztmelts-rtt ) )
382	ze_newice(ji) = MAX( ze_newice(ji) , 0.0 ) + &
383	MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) ) &
384	* rhoic * lfus
385	END DO ! ji
386	!----------------
387	! Age of new ice
388	!----------------
389	DO ji = 1, nbpac
390	zo_newice(ji) = 0.0
391	END DO ! ji
392
393	!--------------------------
394	! Open water energy budget
395	!--------------------------
396	DO ji = 1, nbpac
397	zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) !<0
398	END DO ! ji
399
400	!-------------------
401	! Volume of new ice
402	!-------------------
403	DO ji = 1, nbpac
404	zv_newice(ji) = - zqbgow(ji) / ze_newice(ji)
405
406	! A fraction zfrazb of frazil ice is accreted at the ice bottom
407	zfrazb = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) &
408	+ 1.0 ) / 2.0 * maxfrazb
409	zdh_frazb(ji) = zfrazb*zv_newice(ji)
410	zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
411	END DO
412
413	!---------------------------------
414	! Salt flux due to new ice growth
415	!---------------------------------
416	IF ( ( num_sal .EQ. 4 ) ) THEN
417	DO ji = 1, nbpac
418	zji = MOD( npac(ji) - 1, jpi ) + 1
419	zjj = ( npac(ji) - 1 ) / jpi + 1
420	fseqv_1d(ji) = fseqv_1d(ji) + &
421	( sss_m(zji,zjj) - bulk_sal ) * rhoic * &
422	zv_newice(ji) / rdt_ice
423	END DO
424	ELSE
425	DO ji = 1, nbpac
426	zji = MOD( npac(ji) - 1, jpi ) + 1
427	zjj = ( npac(ji) - 1 ) / jpi + 1
428	fseqv_1d(ji) = fseqv_1d(ji) + &
429	( sss_m(zji,zjj) - zs_newice(ji) ) * rhoic * &
430	zv_newice(ji) / rdt_ice
431	END DO ! ji
432	ENDIF
433
434	!------------------------------------
435	! Diags for energy conservation test
436	!------------------------------------
437	DO ji = 1, nbpac
438	! Volume
439	zji = MOD( npac(ji) - 1, jpi ) + 1
440	zjj = ( npac(ji) - 1 ) / jpi + 1
441	vt_i_init(zji,zjj) = vt_i_init(zji,zjj) + zv_newice(ji)
442	! Energy
443	zde = ze_newice(ji) / unit_fac
444	zde = zde * area(zji,zjj) * zv_newice(ji)
445	et_i_init(zji,zjj) = et_i_init(zji,zjj) + zde
446	END DO
447
448	! keep new ice volume in memory
449	CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , &
450	jpi, jpj )
451
452	!-----------------
453	! Area of new ice
454	!-----------------
455	DO ji = 1, nbpac
456	za_newice(ji) = zv_newice(ji) / zh_newice(ji)
457	! diagnostic
458	zji = MOD( npac(ji) - 1, jpi ) + 1
459	zjj = ( npac(ji) - 1 ) / jpi + 1
460	diag_lat_gr(zji,zjj) = zv_newice(ji) / rdt_ice
461	END DO !ji
462
463	!------------------------------------------------------------------------------!
464	! 6) Redistribute new ice area and volume into ice categories !
465	!------------------------------------------------------------------------------!
466
467	!-----------------------------------------
468	! Keep old ice areas and volume in memory
469	!-----------------------------------------
470	zv_old(:,:) = zv_i_ac(:,:)
471	za_old(:,:) = za_i_ac(:,:)
472
473	!-------------------------------------------
474	! Compute excessive new ice area and volume
475	!-------------------------------------------
476	! If lateral ice growth gives an ice concentration gt 1, then
477	! we keep the excessive volume in memory and attribute it later
478	! to bottom accretion
479	DO ji = 1, nbpac
480	! vectorize
481	IF ( za_newice(ji) .GT. ( 1.0 - zat_i_ac(ji) ) ) THEN
482	zda_res(ji) = za_newice(ji) - (1.0 - zat_i_ac(ji) )
483	zdv_res(ji) = zda_res(ji) * zh_newice(ji)
484	za_newice(ji) = za_newice(ji) - zda_res(ji)
485	zv_newice(ji) = zv_newice(ji) - zdv_res(ji)
486	ELSE
487	zda_res(ji) = 0.0
488	zdv_res(ji) = 0.0
489	ENDIF
490	END DO ! ji
491
492	!------------------------------------------------
493	! Laterally redistribute new ice volume and area
494	!------------------------------------------------
495	zat_i_ac(:) = 0._wp
496	DO jl = 1, jpl
497	DO ji = 1, nbpac
498	IF( hi_max (jl-1) < zh_newice(ji) .AND. &
499	& zh_newice(ji) <= hi_max (jl) ) THEN
500	za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji)
501	zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji)
502	zat_i_ac(ji) = zat_i_ac(ji) + za_i_ac (ji,jl)
503	zcatac (ji) = jl
504	ENDIF
505	END DO ! ji
506	END DO ! jl
507
508	!----------------------------------
509	! Heat content - lateral accretion
510	!----------------------------------
511	DO ji = 1, nbpac
512	jl = zcatac(ji) ! categroy in which new ice is put
513	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) ) ) ! zindb=1 if ice =0 otherwise
514	zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb ! old ice thickness
515	zdhex (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) ) ! difference in thickness
516	zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi11 ) ) ! ice totally new in jl category
517	END DO
518
519	DO jk = 1, nlay_i
520	DO ji = 1, nbpac
521	jl = zcatac(ji)
522	zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ]
523	zalphai = MIN( zhice_old(ji,jl) * jk / nlay_i , zh_newice(ji) ) &
524	& - MIN( zhice_old(ji,jl) * ( jk - 1 ) / nlay_i , zh_newice(ji) )
525	ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji) &
526	+ ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / nlay_i &
527	+ za_newice(ji) * ze_newice(ji) * zalphai &
528	+ za_newice(ji) * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( ( zv_i_ac(ji,jl) ) / nlay_i )
529	END DO
530	END DO
531
532	!-----------------------------------------------
533	! Add excessive volume of new ice at the bottom
534	!-----------------------------------------------
535	! If the ice concentration exceeds 1, the remaining volume of new ice
536	! is equally redistributed among all ice categories in which there is
537	! ice
538
539	! Fraction of level ice
540	jm = 1
541	zat_i_lev(:) = 0._wp
542
543	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
544	DO ji = 1, nbpac
545	zat_i_lev(ji) = zat_i_lev(ji) + za_i_ac(ji,jl)
546	END DO
547	END DO
548
549	IF( ln_nicep ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl)
550	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
551	DO ji = 1, nbpac
552	zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) )
553	zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 )
554	END DO
555	END DO
556	IF( ln_nicep ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindx, 1:jpl)
557
558	!---------------------------------
559	! Heat content - bottom accretion
560	!---------------------------------
561	jm = 1
562	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
563	DO ji = 1, nbpac
564	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) ) ) ! zindb=1 if ice =0 otherwise
565	zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
566	zdhicbot (ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb &
567	& + zindb * zdh_frazb(ji) ! frazil ice may coalesce
568	zdummy(ji,jl) = zv_i_ac(ji,jl)/MAX(za_i_ac(ji,jl),epsi10)*zindb ! thickness of residual ice
569	END DO
570	END DO
571
572	! old layers thicknesses and enthalpies
573	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
574	DO jk = 1, nlay_i
575	DO ji = 1, nbpac
576	zthick0(ji,jk,jl) = zhice_old(ji,jl) / nlay_i
577	zqm0 (ji,jk,jl) = ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)
578	END DO
579	END DO
580	END DO
581	!!gm ??? why the previous do loop if ocerwriten by the following one ?
582	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
583	DO ji = 1, nbpac
584	zthick0(ji,nlay_i+1,jl) = zdhicbot(ji,jl)
585	zqm0 (ji,nlay_i+1,jl) = ze_newice(ji) * zdhicbot(ji,jl)
586	END DO ! ji
587	END DO ! jl
588
589	! Redistributing energy on the new grid
590	ze_i_ac(:,:,:) = 0._wp
591	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
592	DO jk = 1, nlay_i
593	DO layer = 1, nlay_i + 1
594	DO ji = 1, nbpac
595	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) )
596	! Redistributing energy on the new grid
597	zweight = MAX ( MIN( zhice_old(ji,jl) * layer , zdummy(ji,jl) * jk ) &
598	& - MAX( zhice_old(ji,jl) * ( layer - 1 ) , zdummy(ji,jl) * ( jk - 1 ) ) , 0._wp ) &
599	& /( MAX(nlay_i * zthick0(ji,layer,jl),epsi10) ) * zindb
600	ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl)
601	END DO ! ji
602	END DO ! layer
603	END DO ! jk
604	END DO ! jl
605
606	DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
607	DO jk = 1, nlay_i
608	DO ji = 1, nbpac
609	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )
610	ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) &
611	& / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * nlay_i * zindb
612	END DO
613	END DO
614	END DO
615
616	!------------
617	! Update age
618	!------------
619	DO jl = 1, jpl
620	DO ji = 1, nbpac
621	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes
622	zoa_i_ac(ji,jl) = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
623	END DO
624	END DO
625
626	!-----------------
627	! Update salinity
628	!-----------------
629	IF( num_sal == 2 .OR. num_sal == 4 ) THEN
630	DO jl = 1, jpl
631	DO ji = 1, nbpac
632	zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes
633	zdv = zv_i_ac(ji,jl) - zv_old(ji,jl)
634	zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb
635	END DO
636	END DO
637	ENDIF
638
639	!------------------------------------------------------------------------------!
640	! 8) Change 2D vectors to 1D vectors
641	!------------------------------------------------------------------------------!
642	DO jl = 1, jpl
643	CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ac (1:nbpac,jl), jpi, jpj )
644	CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj )
645	CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj )
646	IF ( num_sal == 2 .OR. num_sal == 4 ) &
647	CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj )
648	DO jk = 1, nlay_i
649	CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj )
650	END DO
651	END DO
652	CALL tab_1d_2d( nbpac, fseqv , npac(1:nbpac), fseqv_1d (1:nbpac) , jpi, jpj )
653	!
654	ENDIF ! nbpac > 0
655
656	!------------------------------------------------------------------------------!
657	! 9) Change units for e_i
658	!------------------------------------------------------------------------------!
659	DO jl = 1, jpl
660	DO jk = 1, nlay_i ! heat content in 10^9 Joules
661	e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i / unit_fac
662	END DO
663	END DO
664
665	!------------------------------------------------------------------------------\|
666	! 10) Conservation check and changes in each ice category
667	!------------------------------------------------------------------------------\|
668	IF( con_i ) THEN
669	CALL lim_column_sum (jpl, v_i, vt_i_final)
670	fieldid = 'v_i, limthd_lac'
671	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
672	!
673	CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)
674	fieldid = 'e_i, limthd_lac'
675	CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid)
676	!
677	CALL lim_column_sum (jpl, v_s, vt_s_final)
678	fieldid = 'v_s, limthd_lac'
679	CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)
680	!
681	! CALL lim_column_sum (jpl, e_s(:,:,1,:) , et_s_init)
682	! fieldid = 'e_s, limthd_lac'
683	! CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)
684	IF( ln_nicep ) THEN
685	WRITE(numout,*) ' vt_i_init : ', vt_i_init(jiindx,jjindx)
686	WRITE(numout,*) ' vt_i_final: ', vt_i_final(jiindx,jjindx)
687	WRITE(numout,*) ' et_i_init : ', et_i_init(jiindx,jjindx)
688	WRITE(numout,*) ' et_i_final: ', et_i_final(jiindx,jjindx)
689	ENDIF
690	!
691	ENDIF
692	!
693	IF( wrk_not_released(1, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15) .OR. &
694	wrk_not_released(2, 1,2,3,4,5,6,7,8) ) &
695	CALL ctl_stop( 'lim_thd_lac : failed to release workspace arrays' )
696	!
697	END SUBROUTINE lim_thd_lac
698
699	#else
700	!!----------------------------------------------------------------------
701	!! Default option NO LIM3 sea-ice model
702	!!----------------------------------------------------------------------
703	CONTAINS
704	SUBROUTINE lim_thd_lac ! Empty routine
705	END SUBROUTINE lim_thd_lac
706	#endif
707
708	!!======================================================================
709	END MODULE limthd_lac

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

Download in other formats: