source: LMDZ4/trunk/libf/phy_IPCC_AR4/climb_wind_mod.F90 @ 1068

Last change on this file since 1068 was 868, checked in by Laurent Fairhead, 17 years ago

Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4
par la version de la physique avec thermique. On garde le repertoire phylmd
pour un petit moment pour que les utilisateurs ne soient pas trop perdus ...
phy_IPCC_AR4 = phylmd
LF
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 12.2 KB
Line 
1!
2! $Header$
3!
4MODULE climb_wind_mod
5!
6! Module to solve the verctical diffusion of the wind components "u" and "v".
7!
8  USE dimphy
9
10  IMPLICIT NONE
11
12  SAVE
13  PRIVATE
14 
15  REAL, DIMENSION(:),   ALLOCATABLE  :: alf1, alf2
16  !$OMP THREADPRIVATE(alf1,alf2)
17  REAL, DIMENSION(:,:), ALLOCATABLE  :: Kcoefm
18  !$OMP THREADPRIVATE(Kcoefm)
19  REAL, DIMENSION(:,:), ALLOCATABLE  :: Ccoef_U, Dcoef_U
20  !$OMP THREADPRIVATE(Ccoef_U, Dcoef_U)
21  REAL, DIMENSION(:,:), ALLOCATABLE  :: Ccoef_V, Dcoef_V
22  !$OMP THREADPRIVATE(Ccoef_V, Dcoef_V)
23  LOGICAL                            :: firstcall=.TRUE.
24  !$OMP THREADPRIVATE(firstcall)
25
26 
27  PUBLIC :: climb_wind_down, calcul_wind_flux, climb_wind_up
28
29CONTAINS
30!
31!****************************************************************************************
32!
33  SUBROUTINE climb_wind_init
34
35    INTEGER             :: ierr
36    CHARACTER(len = 20) :: modname = 'climb_wind_init'   
37
38!****************************************************************************************
39! Allocation of global module variables
40!
41!****************************************************************************************
42
43    ALLOCATE(alf1(klon), stat=ierr)
44    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocate alf2',1)
45
46    ALLOCATE(alf2(klon), stat=ierr)
47    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocate alf2',1)
48
49    ALLOCATE(Kcoefm(klon,klev), stat=ierr)
50    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocate Kcoefm',1)
51
52    ALLOCATE(Ccoef_U(klon,klev), stat=ierr)
53    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocate Ccoef_U',1)
54
55    ALLOCATE(Dcoef_U(klon,klev), stat=ierr)
56    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocation Dcoef_U',1)
57
58    ALLOCATE(Ccoef_V(klon,klev), stat=ierr)
59    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocation Ccoef_V',1)
60
61    ALLOCATE(Dcoef_V(klon,klev), stat=ierr)
62    IF (ierr /= 0) CALL abort_gcm(modname,'Pb in allocation Dcoef_V',1)
63
64    firstcall=.FALSE.
65
66  END SUBROUTINE climb_wind_init
67!
68!****************************************************************************************
69!
70  SUBROUTINE climb_wind_down(knon, dtime, coef_in, pplay, paprs, temp, delp, u_old, v_old)
71!
72! This routine calculates for the wind components u and v,
73! recursivly the coefficients C and D in equation
74! X(k) = C(k) + D(k)*X(k-1), X=[u,v], k=[1,klev] is the vertical layer.
75!
76!
77    INCLUDE "YOMCST.h"
78! Input arguments
79!****************************************************************************************
80    INTEGER, INTENT(IN)                      :: knon
81    REAL, INTENT(IN)                         :: dtime
82    REAL, DIMENSION(klon,klev), INTENT(IN)   :: coef_in
83    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay ! pres au milieu de couche (Pa)
84    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression a inter-couche (Pa)
85    REAL, DIMENSION(klon,klev), INTENT(IN)   :: temp  ! temperature
86    REAL, DIMENSION(klon,klev), INTENT(IN)   :: delp
87    REAL, DIMENSION(klon,klev), INTENT(IN)   :: u_old
88    REAL, DIMENSION(klon,klev), INTENT(IN)   :: v_old
89
90! Local variables
91!****************************************************************************************
92    REAL, DIMENSION(klon)                    :: u1lay, v1lay
93    INTEGER                                  :: k, i
94
95
96!****************************************************************************************
97! Initialize module
98    IF (firstcall) CALL climb_wind_init
99
100!****************************************************************************************
101! Calculate the coefficients C and D in : u(k) = C(k) + D(k)*u(k-1)
102!
103!****************************************************************************************
104
105! - Define alpha (alf1 and alf2)
106    alf1(:) = 1.0
107    alf2(:) = 1.0 - alf1(:)
108
109! - Calculte the wind components for the first layer
110    u1lay(1:knon) = u_old(1:knon,1)*alf1(1:knon) + u_old(1:knon,2)*alf2(1:knon)
111    v1lay(1:knon) = v_old(1:knon,1)*alf1(1:knon) + v_old(1:knon,2)*alf2(1:knon)   
112
113! - Calculate K
114    Kcoefm(:,:) = 0.0
115    DO i = 1, knon
116       Kcoefm(i,1) = coef_in(i,1) * (1.0+SQRT(u1lay(i)**2+v1lay(i)**2))* &
117            pplay(i,1)/(RD*temp(i,1))
118       Kcoefm(i,1) = Kcoefm(i,1) * dtime*RG
119    END DO
120
121    DO k = 2, klev
122       DO i=1,knon
123          Kcoefm(i,k) = coef_in(i,k)*RG/(pplay(i,k-1)-pplay(i,k)) &
124               *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
125          Kcoefm(i,k) = Kcoefm(i,k) * dtime*RG
126       END DO
127    END DO
128
129! - Calculate the coefficients C and D, component "u"
130    CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), &
131         u_old(:,:), alf1(:), alf2(:),  &
132         Ccoef_U(:,:), Dcoef_U(:,:))
133
134! - Calculate the coefficients C and D, component "v"
135    CALL calc_coef(knon, Kcoefm(:,:), delp(:,:), &
136         v_old(:,:), alf1(:), alf2(:),  &
137         Ccoef_V(:,:), Dcoef_V(:,:))
138
139  END SUBROUTINE climb_wind_down
140!
141!****************************************************************************************
142!
143  SUBROUTINE calc_coef(knon, Kcoef, dels, X, alfa1, alfa2, Ccoef, Dcoef)
144!
145! Find the coefficients C and D in fonction of alfa, K and dels
146!
147! Input arguments
148!****************************************************************************************
149    INTEGER, INTENT(IN)                      :: knon
150    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, dels
151    REAL, DIMENSION(klon,klev), INTENT(IN)   :: X
152    REAL, DIMENSION(klon), INTENT(IN)        :: alfa1, alfa2
153
154! Output arguments
155!****************************************************************************************
156    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef, Dcoef
157 
158! local variables
159!****************************************************************************************
160    INTEGER                                  :: k, i
161    REAL                                     :: buf
162
163!****************************************************************************************
164!
165! Niveau au sommet, k=klev
166!
167    Ccoef(:,:) = 0.0
168    Dcoef(:,:) = 0.0
169
170    DO i = 1, knon
171       buf = dels(i,klev) + Kcoef(i,klev)
172
173       Ccoef(i,klev) = X(i,klev)*dels(i,klev)/buf
174       Dcoef(i,klev) = Kcoef(i,klev)/buf
175    END DO
176   
177!
178! Niveau  (klev-1) <= k <= 2
179!
180    DO k=(klev-1),2,-1
181       DO i = 1, knon
182          buf = dels(i,k) + Kcoef(i,k) + &
183               Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
184         
185          Ccoef(i,k) = X(i,k)*dels(i,k) + &
186               Kcoef(i,k+1)*Ccoef(i,k+1)
187          Ccoef(i,k) = Ccoef(i,k)/buf         
188          Dcoef(i,k) = Kcoef(i,k)/buf
189       END DO
190    END DO
191
192!
193! Niveau k=1
194!
195    DO i = 1, knon
196       buf = dels(i,1) + &
197            (alfa1(i) + alfa2(i)*Dcoef(i,2)) * Kcoef(i,1) + &
198            (1.-Dcoef(i,2))*Kcoef(i,2)
199       
200       Ccoef(i,1) = X(i,1)*dels(i,1) + &
201            (Kcoef(i,2)-Kcoef(i,1)*alfa2(i)) * Ccoef(i,2)
202       Ccoef(i,1) = Ccoef(i,1)/buf
203       Dcoef(i,1) = Kcoef(i,1)/buf
204    END DO
205
206  END SUBROUTINE calc_coef
207!
208!****************************************************************************************
209!
210  SUBROUTINE calcul_wind_flux(knon, dtime, &
211       flux_u, flux_v)
212
213    INCLUDE "YOMCST.h"
214
215! Input arguments
216!****************************************************************************************
217    INTEGER, INTENT(IN)                  :: knon
218    REAL, INTENT(IN)                     :: dtime
219
220! Output arguments
221!****************************************************************************************
222    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u
223    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_v
224
225! Local variables
226!****************************************************************************************
227    INTEGER                              :: i
228    REAL, DIMENSION(klon)                :: u0, v0  ! u and v at niveau 0
229    REAL, DIMENSION(klon)                :: u1, v1  ! u and v at niveau 1
230    REAL, DIMENSION(klon)                :: u2, v2  ! u and v at niveau 2
231   
232
233!****************************************************************************************
234! Les vents de surface sont supposes nuls
235!
236!****************************************************************************************
237    u0(:) = 0.0
238    v0(:) = 0.0
239
240!****************************************************************************************
241! On calcule les vents du couhes 1 et 2 recurviement
242!
243!****************************************************************************************
244    DO i = 1, knon
245       u1(i) = Ccoef_U(i,1) + Dcoef_U(i,1)*u0(i)
246       v1(i) = Ccoef_V(i,1) + Dcoef_V(i,1)*v0(i)
247       u2(i) = Ccoef_U(i,2) + Dcoef_U(i,2)*u1(i)
248       v2(i) = Ccoef_V(i,2) + Dcoef_V(i,2)*v1(i)
249    END DO
250
251!****************************************************************************************
252! On calcule le flux
253!
254!****************************************************************************************
255    flux_u(:) = 0.0
256    flux_v(:) = 0.0
257
258    DO i=1,knon
259       flux_u(i) = Kcoefm(i,1)/RG/dtime * (u1(i)*alf1(i) + u2(i)*alf2(i) - u0(i))
260       flux_v(i) = Kcoefm(i,1)/RG/dtime * (v1(i)*alf1(i) + v2(i)*alf2(i) - v0(i))
261    END DO
262
263  END SUBROUTINE calcul_wind_flux
264!
265!****************************************************************************************
266!
267  SUBROUTINE climb_wind_up(knon, dtime, u_old, v_old, &
268       flx_u_new, flx_v_new, d_u_new, d_v_new)
269!
270! Diffuse the wind components from the surface and up to the top layer. Coefficents
271! C and D are known from before. The values for U and V at surface are supposed to be
272! zero (this could be modified).
273!
274! u(k) = Cu(k) + Du(k)*u(k-1)
275! v(k) = Cv(k) + Dv(k)*v(k-1)
276! [1 <= k <= klev]
277!
278!****************************************************************************************
279    INCLUDE "YOMCST.h"
280
281! Input arguments
282!****************************************************************************************
283    INTEGER, INTENT(IN)                     :: knon
284    REAL, INTENT(IN)                        :: dtime
285    REAL, DIMENSION(klon,klev), INTENT(IN)  :: u_old
286    REAL, DIMENSION(klon,klev), INTENT(IN)  :: v_old
287
288! Output arguments
289!****************************************************************************************
290    REAL, DIMENSION(klon,klev), INTENT(OUT) :: flx_u_new, flx_v_new
291    REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_u_new, d_v_new
292
293! Local variables
294!****************************************************************************************
295    REAL, DIMENSION(klon,klev)              :: u_new, v_new
296    REAL, DIMENSION(klon)                   :: u0, v0
297    INTEGER                                 :: k, i
298   
299!
300!****************************************************************************************
301
302! Niveau 0
303    u0(1:knon) = 0.0
304    v0(1:knon) = 0.0
305
306! Niveau 1
307    DO i = 1, knon
308       u_new(i,1) = Ccoef_U(i,1) + Dcoef_U(i,1) * u0(i)
309       v_new(i,1) = Ccoef_V(i,1) + Dcoef_V(i,1) * v0(i)
310    END DO
311
312! Niveau 2 jusqu'au sommet klev
313    DO k = 2, klev
314       DO i=1, knon
315          u_new(i,k) = Ccoef_U(i,k) + Dcoef_U(i,k) * u_new(i,k-1)
316          v_new(i,k) = Ccoef_V(i,k) + Dcoef_V(i,k) * v_new(i,k-1)
317       END DO
318    END DO
319
320!****************************************************************************************
321! Calcul flux
322!
323!== flux_u/v est le flux de moment angulaire (positif vers bas)
324!== dont l'unite est: (kg m/s)/(m**2 s)
325!
326!****************************************************************************************
327!
328    flx_u_new(:,:) = 0.0
329    flx_v_new(:,:) = 0.0
330
331! Niveau 1
332    DO i = 1, knon
333       flx_u_new(i,1) = Kcoefm(i,1)/RG/dtime * &
334            (u_new(i,1)*alf1(i)+u_new(i,2)*alf2(i) - u0(i))
335       flx_v_new(i,1) = Kcoefm(i,1)/RG/dtime * &
336            (v_new(i,1)*alf1(i)+v_new(i,2)*alf2(i) - v0(i))
337    END DO
338
339! Niveau 2->klev
340    DO k = 2, klev
341       DO i = 1, knon
342          flx_u_new(i,k) = Kcoefm(i,k)/RG/dtime * &
343               (u_new(i,k)-u_new(i,k-1))
344         
345          flx_v_new(i,k) = Kcoefm(i,k)/RG/dtime * &
346               (v_new(i,k)-v_new(i,k-1))
347       END DO
348    END DO
349
350!****************************************************************************************
351! Calcul tendances
352!
353!****************************************************************************************
354    d_u_new(:,:) = 0.0
355    d_v_new(:,:) = 0.0
356    DO k = 1, klev
357       DO i = 1, knon
358          d_u_new(i,k) = u_new(i,k) - u_old(i,k)
359          d_v_new(i,k) = v_new(i,k) - v_old(i,k)
360       END DO
361    END DO
362
363  END SUBROUTINE climb_wind_up
364!
365!****************************************************************************************
366!
367END MODULE climb_wind_mod
Note: See TracBrowser for help on using the repository browser.