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.

icerdgrft.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/icerdgrft.F90 @ 8498

Last change on this file since 8498 was 8498, checked in by clem, 7 years ago
changes in style - part2 -
File size: 47.8 KB

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

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

Download in other formats: