GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/wx_pbl_mod.F90 Lines: 0 348 0.0 %
Date: 2023-06-30 12:51:15 Branches: 0 176 0.0 %

Line Branch Exec Source
1
MODULE wx_pbl_mod
2
!
3
! Split Planetary Boundary Layer
4
!
5
! This module manages the splitting of the boundary layer between two regions; the (w)
6
! region (inside cold pools) and the (x) region (outside cold pools)
7
!
8
  USE dimphy
9
10
  IMPLICIT NONE
11
12
CONTAINS
13
!
14
!****************************************************************************************
15
!
16
SUBROUTINE wx_pbl0_merge(knon, ypplay, ypaprs,  &
17
                                 sigw, dTs_forcing, dqs_forcing,  &
18
                                 yt_x, yt_w, yq_x, yq_w, &
19
                                 yu_x, yu_w, yv_x, yv_w, &
20
                                 ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w, &
21
                                 ycdragm_x, ycdragm_w, &
22
                                 AcoefT_x, AcoefT_w, AcoefQ_x, AcoefQ_w, &
23
                                 AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w, &
24
                                 BcoefT_x, BcoefT_w, BcoefQ_x, BcoefQ_w, &
25
                                 BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w, &
26
                                 AcoefT, AcoefQ, AcoefU, AcoefV, &
27
                                 BcoefT, BcoefQ, BcoefU, BcoefV, &
28
                                 ycdragh, ycdragq, ycdragm, &
29
                                 yt1, yq1, yu1, yv1 &
30
                                 )
31
!
32
33
    USE wx_pbl_var_mod
34
35
    USE print_control_mod, ONLY: prt_level,lunout
36
    USE indice_sol_mod, ONLY: is_oce
37
!
38
    INCLUDE "YOMCST.h"
39
    INCLUDE "FCTTRE.h"
40
    INCLUDE "YOETHF.h"
41
    INCLUDE "clesphys.h"
42
!
43
    INTEGER,                      INTENT(IN)        :: knon    ! number of grid cells
44
    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypplay  ! mid-layer pressure (Pa)
45
    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypaprs  ! pressure at layer interfaces (pa)
46
    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw ! cold pools fractional area
47
    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs_forcing ! forced temperature difference (w)-(x)
48
    REAL, DIMENSION(knon),        INTENT(IN)        :: dqs_forcing ! forced humidity difference (w)-(x)
49
    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: yt_x, yt_w, yq_x, yq_w
50
    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: yu_x, yu_w, yv_x, yv_w
51
    REAL, DIMENSION(knon),        INTENT(IN)        :: ycdragh_x, ycdragh_w, ycdragq_x, ycdragq_w
52
    REAL, DIMENSION(knon),        INTENT(IN)        :: ycdragm_x, ycdragm_w
53
    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT_x, AcoefT_w, AcoefQ_x, AcoefQ_w
54
    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefU_x, AcoefU_w, AcoefV_x, AcoefV_w
55
    REAL, DIMENSION(knon),        INTENT(IN)        :: BcoefT_x, BcoefT_w, BcoefQ_x, BcoefQ_w
56
    REAL, DIMENSION(knon),        INTENT(IN)        :: BcoefU_x, BcoefU_w, BcoefV_x, BcoefV_w
57
    REAL, DIMENSION(knon),        INTENT(OUT)       :: AcoefT, AcoefQ, AcoefU, AcoefV
58
    REAL, DIMENSION(knon),        INTENT(OUT)       :: BcoefT, BcoefQ, BcoefU, BcoefV
59
    REAL, DIMENSION(knon),        INTENT(OUT)       :: ycdragh, ycdragq, ycdragm
60
    REAL, DIMENSION(knon),        INTENT(OUT)       :: yt1, yq1, yu1, yv1  ! Apparent T, q, u, v at first level, as
61
                                                                           !seen by surface modules
62
!
63
! Local variables
64
    INTEGER                    :: j
65
    REAL                       :: dd_Kh
66
    REAL                       :: dd_Kq
67
    REAL                       :: dd_Km
68
    REAL                       :: dd_u
69
    REAL                       :: dd_v
70
    REAL                       :: dd_t
71
    REAL                       :: dd_q
72
!
73
    REAL                       :: LambdaTs, LambdaQs, LambdaUs, LambdaVs
74
!
75
    REAL, DIMENSION(knon)      :: sigx       ! fractional area of (x) region
76
!
77
!
78
   sigx(1:knon) = 1.-sigw(1:knon)
79
!
80
!
81
        DO j=1,knon
82
!
83
!
84
! Compute w-x differences
85
       dd_t = yt_w(j,1) - yt_x(j,1)
86
       dd_q = yq_w(j,1) - yq_x(j,1)
87
       dd_u = yu_w(j,1) - yu_x(j,1)
88
       dd_v = yv_w(j,1) - yv_x(j,1)
89
!
90
! Merged exchange coefficients
91
         dd_Kh = Kech_h_w(j) - Kech_h_x(j)
92
         dd_Kq = Kech_q_w(j) - Kech_q_x(j)
93
         dd_Km = Kech_m_w(j) - Kech_m_x(j)
94
!
95
       LambdaTs = dd_KTp(j)/Kech_Tp(j)
96
       LambdaQs = dd_KQs(j)/Kech_Qs(j)
97
       LambdaUs = dd_KUp(j)/Kech_Up(j)
98
       LambdaVs = dd_KVp(j)/Kech_Vp(j)
99
!
100
! Calcul des coef A, B \'equivalents dans la couche 1
101
!
102
! The dTs_forcing and dqs_forcing terms are added for diagnostic purpose ; they should be zero in normal operation.
103
       AcoefT(j) = AcoefT_x(j) + sigw(j)*(1.+sigx(j)*LambdaTs)*(dd_AT(j) - C_p(j)*dTs_forcing(j))
104
       AcoefQ(j) = AcoefQ_x(j) + sigw(j)*(1.+sigx(j)*LambdaQs)*(dd_AQ(j) - dqs_forcing(j))
105
       AcoefU(j) = AcoefU_x(j) + sigw(j)*(1.+sigx(j)*LambdaUs)*dd_AU(j)
106
       AcoefV(j) = AcoefV_x(j) + sigw(j)*(1.+sigx(j)*LambdaVs)*dd_AV(j)
107
!
108
!
109
!!       BcoefT(j) = (sigw(j)*Kech_h_w(j)*Kech_T_pw(j)*BcoefT_w(j) + &
110
!!                sigx(j)*Kech_h_x(j)*Kech_T_px(j)*BcoefT_x(j) )/(Kech_h(j)*Kech_Tp(j))
111
!!       BcoefQ(j) = (sigw(j)*Kech_q_w(j)*Kech_Q_pw(j)*BcoefQ_w(j) + &
112
!!                sigx(j)*Kech_q_x(j)*Kech_Q_px(j)*BcoefQ_x(j) )/(Kech_q(j)*Kech_Qp(j))
113
!!       BcoefU(j) = (sigw(j)*Kech_m_w(j)*Kech_U_pw(j)*BcoefU_w(j) + &
114
!!                sigx(j)*Kech_m_x(j)*Kech_U_px(j)*BcoefU_x(j) )/(Kech_m(j)*Kech_Up(j))
115
!!       BcoefV(j) = (sigw(j)*Kech_m_w(j)*Kech_V_pw(j)*BcoefV_w(j) + &
116
!!                sigx(j)*Kech_m_x(j)*Kech_V_px(j)*BcoefV_x(j) )/(Kech_m(j)*Kech_Vp(j))
117
!
118
!!  Print *,'YYYYpbl0: BcoefT_x, sigw, sigx, dd_Kh, dd_KTp, Kech_h_w ', &
119
!!                     BcoefT_x, sigw, sigx, dd_Kh, dd_KTp, Kech_h_w
120
!!  Print *,'YYYYpbl0: Kech_T_pw, dd_BT, Kech_h, Kech_Tp ', &
121
!!                     Kech_T_pw, dd_BT, Kech_h, Kech_Tp
122
       BcoefT(j) = BcoefT_x(j) + sigw(j)*(sigx(j)*dd_Kh*dd_KTp(j)*BcoefT_x(j) + &
123
                                  Kech_h_w(j)*Kech_T_pw(j)*dd_BT(j))/(Kech_h(j)*Kech_Tp(j))
124
       BcoefQ(j) = BcoefQ_x(j) + sigw(j)*(sigx(j)*dd_Kq*dd_KQs(j)*BcoefQ_x(j) + &
125
                                  Kech_q_w(j)*Kech_Q_sw(j)*dd_BQ(j))/(Kech_q(j)*Kech_Qs(j))
126
       BcoefU(j) = BcoefU_x(j) + sigw(j)*(sigx(j)*dd_Km*dd_KUp(j)*BcoefU_x(j) + &
127
                                  Kech_m_w(j)*Kech_U_pw(j)*dd_BU(j))/(Kech_m(j)*Kech_Up(j))
128
       BcoefV(j) = BcoefV_x(j) + sigw(j)*(sigx(j)*dd_Km*dd_KVp(j)*BcoefV_x(j) + &
129
                                  Kech_m_w(j)*Kech_V_pw(j)*dd_BV(j))/(Kech_m(j)*Kech_Vp(j))
130
!>jyg
131
!
132
!
133
! Calcul des cdrag \'equivalents dans la couche
134
!
135
       ycdragm(j) = ycdragm_x(j) + sigw(j)*dd_Cdragm(j)
136
       ycdragh(j) = ycdragh_x(j) + sigw(j)*dd_Cdragh(j)
137
       ycdragq(j) = ycdragq_x(j) + sigw(j)*dd_Cdragq(j)
138
!
139
! Calcul de T, q, u et v \'equivalents dans la couche 1
140
!!       yt1(j) = yt_x(j,1) + sigw(j)*dd_t*(1.+sigx(j)*dd_Kh/KCT)
141
!!       yq1(j) = yq_x(j,1) + sigw(j)*dd_q*(1.+sigx(j)*dd_Kh/KCQ)
142
!!       yu1(j) = yu_x(j,1) + sigw(j)*dd_u*(1.+sigx(j)*dd_Km/KCU)
143
!!       yv1(j) = yv_x(j,1) + sigw(j)*dd_v*(1.+sigx(j)*dd_Km/KCV)
144
       yt1(j) = yt_x(j,1) + sigw(j)*dd_t
145
       yq1(j) = yq_x(j,1) + sigw(j)*dd_q
146
       yu1(j) = yu_x(j,1) + sigw(j)*dd_u
147
       yv1(j) = yv_x(j,1) + sigw(j)*dd_v
148
149
150
        ENDDO
151
152
        RETURN
153
154
END SUBROUTINE wx_pbl0_merge
155
156
SUBROUTINE wx_pbl_dts_merge(knon, dtime, ypplay, ypaprs, &
157
                                 sigw, beta, wcstar, wdens, &
158
                                 AT_x, AT_w, &
159
                                 BT_x, BT_w, &
160
                                 AcoefT0, AcoefQ0, BcoefT0, BcoefQ0, &
161
                                 AcoefT,  AcoefQ,  BcoefT,  BcoefQ, &
162
                                 HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
163
                                 phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0, &
164
                                 g_T, g_Q, &
165
                                 Gamma_phiT, Gamma_phiQ, &
166
                                 dTs_ins, dqsatsrf_ins &
167
                                 )
168
!
169
170
    USE wx_pbl_var_mod
171
172
    USE print_control_mod, ONLY: prt_level,lunout
173
!
174
    INCLUDE "YOMCST.h"
175
    INCLUDE "FCTTRE.h"
176
    INCLUDE "YOETHF.h"
177
!
178
    INTEGER,                      INTENT(IN)        :: knon    ! number of grid cells
179
    REAL,                         INTENT(IN)        :: dtime   ! time step size (s)
180
    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypplay  ! mid-layer pressure (Pa)
181
    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypaprs  ! pressure at layer interfaces (pa)
182
    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw    ! cold pool fractional area
183
    REAL, DIMENSION(knon),        INTENT(IN)        :: beta    ! evaporation by potential evaporation
184
    REAL, DIMENSION(knon),        INTENT(IN)        :: wcstar   ! cold pool gust front speed
185
    REAL, DIMENSION(knon),        INTENT(IN)        :: wdens    ! cold pool number density
186
    REAL, DIMENSION(knon),        INTENT(IN)        :: AT_x, AT_w
187
    REAL, DIMENSION(knon),        INTENT(IN)        :: BT_x, BT_w
188
    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT0, AcoefQ0, BcoefT0, BcoefQ0
189
!
190
    REAL, DIMENSION(knon),        INTENT(OUT)       :: AcoefT, AcoefQ, BcoefT, BcoefQ
191
    REAL, DIMENSION(knon),        INTENT(OUT)       :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn
192
    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiT0_b, dphiT0, phiQ0_b, dphiQ0, Rn0_b, dRn0
193
    REAL, DIMENSION(knon),        INTENT(OUT)       :: g_T, g_Q
194
    REAL, DIMENSION(knon),        INTENT(OUT)       :: Gamma_phiT, Gamma_phiQ
195
    REAL, DIMENSION(knon),        INTENT(OUT)       :: dTs_ins, dqsatsrf_ins
196
!
197
! Local variables
198
    REAL, DIMENSION(knon)      :: qsat_x
199
    REAL, DIMENSION(knon)      :: qsat_w
200
    REAL, DIMENSION(knon)      :: dqsatdT_x
201
    REAL, DIMENSION(knon)      :: dqsatdT_w
202
!
203
    REAL, DIMENSION(knon)      :: T10_x
204
    REAL, DIMENSION(knon)      :: T10_w
205
    REAL, DIMENSION(knon)      :: phiT0_x
206
    REAL, DIMENSION(knon)      :: phiT0_w
207
    REAL, DIMENSION(knon)      :: phiQ0_x
208
    REAL, DIMENSION(knon)      :: phiQ0_w
209
    REAL, DIMENSION(knon)      :: Rn0_x
210
    REAL, DIMENSION(knon)      :: Rn0_w
211
    REAL, DIMENSION(knon)      :: Rp1_x
212
    REAL, DIMENSION(knon)      :: Rp1_w
213
    REAL, DIMENSION(knon)      :: Rps_x
214
    REAL, DIMENSION(knon)      :: Rps_w
215
!
216
    REAL, DIMENSION(knon)      :: HTphiT_x
217
    REAL, DIMENSION(knon)      :: HTphiT_w
218
    REAL, DIMENSION(knon)      :: HTphiQ_x
219
    REAL, DIMENSION(knon)      :: HTphiQ_w
220
    REAL, DIMENSION(knon)      :: HTRn_x
221
    REAL, DIMENSION(knon)      :: HTRn_w
222
!
223
    REAL, DIMENSION(knon)      :: HQphiT_x
224
    REAL, DIMENSION(knon)      :: HQphiT_w
225
    REAL, DIMENSION(knon)      :: HQphiQ_x
226
    REAL, DIMENSION(knon)      :: HQphiQ_w
227
    REAL, DIMENSION(knon)      :: HQRn_x
228
    REAL, DIMENSION(knon)      :: HQRn_w
229
!
230
    REAL, DIMENSION(knon)      :: HQphiT_b
231
    REAL, DIMENSION(knon)      :: dd_HQphiT
232
    REAL, DIMENSION(knon)      :: HQphiQ_b
233
    REAL, DIMENSION(knon)      :: dd_HQphiQ
234
    REAL, DIMENSION(knon)      :: HQRn_b
235
    REAL, DIMENSION(knon)      :: dd_HQRn
236
!
237
238
    REAL, DIMENSION(knon)    :: sigx
239
!
240
    REAL, DIMENSION(knon)    :: Ts, T1
241
!!!    REAL, DIMENSION(knon)    :: qsat, dqsat_dT
242
!!!    REAL, DIMENSION(knon)    :: phiT0
243
!
244
!!!    REAL, DIMENSION(knon)    :: Cp, Lv
245
    REAL, DIMENSION(knon)    :: tau, Inert
246
!
247
    REAL                     :: dd_Kh
248
    REAL                     :: zdelta, zcvm5, zcor
249
    REAL                     :: qsat
250
!
251
    INTEGER                  :: j
252
253
254
!----------------------------------------------------------------------------
255
!  Reference state
256
!  ---------------
257
!   dqsat_dT_w = dqsat_dT(Ts0_w)                          dqsat_dT_x = dqsat_dT(Ts0_x)
258
!   T10_w = (AT_w/Cp - Kech_T_w BT_w dtime Ts0_w)/(1 - Kech_T_w BT_w dtime)
259
!                                                T10_x = (AT_x/Cp - Kech_T_x BT_x dtime Ts0_x)/(1 - Kech_T_x BT_x dtime)
260
!   phiT0_w = Kech_T_pw (AT_w - Cp Ts0_w)                 phiT0_x = Kech_T_px (AT_x - Cp Ts0_x)
261
!   phiQ0_w = Kech_Q_sw (beta AQ_w - qsatsrf0_w)          phiQ0_x = Kech_Q_sx (beta AQ_x - qsatsrf0_x)
262
!   Rn0_w = eps_1 Rsigma T10_w^4 - Rsigma Ts0_w^4         Rn0_x = eps_1 Rsigma T10_x^4 - Rsigma Ts0_x^4
263
!   Rp1_w = 4 eps_1 Rsigma T10_w^3                        Rp1_x = 4 eps_1 Rsigma T10_x^3
264
!   Rps_w = 4 Rsigma Ts0_w^3                              Rps_x = 4 Rsigma Ts0_x^3
265
!
266
!   phiT0_b = sigw phiT0_w + sigx phiT0_x
267
!   dphiT0 = phiT0_w - phiT0_x
268
!   phiQ0_b = sigw phiQ0_w + sigx phiQ0_x
269
!   dphiQ0 = phiQ0_w - phiQ0_x
270
!   Rn0_b = sigw Rn0_w + sigx Rn0_x
271
    dRn0 = Rn0_w - Rn0_x
272
!
273
!
274
!----------------------------------------------------------------------------
275
!  Elementary enthalpy equations
276
!  -----------------------------
277
!   phiT_w = phiT0_w - HTphiT_w (Ts_w-Ts0_w)            phiT_x = phiT0_x - HTphiT_x (Ts_x-Ts0_x)
278
!   phiQ_w = phiQ0_w - HTphiQ_w (Ts_w-Ts0_w)            phiQ_x = phiQ0_x - HTphiQ_x (Ts_x-Ts0_x)
279
!   Rn_w   = Rn0_w   - HTRn_w   (Ts_w-Ts0_w)            Rn_x   = Rn0_x   - HTRn_x   (Ts_x-Ts0_x)
280
!  DFlux_DT coefficients
281
!  ---------------------
282
!   Heat flux equation
283
!     HTphiT_w = Cp Kech_T_pw                            HTphiT_x = Cp Kech_T_px
284
!   Moisture flux equation
285
!     HTphiQ_w = beta Kech_Q_sw dqsat_dT_w               HTphiQ_x = beta Kech_Q_sx dqsat_dT_x
286
!   Radiation equation
287
!     HTRn_w = Rp1_w Kech_T_pw BcoefT_w dtime + Rps_w    HTRn_x = Rp1_x Kech_T_px BcoefT_x dtime + Rps_x
288
!
289
!----------------------------------------------------------------------------
290
!  Elementary moisture equations
291
!  -----------------------------
292
!   beta Ts_w   = beta Ts0_w    + QQ_w     (qsatsrf_w-qsatsrf0_w)    beta Ts_x   = beta Ts0_x    + QQ_x     (qsatsrf_x-qsatsrf0_x)
293
!   beta phiT_w = beta phiT0_w - HQphiT_w (qsatsrf_w-qsatsrf0_w)    beta phiQ_x = beta phiQ0_x - HTphiQ_x (qsatsrf_x-qsatsrf0_x)
294
!   beta phiQ_w = beta phiQ0_w - HQphiQ_w (qsatsrf_w-qsatsrf0_w)    beta phiQ_x = beta phiQ0_x - HTphiQ_x (qsatsrf_x-qsatsrf0_x)
295
!   beta Rn_w   = beta Rn0_w   - HQRn_w   (qsatsrf_w-qsatsrf0_w)    beta Rn_x   = beta Rn0_x   - HTRn_x   (qsatsrf_x-qsatsrf0_x)
296
!  DFluxDQ coefficients
297
!  ---------------------
298
!   dqsat_dT equation
299
!     QQ_w = 1. / dqsat_dT_w                             QQ_x = 1. / dqsat_dT_x
300
!   Heat flux equation
301
!     HQphiT_w = Cp Kech_T_pw QQ_w                       HQphiT_x = Cp Kech_T_px QQ_x
302
!   Moisture flux equation
303
!     HQphiQ_w = beta Kech_Q_sw                          HQphiQ_x = beta Kech_Q_sx
304
!   Radiation equation
305
!     HQRn_w = (Rp1_w Kech_T_pw BcoefT_w dtime + Rps_w) QQ_w
306
!                                         HQRn_x = (Rp1_x Kech_T_px BcoefT_x dtime + Rps_x) QQ_x
307
!
308
!----------------------------------------------------------------------------
309
! Mean values and w-x differences
310
! -------------------------------
311
!  HTphiT_b = sigw HTphiT_w + sigx HTphiT_x               dd_HTphiT = HTphiT_w - HTphiT_x
312
!  HTphiQ_b = sigw HTphiQ_w + sigx HTphiQ_x               dd_HTphiQ = HTphiQ_w - HTphiQ_x
313
!  HTRn_b   = sigw HTRn_w   + sigx HTRn_x                 dd_HTRn   = HTRn_w   - HTRn_x
314
!
315
!  QQ_b     = sigw QQ_w     + sigx QQ_x                   dd_QQ     = QQ_w     - QQ_x
316
!  HQphiT_b = sigw HQphiT_w + sigx HQphiT_x               dd_HQphiT = HQphiT_w - HQphiT_x
317
!  HQphiQ_b = sigw HQphiQ_w + sigx HQphiQ_x               dd_HQphiQ = HQphiQ_w - HQphiQ_x
318
!  HQRn_b   = sigw HQRn_w   + sigx HQRn_x                 dd_HQRn   = HQRn_w   - HQRn_x
319
!
320
!----------------------------------------------------------------------------
321
!  Equations
322
!  ---------
323
! (1 - g_T) dTs    = dTs_ins    + Gamma_phiT phiT
324
! (1 - g_Q) dqsatsrf = dqsatsrf_ins + Gamma_phiQ phiQ
325
!
326
! Feedback Gains
327
! --------------
328
! g_T = - (sqrt(tau)/I) [ HTphiT_b + Lv HTphiQ_b + HTRn_b +  &
329
!                        (dd_HTphiT + Lv dd_HTphiQ + dd_HTRn) (sigx - sigw - sigw sigx dd_HTphiT/HTphiT_b) ]
330
! g_Q = - (sqrt(tau)/(I QQ_b)) ( HQphiT_b + Lv HQphiQ_b + HQRn_b ) -  &
331
!         (sigx - sigw - sigw sigx dd_HQphiQ/HQphiQ_b)   &
332
!                          [ dd_QQ/QQ_b + (sqrt(tau)/(I QQ_b))(dd_HQphiT + Lv dd_HQphiQ + dd_HQRn) ]
333
!
334
!  Ts, qs Coupling coefficients                /
335
!  ----------------------------
336
! Gamma_phiT = (sqrt(tau)/(I HTphiT_b)) (dd_HTphiT + Lv dd_HTphiQ + dd_HTRn)
337
! Gamma_phiQ = (1/(HQphiQ_b QQ_b)) [ dd_QQ +  (sqrt(tau)/(I )) (dd_HQphiT + Lv dd_HQphiQ + dd_HQRn) ]
338
!
339
!  Insensitive changes
340
!  -------------------
341
! dTs_ins    = (1 - g_T) dTs0    - Gamma_phiT phiT0_b
342
! dqsatsrf_ins = (1 - g_Q) dqsatsrf0 - Gamma_phiQ phiQ0_b
343
!
344
!----------------------------------------------------------------------------
345
!  Effective coefficients Acoef and Bcoef
346
!  --------------------------------------
347
!  Equations
348
!  ---------
349
! Cp Ta = AcoefT + BcoefT phiT dtime
350
!    qa = AcoefQ + BcoefQ phiQ dtime
351
!  Coefficients
352
!  ------------
353
! AcoefT = AcoefT0 - sigw sigx (dd_KTp/Kech_Tp) Cp dTs_ins/(1 - g_T)
354
! BcoefT = BcoefT0 - sigw sigx (dd_KTp/Kech_Tp) Cp Gamma_phiT/(1 - g_T)/dtime
355
!
356
! AcoefQ = AcoefQ0 - sigw sigx (dd_KQp/Kech_Qp) dqs_ins/(1 - g_Q)
357
! BcoefQ = BcoefQ0 - sigw sigx (dd_KQp/Kech_Qp) Gamma_phiq/(1 - g_Q)/dtime
358
!
359
!==============================================================================
360
!
361
!
362
!  Parameters
363
!  ----------
364
   Inert(1:knon) = 2000.
365
   tau(1:knon) = sqrt(sigw(1:knon)/max(rpi*wdens(1:knon)*wcstar(1:knon)**2 , &
366
                                       sigw(1:knon)*1.e-12,smallestreal))
367
   sigx(1:knon) = 1.-sigw(1:knon)
368
!! Compute Cp, Lv, qsat, dqsat_dT.
369
!   C_p(1:knon) = RCpd
370
!   L_v(1:knon) = RLvtt
371
!
372
!      print *,' AAAA wx_pbl_dTs, C_p(j), qsat0(j), Ts0(j) : ', C_p(:), qsat0(:), Ts0(:)
373
!
374
!
375
   T10_x(1:knon) = (AT_x(1:knon)/C_p(1:knon) - Kech_h_x(1:knon)*BT_x(1:knon)*dtime*Ts0_x(1:knon))/  &
376
                   (1 - Kech_h_x(1:knon)*BT_x(1:knon)*dtime)
377
   T10_w(1:knon) = (AT_w(1:knon)/C_p(1:knon) - Kech_h_w(1:knon)*BT_w(1:knon)*dtime*Ts0_w(1:knon))/  &
378
                   (1 - Kech_h_w(1:knon)*BT_w(1:knon)*dtime)
379
!
380
   phiT0_x(1:knon) = Kech_T_px(1:knon)*(AT_x(1:knon) - C_p(1:knon)*Ts0_x(1:knon))
381
   phiT0_w(1:knon) = Kech_T_pw(1:knon)*(AT_w(1:knon) - C_p(1:knon)*Ts0_w(1:knon))
382
!
383
   phiQ0_x(1:knon) = Kech_Q_sx(1:knon)*(beta(1:knon)*AQ_x(1:knon) - qsatsrf0_x(1:knon))
384
   phiQ0_w(1:knon) = Kech_Q_sw(1:knon)*(beta(1:knon)*AQ_w(1:knon) - qsatsrf0_w(1:knon))
385
!
386
   Rn0_x(1:knon) = eps_1*Rsigma*T10_x(1:knon)**4 - Rsigma*Ts0_x(1:knon)**4
387
   Rn0_w(1:knon) = eps_1*Rsigma*T10_w(1:knon)**4 - Rsigma*Ts0_w(1:knon)**4
388
!
389
   Rp1_x(1:knon) = 4*eps_1*Rsigma*T10_x(1:knon)**3
390
   Rp1_w(1:knon) = 4*eps_1*Rsigma*T10_w(1:knon)**3
391
!
392
   Rps_x(1:knon) = 4*Rsigma*Ts0_x(1:knon)**3
393
   Rps_w(1:knon) = 4*Rsigma*Ts0_w(1:knon)**3
394
!
395
!  DFlux_DT coefficients
396
!  ---------------------
397
!   Heat flux equation
398
     HTphiT_x(1:knon) = C_p(1:knon)*Kech_T_px(1:knon)
399
     HTphiT_w(1:knon) = C_p(1:knon)*Kech_T_pw(1:knon)
400
!   Moisture flux equation
401
     HTphiQ_x(1:knon) = beta(1:knon)*Kech_Q_sx(1:knon)*dqsatdT0_x(1:knon)
402
     HTphiQ_w(1:knon) = beta(1:knon)*Kech_Q_sw(1:knon)*dqsatdT0_w(1:knon)
403
!   Radiation equation
404
     HTRn_x(1:knon) = Rp1_x(1:knon)*Kech_T_px(1:knon)*BT_x(1:knon)*dtime + Rps_x(1:knon)
405
     HTRn_w(1:knon) = Rp1_w(1:knon)*Kech_T_pw(1:knon)*BT_w(1:knon)*dtime + Rps_w(1:knon)
406
!
407
!  DFluxDQ coefficients
408
!  ---------------------
409
!   Heat flux equation
410
     HQphiT_x(1:knon) = C_p(1:knon)*Kech_T_px(1:knon)*QQ_x(1:knon)
411
     HQphiT_w(1:knon) = C_p(1:knon)*Kech_T_pw(1:knon)*QQ_w(1:knon)
412
!   Moisture flux equation
413
     HQphiQ_x(1:knon) = beta(1:knon)*Kech_Q_sx(1:knon)
414
     HQphiQ_w(1:knon) = beta(1:knon)*Kech_Q_sw(1:knon)
415
!   Radiation equation
416
     HQRn_x(1:knon) = (Rp1_x(1:knon)*Kech_T_px(1:knon)*BT_x(1:knon)*dtime + Rps_x(1:knon))*QQ_x(1:knon)
417
     HQRn_w(1:knon) = (Rp1_w(1:knon)*Kech_T_pw(1:knon)*BT_w(1:knon)*dtime + Rps_w(1:knon))*QQ_w(1:knon)
418
!
419
! Mean values and w-x differences
420
! -------------------------------
421
  phiT0_b(1:knon) = sigw(1:knon)*phiT0_w(1:knon) + sigx(1:knon)*phiT0_x(1:knon)
422
  phiQ0_b(1:knon) = sigw(1:knon)*phiQ0_w(1:knon) + sigx(1:knon)*phiQ0_x(1:knon)
423
  Rn0_b(1:knon)   = sigw(1:knon)*Rn0_w(1:knon)   + sigx(1:knon)*Rn0_x(1:knon)
424
!
425
  dphiT0(1:knon) = phiT0_w(1:knon) - phiT0_x(1:knon)
426
  dphiQ0(1:knon) = phiQ0_w(1:knon) - phiQ0_x(1:knon)
427
  dRn0(1:knon)   = Rn0_w(1:knon)   - Rn0_x(1:knon)
428
!
429
  HTphiT_b(1:knon) = sigw(1:knon)*HTphiT_w(1:knon) + sigx(1:knon)*HTphiT_x(1:knon)
430
  dd_HTphiT(1:knon) = HTphiT_w(1:knon) - HTphiT_x(1:knon)
431
!
432
  HTphiQ_b(1:knon) = sigw(1:knon)*HTphiQ_w(1:knon) + sigx(1:knon)*HTphiQ_x(1:knon)
433
  dd_HTphiQ(1:knon) = HTphiQ_w(1:knon) - HTphiQ_x(1:knon)
434
!
435
  HTRn_b(1:knon)   = sigw(1:knon)*HTRn_w(1:knon)   + sigx(1:knon)*HTRn_x(1:knon)
436
  dd_HTRn(1:knon)   = HTRn_w(1:knon)   - HTRn_x(1:knon)
437
!
438
  HQphiT_b(1:knon) = sigw(1:knon)*HQphiT_w(1:knon) + sigx(1:knon)*HQphiT_x(1:knon)
439
  dd_HQphiT(1:knon) = HQphiT_w(1:knon) - HQphiT_x(1:knon)
440
!
441
  HQphiQ_b(1:knon) = sigw(1:knon)*HQphiQ_w(1:knon) + sigx(1:knon)*HQphiQ_x(1:knon)
442
  dd_HQphiQ(1:knon) = HQphiQ_w - HQphiQ_x(1:knon)
443
!
444
  HQRn_b(1:knon)   = sigw(1:knon)*HQRn_w(1:knon)   + sigx(1:knon)*HQRn_x(1:knon)
445
  dd_HQRn(1:knon)   = HQRn_w(1:knon)   - HQRn_x(1:knon)
446
!
447
! Feedback Gains
448
! --------------
449
 g_T(1:knon) = - (sqrt(tau(1:knon))/Inert(1:knon))  &
450
               * (HTphiT_b(1:knon) + L_v(1:knon)*HTphiQ_b(1:knon) + HTRn_b(1:knon)  &
451
                 + (dd_HTphiT(1:knon) + L_v(1:knon)*dd_HTphiQ(1:knon) + dd_HTRn(1:knon))  &
452
                 * (sigx(1:knon) - sigw(1:knon) - sigw(1:knon)*sigx(1:knon)*dd_HTphiT(1:knon)/HTphiT_b(1:knon)) )
453
!
454
!!!! DO j = 1,knon
455
!!!!  IF (mod(j,20) .eq.0) THEN
456
!!!!   print *, '   j     dd_QQ       QQ_b  dd_HQphiQ  dd_HQphiT   dd_HQRn   HQphiQ_b   HQphiT_b     HQRn_b '
457
!!!!  ENDIF
458
!!!!   print 1789, j, dd_QQ(j), QQ_b(j), dd_HQphiQ(j), dd_HQphiT(j), dd_HQRn(j), HQphiQ_b(j), HQphiT_b(j), HQRn_b(j)
459
!!!! 1789 FORMAT( I4, 10(1X,E10.2))
460
!!!! ENDDO
461
   g_Q(1:knon) = - (dd_QQ(1:knon)/QQ_b(1:knon)) *  &
462
                    (sigx(1:knon)-sigw(1:knon)-sigw(1:knon)*sigx(1:knon)*dd_KQs(1:knon)/Kech_Qs(1:knon)) &
463
                 - sqrt(tau(1:knon))/(Inert(1:knon)*QQ_b(1:knon)) *  &
