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

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

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

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