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.

limthd_dif.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 8420

Last change on this file since 8420 was 8420, checked in by clem, 7 years ago
changing names
Property svn:keywords set to `Id`
File size: 41.1 KB

Rev	Line
[825]	1	MODULE limthd_dif
	2	!!======================================================================
	3	!! * MODULE limthd_dif *
	4	!! heat diffusion in sea ice
	5	!! computation of surface and inner T
	6	!!======================================================================
[2715]	7	!! History : LIM ! 02-2003 (M. Vancoppenolle) original 1D code
	8	!! ! 06-2005 (M. Vancoppenolle) 3d version
	9	!! ! 11-2006 (X Fettweis) Vectorization by Xavier
	10	!! ! 04-2007 (M. Vancoppenolle) Energy conservation
	11	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
[4161]	12	!! - ! 2012-05 (C. Rousset) add penetration solar flux
[825]	13	!!----------------------------------------------------------------------
[2528]	14	#if defined key_lim3
	15	!!----------------------------------------------------------------------
	16	!! 'key_lim3' LIM3 sea-ice model
	17	!!----------------------------------------------------------------------
[3625]	18	USE par_oce ! ocean parameters
	19	USE phycst ! physical constants (ocean directory)
	20	USE ice ! LIM-3 variables
[8420]	21	USE ice1D ! LIM-3: thermodynamics
[8378]	22	!
[3625]	23	USE in_out_manager ! I/O manager
	24	USE lib_mpp ! MPP library
	25	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[921]	26
[825]	27	IMPLICIT NONE
	28	PRIVATE
	29
[2528]	30	PUBLIC lim_thd_dif ! called by lim_thd
[825]	31
	32	!!----------------------------------------------------------------------
[4161]	33	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	34	!! $Id$
[2591]	35	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	36	!!----------------------------------------------------------------------
	37	CONTAINS
	38
[8342]	39	SUBROUTINE lim_thd_dif
[921]	40	!!------------------------------------------------------------------
	41	!! * ROUTINE lim_thd_dif *
	42	!! ** Purpose :
	43	!! This routine determines the time evolution of snow and sea-ice
	44	!! temperature profiles.
	45	!! ** Method :
	46	!! This is done by solving the heat equation diffusion with
	47	!! a Neumann boundary condition at the surface and a Dirichlet one
	48	!! at the bottom. Solar radiation is partially absorbed into the ice.
	49	!! The specific heat and thermal conductivities depend on ice salinity
	50	!! and temperature to take into account brine pocket melting. The
	51	!! numerical
	52	!! scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid
	53	!! in the ice and snow system.
	54	!!
	55	!! The successive steps of this routine are
	56	!! 1. Thermal conductivity at the interfaces of the ice layers
	57	!! 2. Internal absorbed radiation
	58	!! 3. Scale factors due to non-uniform grid
	59	!! 4. Kappa factors
	60	!! Then iterative procedure begins
	61	!! 5. specific heat in the ice
	62	!! 6. eta factors
	63	!! 7. surface flux computation
	64	!! 8. tridiagonal system terms
	65	!! 9. solving the tridiagonal system with Gauss elimination
	66	!! Iterative procedure ends according to a criterion on evolution
	67	!! of temperature
	68	!!
	69	!!
	70	!! ** Inputs / Ouputs : (global commons)
[4872]	71	!! surface temperature : t_su_1d
	72	!! ice/snow temperatures : t_i_1d, t_s_1d
	73	!! ice salinities : s_i_1d
[921]	74	!! number of layers in the ice/snow: nlay_i, nlay_s
	75	!! profile of the ice/snow layers : z_i, z_s
[4872]	76	!! total ice/snow thickness : ht_i_1d, ht_s_1d
[921]	77	!!
	78	!! ** External :
	79	!!
	80	!! ** References :
	81	!!
	82	!! ** History :
	83	!! (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium
	84	!! (06-2005) Martin Vancoppenolle, 3d version
	85	!! (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR)
	86	!! (04-2007) Energy conservation tested by M. Vancoppenolle
	87	!!------------------------------------------------------------------
	88	!! * Local variables
[8313]	89	INTEGER :: ji ! spatial loop index
	90	INTEGER :: ii, ij ! temporary dummy loop index
	91	INTEGER :: numeq ! current reference number of equation
	92	INTEGER :: jk ! vertical dummy loop index
[3294]	93	INTEGER :: minnumeqmin, maxnumeqmax
[8313]	94	INTEGER :: iconv ! number of iterations in iterative procedure
	95	INTEGER :: iconv_max = 50 ! max number of iterations in iterative procedure
[5123]	96
[8373]	97	INTEGER, DIMENSION(jpij) :: numeqmin ! reference number of top equation
	98	INTEGER, DIMENSION(jpij) :: numeqmax ! reference number of bottom equation
[5123]	99
[3294]	100	REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system
	101	REAL(wp) :: zg1 = 2._wp !
	102	REAL(wp) :: zgamma = 18009._wp ! for specific heat
	103	REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13)
[4990]	104	REAL(wp) :: zraext_s = 10._wp ! extinction coefficient of radiation in the snow
[3294]	105	REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity
[8325]	106	REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C
[8313]	107	REAL(wp) :: ztmelt_i ! ice melting temperature
[5123]	108	REAL(wp) :: zhsu
[8313]	109	REAL(wp) :: zdti_max ! current maximal error on temperature
	110	REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature
[5123]	111
[8373]	112	REAL(wp), DIMENSION(jpij) :: isnow ! switch for presence (1) or absence (0) of snow
	113	REAL(wp), DIMENSION(jpij) :: ztsub ! old surface temperature (before the iterative procedure )
	114	REAL(wp), DIMENSION(jpij) :: ztsubit ! surface temperature at previous iteration
	115	REAL(wp), DIMENSION(jpij) :: zh_i ! ice layer thickness
	116	REAL(wp), DIMENSION(jpij) :: zh_s ! snow layer thickness
	117	REAL(wp), DIMENSION(jpij) :: zfsw ! solar radiation absorbed at the surface
	118	REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface
	119	REAL(wp), DIMENSION(jpij) :: zf ! surface flux function
	120	REAL(wp), DIMENSION(jpij) :: dzf ! derivative of the surface flux function
	121	REAL(wp), DIMENSION(jpij) :: zdti ! current error on temperature
	122	REAL(wp), DIMENSION(jpij) :: zdifcase ! case of the equation resolution (1->4)
	123	REAL(wp), DIMENSION(jpij) :: zftrice ! solar radiation transmitted through the ice
	124	REAL(wp), DIMENSION(jpij) :: zihic
[5123]	125
[8373]	126	REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity
	127	REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice
	128	REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice
	129	REAL(wp), DIMENSION(jpij,0:nlay_i) :: zkappa_i ! Kappa factor in the ice
	130	REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztib ! Old temperature in the ice
	131	REAL(wp), DIMENSION(jpij,0:nlay_i) :: zeta_i ! Eta factor in the ice
	132	REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztitemp ! Temporary temperature in the ice to check the convergence
	133	REAL(wp), DIMENSION(jpij,0:nlay_i) :: zspeche_i ! Ice specific heat
	134	REAL(wp), DIMENSION(jpij,0:nlay_i) :: z_i ! Vertical cotes of the layers in the ice
	135	REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow
	136	REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow
	137	REAL(wp), DIMENSION(jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow
	138	REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow
	139	REAL(wp), DIMENSION(jpij,0:nlay_s) :: ztstemp ! Temporary temperature in the snow to check the convergence
	140	REAL(wp), DIMENSION(jpij,0:nlay_s) :: ztsb ! Temporary temperature in the snow
	141	REAL(wp), DIMENSION(jpij,0:nlay_s) :: z_s ! Vertical cotes of the layers in the snow
	142	REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindterm ! 'Ind'ependent term
	143	REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindtbis ! Temporary 'ind'ependent term
	144	REAL(wp), DIMENSION(jpij,nlay_i+3) :: zdiagbis ! Temporary 'dia'gonal term
	145	REAL(wp), DIMENSION(jpij,nlay_i+3,3) :: ztrid ! Tridiagonal system terms
	146	REAL(wp), DIMENSION(jpij) :: zdq, zq_ini, zhfx_err ! diag errors on heat
	147	REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat
[5123]	148
	149	! Mono-category
[8373]	150	REAL(wp) :: zepsilon ! determines thres. above which computation of G(h) is done
	151	REAL(wp) :: zratio_s ! dummy factor
	152	REAL(wp) :: zratio_i ! dummy factor
	153	REAL(wp) :: zh_thres ! thickness thres. for G(h) computation
	154	REAL(wp) :: zhe ! dummy factor
	155	REAL(wp) :: zkimean ! mean sea ice thermal conductivity
	156	REAL(wp) :: zfac ! dummy factor
	157	REAL(wp) :: zihe ! dummy factor
	158	REAL(wp) :: zheshth ! dummy factor
[3625]	159	!!------------------------------------------------------------------
[4688]	160
	161	! --- diag error on heat diffusion - PART 1 --- !
	162	zdq(:) = 0._wp ; zq_ini(:) = 0._wp
[8342]	163	DO ji = 1, nidx
[8325]	164	zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + &
	165	& SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
[4688]	166	END DO
	167
[921]	168	!------------------------------------------------------------------------------!
	169	! 1) Initialization !
	170	!------------------------------------------------------------------------------!
[8342]	171	DO ji = 1 , nidx
[5123]	172	isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) ) ! is there snow or not
[921]	173	! layer thickness
[5123]	174	zh_i(ji) = ht_i_1d(ji) * r1_nlay_i
	175	zh_s(ji) = ht_s_1d(ji) * r1_nlay_s
[921]	176	END DO
[825]	177
[921]	178	!--------------------
	179	! Ice / snow layers
	180	!--------------------
[825]	181
[2715]	182	z_s(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st snow layer
	183	z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer
[825]	184
[4870]	185	DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer
[8342]	186	DO ji = 1 , nidx
[5123]	187	z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s
[921]	188	END DO
	189	END DO
[825]	190
[4870]	191	DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer
[8342]	192	DO ji = 1 , nidx
[5123]	193	z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i
[921]	194	END DO
	195	END DO
	196	!
	197	!------------------------------------------------------------------------------\|
[5123]	198	! 2) Radiation \|
[921]	199	!------------------------------------------------------------------------------\|
	200	!
	201	!-------------------
	202	! Computation of i0
	203	!-------------------
	204	! i0 describes the fraction of solar radiation which does not contribute
	205	! to the surface energy budget but rather penetrates inside the ice.
	206	! We assume that no radiation is transmitted through the snow
	207	! If there is no no snow
	208	! zfsw = (1-i0).qsr_ice is absorbed at the surface
	209	! zftrice = io.qsr_ice is below the surface
[4688]	210	! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
[5123]	211	! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
	212	zhsu = 0.1_wp ! threshold for the computation of i0
[8342]	213	DO ji = 1 , nidx
[921]	214	! switches
[5123]	215	isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )
[921]	216	! hs > 0, isnow = 1
[5123]	217	zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) )
[825]	218
[5123]	219	i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
[921]	220	END DO
[825]	221
[921]	222	!-------------------------------------------------------
	223	! Solar radiation absorbed / transmitted at the surface
	224	! Derivative of the non solar flux
	225	!-------------------------------------------------------
[8342]	226	DO ji = 1 , nidx
[5146]	227	zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface
	228	zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer
	229	dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux
	230	zqns_ice_b(ji) = qns_ice_1d(ji) ! store previous qns_ice_1d value
[921]	231	END DO
[825]	232
[921]	233	!---------------------------------------------------------
	234	! Transmission - absorption of solar radiation in the ice
	235	!---------------------------------------------------------
[825]	236
[8342]	237	DO ji = 1, nidx ! snow initialization
[2715]	238	zradtr_s(ji,0) = zftrice(ji) ! radiation penetrating through snow
[921]	239	END DO
[825]	240
[4870]	241	DO jk = 1, nlay_s ! Radiation through snow
[8342]	242	DO ji = 1, nidx
[2715]	243	! ! radiation transmitted below the layer-th snow layer
[4870]	244	zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) )
[2715]	245	! ! radiation absorbed by the layer-th snow layer
[4870]	246	zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
[921]	247	END DO
	248	END DO
[825]	249
[8342]	250	DO ji = 1, nidx ! ice initialization
[5123]	251	zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) )
[921]	252	END DO
[825]	253
[4870]	254	DO jk = 1, nlay_i ! Radiation through ice
[8342]	255	DO ji = 1, nidx
[2715]	256	! ! radiation transmitted below the layer-th ice layer
[5123]	257	zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) )
[2715]	258	! ! radiation absorbed by the layer-th ice layer
[4870]	259	zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
[921]	260	END DO
	261	END DO
[825]	262
[8342]	263	DO ji = 1, nidx ! Radiation transmitted below the ice
[4688]	264	ftr_ice_1d(ji) = zradtr_i(ji,nlay_i)
[921]	265	END DO
[834]	266
[921]	267	!------------------------------------------------------------------------------\|
	268	! 3) Iterative procedure begins \|
	269	!------------------------------------------------------------------------------\|
	270	!
[8342]	271	DO ji = 1, nidx ! Old surface temperature
[4872]	272	ztsub (ji) = t_su_1d(ji) ! temperature at the beg of iter pr.
	273	ztsubit(ji) = t_su_1d(ji) ! temperature at the previous iter
[5123]	274	t_su_1d(ji) = MIN( t_su_1d(ji), rt0 - ztsu_err ) ! necessary
[8313]	275	zdti (ji) = 1000._wp ! initial value of error
[921]	276	END DO
[825]	277
[4870]	278	DO jk = 1, nlay_s ! Old snow temperature
[8342]	279	DO ji = 1 , nidx
[4872]	280	ztsb(ji,jk) = t_s_1d(ji,jk)
[921]	281	END DO
	282	END DO
[825]	283
[4870]	284	DO jk = 1, nlay_i ! Old ice temperature
[8342]	285	DO ji = 1 , nidx
[4872]	286	ztib(ji,jk) = t_i_1d(ji,jk)
[921]	287	END DO
	288	END DO
[825]	289
[8313]	290	iconv = 0 ! number of iterations
	291	zdti_max = 1000._wp ! maximal value of error on all points
[825]	292
[8313]	293	DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )
[921]	294	!
[8313]	295	iconv = iconv + 1
[2715]	296	!
[921]	297	!------------------------------------------------------------------------------\|
	298	! 4) Sea ice thermal conductivity \|
	299	!------------------------------------------------------------------------------\|
	300	!
[5123]	301	IF( nn_ice_thcon == 0 ) THEN ! Untersteiner (1964) formula
[8342]	302	DO ji = 1 , nidx
[5123]	303	ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
	304	ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
[921]	305	END DO
[4870]	306	DO jk = 1, nlay_i-1
[8342]	307	DO ji = 1 , nidx
[5123]	308	ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / &
	309	MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0)
	310	ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
[921]	311	END DO
	312	END DO
	313	ENDIF
[825]	314
[5123]	315	IF( nn_ice_thcon == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
[8342]	316	DO ji = 1 , nidx
[5123]	317	ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) &
	318	& - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
[2715]	319	ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
[921]	320	END DO
[4870]	321	DO jk = 1, nlay_i-1
[8342]	322	DO ji = 1 , nidx
[5123]	323	ztcond_i(ji,jk) = rcdic + &
	324	& 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) &
	325	& / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 ) &
	326	& - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 )
[4870]	327	ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
[2715]	328	END DO
	329	END DO
[8342]	330	DO ji = 1 , nidx
[5123]	331	ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) &
	332	& - 0.011_wp * ( t_bo_1d(ji) - rt0 )
[2715]	333	ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
	334	END DO
[921]	335	ENDIF
[5123]	336
[921]	337	!
	338	!------------------------------------------------------------------------------\|
[5123]	339	! 5) G(he) - enhancement of thermal conductivity in mono-category case \|
[921]	340	!------------------------------------------------------------------------------\|
	341	!
[5123]	342	! Computation of effective thermal conductivity G(h)
	343	! Used in mono-category case only to simulate an ITD implicitly
	344	! Fichefet and Morales Maqueda, JGR 1997
	345
	346	zghe(:) = 1._wp
	347
	348	SELECT CASE ( nn_monocat )
	349
	350	CASE (1,3) ! LIM3
	351
	352	zepsilon = 0.1_wp
	353	zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp
	354
[8342]	355	DO ji = 1, nidx
[5123]	356
	357	! Mean sea ice thermal conductivity
	358	zkimean = SUM( ztcond_i(ji,0:nlay_i) ) / REAL( nlay_i+1, wp )
	359
	360	! Effective thickness he (zhe)
[7813]	361	zfac = 1._wp / ( rn_cdsn + zkimean )
	362	zratio_s = rn_cdsn * zfac
[5123]	363	zratio_i = zkimean * zfac
	364	zhe = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji)
	365
	366	! G(he)
	367	rswitch = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) ) ! =0 if zhe < zh_thres, if >
	368	zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) )
	369
	370	! Impose G(he) < 2.
	371	zghe(ji) = MIN( zghe(ji), 2._wp )
	372
	373	END DO
	374
	375	END SELECT
	376
	377	!
	378	!------------------------------------------------------------------------------\|
	379	! 6) kappa factors \|
	380	!------------------------------------------------------------------------------\|
	381	!
	382	!--- Snow
[8342]	383	DO ji = 1, nidx
[5123]	384	zfac = 1. / MAX( epsi10 , zh_s(ji) )
[7813]	385	zkappa_s(ji,0) = zghe(ji) * rn_cdsn * zfac
	386	zkappa_s(ji,nlay_s) = zghe(ji) * rn_cdsn * zfac
[921]	387	END DO
[825]	388
[4870]	389	DO jk = 1, nlay_s-1
[8342]	390	DO ji = 1 , nidx
[7813]	391	zkappa_s(ji,jk) = zghe(ji) * 2.0 * rn_cdsn / MAX( epsi10, 2.0 * zh_s(ji) )
[921]	392	END DO
	393	END DO
[825]	394
[5123]	395	!--- Ice
[4870]	396	DO jk = 1, nlay_i-1
[8342]	397	DO ji = 1 , nidx
[5123]	398	zkappa_i(ji,jk) = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0 * zh_i(ji) )
[921]	399	END DO
	400	END DO
[825]	401
[5123]	402	!--- Snow-ice interface
[8342]	403	DO ji = 1 , nidx
[5123]	404	zfac = 1./ MAX( epsi10 , zh_i(ji) )
	405	zkappa_i(ji,0) = zghe(ji) * ztcond_i(ji,0) * zfac
	406	zkappa_i(ji,nlay_i) = zghe(ji) * ztcond_i(ji,nlay_i) * zfac
[7813]	407	zkappa_s(ji,nlay_s) = zghe(ji) * zghe(ji) * 2.0 * rn_cdsn * ztcond_i(ji,0) / &
	408	& MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rn_cdsn * zh_i(ji) ) )
[5123]	409	zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
[921]	410	END DO
[5123]	411
[921]	412	!
	413	!------------------------------------------------------------------------------\|
[5123]	414	! 7) Sea ice specific heat, eta factors \|
[921]	415	!------------------------------------------------------------------------------\|
	416	!
[4870]	417	DO jk = 1, nlay_i
[8342]	418	DO ji = 1 , nidx
[4872]	419	ztitemp(ji,jk) = t_i_1d(ji,jk)
[5123]	420	zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 )
	421	zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 )
[921]	422	END DO
	423	END DO
[825]	424
[4870]	425	DO jk = 1, nlay_s
[8342]	426	DO ji = 1 , nidx
[4872]	427	ztstemp(ji,jk) = t_s_1d(ji,jk)
[5123]	428	zeta_s(ji,jk) = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 )
[921]	429	END DO
	430	END DO
[5123]	431
[921]	432	!
	433	!------------------------------------------------------------------------------\|
[5123]	434	! 8) surface flux computation \|
[921]	435	!------------------------------------------------------------------------------\|
	436	!
[8313]	437	IF ( ln_dqnsice ) THEN
[8342]	438	DO ji = 1 , nidx
[4990]	439	! update of the non solar flux according to the update in T_su
	440	qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) )
	441	END DO
	442	ENDIF
	443
	444	! Update incoming flux
[8342]	445	DO ji = 1 , nidx
[921]	446	! update incoming flux
[5123]	447	zf(ji) = zfsw(ji) & ! net absorbed solar radiation
	448	& + qns_ice_1d(ji) ! non solar total flux (LWup, LWdw, SH, LH)
[921]	449	END DO
[825]	450
[921]	451	!
	452	!------------------------------------------------------------------------------\|
[5123]	453	! 9) tridiagonal system terms \|
[921]	454	!------------------------------------------------------------------------------\|
	455	!
	456	!!layer denotes the number of the layer in the snow or in the ice
	457	!!numeq denotes the reference number of the equation in the tridiagonal
	458	!!system, terms of tridiagonal system are indexed as following :
	459	!!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
[825]	460
[921]	461	!!ice interior terms (top equation has the same form as the others)
	462
[4873]	463	DO numeq=1,nlay_i+3
[8342]	464	DO ji = 1 , nidx
[921]	465	ztrid(ji,numeq,1) = 0.
	466	ztrid(ji,numeq,2) = 0.
	467	ztrid(ji,numeq,3) = 0.
[5123]	468	zindterm(ji,numeq)= 0.
	469	zindtbis(ji,numeq)= 0.
[921]	470	zdiagbis(ji,numeq)= 0.
	471	ENDDO
	472	ENDDO
	473
	474	DO numeq = nlay_s + 2, nlay_s + nlay_i
[8342]	475	DO ji = 1 , nidx
[5123]	476	jk = numeq - nlay_s - 1
	477	ztrid(ji,numeq,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
	478	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
	479	ztrid(ji,numeq,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk)
	480	zindterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
[921]	481	END DO
	482	ENDDO
	483
	484	numeq = nlay_s + nlay_i + 1
[8342]	485	DO ji = 1 , nidx
[825]	486	!!ice bottom term
	487	ztrid(ji,numeq,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)
[5123]	488	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
[825]	489	ztrid(ji,numeq,3) = 0.0
[5123]	490	zindterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i) * &
	491	& ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )
[921]	492	ENDDO
[825]	493
	494
[8342]	495	DO ji = 1 , nidx
[5123]	496	IF ( ht_s_1d(ji) > 0.0 ) THEN
[921]	497	!
	498	!------------------------------------------------------------------------------\|
	499	! snow-covered cells \|
	500	!------------------------------------------------------------------------------\|
	501	!
	502	!!snow interior terms (bottom equation has the same form as the others)
	503	DO numeq = 3, nlay_s + 1
[5123]	504	jk = numeq - 1
	505	ztrid(ji,numeq,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
	506	ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
[4870]	507	ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk)
[5123]	508	zindterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
[921]	509	END DO
[825]	510
[921]	511	!!case of only one layer in the ice (ice equation is altered)
	512	IF ( nlay_i.eq.1 ) THEN
	513	ztrid(ji,nlay_s+2,3) = 0.0
[5123]	514	zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)
[921]	515	ENDIF
[834]	516
[5123]	517	IF ( t_su_1d(ji) < rt0 ) THEN
[825]	518
[921]	519	!------------------------------------------------------------------------------\|
	520	! case 1 : no surface melting - snow present \|
	521	!------------------------------------------------------------------------------\|
	522	zdifcase(ji) = 1.0
	523	numeqmin(ji) = 1
	524	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	525
[921]	526	!!surface equation
[5123]	527	ztrid(ji,1,1) = 0.0
	528	ztrid(ji,1,2) = dzf(ji) - zg1s * zkappa_s(ji,0)
	529	ztrid(ji,1,3) = zg1s * zkappa_s(ji,0)
	530	zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji)
[825]	531
[921]	532	!!first layer of snow equation
[5123]	533	ztrid(ji,2,1) = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
	534	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
[921]	535	ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1)
[5123]	536	zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
[825]	537
[921]	538	ELSE
	539	!
	540	!------------------------------------------------------------------------------\|
	541	! case 2 : surface is melting - snow present \|
	542	!------------------------------------------------------------------------------\|
	543	!
	544	zdifcase(ji) = 2.0
	545	numeqmin(ji) = 2
	546	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	547
[921]	548	!!first layer of snow equation
	549	ztrid(ji,2,1) = 0.0
[5123]	550	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
[921]	551	ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1)
[5123]	552	zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * &
	553	& ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )
[921]	554	ENDIF
	555	ELSE
	556	!
	557	!------------------------------------------------------------------------------\|
	558	! cells without snow \|
	559	!------------------------------------------------------------------------------\|
	560	!
[5123]	561	IF ( t_su_1d(ji) < rt0 ) THEN
[921]	562	!
	563	!------------------------------------------------------------------------------\|
	564	! case 3 : no surface melting - no snow \|
	565	!------------------------------------------------------------------------------\|
	566	!
	567	zdifcase(ji) = 3.0
	568	numeqmin(ji) = nlay_s + 1
	569	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	570
[921]	571	!!surface equation
	572	ztrid(ji,numeqmin(ji),1) = 0.0
	573	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1
	574	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1
[5123]	575	zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji)
[825]	576
[921]	577	!!first layer of ice equation
	578	ztrid(ji,numeqmin(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
[5123]	579	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
	580	ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1)
	581	zindterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)
[825]	582
[921]	583	!!case of only one layer in the ice (surface & ice equations are altered)
[825]	584
[5123]	585	IF ( nlay_i == 1 ) THEN
[921]	586	ztrid(ji,numeqmin(ji),1) = 0.0
[5123]	587	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0) * 2.0
	588	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0) * 2.0
	589	ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
	590	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
[921]	591	ztrid(ji,numeqmin(ji)+1,3) = 0.0
[825]	592
[5123]	593	zindterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1) * &
	594	& ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
[921]	595	ENDIF
[825]	596
[921]	597	ELSE
[825]	598
[921]	599	!
	600	!------------------------------------------------------------------------------\|
	601	! case 4 : surface is melting - no snow \|
	602	!------------------------------------------------------------------------------\|
	603	!
	604	zdifcase(ji) = 4.0
	605	numeqmin(ji) = nlay_s + 2
	606	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	607
[921]	608	!!first layer of ice equation
	609	ztrid(ji,numeqmin(ji),1) = 0.0
[5123]	610	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
[921]	611	ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)
[5123]	612	zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1) * &
	613	& ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) )
[825]	614
[921]	615	!!case of only one layer in the ice (surface & ice equations are altered)
[5123]	616	IF ( nlay_i == 1 ) THEN
[921]	617	ztrid(ji,numeqmin(ji),1) = 0.0
[5123]	618	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
[921]	619	ztrid(ji,numeqmin(ji),3) = 0.0
[5123]	620	zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &
	621	& + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
[921]	622	ENDIF
[825]	623
[921]	624	ENDIF
	625	ENDIF
[825]	626
[921]	627	END DO
[825]	628
[921]	629	!
	630	!------------------------------------------------------------------------------\|
[5123]	631	! 10) tridiagonal system solving \|
[921]	632	!------------------------------------------------------------------------------\|
	633	!
[825]	634
[921]	635	! Solve the tridiagonal system with Gauss elimination method.
	636	! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
	637	! McGraw-Hill 1984.
[825]	638
[921]	639	maxnumeqmax = 0
[4873]	640	minnumeqmin = nlay_i+5
[825]	641
[8342]	642	DO ji = 1 , nidx
[5123]	643	zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji))
[921]	644	zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2)
	645	minnumeqmin = MIN(numeqmin(ji),minnumeqmin)
	646	maxnumeqmax = MAX(numeqmax(ji),maxnumeqmax)
	647	END DO
	648
[4870]	649	DO jk = minnumeqmin+1, maxnumeqmax
[8342]	650	DO ji = 1 , nidx
[4870]	651	numeq = min(max(numeqmin(ji)+1,jk),numeqmax(ji))
[5123]	652	zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3) / zdiagbis(ji,numeq-1)
	653	zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1)
[921]	654	END DO
	655	END DO
	656
[8342]	657	DO ji = 1 , nidx
[921]	658	! ice temperatures
[5123]	659	t_i_1d(ji,nlay_i) = zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji))
[921]	660	END DO
	661
[5180]	662	DO numeq = nlay_i + nlay_s, nlay_s + 2, -1
[8342]	663	DO ji = 1 , nidx
[4870]	664	jk = numeq - nlay_s - 1
[5123]	665	t_i_1d(ji,jk) = ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq)
[921]	666	END DO
	667	END DO
	668
[8342]	669	DO ji = 1 , nidx
[825]	670	! snow temperatures
[5123]	671	IF (ht_s_1d(ji) > 0._wp) &
	672	t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) &
	673	& / zdiagbis(ji,nlay_s+1) * MAX( 0.0, SIGN( 1.0, ht_s_1d(ji) ) )
[825]	674
	675	! surface temperature
[5123]	676	isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) )
[4872]	677	ztsubit(ji) = t_su_1d(ji)
[5123]	678	IF( t_su_1d(ji) < rt0 ) &
	679	t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) * &
	680	& ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))
[921]	681	END DO
	682	!
	683	!--------------------------------------------------------------------------
[5123]	684	! 11) Has the scheme converged ?, end of the iterative procedure \|
[921]	685	!--------------------------------------------------------------------------
	686	!
	687	! check that nowhere it has started to melt
[8313]	688	! zdti(ji) is a measure of error, it has to be under zdti_bnd
[8342]	689	DO ji = 1 , nidx
[5123]	690	t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , 190._wp )
[8313]	691	zdti (ji) = ABS( t_su_1d(ji) - ztsubit(ji) )
[921]	692	END DO
[825]	693
[4870]	694	DO jk = 1, nlay_s
[8342]	695	DO ji = 1 , nidx
[5123]	696	t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), 190._wp )
[8313]	697	zdti (ji) = MAX( zdti(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) )
[921]	698	END DO
	699	END DO
[825]	700
[4870]	701	DO jk = 1, nlay_i
[8342]	702	DO ji = 1 , nidx
[5123]	703	ztmelt_i = -tmut * s_i_1d(ji,jk) + rt0
	704	t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp )
[8313]	705	zdti (ji) = MAX( zdti(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) )
[921]	706	END DO
	707	END DO
[825]	708
[921]	709	! Compute spatial maximum over all errors
[2715]	710	! note that this could be optimized substantially by iterating only the non-converging points
[8313]	711	zdti_max = 0._wp
[8342]	712	DO ji = 1, nidx
[8313]	713	zdti_max = MAX( zdti_max, zdti(ji) )
[921]	714	END DO
[8313]	715	IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
[825]	716
	717	END DO ! End of the do while iterative procedure
	718
[8239]	719	! MV SIMIP 2016
	720	!--- Snow-ice interfacial temperature (diagnostic SIMIP)
[8342]	721	DO ji = 1, nidx
[8239]	722	zfac = 1. / MAX( epsi10 , rn_cdsn * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) )
[8369]	723	IF( zh_s(ji) >= 1.e-3 ) THEN
	724	t_si_1d(ji) = ( rn_cdsn * zh_i(ji) * t_s_1d(ji,1) + &
	725	& ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) * zfac
	726	ELSE
	727	t_si_1d(ji) = t_su_1d(ji)
	728	ENDIF
[8239]	729	END DO
	730	! END MV SIMIP 2016
	731
[7646]	732	IF( ln_limctl .AND. lwp ) THEN
[8313]	733	WRITE(numout,*) ' zdti_max : ', zdti_max
	734	WRITE(numout,*) ' iconv : ', iconv
[1055]	735	ENDIF
[825]	736
[921]	737	!
[2715]	738	!-------------------------------------------------------------------------!
[5123]	739	! 12) Fluxes at the interfaces !
[2715]	740	!-------------------------------------------------------------------------!
[8342]	741	DO ji = 1, nidx
[2715]	742	! ! surface ice conduction flux
[5123]	743	isnow(ji) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) )
	744	fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) &
	745	& - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji))
[2715]	746	! ! bottom ice conduction flux
[4872]	747	fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
[921]	748	END DO
[825]	749
[8239]	750	! MV SIMIP 2016
	751	!--- Conduction fluxes (positive downwards)
	752	diag_fc_bo_1d(:) = diag_fc_bo_1d(:) + fc_bo_i(:) * a_i_1d(:) / at_i_1d(:)
	753	diag_fc_su_1d(:) = diag_fc_su_1d(:) + fc_su(:) * a_i_1d(:) / at_i_1d(:)
	754	! END MV SIMIP 2016
	755
[4688]	756	! --- computes sea ice energy of melting compulsory for limthd_dh --- !
[8342]	757	CALL lim_thd_enmelt
[4688]	758
[5146]	759	! --- diagnose the change in non-solar flux due to surface temperature change --- !
[8313]	760	IF ( ln_dqnsice ) THEN
[8342]	761	DO ji = 1, nidx
[5474]	762	hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji)
	763	END DO
	764	END IF
[5123]	765
[4990]	766	! --- diag conservation imbalance on heat diffusion - PART 2 --- !
[8342]	767	DO ji = 1, nidx
[8325]	768	zdq(ji) = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + &
	769	& SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
[5123]	770	IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC
	771	zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice
	772	ELSE ! case T_su = 0degC
	773	zhfx_err(ji) = fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice
	774	ENDIF
[4990]	775	hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
[5167]	776
	777	! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
	778	hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
[4990]	779	END DO
	780
[5123]	781	!-----------------------------------------
	782	! Heat flux used to warm/cool ice in W.m-2
	783	!-----------------------------------------
[8342]	784	DO ji = 1, nidx
[5123]	785	IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC
	786	hfx_dif_1d(ji) = hfx_dif_1d(ji) + &
	787	& ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
	788	ELSE ! case T_su = 0degC
	789	hfx_dif_1d(ji) = hfx_dif_1d(ji) + &
	790	& ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
	791	ENDIF
[5167]	792	! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
	793	hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji)
	794	END DO
[4688]	795	!
	796	END SUBROUTINE lim_thd_dif
	797
[8342]	798	SUBROUTINE lim_thd_enmelt
[4688]	799	!!-----------------------------------------------------------------------
	800	!! * ROUTINE lim_thd_enmelt *
	801	!!
	802	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature
	803	!!
	804	!! ** Method : Formula (Bitz and Lipscomb, 1999)
	805	!!-------------------------------------------------------------------
	806	INTEGER :: ji, jk ! dummy loop indices
[4990]	807	REAL(wp) :: ztmelts ! local scalar
[4688]	808	!!-------------------------------------------------------------------
	809	!
	810	DO jk = 1, nlay_i ! Sea ice energy of melting
[8342]	811	DO ji = 1, nidx
[5202]	812	ztmelts = - tmut * s_i_1d(ji,jk) + rt0
	813	t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point
	814	! (sometimes dif scheme produces abnormally high temperatures)
[8325]	815	e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) &
[5202]	816	& + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) ) &
	817	& - rcp * ( ztmelts-rt0 ) )
[921]	818	END DO
[4688]	819	END DO
	820	DO jk = 1, nlay_s ! Snow energy of melting
[8342]	821	DO ji = 1, nidx
[8325]	822	e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
[921]	823	END DO
[4688]	824	END DO
[2715]	825	!
[4688]	826	END SUBROUTINE lim_thd_enmelt
[825]	827
	828	#else
[2715]	829	!!----------------------------------------------------------------------
	830	!! Dummy Module No LIM-3 sea-ice model
	831	!!----------------------------------------------------------------------
[825]	832	CONTAINS
	833	SUBROUTINE lim_thd_dif ! Empty routine
	834	END SUBROUTINE lim_thd_dif
	835	#endif
[2528]	836	!!======================================================================
[921]	837	END MODULE limthd_dif

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

Download in other formats: