Context Navigation

wx_pbl_mod.f90 @ 5428

Last change on this file since 5428 was 5285, checked in by abarral, 8 weeks ago
As discussed internally, remove generic ONLY: ... for new _mod_h modules
File size: 60.1 KB

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

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

Download in other formats:

Original Format