My Project
 All Classes Files Functions Variables Macros
pentes_ini.F
Go to the documentation of this file.
1 !
2 ! $Header$
3 !
4  SUBROUTINE pentes_ini (q,w,masse,pbaru,pbarv,mode)
5  IMPLICIT NONE
6 
7 c=======================================================================
8 c Adaptation LMDZ: A.Armengaud (LGGE)
9 c ----------------
10 c
11 c ********************************************************************
12 c Transport des traceurs par la methode des pentes
13 c ********************************************************************
14 c Reference possible : Russel. G.L., Lerner J.A.:
15 c A new Finite-Differencing Scheme for Traceur Transport
16 c Equation , Journal of Applied Meteorology, pp 1483-1498,dec. 81
17 c ********************************************************************
18 c q,w,masse,pbaru et pbarv
19 c sont des arguments d'entree pour le s-pg ....
20 c
21 c=======================================================================
22 
23 
24 #include "dimensions.h"
25 #include "paramet.h"
26 #include "comconst.h"
27 #include "comvert.h"
28 #include "comgeom2.h"
29 
30 c Arguments:
31 c ----------
32  integer mode
33  REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm )
34  REAL q( iip1,jjp1,llm,0:3)
35  REAL w( ip1jmp1,llm )
36  REAL masse( iip1,jjp1,llm)
37 c Local:
38 c ------
39  LOGICAL limit
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 masn,mass,zz
44  INTEGER i,j,l,iq
45 
46 c modif Fred 24 03 96
47 
48  real sinlon(iip1),sinlondlon(iip1)
49  real coslon(iip1),coslondlon(iip1)
50  save sinlon,coslon,sinlondlon,coslondlon
51  real dyn1,dyn2,dys1,dys2
52  real qpn,qps,dqzpn,dqzps
53  real smn,sms,s0n,s0s,sxn(iip1),sxs(iip1)
54  real qmin,zq,pente_max
55 c
56  REAL ssum
57  integer ismax,ismin,lati,latf
58  EXTERNAL ssum, convflu,ismin,ismax
59  logical first
60  save first
61 c fin modif
62 
63 c EXTERNAL masskg
64  EXTERNAL advx
65  EXTERNAL advy
66  EXTERNAL advz
67 
68 c modif Fred 24 03 96
69  data first/.true./
70 
71  limit = .true.
72  pente_max=2
73 c if (mode.eq.1.or.mode.eq.3) then
74 c if (mode.eq.1) then
75  if (mode.ge.1) then
76  lati=2
77  latf=jjm
78  else
79  lati=1
80  latf=jjp1
81  endif
82 
83  qmin=0.4995
84  qmin=0.
85  if(first) then
86  print*,'SCHEMA AMONT NOUVEAU'
87  first=.false.
88  do i=2,iip1
89  coslon(i)=cos(rlonv(i))
90  sinlon(i)=sin(rlonv(i))
91  coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
92  sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
93  print*,coslondlon(i),sinlondlon(i)
94  enddo
95  coslon(1)=coslon(iip1)
96  coslondlon(1)=coslondlon(iip1)
97  sinlon(1)=sinlon(iip1)
98  sinlondlon(1)=sinlondlon(iip1)
99  print*,'sum sinlondlon ',ssum(iim,sinlondlon,1)/sinlondlon(1)
100  print*,'sum coslondlon ',ssum(iim,coslondlon,1)/coslondlon(1)
101  DO l = 1,llm
102  DO j = 1,jjp1
103  DO i = 1,iip1
104  q( i,j,l,1 )=0.
105  q( i,j,l,2 )=0.
106  q( i,j,l,3 )=0.
107  ENDDO
108  ENDDO
109  ENDDO
110 
111  endif
112 c Fin modif Fred
113 
114 c *** q contient les qqtes de traceur avant l'advection
115 
116 c *** Affectation des tableaux S a partir de Q
117 c *** Rem : utilisation de SCOPY ulterieurement
118 
119  DO l = 1,llm
120  DO j = 1,jjp1
121  DO i = 1,iip1
122  s0( i,j,llm+1-l ) = q( i,j,l,0 )
123  sx( i,j,llm+1-l ) = q( i,j,l,1 )
124  sy( i,j,llm+1-l ) = q( i,j,l,2 )
125  sz( i,j,llm+1-l ) = q( i,j,l,3 )
126  ENDDO
127  ENDDO
128  ENDDO
129 
130 c PRINT*,'----- S0 just before conversion -------'
131 c PRINT*,'S0(16,12,1)=',s0(16,12,1)
132 c PRINT*,'Q(16,12,1,4)=',q(16,12,1,4)
133 
134 c *** On calcule la masse d'air en kg
135 
136  DO l = 1,llm
137  DO j = 1,jjp1
138  DO i = 1,iip1
139  sm( i,j,llm+1-l)=masse( i,j,l )
140  ENDDO
141  ENDDO
142  ENDDO
143 
144 c *** On converti les champs S en atome (resp. kg)
145 c *** Les routines d'advection traitent les champs
146 c *** a advecter si ces derniers sont en atome (resp. kg)
147 c *** A optimiser !!!
148 
149  DO l = 1,llm
150  DO j = 1,jjp1
151  DO i = 1,iip1
152  s0(i,j,l) = s0(i,j,l) * sm( i,j,l )
153  sx(i,j,l) = sx(i,j,l) * sm( i,j,l )
154  sy(i,j,l) = sy(i,j,l) * sm( i,j,l )
155  sz(i,j,l) = sz(i,j,l) * sm( i,j,l )
156  ENDDO
157  ENDDO
158  ENDDO
159 
160 c ss0 = 0.
161 c DO l = 1,llm
162 c DO j = 1,jjp1
163 c DO i = 1,iim
164 c ss0 = ss0 + s0 ( i,j,l )
165 c ENDDO
166 c ENDDO
167 c ENDDO
168 c PRINT*, 'valeur tot s0 avant advection=',ss0
169 
170 c *** Appel des subroutines d'advection en X, en Y et en Z
171 c *** Advection avec "time-splitting"
172 
173 c-----------------------------------------------------------
174 c PRINT*,'----- S0 just before ADVX -------'
175 c PRINT*,'S0(16,12,1)=',s0(16,12,1)
176 
177 c-----------------------------------------------------------
178 c do l=1,llm
179 c do j=1,jjp1
180 c do i=1,iip1
181 c zq=s0(i,j,l)/sm(i,j,l)
182 c if(zq.lt.qmin)
183 c , print*,'avant advx1, s0(',i,',',j,',',l,')=',zq
184 c enddo
185 c enddo
186 c enddo
187 CCC
188  if(mode.eq.2) then
189  do l=1,llm
190  s0s=0.
191  s0n=0.
192  dyn1=0.
193  dys1=0.
194  dyn2=0.
195  dys2=0.
196  smn=0.
197  sms=0.
198  do i=1,iim
199  smn=smn+sm(i,1,l)
200  sms=sms+sm(i,jjp1,l)
201  s0n=s0n+s0(i,1,l)
202  s0s=s0s+s0(i,jjp1,l)
203  zz=sy(i,1,l)/sm(i,1,l)
204  dyn1=dyn1+sinlondlon(i)*zz
205  dyn2=dyn2+coslondlon(i)*zz
206  zz=sy(i,jjp1,l)/sm(i,jjp1,l)
207  dys1=dys1+sinlondlon(i)*zz
208  dys2=dys2+coslondlon(i)*zz
209  enddo
210  do i=1,iim
211  sy(i,1,l)=dyn1*sinlon(i)+dyn2*coslon(i)
212  sy(i,jjp1,l)=dys1*sinlon(i)+dys2*coslon(i)
213  enddo
214  do i=1,iim
215  s0(i,1,l)=s0n/smn+sy(i,1,l)
216  s0(i,jjp1,l)=s0s/sms-sy(i,jjp1,l)
217  enddo
218 
219  s0(iip1,1,l)=s0(1,1,l)
220  s0(iip1,jjp1,l)=s0(1,jjp1,l)
221 
222  do i=1,iim
223  sxn(i)=s0(i+1,1,l)-s0(i,1,l)
224  sxs(i)=s0(i+1,jjp1,l)-s0(i,jjp1,l)
225 c on rerentre les masses
226  enddo
227  do i=1,iim
228  sy(i,1,l)=sy(i,1,l)*sm(i,1,l)
229  sy(i,jjp1,l)=sy(i,jjp1,l)*sm(i,jjp1,l)
230  s0(i,1,l)=s0(i,1,l)*sm(i,1,l)
231  s0(i,jjp1,l)=s0(i,jjp1,l)*sm(i,jjp1,l)
232  enddo
233  sxn(iip1)=sxn(1)
234  sxs(iip1)=sxs(1)
235  do i=1,iim
236  sx(i+1,1,l)=0.25*(sxn(i)+sxn(i+1))*sm(i+1,1,l)
237  sx(i+1,jjp1,l)=0.25*(sxs(i)+sxs(i+1))*sm(i+1,jjp1,l)
238  enddo
239  s0(iip1,1,l)=s0(1,1,l)
240  s0(iip1,jjp1,l)=s0(1,jjp1,l)
241  sy(iip1,1,l)=sy(1,1,l)
242  sy(iip1,jjp1,l)=sy(1,jjp1,l)
243  sx(1,1,l)=sx(iip1,1,l)
244  sx(1,jjp1,l)=sx(iip1,jjp1,l)
245  enddo
246  endif
247 
248  if (mode.eq.4) then
249  do l=1,llm
250  do i=1,iip1
251  sx(i,1,l)=0.
252  sx(i,jjp1,l)=0.
253  sy(i,1,l)=0.
254  sy(i,jjp1,l)=0.
255  enddo
256  enddo
257  endif
258  call limx(s0,sx,sm,pente_max)
259 c call minmaxq(zq,1.e33,-1.e33,'avant advx ')
260  call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
261 c call minmaxq(zq,1.e33,-1.e33,'avant advy ')
262  if (mode.eq.4) then
263  do l=1,llm
264  do i=1,iip1
265  sx(i,1,l)=0.
266  sx(i,jjp1,l)=0.
267  sy(i,1,l)=0.
268  sy(i,jjp1,l)=0.
269  enddo
270  enddo
271  endif
272  call limy(s0,sy,sm,pente_max)
273  call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
274 c call minmaxq(zq,1.e33,-1.e33,'avant advz ')
275  do j=1,jjp1
276  do i=1,iip1
277  sz(i,j,1)=0.
278  sz(i,j,llm)=0.
279  enddo
280  enddo
281  call limz(s0,sz,sm,pente_max)
282  call advz( limit,dtvr,w,sm,s0,sx,sy,sz )
283  if (mode.eq.4) then
284  do l=1,llm
285  do i=1,iip1
286  sx(i,1,l)=0.
287  sx(i,jjp1,l)=0.
288  sy(i,1,l)=0.
289  sy(i,jjp1,l)=0.
290  enddo
291  enddo
292  endif
293  call limy(s0,sy,sm,pente_max)
294  call advy( limit,.5*dtvr,pbarv,sm,s0,sx,sy,sz )
295  do l=1,llm
296  do j=1,jjp1
297  sm(iip1,j,l)=sm(1,j,l)
298  s0(iip1,j,l)=s0(1,j,l)
299  sx(iip1,j,l)=sx(1,j,l)
300  sy(iip1,j,l)=sy(1,j,l)
301  sz(iip1,j,l)=sz(1,j,l)
302  enddo
303  enddo
304 
305 
306 c call minmaxq(zq,1.e33,-1.e33,'avant advx ')
307  if (mode.eq.4) then
308  do l=1,llm
309  do i=1,iip1
310  sx(i,1,l)=0.
311  sx(i,jjp1,l)=0.
312  sy(i,1,l)=0.
313  sy(i,jjp1,l)=0.
314  enddo
315  enddo
316  endif
317  call limx(s0,sx,sm,pente_max)
318  call advx( limit,.5*dtvr,pbaru,sm,s0,sx,sy,sz,lati,latf)
319 c call minmaxq(zq,1.e33,-1.e33,'apres advx ')
320 c do l=1,llm
321 c do j=1,jjp1
322 c do i=1,iip1
323 c zq=s0(i,j,l)/sm(i,j,l)
324 c if(zq.lt.qmin)
325 c , print*,'apres advx2, s0(',i,',',j,',',l,')=',zq
326 c enddo
327 c enddo
328 c enddo
329 c *** On repasse les S dans la variable q directement 14/10/94
330 c On revient a des rapports de melange en divisant par la masse
331 
332 c En dehors des poles:
333 
334  DO l = 1,llm
335  DO j = 1,jjp1
336  DO i = 1,iim
337  q(i,j,llm+1-l,0)=s0(i,j,l)/sm(i,j,l)
338  q(i,j,llm+1-l,1)=sx(i,j,l)/sm(i,j,l)
339  q(i,j,llm+1-l,2)=sy(i,j,l)/sm(i,j,l)
340  q(i,j,llm+1-l,3)=sz(i,j,l)/sm(i,j,l)
341  ENDDO
342  ENDDO
343  ENDDO
344 
345 c Traitements specifiques au pole
346 
347  if(mode.ge.1) then
348  DO l=1,llm
349 c filtrages aux poles
350  masn=ssum(iim,sm(1,1,l),1)
351  mass=ssum(iim,sm(1,jjp1,l),1)
352  qpn=ssum(iim,s0(1,1,l),1)/masn
353  qps=ssum(iim,s0(1,jjp1,l),1)/mass
354  dqzpn=ssum(iim,sz(1,1,l),1)/masn
355  dqzps=ssum(iim,sz(1,jjp1,l),1)/mass
356  do i=1,iip1
357  q( i,1,llm+1-l,3)=dqzpn
358  q( i,jjp1,llm+1-l,3)=dqzps
359  q( i,1,llm+1-l,0)=qpn
360  q( i,jjp1,llm+1-l,0)=qps
361  enddo
362  if(mode.eq.3) then
363  dyn1=0.
364  dys1=0.
365  dyn2=0.
366  dys2=0.
367  do i=1,iim
368  dyn1=dyn1+sinlondlon(i)*sy(i,1,l)/sm(i,1,l)
369  dyn2=dyn2+coslondlon(i)*sy(i,1,l)/sm(i,1,l)
370  dys1=dys1+sinlondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
371  dys2=dys2+coslondlon(i)*sy(i,jjp1,l)/sm(i,jjp1,l)
372  enddo
373  do i=1,iim
374  q(i,1,llm+1-l,2)=
375  s(sinlon(i)*dyn1+coslon(i)*dyn2)
376  q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
377  q(i,jjp1,llm+1-l,2)=
378  s(sinlon(i)*dys1+coslon(i)*dys2)
379  q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
380  s -q(i,jjp1,llm+1-l,2)
381  enddo
382  endif
383  if(mode.eq.1) then
384 c on filtre les valeurs au bord de la "grande maille pole"
385  dyn1=0.
386  dys1=0.
387  dyn2=0.
388  dys2=0.
389  do i=1,iim
390  zz=s0(i,2,l)/sm(i,2,l)-q(i,1,llm+1-l,0)
391  dyn1=dyn1+sinlondlon(i)*zz
392  dyn2=dyn2+coslondlon(i)*zz
393  zz=q(i,jjp1,llm+1-l,0)-s0(i,jjm,l)/sm(i,jjm,l)
394  dys1=dys1+sinlondlon(i)*zz
395  dys2=dys2+coslondlon(i)*zz
396  enddo
397  do i=1,iim
398  q(i,1,llm+1-l,2)=
399  s(sinlon(i)*dyn1+coslon(i)*dyn2)/2.
400  q(i,1,llm+1-l,0)=q(i,1,llm+1-l,0)+q(i,1,llm+1-l,2)
401  q(i,jjp1,llm+1-l,2)=
402  s(sinlon(i)*dys1+coslon(i)*dys2)/2.
403  q(i,jjp1,llm+1-l,0)=q(i,jjp1,llm+1-l,0)
404  s -q(i,jjp1,llm+1-l,2)
405  enddo
406  q(iip1,1,llm+1-l,0)=q(1,1,llm+1-l,0)
407  q(iip1,jjp1,llm+1-l,0)=q(1,jjp1,llm+1-l,0)
408 
409  do i=1,iim
410  sxn(i)=q(i+1,1,llm+1-l,0)-q(i,1,llm+1-l,0)
411  sxs(i)=q(i+1,jjp1,llm+1-l,0)-q(i,jjp1,llm+1-l,0)
412  enddo
413  sxn(iip1)=sxn(1)
414  sxs(iip1)=sxs(1)
415  do i=1,iim
416  q(i+1,1,llm+1-l,1)=0.25*(sxn(i)+sxn(i+1))
417  q(i+1,jjp1,llm+1-l,1)=0.25*(sxs(i)+sxs(i+1))
418  enddo
419  q(1,1,llm+1-l,1)=q(iip1,1,llm+1-l,1)
420  q(1,jjp1,llm+1-l,1)=q(iip1,jjp1,llm+1-l,1)
421 
422  endif
423 
424  ENDDO
425  endif
426 
427 c bouclage en longitude
428  do iq=0,3
429  do l=1,llm
430  do j=1,jjp1
431  q(iip1,j,l,iq)=q(1,j,l,iq)
432  enddo
433  enddo
434  enddo
435 
436 c PRINT*, ' SORTIE DE PENTES --- ca peut glisser ....'
437 
438  DO l = 1,llm
439  DO j = 1,jjp1
440  DO i = 1,iip1
441  IF (q(i,j,l,0).lt.0.) THEN
442 c PRINT*,'------------ BIP-----------'
443 c PRINT*,'Q0(',i,j,l,')=',q(i,j,l,0)
444 c PRINT*,'QX(',i,j,l,')=',q(i,j,l,1)
445 c PRINT*,'QY(',i,j,l,')=',q(i,j,l,2)
446 c PRINT*,'QZ(',i,j,l,')=',q(i,j,l,3)
447 c PRINT*,' PBL EN SORTIE DE PENTES'
448  q(i,j,l,0)=0.
449 c STOP
450  ENDIF
451  ENDDO
452  ENDDO
453  ENDDO
454 
455 c PRINT*, '-------------------------------------------'
456 
457  do l=1,llm
458  do j=1,jjp1
459  do i=1,iip1
460  if(q(i,j,l,0).lt.qmin)
461  , print*,'apres pentes, s0(',i,',',j,',',l,')=',q(i,j,l,0)
462  enddo
463  enddo
464  enddo
465  RETURN
466  END
467 
468 
469 
470 
471 
472 
473 
474 
475 
476 
477 
478