source: LMDZ6/branches/Amaury_dev/libf/phylmd/inlandsis/sisvat_ts2.f90 @ 5143

Last change on this file since 5143 was 5143, checked in by abarral, 3 months ago

Put YOEGWD.h, FCTTRE.h into modules
File size: 18.7 KB
Line 
1SUBROUTINE SISVAT_TS2
2  ! #ES.                     (ETSo_0,ETSo_1,ETSo_d)
3
4  ! +------------------------------------------------------------------------+
5  ! | MAR          SISVAT_TS2                            Mon 16-08-2009  MAR |
6  ! |   SubRoutine SISVAT_TS2 computes the Soil/Snow temperature and fluxes  |
7  ! |   using the same method as in LMDZ for consistency.                    |
8  ! |   The corresponding LMDZ routines are soil (soil.F90) and calcul_fluxs |
9  ! |   (calcul_fluxs_mod.F90).                                              |
10  ! +------------------------------------------------------------------------+
11  ! |                                                                        |
12  ! |                                                                        |
13  ! |   PARAMETERS:  klonv: Total Number of columns =                        |
14  ! |   ^^^^^^^^^^        = Total Number of grid boxes of surface type       |
15  ! |                       (land ice for now)                               |
16  ! |                                                                        |
17  ! |   INPUT:   isnoSV   = total Nb of Ice/Snow Layers                      |
18  ! |   ^^^^^    sol_SV   : Downward Solar Radiation                  [W/m2] |
19  ! |            IRd_SV   : Surface Downward  Longwave   Radiation    [W/m2] |
20  ! |            VV__SV   : SBL Top    Wind Speed                      [m/s] |
21  ! |            TaT_SV   : SBL Top    Temperature                       [K] |
22  ! |            QaT_SV   : SBL Top    Specific  Humidity            [kg/kg] |
23  ! |            dzsnSV   : Snow Layer Thickness                         [m] |
24  ! |            dt__SV   : Time Step                                    [s] |
25  ! |                                                                        |
26  ! |            SoSosv   : Absorbed Solar Radiation by Surfac.(Normaliz)[-] |
27  ! |            Eso_sv   : Soil+Snow       Emissivity                   [-] |
28  ! |   ?        rah_sv   : Aerodynamic Resistance for Heat            [s/m] |
29  ! |                                                                        |
30  ! |            dz1_SV    : "inverse" layer thickness (centered)      [1/m] |
31  ! |            dz2_SV    : layer thickness (layer above (?))           [m] |
32  ! |            AcoHSV    : coefficient for Enthalpy evolution, from atm.   |
33  ! |            AcoHSV    : coefficient for Enthalpy evolution, from atm.   |
34  ! |            AcoQSV    : coefficient for Humidity evolution, from atm.   |
35  ! |            BcoQSV    : coefficient for Humidity evolution, from atm.   |
36  ! |            ps__SV    : surface pressure                           [Pa] |
37  ! |            p1l_SV    : 1st atmospheric layer pressure             [Pa] |
38  ! |            cdH_SV    : drag coeff Energy (?)                           |
39  ! |            rsolSV    : Radiation balance surface                [W/m2] |
40  ! |            lambSV    : Coefficient for soil layer geometry         [-] |
41  ! |                                                                        |
42  ! |   INPUT /  TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
43  ! |   OUTPUT:           & Snow     Temperatures (layers  1,2,...,nsno) [K] |
44  ! |   ^^^^^^   rsolSV   : Radiation balance surface                 [W/m2] |
45  ! |                                                                        |
46  ! |   OUTPUT:  IRs_SV   : Soil      IR Radiation                    [W/m2] |
47  ! |   ^^^^^^   HSs_sv   : Sensible  Heat Flux                       [W/m2] |
48  ! |            HLs_sv   : Latent    Heat Flux                       [W/m2] |
49  ! |            TsfnSV   : new surface temperature                      [K] |
50  ! |            Evp_sv   : Evaporation                              [kg/m2] |
51  ! |            dSdTSV   : Sensible Heat Flux temp. derivative     [W/m2/K] |
52  ! |            dLdTSV   : Latent Heat Flux temp. derivative       [W/m2/K] |
53  ! |                                                                        |
54  ! |   ?        ETSo_0   : Snow/Soil Energy Power, before Forcing    [W/m2] |
55  ! |   ?        ETSo_1   : Snow/Soil Energy Power, after  Forcing    [W/m2] |
56  ! |   ?        ETSo_d   : Snow/Soil Energy Power         Forcing    [W/m2] |
57  ! |                                                                        |
58  ! |________________________________________________________________________|
59
60  USE VAR_SV
61  USE VARdSV
62
63  USE VARySV
64  USE VARtSV
65  USE VARxSV
66  USE VARphy
67  USE indice_sol_mod
68  USE lmdz_comsoil, ONLY: inertie_sol, inertie_sno, inertie_sic, inertie_lic, iflag_sic, iflag_inertie
69  USE lmdz_YOETHF
70  USE lmdz_fcttre, ONLY: foeew, foede, qsats, qsatl, dqsats, dqsatl, thermcep
71
72  IMPLICIT NONE
73
74
75  ! +--Global Variables
76  ! +  ================
77
78  INCLUDE "YOMCST.h"
79  ! INCLUDE "indicesol.h"
80  ! include  "LMDZphy.inc"
81
82  ! +--OUTPUT for Stand Alone NetCDF File
83  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
84  ! #NC      real*8        SOsoKL(klonv)             ! Absorbed Solar Radiation
85  ! #NC      real*8        IRsoKL(klonv)             ! Absorbed IR    Radiation
86  ! #NC      real*8        HSsoKL(klonv)             ! Absorbed Sensible Heat Flux
87  ! #NC      real*8        HLsoKL(klonv)             ! Absorbed Latent   Heat Flux
88  ! #NC      real*8        HLs_KL(klonv)             ! Evaporation
89  ! #NC      real*8        HLv_KL(klonv)             ! Transpiration
90  ! #NC      common/DumpNC/SOsoKL,IRsoKL
91  ! #NC     .             ,HSsoKL,HLsoKL
92  ! #NC     .             ,HLs_KL,HLv_KL
93
94  ! +--Internal Variables
95  ! +  ==================
96
97  INTEGER :: ig, jk, isl
98  REAL :: mu
99  REAL :: Tsrf(klonv)               ! surface temperature as extrapolated from soil
100  REAL :: mug(klonv)                 !hj coef top layers
101  REAL :: ztherm_i(klonv), zdz2(klonv, -nsol:nsno), z1s
102  REAL :: pfluxgrd(klonv), pcapcal(klonv), cal(klonv)
103  REAL :: beta(klonv), dif_grnd(klonv)
104  REAL :: C_coef(klonv, -nsol:nsno), D_coef(klonv, -nsol:nsno)
105
106  REAL, DIMENSION(klonv) :: zx_mh, zx_nh, zx_oh
107  REAL, DIMENSION(klonv) :: zx_mq, zx_nq, zx_oq
108  REAL, DIMENSION(klonv) :: zx_pkh, zx_dq_s_dt, zx_qsat, zx_coef
109  REAL, DIMENSION(klonv) :: zx_sl, zx_k1
110  REAL, DIMENSION(klonv) :: d_ts
111  REAL :: zdelta, zcvm5, zx_qs, zcor, zx_dq_s_dh
112  REAL :: qsat_new, q1_new
113  ! REAL, PARAMETER          :: t_grnd = 271.35, t_coup = 273.15
114  ! REAL, PARAMETER          :: max_eau_sol = 150.0
115  REAL, DIMENSION(klonv) :: IRs__D, dIRsdT
116
117  REAL :: t_grnd                      ! not used
118  parameter(t_grnd = 271.35)       !
119  REAL :: t_coup                      ! distinguish evap/sublimation
120  parameter(t_coup = 273.15)       !
121  REAL :: max_eau_sol
122  parameter(max_eau_sol = 150.0)
123
124
125  ! WRITE(*,*)'T check'
126
127  ! DO  ig = 1,knonv
128  !     DO  jk = 1,isnoSV(ig) !nsno
129  !       IF (TsisSV(ig,jk) <= 1.) THEN          !hj check
130  !         TsisSV(ig,jk) = TsisSV(ig,isnoSV(ig))
131  !       ENDIF
132
133  !       IF (TsisSV(ig,jk) <= 1.) THEN          !hj check
134  !         TsisSV(ig,jk) = 273.15
135  !       ENDIF
136  !     END DO
137  !        END DO
138
139  !!=======================================================================
140  !! I. First part: corresponds to soil.F90 in LMDZ
141  !!=======================================================================
142
143  DO ig = 1, knonv
144    DO jk = 1, isnoSV(ig)
145      dz2_SV(ig, jk) = dzsnSV(ig, jk)
146      !! use arithmetic center between layers to derive dz1 for snow layers for simplicity:
147      dz1_SV(ig, jk) = 2. / (dzsnSV(ig, jk) + dzsnSV(ig, jk - 1))
148    ENDDO
149  ENDDO
150
151  DO ig = 1, knonv
152    ztherm_i(ig) = inertie_lic
153    IF (isnoSV(ig) > 0) ztherm_i(ig) = inertie_sno
154  ENDDO
155
156  !!-----------------------------------------------------------------------
157  !! 1)
158  !! Calculation of Cgrf and Dgrd coefficients using soil temperature from
159  !! previous time step.
160  !!
161  !! These variables are recalculated on the local compressed grid instead
162  !! of saved in restart file.
163  !!-----------------------------------------------------------------------
164  DO ig = 1, knonv
165    DO jk = -nsol, nsno
166      zdz2(ig, jk) = dz2_SV(ig, jk) / dt__SV                                !ptimestep
167    ENDDO
168  ENDDO
169
170  DO ig = 1, knonv
171    z1s = zdz2(ig, -nsol) + dz1_SV(ig, -nsol + 1)
172    C_coef(ig, -nsol + 1) = zdz2(ig, -nsol) * TsisSV(ig, -nsol) / z1s
173    D_coef(ig, -nsol + 1) = dz1_SV(ig, -nsol + 1) / z1s
174  ENDDO
175
176  DO ig = 1, knonv
177    DO jk = -nsol + 1, isnoSV(ig) - 1, 1
178      z1s = 1. / (zdz2(ig, jk) + dz1_SV(ig, jk + 1) + dz1_SV(ig, jk)           &
179              * (1. - D_coef(ig, jk)))
180      C_coef(ig, jk + 1) = &
181              (TsisSV(ig, jk) * zdz2(ig, jk)                              &
182                      + dz1_SV(ig, jk) * C_coef(ig, jk)) * z1s
183      D_coef(ig, jk + 1) = dz1_SV(ig, jk + 1) * z1s
184    ENDDO
185  ENDDO
186
187  !!-----------------------------------------------------------------------
188  !! 2)
189  !! Computation of the soil temperatures using the Cgrd and Dgrd
190  !! coefficient computed above
191  !!
192  !!-----------------------------------------------------------------------
193  !! Extrapolate surface Temperature                   !hj check
194  mu = 1. / ((2.**1.5 - 1.) / (2.**(0.5) - 1.) - 1.)
195
196  ! IF (knonv>0) THEN
197  !  DO ig=1,8
198  !    WRITE(*,*)ig,'sisvat: Tsis ',TsisSV(ig,isnoSV(ig))
199  !    WRITE(*,*)'max-1            ',TsisSV(ig,isnoSV(ig)-1)
200  !    WRITE(*,*)'max-2            ',TsisSV(ig,isnoSV(ig)-2)
201  !    WRITE(*,*)'0                ',TsisSV(ig,0)
202  !!        WRITE(*,*)min(max(isnoSV(ig),0),1),max(1-isnoSV(ig),0)
203  !  ENDDO
204  ! END IF
205
206  DO ig = 1, knonv
207    IF (isnoSV(ig)>0) THEN
208      IF (isnoSV(ig)>1) THEN
209        mug(ig) = 1. / (1. + dzsnSV(ig, isnoSV(ig) - 1) / dzsnSV(ig, isnoSV(ig))) !mu
210      ELSE
211        mug(ig) = 1. / (1. + dzsnSV(ig, isnoSV(ig) - 1) / dz_dSV(0)) !mu
212      ENDIF
213    ELSE
214      mug(ig) = lambSV
215    ENDIF
216
217    IF (mug(ig)  <= 0.05) THEN
218      WRITE(*, *)'Attention mu low', mug(ig)
219    ENDIF
220    IF (mug(ig)  >= 0.98) THEN
221      WRITE(*, *)'Attention mu high', mug(ig)
222    ENDIF
223
224    Tsrf(ig) = (1.5 * TsisSV(ig, isnoSV(ig)) - 0.5 * TsisSV(ig, isnoSV(ig) - 1))&
225            * min(max(isnoSV(ig), 0), 1) + &
226            ((mug(ig) + 1) * TsisSV(ig, 0) - mug(ig) * TsisSV(ig, -1))       &
227                    * max(1 - isnoSV(ig), 0)
228  ENDDO
229
230
231
232  !! Surface temperature
233  DO ig = 1, knonv
234    TsisSV(ig, isnoSV(ig)) = (mug(ig) * C_coef(ig, isnoSV(ig)) + Tsf_SV(ig)) / &
235            (mug(ig) * (1. - D_coef(ig, isnoSV(ig))) + 1.)
236  ENDDO
237
238  !! Other temperatures
239  DO ig = 1, knonv
240    DO jk = isnoSV(ig), -nsol + 1, -1
241      TsisSV(ig, jk - 1) = C_coef(ig, jk) + D_coef(ig, jk)                   &
242              * TsisSV(ig, jk)
243    ENDDO
244  ENDDO
245  ! WRITE(*,*)ig,'Tsis',TsisSV(ig,0)
246
247  ! IF (indice == is_sic) THEN
248  !   DO ig = 1,knonv
249  !     TsisSV(ig,-nsol) = RTT - 1.8
250  !   END DO
251  ! ENDIF
252
253  !C      !hj new 11 03 2010
254  DO ig = 1, knonv
255    isl = isnoSV(ig)
256    ! dIRsdT(ig) = Eso_sv(ig)* SteBo   * 4.                        & ! - d(IR)/d(T)
257    ! &                             * Tsf_SV(ig)                         & !T TsisSV(ig,isl)           !
258    ! &                             * Tsf_SV(ig)                         & !TsisSV(ig,isl)           !
259    ! &                             * Tsf_SV(ig) !TsisSV(ig,isl)           !
260    !      IRs__D(ig) = dIRsdT(ig)* Tsf_SV(ig) * 0.75 !TsisSV(ig,isl) * 0.75    !:
261    dIRsdT(ig) = Eso_sv(ig) * StefBo * 4.                        & ! - d(IR)/d(T)
262            * TsisSV(ig, isl)                     & !
263            * TsisSV(ig, isl)                     & !
264            * TsisSV(ig, isl)
265    IRs__D(ig) = dIRsdT(ig) * TsisSV(ig, isl) * 0.75                  !:
266  END DO
267  !hj
268  !!-----------------------------------------------------------------------
269  !! 3)
270  !! Calculate the Cgrd and Dgrd coefficient corresponding to actual soil
271  !! temperature
272  !!-----------------------------------------------------------------------
273  DO ig = 1, knonv
274    z1s = zdz2(ig, -nsol) + dz1_SV(ig, -nsol + 1)
275    C_coef(ig, -nsol + 1) = zdz2(ig, -nsol) * TsisSV(ig, -nsol) / z1s
276    D_coef(ig, -nsol + 1) = dz1_SV(ig, -nsol + 1) / z1s
277  ENDDO
278
279  DO ig = 1, knonv
280    DO jk = -nsol + 1, isnoSV(ig) - 1, 1
281      z1s = 1. / (zdz2(ig, jk) + dz1_SV(ig, jk + 1) + dz1_SV(ig, jk)           &
282              * (1. - D_coef(ig, jk)))
283      C_coef(ig, jk + 1) = (TsisSV(ig, jk) * zdz2(ig, jk) + &
284              dz1_SV(ig, jk) * C_coef(ig, jk)) * z1s
285      D_coef(ig, jk + 1) = dz1_SV(ig, jk + 1) * z1s
286    ENDDO
287  ENDDO
288
289  !!-----------------------------------------------------------------------
290  !! 4)
291  !! Computation of the surface diffusive flux from ground and
292  !! calorific capacity of the ground
293  !!-----------------------------------------------------------------------
294  DO ig = 1, knonv
295    !! (pfluxgrd)
296    pfluxgrd(ig) = ztherm_i(ig) * dz1_SV(ig, isnoSV(ig)) * &
297            (C_coef(ig, isnoSV(ig)) + (D_coef(ig, isnoSV(ig)) - 1.)           &
298                    * TsisSV(ig, isnoSV(ig)))
299    !! (pcapcal)
300    pcapcal(ig) = ztherm_i(ig) * &
301            (dz2_SV(ig, isnoSV(ig)) + dt__SV * (1. - D_coef(ig, isnoSV(ig)))    &
302                    * dz1_SV(ig, isnoSV(ig)))
303    z1s = mug(ig) * (1. - D_coef(ig, isnoSV(ig))) + 1.
304    pcapcal(ig) = pcapcal(ig) / z1s
305    pfluxgrd(ig) = (pfluxgrd(ig)                                   &
306            + pcapcal(ig) * (TsisSV(ig, isnoSV(ig)) * z1s               &
307                    - mug(ig) * C_coef(ig, isnoSV(ig))                           &
308                    - Tsf_SV(ig)) / dt__SV)
309  ENDDO
310
311  cal(1:knonv) = RCPD / pcapcal(1:knonv)
312  rsolSV(1:knonv) = rsolSV(1:knonv) + pfluxgrd(1:knonv)
313  !!=======================================================================
314  !! II. Second part: corresponds to calcul_fluxs_mod.F90 in LMDZ
315  !!=======================================================================
316
317  Evp_sv = 0.
318  ! #NC HSsoKL=0.
319  ! #NC HLsoKL=0.
320  dSdTSV = 0.
321  dLdTSV = 0.
322
323  beta(:) = 1.0
324  dif_grnd(:) = 0.0
325
326  !! zx_qs = qsat en kg/kg
327  !!**********************************************************************x***************
328
329  DO ig = 1, knonv
330    IF (ps__SV(ig)<1.) THEN
331      ! WRITE(*,*)'ig',ig,'ps',ps__SV(ig)
332      ps__SV(ig) = max(ps__SV(ig), 1.e-8)
333    ENDIF
334    IF (p1l_SV(ig)<1.) THEN
335      ! WRITE(*,*)'ig',ig,'p1l',p1l_SV(ig)
336      p1l_SV(ig) = max(p1l_SV(ig), 1.e-8)
337    ENDIF
338    IF (TaT_SV(ig)<180.) THEN
339      ! WRITE(*,*)'ig',ig,'TaT',TaT_SV(ig)
340      TaT_SV(ig) = max(TaT_SV(ig), 180.)
341    ENDIF
342    IF (QaT_SV(ig)<1.e-8) THEN
343      ! WRITE(*,*)'ig',ig,'QaT',QaT_SV(ig)
344      QaT_SV(ig) = max(QaT_SV(ig), 1.e-8)
345    ENDIF
346    IF (Tsf_SV(ig)<100.) THEN
347      ! WRITE(*,*)'ig',ig,'Tsf',Tsf_SV(ig)
348      Tsf_SV(ig) = max(Tsf_SV(ig), 180.)
349    ENDIF
350    IF (Tsf_SV(ig)>500.) THEN
351      ! WRITE(*,*)'ig',ig,'Tsf',Tsf_SV(ig)
352      Tsf_SV(ig) = min(Tsf_SV(ig), 400.)
353    ENDIF
354    ! IF (Tsrf(ig).LT.1.) THEN
355    !!          WRITE(*,*)'ig',ig,'Tsrf',Tsrf(ig)
356    !   Tsrf(ig)=max(Tsrf(ig),TaT_SV(ig)-20.)
357    ! ENDIF
358    IF (cdH_SV(ig)<1.e-10) THEN
359      ! IF (ig.le.3)   WRITE(*,*)'ig',ig,'cdH',cdH_SV(ig)
360      cdH_SV(ig) = .5
361    ENDIF
362  ENDDO
363
364  DO ig = 1, knonv
365    zx_pkh(ig) = 1. ! (ps__SV(ig)/ps__SV(ig))**RKAPPA
366    IF (thermcep) THEN
367      zdelta = MAX(0., SIGN(1., rtt - Tsf_SV(ig)))
368      zcvm5 = R5LES * LhvH2O * (1. - zdelta) + R5IES * LhsH2O * zdelta
369      zcvm5 = zcvm5 / RCPD / (1.0 + RVTMP2 * QaT_SV(ig))
370      zx_qs = r2es * FOEEW(Tsf_SV(ig), zdelta) / ps__SV(ig)
371      zx_qs = MIN(0.5, zx_qs)
372      !WRITE(*,*)'zcor',retv*zx_qs
373      zcor = 1. / (1. - retv * zx_qs)
374      zx_qs = zx_qs * zcor
375      zx_dq_s_dh = FOEDE(Tsf_SV(ig), zdelta, zcvm5, zx_qs, zcor)        &
376              / LhvH2O / zx_pkh(ig)
377    ELSE
378      IF (Tsf_SV(ig)<t_coup) THEN
379        zx_qs = qsats(Tsf_SV(ig)) / ps__SV(ig)
380        zx_dq_s_dh = dqsats(Tsf_SV(ig), zx_qs) / LhvH2O               &
381                / zx_pkh(ig)
382      ELSE
383        zx_qs = qsatl(Tsf_SV(ig)) / ps__SV(ig)
384        zx_dq_s_dh = dqsatl(Tsf_SV(ig), zx_qs) / LhvH2O               &
385                / zx_pkh(ig)
386      ENDIF
387    ENDIF
388    zx_dq_s_dt(ig) = RCPD * zx_pkh(ig) * zx_dq_s_dh
389    zx_qsat(ig) = zx_qs
390    ! zx_coef(ig) = cdH_SV(ig) *                                     &
391    ! &       (1.0+SQRT(u1lay(ig)**2+v1lay(ig)**2)) *                   &
392    ! &       p1l_SV(ig)/(RD*t1lay(ig))
393    zx_coef(ig) = cdH_SV(ig) * &
394            (1.0 + VV__SV(ig)) * &
395            p1l_SV(ig) / (RD * TaT_SV(ig))
396
397  ENDDO
398
399
400  !! === Calcul de la temperature de surface ===
401  !! zx_sl = chaleur latente d'evaporation ou de sublimation
402  !!****************************************************************************
403
404  DO ig = 1, knonv
405    zx_sl(ig) = LhvH2O
406    IF (Tsf_SV(ig) < RTT) zx_sl(ig) = LhsH2O
407    zx_k1(ig) = zx_coef(ig)
408  ENDDO
409
410  DO ig = 1, knonv
411    !! Q
412    zx_oq(ig) = 1. - (beta(ig) * zx_k1(ig) * BcoQSV(ig) * dt__SV)
413    zx_mq(ig) = beta(ig) * zx_k1(ig) * &
414            (AcoQSV(ig) - zx_qsat(ig) + &
415                    zx_dq_s_dt(ig) * Tsf_SV(ig))                               &
416            / zx_oq(ig)
417    zx_nq(ig) = beta(ig) * zx_k1(ig) * (-1. * zx_dq_s_dt(ig))       &
418            / zx_oq(ig)
419
420    !! H
421    zx_oh(ig) = 1. - (zx_k1(ig) * BcoHSV(ig) * dt__SV)
422    zx_mh(ig) = zx_k1(ig) * AcoHSV(ig) / zx_oh(ig)
423    zx_nh(ig) = - (zx_k1(ig) * RCPD * zx_pkh(ig)) / zx_oh(ig)
424
425    !! surface temperature
426    TsfnSV(ig) = (Tsf_SV(ig) + cal(ig) / RCPD * zx_pkh(ig) * dt__SV * &
427            (rsolSV(ig) + zx_mh(ig) + zx_sl(ig) * zx_mq(ig))           &
428            + dif_grnd(ig) * t_grnd * dt__SV) / &
429            (1. - dt__SV * cal(ig) / (RCPD * zx_pkh(ig)) * &
430                    (zx_nh(ig) + zx_sl(ig) * zx_nq(ig))                        &
431                    + dt__SV * dif_grnd(ig))
432
433    !hj rajoute 22 11 2010 tuning...
434    TsfnSV(ig) = min(RTT + 0.02, TsfnSV(ig))
435
436    d_ts(ig) = TsfnSV(ig) - Tsf_SV(ig)
437
438
439    !!== flux_q est le flux de vapeur d'eau: kg/(m**2 s)  positive vers bas
440    !!== flux_t est le flux de cpt (energie sensible): j/(m**2 s)
441    Evp_sv(ig) = - zx_mq(ig) - zx_nq(ig) * TsfnSV(ig)
442    HLs_sv(ig) = - Evp_sv(ig) * zx_sl(ig)
443    HSs_sv(ig) = zx_mh(ig) + zx_nh(ig) * TsfnSV(ig)
444
445    !! Derives des flux dF/dTs (W m-2 K-1):
446    dSdTSV(ig) = zx_nh(ig)
447    dLdTSV(ig) = zx_sl(ig) * zx_nq(ig)
448
449
450    !hj  new 11 03 2010
451    isl = isnoSV(ig)
452    ! TsisSV(ig,isl) = TsfnSV(ig)
453    IRs_SV(ig) = IRs__D(ig)                                         & !
454            - dIRsdT(ig) * TsfnSV(ig) !TsisSV(ig,isl)?      !
455
456    ! hj
457    ! #NC   SOsoKL(ig) = sol_SV(ig) * SoSosv(ig)              ! Absorbed Sol.
458    ! #NC   IRsoKL(ig) =               IRs_SV(ig)                           & !Up Surf. IR
459    ! #NC&        +     tau_sv(ig)      *IRd_SV(ig)*Eso_sv(ig)              & !Down Atm IR
460    ! #NC&        -(1.0-tau_sv(ig)) *0.5*IRv_sv(ig)            ! Down Veg IR
461    ! #NC   HLsoKL(ig) = HLs_sv(ig)
462    ! #NC   HSsoKL(ig) = HSs_sv(ig)
463    ! #NC   HLs_KL(ig) = Evp_sv(ig)
464
465    !! Nouvelle valeure de l'humidite au dessus du sol
466    qsat_new = zx_qsat(ig) + zx_dq_s_dt(ig) * d_ts(ig)
467    q1_new = AcoQSV(ig) - BcoQSV(ig) * Evp_sv(ig) * dt__SV
468    QaT_SV(ig) = q1_new * (1. - beta(ig)) + beta(ig) * qsat_new
469
470  ENDDO
471
472END SUBROUTINE sisvat_ts2
Note: See TracBrowser for help on using the repository browser.