/[lmdze]/trunk/dyn3d/vlspltqs.f
ViewVC logotype

Diff of /trunk/dyn3d/vlspltqs.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/dyn3d/vlspltqs.f revision 39 by guez, Tue Jan 25 15:11:05 2011 UTC trunk/libf/dyn3d/vlspltqs.f90 revision 43 by guez, Fri Apr 8 12:43:31 2011 UTC
# Line 1  Line 1 
1  !  !
2  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/vlspltqs.F,v 1.2 2005/02/24 12:16:57 fairhead Exp $  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/vlspltqs.F,v 1.2 2005/02/24 12:16:57 fairhead Exp $
3  !  !
4         SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,         SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt, &
5       ,                                  p,pk,teta                 )                                          p,pk,teta                 )
6  c  !
7  c     Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron  !     Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron
8  c  !
9  c    ********************************************************************  !    ********************************************************************
10  c          Shema  d'advection " pseudo amont " .  !          Shema  d'advection " pseudo amont " .
11  c      + test sur humidite specifique: Q advecte< Qsat aval  !      + test sur humidite specifique: Q advecte< Qsat aval
12  c                   (F. Codron, 10/99)  !                   (F. Codron, 10/99)
13  c    ********************************************************************  !    ********************************************************************
14  c     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....  !     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
15  c  !
16  c     pente_max facteur de limitation des pentes: 2 en general  !     pente_max facteur de limitation des pentes: 2 en general
17  c                                                0 pour un schema amont  !                                                0 pour un schema amont
18  c     pbaru,pbarv,w flux de masse en u ,v ,w  !     pbaru,pbarv,w flux de masse en u ,v ,w
19  c     pdt pas de temps  !     pdt pas de temps
20  c  !
21  c     teta temperature potentielle, p pression aux interfaces,  !     teta temperature potentielle, p pression aux interfaces,
22  c     pk exner au milieu des couches necessaire pour calculer Qsat  !     pk exner au milieu des couches necessaire pour calculer Qsat
23  c   --------------------------------------------------------------------  !   --------------------------------------------------------------------
24        use dimens_m        use dimens_m
25        use paramet_m        use paramet_m
26        use comconst        use comconst
27        use comvert        use comvert
28        use logic        use logic
29        IMPLICIT NONE        IMPLICIT NONE
30  c  !
31    
32  c  !
33  c   Arguments:  !   Arguments:
34  c   ----------  !   ----------
35        REAL masse(ip1jmp1,llm),pente_max        REAL masse(ip1jmp1,llm),pente_max
36        REAL, intent(in):: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)        REAL, intent(in):: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
37        REAL q(ip1jmp1,llm)        REAL q(ip1jmp1,llm)
# Line 39  c   ---------- Line 39  c   ----------
39        real, intent(in):: pdt        real, intent(in):: pdt
40        REAL, intent(in):: p(ip1jmp1,llmp1)        REAL, intent(in):: p(ip1jmp1,llmp1)
41        real teta(ip1jmp1,llm),pk(ip1jmp1,llm)        real teta(ip1jmp1,llm),pk(ip1jmp1,llm)
42  c  !
43  c      Local  !      Local
44  c   ---------  !   ---------
45  c  !
46        INTEGER i,ij,l,j,ii        INTEGER i,ij,l,j,ii
47  c  !
48        REAL qsat(ip1jmp1,llm)        REAL qsat(ip1jmp1,llm)
49        REAL zm(ip1jmp1,llm)        REAL zm(ip1jmp1,llm)
50        REAL mu(ip1jmp1,llm)        REAL mu(ip1jmp1,llm)
# Line 62  c Line 62  c
62        DATA testcpu/.false./        DATA testcpu/.false./
63        DATA temps1,temps2,temps3/0.,0.,0./        DATA temps1,temps2,temps3/0.,0.,0./
64    
65  c--pour rapport de melange saturant--  !--pour rapport de melange saturant--
66    
67        REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play        REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
68        REAL ptarg,pdelarg,foeew,zdelta        REAL ptarg,pdelarg,foeew,zdelta
69        REAL tempe(ip1jmp1)        REAL tempe(ip1jmp1)
70    
71  c    fonction psat(T)  !    fonction psat(T)
72    
73         FOEEW ( PTARG,PDELARG ) = EXP (         FOEEW ( PTARG,PDELARG ) = EXP ( &
74       *          (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT)                  (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) &
75       * / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )         / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
76    
77          r2es  = 380.11733          r2es  = 380.11733
78          r3les = 17.269          r3les = 17.269
79          r3ies = 21.875          r3ies = 21.875
80          r4les = 35.86          r4les = 35.86
# Line 82  c    fonction psat(T) Line 82  c    fonction psat(T)
82          retv = 0.6077667          retv = 0.6077667
83          rtt  = 273.16          rtt  = 273.16
84    
85  c-- Calcul de Qsat en chaque point  !-- Calcul de Qsat en chaque point
86  c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2  !-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
87  c   pour eviter une exponentielle.  !   pour eviter une exponentielle.
88          DO l = 1, llm          DO l = 1, llm
89           DO ij = 1, ip1jmp1           DO ij = 1, ip1jmp1
90            tempe(ij) = teta(ij,l) * pk(ij,l) /cpp            tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
# Line 118  c   pour eviter une exponentielle. Line 118  c   pour eviter une exponentielle.
118        CALL SCOPY(ijp1llm,q,1,zq,1)        CALL SCOPY(ijp1llm,q,1,zq,1)
119        CALL SCOPY(ijp1llm,masse,1,zm,1)        CALL SCOPY(ijp1llm,masse,1,zm,1)
120    
121  c      call minmaxq(zq,qmin,qmax,'avant vlxqs     ')  !      call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
122        call vlxqs(zq,pente_max,zm,mu,qsat)        call vlxqs(zq,pente_max,zm,mu,qsat)
123    
124    
125  c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')  !     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
126    
127        call vlyqs(zq,pente_max,zm,mv,qsat)        call vlyqs(zq,pente_max,zm,mv,qsat)
128    
129    
130  c      call minmaxq(zq,qmin,qmax,'avant vlz     ')  !      call minmaxq(zq,qmin,qmax,'avant vlz     ')
131    
132        call vlz(zq,pente_max,zm,mw)        call vlz(zq,pente_max,zm,mw)
133    
134    
135  c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')  !     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
136  c     call minmaxq(zm,qmin,qmax,'M avant vlyqs     ')  !     call minmaxq(zm,qmin,qmax,'M avant vlyqs     ')
137    
138        call vlyqs(zq,pente_max,zm,mv,qsat)        call vlyqs(zq,pente_max,zm,mv,qsat)
139    
140    
141  c     call minmaxq(zq,qmin,qmax,'avant vlxqs     ')  !     call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
142  c     call minmaxq(zm,qmin,qmax,'M avant vlxqs     ')  !     call minmaxq(zm,qmin,qmax,'M avant vlxqs     ')
143    
144        call vlxqs(zq,pente_max,zm,mu,qsat)        call vlxqs(zq,pente_max,zm,mu,qsat)
145    
146  c     call minmaxq(zq,qmin,qmax,'apres vlxqs     ')  !     call minmaxq(zq,qmin,qmax,'apres vlxqs     ')
147  c     call minmaxq(zm,qmin,qmax,'M apres vlxqs     ')  !     call minmaxq(zm,qmin,qmax,'M apres vlxqs     ')
148    
149    
150        DO l=1,llm        DO l=1,llm
# Line 157  c     call minmaxq(zm,qmin,qmax,'M apres Line 157  c     call minmaxq(zm,qmin,qmax,'M apres
157        ENDDO        ENDDO
158    
159        RETURN        RETURN
       END  
       SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat)  
 c  
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget  
 c  
 c    ********************************************************************  
 c     Shema  d'advection " pseudo amont " .  
 c    ********************************************************************  
 c  
 c   --------------------------------------------------------------------  
       use dimens_m  
       use paramet_m  
       use comconst  
       use comvert  
       use logic  
       IMPLICIT NONE  
 c  
 c  
 c  
 c   Arguments:  
 c   ----------  
       REAL masse(ip1jmp1,llm),pente_max  
       REAL u_m( ip1jmp1,llm )  
       REAL q(ip1jmp1,llm)  
       REAL qsat(ip1jmp1,llm)  
 c  
 c      Local  
 c   ---------  
 c  
       INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju  
       INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)  
 c  
       REAL new_m,zu_m,zdum(ip1jmp1,llm)  
       REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)  
       REAL zz(ip1jmp1)  
       REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)  
       REAL u_mq(ip1jmp1,llm)  
   
       Logical first,testcpu  
       SAVE first,testcpu  
   
       REAL      SSUM  
       REAL temps0,temps1,temps2,temps3,temps4,temps5  
       SAVE temps0,temps1,temps2,temps3,temps4,temps5  
   
   
       DATA first,testcpu/.true.,.false./  
   
       IF(first) THEN  
          temps1=0.  
          temps2=0.  
          temps3=0.  
          temps4=0.  
          temps5=0.  
          first=.false.  
       ENDIF  
   
 c   calcul de la pente a droite et a gauche de la maille  
   
   
       IF (pente_max.gt.-1.e-5) THEN  
 c     IF (pente_max.gt.10) THEN  
   
 c   calcul des pentes avec limitation, Van Leer scheme I:  
 c   -----------------------------------------------------  
   
 c   calcul de la pente aux points u  
          DO l = 1, llm  
             DO ij=iip2,ip1jm-1  
                dxqu(ij)=q(ij+1,l)-q(ij,l)  
 c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'  
 c              sigu(ij)=u_m(ij,l)/masse(ij,l)  
             ENDDO  
             DO ij=iip1+iip1,ip1jm,iip1  
                dxqu(ij)=dxqu(ij-iim)  
 c              sigu(ij)=sigu(ij-iim)  
             ENDDO  
   
             DO ij=iip2,ip1jm  
                adxqu(ij)=abs(dxqu(ij))  
             ENDDO  
   
 c   calcul de la pente maximum dans la maille en valeur absolue  
   
             DO ij=iip2+1,ip1jm  
                dxqmax(ij,l)=pente_max*  
      ,      min(adxqu(ij-1),adxqu(ij))  
 c limitation subtile  
 c    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))  
             
   
             ENDDO  
   
             DO ij=iip1+iip1,ip1jm,iip1  
                dxqmax(ij-iim,l)=dxqmax(ij,l)  
             ENDDO  
   
             DO ij=iip2+1,ip1jm  
                IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN  
                   dxq(ij,l)=dxqu(ij-1)+dxqu(ij)  
                ELSE  
 c   extremum local  
                   dxq(ij,l)=0.  
                ENDIF  
                dxq(ij,l)=0.5*dxq(ij,l)  
                dxq(ij,l)=  
      ,         sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))  
             ENDDO  
   
          ENDDO ! l=1,llm  
   
       ELSE ! (pente_max.lt.-1.e-5)  
   
 c   Pentes produits:  
 c   ----------------  
   
          DO l = 1, llm  
             DO ij=iip2,ip1jm-1  
                dxqu(ij)=q(ij+1,l)-q(ij,l)  
             ENDDO  
             DO ij=iip1+iip1,ip1jm,iip1  
                dxqu(ij)=dxqu(ij-iim)  
             ENDDO  
   
             DO ij=iip2+1,ip1jm  
                zz(ij)=dxqu(ij-1)*dxqu(ij)  
                zz(ij)=zz(ij)+zz(ij)  
                IF(zz(ij).gt.0) THEN  
                   dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))  
                ELSE  
 c   extremum local  
                   dxq(ij,l)=0.  
                ENDIF  
             ENDDO  
   
          ENDDO  
   
       ENDIF ! (pente_max.lt.-1.e-5)  
   
 c   bouclage de la pente en iip1:  
 c   -----------------------------  
   
       DO l=1,llm  
          DO ij=iip1+iip1,ip1jm,iip1  
             dxq(ij-iim,l)=dxq(ij,l)  
          ENDDO  
   
          DO ij=1,ip1jmp1  
             iadvplus(ij,l)=0  
          ENDDO  
   
       ENDDO  
   
   
 c   calcul des flux a gauche et a droite  
   
 c   on cumule le flux correspondant a toutes les mailles dont la masse  
 c   au travers de la paroi pENDant le pas de temps.  
 c   le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)  
       DO l=1,llm  
        DO ij=iip2,ip1jm-1  
           IF (u_m(ij,l).gt.0.) THEN  
              zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)  
              u_mq(ij,l)=u_m(ij,l)*  
      $         min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))  
           ELSE  
              zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)  
              u_mq(ij,l)=u_m(ij,l)*  
      $         min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))  
           ENDIF  
        ENDDO  
       ENDDO  
   
   
 c   detection des points ou on advecte plus que la masse de la  
 c   maille  
       DO l=1,llm  
          DO ij=iip2,ip1jm-1  
             IF(zdum(ij,l).lt.0) THEN  
                iadvplus(ij,l)=1  
                u_mq(ij,l)=0.  
             ENDIF  
          ENDDO  
       ENDDO  
       DO l=1,llm  
        DO ij=iip1+iip1,ip1jm,iip1  
           iadvplus(ij,l)=iadvplus(ij-iim,l)  
        ENDDO  
       ENDDO  
   
   
   
 c   traitement special pour le cas ou on advecte en longitude plus que le  
 c   contenu de la maille.  
 c   cette partie est mal vectorisee.  
   
 c   pas d'influence de la pression saturante (pour l'instant)  
   
 c  calcul du nombre de maille sur lequel on advecte plus que la maille.  
   
       n0=0  
       DO l=1,llm  
          nl(l)=0  
          DO ij=iip2,ip1jm  
             nl(l)=nl(l)+iadvplus(ij,l)  
          ENDDO  
          n0=n0+nl(l)  
       ENDDO  
   
       IF(n0.gt.0) THEN  
          DO l=1,llm  
             IF(nl(l).gt.0) THEN  
                iju=0  
 c   indicage des mailles concernees par le traitement special  
                DO ij=iip2,ip1jm  
                   IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN  
                      iju=iju+1  
                      indu(iju)=ij  
                   ENDIF  
                ENDDO  
                niju=iju  
   
 c  traitement des mailles  
                DO iju=1,niju  
                   ij=indu(iju)  
                   j=(ij-1)/iip1+1  
                   zu_m=u_m(ij,l)  
                   u_mq(ij,l)=0.  
                   IF(zu_m.gt.0.) THEN  
                      ijq=ij  
                      i=ijq-(j-1)*iip1  
 c   accumulation pour les mailles completements advectees  
                      do while(zu_m.gt.masse(ijq,l))  
                         u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)  
                         zu_m=zu_m-masse(ijq,l)  
                         i=mod(i-2+iim,iim)+1  
                         ijq=(j-1)*iip1+i  
                      ENDDO  
 c   ajout de la maille non completement advectee  
                      u_mq(ij,l)=u_mq(ij,l)+zu_m*  
      &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))  
                   ELSE  
                      ijq=ij+1  
                      i=ijq-(j-1)*iip1  
 c   accumulation pour les mailles completements advectees  
                      do while(-zu_m.gt.masse(ijq,l))  
                         u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)  
                         zu_m=zu_m+masse(ijq,l)  
                         i=mod(i,iim)+1  
                         ijq=(j-1)*iip1+i  
                      ENDDO  
 c   ajout de la maille non completement advectee  
                      u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-  
      &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))  
                   ENDIF  
                ENDDO  
             ENDIF  
          ENDDO  
       ENDIF  ! n0.gt.0  
   
   
   
 c   bouclage en latitude  
   
       DO l=1,llm  
         DO ij=iip1+iip1,ip1jm,iip1  
            u_mq(ij,l)=u_mq(ij-iim,l)  
         ENDDO  
       ENDDO  
   
   
 c   calcul des tendances  
   
       DO l=1,llm  
          DO ij=iip2+1,ip1jm  
             new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)  
             q(ij,l)=(q(ij,l)*masse(ij,l)+  
      &      u_mq(ij-1,l)-u_mq(ij,l))  
      &      /new_m  
             masse(ij,l)=new_m  
          ENDDO  
 c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)  
          DO ij=iip1+iip1,ip1jm,iip1  
             q(ij-iim,l)=q(ij,l)  
             masse(ij-iim,l)=masse(ij,l)  
          ENDDO  
       ENDDO  
   
 c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)  
 c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)  
   
   
       RETURN  
       END  
       SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat)  
 c  
 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget  
 c  
 c    ********************************************************************  
 c     Shema  d'advection " pseudo amont " .  
 c    ********************************************************************  
 c     q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....  
 c     qsat             est   un argument de sortie pour le s-pg ....  
 c  
 c  
 c   --------------------------------------------------------------------  
 c  
       use dimens_m  
       use paramet_m  
       use comconst  
       use comvert  
       use logic  
       use comgeom  
       USE nr_util, ONLY : pi  
       IMPLICIT NONE  
 c  
 c  
 c   Arguments:  
 c   ----------  
       REAL masse(ip1jmp1,llm),pente_max  
       REAL masse_adv_v( ip1jm,llm)  
       REAL q(ip1jmp1,llm)  
       REAL qsat(ip1jmp1,llm)  
 c  
 c      Local  
 c   ---------  
 c  
       INTEGER i,ij,l  
 c  
       REAL airej2,airejjm,airescb(iim),airesch(iim)  
       REAL dyq(ip1jmp1,llm),dyqv(ip1jm)  
       REAL adyqv(ip1jm),dyqmax(ip1jmp1)  
       REAL qbyv(ip1jm,llm)  
   
       REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs  
 c     REAL newq,oldmasse  
       Logical first,testcpu  
       REAL temps0,temps1,temps2,temps3,temps4,temps5  
       SAVE temps0,temps1,temps2,temps3,temps4,temps5  
       SAVE first,testcpu  
   
       REAL convpn,convps,convmpn,convmps  
       REAL sinlon(iip1),sinlondlon(iip1)  
       REAL coslon(iip1),coslondlon(iip1)  
       SAVE sinlon,coslon,sinlondlon,coslondlon  
       SAVE airej2,airejjm  
 c  
 c  
       REAL      SSUM  
   
       DATA first,testcpu/.true.,.false./  
       DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./  
   
       IF(first) THEN  
          PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '  
          first=.false.  
          do i=2,iip1  
             coslon(i)=cos(rlonv(i))  
             sinlon(i)=sin(rlonv(i))  
             coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi  
             sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi  
          ENDDO  
          coslon(1)=coslon(iip1)  
          coslondlon(1)=coslondlon(iip1)  
          sinlon(1)=sinlon(iip1)  
          sinlondlon(1)=sinlondlon(iip1)  
          airej2 = SSUM( iim, aire(iip2), 1 )  
          airejjm= SSUM( iim, aire(ip1jm -iim), 1 )  
       ENDIF  
   
 c  
   
   
       DO l = 1, llm  
 c  
 c   --------------------------------  
 c      CALCUL EN LATITUDE  
 c   --------------------------------  
   
 c   On commence par calculer la valeur du traceur moyenne sur le premier cercle  
 c   de latitude autour du pole (qpns pour le pole nord et qpsn pour  
 c    le pole nord) qui sera utilisee pour evaluer les pentes au pole.  
   
       DO i = 1, iim  
       airescb(i) = aire(i+ iip1) * q(i+ iip1,l)  
       airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)  
       ENDDO  
       qpns   = SSUM( iim,  airescb ,1 ) / airej2  
       qpsn   = SSUM( iim,  airesch ,1 ) / airejjm  
   
 c   calcul des pentes aux points v  
   
       DO ij=1,ip1jm  
          dyqv(ij)=q(ij,l)-q(ij+iip1,l)  
          adyqv(ij)=abs(dyqv(ij))  
       ENDDO  
   
 c   calcul des pentes aux points scalaires  
   
       DO ij=iip2,ip1jm  
          dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))  
          dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))  
          dyqmax(ij)=pente_max*dyqmax(ij)  
       ENDDO  
   
 c   calcul des pentes aux poles  
   
       DO ij=1,iip1  
          dyq(ij,l)=qpns-q(ij+iip1,l)  
          dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn  
       ENDDO  
   
 c   filtrage de la derivee  
       dyn1=0.  
       dys1=0.  
       dyn2=0.  
       dys2=0.  
       DO ij=1,iim  
          dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)  
          dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)  
          dyn2=dyn2+coslondlon(ij)*dyq(ij,l)  
          dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)  
       ENDDO  
       DO ij=1,iip1  
          dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)  
          dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)  
       ENDDO  
   
 c   calcul des pentes limites aux poles  
   
       fn=1.  
       fs=1.  
       DO ij=1,iim  
          IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN  
             fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)  
          ENDIF  
       IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN  
          fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)  
          ENDIF  
       ENDDO  
       DO ij=1,iip1  
          dyq(ij,l)=fn*dyq(ij,l)  
          dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)  
       ENDDO  
   
 c   calcul des pentes limitees  
   
       DO ij=iip2,ip1jm  
          IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN  
             dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))  
          ELSE  
             dyq(ij,l)=0.  
          ENDIF  
       ENDDO  
   
       ENDDO  
   
       DO l=1,llm  
        DO ij=1,ip1jm  
          IF( masse_adv_v(ij,l).GT.0. ) THEN  
            qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l )  +  
      ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))  
          ELSE  
               qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *  
      ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )  
          ENDIF  
           qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)  
        ENDDO  
       ENDDO  
   
   
       DO l=1,llm  
          DO ij=iip2,ip1jm  
             newmasse=masse(ij,l)  
      &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)  
             q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))  
      &         /newmasse  
             masse(ij,l)=newmasse  
          ENDDO  
 c.-. ancienne version  
          convpn=SSUM(iim,qbyv(1,l),1)/apoln  
          convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln  
          DO ij = 1,iip1  
             newmasse=masse(ij,l)+convmpn*aire(ij)  
             q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/  
      &               newmasse  
             masse(ij,l)=newmasse  
          ENDDO  
          convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols  
          convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols  
          DO ij = ip1jm+1,ip1jmp1  
             newmasse=masse(ij,l)+convmps*aire(ij)  
             q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/  
      &               newmasse  
             masse(ij,l)=newmasse  
          ENDDO  
 c.-. fin ancienne version  
   
 c._. nouvelle version  
 c        convpn=SSUM(iim,qbyv(1,l),1)  
 c        convmpn=ssum(iim,masse_adv_v(1,l),1)  
 c        oldmasse=ssum(iim,masse(1,l),1)  
 c        newmasse=oldmasse+convmpn  
 c        newq=(q(1,l)*oldmasse+convpn)/newmasse  
 c        newmasse=newmasse/apoln  
 c        DO ij = 1,iip1  
 c           q(ij,l)=newq  
 c           masse(ij,l)=newmasse*aire(ij)  
 c        ENDDO  
 c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)  
 c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)  
 c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)  
 c        newmasse=oldmasse+convmps  
 c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse  
 c        newmasse=newmasse/apols  
 c        DO ij = ip1jm+1,ip1jmp1  
 c           q(ij,l)=newq  
 c           masse(ij,l)=newmasse*aire(ij)  
 c        ENDDO  
 c._. fin nouvelle version  
       ENDDO  
   
       RETURN  
160        END        END
Legend:
Removed from v.39  
changed lines
  Added in v.43
  ViewVC Help
Powered by ViewVC 1.1.21