source: LMDZ.3.3/trunk/libf/dyn3d/vlspltqs.F @ 1611

Last change on this file since 1611 was 78, checked in by lmdz, 25 years ago

Modification de F.Codron pour le schema d'advection de l'eau

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