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

Last change on this file since 1081 was 1067, checked in by Laurent Fairhead, 16 years ago
  • Modifications lie au premier niveau du modele pour la diffusion turbulent

du vent.

  • Preparation pour un couplage des courrant oceaniques.

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