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

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

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

Context Navigation

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

Last change on this file since 5187 was 5181, checked in by vancop, 9 years ago
Remove out-of-bound exceptions
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
[4688]	7	!! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_dyn
[2715]	8	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
[825]	9	!!----------------------------------------------------------------------
[1572]	10	#if defined key_lim3
	11	!!----------------------------------------------------------------------
[3625]	12	!! 'key_lim3' LIM-3 sea-ice model
[1572]	13	!!----------------------------------------------------------------------
[4161]	14	USE par_oce ! ocean parameters
	15	USE dom_oce ! ocean domain
	16	USE phycst ! physical constants (ocean directory)
	17	USE sbc_oce ! surface boundary condition: ocean fields
	18	USE thd_ice ! LIM thermodynamics
	19	USE ice ! LIM variables
	20	USE dom_ice ! LIM domain
	21	USE limvar ! LIM
	22	USE lbclnk ! lateral boundary condition - MPP exchanges
	23	USE lib_mpp ! MPP library
	24	USE wrk_nemo ! work arrays
	25	USE prtctl ! Print control
[5123]	26
[5134]	27	USE in_out_manager ! I/O manager
[4161]	28	USE iom ! I/O manager
[4688]	29	USE lib_fortran ! glob_sum
[4161]	30	USE limdiahsb
[4688]	31	USE timing ! Timing
	32	USE limcons ! conservation tests
[921]	33
[825]	34	IMPLICIT NONE
	35	PRIVATE
	36
[2715]	37	PUBLIC lim_itd_me ! called by ice_stp
	38	PUBLIC lim_itd_me_icestrength
	39	PUBLIC lim_itd_me_init
[5123]	40	PUBLIC lim_itd_me_alloc ! called by sbc_lim_init
[825]	41
[921]	42	!-----------------------------------------------------------------------
	43	! Variables shared among ridging subroutines
	44	!-----------------------------------------------------------------------
[4333]	45	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: asum ! sum of total ice and open water area
	46	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged
	47	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/
	48	! ! closing associated w/ category n
	49	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness
	50	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness
	51	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice
	52	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: krdg ! mean ridge thickness/thickness of ridging ice
	53	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: aridge ! participating ice ridging
	54	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: araft ! participating ice rafting
[825]	55
[2715]	56	REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier
	57	REAL(wp), PARAMETER :: kraft = 2.0_wp ! rafting multipliyer
[4333]	58	REAL(wp), PARAMETER :: kamax = 1.0_wp ! max of ice area authorized (clem: scheme is not stable if kamax <= 0.99)
[825]	59
[2715]	60	REAL(wp) :: Cp !
[921]	61	!
	62	!-----------------------------------------------------------------------
	63	! Ridging diagnostic arrays for history files
	64	!-----------------------------------------------------------------------
[2715]	65	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dardg1dt ! rate of fractional area loss by ridging ice (1/s)
	66	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dardg2dt ! rate of fractional area gain by new ridges (1/s)
	67	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dvirdgdt ! rate of ice volume ridged (m/s)
	68	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: opening ! rate of opening due to divergence/shear (1/s)
[3625]	69	!
[825]	70	!!----------------------------------------------------------------------
[2528]	71	!! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
[1156]	72	!! $Id$
[2715]	73	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	74	!!----------------------------------------------------------------------
[2715]	75	CONTAINS
[825]	76
[2715]	77	INTEGER FUNCTION lim_itd_me_alloc()
	78	!!---------------------------------------------------------------------!
	79	!! * ROUTINE lim_itd_me_alloc *
	80	!!---------------------------------------------------------------------!
	81	ALLOCATE( &
	82	!* Variables shared among ridging subroutines
	83	& asum (jpi,jpj) , athorn(jpi,jpj,0:jpl) , &
	84	& aksum(jpi,jpj) , &
	85	!
	86	& hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , &
	87	& hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , &
	88	!
	89	!* Ridging diagnostic arrays for history files
	90	& dardg1dt(jpi,jpj) , dardg2dt(jpi,jpj) , &
	91	& dvirdgdt(jpi,jpj) , opening(jpi,jpj) , STAT=lim_itd_me_alloc )
	92	!
	93	IF( lim_itd_me_alloc /= 0 ) CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' )
	94	!
	95	END FUNCTION lim_itd_me_alloc
[1156]	96
[825]	97
[2715]	98	SUBROUTINE lim_itd_me
[921]	99	!!---------------------------------------------------------------------!
	100	!! * ROUTINE lim_itd_me *
	101	!!
[2715]	102	!! ** Purpose : computes the mechanical redistribution of ice thickness
[921]	103	!!
[2715]	104	!! ** Method : Steps :
	105	!! 1) Thickness categories boundaries, ice / o.w. concentrations
	106	!! Ridge preparation
	107	!! 2) Dynamical inputs (closing rate, divu_adv, opning)
	108	!! 3) Ridging iteration
	109	!! 4) Ridging diagnostics
	110	!! 5) Heat, salt and freshwater fluxes
	111	!! 6) Compute increments of tate variables and come back to old values
[921]	112	!!
[2715]	113	!! References : Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
	114	!! Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
	115	!! Rothrock, D. A., 1975: JGR, 80, 4514-4519.
	116	!! Thorndike et al., 1975, JGR, 80, 4501-4513.
	117	!! Bitz et al., JGR, 2001
	118	!! Amundrud and Melling, JGR 2005
	119	!! Babko et al., JGR 2002
[921]	120	!!
[1572]	121	!! This routine is based on CICE code and authors William H. Lipscomb,
	122	!! and Elizabeth C. Hunke, LANL are gratefully acknowledged
[921]	123	!!--------------------------------------------------------------------!
[5123]	124	INTEGER :: ji, jj, jk, jl ! dummy loop index
	125	INTEGER :: niter ! local integer
[3625]	126	INTEGER :: iterate_ridging ! if true, repeat the ridging
[5123]	127	REAL(wp) :: za, zfac ! local scalar
[2715]	128	CHARACTER (len = 15) :: fieldid
[3294]	129	REAL(wp), POINTER, DIMENSION(:,:) :: closing_net ! net rate at which area is removed (1/s)
	130	! (ridging ice area - area of new ridges) / dt
	131	REAL(wp), POINTER, DIMENSION(:,:) :: divu_adv ! divu as implied by transport scheme (1/s)
	132	REAL(wp), POINTER, DIMENSION(:,:) :: opning ! rate of opening due to divergence/shear
	133	REAL(wp), POINTER, DIMENSION(:,:) :: closing_gross ! rate at which area removed, not counting area of new ridges
	134	REAL(wp), POINTER, DIMENSION(:,:) :: msnow_mlt ! mass of snow added to ocean (kg m-2)
	135	REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)
	136	REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories
[4688]	137	!
[5123]	138	INTEGER, PARAMETER :: nitermax = 20
	139	!
[4688]	140	REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
[2715]	141	!!-----------------------------------------------------------------------------
[4161]	142	IF( nn_timing == 1 ) CALL timing_start('limitd_me')
[825]	143
[3294]	144	CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
[825]	145
	146	IF(ln_ctl) THEN
	147	CALL prt_ctl(tab2d_1=ato_i , clinfo1=' lim_itd_me: ato_i : ', tab2d_2=at_i , clinfo2=' at_i : ')
	148	CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ')
	149	ENDIF
	150
[4161]	151	IF( ln_limdyn ) THEN ! Start ridging and rafting !
	152
[4688]	153	! conservation test
	154	IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4161]	155
[5167]	156	CALL lim_var_glo2eqv ! equivalent variables, requested for rafting
[921]	157	!-----------------------------------------------------------------------------!
	158	! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
	159	!-----------------------------------------------------------------------------!
[5134]	160	Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE
[3625]	161	!
	162	CALL lim_itd_me_ridgeprep ! prepare ridging
	163	!
[2715]	164	IF( con_i) CALL lim_column_sum( jpl, v_i, vt_i_init ) ! conservation check
[825]	165
[2715]	166	DO jj = 1, jpj ! Initialize arrays.
[825]	167	DO ji = 1, jpi
[2715]	168	msnow_mlt(ji,jj) = 0._wp
	169	esnow_mlt(ji,jj) = 0._wp
[3625]	170	dardg1dt (ji,jj) = 0._wp
	171	dardg2dt (ji,jj) = 0._wp
	172	dvirdgdt (ji,jj) = 0._wp
	173	opening (ji,jj) = 0._wp
[825]	174
[921]	175	!-----------------------------------------------------------------------------!
	176	! 2) Dynamical inputs (closing rate, divu_adv, opning)
	177	!-----------------------------------------------------------------------------!
	178	!
	179	! 2.1 closing_net
	180	!-----------------
	181	! Compute the net rate of closing due to convergence
	182	! and shear, based on Flato and Hibler (1995).
	183	!
	184	! The energy dissipation rate is equal to the net closing rate
	185	! times the ice strength.
	186	!
	187	! NOTE: The NET closing rate is equal to the rate that open water
	188	! area is removed, plus the rate at which ice area is removed by
	189	! ridging, minus the rate at which area is added in new ridges.
	190	! The GROSS closing rate is equal to the first two terms (open
	191	! water closing and thin ice ridging) without the third term
	192	! (thick, newly ridged ice).
[825]	193
[5123]	194	closing_net(ji,jj) = rn_cs * 0.5 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
[825]	195
[921]	196	! 2.2 divu_adv
	197	!--------------
	198	! Compute divu_adv, the divergence rate given by the transport/
	199	! advection scheme, which may not be equal to divu as computed
	200	! from the velocity field.
	201	!
	202	! If divu_adv < 0, make sure the closing rate is large enough
	203	! to give asum = 1.0 after ridging.
[825]	204
[4161]	205	divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep
[825]	206
[2715]	207	IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
[825]	208
[921]	209	! 2.3 opning
	210	!------------
[3625]	211	! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
[921]	212	opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
[825]	213	END DO
	214	END DO
	215
[921]	216	!-----------------------------------------------------------------------------!
	217	! 3) Ridging iteration
	218	!-----------------------------------------------------------------------------!
[2715]	219	niter = 1 ! iteration counter
[869]	220	iterate_ridging = 1
[825]	221
[869]	222	DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
[825]	223
[921]	224	DO jj = 1, jpj
	225	DO ji = 1, jpi
[825]	226
[921]	227	! 3.2 closing_gross
	228	!-----------------------------------------------------------------------------!
	229	! Based on the ITD of ridging and ridged ice, convert the net
	230	! closing rate to a gross closing rate.
	231	! NOTE: 0 < aksum <= 1
	232	closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj)
[825]	233
[921]	234	! correction to closing rate and opening if closing rate is excessive
	235	!---------------------------------------------------------------------
	236	! Reduce the closing rate if more than 100% of the open water
	237	! would be removed. Reduce the opening rate proportionately.
[5167]	238	za = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice
	239	IF( za > epsi20 ) THEN
	240	zfac = MIN( 1._wp, ato_i(ji,jj) / za )
	241	closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
	242	opning (ji,jj) = opning (ji,jj) * zfac
