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