source: LMDZ6/branches/Amaury_dev/libf/phylmd/climb_wind_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: 13.3 KB
Line 
1MODULE climb_wind_mod
2
3  ! Module to solve the verctical diffusion of the wind components "u" and "v".
4
5  USE dimphy
6  USE lmdz_abort_physic, ONLY: abort_physic
7
8  IMPLICIT NONE
9
10  SAVE
11  PRIVATE
12
13  REAL, DIMENSION(:), ALLOCATABLE :: alf1, alf2
14  !$OMP THREADPRIVATE(alf1,alf2)
15  REAL, DIMENSION(:, :), ALLOCATABLE :: Kcoefm
16  !$OMP THREADPRIVATE(Kcoefm)
17  REAL, DIMENSION(:, :), ALLOCATABLE :: Ccoef_U, Dcoef_U
18  !$OMP THREADPRIVATE(Ccoef_U, Dcoef_U)
19  REAL, DIMENSION(:, :), ALLOCATABLE :: Ccoef_V, Dcoef_V
20  !$OMP THREADPRIVATE(Ccoef_V, Dcoef_V)
21  REAL, DIMENSION(:), ALLOCATABLE :: Acoef_U, Bcoef_U
22  !$OMP THREADPRIVATE(Acoef_U, Bcoef_U)
23  REAL, DIMENSION(:), ALLOCATABLE :: Acoef_V, Bcoef_V
24  !$OMP THREADPRIVATE(Acoef_V, Bcoef_V)
25  LOGICAL :: firstcall = .TRUE.
26  !$OMP THREADPRIVATE(firstcall)
27
28  PUBLIC :: climb_wind_down, climb_wind_up
29
30CONTAINS
31
32  !****************************************************************************************
33
34  SUBROUTINE climb_wind_init
35
36    INTEGER :: ierr
37    CHARACTER(len = 20) :: modname = 'climb_wind_init'
38
39    !****************************************************************************************
40    ! Allocation of global module variables
41
42    !****************************************************************************************
43
44    ALLOCATE(alf1(klon), stat = ierr)
45    IF (ierr /= 0) CALL abort_physic(modname, 'Pb in allocate alf1', 1)
46
47    ALLOCATE(alf2(klon), stat = ierr)
48    IF (ierr /= 0) CALL abort_physic(modname, 'Pb in allocate alf2', 1)
49
50    ALLOCATE(Kcoefm(klon, klev), stat = ierr)
51    IF (ierr /= 0) CALL abort_physic(modname, 'Pb in allocate Kcoefm', 1)
52
53    ALLOCATE(Ccoef_U(klon, klev), stat = ierr)
54    IF (ierr /= 0) CALL abort_physic(modname, 'Pb in allocate Ccoef_U', 1)
55
56    ALLOCATE(Dcoef_U(klon, klev), stat = ierr)
57    IF (ierr /= 0) CALL abort_physic(modname, 'Pb in allocation Dcoef_U', 1)
58
59    ALLOCATE(Ccoef_V(klon, klev), stat = ierr)
60    IF (ierr /= 0) CALL abort_physic(modname, 'Pb in allocation Ccoef_V', 1)
61
62    ALLOCATE(Dcoef_V(klon, klev), stat = ierr)
63    IF (ierr /= 0) CALL abort_physic(modname, 'Pb in allocation Dcoef_V', 1)
64
65    ALLOCATE(Acoef_U(klon), Bcoef_U(klon), Acoef_V(klon), Bcoef_V(klon), STAT = ierr)
66    IF (ierr /= 0)  PRINT*, ' pb in allloc Acoef_U and Bcoef_U, ierr=', ierr
67
68    firstcall = .FALSE.
69
70  END SUBROUTINE climb_wind_init
71
72  !****************************************************************************************
73
74  SUBROUTINE climb_wind_down(knon, dtime, coef_in, pplay, paprs, temp, delp, u_old, v_old, &
75          !!! nrlmd le 02/05/2011
76          Ccoef_U_out, Ccoef_V_out, Dcoef_U_out, Dcoef_V_out, &
77          Kcoef_m_out, alf_1_out, alf_2_out, &
78          !!!
79          Acoef_U_out, Acoef_V_out, Bcoef_U_out, Bcoef_V_out)
80
81    ! This routine calculates for the wind components u and v,
82    ! recursivly the coefficients C and D in equation
83    ! X(k) = C(k) + D(k)*X(k-1), X=[u,v], k=[1,klev] is the vertical layer.
84
85
86    ! Input arguments
87    !****************************************************************************************
88    USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
89    USE lmdz_yomcst
90
91    IMPLICIT NONE
92    INTEGER, INTENT(IN) :: knon
93    REAL, INTENT(IN) :: dtime
94    REAL, DIMENSION(klon, klev), INTENT(IN) :: coef_in
95    REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay ! pres au milieu de couche (Pa)
96    REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs ! pression a inter-couche (Pa)
97    REAL, DIMENSION(klon, klev), INTENT(IN) :: temp  ! temperature
98    REAL, DIMENSION(klon, klev), INTENT(IN) :: delp
99    REAL, DIMENSION(klon, klev), INTENT(IN) :: u_old
100    REAL, DIMENSION(klon, klev), INTENT(IN) :: v_old
101
102    ! Output arguments
103    !****************************************************************************************
104    REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_U_out
105    REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_V_out
106    REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_U_out
107    REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_V_out
108
109    !!! nrlmd le 02/05/2011
110    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Ccoef_U_out
111    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Ccoef_V_out
112    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Dcoef_U_out
113    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Dcoef_V_out
114    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Kcoef_m_out
115    REAL, DIMENSION(klon), INTENT(OUT) :: alf_1_out
116    REAL, DIMENSION(klon), INTENT(OUT) :: alf_2_out
117    !!!
118
119    ! Local variables
120    !****************************************************************************************
121    REAL, DIMENSION(klon) :: u1lay, v1lay
122    INTEGER :: k, i
123
124    !****************************************************************************************
125    ! Initialize module
126    IF (firstcall) CALL climb_wind_init
127
128    !****************************************************************************************
129    ! Calculate the coefficients C and D in : u(k) = C(k) + D(k)*u(k-1)
130
131    !****************************************************************************************
132    ! - Define alpha (alf1 and alf2)
133    alf1(:) = 1.0
134    alf2(:) = 1.0 - alf1(:)
135
136    ! - Calculate the coefficients K
137    Kcoefm(:, :) = 0.0
138    DO k = 2, klev
139      DO i = 1, knon
140        Kcoefm(i, k) = coef_in(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      END DO
143    END DO
144
145    ! - Calculate the coefficients C and D, component "u"
146    CALL calc_coef(knon, Kcoefm(:, :), delp(:, :), &
147            u_old(:, :), alf1(:), alf2(:), &
148            Ccoef_U(:, :), Dcoef_U(:, :), Acoef_U(:), Bcoef_U(:))
149
150    ! - Calculate the coefficients C and D, component "v"
151    CALL calc_coef(knon, Kcoefm(:, :), delp(:, :), &
152            v_old(:, :), alf1(:), alf2(:), &
153            Ccoef_V(:, :), Dcoef_V(:, :), Acoef_V(:), Bcoef_V(:))
154
155    !****************************************************************************************
156    ! 6)
157    ! Return the first layer in output variables
158
159    !****************************************************************************************
160    Acoef_U_out = Acoef_U
161    Bcoef_U_out = Bcoef_U
162    Acoef_V_out = Acoef_V
163    Bcoef_V_out = Bcoef_V
164
165    !****************************************************************************************
166    ! 7)
167    ! If Pbl is split, return also the other layers in output variables
168
169    !****************************************************************************************
170    !!! jyg le 07/02/2012
171    !!jyg       IF (mod(iflag_pbl_split,2) .EQ.1) THEN
172    IF (mod(iflag_pbl_split, 10) >=1) THEN
173      !!! nrlmd le 02/05/2011
174      DO k = 1, klev
175        DO i = 1, klon
176          Ccoef_U_out(i, k) = Ccoef_U(i, k)
177          Ccoef_V_out(i, k) = Ccoef_V(i, k)
178          Dcoef_U_out(i, k) = Dcoef_U(i, k)
179          Dcoef_V_out(i, k) = Dcoef_V(i, k)
180          Kcoef_m_out(i, k) = Kcoefm(i, k)
181        ENDDO
182      ENDDO
183      DO i = 1, klon
184        alf_1_out(i) = alf1(i)
185        alf_2_out(i) = alf2(i)
186      ENDDO
187      !!!
188    ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
189    !!!
190
191  END SUBROUTINE climb_wind_down
192
193  !****************************************************************************************
194
195  SUBROUTINE calc_coef(knon, Kcoef, delp, X, alfa1, alfa2, Ccoef, Dcoef, Acoef, Bcoef)
196    USE lmdz_yomcst
197
198    IMPLICIT NONE
199
200    ! Find the coefficients C and D in fonction of alfa, K and delp
201
202    ! Input arguments
203    !****************************************************************************************
204    INTEGER, INTENT(IN) :: knon
205    REAL, DIMENSION(klon, klev), INTENT(IN) :: Kcoef, delp
206    REAL, DIMENSION(klon, klev), INTENT(IN) :: X
207    REAL, DIMENSION(klon), INTENT(IN) :: alfa1, alfa2
208
209    ! Output arguments
210    !****************************************************************************************
211    REAL, DIMENSION(klon), INTENT(OUT) :: Acoef, Bcoef
212    REAL, DIMENSION(klon, klev), INTENT(OUT) :: Ccoef, Dcoef
213
214    ! local variables
215    !****************************************************************************************
216    INTEGER :: k, i
217    REAL :: buf
218
219    !****************************************************************************************
220
221    ! Calculate coefficients C and D at top level, k=klev
222
223    Ccoef(:, :) = 0.0
224    Dcoef(:, :) = 0.0
225
226    DO i = 1, knon
227      buf = delp(i, klev) + Kcoef(i, klev)
228
229      Ccoef(i, klev) = X(i, klev) * delp(i, klev) / buf
230      Dcoef(i, klev) = Kcoef(i, klev) / buf
231    END DO
232
233    ! Calculate coefficients C and D at top level (klev-1) <= k <= 2
234
235    DO k = (klev - 1), 2, -1
236      DO i = 1, knon
237        buf = delp(i, k) + Kcoef(i, k) + Kcoef(i, k + 1) * (1. - Dcoef(i, k + 1))
238
239        Ccoef(i, k) = (X(i, k) * delp(i, k) + Kcoef(i, k + 1) * Ccoef(i, k + 1)) / buf
240        Dcoef(i, k) = Kcoef(i, k) / buf
241      END DO
242    END DO
243
244    ! Calculate coeffiecent A and B at surface
245
246    DO i = 1, knon
247      buf = delp(i, 1) + Kcoef(i, 2) * (1 - Dcoef(i, 2))
248      Acoef(i) = (X(i, 1) * delp(i, 1) + Kcoef(i, 2) * Ccoef(i, 2)) / buf
249      Bcoef(i) = -RG / buf
250    END DO
251
252  END SUBROUTINE calc_coef
253
254  !****************************************************************************************
255
256  SUBROUTINE climb_wind_up(knon, dtime, u_old, v_old, flx_u1, flx_v1, &
257          !!! nrlmd le 02/05/2011
258          Acoef_U_in, Acoef_V_in, Bcoef_U_in, Bcoef_V_in, &
259          Ccoef_U_in, Ccoef_V_in, Dcoef_U_in, Dcoef_V_in, &
260          Kcoef_m_in, &
261          !!!
262          flx_u_new, flx_v_new, d_u_new, d_v_new)
263
264    ! Diffuse the wind components from the surface layer and up to the top layer.
265    ! Coefficents A, B, C and D are known from before. Start values for the diffusion are the
266    ! momentum fluxes at surface.
267
268    ! u(k=1) = A + B*flx*dtime
269    ! u(k)   = C(k) + D(k)*u(k-1)  [2 <= k <= klev]
270
271    !****************************************************************************************
272
273    ! Input arguments
274    !****************************************************************************************
275    USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
276    USE lmdz_yomcst
277
278    IMPLICIT NONE
279
280    INTEGER, INTENT(IN) :: knon
281    REAL, INTENT(IN) :: dtime
282    REAL, DIMENSION(klon, klev), INTENT(IN) :: u_old
283    REAL, DIMENSION(klon, klev), INTENT(IN) :: v_old
284    REAL, DIMENSION(klon), INTENT(IN) :: flx_u1, flx_v1 ! momentum flux
285
286    !!! nrlmd le 02/05/2011
287    REAL, DIMENSION(klon), INTENT(IN) :: Acoef_U_in, Acoef_V_in, Bcoef_U_in, Bcoef_V_in
288    REAL, DIMENSION(klon, klev), INTENT(IN) :: Ccoef_U_in, Ccoef_V_in, Dcoef_U_in, Dcoef_V_in
289    REAL, DIMENSION(klon, klev), INTENT(IN) :: Kcoef_m_in
290    !!!
291
292    ! Output arguments
293    !****************************************************************************************
294    REAL, DIMENSION(klon, klev), INTENT(OUT) :: flx_u_new, flx_v_new
295    REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_u_new, d_v_new
296
297    ! Local variables
298    !****************************************************************************************
299    REAL, DIMENSION(klon, klev) :: u_new, v_new
300    INTEGER :: k, i
301    !****************************************************************************************
302
303    !!! jyg le 07/02/2012
304    !!jyg       IF (mod(iflag_pbl_split,2) .EQ.1) THEN
305    IF (mod(iflag_pbl_split, 10) >=1) THEN
306      !!! nrlmd le 02/05/2011
307      DO i = 1, knon
308        Acoef_U(i) = Acoef_U_in(i)
309        Acoef_V(i) = Acoef_V_in(i)
310        Bcoef_U(i) = Bcoef_U_in(i)
311        Bcoef_V(i) = Bcoef_V_in(i)
312      ENDDO
313      DO k = 1, klev
314        DO i = 1, knon
315          Ccoef_U(i, k) = Ccoef_U_in(i, k)
316          Ccoef_V(i, k) = Ccoef_V_in(i, k)
317          Dcoef_U(i, k) = Dcoef_U_in(i, k)
318          Dcoef_V(i, k) = Dcoef_V_in(i, k)
319          Kcoefm(i, k) = Kcoef_m_in(i, k)
320        ENDDO
321      ENDDO
322      !!!
323    ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
324    !!!
325
326    ! Niveau 1
327    DO i = 1, knon
328      u_new(i, 1) = Acoef_U(i) + Bcoef_U(i) * flx_u1(i) * dtime
329      v_new(i, 1) = Acoef_V(i) + Bcoef_V(i) * flx_v1(i) * dtime
330    END DO
331
332    ! Niveau 2 jusqu'au sommet klev
333    DO k = 2, klev
334      DO i = 1, knon
335        u_new(i, k) = Ccoef_U(i, k) + Dcoef_U(i, k) * u_new(i, k - 1)
336        v_new(i, k) = Ccoef_V(i, k) + Dcoef_V(i, k) * v_new(i, k - 1)
337      END DO
338    END DO
339
340    !****************************************************************************************
341    ! Calcul flux
342
343    !== flux_u/v est le flux de moment angulaire (positif vers bas)
344    !== dont l'unite est: (kg m/s)/(m**2 s)
345
346    !****************************************************************************************
347
348    flx_u_new(:, :) = 0.0
349    flx_v_new(:, :) = 0.0
350
351    flx_u_new(1:knon, 1) = flx_u1(1:knon)
352    flx_v_new(1:knon, 1) = flx_v1(1:knon)
353
354    ! Niveau 2->klev
355    DO k = 2, klev
356      DO i = 1, knon
357        flx_u_new(i, k) = Kcoefm(i, k) / RG / dtime * &
358                (u_new(i, k) - u_new(i, k - 1))
359
360        flx_v_new(i, k) = Kcoefm(i, k) / RG / dtime * &
361                (v_new(i, k) - v_new(i, k - 1))
362      END DO
363    END DO
364
365    !****************************************************************************************
366    ! Calcul tendances
367
368    !****************************************************************************************
369    d_u_new(:, :) = 0.0
370    d_v_new(:, :) = 0.0
371    DO k = 1, klev
372      DO i = 1, knon
373        d_u_new(i, k) = u_new(i, k) - u_old(i, k)
374        d_v_new(i, k) = v_new(i, k) - v_old(i, k)
375      END DO
376    END DO
377
378  END SUBROUTINE climb_wind_up
379
380  !****************************************************************************************
381
382END MODULE climb_wind_mod
Note: See TracBrowser for help on using the repository browser.