Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

limitd_me.F90 @ 8554

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

Download in other formats: