source: LMDZ6/branches/LMDZ-tracers/libf/dyn3dmem/vlsplt_loc.F @ 4448

Last change on this file since 4448 was 3891, checked in by dcugnet, 3 years ago
  • Bugs corrections:
    • sequential gcm fixed
    • parallel gcm compilation fixed ; to be tested
  • Some generic operations moved from infotrac to readTracFile
  • Fixed algebrical reduction routine, used in the isotopes parameters file.
  • Additional component "comp" in the tracers descriptor derived type "tra",

specifying the model component name(s) (cf. tracers sections) it belongs.

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