Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

sisvat_qso.F @ 3817

Last change on this file since 3817 was 3792, checked in by evignon, 4 years ago
Ajout de INLANDSIS, nouvelle interface entre LMDZ et la neige de SISVAT Etienne, 04/01/2021
File size: 29.2 KB

Rev	Line
[3792]	1
	2
	3	subroutine SISVAT_qSo
	4	! #m0. (Wats_0,Wats_1,Wats_d)
	5
	6	C +------------------------------------------------------------------------+
	7	C \| MAR SISVAT_qSo 6-04-2001 MAR \|
	8	C \| SubRoutine SISVAT_qSo computes the Soil Water Balance \|
	9	C +------------------------------------------------------------------------+
	10	C \| \|
	11	C \| PARAMETERS: knonv: Total Number of columns = \|
	12	C \| ^^^^^^^^^^ = Total Number of continental grid boxes \|
	13	C \| X Number of Mosaic Cell per grid box \|
	14	C \| \|
	15	C \| INPUT: isnoSV = total Nb of Ice/Snow Layers \|
	16	C \| ^^^^^ isotSV = 0,...,11: Soil Type \|
	17	C \| 0: Water, Solid or Liquid \|
	18	C \| \|
	19	C \| INPUT: rhT_SV : SBL Top Air Density [kg/m3] \|
	20	C \| ^^^^^ drr_SV : Rain Intensity [kg/m2/s] \|
	21	C \| LSdzsv : Vertical Discretization Factor [-] \|
	22	C \| = 1. Soil \|
	23	C \| = 1000. Ocean \|
	24	C \| dt__SV : Time Step [s] \|
	25	C \| \|
	26	C \| Lx_H2O : Latent Heat of Vaporization/Sublimation [J/kg] \|
	27	C \| HLs_sv : Latent Heat Flux [W/m2] \|
	28	C \| \|
	29	C \| INPUT / eta_SV : Water Content [m3/m3] \|
	30	C \| OUTPUT: Khydsv : Soil Hydraulic Conductivity [m/s] \|
	31	C \| ^^^^^^ \|
	32	C \| \|
	33	C \| OUTPUT: RnofSV : RunOFF Intensity [kg/m2/s] \|
	34	C \| ^^^^^^ Wats_0 : Soil Water, before Forcing [mm] \|
	35	C \| Wats_1 : Soil Water, after Forcing [mm] \|
	36	C \| Wats_d : Soil Water Forcing [mm] \|
	37	C \| \|
	38	C \| Internal Variables: \|
	39	C \| ^^^^^^^^^^^^^^^^^^ \|
	40	C \| z_Bump : (Partly)Bumpy Layers Height [m] \|
	41	C \| z0Bump : Bumpy Layers Height [m] \|
	42	C \| dzBump : Lowest Bumpy Layer: [m] \|
	43	C \| etBump : Bumps Layer Averaged Humidity [m3/m3] \|
	44	C \| etaMid : Layer Interface's Humidity [m3/m3] \|
	45	C \| eta__f : Layer Humidity (Water Front)[m3/m3] \|
	46	C \| Dhyd_f : Soil Hydraulic Diffusivity (Water Front) [m2/s] \|
	47	C \| Dhydif : Soil Hydraulic Diffusivity [m2/s] \|
	48	C \| WgFlow : Water gravitational Flux [kg/m2/s] \|
	49	C \| Wg_MAX : Water MAXIMUM gravitational Flux [kg/m2/s] \|
	50	C \| SatRat : Water Saturation Flux [kg/m2/s] \|
	51	C \| WExces : Water Saturation Excess Flux [kg/m2/s] \|
	52	C \| Dhydtz : Dhydif * dt / dz [m] \|
	53	C \| FreeDr : Free Drainage Fraction [-] \|
	54	C \| Elem_A : A Diagonal Coefficient \|
	55	C \| Elem_C : C Diagonal Coefficient \|
	56	C \| Diag_A : A Diagonal \|
	57	C \| Diag_B : B Diagonal \|
	58	C \| Diag_C : C Diagonal \|
	59	C \| Term_D : Independant Term \|
	60	C \| Aux__P : P Auxiliary Variable \|
	61	C \| Aux__Q : Q Auxiliary Variable \|
	62	C \| \|
	63	C \| TUNING PARAMETER: \|
	64	C \| ^^^^^^^^^^^^^^^^ \|
	65	C \| z0soil : Soil Surface averaged Bumps Height [m] \|
	66	C \| \|
	67	C \| METHOD: NO Skin Surface Humidity \|
	68	C \| ^^^^^^ Semi-Implicit Crank Nicholson Scheme \|
	69	C \| (Partial) free Drainage, Water Bodies excepted (Lakes, Sea) \|
	70	C \| \|
	71
	72	C \| \|
	73	C \| # OPTIONS: #GF: Saturation Front \|
	74	C \| # ^^^^^^^ #GH: Saturation Front allows Horton Runoff \|
	75	C \| # #GA: Soil Humidity Geometric Average \|
	76	C \| # #BP: Parameterization of Terrain Bumps \|
	77	C \| \|
	78	C \| \|
	79	C +------------------------------------------------------------------------+
	80
	81
	82
	83
	84	C +--Global Variables
	85	C + ================
	86
	87	use VARphy
	88	use VAR_SV
	89	use VARdSV
	90	use VAR0SV
	91	use VARxSV
	92	use VARySV
	93
	94
	95	IMPLICIT NONE
	96
	97
	98	C +--OUTPUT
	99	C + ------
	100
	101	! Water (Mass) Budget
	102	! ~~~~~~~~~~~~~~~~~~~
	103	! #m0 real Wats_0(knonv) ! Soil Water, before forcing
	104	! #m0 real Wats_1(knonv) ! Soil Water, after forcing
	105	! #m0 real Wats_d(knonv) ! Soil Water forcing
	106
	107
	108	C +--Internal Variables
	109	C + ==================
	110
	111	integer isl ,jsl ,ist ,ikl !
	112	integer ikm ,ikp ,ik0 ,ik1 !
	113	integer ist__s,ist__w ! Soil/Water Body Identifier
	114	c #BP real z0soil ! Soil Surface Bumps Height [m]
	115	c #BP real z_Bump !(Partly)Bumpy Layers Height [m]
	116	c #BP real z0Bump ! Bumpy Layers Height [m]
	117	c #BP real dzBump ! Lowest Bumpy Layer:
	118
	119	c #BP real etBump(knonv) ! Bumps Layer Averaged Humidity
	120	real etaMid ! Layer Interface's Humidity
	121	real Dhydif ! Hydraulic Diffusivity [m2/s]
	122	real eta__f ! Water Front Soil Water Content
	123	real Khyd_f ! Water Front Hydraulic Conduct.
	124	real Khydav ! Hydraulic Conductivity [m/s]
	125	real WgFlow ! Water gravitat. Flux [kg/m2/s]
	126	real Wg_MAX ! Water MAX.grav. Flux [kg/m2/s]
	127	real SatRat ! Saturation Flux [kg/m2/s]
	128	real WExces ! Saturat. Excess Flux [kg/m2/s]
	129	real SoRnOF(knonv) ! Soil Run OFF
	130	real Dhydtz(knonv,-nsol:0) ! Dhydif * dt / dz [m]
	131	real Elem_A,Elem_B,Elem_C ! Diagonal Coefficients
	132	real Diag_A(knonv,-nsol:0) ! A Diagonal
	133	real Diag_B(knonv,-nsol:0) ! B Diagonal
	134	real Diag_C(knonv,-nsol:0) ! C Diagonal
	135	real Term_D(knonv,-nsol:0) ! Independant Term
	136	real Aux__P(knonv,-nsol:0) ! P Auxiliary Variable
	137	real Aux__Q(knonv,-nsol:0) ! Q Auxiliary Variable
	138	real etaaux(knonv,-nsol:-nsol+1) ! Soil Water Content [m3/m3]
	139	real FreeDr ! Free Drainage Fraction (actual)
	140	real FreeD0 ! Free Drainage Fraction (1=Full)
	141	real aKdtSV3( 0:nsot, 0:nkhy) ! Khyd=aeta+b: a dt
	142	real bKdtSV3( 0:nsot, 0:nkhy) ! Khyd=aeta+b: b dt
	143
	144	! Water (Mass) Budget
	145	! ~~~~~~~~~~~~~~~~~~~
	146	c #mw logical mwopen ! IO Switch
	147	c #mw common/Sm_qSo_L/mwopen !
	148	c #mw real hourwr,timewr !
	149	c #mw common/Sm_qSo_R/timewr !
	150	c #mw real Evapor(knonv) !
	151
	152
	153	C +--Internal DATA
	154	C + =============
	155
	156	c #BP data z0soil/0.020/ ! Soil Surface Bumps Height [m]
	157	data FreeD0/1.000/ ! Free Drainage Fraction (1=Full)
	158
	159	aKdtSV3=aKdtSV2*dt__SV
	160	bKdtSV3=bKdtSV2*dt__SV
	161
	162	! Water Budget (IN)
	163	! ==================
	164
	165	! #m0 DO ikl=1,knonv
	166	! #m0 Wats_0(ikl) = 0. ! OLD RunOFF Contrib.
	167	! #m0 Wats_d(ikl) = drr_SV(ikl) ! Water Surface Forc.
	168	! #m0 END DO
	169
	170	! #m0 isl= -nsol
	171	! #m0 DO ikl=1,knonv
	172	! #m0 Wats_0(ikl) = Wats_0(ikl)
	173	! #m0. + ro_Wat ( eta_SV(ikl,isl) dz78SV(isl)
	174	! #m0. + eta_SV(ikl,isl+1) dz_8SV(isl) ) LSdzsv(ikl)
	175	! #m0 END DO
	176
	177	! #m0 DO isl= -nsol+1,-1
	178	! #m0 DO ikl=1,knonv
	179	! #m0 Wats_0(ikl) = Wats_0(ikl)
	180	! #m0. + ro_Wat ( eta_SV(ikl,isl) dz34SV(isl)
	181	! #m0. +(eta_SV(ikl,isl-1)
	182	! #m0. +eta_SV(ikl,isl+1))dz_8SV(isl) ) LSdzsv(ikl)
	183	! #m0 END DO
	184	! #m0 END DO
	185
	186	! #m0 isl= 0
	187	! #m0 DO ikl=1,knonv
	188	! #m0 Wats_0(ikl) = Wats_0(ikl)
	189	! #m0. + ro_Wat ( eta_SV(ikl,isl) dz78SV(isl)
	190	! #m0. + eta_SV(ikl,isl-1) dz_8SV(isl) ) LSdzsv(ikl)
	191	! #m0 END DO
	192
	193
	194	C +--Gravitational Flow
	195	C + ==================
	196
	197	C +... METHOD: Surface Water Flux saturates successively the soil layers
	198	C + ^^^^^^ from up to below, but is limited by infiltration capacity.
	199	C + Hydraulic Conductivity again contributes after this step,
	200	C + not redundantly because of a constant (saturated) profile.
	201
	202	C +--Flux Limitor
	203	C + ^^^^^^^^^^^^^
	204	isl=0
	205	DO ikl=1,knonv
	206	ist = isotSV(ikl) ! Soil Type
	207	ist__s = min(ist, 1) ! 1 => Soil
	208	ist__w = 1 - ist__s ! 1 => Water Body
	209	Dhydif = s1__SV(ist)
	210	. *max(epsi,eta_SV(ikl,isl)) ! Hydraulic Diffusivity
	211	. **(bCHdSV(ist)+2.) ! DR97, Eqn.(3.36)
	212	Dhydif = ist__s * Dhydif !
	213	. + ist__w * vK_dSV ! Water Bodies
	214	C +
	215	Khydav = ist__s * Ks_dSV(ist) ! DR97 Assumption
	216	. + ist__w * vK_dSV ! Water Bodies
	217	C +
	218	Wg_MAX = ro_Wat *Dhydif ! MAXimum Infiltration
	219	. *(etadSV(ist)-eta_SV(ikl,isl)) ! Rate
	220	. /(dzAvSV(isl)*LSdzsv(ikl) ) !
	221	. + ro_Wat *Khydav !
	222
	223	C +--Surface Horton RunOFF
	224	C + ^^^^^^^^^^^^^^^^^^^^^
	225	SoRnOF(ikl) =
	226	. max(zero,drr_SV(ikl)-Wg_MAX)
	227	RuofSV(ikl,1) = RuofSV(ikl,1) + SoRnOF(ikl)
	228	drr_SV(ikl) = drr_SV(ikl)-SoRnOF(ikl)
	229	RuofSV(ikl,2) = RuofSV(ikl,2) +max(0.,drr_SV(ikl))
	230	END DO
	231
	232	c #GF DO isl=0,-nsol,-1
	233	c #GF DO ikl=1,knonv
	234	c #GF ist = isotSV(ikl) ! Soil Type
	235	c #GF ist__s = min(ist, 1) ! 1 => Soil
	236	c #GF ist__w = 1 - ist__s ! 1 => Water Body
	237
	238	C +--Water Diffusion
	239	C + ^^^^^^^^^^^^^^^
	240	c #GF Dhydif = s1__SV(ist)
	241	c #GF. *max(epsi,eta_SV(ikl,isl)) ! Hydraulic Diffusivity
	242	c #GF. **(bCHdSV(ist)+2.) ! DR97, Eqn.(3.36)
	243	c #GF Dhydif = ist__s * Dhydif !
	244	c #GF. + ist__w * vK_dSV ! Water Bodies
	245
	246	C +--Water Conduction (without Horton Runoff)
	247	C + ^^^^^^^^^^^^^^^^
	248	c #GF Khyd_f = Ks_dSV(ist)
	249	C +... Uses saturated K ==> Horton Runoff ~0 !
	250
	251	C +--Water Conduction (with Horton Runoff)
	252	C + ^^^^^^^^^^^^^^^^
	253	c #GH ik0 = nkhy *eta_SV(ikl,isl)
	254	c #GH. /etadSV(ist)
	255	c #GH eta__f = 1.
	256	c #GH. -aKdtSV3(ist,ik0)/(2. *dzAvSV(isl)
	257	c #GH. *LSdzsv(ikl))
	258	c #GH eta__f = max(eps_21,eta__f)
	259	c #GH eta__f = min(etadSV(ist),
	260	c #GH. eta_SV(ikl,isl) +
	261	c #GH. (aKdtSV3(ist,ik0) *eta_SV(ikl,isl)
	262	c #GH. +bKdtSV3(ist,ik0)) /(dzAvSV(isl)
	263	c #GH. *LSdzsv(ikl))
	264	c #GH. / eta__f )
	265	c #GH eta__f = .5*(eta_SV(ikl,isl)
	266	c #GH. +eta__f)
	267
	268	c #gh eta__f = eta_SV(ikl,isl)
	269
	270	c #GH ik0 = nkhy *eta__f
	271	c #GH. /etadSV(ist)
	272	c #GH Khyd_f =
	273	c #GH. (aKdtSV3(ist,ik0) *eta__f
	274	c #GH. +bKdtSV3(ist,ik0)) /dt__SV
	275
	276	c #GF Khydav = ist__s * Khyd_f ! DR97 Assumption
	277	c #GF. + ist__w * vK_dSV ! Water Bodies
	278
	279	C +--Gravitational Flow
	280	C + ^^^^^^^^^^^^^^^^^^
	281	c #GF Wg_MAX = ! MAXimum Infiltration
	282	c #GF. ro_Wat *Dhydif ! Rate
	283	c #GF. *(etadSV(ist)-eta_SV(ikl,isl)) !
	284	c #GF. /(dzAvSV(isl)*LSdzsv(ikl) ) !
	285	c #GF. + ro_Wat *Khydav !
	286	c #GF END DO
	287	c #GF END DO
	288	c #GF DO ikl=1,knonv
	289	c #GF SoRnOF(ikl) = SoRnOF(ikl) ! RunOFF Intensity
	290	c #GF. + drr_SV(ikl) ! [kg/m2/s]
	291	C +!!! Inclure la possibilite de creer une mare sur un bedrock impermeable
	292	c #GF drr_SV(ikl) = 0.
	293	c #GF END DO
	294
	295
	296	C +--Temperature Correction due to a changed Soil Energy Content
	297	C + ===========================================================
	298
	299	C +!!! Mettre en oeuvre le couplage humidit?-?nergie
	300
	301
	302	C +--Full Resolution of the Richard's Equation
	303	C + =========================================
	304
	305	C +... METHOD: Water content evolution results from water fluxes
	306	C + ^^^^^^ at the layer boundaries
	307	C + Conductivity is approximated by a piecewise linear profile.
	308	C + Semi-Implicit Crank-Nicholson scheme is used.
	309	C + (Bruen, 1997, Sensitivity of hydrological processes
	310	C + at the land-atmosphere interface.
	311	C + Proc. Royal Irish Academy, IGBP symposium
	312	C + on global change and the Irish Environment.
	313	C + Publ.: Maynooth)
	314
	315	C + - - - - - - - - isl+1/2 - - ^
	316	C + \|
	317	C + eta_SV(isl) --------------- isl ----- +--dz_dSV(isl) ^
	318	C + \| \|
	319	C + Dhydtz(isl) etaMid - - - - - - - - isl-1/2 - - v dzmiSV(isl)--+
	320	C + \|
	321	C + eta_SV(isl-1) --------------- isl-1 ----- v
	322
	323	C +--Transfert Coefficients
	324	C + ----------------------------
	325
	326	DO isl=-nsol+1,0
	327	DO ikl=1,knonv
	328	ist = isotSV(ikl) ! Soil Type
	329	ist__s = min(ist, 1) ! 1 => Soil
	330	ist__w = 1 - ist__s ! 1 => Water Body
	331	etaMid = (dz_dSV(isl) *eta_SV(ikl,isl-1) ! eta at layers
	332	. +dz_dSV(isl-1)*eta_SV(ikl,isl) ) ! interface
	333	. /(2.0* dzmiSV(isl)) ! LSdzsv implicit !
	334	c #GA etaMid = sqrt(dz_dSV(isl) *eta_SV(ikl,isl-1) ! Idem, geometric
	335	c #GA. dz_dSV(isl-1)eta_SV(ikl,isl) ) ! average
	336	c #GA. /(2.0* dzmiSV(isl)) ! (Vauclin&al.1979)
	337	Dhydif = s1__SV(ist) ! Hydraul.Diffusi.
	338	. (etaMid *( bCHdSV(ist)+2.)) ! DR97, Eqn.(3.36)
	339	Dhydtz(ikl,isl) = Dhydif*dt__SV !
	340	. /(dzmiSV(isl) !
	341	. *LSdzsv(ikl)) !
	342	Dhydtz(ikl,isl) = Dhydtz(ikl,isl) * ist__s ! Soil
	343	. +0.5dzmiSV(isl)LSdzsv(ikl) * ist__w ! Water bodies
	344
	345	END DO
	346	END DO
	347	isl=-nsol
	348	DO ikl=1,knonv
	349	Dhydtz(ikl,isl) = 0.0 !
	350	END DO
	351
	352
	353	C +--Tridiagonal Elimination: Set Up
	354	C + -------------------------------
	355
	356	C +--Soil/Snow Interior
	357	C + ^^^^^^^^^^^^^^^^^^
	358
	359	DO isl=0,-nsol,-1
	360	DO ikl=1,knonv
	361	ist = isotSV(ikl)
	362	eta_SV(ikl,isl) = max(epsi, eta_SV(ikl,isl))
	363	END DO
	364	END DO
	365
	366	DO isl=-nsol,-nsol+1
	367	DO ikl=1,knonv
	368	etaaux(ikl,isl) = eta_SV(ikl,isl)
	369	END DO
	370	END DO
	371
	372	DO isl=-nsol+1,-1
	373	DO ikl=1,knonv
	374	ist = isotSV(ikl)
	375	ikm = nkhy * eta_SV(ikl,isl-1) / etadSV(ist)
	376	ik0 = nkhy * eta_SV(ikl,isl) / etadSV(ist)
	377	ikp = nkhy * eta_SV(ikl,isl+1) / etadSV(ist)
	378
	379	if(ikm<0.or.ik0<0.or.ikp<0)then
	380	print *,"CRASH1 in sisvat_qso.f on pixel (i,j,n)",
	381	. ii__SV(ikl),jj__SV(ikl),nn__SV(ikl)
	382	print *,"decrease your time step or increase ntphys "//
	383	. "and ntdiff in time_steps.f"
	384	stop
	385	endif
	386
	387
	388	Elem_A = Dhydtz(ikl,isl)
	389	. - aKdtSV3(ist,ikm)* dziiSV(isl) *LSdzsv(ikl)
	390	Elem_B = - (Dhydtz(ikl,isl)
	391	. +Dhydtz(ikl,isl+1)
	392	. -aKdtSV3(ist,ik0)*(dziiSV(isl+1)
	393	. -dzi_SV(isl) )*LSdzsv(ikl))
	394	Elem_C = Dhydtz(ikl,isl+1)
	395	. + aKdtSV3(ist,ikp)* dzi_SV(isl+1)*LSdzsv(ikl)
	396	Diag_A(ikl,isl) = dz_8SV(isl) *LSdzsv(ikl)
	397	. -Implic * Elem_A
	398	Diag_B(ikl,isl) = dz34SV(isl) *LSdzsv(ikl)
	399	. -Implic * Elem_B
	400	Diag_C(ikl,isl) = dz_8SV(isl) *LSdzsv(ikl)
	401	. -Implic * Elem_C
	402
	403	Term_D(ikl,isl) = (dz_8SV(isl) *LSdzsv(ikl)
	404	. +Explic Elem_A )eta_SV(ikl,isl-1)
	405	. + (dz34SV(isl) *LSdzsv(ikl)
	406	. +Explic Elem_B )eta_SV(ikl,isl)
	407	. + (dz_8SV(isl) *LSdzsv(ikl)
	408	. +Explic Elem_C )eta_SV(ikl,isl+1)
	409	. + (bKdtSV3(ist,ikp)* dzi_SV(isl+1)
	410	. +bKdtSV3(ist,ik0)*(dziiSV(isl+1)
	411	. -dzi_SV(isl) )
	412	. -bKdtSV3(ist,ikm)* dziiSV(isl) )
	413	. * LSdzsv(ikl)
	414	END DO
	415	END DO
	416
	417	isl=-nsol
	418	DO ikl=1,knonv
	419	ist = isotSV(ikl)
	420	c # FreeDr = FreeD0 * min(ist,1)
	421	FreeDr = iWaFSV(ikl) * min(ist,1)
	422	ik0 = nkhy * eta_SV(ikl,isl ) / etadSV(ist)
	423	ikp = nkhy * eta_SV(ikl,isl+1) / etadSV(ist)
	424
	425	if(ik0<0.or.ikp<0)then
	426	print *,"CRASH2 in sisvat_qso.f on pixel (i,j,n)",
	427	. ii__SV(ikl),jj__SV(ikl),nn__SV(ikl)
	428	print *,"decrease your time step or increase ntphys "//
	429	. "and ntdiff in time_steps.f"
	430	stop
	431	endif
	432
	433	Elem_A = 0.
	434	Elem_B = - (Dhydtz(ikl,isl+1)
	435	. -aKdtSV3(ist,ik0)(dziiSV(isl+1)LSdzsv(ikl)
	436	. -FreeDr ))
	437	Elem_C = Dhydtz(ikl,isl+1)
	438	. + aKdtSV3(ist,ikp)* dzi_SV(isl+1)*LSdzsv(ikl)
	439	Diag_A(ikl,isl) = 0.
	440	Diag_B(ikl,isl) = dz78SV(isl) *LSdzsv(ikl)
	441	. -Implic *Elem_B
	442	Diag_C(ikl,isl) = dz_8SV(isl) *LSdzsv(ikl)
	443	. -Implic *Elem_C
	444
	445	Term_D(ikl,isl) = (dz78SV(isl) *LSdzsv(ikl)
	446	. +Explic Elem_B )eta_SV(ikl,isl)
	447	. + (dz_8SV(isl) *LSdzsv(ikl)
	448	. +Explic Elem_C )eta_SV(ikl,isl+1)
	449	. + (bKdtSV3(ist,ikp)* dzi_SV(isl+1)*LSdzsv(ikl)
	450	. +bKdtSV3(ist,ik0)(dziiSV(isl+1)LSdzsv(ikl)
	451	. -FreeDr ))
	452	END DO
	453
	454	isl=0
	455	DO ikl=1,knonv
	456	ist = isotSV(ikl)
	457	ikm = nkhy * eta_SV(ikl,isl-1) / etadSV(ist)
	458	ik0 = nkhy * eta_SV(ikl,isl) / etadSV(ist)
	459	Elem_A = Dhydtz(ikl,isl)
	460	. - aKdtSV3(ist,ikm)* dziiSV(isl)*LSdzsv(ikl)
	461	Elem_B = - (Dhydtz(ikl,isl)
	462	. +aKdtSV3(ist,ik0)* dzi_SV(isl)*LSdzsv(ikl))
	463	Elem_C = 0.
	464	Diag_A(ikl,isl) = dz_8SV(isl) *LSdzsv(ikl)
	465	. - Implic *Elem_A
	466	Diag_B(ikl,isl) = dz78SV(isl) *LSdzsv(ikl)
	467	. - Implic *Elem_B
	468	Diag_C(ikl,isl) = 0.
	469	C +
	470	Term_D(ikl,isl) = (dz_8SV(isl) *LSdzsv(ikl)
	471	. +Explic Elem_A )eta_SV(ikl,isl-1)
	472	. + (dz78SV(isl) *LSdzsv(ikl)
	473	. +Explic Elem_B )eta_SV(ikl,isl)
	474	. - (bKdtSV3(ist,ik0)* dzi_SV(isl)
	475	. +bKdtSV3(ist,ikm)* dziiSV(isl))*LSdzsv(ikl)
	476	. + dt__SV (HLs_sv(ikl) (1-min(1,isnoSV(ikl)))
	477	. / (ro_Wat dz_dSV(0) Lx_H2O(ikl))
	478	cXF bug 17/05/2017
	479	. +drr_SV(ikl))/ro_Wat
	480	END DO
	481
	482	DO ikl=1,knonv
	483	drr_SV(ikl)=0. ! drr is included in the 1st soil layer
	484	ENDDO
	485
	486	C +
	487	C +
	488	C +--Tridiagonal Elimination
	489	C + =======================
	490	C +
	491	C +--Forward Sweep
	492	C + ^^^^^^^^^^^^^^
	493	DO ikl= 1,knonv
	494	Aux__P(ikl,-nsol) = Diag_B(ikl,-nsol)
	495	Aux__Q(ikl,-nsol) =-Diag_C(ikl,-nsol)/Aux__P(ikl,-nsol)
	496	END DO
	497	C +
	498	DO isl=-nsol+1,0
	499	DO ikl= 1,knonv
	500	Aux__P(ikl,isl) = Diag_A(ikl,isl) *Aux__Q(ikl,isl-1)
	501	. +Diag_B(ikl,isl)
	502	Aux__Q(ikl,isl) =-Diag_C(ikl,isl) /Aux__P(ikl,isl)
	503	END DO
	504	END DO
	505	C +
	506	DO ikl= 1,knonv
	507	eta_SV(ikl,-nsol) = Term_D(ikl,-nsol)/Aux__P(ikl,-nsol)
	508	END DO
	509	C +
	510	DO isl=-nsol+1,0
	511	DO ikl= 1,knonv
	512	eta_SV(ikl,isl) =(Term_D(ikl,isl)
	513	. -Diag_A(ikl,isl) *eta_SV(ikl,isl-1))
	514	. /Aux__P(ikl,isl)
	515	END DO
	516	END DO
	517
	518	C +--Backward Sweep
	519	C + ^^^^^^^^^^^^^^
	520	DO isl=-1,-nsol,-1
	521	DO ikl= 1,knonv
	522	eta_SV(ikl,isl) = Aux__Q(ikl,isl) *eta_SV(ikl,isl+1)
	523	. +eta_SV(ikl,isl)
	524	END DO
	525	END DO
	526
	527
	528	C +--Horton RunOFF Intensity
	529	C + =======================
	530
	531	DO isl=0,-nsol,-1
	532	DO ikl=1,knonv
	533	ist = isotSV(ikl) ! Soil Type
	534	SatRat = (eta_SV(ikl,isl)-etadSV(ist)) ! OverSaturation Rate
	535	. ro_Wat dzAvSV(isl) !
	536	. *LSdzsv(ikl) !
	537	. /dt__SV !
	538	SoRnOF(ikl) = SoRnOF(ikl) !
	539	. + max(zero,SatRat) !
	540	RuofSV(ikl,3) = RuofSV(ikl,3) +
	541	. + max(zero,SatRat)
	542	eta_SV(ikl,isl) = max(epsi !
	543	c #ED. +etamSV(isotSV(ikl))!
	544	. ,eta_SV(ikl,isl)) !
	545	eta_SV(ikl,isl) = min(eta_SV(ikl,isl) !
	546	. ,etadSV(ist) ) !
	547	END DO
	548	END DO
	549
	550	C +--IO, for Verification
	551	C + ~~~~~~~~~~~~~~~~~~~~
	552	c #WR write(6,6010)
	553	6010 format(/,1x)
	554	DO isl= 0,-nsol,-1
	555	DO ikl= 1,knonv
	556	ist = isotSV(ikl)
	557	ikp = nkhy * eta_SV(ikl,isl) /etadSV(ist)
	558	Khydsv(ikl,isl) =(aKdtSV3(ist,ikp) *eta_SV(ikl,isl)
	559	. +bKdtSV3(ist,ikp)) *2.0/dt__SV
	560	c #WR write(6,6011) ikl,isl,eta_SV(ikl,isl)*1.e3,
	561	c #WR. ikp, aKdtSV3(ist,ikp),bKdtSV3(ist,ikp),
	562	c #WR. Khydsv(ikl,isl)
	563	6011 format(2i3,f8.1,i3,3e12.3)
	564	END DO
	565	END DO
	566
	567
	568	C +--Additional RunOFF Intensity
	569	C + ===========================
	570
	571	DO ikl=1,knonv
	572	ist = isotSV(ikl)
	573	ik0 = nkhy * etaaux(ikl,-nsol ) /etadSV(ist)
	574	c # FreeDr = FreeD0 * min(ist,1)
	575	FreeDr = iWaFSV(ikl) * min(ist,1)
	576	SoRnOF(ikl) = SoRnOF(ikl)
	577	. + (aKdtSV3(ist,ik0)(etaaux(ikl,-nsol)Explic
	578	. +eta_SV(ikl,-nsol)*Implic)
	579	. + bKdtSV3(ist,ik0) )
	580	. * FreeDr *ro_Wat /dt__SV
	581	RuofSV(ikl,3) = RuofSV(ikl,3)
	582	. + (aKdtSV3(ist,ik0)(etaaux(ikl,-nsol)Explic
	583	. +eta_SV(ikl,-nsol)*Implic)
	584	. + bKdtSV3(ist,ik0) )
	585	. * FreeDr *ro_Wat /dt__SV
	586
	587	C +--Full Run OFF: Update
	588	C + ~~~~~~~~~~~~~~~~~~~~
	589	RnofSV(ikl) = RnofSV(ikl) + SoRnOF(ikl)
	590	RuofSV(ikl,4) = RuofSV(ikl,4) + SoRnOF(ikl)
	591	END DO
	592
	593
	594	C +--Temperature Correction due to a changed Soil Energy Content
	595	C + ===========================================================
	596
	597	C +!!! Mettre en oeuvre le couplage humidit?-?nergie
	598
	599
	600	C +--Bumps/Asperites Treatment
	601	C + =========================
	602
	603	C +--Average over Bump Depth (z0soil)
	604	C + --------------------------------
	605
	606	c #BP z_Bump = 0.
	607	c #BP DO ikl=1,knonv
	608	c #BP etBump(ikl) = 0.
	609	c #BP END DO
	610	C +
	611	c #BP DO isl=0,-nsol,-1
	612	c #BP z0Bump = z_Bump
	613	c #BP z_Bump = z_Bump + dzAvSV(isl)
	614	c #BP IF (z_Bump.lt.z0soil) THEN
	615	c #BP DO ikl=1,knonv
	616	c #BP etBump(ikl) = etBump(ikl) + dzAvSV(isl) *eta_SV(ikl,isl)
	617	c #BP END DO
	618	c #BP END IF
	619	c #BP IF (z_Bump.gt.z0soil.AND.z0Bump.lt.z0soil) THEN
	620	c #BP DO ikl=1,knonv
	621	c #BP etBump(ikl) = etBump(ikl) + (z0soil-z0Bump)*eta_SV(ikl,isl)
	622	c #BP etBump(ikl) = etBump(ikl) / z0soil
	623	c #BP END DO
	624	c #BP END IF
	625	c #BP END DO
	626
	627
	628	C +--Correction
	629	C + ----------
	630
	631	c #BP z_Bump = 0.
	632	c #BP DO isl=0,-nsol,-1
	633	c #BP z0Bump = z_Bump
	634	c #BP z_Bump = z_Bump +dzAvSV(isl)
	635	c #BP IF (z_Bump.lt.z0soil) THEN
	636	c #BP DO ikl=1,knonv
	637	c #BP eta_SV(ikl,isl) = etBump(ikl)
	638	c #BP END DO
	639	c #BP END IF
	640	c #BP IF (z_Bump.gt.z0soil.AND.z0Bump.lt.z0soil) THEN
	641	c #BP dzBump = z_Bump - z0soil
	642	c #BP DO ikl=1,knonv
	643	c #BP eta_SV(ikl,isl) =(etBump(ikl) *(dzAvSV(isl)-dzBump)
	644	c #BP. + eta_SV(ikl,isl)* dzBump)
	645	c #BP. / dzAvSV(isl)
	646	c #BP END DO
	647	c #BP END IF
	648	c #BP END DO
	649
	650
	651	C +--Positive Definite
	652	C + =================
	653
	654	c #BP DO isl= 0,-nsol,-1
	655	c #BP DO ikl= 1,knonv
	656	c #BP eta_SV(ikl,isl) = max(epsi,eta_SV(ikl,isl))
	657	c #BP END DO
	658	c #BP END DO
	659
	660
	661	C +--Water Budget (OUT)
	662	C + ===================
	663
	664	! #m0 DO ikl=1,knonv
	665	! #m0 Wats_d(ikl) = Wats_d(ikl) !
	666	! #m0. + drr_SV(ikl) *zero ! Precipitation is
	667	C + \______________ already included
	668	! #m0. + HLs_sv(ikl)
	669	! #m0. *(1-min(isnoSV(ikl),1)) /Lx_H2O(ikl) ! Evaporation
	670	! #m0. - SoRnOF(ikl) ! Soil RunOFF Contrib.
	671	! #m0 Wats_1(ikl) = 0. !
	672	c #mw Evapor(ikl) = HLs_sv(ikl) *dt__SV !
	673	c #mw. *(1-min(isnoSV(ikl),1)) /Lx_H2O(ikl) !
	674	! #m0 END DO
	675
	676	! #m0 DO isl= -nsol,0
	677	! #m0 DO ikl=1,knonv
	678	! #m0 Wats_d(ikl) = Wats_d(ikl) !
	679	! #m0 END DO
	680	! #m0 END DO
	681	! #m0 DO ikl=1,knonv
	682	! #m0 Wats_d(ikl) = Wats_d(ikl) *dt__SV !
	683	! #m0 END DO
	684
	685	! #m0 isl= -nsol
	686	! #m0 DO ikl=1,knonv
	687	! #m0 Wats_1(ikl) = Wats_1(ikl)
	688	! #m0. + ro_Wat ( eta_SV(ikl,isl) dz78SV(isl)
	689	! #m0. + eta_SV(ikl,isl+1) dz_8SV(isl) ) LSdzsv(ikl)
	690	! #m0 END DO
	691
	692	! #m0 DO isl= -nsol+1,-1
	693	! #m0 DO ikl=1,knonv
	694	! #m0 Wats_1(ikl) = Wats_1(ikl)
	695	! #m0. + ro_Wat ( eta_SV(ikl,isl) dz34SV(isl)
	696	! #m0. +(eta_SV(ikl,isl-1)
	697	! #m0. +eta_SV(ikl,isl+1))dz_8SV(isl) ) LSdzsv(ikl)
	698	! #m0 END DO
	699	! #m0 END DO
	700
	701	! #m0 isl= 0
	702	! #m0 DO ikl=1,knonv
	703	! #m0 Wats_1(ikl) = Wats_1(ikl)
	704	! #m0. + ro_Wat ( eta_SV(ikl,isl) dz78SV(isl)
	705	! #m0. + eta_SV(ikl,isl-1) dz_8SV(isl) ) LSdzsv(ikl)
	706	! #m0 END DO
	707
	708
	709	return
	710	end

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

Download in other formats:

Original Format