Context Navigation

vlsplt.F @ 4050

Last change on this file since 4050 was 4050, checked in by dcugnet, 2 years ago

Second commit for new tracers.

include most of the keys in the tracers descriptor vector "tracers(:)".
fix in phylmdiso/cv3_routines: fq_* variables were used where their fxt_* counterparts were expected.
multiple IF(nqdesc(iq)>0) and IF(nqfils(iq)>0) tests suppressed, because they are not needed: "do ... enddo" loops with 0 upper bound are not executed.
remove French accents from comments (encoding problem) in phylmdiso/cv3_routines and phylmdiso/cv30_routines.
modifications in "isotopes_verif_mod", where the call to function "iso_verif_tag17_q_deltad_chn" in "iso_verif_tag17_q_deltad_chn" was not detected at linking stage, although defined in the same module (?).

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

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

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

Download in other formats:

Original Format