My Project
 All Classes Files Functions Variables Macros
advn.F
Go to the documentation of this file.
1 !
2 ! $Header$
3 !
4  SUBROUTINE advn(q,masse,w,pbaru,pbarv,pdt,mode)
5 c
6 c Auteur : F. Hourdin
7 c
8 c ********************************************************************
9 c Shema d'advection " pseudo amont " .
10 c ********************************************************************
11 c q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
12 c
13 c pbaru,pbarv,w flux de masse en u ,v ,w
14 c pdt pas de temps
15 c
16 c --------------------------------------------------------------------
17  IMPLICIT NONE
18 c
19 #include "dimensions.h"
20 #include "paramet.h"
21 #include "logic.h"
22 #include "comvert.h"
23 #include "comconst.h"
24 #include "comgeom.h"
25 #include "iniprint.h"
26 
27 c
28 c Arguments:
29 c ----------
30  integer mode
31  real masse(ip1jmp1,llm)
32  REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
33  REAL q(ip1jmp1,llm)
34  REAL w(ip1jmp1,llm),pdt
35 c
36 c Local
37 c ---------
38 c
39  INTEGER i,ij,l,j,ii
40  integer ijlqmin,iqmin,jqmin,lqmin
41  integer ismin
42 c
43  real zm(ip1jmp1,llm),newmasse
44  real mu(ip1jmp1,llm)
45  real mv(ip1jm,llm)
46  real mw(ip1jmp1,llm+1)
47  real zq(ip1jmp1,llm),zz,qpn,qps
48  real zqg(ip1jmp1,llm),zqd(ip1jmp1,llm)
49  real zqs(ip1jmp1,llm),zqn(ip1jmp1,llm)
50  real zqh(ip1jmp1,llm),zqb(ip1jmp1,llm)
51  real temps0,temps1,temps2,temps3
52  real ztemps1,ztemps2,ztemps3,ssum
53  logical testcpu
54  save testcpu
55  save temps1,temps2,temps3
56  real zzpbar,zzw
57 
58 #ifdef CRAY
59  real second
60 #endif
61 
62  real qmin,qmax
63  data qmin,qmax/0.,1./
64  data testcpu/.false./
65  data temps1,temps2,temps3/0.,0.,0./
66 
67  zzpbar = 0.5 * pdt
68  zzw = pdt
69 
70  DO l=1,llm
71  DO ij = iip2,ip1jm
72  mu(ij,l)=pbaru(ij,l) * zzpbar
73  ENDDO
74  DO ij=1,ip1jm
75  mv(ij,l)=pbarv(ij,l) * zzpbar
76  ENDDO
77  DO ij=1,ip1jmp1
78  mw(ij,l)=w(ij,l) * zzw
79  ENDDO
80  ENDDO
81 
82  DO ij=1,ip1jmp1
83  mw(ij,llm+1)=0.
84  ENDDO
85 
86  do l=1,llm
87  qpn=0.
88  qps=0.
89  do ij=1,iim
90  qpn=qpn+q(ij,l)*masse(ij,l)
91  qps=qps+q(ip1jm+ij,l)*masse(ip1jm+ij,l)
92  enddo
93  qpn=qpn/ssum(iim,masse(1,l),1)
94  qps=qps/ssum(iim,masse(ip1jm+1,l),1)
95  do ij=1,iip1
96  q(ij,l)=qpn
97  q(ip1jm+ij,l)=qps
98  enddo
99  enddo
100 
101  do ij=1,ip1jmp1
102  mw(ij,llm+1)=0.
103  enddo
104  do l=1,llm
105  do ij=1,ip1jmp1
106  zq(ij,l)=q(ij,l)
107  zm(ij,l)=masse(ij,l)
108  enddo
109  enddo
110 
111 c call minmaxq(zq,qmin,qmax,'avant vlx ')
112  call advnqx(zq,zqg,zqd)
113  call advnx(zq,zqg,zqd,zm,mu,mode)
114  call advnqy(zq,zqs,zqn)
115  call advny(zq,zqs,zqn,zm,mv)
116  call advnqz(zq,zqh,zqb)
117  call advnz(zq,zqh,zqb,zm,mw)
118 c call vlz(zq,0.,zm,mw)
119  call advnqy(zq,zqs,zqn)
120  call advny(zq,zqs,zqn,zm,mv)
121  call advnqx(zq,zqg,zqd)
122  call advnx(zq,zqg,zqd,zm,mu,mode)
123 c call minmaxq(zq,qmin,qmax,'apres vlx ')
124 
125 #ifdef CRAY
126  if(testcpu) then
127  ztemps1=second(0.)
128  temps1=temps1+ztemps1-ztemps2
129  print*,'VLSPLT X:',temps1,' Y:',temps2,' Z:',temps3
130  endif
131 #endif
132  do l=1,llm
133  do ij=1,ip1jmp1
134  q(ij,l)=zq(ij,l)
135  enddo
136  do ij=1,ip1jm+1,iip1
137  q(ij+iim,l)=q(ij,l)
138  enddo
139  enddo
140 
141  RETURN
142  END
143 
144  SUBROUTINE advnqx(q,qg,qd)
145 c
146 c Auteurs: Calcul des valeurs de q aux point u.
147 c
148 c --------------------------------------------------------------------
149  IMPLICIT NONE
150 c
151 #include "dimensions.h"
152 #include "paramet.h"
153 #include "iniprint.h"
154 c
155 c
156 c Arguments:
157 c ----------
158  real q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)
159 c
160 c Local
161 c ---------
162 c
163  INTEGER ij,l
164 c
165  real dxqu(ip1jmp1),zqu(ip1jmp1)
166  real zqmax(ip1jmp1),zqmin(ip1jmp1)
167  logical extremum(ip1jmp1)
168 
169  integer mode
170  save mode
171  data mode/1/
172 
173 c calcul des pentes en u:
174 c -----------------------
175  if (mode.eq.0) then
176  do l=1,llm
177  do ij=1,ip1jm
178  qd(ij,l)=q(ij,l)
179  qg(ij,l)=q(ij,l)
180  enddo
181  enddo
182  else
183  do l = 1, llm
184  do ij=iip2,ip1jm-1
185  dxqu(ij)=q(ij+1,l)-q(ij,l)
186  zqu(ij)=0.5*(q(ij+1,l)+q(ij,l))
187  enddo
188  do ij=iip1+iip1,ip1jm,iip1
189  dxqu(ij)=dxqu(ij-iim)
190  zqu(ij)=zqu(ij-iim)
191  enddo
192  do ij=iip2,ip1jm-1
193  zqu(ij)=zqu(ij)-dxqu(ij+1)/12.
194  enddo
195  do ij=iip1+iip1,ip1jm,iip1
196  zqu(ij)=zqu(ij-iim)
197  enddo
198  do ij=iip2+1,ip1jm
199  zqu(ij)=zqu(ij)+dxqu(ij-1)/12.
200  enddo
201  do ij=iip1+iip1,ip1jm,iip1
202  zqu(ij-iim)=zqu(ij)
203  enddo
204 
205 c calcul des valeurs max et min acceptees aux interfaces
206 
207  do ij=iip2,ip1jm-1
208  zqmax(ij)=max(q(ij+1,l),q(ij,l))
209  zqmin(ij)=min(q(ij+1,l),q(ij,l))
210  enddo
211  do ij=iip1+iip1,ip1jm,iip1
212  zqmax(ij)=zqmax(ij-iim)
213  zqmin(ij)=zqmin(ij-iim)
214  enddo
215  do ij=iip2+1,ip1jm
216  extremum(ij)=dxqu(ij)*dxqu(ij-1).le.0.
217  enddo
218  do ij=iip1+iip1,ip1jm,iip1
219  extremum(ij-iim)=extremum(ij)
220  enddo
221  do ij=iip2,ip1jm
222  zqu(ij)=min(max(zqmin(ij),zqu(ij)),zqmax(ij))
223  enddo
224  do ij=iip2+1,ip1jm
225  if(extremum(ij)) then
226  qg(ij,l)=q(ij,l)
227  qd(ij,l)=q(ij,l)
228  else
229  qd(ij,l)=zqu(ij)
230  qg(ij,l)=zqu(ij-1)
231  endif
232  enddo
233  do ij=iip1+iip1,ip1jm,iip1
234  qd(ij-iim,l)=qd(ij,l)
235  qg(ij-iim,l)=qg(ij,l)
236  enddo
237 
238  goto 8888
239 
240  do ij=iip2+1,ip1jm
241  if(extremum(ij).and..not.extremum(ij-1))
242  s qd(ij-1,l)=q(ij,l)
243  enddo
244 
245  do ij=iip1+iip1,ip1jm,iip1
246  qd(ij-iim,l)=qd(ij,l)
247  enddo
248  do ij=iip2,ip1jm-1
249  if (extremum(ij).and..not.extremum(ij+1))
250  s qg(ij+1,l)=q(ij,l)
251  enddo
252 
253  do ij=iip1+iip1,ip1jm,iip1
254  qg(ij,l)=qg(ij-iim,l)
255  enddo
256 8888 continue
257  enddo
258  endif
259  RETURN
260  END
261  SUBROUTINE advnqy(q,qs,qn)
262 c
263 c Auteurs: Calcul des valeurs de q aux point v.
264 c
265 c --------------------------------------------------------------------
266  IMPLICIT NONE
267 c
268 #include "dimensions.h"
269 #include "paramet.h"
270 #include "iniprint.h"
271 c
272 c
273 c Arguments:
274 c ----------
275  real q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)
276 c
277 c Local
278 c ---------
279 c
280  INTEGER ij,l
281 c
282  real dyqv(ip1jm),zqv(ip1jm,llm)
283  real zqmax(ip1jm),zqmin(ip1jm)
284  logical extremum(ip1jmp1)
285 
286  integer mode
287  save mode
288  data mode/1/
289 
290  if (mode.eq.0) then
291  do l=1,llm
292  do ij=1,ip1jmp1
293  qn(ij,l)=q(ij,l)
294  qs(ij,l)=q(ij,l)
295  enddo
296  enddo
297  else
298 
299 c calcul des pentes en u:
300 c -----------------------
301  do l = 1, llm
302  do ij=1,ip1jm
303  dyqv(ij)=q(ij,l)-q(ij+iip1,l)
304  enddo
305 
306  do ij=iip2,ip1jm-iip1
307  zqv(ij,l)=0.5*(q(ij+iip1,l)+q(ij,l))
308  zqv(ij,l)=zqv(ij,l)+(dyqv(ij+iip1)-dyqv(ij-iip1))/12.
309  enddo
310 
311  do ij=iip2,ip1jm
312  extremum(ij)=dyqv(ij)*dyqv(ij-iip1).le.0.
313  enddo
314 
315 c Pas de pentes aux poles
316  do ij=1,iip1
317  zqv(ij,l)=q(ij,l)
318  zqv(ip1jm-iip1+ij,l)=q(ip1jm+ij,l)
319  extremum(ij)=.true.
320  extremum(ip1jmp1-iip1+ij)=.true.
321  enddo
322 
323 c calcul des valeurs max et min acceptees aux interfaces
324  do ij=1,ip1jm
325  zqmax(ij)=max(q(ij+iip1,l),q(ij,l))
326  zqmin(ij)=min(q(ij+iip1,l),q(ij,l))
327  enddo
328 
329  do ij=1,ip1jm
330  zqv(ij,l)=min(max(zqmin(ij),zqv(ij,l)),zqmax(ij))
331  enddo
332 
333  do ij=iip2,ip1jm
334  if(extremum(ij)) then
335  qs(ij,l)=q(ij,l)
336  qn(ij,l)=q(ij,l)
337 c if (.not.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
338 c if (.not.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l)
339  else
340  qs(ij,l)=zqv(ij,l)
341  qn(ij,l)=zqv(ij-iip1,l)
342  endif
343  enddo
344 
345  do ij=1,iip1
346  qs(ij,l)=q(ij,l)
347  qn(ij,l)=q(ij,l)
348  qs(ip1jm+ij,l)=q(ip1jm+ij,l)
349  qn(ip1jm+ij,l)=q(ip1jm+ij,l)
350  enddo
351 
352  enddo
353  endif
354  RETURN
355  END
356 
357  SUBROUTINE advnqz(q,qh,qb)
358 c
359 c Auteurs: Calcul des valeurs de q aux point v.
360 c
361 c --------------------------------------------------------------------
362  IMPLICIT NONE
363 c
364 #include "dimensions.h"
365 #include "paramet.h"
366 #include "iniprint.h"
367 c
368 c
369 c Arguments:
370 c ----------
371  real q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)
372 c
373 c Local
374 c ---------
375 c
376  INTEGER ij,l
377 c
378  real dzqw(ip1jmp1,llm+1),zqw(ip1jmp1,llm+1)
379  real zqmax(ip1jmp1,llm),zqmin(ip1jmp1,llm)
380  logical extremum(ip1jmp1,llm)
381 
382  integer mode
383  save mode
384 
385  data mode/1/
386 
387 c calcul des pentes en u:
388 c -----------------------
389 
390  if (mode.eq.0) then
391  do l=1,llm
392  do ij=1,ip1jmp1
393  qb(ij,l)=q(ij,l)
394  qh(ij,l)=q(ij,l)
395  enddo
396  enddo
397  else
398  do l = 2, llm
399  do ij=1,ip1jmp1
400  dzqw(ij,l)=q(ij,l-1)-q(ij,l)
401  zqw(ij,l)=0.5*(q(ij,l-1)+q(ij,l))
402  enddo
403  enddo
404  do ij=1,ip1jmp1
405  dzqw(ij,1)=0.
406  dzqw(ij,llm+1)=0.
407  enddo
408  do l=2,llm
409  do ij=1,ip1jmp1
410  zqw(ij,l)=zqw(ij,l)+(dzqw(ij,l+1)-dzqw(ij,l-1))/12.
411  enddo
412  enddo
413  do l=2,llm-1
414  do ij=1,ip1jmp1
415  extremum(ij,l)=dzqw(ij,l)*dzqw(ij,l+1).le.0.
416  enddo
417  enddo
418 
419 c Pas de pentes en bas et en haut
420  do ij=1,ip1jmp1
421  zqw(ij,2)=q(ij,1)
422  zqw(ij,llm)=q(ij,llm)
423  extremum(ij,1)=.true.
424  extremum(ij,llm)=.true.
425  enddo
426 
427 c calcul des valeurs max et min acceptees aux interfaces
428  do l=2,llm
429  do ij=1,ip1jmp1
430  zqmax(ij,l)=max(q(ij,l-1),q(ij,l))
431  zqmin(ij,l)=min(q(ij,l-1),q(ij,l))
432  enddo
433  enddo
434 
435  do l=2,llm
436  do ij=1,ip1jmp1
437  zqw(ij,l)=min(max(zqmin(ij,l),zqw(ij,l)),zqmax(ij,l))
438  enddo
439  enddo
440 
441  do l=2,llm-1
442  do ij=1,ip1jmp1
443  if(extremum(ij,l)) then
444  qh(ij,l)=q(ij,l)
445  qb(ij,l)=q(ij,l)
446  else
447  qh(ij,l)=zqw(ij,l+1)
448  qb(ij,l)=zqw(ij,l)
449  endif
450  enddo
451  enddo
452 c do l=2,llm-1
453 c do ij=1,ip1jmp1
454 c if(extremum(ij,l)) then
455 c if (.not.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
456 c if (.not.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
457 c endif
458 c enddo
459 c enddo
460 
461  do ij=1,ip1jmp1
462  qb(ij,1)=q(ij,1)
463  qh(ij,1)=q(ij,1)
464  qb(ij,llm)=q(ij,llm)
465  qh(ij,llm)=q(ij,llm)
466  enddo
467 
468  endif
469 
470  RETURN
471  END
472 
473  SUBROUTINE advnx(q,qg,qd,masse,u_m,mode)
474 c
475 c Auteur : F. Hourdin
476 c
477 c ********************************************************************
478 c Shema d'advection " pseudo amont " .
479 c ********************************************************************
480 c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
481 c
482 c
483 c --------------------------------------------------------------------
484  IMPLICIT NONE
485 c
486 #include "dimensions.h"
487 #include "paramet.h"
488 #include "logic.h"
489 #include "comvert.h"
490 #include "comconst.h"
491 #include "iniprint.h"
492 c
493 c
494 c Arguments:
495 c ----------
496  integer mode
497  real masse(ip1jmp1,llm)
498  real u_m( ip1jmp1,llm )
499  real q(ip1jmp1,llm),qd(ip1jmp1,llm),qg(ip1jmp1,llm)
500 c
501 c Local
502 c ---------
503 c
504  INTEGER i,j,ij,l,indu(ip1jmp1),niju,iju,ijq
505  integer n0,nl(llm)
506 c
507  real new_m,zu_m,zdq,zz
508  real zsigg(ip1jmp1,llm),zsigd(ip1jmp1,llm),zsig
509  real u_mq(ip1jmp1,llm)
510 
511  real zm,zq,zsigm,zsigp,zqm,zqp,zu
512 
513  logical ladvplus(ip1jmp1,llm)
514 
515  real prec
516  save prec
517 
518 #ifdef CRAY
519  data prec/1.e-24/
520 #else
521  data prec/1.e-15/
522 #endif
523 
524  do l=1,llm
525  do ij=iip2,ip1jm
526  zdq=qd(ij,l)-qg(ij,l)
527 c if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) then
528 c print*,'probleme au point ij=',ij,' l=',l
529 c print*,qd(ij,l),q(ij,l),qg(ij,l)
530 c qd(ij,l)=q(ij,l)
531 c qg(ij,l)=q(ij,l)
532 c endif
533  if(abs(zdq).gt.prec) then
534  zsigd(ij,l)=(q(ij,l)-qg(ij,l))/zdq
535  zsigg(ij,l)=1.-zsigd(ij,l)
536 c if(.not.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .and.
537 c s zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) then
538 c print*,'probleme au point ij=',ij,' l=',l
539 c print*,'sigg=',zsigg(ij,l),' sigd=',zsigd(ij,l)
540 c print*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq
541 c stop
542 c endif
543  else
544  zsigd(ij,l)=0.5
545  zsigg(ij,l)=0.5
546  qd(ij,l)=q(ij,l)
547  qg(ij,l)=q(ij,l)
548  endif
549  enddo
550  enddo
551 
552 c calcul de la pente maximum dans la maille en valeur absolue
553 
554  do l=1,llm
555  do ij=iip2,ip1jm-1
556  if (u_m(ij,l).ge.0.) then
557  zsigp=zsigd(ij,l)
558  zsigm=zsigg(ij,l)
559  zqp=qd(ij,l)
560  zqm=qg(ij,l)
561  zm=masse(ij,l)
562  zq=q(ij,l)
563  else
564  zsigm=zsigd(ij+1,l)
565  zsigp=zsigg(ij+1,l)
566  zqm=qd(ij+1,l)
567  zqp=qg(ij+1,l)
568  zm=masse(ij+1,l)
569  zq=q(ij+1,l)
570  endif
571  zu=abs(u_m(ij,l))
572  ladvplus(ij,l)=zu.gt.zm
573  zsig=zu/zm
574  if(zsig.eq.0.) zsigp=0.1
575  if (mode.eq.1) then
576  if (zsig.le.zsigp) then
577  u_mq(ij,l)=u_m(ij,l)*zqp
578  else if (mode.eq.1) then
579  u_mq(ij,l)=
580  s sign(zm,u_m(ij,l))*(zsigp*zqp+(zsig-zsigp)*zqm)
581  endif
582  else
583  if (zsig.le.zsigp) then
584  u_mq(ij,l)=u_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
585  else
586  zz=0.5*(zsig-zsigp)/zsigm
587  u_mq(ij,l)=sign(zm,u_m(ij,l))*( 0.5*(zq+zqp)*zsigp
588  s +(zsig-zsigp)*(zq+zz*(zqm-zq)) )
589  endif
590  endif
591 c if(zsig.lt.0.) then
592 c print*,'au point ij=',ij,' l=',l,' sig=',zsig
593 c stop
594 c endif
595  enddo
596  enddo
597 
598  do l=1,llm
599  do ij=iip1+iip1,ip1jm,iip1
600  u_mq(ij,l)=u_mq(ij-iim,l)
601  ladvplus(ij,l)=ladvplus(ij-iim,l)
602  enddo
603  enddo
604 
605 c=================================================================
606 C SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES
607 c=================================================================
608 c tris des regions a traiter
609  n0=0
610  do l=1,llm
611  nl(l)=0
612  do ij=iip2,ip1jm
613  if(ladvplus(ij,l)) then
614  nl(l)=nl(l)+1
615  u_mq(ij,l)=0.
616  endif
617  enddo
618  n0=n0+nl(l)
619  enddo
620 
621  if(n0.gt.1) then
622  IF (prt_level > 9) WRITE(lunout,*)
623  & 'Nombre de points pour lesquels on advect plus que le'
624  & ,'contenu de la maille : ',n0
625 
626  do l=1,llm
627  if(nl(l).gt.0) then
628  iju=0
629 c indicage des mailles concernees par le traitement special
630  do ij=iip2,ip1jm
631  if(ladvplus(ij,l).and.mod(ij,iip1).ne.0) then
632  iju=iju+1
633  indu(iju)=ij
634  endif
635  enddo
636  niju=iju
637 c print*,'niju,nl',niju,nl(l)
638 
639 c traitement des mailles
640  do iju=1,niju
641  ij=indu(iju)
642  j=(ij-1)/iip1+1
643  zu_m=u_m(ij,l)
644  u_mq(ij,l)=0.
645  if(zu_m.gt.0.) then
646  ijq=ij
647  i=ijq-(j-1)*iip1
648 c accumulation pour les mailles completements advectees
649  do while(zu_m.gt.masse(ijq,l))
650  u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
651  zu_m=zu_m-masse(ijq,l)
652  i=mod(i-2+iim,iim)+1
653  ijq=(j-1)*iip1+i
654  enddo
655 c MODIFS SPECIFIQUES DU SCHEMA
656 c ajout de la maille non completement advectee
657  zsig=zu_m/masse(ijq,l)
658  if(zsig.le.zsigd(ijq,l)) then
659  u_mq(ij,l)=u_mq(ij,l)+zu_m*(qd(ijq,l)
660  s -0.5*zsig/zsigd(ijq,l)*(qd(ijq,l)-q(ijq,l)))
661  else
662 c u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
663 c goto 8888
664  zz=0.5*(zsig-zsigd(ijq,l))/zsigg(ijq,l)
665  if(.not.(zz.gt.0..and.zz.le.0.5)) then
666  WRITE(lunout,*)'probleme2 au point ij=',ij,
667  s ' l=',l
668  WRITE(lunout,*)'zz=',zz
669  stop
670  endif
671  u_mq(ij,l)=u_mq(ij,l)+masse(ijq,l)*(
672  s 0.5*(q(ijq,l)+qd(ijq,l))*zsigd(ijq,l)
673  s +(zsig-zsigd(ijq,l))*(q(ijq,l)+zz*(qg(ijq,l)-q(ijq,l))) )
674  endif
675  else
676  ijq=ij+1
677  i=ijq-(j-1)*iip1
678 c accumulation pour les mailles completements advectees
679  do while(-zu_m.gt.masse(ijq,l))
680  u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
681  zu_m=zu_m+masse(ijq,l)
682  i=mod(i,iim)+1
683  ijq=(j-1)*iip1+i
684  enddo
685 c ajout de la maille non completement advectee
686 c 2eme MODIF SPECIFIQUE
687  zsig=-zu_m/masse(ij+1,l)
688  if(zsig.le.zsigg(ijq,l)) then
689  u_mq(ij,l)=u_mq(ij,l)+zu_m*(qg(ijq,l)
690  s -0.5*zsig/zsigg(ijq,l)*(qg(ijq,l)-q(ijq,l)))
691  else
692 c u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
693 c goto 9999
694  zz=0.5*(zsig-zsigg(ijq,l))/zsigd(ijq,l)
695  if(.not.(zz.gt.0..and.zz.le.0.5)) then
696  WRITE(lunout,*)'probleme22 au point ij=',ij
697  s ,' l=',l
698  WRITE(lunout,*)'zz=',zz
699  stop
700  endif
701  u_mq(ij,l)=u_mq(ij,l)-masse(ijq,l)*(
702  s 0.5*(q(ijq,l)+qg(ijq,l))*zsigg(ijq,l)
703  s +(zsig-zsigg(ijq,l))*
704  s(q(ijq,l)+zz*(qd(ijq,l)-q(ijq,l))) )
705  endif
706 c fin de la modif
707  endif
708  enddo
709  endif
710  enddo
711  endif ! n0.gt.0
712 
713 c bouclage en latitude
714  do l=1,llm
715  do ij=iip1+iip1,ip1jm,iip1
716  u_mq(ij,l)=u_mq(ij-iim,l)
717  enddo
718  enddo
719 
720 c=================================================================
721 c CALCUL DE LA CONVERGENCE DES FLUX
722 c=================================================================
723 
724  do l=1,llm
725  do ij=iip2+1,ip1jm
726  new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
727  q(ij,l)=(q(ij,l)*masse(ij,l)+
728  & u_mq(ij-1,l)-u_mq(ij,l))
729  & /new_m
730  masse(ij,l)=new_m
731  enddo
732 c Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
733  do ij=iip1+iip1,ip1jm,iip1
734  q(ij-iim,l)=q(ij,l)
735  masse(ij-iim,l)=masse(ij,l)
736  enddo
737  enddo
738 
739  RETURN
740  END
741  SUBROUTINE advny(q,qs,qn,masse,v_m)
742 c
743 c Auteur : F. Hourdin
744 c
745 c ********************************************************************
746 c Shema d'advection " pseudo amont " .
747 c ********************************************************************
748 c nq,iq,q,pbaru,pbarv,w sont des arguments d'entree pour le s-pg ....
749 c
750 c
751 c --------------------------------------------------------------------
752  IMPLICIT NONE
753 c
754 #include "dimensions.h"
755 #include "paramet.h"
756 #include "comgeom.h"
757 #include "iniprint.h"
758 c
759 c
760 c Arguments:
761 c ----------
762  real masse(ip1jmp1,llm)
763  real v_m( ip1jm,llm )
764  real q(ip1jmp1,llm),qn(ip1jmp1,llm),qs(ip1jmp1,llm)
765 c
766 c Local
767 c ---------
768 c
769  INTEGER ij,l
770 c
771  real new_m,zdq,zz
772  real zsigs(ip1jmp1),zsign(ip1jmp1),zsig
773  real v_mq(ip1jm,llm)
774  real convpn,convps,convmpn,convmps,massen,masses
775  real zm,zq,zsigm,zsigp,zqm,zqp
776  real ssum
777  real prec
778  save prec
779 
780 #ifdef CRAY
781  data prec/1.e-24/
782 #else
783  data prec/1.e-15/
784 #endif
785  do l=1,llm
786  do ij=1,ip1jmp1
787  zdq=qn(ij,l)-qs(ij,l)
788 c if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) then
789 c print*,'probleme au point ij=',ij,' l=',l,' advnqx'
790 c print*,qn(ij,l),q(ij,l),qs(ij,l)
791 c qn(ij,l)=q(ij,l)
792 c qs(ij,l)=q(ij,l)
793 c endif
794  if(abs(zdq).gt.prec) then
795  zsign(ij)=(q(ij,l)-qs(ij,l))/zdq
796  zsigs(ij)=1.-zsign(ij)
797 c if(.not.(zsign(ij).ge.0..and.zsign(ij).le.1. .and.
798 c s zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) then
799 c print*,'probleme au point ij=',ij,' l=',l
800 c print*,'sigs=',zsigs(ij),' sign=',zsign(ij)
801 c stop
802 c endif
803  else
804  zsign(ij)=0.5
805  zsigs(ij)=0.5
806  endif
807  enddo
808 
809 c calcul de la pente maximum dans la maille en valeur absolue
810 
811  do ij=1,ip1jm
812  if (v_m(ij,l).ge.0.) then
813  zsigp=zsign(ij+iip1)
814  zsigm=zsigs(ij+iip1)
815  zqp=qn(ij+iip1,l)
816  zqm=qs(ij+iip1,l)
817  zm=masse(ij+iip1,l)
818  zq=q(ij+iip1,l)
819  else
820  zsigm=zsign(ij)
821  zsigp=zsigs(ij)
822  zqm=qn(ij,l)
823  zqp=qs(ij,l)
824  zm=masse(ij,l)
825  zq=q(ij,l)
826  endif
827  zsig=abs(v_m(ij,l))/zm
828  if(zsig.eq.0.) zsigp=0.1
829  if (zsig.le.zsigp) then
830  v_mq(ij,l)=v_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
831  else
832  zz=0.5*(zsig-zsigp)/zsigm
833  v_mq(ij,l)=sign(zm,v_m(ij,l))*( 0.5*(zq+zqp)*zsigp
834  s +(zsig-zsigp)*(zq+zz*(zqm-zq)) )
835  endif
836  enddo
837  enddo
838 
839  do l=1,llm
840  do ij=iip2,ip1jm
841  new_m=masse(ij,l)
842  & +v_m(ij,l)-v_m(ij-iip1,l)
843  q(ij,l)=(q(ij,l)*masse(ij,l)+v_mq(ij,l)-v_mq(ij-iip1,l))
844  & /new_m
845  masse(ij,l)=new_m
846  enddo
847 c.-. ancienne version
848  convpn=ssum(iim,v_mq(1,l),1)
849  convmpn=ssum(iim,v_m(1,l),1)
850  massen=ssum(iim,masse(1,l),1)
851  new_m=massen+convmpn
852  q(1,l)=(q(1,l)*massen+convpn)/new_m
853  do ij = 1,iip1
854  q(ij,l)=q(1,l)
855  masse(ij,l)=new_m*aire(ij)/apoln
856  enddo
857 
858  convps=-ssum(iim,v_mq(ip1jm-iim,l),1)
859  convmps=-ssum(iim,v_m(ip1jm-iim,l),1)
860  masses=ssum(iim,masse(ip1jm+1,l),1)
861  new_m=masses+convmps
862  q(ip1jm+1,l)=(q(ip1jm+1,l)*masses+convps)/new_m
863  do ij = ip1jm+1,ip1jmp1
864  q(ij,l)=q(ip1jm+1,l)
865  masse(ij,l)=new_m*aire(ij)/apols
866  enddo
867  enddo
868 
869  RETURN
870  END
871  SUBROUTINE advnz(q,qh,qb,masse,w_m)
872 c
873 c Auteurs: F.Hourdin
874 c
875 c ********************************************************************
876 c Shema d'advection " pseudo amont " .
877 c b designe le bas et h le haut
878 c il y a une correspondance entre le b en z et le d en x
879 c ********************************************************************
880 c
881 c
882 c --------------------------------------------------------------------
883  IMPLICIT NONE
884 c
885 #include "dimensions.h"
886 #include "paramet.h"
887 #include "comgeom.h"
888 #include "iniprint.h"
889 c
890 c
891 c Arguments:
892 c ----------
893  real masse(ip1jmp1,llm)
894  real w_m( ip1jmp1,llm+1)
895  real q(ip1jmp1,llm),qb(ip1jmp1,llm),qh(ip1jmp1,llm)
896 
897 c
898 c Local
899 c ---------
900 c
901  INTEGER ij,l
902 c
903  real new_m,zdq,zz
904  real zsigh(ip1jmp1,llm),zsigb(ip1jmp1,llm),zsig
905  real w_mq(ip1jmp1,llm+1)
906  real zm,zq,zsigm,zsigp,zqm,zqp
907  real prec
908  save prec
909 
910 #ifdef CRAY
911  data prec/1.e-24/
912 #else
913  data prec/1.e-13/
914 #endif
915 
916  do l=1,llm
917  do ij=1,ip1jmp1
918  zdq=qb(ij,l)-qh(ij,l)
919 c if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) then
920 c print*,'probleme au point ij=',ij,' l=',l
921 c print*,qh(ij,l),q(ij,l),qb(ij,l)
922 c qh(ij,l)=q(ij,l)
923 c qb(ij,l)=q(ij,l)
924 c endif
925 
926  if(abs(zdq).gt.prec) then
927  zsigb(ij,l)=(q(ij,l)-qh(ij,l))/zdq
928  zsigh(ij,l)=1.-zsigb(ij,l)
929  zsigb(ij,l)=min(max(zsigb(ij,l),0.),1.)
930  else
931  zsigb(ij,l)=0.5
932  zsigh(ij,l)=0.5
933  endif
934  enddo
935  enddo
936 
937 c print*,'ok1'
938 c calcul de la pente maximum dans la maille en valeur absolue
939  do l=2,llm
940  do ij=1,ip1jmp1
941  if (w_m(ij,l).ge.0.) then
942  zsigp=zsigb(ij,l)
943  zsigm=zsigh(ij,l)
944  zqp=qb(ij,l)
945  zqm=qh(ij,l)
946  zm=masse(ij,l)
947  zq=q(ij,l)
948  else
949  zsigm=zsigb(ij,l-1)
950  zsigp=zsigh(ij,l-1)
951  zqm=qb(ij,l-1)
952  zqp=qh(ij,l-1)
953  zm=masse(ij,l-1)
954  zq=q(ij,l-1)
955  endif
956  zsig=abs(w_m(ij,l))/zm
957  if(zsig.eq.0.) zsigp=0.1
958  if (zsig.le.zsigp) then
959  w_mq(ij,l)=w_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
960  else
961  zz=0.5*(zsig-zsigp)/zsigm
962  w_mq(ij,l)=sign(zm,w_m(ij,l))*( 0.5*(zq+zqp)*zsigp
963  s +(zsig-zsigp)*(zq+zz*(zqm-zq)) )
964  endif
965  enddo
966  enddo
967 
968  do ij=1,ip1jmp1
969  w_mq(ij,llm+1)=0.
970  w_mq(ij,1)=0.
971  enddo
972 
973  do l=1,llm
974  do ij=1,ip1jmp1
975  new_m=masse(ij,l)+w_m(ij,l+1)-w_m(ij,l)
976  q(ij,l)=(q(ij,l)*masse(ij,l)+w_mq(ij,l+1)-w_mq(ij,l))
977  & /new_m
978  masse(ij,l)=new_m
979  enddo
980  enddo
981 c print*,'ok3'
982  RETURN
983  END