source: LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_tso.F @ 4629

Last change on this file since 4629 was 3900, checked in by evignon, 3 years ago

Commission de la nouvelle interface LMDZ-SISVAT
la prochaine commission consistera a supprimer l'ancien repertoire sisvat
et a faire un peu de nettoyage.

File size: 29.3 KB
Line 
1 
2 
3 
4 
5      subroutine SISVAT_TSo
6! #e1.                     (ETSo_0,ETSo_1,ETSo_d)
7 
8C +------------------------------------------------------------------------+
9C | MAR          SISVAT_TSo                                06-10-2020  MAR |
10C |   SubRoutine SISVAT_TSo computes the Soil/Snow Energy Balance          |
11C +------------------------------------------------------------------------+
12C |                                                                        |
13C |   PARAMETERS:  knonv: Total Number of columns =                        |
14C |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
15C |                     X       Number of Mosaic Cell per grid box         |
16C |                                                                        |
17C |   INPUT:   isotSV   = 0,...,11:   Soil       Type                      |
18C |   ^^^^^               0:          Water, Solid or Liquid               |
19C |            isnoSV   = total Nb of Ice/Snow Layers                      |
20C |            dQa_SV   = Limitation of  Water Vapor  Turbulent Flux       |
21C |                                                                        |
22C |   INPUT:   sol_SV   : Downward Solar Radiation                  [W/m2] |
23C |   ^^^^^    IRd_SV   : Surface Downward  Longwave   Radiation    [W/m2] |
24C |            za__SV   : SBL Top    Height                            [m] |
25C |            VV__SV   : SBL Top    Wind Speed                      [m/s] |
26C |            TaT_SV   : SBL Top    Temperature                       [K] |
27C |            rhT_SV   : SBL Top    Air  Density                  [kg/m3] |
28C |            QaT_SV   : SBL Top    Specific  Humidity            [kg/kg] |
29C |            LSdzsv   : Vertical   Discretization Factor             [-] |
30C |                     =    1. Soil                                       |
31C |                     = 1000. Ocean                                      |
32C |            dzsnSV   : Snow Layer Thickness                         [m] |
33C |            ro__SV   : Snow/Soil  Volumic Mass                  [kg/m3] |
34C |            eta_SV   : Soil Water Content                       [m3/m3] |
35C |            dt__SV   : Time Step                                    [s] |
36C |                                                                        |
37C |            SoSosv   : Absorbed Solar Radiation by Surfac.(Normaliz)[-] |
38C |            Eso_sv   : Soil+Snow       Emissivity                   [-] |
39C |            rah_sv   : Aerodynamic Resistance for Heat            [s/m] |
40C |            Lx_H2O   : Latent Heat of Vaporization/Sublimation   [J/kg] |
41C |            sEX_sv   : Verticaly Integrated Extinction Coefficient  [-] |
42C |                                                                        |
43C |   INPUT /  TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
44C |   OUTPUT:           & Snow     Temperatures (layers  1,2,...,nsno) [K] |
45C |   ^^^^^^                                                               |
46C |                                                                        |
47C |   OUTPUT:  IRs_SV   : Soil      IR Radiation                    [W/m2] |
48C |   ^^^^^^   HSs_sv   : Sensible  Heat Flux                       [W/m2] |
49C |            HLs_sv   : Latent    Heat Flux                       [W/m2] |
50C |            ETSo_0   : Snow/Soil Energy Power, before Forcing    [W/m2] |
51C |            ETSo_1   : Snow/Soil Energy Power, after  Forcing    [W/m2] |
52C |            ETSo_d   : Snow/Soil Energy Power         Forcing    [W/m2] |
53C |                                                                        |
54C |   Internal Variables:                                                  |
55C |   ^^^^^^^^^^^^^^^^^^                                                   |
56C |                                                                        |
57C |   METHOD: NO   Skin Surface Temperature                                |
58C |   ^^^^^^  Semi-Implicit Crank Nicholson Scheme                         |
59C |                                                                        |
60C | # OPTIONS: #E0: Energy Budget Verification                             |
61C | # ^^^^^^^  #kd: KDsvat Option:NO Flux  Limitor     on HL               |
62C | #          #KD: KDsvat Option:Explicit Formulation of HL               |
63C | #          #NC: OUTPUT for Stand Alone NetCDF File                     |
64C |                                                                        |
65C +------------------------------------------------------------------------+
66 
67 
68 
69 
70C +--Global Variables
71C +  ================
72 
73      use VARphy
74      use VAR_SV
75      use VARdSV
76      use VARxSV
77      use VARySV
78      use VARtSV
79      use VAR0SV
80     
81
82      IMPLICIT NONE
83
84 
85C +--OUTPUT
86C +  ------
87 
88! #e1 real      ETSo_0(knonv)                 ! Soil/Snow Power, before Forcing
89! #e1 real      ETSo_1(knonv)                 ! Soil/Snow Power, after  Forcing
90! #e1 real      ETSo_d(knonv)                 ! Soil/Snow Power, Forcing
91 
92 
93C +--Internal Variables
94C +  ==================
95 
96      integer   ikl   ,isl   ,jsl   ,ist      !
97      integer   ist__s,ist__w                 ! Soil/Water  Body Identifier
98      integer   islsgn                        ! Soil/Snow Surfac.Identifier
99      real      eps__3                        ! Arbitrary    Low Number
100      real      etaMid,psiMid                 ! Layer Interface's Humidity
101      real      mu_eta                        !     Soil thermal Conductivity
102      real      mu_exp                        ! arg Soil thermal Conductivity
103      real      mu_min                        ! Min Soil thermal Conductivity
104      real      mu_max                        ! Max Soil thermal Conductivity
105      real      mu_sno(knonv),mu_aux          !     Snow thermal Conductivity
106      real      mu__dz(knonv,-nsol:nsno+1)    ! mu_(eta,sno)   / dz
107      real      dtC_sv(knonv,-nsol:nsno)      ! dt      / C
108      real      IRs__D(knonv)                 ! UpwardIR Previous Iter.Contr.
109      real      dIRsdT(knonv)                 ! UpwardIR           T Derivat.
110      real      f_HSHL(knonv)                 ! Factor common to HS and HL
111      real      dRidTs(knonv)                 ! d(Rib)/d(Ts)
112      real      HS___D(knonv)                 ! Sensible Heat Flux Atm.Contr.
113      real      f___HL(knonv)                 !
114      real      HL___D(knonv)                 ! Latent   Heat Flux Atm.Contr.
115      REAL      TSurf0(knonv),dTSurf          ! Previous Surface Temperature
116      real      qsatsg(knonv) !,den_qs,arg_qs ! Soil   Saturat. Spec. Humidity
117      real      dqs_dT(knonv)                 ! d(qsatsg)/dTv
118      real      Psi(   knonv)                 ! 1st Soil Layer Water Potential
119      real      RHuSol(knonv)                 ! Soil Surface Relative Humidity
120      real      etaSol                        ! Soil Surface          Humidity
121      real      d__eta                        ! Soil Surface Humidity Increm.
122      real      Elem_A,Elem_C                 !   Diagonal Coefficients
123      real      Diag_A(knonv,-nsol:nsno)      ! A Diagonal
124      real      Diag_B(knonv,-nsol:nsno)      ! B Diagonal
125      real      Diag_C(knonv,-nsol:nsno)      ! C Diagonal
126      real      Term_D(knonv,-nsol:nsno)      !   Independant Term
127      real      Aux__P(knonv,-nsol:nsno)      ! P Auxiliary Variable
128      real      Aux__Q(knonv,-nsol:nsno)      ! Q Auxiliary Variable
129      real      Ts_Min,Ts_Max                 ! Temperature Limits
130! #e1 real      Exist0                        ! Existing Layer Switch
131      real      psat_wat, psat_ice, sp        ! computation of qsat
132
133      integer   nt_srf,it_srf,itEuBk          ! HL: Surface Scheme
134      parameter(nt_srf=10)                     ! 10 before
135      real      agpsrf,xgpsrf,dt_srf,dt_ver   !
136      real      etaBAK(knonv)                 !
137      real      etaNEW(knonv)                 !
138      real      etEuBk(knonv)                 !
139      real      fac_dt(knonv),faceta(knonv)   !
140      real      PsiArg(knonv),SHuSol(knonv)   !
141 
142
143 
144C +--Internal DATA
145C +  =============
146 
147      data      eps__3 /   1.e-3   /          ! Arbitrary    Low Number
148      data      mu_exp /  -0.4343  /          ! Soil Thermal Conductivity
149      data      mu_min /   0.172   /          ! Min Soil Thermal Conductivity
150      data      mu_max /   2.000   /          ! Max Soil Thermal Conductivity
151      data      Ts_Min / 175.      /          ! Temperature            Minimum
152      data      Ts_Max / 300.      /          ! Temperature Acceptable Maximum
153C +                                           ! including   Snow Melt  Energy
154 
155C +-- Initilialisation of local arrays
156C +   ================================
157        DO ikl=1,knonv
158
159          mu_sno(ikl)=0.
160          mu__dz(ikl,:)=0.   
161          dtC_sv(ikl,:)=0.     
162          IRs__D(ikl)=0.                 
163          dIRsdT(ikl)=0.                 
164          f_HSHL(ikl)=0.                 
165          dRidTs(ikl)=0.                 
166          HS___D(ikl)=0.               
167          f___HL(ikl)=0.             
168          HL___D(ikl)=0.           
169          TSurf0(ikl)=0.     
170          qsatsg(ikl)=0.
171          dqs_dT(ikl)=0.                 
172          Psi(ikl)=0.               
173          RHuSol(ikl)=0.                 
174          Diag_A(ikl,:)=0.     
175          Diag_B(ikl,:)=0.     
176          Diag_C(ikl,:)=0.     
177          Term_D(ikl,:)=0.     
178          Aux__P(ikl,:)=0.     
179          Aux__Q(ikl,:)=0.     
180          etaBAK(ikl)=0.               
181          etaNEW(ikl)=0.               
182          etEuBk(ikl)=0.                 
183          fac_dt(ikl)=0.
184          faceta(ikl)=0. 
185          PsiArg(ikl)=0.
186          SHuSol(ikl)=0. 
187
188        END DO
189
190       
191
192C +--Heat Conduction Coefficient (zero in the Layers over the highest one)
193C +  ===========================
194C +                             ---------------- isl    eta_SV, rho C (isl)
195C +
196C +--Soil                       ++++++++++++++++        etaMid,    mu (isl)
197C +  ----
198C +                             ---------------- isl-1  eta_SV, rho C (isl-1)
199           isl=-nsol
200        DO ikl=1,knonv
201
202          mu__dz(ikl,isl) = 0.
203 
204          dtC_sv(ikl,isl) = dtz_SV2(isl)  * dt__SV        ! dt / (dz X rho C)
205     .                   /((rocsSV(isotSV(ikl))           ! [s / (m.J/m3/K)]
206     .                     +rcwdSV*eta_SV(ikl,isl))       !
207     .                     *LSdzsv(ikl)            )      !
208        END DO
209      DO   isl=-nsol+1,0
210        DO ikl=1,knonv
211          ist    =      isotSV(ikl)                       ! Soil Type
212          ist__s =  min(ist, 1)                           ! 1 => Soil
213          ist__w =  1 - ist__s                            ! 1 => Water Body
214 
215          etaMid = 0.5*(dz_dSV(isl-1)*eta_SV(ikl,isl-1)   ! eta at layers
216     .                 +dz_dSV(isl)  *eta_SV(ikl,isl)  )  !     interface
217     .                 /dzmiSV(isl)                       ! LSdzsv implicit !
218          etaMid =  max(etaMid,epsi)
219          psiMid =      psidSV(ist)
220     .                *(etadSV(ist)/etaMid)**bCHdSV(ist)
221          mu_eta =      3.82      *(psiMid)**mu_exp       ! Soil Thermal
222          mu_eta =  min(max(mu_eta, mu_min), mu_max)      ! Conductivity
223C +                                                       ! DR97 eq.3.31
224          mu_eta =  ist__s *mu_eta +ist__w * vK_dSV       ! Water Bodies
225C +                                                       ! Correction
226          mu__dz(ikl,isl) = mu_eta/(dzmiSV(isl)           !
227     .                             *LSdzsv(ikl))          !
228 
229          dtC_sv(ikl,isl) = dtz_SV2(isl)* dt__SV          ! dt / (dz X rho C)
230     .                   /((rocsSV(isotSV(ikl))           !
231     .                     +rcwdSV*eta_SV(ikl,isl))       !
232     .                     *LSdzsv(ikl)            )      !
233        END DO
234      END DO
235 
236 
237C +--Soil/Snow Interface
238C +  -------------------
239 
240C +--Soil Contribution
241C +  ^^^^^^^^^^^^^^^^^
242           isl=1
243        DO ikl=1,knonv
244          ist    =      isotSV(ikl)                       ! Soil Type
245          ist__s =  min(ist, 1)                           ! 1 => Soil
246          ist__w =  1 - ist__s                            ! 1 => Water Body
247          psiMid =      psidSV(ist)                       ! Snow => Saturation
248          mu_eta =      3.82      *(psiMid)**mu_exp       ! Soil Thermal
249          mu_eta =  min(max(mu_eta, mu_min), mu_max)      ! Conductivity
250C +                                                       ! DR97 eq.3.31
251          mu_eta =  ist__s *mu_eta +ist__w * vK_dSV       ! Water Bodies
252 
253C +--Snow Contribution
254C +  ^^^^^^^^^^^^^^^^^
255          mu_sno(ikl) =  CdidSV                           !
256     .                 *(ro__SV(ikl,isl) /ro_Wat) ** 1.88 !
257          mu_sno(ikl) =          max(epsi,mu_sno(ikl))    !
258C +...    mu_sno :  Snow Heat Conductivity Coefficient [Wm/K]
259C +                 (Yen 1981, CRREL Rep., 81-10)
260 
261C +--Combined Heat Conductivity
262C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
263          mu__dz(ikl,isl) = 2./(dzsnSV(ikl,isl  )         ! Combined Heat
264     .                         /mu_sno(ikl)               ! Conductivity
265     .                         +LSdzsv(ikl)               !
266     .                         *dz_dSV(    isl-1)/mu_eta) ! Coefficient
267 
268C +--Inverted Heat Capacity
269C +  ^^^^^^^^^^^^^^^^^^^^^^
270          dtC_sv(ikl,isl) = dt__SV/max(epsi,              ! dt / (dz X rho C)
271     .    dzsnSV(ikl,isl) * ro__SV(ikl,isl) *Cn_dSV)      !
272        END DO
273 
274 
275C +--Snow
276C +  ----
277 
278      DO ikl=1,knonv
279      DO   isl=1,min(nsno,isnoSV(ikl)+1)
280          ro__SV(ikl,isl) =                               !
281     .                   ro__SV(ikl ,isl)                 !
282     .       * max(0,min(isnoSV(ikl)-isl+1,1))            !
283
284        END DO
285      END DO
286 
287      DO ikl=1,knonv
288      DO   isl=1,min(nsno,isnoSV(ikl)+1)
289 
290C +--Combined Heat Conductivity
291C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
292          mu_aux      =  CdidSV                           !
293     .                 *(ro__SV(ikl,isl) /ro_Wat) ** 1.88 !
294          mu__dz(ikl,isl) =                               !
295     .      2.                        *mu_aux*mu_sno(ikl) ! Combined Heat
296     .     /max(epsi,dzsnSV(ikl,isl  )*mu_sno(ikl)        ! Conductivity
297     .              +dzsnSV(ikl,isl-1)*mu_aux     )       ! For upper Layer
298          mu_sno(ikl)     =            mu_aux             !
299 
300C +--Inverted Heat Capacity
301C +  ^^^^^^^^^^^^^^^^^^^^^^
302          dtC_sv(ikl,isl) = dt__SV/max(eps__3,            ! dt / (dz X rho C)
303     .    dzsnSV(ikl,isl) * ro__SV(ikl,isl) *Cn_dSV)      !
304        END DO
305      END DO
306 
307 
308C +--Uppermost Effective Layer: NO conduction
309C +  ----------------------------------------
310 
311        DO ikl=1,knonv
312          mu__dz(ikl,isnoSV(ikl)+1) = 0.0
313        END DO
314 
315 
316C +--Energy Budget (IN)
317C +  ==================
318 
319! #e1   DO ikl=1,knonv
320! #e1     ETSo_0(ikl) = 0.
321! #e1   END DO
322! #e1 DO   isl= -nsol,nsno
323! #e1   DO ikl=1,knonv
324! #e1     Exist0      = isl -           isnoSV(ikl)
325! #e1     Exist0      = 1.  - max(zero,min(unun,Exist0))
326! #e1     ETSo_0(ikl) = ETSo_0(ikl)
327! #e1.                +(TsisSV(ikl,isl)-TfSnow)*Exist0
328! #e1.                                 /dtC_sv(ikl,isl)
329! #e1   END DO
330! #e1 END DO
331 
332 
333C +--Tridiagonal Elimination: Set Up
334C +  ===============================
335 
336C +--Soil/Snow Interior
337C +  ^^^^^^^^^^^^^^^^^^
338      DO ikl=1,knonv
339      DO   isl=-nsol+1,min(nsno-1,isnoSV(ikl)+1)
340
341          Elem_A          =  dtC_sv(ikl,isl)         *mu__dz(ikl,isl)
342          Elem_C          =  dtC_sv(ikl,isl)         *mu__dz(ikl,isl+1)
343          Diag_A(ikl,isl) = -Elem_A  *Implic
344          Diag_C(ikl,isl) = -Elem_C  *Implic
345          Diag_B(ikl,isl) =  1.0d+0  -Diag_A(ikl,isl)-Diag_C(ikl,isl)
346          Term_D(ikl,isl) =  Explic *(Elem_A         *TsisSV(ikl,isl-1)
347     .                               +Elem_C         *TsisSV(ikl,isl+1))
348     .             +(1.0d+0 -Explic *(Elem_A+Elem_C))*TsisSV(ikl,isl)
349     .  + dtC_sv(ikl,isl)           * sol_SV(ikl)    *SoSosv(ikl)
350     .                                               *(sEX_sv(ikl,isl+1)
351     .                                               -sEX_sv(ikl,isl  ))
352        END DO
353      END DO
354 
355C +--Soil  lowest Layer
356C +  ^^^^^^^^^^^^^^^^^^
357           isl= -nsol
358        DO ikl=1,knonv
359          Elem_A          =  0.
360          Elem_C          =  dtC_sv(ikl,isl)         *mu__dz(ikl,isl+1)
361          Diag_A(ikl,isl) =  0.
362          Diag_C(ikl,isl) = -Elem_C  *Implic
363          Diag_B(ikl,isl) =  1.0d+0  -Diag_A(ikl,isl)-Diag_C(ikl,isl)
364          Term_D(ikl,isl) =  Explic * Elem_C         *TsisSV(ikl,isl+1)
365     .             +(1.0d+0 -Explic * Elem_C)        *TsisSV(ikl,isl)
366     .  + dtC_sv(ikl,isl)           * sol_SV(ikl)    *SoSosv(ikl)
367     .                                              *(sEX_sv(ikl,isl+1)
368     .                                               -sEX_sv(ikl,isl  ))
369        END DO
370 
371C +--Snow highest Layer (dummy!)
372C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
373
374        !EV!isl=  min(isnoSV(1)+1,nsno)
375
376        DO ikl=1,knonv
377! EV try to calculate isl at the ikl grid point
378          isl=  min(isnoSV(ikl)+1,nsno)
379
380          Elem_A          =  dtC_sv(ikl,isl)  *mu__dz(ikl,isl)
381          Elem_C          =  0.
382          Diag_A(ikl,isl) = -Elem_A  *Implic
383          Diag_C(ikl,isl) =  0.
384          Diag_B(ikl,isl) =  1.0d+0  -Diag_A(ikl,isl)
385          Term_D(ikl,isl) =  Explic * Elem_A  *TsisSV(ikl,isl-1)
386     .             +(1.0d+0 -Explic * Elem_A) *TsisSV(ikl,isl)
387     .  + dtC_sv(ikl,isl) * (sol_SV(ikl)      *SoSosv(ikl)
388     .                                       *(sEX_sv(ikl,isl+1)
389     .                                        -sEX_sv(ikl,isl  )))
390        END DO
391 
392C +--Surface: UPwardIR Heat Flux
393C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
394        DO ikl=1,knonv
395          isl         = isnoSV(ikl)
396          dIRsdT(ikl) = Eso_sv(ikl)* StefBo          * 4.      ! - d(IR)/d(T)
397     .                             * TsisSV(ikl,isl)           !
398     .                             * TsisSV(ikl,isl)           !
399     .                             * TsisSV(ikl,isl)           !
400          IRs__D(ikl) = dIRsdT(ikl)* TsisSV(ikl,isl) * 0.75    !
401 
402C +--Surface: Richardson Number:   T Derivative
403C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
404c #RC     dRidTs(ikl) =-gravit      *    za__SV(ikl)
405c #RC.                /(TaT_SV(ikl) *    VV__SV(ikl)
406c #RC.                              *    VV__SV(ikl))
407 
408C +--Surface: Turbulent Heat Flux: Factors
409C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
410          f_HSHL(ikl) = rhT_SV(ikl)  /    rah_sv(ikl)           ! to  HS, HL
411          f___HL(ikl) = f_HSHL(ikl) *    Lx_H2O(ikl)
412 
413C +--Surface: Sensible  Heat Flux: T Derivative
414C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
415          dSdTSV(ikl) = f_HSHL(ikl) *    Cp                    !#- d(HS)/d(T)
416c #RC.         *(1.0  -(TsisSV(ikl,isl) -TaT_SV(ikl))          !#Richardson
417c #RC.         * dRidTs(ikl)*dFh_sv(ikl)/rah_sv(ikl))          ! Nb. Correct.
418          HS___D(ikl) = dSdTSV(ikl) *    TaT_SV(ikl)           !
419 
420C +--Surface: Latent    Heat Flux: Saturation Specific Humidity
421C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
422c         den_qs      =         TsisSV(ikl,isl)- 35.8          !
423c         arg_qs      = 17.27 *(TsisSV(ikl,isl)-273.16)        !
424c    .                                   / den_qs              !
425c         qsatsg(ikl) = .0038 *        exp(arg_qs)             !
426!          sp = (pst_SV(ikl) + ptopSV) * 10.
427
428          !sp=ps__SV(ikl)
429          ! Etienne: in the formula herebelow sp should be in hPa, not
430          ! in Pa so I divide by 100.
431          sp=ps__SV(ikl)/100.
432          psat_ice = 6.1070 * exp(6150. *(1./273.16 -
433     .                                              1./TsisSV(ikl,isl)))
434
435          psat_wat = 6.1078 * exp (5.138*log(273.16   /TsisSV(ikl,isl)))
436     .             * exp (6827.*(1.         /273.16-1./TsisSV(ikl,isl)))
437
438          if(TsisSV(ikl,isl)<=273.16) then
439            qsatsg(ikl) = 0.622 * psat_ice / (sp - 0.378 * psat_ice)
440          else
441            qsatsg(ikl) = 0.622 * psat_wat / (sp - 0.378 * psat_wat)
442          endif
443          QsT_SV(ikl)=qsatsg(ikl)
444
445c         dqs_dT(ikl) = qsatsg(ikl)* 4099.2   /(den_qs *den_qs)!
446          fac_dt(ikl) = f_HSHL(ikl)/(ro_Wat   * dz_dSV(0))     !
447        END DO
448
449
450 
451C +--Surface: Latent    Heat Flux: Surface    Relative Humidity
452C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
453              xgpsrf       =   1.05                            !
454              agpsrf       = dt__SV*(   1.0-xgpsrf        )    !
455     .                             /(   1.0-xgpsrf**nt_srf)    !
456              dt_srf       = agpsrf                            !
457              dt_ver       = 0.               
458
459            DO ikl=1,knonv
460              isl          =          isnoSV(ikl)             
461              ist          = max(0,isotSV(ikl)-100*isnoSV(ikl))! 0 if    H2O                         
462              ist__s       = min(1,ist)                                                                             
463              etaBAK(ikl)  = max(epsi,eta_SV(ikl ,isl))        !
464              etaNEW(ikl)  =          etaBAK(ikl)              !
465              etEuBk(ikl)  =          etaNEW(ikl)              !
466            END DO     
467
468        if(ist__s==1) then ! to reduce computer time                                                 
469                                          !
470        DO it_srf=1,nt_srf                                     !
471              dt_ver       = dt_ver     +dt_srf                !
472            DO ikl=1,knonv                                     !
473              faceta(ikl)  = fac_dt(ikl)*dt_srf                !
474c #VX         faceta(ikl)  = faceta(ikl)                       !
475c #VX.                  /(1.+faceta(ikl)*dQa_SV(ikl))          !    Limitation
476                                                               ! by Atm.Conten
477c #??.        *max(0,sign(1.,qsatsg(ikl)-QaT_SV(ikl))))        ! NO Limitation
478                                                               ! of Downw.Flux
479            END DO                                             !
480          DO itEuBk=1,2                                        !
481            DO ikl=1,knonv
482              ist    = max(0,isotSV(ikl)-100*isnoSV(ikl))      ! 0 if    H2O
483                                                               !
484              Psi(ikl) =                                       !
485     .                psidSV(ist)                              ! DR97, Eqn 3.34
486     .              *(etadSV(ist)                              !
487     .           /max(etEuBk(ikl),epsi))                       !
488     .              **bCHdSV(ist)                              !
489              PsiArg(ikl) = 7.2E-5*Psi(ikl)                    !
490              RHuSol(ikl) =   exp(-min(0.,PsiArg(ikl)))    !
491              SHuSol(ikl) =     qsatsg(ikl)  *RHuSol(ikl)      ! DR97, Eqn 3.15
492              etEuBk(ikl) =                                    !
493     .       (etaNEW(ikl) + faceta(ikl)*(QaT_SV(ikl)           !
494     .                                  -SHuSol(ikl)           !
495     .                    *(1.          -bCHdSV(ist)           !
496     .                                  *PsiArg(ikl))       )) !
497     .      /(1.          + faceta(ikl)* SHuSol(ikl)           !
498     .                                  *bCHdSV(ist)           !
499     .                                  *PsiArg(ikl)           !
500     .                                  /etaNEW(ikl))          !
501              etEuBk(ikl) = etEuBk(ikl)          !
502c    .                                 /(Ro_Wat*dz_dSV(0))     !
503     .                    * dt_srf     /(Ro_Wat*dz_dSV(0))     !
504cXF 15/05/2017 BUG
505            END DO                                             !
506          END DO                                               !
507            DO ikl=1,knonv                                     !
508              etaNEW(ikl) =  max(etEuBk(ikl),epsi)             !
509            END DO                                             !
510              dt_srf      =      dt_srf         * xgpsrf       !
511        END DO         
512
513
514        endif                                       !
515 
516C +--Surface: Latent    Heat Flux: Soil/Water Surface Contributions
517C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
518        DO ikl=1,knonv                                         !
519          isl        =  isnoSV(ikl)                            !
520          ist   = max(0,isotSV(ikl)-100*isnoSV(ikl))           ! 0 if    H2O
521          ist__s= min(1,ist)                                   ! 1 if no H2O
522          ist__w=     1-ist__s                                 ! 1 if    H2O
523          d__eta     =  eta_SV(ikl,isl)-etaNEW(ikl)            !
524          ! latent heat flux computation
525          HL___D(ikl)=( ist__s *ro_Wat *dz_dSV(0)              ! Soil Contrib.
526     .                *(etaNEW(ikl)    -etaBAK(ikl)) / dt__SV  !
527     .                 +ist__w         *f_HSHL(ikl)            ! H2O  Contrib.
528     .                *(QaT_SV(ikl)    - qsatsg(ikl))        ) !
529     .                * Lx_H2O(ikl)                            ! common factor
530 
531c #DL     RHuSol(ikl) =(QaT_SV(ikl)                            !
532c #DL.                 -HL___D(ikl)    / f___HL(ikl))          !
533c #DL.                / qsatsg(ikl)                            !
534 
535C +--Surface: Latent    Heat Flux: T Derivative
536C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
537          dLdTSV(ikl) = 0.
538c #DL     dLdTSV(ikl) = f___HL(ikl) * RHuSol(ikl) *dqs_dT(ikl) ! - d(HL)/d(T)
539c #DL     HL___D(ikl) = HL___D(ikl)                            !
540c #DL.                 +dLdTSV(ikl) * TsisSV(ikl,isl)          !
541        END DO                                                 !
542 
543C +--Surface: Tridiagonal Matrix Set Up
544C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
545        DO ikl=1,knonv
546          isl             =  isnoSV(ikl)
547          TSurf0(ikl)     =  TsisSV(ikl,isl)
548
549          Elem_A          =  dtC_sv(ikl,isl)*mu__dz(ikl,isl)
550          Elem_C          =  0.
551          Diag_A(ikl,isl) = -Elem_A *Implic
552          Diag_C(ikl,isl) =  0.
553          Diag_B(ikl,isl) =  1.0d+0 -Diag_A(ikl,isl)
554          Diag_B(ikl,isl) =  Diag_B(ikl,isl)
555     .  + dtC_sv(ikl,isl) * (dIRsdT(ikl)                       ! Upw. Sol IR
556     .                      +dSdTSV(ikl)                       ! HS/Surf.Contr.
557     .                      +dLdTSV(ikl))                      ! HL/Surf.Contr.
558
559          Term_D(ikl,isl) =  Explic *Elem_A *TsisSV(ikl,isl-1)
560     .             +(1.0d+0 -Explic *Elem_A)*TsisSV(ikl,isl)
561
562
563
564          Term_D(ikl,isl) =  Term_D(ikl,isl)
565     .  + dtC_sv(ikl,isl) * (sol_SV(ikl)    *SoSosv(ikl)       ! Absorbed
566     .                                     *(sEX_sv(ikl,isl+1) ! Solar
567     .                                      -sEX_sv(ikl,isl  ))!
568     .                              +   IRd_SV(ikl)*Eso_sv(ikl) ! Down Atm IR
569     .                                +IRs__D(ikl)             ! Upw. Sol IR
570     .                                +HS___D(ikl)             ! HS/Atmo.Contr.
571     .                                +HL___D(ikl)            )! HL/Atmo.Contr.
572     
573        END DO
574 
575 
576C +--Tridiagonal Elimination
577C +  =======================
578 
579C +--Forward  Sweep
580C +  ^^^^^^^^^^^^^^
581        DO ikl=  1,knonv
582          Aux__P(ikl,-nsol) = Diag_B(ikl,-nsol)
583          Aux__Q(ikl,-nsol) =-Diag_C(ikl,-nsol)/Aux__P(ikl,-nsol)
584        END DO
585 
586        DO ikl=      1,knonv
587
588        DO   isl=-nsol+1,min(nsno,isnoSV(ikl)+1)
589          Aux__P(ikl,isl)   = Diag_A(ikl,isl)  *Aux__Q(ikl,isl-1)
590     .                       +Diag_B(ikl,isl)
591          Aux__Q(ikl,isl)   =-Diag_C(ikl,isl)  /Aux__P(ikl,isl)
592        END DO
593        END DO
594 
595        DO ikl=      1,knonv
596          TsisSV(ikl,-nsol) = Term_D(ikl,-nsol)/Aux__P(ikl,-nsol)
597        END DO
598 
599        DO ikl=      1,knonv
600        DO   isl=-nsol+1,min(nsno,isnoSV(ikl)+1)
601          TsisSV(ikl,isl)   =(Term_D(ikl,isl)
602     .                       -Diag_A(ikl,isl)  *TsisSV(ikl,isl-1))
603     .                                         /Aux__P(ikl,isl)
604
605 
606        END DO
607        END DO
608 
609C +--Backward Sweep
610C +  ^^^^^^^^^^^^^^
611        DO ikl=     1,knonv
612        DO   isl=min(nsno-1,isnoSV(ikl)+1),-nsol,-1
613           
614
615          TsisSV(ikl,isl)   = Aux__Q(ikl,isl)  *TsisSV(ikl,isl+1)
616     .                                         +TsisSV(ikl,isl)
617          if(isl==0.and.isnoSV(ikl)==0) then
618 
619           TsisSV(ikl,isl)  = min(TaT_SV(ikl)+30,TsisSV(ikl,isl))
620           TsisSV(ikl,isl)  = max(TaT_SV(ikl)-30,TsisSV(ikl,isl))
621
622
623c #EU      TsisSV(ikl,isl)  = max(TaT_SV(ikl)-15.,TsisSV(ikl,isl))
624
625          !XF 18/11/2018 to avoid ST reaching 70�C!!
626          !It is an error compensation but does not work over tundra
627 
628          endif
629
630   
631
632        END DO
633       
634      END DO
635
636
637 
638C +--Temperature Limits (avoids problems in case of no Snow Layers)
639C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
640        DO ikl=     1,knonv
641           isl              = isnoSV(ikl)
642
643           dTSurf            = TsisSV(ikl,isl) -     TSurf0(ikl)
644          TsisSV(ikl,isl)   = TSurf0(ikl) + sign(1.,dTSurf) ! 180.0 dgC/hr
645     .              * min(abs(dTSurf),5.e-2*dt__SV)         ! =0.05 dgC/s
646
647
648
649        END DO
650
651        DO ikl=     1,knonv
652        DO   isl=min(nsno,isnoSV(ikl)+1),1      ,-1
653          TsisSV(ikl,isl)   = max(Ts_Min,       TsisSV(ikl,isl))
654          TsisSV(ikl,isl)   = min(Ts_Max,       TsisSV(ikl,isl))
655        END DO
656
657        END DO
658 
659C +--Update Surface    Fluxes
660C +  ========================
661       
662
663
664        DO ikl=      1,knonv
665          isl         = isnoSV(ikl)
666          IRs_SV(ikl) = IRs__D(ikl)                          !
667     .                - dIRsdT(ikl) * TsisSV(ikl,isl)        !
668          HSs_sv(ikl) = HS___D(ikl)                          ! Sensible Heat
669     .                - dSdTSV(ikl) * TsisSV(ikl,isl)        ! Downward > 0
670          HLs_sv(ikl) = HL___D(ikl)                          ! Latent   Heat
671     .                - dLdTSV(ikl) * TsisSV(ikl,isl)        ! Downward > 0
672        END DO
673
674      return
675      end
Note: See TracBrowser for help on using the repository browser.