Ignore:
Timestamp:
Jul 23, 2024, 7:14:34 PM (8 weeks ago)
Author:
abarral
Message:

Replace 1DUTILS.h by module lmdz_1dutils.f90
Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files)
Convert *.F to *.f90
Fix gradsdef.h formatting
Remove unnecessary "RETURN" at the end of functions/subroutines

File:
1 moved

Legend:

Unmodified
Added
Removed

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

    r5104 r5105  
    22! $Header$
    33
    4       SUBROUTINE advn(q,masse,w,pbaru,pbarv,pdt,mode)
    5 c
    6 c     Auteur : F. Hourdin
    7 c
    8 c    ********************************************************************
    9 c     Shema  d'advection " pseudo amont " .
    10 c    ********************************************************************
    11 c     q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
    12 c
    13 c   pbaru,pbarv,w flux de masse en u ,v ,w
    14 c   pdt pas de temps
    15 c
    16 c   --------------------------------------------------------------------
    17       IMPLICIT NONE
    18 c
    19       include "dimensions.h"
    20       include "paramet.h"
    21       include "comgeom.h"
    22       include "iniprint.h"
    23 
    24 c
    25 c   Arguments:
    26 c   ----------
    27       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
    32 c
    33 c      Local
    34 c   ---------
    35 c
    36       INTEGER i,ij,l,j,ii
    37       integer ijlqmin,iqmin,jqmin,lqmin
    38       integer ismin
    39 c
    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./
    56 
    57       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
    70       ENDDO
    71 
    72       DO ij=1,ip1jmp1
    73          mw(ij,llm+1)=0.
    74       ENDDO
    75 
    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
    99       enddo
    100 
    101 c     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)
    108 c     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)
    113 c     CALL minmaxq(zq,qmin,qmax,'apres vlx     ')
    114 
    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
    122       enddo
    123 
    124       RETURN
    125       END
    126 
    127       SUBROUTINE advnqx(q,qg,qd)
    128 c
    129 c     Auteurs:   Calcul des valeurs de q aux point u.
    130 c
    131 c   --------------------------------------------------------------------
    132       IMPLICIT NONE
    133 c
    134       INCLUDE "dimensions.h"
    135       INCLUDE "paramet.h"
    136       INCLUDE "iniprint.h"
    137 c
    138 c
    139 c   Arguments:
    140 c   ----------
    141       real q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)
    142 c
    143 c      Local
    144 c   ---------
    145 c
    146       INTEGER ij,l
    147 c
    148       real dxqu(ip1jmp1),zqu(ip1jmp1)
    149       real zqmax(ip1jmp1),zqmin(ip1jmp1)
    150       logical extremum(ip1jmp1)
    151 
    152       integer mode
    153       save mode
    154       data mode/1/
    155 
    156 c   calcul des pentes en u:
    157 c   -----------------------
    158       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
     4SUBROUTINE advn(q,masse,w,pbaru,pbarv,pdt,mode)
     5  !
     6  ! Auteur : F. Hourdin
     7  !
     8  !    ********************************************************************
     9  ! Shema  d'advection " pseudo amont " .
     10  !    ********************************************************************
     11  ! q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
     12  !
     13  !   pbaru,pbarv,w flux de masse en u ,v ,w
     14  !   pdt pas de temps
     15  !
     16  !   --------------------------------------------------------------------
     17  IMPLICIT NONE
     18  !
     19  include "dimensions.h"
     20  include "paramet.h"
     21  include "comgeom.h"
     22  include "iniprint.h"
     23
     24  !
     25  !   Arguments:
     26  !   ----------
     27  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
     32  !
     33  !  Local
     34  !   ---------
     35  !
     36  INTEGER :: i,ij,l,j,ii
     37  integer :: ijlqmin,iqmin,jqmin,lqmin
     38  integer :: ismin
     39  !
     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./
     56
     57  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
     70  ENDDO
     71
     72  DO ij=1,ip1jmp1
     73     mw(ij,llm+1)=0.
     74  ENDDO
     75
     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
     99  enddo
     100
     101  ! 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)
     108  ! 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)
     113  ! CALL minmaxq(zq,qmin,qmax,'apres vlx     ')
     114
     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
     122  enddo
     123
     124  RETURN
     125END SUBROUTINE advn
     126
     127SUBROUTINE advnqx(q,qg,qd)
     128  !
     129  ! Auteurs:   Calcul des valeurs de q aux point u.
     130  !
     131  !   --------------------------------------------------------------------
     132  IMPLICIT NONE
     133  !
     134  INCLUDE "dimensions.h"
     135  INCLUDE "paramet.h"
     136  INCLUDE "iniprint.h"
     137  !
     138  !
     139  !   Arguments:
     140  !   ----------
     141  real :: q(ip1jmp1,llm),qg(ip1jmp1,llm),qd(ip1jmp1,llm)
     142  !
     143  !  Local
     144  !   ---------
     145  !
     146  INTEGER :: ij,l
     147  !
     148  real :: dxqu(ip1jmp1),zqu(ip1jmp1)
     149  real :: zqmax(ip1jmp1),zqmin(ip1jmp1)
     150  logical :: extremum(ip1jmp1)
     151
     152  integer :: mode
     153  save mode
     154  data mode/1/
     155
     156  !   calcul des pentes en u:
     157  !   -----------------------
     158  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
     165  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
     208        if(extremum(ij)) then
     209           qg(ij,l)=q(ij,l)
     210           qd(ij,l)=q(ij,l)
     211        else
     212           qd(ij,l)=zqu(ij)
     213           qg(ij,l)=zqu(ij-1)
     214        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
     2398888   continue
     240  enddo
     241  endif
     242  RETURN
     243END SUBROUTINE advnqx
     244SUBROUTINE advnqy(q,qs,qn)
     245  !
     246  ! Auteurs:   Calcul des valeurs de q aux point v.
     247  !
     248  !   --------------------------------------------------------------------
     249  IMPLICIT NONE
     250  !
     251  INCLUDE "dimensions.h"
     252  INCLUDE "paramet.h"
     253  INCLUDE "iniprint.h"
     254  !
     255  !
     256  !   Arguments:
     257  !   ----------
     258  real :: q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)
     259  !
     260  !  Local
     261  !   ---------
     262  !
     263  INTEGER :: ij,l
     264  !
     265  real :: dyqv(ip1jm),zqv(ip1jm,llm)
     266  real :: zqmax(ip1jm),zqmin(ip1jm)
     267  logical :: extremum(ip1jmp1)
     268
     269  integer :: mode
     270  save mode
     271  data mode/1/
     272
     273  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
     280  else
     281
     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
     317        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)
     322        else
     323           qs(ij,l)=zqv(ij,l)
     324           qn(ij,l)=zqv(ij-iip1,l)
     325        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
     336  endif
     337  RETURN
     338END SUBROUTINE advnqy
     339
     340SUBROUTINE advnqz(q,qh,qb)
     341  !
     342  ! Auteurs:   Calcul des valeurs de q aux point v.
     343  !
     344  !   --------------------------------------------------------------------
     345  IMPLICIT NONE
     346  !
     347  INCLUDE "dimensions.h"
     348  INCLUDE "paramet.h"
     349  INCLUDE "iniprint.h"
     350  !
     351  !
     352  !   Arguments:
     353  !   ----------
     354  real :: q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)
     355  !
     356  !  Local
     357  !   ---------
     358  !
     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)
     364
     365  integer :: mode
     366  save mode
     367
     368  data mode/1/
     369
     370  !   calcul des pentes en u:
     371  !   -----------------------
     372
     373  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
     380  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)
     429        else
     430           qh(ij,l)=zqw(ij,l+1)
     431           qb(ij,l)=zqw(ij,l)
     432        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
     450
     451  endif
     452
     453  RETURN
     454END SUBROUTINE advnqz
     455
     456SUBROUTINE advnx(q,qg,qd,masse,u_m,mode)
     457  !
     458  ! Auteur : F. Hourdin
     459  !
     460  !    ********************************************************************
     461  ! Shema  d'advection " pseudo amont " .
     462  !    ********************************************************************
     463  ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
     464  !
     465  !
     466  !   --------------------------------------------------------------------
     467  IMPLICIT NONE
     468  !
     469  include "dimensions.h"
     470  include "paramet.h"
     471  include "iniprint.h"
     472  !
     473  !
     474  !   Arguments:
     475  !   ----------
     476  integer :: mode
     477  real :: masse(ip1jmp1,llm)
     478  real :: u_m( ip1jmp1,llm )
     479  real :: q(ip1jmp1,llm),qd(ip1jmp1,llm),qg(ip1jmp1,llm)
     480  !
     481  !  Local
     482  !   ---------
     483  !
     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)
     494
     495  real :: prec
     496  save prec
     497
     498  data prec/1.e-15/
     499
     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           ! endif
     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              ! endif
     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
     527
     528  !   calcul de la pente maximum dans la maille en valeur absolue
     529
     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)
    165539      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 c   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
    208             if(extremum(ij)) then
    209                qg(ij,l)=q(ij,l)
    210                qd(ij,l)=q(ij,l)
    211             else
    212                qd(ij,l)=zqu(ij)
    213                qg(ij,l)=zqu(ij-1)
     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
     551      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
     558      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
     566      endif
     567      ! if(zsig.lt.0.) then
     568      !    PRINT*,'au point ij=',ij,'  l=',l,'  sig=',zsig
     569      !    stop
     570      ! endif
     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
     579  enddo
     580
     581  !=================================================================
     582  !   SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES
     583  !=================================================================
     584  !   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)
     595  enddo
     596
     597  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
     610              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
    214646            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      s         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      s         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
     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
     676            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
     687  endif  ! n0.gt.0
     688
     689  !   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)
     693    enddo
     694  enddo
     695
     696  !=================================================================
     697  !   CALCUL DE LA CONVERGENCE DES FLUX
     698  !=================================================================
     699
     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
     713  enddo
     714
     715  RETURN
     716END SUBROUTINE advnx
     717SUBROUTINE advny(q,qs,qn,masse,v_m)
     718  !
     719  ! Auteur : F. Hourdin
     720  !
     721  !    ********************************************************************
     722  ! Shema  d'advection " pseudo amont " .
     723  !    ********************************************************************
     724  ! nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
     725  !
     726  !
     727  !   --------------------------------------------------------------------
     728  IMPLICIT NONE
     729  !
     730  INCLUDE "dimensions.h"
     731  INCLUDE "paramet.h"
     732  INCLUDE "comgeom.h"
     733  INCLUDE "iniprint.h"
     734  !
     735  !
     736  !   Arguments:
     737  !   ----------
     738  real :: masse(ip1jmp1,llm)
     739  real :: v_m( ip1jm,llm )
     740  real :: q(ip1jmp1,llm),qn(ip1jmp1,llm),qs(ip1jmp1,llm)
     741  !
     742  !  Local
     743  !   ---------
     744  !
     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
     752  real :: ssum
     753  real :: prec
     754  save prec
     755
     756  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           ! endif
     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              ! endif
     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)
     791      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)
    241798      endif
    242       RETURN
    243       END
    244       SUBROUTINE advnqy(q,qs,qn)
    245 c
    246 c     Auteurs:   Calcul des valeurs de q aux point v.
    247 c
    248 c   --------------------------------------------------------------------
    249       IMPLICIT NONE
    250 c
    251       INCLUDE "dimensions.h"
    252       INCLUDE "paramet.h"
    253       INCLUDE "iniprint.h"
    254 c
    255 c
    256 c   Arguments:
    257 c   ----------
    258       real q(ip1jmp1,llm),qs(ip1jmp1,llm),qn(ip1jmp1,llm)
    259 c
    260 c      Local
    261 c   ---------
    262 c
    263       INTEGER ij,l
    264 c
    265       real dyqv(ip1jm),zqv(ip1jm,llm)
    266       real zqmax(ip1jm),zqmin(ip1jm)
    267       logical extremum(ip1jmp1)
    268 
    269       integer mode
    270       save mode
    271       data mode/1/
    272 
    273       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
     799      zsig=abs(v_m(ij,l))/zm
     800      if(zsig==0.) zsigp=0.1
     801      if (zsig<=zsigp) then
     802          v_mq(ij,l)=v_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
    280803      else
    281 
    282 c   calcul des pentes en u:
    283 c   -----------------------
    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 c 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 c   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
    317             if(extremum(ij)) then
    318                qs(ij,l)=q(ij,l)
    319                qn(ij,l)=q(ij,l)
    320 c              if (.not.extremum(ij-iip1)) qs(ij-iip1,l)=q(ij,l)
    321 c              if (.not.extremum(ij+iip1)) qn(ij+iip1,l)=q(ij,l)
    322             else
    323                qs(ij,l)=zqv(ij,l)
    324                qn(ij,l)=zqv(ij-iip1,l)
    325             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
     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)) )
    336807      endif
    337       RETURN
    338       END
    339 
    340       SUBROUTINE advnqz(q,qh,qb)
    341 c
    342 c     Auteurs:   Calcul des valeurs de q aux point v.
    343 c
    344 c   --------------------------------------------------------------------
    345       IMPLICIT NONE
    346 c
    347       INCLUDE "dimensions.h"
    348       INCLUDE "paramet.h"
    349       INCLUDE "iniprint.h"
    350 c
    351 c
    352 c   Arguments:
    353 c   ----------
    354       real q(ip1jmp1,llm),qh(ip1jmp1,llm),qb(ip1jmp1,llm)
    355 c
    356 c      Local
    357 c   ---------
    358 c
    359       INTEGER ij,l
    360 c
    361       real dzqw(ip1jmp1,llm+1),zqw(ip1jmp1,llm+1)
    362       real zqmax(ip1jmp1,llm),zqmin(ip1jmp1,llm)
    363       logical extremum(ip1jmp1,llm)
    364 
    365       integer mode
    366       save mode
    367 
    368       data mode/1/
    369 
    370 c   calcul des pentes en u:
    371 c   -----------------------
    372 
    373       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
     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
     839  enddo
     840
     841  RETURN
     842END SUBROUTINE advny
     843SUBROUTINE advnz(q,qh,qb,masse,w_m)
     844  !
     845  ! Auteurs:   F.Hourdin
     846  !
     847  !    ********************************************************************
     848  ! Shema  d'advection " pseudo amont " .
     849  ! b designe le bas et h le haut
     850  ! il y a une correspondance entre le b en z et le d en x
     851  !    ********************************************************************
     852  !
     853  !
     854  !   --------------------------------------------------------------------
     855  IMPLICIT NONE
     856  !
     857  INCLUDE "dimensions.h"
     858  INCLUDE "paramet.h"
     859  INCLUDE "comgeom.h"
     860  INCLUDE "iniprint.h"
     861  !
     862  !
     863  !   Arguments:
     864  !   ----------
     865  real :: masse(ip1jmp1,llm)
     866  real :: w_m( ip1jmp1,llm+1)
     867  real :: q(ip1jmp1,llm),qb(ip1jmp1,llm),qh(ip1jmp1,llm)
     868
     869  !
     870  !  Local
     871  !   ---------
     872  !
     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
     879  real :: prec
     880  save prec
     881
     882  data prec/1.e-13/
     883
     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           ! endif
     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'
     906  !   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)
    380916      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 c 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 c   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)
    429             else
    430                qh(ij,l)=zqw(ij,l+1)
    431                qb(ij,l)=zqw(ij,l)
    432             endif
    433          enddo
    434       enddo
    435 c     do l=2,llm-1
    436 c        do ij=1,ip1jmp1
    437 c           if(extremum(ij,l)) then
    438 c              if (.not.extremum(ij,l-1)) qh(ij,l-1)=q(ij,l)
    439 c              if (.not.extremum(ij,l+1)) qb(ij,l+1)=q(ij,l)
    440 c           endif
    441 c        enddo
    442 c     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
    450 
     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)
    451923      endif
    452 
    453       RETURN
    454       END
    455 
    456       SUBROUTINE advnx(q,qg,qd,masse,u_m,mode)
    457 c
    458 c     Auteur : F. Hourdin
    459 c
    460 c    ********************************************************************
    461 c     Shema  d'advection " pseudo amont " .
    462 c    ********************************************************************
    463 c     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
    464 c
    465 c
    466 c   --------------------------------------------------------------------
    467       IMPLICIT NONE
    468 c
    469       include "dimensions.h"
    470       include "paramet.h"
    471       include "iniprint.h"
    472 c
    473 c
    474 c   Arguments:
    475 c   ----------
    476       integer mode
    477       real masse(ip1jmp1,llm)
    478       real u_m( ip1jmp1,llm )
    479       real q(ip1jmp1,llm),qd(ip1jmp1,llm),qg(ip1jmp1,llm)
    480 c
    481 c      Local
    482 c   ---------
    483 c
    484       INTEGER i,j,ij,l,indu(ip1jmp1),niju,iju,ijq
    485       integer n0,nl(llm)
    486 c
    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)
    494 
    495       real prec
    496       save prec
    497 
    498       data prec/1.e-15/
    499 
    500       do l=1,llm
    501             do ij=iip2,ip1jm
    502                zdq=qd(ij,l)-qg(ij,l)
    503 c              if((qd(ij,l)-q(ij,l))*(q(ij,l)-qg(ij,l)).lt.0.) then
    504 c                 PRINT*,'probleme au point ij=',ij,'  l=',l
    505 c                 PRINT*,qd(ij,l),q(ij,l),qg(ij,l)
    506 c                 qd(ij,l)=q(ij,l)
    507 c                 qg(ij,l)=q(ij,l)
    508 c              endif
    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 c                 if(.not.(zsigd(ij,l).ge.0..and.zsigd(ij,l).le.1. .and.
    513 c    s               zsigg(ij,l).ge.0..or.zsigg(ij,l).le.1.) ) then
    514 c                    PRINT*,'probleme au point ij=',ij,'  l=',l
    515 c                    PRINT*,'sigg=',zsigg(ij,l),'  sigd=',zsigd(ij,l)
    516 c                    PRINT*,'q d,c,g ',qd(ij,l),q(ij,l),qg(ij,l),zdq
    517 c                    stop
    518 c                 endif
    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
    527 
    528 c   calcul de la pente maximum dans la maille en valeur absolue
    529 
    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)
    539           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
    551           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      s           sign(zm,u_m(ij,l))*(zsigp*zqp+(zsig-zsigp)*zqm)
    557              endif
    558           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      s          +(zsig-zsigp)*(zq+zz*(zqm-zq)) )
    565              endif
    566           endif
    567 c         if(zsig.lt.0.) then
    568 c            PRINT*,'au point ij=',ij,'  l=',l,'  sig=',zsig
    569 c            stop
    570 c         endif
    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
    579       enddo
    580 
    581 c=================================================================
    582 C   SCHEMA SEMI-LAGRAGIEN EN X DANS LES REGIONS POLAIRES
    583 c=================================================================
    584 c   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)
    595       enddo
    596 
    597       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 c   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
    610                   endif
    611                enddo
    612                niju=iju
    613 c              PRINT*,'niju,nl',niju,nl(l)
    614 
    615 c  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 c   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 c   MODIFS SPECIFIQUES DU SCHEMA
    632 c   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      s          -0.5*zsig/zsigd(ijq,l)*(qd(ijq,l)-q(ijq,l)))
    637              else
    638 c               u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
    639 c         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      s               '  l=',l
    644                      WRITE(lunout,*)'zz=',zz
    645                      stop
    646                 endif
    647                 u_mq(ij,l)=u_mq(ij,l)+masse(ijq,l)*(
    648      s          0.5*(q(ijq,l)+qd(ijq,l))*zsigd(ijq,l)
    649      s        +(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 c   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 c   ajout de la maille non completement advectee
    662 c 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      s          -0.5*zsig/zsigg(ijq,l)*(qg(ijq,l)-q(ijq,l)))
    667              else
    668 c               u_mq(ij,l)=u_mq(ij,l)+zu_m*q(ijq,l)
    669 c           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      s               ,'  l=',l
    674                      WRITE(lunout,*)'zz=',zz
    675                      stop
    676                 endif
    677                 u_mq(ij,l)=u_mq(ij,l)-masse(ijq,l)*(
    678      s          0.5*(q(ijq,l)+qg(ijq,l))*zsigg(ijq,l)
    679      s          +(zsig-zsigg(ijq,l))*
    680      s           (q(ijq,l)+zz*(qd(ijq,l)-q(ijq,l))) )
    681              endif
    682 c   fin de la modif
    683                   endif
    684                enddo
    685             endif
    686          enddo
    687       endif  ! n0.gt.0
    688 
    689 c   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)
    693         enddo
    694       enddo
    695 
    696 c=================================================================
    697 c   CALCUL DE LA CONVERGENCE DES FLUX
    698 c=================================================================
    699 
    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 c   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
    713       enddo
    714 
    715       RETURN
    716       END
    717       SUBROUTINE advny(q,qs,qn,masse,v_m)
    718 c
    719 c     Auteur : F. Hourdin
    720 c
    721 c    ********************************************************************
    722 c     Shema  d'advection " pseudo amont " .
    723 c    ********************************************************************
    724 c     nq,iq,q,pbaru,pbarv,w sont des arguments d'entree  pour le s-pg ....
    725 c
    726 c
    727 c   --------------------------------------------------------------------
    728       IMPLICIT NONE
    729 c
    730       INCLUDE "dimensions.h"
    731       INCLUDE "paramet.h"
    732       INCLUDE "comgeom.h"
    733       INCLUDE "iniprint.h"
    734 c
    735 c
    736 c   Arguments:
    737 c   ----------
    738       real masse(ip1jmp1,llm)
    739       real v_m( ip1jm,llm )
    740       real q(ip1jmp1,llm),qn(ip1jmp1,llm),qs(ip1jmp1,llm)
    741 c
    742 c      Local
    743 c   ---------
    744 c
    745       INTEGER ij,l
    746 c
    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
    752       real ssum
    753       real prec
    754       save prec
    755 
    756       data prec/1.e-15/
    757       do l=1,llm
    758             do ij=1,ip1jmp1
    759                zdq=qn(ij,l)-qs(ij,l)
    760 c              if((qn(ij,l)-q(ij,l))*(q(ij,l)-qs(ij,l)).lt.0.) then
    761 c                 PRINT*,'probleme au point ij=',ij,'  l=',l,'  advnqx'
    762 c                 PRINT*,qn(ij,l),q(ij,l),qs(ij,l)
    763 c                 qn(ij,l)=q(ij,l)
    764 c                 qs(ij,l)=q(ij,l)
    765 c              endif
    766                if(abs(zdq)>prec) then
    767                   zsign(ij)=(q(ij,l)-qs(ij,l))/zdq
    768                   zsigs(ij)=1.-zsign(ij)
    769 c                 if(.not.(zsign(ij).ge.0..and.zsign(ij).le.1. .and.
    770 c    s               zsigs(ij).ge.0..or.zsigs(ij).le.1.) ) then
    771 c                    PRINT*,'probleme au point ij=',ij,'  l=',l
    772 c                    PRINT*,'sigs=',zsigs(ij),'  sign=',zsign(ij)
    773 c                    stop
    774 c                 endif
    775                else
    776                   zsign(ij)=0.5
    777                   zsigs(ij)=0.5
    778                endif
    779             enddo
    780 
    781 c   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)
    791           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
    801           if (zsig<=zsigp) then
    802               v_mq(ij,l)=v_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
    803           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      s        +(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 c.-. 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
    839       enddo
    840 
    841       RETURN
    842       END
    843       SUBROUTINE advnz(q,qh,qb,masse,w_m)
    844 c
    845 c     Auteurs:   F.Hourdin
    846 c
    847 c    ********************************************************************
    848 c     Shema  d'advection " pseudo amont " .
    849 c     b designe le bas et h le haut
    850 c     il y a une correspondance entre le b en z et le d en x
    851 c    ********************************************************************
    852 c
    853 c
    854 c   --------------------------------------------------------------------
    855       IMPLICIT NONE
    856 c
    857       INCLUDE "dimensions.h"
    858       INCLUDE "paramet.h"
    859       INCLUDE "comgeom.h"
    860       INCLUDE "iniprint.h"
    861 c
    862 c
    863 c   Arguments:
    864 c   ----------
    865       real masse(ip1jmp1,llm)
    866       real w_m( ip1jmp1,llm+1)
    867       real q(ip1jmp1,llm),qb(ip1jmp1,llm),qh(ip1jmp1,llm)
    868 
    869 c
    870 c      Local
    871 c   ---------
    872 c
    873       INTEGER ij,l
    874 c
    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
    879       real prec
    880       save prec
    881 
    882       data prec/1.e-13/
    883 
    884       do l=1,llm
    885             do ij=1,ip1jmp1
    886                zdq=qb(ij,l)-qh(ij,l)
    887 c              if((qh(ij,l)-q(ij,l))*(q(ij,l)-qb(ij,l)).lt.0.) then
    888 c                 PRINT*,'probleme au point ij=',ij,'  l=',l
    889 c                 PRINT*,qh(ij,l),q(ij,l),qb(ij,l)
    890 c                 qh(ij,l)=q(ij,l)
    891 c                 qb(ij,l)=q(ij,l)
    892 c              endif
    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 c      PRINT*,'ok1'
    906 c   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)
    916           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
    926           if (zsig<=zsigp) then
    927               w_mq(ij,l)=w_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
    928           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      s        +(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
    948       enddo
    949 c     PRINT*,'ok3'
    950       RETURN
    951       END
     924      zsig=abs(w_m(ij,l))/zm
     925      if(zsig==0.) zsigp=0.1
     926      if (zsig<=zsigp) then
     927          w_mq(ij,l)=w_m(ij,l)*(zqp-0.5*zsig/zsigp*(zqp-zq))
     928      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
     948  enddo
     949  ! PRINT*,'ok3'
     950  RETURN
     951END SUBROUTINE advnz
Note: See TracChangeset for help on using the changeset viewer.