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.

limitd_me.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 5167

Last change on this file since 5167 was 5167, checked in by clem, 9 years ago
LIM3 bug fix. see ticket #1492 (forthcoming update) which also solve ticket #1497
Property svn:keywords set to `Id`
File size: 59.5 KB

Line
1	MODULE limitd_me
2	!!======================================================================
3	!! * MODULE limitd_me *
4	!! LIM-3 : Mechanical impact on ice thickness distribution
5	!!======================================================================
6	!! History : LIM ! 2006-02 (M. Vancoppenolle) Original code
7	!! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_dyn
8	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
9	!!----------------------------------------------------------------------
10	#if defined key_lim3
11	!!----------------------------------------------------------------------
12	!! 'key_lim3' LIM-3 sea-ice model
13	!!----------------------------------------------------------------------
14	USE par_oce ! ocean parameters
15	USE dom_oce ! ocean domain
16	USE phycst ! physical constants (ocean directory)
17	USE sbc_oce ! surface boundary condition: ocean fields
18	USE thd_ice ! LIM thermodynamics
19	USE ice ! LIM variables
20	USE dom_ice ! LIM domain
21	USE limvar ! LIM
22	USE lbclnk ! lateral boundary condition - MPP exchanges
23	USE lib_mpp ! MPP library
24	USE wrk_nemo ! work arrays
25	USE prtctl ! Print control
26
27	USE in_out_manager ! I/O manager
28	USE iom ! I/O manager
29	USE lib_fortran ! glob_sum
30	USE limdiahsb
31	USE timing ! Timing
32	USE limcons ! conservation tests
33
34	IMPLICIT NONE
35	PRIVATE
36
37	PUBLIC lim_itd_me ! called by ice_stp
38	PUBLIC lim_itd_me_icestrength
39	PUBLIC lim_itd_me_init
40	PUBLIC lim_itd_me_alloc ! called by sbc_lim_init
41
42	!-----------------------------------------------------------------------
43	! Variables shared among ridging subroutines
44	!-----------------------------------------------------------------------
45	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: asum ! sum of total ice and open water area
46	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged
47	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/
48	! ! closing associated w/ category n
49	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness
50	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness
51	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice
52	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: krdg ! mean ridge thickness/thickness of ridging ice
53	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: aridge ! participating ice ridging
54	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: araft ! participating ice rafting
55
56	REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier
57	REAL(wp), PARAMETER :: kraft = 2.0_wp ! rafting multipliyer
58	REAL(wp), PARAMETER :: kamax = 1.0_wp ! max of ice area authorized (clem: scheme is not stable if kamax <= 0.99)
59
60	REAL(wp) :: Cp !
61	!
62	!-----------------------------------------------------------------------
63	! Ridging diagnostic arrays for history files
64	!-----------------------------------------------------------------------
65	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dardg1dt ! rate of fractional area loss by ridging ice (1/s)
66	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dardg2dt ! rate of fractional area gain by new ridges (1/s)
67	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dvirdgdt ! rate of ice volume ridged (m/s)
68	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: opening ! rate of opening due to divergence/shear (1/s)
69	!
70	!!----------------------------------------------------------------------
71	!! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
72	!! $Id$
73	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
74	!!----------------------------------------------------------------------
75	CONTAINS
76
77	INTEGER FUNCTION lim_itd_me_alloc()
78	!!---------------------------------------------------------------------!
79	!! * ROUTINE lim_itd_me_alloc *
80	!!---------------------------------------------------------------------!
81	ALLOCATE( &
82	!* Variables shared among ridging subroutines
83	& asum (jpi,jpj) , athorn(jpi,jpj,0:jpl) , &
84	& aksum(jpi,jpj) , &
85	!
86	& hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , &
87	& hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , &
88	!
89	!* Ridging diagnostic arrays for history files
90	& dardg1dt(jpi,jpj) , dardg2dt(jpi,jpj) , &
91	& dvirdgdt(jpi,jpj) , opening(jpi,jpj) , STAT=lim_itd_me_alloc )
92	!
93	IF( lim_itd_me_alloc /= 0 ) CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' )
94	!
95	END FUNCTION lim_itd_me_alloc
96
97
98	SUBROUTINE lim_itd_me
99	!!---------------------------------------------------------------------!
100	!! * ROUTINE lim_itd_me *
101	!!
102	!! ** Purpose : computes the mechanical redistribution of ice thickness
103	!!
104	!! ** Method : Steps :
105	!! 1) Thickness categories boundaries, ice / o.w. concentrations
106	!! Ridge preparation
107	!! 2) Dynamical inputs (closing rate, divu_adv, opning)
108	!! 3) Ridging iteration
109	!! 4) Ridging diagnostics
110	!! 5) Heat, salt and freshwater fluxes
111	!! 6) Compute increments of tate variables and come back to old values
112	!!
113	!! References : Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
114	!! Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
115	!! Rothrock, D. A., 1975: JGR, 80, 4514-4519.
116	!! Thorndike et al., 1975, JGR, 80, 4501-4513.
117	!! Bitz et al., JGR, 2001
118	!! Amundrud and Melling, JGR 2005
119	!! Babko et al., JGR 2002
120	!!
121	!! This routine is based on CICE code and authors William H. Lipscomb,
122	!! and Elizabeth C. Hunke, LANL are gratefully acknowledged
123	!!--------------------------------------------------------------------!
124	INTEGER :: ji, jj, jk, jl ! dummy loop index
125	INTEGER :: niter ! local integer
126	INTEGER :: iterate_ridging ! if true, repeat the ridging
127	REAL(wp) :: za, zfac ! local scalar
128	CHARACTER (len = 15) :: fieldid
129	REAL(wp), POINTER, DIMENSION(:,:) :: closing_net ! net rate at which area is removed (1/s)
130	! (ridging ice area - area of new ridges) / dt
131	REAL(wp), POINTER, DIMENSION(:,:) :: divu_adv ! divu as implied by transport scheme (1/s)
132	REAL(wp), POINTER, DIMENSION(:,:) :: opning ! rate of opening due to divergence/shear
133	REAL(wp), POINTER, DIMENSION(:,:) :: closing_gross ! rate at which area removed, not counting area of new ridges
134	REAL(wp), POINTER, DIMENSION(:,:) :: msnow_mlt ! mass of snow added to ocean (kg m-2)
135	REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)
136	REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories
137	!
138	INTEGER, PARAMETER :: nitermax = 20
139	!
140	REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
141	!!-----------------------------------------------------------------------------
142	IF( nn_timing == 1 ) CALL timing_start('limitd_me')
143
144	CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
145
146	IF(ln_ctl) THEN
147	CALL prt_ctl(tab2d_1=ato_i , clinfo1=' lim_itd_me: ato_i : ', tab2d_2=at_i , clinfo2=' at_i : ')
148	CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ')
149	ENDIF
150
151	IF( ln_limdyn ) THEN ! Start ridging and rafting !
152
153	! conservation test
154	IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
155
156	CALL lim_var_glo2eqv ! equivalent variables, requested for rafting
157	!-----------------------------------------------------------------------------!
158	! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
159	!-----------------------------------------------------------------------------!
160	Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE
161	!
162	CALL lim_itd_me_ridgeprep ! prepare ridging
163	!
164	IF( con_i) CALL lim_column_sum( jpl, v_i, vt_i_init ) ! conservation check
165
166	DO jj = 1, jpj ! Initialize arrays.
167	DO ji = 1, jpi
168	msnow_mlt(ji,jj) = 0._wp
169	esnow_mlt(ji,jj) = 0._wp
170	dardg1dt (ji,jj) = 0._wp
171	dardg2dt (ji,jj) = 0._wp
172	dvirdgdt (ji,jj) = 0._wp
173	opening (ji,jj) = 0._wp
174
175	!-----------------------------------------------------------------------------!
176	! 2) Dynamical inputs (closing rate, divu_adv, opning)
177	!-----------------------------------------------------------------------------!
178	!
179	! 2.1 closing_net
180	!-----------------
181	! Compute the net rate of closing due to convergence
182	! and shear, based on Flato and Hibler (1995).
183	!
184	! The energy dissipation rate is equal to the net closing rate
185	! times the ice strength.
186	!
187	! NOTE: The NET closing rate is equal to the rate that open water
188	! area is removed, plus the rate at which ice area is removed by
189	! ridging, minus the rate at which area is added in new ridges.
190	! The GROSS closing rate is equal to the first two terms (open
191	! water closing and thin ice ridging) without the third term
192	! (thick, newly ridged ice).
193
194	closing_net(ji,jj) = rn_cs * 0.5 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
195
196	! 2.2 divu_adv
197	!--------------
198	! Compute divu_adv, the divergence rate given by the transport/
199	! advection scheme, which may not be equal to divu as computed
200	! from the velocity field.
201	!
202	! If divu_adv < 0, make sure the closing rate is large enough
203	! to give asum = 1.0 after ridging.
204
205	divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep
206
207	IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
208
209	! 2.3 opning
210	!------------
211	! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
212	opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
213	END DO
214	END DO
215
216	!-----------------------------------------------------------------------------!
217	! 3) Ridging iteration
218	!-----------------------------------------------------------------------------!
219	niter = 1 ! iteration counter
220	iterate_ridging = 1
221
222	DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
223
224	DO jj = 1, jpj
225	DO ji = 1, jpi
226
227	! 3.2 closing_gross
228	!-----------------------------------------------------------------------------!
229	! Based on the ITD of ridging and ridged ice, convert the net
230	! closing rate to a gross closing rate.
231	! NOTE: 0 < aksum <= 1
232	closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj)
233
234	! correction to closing rate and opening if closing rate is excessive
235	!---------------------------------------------------------------------
236	! Reduce the closing rate if more than 100% of the open water
237	! would be removed. Reduce the opening rate proportionately.
238	za = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice
239	IF( za > epsi20 ) THEN
240	zfac = MIN( 1._wp, ato_i(ji,jj) / za )
241	closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
242	opning (ji,jj) = opning (ji,jj) * zfac
243	ENDIF
244
245	END DO
246	END DO
247
248	! correction to closing rate / opening if excessive ice removal
249	!---------------------------------------------------------------
250	! Reduce the closing rate if more than 100% of any ice category
251	! would be removed. Reduce the opening rate proportionately.
252	DO jl = 1, jpl
253	DO jj = 1, jpj
254	DO ji = 1, jpi
255	za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
256	IF( za > epsi20 ) THEN
257	zfac = MIN( 1._wp, a_i(ji,jj,jl) / za )
258	closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
259	opning (ji,jj) = opning (ji,jj) * zfac
260	ENDIF
261	END DO
262	END DO
263	END DO
264
265	! 3.3 Redistribute area, volume, and energy.
266	!-----------------------------------------------------------------------------!
267
268	CALL lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
269
270	! 3.4 Compute total area of ice plus open water after ridging.
271	!-----------------------------------------------------------------------------!
272	! This is in general not equal to one because of divergence during transport
273	asum(:,:) = ato_i(:,:)
274	DO jl = 1, jpl
275	asum(:,:) = asum(:,:) + a_i(:,:,jl)
276	END DO
277
278	! 3.5 Do we keep on iterating ???
279	!-----------------------------------------------------------------------------!
280	! Check whether asum = 1. If not (because the closing and opening
281	! rates were reduced above), ridge again with new rates.
282
283	iterate_ridging = 0
284
285	DO jj = 1, jpj
286	DO ji = 1, jpi
287	IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN
288	closing_net(ji,jj) = 0._wp
289	opning (ji,jj) = 0._wp
290	ELSE
291	iterate_ridging = 1
292	divu_adv (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice
293	closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
294	opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) )
295	ENDIF
296	END DO
297	END DO
298
299	IF( lk_mpp ) CALL mpp_max( iterate_ridging )
300
301	! Repeat if necessary.
302	! NOTE: If strength smoothing is turned on, the ridging must be
303	! iterated globally because of the boundary update in the
304	! smoothing.
305
306	niter = niter + 1
307
308	IF( iterate_ridging == 1 ) THEN
309	IF( niter > nitermax ) THEN
310	WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
311	WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
312	ENDIF
313	CALL lim_itd_me_ridgeprep
314	ENDIF
315
316	END DO !! on the do while over iter
317
318	!-----------------------------------------------------------------------------!
319	! 4) Ridging diagnostics
320	!-----------------------------------------------------------------------------!
321	! Convert ridging rate diagnostics to correct units.
322	! Update fresh water and heat fluxes due to snow melt.
323	DO jj = 1, jpj
324	DO ji = 1, jpi
325
326	dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice
327	dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice
328	dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice
329	opening (ji,jj) = opening (ji,jj) * r1_rdtice
330
331	!-----------------------------------------------------------------------------!
332	! 5) Heat, salt and freshwater fluxes
333	!-----------------------------------------------------------------------------!
334	wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean
335	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice ! heat sink for ocean (<0, W.m-2)
336
337	END DO
338	END DO
339
340	! Check if there is a ridging error
341	IF( lwp ) THEN
342	DO jj = 1, jpj
343	DO ji = 1, jpi
344	IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug
345	WRITE(numout,*) ' '
346	WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
347	WRITE(numout,*) ' limitd_me '
348	WRITE(numout,*) ' POINT : ', ji, jj
349	WRITE(numout,*) ' jpl, a_i, athorn '
350	WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
351	DO jl = 1, jpl
352	WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
353	END DO
354	ENDIF
355	END DO
356	END DO
357	END IF
358
359	! Conservation check
360	IF ( con_i ) THEN
361	CALL lim_column_sum (jpl, v_i, vt_i_final)
362	fieldid = ' v_i : limitd_me '
363	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
364	ENDIF
365
366	! updates
367	CALL lim_var_agg( 1 )
368
369	!-----------------------------------------------------------------------------!
370	! control prints
371	!-----------------------------------------------------------------------------!
372	IF(ln_ctl) THEN
373	CALL lim_var_glo2eqv
374
375	CALL prt_ctl_info(' ')
376	CALL prt_ctl_info(' - Cell values : ')
377	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
378	CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_me : cell area :')
379	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :')
380	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :')
381	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me : vt_s :')
382	DO jl = 1, jpl
383	CALL prt_ctl_info(' ')
384	CALL prt_ctl_info(' - Category : ', ivar1=jl)
385	CALL prt_ctl_info(' ~~~~~~~~~~')
386	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_me : a_i : ')
387	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_me : ht_i : ')
388	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_me : ht_s : ')
389	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_me : v_i : ')
390	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_me : v_s : ')
391	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_me : e_s : ')
392	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_me : t_su : ')
393	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_me : t_snow : ')
394	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_me : sm_i : ')
395	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_me : smv_i : ')
396	DO jk = 1, nlay_i
397	CALL prt_ctl_info(' ')
398	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
399	CALL prt_ctl_info(' ~~~~~~~')
400	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : t_i : ')
401	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : e_i : ')
402	END DO
403	END DO
404	ENDIF
405
406	! conservation test
407	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
408
409	ENDIF ! ln_limdyn=.true.
410	!
411	CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
412	!
413	IF( nn_timing == 1 ) CALL timing_stop('limitd_me')
414	END SUBROUTINE lim_itd_me
415
416
417	SUBROUTINE lim_itd_me_icestrength( kstrngth )
418	!!----------------------------------------------------------------------
419	!! * ROUTINE lim_itd_me_icestrength *
420	!!
421	!! ** Purpose : computes ice strength used in dynamics routines of ice thickness
422	!!
423	!! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2)
424	!! dissipated per unit area removed from the ice pack under compression,
425	!! and assumed proportional to the change in potential energy caused
426	!! by ridging. Note that only Hibler's formulation is stable and that
427	!! ice strength has to be smoothed
428	!!
429	!! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
430	!!----------------------------------------------------------------------
431	INTEGER, INTENT(in) :: kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979)
432	INTEGER :: ji,jj, jl ! dummy loop indices
433	INTEGER :: ksmooth ! smoothing the resistance to deformation
434	INTEGER :: numts_rm ! number of time steps for the P smoothing
435	REAL(wp) :: zhi, zp, z1_3 ! local scalars
436	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
437	!!----------------------------------------------------------------------
438
439	CALL wrk_alloc( jpi, jpj, zworka )
440
441	!------------------------------------------------------------------------------!
442	! 1) Initialize
443	!------------------------------------------------------------------------------!
444	strength(:,:) = 0._wp
445
446	!------------------------------------------------------------------------------!
447	! 2) Compute thickness distribution of ridging and ridged ice
448	!------------------------------------------------------------------------------!
449	CALL lim_itd_me_ridgeprep
450
451	!------------------------------------------------------------------------------!
452	! 3) Rothrock(1975)'s method
453	!------------------------------------------------------------------------------!
454	IF( kstrngth == 1 ) THEN
455	z1_3 = 1._wp / 3._wp
456	DO jl = 1, jpl
457	DO jj= 1, jpj
458	DO ji = 1, jpi
459	!
460	IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN
461	zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
462	!----------------------------
463	! PE loss from deforming ice
464	!----------------------------
465	strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * zhi * zhi
466
467	!--------------------------
468	! PE gain from rafting ice
469	!--------------------------
470	strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * zhi * zhi
471
472	!----------------------------
473	! PE gain from ridging ice
474	!----------------------------
475	strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) / krdg(ji,jj,jl) &
476	* z1_3 * ( hrmax(ji,jj,jl)2 + hrmin(ji,jj,jl)2 + hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )
477	!!(a3-b3)/(a-b) = aa+ab+bb
478	ENDIF
479	!
480	END DO
481	END DO
482	END DO
483
484	strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:)
485	! where Cp = (g/2)(rhow-rhoi)(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
486	ksmooth = 1
487
488	!------------------------------------------------------------------------------!
489	! 4) Hibler (1979)' method
490	!------------------------------------------------------------------------------!
491	ELSE ! kstrngth ne 1: Hibler (1979) form
492	!
493	strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) )
494	!
495	ksmooth = 1
496	!
497	ENDIF ! kstrngth
498
499	!
500	!------------------------------------------------------------------------------!
501	! 5) Impact of brine volume
502	!------------------------------------------------------------------------------!
503	! CAN BE REMOVED
504	!
505	IF( ln_icestr_bvf ) THEN
506
507	DO jj = 1, jpj
508	DO ji = 1, jpi
509	strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
510	END DO
511	END DO
512
513	ENDIF
514
515	!
516	!------------------------------------------------------------------------------!
517	! 6) Smoothing ice strength
518	!------------------------------------------------------------------------------!
519	!
520	!-------------------
521	! Spatial smoothing
522	!-------------------
523	IF ( ksmooth == 1 ) THEN
524
525	CALL lbc_lnk( strength, 'T', 1. )
526
527	DO jj = 2, jpjm1
528	DO ji = 2, jpim1
529	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN
530	zworka(ji,jj) = ( 4.0 * strength(ji,jj) &
531	& + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &
532	& + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
533	& ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
534	ELSE
535	zworka(ji,jj) = 0._wp
536	ENDIF
537	END DO
538	END DO
539
540	DO jj = 2, jpjm1
541	DO ji = 2, jpim1
542	strength(ji,jj) = zworka(ji,jj)
543	END DO
544	END DO
545	CALL lbc_lnk( strength, 'T', 1. )
546
547	ENDIF
548
549	!--------------------
550	! Temporal smoothing
551	!--------------------
552	IF ( numit == nit000 + nn_fsbc - 1 ) THEN
553	strp1(:,:) = 0.0
554	strp2(:,:) = 0.0
555	ENDIF
556
557	IF ( ksmooth == 2 ) THEN
558
559
560	CALL lbc_lnk( strength, 'T', 1. )
561
562	DO jj = 1, jpj - 1
563	DO ji = 1, jpi - 1
564	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN
565	numts_rm = 1 ! number of time steps for the running mean
566	IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
567	IF ( strp2(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
568	zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
569	strp2(ji,jj) = strp1(ji,jj)
570	strp1(ji,jj) = strength(ji,jj)
571	strength(ji,jj) = zp
572
573	ENDIF
574	END DO
575	END DO
576
577	ENDIF ! ksmooth
578
579	CALL lbc_lnk( strength, 'T', 1. ) ! Boundary conditions
580
581	CALL wrk_dealloc( jpi, jpj, zworka )
582	!
583	END SUBROUTINE lim_itd_me_icestrength
584
585
586	SUBROUTINE lim_itd_me_ridgeprep
587	!!---------------------------------------------------------------------!
588	!! * ROUTINE lim_itd_me_ridgeprep *
589	!!
590	!! ** Purpose : preparation for ridging and strength calculations
591	!!
592	!! ** Method : Compute the thickness distribution of the ice and open water
593	!! participating in ridging and of the resulting ridges.
594	!!---------------------------------------------------------------------!
595	INTEGER :: ji,jj, jl ! dummy loop indices
596	REAL(wp) :: Gstari, astari, zhi, hrmean, zdummy ! local scalar
597	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
598	REAL(wp), POINTER, DIMENSION(:,:,:) :: Gsum ! Gsum(n) = sum of areas in categories 0 to n
599	!------------------------------------------------------------------------------!
600
601	CALL wrk_alloc( jpi,jpj, zworka )
602	CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
603
604	Gstari = 1.0/rn_gstar
605	astari = 1.0/rn_astar
606	aksum(:,:) = 0.0
607	athorn(:,:,:) = 0.0
608	aridge(:,:,:) = 0.0
609	araft (:,:,:) = 0.0
610	hrmin(:,:,:) = 0.0
611	hrmax(:,:,:) = 0.0
612	hraft(:,:,:) = 0.0
613	krdg (:,:,:) = 1.0
614
615	! ! Zero out categories with very small areas
616	CALL lim_var_zapsmall
617
618	!------------------------------------------------------------------------------!
619	! 1) Participation function
620	!------------------------------------------------------------------------------!
621
622	! Compute total area of ice plus open water.
623	! This is in general not equal to one because of divergence during transport
624	asum(:,:) = ato_i(:,:)
625	DO jl = 1, jpl
626	asum(:,:) = asum(:,:) + a_i(:,:,jl)
627	END DO
628
629	! Compute cumulative thickness distribution function
630	! Compute the cumulative thickness distribution function Gsum,
631	! where Gsum(n) is the fractional area in categories 0 to n.
632	! initial value (in h = 0) equals open water area
633
634	Gsum(:,:,-1) = 0._wp
635	Gsum(:,:,0 ) = ato_i(:,:)
636
637	! for each value of h, you have to add ice concentration then
638	DO jl = 1, jpl
639	Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl)
640	END DO
641
642	! Normalize the cumulative distribution to 1
643	zworka(:,:) = 1._wp / Gsum(:,:,jpl)
644	DO jl = 0, jpl
645	Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:)
646	END DO
647
648	! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
649	!--------------------------------------------------------------------------------------------------
650	! Compute the participation function athorn; this is analogous to
651	! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
652	! area lost from category n due to ridging/closing
653	! athorn(n) = total area lost due to ridging/closing
654	! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
655	!
656	! The expressions for athorn are found by integrating b(h)g(h) between
657	! the category boundaries.
658	!-----------------------------------------------------------------
659
660	IF( nn_partfun == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975)
661	DO jl = 0, jpl
662	DO jj = 1, jpj
663	DO ji = 1, jpi
664	IF( Gsum(ji,jj,jl) < rn_gstar) THEN
665	athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * &
666	& ( 2.0 - (Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari )
667	ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN
668	athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) * &
669	& ( 2.0 - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari )
670	ELSE
671	athorn(ji,jj,jl) = 0.0
672	ENDIF
673	END DO
674	END DO
675	END DO
676
677	ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007)
678	!
679	zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array
680	DO jl = -1, jpl
681	Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
682	END DO
683	DO jl = 0, jpl
684	athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
685	END DO
686	!
687	ENDIF
688
689	IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions
690	!
691	DO jl = 1, jpl
692	DO jj = 1, jpj
693	DO ji = 1, jpi
694	IF ( athorn(ji,jj,jl) > 0._wp ) THEN
695	!!gm TANH( -X ) = - TANH( X ) so can be computed only 1 time....
696	aridge(ji,jj,jl) = ( TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
697	araft (ji,jj,jl) = ( TANH ( -rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
698	IF ( araft(ji,jj,jl) < epsi06 ) araft(ji,jj,jl) = 0._wp
699	aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 )
700	ENDIF
701	END DO
702	END DO
703	END DO
704
705	ELSE
706	!
707	DO jl = 1, jpl
708	aridge(:,:,jl) = athorn(:,:,jl)
709	END DO
710	!
711	ENDIF
712
713	IF( ln_rafting ) THEN
714
715	IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 .AND. lwp ) THEN
716	DO jl = 1, jpl
717	DO jj = 1, jpj
718	DO ji = 1, jpi
719	IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN
720	WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
721	WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
722	WRITE(numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj)
723	WRITE(numout,*) ' aridge : ', aridge(ji,jj,1:jpl)
724	WRITE(numout,*) ' araft : ', araft(ji,jj,1:jpl)
725	WRITE(numout,*) ' athorn : ', athorn(ji,jj,1:jpl)
726	ENDIF
727	END DO
728	END DO
729	END DO
730	ENDIF
731
732	ENDIF
733
734	!-----------------------------------------------------------------
735	! 2) Transfer function
736	!-----------------------------------------------------------------
737	! Compute max and min ridged ice thickness for each ridging category.
738	! Assume ridged ice is uniformly distributed between hrmin and hrmax.
739	!
740	! This parameterization is a modified version of Hibler (1980).
741	! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
742	! and for very thick ridging ice must be >= krdgmin*hi
743	!
744	! The minimum ridging thickness, hrmin, is equal to 2*hi
745	! (i.e., rafting) and for very thick ridging ice is
746	! constrained by hrmin <= (hrmean + hi)/2.
747	!
748	! The maximum ridging thickness, hrmax, is determined by
749	! hrmean and hrmin.
750	!
751	! These modifications have the effect of reducing the ice strength
752	! (relative to the Hibler formulation) when very thick ice is
753	! ridging.
754	!
755	! aksum = net area removed/ total area removed
756	! where total area removed = area of ice that ridges
757	! net area removed = total area removed - area of new ridges
758	!-----------------------------------------------------------------
759
760	! Transfer function
761	DO jl = 1, jpl !all categories have a specific transfer function
762	DO jj = 1, jpj
763	DO ji = 1, jpi
764
765	IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN
766	zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
767	hrmean = MAX(SQRT(rn_hstarzhi), zhikrdgmin)
768	hrmin(ji,jj,jl) = MIN(2.0zhi, 0.5(hrmean + zhi))
769	hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)
770	hraft(ji,jj,jl) = kraft*zhi
771	krdg(ji,jj,jl) = hrmean / zhi
772	ELSE
773	hraft(ji,jj,jl) = 0.0
774	hrmin(ji,jj,jl) = 0.0
775	hrmax(ji,jj,jl) = 0.0
776	krdg (ji,jj,jl) = 1.0
777	ENDIF
778
779	END DO
780	END DO
781	END DO
782
783	! Normalization factor : aksum, ensures mass conservation
784	aksum(:,:) = athorn(:,:,0)
785	DO jl = 1, jpl
786	aksum(:,:) = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) ) &
787	& + araft (:,:,jl) * ( 1._wp - 1._wp / kraft )
788	END DO
789	!
790	CALL wrk_dealloc( jpi,jpj, zworka )
791	CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
792	!
793	END SUBROUTINE lim_itd_me_ridgeprep
794
795
796	SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
797	!!----------------------------------------------------------------------
798	!! * ROUTINE lim_itd_me_icestrength *
799	!!
800	!! ** Purpose : shift ridging ice among thickness categories of ice thickness
801	!!
802	!! ** Method : Remove area, volume, and energy from each ridging category
803	!! and add to thicker ice categories.
804	!!----------------------------------------------------------------------
805	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear
806	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area removed, excluding area of new ridges
807	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: msnow_mlt ! mass of snow added to ocean (kg m-2)
808	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)
809	!
810	CHARACTER (len=80) :: fieldid ! field identifier
811	LOGICAL, PARAMETER :: l_conservation_check = .true. ! if true, check conservation (useful for debugging)
812	!
813	INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices
814	INTEGER :: ij ! horizontal index, combines i and j loops
815	INTEGER :: icells ! number of cells with aicen > puny
816	REAL(wp) :: hL, hR, farea, ztmelts ! left and right limits of integration
817
818	INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices
819
820	REAL(wp), POINTER, DIMENSION(:,:) :: vice_init, vice_final ! ice volume summed over categories
821	REAL(wp), POINTER, DIMENSION(:,:) :: eice_init, eice_final ! ice energy summed over layers
822
823	REAL(wp), POINTER, DIMENSION(:,:,:) :: aicen_init, vicen_init ! ice area & volume before ridging
824	REAL(wp), POINTER, DIMENSION(:,:,:) :: vsnwn_init, esnwn_init ! snow volume & energy before ridging
825	REAL(wp), POINTER, DIMENSION(:,:,:) :: smv_i_init, oa_i_init ! ice salinity & age before ridging
826
827	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: eicen_init ! ice energy before ridging
828
829	REAL(wp), POINTER, DIMENSION(:,:) :: afrac , fvol ! fraction of category area ridged & new ridge volume going to n2
830	REAL(wp), POINTER, DIMENSION(:,:) :: ardg1 , ardg2 ! area of ice ridged & new ridges
831	REAL(wp), POINTER, DIMENSION(:,:) :: vsrdg , esrdg ! snow volume & energy of ridging ice
832	REAL(wp), POINTER, DIMENSION(:,:) :: oirdg1, oirdg2 ! areal age content of ridged & rifging ice
833	REAL(wp), POINTER, DIMENSION(:,:) :: dhr , dhr2 ! hrmax - hrmin & hrmax^2 - hrmin^2
834
835	REAL(wp), POINTER, DIMENSION(:,:) :: vrdg1 ! volume of ice ridged
836	REAL(wp), POINTER, DIMENSION(:,:) :: vrdg2 ! volume of new ridges
837	REAL(wp), POINTER, DIMENSION(:,:) :: vsw ! volume of seawater trapped into ridges
838	REAL(wp), POINTER, DIMENSION(:,:) :: srdg1 ! sal*volume of ice ridged
839	REAL(wp), POINTER, DIMENSION(:,:) :: srdg2 ! sal*volume of new ridges
840	REAL(wp), POINTER, DIMENSION(:,:) :: smsw ! sal*volume of water trapped into ridges
841
842	REAL(wp), POINTER, DIMENSION(:,:) :: afrft ! fraction of category area rafted
843	REAL(wp), POINTER, DIMENSION(:,:) :: arft1 , arft2 ! area of ice rafted and new rafted zone
844	REAL(wp), POINTER, DIMENSION(:,:) :: virft , vsrft ! ice & snow volume of rafting ice
845	REAL(wp), POINTER, DIMENSION(:,:) :: esrft , smrft ! snow energy & salinity of rafting ice
846	REAL(wp), POINTER, DIMENSION(:,:) :: oirft1, oirft2 ! areal age content of rafted ice & rafting ice
847
848	REAL(wp), POINTER, DIMENSION(:,:,:) :: eirft ! ice energy of rafting ice
849	REAL(wp), POINTER, DIMENSION(:,:,:) :: erdg1 ! enth*volume of ice ridged
850	REAL(wp), POINTER, DIMENSION(:,:,:) :: erdg2 ! enth*volume of new ridges
851	REAL(wp), POINTER, DIMENSION(:,:,:) :: ersw ! enth of water trapped into ridges
852	!!----------------------------------------------------------------------
853
854	CALL wrk_alloc( (jpi+1)*(jpj+1), indxi, indxj )
855	CALL wrk_alloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )
856	CALL wrk_alloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
857	CALL wrk_alloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw )
858	CALL wrk_alloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
859	CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
860	CALL wrk_alloc( jpi, jpj, nlay_i+1, eirft, erdg1, erdg2, ersw )
861	CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, eicen_init )
862
863	! Conservation check
864	eice_init(:,:) = 0._wp
865
866	IF( con_i ) THEN
867	CALL lim_column_sum (jpl, v_i, vice_init )
868	CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_init )
869	DO ji = mi0(iiceprt), mi1(iiceprt)
870	DO jj = mj0(jiceprt), mj1(jiceprt)
871	WRITE(numout,*) ' vice_init : ', vice_init(ji,jj)
872	WRITE(numout,*) ' eice_init : ', eice_init(ji,jj)
873	END DO
874	END DO
875	ENDIF
876
877	!-------------------------------------------------------------------------------
878	! 1) Compute change in open water area due to closing and opening.
879	!-------------------------------------------------------------------------------
880	DO jj = 1, jpj
881	DO ji = 1, jpi
882	ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice &
883	& + opning(ji,jj) * rdt_ice
884	IF ( ato_i(ji,jj) < -epsi10 ) THEN ! there is a bug
885	IF(lwp) WRITE(numout,*) 'Ridging error: ato_i < 0 -- ato_i : ',ato_i(ji,jj)
886	ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error
887	ato_i(ji,jj) = 0._wp
888	ENDIF
889	END DO
890	END DO
891
892	!-----------------------------------------------------------------
893	! 2) Save initial state variables
894	!-----------------------------------------------------------------
895	aicen_init(:,:,:) = a_i (:,:,:)
896	vicen_init(:,:,:) = v_i (:,:,:)
897	vsnwn_init(:,:,:) = v_s (:,:,:)
898	smv_i_init(:,:,:) = smv_i(:,:,:)
899	oa_i_init (:,:,:) = oa_i (:,:,:)
900	esnwn_init(:,:,:) = e_s (:,:,1,:)
901	eicen_init(:,:,:,:) = e_i (:,:,:,:)
902
903	!
904	!-----------------------------------------------------------------
905	! 3) Pump everything from ice which is being ridged / rafted
906	!-----------------------------------------------------------------
907	! Compute the area, volume, and energy of ice ridging in each
908	! category, along with the area of the resulting ridge.
909
910	DO jl1 = 1, jpl !jl1 describes the ridging category
911
912	!------------------------------------------------
913	! 3.1) Identify grid cells with nonzero ridging
914	!------------------------------------------------
915
916	icells = 0
917	DO jj = 1, jpj
918	DO ji = 1, jpi
919	IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp &
920	& .AND. closing_gross(ji,jj) > 0._wp ) THEN
921	icells = icells + 1
922	indxi(icells) = ji
923	indxj(icells) = jj
924	ENDIF
925	END DO
926	END DO
927
928	DO ij = 1, icells
929	ji = indxi(ij)
930	jj = indxj(ij)
931
932	!--------------------------------------------------------------------
933	! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
934	!--------------------------------------------------------------------
935
936	ardg1(ji,jj) = aridge(ji,jj,jl1)closing_gross(ji,jj)rdt_ice
937	arft1(ji,jj) = araft (ji,jj,jl1)closing_gross(ji,jj)rdt_ice
938	ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1)
939	arft2(ji,jj) = arft1(ji,jj) / kraft
940
941	oirdg1(ji,jj)= aridge(ji,jj,jl1)closing_gross(ji,jj)rdt_ice
942	oirft1(ji,jj)= araft (ji,jj,jl1)closing_gross(ji,jj)rdt_ice
943	oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1)
944	oirft2(ji,jj)= oirft1(ji,jj) / kraft
945
946	!---------------------------------------------------------------
947	! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
948	!---------------------------------------------------------------
949
950	afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging
951	afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting
952
953	IF( afrac(ji,jj) > kamax + epsi10 ) THEN ! there is a bug
954	IF(lwp) WRITE(numout,*) ' ardg > a_i -- ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1)
955	ELSEIF( afrac(ji,jj) > kamax ) THEN ! roundoff error
956	afrac(ji,jj) = kamax
957	ENDIF
958
959	IF( afrft(ji,jj) > kamax + epsi10 ) THEN ! there is a bug
960	IF(lwp) WRITE(numout,*) ' arft > a_i -- arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)
961	ELSEIF( afrft(ji,jj) > kamax) THEN ! roundoff error
962	afrft(ji,jj) = kamax
963	ENDIF
964
965	!--------------------------------------------------------------------------
966	! 3.4) Subtract area, volume, and energy from ridging
967	! / rafting category n1.
968	!--------------------------------------------------------------------------
969	vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj)
970	vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + rn_por_rdg )
971	vsw (ji,jj) = vrdg1(ji,jj) * rn_por_rdg
972
973	vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj)
974	esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj)
975	srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
976
977	! rafting volumes, heat contents ...
978	virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)
979	vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj)
980	esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj)
981	smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)
982
983	! substract everything
984	a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - ardg1(ji,jj) - arft1(ji,jj)
985	v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - vrdg1(ji,jj) - virft(ji,jj)
986	v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - vsrdg(ji,jj) - vsrft(ji,jj)
987	e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj) - esrft(ji,jj)
988	oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1) - oirdg1(ji,jj) - oirft1(ji,jj)
989	smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj) - smrft(ji,jj)
990
991	!-----------------------------------------------------------------
992	! 3.5) Compute properties of new ridges
993	!-----------------------------------------------------------------
994	!-------------
995	! Salinity
996	!-------------
997	smsw(ji,jj) = vsw(ji,jj) * sss_m(ji,jj) ! salt content of seawater frozen in voids !! MV HC2014
998	srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj) ! salt content of new ridge
999
1000	!srdg2(ji,jj) = MIN( rn_simax * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity
1001
1002	sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice
1003	wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice ! increase in ice volume du to seawater frozen in voids
1004
1005	!------------------------------------
1006	! 3.6 Increment ridging diagnostics
1007	!------------------------------------
1008
1009	! jl1 looping 1-jpl
1010	! ij looping 1-icells
1011
1012	dardg1dt (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)
1013	dardg2dt (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)
1014	opening (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice
1015
1016	IF( con_i ) vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj)
1017
1018	!------------------------------------------
1019	! 3.7 Put the snow somewhere in the ocean
1020	!------------------------------------------
1021	! Place part of the snow lost by ridging into the ocean.
1022	! Note that esnow_mlt < 0; the ocean must cool to melt snow.
1023	! If the ocean temp = Tf already, new ice must grow.
1024	! During the next time step, thermo_rates will determine whether
1025	! the ocean cools or new ice grows.
1026	! jl1 looping 1-jpl
1027	! ij looping 1-icells
1028
1029	msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosnvsrdg(ji,jj)(1.0-rn_fsnowrdg) & ! rafting included
1030	& + rhosnvsrft(ji,jj)(1.0-rn_fsnowrft)
1031
1032	! in J/m2 (same as e_s)
1033	esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-rn_fsnowrdg) & !rafting included
1034	& - esrft(ji,jj)*(1.0-rn_fsnowrft)
1035
1036	!-----------------------------------------------------------------
1037	! 3.8 Compute quantities used to apportion ice among categories
1038	! in the n2 loop below
1039	!-----------------------------------------------------------------
1040
1041	! jl1 looping 1-jpl
1042	! ij looping 1-icells
1043
1044	dhr (ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1)
1045	dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1)
1046
1047	END DO
1048
1049	!--------------------------------------------------------------------
1050	! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
1051	! compute ridged ice enthalpy
1052	!--------------------------------------------------------------------
1053	DO jk = 1, nlay_i
1054	DO ij = 1, icells
1055	ji = indxi(ij)
1056	jj = indxj(ij)
1057	! heat content of ridged ice
1058	erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj)
1059	eirft(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)
1060	e_i (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk)
1061
1062
1063	! enthalpy of the trapped seawater (J/m2, >0)
1064	! clem: if sst>0, then ersw <0 (is that possible?)
1065	ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * sst_m(ji,jj) * r1_nlay_i
1066
1067	! heat flux to the ocean
1068	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux
1069
1070	! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
1071	erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
1072
1073	END DO
1074	END DO
1075
1076
1077	IF( con_i ) THEN
1078	DO jk = 1, nlay_i
1079	DO ij = 1, icells
1080	ji = indxi(ij)
1081	jj = indxj(ij)
1082	eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk)
1083	END DO
1084	END DO
1085	ENDIF
1086
1087	!-------------------------------------------------------------------------------
1088	! 4) Add area, volume, and energy of new ridge to each category jl2
1089	!-------------------------------------------------------------------------------
1090	! jl1 looping 1-jpl
1091	DO jl2 = 1, jpl
1092	! over categories to which ridged ice is transferred
1093	DO ij = 1, icells
1094	ji = indxi(ij)
1095	jj = indxj(ij)
1096
1097	! Compute the fraction of ridged ice area and volume going to
1098	! thickness category jl2.
1099	! Transfer area, volume, and energy accordingly.
1100
1101	IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN
1102	hL = 0._wp
1103	hR = 0._wp
1104	ELSE
1105	hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
1106	hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2) )
1107	ENDIF
1108
1109	! fraction of ridged ice area and volume going to n2
1110	farea = ( hR - hL ) / dhr(ji,jj)
1111	fvol(ji,jj) = ( hRhR - hLhL ) / dhr2(ji,jj)
1112
1113	a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ardg2 (ji,jj) * farea
1114	v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj)
1115	v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg
1116	e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg
1117	smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + srdg2 (ji,jj) * fvol(ji,jj)
1118	oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirdg2(ji,jj) * farea
1119
1120	END DO ! ij
1121
1122	! Transfer ice energy to category jl2 by ridging
1123	DO jk = 1, nlay_i
1124	DO ij = 1, icells
1125	ji = indxi(ij)
1126	jj = indxj(ij)
1127	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) * erdg2(ji,jj,jk)
1128	END DO
1129	END DO
1130	!
1131	END DO ! jl2 (new ridges)
1132
1133	DO jl2 = 1, jpl
1134
1135	DO ij = 1, icells
1136	ji = indxi(ij)
1137	jj = indxj(ij)
1138	! Compute the fraction of rafted ice area and volume going to
1139	! thickness category jl2, transfer area, volume, and energy accordingly.
1140	!
1141	IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN
1142	a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + arft2 (ji,jj)
1143	v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + virft (ji,jj)
1144	v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrft (ji,jj) * rn_fsnowrft
1145	e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrft (ji,jj) * rn_fsnowrft
1146	smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + smrft (ji,jj)
1147	oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirft2(ji,jj)
1148	ENDIF
1149	!
1150	END DO
1151
1152	! Transfer rafted ice energy to category jl2
1153	DO jk = 1, nlay_i
1154	DO ij = 1, icells
1155	ji = indxi(ij)
1156	jj = indxj(ij)
1157	IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN
1158	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk)
1159	ENDIF
1160	END DO
1161	END DO
1162
1163	END DO
1164
1165	END DO ! jl1 (deforming categories)
1166
1167	! Conservation check
1168	IF ( con_i ) THEN
1169	CALL lim_column_sum (jpl, v_i, vice_final)
1170	fieldid = ' v_i : limitd_me '
1171	CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid)
1172
1173	CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_final )
1174	fieldid = ' e_i : limitd_me '
1175	CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid)
1176
1177	DO ji = mi0(iiceprt), mi1(iiceprt)
1178	DO jj = mj0(jiceprt), mj1(jiceprt)
1179	WRITE(numout,*) ' vice_init : ', vice_init (ji,jj)
1180	WRITE(numout,*) ' vice_final : ', vice_final(ji,jj)
1181	WRITE(numout,*) ' eice_init : ', eice_init (ji,jj)
1182	WRITE(numout,*) ' eice_final : ', eice_final(ji,jj)
1183	END DO
1184	END DO
1185	ENDIF
1186	!
1187	CALL wrk_dealloc( (jpi+1)*(jpj+1), indxi, indxj )
1188	CALL wrk_dealloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )
1189	CALL wrk_dealloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
1190	CALL wrk_dealloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw )
1191	CALL wrk_dealloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
1192	CALL wrk_dealloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
1193	CALL wrk_dealloc( jpi, jpj, nlay_i+1, eirft, erdg1, erdg2, ersw )
1194	CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, eicen_init )
1195	!
1196	END SUBROUTINE lim_itd_me_ridgeshift
1197
1198	SUBROUTINE lim_itd_me_init
1199	!!-------------------------------------------------------------------
1200	!! * ROUTINE lim_itd_me_init *
1201	!!
1202	!! ** Purpose : Physical constants and parameters linked
1203	!! to the mechanical ice redistribution
1204	!!
1205	!! ** Method : Read the namiceitdme namelist
1206	!! and check the parameters values
1207	!! called at the first timestep (nit000)
1208	!!
1209	!! ** input : Namelist namiceitdme
1210	!!-------------------------------------------------------------------
1211	INTEGER :: ios ! Local integer output status for namelist read
1212	NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft, &
1213	& rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, &
1214	& nn_partfun
1215	!!-------------------------------------------------------------------
1216	!
1217	REWIND( numnam_ice_ref ) ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
1218	READ ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
1219	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
1220
1221	REWIND( numnam_ice_cfg ) ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
1222	READ ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
1223	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
1224	IF(lwm) WRITE ( numoni, namiceitdme )
1225	!
1226	IF (lwp) THEN ! control print
1227	WRITE(numout,*)
1228	WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
1229	WRITE(numout,*)' ~~~~~~~~~~~~~~~'
1230	WRITE(numout,*)' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs
1231	WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg
1232	WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft
1233	WRITE(numout,*)' Fraction of total ice coverage contributing to ridging rn_gstar = ', rn_gstar
1234	WRITE(numout,)' Equivalent to G for an exponential part function rn_astar = ', rn_astar
1235	WRITE(numout,*)' Quantity playing a role in max ridged ice thickness rn_hstar = ', rn_hstar
1236	WRITE(numout,*)' Rafting of ice sheets or not ln_rafting = ', ln_rafting
1237	WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft
1238	WRITE(numout,*)' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft
1239	WRITE(numout,*)' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg
1240	WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun
1241	ENDIF
1242	!
1243	END SUBROUTINE lim_itd_me_init
1244
1245	#else
1246	!!----------------------------------------------------------------------
1247	!! Default option Empty module NO LIM-3 sea-ice model
1248	!!----------------------------------------------------------------------
1249	CONTAINS
1250	SUBROUTINE lim_itd_me ! Empty routines
1251	END SUBROUTINE lim_itd_me
1252	SUBROUTINE lim_itd_me_icestrength
1253	END SUBROUTINE lim_itd_me_icestrength
1254	SUBROUTINE lim_itd_me_init
1255	END SUBROUTINE lim_itd_me_init
1256	#endif
1257	!!======================================================================
1258	END MODULE limitd_me

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

Download in other formats: