New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

limitd_me.F90 in branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 6584

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

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

Download in other formats: