source: LMDZ4/trunk/libf/dyn3dpar/vlspltqs.F @ 701

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

Import d'une version parallele de la dynamique YM
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 21.1 KB
RevLine 
[630]1!
2! $Header$
3!
4       SUBROUTINE vlspltqs ( q,pente_max,masse,w,pbaru,pbarv,pdt,
5     ,                                  p,pk,teta                 )
6c
7c     Auteurs:   P.Le Van, F.Hourdin, F.Forget, F.Codron
8c
9c    ********************************************************************
10c          Shema  d'advection " pseudo amont " .
11c      + test sur humidite specifique: Q advecte< Qsat aval
12c                   (F. Codron, 10/99)
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     teta temperature potentielle, p pression aux interfaces,
22c     pk exner au milieu des couches necessaire pour calculer Qsat
23c   --------------------------------------------------------------------
24      IMPLICIT NONE
25c
26#include "dimensions.h"
27#include "paramet.h"
28#include "logic.h"
29#include "comvert.h"
30#include "comconst.h"
31
32c
33c   Arguments:
34c   ----------
35      REAL masse(ip1jmp1,llm),pente_max
36      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
37      REAL q(ip1jmp1,llm)
38      REAL w(ip1jmp1,llm),pdt
39      REAL p(ip1jmp1,llmp1),teta(ip1jmp1,llm),pk(ip1jmp1,llm)
40c
41c      Local
42c   ---------
43c
44      INTEGER i,ij,l,j,ii
45c
46      REAL qsat(ip1jmp1,llm)
47      REAL zm(ip1jmp1,llm)
48      REAL mu(ip1jmp1,llm)
49      REAL mv(ip1jm,llm)
50      REAL mw(ip1jmp1,llm+1)
51      REAL zq(ip1jmp1,llm)
52      REAL temps1,temps2,temps3
53      REAL zzpbar, zzw
54      LOGICAL testcpu
55      SAVE testcpu
56      SAVE temps1,temps2,temps3
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
63c--pour rapport de melange saturant--
64
65      REAL rtt,retv,r2es,r3les,r3ies,r4les,r4ies,play
66      REAL ptarg,pdelarg,foeew,zdelta
67      REAL tempe(ip1jmp1)
68
69c    fonction psat(T)
70
71       FOEEW ( PTARG,PDELARG ) = EXP (
72     *          (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT)
73     * / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) )
74
75        r2es  = 380.11733
76        r3les = 17.269
77        r3ies = 21.875
78        r4les = 35.86
79        r4ies = 7.66
80        retv = 0.6077667
81        rtt  = 273.16
82
83c-- Calcul de Qsat en chaque point
84c-- approximation: au milieu des couches play(l)=(p(l)+p(l+1))/2
85c   pour eviter une exponentielle.
86        DO l = 1, llm
87         DO ij = 1, ip1jmp1
88          tempe(ij) = teta(ij,l) * pk(ij,l) /cpp
89         ENDDO
90         DO ij = 1, ip1jmp1
91          zdelta = MAX( 0., SIGN(1., rtt - tempe(ij)) )
92          play   = 0.5*(p(ij,l)+p(ij,l+1))
93          qsat(ij,l) = MIN(0.5, r2es* FOEEW(tempe(ij),zdelta) / play )
94          qsat(ij,l) = qsat(ij,l) / ( 1. - retv * qsat(ij,l) )
95         ENDDO
96        ENDDO
97
98c      PRINT*,'Debut vlsplt version debug sans vlyqs'
99
100        zzpbar = 0.5 * pdt
101        zzw    = pdt
102      DO l=1,llm
103        DO ij = iip2,ip1jm
104            mu(ij,l)=pbaru(ij,l) * zzpbar
105         ENDDO
106         DO ij=1,ip1jm
107            mv(ij,l)=pbarv(ij,l) * zzpbar
108         ENDDO
109         DO ij=1,ip1jmp1
110            mw(ij,l)=w(ij,l) * zzw
111         ENDDO
112      ENDDO
113
114      DO ij=1,ip1jmp1
115         mw(ij,llm+1)=0.
116      ENDDO
117
118      CALL SCOPY(ijp1llm,q,1,zq,1)
119      CALL SCOPY(ijp1llm,masse,1,zm,1)
120
121c      call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
122      call vlxqs(zq,pente_max,zm,mu,qsat)
123
124
125c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
126
127      call vlyqs(zq,pente_max,zm,mv,qsat)
128
129
130c      call minmaxq(zq,qmin,qmax,'avant vlz     ')
131
132      call vlz(zq,pente_max,zm,mw)
133
134
135c     call minmaxq(zq,qmin,qmax,'avant vlyqs     ')
136c     call minmaxq(zm,qmin,qmax,'M avant vlyqs     ')
137
138      call vlyqs(zq,pente_max,zm,mv,qsat)
139
140
141c     call minmaxq(zq,qmin,qmax,'avant vlxqs     ')
142c     call minmaxq(zm,qmin,qmax,'M avant vlxqs     ')
143
144      call vlxqs(zq,pente_max,zm,mu,qsat)
145
146c     call minmaxq(zq,qmin,qmax,'apres vlxqs     ')
147c     call minmaxq(zm,qmin,qmax,'M apres vlxqs     ')
148
149
150      DO l=1,llm
151         DO ij=1,ip1jmp1
152           q(ij,l)=zq(ij,l)
153         ENDDO
154         DO ij=1,ip1jm+1,iip1
155            q(ij+iim,l)=q(ij,l)
156         ENDDO
157      ENDDO
158
159      RETURN
160      END
161      SUBROUTINE vlxqs(q,pente_max,masse,u_m,qsat)
162c
163c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
164c
165c    ********************************************************************
166c     Shema  d'advection " pseudo amont " .
167c    ********************************************************************
168c
169c   --------------------------------------------------------------------
170      IMPLICIT NONE
171c
172#include "dimensions.h"
173#include "paramet.h"
174#include "logic.h"
175#include "comvert.h"
176#include "comconst.h"
177c
178c
179c   Arguments:
180c   ----------
181      REAL masse(ip1jmp1,llm),pente_max
182      REAL u_m( ip1jmp1,llm )
183      REAL q(ip1jmp1,llm)
184      REAL qsat(ip1jmp1,llm)
185c
186c      Local
187c   ---------
188c
189      INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
190      INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
191c
192      REAL new_m,zu_m,zdum(ip1jmp1,llm)
193      REAL dxq(ip1jmp1,llm),dxqu(ip1jmp1)
194      REAL zz(ip1jmp1)
195      REAL adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
196      REAL u_mq(ip1jmp1,llm)
197
198      Logical first,testcpu
199      SAVE first,testcpu
200
201      REAL      SSUM
202      REAL temps0,temps1,temps2,temps3,temps4,temps5
203      SAVE temps0,temps1,temps2,temps3,temps4,temps5
204
205
206      DATA first,testcpu/.true.,.false./
207
208      IF(first) THEN
209         temps1=0.
210         temps2=0.
211         temps3=0.
212         temps4=0.
213         temps5=0.
214         first=.false.
215      ENDIF
216
217c   calcul de la pente a droite et a gauche de la maille
218
219
220      IF (pente_max.gt.-1.e-5) THEN
221c     IF (pente_max.gt.10) THEN
222
223c   calcul des pentes avec limitation, Van Leer scheme I:
224c   -----------------------------------------------------
225
226c   calcul de la pente aux points u
227         DO l = 1, llm
228            DO ij=iip2,ip1jm-1
229               dxqu(ij)=q(ij+1,l)-q(ij,l)
230c              IF(u_m(ij,l).lt.0.) stop'limx n admet pas les U<0'
231c              sigu(ij)=u_m(ij,l)/masse(ij,l)
232            ENDDO
233            DO ij=iip1+iip1,ip1jm,iip1
234               dxqu(ij)=dxqu(ij-iim)
235c              sigu(ij)=sigu(ij-iim)
236            ENDDO
237
238            DO ij=iip2,ip1jm
239               adxqu(ij)=abs(dxqu(ij))
240            ENDDO
241
242c   calcul de la pente maximum dans la maille en valeur absolue
243
244            DO ij=iip2+1,ip1jm
245               dxqmax(ij,l)=pente_max*
246     ,      min(adxqu(ij-1),adxqu(ij))
247c limitation subtile
248c    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
249         
250
251            ENDDO
252
253            DO ij=iip1+iip1,ip1jm,iip1
254               dxqmax(ij-iim,l)=dxqmax(ij,l)
255            ENDDO
256
257            DO ij=iip2+1,ip1jm
258#ifdef CRAY
259               dxq(ij,l)=
260     ,         cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
261#else
262               IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
263                  dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
264               ELSE
265c   extremum local
266                  dxq(ij,l)=0.
267               ENDIF
268#endif
269               dxq(ij,l)=0.5*dxq(ij,l)
270               dxq(ij,l)=
271     ,         sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
272            ENDDO
273
274         ENDDO ! l=1,llm
275
276      ELSE ! (pente_max.lt.-1.e-5)
277
278c   Pentes produits:
279c   ----------------
280
281         DO l = 1, llm
282            DO ij=iip2,ip1jm-1
283               dxqu(ij)=q(ij+1,l)-q(ij,l)
284            ENDDO
285            DO ij=iip1+iip1,ip1jm,iip1
286               dxqu(ij)=dxqu(ij-iim)
287            ENDDO
288
289            DO ij=iip2+1,ip1jm
290               zz(ij)=dxqu(ij-1)*dxqu(ij)
291               zz(ij)=zz(ij)+zz(ij)
292               IF(zz(ij).gt.0) THEN
293                  dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
294               ELSE
295c   extremum local
296                  dxq(ij,l)=0.
297               ENDIF
298            ENDDO
299
300         ENDDO
301
302      ENDIF ! (pente_max.lt.-1.e-5)
303
304c   bouclage de la pente en iip1:
305c   -----------------------------
306
307      DO l=1,llm
308         DO ij=iip1+iip1,ip1jm,iip1
309            dxq(ij-iim,l)=dxq(ij,l)
310         ENDDO
311
312         DO ij=1,ip1jmp1
313            iadvplus(ij,l)=0
314         ENDDO
315
316      ENDDO
317
318
319c   calcul des flux a gauche et a droite
320
321#ifdef CRAY
322c--pas encore modification sur Qsat
323      DO l=1,llm
324       DO ij=iip2,ip1jm-1
325          zdum(ij,l)=cvmgp(1.-u_m(ij,l)/masse(ij,l),
326     ,                     1.+u_m(ij,l)/masse(ij+1,l),
327     ,                     u_m(ij,l))
328          zdum(ij,l)=0.5*zdum(ij,l)
329          u_mq(ij,l)=cvmgp(
330     ,                q(ij,l)+zdum(ij,l)*dxq(ij,l),
331     ,                q(ij+1,l)-zdum(ij,l)*dxq(ij+1,l),
332     ,                u_m(ij,l))
333          u_mq(ij,l)=u_m(ij,l)*u_mq(ij,l)
334       ENDDO
335      ENDDO
336#else
337c   on cumule le flux correspondant a toutes les mailles dont la masse
338c   au travers de la paroi pENDant le pas de temps.
339c   le rapport de melange de l'air advecte est min(q_vanleer, Qsat_downwind)
340      DO l=1,llm
341       DO ij=iip2,ip1jm-1
342          IF (u_m(ij,l).gt.0.) THEN
343             zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l)
344             u_mq(ij,l)=u_m(ij,l)*
345     $         min(q(ij,l)+0.5*zdum(ij,l)*dxq(ij,l),qsat(ij+1,l))
346          ELSE
347             zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l)
348             u_mq(ij,l)=u_m(ij,l)*
349     $         min(q(ij+1,l)-0.5*zdum(ij,l)*dxq(ij+1,l),qsat(ij,l))
350          ENDIF
351       ENDDO
352      ENDDO
353#endif
354
355
356c   detection des points ou on advecte plus que la masse de la
357c   maille
358      DO l=1,llm
359         DO ij=iip2,ip1jm-1
360            IF(zdum(ij,l).lt.0) THEN
361               iadvplus(ij,l)=1
362               u_mq(ij,l)=0.
363            ENDIF
364         ENDDO
365      ENDDO
366      DO l=1,llm
367       DO ij=iip1+iip1,ip1jm,iip1
368          iadvplus(ij,l)=iadvplus(ij-iim,l)
369       ENDDO
370      ENDDO
371
372
373
374c   traitement special pour le cas ou on advecte en longitude plus que le
375c   contenu de la maille.
376c   cette partie est mal vectorisee.
377
378c   pas d'influence de la pression saturante (pour l'instant)
379
380c  calcul du nombre de maille sur lequel on advecte plus que la maille.
381
382      n0=0
383      DO l=1,llm
384         nl(l)=0
385         DO ij=iip2,ip1jm
386            nl(l)=nl(l)+iadvplus(ij,l)
387         ENDDO
388         n0=n0+nl(l)
389      ENDDO
390
391cym      IF(n0.gt.1) THEN
392         IF(n0.gt.0) THEN
393ccc      PRINT*,'Nombre de points pour lesquels on advect plus que le'
394ccc     &       ,'contenu de la maille : ',n0
395
396         DO l=1,llm
397            IF(nl(l).gt.0) THEN
398               iju=0
399c   indicage des mailles concernees par le traitement special
400               DO ij=iip2,ip1jm
401                  IF(iadvplus(ij,l).eq.1.and.mod(ij,iip1).ne.0) THEN
402                     iju=iju+1
403                     indu(iju)=ij
404                  ENDIF
405               ENDDO
406               niju=iju
407c              PRINT*,'niju,nl',niju,nl(l)
408
409c  traitement des mailles
410               DO iju=1,niju
411                  ij=indu(iju)
412                  j=(ij-1)/iip1+1
413                  zu_m=u_m(ij,l)
414                  u_mq(ij,l)=0.
415                  IF(zu_m.gt.0.) THEN
416                     ijq=ij
417                     i=ijq-(j-1)*iip1
418c   accumulation pour les mailles completements advectees
419                     do while(zu_m.gt.masse(ijq,l))
420                        u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
421                        zu_m=zu_m-masse(ijq,l)
422                        i=mod(i-2+iim,iim)+1
423                        ijq=(j-1)*iip1+i
424                     ENDDO
425c   ajout de la maille non completement advectee
426                     u_mq(ij,l)=u_mq(ij,l)+zu_m*
427     &               (q(ijq,l)+0.5*(1.-zu_m/masse(ijq,l))*dxq(ijq,l))
428                  ELSE
429                     ijq=ij+1
430                     i=ijq-(j-1)*iip1
431c   accumulation pour les mailles completements advectees
432                     do while(-zu_m.gt.masse(ijq,l))
433                        u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
434                        zu_m=zu_m+masse(ijq,l)
435                        i=mod(i,iim)+1
436                        ijq=(j-1)*iip1+i
437                     ENDDO
438c   ajout de la maille non completement advectee
439                     u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l)-
440     &               0.5*(1.+zu_m/masse(ijq,l))*dxq(ijq,l))
441                  ENDIF
442               ENDDO
443            ENDIF
444         ENDDO
445      ENDIF  ! n0.gt.0
446
447
448
449c   bouclage en latitude
450
451      DO l=1,llm
452        DO ij=iip1+iip1,ip1jm,iip1
453           u_mq(ij,l)=u_mq(ij-iim,l)
454        ENDDO
455      ENDDO
456
457
458c   calcul des tendances
459
460      DO l=1,llm
461         DO ij=iip2+1,ip1jm
462            new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
463            q(ij,l)=(q(ij,l)*masse(ij,l)+
464     &      u_mq(ij-1,l)-u_mq(ij,l))
465     &      /new_m
466            masse(ij,l)=new_m
467         ENDDO
468c   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
469         DO ij=iip1+iip1,ip1jm,iip1
470            q(ij-iim,l)=q(ij,l)
471            masse(ij-iim,l)=masse(ij,l)
472         ENDDO
473      ENDDO
474
475c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
476c     CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
477
478
479      RETURN
480      END
481      SUBROUTINE vlyqs(q,pente_max,masse,masse_adv_v,qsat)
482c
483c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
484c
485c    ********************************************************************
486c     Shema  d'advection " pseudo amont " .
487c    ********************************************************************
488c     q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
489c     qsat             est   un argument de sortie pour le s-pg ....
490c
491c
492c   --------------------------------------------------------------------
493      IMPLICIT NONE
494c
495#include "dimensions.h"
496#include "paramet.h"
497#include "logic.h"
498#include "comvert.h"
499#include "comconst.h"
500#include "comgeom.h"
501c
502c
503c   Arguments:
504c   ----------
505      REAL masse(ip1jmp1,llm),pente_max
506      REAL masse_adv_v( ip1jm,llm)
507      REAL q(ip1jmp1,llm)
508      REAL qsat(ip1jmp1,llm)
509c
510c      Local
511c   ---------
512c
513      INTEGER i,ij,l
514c
515      REAL airej2,airejjm,airescb(iim),airesch(iim)
516      REAL dyq(ip1jmp1,llm),dyqv(ip1jm)
517      REAL adyqv(ip1jm),dyqmax(ip1jmp1)
518      REAL qbyv(ip1jm,llm)
519
520      REAL qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
521c     REAL newq,oldmasse
522      Logical first,testcpu
523      REAL temps0,temps1,temps2,temps3,temps4,temps5
524      SAVE temps0,temps1,temps2,temps3,temps4,temps5
525      SAVE first,testcpu
526
527      REAL convpn,convps,convmpn,convmps
528      REAL sinlon(iip1),sinlondlon(iip1)
529      REAL coslon(iip1),coslondlon(iip1)
530      SAVE sinlon,coslon,sinlondlon,coslondlon
531      SAVE airej2,airejjm
532c
533c
534      REAL      SSUM
535
536      DATA first,testcpu/.true.,.false./
537      DATA temps0,temps1,temps2,temps3,temps4,temps5/0.,0.,0.,0.,0.,0./
538
539      IF(first) THEN
540         PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
541         first=.false.
542         do i=2,iip1
543            coslon(i)=cos(rlonv(i))
544            sinlon(i)=sin(rlonv(i))
545            coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
546            sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
547         ENDDO
548         coslon(1)=coslon(iip1)
549         coslondlon(1)=coslondlon(iip1)
550         sinlon(1)=sinlon(iip1)
551         sinlondlon(1)=sinlondlon(iip1)
552         airej2 = SSUM( iim, aire(iip2), 1 )
553         airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
554      ENDIF
555
556c
557
558
559      DO l = 1, llm
560c
561c   --------------------------------
562c      CALCUL EN LATITUDE
563c   --------------------------------
564
565c   On commence par calculer la valeur du traceur moyenne sur le premier cercle
566c   de latitude autour du pole (qpns pour le pole nord et qpsn pour
567c    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
568
569      DO i = 1, iim
570      airescb(i) = aire(i+ iip1) * q(i+ iip1,l)
571      airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l)
572      ENDDO
573      qpns   = SSUM( iim,  airescb ,1 ) / airej2
574      qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
575
576c   calcul des pentes aux points v
577
578      DO ij=1,ip1jm
579         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
580         adyqv(ij)=abs(dyqv(ij))
581      ENDDO
582
583c   calcul des pentes aux points scalaires
584
585      DO ij=iip2,ip1jm
586         dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
587         dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
588         dyqmax(ij)=pente_max*dyqmax(ij)
589      ENDDO
590
591c   calcul des pentes aux poles
592
593      DO ij=1,iip1
594         dyq(ij,l)=qpns-q(ij+iip1,l)
595         dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l)-qpsn
596      ENDDO
597
598c   filtrage de la derivee
599      dyn1=0.
600      dys1=0.
601      dyn2=0.
602      dys2=0.
603      DO ij=1,iim
604         dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
605         dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
606         dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
607         dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
608      ENDDO
609      DO ij=1,iip1
610         dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
611         dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
612      ENDDO
613
614c   calcul des pentes limites aux poles
615
616      fn=1.
617      fs=1.
618      DO ij=1,iim
619         IF(pente_max*adyqv(ij).lt.abs(dyq(ij,l))) THEN
620            fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
621         ENDIF
622      IF(pente_max*adyqv(ij+ip1jm-iip1).lt.abs(dyq(ij+ip1jm,l))) THEN
623         fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
624         ENDIF
625      ENDDO
626      DO ij=1,iip1
627         dyq(ij,l)=fn*dyq(ij,l)
628         dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
629      ENDDO
630
631CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
632C  En memoire de dIFferents tests sur la
633C  limitation des pentes aux poles.
634CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
635C     PRINT*,dyq(1)
636C     PRINT*,dyqv(iip1+1)
637C     apn=abs(dyq(1)/dyqv(iip1+1))
638C     PRINT*,dyq(ip1jm+1)
639C     PRINT*,dyqv(ip1jm-iip1+1)
640C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
641C     DO ij=2,iim
642C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
643C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
644C     ENDDO
645C     apn=min(pente_max/apn,1.)
646C     aps=min(pente_max/aps,1.)
647C
648C
649C   cas ou on a un extremum au pole
650C
651C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
652C    &   apn=0.
653C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
654C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
655C    &   aps=0.
656C
657C   limitation des pentes aux poles
658C     DO ij=1,iip1
659C        dyq(ij)=apn*dyq(ij)
660C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
661C     ENDDO
662C
663C   test
664C      DO ij=1,iip1
665C         dyq(iip1+ij)=0.
666C         dyq(ip1jm+ij-iip1)=0.
667C      ENDDO
668C      DO ij=1,ip1jmp1
669C         dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
670C      ENDDO
671C
672C changement 10 07 96
673C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
674C    &   THEN
675C        DO ij=1,iip1
676C           dyqmax(ij)=0.
677C        ENDDO
678C     ELSE
679C        DO ij=1,iip1
680C           dyqmax(ij)=pente_max*abs(dyqv(ij))
681C        ENDDO
682C     ENDIF
683C
684C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
685C    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
686C    &THEN
687C        DO ij=ip1jm+1,ip1jmp1
688C           dyqmax(ij)=0.
689C        ENDDO
690C     ELSE
691C        DO ij=ip1jm+1,ip1jmp1
692C           dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
693C        ENDDO
694C     ENDIF
695C   fin changement 10 07 96
696CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
697
698c   calcul des pentes limitees
699
700      DO ij=iip2,ip1jm
701         IF(dyqv(ij)*dyqv(ij-iip1).gt.0.) THEN
702            dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
703         ELSE
704            dyq(ij,l)=0.
705         ENDIF
706      ENDDO
707
708      ENDDO
709
710      DO l=1,llm
711       DO ij=1,ip1jm
712         IF( masse_adv_v(ij,l).GT.0. ) THEN
713           qbyv(ij,l)= MIN( qsat(ij+iip1,l), q(ij+iip1,l )  +
714     ,      dyq(ij+iip1,l)*0.5*(1.-masse_adv_v(ij,l)/masse(ij+iip1,l)))
715         ELSE
716              qbyv(ij,l)= MIN( qsat(ij,l), q(ij,l) - dyq(ij,l) *
717     ,                   0.5*(1.+masse_adv_v(ij,l)/masse(ij,l)) )
718         ENDIF
719          qbyv(ij,l) = masse_adv_v(ij,l)*qbyv(ij,l)
720       ENDDO
721      ENDDO
722
723
724      DO l=1,llm
725         DO ij=iip2,ip1jm
726            newmasse=masse(ij,l)
727     &      +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
728            q(ij,l)=(q(ij,l)*masse(ij,l)+qbyv(ij,l)-qbyv(ij-iip1,l))
729     &         /newmasse
730            masse(ij,l)=newmasse
731         ENDDO
732c.-. ancienne version
733         convpn=SSUM(iim,qbyv(1,l),1)/apoln
734         convmpn=ssum(iim,masse_adv_v(1,l),1)/apoln
735         DO ij = 1,iip1
736            newmasse=masse(ij,l)+convmpn*aire(ij)
737            q(ij,l)=(q(ij,l)*masse(ij,l)+convpn*aire(ij))/
738     &               newmasse
739            masse(ij,l)=newmasse
740         ENDDO
741         convps  = -SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
742         convmps = -SSUM(iim,masse_adv_v(ip1jm-iim,l),1)/apols
743         DO ij = ip1jm+1,ip1jmp1
744            newmasse=masse(ij,l)+convmps*aire(ij)
745            q(ij,l)=(q(ij,l)*masse(ij,l)+convps*aire(ij))/
746     &               newmasse
747            masse(ij,l)=newmasse
748         ENDDO
749c.-. fin ancienne version
750
751c._. nouvelle version
752c        convpn=SSUM(iim,qbyv(1,l),1)
753c        convmpn=ssum(iim,masse_adv_v(1,l),1)
754c        oldmasse=ssum(iim,masse(1,l),1)
755c        newmasse=oldmasse+convmpn
756c        newq=(q(1,l)*oldmasse+convpn)/newmasse
757c        newmasse=newmasse/apoln
758c        DO ij = 1,iip1
759c           q(ij,l)=newq
760c           masse(ij,l)=newmasse*aire(ij)
761c        ENDDO
762c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
763c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
764c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
765c        newmasse=oldmasse+convmps
766c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
767c        newmasse=newmasse/apols
768c        DO ij = ip1jm+1,ip1jmp1
769c           q(ij,l)=newq
770c           masse(ij,l)=newmasse*aire(ij)
771c        ENDDO
772c._. fin nouvelle version
773      ENDDO
774
775      RETURN
776      END
Note: See TracBrowser for help on using the repository browser.