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

Context Navigation

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 7560

Last change on this file since 7560 was 5080, checked in by clem, 9 years ago
LIM3 cleaning again (almost done)
Property svn:keywords set to `Id`
File size: 62.8 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 limthd_lac ! LIM
	22	USE limvar ! LIM
	23	USE in_out_manager ! I/O manager
	24	USE lbclnk ! lateral boundary condition - MPP exchanges
	25	USE lib_mpp ! MPP library
	26	USE wrk_nemo ! work arrays
	27	USE prtctl ! Print control
[5051]	28
[4161]	29	USE iom ! I/O manager
[4688]	30	USE lib_fortran ! glob_sum
[4161]	31	USE limdiahsb
[4688]	32	USE timing ! Timing
	33	USE limcons ! conservation tests
[921]	34
[825]	35	IMPLICIT NONE
	36	PRIVATE
	37
[2715]	38	PUBLIC lim_itd_me ! called by ice_stp
	39	PUBLIC lim_itd_me_icestrength
	40	PUBLIC lim_itd_me_init
[5051]	41	PUBLIC lim_itd_me_alloc ! called by sbc_lim_init
[825]	42
[921]	43	!-----------------------------------------------------------------------
	44	! Variables shared among ridging subroutines
	45	!-----------------------------------------------------------------------
[4333]	46	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: asum ! sum of total ice and open water area
	47	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged
	48	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/
	49	! ! closing associated w/ category n
	50	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness
	51	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness
	52	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice
	53	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: krdg ! mean ridge thickness/thickness of ridging ice
	54	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: aridge ! participating ice ridging
	55	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: araft ! participating ice rafting
[825]	56
[2715]	57	REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier
	58	REAL(wp), PARAMETER :: kraft = 2.0_wp ! rafting multipliyer
[4333]	59	REAL(wp), PARAMETER :: kamax = 1.0_wp ! max of ice area authorized (clem: scheme is not stable if kamax <= 0.99)
[825]	60
[2715]	61	REAL(wp) :: Cp !
[921]	62	!
	63	!-----------------------------------------------------------------------
	64	! Ridging diagnostic arrays for history files
	65	!-----------------------------------------------------------------------
[2715]	66	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dardg1dt ! rate of fractional area loss by ridging ice (1/s)
	67	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dardg2dt ! rate of fractional area gain by new ridges (1/s)
	68	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dvirdgdt ! rate of ice volume ridged (m/s)
	69	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: opening ! rate of opening due to divergence/shear (1/s)
[3625]	70	!
[825]	71	!!----------------------------------------------------------------------
[2528]	72	!! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
[1156]	73	!! $Id$
[2715]	74	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	75	!!----------------------------------------------------------------------
[2715]	76	CONTAINS
[825]	77
[2715]	78	INTEGER FUNCTION lim_itd_me_alloc()
	79	!!---------------------------------------------------------------------!
	80	!! * ROUTINE lim_itd_me_alloc *
	81	!!---------------------------------------------------------------------!
	82	ALLOCATE( &
	83	!* Variables shared among ridging subroutines
	84	& asum (jpi,jpj) , athorn(jpi,jpj,0:jpl) , &
	85	& aksum(jpi,jpj) , &
	86	!
	87	& hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , &
	88	& hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , &
	89	!
	90	!* Ridging diagnostic arrays for history files
	91	& dardg1dt(jpi,jpj) , dardg2dt(jpi,jpj) , &
	92	& dvirdgdt(jpi,jpj) , opening(jpi,jpj) , STAT=lim_itd_me_alloc )
	93	!
	94	IF( lim_itd_me_alloc /= 0 ) CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' )
	95	!
	96	END FUNCTION lim_itd_me_alloc
[1156]	97
[825]	98
[2715]	99	SUBROUTINE lim_itd_me
[921]	100	!!---------------------------------------------------------------------!
	101	!! * ROUTINE lim_itd_me *
	102	!!
[2715]	103	!! ** Purpose : computes the mechanical redistribution of ice thickness
[921]	104	!!
[2715]	105	!! ** Method : Steps :
	106	!! 1) Thickness categories boundaries, ice / o.w. concentrations
	107	!! Ridge preparation
	108	!! 2) Dynamical inputs (closing rate, divu_adv, opning)
	109	!! 3) Ridging iteration
	110	!! 4) Ridging diagnostics
	111	!! 5) Heat, salt and freshwater fluxes
	112	!! 6) Compute increments of tate variables and come back to old values
[921]	113	!!
[2715]	114	!! References : Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
	115	!! Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
	116	!! Rothrock, D. A., 1975: JGR, 80, 4514-4519.
	117	!! Thorndike et al., 1975, JGR, 80, 4501-4513.
	118	!! Bitz et al., JGR, 2001
	119	!! Amundrud and Melling, JGR 2005
	120	!! Babko et al., JGR 2002
[921]	121	!!
[1572]	122	!! This routine is based on CICE code and authors William H. Lipscomb,
	123	!! and Elizabeth C. Hunke, LANL are gratefully acknowledged
[921]	124	!!--------------------------------------------------------------------!
[5080]	125	INTEGER :: ji, jj, jk, jl ! dummy loop index
	126	INTEGER :: niter ! local integer
[3625]	127	INTEGER :: iterate_ridging ! if true, repeat the ridging
[5080]	128	REAL(wp) :: za, zfac ! local scalar
[2715]	129	CHARACTER (len = 15) :: fieldid
[3294]	130	REAL(wp), POINTER, DIMENSION(:,:) :: closing_net ! net rate at which area is removed (1/s)
	131	! (ridging ice area - area of new ridges) / dt
	132	REAL(wp), POINTER, DIMENSION(:,:) :: divu_adv ! divu as implied by transport scheme (1/s)
	133	REAL(wp), POINTER, DIMENSION(:,:) :: opning ! rate of opening due to divergence/shear
	134	REAL(wp), POINTER, DIMENSION(:,:) :: closing_gross ! rate at which area removed, not counting area of new ridges
	135	REAL(wp), POINTER, DIMENSION(:,:) :: msnow_mlt ! mass of snow added to ocean (kg m-2)
	136	REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)
	137	REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories
[4688]	138	!
[5080]	139	INTEGER, PARAMETER :: nitermax = 20
	140	!
[4688]	141	REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
[2715]	142	!!-----------------------------------------------------------------------------
[4161]	143	IF( nn_timing == 1 ) CALL timing_start('limitd_me')
[825]	144
[3294]	145	CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
[825]	146
	147	IF(ln_ctl) THEN
	148	CALL prt_ctl(tab2d_1=ato_i , clinfo1=' lim_itd_me: ato_i : ', tab2d_2=at_i , clinfo2=' at_i : ')
	149	CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ')
	150	ENDIF
	151
[4161]	152	IF( ln_limdyn ) THEN ! Start ridging and rafting !
	153
[4688]	154	! conservation test
	155	IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4161]	156
[921]	157	!-----------------------------------------------------------------------------!
	158	! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
	159	!-----------------------------------------------------------------------------!
[5070]	160	Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE
[3625]	161	!
	162	CALL lim_itd_me_ridgeprep ! prepare ridging
	163	!
[2715]	164	IF( con_i) CALL lim_column_sum( jpl, v_i, vt_i_init ) ! conservation check
[825]	165
[2715]	166	DO jj = 1, jpj ! Initialize arrays.
[825]	167	DO ji = 1, jpi
[2715]	168	msnow_mlt(ji,jj) = 0._wp
	169	esnow_mlt(ji,jj) = 0._wp
[3625]	170	dardg1dt (ji,jj) = 0._wp
	171	dardg2dt (ji,jj) = 0._wp
	172	dvirdgdt (ji,jj) = 0._wp
	173	opening (ji,jj) = 0._wp
[825]	174
[921]	175	!-----------------------------------------------------------------------------!
	176	! 2) Dynamical inputs (closing rate, divu_adv, opning)
	177	!-----------------------------------------------------------------------------!
	178	!
	179	! 2.1 closing_net
	180	!-----------------
	181	! Compute the net rate of closing due to convergence
	182	! and shear, based on Flato and Hibler (1995).
	183	!
	184	! The energy dissipation rate is equal to the net closing rate
	185	! times the ice strength.
	186	!
	187	! NOTE: The NET closing rate is equal to the rate that open water
	188	! area is removed, plus the rate at which ice area is removed by
	189	! ridging, minus the rate at which area is added in new ridges.
	190	! The GROSS closing rate is equal to the first two terms (open
	191	! water closing and thin ice ridging) without the third term
	192	! (thick, newly ridged ice).
[825]	193
[5070]	194	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]	195
[921]	196	! 2.2 divu_adv
	197	!--------------
	198	! Compute divu_adv, the divergence rate given by the transport/
	199	! advection scheme, which may not be equal to divu as computed
	200	! from the velocity field.
	201	!
	202	! If divu_adv < 0, make sure the closing rate is large enough
	203	! to give asum = 1.0 after ridging.
[825]	204
[4161]	205	divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep
[825]	206
[2715]	207	IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
[825]	208
[921]	209	! 2.3 opning
	210	!------------
[3625]	211	! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
[921]	212	opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
[825]	213	END DO
	214	END DO
	215
[921]	216	!-----------------------------------------------------------------------------!
	217	! 3) Ridging iteration
	218	!-----------------------------------------------------------------------------!
[2715]	219	niter = 1 ! iteration counter
[869]	220	iterate_ridging = 1
[825]	221
[869]	222	DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
[825]	223
[921]	224	DO jj = 1, jpj
	225	DO ji = 1, jpi
[825]	226
[921]	227	! 3.2 closing_gross
	228	!-----------------------------------------------------------------------------!
	229	! Based on the ITD of ridging and ridged ice, convert the net
	230	! closing rate to a gross closing rate.
	231	! NOTE: 0 < aksum <= 1
	232	closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj)
[825]	233
[921]	234	! correction to closing rate and opening if closing rate is excessive
	235	!---------------------------------------------------------------------
	236	! Reduce the closing rate if more than 100% of the open water
	237	! would be removed. Reduce the opening rate proportionately.
[5070]	238	IF ( ato_i(ji,jj) > epsi10 .AND. athorn(ji,jj,0) > 0.0 ) THEN
[5080]	239	za = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice
	240	IF ( za > ato_i(ji,jj)) THEN
	241	zfac = ato_i(ji,jj) / za
	242	closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
	243	opning(ji,jj) = opning(ji,jj) * zfac
	244	ENDIF
	245	ENDIF
[825]	246
[5080]	247	END DO
	248	END DO
[825]	249
[921]	250	! correction to closing rate / opening if excessive ice removal
	251	!---------------------------------------------------------------
	252	! Reduce the closing rate if more than 100% of any ice category
	253	! would be removed. Reduce the opening rate proportionately.
[825]	254
[921]	255	DO jl = 1, jpl
	256	DO jj = 1, jpj
	257	DO ji = 1, jpi
[4333]	258	IF ( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp )THEN
[5080]	259	za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
	260	IF ( za > a_i(ji,jj,jl) ) THEN
	261	zfac = a_i(ji,jj,jl) / za
	262	closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
	263	opning (ji,jj) = opning (ji,jj) * zfac
[921]	264	ENDIF
[825]	265	ENDIF
[5080]	266	END DO
	267	END DO
	268	END DO
[825]	269
[921]	270	! 3.3 Redistribute area, volume, and energy.
	271	!-----------------------------------------------------------------------------!
[825]	272
[2715]	273	CALL lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
[825]	274
[921]	275	! 3.4 Compute total area of ice plus open water after ridging.
	276	!-----------------------------------------------------------------------------!
[825]	277
[921]	278	CALL lim_itd_me_asumr
[825]	279
[921]	280	! 3.5 Do we keep on iterating ???
	281	!-----------------------------------------------------------------------------!
	282	! Check whether asum = 1. If not (because the closing and opening
	283	! rates were reduced above), ridge again with new rates.
[825]	284
[921]	285	iterate_ridging = 0
[825]	286
[921]	287	DO jj = 1, jpj
	288	DO ji = 1, jpi
[5070]	289	IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN
[2715]	290	closing_net(ji,jj) = 0._wp
	291	opning (ji,jj) = 0._wp
[921]	292	ELSE
	293	iterate_ridging = 1
[4161]	294	divu_adv (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice
[2715]	295	closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
	296	opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) )
[921]	297	ENDIF
	298	END DO
[825]	299	END DO
	300
[2715]	301	IF( lk_mpp ) CALL mpp_max( iterate_ridging )
[869]	302
[921]	303	! Repeat if necessary.
	304	! NOTE: If strength smoothing is turned on, the ridging must be
	305	! iterated globally because of the boundary update in the
	306	! smoothing.
[825]	307
[921]	308	niter = niter + 1
[825]	309
[2715]	310	IF( iterate_ridging == 1 ) THEN
[3625]	311	IF( niter > nitermax ) THEN
[921]	312	WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
	313	WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
	314	ENDIF
	315	CALL lim_itd_me_ridgeprep
[825]	316	ENDIF
	317
	318	END DO !! on the do while over iter
	319
[921]	320	!-----------------------------------------------------------------------------!
	321	! 4) Ridging diagnostics
	322	!-----------------------------------------------------------------------------!
	323	! Convert ridging rate diagnostics to correct units.
	324	! Update fresh water and heat fluxes due to snow melt.
[825]	325	DO jj = 1, jpj
	326	DO ji = 1, jpi
	327
[3625]	328	dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice
	329	dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice
	330	dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice
	331	opening (ji,jj) = opening (ji,jj) * r1_rdtice
[825]	332
[921]	333	!-----------------------------------------------------------------------------!
	334	! 5) Heat, salt and freshwater fluxes
	335	!-----------------------------------------------------------------------------!
[4688]	336	wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean
[5064]	337	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice ! heat sink for ocean (<0, W.m-2)
[825]	338
	339	END DO
	340	END DO
	341
	342	! Check if there is a ridging error
	343	DO jj = 1, jpj
	344	DO ji = 1, jpi
[4333]	345	IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug
[825]	346	WRITE(numout,*) ' '
	347	WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
	348	WRITE(numout,*) ' limitd_me '
	349	WRITE(numout,*) ' POINT : ', ji, jj
	350	WRITE(numout,*) ' jpl, a_i, athorn '
	351	WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
	352	DO jl = 1, jpl
	353	WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
	354	END DO
[5070]	355	ENDIF
	356	END DO
	357	END DO
[825]	358
	359	! Conservation check
[834]	360	IF ( con_i ) THEN
	361	CALL lim_column_sum (jpl, v_i, vt_i_final)
	362	fieldid = ' v_i : limitd_me '
	363	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
	364	ENDIF
[825]	365
[921]	366	!-----------------------------------------------------------------------------!
[4333]	367	! 6) Updating state variables and trend terms (done in limupdate)
[921]	368	!-----------------------------------------------------------------------------!
[834]	369	CALL lim_var_glo2eqv
[5051]	370	CALL lim_var_zapsmall
[5053]	371	CALL lim_var_agg( 1 )
[825]	372
[4161]	373
[863]	374	IF(ln_ctl) THEN ! Control print
[867]	375	CALL prt_ctl_info(' ')
	376	CALL prt_ctl_info(' - Cell values : ')
	377	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
[5064]	378	CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_me : cell area :')
[863]	379	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :')
	380	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :')
	381	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me : vt_s :')
	382	DO jl = 1, jpl
[867]	383	CALL prt_ctl_info(' ')
[863]	384	CALL prt_ctl_info(' - Category : ', ivar1=jl)
	385	CALL prt_ctl_info(' ~~~~~~~~~~')
	386	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_me : a_i : ')
	387	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_me : ht_i : ')
	388	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_me : ht_s : ')
	389	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_me : v_i : ')
	390	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_me : v_s : ')
	391	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_me : e_s : ')
	392	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_me : t_su : ')
	393	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_me : t_snow : ')
	394	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_me : sm_i : ')
	395	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_me : smv_i : ')
	396	DO jk = 1, nlay_i
[867]	397	CALL prt_ctl_info(' ')
[863]	398	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
	399	CALL prt_ctl_info(' ~~~~~~~')
	400	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : t_i : ')
	401	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : e_i : ')
	402	END DO
	403	END DO
	404	ENDIF
[825]	405
[4688]	406	! conservation test
	407	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4161]	408
	409	ENDIF ! ln_limdyn=.true.
[3625]	410	!
[3294]	411	CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
[2715]	412	!
[4161]	413	IF( nn_timing == 1 ) CALL timing_stop('limitd_me')
[825]	414	END SUBROUTINE lim_itd_me
	415
	416
[2715]	417	SUBROUTINE lim_itd_me_icestrength( kstrngth )
[921]	418	!!----------------------------------------------------------------------
	419	!! * ROUTINE lim_itd_me_icestrength *
	420	!!
[2715]	421	!! ** Purpose : computes ice strength used in dynamics routines of ice thickness
[921]	422	!!
[2715]	423	!! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2)
	424	!! dissipated per unit area removed from the ice pack under compression,
	425	!! and assumed proportional to the change in potential energy caused
	426	!! by ridging. Note that only Hibler's formulation is stable and that
	427	!! ice strength has to be smoothed
	428	!!
[921]	429	!! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
	430	!!----------------------------------------------------------------------
[2715]	431	INTEGER, INTENT(in) :: kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979)
[5080]	432	INTEGER :: ji,jj, jl ! dummy loop indices
	433	INTEGER :: ksmooth ! smoothing the resistance to deformation
	434	INTEGER :: numts_rm ! number of time steps for the P smoothing
	435	REAL(wp) :: zhi, zp, z1_3 ! local scalars
[3294]	436	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
[2715]	437	!!----------------------------------------------------------------------
[825]	438
[3294]	439	CALL wrk_alloc( jpi, jpj, zworka )
[825]	440
[921]	441	!------------------------------------------------------------------------------!
	442	! 1) Initialize
	443	!------------------------------------------------------------------------------!
[2715]	444	strength(:,:) = 0._wp
[825]	445
[921]	446	!------------------------------------------------------------------------------!
	447	! 2) Compute thickness distribution of ridging and ridged ice
	448	!------------------------------------------------------------------------------!
[825]	449	CALL lim_itd_me_ridgeprep
	450
[921]	451	!------------------------------------------------------------------------------!
	452	! 3) Rothrock(1975)'s method
	453	!------------------------------------------------------------------------------!
[2715]	454	IF( kstrngth == 1 ) THEN
	455	z1_3 = 1._wp / 3._wp
[825]	456	DO jl = 1, jpl
	457	DO jj= 1, jpj
	458	DO ji = 1, jpi
[2715]	459	!
[4333]	460	IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN
[5080]	461	zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
[825]	462	!----------------------------
	463	! PE loss from deforming ice
	464	!----------------------------
[5080]	465	strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * zhi * zhi
[825]	466
	467	!--------------------------
	468	! PE gain from rafting ice
	469	!--------------------------
[5080]	470	strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * zhi * zhi
[825]	471
	472	!----------------------------
	473	! PE gain from ridging ice
	474	!----------------------------
[2715]	475	strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl)/krdg(ji,jj,jl) &
[5080]	476	* z1_3 * ( hrmax(ji,jj,jl)2 + hrmin(ji,jj,jl)2 + hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )
	477	!!(a3-b3)/(a-b) = aa+ab+bb
[5070]	478	ENDIF
[2715]	479	!
[5070]	480	END DO
	481	END DO
	482	END DO
[5080]	483
	484	strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:)
	485	! where Cp = (g/2)(rhow-rhoi)(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
[825]	486	ksmooth = 1
	487
[921]	488	!------------------------------------------------------------------------------!
	489	! 4) Hibler (1979)' method
	490	!------------------------------------------------------------------------------!
[825]	491	ELSE ! kstrngth ne 1: Hibler (1979) form
[2715]	492	!
[5067]	493	strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) )
[2715]	494	!
[825]	495	ksmooth = 1
[2715]	496	!
[825]	497	ENDIF ! kstrngth
	498
[921]	499	!
	500	!------------------------------------------------------------------------------!
	501	! 5) Impact of brine volume
	502	!------------------------------------------------------------------------------!
	503	! CAN BE REMOVED
	504	!
[5070]	505	IF( ln_icestr_bvf ) THEN
[825]	506
	507	DO jj = 1, jpj
	508	DO ji = 1, jpi
	509	strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
[5070]	510	END DO
	511	END DO
[825]	512
	513	ENDIF
	514
[921]	515	!
	516	!------------------------------------------------------------------------------!
	517	! 6) Smoothing ice strength
	518	!------------------------------------------------------------------------------!
	519	!
[825]	520	!-------------------
	521	! Spatial smoothing
	522	!-------------------
[2715]	523	IF ( ksmooth == 1 ) THEN
[825]	524
	525	CALL lbc_lnk( strength, 'T', 1. )
	526
[5080]	527	DO jj = 2, jpjm1
	528	DO ji = 2, jpim1
[5070]	529	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present
[825]	530	zworka(ji,jj) = 4.0 * strength(ji,jj) &
[5067]	531	& + strength(ji-1,jj) * tmask(ji-1,jj,1) &
	532	& + strength(ji+1,jj) * tmask(ji+1,jj,1) &
	533	& + strength(ji,jj-1) * tmask(ji,jj-1,1) &
	534	& + strength(ji,jj+1) * tmask(ji,jj+1,1)
[825]	535
[5080]	536	zworka(ji,jj) = zworka(ji,jj) / &
	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
[5080]	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
[5070]	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
[5070]	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
[5070]	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
[5080]	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
[5067]	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
[5051]	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.
	627	CALL lim_itd_me_asumr
	628	! This is in general not equal to one
	629	! because of divergence during transport
	630
	631	! Compute cumulative thickness distribution function
	632	! Compute the cumulative thickness distribution function Gsum,
	633	! where Gsum(n) is the fractional area in categories 0 to n.
	634	! initial value (in h = 0) equals open water area
	635
[2715]	636	Gsum(:,:,-1) = 0._wp
[825]	637
	638	DO jj = 1, jpj
	639	DO ji = 1, jpi
[4333]	640	IF( ato_i(ji,jj) > epsi10 ) THEN ; Gsum(ji,jj,0) = ato_i(ji,jj)
[2715]	641	ELSE ; Gsum(ji,jj,0) = 0._wp
[825]	642	ENDIF
	643	END DO
	644	END DO
	645
	646	! for each value of h, you have to add ice concentration then
	647	DO jl = 1, jpl
	648	DO jj = 1, jpj
	649	DO ji = 1, jpi
[5070]	650	IF( a_i(ji,jj,jl) > epsi10 ) THEN ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)
	651	ELSE ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1)
[825]	652	ENDIF
	653	END DO
	654	END DO
	655	END DO
	656
	657	! Normalize the cumulative distribution to 1
[2715]	658	zworka(:,:) = 1._wp / Gsum(:,:,jpl)
[825]	659	DO jl = 0, jpl
[2715]	660	Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:)
[825]	661	END DO
	662
[921]	663	! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
	664	!--------------------------------------------------------------------------------------------------
	665	! Compute the participation function athorn; this is analogous to
	666	! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
	667	! area lost from category n due to ridging/closing
	668	! athorn(n) = total area lost due to ridging/closing
	669	! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
	670	!
	671	! The expressions for athorn are found by integrating b(h)g(h) between
	672	! the category boundaries.
	673	!-----------------------------------------------------------------
[825]	674
[5067]	675	IF( nn_partfun == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975)
[4869]	676	DO jl = 0, jpl
[921]	677	DO jj = 1, jpj
	678	DO ji = 1, jpi
[5067]	679	IF( Gsum(ji,jj,jl) < rn_gstar) THEN
[921]	680	athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * &
	681	(2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari)
[5067]	682	ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN
	683	athorn(ji,jj,jl) = Gstari * (rn_gstar-Gsum(ji,jj,jl-1)) * &
	684	(2.0 - (Gsum(ji,jj,jl-1)+rn_gstar)*Gstari)
[921]	685	ELSE
	686	athorn(ji,jj,jl) = 0.0
	687	ENDIF
[5080]	688	END DO
	689	END DO
	690	END DO
[825]	691
[2715]	692	ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007)
	693	!
	694	zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array
[825]	695
[921]	696	DO jl = -1, jpl
[2715]	697	Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
[921]	698	END DO !jl
[4869]	699	DO jl = 0, jpl
[2715]	700	athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
	701	END DO
	702	!
[5067]	703	ENDIF ! nn_partfun
[825]	704
[5070]	705	IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions
[2715]	706	!
[825]	707	DO jl = 1, jpl
	708	DO jj = 1, jpj
	709	DO ji = 1, jpi
