Context Navigation

vlsplt.F @ 4114

Last change on this file since 4114 was 4064, checked in by dcugnet, 3 years ago

minor fixes (unused variables suppressed, comas after a WRITE() statement, etc.)
parser routines taken from version 7 of https://svn.lmd.jussieu.fr/tracers-parser
few changes in infotrac, and few fixes of (at least) the sequential version:
- uadv and vadv were deallocated twice (fix was lost by mistake just before last commit)
- in [( dum(im), im=1, nm)] implicit loops, ifort evaluates "dum(im)" even if nm==0,

resulting in a crash, "im" being unitialized.

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

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

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

Download in other formats:

Original Format