source: LMDZ5/trunk/libf/dyn3d/vlspltqs.F @ 2442

Last change on this file since 2442 was 2286, checked in by Ehouarn Millour, 9 years ago
  • fix in dyn3d the array out of bound issue (in qminimum) that was only corrected in dyn3dmem in rev 2285.
  • comment out many invasive debug writes.

EM

  • 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:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.1 KB
Line 
1c
2c $Id: vlspltqs.F 2286 2015-05-20 13:27:07Z fhourdin $
3c
4       SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
5     ,                                  p,pk,teta,iq             )
6       USE infotrac, ONLY: nqtot,nqdesc,iqfils
7c
8c     Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron
9c
10c    ********************************************************************
11c          Shema  d'advection " pseudo amont " .
12c      + test sur humidite specifique: Q advecte< Qsat aval
13c                   (F. Codron, 10/99)
14c    ********************************************************************
15c     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
16c
17c     pente_max facteur de limitation des pentes: 2 en general
18c                                                0 pour un schema amont
19c     pbaru,pbarv,w flux de masse en u ,v ,w
20c     pdt pas de temps
21c
22c     teta temperature potentielle, p pression aux interfaces,
23c     pk exner au milieu des couches necessaire pour calculer Qsat
24c   --------------------------------------------------------------------
25      IMPLICIT NONE
26c
27#include "dimensions.h"
28#include "paramet.h"
29#include "logic.h"
30#include "comvert.h"
31#include "comconst.h"
32
33c
34c   Arguments:
35c   ----------
36      REAL masse(ip1jmp1,llm),pente_max
37      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
38      REAL q(ip1jmp1,llm,nqtot)
39      REAL w(ip1jmp1,llm),pdt
40      REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
41      INTEGER iq ! CRisi
42c
43c      Local
44c   ---------
45c
46      INTEGER i,ij,l,j,ii
47      INTEGER ifils,iq2 ! CRisi
48c
49      REAL qsat(ip1jmp1,llm)
50      REAL zm(ip1jmp1,llm,nqtot)
51      REAL mu(ip1jmp1,llm)
52      REAL mv(ip1jm,llm)
53      REAL mw(ip1jmp1,llm+1)
54      REAL zq(ip1jmp1,llm,nqtot)
55      REAL temps1,temps2,temps3
56      REAL zzpbar, zzw
57      LOGICAL testcpu
58      SAVE testcpu
59      SAVE temps1,temps2,temps3
60
61      REAL qmin,qmax
62      DATA qmin,qmax/0.,1.e33/
63      DATA testcpu/.false./
64      DATA temps1,temps2,temps3/0.,0.,0./
65
66c--pour rapport de melange saturant--
67
68      REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
69      REAL ptarg,pdelarg,foeew,zdelta
70      REAL tempe(ip1jmp1)
71
72c    fonction psat(T)
73
74       FOEEW ( PTARG,PDELARG ) = EXP (
75     *          (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT)
76     * / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
77
78        r2es  = 380.11733
79        r3les = 17.269
80        r3ies = 21.875
81        r4les = 35.86
82        r4ies = 7.66
83        retv = 0.6077667
84        rtt  = 273.16
85
86c-- Calcul de Qsat en chaque point
87c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
88c   pour eviter une exponentielle.
89        DO l = 1, llm
90         DO ij = 1, ip1jmp1
91          tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
92         ENDDO
93         DO ij = 1, ip1jmp1
94          zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
95          play   = 0.5*(p(ij,l)+p(ij,l+1))
96          qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
97          qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
98         ENDDO
99        ENDDO
100
101c      PRINT*,'Debut vlsplt version debug sans vlyqs'
102
103        zzpbar = 0.5 * pdt
104        zzw    = pdt
105      DO l=1,llm
106        DO ij = iip2,ip1jm
107            mu(ij,l)=pbaru(ij,l) * zzpbar
108         ENDDO
109         DO ij=1,ip1jm
110            mv(ij,l)=pbarv(ij,l) * zzpbar
111         ENDDO
112         DO ij=1,ip1jmp1
113            mw(ij,l)=w(ij,l) * zzw
114         ENDDO
115      ENDDO
116
117      DO ij=1,ip1jmp1
118         mw(ij,llm+1)=0.
119      ENDDO
120
121      CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1)
122      CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
123      if (nqdesc(iq).gt.0) then 
124       do ifils=1,nqdesc(iq)
125        iq2=iqfils(ifils,iq)
126        CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
127       enddo 
128      endif !if (nqfils(iq).gt.0) then
129
130c      call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
131      call vlxqs(zq,pente_max,zm,mu,qsat,iq)
132
133c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
134
135      call vlyqs(zq,pente_max,zm,mv,qsat,iq)
136
137c      call minmaxq(zq,qmin,qmax,'avant vlz     ')
138
139      call vlz(zq,pente_max,zm,mw,iq)
140
141c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
142c     call minmaxq(zm,qmin,qmax,'M avant vlyqs     ')
143
144      call vlyqs(zq,pente_max,zm,mv,qsat,iq)
145
146c     call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
147c     call minmaxq(zm,qmin,qmax,'M avant vlxqs     ')
148
149      call vlxqs(zq,pente_max,zm,mu,qsat,iq)
150
151c     call minmaxq(zq,qmin,qmax,'apres vlxqs     ')
152c     call minmaxq(zm,qmin,qmax,'M apres vlxqs     ')
153
154
155      DO l=1,llm
156         DO ij=1,ip1jmp1
157           q(ij,l,iq)=zq(ij,l,iq)
158         ENDDO
159         DO ij=1,ip1jm+1,iip1
160            q(ij+iim,l,iq)=q(ij,l,iq)
161         ENDDO
162      ENDDO
163      ! CRisi: aussi pour les fils
164      if (nqdesc(iq).gt.0) then
165      do ifils=1,nqdesc(iq)
166        iq2=iqfils(ifils,iq)
167        DO l=1,llm
168         DO ij=1,ip1jmp1
169           q(ij,l,iq2)=zq(ij,l,iq2)
170         ENDDO
171         DO ij=1,ip1jm+1,iip1
172            q(ij+iim,l,iq2)=q(ij,l,iq2)
173         ENDDO
174        ENDDO
175      enddo !do ifils=1,nqdesc(iq) 
176      endif ! if (nqfils(iq).gt.0) then
177      !write(*,*) 'vlspltqs 183: fin de la routine'
178
179      RETURN
180      END
181      SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat,iq)
182      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
183
184c
185c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
186c
187c    ********************************************************************
188c     Shema  d'advection " pseudo amont " .
189c    ********************************************************************
190c
191c   --------------------------------------------------------------------
192      IMPLICIT NONE
193c
194#include "dimensions.h"
195#include "paramet.h"
196#include "logic.h"
197#include "comvert.h"
198#include "comconst.h"
199c
200c
201c   Arguments:
202c   ----------
203      REAL masse(ip1jmp1,llm,nqtot),pente_max
204      REAL u_m( ip1jmp1,llm )
205      REAL q(ip1jmp1,llm,nqtot)
206      REAL qsat(ip1jmp1,llm)
207      INTEGER iq ! CRisi
208c
209c      Local
210c   ---------
211c
212      INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
213      INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
214c
215      REAL new_m,zu_m,zdum(ip1jmp1,llm)
216      REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
217      REAL zz(ip1jmp1)
218      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
219      REAL u_mq(ip1jmp1,llm)
220
221      ! CRisi
222      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
223      INTEGER ifils,iq2 ! CRisi
224
225      Logical first,testcpu
226      SAVE first,testcpu
227
228      REAL      SSUM
229      REAL temps0,temps1,temps2,temps3,temps4,temps5
230      SAVE temps0,temps1,temps2,temps3,temps4,temps5
231
232
233      DATA first,testcpu/.true.,.false./
234
235      IF(first) THEN
236         temps1=0.
237         temps2=0.
238         temps3=0.
239         temps4=0.
240         temps5=0.
241         first=.false.
242      ENDIF
243
244c   calcul de la pente a droite et a gauche de la maille
245
246
247      IF (pente_max.gt.-1.e-5) THEN
248c     IF (pente_max.gt.10) THEN
249
250c   calcul des pentes avec limitation, Van Leer scheme I:
251c   -----------------------------------------------------
252
253c   calcul de la pente aux points u
254         DO l = 1, llm
255            DO ij=iip2,ip1jm-1
256               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
257c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
258c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
259            ENDDO
260            DO ij=iip1+iip1,ip1jm,iip1
261               dxqu(ij)=dxqu(ij-iim)
262c              sigu(ij)=sigu(ij-iim)
263            ENDDO
264
265            DO ij=iip2,ip1jm
266               adxqu(ij)=abs(dxqu(ij))
267            ENDDO
268
269c   calcul de la pente maximum dans la maille en valeur absolue
270
271            DO ij=iip2+1,ip1jm
272               dxqmax(ij,l)=pente_max*
273     ,      min(adxqu(ij-1),adxqu(ij))
274c limitation subtile
275c    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
276         
277
278            ENDDO
279
280            DO ij=iip1+iip1,ip1jm,iip1
281               dxqmax(ij-iim,l)=dxqmax(ij,l)
282            ENDDO
283
284            DO ij=iip2+1,ip1jm
285#ifdef CRAY
286               dxq(ij,l)=
287     ,         cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
288#else
289               IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
290                  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
291               ELSE
292c   extremum local
293                  dxq(ij,l)=0.
294               ENDIF
295#endif
296               dxq(ij,l)=0.5*dxq(ij,l)
297               dxq(ij,l)=
298     ,         sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
299            ENDDO
300
301         ENDDO ! l=1,llm
302
303      ELSE ! (pente_max.lt.-1.e-5)
304
305c   Pentes produits:
306c   ----------------
307
308         DO l = 1, llm
309            DO ij=iip2,ip1jm-1
310               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
311            ENDDO
312            DO ij=iip1+iip1,ip1jm,iip1
313               dxqu(ij)=dxqu(ij-iim)
314            ENDDO
315
316            DO ij=iip2+1,ip1jm
317               zz(ij)=dxqu(ij-1)*dxqu(ij)
318               zz(ij)=zz(ij)+zz(ij)
319               IF(zz(ij).gt.0) THEN
320                  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
321               ELSE
322c   extremum local
323                  dxq(ij,l)=0.
324               ENDIF
325            ENDDO
326
327         ENDDO
328
329      ENDIF ! (pente_max.lt.-1.e-5)
330
331c   bouclage de la pente en iip1:
332c   -----------------------------
333
334      DO l=1,llm
335         DO ij=iip1+iip1,ip1jm,iip1
336            dxq(ij-iim,l)=dxq(ij,l)
337         ENDDO
338
339         DO ij=1,ip1jmp1
340            iadvplus(ij,l)=0
341         ENDDO
342
343      ENDDO
344
345
346c   calcul des flux a gauche et a droite
347
348#ifdef CRAY
349c--pas encore modification sur Qsat
350      DO l=1,llm
351       DO ij=iip2,ip1jm-1
352          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
353     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
354     ,                     u_m(ij,l))
355          zdum(ij,l)=0.5*zdum(ij,l)
356          u_mq(ij,l)=cvmgp(
357     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
358     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
359     ,                u_m(ij,l))
360          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
361       ENDDO
362      ENDDO
363#else
364c   on cumule le flux correspondant a toutes les mailles dont la masse
365c   au travers de la paroi pENDant le pas de temps.
366c   le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
367      DO l=1,llm
368       DO ij=iip2,ip1jm-1
369          IF (u_m(ij,l).gt.0.) THEN
370             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
371             u_mq(ij,l)=u_m(ij,l)*
372     $         min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
373          ELSE
374             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
375             u_mq(ij,l)=u_m(ij,l)*
376     $         min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
377          ENDIF
378       ENDDO
379      ENDDO
380#endif
381
382
383c   detection des points ou on advecte plus que la masse de la
384c   maille
385      DO l=1,llm
386         DO ij=iip2,ip1jm-1
387            IF(zdum(ij,l).lt.0) THEN
388               iadvplus(ij,l)=1
389               u_mq(ij,l)=0.
390            ENDIF
391         ENDDO
392      ENDDO
393      DO l=1,llm
394       DO ij=iip1+iip1,ip1jm,iip1
395          iadvplus(ij,l)=iadvplus(ij-iim,l)
396       ENDDO
397      ENDDO
398
399
400
401c   traitement special pour le cas ou on advecte en longitude plus que le
402c   contenu de la maille.
403c   cette partie est mal vectorisee.
404
405c   pas d'influence de la pression saturante (pour l'instant)
406
407c  calcul du nombre de maille sur lequel on advecte plus que la maille.
408
409      n0=0
410      DO l=1,llm
411         nl(l)=0
412         DO ij=iip2,ip1jm
413            nl(l)=nl(l)+iadvplus(ij,l)
414         ENDDO
415         n0=n0+nl(l)
416      ENDDO
417
418      IF(n0.gt.0) THEN
419ccc      PRINT*,'Nombre de points pour lesquels on advect plus que le'
420ccc     &       ,'contenu de la maille : ',n0
421
422         DO l=1,llm
423            IF(nl(l).gt.0) THEN
424               iju=0
425c   indicage des mailles concernees par le traitement special
426               DO ij=iip2,ip1jm
427                  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
428                     iju=iju+1
429                     indu(iju)=ij
430                  ENDIF
431               ENDDO
432               niju=iju
433c              PRINT*,'niju,nl',niju,nl(l)
434
435c  traitement des mailles
436               DO iju=1,niju
437                  ij=indu(iju)
438                  j=(ij-1)/iip1+1
439                  zu_m=u_m(ij,l)
440                  u_mq(ij,l)=0.
441                  IF(zu_m.gt.0.) THEN
442                     ijq=ij
443                     i=ijq-(j-1)*iip1
444c   accumulation pour les mailles completements advectees
445                     do while(zu_m.gt.masse(ijq,l,iq))
446                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
447     &                          *masse(ijq,l,iq)
448                        zu_m=zu_m-masse(ijq,l,iq)
449                        i=mod(i-2+iim,iim)+1
450                        ijq=(j-1)*iip1+i
451                     ENDDO
452c   ajout de la maille non completement advectee
453                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
454     &                  (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
455     &                  *dxq(ijq,l))
456                  ELSE
457                     ijq=ij+1
458                     i=ijq-(j-1)*iip1
459c   accumulation pour les mailles completements advectees
460                     do while(-zu_m.gt.masse(ijq,l,iq))
461                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
462     &                          *masse(ijq,l,iq)
463                        zu_m=zu_m+masse(ijq,l,iq)
464                        i=mod(i,iim)+1
465                        ijq=(j-1)*iip1+i
466                     ENDDO
467c   ajout de la maille non completement advectee
468                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
469     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
470                  ENDIF
471               ENDDO
472            ENDIF
473         ENDDO
474      ENDIF  ! n0.gt.0
475
476
477
478c   bouclage en latitude
479
480      DO l=1,llm
481        DO ij=iip1+iip1,ip1jm,iip1
482           u_mq(ij,l)=u_mq(ij-iim,l)
483        ENDDO
484      ENDDO
485
486! CRisi: appel récursif de l'advection sur les fils.
487! Il faut faire ça avant d'avoir mis à jour q et masse
488      !write(*,*) 'vlspltqs 326: iq,nqfils(iq)=',iq,nqfils(iq)
489     
490      if (nqfils(iq).gt.0) then 
491       do ifils=1,nqdesc(iq)
492         iq2=iqfils(ifils,iq)
493         DO l=1,llm
494          DO ij=iip2,ip1jm
495           ! On a besoin de q et masse seulement entre iip2 et ip1jm
496           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
497           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
498          enddo   
499         enddo
500        enddo !do ifils=1,nqdesc(iq)
501        do ifils=1,nqfils(iq)
502         iq2=iqfils(ifils,iq)
503         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
504        enddo !do ifils=1,nqfils(iq)
505      endif !if (nqfils(iq).gt.0) then
506! end CRisi
507
508c   calcul des tendances
509
510      DO l=1,llm
511         DO ij=iip2+1,ip1jm
512            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
513            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
514     &      u_mq(ij-1,l)-u_mq(ij,l))
515     &      /new_m
516            masse(ij,l,iq)=new_m
517         ENDDO
518c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
519         DO ij=iip1+iip1,ip1jm,iip1
520            q(ij-iim,l,iq)=q(ij,l,iq)
521            masse(ij-iim,l,iq)=masse(ij,l,iq)
522         ENDDO
523      ENDDO
524
525      ! retablir les fils en rapport de melange par rapport a l'air:
526      ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
527      ! puis on boucle en longitude
528      if (nqdesc(iq).gt.0) then 
529       do ifils=1,nqdesc(iq)
530         iq2=iqfils(ifils,iq) 
531         DO l=1,llm
532          DO ij=iip2+1,ip1jm
533            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
534          enddo
535          DO ij=iip1+iip1,ip1jm,iip1
536             q(ij-iim,l,iq2)=q(ij,l,iq2)
537          enddo ! DO ij=ijb+iip1-1,ije,iip1
538         enddo !DO l=1,llm
539        enddo !do ifils=1,nqdesc(iq)
540      endif !if (nqfils(iq).gt.0) then
541
542c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
543c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
544
545
546      RETURN
547      END
548      SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq)
549      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
550c
551c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
552c
553c    ********************************************************************
554c     Shema  d'advection " pseudo amont " .
555c    ********************************************************************
556c     q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
557c     qsat             est   un argument de sortie pour le s-pg ....
558c
559c
560c   --------------------------------------------------------------------
561      IMPLICIT NONE
562c
563#include "dimensions.h"
564#include "paramet.h"
565#include "logic.h"
566#include "comvert.h"
567#include "comconst.h"
568#include "comgeom.h"
569c
570c
571c   Arguments:
572c   ----------
573      REAL masse(ip1jmp1,llm,nqtot),pente_max
574      REAL masse_adv_v( ip1jm,llm)
575      REAL q(ip1jmp1,llm,nqtot)
576      REAL qsat(ip1jmp1,llm)
577      INTEGER iq ! CRisi
578c
579c      Local
580c   ---------
581c
582      INTEGER i,ij,l
583c
584      REAL airej2,airejjm,airescb(iim),airesch(iim)
585      REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
586      REAL adyqv(ip1jm),dyqmax(ip1jmp1)
587      REAL qbyv(ip1jm,llm)
588
589      REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
590c     REAL newq,oldmasse
591      Logical first,testcpu
592      REAL temps0,temps1,temps2,temps3,temps4,temps5
593      SAVE temps0,temps1,temps2,temps3,temps4,temps5
594      SAVE first,testcpu
595
596      REAL convpn,convps,convmpn,convmps
597      REAL sinlon(iip1),sinlondlon(iip1)
598      REAL coslon(iip1),coslondlon(iip1)
599      SAVE sinlon,coslon,sinlondlon,coslondlon
600      SAVE airej2,airejjm
601
602      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
603      INTEGER ifils,iq2 ! CRisi
604c
605c
606      REAL      SSUM
607
608      DATA first,testcpu/.true.,.false./
609      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
610
611      IF(first) THEN
612         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
613         first=.false.
614         do i=2,iip1
615            coslon(i)=cos(rlonv(i))
616            sinlon(i)=sin(rlonv(i))
617            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
618            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
619         ENDDO
620         coslon(1)=coslon(iip1)
621         coslondlon(1)=coslondlon(iip1)
622         sinlon(1)=sinlon(iip1)
623         sinlondlon(1)=sinlondlon(iip1)
624         airej2 = SSUM( iim, aire(iip2), 1 )
625         airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
626      ENDIF
627
628c
629
630
631      DO l = 1, llm
632c
633c   --------------------------------
634c      CALCUL EN LATITUDE
635c   --------------------------------
636
637c   On commence par calculer la valeur du traceur moyenne sur le premier cercle
638c   de latitude autour du pole (qpns pour le pole nord et qpsn pour
639c    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
640
641      DO i = 1, iim
642      airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
643      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
644      ENDDO
645      qpns   = SSUM( iim,  airescb ,1 ) / airej2
646      qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
647
648c   calcul des pentes aux points v
649
650      DO ij=1,ip1jm
651         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
652         adyqv(ij)=abs(dyqv(ij))
653      ENDDO
654
655c   calcul des pentes aux points scalaires
656
657      DO ij=iip2,ip1jm
658         dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
659         dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
660         dyqmax(ij)=pente_max*dyqmax(ij)
661      ENDDO
662
663c   calcul des pentes aux poles
664
665      DO ij=1,iip1
666         dyq(ij,l)=qpns-q(ij+iip1,l,iq)
667         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
668      ENDDO
669
670c   filtrage de la derivee
671      dyn1=0.
672      dys1=0.
673      dyn2=0.
674      dys2=0.
675      DO ij=1,iim
676         dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
677         dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
678         dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
679         dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
680      ENDDO
681      DO ij=1,iip1
682         dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
683         dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
684      ENDDO
685
686c   calcul des pentes limites aux poles
687
688      fn=1.
689      fs=1.
690      DO ij=1,iim
691         IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
692            fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
693         ENDIF
694      IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
695         fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
696         ENDIF
697      ENDDO
698      DO ij=1,iip1
699         dyq(ij,l)=fn*dyq(ij,l)
700         dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
701      ENDDO
702
703CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
704C  En memoire de dIFferents tests sur la
705C  limitation des pentes aux poles.
706CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
707C     PRINT*,dyq(1)
708C     PRINT*,dyqv(iip1+1)
709C     appn=abs(dyq(1)/dyqv(iip1+1))
710C     PRINT*,dyq(ip1jm+1)
711C     PRINT*,dyqv(ip1jm-iip1+1)
712C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
713C     DO ij=2,iim
714C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
715C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
716C     ENDDO
717C     appn=min(pente_max/appn,1.)
718C     apps=min(pente_max/apps,1.)
719C
720C
721C   cas ou on a un extremum au pole
722C
723C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
724C    &   appn=0.
725C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
726C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
727C    &   apps=0.
728C
729C   limitation des pentes aux poles
730C     DO ij=1,iip1
731C        dyq(ij)=appn*dyq(ij)
732C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
733C     ENDDO
734C
735C   test
736C      DO ij=1,iip1
737C         dyq(iip1+ij)=0.
738C         dyq(ip1jm+ij-iip1)=0.
739C      ENDDO
740C      DO ij=1,ip1jmp1
741C         dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
742C      ENDDO
743C
744C changement 10 07 96
745C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
746C    &   THEN
747C        DO ij=1,iip1
748C           dyqmax(ij)=0.
749C        ENDDO
750C     ELSE
751C        DO ij=1,iip1
752C           dyqmax(ij)=pente_max*abs(dyqv(ij))
753C        ENDDO
754C     ENDIF
755C
756C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
757C    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
758C    &THEN
759C        DO ij=ip1jm+1,ip1jmp1
760C           dyqmax(ij)=0.
761C        ENDDO
762C     ELSE
763C        DO ij=ip1jm+1,ip1jmp1
764C           dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
765C        ENDDO
766C     ENDIF
767C   fin changement 10 07 96
768CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
769
770c   calcul des pentes limitees
771
772      DO ij=iip2,ip1jm
773         IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
774            dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
775         ELSE
776            dyq(ij,l)=0.
777         ENDIF
778      ENDDO
779
780      ENDDO
781
782      DO l=1,llm
783       DO ij=1,ip1jm
784         IF( masse_adv_v(ij,l).GT.0. ) THEN
785           qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq )  +
786     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)
787     ,      /masse(ij+iip1,l,iq)))
788         ELSE
789              qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) *
790     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) )
791         ENDIF
792          qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
793       ENDDO
794      ENDDO
795
796
797! CRisi: appel récursif de l'advection sur les fils.
798! Il faut faire ça avant d'avoir mis à jour q et masse
799      !write(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq)
800   
801      if (nqfils(iq).gt.0) then 
802       do ifils=1,nqdesc(iq)
803         iq2=iqfils(ifils,iq)
804         DO l=1,llm
805         DO ij=1,ip1jmp1
806           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
807           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
808          enddo   
809         enddo
810        enddo !do ifils=1,nqdesc(iq)
811
812        do ifils=1,nqfils(iq)
813         iq2=iqfils(ifils,iq)
814         !write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2
815         call vly(Ratio,pente_max,masseq,qbyv,iq2)
816        enddo !do ifils=1,nqfils(iq)
817      endif !if (nqfils(iq).gt.0) then
818
819      DO l=1,llm
820         DO ij=iip2,ip1jm
821            newmasse=masse(ij,l,iq)
822     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
823            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
824     &         -qbyv(ij-iip1,l))/newmasse
825            masse(ij,l,iq)=newmasse
826         ENDDO
827c.-. ancienne version
828         convpn=SSUM(iim,qbyv(1,l),1)/apoln
829         convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
830         DO ij = 1,iip1
831            newmasse=masse(ij,l,iq)+convmpn*aire(ij)
832            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convpn*aire(ij))/
833     &               newmasse
834            masse(ij,l,iq)=newmasse
835         ENDDO
836         convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
837         convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
838         DO ij = ip1jm+1,ip1jmp1
839            newmasse=masse(ij,l,iq)+convmps*aire(ij)
840            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convps*aire(ij))/
841     &               newmasse
842            masse(ij,l,iq)=newmasse
843         ENDDO
844c.-. fin ancienne version
845
846c._. nouvelle version
847c        convpn=SSUM(iim,qbyv(1,l),1)
848c        convmpn=ssum(iim,masse_adv_v(1,l),1)
849c        oldmasse=ssum(iim,masse(1,l),1)
850c        newmasse=oldmasse+convmpn
851c        newq=(q(1,l)*oldmasse+convpn)/newmasse
852c        newmasse=newmasse/apoln
853c        DO ij = 1,iip1
854c           q(ij,l)=newq
855c           masse(ij,l,iq)=newmasse*aire(ij)
856c        ENDDO
857c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
858c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
859c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
860c        newmasse=oldmasse+convmps
861c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
862c        newmasse=newmasse/apols
863c        DO ij = ip1jm+1,ip1jmp1
864c           q(ij,l)=newq
865c           masse(ij,l,iq)=newmasse*aire(ij)
866c        ENDDO
867c._. fin nouvelle version
868      ENDDO
869
870      !write(*,*) 'vly 866'
871
872! retablir les fils en rapport de melange par rapport a l'air:
873      if (nqdesc(iq).gt.0) then 
874       do ifils=1,nqdesc(iq)
875         iq2=iqfils(ifils,iq) 
876         DO l=1,llm
877          DO ij=1,ip1jmp1
878            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
879          enddo
880         enddo
881        enddo !do ifils=1,nqdesc(iq)
882      endif !if (nqfils(iq).gt.0) then
883      !write(*,*) 'vly 879'
884
885      RETURN
886      END
Note: See TracBrowser for help on using the repository browser.