[5070]	710	IF ( athorn(ji,jj,jl) > 0._wp ) THEN
[2715]	711	!!gm TANH( -X ) = - TANH( X ) so can be computed only 1 time....
[5067]	712	aridge(ji,jj,jl) = ( TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
	713	araft (ji,jj,jl) = ( TANH ( -rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
[2715]	714	IF ( araft(ji,jj,jl) < epsi06 ) araft(ji,jj,jl) = 0._wp
	715	aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 )
[5070]	716	ENDIF
	717	END DO
	718	END DO
	719	END DO
[825]	720
[5070]	721	ELSE
[2715]	722	!
[825]	723	DO jl = 1, jpl
[2715]	724	aridge(:,:,jl) = athorn(:,:,jl)
[825]	725	END DO
[2715]	726	!
[825]	727	ENDIF
	728
[5070]	729	IF( ln_rafting ) THEN
[825]	730
[5070]	731	IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 ) THEN
[921]	732	DO jl = 1, jpl
	733	DO jj = 1, jpj
	734	DO ji = 1, jpi
[5070]	735	IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN
[921]	736	WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
	737	WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
	738	WRITE(numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj)
	739	WRITE(numout,*) ' aridge : ', aridge(ji,jj,1:jpl)
	740	WRITE(numout,*) ' araft : ', araft(ji,jj,1:jpl)
	741	WRITE(numout,*) ' athorn : ', athorn(ji,jj,1:jpl)
	742	ENDIF
	743	END DO
[868]	744	END DO
[825]	745	END DO
[921]	746	ENDIF
[825]	747
	748	ENDIF
	749
[921]	750	!-----------------------------------------------------------------
	751	! 2) Transfer function
	752	!-----------------------------------------------------------------
	753	! Compute max and min ridged ice thickness for each ridging category.
	754	! Assume ridged ice is uniformly distributed between hrmin and hrmax.
	755	!
	756	! This parameterization is a modified version of Hibler (1980).
	757	! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
	758	! and for very thick ridging ice must be >= krdgmin*hi
	759	!
	760	! The minimum ridging thickness, hrmin, is equal to 2*hi
	761	! (i.e., rafting) and for very thick ridging ice is
	762	! constrained by hrmin <= (hrmean + hi)/2.
	763	!
	764	! The maximum ridging thickness, hrmax, is determined by
	765	! hrmean and hrmin.
	766	!
	767	! These modifications have the effect of reducing the ice strength
	768	! (relative to the Hibler formulation) when very thick ice is
	769	! ridging.
	770	!
	771	! aksum = net area removed/ total area removed
	772	! where total area removed = area of ice that ridges
	773	! net area removed = total area removed - area of new ridges
	774	!-----------------------------------------------------------------
[825]	775
	776	! Transfer function
	777	DO jl = 1, jpl !all categories have a specific transfer function
	778	DO jj = 1, jpj
	779	DO ji = 1, jpi
	780
[5070]	781	IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN
[5080]	782	zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
	783	hrmean = MAX(SQRT(rn_hstarzhi), zhikrdgmin)
	784	hrmin(ji,jj,jl) = MIN(2.0zhi, 0.5(hrmean + zhi))
[825]	785	hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)
[5080]	786	hraft(ji,jj,jl) = kraft*zhi
	787	krdg(ji,jj,jl) = hrmean / zhi
[825]	788	ELSE
	789	hraft(ji,jj,jl) = 0.0
	790	hrmin(ji,jj,jl) = 0.0
	791	hrmax(ji,jj,jl) = 0.0
	792	krdg (ji,jj,jl) = 1.0
	793	ENDIF
	794
	795	END DO ! ji
	796	END DO ! jj
	797	END DO ! jl
	798
	799	! Normalization factor : aksum, ensures mass conservation
[2777]	800	aksum(:,:) = athorn(:,:,0)
[825]	801	DO jl = 1, jpl
[2715]	802	aksum(:,:) = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) ) &
	803	& + araft (:,:,jl) * ( 1._wp - 1._wp / kraft )
[825]	804	END DO
[2715]	805	!
[3294]	806	CALL wrk_dealloc( jpi,jpj, zworka )
	807	CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
	808	!
[921]	809	END SUBROUTINE lim_itd_me_ridgeprep
[825]	810
	811
[2715]	812	SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
	813	!!----------------------------------------------------------------------
[921]	814	!! * ROUTINE lim_itd_me_icestrength *
	815	!!
[2715]	816	!! ** Purpose : shift ridging ice among thickness categories of ice thickness
[921]	817	!!
[2715]	818	!! ** Method : Remove area, volume, and energy from each ridging category
	819	!! and add to thicker ice categories.
	820	!!----------------------------------------------------------------------
	821	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear
	822	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area removed, excluding area of new ridges
	823	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: msnow_mlt ! mass of snow added to ocean (kg m-2)
	824	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)
	825	!
	826	CHARACTER (len=80) :: fieldid ! field identifier
	827	LOGICAL, PARAMETER :: l_conservation_check = .true. ! if true, check conservation (useful for debugging)
	828	!
	829	LOGICAL :: neg_ato_i ! flag for ato_i(i,j) < -puny
	830	LOGICAL :: large_afrac ! flag for afrac > 1
	831	LOGICAL :: large_afrft ! flag for afrac > 1
	832	INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices
	833	INTEGER :: ij ! horizontal index, combines i and j loops
	834	INTEGER :: icells ! number of cells with aicen > puny
[5080]	835	REAL(wp) :: hL, hR, farea, ztmelts ! left and right limits of integration
[4688]	836	REAL(wp) :: zsstK ! SST in Kelvin
[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 )
[4333]	889	DO ji = mi0(jiindx), mi1(jiindx)
	890	DO jj = mj0(jjindx), mj1(jjindx)
	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
	912	END DO !jj
	913	END DO !ji
	914
	915	! if negative open water area alert it
[2715]	916	IF( neg_ato_i ) 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)
[4333]	923	ENDIF ! ato_i < -epsi10
[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)
[825]	939	END DO !jl
[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
	970	ENDIF ! test on a_icen_init
	971	END DO ! ji
	972	END DO ! jj
	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)
[5067]	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
[5067]	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
[5067]	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
[5064]	1081	! in J/m2 (same as e_s)
[5067]	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
	1096	END DO ! ij
	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?)
	1114	zsstK = sst_m(ji,jj) + rt0
