source: LMDZ4/trunk/libf/dyn3d/vlsplt.F @ 962

Last change on this file since 962 was 595, checked in by Laurent Fairhead, 20 years ago

Si il n'y avait qu'un point ou on advectait plus que le contenu de la maille,
il n'etait jamais traite YM
LF

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