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/dyn3d/vlspltqs.F

    r1907 r2270  
    33c
    44       SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
    5      ,                                  p,pk,teta                 )
     5     ,                                  p,pk,teta,iq             )
     6       USE infotrac, ONLY: nqtot,nqdesc,iqfils
    67c
    78c     Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron
    … …  
    3536      REAL masse(ip1jmp1,llm),pente_max
    3637      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
    37       REAL q(ip1jmp1,llm)
     38      REAL q(ip1jmp1,llm,nqtot)
    3839      REAL w(ip1jmp1,llm),pdt
    3940      REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
     41      INTEGER iq ! CRisi
    4042c
    4143c      Local
    … …  
    4345c
    4446      INTEGER i,ij,l,j,ii
     47      INTEGER ifils,iq2 ! CRisi
    4548c
    4649      REAL qsat(ip1jmp1,llm)
    47       REAL zm(ip1jmp1,llm)
     50      REAL zm(ip1jmp1,llm,nqtot)
    4851      REAL mu(ip1jmp1,llm)
    4952      REAL mv(ip1jm,llm)
    5053      REAL mw(ip1jmp1,llm+1)
    51       REAL zq(ip1jmp1,llm)
     54      REAL zq(ip1jmp1,llm,nqtot)
    5255      REAL temps1,temps2,temps3
    5356      REAL zzpbar, zzw
    … …  
    116119      ENDDO
    117120
    118       CALL SCOPY(ijp1llm,q,1,zq,1)
    119       CALL SCOPY(ijp1llm,masse,1,zm,1)
     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
    120129
    121130c      call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
    122       call vlxqs(zq,pente_max,zm,mu,qsat)
    123 
     131      call vlxqs(zq,pente_max,zm,mu,qsat,iq)
    124132
    125133c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
    126134
    127       call vlyqs(zq,pente_max,zm,mv,qsat)
    128 
     135      call vlyqs(zq,pente_max,zm,mv,qsat,iq)
    129136
    130137c      call minmaxq(zq,qmin,qmax,'avant vlz     ')
    131138
    132       call vlz(zq,pente_max,zm,mw)
    133 
     139      call vlz(zq,pente_max,zm,mw,iq)
    134140
    135141c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
    136142c     call minmaxq(zm,qmin,qmax,'M avant vlyqs     ')
    137143
    138       call vlyqs(zq,pente_max,zm,mv,qsat)
    139 
     144      call vlyqs(zq,pente_max,zm,mv,qsat,iq)
    140145
    141146c     call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
    142147c     call minmaxq(zm,qmin,qmax,'M avant vlxqs     ')
    143148
    144       call vlxqs(zq,pente_max,zm,mu,qsat)
     149      call vlxqs(zq,pente_max,zm,mu,qsat,iq)
    145150
    146151c     call minmaxq(zq,qmin,qmax,'apres vlxqs     ')
    … …  
    150155      DO l=1,llm
    151156         DO ij=1,ip1jmp1
    152            q(ij,l)=zq(ij,l)
     157           q(ij,l,iq)=zq(ij,l,iq)
    153158         ENDDO
    154159         DO ij=1,ip1jm+1,iip1
    155             q(ij+iim,l)=q(ij,l)
    156          ENDDO
    157       ENDDO
     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'
    158178
    159179      RETURN
    160180      END
    161       SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat)
     181      SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat,iq)
     182      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
     183
    162184c
    163185c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    179201c   Arguments:
    180202c   ----------
    181       REAL masse(ip1jmp1,llm),pente_max
     203      REAL masse(ip1jmp1,llm,nqtot),pente_max
    182204      REAL u_m( ip1jmp1,llm )
    183       REAL q(ip1jmp1,llm)
     205      REAL q(ip1jmp1,llm,nqtot)
    184206      REAL qsat(ip1jmp1,llm)
     207      INTEGER iq ! CRisi
    185208c
    186209c      Local
    … …  
    195218      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
    196219      REAL u_mq(ip1jmp1,llm)
     220
     221      ! CRisi
     222      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
     223      INTEGER ifils,iq2 ! CRisi
    197224
    198225      Logical first,testcpu
    … …  
    227254         DO l = 1, llm
    228255            DO ij=iip2,ip1jm-1
    229                dxqu(ij)=q(ij+1,l)-q(ij,l)
     256               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    230257c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
    231 c              sigu(ij)=u_m(ij,l)/masse(ij,l)
     258c              sigu(ij)=u_m(ij,l)/masse(ij,l,iq)
    232259            ENDDO
    233260            DO ij=iip1+iip1,ip1jm,iip1
    … …  
    281308         DO l = 1, llm
    282309            DO ij=iip2,ip1jm-1
    283                dxqu(ij)=q(ij+1,l)-q(ij,l)
     310               dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    284311            ENDDO
    285312            DO ij=iip1+iip1,ip1jm,iip1
    … …  
    323350      DO l=1,llm
    324351       DO ij=iip2,ip1jm-1
    325           zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
    326      ,                     1.+u_m(ij,l)/masse(ij+1,l),
     352          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l,iq),
     353     ,                     1.+u_m(ij,l)/masse(ij+1,l,iq),
    327354     ,                     u_m(ij,l))
    328355          zdum(ij,l)=0.5*zdum(ij,l)
    329356          u_mq(ij,l)=cvmgp(
    330      ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
    331      ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
     357     ,                q(ij,l,iq)+zdum(ij,l)*dxq(ij,l),
     358     ,                q(ij+1,l,iq)-zdum(ij,l)*dxq(ij+1,l),
    332359     ,                u_m(ij,l))
    333360          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
    … …  
    341368       DO ij=iip2,ip1jm-1
    342369          IF (u_m(ij,l).gt.0.) THEN
    343              zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
     370             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
    344371             u_mq(ij,l)=u_m(ij,l)*
    345      $         min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
     372     $         min(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
    346373          ELSE
    347              zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
     374             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
    348375             u_mq(ij,l)=u_m(ij,l)*
    349      $         min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
     376     $         min(q(ij+1,l,iq)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
    350377          ENDIF
    351378       ENDDO
    … …  
    416443                     i=ijq-(j-1)*iip1
    417444c   accumulation pour les mailles completements advectees
    418                      do while(zu_m.gt.masse(ijq,l))
    419                         u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
    420                         zu_m=zu_m-masse(ijq,l)
     445                     do while(zu_m.gt.masse(ijq,l,iq))
     446                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq)
     447     &                          *masse(ijq,l,iq)
     448                        zu_m=zu_m-masse(ijq,l,iq)
    421449                        i=mod(i-2+iim,iim)+1
    422450                        ijq=(j-1)*iip1+i
    … …  
    424452c   ajout de la maille non completement advectee
    425453                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
    426      &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
     454     &                  (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq))
     455     &                  *dxq(ijq,l))
    427456                  ELSE
    428457                     ijq=ij+1
    429458                     i=ijq-(j-1)*iip1
    430459c   accumulation pour les mailles completements advectees
    431                      do while(-zu_m.gt.masse(ijq,l))
    432                         u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
    433                         zu_m=zu_m+masse(ijq,l)
     460                     do while(-zu_m.gt.masse(ijq,l,iq))
     461                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq)
     462     &                          *masse(ijq,l,iq)
     463                        zu_m=zu_m+masse(ijq,l,iq)
    434464                        i=mod(i,iim)+1
    435465                        ijq=(j-1)*iip1+i
    436466                     ENDDO
    437467c   ajout de la maille non completement advectee
    438                      u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
    439      &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
     468                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)-
     469     &               0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
    440470                  ENDIF
    441471               ENDDO
    … …  
    454484      ENDDO
    455485
     486! CRisi: appel récursif de l'advection sur les fils.
     487! Il faut faire ça avant d'avoir mis à jour q et masse
     488      write(*,*) 'vlspltqs 326: iq,nqfils(iq)=',iq,nqfils(iq)
     489     
     490      if (nqfils(iq).gt.0) then 
     491       do ifils=1,nqdesc(iq)
     492         iq2=iqfils(ifils,iq)
     493         DO l=1,llm
     494          DO ij=iip2,ip1jm
     495           ! On a besoin de q et masse seulement entre iip2 et ip1jm
     496           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     497           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     498          enddo   
     499         enddo
     500        enddo !do ifils=1,nqdesc(iq)
     501        do ifils=1,nqfils(iq)
     502         iq2=iqfils(ifils,iq)
     503         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
     504        enddo !do ifils=1,nqfils(iq)
     505      endif !if (nqfils(iq).gt.0) then
     506! end CRisi
    456507
    457508c   calcul des tendances
    … …  
    459510      DO l=1,llm
    460511         DO ij=iip2+1,ip1jm
    461             new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
    462             q(ij,l)=(q(ij,l)*masse(ij,l)+
     512            new_m=masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l)
     513            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+
    463514     &      u_mq(ij-1,l)-u_mq(ij,l))
    464515     &      /new_m
    465             masse(ij,l)=new_m
     516            masse(ij,l,iq)=new_m
    466517         ENDDO
    467518c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
    468519         DO ij=iip1+iip1,ip1jm,iip1
    469             q(ij-iim,l)=q(ij,l)
    470             masse(ij-iim,l)=masse(ij,l)
    471          ENDDO
    472       ENDDO
     520            q(ij-iim,l,iq)=q(ij,l,iq)
     521            masse(ij-iim,l,iq)=masse(ij,l,iq)
     522         ENDDO
     523      ENDDO
     524
     525      ! retablir les fils en rapport de melange par rapport a l'air:
     526      ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
     527      ! puis on boucle en longitude
     528      if (nqdesc(iq).gt.0) then 
     529       do ifils=1,nqdesc(iq)
     530         iq2=iqfils(ifils,iq) 
     531         DO l=1,llm
     532          DO ij=iip2+1,ip1jm
     533            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     534          enddo
     535          DO ij=iip1+iip1,ip1jm,iip1
     536             q(ij-iim,l,iq2)=q(ij,l,iq2)
     537          enddo ! DO ij=ijb+iip1-1,ije,iip1
     538         enddo !DO l=1,llm
     539        enddo !do ifils=1,nqdesc(iq)
     540      endif !if (nqfils(iq).gt.0) then
    473541
    474542c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
    … …  
    478546      RETURN
    479547      END
    480       SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat)
     548      SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat,iq)
     549      USE infotrac, ONLY : nqtot,nqfils,nqdesc,iqfils ! CRisi
    481550c
    482551c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    502571c   Arguments:
    503572c   ----------
    504       REAL masse(ip1jmp1,llm),pente_max
     573      REAL masse(ip1jmp1,llm,nqtot),pente_max
    505574      REAL masse_adv_v( ip1jm,llm)
    506       REAL q(ip1jmp1,llm)
     575      REAL q(ip1jmp1,llm,nqtot)
    507576      REAL qsat(ip1jmp1,llm)
     577      INTEGER iq ! CRisi
    508578c
    509579c      Local
    … …  
    529599      SAVE sinlon,coslon,sinlondlon,coslondlon
    530600      SAVE airej2,airejjm
     601
     602      REAL masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
     603      INTEGER ifils,iq2 ! CRisi
    531604c
    532605c
    … …  
    567640
    568641      DO i = 1, iim
    569       airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
    570       airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
     642      airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
     643      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
    571644      ENDDO
    572645      qpns   = SSUM( iim,  airescb ,1 ) / airej2
    … …  
    576649
    577650      DO ij=1,ip1jm
    578          dyqv(ij)=q(ij,l)-q(ij+iip1,l)
     651         dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
    579652         adyqv(ij)=abs(dyqv(ij))
    580653      ENDDO
    … …  
    591664
    592665      DO ij=1,iip1
    593          dyq(ij,l)=qpns-q(ij+iip1,l)
    594          dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
     666         dyq(ij,l)=qpns-q(ij+iip1,l,iq)
     667         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
    595668      ENDDO
    596669
    … …  
    710783       DO ij=1,ip1jm
    711784         IF( masse_adv_v(ij,l).GT.0. ) THEN
    712            qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l )  +
    713      ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
     785           qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l,iq )  +
     786     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)
     787     ,      /masse(ij+iip1,l,iq)))
    714788         ELSE
    715               qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
    716      ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
     789              qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l,iq) - dyq(ij,l) *
     790     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l,iq)) )
    717791         ENDIF
    718792          qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
    … …  
    721795
    722796
     797! CRisi: appel récursif de l'advection sur les fils.
     798! Il faut faire ça avant d'avoir mis à jour q et masse
     799      write(*,*) 'vlyqs 689: iq,nqfils(iq)=',iq,nqfils(iq)
     800   
     801      if (nqfils(iq).gt.0) then 
     802       do ifils=1,nqdesc(iq)
     803         iq2=iqfils(ifils,iq)
     804         DO l=1,llm
     805         DO ij=1,ip1jmp1
     806           masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
     807           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)     
     808          enddo   
     809         enddo
     810        enddo !do ifils=1,nqdesc(iq)
     811
     812        do ifils=1,nqfils(iq)
     813         iq2=iqfils(ifils,iq)
     814         write(*,*) 'vlyqs 783: appel rec de vly, iq2=',iq2
     815         call vly(Ratio,pente_max,masseq,qbyv,iq2)
     816        enddo !do ifils=1,nqfils(iq)
     817      endif !if (nqfils(iq).gt.0) then
     818
    723819      DO l=1,llm
    724820         DO ij=iip2,ip1jm
    725             newmasse=masse(ij,l)
     821            newmasse=masse(ij,l,iq)
    726822     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
    727             q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
    728      &         /newmasse
    729             masse(ij,l)=newmasse
     823            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l)
     824     &         -qbyv(ij-iip1,l))/newmasse
     825            masse(ij,l,iq)=newmasse
    730826         ENDDO
    731827c.-. ancienne version
    … …  
    733829         convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
    734830         DO ij = 1,iip1
    735             newmasse=masse(ij,l)+convmpn*aire(ij)
    736             q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
     831            newmasse=masse(ij,l,iq)+convmpn*aire(ij)
     832            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convpn*aire(ij))/
    737833     &               newmasse
    738             masse(ij,l)=newmasse
     834            masse(ij,l,iq)=newmasse
    739835         ENDDO
    740836         convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
    741837         convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
    742838         DO ij = ip1jm+1,ip1jmp1
    743             newmasse=masse(ij,l)+convmps*aire(ij)
    744             q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
     839            newmasse=masse(ij,l,iq)+convmps*aire(ij)
     840            q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+convps*aire(ij))/
    745841     &               newmasse
    746             masse(ij,l)=newmasse
     842            masse(ij,l,iq)=newmasse
    747843         ENDDO
    748844c.-. fin ancienne version
    … …  
    757853c        DO ij = 1,iip1
    758854c           q(ij,l)=newq
    759 c           masse(ij,l)=newmasse*aire(ij)
     855c           masse(ij,l,iq)=newmasse*aire(ij)
    760856c        ENDDO
    761857c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    … …  
    767863c        DO ij = ip1jm+1,ip1jmp1
    768864c           q(ij,l)=newq
    769 c           masse(ij,l)=newmasse*aire(ij)
     865c           masse(ij,l,iq)=newmasse*aire(ij)
    770866c        ENDDO
    771867c._. fin nouvelle version
    772868      ENDDO
    773869
     870      write(*,*) 'vly 866'
     871
     872! retablir les fils en rapport de melange par rapport a l'air:
     873      if (nqdesc(iq).gt.0) then 
     874       do ifils=1,nqdesc(iq)
     875         iq2=iqfils(ifils,iq) 
     876         DO l=1,llm
     877          DO ij=1,ip1jmp1
     878            q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)           
     879          enddo
     880         enddo
     881        enddo !do ifils=1,nqdesc(iq)
     882      endif !if (nqfils(iq).gt.0) then
     883      write(*,*) 'vly 879'
     884
    774885      RETURN
    775886      END
Note: See TracChangeset for help on using the changeset viewer.