Context Navigation

vlspltqs_loc.F @ 2069

Last change on this file since 2069 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

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_lmdz
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_lmdz
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