Ignore:
Timestamp:
Jul 24, 2024, 4:39:59 PM (3 months ago)
Author:
abarral
Message:

Replace iniprint.h by lmdz_iniprint.f90
(lint) along the way

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/dyn3d_common/advn.f90

    r5117 r5118  
    1 
    21! $Header$
    32
    4 SUBROUTINE advn(q,masse,w,pbaru,pbarv,pdt,mode)
     3SUBROUTINE advn(q, masse, w, pbaru, pbarv, pdt, mode)
    54  !
    65  ! Auteur : F. Hourdin
    … …  
    1514  !
    1615  !   --------------------------------------------------------------------
     16  USE lmdz_iniprint, ONLY: lunout, prt_level
    1717  IMPLICIT NONE
    1818  !
    … …  
    2020  include "paramet.h"
    2121  include "comgeom.h"
    22   include "iniprint.h"
    2322
    2423  !
    … …  
    2625  !   ----------
    2726  INTEGER :: mode
    28   REAL :: masse(ip1jmp1,llm)
    29   REAL :: pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm)
    30   REAL :: q(ip1jmp1,llm)
    31   REAL :: w(ip1jmp1,llm),pdt
     27  REAL :: masse(ip1jmp1, llm)
     28  REAL :: pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
     29  REAL :: q(ip1jmp1, llm)
     30  REAL :: w(ip1jmp1, llm), pdt
    3231  !
    3332  !  Local
    3433  !   ---------
    3534  !
    36   INTEGER :: i,ij,l,j,ii
    37   INTEGER :: ijlqmin,iqmin,jqmin,lqmin
     35  INTEGER :: i, ij, l, j, ii
     36  INTEGER :: ijlqmin, iqmin, jqmin, lqmin
    3837  INTEGER :: ismin
    3938  !
    40   REAL :: zm(ip1jmp1,llm),newmasse
    41   REAL :: mu(ip1jmp1,llm)
    42   REAL :: mv(ip1jm,llm)
    43   REAL :: mw(ip1jmp1,llm+1)
    44   REAL :: zq(ip1jmp1,llm),zz,qpn,qps
    45   REAL :: zqg(ip1jmp1,llm),zqd(ip1jmp1,llm)
    46   REAL :: zqs(ip1jmp1,llm),zqn(ip1jmp1,llm)
    47   REAL :: zqh(ip1jmp1,llm),zqb(ip1jmp1,llm)
    48   REAL :: temps0,temps1,temps2,temps3
    49   REAL :: ztemps1,ztemps2,ssum
    50   save temps1,temps2,temps3
    51   REAL :: zzpbar,zzw
    52 
    53   REAL :: qmin,qmax
    54   data qmin,qmax/0.,1./
    55   data temps1,temps2,temps3/0.,0.,0./
     39  REAL :: zm(ip1jmp1, llm), newmasse
     40  REAL :: mu(ip1jmp1, llm)
     41  REAL :: mv(ip1jm, llm)
     42  REAL :: mw(ip1jmp1, llm + 1)
     43  REAL :: zq(ip1jmp1, llm), zz, qpn, qps
     44  REAL :: zqg(ip1jmp1, llm), zqd(ip1jmp1, llm)
     45  REAL :: zqs(ip1jmp1, llm), zqn(ip1jmp1, llm)
     46  REAL :: zqh(ip1jmp1, llm), zqb(ip1jmp1, llm)
     47  REAL :: temps0, temps1, temps2, temps3
     48  REAL :: ztemps1, ztemps2, ssum
     49  save temps1, temps2, temps3
     50  REAL :: zzpbar, zzw
     51
     52  REAL :: qmin, qmax
     53  data qmin, qmax/0., 1./
     54  data temps1, temps2, temps3/0., 0., 0./
    5655
    5756  zzpbar = 0.5 * pdt
    58   zzw    = pdt
    59 
    60   DO l=1,llm
    61     DO ij = iip2,ip1jm
    62         mu(ij,l)=pbaru(ij,l) * zzpbar
    63      ENDDO
    64      DO ij=1,ip1jm
    65         mv(ij,l)=pbarv(ij,l) * zzpbar
    66      ENDDO
    67      DO ij=1,ip1jmp1
    68         mw(ij,l)=w(ij,l) * zzw
    69      ENDDO
     57  zzw = pdt
     58
     59  DO l = 1, llm
     60    DO ij = iip2, ip1jm
     61      mu(ij, l) = pbaru(ij, l) * zzpbar
     62    ENDDO
     63    DO ij = 1, ip1jm
     64      mv(ij, l) = pbarv(ij, l) * zzpbar
     65    ENDDO
     66    DO ij = 1, ip1jmp1
     67      mw(ij, l) = w(ij, l) * zzw
     68    ENDDO
    7069  ENDDO
    7170
    72   DO ij=1,ip1jmp1
    73      mw(ij,llm+1)=0.
     71  DO ij = 1, ip1jmp1
     72    mw(ij, llm + 1) = 0.
    7473  ENDDO
    7574
    76   do l=1,llm
    77      qpn=0.
    78      qps=0.
    79      do ij=1,iim
    80         qpn=qpn+q(ij,l)*masse(ij,l)
    81         qps=qps+q(ip1jm+ij,l)*masse(ip1jm+ij,l)
    82      enddo
    83      qpn=qpn/ssum(iim,masse(1,l),1)
    84      qps=qps/ssum(iim,masse(ip1jm+1,l),1)
    85      do ij=1,iip1
    86         q(ij,l)=qpn
    87         q(ip1jm+ij,l)=qps
    88      enddo
    89   enddo
    90 
    91   do ij=1,ip1jmp1
    92      mw(ij,llm+1)=0.
    93   enddo
    94   do l=1,llm
    95      do ij=1,ip1jmp1
    96         zq(ij,l)=q(ij,l)
    97         zm(ij,l)=masse(ij,l)
    98      enddo
     75  do l = 1, llm
     76    qpn = 0.
     77    qps = 0.
     78    do ij = 1, iim
     79      qpn = qpn + q(ij, l) * masse(ij, l)
     80      qps = qps + q(ip1jm + ij, l) * masse(ip1jm + ij, l)
     81    enddo
     82    qpn = qpn / ssum(iim, masse(1, l), 1)
     83    qps = qps / ssum(iim, masse(ip1jm + 1, l), 1)
     84    do ij = 1, iip1
     85      q(ij, l) = qpn
     86      q(ip1jm + ij, l) = qps
     87    enddo
     88  enddo
     89
     90  do ij = 1, ip1jmp1
     91    mw(ij, llm + 1) = 0.
     92  enddo
     93  do l = 1, llm
     94    do ij = 1, ip1jmp1
     95      zq(ij, l) = q(ij, l)
     96      zm(ij, l) = masse(ij, l)
     97    enddo
    9998  enddo
    10099
    101100  ! CALL minmaxq(zq,qmin,qmax,'avant vlx     ')
    102   CALL advnqx(zq,zqg,zqd)
    103   CALL advnx(zq,zqg,zqd,zm,mu,mode)
    104   CALL advnqy(zq,zqs,zqn)
    105   CALL advny(zq,zqs,zqn,zm,mv)
    106   CALL advnqz(zq,zqh,zqb)
    107   CALL advnz(zq,zqh,zqb,zm,mw)
     101  CALL advnqx(zq, zqg, zqd)
     102  CALL advnx(zq, zqg, zqd, zm, mu, mode)
     103  CALL advnqy(zq, zqs, zqn)
     104  CALL advny(zq, zqs, zqn, zm, mv)
     105  CALL advnqz(zq, zqh, zqb)
     106  CALL advnz(zq, zqh, zqb, zm, mw)
    108107  ! CALL vlz(zq,0.,zm,mw)
    109   CALL advnqy(zq,zqs,zqn)
    110   CALL advny(zq,zqs,zqn,zm,mv)
    111   CALL advnqx(zq,zqg,zqd)
    112   CALL advnx(zq,zqg,zqd,zm,mu,mode)
     108  CALL advnqy(zq, zqs, zqn)
     109  CALL advny(zq, zqs, zqn, zm, mv)
     110  CALL advnqx(zq, zqg, zqd)
     111  CALL advnx(zq, zqg, zqd, zm, mu, mode)
    113112  ! CALL minmaxq(zq,qmin,qmax,'apres vlx     ')
    114113
    115   do l=1,llm
    116      do ij=1,ip1jmp1
    117        q(ij,l)=zq(ij,l)
    118      enddo
    119      do ij=1,ip1jm+1,iip1
    120         q(ij+iim,l)=q(ij,l)
    121      enddo
     114  do l = 1, llm
     115    do ij = 1, ip1jmp1
     116      q(ij, l) = zq(ij, l)
     117    enddo
     118    do ij = 1, ip1jm + 1, iip1
     119      q(ij + iim, l) = q(ij, l)
     120    enddo
    122121  enddo
    123122
    … …  
    125124END SUBROUTINE advn
    126125
    127 SUBROUTINE advnqx(q,qg,qd)
     126SUBROUTINE advnqx(q, qg, qd)
    128127  !
    129128  ! Auteurs:   Calcul des valeurs de q aux point u.
    130129  !
    131130  !   --------------------------------------------------------------------
     131  USE lmdz_iniprint, ONLY: lunout, prt_level
    132132  IMPLICIT NONE
    133133  !
    134134  INCLUDE "dimensions.h"
    135135  INCLUDE "paramet.h"
    136   INCLUDE "iniprint.h"
    137136  !
    138137  !
    139138  !   Arguments:
    140139  !   ----------
    141   REAL :: q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)
     140  REAL :: q(ip1jmp1, llm), qg(ip1jmp1, llm), qd(ip1jmp1, llm)
    142141  !
    143142  !  Local
    144143  !   ---------
    145144  !
    146   INTEGER :: ij,l
    147   !
    148   REAL :: dxqu(ip1jmp1),zqu(ip1jmp1)
    149   REAL :: zqmax(ip1jmp1),zqmin(ip1jmp1)
     145  INTEGER :: ij, l
     146  !
     147  REAL :: dxqu(ip1jmp1), zqu(ip1jmp1)
     148  REAL :: zqmax(ip1jmp1), zqmin(ip1jmp1)
    150149  LOGICAL :: extremum(ip1jmp1)
    151150
    … …  
    157156  !   -----------------------
    158157  IF (mode==0) THEN
    159      do l=1,llm
    160         do ij=1,ip1jm
    161            qd(ij,l)=q(ij,l)
    162            qg(ij,l)=q(ij,l)
    163         enddo
    164      enddo
     158    do l = 1, llm
     159      do ij = 1, ip1jm
     160        qd(ij, l) = q(ij, l)
     161        qg(ij, l) = q(ij, l)
     162      enddo
     163    enddo
    165164  else
    166   do l = 1, llm
    167      do ij=iip2,ip1jm-1
    168         dxqu(ij)=q(ij+1,l)-q(ij,l)
    169         zqu(ij)=0.5*(q(ij+1,l)+q(ij,l))
    170      enddo
    171      do ij=iip1+iip1,ip1jm,iip1
    172         dxqu(ij)=dxqu(ij-iim)
    173         zqu(ij)=zqu(ij-iim)
    174      enddo
    175      do ij=iip2,ip1jm-1
    176         zqu(ij)=zqu(ij)-dxqu(ij+1)/12.
    177      enddo
    178      do ij=iip1+iip1,ip1jm,iip1
    179         zqu(ij)=zqu(ij-iim)
    180      enddo
    181      do ij=iip2+1,ip1jm
    182         zqu(ij)=zqu(ij)+dxqu(ij-1)/12.
    183      enddo
    184      do ij=iip1+iip1,ip1jm,iip1
    185         zqu(ij-iim)=zqu(ij)
    186      enddo
    187 
    188   !   calcul des valeurs max et min acceptees aux interfaces
    189 
    190      do ij=iip2,ip1jm-1
    191         zqmax(ij)=max(q(ij+1,l),q(ij,l))
    192         zqmin(ij)=min(q(ij+1,l),q(ij,l))
    193      enddo
    194      do ij=iip1+iip1,ip1jm,iip1
    195         zqmax(ij)=zqmax(ij-iim)
    196         zqmin(ij)=zqmin(ij-iim)
    197      enddo
    198      do ij=iip2+1,ip1jm
    199         extremum(ij)=dxqu(ij)*dxqu(ij-1)<=0.
    200      enddo
    201      do ij=iip1+iip1,ip1jm,iip1
    202         extremum(ij-iim)=extremum(ij)
    203      enddo
    204      do ij=iip2,ip1jm
    205         zqu(ij)=min(max(zqmin(ij),zqu(ij)),zqmax(ij))
    206      enddo
    207      do ij=iip2+1,ip1jm
     165    do l = 1, llm
     166      do ij = iip2, ip1jm - 1
     167        dxqu(ij) = q(ij + 1, l) - q(ij, l)
     168        zqu(ij) = 0.5 * (q(ij + 1, l) + q(ij, l))
     169      enddo
     170      do ij = iip1 + iip1, ip1jm, iip1
     171        dxqu(ij) = dxqu(ij - iim)
     172        zqu(ij) = zqu(ij - iim)
     173      enddo
     174      do ij = iip2, ip1jm - 1
     175        zqu(ij) = zqu(ij) - dxqu(ij + 1) / 12.
     176      enddo
     177      do ij = iip1 + iip1, ip1jm, iip1
     178        zqu(ij) = zqu(ij - iim)
     179      enddo
     180      do ij = iip2 + 1, ip1jm
     181        zqu(ij) = zqu(ij) + dxqu(ij - 1) / 12.
     182      enddo
     183      do ij = iip1 + iip1, ip1jm, iip1
     184        zqu(ij - iim) = zqu(ij)
     185      enddo
     186
     187      !   calcul des valeurs max et min acceptees aux interfaces
     188
     189      do ij = iip2, ip1jm - 1
     190        zqmax(ij) = max(q(ij + 1, l), q(ij, l))
     191        zqmin(ij) = min(q(ij + 1, l), q(ij, l))
     192      enddo
     193      do ij = iip1 + iip1, ip1jm, iip1
     194        zqmax(ij) = zqmax(ij - iim)
     195        zqmin(ij) = zqmin(ij - iim)
     196      enddo
     197      do ij = iip2 + 1, ip1jm
     198        extremum(ij) = dxqu(ij) * dxqu(ij - 1)<=0.
     199      enddo
     200      do ij = iip1 + iip1, ip1jm, iip1
     201        extremum(ij - iim) = extremum(ij)
     202      enddo
     203      do ij = iip2, ip1jm
     204        zqu(ij) = min(max(zqmin(ij), zqu(ij)), zqmax(ij))
     205      enddo
     206      do ij = iip2 + 1, ip1jm
    208207        IF(extremum(ij)) THEN
    209            qg(ij,l)=q(ij,l)
    210            qd(ij,l)=q(ij,l)
     208          qg(ij, l) = q(ij, l)
     209          qd(ij, l) = q(ij, l)
    211210        else
    212            qd(ij,l)=zqu(ij)
    213            qg(ij,l)=zqu(ij-1)
     211          qd(ij, l) = zqu(ij)
     212          qg(ij, l) = zqu(ij - 1)
    214213        endif
    215      enddo
    216      do ij=iip1+iip1,ip1jm,iip1
    217         qd(ij-iim,l)=qd(ij,l)
    218         qg(ij-iim,l)=qg(ij,l)
    219      enddo
    220 
    221      goto 8888
    222 
    223      do ij=iip2+1,ip1jm
    224         IF(extremum(ij).and..not.extremum(ij-1)) &
    225               qd(ij-1,l)=q(ij,l)
    226      enddo
    227 
    228      do ij=iip1+iip1,ip1jm,iip1
    229         qd(ij-iim,l)=qd(ij,l)
    230      enddo
    231      do ij=iip2,ip1jm-1
    232         IF (extremum(ij).and..not.extremum(ij+1)) &
    233               qg(ij+1,l)=q(ij,l)
    234      enddo
    235 
    236      do ij=iip1+iip1,ip1jm,iip1
    237         qg(ij,l)=qg(ij-iim,l)
    238      enddo
    239 8888   continue
    240   enddo
     214      enddo
     215      do ij = iip1 + iip1, ip1jm, iip1
     216        qd(ij - iim, l) = qd(ij, l)
     217        qg(ij - iim, l) = qg(ij, l)
     218      enddo
     219
     220      goto 8888
     221
     222      do ij = iip2 + 1, ip1jm
     223        IF(extremum(ij).and..not.extremum(ij - 1)) &
     224                qd(ij - 1, l) = q(ij, l)
     225      enddo
     226
     227      do ij = iip1 + iip1, ip1jm, iip1
     228        qd(ij - iim, l) = qd(ij, l)
     229      enddo
     230      do ij = iip2, ip1jm - 1
     231        IF (extremum(ij).and..not.extremum(ij + 1)) &
     232                qg(ij + 1, l) = q(ij, l)
     233      enddo
     234
     235      do ij = iip1 + iip1, ip1jm, iip1
     236        qg(ij, l) = qg(ij - iim, l)
     237      enddo
     238      8888   continue
     239    enddo
    241240  ENDIF
    242241  RETURN
    243242END SUBROUTINE advnqx
    244 SUBROUTINE advnqy(q,qs,qn)
     243SUBROUTINE advnqy(q, qs, qn)
    245244  !
    246245  ! Auteurs:   Calcul des valeurs de q aux point v.
    247246  !
    248247  !   --------------------------------------------------------------------
     248  USE lmdz_iniprint, ONLY: lunout, prt_level
    249249  IMPLICIT NONE
    250250  !
    251251  INCLUDE "dimensions.h"
    252252  INCLUDE "paramet.h"
    253   INCLUDE "iniprint.h"
    254253  !
    255254  !
    256255  !   Arguments:
    257256  !   ----------
    258   REAL :: q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)
     257  REAL :: q(ip1jmp1, llm), qs(ip1jmp1, llm), qn(ip1jmp1, llm)
    259258  !
    260259  !  Local
    261260  !   ---------
    262261  !
    263   INTEGER :: ij,l
    264   !
    265   REAL :: dyqv(ip1jm),zqv(ip1jm,llm)
    266   REAL :: zqmax(ip1jm),zqmin(ip1jm)
     262  INTEGER :: ij, l
     263  !
     264  REAL :: dyqv(ip1jm), zqv(ip1jm, llm)
     265  REAL :: zqmax(ip1jm), zqmin(ip1jm)
    267266  LOGICAL :: extremum(ip1jmp1)
    268267
    … …  
    272271
    273272  IF (mode==0) THEN
    274      do l=1,llm
    275         do ij=1,ip1jmp1
    276            qn(ij,l)=q(ij,l)
    277            qs(ij,l)=q(ij,l)
    278         enddo
    279      enddo
     273    do l = 1, llm
     274      do ij = 1, ip1jmp1
     275        qn(ij, l) = q(ij, l)
     276        qs(ij, l) = q(ij, l)
     277      enddo
     278    enddo
    280279  else
    281280
    282   !   calcul des pentes en u:
    283   !   -----------------------
    284   do l = 1, llm
    285      do ij=1,ip1jm
    286         dyqv(ij)=q(ij,l)-q(ij+iip1,l)
    287      enddo
    288 
    289      do ij=iip2,ip1jm-iip1
    290         zqv(ij,l)=0.5*(q(ij+iip1,l)+q(ij,l))
    291         zqv(ij,l)=zqv(ij,l)+(dyqv(ij+iip1)-dyqv(ij-iip1))/12.
    292      enddo
    293 
    294      do ij=iip2,ip1jm
    295         extremum(ij)=dyqv(ij)*dyqv(ij-iip1)<=0.
    296      enddo
    297 
    298   ! Pas de pentes aux poles
    299      do ij=1,iip1
    300         zqv(ij,l)=q(ij,l)
    301         zqv(ip1jm-iip1+ij,l)=q(ip1jm+ij,l)
    302         extremum(ij)=.TRUE.
    303         extremum(ip1jmp1-iip1+ij)=.TRUE.
    304      enddo
    305 
    306   !   calcul des valeurs max et min acceptees aux interfaces
    307      do ij=1,ip1jm
    308         zqmax(ij)=max(q(ij+iip1,l),q(ij,l))
    309         zqmin(ij)=min(q(ij+iip1,l),q(ij,l))
    310      enddo
    311 
    312      do ij=1,ip1jm
    313         zqv(ij,l)=min(max(zqmin(ij),zqv(ij,l)),zqmax(ij))
    314      enddo
    315 
    316      do ij=iip2,ip1jm
     281    !   calcul des pentes en u:
     282    !   -----------------------
     283    do l = 1, llm
     284      do ij = 1, ip1jm
     285        dyqv(ij) = q(ij, l) - q(ij + iip1, l)
     286      enddo
     287
     288      do ij = iip2, ip1jm - iip1
     289        zqv(ij, l) = 0.5 * (q(ij + iip1, l) + q(ij, l))
     290        zqv(ij, l) = zqv(ij, l) + (dyqv(ij + iip1) - dyqv(ij - iip1)) / 12.
     291      enddo
     292
     293      do ij = iip2, ip1jm
     294        extremum(ij) = dyqv(ij) * dyqv(ij - iip1)<=0.
     295      enddo
     296
     297      ! Pas de pentes aux poles
     298      do ij = 1, iip1
     299        zqv(ij, l) = q(ij, l)
     300        zqv(ip1jm - iip1 + ij, l) = q(ip1jm + ij, l)
     301        extremum(ij) = .TRUE.
     302        extremum(ip1jmp1 - iip1 + ij) = .TRUE.
     303      enddo
     304
     305      !   calcul des valeurs max et min acceptees aux interfaces
     306      do ij = 1, ip1jm
     307        zqmax(ij) = max(q(ij + iip1, l), q(ij, l))
     308        zqmin(ij) = min(q(ij + iip1, l), q(ij, l))
     309      enddo
     310
     311      do ij = 1, ip1jm
     312        zqv(ij, l) = min(max(zqmin(ij), zqv(ij, l)), zqmax(ij))
     313      enddo
     314
     315      do ij = iip2, ip1jm
    317316        IF(extremum(ij)) THEN
    318            qs(ij,l)=q(ij,l)
    319            qn(ij,l)=q(ij,l)
    320            ! if (.NOT.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
    321            ! if (.NOT.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l)
     317          qs(ij, l) = q(ij, l)
     318          qn(ij, l) = q(ij, l)
     319          ! if (.NOT.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
     320          ! if (.NOT.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l)
    322321        else
    323            qs(ij,l)=zqv(ij,l)
    324            qn(ij,l)=zqv(ij-iip1,l)
     322          qs(ij, l) = zqv(ij, l)
     323          qn(ij, l) = zqv(ij - iip1, l)
    325324        endif
    326      enddo
    327 
    328      do ij=1,iip1
    329         qs(ij,l)=q(ij,l)
    330         qn(ij,l)=q(ij,l)
    331         qs(ip1jm+ij,l)=q(ip1jm+ij,l)
    332         qn(ip1jm+ij,l)=q(ip1jm+ij,l)
    333      enddo
    334 
    335   enddo
     325      enddo
     326
     327      do ij = 1, iip1
     328        qs(ij, l) = q(ij, l)
     329        qn(ij, l) = q(ij, l)
     330        qs(ip1jm + ij, l) = q(ip1jm + ij, l)
     331        qn(ip1jm + ij, l) = q(ip1jm + ij, l)
     332      enddo
     333
     334    enddo
    336335  ENDIF
    337336  RETURN
    338337END SUBROUTINE advnqy
    339338
    340 SUBROUTINE advnqz(q,qh,qb)
     339SUBROUTINE advnqz(q, qh, qb)
    341340  !
    342341  ! Auteurs:   Calcul des valeurs de q aux point v.
    343342  !
    344343  !   --------------------------------------------------------------------
     344  USE lmdz_iniprint, ONLY: lunout, prt_level
    345345  IMPLICIT NONE
    346346  !
    347347  INCLUDE "dimensions.h"
    348348  INCLUDE "paramet.h"
    349   INCLUDE "iniprint.h"
    350349  !
    351350  !
    352351  !   Arguments:
    353352  !   ----------
    354   REAL :: q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)
     353  REAL :: q(ip1jmp1, llm), qh(ip1jmp1, llm), qb(ip1jmp1, llm)
    355354  !
    356355  !  Local
    357356  !   ---------
    358357  !
    359   INTEGER :: ij,l
    360   !
    361   REAL :: dzqw(ip1jmp1,llm+1),zqw(ip1jmp1,llm+1)
    362   REAL :: zqmax(ip1jmp1,llm),zqmin(ip1jmp1,llm)
    363   LOGICAL :: extremum(ip1jmp1,llm)
     358  INTEGER :: ij, l
     359  !
     360  REAL :: dzqw(ip1jmp1, llm + 1), zqw(ip1jmp1, llm + 1)
     361  REAL :: zqmax(ip1jmp1, llm), zqmin(ip1jmp1, llm)
     362  LOGICAL :: extremum(ip1jmp1, llm)
    364363
    365364  INTEGER :: mode
    … …  
    372371
    373372  IF (mode==0) THEN
    374      do l=1,llm
    375         do ij=1,ip1jmp1
    376            qb(ij,l)=q(ij,l)
    377            qh(ij,l)=q(ij,l)
    378         enddo
    379      enddo
     373    do l = 1, llm
     374      do ij = 1, ip1jmp1
     375        qb(ij, l) = q(ij, l)
     376        qh(ij, l) = q(ij, l)
     377      enddo
     378    enddo
    380379  else
    381   do l = 2, llm
    382      do ij=1,ip1jmp1
    383         dzqw(ij,l)=q(ij,l-1)-q(ij,l)
    384         zqw(ij,l)=0.5*(q(ij,l-1)+q(ij,l))
    385      enddo
    386   enddo
    387   do ij=1,ip1jmp1
    388      dzqw(ij,1)=0.
    389      dzqw(ij,llm+1)=0.
    390   enddo
    391   do l=2,llm
    392      do ij=1,ip1jmp1
    393         zqw(ij,l)=zqw(ij,l)+(dzqw(ij,l+1)-dzqw(ij,l-1))/12.
    394      enddo
    395   enddo
    396   do l=2,llm-1
    397      do ij=1,ip1jmp1
    398         extremum(ij,l)=dzqw(ij,l)*dzqw(ij,l+1)<=0.
    399      enddo
    400   enddo
    401 
    402   ! Pas de pentes en bas et en haut
    403      do ij=1,ip1jmp1
    404         zqw(ij,2)=q(ij,1)
    405         zqw(ij,llm)=q(ij,llm)
    406         extremum(ij,1)=.TRUE.
    407         extremum(ij,llm)=.TRUE.
    408      enddo
    409 
    410   !   calcul des valeurs max et min acceptees aux interfaces
    411   do l=2,llm
    412      do ij=1,ip1jmp1
    413         zqmax(ij,l)=max(q(ij,l-1),q(ij,l))
    414         zqmin(ij,l)=min(q(ij,l-1),q(ij,l))
    415      enddo
    416   enddo
    417 
    418   do l=2,llm
    419      do ij=1,ip1jmp1
    420         zqw(ij,l)=min(max(zqmin(ij,l),zqw(ij,l)),zqmax(ij,l))
    421      enddo
    422   enddo
    423 
    424   do l=2,llm-1
    425      do ij=1,ip1jmp1
    426         IF(extremum(ij,l)) THEN
    427            qh(ij,l)=q(ij,l)
    428            qb(ij,l)=q(ij,l)
     380    do l = 2, llm
     381      do ij = 1, ip1jmp1
     382        dzqw(ij, l) = q(ij, l - 1) - q(ij, l)
     383        zqw(ij, l) = 0.5 * (q(ij, l - 1) + q(ij, l))
     384      enddo
     385    enddo
     386    do ij = 1, ip1jmp1
     387      dzqw(ij, 1) = 0.
     388      dzqw(ij, llm + 1) = 0.
     389    enddo
     390    do l = 2, llm
     391      do ij = 1, ip1jmp1
     392        zqw(ij, l) = zqw(ij, l) + (dzqw(ij, l + 1) - dzqw(ij, l - 1)) / 12.
     393      enddo
     394    enddo
     395    do l = 2, llm - 1
     396      do ij = 1, ip1jmp1
     397        extremum(ij, l) = dzqw(ij, l) * dzqw(ij, l + 1)<=0.
     398      enddo
     399    enddo
     400
     401    ! Pas de pentes en bas et en haut
     402    do ij = 1, ip1jmp1
     403      zqw(ij, 2) = q(ij, 1)
     404      zqw(ij, llm) = q(ij, llm)
     405      extremum(ij, 1) = .TRUE.
     406      extremum(ij, llm) = .TRUE.
     407    enddo
     408
     409    !   calcul des valeurs max et min acceptees aux interfaces
     410    do l = 2, llm
     411      do ij = 1, ip1jmp1
     412        zqmax(ij, l) = max(q(ij, l - 1), q(ij, l))
     413        zqmin(ij, l) = min(q(ij, l - 1), q(ij, l))
     414      enddo
     415    enddo
     416
     417    do l = 2, llm
     418      do ij = 1, ip1jmp1
     419        zqw(ij, l) = min(max(zqmin(ij, l), zqw(ij, l)), zqmax(ij, l))
     420      enddo
     421    enddo
     422
     423    do l = 2, llm - 1
     424      do ij = 1, ip1jmp1
     425        IF(extremum(ij, l)) THEN
     426          qh(ij, l) = q(ij, l)
     427          qb(ij, l) = q(ij, l)
    429428        else
    430            qh(ij,l)=zqw(ij,l+1)
    431            qb(ij,l)=zqw(ij,l)
     429          qh(ij, l) = zqw(ij, l + 1)
     430          qb(ij, l) = zqw(ij, l)
    432431        endif
    433      enddo
    434   enddo
    435   ! do l=2,llm-1
    436   !    do ij=1,ip1jmp1
    437   !       IF(extremum(ij,l)) THEN
    438   !          if (.NOT.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
    439   !          if (.NOT.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
    440   !       endif
    441   !    enddo
    442   ! enddo
    443 
    444   do ij=1,ip1jmp1
    445      qb(ij,1)=q(ij,1)
    446      qh(ij,1)=q(ij,1)
    447      qb(ij,llm)=q(ij,llm)
    448      qh(ij,llm)=q(ij,llm)
    449   enddo
     432      enddo
     433    enddo
     434    ! do l=2,llm-1
     435    !    do ij=1,ip1jmp1
     436    !       IF(extremum(ij,l)) THEN
     437    !          if (.NOT.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
     438    !          if (.NOT.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
     439    !       endif
     440    !    enddo
     441    ! enddo
     442
     443    do ij = 1, ip1jmp1
     444      qb(ij, 1) = q(ij, 1)
     445      qh(ij, 1) = q(ij, 1)
     446      qb(ij, llm) = q(ij, llm)
     447      qh(ij, llm) = q(ij, llm)
     448    enddo
    450449
    451450  ENDIF
    … …  
    454453END SUBROUTINE advnqz
    455454
    456 SUBROUTINE advnx(q,qg,qd,masse,u_m,mode)
     455SUBROUTINE advnx(q, qg, qd, masse, u_m, mode)
    457456  !
    458457  ! Auteur : F. Hourdin
    … …  
    465464  !
    466465  !   --------------------------------------------------------------------
     466  USE lmdz_iniprint, ONLY: lunout, prt_level
    467467  IMPLICIT NONE
    468468  !
    469469  include "dimensions.h"
    470470  include "paramet.h"
    471   include "iniprint.h"
    472471  !
    473472  !
    … …  
    475474  !   ----------
    476475  INTEGER :: mode
    477   REAL :: masse(ip1jmp1,llm)
    478   REAL :: u_m( ip1jmp1,llm )
    479   REAL :: q(ip1jmp1,llm),qd(ip1jmp1,llm),qg(ip1jmp1,llm)
     476  REAL :: masse(ip1jmp1, llm)
     477  REAL :: u_m(ip1jmp1, llm)
     478  REAL :: q(ip1jmp1, llm), qd(ip1jmp1, llm), qg(ip1jmp1, llm)
    480479  !
    481480  !  Local
    482481  !   ---------
    483482  !
    484   INTEGER :: i,j,ij,l,indu(ip1jmp1),niju,iju,ijq
    485   INTEGER :: n0,nl(llm)
    486   !
    487   REAL :: new_m,zu_m,zdq,zz
    488   REAL :: zsigg(ip1jmp1,llm),zsigd(ip1jmp1,llm),zsig
    489   REAL :: u_mq(ip1jmp1,llm)
    490 
    491   REAL :: zm,zq,zsigm,zsigp,zqm,zqp,zu
    492 
    493   LOGICAL :: ladvplus(ip1jmp1,llm)
     483  INTEGER :: i, j, ij, l, indu(ip1jmp1), niju, iju, ijq
     484  INTEGER :: n0, nl(llm)
     485  !
     486  REAL :: new_m, zu_m, zdq, zz
     487  REAL :: zsigg(ip1jmp1, llm), zsigd(ip1jmp1, llm), zsig
     488  REAL :: u_mq(ip1jmp1, llm)
     489
     490  REAL :: zm, zq, zsigm, zsigp, zqm, zqp, zu
     491
     492  LOGICAL :: ladvplus(ip1jmp1, llm)
    494493
    495494  REAL :: prec
    … …  
    498497  data prec/1.e-15/
    499498
    500   do l=1,llm
    501         do ij=iip2,ip1jm
    502            zdq=qd(ij,l)-qg(ij,l)
    503            ! if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) THEN
    504            !    PRINT*,'probleme au point ij=',ij,'  l=',l
    505            !    PRINT*,qd(ij,l),q(ij,l),qg(ij,l)
    506            !    qd(ij,l)=q(ij,l)
    507            !    qg(ij,l)=q(ij,l)
    508            ! END IF
    509            IF(abs(zdq)>prec) THEN
    510               zsigd(ij,l)=(q(ij,l)-qg(ij,l))/zdq
    511               zsigg(ij,l)=1.-zsigd(ij,l)
    512               ! IF(.NOT.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .AND.
    513   !    s               zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) THEN
    514               !    PRINT*,'probleme au point ij=',ij,'  l=',l
    515               !    PRINT*,'sigg=',zsigg(ij,l),'  sigd=',zsigd(ij,l)
    516               !    PRINT*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq
    517               !    stop
    518               ! END IF
    519            else
    520               zsigd(ij,l)=0.5
    521               zsigg(ij,l)=0.5
    522               qd(ij,l)=q(ij,l)
    523               qg(ij,l)=q(ij,l)
    524            endif
    525         enddo
    526    enddo
     499  do l = 1, llm
     500    do ij = iip2, ip1jm
     501      zdq = qd(ij, l) - qg(ij, l)
     502      ! if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) THEN
     503      !    PRINT*,'probleme au point ij=',ij,'  l=',l
     504      !    PRINT*,qd(ij,l),q(ij,l),qg(ij,l)
     505      !    qd(ij,l)=q(ij,l)
     506      !    qg(ij,l)=q(ij,l)
     507      ! END IF
     508      IF(abs(zdq)>prec) THEN
     509        zsigd(ij, l) = (q(ij, l) - qg(ij, l)) / zdq
     510        zsigg(ij, l) = 1. - zsigd(ij, l)
     511        ! IF(.NOT.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .AND.
     512        !    s               zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) THEN
     513        !    PRINT*,'probleme au point ij=',ij,'  l=',l
     514        !    PRINT*,'sigg=',zsigg(ij,l),'  sigd=',zsigd(ij,l)
     515        !    PRINT*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq
     516        !    stop
     517        ! END IF
     518      else
     519        zsigd(ij, l) = 0.5
     520        zsigg(ij, l) = 0.5
     521        qd(ij, l) = q(ij, l)
     522        qg(ij, l) = q(ij, l)
     523      endif
     524    enddo
     525  enddo
    527526
    528527  !   calcul de la pente maximum dans la maille en valeur absolue
    529528
    530    do l=1,llm
    531    do ij=iip2,ip1jm-1
    532       IF (u_m(ij,l)>=0.) THEN
    533          zsigp=zsigd(ij,l)
    534          zsigm=zsigg(ij,l)
    535          zqp=qd(ij,l)
    536          zqm=qg(ij,l)
    537          zm=masse(ij,l)
    538          zq=q(ij,l)
     529  do l = 1, llm
     530    do ij = iip2, ip1jm - 1
     531      IF (u_m(ij, l)>=0.) THEN
     532        zsigp = zsigd(ij, l)
     533        zsigm = zsigg(ij, l)
     534        zqp = qd(ij, l)
     535        zqm = qg(ij, l)
     536        zm = masse(ij, l)
     537        zq = q(ij, l)
    539538      else
    540          zsigm=zsigd(ij+1,l)
    541          zsigp=zsigg(ij+1,l)
    542          zqm=qd(ij+1,l)
    543          zqp=qg(ij+1,l)
    544          zm=masse(ij+1,l)
    545          zq=q(ij+1,l)
    546       endif
    547       zu=abs(u_m(ij,l))
    548       ladvplus(ij,l)=zu>zm
    549       zsig=zu/zm
    550       IF(zsig==0.) zsigp=0.1
     539        zsigm = zsigd(ij + 1, l)
     540        zsigp = zsigg(ij + 1, l)
     541        zqm = qd(ij + 1, l)
     542        zqp = qg(ij + 1, l)
     543        zm = masse(ij + 1, l)
     544        zq = q(ij + 1, l)
     545      endif
     546      zu = abs(u_m(ij, l))
     547      ladvplus(ij, l) = zu>zm
     548      zsig = zu / zm
     549      IF(zsig==0.) zsigp = 0.1
    551550      IF (mode==1) THEN
    552          IF (zsig<=zsigp) THEN
    553              u_mq(ij,l)=u_m(ij,l)*zqp
    554          ELSE IF (mode==1) THEN
    555              u_mq(ij,l)= &
    556                    sign(zm,u_m(ij,l))*(zsigp*zqp+(zsig-zsigp)*zqm)
    557          endif
     551        IF (zsig<=zsigp) THEN
     552          u_mq(ij, l) = u_m(ij, l) * zqp
     553        ELSE IF (mode==1) THEN
     554          u_mq(ij, l) = &
     555                  sign(zm, u_m(ij, l)) * (zsigp * zqp + (zsig - zsigp) * zqm)
     556        endif
    558557      else
    559          IF (zsig<=zsigp) THEN
    560              u_mq(ij,l)=u_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
    561          else
    562             zz=0.5*(zsig-zsigp)/zsigm
    563             u_mq(ij,l)=sign(zm,u_m(ij,l))*( 0.5*(zq+zqp)*zsigp &
    564                   +(zsig-zsigp)*(zq+zz*(zqm-zq)) )
    565          endif
     558        IF (zsig<=zsigp) THEN
     559          u_mq(ij, l) = u_m(ij, l) * (zqp - 0.5 * zsig / zsigp * (zqp - zq))
     560        else
     561          zz = 0.5 * (zsig - zsigp) / zsigm
     562          u_mq(ij, l) = sign(zm, u_m(ij, l)) * (0.5 * (zq + zqp) * zsigp &
     563                  + (zsig - zsigp) * (zq + zz * (zqm - zq)))
     564        endif
    566565      endif
    567566      ! IF(zsig.lt.0.) THEN
    … …  
    569568      !    stop
    570569      ! END IF
    571   enddo
    572   enddo
    573 
    574   do l=1,llm
    575    do ij=iip1+iip1,ip1jm,iip1
    576       u_mq(ij,l)=u_mq(ij-iim,l)
    577       ladvplus(ij,l)=ladvplus(ij-iim,l)
    578    enddo
     570    enddo
     571  enddo
     572
     573  do l = 1, llm
     574    do ij = iip1 + iip1, ip1jm, iip1
     575      u_mq(ij, l) = u_mq(ij - iim, l)
     576      ladvplus(ij, l) = ladvplus(ij - iim, l)
     577    enddo
    579578  enddo
    580579
    … …  
    583582  !=================================================================
    584583  !   tris des regions a traiter
    585   n0=0
    586   do l=1,llm
    587      nl(l)=0
    588      do ij=iip2,ip1jm
    589         IF(ladvplus(ij,l)) THEN
    590            nl(l)=nl(l)+1
    591            u_mq(ij,l)=0.
    592         endif
    593      enddo
    594      n0=n0+nl(l)
     584  n0 = 0
     585  do l = 1, llm
     586    nl(l) = 0
     587    do ij = iip2, ip1jm
     588      IF(ladvplus(ij, l)) THEN
     589        nl(l) = nl(l) + 1
     590        u_mq(ij, l) = 0.
     591      endif
     592    enddo
     593    n0 = n0 + nl(l)
    595594  enddo
    596595
    597596  IF(n0>1) THEN
    598   IF (prt_level > 9) WRITE(lunout,*) &
    599         'Nombre de points pour lesquels on advect plus que le' &
    600         ,'contenu de la maille : ',n0
    601 
    602      do l=1,llm
    603         IF(nl(l)>0) THEN
    604            iju=0
    605   !   indicage des mailles concernees par le traitement special
    606            do ij=iip2,ip1jm
    607               IF(ladvplus(ij,l).AND.mod(ij,iip1)/=0) THEN
    608                  iju=iju+1
    609                  indu(iju)=ij
     597    IF (prt_level > 9) WRITE(lunout, *) &
     598            'Nombre de points pour lesquels on advect plus que le' &
     599            , 'contenu de la maille : ', n0
     600
     601    do l = 1, llm
     602      IF(nl(l)>0) THEN
     603        iju = 0
     604        !   indicage des mailles concernees par le traitement special
     605        do ij = iip2, ip1jm
     606          IF(ladvplus(ij, l).AND.mod(ij, iip1)/=0) THEN
     607            iju = iju + 1
     608            indu(iju) = ij
     609          endif
     610        enddo
     611        niju = iju
     612        ! PRINT*,'niju,nl',niju,nl(l)
     613
     614        !  traitement des mailles
     615        do iju = 1, niju
     616          ij = indu(iju)
     617          j = (ij - 1) / iip1 + 1
     618          zu_m = u_m(ij, l)
     619          u_mq(ij, l) = 0.
     620          IF(zu_m>0.) THEN
     621            ijq = ij
     622            i = ijq - (j - 1) * iip1
     623            !   accumulation pour les mailles completements advectees
     624            do while(zu_m>masse(ijq, l))
     625              u_mq(ij, l) = u_mq(ij, l) + q(ijq, l) * masse(ijq, l)
     626              zu_m = zu_m - masse(ijq, l)
     627              i = mod(i - 2 + iim, iim) + 1
     628              ijq = (j - 1) * iip1 + i
     629            enddo
     630            !   MODIFS SPECIFIQUES DU SCHEMA
     631            !   ajout de la maille non completement advectee
     632            zsig = zu_m / masse(ijq, l)
     633            IF(zsig<=zsigd(ijq, l)) THEN
     634              u_mq(ij, l) = u_mq(ij, l) + zu_m * (qd(ijq, l) &
     635                      - 0.5 * zsig / zsigd(ijq, l) * (qd(ijq, l) - q(ijq, l)))
     636            else
     637              ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
     638              ! goto 8888
     639              zz = 0.5 * (zsig - zsigd(ijq, l)) / zsigg(ijq, l)
     640              IF(.NOT.(zz>0..and.zz<=0.5)) THEN
     641                WRITE(lunout, *)'probleme2 au point ij=', ij, &
     642                        '  l=', l
     643                WRITE(lunout, *)'zz=', zz
     644                stop
    610645              endif
    611            enddo
    612            niju=iju
    613            ! PRINT*,'niju,nl',niju,nl(l)
    614 
    615   !  traitement des mailles
    616            do iju=1,niju
    617               ij=indu(iju)
    618               j=(ij-1)/iip1+1
    619               zu_m=u_m(ij,l)
    620               u_mq(ij,l)=0.
    621               IF(zu_m>0.) THEN
    622                  ijq=ij
    623                  i=ijq-(j-1)*iip1
    624   !   accumulation pour les mailles completements advectees
    625                  do while(zu_m>masse(ijq,l))
    626                     u_mq(ij,l)=u_mq(ij,l)+q(ijq,l)*masse(ijq,l)
    627                     zu_m=zu_m-masse(ijq,l)
    628                     i=mod(i-2+iim,iim)+1
    629                     ijq=(j-1)*iip1+i
    630                  enddo
    631   !   MODIFS SPECIFIQUES DU SCHEMA
    632   !   ajout de la maille non completement advectee
    633          zsig=zu_m/masse(ijq,l)
    634          IF(zsig<=zsigd(ijq,l)) THEN
    635             u_mq(ij,l)=u_mq(ij,l)+zu_m*(qd(ijq,l) &
    636                   -0.5*zsig/zsigd(ijq,l)*(qd(ijq,l)-q(ijq,l)))
    637          else
    638             ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
    639       ! goto 8888
    640             zz=0.5*(zsig-zsigd(ijq,l))/zsigg(ijq,l)
    641             IF(.NOT.(zz>0..and.zz<=0.5)) THEN
    642                  WRITE(lunout,*)'probleme2 au point ij=',ij, &
    643                        '  l=',l
    644                  WRITE(lunout,*)'zz=',zz
    645                  stop
     646              u_mq(ij, l) = u_mq(ij, l) + masse(ijq, l) * (&
     647                      0.5 * (q(ijq, l) + qd(ijq, l)) * zsigd(ijq, l) &
     648                              + (zsig - zsigd(ijq, l)) * (q(ijq, l) + zz * (qg(ijq, l) - q(ijq, l))))
    646649            endif
    647             u_mq(ij,l)=u_mq(ij,l)+masse(ijq,l)*( &
    648                   0.5*(q(ijq,l)+qd(ijq,l))*zsigd(ijq,l) &
    649                   +(zsig-zsigd(ijq,l))*(q(ijq,l)+zz*(qg(ijq,l)-q(ijq,l))) )
    650          endif
    651               else
    652                  ijq=ij+1
    653                  i=ijq-(j-1)*iip1
    654   !   accumulation pour les mailles completements advectees
    655                  do while(-zu_m>masse(ijq,l))
    656                     u_mq(ij,l)=u_mq(ij,l)-q(ijq,l)*masse(ijq,l)
    657                     zu_m=zu_m+masse(ijq,l)
    658                     i=mod(i,iim)+1
    659                     ijq=(j-1)*iip1+i
    660                  enddo
    661   !   ajout de la maille non completement advectee
    662   ! 2eme MODIF SPECIFIQUE
    663          zsig=-zu_m/masse(ij+1,l)
    664          IF(zsig<=zsigg(ijq,l)) THEN
    665             u_mq(ij,l)=u_mq(ij,l)+zu_m*(qg(ijq,l) &
    666                   -0.5*zsig/zsigg(ijq,l)*(qg(ijq,l)-q(ijq,l)))
    667          else
    668             ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
    669         ! goto 9999
    670             zz=0.5*(zsig-zsigg(ijq,l))/zsigd(ijq,l)
    671             IF(.NOT.(zz>0..and.zz<=0.5)) THEN
    672                  WRITE(lunout,*)'probleme22 au point ij=',ij &
    673                        ,'  l=',l
    674                  WRITE(lunout,*)'zz=',zz
    675                  stop
     650          else
     651            ijq = ij + 1
     652            i = ijq - (j - 1) * iip1
     653            !   accumulation pour les mailles completements advectees
     654            do while(-zu_m>masse(ijq, l))
     655              u_mq(ij, l) = u_mq(ij, l) - q(ijq, l) * masse(ijq, l)
     656              zu_m = zu_m + masse(ijq, l)
     657              i = mod(i, iim) + 1
     658              ijq = (j - 1) * iip1 + i
     659            enddo
     660            !   ajout de la maille non completement advectee
     661            ! 2eme MODIF SPECIFIQUE
     662            zsig = -zu_m / masse(ij + 1, l)
     663            IF(zsig<=zsigg(ijq, l)) THEN
     664              u_mq(ij, l) = u_mq(ij, l) + zu_m * (qg(ijq, l) &
     665                      - 0.5 * zsig / zsigg(ijq, l) * (qg(ijq, l) - q(ijq, l)))
     666            else
     667              ! u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
     668              ! goto 9999
     669              zz = 0.5 * (zsig - zsigg(ijq, l)) / zsigd(ijq, l)
     670              IF(.NOT.(zz>0..and.zz<=0.5)) THEN
     671                WRITE(lunout, *)'probleme22 au point ij=', ij &
     672                        , '  l=', l
     673                WRITE(lunout, *)'zz=', zz
     674                stop
     675              endif
     676              u_mq(ij, l) = u_mq(ij, l) - masse(ijq, l) * (&
     677                      0.5 * (q(ijq, l) + qg(ijq, l)) * zsigg(ijq, l) &
     678                              + (zsig - zsigg(ijq, l)) * &
     679                              (q(ijq, l) + zz * (qd(ijq, l) - q(ijq, l))))
    676680            endif
    677             u_mq(ij,l)=u_mq(ij,l)-masse(ijq,l)*( &
    678                   0.5*(q(ijq,l)+qg(ijq,l))*zsigg(ijq,l) &
    679                   +(zsig-zsigg(ijq,l))* &
    680                   (q(ijq,l)+zz*(qd(ijq,l)-q(ijq,l))) )
    681          endif
    682   !   fin de la modif
    683               endif
    684            enddo
    685         endif
    686      enddo
     681            !   fin de la modif
     682          endif
     683        enddo
     684      endif
     685    enddo
    687686  ENDIF  ! n0.gt.0
    688687
    689688  !   bouclage en latitude
    690   do l=1,llm
    691     do ij=iip1+iip1,ip1jm,iip1
    692        u_mq(ij,l)=u_mq(ij-iim,l)
     689  do l = 1, llm
     690    do ij = iip1 + iip1, ip1jm, iip1
     691      u_mq(ij, l) = u_mq(ij - iim, l)
    693692    enddo
    694693  enddo
    … …  
    698697  !=================================================================
    699698
    700   do l=1,llm
    701      do ij=iip2+1,ip1jm
    702         new_m=masse(ij,l)+u_m(ij-1,l)-u_m(ij,l)
    703         q(ij,l)=(q(ij,l)*masse(ij,l)+ &
    704               u_mq(ij-1,l)-u_mq(ij,l)) &
    705               /new_m
    706         masse(ij,l)=new_m
    707      enddo
    708   !   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
    709      do ij=iip1+iip1,ip1jm,iip1
    710         q(ij-iim,l)=q(ij,l)
    711         masse(ij-iim,l)=masse(ij,l)
    712      enddo
     699  do l = 1, llm
     700    do ij = iip2 + 1, ip1jm
     701      new_m = masse(ij, l) + u_m(ij - 1, l) - u_m(ij, l)
     702      q(ij, l) = (q(ij, l) * masse(ij, l) + &
     703              u_mq(ij - 1, l) - u_mq(ij, l)) &
     704              / new_m
     705      masse(ij, l) = new_m
     706    enddo
     707    !   Modif Fred 22 03 96 correction d'un bug (les scopy ci-dessous)
     708    do ij = iip1 + iip1, ip1jm, iip1
     709      q(ij - iim, l) = q(ij, l)
     710      masse(ij - iim, l) = masse(ij, l)
     711    enddo
    713712  enddo
    714713
    715714  RETURN
    716715END SUBROUTINE advnx
    717 SUBROUTINE advny(q,qs,qn,masse,v_m)
     716SUBROUTINE advny(q, qs, qn, masse, v_m)
    718717  !
    719718  ! Auteur : F. Hourdin
    … …  
    726725  !
    727726  !   --------------------------------------------------------------------
     727  USE lmdz_iniprint, ONLY: lunout, prt_level
    728728  IMPLICIT NONE
    729729  !
    … …  
    731731  INCLUDE "paramet.h"
    732732  INCLUDE "comgeom.h"
    733   INCLUDE "iniprint.h"
    734733  !
    735734  !
    736735  !   Arguments:
    737736  !   ----------
    738   REAL :: masse(ip1jmp1,llm)
    739   REAL :: v_m( ip1jm,llm )
    740   REAL :: q(ip1jmp1,llm),qn(ip1jmp1,llm),qs(ip1jmp1,llm)
     737  REAL :: masse(ip1jmp1, llm)
     738  REAL :: v_m(ip1jm, llm)
     739  REAL :: q(ip1jmp1, llm), qn(ip1jmp1, llm), qs(ip1jmp1, llm)
    741740  !
    742741  !  Local
    743742  !   ---------
    744743  !
    745   INTEGER :: ij,l
    746   !
    747   REAL :: new_m,zdq,zz
    748   REAL :: zsigs(ip1jmp1),zsign(ip1jmp1),zsig
    749   REAL :: v_mq(ip1jm,llm)
    750   REAL :: convpn,convps,convmpn,convmps,massen,masses
    751   REAL :: zm,zq,zsigm,zsigp,zqm,zqp
     744  INTEGER :: ij, l
     745  !
     746  REAL :: new_m, zdq, zz
     747  REAL :: zsigs(ip1jmp1), zsign(ip1jmp1), zsig
     748  REAL :: v_mq(ip1jm, llm)
     749  REAL :: convpn, convps, convmpn, convmps, massen, masses
     750  REAL :: zm, zq, zsigm, zsigp, zqm, zqp
    752751  REAL :: ssum
    753752  REAL :: prec
    … …  
    755754
    756755  data prec/1.e-15/
    757   do l=1,llm
    758         do ij=1,ip1jmp1
    759            zdq=qn(ij,l)-qs(ij,l)
    760            ! if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) THEN
    761            !    PRINT*,'probleme au point ij=',ij,'  l=',l,'  advnqx'
    762            !    PRINT*,qn(ij,l),q(ij,l),qs(ij,l)
    763            !    qn(ij,l)=q(ij,l)
    764            !    qs(ij,l)=q(ij,l)
    765            ! END IF
    766            IF(abs(zdq)>prec) THEN
    767               zsign(ij)=(q(ij,l)-qs(ij,l))/zdq
    768               zsigs(ij)=1.-zsign(ij)
    769               ! IF(.NOT.(zsign(ij).ge.0..and.zsign(ij).le.1. .AND.
    770   !    s               zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) THEN
    771               !    PRINT*,'probleme au point ij=',ij,'  l=',l
    772               !    PRINT*,'sigs=',zsigs(ij),'  sign=',zsign(ij)
    773               !    stop
    774               ! END IF
    775            else
    776               zsign(ij)=0.5
    777               zsigs(ij)=0.5
    778            endif
    779         enddo
    780 
    781   !   calcul de la pente maximum dans la maille en valeur absolue
    782 
    783    do ij=1,ip1jm
    784       IF (v_m(ij,l)>=0.) THEN
    785          zsigp=zsign(ij+iip1)
    786          zsigm=zsigs(ij+iip1)
    787          zqp=qn(ij+iip1,l)
    788          zqm=qs(ij+iip1,l)
    789          zm=masse(ij+iip1,l)
    790          zq=q(ij+iip1,l)
     756  do l = 1, llm
     757    do ij = 1, ip1jmp1
     758      zdq = qn(ij, l) - qs(ij, l)
     759      ! if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) THEN
     760      !    PRINT*,'probleme au point ij=',ij,'  l=',l,'  advnqx'
     761      !    PRINT*,qn(ij,l),q(ij,l),qs(ij,l)
     762      !    qn(ij,l)=q(ij,l)
     763      !    qs(ij,l)=q(ij,l)
     764      ! END IF
     765      IF(abs(zdq)>prec) THEN
     766        zsign(ij) = (q(ij, l) - qs(ij, l)) / zdq
     767        zsigs(ij) = 1. - zsign(ij)
     768        ! IF(.NOT.(zsign(ij).ge.0..and.zsign(ij).le.1. .AND.
     769        !    s               zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) THEN
     770        !    PRINT*,'probleme au point ij=',ij,'  l=',l
     771        !    PRINT*,'sigs=',zsigs(ij),'  sign=',zsign(ij)
     772        !    stop
     773        ! END IF
    791774      else
    792          zsigm=zsign(ij)
    793          zsigp=zsigs(ij)
    794          zqm=qn(ij,l)
    795          zqp=qs(ij,l)
    796          zm=masse(ij,l)
    797          zq=q(ij,l)
    798       endif
    799       zsig=abs(v_m(ij,l))/zm
    800       IF(zsig==0.) zsigp=0.1
     775        zsign(ij) = 0.5
     776        zsigs(ij) = 0.5
     777      endif
     778    enddo
     779
     780    !   calcul de la pente maximum dans la maille en valeur absolue
     781
     782    do ij = 1, ip1jm
     783      IF (v_m(ij, l)>=0.) THEN
     784        zsigp = zsign(ij + iip1)
     785        zsigm = zsigs(ij + iip1)
     786        zqp = qn(ij + iip1, l)
     787        zqm = qs(ij + iip1, l)
     788        zm = masse(ij + iip1, l)
     789        zq = q(ij + iip1, l)
     790      else
     791        zsigm = zsign(ij)
     792        zsigp = zsigs(ij)
     793        zqm = qn(ij, l)
     794        zqp = qs(ij, l)
     795        zm = masse(ij, l)
     796        zq = q(ij, l)
     797      endif
     798      zsig = abs(v_m(ij, l)) / zm
     799      IF(zsig==0.) zsigp = 0.1
    801800      IF (zsig<=zsigp) THEN
    802           v_mq(ij,l)=v_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
     801        v_mq(ij, l) = v_m(ij, l) * (zqp - 0.5 * zsig / zsigp * (zqp - zq))
    803802      else
    804           zz=0.5*(zsig-zsigp)/zsigm
    805           v_mq(ij,l)=sign(zm,v_m(ij,l))*( 0.5*(zq+zqp)*zsigp &
    806                 +(zsig-zsigp)*(zq+zz*(zqm-zq)) )
    807       endif
    808    enddo
    809   enddo
    810 
    811   do l=1,llm
    812      do ij=iip2,ip1jm
    813         new_m=masse(ij,l) &
    814               +v_m(ij,l)-v_m(ij-iip1,l)
    815         q(ij,l)=(q(ij,l)*masse(ij,l)+v_mq(ij,l)-v_mq(ij-iip1,l)) &
    816               /new_m
    817         masse(ij,l)=new_m
    818      enddo
    819   !.-. ancienne version
    820      convpn=SSUM(iim,v_mq(1,l),1)
    821      convmpn=ssum(iim,v_m(1,l),1)
    822      massen=ssum(iim,masse(1,l),1)
    823      new_m=massen+convmpn
    824      q(1,l)=(q(1,l)*massen+convpn)/new_m
    825      do ij = 1,iip1
    826         q(ij,l)=q(1,l)
    827         masse(ij,l)=new_m*aire(ij)/apoln
    828      enddo
    829 
    830      convps=-SSUM(iim,v_mq(ip1jm-iim,l),1)
    831      convmps=-ssum(iim,v_m(ip1jm-iim,l),1)
    832      masses=ssum(iim,masse(ip1jm+1,l),1)
    833      new_m=masses+convmps
    834      q(ip1jm+1,l)=(q(ip1jm+1,l)*masses+convps)/new_m
    835      do ij = ip1jm+1,ip1jmp1
    836         q(ij,l)=q(ip1jm+1,l)
    837         masse(ij,l)=new_m*aire(ij)/apols
    838      enddo
     803        zz = 0.5 * (zsig - zsigp) / zsigm
     804        v_mq(ij, l) = sign(zm, v_m(ij, l)) * (0.5 * (zq + zqp) * zsigp &
     805                + (zsig - zsigp) * (zq + zz * (zqm - zq)))
     806      endif
     807    enddo
     808  enddo
     809
     810  do l = 1, llm
     811    do ij = iip2, ip1jm
     812      new_m = masse(ij, l) &
     813              + v_m(ij, l) - v_m(ij - iip1, l)
     814      q(ij, l) = (q(ij, l) * masse(ij, l) + v_mq(ij, l) - v_mq(ij - iip1, l)) &
     815              / new_m
     816      masse(ij, l) = new_m
     817    enddo
     818    !.-. ancienne version
     819    convpn = SSUM(iim, v_mq(1, l), 1)
     820    convmpn = ssum(iim, v_m(1, l), 1)
     821    massen = ssum(iim, masse(1, l), 1)
     822    new_m = massen + convmpn
     823    q(1, l) = (q(1, l) * massen + convpn) / new_m
     824    do ij = 1, iip1
     825      q(ij, l) = q(1, l)
     826      masse(ij, l) = new_m * aire(ij) / apoln
     827    enddo
     828
     829    convps = -SSUM(iim, v_mq(ip1jm - iim, l), 1)
     830    convmps = -ssum(iim, v_m(ip1jm - iim, l), 1)
     831    masses = ssum(iim, masse(ip1jm + 1, l), 1)
     832    new_m = masses + convmps
     833    q(ip1jm + 1, l) = (q(ip1jm + 1, l) * masses + convps) / new_m
     834    do ij = ip1jm + 1, ip1jmp1
     835      q(ij, l) = q(ip1jm + 1, l)
     836      masse(ij, l) = new_m * aire(ij) / apols
     837    enddo
    839838  enddo
    840839
    841840  RETURN
    842841END SUBROUTINE advny
    843 SUBROUTINE advnz(q,qh,qb,masse,w_m)
     842SUBROUTINE advnz(q, qh, qb, masse, w_m)
    844843  !
    845844  ! Auteurs:   F.Hourdin
    … …  
    853852  !
    854853  !   --------------------------------------------------------------------
     854  USE lmdz_iniprint, ONLY: lunout, prt_level
    855855  IMPLICIT NONE
    856856  !
    … …  
    858858  INCLUDE "paramet.h"
    859859  INCLUDE "comgeom.h"
    860   INCLUDE "iniprint.h"
    861860  !
    862861  !
    863862  !   Arguments:
    864863  !   ----------
    865   REAL :: masse(ip1jmp1,llm)
    866   REAL :: w_m( ip1jmp1,llm+1)
    867   REAL :: q(ip1jmp1,llm),qb(ip1jmp1,llm),qh(ip1jmp1,llm)
     864  REAL :: masse(ip1jmp1, llm)
     865  REAL :: w_m(ip1jmp1, llm + 1)
     866  REAL :: q(ip1jmp1, llm), qb(ip1jmp1, llm), qh(ip1jmp1, llm)
    868867
    869868  !
    … …  
    871870  !   ---------
    872871  !
    873   INTEGER :: ij,l
    874   !
    875   REAL :: new_m,zdq,zz
    876   REAL :: zsigh(ip1jmp1,llm),zsigb(ip1jmp1,llm),zsig
    877   REAL :: w_mq(ip1jmp1,llm+1)
    878   REAL :: zm,zq,zsigm,zsigp,zqm,zqp
     872  INTEGER :: ij, l
     873  !
     874  REAL :: new_m, zdq, zz
     875  REAL :: zsigh(ip1jmp1, llm), zsigb(ip1jmp1, llm), zsig
     876  REAL :: w_mq(ip1jmp1, llm + 1)
     877  REAL :: zm, zq, zsigm, zsigp, zqm, zqp
    879878  REAL :: prec
    880879  save prec
    … …  
    882881  data prec/1.e-13/
    883882
    884   do l=1,llm
    885         do ij=1,ip1jmp1
    886            zdq=qb(ij,l)-qh(ij,l)
    887            ! if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) THEN
    888            !    PRINT*,'probleme au point ij=',ij,'  l=',l
    889            !    PRINT*,qh(ij,l),q(ij,l),qb(ij,l)
    890            !    qh(ij,l)=q(ij,l)
    891            !    qb(ij,l)=q(ij,l)
    892            ! END IF
    893 
    894            IF(abs(zdq)>prec) THEN
    895               zsigb(ij,l)=(q(ij,l)-qh(ij,l))/zdq
    896               zsigh(ij,l)=1.-zsigb(ij,l)
    897               zsigb(ij,l)=min(max(zsigb(ij,l),0.),1.)
    898            else
    899               zsigb(ij,l)=0.5
    900               zsigh(ij,l)=0.5
    901            endif
    902         enddo
    903    enddo
    904 
    905    ! PRINT*,'ok1'
     883  do l = 1, llm
     884    do ij = 1, ip1jmp1
     885      zdq = qb(ij, l) - qh(ij, l)
     886      ! if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) THEN
     887      !    PRINT*,'probleme au point ij=',ij,'  l=',l
     888      !    PRINT*,qh(ij,l),q(ij,l),qb(ij,l)
     889      !    qh(ij,l)=q(ij,l)
     890      !    qb(ij,l)=q(ij,l)
     891      ! END IF
     892
     893      IF(abs(zdq)>prec) THEN
     894        zsigb(ij, l) = (q(ij, l) - qh(ij, l)) / zdq
     895        zsigh(ij, l) = 1. - zsigb(ij, l)
     896        zsigb(ij, l) = min(max(zsigb(ij, l), 0.), 1.)
     897      else
     898        zsigb(ij, l) = 0.5
     899        zsigh(ij, l) = 0.5
     900      endif
     901    enddo
     902  enddo
     903
     904  ! PRINT*,'ok1'
    906905  !   calcul de la pente maximum dans la maille en valeur absolue
    907    do l=2,llm
    908    do ij=1,ip1jmp1
    909       IF (w_m(ij,l)>=0.) THEN
    910          zsigp=zsigb(ij,l)
    911          zsigm=zsigh(ij,l)
    912          zqp=qb(ij,l)
    913          zqm=qh(ij,l)
    914          zm=masse(ij,l)
    915          zq=q(ij,l)
     906  do l = 2, llm
     907    do ij = 1, ip1jmp1
     908      IF (w_m(ij, l)>=0.) THEN
     909        zsigp = zsigb(ij, l)
     910        zsigm = zsigh(ij, l)
     911        zqp = qb(ij, l)
     912        zqm = qh(ij, l)
     913        zm = masse(ij, l)
     914        zq = q(ij, l)
    916915      else
    917          zsigm=zsigb(ij,l-1)
    918          zsigp=zsigh(ij,l-1)
    919          zqm=qb(ij,l-1)
    920          zqp=qh(ij,l-1)
    921          zm=masse(ij,l-1)
    922          zq=q(ij,l-1)
    923       endif
    924       zsig=abs(w_m(ij,l))/zm
    925       IF(zsig==0.) zsigp=0.1
     916        zsigm = zsigb(ij, l - 1)
     917        zsigp = zsigh(ij, l - 1)
     918        zqm = qb(ij, l - 1)
     919        zqp = qh(ij, l - 1)
     920        zm = masse(ij, l - 1)
     921        zq = q(ij, l - 1)
     922      endif
     923      zsig = abs(w_m(ij, l)) / zm
     924      IF(zsig==0.) zsigp = 0.1
    926925      IF (zsig<=zsigp) THEN
    927           w_mq(ij,l)=w_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
     926        w_mq(ij, l) = w_m(ij, l) * (zqp - 0.5 * zsig / zsigp * (zqp - zq))
    928927      else
    929           zz=0.5*(zsig-zsigp)/zsigm
    930           w_mq(ij,l)=sign(zm,w_m(ij,l))*( 0.5*(zq+zqp)*zsigp &
    931                 +(zsig-zsigp)*(zq+zz*(zqm-zq)) )
    932       endif
    933   enddo
    934   enddo
    935 
    936    do ij=1,ip1jmp1
    937       w_mq(ij,llm+1)=0.
    938       w_mq(ij,1)=0.
    939    enddo
    940 
    941   do l=1,llm
    942      do ij=1,ip1jmp1
    943         new_m=masse(ij,l)+w_m(ij,l+1)-w_m(ij,l)
    944         q(ij,l)=(q(ij,l)*masse(ij,l)+w_mq(ij,l+1)-w_mq(ij,l)) &
    945               /new_m
    946         masse(ij,l)=new_m
    947      enddo
     928        zz = 0.5 * (zsig - zsigp) / zsigm
     929        w_mq(ij, l) = sign(zm, w_m(ij, l)) * (0.5 * (zq + zqp) * zsigp &
     930                + (zsig - zsigp) * (zq + zz * (zqm - zq)))
     931      endif
     932    enddo
     933  enddo
     934
     935  do ij = 1, ip1jmp1
     936    w_mq(ij, llm + 1) = 0.
     937    w_mq(ij, 1) = 0.
     938  enddo
     939
     940  do l = 1, llm
     941    do ij = 1, ip1jmp1
     942      new_m = masse(ij, l) + w_m(ij, l + 1) - w_m(ij, l)
     943      q(ij, l) = (q(ij, l) * masse(ij, l) + w_mq(ij, l + 1) - w_mq(ij, l)) &
     944              / new_m
     945      masse(ij, l) = new_m
     946    enddo
    948947  enddo
    949948  ! PRINT*,'ok3'
Note: See TracChangeset for help on using the changeset viewer.