Context Navigation

vlsplt.F @ 1508

Last change on this file since 1508 was 1508, checked in by emillour, 9 years ago

Common dynamics:
Updates in the dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2325):
IMPORTANT: Modifications for isotopes are only done in dyn3d, not in dyn3dpar

as in LMDZ5 these modifications were done in dyn3dmem.
Related LMDZ5 revisions are r2270 and r2281

in dynlonlat_phylonlat:
add module "grid_atob_m.F90" (a regridding utility so far only used by phylmd/ce0l.F90, used to be dyn3d_common/grid_atob.F)

in misc:
follow up updates on wxios.F (add missing_val module variable)

in dyn3d_common:
pression.F => pression.F90
misc_mod.F90: moved from misc to dyn3d_common
added new iso_verif_dyn.F
covcont.F => covcont.F90
infotrac.F90 : add handling of isotopes (reading of corresponding traceur.def for planets not implemented)
dynetat0.F => dynetat0.F90 with some code factorization
dynredem.F => dynredem.F90 with some code factorization
added dynredem_mod.F90: routines used by dynredem
iniacademic.F90 : added isotopes-related initialization for Earth case

in dyn3d:
added check_isotopes.F
modified (isotopes) advtrac.F90, caladvtrac.F
guide_mod.F90: ported updates
leapfrog.F : (isotopes) updates (NB: call integrd with nqtot tracers)
qminimium.F : adaptations for isotopes (copied over, except that #include comvert.h is not needed).
vlsplt.F: adaptations for isotopes (copied over, except than #include logic.h, comvert.h not needed, and replace "include comconst.h" with use comconst_mod, ONLY: pi)
vlspltqs.F : same as vlsplt.F, but also keeping added modification for CP(T)

in dyn3dpar:
leapfrog_p.F: remove unecessary #ifdef CPP_EARTH cpp flag. and call integrd_p with nqtot tracers (only important for Earth)
dynredem_p.F => dynredem_p.F90 and some code factorization
and no isotopes-relates changes in dyn3dpar (since these changes have been made in LMDZ5 dyn3dmem).

File size: 30.6 KB

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

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

Download in other formats:

Original Format