Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sisvat_qso.F @ 4163

Last change on this file since 4163 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

Line
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