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/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 5600

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

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

Download in other formats: