My Project
 All Classes Files Functions Variables Macros
climb_hq_mod.F90
Go to the documentation of this file.
1 MODULE 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 
26 CONTAINS
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 !
373 END MODULE climb_hq_mod
374 
375 
376 
377 
378 
379