Context Navigation

vlsplt.F @ 1549

Last change on this file since 1549 was 1520, checked in by Ehouarn Millour, 13 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 25.1 KB

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

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

Download in other formats:

Original Format