source: LMDZ6/branches/Amaury_dev/libf/phylmd/climb_qbs_mod.F90 @ 5099

Last change on this file since 5099 was 5099, checked in by abarral, 2 months ago

Replace most uses of CPP_DUST by the corresponding logical defined in lmdz_cppkeys_wrapper.F90
Convert several files from .F to .f90 to allow Dust to compile w/o rrtm/ecrad
Create lmdz_yoerad.f90
(lint) Remove "!" on otherwise empty line
File size: 13.5 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) >=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==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! Input arguments
261!****************************************************************************************
262    INTEGER, INTENT(IN)                      :: knon
263    REAL, INTENT(IN)                         :: dtime
264    REAL, DIMENSION(klon,klev), INTENT(IN)   :: qbs_old
265    REAL, DIMENSION(klon), INTENT(IN)        :: flx_qbs1
266    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
267    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
268
269!!! nrlmd le 02/05/2011
270    REAL, DIMENSION(klon), INTENT(IN)        :: Acoef_QBS_in, Bcoef_QBS_in
271    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Ccoef_QBS_in, Dcoef_QBS_in
272    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef_qbs_in, gama_qbs_in
273!!!
274
275! Output arguments
276!****************************************************************************************
277    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: flux_qbs, d_qbs
278
279! Local variables
280!****************************************************************************************
281    LOGICAL, SAVE                            :: last=.FALSE.
282    REAL, DIMENSION(klon,klev)               :: qbs_new
283    REAL, DIMENSION(klon)                    :: psref         
284    INTEGER                                  :: k, i, ierr
285 
286! Include
287!****************************************************************************************
288    INCLUDE "YOMCST.h"
289    INCLUDE "compbl.h"   
290
291!****************************************************************************************
292! 1)
293! Definition of some variables
294    REAL, DIMENSION(klon,klev)               :: zairm
295
296!****************************************************************************************
297    flux_qbs(:,:) = 0.0
298    d_qbs(:,:)    = 0.0
299
300    psref(1:knon) = paprs(1:knon,1) 
301
302       IF (mod(iflag_pbl_split,10) >=1) THEN
303    DO i = 1, knon
304      Acoef_QBS(i)=Acoef_QBS_in(i)
305      Bcoef_QBS(i)=Bcoef_QBS_in(i)
306    ENDDO
307    DO k = 1, klev
308      DO i = 1, knon
309        Ccoef_QBS(i,k)=Ccoef_QBS_in(i,k)
310        Dcoef_QBS(i,k)=Dcoef_QBS_in(i,k)
311        Kcoefqbs(i,k)=Kcoef_qbs_in(i,k)
312          IF (k>1) THEN
313            gamaqbs(i,k)=gama_qbs_in(i,k)
314          ENDIF
315      ENDDO
316    ENDDO
317!!!     
318       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
319!!!
320
321!****************************************************************************************
322! 2)
323! Calculation of QBS
324
325!****************************************************************************************
326
327!- First layer
328    qbs_new(1:knon,1) = Acoef_QBS(1:knon) + Bcoef_QBS(1:knon)*flx_qbs1(1:knon)*dtime
329!- All the other layers
330    DO k = 2, klev
331       DO i = 1, knon
332          qbs_new(i,k) = Ccoef_QBS(i,k) + Dcoef_QBS(i,k)*qbs_new(i,k-1)
333       END DO
334    END DO
335!****************************************************************************************
336! 3)
337! Calculation of the flux for QBS
338
339!****************************************************************************************
340
341!- The flux at first layer, k=1
342    flux_qbs(1:knon,1)=flx_qbs1(1:knon)
343
344!- The flux at all layers above surface
345    DO k = 2, klev
346       DO i = 1, knon
347          flux_qbs(i,k) = (Kcoefqbs(i,k)/RG/dtime) * &
348               (qbs_new(i,k)-qbs_new(i,k-1)+gamaqbs(i,k))
349       END DO
350    END DO
351
352!****************************************************************************************
353! 4)
354! Calculation of tendency for QBS
355
356!****************************************************************************************
357    DO k = 1, klev
358       DO i = 1, knon
359          d_qbs(i,k) = qbs_new(i,k) - qbs_old(i,k)
360          zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
361        END DO
362    END DO
363
364!****************************************************************************************
365! Some deallocations
366
367!****************************************************************************************
368    IF (last) THEN
369       DEALLOCATE(Ccoef_QBS, Dcoef_QBS,stat=ierr)   
370       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_QBS, Dcoef_QBS, ierr=', ierr
371       DEALLOCATE(Acoef_QBS, Bcoef_QBS,stat=ierr)   
372       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Acoef_QBS, Bcoef_QBS, ierr=', ierr
373       DEALLOCATE(gamaqbs,stat=ierr)
374       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaqbs, ierr=', ierr
375       DEALLOCATE(Kcoefqbs,stat=ierr)
376       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefqbs, ierr=', ierr
377    END IF
378  END SUBROUTINE climb_qbs_up
379
380!****************************************************************************************
381
382END MODULE climb_qbs_mod
383
384 
385
386
387
388
Note: See TracBrowser for help on using the repository browser.