Context Navigation

wx_pbl_mod.f90 @ 5884

Last change on this file since 5884 was 5884, checked in by yann meurdesoif, 3 days ago

GPU port of wx_evappot, wx_pbl_split, wx_pbl_check, wx_pbl_dts_check

File size: 61.0 KB

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

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

Download in other formats:

Original Format