Context Navigation

vlsplt.F @ 3863

Last change on this file since 3863 was 3852, checked in by dcugnet, 4 years ago

Extension of the tracers management.

The tracers files can be:

1) "traceur.def": old format, with:

the number of tracers on the first line
one line for each tracer: <tracer name> <hadv> <vadv> [<parent name>]

2) "tracer.def": new format with one section each model component.
3) "tracer_<name>.def": new format with a single section.

The formats 2 and 3 reading is driven by the "type_trac" key, which can be a

coma-separated list of components.

Format 2: read the sections from the "tracer.def" file.
format 3: read one section each "tracer_<section name>.def" file.
the first line of a section is "&<section name>
the other lines start with a tracer name followed by <key>=<val> pairs.
the "default" tracer name is reserved ; the other tracers of the section inherit its <key>=<val>, except for the keys that are redefined locally.

This format helps keeping the tracers files compact, thanks to the "default"
special tracer and the three levels of factorization:

on the tracers names: a tracer name can be a coma-separated list of tracers => all the tracers of the list have the same <key>=<val> properties
on the parents names: the value of the "parent" property can be a coma-separated list of tracers => only possible for geographic tagging tracers
on the phases: the property "phases" is [g](l][s] (gas/liquid/solid)

Read information is stored in the vector "tracers(:)", of derived type "tra".

"isotopes_params.def" is a similar file, with one section each isotopes family.
It contains a database of isotopes properties ; if there are second generation
tracers (isotopes), the corresponding sections are read.

Read information is stored in the vector "isotopes(:)", of derived type "iso".

The "getKey" function helps to get the values of the parameters stored in
"tracers" or "isotopes".

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.7 KB

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