Context Navigation

advn.F @ 5313

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

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

Download in other formats:

Original Format