Context Navigation

vlsplt.F @ 5360

Last change on this file since 5360 was 1907, checked in by lguez, 11 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
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 25.2 KB

Line
1	c
2	c $Id: vlsplt.F 1907 2013-11-26 13:10:46Z abarral $
3	c
4
5	SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt)
6	c
7	c Auteurs: P.Le Van, F.Hourdin, F.Forget
8	c
9	c ********************************************************************
10	c Shema d'advection " pseudo amont " .
11	c ********************************************************************
12	c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
13	c
14	c pente_max facteur de limitation des pentes: 2 en general
15	c 0 pour un schema amont
16	c pbaru,pbarv,w flux de masse en u ,v ,w
17	c pdt pas de temps
18	c
19	c --------------------------------------------------------------------
20	IMPLICIT NONE
21	c
22	#include "dimensions.h"
23	#include "paramet.h"
24	#include "logic.h"
25	#include "comvert.h"
26	#include "comconst.h"
27
28	c
29	c Arguments:
30	c ----------
31	REAL masse(ip1jmp1,llm),pente_max
32	c REAL masse(iip1,jjp1,llm),pente_max
33	REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
34	REAL q(ip1jmp1,llm)
35	c REAL q(iip1,jjp1,llm)
36	REAL w(ip1jmp1,llm),pdt
37	c
38	c Local
39	c ---------
40	c
41	INTEGER i,ij,l,j,ii
42	INTEGER ijlqmin,iqmin,jqmin,lqmin
43	c
44	REAL zm(ip1jmp1,llm),newmasse
45	REAL mu(ip1jmp1,llm)
46	REAL mv(ip1jm,llm)
47	REAL mw(ip1jmp1,llm+1)
48	REAL zq(ip1jmp1,llm),zz
49	REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
50	REAL second,temps0,temps1,temps2,temps3
51	REAL ztemps1,ztemps2,ztemps3
52	REAL zzpbar, zzw
53	LOGICAL testcpu
54	SAVE testcpu
55	SAVE temps1,temps2,temps3
56	INTEGER iminn,imaxx
57
58	REAL qmin,qmax
59	DATA qmin,qmax/0.,1.e33/
60	DATA testcpu/.false./
61	DATA temps1,temps2,temps3/0.,0.,0./
62
63
64	zzpbar = 0.5 * pdt
65	zzw = pdt
66	DO l=1,llm
67	DO ij = iip2,ip1jm
68	mu(ij,l)=pbaru(ij,l) * zzpbar
69	ENDDO
70	DO ij=1,ip1jm
71	mv(ij,l)=pbarv(ij,l) * zzpbar
72	ENDDO
73	DO ij=1,ip1jmp1
74	mw(ij,l)=w(ij,l) * zzw
75	ENDDO
76	ENDDO
77
78	DO ij=1,ip1jmp1
79	mw(ij,llm+1)=0.
80	ENDDO
81
82	CALL SCOPY(ijp1llm,q,1,zq,1)
83	CALL SCOPY(ijp1llm,masse,1,zm,1)
84
85	cprint*,'Entree vlx1'
86	c call minmaxq(zq,qmin,qmax,'avant vlx ')
87	call vlx(zq,pente_max,zm,mu)
88	cprint*,'Sortie vlx1'
89	c call minmaxq(zq,qmin,qmax,'apres vlx1 ')
90
91	c print*,'Entree vly1'
92	call vly(zq,pente_max,zm,mv)
93	c call minmaxq(zq,qmin,qmax,'apres vly1 ')
94	cprint*,'Sortie vly1'
95	call vlz(zq,pente_max,zm,mw)
96	c call minmaxq(zq,qmin,qmax,'apres vlz ')
97
98
99	call vly(zq,pente_max,zm,mv)
100	c call minmaxq(zq,qmin,qmax,'apres vly ')
101
102
103	call vlx(zq,pente_max,zm,mu)
104	c call minmaxq(zq,qmin,qmax,'apres vlx2 ')
105
106
107	DO l=1,llm
108	DO ij=1,ip1jmp1
109	q(ij,l)=zq(ij,l)
110	ENDDO
111	DO ij=1,ip1jm+1,iip1
112	q(ij+iim,l)=q(ij,l)
113	ENDDO
114	ENDDO
115
116	RETURN
117	END
118	SUBROUTINE vlx(q,pente_max,masse,u_m)
119
120	c Auteurs: P.Le Van, F.Hourdin, F.Forget
121	c
122	c ********************************************************************
123	c Shema d'advection " pseudo amont " .
124	c ********************************************************************
125	c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
126	c
127	c
128	c --------------------------------------------------------------------
129	IMPLICIT NONE
130	c
131	include "dimensions.h"
132	include "paramet.h"
133	include "logic.h"
134	include "comvert.h"
135	include "comconst.h"
136	include "iniprint.h"
137	c
138	c
139	c Arguments:
140	c ----------
141	REAL masse(ip1jmp1,llm),pente_max
142	REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
143	REAL q(ip1jmp1,llm)
144	REAL w(ip1jmp1,llm)
145	c
146	c Local
147	c ---------
148	c
149	INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
150	INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
151	c
152	REAL new_m,zu_m,zdum(ip1jmp1,llm)
153	REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
154	REAL zz(ip1jmp1)
155	REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
156	REAL u_mq(ip1jmp1,llm)
157
158	Logical extremum,first,testcpu
159	SAVE first,testcpu
160
161	REAL SSUM
162	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
163	SAVE temps0,temps1,temps2,temps3,temps4,temps5
164
165	REAL z1,z2,z3
166
167	DATA first,testcpu/.true.,.false./
168
169	IF(first) THEN
170	temps1=0.
171	temps2=0.
172	temps3=0.
173	temps4=0.
174	temps5=0.
175	first=.false.
176	ENDIF
177
178	c calcul de la pente a droite et a gauche de la maille
179
180
181	IF (pente_max.gt.-1.e-5) THEN
182	c IF (pente_max.gt.10) THEN
183
184	c calcul des pentes avec limitation, Van Leer scheme I:
185	c -----------------------------------------------------
186
187	c calcul de la pente aux points u
188	DO l = 1, llm
189	DO ij=iip2,ip1jm-1
190	dxqu(ij)=q(ij+1,l)-q(ij,l)
191	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
192	c sigu(ij)=u_m(ij,l)/masse(ij,l)
193	ENDDO
194	DO ij=iip1+iip1,ip1jm,iip1
195	dxqu(ij)=dxqu(ij-iim)
196	c sigu(ij)=sigu(ij-iim)
197	ENDDO
198
199	DO ij=iip2,ip1jm
200	adxqu(ij)=abs(dxqu(ij))
201	ENDDO
202
203	c calcul de la pente maximum dans la maille en valeur absolue
204
205	DO ij=iip2+1,ip1jm
206	dxqmax(ij,l)=pente_max*
207	, min(adxqu(ij-1),adxqu(ij))
208	c limitation subtile
209	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
210
211
212	ENDDO
213
214	DO ij=iip1+iip1,ip1jm,iip1
215	dxqmax(ij-iim,l)=dxqmax(ij,l)
216	ENDDO
217
218	DO ij=iip2+1,ip1jm
219	#ifdef CRAY
220	dxq(ij,l)=
221	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
222	#else
223	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
224	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
225	ELSE
226	c extremum local
227	dxq(ij,l)=0.
228	ENDIF
229	#endif
230	dxq(ij,l)=0.5*dxq(ij,l)
231	dxq(ij,l)=
232	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
233	ENDDO
234
235	ENDDO ! l=1,llm
236	cprint*,'Ok calcul des pentes'
237
238	ELSE ! (pente_max.lt.-1.e-5)
239
240	c Pentes produits:
241	c ----------------
242
243	DO l = 1, llm
244	DO ij=iip2,ip1jm-1
245	dxqu(ij)=q(ij+1,l)-q(ij,l)
246	ENDDO
247	DO ij=iip1+iip1,ip1jm,iip1
248	dxqu(ij)=dxqu(ij-iim)
249	ENDDO
250
251	DO ij=iip2+1,ip1jm
252	zz(ij)=dxqu(ij-1)*dxqu(ij)
253	zz(ij)=zz(ij)+zz(ij)
254	IF(zz(ij).gt.0) THEN
255	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
256	ELSE
257	c extremum local
258	dxq(ij,l)=0.
259	ENDIF
260	ENDDO
261
262	ENDDO
263
264	ENDIF ! (pente_max.lt.-1.e-5)
265
266	c bouclage de la pente en iip1:
267	c -----------------------------
268
269	DO l=1,llm
270	DO ij=iip1+iip1,ip1jm,iip1
271	dxq(ij-iim,l)=dxq(ij,l)
272	ENDDO
273	DO ij=1,ip1jmp1
274	iadvplus(ij,l)=0
275	ENDDO
276
277	ENDDO
278
279	c print*,'Bouclage en iip1'
280
281	c calcul des flux a gauche et a droite
282
283	#ifdef CRAY
284
285	DO l=1,llm
286	DO ij=iip2,ip1jm-1
287	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
288	, 1.+u_m(ij,l)/masse(ij+1,l),
289	, u_m(ij,l))
290	zdum(ij,l)=0.5*zdum(ij,l)
291	u_mq(ij,l)=cvmgp(
292	, q(ij,l)+zdum(ij,l)*dxq(ij,l),
293	, q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
294	, u_m(ij,l))
295	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
296	ENDDO
297	ENDDO
298	#else
299	c on cumule le flux correspondant a toutes les mailles dont la masse
300	c au travers de la paroi pENDant le pas de temps.
301	cprint*,'Cumule ....'
302
303	DO l=1,llm
304	DO ij=iip2,ip1jm-1
305	c print*,'masse(',ij,')=',masse(ij,l)
306	IF (u_m(ij,l).gt.0.) THEN
307	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
308	u_mq(ij,l)=u_m(ij,l)(q(ij,l)+0.5zdum(ij,l)*dxq(ij,l))
309	ELSE
310	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
311	u_mq(ij,l)=u_m(ij,l)(q(ij+1,l)-0.5zdum(ij,l)*dxq(ij+1,l))
312	ENDIF
313	ENDDO
314	ENDDO
315	#endif
316	c stop
317
318	c go to 9999
319	c detection des points ou on advecte plus que la masse de la
320	c maille
321	DO l=1,llm
322	DO ij=iip2,ip1jm-1
323	IF(zdum(ij,l).lt.0) THEN
324	iadvplus(ij,l)=1
325	u_mq(ij,l)=0.
326	ENDIF
327	ENDDO
328	ENDDO
329	cprint*,'Ok test 1'
330	DO l=1,llm
331	DO ij=iip1+iip1,ip1jm,iip1
332	iadvplus(ij,l)=iadvplus(ij-iim,l)
333	ENDDO
334	ENDDO
335	c print*,'Ok test 2'
336
337
338	c traitement special pour le cas ou on advecte en longitude plus que le
339	c contenu de la maille.
340	c cette partie est mal vectorisee.
341
342	c calcul du nombre de maille sur lequel on advecte plus que la maille.
343
344	n0=0
345	DO l=1,llm
346	nl(l)=0
347	DO ij=iip2,ip1jm
348	nl(l)=nl(l)+iadvplus(ij,l)
349	ENDDO
350	n0=n0+nl(l)
351	ENDDO
352
353	IF(n0.gt.0) THEN
354	if (prt_level > 2) PRINT *,
355	$ 'Nombre de points pour lesquels on advect plus que le'
356	& ,'contenu de la maille : ',n0
357
358	DO l=1,llm
359	IF(nl(l).gt.0) THEN
360	iju=0
361	c indicage des mailles concernees par le traitement special
362	DO ij=iip2,ip1jm
363	IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
364	iju=iju+1
365	indu(iju)=ij
366	ENDIF
367	ENDDO
368	niju=iju
369	c PRINT*,'niju,nl',niju,nl(l)
370
371	c traitement des mailles
372	DO iju=1,niju
373	ij=indu(iju)
374	j=(ij-1)/iip1+1
375	zu_m=u_m(ij,l)
376	u_mq(ij,l)=0.
377	IF(zu_m.gt.0.) THEN
378	ijq=ij
379	i=ijq-(j-1)*iip1
380	c accumulation pour les mailles completements advectees
381	do while(zu_m.gt.masse(ijq,l))
382	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
383	zu_m=zu_m-masse(ijq,l)
384	i=mod(i-2+iim,iim)+1
385	ijq=(j-1)*iip1+i
386	ENDDO
387	c ajout de la maille non completement advectee
388	u_mq(ij,l)=u_mq(ij,l)+zu_m*
389	& (q(ijq,l)+0.5(1.-zu_m/masse(ijq,l))dxq(ijq,l))
390	ELSE
391	ijq=ij+1
392	i=ijq-(j-1)*iip1
393	c accumulation pour les mailles completements advectees
394	do while(-zu_m.gt.masse(ijq,l))
395	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
396	zu_m=zu_m+masse(ijq,l)
397	i=mod(i,iim)+1
398	ijq=(j-1)*iip1+i
399	ENDDO
400	c ajout de la maille non completement advectee
401	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
402	& 0.5(1.+zu_m/masse(ijq,l))dxq(ijq,l))
403	ENDIF
404	ENDDO
405	ENDIF
406	ENDDO
407	ENDIF ! n0.gt.0
408	9999 continue
409
410
411	c bouclage en latitude
412	cprint*,'cvant bouclage en latitude'
413	DO l=1,llm
414	DO ij=iip1+iip1,ip1jm,iip1
415	u_mq(ij,l)=u_mq(ij-iim,l)
416	ENDDO
417	ENDDO
418
419
420	c calcul des tENDances
421
422	DO l=1,llm
423	DO ij=iip2+1,ip1jm
424	new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
425	q(ij,l)=(q(ij,l)*masse(ij,l)+
426	& u_mq(ij-1,l)-u_mq(ij,l))
427	& /new_m
428	masse(ij,l)=new_m
429	ENDDO
430	c ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
431	DO ij=iip1+iip1,ip1jm,iip1
432	q(ij-iim,l)=q(ij,l)
433	masse(ij-iim,l)=masse(ij,l)
434	ENDDO
435	ENDDO
436	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
437	c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
438
439
440	RETURN
441	END
442	SUBROUTINE vly(q,pente_max,masse,masse_adv_v)
443	c
444	c Auteurs: P.Le Van, F.Hourdin, F.Forget
445	c
446	c ********************************************************************
447	c Shema d'advection " pseudo amont " .
448	c ********************************************************************
449	c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
450	c dq sont des arguments de sortie pour le s-pg ....
451	c
452	c
453	c --------------------------------------------------------------------
454	IMPLICIT NONE
455	c
456	#include "dimensions.h"
457	#include "paramet.h"
458	#include "logic.h"
459	#include "comvert.h"
460	#include "comconst.h"
461	#include "comgeom.h"
462	c
463	c
464	c Arguments:
465	c ----------
466	REAL masse(ip1jmp1,llm),pente_max
467	REAL masse_adv_v( ip1jm,llm)
468	REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
469	c
470	c Local
471	c ---------
472	c
473	INTEGER i,ij,l
474	c
475	REAL airej2,airejjm,airescb(iim),airesch(iim)
476	REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
477	REAL adyqv(ip1jm),dyqmax(ip1jmp1)
478	REAL qbyv(ip1jm,llm)
479
480	REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
481	c REAL newq,oldmasse
482	Logical extremum,first,testcpu
483	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
484	SAVE temps0,temps1,temps2,temps3,temps4,temps5
485	SAVE first,testcpu
486
487	REAL convpn,convps,convmpn,convmps
488	real massepn,masseps,qpn,qps
489	REAL sinlon(iip1),sinlondlon(iip1)
490	REAL coslon(iip1),coslondlon(iip1)
491	SAVE sinlon,coslon,sinlondlon,coslondlon
492	SAVE airej2,airejjm
493	c
494	c
495	REAL SSUM
496
497	DATA first,testcpu/.true.,.false./
498	DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
499
500	IF(first) THEN
501	PRINT*,'Shema Amont nouveau appele dans Vanleer '
502	first=.false.
503	do i=2,iip1
504	coslon(i)=cos(rlonv(i))
505	sinlon(i)=sin(rlonv(i))
506	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
507	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
508	ENDDO
509	coslon(1)=coslon(iip1)
510	coslondlon(1)=coslondlon(iip1)
511	sinlon(1)=sinlon(iip1)
512	sinlondlon(1)=sinlondlon(iip1)
513	airej2 = SSUM( iim, aire(iip2), 1 )
514	airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
515	ENDIF
516
517	c
518	cPRINT*,'CALCUL EN LATITUDE'
519
520	DO l = 1, llm
521	c
522	c --------------------------------
523	c CALCUL EN LATITUDE
524	c --------------------------------
525
526	c On commence par calculer la valeur du traceur moyenne sur le premier cercle
527	c de latitude autour du pole (qpns pour le pole nord et qpsn pour
528	c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
529
530	DO i = 1, iim
531	airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
532	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
533	ENDDO
534	qpns = SSUM( iim, airescb ,1 ) / airej2
535	qpsn = SSUM( iim, airesch ,1 ) / airejjm
536
537	c calcul des pentes aux points v
538
539	DO ij=1,ip1jm
540	dyqv(ij)=q(ij,l)-q(ij+iip1,l)
541	adyqv(ij)=abs(dyqv(ij))
542	ENDDO
543
544	c calcul des pentes aux points scalaires
545
546	DO ij=iip2,ip1jm
547	dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
548	dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
549	dyqmax(ij)=pente_max*dyqmax(ij)
550	ENDDO
551
552	c calcul des pentes aux poles
553
554	DO ij=1,iip1
555	dyq(ij,l)=qpns-q(ij+iip1,l)
556	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
557	ENDDO
558
559	c filtrage de la derivee
560	dyn1=0.
561	dys1=0.
562	dyn2=0.
563	dys2=0.
564	DO ij=1,iim
565	dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
566	dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
567	dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
568	dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
569	ENDDO
570	DO ij=1,iip1
571	dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
572	dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
573	ENDDO
574
575	c calcul des pentes limites aux poles
576
577	goto 8888
578	fn=1.
579	fs=1.
580	DO ij=1,iim
581	IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
582	fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
583	ENDIF
584	IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
585	fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
586	ENDIF
587	ENDDO
588	DO ij=1,iip1
589	dyq(ij,l)=fn*dyq(ij,l)
590	dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
591	ENDDO
592	8888 continue
593	DO ij=1,iip1
594	dyq(ij,l)=0.
595	dyq(ip1jm+ij,l)=0.
596	ENDDO
597
598	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
599	C En memoire de dIFferents tests sur la
600	C limitation des pentes aux poles.
601	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
602	C PRINT*,dyq(1)
603	C PRINT*,dyqv(iip1+1)
604	C appn=abs(dyq(1)/dyqv(iip1+1))
605	C PRINT*,dyq(ip1jm+1)
606	C PRINT*,dyqv(ip1jm-iip1+1)
607	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
608	C DO ij=2,iim
609	C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
610	C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
611	C ENDDO
612	C appn=min(pente_max/appn,1.)
613	C apps=min(pente_max/apps,1.)
614	C
615	C
616	C cas ou on a un extremum au pole
617	C
618	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
619	C & appn=0.
620	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
621	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
622	C & apps=0.
623	C
624	C limitation des pentes aux poles
625	C DO ij=1,iip1
626	C dyq(ij)=appn*dyq(ij)
627	C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
628	C ENDDO
629	C
630	C test
631	C DO ij=1,iip1
632	C dyq(iip1+ij)=0.
633	C dyq(ip1jm+ij-iip1)=0.
634	C ENDDO
635	C DO ij=1,ip1jmp1
636	C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
637	C ENDDO
638	C
639	C changement 10 07 96
640	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
641	C & THEN
642	C DO ij=1,iip1
643	C dyqmax(ij)=0.
644	C ENDDO
645	C ELSE
646	C DO ij=1,iip1
647	C dyqmax(ij)=pente_max*abs(dyqv(ij))
648	C ENDDO
649	C ENDIF
650	C
651	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
652	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
653	C &THEN
654	C DO ij=ip1jm+1,ip1jmp1
655	C dyqmax(ij)=0.
656	C ENDDO
657	C ELSE
658	C DO ij=ip1jm+1,ip1jmp1
659	C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
660	C ENDDO
661	C ENDIF
662	C fin changement 10 07 96
663	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
664
665	c calcul des pentes limitees
666
667	DO ij=iip2,ip1jm
668	IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
669	dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
670	ELSE
671	dyq(ij,l)=0.
672	ENDIF
673	ENDDO
674
675	ENDDO
676
677	DO l=1,llm
678	DO ij=1,ip1jm
679	IF(masse_adv_v(ij,l).gt.0) THEN
680	qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
681	, 0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
682	ELSE
683	qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
684	, 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
685	ENDIF
686	qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
687	ENDDO
688	ENDDO
689
690
691	DO l=1,llm
692	DO ij=iip2,ip1jm
693	newmasse=masse(ij,l)
694	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
695	q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
696	& /newmasse
697	masse(ij,l)=newmasse
698	ENDDO
699	c.-. ancienne version
700	c convpn=SSUM(iim,qbyv(1,l),1)/apoln
701	c convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
702
703	convpn=SSUM(iim,qbyv(1,l),1)
704	convmpn=ssum(iim,masse_adv_v(1,l),1)
705	massepn=ssum(iim,masse(1,l),1)
706	qpn=0.
707	do ij=1,iim
708	qpn=qpn+masse(ij,l)*q(ij,l)
709	enddo
710	qpn=(qpn+convpn)/(massepn+convmpn)
711	do ij=1,iip1
712	q(ij,l)=qpn
713	enddo
714
715	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
716	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
717
718	convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
719	convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
720	masseps=ssum(iim, masse(ip1jm+1,l),1)
721	qps=0.
722	do ij = ip1jm+1,ip1jmp1-1
723	qps=qps+masse(ij,l)*q(ij,l)
724	enddo
725	qps=(qps+convps)/(masseps+convmps)
726	do ij=ip1jm+1,ip1jmp1
727	q(ij,l)=qps
728	enddo
729
730	c.-. fin ancienne version
731
732	c._. nouvelle version
733	c convpn=SSUM(iim,qbyv(1,l),1)
734	c convmpn=ssum(iim,masse_adv_v(1,l),1)
735	c oldmasse=ssum(iim,masse(1,l),1)
736	c newmasse=oldmasse+convmpn
737	c newq=(q(1,l)*oldmasse+convpn)/newmasse
738	c newmasse=newmasse/apoln
739	c DO ij = 1,iip1
740	c q(ij,l)=newq
741	c masse(ij,l)=newmasse*aire(ij)
742	c ENDDO
743	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
744	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
745	c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
746	c newmasse=oldmasse+convmps
747	c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
748	c newmasse=newmasse/apols
749	c DO ij = ip1jm+1,ip1jmp1
750	c q(ij,l)=newq
751	c masse(ij,l)=newmasse*aire(ij)
752	c ENDDO
753	c._. fin nouvelle version
754	ENDDO
755
756	RETURN
757	END
758	SUBROUTINE vlz(q,pente_max,masse,w)
759	c
760	c Auteurs: P.Le Van, F.Hourdin, F.Forget
761	c
762	c ********************************************************************
763	c Shema d'advection " pseudo amont " .
764	c ********************************************************************
765	c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
766	c dq sont des arguments de sortie pour le s-pg ....
767	c
768	c
769	c --------------------------------------------------------------------
770	IMPLICIT NONE
771	c
772	#include "dimensions.h"
773	#include "paramet.h"
774	#include "logic.h"
775	#include "comvert.h"
776	#include "comconst.h"
777	c
778	c
779	c Arguments:
780	c ----------
781	REAL masse(ip1jmp1,llm),pente_max
782	REAL q(ip1jmp1,llm)
783	REAL w(ip1jmp1,llm+1)
784	c
785	c Local
786	c ---------
787	c
788	INTEGER i,ij,l,j,ii
789	c
790	REAL wq(ip1jmp1,llm+1),newmasse
791
792	REAL dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm),dzqmax
793	REAL sigw
794
795	LOGICAL testcpu
796	SAVE testcpu
797
798	REAL temps0,temps1,temps2,temps3,temps4,temps5,second
799	SAVE temps0,temps1,temps2,temps3,temps4,temps5
800	REAL SSUM
801
802	DATA testcpu/.false./
803	DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
804
805	c On oriente tout dans le sens de la pression c'est a dire dans le
806	c sens de W
807
808	#ifdef BIDON
809	IF(testcpu) THEN
810	temps0=second(0.)
811	ENDIF
812	#endif
813	DO l=2,llm
814	DO ij=1,ip1jmp1
815	dzqw(ij,l)=q(ij,l-1)-q(ij,l)
816	adzqw(ij,l)=abs(dzqw(ij,l))
817	ENDDO
818	ENDDO
819
820	DO l=2,llm-1
821	DO ij=1,ip1jmp1
822	#ifdef CRAY
823	dzq(ij,l)=0.5*
824	, cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1))
825	#else
826	IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
827	dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
828	ELSE
829	dzq(ij,l)=0.
830	ENDIF
831	#endif
832	dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
833	dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
834	ENDDO
835	ENDDO
836
837	DO ij=1,ip1jmp1
838	dzq(ij,1)=0.
839	dzq(ij,llm)=0.
840	ENDDO
841
842	#ifdef BIDON
843	IF(testcpu) THEN
844	temps1=temps1+second(0.)-temps0
845	ENDIF
846	#endif
847	c ---------------------------------------------------------------
848	c .... calcul des termes d'advection verticale .......
849	c ---------------------------------------------------------------
850
851	c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq
852
853	DO l = 1,llm-1
854	do ij = 1,ip1jmp1
855	IF(w(ij,l+1).gt.0.) THEN
856	sigw=w(ij,l+1)/masse(ij,l+1)
857	wq(ij,l+1)=w(ij,l+1)(q(ij,l+1)+0.5(1.-sigw)*dzq(ij,l+1))
858	ELSE
859	sigw=w(ij,l+1)/masse(ij,l)
860	wq(ij,l+1)=w(ij,l+1)(q(ij,l)-0.5(1.+sigw)*dzq(ij,l))
861	ENDIF
862	ENDDO
863	ENDDO
864
865	DO ij=1,ip1jmp1
866	wq(ij,llm+1)=0.
867	wq(ij,1)=0.
868	ENDDO
869
870	DO l=1,llm
871	DO ij=1,ip1jmp1
872	newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
873	q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
874	& /newmasse
875	masse(ij,l)=newmasse
876	ENDDO
877	ENDDO
878
879
880	RETURN
881	END
882	c SUBROUTINE minmaxq(zq,qmin,qmax,comment)
883	c
884	c#include "dimensions.h"
885	c#include "paramet.h"
886
887	c CHARACTER() comment
888	c real qmin,qmax
889	c real zq(ip1jmp1,llm)
890
891	c INTEGER jadrs(ip1jmp1), jbad, k, i
892
893
894	c DO k = 1, llm
895	c jbad = 0
896	c DO i = 1, ip1jmp1
897	c IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
898	c jbad = jbad + 1
899	c jadrs(jbad) = i
900	c ENDIF
901	c ENDDO
902	c IF (jbad.GT.0) THEN
903	c PRINT*, comment
904	c DO i = 1, jbad
905	cc PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
906	c ENDDO
907	c ENDIF
908	c ENDDO
909
910	c return
911	c end
912	subroutine minmaxq(zq,qmin,qmax,comment)
913
914	#include "dimensions.h"
915	#include "paramet.h"
916
917	character*20 comment
918	real qmin,qmax
919	real zq(ip1jmp1,llm)
920	real zzq(iip1,jjp1,llm)
921
922	integer imin,jmin,lmin,ijlmin
923	integer imax,jmax,lmax,ijlmax
924
925	integer ismin,ismax
926
927	#ifdef isminismax
928	call scopy (ip1jmp1*llm,zq,1,zzq,1)
929
930	ijlmin=ismin(ijp1llm,zq,1)
931	lmin=(ijlmin-1)/ip1jmp1+1
932	ijlmin=ijlmin-(lmin-1.)*ip1jmp1
933	jmin=(ijlmin-1)/iip1+1
934	imin=ijlmin-(jmin-1.)*iip1
935	zqmin=zq(ijlmin,lmin)
936
937	ijlmax=ismax(ijp1llm,zq,1)
938	lmax=(ijlmax-1)/ip1jmp1+1
939	ijlmax=ijlmax-(lmax-1.)*ip1jmp1
940	jmax=(ijlmax-1)/iip1+1
941	imax=ijlmax-(jmax-1.)*iip1
942	zqmax=zq(ijlmax,lmax)
943
944	if(zqmin.lt.qmin)
945	c s write(*,9999) comment,
946	s write(,) comment,
947	s imin,jmin,lmin,zqmin,zzq(imin,jmin,lmin)
948	if(zqmax.gt.qmax)
949	c s write(*,9999) comment,
950	s write(,) comment,
951	s imax,jmax,lmax,zqmax,zzq(imax,jmax,lmax)
952
953	#endif
954	return
955	9999 format(a20,' q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
956	end
957
958
959

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

Download in other formats:

Original Format