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

Last change on this file since 996 was 996, checked in by lsce, 16 years ago
  • Modifications liées au calcul des nouveau sous-fractions
  • Nettoyage de ocean slab : il reste uniquement la version avec glace de mer forcé
  • Nouveaux variables pour distiguer la version et type d'ocean : type_ocean=force/slab/couple, version_ocean=opa8/nemo pour couplé ou version_ocean=sicOBS pour slab

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