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