1 |
|
|
MODULE climb_qbs_mod |
2 |
|
|
! |
3 |
|
|
! Module to solve the verctical diffusion of blowing snow; |
4 |
|
|
! |
5 |
|
|
USE dimphy |
6 |
|
|
|
7 |
|
|
IMPLICIT NONE |
8 |
|
|
SAVE |
9 |
|
|
PRIVATE |
10 |
|
|
PUBLIC :: climb_qbs_down, climb_qbs_up |
11 |
|
|
|
12 |
|
|
REAL, DIMENSION(:,:), ALLOCATABLE :: gamaqbs |
13 |
|
|
!$OMP THREADPRIVATE(gamaqbs) |
14 |
|
|
REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_QBS, Dcoef_QBS |
15 |
|
|
!$OMP THREADPRIVATE(Ccoef_QBS, Dcoef_QBS) |
16 |
|
|
REAL, DIMENSION(:), ALLOCATABLE :: Acoef_QBS, Bcoef_QBS |
17 |
|
|
!$OMP THREADPRIVATE(Acoef_QBS, Bcoef_QBS) |
18 |
|
|
REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefqbs |
19 |
|
|
!$OMP THREADPRIVATE(Kcoefqbs) |
20 |
|
|
|
21 |
|
|
CONTAINS |
22 |
|
|
! |
23 |
|
|
!**************************************************************************************** |
24 |
|
|
! |
25 |
|
|
SUBROUTINE climb_qbs_down(knon, coefqbs, paprs, pplay, & |
26 |
|
|
delp, temp, qbs, dtime, & |
27 |
|
|
Ccoef_QBS_out, Dcoef_QBS_out, & |
28 |
|
|
Kcoef_qbs_out, gama_qbs_out, & |
29 |
|
|
Acoef_QBS_out, Bcoef_QBS_out) |
30 |
|
|
|
31 |
|
|
! This routine calculates recursivly the coefficients C and D |
32 |
|
|
! for the quantity X=[QBS] in equation X(k) = C(k) + D(k)*X(k-1), where k is |
33 |
|
|
! the index of the vertical layer. |
34 |
|
|
! |
35 |
|
|
! Input arguments |
36 |
|
|
!**************************************************************************************** |
37 |
|
|
INTEGER, INTENT(IN) :: knon |
38 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: coefqbs |
39 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay |
40 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs |
41 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: delp |
42 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: temp |
43 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: qbs |
44 |
|
|
REAL, INTENT(IN) :: dtime |
45 |
|
|
|
46 |
|
|
! Output arguments |
47 |
|
|
!**************************************************************************************** |
48 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_QBS_out |
49 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_QBS_out |
50 |
|
|
|
51 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_QBS_out |
52 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_QBS_out |
53 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: Kcoef_qbs_out |
54 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: gama_qbs_out |
55 |
|
|
|
56 |
|
|
! Local variables |
57 |
|
|
!**************************************************************************************** |
58 |
|
|
LOGICAL, SAVE :: first=.TRUE. |
59 |
|
|
!$OMP THREADPRIVATE(first) |
60 |
|
|
REAL, DIMENSION(klon) :: psref |
61 |
|
|
REAL :: delz, pkh |
62 |
|
|
INTEGER :: k, i, ierr |
63 |
|
|
! Include |
64 |
|
|
!**************************************************************************************** |
65 |
|
|
INCLUDE "YOMCST.h" |
66 |
|
|
INCLUDE "compbl.h" |
67 |
|
|
|
68 |
|
|
|
69 |
|
|
!**************************************************************************************** |
70 |
|
|
! 1) |
71 |
|
|
! Allocation at first time step only |
72 |
|
|
! |
73 |
|
|
!**************************************************************************************** |
74 |
|
|
|
75 |
|
|
IF (first) THEN |
76 |
|
|
first=.FALSE. |
77 |
|
|
ALLOCATE(Ccoef_QBS(klon,klev), STAT=ierr) |
78 |
|
|
IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_QBS, ierr=', ierr |
79 |
|
|
|
80 |
|
|
ALLOCATE(Dcoef_QBS(klon,klev), STAT=ierr) |
81 |
|
|
IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_QBS, ierr=', ierr |
82 |
|
|
|
83 |
|
|
ALLOCATE(Acoef_QBS(klon), Bcoef_QBS(klon), STAT=ierr) |
84 |
|
|
IF ( ierr /= 0 ) PRINT*,' pb in allloc Acoef_BS and Bcoef_BS, ierr=', ierr |
85 |
|
|
|
86 |
|
|
ALLOCATE(Kcoefqbs(klon,klev), STAT=ierr) |
87 |
|
|
IF ( ierr /= 0 ) PRINT*,' pb in allloc Kcoefqbs, ierr=', ierr |
88 |
|
|
|
89 |
|
|
ALLOCATE(gamaqbs(1:klon,2:klev), STAT=ierr) |
90 |
|
|
IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaqbs, ierr=', ierr |
91 |
|
|
|
92 |
|
|
END IF |
93 |
|
|
|
94 |
|
|
!**************************************************************************************** |
95 |
|
|
! 2) |
96 |
|
|
! Definition of the coeficient K |
97 |
|
|
! |
98 |
|
|
!**************************************************************************************** |
99 |
|
|
Kcoefqbs(:,:) = 0.0 |
100 |
|
|
DO k = 2, klev |
101 |
|
|
DO i = 1, knon |
102 |
|
|
Kcoefqbs(i,k) = & |
103 |
|
|
coefqbs(i,k)*RG*RG*dtime /(pplay(i,k-1)-pplay(i,k)) & |
104 |
|
|
*(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2 |
105 |
|
|
ENDDO |
106 |
|
|
ENDDO |
107 |
|
|
|
108 |
|
|
!**************************************************************************************** |
109 |
|
|
! 3) |
110 |
|
|
! Calculation of gama for "Q" and "H" |
111 |
|
|
! |
112 |
|
|
!**************************************************************************************** |
113 |
|
|
! surface pressure is used as reference |
114 |
|
|
psref(:) = paprs(:,1) |
115 |
|
|
|
116 |
|
|
! definition of gama |
117 |
|
|
IF (iflag_pbl == 1) THEN |
118 |
|
|
gamaqbs(:,:) = 0.0 |
119 |
|
|
|
120 |
|
|
! conversion de gama |
121 |
|
|
DO k = 2, klev |
122 |
|
|
DO i = 1, knon |
123 |
|
|
delz = RD * (temp(i,k-1)+temp(i,k)) / & |
124 |
|
|
2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k)) |
125 |
|
|
pkh = (psref(i)/paprs(i,k))**RKAPPA |
126 |
|
|
|
127 |
|
|
! convertie gradient verticale de contenu en neige soufflee en difference de neige soufflee entre centre de couches |
128 |
|
|
gamaqbs(i,k) = gamaqbs(i,k) * delz |
129 |
|
|
ENDDO |
130 |
|
|
ENDDO |
131 |
|
|
|
132 |
|
|
ELSE |
133 |
|
|
gamaqbs(:,:) = 0.0 |
134 |
|
|
ENDIF |
135 |
|
|
|
136 |
|
|
|
137 |
|
|
!**************************************************************************************** |
138 |
|
|
! 4) |
139 |
|
|
! Calculte the coefficients C and D for specific content of blowing snow, qbs |
140 |
|
|
! |
141 |
|
|
!**************************************************************************************** |
142 |
|
|
|
143 |
|
|
CALL calc_coef_qbs(knon, Kcoefqbs(:,:), gamaqbs(:,:), delp(:,:), qbs(:,:), & |
144 |
|
|
Ccoef_QBS(:,:), Dcoef_QBS(:,:), Acoef_QBS, Bcoef_QBS) |
145 |
|
|
|
146 |
|
|
|
147 |
|
|
!**************************************************************************************** |
148 |
|
|
! 5) |
149 |
|
|
! Return the first layer in output variables |
150 |
|
|
! |
151 |
|
|
!**************************************************************************************** |
152 |
|
|
Acoef_QBS_out = Acoef_QBS |
153 |
|
|
Bcoef_QBS_out = Bcoef_QBS |
154 |
|
|
|
155 |
|
|
!**************************************************************************************** |
156 |
|
|
! 6) |
157 |
|
|
! If Pbl is split, return also the other layers in output variables |
158 |
|
|
! |
159 |
|
|
!**************************************************************************************** |
160 |
|
|
IF (mod(iflag_pbl_split,10) .ge.1) THEN |
161 |
|
|
DO k= 1, klev |
162 |
|
|
DO i= 1, klon |
163 |
|
|
Ccoef_QBS_out(i,k) = Ccoef_QBS(i,k) |
164 |
|
|
Dcoef_QBS_out(i,k) = Dcoef_QBS(i,k) |
165 |
|
|
Kcoef_qbs_out(i,k) = Kcoefqbs(i,k) |
166 |
|
|
IF (k.eq.1) THEN |
167 |
|
|
gama_qbs_out(i,k) = 0. |
168 |
|
|
ELSE |
169 |
|
|
gama_qbs_out(i,k) = gamaqbs(i,k) |
170 |
|
|
ENDIF |
171 |
|
|
ENDDO |
172 |
|
|
ENDDO |
173 |
|
|
!!! |
174 |
|
|
ENDIF ! (mod(iflag_pbl_split,2) .ge.1) |
175 |
|
|
!!! |
176 |
|
|
|
177 |
|
|
END SUBROUTINE climb_qbs_down |
178 |
|
|
! |
179 |
|
|
!**************************************************************************************** |
180 |
|
|
! |
181 |
|
|
SUBROUTINE calc_coef_qbs(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef) |
182 |
|
|
! |
183 |
|
|
! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1) |
184 |
|
|
! where X is QQBS, and k the vertical level k=1,klev |
185 |
|
|
! |
186 |
|
|
INCLUDE "YOMCST.h" |
187 |
|
|
! Input arguments |
188 |
|
|
!**************************************************************************************** |
189 |
|
|
INTEGER, INTENT(IN) :: knon |
190 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef, delp |
191 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: X |
192 |
|
|
REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama |
193 |
|
|
|
194 |
|
|
! Output arguments |
195 |
|
|
!**************************************************************************************** |
196 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: Acoef, Bcoef |
197 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef, Dcoef |
198 |
|
|
|
199 |
|
|
! Local variables |
200 |
|
|
!**************************************************************************************** |
201 |
|
|
INTEGER :: k, i |
202 |
|
|
REAL :: buf |
203 |
|
|
|
204 |
|
|
!**************************************************************************************** |
205 |
|
|
! Niveau au sommet, k=klev |
206 |
|
|
! |
207 |
|
|
!**************************************************************************************** |
208 |
|
|
Ccoef(:,:) = 0.0 |
209 |
|
|
Dcoef(:,:) = 0.0 |
210 |
|
|
|
211 |
|
|
DO i = 1, knon |
212 |
|
|
buf = delp(i,klev) + Kcoef(i,klev) |
213 |
|
|
|
214 |
|
|
Ccoef(i,klev) = (X(i,klev)*delp(i,klev) - Kcoef(i,klev)*gama(i,klev))/buf |
215 |
|
|
Dcoef(i,klev) = Kcoef(i,klev)/buf |
216 |
|
|
END DO |
217 |
|
|
|
218 |
|
|
|
219 |
|
|
!**************************************************************************************** |
220 |
|
|
! Niveau (klev-1) <= k <= 2 |
221 |
|
|
! |
222 |
|
|
!**************************************************************************************** |
223 |
|
|
|
224 |
|
|
DO k=(klev-1),2,-1 |
225 |
|
|
DO i = 1, knon |
226 |
|
|
buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1)) |
227 |
|
|
Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1) + & |
228 |
|
|
Kcoef(i,k+1)*gama(i,k+1) - Kcoef(i,k)*gama(i,k))/buf |
229 |
|
|
Dcoef(i,k) = Kcoef(i,k)/buf |
230 |
|
|
END DO |
231 |
|
|
END DO |
232 |
|
|
|
233 |
|
|
!**************************************************************************************** |
234 |
|
|
! Niveau k=1 |
235 |
|
|
! |
236 |
|
|
!**************************************************************************************** |
237 |
|
|
|
238 |
|
|
DO i = 1, knon |
239 |
|
|
buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2)) |
240 |
|
|
Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*(gama(i,2)+Ccoef(i,2)))/buf |
241 |
|
|
Bcoef(i) = -1. * RG / buf |
242 |
|
|
END DO |
243 |
|
|
|
244 |
|
|
END SUBROUTINE calc_coef_qbs |
245 |
|
|
! |
246 |
|
|
!**************************************************************************************** |
247 |
|
|
! |
248 |
|
|
SUBROUTINE climb_qbs_up(knon, dtime, qbs_old, & |
249 |
|
|
flx_qbs1, paprs, pplay, & |
250 |
|
|
Acoef_QBS_in, Bcoef_QBS_in, & |
251 |
|
|
Ccoef_QBS_in, Dcoef_QBS_in, & |
252 |
|
|
Kcoef_qbs_in, gama_qbs_in, & |
253 |
|
|
flux_qbs, d_qbs) |
254 |
|
|
! |
255 |
|
|
! This routine calculates the flux and tendency of the specific content of blowing snow qbs |
256 |
|
|
! The quantity qbs is calculated according to |
257 |
|
|
! X(k) = C(k) + D(k)*X(k-1) for X=[qbs], where the coefficients |
258 |
|
|
! C and D are known from before and k is index of the vertical layer. |
259 |
|
|
! |
260 |
|
|
|
261 |
|
|
! Input arguments |
262 |
|
|
!**************************************************************************************** |
263 |
|
|
INTEGER, INTENT(IN) :: knon |
264 |
|
|
REAL, INTENT(IN) :: dtime |
265 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: qbs_old |
266 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: flx_qbs1 |
267 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs |
268 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay |
269 |
|
|
|
270 |
|
|
!!! nrlmd le 02/05/2011 |
271 |
|
|
REAL, DIMENSION(klon), INTENT(IN) :: Acoef_QBS_in, Bcoef_QBS_in |
272 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: Ccoef_QBS_in, Dcoef_QBS_in |
273 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef_qbs_in, gama_qbs_in |
274 |
|
|
!!! |
275 |
|
|
|
276 |
|
|
! Output arguments |
277 |
|
|
!**************************************************************************************** |
278 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: flux_qbs, d_qbs |
279 |
|
|
|
280 |
|
|
! Local variables |
281 |
|
|
!**************************************************************************************** |
282 |
|
|
LOGICAL, SAVE :: last=.FALSE. |
283 |
|
|
REAL, DIMENSION(klon,klev) :: qbs_new |
284 |
|
|
REAL, DIMENSION(klon) :: psref |
285 |
|
|
INTEGER :: k, i, ierr |
286 |
|
|
|
287 |
|
|
! Include |
288 |
|
|
!**************************************************************************************** |
289 |
|
|
INCLUDE "YOMCST.h" |
290 |
|
|
INCLUDE "compbl.h" |
291 |
|
|
|
292 |
|
|
!**************************************************************************************** |
293 |
|
|
! 1) |
294 |
|
|
! Definition of some variables |
295 |
|
|
REAL, DIMENSION(klon,klev) :: zairm |
296 |
|
|
! |
297 |
|
|
!**************************************************************************************** |
298 |
|
|
flux_qbs(:,:) = 0.0 |
299 |
|
|
d_qbs(:,:) = 0.0 |
300 |
|
|
|
301 |
|
|
psref(1:knon) = paprs(1:knon,1) |
302 |
|
|
|
303 |
|
|
IF (mod(iflag_pbl_split,10) .ge.1) THEN |
304 |
|
|
DO i = 1, knon |
305 |
|
|
Acoef_QBS(i)=Acoef_QBS_in(i) |
306 |
|
|
Bcoef_QBS(i)=Bcoef_QBS_in(i) |
307 |
|
|
ENDDO |
308 |
|
|
DO k = 1, klev |
309 |
|
|
DO i = 1, knon |
310 |
|
|
Ccoef_QBS(i,k)=Ccoef_QBS_in(i,k) |
311 |
|
|
Dcoef_QBS(i,k)=Dcoef_QBS_in(i,k) |
312 |
|
|
Kcoefqbs(i,k)=Kcoef_qbs_in(i,k) |
313 |
|
|
IF (k.gt.1) THEN |
314 |
|
|
gamaqbs(i,k)=gama_qbs_in(i,k) |
315 |
|
|
ENDIF |
316 |
|
|
ENDDO |
317 |
|
|
ENDDO |
318 |
|
|
!!! |
319 |
|
|
ENDIF ! (mod(iflag_pbl_split,2) .ge.1) |
320 |
|
|
!!! |
321 |
|
|
|
322 |
|
|
!**************************************************************************************** |
323 |
|
|
! 2) |
324 |
|
|
! Calculation of QBS |
325 |
|
|
! |
326 |
|
|
!**************************************************************************************** |
327 |
|
|
|
328 |
|
|
!- First layer |
329 |
|
|
qbs_new(1:knon,1) = Acoef_QBS(1:knon) + Bcoef_QBS(1:knon)*flx_qbs1(1:knon)*dtime |
330 |
|
|
!- All the other layers |
331 |
|
|
DO k = 2, klev |
332 |
|
|
DO i = 1, knon |
333 |
|
|
qbs_new(i,k) = Ccoef_QBS(i,k) + Dcoef_QBS(i,k)*qbs_new(i,k-1) |
334 |
|
|
END DO |
335 |
|
|
END DO |
336 |
|
|
!**************************************************************************************** |
337 |
|
|
! 3) |
338 |
|
|
! Calculation of the flux for QBS |
339 |
|
|
! |
340 |
|
|
!**************************************************************************************** |
341 |
|
|
|
342 |
|
|
!- The flux at first layer, k=1 |
343 |
|
|
flux_qbs(1:knon,1)=flx_qbs1(1:knon) |
344 |
|
|
|
345 |
|
|
!- The flux at all layers above surface |
346 |
|
|
DO k = 2, klev |
347 |
|
|
DO i = 1, knon |
348 |
|
|
flux_qbs(i,k) = (Kcoefqbs(i,k)/RG/dtime) * & |
349 |
|
|
(qbs_new(i,k)-qbs_new(i,k-1)+gamaqbs(i,k)) |
350 |
|
|
END DO |
351 |
|
|
END DO |
352 |
|
|
|
353 |
|
|
!**************************************************************************************** |
354 |
|
|
! 4) |
355 |
|
|
! Calculation of tendency for QBS |
356 |
|
|
! |
357 |
|
|
!**************************************************************************************** |
358 |
|
|
DO k = 1, klev |
359 |
|
|
DO i = 1, knon |
360 |
|
|
d_qbs(i,k) = qbs_new(i,k) - qbs_old(i,k) |
361 |
|
|
zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg |
362 |
|
|
END DO |
363 |
|
|
END DO |
364 |
|
|
|
365 |
|
|
!**************************************************************************************** |
366 |
|
|
! Some deallocations |
367 |
|
|
! |
368 |
|
|
!**************************************************************************************** |
369 |
|
|
IF (last) THEN |
370 |
|
|
DEALLOCATE(Ccoef_QBS, Dcoef_QBS,stat=ierr) |
371 |
|
|
IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Ccoef_QBS, Dcoef_QBS, ierr=', ierr |
372 |
|
|
DEALLOCATE(Acoef_QBS, Bcoef_QBS,stat=ierr) |
373 |
|
|
IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Acoef_QBS, Bcoef_QBS, ierr=', ierr |
374 |
|
|
DEALLOCATE(gamaqbs,stat=ierr) |
375 |
|
|
IF ( ierr /= 0 ) PRINT*,' pb in dealllocate gamaqbs, ierr=', ierr |
376 |
|
|
DEALLOCATE(Kcoefqbs,stat=ierr) |
377 |
|
|
IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Kcoefqbs, ierr=', ierr |
378 |
|
|
END IF |
379 |
|
|
END SUBROUTINE climb_qbs_up |
380 |
|
|
! |
381 |
|
|
!**************************************************************************************** |
382 |
|
|
! |
383 |
|
|
END MODULE climb_qbs_mod |
384 |
|
|
|
385 |
|
|
|
386 |
|
|
|
387 |
|
|
|
388 |
|
|
|
389 |
|
|
|