Context Navigation

climb_wind_mod.F90 @ 814

Last change on this file since 814 was 814, checked in by Laurent Fairhead, 17 years ago
Rajout de la physique utilisant les thermiques FH 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	!
4	MODULE 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
29	CONTAINS
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	!
367	END MODULE climb_wind_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format