source: LMDZ6/trunk/libf/dyn3d/vlsplt.F

Last change on this file was 4470, checked in by Laurent Fairhead, 15 months ago

Replaced STOP instructions by calls to abort_gcm for a cleaner exit

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