Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

advn.F90 @ 5246

Last change on this file since 5246 was 5246, checked in by abarral, 23 hours ago
Convert fixed-form to free-form sources .F -> .{f,F}90 (WIP: some .F remain, will be handled in subsequent commits)
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: 23.8 KB

Line
1	!
2	! $Header$
3	!
4	SUBROUTINE advn(q,masse,w,pbaru,pbarv,pdt,mode)
5	!
6	! Auteur : F. Hourdin
7	!
8	! ********************************************************************
9	! Shema d'advection " pseudo amont " .
10	! ********************************************************************
11	! q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
12	!
13	! pbaru,pbarv,w flux de masse en u ,v ,w
14	! pdt pas de temps
15	!
16	! --------------------------------------------------------------------
17	IMPLICIT NONE
18	!
19	include "dimensions.h"
20	include "paramet.h"
21	include "comgeom.h"
22	include "iniprint.h"
23
24	!
25	! Arguments:
26	! ----------
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	!
33	! Local
34	! ---------
35	!
36	INTEGER :: i,ij,l,j,ii
37	integer :: ijlqmin,iqmin,jqmin,lqmin
38	integer :: ismin
39	!
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	! 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	! 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	! 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 SUBROUTINE advn
140
141	SUBROUTINE advnqx(q,qg,qd)
142	!
143	! Auteurs: Calcul des valeurs de q aux point u.
144	!
145	! --------------------------------------------------------------------
146	IMPLICIT NONE
147	!
148	INCLUDE "dimensions.h"
149	INCLUDE "paramet.h"
150	INCLUDE "iniprint.h"
151	!
152	!
153	! Arguments:
154	! ----------
155	real :: q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)
156	!
157	! Local
158	! ---------
159	!
160	INTEGER :: ij,l
161	!
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	! calcul des pentes en u:
171	! -----------------------
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	! 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	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	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 SUBROUTINE advnqx
258	SUBROUTINE advnqy(q,qs,qn)
259	!
260	! Auteurs: Calcul des valeurs de q aux point v.
261	!
262	! --------------------------------------------------------------------
263	IMPLICIT NONE
264	!
265	INCLUDE "dimensions.h"
266	INCLUDE "paramet.h"
267	INCLUDE "iniprint.h"
268	!
269	!
270	! Arguments:
271	! ----------
272	real :: q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)
273	!
274	! Local
275	! ---------
276	!
277	INTEGER :: ij,l
278	!
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	! calcul des pentes en u:
297	! -----------------------
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	! 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	! 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	! if (.not.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
335	! 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 SUBROUTINE advnqy
353
354	SUBROUTINE advnqz(q,qh,qb)
355	!
356	! Auteurs: Calcul des valeurs de q aux point v.
357	!
358	! --------------------------------------------------------------------
359	IMPLICIT NONE
360	!
361	INCLUDE "dimensions.h"
362	INCLUDE "paramet.h"
363	INCLUDE "iniprint.h"
364	!
365	!
366	! Arguments:
367	! ----------
368	real :: q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)
369	!
370	! Local
371	! ---------
372	!
373	INTEGER :: ij,l
374	!
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	! calcul des pentes en u:
385	! -----------------------
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	! 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	! 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	! do l=2,llm-1
450	! do ij=1,ip1jmp1
451	! if(extremum(ij,l)) then
452	! if (.not.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
453	! if (.not.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
454	! endif
455	! enddo
456	! 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 SUBROUTINE advnqz
469
470	SUBROUTINE advnx(q,qg,qd,masse,u_m,mode)
471	!
472	! Auteur : F. Hourdin
473	!
474	! ********************************************************************
475	! Shema d'advection " pseudo amont " .
476	! ********************************************************************
477	! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
478	!
479	!
480	! --------------------------------------------------------------------
481	IMPLICIT NONE
482	!
483	include "dimensions.h"
484	include "paramet.h"
485	include "iniprint.h"
486	!
487	!
488	! Arguments:
489	! ----------
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	!
495	! Local
496	! ---------
497	!
498	INTEGER :: i,j,ij,l,indu(ip1jmp1),niju,iju,ijq
499	integer :: n0,nl(llm)
500	!
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	! if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) then
522	! print*,'probleme au point ij=',ij,' l=',l
523	! print*,qd(ij,l),q(ij,l),qg(ij,l)
524	! qd(ij,l)=q(ij,l)
525	! qg(ij,l)=q(ij,l)
526	! 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	! if(.not.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .and.
531	! s zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) then
532	! print*,'probleme au point ij=',ij,' l=',l
533	! print*,'sigg=',zsigg(ij,l),' sigd=',zsigd(ij,l)
534	! print*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq
535	! stop
536	! 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	! 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	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	+(zsig-zsigp)(zq+zz(zqm-zq)) )
583	endif
584	endif
585	! if(zsig.lt.0.) then
586	! print*,'au point ij=',ij,' l=',l,' sig=',zsig
587	! stop
588	! 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	!=================================================================
600	! SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES
601	!=================================================================
602	! 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	! 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	! print*,'niju,nl',niju,nl(l)
632
633	! 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	! 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	! MODIFS SPECIFIQUES DU SCHEMA
650	! 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	-0.5zsig/zsigd(ijq,l)(qd(ijq,l)-q(ijq,l)))
655	else
656	! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
657	! 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	' l=',l
662	WRITE(lunout,*)'zz=',zz
663	stop
664	endif
665	u_mq(ij,l)=u_mq(ij,l)+masse(ijq,l)*( &
666	0.5(q(ijq,l)+qd(ijq,l))zsigd(ijq,l) &
667	+(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	! 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	! ajout de la maille non completement advectee
680	! 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	-0.5zsig/zsigg(ijq,l)(qg(ijq,l)-q(ijq,l)))
685	else
686	! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
687	! 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	,' l=',l
692	WRITE(lunout,*)'zz=',zz
693	stop
694	endif
695	u_mq(ij,l)=u_mq(ij,l)-masse(ijq,l)*( &
696	0.5(q(ijq,l)+qg(ijq,l))zsigg(ijq,l) &
697	+(zsig-zsigg(ijq,l))* &
698	(q(ijq,l)+zz*(qd(ijq,l)-q(ijq,l))) )
699	endif
700	! fin de la modif
701	endif
702	enddo
703	endif
704	enddo
705	endif ! n0.gt.0
706
707	! 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	!=================================================================
715	! CALCUL DE LA CONVERGENCE DES FLUX
716	!=================================================================
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	! 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 SUBROUTINE advnx
735	SUBROUTINE advny(q,qs,qn,masse,v_m)
736	!
737	! Auteur : F. Hourdin
738	!
739	! ********************************************************************
740	! Shema d'advection " pseudo amont " .
741	! ********************************************************************
742	! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
743	!
744	!
745	! --------------------------------------------------------------------
746	IMPLICIT NONE
747	!
748	INCLUDE "dimensions.h"
749	INCLUDE "paramet.h"
750	INCLUDE "comgeom.h"
751	INCLUDE "iniprint.h"
752	!
753	!
754	! Arguments:
755	! ----------
756	real :: masse(ip1jmp1,llm)
757	real :: v_m( ip1jm,llm )
758	real :: q(ip1jmp1,llm),qn(ip1jmp1,llm),qs(ip1jmp1,llm)
759	!
760	! Local
761	! ---------
762	!
763	INTEGER :: ij,l
764	!
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	! if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) then
783	! print*,'probleme au point ij=',ij,' l=',l,' advnqx'
784	! print*,qn(ij,l),q(ij,l),qs(ij,l)
785	! qn(ij,l)=q(ij,l)
786	! qs(ij,l)=q(ij,l)
787	! 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	! if(.not.(zsign(ij).ge.0..and.zsign(ij).le.1. .and.
792	! s zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) then
793	! print*,'probleme au point ij=',ij,' l=',l
794	! print*,'sigs=',zsigs(ij),' sign=',zsign(ij)
795	! stop
796	! endif
797	else
798	zsign(ij)=0.5
799	zsigs(ij)=0.5
800	endif
801	enddo
802
803	! 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	+(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	!.-. 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 SUBROUTINE advny
865	SUBROUTINE advnz(q,qh,qb,masse,w_m)
866	!
867	! Auteurs: F.Hourdin
868	!
869	! ********************************************************************
870	! Shema d'advection " pseudo amont " .
871	! b designe le bas et h le haut
872	! il y a une correspondance entre le b en z et le d en x
873	! ********************************************************************
874	!
875	!
876	! --------------------------------------------------------------------
877	IMPLICIT NONE
878	!
879	INCLUDE "dimensions.h"
880	INCLUDE "paramet.h"
881	INCLUDE "comgeom.h"
882	INCLUDE "iniprint.h"
883	!
884	!
885	! Arguments:
886	! ----------
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	!
892	! Local
893	! ---------
894	!
895	INTEGER :: ij,l
896	!
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	! if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) then
914	! print*,'probleme au point ij=',ij,' l=',l
915	! print*,qh(ij,l),q(ij,l),qb(ij,l)
916	! qh(ij,l)=q(ij,l)
917	! qb(ij,l)=q(ij,l)
918	! 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	! print*,'ok1'
932	! 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	+(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	! print*,'ok3'
976	RETURN
977	END SUBROUTINE advnz

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

Download in other formats:

Original Format