/[lmdze]/trunk/Sources/dyn3d/Vlsplt/vlx.f
ViewVC logotype

Diff of /trunk/Sources/dyn3d/Vlsplt/vlx.f

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

revision 156 by guez, Wed Apr 29 15:47:56 2015 UTC revision 157 by guez, Mon Jul 20 16:01:49 2015 UTC
# Line 1  Line 1 
1        SUBROUTINE vlx(q,pente_max,masse,u_m)  module vlx_m
2    
3  !     Auteurs:   P.Le Van, F.Hourdin, F.Forget    IMPLICIT NONE
4  !  
5  !   *******************************************************  contains
6  !     Shema  d'advection " pseudo amont " .  
7  !    ************************************************************    SUBROUTINE vlx(q, pente_max, masse, u_m)
8  !     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le  
9  !     s-pg ....      ! Authors: P. Le Van, F. Hourdin, F. Forget
10  !  
11  !      ! Sch\'ema d'advection "pseudo-amont".
12  !   --------------------------------------------------------------  
13        use dimens_m      use dimens_m, only: iim, llm
14        use paramet_m      use paramet_m, only: ip1jmp1, iip1, iip2, ip1jm
15        use comconst  
16        use disvert_m      REAL, intent(inout):: q(ip1jmp1, llm)
17        use conf_gcm_m      REAL, intent(in):: pente_max
18        IMPLICIT NONE      real, intent(inout):: masse(ip1jmp1, llm)
19  !      REAL, intent(in):: u_m(ip1jmp1, llm)
20  !  
21  !      ! Local:
22  !   Arguments:  
23  !   ----------      INTEGER ij, l, j, i, iju, ijq, indu(ip1jmp1), niju
24        REAL masse(ip1jmp1,llm),pente_max      INTEGER n0, iadvplus(ip1jmp1, llm), nl(llm)
25        REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)      REAL new_m, zu_m, zdum(ip1jmp1, llm)
26        REAL q(ip1jmp1,llm)      REAL dxq(ip1jmp1, llm), dxqu(iip2:ip1jm)
27        REAL w(ip1jmp1,llm)      REAL zz(ip1jmp1)
28  !      REAL adxqu(iip2:ip1jm), dxqmax(ip1jmp1, llm)
29  !      Local      REAL u_mq(ip1jmp1, llm)
30  !   ---------  
31  !      !-----------------------------------------------------
32        INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju  
33        INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)      ! calcul de la pente a droite et a gauche de la maille
34  !  
35        REAL new_m,zu_m,zdum(ip1jmp1,llm)      IF (pente_max > - 1e-5) THEN
36        REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)         ! calcul des pentes avec limitation, Van Leer scheme I:
37        REAL zz(ip1jmp1)  
38        REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)         ! calcul de la pente aux points u
39        REAL u_mq(ip1jmp1,llm)         DO l = 1, llm
40              DO ij = iip2, ip1jm - 1
41        Logical extremum,first,testcpu               dxqu(ij) = q(ij + 1, l) - q(ij, l)
42        SAVE first,testcpu            ENDDO
43              DO ij = iip1 + iip1, ip1jm, iip1
44        REAL      SSUM               dxqu(ij) = dxqu(ij - iim)
45        REAL temps0,temps1,temps2,temps3,temps4,temps5,second            ENDDO
46        SAVE temps0,temps1,temps2,temps3,temps4,temps5  
47              DO ij = iip2, ip1jm
48        REAL z1,z2,z3               adxqu(ij) = abs(dxqu(ij))
49              ENDDO
50        DATA first,testcpu/.true.,.false./  
51              ! calcul de la pente maximum dans la maille en valeur absolue
52        IF(first) THEN  
53           temps1=0.            DO ij = iip2 + 1, ip1jm
54           temps2=0.               dxqmax(ij, l) = pente_max * min(adxqu(ij - 1), adxqu(ij))
55           temps3=0.            ENDDO
56           temps4=0.  
57           temps5=0.            DO ij = iip1 + iip1, ip1jm, iip1
58           first=.false.               dxqmax(ij - iim, l) = dxqmax(ij, l)
59        ENDIF            ENDDO
60    
61  !   calcul de la pente a droite et a gauche de la maille            DO ij = iip2 + 1, ip1jm
62                 IF (dxqu(ij - 1) * dxqu(ij) > 0.) THEN
63                    dxq(ij, l) = dxqu(ij - 1) + dxqu(ij)
64        IF (pente_max.gt.-1.e-5) THEN               ELSE
65  !       IF (pente_max.gt.10) THEN                  ! extremum local
66                    dxq(ij, l) = 0.
67  !   calcul des pentes avec limitation, Van Leer scheme I:               ENDIF
68  !   -----------------------------------------------------               dxq(ij, l) = 0.5 * dxq(ij, l)
69                 dxq(ij, l) =  sign(min(abs(dxq(ij, l)), dxqmax(ij, l)), dxq(ij, l))
70  !   calcul de la pente aux points u            ENDDO
71           DO l = 1, llm         ENDDO
72              DO ij=iip2,ip1jm-1      ELSE
73                 dxqu(ij)=q(ij+1,l)-q(ij,l)         ! Pentes produits:
74  !              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'  
75  !              sigu(ij)=u_m(ij,l)/masse(ij,l)         DO l = 1, llm
76              ENDDO            DO ij = iip2, ip1jm - 1
77              DO ij=iip1+iip1,ip1jm,iip1               dxqu(ij) = q(ij + 1, l) - q(ij, l)
78                 dxqu(ij)=dxqu(ij-iim)            ENDDO
79  !              sigu(ij)=sigu(ij-iim)            DO ij = iip1 + iip1, ip1jm, iip1
80              ENDDO               dxqu(ij) = dxqu(ij - iim)
81              ENDDO
82              DO ij=iip2,ip1jm  
83                 adxqu(ij)=abs(dxqu(ij))            DO ij = iip2 + 1, ip1jm
84              ENDDO               zz(ij) = dxqu(ij - 1) * dxqu(ij)
85                 zz(ij) = zz(ij) + zz(ij)
86  !   calcul de la pente maximum dans la maille en valeur absolue               IF (zz(ij) > 0) THEN
87                    dxq(ij, l) = zz(ij) / (dxqu(ij - 1) + dxqu(ij))
88              DO ij=iip2+1,ip1jm               ELSE
89                 dxqmax(ij,l)=pente_max*                                  &                  ! extremum local
90       &      min(adxqu(ij-1),adxqu(ij))                  dxq(ij, l) = 0.
91  ! limitation subtile               ENDIF
92  !    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))            ENDDO
93           ENDDO
94        ENDIF
95              ENDDO  
96        ! bouclage de la pente en iip1:
97              DO ij=iip1+iip1,ip1jm,iip1  
98                 dxqmax(ij-iim,l)=dxqmax(ij,l)      DO l = 1, llm
99              ENDDO         DO ij = iip1 + iip1, ip1jm, iip1
100              dxq(ij - iim, l) = dxq(ij, l)
101              DO ij=iip2+1,ip1jm         ENDDO
102                 IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN         DO ij = 1, ip1jmp1
103                    dxq(ij,l)=dxqu(ij-1)+dxqu(ij)            iadvplus(ij, l) = 0
104                 ELSE         ENDDO
105  !   extremum local      ENDDO
106                    dxq(ij,l)=0.  
107                 ENDIF      ! calcul des flux a gauche et a droite
108                 dxq(ij,l)=0.5*dxq(ij,l)  
109                 dxq(ij,l)=                                               &      ! on cumule le flux correspondant a toutes les mailles dont la masse
110       &         sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))      ! au travers de la paroi pENDant le pas de temps.
111              ENDDO  
112        DO l = 1, llm
113                 ! l=1,llm         DO ij = iip2, ip1jm - 1
114           ENDDO            IF (u_m(ij, l) > 0.) THEN
115             ! (pente_max.lt.-1.e-5)               zdum(ij, l) = 1. - u_m(ij, l) / masse(ij, l)
116        ELSE               u_mq(ij, l) = u_m(ij, l) &
117                      * (q(ij, l) + 0.5 * zdum(ij, l) * dxq(ij, l))
 !   Pentes produits:  
 !   ----------------  
   
          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  
 !   extremum local  
                   dxq(ij,l)=0.  
                ENDIF  
             ENDDO  
   
          ENDDO  
   
             ! (pente_max.lt.-1.e-5)  
       ENDIF  
   
 !   bouclage de la pente en iip1:  
 !   -----------------------------  
   
       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  
   
 !   calcul des flux a gauche et a droite  
   
 !   on cumule le flux correspondant a toutes les mailles dont la masse  
 !   au travers de la paroi pENDant le pas de temps.  
 !print*,'Cumule ....'  
   
       DO l=1,llm  
        DO ij=iip2,ip1jm-1  
 !      print*,'masse(',ij,')=',masse(ij,l)  
           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)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))  
118            ELSE            ELSE
119               zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)               zdum(ij, l) = 1. + u_m(ij, l) / masse(ij + 1, l)
120               u_mq(ij,l)=u_m(ij,l)                                       &               u_mq(ij, l) = u_m(ij, l) &
121       &            *(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))                    * (q(ij + 1, l) - 0.5 * zdum(ij, l) * dxq(ij + 1, l))
122              ENDIF
123           ENDDO
124        ENDDO
125    
126        ! detection des points ou on advecte plus que la masse de la
127        ! maille
128        DO l = 1, llm
129           DO ij = iip2, ip1jm - 1
130              IF (zdum(ij, l) < 0) THEN
131                 iadvplus(ij, l) = 1
132                 u_mq(ij, l) = 0.
133              ENDIF
134           ENDDO
135        ENDDO
136    
137        DO l = 1, llm
138           DO ij = iip1 + iip1, ip1jm, iip1
139              iadvplus(ij, l) = iadvplus(ij - iim, l)
140           ENDDO
141        ENDDO
142    
143        ! traitement special pour le cas ou on advecte en longitude plus
144        ! que le contenu de la maille.  cette partie est mal vectorisee.
145    
146        ! calcul du nombre de maille sur lequel on advecte plus que la maille.
147    
148        n0 = 0
149        DO l = 1, llm
150           nl(l) = 0
151           DO ij = iip2, ip1jm
152              nl(l) = nl(l) + iadvplus(ij, l)
153           ENDDO
154           n0 = n0 + nl(l)
155        ENDDO
156    
157        IF (n0 > 0) THEN
158           DO l = 1, llm
159              IF (nl(l) > 0) THEN
160                 iju = 0
161                 ! indicage des mailles concernees par le traitement special
162                 DO ij = iip2, ip1jm
163                    IF (iadvplus(ij, l) == 1.and.mod(ij, iip1) /= 0) THEN
164                       iju = iju + 1
165                       indu(iju) = ij
166                    ENDIF
167                 ENDDO
168                 niju = iju
169    
170                 ! traitement des mailles
171                 DO iju = 1, niju
172                    ij = indu(iju)
173                    j = (ij - 1) / iip1 + 1
174                    zu_m = u_m(ij, l)
175                    u_mq(ij, l) = 0.
176                    IF (zu_m > 0.) THEN
177                       ijq = ij
178                       i = ijq - (j - 1) * iip1
179                       ! accumulation pour les mailles completements advectees
180                       do while(zu_m > masse(ijq, l))
181                          u_mq(ij, l) = u_mq(ij, l) + q(ijq, l) * masse(ijq, l)
182                          zu_m = zu_m - masse(ijq, l)
183                          i = mod(i - 2 + iim, iim) + 1
184                          ijq = (j - 1) * iip1 + i
185                       ENDDO
186                       ! ajout de la maille non completement advectee
187                       u_mq(ij, l) = u_mq(ij, l) + zu_m * (q(ijq, l) + 0.5 * (1. &
188                            - zu_m / masse(ijq, l)) * dxq(ijq, l))
189                    ELSE
190                       ijq = ij + 1
191                       i = ijq - (j - 1) * iip1
192                       ! accumulation pour les mailles completements advectees
193                       do while(- zu_m > masse(ijq, l))
194                          u_mq(ij, l) = u_mq(ij, l) - q(ijq, l) * masse(ijq, l)
195                          zu_m = zu_m + masse(ijq, l)
196                          i = mod(i, iim) + 1
197                          ijq = (j - 1) * iip1 + i
198                       ENDDO
199                       ! ajout de la maille non completement advectee
200                       u_mq(ij, l) = u_mq(ij, l) + zu_m * (q(ijq, l) - 0.5 * (1. &
201                            + zu_m / masse(ijq, l)) * dxq(ijq, l))
202                    ENDIF
203                 ENDDO
204            ENDIF            ENDIF
205         ENDDO         ENDDO
206        ENDDO      ENDIF
 !      stop  
207    
208  !      go to 9999      ! bouclage en latitude
209  !   detection des points ou on advecte plus que la masse de la      DO l = 1, llm
210  !   maille         DO ij = iip1 + iip1, ip1jm, iip1
211        DO l=1,llm            u_mq(ij, l) = u_mq(ij - iim, l)
212           DO ij=iip2,ip1jm-1         ENDDO
213              IF(zdum(ij,l).lt.0) THEN      ENDDO
214                 iadvplus(ij,l)=1  
215                 u_mq(ij,l)=0.      ! calcul des tENDances
216              ENDIF  
217           ENDDO      DO l = 1, llm
218        ENDDO         DO ij = iip2 + 1, ip1jm
219  !print*,'Ok test 1'            new_m = masse(ij, l) + u_m(ij - 1, l) - u_m(ij, l)
220        DO l=1,llm            q(ij, l) = (q(ij, l) * masse(ij, l) + u_mq(ij - 1, l) &
221         DO ij=iip1+iip1,ip1jm,iip1                 - u_mq(ij, l)) / new_m
222            iadvplus(ij,l)=iadvplus(ij-iim,l)            masse(ij, l) = new_m
223         ENDDO         ENDDO
224        ENDDO  
225  ! print*,'Ok test 2'         DO ij = iip1 + iip1, ip1jm, iip1
226              q(ij - iim, l) = q(ij, l)
227              masse(ij - iim, l) = masse(ij, l)
228  !   traitement special pour le cas ou on advecte en longitude plus         ENDDO
229  !     que le      ENDDO
 !   contenu de la maille.  
 !   cette partie est mal vectorisee.  
   
 !  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  
 !C      PRINT*,'Nombre de points pour lesquels on advect plus que le'  
 !C     &       ,'contenu de la maille : ',n0  
   
          DO l=1,llm  
             IF(nl(l).gt.0) THEN  
                iju=0  
 !   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  
 !              PRINT*,'niju,nl',niju,nl(l)  
   
 !  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  
 !   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  
 !   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  
 !   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  
 !   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  
              ! n0.gt.0  
       ENDIF  
  9999   continue  
   
   
 !   bouclage en latitude  
 !print*,'cvant bouclage en latitude'  
       DO l=1,llm  
         DO ij=iip1+iip1,ip1jm,iip1  
            u_mq(ij,l)=u_mq(ij-iim,l)  
         ENDDO  
       ENDDO  
   
   
 !   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  
 !   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  
230    
231      END SUBROUTINE vlx
232    
233        RETURN  end module vlx_m
       END  
Legend:
Removed from v.156  
changed lines
  Added in v.157
  ViewVC Help
Powered by ViewVC 1.1.21