source: LMDZ6/branches/Amaury_dev/libf/phylmd/climb_wind_mod.F90 @ 5139

Last change on this file since 5139 was 5139, checked in by abarral, 8 weeks ago

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