Context Navigation

advn.F @ 5099

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

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

Download in other formats:

Original Format