Context Navigation

icedyn_rdgrft.F90

Last change on this file was 8879, checked in by frrh, 6 years ago
Merge in http://fcm3/projects/NEMO.xm/log/branches/UKMO/dev_r8183_ICEMODEL_svn_removed revisions 8738:8847 inclusive.
File size: 50.0 KB

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

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

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

Context Navigation

source: branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90

Download in other formats: