Ignore:
Timestamp:
Jun 14, 2015, 9:13:32 PM (9 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes -r2237:2291 into testing branch

Location:
LMDZ5/branches/testing
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/dyn3d/vlsplt.F

    r1910 r2298  
    33c
    44
    5       SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt)
     5      SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq)
     6      USE infotrac, ONLY: nqtot,nqdesc,iqfils
    67c
    78c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    3233c      REAL masse(iip1,jjp1,llm),pente_max
    3334      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
    34       REAL q(ip1jmp1,llm)
     35      REAL q(ip1jmp1,llm,nqtot)
    3536c      REAL q(iip1,jjp1,llm)
    3637      REAL w(ip1jmp1,llm),pdt
     38      INTEGER iq ! CRisi
    3739c
    3840c      Local
    … …  
    4244      INTEGER ijlqmin,iqmin,jqmin,lqmin
    4345c
    44       REAL zm(ip1jmp1,llm),newmasse
     46      REAL zm(ip1jmp1,llm,nqtot),newmasse
    4547      REAL mu(ip1jmp1,llm)
    4648      REAL mv(ip1jm,llm)
    4749      REAL mw(ip1jmp1,llm+1)
    48       REAL zq(ip1jmp1,llm),zz
     50      REAL zq(ip1jmp1,llm,nqtot),zz
    4951      REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
    5052      REAL second,temps0,temps1,temps2,temps3
    … …  
    5557      SAVE temps1,temps2,temps3
    5658      INTEGER iminn,imaxx
     59      INTEGER ifils,iq2 ! CRisi
    5760
    5861      REAL qmin,qmax
    … …  
    7982         mw(ij,llm+1)=0.
    8083      ENDDO
    81      
    82       CALL SCOPY(ijp1llm,q,1,zq,1)
    83       CALL SCOPY(ijp1llm,masse,1,zm,1)
     84           
     85      CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1)
     86      CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
     87       
     88      if (nqdesc(iq).gt.0) then 
     89        do ifils=1,nqdesc(iq)
     90          iq2=iqfils(ifils,iq)
     91          CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
     92        enddo 
     93      endif !if (nqfils(iq).gt.0) then
    8494
    8595cprint*,'Entree vlx1'
    8696c       call minmaxq(zq,qmin,qmax,'avant vlx     ')
    87       call vlx(zq,pente_max,zm,mu)
     97      call vlx(zq,pente_max,zm,mu,iq)
    8898cprint*,'Sortie vlx1'
    8999c       call minmaxq(zq,qmin,qmax,'apres vlx1    ')
    90100
    91101c print*,'Entree vly1'
    92       call vly(zq,pente_max,zm,mv)
     102
     103      call vly(zq,pente_max,zm,mv,iq)
    93104c       call minmaxq(zq,qmin,qmax,'apres vly1     ')
    94105cprint*,'Sortie vly1'
    95       call vlz(zq,pente_max,zm,mw)
     106      call vlz(zq,pente_max,zm,mw,iq)
    96107c       call minmaxq(zq,qmin,qmax,'apres vlz     ')
    97108
    98109
    99       call vly(zq,pente_max,zm,mv)
     110      call vly(zq,pente_max,zm,mv,iq)
    100111c       call minmaxq(zq,qmin,qmax,'apres vly     ')
    101112
    102113
    103       call vlx(zq,pente_max,zm,mu)
     114      call vlx(zq,pente_max,zm,mu,iq)
    104115c       call minmaxq(zq,qmin,qmax,'apres vlx2    ')
    105116       
    … …  
    107118      DO l=1,llm
    108119         DO ij=1,ip1jmp1
    109            q(ij,l)=zq(ij,l)
     120           q(ij,l,iq)=zq(ij,l,iq)
    110121         ENDDO
    111122         DO ij=1,ip1jm+1,iip1
    112             q(ij+iim,l)=q(ij,l)
    113          ENDDO
    114       ENDDO
     123            q(ij+iim,l,iq)=q(ij,l,iq)
     124         ENDDO
     125      ENDDO
     126      ! CRisi: aussi pour les fils
     127      if (nqdesc(iq).gt.0) then
     128      do ifils=1,nqdesc(iq)
     129        iq2=iqfils(ifils,iq)
     130        DO l=1,llm
     131         DO ij=1,ip1jmp1
     132           q(ij,l,iq2)=zq(ij,l,iq2)
     133         ENDDO
     134         DO ij=1,ip1jm+1,iip1
     135            q(ij+iim,l,iq2)=q(ij,l,iq2)
     136         ENDDO
     137        ENDDO
     138      enddo !do ifils=1,nqdesc(iq)   
     139      endif ! if (nqdesc(iq).gt.0) then   
    115140
    116141      RETURN
    117142      END
    118       SUBROUTINE vlx(q,pente_max,masse,u_m)
     143      RECURSIVE SUBROUTINE vlx(q,pente_max,masse,u_m,iq)
     144      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    119145
    120146c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    139165c   Arguments:
    140166c   ----------
    141       REAL masse(ip1jmp1,llm),pente_max
     167      REAL masse(ip1jmp1,llm,nqtot),pente_max
    142168      REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
    143       REAL q(ip1jmp1,llm)
     169      REAL q(ip1jmp1,llm,nqtot)
    144170      REAL w(ip1jmp1,llm)
     171      INTEGER iq ! CRisi
    145172c
    146173c      Local
    … …  
    155182      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
    156183      REAL u_mq(ip1jmp1,llm)
     184
     185      ! CRisi
     186      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
     187      INTEGER ifils,iq2 ! CRisi
    157188
    158189      Logical extremum,first,testcpu
    … …  
    188219         DO l = 1, llm
    189220            DO ij=iip2,ip1jm-1
    190                dxqu(ij)=q(ij+1,l)-q(ij,l)
     221               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    191222c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
    192 c              sigu(ij)=u_m(ij,l)/masse(ij,l)
     223c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
    193224            ENDDO
    194225            DO ij=iip1+iip1,ip1jm,iip1
    … …  
    243274         DO l = 1, llm
    244275            DO ij=iip2,ip1jm-1
    245                dxqu(ij)=q(ij+1,l)-q(ij,l)
     276               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    246277            ENDDO
    247278            DO ij=iip1+iip1,ip1jm,iip1
    … …  
    285316      DO l=1,llm
    286317       DO ij=iip2,ip1jm-1
    287           zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
    288      ,                     1.+u_m(ij,l)/masse(ij+1,l),
     318          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
     319     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
    289320     ,                     u_m(ij,l))
    290321          zdum(ij,l)=0.5*zdum(ij,l)
    291322          u_mq(ij,l)=cvmgp(
    292      ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
    293      ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
     323     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
     324     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
    294325     ,                u_m(ij,l))
    295326          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
    … …  
    303334      DO l=1,llm
    304335       DO ij=iip2,ip1jm-1
    305 c       print*,'masse(',ij,')=',masse(ij,l)
     336c       print*,'masse(',ij,')=',masse(ij,l,iq)
    306337          IF (u_m(ij,l).gt.0.) THEN
    307              zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
    308              u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
     338             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
     339             u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l))
    309340          ELSE
    310              zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
    311              u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
     341             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
     342             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq)
     343     &           -0.5*zdum(ij,l)*dxq(ij+1,l))
    312344          ENDIF
    313345       ENDDO
    … …  
    379411                     i=ijq-(j-1)*iip1
    380412c   accumulation pour les mailles completements advectees
    381                      do while(zu_m.gt.masse(ijq,l))
    382                         u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
    383                         zu_m=zu_m-masse(ijq,l)
     413                     do while(zu_m.gt.masse(ijq,l,iq))
     414                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
     415     &                          *masse(ijq,l,iq)
     416                        zu_m=zu_m-masse(ijq,l,iq)
    384417                        i=mod(i-2+iim,iim)+1
    385418                        ijq=(j-1)*iip1+i
    … …  
    387420c   ajout de la maille non completement advectee
    388421                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
    389      &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
     422     &                  (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
     423     &                  *dxq(ijq,l))
    390424                  ELSE
    391425                     ijq=ij+1
    392426                     i=ijq-(j-1)*iip1
    393427c   accumulation pour les mailles completements advectees
    394                      do while(-zu_m.gt.masse(ijq,l))
    395                         u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
    396                         zu_m=zu_m+masse(ijq,l)
     428                     do while(-zu_m.gt.masse(ijq,l,iq))
     429                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
     430     &                          *masse(ijq,l,iq)
     431                        zu_m=zu_m+masse(ijq,l,iq)
    397432                        i=mod(i,iim)+1
    398433                        ijq=(j-1)*iip1+i
    399434                     ENDDO
    400435c   ajout de la maille non completement advectee
    401                      u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
    402      &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
     436                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
     437     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
    403438                  ENDIF
    404439               ENDDO
    … …  
    417452      ENDDO
    418453
     454! CRisi: appel récursif de l'advection sur les fils.
     455! Il faut faire ça avant d'avoir mis à jour q et masse
     456      !write(*,*) 'vlsplt 326: iq,nqfils(iq)=',iq,nqfils(iq)
     457     
     458      if (nqdesc(iq).gt.0) then 
     459       do ifils=1,nqdesc(iq)
     460         iq2=iqfils(ifils,iq)
     461         DO l=1,llm
     462          DO ij=iip2,ip1jm
     463           ! On a besoin de q et masse seulement entre iip2 et ip1jm
     464           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     465           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     466          enddo   
     467         enddo
     468        enddo !do ifils=1,nqdesc(iq)
     469        do ifils=1,nqfils(iq)
     470         iq2=iqfils(ifils,iq)
     471         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
     472        enddo !do ifils=1,nqfils(iq)
     473      endif !if (nqfils(iq).gt.0) then
     474! end CRisi
     475
    419476
    420477c   calcul des tENDances
    … …  
    422479      DO l=1,llm
    423480         DO ij=iip2+1,ip1jm
    424             new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
    425             q(ij,l)=(q(ij,l)*masse(ij,l)+
     481            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
     482            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
    426483     &      u_mq(ij-1,l)-u_mq(ij,l))
    427484     &      /new_m
    428             masse(ij,l)=new_m
     485            masse(ij,l,iq)=new_m
    429486         ENDDO
    430487c   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
    431488         DO ij=iip1+iip1,ip1jm,iip1
    432             q(ij-iim,l)=q(ij,l)
    433             masse(ij-iim,l)=masse(ij,l)
    434          ENDDO
    435       ENDDO
     489            q(ij-iim,l,iq)=q(ij,l,iq)
     490            masse(ij-iim,l,iq)=masse(ij,l,iq)
     491         ENDDO
     492      ENDDO
     493
     494      ! retablir les fils en rapport de melange par rapport a l'air:
     495      ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
     496      ! puis on boucle en longitude
     497      if (nqdesc(iq).gt.0) then 
     498       do ifils=1,nqdesc(iq)
     499         iq2=iqfils(ifils,iq) 
     500         DO l=1,llm
     501          DO ij=iip2+1,ip1jm
     502            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     503          enddo
     504          DO ij=iip1+iip1,ip1jm,iip1
     505             q(ij-iim,l,iq2)=q(ij,l,iq2)
     506          enddo ! DO ij=ijb+iip1-1,ije,iip1
     507         enddo !DO l=1,llm
     508        enddo !do ifils=1,nqdesc(iq)
     509      endif !if (nqfils(iq).gt.0) then
     510
    436511c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
    437512c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
    … …  
    440515      RETURN
    441516      END
    442       SUBROUTINE vly(q,pente_max,masse,masse_adv_v)
     517      RECURSIVE SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq)
     518      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    443519c
    444520c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    464540c   Arguments:
    465541c   ----------
    466       REAL masse(ip1jmp1,llm),pente_max
     542      REAL masse(ip1jmp1,llm,nqtot),pente_max
    467543      REAL masse_adv_v( ip1jm,llm)
    468       REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
     544      REAL q(ip1jmp1,llm,nqtot), dq( ip1jmp1,llm)
     545      INTEGER iq ! CRisi
    469546c
    470547c      Local
    … …  
    491568      SAVE sinlon,coslon,sinlondlon,coslondlon
    492569      SAVE airej2,airejjm
     570
     571      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
     572      INTEGER ifils,iq2 ! CRisi
     573
    493574c
    494575c
    … …  
    497578      DATA first,testcpu/.true.,.false./
    498579      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
     580
     581      !write(*,*) 'vly 578: entree, iq=',iq
    499582
    500583      IF(first) THEN
    … …  
    529612
    530613      DO i = 1, iim
    531       airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
    532       airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
     614      airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
     615      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
    533616      ENDDO
    534617      qpns   = SSUM( iim,  airescb ,1 ) / airej2
    … …  
    538621
    539622      DO ij=1,ip1jm
    540          dyqv(ij)=q(ij,l)-q(ij+iip1,l)
     623         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
    541624         adyqv(ij)=abs(dyqv(ij))
    542625      ENDDO
    … …  
    553636
    554637      DO ij=1,iip1
    555          dyq(ij,l)=qpns-q(ij+iip1,l)
    556          dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
     638         dyq(ij,l)=qpns-q(ij+iip1,l,iq)
     639         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
    557640      ENDDO
    558641
    … …  
    675758      ENDDO
    676759
     760      !write(*,*) 'vly 756'
    677761      DO l=1,llm
    678762       DO ij=1,ip1jm
    679763          IF(masse_adv_v(ij,l).gt.0) THEN
    680               qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
    681      ,                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
     764              qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)*
     765     ,                   0.5*(1.-masse_adv_v(ij,l)
     766     ,                   /masse(ij+iip1,l,iq))
    682767          ELSE
    683               qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
    684      ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
     768              qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)*
     769     ,                   0.5*(1.+masse_adv_v(ij,l)
     770     ,                   /masse(ij,l,iq))
    685771          ENDIF
    686772          qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
    … …  
    688774      ENDDO
    689775
     776! CRisi: appel récursif de l'advection sur les fils.
     777! Il faut faire ça avant d'avoir mis à jour q et masse
     778      !write(*,*) 'vly 689: iq,nqfils(iq)=',iq,nqfils(iq)
     779   
     780      if (nqfils(iq).gt.0) then 
     781       do ifils=1,nqdesc(iq)
     782         iq2=iqfils(ifils,iq)
     783         DO l=1,llm
     784         DO ij=1,ip1jmp1
     785           ! attention, chaque fils doit avoir son masseq, sinon, le 1er
     786           ! fils ecrase le masseq de ses freres.
     787           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     788           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
     789          enddo   
     790         enddo
     791        enddo !do ifils=1,nqdesc(iq)
     792
     793        do ifils=1,nqfils(iq)
     794         iq2=iqfils(ifils,iq)
     795         call vly(Ratio,pente_max,masseq,qbyv,iq2)
     796        enddo !do ifils=1,nqfils(iq)
     797      endif !if (nqfils(iq).gt.0) then
    690798
    691799      DO l=1,llm
    692800         DO ij=iip2,ip1jm
    693             newmasse=masse(ij,l)
     801            newmasse=masse(ij,l,iq)
    694802     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
    695             q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
    696      &         /newmasse
    697             masse(ij,l)=newmasse
     803            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
     804     &         -qbyv(ij-iip1,l))/newmasse
     805            masse(ij,l,iq)=newmasse
    698806         ENDDO
    699807c.-. ancienne version
    … …  
    703811         convpn=SSUM(iim,qbyv(1,l),1)
    704812         convmpn=ssum(iim,masse_adv_v(1,l),1)
    705          massepn=ssum(iim,masse(1,l),1)
     813         massepn=ssum(iim,masse(1,l,iq),1)
    706814         qpn=0.
    707815         do ij=1,iim
    708             qpn=qpn+masse(ij,l)*q(ij,l)
     816            qpn=qpn+masse(ij,l,iq)*q(ij,l,iq)
    709817         enddo
    710818         qpn=(qpn+convpn)/(massepn+convmpn)
    711819         do ij=1,iip1
    712             q(ij,l)=qpn
     820            q(ij,l,iq)=qpn
    713821         enddo
    714822
    … …  
    718826         convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    719827         convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
    720          masseps=ssum(iim, masse(ip1jm+1,l),1)
     828         masseps=ssum(iim, masse(ip1jm+1,l,iq),1)
    721829         qps=0.
    722830         do ij = ip1jm+1,ip1jmp1-1
    723             qps=qps+masse(ij,l)*q(ij,l)
     831            qps=qps+masse(ij,l,iq)*q(ij,l,iq)
    724832         enddo
    725833         qps=(qps+convps)/(masseps+convmps)
    726834         do ij=ip1jm+1,ip1jmp1
    727             q(ij,l)=qps
     835            q(ij,l,iq)=qps
    728836         enddo
    729837
    … …  
    739847c        DO ij = 1,iip1
    740848c           q(ij,l)=newq
    741 c           masse(ij,l)=newmasse*aire(ij)
     849c           masse(ij,l,iq)=newmasse*aire(ij)
    742850c        ENDDO
    743851c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    … …  
    749857c        DO ij = ip1jm+1,ip1jmp1
    750858c           q(ij,l)=newq
    751 c           masse(ij,l)=newmasse*aire(ij)
     859c           masse(ij,l,iq)=newmasse*aire(ij)
    752860c        ENDDO
    753861c._. fin nouvelle version
    754862      ENDDO
     863 
     864! retablir les fils en rapport de melange par rapport a l'air:
     865      if (nqfils(iq).gt.0) then 
     866       do ifils=1,nqdesc(iq)
     867         iq2=iqfils(ifils,iq) 
     868         DO l=1,llm
     869          DO ij=1,ip1jmp1
     870            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     871          enddo
     872         enddo
     873        enddo !do ifils=1,nqdesc(iq)
     874      endif !if (nqfils(iq).gt.0) then
     875
     876      !write(*,*) 'vly 853: sortie'
    755877
    756878      RETURN
    757879      END
    758       SUBROUTINE vlz(q,pente_max,masse,w)
     880      RECURSIVE SUBROUTINE vlz(q,pente_max,masse,w,iq)
     881      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    759882c
    760883c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    779902c   Arguments:
    780903c   ----------
    781       REAL masse(ip1jmp1,llm),pente_max
    782       REAL q(ip1jmp1,llm)
     904      REAL masse(ip1jmp1,llm,nqtot),pente_max
     905      REAL q(ip1jmp1,llm,nqtot)
    783906      REAL w(ip1jmp1,llm+1)
     907      INTEGER iq
    784908c
    785909c      Local
    … …  
    792916      REAL dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm),dzqmax
    793917      REAL sigw
     918
     919      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
     920      INTEGER ifils,iq2 ! CRisi
    794921
    795922      LOGICAL testcpu
    … …  
    805932c    On oriente tout dans le sens de la pression c'est a dire dans le
    806933c    sens de W
     934
     935      !write(*,*) 'vlz 923: entree'
    807936
    808937#ifdef BIDON
    … …  
    813942      DO l=2,llm
    814943         DO ij=1,ip1jmp1
    815             dzqw(ij,l)=q(ij,l-1)-q(ij,l)
     944            dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq)
    816945            adzqw(ij,l)=abs(dzqw(ij,l))
    817946         ENDDO
    … …  
    835964      ENDDO
    836965
     966      !write(*,*) 'vlz 954'
    837967      DO ij=1,ip1jmp1
    838968         dzq(ij,1)=0.
    … …  
    851981c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
    852982
     983       !write(*,*) 'vlz 969'
    853984       DO l = 1,llm-1
    854985         do  ij = 1,ip1jmp1
    855986          IF(w(ij,l+1).gt.0.) THEN
    856              sigw=w(ij,l+1)/masse(ij,l+1)
    857              wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
     987             sigw=w(ij,l+1)/masse(ij,l+1,iq)
     988             wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1,iq)
     989     &           +0.5*(1.-sigw)*dzq(ij,l+1))
    858990          ELSE
    859              sigw=w(ij,l+1)/masse(ij,l)
    860              wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
     991             sigw=w(ij,l+1)/masse(ij,l,iq)
     992             wq(ij,l+1)=w(ij,l+1)*(q(ij,l,iq)-0.5*(1.+sigw)*dzq(ij,l))
    861993          ENDIF
    862994         ENDDO
    … …  
    8681000       ENDDO
    8691001
     1002! CRisi: appel récursif de l'advection sur les fils.
     1003! Il faut faire ça avant d'avoir mis à jour q et masse
     1004      !write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq)
     1005      if (nqfils(iq).gt.0) then 
     1006       do ifils=1,nqdesc(iq)
     1007         iq2=iqfils(ifils,iq)
     1008         DO l=1,llm
     1009          DO ij=1,ip1jmp1
     1010           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     1011           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)       
     1012          enddo   
     1013         enddo
     1014        enddo !do ifils=1,nqdesc(iq)
     1015       
     1016        do ifils=1,nqfils(iq)
     1017         iq2=iqfils(ifils,iq)         
     1018         call vlz(Ratio,pente_max,masseq,wq,iq2)
     1019        enddo !do ifils=1,nqfils(iq)
     1020      endif !if (nqfils(iq).gt.0) then
     1021! end CRisi 
     1022
    8701023      DO l=1,llm
    8711024         DO ij=1,ip1jmp1
    872             newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
    873             q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
     1025            newmasse=masse(ij,l,iq)+w(ij,l+1)-w(ij,l)
     1026            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+wq(ij,l+1)-wq(ij,l))
    8741027     &         /newmasse
    875             masse(ij,l)=newmasse
    876          ENDDO
    877       ENDDO
    878 
     1028            masse(ij,l,iq)=newmasse
     1029         ENDDO
     1030      ENDDO
     1031
     1032! retablir les fils en rapport de melange par rapport a l'air:
     1033      if (nqfils(iq).gt.0) then 
     1034       do ifils=1,nqdesc(iq)
     1035         iq2=iqfils(ifils,iq) 
     1036         DO l=1,llm
     1037          DO ij=1,ip1jmp1
     1038            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     1039          enddo
     1040         enddo
     1041        enddo !do ifils=1,nqdesc(iq)
     1042      endif !if (nqfils(iq).gt.0) then
     1043      !write(*,*) 'vlsplt 1032'
    8791044
    8801045      RETURN
Note: See TracChangeset for help on using the changeset viewer.