source: LMDZ6/trunk/libf/phylmd/climb_qbs_mod.F90 @ 4771

Last change on this file since 4771 was 4523, checked in by evignon, 18 months ago

merge de la branche blowing snow vers la trunk
premiere tentative
Etienne
File size: 13.6 KB
RevLine 
[4485]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.