source: LMDZ6/branches/DYNAMICO-conv/libf/dyn3d/vlsplt.F @ 4976

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

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

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