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

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

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

Context Navigation

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

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

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

Download in other formats: