source: LMDZ6/trunk/libf/phylmd/climb_qbs_mod.f90 @ 5296

Last change on this file since 5296 was 5296, checked in by abarral, 29 hours ago

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