Context Navigation

advn.F @ 2597

Last change on this file since 2597 was 2597, checked in by Ehouarn Millour, 8 years ago
Cleanup in the dynamics: get rid of comconst.h, make it a module comconst_mod. EM
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 Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 25.8 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE advn(q,masse,w,pbaru,pbarv,pdt,mode)
5	c
6	c Auteur : F. Hourdin
7	c
8	c ********************************************************************
9	c Shema d'advection " pseudo amont " .
10	c ********************************************************************
11	c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
12	c
13	c pbaru,pbarv,w flux de masse en u ,v ,w
14	c pdt pas de temps
15	c
16	c --------------------------------------------------------------------
17	IMPLICIT NONE
18	c
19	#include "dimensions.h"
20	#include "paramet.h"
21	#include "logic.h"
22	#include "comvert.h"
23	#include "comgeom.h"
24	#include "iniprint.h"
25
26	c
27	c Arguments:
28	c ----------
29	integer mode
30	real masse(ip1jmp1,llm)
31	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
32	REAL q(ip1jmp1,llm)
33	REAL w(ip1jmp1,llm),pdt
34	c
35	c Local
36	c ---------
37	c
38	INTEGER i,ij,l,j,ii
39	integer ijlqmin,iqmin,jqmin,lqmin
40	integer ismin
41	c
42	real zm(ip1jmp1,llm),newmasse
43	real mu(ip1jmp1,llm)
44	real mv(ip1jm,llm)
45	real mw(ip1jmp1,llm+1)
46	real zq(ip1jmp1,llm),zz,qpn,qps
47	real zqg(ip1jmp1,llm),zqd(ip1jmp1,llm)
48	real zqs(ip1jmp1,llm),zqn(ip1jmp1,llm)
49	real zqh(ip1jmp1,llm),zqb(ip1jmp1,llm)
50	real temps0,temps1,temps2,temps3
51	real ztemps1,ztemps2,ztemps3,ssum
52	logical testcpu
53	save testcpu
54	save temps1,temps2,temps3
55	real zzpbar,zzw
56
57	#ifdef CRAY
58	real second
59	#endif
60
61	real qmin,qmax
62	data qmin,qmax/0.,1./
63	data testcpu/.false./
64	data temps1,temps2,temps3/0.,0.,0./
65
66	zzpbar = 0.5 * pdt
67	zzw = pdt
68
69	DO l=1,llm
70	DO ij = iip2,ip1jm
71	mu(ij,l)=pbaru(ij,l) * zzpbar
72	ENDDO
73	DO ij=1,ip1jm
74	mv(ij,l)=pbarv(ij,l) * zzpbar
75	ENDDO
76	DO ij=1,ip1jmp1
77	mw(ij,l)=w(ij,l) * zzw
78	ENDDO
79	ENDDO
80
81	DO ij=1,ip1jmp1
82	mw(ij,llm+1)=0.
83	ENDDO
84
85	do l=1,llm
86	qpn=0.
87	qps=0.
88	do ij=1,iim
89	qpn=qpn+q(ij,l)*masse(ij,l)
90	qps=qps+q(ip1jm+ij,l)*masse(ip1jm+ij,l)
91	enddo
92	qpn=qpn/ssum(iim,masse(1,l),1)
93	qps=qps/ssum(iim,masse(ip1jm+1,l),1)
94	do ij=1,iip1
95	q(ij,l)=qpn
96	q(ip1jm+ij,l)=qps
97	enddo
98	enddo
99
100	do ij=1,ip1jmp1
101	mw(ij,llm+1)=0.
102	enddo
103	do l=1,llm
104	do ij=1,ip1jmp1
105	zq(ij,l)=q(ij,l)
106	zm(ij,l)=masse(ij,l)
107	enddo
108	enddo
109
110	c call minmaxq(zq,qmin,qmax,'avant vlx ')
111	call advnqx(zq,zqg,zqd)
112	call advnx(zq,zqg,zqd,zm,mu,mode)
113	call advnqy(zq,zqs,zqn)
114	call advny(zq,zqs,zqn,zm,mv)
115	call advnqz(zq,zqh,zqb)
116	call advnz(zq,zqh,zqb,zm,mw)
117	c call vlz(zq,0.,zm,mw)
118	call advnqy(zq,zqs,zqn)
119	call advny(zq,zqs,zqn,zm,mv)
120	call advnqx(zq,zqg,zqd)
121	call advnx(zq,zqg,zqd,zm,mu,mode)
122	c call minmaxq(zq,qmin,qmax,'apres vlx ')
123
124	#ifdef CRAY
125	if(testcpu) then
126	ztemps1=second(0.)
127	temps1=temps1+ztemps1-ztemps2
128	print*,'VLSPLT X:',temps1,' Y:',temps2,' Z:',temps3
129	endif
130	#endif
131	do l=1,llm
132	do ij=1,ip1jmp1
133	q(ij,l)=zq(ij,l)
134	enddo
135	do ij=1,ip1jm+1,iip1
136	q(ij+iim,l)=q(ij,l)
137	enddo
138	enddo
139
140	RETURN
141	END
142
143	SUBROUTINE advnqx(q,qg,qd)
144	c
145	c Auteurs: Calcul des valeurs de q aux point u.
146	c
147	c --------------------------------------------------------------------
148	IMPLICIT NONE
149	c
150	#include "dimensions.h"
151	#include "paramet.h"
152	#include "iniprint.h"
153	c
154	c
155	c Arguments:
156	c ----------
157	real q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)
158	c
159	c Local
160	c ---------
161	c
162	INTEGER ij,l
163	c
164	real dxqu(ip1jmp1),zqu(ip1jmp1)
165	real zqmax(ip1jmp1),zqmin(ip1jmp1)
166	logical extremum(ip1jmp1)
167
168	integer mode
169	save mode
170	data mode/1/
171
172	c calcul des pentes en u:
173	c -----------------------
174	if (mode.eq.0) then
175	do l=1,llm
176	do ij=1,ip1jm
177	qd(ij,l)=q(ij,l)
178	qg(ij,l)=q(ij,l)
179	enddo
180	enddo
181	else
182	do l = 1, llm
183	do ij=iip2,ip1jm-1
184	dxqu(ij)=q(ij+1,l)-q(ij,l)
185	zqu(ij)=0.5*(q(ij+1,l)+q(ij,l))
186	enddo
187	do ij=iip1+iip1,ip1jm,iip1
188	dxqu(ij)=dxqu(ij-iim)
189	zqu(ij)=zqu(ij-iim)
190	enddo
191	do ij=iip2,ip1jm-1
192	zqu(ij)=zqu(ij)-dxqu(ij+1)/12.
193	enddo
194	do ij=iip1+iip1,ip1jm,iip1
195	zqu(ij)=zqu(ij-iim)
196	enddo
197	do ij=iip2+1,ip1jm
198	zqu(ij)=zqu(ij)+dxqu(ij-1)/12.
199	enddo
200	do ij=iip1+iip1,ip1jm,iip1
201	zqu(ij-iim)=zqu(ij)
202	enddo
203
204	c calcul des valeurs max et min acceptees aux interfaces
205
206	do ij=iip2,ip1jm-1
207	zqmax(ij)=max(q(ij+1,l),q(ij,l))
208	zqmin(ij)=min(q(ij+1,l),q(ij,l))
209	enddo
210	do ij=iip1+iip1,ip1jm,iip1
211	zqmax(ij)=zqmax(ij-iim)
212	zqmin(ij)=zqmin(ij-iim)
213	enddo
214	do ij=iip2+1,ip1jm
215	extremum(ij)=dxqu(ij)*dxqu(ij-1).le.0.
216	enddo
217	do ij=iip1+iip1,ip1jm,iip1
218	extremum(ij-iim)=extremum(ij)
219	enddo
220	do ij=iip2,ip1jm
221	zqu(ij)=min(max(zqmin(ij),zqu(ij)),zqmax(ij))
222	enddo
223	do ij=iip2+1,ip1jm
224	if(extremum(ij)) then
225	qg(ij,l)=q(ij,l)
226	qd(ij,l)=q(ij,l)
227	else
228	qd(ij,l)=zqu(ij)
229	qg(ij,l)=zqu(ij-1)
230	endif
231	enddo
232	do ij=iip1+iip1,ip1jm,iip1
233	qd(ij-iim,l)=qd(ij,l)
234	qg(ij-iim,l)=qg(ij,l)
235	enddo
236
237	goto 8888
238
239	do ij=iip2+1,ip1jm
240	if(extremum(ij).and..not.extremum(ij-1))
241	s qd(ij-1,l)=q(ij,l)
242	enddo
243
244	do ij=iip1+iip1,ip1jm,iip1
245	qd(ij-iim,l)=qd(ij,l)
246	enddo
247	do ij=iip2,ip1jm-1
248	if (extremum(ij).and..not.extremum(ij+1))
249	s qg(ij+1,l)=q(ij,l)
250	enddo
251
252	do ij=iip1+iip1,ip1jm,iip1
253	qg(ij,l)=qg(ij-iim,l)
254	enddo
255	8888 continue
256	enddo
257	endif
258	RETURN
259	END
260	SUBROUTINE advnqy(q,qs,qn)
261	c
262	c Auteurs: Calcul des valeurs de q aux point v.
263	c
264	c --------------------------------------------------------------------
265	IMPLICIT NONE
266	c
267	#include "dimensions.h"
268	#include "paramet.h"
269	#include "iniprint.h"
270	c
271	c
272	c Arguments:
273	c ----------
274	real q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)
275	c
276	c Local
277	c ---------
278	c
279	INTEGER ij,l
280	c
281	real dyqv(ip1jm),zqv(ip1jm,llm)
282	real zqmax(ip1jm),zqmin(ip1jm)
283	logical extremum(ip1jmp1)
284
285	integer mode
286	save mode
287	data mode/1/
288
289	if (mode.eq.0) then
290	do l=1,llm
291	do ij=1,ip1jmp1
292	qn(ij,l)=q(ij,l)
293	qs(ij,l)=q(ij,l)
294	enddo
295	enddo
296	else
297
298	c calcul des pentes en u:
299	c -----------------------
300	do l = 1, llm
301	do ij=1,ip1jm
302	dyqv(ij)=q(ij,l)-q(ij+iip1,l)
303	enddo
304
305	do ij=iip2,ip1jm-iip1
306	zqv(ij,l)=0.5*(q(ij+iip1,l)+q(ij,l))
307	zqv(ij,l)=zqv(ij,l)+(dyqv(ij+iip1)-dyqv(ij-iip1))/12.
308	enddo
309
310	do ij=iip2,ip1jm
311	extremum(ij)=dyqv(ij)*dyqv(ij-iip1).le.0.
312	enddo
313
314	c Pas de pentes aux poles
315	do ij=1,iip1
316	zqv(ij,l)=q(ij,l)
317	zqv(ip1jm-iip1+ij,l)=q(ip1jm+ij,l)
318	extremum(ij)=.true.
319	extremum(ip1jmp1-iip1+ij)=.true.
320	enddo
321
322	c calcul des valeurs max et min acceptees aux interfaces
323	do ij=1,ip1jm
324	zqmax(ij)=max(q(ij+iip1,l),q(ij,l))
325	zqmin(ij)=min(q(ij+iip1,l),q(ij,l))
326	enddo
327
328	do ij=1,ip1jm
329	zqv(ij,l)=min(max(zqmin(ij),zqv(ij,l)),zqmax(ij))
330	enddo
331
332	do ij=iip2,ip1jm
333	if(extremum(ij)) then
334	qs(ij,l)=q(ij,l)
335	qn(ij,l)=q(ij,l)
336	c if (.not.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
337	c if (.not.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l)
338	else
339	qs(ij,l)=zqv(ij,l)
340	qn(ij,l)=zqv(ij-iip1,l)
341	endif
342	enddo
343
344	do ij=1,iip1
345	qs(ij,l)=q(ij,l)
346	qn(ij,l)=q(ij,l)
347	qs(ip1jm+ij,l)=q(ip1jm+ij,l)
348	qn(ip1jm+ij,l)=q(ip1jm+ij,l)
349	enddo
350
351	enddo
352	endif
353	RETURN
354	END
355
356	SUBROUTINE advnqz(q,qh,qb)
357	c
358	c Auteurs: Calcul des valeurs de q aux point v.
359	c
360	c --------------------------------------------------------------------
361	IMPLICIT NONE
362	c
363	#include "dimensions.h"
364	#include "paramet.h"
365	#include "iniprint.h"
366	c
367	c
368	c Arguments:
369	c ----------
370	real q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)
371	c
372	c Local
373	c ---------
374	c
375	INTEGER ij,l
376	c
377	real dzqw(ip1jmp1,llm+1),zqw(ip1jmp1,llm+1)
378	real zqmax(ip1jmp1,llm),zqmin(ip1jmp1,llm)
379	logical extremum(ip1jmp1,llm)
380
381	integer mode
382	save mode
383
384	data mode/1/
385
386	c calcul des pentes en u:
387	c -----------------------
388
389	if (mode.eq.0) then
390	do l=1,llm
391	do ij=1,ip1jmp1
392	qb(ij,l)=q(ij,l)
393	qh(ij,l)=q(ij,l)
394	enddo
395	enddo
396	else
397	do l = 2, llm
398	do ij=1,ip1jmp1
399	dzqw(ij,l)=q(ij,l-1)-q(ij,l)
400	zqw(ij,l)=0.5*(q(ij,l-1)+q(ij,l))
401	enddo
402	enddo
403	do ij=1,ip1jmp1
404	dzqw(ij,1)=0.
405	dzqw(ij,llm+1)=0.
406	enddo
407	do l=2,llm
408	do ij=1,ip1jmp1
409	zqw(ij,l)=zqw(ij,l)+(dzqw(ij,l+1)-dzqw(ij,l-1))/12.
410	enddo
411	enddo
412	do l=2,llm-1
413	do ij=1,ip1jmp1
414	extremum(ij,l)=dzqw(ij,l)*dzqw(ij,l+1).le.0.
415	enddo
416	enddo
417
418	c Pas de pentes en bas et en haut
419	do ij=1,ip1jmp1
420	zqw(ij,2)=q(ij,1)
421	zqw(ij,llm)=q(ij,llm)
422	extremum(ij,1)=.true.
423	extremum(ij,llm)=.true.
424	enddo
425
426	c calcul des valeurs max et min acceptees aux interfaces
427	do l=2,llm
428	do ij=1,ip1jmp1
429	zqmax(ij,l)=max(q(ij,l-1),q(ij,l))
430	zqmin(ij,l)=min(q(ij,l-1),q(ij,l))
431	enddo
432	enddo
433
434	do l=2,llm
435	do ij=1,ip1jmp1
436	zqw(ij,l)=min(max(zqmin(ij,l),zqw(ij,l)),zqmax(ij,l))
437	enddo
438	enddo
439
440	do l=2,llm-1
441	do ij=1,ip1jmp1
442	if(extremum(ij,l)) then
443	qh(ij,l)=q(ij,l)
444	qb(ij,l)=q(ij,l)
445	else
446	qh(ij,l)=zqw(ij,l+1)
447	qb(ij,l)=zqw(ij,l)
448	endif
449	enddo
450	enddo
451	c do l=2,llm-1
452	c do ij=1,ip1jmp1
453	c if(extremum(ij,l)) then
454	c if (.not.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
455	c if (.not.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
456	c endif
457	c enddo
458	c enddo
459
460	do ij=1,ip1jmp1
461	qb(ij,1)=q(ij,1)
462	qh(ij,1)=q(ij,1)
463	qb(ij,llm)=q(ij,llm)
464	qh(ij,llm)=q(ij,llm)
465	enddo
466
467	endif
468
469	RETURN
470	END
471
472	SUBROUTINE advnx(q,qg,qd,masse,u_m,mode)
473	c
474	c Auteur : F. Hourdin
475	c
476	c ********************************************************************
477	c Shema d'advection " pseudo amont " .
478	c ********************************************************************
479	c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
480	c
481	c
482	c --------------------------------------------------------------------
483	IMPLICIT NONE
484	c
485	#include "dimensions.h"
486	#include "paramet.h"
487	#include "logic.h"
488	#include "comvert.h"
489	#include "iniprint.h"
490	c
491	c
492	c Arguments:
493	c ----------
494	integer mode
495	real masse(ip1jmp1,llm)
496	real u_m( ip1jmp1,llm )
497	real q(ip1jmp1,llm),qd(ip1jmp1,llm),qg(ip1jmp1,llm)
498	c
499	c Local
500	c ---------
501	c
502	INTEGER i,j,ij,l,indu(ip1jmp1),niju,iju,ijq
503	integer n0,nl(llm)
504	c
505	real new_m,zu_m,zdq,zz
506	real zsigg(ip1jmp1,llm),zsigd(ip1jmp1,llm),zsig
507	real u_mq(ip1jmp1,llm)
508
509	real zm,zq,zsigm,zsigp,zqm,zqp,zu
510
511	logical ladvplus(ip1jmp1,llm)
512
513	real prec
514	save prec
515
516	#ifdef CRAY
517	data prec/1.e-24/
518	#else
519	data prec/1.e-15/
520	#endif
521
522	do l=1,llm
523	do ij=iip2,ip1jm
524	zdq=qd(ij,l)-qg(ij,l)
525	c if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) then
526	c print*,'probleme au point ij=',ij,' l=',l
527	c print*,qd(ij,l),q(ij,l),qg(ij,l)
528	c qd(ij,l)=q(ij,l)
529	c qg(ij,l)=q(ij,l)
530	c endif
531	if(abs(zdq).gt.prec) then
532	zsigd(ij,l)=(q(ij,l)-qg(ij,l))/zdq
533	zsigg(ij,l)=1.-zsigd(ij,l)
534	c if(.not.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .and.
535	c s zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) then
536	c print*,'probleme au point ij=',ij,' l=',l
537	c print*,'sigg=',zsigg(ij,l),' sigd=',zsigd(ij,l)
538	c print*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq
539	c stop
540	c endif
541	else
542	zsigd(ij,l)=0.5
543	zsigg(ij,l)=0.5
544	qd(ij,l)=q(ij,l)
545	qg(ij,l)=q(ij,l)
546	endif
547	enddo
548	enddo
549
550	c calcul de la pente maximum dans la maille en valeur absolue
551
552	do l=1,llm
553	do ij=iip2,ip1jm-1
554	if (u_m(ij,l).ge.0.) then
555	zsigp=zsigd(ij,l)
556	zsigm=zsigg(ij,l)
557	zqp=qd(ij,l)
558	zqm=qg(ij,l)
559	zm=masse(ij,l)
560	zq=q(ij,l)
561	else
562	zsigm=zsigd(ij+1,l)
563	zsigp=zsigg(ij+1,l)
564	zqm=qd(ij+1,l)
565	zqp=qg(ij+1,l)
566	zm=masse(ij+1,l)
567	zq=q(ij+1,l)
568	endif
569	zu=abs(u_m(ij,l))
570	ladvplus(ij,l)=zu.gt.zm
571	zsig=zu/zm
572	if(zsig.eq.0.) zsigp=0.1
573	if (mode.eq.1) then
574	if (zsig.le.zsigp) then
575	u_mq(ij,l)=u_m(ij,l)*zqp
576	else if (mode.eq.1) then
577	u_mq(ij,l)=
578	s sign(zm,u_m(ij,l))(zsigpzqp+(zsig-zsigp)*zqm)
579	endif
580	else
581	if (zsig.le.zsigp) then
582	u_mq(ij,l)=u_m(ij,l)(zqp-0.5zsig/zsigp*(zqp-zq))
583	else
584	zz=0.5*(zsig-zsigp)/zsigm
585	u_mq(ij,l)=sign(zm,u_m(ij,l))( 0.5(zq+zqp)*zsigp
586	s +(zsig-zsigp)(zq+zz(zqm-zq)) )
587	endif
588	endif
589	c if(zsig.lt.0.) then
590	c print*,'au point ij=',ij,' l=',l,' sig=',zsig
591	c stop
592	c endif
593	enddo
594	enddo
595
596	do l=1,llm
597	do ij=iip1+iip1,ip1jm,iip1
598	u_mq(ij,l)=u_mq(ij-iim,l)
599	ladvplus(ij,l)=ladvplus(ij-iim,l)
600	enddo
601	enddo
602
603	c=================================================================
604	C SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES
605	c=================================================================
606	c tris des regions a traiter
607	n0=0
608	do l=1,llm
609	nl(l)=0
610	do ij=iip2,ip1jm
611	if(ladvplus(ij,l)) then
612	nl(l)=nl(l)+1
613	u_mq(ij,l)=0.
614	endif
615	enddo
616	n0=n0+nl(l)
617	enddo
618
619	if(n0.gt.1) then
620	IF (prt_level > 9) WRITE(lunout,*)
621	& 'Nombre de points pour lesquels on advect plus que le'
622	& ,'contenu de la maille : ',n0
623
624	do l=1,llm
625	if(nl(l).gt.0) then
626	iju=0
627	c indicage des mailles concernees par le traitement special
628	do ij=iip2,ip1jm
629	if(ladvplus(ij,l).and.mod(ij,iip1).ne.0) then
630	iju=iju+1
631	indu(iju)=ij
632	endif
633	enddo
634	niju=iju
635	c print*,'niju,nl',niju,nl(l)
636
637	c traitement des mailles
638	do iju=1,niju
639	ij=indu(iju)
640	j=(ij-1)/iip1+1
641	zu_m=u_m(ij,l)
642	u_mq(ij,l)=0.
643	if(zu_m.gt.0.) then
644	ijq=ij
645	i=ijq-(j-1)*iip1
646	c accumulation pour les mailles completements advectees
647	do while(zu_m.gt.masse(ijq,l))
648	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
649	zu_m=zu_m-masse(ijq,l)
650	i=mod(i-2+iim,iim)+1
651	ijq=(j-1)*iip1+i
652	enddo
653	c MODIFS SPECIFIQUES DU SCHEMA
654	c ajout de la maille non completement advectee
655	zsig=zu_m/masse(ijq,l)
656	if(zsig.le.zsigd(ijq,l)) then
657	u_mq(ij,l)=u_mq(ij,l)+zu_m*(qd(ijq,l)
658	s -0.5zsig/zsigd(ijq,l)(qd(ijq,l)-q(ijq,l)))
659	else
660	c u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
661	c goto 8888
662	zz=0.5*(zsig-zsigd(ijq,l))/zsigg(ijq,l)
663	if(.not.(zz.gt.0..and.zz.le.0.5)) then
664	WRITE(lunout,*)'probleme2 au point ij=',ij,
665	s ' l=',l
666	WRITE(lunout,*)'zz=',zz
667	stop
668	endif
669	u_mq(ij,l)=u_mq(ij,l)+masse(ijq,l)*(
670	s 0.5(q(ijq,l)+qd(ijq,l))zsigd(ijq,l)
671	s +(zsig-zsigd(ijq,l))(q(ijq,l)+zz(qg(ijq,l)-q(ijq,l))) )
672	endif
673	else
674	ijq=ij+1
675	i=ijq-(j-1)*iip1
676	c accumulation pour les mailles completements advectees
677	do while(-zu_m.gt.masse(ijq,l))
678	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
679	zu_m=zu_m+masse(ijq,l)
680	i=mod(i,iim)+1
681	ijq=(j-1)*iip1+i
682	enddo
683	c ajout de la maille non completement advectee
684	c 2eme MODIF SPECIFIQUE
685	zsig=-zu_m/masse(ij+1,l)
686	if(zsig.le.zsigg(ijq,l)) then
687	u_mq(ij,l)=u_mq(ij,l)+zu_m*(qg(ijq,l)
688	s -0.5zsig/zsigg(ijq,l)(qg(ijq,l)-q(ijq,l)))
689	else
690	c u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
691	c goto 9999
692	zz=0.5*(zsig-zsigg(ijq,l))/zsigd(ijq,l)
693	if(.not.(zz.gt.0..and.zz.le.0.5)) then
694	WRITE(lunout,*)'probleme22 au point ij=',ij
695	s ,' l=',l
696	WRITE(lunout,*)'zz=',zz
697	stop
698	endif
699	u_mq(ij,l)=u_mq(ij,l)-masse(ijq,l)*(
700	s 0.5(q(ijq,l)+qg(ijq,l))zsigg(ijq,l)
701	s +(zsig-zsigg(ijq,l))*
702	s (q(ijq,l)+zz*(qd(ijq,l)-q(ijq,l))) )
703	endif
704	c fin de la modif
705	endif
706	enddo
707	endif
708	enddo
709	endif ! n0.gt.0
710
711	c bouclage en latitude
712	do l=1,llm
713	do ij=iip1+iip1,ip1jm,iip1
714	u_mq(ij,l)=u_mq(ij-iim,l)
715	enddo
716	enddo
717
718	c=================================================================
719	c CALCUL DE LA CONVERGENCE DES FLUX
720	c=================================================================
721
722	do l=1,llm
723	do ij=iip2+1,ip1jm
724	new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
725	q(ij,l)=(q(ij,l)*masse(ij,l)+
726	& u_mq(ij-1,l)-u_mq(ij,l))
727	& /new_m
728	masse(ij,l)=new_m
729	enddo
730	c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
731	do ij=iip1+iip1,ip1jm,iip1
732	q(ij-iim,l)=q(ij,l)
733	masse(ij-iim,l)=masse(ij,l)
734	enddo
735	enddo
736
737	RETURN
738	END
739	SUBROUTINE advny(q,qs,qn,masse,v_m)
740	c
741	c Auteur : F. Hourdin
742	c
743	c ********************************************************************
744	c Shema d'advection " pseudo amont " .
745	c ********************************************************************
746	c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
747	c
748	c
749	c --------------------------------------------------------------------
750	IMPLICIT NONE
751	c
752	#include "dimensions.h"
753	#include "paramet.h"
754	#include "comgeom.h"
755	#include "iniprint.h"
756	c
757	c
758	c Arguments:
759	c ----------
760	real masse(ip1jmp1,llm)
761	real v_m( ip1jm,llm )
762	real q(ip1jmp1,llm),qn(ip1jmp1,llm),qs(ip1jmp1,llm)
763	c
764	c Local
765	c ---------
766	c
767	INTEGER ij,l
768	c
769	real new_m,zdq,zz
770	real zsigs(ip1jmp1),zsign(ip1jmp1),zsig
771	real v_mq(ip1jm,llm)
772	real convpn,convps,convmpn,convmps,massen,masses
773	real zm,zq,zsigm,zsigp,zqm,zqp
774	real ssum
775	real prec
776	save prec
777
778	#ifdef CRAY
779	data prec/1.e-24/
780	#else
781	data prec/1.e-15/
782	#endif
783	do l=1,llm
784	do ij=1,ip1jmp1
785	zdq=qn(ij,l)-qs(ij,l)
786	c if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) then
787	c print*,'probleme au point ij=',ij,' l=',l,' advnqx'
788	c print*,qn(ij,l),q(ij,l),qs(ij,l)
789	c qn(ij,l)=q(ij,l)
790	c qs(ij,l)=q(ij,l)
791	c endif
792	if(abs(zdq).gt.prec) then
793	zsign(ij)=(q(ij,l)-qs(ij,l))/zdq
794	zsigs(ij)=1.-zsign(ij)
795	c if(.not.(zsign(ij).ge.0..and.zsign(ij).le.1. .and.
796	c s zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) then
797	c print*,'probleme au point ij=',ij,' l=',l
798	c print*,'sigs=',zsigs(ij),' sign=',zsign(ij)
799	c stop
800	c endif
801	else
802	zsign(ij)=0.5
803	zsigs(ij)=0.5
804	endif
805	enddo
806
807	c calcul de la pente maximum dans la maille en valeur absolue
808
809	do ij=1,ip1jm
810	if (v_m(ij,l).ge.0.) then
811	zsigp=zsign(ij+iip1)
812	zsigm=zsigs(ij+iip1)
813	zqp=qn(ij+iip1,l)
814	zqm=qs(ij+iip1,l)
815	zm=masse(ij+iip1,l)
816	zq=q(ij+iip1,l)
817	else
818	zsigm=zsign(ij)
819	zsigp=zsigs(ij)
820	zqm=qn(ij,l)
821	zqp=qs(ij,l)
822	zm=masse(ij,l)
823	zq=q(ij,l)
824	endif
825	zsig=abs(v_m(ij,l))/zm
826	if(zsig.eq.0.) zsigp=0.1
827	if (zsig.le.zsigp) then
828	v_mq(ij,l)=v_m(ij,l)(zqp-0.5zsig/zsigp*(zqp-zq))
829	else
830	zz=0.5*(zsig-zsigp)/zsigm
831	v_mq(ij,l)=sign(zm,v_m(ij,l))( 0.5(zq+zqp)*zsigp
832	s +(zsig-zsigp)(zq+zz(zqm-zq)) )
833	endif
834	enddo
835	enddo
836
837	do l=1,llm
838	do ij=iip2,ip1jm
839	new_m=masse(ij,l)
840	& +v_m(ij,l)-v_m(ij-iip1,l)
841	q(ij,l)=(q(ij,l)*masse(ij,l)+v_mq(ij,l)-v_mq(ij-iip1,l))
842	& /new_m
843	masse(ij,l)=new_m
844	enddo
845	c.-. ancienne version
846	convpn=SSUM(iim,v_mq(1,l),1)
847	convmpn=ssum(iim,v_m(1,l),1)
848	massen=ssum(iim,masse(1,l),1)
849	new_m=massen+convmpn
850	q(1,l)=(q(1,l)*massen+convpn)/new_m
851	do ij = 1,iip1
852	q(ij,l)=q(1,l)
853	masse(ij,l)=new_m*aire(ij)/apoln
854	enddo
855
856	convps=-SSUM(iim,v_mq(ip1jm-iim,l),1)
857	convmps=-ssum(iim,v_m(ip1jm-iim,l),1)
858	masses=ssum(iim,masse(ip1jm+1,l),1)
859	new_m=masses+convmps
860	q(ip1jm+1,l)=(q(ip1jm+1,l)*masses+convps)/new_m
861	do ij = ip1jm+1,ip1jmp1
862	q(ij,l)=q(ip1jm+1,l)
863	masse(ij,l)=new_m*aire(ij)/apols
864	enddo
865	enddo
866
867	RETURN
868	END
869	SUBROUTINE advnz(q,qh,qb,masse,w_m)
870	c
871	c Auteurs: F.Hourdin
872	c
873	c ********************************************************************
874	c Shema d'advection " pseudo amont " .
875	c b designe le bas et h le haut
876	c il y a une correspondance entre le b en z et le d en x
877	c ********************************************************************
878	c
879	c
880	c --------------------------------------------------------------------
881	IMPLICIT NONE
882	c
883	#include "dimensions.h"
884	#include "paramet.h"
885	#include "comgeom.h"
886	#include "iniprint.h"
887	c
888	c
889	c Arguments:
890	c ----------
891	real masse(ip1jmp1,llm)
892	real w_m( ip1jmp1,llm+1)
893	real q(ip1jmp1,llm),qb(ip1jmp1,llm),qh(ip1jmp1,llm)
894
895	c
896	c Local
897	c ---------
898	c
899	INTEGER ij,l
900	c
901	real new_m,zdq,zz
902	real zsigh(ip1jmp1,llm),zsigb(ip1jmp1,llm),zsig
903	real w_mq(ip1jmp1,llm+1)
904	real zm,zq,zsigm,zsigp,zqm,zqp
905	real prec
906	save prec
907
908	#ifdef CRAY
909	data prec/1.e-24/
910	#else
911	data prec/1.e-13/
912	#endif
913
914	do l=1,llm
915	do ij=1,ip1jmp1
916	zdq=qb(ij,l)-qh(ij,l)
917	c if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) then
918	c print*,'probleme au point ij=',ij,' l=',l
919	c print*,qh(ij,l),q(ij,l),qb(ij,l)
920	c qh(ij,l)=q(ij,l)
921	c qb(ij,l)=q(ij,l)
922	c endif
923
924	if(abs(zdq).gt.prec) then
925	zsigb(ij,l)=(q(ij,l)-qh(ij,l))/zdq
926	zsigh(ij,l)=1.-zsigb(ij,l)
927	zsigb(ij,l)=min(max(zsigb(ij,l),0.),1.)
928	else
929	zsigb(ij,l)=0.5
930	zsigh(ij,l)=0.5
931	endif
932	enddo
933	enddo
934
935	c print*,'ok1'
936	c calcul de la pente maximum dans la maille en valeur absolue
937	do l=2,llm
938	do ij=1,ip1jmp1
939	if (w_m(ij,l).ge.0.) then
940	zsigp=zsigb(ij,l)
941	zsigm=zsigh(ij,l)
942	zqp=qb(ij,l)
943	zqm=qh(ij,l)
944	zm=masse(ij,l)
945	zq=q(ij,l)
946	else
947	zsigm=zsigb(ij,l-1)
948	zsigp=zsigh(ij,l-1)
949	zqm=qb(ij,l-1)
950	zqp=qh(ij,l-1)
951	zm=masse(ij,l-1)
952	zq=q(ij,l-1)
953	endif
954	zsig=abs(w_m(ij,l))/zm
955	if(zsig.eq.0.) zsigp=0.1
956	if (zsig.le.zsigp) then
957	w_mq(ij,l)=w_m(ij,l)(zqp-0.5zsig/zsigp*(zqp-zq))
958	else
959	zz=0.5*(zsig-zsigp)/zsigm
960	w_mq(ij,l)=sign(zm,w_m(ij,l))( 0.5(zq+zqp)*zsigp
961	s +(zsig-zsigp)(zq+zz(zqm-zq)) )
962	endif
963	enddo
964	enddo
965
966	do ij=1,ip1jmp1
967	w_mq(ij,llm+1)=0.
968	w_mq(ij,1)=0.
969	enddo
970
971	do l=1,llm
972	do ij=1,ip1jmp1
973	new_m=masse(ij,l)+w_m(ij,l+1)-w_m(ij,l)
974	q(ij,l)=(q(ij,l)*masse(ij,l)+w_mq(ij,l+1)-w_mq(ij,l))
975	& /new_m
976	masse(ij,l)=new_m
977	enddo
978	enddo
979	c print*,'ok3'
980	RETURN
981	END

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

Download in other formats:

Original Format