source: LMDZ5/branches/LMDZ5_AR5/libf/dyn3d/vlsplt.F @ 4116

Last change on this file since 4116 was 1550, checked in by lguez, 13 years ago

Bug fix in "bilan_dyn_p". The index was out of bounds in the removed
assignment . Also, the removed assignment was useless.

Bug fix in "coefkzmin". The size of a dummy array cannot exceed the
size of the associated actual array. ("coefkzmin" is called by
"coef_diff_turb".) "km(:, klev+1)" and "kn(:, klev+1)" were not
defined in "coefkzmin" so this was maybe an innocuous bug.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.2 KB
Line 
1c
2c $Id: vlsplt.F 1550 2011-07-05 09:44:55Z idelkadi $
3c
4
5      SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt)
6c
7c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
8c
9c    ********************************************************************
10c     Shema  d'advection " pseudo amont " .
11c    ********************************************************************
12c     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
13c
14c   pente_max facteur de limitation des pentes: 2 en general
15c                                               0 pour un schema amont
16c   pbaru,pbarv,w flux de masse en u ,v ,w
17c   pdt pas de temps
18c
19c   --------------------------------------------------------------------
20      IMPLICIT NONE
21c
22#include "dimensions.h"
23#include "paramet.h"
24#include "logic.h"
25#include "comvert.h"
26#include "comconst.h"
27
28c
29c   Arguments:
30c   ----------
31      REAL masse(ip1jmp1,llm),pente_max
32c      REAL masse(iip1,jjp1,llm),pente_max
33      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
34      REAL q(ip1jmp1,llm)
35c      REAL q(iip1,jjp1,llm)
36      REAL w(ip1jmp1,llm),pdt
37c
38c      Local
39c   ---------
40c
41      INTEGER i,ij,l,j,ii
42      INTEGER ijlqmin,iqmin,jqmin,lqmin
43c
44      REAL zm(ip1jmp1,llm),newmasse
45      REAL mu(ip1jmp1,llm)
46      REAL mv(ip1jm,llm)
47      REAL mw(ip1jmp1,llm+1)
48      REAL zq(ip1jmp1,llm),zz
49      REAL dqx(ip1jmp1,llm),dqy(ip1jmp1,llm),dqz(ip1jmp1,llm)
50      REAL second,temps0,temps1,temps2,temps3
51      REAL ztemps1,ztemps2,ztemps3
52      REAL zzpbar, zzw
53      LOGICAL testcpu
54      SAVE testcpu
55      SAVE temps1,temps2,temps3
56      INTEGER iminn,imaxx
57
58      REAL qmin,qmax
59      DATA qmin,qmax/0.,1.e33/
60      DATA testcpu/.false./
61      DATA temps1,temps2,temps3/0.,0.,0./
62
63
64        zzpbar = 0.5 * pdt
65        zzw    = pdt
66      DO l=1,llm
67        DO ij = iip2,ip1jm
68            mu(ij,l)=pbaru(ij,l) * zzpbar
69         ENDDO
70         DO ij=1,ip1jm
71            mv(ij,l)=pbarv(ij,l) * zzpbar
72         ENDDO
73         DO ij=1,ip1jmp1
74            mw(ij,l)=w(ij,l) * zzw
75         ENDDO
76      ENDDO
77
78      DO ij=1,ip1jmp1
79         mw(ij,llm+1)=0.
80      ENDDO
81     
82      CALL SCOPY(ijp1llm,q,1,zq,1)
83      CALL SCOPY(ijp1llm,masse,1,zm,1)
84
85cprint*,'Entree vlx1'
86c       call minmaxq(zq,qmin,qmax,'avant vlx     ')
87      call vlx(zq,pente_max,zm,mu)
88cprint*,'Sortie vlx1'
89c       call minmaxq(zq,qmin,qmax,'apres vlx1    ')
90
91c print*,'Entree vly1'
92      call vly(zq,pente_max,zm,mv)
93c       call minmaxq(zq,qmin,qmax,'apres vly1     ')
94cprint*,'Sortie vly1'
95      call vlz(zq,pente_max,zm,mw)
96c       call minmaxq(zq,qmin,qmax,'apres vlz     ')
97
98
99      call vly(zq,pente_max,zm,mv)
100c       call minmaxq(zq,qmin,qmax,'apres vly     ')
101
102
103      call vlx(zq,pente_max,zm,mu)
104c       call minmaxq(zq,qmin,qmax,'apres vlx2    ')
105       
106
107      DO l=1,llm
108         DO ij=1,ip1jmp1
109           q(ij,l)=zq(ij,l)
110         ENDDO
111         DO ij=1,ip1jm+1,iip1
112            q(ij+iim,l)=q(ij,l)
113         ENDDO
114      ENDDO
115
116      RETURN
117      END
118      SUBROUTINE vlx(q,pente_max,masse,u_m)
119
120c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
121c
122c    ********************************************************************
123c     Shema  d'advection " pseudo amont " .
124c    ********************************************************************
125c     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
126c
127c
128c   --------------------------------------------------------------------
129      IMPLICIT NONE
130c
131      include "dimensions.h"
132      include "paramet.h"
133      include "logic.h"
134      include "comvert.h"
135      include "comconst.h"
136      include "iniprint.h"
137c
138c
139c   Arguments:
140c   ----------
141      REAL masse(ip1jmp1,llm),pente_max
142      REAL u_m( ip1jmp1,llm ),pbarv( iip1,jjm,llm)
143      REAL q(ip1jmp1,llm)
144      REAL w(ip1jmp1,llm)
145c
146c      Local
147c   ---------
148c
149      INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
150      INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
151c
152      REAL new_m,zu_m,zdum(ip1jmp1,llm)
153      REAL sigu(ip1jmp1),dxq(ip1jmp1,llm),dxqu(ip1jmp1)
154      REAL zz(ip1jmp1)
155      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
156      REAL u_mq(ip1jmp1,llm)
157
158      Logical extremum,first,testcpu
159      SAVE first,testcpu
160
161      REAL      SSUM
162      REAL temps0,temps1,temps2,temps3,temps4,temps5,second
163      SAVE temps0,temps1,temps2,temps3,temps4,temps5
164
165      REAL z1,z2,z3
166
167      DATA first,testcpu/.true.,.false./
168
169      IF(first) THEN
170         temps1=0.
171         temps2=0.
172         temps3=0.
173         temps4=0.
174         temps5=0.
175         first=.false.
176      ENDIF
177
178c   calcul de la pente a droite et a gauche de la maille
179
180
181      IF (pente_max.gt.-1.e-5) THEN
182c       IF (pente_max.gt.10) THEN
183
184c   calcul des pentes avec limitation, Van Leer scheme I:
185c   -----------------------------------------------------
186
187c   calcul de la pente aux points u
188         DO l = 1, llm
189            DO ij=iip2,ip1jm-1
190               dxqu(ij)=q(ij+1,l)-q(ij,l)
191c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
192c              sigu(ij)=u_m(ij,l)/masse(ij,l)
193            ENDDO
194            DO ij=iip1+iip1,ip1jm,iip1
195               dxqu(ij)=dxqu(ij-iim)
196c              sigu(ij)=sigu(ij-iim)
197            ENDDO
198
199            DO ij=iip2,ip1jm
200               adxqu(ij)=abs(dxqu(ij))
201            ENDDO
202
203c   calcul de la pente maximum dans la maille en valeur absolue
204
205            DO ij=iip2+1,ip1jm
206               dxqmax(ij,l)=pente_max*
207     ,      min(adxqu(ij-1),adxqu(ij))
208c limitation subtile
209c    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
210         
211
212            ENDDO
213
214            DO ij=iip1+iip1,ip1jm,iip1
215               dxqmax(ij-iim,l)=dxqmax(ij,l)
216            ENDDO
217
218            DO ij=iip2+1,ip1jm
219#ifdef CRAY
220               dxq(ij,l)=
221     ,         cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
222#else
223               IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
224                  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
225               ELSE
226c   extremum local
227                  dxq(ij,l)=0.
228               ENDIF
229#endif
230               dxq(ij,l)=0.5*dxq(ij,l)
231               dxq(ij,l)=
232     ,         sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
233            ENDDO
234
235         ENDDO ! l=1,llm
236cprint*,'Ok calcul des pentes'
237
238      ELSE ! (pente_max.lt.-1.e-5)
239
240c   Pentes produits:
241c   ----------------
242
243         DO l = 1, llm
244            DO ij=iip2,ip1jm-1
245               dxqu(ij)=q(ij+1,l)-q(ij,l)
246            ENDDO
247            DO ij=iip1+iip1,ip1jm,iip1
248               dxqu(ij)=dxqu(ij-iim)
249            ENDDO
250
251            DO ij=iip2+1,ip1jm
252               zz(ij)=dxqu(ij-1)*dxqu(ij)
253               zz(ij)=zz(ij)+zz(ij)
254               IF(zz(ij).gt.0) THEN
255                  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
256               ELSE
257c   extremum local
258                  dxq(ij,l)=0.
259               ENDIF
260            ENDDO
261
262         ENDDO
263
264      ENDIF ! (pente_max.lt.-1.e-5)
265
266c   bouclage de la pente en iip1:
267c   -----------------------------
268
269      DO l=1,llm
270         DO ij=iip1+iip1,ip1jm,iip1
271            dxq(ij-iim,l)=dxq(ij,l)
272         ENDDO
273         DO ij=1,ip1jmp1
274            iadvplus(ij,l)=0
275         ENDDO
276
277      ENDDO
278
279c print*,'Bouclage en iip1'
280
281c   calcul des flux a gauche et a droite
282
283#ifdef CRAY
284
285      DO l=1,llm
286       DO ij=iip2,ip1jm-1
287          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
288     ,                     1.+u_m(ij,l)/masse(ij+1,l),
289     ,                     u_m(ij,l))
290          zdum(ij,l)=0.5*zdum(ij,l)
291          u_mq(ij,l)=cvmgp(
292     ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
293     ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
294     ,                u_m(ij,l))
295          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
296       ENDDO
297      ENDDO
298#else
299c   on cumule le flux correspondant a toutes les mailles dont la masse
300c   au travers de la paroi pENDant le pas de temps.
301cprint*,'Cumule ....'
302
303      DO l=1,llm
304       DO ij=iip2,ip1jm-1
305c       print*,'masse(',ij,')=',masse(ij,l)
306          IF (u_m(ij,l).gt.0.) THEN
307             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
308             u_mq(ij,l)=u_m(ij,l)*(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l))
309          ELSE
310             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
311             u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l))
312          ENDIF
313       ENDDO
314      ENDDO
315#endif
316c       stop
317
318c       go to 9999
319c   detection des points ou on advecte plus que la masse de la
320c   maille
321      DO l=1,llm
322         DO ij=iip2,ip1jm-1
323            IF(zdum(ij,l).lt.0) THEN
324               iadvplus(ij,l)=1
325               u_mq(ij,l)=0.
326            ENDIF
327         ENDDO
328      ENDDO
329cprint*,'Ok test 1'
330      DO l=1,llm
331       DO ij=iip1+iip1,ip1jm,iip1
332          iadvplus(ij,l)=iadvplus(ij-iim,l)
333       ENDDO
334      ENDDO
335c print*,'Ok test 2'
336
337
338c   traitement special pour le cas ou on advecte en longitude plus que le
339c   contenu de la maille.
340c   cette partie est mal vectorisee.
341
342c  calcul du nombre de maille sur lequel on advecte plus que la maille.
343
344      n0=0
345      DO l=1,llm
346         nl(l)=0
347         DO ij=iip2,ip1jm
348            nl(l)=nl(l)+iadvplus(ij,l)
349         ENDDO
350         n0=n0+nl(l)
351      ENDDO
352
353      IF(n0.gt.0) THEN
354      if (prt_level > 2) PRINT *,
355     $        '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,appn,apps,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     appn=abs(dyq(1)/dyqv(iip1+1))
605C     PRINT*,dyq(ip1jm+1)
606C     PRINT*,dyqv(ip1jm-iip1+1)
607C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
608C     DO ij=2,iim
609C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
610C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
611C     ENDDO
612C     appn=min(pente_max/appn,1.)
613C     apps=min(pente_max/apps,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    &   appn=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    &   apps=0.
623C
624C   limitation des pentes aux poles
625C     DO ij=1,iip1
626C        dyq(ij)=appn*dyq(ij)
627C        dyq(ip1jm+ij)=apps*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.