source: LMDZ5/branches/IPSLCM5A2.1/libf/phylmd/climb_hq_mod.F90 @ 5421

Last change on this file since 5421 was 2159, checked in by jyg, 10 years ago

1/ Splitting of the boundary layer : the climbing down and up of Pbl_surface is
split between the off-wake and wake regions ; the thermal scheme is applied
only to the off-wake region.
2/ Elimination of wake_scal and calwake_scal.

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 16.3 KB
Line 
1MODULE climb_hq_mod
2!
3! Module to solve the verctical diffusion of "q" and "H";
4! specific humidity and potential energi.
5!
6  USE dimphy
7
8  IMPLICIT NONE
9  SAVE
10  PRIVATE
11  PUBLIC :: climb_hq_down, climb_hq_up
12
13  REAL, DIMENSION(:,:), ALLOCATABLE :: gamaq, gamah
14  !$OMP THREADPRIVATE(gamaq,gamah)
15  REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_Q, Dcoef_Q
16  !$OMP THREADPRIVATE(Ccoef_Q, Dcoef_Q)
17  REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_H, Dcoef_H
18  !$OMP THREADPRIVATE(Ccoef_H, Dcoef_H)
19  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_Q, Bcoef_Q
20  !$OMP THREADPRIVATE(Acoef_Q, Bcoef_Q)
21  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_H, Bcoef_H
22  !$OMP THREADPRIVATE(Acoef_H, Bcoef_H)
23  REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefhq
24  !$OMP THREADPRIVATE(Kcoefhq)
25
26CONTAINS
27!
28!****************************************************************************************
29!
30  SUBROUTINE climb_hq_down(knon, coefhq, paprs, pplay, &
31       delp, temp, q, dtime, &
32!!! nrlmd le 02/05/2011
33       Ccoef_H_out, Ccoef_Q_out, Dcoef_H_out, Dcoef_Q_out, &
34       Kcoef_hq_out, gama_q_out, gama_h_out, &
35!!!
36       Acoef_H_out, Acoef_Q_out, Bcoef_H_out, Bcoef_Q_out)
37
38! This routine calculates recursivly the coefficients C and D
39! for the quantity X=[Q,H] in equation X(k) = C(k) + D(k)*X(k-1), where k is
40! the index of the vertical layer.
41!
42! Input arguments
43!****************************************************************************************
44    INTEGER, INTENT(IN)                      :: knon
45    REAL, DIMENSION(klon,klev), INTENT(IN)   :: coefhq
46    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
47    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
48    REAL, DIMENSION(klon,klev), INTENT(IN)   :: temp, delp  ! temperature
49    REAL, DIMENSION(klon,klev), INTENT(IN)   :: q
50    REAL, INTENT(IN)                         :: dtime
51
52! Output arguments
53!****************************************************************************************
54    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_H_out
55    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_Q_out
56    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_H_out
57    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_Q_out
58
59!!! nrlmd le 02/05/2011
60    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_H_out
61    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_Q_out
62    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_H_out
63    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_Q_out
64    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Kcoef_hq_out
65    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: gama_q_out
66    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: gama_h_out
67!!!
68
69! Local variables
70!****************************************************************************************
71    LOGICAL, SAVE                            :: first=.TRUE.
72    !$OMP THREADPRIVATE(first)
73    REAL, DIMENSION(klon,klev)               :: local_H
74    REAL, DIMENSION(klon)                    :: psref
75    REAL                                     :: delz, pkh
76    INTEGER                                  :: k, i, ierr
77
78! Include
79!****************************************************************************************
80    INCLUDE "YOMCST.h"
81    INCLUDE "compbl.h"   
82
83
84!****************************************************************************************
85! 1)
86! Allocation at first time step only
87!   
88!****************************************************************************************
89
90    IF (first) THEN
91       first=.FALSE.
92       ALLOCATE(Ccoef_Q(klon,klev), STAT=ierr)
93       IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_Q, ierr=', ierr
94       
95       ALLOCATE(Dcoef_Q(klon,klev), STAT=ierr)
96       IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_Q, ierr=', ierr
97       
98       ALLOCATE(Ccoef_H(klon,klev), STAT=ierr)
99       IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_H, ierr=', ierr
100       
101       ALLOCATE(Dcoef_H(klon,klev), STAT=ierr)
102       IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_H, ierr=', ierr
103       
104       ALLOCATE(Acoef_Q(klon), Bcoef_Q(klon), Acoef_H(klon), Bcoef_H(klon), STAT=ierr)
105       IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_X and Bcoef_X, ierr=', ierr
106       
107       ALLOCATE(Kcoefhq(klon,klev), STAT=ierr)
108       IF ( ierr /= 0 )  PRINT*,' pb in allloc Kcoefhq, ierr=', ierr
109       
110       ALLOCATE(gamaq(1:klon,2:klev), STAT=ierr)
111       IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaq, ierr=', ierr
112       
113       ALLOCATE(gamah(1:klon,2:klev), STAT=ierr)
114       IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr
115    END IF
116
117!****************************************************************************************
118! 2)
119! Definition of the coeficient K
120!
121!****************************************************************************************
122    Kcoefhq(:,:) = 0.0
123    DO k = 2, klev
124       DO i = 1, knon
125          Kcoefhq(i,k) = &
126               coefhq(i,k)*RG*RG*dtime /(pplay(i,k-1)-pplay(i,k)) &
127               *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
128       ENDDO
129    ENDDO
130
131!****************************************************************************************
132! 3)
133! Calculation of gama for "Q" and "H"
134!
135!****************************************************************************************
136!   surface pressure is used as reference
137    psref(:) = paprs(:,1)
138
139!   definition of gama
140    IF (iflag_pbl == 1) THEN
141       gamaq(:,:) = 0.0
142       gamah(:,:) = -1.0e-03
143       gamah(:,2) = -2.5e-03
144 
145! conversion de gama
146       DO k = 2, klev
147          DO i = 1, knon
148             delz = RD * (temp(i,k-1)+temp(i,k)) / &
149                    2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
150             pkh  = (psref(i)/paprs(i,k))**RKAPPA
151         
152! convertie gradient verticale d'humidite specifique en difference d'humidite specifique entre centre de couches
153             gamaq(i,k) = gamaq(i,k) * delz   
154! convertie gradient verticale de temperature en difference de temperature potentielle entre centre de couches
155             gamah(i,k) = gamah(i,k) * delz * RCPD * pkh
156          ENDDO
157       ENDDO
158
159    ELSE
160       gamaq(:,:) = 0.0
161       gamah(:,:) = 0.0
162    ENDIF
163   
164
165!****************************************************************************************   
166! 4)
167! Calculte the coefficients C and D for specific humidity, q
168!
169!****************************************************************************************
170   
171    CALL calc_coef(knon, Kcoefhq(:,:), gamaq(:,:), delp(:,:), q(:,:), &
172         Ccoef_Q(:,:), Dcoef_Q(:,:), Acoef_Q, Bcoef_Q)
173
174!****************************************************************************************
175! 5)
176! Calculte the coefficients C and D for potentiel entalpie, H
177!
178!****************************************************************************************
179    local_H(:,:) = 0.0
180
181    DO k=1,klev
182       DO i = 1, knon
183          ! convertie la temperature en entalpie potentielle
184          local_H(i,k) = RCPD * temp(i,k) * &
185               (psref(i)/pplay(i,k))**RKAPPA
186       ENDDO
187    ENDDO
188
189    CALL calc_coef(knon, Kcoefhq(:,:), gamah(:,:), delp(:,:), local_H(:,:), &
190         Ccoef_H(:,:), Dcoef_H(:,:), Acoef_H, Bcoef_H)
191 
192!****************************************************************************************
193! 6)
194! Return the first layer in output variables
195!
196!****************************************************************************************
197    Acoef_H_out = Acoef_H
198    Bcoef_H_out = Bcoef_H
199    Acoef_Q_out = Acoef_Q
200    Bcoef_Q_out = Bcoef_Q
201
202!****************************************************************************************
203! 7)
204! If Pbl is split, return also the other layers in output variables
205!
206!****************************************************************************************
207!!! jyg le 07/02/2012
208       IF (mod(iflag_pbl_split,2) .eq.1) THEN
209!!! nrlmd le 02/05/2011
210    DO k= 1, klev
211      DO i= 1, klon
212        Ccoef_H_out(i,k) = Ccoef_H(i,k)
213        Dcoef_H_out(i,k) = Dcoef_H(i,k)
214        Ccoef_Q_out(i,k) = Ccoef_Q(i,k)
215        Dcoef_Q_out(i,k) = Dcoef_Q(i,k)
216        Kcoef_hq_out(i,k) = Kcoefhq(i,k)
217          IF (k.eq.1) THEN
218            gama_h_out(i,k)  = 0.
219            gama_q_out(i,k)  = 0.
220          ELSE
221            gama_h_out(i,k)  = gamah(i,k)
222            gama_q_out(i,k)  = gamaq(i,k)
223          ENDIF
224      ENDDO
225    ENDDO
226!!!     
227       ENDIF  ! (mod(iflag_pbl_split,2) .eq.1)
228!!!
229
230  END SUBROUTINE climb_hq_down
231!
232!****************************************************************************************
233!
234  SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
235!
236! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
237! where X is H or Q, and k the vertical level k=1,klev
238!
239    INCLUDE "YOMCST.h"
240! Input arguments
241!****************************************************************************************
242    INTEGER, INTENT(IN)                      :: knon
243    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, delp
244    REAL, DIMENSION(klon,klev), INTENT(IN)   :: X
245    REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama
246
247! Output arguments
248!****************************************************************************************
249    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef, Bcoef
250    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef, Dcoef
251
252! Local variables
253!****************************************************************************************
254    INTEGER                                  :: k, i
255    REAL                                     :: buf
256
257!****************************************************************************************
258! Niveau au sommet, k=klev
259!
260!****************************************************************************************
261    Ccoef(:,:) = 0.0
262    Dcoef(:,:) = 0.0
263
264    DO i = 1, knon
265       buf = delp(i,klev) + Kcoef(i,klev)
266       
267       Ccoef(i,klev) = (X(i,klev)*delp(i,klev) - Kcoef(i,klev)*gama(i,klev))/buf
268       Dcoef(i,klev) = Kcoef(i,klev)/buf
269    END DO
270
271
272!****************************************************************************************
273! Niveau  (klev-1) <= k <= 2
274!
275!****************************************************************************************
276
277    DO k=(klev-1),2,-1
278       DO i = 1, knon
279          buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
280          Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1) + &
281               Kcoef(i,k+1)*gama(i,k+1) - Kcoef(i,k)*gama(i,k))/buf
282          Dcoef(i,k) = Kcoef(i,k)/buf
283       END DO
284    END DO
285
286!****************************************************************************************
287! Niveau k=1
288!
289!****************************************************************************************
290
291    DO i = 1, knon
292       buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
293       Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*(gama(i,2)+Ccoef(i,2)))/buf
294       Bcoef(i) = -1. * RG / buf
295    END DO
296
297  END SUBROUTINE calc_coef
298!
299!****************************************************************************************
300!
301  SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, &
302       flx_q1, flx_h1, paprs, pplay, &
303!!! nrlmd le 02/05/2011
304       Acoef_H_in, Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in, &
305       Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in, &
306       Kcoef_hq_in, gama_q_in, gama_h_in, &
307!!!
308       flux_q, flux_h, d_q, d_t)
309!
310! This routine calculates the flux and tendency of the specific humidity q and
311! the potential engergi H.
312! The quantities q and H are calculated according to
313! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
314! C and D are known from before and k is index of the vertical layer.
315!   
316
317! Input arguments
318!****************************************************************************************
319    INTEGER, INTENT(IN)                      :: knon
320    REAL, INTENT(IN)                         :: dtime
321    REAL, DIMENSION(klon,klev), INTENT(IN)   :: t_old, q_old
322    REAL, DIMENSION(klon), INTENT(IN)        :: flx_q1, flx_h1
323    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
324    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
325
326!!! nrlmd le 02/05/2011
327    REAL, DIMENSION(klon), INTENT(IN)        :: Acoef_H_in,Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in
328    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in
329    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef_hq_in, gama_q_in, gama_h_in
330!!!
331
332! Output arguments
333!****************************************************************************************
334    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: flux_q, flux_h, d_q, d_t
335
336! Local variables
337!****************************************************************************************
338    LOGICAL, SAVE                            :: last=.FALSE.
339    REAL, DIMENSION(klon,klev)               :: h_new, q_new
340    REAL, DIMENSION(klon)                    :: psref         
341    INTEGER                                  :: k, i, ierr
342 
343! Include
344!****************************************************************************************
345    INCLUDE "YOMCST.h"
346    INCLUDE "compbl.h"   
347
348!****************************************************************************************
349! 1)
350! Definition of some variables
351!
352!****************************************************************************************
353    flux_q(:,:) = 0.0
354    flux_h(:,:) = 0.0
355    d_q(:,:)    = 0.0
356    d_t(:,:)    = 0.0
357
358    psref(1:knon) = paprs(1:knon,1) 
359
360!!! jyg le 07/02/2012
361       IF (mod(iflag_pbl_split,2) .eq.1) THEN
362!!! nrlmd le 02/05/2011
363    DO i = 1, knon
364      Acoef_H(i)=Acoef_H_in(i)
365      Acoef_Q(i)=Acoef_Q_in(i)
366      Bcoef_H(i)=Bcoef_H_in(i)
367      Bcoef_Q(i)=Bcoef_Q_in(i)
368    ENDDO
369    DO k = 1, klev
370      DO i = 1, knon
371        Ccoef_H(i,k)=Ccoef_H_in(i,k)
372        Ccoef_Q(i,k)=Ccoef_Q_in(i,k)
373        Dcoef_H(i,k)=Dcoef_H_in(i,k)
374        Dcoef_Q(i,k)=Dcoef_Q_in(i,k)
375        Kcoefhq(i,k)=Kcoef_hq_in(i,k)
376          IF (k.gt.1) THEN
377            gamah(i,k)=gama_h_in(i,k)
378            gamaq(i,k)=gama_q_in(i,k)
379          ENDIF
380      ENDDO
381    ENDDO
382!!!     
383       ENDIF  ! (mod(iflag_pbl_split,2) .eq.1)
384!!!
385
386!****************************************************************************************
387! 2)
388! Calculation of Q and H
389!
390!****************************************************************************************
391
392!- First layer
393    q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)*flx_q1(1:knon)*dtime
394    h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)*flx_h1(1:knon)*dtime
395   
396!- All the other layers
397    DO k = 2, klev
398       DO i = 1, knon
399          q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1)
400          h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1)
401       END DO
402    END DO
403!****************************************************************************************
404! 3)
405! Calculation of the flux for Q and H
406!
407!****************************************************************************************
408
409!- The flux at first layer, k=1
410    flux_q(1:knon,1)=flx_q1(1:knon)
411    flux_h(1:knon,1)=flx_h1(1:knon)
412
413!- The flux at all layers above surface
414    DO k = 2, klev
415       DO i = 1, knon
416          flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
417               (q_new(i,k)-q_new(i,k-1)+gamaq(i,k))
418
419          flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
420               (h_new(i,k)-h_new(i,k-1)+gamah(i,k))
421       END DO
422    END DO
423
424!****************************************************************************************
425! 4)
426! Calculation of tendency for Q and H
427!
428!****************************************************************************************
429
430    DO k = 1, klev
431       DO i = 1, knon
432          d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
433          d_q(i,k) = q_new(i,k) - q_old(i,k)
434       END DO
435    END DO
436
437!****************************************************************************************
438! Some deallocations
439!
440!****************************************************************************************
441    IF (last) THEN
442       DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)   
443       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
444       DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)   
445       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
446       DEALLOCATE(gamaq, gamah,stat=ierr)
447       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
448       DEALLOCATE(Kcoefhq,stat=ierr)
449       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
450    END IF
451  END SUBROUTINE climb_hq_up
452!
453!****************************************************************************************
454!
455END MODULE climb_hq_mod
456
457 
458
459
460
461
Note: See TracBrowser for help on using the repository browser.