Context Navigation

wx_pbl_mod.f90 @ 5274

Last change on this file since 5274 was 5274, checked in by abarral, 9 hours ago
Replace yomcst.h by existing module
File size: 64.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format