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.

icevar.F90 in NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icevar.F90 @ 11352

Last change on this file since 11352 was 11352, checked in by smasson, 5 years ago
dev_r10984_HPC-13 : small bugfixes
Property svn:keywords set to `Id`
File size: 57.8 KB

Rev	Line
[8586]	1	MODULE icevar
	2	!!======================================================================
	3	!! * MODULE icevar *
[9604]	4	!! sea-ice: series of functions to transform or compute ice variables
	5	!!======================================================================
	6	!! History : - ! 2006-01 (M. Vancoppenolle) Original code
	7	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
	8	!!----------------------------------------------------------------------
	9	#if defined key_si3
	10	!!----------------------------------------------------------------------
	11	!! 'key_si3' SI3 sea-ice model
	12	!!----------------------------------------------------------------------
[8586]	13	!!
	14	!! There are three sets of variables
	15	!! VGLO : global variables of the model
	16	!! - v_i (jpi,jpj,jpl)
	17	!! - v_s (jpi,jpj,jpl)
	18	!! - a_i (jpi,jpj,jpl)
	19	!! - t_s (jpi,jpj,jpl)
	20	!! - e_i (jpi,jpj,nlay_i,jpl)
[9604]	21	!! - e_s (jpi,jpj,nlay_s,jpl)
[8586]	22	!! - sv_i(jpi,jpj,jpl)
	23	!! - oa_i(jpi,jpj,jpl)
	24	!! VEQV : equivalent variables sometimes used in the model
	25	!! - h_i(jpi,jpj,jpl)
	26	!! - h_s(jpi,jpj,jpl)
	27	!! - t_i(jpi,jpj,nlay_i,jpl)
	28	!! ...
	29	!! VAGG : aggregate variables, averaged/summed over all
	30	!! thickness categories
	31	!! - vt_i(jpi,jpj)
	32	!! - vt_s(jpi,jpj)
	33	!! - at_i(jpi,jpj)
[9433]	34	!! - et_s(jpi,jpj) total snow heat content
	35	!! - et_i(jpi,jpj) total ice thermal content
	36	!! - sm_i(jpi,jpj) mean ice salinity
	37	!! - tm_i(jpi,jpj) mean ice temperature
	38	!! - tm_s(jpi,jpj) mean snw temperature
[8586]	39	!!----------------------------------------------------------------------
	40	!! ice_var_agg : integrate variables over layers and categories
	41	!! ice_var_glo2eqv : transform from VGLO to VEQV
	42	!! ice_var_eqv2glo : transform from VEQV to VGLO
	43	!! ice_var_salprof : salinity profile in the ice
	44	!! ice_var_salprof1d : salinity profile in the ice 1D
	45	!! ice_var_zapsmall : remove very small area and volume
[10994]	46	!! ice_var_zapneg : remove negative ice fields
	47	!! ice_var_roundoff : remove negative values arising from roundoff erros
[8813]	48	!! ice_var_itd : convert 1-cat to jpl-cat
	49	!! ice_var_itd2 : convert N-cat to jpl-cat
[8586]	50	!! ice_var_bv : brine volume
[8984]	51	!! ice_var_enthalpy : compute ice and snow enthalpies from temperature
[10332]	52	!! ice_var_sshdyn : compute equivalent ssh in lead
[8586]	53	!!----------------------------------------------------------------------
	54	USE dom_oce ! ocean space and time domain
	55	USE phycst ! physical constants (ocean directory)
[10332]	56	USE sbc_oce , ONLY : sss_m, ln_ice_embd, nn_fsbc
[8586]	57	USE ice ! sea-ice: variables
	58	USE ice1D ! sea-ice: thermodynamics variables
	59	!
	60	USE in_out_manager ! I/O manager
	61	USE lib_mpp ! MPP library
	62	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	63
	64	IMPLICIT NONE
	65	PRIVATE
	66
	67	PUBLIC ice_var_agg
	68	PUBLIC ice_var_glo2eqv
	69	PUBLIC ice_var_eqv2glo
	70	PUBLIC ice_var_salprof
	71	PUBLIC ice_var_salprof1d
	72	PUBLIC ice_var_zapsmall
[9943]	73	PUBLIC ice_var_zapneg
[10994]	74	PUBLIC ice_var_roundoff
[8586]	75	PUBLIC ice_var_bv
[8984]	76	PUBLIC ice_var_enthalpy
[10332]	77	PUBLIC ice_var_sshdyn
[11223]	78	PUBLIC ice_var_itd
[8586]	79
[11223]	80	INTERFACE ice_var_itd
	81	MODULE PROCEDURE ice_var_itd_1c1c, ice_var_itd_Nc1c, ice_var_itd_1cMc, ice_var_itd_NcMc
	82	END INTERFACE
	83
[8586]	84	!!----------------------------------------------------------------------
[9598]	85	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	86	!! $Id$
[10068]	87	!! Software governed by the CeCILL license (see ./LICENSE)
[8586]	88	!!----------------------------------------------------------------------
	89	CONTAINS
	90
	91	SUBROUTINE ice_var_agg( kn )
	92	!!-------------------------------------------------------------------
	93	!! * ROUTINE ice_var_agg *
	94	!!
	95	!! ** Purpose : aggregates ice-thickness-category variables to
	96	!! all-ice variables, i.e. it turns VGLO into VAGG
	97	!!-------------------------------------------------------------------
	98	INTEGER, INTENT( in ) :: kn ! =1 state variables only
	99	! ! >1 state variables + others
	100	!
	101	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[9433]	102	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z1_at_i, z1_vt_i, z1_vt_s
[8586]	103	!!-------------------------------------------------------------------
	104	!
	105	! ! integrated values
	106	vt_i(:,:) = SUM( v_i(:,:,:) , dim=3 )
	107	vt_s(:,:) = SUM( v_s(:,:,:) , dim=3 )
	108	at_i(:,:) = SUM( a_i(:,:,:) , dim=3 )
	109	et_s(:,:) = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
	110	et_i(:,:) = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
[8813]	111	!
[8637]	112	at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
	113	vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
[8813]	114	!
	115	ato_i(:,:) = 1._wp - at_i(:,:) ! open water fraction
[8586]	116
[9725]	117	! The following fields are calculated for diagnostics and outputs only
	118	! ==> Do not use them for other purposes
[8586]	119	IF( kn > 1 ) THEN
	120	!
[9433]	121	ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) )
[8586]	122	WHERE( at_i(:,:) > epsi20 ) ; z1_at_i(:,:) = 1._wp / at_i(:,:)
	123	ELSEWHERE ; z1_at_i(:,:) = 0._wp
	124	END WHERE
	125	WHERE( vt_i(:,:) > epsi20 ) ; z1_vt_i(:,:) = 1._wp / vt_i(:,:)
	126	ELSEWHERE ; z1_vt_i(:,:) = 0._wp
	127	END WHERE
[9433]	128	WHERE( vt_s(:,:) > epsi20 ) ; z1_vt_s(:,:) = 1._wp / vt_s(:,:)
	129	ELSEWHERE ; z1_vt_s(:,:) = 0._wp
	130	END WHERE
[8586]	131	!
	132	! ! mean ice/snow thickness
	133	hm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
	134	hm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
	135	!
	136	! ! mean temperature (K), salinity and age
	137	tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
	138	tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
	139	om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
[9118]	140	sm_i (:,:) = SUM( sv_i(:,:,:) , dim=3 ) * z1_vt_i(:,:)
[8586]	141	!
	142	tm_i(:,:) = 0._wp
[9433]	143	tm_s(:,:) = 0._wp
[8586]	144	DO jl = 1, jpl
	145	DO jk = 1, nlay_i
	146	tm_i(:,:) = tm_i(:,:) + r1_nlay_i * t_i (:,:,jk,jl) * v_i(:,:,jl) * z1_vt_i(:,:)
	147	END DO
[9433]	148	DO jk = 1, nlay_s
	149	tm_s(:,:) = tm_s(:,:) + r1_nlay_s * t_s (:,:,jk,jl) * v_s(:,:,jl) * z1_vt_s(:,:)
	150	END DO
[8586]	151	END DO
	152	!
[8637]	153	! ! put rt0 where there is no ice
	154	WHERE( at_i(:,:)<=epsi20 )
	155	tm_su(:,:) = rt0
	156	tm_si(:,:) = rt0
	157	tm_i (:,:) = rt0
[9433]	158	tm_s (:,:) = rt0
[8637]	159	END WHERE
	160
[9433]	161	DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s )
[8586]	162	ENDIF
	163	!
	164	END SUBROUTINE ice_var_agg
	165
	166
	167	SUBROUTINE ice_var_glo2eqv
	168	!!-------------------------------------------------------------------
	169	!! * ROUTINE ice_var_glo2eqv *
	170	!!
	171	!! ** Purpose : computes equivalent variables as function of
	172	!! global variables, i.e. it turns VGLO into VEQV
	173	!!-------------------------------------------------------------------
	174	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	175	REAL(wp) :: ze_i ! local scalars
	176	REAL(wp) :: ze_s, ztmelts, zbbb, zccc ! - -
	177	REAL(wp) :: zhmax, z1_zhmax ! - -
	178	REAL(wp) :: zlay_i, zlay_s ! - -
	179	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i
	180	!!-------------------------------------------------------------------
	181
	182	!!gm Question 2: It is possible to define existence of sea-ice in a common way between
	183	!! ice area and ice volume ?
	184	!! the idea is to be able to define one for all at the begining of this routine
	185	!! a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
	186
[8637]	187	!---------------------------------------------------------------
	188	! Ice thickness, snow thickness, ice salinity, ice age and ponds
	189	!---------------------------------------------------------------
[8586]	190	! !--- inverse of the ice area
	191	WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
	192	ELSEWHERE ; z1_a_i(:,:,:) = 0._wp
	193	END WHERE
	194	!
	195	WHERE( v_i(:,:,:) > epsi20 ) ; z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
	196	ELSEWHERE ; z1_v_i(:,:,:) = 0._wp
	197	END WHERE
[8637]	198	! !--- ice thickness
	199	h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
[8586]	200
	201	zhmax = hi_max(jpl)
	202	z1_zhmax = 1._wp / hi_max(jpl)
[8637]	203	WHERE( h_i(:,:,jpl) > zhmax ) ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area
[8586]	204	h_i (:,:,jpl) = zhmax
	205	a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax
[8637]	206	z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)
[8586]	207	END WHERE
[8637]	208	! !--- snow thickness
	209	h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)
	210	! !--- ice age
	211	o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:)
	212	! !--- pond fraction and thickness
	213	a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:)
	214	WHERE( a_ip_frac(:,:,:) > epsi20 ) ; h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:)
	215	ELSEWHERE ; h_ip(:,:,:) = 0._wp
	216	END WHERE
	217	!
	218	! !--- salinity (with a minimum value imposed everywhere)
	219	IF( nn_icesal == 2 ) THEN
[8586]	220	WHERE( v_i(:,:,:) > epsi20 ) ; s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) )
	221	ELSEWHERE ; s_i(:,:,:) = rn_simin
	222	END WHERE
	223	ENDIF
[8637]	224	CALL ice_var_salprof ! salinity profile
[8586]	225
	226	!-------------------
	227	! Ice temperature [K] (with a minimum value (rt0 - 100.))
	228	!-------------------
	229	zlay_i = REAL( nlay_i , wp ) ! number of layers
	230	DO jl = 1, jpl
	231	DO jk = 1, nlay_i
	232	DO jj = 1, jpj
	233	DO ji = 1, jpi
	234	IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area
	235	!
[10930]	236	ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3]
[9935]	237	ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C]
[8586]	238	! Conversion q(S,T) -> T (second order equation)
[9935]	239	zbbb = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus
	240	zccc = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) )
	241	t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts
[8586]	242	!
[10930]	243	ELSE !--- no ice
[8586]	244	t_i(ji,jj,jk,jl) = rt0
	245	ENDIF
	246	END DO
	247	END DO
	248	END DO
	249	END DO
	250
	251	!--------------------
	252	! Snow temperature [K] (with a minimum value (rt0 - 100.))
	253	!--------------------
	254	zlay_s = REAL( nlay_s , wp )
	255	DO jk = 1, nlay_s
	256	WHERE( v_s(:,:,:) > epsi20 ) !--- icy area
[9121]	257	t_s(:,:,jk,:) = rt0 + MAX( -100._wp , &
[9935]	258	& MIN( r1_rcpi * ( -r1_rhos * ( e_s(:,:,jk,:) / v_s(:,:,:) * zlay_s ) + rLfus ) , 0._wp ) )
[8586]	259	ELSEWHERE !--- no ice
	260	t_s(:,:,jk,:) = rt0
	261	END WHERE
	262	END DO
[8813]	263	!
[8586]	264	! integrated values
	265	vt_i (:,:) = SUM( v_i, dim=3 )
	266	vt_s (:,:) = SUM( v_s, dim=3 )
	267	at_i (:,:) = SUM( a_i, dim=3 )
	268	!
	269	END SUBROUTINE ice_var_glo2eqv
	270
	271
	272	SUBROUTINE ice_var_eqv2glo
	273	!!-------------------------------------------------------------------
	274	!! * ROUTINE ice_var_eqv2glo *
	275	!!
	276	!! ** Purpose : computes global variables as function of
	277	!! equivalent variables, i.e. it turns VEQV into VGLO
	278	!!-------------------------------------------------------------------
	279	!
[8637]	280	v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:)
	281	v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:)
	282	sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
	283	v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
[8586]	284	!
	285	END SUBROUTINE ice_var_eqv2glo
	286
	287
	288	SUBROUTINE ice_var_salprof
	289	!!-------------------------------------------------------------------
	290	!! * ROUTINE ice_var_salprof *
	291	!!
	292	!! ** Purpose : computes salinity profile in function of bulk salinity
	293	!!
	294	!! ** Method : If bulk salinity greater than zsi1,
	295	!! the profile is assumed to be constant (S_inf)
	296	!! If bulk salinity lower than zsi0,
	297	!! the profile is linear with 0 at the surface (S_zero)
	298	!! If it is between zsi0 and zsi1, it is a
	299	!! alpha-weighted linear combination of s_inf and s_zero
	300	!!
	301	!! ** References : Vancoppenolle et al., 2007
	302	!!-------------------------------------------------------------------
	303	INTEGER :: ji, jj, jk, jl ! dummy loop index
	304	REAL(wp) :: zsal, z1_dS
	305	REAL(wp) :: zargtemp , zs0, zs
	306	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_slope_s, zalpha ! case 2 only
	307	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
	308	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
	309	!!-------------------------------------------------------------------
	310
	311	!!gm Question: Remove the option 3 ? How many years since it last use ?
	312
	313	SELECT CASE ( nn_icesal )
	314	!
	315	! !---------------------------------------!
	316	CASE( 1 ) ! constant salinity in time and space !
	317	! !---------------------------------------!
	318	sz_i(:,:,:,:) = rn_icesal
[9118]	319	s_i (:,:,:) = rn_icesal
[8586]	320	!
	321	! !---------------------------------------------!
	322	CASE( 2 ) ! time varying salinity with linear profile !
	323	! !---------------------------------------------!
	324	!
	325	ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
	326	!
	327	DO jl = 1, jpl
	328	DO jk = 1, nlay_i
	329	sz_i(:,:,jk,jl) = s_i(:,:,jl)
	330	END DO
	331	END DO
	332	! ! Slope of the linear profile
	333	WHERE( h_i(:,:,:) > epsi20 ) ; z_slope_s(:,:,:) = 2._wp * s_i(:,:,:) / h_i(:,:,:)
[9118]	334	ELSEWHERE ; z_slope_s(:,:,:) = 0._wp
[8586]	335	END WHERE
	336	!
	337	z1_dS = 1._wp / ( zsi1 - zsi0 )
	338	DO jl = 1, jpl
	339	DO jj = 1, jpj
	340	DO ji = 1, jpi
	341	zalpha(ji,jj,jl) = MAX( 0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp ) )
	342	! ! force a constant profile when SSS too low (Baltic Sea)
	343	IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) ) zalpha(ji,jj,jl) = 0._wp
	344	END DO
	345	END DO
	346	END DO
[8813]	347	!
[8586]	348	! Computation of the profile
	349	DO jl = 1, jpl
	350	DO jk = 1, nlay_i
	351	DO jj = 1, jpj
	352	DO ji = 1, jpi
	353	! ! linear profile with 0 surface value
	354	zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
	355	zs = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl) ! weighting the profile
	356	sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
	357	END DO
	358	END DO
	359	END DO
	360	END DO
	361	!
	362	DEALLOCATE( z_slope_s , zalpha )
	363	!
	364	! !-------------------------------------------!
	365	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
	366	! !-------------------------------------------! (mean = 2.30)
	367	!
	368	s_i(:,:,:) = 2.30_wp
	369	!!gm Remark: if we keep the case 3, then compute an store one for all time-step
	370	!! a array S_prof(1:nlay_i) containing the calculation and just do:
	371	! DO jk = 1, nlay_i
	372	! sz_i(:,:,jk,:) = S_prof(jk)
	373	! END DO
	374	!!gm end
	375	!
	376	DO jl = 1, jpl
	377	DO jk = 1, nlay_i
	378	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
	379	sz_i(:,:,jk,jl) = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) )
	380	END DO
	381	END DO
	382	!
	383	END SELECT
	384	!
	385	END SUBROUTINE ice_var_salprof
	386
[8813]	387
[8586]	388	SUBROUTINE ice_var_salprof1d
	389	!!-------------------------------------------------------------------
	390	!! * ROUTINE ice_var_salprof1d *
	391	!!
	392	!! ** Purpose : 1d computation of the sea ice salinity profile
	393	!! Works with 1d vectors and is used by thermodynamic modules
	394	!!-------------------------------------------------------------------
	395	INTEGER :: ji, jk ! dummy loop indices
	396	REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars
	397	REAL(wp) :: zs, zs0 ! - -
	398	!
	399	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s, zalpha !
	400	REAL(wp), PARAMETER :: zsi0 = 3.5_wp
	401	REAL(wp), PARAMETER :: zsi1 = 4.5_wp
	402	!!-------------------------------------------------------------------
	403	!
	404	SELECT CASE ( nn_icesal )
	405	!
	406	! !---------------------------------------!
	407	CASE( 1 ) ! constant salinity in time and space !
	408	! !---------------------------------------!
	409	sz_i_1d(1:npti,:) = rn_icesal
	410	!
	411	! !---------------------------------------------!
	412	CASE( 2 ) ! time varying salinity with linear profile !
	413	! !---------------------------------------------!
	414	!
	415	ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
	416	!
	417	! ! Slope of the linear profile
	418	WHERE( h_i_1d(1:npti) > epsi20 ) ; z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
[8984]	419	ELSEWHERE ; z_slope_s(1:npti) = 0._wp
[8586]	420	END WHERE
	421
	422	z1_dS = 1._wp / ( zsi1 - zsi0 )
	423	DO ji = 1, npti
	424	zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp ) )
	425	! ! force a constant profile when SSS too low (Baltic Sea)
	426	IF( 2._wp * s_i_1d(ji) >= sss_1d(ji) ) zalpha(ji) = 0._wp
	427	END DO
	428	!
	429	! Computation of the profile
	430	DO jk = 1, nlay_i
	431	DO ji = 1, npti
	432	! ! linear profile with 0 surface value
	433	zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
	434	zs = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * s_i_1d(ji)
	435	sz_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
	436	END DO
	437	END DO
	438	!
	439	DEALLOCATE( z_slope_s, zalpha )
	440
	441	! !-------------------------------------------!
	442	CASE( 3 ) ! constant salinity with a fix profile ! (Schwarzacher (1959) multiyear salinity profile
	443	! !-------------------------------------------! (mean = 2.30)
	444	!
	445	s_i_1d(1:npti) = 2.30_wp
	446	!
	447	!!gm cf remark in ice_var_salprof routine, CASE( 3 )
	448	DO jk = 1, nlay_i
	449	zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
	450	zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
	451	DO ji = 1, npti
	452	sz_i_1d(ji,jk) = zsal
	453	END DO
	454	END DO
	455	!
	456	END SELECT
	457	!
	458	END SUBROUTINE ice_var_salprof1d
	459
	460
	461	SUBROUTINE ice_var_zapsmall
	462	!!-------------------------------------------------------------------
	463	!! * ROUTINE ice_var_zapsmall *
	464	!!
	465	!! ** Purpose : Remove too small sea ice areas and correct fluxes
	466	!!-------------------------------------------------------------------
	467	INTEGER :: ji, jj, jl, jk ! dummy loop indices
	468	REAL(wp), DIMENSION(jpi,jpj) :: zswitch
	469	!!-------------------------------------------------------------------
	470	!
	471	DO jl = 1, jpl !== loop over the categories ==!
	472	!
[9448]	473	WHERE( a_i(:,:,jl) > epsi10 ) ; h_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
[8586]	474	ELSEWHERE ; h_i(:,:,jl) = 0._wp
	475	END WHERE
	476	!
	477	WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. h_i(:,:,jl) < epsi10 ) ; zswitch(:,:) = 0._wp
	478	ELSEWHERE ; zswitch(:,:) = 1._wp
	479	END WHERE
[8813]	480	!
[9943]	481	!-----------------------------------------------------------------
	482	! Zap ice energy and use ocean heat to melt ice
	483	!-----------------------------------------------------------------
[8586]	484	DO jk = 1, nlay_i
	485	DO jj = 1 , jpj
	486	DO ji = 1 , jpi
	487	! update exchanges with ocean
	488	hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
	489	e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
	490	t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
	491	END DO
	492	END DO
	493	END DO
[8813]	494	!
[9271]	495	DO jk = 1, nlay_s
	496	DO jj = 1 , jpj
	497	DO ji = 1 , jpi
	498	! update exchanges with ocean
	499	hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
	500	e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
	501	t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
	502	END DO
	503	END DO
	504	END DO
	505	!
[9943]	506	!-----------------------------------------------------------------
	507	! zap ice and snow volume, add water and salt to ocean
	508	!-----------------------------------------------------------------
[8586]	509	DO jj = 1 , jpj
	510	DO ji = 1 , jpi
	511	! update exchanges with ocean
[9935]	512	sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_rdtice
	513	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_rdtice
	514	wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_rdtice
[8813]	515	!
[8586]	516	a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
	517	v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
	518	v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
	519	t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
	520	oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
	521	sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
[8813]	522	!
[8586]	523	h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
	524	h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
[8813]	525	!
[8637]	526	a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
	527	v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
[8813]	528	!
[8586]	529	END DO
	530	END DO
[8813]	531	!
[8586]	532	END DO
	533
	534	! to be sure that at_i is the sum of a_i(jl)
	535	at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
	536	vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
	537
	538	! open water = 1 if at_i=0
	539	WHERE( at_i(:,:) == 0._wp ) ato_i(:,:) = 1._wp
	540	!
	541	END SUBROUTINE ice_var_zapsmall
	542
	543
[10930]	544	SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
[9943]	545	!!-------------------------------------------------------------------
	546	!! * ROUTINE ice_var_zapneg *
	547	!!
	548	!! ** Purpose : Remove negative sea ice fields and correct fluxes
	549	!!-------------------------------------------------------------------
[10930]	550	REAL(wp) , INTENT(in ) :: pdt ! tracer time-step
[9943]	551	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area
	552	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume
	553	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_s ! snw volume
	554	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: psv_i ! salt content
	555	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: poa_i ! age content
	556	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_i ! ice concentration
	557	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
	558	REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_ip ! melt pond volume
	559	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content
	560	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content
[10930]	561	!
	562	INTEGER :: ji, jj, jl, jk ! dummy loop indices
	563	REAL(wp) :: z1_dt
[9943]	564	!!-------------------------------------------------------------------
	565	!
[10930]	566	z1_dt = 1._wp / pdt
[9943]	567	!
	568	DO jl = 1, jpl !== loop over the categories ==!
	569	!
[10945]	570	! make sure a_i=0 where v_i<=0
	571	WHERE( pv_i(:,:,:) <= 0._wp ) pa_i(:,:,:) = 0._wp
	572
[9943]	573	!----------------------------------------
	574	! zap ice energy and send it to the ocean
	575	!----------------------------------------
	576	DO jk = 1, nlay_i
	577	DO jj = 1 , jpj
	578	DO ji = 1 , jpi
[10945]	579	IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
[10930]	580	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
[9943]	581	pe_i(ji,jj,jk,jl) = 0._wp
	582	ENDIF
	583	END DO
	584	END DO
	585	END DO
	586	!
	587	DO jk = 1, nlay_s
	588	DO jj = 1 , jpj
	589	DO ji = 1 , jpi
[10945]	590	IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
[10930]	591	hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
[9943]	592	pe_s(ji,jj,jk,jl) = 0._wp
	593	ENDIF
	594	END DO
	595	END DO
	596	END DO
	597	!
	598	!-----------------------------------------------------
	599	! zap ice and snow volume, add water and salt to ocean
	600	!-----------------------------------------------------
	601	DO jj = 1 , jpj
	602	DO ji = 1 , jpi
[10945]	603	IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
[10930]	604	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
[9943]	605	pv_i (ji,jj,jl) = 0._wp
	606	ENDIF
[10945]	607	IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
[10930]	608	wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
[9943]	609	pv_s (ji,jj,jl) = 0._wp
	610	ENDIF
[10945]	611	IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
[10930]	612	sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
[9943]	613	psv_i (ji,jj,jl) = 0._wp
	614	ENDIF
	615	END DO
	616	END DO
	617	!
	618	END DO
	619	!
[10413]	620	WHERE( pato_i(:,:) < 0._wp ) pato_i(:,:) = 0._wp
	621	WHERE( poa_i (:,:,:) < 0._wp ) poa_i (:,:,:) = 0._wp
	622	WHERE( pa_i (:,:,:) < 0._wp ) pa_i (:,:,:) = 0._wp
	623	WHERE( pa_ip (:,:,:) < 0._wp ) pa_ip (:,:,:) = 0._wp
	624	WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
	625	! but it does not change conservation, so keep it this way is ok
	626	!
[9943]	627	END SUBROUTINE ice_var_zapneg
	628
[10994]	629
	630	SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i )
	631	!!-------------------------------------------------------------------
	632	!! * ROUTINE ice_var_roundoff *
	633	!!
	634	!! ** Purpose : Remove negative sea ice values arising from roundoff errors
	635	!!-------------------------------------------------------------------
	636	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_i ! ice concentration
	637	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_i ! ice volume
	638	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_s ! snw volume
	639	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: psv_i ! salt content
	640	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: poa_i ! age content
	641	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction
	642	REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume
	643	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content
	644	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content
	645	!!-------------------------------------------------------------------
	646	!
	647	WHERE( pa_i (1:npti,:) < 0._wp .AND. pa_i (1:npti,:) > -epsi10 ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0
	648	WHERE( pv_i (1:npti,:) < 0._wp .AND. pv_i (1:npti,:) > -epsi10 ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0
	649	WHERE( pv_s (1:npti,:) < 0._wp .AND. pv_s (1:npti,:) > -epsi10 ) pv_s (1:npti,:) = 0._wp ! v_s must be >= 0
	650	WHERE( psv_i(1:npti,:) < 0._wp .AND. psv_i(1:npti,:) > -epsi10 ) psv_i(1:npti,:) = 0._wp ! sv_i must be >= 0
	651	WHERE( poa_i(1:npti,:) < 0._wp .AND. poa_i(1:npti,:) > -epsi10 ) poa_i(1:npti,:) = 0._wp ! oa_i must be >= 0
	652	WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0
	653	WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0
	654	IF ( ln_pnd_H12 ) THEN
	655	WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0
	656	WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0
	657	ENDIF
	658	!
	659	END SUBROUTINE ice_var_roundoff
[9943]	660
[11223]	661
	662	SUBROUTINE ice_var_bv
	663	!!-------------------------------------------------------------------
	664	!! * ROUTINE ice_var_bv *
	665	!!
	666	!! ** Purpose : computes mean brine volume (%) in sea ice
	667	!!
	668	!! ** Method : e = - 0.054 * S (ppt) / T (C)
	669	!!
	670	!! References : Vancoppenolle et al., JGR, 2007
	671	!!-------------------------------------------------------------------
	672	INTEGER :: ji, jj, jk, jl ! dummy loop indices
	673	!!-------------------------------------------------------------------
	674	!
	675	!!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done
	676	!! instead of setting everything to zero as just below
	677	bv_i (:,:,:) = 0._wp
	678	DO jl = 1, jpl
	679	DO jk = 1, nlay_i
	680	WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )
	681	bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
	682	END WHERE
	683	END DO
	684	END DO
	685	WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
	686	ELSEWHERE ; bvm_i(:,:) = 0._wp
	687	END WHERE
	688	!
	689	END SUBROUTINE ice_var_bv
	690
	691
	692	SUBROUTINE ice_var_enthalpy
	693	!!-------------------------------------------------------------------
	694	!! * ROUTINE ice_var_enthalpy *
	695	!!
	696	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
	697	!!
	698	!! ** Method : Formula (Bitz and Lipscomb, 1999)
	699	!!-------------------------------------------------------------------
	700	INTEGER :: ji, jk ! dummy loop indices
	701	REAL(wp) :: ztmelts ! local scalar
	702	!!-------------------------------------------------------------------
	703	!
	704	DO jk = 1, nlay_i ! Sea ice energy of melting
	705	DO ji = 1, npti
	706	ztmelts = - rTmlt * sz_i_1d(ji,jk)
	707	t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue
	708	! (sometimes zdf scheme produces abnormally high temperatures)
	709	e_i_1d(ji,jk) = rhoi * ( rcpi * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) ) &
	710	& + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) ) &
	711	& - rcp * ztmelts )
	712	END DO
	713	END DO
	714	DO jk = 1, nlay_s ! Snow energy of melting
	715	DO ji = 1, npti
	716	e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
	717	END DO
	718	END DO
	719	!
	720	END SUBROUTINE ice_var_enthalpy
	721
[10994]	722
[11223]	723	FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
	724	!!---------------------------------------------------------------------
	725	!! * ROUTINE ice_var_sshdyn *
	726	!!
	727	!! ** Purpose : compute the equivalent ssh in lead when sea ice is embedded
	728	!!
	729	!! ** Method : ssh_lead = ssh + (Mice + Msnow) / rau0
	730	!!
	731	!! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
	732	!! Sea ice-ocean coupling using a rescaled vertical coordinate z*,
	733	!! Ocean Modelling, Volume 24, Issues 1-2, 2008
	734	!!----------------------------------------------------------------------
	735	!
	736	! input
	737	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh !: ssh [m]
	738	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass !: mass of snow and ice at current ice time step [Kg/m2]
	739	REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b !: mass of snow and ice at previous ice time step [Kg/m2]
	740	!
	741	! output
	742	REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn ! equivalent ssh in lead [m]
	743	!
	744	! temporary
	745	REAL(wp) :: zintn, zintb ! time interpolation weights []
	746	REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload ! snow and ice load [m]
	747	!
	748	! compute ice load used to define the equivalent ssh in lead
	749	IF( ln_ice_embd ) THEN
	750	!
	751	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
	752	! = (1/nn_fsbc)^2 * {SUM[n] , n=0,nn_fsbc-1}
	753	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
	754	!
	755	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
	756	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
	757	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
	758	!
	759	zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
	760	!
	761	ELSE
	762	zsnwiceload(:,:) = 0.0_wp
	763	ENDIF
	764	! compute equivalent ssh in lead
	765	ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
	766	!
	767	END FUNCTION ice_var_sshdyn
	768
	769
	770	!!-------------------------------------------------------------------
	771	!! * INTERFACE ice_var_itd *
	772	!!
	773	!! ** Purpose : converting N-cat ice to jpl ice categories
	774	!!-------------------------------------------------------------------
[11332]	775	SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, &
	776	& ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i )
[8586]	777	!!-------------------------------------------------------------------
[11223]	778	!! ** Purpose : converting 1-cat ice to 1 ice category
	779	!!-------------------------------------------------------------------
[11332]	780	REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
	781	REAL(wp), DIMENSION(:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
	782	REAL(wp), DIMENSION(:), INTENT(in) , OPTIONAL :: ptmi, ptms, ptmsu, psmi ! input ice/snow temp & sal
	783	REAL(wp), DIMENSION(:), INTENT(inout), OPTIONAL :: pt_i, pt_s, pt_su, ps_i ! output ice/snow temp & sal
[11223]	784	!!-------------------------------------------------------------------
[11332]	785	! == thickness and concentration == !
	786	ph_i(:) = phti(:)
	787	ph_s(:) = phts(:)
	788	pa_i(:) = pati(:)
	789	!
	790	! == temperature and salinity == !
	791	IF( PRESENT( pt_i ) ) pt_i (:) = ptmi (:)
	792	IF( PRESENT( pt_s ) ) pt_s (:) = ptms (:)
	793	IF( PRESENT( pt_su ) ) pt_su(:) = ptmsu(:)
	794	IF( PRESENT( ps_i ) ) ps_i (:) = psmi (:)
	795
[11223]	796	END SUBROUTINE ice_var_itd_1c1c
	797
[11332]	798	SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, &
	799	& ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i )
[11223]	800	!!-------------------------------------------------------------------
	801	!! ** Purpose : converting N-cat ice to 1 ice category
	802	!!-------------------------------------------------------------------
[11332]	803	REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
	804	REAL(wp), DIMENSION(:) , INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
	805	REAL(wp), DIMENSION(:,:), INTENT(in) , OPTIONAL :: ptmi, ptms, ptmsu, psmi ! input ice/snow temp & sal
	806	REAL(wp), DIMENSION(:) , INTENT(inout), OPTIONAL :: pt_i, pt_s, pt_su, ps_i ! output ice/snow temp & sal
	807	!
	808	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs
	809	!
	810	INTEGER :: idim
[11223]	811	!!-------------------------------------------------------------------
	812	!
[11332]	813	idim = SIZE( phti, 1 )
[11223]	814	!
[11332]	815	! == thickness and concentration == !
	816	ALLOCATE( z1_ai(idim) )
	817	!
	818	pa_i(:) = SUM( pati(:,:), dim=2 )
	819
	820	WHERE( ( pa_i(:) ) /= 0._wp ) ; z1_ai(:) = 1._wp / pa_i(:)
	821	ELSEWHERE ; z1_ai(:) = 0._wp
[11223]	822	END WHERE
[11332]	823
	824	ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
	825	ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:)
[11223]	826	!
[11332]	827	! == temperature and salinity == !
	828	IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN
	829	!
	830	ALLOCATE( z1_vi(idim), z1_vs(idim) )
	831	!
	832	WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp ) ; z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) )
	833	ELSEWHERE ; z1_vi(:) = 0._wp
	834	END WHERE
	835	WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp ) ; z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) )
	836	ELSEWHERE ; z1_vs(:) = 0._wp
	837	END WHERE
	838	!
	839	IF( PRESENT( pt_i ) ) pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
	840	IF( PRESENT( pt_s ) ) pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
	841	IF( PRESENT( pt_su ) ) pt_su(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:)
	842	IF( PRESENT( ps_i ) ) ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
	843	!
	844	DEALLOCATE( z1_vi, z1_vs )
	845	!
	846	ENDIF
	847	!
	848	DEALLOCATE( z1_ai )
	849	!
[11223]	850	END SUBROUTINE ice_var_itd_Nc1c
	851
[11332]	852	SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, &
	853	& ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i )
[11223]	854	!!-------------------------------------------------------------------
[8586]	855	!!
[11223]	856	!! ** Purpose : converting 1-cat ice to jpl ice categories
[8586]	857	!!
	858	!! ice thickness distribution follows a gaussian law
	859	!! around the concentration of the most likely ice thickness
	860	!! (similar as iceistate.F90)
	861	!!
	862	!! ** Method: Iterative procedure
	863	!!
	864	!! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
	865	!!
	866	!! 2) Check whether the distribution conserves area and volume, positivity and
	867	!! category boundaries
	868	!!
	869	!! 3) If not (input ice is too thin), the last category is empty and
	870	!! the number of categories is reduced (jpl-1)
	871	!!
	872	!! 4) Iterate until ok (SUM(itest(:) = 4)
	873	!!
[11332]	874	!! ** Arguments : phti: 1-cat ice thickness
	875	!! phts: 1-cat snow depth
	876	!! pati: 1-cat ice concentration
[8586]	877	!!
	878	!! ** Output : jpl-cat
	879	!!
	880	!! (Example of application: BDY forcings when input are cell averaged)
	881	!!-------------------------------------------------------------------
[11332]	882	REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
	883	REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
	884	REAL(wp), DIMENSION(:) , INTENT(in) , OPTIONAL :: ptmi, ptms, ptmsu, psmi ! input ice/snow temp & sal
	885	REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL :: pt_i, pt_s, pt_su, ps_i ! output ice/snow temp & sal
	886	!
	887	INTEGER , DIMENSION(4) :: itest
	888	INTEGER :: ji, jk, jl
	889	INTEGER :: idim, i_fill, jl0
	890	REAL(wp) :: zarg, zV, zconv, zdh, zdv
[8586]	891	!!-------------------------------------------------------------------
	892	!
[11332]	893	! == thickness and concentration == !
[8586]	894	! distribution over the jpl ice categories
[11332]	895	! a gaussian distribution for ice concentration is used
	896	! then we check whether the distribution fullfills
	897	! volume and area conservation, positivity and ice categories bounds
	898	idim = SIZE( phti , 1 )
[8813]	899	!
[11332]	900	ph_i(1:idim,1:jpl) = 0._wp
	901	ph_s(1:idim,1:jpl) = 0._wp
	902	pa_i(1:idim,1:jpl) = 0._wp
[11223]	903	!
[8813]	904	DO ji = 1, idim
	905	!
[11332]	906	IF( phti(ji) > 0._wp ) THEN
[8813]	907	!
[8586]	908	! find which category (jl0) the input ice thickness falls into
	909	jl0 = jpl
	910	DO jl = 1, jpl
[11332]	911	IF ( ( phti(ji) >= hi_max(jl-1) ) .AND. ( phti(ji) < hi_max(jl) ) ) THEN
[8586]	912	jl0 = jl
	913	CYCLE
	914	ENDIF
	915	END DO
[8813]	916	!
[8586]	917	itest(:) = 0
	918	i_fill = jpl + 1 !------------------------------------
	919	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
	920	! !------------------------------------
	921	i_fill = i_fill - 1
	922	!
[11332]	923	ph_i(ji,1:jpl) = 0._wp
	924	pa_i(ji,1:jpl) = 0._wp
[8813]	925	itest(:) = 0
	926	!
[8586]	927	IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1
[11332]	928	ph_i(ji,1) = phti(ji)
	929	pa_i (ji,1) = pati (ji)
[8586]	930	ELSE !-- case ice is thicker: fill categories >1
	931	! thickness
	932	DO jl = 1, i_fill - 1
[11332]	933	ph_i(ji,jl) = hi_mean(jl)
[8586]	934	END DO
[8813]	935	!
[8586]	936	! concentration
[11332]	937	pa_i(ji,jl0) = pati(ji) / SQRT(REAL(jpl))
[8586]	938	DO jl = 1, i_fill - 1
	939	IF ( jl /= jl0 ) THEN
[11332]	940	zarg = ( ph_i(ji,jl) - phti(ji) ) / ( phti(ji) * 0.5_wp )
	941	pa_i(ji,jl) = pa_i (ji,jl0) * EXP(-zarg**2)
[8586]	942	ENDIF
	943	END DO
[8813]	944	!
[8586]	945	! last category
[11332]	946	pa_i(ji,i_fill) = pati(ji) - SUM( pa_i(ji,1:i_fill-1) )
	947	zV = SUM( pa_i(ji,1:i_fill-1) * ph_i(ji,1:i_fill-1) )
	948	ph_i(ji,i_fill) = ( phti(ji) * pati(ji) - zV ) / MAX( pa_i(ji,i_fill), epsi10 )
[8813]	949	!
[8885]	950	! correction if concentration of upper cat is greater than lower cat
	951	! (it should be a gaussian around jl0 but sometimes it is not)
[8586]	952	IF ( jl0 /= jpl ) THEN
	953	DO jl = jpl, jl0+1, -1
[11332]	954	IF ( pa_i(ji,jl) > pa_i(ji,jl-1) ) THEN
	955	zdv = ph_i(ji,jl) * pa_i(ji,jl)
	956	ph_i(ji,jl ) = 0._wp
	957	pa_i (ji,jl ) = 0._wp
	958	pa_i (ji,1:jl-1) = pa_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * phti(ji), epsi10 )
[8586]	959	END IF
[8813]	960	END DO
[8586]	961	ENDIF
[8813]	962	!
[8586]	963	ENDIF
[8813]	964	!
[8586]	965	! Compatibility tests
[11332]	966	zconv = ABS( pati(ji) - SUM( pa_i(ji,1:jpl) ) )
[8813]	967	IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation
	968	!
[11332]	969	zconv = ABS( phti(ji)pati(ji) - SUM( pa_i(ji,1:jpl)ph_i(ji,1:jpl) ) )
[8813]	970	IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation
	971	!
[11332]	972	IF ( ph_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?
[8813]	973	!
[8586]	974	itest(4) = 1
	975	DO jl = 1, i_fill
[11332]	976	IF ( pa_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations
[8586]	977	END DO
	978	! !----------------------------
	979	END DO ! end iteration on categories
[8813]	980	! !----------------------------
[8586]	981	ENDIF
	982	END DO
	983
[11332]	984	! Add Snow in each category where pa_i is not 0
[8586]	985	DO jl = 1, jpl
[8813]	986	DO ji = 1, idim
[11332]	987	IF( pa_i(ji,jl) > 0._wp ) THEN
	988	ph_s(ji,jl) = ph_i(ji,jl) * ( phts(ji) / phti(ji) )
[8586]	989	! In case snow load is in excess that would lead to transformation from snow to ice
	990	! Then, transfer the snow excess into the ice (different from icethd_dh)
[11332]	991	zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 )
[8586]	992	! recompute h_i, h_s avoiding out of bounds values
[11332]	993	ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh )
	994	ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos )
[8586]	995	ENDIF
	996	END DO
	997	END DO
	998	!
[11332]	999	! == temperature and salinity == !
	1000	IF( PRESENT( pt_i ) ) THEN
	1001	DO jl = 1, jpl
	1002	pt_i(:,jl) = ptmi (:)
	1003	END DO
	1004	ENDIF
	1005	IF( PRESENT( pt_s ) ) THEN
	1006	DO jl = 1, jpl
	1007	pt_s (:,jl) = ptms (:)
	1008	END DO
	1009	ENDIF
	1010	IF( PRESENT( pt_su ) ) THEN
	1011	DO jl = 1, jpl
	1012	pt_su(:,jl) = ptmsu(:)
	1013	END DO
	1014	ENDIF
	1015	IF( PRESENT( ps_i ) ) THEN
	1016	DO jl = 1, jpl
	1017	ps_i (:,jl) = psmi (:)
	1018	END DO
	1019	ENDIF
	1020	!
[11223]	1021	END SUBROUTINE ice_var_itd_1cMc
[8586]	1022
[11332]	1023	SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, &
	1024	& ptmi, ptms, ptmsu, psmi, pt_i, pt_s, pt_su, ps_i )
[8586]	1025	!!-------------------------------------------------------------------
[8813]	1026	!!
	1027	!! ** Purpose : converting N-cat ice to jpl ice categories
	1028	!!
	1029	!! ice thickness distribution follows a gaussian law
	1030	!! around the concentration of the most likely ice thickness
	1031	!! (similar as iceistate.F90)
	1032	!!
	1033	!! ** Method: Iterative procedure
	1034	!!
	1035	!! 1) Fill ice cat that correspond to input thicknesses
	1036	!! Find the lowest(jlmin) and highest(jlmax) cat that are filled
	1037	!!
	1038	!! 2) Expand the filling to the cat jlmin-1 and jlmax+1
[10994]	1039	!! by removing 25% ice area from jlmin and jlmax (resp.)
[8813]	1040	!!
	1041	!! 3) Expand the filling to the empty cat between jlmin and jlmax
	1042	!! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
	1043	!! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
	1044	!!
[11332]	1045	!! ** Arguments : phti: N-cat ice thickness
	1046	!! phts: N-cat snow depth
	1047	!! pati: N-cat ice concentration
[8813]	1048	!!
	1049	!! ** Output : jpl-cat
	1050	!!
	1051	!! (Example of application: BDY forcings when inputs have N-cat /= jpl)
	1052	!!-------------------------------------------------------------------
[11332]	1053	REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables
	1054	REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables
	1055	REAL(wp), DIMENSION(:,:), INTENT(in) , OPTIONAL :: ptmi, ptms, ptmsu, psmi ! input ice/snow temp & sal
	1056	REAL(wp), DIMENSION(:,:), INTENT(inout), OPTIONAL :: pt_i, pt_s, pt_su, ps_i ! output ice/snow temp & sal
	1057	!
	1058	INTEGER , ALLOCATABLE, DIMENSION(:,:) :: jlfil, jlfil2
	1059	INTEGER , ALLOCATABLE, DIMENSION(:) :: jlmax, jlmin
	1060	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs, ztmp
	1061	!
	1062	REAL(wp), PARAMETER :: ztrans = 0.25_wp
	1063	INTEGER :: ji, jl, jl1, jl2
[8813]	1064	INTEGER :: idim, icat
	1065	!!-------------------------------------------------------------------
	1066	!
[11332]	1067	idim = SIZE( phti, 1 )
	1068	icat = SIZE( phti, 2 )
	1069	!
	1070	! == thickness and concentration == !
[11223]	1071	! ! ---------------------- !
	1072	IF( icat == jpl ) THEN ! input cat = output cat !
	1073	! ! ---------------------- !
[11332]	1074	ph_i(:,:) = phti(:,:)
	1075	ph_s(:,:) = phts(:,:)
	1076	pa_i(:,:) = pati(:,:)
	1077	!
	1078	! == temperature and salinity == !
	1079	IF( PRESENT( pt_i ) ) pt_i (:,:) = ptmi (:,:)
	1080	IF( PRESENT( pt_s ) ) pt_s (:,:) = ptms (:,:)
	1081	IF( PRESENT( pt_su ) ) pt_su(:,:) = ptmsu(:,:)
	1082	IF( PRESENT( ps_i ) ) ps_i (:,:) = psmi (:,:)
[11223]	1083	! ! ---------------------- !
[11332]	1084	ELSEIF( icat == 1 ) THEN ! input cat = 1 !
[11223]	1085	! ! ---------------------- !
[11352]	1086	CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), ph_i(:,:), ph_s(:,:), pa_i (:,:) )
	1087	!!$ CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), ph_i(:,:), ph_s(:,:), pa_i (:,:), &
	1088	!!$ & ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:) )
[11223]	1089	! ! ---------------------- !
[11332]	1090	ELSEIF( jpl == 1 ) THEN ! output cat = 1 !
[11223]	1091	! ! ---------------------- !
[11352]	1092	CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), ph_i(:,1), ph_s(:,1), pa_i (:,1) )
	1093	!!$ CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), ph_i(:,1), ph_s(:,1), pa_i (:,1), &
	1094	!!$ & ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1) )
[11223]	1095	! ! ----------------------- !
	1096	ELSE ! input cat /= output cat !
	1097	! ! ----------------------- !
	1098
	1099	ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays
	1100	ALLOCATE( jlmin(idim), jlmax(idim) )
[8813]	1101
[11223]	1102	! --- initialize output fields to 0 --- !
[11332]	1103	ph_i(1:idim,1:jpl) = 0._wp
	1104	ph_s(1:idim,1:jpl) = 0._wp
	1105	pa_i(1:idim,1:jpl) = 0._wp
[11223]	1106	!
	1107	! --- fill the categories --- !
	1108	! find where cat-input = cat-output and fill cat-output fields
	1109	jlmax(:) = 0
	1110	jlmin(:) = 999
	1111	jlfil(:,:) = 0
	1112	DO jl1 = 1, jpl
	1113	DO jl2 = 1, icat
	1114	DO ji = 1, idim
[11332]	1115	IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN
[11223]	1116	! fill the right category
[11332]	1117	ph_i(ji,jl1) = phti(ji,jl2)
	1118	ph_s(ji,jl1) = phts(ji,jl2)
	1119	pa_i(ji,jl1) = pati(ji,jl2)
[11223]	1120	! record categories that are filled
	1121	jlmax(ji) = MAX( jlmax(ji), jl1 )
	1122	jlmin(ji) = MIN( jlmin(ji), jl1 )
	1123	jlfil(ji,jl1) = jl1
	1124	ENDIF
	1125	END DO
[8813]	1126	END DO
	1127	END DO
[11223]	1128	!
	1129	! --- fill the gaps between categories --- !
	1130	! transfer from categories filled at the previous step to the empty ones in between
[8813]	1131	DO ji = 1, idim
[11223]	1132	jl1 = jlmin(ji)
	1133	jl2 = jlmax(ji)
	1134	IF( jl1 > 1 ) THEN
	1135	! fill the lower cat (jl1-1)
[11332]	1136	pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1)
	1137	ph_i(ji,jl1-1) = hi_mean(jl1-1)
[11223]	1138	! remove from cat jl1
[11332]	1139	pa_i(ji,jl1 ) = ( 1._wp - ztrans ) * pa_i(ji,jl1)
[8813]	1140	ENDIF
[11223]	1141	IF( jl2 < jpl ) THEN
	1142	! fill the upper cat (jl2+1)
[11332]	1143	pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2)
	1144	ph_i(ji,jl2+1) = hi_mean(jl2+1)
[11223]	1145	! remove from cat jl2
[11332]	1146	pa_i(ji,jl2 ) = ( 1._wp - ztrans ) * pa_i(ji,jl2)
[8813]	1147	ENDIF
	1148	END DO
[11223]	1149	!
	1150	jlfil2(:,:) = jlfil(:,:)
	1151	! fill categories from low to high
	1152	DO jl = 2, jpl-1
	1153	DO ji = 1, idim
	1154	IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
	1155	! fill high
[11332]	1156	pa_i(ji,jl) = ztrans * pa_i(ji,jl-1)
	1157	ph_i(ji,jl) = hi_mean(jl)
[11223]	1158	jlfil(ji,jl) = jl
	1159	! remove low
[11332]	1160	pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1)
[11223]	1161	ENDIF
	1162	END DO
[8586]	1163	END DO
[11223]	1164	!
	1165	! fill categories from high to low
	1166	DO jl = jpl-1, 2, -1
	1167	DO ji = 1, idim
	1168	IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
	1169	! fill low
[11332]	1170	pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1)
	1171	ph_i(ji,jl) = hi_mean(jl)
[11223]	1172	jlfil2(ji,jl) = jl
	1173	! remove high
[11332]	1174	pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1)
[11223]	1175	ENDIF
	1176	END DO
[8984]	1177	END DO
[10332]	1178	!
[11223]	1179	DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays
	1180	DEALLOCATE( jlmin, jlmax )
[10332]	1181	!
[11332]	1182	! == temperature and salinity == !
	1183	!
	1184	IF( PRESENT( pt_i ) .OR. PRESENT( pt_s ) .OR. PRESENT( pt_su ) .OR. PRESENT( ps_i ) ) THEN
	1185	!
	1186	ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) )
	1187	!
	1188	WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp ) ; z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 )
	1189	ELSEWHERE ; z1_ai(:) = 0._wp
	1190	END WHERE
	1191	WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp ) ; z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 )
	1192	ELSEWHERE ; z1_vi(:) = 0._wp
	1193	END WHERE
	1194	WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp ) ; z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 )
	1195	ELSEWHERE ; z1_vs(:) = 0._wp
	1196	END WHERE
	1197	!
	1198	! fill all the categories with the same value
	1199	IF( PRESENT( pt_i ) ) THEN
	1200	ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
	1201	DO jl = 1, jpl
	1202	pt_i (:,jl) = ztmp(:)
	1203	END DO
	1204	ENDIF
	1205	IF( PRESENT( pt_s ) ) THEN
	1206	ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
	1207	DO jl = 1, jpl
	1208	pt_s (:,jl) = ztmp(:)
	1209	END DO
	1210	ENDIF
	1211	IF( PRESENT( pt_su ) ) THEN
	1212	ztmp(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:)
	1213	DO jl = 1, jpl
	1214	pt_su(:,jl) = ztmp(:)
	1215	END DO
	1216	ENDIF
	1217	IF( PRESENT( ps_i ) ) THEN
	1218	ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
	1219	DO jl = 1, jpl
	1220	ps_i (:,jl) = ztmp(:)
	1221	END DO
	1222	ENDIF
	1223	!
	1224	DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp )
	1225	!
	1226	ENDIF
	1227	!
[10332]	1228	ENDIF
	1229	!
[11223]	1230	END SUBROUTINE ice_var_itd_NcMc
[10332]	1231
[8586]	1232	#else
	1233	!!----------------------------------------------------------------------
[9570]	1234	!! Default option Dummy module NO SI3 sea-ice model
[8586]	1235	!!----------------------------------------------------------------------
	1236	#endif
	1237
	1238	!!======================================================================
	1239	END MODULE icevar

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

Download in other formats: