Context Navigation

icethd_zdf_bl99.F90 @ 14103

Last change on this file since 14103 was 14072, checked in by laurent, 3 years ago
Merging branch "2020/dev_r13648_ASINTER-04_laurent_bulk_ice", ticket #2369
Property svn:keywords set to `Id`
File size: 48.2 KB

Rev	Line
[8984]	1	MODULE icethd_zdf_BL99
	2	!!======================================================================
	3	!! * MODULE icethd_zdf_BL99 *
[14072]	4	!! sea-ice: vertical heat diffusion in sea ice (computation of temperatures)
[8984]	5	!!======================================================================
[9656]	6	!! History : ! 2003-02 (M. Vancoppenolle) original 1D code
[9604]	7	!! ! 2005-06 (M. Vancoppenolle) 3d version
	8	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
[8984]	9	!!----------------------------------------------------------------------
[9570]	10	#if defined key_si3
[8984]	11	!!----------------------------------------------------------------------
[9570]	12	!! 'key_si3' SI3 sea-ice model
[8984]	13	!!----------------------------------------------------------------------
	14	!! ice_thd_zdf_BL99 : vertical diffusion computation
	15	!!----------------------------------------------------------------------
	16	USE dom_oce ! ocean space and time domain
[14072]	17	USE phycst ! physical constants (ocean directory)
[8984]	18	USE ice ! sea-ice: variables
	19	USE ice1D ! sea-ice: thermodynamics variables
	20	USE icevar ! sea-ice: operations
	21	!
	22	USE in_out_manager ! I/O manager
	23	USE lib_mpp ! MPP library
	24	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	25
	26	IMPLICIT NONE
	27	PRIVATE
	28
	29	PUBLIC ice_thd_zdf_BL99 ! called by icethd_zdf
	30
	31	!!----------------------------------------------------------------------
[9598]	32	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	33	!! $Id$
[10068]	34	!! Software governed by the CeCILL license (see ./LICENSE)
[8984]	35	!!----------------------------------------------------------------------
	36	CONTAINS
	37
[10534]	38	SUBROUTINE ice_thd_zdf_BL99( k_cnd )
[8984]	39	!!-------------------------------------------------------------------
	40	!! * ROUTINE ice_thd_zdf_BL99 *
	41	!!
	42	!! ** Purpose : computes the time evolution of snow and sea-ice temperature
	43	!! profiles, using the original Bitz and Lipscomb (1999) algorithm
	44	!!
	45	!! ** Method : solves the heat equation diffusion with a Neumann boundary
[14072]	46	!! condition at the surface and a Dirichlet one at the bottom.
[8984]	47	!! Solar radiation is partially absorbed into the ice.
[14072]	48	!! The specific heat and thermal conductivities depend on ice
	49	!! salinity and temperature to take into account brine pocket
[8984]	50	!! melting. The numerical scheme is an iterative Crank-Nicolson
	51	!! on a non-uniform multilayer grid in the ice and snow system.
	52	!!
	53	!! The successive steps of this routine are
	54	!! 1. initialization of ice-snow layers thicknesses
	55	!! 2. Internal absorbed and transmitted radiation
	56	!! Then iterative procedure begins
	57	!! 3. Thermal conductivity
	58	!! 4. Kappa factors
	59	!! 5. specific heat in the ice
	60	!! 6. eta factors
	61	!! 7. surface flux computation
	62	!! 8. tridiagonal system terms
	63	!! 9. solving the tridiagonal system with Gauss elimination
	64	!! Iterative procedure ends according to a criterion on evolution
	65	!! of temperature
	66	!! 10. Fluxes at the interfaces
	67	!!
	68	!! ** Inputs / Ouputs : (global commons)
	69	!! surface temperature : t_su_1d
	70	!! ice/snow temperatures : t_i_1d, t_s_1d
	71	!! ice salinities : sz_i_1d
	72	!! number of layers in the ice/snow : nlay_i, nlay_s
	73	!! total ice/snow thickness : h_i_1d, h_s_1d
	74	!!-------------------------------------------------------------------
[10534]	75	INTEGER, INTENT(in) :: k_cnd ! conduction flux (off, on, emulated)
[8984]	76	!
	77	INTEGER :: ji, jk ! spatial loop index
	78	INTEGER :: jm ! current reference number of equation
	79	INTEGER :: jm_mint, jm_maxt
	80	INTEGER :: iconv ! number of iterations in iterative procedure
	81	INTEGER :: iconv_max = 50 ! max number of iterations in iterative procedure
	82	!
	83	INTEGER, DIMENSION(jpij) :: jm_min ! reference number of top equation
	84	INTEGER, DIMENSION(jpij) :: jm_max ! reference number of bottom equation
[10425]	85
	86	LOGICAL, DIMENSION(jpij) :: l_T_converged ! true when T converges (per grid point)
[13472]	87	!
[8984]	88	REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system
	89	REAL(wp) :: zg1 = 2._wp !
	90	REAL(wp) :: zgamma = 18009._wp ! for specific heat
	91	REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13)
	92	REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity
[14072]	93	REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C
	94	REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature
	95	REAL(wp) :: zhs_ssl = 0.03_wp ! surface scattering layer in the snow
[13472]	96	REAL(wp) :: zhi_ssl = 0.10_wp ! surface scattering layer in the ice
	97	REAL(wp) :: zh_min = 1.e-3_wp ! minimum ice/snow thickness for conduction
[9935]	98	REAL(wp) :: ztmelts ! ice melting temperature
[14072]	99	REAL(wp) :: zdti_max ! current maximal error on temperature
[8984]	100	REAL(wp) :: zcpi ! Ice specific heat
	101	REAL(wp) :: zhfx_err, zdq ! diag errors on heat
	102	!
[13472]	103	REAL(wp), DIMENSION(jpij) :: zraext_s ! extinction coefficient of radiation in the snow
[8984]	104	REAL(wp), DIMENSION(jpij) :: ztsub ! surface temperature at previous iteration
	105	REAL(wp), DIMENSION(jpij) :: zh_i, z1_h_i ! ice layer thickness
	106	REAL(wp), DIMENSION(jpij) :: zh_s, z1_h_s ! snow layer thickness
	107	REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface
	108	REAL(wp), DIMENSION(jpij) :: zfnet ! surface flux function
	109	REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function
	110	!
[14005]	111	REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice
	112	REAL(wp), DIMENSION(jpij,nlay_i) :: ztiold ! Old temperature in the ice
	113	REAL(wp), DIMENSION(jpij,nlay_s) :: ztsold ! Old temperature in the snow
	114	REAL(wp), DIMENSION(jpij,nlay_i) :: ztib ! Temporary temperature in the ice to check the convergence
	115	REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence
	116	REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity
	117	REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i_cp ! copy
	118	REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice
	119	REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice
	120	REAL(wp), DIMENSION(jpij,0:nlay_i) :: zkappa_i ! Kappa factor in the ice
	121	REAL(wp), DIMENSION(jpij,0:nlay_i) :: zeta_i ! Eta factor in the ice
	122	REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow
	123	REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow
	124	REAL(wp), DIMENSION(jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow
	125	REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow
	126	REAL(wp), DIMENSION(jpij) :: zkappa_comb ! Combined snow and ice surface conductivity
	127	REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat
	128	REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat
[14072]	129	REAL(wp), DIMENSION(jpij) :: za_s_fra ! ice fraction covered by snow
	130	REAL(wp), DIMENSION(jpij) :: isnow ! snow presence (1) or not (0)
	131	REAL(wp), DIMENSION(jpij) :: isnow_comb ! snow presence for met-office
[14005]	132	REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zindterm ! 'Ind'ependent term
	133	REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zindtbis ! Temporary 'ind'ependent term
	134	REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zdiagbis ! Temporary 'dia'gonal term
	135	REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1,3) :: ztrid ! Tridiagonal system terms
[8984]	136	!
	137	! Mono-category
	138	REAL(wp) :: zepsilon ! determines thres. above which computation of G(h) is done
	139	REAL(wp) :: zhe ! dummy factor
	140	REAL(wp) :: zcnd_i ! mean sea ice thermal conductivity
[14072]	141	!!------------------------------------------------------------------
[8984]	142
	143	! --- diag error on heat diffusion - PART 1 --- !
	144	DO ji = 1, npti
	145	zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + &
[14072]	146	& SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
[8984]	147	END DO
	148
[13472]	149	! calculate ice fraction covered by snow for radiation
	150	CALL ice_var_snwfra( h_s_1d(1:npti), za_s_fra(1:npti) )
[14072]	151
[8984]	152	!------------------
	153	! 1) Initialization
	154	!------------------
[13472]	155	!
	156	! extinction radiation in the snow
[14072]	157	IF ( nn_qtrice == 0 ) THEN ! constant
[13472]	158	zraext_s(1:npti) = rn_kappa_s
	159	ELSEIF( nn_qtrice == 1 ) THEN ! depends on melting/freezing conditions
	160	WHERE( t_su_1d(1:npti) < rt0 ) ; zraext_s(1:npti) = rn_kappa_sdry ! no surface melting
	161	ELSEWHERE ; zraext_s(1:npti) = rn_kappa_smlt ! surface melting
	162	END WHERE
	163	ENDIF
	164	!
	165	! thicknesses
[8984]	166	DO ji = 1, npti
[13472]	167	! ice thickness
[14072]	168	IF( h_i_1d(ji) > 0._wp ) THEN
[13472]	169	zh_i (ji) = MAX( zh_min , h_i_1d(ji) ) * r1_nlay_i ! set a minimum thickness for conduction
	170	z1_h_i(ji) = 1._wp / zh_i(ji) ! it must be very small
	171	ELSE
	172	zh_i (ji) = 0._wp
	173	z1_h_i(ji) = 0._wp
	174	ENDIF
	175	! snow thickness
	176	IF( h_s_1d(ji) > 0._wp ) THEN
	177	zh_s (ji) = MAX( zh_min , h_s_1d(ji) ) * r1_nlay_s ! set a minimum thickness for conduction
	178	z1_h_s(ji) = 1._wp / zh_s(ji) ! it must be very small
	179	isnow (ji) = 1._wp
	180	ELSE
	181	zh_s (ji) = 0._wp
	182	z1_h_s(ji) = 0._wp
	183	isnow (ji) = 0._wp
	184	ENDIF
	185	! for Met-Office
	186	IF( h_s_1d(ji) < zh_min ) THEN
	187	isnow_comb(ji) = h_s_1d(ji) / zh_min
	188	ELSE
	189	isnow_comb(ji) = 1._wp
	190	ENDIF
[8984]	191	END DO
[13472]	192	! clem: we should apply correction on snow thickness to take into account snow fraction
	193	! it must be a distribution, so it is a bit complicated
[8984]	194	!
	195	! Store initial temperatures and non solar heat fluxes
[10534]	196	IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
[8984]	197	ztsub (1:npti) = t_su_1d(1:npti) ! surface temperature at iteration n-1
	198	ztsuold (1:npti) = t_su_1d(1:npti) ! surface temperature initial value
	199	t_su_1d (1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err ) ! required to leave the choice between melting or not
[14072]	200	zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti) ! derivative of incoming nonsolar flux
[8984]	201	zqns_ice_b (1:npti) = qns_ice_1d(1:npti) ! store previous qns_ice_1d value
	202	!
	203	ENDIF
	204	!
	205	ztsold (1:npti,:) = t_s_1d(1:npti,:) ! Old snow temperature
	206	ztiold (1:npti,:) = t_i_1d(1:npti,:) ! Old ice temperature
	207
	208	!-------------
	209	! 2) Radiation
	210	!-------------
	211	! --- Transmission/absorption of solar radiation in the ice --- !
[9910]	212	zradtr_s(1:npti,0) = qtr_ice_top_1d(1:npti)
[8984]	213	DO jk = 1, nlay_s
	214	DO ji = 1, npti
	215	! ! radiation transmitted below the layer-th snow layer
[13472]	216	zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s(ji) * MAX( 0._wp, zh_s(ji) * REAL(jk) - zhs_ssl ) )
[8984]	217	! ! radiation absorbed by the layer-th snow layer
	218	zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
	219	END DO
	220	END DO
	221	!
[13472]	222	zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * za_s_fra(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - za_s_fra(1:npti) )
[14072]	223	DO jk = 1, nlay_i
[8984]	224	DO ji = 1, npti
	225	! ! radiation transmitted below the layer-th ice layer
[13472]	226	zradtr_i(ji,jk) = za_s_fra(ji) * zradtr_s(ji,nlay_s) & ! part covered by snow
	227	& * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zh_min ) ) &
	228	& + ( 1._wp - za_s_fra(ji) ) * qtr_ice_top_1d(ji) & ! part snow free
[14072]	229	& * EXP( - rn_kappa_i * MAX( 0._wp, zh_i(ji) * REAL(jk) - zhi_ssl ) )
[8984]	230	! ! radiation absorbed by the layer-th ice layer
	231	zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
	232	END DO
	233	END DO
	234	!
[9910]	235	qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti,nlay_i) ! record radiation transmitted below the ice
[8984]	236	!
[13472]	237	iconv = 0 ! number of iterations
[8984]	238	!
[10425]	239	l_T_converged(:) = .FALSE.
	240	! Convergence calculated until all sub-domain grid points have converged
	241	! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation)
	242	! but values are not taken into account (results independant of MPI partitioning)
	243	!
[10926]	244	! !============================!
[10425]	245	DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max ) ! Iterative procedure begins !
[10926]	246	! !============================!
[8984]	247	iconv = iconv + 1
	248	!
	249	ztib(1:npti,:) = t_i_1d(1:npti,:)
	250	ztsb(1:npti,:) = t_s_1d(1:npti,:)
	251	!
	252	!--------------------------------
	253	! 3) Sea ice thermal conductivity
	254	!--------------------------------
	255	IF( ln_cndi_U64 ) THEN !-- Untersteiner (1964) formula: k = k0 + beta.S/T
	256	!
	257	DO ji = 1, npti
[10425]	258	ztcond_i_cp(ji,0) = rcnd_i + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
	259	ztcond_i_cp(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )
[8984]	260	END DO
	261	DO jk = 1, nlay_i-1
	262	DO ji = 1, npti
[10425]	263	ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &
[13472]	264	& MIN( -epsi10, 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 )
[8984]	265	END DO
	266	END DO
	267	!
	268	ELSEIF( ln_cndi_P07 ) THEN !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
	269	!
	270	DO ji = 1, npti
[10425]	271	ztcond_i_cp(ji,0) = rcnd_i + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) &
[13472]	272	& - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
[10425]	273	ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) &
[13472]	274	& - 0.011_wp * ( t_bo_1d(ji) - rt0 )
[8984]	275	END DO
	276	DO jk = 1, nlay_i-1
	277	DO ji = 1, npti
[13472]	278	ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &
	279	& MIN( -epsi10, 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 ) &
	280	& - 0.011_wp * ( 0.5_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) ) - rt0 )
[8984]	281	END DO
	282	END DO
	283	!
	284	ENDIF
[10425]	285
	286	! Variable used after iterations
	287	! Value must be frozen after convergence for MPP independance reason
	288	DO ji = 1, npti
	289	IF ( .NOT. l_T_converged(ji) ) &
[14072]	290	ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) )
[10425]	291	END DO
[8984]	292	!
	293	!--- G(he) : enhancement of thermal conductivity in mono-category case
	294	! Computation of effective thermal conductivity G(h)
	295	! Used in mono-category case only to simulate an ITD implicitly
	296	! Fichefet and Morales Maqueda, JGR 1997
	297	zghe(1:npti) = 1._wp
	298	!
[10531]	299	IF( ln_virtual_itd ) THEN
[8984]	300	!
	301	zepsilon = 0.1_wp
	302	DO ji = 1, npti
	303	zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp ) ! Mean sea ice thermal conductivity
	304	zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i ) ! Effective thickness he (zhe)
	305	IF( zhe >= zepsilon * 0.5_wp * EXP(1._wp) ) &
	306	& zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) ) ! G(he)
	307	END DO
	308	!
[10531]	309	ENDIF
[8984]	310	!
	311	!-----------------
	312	! 4) kappa factors
	313	!-----------------
	314	!--- Snow
[10425]	315	! Variable used after iterations
	316	! Value must be frozen after convergence for MPP independance reason
[8984]	317	DO jk = 0, nlay_s-1
	318	DO ji = 1, npti
[10425]	319	IF ( .NOT. l_T_converged(ji) ) &
	320	zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
[8984]	321	END DO
	322	END DO
	323	DO ji = 1, npti ! Snow-ice interface
[13472]	324	IF ( .NOT. l_T_converged(ji) ) &
	325	zkappa_s(ji,nlay_s) = isnow(ji) * zghe(ji) * rn_cnd_s * ztcond_i(ji,0) &
	326	& / ( 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) )
[8984]	327	END DO
	328
	329	!--- Ice
[10425]	330	! Variable used after iterations
	331	! Value must be frozen after convergence for MPP independance reason
[8984]	332	DO jk = 0, nlay_i
	333	DO ji = 1, npti
[10425]	334	IF ( .NOT. l_T_converged(ji) ) &
	335	zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
[8984]	336	END DO
	337	END DO
	338	DO ji = 1, npti ! Snow-ice interface
[13472]	339	IF ( .NOT. l_T_converged(ji) ) THEN
	340	! Calculate combined surface snow and ice conductivity to pass through the coupler (met-office)
	341	zkappa_comb(ji) = isnow_comb(ji) * zkappa_s(ji,0) + ( 1._wp - isnow_comb(ji) ) * zkappa_i(ji,0)
	342	! If there is snow then use the same snow-ice interface conductivity for the top layer of ice
	343	IF( h_s_1d(ji) > 0._wp ) zkappa_i(ji,0) = zkappa_s(ji,nlay_s)
	344	ENDIF
[8984]	345	END DO
	346	!
	347	!--------------------------------------
	348	! 5) Sea ice specific heat, eta factors
	349	!--------------------------------------
	350	DO jk = 1, nlay_i
	351	DO ji = 1, npti
[9935]	352	zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
[13472]	353	zeta_i(ji,jk) = rDt_ice * r1_rhoi * z1_h_i(ji) / zcpi
[8984]	354	END DO
	355	END DO
	356
	357	DO jk = 1, nlay_s
	358	DO ji = 1, npti
[12489]	359	zeta_s(ji,jk) = rDt_ice * r1_rhos * r1_rcpi * z1_h_s(ji)
[8984]	360	END DO
	361	END DO
	362	!
	363	!----------------------------------------!
	364	! !
[10534]	365	! Conduction flux is off or emulated !
[8984]	366	! !
	367	!----------------------------------------!
	368	!
[10534]	369	IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
[8984]	370	!
	371	! ==> The original BL99 temperature computation is used
	372	! (with qsr_ice, qns_ice and dqns_ice as inputs)
	373	!
	374	!----------------------------
	375	! 6) surface flux computation
	376	!----------------------------
	377	! update of the non solar flux according to the update in T_su
	378	DO ji = 1, npti
[10425]	379	! Variable used after iterations
	380	! Value must be frozen after convergence for MPP independance reason
	381	IF ( .NOT. l_T_converged(ji) ) &
	382	qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
[8984]	383	END DO
	384
	385	DO ji = 1, npti
[9910]	386	zfnet(ji) = qsr_ice_1d(ji) - qtr_ice_top_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar
[8984]	387	END DO
	388	!
	389	!----------------------------
	390	! 7) tridiagonal system terms
	391	!----------------------------
	392	! layer denotes the number of the layer in the snow or in the ice
	393	! jm denotes the reference number of the equation in the tridiagonal
	394	! system, terms of tridiagonal system are indexed as following :
	395	! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
	396
	397	! ice interior terms (top equation has the same form as the others)
	398	ztrid (1:npti,:,:) = 0._wp
	399	zindterm(1:npti,:) = 0._wp
	400	zindtbis(1:npti,:) = 0._wp
	401	zdiagbis(1:npti,:) = 0._wp
	402
[14072]	403	DO jm = nlay_s + 2, nlay_s + nlay_i
[8984]	404	DO ji = 1, npti
	405	jk = jm - nlay_s - 1
	406	ztrid (ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
	407	ztrid (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
	408	ztrid (ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk)
	409	zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
	410	END DO
	411	END DO
	412
	413	jm = nlay_s + nlay_i + 1
	414	DO ji = 1, npti
	415	! ice bottom term
[14072]	416	ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1)
[8984]	417	ztrid (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
	418	ztrid (ji,jm,3) = 0._wp
	419	zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * &
[14072]	420	& ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )
[8984]	421	END DO
	422
	423	DO ji = 1, npti
	424	! !---------------------!
	425	IF( h_s_1d(ji) > 0._wp ) THEN ! snow-covered cells !
	426	! !---------------------!
	427	! snow interior terms (bottom equation has the same form as the others)
	428	DO jm = 3, nlay_s + 1
	429	jk = jm - 1
	430	ztrid (ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
	431	ztrid (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
	432	ztrid (ji,jm,3) = - zeta_s(ji,jk) * zkappa_s(ji,jk)
	433	zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
	434	END DO
[14072]	435
[8984]	436	! case of only one layer in the ice (ice equation is altered)
	437	IF( nlay_i == 1 ) THEN
	438	ztrid (ji,nlay_s+2,3) = 0._wp
[14072]	439	zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji)
[8984]	440	ENDIF
[14072]	441
[8984]	442	IF( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting
[14072]	443
[8984]	444	jm_min(ji) = 1
	445	jm_max(ji) = nlay_i + nlay_s + 1
[14072]	446
[8984]	447	! surface equation
	448	ztrid (ji,1,1) = 0._wp
	449	ztrid (ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
	450	ztrid (ji,1,3) = zg1s * zkappa_s(ji,0)
	451	zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
[14072]	452
[8984]	453	! first layer of snow equation
	454	ztrid (ji,2,1) = - zeta_s(ji,1) * zkappa_s(ji,0) * zg1s
	455	ztrid (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
	456	ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1)
	457	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
[14072]	458
[8984]	459	ELSE !-- case 2 : surface is melting
	460	!
	461	jm_min(ji) = 2
	462	jm_max(ji) = nlay_i + nlay_s + 1
[14072]	463
[8984]	464	! first layer of snow equation
	465	ztrid (ji,2,1) = 0._wp
	466	ztrid (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
[14072]	467	ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1)
	468	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )
[8984]	469	ENDIF
	470	! !---------------------!
	471	ELSE ! cells without snow !
	472	! !---------------------!
	473	!
	474	IF( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting
	475	!
	476	jm_min(ji) = nlay_s + 1
	477	jm_max(ji) = nlay_i + nlay_s + 1
[14072]	478
	479	! surface equation
[8984]	480	ztrid (ji,jm_min(ji),1) = 0._wp
[14072]	481	ztrid (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1
[8984]	482	ztrid (ji,jm_min(ji),3) = zkappa_i(ji,0) * zg1
	483	zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
[14072]	484
[8984]	485	! first layer of ice equation
	486	ztrid (ji,jm_min(ji)+1,1) = - zeta_i(ji,1) * zkappa_i(ji,0) * zg1
	487	ztrid (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
[14072]	488	ztrid (ji,jm_min(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1)
	489	zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)
	490
[8984]	491	! case of only one layer in the ice (surface & ice equations are altered)
	492	IF( nlay_i == 1 ) THEN
	493	ztrid (ji,jm_min(ji),1) = 0._wp
	494	ztrid (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2._wp
	495	ztrid (ji,jm_min(ji),3) = zkappa_i(ji,0) * 2._wp
	496	ztrid (ji,jm_min(ji)+1,1) = - zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp
	497	ztrid (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
	498	ztrid (ji,jm_min(ji)+1,3) = 0._wp
	499	zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji))
	500	ENDIF
[14072]	501
[8984]	502	ELSE !-- case 2 : surface is melting
[14072]	503
[8984]	504	jm_min(ji) = nlay_s + 2
	505	jm_max(ji) = nlay_i + nlay_s + 1
[14072]	506
[8984]	507	! first layer of ice equation
	508	ztrid (ji,jm_min(ji),1) = 0._wp
[14072]	509	ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
[8984]	510	ztrid (ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)
[14072]	511	zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji))
	512
[8984]	513	! case of only one layer in the ice (surface & ice equations are altered)
	514	IF( nlay_i == 1 ) THEN
	515	ztrid (ji,jm_min(ji),1) = 0._wp
	516	ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
	517	ztrid (ji,jm_min(ji),3) = 0._wp
	518	zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &
	519	& + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp
	520	ENDIF
[14072]	521
[8984]	522	ENDIF
	523	ENDIF
	524	!
	525	zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
	526	zdiagbis(ji,jm_min(ji)) = ztrid (ji,jm_min(ji),2)
	527	!
	528	END DO
	529	!
	530	!------------------------------
	531	! 8) tridiagonal system solving
	532	!------------------------------
	533	! Solve the tridiagonal system with Gauss elimination method.
	534	! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
[14005]	535	!!$ jm_maxt = 0
	536	!!$ jm_mint = nlay_i+5
	537	!!$ DO ji = 1, npti
	538	!!$ jm_mint = MIN(jm_min(ji),jm_mint)
	539	!!$ jm_maxt = MAX(jm_max(ji),jm_maxt)
	540	!!$ END DO
	541	!!$ !!clem SNWLAY => check why LIM1D does not get this loop. Is nlay_i+5 correct?
[14072]	542	!!$
[14005]	543	!!$ DO jk = jm_mint+1, jm_maxt
	544	!!$ DO ji = 1, npti
	545	!!$ jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
	546	!!$ zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1)
	547	!!$ zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1)
	548	!!$ END DO
	549	!!$ END DO
	550	! clem: maybe one should find a way to reverse this loop for mpi performance
[8984]	551	DO ji = 1, npti
[14005]	552	jm_mint = jm_min(ji)
	553	jm_maxt = jm_max(ji)
	554	DO jm = jm_mint+1, jm_maxt
[8984]	555	zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1)
	556	zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1)
	557	END DO
	558	END DO
	559
	560	! ice temperatures
	561	DO ji = 1, npti
[10425]	562	! Variable used after iterations
	563	! Value must be frozen after convergence for MPP independance reason
	564	IF ( .NOT. l_T_converged(ji) ) &
	565	t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
[8984]	566	END DO
	567
	568	DO jm = nlay_i + nlay_s, nlay_s + 2, -1
	569	DO ji = 1, npti
	570	jk = jm - nlay_s - 1
[10425]	571	IF ( .NOT. l_T_converged(ji) ) &
	572	t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
[8984]	573	END DO
	574	END DO
	575
[14072]	576	! snow temperatures
[8984]	577	DO ji = 1, npti
[10425]	578	! Variables used after iterations
	579	! Value must be frozen after convergence for MPP independance reason
[14005]	580	IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) &
	581	& t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
	582	END DO
	583	!!clem SNWLAY
	584	DO jm = nlay_s, 2, -1
	585	DO ji = 1, npti
	586	jk = jm - 1
	587	IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) &
	588	& t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm)
	589	END DO
	590	END DO
[14072]	591
[14005]	592	! surface temperature
	593	DO ji = 1, npti
	594	IF( .NOT. l_T_converged(ji) ) THEN
[10425]	595	ztsub(ji) = t_su_1d(ji)
	596	IF( t_su_1d(ji) < rt0 ) THEN
[14005]	597	t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * &
	598	& ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
[10425]	599	ENDIF
[8984]	600	ENDIF
	601	END DO
	602	!
	603	!--------------------------------------------------------------
	604	! 9) Has the scheme converged?, end of the iterative procedure
	605	!--------------------------------------------------------------
	606	! check that nowhere it has started to melt
	607	! zdti_max is a measure of error, it has to be under zdti_bnd
[10425]	608
[8984]	609	DO ji = 1, npti
	610
[10425]	611	zdti_max = 0._wp
[8984]	612
[10425]	613	IF ( .NOT. l_T_converged(ji) ) THEN
[13472]	614
[10425]	615	t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
	616	zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )
	617
[13472]	618	IF( h_s_1d(ji) > 0._wp ) THEN
	619	DO jk = 1, nlay_s
	620	t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
	621	zdti_max = MAX ( zdti_max , ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
	622	END DO
	623	ENDIF
[10425]	624
	625	DO jk = 1, nlay_i
	626	ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0
	627	t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
	628	zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
	629	END DO
[14072]	630
[13472]	631	! convergence test
	632	IF( ln_zdf_chkcvg ) THEN
	633	tice_cvgerr_1d(ji) = zdti_max
	634	tice_cvgstp_1d(ji) = REAL(iconv)
	635	ENDIF
[10425]	636
[13472]	637	IF( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE.
[10425]	638
	639	ENDIF
	640
[8984]	641	END DO
	642
	643	!----------------------------------------!
	644	! !
[10534]	645	! Conduction flux is on !
[8984]	646	! !
	647	!----------------------------------------!
	648	!
[10534]	649	ELSEIF( k_cnd == np_cnd_ON ) THEN
[8984]	650	!
	651	! ==> we use a modified BL99 solver with conduction flux (qcn_ice) as forcing term
	652	!
	653	!----------------------------
	654	! 7) tridiagonal system terms
	655	!----------------------------
	656	! layer denotes the number of the layer in the snow or in the ice
	657	! jm denotes the reference number of the equation in the tridiagonal
	658	! system, terms of tridiagonal system are indexed as following :
	659	! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
	660
	661	! ice interior terms (top equation has the same form as the others)
	662	ztrid (1:npti,:,:) = 0._wp
	663	zindterm(1:npti,:) = 0._wp
	664	zindtbis(1:npti,:) = 0._wp
	665	zdiagbis(1:npti,:) = 0._wp
	666
[14072]	667	DO jm = nlay_s + 2, nlay_s + nlay_i
[8984]	668	DO ji = 1, npti
	669	jk = jm - nlay_s - 1
	670	ztrid (ji,jm,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
	671	ztrid (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
	672	ztrid (ji,jm,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk)
	673	zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
	674	END DO
	675	ENDDO
	676
	677	jm = nlay_s + nlay_i + 1
	678	DO ji = 1, npti
	679	! ice bottom term
[14072]	680	ztrid (ji,jm,1) = - zeta_i(ji,nlay_i) * zkappa_i(ji,nlay_i-1)
[8984]	681	ztrid (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
	682	ztrid (ji,jm,3) = 0._wp
	683	zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) * &
[14072]	684	& ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )
[8984]	685	ENDDO
	686
	687	DO ji = 1, npti
	688	! !---------------------!
	689	IF( h_s_1d(ji) > 0._wp ) THEN ! snow-covered cells !
	690	! !---------------------!
	691	! snow interior terms (bottom equation has the same form as the others)
	692	DO jm = 3, nlay_s + 1
	693	jk = jm - 1
	694	ztrid (ji,jm,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
	695	ztrid (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
	696	ztrid (ji,jm,3) = - zeta_s(ji,jk) * zkappa_s(ji,jk)
	697	zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
	698	END DO
[14072]	699
[8984]	700	! case of only one layer in the ice (ice equation is altered)
	701	IF ( nlay_i == 1 ) THEN
	702	ztrid (ji,nlay_s+2,3) = 0._wp
[14072]	703	zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji)
[8984]	704	ENDIF
[14072]	705
[8984]	706	jm_min(ji) = 2
	707	jm_max(ji) = nlay_i + nlay_s + 1
[14072]	708
[8984]	709	! first layer of snow equation
	710	ztrid (ji,2,1) = 0._wp
	711	ztrid (ji,2,2) = 1._wp + zeta_s(ji,1) * zkappa_s(ji,1)
[14072]	712	ztrid (ji,2,3) = - zeta_s(ji,1) * zkappa_s(ji,1)
	713	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) )
	714
[8984]	715	! !---------------------!
	716	ELSE ! cells without snow !
	717	! !---------------------!
	718	jm_min(ji) = nlay_s + 2
	719	jm_max(ji) = nlay_i + nlay_s + 1
[14072]	720
[8984]	721	! first layer of ice equation
	722	ztrid (ji,jm_min(ji),1) = 0._wp
[14072]	723	ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)
[8984]	724	ztrid (ji,jm_min(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)
	725	zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )
[14072]	726
[8984]	727	! case of only one layer in the ice (surface & ice equations are altered)
	728	IF( nlay_i == 1 ) THEN
	729	ztrid (ji,jm_min(ji),1) = 0._wp
	730	ztrid (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)
	731	ztrid (ji,jm_min(ji),3) = 0._wp
	732	zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * &
	733	& ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) )
	734	ENDIF
[14072]	735
[8984]	736	ENDIF
	737	!
	738	zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
	739	zdiagbis(ji,jm_min(ji)) = ztrid (ji,jm_min(ji),2)
	740	!
	741	END DO
	742	!
	743	!------------------------------
	744	! 8) tridiagonal system solving
	745	!------------------------------
	746	! Solve the tridiagonal system with Gauss elimination method.
	747	! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
[14005]	748	!!$ jm_maxt = 0
	749	!!$ jm_mint = nlay_i+5
	750	!!$ DO ji = 1, npti
	751	!!$ jm_mint = MIN(jm_min(ji),jm_mint)
	752	!!$ jm_maxt = MAX(jm_max(ji),jm_maxt)
	753	!!$ END DO
[14072]	754	!!$
[14005]	755	!!$ DO jk = jm_mint+1, jm_maxt
	756	!!$ DO ji = 1, npti
	757	!!$ jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
	758	!!$ zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1)
	759	!!$ zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
	760	!!$ END DO
	761	!!$ END DO
	762	! clem: maybe one should find a way to reverse this loop for mpi performance
[8984]	763	DO ji = 1, npti
[14005]	764	jm_mint = jm_min(ji)
	765	jm_maxt = jm_max(ji)
	766	DO jm = jm_mint+1, jm_maxt
[8984]	767	zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1)
[14005]	768	zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1)
[8984]	769	END DO
	770	END DO
[14005]	771
[8984]	772	! ice temperatures
[10425]	773	DO ji = 1, npti
	774	! Variable used after iterations
	775	! Value must be frozen after convergence for MPP independance reason
	776	IF ( .NOT. l_T_converged(ji) ) &
	777	t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
[8984]	778	END DO
	779
	780	DO jm = nlay_i + nlay_s, nlay_s + 2, -1
	781	DO ji = 1, npti
[10425]	782	IF ( .NOT. l_T_converged(ji) ) THEN
	783	jk = jm - nlay_s - 1
	784	t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
	785	ENDIF
[8984]	786	END DO
	787	END DO
[14072]	788
	789	! snow temperatures
[8984]	790	DO ji = 1, npti
[14005]	791	! Variables used after iterations
[10425]	792	! Value must be frozen after convergence for MPP independance reason
[14005]	793	IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) &
	794	& t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
[8984]	795	END DO
[14005]	796	!!clem SNWLAY
	797	DO jm = nlay_s, 2, -1
	798	DO ji = 1, npti
	799	jk = jm - 1
	800	IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) &
	801	& t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm)
	802	END DO
	803	END DO
[8984]	804	!
	805	!--------------------------------------------------------------
	806	! 9) Has the scheme converged?, end of the iterative procedure
	807	!--------------------------------------------------------------
	808	! check that nowhere it has started to melt
	809	! zdti_max is a measure of error, it has to be under zdti_bnd
	810
[10425]	811	DO ji = 1, npti
	812
	813	zdti_max = 0._wp
	814
	815	IF ( .NOT. l_T_converged(ji) ) THEN
[13472]	816
	817	IF( h_s_1d(ji) > 0._wp ) THEN
	818	DO jk = 1, nlay_s
	819	t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
	820	zdti_max = MAX ( zdti_max , ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
	821	END DO
	822	ENDIF
	823
[10924]	824	DO jk = 1, nlay_i
[14072]	825	ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0
[10425]	826	t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
	827	zdti_max = MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
	828	END DO
	829
[13472]	830	! convergence test
	831	IF( ln_zdf_chkcvg ) THEN
	832	tice_cvgerr_1d(ji) = zdti_max
	833	tice_cvgstp_1d(ji) = REAL(iconv)
	834	ENDIF
[10425]	835
[13472]	836	IF( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE.
	837
[10425]	838	ENDIF
	839
[8984]	840	END DO
	841
[10534]	842	ENDIF ! k_cnd
[13472]	843
[8984]	844	END DO ! End of the do while iterative procedure
	845	!
	846	!-----------------------------
	847	! 10) Fluxes at the interfaces
	848	!-----------------------------
	849	!
[9916]	850	! --- calculate conduction fluxes (positive downward)
[12396]	851	! bottom ice conduction flux
[8984]	852	DO ji = 1, npti
[13472]	853	qcn_ice_bot_1d(ji) = - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )
[8984]	854	END DO
[12396]	855	! surface ice conduction flux
[10534]	856	IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
[8984]	857	!
	858	DO ji = 1, npti
[13472]	859	qcn_ice_top_1d(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) ) &
	860	& - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * ( t_i_1d(ji,1) - t_su_1d(ji) )
[8984]	861	END DO
	862	!
[10534]	863	ELSEIF( k_cnd == np_cnd_ON ) THEN
[8984]	864	!
	865	DO ji = 1, npti
[12396]	866	qcn_ice_top_1d(ji) = qcn_ice_1d(ji)
[8984]	867	END DO
	868	!
	869	ENDIF
[12396]	870	! surface ice temperature
	871	IF( k_cnd == np_cnd_ON .AND. ln_cndemulate ) THEN
	872	!
	873	DO ji = 1, npti
[13472]	874	t_su_1d(ji) = ( qcn_ice_top_1d(ji) + isnow(ji) * zkappa_s(ji,0) * zg1s * t_s_1d(ji,1) + &
	875	& ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * t_i_1d(ji,1) ) &
	876	& / MAX( epsi10, isnow(ji) * zkappa_s(ji,0) * zg1s + ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 )
[12396]	877	t_su_1d(ji) = MAX( MIN( t_su_1d(ji), rt0 ), rt0 - 100._wp ) ! cap t_su
	878	END DO
	879	!
	880	ENDIF
[8984]	881	!
[12396]	882	! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !
	883	!
	884	IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_EMU ) THEN
	885	!
	886	DO ji = 1, npti
[14072]	887	hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji)
[12396]	888	END DO
	889	!
	890	ENDIF
	891	!
[8984]	892	! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !
	893	!
[10534]	894	IF( k_cnd == np_cnd_OFF .OR. k_cnd == np_cnd_ON ) THEN
[14072]	895
	896	CALL ice_var_enthalpy
	897
[8984]	898	! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
	899	DO ji = 1, npti
	900	zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i + &
	901	& SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
[14072]	902
[10534]	903	IF( k_cnd == np_cnd_OFF ) THEN
[14072]	904
[8984]	905	IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC
[9916]	906	zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji) &
[12489]	907	& + zdq * r1_Dt_ice ) * a_i_1d(ji)
[8984]	908	ELSE ! case T_su = 0degC
[9916]	909	zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji) &
[12489]	910	& + zdq * r1_Dt_ice ) * a_i_1d(ji)
[8984]	911	ENDIF
[14072]	912
[10534]	913	ELSEIF( k_cnd == np_cnd_ON ) THEN
[14072]	914
[9916]	915	zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji) &
[12489]	916	& + zdq * r1_Dt_ice ) * a_i_1d(ji)
[14072]	917
[8984]	918	ENDIF
	919	!
	920	! total heat sink to be sent to the ocean
	921	hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
	922	!
[14072]	923	! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2
[12489]	924	hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_Dt_ice * a_i_1d(ji)
[8984]	925	!
	926	END DO
	927	!
	928	ENDIF
	929	!
[10534]	930	!--------------------------------------------------------------------
	931	! 11) reset inner snow and ice temperatures, update conduction fluxes
	932	!--------------------------------------------------------------------
[8984]	933	! effective conductivity and 1st layer temperature (needed by Met Office)
[13472]	934	! this is a conductivity at mid-layer, hence the factor 2
[8984]	935	DO ji = 1, npti
[13472]	936	IF( h_i_1d(ji) >= zhi_ssl ) THEN
	937	cnd_ice_1d(ji) = 2._wp * zkappa_comb(ji)
	938	!!cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
[8984]	939	ELSE
[13472]	940	cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / zhi_ssl ! cnd_ice is capped by: cond_i/zhi_ssl
[8984]	941	ENDIF
	942	t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
	943	END DO
	944	!
[10534]	945	IF( k_cnd == np_cnd_EMU ) THEN
[8984]	946	! Restore temperatures to their initial values
[9916]	947	t_s_1d (1:npti,:) = ztsold (1:npti,:)
	948	t_i_1d (1:npti,:) = ztiold (1:npti,:)
	949	qcn_ice_1d(1:npti) = qcn_ice_top_1d(1:npti)
[8984]	950	ENDIF
	951	!
[9916]	952	! --- SIMIP diagnostics
	953	!
[14072]	954	DO ji = 1, npti
[9916]	955	!--- Snow-ice interfacial temperature (diagnostic SIMIP)
[13472]	956	IF( h_s_1d(ji) >= zhs_ssl ) THEN
[14005]	957	t_si_1d(ji) = ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,nlay_s) &
	958	& + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) &
[13472]	959	& / ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i &
	960	& + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s )
[9916]	961	ELSE
	962	t_si_1d(ji) = t_su_1d(ji)
	963	ENDIF
	964	END DO
	965	!
[8984]	966	END SUBROUTINE ice_thd_zdf_BL99
	967
	968	#else
	969	!!----------------------------------------------------------------------
[9570]	970	!! Default option Dummy Module No SI3 sea-ice model
[8984]	971	!!----------------------------------------------------------------------
	972	#endif
	973
	974	!!======================================================================
	975	END MODULE icethd_zdf_BL99

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/trunk/src/ICE/icethd_zdf_bl99.F90 @ 14103

Download in other formats: