Context Navigation

sisvat_qso.f90 @ 5300

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

Line
1
2
3	subroutine SISVAT_qSo
4	! #m0. (Wats_0,Wats_1,Wats_d)
5
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) + &
541	max(zero,SatRat)
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