464
                   ( HQphiT_b(1:knon) + L_v(1:knon)*HQphiQ_b(1:knon) + HQRn_b(1:knon) +  &
465
                      (sigx(1:knon) - sigw(1:knon) - sigw(1:knon)*sigx(1:knon)*dd_KQs(1:knon)/Kech_Qs(1:knon)) *  &
466
                       (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
467
468
!!   g_Q(1:knon) = - (dd_QQ(1:knon)/QQ_b(1:knon)) *  &
469
!!                    (sigx(1:knon)-sigw(1:knon)-sigw(1:knon)*sigx(1:knon)*dd_HQphiQ(1:knon)/HQphiQ_b(1:knon)) &
470
!!                 - sqrt(tau(1:knon))/(Inert(1:knon)*QQ_b(1:knon)) *  &
471
!!                   ( HQphiT_b(1:knon) + L_v(1:knon)*HQphiQ_b(1:knon) + HQRn_b(1:knon) +  &
472
!!                      (sigx(1:knon) - sigw(1:knon) - sigw(1:knon)*sigx(1:knon)*dd_HQphiQ(1:knon)/HQphiQ_b(1:knon)) *  &
473
!!                       (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
474
475
!! g_Q(1:knon) = - (sqrt(tau(1:knon))/(Inert(1:knon)*QQ_b(1:knon))) *  &
476
!!                 ( HQphiT_b(1:knon) + L_v(1:knon)*HQphiQ_b(1:knon) + HQRn_b(1:knon) )  &
477
!!               - (sigx(1:knon) - sigw(1:knon) - sigw(1:knon)*sigx(1:knon)*dd_HQphiQ(1:knon)/HQphiQ_b(1:knon)) *   &
478
!!                 ( dd_QQ(1:knon)/QQ_b(1:knon)   &
479
!!                 + (sqrt(tau(1:knon))/(Inert(1:knon)*QQ_b(1:knon)))  &
480
!!                 * (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
481
482
!  Ts, qs Coupling coefficients                /
483
!  ----------------------------
484
  Gamma_phiT(1:knon) = (sqrt(tau(1:knon))/(Inert(1:knon)*HTphiT_b(1:knon)))  &
485
                     * (dd_HTphiT(1:knon) + L_v(1:knon)*dd_HTphiQ(1:knon) + dd_HTRn(1:knon))
486
!
487
  Gamma_phiQ(1:knon) = (1./(Kech_Qs(1:knon)*QQ_b(1:knon))) * &
488
                        ( dd_QQ(1:knon)   &
489
                         + (sqrt(tau(1:knon))/(Inert(1:knon))) *  &
490
                          (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
491
492
!!  Gamma_phiQ(1:knon) = (beta(1:knon)/(HQphiQ_b(1:knon)*QQ_b(1:knon))) * &
493
!!                        ( dd_QQ(1:knon)   &
494
!!                         + (sqrt(tau(1:knon))/(Inert(1:knon))) *  &
495
!!                          (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
496
497
!!  Gamma_phiQ(1:knon) = (1/(HQphiQ_b(1:knon)*QQ_b(1:knon)))   &
498
!!                     * ( dd_QQ(1:knon)   &
499
!!                       + (sqrt(tau(1:knon))/(Inert(1:knon)))  &
500
!!                       * (dd_HQphiT(1:knon) + L_v(1:knon)*dd_HQphiQ(1:knon) + dd_HQRn(1:knon)) )
501
!
502
!  Insensitive changes
503
!  -------------------
504
  dTs_ins(1:knon)    = (sqrt(tau(1:knon))/Inert(1:knon))*  &
505
                       (dphiT0(1:knon) + L_v(1:knon)*dphiQ0(1:knon) + dRn0(1:knon))
506
!
507
  dqsatsrf_ins(1:knon) = (beta(1:knon)/QQ_b(1:knon))*dTs_ins(1:knon)
508
!
509
   IF (prt_level .Ge. 10) THEN
510
      print *,'wx_pbl_merge, tau         ', tau
511
      print *,'wx_pbl_merge, AcoefT0     ', AcoefT0
512
      print *,'wx_pbl_merge, AcoefQ0     ', AcoefQ0
513
      print *,'wx_pbl_merge, BcoefT0     ', BcoefT0
514
      print *,'wx_pbl_merge, BcoefQ0     ', BcoefQ0
515
      print *,'wx_pbl_merge, qsat0_w, qsat0_x ', (qsat0_w(j), qsat0_x(j),j=1,knon)
516
      print *,'wx_pbl_merge, dqsatdT0_w, dqsatdT0_x ', (dqsatdT0_w(j), dqsatdT0_x(j),j=1,knon)
517
   ENDIF
518
!
519
!----------------------------------------------------------------------------
520
!
521
!------------------------------------------------------------------------------
522
!
523
!    Effective coefficients Acoef and Bcoef
524
!    --------------------------------------
525
   DO j = 1,knon
526
     AcoefT(j) = AcoefT0(j) - sigw(j)*sigx(j)*(dd_KTp(j)/Kech_Tp(j))*C_p(j)*   &
527
                 (dTs0(j) + (dTs_ins(j)-dTs0(j)-Gamma_phiT(j)*phiT0_b(j))/(1. - g_T(j)))
528
     BcoefT(j) = BcoefT0(j) - sigw(j)*sigx(j)*(dd_KTp(j)/Kech_Tp(j))*C_p(j)*Gamma_phiT(j)/(1. - g_T(j))/dtime
529
530
     AcoefQ(j) = AcoefQ0(j) - sigw(j)*sigx(j)*(dd_KQs(j)/Kech_Qs(j))*    &
531
                 (dqsatsrf0(j) + (dqsatsrf_ins(j)-(beta(j)/QQ_b(j))*dTs0(j)-Gamma_phiQ(j)*phiQ0_b(j))/(1 - g_Q(j)))/ &
532
                 max(beta(j),1.e-4)
533
     BcoefQ(j) = BcoefQ0(j) - sigw(j)*sigx(j)*(dd_KQs(j)/Kech_Qs(j))*Gamma_phiQ(j)/(1 - g_Q(j))/ &
534
                 (max(beta(j),1.e-4)*dtime)
535
!!     AcoefQ(j) = AcoefQ0(j) - sigw(j)*sigx(j)*(dd_KQs(j)/Kech_Qs(j))*    &
536
!!                 (dqsatsrf0(j) + (dqsatsrf_ins(j)-(beta(j)/QQ_b(j))*dTs0(j)-Gamma_phiQ(j)*phiQ0_b(j))/(1 - g_Q(j)))/ &
537
!!                 beta(j)
538
!!     BcoefQ(j) = BcoefQ0(j) - sigw(j)*sigx(j)*(dd_KQs(j)/Kech_Qs(j))*Gamma_phiQ(j)/(1 - g_Q(j))/(beta(j)*dtime)
539
   ENDDO ! j = 1,knon
540
541
   IF (prt_level .Ge. 10) THEN
542
   print *,'wx_pbl_dts AAAA BcoefQ, BcoefQ0, sigw ', &
543
                            BcoefQ, BcoefQ0, sigw
544
      print *,'wx_pbl_dts_merge, dTs_ins      ', dTs_ins
545
      print *,'wx_pbl_dts_merge, dqs_ins      ', dqsatsrf_ins
546
   ENDIF
547
548
     RETURN
549
550
END SUBROUTINE wx_pbl_dts_merge
551
552
SUBROUTINE wx_pbl_split(knon, nsrf, dtime, sigw, beta, iflag_split, &
553
                       g_T, g_Q, &
554
                       Gamma_phiT, Gamma_phiQ, &
555
                       dTs_ins, dqsatsrf_ins, &
556
                       phiT, phiQ, phiU, phiV, &
557
!!!!                       HTRn_b, dd_HTRn, HTphiT_b, dd_HTphiT, &
558
                       phiQ0_b, phiT0_b, &
559
                       phiT_x, phiT_w, &
560
                       phiQ_x, phiQ_w, &
561
                       phiU_x, phiU_w, &
562
                       phiV_x, phiV_w, &
563
                       philat_x, philat_w, &
564
!!!!                       Rn_b, dRn, &
565
                       dqsatsrf, &
566
                       dTs, delta_qsurf &
567
                       )
568
!
569
570
    USE wx_pbl_var_mod
571
572
    USE print_control_mod, ONLY: prt_level,lunout
573
    USE indice_sol_mod, ONLY: is_oce
574
!
575
    INCLUDE "YOMCST.h"
576
!
577
    INTEGER,                      INTENT(IN)        :: knon    ! number of grid cells
578
    INTEGER,                      INTENT(IN)        :: nsrf    ! surface type
579
    REAL,                         INTENT(IN)        :: dtime   ! time step size (s)
580
    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw ! cold pools fractional area
581
    REAL, DIMENSION(knon),        INTENT(IN)        :: beta ! aridity factor
582
    INTEGER,                      INTENT(IN)        :: iflag_split
583
    REAL, DIMENSION(knon),        INTENT(IN)        :: g_T, g_Q
584
    REAL, DIMENSION(knon),        INTENT(IN)        :: Gamma_phiT, Gamma_phiQ
585
    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs_ins, dqsatsrf_ins
586
    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT, phiQ, phiU, phiV
587
    REAL, DIMENSION(knon),        INTENT(IN)        :: phiQ0_b, phiT0_b
588
!
589
    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiT_x, phiT_w
590
    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiQ_x, phiQ_w
591
    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiU_x, phiU_w
592
    REAL, DIMENSION(knon),        INTENT(OUT)       :: phiV_x, phiV_w
593
    REAL, DIMENSION(knon),        INTENT(OUT)       :: philat_x, philat_w
594
    REAL, DIMENSION(knon),        INTENT(OUT)       :: dqsatsrf      ! beta delta(qsat(Ts))
595
    REAL, DIMENSION(knon),        INTENT(OUT)       :: dTs           ! Temperature difference at surface
596
    REAL, DIMENSION(knon),        INTENT(OUT)       :: delta_qsurf
597
!
598
!! Local variables
599
    INTEGER                    :: j
600
    REAL, DIMENSION(knon)      :: dphiT, dphiQ, dphiU, dphiV
601
    REAL, DIMENSION(knon)      :: q1_x, q1_w
602
!
603
    REAL, DIMENSION(knon)      :: sigx       ! fractional area of (x) region
604
605
!----------------------------------------------------------------------------
606
!  Equations
607
!  ---------
608
!!!!!! (1 - g_T) dTs    = dTs_ins    + Gamma_phiT phiT
609
!!!!!! (1 - g_Q) dqsatsrf = dqsatsrf_ins + Gamma_phiQ phiQ
610
!!!!!! dphiT = (dd_KTp/KTp) phiT + (     dd_AT - C_p dTs)*KxKwTp/KTp
611
!!!!!! dphiQ = (dd_KQs/KQs) phiQ + (beta dd_AQ - dqsatsrf )*KxKwQs/KQs
612
!!!!!! dphiU = (dd_KUp/KUp) phiU + (     dd_AU          )*KxKwUp/KUp
613
!!!!!! dphiV = (dd_KVp/KVp) phiV + (     dd_AV          )*KxKwVp/KVp
614
!
615
! (1 - g_T) (dTs-dTs0)    = dTs_ins-dTs0    + Gamma_phiT (phiT-phiT0)
616
! (1 - g_Q) dqsatsrf = dqsatsrf_ins + Gamma_phiQ phiQ
617
! dphiT = (dd_KTp/KTp) phiT + (     dd_AT - C_p dTs)*KxKwTp/KTp
618
! dphiQ = (dd_KQs/KQs) phiQ + (beta dd_AQ - dqsatsrf )*KxKwQs/KQs
619
! dphiU = (dd_KUp/KUp) phiU + (     dd_AU          )*KxKwUp/KUp
620
! dphiV = (dd_KVp/KVp) phiV + (     dd_AV          )*KxKwVp/KVp
621
!
622
!!
623
        sigx(:) = 1.-sigw(:)
624
!
625
!      print *,' AAAA wx_pbl_split, C_p(j), qsat0(j), Ts0(j) : ', C_p(:), qsat0(:), Ts0(:)
626
!
627
   IF (iflag_split .EQ. 2 .AND. nsrf .NE. is_oce) THEN
628
!
629
!   Delta_tsurf and  Delta_qsurf computation
630
!   -----------------------------------------
631
      IF (prt_level >=10 ) THEN
632
        print *,' wx_pbl_split, dTs_ins, dTs0 , Gamma_phiT, g_T ', dTs_ins, dTs0, Gamma_phiT, g_T
633
        print *,' wx_pbl_split, dqsatsrf_ins, Gamma_phiQ, g_q ', dqsatsrf_ins, Gamma_phiQ, g_q
634
      ENDIF
635
!
636
      DO j = 1,knon
637
        dTs(j)    = dTs0(j) + (dTs_ins(j) - dTs0(j) + Gamma_phiT(j)*(phiT(j)-phiT0_b(j)) )/(1 - g_T(j))
638
        dqsatsrf(j) = dqsatsrf0(j) + (dqsatsrf_ins(j) - (beta(j)/QQ_b(j))*dTs0(j) + &
639
                       Gamma_phiQ(j)*(phiQ(j)-phiQ0_b(j)) )/(1 - g_Q(j))
640
      ENDDO ! j = 1,knon
641
!
642
        IF (prt_level >=10 ) THEN
643
          print *,' wx_pbl_split, dqsatsrf0, QQ_b ', dqsatsrf0, QQ_b
644
          print *,' wx_pbl_split, phiT0_b, phiT, dTs ', phiT0_b, phiT, dTs
645
          print *,' wx_pbl_split, phiQ0_b, phiQ, dqsatsrf ', phiQ0_b, phiQ, dqsatsrf
646
        ENDIF
647
   ELSE
648
        dTs(:) = 0.
649
        dqsatsrf(:) = 0.
650
   ENDIF ! (iflag_split .EQ. 2 .AND. nsrf .NE. is_oce)
651
!
652
     DO j = 1,knon
653
       dphiT(j) = (phiT(j)*dd_KTp(j) + (        dd_AT(j) - C_p(j)*dTs(j))*KxKwTp(j))/Kech_Tp(j)
654
       dphiQ(j) = (phiQ(j)*dd_KQs(j) + (beta(j)*dd_AQ(j) -        dqsatsrf(j))*KxKwQs(j))/Kech_Qs(j)
655
       dphiU(j) = (phiU(j)*dd_KUp(j) +          dd_AU(j)                 *KxKwUp(j))/Kech_Up(j)
656
       dphiV(j) = (phiV(j)*dd_KVp(j) +          dd_AV(j)                 *KxKwVp(j))/Kech_Vp(j)
657
!
658
       phiT_x(j)=phiT(j) - sigw(j)*dphiT(j)
659
       phiT_w(j)=phiT(j) + sigx(j)*dphiT(j)
660
       phiQ_x(j)=phiQ(j) - sigw(j)*dphiQ(j)
661
       phiQ_w(j)=phiQ(j) + sigx(j)*dphiQ(j)
662
       phiU_x(j)=phiU(j) - sigw(j)*dphiU(j)
663
       phiU_w(j)=phiU(j) + sigx(j)*dphiU(j)
664
       phiV_x(j)=phiV(j) - sigw(j)*dphiV(j)
665
       phiV_w(j)=phiV(j) + sigx(j)*dphiV(j)
666
!
667
       philat_x(j)=phiQ_x(j)*RLVTT
668
       philat_w(j)=phiQ_w(j)*RLVTT
669
     ENDDO ! j = 1,knon
670
!
671
     DO j = 1,knon
672
       q1_x(j) = AQ_x(j) + BQ_x(j)*phiQ_x(j)*dtime
673
       q1_w(j) = AQ_w(j) + BQ_w(j)*phiQ_w(j)*dtime
674
     ENDDO ! j = 1,knon
675
     DO j = 1,knon
676
       delta_qsurf(j) = (1.-beta(j))*(q1_w(j) - q1_x(j)) + dqsatsrf(j)
677
     ENDDO ! j = 1,knon
678
!
679
!!  Do j = 1,knon
680
!!     print *,'XXXsplit : j, q1_x(j), AQ_x(j), BQ_x(j), phiQ_x(j) ', j, q1_x(j), AQ_x(j), BQ_x(j), phiQ_x(j)
681
!!     print *,'XXXsplit : j, q1_w(j), AQ_w(j), BQ_w(j), phiQ_w(j) ', j, q1_w(j), AQ_w(j), BQ_w(j), phiQ_w(j)
682
!!  ENDDO
683
!
684
        IF (prt_level >=10 ) THEN
685
          print *,' wx_pbl_split, phiT, dphiT, dTs ', phiT, dphiT, dTs
686
          print *,' wx_pbl_split, phiQ, dphiQ, dqsatsrf ', phiQ, dphiQ, dqsatsrf
687
        ENDIF
688
!
689
        IF (prt_level >=10 ) THEN
690
!!          print *,' wx_pbl_split, verif dqsatsrf = beta dqsatdT0 dTs '
691
!!          print *,' wx_pbl_split, dqsatsrf, dqsatdT0*dTs ', dqsatsrf, dqsatdT0*dTs
692
        ENDIF
693
!
694
!!    IF (knon .NE. 0) THEN
695
!!       call  iophys_ecrit('sigw', 1,'sigw', '.',sigw)
696
!!       call  iophys_ecrit('phit', 1,'phit', 'W/m2',phit)
697
!!       call  iophys_ecrit('phit_w', 1,'phit_w', 'W/m2',phit_w)
698
!!       call  iophys_ecrit('phit_x', 1,'phit_x', 'W/m2',phit_x)
699
!!       call  iophys_ecrit('phiq', 1,'phiq', 'kg/m2/s',phiq)
700
!!       call  iophys_ecrit('phiq_w', 1,'phiq_w', 'kg/m2/s',phiq_w)
701
!!       call  iophys_ecrit('phiq_x', 1,'phiq_x', 'kg/m2/s',phiq_x)
702
!!       call  iophys_ecrit('q1_w', 1,'q1_w', '.',q1_w)
703
!!       call  iophys_ecrit('q1_x', 1,'q1_x', '.',q1_x)
704
!!    ENDIF  ! (knon .NE. 0)
705
!
706
        RETURN
707
708
END SUBROUTINE wx_pbl_split
709
710
SUBROUTINE wx_pbl_check( knon, dtime, ypplay, ypaprs, &
711
                               sigw, beta, iflag_split,   &
712
                               Ts0_b9, dTs09,   &
713
                               qs_b9, Ts_b9,  &                         ! yqsurf, Tsurf_new
714
                               dTs9, dqsatsrf9,   &
715
                               AcoefT_x, AcoefT_w, &
716
                               BcoefT_x, BcoefT_w, &
717
                               AcoefT0, AcoefQ0, BcoefT0, BcoefQ0, &
718
                               AcoefT,  AcoefQ,  BcoefT,  BcoefQ, &
719
                               phiT_b9, phiQ_b9,  &
720
                               phiT_x9, phiT_w9, &
721
                               phiQ_x9, phiQ_w9 &
722
                               )
723
!
724
725
    USE wx_pbl_var_mod
726
727
    USE print_control_mod, ONLY: prt_level,lunout
728
!
729
    INCLUDE "YOMCST.h"
730
    INCLUDE "FCTTRE.h"
731
    INCLUDE "YOETHF.h"
732
!
733
    INTEGER,                      INTENT(IN)        :: knon         ! number of grid cells
734
    REAL,                         INTENT(IN)        :: dtime        ! time step size (s)
735
    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypplay       ! mid-layer pressure (Pa)
736
    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypaprs       ! pressure at layer interfaces (pa)
737
    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw         ! cold pools fractional area
738
    REAL, DIMENSION(knon),        INTENT(IN)        :: beta         ! aridity factor
739
    INTEGER,                      INTENT(IN)        :: iflag_split
740
    REAL, DIMENSION(knon),        INTENT(IN)        :: Ts0_b9, dTs09
741
    REAL, DIMENSION(knon),        INTENT(IN)        :: qs_b9, Ts_b9         ! yqsurf, Tsurf_new
742
    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs9, dqsatsrf9
743
    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT_x, AcoefT_w
744
    REAL, DIMENSION(knon),        INTENT(IN)        :: BcoefT_x, BcoefT_w
745
    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT0, AcoefQ0, BcoefT0, BcoefQ0
746
!
747
    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT, AcoefQ, BcoefT, BcoefQ
748
    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT_b9, phiQ_b9
749
    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT_x9, phiT_w9
750
    REAL, DIMENSION(knon),        INTENT(IN)        :: phiQ_x9, phiQ_w9
751
!
752
!! Local variables
753
    INTEGER                    :: j
754
    REAL, DIMENSION(knon)      :: sigx                 ! fractional area of (x) region
755
    REAL, DIMENSION(knon)      :: AcoefT_b, AcoefQ_b   ! mean values of AcoefT and AcoefQ
756
    REAL                       :: zzt, zzq, zzqsat
757
    REAL                       :: zdelta, zcvm5, zcor, qsat
758
    REAL, DIMENSION(knon)      :: qsat_w, qsat_x
759
    REAL, DIMENSION(knon)      :: dqsatdT_w, dqsatdT_x
760
    REAL, DIMENSION(knon)      :: qsat_bs              ! qsat(Ts_b)
761
    REAL, DIMENSION(knon)      :: qsat01, dqsatdT01
762
    REAL, DIMENSION(knon)      :: Ts_x, Ts_w, qs_x, qs_w
763
    REAL, DIMENSION(knon)      :: T1_x, T1_w, q1_x, q1_w
764
    REAL, DIMENSION(knon)      :: Rn_x, Rn_w
765
    REAL, DIMENSION(knon)      :: phiQ0_x, phiQ0_w
766
    REAL, DIMENSION(knon)      :: Ta, qa
767
    REAL, DIMENSION(knon)      :: qsatsrf_w, qsatsrf_x, qsatsrf_b
768
    REAL, DIMENSION(knon)      :: qsurf_w, qsurf_x
769
    REAL                       :: dphiT, dphiQ
770
    REAL                       :: dqsatsrf1
771
    REAL                       :: phiT_w1, phiT_w2
772
    REAL                       :: phiT_x1, phiT_x2
773
    REAL                       :: phiQ_w1, phiQ_w2, phiQ_w3
774
    REAL                       :: phiQ_x1, phiQ_x2, phiQ_x3
775
    REAL                       :: phiT_b1, phiQ_b1
776
    REAL                       :: Kech_Q_sw1, Kech_Q_sx1
777
    REAL                       :: evap_pot
778
779
!----------------------------------------------------------------------------
780
! Equations to be checked:
781
! -----------------------
782
!  Input : Ts0_b, dTs0, Ts_b, dTs, qsatsrf_b, dqsatsrf,
783
!          phiT_b, phiQ_b, phiT_w, phiT_x, phiQ_w, phiQ_x,
784
!
785
!          AcoefT, AcoefQ, AcoefT_w, AcoefQ_w, AcoefT_x, AcoefQ_x,
786
!          BcoefT, BcoefQ, BcoefT_w, BcoefQ_w, BcoefT_x, BcoefQ_x
787
!
788
!  C_p T1_w = AcoefT_w + BcoefT_w phiT_w Delta t          C_p T1_x = AcoefT_x + BcoefT_x phiT_x Delta t
789
!  q1_w = AQ_w + BQ_w phiQ_w Delta t                      q1_x = AQ_x + BQ_x phiQ_x Delta t
790
!  qsatsrf_w = beta qsat(Ts_w)                            qsatsrf_x = beta qsat(Ts_x)
791
!  qsurf_w = (1-beta) q1_w + qsatsrf_w                    qsurf_x = (1-beta) q1_x + qsatsrf_x
792
!  phiT_w = Kech_h_w C_p ( T1_w - Ts_w)                   phiT_x = Kech_h_x C_p ( T1_x - Ts_x)
793
!  phiT_w = Kech_T_pw ( AcoefT_w - C_p Ts_w)              phiT_x = Kech_T_px ( AcoefT_x - C_p Ts_x)
794
!  phiq_w = Kech_h_w ( beta q1_w - qsatsrf_w)             phiq_x = Kech_h_x ( beta q1_x - qsatsrf_x))
795
!  phiq_w = Kech_Q_sw (beta AQ_w -qsatsrf_w)              phiq_x = Kech_Q_sx (beta AQ_x -qsatsrf_x)
796
!  phiq_w = Kech_h_w (q1_w - qsurf_w)                     phiq_x = Kech_h_x (q1_x - qsurf_x)
797
!  phiT_b = sigw phiT_w + sigx phiT_x                     dphiT = phiT_w - phiT_x
798
!  phiQ_b = sigw phiQ_w + sigx phiQ_x                     dphiQ = phiQ_w - phiQ_x
799
!  Ts_b = sigw Ts_w + sigx Ts_x                           dTs = Ts_w - Ts_x
800
!  qsatsrf_b = sigw qsatsrf_w + sigx qsatsrf_x
801
!  C_p Ta = AcoefT + BcoefT phiT_b Delta t
802
!  qa = AcoefQ + BcoefQ phiQ_b Delta t
803
!  phiT_b = Kech_h C_p (Ta - Ts_b)
804
!  phiQ_b = beta Kech_h (qa - qsatsrf_b)
805
!  dTs = sqrt(tau)/I (dphit + L_v dphiq + dR)
806
807
!----------------------------------------------------------------------------
808
!
809
!!
810
        sigx(:) = 1.-sigw(:)
811
        AcoefT_b(1:knon) = AcoefT_x(1:knon) + sigw(1:knon)*dd_AT(1:knon)
812
        AcoefQ_b(1:knon) = AQ_x(1:knon) + sigw(1:knon)*dd_AQ(1:knon)
813
814
! Compute the three qsat and dqsatdTs
815
! ---------------------------------------------
816
!!   C_p(1:knon) = RCpd
817
!!   L_v(1:knon) = RLvtt
818
    IF (prt_level >=10 ) THEN
819
      print *,' AAAA wx_pbl_check, C_p(j), qsat0(j), Ts0(j) : ', C_p(:), qsat0(:), Ts0(:)
820
    ENDIF ! (prt_level >=10 )
821
!
822
   DO j = 1, knon
823
      zdelta = MAX(0.,SIGN(1.,RTT-Ts0_b9(j)))
824
      zcvm5 = R5LES*(1.-zdelta) + R5IES*zdelta
825
      qsat = R2ES*FOEEW(Ts0_b9(j),zdelta)/ypaprs(j,1)
826
      qsat = MIN(0.5,qsat)
827
      zcor = 1./(1.-RETV*qsat)
828
      qsat01(j) = fqsat*qsat*zcor
829
!!      dqsatdT0(j) = FOEDE(Ts0_b(j),zdelta,zcvm5,qsat0(j),zcor)/RLVTT    ! jyg 20210116
830
!!      dqsatdT0(j) = (RLvtt*(1.-zdelta)+RLSTT*zdelta)*qsat0(j)/(Rv*Ts0_b(j)*Ts0_b(j))
831
      dqsatdT01(j) = fqsat*FOEDE(Ts0_b9(j),zdelta,zcvm5,qsat01(j),zcor)
832
   ENDDO
833
!
834
!--------------------------------------------------------------------------------------------------
835
        IF (prt_level >=10 ) THEN
836
!
837
          DO j = 1, knon
838
!
839
   print *,'wx_pbl_check: Kech_h, Kech_q ', Kech_h(j), Kech_q(j)
840
!
841
  Ta(j) = (AcoefT(j) + BcoefT(j)*phiT_b9(j)*dtime)/C_p(j)
842
  qa(j) = AcoefQ(j) + BcoefQ(j)*phiQ_b9(j)*dtime
843
    print *, 'wx_pbl_check: j, Ta, qa ', Ta(j), qa(j)
844
!
845
  qsat_bs(j) = qsat01(j) + dqsatdT01(j)*(Ts_b9(j)-Ts0_b9(j))
846
!
847
   print *,'wx_pbl_check: qsat01, qsat_bs ', j,qsat01(j), qsat_bs(j)
848
!
849
  Ts_x(j) = Ts_b9(j) - sigw(j)*dTs9(j)
850
  Ts_w(j) = Ts_b9(j) + sigx(j)*dTs9(j)
851
    print *, 'wx_pbl_check: j, Ts_b9, Ts_w, Ts_x ', j, Ts_b9(j), Ts_w(j), Ts_x(j)
852
!
853
  qsat_x(j) = qsat0_x(j) + dqsatdT0_x(j)*(Ts_x(j)-Ts0_x(j))
854
  qsat_w(j) = qsat0_w(j) + dqsatdT0_w(j)*(Ts_w(j)-Ts0_w(j))
855
!
856
   print *,'wx_pbl_check: qsat0_w, qsat0_x, qsat_w, qsat_x ', qsat0_w(j), qsat0_x(j), qsat_w(j), qsat_x(j)
857
!
858
  T1_x(j) = (AcoefT_x(j) + BcoefT_x(j)*phiT_x9(j)*dtime) / C_p(j)
859
  T1_w(j) = (AcoefT_w(j) + BcoefT_w(j)*phiT_w9(j)*dtime) / C_p(j)
860
    print *, 'wx_pbl_check: j, T1_w, T1_x ', j, T1_w(j), T1_x(j)
861
!
862
  q1_x(j) = AQ_x(j) + BQ_x(j)*phiQ_x9(j)*dtime
863
  q1_w(j) = AQ_w(j) + BQ_w(j)*phiQ_w9(j)*dtime
864
    print *, 'wx_pbl_check: j, q1_w, q1_x ', j, q1_w(j), q1_x(j)
865
!
866
   qsatsrf_x(j) = beta(j)*qsat_x(j)
867
   qsatsrf_w(j) = beta(j)*qsat_w(j)
868
   qsatsrf_b(j) = sigw(j)*qsatsrf_w(j) + sigx(j)*qsatsrf_x(j)
869
!
870
   dqsatsrf1 = qsatsrf_w(j) - qsatsrf_x(j)
871
    print *, 'wx_pbl_check: j, qsatsrf_w, qsatsrf_x, dqsatsrf1, dqsatsrf9 ', &
872
                          qsatsrf_w(j), qsatsrf_x(j), dqsatsrf1, dqsatsrf9(j)
873
!
874
   qsurf_x(j) = (1-beta(j))*q1_x(j) + qsatsrf_x(j)
875
   qsurf_w(j) = (1-beta(j))*q1_w(j) + qsatsrf_w(j)
876
    print *, 'wx_pbl_check: j, qsurf_w, qsurf_x ', j, qsurf_w(j), qsurf_x(j)
877
!
878
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
879
!  Test qsat01 = qsat0    et   dqsatdT01 = dqsatdT0
880
!------------------------------------------------------------------------------------------------------
881
   print *, 'wx_pbl_check: j, qsat01(j), qsat0(j) ', j, qsat01(j), qsat0(j)
882
   print *, 'wx_pbl_check: j, dqsatdT01(j), dqsatdT0(j) ', j, dqsatdT01(j), dqsatdT0(j)
883
!
884
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
885
!  Test Kexh_Q_sw = Kech_q_w/(1.-beta*Kech_q_w*BcoefQ)   Kexh_Q_sx = Kech_q_x/(1.-beta*Kech_q_x*BcoefQ)
886
!------------------------------------------------------------------------------------------------------
887
  Kech_Q_sx1 = Kech_q_x(j)/(1.-beta(j)*Kech_q_x(j)*BQ_x(j)*dtime)
888
  Kech_Q_sw1 = Kech_q_w(j)/(1.-beta(j)*Kech_q_w(j)*BQ_w(j)*dtime)
889
    print *, 'wx_pbl_check: j, Kech_Q_sx1, Kech_Q_sx(j)', j, Kech_Q_sx1, Kech_Q_sx(j)
890
    print *, 'wx_pbl_check: j, Kech_Q_sw1, Kech_Q_sw(j)', j, Kech_Q_sw1, Kech_Q_sw(j)
891
!
892
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
893
!  Test phiT_w = Kech_h_w*C_p(j)*(T1_w(j)-Ts_w(j))        phiT_x = Kech_h_x*C_p(j)*(T1_x(j)-Ts_x(j))
894
!-----------------------------------------------------
895
    phiT_x1 = Kech_h_x(j)*C_p(j)*(T1_x(j)-Ts_x(j))
896
    phiT_w1 = Kech_h_w(j)*C_p(j)*(T1_w(j)-Ts_w(j))
897
!
898
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
899
!  Test phiT_w = Kech_T_pw*(AcoefT_w(j)-C_p(j)*Ts_w(j))   phiT_x = Kech_T_px*(AcoefT_x(j)-C_p(j)*Ts_x(j))
900
!-----------------------------------------------------
901
    phiT_x2 = Kech_T_px(j)*(AcoefT_x(j)-C_p(j)*Ts_x(j))
902
    phiT_w2 = Kech_T_pw(j)*(AcoefT_w(j)-C_p(j)*Ts_w(j))
903
    print *, 'wx_pbl_check: j, phiT_w1, phiT_w2, phiT_w9 ', j, phiT_w1, phiT_w2, phiT_w9(j)
904
    print *, 'wx_pbl_check: j, phiT_x1, phiT_x2, phiT_x9 ', j, phiT_x1, phiT_x2, phiT_x9(j)
905
!
906
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
907
!  Test phiq_w = Kech_q_w ( beta q1_w - qsatsrf_w)    phiq_x = Kech_q_x ( beta q1_x - qsatsrf_x))
908
!--------------------------------------------------------------
909
    phiq_x1 = Kech_q_x(j)*( beta(j)*q1_x(j) - qsatsrf_x(j))
910
    phiq_w1 = Kech_q_w(j)*( beta(j)*q1_w(j) - qsatsrf_w(j))
911
!
912
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
913
!  Test  phiq_w = Kech_Q_sw (beta AQ_w -qsatsrf_w)     phiq_x = Kech_Q_sx (beta AQ_x -qsatsrf_x)
914
!--------------------------------------------------------------
915
    phiq_x2 = Kech_Q_sx(j)*(beta(j)*AQ_x(j) -qsatsrf_x(j))
916
    phiq_w2 = Kech_Q_sw(j)*(beta(j)*AQ_w(j) -qsatsrf_w(j))
917
!
918
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
919
!  Test phiq_w = Kech_q_w ( q1_w - qsurf_w)    phiq_x = Kech_q_x ( q1_x - qsurf_x))
920
!--------------------------------------------------------------
921
    phiq_x3 = Kech_q_x(j)*( q1_x(j) - qsurf_x(j))
922
    phiq_w3 = Kech_q_w(j)*( q1_w(j) - qsurf_w(j))
923
    print *, 'wx_pbl_check: j, phiQ_w1, phiQ_w2, phiQ_w3, phiQ_w9 ', j, phiQ_w1, phiQ_w2, phiQ_w3, phiQ_w9(j)
924
    print *, 'wx_pbl_check: j, phiQ_x1, phiQ_x2, phiQ_x3, phiQ_x9 ', j, phiQ_x1, phiQ_x2, phiQ_x3, phiQ_x9(j)
925
!
926
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
927
!  Test phiT_b = Kech_h C_p (Ta - Ts_b)
928
!--------------------------------------------------------------
929
   phiT_b1 = Kech_h(j)*C_p(j)*(Ta(j) - Ts_b9(j))
930
   print *, 'wx_pbl_check: j, phiT_b1, PhiT_b9 ', j, phiT_b1, PhiT_b9(j)
931
!
932
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
933
!  Test phiQ_b = beta Kech_q (qa - qsat_bs)
934
!--------------------------------------------------------------
935
   evap_pot = Kech_q(j)*(qa(j) - qsat_bs(j))
936
   phiQ_b1 = beta(j)*Kech_q(j)*(qa(j) - qsat_bs(j))
937
   print *, 'wx_pbl_check: j, beta, evap_pot, phiQ_b1, PhiQ_b9 ', j, beta(j), evap_pot, phiQ_b1, PhiQ_b9(j)
938
!
939
!
940
          ENDDO  ! j = 1, knon
941
942
        ENDIF   ! (prt_level >=10 )
943
!--------------------------------------------------------------------------------------------------
944
945
        RETURN
946
947
END SUBROUTINE wx_pbl_check
948
949
SUBROUTINE wx_pbl_dts_check( knon, dtime, ypplay, ypaprs, &
950
                               sigw, beta, iflag_split,   &
951
                               Ts0_b9, dTs09,   &
952
                               qs_b9, Ts_b9,  &                         ! yqsurf, Tsurf_new
953
                               dqsatsrf9, dTs9, delta_qsurf9,   &
954
                               AcoefT_x, AcoefT_w, &
955
                               BcoefT_x, BcoefT_w, &
956
                               AcoefT0, AcoefQ0, BcoefT0, BcoefQ0, &
957
                               AcoefT,  AcoefQ,  BcoefT,  BcoefQ, &
958
                               HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn, &
959
                               phiT0_b9, dphiT09, phiQ0_b9, dphiQ09, Rn0_b9, dRn09, &
960
                               g_T, g_Q, &
961
                               Gamma_phiT, Gamma_phiQ, &
962
                               dTs_ins, dqsatsrf_ins, &
963
                               phiT_b9, phiQ_b9,  &
964
                               phiT_x9, phiT_w9, &
965
                               phiQ_x9, phiQ_w9  &
966
                               )
967
!
968
969
    USE wx_pbl_var_mod
970
971
    USE print_control_mod, ONLY: prt_level,lunout
972
!
973
    INCLUDE "YOMCST.h"
974
    INCLUDE "FCTTRE.h"
975
    INCLUDE "YOETHF.h"
976
!
977
    INTEGER,                      INTENT(IN)        :: knon         ! number of grid cells
978
    REAL,                         INTENT(IN)        :: dtime        ! time step size (s)
979
    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypplay       ! mid-layer pressure (Pa)
980
    REAL, DIMENSION(knon,klev),   INTENT(IN)        :: ypaprs       ! pressure at layer interfaces (pa)
981
    REAL, DIMENSION(knon),        INTENT(IN)        :: sigw         ! cold pools fractional area
982
    REAL, DIMENSION(knon),        INTENT(IN)        :: beta         ! aridity factor
983
    INTEGER,                      INTENT(IN)        :: iflag_split
984
    REAL, DIMENSION(knon),        INTENT(IN)        :: Ts0_b9, dTs09
985
    REAL, DIMENSION(knon),        INTENT(IN)        :: qs_b9, Ts_b9         ! yqsurf, Tsurf_new
986
    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs9, dqsatsrf9
987
    REAL, DIMENSION(knon),        INTENT(IN)        :: delta_qsurf9
988
    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT_x, AcoefT_w
989
    REAL, DIMENSION(knon),        INTENT(IN)        :: BcoefT_x, BcoefT_w
990
    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT0, AcoefQ0, BcoefT0, BcoefQ0
991
!
992
    REAL, DIMENSION(knon),        INTENT(IN)        :: AcoefT, AcoefQ, BcoefT, BcoefQ
993
    REAL, DIMENSION(knon),        INTENT(IN)        :: HTphiT_b, dd_HTphiT, HTphiQ_b, dd_HTphiQ, HTRn_b, dd_HTRn
994
    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT0_b9, dphiT09, phiQ0_b9, dphiQ09, Rn0_b9, dRn09
995
    REAL, DIMENSION(knon),        INTENT(IN)        :: g_T, g_Q
996
    REAL, DIMENSION(knon),        INTENT(IN)        :: Gamma_phiT, Gamma_phiQ
997
    REAL, DIMENSION(knon),        INTENT(IN)        :: dTs_ins, dqsatsrf_ins
998
    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT_b9, phiQ_b9
999
    REAL, DIMENSION(knon),        INTENT(IN)        :: phiT_x9, phiT_w9
1000
    REAL, DIMENSION(knon),        INTENT(IN)        :: phiQ_x9, phiQ_w9
1001
!
1002
!! Local variables
1003
    INTEGER                    :: j
1004
    REAL, DIMENSION(knon)      :: sigx       ! fractional area of (x) region
1005
    REAL, DIMENSION(knon)      :: AcoefT_b, AcoefQ_b   ! mean values of AcoefT and AcoefQ
1006
    REAL                       :: zzt, zzq, zzqsat
1007
    REAL                       :: zdelta, zcvm5, zcor, qsat
1008
    REAL, DIMENSION(knon)      :: qsat_w, qsat_x
1009
    REAL, DIMENSION(knon)      :: Ts_x, Ts_w, qs_x, qs_w
1010
    REAL, DIMENSION(knon)      :: T1_x, T1_w, q1_x, q1_w
1011
    REAL, DIMENSION(knon)      :: Rn_x, Rn_w
1012
    REAL, DIMENSION(knon)      :: Rn_b, dRn
1013
    REAL, DIMENSION(knon)      :: phiQ0_x, phiQ0_w
1014
    REAL, DIMENSION(knon)      :: Ta, qa
1015
    REAL, DIMENSION(knon)      :: err_phiT_w, err_phiT_x
1016
    REAL, DIMENSION(knon)      :: err_phiq_w, err_phiq_x
1017
    REAL, DIMENSION(knon)      :: err_phiT_b
1018
    REAL, DIMENSION(knon)      :: err_phiQ_b
1019
    REAL, DIMENSION(knon)      :: err2_phiT_b
1020
    REAL                       :: T1A_x, T1A_w, q1A_x, q1A_w
1021
    REAL                       :: qsatsrf_w, qsatsrf_x, qsatsrfb, qsbA
1022
    REAL                       :: dphiT, dphiQ
1023
    REAL                       :: dphiT_H, dphiQ_H
1024
    REAL                       :: phiQ_pot
1025
    REAL                       :: phiQ_w_m_phiQ0_w
1026
    REAL                       :: phiQ_x_m_phiQ0_x
1027
    REAL                       :: dphiQ_m_dphiQ0
1028
    REAL                       :: dphiT_m_dphiT0
1029
    REAL                       :: dRN_m_dRn0
1030
    REAL                       :: phiTb_m_phiT0b
1031
1032
!----------------------------------------------------------------------------
1033
! Equations to be checked:
1034
! -----------------------
1035
!  Input : Ts0_b, dTs0, Ts_b, dTs, qsatsrf_b, dqsatsrf,
1036
!          phiT_b, phiQ_b, phiT_w, phiT_x, phiQ_w, phiQ_x,
1037
!
1038
!          AcoefT, AcoefQ, AcoefT_w, AcoefQ_w, AcoefT_x, AcoefQ_x,
1039
!          BcoefT, BcoefQ, BcoefT_w, BcoefQ_w, BcoefT_x, BcoefQ_x
1040
!
1041
!  Ts_w = Ts_b + sigx dTs                                 Ts_x = Ts_b - sigw dTs
1042
!  T1_w = AcoefT_w + BcoefT_w phiT_w Delta t              T1_x = AcoefT_x + BcoefT_x phiT_x Delta t
1043
!  q1_w = AcoefQ_w + BcoefQ_w phiQ_w Delta t              q1_x = AcoefQ_x + BcoefQ_x phiQ_x Delta t
1044
!  phiT_w = Kech_h_w ( T1_w - Ts_w)                       phiT_x = Kech_h_x ( T1_x - Ts_x)
1045
!  phiq_w = beta Kech_h_w ( q1_w - qsat(Ts_w))            phiq_x = beta Kech_h_x ( q1_x - qsat(Ts_x))
1046
!  phiT_b = sigw phiT_w + sigx phiT_x                     dphiT = phiT_w - phiT_x
1047
!  phiQ_b = sigw phiQ_w + sigx phiQ_x                     dphiQ = phiQ_w - phiQ_x
1048
!  Ts_b = sigw Ts_w + sigx Ts_x                           dTs = Ts_w - Ts_x
1049
!  Ta = AcoefT + BcoefT phiT_b Delta t
1050
!  qa = AcoefQ + BcoefQ phiQ_b Delta t
1051
!  phiT_b = Kech_h (Ta - Ts_b)
1052
!  phiQ_b = beta Kech_h (qa - qsat(Ts_b))
1053
!  dTs = sqrt(tau)/I (dphit + L_v dphiq + dR)
1054
1055
!----------------------------------------------------------------------------
1056
!
1057
!!
1058
        sigx(:) = 1.-sigw(:)
1059
        AcoefT_b(1:knon) = AcoefT_x(1:knon) + sigw(1:knon)*dd_AT(1:knon)
1060
        AcoefQ_b(1:knon) = AQ_x(1:knon) + sigw(1:knon)*dd_AQ(1:knon)
1061
1062
   IF (prt_level >=10 ) THEN
1063
    print *,'->wx_pbl_dts_check, HTphiT_b, HTphiQ_b, HTRn_b ', &
1064
                             HTphiT_b, HTphiQ_b, HTRn_b
1065
    print *,'->wx_pbl_dts_check, dd_HTphiT, dd_HTphiQ, dd_HTRn ', &
1066
                             dd_HTphiT, dd_HTphiQ, dd_HTRn
1067
   ENDIF ! (prt_level >=10 )
1068
!
1069
! Compute the three qsat and dqsatdTs
1070
! ---------------------------------------------
1071
!!      print *,' AAAA wx_pbl_dts_check, C_p(j), qsat0(j), Ts0(j) : ',  &
1072
!!                                      (C_p(j), qsat0(j), Ts0(j), j = 1,knon)
1073
!
1074
!
1075
!--------------------------------------------------------------------------------------------------
1076
        IF (prt_level >=10 ) THEN
1077
!
1078
          DO j = 1, knon
1079
  Ts_x(j) = Ts_b9(j) - sigw(j)*dTs9(j)
1080
  Ts_w(j) = Ts_b9(j) + sigx(j)*dTs9(j)
1081
    print *, 'wx_pbl_dts_check: j, Ts_b9, Ts_w, Ts_x ', j, Ts_b9(j), Ts_w(j), Ts_x(j)
1082
!
1083
  qsat_x(j) = qsat0_x(j) + dqsatdT0_x(j)*(Ts_x(j)-Ts0_x(j))
1084
  qsat_w(j) = qsat0_w(j) + dqsatdT0_w(j)*(Ts_w(j)-Ts0_w(j))
1085
!
1086
  T1_x(j) = (AcoefT_x(j) + BcoefT_x(j)*phiT_x9(j)*dtime) / C_p(j)
1087
  T1_w(j) = (AcoefT_w(j) + BcoefT_w(j)*phiT_w9(j)*dtime) / C_p(j)
1088
    print *, 'wx_pbl_dts_check: j, T1_w, T1_x ', j, T1_w(j), T1_x(j)
1089
!
1090
  q1_x(j) = AQ_x(j) + BQ_x(j)*phiQ_x9(j)*dtime
1091
  q1_w(j) = AQ_w(j) + BQ_w(j)*phiQ_w9(j)*dtime
1092
    print *, 'wx_pbl_dts_check: j, q1_w, q1_x ', j, q1_w(j), q1_x(j)
1093
!
1094
    Rn_x(j) = eps_1*Rsigma*T1_x(j)**4 - Rsigma*Ts_x(j)**4
1095
    Rn_w(j) = eps_1*Rsigma*T1_w(j)**4 - Rsigma*Ts_w(j)**4
1096
    Rn_b(j) = sigw(j)*Rn_w(j) + sigx(j)*Rn_x(j)
1097
    dRn(j) = dRn09(j) - ( HTRn_b(j) &
1098
                        +(sigx(j)-sigw(j))*dd_HTRn(j) &
1099
                        -sigw(j)*sigx(j)*dd_HTRn(j)*dd_HTphiT(j)/HTphiT_b(j) &
1100
                       )*(dTs9(j)-dTs09(j)) &
1101
                    + dd_HTRn(j)/HTphiT_b(j)*(phiT_b9(j)-phiT0_b9(j))
1102
!
1103
          print *,'wx_pbl_dts_check, dphiT, L_v*dphiQ, dRn, dTs ', &
1104
              phiT_w9(j)-phiT_x9(j), L_v(j)*(phiQ_w9(j)-phiQ_x9(j)), dRn(j), dTs9(j)
1105
!
1106
  phiQ0_x(j) = PhiQ0_b9(j) - sigw(j)*dphiQ09(j)
1107
  phiQ0_w(j) = PhiQ0_b9(j) + sigx(j)*dphiQ09(j)
1108
!
1109
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1110
!  Test phiQ_w-phiQ0_w = -beta*Kech_Q_sw*dqsatdT_w*(Ts_w-Ts0_w)
1111
!--------------------------------------------------------------
1112
  print *,'wx_pbl_dts_check: beta(j), Kech_Q_sw(j), dqsatdT0_w(j), Ts_w(j), Ts0_w(j) ', &
1113
                         beta(j), Kech_Q_sw(j), dqsatdT0_w(j), Ts_w(j), Ts0_w(j)
1114
  phiQ_w_m_phiQ0_w = -beta(j)*Kech_Q_sw(j)*dqsatdT0_w(j)*(Ts_w(j)-Ts0_w(j))
1115
    print *,'wx_pbl_dts_check: j, phiQ_w9-phiQ0_w, phiQ_w_m_phiQ0_w ', &
1116
                           j, phiQ_w9(j)-phiQ0_w(j), phiQ_w_m_phiQ0_w
1117
!
1118
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1119
!  Test phiQ_x-phiQ0_x = -beta*Kech_Q_sx*dqsatdT_x*(Ts_x-Ts0_x)
1120
!--------------------------------------------------------------
1121
  phiQ_x_m_phiQ0_x = -beta(j)*Kech_Q_sx(j)*dqsatdT0_x(j)*(Ts_x(j)-Ts0_x(j))
1122
    print *,'wx_pbl_dts_check: j, phiQ_x9-phiQ0_x, phiQ_x_m_phiQ0_x ', &
1123
                           j, phiQ_x9(j)-phiQ0_x(j), phiQ_x_m_phiQ0_x
1124
!
1125
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1126
!  Test dphiT-dphiT0 = -(HTphiT_b+(sigx-sigw)*dd_HTphiT)*(dTs-dTs0) - dd_HTphiT*(Ts_b-Ts0_b)
1127
!-------------------------------------------------------------------------------------------
1128
 dphiT = phiT_w9(j) - phiT_x9(j)
1129
 dphiT_m_dphiT0 = -(HTphiT_b(j)+(sigx(j)-sigw(j))*dd_HTphiT(j))*(dTs9(j)-dTs09(j)) &
1130
                  - dd_HTphiT(j)*(Ts_b9(j)-Ts0_b9(j))
1131
 print *,'wx_pbl_dts_check: j, dphiT-dphiT09, dphiT_m_dphiT0 ',j, dphiT-dphiT09(j), dphiT_m_dphiT0
1132
!
1133
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1134
!  Test dphiQ-dphiQ0 = -(HTphiQ_b+(sigx-sigw)*dd_HTphiQ)*(dTs-dTs0) - dd_HTphiQ*(Ts_b-Ts0_b)
1135
!-------------------------------------------------------------------------------------------
1136
 dphiQ = phiQ_w9(j) - phiQ_x9(j)
1137
 dphiQ_m_dphiQ0 = -(HTphiQ_b(j)+(sigx(j)-sigw(j))*dd_HTphiQ(j))*(dTs9(j)-dTs09(j)) &
1138
                  - dd_HTphiQ(j)*(Ts_b9(j)-Ts0_b9(j))
1139
 print *,'wx_pbl_dts_check: j, dphiQ-dphiQ09, dphiQ_m_dphiQ0 ',j, dphiQ-dphiQ09(j), dphiQ_m_dphiQ0
1140
!
1141
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1142
!  Test dRn-dRn0 = -(HTRn_b+(sigx-sigw)*dd_HTRn)*(dTs-dTs0) - dd_HTRn*(Ts_b-Ts0_b)
1143
!-------------------------------------------------------------------------------------------
1144
 dRn_m_dRn0 = -(HTRn_b(j)+(sigx(j)-sigw(j))*dd_HTRn(j))*(dTs9(j)-dTs09(j)) &
1145
                  - dd_HTRn(j)*(Ts_b9(j)-Ts0_b9(j))
1146
 print *,'wx_pbl_dts_check: j, dRn-dRn09, dRn_m_dRn0 ',j, dRn-dRn09(j), dRn_m_dRn0
1147
!
1148
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1149
!  Test phiT_b-phiT0_b = -sigx*sigw*dd_HTphiT*(dTs-dTs0) - HTphiT_b*(Ts_b-Ts0_b)
1150
!-------------------------------------------------------------------------------
1151
   phiTb_m_phiT0b = -sigx(j)*sigw(j)*dd_HTphiT(j)*(dTs9(j)-dTs09(j)) - HTphiT_b(j)*(Ts_b9(j)-Ts0_b9(j))
1152
   print *,'wx_pbl_dts_check: j, phiT_b9-phiT0_b9, phiTb_m_phiT0b ',j ,phiT_b9(j)-phiT0_b9(j), phiTb_m_phiT0b
1153
!
1154
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1155
!  Test phiT_w, phiT_x, dphiT from HTphiT
1156
!------------------------------------------
1157
!  phiT_w = Kech_h_w C_p ( T1_w - Ts_w)                   phiT_x = Kech_h_x C_p ( T1_x - Ts_x)
1158
  err_phiT_x(j) = Kech_h_x(j)*C_p(j)*(T1_x(j) - Ts_x(j)) - phiT_x9(j)
1159
  err_phiT_w(j) = Kech_h_w(j)*C_p(j)*(T1_w(j) - Ts_w(j)) - phiT_w9(j)
1160
    print *, 'wx_pbl_dts_check: j, phiT_w9, phiT_x9, err_phiT_w, err_phiT_x ',   &
1161
                            j, phiT_w9(j), phiT_x9(j), err_phiT_w(j), err_phiT_x(j)
1162
  dphiT = phiT_w9(j) - phiT_x9(j)
1163
  dphiT_H = dphiT09(j) - ( HTphiT_b(j) &
1164
                            +(sigx(j)-sigw(j))*dd_HTphiT(j) &
1165
                            -sigw(j)*sigx(j)*dd_HTphiT(j)*dd_HTphiT(j)/HTphiT_b(j) &
1166
                           )*(dTs9(j)-dTs09(j)) &
1167
                         + dd_HTphiT(j)/HTphiT_b(j)*(phiT_b9(j)-phiT0_b9(j))
1168
    print *,'wx_pbl_dts_check: j, dphiT, dphiT_H ', j, dphiT, dphiT_H
1169
1170
!
1171
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1172
!  Test phiq_w, phiq_x, dphiq from HTphiq
1173
!------------------------------------------
1174
!
1175
!  phiq_w = beta Kech_q_w ( q1_w - qsat(Ts_w))            phiq_x = beta Kech_q_x ( q1_x - qsat(Ts_x))
1176
  err_phiq_x(j) = beta(j)*Kech_q_x(j)*( q1_x(j) - qsat_x(j)) - phiq_x9(j)
1177
  err_phiq_w(j) = beta(j)*Kech_q_w(j)*( q1_w(j) - qsat_w(j)) - phiq_w9(j)
1178
  dphiQ = phiQ_w9(j) - phiQ_x9(j)
1179
  dphiQ_H = dphiQ09(j) - ( HTphiQ_b(j) &
1180
                            +(sigx(j)-sigw(j))*dd_HTphiQ(j) &
1181
                            -sigw(j)*sigx(j)*dd_HTphiQ(j)*dd_HTphiT(j)/HTphiT_b(j) &
1182
                           )*(dTs9(j)-dTs09(j)) &
1183
                         + dd_HTphiQ(j)/HTphiT_b(j)*(phiT_b9(j)-phiT0_b9(j))
1184
    print *,'wx_pbl_dts_check: j, dphiQ, dphiQ_H ', j, dphiQ, dphiQ_H
1185
!
1186
!  phiT_b = sigw phiT_w + sigx phiT_x                     dphiT = phiT_w - phiT_x
1187
  err_phiT_b(j) = sigw(j)*phiT_w9(j) + sigx(j)*phiT_x9(j) - phiT_b9(j)
1188
!
1189
!  phiQ_b = sigw phiQ_w + sigx phiQ_x                     dphiQ = phiQ_w - phiQ_x
1190
  err_phiQ_b(j) = sigw(j)*phiQ_w9(j) + sigx(j)*phiQ_x9(j) - phiQ_b9(j)
1191
!
1192
!  Ta = AcoefT + BcoefT phiT_b Delta t
1193
!  phiT_b = Kech_h C_p (Ta - Ts_b)
1194
  Ta(j) = (AcoefT(j) + BcoefT(j)*phiT_b9(j)*dtime) / C_p(j)
1195
  err2_phiT_b(j) = Kech_h(j)*C_p(j)*(Ta(j) - Ts_b9(j)) - phiT_b9(j)
1196
    print *, 'wx_pbl_dts_check: j, Ta, phiT_b9, err2_phiT_b ',   &
1197
                            j, Ta(j), phiT_b9(j), err2_phiT_b(j)
1198
!
1199
          ENDDO  ! j = 1, knon
1200
1201
        ENDIF   ! (prt_level >=10 )
1202
!--------------------------------------------------------------------------------------------------
1203
        RETURN
1204
1205
END SUBROUTINE wx_pbl_dts_check
1206
1207
SUBROUTINE wx_evappot(knon, q1, Ts, evap_pot)
1208
1209
    USE wx_pbl_var_mod
1210
1211
    INTEGER,                      INTENT(IN)        :: knon     ! number of grid cells
1212
    REAL, DIMENSION(knon),        INTENT(IN)        :: q1       ! specific humidity in layer 1
1213
    REAL, DIMENSION(knon),        INTENT(IN)        :: Ts       ! surface temperature
1214
!
1215
    REAL, DIMENSION(knon),        INTENT(OUT)       :: evap_pot ! potential evaporation
1216
!
1217
    INTEGER                   :: j
1218
    REAL                      :: qsat_bs
1219
!
1220
 DO j = 1,knon
1221
   evap_pot(j) = Kech_q(j)*(qsat0(j)+dqsatdT0(j)*(Ts(j)-Ts0(j))-q1(j))
1222
!
1223
  qsat_bs = qsat0(j)+dqsatdT0(j)*(Ts(j)-Ts0(j))
1224
!!  print *,'wx_evappot : Kech_q, qsat_bs, qa, evap_pot ', Kech_q(j), qsat_bs, q1(j), evap_pot(j)
1225
 ENDDO
1226
!
1227
 RETURN
1228
END SUBROUTINE wx_evappot
1229
1230
END MODULE wx_pbl_mod
1231