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