source: LMDZ6/trunk/libf/phylmd/wx_pbl_mod.F90 @ 5209

Last change on this file since 5209 was 3906, checked in by jyg, 3 years ago

1/ Bug fix in wx_pbl_mod.F90

2/ New pbl_split option:
iflag_pbl_split=3 ==> iflag_split=1 over ocean,

iflag_split=0 everywhere else.

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