source: LMDZ5/branches/LMDZ5-DOFOCO/libf/phylmd/climb_hq_mod.F90 @ 2265

Last change on this file since 2265 was 1618, checked in by lguez, 13 years ago

In module "climb_wind_mod", in procedure "calc_coef", "bcoef" is
declared with dimension "klon" but was only defined up to index
"knon". In procedure "climb_wind_down", "bcoef_v" is defined from
argument "bcoef" of "calc_coef". The whole array "bcoef_v" is
used. This produced a run-time error with the XLF compiler and
debugging options because undefined elements of "bcoef_v" contained
NaN. Fixed by filling the useless part of "bcoef" with zeros. Same fix
for "acoef". Same fix in "climb_hq_mod".
  • 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    acoef(knon+1: klon) = 0.
255    bcoef(knon+1: klon) = 0.
256
257  END SUBROUTINE calc_coef
258!
259!****************************************************************************************
260!
261  SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, &
262       flx_q1, flx_h1, paprs, pplay, &
263       flux_q, flux_h, d_q, d_t)
264!
265! This routine calculates the flux and tendency of the specific humidity q and
266! the potential engergi H.
267! The quantities q and H are calculated according to
268! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
269! C and D are known from before and k is index of the vertical layer.
270!   
271    INCLUDE "YOMCST.h"
272! Input arguments
273!****************************************************************************************
274    INTEGER, INTENT(IN)                      :: knon
275    REAL, INTENT(IN)                         :: dtime
276    REAL, DIMENSION(klon,klev), INTENT(IN)   :: t_old, q_old
277    REAL, DIMENSION(klon), INTENT(IN)        :: flx_q1, flx_h1
278    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
279    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
280
281! Output arguments
282!****************************************************************************************
283    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: flux_q, flux_h, d_q, d_t
284
285! Local variables
286!****************************************************************************************
287    LOGICAL, SAVE                            :: last=.FALSE.
288    REAL, DIMENSION(klon,klev)               :: h_new, q_new
289    REAL, DIMENSION(klon)                    :: psref         
290    INTEGER                                  :: k, i, ierr
291
292!****************************************************************************************
293! 1)
294! Definition of some variables
295!
296!****************************************************************************************
297    flux_q(:,:) = 0.0
298    flux_h(:,:) = 0.0
299    d_q(:,:)    = 0.0
300    d_t(:,:)    = 0.0
301
302    psref(1:knon) = paprs(1:knon,1) 
303
304!****************************************************************************************
305! 2)
306! Calculation of Q and H
307!
308!****************************************************************************************
309
310!- First layer
311    q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)*flx_q1(1:knon)*dtime
312    h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)*flx_h1(1:knon)*dtime
313   
314!- All the other layers
315    DO k = 2, klev
316       DO i = 1, knon
317          q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1)
318          h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1)
319       END DO
320    END DO
321!****************************************************************************************
322! 3)
323! Calculation of the flux for Q and H
324!
325!****************************************************************************************
326
327!- The flux at first layer, k=1
328    flux_q(1:knon,1)=flx_q1(1:knon)
329    flux_h(1:knon,1)=flx_h1(1:knon)
330
331!- The flux at all layers above surface
332    DO k = 2, klev
333       DO i = 1, knon
334          flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
335               (q_new(i,k)-q_new(i,k-1)+gamaq(i,k))
336
337          flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
338               (h_new(i,k)-h_new(i,k-1)+gamah(i,k))
339       END DO
340    END DO
341
342!****************************************************************************************
343! 4)
344! Calculation of tendency for Q and H
345!
346!****************************************************************************************
347
348    DO k = 1, klev
349       DO i = 1, knon
350          d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
351          d_q(i,k) = q_new(i,k) - q_old(i,k)
352       END DO
353    END DO
354
355!****************************************************************************************
356! Some deallocations
357!
358!****************************************************************************************
359    IF (last) THEN
360       DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)   
361       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
362       DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)   
363       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
364       DEALLOCATE(gamaq, gamah,stat=ierr)
365       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
366       DEALLOCATE(Kcoefhq,stat=ierr)
367       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
368    END IF
369  END SUBROUTINE climb_hq_up
370!
371!****************************************************************************************
372!
373END MODULE climb_hq_mod
374
375 
376
377
378
379
Note: See TracBrowser for help on using the repository browser.