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

Last change on this file since 5153 was 5136, checked in by abarral, 4 months ago

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