Context Navigation

wx_pbl_mod.F90 @ 3893

Last change on this file since 3893 was 3888, checked in by jyg, 4 years ago

New provisional version of the splitting of the
diffusive boundary layer into inwake and offwake
PBLs. The splitting of the diffuse BL should NOT
be activated yet for general purpose simulations.

The splitting is activated by:
mod(iflag_pbl_split,10)=1 for the option with
fixed surface temperature and
mod(iflag_pbl_split,10)=2 for the option with
coupled surface temperature.

iflag_pbl_split=0 ==> no splittingat all.
iflag_pbl_split=10 ==> splitting of thermals.
iflag_pbl_split=11 ==> splitting of thermals and
of vertical diffusion (fixed surf. temp.).
iflag_pbl_split=12 ==> splitting of thermals and
of vertical diffusion (coupled surf. temp.).

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

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

Download in other formats:

Original Format