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/fix_vvl_ticket1791/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2019/fix_vvl_ticket1791/src/ICE/icevar.F90 @ 11422

Last change on this file since 11422 was 11422, checked in by jchanut, 5 years ago
#1791, merge with trunk
Property svn:keywords set to `Id`
File size: 51.5 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
[11422]	78	PUBLIC ice_var_itd
[8586]	79
[11422]	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
[11422]	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
[11422]	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	!!-------------------------------------------------------------------
	775	SUBROUTINE ice_var_itd_1c1c( zhti, zhts, zati, zh_i, zh_s, za_i )
[8586]	776	!!-------------------------------------------------------------------
[11422]	777	!! ** Purpose : converting 1-cat ice to 1 ice category
	778	!!-------------------------------------------------------------------
	779	REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
	780	REAL(wp), DIMENSION(:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
	781	!!-------------------------------------------------------------------
	782	zh_i(:) = zhti(:)
	783	zh_s(:) = zhts(:)
	784	za_i(:) = zati(:)
	785	END SUBROUTINE ice_var_itd_1c1c
	786
	787	SUBROUTINE ice_var_itd_Nc1c( zhti, zhts, zati, zh_i, zh_s, za_i )
	788	!!-------------------------------------------------------------------
	789	!! ** Purpose : converting N-cat ice to 1 ice category
	790	!!-------------------------------------------------------------------
	791	REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
	792	REAL(wp), DIMENSION(:) , INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
	793	!!-------------------------------------------------------------------
	794	!
	795	za_i(:) = SUM( zati(:,:), dim=2 )
	796	!
	797	WHERE( za_i(:) /= 0._wp )
	798	zh_i(:) = SUM( zhti(:,:) * zati(:,:), dim=2 ) / za_i(:)
	799	zh_s(:) = SUM( zhts(:,:) * zati(:,:), dim=2 ) / za_i(:)
	800	ELSEWHERE
	801	zh_i(:) = 0._wp
	802	zh_s(:) = 0._wp
	803	END WHERE
	804	!
	805	END SUBROUTINE ice_var_itd_Nc1c
	806
	807	SUBROUTINE ice_var_itd_1cMc( zhti, zhts, zati, zh_i, zh_s, za_i )
	808	!!-------------------------------------------------------------------
[8586]	809	!!
[11422]	810	!! ** Purpose : converting 1-cat ice to jpl ice categories
[8586]	811	!!
	812	!! ice thickness distribution follows a gaussian law
	813	!! around the concentration of the most likely ice thickness
	814	!! (similar as iceistate.F90)
	815	!!
	816	!! ** Method: Iterative procedure
	817	!!
	818	!! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
	819	!!
	820	!! 2) Check whether the distribution conserves area and volume, positivity and
	821	!! category boundaries
	822	!!
	823	!! 3) If not (input ice is too thin), the last category is empty and
	824	!! the number of categories is reduced (jpl-1)
	825	!!
	826	!! 4) Iterate until ok (SUM(itest(:) = 4)
	827	!!
	828	!! ** Arguments : zhti: 1-cat ice thickness
	829	!! zhts: 1-cat snow depth
	830	!! zati: 1-cat ice concentration
	831	!!
	832	!! ** Output : jpl-cat
	833	!!
	834	!! (Example of application: BDY forcings when input are cell averaged)
	835	!!-------------------------------------------------------------------
	836	INTEGER :: ji, jk, jl ! dummy loop indices
[8813]	837	INTEGER :: idim, i_fill, jl0
[8586]	838	REAL(wp) :: zarg, zV, zconv, zdh, zdv
[10994]	839	REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
	840	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
[8586]	841	INTEGER , DIMENSION(4) :: itest
	842	!!-------------------------------------------------------------------
	843	!
	844	! ----------------------------------------
	845	! distribution over the jpl ice categories
	846	! ----------------------------------------
	847	! a gaussian distribution for ice concentration is used
	848	! then we check whether the distribution fullfills
	849	! volume and area conservation, positivity and ice categories bounds
[8813]	850	idim = SIZE( zhti , 1 )
	851	zh_i(1:idim,1:jpl) = 0._wp
	852	zh_s(1:idim,1:jpl) = 0._wp
	853	za_i(1:idim,1:jpl) = 0._wp
	854	!
	855	DO ji = 1, idim
	856	!
[8586]	857	IF( zhti(ji) > 0._wp ) THEN
[8813]	858	!
[8586]	859	! find which category (jl0) the input ice thickness falls into
	860	jl0 = jpl
	861	DO jl = 1, jpl
	862	IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
	863	jl0 = jl
	864	CYCLE
	865	ENDIF
	866	END DO
[8813]	867	!
[8586]	868	itest(:) = 0
	869	i_fill = jpl + 1 !------------------------------------
	870	DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories
	871	! !------------------------------------
	872	i_fill = i_fill - 1
	873	!
	874	zh_i(ji,1:jpl) = 0._wp
[8813]	875	za_i(ji,1:jpl) = 0._wp
	876	itest(:) = 0
	877	!
[8586]	878	IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1
	879	zh_i(ji,1) = zhti(ji)
	880	za_i (ji,1) = zati (ji)
	881	ELSE !-- case ice is thicker: fill categories >1
	882	! thickness
	883	DO jl = 1, i_fill - 1
	884	zh_i(ji,jl) = hi_mean(jl)
	885	END DO
[8813]	886	!
[8586]	887	! concentration
	888	za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))
	889	DO jl = 1, i_fill - 1
	890	IF ( jl /= jl0 ) THEN
	891	zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
	892	za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)
	893	ENDIF
	894	END DO
[8813]	895	!
[8586]	896	! last category
	897	za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )
	898	zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )
	899	zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )
[8813]	900	!
[8885]	901	! correction if concentration of upper cat is greater than lower cat
	902	! (it should be a gaussian around jl0 but sometimes it is not)
[8586]	903	IF ( jl0 /= jpl ) THEN
	904	DO jl = jpl, jl0+1, -1
	905	IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
	906	zdv = zh_i(ji,jl) * za_i(ji,jl)
	907	zh_i(ji,jl ) = 0._wp
	908	za_i (ji,jl ) = 0._wp
	909	za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
	910	END IF
[8813]	911	END DO
[8586]	912	ENDIF
[8813]	913	!
[8586]	914	ENDIF
[8813]	915	!
[8586]	916	! Compatibility tests
	917	zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )
[8813]	918	IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation
	919	!
[8586]	920	zconv = ABS( zhti(ji)zati(ji) - SUM( za_i(ji,1:jpl)zh_i(ji,1:jpl) ) )
[8813]	921	IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation
	922	!
	923	IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?
	924	!
[8586]	925	itest(4) = 1
	926	DO jl = 1, i_fill
	927	IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations
	928	END DO
	929	! !----------------------------
	930	END DO ! end iteration on categories
[8813]	931	! !----------------------------
[8586]	932	ENDIF
	933	END DO
	934
	935	! Add Snow in each category where za_i is not 0
	936	DO jl = 1, jpl
[8813]	937	DO ji = 1, idim
[8586]	938	IF( za_i(ji,jl) > 0._wp ) THEN
	939	zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )
	940	! In case snow load is in excess that would lead to transformation from snow to ice
	941	! Then, transfer the snow excess into the ice (different from icethd_dh)
[9935]	942	zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )
[8586]	943	! recompute h_i, h_s avoiding out of bounds values
	944	zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )
[9935]	945	zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )
[8586]	946	ENDIF
	947	END DO
	948	END DO
	949	!
[11422]	950	END SUBROUTINE ice_var_itd_1cMc
[8586]	951
[11422]	952	SUBROUTINE ice_var_itd_NcMc( zhti, zhts, zati, zh_i, zh_s, za_i )
[8586]	953	!!-------------------------------------------------------------------
[8813]	954	!!
	955	!! ** Purpose : converting N-cat ice to jpl ice categories
	956	!!
	957	!! ice thickness distribution follows a gaussian law
	958	!! around the concentration of the most likely ice thickness
	959	!! (similar as iceistate.F90)
	960	!!
	961	!! ** Method: Iterative procedure
	962	!!
	963	!! 1) Fill ice cat that correspond to input thicknesses
	964	!! Find the lowest(jlmin) and highest(jlmax) cat that are filled
	965	!!
	966	!! 2) Expand the filling to the cat jlmin-1 and jlmax+1
[10994]	967	!! by removing 25% ice area from jlmin and jlmax (resp.)
[8813]	968	!!
	969	!! 3) Expand the filling to the empty cat between jlmin and jlmax
	970	!! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
	971	!! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
	972	!!
	973	!! ** Arguments : zhti: N-cat ice thickness
	974	!! zhts: N-cat snow depth
	975	!! zati: N-cat ice concentration
	976	!!
	977	!! ** Output : jpl-cat
	978	!!
	979	!! (Example of application: BDY forcings when inputs have N-cat /= jpl)
	980	!!-------------------------------------------------------------------
	981	INTEGER :: ji, jl, jl1, jl2 ! dummy loop indices
	982	INTEGER :: idim, icat
[10589]	983	REAL(wp), PARAMETER :: ztrans = 0.25_wp
[8813]	984	REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables
	985	REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables
	986	INTEGER , DIMENSION(:,:), ALLOCATABLE :: jlfil, jlfil2
	987	INTEGER , DIMENSION(:) , ALLOCATABLE :: jlmax, jlmin
	988	!!-------------------------------------------------------------------
	989	!
	990	idim = SIZE( zhti, 1 )
	991	icat = SIZE( zhti, 2 )
[11422]	992	! ! ---------------------- !
	993	IF( icat == jpl ) THEN ! input cat = output cat !
	994	! ! ---------------------- !
	995	zh_i(:,:) = zhti(:,:)
	996	zh_s(:,:) = zhts(:,:)
	997	za_i(:,:) = zati(:,:)
	998	! ! ---------------------- !
	999	ELSEIF( icat == 1 ) THEN ! input cat = 1 !
	1000	! ! ---------------------- !
	1001	CALL ice_var_itd_1cMc( zhti(:,1), zhts(:,1), zati(:,1), zh_i(:,:), zh_s(:,:), za_i(:,:) )
	1002	! ! ---------------------- !
	1003	ELSEIF( jpl == 1 ) THEN ! output cat = 1 !
	1004	! ! ---------------------- !
	1005	CALL ice_var_itd_Nc1c( zhti(:,:), zhts(:,:), zati(:,:), zh_i(:,1), zh_s(:,1), za_i(:,1) )
	1006	! ! ----------------------- !
	1007	ELSE ! input cat /= output cat !
	1008	! ! ----------------------- !
	1009
	1010	ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays
	1011	ALLOCATE( jlmin(idim), jlmax(idim) )
[8813]	1012
[11422]	1013	! --- initialize output fields to 0 --- !
	1014	zh_i(1:idim,1:jpl) = 0._wp
	1015	zh_s(1:idim,1:jpl) = 0._wp
	1016	za_i(1:idim,1:jpl) = 0._wp
	1017	!
	1018	! --- fill the categories --- !
	1019	! find where cat-input = cat-output and fill cat-output fields
	1020	jlmax(:) = 0
	1021	jlmin(:) = 999
	1022	jlfil(:,:) = 0
	1023	DO jl1 = 1, jpl
	1024	DO jl2 = 1, icat
	1025	DO ji = 1, idim
	1026	IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
	1027	! fill the right category
	1028	zh_i(ji,jl1) = zhti(ji,jl2)
	1029	zh_s(ji,jl1) = zhts(ji,jl2)
	1030	za_i(ji,jl1) = zati(ji,jl2)
	1031	! record categories that are filled
	1032	jlmax(ji) = MAX( jlmax(ji), jl1 )
	1033	jlmin(ji) = MIN( jlmin(ji), jl1 )
	1034	jlfil(ji,jl1) = jl1
	1035	ENDIF
	1036	END DO
[8813]	1037	END DO
	1038	END DO
[11422]	1039	!
	1040	! --- fill the gaps between categories --- !
	1041	! transfer from categories filled at the previous step to the empty ones in between
[8813]	1042	DO ji = 1, idim
[11422]	1043	jl1 = jlmin(ji)
	1044	jl2 = jlmax(ji)
	1045	IF( jl1 > 1 ) THEN
	1046	! fill the lower cat (jl1-1)
	1047	za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
	1048	zh_i(ji,jl1-1) = hi_mean(jl1-1)
	1049	! remove from cat jl1
	1050	za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
[8813]	1051	ENDIF
[11422]	1052	IF( jl2 < jpl ) THEN
	1053	! fill the upper cat (jl2+1)
	1054	za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
	1055	zh_i(ji,jl2+1) = hi_mean(jl2+1)
	1056	! remove from cat jl2
	1057	za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
[8813]	1058	ENDIF
	1059	END DO
[11422]	1060	!
	1061	jlfil2(:,:) = jlfil(:,:)
	1062	! fill categories from low to high
	1063	DO jl = 2, jpl-1
	1064	DO ji = 1, idim
	1065	IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
	1066	! fill high
	1067	za_i(ji,jl) = ztrans * za_i(ji,jl-1)
	1068	zh_i(ji,jl) = hi_mean(jl)
	1069	jlfil(ji,jl) = jl
	1070	! remove low
	1071	za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
	1072	ENDIF
	1073	END DO
[8586]	1074	END DO
[11422]	1075	!
	1076	! fill categories from high to low
	1077	DO jl = jpl-1, 2, -1
	1078	DO ji = 1, idim
	1079	IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
	1080	! fill low
	1081	za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
	1082	zh_i(ji,jl) = hi_mean(jl)
	1083	jlfil2(ji,jl) = jl
	1084	! remove high
	1085	za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
	1086	ENDIF
	1087	END DO
[8984]	1088	END DO
[10332]	1089	!
[11422]	1090	DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays
	1091	DEALLOCATE( jlmin, jlmax )
[10332]	1092	!
	1093	ENDIF
	1094	!
[11422]	1095	END SUBROUTINE ice_var_itd_NcMc
[10332]	1096
[8586]	1097	#else
	1098	!!----------------------------------------------------------------------
[9570]	1099	!! Default option Dummy module NO SI3 sea-ice model
[8586]	1100	!!----------------------------------------------------------------------
	1101	#endif
	1102
	1103	!!======================================================================
	1104	END MODULE icevar

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

Download in other formats: