source: LMDZ6/trunk/libf/dyn3d_common/prather.f90 @ 5281

Last change on this file since 5281 was 5281, checked in by abarral, 4 days ago

Turn 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.9 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, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
9          ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm
10IMPLICIT 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  !   Arguments:
26  !   ----------
27  INTEGER :: iq,nt
28  REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
29  REAL :: masse(iip1,jjp1,llm)
30  REAL :: q( iip1,jjp1,llm,0:9)
31  REAL :: w( ip1jmp1,llm )
32  integer :: ordre,ilim
33
34  !   Local:
35  !   ------
36  LOGICAL :: limit
37  real :: zq(iip1,jjp1,llm)
38  REAL :: sm ( iip1,jjp1, llm )
39  REAL :: s0( iip1,jjp1,llm ),  sx( iip1,jjp1,llm )
40  REAL :: sy( iip1,jjp1,llm ),  sz( iip1,jjp1,llm )
41  REAL :: sxx( iip1,jjp1,llm)
42  REAL :: sxy( iip1,jjp1,llm)
43  REAL :: sxz( iip1,jjp1,llm)
44  REAL :: syy( iip1,jjp1,llm )
45  REAL :: syz( iip1,jjp1,llm )
46  REAL :: szz( iip1,jjp1,llm ),zz
47  INTEGER :: i,j,l,indice
48  real :: sxn(iip1),sxs(iip1)
49
50  real :: sinlon(iip1),sinlondlon(iip1)
51  real :: coslon(iip1),coslondlon(iip1)
52  real :: qmin,qmax
53  save qmin,qmax
54  save sinlon,coslon,sinlondlon,coslondlon
55  real :: dyn1,dyn2,dys1,dys2,qpn,qps,dqzpn,dqzps
56  real :: masn,mass
57  !
58  REAL :: SSUM
59  integer :: ismax,ismin
60  EXTERNAL  SSUM, ismin,ismax
61  logical :: first
62  save first
63  EXTERNAL advxp,advyp,advzp
64
65
66  data first/.true./
67  data qmin,qmax/-1.e33,1.e33/
68
69
70  !==========================================================================
71  !==========================================================================
72  ! MODIFICATION POUR PAS DE TEMPS ADAPTATIF, dtvr remplace par dt
73  !==========================================================================
74  !==========================================================================
75  REAL :: dt
76  !==========================================================================
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
109  !   Fin modif Fred
110
111  ! *** 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
121  ! *** q contient les qqtes de traceur avant l'advection
122
123  ! *** 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
141  ! *** Appel des subroutines d'advection en X, en Y et en Z
142  ! *** Advection avec "time-splitting"
143
144  !-----------------------------------------------------------
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
155  !---------------------------------------------------------
156   call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz &
157         ,sxx,sxy,sxz,syy,syz,szz,1 )
158  !---------------------------------------------------------
159
160  !---------------------------------------------------------
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
182  !---------------------------------------------------------
183
184  !---------------------------------------------------------
185   call advyp( limit,.5*dt*nt,pbarv,sm,s0,sx,sy,sz &
186         ,sxx,sxy,sxz,syy,syz,szz,1 )
187  !---------------------------------------------------------
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
206  !---------------------------------------------------------
207  !---------------------------------------------------------
208  ! ***   On repasse les S dans la variable qpr
209  ! ***   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
228  !---------------------------------------------------------
229   ! go to  777
230  !   filtrages aux poles
231
232  ! Traitements specifiques au pole
233
234  !   filtrages aux poles
235     DO l=1,llm
236  !   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
249    ! enddo
250    !   print*,'qpn',qpn,'qps',qps
251    !    print*,'dqzpn',dqzpn,'dqzps',dqzps
252    ! 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
309  !
310  !   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)
337  !                  PRINT*,' PBL EN SORTIE D'' ADVZP'
338                 q(i,j,l,0)=0.
339               ! 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.
353               ! STOP
354           ENDIF
355     enddo
356     enddo
357    ENDDO
358  RETURN
359END SUBROUTINE prather
Note: See TracBrowser for help on using the repository browser.