Context Navigation

vlspltqs.F @ 4350

Last change on this file since 4350 was 3852, checked in by dcugnet, 3 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: 25.0 KB

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

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

Download in other formats:

Original Format