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

Diff of /trunk/dyn3d/advn.f

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

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