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

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 7698

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

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

Download in other formats: