Context Navigation

vlspltqs_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: 24.7 KB

Line
1	c
2	c $Id: vlspltqs_p.F 2597 2016-07-22 06:44:47Z emillour $
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	USE comconst_mod, ONLY: cpp
28	IMPLICIT NONE
29
30	c
31	#include "dimensions.h"
32	#include "paramet.h"
33	#include "logic.h"
34	#include "comvert.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	c
241	c
242	c Arguments:
243	c ----------
244	REAL masse(ip1jmp1,llm),pente_max
245	REAL u_m( ip1jmp1,llm )
246	REAL q(ip1jmp1,llm)
247	REAL qsat(ip1jmp1,llm)
248	c
249	c Local
250	c ---------
251	c
252	INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
253	INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
254	c
255	REAL new_m,zu_m,zdum(ip1jmp1,llm)
256	REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
257	REAL zz(ip1jmp1)
258	REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
259	REAL u_mq(ip1jmp1,llm)
260
261	REAL SSUM
262
263
264	INTEGER ijb,ije,ijb_x,ije_x
265
266
267	c calcul de la pente a droite et a gauche de la maille
268
269	c ijb=ij_begin
270	c ije=ij_end
271
272	ijb=ijb_x
273	ije=ije_x
274
275	if (pole_nord.and.ijb==1) ijb=ijb+iip1
276	if (pole_sud.and.ije==ip1jmp1) ije=ije-iip1
277
278	IF (pente_max.gt.-1.e-5) THEN
279	c IF (pente_max.gt.10) THEN
280
281	c calcul des pentes avec limitation, Van Leer scheme I:
282	c -----------------------------------------------------
283
284	c calcul de la pente aux points u
285
286	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
287	DO l = 1, llm
288	DO ij=ijb,ije-1
289	dxqu(ij)=q(ij+1,l)-q(ij,l)
290	c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
291	c sigu(ij)=u_m(ij,l)/masse(ij,l)
292	ENDDO
293	DO ij=ijb+iip1-1,ije,iip1
294	dxqu(ij)=dxqu(ij-iim)
295	c sigu(ij)=sigu(ij-iim)
296	ENDDO
297
298	DO ij=ijb,ije
299	adxqu(ij)=abs(dxqu(ij))
300	ENDDO
301
302	c calcul de la pente maximum dans la maille en valeur absolue
303
304	DO ij=ijb+1,ije
305	dxqmax(ij,l)=pente_max*
306	, min(adxqu(ij-1),adxqu(ij))
307	c limitation subtile
308	c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
309
310
311	ENDDO
312
313	DO ij=ijb+iip1-1,ije,iip1
314	dxqmax(ij-iim,l)=dxqmax(ij,l)
315	ENDDO
316
317	DO ij=ijb+1,ije
318	#ifdef CRAY
319	dxq(ij,l)=
320	, cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
321	#else
322	IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
323	dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
324	ELSE
325	c extremum local
326	dxq(ij,l)=0.
327	ENDIF
328	#endif
329	dxq(ij,l)=0.5*dxq(ij,l)
330	dxq(ij,l)=
331	, sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
332	ENDDO
333
334	ENDDO ! l=1,llm
335	c$OMP END DO NOWAIT
336
337	ELSE ! (pente_max.lt.-1.e-5)
338
339	c Pentes produits:
340	c ----------------
341	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
342	DO l = 1, llm
343	DO ij=ijb,ije-1
344	dxqu(ij)=q(ij+1,l)-q(ij,l)
345	ENDDO
346	DO ij=ijb+iip1-1,ije,iip1
347	dxqu(ij)=dxqu(ij-iim)
348	ENDDO
349
350	DO ij=ijb+1,ije
351	zz(ij)=dxqu(ij-1)*dxqu(ij)
352	zz(ij)=zz(ij)+zz(ij)
353	IF(zz(ij).gt.0) THEN
354	dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
355	ELSE
356	c extremum local
357	dxq(ij,l)=0.
358	ENDIF
359	ENDDO
360
361	ENDDO
362	c$OMP END DO NOWAIT
363	ENDIF ! (pente_max.lt.-1.e-5)
364
365	c bouclage de la pente en iip1:
366	c -----------------------------
367	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
368	DO l=1,llm
369	DO ij=ijb+iip1-1,ije,iip1
370	dxq(ij-iim,l)=dxq(ij,l)
371	ENDDO
372
373	DO ij=ijb,ije
374	iadvplus(ij,l)=0
375	ENDDO
376
377	ENDDO
378	c$OMP END DO NOWAIT
379
380	if (pole_nord) THEN
381	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
382	DO l=1,llm
383	iadvplus(1:iip1,l)=0
384	ENDDO
385	c$OMP END DO NOWAIT
386	endif
387
388	if (pole_sud) THEN
389	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
390	DO l=1,llm
391	iadvplus(ip1jm+1:ip1jmp1,l)=0
392	ENDDO
393	c$OMP END DO NOWAIT
394	endif
395
396	c calcul des flux a gauche et a droite
397
398	#ifdef CRAY
399	c--pas encore modification sur Qsat
400	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
401	DO l=1,llm
402	DO ij=ijb,ije-1
403	zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
404	, 1.+u_m(ij,l)/masse(ij+1,l),
405	, u_m(ij,l))
406	zdum(ij,l)=0.5*zdum(ij,l)
407	u_mq(ij,l)=cvmgp(
408	, q(ij,l)+zdum(ij,l)*dxq(ij,l),
409	, q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
410	, u_m(ij,l))
411	u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
412	ENDDO
413	ENDDO
414	c$OMP END DO NOWAIT
415
416	#else
417	c on cumule le flux correspondant a toutes les mailles dont la masse
418	c au travers de la paroi pENDant le pas de temps.
419	c le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
420	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
421	DO l=1,llm
422	DO ij=ijb,ije-1
423	IF (u_m(ij,l).gt.0.) THEN
424	zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
425	u_mq(ij,l)=u_m(ij,l)*
426	$ min(q(ij,l)+0.5zdum(ij,l)dxq(ij,l),qsat(ij+1,l))
427	ELSE
428	zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
429	u_mq(ij,l)=u_m(ij,l)*
430	$ min(q(ij+1,l)-0.5zdum(ij,l)dxq(ij+1,l),qsat(ij,l))
431	ENDIF
432	ENDDO
433	ENDDO
434	c$OMP END DO NOWAIT
435	#endif
436
437
438	c detection des points ou on advecte plus que la masse de la
439	c maille
440	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
441	DO l=1,llm
442	DO ij=ijb,ije-1
443	IF(zdum(ij,l).lt.0) THEN
444	iadvplus(ij,l)=1
445	u_mq(ij,l)=0.
446	ENDIF
447	ENDDO
448	ENDDO
449	c$OMP END DO NOWAIT
450
451	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
452	DO l=1,llm
453	DO ij=ijb+iip1-1,ije,iip1
454	iadvplus(ij,l)=iadvplus(ij-iim,l)
455	ENDDO
456	ENDDO
457	c$OMP END DO NOWAIT
458
459
460
461	c traitement special pour le cas ou on advecte en longitude plus que le
462	c contenu de la maille.
463	c cette partie est mal vectorisee.
464
465	c pas d'influence de la pression saturante (pour l'instant)
466
467	c calcul du nombre de maille sur lequel on advecte plus que la maille.
468
469	n0=0
470	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
471	DO l=1,llm
472	nl(l)=0
473	DO ij=ijb,ije
474	nl(l)=nl(l)+iadvplus(ij,l)
475	ENDDO
476	n0=n0+nl(l)
477	ENDDO
478	c$OMP END DO NOWAIT
479
480	cym ATTENTION ICI en OpenMP reduction pas forcement nécessaire
481	cym IF(n0.gt.1) THEN
482	cym IF(n0.gt.0) THEN
483	ccc PRINT*,'Nombre de points pour lesquels on advect plus que le'
484	ccc & ,'contenu de la maille : ',n0
485	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
486	DO l=1,llm
487	IF(nl(l).gt.0) THEN
488	iju=0
489	c indicage des mailles concernees par le traitement special
490	DO ij=ijb,ije
491	IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
492	iju=iju+1
493	indu(iju)=ij
494	ENDIF
495	ENDDO
496	niju=iju
497	c PRINT*,'niju,nl',niju,nl(l)
498
499	c traitement des mailles
500	DO iju=1,niju
501	ij=indu(iju)
502	j=(ij-1)/iip1+1
503	zu_m=u_m(ij,l)
504	u_mq(ij,l)=0.
505	IF(zu_m.gt.0.) THEN
506	ijq=ij
507	i=ijq-(j-1)*iip1
508	c accumulation pour les mailles completements advectees
509	do while(zu_m.gt.masse(ijq,l))
510	u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
511	zu_m=zu_m-masse(ijq,l)
512	i=mod(i-2+iim,iim)+1
513	ijq=(j-1)*iip1+i
514	ENDDO
515	c ajout de la maille non completement advectee
516	u_mq(ij,l)=u_mq(ij,l)+zu_m*
517	& (q(ijq,l)+0.5(1.-zu_m/masse(ijq,l))dxq(ijq,l))
518	ELSE
519	ijq=ij+1
520	i=ijq-(j-1)*iip1
521	c accumulation pour les mailles completements advectees
522	do while(-zu_m.gt.masse(ijq,l))
523	u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
524	zu_m=zu_m+masse(ijq,l)
525	i=mod(i,iim)+1
526	ijq=(j-1)*iip1+i
527	ENDDO
528	c ajout de la maille non completement advectee
529	u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
530	& 0.5(1.+zu_m/masse(ijq,l))dxq(ijq,l))
531	ENDIF
532	ENDDO
533	ENDIF
534	ENDDO
535	c$OMP END DO NOWAIT
536	cym ENDIF ! n0.gt.0
537
538
539
540	c bouclage en latitude
541	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
542	DO l=1,llm
543	DO ij=ijb+iip1-1,ije,iip1
544	u_mq(ij,l)=u_mq(ij-iim,l)
545	ENDDO
546	ENDDO
547	c$OMP END DO NOWAIT
548
549	c calcul des tendances
550	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
551	DO l=1,llm
552	DO ij=ijb+1,ije
553	new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
554	q(ij,l)=(q(ij,l)*masse(ij,l)+
555	& u_mq(ij-1,l)-u_mq(ij,l))
556	& /new_m
557	masse(ij,l)=new_m
558	ENDDO
559	c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
560	DO ij=ijb+iip1-1,ije,iip1
561	q(ij-iim,l)=q(ij,l)
562	masse(ij-iim,l)=masse(ij,l)
563	ENDDO
564	ENDDO
565	c$OMP END DO NOWAIT
566	c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
567	c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
568
569
570	RETURN
571	END
572	SUBROUTINE vlyqs_p(q,pente_max,masse,masse_adv_v,qsat)
573	c
574	c Auteurs: P.Le Van, F.Hourdin, F.Forget
575	c
576	c ********************************************************************
577	c Shema d'advection " pseudo amont " .
578	c ********************************************************************
579	c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
580	c qsat est un argument de sortie pour le s-pg ....
581	c
582	c
583	c --------------------------------------------------------------------
584	USE parallel_lmdz
585	USE comconst_mod, ONLY: pi
586	IMPLICIT NONE
587	c
588	#include "dimensions.h"
589	#include "paramet.h"
590	#include "logic.h"
591	#include "comvert.h"
592	#include "comgeom.h"
593	c
594	c
595	c Arguments:
596	c ----------
597	REAL masse(ip1jmp1,llm),pente_max
598	REAL masse_adv_v( ip1jm,llm)
599	REAL q(ip1jmp1,llm)
600	REAL qsat(ip1jmp1,llm)
601	c
602	c Local
603	c ---------
604	c
605	INTEGER i,ij,l
606	c
607	REAL airej2,airejjm,airescb(iim),airesch(iim)
608	REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
609	REAL adyqv(ip1jm),dyqmax(ip1jmp1)
610	REAL qbyv(ip1jm,llm)
611
612	REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
613	c REAL newq,oldmasse
614	Logical first
615	SAVE first
616	c$OMP THREADPRIVATE(first)
617	REAL convpn,convps,convmpn,convmps
618	REAL sinlon(iip1),sinlondlon(iip1)
619	REAL coslon(iip1),coslondlon(iip1)
620	SAVE sinlon,coslon,sinlondlon,coslondlon
621	SAVE airej2,airejjm
622	c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
623	c$OMP THREADPRIVATE(airej2,airejjm)
624	c
625	c
626	REAL SSUM
627
628	DATA first/.true./
629	INTEGER ijb,ije
630
631	IF(first) THEN
632	PRINT*,'Shema Amont nouveau appele dans Vanleer '
633	first=.false.
634	do i=2,iip1
635	coslon(i)=cos(rlonv(i))
636	sinlon(i)=sin(rlonv(i))
637	coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
638	sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
639	ENDDO
640	coslon(1)=coslon(iip1)
641	coslondlon(1)=coslondlon(iip1)
642	sinlon(1)=sinlon(iip1)
643	sinlondlon(1)=sinlondlon(iip1)
644	airej2 = SSUM( iim, aire(iip2), 1 )
645	airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
646	ENDIF
647
648	c
649
650	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
651	DO l = 1, llm
652	c
653	c --------------------------------
654	c CALCUL EN LATITUDE
655	c --------------------------------
656
657	c On commence par calculer la valeur du traceur moyenne sur le premier cercle
658	c de latitude autour du pole (qpns pour le pole nord et qpsn pour
659	c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
660
661	if (pole_nord) then
662	DO i = 1, iim
663	airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
664	ENDDO
665	qpns = SSUM( iim, airescb ,1 ) / airej2
666	endif
667
668	if (pole_sud) then
669	DO i = 1, iim
670	airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
671	ENDDO
672	qpsn = SSUM( iim, airesch ,1 ) / airejjm
673	endif
674
675
676	c calcul des pentes aux points v
677
678	ijb=ij_begin-2*iip1
679	ije=ij_end+iip1
680	if (pole_nord) ijb=ij_begin
681	if (pole_sud) ije=ij_end-iip1
682
683	DO ij=ijb,ije
684	dyqv(ij)=q(ij,l)-q(ij+iip1,l)
685	adyqv(ij)=abs(dyqv(ij))
686	ENDDO
687
688
689	c calcul des pentes aux points scalaires
690
691	ijb=ij_begin-iip1
692	ije=ij_end+iip1
693	if (pole_nord) ijb=ij_begin+iip1
694	if (pole_sud) ije=ij_end-iip1
695
696	DO ij=ijb,ije
697	dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
698	dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
699	dyqmax(ij)=pente_max*dyqmax(ij)
700	ENDDO
701
702	IF (pole_nord) THEN
703
704	c calcul des pentes aux poles
705	DO ij=1,iip1
706	dyq(ij,l)=qpns-q(ij+iip1,l)
707	ENDDO
708
709	c filtrage de la derivee
710	dyn1=0.
711	dyn2=0.
712	DO ij=1,iim
713	dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
714	dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
715	ENDDO
716	DO ij=1,iip1
717	dyq(ij,l)=dyn1sinlon(ij)+dyn2coslon(ij)
718	ENDDO
719
720	c calcul des pentes limites aux poles
721	fn=1.
722	DO ij=1,iim
723	IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
724	fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
725	ENDIF
726	ENDDO
727
728	DO ij=1,iip1
729	dyq(ij,l)=fn*dyq(ij,l)
730	ENDDO
731
732	ENDIF
733
734	IF (pole_sud) THEN
735
736	DO ij=1,iip1
737	dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
738	ENDDO
739
740	dys1=0.
741	dys2=0.
742
743	DO ij=1,iim
744	dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
745	dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
746	ENDDO
747
748	DO ij=1,iip1
749	dyq(ip1jm+ij,l)=dys1sinlon(ij)+dys2coslon(ij)
750	ENDDO
751
752	c calcul des pentes limites aux poles
753	fs=1.
754	DO ij=1,iim
755	IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
756	fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
757	ENDIF
758	ENDDO
759
760	DO ij=1,iip1
761	dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
762	ENDDO
763
764	ENDIF
765
766
767	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
768	C En memoire de dIFferents tests sur la
769	C limitation des pentes aux poles.
770	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
771	C PRINT*,dyq(1)
772	C PRINT*,dyqv(iip1+1)
773	C appn=abs(dyq(1)/dyqv(iip1+1))
774	C PRINT*,dyq(ip1jm+1)
775	C PRINT*,dyqv(ip1jm-iip1+1)
776	C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
777	C DO ij=2,iim
778	C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
779	C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
780	C ENDDO
781	C appn=min(pente_max/appn,1.)
782	C apps=min(pente_max/apps,1.)
783	C
784	C
785	C cas ou on a un extremum au pole
786	C
787	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
788	C & appn=0.
789	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
790	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
791	C & apps=0.
792	C
793	C limitation des pentes aux poles
794	C DO ij=1,iip1
795	C dyq(ij)=appn*dyq(ij)
796	C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
797	C ENDDO
798	C
799	C test
800	C DO ij=1,iip1
801	C dyq(iip1+ij)=0.
802	C dyq(ip1jm+ij-iip1)=0.
803	C ENDDO
804	C DO ij=1,ip1jmp1
805	C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
806	C ENDDO
807	C
808	C changement 10 07 96
809	C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
810	C & THEN
811	C DO ij=1,iip1
812	C dyqmax(ij)=0.
813	C ENDDO
814	C ELSE
815	C DO ij=1,iip1
816	C dyqmax(ij)=pente_max*abs(dyqv(ij))
817	C ENDDO
818	C ENDIF
819	C
820	C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
821	C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
822	C &THEN
823	C DO ij=ip1jm+1,ip1jmp1
824	C dyqmax(ij)=0.
825	C ENDDO
826	C ELSE
827	C DO ij=ip1jm+1,ip1jmp1
828	C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
829	C ENDDO
830	C ENDIF
831	C fin changement 10 07 96
832	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
833
834	c calcul des pentes limitees
835	ijb=ij_begin-iip1
836	ije=ij_end+iip1
837	if (pole_nord) ijb=ij_begin+iip1
838	if (pole_sud) ije=ij_end-iip1
839
840	DO ij=ijb,ije
841	IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
842	dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
843	ELSE
844	dyq(ij,l)=0.
845	ENDIF
846	ENDDO
847
848	ENDDO
849	c$OMP END DO NOWAIT
850
851	ijb=ij_begin-iip1
852	ije=ij_end
853	if (pole_nord) ijb=ij_begin
854	if (pole_sud) ije=ij_end-iip1
855
856	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
857	DO l=1,llm
858	DO ij=ijb,ije
859	IF( masse_adv_v(ij,l).GT.0. ) THEN
860	qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l ) +
861	, dyq(ij+iip1,l)0.5(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
862	ELSE
863	qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
864	, 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
865	ENDIF
866	qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
867	ENDDO
868	ENDDO
869	c$OMP END DO NOWAIT
870
871	ijb=ij_begin
872	ije=ij_end
873	if (pole_nord) ijb=ij_begin+iip1
874	if (pole_sud) ije=ij_end-iip1
875
876	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
877	DO l=1,llm
878	DO ij=ijb,ije
879	newmasse=masse(ij,l)
880	& +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
881	q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
882	& /newmasse
883	masse(ij,l)=newmasse
884	ENDDO
885	c.-. ancienne version
886
887	IF (pole_nord) THEN
888
889	convpn=SSUM(iim,qbyv(1,l),1)/apoln
890	convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
891	DO ij = 1,iip1
892	newmasse=masse(ij,l)+convmpn*aire(ij)
893	q(ij,l)=(q(ij,l)masse(ij,l)+convpnaire(ij))/
894	& newmasse
895	masse(ij,l)=newmasse
896	ENDDO
897
898	ENDIF
899
900	IF (pole_sud) THEN
901
902	convps = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
903	convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
904	DO ij = ip1jm+1,ip1jmp1
905	newmasse=masse(ij,l)+convmps*aire(ij)
906	q(ij,l)=(q(ij,l)masse(ij,l)+convpsaire(ij))/
907	& newmasse
908	masse(ij,l)=newmasse
909	ENDDO
910
911	ENDIF
912	c.-. fin ancienne version
913
914	c._. nouvelle version
915	c convpn=SSUM(iim,qbyv(1,l),1)
916	c convmpn=ssum(iim,masse_adv_v(1,l),1)
917	c oldmasse=ssum(iim,masse(1,l),1)
918	c newmasse=oldmasse+convmpn
919	c newq=(q(1,l)*oldmasse+convpn)/newmasse
920	c newmasse=newmasse/apoln
921	c DO ij = 1,iip1
922	c q(ij,l)=newq
923	c masse(ij,l)=newmasse*aire(ij)
924	c ENDDO
925	c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
926	c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
927	c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
928	c newmasse=oldmasse+convmps
929	c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
930	c newmasse=newmasse/apols
931	c DO ij = ip1jm+1,ip1jmp1
932	c q(ij,l)=newq
933	c masse(ij,l)=newmasse*aire(ij)
934	c ENDDO
935	c._. fin nouvelle version
936	ENDDO
937	c$OMP END DO NOWAIT
938	RETURN
939	END

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

Download in other formats:

Original Format