/[lmdze]/trunk/dyn3d/prather.f90
ViewVC logotype

Diff of /trunk/dyn3d/prather.f90

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

trunk/dyn3d/prather.f revision 80 by guez, Fri Nov 15 18:45:49 2013 UTC trunk/dyn3d/prather.f90 revision 81 by guez, Wed Mar 5 14:38:41 2014 UTC
# Line 1  Line 1 
1  !  
2  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/prather.F,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/prather.F,v 1.1.1.1 2004/05/19
3  !  ! 12:53:07 lmdzadmin Exp $
4        SUBROUTINE prather (q,w,masse,pbaru,pbarv,nt,dt)  
5        use dimens_m  SUBROUTINE prather(q, w, masse, pbaru, pbarv, nt, dt)
6        use paramet_m    USE dimens_m
7        use comconst    USE paramet_m
8        use disvert_m    USE comconst
9        use comgeom    USE disvert_m
10        USE nr_util, ONLY : pi    USE comgeom
11        IMPLICIT NONE    USE nr_util, ONLY: pi
12      IMPLICIT NONE
13  c=======================================================================  
14  c   Adaptation LMDZ:  A.Armengaud (LGGE)    ! =======================================================================
15  c   ----------------    ! Adaptation LMDZ:  A.Armengaud (LGGE)
16  c    ! ----------------
17  c   ************************************************  
18  c   Transport des traceurs par la methode de prather    ! ************************************************
19  c   Ref :    ! Transport des traceurs par la methode de prather
20  c    ! Ref :
21  c   ************************************************  
22  c   q,w,pext,pbaru et pbarv : arguments d'entree  pour le s-pg    ! ************************************************
23  c    ! q,w,pext,pbaru et pbarv : arguments d'entree  pour le s-pg
24  c=======================================================================  
25      ! =======================================================================
26    
27    
28  c   Arguments:  
29  c   ----------    ! Arguments:
30        INTEGER iq,nt    ! ----------
31        REAL, intent(in):: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )    INTEGER iq, nt
32        REAL masse(iip1,jjp1,llm)    REAL, INTENT (IN) :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
33        REAL q( iip1,jjp1,llm,0:9)    REAL masse(iip1, jjp1, llm)
34        REAL w( ip1jmp1,llm )    REAL q(iip1, jjp1, llm, 0:9)
35        integer ordre,ilim    REAL w(ip1jmp1, llm)
36      INTEGER ordre, ilim
37  c   Local:  
38  c   ------    ! Local:
39        LOGICAL limit    ! ------
40        real zq(iip1,jjp1,llm)    LOGICAL limit
41        REAL sm ( iip1,jjp1, llm )    REAL zq(iip1, jjp1, llm)
42        REAL s0( iip1,jjp1,llm ),  sx( iip1,jjp1,llm )    REAL sm(iip1, jjp1, llm)
43        REAL sy( iip1,jjp1,llm ),  sz( iip1,jjp1,llm )    REAL s0(iip1, jjp1, llm), sx(iip1, jjp1, llm)
44        REAL sxx( iip1,jjp1,llm)    REAL sy(iip1, jjp1, llm), sz(iip1, jjp1, llm)
45        REAL sxy( iip1,jjp1,llm)    REAL sxx(iip1, jjp1, llm)
46        REAL sxz( iip1,jjp1,llm)    REAL sxy(iip1, jjp1, llm)
47        REAL syy( iip1,jjp1,llm )    REAL sxz(iip1, jjp1, llm)
48        REAL syz( iip1,jjp1,llm )    REAL syy(iip1, jjp1, llm)
49        REAL szz( iip1,jjp1,llm ),zz    REAL syz(iip1, jjp1, llm)
50        INTEGER i,j,l,indice    REAL szz(iip1, jjp1, llm), zz
51        real sxn(iip1),sxs(iip1)    INTEGER i, j, l, indice
52      REAL sxn(iip1), sxs(iip1)
53        real sinlon(iip1),sinlondlon(iip1)  
54        real coslon(iip1),coslondlon(iip1)    REAL sinlon(iip1), sinlondlon(iip1)
55        real qmin,qmax    REAL coslon(iip1), coslondlon(iip1)
56        save qmin,qmax    REAL qmin, qmax
57        save sinlon,coslon,sinlondlon,coslondlon    SAVE qmin, qmax
58        real dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps    SAVE sinlon, coslon, sinlondlon, coslondlon
59        real masn,mass    REAL dyn1, dyn2, dys1, dys2, qpn, qps, dqzpn, dqzps
60  c    REAL masn, mass
61        REAL      SSUM  
62        integer ismax,ismin    REAL ssum
63        EXTERNAL  SSUM, convflu,ismin,ismax    INTEGER ismax, ismin
64        logical first    EXTERNAL ssum, convflu, ismin, ismax
65        save first    LOGICAL first
66        EXTERNAL advxp,advyp,advzp    SAVE first
67      EXTERNAL advxp, advyp, advzp
68    
69        data first/.true./  
70        data qmin,qmax/-1.e33,1.e33/    DATA first/.TRUE./
71      DATA qmin, qmax/ -1.E33, 1.E33/
72    
73  c==========================================================================  
74  c==========================================================================    ! ==========================================================================
75  c     MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt    ! ==========================================================================
76  c==========================================================================    ! MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
77  c==========================================================================    ! ==========================================================================
78        REAL dt    ! ==========================================================================
79  c==========================================================================    REAL dt
80        limit = .TRUE.    ! ==========================================================================
81      limit = .TRUE.
82        if(first) then  
83           print*,'SCHEMA PRATHER'    IF (first) THEN
84           first=.false.      PRINT *, 'SCHEMA PRATHER'
85           do i=2,iip1      first = .FALSE.
86              coslon(i)=cos(rlonv(i))      DO i = 2, iip1
87              sinlon(i)=sin(rlonv(i))        coslon(i) = cos(rlonv(i))
88              coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi        sinlon(i) = sin(rlonv(i))
89              sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi        coslondlon(i) = coslon(i)*(rlonu(i)-rlonu(i-1))/pi
90           enddo        sinlondlon(i) = sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
91           coslon(1)=coslon(iip1)      END DO
92           coslondlon(1)=coslondlon(iip1)      coslon(1) = coslon(iip1)
93           sinlon(1)=sinlon(iip1)      coslondlon(1) = coslondlon(iip1)
94           sinlondlon(1)=sinlondlon(iip1)      sinlon(1) = sinlon(iip1)
95        sinlondlon(1) = sinlondlon(iip1)
96          DO l = 1,llm  
97          DO j = 1,jjp1      DO l = 1, llm
98          DO i = 1,iip1        DO j = 1, jjp1
99          q( i,j,l,1 )=0.          DO i = 1, iip1
100          q( i,j,l,2)=0.            q(i, j, l, 1) = 0.
101          q( i,j,l,3)=0.            q(i, j, l, 2) = 0.
102          q( i,j,l,4)=0.            q(i, j, l, 3) = 0.
103          q( i,j,l,5)=0.            q(i, j, l, 4) = 0.
104          q( i,j,l,6)=0.            q(i, j, l, 5) = 0.
105          q( i,j,l,7)=0.            q(i, j, l, 6) = 0.
106          q( i,j,l,8)=0.            q(i, j, l, 7) = 0.
107          q( i,j,l,9)=0.            q(i, j, l, 8) = 0.
108          ENDDO            q(i, j, l, 9) = 0.
109          ENDDO          END DO
110          ENDDO        END DO
111        endif      END DO
112  c   Fin modif Fred    END IF
113      ! Fin modif Fred
114  c *** On calcule la masse d'air en kg  
115      ! *** On calcule la masse d'air en kg
116         DO l = 1,llm  
117          DO j = 1,jjp1    DO l = 1, llm
118           DO i = 1,iip1      DO j = 1, jjp1
119           sm( i,j,llm+1-l ) =masse(i,j,l)        DO i = 1, iip1
120           ENDDO          sm(i, j, llm+1-l) = masse(i, j, l)
121          ENDDO        END DO
122         ENDDO      END DO
123      END DO
124  c *** q contient les qqtes de traceur avant l'advection  
125      ! *** q contient les qqtes de traceur avant l'advection
126  c *** Affectation des tableaux S a partir de Q  
127      ! *** Affectation des tableaux S a partir de Q
128         DO l = 1,llm  
129          DO j = 1,jjp1    DO l = 1, llm
130           DO i = 1,iip1      DO j = 1, jjp1
131         s0( i,j,l) = q ( i,j,llm+1-l,0 )*sm(i,j,l)        DO i = 1, iip1
132         sx( i,j,l) = q( i,j,llm+1-l,1 )*sm(i,j,l)          s0(i, j, l) = q(i, j, llm+1-l, 0)*sm(i, j, l)
133         sy( i,j,l) = q( i,j,llm+1-l,2)*sm(i,j,l)          sx(i, j, l) = q(i, j, llm+1-l, 1)*sm(i, j, l)
134         sz( i,j,l) = q( i,j,llm+1-l,3)*sm(i,j,l)          sy(i, j, l) = q(i, j, llm+1-l, 2)*sm(i, j, l)
135         sxx( i,j,l) = q( i,j,llm+1-l,4)*sm(i,j,l)          sz(i, j, l) = q(i, j, llm+1-l, 3)*sm(i, j, l)
136         sxy( i,j,l) = q( i,j,llm+1-l,5)*sm(i,j,l)          sxx(i, j, l) = q(i, j, llm+1-l, 4)*sm(i, j, l)
137         sxz( i,j,l) = q( i,j,llm+1-l,6)*sm(i,j,l)          sxy(i, j, l) = q(i, j, llm+1-l, 5)*sm(i, j, l)
138         syy( i,j,l) = q( i,j,llm+1-l,7)*sm(i,j,l)          sxz(i, j, l) = q(i, j, llm+1-l, 6)*sm(i, j, l)
139         syz( i,j,l) = q( i,j,llm+1-l,8)*sm(i,j,l)          syy(i, j, l) = q(i, j, llm+1-l, 7)*sm(i, j, l)
140         szz( i,j,l) = q( i,j,llm+1-l,9)*sm(i,j,l)          syz(i, j, l) = q(i, j, llm+1-l, 8)*sm(i, j, l)
141           ENDDO          szz(i, j, l) = q(i, j, llm+1-l, 9)*sm(i, j, l)
142          ENDDO        END DO
143         ENDDO      END DO
144  c *** Appel des subroutines d'advection en X, en Y et en Z    END DO
145  c *** Advection avec "time-splitting"    ! *** Appel des subroutines d'advection en X, en Y et en Z
146            ! *** Advection avec "time-splitting"
147  c-----------------------------------------------------------  
148         do indice =1,nt    ! -----------------------------------------------------------
149         call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz    DO indice = 1, nt
150       .             ,sxx,sxy,sxz,syy,syz,szz,1 )      CALL advxp(limit, 0.5*dt, pbaru, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
151          end do        syz, szz, 1)
152          do l=1,llm    END DO
153          do i=1,iip1    DO l = 1, llm
154          sy(i,1,l)=0.      DO i = 1, iip1
155          sy(i,jjp1,l)=0.        sy(i, 1, l) = 0.
156          enddo        sy(i, jjp1, l) = 0.
157          enddo      END DO
158  c---------------------------------------------------------    END DO
159         call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz    ! ---------------------------------------------------------
160       .             ,sxx,sxy,sxz,syy,syz,szz,1 )    CALL advyp(limit, .5*dt*nt, pbarv, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
161  c---------------------------------------------------------      syz, szz, 1)
162      ! ---------------------------------------------------------
163  c---------------------------------------------------------  
164         do j=1,jjp1    ! ---------------------------------------------------------
165            do i=1,iip1    DO j = 1, jjp1
166               sz(i,j,1)=0.      DO i = 1, iip1
167               sz(i,j,llm)=0.        sz(i, j, 1) = 0.
168               sxz(i,j,1)=0.        sz(i, j, llm) = 0.
169               sxz(i,j,llm)=0.        sxz(i, j, 1) = 0.
170               syz(i,j,1)=0.        sxz(i, j, llm) = 0.
171               syz(i,j,llm)=0.        syz(i, j, 1) = 0.
172               szz(i,j,1)=0.        syz(i, j, llm) = 0.
173               szz(i,j,llm)=0.        szz(i, j, 1) = 0.
174            enddo        szz(i, j, llm) = 0.
175         enddo      END DO
176         call advzp( limit,dt*nt,w,sm,s0,sx,sy,sz    END DO
177       .             ,sxx,sxy,sxz,syy,syz,szz,1 )    CALL advzp(limit, dt*nt, w, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, syz, &
178          do l=1,llm      szz, 1)
179          do i=1,iip1    DO l = 1, llm
180          sy(i,1,l)=0.      DO i = 1, iip1
181          sy(i,jjp1,l)=0.        sy(i, 1, l) = 0.
182          enddo        sy(i, jjp1, l) = 0.
183          enddo      END DO
184      END DO
185  c---------------------------------------------------------  
186      ! ---------------------------------------------------------
187  c---------------------------------------------------------  
188         call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz    ! ---------------------------------------------------------
189       .             ,sxx,sxy,sxz,syy,syz,szz,1 )    CALL advyp(limit, .5*dt*nt, pbarv, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
190  c---------------------------------------------------------      syz, szz, 1)
191         DO l = 1,llm    ! ---------------------------------------------------------
192          DO j = 1,jjp1    DO l = 1, llm
193               s0( iip1,j,l)=s0( 1,j,l )      DO j = 1, jjp1
194               sx( iip1,j,l)=sx( 1,j,l )        s0(iip1, j, l) = s0(1, j, l)
195               sy( iip1,j,l)=sy( 1,j,l )        sx(iip1, j, l) = sx(1, j, l)
196               sz( iip1,j,l)=sz( 1,j,l )        sy(iip1, j, l) = sy(1, j, l)
197               sxx( iip1,j,l)=sxx( 1,j,l )        sz(iip1, j, l) = sz(1, j, l)
198               sxy( iip1,j,l)=sxy( 1,j,l)        sxx(iip1, j, l) = sxx(1, j, l)
199               sxz( iip1,j,l)=sxz( 1,j,l )        sxy(iip1, j, l) = sxy(1, j, l)
200               syy( iip1,j,l)=syy( 1,j,l )        sxz(iip1, j, l) = sxz(1, j, l)
201               syz( iip1,j,l)=syz( 1,j,l)        syy(iip1, j, l) = syy(1, j, l)
202               szz( iip1,j,l)=szz( 1,j,l )        syz(iip1, j, l) = syz(1, j, l)
203          ENDDO        szz(iip1, j, l) = szz(1, j, l)
204         ENDDO      END DO
205         do indice=1,nt    END DO
206         call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz    DO indice = 1, nt
207       .             ,sxx,sxy,sxz,syy,syz,szz,1 )      CALL advxp(limit, 0.5*dt, pbaru, sm, s0, sx, sy, sz, sxx, sxy, sxz, syy, &
208          end do        syz, szz, 1)
209  c---------------------------------------------------------    END DO
210  c---------------------------------------------------------    ! ---------------------------------------------------------
211  c ***   On repasse les S dans la variable qpr    ! ---------------------------------------------------------
212  c ***   On repasse les S dans la variable q directement 14/10/94    ! ***   On repasse les S dans la variable qpr
213      ! ***   On repasse les S dans la variable q directement 14/10/94
214         DO  l = 1,llm  
215          DO  j = 1,jjp1    DO l = 1, llm
216           DO  i = 1,iip1      DO j = 1, jjp1
217        q( i,j,llm+1-l,0 )=s0( i,j,l )/sm(i,j,l)        DO i = 1, iip1
218        q( i,j,llm+1-l,1 ) = sx( i,j,l )/sm(i,j,l)          q(i, j, llm+1-l, 0) = s0(i, j, l)/sm(i, j, l)
219        q( i,j,llm+1-l,2 ) = sy( i,j,l )/sm(i,j,l)          q(i, j, llm+1-l, 1) = sx(i, j, l)/sm(i, j, l)
220        q( i,j,llm+1-l,3 ) = sz( i,j,l )/sm(i,j,l)          q(i, j, llm+1-l, 2) = sy(i, j, l)/sm(i, j, l)
221        q( i,j,llm+1-l,4 ) = sxx( i,j,l )/sm(i,j,l)          q(i, j, llm+1-l, 3) = sz(i, j, l)/sm(i, j, l)
222        q( i,j,llm+1-l,5 ) = sxy( i,j,l )/sm(i,j,l)          q(i, j, llm+1-l, 4) = sxx(i, j, l)/sm(i, j, l)
223        q( i,j,llm+1-l,6 ) = sxz( i,j,l )/sm(i,j,l)          q(i, j, llm+1-l, 5) = sxy(i, j, l)/sm(i, j, l)
224        q( i,j,llm+1-l,7 ) = syy( i,j,l )/sm(i,j,l)          q(i, j, llm+1-l, 6) = sxz(i, j, l)/sm(i, j, l)
225        q( i,j,llm+1-l,8 ) = syz( i,j,l )/sm(i,j,l)          q(i, j, llm+1-l, 7) = syy(i, j, l)/sm(i, j, l)
226        q( i,j,llm+1-l,9 ) = szz( i,j,l )/sm(i,j,l)          q(i, j, llm+1-l, 8) = syz(i, j, l)/sm(i, j, l)
227        ENDDO          q(i, j, llm+1-l, 9) = szz(i, j, l)/sm(i, j, l)
228        ENDDO        END DO
229        ENDDO      END DO
230      END DO
231  c---------------------------------------------------------  
232  c      go to  777    ! ---------------------------------------------------------
233  c   filtrages aux poles    ! go to  777
234      ! filtrages aux poles
235  c Traitements specifiques au pole  
236      ! Traitements specifiques au pole
237  c   filtrages aux poles  
238           DO l=1,llm    ! filtrages aux poles
239  c   filtrages aux poles    DO l = 1, llm
240           masn=ssum(iim,sm(1,1,l),1)      ! filtrages aux poles
241           mass=ssum(iim,sm(1,jjp1,l),1)      masn = ssum(iim, sm(1,1,l), 1)
242           qpn=ssum(iim,s0(1,1,l),1)/masn      mass = ssum(iim, sm(1,jjp1,l), 1)
243           qps=ssum(iim,s0(1,jjp1,l),1)/mass      qpn = ssum(iim, s0(1,1,l), 1)/masn
244           dqzpn=ssum(iim,sz(1,1,l),1)/masn      qps = ssum(iim, s0(1,jjp1,l), 1)/mass
245           dqzps=ssum(iim,sz(1,jjp1,l),1)/mass      dqzpn = ssum(iim, sz(1,1,l), 1)/masn
246           do i=1,iip1      dqzps = ssum(iim, sz(1,jjp1,l), 1)/mass
247            q( i,1,llm+1-l,3)=dqzpn      DO i = 1, iip1
248            q( i,jjp1,llm+1-l,3)=dqzps        q(i, 1, llm+1-l, 3) = dqzpn
249            q( i,1,llm+1-l,0)=qpn        q(i, jjp1, llm+1-l, 3) = dqzps
250            q( i,jjp1,llm+1-l,0)=qps        q(i, 1, llm+1-l, 0) = qpn
251           enddo        q(i, jjp1, llm+1-l, 0) = qps
252             dyn1=0.      END DO
253             dys1=0.      dyn1 = 0.
254             dyn2=0.      dys1 = 0.
255             dys2=0.      dyn2 = 0.
256          do i=1,iim      dys2 = 0.
257          zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)      DO i = 1, iim
258          dyn1=dyn1+sinlondlon(i)*zz        zz = s0(i, 2, l)/sm(i, 2, l) - q(i, 1, llm+1-l, 0)
259          dyn2=dyn2+coslondlon(i)*zz        dyn1 = dyn1 + sinlondlon(i)*zz
260          zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)        dyn2 = dyn2 + coslondlon(i)*zz
261          dys1=dys1+sinlondlon(i)*zz        zz = q(i, jjp1, llm+1-l, 0) - s0(i, jjm, l)/sm(i, jjm, l)
262          dys2=dys2+coslondlon(i)*zz        dys1 = dys1 + sinlondlon(i)*zz
263          enddo        dys2 = dys2 + coslondlon(i)*zz
264           do i=1,iim      END DO
265           q(i,1,llm+1-l,2)=      DO i = 1, iim
266       $   (sinlon(i)*dyn1+coslon(i)*dyn2)/2.        q(i, 1, llm+1-l, 2) = (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
267           q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)        q(i, 1, llm+1-l, 0) = q(i, 1, llm+1-l, 0) + q(i, 1, llm+1-l, 2)
268       $          +q(i,1,llm+1-l,2)        q(i, jjp1, llm+1-l, 2) = (sinlon(i)*dys1+coslon(i)*dys2)/2.
269           q(i,jjp1,llm+1-l,2)=        q(i, jjp1, llm+1-l, 0) = q(i, jjp1, llm+1-l, 0) - &
270       $   (sinlon(i)*dys1+coslon(i)*dys2)/2.          q(i, jjp1, llm+1-l, 2)
271           q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)      END DO
272       $      -q(i,jjp1,llm+1-l,2)      q(iip1, 1, llm+1-l, 0) = q(1, 1, llm+1-l, 0)
273           enddo      q(iip1, jjp1, llm+1-l, 0) = q(1, jjp1, llm+1-l, 0)
274        q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)      DO i = 1, iim
275        q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)        sxn(i) = q(i+1, 1, llm+1-l, 0) - q(i, 1, llm+1-l, 0)
276        do i=1,iim        sxs(i) = q(i+1, jjp1, llm+1-l, 0) - q(i, jjp1, llm+1-l, 0)
277        sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)      END DO
278        sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)      sxn(iip1) = sxn(1)
279        enddo      sxs(iip1) = sxs(1)
280        sxn(iip1)=sxn(1)      DO i = 1, iim
281        sxs(iip1)=sxs(1)        q(i+1, 1, llm+1-l, 1) = 0.25*(sxn(i)+sxn(i+1))
282        do i=1,iim        q(i+1, jjp1, llm+1-l, 1) = 0.25*(sxs(i)+sxs(i+1))
283        q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))      END DO
284        q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))      q(1, 1, llm+1-l, 1) = q(iip1, 1, llm+1-l, 1)
285        q(1, jjp1, llm+1-l, 1) = q(iip1, jjp1, llm+1-l, 1)
286      END DO
287      DO l = 1, llm
288        DO i = 1, iim
289          q(i, 1, llm+1-l, 4) = 0.
290          q(i, jjp1, llm+1-l, 4) = 0.
291          q(i, 1, llm+1-l, 5) = 0.
292          q(i, jjp1, llm+1-l, 5) = 0.
293          q(i, 1, llm+1-l, 6) = 0.
294          q(i, jjp1, llm+1-l, 6) = 0.
295          q(i, 1, llm+1-l, 7) = 0.
296          q(i, jjp1, llm+1-l, 7) = 0.
297          q(i, 1, llm+1-l, 8) = 0.
298          q(i, jjp1, llm+1-l, 8) = 0.
299          q(i, 1, llm+1-l, 9) = 0.
300          q(i, jjp1, llm+1-l, 9) = 0.
301        END DO
302      END DO
303    
304      ! bouclage en longitude
305      DO l = 1, llm
306        DO j = 1, jjp1
307          q(iip1, j, l, 0) = q(1, j, l, 0)
308          q(iip1, j, llm+1-l, 0) = q(1, j, llm+1-l, 0)
309          q(iip1, j, llm+1-l, 1) = q(1, j, llm+1-l, 1)
310          q(iip1, j, llm+1-l, 2) = q(1, j, llm+1-l, 2)
311          q(iip1, j, llm+1-l, 3) = q(1, j, llm+1-l, 3)
312          q(iip1, j, llm+1-l, 4) = q(1, j, llm+1-l, 4)
313          q(iip1, j, llm+1-l, 5) = q(1, j, llm+1-l, 5)
314          q(iip1, j, llm+1-l, 6) = q(1, j, llm+1-l, 6)
315          q(iip1, j, llm+1-l, 7) = q(1, j, llm+1-l, 7)
316          q(iip1, j, llm+1-l, 8) = q(1, j, llm+1-l, 8)
317          q(iip1, j, llm+1-l, 9) = q(1, j, llm+1-l, 9)
318        END DO
319      END DO
320      DO l = 1, llm
321        DO j = 2, jjm
322          DO i = 1, iip1
323            IF (q(i,j,l,0)<0.) THEN
324              PRINT *, '------------ BIP-----------'
325              PRINT *, 'S0(', i, j, l, ')=', q(i, j, l, 0), q(i, j-1, l, 0)
326              PRINT *, 'SX(', i, j, l, ')=', q(i, j, l, 1)
327              PRINT *, 'SY(', i, j, l, ')=', q(i, j, l, 2), q(i, j-1, l, 2)
328              PRINT *, 'SZ(', i, j, l, ')=', q(i, j, l, 3)
329              q(i, j, l, 0) = 0.
330            END IF
331          END DO
332        END DO
333        DO j = 1, jjp1, jjm
334          DO i = 1, iip1
335            IF (q(i,j,l,0)<0.) THEN
336              PRINT *, '------------ BIP 2-----------'
337              PRINT *, 'S0(', i, j, l, ')=', q(i, j, l, 0)
338              PRINT *, 'SX(', i, j, l, ')=', q(i, j, l, 1)
339              PRINT *, 'SY(', i, j, l, ')=', q(i, j, l, 2)
340              PRINT *, 'SZ(', i, j, l, ')=', q(i, j, l, 3)
341    
342              q(i, j, l, 0) = 0.
343              ! STOP
344            END IF
345        END DO        END DO
346        q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)      END DO
347        q(1,jjp1,llm+1-l,1)=    END DO
348       $   q(iip1,jjp1,llm+1-l,1)    RETURN
349          enddo  END SUBROUTINE prather
          do l=1,llm  
            do i=1,iim  
             q( i,1,llm+1-l,4)=0.  
             q( i,jjp1,llm+1-l,4)=0.  
             q( i,1,llm+1-l,5)=0.  
             q( i,jjp1,llm+1-l,5)=0.  
             q( i,1,llm+1-l,6)=0.  
             q( i,jjp1,llm+1-l,6)=0.  
             q( i,1,llm+1-l,7)=0.  
             q( i,jjp1,llm+1-l,7)=0.  
             q( i,1,llm+1-l,8)=0.  
             q( i,jjp1,llm+1-l,8)=0.  
             q( i,1,llm+1-l,9)=0.  
             q( i,jjp1,llm+1-l,9)=0.  
           enddo  
          ENDDO  
   
 777      continue  
 c  
 c   bouclage en longitude  
       do l=1,llm  
       do j=1,jjp1  
       q(iip1,j,l,0)=q(1,j,l,0)  
       q(iip1,j,llm+1-l,0)=q(1,j,llm+1-l,0)  
       q(iip1,j,llm+1-l,1)=q(1,j,llm+1-l,1)  
       q(iip1,j,llm+1-l,2)=q(1,j,llm+1-l,2)  
       q(iip1,j,llm+1-l,3)=q(1,j,llm+1-l,3)  
       q(iip1,j,llm+1-l,4)=q(1,j,llm+1-l,4)  
       q(iip1,j,llm+1-l,5)=q(1,j,llm+1-l,5)  
       q(iip1,j,llm+1-l,6)=q(1,j,llm+1-l,6)  
       q(iip1,j,llm+1-l,7)=q(1,j,llm+1-l,7)  
       q(iip1,j,llm+1-l,8)=q(1,j,llm+1-l,8)  
       q(iip1,j,llm+1-l,9)=q(1,j,llm+1-l,9)  
       enddo  
       enddo  
         DO l = 1,llm  
          DO j = 2,jjm  
            DO i = 1,iip1  
          IF (q(i,j,l,0).lt.0.)  THEN  
          PRINT*,'------------ BIP-----------'  
          PRINT*,'S0(',i,j,l,')=',q(i,j,l,0),  
      $          q(i,j-1,l,0)  
          PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)  
          PRINT*,'SY(',i,j,l,')=',q(i,j,l,2),  
      $   q(i,j-1,l,2)    
          PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)  
                      q(i,j,l,0)=0.  
                ENDIF  
            ENDDO  
          ENDDO  
          do j=1,jjp1,jjm  
          do i=1,iip1  
                IF (q(i,j,l,0).lt.0.)  THEN  
                PRINT*,'------------ BIP 2-----------'  
          PRINT*,'S0(',i,j,l,')=',q(i,j,l,0)  
          PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)  
          PRINT*,'SY(',i,j,l,')=',q(i,j,l,2)  
          PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)  
   
                      q(i,j,l,0)=0.  
 c                  STOP  
                ENDIF  
          enddo  
          enddo  
         ENDDO  
       RETURN  
       END  
Legend:
Removed from v.80  
changed lines
  Added in v.81
  ViewVC Help
Powered by ViewVC 1.1.21