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

Context Navigation

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 @ 3523

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

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

Download in other formats: