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

Diff of /trunk/dyn3d/pentes_ini.f

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

trunk/dyn3d/pentes_ini.f revision 76 by guez, Fri Nov 15 18:45:49 2013 UTC trunk/dyn3d/pentes_ini.f90 revision 81 by guez, Wed Mar 5 14:38:41 2014 UTC
# Line 1  Line 1 
1  !  
2  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/pentes_ini.F,v 1.1.1.1 2004/05/19 12:53:07 lmdzadmin Exp $  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/pentes_ini.F,v 1.1.1.1 2004/05/19
3  !  ! 12:53:07 lmdzadmin Exp $
4        SUBROUTINE pentes_ini (q,w,masse,pbaru,pbarv,mode)  
5        use dimens_m  SUBROUTINE pentes_ini(q, w, masse, pbaru, pbarv, mode)
6        use paramet_m    USE dimens_m
7        use comconst    USE paramet_m
8        use disvert_m    USE comconst
9        use comgeom    USE disvert_m
10        USE nr_util, ONLY : pi    USE comgeom
11        IMPLICIT NONE    USE nr_util, ONLY: pi
12      IMPLICIT NONE
13  c=======================================================================  
14  c   Adaptation LMDZ:  A.Armengaud (LGGE)    ! =======================================================================
15  c   ----------------    ! Adaptation LMDZ:  A.Armengaud (LGGE)
16  c    ! ----------------
17  c   ********************************************************************  
18  c   Transport des traceurs par la methode des pentes    ! ********************************************************************
19  c   ********************************************************************    ! Transport des traceurs par la methode des pentes
20  c   Reference possible : Russel. G.L., Lerner J.A.:    ! ********************************************************************
21  c         A new Finite-Differencing Scheme for Traceur Transport    ! Reference possible : Russel. G.L., Lerner J.A.:
22  c         Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81    ! A new Finite-Differencing Scheme for Traceur Transport
23  c   ********************************************************************    ! Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81
24  c   q,w,masse,pbaru et pbarv    ! ********************************************************************
25  c                      sont des arguments d'entree  pour le s-pg ....    ! q,w,masse,pbaru et pbarv
26  c    ! sont des arguments d'entree  pour le s-pg ....
27  c=======================================================================  
28      ! =======================================================================
29    
30    
31  c   Arguments:  
32  c   ----------    ! Arguments:
33        integer mode    ! ----------
34        REAL, intent(in):: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )    INTEGER mode
35        REAL q( iip1,jjp1,llm,0:3)    REAL, INTENT (IN) :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
36        REAL w( ip1jmp1,llm )    REAL q(iip1, jjp1, llm, 0:3)
37        REAL masse( iip1,jjp1,llm)    REAL w(ip1jmp1, llm)
38  c   Local:    REAL masse(iip1, jjp1, llm)
39  c   ------    ! Local:
40        LOGICAL limit    ! ------
41        REAL sm ( iip1,jjp1, llm )    LOGICAL limit
42        REAL s0( iip1,jjp1,llm ),  sx( iip1,jjp1,llm )    REAL sm(iip1, jjp1, llm)
43        REAL sy( iip1,jjp1,llm ),  sz( iip1,jjp1,llm )    REAL s0(iip1, jjp1, llm), sx(iip1, jjp1, llm)
44        real masn,mass,zz    REAL sy(iip1, jjp1, llm), sz(iip1, jjp1, llm)
45        INTEGER i,j,l,iq    REAL masn, mass, zz
46      INTEGER i, j, l, iq
47  c  modif Fred 24 03 96  
48      ! modif Fred 24 03 96
49        real sinlon(iip1),sinlondlon(iip1)  
50        real coslon(iip1),coslondlon(iip1)    REAL sinlon(iip1), sinlondlon(iip1)
51        save sinlon,coslon,sinlondlon,coslondlon    REAL coslon(iip1), coslondlon(iip1)
52        real dyn1,dyn2,dys1,dys2    SAVE sinlon, coslon, sinlondlon, coslondlon
53        real qpn,qps,dqzpn,dqzps    REAL dyn1, dyn2, dys1, dys2
54        real smn,sms,s0n,s0s,sxn(iip1),sxs(iip1)    REAL qpn, qps, dqzpn, dqzps
55        real qmin,zq,pente_max    REAL smn, sms, s0n, s0s, sxn(iip1), sxs(iip1)
56  c    REAL qmin, zq, pente_max
57        REAL      SSUM  
58        integer ismax,ismin,lati,latf    REAL ssum
59        EXTERNAL  SSUM, convflu,ismin,ismax    INTEGER ismax, ismin, lati, latf
60        logical first    EXTERNAL ssum, convflu, ismin, ismax
61        save first    LOGICAL first
62  c   fin modif    SAVE first
63      ! fin modif
64  c      EXTERNAL masskg  
65        EXTERNAL advx    ! EXTERNAL masskg
66        EXTERNAL advy    EXTERNAL advx
67        EXTERNAL advz    EXTERNAL advy
68      EXTERNAL advz
69  c  modif Fred 24 03 96  
70        data first/.true./    ! modif Fred 24 03 96
71      DATA first/.TRUE./
72        limit = .TRUE.  
73        pente_max=2    limit = .TRUE.
74  c     if (mode.eq.1.or.mode.eq.3) then    pente_max = 2
75  c     if (mode.eq.1) then    ! if (mode.eq.1.or.mode.eq.3) then
76        if (mode.ge.1) then    ! if (mode.eq.1) then
77          lati=2    IF (mode>=1) THEN
78          latf=jjm      lati = 2
79        else      latf = jjm
80          lati=1    ELSE
81          latf=jjp1      lati = 1
82        endif      latf = jjp1
83      END IF
84        qmin=0.4995  
85        qmin=0.    qmin = 0.4995
86        if(first) then    qmin = 0.
87           print*,'SCHEMA AMONT NOUVEAU'    IF (first) THEN
88           first=.false.      PRINT *, 'SCHEMA AMONT NOUVEAU'
89           do i=2,iip1      first = .FALSE.
90              coslon(i)=cos(rlonv(i))      DO i = 2, iip1
91              sinlon(i)=sin(rlonv(i))        coslon(i) = cos(rlonv(i))
92              coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi        sinlon(i) = sin(rlonv(i))
93              sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi        coslondlon(i) = coslon(i)*(rlonu(i)-rlonu(i-1))/pi
94              print*,coslondlon(i),sinlondlon(i)        sinlondlon(i) = sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
95           enddo        PRINT *, coslondlon(i), sinlondlon(i)
96           coslon(1)=coslon(iip1)      END DO
97           coslondlon(1)=coslondlon(iip1)      coslon(1) = coslon(iip1)
98           sinlon(1)=sinlon(iip1)      coslondlon(1) = coslondlon(iip1)
99           sinlondlon(1)=sinlondlon(iip1)      sinlon(1) = sinlon(iip1)
100           print*,'sum sinlondlon ',ssum(iim,sinlondlon,1)/sinlondlon(1)      sinlondlon(1) = sinlondlon(iip1)
101           print*,'sum coslondlon ',ssum(iim,coslondlon,1)/coslondlon(1)      PRINT *, 'sum sinlondlon ', ssum(iim, sinlondlon, 1)/sinlondlon(1)
102          DO l = 1,llm      PRINT *, 'sum coslondlon ', ssum(iim, coslondlon, 1)/coslondlon(1)
103          DO j = 1,jjp1      DO l = 1, llm
104           DO i = 1,iip1        DO j = 1, jjp1
105           q ( i,j,l,1 )=0.          DO i = 1, iip1
106           q ( i,j,l,2 )=0.            q(i, j, l, 1) = 0.
107           q ( i,j,l,3 )=0.              q(i, j, l, 2) = 0.
108           ENDDO            q(i, j, l, 3) = 0.
109           ENDDO          END DO
110          ENDDO        END DO
111                END DO
112        endif  
113      END IF
114  c *** q contient les qqtes de traceur avant l'advection  
115      ! *** q contient les qqtes de traceur avant l'advection
116  c *** Affectation des tableaux S a partir de Q  
117  c *** Rem : utilisation de SCOPY ulterieurement    ! *** Affectation des tableaux S a partir de Q
118      ! *** Rem : utilisation de SCOPY ulterieurement
119         DO l = 1,llm  
120          DO j = 1,jjp1    DO l = 1, llm
121           DO i = 1,iip1      DO j = 1, jjp1
122               s0( i,j,llm+1-l ) = q ( i,j,l,0 )        DO i = 1, iip1
123               sx( i,j,llm+1-l ) = q ( i,j,l,1 )          s0(i, j, llm+1-l) = q(i, j, l, 0)
124               sy( i,j,llm+1-l ) = q ( i,j,l,2 )          sx(i, j, llm+1-l) = q(i, j, l, 1)
125               sz( i,j,llm+1-l ) = q ( i,j,l,3 )          sy(i, j, llm+1-l) = q(i, j, l, 2)
126           ENDDO          sz(i, j, llm+1-l) = q(i, j, l, 3)
127          ENDDO        END DO
128         ENDDO      END DO
129      END DO
130  c *** On calcule la masse d'air en kg  
131      ! *** On calcule la masse d'air en kg
132         DO  l = 1,llm  
133           DO  j = 1,jjp1    DO l = 1, llm
134             DO  i = 1,iip1      DO j = 1, jjp1
135              sm ( i,j,llm+1-l)=masse( i,j,l )        DO i = 1, iip1
136            ENDDO          sm(i, j, llm+1-l) = masse(i, j, l)
137           ENDDO        END DO
138         ENDDO      END DO
139      END DO
140  c *** On converti les champs S en atome (resp. kg)  
141  c *** Les routines d'advection traitent les champs    ! *** On converti les champs S en atome (resp. kg)
142  c *** a advecter si ces derniers sont en atome (resp. kg)    ! *** Les routines d'advection traitent les champs
143  c *** A optimiser !!!    ! *** a advecter si ces derniers sont en atome (resp. kg)
144      ! *** A optimiser !!!
145         DO  l = 1,llm  
146           DO  j = 1,jjp1    DO l = 1, llm
147             DO  i = 1,iip1      DO j = 1, jjp1
148                 s0(i,j,l) = s0(i,j,l) * sm ( i,j,l )        DO i = 1, iip1
149                 sx(i,j,l) = sx(i,j,l) * sm ( i,j,l )          s0(i, j, l) = s0(i, j, l)*sm(i, j, l)
150                 sy(i,j,l) = sy(i,j,l) * sm ( i,j,l )          sx(i, j, l) = sx(i, j, l)*sm(i, j, l)
151                 sz(i,j,l) = sz(i,j,l) * sm ( i,j,l )          sy(i, j, l) = sy(i, j, l)*sm(i, j, l)
152             ENDDO          sz(i, j, l) = sz(i, j, l)*sm(i, j, l)
153           ENDDO        END DO
154         ENDDO      END DO
155      END DO
156  c *** Appel des subroutines d'advection en X, en Y et en Z  
157  c *** Advection avec "time-splitting"    ! *** Appel des subroutines d'advection en X, en Y et en Z
158            ! *** Advection avec "time-splitting"
159         if(mode.eq.2) then  
160            do l=1,llm    IF (mode==2) THEN
161              s0s=0.      DO l = 1, llm
162              s0n=0.        s0s = 0.
163              dyn1=0.        s0n = 0.
164              dys1=0.        dyn1 = 0.
165              dyn2=0.        dys1 = 0.
166              dys2=0.        dyn2 = 0.
167              smn=0.        dys2 = 0.
168              sms=0.        smn = 0.
169              do i=1,iim        sms = 0.
170                 smn=smn+sm(i,1,l)        DO i = 1, iim
171                 sms=sms+sm(i,jjp1,l)          smn = smn + sm(i, 1, l)
172                 s0n=s0n+s0(i,1,l)          sms = sms + sm(i, jjp1, l)
173                 s0s=s0s+s0(i,jjp1,l)          s0n = s0n + s0(i, 1, l)
174                 zz=sy(i,1,l)/sm(i,1,l)          s0s = s0s + s0(i, jjp1, l)
175                 dyn1=dyn1+sinlondlon(i)*zz          zz = sy(i, 1, l)/sm(i, 1, l)
176                 dyn2=dyn2+coslondlon(i)*zz          dyn1 = dyn1 + sinlondlon(i)*zz
177                 zz=sy(i,jjp1,l)/sm(i,jjp1,l)          dyn2 = dyn2 + coslondlon(i)*zz
178                 dys1=dys1+sinlondlon(i)*zz          zz = sy(i, jjp1, l)/sm(i, jjp1, l)
179                 dys2=dys2+coslondlon(i)*zz          dys1 = dys1 + sinlondlon(i)*zz
180              enddo          dys2 = dys2 + coslondlon(i)*zz
181              do i=1,iim        END DO
182                 sy(i,1,l)=dyn1*sinlon(i)+dyn2*coslon(i)        DO i = 1, iim
183                 sy(i,jjp1,l)=dys1*sinlon(i)+dys2*coslon(i)          sy(i, 1, l) = dyn1*sinlon(i) + dyn2*coslon(i)
184              enddo          sy(i, jjp1, l) = dys1*sinlon(i) + dys2*coslon(i)
185              do i=1,iim        END DO
186                 s0(i,1,l)=s0n/smn+sy(i,1,l)        DO i = 1, iim
187                 s0(i,jjp1,l)=s0s/sms-sy(i,jjp1,l)          s0(i, 1, l) = s0n/smn + sy(i, 1, l)
188              enddo          s0(i, jjp1, l) = s0s/sms - sy(i, jjp1, l)
189          END DO
190              s0(iip1,1,l)=s0(1,1,l)  
191              s0(iip1,jjp1,l)=s0(1,jjp1,l)        s0(iip1, 1, l) = s0(1, 1, l)
192          s0(iip1, jjp1, l) = s0(1, jjp1, l)
193              do i=1,iim  
194                 sxn(i)=s0(i+1,1,l)-s0(i,1,l)        DO i = 1, iim
195                 sxs(i)=s0(i+1,jjp1,l)-s0(i,jjp1,l)          sxn(i) = s0(i+1, 1, l) - s0(i, 1, l)
196  c   on rerentre les masses          sxs(i) = s0(i+1, jjp1, l) - s0(i, jjp1, l)
197              enddo          ! on rerentre les masses
198              do i=1,iim        END DO
199                 sy(i,1,l)=sy(i,1,l)*sm(i,1,l)        DO i = 1, iim
200                 sy(i,jjp1,l)=sy(i,jjp1,l)*sm(i,jjp1,l)          sy(i, 1, l) = sy(i, 1, l)*sm(i, 1, l)
201                 s0(i,1,l)=s0(i,1,l)*sm(i,1,l)          sy(i, jjp1, l) = sy(i, jjp1, l)*sm(i, jjp1, l)
202                 s0(i,jjp1,l)=s0(i,jjp1,l)*sm(i,jjp1,l)          s0(i, 1, l) = s0(i, 1, l)*sm(i, 1, l)
203              enddo          s0(i, jjp1, l) = s0(i, jjp1, l)*sm(i, jjp1, l)
204              sxn(iip1)=sxn(1)        END DO
205              sxs(iip1)=sxs(1)        sxn(iip1) = sxn(1)
206              do i=1,iim        sxs(iip1) = sxs(1)
207                 sx(i+1,1,l)=0.25*(sxn(i)+sxn(i+1))*sm(i+1,1,l)        DO i = 1, iim
208                 sx(i+1,jjp1,l)=0.25*(sxs(i)+sxs(i+1))*sm(i+1,jjp1,l)          sx(i+1, 1, l) = 0.25*(sxn(i)+sxn(i+1))*sm(i+1, 1, l)
209              enddo          sx(i+1, jjp1, l) = 0.25*(sxs(i)+sxs(i+1))*sm(i+1, jjp1, l)
210              s0(iip1,1,l)=s0(1,1,l)        END DO
211              s0(iip1,jjp1,l)=s0(1,jjp1,l)        s0(iip1, 1, l) = s0(1, 1, l)
212              sy(iip1,1,l)=sy(1,1,l)        s0(iip1, jjp1, l) = s0(1, jjp1, l)
213              sy(iip1,jjp1,l)=sy(1,jjp1,l)        sy(iip1, 1, l) = sy(1, 1, l)
214              sx(1,1,l)=sx(iip1,1,l)        sy(iip1, jjp1, l) = sy(1, jjp1, l)
215              sx(1,jjp1,l)=sx(iip1,jjp1,l)        sx(1, 1, l) = sx(iip1, 1, l)
216            enddo        sx(1, jjp1, l) = sx(iip1, jjp1, l)
217        endif      END DO
218      END IF
219        if (mode.eq.4) then  
220           do l=1,llm    IF (mode==4) THEN
221              do i=1,iip1      DO l = 1, llm
222                 sx(i,1,l)=0.        DO i = 1, iip1
223                 sx(i,jjp1,l)=0.          sx(i, 1, l) = 0.
224                 sy(i,1,l)=0.          sx(i, jjp1, l) = 0.
225                 sy(i,jjp1,l)=0.          sy(i, 1, l) = 0.
226              enddo          sy(i, jjp1, l) = 0.
227           enddo        END DO
228        endif      END DO
229        call limx(s0,sx,sm,pente_max)    END IF
230         call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)    CALL limx(s0, sx, sm, pente_max)
231        if (mode.eq.4) then    CALL advx(limit, .5*dtvr, pbaru, sm, s0, sx, sy, sz, lati, latf)
232           do l=1,llm    IF (mode==4) THEN
233              do i=1,iip1      DO l = 1, llm
234                 sx(i,1,l)=0.        DO i = 1, iip1
235                 sx(i,jjp1,l)=0.          sx(i, 1, l) = 0.
236                 sy(i,1,l)=0.          sx(i, jjp1, l) = 0.
237                 sy(i,jjp1,l)=0.          sy(i, 1, l) = 0.
238              enddo          sy(i, jjp1, l) = 0.
239           enddo        END DO
240        endif      END DO
241         call   limy(s0,sy,sm,pente_max)    END IF
242         call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )    CALL limy(s0, sy, sm, pente_max)
243         do j=1,jjp1    CALL advy(limit, .5*dtvr, pbarv, sm, s0, sx, sy, sz)
244            do i=1,iip1    DO j = 1, jjp1
245               sz(i,j,1)=0.      DO i = 1, iip1
246               sz(i,j,llm)=0.        sz(i, j, 1) = 0.
247            enddo        sz(i, j, llm) = 0.
248         enddo      END DO
249         call limz(s0,sz,sm,pente_max)    END DO
250         call advz( limit,dtvr,w,sm,s0,sx,sy,sz )    CALL limz(s0, sz, sm, pente_max)
251        if (mode.eq.4) then    CALL advz(limit, dtvr, w, sm, s0, sx, sy, sz)
252           do l=1,llm    IF (mode==4) THEN
253              do i=1,iip1      DO l = 1, llm
254                 sx(i,1,l)=0.        DO i = 1, iip1
255                 sx(i,jjp1,l)=0.          sx(i, 1, l) = 0.
256                 sy(i,1,l)=0.          sx(i, jjp1, l) = 0.
257                 sy(i,jjp1,l)=0.          sy(i, 1, l) = 0.
258              enddo          sy(i, jjp1, l) = 0.
259           enddo        END DO
260        endif      END DO
261          call limy(s0,sy,sm,pente_max)    END IF
262         call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )    CALL limy(s0, sy, sm, pente_max)
263         do l=1,llm    CALL advy(limit, .5*dtvr, pbarv, sm, s0, sx, sy, sz)
264            do j=1,jjp1    DO l = 1, llm
265               sm(iip1,j,l)=sm(1,j,l)      DO j = 1, jjp1
266               s0(iip1,j,l)=s0(1,j,l)        sm(iip1, j, l) = sm(1, j, l)
267               sx(iip1,j,l)=sx(1,j,l)        s0(iip1, j, l) = s0(1, j, l)
268               sy(iip1,j,l)=sy(1,j,l)        sx(iip1, j, l) = sx(1, j, l)
269               sz(iip1,j,l)=sz(1,j,l)        sy(iip1, j, l) = sy(1, j, l)
270            enddo        sz(iip1, j, l) = sz(1, j, l)
271         enddo      END DO
272      END DO
273    
274        if (mode.eq.4) then  
275           do l=1,llm    IF (mode==4) THEN
276              do i=1,iip1      DO l = 1, llm
277                 sx(i,1,l)=0.        DO i = 1, iip1
278                 sx(i,jjp1,l)=0.          sx(i, 1, l) = 0.
279                 sy(i,1,l)=0.          sx(i, jjp1, l) = 0.
280                 sy(i,jjp1,l)=0.          sy(i, 1, l) = 0.
281              enddo          sy(i, jjp1, l) = 0.
282           enddo        END DO
283        endif      END DO
284         call limx(s0,sx,sm,pente_max)    END IF
285         call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)    CALL limx(s0, sx, sm, pente_max)
286  c ***   On repasse les S dans la variable q directement 14/10/94    CALL advx(limit, .5*dtvr, pbaru, sm, s0, sx, sy, sz, lati, latf)
287  c   On revient a des rapports de melange en divisant par la masse    ! ***   On repasse les S dans la variable q directement 14/10/94
288      ! On revient a des rapports de melange en divisant par la masse
289  c En dehors des poles:  
290      ! En dehors des poles:
291         DO  l = 1,llm  
292          DO  j = 1,jjp1    DO l = 1, llm
293           DO  i = 1,iim      DO j = 1, jjp1
294               q(i,j,llm+1-l,0)=s0(i,j,l)/sm(i,j,l)        DO i = 1, iim
295               q(i,j,llm+1-l,1)=sx(i,j,l)/sm(i,j,l)          q(i, j, llm+1-l, 0) = s0(i, j, l)/sm(i, j, l)
296               q(i,j,llm+1-l,2)=sy(i,j,l)/sm(i,j,l)          q(i, j, llm+1-l, 1) = sx(i, j, l)/sm(i, j, l)
297               q(i,j,llm+1-l,3)=sz(i,j,l)/sm(i,j,l)          q(i, j, llm+1-l, 2) = sy(i, j, l)/sm(i, j, l)
298           ENDDO          q(i, j, llm+1-l, 3) = sz(i, j, l)/sm(i, j, l)
299          ENDDO        END DO
300        ENDDO      END DO
301      END DO
302  c Traitements specifiques au pole  
303      ! Traitements specifiques au pole
304        if(mode.ge.1) then  
305        DO l=1,llm    IF (mode>=1) THEN
306  c   filtrages aux poles      DO l = 1, llm
307           masn=ssum(iim,sm(1,1,l),1)        ! filtrages aux poles
308           mass=ssum(iim,sm(1,jjp1,l),1)        masn = ssum(iim, sm(1,1,l), 1)
309           qpn=ssum(iim,s0(1,1,l),1)/masn        mass = ssum(iim, sm(1,jjp1,l), 1)
310           qps=ssum(iim,s0(1,jjp1,l),1)/mass        qpn = ssum(iim, s0(1,1,l), 1)/masn
311           dqzpn=ssum(iim,sz(1,1,l),1)/masn        qps = ssum(iim, s0(1,jjp1,l), 1)/mass
312           dqzps=ssum(iim,sz(1,jjp1,l),1)/mass        dqzpn = ssum(iim, sz(1,1,l), 1)/masn
313           do i=1,iip1        dqzps = ssum(iim, sz(1,jjp1,l), 1)/mass
314              q( i,1,llm+1-l,3)=dqzpn        DO i = 1, iip1
315              q( i,jjp1,llm+1-l,3)=dqzps          q(i, 1, llm+1-l, 3) = dqzpn
316              q( i,1,llm+1-l,0)=qpn          q(i, jjp1, llm+1-l, 3) = dqzps
317              q( i,jjp1,llm+1-l,0)=qps          q(i, 1, llm+1-l, 0) = qpn
318           enddo          q(i, jjp1, llm+1-l, 0) = qps
319           if(mode.eq.3) then        END DO
320              dyn1=0.        IF (mode==3) THEN
321              dys1=0.          dyn1 = 0.
322              dyn2=0.          dys1 = 0.
323              dys2=0.          dyn2 = 0.
324              do i=1,iim          dys2 = 0.
325                 dyn1=dyn1+sinlondlon(i)*sy(i,1,l)/sm(i,1,l)          DO i = 1, iim
326                 dyn2=dyn2+coslondlon(i)*sy(i,1,l)/sm(i,1,l)            dyn1 = dyn1 + sinlondlon(i)*sy(i, 1, l)/sm(i, 1, l)
327                 dys1=dys1+sinlondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)            dyn2 = dyn2 + coslondlon(i)*sy(i, 1, l)/sm(i, 1, l)
328                 dys2=dys2+coslondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)            dys1 = dys1 + sinlondlon(i)*sy(i, jjp1, l)/sm(i, jjp1, l)
329              enddo            dys2 = dys2 + coslondlon(i)*sy(i, jjp1, l)/sm(i, jjp1, l)
330              do i=1,iim          END DO
331                 q(i,1,llm+1-l,2)=          DO i = 1, iim
332       s          (sinlon(i)*dyn1+coslon(i)*dyn2)            q(i, 1, llm+1-l, 2) = (sinlon(i)*dyn1+coslon(i)*dyn2)
333                 q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)            q(i, 1, llm+1-l, 0) = q(i, 1, llm+1-l, 0) + q(i, 1, llm+1-l, 2)
334                 q(i,jjp1,llm+1-l,2)=            q(i, jjp1, llm+1-l, 2) = (sinlon(i)*dys1+coslon(i)*dys2)
335       s          (sinlon(i)*dys1+coslon(i)*dys2)            q(i, jjp1, llm+1-l, 0) = q(i, jjp1, llm+1-l, 0) - &
336                 q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)              q(i, jjp1, llm+1-l, 2)
337       s         -q(i,jjp1,llm+1-l,2)          END DO
338              enddo        END IF
339           endif        IF (mode==1) THEN
340           if(mode.eq.1) then          ! on filtre les valeurs au bord de la "grande maille pole"
341  c   on filtre les valeurs au bord de la "grande maille pole"          dyn1 = 0.
342              dyn1=0.          dys1 = 0.
343              dys1=0.          dyn2 = 0.
344              dyn2=0.          dys2 = 0.
345              dys2=0.          DO i = 1, iim
346              do i=1,iim            zz = s0(i, 2, l)/sm(i, 2, l) - q(i, 1, llm+1-l, 0)
347                 zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)            dyn1 = dyn1 + sinlondlon(i)*zz
348                 dyn1=dyn1+sinlondlon(i)*zz            dyn2 = dyn2 + coslondlon(i)*zz
349                 dyn2=dyn2+coslondlon(i)*zz            zz = q(i, jjp1, llm+1-l, 0) - s0(i, jjm, l)/sm(i, jjm, l)
350                 zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)            dys1 = dys1 + sinlondlon(i)*zz
351                 dys1=dys1+sinlondlon(i)*zz            dys2 = dys2 + coslondlon(i)*zz
352                 dys2=dys2+coslondlon(i)*zz          END DO
353              enddo          DO i = 1, iim
354              do i=1,iim            q(i, 1, llm+1-l, 2) = (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
355                 q(i,1,llm+1-l,2)=            q(i, 1, llm+1-l, 0) = q(i, 1, llm+1-l, 0) + q(i, 1, llm+1-l, 2)
356       s          (sinlon(i)*dyn1+coslon(i)*dyn2)/2.            q(i, jjp1, llm+1-l, 2) = (sinlon(i)*dys1+coslon(i)*dys2)/2.
357                 q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)            q(i, jjp1, llm+1-l, 0) = q(i, jjp1, llm+1-l, 0) - &
358                 q(i,jjp1,llm+1-l,2)=              q(i, jjp1, llm+1-l, 2)
359       s          (sinlon(i)*dys1+coslon(i)*dys2)/2.          END DO
360                 q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)          q(iip1, 1, llm+1-l, 0) = q(1, 1, llm+1-l, 0)
361       s         -q(i,jjp1,llm+1-l,2)          q(iip1, jjp1, llm+1-l, 0) = q(1, jjp1, llm+1-l, 0)
362              enddo  
363              q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)          DO i = 1, iim
364              q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)            sxn(i) = q(i+1, 1, llm+1-l, 0) - q(i, 1, llm+1-l, 0)
365              sxs(i) = q(i+1, jjp1, llm+1-l, 0) - q(i, jjp1, llm+1-l, 0)
366              do i=1,iim          END DO
367                 sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)          sxn(iip1) = sxn(1)
368                 sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)          sxs(iip1) = sxs(1)
369              enddo          DO i = 1, iim
370              sxn(iip1)=sxn(1)            q(i+1, 1, llm+1-l, 1) = 0.25*(sxn(i)+sxn(i+1))
371              sxs(iip1)=sxs(1)            q(i+1, jjp1, llm+1-l, 1) = 0.25*(sxs(i)+sxs(i+1))
372              do i=1,iim          END DO
373                 q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))          q(1, 1, llm+1-l, 1) = q(iip1, 1, llm+1-l, 1)
374                 q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))          q(1, jjp1, llm+1-l, 1) = q(iip1, jjp1, llm+1-l, 1)
375              enddo  
376              q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)        END IF
377              q(1,jjp1,llm+1-l,1)=q(iip1,jjp1,llm+1-l,1)  
378        END DO
379           endif    END IF
380    
381         ENDDO    ! bouclage en longitude
382         endif    DO iq = 0, 3
383        DO l = 1, llm
384  c bouclage en longitude        DO j = 1, jjp1
385        do iq=0,3          q(iip1, j, l, iq) = q(1, j, l, iq)
386           do l=1,llm        END DO
387              do j=1,jjp1      END DO
388                 q(iip1,j,l,iq)=q(1,j,l,iq)    END DO
389              enddo  
390           enddo    DO l = 1, llm
391        enddo      DO j = 1, jjp1
392          DO i = 1, iip1
393          DO l = 1,llm          IF (q(i,j,l,0)<0.) THEN
394           DO j = 1,jjp1            q(i, j, l, 0) = 0.
395            DO i = 1,iip1          END IF
396                  IF (q(i,j,l,0).lt.0.)  THEN        END DO
397                       q(i,j,l,0)=0.      END DO
398                   ENDIF    END DO
399            ENDDO  
400           ENDDO    DO l = 1, llm
401          ENDDO      DO j = 1, jjp1
402          DO i = 1, iip1
403         do l=1,llm          IF (q(i,j,l,0)<qmin) PRINT *, 'apres pentes, s0(', i, ',', j, ',', l, &
404            do j=1,jjp1            ')=', q(i, j, l, 0)
405             do i=1,iip1        END DO
406               if(q(i,j,l,0).lt.qmin)      END DO
407       ,       print*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0)    END DO
408             enddo    RETURN
409            enddo  END SUBROUTINE pentes_ini
        enddo  
       RETURN  
       END  
Legend:
Removed from v.76  
changed lines
  Added in v.81
  ViewVC Help
Powered by ViewVC 1.1.21