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

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

limitd_me.F90 in branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 4635

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

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

Download in other formats: