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