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/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4899

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

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

Download in other formats: