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_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/ICE/icevar.F90 @ 11413

Last change on this file since 11413 was 11413, checked in by gsamson, 5 years ago

dev_r11265_ABL : see #2131

merge src and cfgs from HPC-13_IRRMANN_BDY_optimization branch @ r11402 with dev_r11265_ABL branch @ r11363
change ORCA2 results due to ice rheology "cleaning" (see commit r11377)

Property svn:keywords set to Id

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

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

Download in other formats: