Ignore:
Timestamp:
Oct 21, 2024, 7:05:31 PM (17 hours ago)
Author:
abarral
Message:

(cont.) Convert fixed-form to free-form sources .F -> .{f,F}90

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/dyn3d/vlsplt.F90

    r5247 r5248  
    1 c
    2 c $Id$
    3 c
    4 
    5       SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq)
    6       USE infotrac, ONLY: nqtot,tracers
    7 c
    8 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    9 c
    10 c    ********************************************************************
    11 c     Shema  d'advection " pseudo amont " .
    12 c    ********************************************************************
    13 c     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
    14 c
    15 c   pente_max facteur de limitation des pentes: 2 en general
    16 c                                               0 pour un schema amont
    17 c   pbaru,pbarv,w flux de masse en u ,v ,w
    18 c   pdt pas de temps
    19 c
    20 c   --------------------------------------------------------------------
    21       IMPLICIT NONE
    22 c
    23       include "dimensions.h"
    24       include "paramet.h"
    25 
    26 c
    27 c   Arguments:
    28 c   ----------
    29       REAL masse(ip1jmp1,llm),pente_max
    30 c      REAL masse(iip1,jjp1,llm),pente_max
    31       REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
    32       REAL q(ip1jmp1,llm,nqtot)
    33 c      REAL q(iip1,jjp1,llm)
    34       REAL w(ip1jmp1,llm),pdt
    35       INTEGER iq ! CRisi
    36 c
    37 c      Local
    38 c   ---------
    39 c
    40       INTEGER ij,l
    41 c
    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
     1!
     2! $Id$
     3!
     4
     5SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq)
     6  USE infotrac, ONLY: nqtot,tracers
     7  !
     8  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
     9  !
     10  !    ********************************************************************
     11  ! Shema  d'advection " pseudo amont " .
     12  !    ********************************************************************
     13  ! q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
     14  !
     15  !   pente_max facteur de limitation des pentes: 2 en general
     16  !                                           0 pour un schema amont
     17  !   pbaru,pbarv,w flux de masse en u ,v ,w
     18  !   pdt pas de temps
     19  !
     20  !   --------------------------------------------------------------------
     21  IMPLICIT NONE
     22  !
     23  include "dimensions.h"
     24  include "paramet.h"
     25
     26  !
     27  !   Arguments:
     28  !   ----------
     29  REAL :: masse(ip1jmp1,llm),pente_max
     30   ! REAL masse(iip1,jjp1,llm),pente_max
     31  REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
     32  REAL :: q(ip1jmp1,llm,nqtot)
     33   ! REAL q(iip1,jjp1,llm)
     34  REAL :: w(ip1jmp1,llm),pdt
     35  INTEGER :: iq ! CRisi
     36  !
     37  !  Local
     38  !   ---------
     39  !
     40  INTEGER :: ij,l
     41  !
     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
     79  !print*,'Entree vlx1'
     80  ! call minmaxq(zq,qmin,qmax,'avant vlx     ')
     81  call vlx(zq,pente_max,zm,mu,iq)
     82  !print*,'Sortie vlx1'
     83  ! call minmaxq(zq,qmin,qmax,'apres vlx1    ')
     84
     85  ! print*,'Entree vly1'
     86
     87  call vly(zq,pente_max,zm,mv,iq)
     88  ! call minmaxq(zq,qmin,qmax,'apres vly1     ')
     89  !print*,'Sortie vly1'
     90  call vlz(zq,pente_max,zm,mw,iq)
     91  ! call minmaxq(zq,qmin,qmax,'apres vlz     ')
     92
     93
     94  call vly(zq,pente_max,zm,mv,iq)
     95  ! call minmaxq(zq,qmin,qmax,'apres vly     ')
     96
     97
     98  call vlx(zq,pente_max,zm,mu,iq)
     99  ! 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)
    65116      ENDDO
    66 
     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
     124END SUBROUTINE vlsplt
     125RECURSIVE 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
     129  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
     130  !
     131  !    ********************************************************************
     132  ! Shema  d'advection " pseudo amont " .
     133  !    ********************************************************************
     134  ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
     135  !
     136  !
     137  !   --------------------------------------------------------------------
     138  IMPLICIT NONE
     139  !
     140  include "dimensions.h"
     141  include "paramet.h"
     142  include "iniprint.h"
     143  !
     144  !
     145  !   Arguments:
     146  !   ----------
     147  REAL :: masse(ip1jmp1,llm,nqtot),pente_max
     148  REAL :: u_m( ip1jmp1,llm )
     149  REAL :: q(ip1jmp1,llm,nqtot)
     150  INTEGER :: iq ! CRisi
     151  !
     152  !  Local
     153  !   ---------
     154  !
     155  INTEGER :: ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
     156  INTEGER :: n0,iadvplus(ip1jmp1,llm),nl(llm)
     157  !
     158  REAL :: new_m,zu_m,zdum(ip1jmp1,llm)
     159   ! 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
     173  !   calcul de la pente a droite et a gauche de la maille
     174
     175
     176  IF (pente_max.gt.-1.e-5) THEN
     177    ! IF (pente_max.gt.10) THEN
     178
     179  !   calcul des pentes avec limitation, Van Leer scheme I:
     180  !   -----------------------------------------------------
     181
     182  !   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)
     189           ! sigu(ij)=sigu(ij-iim)
     190        ENDDO
     191
     192        DO ij=iip2,ip1jm
     193           adxqu(ij)=abs(dxqu(ij))
     194        ENDDO
     195
     196  !   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))
     201  ! limitation subtile
     202  !    ,      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
     219  !   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
     229  !print*,'Ok calcul des pentes'
     230
     231  ELSE ! (pente_max.lt.-1.e-5)
     232
     233  !   Pentes produits:
     234  !   ----------------
     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
     250  !   extremum local
     251              dxq(ij,l)=0.
     252           ENDIF
     253        ENDDO
     254
     255     ENDDO
     256
     257  ENDIF ! (pente_max.lt.-1.e-5)
     258
     259  !   bouclage de la pente en iip1:
     260  !   -----------------------------
     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
     272  ! print*,'Bouclage en iip1'
     273
     274  !   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
     292  !   on cumule le flux correspondant a toutes les mailles dont la masse
     293  !   au travers de la paroi pENDant le pas de temps.
     294  !print*,'Cumule ....'
     295
     296  DO l=1,llm
     297   DO ij=iip2,ip1jm-1
     298  ! 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
     311  ! go to 9999
     312  !   detection des points ou on advecte plus que la masse de la
     313  !   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
     322  !print*,'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
     328  ! print*,'Ok test 2'
     329
     330
     331  !   traitement special pour le cas ou on advecte en longitude plus que le
     332  !   contenu de la maille.
     333  !   cette partie est mal vectorisee.
     334
     335  !  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
     354  !   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
     362           ! PRINT*,'niju,nl',niju,nl(l)
     363
     364  !  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
     373  !   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
     381  !   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
     388  !   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
     396  !   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
     404  !9999    continue
     405
     406
     407  !   bouclage en latitude
     408  !print*,'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
     443  !   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
     454  !   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
     476  ! CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
     477  ! CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
     478
     479
     480  RETURN
     481END SUBROUTINE vlx
     482RECURSIVE 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
     485  !
     486  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
     487  !
     488  !    ********************************************************************
     489  ! Shema  d'advection " pseudo amont " .
     490  !    ********************************************************************
     491  ! q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
     492  ! dq         sont des arguments de sortie pour le s-pg ....
     493  !
     494  !
     495  !   --------------------------------------------------------------------
     496  USE comconst_mod, ONLY: pi
     497  IMPLICIT NONE
     498  !
     499  include "dimensions.h"
     500  include "paramet.h"
     501  include "comgeom.h"
     502  !
     503  !
     504  !   Arguments:
     505  !   ----------
     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
     510  !
     511  !  Local
     512  !   ---------
     513  !
     514  INTEGER :: i,ij,l
     515  !
     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
     522  ! REAL appn apps
     523  ! 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
     537  !
     538  !
     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
     562  !
     563  !PRINT*,'CALCUL EN LATITUDE'
     564
     565  DO l = 1, llm
     566  !
     567  !   --------------------------------
     568  !  CALCUL EN LATITUDE
     569  !   --------------------------------
     570
     571  !   On commence par calculer la valeur du traceur moyenne sur le premier cercle
     572  !   de latitude autour du pole (qpns pour le pole nord et qpsn pour
     573  !    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
     582  !   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
     589  !   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
     597  !   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
     604  !   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
     620  !   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
     643  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     644  !  En memoire de dIFferents tests sur la
     645  !  limitation des pentes aux poles.
     646  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     647  ! PRINT*,dyq(1)
     648  ! PRINT*,dyqv(iip1+1)
     649  ! appn=abs(dyq(1)/dyqv(iip1+1))
     650  ! PRINT*,dyq(ip1jm+1)
     651  ! PRINT*,dyqv(ip1jm-iip1+1)
     652  ! apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     653  ! DO ij=2,iim
     654  !    appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
     655  !    apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
     656  ! ENDDO
     657  ! appn=min(pente_max/appn,1.)
     658  ! apps=min(pente_max/apps,1.)
     659  !
     660  !
     661  !   cas ou on a un extremum au pole
     662  !
     663  ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
     664  !    &   appn=0.
     665  ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
     666  !    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
     667  !    &   apps=0.
     668  !
     669  !   limitation des pentes aux poles
     670  ! DO ij=1,iip1
     671  !    dyq(ij)=appn*dyq(ij)
     672  !    dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
     673  ! ENDDO
     674  !
     675  !   test
     676  !  DO ij=1,iip1
     677  !     dyq(iip1+ij)=0.
     678  !     dyq(ip1jm+ij-iip1)=0.
     679  !  ENDDO
     680  !  DO ij=1,ip1jmp1
     681  !     dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
     682  !  ENDDO
     683  !
     684  ! changement 10 07 96
     685  ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
     686  !    &   THEN
     687  !    DO ij=1,iip1
     688  !       dyqmax(ij)=0.
     689  !    ENDDO
     690  ! ELSE
     691  !    DO ij=1,iip1
     692  !       dyqmax(ij)=pente_max*abs(dyqv(ij))
     693  !    ENDDO
     694  ! ENDIF
     695  !
     696  ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
     697  !    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
     698  !    &THEN
     699  !    DO ij=ip1jm+1,ip1jmp1
     700  !       dyqmax(ij)=0.
     701  !    ENDDO
     702  ! ELSE
     703  !    DO ij=ip1jm+1,ip1jmp1
     704  !       dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
     705  !    ENDDO
     706  ! ENDIF
     707  !   fin changement 10 07 96
     708  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     709
     710  !   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
    67745      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 
    79 cprint*,'Entree vlx1'
    80 c       call minmaxq(zq,qmin,qmax,'avant vlx     ')
    81       call vlx(zq,pente_max,zm,mu,iq)
    82 cprint*,'Sortie vlx1'
    83 c       call minmaxq(zq,qmin,qmax,'apres vlx1    ')
    84 
    85 c print*,'Entree vly1'
    86 
    87       call vly(zq,pente_max,zm,mv,iq)
    88 c       call minmaxq(zq,qmin,qmax,'apres vly1     ')
    89 cprint*,'Sortie vly1'
    90       call vlz(zq,pente_max,zm,mw,iq)
    91 c       call minmaxq(zq,qmin,qmax,'apres vlz     ')
    92 
    93 
    94       call vly(zq,pente_max,zm,mv,iq)
    95 c       call minmaxq(zq,qmin,qmax,'apres vly     ')
    96 
    97 
    98       call vlx(zq,pente_max,zm,mu,iq)
    99 c       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
     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
    121757      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 
    129 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    130 c
    131 c    ********************************************************************
    132 c     Shema  d'advection " pseudo amont " .
    133 c    ********************************************************************
    134 c     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
    135 c
    136 c
    137 c   --------------------------------------------------------------------
    138       IMPLICIT NONE
    139 c
    140       include "dimensions.h"
    141       include "paramet.h"
    142       include "iniprint.h"
    143 c
    144 c
    145 c   Arguments:
    146 c   ----------
    147       REAL masse(ip1jmp1,llm,nqtot),pente_max
    148       REAL u_m( ip1jmp1,llm )
    149       REAL q(ip1jmp1,llm,nqtot)
    150       INTEGER iq ! CRisi
    151 c
    152 c      Local
    153 c   ---------
    154 c
    155       INTEGER ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
    156       INTEGER n0,iadvplus(ip1jmp1,llm),nl(llm)
    157 c
    158       REAL new_m,zu_m,zdum(ip1jmp1,llm)
    159 c      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 
    173 c   calcul de la pente a droite et a gauche de la maille
    174 
    175 
    176       IF (pente_max.gt.-1.e-5) THEN
    177 c       IF (pente_max.gt.10) THEN
    178 
    179 c   calcul des pentes avec limitation, Van Leer scheme I:
    180 c   -----------------------------------------------------
    181 
    182 c   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)
    189 c              sigu(ij)=sigu(ij-iim)
    190             ENDDO
    191 
    192             DO ij=iip2,ip1jm
    193                adxqu(ij)=abs(dxqu(ij))
    194             ENDDO
    195 
    196 c   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))
    201 c limitation subtile
    202 c    ,      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
     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
     774  !.-. ancienne version
     775     ! convpn=SSUM(iim,qbyv(1,l),1)/apoln
     776     ! 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
     790     ! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
     791     ! 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
     805  !.-. fin ancienne version
     806
     807  !._. nouvelle version
     808     ! convpn=SSUM(iim,qbyv(1,l),1)
     809     ! convmpn=ssum(iim,masse_adv_v(1,l),1)
     810     ! oldmasse=ssum(iim,masse(1,l),1)
     811     ! newmasse=oldmasse+convmpn
     812     ! newq=(q(1,l)*oldmasse+convpn)/newmasse
     813     ! newmasse=newmasse/apoln
     814     ! DO ij = 1,iip1
     815     !    q(ij,l)=newq
     816     !    masse(ij,l,iq)=newmasse*aire(ij)
     817     ! ENDDO
     818     ! convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
     819     ! convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
     820     ! oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
     821     ! newmasse=oldmasse+convmps
     822     ! newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
     823     ! newmasse=newmasse/apols
     824     ! DO ij = ip1jm+1,ip1jmp1
     825     !    q(ij,l)=newq
     826     !    masse(ij,l,iq)=newmasse*aire(ij)
     827     ! ENDDO
     828  !._. 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
     844END SUBROUTINE vly
     845RECURSIVE 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
     848  !
     849  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
     850  !
     851  !    ********************************************************************
     852  ! Shema  d'advection " pseudo amont " .
     853  !    ********************************************************************
     854  !    q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
     855  ! dq         sont des arguments de sortie pour le s-pg ....
     856  !
     857  !
     858  !   --------------------------------------------------------------------
     859  IMPLICIT NONE
     860  !
     861  include "dimensions.h"
     862  include "paramet.h"
     863  !
     864  !
     865  !   Arguments:
     866  !   ----------
     867  REAL :: masse(ip1jmp1,llm,nqtot),pente_max
     868  REAL :: q(ip1jmp1,llm,nqtot)
     869  REAL :: w(ip1jmp1,llm+1)
     870  INTEGER :: iq
     871  !
     872  !  Local
     873  !   ---------
     874  !
     875  INTEGER :: ij,l
     876  !
     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
     896  !    On oriente tout dans le sens de la pression c'est a dire dans le
     897  !    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
    212915#ifdef CRAY
    213                dxq(ij,l)=
    214      ,         cvmgp(dxqu(ij-1)+dxqu(ij),0.,dxqu(ij-1)*dxqu(ij))
     916        dzq(ij,l)=0.5* &
     917              cvmgp(dzqw(ij,l)+dzqw(ij,l+1),0.,dzqw(ij,l)*dzqw(ij,l+1))
    215918#else
    216                IF(dxqu(ij-1)*dxqu(ij).gt.0) THEN
    217                   dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
    218                ELSE
    219 c   extremum local
    220                   dxq(ij,l)=0.
    221                ENDIF
     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
    222924#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
    229 cprint*,'Ok calcul des pentes'
    230 
    231       ELSE ! (pente_max.lt.-1.e-5)
    232 
    233 c   Pentes produits:
    234 c   ----------------
    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
    250 c   extremum local
    251                   dxq(ij,l)=0.
    252                ENDIF
    253             ENDDO
    254 
    255          ENDDO
    256 
    257       ENDIF ! (pente_max.lt.-1.e-5)
    258 
    259 c   bouclage de la pente en iip1:
    260 c   -----------------------------
    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 
    272 c print*,'Bouclage en iip1'
    273 
    274 c   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
    292 c   on cumule le flux correspondant a toutes les mailles dont la masse
    293 c   au travers de la paroi pENDant le pas de temps.
    294 cprint*,'Cumule ....'
    295 
    296       DO l=1,llm
    297        DO ij=iip2,ip1jm-1
    298 c       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
     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
    309940#endif
    310 
    311 c       go to 9999
    312 c   detection des points ou on advecte plus que la masse de la
    313 c   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
    322 cprint*,'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
    328 c print*,'Ok test 2'
    329 
    330 
    331 c   traitement special pour le cas ou on advecte en longitude plus que le
    332 c   contenu de la maille.
    333 c   cette partie est mal vectorisee.
    334 
    335 c  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
    354 c   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
    362 c              PRINT*,'niju,nl',niju,nl(l)
    363 
    364 c  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
    373 c   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
    381 c   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
    388 c   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
    396 c   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
    404 c9999    continue
    405 
    406 
    407 c   bouclage en latitude
    408 cprint*,'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
     941  ! ---------------------------------------------------------------
     942  !   .... calcul des termes d'advection verticale  .......
     943  ! ---------------------------------------------------------------
     944
     945  ! 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
    435982      enddo
    436       do ifils=1,tracers(iq)%nqChildren
    437         iq2=tracers(iq)%iqDescen(ifils)
    438         call vlx(Ratio,pente_max,masseq,u_mq,iq2)
     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)
    4391007      enddo
    440 ! end CRisi
    441 
    442 
    443 c   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
    454 c   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 
    476 c     CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
    477 c     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
    485 c
    486 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    487 c
    488 c    ********************************************************************
    489 c     Shema  d'advection " pseudo amont " .
    490 c    ********************************************************************
    491 c     q,masse_adv_v,w sont des arguments d'entree  pour le s-pg ....
    492 c     dq               sont des arguments de sortie pour le s-pg ....
    493 c
    494 c
    495 c   --------------------------------------------------------------------
    496       USE comconst_mod, ONLY: pi
    497       IMPLICIT NONE
    498 c
    499       include "dimensions.h"
    500       include "paramet.h"
    501       include "comgeom.h"
    502 c
    503 c
    504 c   Arguments:
    505 c   ----------
    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
    510 c
    511 c      Local
    512 c   ---------
    513 c
    514       INTEGER i,ij,l
    515 c
    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
    522 c     REAL appn apps
    523 c     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 
    537 c
    538 c
    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 
    562 c
    563 cPRINT*,'CALCUL EN LATITUDE'
    564 
    565       DO l = 1, llm
    566 c
    567 c   --------------------------------
    568 c      CALCUL EN LATITUDE
    569 c   --------------------------------
    570 
    571 c   On commence par calculer la valeur du traceur moyenne sur le premier cercle
    572 c   de latitude autour du pole (qpns pour le pole nord et qpsn pour
    573 c    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 
    582 c   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 
    589 c   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 
    597 c   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 
    604 c   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 
    620 c   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
    637 8888    continue
    638       DO ij=1,iip1
    639          dyq(ij,l)=0.
    640          dyq(ip1jm+ij,l)=0.
    641       ENDDO
    642 
    643 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    644 C  En memoire de dIFferents tests sur la
    645 C  limitation des pentes aux poles.
    646 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    647 C     PRINT*,dyq(1)
    648 C     PRINT*,dyqv(iip1+1)
    649 C     appn=abs(dyq(1)/dyqv(iip1+1))
    650 C     PRINT*,dyq(ip1jm+1)
    651 C     PRINT*,dyqv(ip1jm-iip1+1)
    652 C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    653 C     DO ij=2,iim
    654 C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
    655 C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    656 C     ENDDO
    657 C     appn=min(pente_max/appn,1.)
    658 C     apps=min(pente_max/apps,1.)
    659 C
    660 C
    661 C   cas ou on a un extremum au pole
    662 C
    663 C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    664 C    &   appn=0.
    665 C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    666 C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    667 C    &   apps=0.
    668 C
    669 C   limitation des pentes aux poles
    670 C     DO ij=1,iip1
    671 C        dyq(ij)=appn*dyq(ij)
    672 C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    673 C     ENDDO
    674 C
    675 C   test
    676 C      DO ij=1,iip1
    677 C         dyq(iip1+ij)=0.
    678 C         dyq(ip1jm+ij-iip1)=0.
    679 C      ENDDO
    680 C      DO ij=1,ip1jmp1
    681 C         dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
    682 C      ENDDO
    683 C
    684 C changement 10 07 96
    685 C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    686 C    &   THEN
    687 C        DO ij=1,iip1
    688 C           dyqmax(ij)=0.
    689 C        ENDDO
    690 C     ELSE
    691 C        DO ij=1,iip1
    692 C           dyqmax(ij)=pente_max*abs(dyqv(ij))
    693 C        ENDDO
    694 C     ENDIF
    695 C
    696 C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    697 C    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    698 C    &THEN
    699 C        DO ij=ip1jm+1,ip1jmp1
    700 C           dyqmax(ij)=0.
    701 C        ENDDO
    702 C     ELSE
    703 C        DO ij=ip1jm+1,ip1jmp1
    704 C           dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
    705 C        ENDDO
    706 C     ENDIF
    707 C   fin changement 10 07 96
    708 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    709 
    710 c   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
    774 c.-. ancienne version
    775 c        convpn=SSUM(iim,qbyv(1,l),1)/apoln
    776 c        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 
    790 c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)/apols
    791 c        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 
    805 c.-. fin ancienne version
    806 
    807 c._. nouvelle version
    808 c        convpn=SSUM(iim,qbyv(1,l),1)
    809 c        convmpn=ssum(iim,masse_adv_v(1,l),1)
    810 c        oldmasse=ssum(iim,masse(1,l),1)
    811 c        newmasse=oldmasse+convmpn
    812 c        newq=(q(1,l)*oldmasse+convpn)/newmasse
    813 c        newmasse=newmasse/apoln
    814 c        DO ij = 1,iip1
    815 c           q(ij,l)=newq
    816 c           masse(ij,l,iq)=newmasse*aire(ij)
    817 c        ENDDO
    818 c        convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    819 c        convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
    820 c        oldmasse=ssum(iim,masse(ip1jm-iim,l),1)
    821 c        newmasse=oldmasse+convmps
    822 c        newq=(q(ip1jmp1,l)*oldmasse+convps)/newmasse
    823 c        newmasse=newmasse/apols
    824 c        DO ij = ip1jm+1,ip1jmp1
    825 c           q(ij,l)=newq
    826 c           masse(ij,l,iq)=newmasse*aire(ij)
    827 c        ENDDO
    828 c._. 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
    848 c
    849 c     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    850 c
    851 c    ********************************************************************
    852 c     Shema  d'advection " pseudo amont " .
    853 c    ********************************************************************
    854 c    q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
    855 c     dq               sont des arguments de sortie pour le s-pg ....
    856 c
    857 c
    858 c   --------------------------------------------------------------------
    859       IMPLICIT NONE
    860 c
    861       include "dimensions.h"
    862       include "paramet.h"
    863 c
    864 c
    865 c   Arguments:
    866 c   ----------
    867       REAL masse(ip1jmp1,llm,nqtot),pente_max
    868       REAL q(ip1jmp1,llm,nqtot)
    869       REAL w(ip1jmp1,llm+1)
    870       INTEGER iq
    871 c
    872 c      Local
    873 c   ---------
    874 c
    875       INTEGER ij,l
    876 c
    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 
    896 c    On oriente tout dans le sens de la pression c'est a dire dans le
    897 c    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
    941 c ---------------------------------------------------------------
    942 c   .... calcul des termes d'advection verticale  .......
    943 c ---------------------------------------------------------------
    944 
    945 c 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
    1014 c      SUBROUTINE minmaxq(zq,qmin,qmax,comment)
    1015 c
    1016 c#include "dimensions.h"
    1017 c#include "paramet.h"
    1018 
    1019 c      CHARACTER*(*) comment
    1020 c      real qmin,qmax
    1021 c      real zq(ip1jmp1,llm)
    1022 
    1023 c      INTEGER jadrs(ip1jmp1), jbad, k, i
    1024 
    1025 
    1026 c      DO k = 1, llm
    1027 c         jbad = 0
    1028 c         DO i = 1, ip1jmp1
    1029 c         IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
    1030 c            jbad = jbad + 1
    1031 c            jadrs(jbad) = i
    1032 c         ENDIF
    1033 c         ENDDO
    1034 c         IF (jbad.GT.0) THEN
    1035 c         PRINT*, comment
    1036 c         DO i = 1, jbad
    1037 cc            PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
    1038 c         ENDDO
    1039 c         ENDIF
    1040 c      ENDDO
    1041 
    1042 c      return
    1043 c      end
    1044       subroutine minmaxq(zq,qmin,qmax,comment)
     1008    enddo
     1009  enddo
     1010  ! !write(*,*) 'vlsplt 1032'
     1011
     1012  RETURN
     1013END SUBROUTINE vlz
     1014 ! SUBROUTINE minmaxq(zq,qmin,qmax,comment)
     1015!
     1016!#include "dimensions.h"
     1017!#include "paramet.h"
     1018
     1019!  CHARACTER*(*) comment
     1020!  real qmin,qmax
     1021!  real zq(ip1jmp1,llm)
     1022
     1023!  INTEGER jadrs(ip1jmp1), jbad, k, i
     1024
     1025
     1026!  DO k = 1, llm
     1027!     jbad = 0
     1028!     DO i = 1, ip1jmp1
     1029!     IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN
     1030!        jbad = jbad + 1
     1031!        jadrs(jbad) = i
     1032!     ENDIF
     1033!     ENDDO
     1034!     IF (jbad.GT.0) THEN
     1035!     PRINT*, comment
     1036!     DO i = 1, jbad
     1037!c            PRINT*, "i,k,zq=", jadrs(i),k,zq(jadrs(i),k)
     1038!     ENDDO
     1039!     ENDIF
     1040!  ENDDO
     1041
     1042!  return
     1043!  end
     1044subroutine minmaxq(zq,qmin,qmax,comment)
    10451045
    10461046#include "dimensions.h"
    10471047#include "paramet.h"
    10481048
    1049       character*20 comment
    1050       real qmin,qmax
    1051       real zq(ip1jmp1,llm)
    1052       real zzq(iip1,jjp1,llm)
     1049  character(len=20) :: comment
     1050  real :: qmin,qmax
     1051  real :: zq(ip1jmp1,llm)
     1052  real :: zzq(iip1,jjp1,llm)
    10531053
    10541054#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)
    1077 c    s     write(*,9999) comment,
    1078      s     write(*,*) comment,
    1079      s     imin,jmin,lmin,zqmin,zzq(imin,jmin,lmin)
    1080        if(zqmax.gt.qmax)
    1081 c    s     write(*,9999) comment,
    1082      s     write(*,*) comment,
    1083      s     imax,jmax,lmax,zqmax,zzq(imax,jmax,lmax)
     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) &
     1077  ! s     write(*,9999) comment,
     1078         write(*,*) comment, &
     1079         imin,jmin,lmin,zqmin,zzq(imin,jmin,lmin)
     1080   if(zqmax.gt.qmax) &
     1081  ! s     write(*,9999) comment,
     1082         write(*,*) comment, &
     1083         imax,jmax,lmax,zqmax,zzq(imax,jmax,lmax)
    10841084
    10851085#endif
    1086       return
    1087 c9999  format(a20,'  q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
    1088       end
    1089 
    1090 
    1091 
     1086  return
     1087  !9999  format(a20,'  q(',i3,',',i2,',',i2,')=',e12.5,e12.5)
     1088end subroutine minmaxq
     1089
     1090
     1091
Note: See TracChangeset for help on using the changeset viewer.