Context Navigation

vlspltqs_p.F @ 2146

Last change on this file since 2146 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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

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

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

Download in other formats:

Original Format