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

Last change on this file since 5725 was 5144, checked in by abarral, 16 months ago

Put YOMCST.h into modules
  • 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.9 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    USE lmdz_yomcst
51
52    IMPLICIT NONE
53    INTEGER, INTENT(IN) :: knon
54    REAL, DIMENSION(klon, klev), INTENT(IN) :: coefhq
55    REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
56    REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs
57    REAL, DIMENSION(klon, klev), INTENT(IN) :: temp, delp  ! temperature
58    REAL, DIMENSION(klon, klev), INTENT(IN) :: q
59    REAL, INTENT(IN) :: dtime
60
61    ! Output arguments
62    !****************************************************************************************
63    REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_H_out
64    REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_Q_out
65    REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_H_out
66    REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_Q_out
67
68    !!! nrlmd le 02/05/2011
69    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Ccoef_H_out
70    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Ccoef_Q_out
71    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Dcoef_H_out
72    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Dcoef_Q_out
73    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Kcoef_hq_out
74    REAL, DIMENSION(klon, klev), INTENT(OUT) :: gama_q_out
75    REAL, DIMENSION(klon, klev), INTENT(OUT) :: gama_h_out
76    !!!
77
78    ! Local variables
79    !****************************************************************************************
80    LOGICAL, SAVE :: first = .TRUE.
81    !$OMP THREADPRIVATE(first)
82    ! JLD now renamed h_old and declared in module
83    !    REAL, DIMENSION(klon,klev)               :: local_H
84    REAL, DIMENSION(klon) :: psref
85    REAL :: delz, pkh
86    INTEGER :: k, i, ierr
87    ! Include
88    !****************************************************************************************
89    ! 1)
90    ! Allocation at first time step only
91
92    !****************************************************************************************
93
94    IF (first) THEN
95      first = .FALSE.
96      ALLOCATE(Ccoef_Q(klon, klev), STAT = ierr)
97      IF (ierr /= 0)  PRINT*, ' pb in allloc Ccoef_Q, ierr=', ierr
98
99      ALLOCATE(Dcoef_Q(klon, klev), STAT = ierr)
100      IF (ierr /= 0)  PRINT*, ' pb in allloc Dcoef_Q, ierr=', ierr
101
102      ALLOCATE(Ccoef_H(klon, klev), STAT = ierr)
103      IF (ierr /= 0)  PRINT*, ' pb in allloc Ccoef_H, ierr=', ierr
104
105      ALLOCATE(Dcoef_H(klon, klev), STAT = ierr)
106      IF (ierr /= 0)  PRINT*, ' pb in allloc Dcoef_H, ierr=', ierr
107
108      ALLOCATE(Acoef_Q(klon), Bcoef_Q(klon), Acoef_H(klon), Bcoef_H(klon), STAT = ierr)
109      IF (ierr /= 0)  PRINT*, ' pb in allloc Acoef_X and Bcoef_X, ierr=', ierr
110
111      ALLOCATE(Kcoefhq(klon, klev), STAT = ierr)
112      IF (ierr /= 0)  PRINT*, ' pb in allloc Kcoefhq, ierr=', ierr
113
114      ALLOCATE(gamaq(1:klon, 2:klev), STAT = ierr)
115      IF (ierr /= 0) PRINT*, ' pb in allloc gamaq, ierr=', ierr
116
117      ALLOCATE(gamah(1:klon, 2:klev), STAT = ierr)
118      IF (ierr /= 0) PRINT*, ' pb in allloc gamah, ierr=', ierr
119
120      ALLOCATE(h_old(klon, klev), STAT = ierr)
121      IF (ierr /= 0)  PRINT*, ' pb in allloc h_old, ierr=', ierr
122
123      ALLOCATE(d_h_col_vdf(klon), STAT = ierr)
124      IF (ierr /= 0)  PRINT*, ' pb in allloc d_h_col_vdf, ierr=', ierr
125
126      ALLOCATE(f_h_bnd(klon), STAT = ierr)
127      IF (ierr /= 0)  PRINT*, ' pb in allloc f_h_bnd, ierr=', ierr
128    END IF
129
130    !****************************************************************************************
131    ! 2)
132    ! Definition of the coeficient K
133
134    !****************************************************************************************
135    Kcoefhq(:, :) = 0.0
136    DO k = 2, klev
137      DO i = 1, knon
138        Kcoefhq(i, k) = &
139                coefhq(i, k) * RG * RG * dtime / (pplay(i, k - 1) - pplay(i, k)) &
140                        * (paprs(i, k) * 2 / (temp(i, k) + temp(i, k - 1)) / RD)**2
141      ENDDO
142    ENDDO
143
144    !****************************************************************************************
145    ! 3)
146    ! Calculation of gama for "Q" and "H"
147
148    !****************************************************************************************
149    !   surface pressure is used as reference
150    psref(:) = paprs(:, 1)
151
152    !   definition of gama
153    IF (iflag_pbl == 1) THEN
154      gamaq(:, :) = 0.0
155      gamah(:, :) = -1.0e-03
156      gamah(:, 2) = -2.5e-03
157
158      ! conversion de gama
159      DO k = 2, klev
160        DO i = 1, knon
161          delz = RD * (temp(i, k - 1) + temp(i, k)) / &
162                  2.0 / RG / paprs(i, k) * (pplay(i, k - 1) - pplay(i, k))
163          pkh = (psref(i) / paprs(i, k))**RKAPPA
164
165          ! convertie gradient verticale d'humidite specifique en difference d'humidite specifique entre centre de couches
166          gamaq(i, k) = gamaq(i, k) * delz
167          ! convertie gradient verticale de temperature en difference de temperature potentielle entre centre de couches
168          gamah(i, k) = gamah(i, k) * delz * RCPD * pkh
169        ENDDO
170      ENDDO
171
172    ELSE
173      gamaq(:, :) = 0.0
174      gamah(:, :) = 0.0
175    ENDIF
176
177
178    !****************************************************************************************
179    ! 4)
180    ! Calculte the coefficients C and D for specific humidity, q
181
182    !****************************************************************************************
183
184    CALL calc_coef(knon, Kcoefhq(:, :), gamaq(:, :), delp(:, :), q(:, :), &
185            Ccoef_Q(:, :), Dcoef_Q(:, :), Acoef_Q, Bcoef_Q)
186
187    !****************************************************************************************
188    ! 5)
189    ! Calculte the coefficients C and D for potentiel entalpie, H
190
191    !****************************************************************************************
192    h_old(:, :) = 0.0
193
194    DO k = 1, klev
195      DO i = 1, knon
196        ! convertie la temperature en entalpie potentielle
197        h_old(i, k) = RCPD * temp(i, k) * &
198                (psref(i) / pplay(i, k))**RKAPPA
199      ENDDO
200    ENDDO
201
202    CALL calc_coef(knon, Kcoefhq(:, :), gamah(:, :), delp(:, :), h_old(:, :), &
203            Ccoef_H(:, :), Dcoef_H(:, :), Acoef_H, Bcoef_H)
204
205    !****************************************************************************************
206    ! 6)
207    ! Return the first layer in output variables
208
209    !****************************************************************************************
210    Acoef_H_out = Acoef_H
211    Bcoef_H_out = Bcoef_H
212    Acoef_Q_out = Acoef_Q
213    Bcoef_Q_out = Bcoef_Q
214
215    !****************************************************************************************
216    ! 7)
217    ! If Pbl is split, return also the other layers in output variables
218
219    !****************************************************************************************
220    !!! jyg le 07/02/2012
221    !!jyg       IF (mod(iflag_pbl_split,2) .EQ.1) THEN
222    IF (mod(iflag_pbl_split, 10) >=1) THEN
223      !!! nrlmd le 02/05/2011
224      DO k = 1, klev
225        DO i = 1, klon
226          Ccoef_H_out(i, k) = Ccoef_H(i, k)
227          Dcoef_H_out(i, k) = Dcoef_H(i, k)
228          Ccoef_Q_out(i, k) = Ccoef_Q(i, k)
229          Dcoef_Q_out(i, k) = Dcoef_Q(i, k)
230          Kcoef_hq_out(i, k) = Kcoefhq(i, k)
231          IF (k==1) THEN
232            gama_h_out(i, k) = 0.
233            gama_q_out(i, k) = 0.
234          ELSE
235            gama_h_out(i, k) = gamah(i, k)
236            gama_q_out(i, k) = gamaq(i, k)
237          ENDIF
238        ENDDO
239      ENDDO
240      !!!
241    ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
242    !!!
243
244  END SUBROUTINE climb_hq_down
245
246  !****************************************************************************************
247
248  SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
249    USE lmdz_yomcst
250
251    IMPLICIT NONE
252
253    ! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
254    ! where X is H or Q, and k the vertical level k=1,klev
255
256    ! Input arguments
257    !****************************************************************************************
258    INTEGER, INTENT(IN) :: knon
259    REAL, DIMENSION(klon, klev), INTENT(IN) :: Kcoef, delp
260    REAL, DIMENSION(klon, klev), INTENT(IN) :: X
261    REAL, DIMENSION(klon, 2:klev), INTENT(IN) :: gama
262
263    ! Output arguments
264    !****************************************************************************************
265    REAL, DIMENSION(klon), INTENT(OUT) :: Acoef, Bcoef
266    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Ccoef, Dcoef
267
268    ! Local variables
269    !****************************************************************************************
270    INTEGER :: k, i
271    REAL :: buf
272
273    !****************************************************************************************
274    ! Niveau au sommet, k=klev
275
276    !****************************************************************************************
277    Ccoef(:, :) = 0.0
278    Dcoef(:, :) = 0.0
279
280    DO i = 1, knon
281      buf = delp(i, klev) + Kcoef(i, klev)
282
283      Ccoef(i, klev) = (X(i, klev) * delp(i, klev) - Kcoef(i, klev) * gama(i, klev)) / buf
284      Dcoef(i, klev) = Kcoef(i, klev) / buf
285    END DO
286
287
288    !****************************************************************************************
289    ! Niveau  (klev-1) <= k <= 2
290
291    !****************************************************************************************
292
293    DO k = (klev - 1), 2, -1
294      DO i = 1, knon
295        buf = delp(i, k) + Kcoef(i, k) + Kcoef(i, k + 1) * (1. - Dcoef(i, k + 1))
296        Ccoef(i, k) = (X(i, k) * delp(i, k) + Kcoef(i, k + 1) * Ccoef(i, k + 1) + &
297                Kcoef(i, k + 1) * gama(i, k + 1) - Kcoef(i, k) * gama(i, k)) / buf
298        Dcoef(i, k) = Kcoef(i, k) / buf
299      END DO
300    END DO
301
302    !****************************************************************************************
303    ! Niveau k=1
304
305    !****************************************************************************************
306
307    DO i = 1, knon
308      buf = delp(i, 1) + Kcoef(i, 2) * (1. - Dcoef(i, 2))
309      Acoef(i) = (X(i, 1) * delp(i, 1) + Kcoef(i, 2) * (gama(i, 2) + Ccoef(i, 2))) / buf
310      Bcoef(i) = -1. * RG / buf
311    END DO
312
313  END SUBROUTINE calc_coef
314
315  !****************************************************************************************
316
317  SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, &
318          flx_q1, flx_h1, paprs, pplay, &
319          !!! nrlmd le 02/05/2011
320          Acoef_H_in, Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in, &
321          Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in, &
322          Kcoef_hq_in, gama_q_in, gama_h_in, &
323          !!!
324          flux_q, flux_h, d_q, d_t)
325
326    ! This routine calculates the flux and tendency of the specific humidity q and
327    ! the potential engergi H.
328    ! The quantities q and H are calculated according to
329    ! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
330    ! C and D are known from before and k is index of the vertical layer.
331
332    ! Input arguments
333    !****************************************************************************************
334    USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
335    USE lmdz_yomcst
336
337    IMPLICIT NONE
338    INTEGER, INTENT(IN) :: knon
339    REAL, INTENT(IN) :: dtime
340    REAL, DIMENSION(klon, klev), INTENT(IN) :: t_old, q_old
341    REAL, DIMENSION(klon), INTENT(IN) :: flx_q1, flx_h1
342    REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs
343    REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
344
345    !!! nrlmd le 02/05/2011
346    REAL, DIMENSION(klon), INTENT(IN) :: Acoef_H_in, Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in
347    REAL, DIMENSION(klon, klev), INTENT(IN) :: Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in
348    REAL, DIMENSION(klon, klev), INTENT(IN) :: Kcoef_hq_in, gama_q_in, gama_h_in
349    !!!
350
351    ! Output arguments
352    !****************************************************************************************
353    REAL, DIMENSION(klon, klev), INTENT(OUT) :: flux_q, flux_h, d_q, d_t
354
355    ! Local variables
356    !****************************************************************************************
357    LOGICAL, SAVE :: last = .FALSE.
358    REAL, DIMENSION(klon, klev) :: h_new, q_new
359    REAL, DIMENSION(klon) :: psref
360    INTEGER :: k, i, ierr
361
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.