source: LMDZ4/trunk/libf/dyn3dpar/advn.F @ 1345

Last change on this file since 1345 was 774, checked in by Laurent Fairhead, 17 years ago

Suite du merge entre la version et la HEAD: quelques modifications de
Yann sur le

LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.9 KB
Line 
1!
2! $Header$
3!
4      SUBROUTINE advn(q,masse,w,pbaru,pbarv,pdt,mode)
5c
6c     Auteur : F. Hourdin
7c
8c    ********************************************************************
9c     Shema  d'advection " pseudo amont " .
10c    ********************************************************************
11c     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
12c
13c   pbaru,pbarv,w flux de masse en u ,v ,w
14c   pdt pas de temps
15c
16c   --------------------------------------------------------------------
17      IMPLICIT NONE
18c
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
27c
28c   Arguments:
29c   ----------
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
35c
36c      Local
37c   ---------
38c
39      INTEGER i,ij,l,j,ii
40      integer ijlqmin,iqmin,jqmin,lqmin
41      integer ismin
42c
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
111c     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)
118c     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)
123c     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)
145c
146c     Auteurs:   Calcul des valeurs de q aux point u.
147c
148c   --------------------------------------------------------------------
149      IMPLICIT NONE
150c
151#include "dimensions.h"
152#include "paramet.h"
153#include "iniprint.h"
154c
155c
156c   Arguments:
157c   ----------
158      real q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)
159c
160c      Local
161c   ---------
162c
163      INTEGER ij,l
164c
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
173c   calcul des pentes en u:
174c   -----------------------
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
205c   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
2568888     continue
257      enddo
258      endif
259      RETURN
260      END
261      SUBROUTINE advnqy(q,qs,qn)
262c
263c     Auteurs:   Calcul des valeurs de q aux point v.
264c
265c   --------------------------------------------------------------------
266      IMPLICIT NONE
267c
268#include "dimensions.h"
269#include "paramet.h"
270#include "iniprint.h"
271c
272c
273c   Arguments:
274c   ----------
275      real q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)
276c
277c      Local
278c   ---------
279c
280      INTEGER ij,l
281c
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
299c   calcul des pentes en u:
300c   -----------------------
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
315c 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
323c   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)
337c              if (.not.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
338c              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)
358c
359c     Auteurs:   Calcul des valeurs de q aux point v.
360c
361c   --------------------------------------------------------------------
362      IMPLICIT NONE
363c
364#include "dimensions.h"
365#include "paramet.h"
366#include "iniprint.h"
367c
368c
369c   Arguments:
370c   ----------
371      real q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)
372c
373c      Local
374c   ---------
375c
376      INTEGER ij,l
377c
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
387c   calcul des pentes en u:
388c   -----------------------
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
419c 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
427c   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
452c     do l=2,llm-1
453c        do ij=1,ip1jmp1
454c           if(extremum(ij,l)) then
455c              if (.not.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
456c              if (.not.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
457c           endif
458c        enddo
459c     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)
474c
475c     Auteur : F. Hourdin
476c
477c    ********************************************************************
478c     Shema  d'advection " pseudo amont " .
479c    ********************************************************************
480c     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
481c
482c
483c   --------------------------------------------------------------------
484      IMPLICIT NONE
485c
486#include "dimensions.h"
487#include "paramet.h"
488#include "logic.h"
489#include "comvert.h"
490#include "comconst.h"
491#include "iniprint.h"
492c
493c
494c   Arguments:
495c   ----------
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)
500c
501c      Local
502c   ---------
503c
504      INTEGER i,j,ij,l,indu(ip1jmp1),niju,iju,ijq
505      integer n0,nl(llm)
506c
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)
527c              if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) then
528c                 print*,'probleme au point ij=',ij,'  l=',l
529c                 print*,qd(ij,l),q(ij,l),qg(ij,l)
530c                 qd(ij,l)=q(ij,l)
531c                 qg(ij,l)=q(ij,l)
532c              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)
536c                 if(.not.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .and.
537c    s               zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) then
538c                    print*,'probleme au point ij=',ij,'  l=',l
539c                    print*,'sigg=',zsigg(ij,l),'  sigd=',zsigd(ij,l)
540c                    print*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq
541c                    stop
542c                 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
552c   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
591c         if(zsig.lt.0.) then
592c            print*,'au point ij=',ij,'  l=',l,'  sig=',zsig
593c            stop
594c         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
605c=================================================================
606C   SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES
607c=================================================================
608c   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
629c   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
637c              print*,'niju,nl',niju,nl(l)
638
639c  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
648c   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
655c   MODIFS SPECIFIQUES DU SCHEMA
656c   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
662c               u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
663c         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
678c   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
685c   ajout de la maille non completement advectee
686c 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
692c               u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
693c           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
706c   fin de la modif
707                  endif
708               enddo
709            endif
710         enddo
711      endif  ! n0.gt.0
712
713c   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
720c=================================================================
721c   CALCUL DE LA CONVERGENCE DES FLUX
722c=================================================================
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
732c   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)
742c
743c     Auteur : F. Hourdin
744c
745c    ********************************************************************
746c     Shema  d'advection " pseudo amont " .
747c    ********************************************************************
748c     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
749c
750c
751c   --------------------------------------------------------------------
752      IMPLICIT NONE
753c
754#include "dimensions.h"
755#include "paramet.h"
756#include "comgeom.h"
757#include "iniprint.h"
758c
759c
760c   Arguments:
761c   ----------
762      real masse(ip1jmp1,llm)
763      real v_m( ip1jm,llm )
764      real q(ip1jmp1,llm),qn(ip1jmp1,llm),qs(ip1jmp1,llm)
765c
766c      Local
767c   ---------
768c
769      INTEGER ij,l
770c
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)
788c              if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) then
789c                 print*,'probleme au point ij=',ij,'  l=',l,'  advnqx'
790c                 print*,qn(ij,l),q(ij,l),qs(ij,l)
791c                 qn(ij,l)=q(ij,l)
792c                 qs(ij,l)=q(ij,l)
793c              endif
794               if(abs(zdq).gt.prec) then
795                  zsign(ij)=(q(ij,l)-qs(ij,l))/zdq
796                  zsigs(ij)=1.-zsign(ij)
797c                 if(.not.(zsign(ij).ge.0..and.zsign(ij).le.1. .and.
798c    s               zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) then
799c                    print*,'probleme au point ij=',ij,'  l=',l
800c                    print*,'sigs=',zsigs(ij),'  sign=',zsign(ij)
801c                    stop
802c                 endif
803               else
804                  zsign(ij)=0.5
805                  zsigs(ij)=0.5
806               endif
807            enddo
808
809c   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
847c.-. 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)
872c
873c     Auteurs:   F.Hourdin
874c
875c    ********************************************************************
876c     Shema  d'advection " pseudo amont " .
877c     b designe le bas et h le haut
878c     il y a une correspondance entre le b en z et le d en x
879c    ********************************************************************
880c
881c
882c   --------------------------------------------------------------------
883      IMPLICIT NONE
884c
885#include "dimensions.h"
886#include "paramet.h"
887#include "comgeom.h"
888#include "iniprint.h"
889c
890c
891c   Arguments:
892c   ----------
893      real masse(ip1jmp1,llm)
894      real w_m( ip1jmp1,llm+1)
895      real q(ip1jmp1,llm),qb(ip1jmp1,llm),qh(ip1jmp1,llm)
896
897c
898c      Local
899c   ---------
900c
901      INTEGER ij,l
902c
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)
919c              if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) then
920c                 print*,'probleme au point ij=',ij,'  l=',l
921c                 print*,qh(ij,l),q(ij,l),qb(ij,l)
922c                 qh(ij,l)=q(ij,l)
923c                 qb(ij,l)=q(ij,l)
924c              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
937c      print*,'ok1'
938c   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
981c     print*,'ok3'
982      RETURN
983      END
Note: See TracBrowser for help on using the repository browser.