source: LMDZ6/branches/contrails/libf/phylmd/wx_pbl_mod.f90 @ 5428

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

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