My Project
 All Classes Files Functions Variables Macros
vlsplt.F
Go to the documentation of this file.
1 c
2 c $Id: vlsplt.F 1550 2011-07-05 09:44:55Z lguez $
3 c
4 
5  SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt)
6 c
7 c Auteurs: P.Le Van, F.Hourdin, F.Forget
8 c
9 c ********************************************************************
10 c Shema d'advection " pseudo amont " .
11 c ********************************************************************
12 c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
13 c
14 c pente_max facteur de limitation des pentes: 2 en general
15 c 0 pour un schema amont
16 c pbaru,pbarv,w flux de masse en u ,v ,w
17 c pdt pas de temps
18 c
19 c --------------------------------------------------------------------
20  IMPLICIT NONE
21 c
22 #include "dimensions.h"
23 #include "paramet.h"
24 #include "logic.h"
25 #include "comvert.h"
26 #include "comconst.h"
27 
28 c
29 c Arguments:
30 c ----------
31  REAL masse(ip1jmp1,llm),pente_max
32 c REAL masse(iip1,jjp1,llm),pente_max
33  REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
34  REAL q(ip1jmp1,llm)
35 c REAL q(iip1,jjp1,llm)
36  REAL w(ip1jmp1,llm),pdt
37 c
38 c Local
39 c ---------
40 c
41  INTEGER i,ij,l,j,ii
42  INTEGER ijlqmin,iqmin,jqmin,lqmin
43 c
44  REAL zm(ip1jmp1,llm),newmasse
45  REAL mu(ip1jmp1,llm)
46  REAL mv(ip1jm,llm)
47  REAL mw(ip1jmp1,llm+1)
48  REAL zq(ip1jmp1,llm),zz
49  REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
50  REAL second,temps0,temps1,temps2,temps3
51  REAL ztemps1,ztemps2,ztemps3
52  REAL zzpbar, zzw
53  LOGICAL testcpu
54  SAVE testcpu
55  SAVE temps1,temps2,temps3
56  INTEGER iminn,imaxx
57 
58  REAL qmin,qmax
59  DATA qmin,qmax/0.,1.e33/
60  DATA testcpu/.false./
61  DATA temps1,temps2,temps3/0.,0.,0./
62 
63 
64  zzpbar = 0.5 * pdt
65  zzw = pdt
66  DO l=1,llm
67  DO ij = iip2,ip1jm
68  mu(ij,l)=pbaru(ij,l) * zzpbar
69  ENDDO
70  DO ij=1,ip1jm
71  mv(ij,l)=pbarv(ij,l) * zzpbar
72  ENDDO
73  DO ij=1,ip1jmp1
74  mw(ij,l)=w(ij,l) * zzw
75  ENDDO
76  ENDDO
77 
78  DO ij=1,ip1jmp1
79  mw(ij,llm+1)=0.
80  ENDDO
81 
82  CALL scopy(ijp1llm,q,1,zq,1)
83  CALL scopy(ijp1llm,masse,1,zm,1)
84 
85 cprint*,'Entree vlx1'
86 c call minmaxq(zq,qmin,qmax,'avant vlx ')
87  call vlx(zq,pente_max,zm,mu)
88 cprint*,'Sortie vlx1'
89 c call minmaxq(zq,qmin,qmax,'apres vlx1 ')
90 
91 c print*,'Entree vly1'
92  call vly(zq,pente_max,zm,mv)
93 c call minmaxq(zq,qmin,qmax,'apres vly1 ')
94 cprint*,'Sortie vly1'
95  call vlz(zq,pente_max,zm,mw)
96 c call minmaxq(zq,qmin,qmax,'apres vlz ')
97 
98 
99  call vly(zq,pente_max,zm,mv)
100 c call minmaxq(zq,qmin,qmax,'apres vly ')
101 
102 
103  call vlx(zq,pente_max,zm,mu)
104 c call minmaxq(zq,qmin,qmax,'apres vlx2 ')
105 
106 
107  DO l=1,llm
108  DO ij=1,ip1jmp1
109  q(ij,l)=zq(ij,l)
110  ENDDO
111  DO ij=1,ip1jm+1,iip1
112  q(ij+iim,l)=q(ij,l)
113  ENDDO
114  ENDDO
115 
116  RETURN
117  END
118  SUBROUTINE vlx(q,pente_max,masse,u_m)
119 
120 c Auteurs: P.Le Van, F.Hourdin, F.Forget
121 c
122 c ********************************************************************
123 c Shema d'advection " pseudo amont " .
124 c ********************************************************************
125 c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
126 c
127 c
128 c --------------------------------------------------------------------
129  IMPLICIT NONE
130 c
131  include "dimensions.h"
132  include "paramet.h"
133  include "logic.h"
134  include "comvert.h"
135  include "comconst.h"
136  include "iniprint.h"
137 c
138 c
139 c Arguments:
140 c ----------
141  REAL masse(ip1jmp1,llm),pente_max
142  REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
143  REAL q(ip1jmp1,llm)
144  REAL w(ip1jmp1,llm)
145 c
146 c Local
147 c ---------
148 c
149  INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
150  INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
151 c
152  REAL new_m,zu_m,zdum(ip1jmp1,llm)
153  REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
154  REAL zz(ip1jmp1)
155  REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
156  REAL u_mq(ip1jmp1,llm)
157 
158  Logical extremum,first,testcpu
159  SAVE first,testcpu
160 
161  REAL ssum
162  REAL temps0,temps1,temps2,temps3,temps4,temps5,second
163  SAVE temps0,temps1,temps2,temps3,temps4,temps5
164 
165  REAL z1,z2,z3
166 
167  DATA first,testcpu/.true.,.false./
168 
169  IF(first) THEN
170  temps1=0.
171  temps2=0.
172  temps3=0.
173  temps4=0.
174  temps5=0.
175  first=.false.
176  ENDIF
177 
178 c calcul de la pente a droite et a gauche de la maille
179 
180 
181  IF (pente_max.gt.-1.e-5) THEN
182 c IF (pente_max.gt.10) THEN
183 
184 c calcul des pentes avec limitation, Van Leer scheme I:
185 c -----------------------------------------------------
186 
187 c calcul de la pente aux points u
188  DO l = 1, llm
189  DO ij=iip2,ip1jm-1
190  dxqu(ij)=q(ij+1,l)-q(ij,l)
191 c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
192 c sigu(ij)=u_m(ij,l)/masse(ij,l)
193  ENDDO
194  DO ij=iip1+iip1,ip1jm,iip1
195  dxqu(ij)=dxqu(ij-iim)
196 c sigu(ij)=sigu(ij-iim)
197  ENDDO
198 
199  DO ij=iip2,ip1jm
200  adxqu(ij)=abs(dxqu(ij))
201  ENDDO
202 
203 c calcul de la pente maximum dans la maille en valeur absolue
204 
205  DO ij=iip2+1,ip1jm
206  dxqmax(ij,l)=pente_max*
207  , min(adxqu(ij-1),adxqu(ij))
208 c limitation subtile
209 c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
210 
211 
212  ENDDO
213 
214  DO ij=iip1+iip1,ip1jm,iip1
215  dxqmax(ij-iim,l)=dxqmax(ij,l)
216  ENDDO
217 
218  DO ij=iip2+1,ip1jm
219 #ifdef CRAY
220  dxq(ij,l)=
221  , cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
222 #else
223  IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
224  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
225  ELSE
226 c extremum local
227  dxq(ij,l)=0.
228  ENDIF
229 #endif
230  dxq(ij,l)=0.5*dxq(ij,l)
231  dxq(ij,l)=
232  , sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
233  ENDDO
234 
235  ENDDO ! l=1,llm
236 cprint*,'Ok calcul des pentes'
237 
238  ELSE ! (pente_max.lt.-1.e-5)
239 
240 c Pentes produits:
241 c ----------------
242 
243  DO l = 1, llm
244  DO ij=iip2,ip1jm-1
245  dxqu(ij)=q(ij+1,l)-q(ij,l)
246  ENDDO
247  DO ij=iip1+iip1,ip1jm,iip1
248  dxqu(ij)=dxqu(ij-iim)
249  ENDDO
250 
251  DO ij=iip2+1,ip1jm
252  zz(ij)=dxqu(ij-1)*dxqu(ij)
253  zz(ij)=zz(ij)+zz(ij)
254  IF(zz(ij).gt.0) THEN
255  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
256  ELSE
257 c extremum local
258  dxq(ij,l)=0.
259  ENDIF
260  ENDDO
261 
262  ENDDO
263 
264  ENDIF ! (pente_max.lt.-1.e-5)
265 
266 c bouclage de la pente en iip1:
267 c -----------------------------
268 
269  DO l=1,llm
270  DO ij=iip1+iip1,ip1jm,iip1
271  dxq(ij-iim,l)=dxq(ij,l)
272  ENDDO
273  DO ij=1,ip1jmp1
274  iadvplus(ij,l)=0
275  ENDDO
276 
277  ENDDO
278 
279 c print*,'Bouclage en iip1'
280 
281 c calcul des flux a gauche et a droite
282 
283 #ifdef CRAY
284 
285  DO l=1,llm
286  DO ij=iip2,ip1jm-1
287  zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
288  , 1.+u_m(ij,l)/masse(ij+1,l),
289  , u_m(ij,l))
290  zdum(ij,l)=0.5*zdum(ij,l)
291  u_mq(ij,l)=cvmgp(
292  , q(ij,l)+zdum(ij,l)*dxq(ij,l),
293  , q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
294  , u_m(ij,l))
295  u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
296  ENDDO
297  ENDDO
298 #else
299 c on cumule le flux correspondant a toutes les mailles dont la masse
300 c au travers de la paroi pENDant le pas de temps.
301 cprint*,'Cumule ....'
302 
303  DO l=1,llm
304  DO ij=iip2,ip1jm-1
305 c print*,'masse(',ij,')=',masse(ij,l)
306  IF (u_m(ij,l).gt.0.) THEN
307  zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
308  u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
309  ELSE
310  zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
311  u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
312  ENDIF
313  ENDDO
314  ENDDO
315 #endif
316 c stop
317 
318 c go to 9999
319 c detection des points ou on advecte plus que la masse de la
320 c maille
321  DO l=1,llm
322  DO ij=iip2,ip1jm-1
323  IF(zdum(ij,l).lt.0) THEN
324  iadvplus(ij,l)=1
325  u_mq(ij,l)=0.
326  ENDIF
327  ENDDO
328  ENDDO
329 cprint*,'Ok test 1'
330  DO l=1,llm
331  DO ij=iip1+iip1,ip1jm,iip1
332  iadvplus(ij,l)=iadvplus(ij-iim,l)
333  ENDDO
334  ENDDO
335 c print*,'Ok test 2'
336 
337 
338 c traitement special pour le cas ou on advecte en longitude plus que le
339 c contenu de la maille.
340 c cette partie est mal vectorisee.
341 
342 c calcul du nombre de maille sur lequel on advecte plus que la maille.
343 
344  n0=0
345  DO l=1,llm
346  nl(l)=0
347  DO ij=iip2,ip1jm
348  nl(l)=nl(l)+iadvplus(ij,l)
349  ENDDO
350  n0=n0+nl(l)
351  ENDDO
352 
353  IF(n0.gt.0) THEN
354  if (prt_level > 2) print *,
355  $ 'Nombre de points pour lesquels on advect plus que le'
356  & ,'contenu de la maille : ',n0
357 
358  DO l=1,llm
359  IF(nl(l).gt.0) THEN
360  iju=0
361 c indicage des mailles concernees par le traitement special
362  DO ij=iip2,ip1jm
363  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
364  iju=iju+1
365  indu(iju)=ij
366  ENDIF
367  ENDDO
368  niju=iju
369 c PRINT*,'niju,nl',niju,nl(l)
370 
371 c traitement des mailles
372  DO iju=1,niju
373  ij=indu(iju)
374  j=(ij-1)/iip1+1
375  zu_m=u_m(ij,l)
376  u_mq(ij,l)=0.
377  IF(zu_m.gt.0.) THEN
378  ijq=ij
379  i=ijq-(j-1)*iip1
380 c accumulation pour les mailles completements advectees
381  do while(zu_m.gt.masse(ijq,l))
382  u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
383  zu_m=zu_m-masse(ijq,l)
384  i=mod(i-2+iim,iim)+1
385  ijq=(j-1)*iip1+i
386  ENDDO
387 c ajout de la maille non completement advectee
388  u_mq(ij,l)=u_mq(ij,l)+zu_m*
389  & (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
390  ELSE
391  ijq=ij+1
392  i=ijq-(j-1)*iip1
393 c accumulation pour les mailles completements advectees
394  do while(-zu_m.gt.masse(ijq,l))
395  u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
396  zu_m=zu_m+masse(ijq,l)
397  i=mod(i,iim)+1
398  ijq=(j-1)*iip1+i
399  ENDDO
400 c ajout de la maille non completement advectee
401  u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
402  & 0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
403  ENDIF
404  ENDDO
405  ENDIF
406  ENDDO
407  ENDIF ! n0.gt.0
408 9999 continue
409 
410 
411 c bouclage en latitude
412 cprint*,'cvant bouclage en latitude'
413  DO l=1,llm
414  DO ij=iip1+iip1,ip1jm,iip1
415  u_mq(ij,l)=u_mq(ij-iim,l)
416  ENDDO
417  ENDDO
418 
419 
420 c calcul des tENDances
421 
422  DO l=1,llm
423  DO ij=iip2+1,ip1jm
424  new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
425  q(ij,l)=(q(ij,l)*masse(ij,l)+
426  & u_mq(ij-1,l)-u_mq(ij,l))
427  & /new_m
428  masse(ij,l)=new_m
429  ENDDO
430 c ModIF Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
431  DO ij=iip1+iip1,ip1jm,iip1
432  q(ij-iim,l)=q(ij,l)
433  masse(ij-iim,l)=masse(ij,l)
434  ENDDO
435  ENDDO
436 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
437 c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
438 
439 
440  RETURN
441  END
442  SUBROUTINE vly(q,pente_max,masse,masse_adv_v)
443 c
444 c Auteurs: P.Le Van, F.Hourdin, F.Forget
445 c
446 c ********************************************************************
447 c Shema d'advection " pseudo amont " .
448 c ********************************************************************
449 c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
450 c dq sont des arguments de sortie pour le s-pg ....
451 c
452 c
453 c --------------------------------------------------------------------
454  IMPLICIT NONE
455 c
456 #include "dimensions.h"
457 #include "paramet.h"
458 #include "logic.h"
459 #include "comvert.h"
460 #include "comconst.h"
461 #include "comgeom.h"
462 c
463 c
464 c Arguments:
465 c ----------
466  REAL masse(ip1jmp1,llm),pente_max
467  REAL masse_adv_v( ip1jm,llm)
468  REAL q(ip1jmp1,llm), dq( ip1jmp1,llm)
469 c
470 c Local
471 c ---------
472 c
473  INTEGER i,ij,l
474 c
475  REAL airej2,airejjm,airescb(iim),airesch(iim)
476  REAL dyq(ip1jmp1,llm),dyqv(ip1jm),zdvm(ip1jmp1,llm)
477  REAL adyqv(ip1jm),dyqmax(ip1jmp1)
478  REAL qbyv(ip1jm,llm)
479 
480  REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
481 c REAL newq,oldmasse
482  Logical extremum,first,testcpu
483  REAL temps0,temps1,temps2,temps3,temps4,temps5,second
484  SAVE temps0,temps1,temps2,temps3,temps4,temps5
485  SAVE first,testcpu
486 
487  REAL convpn,convps,convmpn,convmps
488  real massepn,masseps,qpn,qps
489  REAL sinlon(iip1),sinlondlon(iip1)
490  REAL coslon(iip1),coslondlon(iip1)
491  SAVE sinlon,coslon,sinlondlon,coslondlon
492  SAVE airej2,airejjm
493 c
494 c
495  REAL ssum
496 
497  DATA first,testcpu/.true.,.false./
498  DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
499 
500  IF(first) THEN
501  print*,'Shema Amont nouveau appele dans Vanleer '
502  first=.false.
503  do i=2,iip1
504  coslon(i)=cos(rlonv(i))
505  sinlon(i)=sin(rlonv(i))
506  coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
507  sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
508  ENDDO
509  coslon(1)=coslon(iip1)
510  coslondlon(1)=coslondlon(iip1)
511  sinlon(1)=sinlon(iip1)
512  sinlondlon(1)=sinlondlon(iip1)
513  airej2 = ssum( iim, aire(iip2), 1 )
514  airejjm= ssum( iim, aire(ip1jm -iim), 1 )
515  ENDIF
516 
517 c
518 cPRINT*,'CALCUL EN LATITUDE'
519 
520  DO l = 1, llm
521 c
522 c --------------------------------
523 c CALCUL EN LATITUDE
524 c --------------------------------
525 
526 c On commence par calculer la valeur du traceur moyenne sur le premier cercle
527 c de latitude autour du pole (qpns pour le pole nord et qpsn pour
528 c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
529 
530  DO i = 1, iim
531  airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
532  airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
533  ENDDO
534  qpns = ssum( iim, airescb ,1 ) / airej2
535  qpsn = ssum( iim, airesch ,1 ) / airejjm
536 
537 c calcul des pentes aux points v
538 
539  DO ij=1,ip1jm
540  dyqv(ij)=q(ij,l)-q(ij+iip1,l)
541  adyqv(ij)=abs(dyqv(ij))
542  ENDDO
543 
544 c calcul des pentes aux points scalaires
545 
546  DO ij=iip2,ip1jm
547  dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
548  dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
549  dyqmax(ij)=pente_max*dyqmax(ij)
550  ENDDO
551 
552 c calcul des pentes aux poles
553 
554  DO ij=1,iip1
555  dyq(ij,l)=qpns-q(ij+iip1,l)
556  dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
557  ENDDO
558 
559 c filtrage de la derivee
560  dyn1=0.
561  dys1=0.
562  dyn2=0.
563  dys2=0.
564  DO ij=1,iim
565  dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
566  dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
567  dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
568  dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
569  ENDDO
570  DO ij=1,iip1
571  dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
572  dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
573  ENDDO
574 
575 c calcul des pentes limites aux poles
576 
577  goto 8888
578  fn=1.
579  fs=1.
580  DO ij=1,iim
581  IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
582  fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
583  ENDIF
584  IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
585  fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
586  ENDIF
587  ENDDO
588  DO ij=1,iip1
589  dyq(ij,l)=fn*dyq(ij,l)
590  dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
591  ENDDO
592 8888 continue
593  DO ij=1,iip1
594  dyq(ij,l)=0.
595  dyq(ip1jm+ij,l)=0.
596  ENDDO
597 
598 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
599 C En memoire de dIFferents tests sur la
600 C limitation des pentes aux poles.
601 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
602 C PRINT*,dyq(1)
603 C PRINT*,dyqv(iip1+1)
604 C appn=abs(dyq(1)/dyqv(iip1+1))
605 C PRINT*,dyq(ip1jm+1)
606 C PRINT*,dyqv(ip1jm-iip1+1)
607 C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
608 C DO ij=2,iim
609 C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
610 C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
611 C ENDDO
612 C appn=min(pente_max/appn,1.)
613 C apps=min(pente_max/apps,1.)
614 C
615 C
616 C cas ou on a un extremum au pole
617 C
618 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
619 C & appn=0.
620 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
621 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
622 C & apps=0.
623 C
624 C limitation des pentes aux poles
625 C DO ij=1,iip1
626 C dyq(ij)=appn*dyq(ij)
627 C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
628 C ENDDO
629 C
630 C test
631 C DO ij=1,iip1
632 C dyq(iip1+ij)=0.
633 C dyq(ip1jm+ij-iip1)=0.
634 C ENDDO
635 C DO ij=1,ip1jmp1
636 C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
637 C ENDDO
638 C
639 C changement 10 07 96
640 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
641 C & THEN
642 C DO ij=1,iip1
643 C dyqmax(ij)=0.
644 C ENDDO
645 C ELSE
646 C DO ij=1,iip1
647 C dyqmax(ij)=pente_max*abs(dyqv(ij))
648 C ENDDO
649 C ENDIF
650 C
651 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
652 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
653 C &THEN
654 C DO ij=ip1jm+1,ip1jmp1
655 C dyqmax(ij)=0.
656 C ENDDO
657 C ELSE
658 C DO ij=ip1jm+1,ip1jmp1
659 C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
660 C ENDDO
661 C ENDIF
662 C fin changement 10 07 96
663 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
664 
665 c calcul des pentes limitees
666 
667  DO ij=iip2,ip1jm
668  IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
669  dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
670  ELSE
671  dyq(ij,l)=0.
672  ENDIF
673  ENDDO
674 
675  ENDDO
676 
677  DO l=1,llm
678  DO ij=1,ip1jm
679  IF(masse_adv_v(ij,l).gt.0) THEN
680  qbyv(ij,l)=q(ij+iip1,l)+dyq(ij+iip1,l)*
681  , 0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l))
682  ELSE
683  qbyv(ij,l)=q(ij,l)-dyq(ij,l)*
684  , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l))
685  ENDIF
686  qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
687  ENDDO
688  ENDDO
689 
690 
691  DO l=1,llm
692  DO ij=iip2,ip1jm
693  newmasse=masse(ij,l)
694  & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
695  q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
696  & /newmasse
697  masse(ij,l)=newmasse
698  ENDDO
699 c.-. ancienne version
700 c convpn=SSUM(iim,qbyv(1,l),1)/apoln
701 c convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
702 
703  convpn=ssum(iim,qbyv(1,l),1)
704  convmpn=ssum(iim,masse_adv_v(1,l),1)
705  massepn=ssum(iim,masse(1,l),1)
706  qpn=0.
707  do ij=1,iim
708  qpn=qpn+masse(ij,l)*q(ij,l)
709  enddo
710  qpn=(qpn+convpn)/(massepn+convmpn)
711  do ij=1,iip1
712  q(ij,l)=qpn
713  enddo
714 
715 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
716 c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
717 
718  convps=-ssum(iim,qbyv(ip1jm-iim,l),1)
719  convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
720  masseps=ssum(iim, masse(ip1jm+1,l),1)
721  qps=0.
722  do ij = ip1jm+1,ip1jmp1-1
723  qps=qps+masse(ij,l)*q(ij,l)
724  enddo
725  qps=(qps+convps)/(masseps+convmps)
726  do ij=ip1jm+1,ip1jmp1
727  q(ij,l)=qps
728  enddo
729 
730 c.-. fin ancienne version
731 
732 c._. nouvelle version
733 c convpn=SSUM(iim,qbyv(1,l),1)
734 c convmpn=ssum(iim,masse_adv_v(1,l),1)
735 c oldmasse=ssum(iim,masse(1,l),1)
736 c newmasse=oldmasse+convmpn
737 c newq=(q(1,l)*oldmasse+convpn)/newmasse
738 c newmasse=newmasse/apoln
739 c DO ij = 1,iip1
740 c q(ij,l)=newq
741 c masse(ij,l)=newmasse*aire(ij)
742 c ENDDO
743 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
744 c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
745 c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
746 c newmasse=oldmasse+convmps
747 c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
748 c newmasse=newmasse/apols
749 c DO ij = ip1jm+1,ip1jmp1
750 c q(ij,l)=newq
751 c masse(ij,l)=newmasse*aire(ij)
752 c ENDDO
753 c._. fin nouvelle version
754  ENDDO
755 
756  RETURN
757  END
758  SUBROUTINE vlz(q,pente_max,masse,w)
759 c
760 c Auteurs: P.Le Van, F.Hourdin, F.Forget
761 c
762 c ********************************************************************
763 c Shema d'advection " pseudo amont " .
764 c ********************************************************************
765 c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
766 c dq sont des arguments de sortie pour le s-pg ....
767 c
768 c
769 c --------------------------------------------------------------------
770  IMPLICIT NONE
771 c
772 #include "dimensions.h"
773 #include "paramet.h"
774 #include "logic.h"
775 #include "comvert.h"
776 #include "comconst.h"
777 c
778 c
779 c Arguments:
780 c ----------
781  REAL masse(ip1jmp1,llm),pente_max
782  REAL q(ip1jmp1,llm)
783  REAL w(ip1jmp1,llm+1)
784 c
785 c Local
786 c ---------
787 c
788  INTEGER i,ij,l,j,ii
789 c
790  REAL wq(ip1jmp1,llm+1),newmasse
791 
792  REAL dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm),dzqmax
793  REAL sigw
794 
795  LOGICAL testcpu
796  SAVE testcpu
797 
798  REAL temps0,temps1,temps2,temps3,temps4,temps5,second
799  SAVE temps0,temps1,temps2,temps3,temps4,temps5
800  REAL ssum
801 
802  DATA testcpu/.false./
803  DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
804 
805 c On oriente tout dans le sens de la pression c'est a dire dans le
806 c sens de W
807 
808 #ifdef BIDON
809  IF(testcpu) THEN
810  temps0=second(0.)
811  ENDIF
812 #endif
813  DO l=2,llm
814  DO ij=1,ip1jmp1
815  dzqw(ij,l)=q(ij,l-1)-q(ij,l)
816  adzqw(ij,l)=abs(dzqw(ij,l))
817  ENDDO
818  ENDDO
819 
820  DO l=2,llm-1
821  DO ij=1,ip1jmp1
822 #ifdef CRAY
823  dzq(ij,l)=0.5*
824  , cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1))
825 #else
826  IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
827  dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
828  ELSE
829  dzq(ij,l)=0.
830  ENDIF
831 #endif
832  dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
833  dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
834  ENDDO
835  ENDDO
836 
837  DO ij=1,ip1jmp1
838  dzq(ij,1)=0.
839  dzq(ij,llm)=0.
840  ENDDO
841 
842 #ifdef BIDON
843  IF(testcpu) THEN
844  temps1=temps1+second(0.)-temps0
845  ENDIF
846 #endif
847 c ---------------------------------------------------------------
848 c .... calcul des termes d'advection verticale .......
849 c ---------------------------------------------------------------
850 
851 c calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq
852 
853  DO l = 1,llm-1
854  do ij = 1,ip1jmp1
855  IF(w(ij,l+1).gt.0.) THEN
856  sigw=w(ij,l+1)/masse(ij,l+1)
857  wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
858  ELSE
859  sigw=w(ij,l+1)/masse(ij,l)
860  wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
861  ENDIF
862  ENDDO
863  ENDDO
864 
865  DO ij=1,ip1jmp1
866  wq(ij,llm+1)=0.
867  wq(ij,1)=0.
868  ENDDO
869 
870  DO l=1,llm
871  DO ij=1,ip1jmp1
872  newmasse=masse(ij,l)+w(ij,l+1)-w(ij,l)
873  q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))
874  & /newmasse
875  masse(ij,l)=newmasse
876  ENDDO
877  ENDDO
878 
879 
880  RETURN
881  END
882 c SUBROUTINE minmaxq(zq,qmin,qmax,comment)
883 c
884 c#include "dimensions.h"
885 c#include "paramet.h"
886 
887 c CHARACTER*(*) comment
888 c real qmin,qmax
889 c real zq(ip1jmp1,llm)
890 
891 c INTEGER jadrs(ip1jmp1), jbad, k, i
892 
893 
894 c DO k = 1, llm
895 c jbad = 0
896 c DO i = 1, ip1jmp1
897 c IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
898 c jbad = jbad + 1
899 c jadrs(jbad) = i
900 c ENDIF
901 c ENDDO
902 c IF (jbad.GT.0) THEN
903 c PRINT*, comment
904 c DO i = 1, jbad
905 cc PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
906 c ENDDO
907 c ENDIF
908 c ENDDO
909 
910 c return
911 c end
912  subroutine minmaxq(zq,qmin,qmax,comment)
913 
914 #include "dimensions.h"
915 #include "paramet.h"
916 
917  character*20 comment
918  real qmin,qmax
919  real zq(ip1jmp1,llm)
920  real zzq(iip1,jjp1,llm)
921 
922  integer imin,jmin,lmin,ijlmin
923  integer imax,jmax,lmax,ijlmax
924 
925  integer ismin,ismax
926 
927 #ifdef isminismax
928  call scopy(ip1jmp1*llm,zq,1,zzq,1)
929 
930  ijlmin=ismin(ijp1llm,zq,1)
931  lmin=(ijlmin-1)/ip1jmp1+1
932  ijlmin=ijlmin-(lmin-1.)*ip1jmp1
933  jmin=(ijlmin-1)/iip1+1
934  imin=ijlmin-(jmin-1.)*iip1
935  zqmin=zq(ijlmin,lmin)
936 
937  ijlmax=ismax(ijp1llm,zq,1)
938  lmax=(ijlmax-1)/ip1jmp1+1
939  ijlmax=ijlmax-(lmax-1.)*ip1jmp1
940  jmax=(ijlmax-1)/iip1+1
941  imax=ijlmax-(jmax-1.)*iip1
942  zqmax=zq(ijlmax,lmax)
943 
944  if(zqmin.lt.qmin)
945 c s write(*,9999) comment,
946  s write(*,*) comment,
947  s imin,jmin,lmin,zqmin,zzq(imin,jmin,lmin)
948  if(zqmax.gt.qmax)
949 c s write(*,9999) comment,
950  s write(*,*) comment,
951  s imax,jmax,lmax,zqmax,zzq(imax,jmax,lmax)
952 
953 #endif
954  return
955 9999 format(a20,' q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
956  end
957 
958 
959