source: LMDZ4/trunk/libf/dyn3dpar/prather.F @ 5440

Last change on this file since 5440 was 764, checked in by Laurent Fairhead, 18 years ago

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