Context Navigation

vlsplt_p.F @ 2597

Last change on this file since 2597 was 2597, checked in by Ehouarn Millour, 8 years ago
Cleanup in the dynamics: get rid of comconst.h, make it a module comconst_mod. EM
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: 29.8 KB

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

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

Download in other formats:

Original Format