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/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 8142

Last change on this file since 8142 was 8142, checked in by vancop, 7 years ago
Melt pond interfaces practically operational
Property svn:keywords set to `Id`
File size: 50.6 KB

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

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

Download in other formats: