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

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

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

Download in other formats: