Context Navigation

lmdz_wake_dadv.f90 @ 5833

Last change on this file since 5833 was 5805, checked in by fhourdin, 2 months ago
Details poches
File size: 21.1 KB

Line
1	MODULE lmdz_wake_dadv
2	PUBLIC wake_dadv
3	CONTAINS
4
5	SUBROUTINE wake_dadv(klon, klev, dtime, ph, ppi, wk_adv, kupper, &
6	deltomg, dp_deltomg, sigmaw, dsigspread, &
7	thw, thx, qw, qx, &
8	d_deltat_dadv, d_deltaq_dadv, d_tb_dadv, d_qb_dadv)
9
10	USE lmdz_wake_ini , ONLY : flag_dadv_implicit
11
12	IMPLICIT NONE
13
14	INTEGER, INTENT(IN) :: klon, klev
15	REAL, INTENT(IN) :: dtime
16	REAL, DIMENSION (klon, klev+1), INTENT(IN) :: ph
17	REAL, DIMENSION (klon, klev), INTENT(IN) :: ppi
18	LOGICAL, DIMENSION (klon), INTENT(IN) :: wk_adv
19	INTEGER, DIMENSION (klon), INTENT(IN) :: kupper
20	REAL, DIMENSION (klon, klev), INTENT(IN) :: deltomg
21	REAL, DIMENSION (klon, klev), INTENT(IN) :: dp_deltomg
22	REAL, DIMENSION (klon), INTENT(IN) :: sigmaw
23	REAL, DIMENSION (klon), INTENT(IN) :: dsigspread
24	REAL, DIMENSION (klon, klev), INTENT(IN) :: thw ! component # 1
25	REAL, DIMENSION (klon, klev), INTENT(IN) :: thx ! component # 2
26	REAL, DIMENSION (klon, klev), INTENT(IN) :: qw ! component # 1
27	REAL, DIMENSION (klon, klev), INTENT(IN) :: qx ! component # 2
28
29	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltat_dadv
30	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_deltaq_dadv
31	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_tb_dadv
32	REAL, DIMENSION (klon, klev), INTENT(OUT) :: d_qb_dadv
33
34	! Internal variables
35	INTEGER :: i, k
36	REAL, DIMENSION (klon, klev) :: entr_s ! entrainment into wakes due to spread
37	REAL, DIMENSION (klon, klev) :: thb, qb
38	REAL, DIMENSION (klon, klev) :: delta_th, delta_q
39
40	! Tests
41
42	! Arrays used in the implicit scheme
43	REAL, DIMENSION (klon) :: rr11, rr12, rr21, rr22
44
45	REAL, DIMENSION (klon, klev) :: aa11, aa12, aa21, aa22
46	REAL, DIMENSION (klon, klev) :: bb11, bb12, bb21, bb22
47	REAL, DIMENSION (klon, klev) :: cc11, cc12, cc21, cc22
48
49	REAL, DIMENSION (klon, klev) :: alpha11, alpha12, alpha21, alpha22
50	REAL, DIMENSION (klon, klev) :: beta11, beta12, beta21, beta22
51	REAL, DIMENSION (klon, klev) :: gamma11, gamma12, gamma21, gamma22
52	REAL, DIMENSION (klon, klev) :: ai11, ai12, ai21, ai22 ! inverse of alpha
53
54	REAL, DIMENSION (klon, klev) :: xt1, xt2, xq1, xq2
55	REAL, DIMENSION (klon, klev) :: yt1, yt2, yq1, yq2
56	REAL, DIMENSION (klon, klev) :: zt1, zt2, zq1, zq2
57	REAL, DIMENSION (klon, klev) :: th1, th2, q1, q2
58
59	REAL :: coef, det
60
61	REAL, DIMENSION (klon,klev) :: xt1inv, xt2inv, xq1inv, xq2inv
62
63	! Arrays used in the explicit scheme (vertical gradients)
64	REAL, DIMENSION (klon, klev) :: d_thx, d_qx
65	REAL, DIMENSION (klon, klev) :: d_thw, d_qw
66	REAL, DIMENSION (klon, klev) :: d_dth, d_dq
67
68	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
69
70	! print *,'ZZZwake_dadv_IN wk_adv(1) ', wk_adv(1)
71	! print *,'ZZZwake_dadv_IN kupper(1) ', kupper(1)
72	! print *,'ZZZwake_dadv_IN k, thw(1,k), thx(1,k) ', (k, thw(1,k), thx(1,k), k = 1,3)
73	! print *,'ZZZwake_dadv_IN k, deltomg(1,k) ', (k, deltomg(1,k), k = 1,3)
74	! print *,'ZZZwake_dadv_IN k, dp_deltomg(1,k) ', (k, dp_deltomg(1,k), k = 1,3)
75	! print *,'ZZZwake_dadv_IN sigmaw(1) ', sigmaw(1)
76	! print *,'ZZZwake_dadv_IN dsigspread(1) ', dsigspread(1)
77
78	entr_s(:,:) = 0.
79	delta_th(:,:) = 0.
80
81	d_deltat_dadv(:,:) = 0.
82	d_deltaq_dadv(:,:) = 0.
83	d_tb_dadv(:,:) = 0.
84	d_qb_dadv(:,:) = 0.
85
86
87	rr11(:) = sigmaw(:)
88	rr12(:) = 1.-sigmaw(:)
89	rr21(:) = 1.
90	rr22(:) = -1.
91
92	DO k = 1, klev
93	DO i = 1,klon
94	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
95	thb(i,k) = rr11(i)thw(i,k)+rr12(i)thx(i,k)
96	delta_th(i,k) = rr21(i)thw(i,k)+rr22(i)thx(i,k)
97
98	qb(i,k) = rr11(i)qw(i,k)+rr12(i)qx(i,k)
99	delta_q(i,k) = rr21(i)qw(i,k)+rr22(i)qx(i,k)
100	ENDIF
101	ENDDO
102	ENDDO
103
104	DO i = 1, klon
105	entr_s(i,klev) = 0.
106	ENDDO
107
108	DO k = 1, klev-1
109	DO i = 1,klon
110	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
111	!! entr_s(i,k) = dsigspread(i) - sigmaw(i)(1.-sigmaw(i))(deltomg(i,k+1)-deltomg(i,k)) / &
112	!! (ph(i,k)-ph(i,k+1))
113
114	entr_s(i,k) = dsigspread(i) + sigmaw(i)(1.-sigmaw(i))dp_deltomg(i,k)
115	!! print *,'dadv, k, dp_deltomg(i,k), (deltomg(i,k)-deltomg(i,k+1))/(ph(i,k)-ph(i,k+1)) ', &
116	!! k, dp_deltomg(i,k), (deltomg(i,k)-deltomg(i,k+1))/(ph(i,k)-ph(i,k+1))
117
118	ENDIF
119	ENDDO
120	ENDDO
121
122	! -------------------------------------------------------------------------------------------
123	! Depending on flag_dadv_implicit, use implicit upstream scheme or explicit upstream scheme
124	! -------------------------------------------------------------------------------------------
125
126	IF (flag_dadv_implicit) THEN
127
128	! Implicit scheme : solve for d_deltat_dadv and d_tb_dadv
129	! (and similarly for d_deltaq_dadv and d_qb_dadv).
130	! The system to be inverted is block-tridiagonal with 2x2 blocks.
131	! -----------------------------------------------------------------------------------------
132
133	! Matrix indexing: Theta_w Theta_x
134	!
135	! /
136	! Theta_b \| A11 A12 \|
137	! \| \|
138	! delta_theta \| A21 A22 \|
139	! /
140	! Tridiagonal matrix
141	! /
142	! \| aa(1) bb(1) 0 \|
143	! \| cc(2) aa(2) bb(2) 0 \|
144	! \| 0 cc(3) aa(3) bb(3) \|
145	! \| \|
146	! . . . .
147	!
148	! . . . .
149	! \| \|
150	! \| cc(n-2) aa(n-2) bb(n-2) 0 \|
151	! \| 0 cc(n-1) aa(n-1) bb(n-1) \|
152	! 0 cc(n) aa(n) /
153	! -----------------------------------------------------------------------------------------
154
155	!! Building the tridiagonal matrix
156	DO i = 1,klon
157	IF (wk_adv(i)) THEN
158	k = kupper(i)
159	coef = dtime/(ph(i,k)-ph(i,k+1))
160	aa11(i,k) = rr11(i)+coefsigmaw(i)(1.-sigmaw(i))*deltomg(i,k)
161	aa12(i,k) = rr12(i)
162	aa21(i,k) = rr21(i)+coef( dsigspread(i)/sigmaw(i)(ph(i,k)-ph(i,k+1)) + &
163	(1.-sigmaw(i))*deltomg(i,k) )
164	aa22(i,k) = rr22(i)+coef(-dsigspread(i)/sigmaw(i)(ph(i,k)-ph(i,k+1)) - &
165	deltomg(i,k) )
166
167	cc11(i,k) = 0.
168	cc12(i,k) = -coefsigmaw(i)(1.-sigmaw(i))*deltomg(i,k)
169	cc21(i,k) = 0.
170	cc22(i,k) = coefsigmaw(i)deltomg(i,k)
171	ENDIF ! (wk_adv(i))
172	ENDDO
173	DO k = 2, klev-1
174	DO i = 1,klon
175	IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
176	coef = dtime/(ph(i,k)-ph(i,k+1))
177	aa11(i,k) = rr11(i)+coefsigmaw(i)(1.-sigmaw(i))*deltomg(i,k)
178	aa12(i,k) = rr12(i)+coefsigmaw(i)(1.-sigmaw(i))*deltomg(i,k+1)
179	aa21(i,k) = rr21(i)+coef( dsigspread(i)/sigmaw(i)(ph(i,k)-ph(i,k+1)) + &
180	(1.-sigmaw(i))*deltomg(i,k) )
181	aa22(i,k) = rr22(i)+coef(-dsigspread(i)/sigmaw(i)(ph(i,k)-ph(i,k+1)) + &
182	(1.-sigmaw(i))*(deltomg(i,k+1)-deltomg(i,k)) - &
183	sigmaw(i)*deltomg(i,k) )
184
185	bb11(i,k) = -coefsigmaw(i)(1.-sigmaw(i))*deltomg(i,k+1)
186	bb12(i,k) = 0.
187	bb21(i,k) = -coef(1.-sigmaw(i))deltomg(i,k+1)
188	bb22(i,k) = 0.
189
190	cc11(i,k) = 0.
191	cc12(i,k) = -coefsigmaw(i)(1.-sigmaw(i))*deltomg(i,k)
192	cc21(i,k) = 0.
193	cc22(i,k) = coefsigmaw(i)deltomg(i,k)
194	ENDIF ! (wk_adv(i) .AND. k<=kupper(i))
195	ENDDO
196	ENDDO
197	DO i = 1,klon
198	IF (wk_adv(i)) THEN
199	coef = dtime/(ph(i,1)-ph(i,2))
200	aa11(i,1) = rr11(i)+coefsigmaw(i)(1.-sigmaw(i))*deltomg(i,1)
201	aa12(i,1) = rr12(i)+coefsigmaw(i)(1.-sigmaw(i))*deltomg(i,2)
202	aa21(i,1) = rr21(i)+coef( dsigspread(i)/sigmaw(i)(ph(i,1)-ph(i,2)) + &
203	(1.-sigmaw(i))*deltomg(i,1) )
204	aa22(i,1) = rr22(i)+coef(-dsigspread(i)/sigmaw(i)(ph(i,1)-ph(i,2)) + &
205	(1.-sigmaw(i))*(deltomg(i,2)-deltomg(i,1)) - &
206	sigmaw(i)*deltomg(i,1) )
207
208	bb11(i,1) = -coefsigmaw(i)(1.-sigmaw(i))*deltomg(i,2)
209	bb12(i,1) = 0.
210	bb21(i,1) = -coef(1.-sigmaw(i))deltomg(i,2)
211	bb22(i,1) = 0.
212	ENDIF ! (wk_adv(i))
213	ENDDO
214
215	!! printing the tridiagonal matrix
216	!!! First row
217	!! k = 1
218	!! print 1789, k, aa11(1,1), aa12(1,1), bb11(1,1), bb12(1,1)
219	!! print 1789, k, aa21(1,1), aa22(1,1), bb21(1,1), bb22(1,1)
220	!!1789 FORMAT(1X, I3, 3(4X, 2E13.5))
221	!! coef = dtime/(ph(1,k)-ph(1,k+1))
222	!! print *,'rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,1), deltomg(1,2) ', &
223	!! rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,1), deltomg(1,2)
224	!!
225	!!! Rows 2 to klev-1
226	!! DO k = 2, klev-1
227	!! print 1789, k, cc11(1,k), cc12(1,k), aa11(1,k), aa12(1,k), bb11(1,k), bb12(1,k)
228	!! print 1789, k, cc21(1,k), cc22(1,k), aa21(1,k), aa22(1,k), bb21(1,k), bb22(1,k)
229	!! coef = dtime/(ph(1,k)-ph(1,k+1))
230	!! print *,'rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,k), deltomg(1,k+1) ', &
231	!! rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,k), deltomg(1,k+1)
232	!! ENDDO
233	!!
234	!!! Row klev
235	!! print 1789, klev, cc11(1,klev), cc12(1,klev), aa11(1,klev), aa12(1,klev)
236	!! print 1789, klev, cc21(1,klev), cc22(1,klev), aa21(1,klev), aa22(1,klev)
237	!! coef = dtime/(ph(1,klev)-ph(1,klev+1))
238	!! print *,'rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,klev) ', &
239	!! rr22(1), coef, dsigspread(1), sigmaw(1), deltomg(1,klev)
240
241
242	!! Downward loop
243
244	xt1(:,:) = thb(:,:)
245	xt2(:,:) = delta_th(:,:)
246	xq1(:,:) = qb(:,:)
247	xq2(:,:) = delta_q(:,:)
248
249	DO i = 1,klon
250	IF (wk_adv(i)) THEN
251	k = kupper(i)
252	alpha11(:,k)=aa11(:,k)
253	alpha12(:,k)=aa12(:,k)
254	alpha21(:,k)=aa21(:,k)
255	alpha22(:,k)=aa22(:,k)
256	beta11(:,k)=0.
257	beta12(:,k)=0.
258	beta21(:,k)=0.
259	beta22(:,k)=0.
260	yt1(i,k) = xt1(i,k)
261	yt2(i,k) = xt2(i,k)
262	yq1(i,k) = xq1(i,k)
263	yq2(i,k) = xq2(i,k)
264	ENDIF ! (wk_adv(i))
265	ENDDO
266	DO i = 1,klon
267	IF (wk_adv(i)) THEN
268	k = kupper(i)
269	det=alpha11(i,k)alpha22(i,k) - alpha12(i,k)alpha21(i,k)
270	ai11(i,k)= alpha22(i,k)/det
271	ai12(i,k)=-alpha12(i,k)/det
272	ai21(i,k)=-alpha21(i,k)/det
273	ai22(i,k)= alpha11(i,k)/det
274	zt1(i,k) = ai11(i,k)yt1(i,k) + ai12(i,k)yt2(i,k)
275	zt2(i,k) = ai21(i,k)yt1(i,k) + ai22(i,k)yt2(i,k)
276	zq1(i,k) = ai11(i,k)yq1(i,k) + ai12(i,k)yq2(i,k)
277	zq2(i,k) = ai21(i,k)yq1(i,k) + ai22(i,k)yq2(i,k)
278	ENDIF ! (wk_adv(i))
279	ENDDO
280
281	DO k = klev, 2, -1
282	DO i = 1,klon
283	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
284	gamma11(i,k) = ai11(i,k)cc11(i,k) + ai12(i,k)cc21(i,k)
285	gamma12(i,k) = ai11(i,k)cc12(i,k) + ai12(i,k)cc22(i,k)
286	gamma21(i,k) = ai21(i,k)cc11(i,k) + ai22(i,k)cc21(i,k)
287	gamma22(i,k) = ai21(i,k)cc12(i,k) + ai22(i,k)cc22(i,k)
288
289	alpha11(i,k-1) = aa11(i,k-1) - ( bb11(i,k-1)gamma11(i,k)+bb12(i,k-1)gamma21(i,k) )
290	alpha12(i,k-1) = aa12(i,k-1) - ( bb11(i,k-1)gamma12(i,k)+bb12(i,k-1)gamma22(i,k) )
291	alpha21(i,k-1) = aa21(i,k-1) - ( bb21(i,k-1)gamma11(i,k)+bb22(i,k-1)gamma21(i,k) )
292	alpha22(i,k-1) = aa22(i,k-1) - ( bb21(i,k-1)gamma12(i,k)+bb22(i,k-1)gamma22(i,k) )
293
294	beta11(i,k-1) = bb11(i,k-1)ai11(i,k)+bb12(i,k-1)ai21(i,k)
295	beta12(i,k-1) = bb11(i,k-1)ai12(i,k)+bb12(i,k-1)ai22(i,k)
296	beta21(i,k-1) = bb21(i,k-1)ai11(i,k)+bb22(i,k-1)ai21(i,k)
297	beta22(i,k-1) = bb21(i,k-1)ai12(i,k)+bb22(i,k-1)ai22(i,k)
298
299	yt1(i,k-1) = xt1(i,k-1) - ( beta11(i,k-1)yt1(i,k) +beta12(i,k-1)yt2(i,k) )
300	yt2(i,k-1) = xt2(i,k-1) - ( beta21(i,k-1)yt1(i,k) +beta22(i,k-1)yt2(i,k) )
301	yq1(i,k-1) = xq1(i,k-1) - ( beta11(i,k-1)yq1(i,k) +beta12(i,k-1)yq2(i,k) )
302	yq2(i,k-1) = xq2(i,k-1) - ( beta21(i,k-1)yq1(i,k) +beta22(i,k-1)yq2(i,k) )
303
304	det=alpha11(i,k-1)alpha22(i,k-1) - alpha12(i,k-1)alpha21(i,k-1)
305	ai11(i,k-1)= alpha22(i,k-1)/det
306	ai12(i,k-1)=-alpha12(i,k-1)/det
307	ai21(i,k-1)=-alpha21(i,k-1)/det
308	ai22(i,k-1)= alpha11(i,k-1)/det
309
310	zt1(i,k-1) = ai11(i,k-1)yt1(i,k-1)+ai12(i,k-1)yt2(i,k-1)
311	zt2(i,k-1) = ai21(i,k-1)yt1(i,k-1)+ai22(i,k-1)yt2(i,k-1)
312	zq1(i,k-1) = ai11(i,k-1)yq1(i,k-1)+ai12(i,k-1)yq2(i,k-1)
313	zq2(i,k-1) = ai21(i,k-1)yq1(i,k-1)+ai22(i,k-1)yq2(i,k-1)
314	ENDIF ! (wk_adv(i) .AND. k<=kupper(i))
315	ENDDO
316	ENDDO
317
318	!! Upward loop
319
320	DO i = 1,klon
321	IF (wk_adv(i)) THEN
322	th1(i,1) = zt1(i,1)
323	th2(i,1) = zt2(i,1)
324	q1(i,1) = zq1(i,1)
325	q2(i,1) = zq2(i,1)
326
327	d_tb_dadv(i,1) = ( rr11(i)(th1(i,1)-thw(i,1))+rr12(i)(th2(i,1)-thx(i,1)) )*ppi(i,1)
328	d_deltat_dadv(i,1) = ( rr21(i)(th1(i,1)-thw(i,1))+rr22(i)(th2(i,1)-thx(i,1)) )*ppi(i,1)
329	d_qb_dadv(i,1) = rr11(i)(q1(i,1) -qw(i,1)) +rr12(i)(q2(i,1)-qx(i,1))
330	d_deltaq_dadv(i,1) = rr21(i)(q1(i,1) -qw(i,1)) +rr22(i)(q2(i,1)-qx(i,1))
331	ENDIF ! (wk_adv(i))
332	ENDDO
333
334	DO k = 2, klev
335	DO i = 1,klon
336	IF (wk_adv(i) .AND. k<=kupper(i)) THEN
337	th1(i,k) = zt1(i,k) - ( gamma11(i,k)th1(i,k-1)+gamma12(i,k)th2(i,k-1) )
338	th2(i,k) = zt2(i,k) - ( gamma21(i,k)th1(i,k-1)+gamma22(i,k)th2(i,k-1) )
339	q1(i,k) = zq1(i,k) - ( gamma11(i,k)q1(i,k-1) +gamma12(i,k)q2(i,k-1) )
340	q2(i,k) = zq2(i,k) - ( gamma21(i,k)q1(i,k-1) +gamma22(i,k)q2(i,k-1) )
341
342	d_tb_dadv(i,k) = ( rr11(i)(th1(i,k)-thw(i,k))+rr12(i)(th2(i,k)-thx(i,k)) )*ppi(i,k)
343	d_deltat_dadv(i,k) = ( rr21(i)(th1(i,k)-thw(i,k))+rr22(i)(th2(i,k)-thx(i,k)) )*ppi(i,k)
344	d_qb_dadv(i,k) = rr11(i)(q1(i,k)-qw(i,k)) +rr12(i)(q2(i,k)-qx(i,k))
345	d_deltaq_dadv(i,k) = rr21(i)(q1(i,k)-qw(i,k)) +rr22(i)(q2(i,k)-qx(i,k))
346	ENDIF ! (wk_adv(i) .AND. k<=kupper(i))
347	ENDDO
348	ENDDO
349
350	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
351	!!!!!!!!!!!!!!!!!! Verification de l'inversion !!!!!!!!!!!!!!!!!!!!!!!
352	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
353	!!
354	!! DO i = 1,klon
355	!! xt1inv(i,1) = aa11(i,1)th1(i,1) + aa12(i,1)th2(i,1) + bb11(i,1)th1(i,2) + bb12(i,1)th2(i,2)
356	!! xt2inv(i,1) = aa21(i,1)th1(i,1) + aa22(i,1)th2(i,1) + bb21(i,1)th1(i,2) + bb22(i,1)th2(i,2)
357	!! xq1inv(i,1) = aa11(i,1)q1(i,1) + aa12(i,1)q2(i,1) + bb11(i,1)q1(i,2) + bb12(i,1)q2(i,2)
358	!! xq2inv(i,1) = aa21(i,1)q1(i,1) + aa22(i,1)q2(i,1) + bb21(i,1)q1(i,2) + bb22(i,1)q2(i,2)
359	!! ENDDO
360	!!
361	!! DO k = 2, klev-1
362	!! DO i = 1,klon
363	!! xt1inv(i,k) = aa11(i,k)th1(i,k) + aa12(i,k)th2(i,k) + bb11(i,k)th1(i,k+1) + bb12(i,k)th2(i,k+1) &
364	!! + cc11(i,k)th1(i,k-1) + cc12(i,k)th2(i,k-1)
365	!! xt2inv(i,k) = aa21(i,k)th1(i,k) + aa22(i,k)th2(i,k) + bb21(i,k)th1(i,k+1) + bb22(i,k)th2(i,k+1) &
366	!! + cc21(i,k)th1(i,k-1) + cc22(i,k)th2(i,k-1)
367	!! xq1inv(i,k) = aa11(i,k)q1(i,k) + aa12(i,k)q2(i,k) + bb11(i,k)q1(i,k+1) + bb12(i,k)q2(i,k+1) &
368	!! + cc11(i,k)q1(i,k-1) + cc12(i,k)q2(i,k-1)
369	!! xq2inv(i,k) = aa21(i,k)q1(i,k) + aa22(i,k)q2(i,k) + bb21(i,k)q1(i,k+1) + bb22(i,k)q2(i,k+1) &
370	!! + cc21(i,k)q1(i,k-1) + cc22(i,k)q2(i,k-1)
371	!! ENDDO
372	!! ENDDO
373	!!
374	!! DO i = 1,klon
375	!! xt1inv(i,klev) = aa11(i,klev)th1(i,klev) + aa12(i,klev)th2(i,klev) + cc11(i,klev)th1(i,klev-1) + cc12(i,klev)th2(i,klev-1)
376	!! xt2inv(i,klev) = aa21(i,klev)th1(i,klev) + aa22(i,klev)th2(i,klev) + cc21(i,klev)th1(i,klev-1) + cc22(i,klev)th2(i,klev-1)
377	!! xq1inv(i,klev) = aa11(i,klev)q1(i,klev) + aa12(i,klev)q2(i,klev) + cc11(i,klev)q1(i,klev-1) + cc12(i,klev)q2(i,klev-1)
378	!! xq2inv(i,klev) = aa21(i,klev)q1(i,klev) + aa22(i,klev)q2(i,klev) + cc21(i,klev)q1(i,klev-1) + cc22(i,klev)q2(i,klev-1)
379	!! ENDDO
380	!!
381	!! DO k = 1, 20
382	!! IF (abs(xt1inv(1,k)-xt1(1,k)) .GT. 1.e-15*xt1(1,k) ) THEN
383	!! print *,'wake_dadv, k, xt1inv(1,k), xt1(1,k), xt1inv(1,k)-xt1(1,k) ', &
384	!! k, xt1inv(1,k), xt1(1,k), xt1inv(1,k)-xt1(1,k)
385	!! ENDIF
386	!! IF (abs(xt2inv(1,k)-xt2(1,k)) .GT. 1.e-15*xt2(1,k) ) THEN
387	!! print *,'wake_dadv, k, xt2inv(1,k), xt2(1,k), xt2inv(1,k)-xt2(1,k) ', &
388	!! k, xt2inv(1,k), xt2(1,k), xt2inv(1,k)-xt2(1,k)
389	!! ENDIF
390	!! IF (abs(xq1inv(1,k)-xq1(1,k)) .GT. 1.e-15*xq1(1,k) ) THEN
391	!! print *,'wake_dadv, k, xq1inv(1,k), xq1(1,k), xq1inv(1,k)-xq1(1,k) ', &
392	!! k, xq1inv(1,k), xq1(1,k), xq1inv(1,k)-xq1(1,k)
393	!! ENDIF
394	!! IF (abs(xq2inv(1,k)-xq2(1,k)) .GT. 1.e-15*xq2(1,k) ) THEN
395	!! print *,'wake_dadv, k, xq2inv(1,k), xq2(1,k), xq2inv(1,k)-xq2(1,k) ', &
396	!! k, xq2inv(1,k), xq2(1,k), xq2inv(1,k)-xq2(1,k)
397	!! ENDIF
398	!! ENDDO
399
400	ELSE ! (flag_dadv_implicit)
401
402	! Explicit scheme : compute directly d_deltat_dadv and d_tb_dadv
403	! (and similarly for d_deltaq_dadv and d_qb_dadv).
404	! -----------------------------------------------------------------------------------------
405
406	DO i = 1, klon
407	IF (wk_adv(i)) THEN !!! nrlmd
408	d_thx(i, 1) = 0.
409	d_thw(i, 1) = 0.
410	d_dth(i, 1) = 0.
411	d_qx(i, 1) = 0.
412	d_qw(i, 1) = 0.
413	d_dq(i, 1) = 0.
414	END IF
415	END DO
416
417	DO k = 2, klev
418	DO i = 1, klon
419	IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN
420	d_thx(i, k) = thx(i, k-1) - thx(i, k)
421	d_thw(i, k) = thw(i, k-1) - thw(i, k)
422	d_dth(i, k) = delta_th(i, k-1) - delta_th(i, k)
423	d_qx(i, k) = qx(i, k-1) - qx(i, k)
424	d_qw(i, k) = qw(i, k-1) - qw(i, k)
425	d_dq(i, k) = delta_q(i, k-1) - delta_q(i, k)
426	END IF ! (wk_adv(i) .AND. k<=kupper(i)+1)
427	END DO
428	END DO
429
430	DO k = 1, klev-1
431	DO i = 1, klon
432	IF (wk_adv(i) .AND. k<=kupper(i)-1) THEN
433	d_deltat_dadv(i, k) = dtime/(ph(i,k)-ph(i,k+1))* &
434	(rr22(i)deltomg(i,k)sigmaw(i)*d_thx(i,k) - &
435	rr21(i)deltomg(i,k+1)(1.-sigmaw(i))d_thw(i,k+1) )ppi(i, k) - &
436	dtimeentr_s(i,k)delta_th(i,k)/sigmaw(i)*ppi(i, k)
437	!
438	d_deltaq_dadv(i, k) = dtime/(ph(i,k)-ph(i,k+1))* &
439	(rr22(i)deltomg(i,k)sigmaw(i)*d_qx(i,k)- &
440	rr21(i)deltomg(i,k+1)(1.-sigmaw(i))*d_qw(i,k+1) ) - &
441	dtimeentr_s(i,k)delta_q(i,k)/sigmaw(i)
442
443	! and increment large scale tendencies
444	d_tb_dadv(i, k) = dtime((rr12(i)deltomg(i,k)sigmaw(i)d_thx(i,k)- &
445	rr11(i)deltomg(i,k+1)(1.-sigmaw(i))*d_thw(i,k+1))/ &
446	(ph(i,k)-ph(i,k+1)) &
447	-sigmaw(i)(1.-sigmaw(i))delta_th(i,k)*(deltomg(i,k)-deltomg(i,k+1))/ &
448	(ph(i,k)-ph(i,k+1)) )*ppi(i, k)
449
450	d_qb_dadv(i, k) = dtime((rr12(i)deltomg(i,k)sigmaw(i)d_qx(i,k)- &
451	rr11(i)deltomg(i,k+1)(1.-sigmaw(i))*d_qw(i,k+1))/ &
452	(ph(i,k)-ph(i,k+1)) &
453	-sigmaw(i)(1.-sigmaw(i))delta_q(i,k)*(deltomg(i,k)-deltomg(i,k+1))/ &
454	(ph(i,k)-ph(i,k+1)) )
455	ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN
456	d_tb_dadv(i, k) = dtime(rr12(i)deltomg(i, k)sigmaw(i)d_thx(i, k)/(ph(i, k)-ph(i, k+1)))*ppi(i, k)
457
458	d_qb_dadv(i, k) = dtime(rr12(i)deltomg(i, k)sigmaw(i)d_qx(i, k)/(ph(i, k)-ph(i, k+1)))
459	END IF ! (wk_adv(i) .AND. k<=kupper(i)-1)
460	END DO
461	END DO
462
463	ENDIF! (flag_dadv_implicit)
464
465	!print *,'ZZZwake_dadv k, d_deltat_dadv(1,k) ', (k, d_deltat_dadv(1,k), k = 1,3)
466
467	END SUBROUTINE wake_dadv
468	END MODULE lmdz_wake_dadv

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

Download in other formats:

Original Format