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/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 8407

Last change on this file since 8407 was 8378, checked in by clem, 7 years ago
some cleaning
Property svn:keywords set to `Id`
File size: 48.7 KB

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

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

Download in other formats: