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

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

Location:
LMDZ6/branches/Amaury_dev/libf/dyn3d_common
Files:
13 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'

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

    r5117 r5118  
    5353
    5454  USE control_mod, ONLY: planet_type
     55  USE lmdz_iniprint, ONLY: lunout, prt_level
    5556
    5657  IMPLICIT NONE
    … …  
    5960  INCLUDE "paramet.h"
    6061  INCLUDE "comgeom.h"
    61   INCLUDE "iniprint.h"
    6262
    6363  ! Ehouarn: for now set these parameters to what is in Earth physics...

  • LMDZ6/branches/Amaury_dev/libf/dyn3d_common/disvert.F90

    r5117 r5118  
    99                         pseudoalt, pa, preff, scaleheight, presinter
    1010  USE logic_mod, ONLY: ok_strato
     11  USE lmdz_iniprint, ONLY: lunout, prt_level
    1112
    1213  IMPLICIT NONE
    … …  
    1415  include "dimensions.h"
    1516  include "paramet.h"
    16   include "iniprint.h"
    1717
    1818!-------------------------------------------------------------------------------

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

    r5117 r5118  
    1111  USE comconst_mod, ONLY: kappa
    1212  USE logic_mod, ONLY: hybrid
     13  USE lmdz_iniprint, ONLY: lunout, prt_level
    1314
    1415  IMPLICIT NONE
    … …  
    1617  include "dimensions.h"
    1718  include "paramet.h"
    18   include "iniprint.h"
    1919  !
    2020  !=======================================================================

  • LMDZ6/branches/Amaury_dev/libf/dyn3d_common/infotrac.F90

    r5117 r5118  
    118118#endif
    119119    USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_INCA, CPPKEY_STRATAER
     120    USE lmdz_iniprint, ONLY: lunout, prt_level
    120121    IMPLICIT NONE
    121122    !==============================================================================================================================
    … …  
    139140    ! Declarations:
    140141    INCLUDE "dimensions.h"
    141     INCLUDE "iniprint.h"
    142142
    143143    !------------------------------------------------------------------------------------------------------------------------------
    … …  
    228228    IF(fType == 1 .AND. ANY(['inca', 'inco']==type_trac)) THEN         !=== FOUND OLD STYLE INCA "traceur.def"
    229229      !---------------------------------------------------------------------------------------------------------------------------
    230         nqo = SIZE(tracers) - nqCO2
    231         CALL Init_chem_inca_trac(nqINCA)                               !--- Get nqINCA from INCA
    232         nbtr = nqINCA + nqCO2                                          !--- Number of tracers passed to phytrac
    233         nqtrue = nbtr + nqo                                            !--- Total number of "true" tracers
    234         IF(ALL([2, 3] /= nqo)) CALL abort_gcm(modname, 'Only 2 or 3 water phases allowed ; found nqo=' // TRIM(int2str(nqo)), 1)
    235         ALLOCATE(hadv(nqtrue), hadv_inca(nqINCA), conv_flg_inca(nqINCA), solsym_inca(nqINCA))
    236         ALLOCATE(vadv(nqtrue), vadv_inca(nqINCA), pbl_flg_inca(nqINCA))
    237         CALL init_transport(solsym_inca, conv_flg_inca, pbl_flg_inca, hadv_inca, vadv_inca)
    238         ALLOCATE(ttr(nqtrue))
    239         ttr(1:nqo + nqCO2) = tracers
    240         ttr(1:nqo)%component = 'lmdz'
    241         ttr(1 + nqo:nqCO2 + nqo)%component = 'co2i'
    242         ttr(1 + nqo + nqCO2:nqtrue)%component = 'inca'
    243         ttr(1 + nqo:nqtrue)%name = [('CO2     ', iq = 1, nqCO2), solsym_inca]
    244         ttr(1 + nqo + nqCO2:nqtrue)%parent = tran0
    245         ttr(1 + nqo + nqCO2:nqtrue)%phase = 'g'
    246         lerr = getKey('hadv', had, ky = tracers(:)%keys)
    247         lerr = getKey('vadv', vad, ky = tracers(:)%keys)
    248         hadv(1:nqo + nqCO2) = had(:); hadv(1 + nqo + nqCO2:nqtrue) = hadv_inca
    249         vadv(1:nqo + nqCO2) = vad(:); vadv(1 + nqo + nqCO2:nqtrue) = vadv_inca
    250         CALL MOVE_ALLOC(FROM = ttr, TO = tracers)
    251         DO iq = 1, nqtrue
    252           t1 => tracers(iq)
    253           CALL addKey('name', t1%name, t1%keys)
    254           CALL addKey('component', t1%component, t1%keys)
    255           CALL addKey('parent', t1%parent, t1%keys)
    256           CALL addKey('phase', t1%phase, t1%keys)
    257         END DO
    258         IF(setGeneration(tracers)) CALL abort_gcm(modname, 'See above', 1) !- SET FIELDS %iGeneration, %gen0Name
    259         DEALLOCATE(had, hadv_inca, vad, vadv_inca, conv_flg_inca, pbl_flg_inca, solsym_inca)
     230      nqo = SIZE(tracers) - nqCO2
     231      CALL Init_chem_inca_trac(nqINCA)                               !--- Get nqINCA from INCA
     232      nbtr = nqINCA + nqCO2                                          !--- Number of tracers passed to phytrac
     233      nqtrue = nbtr + nqo                                            !--- Total number of "true" tracers
     234      IF(ALL([2, 3] /= nqo)) CALL abort_gcm(modname, 'Only 2 or 3 water phases allowed ; found nqo=' // TRIM(int2str(nqo)), 1)
     235      ALLOCATE(hadv(nqtrue), hadv_inca(nqINCA), conv_flg_inca(nqINCA), solsym_inca(nqINCA))
     236      ALLOCATE(vadv(nqtrue), vadv_inca(nqINCA), pbl_flg_inca(nqINCA))
     237      CALL init_transport(solsym_inca, conv_flg_inca, pbl_flg_inca, hadv_inca, vadv_inca)
     238      ALLOCATE(ttr(nqtrue))
     239      ttr(1:nqo + nqCO2) = tracers
     240      ttr(1:nqo)%component = 'lmdz'
     241      ttr(1 + nqo:nqCO2 + nqo)%component = 'co2i'
     242      ttr(1 + nqo + nqCO2:nqtrue)%component = 'inca'
     243      ttr(1 + nqo:nqtrue)%name = [('CO2     ', iq = 1, nqCO2), solsym_inca]
     244      ttr(1 + nqo + nqCO2:nqtrue)%parent = tran0
     245      ttr(1 + nqo + nqCO2:nqtrue)%phase = 'g'
     246      lerr = getKey('hadv', had, ky = tracers(:)%keys)
     247      lerr = getKey('vadv', vad, ky = tracers(:)%keys)
     248      hadv(1:nqo + nqCO2) = had(:); hadv(1 + nqo + nqCO2:nqtrue) = hadv_inca
     249      vadv(1:nqo + nqCO2) = vad(:); vadv(1 + nqo + nqCO2:nqtrue) = vadv_inca
     250      CALL MOVE_ALLOC(FROM = ttr, TO = tracers)
     251      DO iq = 1, nqtrue
     252        t1 => tracers(iq)
     253        CALL addKey('name', t1%name, t1%keys)
     254        CALL addKey('component', t1%component, t1%keys)
     255        CALL addKey('parent', t1%parent, t1%keys)
     256        CALL addKey('phase', t1%phase, t1%keys)
     257      END DO
     258      IF(setGeneration(tracers)) CALL abort_gcm(modname, 'See above', 1) !- SET FIELDS %iGeneration, %gen0Name
     259      DEALLOCATE(had, hadv_inca, vad, vadv_inca, conv_flg_inca, pbl_flg_inca, solsym_inca)
    260260      !---------------------------------------------------------------------------------------------------------------------------
    261261    ELSE                                                              !=== OTHER CASES (OLD OR NEW FORMAT, NO INCA MODULE)
    … …  
    266266      nbtr = nqtrue - COUNT(delPhase(tracers(:)%gen0Name) == 'H2O' &
    267267              .AND. tracers(:)%component == 'lmdz') !--- Number of tracers passed to phytrac
    268         nqINCA = COUNT(tracers(:)%component == 'inca')
     268      nqINCA = COUNT(tracers(:)%component == 'inca')
    269269      lerr = getKey('hadv', hadv, ky = tracers(:)%keys)
    270270      lerr = getKey('vadv', vadv, ky = tracers(:)%keys)

  • LMDZ6/branches/Amaury_dev/libf/dyn3d_common/iniconst.F90

    r5117 r5118  
    1 
    21! $Id$
    32
    … …  
    76  USE IOIPSL
    87  USE comconst_mod, ONLY: im, imp1, jm, jmp1, lllm, lllmm1, lllmp1, &
    9                           unsim, pi, r, kappa, cpp, dtvr, dtphys
     8          unsim, pi, r, kappa, cpp, dtvr, dtphys
    109  USE comvert_mod, ONLY: disvert_type, pressure_exner
    11  
     10  USE lmdz_iniprint, ONLY: lunout, prt_level
     11
    1212  IMPLICIT NONE
    1313
    … …  
    1919  include "dimensions.h"
    2020  include "paramet.h"
    21   include "iniprint.h"
    2221
    23   CHARACTER(LEN=*),parameter :: modname="iniconst"
    24   CHARACTER(LEN=80) :: abort_message
     22  CHARACTER(LEN = *), parameter :: modname = "iniconst"
     23  CHARACTER(LEN = 80) :: abort_message
    2524
    2625
    … …  
    3029  !   ----------------------
    3130
    32   im      = iim
    33   jm      = jjm
    34   lllm    = llm
    35   imp1    = iim
    36   jmp1    = jjm + 1
    37   lllmm1  = llm - 1
    38   lllmp1  = llm + 1
     31  im = iim
     32  jm = jjm
     33  lllm = llm
     34  imp1 = iim
     35  jmp1 = jjm + 1
     36  lllmm1 = llm - 1
     37  lllmp1 = llm + 1
    3938
    4039  !-----------------------------------------------------------------------
    4140
    42   dtphys  = iphysiq * dtvr
    43   unsim   = 1./iim
    44   pi      = 2.*ASIN( 1. )
     41  dtphys = iphysiq * dtvr
     42  unsim = 1. / iim
     43  pi = 2. * ASIN(1.)
    4544
    4645  !-----------------------------------------------------------------------
    4746
    48   r       = cpp * kappa
     47  r = cpp * kappa
    4948
    50   WRITE(lunout,*) trim(modname),': R  CP  Kappa ',r,cpp,kappa
     49  WRITE(lunout, *) trim(modname), ': R  CP  Kappa ', r, cpp, kappa
    5150
    5251  !-----------------------------------------------------------------------
    … …  
    5453  ! vertical discretization: default behavior depends on planet_type flag
    5554  IF (planet_type=="earth") THEN
    56      disvert_type=1
     55    disvert_type = 1
    5756  else
    58      disvert_type=2
     57    disvert_type = 2
    5958  ENDIF
    6059  ! but user can also specify using one or the other in run.def:
    61   CALL getin('disvert_type',disvert_type)
    62   WRITE(lunout,*) trim(modname),': disvert_type=',disvert_type
     60  CALL getin('disvert_type', disvert_type)
     61  WRITE(lunout, *) trim(modname), ': disvert_type=', disvert_type
    6362
    6463  pressure_exner = disvert_type == 1 ! default value
    … …  
    6665
    6766  IF (disvert_type==1) THEN
    68      ! standard case for Earth (automatic generation of levels)
    69      CALL disvert()
     67    ! standard case for Earth (automatic generation of levels)
     68    CALL disvert()
    7069  ELSE IF (disvert_type==2) THEN
    71      ! standard case for planets (levels generated using z2sig.def file)
    72      CALL disvert_noterre
     70    ! standard case for planets (levels generated using z2sig.def file)
     71    CALL disvert_noterre
    7372  else
    74      WRITE(abort_message,*) "Wrong value for disvert_type: ", disvert_type
    75      CALL abort_gcm(modname,abort_message,0)
     73    WRITE(abort_message, *) "Wrong value for disvert_type: ", disvert_type
     74    CALL abort_gcm(modname, abort_message, 0)
    7675  ENDIF
    7776

  • LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inidissip.F90

    r5117 r5118  
    1717  USE lmdz_libmath, ONLY: minmax
    1818  USE lmdz_ran1, ONLY: ran1
     19  USE lmdz_iniprint, ONLY: lunout, prt_level
    1920
    2021  IMPLICIT NONE
    … …  
    2223  include "paramet.h"
    2324  include "comdissipn.h"
    24   include "iniprint.h"
    2525
    2626  LOGICAL, INTENT(IN) :: lstardis

  • LMDZ6/branches/Amaury_dev/libf/dyn3d_common/initdynav.F90

    r5117 r5118  
    11! $Id$
    22
    3 SUBROUTINE initdynav(day0,anne0,tstep,t_ops,t_wrt)
     3SUBROUTINE initdynav(day0, anne0, tstep, t_ops, t_wrt)
    44
    55  USE IOIPSL
    66  USE infotrac, ONLY: nqtot
    7   USE com_io_dyn_mod, ONLY: histaveid,histvaveid,histuaveid, &
    8        dynhistave_file,dynhistvave_file,dynhistuave_file
     7  USE com_io_dyn_mod, ONLY: histaveid, histvaveid, histuaveid, &
     8          dynhistave_file, dynhistvave_file, dynhistuave_file
    99  USE comconst_mod, ONLY: pi
    1010  USE comvert_mod, ONLY: presnivs
    1111  USE temps_mod, ONLY: itau_dyn
    1212  USE lmdz_description, ONLY: descript
    13  
     13  USE lmdz_iniprint, ONLY: lunout, prt_level
     14
    1415  IMPLICIT NONE
    1516
    … …  
    3839  include "paramet.h"
    3940  include "comgeom.h"
    40   include "iniprint.h"
    4141
    4242  !   Arguments
    … …  
    5151  REAL zjulian
    5252  INTEGER iq
    53   REAL rlong(iip1,jjp1), rlat(iip1,jjp1)
     53  REAL rlong(iip1, jjp1), rlat(iip1, jjp1)
    5454  INTEGER uhoriid, vhoriid, thoriid, zvertiid
    55   INTEGER ii,jj
     55  INTEGER ii, jj
    5656  INTEGER zan, dayref
    5757
    … …  
    6464  !  Appel a histbeg: creation du fichier netcdf et initialisations diverses
    6565
    66 
    6766  zan = anne0
    6867  dayref = day0
    … …  
    7170
    7271  do jj = 1, jjp1
    73      do ii = 1, iip1
    74         rlong(ii,jj) = rlonv(ii) * 180. / pi
    75         rlat(ii,jj) = rlatu(jj) * 180. / pi
    76      enddo
     72    do ii = 1, iip1
     73      rlong(ii, jj) = rlonv(ii) * 180. / pi
     74      rlat(ii, jj) = rlatu(jj) * 180. / pi
     75    enddo
    7776  enddo
    7877
    … …  
    8079  ! Restriction de IOIPSL: seulement 2 coordonnees dans le meme fichier
    8180  ! Grille Scalaire
    82   CALL histbeg(dynhistave_file, iip1, rlong(:,1), jjp1, rlat(1,:), &
    83        1, iip1, 1, jjp1, &
    84        tau0, zjulian, tstep, thoriid,histaveid)
     81  CALL histbeg(dynhistave_file, iip1, rlong(:, 1), jjp1, rlat(1, :), &
     82          1, iip1, 1, jjp1, &
     83          tau0, zjulian, tstep, thoriid, histaveid)
    8584
    8685  ! Creation du fichier histoire pour les grilles en V et U (oblige
    … …  
    8988  ! Grille V
    9089  do jj = 1, jjm
    91      do ii = 1, iip1
    92         rlong(ii,jj) = rlonv(ii) * 180. / pi
    93         rlat(ii,jj) = rlatv(jj) * 180. / pi
    94      enddo
     90    do ii = 1, iip1
     91      rlong(ii, jj) = rlonv(ii) * 180. / pi
     92      rlat(ii, jj) = rlatv(jj) * 180. / pi
     93    enddo
    9594  enddo
    9695
    97   CALL histbeg(dynhistvave_file, iip1, rlong(:,1), jjm, rlat(1,:), &
    98        1, iip1, 1, jjm, &
    99        tau0, zjulian, tstep, vhoriid,histvaveid)
     96  CALL histbeg(dynhistvave_file, iip1, rlong(:, 1), jjm, rlat(1, :), &
     97          1, iip1, 1, jjm, &
     98          tau0, zjulian, tstep, vhoriid, histvaveid)
    10099  ! Grille U
    101100  do jj = 1, jjp1
    102      do ii = 1, iip1
    103         rlong(ii,jj) = rlonu(ii) * 180. / pi
    104         rlat(ii,jj) = rlatu(jj) * 180. / pi
    105      enddo
     101    do ii = 1, iip1
     102      rlong(ii, jj) = rlonu(ii) * 180. / pi
     103      rlat(ii, jj) = rlatu(jj) * 180. / pi
     104    enddo
    106105  enddo
    107106
    108   CALL histbeg(dynhistuave_file, iip1, rlong(:,1),jjp1, rlat(1,:), &
    109        1, iip1, 1, jjp1, &
    110        tau0, zjulian, tstep, uhoriid,histuaveid)
     107  CALL histbeg(dynhistuave_file, iip1, rlong(:, 1), jjp1, rlat(1, :), &
     108          1, iip1, 1, jjp1, &
     109          tau0, zjulian, tstep, uhoriid, histuaveid)
    111110
    112111  !  Appel a histvert pour la grille verticale
    113112
    114   CALL histvert(histaveid,'presnivs','Niveaux Pression approximatifs','mb', &
    115        llm, presnivs/100., zvertiid,'down')
    116   CALL histvert(histuaveid,'presnivs','Niveaux Pression approximatifs','mb', &
    117        llm, presnivs/100., zvertiid,'down')
    118   CALL histvert(histvaveid,'presnivs','Niveaux Pression approximatifs','mb', &
    119        llm, presnivs/100., zvertiid,'down')
     113  CALL histvert(histaveid, 'presnivs', 'Niveaux Pression approximatifs', 'mb', &
     114          llm, presnivs / 100., zvertiid, 'down')
     115  CALL histvert(histuaveid, 'presnivs', 'Niveaux Pression approximatifs', 'mb', &
     116          llm, presnivs / 100., zvertiid, 'down')
     117  CALL histvert(histvaveid, 'presnivs', 'Niveaux Pression approximatifs', 'mb', &
     118          llm, presnivs / 100., zvertiid, 'down')
    120119
    121120  !  Appels a histdef pour la definition des variables a sauvegarder
    … …  
    124123
    125124  CALL histdef(histuaveid, 'u', 'vent u moyen ', &
    126        'm/s', iip1, jjp1, uhoriid, llm, 1, llm, zvertiid, &
    127        32, 'ave(X)', t_ops, t_wrt)
     125          'm/s', iip1, jjp1, uhoriid, llm, 1, llm, zvertiid, &
     126          32, 'ave(X)', t_ops, t_wrt)
    128127
    129128  !  Vents V
    130129
    131130  CALL histdef(histvaveid, 'v', 'vent v moyen', &
    132        'm/s', iip1, jjm, vhoriid, llm, 1, llm, zvertiid, &
    133        32, 'ave(X)', t_ops, t_wrt)
     131          'm/s', iip1, jjm, vhoriid, llm, 1, llm, zvertiid, &
     132          32, 'ave(X)', t_ops, t_wrt)
    134133
    135134
    … …  
    137136
    138137  CALL histdef(histaveid, 'temp', 'temperature moyenne', 'K', &
    139        iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
    140        32, 'ave(X)', t_ops, t_wrt)
     138          iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
     139          32, 'ave(X)', t_ops, t_wrt)
    141140
    142141  !  Temperature potentielle
    143142
    144143  CALL histdef(histaveid, 'theta', 'temperature potentielle', 'K', &
    145        iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
    146        32, 'ave(X)', t_ops, t_wrt)
     144          iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
     145          32, 'ave(X)', t_ops, t_wrt)
    147146
    148147  !  Geopotentiel
    149148
    150149  CALL histdef(histaveid, 'phi', 'geopotentiel moyen', '-', &
    151        iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
    152        32, 'ave(X)', t_ops, t_wrt)
     150          iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
     151          32, 'ave(X)', t_ops, t_wrt)
    153152
    154153  !  Traceurs
    … …  
    164163
    165164  CALL histdef(histaveid, 'masse', 'masse', 'kg', &
    166        iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
    167        32, 'ave(X)', t_ops, t_wrt)
     165          iip1, jjp1, thoriid, llm, 1, llm, zvertiid, &
     166          32, 'ave(X)', t_ops, t_wrt)
    168167
    169168  !  Pression au sol
    170169
    171170  CALL histdef(histaveid, 'ps', 'pression naturelle au sol', 'Pa', &
    172        iip1, jjp1, thoriid, 1, 1, 1, -99, &
    173        32, 'ave(X)', t_ops, t_wrt)
     171          iip1, jjp1, thoriid, 1, 1, 1, -99, &
     172          32, 'ave(X)', t_ops, t_wrt)
    174173
    175174  !  Geopotentiel au sol

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

    r5117 r5118  
    1010  USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn
    1111  USE lmdz_description, ONLY: descript
     12  USE lmdz_iniprint, ONLY: lunout, prt_level
    1213
    1314  IMPLICIT NONE
    … …  
    4344  include "paramet.h"
    4445  include "comgeom.h"
    45   include "iniprint.h"
    4646
    4747  !   Arguments

  • LMDZ6/branches/Amaury_dev/libf/dyn3d_common/inithist.F90

    r5116 r5118  
    1111  USE temps_mod, ONLY: itau_dyn
    1212  USE lmdz_description, ONLY: descript
     13  USE lmdz_iniprint, ONLY: lunout, prt_level
    1314
    1415  IMPLICIT NONE
    … …  
    4243  include "paramet.h"
    4344  include "comgeom.h"
    44   include "iniprint.h"
    4545
    4646  !   Arguments

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

    r5117 r5118  
    1 
    21! $Id$
    32
    43SUBROUTINE sortvarc &
    5         (itau,ucov,teta,ps,masse,pk,phis,vorpot,phi,bern,dp,time , &
    6         vcov )
     4        (itau, ucov, teta, ps, masse, pk, phis, vorpot, phi, bern, dp, time, &
     5        vcov)
    76
    87  USE control_mod, ONLY: resetvarc
    98  USE comconst_mod, ONLY: dtvr, daysec, g, rad, omeg
    109  USE logic_mod, ONLY: read_start
    11   USE ener_mod, ONLY: etot,ptot,ztot,stot,ang, &
    12         etot0,ptot0,ztot0,stot0,ang0, &
    13         rmsdpdt,rmsv
     10  USE ener_mod, ONLY: etot, ptot, ztot, stot, ang, &
     11          etot0, ptot0, ztot0, stot0, ang0, &
     12          rmsdpdt, rmsv
    1413  USE lmdz_filtreg, ONLY: filtreg
     14  USE lmdz_iniprint, ONLY: lunout, prt_level
    1515  IMPLICIT NONE
    1616
    … …  
    3434  INCLUDE "paramet.h"
    3535  INCLUDE "comgeom.h"
    36   INCLUDE "iniprint.h"
    3736
    3837  !   Arguments:
    3938  !   ----------
    4039
    41   INTEGER,INTENT(IN) :: itau
    42   REAL,INTENT(IN) :: ucov(ip1jmp1,llm)
    43   REAL,INTENT(IN) :: teta(ip1jmp1,llm)
    44   REAL,INTENT(IN) :: masse(ip1jmp1,llm)
    45   REAL,INTENT(IN) :: vcov(ip1jm,llm)
    46   REAL,INTENT(IN) :: ps(ip1jmp1)
    47   REAL,INTENT(IN) :: phis(ip1jmp1)
    48   REAL,INTENT(IN) :: vorpot(ip1jm,llm)
    49   REAL,INTENT(IN) :: phi(ip1jmp1,llm)
    50   REAL,INTENT(IN) :: bern(ip1jmp1,llm)
    51   REAL,INTENT(IN) :: dp(ip1jmp1)
    52   REAL,INTENT(IN) :: time
    53   REAL,INTENT(IN) :: pk(ip1jmp1,llm)
     40  INTEGER, INTENT(IN) :: itau
     41  REAL, INTENT(IN) :: ucov(ip1jmp1, llm)
     42  REAL, INTENT(IN) :: teta(ip1jmp1, llm)
     43  REAL, INTENT(IN) :: masse(ip1jmp1, llm)
     44  REAL, INTENT(IN) :: vcov(ip1jm, llm)
     45  REAL, INTENT(IN) :: ps(ip1jmp1)
     46  REAL, INTENT(IN) :: phis(ip1jmp1)
     47  REAL, INTENT(IN) :: vorpot(ip1jm, llm)
     48  REAL, INTENT(IN) :: phi(ip1jmp1, llm)
     49  REAL, INTENT(IN) :: bern(ip1jmp1, llm)
     50  REAL, INTENT(IN) :: dp(ip1jmp1)
     51  REAL, INTENT(IN) :: time
     52  REAL, INTENT(IN) :: pk(ip1jmp1, llm)
    5453
    5554  !   Local:
    5655  !   ------
    5756
    58   REAL :: vor(ip1jm),bernf(ip1jmp1,llm),ztotl(llm)
    59   REAL :: etotl(llm),stotl(llm),rmsvl(llm),angl(llm),ge(ip1jmp1)
    60   REAL :: cosphi(ip1jm),omegcosp(ip1jm)
    61   REAL :: dtvrs1j,rjour,heure,radsg,radomeg
    62   REAL :: massebxy(ip1jm,llm)
     57  REAL :: vor(ip1jm), bernf(ip1jmp1, llm), ztotl(llm)
     58  REAL :: etotl(llm), stotl(llm), rmsvl(llm), angl(llm), ge(ip1jmp1)
     59  REAL :: cosphi(ip1jm), omegcosp(ip1jm)
     60  REAL :: dtvrs1j, rjour, heure, radsg, radomeg
     61  REAL :: massebxy(ip1jm, llm)
    6362  INTEGER :: l, ij, imjmp1
    6463
    6564  REAL :: SSUM
    66   LOGICAL,SAVE :: firstcal=.TRUE.
    67   CHARACTER(LEN=*),PARAMETER :: modname="sortvarc"
     65  LOGICAL, SAVE :: firstcal = .TRUE.
     66  CHARACTER(LEN = *), PARAMETER :: modname = "sortvarc"
    6867
    6968  !-----------------------------------------------------------------------
    7069  ! Ehouarn: when no initialization fields from file, resetvarc should be
    71        ! set to false
    72    IF (firstcal) THEN
    73      IF (.NOT.read_start) THEN
    74        resetvarc=.TRUE.
    75      endif
    76    endif
    77 
    78    dtvrs1j   = dtvr/daysec
    79    rjour     = REAL( INT( itau * dtvrs1j ))
    80    heure     = ( itau*dtvrs1j-rjour ) * 24.
    81    imjmp1    = iim * jjp1
    82    IF(ABS(heure - 24.)<=0.0001 ) heure = 0.
    83   !
    84    CALL massbarxy ( masse, massebxy )
     70  ! set to false
     71  IF (firstcal) THEN
     72    IF (.NOT.read_start) THEN
     73      resetvarc = .TRUE.
     74    endif
     75  endif
     76
     77  dtvrs1j = dtvr / daysec
     78  rjour = REAL(INT(itau * dtvrs1j))
     79  heure = (itau * dtvrs1j - rjour) * 24.
     80  imjmp1 = iim * jjp1
     81  IF(ABS(heure - 24.)<=0.0001) heure = 0.
     82  !
     83  CALL massbarxy (masse, massebxy)
    8584
    8685  !   .....  Calcul  de  rmsdpdt  .....
    8786
    88    ge(:)=dp(:)*dp(:)
    89 
    90    rmsdpdt = SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
    91   !
    92    rmsdpdt = daysec* 1.e-2 * SQRT(rmsdpdt/imjmp1)
    93 
    94    CALL SCOPY( ijp1llm,bern,1,bernf,1 )
    95    CALL filtreg(bernf,jjp1,llm,-2,2,.TRUE.,1)
     87  ge(:) = dp(:) * dp(:)
     88
     89  rmsdpdt = SSUM(ip1jmp1, ge, 1) - SSUM(jjp1, ge, iip1)
     90  !
     91  rmsdpdt = daysec * 1.e-2 * SQRT(rmsdpdt / imjmp1)
     92
     93  CALL SCOPY(ijp1llm, bern, 1, bernf, 1)
     94  CALL filtreg(bernf, jjp1, llm, -2, 2, .TRUE., 1)
    9695
    9796  !   .....  Calcul du moment  angulaire   .....
    9897
    99    radsg    = rad /g
    100    radomeg = rad * omeg
    101   !
    102    DO ij=iip2,ip1jm
    103       cosphi( ij ) = COS(rlatu((ij-1)/iip1+1))
    104       omegcosp(ij) = radomeg  * cosphi(ij)
    105    ENDDO
     98  radsg = rad / g
     99  radomeg = rad * omeg
     100  !
     101  DO ij = iip2, ip1jm
     102    cosphi(ij) = COS(rlatu((ij - 1) / iip1 + 1))
     103    omegcosp(ij) = radomeg * cosphi(ij)
     104  ENDDO
    106105
    107106  !  ...  Calcul  de l'energie,de l'enstrophie,de l'entropie et de rmsv  .
    108107
    109    DO l=1,llm
    110       DO ij = 1,ip1jm
    111          vor(ij)=vorpot(ij,l)*vorpot(ij,l)*massebxy(ij,l)
    112       ENDDO
    113       ztotl(l)=(SSUM(ip1jm,vor,1)-SSUM(jjm,vor,iip1))
    114 
    115       DO ij = 1,ip1jmp1
    116          ge(ij)= masse(ij,l)*(phis(ij)+teta(ij,l)*pk(ij,l) + &
    117                bernf(ij,l)-phi(ij,l))
    118       ENDDO
    119       etotl(l) = SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
    120 
    121       DO   ij  = 1, ip1jmp1
    122          ge(ij) = masse(ij,l)*teta(ij,l)
    123       ENDDO
    124       stotl(l)= SSUM(ip1jmp1,ge,1) - SSUM(jjp1,ge,iip1)
    125 
    126       DO ij=1,ip1jmp1
    127          ge(ij)=masse(ij,l)*AMAX1(bernf(ij,l)-phi(ij,l),0.)
    128       ENDDO
    129       rmsvl(l)=2.*(SSUM(ip1jmp1,ge,1)-SSUM(jjp1,ge,iip1))
    130 
    131       DO ij =iip2,ip1jm
    132          ge(ij)=(ucov(ij,l)/cu(ij)+omegcosp(ij))*masse(ij,l) * &
    133                cosphi(ij)
    134       ENDDO
    135       angl(l) = rad * &
    136             (SSUM(ip1jm-iip1,ge(iip2),1)-SSUM(jjm-1,ge(iip2),iip1))
     108  DO l = 1, llm
     109    DO ij = 1, ip1jm
     110      vor(ij) = vorpot(ij, l) * vorpot(ij, l) * massebxy(ij, l)
     111    ENDDO
     112    ztotl(l) = (SSUM(ip1jm, vor, 1) - SSUM(jjm, vor, iip1))
     113
     114    DO ij = 1, ip1jmp1
     115      ge(ij) = masse(ij, l) * (phis(ij) + teta(ij, l) * pk(ij, l) + &
     116              bernf(ij, l) - phi(ij, l))
     117    ENDDO
     118    etotl(l) = SSUM(ip1jmp1, ge, 1) - SSUM(jjp1, ge, iip1)
     119
     120    DO   ij = 1, ip1jmp1
     121      ge(ij) = masse(ij, l) * teta(ij, l)
     122    ENDDO
     123    stotl(l) = SSUM(ip1jmp1, ge, 1) - SSUM(jjp1, ge, iip1)
     124
     125    DO ij = 1, ip1jmp1
     126      ge(ij) = masse(ij, l) * AMAX1(bernf(ij, l) - phi(ij, l), 0.)
     127    ENDDO
     128    rmsvl(l) = 2. * (SSUM(ip1jmp1, ge, 1) - SSUM(jjp1, ge, iip1))
     129
     130    DO ij = iip2, ip1jm
     131      ge(ij) = (ucov(ij, l) / cu(ij) + omegcosp(ij)) * masse(ij, l) * &
     132              cosphi(ij)
     133    ENDDO
     134    angl(l) = rad * &
     135            (SSUM(ip1jm - iip1, ge(iip2), 1) - SSUM(jjm - 1, ge(iip2), iip1))
    137136  ENDDO
    138137
    139       DO ij=1,ip1jmp1
    140         ge(ij)= ps(ij)*aire(ij)
    141       ENDDO
    142   ptot  = SSUM(ip1jmp1,ge,1)-SSUM(jjp1,ge,iip1)
    143   etot  = SSUM(     llm, etotl, 1 )
    144   ztot  = SSUM(     llm, ztotl, 1 )
    145   stot  = SSUM(     llm, stotl, 1 )
    146   rmsv  = SSUM(     llm, rmsvl, 1 )
    147   ang   = SSUM(     llm,  angl, 1 )
     138  DO ij = 1, ip1jmp1
     139    ge(ij) = ps(ij) * aire(ij)
     140  ENDDO
     141  ptot = SSUM(ip1jmp1, ge, 1) - SSUM(jjp1, ge, iip1)
     142  etot = SSUM(llm, etotl, 1)
     143  ztot = SSUM(llm, ztotl, 1)
     144  stot = SSUM(llm, stotl, 1)
     145  rmsv = SSUM(llm, rmsvl, 1)
     146  ang = SSUM(llm, angl, 1)
    148147
    149148  IF (firstcal.AND.resetvarc) THEN
    150      WRITE(lunout,3500) itau, rjour, heure, time
    151      WRITE(lunout,*) trim(modname), &
    152            ' WARNING!!! Recomputing initial values of : '
    153      WRITE(lunout,*) 'ptot,rmsdpdt,etot,ztot,stot,rmsv,ang'
    154      WRITE(lunout,*) ptot,rmsdpdt,etot,ztot,stot,rmsv,ang
    155      etot0 = etot
    156      ptot0 = ptot
    157      ztot0 = ztot
    158      stot0 = stot
    159      ang0 = ang
     149    WRITE(lunout, 3500) itau, rjour, heure, time
     150    WRITE(lunout, *) trim(modname), &
     151            ' WARNING!!! Recomputing initial values of : '
     152    WRITE(lunout, *) 'ptot,rmsdpdt,etot,ztot,stot,rmsv,ang'
     153    WRITE(lunout, *) ptot, rmsdpdt, etot, ztot, stot, rmsv, ang
     154    etot0 = etot
     155    ptot0 = ptot
     156    ztot0 = ztot
     157    stot0 = stot
     158    ang0 = ang
    160159  END IF
    161160
    … …  
    163162  ! are zero, which can happen when using iniacademic)
    164163  IF (etot0/=0) THEN
    165     etot= etot/etot0
    166   else
    167     etot=1.
    168   ENDIF
    169   rmsv= SQRT(rmsv/ptot)
     164    etot = etot / etot0
     165  else
     166    etot = 1.
     167  ENDIF
     168  rmsv = SQRT(rmsv / ptot)
    170169  IF (ptot0/=0) THEN
    171     ptot= ptot/ptot0
    172   else
    173     ptot=1.
     170    ptot = ptot / ptot0
     171  else
     172    ptot = 1.
    174173  ENDIF
    175174  IF (ztot0/=0) THEN
    176     ztot= ztot/ztot0
    177   else
    178     ztot=1.
     175    ztot = ztot / ztot0
     176  else
     177    ztot = 1.
    179178  ENDIF
    180179  IF (stot0/=0) THEN
    181     stot= stot/stot0
    182   else
    183     stot=1.
     180    stot = stot / stot0
     181  else
     182    stot = 1.
    184183  ENDIF
    185184  IF (ang0/=0) THEN
    186     ang = ang /ang0
    187   else
    188     ang=1.
     185    ang = ang / ang0
     186  else
     187    ang = 1.
    189188  ENDIF
    190189
    191190  firstcal = .FALSE.
    192191
    193   WRITE(lunout,3500) itau, rjour, heure, time
    194   WRITE(lunout,4000) ptot,rmsdpdt,etot,ztot,stot,rmsv,ang
    195 
    196 3500   FORMAT(10("*"),4x,'pas',i7,5x,'jour',f9.0,'heure',f5.1,4x &
    197              ,'date',f14.4,4x,10("*"))
    198 4000   FORMAT(10x,'masse',4x,'rmsdpdt',7x,'energie',2x,'enstrophie' &
    199              ,2x,'entropie',3x,'rmsv',4x,'mt.ang',/,'GLOB  ' &
    200              ,f10.6,e13.6,5f10.3/ &
    201              )
     192  WRITE(lunout, 3500) itau, rjour, heure, time
     193  WRITE(lunout, 4000) ptot, rmsdpdt, etot, ztot, stot, rmsv, ang
     194
     195  3500   FORMAT(10("*"), 4x, 'pas', i7, 5x, 'jour', f9.0, 'heure', f5.1, 4x &
     196          , 'date', f14.4, 4x, 10("*"))
     197  4000   FORMAT(10x, 'masse', 4x, 'rmsdpdt', 7x, 'energie', 2x, 'enstrophie' &
     198          , 2x, 'entropie', 3x, 'rmsv', 4x, 'mt.ang', /, 'GLOB  ' &
     199          , f10.6, e13.6, 5f10.3/ &
     200          )
    202201END SUBROUTINE sortvarc
    203202

  • LMDZ6/branches/Amaury_dev/libf/dyn3d_common/writedynav.F90

    r5117 r5118  
    99  USE temps_mod, ONLY: itau_dyn
    1010  USE lmdz_description, ONLY: descript
     11  USE lmdz_iniprint, ONLY: lunout, prt_level
    1112
    1213  IMPLICIT NONE
    … …  
    3334  include "paramet.h"
    3435  include "comgeom.h"
    35   include "iniprint.h"
    3636
    3737  !   Arguments
    3838
    39   REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) 
    40   REAL teta(ip1jmp1*llm), phi(ip1jmp1, llm), ppk(ip1jmp1*llm)     
    41   REAL ps(ip1jmp1), masse(ip1jmp1, llm)                   
    42   REAL phis(ip1jmp1)                 
     39  REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm)
     40  REAL teta(ip1jmp1 * llm), phi(ip1jmp1, llm), ppk(ip1jmp1 * llm)
     41  REAL ps(ip1jmp1), masse(ip1jmp1, llm)
     42  REAL phis(ip1jmp1)
    4343  REAL q(ip1jmp1, llm, nqtot)
    4444  INTEGER time
    … …  
    4747  !   Variables locales
    4848
    49   INTEGER ndex2d(ip1jmp1), ndexu(ip1jmp1*llm), ndexv(ip1jm*llm)
     49  INTEGER ndex2d(ip1jmp1), ndexu(ip1jmp1 * llm), ndexv(ip1jm * llm)
    5050  INTEGER iq, ii, ll
    51   REAL tm(ip1jmp1*llm)
     51  REAL tm(ip1jmp1 * llm)
    5252  REAL vnat(ip1jm, llm), unat(ip1jmp1, llm)
    5353  LOGICAL ok_sync
    … …  
    7474  !  Vents U
    7575
    76   CALL histwrite(histuaveid, 'u', itau_w, unat,  &
    77        iip1*jjp1*llm, ndexu)
     76  CALL histwrite(histuaveid, 'u', itau_w, unat, &
     77          iip1 * jjp1 * llm, ndexu)
    7878
    7979  !  Vents V
    8080
    81   CALL histwrite(histvaveid, 'v', itau_w, vnat,  &
    82        iip1*jjm*llm, ndexv)
     81  CALL histwrite(histvaveid, 'v', itau_w, vnat, &
     82          iip1 * jjm * llm, ndexv)
    8383
    8484  !  Temperature potentielle moyennee
    8585
    86   CALL histwrite(histaveid, 'theta', itau_w, teta,  &
    87        iip1*jjp1*llm, ndexu)
     86  CALL histwrite(histaveid, 'theta', itau_w, teta, &
     87          iip1 * jjp1 * llm, ndexu)
    8888
    8989  !  Temperature moyennee
    9090
    9191  do ii = 1, ijp1llm
    92      tm(ii) = teta(ii) * ppk(ii)/cpp
     92    tm(ii) = teta(ii) * ppk(ii) / cpp
    9393  enddo
    94   CALL histwrite(histaveid, 'temp', itau_w, tm,  &
    95        iip1*jjp1*llm, ndexu)
     94  CALL histwrite(histaveid, 'temp', itau_w, tm, &
     95          iip1 * jjp1 * llm, ndexu)
    9696
    9797  !  Geopotentiel
    9898
    99   CALL histwrite(histaveid, 'phi', itau_w, phi,  &
    100        iip1*jjp1*llm, ndexu)
     99  CALL histwrite(histaveid, 'phi', itau_w, phi, &
     100          iip1 * jjp1 * llm, ndexu)
    101101
    102102  !  Traceurs
    … …  
    109109  !  Masse
    110110
    111   CALL histwrite(histaveid, 'masse', itau_w, masse,  &
    112        iip1*jjp1*llm, ndexu)
     111  CALL histwrite(histaveid, 'masse', itau_w, masse, &
     112          iip1 * jjp1 * llm, ndexu)
    113113
    114114  !  Pression au sol
    115115
    116   CALL histwrite(histaveid, 'ps', itau_w, ps, iip1*jjp1, ndex2d)
     116  CALL histwrite(histaveid, 'ps', itau_w, ps, iip1 * jjp1, ndex2d)
    117117
    118118  ! Geopotentiel au sol
    … …  
    121121
    122122  IF (ok_sync) THEN
    123      CALL histsync(histaveid)
    124      CALL histsync(histvaveid)
    125      CALL histsync(histuaveid)
     123    CALL histsync(histaveid)
     124    CALL histsync(histvaveid)
     125    CALL histsync(histuaveid)
    126126  ENDIF
    127127

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

    r5117 r5118  
    88  USE temps_mod, ONLY: itau_dyn
    99  USE lmdz_description, ONLY: descript
     10  USE lmdz_iniprint, ONLY: lunout, prt_level
    1011
    1112  IMPLICIT NONE
    … …  
    3637  include "paramet.h"
    3738  include "comgeom.h"
    38   include "iniprint.h"
    3939
    4040  !
Note: See TracChangeset for help on using the changeset viewer.