My Project
 All Classes Files Functions Variables Macros
vlspltqs_loc.F
Go to the documentation of this file.
1  SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x)
2 c
3 c Auteurs: P.Le Van, F.Hourdin, F.Forget
4 c
5 c ********************************************************************
6 c Shema d'advection " pseudo amont " .
7 c ********************************************************************
8 c
9 c --------------------------------------------------------------------
10  USE parallel
11  IMPLICIT NONE
12 c
13 #include "dimensions.h"
14 #include "paramet.h"
15 #include "logic.h"
16 #include "comvert.h"
17 #include "comconst.h"
18 c
19 c
20 c Arguments:
21 c ----------
22  REAL masse(ijb_u:ije_u,llm),pente_max
23  REAL u_m( ijb_u:ije_u,llm )
24  REAL q(ijb_u:ije_u,llm)
25  REAL qsat(ijb_u:ije_u,llm)
26 c
27 c Local
28 c ---------
29 c
30  INTEGER ij,l,j,i,iju,ijq,indu(ijnb_u),niju
31  INTEGER n0,iadvplus(ijb_u:ije_u,llm),nl(llm)
32 c
33  REAL new_m,zu_m,zdum(ijb_u:ije_u,llm)
34  REAL dxq(ijb_u:ije_u,llm),dxqu(ijb_u:ije_u)
35  REAL zz(ijb_u:ije_u)
36  REAL adxqu(ijb_u:ije_u),dxqmax(ijb_u:ije_u,llm)
37  REAL u_mq(ijb_u:ije_u,llm)
38 
39  REAL ssum
40 
41 
42  INTEGER ijb,ije,ijb_x,ije_x
43 
44 
45 c calcul de la pente a droite et a gauche de la maille
46 
47 c ijb=ij_begin
48 c ije=ij_end
49 
50  ijb=ijb_x
51  ije=ije_x
52 
53  if (pole_nord.and.ijb==1) ijb=ijb+iip1
54  if (pole_sud.and.ije==ip1jmp1) ije=ije-iip1
55 
56  IF (pente_max.gt.-1.e-5) THEN
57 c IF (pente_max.gt.10) THEN
58 
59 c calcul des pentes avec limitation, Van Leer scheme I:
60 c -----------------------------------------------------
61 
62 c calcul de la pente aux points u
63 
64 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
65  DO l = 1, llm
66  DO ij=ijb,ije-1
67  dxqu(ij)=q(ij+1,l)-q(ij,l)
68 c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
69 c sigu(ij)=u_m(ij,l)/masse(ij,l)
70  ENDDO
71  DO ij=ijb+iip1-1,ije,iip1
72  dxqu(ij)=dxqu(ij-iim)
73 c sigu(ij)=sigu(ij-iim)
74  ENDDO
75 
76  DO ij=ijb,ije
77  adxqu(ij)=abs(dxqu(ij))
78  ENDDO
79 
80 c calcul de la pente maximum dans la maille en valeur absolue
81 
82  DO ij=ijb+1,ije
83  dxqmax(ij,l)=pente_max*
84  , min(adxqu(ij-1),adxqu(ij))
85 c limitation subtile
86 c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
87 
88 
89  ENDDO
90 
91  DO ij=ijb+iip1-1,ije,iip1
92  dxqmax(ij-iim,l)=dxqmax(ij,l)
93  ENDDO
94 
95  DO ij=ijb+1,ije
96 #ifdef CRAY
97  dxq(ij,l)=
98  , cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
99 #else
100  IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
101  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
102  ELSE
103 c extremum local
104  dxq(ij,l)=0.
105  ENDIF
106 #endif
107  dxq(ij,l)=0.5*dxq(ij,l)
108  dxq(ij,l)=
109  , sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
110  ENDDO
111 
112  ENDDO ! l=1,llm
113 c$OMP END DO NOWAIT
114 
115  ELSE ! (pente_max.lt.-1.e-5)
116 
117 c Pentes produits:
118 c ----------------
119 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
120  DO l = 1, llm
121  DO ij=ijb,ije-1
122  dxqu(ij)=q(ij+1,l)-q(ij,l)
123  ENDDO
124  DO ij=ijb+iip1-1,ije,iip1
125  dxqu(ij)=dxqu(ij-iim)
126  ENDDO
127 
128  DO ij=ijb+1,ije
129  zz(ij)=dxqu(ij-1)*dxqu(ij)
130  zz(ij)=zz(ij)+zz(ij)
131  IF(zz(ij).gt.0) THEN
132  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
133  ELSE
134 c extremum local
135  dxq(ij,l)=0.
136  ENDIF
137  ENDDO
138 
139  ENDDO
140 c$OMP END DO NOWAIT
141  ENDIF ! (pente_max.lt.-1.e-5)
142 
143 c bouclage de la pente en iip1:
144 c -----------------------------
145 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
146  DO l=1,llm
147  DO ij=ijb+iip1-1,ije,iip1
148  dxq(ij-iim,l)=dxq(ij,l)
149  ENDDO
150 
151  DO ij=ijb,ije
152  iadvplus(ij,l)=0
153  ENDDO
154 
155  ENDDO
156 c$OMP END DO NOWAIT
157 
158  if (pole_nord) THEN
159 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
160  DO l=1,llm
161  iadvplus(1:iip1,l)=0
162  ENDDO
163 c$OMP END DO NOWAIT
164  endif
165 
166  if (pole_sud) THEN
167 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
168  DO l=1,llm
169  iadvplus(ip1jm+1:ip1jmp1,l)=0
170  ENDDO
171 c$OMP END DO NOWAIT
172  endif
173 
174 c calcul des flux a gauche et a droite
175 
176 #ifdef CRAY
177 c--pas encore modification sur Qsat
178 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
179  DO l=1,llm
180  DO ij=ijb,ije-1
181  zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
182  , 1.+u_m(ij,l)/masse(ij+1,l),
183  , u_m(ij,l))
184  zdum(ij,l)=0.5*zdum(ij,l)
185  u_mq(ij,l)=cvmgp(
186  , q(ij,l)+zdum(ij,l)*dxq(ij,l),
187  , q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
188  , u_m(ij,l))
189  u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
190  ENDDO
191  ENDDO
192 c$OMP END DO NOWAIT
193 
194 #else
195 c on cumule le flux correspondant a toutes les mailles dont la masse
196 c au travers de la paroi pENDant le pas de temps.
197 c le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
198 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
199  DO l=1,llm
200  DO ij=ijb,ije-1
201  IF (u_m(ij,l).gt.0.) THEN
202  zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
203  u_mq(ij,l)=u_m(ij,l)*
204  $ min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
205  ELSE
206  zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
207  u_mq(ij,l)=u_m(ij,l)*
208  $ min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
209  ENDIF
210  ENDDO
211  ENDDO
212 c$OMP END DO NOWAIT
213 #endif
214 
215 
216 c detection des points ou on advecte plus que la masse de la
217 c maille
218 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
219  DO l=1,llm
220  DO ij=ijb,ije-1
221  IF(zdum(ij,l).lt.0) THEN
222  iadvplus(ij,l)=1
223  u_mq(ij,l)=0.
224  ENDIF
225  ENDDO
226  ENDDO
227 c$OMP END DO NOWAIT
228 
229 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
230  DO l=1,llm
231  DO ij=ijb+iip1-1,ije,iip1
232  iadvplus(ij,l)=iadvplus(ij-iim,l)
233  ENDDO
234  ENDDO
235 c$OMP END DO NOWAIT
236 
237 
238 
239 c traitement special pour le cas ou on advecte en longitude plus que le
240 c contenu de la maille.
241 c cette partie est mal vectorisee.
242 
243 c pas d'influence de la pression saturante (pour l'instant)
244 
245 c calcul du nombre de maille sur lequel on advecte plus que la maille.
246 
247  n0=0
248 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
249  DO l=1,llm
250  nl(l)=0
251  DO ij=ijb,ije
252  nl(l)=nl(l)+iadvplus(ij,l)
253  ENDDO
254  n0=n0+nl(l)
255  ENDDO
256 c$OMP END DO NOWAIT
257 
258 cym ATTENTION ICI en OpenMP reduction pas forcement nécessaire
259 cym IF(n0.gt.1) THEN
260 cym IF(n0.gt.0) THEN
261 ccc PRINT*,'Nombre de points pour lesquels on advect plus que le'
262 ccc & ,'contenu de la maille : ',n0
263 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
264  DO l=1,llm
265  IF(nl(l).gt.0) THEN
266  iju=0
267 c indicage des mailles concernees par le traitement special
268  DO ij=ijb,ije
269  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
270  iju=iju+1
271  indu(iju)=ij
272  ENDIF
273  ENDDO
274  niju=iju
275 c PRINT*,'niju,nl',niju,nl(l)
276 
277 c traitement des mailles
278  DO iju=1,niju
279  ij=indu(iju)
280  j=(ij-1)/iip1+1
281  zu_m=u_m(ij,l)
282  u_mq(ij,l)=0.
283  IF(zu_m.gt.0.) THEN
284  ijq=ij
285  i=ijq-(j-1)*iip1
286 c accumulation pour les mailles completements advectees
287  do while(zu_m.gt.masse(ijq,l))
288  u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
289  zu_m=zu_m-masse(ijq,l)
290  i=mod(i-2+iim,iim)+1
291  ijq=(j-1)*iip1+i
292  ENDDO
293 c ajout de la maille non completement advectee
294  u_mq(ij,l)=u_mq(ij,l)+zu_m*
295  & (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
296  ELSE
297  ijq=ij+1
298  i=ijq-(j-1)*iip1
299 c accumulation pour les mailles completements advectees
300  do while(-zu_m.gt.masse(ijq,l))
301  u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
302  zu_m=zu_m+masse(ijq,l)
303  i=mod(i,iim)+1
304  ijq=(j-1)*iip1+i
305  ENDDO
306 c ajout de la maille non completement advectee
307  u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
308  & 0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
309  ENDIF
310  ENDDO
311  ENDIF
312  ENDDO
313 c$OMP END DO NOWAIT
314 cym ENDIF ! n0.gt.0
315 
316 
317 
318 c bouclage en latitude
319 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
320  DO l=1,llm
321  DO ij=ijb+iip1-1,ije,iip1
322  u_mq(ij,l)=u_mq(ij-iim,l)
323  ENDDO
324  ENDDO
325 c$OMP END DO NOWAIT
326 
327 c calcul des tendances
328 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
329  DO l=1,llm
330  DO ij=ijb+1,ije
331  new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
332  q(ij,l)=(q(ij,l)*masse(ij,l)+
333  & u_mq(ij-1,l)-u_mq(ij,l))
334  & /new_m
335  masse(ij,l)=new_m
336  ENDDO
337 c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
338  DO ij=ijb+iip1-1,ije,iip1
339  q(ij-iim,l)=q(ij,l)
340  masse(ij-iim,l)=masse(ij,l)
341  ENDDO
342  ENDDO
343 c$OMP END DO NOWAIT
344 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
345 c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
346 
347 
348  RETURN
349  END
350  SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat)
351 c
352 c Auteurs: P.Le Van, F.Hourdin, F.Forget
353 c
354 c ********************************************************************
355 c Shema d'advection " pseudo amont " .
356 c ********************************************************************
357 c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
358 c qsat est un argument de sortie pour le s-pg ....
359 c
360 c
361 c --------------------------------------------------------------------
362  USE parallel
363  IMPLICIT NONE
364 c
365 #include "dimensions.h"
366 #include "paramet.h"
367 #include "logic.h"
368 #include "comvert.h"
369 #include "comconst.h"
370 #include "comgeom.h"
371 c
372 c
373 c Arguments:
374 c ----------
375  REAL masse(ijb_u:ije_u,llm),pente_max
376  REAL masse_adv_v( ijb_v:ije_v,llm)
377  REAL q(ijb_u:ije_u,llm)
378  REAL qsat(ijb_u:ije_u,llm)
379 c
380 c Local
381 c ---------
382 c
383  INTEGER i,ij,l
384 c
385  REAL airej2,airejjm,airescb(iim),airesch(iim)
386  REAL dyq(ijb_u:ije_u,llm),dyqv(ijb_v:ije_v)
387  REAL adyqv(ijb_v:ije_v),dyqmax(ijb_u:ije_u)
388  REAL qbyv(ijb_v:ije_v,llm)
389 
390  REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
391 c REAL newq,oldmasse
392  Logical first
393  SAVE first
394 c$OMP THREADPRIVATE(first)
395  REAL convpn,convps,convmpn,convmps
396  REAL sinlon(iip1),sinlondlon(iip1)
397  REAL coslon(iip1),coslondlon(iip1)
398  SAVE sinlon,coslon,sinlondlon,coslondlon
399  SAVE airej2,airejjm
400 c$OMP THREADPRIVATE(sinlon,coslon,sinlondlon,coslondlon)
401 c$OMP THREADPRIVATE(airej2,airejjm)
402 c
403 c
404  REAL ssum
405 
406  DATA first/.true./
407  INTEGER ijb,ije
408 
409  IF(first) THEN
410  print*,'Shema Amont nouveau appele dans Vanleer '
411  first=.false.
412  do i=2,iip1
413  coslon(i)=cos(rlonv(i))
414  sinlon(i)=sin(rlonv(i))
415  coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
416  sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
417  ENDDO
418  coslon(1)=coslon(iip1)
419  coslondlon(1)=coslondlon(iip1)
420  sinlon(1)=sinlon(iip1)
421  sinlondlon(1)=sinlondlon(iip1)
422  airej2 = ssum( iim, aire(iip2), 1 )
423  airejjm= ssum( iim, aire(ip1jm -iim), 1 )
424  ENDIF
425 
426 c
427 
428 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
429  DO l = 1, llm
430 c
431 c --------------------------------
432 c CALCUL EN LATITUDE
433 c --------------------------------
434 
435 c On commence par calculer la valeur du traceur moyenne sur le premier cercle
436 c de latitude autour du pole (qpns pour le pole nord et qpsn pour
437 c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
438 
439  if (pole_nord) then
440  DO i = 1, iim
441  airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
442  ENDDO
443  qpns = ssum( iim, airescb ,1 ) / airej2
444  endif
445 
446  if (pole_sud) then
447  DO i = 1, iim
448  airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
449  ENDDO
450  qpsn = ssum( iim, airesch ,1 ) / airejjm
451  endif
452 
453 
454 c calcul des pentes aux points v
455 
456  ijb=ij_begin-2*iip1
457  ije=ij_end+iip1
458  if (pole_nord) ijb=ij_begin
459  if (pole_sud) ije=ij_end-iip1
460 
461  DO ij=ijb,ije
462  dyqv(ij)=q(ij,l)-q(ij+iip1,l)
463  adyqv(ij)=abs(dyqv(ij))
464  ENDDO
465 
466 
467 c calcul des pentes aux points scalaires
468 
469  ijb=ij_begin-iip1
470  ije=ij_end+iip1
471  if (pole_nord) ijb=ij_begin+iip1
472  if (pole_sud) ije=ij_end-iip1
473 
474  DO ij=ijb,ije
475  dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
476  dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
477  dyqmax(ij)=pente_max*dyqmax(ij)
478  ENDDO
479 
480  IF (pole_nord) THEN
481 
482 c calcul des pentes aux poles
483  DO ij=1,iip1
484  dyq(ij,l)=qpns-q(ij+iip1,l)
485  ENDDO
486 
487 c filtrage de la derivee
488  dyn1=0.
489  dyn2=0.
490  DO ij=1,iim
491  dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
492  dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
493  ENDDO
494  DO ij=1,iip1
495  dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
496  ENDDO
497 
498 c calcul des pentes limites aux poles
499  fn=1.
500  DO ij=1,iim
501  IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
502  fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
503  ENDIF
504  ENDDO
505 
506  DO ij=1,iip1
507  dyq(ij,l)=fn*dyq(ij,l)
508  ENDDO
509 
510  ENDIF
511 
512  IF (pole_sud) THEN
513 
514  DO ij=1,iip1
515  dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
516  ENDDO
517 
518  dys1=0.
519  dys2=0.
520 
521  DO ij=1,iim
522  dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
523  dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
524  ENDDO
525 
526  DO ij=1,iip1
527  dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
528  ENDDO
529 
530 c calcul des pentes limites aux poles
531  fs=1.
532  DO ij=1,iim
533  IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
534  fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
535  ENDIF
536  ENDDO
537 
538  DO ij=1,iip1
539  dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
540  ENDDO
541 
542  ENDIF
543 
544 
545 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
546 C En memoire de dIFferents tests sur la
547 C limitation des pentes aux poles.
548 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
549 C PRINT*,dyq(1)
550 C PRINT*,dyqv(iip1+1)
551 C appn=abs(dyq(1)/dyqv(iip1+1))
552 C PRINT*,dyq(ip1jm+1)
553 C PRINT*,dyqv(ip1jm-iip1+1)
554 C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
555 C DO ij=2,iim
556 C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
557 C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
558 C ENDDO
559 C appn=min(pente_max/appn,1.)
560 C apps=min(pente_max/apps,1.)
561 C
562 C
563 C cas ou on a un extremum au pole
564 C
565 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
566 C & appn=0.
567 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
568 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
569 C & apps=0.
570 C
571 C limitation des pentes aux poles
572 C DO ij=1,iip1
573 C dyq(ij)=appn*dyq(ij)
574 C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
575 C ENDDO
576 C
577 C test
578 C DO ij=1,iip1
579 C dyq(iip1+ij)=0.
580 C dyq(ip1jm+ij-iip1)=0.
581 C ENDDO
582 C DO ij=1,ip1jmp1
583 C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
584 C ENDDO
585 C
586 C changement 10 07 96
587 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
588 C & THEN
589 C DO ij=1,iip1
590 C dyqmax(ij)=0.
591 C ENDDO
592 C ELSE
593 C DO ij=1,iip1
594 C dyqmax(ij)=pente_max*abs(dyqv(ij))
595 C ENDDO
596 C ENDIF
597 C
598 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
599 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
600 C &THEN
601 C DO ij=ip1jm+1,ip1jmp1
602 C dyqmax(ij)=0.
603 C ENDDO
604 C ELSE
605 C DO ij=ip1jm+1,ip1jmp1
606 C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
607 C ENDDO
608 C ENDIF
609 C fin changement 10 07 96
610 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
611 
612 c calcul des pentes limitees
613  ijb=ij_begin-iip1
614  ije=ij_end+iip1
615  if (pole_nord) ijb=ij_begin+iip1
616  if (pole_sud) ije=ij_end-iip1
617 
618  DO ij=ijb,ije
619  IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
620  dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
621  ELSE
622  dyq(ij,l)=0.
623  ENDIF
624  ENDDO
625 
626  ENDDO
627 c$OMP END DO NOWAIT
628 
629  ijb=ij_begin-iip1
630  ije=ij_end
631  if (pole_nord) ijb=ij_begin
632  if (pole_sud) ije=ij_end-iip1
633 
634 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
635  DO l=1,llm
636  DO ij=ijb,ije
637  IF( masse_adv_v(ij,l).GT.0. ) THEN
638  qbyv(ij,l)= min( qsat(ij+iip1,l), q(ij+iip1,l ) +
639  , dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
640  ELSE
641  qbyv(ij,l)= min( qsat(ij,l), q(ij,l) - dyq(ij,l) *
642  , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
643  ENDIF
644  qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
645  ENDDO
646  ENDDO
647 c$OMP END DO NOWAIT
648 
649  ijb=ij_begin
650  ije=ij_end
651  if (pole_nord) ijb=ij_begin+iip1
652  if (pole_sud) ije=ij_end-iip1
653 
654 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
655  DO l=1,llm
656  DO ij=ijb,ije
657  newmasse=masse(ij,l)
658  & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
659  q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
660  & /newmasse
661  masse(ij,l)=newmasse
662  ENDDO
663 c.-. ancienne version
664 
665  IF (pole_nord) THEN
666 
667  convpn=ssum(iim,qbyv(1,l),1)/apoln
668  convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
669  DO ij = 1,iip1
670  newmasse=masse(ij,l)+convmpn*aire(ij)
671  q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
672  & newmasse
673  masse(ij,l)=newmasse
674  ENDDO
675 
676  ENDIF
677 
678  IF (pole_sud) THEN
679 
680  convps = -ssum(iim,qbyv(ip1jm-iim,l),1)/apols
681  convmps = -ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
682  DO ij = ip1jm+1,ip1jmp1
683  newmasse=masse(ij,l)+convmps*aire(ij)
684  q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
685  & newmasse
686  masse(ij,l)=newmasse
687  ENDDO
688 
689  ENDIF
690 c.-. fin ancienne version
691 
692 c._. nouvelle version
693 c convpn=SSUM(iim,qbyv(1,l),1)
694 c convmpn=ssum(iim,masse_adv_v(1,l),1)
695 c oldmasse=ssum(iim,masse(1,l),1)
696 c newmasse=oldmasse+convmpn
697 c newq=(q(1,l)*oldmasse+convpn)/newmasse
698 c newmasse=newmasse/apoln
699 c DO ij = 1,iip1
700 c q(ij,l)=newq
701 c masse(ij,l)=newmasse*aire(ij)
702 c ENDDO
703 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
704 c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
705 c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
706 c newmasse=oldmasse+convmps
707 c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
708 c newmasse=newmasse/apols
709 c DO ij = ip1jm+1,ip1jmp1
710 c q(ij,l)=newq
711 c masse(ij,l)=newmasse*aire(ij)
712 c ENDDO
713 c._. fin nouvelle version
714  ENDDO
715 c$OMP END DO NOWAIT
716  RETURN
717  END