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

Last change on this file since 3913 was 2603, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn logic.h into module logic_mod.F90
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.0 KB
Line 
1c
2c $Id: vlspltqs.F 2603 2016-07-25 09:31:56Z idelkadi $
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     
26      USE comconst_mod, ONLY: cpp
27     
28      IMPLICIT NONE
29c
30      include "dimensions.h"
31      include "paramet.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"
196c
197c
198c   Arguments:
199c   ----------
200      REAL masse(ip1jmp1,llm,nqtot),pente_max
201      REAL u_m( ip1jmp1,llm )
202      REAL q(ip1jmp1,llm,nqtot)
203      REAL qsat(ip1jmp1,llm)
204      INTEGER iq ! CRisi
205c
206c      Local
207c   ---------
208c
209      INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
210      INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
211c
212      REAL new_m,zu_m,zdum(ip1jmp1,llm)
213      REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
214      REAL zz(ip1jmp1)
215      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
216      REAL u_mq(ip1jmp1,llm)
217
218      ! CRisi
219      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
220      INTEGER ifils,iq2 ! CRisi
221
222      Logical first,testcpu
223      SAVE first,testcpu
224
225      REAL      SSUM
226      REAL temps0,temps1,temps2,temps3,temps4,temps5
227      SAVE temps0,temps1,temps2,temps3,temps4,temps5
228
229
230      DATA first,testcpu/.true.,.false./
231
232      IF(first) THEN
233         temps1=0.
234         temps2=0.
235         temps3=0.
236         temps4=0.
237         temps5=0.
238         first=.false.
239      ENDIF
240
241c   calcul de la pente a droite et a gauche de la maille
242
243
244      IF (pente_max.gt.-1.e-5) THEN
245c     IF (pente_max.gt.10) THEN
246
247c   calcul des pentes avec limitation, Van Leer scheme I:
248c   -----------------------------------------------------
249
250c   calcul de la pente aux points u
251         DO l = 1, llm
252            DO ij=iip2,ip1jm-1
253               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
254c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
255c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
256            ENDDO
257            DO ij=iip1+iip1,ip1jm,iip1
258               dxqu(ij)=dxqu(ij-iim)
259c              sigu(ij)=sigu(ij-iim)
260            ENDDO
261
262            DO ij=iip2,ip1jm
263               adxqu(ij)=abs(dxqu(ij))
264            ENDDO
265
266c   calcul de la pente maximum dans la maille en valeur absolue
267
268            DO ij=iip2+1,ip1jm
269               dxqmax(ij,l)=pente_max*
270     ,      min(adxqu(ij-1),adxqu(ij))
271c limitation subtile
272c    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
273         
274
275            ENDDO
276
277            DO ij=iip1+iip1,ip1jm,iip1
278               dxqmax(ij-iim,l)=dxqmax(ij,l)
279            ENDDO
280
281            DO ij=iip2+1,ip1jm
282#ifdef CRAY
283               dxq(ij,l)=
284     ,         cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
285#else
286               IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
287                  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
288               ELSE
289c   extremum local
290                  dxq(ij,l)=0.
291               ENDIF
292#endif
293               dxq(ij,l)=0.5*dxq(ij,l)
294               dxq(ij,l)=
295     ,         sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
296            ENDDO
297
298         ENDDO ! l=1,llm
299
300      ELSE ! (pente_max.lt.-1.e-5)
301
302c   Pentes produits:
303c   ----------------
304
305         DO l = 1, llm
306            DO ij=iip2,ip1jm-1
307               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
308            ENDDO
309            DO ij=iip1+iip1,ip1jm,iip1
310               dxqu(ij)=dxqu(ij-iim)
311            ENDDO
312
313            DO ij=iip2+1,ip1jm
314               zz(ij)=dxqu(ij-1)*dxqu(ij)
315               zz(ij)=zz(ij)+zz(ij)
316               IF(zz(ij).gt.0) THEN
317                  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
318               ELSE
319c   extremum local
320                  dxq(ij,l)=0.
321               ENDIF
322            ENDDO
323
324         ENDDO
325
326      ENDIF ! (pente_max.lt.-1.e-5)
327
328c   bouclage de la pente en iip1:
329c   -----------------------------
330
331      DO l=1,llm
332         DO ij=iip1+iip1,ip1jm,iip1
333            dxq(ij-iim,l)=dxq(ij,l)
334         ENDDO
335
336         DO ij=1,ip1jmp1
337            iadvplus(ij,l)=0
338         ENDDO
339
340      ENDDO
341
342
343c   calcul des flux a gauche et a droite
344
345#ifdef CRAY
346c--pas encore modification sur Qsat
347      DO l=1,llm
348       DO ij=iip2,ip1jm-1
349          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
350     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
351     ,                     u_m(ij,l))
352          zdum(ij,l)=0.5*zdum(ij,l)
353          u_mq(ij,l)=cvmgp(
354     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
355     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
356     ,                u_m(ij,l))
357          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
358       ENDDO
359      ENDDO
360#else
361c   on cumule le flux correspondant a toutes les mailles dont la masse
362c   au travers de la paroi pENDant le pas de temps.
363c   le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
364      DO l=1,llm
365       DO ij=iip2,ip1jm-1
366          IF (u_m(ij,l).gt.0.) THEN
367             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
368             u_mq(ij,l)=u_m(ij,l)*
369     $         min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
370          ELSE
371             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
372             u_mq(ij,l)=u_m(ij,l)*
373     $         min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
374          ENDIF
375       ENDDO
376      ENDDO
377#endif
378
379
380c   detection des points ou on advecte plus que la masse de la
381c   maille
382      DO l=1,llm
383         DO ij=iip2,ip1jm-1
384            IF(zdum(ij,l).lt.0) THEN
385               iadvplus(ij,l)=1
386               u_mq(ij,l)=0.
387            ENDIF
388         ENDDO
389      ENDDO
390      DO l=1,llm
391       DO ij=iip1+iip1,ip1jm,iip1
392          iadvplus(ij,l)=iadvplus(ij-iim,l)
393       ENDDO
394      ENDDO
395
396
397
398c   traitement special pour le cas ou on advecte en longitude plus que le
399c   contenu de la maille.
400c   cette partie est mal vectorisee.
401
402c   pas d'influence de la pression saturante (pour l'instant)
403
404c  calcul du nombre de maille sur lequel on advecte plus que la maille.
405
406      n0=0
407      DO l=1,llm
408         nl(l)=0
409         DO ij=iip2,ip1jm
410            nl(l)=nl(l)+iadvplus(ij,l)
411         ENDDO
412         n0=n0+nl(l)
413      ENDDO
414
415      IF(n0.gt.0) THEN
416ccc      PRINT*,'Nombre de points pour lesquels on advect plus que le'
417ccc     &       ,'contenu de la maille : ',n0
418
419         DO l=1,llm
420            IF(nl(l).gt.0) THEN
421               iju=0
422c   indicage des mailles concernees par le traitement special
423               DO ij=iip2,ip1jm
424                  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
425                     iju=iju+1
426                     indu(iju)=ij
427                  ENDIF
428               ENDDO
429               niju=iju
430c              PRINT*,'niju,nl',niju,nl(l)
431
432c  traitement des mailles
433               DO iju=1,niju
434                  ij=indu(iju)
435                  j=(ij-1)/iip1+1
436                  zu_m=u_m(ij,l)
437                  u_mq(ij,l)=0.
438                  IF(zu_m.gt.0.) THEN
439                     ijq=ij
440                     i=ijq-(j-1)*iip1
441c   accumulation pour les mailles completements advectees
442                     do while(zu_m.gt.masse(ijq,l,iq))
443                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
444     &                          *masse(ijq,l,iq)
445                        zu_m=zu_m-masse(ijq,l,iq)
446                        i=mod(i-2+iim,iim)+1
447                        ijq=(j-1)*iip1+i
448                     ENDDO
449c   ajout de la maille non completement advectee
450                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
451     &                  (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
452     &                  *dxq(ijq,l))
453                  ELSE
454                     ijq=ij+1
455                     i=ijq-(j-1)*iip1
456c   accumulation pour les mailles completements advectees
457                     do while(-zu_m.gt.masse(ijq,l,iq))
458                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
459     &                          *masse(ijq,l,iq)
460                        zu_m=zu_m+masse(ijq,l,iq)
461                        i=mod(i,iim)+1
462                        ijq=(j-1)*iip1+i
463                     ENDDO
464c   ajout de la maille non completement advectee
465                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
466     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
467                  ENDIF
468               ENDDO
469            ENDIF
470         ENDDO
471      ENDIF  ! n0.gt.0
472
473
474
475c   bouclage en latitude
476
477      DO l=1,llm
478        DO ij=iip1+iip1,ip1jm,iip1
479           u_mq(ij,l)=u_mq(ij-iim,l)
480        ENDDO
481      ENDDO
482
483! CRisi: appel récursif de l'advection sur les fils.
484! Il faut faire ça avant d'avoir mis à jour q et masse
485      !write(*,*) 'vlspltqs 326: iq,nqfils(iq)=',iq,nqfils(iq)
486     
487      if (nqfils(iq).gt.0) then 
488       do ifils=1,nqdesc(iq)
489         iq2=iqfils(ifils,iq)
490         DO l=1,llm
491          DO ij=iip2,ip1jm
492           ! On a besoin de q et masse seulement entre iip2 et ip1jm
493           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
494           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
495          enddo   
496         enddo
497        enddo !do ifils=1,nqdesc(iq)
498        do ifils=1,nqfils(iq)
499         iq2=iqfils(ifils,iq)
500         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
501        enddo !do ifils=1,nqfils(iq)
502      endif !if (nqfils(iq).gt.0) then
503! end CRisi
504
505c   calcul des tendances
506
507      DO l=1,llm
508         DO ij=iip2+1,ip1jm
509            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
510            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
511     &      u_mq(ij-1,l)-u_mq(ij,l))
512     &      /new_m
513            masse(ij,l,iq)=new_m
514         ENDDO
515c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
516         DO ij=iip1+iip1,ip1jm,iip1
517            q(ij-iim,l,iq)=q(ij,l,iq)
518            masse(ij-iim,l,iq)=masse(ij,l,iq)
519         ENDDO
520      ENDDO
521
522      ! retablir les fils en rapport de melange par rapport a l'air:
523      ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
524      ! puis on boucle en longitude
525      if (nqdesc(iq).gt.0) then 
526       do ifils=1,nqdesc(iq)
527         iq2=iqfils(ifils,iq) 
528         DO l=1,llm
529          DO ij=iip2+1,ip1jm
530            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
531          enddo
532          DO ij=iip1+iip1,ip1jm,iip1
533             q(ij-iim,l,iq2)=q(ij,l,iq2)
534          enddo ! DO ij=ijb+iip1-1,ije,iip1
535         enddo !DO l=1,llm
536        enddo !do ifils=1,nqdesc(iq)
537      endif !if (nqfils(iq).gt.0) then
538
539c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
540c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
541
542
543      RETURN
544      END
545      SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq)
546      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
547c
548c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
549c
550c    ********************************************************************
551c     Shema  d'advection " pseudo amont " .
552c    ********************************************************************
553c     q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
554c     qsat             est   un argument de sortie pour le s-pg ....
555c
556c
557c   --------------------------------------------------------------------
558     
559      USE comconst_mod, ONLY: pi
560     
561      IMPLICIT NONE
562c
563      include "dimensions.h"
564      include "paramet.h"
565      include "comgeom.h"
566c
567c
568c   Arguments:
569c   ----------
570      REAL masse(ip1jmp1,llm,nqtot),pente_max
571      REAL masse_adv_v( ip1jm,llm)
572      REAL q(ip1jmp1,llm,nqtot)
573      REAL qsat(ip1jmp1,llm)
574      INTEGER iq ! CRisi
575c
576c      Local
577c   ---------
578c
579      INTEGER i,ij,l
580c
581      REAL airej2,airejjm,airescb(iim),airesch(iim)
582      REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
583      REAL adyqv(ip1jm),dyqmax(ip1jmp1)
584      REAL qbyv(ip1jm,llm)
585
586      REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
587c     REAL newq,oldmasse
588      Logical first,testcpu
589      REAL temps0,temps1,temps2,temps3,temps4,temps5
590      SAVE temps0,temps1,temps2,temps3,temps4,temps5
591      SAVE first,testcpu
592
593      REAL convpn,convps,convmpn,convmps
594      REAL sinlon(iip1),sinlondlon(iip1)
595      REAL coslon(iip1),coslondlon(iip1)
596      SAVE sinlon,coslon,sinlondlon,coslondlon
597      SAVE airej2,airejjm
598
599      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
600      INTEGER ifils,iq2 ! CRisi
601c
602c
603      REAL      SSUM
604
605      DATA first,testcpu/.true.,.false./
606      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
607
608      IF(first) THEN
609         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
610         first=.false.
611         do i=2,iip1
612            coslon(i)=cos(rlonv(i))
613            sinlon(i)=sin(rlonv(i))
614            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
615            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
616         ENDDO
617         coslon(1)=coslon(iip1)
618         coslondlon(1)=coslondlon(iip1)
619         sinlon(1)=sinlon(iip1)
620         sinlondlon(1)=sinlondlon(iip1)
621         airej2 = SSUM( iim, aire(iip2), 1 )
622         airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
623      ENDIF
624
625c
626
627
628      DO l = 1, llm
629c
630c   --------------------------------
631c      CALCUL EN LATITUDE
632c   --------------------------------
633
634c   On commence par calculer la valeur du traceur moyenne sur le premier cercle
635c   de latitude autour du pole (qpns pour le pole nord et qpsn pour
636c    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
637
638      DO i = 1, iim
639      airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
640      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
641      ENDDO
642      qpns   = SSUM( iim,  airescb ,1 ) / airej2
643      qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
644
645c   calcul des pentes aux points v
646
647      DO ij=1,ip1jm
648         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
649         adyqv(ij)=abs(dyqv(ij))
650      ENDDO
651
652c   calcul des pentes aux points scalaires
653
654      DO ij=iip2,ip1jm
655         dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
656         dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
657         dyqmax(ij)=pente_max*dyqmax(ij)
658      ENDDO
659
660c   calcul des pentes aux poles
661
662      DO ij=1,iip1
663         dyq(ij,l)=qpns-q(ij+iip1,l,iq)
664         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
665      ENDDO
666
667c   filtrage de la derivee
668      dyn1=0.
669      dys1=0.
670      dyn2=0.
671      dys2=0.
672      DO ij=1,iim
673         dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
674         dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
675         dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
676         dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
677      ENDDO
678      DO ij=1,iip1
679         dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
680         dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
681      ENDDO
682
683c   calcul des pentes limites aux poles
684
685      fn=1.
686      fs=1.
687      DO ij=1,iim
688         IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
689            fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
690         ENDIF
691      IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
692         fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
693         ENDIF
694      ENDDO
695      DO ij=1,iip1
696         dyq(ij,l)=fn*dyq(ij,l)
697         dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
698      ENDDO
699
700CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
701C  En memoire de dIFferents tests sur la
702C  limitation des pentes aux poles.
703CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
704C     PRINT*,dyq(1)
705C     PRINT*,dyqv(iip1+1)
706C     appn=abs(dyq(1)/dyqv(iip1+1))
707C     PRINT*,dyq(ip1jm+1)
708C     PRINT*,dyqv(ip1jm-iip1+1)
709C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
710C     DO ij=2,iim
711C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
712C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
713C     ENDDO
714C     appn=min(pente_max/appn,1.)
715C     apps=min(pente_max/apps,1.)
716C
717C
718C   cas ou on a un extremum au pole
719C
720C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
721C    &   appn=0.
722C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
723C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
724C    &   apps=0.
725C
726C   limitation des pentes aux poles
727C     DO ij=1,iip1
728C        dyq(ij)=appn*dyq(ij)
729C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
730C     ENDDO
731C
732C   test
733C      DO ij=1,iip1
734C         dyq(iip1+ij)=0.
735C         dyq(ip1jm+ij-iip1)=0.
736C      ENDDO
737C      DO ij=1,ip1jmp1
738C         dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
739C      ENDDO
740C
741C changement 10 07 96
742C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
743C    &   THEN
744C        DO ij=1,iip1
745C           dyqmax(ij)=0.
746C        ENDDO
747C     ELSE
748C        DO ij=1,iip1
749C           dyqmax(ij)=pente_max*abs(dyqv(ij))
750C        ENDDO
751C     ENDIF
752C
753C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
754C    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
755C    &THEN
756C        DO ij=ip1jm+1,ip1jmp1
757C           dyqmax(ij)=0.
758C        ENDDO
759C     ELSE
760C        DO ij=ip1jm+1,ip1jmp1
761C           dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
762C        ENDDO
763C     ENDIF
764C   fin changement 10 07 96
765CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
766
767c   calcul des pentes limitees
768
769      DO ij=iip2,ip1jm
770         IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
771            dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
772         ELSE
773            dyq(ij,l)=0.
774         ENDIF
775      ENDDO
776
777      ENDDO
778
779      DO l=1,llm
780       DO ij=1,ip1jm
781         IF( masse_adv_v(ij,l).GT.0. ) THEN
782           qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq )  +
783     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)
784     ,      /masse(ij+iip1,l,iq)))
785         ELSE
786              qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) *
787     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) )
788         ENDIF
789          qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
790       ENDDO
791      ENDDO
792
793
794! CRisi: appel récursif de l'advection sur les fils.
795! Il faut faire ça avant d'avoir mis à jour q et masse
796      !write(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq)
797   
798      if (nqfils(iq).gt.0) then 
799       do ifils=1,nqdesc(iq)
800         iq2=iqfils(ifils,iq)
801         DO l=1,llm
802         DO ij=1,ip1jmp1
803           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
804           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
805          enddo   
806         enddo
807        enddo !do ifils=1,nqdesc(iq)
808
809        do ifils=1,nqfils(iq)
810         iq2=iqfils(ifils,iq)
811         !write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2
812         call vly(Ratio,pente_max,masseq,qbyv,iq2)
813        enddo !do ifils=1,nqfils(iq)
814      endif !if (nqfils(iq).gt.0) then
815
816      DO l=1,llm
817         DO ij=iip2,ip1jm
818            newmasse=masse(ij,l,iq)
819     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
820            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
821     &         -qbyv(ij-iip1,l))/newmasse
822            masse(ij,l,iq)=newmasse
823         ENDDO
824c.-. ancienne version
825         convpn=SSUM(iim,qbyv(1,l),1)/apoln
826         convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
827         DO ij = 1,iip1
828            newmasse=masse(ij,l,iq)+convmpn*aire(ij)
829            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convpn*aire(ij))/
830     &               newmasse
831            masse(ij,l,iq)=newmasse
832         ENDDO
833         convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
834         convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
835         DO ij = ip1jm+1,ip1jmp1
836            newmasse=masse(ij,l,iq)+convmps*aire(ij)
837            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convps*aire(ij))/
838     &               newmasse
839            masse(ij,l,iq)=newmasse
840         ENDDO
841c.-. fin ancienne version
842
843c._. nouvelle version
844c        convpn=SSUM(iim,qbyv(1,l),1)
845c        convmpn=ssum(iim,masse_adv_v(1,l),1)
846c        oldmasse=ssum(iim,masse(1,l),1)
847c        newmasse=oldmasse+convmpn
848c        newq=(q(1,l)*oldmasse+convpn)/newmasse
849c        newmasse=newmasse/apoln
850c        DO ij = 1,iip1
851c           q(ij,l)=newq
852c           masse(ij,l,iq)=newmasse*aire(ij)
853c        ENDDO
854c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
855c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
856c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
857c        newmasse=oldmasse+convmps
858c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
859c        newmasse=newmasse/apols
860c        DO ij = ip1jm+1,ip1jmp1
861c           q(ij,l)=newq
862c           masse(ij,l,iq)=newmasse*aire(ij)
863c        ENDDO
864c._. fin nouvelle version
865      ENDDO
866
867      !write(*,*) 'vly 866'
868
869! retablir les fils en rapport de melange par rapport a l'air:
870      if (nqdesc(iq).gt.0) then 
871       do ifils=1,nqdesc(iq)
872         iq2=iqfils(ifils,iq) 
873         DO l=1,llm
874          DO ij=1,ip1jmp1
875            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
876          enddo
877         enddo
878        enddo !do ifils=1,nqdesc(iq)
879      endif !if (nqfils(iq).gt.0) then
880      !write(*,*) 'vly 879'
881
882      RETURN
883      END
Note: See TracBrowser for help on using the repository browser.