source: LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_qso.F @ 4744

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