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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

limthd_dif.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

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

Last change on this file since 4416 was 3808, checked in by gm, 11 years ago
dev_MERGE_2012: #997 and #898: in LIM3, no use of qla_ice in coupled mode and its capping in forced mode
Property svn:keywords set to `Id`
File size: 36.9 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
[825]	12	!!----------------------------------------------------------------------
[2528]	13	#if defined key_lim3
	14	!!----------------------------------------------------------------------
	15	!! 'key_lim3' LIM3 sea-ice model
	16	!!----------------------------------------------------------------------
[3625]	17	USE par_oce ! ocean parameters
	18	USE phycst ! physical constants (ocean directory)
	19	USE ice ! LIM-3 variables
	20	USE par_ice ! LIM-3 parameters
	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)
[921]	26
[825]	27	IMPLICIT NONE
	28	PRIVATE
	29
[2528]	30	PUBLIC lim_thd_dif ! called by lim_thd
[825]	31
[2715]	32	REAL(wp) :: epsi20 = 1e-20 ! constant values
	33	REAL(wp) :: epsi13 = 1e-13 ! constant values
[825]	34
	35	!!----------------------------------------------------------------------
[3625]	36	!! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)
[1156]	37	!! $Id$
[2591]	38	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	39	!!----------------------------------------------------------------------
	40	CONTAINS
	41
	42	SUBROUTINE lim_thd_dif( kideb , kiut , jl )
[921]	43	!!------------------------------------------------------------------
	44	!! * ROUTINE lim_thd_dif *
	45	!! ** Purpose :
	46	!! This routine determines the time evolution of snow and sea-ice
	47	!! temperature profiles.
	48	!! ** Method :
	49	!! This is done by solving the heat equation diffusion with
	50	!! a Neumann boundary condition at the surface and a Dirichlet one
	51	!! at the bottom. Solar radiation is partially absorbed into the ice.
	52	!! The specific heat and thermal conductivities depend on ice salinity
	53	!! and temperature to take into account brine pocket melting. The
	54	!! numerical
	55	!! scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid
	56	!! in the ice and snow system.
	57	!!
	58	!! The successive steps of this routine are
	59	!! 1. Thermal conductivity at the interfaces of the ice layers
	60	!! 2. Internal absorbed radiation
	61	!! 3. Scale factors due to non-uniform grid
	62	!! 4. Kappa factors
	63	!! Then iterative procedure begins
	64	!! 5. specific heat in the ice
	65	!! 6. eta factors
	66	!! 7. surface flux computation
	67	!! 8. tridiagonal system terms
	68	!! 9. solving the tridiagonal system with Gauss elimination
	69	!! Iterative procedure ends according to a criterion on evolution
	70	!! of temperature
	71	!!
	72	!! ** Arguments :
	73	!! kideb , kiut : Starting and ending points on which the
	74	!! the computation is applied
	75	!!
	76	!! ** Inputs / Ouputs : (global commons)
	77	!! surface temperature : t_su_b
	78	!! ice/snow temperatures : t_i_b, t_s_b
	79	!! ice salinities : s_i_b
	80	!! number of layers in the ice/snow: nlay_i, nlay_s
	81	!! profile of the ice/snow layers : z_i, z_s
	82	!! total ice/snow thickness : ht_i_b, ht_s_b
	83	!!
	84	!! ** External :
	85	!!
	86	!! ** References :
	87	!!
	88	!! ** History :
	89	!! (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium
	90	!! (06-2005) Martin Vancoppenolle, 3d version
	91	!! (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR)
	92	!! (04-2007) Energy conservation tested by M. Vancoppenolle
	93	!!------------------------------------------------------------------
[3294]	94	INTEGER , INTENT (in) :: kideb ! Start point on which the the computation is applied
	95	INTEGER , INTENT (in) :: kiut ! End point on which the the computation is applied
	96	INTEGER , INTENT (in) :: jl ! Category number
[825]	97
[921]	98	!! * Local variables
[3294]	99	INTEGER :: ji ! spatial loop index
	100	INTEGER :: ii, ij ! temporary dummy loop index
	101	INTEGER :: numeq ! current reference number of equation
	102	INTEGER :: layer ! vertical dummy loop index
	103	INTEGER :: nconv ! number of iterations in iterative procedure
	104	INTEGER :: minnumeqmin, maxnumeqmax
[3610]	105	INTEGER, DIMENSION(kiut) :: numeqmin ! reference number of top equation
	106	INTEGER, DIMENSION(kiut) :: numeqmax ! reference number of bottom equation
	107	INTEGER, DIMENSION(kiut) :: isnow ! switch for presence (1) or absence (0) of snow
[3294]	108	REAL(wp) :: zeps = 1.e-10_wp !
	109	REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system
	110	REAL(wp) :: zg1 = 2._wp !
	111	REAL(wp) :: zgamma = 18009._wp ! for specific heat
	112	REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13)
	113	REAL(wp) :: zraext_s = 1.e+8_wp ! extinction coefficient of radiation in the snow
	114	REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity
	115	REAL(wp) :: zht_smin = 1.e-4_wp ! minimum snow depth
[2715]	116	REAL(wp) :: ztmelt_i ! ice melting temperature
	117	REAL(wp) :: zerritmax ! current maximal error on temperature
[3610]	118	REAL(wp), DIMENSION(kiut) :: ztfs ! ice melting point
	119	REAL(wp), DIMENSION(kiut) :: ztsuold ! old surface temperature (before the iterative procedure )
	120	REAL(wp), DIMENSION(kiut) :: ztsuoldit ! surface temperature at previous iteration
	121	REAL(wp), DIMENSION(kiut) :: zh_i ! ice layer thickness
	122	REAL(wp), DIMENSION(kiut) :: zh_s ! snow layer thickness
	123	REAL(wp), DIMENSION(kiut) :: zfsw ! solar radiation absorbed at the surface
	124	REAL(wp), DIMENSION(kiut) :: zf ! surface flux function
	125	REAL(wp), DIMENSION(kiut) :: dzf ! derivative of the surface flux function
	126	REAL(wp), DIMENSION(kiut) :: zerrit ! current error on temperature
	127	REAL(wp), DIMENSION(kiut) :: zdifcase ! case of the equation resolution (1->4)
	128	REAL(wp), DIMENSION(kiut) :: zftrice ! solar radiation transmitted through the ice
	129	REAL(wp), DIMENSION(kiut) :: zihic, zhsu
	130	REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztcond_i ! Ice thermal conductivity
	131	REAL(wp), DIMENSION(kiut,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice
	132	REAL(wp), DIMENSION(kiut,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice
	133	REAL(wp), DIMENSION(kiut,0:nlay_i) :: zkappa_i ! Kappa factor in the ice
	134	REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztiold ! Old temperature in the ice
	135	REAL(wp), DIMENSION(kiut,0:nlay_i) :: zeta_i ! Eta factor in the ice
	136	REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztitemp ! Temporary temperature in the ice to check the convergence
	137	REAL(wp), DIMENSION(kiut,0:nlay_i) :: zspeche_i ! Ice specific heat
	138	REAL(wp), DIMENSION(kiut,0:nlay_i) :: z_i ! Vertical cotes of the layers in the ice
	139	REAL(wp), DIMENSION(kiut,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow
	140	REAL(wp), DIMENSION(kiut,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow
	141	REAL(wp), DIMENSION(kiut,0:nlay_s) :: zkappa_s ! Kappa factor in the snow
	142	REAL(wp), DIMENSION(kiut,0:nlay_s) :: zeta_s ! Eta factor in the snow
	143	REAL(wp), DIMENSION(kiut,0:nlay_s) :: ztstemp ! Temporary temperature in the snow to check the convergence
	144	REAL(wp), DIMENSION(kiut,0:nlay_s) :: ztsold ! Temporary temperature in the snow
	145	REAL(wp), DIMENSION(kiut,0:nlay_s) :: z_s ! Vertical cotes of the layers in the snow
	146	REAL(wp), DIMENSION(kiut,jkmax+2) :: zindterm ! Independent term
	147	REAL(wp), DIMENSION(kiut,jkmax+2) :: zindtbis ! temporary independent term
	148	REAL(wp), DIMENSION(kiut,jkmax+2) :: zdiagbis
	149	REAL(wp), DIMENSION(kiut,jkmax+2,3) :: ztrid ! tridiagonal system terms
[3625]	150	!!------------------------------------------------------------------
[3610]	151	!
[921]	152	!------------------------------------------------------------------------------!
	153	! 1) Initialization !
	154	!------------------------------------------------------------------------------!
	155	!
	156	DO ji = kideb , kiut
	157	! is there snow or not
[3625]	158	isnow(ji)= INT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) ) )
[921]	159	! surface temperature of fusion
[2715]	160	!!gm ??? ztfs(ji) = rtt !!!????
[921]	161	ztfs(ji) = isnow(ji) * rtt + (1.0-isnow(ji)) * rtt
	162	! layer thickness
[2715]	163	zh_i(ji) = ht_i_b(ji) / nlay_i
	164	zh_s(ji) = ht_s_b(ji) / nlay_s
[921]	165	END DO
[825]	166
[921]	167	!--------------------
	168	! Ice / snow layers
	169	!--------------------
[825]	170
[2715]	171	z_s(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st snow layer
	172	z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer
[825]	173
[2715]	174	DO layer = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer
[921]	175	DO ji = kideb , kiut
[2715]	176	z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / nlay_s
[921]	177	END DO
	178	END DO
[825]	179
[2715]	180	DO layer = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer
[921]	181	DO ji = kideb , kiut
[2715]	182	z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / nlay_i
[921]	183	END DO
	184	END DO
	185	!
	186	!------------------------------------------------------------------------------\|
	187	! 2) Radiations \|
	188	!------------------------------------------------------------------------------\|
	189	!
	190	!-------------------
	191	! Computation of i0
	192	!-------------------
	193	! i0 describes the fraction of solar radiation which does not contribute
	194	! to the surface energy budget but rather penetrates inside the ice.
	195	! We assume that no radiation is transmitted through the snow
	196	! If there is no no snow
	197	! zfsw = (1-i0).qsr_ice is absorbed at the surface
	198	! zftrice = io.qsr_ice is below the surface
	199	! fstbif = io.qsr_ice.exp(-k(h_i)) transmitted below the ice
[825]	200
[921]	201	DO ji = kideb , kiut
	202	! switches
[3625]	203	isnow(ji) = INT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) ) )
[921]	204	! hs > 0, isnow = 1
[2715]	205	zhsu (ji) = hnzst ! threshold for the computation of i0
	206	zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) )
[825]	207
[2715]	208	i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
[921]	209	!fr1_i0_1d = i0 for a thin ice surface
	210	!fr1_i0_2d = i0 for a thick ice surface
	211	! a function of the cloud cover
	212	!
	213	!i0(ji) = (1.0-FLOAT(isnow(ji)))3.0/(100ht_s_b(ji)+10.0)
	214	!formula used in Cice
	215	END DO
[825]	216
[921]	217	!-------------------------------------------------------
	218	! Solar radiation absorbed / transmitted at the surface
	219	! Derivative of the non solar flux
	220	!-------------------------------------------------------
	221	DO ji = kideb , kiut
[2715]	222	zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface
	223	zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer
	224	dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux
[921]	225	END DO
[825]	226
[921]	227	!---------------------------------------------------------
	228	! Transmission - absorption of solar radiation in the ice
	229	!---------------------------------------------------------
[825]	230
[2715]	231	DO ji = kideb, kiut ! snow initialization
	232	zradtr_s(ji,0) = zftrice(ji) ! radiation penetrating through snow
[921]	233	END DO
[825]	234
[2715]	235	DO layer = 1, nlay_s ! Radiation through snow
	236	DO ji = kideb, kiut
	237	! ! radiation transmitted below the layer-th snow layer
	238	zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) )
	239	! ! radiation absorbed by the layer-th snow layer
[921]	240	zradab_s(ji,layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)
	241	END DO
	242	END DO
[825]	243
[2715]	244	DO ji = kideb, kiut ! ice initialization
	245	zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) )
[921]	246	END DO
[825]	247
[2715]	248	DO layer = 1, nlay_i ! Radiation through ice
	249	DO ji = kideb, kiut
	250	! ! radiation transmitted below the layer-th ice layer
	251	zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) )
	252	! ! radiation absorbed by the layer-th ice layer
[921]	253	zradab_i(ji,layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)
	254	END DO
	255	END DO
[825]	256
[2715]	257	DO ji = kideb, kiut ! Radiation transmitted below the ice
	258	fstbif_1d(ji) = fstbif_1d(ji) + zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji)
[921]	259	END DO
[834]	260
[921]	261	! +++++
	262	! just to check energy conservation
[2715]	263	DO ji = kideb, kiut
[3625]	264	ii = MOD( npb(ji) - 1 , jpi ) + 1
	265	ij = ( npb(ji) - 1 ) / jpi + 1
[2715]	266	fstroc(ii,ij,jl) = zradtr_i(ji,nlay_i)
[921]	267	END DO
	268	! +++++
[825]	269
[921]	270	DO layer = 1, nlay_i
[2715]	271	DO ji = kideb, kiut
[921]	272	radab(ji,layer) = zradab_i(ji,layer)
	273	END DO
	274	END DO
[825]	275
[921]	276	!
	277	!------------------------------------------------------------------------------\|
	278	! 3) Iterative procedure begins \|
	279	!------------------------------------------------------------------------------\|
	280	!
[2715]	281	DO ji = kideb, kiut ! Old surface temperature
	282	ztsuold (ji) = t_su_b(ji) ! temperature at the beg of iter pr.
	283	ztsuoldit(ji) = t_su_b(ji) ! temperature at the previous iter
	284	t_su_b (ji) = MIN( t_su_b(ji), ztfs(ji)-0.00001 ) ! necessary
	285	zerrit (ji) = 1000._wp ! initial value of error
[921]	286	END DO
[825]	287
[2715]	288	DO layer = 1, nlay_s ! Old snow temperature
[921]	289	DO ji = kideb , kiut
[2715]	290	ztsold(ji,layer) = t_s_b(ji,layer)
[921]	291	END DO
	292	END DO
[825]	293
[2715]	294	DO layer = 1, nlay_i ! Old ice temperature
[921]	295	DO ji = kideb , kiut
[2715]	296	ztiold(ji,layer) = t_i_b(ji,layer)
[921]	297	END DO
	298	END DO
[825]	299
[2715]	300	nconv = 0 ! number of iterations
	301	zerritmax = 1000._wp ! maximal value of error on all points
[825]	302
[2715]	303	DO WHILE ( zerritmax > maxer_i_thd .AND. nconv < nconv_i_thd )
[921]	304	!
[2715]	305	nconv = nconv + 1
	306	!
[921]	307	!------------------------------------------------------------------------------\|
	308	! 4) Sea ice thermal conductivity \|
	309	!------------------------------------------------------------------------------\|
	310	!
[2715]	311	IF( thcon_i_swi == 0 ) THEN ! Untersteiner (1964) formula
[921]	312	DO ji = kideb , kiut
[2777]	313	ztcond_i(ji,0) = rcdic + zbeta*s_i_b(ji,1) / MIN(-zeps,t_i_b(ji,1)-rtt)
[921]	314	ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin)
	315	END DO
	316	DO layer = 1, nlay_i-1
	317	DO ji = kideb , kiut
[2777]	318	ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) / &
	319	MIN(-2.0_wp * zeps, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)
[921]	320	ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin)
	321	END DO
	322	END DO
	323	ENDIF
[825]	324
[2715]	325	IF( thcon_i_swi == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
[921]	326	DO ji = kideb , kiut
[2715]	327	ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -zeps, t_i_b(ji,1)-rtt ) &
	328	& - 0.011_wp * ( t_i_b(ji,1) - rtt )
	329	ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
[921]	330	END DO
[2715]	331	DO layer = 1, nlay_i-1
	332	DO ji = kideb , kiut
	333	ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) &
[2777]	334	& / MIN(-2.0_wp * zeps, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) &
[2715]	335	& - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )
	336	ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin )
	337	END DO
	338	END DO
	339	DO ji = kideb , kiut
	340	ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-zeps,t_bo_b(ji)-rtt) &
	341	& - 0.011_wp * ( t_bo_b(ji) - rtt )
	342	ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
	343	END DO
[921]	344	ENDIF
	345	!
	346	!------------------------------------------------------------------------------\|
	347	! 5) kappa factors \|
	348	!------------------------------------------------------------------------------\|
	349	!
	350	DO ji = kideb, kiut
[825]	351
[921]	352	!-- Snow kappa factors
	353	zkappa_s(ji,0) = rcdsn / MAX(zeps,zh_s(ji))
	354	zkappa_s(ji,nlay_s) = rcdsn / MAX(zeps,zh_s(ji))
	355	END DO
[825]	356
[921]	357	DO layer = 1, nlay_s-1
	358	DO ji = kideb , kiut
	359	zkappa_s(ji,layer) = 2.0 * rcdsn / &
	360	MAX(zeps,2.0*zh_s(ji))
	361	END DO
	362	END DO
[825]	363
[921]	364	DO layer = 1, nlay_i-1
	365	DO ji = kideb , kiut
	366	!-- Ice kappa factors
	367	zkappa_i(ji,layer) = 2.0*ztcond_i(ji,layer)/ &
	368	MAX(zeps,2.0*zh_i(ji))
	369	END DO
	370	END DO
[825]	371
[921]	372	DO ji = kideb , kiut
	373	zkappa_i(ji,0) = ztcond_i(ji,0)/MAX(zeps,zh_i(ji))
	374	zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(zeps,zh_i(ji))
	375	!-- Interface
	376	zkappa_s(ji,nlay_s) = 2.0rcdsnztcond_i(ji,0)/MAX(zeps, &
	377	(ztcond_i(ji,0)zh_s(ji) + rcdsnzh_i(ji)))
	378	zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*isnow(ji) &
	379	+ zkappa_i(ji,0)*(1.0-isnow(ji))
	380	END DO
	381	!
	382	!------------------------------------------------------------------------------\|
	383	! 6) Sea ice specific heat, eta factors \|
	384	!------------------------------------------------------------------------------\|
	385	!
	386	DO layer = 1, nlay_i
	387	DO ji = kideb , kiut
	388	ztitemp(ji,layer) = t_i_b(ji,layer)
	389	zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ &
	390	MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),zeps)
	391	zeta_i(ji,layer) = rdt_ice / MAX(rhoiczspeche_i(ji,layer)zh_i(ji), &
	392	zeps)
	393	END DO
	394	END DO
[825]	395
[921]	396	DO layer = 1, nlay_s
	397	DO ji = kideb , kiut
	398	ztstemp(ji,layer) = t_s_b(ji,layer)
	399	zeta_s(ji,layer) = rdt_ice / MAX(rhosncpiczh_s(ji),zeps)
	400	END DO
	401	END DO
	402	!
	403	!------------------------------------------------------------------------------\|
	404	! 7) surface flux computation \|
	405	!------------------------------------------------------------------------------\|
	406	!
	407	DO ji = kideb , kiut
[825]	408
[921]	409	! update of the non solar flux according to the update in T_su
	410	qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * &
	411	( t_su_b(ji) - ztsuoldit(ji) )
[825]	412
[921]	413	! update incoming flux
	414	zf(ji) = zfsw(ji) & ! net absorbed solar radiation
	415	+ qnsr_ice_1d(ji) ! non solar total flux
	416	! (LWup, LWdw, SH, LH)
[825]	417
[921]	418	END DO
[825]	419
[921]	420	!
	421	!------------------------------------------------------------------------------\|
	422	! 8) tridiagonal system terms \|
	423	!------------------------------------------------------------------------------\|
	424	!
	425	!!layer denotes the number of the layer in the snow or in the ice
	426	!!numeq denotes the reference number of the equation in the tridiagonal
	427	!!system, terms of tridiagonal system are indexed as following :
	428	!!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
[825]	429
[921]	430	!!ice interior terms (top equation has the same form as the others)
	431
	432	DO numeq=1,jkmax+2
	433	DO ji = kideb , kiut
	434	ztrid(ji,numeq,1) = 0.
	435	ztrid(ji,numeq,2) = 0.
	436	ztrid(ji,numeq,3) = 0.
	437	zindterm(ji,numeq)= 0.
	438	zindtbis(ji,numeq)= 0.
	439	zdiagbis(ji,numeq)= 0.
	440	ENDDO
	441	ENDDO
	442
	443	DO numeq = nlay_s + 2, nlay_s + nlay_i
	444	DO ji = kideb , kiut
	445	layer = numeq - nlay_s - 1
	446	ztrid(ji,numeq,1) = - zeta_i(ji,layer)*zkappa_i(ji,layer-1)
	447	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,layer)*(zkappa_i(ji,layer-1) + &
	448	zkappa_i(ji,layer))
	449	ztrid(ji,numeq,3) = - zeta_i(ji,layer)*zkappa_i(ji,layer)
	450	zindterm(ji,numeq) = ztiold(ji,layer) + zeta_i(ji,layer)* &
	451	zradab_i(ji,layer)
	452	END DO
	453	ENDDO
	454
	455	numeq = nlay_s + nlay_i + 1
	456	DO ji = kideb , kiut
[825]	457	!!ice bottom term
	458	ztrid(ji,numeq,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)
	459	ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i)( zkappa_i(ji,nlay_i)zg1 &
[921]	460	+ zkappa_i(ji,nlay_i-1) )
[825]	461	ztrid(ji,numeq,3) = 0.0
	462	zindterm(ji,numeq) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* &
[921]	463	( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 &
	464	* t_bo_b(ji) )
	465	ENDDO
[825]	466
	467
[921]	468	DO ji = kideb , kiut
[825]	469	IF ( ht_s_b(ji).gt.0.0 ) THEN
[921]	470	!
	471	!------------------------------------------------------------------------------\|
	472	! snow-covered cells \|
	473	!------------------------------------------------------------------------------\|
	474	!
	475	!!snow interior terms (bottom equation has the same form as the others)
	476	DO numeq = 3, nlay_s + 1
	477	layer = numeq - 1
	478	ztrid(ji,numeq,1) = - zeta_s(ji,layer)*zkappa_s(ji,layer-1)
	479	ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,layer)*( zkappa_s(ji,layer-1) + &
	480	zkappa_s(ji,layer) )
	481	ztrid(ji,numeq,3) = - zeta_s(ji,layer)*zkappa_s(ji,layer)
	482	zindterm(ji,numeq) = ztsold(ji,layer) + zeta_s(ji,layer)* &
	483	zradab_s(ji,layer)
	484	END DO
[825]	485
[921]	486	!!case of only one layer in the ice (ice equation is altered)
	487	IF ( nlay_i.eq.1 ) THEN
	488	ztrid(ji,nlay_s+2,3) = 0.0
	489	zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* &
	490	t_bo_b(ji)
	491	ENDIF
[834]	492
[921]	493	IF ( t_su_b(ji) .LT. rtt ) THEN
[825]	494
[921]	495	!------------------------------------------------------------------------------\|
	496	! case 1 : no surface melting - snow present \|
	497	!------------------------------------------------------------------------------\|
	498	zdifcase(ji) = 1.0
	499	numeqmin(ji) = 1
	500	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	501
[921]	502	!!surface equation
	503	ztrid(ji,1,1) = 0.0
	504	ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0)
	505	ztrid(ji,1,3) = zg1s*zkappa_s(ji,0)
	506	zindterm(ji,1) = dzf(ji)*t_su_b(ji) - zf(ji)
[825]	507
[921]	508	!!first layer of snow equation
	509	ztrid(ji,2,1) = - zkappa_s(ji,0)zg1szeta_s(ji,1)
	510	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)(zkappa_s(ji,1) + zkappa_s(ji,0)zg1s)
	511	ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1)
	512	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)
[825]	513
[921]	514	ELSE
	515	!
	516	!------------------------------------------------------------------------------\|
	517	! case 2 : surface is melting - snow present \|
	518	!------------------------------------------------------------------------------\|
	519	!
	520	zdifcase(ji) = 2.0
	521	numeqmin(ji) = 2
	522	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	523
[921]	524	!!first layer of snow equation
	525	ztrid(ji,2,1) = 0.0
	526	ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + &
	527	zkappa_s(ji,0) * zg1s )
	528	ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1)
	529	zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * &
	530	( zradab_s(ji,1) + &
	531	zkappa_s(ji,0) * zg1s * t_su_b(ji) )
	532	ENDIF
	533	ELSE
	534	!
	535	!------------------------------------------------------------------------------\|
	536	! cells without snow \|
	537	!------------------------------------------------------------------------------\|
	538	!
	539	IF (t_su_b(ji) .LT. rtt) THEN
	540	!
	541	!------------------------------------------------------------------------------\|
	542	! case 3 : no surface melting - no snow \|
	543	!------------------------------------------------------------------------------\|
	544	!
	545	zdifcase(ji) = 3.0
	546	numeqmin(ji) = nlay_s + 1
	547	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	548
[921]	549	!!surface equation
	550	ztrid(ji,numeqmin(ji),1) = 0.0
	551	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1
	552	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1
	553	zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_b(ji) - zf(ji)
[825]	554
[921]	555	!!first layer of ice equation
	556	ztrid(ji,numeqmin(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
	557	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) &
	558	+ zkappa_i(ji,0) * zg1 )
	559	ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1)
	560	zindterm(ji,numeqmin(ji)+1)= ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)
[825]	561
[921]	562	!!case of only one layer in the ice (surface & ice equations are altered)
[825]	563
[921]	564	IF (nlay_i.eq.1) THEN
	565	ztrid(ji,numeqmin(ji),1) = 0.0
	566	ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*2.0
	567	ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*2.0
	568	ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0)2.0zeta_i(ji,1)
	569	ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,0)2.0 + &
	570	zkappa_i(ji,1))
	571	ztrid(ji,numeqmin(ji)+1,3) = 0.0
[825]	572
[921]	573	zindterm(ji,numeqmin(ji)+1) = ztiold(ji,1) + zeta_i(ji,1)* &
	574	( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji) )
	575	ENDIF
[825]	576
[921]	577	ELSE
[825]	578
[921]	579	!
	580	!------------------------------------------------------------------------------\|
	581	! case 4 : surface is melting - no snow \|
	582	!------------------------------------------------------------------------------\|
	583	!
	584	zdifcase(ji) = 4.0
	585	numeqmin(ji) = nlay_s + 2
	586	numeqmax(ji) = nlay_i + nlay_s + 1
[825]	587
[921]	588	!!first layer of ice equation
	589	ztrid(ji,numeqmin(ji),1) = 0.0
	590	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,1) + zkappa_i(ji,0) &
	591	zg1)
	592	ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1)
	593	zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &
	594	zkappa_i(ji,0) * zg1 * t_su_b(ji) )
[825]	595
[921]	596	!!case of only one layer in the ice (surface & ice equations are altered)
	597	IF (nlay_i.eq.1) THEN
	598	ztrid(ji,numeqmin(ji),1) = 0.0
	599	ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)(zkappa_i(ji,0)2.0 + &
	600	zkappa_i(ji,1))
	601	ztrid(ji,numeqmin(ji),3) = 0.0
	602	zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)* &
	603	(zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji)) &
	604	+ t_su_b(ji)zeta_i(ji,1)zkappa_i(ji,0)*2.0
	605	ENDIF
[825]	606
[921]	607	ENDIF
	608	ENDIF
[825]	609
[921]	610	END DO
[825]	611
[921]	612	!
	613	!------------------------------------------------------------------------------\|
	614	! 9) tridiagonal system solving \|
	615	!------------------------------------------------------------------------------\|
	616	!
[825]	617
[921]	618	! Solve the tridiagonal system with Gauss elimination method.
	619	! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,
	620	! McGraw-Hill 1984.
[825]	621
[921]	622	maxnumeqmax = 0
	623	minnumeqmin = jkmax+4
[825]	624
[921]	625	DO ji = kideb , kiut
	626	zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji))
	627	zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2)
	628	minnumeqmin = MIN(numeqmin(ji),minnumeqmin)
	629	maxnumeqmax = MAX(numeqmax(ji),maxnumeqmax)
	630	END DO
	631
	632	DO layer = minnumeqmin+1, maxnumeqmax
	633	DO ji = kideb , kiut
	634	numeq = min(max(numeqmin(ji)+1,layer),numeqmax(ji))
	635	zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* &
	636	ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1)
	637	zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* &
	638	zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1)
	639	END DO
	640	END DO
	641
	642	DO ji = kideb , kiut
	643	! ice temperatures
	644	t_i_b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))
	645	END DO
	646
	647	DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1
	648	DO ji = kideb , kiut
	649	layer = numeq - nlay_s - 1
	650	t_i_b(ji,layer) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &
	651	t_i_b(ji,layer+1))/zdiagbis(ji,numeq)
	652	END DO
	653	END DO
	654
	655	DO ji = kideb , kiut
[825]	656	! snow temperatures
	657	IF (ht_s_b(ji).GT.0) &
[921]	658	t_s_b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
	659	* t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) &
	660	* MAX(0.0,SIGN(1.0,ht_s_b(ji)-zeps))
[825]	661
	662	! surface temperature
[3625]	663	isnow(ji) = INT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) ) ) )
[2715]	664	ztsuoldit(ji) = t_su_b(ji)
[3625]	665	IF( t_su_b(ji) < ztfs(ji) ) &
[2715]	666	t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( isnow(ji)*t_s_b(ji,1) &
	667	& + (1.0-isnow(ji))*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))
[921]	668	END DO
	669	!
	670	!--------------------------------------------------------------------------
	671	! 10) Has the scheme converged ?, end of the iterative procedure \|
	672	!--------------------------------------------------------------------------
	673	!
	674	! check that nowhere it has started to melt
	675	! zerrit(ji) is a measure of error, it has to be under maxer_i_thd
	676	DO ji = kideb , kiut
[2715]	677	t_su_b(ji) = MAX( MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp )
	678	zerrit(ji) = ABS( t_su_b(ji) - ztsuoldit(ji) )
[921]	679	END DO
[825]	680
[921]	681	DO layer = 1, nlay_s
	682	DO ji = kideb , kiut
[2715]	683	ii = MOD( npb(ji) - 1, jpi ) + 1
	684	ij = ( npb(ji) - 1 ) / jpi + 1
	685	t_s_b(ji,layer) = MAX( MIN( t_s_b(ji,layer), rtt ), 190._wp )
	686	zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer)))
[921]	687	END DO
	688	END DO
[825]	689
[921]	690	DO layer = 1, nlay_i
	691	DO ji = kideb , kiut
[2715]	692	ztmelt_i = -tmut * s_i_b(ji,layer) + rtt
	693	t_i_b(ji,layer) = MAX(MIN(t_i_b(ji,layer),ztmelt_i),190.0)
	694	zerrit(ji) = MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer)))
[921]	695	END DO
	696	END DO
[825]	697
[921]	698	! Compute spatial maximum over all errors
[2715]	699	! note that this could be optimized substantially by iterating only the non-converging points
	700	zerritmax = 0._wp
	701	DO ji = kideb, kiut
	702	zerritmax = MAX( zerritmax, zerrit(ji) )
[921]	703	END DO
[2715]	704	IF( lk_mpp ) CALL mpp_max( zerritmax, kcom=ncomm_ice )
[825]	705
	706	END DO ! End of the do while iterative procedure
	707
[1055]	708	IF( ln_nicep ) THEN
	709	WRITE(numout,*) ' zerritmax : ', zerritmax
	710	WRITE(numout,*) ' nconv : ', nconv
	711	ENDIF
[825]	712
[921]	713	!
[2715]	714	!-------------------------------------------------------------------------!
	715	! 11) Fluxes at the interfaces !
	716	!-------------------------------------------------------------------------!
[921]	717	DO ji = kideb, kiut
[3808]	718	#if ! defined key_coupled
	719	! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)
	720	qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) )
	721	#endif
[2715]	722	! ! surface ice conduction flux
	723	isnow(ji) = INT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ) )
	724	fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji)) &
	725	& - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_b(ji,1) - t_su_b(ji))
	726	! ! bottom ice conduction flux
	727	fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
[921]	728	END DO
[825]	729
[921]	730	!-------------------------!
	731	! Heat conservation !
	732	!-------------------------!
[2715]	733	IF( con_i ) THEN
[921]	734	DO ji = kideb, kiut
	735	! Upper snow value
[2715]	736	fc_s(ji,0) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) )
[921]	737	! Bott. snow value
[2715]	738	fc_s(ji,1) = - isnow(ji)* zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) )
[921]	739	END DO
[2715]	740	DO ji = kideb, kiut ! Upper ice layer
[921]	741	fc_i(ji,0) = - isnow(ji) * & ! interface flux if there is snow
	742	( zkappa_i(ji,0) * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) &
	743	- ( 1.0 - isnow(ji) ) * ( zkappa_i(ji,0) * &
	744	zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not
	745	END DO
[2715]	746	DO layer = 1, nlay_i - 1 ! Internal ice layers
[921]	747	DO ji = kideb, kiut
[2715]	748	fc_i(ji,layer) = - zkappa_i(ji,layer) * ( t_i_b(ji,layer+1) - t_i_b(ji,layer) )
	749	ii = MOD( npb(ji) - 1, jpi ) + 1
	750	ij = ( npb(ji) - 1 ) / jpi + 1
[921]	751	END DO
	752	END DO
[2715]	753	DO ji = kideb, kiut ! Bottom ice layers
	754	fc_i(ji,nlay_i) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
[921]	755	END DO
	756	ENDIF
[2715]	757	!
[921]	758	END SUBROUTINE lim_thd_dif
[825]	759
	760	#else
[2715]	761	!!----------------------------------------------------------------------
	762	!! Dummy Module No LIM-3 sea-ice model
	763	!!----------------------------------------------------------------------
[825]	764	CONTAINS
	765	SUBROUTINE lim_thd_dif ! Empty routine
	766	END SUBROUTINE lim_thd_dif
	767	#endif
[2528]	768	!!======================================================================
[921]	769	END MODULE limthd_dif

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

Download in other formats: