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 @ 5076

Last change on this file since 5076 was 5070, checked in by clem, 9 years ago
LIM3 cosmetic changes
Property svn:keywords set to `Id`
File size: 63.2 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	!!--------------------------------------------------------------------!
[2715]	125	INTEGER :: ji, jj, jk, jl ! dummy loop index
	126	INTEGER :: niter, nitermax = 20 ! local integer
[3625]	127	LOGICAL :: asum_error ! flag for asum .ne. 1
	128	INTEGER :: iterate_ridging ! if true, repeat the ridging
	129	REAL(wp) :: w1, tmpfac ! local scalar
[2715]	130	CHARACTER (len = 15) :: fieldid
[3294]	131	REAL(wp), POINTER, DIMENSION(:,:) :: closing_net ! net rate at which area is removed (1/s)
	132	! (ridging ice area - area of new ridges) / dt
	133	REAL(wp), POINTER, DIMENSION(:,:) :: divu_adv ! divu as implied by transport scheme (1/s)
	134	REAL(wp), POINTER, DIMENSION(:,:) :: opning ! rate of opening due to divergence/shear
	135	REAL(wp), POINTER, DIMENSION(:,:) :: closing_gross ! rate at which area removed, not counting area of new ridges
	136	REAL(wp), POINTER, DIMENSION(:,:) :: msnow_mlt ! mass of snow added to ocean (kg m-2)
	137	REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)
	138	REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories
[4688]	139	!
	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	!-----------------------------------------------------------------------------!
[5070]	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
[5070]	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.
[5070]	237	IF ( ato_i(ji,jj) > epsi10 .AND. athorn(ji,jj,0) > 0.0 ) THEN
[921]	238	w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice
[5070]	239	IF ( w1 > ato_i(ji,jj)) THEN
[921]	240	tmpfac = ato_i(ji,jj) / w1
	241	closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac
	242	opning(ji,jj) = opning(ji,jj) * tmpfac
	243	ENDIF !w1
	244	ENDIF !at0i and athorn
[825]	245
[921]	246	END DO ! ji
	247	END DO ! jj
[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
[921]	258	w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
[3625]	259	IF ( w1 > a_i(ji,jj,jl) ) THEN
[921]	260	tmpfac = a_i(ji,jj,jl) / w1
	261	closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac
[2715]	262	opning (ji,jj) = opning (ji,jj) * tmpfac
[921]	263	ENDIF
[825]	264	ENDIF
[921]	265	END DO !ji
	266	END DO ! jj
	267	END DO !jl
[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	!-----------------------------------------------------------------------------!
[825]	276
[921]	277	CALL lim_itd_me_asumr
[825]	278
[921]	279	! 3.5 Do we keep on iterating ???
	280	!-----------------------------------------------------------------------------!
	281	! Check whether asum = 1. If not (because the closing and opening
	282	! rates were reduced above), ridge again with new rates.
[825]	283
[921]	284	iterate_ridging = 0
[825]	285
[921]	286	DO jj = 1, jpj
	287	DO ji = 1, jpi
[5070]	288	IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN
[2715]	289	closing_net(ji,jj) = 0._wp
	290	opning (ji,jj) = 0._wp
[921]	291	ELSE
	292	iterate_ridging = 1
[4161]	293	divu_adv (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice
[2715]	294	closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
	295	opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) )
[921]	296	ENDIF
	297	END DO
[825]	298	END DO
	299
[2715]	300	IF( lk_mpp ) CALL mpp_max( iterate_ridging )
[869]	301
[921]	302	! Repeat if necessary.
	303	! NOTE: If strength smoothing is turned on, the ridging must be
	304	! iterated globally because of the boundary update in the
	305	! smoothing.
[825]	306
[921]	307	niter = niter + 1
[825]	308
[2715]	309	IF( iterate_ridging == 1 ) THEN
[3625]	310	IF( niter > nitermax ) THEN
[921]	311	WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
	312	WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
	313	ENDIF
	314	CALL lim_itd_me_ridgeprep
[825]	315	ENDIF
	316
	317	END DO !! on the do while over iter
	318
[921]	319	!-----------------------------------------------------------------------------!
	320	! 4) Ridging diagnostics
	321	!-----------------------------------------------------------------------------!
	322	! Convert ridging rate diagnostics to correct units.
	323	! Update fresh water and heat fluxes due to snow melt.
[825]	324
	325	asum_error = .false.
	326
	327	DO jj = 1, jpj
	328	DO ji = 1, jpi
	329
[4333]	330	IF(ABS(asum(ji,jj) - kamax) > epsi10 ) asum_error = .true.
[825]	331
[3625]	332	dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice
	333	dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice
	334	dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice
	335	opening (ji,jj) = opening (ji,jj) * r1_rdtice
[825]	336
[921]	337	!-----------------------------------------------------------------------------!
	338	! 5) Heat, salt and freshwater fluxes
	339	!-----------------------------------------------------------------------------!
[4688]	340	wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean
[5064]	341	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice ! heat sink for ocean (<0, W.m-2)
[825]	342
	343	END DO
	344	END DO
	345
	346	! Check if there is a ridging error
	347	DO jj = 1, jpj
	348	DO ji = 1, jpi
[4333]	349	IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug
[825]	350	WRITE(numout,*) ' '
	351	WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
	352	WRITE(numout,*) ' limitd_me '
	353	WRITE(numout,*) ' POINT : ', ji, jj
	354	WRITE(numout,*) ' jpl, a_i, athorn '
	355	WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
	356	DO jl = 1, jpl
	357	WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
	358	END DO
[5070]	359	ENDIF
	360	END DO
	361	END DO
[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
[921]	370	!-----------------------------------------------------------------------------!
[4333]	371	! 6) Updating state variables and trend terms (done in limupdate)
[921]	372	!-----------------------------------------------------------------------------!
[834]	373	CALL lim_var_glo2eqv
[5051]	374	CALL lim_var_zapsmall
[5053]	375	CALL lim_var_agg( 1 )
[825]	376
[4161]	377
[863]	378	IF(ln_ctl) THEN ! Control print
[867]	379	CALL prt_ctl_info(' ')
	380	CALL prt_ctl_info(' - Cell values : ')
	381	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
[5064]	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)
[921]	436
[2715]	437	INTEGER :: ji,jj, jl ! dummy loop indices
	438	INTEGER :: ksmooth ! smoothing the resistance to deformation
	439	INTEGER :: numts_rm ! number of time steps for the P smoothing
	440	REAL(wp) :: hi, zw1, zp, zdummy, zzc, z1_3 ! local scalars
[3294]	441	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
[2715]	442	!!----------------------------------------------------------------------
[825]	443
[3294]	444	CALL wrk_alloc( jpi, jpj, zworka )
[825]	445
[921]	446	!------------------------------------------------------------------------------!
	447	! 1) Initialize
	448	!------------------------------------------------------------------------------!
[2715]	449	strength(:,:) = 0._wp
[825]	450
[921]	451	!------------------------------------------------------------------------------!
	452	! 2) Compute thickness distribution of ridging and ridged ice
	453	!------------------------------------------------------------------------------!
[825]	454	CALL lim_itd_me_ridgeprep
	455
[921]	456	!------------------------------------------------------------------------------!
	457	! 3) Rothrock(1975)'s method
	458	!------------------------------------------------------------------------------!
[2715]	459	IF( kstrngth == 1 ) THEN
	460	z1_3 = 1._wp / 3._wp
[825]	461	DO jl = 1, jpl
	462	DO jj= 1, jpj
	463	DO ji = 1, jpi
[2715]	464	!
[4333]	465	IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN
[825]	466	hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
	467	!----------------------------
	468	! PE loss from deforming ice
	469	!----------------------------
[2715]	470	strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * hi * hi
[825]	471
	472	!--------------------------
	473	! PE gain from rafting ice
	474	!--------------------------
[2715]	475	strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * hi * hi
[825]	476
	477	!----------------------------
	478	! PE gain from ridging ice
	479	!----------------------------
[2715]	480	strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl)/krdg(ji,jj,jl) &
	481	* z1_3 * (hrmax(ji,jj,jl)3 - hrmin(ji,jj,jl)3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) )
	482	!!gm Optimization: (a3-b3)/(a-b) = aa+ab+bb ==> less costly operations even if a*3 is replaced by aa*a...
[5070]	483	ENDIF
[2715]	484	!
[5070]	485	END DO
	486	END DO
	487	END DO
[825]	488
[5067]	489	zzc = rn_pe_rdg * Cp ! where Cp = (g/2)(rhow-rhoi)(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
[2715]	490	strength(:,:) = zzc * strength(:,:) / aksum(:,:)
[921]	491
[825]	492	ksmooth = 1
	493
[921]	494	!------------------------------------------------------------------------------!
	495	! 4) Hibler (1979)' method
	496	!------------------------------------------------------------------------------!
[825]	497	ELSE ! kstrngth ne 1: Hibler (1979) form
[2715]	498	!
[5067]	499	strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) )
[2715]	500	!
[825]	501	ksmooth = 1
[2715]	502	!
[825]	503	ENDIF ! kstrngth
	504
[921]	505	!
	506	!------------------------------------------------------------------------------!
	507	! 5) Impact of brine volume
	508	!------------------------------------------------------------------------------!
	509	! CAN BE REMOVED
	510	!
[5070]	511	IF( ln_icestr_bvf ) THEN
[825]	512
	513	DO jj = 1, jpj
	514	DO ji = 1, jpi
[5070]	515	IF ( bv_i(ji,jj) > 0.0 ) THEN
[825]	516	zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 )
	517	ELSE
	518	zdummy = 0.0
	519	ENDIF
	520	strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
[5070]	521	END DO
	522	END DO
[825]	523
	524	ENDIF
	525
[921]	526	!
	527	!------------------------------------------------------------------------------!
	528	! 6) Smoothing ice strength
	529	!------------------------------------------------------------------------------!
	530	!
[825]	531	!-------------------
	532	! Spatial smoothing
	533	!-------------------
[2715]	534	IF ( ksmooth == 1 ) THEN
[825]	535
	536	CALL lbc_lnk( strength, 'T', 1. )
	537
[869]	538	DO jj = 2, jpj - 1
	539	DO ji = 2, jpi - 1
[5070]	540	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present
[825]	541	zworka(ji,jj) = 4.0 * strength(ji,jj) &
[5067]	542	& + strength(ji-1,jj) * tmask(ji-1,jj,1) &
	543	& + strength(ji+1,jj) * tmask(ji+1,jj,1) &
	544	& + strength(ji,jj-1) * tmask(ji,jj-1,1) &
	545	& + strength(ji,jj+1) * tmask(ji,jj+1,1)
[825]	546
[5067]	547	zw1 = 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1)
[825]	548	zworka(ji,jj) = zworka(ji,jj) / zw1
	549	ELSE
[3625]	550	zworka(ji,jj) = 0._wp
[825]	551	ENDIF
	552	END DO
	553	END DO
	554
[869]	555	DO jj = 2, jpj - 1
	556	DO ji = 2, jpi - 1
[825]	557	strength(ji,jj) = zworka(ji,jj)
	558	END DO
	559	END DO
[869]	560	CALL lbc_lnk( strength, 'T', 1. )
[825]	561
[5070]	562	ENDIF
[825]	563
	564	!--------------------
	565	! Temporal smoothing
	566	!--------------------
[2715]	567	IF ( numit == nit000 + nn_fsbc - 1 ) THEN
[825]	568	strp1(:,:) = 0.0
	569	strp2(:,:) = 0.0
	570	ENDIF
	571
[2715]	572	IF ( ksmooth == 2 ) THEN
[921]	573
	574
[825]	575	CALL lbc_lnk( strength, 'T', 1. )
[921]	576
[825]	577	DO jj = 1, jpj - 1
	578	DO ji = 1, jpi - 1
[5070]	579	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present
[825]	580	numts_rm = 1 ! number of time steps for the running mean
[5070]	581	IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
	582	IF ( strp2(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
[2715]	583	zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
[825]	584	strp2(ji,jj) = strp1(ji,jj)
	585	strp1(ji,jj) = strength(ji,jj)
	586	strength(ji,jj) = zp
	587
	588	ENDIF
	589	END DO
	590	END DO
	591
	592	ENDIF ! ksmooth
[921]	593
[2715]	594	CALL lbc_lnk( strength, 'T', 1. ) ! Boundary conditions
[825]	595
[3294]	596	CALL wrk_dealloc( jpi, jpj, zworka )
[2715]	597	!
[921]	598	END SUBROUTINE lim_itd_me_icestrength
[825]	599
	600
[2715]	601	SUBROUTINE lim_itd_me_ridgeprep
[921]	602	!!---------------------------------------------------------------------!
	603	!! * ROUTINE lim_itd_me_ridgeprep *
	604	!!
[2715]	605	!! ** Purpose : preparation for ridging and strength calculations
[921]	606	!!
[2715]	607	!! ** Method : Compute the thickness distribution of the ice and open water
	608	!! participating in ridging and of the resulting ridges.
[921]	609	!!---------------------------------------------------------------------!
[2715]	610	INTEGER :: ji,jj, jl ! dummy loop indices
	611	REAL(wp) :: Gstari, astari, hi, hrmean, zdummy ! local scalar
[3294]	612	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
	613	REAL(wp), POINTER, DIMENSION(:,:,:) :: Gsum ! Gsum(n) = sum of areas in categories 0 to n
	614	!------------------------------------------------------------------------------!
[825]	615
[3294]	616	CALL wrk_alloc( jpi,jpj, zworka )
	617	CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
[825]	618
[5067]	619	Gstari = 1.0/rn_gstar
	620	astari = 1.0/rn_astar
[825]	621	aksum(:,:) = 0.0
	622	athorn(:,:,:) = 0.0
	623	aridge(:,:,:) = 0.0
	624	araft (:,:,:) = 0.0
	625	hrmin(:,:,:) = 0.0
	626	hrmax(:,:,:) = 0.0
	627	hraft(:,:,:) = 0.0
	628	krdg (:,:,:) = 1.0
	629
[921]	630	! ! Zero out categories with very small areas
[5051]	631	CALL lim_var_zapsmall
[825]	632
[921]	633	!------------------------------------------------------------------------------!
	634	! 1) Participation function
	635	!------------------------------------------------------------------------------!
[825]	636
	637	! Compute total area of ice plus open water.
	638	CALL lim_itd_me_asumr
	639	! This is in general not equal to one
	640	! because of divergence during transport
	641
	642	! Compute cumulative thickness distribution function
	643	! Compute the cumulative thickness distribution function Gsum,
	644	! where Gsum(n) is the fractional area in categories 0 to n.
	645	! initial value (in h = 0) equals open water area
	646
[2715]	647	Gsum(:,:,-1) = 0._wp
[825]	648
	649	DO jj = 1, jpj
	650	DO ji = 1, jpi
[4333]	651	IF( ato_i(ji,jj) > epsi10 ) THEN ; Gsum(ji,jj,0) = ato_i(ji,jj)
[2715]	652	ELSE ; Gsum(ji,jj,0) = 0._wp
[825]	653	ENDIF
	654	END DO
	655	END DO
	656
	657	! for each value of h, you have to add ice concentration then
	658	DO jl = 1, jpl
	659	DO jj = 1, jpj
	660	DO ji = 1, jpi
[5070]	661	IF( a_i(ji,jj,jl) > epsi10 ) THEN ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)
	662	ELSE ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1)
[825]	663	ENDIF
	664	END DO
	665	END DO
	666	END DO
	667
	668	! Normalize the cumulative distribution to 1
[2715]	669	zworka(:,:) = 1._wp / Gsum(:,:,jpl)
[825]	670	DO jl = 0, jpl
[2715]	671	Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:)
[825]	672	END DO
	673
[921]	674	! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
	675	!--------------------------------------------------------------------------------------------------
	676	! Compute the participation function athorn; this is analogous to
	677	! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
	678	! area lost from category n due to ridging/closing
	679	! athorn(n) = total area lost due to ridging/closing
	680	! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
	681	!
	682	! The expressions for athorn are found by integrating b(h)g(h) between
	683	! the category boundaries.
	684	!-----------------------------------------------------------------
[825]	685
[5067]	686	IF( nn_partfun == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975)
[4869]	687	DO jl = 0, jpl
[921]	688	DO jj = 1, jpj
	689	DO ji = 1, jpi
[5067]	690	IF( Gsum(ji,jj,jl) < rn_gstar) THEN
[921]	691	athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * &
	692	(2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari)
[5067]	693	ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN
	694	athorn(ji,jj,jl) = Gstari * (rn_gstar-Gsum(ji,jj,jl-1)) * &
	695	(2.0 - (Gsum(ji,jj,jl-1)+rn_gstar)*Gstari)
[921]	696	ELSE
	697	athorn(ji,jj,jl) = 0.0
	698	ENDIF
	699	END DO ! ji
	700	END DO ! jj
	701	END DO ! jl
[825]	702
[2715]	703	ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007)
	704	!
	705	zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array
[825]	706
[921]	707	DO jl = -1, jpl
[2715]	708	Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
[921]	709	END DO !jl
[4869]	710	DO jl = 0, jpl
[2715]	711	athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
	712	END DO
	713	!
[5067]	714	ENDIF ! nn_partfun
[825]	715
[5070]	716	IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions
[2715]	717	!
[825]	718	DO jl = 1, jpl
	719	DO jj = 1, jpj
	720	DO ji = 1, jpi
[5070]	721	IF ( athorn(ji,jj,jl) > 0._wp ) THEN
[2715]	722	!!gm TANH( -X ) = - TANH( X ) so can be computed only 1 time....
[5067]	723	aridge(ji,jj,jl) = ( TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
	724	araft (ji,jj,jl) = ( TANH ( -rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
[2715]	725	IF ( araft(ji,jj,jl) < epsi06 ) araft(ji,jj,jl) = 0._wp
	726	aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 )
[5070]	727	ENDIF
	728	END DO
	729	END DO
	730	END DO
[825]	731
[5070]	732	ELSE
[2715]	733	!
[825]	734	DO jl = 1, jpl
[2715]	735	aridge(:,:,jl) = athorn(:,:,jl)
[825]	736	END DO
[2715]	737	!
[825]	738	ENDIF
	739
[5070]	740	IF( ln_rafting ) THEN
[825]	741
[5070]	742	IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 ) THEN
[921]	743	DO jl = 1, jpl
	744	DO jj = 1, jpj
	745	DO ji = 1, jpi
[5070]	746	IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN
[921]	747	WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
	748	WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
	749	WRITE(numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj)
	750	WRITE(numout,*) ' aridge : ', aridge(ji,jj,1:jpl)
	751	WRITE(numout,*) ' araft : ', araft(ji,jj,1:jpl)
	752	WRITE(numout,*) ' athorn : ', athorn(ji,jj,1:jpl)
	753	ENDIF
	754	END DO
[868]	755	END DO
[825]	756	END DO
[921]	757	ENDIF
[825]	758
	759	ENDIF
	760
[921]	761	!-----------------------------------------------------------------
	762	! 2) Transfer function
	763	!-----------------------------------------------------------------
	764	! Compute max and min ridged ice thickness for each ridging category.
	765	! Assume ridged ice is uniformly distributed between hrmin and hrmax.
	766	!
	767	! This parameterization is a modified version of Hibler (1980).
	768	! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
	769	! and for very thick ridging ice must be >= krdgmin*hi
	770	!
	771	! The minimum ridging thickness, hrmin, is equal to 2*hi
	772	! (i.e., rafting) and for very thick ridging ice is
	773	! constrained by hrmin <= (hrmean + hi)/2.
	774	!
	775	! The maximum ridging thickness, hrmax, is determined by
	776	! hrmean and hrmin.
	777	!
	778	! These modifications have the effect of reducing the ice strength
	779	! (relative to the Hibler formulation) when very thick ice is
	780	! ridging.
	781	!
	782	! aksum = net area removed/ total area removed
	783	! where total area removed = area of ice that ridges
	784	! net area removed = total area removed - area of new ridges
	785	!-----------------------------------------------------------------
[825]	786
	787	! Transfer function
	788	DO jl = 1, jpl !all categories have a specific transfer function
	789	DO jj = 1, jpj
	790	DO ji = 1, jpi
	791
[5070]	792	IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN
[825]	793	hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
[5067]	794	hrmean = MAX(SQRT(rn_hstarhi), hikrdgmin)
[825]	795	hrmin(ji,jj,jl) = MIN(2.0hi, 0.5(hrmean + hi))
	796	hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)
	797	hraft(ji,jj,jl) = kraft*hi
	798	krdg(ji,jj,jl) = hrmean / hi
	799	ELSE
	800	hraft(ji,jj,jl) = 0.0
	801	hrmin(ji,jj,jl) = 0.0
	802	hrmax(ji,jj,jl) = 0.0
	803	krdg (ji,jj,jl) = 1.0
	804	ENDIF
	805
	806	END DO ! ji
	807	END DO ! jj
	808	END DO ! jl
	809
	810	! Normalization factor : aksum, ensures mass conservation
[2777]	811	aksum(:,:) = athorn(:,:,0)
[825]	812	DO jl = 1, jpl
[2715]	813	aksum(:,:) = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) ) &
	814	& + araft (:,:,jl) * ( 1._wp - 1._wp / kraft )
[825]	815	END DO
[2715]	816	!
[3294]	817	CALL wrk_dealloc( jpi,jpj, zworka )
	818	CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
	819	!
[921]	820	END SUBROUTINE lim_itd_me_ridgeprep
[825]	821
	822
[2715]	823	SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
	824	!!----------------------------------------------------------------------
[921]	825	!! * ROUTINE lim_itd_me_icestrength *
	826	!!
[2715]	827	!! ** Purpose : shift ridging ice among thickness categories of ice thickness
[921]	828	!!
[2715]	829	!! ** Method : Remove area, volume, and energy from each ridging category
	830	!! and add to thicker ice categories.
	831	!!----------------------------------------------------------------------
	832	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear
	833	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area removed, excluding area of new ridges
	834	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: msnow_mlt ! mass of snow added to ocean (kg m-2)
	835	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)
	836	!
	837	CHARACTER (len=80) :: fieldid ! field identifier
	838	LOGICAL, PARAMETER :: l_conservation_check = .true. ! if true, check conservation (useful for debugging)
	839	!
	840	LOGICAL :: neg_ato_i ! flag for ato_i(i,j) < -puny
	841	LOGICAL :: large_afrac ! flag for afrac > 1
	842	LOGICAL :: large_afrft ! flag for afrac > 1
	843	INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices
	844	INTEGER :: ij ! horizontal index, combines i and j loops
	845	INTEGER :: icells ! number of cells with aicen > puny
	846	REAL(wp) :: hL, hR, farea, zdummy, zdummy0, ztmelts ! left and right limits of integration
[4688]	847	REAL(wp) :: zsstK ! SST in Kelvin
[825]	848
[3294]	849	INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices
[825]	850
[3294]	851	REAL(wp), POINTER, DIMENSION(:,:) :: vice_init, vice_final ! ice volume summed over categories
	852	REAL(wp), POINTER, DIMENSION(:,:) :: eice_init, eice_final ! ice energy summed over layers
[825]	853
[3294]	854	REAL(wp), POINTER, DIMENSION(:,:,:) :: aicen_init, vicen_init ! ice area & volume before ridging
[4688]	855	REAL(wp), POINTER, DIMENSION(:,:,:) :: vsnwn_init, esnwn_init ! snow volume & energy before ridging
[3294]	856	REAL(wp), POINTER, DIMENSION(:,:,:) :: smv_i_init, oa_i_init ! ice salinity & age before ridging
[825]	857
[3294]	858	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: eicen_init ! ice energy before ridging
[825]	859
[3294]	860	REAL(wp), POINTER, DIMENSION(:,:) :: afrac , fvol ! fraction of category area ridged & new ridge volume going to n2
	861	REAL(wp), POINTER, DIMENSION(:,:) :: ardg1 , ardg2 ! area of ice ridged & new ridges
	862	REAL(wp), POINTER, DIMENSION(:,:) :: vsrdg , esrdg ! snow volume & energy of ridging ice
	863	REAL(wp), POINTER, DIMENSION(:,:) :: oirdg1, oirdg2 ! areal age content of ridged & rifging ice
	864	REAL(wp), POINTER, DIMENSION(:,:) :: dhr , dhr2 ! hrmax - hrmin & hrmax^2 - hrmin^2
[825]	865
[3294]	866	REAL(wp), POINTER, DIMENSION(:,:) :: vrdg1 ! volume of ice ridged
	867	REAL(wp), POINTER, DIMENSION(:,:) :: vrdg2 ! volume of new ridges
	868	REAL(wp), POINTER, DIMENSION(:,:) :: vsw ! volume of seawater trapped into ridges
	869	REAL(wp), POINTER, DIMENSION(:,:) :: srdg1 ! sal*volume of ice ridged
	870	REAL(wp), POINTER, DIMENSION(:,:) :: srdg2 ! sal*volume of new ridges
	871	REAL(wp), POINTER, DIMENSION(:,:) :: smsw ! sal*volume of water trapped into ridges
[825]	872
[3294]	873	REAL(wp), POINTER, DIMENSION(:,:) :: afrft ! fraction of category area rafted
	874	REAL(wp), POINTER, DIMENSION(:,:) :: arft1 , arft2 ! area of ice rafted and new rafted zone
	875	REAL(wp), POINTER, DIMENSION(:,:) :: virft , vsrft ! ice & snow volume of rafting ice
	876	REAL(wp), POINTER, DIMENSION(:,:) :: esrft , smrft ! snow energy & salinity of rafting ice
	877	REAL(wp), POINTER, DIMENSION(:,:) :: oirft1, oirft2 ! areal age content of rafted ice & rafting ice
[825]	878
[3294]	879	REAL(wp), POINTER, DIMENSION(:,:,:) :: eirft ! ice energy of rafting ice
	880	REAL(wp), POINTER, DIMENSION(:,:,:) :: erdg1 ! enth*volume of ice ridged
	881	REAL(wp), POINTER, DIMENSION(:,:,:) :: erdg2 ! enth*volume of new ridges
	882	REAL(wp), POINTER, DIMENSION(:,:,:) :: ersw ! enth of water trapped into ridges
	883	!!----------------------------------------------------------------------
[825]	884
[3294]	885	CALL wrk_alloc( (jpi+1)*(jpj+1), indxi, indxj )
	886	CALL wrk_alloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )
	887	CALL wrk_alloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
	888	CALL wrk_alloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw )
	889	CALL wrk_alloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
[4688]	890	CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
[4873]	891	CALL wrk_alloc( jpi, jpj, nlay_i+1, eirft, erdg1, erdg2, ersw )
	892	CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, eicen_init )
[3294]	893
[825]	894	! Conservation check
[2715]	895	eice_init(:,:) = 0._wp
[825]	896
[2715]	897	IF( con_i ) THEN
[4333]	898	CALL lim_column_sum (jpl, v_i, vice_init )
[825]	899	CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_init )
[4333]	900	DO ji = mi0(jiindx), mi1(jiindx)
	901	DO jj = mj0(jjindx), mj1(jjindx)
	902	WRITE(numout,*) ' vice_init : ', vice_init(ji,jj)
	903	WRITE(numout,*) ' eice_init : ', eice_init(ji,jj)
	904	END DO
	905	END DO
[825]	906	ENDIF
	907
[921]	908	!-------------------------------------------------------------------------------
	909	! 1) Compute change in open water area due to closing and opening.
	910	!-------------------------------------------------------------------------------
[825]	911
	912	neg_ato_i = .false.
	913
	914	DO jj = 1, jpj
	915	DO ji = 1, jpi
[2715]	916	ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice &
	917	& + opning(ji,jj) * rdt_ice
[4333]	918	IF( ato_i(ji,jj) < -epsi10 ) THEN
[2715]	919	neg_ato_i = .TRUE.
	920	ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error
	921	ato_i(ji,jj) = 0._wp
[825]	922	ENDIF
	923	END DO !jj
	924	END DO !ji
	925
	926	! if negative open water area alert it
[2715]	927	IF( neg_ato_i ) THEN ! there is a bug
[825]	928	DO jj = 1, jpj
	929	DO ji = 1, jpi
[4333]	930	IF( ato_i(ji,jj) < -epsi10 ) THEN
[825]	931	WRITE(numout,*) ''
	932	WRITE(numout,*) 'Ridging error: ato_i < 0'
	933	WRITE(numout,*) 'ato_i : ', ato_i(ji,jj)
[4333]	934	ENDIF ! ato_i < -epsi10
[2715]	935	END DO
	936	END DO
	937	ENDIF
[825]	938
[921]	939	!-----------------------------------------------------------------
	940	! 2) Save initial state variables
	941	!-----------------------------------------------------------------
[825]	942
	943	DO jl = 1, jpl
[2715]	944	aicen_init(:,:,jl) = a_i(:,:,jl)
	945	vicen_init(:,:,jl) = v_i(:,:,jl)
[4688]	946	vsnwn_init(:,:,jl) = v_s(:,:,jl)
[2715]	947	!
	948	smv_i_init(:,:,jl) = smv_i(:,:,jl)
	949	oa_i_init (:,:,jl) = oa_i (:,:,jl)
[825]	950	END DO !jl
[868]	951
[4688]	952	esnwn_init(:,:,:) = e_s(:,:,1,:)
[921]	953
[825]	954	DO jl = 1, jpl
	955	DO jk = 1, nlay_i
[2715]	956	eicen_init(:,:,jk,jl) = e_i(:,:,jk,jl)
	957	END DO
	958	END DO
[825]	959
[921]	960	!
	961	!-----------------------------------------------------------------
	962	! 3) Pump everything from ice which is being ridged / rafted
	963	!-----------------------------------------------------------------
	964	! Compute the area, volume, and energy of ice ridging in each
	965	! category, along with the area of the resulting ridge.
[825]	966
	967	DO jl1 = 1, jpl !jl1 describes the ridging category
	968
[921]	969	!------------------------------------------------
	970	! 3.1) Identify grid cells with nonzero ridging
	971	!------------------------------------------------
[825]	972
	973	icells = 0
	974	DO jj = 1, jpj
	975	DO ji = 1, jpi
[4333]	976	IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp &
	977	& .AND. closing_gross(ji,jj) > 0._wp ) THEN
[825]	978	icells = icells + 1
	979	indxi(icells) = ji
	980	indxj(icells) = jj
	981	ENDIF ! test on a_icen_init
	982	END DO ! ji
	983	END DO ! jj
	984
	985	large_afrac = .false.
	986	large_afrft = .false.
	987
	988	DO ij = 1, icells
	989	ji = indxi(ij)
	990	jj = indxj(ij)
	991
[921]	992	!--------------------------------------------------------------------
	993	! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
	994	!--------------------------------------------------------------------
[825]	995
	996	ardg1(ji,jj) = aridge(ji,jj,jl1)closing_gross(ji,jj)rdt_ice
	997	arft1(ji,jj) = araft (ji,jj,jl1)closing_gross(ji,jj)rdt_ice
	998	ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1)
	999	arft2(ji,jj) = arft1(ji,jj) / kraft
	1000
	1001	oirdg1(ji,jj)= aridge(ji,jj,jl1)closing_gross(ji,jj)rdt_ice
	1002	oirft1(ji,jj)= araft (ji,jj,jl1)closing_gross(ji,jj)rdt_ice
	1003	oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1)
	1004	oirft2(ji,jj)= oirft1(ji,jj) / kraft
	1005
[921]	1006	!---------------------------------------------------------------
	1007	! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
	1008	!---------------------------------------------------------------
[825]	1009
	1010	afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging
	1011	afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting
	1012
[4333]	1013	IF (afrac(ji,jj) > kamax + epsi10) THEN !riging
[825]	1014	large_afrac = .true.
[4161]	1015	ELSEIF (afrac(ji,jj) > kamax) THEN ! roundoff error
	1016	afrac(ji,jj) = kamax
[825]	1017	ENDIF
[4333]	1018	IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting
[825]	1019	large_afrft = .true.
[4161]	1020	ELSEIF (afrft(ji,jj) > kamax) THEN ! roundoff error
	1021	afrft(ji,jj) = kamax
[825]	1022	ENDIF
	1023
[921]	1024	!--------------------------------------------------------------------------
	1025	! 3.4) Subtract area, volume, and energy from ridging
	1026	! / rafting category n1.
	1027	!--------------------------------------------------------------------------
[4788]	1028	vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj)
[5067]	1029	vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + rn_por_rdg )
	1030	vsw (ji,jj) = vrdg1(ji,jj) * rn_por_rdg
[825]	1031
[4688]	1032	vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj)
	1033	esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj)
[4788]	1034	srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
[4688]	1035	srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless
[825]	1036
	1037	! rafting volumes, heat contents ...
	1038	virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)
[4688]	1039	vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj)
	1040	esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj)
[825]	1041	smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)
	1042
	1043	! substract everything
	1044	a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - ardg1(ji,jj) - arft1(ji,jj)
	1045	v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - vrdg1(ji,jj) - virft(ji,jj)
	1046	v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - vsrdg(ji,jj) - vsrft(ji,jj)
	1047	e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj) - esrft(ji,jj)
	1048	oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1) - oirdg1(ji,jj) - oirft1(ji,jj)
	1049	smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj) - smrft(ji,jj)
	1050
[921]	1051	!-----------------------------------------------------------------
	1052	! 3.5) Compute properties of new ridges
	1053	!-----------------------------------------------------------------
[825]	1054	!-------------
	1055	! Salinity
	1056	!-------------
[4688]	1057	smsw(ji,jj) = vsw(ji,jj) * sss_m(ji,jj) ! salt content of seawater frozen in voids !! MV HC2014
	1058	srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj) ! salt content of new ridge
[1572]	1059
[5067]	1060	!srdg2(ji,jj) = MIN( rn_simax * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity
[1572]	1061
[4688]	1062	sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice
	1063	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]	1064
[921]	1065	!------------------------------------
	1066	! 3.6 Increment ridging diagnostics
	1067	!------------------------------------
[825]	1068
[921]	1069	! jl1 looping 1-jpl
	1070	! ij looping 1-icells
[825]	1071
[2715]	1072	dardg1dt (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)
	1073	dardg2dt (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)
[3625]	1074	opening (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice
[825]	1075
[2715]	1076	IF( con_i ) vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj)
[825]	1077
[921]	1078	!------------------------------------------
	1079	! 3.7 Put the snow somewhere in the ocean
	1080	!------------------------------------------
	1081	! Place part of the snow lost by ridging into the ocean.
	1082	! Note that esnow_mlt < 0; the ocean must cool to melt snow.
	1083	! If the ocean temp = Tf already, new ice must grow.
	1084	! During the next time step, thermo_rates will determine whether
	1085	! the ocean cools or new ice grows.
	1086	! jl1 looping 1-jpl
	1087	! ij looping 1-icells
	1088
[5067]	1089	msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosnvsrdg(ji,jj)(1.0-rn_fsnowrdg) & ! rafting included
	1090	& + rhosnvsrft(ji,jj)(1.0-rn_fsnowrft)
[825]	1091
[5064]	1092	! in J/m2 (same as e_s)
[5067]	1093	esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-rn_fsnowrdg) & !rafting included
	1094	& - esrft(ji,jj)*(1.0-rn_fsnowrft)
[825]	1095
[921]	1096	!-----------------------------------------------------------------
	1097	! 3.8 Compute quantities used to apportion ice among categories
	1098	! in the n2 loop below
	1099	!-----------------------------------------------------------------
[825]	1100
[921]	1101	! jl1 looping 1-jpl
	1102	! ij looping 1-icells
[825]	1103
[3625]	1104	dhr (ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1)
[2715]	1105	dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1)
[825]	1106
	1107	END DO ! ij
	1108
[921]	1109	!--------------------------------------------------------------------
	1110	! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
	1111	! compute ridged ice enthalpy
	1112	!--------------------------------------------------------------------
[825]	1113	DO jk = 1, nlay_i
	1114	DO ij = 1, icells
[921]	1115	ji = indxi(ij)
	1116	jj = indxj(ij)
	1117	! heat content of ridged ice
[4788]	1118	erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj)
[921]	1119	eirft(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)
[2715]	1120	e_i (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk)
[4688]	1121
	1122
	1123	! enthalpy of the trapped seawater (J/m2, >0)
	1124	! clem: if sst>0, then ersw <0 (is that possible?)
	1125	zsstK = sst_m(ji,jj) + rt0
	1126	ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) / REAL( nlay_i )
[825]	1127
[4688]	1128	! heat flux to the ocean
	1129	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux
[825]	1130
[5064]	1131	! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
[921]	1132	erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
[4688]	1133
[825]	1134	END DO ! ij
	1135	END DO !jk
	1136
	1137
[2715]	1138	IF( con_i ) THEN
[825]	1139	DO jk = 1, nlay_i
	1140	DO ij = 1, icells
	1141	ji = indxi(ij)
	1142	jj = indxj(ij)
[2715]	1143	eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk)
[825]	1144	END DO ! ij
	1145	END DO !jk
	1146	ENDIF
	1147
[2715]	1148	IF( large_afrac ) THEN ! there is a bug
[825]	1149	DO ij = 1, icells
	1150	ji = indxi(ij)
	1151	jj = indxj(ij)
[4333]	1152	IF( afrac(ji,jj) > kamax + epsi10 ) THEN
[825]	1153	WRITE(numout,*) ''
	1154	WRITE(numout,*) ' ardg > a_i'
[2715]	1155	WRITE(numout,*) ' ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1)
	1156	ENDIF
	1157	END DO
	1158	ENDIF
	1159	IF( large_afrft ) THEN ! there is a bug
[825]	1160	DO ij = 1, icells
	1161	ji = indxi(ij)
	1162	jj = indxj(ij)
[4333]	1163	IF( afrft(ji,jj) > kamax + epsi10 ) THEN
[825]	1164	WRITE(numout,*) ''
	1165	WRITE(numout,*) ' arft > a_i'
[2715]	1166	WRITE(numout,*) ' arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)
	1167	ENDIF
	1168	END DO
	1169	ENDIF
[825]	1170
[921]	1171	!-------------------------------------------------------------------------------
	1172	! 4) Add area, volume, and energy of new ridge to each category jl2
	1173	!-------------------------------------------------------------------------------
	1174	! jl1 looping 1-jpl
[4869]	1175	DO jl2 = 1, jpl
[921]	1176	! over categories to which ridged ice is transferred
[825]	1177	DO ij = 1, icells
	1178	ji = indxi(ij)
	1179	jj = indxj(ij)
	1180
	1181	! Compute the fraction of ridged ice area and volume going to
	1182	! thickness category jl2.
	1183	! Transfer area, volume, and energy accordingly.
	1184
[3625]	1185	IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. &
	1186	hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN
	1187	hL = 0._wp
	1188	hR = 0._wp
[825]	1189	ELSE
[3625]	1190	hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
	1191	hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2) )
[825]	1192	ENDIF
	1193
	1194	! fraction of ridged ice area and volume going to n2
[3625]	1195	farea = ( hR - hL ) / dhr(ji,jj)
	1196	fvol(ji,jj) = ( hRhR - hLhL ) / dhr2(ji,jj)
[825]	1197
[3625]	1198	a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ardg2 (ji,jj) * farea
	1199	v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj)
[5067]	1200	v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg
	1201	e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg
[3625]	1202	smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + srdg2 (ji,jj) * fvol(ji,jj)
	1203	oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirdg2(ji,jj) * farea
[825]	1204
	1205	END DO ! ij
	1206
	1207	! Transfer ice energy to category jl2 by ridging
	1208	DO jk = 1, nlay_i
	1209	DO ij = 1, icells
	1210	ji = indxi(ij)
	1211	jj = indxj(ij)
[2715]	1212	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj)*erdg2(ji,jj,jk)
	1213	END DO
	1214	END DO
	1215	!
[825]	1216	END DO ! jl2 (new ridges)
	1217
[4869]	1218	DO jl2 = 1, jpl
[825]	1219
	1220	DO ij = 1, icells
	1221	ji = indxi(ij)
	1222	jj = indxj(ij)
	1223	! Compute the fraction of rafted ice area and volume going to
	1224	! thickness category jl2, transfer area, volume, and energy accordingly.
[3625]	1225	!
	1226	IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. &
	1227	hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN
	1228	a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + arft2 (ji,jj)
	1229	v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + virft (ji,jj)
[5067]	1230	v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrft (ji,jj) * rn_fsnowrft
	1231	e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrft (ji,jj) * rn_fsnowrft
[3625]	1232	smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + smrft (ji,jj)
	1233	oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirft2(ji,jj)
[5070]	1234	ENDIF
[3625]	1235	!
[5070]	1236	END DO
[825]	1237
	1238	! Transfer rafted ice energy to category jl2
	1239	DO jk = 1, nlay_i
	1240	DO ij = 1, icells
	1241	ji = indxi(ij)
	1242	jj = indxj(ij)
[3625]	1243	IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. &
	1244	hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN
[2715]	1245	e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk)
[825]	1246	ENDIF
[5070]	1247	END DO
	1248	END DO
