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

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

limitd_me.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4313

Last change on this file since 4313 was 4298, checked in by clem, 11 years ago

Property svn:keywords set to `Id`
File size: 76.3 KB

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

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

Download in other formats: