Context Navigation

vlsplt.F @ 5415

Last change on this file since 5415 was 595, checked in by Laurent Fairhead, 20 years ago
Si il n'y avait qu'un point ou on advectait plus que le contenu de la maille, il n'etait jamais traite YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 25.1 KB

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

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

Download in other formats:

Original Format