Context Navigation

sisvat_qso.f90 @ 5441

Last change on this file since 5441 was 5277, checked in by abarral, 2 months ago
Fix: unary operator following arithmetic operator
File size: 27.9 KB

Rev	Line
[3792]	1
	2
[5246]	3	subroutine SISVAT_qSo
	4	! #m0. (Wats_0,Wats_1,Wats_d)
[3792]	5
[5246]	6	! +------------------------------------------------------------------------+
	7	! \| MAR SISVAT_qSo 6-04-2001 MAR \|
	8	! \| SubRoutine SISVAT_qSo computes the Soil Water Balance \|
	9	! +------------------------------------------------------------------------+
	10	! \| \|
	11	! \| PARAMETERS: knonv: Total Number of columns = \|
	12	! \| ^^^^^^^^^^ = Total Number of continental grid boxes \|
	13	! \| X Number of Mosaic Cell per grid box \|
	14	! \| \|
	15	! \| INPUT: isnoSV = total Nb of Ice/Snow Layers \|
	16	! \| ^^^^^ isotSV = 0,...,11: Soil Type \|
	17	! \| 0: Water, Solid or Liquid \|
	18	! \| \|
	19	! \| INPUT: rhT_SV : SBL Top Air Density [kg/m3] \|
	20	! \| ^^^^^ drr_SV : Rain Intensity [kg/m2/s] \|
	21	! \| LSdzsv : Vertical Discretization Factor [-] \|
	22	! \| = 1. Soil \|
	23	! \| = 1000. Ocean \|
	24	! \| dt__SV : Time Step [s] \|
	25	! \| \|
	26	! \| Lx_H2O : Latent Heat of Vaporization/Sublimation [J/kg] \|
	27	! \| HLs_sv : Latent Heat Flux [W/m2] \|
	28	! \| \|
	29	! \| INPUT / eta_SV : Water Content [m3/m3] \|
	30	! \| OUTPUT: Khydsv : Soil Hydraulic Conductivity [m/s] \|
	31	! \| ^^^^^^ \|
	32	! \| \|
	33	! \| OUTPUT: RnofSV : RunOFF Intensity [kg/m2/s] \|
	34	! \| ^^^^^^ Wats_0 : Soil Water, before Forcing [mm] \|
	35	! \| Wats_1 : Soil Water, after Forcing [mm] \|
	36	! \| Wats_d : Soil Water Forcing [mm] \|
	37	! \| \|
	38	! \| Internal Variables: \|
	39	! \| ^^^^^^^^^^^^^^^^^^ \|
	40	! \| z_Bump : (Partly)Bumpy Layers Height [m] \|
	41	! \| z0Bump : Bumpy Layers Height [m] \|
	42	! \| dzBump : Lowest Bumpy Layer: [m] \|
	43	! \| etBump : Bumps Layer Averaged Humidity [m3/m3] \|
	44	! \| etaMid : Layer Interface's Humidity [m3/m3] \|
	45	! \| eta__f : Layer Humidity (Water Front)[m3/m3] \|
	46	! \| Dhyd_f : Soil Hydraulic Diffusivity (Water Front) [m2/s] \|
	47	! \| Dhydif : Soil Hydraulic Diffusivity [m2/s] \|
	48	! \| WgFlow : Water gravitational Flux [kg/m2/s] \|
	49	! \| Wg_MAX : Water MAXIMUM gravitational Flux [kg/m2/s] \|
	50	! \| SatRat : Water Saturation Flux [kg/m2/s] \|
	51	! \| WExces : Water Saturation Excess Flux [kg/m2/s] \|
	52	! \| Dhydtz : Dhydif * dt / dz [m] \|
	53	! \| FreeDr : Free Drainage Fraction [-] \|
	54	! \| Elem_A : A Diagonal Coefficient \|
	55	! \| Elem_C : C Diagonal Coefficient \|
	56	! \| Diag_A : A Diagonal \|
	57	! \| Diag_B : B Diagonal \|
	58	! \| Diag_C : C Diagonal \|
	59	! \| Term_D : Independant Term \|
	60	! \| Aux__P : P Auxiliary Variable \|
	61	! \| Aux__Q : Q Auxiliary Variable \|
	62	! \| \|
	63	! \| TUNING PARAMETER: \|
	64	! \| ^^^^^^^^^^^^^^^^ \|
	65	! \| z0soil : Soil Surface averaged Bumps Height [m] \|
	66	! \| \|
	67	! \| METHOD: NO Skin Surface Humidity \|
	68	! \| ^^^^^^ Semi-Implicit Crank Nicholson Scheme \|
	69	! \| (Partial) free Drainage, Water Bodies excepted (Lakes, Sea) \|
	70	! \| \|
	71
	72	! \| \|
	73	! \| # OPTIONS: #GF: Saturation Front \|
	74	! \| # ^^^^^^^ #GH: Saturation Front allows Horton Runoff \|
	75	! \| # #GA: Soil Humidity Geometric Average \|
	76	! \| # #BP: Parameterization of Terrain Bumps \|
	77	! \| \|
	78	! \| \|
	79	! +------------------------------------------------------------------------+
	80
	81
	82
	83
	84	! +--Global Variables
	85	! + ================
	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	! +--OUTPUT
	99	! + ------
	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	! +--Internal Variables
	109	! + ==================
	110
	111	integer :: isl ,jsl ,ist ,ikl !
	112	integer :: ikm ,ikp ,ik0 ,ik1 !
	113	integer :: ist__s,ist__w ! Soil/Water Body Identifier
	114	! #BP real z0soil ! Soil Surface Bumps Height [m]
	115	! #BP real z_Bump !(Partly)Bumpy Layers Height [m]
	116	! #BP real z0Bump ! Bumpy Layers Height [m]
	117	! #BP real dzBump ! Lowest Bumpy Layer:
	118
	119	! #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	! #mw logical mwopen ! IO Switch
	147	! #mw common/Sm_qSo_L/mwopen !
	148	! #mw real hourwr,timewr !
	149	! #mw common/Sm_qSo_R/timewr !
	150	! #mw real Evapor(knonv) !
	151
	152
	153	! +--Internal DATA
	154	! + =============
	155
	156	! #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	! +--Gravitational Flow
	195	! + ==================
	196
	197	! +... METHOD: Surface Water Flux saturates successively the soil layers
	198	! + ^^^^^^ from up to below, but is limited by infiltration capacity.
	199	! + Hydraulic Conductivity again contributes after this step,
	200	! + not redundantly because of a constant (saturated) profile.
	201
	202	! +--Flux Limitor
	203	! + ^^^^^^^^^^^^^
	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	! +
	215	Khydav = ist__s * Ks_dSV(ist) & ! DR97 Assumption
	216	+ ist__w * vK_dSV ! Water Bodies
	217	! +
	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	! +--Surface Horton RunOFF
	224	! + ^^^^^^^^^^^^^^^^^^^^^
	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	! #GF DO isl=0,-nsol,-1
	233	! #GF DO ikl=1,knonv
	234	! #GF ist = isotSV(ikl) ! Soil Type
	235	! #GF ist__s = min(ist, 1) ! 1 => Soil
	236	! #GF ist__w = 1 - ist__s ! 1 => Water Body
	237
	238	! +--Water Diffusion
	239	! + ^^^^^^^^^^^^^^^
	240	! #GF Dhydif = s1__SV(ist)
	241	! #GF. *max(epsi,eta_SV(ikl,isl)) ! Hydraulic Diffusivity
	242	! #GF. **(bCHdSV(ist)+2.) ! DR97, Eqn.(3.36)
	243	! #GF Dhydif = ist__s * Dhydif !
	244	! #GF. + ist__w * vK_dSV ! Water Bodies
	245
	246	! +--Water Conduction (without Horton Runoff)
	247	! + ^^^^^^^^^^^^^^^^
	248	! #GF Khyd_f = Ks_dSV(ist)
	249	! +... Uses saturated K ==> Horton Runoff ~0 !
	250
	251	! +--Water Conduction (with Horton Runoff)
	252	! + ^^^^^^^^^^^^^^^^
	253	! #GH ik0 = nkhy *eta_SV(ikl,isl)
	254	! #GH. /etadSV(ist)
	255	! #GH eta__f = 1.
	256	! #GH. -aKdtSV3(ist,ik0)/(2. *dzAvSV(isl)
	257	! #GH. *LSdzsv(ikl))
	258	! #GH eta__f = max(eps_21,eta__f)
	259	! #GH eta__f = min(etadSV(ist),
	260	! #GH. eta_SV(ikl,isl) +
	261	! #GH. (aKdtSV3(ist,ik0) *eta_SV(ikl,isl)
	262	! #GH. +bKdtSV3(ist,ik0)) /(dzAvSV(isl)
	263	! #GH. *LSdzsv(ikl))
	264	! #GH. / eta__f )
	265	! #GH eta__f = .5*(eta_SV(ikl,isl)
	266	! #GH. +eta__f)
	267
	268	! #gh eta__f = eta_SV(ikl,isl)
	269
	270	! #GH ik0 = nkhy *eta__f
	271	! #GH. /etadSV(ist)
	272	! #GH Khyd_f =
	273	! #GH. (aKdtSV3(ist,ik0) *eta__f
	274	! #GH. +bKdtSV3(ist,ik0)) /dt__SV
	275
	276	! #GF Khydav = ist__s * Khyd_f ! DR97 Assumption
	277	! #GF. + ist__w * vK_dSV ! Water Bodies
	278
	279	! +--Gravitational Flow
	280	! + ^^^^^^^^^^^^^^^^^^
	281	! #GF Wg_MAX = ! MAXimum Infiltration
	282	! #GF. ro_Wat *Dhydif ! Rate
	283	! #GF. *(etadSV(ist)-eta_SV(ikl,isl)) !
	284	! #GF. /(dzAvSV(isl)*LSdzsv(ikl) ) !
	285	! #GF. + ro_Wat *Khydav !
	286	! #GF END DO
	287	! #GF END DO
	288	! #GF DO ikl=1,knonv
	289	! #GF SoRnOF(ikl) = SoRnOF(ikl) ! RunOFF Intensity
	290	! #GF. + drr_SV(ikl) ! [kg/m2/s]
	291	! +!!! Inclure la possibilite de creer une mare sur un bedrock impermeable
	292	! #GF drr_SV(ikl) = 0.
	293	! #GF END DO
	294
	295
	296	! +--Temperature Correction due to a changed Soil Energy Content
	297	! + ===========================================================
	298
	299	! +!!! Mettre en oeuvre le couplage humidit?-?nergie
	300
	301
	302	! +--Full Resolution of the Richard's Equation
	303	! + =========================================
	304
	305	! +... METHOD: Water content evolution results from water fluxes
	306	! + ^^^^^^ at the layer boundaries
	307	! + Conductivity is approximated by a piecewise linear profile.
	308	! + Semi-Implicit Crank-Nicholson scheme is used.
	309	! + (Bruen, 1997, Sensitivity of hydrological processes
	310	! + at the land-atmosphere interface.
	311	! + Proc. Royal Irish Academy, IGBP symposium
	312	! + on global change and the Irish Environment.
	313	! + Publ.: Maynooth)
	314
	315	! + - - - - - - - - isl+1/2 - - ^
	316	! + \|
	317	! + eta_SV(isl) --------------- isl ----- +--dz_dSV(isl) ^
	318	! + \| \|
	319	! + Dhydtz(isl) etaMid - - - - - - - - isl-1/2 - - v dzmiSV(isl)--+
	320	! + \|
	321	! + eta_SV(isl-1) --------------- isl-1 ----- v
	322
	323	! +--Transfert Coefficients
	324	! + ----------------------------
	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	! #GA etaMid = sqrt(dz_dSV(isl) *eta_SV(ikl,isl-1) ! Idem, geometric
	335	! #GA. dz_dSV(isl-1)eta_SV(ikl,isl) ) ! average
	336	! #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	! +--Tridiagonal Elimination: Set Up
	354	! + -------------------------------
	355
	356	! +--Soil/Snow Interior
	357	! + ^^^^^^^^^^^^^^^^^^
	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	! # 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	! +
	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	!XF 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	! +
	487	! +
	488	! +--Tridiagonal Elimination
	489	! + =======================
	490	! +
	491	! +--Forward Sweep
	492	! + ^^^^^^^^^^^^^^
	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	! +
	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	! +
	506	DO ikl= 1,knonv
	507	eta_SV(ikl,-nsol) = Term_D(ikl,-nsol)/Aux__P(ikl,-nsol)
	508	END DO
	509	! +
	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	! +--Backward Sweep
	519	! + ^^^^^^^^^^^^^^
	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	! +--Horton RunOFF Intensity
	529	! + =======================
	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) + &
