source: LMDZ5/branches/testing/libf/dyn3d_common/prather.F @ 3946

Last change on this file since 3946 was 2641, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2593:2640 into testing branch
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.7 KB
RevLine 
[524]1!
2! $Header$
3!
4      SUBROUTINE prather (q,w,masse,pbaru,pbarv,nt,dt)
[2641]5     
6      USE comconst_mod, ONLY: pi
7     
[524]8      IMPLICIT NONE
9
10c=======================================================================
11c   Adaptation LMDZ:  A.Armengaud (LGGE)
12c   ----------------
13c
14c   ************************************************
15c   Transport des traceurs par la methode de prather
16c   Ref :
17c
18c   ************************************************
19c   q,w,pext,pbaru et pbarv : arguments d'entree  pour le s-pg
20c
21c=======================================================================
22
23
[2641]24      include "dimensions.h"
25      include "paramet.h"
26      include "comgeom2.h"
[524]27
28c   Arguments:
29c   ----------
30      INTEGER iq,nt
31      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
32      REAL masse(iip1,jjp1,llm)
33      REAL q( iip1,jjp1,llm,0:9)
34      REAL w( ip1jmp1,llm )
35      integer ordre,ilim
36
37c   Local:
38c   ------
39      LOGICAL limit
40      real zq(iip1,jjp1,llm)
41      REAL sm ( iip1,jjp1, llm )
42      REAL s0( iip1,jjp1,llm ),  sx( iip1,jjp1,llm )
43      REAL sy( iip1,jjp1,llm ),  sz( iip1,jjp1,llm )
44      REAL sxx( iip1,jjp1,llm)
45      REAL sxy( iip1,jjp1,llm)
46      REAL sxz( iip1,jjp1,llm)
47      REAL syy( iip1,jjp1,llm )
48      REAL syz( iip1,jjp1,llm )
49      REAL szz( iip1,jjp1,llm ),zz
50      INTEGER i,j,l,indice
51      real sxn(iip1),sxs(iip1)
52
53      real sinlon(iip1),sinlondlon(iip1)
54      real coslon(iip1),coslondlon(iip1)
55      real qmin,qmax
56      save qmin,qmax
57      save sinlon,coslon,sinlondlon,coslondlon
58      real dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps
59      real masn,mass
60c
61      REAL      SSUM
62      integer ismax,ismin
[1952]63      EXTERNAL  SSUM, ismin,ismax
[524]64      logical first
65      save first
66      EXTERNAL advxp,advyp,advzp
67
68
69      data first/.true./
70      data qmin,qmax/-1.e33,1.e33/
71
72
73c==========================================================================
74c==========================================================================
75c     MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
76c==========================================================================
77c==========================================================================
78      REAL dt
79c==========================================================================
80      limit = .TRUE.
81 
82      if(first) then
83         print*,'SCHEMA PRATHER'
84         first=.false.
85         do i=2,iip1
86            coslon(i)=cos(rlonv(i))
87            sinlon(i)=sin(rlonv(i))
88            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
89            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
90         enddo
91         coslon(1)=coslon(iip1)
92         coslondlon(1)=coslondlon(iip1)
93         sinlon(1)=sinlon(iip1)
94         sinlondlon(1)=sinlondlon(iip1)
95
96        DO l = 1,llm
97        DO j = 1,jjp1
98        DO i = 1,iip1
99        q( i,j,l,1 )=0.
100        q( i,j,l,2)=0.
101        q( i,j,l,3)=0.
102        q( i,j,l,4)=0.
103        q( i,j,l,5)=0.
104        q( i,j,l,6)=0.
105        q( i,j,l,7)=0.
106        q( i,j,l,8)=0.
107        q( i,j,l,9)=0.
108        ENDDO
109        ENDDO
110        ENDDO
111      endif
112c   Fin modif Fred
113
114c *** On calcule la masse d'air en kg
115
116       DO l = 1,llm
117        DO j = 1,jjp1
118         DO i = 1,iip1
119         sm( i,j,llm+1-l ) =masse(i,j,l)
120         ENDDO
121        ENDDO
122       ENDDO
123
124c *** q contient les qqtes de traceur avant l'advection
125
126c *** Affectation des tableaux S a partir de Q
127 
128       DO l = 1,llm
129        DO j = 1,jjp1
130         DO i = 1,iip1
131       s0( i,j,l) = q ( i,j,llm+1-l,0 )*sm(i,j,l)
132       sx( i,j,l) = q( i,j,llm+1-l,1 )*sm(i,j,l)
133       sy( i,j,l) = q( i,j,llm+1-l,2)*sm(i,j,l)
134       sz( i,j,l) = q( i,j,llm+1-l,3)*sm(i,j,l)
135       sxx( i,j,l) = q( i,j,llm+1-l,4)*sm(i,j,l)
136       sxy( i,j,l) = q( i,j,llm+1-l,5)*sm(i,j,l)
137       sxz( i,j,l) = q( i,j,llm+1-l,6)*sm(i,j,l)
138       syy( i,j,l) = q( i,j,llm+1-l,7)*sm(i,j,l)
139       syz( i,j,l) = q( i,j,llm+1-l,8)*sm(i,j,l)
140       szz( i,j,l) = q( i,j,llm+1-l,9)*sm(i,j,l)
141         ENDDO
142        ENDDO
143       ENDDO
144c *** Appel des subroutines d'advection en X, en Y et en Z
145c *** Advection avec "time-splitting"
146     
147c-----------------------------------------------------------
148       do indice =1,nt
149       call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz
150     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
151        end do
152        do l=1,llm
153        do i=1,iip1
154        sy(i,1,l)=0.
155        sy(i,jjp1,l)=0.
156        enddo
157        enddo
158c---------------------------------------------------------
159       call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz
160     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
161c---------------------------------------------------------
162
163c---------------------------------------------------------
164       do j=1,jjp1
165          do i=1,iip1
166             sz(i,j,1)=0.
167             sz(i,j,llm)=0.
168             sxz(i,j,1)=0.
169             sxz(i,j,llm)=0.
170             syz(i,j,1)=0.
171             syz(i,j,llm)=0.
172             szz(i,j,1)=0.
173             szz(i,j,llm)=0.
174          enddo
175       enddo
176       call advzp( limit,dt*nt,w,sm,s0,sx,sy,sz
177     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
178        do l=1,llm
179        do i=1,iip1
180        sy(i,1,l)=0.
181        sy(i,jjp1,l)=0.
182        enddo
183        enddo
184
185c---------------------------------------------------------
186
187c---------------------------------------------------------
188       call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz
189     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
190c---------------------------------------------------------
191       DO l = 1,llm
192        DO j = 1,jjp1
193             s0( iip1,j,l)=s0( 1,j,l )
194             sx( iip1,j,l)=sx( 1,j,l )
195             sy( iip1,j,l)=sy( 1,j,l )
196             sz( iip1,j,l)=sz( 1,j,l )
197             sxx( iip1,j,l)=sxx( 1,j,l )
198             sxy( iip1,j,l)=sxy( 1,j,l)
199             sxz( iip1,j,l)=sxz( 1,j,l )
200             syy( iip1,j,l)=syy( 1,j,l )
201             syz( iip1,j,l)=syz( 1,j,l)
202             szz( iip1,j,l)=szz( 1,j,l )
203        ENDDO
204       ENDDO
205       do indice=1,nt
206       call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz
207     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
208        end do
209c---------------------------------------------------------
210c---------------------------------------------------------
211c ***   On repasse les S dans la variable qpr
212c ***   On repasse les S dans la variable q directement 14/10/94
213
214       DO  l = 1,llm
215        DO  j = 1,jjp1
216         DO  i = 1,iip1
217      q( i,j,llm+1-l,0 )=s0( i,j,l )/sm(i,j,l)
218      q( i,j,llm+1-l,1 ) = sx( i,j,l )/sm(i,j,l)
219      q( i,j,llm+1-l,2 ) = sy( i,j,l )/sm(i,j,l)
220      q( i,j,llm+1-l,3 ) = sz( i,j,l )/sm(i,j,l)
221      q( i,j,llm+1-l,4 ) = sxx( i,j,l )/sm(i,j,l)
222      q( i,j,llm+1-l,5 ) = sxy( i,j,l )/sm(i,j,l)
223      q( i,j,llm+1-l,6 ) = sxz( i,j,l )/sm(i,j,l)
224      q( i,j,llm+1-l,7 ) = syy( i,j,l )/sm(i,j,l)
225      q( i,j,llm+1-l,8 ) = syz( i,j,l )/sm(i,j,l)
226      q( i,j,llm+1-l,9 ) = szz( i,j,l )/sm(i,j,l)
227      ENDDO
228      ENDDO
229      ENDDO
230
231c---------------------------------------------------------
232c      go to  777
233c   filtrages aux poles
234
235c Traitements specifiques au pole
236
237c   filtrages aux poles
238         DO l=1,llm
239c   filtrages aux poles
240         masn=ssum(iim,sm(1,1,l),1)
241         mass=ssum(iim,sm(1,jjp1,l),1)
242         qpn=ssum(iim,s0(1,1,l),1)/masn
243         qps=ssum(iim,s0(1,jjp1,l),1)/mass
244         dqzpn=ssum(iim,sz(1,1,l),1)/masn
245         dqzps=ssum(iim,sz(1,jjp1,l),1)/mass
246         do i=1,iip1
247          q( i,1,llm+1-l,3)=dqzpn
248          q( i,jjp1,llm+1-l,3)=dqzps
249          q( i,1,llm+1-l,0)=qpn
250          q( i,jjp1,llm+1-l,0)=qps
251         enddo
252c       enddo
253c         print*,'qpn',qpn,'qps',qps
254c          print*,'dqzpn',dqzpn,'dqzps',dqzps
255c       enddo
256           dyn1=0.
257           dys1=0.
258           dyn2=0.
259           dys2=0.
260        do i=1,iim
261        zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
262        dyn1=dyn1+sinlondlon(i)*zz
263        dyn2=dyn2+coslondlon(i)*zz
264        zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
265        dys1=dys1+sinlondlon(i)*zz
266        dys2=dys2+coslondlon(i)*zz
267        enddo
268         do i=1,iim
269         q(i,1,llm+1-l,2)=
270     $   (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
271         q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)
272     $          +q(i,1,llm+1-l,2)
273         q(i,jjp1,llm+1-l,2)=
274     $   (sinlon(i)*dys1+coslon(i)*dys2)/2.
275         q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
276     $      -q(i,jjp1,llm+1-l,2)
277         enddo
278      q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
279      q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
280      do i=1,iim
281      sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
282      sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
283      enddo
284      sxn(iip1)=sxn(1)
285      sxs(iip1)=sxs(1)
286      do i=1,iim
287      q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
288      q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
289      END DO
290      q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
291      q(1,jjp1,llm+1-l,1)=
292     $   q(iip1,jjp1,llm+1-l,1)
293        enddo
294         do l=1,llm
295           do i=1,iim
296            q( i,1,llm+1-l,4)=0.
297            q( i,jjp1,llm+1-l,4)=0.
298            q( i,1,llm+1-l,5)=0.
299            q( i,jjp1,llm+1-l,5)=0.
300            q( i,1,llm+1-l,6)=0.
301            q( i,jjp1,llm+1-l,6)=0.
302            q( i,1,llm+1-l,7)=0.
303            q( i,jjp1,llm+1-l,7)=0.
304            q( i,1,llm+1-l,8)=0.
305            q( i,jjp1,llm+1-l,8)=0.
306            q( i,1,llm+1-l,9)=0.
307            q( i,jjp1,llm+1-l,9)=0.
308          enddo
309         ENDDO
310
311777      continue
312c
313c   bouclage en longitude
314      do l=1,llm
315      do j=1,jjp1
316      q(iip1,j,l,0)=q(1,j,l,0)
317      q(iip1,j,llm+1-l,0)=q(1,j,llm+1-l,0)
318      q(iip1,j,llm+1-l,1)=q(1,j,llm+1-l,1)
319      q(iip1,j,llm+1-l,2)=q(1,j,llm+1-l,2)
320      q(iip1,j,llm+1-l,3)=q(1,j,llm+1-l,3)
321      q(iip1,j,llm+1-l,4)=q(1,j,llm+1-l,4)
322      q(iip1,j,llm+1-l,5)=q(1,j,llm+1-l,5)
323      q(iip1,j,llm+1-l,6)=q(1,j,llm+1-l,6)
324      q(iip1,j,llm+1-l,7)=q(1,j,llm+1-l,7)
325      q(iip1,j,llm+1-l,8)=q(1,j,llm+1-l,8)
326      q(iip1,j,llm+1-l,9)=q(1,j,llm+1-l,9)
327      enddo
328      enddo
329        DO l = 1,llm
330         DO j = 2,jjm
331           DO i = 1,iip1
332         IF (q(i,j,l,0).lt.0.)  THEN
333         PRINT*,'------------ BIP-----------'
334         PRINT*,'S0(',i,j,l,')=',q(i,j,l,0),
335     $          q(i,j-1,l,0)
336         PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
337         PRINT*,'SY(',i,j,l,')=',q(i,j,l,2),
338     $   q(i,j-1,l,2)   
339         PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
340c                    PRINT*,' PBL EN SORTIE D'' ADVZP'
341                     q(i,j,l,0)=0.
342c                  STOP
343               ENDIF
344           ENDDO
345         ENDDO
346         do j=1,jjp1,jjm
347         do i=1,iip1
348               IF (q(i,j,l,0).lt.0.)  THEN
349               PRINT*,'------------ BIP 2-----------'
350         PRINT*,'S0(',i,j,l,')=',q(i,j,l,0)
351         PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
352         PRINT*,'SY(',i,j,l,')=',q(i,j,l,2)
353         PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
354
355                     q(i,j,l,0)=0.
356c                  STOP
357               ENDIF
358         enddo
359         enddo
360        ENDDO
361      RETURN
362      END
Note: See TracBrowser for help on using the repository browser.