My Project
 All Classes Files Functions Variables Macros
vlspltqs.F
Go to the documentation of this file.
1 c
2 c $Id: vlspltqs.F 1520 2011-05-23 11:37:09Z emillour $
3 c
4  SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
5  , p,pk,teta )
6 c
7 c Auteurs: P.Le Van, F.Hourdin, F.Forget, F.Codron
8 c
9 c ********************************************************************
10 c Shema d'advection " pseudo amont " .
11 c + test sur humidite specifique: Q advecte< Qsat aval
12 c (F. Codron, 10/99)
13 c ********************************************************************
14 c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
15 c
16 c pente_max facteur de limitation des pentes: 2 en general
17 c 0 pour un schema amont
18 c pbaru,pbarv,w flux de masse en u ,v ,w
19 c pdt pas de temps
20 c
21 c teta temperature potentielle, p pression aux interfaces,
22 c pk exner au milieu des couches necessaire pour calculer Qsat
23 c --------------------------------------------------------------------
24  IMPLICIT NONE
25 c
26 #include "dimensions.h"
27 #include "paramet.h"
28 #include "logic.h"
29 #include "comvert.h"
30 #include "comconst.h"
31 
32 c
33 c Arguments:
34 c ----------
35  REAL masse(ip1jmp1,llm),pente_max
36  REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
37  REAL q(ip1jmp1,llm)
38  REAL w(ip1jmp1,llm),pdt
39  REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
40 c
41 c Local
42 c ---------
43 c
44  INTEGER i,ij,l,j,ii
45 c
46  REAL qsat(ip1jmp1,llm)
47  REAL zm(ip1jmp1,llm)
48  REAL mu(ip1jmp1,llm)
49  REAL mv(ip1jm,llm)
50  REAL mw(ip1jmp1,llm+1)
51  REAL zq(ip1jmp1,llm)
52  REAL temps1,temps2,temps3
53  REAL zzpbar, zzw
54  LOGICAL testcpu
55  SAVE testcpu
56  SAVE temps1,temps2,temps3
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 c--pour rapport de melange saturant--
64 
65  REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
66  REAL ptarg,pdelarg,foeew,zdelta
67  REAL tempe(ip1jmp1)
68 
69 c fonction psat(T)
70 
71  foeew( ptarg,pdelarg ) = exp(
72  * (r3les*(1.-pdelarg)+r3ies*pdelarg) * (ptarg-rtt)
73  * / (ptarg-(r4les*(1.-pdelarg)+r4ies*pdelarg)) )
74 
75  r2es = 380.11733
76  r3les = 17.269
77  r3ies = 21.875
78  r4les = 35.86
79  r4ies = 7.66
80  retv = 0.6077667
81  rtt = 273.16
82 
83 c-- Calcul de Qsat en chaque point
84 c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
85 c pour eviter une exponentielle.
86  DO l = 1, llm
87  DO ij = 1, ip1jmp1
88  tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
89  ENDDO
90  DO ij = 1, ip1jmp1
91  zdelta = max( 0., sign(1., rtt - tempe(ij)) )
92  play = 0.5*(p(ij,l)+p(ij,l+1))
93  qsat(ij,l) = min(0.5, r2es* foeew(tempe(ij),zdelta) / play )
94  qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
95  ENDDO
96  ENDDO
97 
98 c PRINT*,'Debut vlsplt version debug sans vlyqs'
99 
100  zzpbar = 0.5 * pdt
101  zzw = pdt
102  DO l=1,llm
103  DO ij = iip2,ip1jm
104  mu(ij,l)=pbaru(ij,l) * zzpbar
105  ENDDO
106  DO ij=1,ip1jm
107  mv(ij,l)=pbarv(ij,l) * zzpbar
108  ENDDO
109  DO ij=1,ip1jmp1
110  mw(ij,l)=w(ij,l) * zzw
111  ENDDO
112  ENDDO
113 
114  DO ij=1,ip1jmp1
115  mw(ij,llm+1)=0.
116  ENDDO
117 
118  CALL scopy(ijp1llm,q,1,zq,1)
119  CALL scopy(ijp1llm,masse,1,zm,1)
120 
121 c call minmaxq(zq,qmin,qmax,'avant vlxqs ')
122  call vlxqs(zq,pente_max,zm,mu,qsat)
123 
124 
125 c call minmaxq(zq,qmin,qmax,'avant vlyqs ')
126 
127  call vlyqs(zq,pente_max,zm,mv,qsat)
128 
129 
130 c call minmaxq(zq,qmin,qmax,'avant vlz ')
131 
132  call vlz(zq,pente_max,zm,mw)
133 
134 
135 c call minmaxq(zq,qmin,qmax,'avant vlyqs ')
136 c call minmaxq(zm,qmin,qmax,'M avant vlyqs ')
137 
138  call vlyqs(zq,pente_max,zm,mv,qsat)
139 
140 
141 c call minmaxq(zq,qmin,qmax,'avant vlxqs ')
142 c call minmaxq(zm,qmin,qmax,'M avant vlxqs ')
143 
144  call vlxqs(zq,pente_max,zm,mu,qsat)
145 
146 c call minmaxq(zq,qmin,qmax,'apres vlxqs ')
147 c call minmaxq(zm,qmin,qmax,'M apres vlxqs ')
148 
149 
150  DO l=1,llm
151  DO ij=1,ip1jmp1
152  q(ij,l)=zq(ij,l)
153  ENDDO
154  DO ij=1,ip1jm+1,iip1
155  q(ij+iim,l)=q(ij,l)
156  ENDDO
157  ENDDO
158 
159  RETURN
160  END
161  SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat)
162 c
163 c Auteurs: P.Le Van, F.Hourdin, F.Forget
164 c
165 c ********************************************************************
166 c Shema d'advection " pseudo amont " .
167 c ********************************************************************
168 c
169 c --------------------------------------------------------------------
170  IMPLICIT NONE
171 c
172 #include "dimensions.h"
173 #include "paramet.h"
174 #include "logic.h"
175 #include "comvert.h"
176 #include "comconst.h"
177 c
178 c
179 c Arguments:
180 c ----------
181  REAL masse(ip1jmp1,llm),pente_max
182  REAL u_m( ip1jmp1,llm )
183  REAL q(ip1jmp1,llm)
184  REAL qsat(ip1jmp1,llm)
185 c
186 c Local
187 c ---------
188 c
189  INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
190  INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
191 c
192  REAL new_m,zu_m,zdum(ip1jmp1,llm)
193  REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
194  REAL zz(ip1jmp1)
195  REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
196  REAL u_mq(ip1jmp1,llm)
197 
198  Logical first,testcpu
199  SAVE first,testcpu
200 
201  REAL ssum
202  REAL temps0,temps1,temps2,temps3,temps4,temps5
203  SAVE temps0,temps1,temps2,temps3,temps4,temps5
204 
205 
206  DATA first,testcpu/.true.,.false./
207 
208  IF(first) THEN
209  temps1=0.
210  temps2=0.
211  temps3=0.
212  temps4=0.
213  temps5=0.
214  first=.false.
215  ENDIF
216 
217 c calcul de la pente a droite et a gauche de la maille
218 
219 
220  IF (pente_max.gt.-1.e-5) THEN
221 c IF (pente_max.gt.10) THEN
222 
223 c calcul des pentes avec limitation, Van Leer scheme I:
224 c -----------------------------------------------------
225 
226 c calcul de la pente aux points u
227  DO l = 1, llm
228  DO ij=iip2,ip1jm-1
229  dxqu(ij)=q(ij+1,l)-q(ij,l)
230 c IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
231 c sigu(ij)=u_m(ij,l)/masse(ij,l)
232  ENDDO
233  DO ij=iip1+iip1,ip1jm,iip1
234  dxqu(ij)=dxqu(ij-iim)
235 c sigu(ij)=sigu(ij-iim)
236  ENDDO
237 
238  DO ij=iip2,ip1jm
239  adxqu(ij)=abs(dxqu(ij))
240  ENDDO
241 
242 c calcul de la pente maximum dans la maille en valeur absolue
243 
244  DO ij=iip2+1,ip1jm
245  dxqmax(ij,l)=pente_max*
246  , min(adxqu(ij-1),adxqu(ij))
247 c limitation subtile
248 c , min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
249 
250 
251  ENDDO
252 
253  DO ij=iip1+iip1,ip1jm,iip1
254  dxqmax(ij-iim,l)=dxqmax(ij,l)
255  ENDDO
256 
257  DO ij=iip2+1,ip1jm
258 #ifdef CRAY
259  dxq(ij,l)=
260  , cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
261 #else
262  IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
263  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
264  ELSE
265 c extremum local
266  dxq(ij,l)=0.
267  ENDIF
268 #endif
269  dxq(ij,l)=0.5*dxq(ij,l)
270  dxq(ij,l)=
271  , sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
272  ENDDO
273 
274  ENDDO ! l=1,llm
275 
276  ELSE ! (pente_max.lt.-1.e-5)
277 
278 c Pentes produits:
279 c ----------------
280 
281  DO l = 1, llm
282  DO ij=iip2,ip1jm-1
283  dxqu(ij)=q(ij+1,l)-q(ij,l)
284  ENDDO
285  DO ij=iip1+iip1,ip1jm,iip1
286  dxqu(ij)=dxqu(ij-iim)
287  ENDDO
288 
289  DO ij=iip2+1,ip1jm
290  zz(ij)=dxqu(ij-1)*dxqu(ij)
291  zz(ij)=zz(ij)+zz(ij)
292  IF(zz(ij).gt.0) THEN
293  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
294  ELSE
295 c extremum local
296  dxq(ij,l)=0.
297  ENDIF
298  ENDDO
299 
300  ENDDO
301 
302  ENDIF ! (pente_max.lt.-1.e-5)
303 
304 c bouclage de la pente en iip1:
305 c -----------------------------
306 
307  DO l=1,llm
308  DO ij=iip1+iip1,ip1jm,iip1
309  dxq(ij-iim,l)=dxq(ij,l)
310  ENDDO
311 
312  DO ij=1,ip1jmp1
313  iadvplus(ij,l)=0
314  ENDDO
315 
316  ENDDO
317 
318 
319 c calcul des flux a gauche et a droite
320 
321 #ifdef CRAY
322 c--pas encore modification sur Qsat
323  DO l=1,llm
324  DO ij=iip2,ip1jm-1
325  zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
326  , 1.+u_m(ij,l)/masse(ij+1,l),
327  , u_m(ij,l))
328  zdum(ij,l)=0.5*zdum(ij,l)
329  u_mq(ij,l)=cvmgp(
330  , q(ij,l)+zdum(ij,l)*dxq(ij,l),
331  , q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
332  , u_m(ij,l))
333  u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
334  ENDDO
335  ENDDO
336 #else
337 c on cumule le flux correspondant a toutes les mailles dont la masse
338 c au travers de la paroi pENDant le pas de temps.
339 c le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
340  DO l=1,llm
341  DO ij=iip2,ip1jm-1
342  IF (u_m(ij,l).gt.0.) THEN
343  zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
344  u_mq(ij,l)=u_m(ij,l)*
345  $ min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
346  ELSE
347  zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
348  u_mq(ij,l)=u_m(ij,l)*
349  $ min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
350  ENDIF
351  ENDDO
352  ENDDO
353 #endif
354 
355 
356 c detection des points ou on advecte plus que la masse de la
357 c maille
358  DO l=1,llm
359  DO ij=iip2,ip1jm-1
360  IF(zdum(ij,l).lt.0) THEN
361  iadvplus(ij,l)=1
362  u_mq(ij,l)=0.
363  ENDIF
364  ENDDO
365  ENDDO
366  DO l=1,llm
367  DO ij=iip1+iip1,ip1jm,iip1
368  iadvplus(ij,l)=iadvplus(ij-iim,l)
369  ENDDO
370  ENDDO
371 
372 
373 
374 c traitement special pour le cas ou on advecte en longitude plus que le
375 c contenu de la maille.
376 c cette partie est mal vectorisee.
377 
378 c pas d'influence de la pression saturante (pour l'instant)
379 
380 c calcul du nombre de maille sur lequel on advecte plus que la maille.
381 
382  n0=0
383  DO l=1,llm
384  nl(l)=0
385  DO ij=iip2,ip1jm
386  nl(l)=nl(l)+iadvplus(ij,l)
387  ENDDO
388  n0=n0+nl(l)
389  ENDDO
390 
391  IF(n0.gt.0) THEN
392 ccc PRINT*,'Nombre de points pour lesquels on advect plus que le'
393 ccc & ,'contenu de la maille : ',n0
394 
395  DO l=1,llm
396  IF(nl(l).gt.0) THEN
397  iju=0
398 c indicage des mailles concernees par le traitement special
399  DO ij=iip2,ip1jm
400  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
401  iju=iju+1
402  indu(iju)=ij
403  ENDIF
404  ENDDO
405  niju=iju
406 c PRINT*,'niju,nl',niju,nl(l)
407 
408 c traitement des mailles
409  DO iju=1,niju
410  ij=indu(iju)
411  j=(ij-1)/iip1+1
412  zu_m=u_m(ij,l)
413  u_mq(ij,l)=0.
414  IF(zu_m.gt.0.) THEN
415  ijq=ij
416  i=ijq-(j-1)*iip1
417 c accumulation pour les mailles completements advectees
418  do while(zu_m.gt.masse(ijq,l))
419  u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
420  zu_m=zu_m-masse(ijq,l)
421  i=mod(i-2+iim,iim)+1
422  ijq=(j-1)*iip1+i
423  ENDDO
424 c ajout de la maille non completement advectee
425  u_mq(ij,l)=u_mq(ij,l)+zu_m*
426  & (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
427  ELSE
428  ijq=ij+1
429  i=ijq-(j-1)*iip1
430 c accumulation pour les mailles completements advectees
431  do while(-zu_m.gt.masse(ijq,l))
432  u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
433  zu_m=zu_m+masse(ijq,l)
434  i=mod(i,iim)+1
435  ijq=(j-1)*iip1+i
436  ENDDO
437 c ajout de la maille non completement advectee
438  u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
439  & 0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
440  ENDIF
441  ENDDO
442  ENDIF
443  ENDDO
444  ENDIF ! n0.gt.0
445 
446 
447 
448 c bouclage en latitude
449 
450  DO l=1,llm
451  DO ij=iip1+iip1,ip1jm,iip1
452  u_mq(ij,l)=u_mq(ij-iim,l)
453  ENDDO
454  ENDDO
455 
456 
457 c calcul des tendances
458 
459  DO l=1,llm
460  DO ij=iip2+1,ip1jm
461  new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
462  q(ij,l)=(q(ij,l)*masse(ij,l)+
463  & u_mq(ij-1,l)-u_mq(ij,l))
464  & /new_m
465  masse(ij,l)=new_m
466  ENDDO
467 c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
468  DO ij=iip1+iip1,ip1jm,iip1
469  q(ij-iim,l)=q(ij,l)
470  masse(ij-iim,l)=masse(ij,l)
471  ENDDO
472  ENDDO
473 
474 c CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
475 c CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
476 
477 
478  RETURN
479  END
480  SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat)
481 c
482 c Auteurs: P.Le Van, F.Hourdin, F.Forget
483 c
484 c ********************************************************************
485 c Shema d'advection " pseudo amont " .
486 c ********************************************************************
487 c q,masse_adv_v,w sont des arguments d'entree pour le s-pg ....
488 c qsat est un argument de sortie pour le s-pg ....
489 c
490 c
491 c --------------------------------------------------------------------
492  IMPLICIT NONE
493 c
494 #include "dimensions.h"
495 #include "paramet.h"
496 #include "logic.h"
497 #include "comvert.h"
498 #include "comconst.h"
499 #include "comgeom.h"
500 c
501 c
502 c Arguments:
503 c ----------
504  REAL masse(ip1jmp1,llm),pente_max
505  REAL masse_adv_v( ip1jm,llm)
506  REAL q(ip1jmp1,llm)
507  REAL qsat(ip1jmp1,llm)
508 c
509 c Local
510 c ---------
511 c
512  INTEGER i,ij,l
513 c
514  REAL airej2,airejjm,airescb(iim),airesch(iim)
515  REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
516  REAL adyqv(ip1jm),dyqmax(ip1jmp1)
517  REAL qbyv(ip1jm,llm)
518 
519  REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
520 c REAL newq,oldmasse
521  Logical first,testcpu
522  REAL temps0,temps1,temps2,temps3,temps4,temps5
523  SAVE temps0,temps1,temps2,temps3,temps4,temps5
524  SAVE first,testcpu
525 
526  REAL convpn,convps,convmpn,convmps
527  REAL sinlon(iip1),sinlondlon(iip1)
528  REAL coslon(iip1),coslondlon(iip1)
529  SAVE sinlon,coslon,sinlondlon,coslondlon
530  SAVE airej2,airejjm
531 c
532 c
533  REAL ssum
534 
535  DATA first,testcpu/.true.,.false./
536  DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
537 
538  IF(first) THEN
539  print*,'Shema Amont nouveau appele dans Vanleer '
540  first=.false.
541  do i=2,iip1
542  coslon(i)=cos(rlonv(i))
543  sinlon(i)=sin(rlonv(i))
544  coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
545  sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
546  ENDDO
547  coslon(1)=coslon(iip1)
548  coslondlon(1)=coslondlon(iip1)
549  sinlon(1)=sinlon(iip1)
550  sinlondlon(1)=sinlondlon(iip1)
551  airej2 = ssum( iim, aire(iip2), 1 )
552  airejjm= ssum( iim, aire(ip1jm -iim), 1 )
553  ENDIF
554 
555 c
556 
557 
558  DO l = 1, llm
559 c
560 c --------------------------------
561 c CALCUL EN LATITUDE
562 c --------------------------------
563 
564 c On commence par calculer la valeur du traceur moyenne sur le premier cercle
565 c de latitude autour du pole (qpns pour le pole nord et qpsn pour
566 c le pole nord) qui sera utilisee pour evaluer les pentes au pole.
567 
568  DO i = 1, iim
569  airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
570  airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
571  ENDDO
572  qpns = ssum( iim, airescb ,1 ) / airej2
573  qpsn = ssum( iim, airesch ,1 ) / airejjm
574 
575 c calcul des pentes aux points v
576 
577  DO ij=1,ip1jm
578  dyqv(ij)=q(ij,l)-q(ij+iip1,l)
579  adyqv(ij)=abs(dyqv(ij))
580  ENDDO
581 
582 c calcul des pentes aux points scalaires
583 
584  DO ij=iip2,ip1jm
585  dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
586  dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
587  dyqmax(ij)=pente_max*dyqmax(ij)
588  ENDDO
589 
590 c calcul des pentes aux poles
591 
592  DO ij=1,iip1
593  dyq(ij,l)=qpns-q(ij+iip1,l)
594  dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
595  ENDDO
596 
597 c filtrage de la derivee
598  dyn1=0.
599  dys1=0.
600  dyn2=0.
601  dys2=0.
602  DO ij=1,iim
603  dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
604  dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
605  dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
606  dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
607  ENDDO
608  DO ij=1,iip1
609  dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
610  dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
611  ENDDO
612 
613 c calcul des pentes limites aux poles
614 
615  fn=1.
616  fs=1.
617  DO ij=1,iim
618  IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
619  fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
620  ENDIF
621  IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
622  fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
623  ENDIF
624  ENDDO
625  DO ij=1,iip1
626  dyq(ij,l)=fn*dyq(ij,l)
627  dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
628  ENDDO
629 
630 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
631 C En memoire de dIFferents tests sur la
632 C limitation des pentes aux poles.
633 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
634 C PRINT*,dyq(1)
635 C PRINT*,dyqv(iip1+1)
636 C appn=abs(dyq(1)/dyqv(iip1+1))
637 C PRINT*,dyq(ip1jm+1)
638 C PRINT*,dyqv(ip1jm-iip1+1)
639 C apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
640 C DO ij=2,iim
641 C appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
642 C apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
643 C ENDDO
644 C appn=min(pente_max/appn,1.)
645 C apps=min(pente_max/apps,1.)
646 C
647 C
648 C cas ou on a un extremum au pole
649 C
650 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
651 C & appn=0.
652 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
653 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
654 C & apps=0.
655 C
656 C limitation des pentes aux poles
657 C DO ij=1,iip1
658 C dyq(ij)=appn*dyq(ij)
659 C dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
660 C ENDDO
661 C
662 C test
663 C DO ij=1,iip1
664 C dyq(iip1+ij)=0.
665 C dyq(ip1jm+ij-iip1)=0.
666 C ENDDO
667 C DO ij=1,ip1jmp1
668 C dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
669 C ENDDO
670 C
671 C changement 10 07 96
672 C IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
673 C & THEN
674 C DO ij=1,iip1
675 C dyqmax(ij)=0.
676 C ENDDO
677 C ELSE
678 C DO ij=1,iip1
679 C dyqmax(ij)=pente_max*abs(dyqv(ij))
680 C ENDDO
681 C ENDIF
682 C
683 C IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
684 C & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
685 C &THEN
686 C DO ij=ip1jm+1,ip1jmp1
687 C dyqmax(ij)=0.
688 C ENDDO
689 C ELSE
690 C DO ij=ip1jm+1,ip1jmp1
691 C dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
692 C ENDDO
693 C ENDIF
694 C fin changement 10 07 96
695 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
696 
697 c calcul des pentes limitees
698 
699  DO ij=iip2,ip1jm
700  IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
701  dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
702  ELSE
703  dyq(ij,l)=0.
704  ENDIF
705  ENDDO
706 
707  ENDDO
708 
709  DO l=1,llm
710  DO ij=1,ip1jm
711  IF( masse_adv_v(ij,l).GT.0. ) THEN
712  qbyv(ij,l)= min( qsat(ij+iip1,l), q(ij+iip1,l ) +
713  , dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
714  ELSE
715  qbyv(ij,l)= min( qsat(ij,l), q(ij,l) - dyq(ij,l) *
716  , 0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
717  ENDIF
718  qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
719  ENDDO
720  ENDDO
721 
722 
723  DO l=1,llm
724  DO ij=iip2,ip1jm
725  newmasse=masse(ij,l)
726  & +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
727  q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
728  & /newmasse
729  masse(ij,l)=newmasse
730  ENDDO
731 c.-. ancienne version
732  convpn=ssum(iim,qbyv(1,l),1)/apoln
733  convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
734  DO ij = 1,iip1
735  newmasse=masse(ij,l)+convmpn*aire(ij)
736  q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
737  & newmasse
738  masse(ij,l)=newmasse
739  ENDDO
740  convps = -ssum(iim,qbyv(ip1jm-iim,l),1)/apols
741  convmps = -ssum(iim,masse_adv_v(ip1jm-iim,l),1)/apols
742  DO ij = ip1jm+1,ip1jmp1
743  newmasse=masse(ij,l)+convmps*aire(ij)
744  q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
745  & newmasse
746  masse(ij,l)=newmasse
747  ENDDO
748 c.-. fin ancienne version
749 
750 c._. nouvelle version
751 c convpn=SSUM(iim,qbyv(1,l),1)
752 c convmpn=ssum(iim,masse_adv_v(1,l),1)
753 c oldmasse=ssum(iim,masse(1,l),1)
754 c newmasse=oldmasse+convmpn
755 c newq=(q(1,l)*oldmasse+convpn)/newmasse
756 c newmasse=newmasse/apoln
757 c DO ij = 1,iip1
758 c q(ij,l)=newq
759 c masse(ij,l)=newmasse*aire(ij)
760 c ENDDO
761 c convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
762 c convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
763 c oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
764 c newmasse=oldmasse+convmps
765 c newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
766 c newmasse=newmasse/apols
767 c DO ij = ip1jm+1,ip1jmp1
768 c q(ij,l)=newq
769 c masse(ij,l)=newmasse*aire(ij)
770 c ENDDO
771 c._. fin nouvelle version
772  ENDDO
773 
774  RETURN
775  END