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

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 2633

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

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

Download in other formats: