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 branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4656

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

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

Download in other formats: