Context Navigation

vlsplt_p.F @ 1024

Last change on this file since 1024 was 630, checked in by Laurent Fairhead, 20 years ago
Import d'une version parallele de la dynamique YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 28.8 KB

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

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

Download in other formats:

Original Format