[5277]	541	max(zero,SatRat)
[5246]	542	eta_SV(ikl,isl) = max(epsi & !
	543	! #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	! +--IO, for Verification
	551	! + ~~~~~~~~~~~~~~~~~~~~
	552	! #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	! #WR write(6,6011) ikl,isl,eta_SV(ikl,isl)*1.e3,
	561	! #WR. ikp, aKdtSV3(ist,ikp),bKdtSV3(ist,ikp),
	562	! #WR. Khydsv(ikl,isl)
	563	6011 format(2i3,f8.1,i3,3e12.3)
	564	END DO
	565	END DO
	566
	567
	568	! +--Additional RunOFF Intensity
	569	! + ===========================
	570
	571	DO ikl=1,knonv
	572	ist = isotSV(ikl)
	573	ik0 = nkhy * etaaux(ikl,-nsol ) /etadSV(ist)
	574	! # 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	! +--Full Run OFF: Update
	588	! + ~~~~~~~~~~~~~~~~~~~~
	589	RnofSV(ikl) = RnofSV(ikl) + SoRnOF(ikl)
	590	RuofSV(ikl,4) = RuofSV(ikl,4) + SoRnOF(ikl)
	591	END DO
	592
	593
	594	! +--Temperature Correction due to a changed Soil Energy Content
	595	! + ===========================================================
	596
	597	! +!!! Mettre en oeuvre le couplage humidit?-?nergie
	598
	599
	600	! +--Bumps/Asperites Treatment
	601	! + =========================
	602
	603	! +--Average over Bump Depth (z0soil)
	604	! + --------------------------------
	605
	606	! #BP z_Bump = 0.
	607	! #BP DO ikl=1,knonv
	608	! #BP etBump(ikl) = 0.
	609	! #BP END DO
	610	! +
	611	! #BP DO isl=0,-nsol,-1
	612	! #BP z0Bump = z_Bump
	613	! #BP z_Bump = z_Bump + dzAvSV(isl)
	614	! #BP IF (z_Bump.lt.z0soil) THEN
	615	! #BP DO ikl=1,knonv
	616	! #BP etBump(ikl) = etBump(ikl) + dzAvSV(isl) *eta_SV(ikl,isl)
	617	! #BP END DO
	618	! #BP END IF
	619	! #BP IF (z_Bump.gt.z0soil.AND.z0Bump.lt.z0soil) THEN
	620	! #BP DO ikl=1,knonv
	621	! #BP etBump(ikl) = etBump(ikl) + (z0soil-z0Bump)*eta_SV(ikl,isl)
	622	! #BP etBump(ikl) = etBump(ikl) / z0soil
	623	! #BP END DO
	624	! #BP END IF
	625	! #BP END DO
	626
	627
	628	! +--Correction
	629	! + ----------
	630
	631	! #BP z_Bump = 0.
	632	! #BP DO isl=0,-nsol,-1
	633	! #BP z0Bump = z_Bump
	634	! #BP z_Bump = z_Bump +dzAvSV(isl)
	635	! #BP IF (z_Bump.lt.z0soil) THEN
	636	! #BP DO ikl=1,knonv
	637	! #BP eta_SV(ikl,isl) = etBump(ikl)
	638	! #BP END DO
	639	! #BP END IF
	640	! #BP IF (z_Bump.gt.z0soil.AND.z0Bump.lt.z0soil) THEN
	641	! #BP dzBump = z_Bump - z0soil
	642	! #BP DO ikl=1,knonv
	643	! #BP eta_SV(ikl,isl) =(etBump(ikl) *(dzAvSV(isl)-dzBump)
	644	! #BP. + eta_SV(ikl,isl)* dzBump)
	645	! #BP. / dzAvSV(isl)
	646	! #BP END DO
	647	! #BP END IF
	648	! #BP END DO
	649
	650
	651	! +--Positive Definite
	652	! + =================
	653
	654	! #BP DO isl= 0,-nsol,-1
	655	! #BP DO ikl= 1,knonv
	656	! #BP eta_SV(ikl,isl) = max(epsi,eta_SV(ikl,isl))
	657	! #BP END DO
	658	! #BP END DO
	659
	660
	661	! +--Water Budget (OUT)
	662	! + ===================
	663
	664	! #m0 DO ikl=1,knonv
	665	! #m0 Wats_d(ikl) = Wats_d(ikl) !
	666	! #m0. + drr_SV(ikl) *zero ! Precipitation is
	667	! + \______________ 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	! #mw Evapor(ikl) = HLs_sv(ikl) *dt__SV !
	673	! #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 subroutine sisvat_qso

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

Download in other formats:

Original Format