source: LMDZ6/branches/LMDZ_cdrag_LSCE/libf/dyn3dmem/vlsplt_loc.F @ 5353

Last change on this file since 5353 was 4469, checked in by Laurent Fairhead, 22 months ago

Replaced STOP instructions by calls to abort_gcm for a cleaner exit

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