source: LMDZ4/trunk/libf/phylmd/climb_hq_mod.F90 @ 970

Last change on this file since 970 was 793, checked in by Laurent Fairhead, 17 years ago

Modifications suite a la transformation des fichiers include pour
qu'ils soient compatibles a la fois au format fixe et au format libre
Un bon nombre de fichiers *.inc du coup disparaissent
LF
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 12.5 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 :: Kcoefhq
20  !$OMP THREADPRIVATE(Kcoefhq)
21
22CONTAINS
23!
24!****************************************************************************************
25!
26  SUBROUTINE climb_hq_down(knon, coefhq, paprs, pplay, &
27       delp, temp, q, dtime, &
28       petAcoef, peqAcoef, petBcoef, peqBcoef)
29
30    INCLUDE "YOMCST.h"
31! This routine calculates recursivly the coefficients C and D
32! for the quantity X=[Q,H] 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)   :: coefhq
39    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
40    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
41    REAL, DIMENSION(klon,klev), INTENT(IN)   :: temp, delp  ! temperature
42    REAL, DIMENSION(klon,klev), INTENT(IN)   :: q
43    REAL, INTENT(IN)                         :: dtime
44
45! Output arguments
46!****************************************************************************************
47    REAL, DIMENSION(klon), INTENT(OUT)       :: petAcoef
48    REAL, DIMENSION(klon), INTENT(OUT)       :: peqAcoef
49    REAL, DIMENSION(klon), INTENT(OUT)       :: petBcoef
50    REAL, DIMENSION(klon), INTENT(OUT)       :: peqBcoef
51
52! Local variables
53!****************************************************************************************
54    REAL, DIMENSION(klon,klev)               :: dels, local_H
55    REAL, DIMENSION(klon)                    :: psref, delz, pkh
56    INTEGER                                  :: k, i, ierr
57
58! Include
59!****************************************************************************************
60  INCLUDE "compbl.h"   
61
62
63!****************************************************************************************
64! 1)
65! Allocation
66!   
67!****************************************************************************************
68
69    ALLOCATE(Ccoef_Q(klon,klev), STAT=ierr)
70    IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_Q, ierr=', ierr
71   
72    ALLOCATE(Dcoef_Q(klon,klev), STAT=ierr)
73    IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_Q, ierr=', ierr
74   
75    ALLOCATE(Ccoef_H(klon,klev), STAT=ierr)
76    IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_H, ierr=', ierr
77   
78    ALLOCATE(Dcoef_H(klon,klev), STAT=ierr)
79    IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_H, ierr=', ierr
80   
81    ALLOCATE(Kcoefhq(klon,klev), STAT=ierr)
82    IF ( ierr /= 0 )  PRINT*,' pb in allloc Kcoefhq, ierr=', ierr
83   
84    ALLOCATE(gamaq(1:klon,2:klev), STAT=ierr)
85    IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaq, ierr=', ierr
86   
87    ALLOCATE(gamah(1:klon,2:klev), STAT=ierr)
88    IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr
89       
90
91!****************************************************************************************
92! 2)
93! Definition of the coeficient K
94!
95!****************************************************************************************
96    Kcoefhq(:,:) = 0.0
97    DO k = 2, klev
98       DO i = 1, knon
99          Kcoefhq(i,k) = &
100               coefhq(i,k)*RG /(pplay(i,k-1)-pplay(i,k)) &
101               *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
102          Kcoefhq(i,k) = Kcoefhq(i,k) * dtime*RG
103       ENDDO
104    ENDDO
105
106!****************************************************************************************
107! 3)
108! Calculation of gama for "Q" and "H"
109!
110!****************************************************************************************
111!   surface pressure is used as reference
112    psref(:) = paprs(:,1)
113
114!   definition of start value for gama
115    IF (iflag_pbl == 1) THEN
116       gamaq(:,:) = 0.0
117       gamah(:,:) = -1.0e-03
118       gamah(:,2) = -2.5e-03
119    ELSE
120       gamaq(:,:) = 0.0
121       gamah(:,:) = 0.0
122    ENDIF
123   
124!   calculation of gama
125    DO k = 2, klev
126       DO i = 1, knon
127          delz(i)  = RD * (temp(i,k-1)+temp(i,k)) / &
128               2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
129          pkh(i) = (psref(i)/paprs(i,k))**RKAPPA
130         
131          gamaq(i,k) = gamaq(i,k) * delz(i)   
132          gamah(i,k) = gamah(i,k) * delz(i) * RCPD * pkh(i)
133       ENDDO
134    ENDDO
135
136!****************************************************************************************   
137! 4)
138! Calculte the coefficients C and D for specific humidity, q
139!
140!****************************************************************************************
141    dels(:,:) = delp(:,:)
142   
143    CALL calc_coef(knon, Kcoefhq(:,:), gamaq(:,:), dels(:,:), q(:,:), &
144         Ccoef_Q(:,:), Dcoef_Q(:,:))
145
146!****************************************************************************************
147! 5)
148! Calculte the coefficients C and D for potentiel entalpie, H
149!
150!****************************************************************************************
151    dels(:,:)    = 0.0
152    local_H(:,:) = 0.0
153
154    DO k=1,klev
155       DO i = 1, knon
156          dels(i,k) = (pplay(i,k)/psref(i))**RKAPPA * delp(i,k)               
157          local_H(i,k) = RCPD * temp(i,k) * &
158               (psref(i)/pplay(i,k))**RKAPPA
159       ENDDO
160    ENDDO
161
162
163    CALL calc_coef(knon, Kcoefhq(:,:), gamah(:,:), dels(:,:), local_H(:,:), &
164         Ccoef_H(:,:), Dcoef_H(:,:))
165 
166!****************************************************************************************
167! 6)
168! Return the first layer in output variables
169!
170!****************************************************************************************
171    petAcoef(1:knon) = Ccoef_H(1:knon,1)
172    peqAcoef(1:knon) = Ccoef_Q(1:knon,1)
173    petBcoef(1:knon) = Dcoef_H(1:knon,1)
174    peqBcoef(1:knon) = Dcoef_Q(1:knon,1)
175
176  END SUBROUTINE climb_hq_down
177!
178!****************************************************************************************
179!
180  SUBROUTINE calc_coef(knon, Kcoef, gama, dels, X, Ccoef, Dcoef)
181!
182! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
183! where X is H or Q, and k the vertical level k=1,klev
184!
185    INCLUDE "YOMCST.h"
186! Input arguments
187!****************************************************************************************
188    INTEGER, INTENT(IN)                      :: knon
189    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, dels
190    REAL, DIMENSION(klon,klev), INTENT(IN)   :: X
191    REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama
192
193! Output arguments
194!****************************************************************************************
195    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef
196    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef
197
198! Local variables
199!****************************************************************************************
200    INTEGER                                  :: k, i
201    REAL                                     :: buf
202
203!****************************************************************************************
204! Niveau au sommet, k=klev
205!
206!****************************************************************************************
207    Ccoef(:,:) = 0.0
208    Dcoef(:,:) = 0.0
209
210    DO i = 1, knon
211       buf = dels(i,klev) + Kcoef(i,klev)
212       
213       Ccoef(i,klev) = (X(i,klev)*dels(i,klev) - Kcoef(i,klev)*gama(i,klev))/buf
214       Dcoef(i,klev) = Kcoef(i,klev)/buf
215    END DO
216
217
218!****************************************************************************************
219! Niveau  (klev-1) <= k <= 2
220!
221!****************************************************************************************
222
223    DO k=(klev-1),2,-1
224       DO i = 1, knon
225          buf = dels(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
226          Ccoef(i,k) = (X(i,k)*dels(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1) + &
227               Kcoef(i,k+1)*gama(i,k+1) - Kcoef(i,k)*gama(i,k))/buf
228          Dcoef(i,k) = Kcoef(i,k)/buf
229       END DO
230    END DO
231
232!****************************************************************************************
233! Niveau k=1
234!
235!****************************************************************************************
236
237    DO i = 1, knon
238       buf = dels(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
239       Ccoef(i,1) = (X(i,1)*dels(i,1) + Kcoef(i,2)*(gama(i,2)+Ccoef(i,2)))/buf
240       Dcoef(i,1) = -1. * RG / buf
241    END DO
242
243  END SUBROUTINE calc_coef
244!
245!****************************************************************************************
246!
247  SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, &
248       flx_q1, flx_t1, paprs, pplay, &
249       flux_q, flux_h, d_q, d_t)
250!
251! This routine calculates the flux and tendency of the specific humidity q and
252! the potential engergi H.
253! The quantities q and H are calculated according to
254! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
255! C and D are known from before and k is index of the vertical layer.
256!   
257    INCLUDE "YOMCST.h"
258! Input arguments
259!****************************************************************************************
260    INTEGER, INTENT(IN)                      :: knon
261    REAL, INTENT(IN)                         :: dtime
262    REAL, DIMENSION(klon,klev), INTENT(IN)   :: t_old, q_old
263    REAL, DIMENSION(klon), INTENT(IN)        :: flx_q1, flx_t1
264    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
265    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
266
267! Output arguments
268!****************************************************************************************
269    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: flux_q, flux_h, d_q, d_t
270
271! Local variables
272!****************************************************************************************
273    REAL, DIMENSION(klon,klev)               :: zx_pkh, zx_pkf
274    REAL, DIMENSION(klon,klev)               :: h_new, q_new
275    REAL, DIMENSION(klon)                    :: psref         
276    INTEGER                                  :: k, i
277
278!****************************************************************************************
279! 1)
280! Definition of some variables
281!
282!****************************************************************************************
283    flux_q(:,:) = 0.0
284    flux_h(:,:) = 0.0
285    d_q(:,:)    = 0.0
286    d_t(:,:)    = 0.0
287
288    psref(1:knon) = paprs(1:knon,1) 
289
290    DO k = 1, klev
291       DO i = 1, knon
292          zx_pkh(i,k) = (psref(i)/paprs(i,k))**RKAPPA
293          zx_pkf(i,k) = (psref(i)/pplay(i,k))**RKAPPA
294       END DO
295    END DO
296!****************************************************************************************
297! 2)
298! Calculation of Q and H
299!
300!****************************************************************************************
301
302!- First layer
303    q_new(1:knon,1) = Ccoef_Q(1:knon,1) + Dcoef_Q(1:knon,1)*flx_q1(1:knon)*dtime
304    h_new(1:knon,1) = Ccoef_H(1:knon,1) + Dcoef_H(1:knon,1)*flx_t1(1:knon)*dtime
305   
306!- All the rest layers
307    DO k = 2, klev
308       DO i = 1, knon
309          q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1)
310          h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1)
311       END DO
312    END DO
313!****************************************************************************************
314! 3)
315! Calculation of the flux for Q and H
316!
317!****************************************************************************************
318
319!- The flux at first layer, k=1
320    flux_q(1:knon,1)=flx_q1(1:knon)
321    flux_h(1:knon,1)=flx_t1(1:knon)
322
323!- The flux at all layers above surface
324    DO k = 2, klev
325       DO i = 1, knon
326          flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
327               (q_new(i,k)-q_new(i,k-1)+gamaq(i,k))
328
329          flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
330               (h_new(i,k)-h_new(i,k-1)+gamah(i,k)) / &
331               zx_pkh(i,k)
332       END DO
333    END DO
334
335!****************************************************************************************
336! 4)
337! Calculation of tendency for Q and H
338!
339!****************************************************************************************
340
341    DO k = 1, klev
342       DO i = 1, knon
343          d_t(i,k) = h_new(i,k)/zx_pkf(i,k)/RCPD - t_old(i,k)
344          d_q(i,k) = q_new(i,k) - q_old(i,k)
345       END DO
346    END DO
347
348!****************************************************************************************
349! Some deallocations
350!
351!****************************************************************************************
352    DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H)   
353    DEALLOCATE(gamaq, gamah)
354    DEALLOCATE(Kcoefhq)
355
356   
357  END SUBROUTINE climb_hq_up
358!
359!****************************************************************************************
360!
361END MODULE climb_hq_mod
362
363
364
365
366
Note: See TracBrowser for help on using the repository browser.