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

Context Navigation

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90 @ 6400

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

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

Download in other formats: