source: LMDZ6/branches/PBLSURF_GPUPORT/libf/phylmd/climb_hq_mod.f90

Last change on this file was 5878, checked in by yann meurdesoif, 3 days ago

GPU port of climb_hq_up +climb_hq_down
YM
  • 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: 20.5 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  USE compbl_mod_h
8  IMPLICIT NONE
9  PRIVATE
10  PUBLIC :: climb_hq_down, climb_hq_up, climb_hq_init, climb_hq_finalize
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_init
36  USE dimphy, ONLY : klon, klev
37  IMPLICIT NONE
38    INTEGER :: ierr
39 
40     ALLOCATE(Ccoef_Q(klon,klev), STAT=ierr)
41     IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_Q, ierr=', ierr
42     Ccoef_Q(:,:) = 0.
43
44     ALLOCATE(Dcoef_Q(klon,klev), STAT=ierr)
45     IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_Q, ierr=', ierr
46     Dcoef_Q(:,:) = 0.
47
48     ALLOCATE(Ccoef_H(klon,klev), STAT=ierr)
49     IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_H, ierr=', ierr
50     Ccoef_H(:,:) = 0.
51
52     ALLOCATE(Dcoef_H(klon,klev), STAT=ierr)
53     IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_H, ierr=', ierr
54     Dcoef_H(:,:) = 0.
55
56     ALLOCATE(Acoef_Q(klon), Bcoef_Q(klon), Acoef_H(klon), Bcoef_H(klon), STAT=ierr)
57     IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_X and Bcoef_X, ierr=', ierr
58     Acoef_Q(:)=0. ; Bcoef_Q(:)=0. ; Acoef_H(:)=0. ; Bcoef_H(:)=0. ;
59
60     ALLOCATE(Kcoefhq(klon,klev), STAT=ierr)
61     IF ( ierr /= 0 )  PRINT*,' pb in allloc Kcoefhq, ierr=', ierr
62     Kcoefhq(:,:) = 0.
63
64     ALLOCATE(gamaq(1:klon,2:klev), STAT=ierr)
65     IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaq, ierr=', ierr
66     gamaq(:,:) = 0.
67
68     ALLOCATE(gamah(1:klon,2:klev), STAT=ierr)
69     IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr
70     gamah(:,:)=0.
71
72     ALLOCATE(h_old(klon,klev), STAT=ierr)
73     IF ( ierr /= 0 )  PRINT*,' pb in allloc h_old, ierr=', ierr
74     h_old(:,:) = 0.
75
76     ALLOCATE(d_h_col_vdf(klon), STAT=ierr)
77     IF ( ierr /= 0 )  PRINT*,' pb in allloc d_h_col_vdf, ierr=', ierr
78     d_h_col_vdf(:) = 0.
79
80     ALLOCATE(f_h_bnd(klon), STAT=ierr)
81     IF ( ierr /= 0 )  PRINT*,' pb in allloc f_h_bnd, ierr=', ierr
82     f_h_bnd(:) = 0.
83
84  END SUBROUTINE climb_hq_init
85
86  SUBROUTINE climb_hq_finalize
87  IMPLICIT NONE
88    INTEGER :: ierr
89
90    !****************************************************************************************
91    ! Some deallocations
92    !
93    !****************************************************************************************
94   
95    DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)   
96    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
97    DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)   
98    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
99    DEALLOCATE(gamaq, gamah,stat=ierr)
100    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
101    DEALLOCATE(Kcoefhq,stat=ierr)
102    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
103    DEALLOCATE(h_old, d_h_col_vdf, f_h_bnd, stat=ierr)
104    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate h_old, d_h_col_vdf, f_h_bnd, ierr=', ierr
105
106  END  SUBROUTINE climb_hq_finalize
107
108  SUBROUTINE climb_hq_down(knon, ni, coefhq, paprs, pplay, &
109       delp, temp, q, dtime, &
110!!! nrlmd le 02/05/2011
111       Ccoef_H_out, Ccoef_Q_out, Dcoef_H_out, Dcoef_Q_out, &
112       Kcoef_hq_out, gama_q_out, gama_h_out, &
113!!!
114       Acoef_H_out, Acoef_Q_out, Bcoef_H_out, Bcoef_Q_out)
115!$gpum horizontal knon
116
117! This routine calculates recursivly the coefficients C and D
118! for the quantity X=[Q,H] in equation X(k) = C(k) + D(k)*X(k-1), where k is
119! the index of the vertical layer.
120   USE yomcst_mod_h
121! Input arguments
122!****************************************************************************************
123    INTEGER, INTENT(IN)                      :: knon
124    INTEGER, INTENT(IN)                      :: ni(knon)
125    REAL, DIMENSION(knon,klev), INTENT(IN)   :: coefhq
126    REAL, DIMENSION(knon,klev), INTENT(IN)   :: pplay
127    REAL, DIMENSION(knon,klev+1), INTENT(IN) :: paprs
128    REAL, DIMENSION(knon,klev), INTENT(IN)   :: temp, delp  ! temperature
129    REAL, DIMENSION(knon,klev), INTENT(IN)   :: q
130    REAL, INTENT(IN)                         :: dtime
131
132! Output arguments
133!****************************************************************************************
134    REAL, DIMENSION(knon), INTENT(OUT)       :: Acoef_H_out
135    REAL, DIMENSION(knon), INTENT(OUT)       :: Acoef_Q_out
136    REAL, DIMENSION(knon), INTENT(OUT)       :: Bcoef_H_out
137    REAL, DIMENSION(knon), INTENT(OUT)       :: Bcoef_Q_out
138
139!!! nrlmd le 02/05/2011
140    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Ccoef_H_out
141    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Ccoef_Q_out
142    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Dcoef_H_out
143    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Dcoef_Q_out
144    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Kcoef_hq_out
145    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: gama_q_out
146    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: gama_h_out
147
148    REAL :: yCcoef_Q(knon,klev)
149    REAL :: yDcoef_Q(knon,klev)
150    REAL :: yCcoef_H(knon,klev)
151    REAL :: yDcoef_H(knon,klev)
152    REAL :: yAcoef_Q(knon), yBcoef_Q(knon), yAcoef_H(knon), yBcoef_H(knon)
153    REAL :: yKcoefhq(knon,klev)
154    REAL :: ygamaq(knon,2:klev)
155    REAL :: ygamah(knon,2:klev)
156    REAL :: yh_old(knon,klev) ! for diagnostics, h before solving diffusion
157    REAL :: yd_h_col_vdf(knon) ! for diagnostics, vertical integral of enthalpy change
158    REAL :: yf_h_bnd(knon) ! for diagnostics, enthalpy flux at surface
159!!!
160
161! Local variables
162!****************************************************************************************
163! JLD now renamed h_old and declared in module
164!    REAL, DIMENSION(knon,klev)               :: local_H
165    REAL, DIMENSION(knon)                    :: psref
166    REAL                                     :: delz, pkh
167    INTEGER                                  :: k, i, j, ierr
168
169
170    yd_h_col_vdf(:) =0.
171    yf_h_bnd(:) = 0.
172!****************************************************************************************
173! 2)
174! Definition of the coeficient K
175!
176!****************************************************************************************
177    yKcoefhq(:,:) = 0.0
178    DO k = 2, klev
179       DO i = 1, knon
180          yKcoefhq(i,k) = &
181               coefhq(i,k)*RG*RG*dtime /(pplay(i,k-1)-pplay(i,k)) &
182               *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
183       ENDDO
184    ENDDO
185
186!****************************************************************************************
187! 3)
188! Calculation of gama for "Q" and "H"
189!
190!****************************************************************************************
191!   surface pressure is used as reference
192    psref(:) = paprs(:,1)
193
194!   definition of gama
195    IF (iflag_pbl == 1) THEN
196       ygamaq(:,:) = 0.0
197       ygamah(:,:) = -1.0e-03
198       ygamah(:,2) = -2.5e-03
199 
200! conversion de gama
201       DO k = 2, klev
202          DO i = 1, knon
203             delz = RD * (temp(i,k-1)+temp(i,k)) / &
204                    2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
205             pkh  = (psref(i)/paprs(i,k))**RKAPPA
206         
207! convertie gradient verticale d'humidite specifique en difference d'humidite specifique entre centre de couches
208             ygamaq(i,k) = ygamaq(i,k) * delz   
209! convertie gradient verticale de temperature en difference de temperature potentielle entre centre de couches
210             ygamah(i,k) = ygamah(i,k) * delz * RCPD * pkh
211          ENDDO
212       ENDDO
213
214    ELSE
215       ygamaq(:,:) = 0.0
216       ygamah(:,:) = 0.0
217    ENDIF
218   
219
220!****************************************************************************************   
221! 4)
222! Calculte the coefficients C and D for specific humidity, q
223!
224!****************************************************************************************
225   
226    CALL calc_coef(knon, yKcoefhq, ygamaq, delp, q, &
227         yCcoef_Q, yDcoef_Q, yAcoef_Q, yBcoef_Q)
228
229!****************************************************************************************
230! 5)
231! Calculte the coefficients C and D for potentiel entalpie, H
232!
233!****************************************************************************************
234    yh_old(:,:) = 0.0
235
236    DO k=1,klev
237       DO i = 1, knon
238          ! convertie la temperature en entalpie potentielle
239          yh_old(i,k) = RCPD * temp(i,k) * &
240               (psref(i)/pplay(i,k))**RKAPPA
241       ENDDO
242    ENDDO
243
244    CALL calc_coef(knon, yKcoefhq, ygamah, delp, yh_old, &
245         yCcoef_H, yDcoef_H, yAcoef_H, yBcoef_H)
246 
247!****************************************************************************************
248! 6)
249! Return the first layer in output variables
250!
251!****************************************************************************************
252    Acoef_H_out = yAcoef_H
253    Bcoef_H_out = yBcoef_H
254    Acoef_Q_out = yAcoef_Q
255    Bcoef_Q_out = yBcoef_Q
256
257!****************************************************************************************
258! 7)
259! If Pbl is split, return also the other layers in output variables
260!
261!****************************************************************************************
262!!! jyg le 07/02/2012
263!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
264       IF (mod(iflag_pbl_split,10) .ge.1) THEN
265!!! nrlmd le 02/05/2011
266    DO k= 1, klev
267      DO i= 1, knon
268        Ccoef_H_out(i,k) = yCcoef_H(i,k)
269        Dcoef_H_out(i,k) = yDcoef_H(i,k)
270        Ccoef_Q_out(i,k) = yCcoef_Q(i,k)
271        Dcoef_Q_out(i,k) = yDcoef_Q(i,k)
272        Kcoef_hq_out(i,k) = yKcoefhq(i,k)
273          IF (k.eq.1) THEN
274            gama_h_out(i,k)  = 0.
275            gama_q_out(i,k)  = 0.
276          ELSE
277            gama_h_out(i,k)  = ygamah(i,k)
278            gama_q_out(i,k)  = ygamaq(i,k)
279          ENDIF
280      ENDDO
281    ENDDO
282!!!     
283       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
284!!!
285    DO k=1,klev
286      DO j=1,knon
287        i=ni(j)
288        IF (k==1) THEN
289          Acoef_Q(i) = yAcoef_Q(j)
290          Bcoef_Q(i) = yBcoef_Q(j)
291          Acoef_H(i) = yAcoef_H(j)
292          Bcoef_H(i) = yBcoef_H(j)
293          d_h_col_vdf(i)= yd_h_col_vdf(j)
294          f_h_bnd(i)= yf_h_bnd(j)
295        ENDIF
296       
297        IF (k>=2) THEN
298          gamaq(i,k)=ygamaq(j,k)
299          gamah(i,k)=ygamah(j,k)
300        ENDIF
301
302        Ccoef_Q(i,k) = yCcoef_Q(j,k)
303        Dcoef_Q(i,k) = yDcoef_Q(j,k)
304        Ccoef_H(i,k) = yCcoef_H(j,k)
305        Dcoef_H(i,k) = yDcoef_H(j,k)
306        Kcoefhq(i,k) = yKcoefhq(j,k)
307        h_old(i,k) = yh_old(j,k)
308      ENDDO
309    ENDDO
310
311  END SUBROUTINE climb_hq_down
312!
313!****************************************************************************************
314!
315  SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
316!$gpum horizontal knon
317!
318! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
319! where X is H or Q, and k the vertical level k=1,klev
320!
321USE yomcst_mod_h! Input arguments
322!****************************************************************************************
323    INTEGER, INTENT(IN)                      :: knon
324    REAL, DIMENSION(knon,klev), INTENT(IN)   :: Kcoef, delp
325    REAL, DIMENSION(knon,klev), INTENT(IN)   :: X
326    REAL, DIMENSION(knon,2:klev), INTENT(IN) :: gama
327
328! Output arguments
329!****************************************************************************************
330    REAL, DIMENSION(knon), INTENT(OUT)       :: Acoef, Bcoef
331    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Ccoef, Dcoef
332
333! Local variables
334!****************************************************************************************
335    INTEGER                                  :: k, i
336    REAL                                     :: buf
337
338!****************************************************************************************
339! Niveau au sommet, k=klev
340!
341!****************************************************************************************
342    Ccoef(:,:) = 0.0
343    Dcoef(:,:) = 0.0
344
345    DO i = 1, knon
346       buf = delp(i,klev) + Kcoef(i,klev)
347       
348       Ccoef(i,klev) = (X(i,klev)*delp(i,klev) - Kcoef(i,klev)*gama(i,klev))/buf
349       Dcoef(i,klev) = Kcoef(i,klev)/buf
350    END DO
351
352
353!****************************************************************************************
354! Niveau  (klev-1) <= k <= 2
355!
356!****************************************************************************************
357
358    DO k=(klev-1),2,-1
359       DO i = 1, knon
360          buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
361          Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1) + &
362               Kcoef(i,k+1)*gama(i,k+1) - Kcoef(i,k)*gama(i,k))/buf
363          Dcoef(i,k) = Kcoef(i,k)/buf
364       END DO
365    END DO
366
367!****************************************************************************************
368! Niveau k=1
369!
370!****************************************************************************************
371
372    DO i = 1, knon
373       buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
374       Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*(gama(i,2)+Ccoef(i,2)))/buf
375       Bcoef(i) = -1. * RG / buf
376    END DO
377
378  END SUBROUTINE calc_coef
379!
380!****************************************************************************************
381!
382  SUBROUTINE climb_hq_up(knon, ni, dtime, t_old, q_old, &
383       flx_q1, flx_h1, paprs, pplay, &
384!!! nrlmd le 02/05/2011
385       Acoef_H_in, Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in, &
386       Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in, &
387       Kcoef_hq_in, gama_q_in, gama_h_in, &
388!!!
389       flux_q, flux_h, d_q, d_t)
390!$gpum horizontal knon
391
392!
393! This routine calculates the flux and tendency of the specific humidity q and
394! the potential engergi H.
395! The quantities q and H are calculated according to
396! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
397! C and D are known from before and k is index of the vertical layer.
398!   
399USE yomcst_mod_h
400USE compbl_mod_h
401! Input arguments
402!****************************************************************************************
403    INTEGER, INTENT(IN)                      :: knon
404    INTEGER, INTENT(IN)                      :: ni(knon)
405    REAL, INTENT(IN)                         :: dtime
406    REAL, DIMENSION(knon,klev), INTENT(IN)   :: t_old, q_old
407    REAL, DIMENSION(knon), INTENT(IN)        :: flx_q1, flx_h1
408    REAL, DIMENSION(knon,klev+1), INTENT(IN) :: paprs
409    REAL, DIMENSION(knon,klev), INTENT(IN)   :: pplay
410
411!!! nrlmd le 02/05/2011
412    REAL, DIMENSION(knon), INTENT(IN)        :: Acoef_H_in,Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in
413    REAL, DIMENSION(knon,klev), INTENT(IN)   :: Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in
414    REAL, DIMENSION(knon,klev), INTENT(IN)   :: Kcoef_hq_in, gama_q_in, gama_h_in
415!!!
416
417! Output arguments
418!****************************************************************************************
419    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: flux_q, flux_h, d_q, d_t
420
421! Local variables
422!****************************************************************************************
423    REAL :: yCcoef_Q(knon,klev)
424    REAL :: yDcoef_Q(knon,klev)
425    REAL :: yCcoef_H(knon,klev)
426    REAL :: yDcoef_H(knon,klev)
427    REAL :: yAcoef_Q(knon), yBcoef_Q(knon), yAcoef_H(knon), yBcoef_H(knon)
428    REAL :: yKcoefhq(knon,klev)
429    REAL :: ygamaq(knon,2:klev)
430    REAL :: ygamah(knon,2:klev)
431    REAL :: yh_old(knon,klev)
432    REAL :: yd_h_col_vdf(knon)
433    REAL :: yf_h_bnd(knon)
434
435    REAL, DIMENSION(knon,klev)               :: h_new, q_new
436    REAL, DIMENSION(knon)                    :: psref         
437    INTEGER                                  :: k, i, j, ierr
438!****************************************************************************************
439! 1)
440! Definition of some variables
441    REAL, DIMENSION(knon,klev)               :: d_h, zairm
442!
443!****************************************************************************************
444    flux_q(:,:) = 0.0
445    flux_h(:,:) = 0.0
446    d_q(:,:)    = 0.0
447    d_t(:,:)    = 0.0
448    d_h(:,:)    = 0.0
449
450    psref(1:knon) = paprs(1:knon,1) 
451
452    DO k=1,klev
453      DO j=1,knon
454        i=ni(j)
455        IF (k==1) THEN
456          yAcoef_Q(j) = Acoef_Q(i)
457          yBcoef_Q(j) = Bcoef_Q(i)
458          yAcoef_H(j) = Acoef_H(i)
459          yBcoef_H(j) = Bcoef_H(i)
460          yd_h_col_vdf(j)= d_h_col_vdf(i)
461          yf_h_bnd(j)= f_h_bnd(i)
462        ENDIF
463       
464        IF (k>=2) THEN
465          ygamaq(j,k)=gamaq(i,k)
466          ygamah(j,k)=gamah(i,k)
467        ENDIF
468        yCcoef_Q(j,k) = Ccoef_Q(i,k)
469        yDcoef_Q(j,k) = Dcoef_Q(i,k)
470        yCcoef_H(j,k) = Ccoef_H(i,k)
471        yDcoef_H(j,k) = Dcoef_H(i,k)
472        yKcoefhq(j,k) = Kcoefhq(i,k)
473        yh_old(j,k) = h_old(i,k)
474      ENDDO
475    ENDDO     
476    yf_h_bnd(:)= 0.0
477!!! jyg le 07/02/2012
478!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
479       IF (mod(iflag_pbl_split,10) .ge.1) THEN
480!!! nrlmd le 02/05/2011
481    DO i = 1, knon
482      yAcoef_H(i)=Acoef_H_in(i)
483      yAcoef_Q(i)=Acoef_Q_in(i)
484      yBcoef_H(i)=Bcoef_H_in(i)
485      yBcoef_Q(i)=Bcoef_Q_in(i)
486    ENDDO
487    DO k = 1, klev
488      DO i = 1, knon
489        yCcoef_H(i,k)=Ccoef_H_in(i,k)
490        yCcoef_Q(i,k)=Ccoef_Q_in(i,k)
491        yDcoef_H(i,k)=Dcoef_H_in(i,k)
492        yDcoef_Q(i,k)=Dcoef_Q_in(i,k)
493        yKcoefhq(i,k)=Kcoef_hq_in(i,k)
494          IF (k.gt.1) THEN
495            ygamah(i,k)=gama_h_in(i,k)
496            ygamaq(i,k)=gama_q_in(i,k)
497          ENDIF
498      ENDDO
499    ENDDO
500   
501!!!     
502       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
503!!!
504
505!****************************************************************************************
506! 2)
507! Calculation of Q and H
508!
509!****************************************************************************************
510
511!- First layer
512    q_new(1:knon,1) = yAcoef_Q(1:knon) + yBcoef_Q(1:knon)*flx_q1(1:knon)*dtime
513    h_new(1:knon,1) = yAcoef_H(1:knon) + yBcoef_H(1:knon)*flx_h1(1:knon)*dtime
514    yf_h_bnd(1:knon) = flx_h1(1:knon)
515!- All the other layers
516    DO k = 2, klev
517       DO i = 1, knon
518          q_new(i,k) = yCcoef_Q(i,k) + yDcoef_Q(i,k)*q_new(i,k-1)
519          h_new(i,k) = yCcoef_H(i,k) + yDcoef_H(i,k)*h_new(i,k-1)
520       END DO
521    END DO
522!****************************************************************************************
523! 3)
524! Calculation of the flux for Q and H
525!
526!****************************************************************************************
527
528!- The flux at first layer, k=1
529    flux_q(1:knon,1)=flx_q1(1:knon)
530    flux_h(1:knon,1)=flx_h1(1:knon)
531
532!- The flux at all layers above surface
533    DO k = 2, klev
534       DO i = 1, knon
535          flux_q(i,k) = (yKcoefhq(i,k)/RG/dtime) * &
536               (q_new(i,k)-q_new(i,k-1)+ygamaq(i,k))
537
538          flux_h(i,k) = (yKcoefhq(i,k)/RG/dtime) * &
539               (h_new(i,k)-h_new(i,k-1)+ygamah(i,k))
540       END DO
541    END DO
542
543!****************************************************************************************
544! 4)
545! Calculation of tendency for Q and H
546!
547!****************************************************************************************
548    yd_h_col_vdf(:) = 0.0
549    DO k = 1, klev
550       DO i = 1, knon
551          d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
552          d_q(i,k) = q_new(i,k) - q_old(i,k)
553          d_h(i,k) = h_new(i,k) - yh_old(i,k)
554!JLD          d_t(i,k) = d_h(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD    !correction a venir
555    ! layer air mass
556          zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
557          yd_h_col_vdf(i) = yd_h_col_vdf(i) + d_h(i,k)*zairm(i,k)
558        END DO
559    END DO
560
561    DO k=1,klev
562      DO j=1,knon
563        i=ni(j)
564        IF (k==1) THEN
565          Acoef_Q(i) = yAcoef_Q(j)
566          Bcoef_Q(i) = yBcoef_Q(j)
567          Acoef_H(i) = yAcoef_H(j)
568          Bcoef_H(i) = yBcoef_H(j)
569          d_h_col_vdf(i)= yd_h_col_vdf(j)
570          f_h_bnd(i)= yf_h_bnd(j)
571        ENDIF
572       
573        IF (k>=2) THEN
574          gamaq(i,k)=ygamaq(j,k)
575          gamah(i,k)=ygamah(j,k)
576        ENDIF
577
578        Ccoef_Q(i,k) = yCcoef_Q(j,k)
579        Dcoef_Q(i,k) = yDcoef_Q(j,k)
580        Ccoef_H(i,k) = yCcoef_H(j,k)
581        Dcoef_H(i,k) = yDcoef_H(j,k)
582        Kcoefhq(i,k) = yKcoefhq(j,k)
583        h_old(i,k) = yh_old(j,k)
584      ENDDO
585    ENDDO
586
587  END SUBROUTINE climb_hq_up
588!
589!****************************************************************************************
590!
591END MODULE climb_hq_mod
592
593 
594
595
596
597
Note: See TracBrowser for help on using the repository browser.