source: LMDZ6/branches/Amaury_dev/libf/phylmd/inlandsis/sisvat_tso.f90 @ 5116

Last change on this file since 5116 was 5116, checked in by abarral, 2 months ago

rename modules properly lmdz_*
move ismin, ismax, minmax into new lmdz_libmath.f90
(lint) uppercase fortran keywords
File size: 26.1 KB
Line 
1
2
3
4
5SUBROUTINE SISVAT_TSo
6  ! #e1.                     (ETSo_0,ETSo_1,ETSo_d)
7
8  ! +------------------------------------------------------------------------+
9  ! | MAR          SISVAT_TSo                                06-10-2020  MAR |
10  ! |   SubRoutine SISVAT_TSo computes the Soil/Snow Energy Balance          |
11  ! +------------------------------------------------------------------------+
12  ! |                                                                        |
13  ! |   PARAMETERS:  knonv: Total Number of columns =                        |
14  ! |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
15  ! |                     X       Number of Mosaic Cell per grid box         |
16  ! |                                                                        |
17  ! |   INPUT:   isotSV   = 0,...,11:   Soil       Type                      |
18  ! |   ^^^^^               0:          Water, Solid or Liquid               |
19  ! |            isnoSV   = total Nb of Ice/Snow Layers                      |
20  ! |            dQa_SV   = Limitation of  Water Vapor  Turbulent Flux       |
21  ! |                                                                        |
22  ! |   INPUT:   sol_SV   : Downward Solar Radiation                  [W/m2] |
23  ! |   ^^^^^    IRd_SV   : Surface Downward  Longwave   Radiation    [W/m2] |
24  ! |            za__SV   : SBL Top    Height                            [m] |
25  ! |            VV__SV   : SBL Top    Wind Speed                      [m/s] |
26  ! |            TaT_SV   : SBL Top    Temperature                       [K] |
27  ! |            rhT_SV   : SBL Top    Air  Density                  [kg/m3] |
28  ! |            QaT_SV   : SBL Top    Specific  Humidity            [kg/kg] |
29  ! |            LSdzsv   : Vertical   Discretization Factor             [-] |
30  ! |                     =    1. Soil                                       |
31  ! |                     = 1000. Ocean                                      |
32  ! |            dzsnSV   : Snow Layer Thickness                         [m] |
33  ! |            ro__SV   : Snow/Soil  Volumic Mass                  [kg/m3] |
34  ! |            eta_SV   : Soil Water Content                       [m3/m3] |
35  ! |            dt__SV   : Time Step                                    [s] |
36  ! |                                                                        |
37  ! |            SoSosv   : Absorbed Solar Radiation by Surfac.(Normaliz)[-] |
38  ! |            Eso_sv   : Soil+Snow       Emissivity                   [-] |
39  ! |            rah_sv   : Aerodynamic Resistance for Heat            [s/m] |
40  ! |            Lx_H2O   : Latent Heat of Vaporization/Sublimation   [J/kg] |
41  ! |            sEX_sv   : Verticaly Integrated Extinction Coefficient  [-] |
42  ! |                                                                        |
43  ! |   INPUT /  TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
44  ! |   OUTPUT:           & Snow     Temperatures (layers  1,2,...,nsno) [K] |
45  ! |   ^^^^^^                                                               |
46  ! |                                                                        |
47  ! |   OUTPUT:  IRs_SV   : Soil      IR Radiation                    [W/m2] |
48  ! |   ^^^^^^   HSs_sv   : Sensible  Heat Flux                       [W/m2] |
49  ! |            HLs_sv   : Latent    Heat Flux                       [W/m2] |
50  ! |            ETSo_0   : Snow/Soil Energy Power, before Forcing    [W/m2] |
51  ! |            ETSo_1   : Snow/Soil Energy Power, after  Forcing    [W/m2] |
52  ! |            ETSo_d   : Snow/Soil Energy Power         Forcing    [W/m2] |
53  ! |                                                                        |
54  ! |   Internal Variables:                                                  |
55  ! |   ^^^^^^^^^^^^^^^^^^                                                   |
56  ! |                                                                        |
57  ! |   METHOD: NO   Skin Surface Temperature                                |
58  ! |   ^^^^^^  Semi-Implicit Crank Nicholson Scheme                         |
59  ! |                                                                        |
60  ! | # OPTIONS: #E0: Energy Budget Verification                             |
61  ! | # ^^^^^^^  #kd: KDsvat Option:NO Flux  Limitor     on HL               |
62  ! | #          #KD: KDsvat Option:Explicit Formulation of HL               |
63  ! | #          #NC: OUTPUT for Stand Alone NetCDF File                     |
64  ! |                                                                        |
65  ! +------------------------------------------------------------------------+
66
67
68
69
70  ! +--Global Variables
71  ! +  ================
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
85  ! +--OUTPUT
86  ! +  ------
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
93  ! +--Internal Variables
94  ! +  ==================
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
144  ! +--Internal DATA
145  ! +  =============
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
153  ! +                                           ! including   Snow Melt  Energy
154
155  ! +-- Initilialisation of local arrays
156  ! +   ================================
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
192  ! +--Heat Conduction Coefficient (zero in the Layers over the highest one)
193  ! +  ===========================
194  ! +                             ---------------- isl    eta_SV, rho C (isl)
195  ! +
196  ! +--Soil                       ++++++++++++++++        etaMid,    mu (isl)
197  ! +  ----
198  ! +                             ---------------- 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
223  ! +                                                       ! DR97 eq.3.31
224      mu_eta =  ist__s *mu_eta +ist__w * vK_dSV       ! Water Bodies
225  ! +                                                       ! 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
237  ! +--Soil/Snow Interface
238  ! +  -------------------
239
240  ! +--Soil Contribution
241  ! +  ^^^^^^^^^^^^^^^^^
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
250  ! +                                                       ! DR97 eq.3.31
251      mu_eta =  ist__s *mu_eta +ist__w * vK_dSV       ! Water Bodies
252
253  ! +--Snow Contribution
254  ! +  ^^^^^^^^^^^^^^^^^
255      mu_sno(ikl) =  CdidSV & !
256            *(ro__SV(ikl,isl) /ro_Wat) ** 1.88 !
257      mu_sno(ikl) =          max(epsi,mu_sno(ikl))    !
258  ! +...    mu_sno :  Snow Heat Conductivity Coefficient [Wm/K]
259  ! +                 (Yen 1981, CRREL Rep., 81-10)
260
261  ! +--Combined Heat Conductivity
262  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
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
268  ! +--Inverted Heat Capacity
269  ! +  ^^^^^^^^^^^^^^^^^^^^^^
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
275  ! +--Snow
276  ! +  ----
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
290  ! +--Combined Heat Conductivity
291  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
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
300  ! +--Inverted Heat Capacity
301  ! +  ^^^^^^^^^^^^^^^^^^^^^^
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
308  ! +--Uppermost Effective Layer: NO conduction
309  ! +  ----------------------------------------
310
311    DO ikl=1,knonv
312      mu__dz(ikl,isnoSV(ikl)+1) = 0.0
313    END DO
314
315
316  ! +--Energy Budget (IN)
317  ! +  ==================
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
333  ! +--Tridiagonal Elimination: Set Up
334  ! +  ===============================
335
336  ! +--Soil/Snow Interior
337  ! +  ^^^^^^^^^^^^^^^^^^
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
355  ! +--Soil  lowest Layer
356  ! +  ^^^^^^^^^^^^^^^^^^
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
371  ! +--Snow highest Layer (dummy!)
372  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
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
392  ! +--Surface: UPwardIR Heat Flux
393  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^
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
402  ! +--Surface: Richardson Number:   T Derivative
403  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
404  ! #RC     dRidTs(ikl) =-gravit      *    za__SV(ikl)
405  ! #RC.                /(TaT_SV(ikl) *    VV__SV(ikl)
406  ! #RC.                              *    VV__SV(ikl))
407
408  ! +--Surface: Turbulent Heat Flux: Factors
409  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
410      f_HSHL(ikl) = rhT_SV(ikl)  /    rah_sv(ikl)           ! to  HS, HL
411      f___HL(ikl) = f_HSHL(ikl) *    Lx_H2O(ikl)
412
413  ! +--Surface: Sensible  Heat Flux: T Derivative
414  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
415      dSdTSV(ikl) = f_HSHL(ikl) *    Cp                    !#- d(HS)/d(T)
416  ! #RC.         *(1.0  -(TsisSV(ikl,isl) -TaT_SV(ikl))          !#Richardson
417  ! #RC.         * dRidTs(ikl)*dFh_sv(ikl)/rah_sv(ikl))          ! Nb. Correct.
418      HS___D(ikl) = dSdTSV(ikl) *    TaT_SV(ikl)           !
419
420  ! +--Surface: Latent    Heat Flux: Saturation Specific Humidity
421  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
422      ! den_qs      =         TsisSV(ikl,isl)- 35.8          !
423      ! arg_qs      = 17.27 *(TsisSV(ikl,isl)-273.16)        !
424  !    .                                   / den_qs              !
425  !         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
445      ! 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
451  ! +--Surface: Latent    Heat Flux: Surface    Relative Humidity
452  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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                !
474  ! #VX         faceta(ikl)  = faceta(ikl)                       !
475  ! #VX.                  /(1.+faceta(ikl)*dQa_SV(ikl))          !    Limitation
476  !                                                              ! by Atm.Conten
477  ! #??.        *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) & !
502  !    .                                 /(Ro_Wat*dz_dSV(0))     !
503                * dt_srf     /(Ro_Wat*dz_dSV(0))     !
504  !XF 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
516  ! +--Surface: Latent    Heat Flux: Soil/Water Surface Contributions
517  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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
531  ! #DL     RHuSol(ikl) =(QaT_SV(ikl)                            !
532  ! #DL.                 -HL___D(ikl)    / f___HL(ikl))          !
533  ! #DL.                / qsatsg(ikl)                            !
534
535  ! +--Surface: Latent    Heat Flux: T Derivative
536  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
537      dLdTSV(ikl) = 0.
538  ! #DL     dLdTSV(ikl) = f___HL(ikl) * RHuSol(ikl) *dqs_dT(ikl) ! - d(HL)/d(T)
539  ! #DL     HL___D(ikl) = HL___D(ikl)                            !
540  ! #DL.                 +dLdTSV(ikl) * TsisSV(ikl,isl)          !
541    END DO                                                 !
542
543  ! +--Surface: Tridiagonal Matrix Set Up
544  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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
576  ! +--Tridiagonal Elimination
577  ! +  =======================
578
579  ! +--Forward  Sweep
580  ! +  ^^^^^^^^^^^^^^
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
609  ! +--Backward Sweep
610  ! +  ^^^^^^^^^^^^^^
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       TsisSV(ikl,isl)  = min(TaT_SV(ikl)+30,TsisSV(ikl,isl))
619       TsisSV(ikl,isl)  = max(TaT_SV(ikl)-30,TsisSV(ikl,isl))
620
621
622  ! #EU      TsisSV(ikl,isl)  = max(TaT_SV(ikl)-15.,TsisSV(ikl,isl))
623
624      !XF 18/11/2018 to avoid ST reaching 70°C!!
625      !It is an error compensation but does not work over tundra
626
627      endif
628
629
630
631    END DO
632
633  END DO
634
635
636
637  ! +--Temperature Limits (avoids problems in case of no Snow Layers)
638  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
639    DO ikl=     1,knonv
640       isl              = isnoSV(ikl)
641
642       dTSurf            = TsisSV(ikl,isl) -     TSurf0(ikl)
643      TsisSV(ikl,isl)   = TSurf0(ikl) + sign(1.,dTSurf) & ! 180.0 dgC/hr
644            * min(abs(dTSurf),5.e-2*dt__SV)         ! =0.05 dgC/s
645
646
647
648    END DO
649
650    DO ikl=     1,knonv
651    DO   isl=min(nsno,isnoSV(ikl)+1),1      ,-1
652      TsisSV(ikl,isl)   = max(Ts_Min,       TsisSV(ikl,isl))
653      TsisSV(ikl,isl)   = min(Ts_Max,       TsisSV(ikl,isl))
654    END DO
655
656    END DO
657
658  ! +--Update Surface    Fluxes
659  ! +  ========================
660
661
662
663    DO ikl=      1,knonv
664      isl         = isnoSV(ikl)
665      IRs_SV(ikl) = IRs__D(ikl) & !
666            - dIRsdT(ikl) * TsisSV(ikl,isl)        !
667      HSs_sv(ikl) = HS___D(ikl) & ! Sensible Heat
668            - dSdTSV(ikl) * TsisSV(ikl,isl)        ! Downward > 0
669      HLs_sv(ikl) = HL___D(ikl) & ! Latent   Heat
670            - dLdTSV(ikl) * TsisSV(ikl,isl)        ! Downward > 0
671    END DO
672
673
674END SUBROUTINE sisvat_tso
Note: See TracBrowser for help on using the repository browser.