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