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

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

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

Download in other formats: