source: LMDZ6/branches/contrails/libf/dyn3d_common/prather.f90 @ 5456

Last change on this file since 5456 was 5285, checked in by abarral, 2 months ago

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