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 @ 5078

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

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

Download in other formats: