Context Navigation

advect_new_p.F @ 802

Last change on this file since 802 was 774, checked in by Laurent Fairhead, 17 years ago
Suite du merge entre la version et la HEAD: quelques modifications de Yann sur le LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 7.0 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE advect_new_p(ucov,vcov,teta,w,massebx,masseby,
5	& du,dv,dteta)
6	USE parallel
7	USE write_field_p
8	IMPLICIT NONE
9	c=======================================================================
10	c
11	c Auteurs: P. Le Van , Fr. Hourdin .
12	c -------
13	c
14	c Objet:
15	c ------
16	c
17	c *************************************************************
18	c .... calcul des termes d'advection vertic.pour u,v,teta,q ...
19	c *************************************************************
20	c ces termes sont ajoutes a du,dv,dteta et dq .
21	c Modif F.Forget 03/94 : on retire q de advect
22	c
23	c=======================================================================
24	c-----------------------------------------------------------------------
25	c Declarations:
26	c -------------
27
28	#include "dimensions.h"
29	#include "paramet.h"
30	#include "comconst.h"
31	#include "comvert.h"
32	#include "comgeom.h"
33	#include "logic.h"
34	#include "ener.h"
35
36	c Arguments:
37	c ----------
38
39	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm),teta(ip1jmp1,llm)
40	REAL massebx(ip1jmp1,llm),masseby(ip1jm,llm),w(ip1jmp1,llm)
41	REAL dv(ip1jm,llm),du(ip1jmp1,llm),dteta(ip1jmp1,llm)
42	REAL,SAVE :: dv1(ip1jm,llm),du1(ip1jmp1,llm),dteta1(ip1jmp1,llm)
43	REAL,SAVE :: dv2(ip1jm,llm),du2(ip1jmp1,llm),dteta2(ip1jmp1,llm)
44	c Local:
45	c ------
46
47	REAL,SAVE :: uav(ip1jmp1,llm),vav(ip1jm,llm)
48	REAL wsur2(ip1jmp1)
49	REAL unsaire2(ip1jmp1), ge(ip1jmp1)
50	REAL deuxjour, ww, gt, uu, vv
51
52	INTEGER ij,l,ijb,ije
53
54	EXTERNAL SSUM
55	REAL SSUM
56
57	c-----------------------------------------------------------------------
58	c 2. Calculs preliminaires:
59	c -------------------------
60
61	IF (conser) THEN
62	deuxjour = 2. * daysec
63
64	DO 1 ij = 1, ip1jmp1
65	unsaire2(ij) = unsaire(ij) * unsaire(ij)
66	1 CONTINUE
67	END IF
68
69
70	c------------------ -yy ----------------------------------------------
71	c . Calcul de u
72
73	c$OMP MASTER
74	ijb=ij_begin
75	ije=ij_end
76	if (pole_nord) ijb=ijb+iip1
77	if (pole_sud) ije=ije-iip1
78
79	DO ij=ijb,ije
80	du2(ij,1)=0.
81	ENDDO
82
83	ijb=ij_begin
84	ije=ij_end
85	if (pole_sud) ije=ij_end-iip1
86
87	DO ij=ijb,ije
88	dv2(ij,1)=0.
89	ENDDO
90
91	ijb=ij_begin
92	ije=ij_end
93
94	DO ij=ijb,ije
95	dteta2(ij,1)=0.
96	ENDDO
97	c$OMP END MASTER
98
99
100	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
101	DO l=1,llm
102
103	ijb=ij_begin
104	ije=ij_end
105	if (pole_nord) ijb=ijb+iip1
106	if (pole_sud) ije=ije-iip1
107
108	c DO ij = iip2, ip1jmp1
109	c uav(ij,l) = 0.25 * ( ucov(ij,l) + ucov(ij-iip1,l) )
110	c ENDDO
111
112	c DO ij = iip2, ip1jm
113	c uav(ij,l) = uav(ij,l) + uav(ij+iip1,l)
114	c ENDDO
115
116	DO ij = ijb, ije
117
118	uav(ij,l)=0.25*(ucov(ij,l)+ucov(ij-iip1,l))
119	. +0.25*(ucov(ij+iip1,l)+ucov(ij,l))
120	ENDDO
121
122	if (pole_nord) then
123	DO ij = 1, iip1
124	uav(ij ,l) = 0.
125	ENDDO
126	endif
127
128	if (pole_sud) then
129	DO ij = 1, iip1
130	uav(ip1jm+ij,l) = 0.
131	ENDDO
132	endif
133
134	ENDDO
135	c$OMP END DO
136	c call write_field3d_p('uav',reshape(uav,(/iip1,jjp1,llm/)))
137
138	c------------------ -xx ----------------------------------------------
139	c . Calcul de v
140
141	ijb=ij_begin
142	ije=ij_end
143	if (pole_sud) ije=ij_end-iip1
144
145	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
146	DO l=1,llm
147
148	DO ij = ijb+1, ije
149	vav(ij,l) = 0.25 * ( vcov(ij,l) + vcov(ij-1,l) )
150	ENDDO
151
152	DO ij = ijb,ije,iip1
153	vav(ij,l) = vav(ij+iim,l)
154	ENDDO
155
156
157	DO ij = ijb, ije-1
158	vav(ij,l) = vav(ij,l) + vav(ij+1,l)
159	ENDDO
160
161	DO ij = ijb, ije, iip1
162	vav(ij+iim,l) = vav(ij,l)
163	ENDDO
164
165	ENDDO
166	c$OMP END DO
167	c call write_field3d_p('vav',reshape(vav,(/iip1,jjm,llm/)))
168
169	c-----------------------------------------------------------------------
170	c$OMP BARRIER
171
172	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
173	DO 20 l = 1, llmm1
174
175
176	c ...... calcul de - w/2. au niveau l+1 .......
177	ijb=ij_begin
178	ije=ij_end+iip1
179	if (pole_sud) ije=ij_end
180
181	DO 5 ij = ijb, ije
182	wsur2( ij ) = - 0.5 * w( ij,l+1 )
183	5 CONTINUE
184
185
186	c ..................... calcul pour du ..................
187
188	ijb=ij_begin
189	ije=ij_end
190	if (pole_nord) ijb=ijb+iip1
191	if (pole_sud) ije=ije-iip1
192
193	DO 6 ij = ijb ,ije-1
194	ww = wsur2 ( ij ) + wsur2( ij+1 )
195	uu = 0.5 * ( ucov(ij,l) + ucov(ij,l+1) )
196	du1(ij,l) = ww * ( uu - uav(ij, l ) )/massebx(ij, l )
197	du2(ij,l+1)= ww * ( uu - uav(ij,l+1) )/massebx(ij,l+1)
198	6 CONTINUE
199
200	c ................. calcul pour dv .....................
201	ijb=ij_begin
202	ije=ij_end
203	if (pole_sud) ije=ij_end-iip1
204
205	DO 8 ij = ijb, ije
206	ww = wsur2( ij+iip1 ) + wsur2( ij )
207	vv = 0.5 * ( vcov(ij,l) + vcov(ij,l+1) )
208	dv1(ij,l) = ww * (vv - vav(ij, l ) )/masseby(ij, l )
209	dv2(ij,l+1)= ww * (vv - vav(ij,l+1) )/masseby(ij,l+1)
210	8 CONTINUE
211
212	c
213
214	c ............................................................
215	c ............... calcul pour dh ...................
216	c ............................................................
217
218	c ---z
219	c calcul de - d( teta * w ) qu'on ajoute a dh
220	c ...............
221	ijb=ij_begin
222	ije=ij_end
223
224	DO 15 ij = ijb, ije
225	ww = wsur2(ij) * (teta(ij,l) + teta(ij,l+1) )
226	dteta1(ij, l ) = ww
227	dteta2(ij,l+1) = ww
228	15 CONTINUE
229
230	c ym ---> conser a voir plus tard
231
232	c IF( conser) THEN
233	c
234	c DO 17 ij = 1,ip1jmp1
235	c ge(ij) = wsur2(ij) * wsur2(ij) * unsaire2(ij)
236	c 17 CONTINUE
237	c gt = SSUM( ip1jmp1,ge,1 )
238	c gtot(l) = deuxjour * SQRT( gt/ip1jmp1 )
239	c END IF
240
241	20 CONTINUE
242	c$OMP END DO
243
244	ijb=ij_begin
245	ije=ij_end
246	if (pole_nord) ijb=ijb+iip1
247	if (pole_sud) ije=ije-iip1
248
249	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
250	DO l=1,llm
251	DO ij=ijb,ije-1
252	du(ij,l)=du(ij,l)+du2(ij,l)-du1(ij,l)
253	ENDDO
254
255	DO ij = ijb+iip1-1, ije, iip1
256	du( ij, l ) = du( ij -iim, l )
257	ENDDO
258	ENDDO
259	c$OMP END DO NOWAIT
260
261	ijb=ij_begin
262	ije=ij_end
263	if (pole_sud) ije=ij_end-iip1
264
265	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
266	DO l=1,llm
267	DO ij=ijb,ije
268	dv(ij,l)=dv(ij,l)+dv2(ij,l)-dv1(ij,l)
269	ENDDO
270	ENDDO
271	c$OMP END DO NOWAIT
272	ijb=ij_begin
273	ije=ij_end
274
275	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
276	DO l=1,llm
277	DO ij=ijb,ije
278	dteta(ij,l)=dteta(ij,l)+dteta2(ij,l)-dteta1(ij,l)
279	ENDDO
280	ENDDO
281	c$OMP END DO NOWAIT
282
283	RETURN
284	END

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

Download in other formats:

Original Format