Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

limitd_me.F90 @ 5715

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

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

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

Context Navigation

source: branches/NERC/dev_r5107_NOC_MEDUSA/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 5715

Download in other formats: