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_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 7910

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

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

Download in other formats: