source: lmdz_wrf/WRFV3/lmdz/climb_hq_mod.F90 @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

  • Property svn:executable set to *
File size: 14.0 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    PRINT *,'    Lluis in climb_hq_down 1 k dplay dpaprs Kcoefhq gamaq gamah '// &
189      'delp q Ccoef_Q Dcoef_Q Acoef_Q Bcoef_Q local_H ________'
190
191    DO k=2,klev
192      PRINT *,k,pplay(1,k-1)-pplay(1,k),paprs(1,k-1)-paprs(1,k),Kcoefhq(1,k), &
193        gamaq(1,k),gamah(1,k),delp(1,k),q(1,k),Ccoef_Q(1,k),Dcoef_Q(1,k),     &
194        Acoef_Q(1),Bcoef_Q(1), local_H(1,k)
195    END DO
196
197  END SUBROUTINE climb_hq_down
198!
199!****************************************************************************************
200!
201  SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
202!
203! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
204! where X is H or Q, and k the vertical level k=1,klev
205!
206    INCLUDE "YOMCST.h"
207! Input arguments
208!****************************************************************************************
209    INTEGER, INTENT(IN)                      :: knon
210    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, delp
211    REAL, DIMENSION(klon,klev), INTENT(IN)   :: X
212    REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama
213
214! Output arguments
215!****************************************************************************************
216    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef, Bcoef
217    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef, Dcoef
218
219! Local variables
220!****************************************************************************************
221    INTEGER                                  :: k, i
222    REAL                                     :: buf
223
224!****************************************************************************************
225! Niveau au sommet, k=klev
226!
227!****************************************************************************************
228    Ccoef(:,:) = 0.0
229    Dcoef(:,:) = 0.0
230
231    DO i = 1, knon
232       buf = delp(i,klev) + Kcoef(i,klev)
233       
234       Ccoef(i,klev) = (X(i,klev)*delp(i,klev) - Kcoef(i,klev)*gama(i,klev))/buf
235       Dcoef(i,klev) = Kcoef(i,klev)/buf
236    END DO
237
238
239!****************************************************************************************
240! Niveau  (klev-1) <= k <= 2
241!
242!****************************************************************************************
243
244    DO k=(klev-1),2,-1
245       DO i = 1, knon
246          buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
247          Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1) + &
248               Kcoef(i,k+1)*gama(i,k+1) - Kcoef(i,k)*gama(i,k))/buf
249          Dcoef(i,k) = Kcoef(i,k)/buf
250       END DO
251    END DO
252
253!****************************************************************************************
254! Niveau k=1
255!
256!****************************************************************************************
257
258    DO i = 1, knon
259       buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
260       Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*(gama(i,2)+Ccoef(i,2)))/buf
261       Bcoef(i) = -1. * RG / buf
262    END DO
263    acoef(knon+1: klon) = 0.
264    bcoef(knon+1: klon) = 0.
265
266  END SUBROUTINE calc_coef
267!
268!****************************************************************************************
269!
270  SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, &
271       flx_q1, flx_h1, paprs, pplay, &
272       flux_q, flux_h, d_q, d_t)
273!
274! This routine calculates the flux and tendency of the specific humidity q and
275! the potential engergi H.
276! The quantities q and H are calculated according to
277! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
278! C and D are known from before and k is index of the vertical layer.
279!   
280    INCLUDE "YOMCST.h"
281! Input arguments
282!****************************************************************************************
283    INTEGER, INTENT(IN)                      :: knon
284    REAL, INTENT(IN)                         :: dtime
285    REAL, DIMENSION(klon,klev), INTENT(IN)   :: t_old, q_old
286    REAL, DIMENSION(klon), INTENT(IN)        :: flx_q1, flx_h1
287    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
288    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
289
290! Output arguments
291!****************************************************************************************
292    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: flux_q, flux_h, d_q, d_t
293
294! Local variables
295!****************************************************************************************
296    LOGICAL, SAVE                            :: last=.FALSE.
297    REAL, DIMENSION(klon,klev)               :: h_new, q_new
298    REAL, DIMENSION(klon)                    :: psref         
299    INTEGER                                  :: k, i, ierr
300
301!****************************************************************************************
302! 1)
303! Definition of some variables
304!
305!****************************************************************************************
306    flux_q(:,:) = 0.0
307    flux_h(:,:) = 0.0
308    d_q(:,:)    = 0.0
309    d_t(:,:)    = 0.0
310
311    psref(1:knon) = paprs(1:knon,1) 
312
313!****************************************************************************************
314! 2)
315! Calculation of Q and H
316!
317!****************************************************************************************
318
319!- First layer
320    q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)*flx_q1(1:knon)*dtime
321    h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)*flx_h1(1:knon)*dtime
322   
323!- All the other layers
324    DO k = 2, klev
325       DO i = 1, knon
326          q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1)
327          h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1)
328       END DO
329    END DO
330!****************************************************************************************
331! 3)
332! Calculation of the flux for Q and H
333!
334!****************************************************************************************
335
336!- The flux at first layer, k=1
337    flux_q(1:knon,1)=flx_q1(1:knon)
338    flux_h(1:knon,1)=flx_h1(1:knon)
339
340!- The flux at all layers above surface
341    DO k = 2, klev
342       DO i = 1, knon
343          flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
344               (q_new(i,k)-q_new(i,k-1)+gamaq(i,k))
345
346          flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
347               (h_new(i,k)-h_new(i,k-1)+gamah(i,k))
348       END DO
349    END DO
350
351!****************************************************************************************
352! 4)
353! Calculation of tendency for Q and H
354!
355!****************************************************************************************
356
357    DO k = 1, klev
358       DO i = 1, knon
359          d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
360          d_q(i,k) = q_new(i,k) - q_old(i,k)
361       END DO
362    END DO
363
364!****************************************************************************************
365! Some deallocations
366!
367!****************************************************************************************
368    IF (last) THEN
369       DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)   
370       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
371       DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)   
372       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
373       DEALLOCATE(gamaq, gamah,stat=ierr)
374       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
375       DEALLOCATE(Kcoefhq,stat=ierr)
376       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
377    END IF
378  END SUBROUTINE climb_hq_up
379!
380!****************************************************************************************
381!
382END MODULE climb_hq_mod
383
384 
385
386
387
388
Note: See TracBrowser for help on using the repository browser.