source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/climb_qbs_mod.F90 @ 4771

Last change on this file since 4771 was 4727, checked in by idelkadi, 9 months ago

Merged trunk changes -r4488:4726 LMDZ_ECRad branch
File size: 13.6 KB
Line 
1MODULE 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
21CONTAINS
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!
383END MODULE climb_qbs_mod
384
385 
386
387
388
389
Note: See TracBrowser for help on using the repository browser.