Context Navigation

vlspltqs_loc.F @ 5018

Last change on this file since 5018 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: 19.2 KB

Line
1	SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x)
2	c
3	c Auteurs: P.Le Van, F.Hourdin, F.Forget
4	c
5	c ********************************************************************
6	c Shema d'advection " pseudo amont " .
7	c ********************************************************************
8	c
9	c --------------------------------------------------------------------
10	USE parallel
11	IMPLICIT NONE
12	c
13	#include "dimensions.h"
14	#include "paramet.h"
15	#include "logic.h"
16	#include "comvert.h"
17	#include "comconst.h"
18	c
19	c
20	c Arguments:
21	c ----------
22	REAL masse(ijb_u:ije_u,llm),pente_max
23	REAL u_m( ijb_u:ije_u,llm )
24	REAL q(ijb_u:ije_u,llm)
25	REAL qsat(ijb_u:ije_u,llm)
26	c
27	c Local
28	c ---------
29	c
30	INTEGER ij,l,j,i,iju,ijq,indu(ijnb_u),niju
31	INTEGER n0,iadvplus(ijb_u:ije_u,llm),nl(llm)
32	c
33	REAL new_m,zu_m,zdum(ijb_u:ije_u,llm)
34	REAL dxq(ijb_u:ije_u,llm),dxqu(ijb_u:ije_u)
35	REAL zz(ijb_u:ije_u)
36	REAL adxqu(ijb_u:ije_u),dxqmax(ijb_u:ije_u,llm)
37	REAL u_mq(ijb_u:ije_u,llm)
38
39	REAL SSUM
40
41
42	INTEGER ijb,ije,ijb_x,ije_x
43
44
45	c calcul de la pente a droite et a gauche de la maille
46
47	c ijb=ij_begin
48	c ije=ij_end
49
50	ijb=ijb_x
51	ije=ije_x
52
53	if (pole_nord.and.ijb==1) ijb=ijb+iip1
54	if (pole_sud.and.ije==ip1jmp1) ije=ije-iip1
55
56	IF (pente_max.gt.-1.e-5) THEN
57	c IF (pente_max.gt.10) THEN
58
59	c calcul des pentes avec limitation, Van Leer scheme I:
60	c -----------------------------------------------------
61
62	c calcul de la pente aux points u
63
64	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
65	DO l = 1, llm
66	DO ij=ijb,ije-1
67	dxqu(ij)=q(ij+1,l)-q(ij,l)
68	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
69	c sigu(ij)=u_m(ij,l)/masse(ij,l)
70	ENDDO
71	DO ij=ijb+iip1-1,ije,iip1
72	dxqu(ij)=dxqu(ij-iim)
73	c sigu(ij)=sigu(ij-iim)
74	ENDDO
75
76	DO ij=ijb,ije
77	adxqu(ij)=abs(dxqu(ij))
78	ENDDO
79
80	c calcul de la pente maximum dans la maille en valeur absolue
81
82	DO ij=ijb+1,ije
83	dxqmax(ij,l)=pente_max*
84	, min(adxqu(ij-1),adxqu(ij))
85	c limitation subtile
86	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
87
88
89	ENDDO
90
91	DO ij=ijb+iip1-1,ije,iip1
92	dxqmax(ij-iim,l)=dxqmax(ij,l)
93	ENDDO
94
95	DO ij=ijb+1,ije
96	#ifdef CRAY
97	dxq(ij,l)=
98	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
99	#else
100	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
101	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
102	ELSE
103	c extremum local
104	dxq(ij,l)=0.
105	ENDIF
106	#endif
107	dxq(ij,l)=0.5*dxq(ij,l)
108	dxq(ij,l)=
109	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
110	ENDDO
111
112	ENDDO ! l=1,llm
113	c$OMP END DO NOWAIT
114
115	ELSE ! (pente_max.lt.-1.e-5)
116
117	c Pentes produits:
118	c ----------------
119	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
120	DO l = 1, llm
121	DO ij=ijb,ije-1
122	dxqu(ij)=q(ij+1,l)-q(ij,l)
123	ENDDO
124	DO ij=ijb+iip1-1,ije,iip1
125	dxqu(ij)=dxqu(ij-iim)
126	ENDDO
127
128	DO ij=ijb+1,ije
129	zz(ij)=dxqu(ij-1)*dxqu(ij)
130	zz(ij)=zz(ij)+zz(ij)
131	IF(zz(ij).gt.0) THEN
132	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
133	ELSE
134	c extremum local
135	dxq(ij,l)=0.
136	ENDIF
137	ENDDO
138
139	ENDDO
140	c$OMP END DO NOWAIT
141	ENDIF ! (pente_max.lt.-1.e-5)
142
143	c bouclage de la pente en iip1:
144	c -----------------------------
145	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
146	DO l=1,llm
147	DO ij=ijb+iip1-1,ije,iip1
148	dxq(ij-iim,l)=dxq(ij,l)
149	ENDDO
150
151	DO ij=ijb,ije
152	iadvplus(ij,l)=0
153	ENDDO
154
155	ENDDO
156	c$OMP END DO NOWAIT
157
158	if (pole_nord) THEN
159	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
160	DO l=1,llm
161	iadvplus(1:iip1,l)=0
162	ENDDO
163	c$OMP END DO NOWAIT
164	endif
165
166	if (pole_sud) THEN
167	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
168	DO l=1,llm
169	iadvplus(ip1jm+1:ip1jmp1,l)=0
170	ENDDO
171	c$OMP END DO NOWAIT
172	endif
173
174	c calcul des flux a gauche et a droite
175
176	#ifdef CRAY
177	c--pas encore modification sur Qsat
178	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
179	DO l=1,llm
180	DO ij=ijb,ije-1
181	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
182	, 1.+u_m(ij,l)/masse(ij+1,l),
183	, u_m(ij,l))
184	zdum(ij,l)=0.5*zdum(ij,l)
185	u_mq(ij,l)=cvmgp(
186	, q(ij,l)+zdum(ij,l)*dxq(ij,l),
187	, q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
188	, u_m(ij,l))
189	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
190	ENDDO
191	ENDDO
192	c$OMP END DO NOWAIT
193
194	#else
195	c on cumule le flux correspondant a toutes les mailles dont la masse
196	c au travers de la paroi pENDant le pas de temps.
197	c le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
198	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
199	DO l=1,llm
200	DO ij=ijb,ije-1
201	IF (u_m(ij,l).gt.0.) THEN
202	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
203	u_mq(ij,l)=u_m(ij,l)*
204	$ min(q(ij,l)+0.5zdum(ij,l)dxq(ij,l),qsat(ij+1,l))
205	ELSE
206	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
207	u_mq(ij,l)=u_m(ij,l)*
208	$ min(q(ij+1,l)-0.5zdum(ij,l)dxq(ij+1,l),qsat(ij,l))
209	ENDIF
210	ENDDO
211	ENDDO
212	c$OMP END DO NOWAIT
213	#endif
214
215
216	c detection des points ou on advecte plus que la masse de la
217	c maille
218	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
219	DO l=1,llm
220	DO ij=ijb,ije-1
221	IF(zdum(ij,l).lt.0) THEN
222	iadvplus(ij,l)=1
223	u_mq(ij,l)=0.
224	ENDIF
225	ENDDO
226	ENDDO
227	c$OMP END DO NOWAIT
228
229	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
230	DO l=1,llm
231	DO ij=ijb+iip1-1,ije,iip1
232	iadvplus(ij,l)=iadvplus(ij-iim,l)
233	ENDDO
234	ENDDO
235	c$OMP END DO NOWAIT
236
237
238
239	c traitement special pour le cas ou on advecte en longitude plus que le
240	c contenu de la maille.
241	c cette partie est mal vectorisee.
242
243	c pas d'influence de la pression saturante (pour l'instant)
244
245	c calcul du nombre de maille sur lequel on advecte plus que la maille.
246
247	n0=0
248	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
249	DO l=1,llm
250	nl(l)=0
251	DO ij=ijb,ije
252	nl(l)=nl(l)+iadvplus(ij,l)
253	ENDDO
254	n0=n0+nl(l)
255	ENDDO
256	c$OMP END DO NOWAIT
257
258	cym ATTENTION ICI en OpenMP reduction pas forcement nécessaire
259	cym IF(n0.gt.1) THEN
260	cym IF(n0.gt.0) THEN
261	ccc PRINT*,'Nombre de points pour lesquels on advect plus que le'
262	ccc & ,'contenu de la maille : ',n0
263	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
264	DO l=1,llm
265	IF(nl(l).gt.0) THEN
266	iju=0
267	c indicage des mailles concernees par le traitement special
268	DO ij=ijb,ije
269	IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
270	iju=iju+1
271	indu(iju)=ij
272	ENDIF
273	ENDDO
274	niju=iju
275	c PRINT*,'niju,nl',niju,nl(l)
276
277	c traitement des mailles
278	DO iju=1,niju
279	ij=indu(iju)
280	j=(ij-1)/iip1+1
281	zu_m=u_m(ij,l)
282	u_mq(ij,l)=0.
283	IF(zu_m.gt.0.) THEN
284	ijq=ij
285	i=ijq-(j-1)*iip1
286	c accumulation pour les mailles completements advectees
287	do while(zu_m.gt.masse(ijq,l))
288	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
289	zu_m=zu_m-masse(ijq,l)
290	i=mod(i-2+iim,iim)+1
291	ijq=(j-1)*iip1+i
292	ENDDO
293	c ajout de la maille non completement advectee
294	u_mq(ij,l)=u_mq(ij,l)+zu_m*
295	& (q(ijq,l)+0.5(1.-zu_m/masse(ijq,l))dxq(ijq,l))
296	ELSE
297	ijq=ij+1
298	i=ijq-(j-1)*iip1
299	c accumulation pour les mailles completements advectees
300	do while(-zu_m.gt.masse(ijq,l))
301	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
302	zu_m=zu_m+masse(ijq,l)
303	i=mod(i,iim)+1
304	ijq=(j-1)*iip1+i
305	ENDDO
306	c ajout de la maille non completement advectee
307	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
308	& 0.5(1.+zu_m/masse(ijq,l))dxq(ijq,l))
309	ENDIF
310	ENDDO
311	ENDIF
312	ENDDO
313	c$OMP END DO NOWAIT
314	cym ENDIF ! n0.gt.0
315
316
317
318	c bouclage en latitude
319	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
320	DO l=1,llm
321	DO ij=ijb+iip1-1,ije,iip1
322	u_mq(ij,l)=u_mq(ij-iim,l)
323	ENDDO
324	ENDDO
325	c$OMP END DO NOWAIT
326
327	c calcul des tendances
328	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
329	DO l=1,llm
330	DO ij=ijb+1,ije
331	new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
332	q(ij,l)=(q(ij,l)*masse(ij,l)+
333	& u_mq(ij-1,l)-u_mq(ij,l))
334	& /new_m
335	masse(ij,l)=new_m
336	ENDDO
337	c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
338	DO ij=ijb+iip1-1,ije,iip1
339	q(ij-iim,l)=q(ij,l)
340	masse(ij-iim,l)=masse(ij,l)
341	ENDDO
342	ENDDO
343	c$OMP END DO NOWAIT
344	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
345	c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
346
347
348	RETURN
349	END
350	SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat)
351	c
352	c Auteurs: P.Le Van, F.Hourdin, F.Forget
353	c
354	c ********************************************************************
355	c Shema d'advection " pseudo amont " .
356	c ********************************************************************
357	c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
358	c qsat est un argument de sortie pour le s-pg ....
359	c
360	c
361	c --------------------------------------------------------------------
362	USE parallel
363	IMPLICIT NONE
364	c
365	#include "dimensions.h"
366	#include "paramet.h"
367	#include "logic.h"
368	#include "comvert.h"
369	#include "comconst.h"
370	#include "comgeom.h"
371	c
372	c
373	c Arguments:
374	c ----------
375	REAL masse(ijb_u:ije_u,llm),pente_max
376	REAL masse_adv_v( ijb_v:ije_v,llm)
377	REAL q(ijb_u:ije_u,llm)
378	REAL qsat(ijb_u:ije_u,llm)
379	c
380	c Local
381	c ---------
382	c
383	INTEGER i,ij,l
384	c
385	REAL airej2,airejjm,airescb(iim),airesch(iim)
386	REAL dyq(ijb_u:ije_u,llm),dyqv(ijb_v:ije_v)
387	REAL adyqv(ijb_v:ije_v),dyqmax(ijb_u:ije_u)
388	REAL qbyv(ijb_v:ije_v,llm)
389
390	REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
391	c REAL newq,oldmasse
392	Logical first
393	SAVE first
394	c$OMP THREADPRIVATE(first)
395	REAL convpn,convps,convmpn,convmps
396	REAL sinlon(iip1),sinlondlon(iip1)
397	REAL coslon(iip1),coslondlon(iip1)
398	SAVE sinlon,coslon,sinlondlon,coslondlon
399	SAVE airej2,airejjm
400	c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
401	c$OMP THREADPRIVATE(airej2,airejjm)
402	c
403	c
404	REAL SSUM
405
406	DATA first/.true./
407	INTEGER ijb,ije
408
409	IF(first) THEN
410	PRINT*,'Shema Amont nouveau appele dans Vanleer '
411	first=.false.
412	do i=2,iip1
413	coslon(i)=cos(rlonv(i))
414	sinlon(i)=sin(rlonv(i))
415	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
416	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
417	ENDDO
418	coslon(1)=coslon(iip1)
419	coslondlon(1)=coslondlon(iip1)
420	sinlon(1)=sinlon(iip1)
421	sinlondlon(1)=sinlondlon(iip1)
422	airej2 = SSUM( iim, aire(iip2), 1 )
423	airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
424	ENDIF
425
426	c
427
428	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
429	DO l = 1, llm
430	c
431	c --------------------------------
432	c CALCUL EN LATITUDE
433	c --------------------------------
434
435	c On commence par calculer la valeur du traceur moyenne sur le premier cercle
436	c de latitude autour du pole (qpns pour le pole nord et qpsn pour
437	c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
438
439	if (pole_nord) then
440	DO i = 1, iim
441	airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
442	ENDDO
443	qpns = SSUM( iim, airescb ,1 ) / airej2
444	endif
445
446	if (pole_sud) then
447	DO i = 1, iim
448	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
449	ENDDO
450	qpsn = SSUM( iim, airesch ,1 ) / airejjm
451	endif
452
453
454	c calcul des pentes aux points v
455
456	ijb=ij_begin-2*iip1
457	ije=ij_end+iip1
458	if (pole_nord) ijb=ij_begin
459	if (pole_sud) ije=ij_end-iip1
460
461	DO ij=ijb,ije
462	dyqv(ij)=q(ij,l)-q(ij+iip1,l)
463	adyqv(ij)=abs(dyqv(ij))
464	ENDDO
465
466
467	c calcul des pentes aux points scalaires
468
469	ijb=ij_begin-iip1
470	ije=ij_end+iip1
471	if (pole_nord) ijb=ij_begin+iip1
472	if (pole_sud) ije=ij_end-iip1
473
474	DO ij=ijb,ije
475	dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
476	dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
477	dyqmax(ij)=pente_max*dyqmax(ij)
478	ENDDO
479
480	IF (pole_nord) THEN
481
482	c calcul des pentes aux poles
483	DO ij=1,iip1
484	dyq(ij,l)=qpns-q(ij+iip1,l)
485	ENDDO
486
487	c filtrage de la derivee
488	dyn1=0.
489	dyn2=0.
490	DO ij=1,iim
491	dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
492	dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
493	ENDDO
494	DO ij=1,iip1
495	dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
496	ENDDO
497
498	c calcul des pentes limites aux poles
499	fn=1.
500	DO ij=1,iim
501	IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
502	fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
503	ENDIF
504	ENDDO
505
506	DO ij=1,iip1
507	dyq(ij,l)=fn*dyq(ij,l)
508	ENDDO
509
510	ENDIF
511
512	IF (pole_sud) THEN
513
514	DO ij=1,iip1
515	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
516	ENDDO
517
518	dys1=0.
519	dys2=0.
520
521	DO ij=1,iim
522	dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
523	dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
524	ENDDO
525
526	DO ij=1,iip1
527	dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
528	ENDDO
529
530	c calcul des pentes limites aux poles
531	fs=1.
532	DO ij=1,iim
533	IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
534	fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
535	ENDIF
536	ENDDO
537
538	DO ij=1,iip1
539	dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
540	ENDDO
541
542	ENDIF
543
544
545	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
546	C En memoire de dIFferents tests sur la
547	C limitation des pentes aux poles.
548	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
549	C PRINT*,dyq(1)
550	C PRINT*,dyqv(iip1+1)
551	C appn=abs(dyq(1)/dyqv(iip1+1))
552	C PRINT*,dyq(ip1jm+1)
553	C PRINT*,dyqv(ip1jm-iip1+1)
554	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
555	C DO ij=2,iim
556	C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
557	C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
558	C ENDDO
559	C appn=min(pente_max/appn,1.)
560	C apps=min(pente_max/apps,1.)
561	C
562	C
563	C cas ou on a un extremum au pole
564	C
565	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
566	C & appn=0.
567	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
568	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
569	C & apps=0.
570	C
571	C limitation des pentes aux poles
572	C DO ij=1,iip1
573	C dyq(ij)=appn*dyq(ij)
574	C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
575	C ENDDO
576	C
577	C test
578	C DO ij=1,iip1
579	C dyq(iip1+ij)=0.
580	C dyq(ip1jm+ij-iip1)=0.
581	C ENDDO
582	C DO ij=1,ip1jmp1
583	C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
584	C ENDDO
585	C
586	C changement 10 07 96
587	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
588	C & THEN
589	C DO ij=1,iip1
590	C dyqmax(ij)=0.
591	C ENDDO
592	C ELSE
593	C DO ij=1,iip1
594	C dyqmax(ij)=pente_max*abs(dyqv(ij))
595	C ENDDO
596	C ENDIF
597	C
598	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
599	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
600	C &THEN
601	C DO ij=ip1jm+1,ip1jmp1
602	C dyqmax(ij)=0.
603	C ENDDO
604	C ELSE
605	C DO ij=ip1jm+1,ip1jmp1
606	C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
607	C ENDDO
608	C ENDIF
609	C fin changement 10 07 96
610	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
611
612	c calcul des pentes limitees
613	ijb=ij_begin-iip1
614	ije=ij_end+iip1
615	if (pole_nord) ijb=ij_begin+iip1
616	if (pole_sud) ije=ij_end-iip1
617
618	DO ij=ijb,ije
619	IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
620	dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
621	ELSE
622	dyq(ij,l)=0.
623	ENDIF
624	ENDDO
625
626	ENDDO
627	c$OMP END DO NOWAIT
628
629	ijb=ij_begin-iip1
630	ije=ij_end
631	if (pole_nord) ijb=ij_begin
632	if (pole_sud) ije=ij_end-iip1
633
634	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
635	DO l=1,llm
636	DO ij=ijb,ije
637	IF( masse_adv_v(ij,l).GT.0. ) THEN
638	qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l ) +
639	, dyq(ij+iip1,l)0.5(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
640	ELSE
641	qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
642	, 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
643	ENDIF
644	qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
645	ENDDO
646	ENDDO
647	c$OMP END DO NOWAIT
648
649	ijb=ij_begin
650	ije=ij_end
651	if (pole_nord) ijb=ij_begin+iip1
652	if (pole_sud) ije=ij_end-iip1
653
654	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
655	DO l=1,llm
656	DO ij=ijb,ije
657	newmasse=masse(ij,l)
658	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
659	q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
660	& /newmasse
661	masse(ij,l)=newmasse
662	ENDDO
663	c.-. ancienne version
664
665	IF (pole_nord) THEN
666
667	convpn=SSUM(iim,qbyv(1,l),1)/apoln
668	convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
669	DO ij = 1,iip1
670	newmasse=masse(ij,l)+convmpn*aire(ij)
671	q(ij,l)=(q(ij,l)masse(ij,l)+convpnaire(ij))/
672	& newmasse
673	masse(ij,l)=newmasse
674	ENDDO
675
676	ENDIF
677
678	IF (pole_sud) THEN
679
680	convps = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
681	convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
682	DO ij = ip1jm+1,ip1jmp1
683	newmasse=masse(ij,l)+convmps*aire(ij)
684	q(ij,l)=(q(ij,l)masse(ij,l)+convpsaire(ij))/
685	& newmasse
686	masse(ij,l)=newmasse
687	ENDDO
688
689	ENDIF
690	c.-. fin ancienne version
691
692	c._. nouvelle version
693	c convpn=SSUM(iim,qbyv(1,l),1)
694	c convmpn=ssum(iim,masse_adv_v(1,l),1)
695	c oldmasse=ssum(iim,masse(1,l),1)
696	c newmasse=oldmasse+convmpn
697	c newq=(q(1,l)*oldmasse+convpn)/newmasse
698	c newmasse=newmasse/apoln
699	c DO ij = 1,iip1
700	c q(ij,l)=newq
701	c masse(ij,l)=newmasse*aire(ij)
702	c ENDDO
703	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
704	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
705	c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
706	c newmasse=oldmasse+convmps
707	c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
708	c newmasse=newmasse/apols
709	c DO ij = ip1jm+1,ip1jmp1
710	c q(ij,l)=newq
711	c masse(ij,l)=newmasse*aire(ij)
712	c ENDDO
713	c._. fin nouvelle version
714	ENDDO
715	c$OMP END DO NOWAIT
716	RETURN
717	END

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

Download in other formats:

Original Format