source: LMDZ5/trunk/libf/dyn3dmem/vlsplt_loc.F @ 4095

Last change on this file since 4095 was 2765, checked in by fhourdin, 8 years ago

Inclusion d'un traitement spécifique dans l'advection verticale
pour les cas de nombre de courant w/masse > 1
ou w est le transfert de masse entre deux couches.

ce traitement garanti la convergence numérique quand le critère
n'est jamais violé.
Dans les autres cas, il fait proprement le transport vertical
qui est autrement hors de contrôle.
Pourrait contribuer a resoudre les plantages himalayens.

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