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

Context Navigation

source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 @ 5048

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

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

Download in other formats: