Context Navigation

wx_pbl_mod.F90 @ 5151

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

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

Download in other formats:

Original Format