source: LMDZ6/branches/Amaury_dev/libf/phylmd/wx_pbl_mod.F90 @ 5725

Last change on this file since 5725 was 5160, checked in by abarral, 16 months ago

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