Ignore:
Timestamp:
Jul 24, 2024, 6:46:45 PM (2 months ago)
Author:
abarral
Message:

enforce PRIVATE by default in several modules, expose PUBLIC as needed
move eigen.f90 to obsolete/
(lint) aslong the way

Location:
LMDZ6/branches/Amaury_dev/libf/dyn3d
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/advtrac.f90

    r5118 r5119  
    1717  USE lmdz_libmath, ONLY: minmax
    1818  USE lmdz_iniprint, ONLY: lunout, prt_level
     19  USE lmdz_ssum_scopy, ONLY: scopy
    1920
    2021  IMPLICIT NONE

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/caladvtrac.F90

    r5117 r5119  
    1111  USE comconst_mod, ONLY: dtvr
    1212  USE lmdz_filtreg, ONLY: filtreg
     13  USE lmdz_ssum_scopy, ONLY: scopy
    1314
    1415  IMPLICIT NONE

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/fluxstokenc.F90

    r5118 r5119  
    77  USE IOIPSL
    88  USE lmdz_iniprint, ONLY: lunout, prt_level
     9  USE lmdz_ssum_scopy, ONLY: scopy
    910  !
    1011  ! Auteur :  F. Hourdin

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/groupe.F90

    r5117 r5119  
    44
    55  USE comconst_mod, ONLY: ngroup
     6  USE lmdz_ssum_scopy, ONLY: scopy
    67
    78  IMPLICIT NONE

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/integrd.F90

    r5118 r5119  
    1212  USE temps_mod, ONLY: dt
    1313  USE lmdz_iniprint, ONLY: lunout, prt_level
     14  USE lmdz_ssum_scopy, ONLY: scopy
    1415
    1516  IMPLICIT NONE

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/leapfrog.F90

    r5118 r5119  
    2727  USE lmdz_description, ONLY: descript
    2828  USE lmdz_iniprint, ONLY: lunout, prt_level
     29  USE lmdz_ssum_scopy, ONLY: scopy
    2930
    3031  IMPLICIT NONE

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/vlsplt.F90

    r5118 r5119  
    33!
    44
    5 SUBROUTINE vlsplt(q,pente_max,masse,w,pbaru,pbarv,pdt,iq)
    6   USE infotrac, ONLY: nqtot,tracers
     5SUBROUTINE vlsplt(q, pente_max, masse, w, pbaru, pbarv, pdt, iq)
     6  USE infotrac, ONLY: nqtot, tracers
     7  USE lmdz_ssum_scopy, ONLY: scopy
    78  !
    89  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    2728  !   Arguments:
    2829  !   ----------
    29   REAL :: masse(ip1jmp1,llm),pente_max
    30   REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
    31   REAL :: q(ip1jmp1,llm,nqtot)
    32   REAL :: w(ip1jmp1,llm),pdt
     30  REAL :: masse(ip1jmp1, llm), pente_max
     31  REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
     32  REAL :: q(ip1jmp1, llm, nqtot)
     33  REAL :: w(ip1jmp1, llm), pdt
    3334  INTEGER :: iq ! CRisi
    3435  !
    … …  
    3637  !   ---------
    3738  !
    38   INTEGER :: ij,l
    39   !
    40   REAL :: zm(ip1jmp1,llm,nqtot)
    41   REAL :: mu(ip1jmp1,llm)
    42   REAL :: mv(ip1jm,llm)
    43   REAL :: mw(ip1jmp1,llm+1)
    44   REAL :: zq(ip1jmp1,llm,nqtot)
     39  INTEGER :: ij, l
     40  !
     41  REAL :: zm(ip1jmp1, llm, nqtot)
     42  REAL :: mu(ip1jmp1, llm)
     43  REAL :: mv(ip1jm, llm)
     44  REAL :: mw(ip1jmp1, llm + 1)
     45  REAL :: zq(ip1jmp1, llm, nqtot)
    4546  REAL :: zzpbar, zzw
    46   INTEGER :: ifils,iq2 ! CRisi
    47 
    48   REAL :: qmin,qmax
    49   DATA qmin,qmax/0.,1.e33/
    50 
    51     zzpbar = 0.5 * pdt
    52     zzw    = pdt
    53   DO l=1,llm
    54     DO ij = iip2,ip1jm
    55         mu(ij,l)=pbaru(ij,l) * zzpbar
    56      ENDDO
    57      DO ij=1,ip1jm
    58         mv(ij,l)=pbarv(ij,l) * zzpbar
    59      ENDDO
    60      DO ij=1,ip1jmp1
    61         mw(ij,l)=w(ij,l) * zzw
    62      ENDDO
    63   ENDDO
    64 
    65   DO ij=1,ip1jmp1
    66      mw(ij,llm+1)=0.
    67   ENDDO
    68 
    69   CALL SCOPY(ijp1llm,q(1,1,iq),1,zq(1,1,iq),1)
    70   CALL SCOPY(ijp1llm,masse,1,zm(1,1,iq),1)
    71 
    72   do ifils=1,tracers(iq)%nqDescen
    73     iq2=tracers(iq)%iqDescen(ifils)
    74     CALL SCOPY(ijp1llm,q(1,1,iq2),1,zq(1,1,iq2),1)
    75   enddo
    76 
    77   CALL vlx(zq,pente_max,zm,mu,iq)
    78   CALL vly(zq,pente_max,zm,mv,iq)
    79   CALL vlz(zq,pente_max,zm,mw,iq)
    80   CALL vly(zq,pente_max,zm,mv,iq)
    81   CALL vlx(zq,pente_max,zm,mu,iq)
    82 
    83   DO l=1,llm
    84      DO ij=1,ip1jmp1
    85        q(ij,l,iq)=zq(ij,l,iq)
    86      ENDDO
    87      DO ij=1,ip1jm+1,iip1
    88         q(ij+iim,l,iq)=q(ij,l,iq)
    89      ENDDO
     47  INTEGER :: ifils, iq2 ! CRisi
     48
     49  REAL :: qmin, qmax
     50  DATA qmin, qmax/0., 1.e33/
     51
     52  zzpbar = 0.5 * pdt
     53  zzw = pdt
     54  DO l = 1, llm
     55    DO ij = iip2, ip1jm
     56      mu(ij, l) = pbaru(ij, l) * zzpbar
     57    ENDDO
     58    DO ij = 1, ip1jm
     59      mv(ij, l) = pbarv(ij, l) * zzpbar
     60    ENDDO
     61    DO ij = 1, ip1jmp1
     62      mw(ij, l) = w(ij, l) * zzw
     63    ENDDO
     64  ENDDO
     65
     66  DO ij = 1, ip1jmp1
     67    mw(ij, llm + 1) = 0.
     68  ENDDO
     69
     70  CALL SCOPY(ijp1llm, q(1, 1, iq), 1, zq(1, 1, iq), 1)
     71  CALL SCOPY(ijp1llm, masse, 1, zm(1, 1, iq), 1)
     72
     73  do ifils = 1, tracers(iq)%nqDescen
     74    iq2 = tracers(iq)%iqDescen(ifils)
     75    CALL SCOPY(ijp1llm, q(1, 1, iq2), 1, zq(1, 1, iq2), 1)
     76  enddo
     77
     78  CALL vlx(zq, pente_max, zm, mu, iq)
     79  CALL vly(zq, pente_max, zm, mv, iq)
     80  CALL vlz(zq, pente_max, zm, mw, iq)
     81  CALL vly(zq, pente_max, zm, mv, iq)
     82  CALL vlx(zq, pente_max, zm, mu, iq)
     83
     84  DO l = 1, llm
     85    DO ij = 1, ip1jmp1
     86      q(ij, l, iq) = zq(ij, l, iq)
     87    ENDDO
     88    DO ij = 1, ip1jm + 1, iip1
     89      q(ij + iim, l, iq) = q(ij, l, iq)
     90    ENDDO
    9091  ENDDO
    9192  ! CRisi: aussi pour les fils
    92   do ifils=1,tracers(iq)%nqDescen
    93     iq2=tracers(iq)%iqDescen(ifils)
    94     DO l=1,llm
    95       DO ij=1,ip1jmp1
    96         q(ij,l,iq2)=zq(ij,l,iq2)
    97       ENDDO
    98       DO ij=1,ip1jm+1,iip1
    99         q(ij+iim,l,iq2)=q(ij,l,iq2)
    100       ENDDO
    101     ENDDO
    102   enddo
    103 
     93  do ifils = 1, tracers(iq)%nqDescen
     94    iq2 = tracers(iq)%iqDescen(ifils)
     95    DO l = 1, llm
     96      DO ij = 1, ip1jmp1
     97        q(ij, l, iq2) = zq(ij, l, iq2)
     98      ENDDO
     99      DO ij = 1, ip1jm + 1, iip1
     100        q(ij + iim, l, iq2) = q(ij, l, iq2)
     101      ENDDO
     102    ENDDO
     103  enddo
    104104
    105105END SUBROUTINE vlsplt
    106 RECURSIVE SUBROUTINE vlx(q,pente_max,masse,u_m,iq)
    107   USE infotrac, ONLY: nqtot,tracers, & ! CRisi
    108         min_qParent,min_qMass,min_ratio ! MVals et CRisi
     106RECURSIVE SUBROUTINE vlx(q, pente_max, masse, u_m, iq)
     107  USE infotrac, ONLY: nqtot, tracers, & ! CRisi
     108          min_qParent, min_qMass, min_ratio ! MVals et CRisi
    109109  USE lmdz_iniprint, ONLY: lunout, prt_level
    110110
    … …  
    126126  !   Arguments:
    127127  !   ----------
    128   REAL :: masse(ip1jmp1,llm,nqtot),pente_max
    129   REAL :: u_m( ip1jmp1,llm )
    130   REAL :: q(ip1jmp1,llm,nqtot)
     128  REAL :: masse(ip1jmp1, llm, nqtot), pente_max
     129  REAL :: u_m(ip1jmp1, llm)
     130  REAL :: q(ip1jmp1, llm, nqtot)
    131131  INTEGER :: iq ! CRisi
    132132  !
    … …  
    134134  !   ---------
    135135  !
    136   INTEGER :: ij,l,j,i,iju,ijq,indu(ip1jmp1),niju
    137   INTEGER :: n0,iadvplus(ip1jmp1,llm),nl(llm)
    138   !
    139   REAL :: new_m,zu_m,zdum(ip1jmp1,llm)
    140   REAL :: dxq(ip1jmp1,llm),dxqu(ip1jmp1)
     136  INTEGER :: ij, l, j, i, iju, ijq, indu(ip1jmp1), niju
     137  INTEGER :: n0, iadvplus(ip1jmp1, llm), nl(llm)
     138  !
     139  REAL :: new_m, zu_m, zdum(ip1jmp1, llm)
     140  REAL :: dxq(ip1jmp1, llm), dxqu(ip1jmp1)
    141141  REAL :: zz(ip1jmp1)
    142   REAL :: adxqu(ip1jmp1),dxqmax(ip1jmp1,llm)
    143   REAL :: u_mq(ip1jmp1,llm)
     142  REAL :: adxqu(ip1jmp1), dxqmax(ip1jmp1, llm)
     143  REAL :: u_mq(ip1jmp1, llm)
    144144
    145145  ! CRisi
    146   REAL :: masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot)
    147   INTEGER :: ifils,iq2 ! CRisi
     146  REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot)
     147  INTEGER :: ifils, iq2 ! CRisi
    148148
    149149  LOGICAL, SAVE :: first
    … …  
    152152  !   calcul de la pente a droite et a gauche de la maille
    153153
    154 
    155154  IF (pente_max>-1.e-5) THEN
    156155    ! IF (pente_max.gt.10) THEN
    157156
    158   !   calcul des pentes avec limitation, Van Leer scheme I:
    159   !   -----------------------------------------------------
    160 
    161   !   calcul de la pente aux points u
    162      DO l = 1, llm
    163         DO ij=iip2,ip1jm-1
    164            dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    165         ENDDO
    166         DO ij=iip1+iip1,ip1jm,iip1
    167            dxqu(ij)=dxqu(ij-iim)
    168            ! sigu(ij)=sigu(ij-iim)
    169         ENDDO
    170 
    171         DO ij=iip2,ip1jm
    172            adxqu(ij)=abs(dxqu(ij))
    173         ENDDO
    174 
    175   !   calcul de la pente maximum dans la maille en valeur absolue
    176 
    177         DO ij=iip2+1,ip1jm
    178            dxqmax(ij,l)=pente_max* &
    179                  min(adxqu(ij-1),adxqu(ij))
    180   ! limitation subtile
    181   !    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
    182 
    183 
    184         ENDDO
    185 
    186         DO ij=iip1+iip1,ip1jm,iip1
    187            dxqmax(ij-iim,l)=dxqmax(ij,l)
    188         ENDDO
    189 
    190         DO ij=iip2+1,ip1jm
    191            IF(dxqu(ij-1)*dxqu(ij)>0) THEN
    192               dxq(ij,l)=dxqu(ij-1)+dxqu(ij)
    193            ELSE
    194   !   extremum local
    195               dxq(ij,l)=0.
    196            ENDIF
    197            dxq(ij,l)=0.5*dxq(ij,l)
    198            dxq(ij,l)= &
    199                  sign(min(abs(dxq(ij,l)),dxqmax(ij,l)),dxq(ij,l))
    200         ENDDO
    201 
    202      ENDDO ! l=1,llm
    203   !print*,'Ok calcul des pentes'
     157    !   calcul des pentes avec limitation, Van Leer scheme I:
     158    !   -----------------------------------------------------
     159
     160    !   calcul de la pente aux points u
     161    DO l = 1, llm
     162      DO ij = iip2, ip1jm - 1
     163        dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq)
     164      ENDDO
     165      DO ij = iip1 + iip1, ip1jm, iip1
     166        dxqu(ij) = dxqu(ij - iim)
     167        ! sigu(ij)=sigu(ij-iim)
     168      ENDDO
     169
     170      DO ij = iip2, ip1jm
     171        adxqu(ij) = abs(dxqu(ij))
     172      ENDDO
     173
     174      !   calcul de la pente maximum dans la maille en valeur absolue
     175
     176      DO ij = iip2 + 1, ip1jm
     177        dxqmax(ij, l) = pente_max * &
     178                min(adxqu(ij - 1), adxqu(ij))
     179        ! limitation subtile
     180        !    ,      min(adxqu(ij-1)/sigu(ij-1),adxqu(ij)/(1.-sigu(ij)))
     181
     182      ENDDO
     183
     184      DO ij = iip1 + iip1, ip1jm, iip1
     185        dxqmax(ij - iim, l) = dxqmax(ij, l)
     186      ENDDO
     187
     188      DO ij = iip2 + 1, ip1jm
     189        IF(dxqu(ij - 1) * dxqu(ij)>0) THEN
     190          dxq(ij, l) = dxqu(ij - 1) + dxqu(ij)
     191        ELSE
     192          !   extremum local
     193          dxq(ij, l) = 0.
     194        ENDIF
     195        dxq(ij, l) = 0.5 * dxq(ij, l)
     196        dxq(ij, l) = &
     197                sign(min(abs(dxq(ij, l)), dxqmax(ij, l)), dxq(ij, l))
     198      ENDDO
     199
     200    ENDDO ! l=1,llm
     201    !print*,'Ok calcul des pentes'
    204202
    205203  ELSE ! (pente_max.lt.-1.e-5)
    206204
    207   !   Pentes produits:
    208   !   ----------------
    209 
    210      DO l = 1, llm
    211         DO ij=iip2,ip1jm-1
    212            dxqu(ij)=q(ij+1,l,iq)-q(ij,l,iq)
    213         ENDDO
    214         DO ij=iip1+iip1,ip1jm,iip1
    215            dxqu(ij)=dxqu(ij-iim)
    216         ENDDO
    217 
    218         DO ij=iip2+1,ip1jm
    219            zz(ij)=dxqu(ij-1)*dxqu(ij)
    220            zz(ij)=zz(ij)+zz(ij)
    221            IF(zz(ij)>0) THEN
    222               dxq(ij,l)=zz(ij)/(dxqu(ij-1)+dxqu(ij))
    223            ELSE
    224   !   extremum local
    225               dxq(ij,l)=0.
    226            ENDIF
    227         ENDDO
    228 
    229      ENDDO
     205    !   Pentes produits:
     206    !   ----------------
     207
     208    DO l = 1, llm
     209      DO ij = iip2, ip1jm - 1
     210        dxqu(ij) = q(ij + 1, l, iq) - q(ij, l, iq)
     211      ENDDO
     212      DO ij = iip1 + iip1, ip1jm, iip1
     213        dxqu(ij) = dxqu(ij - iim)
     214      ENDDO
     215
     216      DO ij = iip2 + 1, ip1jm
     217        zz(ij) = dxqu(ij - 1) * dxqu(ij)
     218        zz(ij) = zz(ij) + zz(ij)
     219        IF(zz(ij)>0) THEN
     220          dxq(ij, l) = zz(ij) / (dxqu(ij - 1) + dxqu(ij))
     221        ELSE
     222          !   extremum local
     223          dxq(ij, l) = 0.
     224        ENDIF
     225      ENDDO
     226
     227    ENDDO
    230228
    231229  ENDIF ! (pente_max.lt.-1.e-5)
    … …  
    234232  !   -----------------------------
    235233
    236   DO l=1,llm
    237      DO ij=iip1+iip1,ip1jm,iip1
    238         dxq(ij-iim,l)=dxq(ij,l)
    239      ENDDO
    240      DO ij=1,ip1jmp1
    241         iadvplus(ij,l)=0
    242      ENDDO
     234  DO l = 1, llm
     235    DO ij = iip1 + iip1, ip1jm, iip1
     236      dxq(ij - iim, l) = dxq(ij, l)
     237    ENDDO
     238    DO ij = 1, ip1jmp1
     239      iadvplus(ij, l) = 0
     240    ENDDO
    243241
    244242  ENDDO
    … …  
    247245  !   on cumule le flux correspondant a toutes les mailles dont la masse
    248246  !   au travers de la paroi pENDant le pas de temps.
    249   DO l=1,llm
    250    DO ij=iip2,ip1jm-1
    251       IF (u_m(ij,l)>0.) THEN
    252          zdum(ij,l)=1.-u_m(ij,l)/masse(ij,l,iq)
    253          u_mq(ij,l)=u_m(ij,l)*(q(ij,l,iq)+0.5*zdum(ij,l)*dxq(ij,l))
     247  DO l = 1, llm
     248    DO ij = iip2, ip1jm - 1
     249      IF (u_m(ij, l)>0.) THEN
     250        zdum(ij, l) = 1. - u_m(ij, l) / masse(ij, l, iq)
     251        u_mq(ij, l) = u_m(ij, l) * (q(ij, l, iq) + 0.5 * zdum(ij, l) * dxq(ij, l))
    254252      ELSE
    255          zdum(ij,l)=1.+u_m(ij,l)/masse(ij+1,l,iq)
    256          u_mq(ij,l)=u_m(ij,l)*(q(ij+1,l,iq) &
    257                -0.5*zdum(ij,l)*dxq(ij+1,l))
     253        zdum(ij, l) = 1. + u_m(ij, l) / masse(ij + 1, l, iq)
     254        u_mq(ij, l) = u_m(ij, l) * (q(ij + 1, l, iq) &
     255                - 0.5 * zdum(ij, l) * dxq(ij + 1, l))
    258256      ENDIF
    259    ENDDO
     257    ENDDO
    260258  ENDDO
    261259
    262260  !   detection des points ou on advecte plus que la masse de la
    263261  !   maille
    264   DO l=1,llm
    265      DO ij=iip2,ip1jm-1
    266         IF(zdum(ij,l)<0) THEN
    267            iadvplus(ij,l)=1
    268            u_mq(ij,l)=0.
    269         ENDIF
    270      ENDDO
    271   ENDDO
    272   DO l=1,llm
    273    DO ij=iip1+iip1,ip1jm,iip1
    274       iadvplus(ij,l)=iadvplus(ij-iim,l)
    275    ENDDO
     262  DO l = 1, llm
     263    DO ij = iip2, ip1jm - 1
     264      IF(zdum(ij, l)<0) THEN
     265        iadvplus(ij, l) = 1
     266        u_mq(ij, l) = 0.
     267      ENDIF
     268    ENDDO
     269  ENDDO
     270  DO l = 1, llm
     271    DO ij = iip1 + iip1, ip1jm, iip1
     272      iadvplus(ij, l) = iadvplus(ij - iim, l)
     273    ENDDO
    276274  ENDDO
    277275
    … …  
    283281  !  calcul du nombre de maille sur lequel on advecte plus que la maille.
    284282
    285   n0=0
    286   DO l=1,llm
    287      nl(l)=0
    288      DO ij=iip2,ip1jm
    289         nl(l)=nl(l)+iadvplus(ij,l)
    290      ENDDO
    291      n0=n0+nl(l)
     283  n0 = 0
     284  DO l = 1, llm
     285    nl(l) = 0
     286    DO ij = iip2, ip1jm
     287      nl(l) = nl(l) + iadvplus(ij, l)
     288    ENDDO
     289    n0 = n0 + nl(l)
    292290  ENDDO
    293291
    294292  IF(n0>0) THEN
    295   IF (prt_level > 2) PRINT *, &
    296         'Nombre de points pour lesquels on advect plus que le' &
    297         ,'contenu de la maille : ',n0
    298 
    299      DO l=1,llm
    300         IF(nl(l)>0) THEN
    301            iju=0
    302   !   indicage des mailles concernees par le traitement special
    303            DO ij=iip2,ip1jm
    304               IF(iadvplus(ij,l)==1.AND.mod(ij,iip1)/=0) THEN
    305                  iju=iju+1
    306                  indu(iju)=ij
    307               ENDIF
    308            ENDDO
    309            niju=iju
    310 
    311   !  traitement des mailles
    312            DO iju=1,niju
    313               ij=indu(iju)
    314               j=(ij-1)/iip1+1
    315               zu_m=u_m(ij,l)
    316               u_mq(ij,l)=0.
    317               IF(zu_m>0.) THEN
    318                  ijq=ij
    319                  i=ijq-(j-1)*iip1
    320   !   accumulation pour les mailles completements advectees
    321                  do while(zu_m>masse(ijq,l,iq))
    322                     u_mq(ij,l)=u_mq(ij,l)+q(ijq,l,iq) &
    323                           *masse(ijq,l,iq)
    324                     zu_m=zu_m-masse(ijq,l,iq)
    325                     i=mod(i-2+iim,iim)+1
    326                     ijq=(j-1)*iip1+i
    327                  ENDDO
    328   !   ajout de la maille non completement advectee
    329                  u_mq(ij,l)=u_mq(ij,l)+zu_m* &
    330                        (q(ijq,l,iq)+0.5*(1.-zu_m/masse(ijq,l,iq)) &
    331                        *dxq(ijq,l))
    332               ELSE
    333                  ijq=ij+1
    334                  i=ijq-(j-1)*iip1
    335   !   accumulation pour les mailles completements advectees
    336                  do while(-zu_m>masse(ijq,l,iq))
    337                     u_mq(ij,l)=u_mq(ij,l)-q(ijq,l,iq) &
    338                           *masse(ijq,l,iq)
    339                     zu_m=zu_m+masse(ijq,l,iq)
    340                     i=mod(i,iim)+1
    341                     ijq=(j-1)*iip1+i
    342                  ENDDO
    343   !   ajout de la maille non completement advectee
    344                  u_mq(ij,l)=u_mq(ij,l)+zu_m*(q(ijq,l,iq)- &
    345                        0.5*(1.+zu_m/masse(ijq,l,iq))*dxq(ijq,l))
    346               ENDIF
    347            ENDDO
    348         ENDIF
    349      ENDDO
     293    IF (prt_level > 2) PRINT *, &
     294            'Nombre de points pour lesquels on advect plus que le' &
     295            , 'contenu de la maille : ', n0
     296
     297    DO l = 1, llm
     298      IF(nl(l)>0) THEN
     299        iju = 0
     300        !   indicage des mailles concernees par le traitement special
     301        DO ij = iip2, ip1jm
     302          IF(iadvplus(ij, l)==1.AND.mod(ij, iip1)/=0) THEN
     303            iju = iju + 1
     304            indu(iju) = ij
     305          ENDIF
     306        ENDDO
     307        niju = iju
     308
     309        !  traitement des mailles
     310        DO iju = 1, niju
     311          ij = indu(iju)
     312          j = (ij - 1) / iip1 + 1
     313          zu_m = u_m(ij, l)
     314          u_mq(ij, l) = 0.
     315          IF(zu_m>0.) THEN
     316            ijq = ij
     317            i = ijq - (j - 1) * iip1
     318            !   accumulation pour les mailles completements advectees
     319            do while(zu_m>masse(ijq, l, iq))
     320              u_mq(ij, l) = u_mq(ij, l) + q(ijq, l, iq) &
     321                      * masse(ijq, l, iq)
     322              zu_m = zu_m - masse(ijq, l, iq)
     323              i = mod(i - 2 + iim, iim) + 1
     324              ijq = (j - 1) * iip1 + i
     325            ENDDO
     326            !   ajout de la maille non completement advectee
     327            u_mq(ij, l) = u_mq(ij, l) + zu_m * &
     328                    (q(ijq, l, iq) + 0.5 * (1. - zu_m / masse(ijq, l, iq)) &
     329                            * dxq(ijq, l))
     330          ELSE
     331            ijq = ij + 1
     332            i = ijq - (j - 1) * iip1
     333            !   accumulation pour les mailles completements advectees
     334            do while(-zu_m>masse(ijq, l, iq))
     335              u_mq(ij, l) = u_mq(ij, l) - q(ijq, l, iq) &
     336                      * masse(ijq, l, iq)
     337              zu_m = zu_m + masse(ijq, l, iq)
     338              i = mod(i, iim) + 1
     339              ijq = (j - 1) * iip1 + i
     340            ENDDO
     341            !   ajout de la maille non completement advectee
     342            u_mq(ij, l) = u_mq(ij, l) + zu_m * (q(ijq, l, iq) - &
     343                    0.5 * (1. + zu_m / masse(ijq, l, iq)) * dxq(ijq, l))
     344          ENDIF
     345        ENDDO
     346      ENDIF
     347    ENDDO
    350348  ENDIF  ! n0.gt.0
    351349
    … …  
    353351  !   bouclage en latitude
    354352  !print*,'cvant bouclage en latitude'
    355   DO l=1,llm
    356     DO ij=iip1+iip1,ip1jm,iip1
    357        u_mq(ij,l)=u_mq(ij-iim,l)
     353  DO l = 1, llm
     354    DO ij = iip1 + iip1, ip1jm, iip1
     355      u_mq(ij, l) = u_mq(ij - iim, l)
    358356    ENDDO
    359357  ENDDO
    … …  
    363361  !WRITE(*,*) 'vlsplt 326: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen
    364362
    365   do ifils=1,tracers(iq)%nqDescen
    366     iq2=tracers(iq)%iqDescen(ifils)
    367     DO l=1,llm
    368       DO ij=iip2,ip1jm
     363  do ifils = 1, tracers(iq)%nqDescen
     364    iq2 = tracers(iq)%iqDescen(ifils)
     365    DO l = 1, llm
     366      DO ij = iip2, ip1jm
    369367        ! On a besoin de q et masse seulement entre iip2 et ip1jm
    370368        !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
    371   !           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     369        !           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
    372370        !Mvals: veiller a ce qu'on n'ait pas de denominateur nul
    373         masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
    374         IF (q(ij,l,iq)>min_qParent) THEN
    375           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     371        masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
     372        IF (q(ij, l, iq)>min_qParent) THEN
     373          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
    376374        else
    377           Ratio(ij,l,iq2)=min_ratio
     375          Ratio(ij, l, iq2) = min_ratio
    378376        endif
    379377      enddo
    380378    enddo
    381379  enddo
    382   do ifils=1,tracers(iq)%nqChildren
    383     iq2=tracers(iq)%iqDescen(ifils)
    384     CALL vlx(Ratio,pente_max,masseq,u_mq,iq2)
     380  do ifils = 1, tracers(iq)%nqChildren
     381    iq2 = tracers(iq)%iqDescen(ifils)
     382    CALL vlx(Ratio, pente_max, masseq, u_mq, iq2)
    385383  enddo
    386384  ! end CRisi
    … …  
    389387  !   calcul des tENDances
    390388
    391   DO l=1,llm
    392      DO ij=iip2+1,ip1jm
    393         !MVals: veiller a ce qu'on ait pas de denominateur nul
    394         new_m=max(masse(ij,l,iq)+u_m(ij-1,l)-u_m(ij,l),min_qMass)
    395         q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+ &
    396               u_mq(ij-1,l)-u_mq(ij,l)) &
    397               /new_m
    398         masse(ij,l,iq)=new_m
    399      ENDDO
    400      DO ij=iip1+iip1,ip1jm,iip1
    401         q(ij-iim,l,iq)=q(ij,l,iq)
    402         masse(ij-iim,l,iq)=masse(ij,l,iq)
    403      ENDDO
     389  DO l = 1, llm
     390    DO ij = iip2 + 1, ip1jm
     391      !MVals: veiller a ce qu'on ait pas de denominateur nul
     392      new_m = max(masse(ij, l, iq) + u_m(ij - 1, l) - u_m(ij, l), min_qMass)
     393      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + &
     394              u_mq(ij - 1, l) - u_mq(ij, l)) &
     395              / new_m
     396      masse(ij, l, iq) = new_m
     397    ENDDO
     398    DO ij = iip1 + iip1, ip1jm, iip1
     399      q(ij - iim, l, iq) = q(ij, l, iq)
     400      masse(ij - iim, l, iq) = masse(ij, l, iq)
     401    ENDDO
    404402  ENDDO
    405403
    … …  
    407405  ! On calcule q entre iip2+1,ip1jm -> on fait pareil pour ratio
    408406  ! puis on boucle en longitude
    409   do ifils=1,tracers(iq)%nqDescen
    410     iq2=tracers(iq)%iqDescen(ifils)
    411     DO l=1,llm
    412       DO ij=iip2+1,ip1jm
    413         q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
     407  do ifils = 1, tracers(iq)%nqDescen
     408    iq2 = tracers(iq)%iqDescen(ifils)
     409    DO l = 1, llm
     410      DO ij = iip2 + 1, ip1jm
     411        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
    414412      enddo
    415       DO ij=iip1+iip1,ip1jm,iip1
    416          q(ij-iim,l,iq2)=q(ij,l,iq2)
     413      DO ij = iip1 + iip1, ip1jm, iip1
     414        q(ij - iim, l, iq2) = q(ij, l, iq2)
    417415      enddo ! DO ij=ijb+iip1-1,ije,iip1
    418416    enddo !DO l=1,llm
    419417  enddo
    420418
    421 
    422419END SUBROUTINE vlx
    423 RECURSIVE SUBROUTINE vly(q,pente_max,masse,masse_adv_v,iq)
    424   USE infotrac, ONLY: nqtot,tracers, & ! CRisi
    425         min_qParent,min_qMass,min_ratio ! MVals et CRisi
     420RECURSIVE SUBROUTINE vly(q, pente_max, masse, masse_adv_v, iq)
     421  USE infotrac, ONLY: nqtot, tracers, & ! CRisi
     422          min_qParent, min_qMass, min_ratio ! MVals et CRisi
    426423  !
    427424  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    445442  !   Arguments:
    446443  !   ----------
    447   REAL :: masse(ip1jmp1,llm,nqtot),pente_max
    448   REAL :: masse_adv_v( ip1jm,llm)
    449   REAL :: q(ip1jmp1,llm,nqtot)
     444  REAL :: masse(ip1jmp1, llm, nqtot), pente_max
     445  REAL :: masse_adv_v(ip1jm, llm)
     446  REAL :: q(ip1jmp1, llm, nqtot)
    450447  INTEGER :: iq ! CRisi
    451448  !
    … …  
    453450  !   ---------
    454451  !
    455   INTEGER :: i,ij,l
    456   !
    457   REAL :: airej2,airejjm,airescb(iim),airesch(iim)
    458   REAL :: dyq(ip1jmp1,llm),dyqv(ip1jm)
    459   REAL :: adyqv(ip1jm),dyqmax(ip1jmp1)
    460   REAL :: qbyv(ip1jm,llm)
    461 
    462   REAL :: qpns,qpsn,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
     452  INTEGER :: i, ij, l
     453  !
     454  REAL :: airej2, airejjm, airescb(iim), airesch(iim)
     455  REAL :: dyq(ip1jmp1, llm), dyqv(ip1jm)
     456  REAL :: adyqv(ip1jm), dyqmax(ip1jmp1)
     457  REAL :: qbyv(ip1jm, llm)
     458
     459  REAL :: qpns, qpsn, dyn1, dys1, dyn2, dys2, newmasse, fn, fs
    463460  LOGICAL, SAVE :: first
    464461
    465   REAL :: convpn,convps,convmpn,convmps
    466   REAL :: massepn,masseps,qpn,qps
    467   REAL :: sinlon(iip1),sinlondlon(iip1)
    468   REAL :: coslon(iip1),coslondlon(iip1)
    469   SAVE sinlon,coslon,sinlondlon,coslondlon
    470   SAVE airej2,airejjm
    471 
    472   REAL :: masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
    473   INTEGER :: ifils,iq2 ! CRisi
     462  REAL :: convpn, convps, convmpn, convmps
     463  REAL :: massepn, masseps, qpn, qps
     464  REAL :: sinlon(iip1), sinlondlon(iip1)
     465  REAL :: coslon(iip1), coslondlon(iip1)
     466  SAVE sinlon, coslon, sinlondlon, coslondlon
     467  SAVE airej2, airejjm
     468
     469  REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot) ! CRisi
     470  INTEGER :: ifils, iq2 ! CRisi
    474471
    475472  !
    … …  
    482479
    483480  IF(first) THEN
    484      PRINT*,'Shema  Amont nouveau  appele dans  Vanleer   '
    485      first=.FALSE.
    486      do i=2,iip1
    487         coslon(i)=cos(rlonv(i))
    488         sinlon(i)=sin(rlonv(i))
    489         coslondlon(i)=coslon(i)*(rlonu(i)-rlonu(i-1))/pi
    490         sinlondlon(i)=sinlon(i)*(rlonu(i)-rlonu(i-1))/pi
    491      ENDDO
    492      coslon(1)=coslon(iip1)
    493      coslondlon(1)=coslondlon(iip1)
    494      sinlon(1)=sinlon(iip1)
    495      sinlondlon(1)=sinlondlon(iip1)
    496      airej2 = SSUM( iim, aire(iip2), 1 )
    497      airejjm= SSUM( iim, aire(ip1jm -iim), 1 )
     481    PRINT*, 'Shema  Amont nouveau  appele dans  Vanleer   '
     482    first = .FALSE.
     483    do i = 2, iip1
     484      coslon(i) = cos(rlonv(i))
     485      sinlon(i) = sin(rlonv(i))
     486      coslondlon(i) = coslon(i) * (rlonu(i) - rlonu(i - 1)) / pi
     487      sinlondlon(i) = sinlon(i) * (rlonu(i) - rlonu(i - 1)) / pi
     488    ENDDO
     489    coslon(1) = coslon(iip1)
     490    coslondlon(1) = coslondlon(iip1)
     491    sinlon(1) = sinlon(iip1)
     492    sinlondlon(1) = sinlondlon(iip1)
     493    airej2 = SSUM(iim, aire(iip2), 1)
     494    airejjm = SSUM(iim, aire(ip1jm - iim), 1)
    498495  ENDIF
    499496
    … …  
    502499
    503500  DO l = 1, llm
    504   !
    505   !   --------------------------------
    506   !  CALCUL EN LATITUDE
    507   !   --------------------------------
    508 
    509   !   On commence par calculer la valeur du traceur moyenne sur le premier cercle
    510   !   de latitude autour du pole (qpns pour le pole nord et qpsn pour
    511   !    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
    512 
    513   DO i = 1, iim
    514   airescb(i) = aire(i+ iip1) * q(i+ iip1,l,iq)
    515   airesch(i) = aire(i+ ip1jm- iip1) * q(i+ ip1jm- iip1,l,iq)
    516   ENDDO
    517   qpns   = SSUM( iim,  airescb ,1 ) / airej2
    518   qpsn   = SSUM( iim,  airesch ,1 ) / airejjm
    519 
    520   !   calcul des pentes aux points v
    521 
    522   DO ij=1,ip1jm
    523      dyqv(ij)=q(ij,l,iq)-q(ij+iip1,l,iq)
    524      adyqv(ij)=abs(dyqv(ij))
    525   ENDDO
    526 
    527   !   calcul des pentes aux points scalaires
    528 
    529   DO ij=iip2,ip1jm
    530      dyq(ij,l)=.5*(dyqv(ij-iip1)+dyqv(ij))
    531      dyqmax(ij)=min(adyqv(ij-iip1),adyqv(ij))
    532      dyqmax(ij)=pente_max*dyqmax(ij)
    533   ENDDO
    534 
    535   !   calcul des pentes aux poles
    536 
    537   DO ij=1,iip1
    538      dyq(ij,l)=qpns-q(ij+iip1,l,iq)
    539      dyq(ip1jm+ij,l)=q(ip1jm+ij-iip1,l,iq)-qpsn
    540   ENDDO
    541 
    542   !   filtrage de la derivee
    543   dyn1=0.
    544   dys1=0.
    545   dyn2=0.
    546   dys2=0.
    547   DO ij=1,iim
    548      dyn1=dyn1+sinlondlon(ij)*dyq(ij,l)
    549      dys1=dys1+sinlondlon(ij)*dyq(ip1jm+ij,l)
    550      dyn2=dyn2+coslondlon(ij)*dyq(ij,l)
    551      dys2=dys2+coslondlon(ij)*dyq(ip1jm+ij,l)
    552   ENDDO
    553   DO ij=1,iip1
    554      dyq(ij,l)=dyn1*sinlon(ij)+dyn2*coslon(ij)
    555      dyq(ip1jm+ij,l)=dys1*sinlon(ij)+dys2*coslon(ij)
    556   ENDDO
    557 
    558   !   calcul des pentes limites aux poles
    559 
    560   goto 8888
    561   fn=1.
    562   fs=1.
    563   DO ij=1,iim
    564      IF(pente_max*adyqv(ij)<abs(dyq(ij,l))) THEN
    565         fn=min(pente_max*adyqv(ij)/abs(dyq(ij,l)),fn)
    566      ENDIF
    567   IF(pente_max*adyqv(ij+ip1jm-iip1)<abs(dyq(ij+ip1jm,l))) THEN
    568      fs=min(pente_max*adyqv(ij+ip1jm-iip1)/abs(dyq(ij+ip1jm,l)),fs)
    569      ENDIF
    570   ENDDO
    571   DO ij=1,iip1
    572      dyq(ij,l)=fn*dyq(ij,l)
    573      dyq(ip1jm+ij,l)=fs*dyq(ip1jm+ij,l)
    574   ENDDO
    575 8888   continue
    576   DO ij=1,iip1
    577      dyq(ij,l)=0.
    578      dyq(ip1jm+ij,l)=0.
    579   ENDDO
    580 
    581   !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    582   !  En memoire de dIFferents tests sur la
    583   !  limitation des pentes aux poles.
    584   !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    585   ! PRINT*,dyq(1)
    586   ! PRINT*,dyqv(iip1+1)
    587   ! appn=abs(dyq(1)/dyqv(iip1+1))
    588   ! PRINT*,dyq(ip1jm+1)
    589   ! PRINT*,dyqv(ip1jm-iip1+1)
    590   ! apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    591   ! DO ij=2,iim
    592   !    appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
    593   !    apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    594   ! ENDDO
    595   ! appn=min(pente_max/appn,1.)
    596   ! apps=min(pente_max/apps,1.)
    597   !
    598   !
    599   !   cas ou on a un extremum au pole
    600   !
    601   ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    602   !    &   appn=0.
    603   ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    604   !    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    605   !    &   apps=0.
    606   !
    607   !   limitation des pentes aux poles
    608   ! DO ij=1,iip1
    609   !    dyq(ij)=appn*dyq(ij)
    610   !    dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    611   ! ENDDO
    612   !
    613   !   test
    614   !  DO ij=1,iip1
    615   !     dyq(iip1+ij)=0.
    616   !     dyq(ip1jm+ij-iip1)=0.
    617   !  ENDDO
    618   !  DO ij=1,ip1jmp1
    619   !     dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
    620   !  ENDDO
    621   !
    622   ! changement 10 07 96
    623   ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    624   !    &   THEN
    625   !    DO ij=1,iip1
    626   !       dyqmax(ij)=0.
    627   !    ENDDO
    628   ! ELSE
    629   !    DO ij=1,iip1
    630   !       dyqmax(ij)=pente_max*abs(dyqv(ij))
    631   !    ENDDO
    632   ! ENDIF
    633   !
    634   ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    635   !    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    636   !    &THEN
    637   !    DO ij=ip1jm+1,ip1jmp1
    638   !       dyqmax(ij)=0.
    639   !    ENDDO
    640   ! ELSE
    641   !    DO ij=ip1jm+1,ip1jmp1
    642   !       dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
    643   !    ENDDO
    644   ! ENDIF
    645   !   fin changement 10 07 96
    646   !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    647 
    648   !   calcul des pentes limitees
    649 
    650   DO ij=iip2,ip1jm
    651      IF(dyqv(ij)*dyqv(ij-iip1)>0.) THEN
    652         dyq(ij,l)=sign(min(abs(dyq(ij,l)),dyqmax(ij)),dyq(ij,l))
    653      ELSE
    654         dyq(ij,l)=0.
    655      ENDIF
    656   ENDDO
     501    !
     502    !   --------------------------------
     503    !  CALCUL EN LATITUDE
     504    !   --------------------------------
     505
     506    !   On commence par calculer la valeur du traceur moyenne sur le premier cercle
     507    !   de latitude autour du pole (qpns pour le pole nord et qpsn pour
     508    !    le pole nord) qui sera utilisee pour evaluer les pentes au pole.
     509
     510    DO i = 1, iim
     511      airescb(i) = aire(i + iip1) * q(i + iip1, l, iq)
     512      airesch(i) = aire(i + ip1jm - iip1) * q(i + ip1jm - iip1, l, iq)
     513    ENDDO
     514    qpns = SSUM(iim, airescb, 1) / airej2
     515    qpsn = SSUM(iim, airesch, 1) / airejjm
     516
     517    !   calcul des pentes aux points v
     518
     519    DO ij = 1, ip1jm
     520      dyqv(ij) = q(ij, l, iq) - q(ij + iip1, l, iq)
     521      adyqv(ij) = abs(dyqv(ij))
     522    ENDDO
     523
     524    !   calcul des pentes aux points scalaires
     525
     526    DO ij = iip2, ip1jm
     527      dyq(ij, l) = .5 * (dyqv(ij - iip1) + dyqv(ij))
     528      dyqmax(ij) = min(adyqv(ij - iip1), adyqv(ij))
     529      dyqmax(ij) = pente_max * dyqmax(ij)
     530    ENDDO
     531
     532    !   calcul des pentes aux poles
     533
     534    DO ij = 1, iip1
     535      dyq(ij, l) = qpns - q(ij + iip1, l, iq)
     536      dyq(ip1jm + ij, l) = q(ip1jm + ij - iip1, l, iq) - qpsn
     537    ENDDO
     538
     539    !   filtrage de la derivee
     540    dyn1 = 0.
     541    dys1 = 0.
     542    dyn2 = 0.
     543    dys2 = 0.
     544    DO ij = 1, iim
     545      dyn1 = dyn1 + sinlondlon(ij) * dyq(ij, l)
     546      dys1 = dys1 + sinlondlon(ij) * dyq(ip1jm + ij, l)
     547      dyn2 = dyn2 + coslondlon(ij) * dyq(ij, l)
     548      dys2 = dys2 + coslondlon(ij) * dyq(ip1jm + ij, l)
     549    ENDDO
     550    DO ij = 1, iip1
     551      dyq(ij, l) = dyn1 * sinlon(ij) + dyn2 * coslon(ij)
     552      dyq(ip1jm + ij, l) = dys1 * sinlon(ij) + dys2 * coslon(ij)
     553    ENDDO
     554
     555    !   calcul des pentes limites aux poles
     556
     557    goto 8888
     558    fn = 1.
     559    fs = 1.
     560    DO ij = 1, iim
     561      IF(pente_max * adyqv(ij)<abs(dyq(ij, l))) THEN
     562        fn = min(pente_max * adyqv(ij) / abs(dyq(ij, l)), fn)
     563      ENDIF
     564      IF(pente_max * adyqv(ij + ip1jm - iip1)<abs(dyq(ij + ip1jm, l))) THEN
     565        fs = min(pente_max * adyqv(ij + ip1jm - iip1) / abs(dyq(ij + ip1jm, l)), fs)
     566      ENDIF
     567    ENDDO
     568    DO ij = 1, iip1
     569      dyq(ij, l) = fn * dyq(ij, l)
     570      dyq(ip1jm + ij, l) = fs * dyq(ip1jm + ij, l)
     571    ENDDO
     572    8888   continue
     573    DO ij = 1, iip1
     574      dyq(ij, l) = 0.
     575      dyq(ip1jm + ij, l) = 0.
     576    ENDDO
     577
     578    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     579    !  En memoire de dIFferents tests sur la
     580    !  limitation des pentes aux poles.
     581    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     582    ! PRINT*,dyq(1)
     583    ! PRINT*,dyqv(iip1+1)
     584    ! appn=abs(dyq(1)/dyqv(iip1+1))
     585    ! PRINT*,dyq(ip1jm+1)
     586    ! PRINT*,dyqv(ip1jm-iip1+1)
     587    ! apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     588    ! DO ij=2,iim
     589    !    appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
     590    !    apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
     591    ! ENDDO
     592    ! appn=min(pente_max/appn,1.)
     593    ! apps=min(pente_max/apps,1.)
     594    !
     595    !
     596    !   cas ou on a un extremum au pole
     597    !
     598    ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
     599    !    &   appn=0.
     600    ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
     601    !    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
     602    !    &   apps=0.
     603    !
     604    !   limitation des pentes aux poles
     605    ! DO ij=1,iip1
     606    !    dyq(ij)=appn*dyq(ij)
     607    !    dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
     608    ! ENDDO
     609    !
     610    !   test
     611    !  DO ij=1,iip1
     612    !     dyq(iip1+ij)=0.
     613    !     dyq(ip1jm+ij-iip1)=0.
     614    !  ENDDO
     615    !  DO ij=1,ip1jmp1
     616    !     dyq(ij)=dyq(ij)*cos(rlatu((ij-1)/iip1+1))
     617    !  ENDDO
     618    !
     619    ! changement 10 07 96
     620    ! IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
     621    !    &   THEN
     622    !    DO ij=1,iip1
     623    !       dyqmax(ij)=0.
     624    !    ENDDO
     625    ! ELSE
     626    !    DO ij=1,iip1
     627    !       dyqmax(ij)=pente_max*abs(dyqv(ij))
     628    !    ENDDO
     629    ! ENDIF
     630    !
     631    ! IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
     632    !    & dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
     633    !    &THEN
     634    !    DO ij=ip1jm+1,ip1jmp1
     635    !       dyqmax(ij)=0.
     636    !    ENDDO
     637    ! ELSE
     638    !    DO ij=ip1jm+1,ip1jmp1
     639    !       dyqmax(ij)=pente_max*abs(dyqv(ij-iip1))
     640    !    ENDDO
     641    ! ENDIF
     642    !   fin changement 10 07 96
     643    !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     644
     645    !   calcul des pentes limitees
     646
     647    DO ij = iip2, ip1jm
     648      IF(dyqv(ij) * dyqv(ij - iip1)>0.) THEN
     649        dyq(ij, l) = sign(min(abs(dyq(ij, l)), dyqmax(ij)), dyq(ij, l))
     650      ELSE
     651        dyq(ij, l) = 0.
     652      ENDIF
     653    ENDDO
    657654
    658655  ENDDO
    659656
    660657  ! !WRITE(*,*) 'vly 756'
    661   DO l=1,llm
    662    DO ij=1,ip1jm
    663       IF(masse_adv_v(ij,l)>0) THEN
    664           qbyv(ij,l)=q(ij+iip1,l,iq)+dyq(ij+iip1,l)* &
    665                 0.5*(1.-masse_adv_v(ij,l) &
    666                 /masse(ij+iip1,l,iq))
     658  DO l = 1, llm
     659    DO ij = 1, ip1jm
     660      IF(masse_adv_v(ij, l)>0) THEN
     661        qbyv(ij, l) = q(ij + iip1, l, iq) + dyq(ij + iip1, l) * &
     662                0.5 * (1. - masse_adv_v(ij, l) &
     663                / masse(ij + iip1, l, iq))
    667664      ELSE
    668           qbyv(ij,l)=q(ij,l,iq)-dyq(ij,l)* &
    669                 0.5*(1.+masse_adv_v(ij,l) &
    670                 /masse(ij,l,iq))
     665        qbyv(ij, l) = q(ij, l, iq) - dyq(ij, l) * &
     666                0.5 * (1. + masse_adv_v(ij, l) &
     667                / masse(ij, l, iq))
    671668      ENDIF
    672       qbyv(ij,l)=masse_adv_v(ij,l)*qbyv(ij,l)
    673    ENDDO
     669      qbyv(ij, l) = masse_adv_v(ij, l) * qbyv(ij, l)
     670    ENDDO
    674671  ENDDO
    675672
    676673  ! CRisi: appel récursif de l'advection sur les fils.
    677674  ! Il faut faire ça avant d'avoir mis à jour q et masse
    678    ! WRITE(*,*) 'vly 689: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen
    679 
    680   do ifils=1,tracers(iq)%nqDescen
    681     iq2=tracers(iq)%iqDescen(ifils)
    682     DO l=1,llm
    683       DO ij=1,ip1jmp1
     675  ! WRITE(*,*) 'vly 689: iq,nqDesc(iq)=',iq,tracers(iq)%nqDescen
     676
     677  do ifils = 1, tracers(iq)%nqDescen
     678    iq2 = tracers(iq)%iqDescen(ifils)
     679    DO l = 1, llm
     680      DO ij = 1, ip1jmp1
    684681        ! ! attention, chaque fils doit avoir son masseq, sinon, le 1er
    685682        ! ! fils ecrase le masseq de ses freres.
    686683        ! !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
    687   !           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     684        !           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
    688685        ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul
    689         masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
    690         IF (q(ij,l,iq)>min_qParent) THEN
    691           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     686        masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
     687        IF (q(ij, l, iq)>min_qParent) THEN
     688          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
    692689        else
    693           Ratio(ij,l,iq2)=min_ratio
     690          Ratio(ij, l, iq2) = min_ratio
    694691        endif
    695692      enddo
    … …  
    697694  enddo
    698695
    699   do ifils=1,tracers(iq)%nqDescen
    700     iq2=tracers(iq)%iqDescen(ifils)
    701     CALL vly(Ratio,pente_max,masseq,qbyv,iq2)
    702   enddo
    703 
    704   DO l=1,llm
    705      DO ij=iip2,ip1jm
    706         newmasse=masse(ij,l,iq) &
    707               +masse_adv_v(ij,l)-masse_adv_v(ij-iip1,l)
    708         q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+qbyv(ij,l) &
    709               -qbyv(ij-iip1,l))/newmasse
    710         masse(ij,l,iq)=newmasse
    711      ENDDO
    712      convpn=SSUM(iim,qbyv(1,l),1)
    713      convmpn=ssum(iim,masse_adv_v(1,l),1)
    714      massepn=ssum(iim,masse(1,l,iq),1)
    715      qpn=0.
    716      do ij=1,iim
    717         qpn=qpn+masse(ij,l,iq)*q(ij,l,iq)
    718      enddo
    719      qpn=(qpn+convpn)/(massepn+convmpn)
    720      do ij=1,iip1
    721         q(ij,l,iq)=qpn
    722      enddo
    723      convps=-SSUM(iim,qbyv(ip1jm-iim,l),1)
    724      convmps=-ssum(iim,masse_adv_v(ip1jm-iim,l),1)
    725      masseps=ssum(iim, masse(ip1jm+1,l,iq),1)
    726      qps=0.
    727      do ij = ip1jm+1,ip1jmp1-1
    728         qps=qps+masse(ij,l,iq)*q(ij,l,iq)
    729      enddo
    730      qps=(qps+convps)/(masseps+convmps)
    731      do ij=ip1jm+1,ip1jmp1
    732         q(ij,l,iq)=qps
    733      enddo
     696  do ifils = 1, tracers(iq)%nqDescen
     697    iq2 = tracers(iq)%iqDescen(ifils)
     698    CALL vly(Ratio, pente_max, masseq, qbyv, iq2)
     699  enddo
     700
     701  DO l = 1, llm
     702    DO ij = iip2, ip1jm
     703      newmasse = masse(ij, l, iq) &
     704              + masse_adv_v(ij, l) - masse_adv_v(ij - iip1, l)
     705      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + qbyv(ij, l) &
     706              - qbyv(ij - iip1, l)) / newmasse
     707      masse(ij, l, iq) = newmasse
     708    ENDDO
     709    convpn = SSUM(iim, qbyv(1, l), 1)
     710    convmpn = ssum(iim, masse_adv_v(1, l), 1)
     711    massepn = ssum(iim, masse(1, l, iq), 1)
     712    qpn = 0.
     713    do ij = 1, iim
     714      qpn = qpn + masse(ij, l, iq) * q(ij, l, iq)
     715    enddo
     716    qpn = (qpn + convpn) / (massepn + convmpn)
     717    do ij = 1, iip1
     718      q(ij, l, iq) = qpn
     719    enddo
     720    convps = -SSUM(iim, qbyv(ip1jm - iim, l), 1)
     721    convmps = -ssum(iim, masse_adv_v(ip1jm - iim, l), 1)
     722    masseps = ssum(iim, masse(ip1jm + 1, l, iq), 1)
     723    qps = 0.
     724    do ij = ip1jm + 1, ip1jmp1 - 1
     725      qps = qps + masse(ij, l, iq) * q(ij, l, iq)
     726    enddo
     727    qps = (qps + convps) / (masseps + convmps)
     728    do ij = ip1jm + 1, ip1jmp1
     729      q(ij, l, iq) = qps
     730    enddo
    734731  ENDDO
    735732
    736733  ! retablir les fils en rapport de melange par rapport a l'air:
    737   do ifils=1,tracers(iq)%nqDescen
    738     iq2=tracers(iq)%iqDescen(ifils)
    739     DO l=1,llm
    740       DO ij=1,ip1jmp1
    741         q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
     734  do ifils = 1, tracers(iq)%nqDescen
     735    iq2 = tracers(iq)%iqDescen(ifils)
     736    DO l = 1, llm
     737      DO ij = 1, ip1jmp1
     738        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
    742739      enddo
    743740    enddo
    … …  
    746743  ! !WRITE(*,*) 'vly 853: sortie'
    747744
    748 
    749745END SUBROUTINE vly
    750 RECURSIVE SUBROUTINE vlz(q,pente_max,masse,w,iq)
    751   USE infotrac, ONLY: nqtot,tracers, & ! CRisi
    752         min_qParent,min_qMass,min_ratio ! MVals et CRisi
     746RECURSIVE SUBROUTINE vlz(q, pente_max, masse, w, iq)
     747  USE infotrac, ONLY: nqtot, tracers, & ! CRisi
     748          min_qParent, min_qMass, min_ratio ! MVals et CRisi
    753749  !
    754750  ! Auteurs:   P.Le Van, F.Hourdin, F.Forget
    … …  
    768764  !   Arguments:
    769765  !   ----------
    770   REAL :: masse(ip1jmp1,llm,nqtot),pente_max
    771   REAL :: q(ip1jmp1,llm,nqtot)
    772   REAL :: w(ip1jmp1,llm+1)
     766  REAL :: masse(ip1jmp1, llm, nqtot), pente_max
     767  REAL :: q(ip1jmp1, llm, nqtot)
     768  REAL :: w(ip1jmp1, llm + 1)
    773769  INTEGER :: iq
    774770  !
    … …  
    776772  !   ---------
    777773  !
    778   INTEGER :: ij,l
    779   !
    780   REAL :: wq(ip1jmp1,llm+1),newmasse
    781 
    782   REAL :: dzq(ip1jmp1,llm),dzqw(ip1jmp1,llm),adzqw(ip1jmp1,llm),dzqmax
     774  INTEGER :: ij, l
     775  !
     776  REAL :: wq(ip1jmp1, llm + 1), newmasse
     777
     778  REAL :: dzq(ip1jmp1, llm), dzqw(ip1jmp1, llm), adzqw(ip1jmp1, llm), dzqmax
    783779  REAL :: sigw
    784780
    785   REAL :: masseq(ip1jmp1,llm,nqtot),Ratio(ip1jmp1,llm,nqtot) ! CRisi
    786   INTEGER :: ifils,iq2 ! CRisi
     781  REAL :: masseq(ip1jmp1, llm, nqtot), Ratio(ip1jmp1, llm, nqtot) ! CRisi
     782  INTEGER :: ifils, iq2 ! CRisi
    787783
    788784#ifdef BIDON
    … …  
    795791  !    On oriente tout dans le sens de la pression c'est a dire dans le
    796792  !    sens de W
    797   DO l=2,llm
    798      DO ij=1,ip1jmp1
    799         dzqw(ij,l)=q(ij,l-1,iq)-q(ij,l,iq)
    800         adzqw(ij,l)=abs(dzqw(ij,l))
    801      ENDDO
    802   ENDDO
    803 
    804   DO l=2,llm-1
    805      DO ij=1,ip1jmp1
    806         IF(dzqw(ij,l)*dzqw(ij,l+1)>0.) THEN
    807             dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
    808         ELSE
    809             dzq(ij,l)=0.
    810         ENDIF
    811         dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
    812         dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
    813      ENDDO
     793  DO l = 2, llm
     794    DO ij = 1, ip1jmp1
     795      dzqw(ij, l) = q(ij, l - 1, iq) - q(ij, l, iq)
     796      adzqw(ij, l) = abs(dzqw(ij, l))
     797    ENDDO
     798  ENDDO
     799
     800  DO l = 2, llm - 1
     801    DO ij = 1, ip1jmp1
     802      IF(dzqw(ij, l) * dzqw(ij, l + 1)>0.) THEN
     803        dzq(ij, l) = 0.5 * (dzqw(ij, l) + dzqw(ij, l + 1))
     804      ELSE
     805        dzq(ij, l) = 0.
     806      ENDIF
     807      dzqmax = pente_max * min(adzqw(ij, l), adzqw(ij, l + 1))
     808      dzq(ij, l) = sign(min(abs(dzq(ij, l)), dzqmax), dzq(ij, l))
     809    ENDDO
    814810  ENDDO
    815811
    816812  ! !WRITE(*,*) 'vlz 954'
    817   DO ij=1,ip1jmp1
    818      dzq(ij,1)=0.
    819      dzq(ij,llm)=0.
     813  DO ij = 1, ip1jmp1
     814    dzq(ij, 1) = 0.
     815    dzq(ij, llm) = 0.
    820816  ENDDO
    821817
    … …  
    826822  ! calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
    827823
    828    DO l = 1,llm-1
    829      do  ij = 1,ip1jmp1
    830       IF(w(ij,l+1)>0.) THEN
    831          sigw=w(ij,l+1)/masse(ij,l+1,iq)
    832          wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1,iq) &
    833                +0.5*(1.-sigw)*dzq(ij,l+1))
     824  DO l = 1, llm - 1
     825    do  ij = 1, ip1jmp1
     826      IF(w(ij, l + 1)>0.) THEN
     827        sigw = w(ij, l + 1) / masse(ij, l + 1, iq)
     828        wq(ij, l + 1) = w(ij, l + 1) * (q(ij, l + 1, iq) &
     829                + 0.5 * (1. - sigw) * dzq(ij, l + 1))
    834830      ELSE
    835          sigw=w(ij,l+1)/masse(ij,l,iq)
    836          wq(ij,l+1)=w(ij,l+1)*(q(ij,l,iq)-0.5*(1.+sigw)*dzq(ij,l))
     831        sigw = w(ij, l + 1) / masse(ij, l, iq)
     832        wq(ij, l + 1) = w(ij, l + 1) * (q(ij, l, iq) - 0.5 * (1. + sigw) * dzq(ij, l))
    837833      ENDIF
    838      ENDDO
    839    ENDDO
    840 
    841    DO ij=1,ip1jmp1
    842       wq(ij,llm+1)=0.
    843       wq(ij,1)=0.
    844    ENDDO
     834    ENDDO
     835  ENDDO
     836
     837  DO ij = 1, ip1jmp1
     838    wq(ij, llm + 1) = 0.
     839    wq(ij, 1) = 0.
     840  ENDDO
    845841
    846842  ! CRisi: appel récursif de l'advection sur les fils.
    847843  ! Il faut faire ça avant d'avoir mis à jour q et masse
    848844  ! !WRITE(*,*) 'vlsplt 942: iq,nqChildren(iq)=',iq,nqChildren(iq)
    849   do ifils=1,tracers(iq)%nqDescen
    850     iq2=tracers(iq)%iqDescen(ifils)
    851     DO l=1,llm
    852       DO ij=1,ip1jmp1
     845  do ifils = 1, tracers(iq)%nqDescen
     846    iq2 = tracers(iq)%iqDescen(ifils)
     847    DO l = 1, llm
     848      DO ij = 1, ip1jmp1
    853849        ! !masseq(ij,l,iq2)=masse(ij,l,iq)*q(ij,l,iq)
    854   !           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     850        !           !Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
    855851        ! !MVals: veiller a ce qu'on n'ait pas de denominateur nul
    856         masseq(ij,l,iq2)=max(masse(ij,l,iq)*q(ij,l,iq),min_qMass)
    857         IF (q(ij,l,iq)>min_qParent) THEN
    858           Ratio(ij,l,iq2)=q(ij,l,iq2)/q(ij,l,iq)
     852        masseq(ij, l, iq2) = max(masse(ij, l, iq) * q(ij, l, iq), min_qMass)
     853        IF (q(ij, l, iq)>min_qParent) THEN
     854          Ratio(ij, l, iq2) = q(ij, l, iq2) / q(ij, l, iq)
    859855        else
    860           Ratio(ij,l,iq2)=min_ratio
     856          Ratio(ij, l, iq2) = min_ratio
    861857        endif
    862858      enddo
    … …  
    864860  enddo
    865861
    866   do ifils=1,tracers(iq)%nqChildren
    867     iq2=tracers(iq)%iqDescen(ifils)
    868     CALL vlz(Ratio,pente_max,masseq,wq,iq2)
     862  do ifils = 1, tracers(iq)%nqChildren
     863    iq2 = tracers(iq)%iqDescen(ifils)
     864    CALL vlz(Ratio, pente_max, masseq, wq, iq2)
    869865  enddo
    870866  ! end CRisi
    871867
    872   DO l=1,llm
    873      DO ij=1,ip1jmp1
    874         newmasse=masse(ij,l,iq)+w(ij,l+1)-w(ij,l)
    875         q(ij,l,iq)=(q(ij,l,iq)*masse(ij,l,iq)+wq(ij,l+1)-wq(ij,l)) &
    876               /newmasse
    877         masse(ij,l,iq)=newmasse
    878      ENDDO
     868  DO l = 1, llm
     869    DO ij = 1, ip1jmp1
     870      newmasse = masse(ij, l, iq) + w(ij, l + 1) - w(ij, l)
     871      q(ij, l, iq) = (q(ij, l, iq) * masse(ij, l, iq) + wq(ij, l + 1) - wq(ij, l)) &
     872              / newmasse
     873      masse(ij, l, iq) = newmasse
     874    ENDDO
    879875  ENDDO
    880876
    881877  ! retablir les fils en rapport de melange par rapport a l'air:
    882   do ifils=1,tracers(iq)%nqDescen
    883     iq2=tracers(iq)%iqDescen(ifils)
    884     DO l=1,llm
    885       DO ij=1,ip1jmp1
    886         q(ij,l,iq2)=q(ij,l,iq)*Ratio(ij,l,iq2)
     878  do ifils = 1, tracers(iq)%nqDescen
     879    iq2 = tracers(iq)%iqDescen(ifils)
     880    DO l = 1, llm
     881      DO ij = 1, ip1jmp1
     882        q(ij, l, iq2) = q(ij, l, iq) * Ratio(ij, l, iq2)
    887883      enddo
    888884    enddo
    889885  enddo
    890886
    891 
    892887END SUBROUTINE vlz
    893888
    894 SUBROUTINE minmaxq(zq,qmin,qmax,comment)
     889SUBROUTINE minmaxq(zq, qmin, qmax, comment)
    895890
    896891  INCLUDE "dimensions.h"
    897892  INCLUDE "paramet.h"
    898893
    899   CHARACTER(LEN=20) :: comment
    900   REAL :: qmin,qmax
    901   REAL :: zq(ip1jmp1,llm)
    902   REAL :: zzq(iip1,jjp1,llm)
     894  CHARACTER(LEN = 20) :: comment
     895  REAL :: qmin, qmax
     896  REAL :: zq(ip1jmp1, llm)
     897  REAL :: zzq(iip1, jjp1, llm)
    903898
    904899END SUBROUTINE  minmaxq

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/vlspltqs.F90

    r5117 r5119  
    2626  USE comconst_mod, ONLY: cpp
    2727  USE logic_mod, ONLY: adv_qsat_liq
     28  USE lmdz_ssum_scopy, ONLY: scopy
    2829  IMPLICIT NONE
    2930  !
    … …  
    172173  enddo
    173174  !WRITE(*,*) 'vlspltqs 183: fin de la routine'
    174 
    175175
    176176END SUBROUTINE vlspltqs
    … …  
    505505  ! CALL SCOPY((jjm-1)*llm,q(iip1+iip1,1),iip1,q(iip2,1),iip1)
    506506  ! CALL SCOPY((jjm-1)*llm,masse(iip1+iip1,1),iip1,masse(iip2,1),iip1)
    507 
    508507
    509508END SUBROUTINE vlxqs
    … …  
    840839  ! !WRITE(*,*) 'vly 879'
    841840
    842 
    843841END SUBROUTINE vlyqs
Note: See TracChangeset for help on using the changeset viewer.