[5078]	1115	ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) * r1_nlay_i
[825]	1116
[4688]	1117	! heat flux to the ocean
	1118	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux
[825]	1119
[5064]	1120	! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
[921]	1121	erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
[4688]	1122
[825]	1123	END DO ! ij
	1124	END DO !jk
	1125
	1126
[2715]	1127	IF( con_i ) THEN
[825]	1128	DO jk = 1, nlay_i
	1129	DO ij = 1, icells
	1130	ji = indxi(ij)
	1131	jj = indxj(ij)
[2715]	1132	eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk)
[825]	1133	END DO ! ij
	1134	END DO !jk
	1135	ENDIF
	1136
[2715]	1137	IF( large_afrac ) THEN ! there is a bug
[825]	1138	DO ij = 1, icells
	1139	ji = indxi(ij)
	1140	jj = indxj(ij)
[4333]	1141	IF( afrac(ji,jj) > kamax + epsi10 ) THEN
[825]	1142	WRITE(numout,*) ''
	1143	WRITE(numout,*) ' ardg > a_i'
[2715]	1144	WRITE(numout,*) ' ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1)
	1145	ENDIF
	1146	END DO
	1147	ENDIF
	1148	IF( large_afrft ) THEN ! there is a bug
[825]	1149	DO ij = 1, icells
	1150	ji = indxi(ij)
	1151	jj = indxj(ij)
[4333]	1152	IF( afrft(ji,jj) > kamax + epsi10 ) THEN
[825]	1153	WRITE(numout,*) ''
	1154	WRITE(numout,*) ' arft > a_i'
[2715]	1155	WRITE(numout,*) ' arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)
	1156	ENDIF
	1157	END DO
	1158	ENDIF
[825]	1159
[921]	1160	!-------------------------------------------------------------------------------
	1161	! 4) Add area, volume, and energy of new ridge to each category jl2
	1162	!-------------------------------------------------------------------------------
	1163	! jl1 looping 1-jpl
[4869]	1164	DO jl2 = 1, jpl
[921]	1165	! over categories to which ridged ice is transferred
[825]	1166	DO ij = 1, icells
	1167	ji = indxi(ij)
	1168	jj = indxj(ij)
	1169
	1170	! Compute the fraction of ridged ice area and volume going to
	1171	! thickness category jl2.
	1172	! Transfer area, volume, and energy accordingly.
	1173
[3625]	1174	IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. &
	1175	hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN
	1176	hL = 0._wp
	1177	hR = 0._wp
[825]	1178	ELSE
[3625]	1179	hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
	1180	hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2) )
[825]	1181	ENDIF
	1182
	1183	! fraction of ridged ice area and volume going to n2
[3625]	1184	farea = ( hR - hL ) / dhr(ji,jj)
	1185	fvol(ji,jj) = ( hRhR - hLhL ) / dhr2(ji,jj)
[825]	1186
[3625]	1187	a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ardg2 (ji,jj) * farea
	1188	v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj)
[5067]	1189	v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg
	1190	e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg
[3625]	1191	smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + srdg2 (ji,jj) * fvol(ji,jj)
	1192	oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirdg2(ji,jj) * farea
[825]	1193
	1194	END DO ! ij
	1195
	1196	! Transfer ice energy to category jl2 by ridging
	1197	DO jk = 1, nlay_i
	1198	DO ij = 1, icells
	1199	ji = indxi(ij)
	1200	jj = indxj(ij)
[2715]	1201	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj)*erdg2(ji,jj,jk)
	1202	END DO
	1203	END DO
	1204	!
[825]	1205	END DO ! jl2 (new ridges)
	1206
[4869]	1207	DO jl2 = 1, jpl
[825]	1208
	1209	DO ij = 1, icells
	1210	ji = indxi(ij)
	1211	jj = indxj(ij)
	1212	! Compute the fraction of rafted ice area and volume going to
	1213	! thickness category jl2, transfer area, volume, and energy accordingly.
[3625]	1214	!
	1215	IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. &
	1216	hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN
	1217	a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + arft2 (ji,jj)
	1218	v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + virft (ji,jj)
[5067]	1219	v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrft (ji,jj) * rn_fsnowrft
	1220	e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrft (ji,jj) * rn_fsnowrft
[3625]	1221	smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + smrft (ji,jj)
	1222	oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirft2(ji,jj)
[5070]	1223	ENDIF
[3625]	1224	!
[5070]	1225	END DO
[825]	1226
	1227	! Transfer rafted ice energy to category jl2
	1228	DO jk = 1, nlay_i
	1229	DO ij = 1, icells
	1230	ji = indxi(ij)
	1231	jj = indxj(ij)
[3625]	1232	IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. &
	1233	hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN
[2715]	1234	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk)
[825]	1235	ENDIF
[5070]	1236	END DO
	1237	END DO
[825]	1238
[5070]	1239	END DO
[825]	1240
	1241	END DO ! jl1 (deforming categories)
	1242
	1243	! Conservation check
	1244	IF ( con_i ) THEN
	1245	CALL lim_column_sum (jpl, v_i, vice_final)
	1246	fieldid = ' v_i : limitd_me '
	1247	CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid)
	1248
	1249	CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_final )
	1250	fieldid = ' e_i : limitd_me '
	1251	CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid)
