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 @ 6796

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

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

Download in other formats: