source: LMDZ6/branches/LMDZ-INCA-Dyn/libf/dyn3dmem/vlsplt_loc.F @ 4270

Last change on this file since 4270 was 3800, checked in by Laurent Fairhead, 4 years ago

Modifications nécessaires pour les isotopes
CRisi

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