[825]	1249
[5070]	1250	END DO
[825]	1251
	1252	END DO ! jl1 (deforming categories)
	1253
	1254	! Conservation check
	1255	IF ( con_i ) THEN
	1256	CALL lim_column_sum (jpl, v_i, vice_final)
	1257	fieldid = ' v_i : limitd_me '
	1258	CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid)
	1259
	1260	CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_final )
	1261	fieldid = ' e_i : limitd_me '
	1262	CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid)
[4333]	1263
	1264	DO ji = mi0(jiindx), mi1(jiindx)
	1265	DO jj = mj0(jjindx), mj1(jjindx)
	1266	WRITE(numout,*) ' vice_init : ', vice_init (ji,jj)
	1267	WRITE(numout,*) ' vice_final : ', vice_final(ji,jj)
	1268	WRITE(numout,*) ' eice_init : ', eice_init (ji,jj)
	1269	WRITE(numout,*) ' eice_final : ', eice_final(ji,jj)
	1270	END DO
	1271	END DO
[825]	1272	ENDIF
[2715]	1273	!
[3294]	1274	CALL wrk_dealloc( (jpi+1)*(jpj+1), indxi, indxj )
	1275	CALL wrk_dealloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )
	1276	CALL wrk_dealloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
	1277	CALL wrk_dealloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw )
	1278	CALL wrk_dealloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
[4688]	1279	CALL wrk_dealloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
[4873]	1280	CALL wrk_dealloc( jpi, jpj, nlay_i+1, eirft, erdg1, erdg2, ersw )
	1281	CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, eicen_init )
[3294]	1282	!
[825]	1283	END SUBROUTINE lim_itd_me_ridgeshift
	1284
	1285
[2715]	1286	SUBROUTINE lim_itd_me_asumr
[921]	1287	!!-----------------------------------------------------------------------------
	1288	!! * ROUTINE lim_itd_me_asumr *
	1289	!!
[2715]	1290	!! ** Purpose : finds total fractional area
[921]	1291	!!
[2715]	1292	!! ** Method : Find the total area of ice plus open water in each grid cell.
	1293	!! This is similar to the aggregate_area subroutine except that the
	1294	!! total area can be greater than 1, so the open water area is
	1295	!! included in the sum instead of being computed as a residual.
	1296	!!-----------------------------------------------------------------------------
	1297	INTEGER :: jl ! dummy loop index
	1298	!!-----------------------------------------------------------------------------
	1299	!
	1300	asum(:,:) = ato_i(:,:) ! open water
	1301	DO jl = 1, jpl ! ice categories
	1302	asum(:,:) = asum(:,:) + a_i(:,:,jl)
[825]	1303	END DO
[2715]	1304	!
[825]	1305	END SUBROUTINE lim_itd_me_asumr
	1306
	1307
[2715]	1308	SUBROUTINE lim_itd_me_init
[825]	1309	!!-------------------------------------------------------------------
	1310	!! * ROUTINE lim_itd_me_init *
	1311	!!
	1312	!! ** Purpose : Physical constants and parameters linked
	1313	!! to the mechanical ice redistribution
	1314	!!
	1315	!! ** Method : Read the namiceitdme namelist
	1316	!! and check the parameters values
	1317	!! called at the first timestep (nit000)
	1318	!!
	1319	!! ** input : Namelist namiceitdme
	1320	!!-------------------------------------------------------------------
[4298]	1321	INTEGER :: ios ! Local integer output status for namelist read
[5070]	1322	NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft, &
	1323	& rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, &
[5067]	1324	& nn_partfun
[825]	1325	!!-------------------------------------------------------------------
[2715]	1326	!
[4298]	1327	REWIND( numnam_ice_ref ) ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
	1328	READ ( numnam_ice_ref, namiceitdme, IOSTAT = ios, ERR = 901)
	1329	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in reference namelist', lwp )
	1330
	1331	REWIND( numnam_ice_cfg ) ! Namelist namiceitdme in configuration namelist : Ice mechanical ice redistribution
	1332	READ ( numnam_ice_cfg, namiceitdme, IOSTAT = ios, ERR = 902 )
	1333	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitdme in configuration namelist', lwp )
[4624]	1334	IF(lwm) WRITE ( numoni, namiceitdme )
[2715]	1335	!
	1336	IF (lwp) THEN ! control print
[825]	1337	WRITE(numout,*)
	1338	WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
	1339	WRITE(numout,*)' ~~~~~~~~~~~~~~~'
[5070]	1340	WRITE(numout,*)' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs
	1341	WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg
	1342	WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft
	1343	WRITE(numout,*)' Fraction of total ice coverage contributing to ridging rn_gstar = ', rn_gstar
	1344	WRITE(numout,)' Equivalent to G for an exponential part function rn_astar = ', rn_astar
	1345	WRITE(numout,*)' Quantity playing a role in max ridged ice thickness rn_hstar = ', rn_hstar
	1346	WRITE(numout,*)' Rafting of ice sheets or not ln_rafting = ', ln_rafting
	1347	WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft
	1348	WRITE(numout,*)' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft
	1349	WRITE(numout,*)' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg
	1350	WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun
[825]	1351	ENDIF
[2715]	1352	!
[825]	1353	END SUBROUTINE lim_itd_me_init
	1354
	1355	#else
[2715]	1356	!!----------------------------------------------------------------------
	1357	!! Default option Empty module NO LIM-3 sea-ice model
	1358	!!----------------------------------------------------------------------
[825]	1359	CONTAINS
	1360	SUBROUTINE lim_itd_me ! Empty routines
	1361	END SUBROUTINE lim_itd_me
	1362	SUBROUTINE lim_itd_me_icestrength
	1363	END SUBROUTINE lim_itd_me_icestrength
	1364	SUBROUTINE lim_itd_me_sort
	1365	END SUBROUTINE lim_itd_me_sort
	1366	SUBROUTINE lim_itd_me_init
	1367	END SUBROUTINE lim_itd_me_init
	1368	#endif
[2715]	1369	!!======================================================================
[825]	1370	END MODULE limitd_me

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

Download in other formats: