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

Context Navigation

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 5070

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

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

Download in other formats: