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

Context Navigation

source: trunk/NEMO/LIM_SRC_3/limitd_me.F90 @ 1410

Last change on this file since 1410 was 1156, checked in by rblod, 16 years ago
Update Id and licence information, see ticket #210
Property svn:keywords set to `Id`
File size: 74.8 KB

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

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

Download in other formats: