source: LMDZ6/branches/Amaury_dev/libf/phylmd/climb_hq_mod.F90 @ 5139

Last change on this file since 5139 was 5139, checked in by abarral, 8 weeks ago

Put nuage.h, flux_arp.h, compbl.h into modules
Move unused phylmd/ini_hist* to obsolete
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 17.8 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  PRIVATE
10  PUBLIC :: climb_hq_down, climb_hq_up, d_h_col_vdf, f_h_bnd
11
12  REAL, DIMENSION(:,:), ALLOCATABLE :: gamaq, gamah
13  !$OMP THREADPRIVATE(gamaq,gamah)
14  REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_Q, Dcoef_Q
15  !$OMP THREADPRIVATE(Ccoef_Q, Dcoef_Q)
16  REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_H, Dcoef_H
17  !$OMP THREADPRIVATE(Ccoef_H, Dcoef_H)
18  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_Q, Bcoef_Q
19  !$OMP THREADPRIVATE(Acoef_Q, Bcoef_Q)
20  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_H, Bcoef_H
21  !$OMP THREADPRIVATE(Acoef_H, Bcoef_H)
22  REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefhq
23  !$OMP THREADPRIVATE(Kcoefhq)
24  REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: h_old ! for diagnostics, h before solving diffusion
25  !$OMP THREADPRIVATE(h_old)
26  REAL, SAVE, DIMENSION(:), ALLOCATABLE :: d_h_col_vdf ! for diagnostics, vertical integral of enthalpy change
27  !$OMP THREADPRIVATE(d_h_col_vdf)
28  REAL, SAVE, DIMENSION(:), ALLOCATABLE :: f_h_bnd ! for diagnostics, enthalpy flux at surface
29  !$OMP THREADPRIVATE(f_h_bnd)
30
31CONTAINS
32
33!****************************************************************************************
34
35  SUBROUTINE climb_hq_down(knon, coefhq, paprs, pplay, &
36       delp, temp, q, dtime, &
37!!! nrlmd le 02/05/2011
38       Ccoef_H_out, Ccoef_Q_out, Dcoef_H_out, Dcoef_Q_out, &
39       Kcoef_hq_out, gama_q_out, gama_h_out, &
40!!!
41       Acoef_H_out, Acoef_Q_out, Bcoef_H_out, Bcoef_Q_out)
42
43! This routine calculates recursivly the coefficients C and D
44! for the quantity X=[Q,H] in equation X(k) = C(k) + D(k)*X(k-1), where k is
45! the index of the vertical layer.
46
47! Input arguments
48!****************************************************************************************
49    USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
50
51    IMPLICIT NONE
52    INTEGER, INTENT(IN)                      :: knon
53    REAL, DIMENSION(klon,klev), INTENT(IN)   :: coefhq
54    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
55    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
56    REAL, DIMENSION(klon,klev), INTENT(IN)   :: temp, delp  ! temperature
57    REAL, DIMENSION(klon,klev), INTENT(IN)   :: q
58    REAL, INTENT(IN)                         :: dtime
59
60! Output arguments
61!****************************************************************************************
62    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_H_out
63    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_Q_out
64    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_H_out
65    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_Q_out
66
67!!! nrlmd le 02/05/2011
68    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_H_out
69    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_Q_out
70    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_H_out
71    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_Q_out
72    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Kcoef_hq_out
73    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: gama_q_out
74    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: gama_h_out
75!!!
76
77! Local variables
78!****************************************************************************************
79    LOGICAL, SAVE                            :: first=.TRUE.
80    !$OMP THREADPRIVATE(first)
81! JLD now renamed h_old and declared in module
82!    REAL, DIMENSION(klon,klev)               :: local_H
83    REAL, DIMENSION(klon)                    :: psref
84    REAL                                     :: delz, pkh
85    INTEGER                                  :: k, i, ierr
86! Include
87!****************************************************************************************
88    INCLUDE "YOMCST.h"
89!****************************************************************************************
90! 1)
91! Allocation at first time step only
92
93!****************************************************************************************
94
95    IF (first) THEN
96       first=.FALSE.
97       ALLOCATE(Ccoef_Q(klon,klev), STAT=ierr)
98       IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_Q, ierr=', ierr
99       
100       ALLOCATE(Dcoef_Q(klon,klev), STAT=ierr)
101       IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_Q, ierr=', ierr
102       
103       ALLOCATE(Ccoef_H(klon,klev), STAT=ierr)
104       IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_H, ierr=', ierr
105       
106       ALLOCATE(Dcoef_H(klon,klev), STAT=ierr)
107       IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_H, ierr=', ierr
108       
109       ALLOCATE(Acoef_Q(klon), Bcoef_Q(klon), Acoef_H(klon), Bcoef_H(klon), STAT=ierr)
110       IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_X and Bcoef_X, ierr=', ierr
111       
112       ALLOCATE(Kcoefhq(klon,klev), STAT=ierr)
113       IF ( ierr /= 0 )  PRINT*,' pb in allloc Kcoefhq, ierr=', ierr
114       
115       ALLOCATE(gamaq(1:klon,2:klev), STAT=ierr)
116       IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaq, ierr=', ierr
117       
118       ALLOCATE(gamah(1:klon,2:klev), STAT=ierr)
119       IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr
120       
121       ALLOCATE(h_old(klon,klev), STAT=ierr)
122       IF ( ierr /= 0 )  PRINT*,' pb in allloc h_old, ierr=', ierr
123       
124       ALLOCATE(d_h_col_vdf(klon), STAT=ierr)
125       IF ( ierr /= 0 )  PRINT*,' pb in allloc d_h_col_vdf, ierr=', ierr
126       
127       ALLOCATE(f_h_bnd(klon), STAT=ierr)
128       IF ( ierr /= 0 )  PRINT*,' pb in allloc f_h_bnd, ierr=', ierr
129    END IF
130
131!****************************************************************************************
132! 2)
133! Definition of the coeficient K
134
135!****************************************************************************************
136    Kcoefhq(:,:) = 0.0
137    DO k = 2, klev
138       DO i = 1, knon
139          Kcoefhq(i,k) = &
140               coefhq(i,k)*RG*RG*dtime /(pplay(i,k-1)-pplay(i,k)) &
141               *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
142       ENDDO
143    ENDDO
144
145!****************************************************************************************
146! 3)
147! Calculation of gama for "Q" and "H"
148
149!****************************************************************************************
150!   surface pressure is used as reference
151    psref(:) = paprs(:,1)
152
153!   definition of gama
154    IF (iflag_pbl == 1) THEN
155       gamaq(:,:) = 0.0
156       gamah(:,:) = -1.0e-03
157       gamah(:,2) = -2.5e-03
158 
159! conversion de gama
160       DO k = 2, klev
161          DO i = 1, knon
162             delz = RD * (temp(i,k-1)+temp(i,k)) / &
163                    2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
164             pkh  = (psref(i)/paprs(i,k))**RKAPPA
165         
166! convertie gradient verticale d'humidite specifique en difference d'humidite specifique entre centre de couches
167             gamaq(i,k) = gamaq(i,k) * delz   
168! convertie gradient verticale de temperature en difference de temperature potentielle entre centre de couches
169             gamah(i,k) = gamah(i,k) * delz * RCPD * pkh
170          ENDDO
171       ENDDO
172
173    ELSE
174       gamaq(:,:) = 0.0
175       gamah(:,:) = 0.0
176    ENDIF
177   
178
179!****************************************************************************************   
180! 4)
181! Calculte the coefficients C and D for specific humidity, q
182
183!****************************************************************************************
184   
185    CALL calc_coef(knon, Kcoefhq(:,:), gamaq(:,:), delp(:,:), q(:,:), &
186         Ccoef_Q(:,:), Dcoef_Q(:,:), Acoef_Q, Bcoef_Q)
187
188!****************************************************************************************
189! 5)
190! Calculte the coefficients C and D for potentiel entalpie, H
191
192!****************************************************************************************
193    h_old(:,:) = 0.0
194
195    DO k=1,klev
196       DO i = 1, knon
197          ! convertie la temperature en entalpie potentielle
198          h_old(i,k) = RCPD * temp(i,k) * &
199               (psref(i)/pplay(i,k))**RKAPPA
200       ENDDO
201    ENDDO
202
203    CALL calc_coef(knon, Kcoefhq(:,:), gamah(:,:), delp(:,:), h_old(:,:), &
204         Ccoef_H(:,:), Dcoef_H(:,:), Acoef_H, Bcoef_H)
205 
206!****************************************************************************************
207! 6)
208! Return the first layer in output variables
209
210!****************************************************************************************
211    Acoef_H_out = Acoef_H
212    Bcoef_H_out = Bcoef_H
213    Acoef_Q_out = Acoef_Q
214    Bcoef_Q_out = Bcoef_Q
215
216!****************************************************************************************
217! 7)
218! If Pbl is split, return also the other layers in output variables
219
220!****************************************************************************************
221!!! jyg le 07/02/2012
222!!jyg       IF (mod(iflag_pbl_split,2) .EQ.1) THEN
223       IF (mod(iflag_pbl_split,10) >=1) THEN
224!!! nrlmd le 02/05/2011
225    DO k= 1, klev
226      DO i= 1, klon
227        Ccoef_H_out(i,k) = Ccoef_H(i,k)
228        Dcoef_H_out(i,k) = Dcoef_H(i,k)
229        Ccoef_Q_out(i,k) = Ccoef_Q(i,k)
230        Dcoef_Q_out(i,k) = Dcoef_Q(i,k)
231        Kcoef_hq_out(i,k) = Kcoefhq(i,k)
232          IF (k==1) THEN
233            gama_h_out(i,k)  = 0.
234            gama_q_out(i,k)  = 0.
235          ELSE
236            gama_h_out(i,k)  = gamah(i,k)
237            gama_q_out(i,k)  = gamaq(i,k)
238          ENDIF
239      ENDDO
240    ENDDO
241!!!     
242       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
243!!!
244
245  END SUBROUTINE climb_hq_down
246
247!****************************************************************************************
248
249  SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
250
251! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
252! where X is H or Q, and k the vertical level k=1,klev
253
254    INCLUDE "YOMCST.h"
255! Input arguments
256!****************************************************************************************
257    INTEGER, INTENT(IN)                      :: knon
258    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, delp
259    REAL, DIMENSION(klon,klev), INTENT(IN)   :: X
260    REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama
261
262! Output arguments
263!****************************************************************************************
264    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef, Bcoef
265    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef, Dcoef
266
267! Local variables
268!****************************************************************************************
269    INTEGER                                  :: k, i
270    REAL                                     :: buf
271
272!****************************************************************************************
273! Niveau au sommet, k=klev
274
275!****************************************************************************************
276    Ccoef(:,:) = 0.0
277    Dcoef(:,:) = 0.0
278
279    DO i = 1, knon
280       buf = delp(i,klev) + Kcoef(i,klev)
281       
282       Ccoef(i,klev) = (X(i,klev)*delp(i,klev) - Kcoef(i,klev)*gama(i,klev))/buf
283       Dcoef(i,klev) = Kcoef(i,klev)/buf
284    END DO
285
286
287!****************************************************************************************
288! Niveau  (klev-1) <= k <= 2
289
290!****************************************************************************************
291
292    DO k=(klev-1),2,-1
293       DO i = 1, knon
294          buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
295          Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1) + &
296               Kcoef(i,k+1)*gama(i,k+1) - Kcoef(i,k)*gama(i,k))/buf
297          Dcoef(i,k) = Kcoef(i,k)/buf
298       END DO
299    END DO
300
301!****************************************************************************************
302! Niveau k=1
303
304!****************************************************************************************
305
306    DO i = 1, knon
307       buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
308       Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*(gama(i,2)+Ccoef(i,2)))/buf
309       Bcoef(i) = -1. * RG / buf
310    END DO
311
312  END SUBROUTINE calc_coef
313
314!****************************************************************************************
315
316  SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, &
317       flx_q1, flx_h1, paprs, pplay, &
318!!! nrlmd le 02/05/2011
319       Acoef_H_in, Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in, &
320       Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in, &
321       Kcoef_hq_in, gama_q_in, gama_h_in, &
322!!!
323       flux_q, flux_h, d_q, d_t)
324
325! This routine calculates the flux and tendency of the specific humidity q and
326! the potential engergi H.
327! The quantities q and H are calculated according to
328! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
329! C and D are known from before and k is index of the vertical layer.
330
331! Input arguments
332!****************************************************************************************
333    USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
334    IMPLICIT NONE
335    INTEGER, INTENT(IN)                      :: knon
336    REAL, INTENT(IN)                         :: dtime
337    REAL, DIMENSION(klon,klev), INTENT(IN)   :: t_old, q_old
338    REAL, DIMENSION(klon), INTENT(IN)        :: flx_q1, flx_h1
339    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
340    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
341
342!!! nrlmd le 02/05/2011
343    REAL, DIMENSION(klon), INTENT(IN)        :: Acoef_H_in,Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in
344    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in
345    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef_hq_in, gama_q_in, gama_h_in
346!!!
347
348! Output arguments
349!****************************************************************************************
350    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: flux_q, flux_h, d_q, d_t
351
352! Local variables
353!****************************************************************************************
354    LOGICAL, SAVE                            :: last=.FALSE.
355    REAL, DIMENSION(klon,klev)               :: h_new, q_new
356    REAL, DIMENSION(klon)                    :: psref         
357    INTEGER                                  :: k, i, ierr
358 
359! Include
360!****************************************************************************************
361    INCLUDE "YOMCST.h"
362!****************************************************************************************
363! 1)
364! Definition of some variables
365    REAL, DIMENSION(klon,klev)               :: d_h, zairm
366
367!****************************************************************************************
368    flux_q(:,:) = 0.0
369    flux_h(:,:) = 0.0
370    d_q(:,:)    = 0.0
371    d_t(:,:)    = 0.0
372    d_h(:,:)    = 0.0
373    f_h_bnd(:)= 0.0
374
375    psref(1:knon) = paprs(1:knon,1) 
376
377!!! jyg le 07/02/2012
378!!jyg       IF (mod(iflag_pbl_split,2) .EQ.1) THEN
379       IF (mod(iflag_pbl_split,10) >=1) THEN
380!!! nrlmd le 02/05/2011
381    DO i = 1, knon
382      Acoef_H(i)=Acoef_H_in(i)
383      Acoef_Q(i)=Acoef_Q_in(i)
384      Bcoef_H(i)=Bcoef_H_in(i)
385      Bcoef_Q(i)=Bcoef_Q_in(i)
386    ENDDO
387    DO k = 1, klev
388      DO i = 1, knon
389        Ccoef_H(i,k)=Ccoef_H_in(i,k)
390        Ccoef_Q(i,k)=Ccoef_Q_in(i,k)
391        Dcoef_H(i,k)=Dcoef_H_in(i,k)
392        Dcoef_Q(i,k)=Dcoef_Q_in(i,k)
393        Kcoefhq(i,k)=Kcoef_hq_in(i,k)
394          IF (k>1) THEN
395            gamah(i,k)=gama_h_in(i,k)
396            gamaq(i,k)=gama_q_in(i,k)
397          ENDIF
398      ENDDO
399    ENDDO
400!!!     
401       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
402!!!
403
404!****************************************************************************************
405! 2)
406! Calculation of Q and H
407
408!****************************************************************************************
409
410!- First layer
411    q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)*flx_q1(1:knon)*dtime
412    h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)*flx_h1(1:knon)*dtime
413    f_h_bnd(1:knon) = flx_h1(1:knon)
414!- All the other layers
415    DO k = 2, klev
416       DO i = 1, knon
417          q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1)
418          h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1)
419       END DO
420    END DO
421!****************************************************************************************
422! 3)
423! Calculation of the flux for Q and H
424
425!****************************************************************************************
426
427!- The flux at first layer, k=1
428    flux_q(1:knon,1)=flx_q1(1:knon)
429    flux_h(1:knon,1)=flx_h1(1:knon)
430
431!- The flux at all layers above surface
432    DO k = 2, klev
433       DO i = 1, knon
434          flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
435               (q_new(i,k)-q_new(i,k-1)+gamaq(i,k))
436
437          flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
438               (h_new(i,k)-h_new(i,k-1)+gamah(i,k))
439       END DO
440    END DO
441
442!****************************************************************************************
443! 4)
444! Calculation of tendency for Q and H
445
446!****************************************************************************************
447    d_h_col_vdf(:) = 0.0
448    DO k = 1, klev
449       DO i = 1, knon
450          d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
451          d_q(i,k) = q_new(i,k) - q_old(i,k)
452          d_h(i,k) = h_new(i,k) - h_old(i,k)
453!JLD          d_t(i,k) = d_h(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD    !correction a venir
454    ! layer air mass
455          zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
456          d_h_col_vdf(i) = d_h_col_vdf(i) + d_h(i,k)*zairm(i,k)
457        END DO
458    END DO
459
460!****************************************************************************************
461! Some deallocations
462
463!****************************************************************************************
464    IF (last) THEN
465       DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)   
466       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
467       DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)   
468       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
469       DEALLOCATE(gamaq, gamah,stat=ierr)
470       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
471       DEALLOCATE(Kcoefhq,stat=ierr)
472       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
473       DEALLOCATE(h_old, d_h_col_vdf, f_h_bnd, stat=ierr)
474       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate h_old, d_h_col_vdf, f_h_bnd, ierr=', ierr
475    END IF
476  END SUBROUTINE climb_hq_up
477
478!****************************************************************************************
479
480END MODULE climb_hq_mod
481
482 
483
484
485
486
Note: See TracBrowser for help on using the repository browser.