source: trunk/libf/dyn3d/prather.F @ 1

Last change on this file since 1 was 1, checked in by emillour, 14 years ago

Import initial LMDZ5
File size: 10.6 KB
RevLine 
[1]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
62      EXTERNAL  SSUM, convflu,ismin,ismax
63      logical first
64      save first
65      EXTERNAL advxp,advyp,advzp
66
67
68      data first/.true./
69      data qmin,qmax/-1.e33,1.e33/
70
71
72c==========================================================================
73c==========================================================================
74c     MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
75c==========================================================================
76c==========================================================================
77      REAL dt
78c==========================================================================
79      limit = .TRUE.
80 
81      if(first) then
82         print*,'SCHEMA PRATHER'
83         first=.false.
84         do i=2,iip1
85            coslon(i)=cos(rlonv(i))
86            sinlon(i)=sin(rlonv(i))
87            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
88            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
89         enddo
90         coslon(1)=coslon(iip1)
91         coslondlon(1)=coslondlon(iip1)
92         sinlon(1)=sinlon(iip1)
93         sinlondlon(1)=sinlondlon(iip1)
94
95        DO l = 1,llm
96        DO j = 1,jjp1
97        DO i = 1,iip1
98        q( i,j,l,1 )=0.
99        q( i,j,l,2)=0.
100        q( i,j,l,3)=0.
101        q( i,j,l,4)=0.
102        q( i,j,l,5)=0.
103        q( i,j,l,6)=0.
104        q( i,j,l,7)=0.
105        q( i,j,l,8)=0.
106        q( i,j,l,9)=0.
107        ENDDO
108        ENDDO
109        ENDDO
110      endif
111c   Fin modif Fred
112
113c *** On calcule la masse d'air en kg
114
115       DO l = 1,llm
116        DO j = 1,jjp1
117         DO i = 1,iip1
118         sm( i,j,llm+1-l ) =masse(i,j,l)
119         ENDDO
120        ENDDO
121       ENDDO
122
123c *** q contient les qqtes de traceur avant l'advection
124
125c *** Affectation des tableaux S a partir de Q
126 
127       DO l = 1,llm
128        DO j = 1,jjp1
129         DO i = 1,iip1
130       s0( i,j,l) = q ( i,j,llm+1-l,0 )*sm(i,j,l)
131       sx( i,j,l) = q( i,j,llm+1-l,1 )*sm(i,j,l)
132       sy( i,j,l) = q( i,j,llm+1-l,2)*sm(i,j,l)
133       sz( i,j,l) = q( i,j,llm+1-l,3)*sm(i,j,l)
134       sxx( i,j,l) = q( i,j,llm+1-l,4)*sm(i,j,l)
135       sxy( i,j,l) = q( i,j,llm+1-l,5)*sm(i,j,l)
136       sxz( i,j,l) = q( i,j,llm+1-l,6)*sm(i,j,l)
137       syy( i,j,l) = q( i,j,llm+1-l,7)*sm(i,j,l)
138       syz( i,j,l) = q( i,j,llm+1-l,8)*sm(i,j,l)
139       szz( i,j,l) = q( i,j,llm+1-l,9)*sm(i,j,l)
140         ENDDO
141        ENDDO
142       ENDDO
143c *** Appel des subroutines d'advection en X, en Y et en Z
144c *** Advection avec "time-splitting"
145     
146c-----------------------------------------------------------
147       do indice =1,nt
148       call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz
149     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
150        end do
151        do l=1,llm
152        do i=1,iip1
153        sy(i,1,l)=0.
154        sy(i,jjp1,l)=0.
155        enddo
156        enddo
157c---------------------------------------------------------
158       call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz
159     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
160c---------------------------------------------------------
161
162c---------------------------------------------------------
163       do j=1,jjp1
164          do i=1,iip1
165             sz(i,j,1)=0.
166             sz(i,j,llm)=0.
167             sxz(i,j,1)=0.
168             sxz(i,j,llm)=0.
169             syz(i,j,1)=0.
170             syz(i,j,llm)=0.
171             szz(i,j,1)=0.
172             szz(i,j,llm)=0.
173          enddo
174       enddo
175       call advzp( limit,dt*nt,w,sm,s0,sx,sy,sz
176     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
177        do l=1,llm
178        do i=1,iip1
179        sy(i,1,l)=0.
180        sy(i,jjp1,l)=0.
181        enddo
182        enddo
183
184c---------------------------------------------------------
185
186c---------------------------------------------------------
187       call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz
188     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
189c---------------------------------------------------------
190       DO l = 1,llm
191        DO j = 1,jjp1
192             s0( iip1,j,l)=s0( 1,j,l )
193             sx( iip1,j,l)=sx( 1,j,l )
194             sy( iip1,j,l)=sy( 1,j,l )
195             sz( iip1,j,l)=sz( 1,j,l )
196             sxx( iip1,j,l)=sxx( 1,j,l )
197             sxy( iip1,j,l)=sxy( 1,j,l)
198             sxz( iip1,j,l)=sxz( 1,j,l )
199             syy( iip1,j,l)=syy( 1,j,l )
200             syz( iip1,j,l)=syz( 1,j,l)
201             szz( iip1,j,l)=szz( 1,j,l )
202        ENDDO
203       ENDDO
204       do indice=1,nt
205       call advxp( limit,0.5*dt,pbaru,sm,s0,sx,sy,sz
206     .             ,sxx,sxy,sxz,syy,syz,szz,1 )
207        end do
208c---------------------------------------------------------
209c---------------------------------------------------------
210c ***   On repasse les S dans la variable qpr
211c ***   On repasse les S dans la variable q directement 14/10/94
212
213       DO  l = 1,llm
214        DO  j = 1,jjp1
215         DO  i = 1,iip1
216      q( i,j,llm+1-l,0 )=s0( i,j,l )/sm(i,j,l)
217      q( i,j,llm+1-l,1 ) = sx( i,j,l )/sm(i,j,l)
218      q( i,j,llm+1-l,2 ) = sy( i,j,l )/sm(i,j,l)
219      q( i,j,llm+1-l,3 ) = sz( i,j,l )/sm(i,j,l)
220      q( i,j,llm+1-l,4 ) = sxx( i,j,l )/sm(i,j,l)
221      q( i,j,llm+1-l,5 ) = sxy( i,j,l )/sm(i,j,l)
222      q( i,j,llm+1-l,6 ) = sxz( i,j,l )/sm(i,j,l)
223      q( i,j,llm+1-l,7 ) = syy( i,j,l )/sm(i,j,l)
224      q( i,j,llm+1-l,8 ) = syz( i,j,l )/sm(i,j,l)
225      q( i,j,llm+1-l,9 ) = szz( i,j,l )/sm(i,j,l)
226      ENDDO
227      ENDDO
228      ENDDO
229
230c---------------------------------------------------------
231c      go to  777
232c   filtrages aux poles
233
234c Traitements specifiques au pole
235
236c   filtrages aux poles
237         DO l=1,llm
238c   filtrages aux poles
239         masn=ssum(iim,sm(1,1,l),1)
240         mass=ssum(iim,sm(1,jjp1,l),1)
241         qpn=ssum(iim,s0(1,1,l),1)/masn
242         qps=ssum(iim,s0(1,jjp1,l),1)/mass
243         dqzpn=ssum(iim,sz(1,1,l),1)/masn
244         dqzps=ssum(iim,sz(1,jjp1,l),1)/mass
245         do i=1,iip1
246          q( i,1,llm+1-l,3)=dqzpn
247          q( i,jjp1,llm+1-l,3)=dqzps
248          q( i,1,llm+1-l,0)=qpn
249          q( i,jjp1,llm+1-l,0)=qps
250         enddo
251c       enddo
252c         print*,'qpn',qpn,'qps',qps
253c          print*,'dqzpn',dqzpn,'dqzps',dqzps
254c       enddo
255           dyn1=0.
256           dys1=0.
257           dyn2=0.
258           dys2=0.
259        do i=1,iim
260        zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
261        dyn1=dyn1+sinlondlon(i)*zz
262        dyn2=dyn2+coslondlon(i)*zz
263        zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
264        dys1=dys1+sinlondlon(i)*zz
265        dys2=dys2+coslondlon(i)*zz
266        enddo
267         do i=1,iim
268         q(i,1,llm+1-l,2)=
269     $   (sinlon(i)*dyn1+coslon(i)*dyn2)/2.
270         q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)
271     $          +q(i,1,llm+1-l,2)
272         q(i,jjp1,llm+1-l,2)=
273     $   (sinlon(i)*dys1+coslon(i)*dys2)/2.
274         q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
275     $      -q(i,jjp1,llm+1-l,2)
276         enddo
277      q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
278      q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
279      do i=1,iim
280      sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
281      sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
282      enddo
283      sxn(iip1)=sxn(1)
284      sxs(iip1)=sxs(1)
285      do i=1,iim
286      q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
287      q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
288      END DO
289      q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
290      q(1,jjp1,llm+1-l,1)=
291     $   q(iip1,jjp1,llm+1-l,1)
292        enddo
293         do l=1,llm
294           do i=1,iim
295            q( i,1,llm+1-l,4)=0.
296            q( i,jjp1,llm+1-l,4)=0.
297            q( i,1,llm+1-l,5)=0.
298            q( i,jjp1,llm+1-l,5)=0.
299            q( i,1,llm+1-l,6)=0.
300            q( i,jjp1,llm+1-l,6)=0.
301            q( i,1,llm+1-l,7)=0.
302            q( i,jjp1,llm+1-l,7)=0.
303            q( i,1,llm+1-l,8)=0.
304            q( i,jjp1,llm+1-l,8)=0.
305            q( i,1,llm+1-l,9)=0.
306            q( i,jjp1,llm+1-l,9)=0.
307          enddo
308         ENDDO
309
310777      continue
311c
312c   bouclage en longitude
313      do l=1,llm
314      do j=1,jjp1
315      q(iip1,j,l,0)=q(1,j,l,0)
316      q(iip1,j,llm+1-l,0)=q(1,j,llm+1-l,0)
317      q(iip1,j,llm+1-l,1)=q(1,j,llm+1-l,1)
318      q(iip1,j,llm+1-l,2)=q(1,j,llm+1-l,2)
319      q(iip1,j,llm+1-l,3)=q(1,j,llm+1-l,3)
320      q(iip1,j,llm+1-l,4)=q(1,j,llm+1-l,4)
321      q(iip1,j,llm+1-l,5)=q(1,j,llm+1-l,5)
322      q(iip1,j,llm+1-l,6)=q(1,j,llm+1-l,6)
323      q(iip1,j,llm+1-l,7)=q(1,j,llm+1-l,7)
324      q(iip1,j,llm+1-l,8)=q(1,j,llm+1-l,8)
325      q(iip1,j,llm+1-l,9)=q(1,j,llm+1-l,9)
326      enddo
327      enddo
328        DO l = 1,llm
329         DO j = 2,jjm
330           DO i = 1,iip1
331         IF (q(i,j,l,0).lt.0.)  THEN
332         PRINT*,'------------ BIP-----------'
333         PRINT*,'S0(',i,j,l,')=',q(i,j,l,0),
334     $          q(i,j-1,l,0)
335         PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
336         PRINT*,'SY(',i,j,l,')=',q(i,j,l,2),
337     $   q(i,j-1,l,2)   
338         PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
339c                    PRINT*,' PBL EN SORTIE D'' ADVZP'
340                     q(i,j,l,0)=0.
341c                  STOP
342               ENDIF
343           ENDDO
344         ENDDO
345         do j=1,jjp1,jjm
346         do i=1,iip1
347               IF (q(i,j,l,0).lt.0.)  THEN
348               PRINT*,'------------ BIP 2-----------'
349         PRINT*,'S0(',i,j,l,')=',q(i,j,l,0)
350         PRINT*,'SX(',i,j,l,')=',q(i,j,l,1)
351         PRINT*,'SY(',i,j,l,')=',q(i,j,l,2)
352         PRINT*,'SZ(',i,j,l,')=',q(i,j,l,3)
353
354                     q(i,j,l,0)=0.
355c                  STOP
356               ENDIF
357         enddo
358         enddo
359        ENDDO
360      RETURN
361      END
Note: See TracBrowser for help on using the repository browser.