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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 2392

Last change on this file since 2392 was 2287, checked in by smasson, 14 years ago
update licence of all NEMO files...
Property svn:keywords set to `Id`
File size: 74.7 KB

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

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

Download in other formats: