Context Navigation

vlsplt_loc.F @ 5440

Last change on this file since 5440 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: 24.2 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_lmdz
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_lmdz
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_lmdz
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