source: LMDZ5/trunk/libf/phylmd/climb_wind_mod.F90 @ 1992

Last change on this file since 1992 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

  • 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:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 10.6 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
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 
29  PUBLIC :: climb_wind_down, climb_wind_up
30
31CONTAINS
32!
33!****************************************************************************************
34!
35  SUBROUTINE climb_wind_init
36
37    INTEGER             :: ierr
38    CHARACTER(len = 20) :: modname = 'climb_wind_init'   
39
40!****************************************************************************************
41! Allocation of global module variables
42!
43!****************************************************************************************
44
45    ALLOCATE(alf1(klon), stat=ierr)
46    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocate alf2',1)
47
48    ALLOCATE(alf2(klon), stat=ierr)
49    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocate alf2',1)
50
51    ALLOCATE(Kcoefm(klon,klev), stat=ierr)
52    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocate Kcoefm',1)
53
54    ALLOCATE(Ccoef_U(klon,klev), stat=ierr)
55    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocate Ccoef_U',1)
56
57    ALLOCATE(Dcoef_U(klon,klev), stat=ierr)
58    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocation Dcoef_U',1)
59
60    ALLOCATE(Ccoef_V(klon,klev), stat=ierr)
61    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocation Ccoef_V',1)
62
63    ALLOCATE(Dcoef_V(klon,klev), stat=ierr)
64    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocation Dcoef_V',1)
65
66    ALLOCATE(Acoef_U(klon), Bcoef_U(klon), Acoef_V(klon), Bcoef_V(klon), STAT=ierr)
67    IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_U and Bcoef_U, ierr=', ierr
68
69    firstcall=.FALSE.
70
71  END SUBROUTINE climb_wind_init
72!
73!****************************************************************************************
74!
75  SUBROUTINE climb_wind_down(knon, dtime, coef_in, pplay, paprs, temp, delp, u_old, v_old, &
76       Acoef_U_out, Acoef_V_out, Bcoef_U_out, Bcoef_V_out)
77!
78! This routine calculates for the wind components u and v,
79! recursivly the coefficients C and D in equation
80! X(k) = C(k) + D(k)*X(k-1), X=[u,v], k=[1,klev] is the vertical layer.
81!
82!
83    INCLUDE "YOMCST.h"
84! Input arguments
85!****************************************************************************************
86    INTEGER, INTENT(IN)                      :: knon
87    REAL, INTENT(IN)                         :: dtime
88    REAL, DIMENSION(klon,klev), INTENT(IN)   :: coef_in
89    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay ! pres au milieu de couche (Pa)
90    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression a inter-couche (Pa)
91    REAL, DIMENSION(klon,klev), INTENT(IN)   :: temp  ! temperature
92    REAL, DIMENSION(klon,klev), INTENT(IN)   :: delp
93    REAL, DIMENSION(klon,klev), INTENT(IN)   :: u_old
94    REAL, DIMENSION(klon,klev), INTENT(IN)   :: v_old
95
96! Output arguments
97!****************************************************************************************
98    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_U_out
99    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_V_out
100    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_U_out
101    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_V_out
102
103! Local variables
104!****************************************************************************************
105    REAL, DIMENSION(klon)                    :: u1lay, v1lay
106    INTEGER                                  :: k, i
107
108
109!****************************************************************************************
110! Initialize module
111    IF (firstcall) CALL climb_wind_init
112
113!****************************************************************************************
114! Calculate the coefficients C and D in : u(k) = C(k) + D(k)*u(k-1)
115!
116!****************************************************************************************
117! - Define alpha (alf1 and alf2)
118    alf1(:) = 1.0
119    alf2(:) = 1.0 - alf1(:)
120
121! - Calculate the coefficients K
122    Kcoefm(:,:) = 0.0
123    DO k = 2, klev
124       DO i=1,knon
125          Kcoefm(i,k) = coef_in(i,k)*RG*RG*dtime/(pplay(i,k-1)-pplay(i,k)) &
126               *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
127       END DO
128    END DO
129
130! - Calculate the coefficients C and D, component "u"
131    CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), &
132         u_old(:,:), alf1(:), alf2(:),  &
133         Ccoef_U(:,:), Dcoef_U(:,:), Acoef_U(:), Bcoef_U(:))
134
135! - Calculate the coefficients C and D, component "v"
136    CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), &
137         v_old(:,:), alf1(:), alf2(:),  &
138         Ccoef_V(:,:), Dcoef_V(:,:), Acoef_V(:), Bcoef_V(:))
139
140!****************************************************************************************
141! 6)
142! Return the first layer in output variables
143!
144!****************************************************************************************
145    Acoef_U_out = Acoef_U
146    Bcoef_U_out = Bcoef_U
147    Acoef_V_out = Acoef_V
148    Bcoef_V_out = Bcoef_V
149
150  END SUBROUTINE climb_wind_down
151!
152!****************************************************************************************
153!
154  SUBROUTINE calc_coef(knon, Kcoef, delp, X, alfa1, alfa2, Ccoef, Dcoef, Acoef, Bcoef)
155!
156! Find the coefficients C and D in fonction of alfa, K and delp
157!
158! Input arguments
159!****************************************************************************************
160    INTEGER, INTENT(IN)                      :: knon
161    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, delp
162    REAL, DIMENSION(klon,klev), INTENT(IN)   :: X
163    REAL, DIMENSION(klon), INTENT(IN)        :: alfa1, alfa2
164
165! Output arguments
166!****************************************************************************************
167    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef, Bcoef
168    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef, Dcoef
169 
170! local variables
171!****************************************************************************************
172    INTEGER                                  :: k, i
173    REAL                                     :: buf
174
175    INCLUDE "YOMCST.h"
176!****************************************************************************************
177!
178
179! Calculate coefficients C and D at top level, k=klev
180!
181    Ccoef(:,:) = 0.0
182    Dcoef(:,:) = 0.0
183
184    DO i = 1, knon
185       buf = delp(i,klev) + Kcoef(i,klev)
186
187       Ccoef(i,klev) = X(i,klev)*delp(i,klev)/buf
188       Dcoef(i,klev) = Kcoef(i,klev)/buf
189    END DO
190   
191!
192! Calculate coefficients C and D at top level (klev-1) <= k <= 2
193!
194    DO k=(klev-1),2,-1
195       DO i = 1, knon
196          buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
197         
198          Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1))/buf
199          Dcoef(i,k) = Kcoef(i,k)/buf
200       END DO
201    END DO
202
203!
204! Calculate coeffiecent A and B at surface
205!
206    DO i = 1, knon
207       buf = delp(i,1) + Kcoef(i,2)*(1-Dcoef(i,2))
208       Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*Ccoef(i,2))/buf
209       Bcoef(i) = -RG/buf
210    END DO
211    acoef(knon+1: klon) = 0.
212    bcoef(knon+1: klon) = 0.
213
214  END SUBROUTINE calc_coef
215!
216!****************************************************************************************
217!
218
219  SUBROUTINE climb_wind_up(knon, dtime, u_old, v_old, flx_u1, flx_v1,  &
220       flx_u_new, flx_v_new, d_u_new, d_v_new)
221!
222! Diffuse the wind components from the surface layer and up to the top layer.
223! Coefficents A, B, C and D are known from before. Start values for the diffusion are the
224! momentum fluxes at surface.
225!
226! u(k=1) = A + B*flx*dtime
227! u(k)   = C(k) + D(k)*u(k-1)  [2 <= k <= klev]
228!
229!****************************************************************************************
230    INCLUDE "YOMCST.h"
231
232! Input arguments
233!****************************************************************************************
234    INTEGER, INTENT(IN)                     :: knon
235    REAL, INTENT(IN)                        :: dtime
236    REAL, DIMENSION(klon,klev), INTENT(IN)  :: u_old
237    REAL, DIMENSION(klon,klev), INTENT(IN)  :: v_old
238    REAL, DIMENSION(klon), INTENT(IN)       :: flx_u1, flx_v1 ! momentum flux
239
240! Output arguments
241!****************************************************************************************
242    REAL, DIMENSION(klon,klev), INTENT(OUT) :: flx_u_new, flx_v_new
243    REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_u_new, d_v_new
244
245! Local variables
246!****************************************************************************************
247    REAL, DIMENSION(klon,klev)              :: u_new, v_new
248    INTEGER                                 :: k, i
249   
250!
251!****************************************************************************************
252
253! Niveau 1
254    DO i = 1, knon
255       u_new(i,1) = Acoef_U(i) + Bcoef_U(i)*flx_u1(i)*dtime
256       v_new(i,1) = Acoef_V(i) + Bcoef_V(i)*flx_v1(i)*dtime
257    END DO
258
259! Niveau 2 jusqu'au sommet klev
260    DO k = 2, klev
261       DO i=1, knon
262          u_new(i,k) = Ccoef_U(i,k) + Dcoef_U(i,k) * u_new(i,k-1)
263          v_new(i,k) = Ccoef_V(i,k) + Dcoef_V(i,k) * v_new(i,k-1)
264       END DO
265    END DO
266
267!****************************************************************************************
268! Calcul flux
269!
270!== flux_u/v est le flux de moment angulaire (positif vers bas)
271!== dont l'unite est: (kg m/s)/(m**2 s)
272!
273!****************************************************************************************
274!
275    flx_u_new(:,:) = 0.0
276    flx_v_new(:,:) = 0.0
277
278    flx_u_new(1:knon,1)=flx_u1(1:knon)
279    flx_v_new(1:knon,1)=flx_v1(1:knon)
280
281! Niveau 2->klev
282    DO k = 2, klev
283       DO i = 1, knon
284          flx_u_new(i,k) = Kcoefm(i,k)/RG/dtime * &
285               (u_new(i,k)-u_new(i,k-1))
286         
287          flx_v_new(i,k) = Kcoefm(i,k)/RG/dtime * &
288               (v_new(i,k)-v_new(i,k-1))
289       END DO
290    END DO
291
292!****************************************************************************************
293! Calcul tendances
294!
295!****************************************************************************************
296    d_u_new(:,:) = 0.0
297    d_v_new(:,:) = 0.0
298    DO k = 1, klev
299       DO i = 1, knon
300          d_u_new(i,k) = u_new(i,k) - u_old(i,k)
301          d_v_new(i,k) = v_new(i,k) - v_old(i,k)
302       END DO
303    END DO
304
305  END SUBROUTINE climb_wind_up
306!
307!****************************************************************************************
308!
309END MODULE climb_wind_mod
Note: See TracBrowser for help on using the repository browser.