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

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 5180

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

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

Download in other formats: