/[lmdze]/trunk/Sources/dyn3d/advn.f
ViewVC logotype

Diff of /trunk/Sources/dyn3d/advn.f

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

trunk/libf/dyn3d/advn.f revision 31 by guez, Thu Apr 1 14:59:19 2010 UTC trunk/Sources/dyn3d/advn.f revision 157 by guez, Mon Jul 20 16:01:49 2015 UTC
# Line 1  Line 1 
1  !  
2  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advn.F,v 1.1.1.1 2004/05/19 12:53:06 lmdzadmin Exp $  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/advn.F,v 1.1.1.1 2004/05/19
3  !  ! 12:53:06 lmdzadmin Exp $
4        SUBROUTINE advn(q,masse,w,pbaru,pbarv,pdt,mode)  
5  c  SUBROUTINE advn(q, masse, w, pbaru, pbarv, pdt, mode)
6  c     Auteur : F. Hourdin  
7  c    ! Auteur : F. Hourdin
8  c    ********************************************************************  
9  c     Shema  d'advection " pseudo amont " .    ! ********************************************************************
10  c    ********************************************************************    ! Shema  d'advection " pseudo amont " .
11  c     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....    ! ********************************************************************
12  c    ! q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
13  c   pbaru,pbarv,w flux de masse en u ,v ,w  
14  c   pdt pas de temps    ! pbaru,pbarv,w flux de masse en u ,v ,w
15  c    ! pdt pas de temps
16  c   --------------------------------------------------------------------  
17        use dimens_m    ! --------------------------------------------------------------------
18        use paramet_m    USE dimens_m
19        use comconst    USE paramet_m
20        use comvert    USE comconst
21        use logic    USE disvert_m
22        use comgeom    USE conf_gcm_m
23        use iniprint    USE comgeom
24        IMPLICIT NONE    IMPLICIT NONE
25  c  
26    
27  c  
28  c   Arguments:    ! Arguments:
29  c   ----------    ! ----------
30        integer mode    INTEGER mode
31        real masse(ip1jmp1,llm)    REAL masse(ip1jmp1, llm)
32        REAL, intent(in):: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)    REAL, INTENT (IN) :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
33        REAL q(ip1jmp1,llm)    REAL q(ip1jmp1, llm)
34        REAL w(ip1jmp1,llm),pdt    REAL w(ip1jmp1, llm), pdt
35  c  
36  c      Local    ! Local
37  c   ---------    ! ---------
38  c  
39        INTEGER i,ij,l,j,ii    INTEGER ij, l
40        integer ijlqmin,iqmin,jqmin,lqmin    REAL zm(ip1jmp1, llm)
41        integer ismin    REAL mu(ip1jmp1, llm)
42  c    REAL mv(ip1jm, llm)
43        real zm(ip1jmp1,llm),newmasse    REAL mw(ip1jmp1, llm+1)
44        real mu(ip1jmp1,llm)    REAL zq(ip1jmp1, llm), qpn, qps
45        real mv(ip1jm,llm)    REAL zqg(ip1jmp1, llm), zqd(ip1jmp1, llm)
46        real mw(ip1jmp1,llm+1)    REAL zqs(ip1jmp1, llm), zqn(ip1jmp1, llm)
47        real zq(ip1jmp1,llm),zz,qpn,qps    REAL zqh(ip1jmp1, llm), zqb(ip1jmp1, llm)
48        real zqg(ip1jmp1,llm),zqd(ip1jmp1,llm)    REAL ssum
49        real zqs(ip1jmp1,llm),zqn(ip1jmp1,llm)    REAL zzpbar, zzw
50        real zqh(ip1jmp1,llm),zqb(ip1jmp1,llm)  
51        real temps0,temps1,temps2,temps3    zzpbar = 0.5*pdt
52        real ztemps1,ztemps2,ztemps3,ssum    zzw = pdt
53        logical testcpu  
54        save testcpu    DO l = 1, llm
55        save temps1,temps2,temps3      DO ij = iip2, ip1jm
56        real zzpbar,zzw        mu(ij, l) = pbaru(ij, l)*zzpbar
57        END DO
58        real qmin,qmax      DO ij = 1, ip1jm
59        data qmin,qmax/0.,1./        mv(ij, l) = pbarv(ij, l)*zzpbar
60        data testcpu/.false./      END DO
61        data temps1,temps2,temps3/0.,0.,0./      DO ij = 1, ip1jmp1
62          mw(ij, l) = w(ij, l)*zzw
63        zzpbar = 0.5 * pdt      END DO
64        zzw    = pdt    END DO
65    
66        DO l=1,llm    DO ij = 1, ip1jmp1
67          DO ij = iip2,ip1jm      mw(ij, llm+1) = 0.
68              mu(ij,l)=pbaru(ij,l) * zzpbar    END DO
69           ENDDO  
70           DO ij=1,ip1jm    DO l = 1, llm
71              mv(ij,l)=pbarv(ij,l) * zzpbar      qpn = 0.
72           ENDDO      qps = 0.
73           DO ij=1,ip1jmp1      DO ij = 1, iim
74              mw(ij,l)=w(ij,l) * zzw        qpn = qpn + q(ij, l)*masse(ij, l)
75           ENDDO        qps = qps + q(ip1jm+ij, l)*masse(ip1jm+ij, l)
76        ENDDO      END DO
77        qpn = qpn/ssum(iim, masse(1,l), 1)
78        DO ij=1,ip1jmp1      qps = qps/ssum(iim, masse(ip1jm+1,l), 1)
79           mw(ij,llm+1)=0.      DO ij = 1, iip1
80        ENDDO        q(ij, l) = qpn
81          q(ip1jm+ij, l) = qps
82        do l=1,llm      END DO
83           qpn=0.    END DO
84           qps=0.  
85           do ij=1,iim    DO ij = 1, ip1jmp1
86              qpn=qpn+q(ij,l)*masse(ij,l)      mw(ij, llm+1) = 0.
87              qps=qps+q(ip1jm+ij,l)*masse(ip1jm+ij,l)    END DO
88           enddo    DO l = 1, llm
89           qpn=qpn/ssum(iim,masse(1,l),1)      DO ij = 1, ip1jmp1
90           qps=qps/ssum(iim,masse(ip1jm+1,l),1)        zq(ij, l) = q(ij, l)
91           do ij=1,iip1        zm(ij, l) = masse(ij, l)
92              q(ij,l)=qpn      END DO
93              q(ip1jm+ij,l)=qps    END DO
94           enddo  
95        enddo    ! call minmaxq(zq,qmin,qmax,'avant vlx     ')
96      CALL advnqx(zq, zqg, zqd)
97        do ij=1,ip1jmp1    CALL advnx(zq, zqg, zqd, zm, mu, mode)
98           mw(ij,llm+1)=0.    CALL advnqy(zq, zqs, zqn)
99        enddo    CALL advny(zq, zqs, zqn, zm, mv)
100        do l=1,llm    CALL advnqz(zq, zqh, zqb)
101           do ij=1,ip1jmp1    CALL advnz(zq, zqh, zqb, zm, mw)
102              zq(ij,l)=q(ij,l)    ! call vlz(zq,0.,zm,mw)
103              zm(ij,l)=masse(ij,l)    CALL advnqy(zq, zqs, zqn)
104           enddo    CALL advny(zq, zqs, zqn, zm, mv)
105        enddo    CALL advnqx(zq, zqg, zqd)
106      CALL advnx(zq, zqg, zqd, zm, mu, mode)
107  c     call minmaxq(zq,qmin,qmax,'avant vlx     ')    ! call minmaxq(zq,qmin,qmax,'apres vlx     ')
108        call advnqx(zq,zqg,zqd)  
109        call advnx(zq,zqg,zqd,zm,mu,mode)    DO l = 1, llm
110        call advnqy(zq,zqs,zqn)      DO ij = 1, ip1jmp1
111        call advny(zq,zqs,zqn,zm,mv)        q(ij, l) = zq(ij, l)
112        call advnqz(zq,zqh,zqb)      END DO
113        call advnz(zq,zqh,zqb,zm,mw)      DO ij = 1, ip1jm + 1, iip1
114  c     call vlz(zq,0.,zm,mw)        q(ij+iim, l) = q(ij, l)
115        call advnqy(zq,zqs,zqn)      END DO
116        call advny(zq,zqs,zqn,zm,mv)    END DO
117        call advnqx(zq,zqg,zqd)  
118        call advnx(zq,zqg,zqd,zm,mu,mode)    RETURN
119  c     call minmaxq(zq,qmin,qmax,'apres vlx     ')  END SUBROUTINE advn
120    
121        do l=1,llm  SUBROUTINE advnqx(q, qg, qd)
122           do ij=1,ip1jmp1  
123             q(ij,l)=zq(ij,l)    ! Auteurs:   Calcul des valeurs de q aux point u.
124           enddo  
125           do ij=1,ip1jm+1,iip1    ! --------------------------------------------------------------------
126              q(ij+iim,l)=q(ij,l)    USE dimens_m
127           enddo    USE paramet_m
128        enddo    USE conf_gcm_m
129      IMPLICIT NONE
130        RETURN  
131        END  
132    
133        SUBROUTINE advnqx(q,qg,qd)    ! Arguments:
134  c    ! ----------
135  c     Auteurs:   Calcul des valeurs de q aux point u.    REAL q(ip1jmp1, llm), qg(ip1jmp1, llm), qd(ip1jmp1, llm)
136  c  
137  c   --------------------------------------------------------------------    ! Local
138        use dimens_m    ! ---------
139        use paramet_m  
140        use iniprint    INTEGER ij, l
141        IMPLICIT NONE  
142  c    REAL dxqu(ip1jmp1), zqu(ip1jmp1)
143  c    REAL zqmax(ip1jmp1), zqmin(ip1jmp1)
144  c    LOGICAL extremum(ip1jmp1)
145  c   Arguments:  
146  c   ----------    INTEGER mode
147        real q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)    SAVE mode
148  c    DATA mode/1/
149  c      Local  
150  c   ---------    ! calcul des pentes en u:
151  c    ! -----------------------
152        INTEGER ij,l    IF (mode==0) THEN
153  c      DO l = 1, llm
154        real dxqu(ip1jmp1),zqu(ip1jmp1)        DO ij = 1, ip1jm
155        real zqmax(ip1jmp1),zqmin(ip1jmp1)          qd(ij, l) = q(ij, l)
156        logical extremum(ip1jmp1)          qg(ij, l) = q(ij, l)
157          END DO
158        integer mode      END DO
159        save mode    ELSE
160        data mode/1/      DO l = 1, llm
161          DO ij = iip2, ip1jm - 1
162  c   calcul des pentes en u:          dxqu(ij) = q(ij+1, l) - q(ij, l)
163  c   -----------------------          zqu(ij) = 0.5*(q(ij+1,l)+q(ij,l))
164        if (mode.eq.0) then        END DO
165           do l=1,llm        DO ij = iip1 + iip1, ip1jm, iip1
166              do ij=1,ip1jm          dxqu(ij) = dxqu(ij-iim)
167                 qd(ij,l)=q(ij,l)          zqu(ij) = zqu(ij-iim)
168                 qg(ij,l)=q(ij,l)        END DO
169              enddo        DO ij = iip2, ip1jm - 1
170           enddo          zqu(ij) = zqu(ij) - dxqu(ij+1)/12.
171        else        END DO
172        do l = 1, llm        DO ij = iip1 + iip1, ip1jm, iip1
173           do ij=iip2,ip1jm-1          zqu(ij) = zqu(ij-iim)
174              dxqu(ij)=q(ij+1,l)-q(ij,l)        END DO
175              zqu(ij)=0.5*(q(ij+1,l)+q(ij,l))        DO ij = iip2 + 1, ip1jm
176           enddo          zqu(ij) = zqu(ij) + dxqu(ij-1)/12.
177           do ij=iip1+iip1,ip1jm,iip1        END DO
178              dxqu(ij)=dxqu(ij-iim)        DO ij = iip1 + iip1, ip1jm, iip1
179              zqu(ij)=zqu(ij-iim)          zqu(ij-iim) = zqu(ij)
180           enddo        END DO
181           do ij=iip2,ip1jm-1  
182              zqu(ij)=zqu(ij)-dxqu(ij+1)/12.        ! calcul des valeurs max et min acceptees aux interfaces
183           enddo  
184           do ij=iip1+iip1,ip1jm,iip1        DO ij = iip2, ip1jm - 1
185              zqu(ij)=zqu(ij-iim)          zqmax(ij) = max(q(ij+1,l), q(ij,l))
186           enddo          zqmin(ij) = min(q(ij+1,l), q(ij,l))
187           do ij=iip2+1,ip1jm        END DO
188              zqu(ij)=zqu(ij)+dxqu(ij-1)/12.        DO ij = iip1 + iip1, ip1jm, iip1
189           enddo          zqmax(ij) = zqmax(ij-iim)
190           do ij=iip1+iip1,ip1jm,iip1          zqmin(ij) = zqmin(ij-iim)
191              zqu(ij-iim)=zqu(ij)        END DO
192           enddo        DO ij = iip2 + 1, ip1jm
193            extremum(ij) = dxqu(ij)*dxqu(ij-1) <= 0.
194  c   calcul des valeurs max et min acceptees aux interfaces        END DO
195          DO ij = iip1 + iip1, ip1jm, iip1
196           do ij=iip2,ip1jm-1          extremum(ij-iim) = extremum(ij)
197              zqmax(ij)=max(q(ij+1,l),q(ij,l))        END DO
198              zqmin(ij)=min(q(ij+1,l),q(ij,l))        DO ij = iip2, ip1jm
199           enddo          zqu(ij) = min(max(zqmin(ij),zqu(ij)), zqmax(ij))
200           do ij=iip1+iip1,ip1jm,iip1        END DO
201              zqmax(ij)=zqmax(ij-iim)        DO ij = iip2 + 1, ip1jm
202              zqmin(ij)=zqmin(ij-iim)          IF (extremum(ij)) THEN
203           enddo            qg(ij, l) = q(ij, l)
204           do ij=iip2+1,ip1jm            qd(ij, l) = q(ij, l)
205              extremum(ij)=dxqu(ij)*dxqu(ij-1).le.0.          ELSE
206           enddo            qd(ij, l) = zqu(ij)
207           do ij=iip1+iip1,ip1jm,iip1            qg(ij, l) = zqu(ij-1)
208              extremum(ij-iim)=extremum(ij)          END IF
209           enddo        END DO
210           do ij=iip2,ip1jm        DO ij = iip1 + iip1, ip1jm, iip1
211              zqu(ij)=min(max(zqmin(ij),zqu(ij)),zqmax(ij))          qd(ij-iim, l) = qd(ij, l)
212           enddo          qg(ij-iim, l) = qg(ij, l)
213           do ij=iip2+1,ip1jm        END DO
214              if(extremum(ij)) then  
215                 qg(ij,l)=q(ij,l)        GO TO 8888
216                 qd(ij,l)=q(ij,l)  
217              else        DO ij = iip2 + 1, ip1jm
218                 qd(ij,l)=zqu(ij)          IF (extremum(ij) .AND. .NOT. extremum(ij-1)) qd(ij-1, l) = q(ij, l)
219                 qg(ij,l)=zqu(ij-1)        END DO
220              endif  
221           enddo        DO ij = iip1 + iip1, ip1jm, iip1
222           do ij=iip1+iip1,ip1jm,iip1          qd(ij-iim, l) = qd(ij, l)
223              qd(ij-iim,l)=qd(ij,l)        END DO
224              qg(ij-iim,l)=qg(ij,l)        DO ij = iip2, ip1jm - 1
225           enddo          IF (extremum(ij) .AND. .NOT. extremum(ij+1)) qg(ij+1, l) = q(ij, l)
226          END DO
227           goto 8888  
228          DO ij = iip1 + iip1, ip1jm, iip1
229           do ij=iip2+1,ip1jm          qg(ij, l) = qg(ij-iim, l)
230              if(extremum(ij).and..not.extremum(ij-1))        END DO
231       s         qd(ij-1,l)=q(ij,l)  8888  CONTINUE
232           enddo      END DO
233      END IF
234           do ij=iip1+iip1,ip1jm,iip1    RETURN
235              qd(ij-iim,l)=qd(ij,l)  END SUBROUTINE advnqx
236           enddo  SUBROUTINE advnqy(q, qs, qn)
237           do ij=iip2,ip1jm-1  
238              if (extremum(ij).and..not.extremum(ij+1))    ! Auteurs:   Calcul des valeurs de q aux point v.
239       s         qg(ij+1,l)=q(ij,l)  
240           enddo    ! --------------------------------------------------------------------
241      USE dimens_m
242           do ij=iip1+iip1,ip1jm,iip1    USE paramet_m
243              qg(ij,l)=qg(ij-iim,l)    USE conf_gcm_m
244           enddo    IMPLICIT NONE
245  8888     continue  
246        enddo  
247        endif  
248        RETURN    ! Arguments:
249        END    ! ----------
250        SUBROUTINE advnqy(q,qs,qn)    REAL q(ip1jmp1, llm), qs(ip1jmp1, llm), qn(ip1jmp1, llm)
251  c  
252  c     Auteurs:   Calcul des valeurs de q aux point v.    ! Local
253  c    ! ---------
254  c   --------------------------------------------------------------------  
255        use dimens_m    INTEGER ij, l
256        use paramet_m  
257        use iniprint    REAL dyqv(ip1jm), zqv(ip1jm, llm)
258        IMPLICIT NONE    REAL zqmax(ip1jm), zqmin(ip1jm)
259  c    LOGICAL extremum(ip1jmp1)
260  c  
261  c    INTEGER mode
262  c   Arguments:    SAVE mode
263  c   ----------    DATA mode/1/
264        real q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)  
265  c    IF (mode==0) THEN
266  c      Local      DO l = 1, llm
267  c   ---------        DO ij = 1, ip1jmp1
268  c          qn(ij, l) = q(ij, l)
269        INTEGER ij,l          qs(ij, l) = q(ij, l)
270  c        END DO
271        real dyqv(ip1jm),zqv(ip1jm,llm)      END DO
272        real zqmax(ip1jm),zqmin(ip1jm)    ELSE
273        logical extremum(ip1jmp1)  
274        ! calcul des pentes en u:
275        integer mode      ! -----------------------
276        save mode      DO l = 1, llm
277        data mode/1/        DO ij = 1, ip1jm
278            dyqv(ij) = q(ij, l) - q(ij+iip1, l)
279        if (mode.eq.0) then        END DO
280           do l=1,llm  
281              do ij=1,ip1jmp1        DO ij = iip2, ip1jm - iip1
282                 qn(ij,l)=q(ij,l)          zqv(ij, l) = 0.5*(q(ij+iip1,l)+q(ij,l))
283                 qs(ij,l)=q(ij,l)          zqv(ij, l) = zqv(ij, l) + (dyqv(ij+iip1)-dyqv(ij-iip1))/12.
284              enddo        END DO
285           enddo  
286        else        DO ij = iip2, ip1jm
287            extremum(ij) = dyqv(ij)*dyqv(ij-iip1) <= 0.
288  c   calcul des pentes en u:        END DO
289  c   -----------------------  
290        do l = 1, llm        ! Pas de pentes aux poles
291           do ij=1,ip1jm        DO ij = 1, iip1
292              dyqv(ij)=q(ij,l)-q(ij+iip1,l)          zqv(ij, l) = q(ij, l)
293           enddo          zqv(ip1jm-iip1+ij, l) = q(ip1jm+ij, l)
294            extremum(ij) = .TRUE.
295           do ij=iip2,ip1jm-iip1          extremum(ip1jmp1-iip1+ij) = .TRUE.
296              zqv(ij,l)=0.5*(q(ij+iip1,l)+q(ij,l))        END DO
297              zqv(ij,l)=zqv(ij,l)+(dyqv(ij+iip1)-dyqv(ij-iip1))/12.  
298           enddo        ! calcul des valeurs max et min acceptees aux interfaces
299          DO ij = 1, ip1jm
300           do ij=iip2,ip1jm          zqmax(ij) = max(q(ij+iip1,l), q(ij,l))
301              extremum(ij)=dyqv(ij)*dyqv(ij-iip1).le.0.          zqmin(ij) = min(q(ij+iip1,l), q(ij,l))
302           enddo        END DO
303    
304  c Pas de pentes aux poles        DO ij = 1, ip1jm
305           do ij=1,iip1          zqv(ij, l) = min(max(zqmin(ij),zqv(ij,l)), zqmax(ij))
306              zqv(ij,l)=q(ij,l)        END DO
307              zqv(ip1jm-iip1+ij,l)=q(ip1jm+ij,l)  
308              extremum(ij)=.true.        DO ij = iip2, ip1jm
309              extremum(ip1jmp1-iip1+ij)=.true.          IF (extremum(ij)) THEN
310           enddo            qs(ij, l) = q(ij, l)
311              qn(ij, l) = q(ij, l)
312  c   calcul des valeurs max et min acceptees aux interfaces            ! if (.not.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
313           do ij=1,ip1jm            ! if (.not.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l)
314              zqmax(ij)=max(q(ij+iip1,l),q(ij,l))          ELSE
315              zqmin(ij)=min(q(ij+iip1,l),q(ij,l))            qs(ij, l) = zqv(ij, l)
316           enddo            qn(ij, l) = zqv(ij-iip1, l)
317            END IF
318           do ij=1,ip1jm        END DO
319              zqv(ij,l)=min(max(zqmin(ij),zqv(ij,l)),zqmax(ij))  
320           enddo        DO ij = 1, iip1
321            qs(ij, l) = q(ij, l)
322           do ij=iip2,ip1jm          qn(ij, l) = q(ij, l)
323              if(extremum(ij)) then          qs(ip1jm+ij, l) = q(ip1jm+ij, l)
324                 qs(ij,l)=q(ij,l)          qn(ip1jm+ij, l) = q(ip1jm+ij, l)
325                 qn(ij,l)=q(ij,l)        END DO
326  c              if (.not.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)  
327  c              if (.not.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l)      END DO
328              else    END IF
329                 qs(ij,l)=zqv(ij,l)    RETURN
330                 qn(ij,l)=zqv(ij-iip1,l)  END SUBROUTINE advnqy
331              endif  
332           enddo  SUBROUTINE advnqz(q, qh, qb)
333    
334           do ij=1,iip1    ! Auteurs:   Calcul des valeurs de q aux point v.
335              qs(ij,l)=q(ij,l)  
336              qn(ij,l)=q(ij,l)    ! --------------------------------------------------------------------
337              qs(ip1jm+ij,l)=q(ip1jm+ij,l)    USE dimens_m
338              qn(ip1jm+ij,l)=q(ip1jm+ij,l)    USE paramet_m
339           enddo    USE conf_gcm_m
340      IMPLICIT NONE
341        enddo  
342        endif  
343        RETURN  
344        END    ! Arguments:
345      ! ----------
346        SUBROUTINE advnqz(q,qh,qb)    REAL q(ip1jmp1, llm), qh(ip1jmp1, llm), qb(ip1jmp1, llm)
347  c  
348  c     Auteurs:   Calcul des valeurs de q aux point v.    ! Local
349  c    ! ---------
350  c   --------------------------------------------------------------------  
351        use dimens_m    INTEGER ij, l
352        use paramet_m  
353        use iniprint    REAL dzqw(ip1jmp1, llm+1), zqw(ip1jmp1, llm+1)
354        IMPLICIT NONE    REAL zqmax(ip1jmp1, llm), zqmin(ip1jmp1, llm)
355  c    LOGICAL extremum(ip1jmp1, llm)
356  c  
357  c    INTEGER mode
358  c   Arguments:    SAVE mode
359  c   ----------  
360        real q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)    DATA mode/1/
361  c  
362  c      Local    ! calcul des pentes en u:
363  c   ---------    ! -----------------------
364  c  
365        INTEGER ij,l    IF (mode==0) THEN
366  c      DO l = 1, llm
367        real dzqw(ip1jmp1,llm+1),zqw(ip1jmp1,llm+1)        DO ij = 1, ip1jmp1
368        real zqmax(ip1jmp1,llm),zqmin(ip1jmp1,llm)          qb(ij, l) = q(ij, l)
369        logical extremum(ip1jmp1,llm)          qh(ij, l) = q(ij, l)
370          END DO
371        integer mode      END DO
372        save mode    ELSE
373        DO l = 2, llm
374        data mode/1/        DO ij = 1, ip1jmp1
375            dzqw(ij, l) = q(ij, l-1) - q(ij, l)
376  c   calcul des pentes en u:          zqw(ij, l) = 0.5*(q(ij,l-1)+q(ij,l))
377  c   -----------------------        END DO
378        END DO
379        if (mode.eq.0) then      DO ij = 1, ip1jmp1
380           do l=1,llm        dzqw(ij, 1) = 0.
381              do ij=1,ip1jmp1        dzqw(ij, llm+1) = 0.
382                 qb(ij,l)=q(ij,l)      END DO
383                 qh(ij,l)=q(ij,l)      DO l = 2, llm
384              enddo        DO ij = 1, ip1jmp1
385           enddo          zqw(ij, l) = zqw(ij, l) + (dzqw(ij,l+1)-dzqw(ij,l-1))/12.
386        else        END DO
387        do l = 2, llm      END DO
388           do ij=1,ip1jmp1      DO l = 2, llm - 1
389              dzqw(ij,l)=q(ij,l-1)-q(ij,l)        DO ij = 1, ip1jmp1
390              zqw(ij,l)=0.5*(q(ij,l-1)+q(ij,l))          extremum(ij, l) = dzqw(ij, l)*dzqw(ij, l+1) <= 0.
391           enddo        END DO
392        enddo      END DO
393        do ij=1,ip1jmp1  
394           dzqw(ij,1)=0.      ! Pas de pentes en bas et en haut
395           dzqw(ij,llm+1)=0.      DO ij = 1, ip1jmp1
396        enddo        zqw(ij, 2) = q(ij, 1)
397        do l=2,llm        zqw(ij, llm) = q(ij, llm)
398           do ij=1,ip1jmp1        extremum(ij, 1) = .TRUE.
399              zqw(ij,l)=zqw(ij,l)+(dzqw(ij,l+1)-dzqw(ij,l-1))/12.        extremum(ij, llm) = .TRUE.
400           enddo      END DO
401        enddo  
402        do l=2,llm-1      ! calcul des valeurs max et min acceptees aux interfaces
403           do ij=1,ip1jmp1      DO l = 2, llm
404              extremum(ij,l)=dzqw(ij,l)*dzqw(ij,l+1).le.0.        DO ij = 1, ip1jmp1
405           enddo          zqmax(ij, l) = max(q(ij,l-1), q(ij,l))
406        enddo          zqmin(ij, l) = min(q(ij,l-1), q(ij,l))
407          END DO
408  c Pas de pentes en bas et en haut      END DO
409           do ij=1,ip1jmp1  
410              zqw(ij,2)=q(ij,1)      DO l = 2, llm
411              zqw(ij,llm)=q(ij,llm)        DO ij = 1, ip1jmp1
412              extremum(ij,1)=.true.          zqw(ij, l) = min(max(zqmin(ij,l),zqw(ij,l)), zqmax(ij,l))
413              extremum(ij,llm)=.true.        END DO
414           enddo      END DO
415    
416  c   calcul des valeurs max et min acceptees aux interfaces      DO l = 2, llm - 1
417        do l=2,llm        DO ij = 1, ip1jmp1
418           do ij=1,ip1jmp1          IF (extremum(ij,l)) THEN
419              zqmax(ij,l)=max(q(ij,l-1),q(ij,l))            qh(ij, l) = q(ij, l)
420              zqmin(ij,l)=min(q(ij,l-1),q(ij,l))            qb(ij, l) = q(ij, l)
421           enddo          ELSE
422        enddo            qh(ij, l) = zqw(ij, l+1)
423              qb(ij, l) = zqw(ij, l)
424        do l=2,llm          END IF
425           do ij=1,ip1jmp1        END DO
426              zqw(ij,l)=min(max(zqmin(ij,l),zqw(ij,l)),zqmax(ij,l))      END DO
427           enddo      ! do l=2,llm-1
428        enddo      ! do ij=1,ip1jmp1
429        ! if(extremum(ij,l)) then
430        do l=2,llm-1      ! if (.not.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
431           do ij=1,ip1jmp1      ! if (.not.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
432              if(extremum(ij,l)) then      ! endif
433                 qh(ij,l)=q(ij,l)      ! enddo
434                 qb(ij,l)=q(ij,l)      ! enddo
435              else  
436                 qh(ij,l)=zqw(ij,l+1)      DO ij = 1, ip1jmp1
437                 qb(ij,l)=zqw(ij,l)        qb(ij, 1) = q(ij, 1)
438              endif        qh(ij, 1) = q(ij, 1)
439           enddo        qb(ij, llm) = q(ij, llm)
440        enddo        qh(ij, llm) = q(ij, llm)
441  c     do l=2,llm-1      END DO
442  c        do ij=1,ip1jmp1  
443  c           if(extremum(ij,l)) then    END IF
444  c              if (.not.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)  
445  c              if (.not.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)    RETURN
446  c           endif  END SUBROUTINE advnqz
447  c        enddo  
448  c     enddo  SUBROUTINE advnx(q, qg, qd, masse, u_m, mode)
449    
450        do ij=1,ip1jmp1    ! Auteur : F. Hourdin
451           qb(ij,1)=q(ij,1)  
452           qh(ij,1)=q(ij,1)    ! ********************************************************************
453           qb(ij,llm)=q(ij,llm)    ! Shema  d'advection " pseudo amont " .
454           qh(ij,llm)=q(ij,llm)    ! ********************************************************************
455        enddo    ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
456    
457        endif  
458      ! --------------------------------------------------------------------
459        RETURN    USE dimens_m
460        END    USE paramet_m
461      USE comconst
462        SUBROUTINE advnx(q,qg,qd,masse,u_m,mode)    USE disvert_m
463  c    USE conf_gcm_m
464  c     Auteur : F. Hourdin    IMPLICIT NONE
465  c  
466  c    ********************************************************************  
467  c     Shema  d'advection " pseudo amont " .  
468  c    ********************************************************************    ! Arguments:
469  c     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....    ! ----------
470  c    INTEGER mode
471  c    REAL masse(ip1jmp1, llm)
472  c   --------------------------------------------------------------------    REAL u_m(ip1jmp1, llm)
473        use dimens_m    REAL q(ip1jmp1, llm), qd(ip1jmp1, llm), qg(ip1jmp1, llm)
474        use paramet_m  
475        use comconst    ! Local
476        use comvert    ! ---------
477        use logic  
478        use iniprint    INTEGER i, j, ij, l, indu(ip1jmp1), niju, iju, ijq
479        IMPLICIT NONE    INTEGER n0, nl(llm)
480  c  
481  c    REAL new_m, zu_m, zdq, zz
482  c    REAL zsigg(ip1jmp1, llm), zsigd(ip1jmp1, llm), zsig
483  c   Arguments:    REAL u_mq(ip1jmp1, llm)
484  c   ----------  
485        integer mode    REAL zm, zq, zsigm, zsigp, zqm, zqp, zu
486        real masse(ip1jmp1,llm)  
487        real u_m( ip1jmp1,llm )    LOGICAL ladvplus(ip1jmp1, llm)
488        real q(ip1jmp1,llm),qd(ip1jmp1,llm),qg(ip1jmp1,llm)  
489  c    REAL prec
490  c      Local    SAVE prec
491  c   ---------  
492  c    DATA prec/1.E-15/
493        INTEGER i,j,ij,l,indu(ip1jmp1),niju,iju,ijq  
494        integer n0,nl(llm)    DO l = 1, llm
495  c      DO ij = iip2, ip1jm
496        real new_m,zu_m,zdq,zz        zdq = qd(ij, l) - qg(ij, l)
497        real zsigg(ip1jmp1,llm),zsigd(ip1jmp1,llm),zsig        IF (abs(zdq)>prec) THEN
498        real u_mq(ip1jmp1,llm)          zsigd(ij, l) = (q(ij,l)-qg(ij,l))/zdq
499            zsigg(ij, l) = 1. - zsigd(ij, l)
500        real zm,zq,zsigm,zsigp,zqm,zqp,zu        ELSE
501            zsigd(ij, l) = 0.5
502        logical ladvplus(ip1jmp1,llm)          zsigg(ij, l) = 0.5
503            qd(ij, l) = q(ij, l)
504        real prec          qg(ij, l) = q(ij, l)
505        save prec        END IF
506        END DO
507        data prec/1.e-15/    END DO
508    
509        do l=1,llm    ! calcul de la pente maximum dans la maille en valeur absolue
510              do ij=iip2,ip1jm  
511                 zdq=qd(ij,l)-qg(ij,l)    DO l = 1, llm
512  c              if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) then      DO ij = iip2, ip1jm - 1
513  c                 print*,'probleme au point ij=',ij,'  l=',l        IF (u_m(ij,l)>=0.) THEN
514  c                 print*,qd(ij,l),q(ij,l),qg(ij,l)          zsigp = zsigd(ij, l)
515  c                 qd(ij,l)=q(ij,l)          zsigm = zsigg(ij, l)
516  c                 qg(ij,l)=q(ij,l)          zqp = qd(ij, l)
517  c              endif          zqm = qg(ij, l)
518                 if(abs(zdq).gt.prec) then          zm = masse(ij, l)
519                    zsigd(ij,l)=(q(ij,l)-qg(ij,l))/zdq          zq = q(ij, l)
520                    zsigg(ij,l)=1.-zsigd(ij,l)        ELSE
521  c                 if(.not.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .and.          zsigm = zsigd(ij+1, l)
522  c    s               zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) then          zsigp = zsigg(ij+1, l)
523  c                    print*,'probleme au point ij=',ij,'  l=',l          zqm = qd(ij+1, l)
524  c                    print*,'sigg=',zsigg(ij,l),'  sigd=',zsigd(ij,l)          zqp = qg(ij+1, l)
525  c                    print*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq          zm = masse(ij+1, l)
526  c                    stop          zq = q(ij+1, l)
527  c                 endif        END IF
528                 else        zu = abs(u_m(ij,l))
529                    zsigd(ij,l)=0.5        ladvplus(ij, l) = zu > zm
530                    zsigg(ij,l)=0.5        zsig = zu/zm
531                    qd(ij,l)=q(ij,l)        IF (zsig==0.) zsigp = 0.1
532                    qg(ij,l)=q(ij,l)        IF (mode==1) THEN
533                 endif          IF (zsig<=zsigp) THEN
534              enddo            u_mq(ij, l) = u_m(ij, l)*zqp
535         enddo          ELSE IF (mode==1) THEN
536              u_mq(ij, l) = sign(zm, u_m(ij,l))*(zsigp*zqp+(zsig-zsigp)*zqm)
537  c   calcul de la pente maximum dans la maille en valeur absolue          END IF
538          ELSE
539         do l=1,llm          IF (zsig<=zsigp) THEN
540         do ij=iip2,ip1jm-1            u_mq(ij, l) = u_m(ij, l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
541            if (u_m(ij,l).ge.0.) then          ELSE
542               zsigp=zsigd(ij,l)            zz = 0.5*(zsig-zsigp)/zsigm
543               zsigm=zsigg(ij,l)            u_mq(ij, l) = sign(zm, u_m(ij,l))*(0.5*(zq+zqp)*zsigp+(zsig-zsigp)* &
544               zqp=qd(ij,l)              (zq+zz*(zqm-zq)))
545               zqm=qg(ij,l)          END IF
546               zm=masse(ij,l)        END IF
547               zq=q(ij,l)      END DO
548            else    END DO
549               zsigm=zsigd(ij+1,l)  
550               zsigp=zsigg(ij+1,l)    DO l = 1, llm
551               zqm=qd(ij+1,l)      DO ij = iip1 + iip1, ip1jm, iip1
552               zqp=qg(ij+1,l)        u_mq(ij, l) = u_mq(ij-iim, l)
553               zm=masse(ij+1,l)        ladvplus(ij, l) = ladvplus(ij-iim, l)
554               zq=q(ij+1,l)      END DO
555            endif    END DO
556            zu=abs(u_m(ij,l))  
557            ladvplus(ij,l)=zu.gt.zm    ! =================================================================
558            zsig=zu/zm    ! SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES
559            if(zsig.eq.0.) zsigp=0.1    ! =================================================================
560            if (mode.eq.1) then    ! tris des regions a traiter
561               if (zsig.le.zsigp) then    n0 = 0
562                   u_mq(ij,l)=u_m(ij,l)*zqp    DO l = 1, llm
563               else if (mode.eq.1) then      nl(l) = 0
564                   u_mq(ij,l)=      DO ij = iip2, ip1jm
565       s           sign(zm,u_m(ij,l))*(zsigp*zqp+(zsig-zsigp)*zqm)        IF (ladvplus(ij,l)) THEN
566               endif          nl(l) = nl(l) + 1
567            else          u_mq(ij, l) = 0.
568               if (zsig.le.zsigp) then        END IF
569                   u_mq(ij,l)=u_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))      END DO
570               else      n0 = n0 + nl(l)
571                  zz=0.5*(zsig-zsigp)/zsigm    END DO
572                  u_mq(ij,l)=sign(zm,u_m(ij,l))*( 0.5*(zq+zqp)*zsigp  
573       s          +(zsig-zsigp)*(zq+zz*(zqm-zq)) )    IF (n0>1) THEN
574               endif      IF (prt_level>9) PRINT *, &
575            endif        'Nombre de points pour lesquels on advect plus que le', &
576  c         if(zsig.lt.0.) then        'contenu de la maille : ', n0
577  c            print*,'au point ij=',ij,'  l=',l,'  sig=',zsig  
578  c            stop      DO l = 1, llm
579  c         endif        IF (nl(l)>0) THEN
580        enddo          iju = 0
581        enddo          ! indicage des mailles concernees par le traitement special
582            DO ij = iip2, ip1jm
583        do l=1,llm            IF (ladvplus(ij,l) .AND. mod(ij,iip1)/=0) THEN
584         do ij=iip1+iip1,ip1jm,iip1              iju = iju + 1
585            u_mq(ij,l)=u_mq(ij-iim,l)              indu(iju) = ij
586            ladvplus(ij,l)=ladvplus(ij-iim,l)            END IF
587         enddo          END DO
588        enddo          niju = iju
589    
590  c=================================================================          ! traitement des mailles
591  C   SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES          DO iju = 1, niju
592  c=================================================================            ij = indu(iju)
593  c   tris des regions a traiter            j = (ij-1)/iip1 + 1
594        n0=0            zu_m = u_m(ij, l)
595        do l=1,llm            u_mq(ij, l) = 0.
596           nl(l)=0            IF (zu_m>0.) THEN
597           do ij=iip2,ip1jm              ijq = ij
598              if(ladvplus(ij,l)) then              i = ijq - (j-1)*iip1
599                 nl(l)=nl(l)+1              ! accumulation pour les mailles completements advectees
600                 u_mq(ij,l)=0.              DO WHILE (zu_m>masse(ijq,l))
601              endif                u_mq(ij, l) = u_mq(ij, l) + q(ijq, l)*masse(ijq, l)
602           enddo                zu_m = zu_m - masse(ijq, l)
603           n0=n0+nl(l)                i = mod(i-2+iim, iim) + 1
604        enddo                ijq = (j-1)*iip1 + i
605                END DO
606        if(n0.gt.1) then              ! MODIFS SPECIFIQUES DU SCHEMA
607        IF (prt_level > 9) print *,              ! ajout de la maille non completement advectee
608       & 'Nombre de points pour lesquels on advect plus que le'              zsig = zu_m/masse(ijq, l)
609       &       ,'contenu de la maille : ',n0              IF (zsig<=zsigd(ijq,l)) THEN
610                  u_mq(ij, l) = u_mq(ij, l) + zu_m*(qd(ijq,l)-0.5*zsig/zsigd(ijq, &
611           do l=1,llm                  l)*(qd(ijq,l)-q(ijq,l)))
612              if(nl(l).gt.0) then              ELSE
613                 iju=0                ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
614  c   indicage des mailles concernees par le traitement special                ! goto 8888
615                 do ij=iip2,ip1jm                zz = 0.5*(zsig-zsigd(ijq,l))/zsigg(ijq, l)
616                    if(ladvplus(ij,l).and.mod(ij,iip1).ne.0) then                IF (.NOT. (zz>0. .AND. zz<=0.5)) THEN
617                       iju=iju+1                  PRINT *, 'probleme2 au point ij=', ij, '  l=', l
618                       indu(iju)=ij                  PRINT *, 'zz=', zz
619                    endif                  STOP
620                 enddo                END IF
621                 niju=iju                u_mq(ij, l) = u_mq(ij, l) + masse(ijq, l)*(0.5*(q(ijq, &
622  c              print*,'niju,nl',niju,nl(l)                  l)+qd(ijq,l))*zsigd(ijq,l)+(zsig-zsigd(ijq,l))*(q(ijq, &
623                    l)+zz*(qg(ijq,l)-q(ijq,l))))
624  c  traitement des mailles              END IF
625                 do iju=1,niju            ELSE
626                    ij=indu(iju)              ijq = ij + 1
627                    j=(ij-1)/iip1+1              i = ijq - (j-1)*iip1
628                    zu_m=u_m(ij,l)              ! accumulation pour les mailles completements advectees
629                    u_mq(ij,l)=0.              DO WHILE (-zu_m>masse(ijq,l))
630                    if(zu_m.gt.0.) then                u_mq(ij, l) = u_mq(ij, l) - q(ijq, l)*masse(ijq, l)
631                       ijq=ij                zu_m = zu_m + masse(ijq, l)
632                       i=ijq-(j-1)*iip1                i = mod(i, iim) + 1
633  c   accumulation pour les mailles completements advectees                ijq = (j-1)*iip1 + i
634                       do while(zu_m.gt.masse(ijq,l))              END DO
635                          u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)              ! ajout de la maille non completement advectee
636                          zu_m=zu_m-masse(ijq,l)              ! 2eme MODIF SPECIFIQUE
637                          i=mod(i-2+iim,iim)+1              zsig = -zu_m/masse(ij+1, l)
638                          ijq=(j-1)*iip1+i              IF (zsig<=zsigg(ijq,l)) THEN
639                       enddo                u_mq(ij, l) = u_mq(ij, l) + zu_m*(qg(ijq,l)-0.5*zsig/zsigg(ijq, &
640  c   MODIFS SPECIFIQUES DU SCHEMA                  l)*(qg(ijq,l)-q(ijq,l)))
641  c   ajout de la maille non completement advectee              ELSE
642               zsig=zu_m/masse(ijq,l)                ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
643               if(zsig.le.zsigd(ijq,l)) then                ! goto 9999
644                  u_mq(ij,l)=u_mq(ij,l)+zu_m*(qd(ijq,l)                zz = 0.5*(zsig-zsigg(ijq,l))/zsigd(ijq, l)
645       s          -0.5*zsig/zsigd(ijq,l)*(qd(ijq,l)-q(ijq,l)))                IF (.NOT. (zz>0. .AND. zz<=0.5)) THEN
646               else                  PRINT *, 'probleme22 au point ij=', ij, '  l=', l
647  c               u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)                  PRINT *, 'zz=', zz
648  c         goto 8888                  STOP
649                  zz=0.5*(zsig-zsigd(ijq,l))/zsigg(ijq,l)                END IF
650                  if(.not.(zz.gt.0..and.zz.le.0.5)) then                u_mq(ij, l) = u_mq(ij, l) - masse(ijq, l)*(0.5*(q(ijq, &
651                       print *,'probleme2 au point ij=',ij,                  l)+qg(ijq,l))*zsigg(ijq,l)+(zsig-zsigg(ijq,l))*(q(ijq, &
652       s               '  l=',l                  l)+zz*(qd(ijq,l)-q(ijq,l))))
653                       print *,'zz=',zz              END IF
654                       stop              ! fin de la modif
655                  endif            END IF
656                  u_mq(ij,l)=u_mq(ij,l)+masse(ijq,l)*(          END DO
657       s          0.5*(q(ijq,l)+qd(ijq,l))*zsigd(ijq,l)        END IF
658       s        +(zsig-zsigd(ijq,l))*(q(ijq,l)+zz*(qg(ijq,l)-q(ijq,l))) )      END DO
659               endif    END IF ! n0.gt.0
660                    else  
661                       ijq=ij+1    ! bouclage en latitude
662                       i=ijq-(j-1)*iip1    DO l = 1, llm
663  c   accumulation pour les mailles completements advectees      DO ij = iip1 + iip1, ip1jm, iip1
664                       do while(-zu_m.gt.masse(ijq,l))        u_mq(ij, l) = u_mq(ij-iim, l)
665                          u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)      END DO
666                          zu_m=zu_m+masse(ijq,l)    END DO
667                          i=mod(i,iim)+1  
668                          ijq=(j-1)*iip1+i    ! =================================================================
669                       enddo    ! CALCUL DE LA CONVERGENCE DES FLUX
670  c   ajout de la maille non completement advectee    ! =================================================================
671  c 2eme MODIF SPECIFIQUE  
672               zsig=-zu_m/masse(ij+1,l)    DO l = 1, llm
673               if(zsig.le.zsigg(ijq,l)) then      DO ij = iip2 + 1, ip1jm
674                  u_mq(ij,l)=u_mq(ij,l)+zu_m*(qg(ijq,l)        new_m = masse(ij, l) + u_m(ij-1, l) - u_m(ij, l)
675       s          -0.5*zsig/zsigg(ijq,l)*(qg(ijq,l)-q(ijq,l)))        q(ij, l) = (q(ij,l)*masse(ij,l)+u_mq(ij-1,l)-u_mq(ij,l))/new_m
676               else        masse(ij, l) = new_m
677  c               u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)      END DO
678  c           goto 9999      ! Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
679                  zz=0.5*(zsig-zsigg(ijq,l))/zsigd(ijq,l)      DO ij = iip1 + iip1, ip1jm, iip1
680                  if(.not.(zz.gt.0..and.zz.le.0.5)) then        q(ij-iim, l) = q(ij, l)
681                       print *,'probleme22 au point ij=',ij        masse(ij-iim, l) = masse(ij, l)
682       s               ,'  l=',l      END DO
683                       print *,'zz=',zz    END DO
684                       stop  
685                  endif    RETURN
686                  u_mq(ij,l)=u_mq(ij,l)-masse(ijq,l)*(  END SUBROUTINE advnx
687       s          0.5*(q(ijq,l)+qg(ijq,l))*zsigg(ijq,l)  SUBROUTINE advny(q, qs, qn, masse, v_m)
688       s          +(zsig-zsigg(ijq,l))*  
689       s           (q(ijq,l)+zz*(qd(ijq,l)-q(ijq,l))) )    ! Auteur : F. Hourdin
690               endif  
691  c   fin de la modif    ! ********************************************************************
692                    endif    ! Shema  d'advection " pseudo amont " .
693                 enddo    ! ********************************************************************
694              endif    ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
695           enddo  
696        endif  ! n0.gt.0  
697      ! --------------------------------------------------------------------
698  c   bouclage en latitude    USE dimens_m
699        do l=1,llm    USE paramet_m
700          do ij=iip1+iip1,ip1jm,iip1    USE comgeom
701             u_mq(ij,l)=u_mq(ij-iim,l)    USE conf_gcm_m
702          enddo    IMPLICIT NONE
703        enddo  
704    
705  c=================================================================  
706  c   CALCUL DE LA CONVERGENCE DES FLUX    ! Arguments:
707  c=================================================================    ! ----------
708      REAL masse(ip1jmp1, llm)
709        do l=1,llm    REAL v_m(ip1jm, llm)
710           do ij=iip2+1,ip1jm    REAL q(ip1jmp1, llm), qn(ip1jmp1, llm), qs(ip1jmp1, llm)
711              new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)  
712              q(ij,l)=(q(ij,l)*masse(ij,l)+    ! Local
713       &      u_mq(ij-1,l)-u_mq(ij,l))    ! ---------
714       &      /new_m  
715              masse(ij,l)=new_m    INTEGER ij, l
716           enddo  
717  c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)    REAL new_m, zdq, zz
718           do ij=iip1+iip1,ip1jm,iip1    REAL zsigs(ip1jmp1), zsign(ip1jmp1), zsig
719              q(ij-iim,l)=q(ij,l)    REAL v_mq(ip1jm, llm)
720              masse(ij-iim,l)=masse(ij,l)    REAL convpn, convps, convmpn, convmps, massen, masses
721           enddo    REAL zm, zq, zsigm, zsigp, zqm, zqp
722        enddo    REAL ssum
723      REAL prec
724        RETURN    SAVE prec
725        END  
726        SUBROUTINE advny(q,qs,qn,masse,v_m)    DATA prec/1.E-15/
727  c  
728  c     Auteur : F. Hourdin    DO l = 1, llm
729  c      DO ij = 1, ip1jmp1
730  c    ********************************************************************        zdq = qn(ij, l) - qs(ij, l)
731  c     Shema  d'advection " pseudo amont " .        IF (abs(zdq)>prec) THEN
732  c    ********************************************************************          zsign(ij) = (q(ij,l)-qs(ij,l))/zdq
733  c     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....          zsigs(ij) = 1. - zsign(ij)
734  c        ELSE
735  c          zsign(ij) = 0.5
736  c   --------------------------------------------------------------------          zsigs(ij) = 0.5
737        use dimens_m        END IF
738        use paramet_m      END DO
739        use comgeom  
740        use iniprint      ! calcul de la pente maximum dans la maille en valeur absolue
741        IMPLICIT NONE  
742  c      DO ij = 1, ip1jm
743  c        IF (v_m(ij,l)>=0.) THEN
744  c          zsigp = zsign(ij+iip1)
745  c   Arguments:          zsigm = zsigs(ij+iip1)
746  c   ----------          zqp = qn(ij+iip1, l)
747        real masse(ip1jmp1,llm)          zqm = qs(ij+iip1, l)
748        real v_m( ip1jm,llm )          zm = masse(ij+iip1, l)
749        real q(ip1jmp1,llm),qn(ip1jmp1,llm),qs(ip1jmp1,llm)          zq = q(ij+iip1, l)
750  c        ELSE
751  c      Local          zsigm = zsign(ij)
752  c   ---------          zsigp = zsigs(ij)
753  c          zqm = qn(ij, l)
754        INTEGER ij,l          zqp = qs(ij, l)
755  c          zm = masse(ij, l)
756        real new_m,zdq,zz          zq = q(ij, l)
757        real zsigs(ip1jmp1),zsign(ip1jmp1),zsig        END IF
758        real v_mq(ip1jm,llm)        zsig = abs(v_m(ij,l))/zm
759        real convpn,convps,convmpn,convmps,massen,masses        IF (zsig==0.) zsigp = 0.1
760        real zm,zq,zsigm,zsigp,zqm,zqp        IF (zsig<=zsigp) THEN
761        real ssum          v_mq(ij, l) = v_m(ij, l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
762        real prec        ELSE
763        save prec          zz = 0.5*(zsig-zsigp)/zsigm
764            v_mq(ij, l) = sign(zm, v_m(ij,l))*(0.5*(zq+zqp)*zsigp+(zsig-zsigp)*( &
765        data prec/1.e-15/            zq+zz*(zqm-zq)))
766        do l=1,llm        END IF
767              do ij=1,ip1jmp1      END DO
768                 zdq=qn(ij,l)-qs(ij,l)    END DO
769  c              if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) then  
770  c                 print*,'probleme au point ij=',ij,'  l=',l,'  advnqx'    DO l = 1, llm
771  c                 print*,qn(ij,l),q(ij,l),qs(ij,l)      DO ij = iip2, ip1jm
772  c                 qn(ij,l)=q(ij,l)        new_m = masse(ij, l) + v_m(ij, l) - v_m(ij-iip1, l)
773  c                 qs(ij,l)=q(ij,l)        q(ij, l) = (q(ij,l)*masse(ij,l)+v_mq(ij,l)-v_mq(ij-iip1,l))/new_m
774  c              endif        masse(ij, l) = new_m
775                 if(abs(zdq).gt.prec) then      END DO
776                    zsign(ij)=(q(ij,l)-qs(ij,l))/zdq      ! .-. ancienne version
777                    zsigs(ij)=1.-zsign(ij)      convpn = ssum(iim, v_mq(1,l), 1)
778  c                 if(.not.(zsign(ij).ge.0..and.zsign(ij).le.1. .and.      convmpn = ssum(iim, v_m(1,l), 1)
779  c    s               zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) then      massen = ssum(iim, masse(1,l), 1)
780  c                    print*,'probleme au point ij=',ij,'  l=',l      new_m = massen + convmpn
781  c                    print*,'sigs=',zsigs(ij),'  sign=',zsign(ij)      q(1, l) = (q(1,l)*massen+convpn)/new_m
782  c                    stop      DO ij = 1, iip1
783  c                 endif        q(ij, l) = q(1, l)
784                 else        masse(ij, l) = new_m*aire(ij)/apoln
785                    zsign(ij)=0.5      END DO
786                    zsigs(ij)=0.5  
787                 endif      convps = -ssum(iim, v_mq(ip1jm-iim,l), 1)
788              enddo      convmps = -ssum(iim, v_m(ip1jm-iim,l), 1)
789        masses = ssum(iim, masse(ip1jm+1,l), 1)
790  c   calcul de la pente maximum dans la maille en valeur absolue      new_m = masses + convmps
791        q(ip1jm+1, l) = (q(ip1jm+1,l)*masses+convps)/new_m
792         do ij=1,ip1jm      DO ij = ip1jm + 1, ip1jmp1
793            if (v_m(ij,l).ge.0.) then        q(ij, l) = q(ip1jm+1, l)
794               zsigp=zsign(ij+iip1)        masse(ij, l) = new_m*aire(ij)/apols
795               zsigm=zsigs(ij+iip1)      END DO
796               zqp=qn(ij+iip1,l)    END DO
797               zqm=qs(ij+iip1,l)  
798               zm=masse(ij+iip1,l)    RETURN
799               zq=q(ij+iip1,l)  END SUBROUTINE advny
800            else  SUBROUTINE advnz(q, qh, qb, masse, w_m)
801               zsigm=zsign(ij)  
802               zsigp=zsigs(ij)    ! Auteurs:   F.Hourdin
803               zqm=qn(ij,l)  
804               zqp=qs(ij,l)    ! ********************************************************************
805               zm=masse(ij,l)    ! Shema  d'advection " pseudo amont " .
806               zq=q(ij,l)    ! b designe le bas et h le haut
807            endif    ! il y a une correspondance entre le b en z et le d en x
808            zsig=abs(v_m(ij,l))/zm    ! ********************************************************************
809            if(zsig.eq.0.) zsigp=0.1  
810            if (zsig.le.zsigp) then  
811                v_mq(ij,l)=v_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))    ! --------------------------------------------------------------------
812            else    USE dimens_m
813                zz=0.5*(zsig-zsigp)/zsigm    USE paramet_m
814                v_mq(ij,l)=sign(zm,v_m(ij,l))*( 0.5*(zq+zqp)*zsigp    USE comgeom
815       s        +(zsig-zsigp)*(zq+zz*(zqm-zq)) )    USE conf_gcm_m
816            endif    IMPLICIT NONE
817         enddo  
818        enddo  
819    
820        do l=1,llm    ! Arguments:
821           do ij=iip2,ip1jm    ! ----------
822              new_m=masse(ij,l)    REAL masse(ip1jmp1, llm)
823       &      +v_m(ij,l)-v_m(ij-iip1,l)    REAL w_m(ip1jmp1, llm+1)
824              q(ij,l)=(q(ij,l)*masse(ij,l)+v_mq(ij,l)-v_mq(ij-iip1,l))    REAL q(ip1jmp1, llm), qb(ip1jmp1, llm), qh(ip1jmp1, llm)
825       &         /new_m  
826              masse(ij,l)=new_m  
827           enddo    ! Local
828  c.-. ancienne version    ! ---------
829           convpn=SSUM(iim,v_mq(1,l),1)  
830           convmpn=ssum(iim,v_m(1,l),1)    INTEGER ij, l
831           massen=ssum(iim,masse(1,l),1)  
832           new_m=massen+convmpn    REAL new_m, zdq, zz
833           q(1,l)=(q(1,l)*massen+convpn)/new_m    REAL zsigh(ip1jmp1, llm), zsigb(ip1jmp1, llm), zsig
834           do ij = 1,iip1    REAL w_mq(ip1jmp1, llm+1)
835              q(ij,l)=q(1,l)    REAL zm, zq, zsigm, zsigp, zqm, zqp
836              masse(ij,l)=new_m*aire(ij)/apoln    REAL prec
837           enddo    SAVE prec
838    
839           convps=-SSUM(iim,v_mq(ip1jm-iim,l),1)    DATA prec/1.E-13/
840           convmps=-ssum(iim,v_m(ip1jm-iim,l),1)  
841           masses=ssum(iim,masse(ip1jm+1,l),1)    DO l = 1, llm
842           new_m=masses+convmps      DO ij = 1, ip1jmp1
843           q(ip1jm+1,l)=(q(ip1jm+1,l)*masses+convps)/new_m        zdq = qb(ij, l) - qh(ij, l)
844           do ij = ip1jm+1,ip1jmp1        IF (abs(zdq)>prec) THEN
845              q(ij,l)=q(ip1jm+1,l)          zsigb(ij, l) = (q(ij,l)-qh(ij,l))/zdq
846              masse(ij,l)=new_m*aire(ij)/apols          zsigh(ij, l) = 1. - zsigb(ij, l)
847           enddo          zsigb(ij, l) = min(max(zsigb(ij,l),0.), 1.)
848        enddo        ELSE
849            zsigb(ij, l) = 0.5
850        RETURN          zsigh(ij, l) = 0.5
851        END        END IF
852        SUBROUTINE advnz(q,qh,qb,masse,w_m)      END DO
853  c    END DO
854  c     Auteurs:   F.Hourdin  
855  c    ! calcul de la pente maximum dans la maille en valeur absolue
856  c    ********************************************************************    DO l = 2, llm
857  c     Shema  d'advection " pseudo amont " .      DO ij = 1, ip1jmp1
858  c     b designe le bas et h le haut        IF (w_m(ij,l)>=0.) THEN
859  c     il y a une correspondance entre le b en z et le d en x          zsigp = zsigb(ij, l)
860  c    ********************************************************************          zsigm = zsigh(ij, l)
861  c          zqp = qb(ij, l)
862  c          zqm = qh(ij, l)
863  c   --------------------------------------------------------------------          zm = masse(ij, l)
864        use dimens_m          zq = q(ij, l)
865        use paramet_m        ELSE
866        use comgeom          zsigm = zsigb(ij, l-1)
867        use iniprint          zsigp = zsigh(ij, l-1)
868        IMPLICIT NONE          zqm = qb(ij, l-1)
869  c          zqp = qh(ij, l-1)
870  c          zm = masse(ij, l-1)
871  c          zq = q(ij, l-1)
872  c   Arguments:        END IF
873  c   ----------        zsig = abs(w_m(ij,l))/zm
874        real masse(ip1jmp1,llm)        IF (zsig==0.) zsigp = 0.1
875        real w_m( ip1jmp1,llm+1)        IF (zsig<=zsigp) THEN
876        real q(ip1jmp1,llm),qb(ip1jmp1,llm),qh(ip1jmp1,llm)          w_mq(ij, l) = w_m(ij, l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
877          ELSE
878  c          zz = 0.5*(zsig-zsigp)/zsigm
879  c      Local          w_mq(ij, l) = sign(zm, w_m(ij,l))*(0.5*(zq+zqp)*zsigp+(zsig-zsigp)*( &
880  c   ---------            zq+zz*(zqm-zq)))
881  c        END IF
882        INTEGER ij,l      END DO
883  c    END DO
884        real new_m,zdq,zz  
885        real zsigh(ip1jmp1,llm),zsigb(ip1jmp1,llm),zsig    DO ij = 1, ip1jmp1
886        real w_mq(ip1jmp1,llm+1)      w_mq(ij, llm+1) = 0.
887        real zm,zq,zsigm,zsigp,zqm,zqp      w_mq(ij, 1) = 0.
888        real prec    END DO
889        save prec  
890      DO l = 1, llm
891        data prec/1.e-13/      DO ij = 1, ip1jmp1
892          new_m = masse(ij, l) + w_m(ij, l+1) - w_m(ij, l)
893        do l=1,llm        q(ij, l) = (q(ij,l)*masse(ij,l)+w_mq(ij,l+1)-w_mq(ij,l))/new_m
894              do ij=1,ip1jmp1        masse(ij, l) = new_m
895                 zdq=qb(ij,l)-qh(ij,l)      END DO
896  c              if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) then    END DO
897  c                 print*,'probleme au point ij=',ij,'  l=',l  
898  c                 print*,qh(ij,l),q(ij,l),qb(ij,l)  END SUBROUTINE advnz
 c                 qh(ij,l)=q(ij,l)  
 c                 qb(ij,l)=q(ij,l)  
 c              endif  
   
                if(abs(zdq).gt.prec) then  
                   zsigb(ij,l)=(q(ij,l)-qh(ij,l))/zdq  
                   zsigh(ij,l)=1.-zsigb(ij,l)  
                   zsigb(ij,l)=min(max(zsigb(ij,l),0.),1.)  
                else  
                   zsigb(ij,l)=0.5  
                   zsigh(ij,l)=0.5  
                endif  
             enddo  
        enddo  
   
 c      print*,'ok1'  
 c   calcul de la pente maximum dans la maille en valeur absolue  
        do l=2,llm  
        do ij=1,ip1jmp1  
           if (w_m(ij,l).ge.0.) then  
              zsigp=zsigb(ij,l)  
              zsigm=zsigh(ij,l)  
              zqp=qb(ij,l)  
              zqm=qh(ij,l)  
              zm=masse(ij,l)  
              zq=q(ij,l)  
           else  
              zsigm=zsigb(ij,l-1)  
              zsigp=zsigh(ij,l-1)  
              zqm=qb(ij,l-1)  
              zqp=qh(ij,l-1)  
              zm=masse(ij,l-1)  
              zq=q(ij,l-1)  
           endif  
           zsig=abs(w_m(ij,l))/zm  
           if(zsig.eq.0.) zsigp=0.1  
           if (zsig.le.zsigp) then  
               w_mq(ij,l)=w_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))  
           else  
               zz=0.5*(zsig-zsigp)/zsigm  
               w_mq(ij,l)=sign(zm,w_m(ij,l))*( 0.5*(zq+zqp)*zsigp  
      s        +(zsig-zsigp)*(zq+zz*(zqm-zq)) )  
           endif  
       enddo  
       enddo  
   
        do ij=1,ip1jmp1  
           w_mq(ij,llm+1)=0.  
           w_mq(ij,1)=0.  
        enddo  
   
       do l=1,llm  
          do ij=1,ip1jmp1  
             new_m=masse(ij,l)+w_m(ij,l+1)-w_m(ij,l)  
             q(ij,l)=(q(ij,l)*masse(ij,l)+w_mq(ij,l+1)-w_mq(ij,l))  
      &         /new_m  
             masse(ij,l)=new_m  
          enddo  
       enddo  
 c     print*,'ok3'  
       RETURN  
       END  
Legend:
Removed from v.31  
changed lines
  Added in v.157
  ViewVC Help
Powered by ViewVC 1.1.21