source: LMDZ6/branches/Amaury_dev/libf/phylmd/climb_qbs_mod.F90 @ 5877

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

Put YOMCST.h into modules
File size: 13.6 KB
Line 
1MODULE climb_qbs_mod
2
3  ! Module to solve the verctical diffusion of blowing snow;
4
5  USE dimphy
6
7  IMPLICIT NONE
8  SAVE
9  PRIVATE
10  PUBLIC :: climb_qbs_down, climb_qbs_up
11
12  REAL, DIMENSION(:, :), ALLOCATABLE :: gamaqbs
13  !$OMP THREADPRIVATE(gamaqbs)
14  REAL, DIMENSION(:, :), ALLOCATABLE :: Ccoef_QBS, Dcoef_QBS
15  !$OMP THREADPRIVATE(Ccoef_QBS, Dcoef_QBS)
16  REAL, DIMENSION(:), ALLOCATABLE :: Acoef_QBS, Bcoef_QBS
17  !$OMP THREADPRIVATE(Acoef_QBS, Bcoef_QBS)
18  REAL, DIMENSION(:, :), ALLOCATABLE :: Kcoefqbs
19  !$OMP THREADPRIVATE(Kcoefqbs)
20
21CONTAINS
22
23  !****************************************************************************************
24
25  SUBROUTINE climb_qbs_down(knon, coefqbs, paprs, pplay, &
26          delp, temp, qbs, dtime, &
27          Ccoef_QBS_out, Dcoef_QBS_out, &
28          Kcoef_qbs_out, gama_qbs_out, &
29          Acoef_QBS_out, Bcoef_QBS_out)
30
31    ! This routine calculates recursivly the coefficients C and D
32    ! for the quantity X=[QBS] in equation X(k) = C(k) + D(k)*X(k-1), where k is
33    ! the index of the vertical layer.
34
35    ! Input arguments
36    !****************************************************************************************
37    USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
38    USE lmdz_yomcst
39
40    IMPLICIT NONE
41    INTEGER, INTENT(IN) :: knon
42    REAL, DIMENSION(klon, klev), INTENT(IN) :: coefqbs
43    REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
44    REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs
45    REAL, DIMENSION(klon, klev), INTENT(IN) :: delp
46    REAL, DIMENSION(klon, klev), INTENT(IN) :: temp
47    REAL, DIMENSION(klon, klev), INTENT(IN) :: qbs
48    REAL, INTENT(IN) :: dtime
49
50    ! Output arguments
51    !****************************************************************************************
52    REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_QBS_out
53    REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_QBS_out
54
55    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Ccoef_QBS_out
56    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Dcoef_QBS_out
57    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Kcoef_qbs_out
58    REAL, DIMENSION(klon, klev), INTENT(OUT) :: gama_qbs_out
59
60    ! Local variables
61    !****************************************************************************************
62    LOGICAL, SAVE :: first = .TRUE.
63    !$OMP THREADPRIVATE(first)
64    REAL, DIMENSION(klon) :: psref
65    REAL :: delz, pkh
66    INTEGER :: k, i, ierr
67    !****************************************************************************************
68    ! 1)
69    ! Allocation at first time step only
70
71    !****************************************************************************************
72
73    IF (first) THEN
74      first = .FALSE.
75      ALLOCATE(Ccoef_QBS(klon, klev), STAT = ierr)
76      IF (ierr /= 0)  PRINT*, ' pb in allloc Ccoef_QBS, ierr=', ierr
77
78      ALLOCATE(Dcoef_QBS(klon, klev), STAT = ierr)
79      IF (ierr /= 0)  PRINT*, ' pb in allloc Dcoef_QBS, ierr=', ierr
80
81      ALLOCATE(Acoef_QBS(klon), Bcoef_QBS(klon), STAT = ierr)
82      IF (ierr /= 0)  PRINT*, ' pb in allloc Acoef_BS and Bcoef_BS, ierr=', ierr
83
84      ALLOCATE(Kcoefqbs(klon, klev), STAT = ierr)
85      IF (ierr /= 0)  PRINT*, ' pb in allloc Kcoefqbs, ierr=', ierr
86
87      ALLOCATE(gamaqbs(1:klon, 2:klev), STAT = ierr)
88      IF (ierr /= 0) PRINT*, ' pb in allloc gamaqbs, ierr=', ierr
89
90    END IF
91
92    !****************************************************************************************
93    ! 2)
94    ! Definition of the coeficient K
95
96    !****************************************************************************************
97    Kcoefqbs(:, :) = 0.0
98    DO k = 2, klev
99      DO i = 1, knon
100        Kcoefqbs(i, k) = &
101                coefqbs(i, k) * RG * RG * dtime / (pplay(i, k - 1) - pplay(i, k)) &
102                        * (paprs(i, k) * 2 / (temp(i, k) + temp(i, k - 1)) / RD)**2
103      ENDDO
104    ENDDO
105
106    !****************************************************************************************
107    ! 3)
108    ! Calculation of gama for "Q" and "H"
109
110    !****************************************************************************************
111    !   surface pressure is used as reference
112    psref(:) = paprs(:, 1)
113
114    !   definition of gama
115    IF (iflag_pbl == 1) THEN
116      gamaqbs(:, :) = 0.0
117
118      ! conversion de gama
119      DO k = 2, klev
120        DO i = 1, knon
121          delz = RD * (temp(i, k - 1) + temp(i, k)) / &
122                  2.0 / RG / paprs(i, k) * (pplay(i, k - 1) - pplay(i, k))
123          pkh = (psref(i) / paprs(i, k))**RKAPPA
124
125          ! convertie gradient verticale de contenu en neige soufflee en difference de neige soufflee entre centre de couches
126          gamaqbs(i, k) = gamaqbs(i, k) * delz
127        ENDDO
128      ENDDO
129
130    ELSE
131      gamaqbs(:, :) = 0.0
132    ENDIF
133
134
135    !****************************************************************************************
136    ! 4)
137    ! Calculte the coefficients C and D for specific content of blowing snow, qbs
138
139    !****************************************************************************************
140
141    CALL calc_coef_qbs(knon, Kcoefqbs(:, :), gamaqbs(:, :), delp(:, :), qbs(:, :), &
142            Ccoef_QBS(:, :), Dcoef_QBS(:, :), Acoef_QBS, Bcoef_QBS)
143
144
145    !****************************************************************************************
146    ! 5)
147    ! Return the first layer in output variables
148
149    !****************************************************************************************
150    Acoef_QBS_out = Acoef_QBS
151    Bcoef_QBS_out = Bcoef_QBS
152
153    !****************************************************************************************
154    ! 6)
155    ! If Pbl is split, return also the other layers in output variables
156
157    !****************************************************************************************
158    IF (mod(iflag_pbl_split, 10) >=1) THEN
159      DO k = 1, klev
160        DO i = 1, klon
161          Ccoef_QBS_out(i, k) = Ccoef_QBS(i, k)
162          Dcoef_QBS_out(i, k) = Dcoef_QBS(i, k)
163          Kcoef_qbs_out(i, k) = Kcoefqbs(i, k)
164          IF (k==1) THEN
165            gama_qbs_out(i, k) = 0.
166          ELSE
167            gama_qbs_out(i, k) = gamaqbs(i, k)
168          ENDIF
169        ENDDO
170      ENDDO
171      !!!
172    ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
173    !!!
174
175  END SUBROUTINE climb_qbs_down
176
177  !****************************************************************************************
178
179  SUBROUTINE calc_coef_qbs(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
180
181    ! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
182    ! where X is QQBS, and k the vertical level k=1,klev
183    USE lmdz_yomcst
184    IMPLICIT NONE
185    ! Input arguments
186    !****************************************************************************************
187    INTEGER, INTENT(IN) :: knon
188    REAL, DIMENSION(klon, klev), INTENT(IN) :: Kcoef, delp
189    REAL, DIMENSION(klon, klev), INTENT(IN) :: X
190    REAL, DIMENSION(klon, 2:klev), INTENT(IN) :: gama
191
192    ! Output arguments
193    !****************************************************************************************
194    REAL, DIMENSION(klon), INTENT(OUT) :: Acoef, Bcoef
195    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Ccoef, Dcoef
196
197    ! Local variables
198    !****************************************************************************************
199    INTEGER :: k, i
200    REAL :: buf
201
202    !****************************************************************************************
203    ! Niveau au sommet, k=klev
204
205    !****************************************************************************************
206    Ccoef(:, :) = 0.0
207    Dcoef(:, :) = 0.0
208
209    DO i = 1, knon
210      buf = delp(i, klev) + Kcoef(i, klev)
211
212      Ccoef(i, klev) = (X(i, klev) * delp(i, klev) - Kcoef(i, klev) * gama(i, klev)) / buf
213      Dcoef(i, klev) = Kcoef(i, klev) / buf
214    END DO
215
216
217    !****************************************************************************************
218    ! Niveau  (klev-1) <= k <= 2
219
220    !****************************************************************************************
221
222    DO k = (klev - 1), 2, -1
223      DO i = 1, knon
224        buf = delp(i, k) + Kcoef(i, k) + Kcoef(i, k + 1) * (1. - Dcoef(i, k + 1))
225        Ccoef(i, k) = (X(i, k) * delp(i, k) + Kcoef(i, k + 1) * Ccoef(i, k + 1) + &
226                Kcoef(i, k + 1) * gama(i, k + 1) - Kcoef(i, k) * gama(i, k)) / buf
227        Dcoef(i, k) = Kcoef(i, k) / buf
228      END DO
229    END DO
230
231    !****************************************************************************************
232    ! Niveau k=1
233
234    !****************************************************************************************
235
236    DO i = 1, knon
237      buf = delp(i, 1) + Kcoef(i, 2) * (1. - Dcoef(i, 2))
238      Acoef(i) = (X(i, 1) * delp(i, 1) + Kcoef(i, 2) * (gama(i, 2) + Ccoef(i, 2))) / buf
239      Bcoef(i) = -1. * RG / buf
240    END DO
241
242  END SUBROUTINE calc_coef_qbs
243
244  !****************************************************************************************
245
246  SUBROUTINE climb_qbs_up(knon, dtime, qbs_old, &
247          flx_qbs1, paprs, pplay, &
248          Acoef_QBS_in, Bcoef_QBS_in, &
249          Ccoef_QBS_in, Dcoef_QBS_in, &
250          Kcoef_qbs_in, gama_qbs_in, &
251          flux_qbs, d_qbs)
252
253    ! This routine calculates the flux and tendency of the specific content of blowing snow qbs
254    ! The quantity qbs is calculated according to
255    ! X(k) = C(k) + D(k)*X(k-1) for X=[qbs], where the coefficients
256    ! C and D are known from before and k is index of the vertical layer.
257
258    ! Input arguments
259    !****************************************************************************************
260    USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
261    USE lmdz_yomcst
262
263    IMPLICIT NONE
264
265    INTEGER, INTENT(IN) :: knon
266    REAL, INTENT(IN) :: dtime
267    REAL, DIMENSION(klon, klev), INTENT(IN) :: qbs_old
268    REAL, DIMENSION(klon), INTENT(IN) :: flx_qbs1
269    REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs
270    REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay
271
272    !!! nrlmd le 02/05/2011
273    REAL, DIMENSION(klon), INTENT(IN) :: Acoef_QBS_in, Bcoef_QBS_in
274    REAL, DIMENSION(klon, klev), INTENT(IN) :: Ccoef_QBS_in, Dcoef_QBS_in
275    REAL, DIMENSION(klon, klev), INTENT(IN) :: Kcoef_qbs_in, gama_qbs_in
276    !!!
277
278    ! Output arguments
279    !****************************************************************************************
280    REAL, DIMENSION(klon, klev), INTENT(OUT) :: flux_qbs, d_qbs
281
282    ! Local variables
283    !****************************************************************************************
284    LOGICAL, SAVE :: last = .FALSE.
285    REAL, DIMENSION(klon, klev) :: qbs_new
286    REAL, DIMENSION(klon) :: psref
287    INTEGER :: k, i, ierr
288
289    !****************************************************************************************
290    ! 1)
291    ! Definition of some variables
292    REAL, DIMENSION(klon, klev) :: zairm
293
294    !****************************************************************************************
295    flux_qbs(:, :) = 0.0
296    d_qbs(:, :) = 0.0
297
298    psref(1:knon) = paprs(1:knon, 1)
299
300    IF (mod(iflag_pbl_split, 10) >=1) THEN
301      DO i = 1, knon
302        Acoef_QBS(i) = Acoef_QBS_in(i)
303        Bcoef_QBS(i) = Bcoef_QBS_in(i)
304      ENDDO
305      DO k = 1, klev
306        DO i = 1, knon
307          Ccoef_QBS(i, k) = Ccoef_QBS_in(i, k)
308          Dcoef_QBS(i, k) = Dcoef_QBS_in(i, k)
309          Kcoefqbs(i, k) = Kcoef_qbs_in(i, k)
310          IF (k>1) THEN
311            gamaqbs(i, k) = gama_qbs_in(i, k)
312          ENDIF
313        ENDDO
314      ENDDO
315      !!!
316    ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
317    !!!
318
319    !****************************************************************************************
320    ! 2)
321    ! Calculation of QBS
322
323    !****************************************************************************************
324
325    !- First layer
326    qbs_new(1:knon, 1) = Acoef_QBS(1:knon) + Bcoef_QBS(1:knon) * flx_qbs1(1:knon) * dtime
327    !- All the other layers
328    DO k = 2, klev
329      DO i = 1, knon
330        qbs_new(i, k) = Ccoef_QBS(i, k) + Dcoef_QBS(i, k) * qbs_new(i, k - 1)
331      END DO
332    END DO
333    !****************************************************************************************
334    ! 3)
335    ! Calculation of the flux for QBS
336
337    !****************************************************************************************
338
339    !- The flux at first layer, k=1
340    flux_qbs(1:knon, 1) = flx_qbs1(1:knon)
341
342    !- The flux at all layers above surface
343    DO k = 2, klev
344      DO i = 1, knon
345        flux_qbs(i, k) = (Kcoefqbs(i, k) / RG / dtime) * &
346                (qbs_new(i, k) - qbs_new(i, k - 1) + gamaqbs(i, k))
347      END DO
348    END DO
349
350    !****************************************************************************************
351    ! 4)
352    ! Calculation of tendency for QBS
353
354    !****************************************************************************************
355    DO k = 1, klev
356      DO i = 1, knon
357        d_qbs(i, k) = qbs_new(i, k) - qbs_old(i, k)
358        zairm(i, k) = (paprs(i, k) - paprs(i, k + 1)) / rg
359      END DO
360    END DO
361
362    !****************************************************************************************
363    ! Some deallocations
364
365    !****************************************************************************************
366    IF (last) THEN
367      DEALLOCATE(Ccoef_QBS, Dcoef_QBS, stat = ierr)
368      IF (ierr /= 0)  PRINT*, ' pb in dealllocate Ccoef_QBS, Dcoef_QBS, ierr=', ierr
369      DEALLOCATE(Acoef_QBS, Bcoef_QBS, stat = ierr)
370      IF (ierr /= 0)  PRINT*, ' pb in dealllocate Acoef_QBS, Bcoef_QBS, ierr=', ierr
371      DEALLOCATE(gamaqbs, stat = ierr)
372      IF (ierr /= 0)  PRINT*, ' pb in dealllocate gamaqbs, ierr=', ierr
373      DEALLOCATE(Kcoefqbs, stat = ierr)
374      IF (ierr /= 0)  PRINT*, ' pb in dealllocate Kcoefqbs, ierr=', ierr
375    END IF
376  END SUBROUTINE climb_qbs_up
377
378  !****************************************************************************************
379
380END MODULE climb_qbs_mod
381
382 
383
384
385
386
Note: See TracBrowser for help on using the repository browser.