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.

icedyn_rdgrft.F90 in NEMO/branches/UKMO/NEMO_4.0.1_mirror/src/ICE – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.1_mirror/src/ICE/icedyn_rdgrft.F90 @ 11715

Last change on this file since 11715 was 11715, checked in by davestorkey, 4 years ago
UKMO/NEMO_4.0.1_mirror : Remove SVN keywords.
File size: 52.0 KB

Rev	Line
[8586]	1	MODULE icedyn_rdgrft
	2	!!======================================================================
	3	!! * MODULE icedyn_rdgrft *
	4	!! sea-ice : Mechanical impact on ice thickness distribution
	5	!!======================================================================
[9604]	6	!! History : ! 2006-02 (M. Vancoppenolle) Original code
	7	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
[8586]	8	!!----------------------------------------------------------------------
[9570]	9	#if defined key_si3
[8586]	10	!!----------------------------------------------------------------------
[9570]	11	!! 'key_si3' SI3 sea-ice model
[8586]	12	!!----------------------------------------------------------------------
	13	!! ice_dyn_rdgrft : ridging/rafting of sea ice
	14	!! ice_dyn_rdgrft_init : initialization of ridging/rafting of sea ice
	15	!! ice_strength : ice strength calculation
	16	!!----------------------------------------------------------------------
	17	USE dom_oce ! ocean domain
	18	USE phycst ! physical constants (ocean directory)
	19	USE sbc_oce , ONLY : sss_m, sst_m ! surface boundary condition: ocean fields
	20	USE ice1D ! sea-ice: thermodynamics
	21	USE ice ! sea-ice: variables
	22	USE icetab ! sea-ice: 1D <==> 2D transformation
	23	USE icevar ! sea-ice: operations
	24	USE icectl ! sea-ice: control prints
	25	!
	26	USE in_out_manager ! I/O manager
	27	USE iom ! I/O manager library
	28	USE lib_mpp ! MPP library
	29	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	30	USE lbclnk ! lateral boundary conditions (or mpp links)
	31	USE timing ! Timing
	32
	33	IMPLICIT NONE
	34	PRIVATE
	35
	36	PUBLIC ice_dyn_rdgrft ! called by icestp
	37	PUBLIC ice_dyn_rdgrft_init ! called by icedyn
	38	PUBLIC ice_strength ! called by icedyn_rhg_evp
	39
	40	! Variables shared among ridging subroutines
[10413]	41	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: closing_net ! net rate at which area is removed (1/s)
	42	! ! (ridging ice area - area of new ridges) / dt
	43	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: opning ! rate of opening due to divergence/shear
	44	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: closing_gross ! rate at which area removed, not counting area of new ridges
	45	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: apartf ! participation function; fraction of ridging/closing associated w/ category n
	46	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmin ! minimum ridge thickness
	47	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmax ! maximum ridge thickness
	48	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hraft ! thickness of rafted ice
	49	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hi_hrdg ! thickness of ridging ice / mean ridge thickness
	50	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: aridge ! participating ice ridging
	51	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: araft ! participating ice rafting
[8586]	52	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ze_i_2d
	53	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ze_s_2d
	54	!
	55	REAL(wp), PARAMETER :: hrdg_hi_min = 1.1_wp ! min ridge thickness multiplier: min(hrdg/hi)
	56	REAL(wp), PARAMETER :: hi_hrft = 0.5_wp ! rafting multipliyer: (hi/hraft)
	57	!
	58	! namelist (namdyn_rdgrft)
	59	LOGICAL :: ln_str_H79 ! ice strength parameterization (Hibler79)
	60	REAL(wp) :: rn_pstar ! determines ice strength, Hibler JPO79
	61	REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging
	62	LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975))
	63	REAL(wp) :: rn_gstar ! fractional area of young ice contributing to ridging
	64	LOGICAL :: ln_partf_exp ! participation function exponential (Lipscomb et al. (2007))
	65	REAL(wp) :: rn_astar ! equivalent of G* for an exponential participation function
	66	LOGICAL :: ln_ridging ! ridging of ice or not
	67	REAL(wp) :: rn_hstar ! thickness that determines the maximal thickness of ridged ice
	68	REAL(wp) :: rn_porordg ! initial porosity of ridges (0.3 regular value)
	69	REAL(wp) :: rn_fsnwrdg ! fractional snow loss to the ocean during ridging
	70	REAL(wp) :: rn_fpndrdg ! fractional pond loss to the ocean during ridging
	71	LOGICAL :: ln_rafting ! rafting of ice or not
	72	REAL(wp) :: rn_hraft ! threshold thickness (m) for rafting / ridging
	73	REAL(wp) :: rn_craft ! coefficient for smoothness of the hyperbolic tangent in rafting
	74	REAL(wp) :: rn_fsnwrft ! fractional snow loss to the ocean during rafting
	75	REAL(wp) :: rn_fpndrft ! fractional pond loss to the ocean during rafting
	76	!
	77	!!----------------------------------------------------------------------
[9598]	78	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	79	!! $Id$
[10413]	80	!! Software governed by the CeCILL licence (./LICENSE)
[8586]	81	!!----------------------------------------------------------------------
	82	CONTAINS
	83
	84	INTEGER FUNCTION ice_dyn_rdgrft_alloc()
	85	!!-------------------------------------------------------------------
	86	!! * ROUTINE ice_dyn_rdgrft_alloc *
	87	!!-------------------------------------------------------------------
	88	ALLOCATE( closing_net(jpij), opning(jpij) , closing_gross(jpij), &
	89	& apartf(jpij,0:jpl), hrmin(jpij,jpl), hraft(jpij,jpl) , aridge(jpij,jpl), &
	90	& hrmax(jpij,jpl), hi_hrdg(jpij,jpl) , araft (jpij,jpl), &
	91	& ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc )
	92
[10425]	93	CALL mpp_sum ( 'icedyn_rdgrft', ice_dyn_rdgrft_alloc )
	94	IF( ice_dyn_rdgrft_alloc /= 0 ) CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_alloc: failed to allocate arrays' )
[8586]	95	!
	96	END FUNCTION ice_dyn_rdgrft_alloc
	97
	98
	99	SUBROUTINE ice_dyn_rdgrft( kt )
	100	!!-------------------------------------------------------------------
	101	!! * ROUTINE ice_dyn_rdgrft *
	102	!!
	103	!! ** Purpose : computes the mechanical redistribution of ice thickness
	104	!!
	105	!! ** Method : Steps :
[9452]	106	!! 0) Identify grid cells with ice
	107	!! 1) Calculate closing rate, divergence and opening
	108	!! 2) Identify grid cells with ridging
	109	!! 3) Start ridging iterations
	110	!! - prep = ridged and rafted ice + closing_gross
	111	!! - shift = move ice from one category to another
[8586]	112	!!
[9452]	113	!! ** Details
	114	!! step1: The net rate of closing is due to convergence and shear, based on Flato and Hibler (1995).
	115	!! The energy dissipation rate is equal to the net closing rate times the ice strength.
	116	!!
	117	!! step3: The gross closing rate is equal to the first two terms (open
	118	!! water closing and thin ice ridging) without the third term
	119	!! (thick, newly ridged ice).
	120	!!
[8586]	121	!! References : Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626.
	122	!! Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980.
	123	!! Rothrock, D. A., 1975: JGR, 80, 4514-4519.
	124	!! Thorndike et al., 1975, JGR, 80, 4501-4513.
	125	!! Bitz et al., JGR, 2001
	126	!! Amundrud and Melling, JGR 2005
	127	!! Babko et al., JGR 2002
	128	!!
	129	!! This routine is based on CICE code and authors William H. Lipscomb,
[9452]	130	!! and Elizabeth C. Hunke, LANL are gratefully acknowledged
[8586]	131	!!-------------------------------------------------------------------
	132	INTEGER, INTENT(in) :: kt ! number of iteration
	133	!!
	134	INTEGER :: ji, jj, jk, jl ! dummy loop index
	135	INTEGER :: iter, iterate_ridging ! local integer
	136	INTEGER :: ipti ! local integer
	137	REAL(wp) :: zfac ! local scalar
	138	INTEGER , DIMENSION(jpij) :: iptidx ! compute ridge/raft or not
	139	REAL(wp), DIMENSION(jpij) :: zdivu_adv ! divu as implied by transport scheme (1/s)
	140	REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i
	141	!
	142	INTEGER, PARAMETER :: jp_itermax = 20
	143	!!-------------------------------------------------------------------
	144	! controls
[9124]	145	IF( ln_timing ) CALL timing_start('icedyn_rdgrft') ! timing
	146	IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
[11536]	147	IF( ln_icediachk ) CALL ice_cons2D (0, 'icedyn_rdgrft', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation
[8586]	148
	149	IF( kt == nit000 ) THEN
	150	IF(lwp) WRITE(numout,*)
	151	IF(lwp) WRITE(numout,*)'ice_dyn_rdgrft: ice ridging and rafting'
	152	IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~~'
	153	ENDIF
	154
	155	!--------------------------------
	156	! 0) Identify grid cells with ice
	157	!--------------------------------
[10994]	158	at_i(:,:) = SUM( a_i, dim=3 )
	159	!
[8586]	160	npti = 0 ; nptidx(:) = 0
[9452]	161	ipti = 0 ; iptidx(:) = 0
[8586]	162	DO jj = 1, jpj
	163	DO ji = 1, jpi
[10994]	164	IF ( at_i(ji,jj) > epsi10 ) THEN
[8586]	165	npti = npti + 1
	166	nptidx( npti ) = (jj - 1) * jpi + ji
	167	ENDIF
	168	END DO
	169	END DO
[9452]	170
	171	!--------------------------------------------------------
	172	! 1) Dynamical inputs (closing rate, divergence, opening)
	173	!--------------------------------------------------------
[8586]	174	IF( npti > 0 ) THEN
	175
	176	! just needed here
[10994]	177	CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i )
	178	CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i )
[8586]	179	! needed here and in the iteration loop
[10994]	180	CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i )
	181	CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i )
	182	CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i )
[8586]	183
	184	DO ji = 1, npti
[9452]	185	! closing_net = rate at which open water area is removed + ice area removed by ridging
	186	! - ice area added in new ridges
[8586]	187	closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
	188	!
[9452]	189	! divergence given by the advection scheme
	190	! (which may not be equal to divu as computed from the velocity field)
[10413]	191	IF ( ln_adv_Pra ) THEN
	192	zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice
	193	ELSEIF( ln_adv_UMx ) THEN
	194	zdivu_adv(ji) = zdivu(ji)
	195	ENDIF
[9452]	196	!
	197	IF( zdivu_adv(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) ) ! make sure the closing rate is large enough
	198	! ! to give asum = 1.0 after ridging
	199	! Opening rate (non-negative) that will give asum = 1.0 after ridging.
[8586]	200	opning(ji) = closing_net(ji) + zdivu_adv(ji)
	201	END DO
[9452]	202	!
	203	!------------------------------------
	204	! 2) Identify grid cells with ridging
	205	!------------------------------------
	206	CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
[8586]	207
	208	DO ji = 1, npti
[9452]	209	IF( SUM( apartf(ji,1:jpl) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
[8586]	210	ipti = ipti + 1
	211	iptidx (ipti) = nptidx (ji)
[9452]	212	! adjust to new indices
	213	a_i_2d (ipti,:) = a_i_2d (ji,:)
	214	v_i_2d (ipti,:) = v_i_2d (ji,:)
	215	ato_i_1d (ipti) = ato_i_1d (ji)
	216	closing_net(ipti) = closing_net(ji)
	217	zdivu_adv (ipti) = zdivu_adv (ji)
	218	opning (ipti) = opning (ji)
[8586]	219	ENDIF
	220	END DO
	221
	222	ENDIF
	223
[9452]	224	! grid cells with ridging
	225	nptidx(:) = iptidx(:)
	226	npti = ipti
	227
	228	!-----------------
	229	! 3) Start ridging
	230	!-----------------
[8586]	231	IF( npti > 0 ) THEN
	232
[9452]	233	CALL ice_dyn_1d2d( 1 ) ! --- Move to 1D arrays --- !
	234
[8586]	235	iter = 1
	236	iterate_ridging = 1
[9452]	237	! !----------------------!
	238	DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax ) ! ridging iterations !
	239	! !----------------------!
	240	! Calculate participation function (apartf)
	241	! and transfer function
	242	! and closing_gross (+correction on opening)
	243	CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net )
[8586]	244
[9452]	245	! Redistribute area, volume, and energy between categories
[8586]	246	CALL rdgrft_shift
	247
[9452]	248	! Do we keep on iterating?
	249	!-------------------------
	250	! Check whether a_i + ato_i = 0
	251	! If not, because the closing and opening rates were reduced above, ridge again with new rates
[8586]	252	iterate_ridging = 0
	253	DO ji = 1, npti
	254	zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) )
	255	IF( ABS( zfac ) < epsi10 ) THEN
	256	closing_net(ji) = 0._wp
	257	opning (ji) = 0._wp
	258	ato_i_1d (ji) = MAX( 0._wp, 1._wp - SUM( a_i_2d(ji,:) ) )
	259	ELSE
	260	iterate_ridging = 1
	261	zdivu_adv (ji) = zfac * r1_rdtice
	262	closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) )
	263	opning (ji) = MAX( 0._wp, zdivu_adv(ji) )
	264	ENDIF
	265	END DO
	266	!
	267	iter = iter + 1
[10425]	268	IF( iter > jp_itermax ) CALL ctl_stop( 'STOP', 'icedyn_rdgrft: non-converging ridging scheme' )
[8586]	269	!
	270	END DO
	271
[9452]	272	CALL ice_dyn_1d2d( 2 ) ! --- Move to 2D arrays --- !
[8586]	273
[9452]	274	ENDIF
[8586]	275
	276	CALL ice_var_agg( 1 )
	277
	278	! controls
[11536]	279	IF( ln_ctl ) CALL ice_prt3D ('icedyn_rdgrft') ! prints
[11560]	280	IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt,-1, ' - ice dyn rdgrft - ') ! prints
[9124]	281	IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
[11536]	282	IF( ln_icediachk ) CALL ice_cons2D (1, 'icedyn_rdgrft', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation
[9124]	283	IF( ln_timing ) CALL timing_stop ('icedyn_rdgrft') ! timing
[8586]	284	!
	285	END SUBROUTINE ice_dyn_rdgrft
	286
	287
[9452]	288	SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i, pclosing_net )
[8586]	289	!!-------------------------------------------------------------------
	290	!! * ROUTINE rdgrft_prep *
	291	!!
	292	!! ** Purpose : preparation for ridging calculations
	293	!!
	294	!! ** Method : Compute the thickness distribution of the ice and open water
[9452]	295	!! participating in ridging and of the resulting ridges.
[8586]	296	!!-------------------------------------------------------------------
[9452]	297	REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i, pclosing_net
[8586]	298	REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i, pv_i
	299	!!
	300	INTEGER :: ji, jl ! dummy loop indices
	301	REAL(wp) :: z1_gstar, z1_astar, zhmean, zfac ! local scalar
	302	REAL(wp), DIMENSION(jpij) :: zasum, z1_asum, zaksum ! sum of a_i+ato_i and reverse
	303	REAL(wp), DIMENSION(jpij,jpl) :: zhi ! ice thickness
	304	REAL(wp), DIMENSION(jpij,-1:jpl) :: zGsum ! zGsum(n) = sum of areas in categories 0 to n
	305	!--------------------------------------------------------------------
	306
	307	z1_gstar = 1._wp / rn_gstar
	308	z1_astar = 1._wp / rn_astar
	309
	310	! ! Ice thickness needed for rafting
[10994]	311	WHERE( pa_i(1:npti,:) > epsi20 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
	312	ELSEWHERE ; zhi(1:npti,:) = 0._wp
[8586]	313	END WHERE
	314
[9452]	315	! 1) Participation function (apartf): a(h) = b(h).g(h)
[8586]	316	!-----------------------------------------------------------------
[9452]	317	! Compute the participation function = total area lost due to ridging/closing
	318	! This is analogous to
	319	! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
	320	! assuming b(h) = (2/Gstar) * (1 - G(h)/Gstar).
[8586]	321	!
[9452]	322	! apartf = integrating b(h)g(h) between the category boundaries
[8586]	323	! apartf is always >= 0 and SUM(apartf(0:jpl))=1
	324	!-----------------------------------------------------------------
	325	!
	326	! Compute total area of ice plus open water.
	327	! This is in general not equal to one because of divergence during transport
	328	zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
	329	!
[10994]	330	WHERE( zasum(1:npti) > epsi20 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti)
	331	ELSEWHERE ; z1_asum(1:npti) = 0._wp
[8586]	332	END WHERE
[9452]	333	!
[8586]	334	! Compute cumulative thickness distribution function
	335	! Compute the cumulative thickness distribution function zGsum,
	336	! where zGsum(n) is the fractional area in categories 0 to n.
[9452]	337	! initial value (in h = 0) = open water area
[8586]	338	zGsum(1:npti,-1) = 0._wp
	339	zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti)
	340	DO jl = 1, jpl
	341	zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti) ! sum(1:jl) is ok (and not jpl)
	342	END DO
	343	!
	344	IF( ln_partf_lin ) THEN !--- Linear formulation (Thorndike et al., 1975)
	345	DO jl = 0, jpl
	346	DO ji = 1, npti
	347	IF ( zGsum(ji,jl) < rn_gstar ) THEN
	348	apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * &
	349	& ( 2._wp - ( zGsum(ji,jl-1) + zGsum(ji,jl) ) * z1_gstar )
	350	ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN
	351	apartf(ji,jl) = z1_gstar * ( rn_gstar - zGsum(ji,jl-1) ) * &
	352	& ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar ) * z1_gstar )
	353	ELSE
	354	apartf(ji,jl) = 0._wp
	355	ENDIF
	356	END DO
	357	END DO
	358	!
	359	ELSEIF( ln_partf_exp ) THEN !--- Exponential, more stable formulation (Lipscomb et al, 2007)
	360	!
	361	zfac = 1._wp / ( 1._wp - EXP(-z1_astar) )
	362	DO jl = -1, jpl
	363	DO ji = 1, npti
	364	zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac
	365	END DO
	366	END DO
	367	DO jl = 0, jpl
	368	DO ji = 1, npti
	369	apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl)
	370	END DO
	371	END DO
	372	!
	373	ENDIF
	374
	375	! !--- Ridging and rafting participation concentrations
	376	IF( ln_rafting .AND. ln_ridging ) THEN !- ridging & rafting
	377	DO jl = 1, jpl
	378	DO ji = 1, npti
	379	aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl)
	380	araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl)
	381	END DO
	382	END DO
	383	ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN !- ridging alone
	384	DO jl = 1, jpl
	385	DO ji = 1, npti
	386	aridge(ji,jl) = apartf(ji,jl)
	387	araft (ji,jl) = 0._wp
	388	END DO
	389	END DO
	390	ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN !- rafting alone
	391	DO jl = 1, jpl
	392	DO ji = 1, npti
	393	aridge(ji,jl) = 0._wp
	394	araft (ji,jl) = apartf(ji,jl)
	395	END DO
	396	END DO
	397	ELSE !- no ridging & no rafting
	398	DO jl = 1, jpl
	399	DO ji = 1, npti
	400	aridge(ji,jl) = 0._wp
	401	araft (ji,jl) = 0._wp
	402	END DO
	403	END DO
	404	ENDIF
	405
	406	! 2) Transfer function
	407	!-----------------------------------------------------------------
	408	! Compute max and min ridged ice thickness for each ridging category.
	409	! Assume ridged ice is uniformly distributed between hrmin and hrmax.
	410	!
	411	! This parameterization is a modified version of Hibler (1980).
	412	! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
	413	! and for very thick ridging ice must be >= hrdg_hi_min*hi
	414	!
	415	! The minimum ridging thickness, hrmin, is equal to 2*hi
	416	! (i.e., rafting) and for very thick ridging ice is
	417	! constrained by hrmin <= (zhmean + hi)/2.
	418	!
[9452]	419	! The maximum ridging thickness, hrmax, is determined by zhmean and hrmin.
[8586]	420	!
	421	! These modifications have the effect of reducing the ice strength
[9452]	422	! (relative to the Hibler formulation) when very thick ice is ridging.
[8586]	423	!
	424	! zaksum = net area removed/ total area removed
	425	! where total area removed = area of ice that ridges
	426	! net area removed = total area removed - area of new ridges
	427	!-----------------------------------------------------------------
	428	zfac = 1._wp / hi_hrft
	429	zaksum(1:npti) = apartf(1:npti,0)
[9452]	430	!
	431	DO jl = 1, jpl
[8586]	432	DO ji = 1, npti
	433	IF ( apartf(ji,jl) > 0._wp ) THEN
	434	zhmean = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min )
	435	hrmin (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) )
	436	hrmax (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)
	437	hraft (ji,jl) = zhi(ji,jl) * zfac
	438	hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 )
	439	!
	440	! Normalization factor : zaksum, ensures mass conservation
	441	zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) ) &
	442	& + araft (ji,jl) * ( 1._wp - hi_hrft )
	443	ELSE
	444	hrmin (ji,jl) = 0._wp
	445	hrmax (ji,jl) = 0._wp
	446	hraft (ji,jl) = 0._wp
	447	hi_hrdg(ji,jl) = 1._wp
	448	ENDIF
	449	END DO
	450	END DO
	451	!
[9452]	452	! 3) closing_gross
	453	!-----------------
	454	! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.
[8586]	455	! NOTE: 0 < aksum <= 1
[10994]	456	WHERE( zaksum(1:npti) > epsi20 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
	457	ELSEWHERE ; closing_gross(1:npti) = 0._wp
[8586]	458	END WHERE
	459
[9452]	460	! correction to closing rate if excessive ice removal
	461	!----------------------------------------------------
	462	! Reduce the closing rate if more than 100% of any ice category would be removed
	463	! Reduce the opening rate in proportion
[8586]	464	DO jl = 1, jpl
	465	DO ji = 1, npti
	466	zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice
[10994]	467	IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN
[9452]	468	closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice
[8586]	469	ENDIF
	470	END DO
	471	END DO
[9452]	472
	473	! 4) correction to opening if excessive open water removal
	474	!---------------------------------------------------------
	475	! Reduce the closing rate if more than 100% of the open water would be removed
	476	! Reduce the opening rate in proportion
	477	DO ji = 1, npti
	478	zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice
	479	IF( zfac < 0._wp ) THEN ! would lead to negative ato_i
	480	opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_rdtice
	481	ELSEIF( zfac > zasum(ji) ) THEN ! would lead to ato_i > asum
	482	opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_rdtice
	483	ENDIF
	484	END DO
[8586]	485	!
	486	END SUBROUTINE rdgrft_prep
	487
	488
	489	SUBROUTINE rdgrft_shift
	490	!!-------------------------------------------------------------------
	491	!! * ROUTINE rdgrft_shift *
	492	!!
	493	!! ** Purpose : shift ridging ice among thickness categories of ice thickness
	494	!!
	495	!! ** Method : Remove area, volume, and energy from each ridging category
[9452]	496	!! and add to thicker ice categories.
[8586]	497	!!-------------------------------------------------------------------
	498	!
	499	INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices
	500	REAL(wp) :: hL, hR, farea ! left and right limits of integration and new area going to jl2
	501	REAL(wp) :: vsw ! vol of water trapped into ridges
	502	REAL(wp) :: afrdg, afrft ! fraction of category area ridged/rafted
	503	REAL(wp) :: airdg1, oirdg1, aprdg1, virdg1, sirdg1
	504	REAL(wp) :: airft1, oirft1, aprft1
[9271]	505	REAL(wp), DIMENSION(jpij) :: airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg ! area etc of new ridges
	506	REAL(wp), DIMENSION(jpij) :: airft2, oirft2, aprft2, virft , sirft , vsrft, vprft ! area etc of rafted ice
[8586]	507	!
	508	REAL(wp), DIMENSION(jpij) :: ersw ! enth of water trapped into ridges
	509	REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to jl2
	510	REAL(wp), DIMENSION(jpij) :: z1_ai ! 1 / a
[10994]	511	REAL(wp), DIMENSION(jpij) :: zvti ! sum(v_i)
[8586]	512	!
[9271]	513	REAL(wp), DIMENSION(jpij,nlay_s) :: esrft ! snow energy of rafting ice
	514	REAL(wp), DIMENSION(jpij,nlay_i) :: eirft ! ice energy of rafting ice
	515	REAL(wp), DIMENSION(jpij,nlay_s) :: esrdg ! enth*volume of new ridges
[8586]	516	REAL(wp), DIMENSION(jpij,nlay_i) :: eirdg ! enth*volume of new ridges
[9452]	517	!
	518	INTEGER , DIMENSION(jpij) :: itest_rdg, itest_rft ! test for conservation
[8586]	519	!!-------------------------------------------------------------------
[10994]	520	!
	521	zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 ) ! total ice volume
	522	!
[9452]	523	! 1) Change in open water area due to closing and opening
	524	!--------------------------------------------------------
[8586]	525	DO ji = 1, npti
	526	ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice )
	527	END DO
	528
[9452]	529	! 2) compute categories in which ice is removed (jl1)
	530	!----------------------------------------------------
[8586]	531	DO jl1 = 1, jpl
	532
	533	CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) )
	534
	535	DO ji = 1, npti
	536
[9452]	537	IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ! only if ice is ridging
[8586]	538
[10994]	539	IF( a_i_2d(ji,jl1) > epsi20 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
	540	ELSE ; z1_ai(ji) = 0._wp
	541	ENDIF
	542
[9452]	543	! area of ridging / rafting ice (airdg1) and of new ridge (airdg2)
[8586]	544	airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice
	545	airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice
	546
	547	airdg2(ji) = airdg1 * hi_hrdg(ji,jl1)
	548	airft2(ji) = airft1 * hi_hrft
	549
[9452]	550	! ridging /rafting fractions
[8586]	551	afrdg = airdg1 * z1_ai(ji)
	552	afrft = airft1 * z1_ai(ji)
	553
	554	! volume and enthalpy (J/m2, >0) of seawater trapped into ridges
[10994]	555	IF ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg ! v <= 10m then porosity = rn_porordg
	556	ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp ! v >= 20m then porosity = 0
	557	ELSE ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg * MAX( 0._wp, 2._wp - 0.1_wp * zvti(ji) ) ! v > 10m and v < 20m then porosity = linear transition to 0
	558	ENDIF
[9935]	559	ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji) ! clem: if sst>0, then ersw <0 (is that possible?)
[8586]	560
[9452]	561	! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei)
[8586]	562	virdg1 = v_i_2d (ji,jl1) * afrdg
[10994]	563	virdg2(ji) = v_i_2d (ji,jl1) * afrdg + vsw
[8586]	564	vsrdg(ji) = v_s_2d (ji,jl1) * afrdg
	565	sirdg1 = sv_i_2d(ji,jl1) * afrdg
	566	sirdg2(ji) = sv_i_2d(ji,jl1) * afrdg + vsw * sss_1d(ji)
	567	oirdg1 = oa_i_2d(ji,jl1) * afrdg
	568	oirdg2(ji) = oa_i_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
	569
	570	virft(ji) = v_i_2d (ji,jl1) * afrft
	571	vsrft(ji) = v_s_2d (ji,jl1) * afrft
	572	sirft(ji) = sv_i_2d(ji,jl1) * afrft
	573	oirft1 = oa_i_2d(ji,jl1) * afrft
	574	oirft2(ji) = oa_i_2d(ji,jl1) * afrft * hi_hrft
	575
