subroutine SISVAT_qSo ! #m0. (Wats_0,Wats_1,Wats_d) ! +------------------------------------------------------------------------+ ! | MAR SISVAT_qSo 6-04-2001 MAR | ! | SubRoutine SISVAT_qSo computes the Soil Water Balance | ! +------------------------------------------------------------------------+ ! | | ! | PARAMETERS: knonv: Total Number of columns = | ! | ^^^^^^^^^^ = Total Number of continental grid boxes | ! | X Number of Mosaic Cell per grid box | ! | | ! | INPUT: isnoSV = total Nb of Ice/Snow Layers | ! | ^^^^^ isotSV = 0,...,11: Soil Type | ! | 0: Water, Solid or Liquid | ! | | ! | INPUT: rhT_SV : SBL Top Air Density [kg/m3] | ! | ^^^^^ drr_SV : Rain Intensity [kg/m2/s] | ! | LSdzsv : Vertical Discretization Factor [-] | ! | = 1. Soil | ! | = 1000. Ocean | ! | dt__SV : Time Step [s] | ! | | ! | Lx_H2O : Latent Heat of Vaporization/Sublimation [J/kg] | ! | HLs_sv : Latent Heat Flux [W/m2] | ! | | ! | INPUT / eta_SV : Water Content [m3/m3] | ! | OUTPUT: Khydsv : Soil Hydraulic Conductivity [m/s] | ! | ^^^^^^ | ! | | ! | OUTPUT: RnofSV : RunOFF Intensity [kg/m2/s] | ! | ^^^^^^ Wats_0 : Soil Water, before Forcing [mm] | ! | Wats_1 : Soil Water, after Forcing [mm] | ! | Wats_d : Soil Water Forcing [mm] | ! | | ! | Internal Variables: | ! | ^^^^^^^^^^^^^^^^^^ | ! | z_Bump : (Partly)Bumpy Layers Height [m] | ! | z0Bump : Bumpy Layers Height [m] | ! | dzBump : Lowest Bumpy Layer: [m] | ! | etBump : Bumps Layer Averaged Humidity [m3/m3] | ! | etaMid : Layer Interface's Humidity [m3/m3] | ! | eta__f : Layer Humidity (Water Front)[m3/m3] | ! | Dhyd_f : Soil Hydraulic Diffusivity (Water Front) [m2/s] | ! | Dhydif : Soil Hydraulic Diffusivity [m2/s] | ! | WgFlow : Water gravitational Flux [kg/m2/s] | ! | Wg_MAX : Water MAXIMUM gravitational Flux [kg/m2/s] | ! | SatRat : Water Saturation Flux [kg/m2/s] | ! | WExces : Water Saturation Excess Flux [kg/m2/s] | ! | Dhydtz : Dhydif * dt / dz [m] | ! | FreeDr : Free Drainage Fraction [-] | ! | Elem_A : A Diagonal Coefficient | ! | Elem_C : C Diagonal Coefficient | ! | Diag_A : A Diagonal | ! | Diag_B : B Diagonal | ! | Diag_C : C Diagonal | ! | Term_D : Independant Term | ! | Aux__P : P Auxiliary Variable | ! | Aux__Q : Q Auxiliary Variable | ! | | ! | TUNING PARAMETER: | ! | ^^^^^^^^^^^^^^^^ | ! | z0soil : Soil Surface averaged Bumps Height [m] | ! | | ! | METHOD: NO Skin Surface Humidity | ! | ^^^^^^ Semi-Implicit Crank Nicholson Scheme | ! | (Partial) free Drainage, Water Bodies excepted (Lakes, Sea) | ! | | ! | | ! | # OPTIONS: #GF: Saturation Front | ! | # ^^^^^^^ #GH: Saturation Front allows Horton Runoff | ! | # #GA: Soil Humidity Geometric Average | ! | # #BP: Parameterization of Terrain Bumps | ! | | ! | | ! +------------------------------------------------------------------------+ ! +--Global Variables ! + ================ use VARphy use VAR_SV use VARdSV use VAR0SV use VARxSV use VARySV IMPLICIT NONE ! +--OUTPUT ! + ------ ! Water (Mass) Budget ! ~~~~~~~~~~~~~~~~~~~ ! #m0 real Wats_0(knonv) ! Soil Water, before forcing ! #m0 real Wats_1(knonv) ! Soil Water, after forcing ! #m0 real Wats_d(knonv) ! Soil Water forcing ! +--Internal Variables ! + ================== integer :: isl ,jsl ,ist ,ikl ! integer :: ikm ,ikp ,ik0 ,ik1 ! integer :: ist__s,ist__w ! Soil/Water Body Identifier ! #BP real z0soil ! Soil Surface Bumps Height [m] ! #BP real z_Bump !(Partly)Bumpy Layers Height [m] ! #BP real z0Bump ! Bumpy Layers Height [m] ! #BP real dzBump ! Lowest Bumpy Layer: ! #BP real etBump(knonv) ! Bumps Layer Averaged Humidity real :: etaMid ! Layer Interface's Humidity real :: Dhydif ! Hydraulic Diffusivity [m2/s] real :: eta__f ! Water Front Soil Water Content real :: Khyd_f ! Water Front Hydraulic Conduct. real :: Khydav ! Hydraulic Conductivity [m/s] real :: WgFlow ! Water gravitat. Flux [kg/m2/s] real :: Wg_MAX ! Water MAX.grav. Flux [kg/m2/s] real :: SatRat ! Saturation Flux [kg/m2/s] real :: WExces ! Saturat. Excess Flux [kg/m2/s] real :: SoRnOF(knonv) ! Soil Run OFF real :: Dhydtz(knonv,-nsol:0) ! Dhydif * dt / dz [m] real :: Elem_A,Elem_B,Elem_C ! Diagonal Coefficients real :: Diag_A(knonv,-nsol:0) ! A Diagonal real :: Diag_B(knonv,-nsol:0) ! B Diagonal real :: Diag_C(knonv,-nsol:0) ! C Diagonal real :: Term_D(knonv,-nsol:0) ! Independant Term real :: Aux__P(knonv,-nsol:0) ! P Auxiliary Variable real :: Aux__Q(knonv,-nsol:0) ! Q Auxiliary Variable real :: etaaux(knonv,-nsol:-nsol+1) ! Soil Water Content [m3/m3] real :: FreeDr ! Free Drainage Fraction (actual) real :: FreeD0 ! Free Drainage Fraction (1=Full) real :: aKdtSV3( 0:nsot, 0:nkhy) ! Khyd=a*eta+b: a * dt real :: bKdtSV3( 0:nsot, 0:nkhy) ! Khyd=a*eta+b: b * dt ! Water (Mass) Budget ! ~~~~~~~~~~~~~~~~~~~ ! #mw logical mwopen ! IO Switch ! #mw common/Sm_qSo_L/mwopen ! ! #mw real hourwr,timewr ! ! #mw common/Sm_qSo_R/timewr ! ! #mw real Evapor(knonv) ! ! +--Internal DATA ! + ============= ! #BP data z0soil/0.020/ ! Soil Surface Bumps Height [m] data FreeD0/1.000/ ! Free Drainage Fraction (1=Full) aKdtSV3=aKdtSV2*dt__SV bKdtSV3=bKdtSV2*dt__SV ! Water Budget (IN) ! ================== ! #m0 DO ikl=1,knonv ! #m0 Wats_0(ikl) = 0. ! OLD RunOFF Contrib. ! #m0 Wats_d(ikl) = drr_SV(ikl) ! Water Surface Forc. ! #m0 END DO ! #m0 isl= -nsol ! #m0 DO ikl=1,knonv ! #m0 Wats_0(ikl) = Wats_0(ikl) ! #m0. + ro_Wat *( eta_SV(ikl,isl) *dz78SV(isl) ! #m0. + eta_SV(ikl,isl+1) *dz_8SV(isl) ) * LSdzsv(ikl) ! #m0 END DO ! #m0 DO isl= -nsol+1,-1 ! #m0 DO ikl=1,knonv ! #m0 Wats_0(ikl) = Wats_0(ikl) ! #m0. + ro_Wat *( eta_SV(ikl,isl) *dz34SV(isl) ! #m0. +(eta_SV(ikl,isl-1) ! #m0. +eta_SV(ikl,isl+1))*dz_8SV(isl) ) * LSdzsv(ikl) ! #m0 END DO ! #m0 END DO ! #m0 isl= 0 ! #m0 DO ikl=1,knonv ! #m0 Wats_0(ikl) = Wats_0(ikl) ! #m0. + ro_Wat *( eta_SV(ikl,isl) *dz78SV(isl) ! #m0. + eta_SV(ikl,isl-1) *dz_8SV(isl) ) * LSdzsv(ikl) ! #m0 END DO ! +--Gravitational Flow ! + ================== ! +... METHOD: Surface Water Flux saturates successively the soil layers ! + ^^^^^^ from up to below, but is limited by infiltration capacity. ! + Hydraulic Conductivity again contributes after this step, ! + not redundantly because of a constant (saturated) profile. ! +--Flux Limitor ! + ^^^^^^^^^^^^^ isl=0 DO ikl=1,knonv ist = isotSV(ikl) ! Soil Type ist__s = min(ist, 1) ! 1 => Soil ist__w = 1 - ist__s ! 1 => Water Body Dhydif = s1__SV(ist) & *max(epsi,eta_SV(ikl,isl)) & ! Hydraulic Diffusivity **(bCHdSV(ist)+2.) ! DR97, Eqn.(3.36) Dhydif = ist__s * Dhydif & ! + ist__w * vK_dSV ! Water Bodies ! + Khydav = ist__s * Ks_dSV(ist) & ! DR97 Assumption + ist__w * vK_dSV ! Water Bodies ! + Wg_MAX = ro_Wat *Dhydif & ! MAXimum Infiltration *(etadSV(ist)-eta_SV(ikl,isl)) & ! Rate /(dzAvSV(isl)*LSdzsv(ikl) ) & ! + ro_Wat *Khydav ! ! +--Surface Horton RunOFF ! + ^^^^^^^^^^^^^^^^^^^^^ SoRnOF(ikl) = & max(zero,drr_SV(ikl)-Wg_MAX) RuofSV(ikl,1) = RuofSV(ikl,1) + SoRnOF(ikl) drr_SV(ikl) = drr_SV(ikl)-SoRnOF(ikl) RuofSV(ikl,2) = RuofSV(ikl,2) +max(0.,drr_SV(ikl)) END DO ! #GF DO isl=0,-nsol,-1 ! #GF DO ikl=1,knonv ! #GF ist = isotSV(ikl) ! Soil Type ! #GF ist__s = min(ist, 1) ! 1 => Soil ! #GF ist__w = 1 - ist__s ! 1 => Water Body ! +--Water Diffusion ! + ^^^^^^^^^^^^^^^ ! #GF Dhydif = s1__SV(ist) ! #GF. *max(epsi,eta_SV(ikl,isl)) ! Hydraulic Diffusivity ! #GF. **(bCHdSV(ist)+2.) ! DR97, Eqn.(3.36) ! #GF Dhydif = ist__s * Dhydif ! ! #GF. + ist__w * vK_dSV ! Water Bodies ! +--Water Conduction (without Horton Runoff) ! + ^^^^^^^^^^^^^^^^ ! #GF Khyd_f = Ks_dSV(ist) ! +... Uses saturated K ==> Horton Runoff ~0 ! ! +--Water Conduction (with Horton Runoff) ! + ^^^^^^^^^^^^^^^^ ! #GH ik0 = nkhy *eta_SV(ikl,isl) ! #GH. /etadSV(ist) ! #GH eta__f = 1. ! #GH. -aKdtSV3(ist,ik0)/(2. *dzAvSV(isl) ! #GH. *LSdzsv(ikl)) ! #GH eta__f = max(eps_21,eta__f) ! #GH eta__f = min(etadSV(ist), ! #GH. eta_SV(ikl,isl) + ! #GH. (aKdtSV3(ist,ik0) *eta_SV(ikl,isl) ! #GH. +bKdtSV3(ist,ik0)) /(dzAvSV(isl) ! #GH. *LSdzsv(ikl)) ! #GH. / eta__f ) ! #GH eta__f = .5*(eta_SV(ikl,isl) ! #GH. +eta__f) ! #gh eta__f = eta_SV(ikl,isl) ! #GH ik0 = nkhy *eta__f ! #GH. /etadSV(ist) ! #GH Khyd_f = ! #GH. (aKdtSV3(ist,ik0) *eta__f ! #GH. +bKdtSV3(ist,ik0)) /dt__SV ! #GF Khydav = ist__s * Khyd_f ! DR97 Assumption ! #GF. + ist__w * vK_dSV ! Water Bodies ! +--Gravitational Flow ! + ^^^^^^^^^^^^^^^^^^ ! #GF Wg_MAX = ! MAXimum Infiltration ! #GF. ro_Wat *Dhydif ! Rate ! #GF. *(etadSV(ist)-eta_SV(ikl,isl)) ! ! #GF. /(dzAvSV(isl)*LSdzsv(ikl) ) ! ! #GF. + ro_Wat *Khydav ! ! #GF END DO ! #GF END DO ! #GF DO ikl=1,knonv ! #GF SoRnOF(ikl) = SoRnOF(ikl) ! RunOFF Intensity ! #GF. + drr_SV(ikl) ! [kg/m2/s] ! +!!! Inclure la possibilite de creer une mare sur un bedrock impermeable ! #GF drr_SV(ikl) = 0. ! #GF END DO ! +--Temperature Correction due to a changed Soil Energy Content ! + =========================================================== ! +!!! Mettre en oeuvre le couplage humidit?-?nergie ! +--Full Resolution of the Richard's Equation ! + ========================================= ! +... METHOD: Water content evolution results from water fluxes ! + ^^^^^^ at the layer boundaries ! + Conductivity is approximated by a piecewise linear profile. ! + Semi-Implicit Crank-Nicholson scheme is used. ! + (Bruen, 1997, Sensitivity of hydrological processes ! + at the land-atmosphere interface. ! + Proc. Royal Irish Academy, IGBP symposium ! + on global change and the Irish Environment. ! + Publ.: Maynooth) ! + - - - - - - - - isl+1/2 - - ^ ! + | ! + eta_SV(isl) --------------- isl ----- +--dz_dSV(isl) ^ ! + | | ! + Dhydtz(isl) etaMid - - - - - - - - isl-1/2 - - v dzmiSV(isl)--+ ! + | ! + eta_SV(isl-1) --------------- isl-1 ----- v ! +--Transfert Coefficients ! + ---------------------------- DO isl=-nsol+1,0 DO ikl=1,knonv ist = isotSV(ikl) ! Soil Type ist__s = min(ist, 1) ! 1 => Soil ist__w = 1 - ist__s ! 1 => Water Body etaMid = (dz_dSV(isl) *eta_SV(ikl,isl-1) & ! eta at layers +dz_dSV(isl-1)*eta_SV(ikl,isl) ) & ! interface /(2.0* dzmiSV(isl)) ! LSdzsv implicit ! ! #GA etaMid = sqrt(dz_dSV(isl) *eta_SV(ikl,isl-1) ! Idem, geometric ! #GA. *dz_dSV(isl-1)*eta_SV(ikl,isl) ) ! average ! #GA. /(2.0* dzmiSV(isl)) ! (Vauclin&al.1979) Dhydif = s1__SV(ist) & ! Hydraul.Diffusi. *(etaMid **( bCHdSV(ist)+2.)) ! DR97, Eqn.(3.36) Dhydtz(ikl,isl) = Dhydif*dt__SV & ! /(dzmiSV(isl) & ! *LSdzsv(ikl)) ! Dhydtz(ikl,isl) = Dhydtz(ikl,isl) * ist__s & ! Soil +0.5*dzmiSV(isl)*LSdzsv(ikl) * ist__w ! Water bodies END DO END DO isl=-nsol DO ikl=1,knonv Dhydtz(ikl,isl) = 0.0 ! END DO ! +--Tridiagonal Elimination: Set Up ! + ------------------------------- ! +--Soil/Snow Interior ! + ^^^^^^^^^^^^^^^^^^ DO isl=0,-nsol,-1 DO ikl=1,knonv ist = isotSV(ikl) eta_SV(ikl,isl) = max(epsi, eta_SV(ikl,isl)) END DO END DO DO isl=-nsol,-nsol+1 DO ikl=1,knonv etaaux(ikl,isl) = eta_SV(ikl,isl) END DO END DO DO isl=-nsol+1,-1 DO ikl=1,knonv ist = isotSV(ikl) ikm = nkhy * eta_SV(ikl,isl-1) / etadSV(ist) ik0 = nkhy * eta_SV(ikl,isl) / etadSV(ist) ikp = nkhy * eta_SV(ikl,isl+1) / etadSV(ist) if(ikm<0.or.ik0<0.or.ikp<0)then print *,"CRASH1 in sisvat_qso.f on pixel (i,j,n)", & ii__SV(ikl),jj__SV(ikl),nn__SV(ikl) print *,"decrease your time step or increase ntphys "// & "and ntdiff in time_steps.f" stop endif Elem_A = Dhydtz(ikl,isl) & - aKdtSV3(ist,ikm)* dziiSV(isl) *LSdzsv(ikl) Elem_B = - (Dhydtz(ikl,isl) & +Dhydtz(ikl,isl+1) & -aKdtSV3(ist,ik0)*(dziiSV(isl+1) & -dzi_SV(isl) )*LSdzsv(ikl)) Elem_C = Dhydtz(ikl,isl+1) & + aKdtSV3(ist,ikp)* dzi_SV(isl+1)*LSdzsv(ikl) Diag_A(ikl,isl) = dz_8SV(isl) *LSdzsv(ikl) & -Implic * Elem_A Diag_B(ikl,isl) = dz34SV(isl) *LSdzsv(ikl) & -Implic * Elem_B Diag_C(ikl,isl) = dz_8SV(isl) *LSdzsv(ikl) & -Implic * Elem_C Term_D(ikl,isl) = (dz_8SV(isl) *LSdzsv(ikl) & +Explic *Elem_A )*eta_SV(ikl,isl-1) & + (dz34SV(isl) *LSdzsv(ikl) & +Explic *Elem_B )*eta_SV(ikl,isl) & + (dz_8SV(isl) *LSdzsv(ikl) & +Explic *Elem_C )*eta_SV(ikl,isl+1) & + (bKdtSV3(ist,ikp)* dzi_SV(isl+1) & +bKdtSV3(ist,ik0)*(dziiSV(isl+1) & -dzi_SV(isl) ) & -bKdtSV3(ist,ikm)* dziiSV(isl) ) & * LSdzsv(ikl) END DO END DO isl=-nsol DO ikl=1,knonv ist = isotSV(ikl) ! # FreeDr = FreeD0 * min(ist,1) FreeDr = iWaFSV(ikl) * min(ist,1) ik0 = nkhy * eta_SV(ikl,isl ) / etadSV(ist) ikp = nkhy * eta_SV(ikl,isl+1) / etadSV(ist) if(ik0<0.or.ikp<0)then print *,"CRASH2 in sisvat_qso.f on pixel (i,j,n)", & ii__SV(ikl),jj__SV(ikl),nn__SV(ikl) print *,"decrease your time step or increase ntphys "// & "and ntdiff in time_steps.f" stop endif Elem_A = 0. Elem_B = - (Dhydtz(ikl,isl+1) & -aKdtSV3(ist,ik0)*(dziiSV(isl+1)*LSdzsv(ikl) & -FreeDr )) Elem_C = Dhydtz(ikl,isl+1) & + aKdtSV3(ist,ikp)* dzi_SV(isl+1)*LSdzsv(ikl) Diag_A(ikl,isl) = 0. Diag_B(ikl,isl) = dz78SV(isl) *LSdzsv(ikl) & -Implic *Elem_B Diag_C(ikl,isl) = dz_8SV(isl) *LSdzsv(ikl) & -Implic *Elem_C Term_D(ikl,isl) = (dz78SV(isl) *LSdzsv(ikl) & +Explic *Elem_B )*eta_SV(ikl,isl) & + (dz_8SV(isl) *LSdzsv(ikl) & +Explic *Elem_C )*eta_SV(ikl,isl+1) & + (bKdtSV3(ist,ikp)* dzi_SV(isl+1)*LSdzsv(ikl) & +bKdtSV3(ist,ik0)*(dziiSV(isl+1)*LSdzsv(ikl) & -FreeDr )) END DO isl=0 DO ikl=1,knonv ist = isotSV(ikl) ikm = nkhy * eta_SV(ikl,isl-1) / etadSV(ist) ik0 = nkhy * eta_SV(ikl,isl) / etadSV(ist) Elem_A = Dhydtz(ikl,isl) & - aKdtSV3(ist,ikm)* dziiSV(isl)*LSdzsv(ikl) Elem_B = - (Dhydtz(ikl,isl) & +aKdtSV3(ist,ik0)* dzi_SV(isl)*LSdzsv(ikl)) Elem_C = 0. Diag_A(ikl,isl) = dz_8SV(isl) *LSdzsv(ikl) & - Implic *Elem_A Diag_B(ikl,isl) = dz78SV(isl) *LSdzsv(ikl) & - Implic *Elem_B Diag_C(ikl,isl) = 0. ! + Term_D(ikl,isl) = (dz_8SV(isl) *LSdzsv(ikl) & +Explic *Elem_A )*eta_SV(ikl,isl-1) & + (dz78SV(isl) *LSdzsv(ikl) & +Explic *Elem_B )*eta_SV(ikl,isl) & - (bKdtSV3(ist,ik0)* dzi_SV(isl) & +bKdtSV3(ist,ikm)* dziiSV(isl))*LSdzsv(ikl) & + dt__SV *(HLs_sv(ikl) * (1-min(1,isnoSV(ikl))) & / (ro_Wat *dz_dSV(0) * Lx_H2O(ikl)) & !XF bug 17/05/2017 +drr_SV(ikl))/ro_Wat END DO DO ikl=1,knonv drr_SV(ikl)=0. ! drr is included in the 1st soil layer ENDDO ! + ! + ! +--Tridiagonal Elimination ! + ======================= ! + ! +--Forward Sweep ! + ^^^^^^^^^^^^^^ DO ikl= 1,knonv Aux__P(ikl,-nsol) = Diag_B(ikl,-nsol) Aux__Q(ikl,-nsol) =-Diag_C(ikl,-nsol)/Aux__P(ikl,-nsol) END DO ! + DO isl=-nsol+1,0 DO ikl= 1,knonv Aux__P(ikl,isl) = Diag_A(ikl,isl) *Aux__Q(ikl,isl-1) & +Diag_B(ikl,isl) Aux__Q(ikl,isl) =-Diag_C(ikl,isl) /Aux__P(ikl,isl) END DO END DO ! + DO ikl= 1,knonv eta_SV(ikl,-nsol) = Term_D(ikl,-nsol)/Aux__P(ikl,-nsol) END DO ! + DO isl=-nsol+1,0 DO ikl= 1,knonv eta_SV(ikl,isl) =(Term_D(ikl,isl) & -Diag_A(ikl,isl) *eta_SV(ikl,isl-1)) & /Aux__P(ikl,isl) END DO END DO ! +--Backward Sweep ! + ^^^^^^^^^^^^^^ DO isl=-1,-nsol,-1 DO ikl= 1,knonv eta_SV(ikl,isl) = Aux__Q(ikl,isl) *eta_SV(ikl,isl+1) & +eta_SV(ikl,isl) END DO END DO ! +--Horton RunOFF Intensity ! + ======================= DO isl=0,-nsol,-1 DO ikl=1,knonv ist = isotSV(ikl) ! Soil Type SatRat = (eta_SV(ikl,isl)-etadSV(ist)) & ! OverSaturation Rate *ro_Wat *dzAvSV(isl) & ! *LSdzsv(ikl) & ! /dt__SV ! SoRnOF(ikl) = SoRnOF(ikl) & ! + max(zero,SatRat) ! RuofSV(ikl,3) = RuofSV(ikl,3) + & max(zero,SatRat) eta_SV(ikl,isl) = max(epsi & ! ! #ED. +etamSV(isotSV(ikl))! ,eta_SV(ikl,isl)) ! eta_SV(ikl,isl) = min(eta_SV(ikl,isl) & ! ,etadSV(ist) ) ! END DO END DO ! +--IO, for Verification ! + ~~~~~~~~~~~~~~~~~~~~ ! #WR write(6,6010) 6010 format(/,1x) DO isl= 0,-nsol,-1 DO ikl= 1,knonv ist = isotSV(ikl) ikp = nkhy * eta_SV(ikl,isl) /etadSV(ist) Khydsv(ikl,isl) =(aKdtSV3(ist,ikp) *eta_SV(ikl,isl) & +bKdtSV3(ist,ikp)) *2.0/dt__SV ! #WR write(6,6011) ikl,isl,eta_SV(ikl,isl)*1.e3, ! #WR. ikp, aKdtSV3(ist,ikp),bKdtSV3(ist,ikp), ! #WR. Khydsv(ikl,isl) 6011 format(2i3,f8.1,i3,3e12.3) END DO END DO ! +--Additional RunOFF Intensity ! + =========================== DO ikl=1,knonv ist = isotSV(ikl) ik0 = nkhy * etaaux(ikl,-nsol ) /etadSV(ist) ! # FreeDr = FreeD0 * min(ist,1) FreeDr = iWaFSV(ikl) * min(ist,1) SoRnOF(ikl) = SoRnOF(ikl) & + (aKdtSV3(ist,ik0)*(etaaux(ikl,-nsol)*Explic & +eta_SV(ikl,-nsol)*Implic) & + bKdtSV3(ist,ik0) ) & * FreeDr *ro_Wat /dt__SV RuofSV(ikl,3) = RuofSV(ikl,3) & + (aKdtSV3(ist,ik0)*(etaaux(ikl,-nsol)*Explic & +eta_SV(ikl,-nsol)*Implic) & + bKdtSV3(ist,ik0) ) & * FreeDr *ro_Wat /dt__SV ! +--Full Run OFF: Update ! + ~~~~~~~~~~~~~~~~~~~~ RnofSV(ikl) = RnofSV(ikl) + SoRnOF(ikl) RuofSV(ikl,4) = RuofSV(ikl,4) + SoRnOF(ikl) END DO ! +--Temperature Correction due to a changed Soil Energy Content ! + =========================================================== ! +!!! Mettre en oeuvre le couplage humidit?-?nergie ! +--Bumps/Asperites Treatment ! + ========================= ! +--Average over Bump Depth (z0soil) ! + -------------------------------- ! #BP z_Bump = 0. ! #BP DO ikl=1,knonv ! #BP etBump(ikl) = 0. ! #BP END DO ! + ! #BP DO isl=0,-nsol,-1 ! #BP z0Bump = z_Bump ! #BP z_Bump = z_Bump + dzAvSV(isl) ! #BP IF (z_Bump.lt.z0soil) THEN ! #BP DO ikl=1,knonv ! #BP etBump(ikl) = etBump(ikl) + dzAvSV(isl) *eta_SV(ikl,isl) ! #BP END DO ! #BP END IF ! #BP IF (z_Bump.gt.z0soil.AND.z0Bump.lt.z0soil) THEN ! #BP DO ikl=1,knonv ! #BP etBump(ikl) = etBump(ikl) + (z0soil-z0Bump)*eta_SV(ikl,isl) ! #BP etBump(ikl) = etBump(ikl) / z0soil ! #BP END DO ! #BP END IF ! #BP END DO ! +--Correction ! + ---------- ! #BP z_Bump = 0. ! #BP DO isl=0,-nsol,-1 ! #BP z0Bump = z_Bump ! #BP z_Bump = z_Bump +dzAvSV(isl) ! #BP IF (z_Bump.lt.z0soil) THEN ! #BP DO ikl=1,knonv ! #BP eta_SV(ikl,isl) = etBump(ikl) ! #BP END DO ! #BP END IF ! #BP IF (z_Bump.gt.z0soil.AND.z0Bump.lt.z0soil) THEN ! #BP dzBump = z_Bump - z0soil ! #BP DO ikl=1,knonv ! #BP eta_SV(ikl,isl) =(etBump(ikl) *(dzAvSV(isl)-dzBump) ! #BP. + eta_SV(ikl,isl)* dzBump) ! #BP. / dzAvSV(isl) ! #BP END DO ! #BP END IF ! #BP END DO ! +--Positive Definite ! + ================= ! #BP DO isl= 0,-nsol,-1 ! #BP DO ikl= 1,knonv ! #BP eta_SV(ikl,isl) = max(epsi,eta_SV(ikl,isl)) ! #BP END DO ! #BP END DO ! +--Water Budget (OUT) ! + =================== ! #m0 DO ikl=1,knonv ! #m0 Wats_d(ikl) = Wats_d(ikl) ! ! #m0. + drr_SV(ikl) *zero ! Precipitation is ! + \______________ already included ! #m0. + HLs_sv(ikl) ! #m0. *(1-min(isnoSV(ikl),1)) /Lx_H2O(ikl) ! Evaporation ! #m0. - SoRnOF(ikl) ! Soil RunOFF Contrib. ! #m0 Wats_1(ikl) = 0. ! ! #mw Evapor(ikl) = HLs_sv(ikl) *dt__SV ! ! #mw. *(1-min(isnoSV(ikl),1)) /Lx_H2O(ikl) ! ! #m0 END DO ! #m0 DO isl= -nsol,0 ! #m0 DO ikl=1,knonv ! #m0 Wats_d(ikl) = Wats_d(ikl) ! ! #m0 END DO ! #m0 END DO ! #m0 DO ikl=1,knonv ! #m0 Wats_d(ikl) = Wats_d(ikl) *dt__SV ! ! #m0 END DO ! #m0 isl= -nsol ! #m0 DO ikl=1,knonv ! #m0 Wats_1(ikl) = Wats_1(ikl) ! #m0. + ro_Wat *( eta_SV(ikl,isl) *dz78SV(isl) ! #m0. + eta_SV(ikl,isl+1) *dz_8SV(isl) ) *LSdzsv(ikl) ! #m0 END DO ! #m0 DO isl= -nsol+1,-1 ! #m0 DO ikl=1,knonv ! #m0 Wats_1(ikl) = Wats_1(ikl) ! #m0. + ro_Wat *( eta_SV(ikl,isl) *dz34SV(isl) ! #m0. +(eta_SV(ikl,isl-1) ! #m0. +eta_SV(ikl,isl+1))*dz_8SV(isl) ) *LSdzsv(ikl) ! #m0 END DO ! #m0 END DO ! #m0 isl= 0 ! #m0 DO ikl=1,knonv ! #m0 Wats_1(ikl) = Wats_1(ikl) ! #m0. + ro_Wat *( eta_SV(ikl,isl) *dz78SV(isl) ! #m0. + eta_SV(ikl,isl-1) *dz_8SV(isl) ) *LSdzsv(ikl) ! #m0 END DO return end subroutine sisvat_qso