New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

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

Context Navigation

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

Last change on this file since 5123 was 5123, checked in by clem, 9 years ago
major LIM3 cleaning + monocat capabilities + NEMO namelist-consistency; sette to follow
Property svn:keywords set to `Id`
File size: 62.8 KB

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

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

Download in other formats: