Ignore:
Timestamp:
May 7, 2015, 5:45:04 PM (9 years ago)
Author:
crisi
Message:

Adding isotopes in the dynamics and more generally tracers of tracers.
CRisi

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3dmem/vlspltqs_loc.F

    r1907 r2270  
    1       SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x)
     1      SUBROUTINE vlxqs_loc(q,pente_max,masse,u_m,qsat,ijb_x,ije_x,iq)
    22c
    33c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    44c
    55c    ********************************************************************
    6 c     Shema  d'advection " pseudo amont " .
     6c     Shema  d''advection " pseudo amont " .
    77c    ********************************************************************
    88c
    99c   --------------------------------------------------------------------
    1010      USE parallel_lmdz
     11      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 
    1112      IMPLICIT NONE
    1213c
    … …  
    2021c   Arguments:
    2122c   ----------
    22       REAL masse(ijb_u:ije_u,llm),pente_max
     23      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
    2324      REAL u_m( ijb_u:ije_u,llm )
    24       REAL q(ijb_u:ije_u,llm)
     25      REAL q(ijb_u:ije_u,llm,nqtot)
    2526      REAL qsat(ijb_u:ije_u,llm)
     27      INTEGER iq ! CRisi
    2628c
    2729c      Local
    … …  
    3638      REAL adxqu(ijb_u:ije_u),dxqmax(ijb_u:ije_u,llm)
    3739      REAL u_mq(ijb_u:ije_u,llm)
     40      REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
     41      INTEGER ifils,iq2 ! CRisi
     42
    3843
    3944      REAL      SSUM
    … …  
    4247      INTEGER ijb,ije,ijb_x,ije_x
    4348     
     49      write(*,*) 'vlspltqs 58: entree vlxqs_loc, iq,ijb_x=',
     50     &   iq,ijb_x
    4451
    4552c   calcul de la pente a droite et a gauche de la maille
    … …  
    6572         DO l = 1, llm
    6673            DO ij=ijb,ije-1
    67                dxqu(ij)=q(ij+1,l)-q(ij,l)
     74               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    6875c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
    69 c              sigu(ij)=u_m(ij,l)/masse(ij,l)
     76c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
    7077            ENDDO
    7178            DO ij=ijb+iip1-1,ije,iip1
    … …  
    120127         DO l = 1, llm
    121128            DO ij=ijb,ije-1
    122                dxqu(ij)=q(ij+1,l)-q(ij,l)
     129               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    123130            ENDDO
    124131            DO ij=ijb+iip1-1,ije,iip1
    … …  
    179186      DO l=1,llm
    180187       DO ij=ijb,ije-1
    181           zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
    182      ,                     1.+u_m(ij,l)/masse(ij+1,l),
     188          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
     189     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
    183190     ,                     u_m(ij,l))
    184191          zdum(ij,l)=0.5*zdum(ij,l)
    185192          u_mq(ij,l)=cvmgp(
    186      ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
    187      ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
     193     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
     194     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
    188195     ,                u_m(ij,l))
    189196          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
    … …  
    195202c   on cumule le flux correspondant a toutes les mailles dont la masse
    196203c   au travers de la paroi pENDant le pas de temps.
    197 c   le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
     204c   le rapport de melange de l''air advecte est min(q_vanleer, Qsat_downwind)
    198205c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    199206      DO l=1,llm
    200207       DO ij=ijb,ije-1
    201208          IF (u_m(ij,l).gt.0.) THEN
    202              zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
     209             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
    203210             u_mq(ij,l)=u_m(ij,l)*
    204      $         min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
     211     $         min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
    205212          ELSE
    206              zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
     213             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
    207214             u_mq(ij,l)=u_m(ij,l)*
    208      $         min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
     215     $         min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
    209216          ENDIF
    210217       ENDDO
    … …  
    273280               ENDDO
    274281               niju=iju
    275 c              PRINT*,'niju,nl',niju,nl(l)
     282               PRINT*,'vlxqs 280: niju,nl',niju,nl(l)
    276283
    277284c  traitement des mailles
    … …  
    285292                     i=ijq-(j-1)*iip1
    286293c   accumulation pour les mailles completements advectees
    287                      do while(zu_m.gt.masse(ijq,l))
    288                         u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
    289                         zu_m=zu_m-masse(ijq,l)
     294                     do while(zu_m.gt.masse(ijq,l,iq))
     295                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
     296     &                     *masse(ijq,l,iq)
     297                        zu_m=zu_m-masse(ijq,l,iq)
    290298                        i=mod(i-2+iim,iim)+1
    291299                        ijq=(j-1)*iip1+i
    292300                     ENDDO
    293301c   ajout de la maille non completement advectee
    294                      u_mq(ij,l)=u_mq(ij,l)+zu_m*
    295      &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
     302                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)
     303     &                 +0.5*(1.-zu_m/masse(ijq,l,iq))*dxq(ijq,l))
    296304                  ELSE
    297305                     ijq=ij+1
    298306                     i=ijq-(j-1)*iip1
    299307c   accumulation pour les mailles completements advectees
    300                      do while(-zu_m.gt.masse(ijq,l))
    301                         u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
    302                         zu_m=zu_m+masse(ijq,l)
     308                     do while(-zu_m.gt.masse(ijq,l,iq))
     309                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
     310     &                   *masse(ijq,l,iq)
     311                        zu_m=zu_m+masse(ijq,l,iq)
    303312                        i=mod(i,iim)+1
    304313                        ijq=(j-1)*iip1+i
    305314                     ENDDO
    306315c   ajout de la maille non completement advectee
    307                      u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
    308      &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
     316                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
     317     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
    309318                  ENDIF
    310319               ENDDO
    … …  
    325334c$OMP END DO NOWAIT
    326335
     336! CRisi: appel récursif de l'advection sur les fils.
     337! Il faut faire ça avant d'avoir mis à jour q et masse
     338      write(*,*) 'vlspltqs 336: iq,ijb_x,nqfils(iq)=',
     339     &     iq,ijb_x,nqfils(iq) 
     340
     341      if (nqfils(iq).gt.0) then 
     342       do ifils=1,nqdesc(iq)
     343         iq2=iqfils(ifils,iq)
     344c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
     345         DO l=1,llm
     346          DO ij=ijb,ije
     347           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     348           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
     349          enddo   
     350         enddo
     351c$OMP END DO NOWAIT
     352        enddo !do ifils=1,nqfils(iq)
     353        do ifils=1,nqfils(iq)
     354         iq2=iqfils(ifils,iq)
     355         write(*,*) 'vlxqs 349: on appelle vlx pour iq2=',iq2
     356         call vlx_loc(Ratio,pente_max,masseq,u_mq,ijb_x,ije_x,iq2)
     357        enddo !do ifils=1,nqfils(iq)
     358      endif !if (nqfils(iq).gt.0) then
     359! end CRisi
     360
     361      write(*,*) 'vlspltqs 360: iq,ijb_x=',iq,ijb_x   
     362
    327363c   calcul des tendances
    328364c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    329365      DO l=1,llm
    330366         DO ij=ijb+1,ije
    331             new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
    332             q(ij,l)=(q(ij,l)*masse(ij,l)+
     367            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
     368            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
    333369     &      u_mq(ij-1,l)-u_mq(ij,l))
    334370     &      /new_m
    335             masse(ij,l)=new_m
    336          ENDDO
    337 c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
     371            masse(ij,l,iq)=new_m
     372         ENDDO
     373c   Modif Fred 22 03 96 correction d''un bug (les scopy ci-dessous)
    338374         DO ij=ijb+iip1-1,ije,iip1
    339             q(ij-iim,l)=q(ij,l)
    340             masse(ij-iim,l)=masse(ij,l)
    341          ENDDO
    342       ENDDO
    343 c$OMP END DO NOWAIT
     375            q(ij-iim,l,iq)=q(ij,l,iq)
     376            masse(ij-iim,l,iq)=masse(ij,l,iq)
     377         ENDDO
     378      ENDDO
     379c$OMP END DO NOWAIT
     380
     381      write(*,*) 'vlspltqs 380: iq,ijb_x=',iq,ijb_x
     382
     383! retablir les fils en rapport de melange par rapport a l'air:
     384      if (nqfils(iq).gt.0) then 
     385       do ifils=1,nqdesc(iq)
     386         iq2=iqfils(ifils,iq) 
     387c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
     388         DO l=1,llm
     389          DO ij=ijb+1,ije
     390            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     391          enddo
     392          DO ij=ijb+iip1-1,ije,iip1
     393             q(ij-iim,l,iq2)=q(ij,l,iq2)
     394          enddo ! DO ij=ijb+iip1-1,ije,iip1
     395         enddo
     396c$OMP END DO NOWAIT
     397        enddo !do ifils=1,nqdesc(iq)
     398      endif !if (nqfils(iq).gt.0) then
     399
     400      write(*,*) 'vlspltqs 399: iq,ijb_x=',iq,ijb_x
     401
    344402c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
    345 c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
     403c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1,iq),iip1,masse(iip2,1,iq),iip1)
    346404
    347405
    348406      RETURN
    349407      END
    350       SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat)
     408      SUBROUTINE vlyqs_loc(q,pente_max,masse,masse_adv_v,qsat,iq)
    351409c
    352410c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    361419c   --------------------------------------------------------------------
    362420      USE parallel_lmdz
     421      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 
    363422      IMPLICIT NONE
    364423c
    … …  
    373432c   Arguments:
    374433c   ----------
    375       REAL masse(ijb_u:ije_u,llm),pente_max
     434      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
    376435      REAL masse_adv_v( ijb_v:ije_v,llm)
    377       REAL q(ijb_u:ije_u,llm)
     436      REAL q(ijb_u:ije_u,llm,nqtot)
    378437      REAL qsat(ijb_u:ije_u,llm)
     438      INTEGER iq ! CRisi
    379439c
    380440c      Local
    … …  
    386446      REAL dyq(ijb_u:ije_u,llm),dyqv(ijb_v:ije_v)
    387447      REAL adyqv(ijb_v:ije_v),dyqmax(ijb_u:ije_u)
    388       REAL qbyv(ijb_v:ije_v,llm)
     448      REAL qbyv(ijb_v:ije_v,llm,nqtot)
    389449
    390450      REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
    … …  
    402462c
    403463c
     464      REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
     465      INTEGER ifils,iq2 ! CRisi
     466
    404467      REAL      SSUM
    405468
    … …  
    407470      INTEGER ijb,ije
    408471
     472      ijb=ij_begin-2*iip1
     473      ije=ij_end+2*iip1 
     474      if (pole_nord) ijb=ij_begin
     475      if (pole_sud)  ije=ij_end
     476      ij=3525
     477      l=3
     478      if ((ij.ge.ijb).and.(ij.le.ije)) then
     479        write(*,*) 'vlyqs 480: ij,l,iq,ijb,q(ij,l,:)=',
     480     &             ij,l,iq,ijb,q(ij,l,:)
     481      endif 
     482
    409483      IF(first) THEN
    410484         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
     485         PRINT*,'vlyqs_loc, iq=',iq
    411486         first=.false.
    412487         do i=2,iip1
    … …  
    439514      if (pole_nord) then
    440515        DO i = 1, iim
    441           airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
     516          airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
    442517        ENDDO
    443518        qpns   = SSUM( iim,  airescb ,1 ) / airej2
    … …  
    446521      if (pole_sud) then
    447522        DO i = 1, iim
    448           airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
     523          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
    449524        ENDDO
    450525        qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
    … …  
    460535     
    461536      DO ij=ijb,ije
    462          dyqv(ij)=q(ij,l)-q(ij+iip1,l)
     537         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
    463538         adyqv(ij)=abs(dyqv(ij))
    464539      ENDDO
    … …  
    482557c   calcul des pentes aux poles
    483558        DO ij=1,iip1
    484            dyq(ij,l)=qpns-q(ij+iip1,l)
     559           dyq(ij,l)=qpns-q(ij+iip1,l,iq)
    485560        ENDDO
    486561
    … …  
    513588
    514589        DO ij=1,iip1
    515            dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
     590           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
    516591        ENDDO
    517592
    … …  
    636711       DO ij=ijb,ije
    637712         IF( masse_adv_v(ij,l).GT.0. ) THEN
    638            qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l )  +
    639      ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
     713           qbyv(ij,l,iq)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq )  +
     714     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)
     715     ,      /masse(ij+iip1,l,iq)))
    640716         ELSE
    641               qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
    642      ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
     717              qbyv(ij,l,iq)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) *
     718     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) )
    643719         ENDIF
    644           qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
     720          qbyv(ij,l,iq) = masse_adv_v(ij,l)*qbyv(ij,l,iq)
    645721       ENDDO
    646722      ENDDO
    647723c$OMP END DO NOWAIT
     724
     725! CRisi: appel récursif de l'advection sur les fils.
     726! Il faut faire ça avant d'avoir mis à jour q et masse
     727      write(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq)
     728     
     729      ijb=ij_begin-2*iip1
     730      ije=ij_end+2*iip1
     731      if (pole_nord) ijb=ij_begin
     732      if (pole_sud)  ije=ij_end 
     733
     734      if (nqfils(iq).gt.0) then 
     735       do ifils=1,nqdesc(iq)
     736         iq2=iqfils(ifils,iq)
     737c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
     738         DO l=1,llm
     739          DO ij=ijb,ije
     740           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     741           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)   
     742          enddo   
     743         enddo
     744c$OMP END DO NOWAIT
     745        enddo !do ifils=1,nqdesc(iq)
     746        do ifils=1,nqfils(iq)
     747         iq2=iqfils(ifils,iq)
     748         call vly_loc(Ratio,pente_max,masseq,qbyv,iq2)
     749        enddo !do ifils=1,nqfils(iq)
     750      endif !if (nqfils(iq).gt.0) then
     751
     752       
     753! end CRisi
    648754
    649755      ijb=ij_begin
    … …  
    655761      DO l=1,llm
    656762         DO ij=ijb,ije
    657             newmasse=masse(ij,l)
     763            newmasse=masse(ij,l,iq)
    658764     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
    659             q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
    660      &         /newmasse
    661             masse(ij,l)=newmasse
     765            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l,iq)
     766     &         -qbyv(ij-iip1,l,iq))/newmasse
     767            masse(ij,l,iq)=newmasse
    662768         ENDDO
    663769c.-. ancienne version
    … …  
    665771         IF (pole_nord) THEN
    666772
    667            convpn=SSUM(iim,qbyv(1,l),1)/apoln
     773           convpn=SSUM(iim,qbyv(1,l,iq),1)/apoln
    668774           convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
    669775           DO ij = 1,iip1
    670               newmasse=masse(ij,l)+convmpn*aire(ij)
    671               q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
     776              newmasse=masse(ij,l,iq)+convmpn*aire(ij)
     777              q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convpn*aire(ij))/
    672778     &                 newmasse
    673               masse(ij,l)=newmasse
     779              masse(ij,l,iq)=newmasse
    674780           ENDDO
    675781         
    … …  
    678784         IF (pole_sud) THEN
    679785         
    680            convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
     786           convps  = -SSUM(iim,qbyv(ip1jm-iim,l,iq),iq,1)/apols
    681787           convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
    682788           DO ij = ip1jm+1,ip1jmp1
    683               newmasse=masse(ij,l)+convmps*aire(ij)
    684               q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
     789              newmasse=masse(ij,l,iq)+convmps*aire(ij)
     790              q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convps*aire(ij))/
    685791     &                 newmasse
    686               masse(ij,l)=newmasse
     792              masse(ij,l,iq)=newmasse
    687793           ENDDO
    688794         
    … …  
    691797
    692798c._. nouvelle version
    693 c        convpn=SSUM(iim,qbyv(1,l),1)
     799c        convpn=SSUM(iim,qbyv(1,l,iq),1)
    694800c        convmpn=ssum(iim,masse_adv_v(1,l),1)
    695 c        oldmasse=ssum(iim,masse(1,l),1)
     801c        oldmasse=ssum(iim,masse(1,l,iq),1)
    696802c        newmasse=oldmasse+convmpn
    697 c        newq=(q(1,l)*oldmasse+convpn)/newmasse
     803c        newq=(q(1,l,iq)*oldmasse+convpn)/newmasse
    698804c        newmasse=newmasse/apoln
    699805c        DO ij = 1,iip1
    700 c           q(ij,l)=newq
    701 c           masse(ij,l)=newmasse*aire(ij)
     806c           q(ij,l,iq)=newq
     807c           masse(ij,l,iq)=newmasse*aire(ij)
    702808c        ENDDO
    703 c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
     809c        convps=-SSUM(iim,qbyv(ip1jm-iim,l,iq),1)
    704810c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
    705 c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
     811c        oldmasse=ssum(iim,masse(ip1jm-iim,l,iq),1)
    706812c        newmasse=oldmasse+convmps
    707 c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
     813c        newq=(q(ip1jmp1,l,iq)*oldmasse+convps)/newmasse
    708814c        newmasse=newmasse/apols
    709815c        DO ij = ip1jm+1,ip1jmp1
    710 c           q(ij,l)=newq
    711 c           masse(ij,l)=newmasse*aire(ij)
     816c           q(ij,l,iq)=newq
     817c           masse(ij,l,iq)=newmasse*aire(ij)
    712818c        ENDDO
    713819c._. fin nouvelle version
    714820      ENDDO
    715821c$OMP END DO NOWAIT
     822
     823! retablir les fils en rapport de melange par rapport a l'air:
     824      ijb=ij_begin
     825      ije=ij_end
     826!      if (pole_nord) ijb=ij_begin+iip1
     827!      if (pole_sud)  ije=ij_end-iip1
     828 
     829      if (nqfils(iq).gt.0) then 
     830       do ifils=1,nqdesc(iq)
     831         iq2=iqfils(ifils,iq) 
     832c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
     833         DO l=1,llm
     834          DO ij=ijb,ije
     835            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     836          enddo
     837         enddo
     838c$OMP END DO NOWAIT
     839        enddo !do ifils=1,nqdesc(iq)
     840      endif !if (nqfils(iq).gt.0) then
     841
     842
    716843      RETURN
    717844      END
Note: See TracChangeset for help on using the changeset viewer.