[8637]	576	IF ( ln_pnd_H12 ) THEN
[8586]	577	aprdg1 = a_ip_2d(ji,jl1) * afrdg
	578	aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
	579	vprdg (ji) = v_ip_2d(ji,jl1) * afrdg
	580	aprft1 = a_ip_2d(ji,jl1) * afrft
	581	aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft
	582	vprft (ji) = v_ip_2d(ji,jl1) * afrft
	583	ENDIF
	584
[9452]	585	! Ice-ocean exchanges associated with ice porosity
[9935]	586	wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoi * r1_rdtice ! increase in ice volume due to seawater frozen in voids
	587	sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoi * r1_rdtice
[8586]	588	hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_rdtice ! > 0 [W.m-2]
	589
	590	! Put the snow lost by ridging into the ocean
	591	! Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow.
[9935]	592	wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhos * vsrdg(ji) * ( 1._wp - rn_fsnwrdg ) & ! fresh water source for ocean
	593	& + rhos * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
[8586]	594
	595	! virtual salt flux to keep salinity constant
	596	IF( nn_icesal /= 2 ) THEN
	597	sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - s_i_1d(ji) ) ! ridge salinity = s_i
[9935]	598	sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_rdtice & ! put back sss_m into the ocean
	599	& - s_i_1d(ji) * vsw * rhoi * r1_rdtice ! and get s_i from the ocean
[8586]	600	ENDIF
	601
[9452]	602	! Remove area, volume of new ridge to each category jl1
	603	!------------------------------------------------------
[8586]	604	a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1 - airft1
	605	v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1 - virft(ji)
	606	v_s_2d (ji,jl1) = v_s_2d (ji,jl1) - vsrdg(ji) - vsrft(ji)
	607	sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1 - sirft(ji)
	608	oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1 - oirft1
[8637]	609	IF ( ln_pnd_H12 ) THEN
[8586]	610	a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1 - aprft1
	611	v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
	612	ENDIF
	613	ENDIF
	614
	615	END DO ! ji
	616
[9271]	617	! special loop for e_s because of layers jk
	618	DO jk = 1, nlay_s
	619	DO ji = 1, npti
	620	IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
	621	! Compute ridging /rafting fractions
	622	afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
	623	afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
[9452]	624	! Compute ridging /rafting ice and new ridges for es
[9271]	625	esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg
	626	esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft
	627	! Put the snow lost by ridging into the ocean
	628	hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg ) & ! heat sink for ocean (<0, W.m-2)
	629	& - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice
[9452]	630	!
	631	! Remove energy of new ridge to each category jl1
	632	!-------------------------------------------------
[9271]	633	ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft )
	634	ENDIF
	635	END DO
	636	END DO
	637
[8586]	638	! special loop for e_i because of layers jk
	639	DO jk = 1, nlay_i
	640	DO ji = 1, npti
	641	IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
	642	! Compute ridging /rafting fractions
	643	afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
	644	afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji)
	645	! Compute ridging ice and new ridges for ei
	646	eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i
	647	eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft
[9452]	648	!
	649	! Remove energy of new ridge to each category jl1
	650	!-------------------------------------------------
[8586]	651	ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft )
	652	ENDIF
	653	END DO
	654	END DO
	655
[9452]	656	! 3) compute categories in which ice is added (jl2)
	657	!--------------------------------------------------
	658	itest_rdg(1:npti) = 0
	659	itest_rft(1:npti) = 0
[8586]	660	DO jl2 = 1, jpl
	661	!
	662	DO ji = 1, npti
	663
	664	IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
	665
[9452]	666	! Compute the fraction of ridged ice area and volume going to thickness category jl2
[8586]	667	IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN
	668	hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) )
	669	hR = MIN( hrmax(ji,jl1), hi_max(jl2) )
	670	farea = ( hR - hL ) / ( hrmax(ji,jl1) - hrmin(ji,jl1) )
	671	fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) )
[9452]	672	!
	673	itest_rdg(ji) = 1 ! test for conservation
[8586]	674	ELSE
	675	farea = 0._wp
	676	fvol(ji) = 0._wp
	677	ENDIF
	678
	679	! Compute the fraction of rafted ice area and volume going to thickness category jl2
[9452]	680	IF( hraft(ji,jl1) <= hi_max(jl2) .AND. hraft(ji,jl1) > hi_max(jl2-1) ) THEN
	681	zswitch(ji) = 1._wp
	682	!
	683	itest_rft(ji) = 1 ! test for conservation
	684	ELSE
	685	zswitch(ji) = 0._wp
[8586]	686	ENDIF
	687	!
[9452]	688	! Patch to ensure perfect conservation if ice thickness goes mad
	689	! Sometimes thickness is larger than hi_max(jpl) because of advection scheme (for very small areas)
	690	! Then ice volume is removed from one category but the ridging/rafting scheme
	691	! does not know where to move it, leading to a conservation issue.
	692	IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN ; farea = 1._wp ; fvol(ji) = 1._wp ; ENDIF
	693	IF( itest_rft(ji) == 0 .AND. jl2 == jpl ) zswitch(ji) = 1._wp
	694	!
	695	! Add area, volume of new ridge to category jl2
	696	!----------------------------------------------
[8586]	697	a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea + airft2(ji) * zswitch(ji) )
	698	oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea + oirft2(ji) * zswitch(ji) )
	699	v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) )
	700	sv_i_2d(ji,jl2) = sv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) )
	701	v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji) + &
[9271]	702	& vsrft (ji) * rn_fsnwrft * zswitch(ji) )
[8637]	703	IF ( ln_pnd_H12 ) THEN
[8586]	704	v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + ( vprdg (ji) * rn_fpndrdg * fvol (ji) &
	705	& + vprft (ji) * rn_fpndrft * zswitch(ji) )
	706	a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + ( aprdg2(ji) * rn_fpndrdg * farea &
	707	& + aprft2(ji) * rn_fpndrft * zswitch(ji) )
	708	ENDIF
	709
	710	ENDIF
	711
	712	END DO
[9452]	713	! Add snow energy of new ridge to category jl2
	714	!---------------------------------------------
[9271]	715	DO jk = 1, nlay_s
	716	DO ji = 1, npti
	717	IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) &
	718	& ze_s_2d(ji,jk,jl2) = ze_s_2d(ji,jk,jl2) + ( esrdg(ji,jk) * rn_fsnwrdg * fvol(ji) + &
	719	& esrft(ji,jk) * rn_fsnwrft * zswitch(ji) )
	720	END DO
	721	END DO
[9452]	722	! Add ice energy of new ridge to category jl2
	723	!--------------------------------------------
[8586]	724	DO jk = 1, nlay_i
	725	DO ji = 1, npti
	726	IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) &
	727	& ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)
	728	END DO
	729	END DO
	730	!
	731	END DO ! jl2
	732	!
	733	END DO ! jl1
	734	!
[10994]	735	! roundoff errors
	736	!----------------
[8586]	737	! In case ridging/rafting lead to very small negative values (sometimes it happens)
[10994]	738	CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d )
[8586]	739	!
	740	END SUBROUTINE rdgrft_shift
	741
	742
	743	SUBROUTINE ice_strength
	744	!!----------------------------------------------------------------------
	745	!! * ROUTINE ice_strength *
	746	!!
	747	!! ** Purpose : computes ice strength used in dynamics routines of ice thickness
	748	!!
	749	!! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2)
	750	!! dissipated per unit area removed from the ice pack under compression,
	751	!! and assumed proportional to the change in potential energy caused
	752	!! by ridging. Note that only Hibler's formulation is stable and that
	753	!! ice strength has to be smoothed
	754	!!----------------------------------------------------------------------
	755	INTEGER :: ji, jj, jl ! dummy loop indices
	756	INTEGER :: ismooth ! smoothing the resistance to deformation
	757	INTEGER :: itframe ! number of time steps for the P smoothing
	758	REAL(wp) :: zp, z1_3 ! local scalars
	759	REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here
	760	REAL(wp), DIMENSION(jpi,jpj) :: zstrp1, zstrp2 ! strength at previous time steps
	761	!!----------------------------------------------------------------------
	762	! !--------------------------------------------------!
	763	IF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method !
	764	! !--------------------------------------------------!
	765	strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
	766	ismooth = 1
	767	! !--------------------------------------------------!
	768	ELSE ! Zero strength !
	769	! !--------------------------------------------------!
	770	strength(:,:) = 0._wp
	771	ismooth = 0
	772	ENDIF
	773	! !--------------------------------------------------!
	774	SELECT CASE( ismooth ) ! Smoothing ice strength !
	775	! !--------------------------------------------------!
	776	CASE( 1 ) !--- Spatial smoothing
	777	DO jj = 2, jpjm1
	778	DO ji = 2, jpim1
	779	IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
	780	zworka(ji,jj) = ( 4.0 * strength(ji,jj) &
	781	& + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &
	782	& + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) &
	783	& ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) )
	784	ELSE
	785	zworka(ji,jj) = 0._wp
	786	ENDIF
	787	END DO
	788	END DO
	789
	790	DO jj = 2, jpjm1
	791	DO ji = 2, jpim1
	792	strength(ji,jj) = zworka(ji,jj)
	793	END DO
	794	END DO
[10425]	795	CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
[8586]	796	!
	797	CASE( 2 ) !--- Temporal smoothing
	798	IF ( kt_ice == nit000 ) THEN
	799	zstrp1(:,:) = 0._wp
	800	zstrp2(:,:) = 0._wp
	801	ENDIF
	802	!
	803	DO jj = 2, jpjm1
	804	DO ji = 2, jpim1
	805	IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN
	806	itframe = 1 ! number of time steps for the running mean
	807	IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1
	808	IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1
	809	zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe
	810	zstrp2 (ji,jj) = zstrp1 (ji,jj)
	811	zstrp1 (ji,jj) = strength(ji,jj)
	812	strength(ji,jj) = zp
	813	ENDIF
	814	END DO
	815	END DO
[10425]	816	CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. )
[8586]	817	!
	818	END SELECT
	819	!
	820	END SUBROUTINE ice_strength
	821
[9452]	822
	823	SUBROUTINE ice_dyn_1d2d( kn )
	824	!!-----------------------------------------------------------------------
	825	!! * ROUTINE ice_dyn_1d2d *
	826	!!
	827	!! ** Purpose : move arrays from 1d to 2d and the reverse
	828	!!-----------------------------------------------------------------------
	829	INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D ; 2= from 1D to 2D
	830	!
	831	INTEGER :: jl, jk ! dummy loop indices
	832	!!-----------------------------------------------------------------------
	833	!
	834	SELECT CASE( kn )
	835	! !---------------------!
	836	CASE( 1 ) !== from 2D to 1D ==!
	837	! !---------------------!
	838	! fields used but not modified
	839	CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
	840	CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
	841	! the following fields are modified in this routine
	842	!!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
	843	!!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
	844	!!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
	845	CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
	846	CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
	847	CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
	848	CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
	849	CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
	850	DO jl = 1, jpl
	851	DO jk = 1, nlay_s
	852	CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
	853	END DO
	854	DO jk = 1, nlay_i
	855	CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
	856	END DO
	857	END DO
	858	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) )
	859	CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) )
	860	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) )
	861	CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d (1:npti), hfx_dyn (:,:) )
	862	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
[10511]	863	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) )
[10994]	864	!
[9452]	865	! !---------------------!
	866	CASE( 2 ) !== from 1D to 2D ==!
	867	! !---------------------!
	868	CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
	869	CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
	870	CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
	871	CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
	872	CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
	873	CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
	874	CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
	875	CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
	876	DO jl = 1, jpl
	877	DO jk = 1, nlay_s
	878	CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
	879	END DO
	880	DO jk = 1, nlay_i
	881	CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
	882	END DO
	883	END DO
	884	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) )
	885	CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) )
	886	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) )
	887	CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d (1:npti), hfx_dyn (:,:) )
	888	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
[10511]	889	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) )
[9452]	890	!
	891	END SELECT
	892	!
	893	END SUBROUTINE ice_dyn_1d2d
	894
	895
[8586]	896	SUBROUTINE ice_dyn_rdgrft_init
	897	!!-------------------------------------------------------------------
	898	!! * ROUTINE ice_dyn_rdgrft_init *
	899	!!
	900	!! ** Purpose : Physical constants and parameters linked
	901	!! to the mechanical ice redistribution
	902	!!
	903	!! ** Method : Read the namdyn_rdgrft namelist
	904	!! and check the parameters values
	905	!! called at the first timestep (nit000)
	906	!!
	907	!! ** input : Namelist namdyn_rdgrft
	908	!!-------------------------------------------------------------------
	909	INTEGER :: ios ! Local integer output status for namelist read
	910	!!
	911	NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
	912	& rn_csrdg , &
	913	& ln_partf_lin, rn_gstar, &
	914	& ln_partf_exp, rn_astar, &
	915	& ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg, &
	916	& ln_rafting, rn_hraft, rn_craft , rn_fsnwrft, rn_fpndrft
	917	!!-------------------------------------------------------------------
	918	!
	919	REWIND( numnam_ice_ref ) ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
	920	READ ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901)
[11536]	921	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist' )
[8586]	922	REWIND( numnam_ice_cfg ) ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution
	923	READ ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 )
[11536]	924	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist' )
[8586]	925	IF(lwm) WRITE ( numoni, namdyn_rdgrft )
	926	!
	927	IF (lwp) THEN ! control print
	928	WRITE(numout,*)
	929	WRITE(numout,*) 'ice_dyn_rdgrft_init: ice parameters for ridging/rafting '
	930	WRITE(numout,*) '~~~~~~~~~~~~~~~~~~'
	931	WRITE(numout,*) ' Namelist namdyn_rdgrft:'
	932	WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_H79 = ', ln_str_H79
	933	WRITE(numout,*) ' 1st bulk-rheology parameter rn_pstar = ', rn_pstar
	934	WRITE(numout,*) ' 2nd bulk-rhelogy parameter rn_crhg = ', rn_crhg
	935	WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg
	936	WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin
	937	WRITE(numout,*) ' Fraction of ice coverage contributing to ridging rn_gstar = ', rn_gstar
	938	WRITE(numout,*) ' Exponential ridging participation function ln_partf_exp = ', ln_partf_exp
	939	WRITE(numout,) ' Equivalent to G for an exponential function rn_astar = ', rn_astar
	940	WRITE(numout,*) ' Ridging of ice sheets or not ln_ridging = ', ln_ridging
	941	WRITE(numout,*) ' max ridged ice thickness rn_hstar = ', rn_hstar
	942	WRITE(numout,*) ' Initial porosity of ridges rn_porordg = ', rn_porordg
	943	WRITE(numout,*) ' Fraction of snow volume conserved during ridging rn_fsnwrdg = ', rn_fsnwrdg
	944	WRITE(numout,*) ' Fraction of pond volume conserved during ridging rn_fpndrdg = ', rn_fpndrdg
	945	WRITE(numout,*) ' Rafting of ice sheets or not ln_rafting = ', ln_rafting
	946	WRITE(numout,*) ' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft
	947	WRITE(numout,*) ' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft
	948	WRITE(numout,*) ' Fraction of snow volume conserved during rafting rn_fsnwrft = ', rn_fsnwrft
	949	WRITE(numout,*) ' Fraction of pond volume conserved during rafting rn_fpndrft = ', rn_fpndrft
	950	ENDIF
	951	!
	952	IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
	953	CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
	954	ENDIF
[10945]	955	!
	956	IF( .NOT. ln_icethd ) THEN
	957	rn_porordg = 0._wp
	958	rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp
	959	rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp
	960	IF( lwp ) THEN
	961	WRITE(numout,*) ' ==> only ice dynamics is activated, thus some parameters must be changed'
	962	WRITE(numout,*) ' rn_porordg = ', rn_porordg
	963	WRITE(numout,*) ' rn_fsnwrdg = ', rn_fsnwrdg
	964	WRITE(numout,*) ' rn_fpndrdg = ', rn_fpndrdg
	965	WRITE(numout,*) ' rn_fsnwrft = ', rn_fsnwrft
	966	WRITE(numout,*) ' rn_fpndrft = ', rn_fpndrft
	967	ENDIF
	968	ENDIF
	969	! ! allocate arrays
[8586]	970	IF( ice_dyn_rdgrft_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' )
	971	!
	972	END SUBROUTINE ice_dyn_rdgrft_init
	973
	974	#else
	975	!!----------------------------------------------------------------------
[9570]	976	!! Default option Empty module NO SI3 sea-ice model
[8586]	977	!!----------------------------------------------------------------------
	978	#endif
	979
	980	!!======================================================================
	981	END MODULE icedyn_rdgrft

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

Download in other formats: