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