source: LMDZ6/branches/Amaury_dev/libf/dyn3d_common/prather.f90 @ 5209

Last change on this file since 5209 was 5159, checked in by abarral, 7 weeks ago

Put dimensions.h and paramet.h into modules
  • 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: 9.7 KB
Line 
1
2! $Header$
3
4SUBROUTINE prather(q,w,masse,pbaru,pbarv,nt,dt)
5
6  USE comconst_mod, ONLY: pi
7  USE lmdz_ssum_scopy, ONLY: ssum
8  USE lmdz_comgeom2
9
10  USE lmdz_dimensions, ONLY: iim, jjm, llm, ndm
11  USE lmdz_paramet
12  IMPLICIT NONE
13
14  !=======================================================================
15  !   Adaptation LMDZ:  A.Armengaud (LGGE)
16  !   ----------------
17
18  !   ************************************************
19  !   Transport des traceurs par la methode de prather
20  !   Ref :
21
22  !   ************************************************
23  !   q,w,pext,pbaru et pbarv : arguments d'entree  pour le s-pg
24
25  !=======================================================================
26
27
28
29
30
31  !   Arguments:
32  !   ----------
33  INTEGER :: iq,nt
34  REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
35  REAL :: masse(iip1,jjp1,llm)
36  REAL :: q( iip1,jjp1,llm,0:9)
37  REAL :: w( ip1jmp1,llm )
38  INTEGER :: ordre,ilim
39
40  !   Local:
41  !   ------
42  LOGICAL :: limit
43  REAL :: zq(iip1,jjp1,llm)
44  REAL :: sm ( iip1,jjp1, llm )
45  REAL :: s0( iip1,jjp1,llm ),  sx( iip1,jjp1,llm )
46  REAL :: sy( iip1,jjp1,llm ),  sz( iip1,jjp1,llm )
47  REAL :: sxx( iip1,jjp1,llm)
48  REAL :: sxy( iip1,jjp1,llm)
49  REAL :: sxz( iip1,jjp1,llm)
50  REAL :: syy( iip1,jjp1,llm )
51  REAL :: syz( iip1,jjp1,llm )
52  REAL :: szz( iip1,jjp1,llm ),zz
53  INTEGER :: i,j,l,indice
54  REAL :: sxn(iip1),sxs(iip1)
55
56  REAL :: sinlon(iip1),sinlondlon(iip1)
57  REAL :: coslon(iip1),coslondlon(iip1)
58  REAL :: qmin,qmax
59  save qmin,qmax
60  save sinlon,coslon,sinlondlon,coslondlon
61  REAL :: dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps
62  REAL :: masn,mass
63
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
73  !==========================================================================
74  !==========================================================================
75  ! MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
76  !==========================================================================
77  !==========================================================================
78  REAL :: dt
79  !==========================================================================
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
112  !   Fin modif Fred
113
114  ! *** 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
124  ! *** q contient les qqtes de traceur avant l'advection
125
126  ! *** 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
144  ! *** Appel des subroutines d'advection en X, en Y et en Z
145  ! *** Advection avec "time-splitting"
146
147  !-----------------------------------------------------------
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
158  !---------------------------------------------------------
159   CALL advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz &
160         ,sxx,sxy,sxz,syy,syz,szz,1 )
161  !---------------------------------------------------------
162
163  !---------------------------------------------------------
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
185  !---------------------------------------------------------
186
187  !---------------------------------------------------------
188   CALL advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz &
189         ,sxx,sxy,sxz,syy,syz,szz,1 )
190  !---------------------------------------------------------
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
209  !---------------------------------------------------------
210  !---------------------------------------------------------
211  ! ***   On repasse les S dans la variable qpr
212  ! ***   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
231  !---------------------------------------------------------
232   ! go to  777
233  !   filtrages aux poles
234
235  ! Traitements specifiques au pole
236
237  !   filtrages aux poles
238     DO l=1,llm
239  !   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
252    ! enddo
253    !   PRINT*,'qpn',qpn,'qps',qps
254    !    PRINT*,'dqzpn',dqzpn,'dqzps',dqzps
255    ! 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
312
313  !   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)<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)
340  !                 PRINT*,' PBL EN SORTIE D'' ADVZP'
341                 q(i,j,l,0)=0.
342               ! 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)<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.
356               ! STOP
357           ENDIF
358     enddo
359     enddo
360    ENDDO
361  RETURN
362END SUBROUTINE prather
Note: See TracBrowser for help on using the repository browser.