Context Navigation

vlspltqs.F @ 4402

Last change on this file since 4402 was 4325, checked in by dcugnet, 2 years ago
simplify the parser usage: the getKey_init routine is now embedded in the readTracersFile routine. the initIsotopes routine is now embedded in the readIsotopesFile routine. the database is now unique, but can be changed using the get/setKeysDBase. the derived types descriptions, originally located in trac_types_mod, are moved to readTracFiles_mod. few checkings moved from infotrac to the routine testIsotopes, contained in the readIsotopesFile function from readTracFiles_mod. the readTracersFiles and readIsotopesFile routines no longer use a tracers/isotopes argument. remove tnat and alpha_ideal from infotrac ; use instead getKey to get them where they are used (check_isotopes, dynetat0, iniacademic) the trac_type field %Childs is renamed %Children move the isoSelect routine and the corresponding variables routine from infotrac and infotrac_phy to readTracFiles_mod infotrac_phy routine is now fully independant of the (very similar) routine infotrac (init_infotrac_phy has no arguments left). all the explicit keys of the trac_type are now included in the embedded keys database, accessible using the getKey function. the getKey/addKey routines are expanded to handle vectors of integers, reals, logicals or strings. few subroutines converted into functions with error return value. corrections for isotopic tagging tracers mode (to be continued).
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: 24.6 KB

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

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

Download in other formats:

Original Format