LMDZ
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, 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 
pentes_ini
subroutine pentes_ini(q, w, masse, pbaru, pbarv, mode)
Definition: pentes_ini.F:5
ip1jmp1
!$Header llmm1 INTEGER ip1jmp1
Definition: paramet.h:14
limy
subroutine limy(s0, sy, sm, pente_max)
Definition: limy.F:5
advz
subroutine advz(limit, dtz, w, sm, s0, sx, sy, sz)
Definition: advz.F:5
llm
!$Id Turb_fcg_gcssold get_uvd hqturb_gcssold endif!large scale llm day day1 day day1 *dt_toga endif!time annee_ref dt_toga u_toga vq_toga w_prof vq_prof llm day day1 day day1 *dt_dice endif!time annee_ref dt_dice swup_dice vg_dice omega_dice tg_prof vg_profd w_profd omega_profd!do llm!print llm l llm
Definition: 1D_interp_cases.h:103
pi
!$Id mode_top_bound COMMON comconstr && pi
Definition: comconst.h:7
limx
subroutine limx(s0, sx, sm, pente_max)
Definition: limx.F:5
advx
subroutine advx(limit, dtx, pbaru, sm, s0, sx, sy, sz, lati, latf)
Definition: advx.F:6
ip1jm
!$Header llmm1 INTEGER ip1jm
Definition: paramet.h:14
false
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
jjp1
!$Header jjp1
Definition: paramet.h:14
rlonu
!$Header!CDK comgeom COMMON comgeom rlonu
Definition: comgeom.h:25
limz
subroutine limz(s0, sz, sm, pente_max)
Definition: limz.F:5
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
advy
subroutine advy(limit, dty, pbarv, sm, s0, sx, sy, sz)
Definition: advy.F:5
dtvr
!$Id mode_top_bound COMMON comconstr dtvr
Definition: comconst.h:7
iim
c c zjulian c cym CALL iim cym klev iim
Definition: ini_bilKP_ave.h:24
rlonv
!$Header!CDK comgeom COMMON comgeom rlonv
Definition: comgeom.h:25