subroutine SISVAT_TSo ! #e1. (ETSo_0,ETSo_1,ETSo_d) ! +------------------------------------------------------------------------+ ! | MAR SISVAT_TSo 06-10-2020 MAR | ! | SubRoutine SISVAT_TSo computes the Soil/Snow Energy Balance | ! +------------------------------------------------------------------------+ ! | | ! | PARAMETERS: knonv: Total Number of columns = | ! | ^^^^^^^^^^ = Total Number of continental grid boxes | ! | X Number of Mosaic Cell per grid box | ! | | ! | INPUT: isotSV = 0,...,11: Soil Type | ! | ^^^^^ 0: Water, Solid or Liquid | ! | isnoSV = total Nb of Ice/Snow Layers | ! | dQa_SV = Limitation of Water Vapor Turbulent Flux | ! | | ! | INPUT: sol_SV : Downward Solar Radiation [W/m2] | ! | ^^^^^ IRd_SV : Surface Downward Longwave Radiation [W/m2] | ! | za__SV : SBL Top Height [m] | ! | VV__SV : SBL Top Wind Speed [m/s] | ! | TaT_SV : SBL Top Temperature [K] | ! | rhT_SV : SBL Top Air Density [kg/m3] | ! | QaT_SV : SBL Top Specific Humidity [kg/kg] | ! | LSdzsv : Vertical Discretization Factor [-] | ! | = 1. Soil | ! | = 1000. Ocean | ! | dzsnSV : Snow Layer Thickness [m] | ! | ro__SV : Snow/Soil Volumic Mass [kg/m3] | ! | eta_SV : Soil Water Content [m3/m3] | ! | dt__SV : Time Step [s] | ! | | ! | SoSosv : Absorbed Solar Radiation by Surfac.(Normaliz)[-] | ! | Eso_sv : Soil+Snow Emissivity [-] | ! | rah_sv : Aerodynamic Resistance for Heat [s/m] | ! | Lx_H2O : Latent Heat of Vaporization/Sublimation [J/kg] | ! | sEX_sv : Verticaly Integrated Extinction Coefficient [-] | ! | | ! | INPUT / TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)| ! | OUTPUT: & Snow Temperatures (layers 1,2,...,nsno) [K] | ! | ^^^^^^ | ! | | ! | OUTPUT: IRs_SV : Soil IR Radiation [W/m2] | ! | ^^^^^^ HSs_sv : Sensible Heat Flux [W/m2] | ! | HLs_sv : Latent Heat Flux [W/m2] | ! | ETSo_0 : Snow/Soil Energy Power, before Forcing [W/m2] | ! | ETSo_1 : Snow/Soil Energy Power, after Forcing [W/m2] | ! | ETSo_d : Snow/Soil Energy Power Forcing [W/m2] | ! | | ! | Internal Variables: | ! | ^^^^^^^^^^^^^^^^^^ | ! | | ! | METHOD: NO Skin Surface Temperature | ! | ^^^^^^ Semi-Implicit Crank Nicholson Scheme | ! | | ! | # OPTIONS: #E0: Energy Budget Verification | ! | # ^^^^^^^ #kd: KDsvat Option:NO Flux Limitor on HL | ! | # #KD: KDsvat Option:Explicit Formulation of HL | ! | # #NC: OUTPUT for Stand Alone NetCDF File | ! | | ! +------------------------------------------------------------------------+ ! +--Global Variables ! + ================ use VARphy use VAR_SV use VARdSV use VARxSV use VARySV use VARtSV use VAR0SV IMPLICIT NONE ! +--OUTPUT ! + ------ ! #e1 real ETSo_0(knonv) ! Soil/Snow Power, before Forcing ! #e1 real ETSo_1(knonv) ! Soil/Snow Power, after Forcing ! #e1 real ETSo_d(knonv) ! Soil/Snow Power, Forcing ! +--Internal Variables ! + ================== integer :: ikl ,isl ,jsl ,ist ! integer :: ist__s,ist__w ! Soil/Water Body Identifier integer :: islsgn ! Soil/Snow Surfac.Identifier real :: eps__3 ! Arbitrary Low Number real :: etaMid,psiMid ! Layer Interface's Humidity real :: mu_eta ! Soil thermal Conductivity real :: mu_exp ! arg Soil thermal Conductivity real :: mu_min ! Min Soil thermal Conductivity real :: mu_max ! Max Soil thermal Conductivity real :: mu_sno(knonv),mu_aux ! Snow thermal Conductivity real :: mu__dz(knonv,-nsol:nsno+1) ! mu_(eta,sno) / dz real :: dtC_sv(knonv,-nsol:nsno) ! dt / C real :: IRs__D(knonv) ! UpwardIR Previous Iter.Contr. real :: dIRsdT(knonv) ! UpwardIR T Derivat. real :: f_HSHL(knonv) ! Factor common to HS and HL real :: dRidTs(knonv) ! d(Rib)/d(Ts) real :: HS___D(knonv) ! Sensible Heat Flux Atm.Contr. real :: f___HL(knonv) ! real :: HL___D(knonv) ! Latent Heat Flux Atm.Contr. REAL :: TSurf0(knonv),dTSurf ! Previous Surface Temperature real :: qsatsg(knonv) !,den_qs,arg_qs ! Soil Saturat. Spec. Humidity real :: dqs_dT(knonv) ! d(qsatsg)/dTv real :: Psi( knonv) ! 1st Soil Layer Water Potential real :: RHuSol(knonv) ! Soil Surface Relative Humidity real :: etaSol ! Soil Surface Humidity real :: d__eta ! Soil Surface Humidity Increm. real :: Elem_A,Elem_C ! Diagonal Coefficients real :: Diag_A(knonv,-nsol:nsno) ! A Diagonal real :: Diag_B(knonv,-nsol:nsno) ! B Diagonal real :: Diag_C(knonv,-nsol:nsno) ! C Diagonal real :: Term_D(knonv,-nsol:nsno) ! Independant Term real :: Aux__P(knonv,-nsol:nsno) ! P Auxiliary Variable real :: Aux__Q(knonv,-nsol:nsno) ! Q Auxiliary Variable real :: Ts_Min,Ts_Max ! Temperature Limits ! #e1 real Exist0 ! Existing Layer Switch real :: psat_wat, psat_ice, sp ! computation of qsat integer :: nt_srf,it_srf,itEuBk ! HL: Surface Scheme parameter(nt_srf=10) ! 10 before real :: agpsrf,xgpsrf,dt_srf,dt_ver ! real :: etaBAK(knonv) ! real :: etaNEW(knonv) ! real :: etEuBk(knonv) ! real :: fac_dt(knonv),faceta(knonv) ! real :: PsiArg(knonv),SHuSol(knonv) ! ! +--Internal DATA ! + ============= data eps__3 / 1.e-3 / ! Arbitrary Low Number data mu_exp / -0.4343 / ! Soil Thermal Conductivity data mu_min / 0.172 / ! Min Soil Thermal Conductivity data mu_max / 2.000 / ! Max Soil Thermal Conductivity data Ts_Min / 175. / ! Temperature Minimum data Ts_Max / 300. / ! Temperature Acceptable Maximum ! + ! including Snow Melt Energy ! +-- Initilialisation of local arrays ! + ================================ DO ikl=1,knonv mu_sno(ikl)=0. mu__dz(ikl,:)=0. dtC_sv(ikl,:)=0. IRs__D(ikl)=0. dIRsdT(ikl)=0. f_HSHL(ikl)=0. dRidTs(ikl)=0. HS___D(ikl)=0. f___HL(ikl)=0. HL___D(ikl)=0. TSurf0(ikl)=0. qsatsg(ikl)=0. dqs_dT(ikl)=0. Psi(ikl)=0. RHuSol(ikl)=0. Diag_A(ikl,:)=0. Diag_B(ikl,:)=0. Diag_C(ikl,:)=0. Term_D(ikl,:)=0. Aux__P(ikl,:)=0. Aux__Q(ikl,:)=0. etaBAK(ikl)=0. etaNEW(ikl)=0. etEuBk(ikl)=0. fac_dt(ikl)=0. faceta(ikl)=0. PsiArg(ikl)=0. SHuSol(ikl)=0. END DO ! +--Heat Conduction Coefficient (zero in the Layers over the highest one) ! + =========================== ! + ---------------- isl eta_SV, rho C (isl) ! + ! +--Soil ++++++++++++++++ etaMid, mu (isl) ! + ---- ! + ---------------- isl-1 eta_SV, rho C (isl-1) isl=-nsol DO ikl=1,knonv mu__dz(ikl,isl) = 0. dtC_sv(ikl,isl) = dtz_SV2(isl) * dt__SV & ! dt / (dz X rho C) /((rocsSV(isotSV(ikl)) & ! [s / (m.J/m3/K)] +rcwdSV*eta_SV(ikl,isl)) & ! *LSdzsv(ikl) ) ! END DO 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 = 0.5*(dz_dSV(isl-1)*eta_SV(ikl,isl-1) & ! eta at layers +dz_dSV(isl) *eta_SV(ikl,isl) ) & ! interface /dzmiSV(isl) ! LSdzsv implicit ! etaMid = max(etaMid,epsi) psiMid = psidSV(ist) & *(etadSV(ist)/etaMid)**bCHdSV(ist) mu_eta = 3.82 *(psiMid)**mu_exp ! Soil Thermal mu_eta = min(max(mu_eta, mu_min), mu_max) ! Conductivity ! + ! DR97 eq.3.31 mu_eta = ist__s *mu_eta +ist__w * vK_dSV ! Water Bodies ! + ! Correction mu__dz(ikl,isl) = mu_eta/(dzmiSV(isl) & ! *LSdzsv(ikl)) ! dtC_sv(ikl,isl) = dtz_SV2(isl)* dt__SV & ! dt / (dz X rho C) /((rocsSV(isotSV(ikl)) & ! +rcwdSV*eta_SV(ikl,isl)) & ! *LSdzsv(ikl) ) ! END DO END DO ! +--Soil/Snow Interface ! + ------------------- ! +--Soil Contribution ! + ^^^^^^^^^^^^^^^^^ isl=1 DO ikl=1,knonv ist = isotSV(ikl) ! Soil Type ist__s = min(ist, 1) ! 1 => Soil ist__w = 1 - ist__s ! 1 => Water Body psiMid = psidSV(ist) ! Snow => Saturation mu_eta = 3.82 *(psiMid)**mu_exp ! Soil Thermal mu_eta = min(max(mu_eta, mu_min), mu_max) ! Conductivity ! + ! DR97 eq.3.31 mu_eta = ist__s *mu_eta +ist__w * vK_dSV ! Water Bodies ! +--Snow Contribution ! + ^^^^^^^^^^^^^^^^^ mu_sno(ikl) = CdidSV & ! *(ro__SV(ikl,isl) /ro_Wat) ** 1.88 ! mu_sno(ikl) = max(epsi,mu_sno(ikl)) ! ! +... mu_sno : Snow Heat Conductivity Coefficient [Wm/K] ! + (Yen 1981, CRREL Rep., 81-10) ! +--Combined Heat Conductivity ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^ mu__dz(ikl,isl) = 2./(dzsnSV(ikl,isl ) & ! Combined Heat /mu_sno(ikl) & ! Conductivity +LSdzsv(ikl) & ! *dz_dSV( isl-1)/mu_eta) ! Coefficient ! +--Inverted Heat Capacity ! + ^^^^^^^^^^^^^^^^^^^^^^ dtC_sv(ikl,isl) = dt__SV/max(epsi, & ! dt / (dz X rho C) dzsnSV(ikl,isl) * ro__SV(ikl,isl) *Cn_dSV) ! END DO ! +--Snow ! + ---- DO ikl=1,knonv DO isl=1,min(nsno,isnoSV(ikl)+1) ro__SV(ikl,isl) = & ! ro__SV(ikl ,isl) & ! * max(0,min(isnoSV(ikl)-isl+1,1)) ! END DO END DO DO ikl=1,knonv DO isl=1,min(nsno,isnoSV(ikl)+1) ! +--Combined Heat Conductivity ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^ mu_aux = CdidSV & ! *(ro__SV(ikl,isl) /ro_Wat) ** 1.88 ! mu__dz(ikl,isl) = & ! 2. *mu_aux*mu_sno(ikl) & ! Combined Heat /max(epsi,dzsnSV(ikl,isl )*mu_sno(ikl) & ! Conductivity +dzsnSV(ikl,isl-1)*mu_aux ) ! For upper Layer mu_sno(ikl) = mu_aux ! ! +--Inverted Heat Capacity ! + ^^^^^^^^^^^^^^^^^^^^^^ dtC_sv(ikl,isl) = dt__SV/max(eps__3, & ! dt / (dz X rho C) dzsnSV(ikl,isl) * ro__SV(ikl,isl) *Cn_dSV) ! END DO END DO ! +--Uppermost Effective Layer: NO conduction ! + ---------------------------------------- DO ikl=1,knonv mu__dz(ikl,isnoSV(ikl)+1) = 0.0 END DO ! +--Energy Budget (IN) ! + ================== ! #e1 DO ikl=1,knonv ! #e1 ETSo_0(ikl) = 0. ! #e1 END DO ! #e1 DO isl= -nsol,nsno ! #e1 DO ikl=1,knonv ! #e1 Exist0 = isl - isnoSV(ikl) ! #e1 Exist0 = 1. - max(zero,min(unun,Exist0)) ! #e1 ETSo_0(ikl) = ETSo_0(ikl) ! #e1. +(TsisSV(ikl,isl)-TfSnow)*Exist0 ! #e1. /dtC_sv(ikl,isl) ! #e1 END DO ! #e1 END DO ! +--Tridiagonal Elimination: Set Up ! + =============================== ! +--Soil/Snow Interior ! + ^^^^^^^^^^^^^^^^^^ DO ikl=1,knonv DO isl=-nsol+1,min(nsno-1,isnoSV(ikl)+1) Elem_A = dtC_sv(ikl,isl) *mu__dz(ikl,isl) Elem_C = dtC_sv(ikl,isl) *mu__dz(ikl,isl+1) Diag_A(ikl,isl) = -Elem_A *Implic Diag_C(ikl,isl) = -Elem_C *Implic Diag_B(ikl,isl) = 1.0d+0 -Diag_A(ikl,isl)-Diag_C(ikl,isl) Term_D(ikl,isl) = Explic *(Elem_A *TsisSV(ikl,isl-1) & +Elem_C *TsisSV(ikl,isl+1)) & +(1.0d+0 -Explic *(Elem_A+Elem_C))*TsisSV(ikl,isl) & + dtC_sv(ikl,isl) * sol_SV(ikl) *SoSosv(ikl) & *(sEX_sv(ikl,isl+1) & -sEX_sv(ikl,isl )) END DO END DO ! +--Soil lowest Layer ! + ^^^^^^^^^^^^^^^^^^ isl= -nsol DO ikl=1,knonv Elem_A = 0. Elem_C = dtC_sv(ikl,isl) *mu__dz(ikl,isl+1) Diag_A(ikl,isl) = 0. Diag_C(ikl,isl) = -Elem_C *Implic Diag_B(ikl,isl) = 1.0d+0 -Diag_A(ikl,isl)-Diag_C(ikl,isl) Term_D(ikl,isl) = Explic * Elem_C *TsisSV(ikl,isl+1) & +(1.0d+0 -Explic * Elem_C) *TsisSV(ikl,isl) & + dtC_sv(ikl,isl) * sol_SV(ikl) *SoSosv(ikl) & *(sEX_sv(ikl,isl+1) & -sEX_sv(ikl,isl )) END DO ! +--Snow highest Layer (dummy!) ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! !EV!isl= min(isnoSV(1)+1,nsno) DO ikl=1,knonv ! EV try to calculate isl at the ikl grid point isl= min(isnoSV(ikl)+1,nsno) Elem_A = dtC_sv(ikl,isl) *mu__dz(ikl,isl) Elem_C = 0. Diag_A(ikl,isl) = -Elem_A *Implic Diag_C(ikl,isl) = 0. Diag_B(ikl,isl) = 1.0d+0 -Diag_A(ikl,isl) Term_D(ikl,isl) = Explic * Elem_A *TsisSV(ikl,isl-1) & +(1.0d+0 -Explic * Elem_A) *TsisSV(ikl,isl) & + dtC_sv(ikl,isl) * (sol_SV(ikl) *SoSosv(ikl) & *(sEX_sv(ikl,isl+1) & -sEX_sv(ikl,isl ))) END DO ! +--Surface: UPwardIR Heat Flux ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^ DO ikl=1,knonv isl = isnoSV(ikl) dIRsdT(ikl) = Eso_sv(ikl)* StefBo * 4. & ! - d(IR)/d(T) * TsisSV(ikl,isl) & ! * TsisSV(ikl,isl) & ! * TsisSV(ikl,isl) ! IRs__D(ikl) = dIRsdT(ikl)* TsisSV(ikl,isl) * 0.75 ! ! +--Surface: Richardson Number: T Derivative ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! #RC dRidTs(ikl) =-gravit * za__SV(ikl) ! #RC. /(TaT_SV(ikl) * VV__SV(ikl) ! #RC. * VV__SV(ikl)) ! +--Surface: Turbulent Heat Flux: Factors ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ f_HSHL(ikl) = rhT_SV(ikl) / rah_sv(ikl) ! to HS, HL f___HL(ikl) = f_HSHL(ikl) * Lx_H2O(ikl) ! +--Surface: Sensible Heat Flux: T Derivative ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ dSdTSV(ikl) = f_HSHL(ikl) * Cp !#- d(HS)/d(T) ! #RC. *(1.0 -(TsisSV(ikl,isl) -TaT_SV(ikl)) !#Richardson ! #RC. * dRidTs(ikl)*dFh_sv(ikl)/rah_sv(ikl)) ! Nb. Correct. HS___D(ikl) = dSdTSV(ikl) * TaT_SV(ikl) ! ! +--Surface: Latent Heat Flux: Saturation Specific Humidity ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! den_qs = TsisSV(ikl,isl)- 35.8 ! ! arg_qs = 17.27 *(TsisSV(ikl,isl)-273.16) ! ! . / den_qs ! ! qsatsg(ikl) = .0038 * exp(arg_qs) ! ! sp = (pst_SV(ikl) + ptopSV) * 10. ! !sp=ps__SV(ikl) ! ! Etienne: in the formula herebelow sp should be in hPa, not ! ! in Pa so I divide by 100. sp=ps__SV(ikl)/100. psat_ice = 6.1070 * exp(6150. *(1./273.16 - & 1./TsisSV(ikl,isl))) psat_wat = 6.1078 * exp (5.138*log(273.16 /TsisSV(ikl,isl))) & * exp (6827.*(1. /273.16-1./TsisSV(ikl,isl))) if(TsisSV(ikl,isl)<=273.16) then qsatsg(ikl) = 0.622 * psat_ice / (sp - 0.378 * psat_ice) else qsatsg(ikl) = 0.622 * psat_wat / (sp - 0.378 * psat_wat) endif QsT_SV(ikl)=qsatsg(ikl) ! dqs_dT(ikl) = qsatsg(ikl)* 4099.2 /(den_qs *den_qs)! fac_dt(ikl) = f_HSHL(ikl)/(ro_Wat * dz_dSV(0)) ! END DO ! +--Surface: Latent Heat Flux: Surface Relative Humidity ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ xgpsrf = 1.05 ! agpsrf = dt__SV*( 1.0-xgpsrf ) & ! /( 1.0-xgpsrf**nt_srf) ! dt_srf = agpsrf ! dt_ver = 0. DO ikl=1,knonv isl = isnoSV(ikl) ist = max(0,isotSV(ikl)-100*isnoSV(ikl))! 0 if H2O ist__s = min(1,ist) etaBAK(ikl) = max(epsi,eta_SV(ikl ,isl)) ! etaNEW(ikl) = etaBAK(ikl) ! etEuBk(ikl) = etaNEW(ikl) ! END DO if(ist__s==1) then ! to reduce computer time ! ! DO it_srf=1,nt_srf ! dt_ver = dt_ver +dt_srf ! DO ikl=1,knonv ! faceta(ikl) = fac_dt(ikl)*dt_srf ! ! #VX faceta(ikl) = faceta(ikl) ! ! #VX. /(1.+faceta(ikl)*dQa_SV(ikl)) ! Limitation ! ! by Atm.Conten ! #??. *max(0,sign(1.,qsatsg(ikl)-QaT_SV(ikl)))) ! NO Limitation ! ! of Downw.Flux END DO ! DO itEuBk=1,2 ! DO ikl=1,knonv ist = max(0,isotSV(ikl)-100*isnoSV(ikl)) ! 0 if H2O ! ! Psi(ikl) = & ! psidSV(ist) & ! DR97, Eqn 3.34 *(etadSV(ist) & ! /max(etEuBk(ikl),epsi)) & ! **bCHdSV(ist) ! PsiArg(ikl) = 7.2E-5*Psi(ikl) ! RHuSol(ikl) = exp(-min(0.,PsiArg(ikl))) ! SHuSol(ikl) = qsatsg(ikl) *RHuSol(ikl) ! DR97, Eqn 3.15 etEuBk(ikl) = & ! (etaNEW(ikl) + faceta(ikl)*(QaT_SV(ikl) & ! -SHuSol(ikl) & ! *(1. -bCHdSV(ist) & ! *PsiArg(ikl)) )) & ! /(1. + faceta(ikl)* SHuSol(ikl) & ! *bCHdSV(ist) & ! *PsiArg(ikl) & ! /etaNEW(ikl)) ! etEuBk(ikl) = etEuBk(ikl) & ! ! . /(Ro_Wat*dz_dSV(0)) ! * dt_srf /(Ro_Wat*dz_dSV(0)) ! !XF 15/05/2017 BUG END DO ! END DO ! DO ikl=1,knonv ! etaNEW(ikl) = max(etEuBk(ikl),epsi) ! END DO ! dt_srf = dt_srf * xgpsrf ! END DO endif ! ! +--Surface: Latent Heat Flux: Soil/Water Surface Contributions ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ DO ikl=1,knonv ! isl = isnoSV(ikl) ! ist = max(0,isotSV(ikl)-100*isnoSV(ikl)) ! 0 if H2O ist__s= min(1,ist) ! 1 if no H2O ist__w= 1-ist__s ! 1 if H2O d__eta = eta_SV(ikl,isl)-etaNEW(ikl) ! ! ! latent heat flux computation HL___D(ikl)=( ist__s *ro_Wat *dz_dSV(0) & ! Soil Contrib. *(etaNEW(ikl) -etaBAK(ikl)) / dt__SV & ! +ist__w *f_HSHL(ikl) & ! H2O Contrib. *(QaT_SV(ikl) - qsatsg(ikl)) ) & ! * Lx_H2O(ikl) ! common factor ! #DL RHuSol(ikl) =(QaT_SV(ikl) ! ! #DL. -HL___D(ikl) / f___HL(ikl)) ! ! #DL. / qsatsg(ikl) ! ! +--Surface: Latent Heat Flux: T Derivative ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ dLdTSV(ikl) = 0. ! #DL dLdTSV(ikl) = f___HL(ikl) * RHuSol(ikl) *dqs_dT(ikl) ! - d(HL)/d(T) ! #DL HL___D(ikl) = HL___D(ikl) ! ! #DL. +dLdTSV(ikl) * TsisSV(ikl,isl) ! END DO ! ! +--Surface: Tridiagonal Matrix Set Up ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ DO ikl=1,knonv isl = isnoSV(ikl) TSurf0(ikl) = TsisSV(ikl,isl) Elem_A = dtC_sv(ikl,isl)*mu__dz(ikl,isl) Elem_C = 0. Diag_A(ikl,isl) = -Elem_A *Implic Diag_C(ikl,isl) = 0. Diag_B(ikl,isl) = 1.0d+0 -Diag_A(ikl,isl) Diag_B(ikl,isl) = Diag_B(ikl,isl) & + dtC_sv(ikl,isl) * (dIRsdT(ikl) & ! Upw. Sol IR +dSdTSV(ikl) & ! HS/Surf.Contr. +dLdTSV(ikl)) ! HL/Surf.Contr. Term_D(ikl,isl) = Explic *Elem_A *TsisSV(ikl,isl-1) & +(1.0d+0 -Explic *Elem_A)*TsisSV(ikl,isl) Term_D(ikl,isl) = Term_D(ikl,isl) & + dtC_sv(ikl,isl) * (sol_SV(ikl) *SoSosv(ikl) & ! Absorbed *(sEX_sv(ikl,isl+1) & ! Solar -sEX_sv(ikl,isl )) & ! + IRd_SV(ikl)*Eso_sv(ikl) & ! Down Atm IR +IRs__D(ikl) & ! Upw. Sol IR +HS___D(ikl) & ! HS/Atmo.Contr. +HL___D(ikl) )! HL/Atmo.Contr. END DO ! +--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 ikl= 1,knonv DO isl=-nsol+1,min(nsno,isnoSV(ikl)+1) 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 TsisSV(ikl,-nsol) = Term_D(ikl,-nsol)/Aux__P(ikl,-nsol) END DO DO ikl= 1,knonv DO isl=-nsol+1,min(nsno,isnoSV(ikl)+1) TsisSV(ikl,isl) =(Term_D(ikl,isl) & -Diag_A(ikl,isl) *TsisSV(ikl,isl-1)) & /Aux__P(ikl,isl) END DO END DO ! +--Backward Sweep ! + ^^^^^^^^^^^^^^ DO ikl= 1,knonv DO isl=min(nsno-1,isnoSV(ikl)+1),-nsol,-1 TsisSV(ikl,isl) = Aux__Q(ikl,isl) *TsisSV(ikl,isl+1) & +TsisSV(ikl,isl) if(isl==0.and.isnoSV(ikl)==0) then TsisSV(ikl,isl) = min(TaT_SV(ikl)+30,TsisSV(ikl,isl)) TsisSV(ikl,isl) = max(TaT_SV(ikl)-30,TsisSV(ikl,isl)) ! #EU TsisSV(ikl,isl) = max(TaT_SV(ikl)-15.,TsisSV(ikl,isl)) ! !XF 18/11/2018 to avoid ST reaching 70�C!! ! !It is an error compensation but does not work over tundra endif END DO END DO ! +--Temperature Limits (avoids problems in case of no Snow Layers) ! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ DO ikl= 1,knonv isl = isnoSV(ikl) dTSurf = TsisSV(ikl,isl) - TSurf0(ikl) TsisSV(ikl,isl) = TSurf0(ikl) + sign(1.,dTSurf) & ! 180.0 dgC/hr * min(abs(dTSurf),5.e-2*dt__SV) ! =0.05 dgC/s END DO DO ikl= 1,knonv DO isl=min(nsno,isnoSV(ikl)+1),1 ,-1 TsisSV(ikl,isl) = max(Ts_Min, TsisSV(ikl,isl)) TsisSV(ikl,isl) = min(Ts_Max, TsisSV(ikl,isl)) END DO END DO ! +--Update Surface Fluxes ! + ======================== DO ikl= 1,knonv isl = isnoSV(ikl) IRs_SV(ikl) = IRs__D(ikl) & ! - dIRsdT(ikl) * TsisSV(ikl,isl) ! HSs_sv(ikl) = HS___D(ikl) & ! Sensible Heat - dSdTSV(ikl) * TsisSV(ikl,isl) ! Downward > 0 HLs_sv(ikl) = HL___D(ikl) & ! Latent Heat - dLdTSV(ikl) * TsisSV(ikl,isl) ! Downward > 0 END DO return end subroutine sisvat_tso