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/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 3649

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

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

Download in other formats: