Context Navigation

vlsplt.F @ 2720

Last change on this file since 2720 was 2286, checked in by Ehouarn Millour, 10 years ago

fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
comment out many invasive debug writes.

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: 30.8 KB

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

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

Download in other formats:

Original Format