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

Last change on this file since 1408 was 1084, checked in by jghattas, 16 years ago

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