[5123]	243	ENDIF
[825]	244
[5123]	245	END DO
	246	END DO
[825]	247
[921]	248	! correction to closing rate / opening if excessive ice removal
	249	!---------------------------------------------------------------
	250	! Reduce the closing rate if more than 100% of any ice category
	251	! would be removed. Reduce the opening rate proportionately.
	252	DO jl = 1, jpl
	253	DO jj = 1, jpj
	254	DO ji = 1, jpi
[5167]	255	za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
	256	IF( za > epsi20 ) THEN
	257	zfac = MIN( 1._wp, a_i(ji,jj,jl) / za )
	258	closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
	259	opning (ji,jj) = opning (ji,jj) * zfac
[825]	260	ENDIF
[5123]	261	END DO
	262	END DO
	263	END DO
[825]	264
[921]	265	! 3.3 Redistribute area, volume, and energy.
	266	!-----------------------------------------------------------------------------!
[825]	267
[2715]	268	CALL lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
[825]	269
[921]	270	! 3.4 Compute total area of ice plus open water after ridging.
	271	!-----------------------------------------------------------------------------!
[5134]	272	! This is in general not equal to one because of divergence during transport
	273	asum(:,:) = ato_i(:,:)
	274	DO jl = 1, jpl
	275	asum(:,:) = asum(:,:) + a_i(:,:,jl)
	276	END DO
[825]	277
[921]	278	! 3.5 Do we keep on iterating ???
	279	!-----------------------------------------------------------------------------!
	280	! Check whether asum = 1. If not (because the closing and opening
	281	! rates were reduced above), ridge again with new rates.
[825]	282
[921]	283	iterate_ridging = 0
[825]	284
[921]	285	DO jj = 1, jpj
	286	DO ji = 1, jpi
[5123]	287	IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN
[2715]	288	closing_net(ji,jj) = 0._wp
	289	opning (ji,jj) = 0._wp
[921]	290	ELSE
	291	iterate_ridging = 1
[4161]	292	divu_adv (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice
[2715]	293	closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
	294	opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) )
[921]	295	ENDIF
	296	END DO
[825]	297	END DO
	298
[2715]	299	IF( lk_mpp ) CALL mpp_max( iterate_ridging )
[869]	300
[921]	301	! Repeat if necessary.
	302	! NOTE: If strength smoothing is turned on, the ridging must be
	303	! iterated globally because of the boundary update in the
	304	! smoothing.
[825]	305
[921]	306	niter = niter + 1
[825]	307
[2715]	308	IF( iterate_ridging == 1 ) THEN
[3625]	309	IF( niter > nitermax ) THEN
[921]	310	WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
	311	WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
	312	ENDIF
	313	CALL lim_itd_me_ridgeprep
[825]	314	ENDIF
	315
	316	END DO !! on the do while over iter
	317
[921]	318	!-----------------------------------------------------------------------------!
	319	! 4) Ridging diagnostics
	320	!-----------------------------------------------------------------------------!
	321	! Convert ridging rate diagnostics to correct units.
	322	! Update fresh water and heat fluxes due to snow melt.
[825]	323	DO jj = 1, jpj
	324	DO ji = 1, jpi
	325
[3625]	326	dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice
	327	dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice
	328	dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice
	329	opening (ji,jj) = opening (ji,jj) * r1_rdtice
[825]	330
[921]	331	!-----------------------------------------------------------------------------!
	332	! 5) Heat, salt and freshwater fluxes
	333	!-----------------------------------------------------------------------------!
[4688]	334	wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean
[5123]	335	hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice ! heat sink for ocean (<0, W.m-2)
[825]	336
	337	END DO
	338	END DO
	339
	340	! Check if there is a ridging error
[5134]	341	IF( lwp ) THEN
	342	DO jj = 1, jpj
	343	DO ji = 1, jpi
	344	IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug
	345	WRITE(numout,*) ' '
	346	WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
	347	WRITE(numout,*) ' limitd_me '
	348	WRITE(numout,*) ' POINT : ', ji, jj
	349	WRITE(numout,*) ' jpl, a_i, athorn '
	350	WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
	351	DO jl = 1, jpl
	352	WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
	353	END DO
	354	ENDIF
	355	END DO
[5123]	356	END DO
[5134]	357	END IF
[825]	358
	359	! Conservation check
[834]	360	IF ( con_i ) THEN
	361	CALL lim_column_sum (jpl, v_i, vt_i_final)
	362	fieldid = ' v_i : limitd_me '
	363	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
	364	ENDIF
[825]	365
[5134]	366	! updates
[5123]	367	CALL lim_var_agg( 1 )
[825]	368
[5134]	369	!-----------------------------------------------------------------------------!
	370	! control prints
	371	!-----------------------------------------------------------------------------!
	372	IF(ln_ctl) THEN
[5167]	373	CALL lim_var_glo2eqv
	374
[867]	375	CALL prt_ctl_info(' ')
	376	CALL prt_ctl_info(' - Cell values : ')
	377	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
[5123]	378	CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_me : cell area :')
[863]	379	CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :')
	380	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :')
	381	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me : vt_s :')
	382	DO jl = 1, jpl
[867]	383	CALL prt_ctl_info(' ')
[863]	384	CALL prt_ctl_info(' - Category : ', ivar1=jl)
	385	CALL prt_ctl_info(' ~~~~~~~~~~')
	386	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_me : a_i : ')
	387	CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_me : ht_i : ')
	388	CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_me : ht_s : ')
	389	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_me : v_i : ')
	390	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_me : v_s : ')
	391	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_me : e_s : ')
	392	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_me : t_su : ')
	393	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_me : t_snow : ')
	394	CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_me : sm_i : ')
	395	CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_me : smv_i : ')
	396	DO jk = 1, nlay_i
[867]	397	CALL prt_ctl_info(' ')
[863]	398	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
	399	CALL prt_ctl_info(' ~~~~~~~')
	400	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : t_i : ')
	401	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : e_i : ')
	402	END DO
	403	END DO
	404	ENDIF
[825]	405
[4688]	406	! conservation test
	407	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4161]	408
	409	ENDIF ! ln_limdyn=.true.
[3625]	410	!
[3294]	411	CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
[2715]	412	!
[4161]	413	IF( nn_timing == 1 ) CALL timing_stop('limitd_me')
[825]	414	END SUBROUTINE lim_itd_me
	415
	416
[2715]	417	SUBROUTINE lim_itd_me_icestrength( kstrngth )
[921]	418	!!----------------------------------------------------------------------
	419	!! * ROUTINE lim_itd_me_icestrength *
	420	!!
[2715]	421	!! ** Purpose : computes ice strength used in dynamics routines of ice thickness
[921]	422	!!
[2715]	423	!! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2)
	424	!! dissipated per unit area removed from the ice pack under compression,
	425	!! and assumed proportional to the change in potential energy caused
	426	!! by ridging. Note that only Hibler's formulation is stable and that
	427	!! ice strength has to be smoothed
	428	!!
[921]	429	!! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using)
	430	!!----------------------------------------------------------------------
[2715]	431	INTEGER, INTENT(in) :: kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979)
[5123]	432	INTEGER :: ji,jj, jl ! dummy loop indices
	433	INTEGER :: ksmooth ! smoothing the resistance to deformation
	434	INTEGER :: numts_rm ! number of time steps for the P smoothing
	435	REAL(wp) :: zhi, zp, z1_3 ! local scalars
[3294]	436	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
[2715]	437	!!----------------------------------------------------------------------
[825]	438
[3294]	439	CALL wrk_alloc( jpi, jpj, zworka )
[825]	440
[921]	441	!------------------------------------------------------------------------------!
	442	! 1) Initialize
	443	!------------------------------------------------------------------------------!
[2715]	444	strength(:,:) = 0._wp
[825]	445
[921]	446	!------------------------------------------------------------------------------!
	447	! 2) Compute thickness distribution of ridging and ridged ice
	448	!------------------------------------------------------------------------------!
[825]	449	CALL lim_itd_me_ridgeprep
	450
[921]	451	!------------------------------------------------------------------------------!
	452	! 3) Rothrock(1975)'s method
	453	!------------------------------------------------------------------------------!
[2715]	454	IF( kstrngth == 1 ) THEN
	455	z1_3 = 1._wp / 3._wp
[825]	456	DO jl = 1, jpl
	457	DO jj= 1, jpj
	458	DO ji = 1, jpi
[2715]	459	!
[4333]	460	IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN
[5123]	461	zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
[825]	462	!----------------------------
	463	! PE loss from deforming ice
	464	!----------------------------
[5123]	465	strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * zhi * zhi
[825]	466
	467	!--------------------------
	468	! PE gain from rafting ice
	469	!--------------------------
[5123]	470	strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * zhi * zhi
[825]	471
	472	!----------------------------
	473	! PE gain from ridging ice
	474	!----------------------------
[5134]	475	strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) / krdg(ji,jj,jl) &
[5123]	476	* z1_3 * ( hrmax(ji,jj,jl)2 + hrmin(ji,jj,jl)2 + hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )
	477	!!(a3-b3)/(a-b) = aa+ab+bb
	478	ENDIF
[2715]	479	!
[5123]	480	END DO
	481	END DO
	482	END DO
	483
	484	strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:)
	485	! where Cp = (g/2)(rhow-rhoi)(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
[825]	486	ksmooth = 1
	487
[921]	488	!------------------------------------------------------------------------------!
	489	! 4) Hibler (1979)' method
	490	!------------------------------------------------------------------------------!
[825]	491	ELSE ! kstrngth ne 1: Hibler (1979) form
[2715]	492	!
[5123]	493	strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) )
[2715]	494	!
[825]	495	ksmooth = 1
[2715]	496	!
[825]	497	ENDIF ! kstrngth
	498
[921]	499	!
	500	!------------------------------------------------------------------------------!
	501	! 5) Impact of brine volume
	502	!------------------------------------------------------------------------------!
	503	! CAN BE REMOVED
	504	!
[5123]	505	IF( ln_icestr_bvf ) THEN
[825]	506
	507	DO jj = 1, jpj
	508	DO ji = 1, jpi
	509	strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
[5123]	510	END DO
	511	END DO
[825]	512
	513	ENDIF
	514
[921]	515	!
	516	!------------------------------------------------------------------------------!
	517	! 6) Smoothing ice strength
	518	!------------------------------------------------------------------------------!
	519	!
[825]	520	!-------------------
	521	! Spatial smoothing
	522	!-------------------
[2715]	523	IF ( ksmooth == 1 ) THEN
[825]	524
	525	CALL lbc_lnk( strength, 'T', 1. )
	526
[5123]	527	DO jj = 2, jpjm1
	528	DO ji = 2, jpim1
[5167]	529	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN
[5134]	530	zworka(ji,jj) = ( 4.0 * strength(ji,jj) &
	531	& + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &
	532	& + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
	533	& ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
[825]	534	ELSE
[3625]	535	zworka(ji,jj) = 0._wp
[825]	536	ENDIF
	537	END DO
	538	END DO
	539
[5123]	540	DO jj = 2, jpjm1
	541	DO ji = 2, jpim1
[825]	542	strength(ji,jj) = zworka(ji,jj)
	543	END DO
	544	END DO
[869]	545	CALL lbc_lnk( strength, 'T', 1. )
[825]	546
[5123]	547	ENDIF
[825]	548
	549	!--------------------
	550	! Temporal smoothing
	551	!--------------------
[2715]	552	IF ( numit == nit000 + nn_fsbc - 1 ) THEN
[825]	553	strp1(:,:) = 0.0
	554	strp2(:,:) = 0.0
	555	ENDIF
	556
[2715]	557	IF ( ksmooth == 2 ) THEN
[921]	558
	559
[825]	560	CALL lbc_lnk( strength, 'T', 1. )
[921]	561
[825]	562	DO jj = 1, jpj - 1
	563	DO ji = 1, jpi - 1
[5167]	564	IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN
[825]	565	numts_rm = 1 ! number of time steps for the running mean
[5123]	566	IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
	567	IF ( strp2(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
[2715]	568	zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
[825]	569	strp2(ji,jj) = strp1(ji,jj)
	570	strp1(ji,jj) = strength(ji,jj)
	571	strength(ji,jj) = zp
	572
	573	ENDIF
	574	END DO
	575	END DO
	576
	577	ENDIF ! ksmooth
[921]	578
[2715]	579	CALL lbc_lnk( strength, 'T', 1. ) ! Boundary conditions
[825]	580
[3294]	581	CALL wrk_dealloc( jpi, jpj, zworka )
[2715]	582	!
[921]	583	END SUBROUTINE lim_itd_me_icestrength
[825]	584
	585
[2715]	586	SUBROUTINE lim_itd_me_ridgeprep
[921]	587	!!---------------------------------------------------------------------!
	588	!! * ROUTINE lim_itd_me_ridgeprep *
	589	!!
[2715]	590	!! ** Purpose : preparation for ridging and strength calculations
[921]	591	!!
[2715]	592	!! ** Method : Compute the thickness distribution of the ice and open water
	593	!! participating in ridging and of the resulting ridges.
[921]	594	!!---------------------------------------------------------------------!
[2715]	595	INTEGER :: ji,jj, jl ! dummy loop indices
[5123]	596	REAL(wp) :: Gstari, astari, zhi, hrmean, zdummy ! local scalar
[3294]	597	REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here
	598	REAL(wp), POINTER, DIMENSION(:,:,:) :: Gsum ! Gsum(n) = sum of areas in categories 0 to n
	599	!------------------------------------------------------------------------------!
[825]	600
[3294]	601	CALL wrk_alloc( jpi,jpj, zworka )
	602	CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
[825]	603
[5123]	604	Gstari = 1.0/rn_gstar
	605	astari = 1.0/rn_astar
[825]	606	aksum(:,:) = 0.0
	607	athorn(:,:,:) = 0.0
	608	aridge(:,:,:) = 0.0
	609	araft (:,:,:) = 0.0
	610	hrmin(:,:,:) = 0.0
	611	hrmax(:,:,:) = 0.0
	612	hraft(:,:,:) = 0.0
	613	krdg (:,:,:) = 1.0
	614
[921]	615	! ! Zero out categories with very small areas
[5123]	616	CALL lim_var_zapsmall
[825]	617
[921]	618	!------------------------------------------------------------------------------!
	619	! 1) Participation function
	620	!------------------------------------------------------------------------------!
[825]	621
	622	! Compute total area of ice plus open water.
[5134]	623	! This is in general not equal to one because of divergence during transport
	624	asum(:,:) = ato_i(:,:)
	625	DO jl = 1, jpl
	626	asum(:,:) = asum(:,:) + a_i(:,:,jl)
	627	END DO
[825]	628
	629	! Compute cumulative thickness distribution function
	630	! Compute the cumulative thickness distribution function Gsum,
	631	! where Gsum(n) is the fractional area in categories 0 to n.
	632	! initial value (in h = 0) equals open water area
	633
[2715]	634	Gsum(:,:,-1) = 0._wp
[5167]	635	Gsum(:,:,0 ) = ato_i(:,:)
[825]	636
	637	! for each value of h, you have to add ice concentration then
	638	DO jl = 1, jpl
[5167]	639	Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl)
[825]	640	END DO
	641
	642	! Normalize the cumulative distribution to 1
[2715]	643	zworka(:,:) = 1._wp / Gsum(:,:,jpl)
[825]	644	DO jl = 0, jpl
[2715]	645	Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:)
[825]	646	END DO
	647
[921]	648	! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
	649	!--------------------------------------------------------------------------------------------------
	650	! Compute the participation function athorn; this is analogous to
	651	! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
	652	! area lost from category n due to ridging/closing
	653	! athorn(n) = total area lost due to ridging/closing
	654	! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).
	655	!
	656	! The expressions for athorn are found by integrating b(h)g(h) between
	657	! the category boundaries.
	658	!-----------------------------------------------------------------
[825]	659
[5123]	660	IF( nn_partfun == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975)
[4869]	661	DO jl = 0, jpl
[921]	662	DO jj = 1, jpj
	663	DO ji = 1, jpi
[5123]	664	IF( Gsum(ji,jj,jl) < rn_gstar) THEN
[5134]	665	athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * &
	666	& ( 2.0 - (Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari )
[5123]	667	ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN
[5134]	668	athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) * &
	669	& ( 2.0 - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari )
[921]	670	ELSE
	671	athorn(ji,jj,jl) = 0.0
	672	ENDIF
[5123]	673	END DO
	674	END DO
	675	END DO
[825]	676
[2715]	677	ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007)
	678	!
	679	zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array
[921]	680	DO jl = -1, jpl
[2715]	681	Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
[5134]	682	END DO
[4869]	683	DO jl = 0, jpl
[2715]	684	athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
	685	END DO
	686	!
[5134]	687	ENDIF
[825]	688
[5123]	689	IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions
[2715]	690	!
[825]	691	DO jl = 1, jpl
	692	DO jj = 1, jpj
	693	DO ji = 1, jpi
[5123]	694	IF ( athorn(ji,jj,jl) > 0._wp ) THEN
[2715]	695	!!gm TANH( -X ) = - TANH( X ) so can be computed only 1 time....
[5123]	696	aridge(ji,jj,jl) = ( TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
	697	araft (ji,jj,jl) = ( TANH ( -rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
[2715]	698	IF ( araft(ji,jj,jl) < epsi06 ) araft(ji,jj,jl) = 0._wp
	699	aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 )
[5123]	700	ENDIF
	701	END DO
	702	END DO
	703	END DO
[825]	704
[5123]	705	ELSE
[2715]	706	!
[825]	707	DO jl = 1, jpl
[2715]	708	aridge(:,:,jl) = athorn(:,:,jl)
[825]	709	END DO
[2715]	710	!
[825]	711	ENDIF
	712
[5123]	713	IF( ln_rafting ) THEN
[825]	714
[5134]	715	IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 .AND. lwp ) THEN
[921]	716	DO jl = 1, jpl
	717	DO jj = 1, jpj
	718	DO ji = 1, jpi
[5123]	719	IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN
[921]	720	WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
	721	WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
	722	WRITE(numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj)
	723	WRITE(numout,*) ' aridge : ', aridge(ji,jj,1:jpl)
	724	WRITE(numout,*) ' araft : ', araft(ji,jj,1:jpl)
	725	WRITE(numout,*) ' athorn : ', athorn(ji,jj,1:jpl)
	726	ENDIF
	727	END DO
[868]	728	END DO
[825]	729	END DO
[921]	730	ENDIF
[825]	731
	732	ENDIF
	733
[921]	734	!-----------------------------------------------------------------
	735	! 2) Transfer function
	736	!-----------------------------------------------------------------
	737	! Compute max and min ridged ice thickness for each ridging category.
	738	! Assume ridged ice is uniformly distributed between hrmin and hrmax.
	739	!
	740	! This parameterization is a modified version of Hibler (1980).
	741	! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
	742	! and for very thick ridging ice must be >= krdgmin*hi
	743	!
	744	! The minimum ridging thickness, hrmin, is equal to 2*hi
	745	! (i.e., rafting) and for very thick ridging ice is
	746	! constrained by hrmin <= (hrmean + hi)/2.
	747	!
	748	! The maximum ridging thickness, hrmax, is determined by
	749	! hrmean and hrmin.
	750	!
	751	! These modifications have the effect of reducing the ice strength
	752	! (relative to the Hibler formulation) when very thick ice is
	753	! ridging.
	754	!
	755	! aksum = net area removed/ total area removed
	756	! where total area removed = area of ice that ridges
	757	! net area removed = total area removed - area of new ridges
	758	!-----------------------------------------------------------------
[825]	759
	760	! Transfer function
	761	DO jl = 1, jpl !all categories have a specific transfer function
	762	DO jj = 1, jpj
	763	DO ji = 1, jpi
	764
[5123]	765	IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN
	766	zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
	767	hrmean = MAX(SQRT(rn_hstarzhi), zhikrdgmin)
	768	hrmin(ji,jj,jl) = MIN(2.0zhi, 0.5(hrmean + zhi))
[825]	769	hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)
[5123]	770	hraft(ji,jj,jl) = kraft*zhi
	771	krdg(ji,jj,jl) = hrmean / zhi
[825]	772	ELSE
	773	hraft(ji,jj,jl) = 0.0
	774	hrmin(ji,jj,jl) = 0.0
	775	hrmax(ji,jj,jl) = 0.0
	776	krdg (ji,jj,jl) = 1.0
	777	ENDIF
	778
[5134]	779	END DO
	780	END DO
	781	END DO
[825]	782
	783	! Normalization factor : aksum, ensures mass conservation
[2777]	784	aksum(:,:) = athorn(:,:,0)
[825]	785	DO jl = 1, jpl
[2715]	786	aksum(:,:) = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) ) &
	787	& + araft (:,:,jl) * ( 1._wp - 1._wp / kraft )
[825]	788	END DO
[2715]	789	!
[3294]	790	CALL wrk_dealloc( jpi,jpj, zworka )
	791	CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
	792	!
[921]	793	END SUBROUTINE lim_itd_me_ridgeprep
[825]	794
	795
[2715]	796	SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
	797	!!----------------------------------------------------------------------
[921]	798	!! * ROUTINE lim_itd_me_icestrength *
	799	!!
[2715]	800	!! ** Purpose : shift ridging ice among thickness categories of ice thickness
[921]	801	!!
[2715]	802	!! ** Method : Remove area, volume, and energy from each ridging category
	803	!! and add to thicker ice categories.
	804	!!----------------------------------------------------------------------
	805	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear
	806	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area removed, excluding area of new ridges
	807	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: msnow_mlt ! mass of snow added to ocean (kg m-2)
	808	REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)
	809	!
	810	CHARACTER (len=80) :: fieldid ! field identifier
	811	LOGICAL, PARAMETER :: l_conservation_check = .true. ! if true, check conservation (useful for debugging)
	812	!
	813	INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices
	814	INTEGER :: ij ! horizontal index, combines i and j loops
	815	INTEGER :: icells ! number of cells with aicen > puny
[5123]	816	REAL(wp) :: hL, hR, farea, ztmelts ! left and right limits of integration
[825]	817
[3294]	818	INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices
[825]	819
[3294]	820	REAL(wp), POINTER, DIMENSION(:,:) :: vice_init, vice_final ! ice volume summed over categories
	821	REAL(wp), POINTER, DIMENSION(:,:) :: eice_init, eice_final ! ice energy summed over layers
[825]	822
[3294]	823	REAL(wp), POINTER, DIMENSION(:,:,:) :: aicen_init, vicen_init ! ice area & volume before ridging
[4688]	824	REAL(wp), POINTER, DIMENSION(:,:,:) :: vsnwn_init, esnwn_init ! snow volume & energy before ridging
[3294]	825	REAL(wp), POINTER, DIMENSION(:,:,:) :: smv_i_init, oa_i_init ! ice salinity & age before ridging
[825]	826
[3294]	827	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: eicen_init ! ice energy before ridging
[825]	828
[3294]	829	REAL(wp), POINTER, DIMENSION(:,:) :: afrac , fvol ! fraction of category area ridged & new ridge volume going to n2
	830	REAL(wp), POINTER, DIMENSION(:,:) :: ardg1 , ardg2 ! area of ice ridged & new ridges
	831	REAL(wp), POINTER, DIMENSION(:,:) :: vsrdg , esrdg ! snow volume & energy of ridging ice
	832	REAL(wp), POINTER, DIMENSION(:,:) :: oirdg1, oirdg2 ! areal age content of ridged & rifging 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
[825]	841
[3294]	842	REAL(wp), POINTER, DIMENSION(:,:) :: afrft ! fraction of category area rafted
	843	REAL(wp), POINTER, DIMENSION(:,:) :: arft1 , arft2 ! area of ice rafted and new rafted zone
	844	REAL(wp), POINTER, DIMENSION(:,:) :: virft , vsrft ! ice & snow volume of rafting ice
	845	REAL(wp), POINTER, DIMENSION(:,:) :: esrft , smrft ! snow energy & salinity of rafting ice
	846	REAL(wp), POINTER, DIMENSION(:,:) :: oirft1, oirft2 ! areal age content of rafted ice & rafting ice
[825]	847
[3294]	848	REAL(wp), POINTER, DIMENSION(:,:,:) :: eirft ! ice energy of rafting ice
	849	REAL(wp), POINTER, DIMENSION(:,:,:) :: erdg1 ! enth*volume of ice ridged
	850	REAL(wp), POINTER, DIMENSION(:,:,:) :: erdg2 ! enth*volume of new ridges
	851	REAL(wp), POINTER, DIMENSION(:,:,:) :: ersw ! enth of water trapped into ridges
	852	!!----------------------------------------------------------------------
[825]	853
[5181]	854	CALL wrk_alloc( (jpi+1)*(jpj+1), indxi, indxj )
	855	CALL wrk_alloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )
	856	CALL wrk_alloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )
	857	CALL wrk_alloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw )
	858	CALL wrk_alloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
	859	CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
	860	CALL wrk_alloc( jpi, jpj, nlay_i, eirft, erdg1, erdg2, ersw )
	861	CALL wrk_alloc( jpi, jpj, nlay_i, jpl, eicen_init )
[3294]	862
[825]	863	! Conservation check
[2715]	864	eice_init(:,:) = 0._wp
[825]	865
[2715]	866	IF( con_i ) THEN
[4333]	867	CALL lim_column_sum (jpl, v_i, vice_init )
[825]	868	CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_init )
[5128]	869	DO ji = mi0(iiceprt), mi1(iiceprt)
	870	DO jj = mj0(jiceprt), mj1(jiceprt)
[4333]	871	WRITE(numout,*) ' vice_init : ', vice_init(ji,jj)
	872	WRITE(numout,*) ' eice_init : ', eice_init(ji,jj)
	873	END DO
	874	END DO
[825]	875	ENDIF
	876
[921]	877	!-------------------------------------------------------------------------------
	878	! 1) Compute change in open water area due to closing and opening.
	879	!-------------------------------------------------------------------------------
[825]	880	DO jj = 1, jpj
	881	DO ji = 1, jpi
[2715]	882	ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice &
	883	& + opning(ji,jj) * rdt_ice
[5167]	884	IF ( ato_i(ji,jj) < -epsi10 ) THEN ! there is a bug
	885	IF(lwp) WRITE(numout,*) 'Ridging error: ato_i < 0 -- ato_i : ',ato_i(ji,jj)
	886	ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error
[2715]	887	ato_i(ji,jj) = 0._wp
[825]	888	ENDIF
[5134]	889	END DO
	890	END DO
[825]	891
[921]	892	!-----------------------------------------------------------------
	893	! 2) Save initial state variables
	894	!-----------------------------------------------------------------
[5167]	895	aicen_init(:,:,:) = a_i (:,:,:)
	896	vicen_init(:,:,:) = v_i (:,:,:)
	897	vsnwn_init(:,:,:) = v_s (:,:,:)
	898	smv_i_init(:,:,:) = smv_i(:,:,:)
	899	oa_i_init (:,:,:) = oa_i (:,:,:)
	900	esnwn_init(:,:,:) = e_s (:,:,1,:)
	901	eicen_init(:,:,:,:) = e_i (:,:,:,:)
[825]	902
[921]	903	!
	904	!-----------------------------------------------------------------
	905	! 3) Pump everything from ice which is being ridged / rafted
	906	!-----------------------------------------------------------------
	907	! Compute the area, volume, and energy of ice ridging in each
	908	! category, along with the area of the resulting ridge.
[825]	909
	910	DO jl1 = 1, jpl !jl1 describes the ridging category
	911
[921]	912	!------------------------------------------------
	913	! 3.1) Identify grid cells with nonzero ridging
	914	!------------------------------------------------
[825]	915
	916	icells = 0
	917	DO jj = 1, jpj
	918	DO ji = 1, jpi
[4333]	919	IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp &
	920	& .AND. closing_gross(ji,jj) > 0._wp ) THEN
[825]	921	icells = icells + 1
	922	indxi(icells) = ji
	923	indxj(icells) = jj
[5134]	924	ENDIF
	925	END DO
	926	END DO
[825]	927
	928	DO ij = 1, icells
	929	ji = indxi(ij)
	930	jj = indxj(ij)
	931
[921]	932	!--------------------------------------------------------------------
	933	! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
	934	!--------------------------------------------------------------------
[825]	935
	936	ardg1(ji,jj) = aridge(ji,jj,jl1)closing_gross(ji,jj)rdt_ice
	937	arft1(ji,jj) = araft (ji,jj,jl1)closing_gross(ji,jj)rdt_ice
	938	ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1)
	939	arft2(ji,jj) = arft1(ji,jj) / kraft
	940
	941	oirdg1(ji,jj)= aridge(ji,jj,jl1)closing_gross(ji,jj)rdt_ice
	942	oirft1(ji,jj)= araft (ji,jj,jl1)closing_gross(ji,jj)rdt_ice
	943	oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1)
	944	oirft2(ji,jj)= oirft1(ji,jj) / kraft
	945
[921]	946	!---------------------------------------------------------------
	947	! 3.3) Compute ridging /rafting fractions, make sure afrac <=1
	948	!---------------------------------------------------------------
[825]	949
	950	afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging
	951	afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting
	952
[5167]	953	IF( afrac(ji,jj) > kamax + epsi10 ) THEN ! there is a bug
	954	IF(lwp) WRITE(numout,*) ' ardg > a_i -- ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1)
	955	ELSEIF( afrac(ji,jj) > kamax ) THEN ! roundoff error
[4161]	956	afrac(ji,jj) = kamax
[825]	957	ENDIF
[5167]	958
	959	IF( afrft(ji,jj) > kamax + epsi10 ) THEN ! there is a bug
	960	IF(lwp) WRITE(numout,*) ' arft > a_i -- arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)
	961	ELSEIF( afrft(ji,jj) > kamax) THEN ! roundoff error
[4161]	962	afrft(ji,jj) = kamax
[825]	963	ENDIF
	964
[921]	965	!--------------------------------------------------------------------------
	966	! 3.4) Subtract area, volume, and energy from ridging
	967	! / rafting category n1.
	968	!--------------------------------------------------------------------------
[4788]	969	vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj)
[5123]	970	vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + rn_por_rdg )
	971	vsw (ji,jj) = vrdg1(ji,jj) * rn_por_rdg
[825]	972
[4688]	973	vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj)
	974	esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj)
[4788]	975	srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
[825]	976
	977	! rafting volumes, heat contents ...
	978	virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)
[4688]	979	vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj)
	980	esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj)
[825]	981	smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)
	982
	983	! substract everything
	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	oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1) - oirdg1(ji,jj) - oirft1(ji,jj)
	989	smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj) - smrft(ji,jj)
	990
[921]	991	!-----------------------------------------------------------------
	992	! 3.5) Compute properties of new ridges
	993	!-----------------------------------------------------------------
[825]	994	!-------------
	995	! Salinity
	996	!-------------
[4688]	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
[5123]	1000	!srdg2(ji,jj) = MIN( rn_simax * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity
[1572]	1001
[4688]	1002	sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice
[5167]	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
[5123]	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
[5123]	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
[5134]	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
[4788]	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)
[4688]	1061
	1062
	1063	! enthalpy of the trapped seawater (J/m2, >0)
	1064	! clem: if sst>0, then ersw <0 (is that possible?)
[5134]	1065	ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * sst_m(ji,jj) * r1_nlay_i
[825]	1066
[4688]	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
[5123]	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)
[4688]	1072
[5134]	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)
[5134]	1083	END DO
	1084	END DO
[825]	1085	ENDIF
	1086
[921]	1087	!-------------------------------------------------------------------------------
	1088	! 4) Add area, volume, and energy of new ridge to each category jl2
	1089	!-------------------------------------------------------------------------------
	1090	! jl1 looping 1-jpl
[4869]	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
[5134]	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)
[5123]	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
	1120	END DO ! ij
	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)
[5134]	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
[4869]	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	!
[5134]	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)
[5123]	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)
	1147	oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirft2(ji,jj)
[5123]	1148	ENDIF
[3625]	1149	!
[5123]	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)
[5134]	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
[5123]	1160	END DO
	1161	END DO
[825]	1162
[5123]	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
[5128]	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	!
[5181]	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, oirdg1, oirdg2, dhr, dhr2 )
	1190	CALL wrk_dealloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw )
	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
[5123]	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,*)' ~~~~~~~~~~~~~~~'
[5123]	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: