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