source: LMDZ6/branches/DYNAMICO-conv/libf/dyn3d_common/advn.F @ 3024

Last change on this file since 3024 was 2603, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn logic.h into module logic_mod.F90
EM

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