[4333]	1252
	1253	DO ji = mi0(jiindx), mi1(jiindx)
	1254	DO jj = mj0(jjindx), mj1(jjindx)
	1255	WRITE(numout,*) ' vice_init : ', vice_init (ji,jj)
	1256	WRITE(numout,*) ' vice_final : ', vice_final(ji,jj)
	1257	WRITE(numout,*) ' eice_init : ', eice_init (ji,jj)
	1258	WRITE(numout,*) ' eice_final : ', eice_final(ji,jj)
	1259	END DO
	1260	END DO
[825]	1261	ENDIF
[2715]	1262	!
[3294]	1263	CALL wrk_dealloc( (jpi+1)*(jpj+1), indxi, indxj )
	1264	CALL wrk_dealloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )
	1265	CALL wrk_dealloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
	1266	CALL wrk_dealloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw )
	1267	CALL wrk_dealloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
[4688]	1268	CALL wrk_dealloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
[4873]	1269	CALL wrk_dealloc( jpi, jpj, nlay_i+1, eirft, erdg1, erdg2, ersw )
	1270	CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, eicen_init )
[3294]	1271	!
[825]	1272	END SUBROUTINE lim_itd_me_ridgeshift
	1273
	1274
[2715]	1275	SUBROUTINE lim_itd_me_asumr
[921]	1276	!!-----------------------------------------------------------------------------
	1277	!! * ROUTINE lim_itd_me_asumr *
	1278	!!
[2715]	1279	!! ** Purpose : finds total fractional area
[921]	1280	!!
[2715]	1281	!! ** Method : Find the total area of ice plus open water in each grid cell.
	1282	!! This is similar to the aggregate_area subroutine except that the
	1283	!! total area can be greater than 1, so the open water area is
	1284	!! included in the sum instead of being computed as a residual.
	1285	!!-----------------------------------------------------------------------------
	1286	INTEGER :: jl ! dummy loop index
	1287	!!-----------------------------------------------------------------------------
	1288	!
	1289	asum(:,:) = ato_i(:,:) ! open water
	1290	DO jl = 1, jpl ! ice categories
	1291	asum(:,:) = asum(:,:) + a_i(:,:,jl)
[825]	1292	END DO
[2715]	1293	!
[825]	1294	END SUBROUTINE lim_itd_me_asumr
	1295
	1296
[2715]	1297	SUBROUTINE lim_itd_me_init
[825]	1298	!!-------------------------------------------------------------------
	1299	!! * ROUTINE lim_itd_me_init *
	1300	!!
	1301	!! ** Purpose : Physical constants and parameters linked
	1302	!! to the mechanical ice redistribution
	1303	!!
	1304	!! ** Method : Read the namiceitdme namelist
	1305	!! and check the parameters values
	1306	!! called at the first timestep (nit000)
	1307	!!
	1308	!! ** input : Namelist namiceitdme
	1309	!!-------------------------------------------------------------------
[4298]	1310	INTEGER :: ios ! Local integer output status for namelist read
[5070]	1311	NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft, &
	1312	& rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, &
[5067]	1313	& nn_partfun
[825]	1314	!!-------------------------------------------------------------------
[2715]	1315	!
[4298]	1316	REWIND( numnam_ice_ref ) ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
	1317	READ ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
	1318	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
	1319
	1320	REWIND( numnam_ice_cfg ) ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
	1321	READ ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
	1322	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
[4624]	1323	IF(lwm) WRITE ( numoni, namiceitdme )
[2715]	1324	!
	1325	IF (lwp) THEN ! control print
[825]	1326	WRITE(numout,*)
	1327	WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
	1328	WRITE(numout,*)' ~~~~~~~~~~~~~~~'
[5070]	1329	WRITE(numout,*)' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs
	1330	WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg
	1331	WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft
	1332	WRITE(numout,*)' Fraction of total ice coverage contributing to ridging rn_gstar = ', rn_gstar
	1333	WRITE(numout,)' Equivalent to G for an exponential part function rn_astar = ', rn_astar
	1334	WRITE(numout,*)' Quantity playing a role in max ridged ice thickness rn_hstar = ', rn_hstar
	1335	WRITE(numout,*)' Rafting of ice sheets or not ln_rafting = ', ln_rafting
	1336	WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft
	1337	WRITE(numout,*)' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft
	1338	WRITE(numout,*)' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg
	1339	WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun
[825]	1340	ENDIF
[2715]	1341	!
[825]	1342	END SUBROUTINE lim_itd_me_init
	1343
	1344	#else
[2715]	1345	!!----------------------------------------------------------------------
	1346	!! Default option Empty module NO LIM-3 sea-ice model
	1347	!!----------------------------------------------------------------------
[825]	1348	CONTAINS
	1349	SUBROUTINE lim_itd_me ! Empty routines
	1350	END SUBROUTINE lim_itd_me
	1351	SUBROUTINE lim_itd_me_icestrength
	1352	END SUBROUTINE lim_itd_me_icestrength
	1353	SUBROUTINE lim_itd_me_sort
	1354	END SUBROUTINE lim_itd_me_sort
	1355	SUBROUTINE lim_itd_me_init
	1356	END SUBROUTINE lim_itd_me_init
	1357	#endif
[2715]	1358	!!======================================================================
[825]	1359	END MODULE limitd_me

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

Download in other formats: