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

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

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

Context Navigation

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

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

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

Download in other formats: