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/dyn3dmem/vlsplt_loc.F

    r1910 r2298  
    22! $Id$
    33!
    4       SUBROUTINE vlx_loc(q,pente_max,masse,u_m,ijb_x,ije_x)
     4      RECURSIVE SUBROUTINE vlx_loc(q,pente_max,masse,u_m,ijb_x,ije_x,iq)
    55
    66c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    1414c   --------------------------------------------------------------------
    1515      USE parallel_lmdz
     16      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    1617      IMPLICIT NONE
    1718c
    … …  
    2526c   Arguments:
    2627c   ----------
    27       REAL masse(ijb_u:ije_u,llm),pente_max
    28       REAL u_m( ijb_u:ije_u,llm ),pbarv( iip1,jjb_v:jje_v,llm)
    29       REAL q(ijb_u:ije_u,llm)
    30       REAL w(ijb_u:ije_u,llm)
     28      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
     29      REAL u_m( ijb_u:ije_u,llm),pbarv( iip1,jjb_v:jje_v,llm)
     30      REAL q(ijb_u:ije_u,llm,nqtot) ! CRisi: ajout dimension nqtot
     31      REAL w(ijb_u:ije_u,llm)
     32      INTEGER iq ! CRisi
    3133c
    3234c      Local
    … …  
    4244      REAL u_mq(ijb_u:ije_u,llm)
    4345
     46      REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
     47      INTEGER ifils,iq2 ! CRisi
     48
    4449      Logical extremum
    4550
    … …  
    5156      INTEGER ijb,ije,ijb_x,ije_x
    5257     
     58      !write(*,*) 'vlsplt 58: entree dans vlx_loc, iq,ijb_x=',
     59!     &   iq,ijb_x
    5360c   calcul de la pente a droite et a gauche de la maille
    5461
    … …  
    6471c   calcul des pentes avec limitation, Van Leer scheme I:
    6572c   -----------------------------------------------------
    66 
     73      ! on a besoin de q entre ijb et ije
    6774c   calcul de la pente aux points u
    6875c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
    … …  
    7077           
    7178            DO ij=ijb,ije-1
    72                dxqu(ij)=q(ij+1,l)-q(ij,l)
     79               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    7380c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
    74 c              sigu(ij)=u_m(ij,l)/masse(ij,l)
     81c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
    7582            ENDDO
    7683            DO ij=ijb+iip1-1,ije,iip1
    … …  
    126133         DO l = 1, llm
    127134            DO ij=ijb,ije-1
    128                dxqu(ij)=q(ij+1,l)-q(ij,l)
     135               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    129136            ENDDO
    130137            DO ij=ijb+iip1-1,ije,iip1
    … …  
    147154      ENDIF ! (pente_max.lt.-1.e-5)
    148155
     156      !write(*,*) 'vlx 156: iq,ijb_x=',iq,ijb_x
     157
    149158c   bouclage de la pente en iip1:
    150159c   -----------------------------
    … …  
    168177      DO l=1,llm
    169178       DO ij=ijb,ije-1
    170           zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
    171      ,                     1.+u_m(ij,l)/masse(ij+1,l),
    172      ,                     u_m(ij,l))
     179          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
     180     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
     181     ,                     u_m(ij,l,iq))
    173182          zdum(ij,l)=0.5*zdum(ij,l)
    174183          u_mq(ij,l)=cvmgp(
    175      ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
    176      ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
     184     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
     185     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
    177186     ,                u_m(ij,l))
    178187          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
    … …  
    185194c       print*,'Cumule ....'
    186195c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     196        ! on a besoin de masse entre ijb et ije
    187197      DO l=1,llm
    188198       DO ij=ijb,ije-1
    189 c       print*,'masse(',ij,')=',masse(ij,l)
     199c       print*,'masse(',ij,')=',masse(ij,l,iq)
    190200          IF (u_m(ij,l).gt.0.) THEN
    191              zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
    192              u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
     201             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
     202             u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq)
     203     :           +0.5*zdum(ij,l)*dxq(ij,l))
    193204          ELSE
    194              zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
    195              u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
     205             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
     206             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq)
     207     :           -0.5*zdum(ij,l)*dxq(ij+1,l))
    196208          ENDIF
    197209       ENDDO
    … …  
    215227c$OMP END DO NOWAIT
    216228c       print*,'Ok test 1'
     229
    217230c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    218231      DO l=1,llm
    … …  
    223236c$OMP END DO NOWAIT
    224237c        print*,'Ok test 2'
    225 
     238       
    226239
    227240c   traitement special pour le cas ou on advecte en longitude plus que le
    … …  
    247260c     &       ,'contenu de la maille : ',n0
    248261c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     262
     263
    249264         DO l=1,llm
    250265            IF(nl(l).gt.0) THEN
    … …  
    258273               ENDDO
    259274               niju=iju
    260 c              PRINT*,'niju,nl',niju,nl(l)
     275               !PRINT*,'vlx 278, niju,nl',niju,nl(l)
    261276
    262277c  traitement des mailles
    … …  
    270285                     i=ijq-(j-1)*iip1
    271286c   accumulation pour les mailles completements advectees
    272                      do while(zu_m.gt.masse(ijq,l))
    273                         u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
    274                         zu_m=zu_m-masse(ijq,l)
     287                     do while(zu_m.gt.masse(ijq,l,iq))
     288                        u_mq(ij,l)=u_mq(ij,l)
     289     &                          +q(ijq,l,iq)*masse(ijq,l,iq)
     290                        zu_m=zu_m-masse(ijq,l,iq)
    275291                        i=mod(i-2+iim,iim)+1
    276292                        ijq=(j-1)*iip1+i
    … …  
    278294c   ajout de la maille non completement advectee
    279295                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
    280      &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
     296     &               (q(ijq,l,iq)+0.5*
     297     &               (1.-zu_m/masse(ijq,l,iq))*dxq(ijq,l))
    281298                  ELSE
    282299                     ijq=ij+1
    283300                     i=ijq-(j-1)*iip1
    284301c   accumulation pour les mailles completements advectees
    285                      do while(-zu_m.gt.masse(ijq,l))
    286                         u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
    287                         zu_m=zu_m+masse(ijq,l)
     302                     do while(-zu_m.gt.masse(ijq,l,iq))
     303                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
     304     &                           *masse(ijq,l,iq)
     305                        zu_m=zu_m+masse(ijq,l,iq)
    288306                        i=mod(i,iim)+1
    289307                        ijq=(j-1)*iip1+i
    290308                     ENDDO
    291309c   ajout de la maille non completement advectee
    292                      u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
    293      &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
     310                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
     311     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
    294312                  ENDIF
    295313               ENDDO
    … …  
    299317cym      ENDIF  ! n0.gt.0
    3003189999    continue
    301 
    302319
    303320c   bouclage en latitude
    … …  
    311328c$OMP END DO NOWAIT
    312329
     330! CRisi: appel récursif de l'advection sur les fils.
     331! Il faut faire ça avant d'avoir mis à jour q et masse
     332
     333      !write(*,*) 'vlsplt 326: iq,ijb_x,nqfils(iq)=',iq,ijb_x,nqfils(iq)
     334
     335      if (nqfils(iq).gt.0) then 
     336       do ifils=1,nqdesc(iq)
     337         iq2=iqfils(ifils,iq)
     338c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     339         DO l=1,llm
     340          DO ij=ijb,ije
     341           ! On a besoin de q et masse seulement entre ijb et ije. On ne
     342           ! les calcule donc que de ijb à ije
     343           masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     344           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     345          enddo   
     346         enddo
     347c$OMP END DO NOWAIT
     348        enddo !do ifils=1,nqdesc(iq)
     349        do ifils=1,nqfils(iq)
     350         iq2=iqfils(ifils,iq)
     351         call vlx_loc(Ratio,pente_max,masse,u_mq,ijb_x,ije_x,iq2)
     352        enddo !do ifils=1,nqfils(iq)
     353      endif !if (nqfils(iq).gt.0) then
     354! end CRisi
     355
     356      !write(*,*) 'vlsplt 360: iq,ijb_x=',iq,ijb_x
     357
    313358c   calcul des tENDances
    314359c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    315360      DO l=1,llm
    316361         DO ij=ijb+1,ije
    317             new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
    318             q(ij,l)=(q(ij,l)*masse(ij,l)+
    319      &      u_mq(ij-1,l)-u_mq(ij,l))
    320      &      /new_m
    321             masse(ij,l)=new_m
     362            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
     363            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
     364     &        u_mq(ij-1,l)-u_mq(ij,l))
     365     &        /new_m
     366            masse(ij,l,iq)=new_m
    322367         ENDDO
    323368c   ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
    324369         DO ij=ijb+iip1-1,ije,iip1
    325             q(ij-iim,l)=q(ij,l)
    326             masse(ij-iim,l)=masse(ij,l)
    327          ENDDO
    328       ENDDO
    329 c$OMP END DO NOWAIT
     370            q(ij-iim,l,iq)=q(ij,l,iq)
     371            masse(ij-iim,l,iq)=masse(ij,l,iq)
     372         ENDDO
     373      ENDDO
     374c$OMP END DO NOWAIT
     375      !write(*,*) 'vlsplt 380: iq,ijb_x=',iq,ijb_x
     376
     377! retablir les fils en rapport de melange par rapport a l'air:
     378      ! On calcule q entre ijb+1 et ije -> on fait pareil pour ratio
     379      ! puis on boucle en longitude
     380      if (nqfils(iq).gt.0) then 
     381       do ifils=1,nqdesc(iq)
     382         iq2=iqfils(ifils,iq) 
     383c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
     384         DO l=1,llm
     385          DO ij=ijb+1,ije
     386            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     387          enddo
     388          DO ij=ijb+iip1-1,ije,iip1
     389             q(ij-iim,l,iq2)=q(ij,l,iq2)
     390          enddo ! DO ij=ijb+iip1-1,ije,iip1
     391         enddo !DO l=1,llm
     392c$OMP END DO NOWAIT
     393        enddo !do ifils=1,nqdesc(iq)
     394      endif !if (nqfils(iq).gt.0) then
     395
     396      !write(*,*) 'vlsplt 399: iq,ijb_x=',iq,ijb_x
    330397c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
    331398c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
    … …  
    336403
    337404
    338       SUBROUTINE vly_loc(q,pente_max,masse,masse_adv_v)
     405      RECURSIVE SUBROUTINE vly_loc(q,pente_max,masse,masse_adv_v,iq)
    339406c
    340407c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    349416c   --------------------------------------------------------------------
    350417      USE parallel_lmdz
     418      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    351419      IMPLICIT NONE
    352420c
    … …  
    361429c   Arguments:
    362430c   ----------
    363       REAL masse(ijb_u:ije_u,llm),pente_max
     431      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
    364432      REAL masse_adv_v( ijb_v:ije_v,llm)
    365       REAL q(ijb_u:ije_u,llm), dq( ijb_u:ije_u,llm)
     433      REAL q(ijb_u:ije_u,llm,nqtot), dq( ijb_u:ije_u,llm)
     434      INTEGER iq ! CRisi
    366435c
    367436c      Local
    … …  
    392461      SAVE airej2,airejjm
    393462c$OMP THREADPRIVATE(airej2,airejjm)
     463
     464      REAL Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
     465      INTEGER ifils,iq2 ! CRisi
    394466c
    395467c
    … …  
    401473      INTEGER ijb,ije
    402474
     475      ijb=ij_begin-2*iip1
     476      ije=ij_end+2*iip1 
     477      if (pole_nord) ijb=ij_begin
     478      if (pole_sud)  ije=ij_end
     479
    403480      IF(first) THEN
    404 c         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
     481         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
    405482         first=.false.
    406483         do i=2,iip1
    … …  
    434511      if (pole_nord) then
    435512        DO i = 1, iim
    436           airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
     513          airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
    437514        ENDDO
    438515        qpns   = SSUM( iim,  airescb ,1 ) / airej2
    … …  
    441518      if (pole_sud) then
    442519        DO i = 1, iim
    443           airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
     520          airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
    444521        ENDDO
    445522        qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
    446523      endif
    447524     
    448      
    449 
    450525c   calcul des pentes aux points v
    451526
    … …  
    455530      if (pole_sud)  ije=ij_end-iip1
    456531     
     532      ! on a besoin de q entre ij_begin-2*iip1 et ij_end+2*iip1
     533      ! Si pole sud, entre ij_begin-2*iip1 et ij_end
     534      ! Si pole Nord, entre ij_begin et ij_end+2*iip1
    457535      DO ij=ijb,ije
    458          dyqv(ij)=q(ij,l)-q(ij+iip1,l)
     536         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
    459537         adyqv(ij)=abs(dyqv(ij))
    460538      ENDDO
     539 
    461540
    462541c   calcul des pentes aux points scalaires
    … …  
    475554      IF (pole_nord) THEN
    476555        DO ij=1,iip1
    477            dyq(ij,l)=qpns-q(ij+iip1,l)
     556           dyq(ij,l)=qpns-q(ij+iip1,l,iq)
    478557        ENDDO
    479558       
    … …  
    497576
    498577        DO ij=1,iip1
    499            dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
     578           dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
    500579        ENDDO
    501580
    … …  
    633712       DO ij=ijb,ije
    634713          IF(masse_adv_v(ij,l).gt.0) THEN
    635               qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
    636      ,                   0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
     714              qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)*
     715     ,                   0.5*(1.-masse_adv_v(ij,l)
     716     ,                   /masse(ij+iip1,l,iq))
    637717          ELSE
    638               qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
    639      ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
     718              qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)*
     719     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq))
    640720          ENDIF
    641721          qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
    … …  
    643723      ENDDO
    644724c$OMP END DO NOWAIT
     725
     726! CRisi: appel récursif de l'advection sur les fils.
     727! Il faut faire ça avant d'avoir mis à jour q et masse
     728      !write(*,*) 'vly 689: iq,nqfils(iq)=',iq,nqfils(iq)
     729
     730      ijb=ij_begin-2*iip1
     731      ije=ij_end+2*iip1
     732      if (pole_nord) ijb=ij_begin
     733      if (pole_sud)  ije=ij_end
     734   
     735      if (nqfils(iq).gt.0) then 
     736       do ifils=1,nqdesc(iq)
     737         iq2=iqfils(ifils,iq)
     738c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     739         DO l=1,llm
     740         DO ij=ijb,ije
     741           masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     742           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
     743          enddo   
     744         enddo
     745c$OMP END DO NOWAIT
     746        enddo !do ifils=1,nqdesc(iq)
     747
     748        do ifils=1,nqfils(iq)
     749         iq2=iqfils(ifils,iq)
     750         call vly_loc(Ratio,pente_max,masse,qbyv,iq2)
     751        enddo !do ifils=1,nqfils(iq)
     752      endif !if (nqfils(iq).gt.0) then
     753! end CRisi
    645754     
    646755      ijb=ij_begin
    … …  
    652761      DO l=1,llm
    653762         DO ij=ijb,ije
    654             newmasse=masse(ij,l)
    655      &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
    656      
    657             q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
    658      &         /newmasse
    659             masse(ij,l)=newmasse
    660          ENDDO
     763            newmasse=masse(ij,l,iq)
     764     &         +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
     765
     766            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
     767     &         -qbyv(ij-iip1,l))/newmasse
     768
     769            masse(ij,l,iq)=newmasse
     770
     771         ENDDO
     772
     773
    661774c.-. ancienne version
    662775c        convpn=SSUM(iim,qbyv(1,l),1)/apoln
    … …  
    665778           convpn=SSUM(iim,qbyv(1,l),1)
    666779           convmpn=ssum(iim,masse_adv_v(1,l),1)
    667            massepn=ssum(iim,masse(1,l),1)
     780           massepn=ssum(iim,masse(1,l,iq),1)
    668781           qpn=0.
    669782           do ij=1,iim
    670               qpn=qpn+masse(ij,l)*q(ij,l)
     783              qpn=qpn+masse(ij,l,iq)*q(ij,l,iq)
    671784           enddo
    672785           qpn=(qpn+convpn)/(massepn+convmpn)
    673786           do ij=1,iip1
    674               q(ij,l)=qpn
     787              q(ij,l,iq)=qpn
    675788           enddo
    676789         endif
    … …  
    683796           convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    684797           convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
    685            masseps=ssum(iim, masse(ip1jm+1,l),1)
     798           masseps=ssum(iim, masse(ip1jm+1,l,iq),1)
    686799           qps=0.
    687800           do ij = ip1jm+1,ip1jmp1-1
    688               qps=qps+masse(ij,l)*q(ij,l)
     801              qps=qps+masse(ij,l,iq)*q(ij,l,iq)
    689802           enddo
    690803           qps=(qps+convps)/(masseps+convmps)
    691804           do ij=ip1jm+1,ip1jmp1
    692               q(ij,l)=qps
     805              q(ij,l,iq)=qps
    693806           enddo
    694807         endif
    … …  
    704817c        DO ij = 1,iip1
    705818c           q(ij,l)=newq
    706 c           masse(ij,l)=newmasse*aire(ij)
     819c           masse(ij,l,iq)=newmasse*aire(ij)
    707820c        ENDDO
    708821c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    … …  
    714827c        DO ij = ip1jm+1,ip1jmp1
    715828c           q(ij,l)=newq
    716 c           masse(ij,l)=newmasse*aire(ij)
     829c           masse(ij,l,iq)=newmasse*aire(ij)
    717830c        ENDDO
    718831c._. fin nouvelle version
    … …  
    720833c$OMP END DO NOWAIT
    721834
     835! retablir les fils en rapport de melange par rapport a l'air:
     836      ijb=ij_begin
     837      ije=ij_end
     838!      if (pole_nord) ijb=ij_begin
     839!      if (pole_sud)  ije=ij_end
     840
     841      if (nqfils(iq).gt.0) then 
     842       do ifils=1,nqdesc(iq)
     843         iq2=iqfils(ifils,iq) 
     844c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
     845         DO l=1,llm
     846          DO ij=ijb,ije
     847            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     848          enddo
     849         enddo
     850c$OMP END DO NOWAIT
     851        enddo !do ifils=1,nqdesc(iq)
     852      endif !if (nqfils(iq).gt.0) then
     853
     854
    722855      RETURN
    723856      END
    … …  
    725858     
    726859     
    727       SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x)
     860      RECURSIVE SUBROUTINE vlz_loc(q,pente_max,masse,w,ijb_x,ije_x,iq)
    728861c
    729862c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    739872      USE parallel_lmdz
    740873      USE vlz_mod
     874      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi 
    741875      IMPLICIT NONE
    742876c
    … …  
    750884c   Arguments:
    751885c   ----------
    752       REAL masse(ijb_u:ije_u,llm),pente_max
    753       REAL q(ijb_u:ije_u,llm)
    754       REAL w(ijb_u:ije_u,llm+1)
     886      REAL masse(ijb_u:ije_u,llm,nqtot),pente_max
     887      REAL q(ijb_u:ije_u,llm,nqtot)
     888      REAL w(ijb_u:ije_u,llm+1,nqtot)
     889      INTEGER iq
    755890c
    756891c      Local
    … …  
    779914      LOGICAL,SAVE :: first=.TRUE.
    780915!$OMP THREADPRIVATE(first)
    781      
     916
     917      !REAL masseq(ijb_u:ije_u,llm,nqtot),Ratio(ijb_u:ije_u,llm,nqtot) ! CRisi
     918      ! Ces varibles doivent être déclarées en pointer et en save dans
     919      ! vlz_loc si on veut qu'elles soient vues par tous les threads. 
     920      INTEGER ifils,iq2 ! CRisi
    782921
    783922      IF (first) THEN
    … …  
    787926c    sens de W
    788927
     928      !write(*,*) 'vlsplt 926: entree dans vlz_loc, iq=',iq
    789929#ifdef BIDON
    790930      IF(testcpu) THEN
    … …  
    799939      DO l=2,llm
    800940         DO ij=ijb,ije
    801             dzqw(ij,l)=q(ij,l-1)-q(ij,l)
     941            dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq)
    802942            adzqw(ij,l)=abs(dzqw(ij,l))
    803943         ENDDO
    … …  
    842982c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
    843983
     984       !write(*,*) 'vlz 982,ijb,ije=',ijb,ije
    844985c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    845986       DO l = 1,llm-1
    846987         do  ij = ijb,ije
    847           IF(w(ij,l+1).gt.0.) THEN
    848              sigw=w(ij,l+1)/masse(ij,l+1)
    849              wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
     988          IF(w(ij,l+1,iq).gt.0.) THEN
     989             sigw=w(ij,l+1,iq)/masse(ij,l+1,iq)
     990             wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l+1,iq)
     991     :           +0.5*(1.-sigw)*dzq(ij,l+1))
    850992          ELSE
    851              sigw=w(ij,l+1)/masse(ij,l)
    852              wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
     993             sigw=w(ij,l+1,iq)/masse(ij,l,iq)
     994             wq(ij,l+1,iq)=w(ij,l+1,iq)*(q(ij,l,iq)
     995     :           -0.5*(1.+sigw)*dzq(ij,l))
    853996          ENDIF
    854997         ENDDO
    855998       ENDDO
    856 c$OMP END DO NOWAIT
     999c$OMP END DO NOWAIT   
     1000       !write(*,*) 'vlz 1001'   
    8571001
    8581002c$OMP MASTER
    8591003       DO ij=ijb,ije
    860           wq(ij,llm+1)=0.
    861           wq(ij,1)=0.
     1004          wq(ij,llm+1,iq)=0.
     1005          wq(ij,1,iq)=0.
    8621006       ENDDO
    8631007c$OMP END MASTER
    8641008c$OMP BARRIER
    8651009
     1010! CRisi: appel récursif de l'advection sur les fils.
     1011! Il faut faire ça avant d'avoir mis à jour q et masse
     1012      !write(*,*) 'vlsplt 942: iq,nqfils(iq)=',iq,nqfils(iq)
     1013      if (nqfils(iq).gt.0) then 
     1014       do ifils=1,nqdesc(iq)
     1015         iq2=iqfils(ifils,iq)
     1016c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     1017         DO l=1,llm
     1018          DO ij=ijb,ije
     1019           masse(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     1020           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     1021           !wq(ij,l,iq2)=wq(ij,l,iq) ! correction bug le 15mai2015
     1022           w(ij,l,iq2)=wq(ij,l,iq)
     1023          enddo   
     1024         enddo
     1025c$OMP END DO NOWAIT
     1026        enddo !do ifils=1,nqdesc(iq)
     1027c$OMP BARRIER
     1028
     1029        do ifils=1,nqfils(iq)
     1030         iq2=iqfils(ifils,iq)
     1031         call vlz_loc(Ratio,pente_max,masse,w,ijb_x,ije_x,iq2)
     1032        enddo !do ifils=1,nqfils(iq)
     1033      endif !if (nqfils(iq).gt.0) then
     1034! end CRisi 
     1035
     1036! CRisi: On rajoute ici une barrière car on veut être sur que tous les
     1037! wq soient synchronisés
     1038
     1039c$OMP BARRIER
    8661040c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    8671041      DO l=1,llm
    8681042         DO ij=ijb,ije
    869             newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
    870             q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
     1043            newmasse=masse(ij,l,iq)+w(ij,l+1,iq)-w(ij,l,iq)
     1044            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)
     1045     &         +wq(ij,l+1,iq)-wq(ij,l,iq))
    8711046     &         /newmasse
    872             masse(ij,l)=newmasse
    873          ENDDO
    874       ENDDO
    875 c$OMP END DO NOWAIT
    876 
     1047            masse(ij,l,iq)=newmasse
     1048         ENDDO
     1049      ENDDO
     1050c$OMP END DO NOWAIT
     1051
     1052     
     1053! retablir les fils en rapport de melange par rapport a l'air:
     1054      if (nqfils(iq).gt.0) then 
     1055       do ifils=1,nqdesc(iq)
     1056         iq2=iqfils(ifils,iq) 
     1057c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
     1058         DO l=1,llm
     1059          DO ij=ijb,ije
     1060            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     1061          enddo
     1062         enddo
     1063c$OMP END DO NOWAIT
     1064        enddo !do ifils=1,nqdesc(iq)
     1065      endif !if (nqfils(iq).gt.0) then
    8771066
    8781067      RETURN
Note: See TracChangeset for help on using the changeset viewer.