Context Navigation

vlsplt_loc.F @ 4106

Last change on this file since 4106 was 1673, checked in by Laurent Fairhead, 12 years ago

Fin du phasage de la dynamique parallele localisee (petite memoire) avec le tronc LMDZ5 r1671
Il reste quelques routines a verifier (en particulier ce qui touche a l'etude des cas academiques)
et la validation a effectuer

End of the phasing of the localised (low memory) parallel dynamics package with the
LMDZ5 trunk (r1671)
Some routines still need some checking (in particular the academic cases) and some
validation is still required

File size: 24.1 KB

Line
1	!
2	! $Id$
3	!
4	SUBROUTINE vlx_loc(q,pente_max,masse,u_m,ijb_x,ije_x)
5
6	c Auteurs: P.Le Van, F.Hourdin, F.Forget
7	c
8	c ********************************************************************
9	c Shema d'advection " pseudo amont " .
10	c ********************************************************************
11	c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
12	c
13	c
14	c --------------------------------------------------------------------
15	USE Parallel
16	IMPLICIT NONE
17	c
18	#include "dimensions.h"
19	#include "paramet.h"
20	#include "logic.h"
21	#include "comvert.h"
22	#include "comconst.h"
23	c
24	c
25	c Arguments:
26	c ----------
27	REAL masse(ijb_u:ije_u,llm),pente_max
28	REAL u_m( ijb_u:ije_u,llm ),pbarv( iip1,jjb_v:jje_v,llm)
29	REAL q(ijb_u:ije_u,llm)
30	REAL w(ijb_u:ije_u,llm)
31	c
32	c Local
33	c ---------
34	c
35	INTEGER ij,l,j,i,iju,ijq,indu(ijnb_u),niju
36	INTEGER n0,iadvplus(ijb_u:ije_u,llm),nl(llm)
37	c
38	REAL new_m,zu_m,zdum(ijb_u:ije_u,llm)
39	REAL sigu(ijb_u:ije_u),dxq(ijb_u:ije_u,llm),dxqu(ijb_u:ije_u)
40	REAL zz(ijb_u:ije_u)
41	REAL adxqu(ijb_u:ije_u),dxqmax(ijb_u:ije_u,llm)
42	REAL u_mq(ijb_u:ije_u,llm)
43
44	Logical extremum
45
46	REAL SSUM
47	EXTERNAL SSUM
48
49	REAL z1,z2,z3
50
51	INTEGER ijb,ije,ijb_x,ije_x
52
53	c calcul de la pente a droite et a gauche de la maille
54
55	ijb=ijb_x
56	ije=ije_x
57
58	if (pole_nord.and.ijb==1) ijb=ijb+iip1
59	if (pole_sud.and.ije==ip1jmp1) ije=ije-iip1
60
61	IF (pente_max.gt.-1.e-5) THEN
62	c IF (pente_max.gt.10) THEN
63
64	c calcul des pentes avec limitation, Van Leer scheme I:
65	c -----------------------------------------------------
66
67	c calcul de la pente aux points u
68	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
69	DO l = 1, llm
70
71	DO ij=ijb,ije-1
72	dxqu(ij)=q(ij+1,l)-q(ij,l)
73	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
74	c sigu(ij)=u_m(ij,l)/masse(ij,l)
75	ENDDO
76	DO ij=ijb+iip1-1,ije,iip1
77	dxqu(ij)=dxqu(ij-iim)
78	c sigu(ij)=sigu(ij-iim)
79	ENDDO
80
81	DO ij=ijb,ije
82	adxqu(ij)=abs(dxqu(ij))
83	ENDDO
84
85	c calcul de la pente maximum dans la maille en valeur absolue
86
87	DO ij=ijb+1,ije
88	dxqmax(ij,l)=pente_max*
89	, min(adxqu(ij-1),adxqu(ij))
90	c limitation subtile
91	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
92
93
94	ENDDO
95
96	DO ij=ijb+iip1-1,ije,iip1
97	dxqmax(ij-iim,l)=dxqmax(ij,l)
98	ENDDO
99
100	DO ij=ijb+1,ije
101	#ifdef CRAY
102	dxq(ij,l)=
103	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
104	#else
105	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
106	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
107	ELSE
108	c extremum local
109	dxq(ij,l)=0.
110	ENDIF
111	#endif
112	dxq(ij,l)=0.5*dxq(ij,l)
113	dxq(ij,l)=
114	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
115	ENDDO
116
117	ENDDO ! l=1,llm
118	c$OMP END DO NOWAIT
119	c print*,'Ok calcul des pentes'
120
121	ELSE ! (pente_max.lt.-1.e-5)
122
123	c Pentes produits:
124	c ----------------
125	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
126	DO l = 1, llm
127	DO ij=ijb,ije-1
128	dxqu(ij)=q(ij+1,l)-q(ij,l)
129	ENDDO
130	DO ij=ijb+iip1-1,ije,iip1
131	dxqu(ij)=dxqu(ij-iim)
132	ENDDO
133
134	DO ij=ijb+1,ije
135	zz(ij)=dxqu(ij-1)*dxqu(ij)
136	zz(ij)=zz(ij)+zz(ij)
137	IF(zz(ij).gt.0) THEN
138	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
139	ELSE
140	c extremum local
141	dxq(ij,l)=0.
142	ENDIF
143	ENDDO
144
145	ENDDO
146	c$OMP END DO NOWAIT
147	ENDIF ! (pente_max.lt.-1.e-5)
148
149	c bouclage de la pente en iip1:
150	c -----------------------------
151	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
152	DO l=1,llm
153	DO ij=ijb+iip1-1,ije,iip1
154	dxq(ij-iim,l)=dxq(ij,l)
155	ENDDO
156	DO ij=ijb,ije
157	iadvplus(ij,l)=0
158	ENDDO
159
160	ENDDO
161	c$OMP END DO NOWAIT
162	c print*,'Bouclage en iip1'
163
164	c calcul des flux a gauche et a droite
165
166	#ifdef CRAY
167	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
168	DO l=1,llm
169	DO ij=ijb,ije-1
170	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
171	, 1.+u_m(ij,l)/masse(ij+1,l),
172	, u_m(ij,l))
173	zdum(ij,l)=0.5*zdum(ij,l)
174	u_mq(ij,l)=cvmgp(
175	, q(ij,l)+zdum(ij,l)*dxq(ij,l),
176	, q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
177	, u_m(ij,l))
178	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
179	ENDDO
180	ENDDO
181	c$OMP END DO NOWAIT
182	#else
183	c on cumule le flux correspondant a toutes les mailles dont la masse
184	c au travers de la paroi pENDant le pas de temps.
185	c print*,'Cumule ....'
186	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
187	DO l=1,llm
188	DO ij=ijb,ije-1
189	c print*,'masse(',ij,')=',masse(ij,l)
190	IF (u_m(ij,l).gt.0.) THEN
191	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
192	u_mq(ij,l)=u_m(ij,l)(q(ij,l)+0.5zdum(ij,l)*dxq(ij,l))
193	ELSE
194	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
195	u_mq(ij,l)=u_m(ij,l)(q(ij+1,l)-0.5zdum(ij,l)*dxq(ij+1,l))
196	ENDIF
197	ENDDO
198	ENDDO
199	c$OMP END DO NOWAIT
200	#endif
201	c stop
202
203	c go to 9999
204	c detection des points ou on advecte plus que la masse de la
205	c maille
206	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
207	DO l=1,llm
208	DO ij=ijb,ije-1
209	IF(zdum(ij,l).lt.0) THEN
210	iadvplus(ij,l)=1
211	u_mq(ij,l)=0.
212	ENDIF
213	ENDDO
214	ENDDO
215	c$OMP END DO NOWAIT
216	c print*,'Ok test 1'
217	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
218	DO l=1,llm
219	DO ij=ijb+iip1-1,ije,iip1
220	iadvplus(ij,l)=iadvplus(ij-iim,l)
221	ENDDO
222	ENDDO
223	c$OMP END DO NOWAIT
224	c print*,'Ok test 2'
225
226
227	c traitement special pour le cas ou on advecte en longitude plus que le
228	c contenu de la maille.
229	c cette partie est mal vectorisee.
230
231	c calcul du nombre de maille sur lequel on advecte plus que la maille.
232
233	n0=0
234	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
235	DO l=1,llm
236	nl(l)=0
237	DO ij=ijb,ije
238	nl(l)=nl(l)+iadvplus(ij,l)
239	ENDDO
240	n0=n0+nl(l)
241	ENDDO
242	c$OMP END DO NOWAIT
243	cym IF(n0.gt.1) THEN
244	cym IF(n0.gt.0) THEN
245
246	c PRINT*,'Nombre de points pour lesquels on advect plus que le'
247	c & ,'contenu de la maille : ',n0
248	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
249	DO l=1,llm
250	IF(nl(l).gt.0) THEN
251	iju=0
252	c indicage des mailles concernees par le traitement special
253	DO ij=ijb,ije
254	IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
255	iju=iju+1
256	indu(iju)=ij
257	ENDIF
258	ENDDO
259	niju=iju
260	c PRINT*,'niju,nl',niju,nl(l)
261
262	c traitement des mailles
263	DO iju=1,niju
264	ij=indu(iju)
265	j=(ij-1)/iip1+1
266	zu_m=u_m(ij,l)
267	u_mq(ij,l)=0.
268	IF(zu_m.gt.0.) THEN
269	ijq=ij
270	i=ijq-(j-1)*iip1
271	c accumulation pour les mailles completements advectees
272	do while(zu_m.gt.masse(ijq,l))
273	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
274	zu_m=zu_m-masse(ijq,l)
275	i=mod(i-2+iim,iim)+1
276	ijq=(j-1)*iip1+i
277	ENDDO
278	c ajout de la maille non completement advectee
279	u_mq(ij,l)=u_mq(ij,l)+zu_m*
280	& (q(ijq,l)+0.5(1.-zu_m/masse(ijq,l))dxq(ijq,l))
281	ELSE
282	ijq=ij+1
283	i=ijq-(j-1)*iip1
284	c accumulation pour les mailles completements advectees
285	do while(-zu_m.gt.masse(ijq,l))
286	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
287	zu_m=zu_m+masse(ijq,l)
288	i=mod(i,iim)+1
289	ijq=(j-1)*iip1+i
290	ENDDO
291	c ajout de la maille non completement advectee
292	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
293	& 0.5(1.+zu_m/masse(ijq,l))dxq(ijq,l))
294	ENDIF
295	ENDDO
296	ENDIF
297	ENDDO
298	c$OMP END DO NOWAIT
299	cym ENDIF ! n0.gt.0
300	9999 continue
301
302
303	c bouclage en latitude
304	c print*,'Avant bouclage en latitude'
305	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
306	DO l=1,llm
307	DO ij=ijb+iip1-1,ije,iip1
308	u_mq(ij,l)=u_mq(ij-iim,l)
309	ENDDO
310	ENDDO
311	c$OMP END DO NOWAIT
312
313	c calcul des tENDances
314	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
315	DO l=1,llm
316	DO ij=ijb+1,ije
317	new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
318	q(ij,l)=(q(ij,l)*masse(ij,l)+
319	& u_mq(ij-1,l)-u_mq(ij,l))
320	& /new_m
321	masse(ij,l)=new_m
322	ENDDO
323	c ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
324	DO ij=ijb+iip1-1,ije,iip1
325	q(ij-iim,l)=q(ij,l)
326	masse(ij-iim,l)=masse(ij,l)
327	ENDDO
328	ENDDO
329	c$OMP END DO NOWAIT
330	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
331	c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
332
333
334	RETURN
335	END
336
337
338	SUBROUTINE vly_loc(q,pente_max,masse,masse_adv_v)
339	c
340	c Auteurs: P.Le Van, F.Hourdin, F.Forget
341	c
342	c ********************************************************************
343	c Shema d'advection " pseudo amont " .
344	c ********************************************************************
345	c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
346	c dq sont des arguments de sortie pour le s-pg ....
347	c
348	c
349	c --------------------------------------------------------------------
350	USE parallel
351	IMPLICIT NONE
352	c
353	#include "dimensions.h"
354	#include "paramet.h"
355	#include "logic.h"
356	#include "comvert.h"
357	#include "comconst.h"
358	#include "comgeom.h"
359	c
360	c
361	c Arguments:
362	c ----------
363	REAL masse(ijb_u:ije_u,llm),pente_max
364	REAL masse_adv_v( ijb_v:ije_v,llm)
365	REAL q(ijb_u:ije_u,llm), dq( ijb_u:ije_u,llm)
366	c
367	c Local
368	c ---------
369	c
370	INTEGER i,ij,l
371	c
372	REAL airej2,airejjm,airescb(iim),airesch(iim)
373	REAL dyq(ijb_u:ije_u,llm),dyqv(ijb_v:ije_v),zdvm(ijb_u:ije_u,llm)
374	REAL adyqv(ijb_v:ije_v),dyqmax(ijb_u:ije_u)
375	REAL qbyv(ijb_v:ije_v,llm)
376
377	REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
378	c REAL newq,oldmasse
379	Logical extremum,first,testcpu
380	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
381	SAVE temps0,temps1,temps2,temps3,temps4,temps5
382	c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
383	SAVE first,testcpu
384	c$OMP THREADPRIVATE(first,testcpu)
385
386	REAL convpn,convps,convmpn,convmps
387	real massepn,masseps,qpn,qps
388	REAL sinlon(iip1),sinlondlon(iip1)
389	REAL coslon(iip1),coslondlon(iip1)
390	SAVE sinlon,coslon,sinlondlon,coslondlon
391	c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
392	SAVE airej2,airejjm
393	c$OMP THREADPRIVATE(airej2,airejjm)
394	c
395	c
396	REAL SSUM
397	EXTERNAL SSUM
398
399	DATA first,testcpu/.true.,.false./
400	DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
401	INTEGER ijb,ije
402
403	IF(first) THEN
404	c PRINT*,'Shema Amont nouveau appele dans Vanleer '
405	first=.false.
406	do i=2,iip1
407	coslon(i)=cos(rlonv(i))
408	sinlon(i)=sin(rlonv(i))
409	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
410	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
411	ENDDO
412	coslon(1)=coslon(iip1)
413	coslondlon(1)=coslondlon(iip1)
414	sinlon(1)=sinlon(iip1)
415	sinlondlon(1)=sinlondlon(iip1)
416	airej2 = SSUM( iim, aire(iip2), 1 )
417	airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
418	ENDIF
419
420	c
421	c PRINT*,'CALCUL EN LATITUDE'
422
423	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
424	DO l = 1, llm
425	c
426	c --------------------------------
427	c CALCUL EN LATITUDE
428	c --------------------------------
429
430	c On commence par calculer la valeur du traceur moyenne sur le premier cercle
431	c de latitude autour du pole (qpns pour le pole nord et qpsn pour
432	c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
433
434	if (pole_nord) then
435	DO i = 1, iim
436	airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
437	ENDDO
438	qpns = SSUM( iim, airescb ,1 ) / airej2
439	endif
440
441	if (pole_sud) then
442	DO i = 1, iim
443	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
444	ENDDO
445	qpsn = SSUM( iim, airesch ,1 ) / airejjm
446	endif
447
448
449
450	c calcul des pentes aux points v
451
452	ijb=ij_begin-2*iip1
453	ije=ij_end+iip1
454	if (pole_nord) ijb=ij_begin
455	if (pole_sud) ije=ij_end-iip1
456
457	DO ij=ijb,ije
458	dyqv(ij)=q(ij,l)-q(ij+iip1,l)
459	adyqv(ij)=abs(dyqv(ij))
460	ENDDO
461
462	c calcul des pentes aux points scalaires
463	ijb=ij_begin-iip1
464	ije=ij_end+iip1
465	if (pole_nord) ijb=ij_begin+iip1
466	if (pole_sud) ije=ij_end-iip1
467
468	DO ij=ijb,ije
469	dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
470	dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
471	dyqmax(ij)=pente_max*dyqmax(ij)
472	ENDDO
473
474	c calcul des pentes aux poles
475	IF (pole_nord) THEN
476	DO ij=1,iip1
477	dyq(ij,l)=qpns-q(ij+iip1,l)
478	ENDDO
479
480	dyn1=0.
481	dyn2=0.
482	DO ij=1,iim
483	dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
484	dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
485	ENDDO
486	DO ij=1,iip1
487	dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
488	ENDDO
489
490	DO ij=1,iip1
491	dyq(ij,l)=0.
492	ENDDO
493	c ym tout cela ne sert pas a grand chose
494	ENDIF
495
496	IF (pole_sud) THEN
497
498	DO ij=1,iip1
499	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
500	ENDDO
501
502	dys1=0.
503	dys2=0.
504
505	DO ij=1,iim
506	dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
507	dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
508	ENDDO
509
510	DO ij=1,iip1
511	dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
512	ENDDO
513
514	DO ij=1,iip1
515	dyq(ip1jm+ij,l)=0.
516	ENDDO
517	c ym tout cela ne sert pas a grand chose
518	ENDIF
519
520	c filtrage de la derivee
521
522	c calcul des pentes limites aux poles
523	c ym partie inutile
524	c goto 8888
525	c fn=1.
526	c fs=1.
527	c DO ij=1,iim
528	c IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
529	c fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
530	c ENDIF
531	c IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
532	c fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
533	c ENDIF
534	c ENDDO
535	c DO ij=1,iip1
536	c dyq(ij,l)=fn*dyq(ij,l)
537	c dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
538	c ENDDO
539	c 8888 continue
540
541
542	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
543	C En memoire de dIFferents tests sur la
544	C limitation des pentes aux poles.
545	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
546	C PRINT*,dyq(1)
547	C PRINT*,dyqv(iip1+1)
548	C appn=abs(dyq(1)/dyqv(iip1+1))
549	C PRINT*,dyq(ip1jm+1)
550	C PRINT*,dyqv(ip1jm-iip1+1)
551	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
552	C DO ij=2,iim
553	C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
554	C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
555	C ENDDO
556	C appn=min(pente_max/appn,1.)
557	C apps=min(pente_max/apps,1.)
558	C
559	C
560	C cas ou on a un extremum au pole
561	C
562	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
563	C & appn=0.
564	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
565	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
566	C & apps=0.
567	C
568	C limitation des pentes aux poles
569	C DO ij=1,iip1
570	C dyq(ij)=appn*dyq(ij)
571	C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
572	C ENDDO
573	C
574	C test
575	C DO ij=1,iip1
576	C dyq(iip1+ij)=0.
577	C dyq(ip1jm+ij-iip1)=0.
578	C ENDDO
579	C DO ij=1,ip1jmp1
580	C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
581	C ENDDO
582	C
583	C changement 10 07 96
584	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
585	C & THEN
586	C DO ij=1,iip1
587	C dyqmax(ij)=0.
588	C ENDDO
589	C ELSE
590	C DO ij=1,iip1
591	C dyqmax(ij)=pente_max*abs(dyqv(ij))
592	C ENDDO
593	C ENDIF
594	C
595	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
596	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
597	C &THEN
598	C DO ij=ip1jm+1,ip1jmp1
599	C dyqmax(ij)=0.
600	C ENDDO
601	C ELSE
602	C DO ij=ip1jm+1,ip1jmp1
603	C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
604	C ENDDO
605	C ENDIF
606	C fin changement 10 07 96
607	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
608
609	c calcul des pentes limitees
610	ijb=ij_begin-iip1
611	ije=ij_end+iip1
612	if (pole_nord) ijb=ij_begin+iip1
613	if (pole_sud) ije=ij_end-iip1
614
615	DO ij=ijb,ije
616	IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
617	dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
618	ELSE
619	dyq(ij,l)=0.
620	ENDIF
621	ENDDO
622
623	ENDDO
624	c$OMP END DO NOWAIT
625
626	ijb=ij_begin-iip1
627	ije=ij_end
628	if (pole_nord) ijb=ij_begin
629	if (pole_sud) ije=ij_end-iip1
630
631	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
632	DO l=1,llm
633	DO ij=ijb,ije
634	IF(masse_adv_v(ij,l).gt.0) THEN
635	qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
636	, 0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
637	ELSE
638	qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
639	, 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
640	ENDIF
641	qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
642	ENDDO
643	ENDDO
644	c$OMP END DO NOWAIT
645
646	ijb=ij_begin
647	ije=ij_end
648	if (pole_nord) ijb=ij_begin+iip1
649	if (pole_sud) ije=ij_end-iip1
650
651	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
652	DO l=1,llm
653	DO ij=ijb,ije
654	newmasse=masse(ij,l)
655	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
656
657	q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
658	& /newmasse
659	masse(ij,l)=newmasse
660	ENDDO
661	c.-. ancienne version
662	c convpn=SSUM(iim,qbyv(1,l),1)/apoln
663	c convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
664	if (pole_nord) then
665	convpn=SSUM(iim,qbyv(1,l),1)
666	convmpn=ssum(iim,masse_adv_v(1,l),1)
667	massepn=ssum(iim,masse(1,l),1)
668	qpn=0.
669	do ij=1,iim
670	qpn=qpn+masse(ij,l)*q(ij,l)
671	enddo
672	qpn=(qpn+convpn)/(massepn+convmpn)
673	do ij=1,iip1
674	q(ij,l)=qpn
675	enddo
676	endif
677
678	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
679	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
680
681	if (pole_sud) then
682
683	convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
684	convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
685	masseps=ssum(iim, masse(ip1jm+1,l),1)
686	qps=0.
687	do ij = ip1jm+1,ip1jmp1-1
688	qps=qps+masse(ij,l)*q(ij,l)
689	enddo
690	qps=(qps+convps)/(masseps+convmps)
691	do ij=ip1jm+1,ip1jmp1
692	q(ij,l)=qps
693	enddo
694	endif
695	c.-. fin ancienne version
696
697	c._. nouvelle version
698	c convpn=SSUM(iim,qbyv(1,l),1)
699	c convmpn=ssum(iim,masse_adv_v(1,l),1)
700	c oldmasse=ssum(iim,masse(1,l),1)
701	c newmasse=oldmasse+convmpn
702	c newq=(q(1,l)*oldmasse+convpn)/newmasse
703	c newmasse=newmasse/apoln
704	c DO ij = 1,iip1
705	c q(ij,l)=newq
706	c masse(ij,l)=newmasse*aire(ij)
707	c ENDDO
708	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
709	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
710	c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
711	c newmasse=oldmasse+convmps
712	c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
713	c newmasse=newmasse/apols
714	c DO ij = ip1jm+1,ip1jmp1
715	c q(ij,l)=newq
716	c masse(ij,l)=newmasse*aire(ij)
717	c ENDDO
718	c._. fin nouvelle version
719	ENDDO
720	c$OMP END DO NOWAIT
721
722	RETURN
723	END
724
725
726
727	SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x)
728	c
729	c Auteurs: P.Le Van, F.Hourdin, F.Forget
730	c
731	c ********************************************************************
732	c Shema d'advection " pseudo amont " .
733	c ********************************************************************
734	c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
735	c dq sont des arguments de sortie pour le s-pg ....
736	c
737	c
738	c --------------------------------------------------------------------
739	USE Parallel
740	USE vlz_mod
741	IMPLICIT NONE
742	c
743	#include "dimensions.h"
744	#include "paramet.h"
745	#include "logic.h"
746	#include "comvert.h"
747	#include "comconst.h"
748	c
749	c
750	c Arguments:
751	c ----------
752	REAL masse(ijb_u:ije_u,llm),pente_max
753	REAL q(ijb_u:ije_u,llm)
754	REAL w(ijb_u:ije_u,llm+1)
755	c
756	c Local
757	c ---------
758	c
759	INTEGER i,ij,l,j,ii
760	c
761	REAL newmasse
762
763	REAL dzqmax
764	REAL sigw
765
766	LOGICAL testcpu
767	SAVE testcpu
768	c$OMP THREADPRIVATE(testcpu)
769	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
770	SAVE temps0,temps1,temps2,temps3,temps4,temps5
771	c$OMP THREADPRIVATE(temps0,temps1,temps2,temps3,temps4,temps5)
772
773	REAL SSUM
774	EXTERNAL SSUM
775
776	DATA testcpu/.false./
777	DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
778	INTEGER ijb,ije,ijb_x,ije_x
779	LOGICAL,SAVE :: first=.TRUE.
780	!$OMP THREADPRIVATE(first)
781
782
783	IF (first) THEN
784	first=.FALSE.
785	ENDIF
786	c On oriente tout dans le sens de la pression c'est a dire dans le
787	c sens de W
788
789	#ifdef BIDON
790	IF(testcpu) THEN
791	temps0=second(0.)
792	ENDIF
793	#endif
794
795	ijb=ijb_x
796	ije=ije_x
797
798	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
799	DO l=2,llm
800	DO ij=ijb,ije
801	dzqw(ij,l)=q(ij,l-1)-q(ij,l)
802	adzqw(ij,l)=abs(dzqw(ij,l))
803	ENDDO
804	ENDDO
805	c$OMP END DO
806
807	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
808	DO l=2,llm-1
809	DO ij=ijb,ije
810	#ifdef CRAY
811	dzq(ij,l)=0.5*
812	, cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1))
813	#else
814	IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
815	dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
816	ELSE
817	dzq(ij,l)=0.
818	ENDIF
819	#endif
820	dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
821	dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
822	ENDDO
823	ENDDO
824	c$OMP END DO NOWAIT
825
826	c$OMP MASTER
827	DO ij=ijb,ije
828	dzq(ij,1)=0.
829	dzq(ij,llm)=0.
830	ENDDO
831	c$OMP END MASTER
832	c$OMP BARRIER
833	#ifdef BIDON
834	IF(testcpu) THEN
835	temps1=temps1+second(0.)-temps0
836	ENDIF
837	#endif
838	c ---------------------------------------------------------------
839	c .... calcul des termes d'advection verticale .......
840	c ---------------------------------------------------------------
841
842	c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq
843
844	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
845	DO l = 1,llm-1
846	do ij = ijb,ije
847	IF(w(ij,l+1).gt.0.) THEN
848	sigw=w(ij,l+1)/masse(ij,l+1)
849	wq(ij,l+1)=w(ij,l+1)(q(ij,l+1)+0.5(1.-sigw)*dzq(ij,l+1))
850	ELSE
851	sigw=w(ij,l+1)/masse(ij,l)
852	wq(ij,l+1)=w(ij,l+1)(q(ij,l)-0.5(1.+sigw)*dzq(ij,l))
853	ENDIF
854	ENDDO
855	ENDDO
856	c$OMP END DO NOWAIT
857
858	c$OMP MASTER
859	DO ij=ijb,ije
860	wq(ij,llm+1)=0.
861	wq(ij,1)=0.
862	ENDDO
863	c$OMP END MASTER
864	c$OMP BARRIER
865
866	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
867	DO l=1,llm
868	DO ij=ijb,ije
869	newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
870	q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
871	& /newmasse
872	masse(ij,l)=newmasse
873	ENDDO
874	ENDDO
875	c$OMP END DO NOWAIT
876
877
878	RETURN
879	END
880	c SUBROUTINE minmaxq(zq,qmin,qmax,comment)
881	c
882	c#include "dimensions.h"
883	c#include "paramet.h"
884
885	c CHARACTER() comment
886	c real qmin,qmax
887	c real zq(ip1jmp1,llm)
888
889	c INTEGER jadrs(ip1jmp1), jbad, k, i
890
891
892	c DO k = 1, llm
893	c jbad = 0
894	c DO i = 1, ip1jmp1
895	c IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
896	c jbad = jbad + 1
897	c jadrs(jbad) = i
898	c ENDIF
899	c ENDDO
900	c IF (jbad.GT.0) THEN
901	c PRINT*, comment
902	c DO i = 1, jbad
903	cc PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
904	c ENDDO
905	c ENDIF
906	c ENDDO
907
908	c return
909	c end
910
911
912
913

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

Download in other formats:

Original Format