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

Last change on this file since 5836 was 5836, checked in by cetlod, 9 years ago
merge the simplification branch onto the trunk, see ticket #1612
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_zapsmall
157	CALL lim_var_glo2eqv ! equivalent variables, requested for rafting
158
159	!-----------------------------------------------------------------------------!
160	! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
161	!-----------------------------------------------------------------------------!
162	Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE
163	!
164	CALL lim_itd_me_ridgeprep ! prepare ridging
165	!
166	IF( con_i) CALL lim_column_sum( jpl, v_i, vt_i_init ) ! conservation check
167
168	DO jj = 1, jpj ! Initialize arrays.
169	DO ji = 1, jpi
170	msnow_mlt(ji,jj) = 0._wp
171	esnow_mlt(ji,jj) = 0._wp
172	dardg1dt (ji,jj) = 0._wp
173	dardg2dt (ji,jj) = 0._wp
174	dvirdgdt (ji,jj) = 0._wp
175	opening (ji,jj) = 0._wp
176
177	!-----------------------------------------------------------------------------!
178	! 2) Dynamical inputs (closing rate, divu_adv, opning)
179	!-----------------------------------------------------------------------------!
180	!
181	! 2.1 closing_net
182	!-----------------
183	! Compute the net rate of closing due to convergence
184	! and shear, based on Flato and Hibler (1995).
185	!
186	! The energy dissipation rate is equal to the net closing rate
187	! times the ice strength.
188	!
189	! NOTE: The NET closing rate is equal to the rate that open water
190	! area is removed, plus the rate at which ice area is removed by
191	! ridging, minus the rate at which area is added in new ridges.
192	! The GROSS closing rate is equal to the first two terms (open
193	! water closing and thin ice ridging) without the third term
194	! (thick, newly ridged ice).
195
196	closing_net(ji,jj) = rn_cs * 0.5 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
197
198	! 2.2 divu_adv
199	!--------------
200	! Compute divu_adv, the divergence rate given by the transport/
201	! advection scheme, which may not be equal to divu as computed
202	! from the velocity field.
203	!
204	! If divu_adv < 0, make sure the closing rate is large enough
205	! to give asum = 1.0 after ridging.
206
207	divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep
208
209	IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
210
211	! 2.3 opning
212	!------------
213	! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
214	opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
215	END DO
216	END DO
217
218	!-----------------------------------------------------------------------------!
219	! 3) Ridging iteration
220	!-----------------------------------------------------------------------------!
221	niter = 1 ! iteration counter
222	iterate_ridging = 1
223
224	DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
225
226	DO jj = 1, jpj
227	DO ji = 1, jpi
228
229	! 3.2 closing_gross
230	!-----------------------------------------------------------------------------!
231	! Based on the ITD of ridging and ridged ice, convert the net
232	! closing rate to a gross closing rate.
233	! NOTE: 0 < aksum <= 1
234	closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj)
235
236	! correction to closing rate and opening if closing rate is excessive
237	!---------------------------------------------------------------------
238	! Reduce the closing rate if more than 100% of the open water
239	! would be removed. Reduce the opening rate proportionately.
240	za = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice
241	IF( za > epsi20 ) THEN
242	zfac = MIN( 1._wp, ato_i(ji,jj) / za )
243	closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
244	opning (ji,jj) = opning (ji,jj) * zfac
245	ENDIF
246
247	END DO
248	END DO
249
250	! correction to closing rate / opening if excessive ice removal
251	!---------------------------------------------------------------
252	! Reduce the closing rate if more than 100% of any ice category
253	! would be removed. Reduce the opening rate proportionately.
254	DO jl = 1, jpl
255	DO jj = 1, jpj
256	DO ji = 1, jpi
257	za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
258	IF( za > epsi20 ) THEN
259	zfac = MIN( 1._wp, a_i(ji,jj,jl) / za )
260	closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
261	opning (ji,jj) = opning (ji,jj) * zfac
262	ENDIF
263	END DO
264	END DO
265	END DO
266
267	! 3.3 Redistribute area, volume, and energy.
268	!-----------------------------------------------------------------------------!
269
270	CALL lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
271
272	! 3.4 Compute total area of ice plus open water after ridging.
273	!-----------------------------------------------------------------------------!
274	! This is in general not equal to one because of divergence during transport
275	asum(:,:) = ato_i(:,:)
276	DO jl = 1, jpl
277	asum(:,:) = asum(:,:) + a_i(:,:,jl)
278	END DO
279
280	! 3.5 Do we keep on iterating ???
281	!-----------------------------------------------------------------------------!
282	! Check whether asum = 1. If not (because the closing and opening
283	! rates were reduced above), ridge again with new rates.
284
285	iterate_ridging = 0
286
287	DO jj = 1, jpj
288	DO ji = 1, jpi
289	IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN
290	closing_net(ji,jj) = 0._wp
291	opning (ji,jj) = 0._wp
292	ELSE
293	iterate_ridging = 1
294	divu_adv (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice
295	closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
296	opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) )
297	ENDIF
298	END DO
299	END DO
300
301	IF( lk_mpp ) CALL mpp_max( iterate_ridging )
302
303	! Repeat if necessary.
304	! NOTE: If strength smoothing is turned on, the ridging must be
305	! iterated globally because of the boundary update in the
306	! smoothing.
307
308	niter = niter + 1
309
310	IF( iterate_ridging == 1 ) THEN
311	IF( niter > nitermax ) THEN
312	WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
313	WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
314	ENDIF
315	CALL lim_itd_me_ridgeprep
316	ENDIF
317
318	END DO !! on the do while over iter
319
320	!-----------------------------------------------------------------------------!
321	! 4) Ridging diagnostics
322	!-----------------------------------------------------------------------------!
323	! Convert ridging rate diagnostics to correct units.
324	! Update fresh water and heat fluxes due to snow melt.
325	DO jj = 1, jpj
326	DO ji = 1, jpi
327
328	dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice
329	dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice
330	dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice
331	opening (ji,jj) = opening (ji,jj) * r1_rdtice
332
333	!-----------------------------------------------------------------------------!
334	! 5) Heat, salt and freshwater fluxes
335	!-----------------------------------------------------------------------------!
336	wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean
337	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice ! heat sink for ocean (<0, W.m-2)
338
339	END DO
340	END DO
341
342	! Check if there is a ridging error
343	IF( lwp ) THEN
344	DO jj = 1, jpj
345	DO ji = 1, jpi
346	IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug
347	WRITE(numout,*) ' '
348	WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
349	WRITE(numout,*) ' limitd_me '
350	WRITE(numout,*) ' POINT : ', ji, jj
351	WRITE(numout,*) ' jpl, a_i, athorn '
352	WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
353	DO jl = 1, jpl
354	WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
355	END DO
356	ENDIF
357	END DO
358	END DO
359	END IF
360
361	! Conservation check
362	IF ( con_i ) THEN
363	CALL lim_column_sum (jpl, v_i, vt_i_final)
364	fieldid = ' v_i : limitd_me '
365	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
366	ENDIF
367
368	CALL lim_var_agg( 1 )
369
370	!-----------------------------------------------------------------------------!
371	! control prints
372	!-----------------------------------------------------------------------------!
373	IF(ln_ctl) THEN
374	CALL lim_var_glo2eqv
375
376	CALL prt_ctl_info(' ')
377	CALL prt_ctl_info(' - Cell values : ')
378	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
379	CALL prt_ctl(tab2d_1=e1e2t, clinfo1=' lim_itd_me : cell area :')
380	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :')
381	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :')
382	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me : vt_s :')
383	DO jl = 1, jpl
384	CALL prt_ctl_info(' ')
385	CALL prt_ctl_info(' - Category : ', ivar1=jl)
386	CALL prt_ctl_info(' ~~~~~~~~~~')
387	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_me : a_i : ')
388	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_me : ht_i : ')
389	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_me : ht_s : ')
390	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_me : v_i : ')
391	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_me : v_s : ')
392	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_me : e_s : ')
393	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_me : t_su : ')
394	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_me : t_snow : ')
395	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_me : sm_i : ')
396	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_me : smv_i : ')
397	DO jk = 1, nlay_i
398	CALL prt_ctl_info(' ')
399	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
400	CALL prt_ctl_info(' ~~~~~~~')
401	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : t_i : ')
402	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : e_i : ')
403	END DO
404	END DO
405	ENDIF
406
407	! conservation test
408	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
409
410	ENDIF ! ln_limdyn=.true.
411	!
412	CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
413	!
414	IF( nn_timing == 1 ) CALL timing_stop('limitd_me')
415	END SUBROUTINE lim_itd_me
416
417
418	SUBROUTINE lim_itd_me_icestrength( kstrngth )
419	!!----------------------------------------------------------------------
420	!! * ROUTINE lim_itd_me_icestrength *
421	!!
422	!! ** Purpose : computes ice strength used in dynamics routines of ice thickness
423	!!
424	!! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2)
425	!! dissipated per unit area removed from the ice pack under compression,
426	!! and assumed proportional to the change in potential energy caused
427	!! by ridging. Note that only Hibler's formulation is stable and that
428	!! ice strength has to be smoothed
429	!!
430	!! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
431	!!----------------------------------------------------------------------
432	INTEGER, INTENT(in) :: kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979)
433	INTEGER :: ji,jj, jl ! dummy loop indices
434	INTEGER :: ksmooth ! smoothing the resistance to deformation
435	INTEGER :: numts_rm ! number of time steps for the P smoothing
436	REAL(wp) :: zhi, zp, z1_3 ! local scalars
437	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
438	!!----------------------------------------------------------------------
439
440	CALL wrk_alloc( jpi, jpj, zworka )
441
442	!------------------------------------------------------------------------------!
443	! 1) Initialize
444	!------------------------------------------------------------------------------!
445	strength(:,:) = 0._wp
446
447	!------------------------------------------------------------------------------!
448	! 2) Compute thickness distribution of ridging and ridged ice
449	!------------------------------------------------------------------------------!
450	CALL lim_itd_me_ridgeprep
451
452	!------------------------------------------------------------------------------!
453	! 3) Rothrock(1975)'s method
454	!------------------------------------------------------------------------------!
455	IF( kstrngth == 1 ) THEN
456	z1_3 = 1._wp / 3._wp
457	DO jl = 1, jpl
458	DO jj= 1, jpj
459	DO ji = 1, jpi
460	!
461	IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN
462	zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
463	!----------------------------
464	! PE loss from deforming ice
465	!----------------------------
466	strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * zhi * zhi
467
468	!--------------------------
469	! PE gain from rafting ice
470	!--------------------------
471	strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * zhi * zhi
472
473	!----------------------------
474	! PE gain from ridging ice
475	!----------------------------
476	strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) / krdg(ji,jj,jl) &
477	* z1_3 * ( hrmax(ji,jj,jl)2 + hrmin(ji,jj,jl)2 + hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )
478	!!(a3-b3)/(a-b) = aa+ab+bb
479	ENDIF
480	!
481	END DO
482	END DO
483	END DO
484
485	strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:)
486	! where Cp = (g/2)(rhow-rhoi)(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
487	ksmooth = 1
488
489	!------------------------------------------------------------------------------!
490	! 4) Hibler (1979)' method
491	!------------------------------------------------------------------------------!
492	ELSE ! kstrngth ne 1: Hibler (1979) form
493	!
494	strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) )
495	!
496	ksmooth = 1
497	!
498	ENDIF ! kstrngth
499
500	!
501	!------------------------------------------------------------------------------!
502	! 5) Impact of brine volume
503	!------------------------------------------------------------------------------!
504	! CAN BE REMOVED
505	!
506	IF( ln_icestr_bvf ) THEN
507
508	DO jj = 1, jpj
509	DO ji = 1, jpi
510	strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
511	END DO
512	END DO
513
514	ENDIF
515
516	!
517	!------------------------------------------------------------------------------!
518	! 6) Smoothing ice strength
519	!------------------------------------------------------------------------------!
520	!
521	!-------------------
522	! Spatial smoothing
523	!-------------------
524	IF ( ksmooth == 1 ) THEN
525
526	CALL lbc_lnk( strength, 'T', 1. )
527
528	DO jj = 2, jpjm1
529	DO ji = 2, jpim1
530	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN
531	zworka(ji,jj) = ( 4.0 * strength(ji,jj) &
532	& + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &
533	& + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
534	& ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
535	ELSE
536	zworka(ji,jj) = 0._wp
537	ENDIF
538	END DO
539	END DO
540
541	DO jj = 2, jpjm1
542	DO ji = 2, jpim1
543	strength(ji,jj) = zworka(ji,jj)
544	END DO
545	END DO
546	CALL lbc_lnk( strength, 'T', 1. )
547
548	ENDIF
549
550	!--------------------
551	! Temporal smoothing
552	!--------------------
553	IF ( numit == nit000 + nn_fsbc - 1 ) THEN
554	strp1(:,:) = 0.0
555	strp2(:,:) = 0.0
556	ENDIF
557
558	IF ( ksmooth == 2 ) THEN
559
560
561	CALL lbc_lnk( strength, 'T', 1. )
562
563	DO jj = 1, jpj - 1
564	DO ji = 1, jpi - 1
565	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN
566	numts_rm = 1 ! number of time steps for the running mean
567	IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
568	IF ( strp2(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
569	zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
570	strp2(ji,jj) = strp1(ji,jj)
571	strp1(ji,jj) = strength(ji,jj)
572	strength(ji,jj) = zp
573
574	ENDIF
575	END DO
576	END DO
577
578	ENDIF ! ksmooth
579
580	CALL lbc_lnk( strength, 'T', 1. ) ! Boundary conditions
581
582	CALL wrk_dealloc( jpi, jpj, zworka )
583	!
584	END SUBROUTINE lim_itd_me_icestrength
585
586
587	SUBROUTINE lim_itd_me_ridgeprep
588	!!---------------------------------------------------------------------!
589	!! * ROUTINE lim_itd_me_ridgeprep *
590	!!
591	!! ** Purpose : preparation for ridging and strength calculations
592	!!
593	!! ** Method : Compute the thickness distribution of the ice and open water
594	!! participating in ridging and of the resulting ridges.
595	!!---------------------------------------------------------------------!
596	INTEGER :: ji,jj, jl ! dummy loop indices
597	REAL(wp) :: Gstari, astari, zhi, hrmean, zdummy ! local scalar
598	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
599	REAL(wp), POINTER, DIMENSION(:,:,:) :: Gsum ! Gsum(n) = sum of areas in categories 0 to n
600	!------------------------------------------------------------------------------!
601
602	CALL wrk_alloc( jpi,jpj, zworka )
603	CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
604
605	Gstari = 1.0/rn_gstar
606	astari = 1.0/rn_astar
607	aksum(:,:) = 0.0
608	athorn(:,:,:) = 0.0
609	aridge(:,:,:) = 0.0
610	araft (:,:,:) = 0.0
611	hrmin(:,:,:) = 0.0
612	hrmax(:,:,:) = 0.0
613	hraft(:,:,:) = 0.0
614	krdg (:,:,:) = 1.0
615
616	! ! Zero out categories with very small areas
617	CALL lim_var_zapsmall
618
619	!------------------------------------------------------------------------------!
620	! 1) Participation function
621	!------------------------------------------------------------------------------!
622
623	! Compute total area of ice plus open water.
624	! This is in general not equal to one because of divergence during transport
625	asum(:,:) = ato_i(:,:)
626	DO jl = 1, jpl
627	asum(:,:) = asum(:,:) + a_i(:,:,jl)
628	END DO
629
630	! Compute cumulative thickness distribution function
631	! Compute the cumulative thickness distribution function Gsum,
632	! where Gsum(n) is the fractional area in categories 0 to n.
633	! initial value (in h = 0) equals open water area
634
635	Gsum(:,:,-1) = 0._wp
636	Gsum(:,:,0 ) = ato_i(:,:)
637
638	! for each value of h, you have to add ice concentration then
639	DO jl = 1, jpl
640	Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl)
641	END DO
642
643	! Normalize the cumulative distribution to 1
644	zworka(:,:) = 1._wp / Gsum(:,:,jpl)
645	DO jl = 0, jpl
646	Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:)
647	END DO
648
649	! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
650	!--------------------------------------------------------------------------------------------------
651	! Compute the participation function athorn; this is analogous to
652	! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
653	! area lost from category n due to ridging/closing
654	! athorn(n) = total area lost due to ridging/closing
655	! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
656	!
657	! The expressions for athorn are found by integrating b(h)g(h) between
658	! the category boundaries.
659	!-----------------------------------------------------------------
660
661	IF( nn_partfun == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975)
662	DO jl = 0, jpl
663	DO jj = 1, jpj
664	DO ji = 1, jpi
665	IF( Gsum(ji,jj,jl) < rn_gstar) THEN
666	athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * &
667	& ( 2.0 - (Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari )
668	ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN
669	athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) * &
670	& ( 2.0 - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari )
671	ELSE
672	athorn(ji,jj,jl) = 0.0
673	ENDIF
674	END DO
675	END DO
676	END DO
677
678	ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007)
679	!
680	zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array
681	DO jl = -1, jpl
682	Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
683	END DO
684	DO jl = 0, jpl
685	athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
686	END DO
687	!
688	ENDIF
689
690	IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions
691	!
692	DO jl = 1, jpl
693	DO jj = 1, jpj
694	DO ji = 1, jpi
695	IF ( athorn(ji,jj,jl) > 0._wp ) THEN
696	!!gm TANH( -X ) = - TANH( X ) so can be computed only 1 time....
697	aridge(ji,jj,jl) = ( TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
698	araft (ji,jj,jl) = ( TANH ( -rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
699	IF ( araft(ji,jj,jl) < epsi06 ) araft(ji,jj,jl) = 0._wp
700	aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 )
701	ENDIF
702	END DO
703	END DO
704	END DO
705
706	ELSE
707	!
708	DO jl = 1, jpl
709	aridge(:,:,jl) = athorn(:,:,jl)
710	END DO
711	!
712	ENDIF
713
714	IF( ln_rafting ) THEN
715
716	IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 .AND. lwp ) THEN
717	DO jl = 1, jpl
718	DO jj = 1, jpj
719	DO ji = 1, jpi
720	IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN
721	WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
722	WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
723	WRITE(numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj)
724	WRITE(numout,*) ' aridge : ', aridge(ji,jj,1:jpl)
725	WRITE(numout,*) ' araft : ', araft(ji,jj,1:jpl)
726	WRITE(numout,*) ' athorn : ', athorn(ji,jj,1:jpl)
727	ENDIF
728	END DO
729	END DO
730	END DO
731	ENDIF
732
733	ENDIF
734
735	!-----------------------------------------------------------------
736	! 2) Transfer function
737	!-----------------------------------------------------------------
738	! Compute max and min ridged ice thickness for each ridging category.
739	! Assume ridged ice is uniformly distributed between hrmin and hrmax.
740	!
741	! This parameterization is a modified version of Hibler (1980).
742	! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
743	! and for very thick ridging ice must be >= krdgmin*hi
744	!
745	! The minimum ridging thickness, hrmin, is equal to 2*hi
746	! (i.e., rafting) and for very thick ridging ice is
747	! constrained by hrmin <= (hrmean + hi)/2.
748	!
749	! The maximum ridging thickness, hrmax, is determined by
750	! hrmean and hrmin.
751	!
752	! These modifications have the effect of reducing the ice strength
753	! (relative to the Hibler formulation) when very thick ice is
754	! ridging.
755	!
756	! aksum = net area removed/ total area removed
757	! where total area removed = area of ice that ridges
758	! net area removed = total area removed - area of new ridges
759	!-----------------------------------------------------------------
760
761	! Transfer function
762	DO jl = 1, jpl !all categories have a specific transfer function
763	DO jj = 1, jpj
764	DO ji = 1, jpi
765
766	IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN
767	zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
768	hrmean = MAX(SQRT(rn_hstarzhi), zhikrdgmin)
769	hrmin(ji,jj,jl) = MIN(2.0zhi, 0.5(hrmean + zhi))
770	hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)
771	hraft(ji,jj,jl) = kraft*zhi
772	krdg(ji,jj,jl) = hrmean / zhi
773	ELSE
774	hraft(ji,jj,jl) = 0.0
775	hrmin(ji,jj,jl) = 0.0
776	hrmax(ji,jj,jl) = 0.0
777	krdg (ji,jj,jl) = 1.0
778	ENDIF
779
780	END DO
781	END DO
782	END DO
783
784	! Normalization factor : aksum, ensures mass conservation
785	aksum(:,:) = athorn(:,:,0)
786	DO jl = 1, jpl
787	aksum(:,:) = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) ) &
788	& + araft (:,:,jl) * ( 1._wp - 1._wp / kraft )
789	END DO
790	!
791	CALL wrk_dealloc( jpi,jpj, zworka )
792	CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
793	!
794	END SUBROUTINE lim_itd_me_ridgeprep
795
796
797	SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
798	!!----------------------------------------------------------------------
799	!! * ROUTINE lim_itd_me_icestrength *
800	!!
801	!! ** Purpose : shift ridging ice among thickness categories of ice thickness
802	!!
803	!! ** Method : Remove area, volume, and energy from each ridging category
804	!! and add to thicker ice categories.
805	!!----------------------------------------------------------------------
806	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear
807	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area removed, excluding area of new ridges
808	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: msnow_mlt ! mass of snow added to ocean (kg m-2)
809	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)
810	!
811	CHARACTER (len=80) :: fieldid ! field identifier
812	LOGICAL, PARAMETER :: l_conservation_check = .true. ! if true, check conservation (useful for debugging)
813	!
814	INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices
815	INTEGER :: ij ! horizontal index, combines i and j loops
816	INTEGER :: icells ! number of cells with aicen > puny
817	REAL(wp) :: hL, hR, farea, ztmelts ! left and right limits of integration
818
819	INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices
820
821	REAL(wp), POINTER, DIMENSION(:,:) :: vice_init, vice_final ! ice volume summed over categories
822	REAL(wp), POINTER, DIMENSION(:,:) :: eice_init, eice_final ! ice energy summed over layers
823
824	REAL(wp), POINTER, DIMENSION(:,:,:) :: aicen_init, vicen_init ! ice area & volume before ridging
825	REAL(wp), POINTER, DIMENSION(:,:,:) :: vsnwn_init, esnwn_init ! snow volume & energy before ridging
826	REAL(wp), POINTER, DIMENSION(:,:,:) :: smv_i_init, oa_i_init ! ice salinity & age before ridging
827
828	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: eicen_init ! ice energy before ridging
829
830	REAL(wp), POINTER, DIMENSION(:,:) :: afrac , fvol ! fraction of category area ridged & new ridge volume going to n2
831	REAL(wp), POINTER, DIMENSION(:,:) :: ardg1 , ardg2 ! area of ice ridged & new ridges
832	REAL(wp), POINTER, DIMENSION(:,:) :: vsrdg , esrdg ! snow volume & energy of ridging 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	REAL(wp), POINTER, DIMENSION(:,:) :: oirdg1, oirdg2 ! ice age of ice ridged
842
843	REAL(wp), POINTER, DIMENSION(:,:) :: afrft ! fraction of category area rafted
844	REAL(wp), POINTER, DIMENSION(:,:) :: arft1 , arft2 ! area of ice rafted and new rafted zone
845	REAL(wp), POINTER, DIMENSION(:,:) :: virft , vsrft ! ice & snow volume of rafting ice
846	REAL(wp), POINTER, DIMENSION(:,:) :: esrft , smrft ! snow energy & salinity of rafting ice
847	REAL(wp), POINTER, DIMENSION(:,:) :: oirft1, oirft2 ! ice age of ice rafted
848
849	REAL(wp), POINTER, DIMENSION(:,:,:) :: eirft ! ice energy of rafting ice
850	REAL(wp), POINTER, DIMENSION(:,:,:) :: erdg1 ! enth*volume of ice ridged
851	REAL(wp), POINTER, DIMENSION(:,:,:) :: erdg2 ! enth*volume of new ridges
852	REAL(wp), POINTER, DIMENSION(:,:,:) :: ersw ! enth of water trapped into ridges
853	!!----------------------------------------------------------------------
854
855	CALL wrk_alloc( (jpi+1)*(jpj+1), indxi, indxj )
856	CALL wrk_alloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )
857	CALL wrk_alloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
858	CALL wrk_alloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 )
859	CALL wrk_alloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
860	CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
861	CALL wrk_alloc( jpi, jpj, nlay_i, eirft, erdg1, erdg2, ersw )
862	CALL wrk_alloc( jpi, jpj, nlay_i, jpl, eicen_init )
863
864	! Conservation check
865	eice_init(:,:) = 0._wp
866
867	IF( con_i ) THEN
868	CALL lim_column_sum (jpl, v_i, vice_init )
869	CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_init )
870	DO ji = mi0(iiceprt), mi1(iiceprt)
871	DO jj = mj0(jiceprt), mj1(jiceprt)
872	WRITE(numout,*) ' vice_init : ', vice_init(ji,jj)
873	WRITE(numout,*) ' eice_init : ', eice_init(ji,jj)
874	END DO
875	END DO
876	ENDIF
877
878	!-------------------------------------------------------------------------------
879	! 1) Compute change in open water area due to closing and opening.
880	!-------------------------------------------------------------------------------
881	DO jj = 1, jpj
882	DO ji = 1, jpi
883	ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice &
884	& + opning(ji,jj) * rdt_ice
885	IF ( ato_i(ji,jj) < -epsi10 ) THEN ! there is a bug
886	IF(lwp) WRITE(numout,*) 'Ridging error: ato_i < 0 -- ato_i : ',ato_i(ji,jj)
887	ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error
888	ato_i(ji,jj) = 0._wp
889	ENDIF
890	END DO
891	END DO
892
893	!-----------------------------------------------------------------
894	! 2) Save initial state variables
895	!-----------------------------------------------------------------
896	aicen_init(:,:,:) = a_i (:,:,:)
897	vicen_init(:,:,:) = v_i (:,:,:)
898	vsnwn_init(:,:,:) = v_s (:,:,:)
899	smv_i_init(:,:,:) = smv_i(:,:,:)
900	esnwn_init(:,:,:) = e_s (:,:,1,:)
901	eicen_init(:,:,:,:) = e_i (:,:,:,:)
902	oa_i_init (:,:,:) = oa_i (:,:,:)
903
904	!
905	!-----------------------------------------------------------------
906	! 3) Pump everything from ice which is being ridged / rafted
907	!-----------------------------------------------------------------
908	! Compute the area, volume, and energy of ice ridging in each
909	! category, along with the area of the resulting ridge.
910
911	DO jl1 = 1, jpl !jl1 describes the ridging category
912
913	!------------------------------------------------
914	! 3.1) Identify grid cells with nonzero ridging
915	!------------------------------------------------
916
917	icells = 0
918	DO jj = 1, jpj
919	DO ji = 1, jpi
920	IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp &
921	& .AND. closing_gross(ji,jj) > 0._wp ) THEN
922	icells = icells + 1
923	indxi(icells) = ji
924	indxj(icells) = jj
925	ENDIF
926	END DO
927	END DO
928
929	DO ij = 1, icells
930	ji = indxi(ij)
931	jj = indxj(ij)
932
933	!--------------------------------------------------------------------
934	! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
935	!--------------------------------------------------------------------
936
937	ardg1(ji,jj) = aridge(ji,jj,jl1)closing_gross(ji,jj)rdt_ice
938	arft1(ji,jj) = araft (ji,jj,jl1)closing_gross(ji,jj)rdt_ice
939	ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1)
940	arft2(ji,jj) = arft1(ji,jj) / kraft
941
942	!---------------------------------------------------------------
943	! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
944	!---------------------------------------------------------------
945
946	afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging
947	afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting
948
949	IF( afrac(ji,jj) > kamax + epsi10 ) THEN ! there is a bug
950	IF(lwp) WRITE(numout,*) ' ardg > a_i -- ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1)
951	ELSEIF( afrac(ji,jj) > kamax ) THEN ! roundoff error
952	afrac(ji,jj) = kamax
953	ENDIF
954
955	IF( afrft(ji,jj) > kamax + epsi10 ) THEN ! there is a bug
956	IF(lwp) WRITE(numout,*) ' arft > a_i -- arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)
957	ELSEIF( afrft(ji,jj) > kamax) THEN ! roundoff error
958	afrft(ji,jj) = kamax
959	ENDIF
960
961	!--------------------------------------------------------------------------
962	! 3.4) Subtract area, volume, and energy from ridging
963	! / rafting category n1.
964	!--------------------------------------------------------------------------
965	vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj)
966	vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + rn_por_rdg )
967	vsw (ji,jj) = vrdg1(ji,jj) * rn_por_rdg
968
969	vsrdg (ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj)
970	esrdg (ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj)
971	srdg1 (ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
972	oirdg1(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj)
973	oirdg2(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj) / krdg(ji,jj,jl1)
974
975	! rafting volumes, heat contents ...
976	virft (ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)
977	vsrft (ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj)
978	esrft (ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj)
979	smrft (ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)
980	oirft1(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj)
981	oirft2(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj) / kraft
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	smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1 (ji,jj) - smrft (ji,jj)
989	oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1) - oirdg1(ji,jj) - oirft1(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
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, dhr, dhr2 )
1190	CALL wrk_dealloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 )
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, eirft, erdg1, erdg2, ersw )
1194	CALL wrk_dealloc( jpi, jpj, nlay_i, 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: