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