Context Navigation

advn.F @ 1907

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

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

Download in other formats:

Original Format