source: LMDZ6/trunk/libf/phylmd/climb_hq_mod.f90 @ 5408

Last change on this file since 5408 was 5400, checked in by evignon, 5 days ago

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