Ignore:
Timestamp:
Jul 19, 2024, 6:40:44 PM (6 months ago)
Author:
Laurent Fairhead
Message:

Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks
AB, LF

Location:
LMDZ6/trunk/libf/dyn3d_common
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/dyn3d_common/advx.F

    r5077 r5084  
    9595
    9696C  -------------------------------------
    97       DO j = 1,jjp1
     97      DO 300 j = 1,jjp1
    9898         NUM(j) = 1
    99       END DO
     99  300 CONTINUE
    100100      sqi = 0.
    101101      sqf = 0.
    … …  
    121121C  ugri est en kg/s
    122122
    123       DO l = 1,llm
    124          DO j = 1,jjm+1
    125             DO i = 1,iip1
     123      DO 500 l = 1,llm
     124         DO 500 j = 1,jjm+1
     125            DO 500 i = 1,iip1 
    126126C            ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
    127127             ugri (i,j,llm+1-l) = pbaru (i,j,l)
    128       END DO
    129       END DO
    130       END DO
     128  500 CONTINUE
    131129
    132130
    … …  
    139137C  boucle principale sur les niveaux et les latitudes
    140138C
    141       DO L=1,NIV
    142       DO K=lati,latf
     139      DO 1 L=1,NIV
     140      DO 1 K=lati,latf
    143141C
    144142C  initialisation
    … …  
    146144C  program assumes periodic boundaries in X
    147145C
    148       DO I=2,LON
     146      DO 10 I=2,LON
    149147         SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
    150       END DO
     148 10   CONTINUE
    151149      SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
    152150C
    … …  
    156154      LONK=LON/NUMK
    157155C
    158       IF(NUMK>1) THEN
    159 C
    160       DO I=1,LON
     156      IF(NUMK.GT.1) THEN
     157C
     158      DO 111 I=1,LON
    161159         TM(I)=0.
    162       END DO
    163       DO JV=1,NTRA
    164       DO I=1,LON
     160 111  CONTINUE
     161      DO 112 JV=1,NTRA
     162      DO 1120 I=1,LON
    165163         T0(I,JV)=0.
    166164         TX(I,JV)=0.
    167165         TY(I,JV)=0.
    168166         TZ(I,JV)=0.
    169       END DO
    170       END DO
    171 C
    172       DO I2=1,NUMK
    173 C
    174          DO I=1,LONK
     167 1120 CONTINUE
     168 112  CONTINUE
     169C
     170      DO 11 I2=1,NUMK
     171C
     172         DO 113 I=1,LONK
    175173            I3=(I-1)*NUMK+I2
    176174            TM(I)=TM(I)+SM(I3,K,L)
    177175            ALF(I)=SM(I3,K,L)/TM(I)
    178176            ALF1(I)=1.-ALF(I)
    179       END DO
     177 113     CONTINUE
    180178C
    181179         DO  JV=1,NTRA
    … …  
    192190         ENDDO
    193191C
    194       END DO
     192 11   CONTINUE
    195193C
    196194      ELSE
    197195C
    198       DO I=1,LON
     196      DO 115 I=1,LON
    199197         TM(I)=SM(I,K,L)
    200       END DO
    201       DO JV=1,NTRA
    202       DO I=1,LON
     198 115  CONTINUE
     199      DO 116 JV=1,NTRA
     200      DO 1160 I=1,LON
    203201         T0(I,JV)=S0(I,K,L,JV)
    204202         TX(I,JV)=sx(I,K,L,JV)
    205203         TY(I,JV)=sy(I,K,L,JV)
    206204         TZ(I,JV)=sz(I,K,L,JV)
    207       END DO
    208       END DO
    209 C
    210       ENDIF
    211 C
    212       DO I=1,LONK
     205 1160 CONTINUE
     206 116  CONTINUE
     207C
     208      ENDIF
     209C
     210      DO 117 I=1,LONK
    213211         UEXT(I)=UGRI(I*NUMK,K,L)
    214       END DO
     212 117  CONTINUE
    215213C
    216214C  place limits on appropriate moments before transport
    … …  
    219217      IF(.NOT.LIMIT) GO TO 13
    220218C
    221       DO JV=1,NTRA
    222       DO I=1,LONK
     219      DO 12 JV=1,NTRA
     220      DO 120 I=1,LONK
    223221        TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV))
    224       END DO
    225       END DO
     222 120  CONTINUE
     223 12   CONTINUE
    226224C
    227225 13   CONTINUE
    … …  
    233231C  flux from IP to I if U(I).lt.0
    234232C
    235       DO I=1,LONK-1
    236          IF(UEXT(I)<0.) THEN
     233      DO 140 I=1,LONK-1
     234         IF(UEXT(I).LT.0.) THEN
    237235           FM(I)=-UEXT(I)*DTX
    238236           ALF(I)=FM(I)/TM(I+1)
    239237           TM(I+1)=TM(I+1)-FM(I)
    240238         ENDIF
    241       END DO
     239 140  CONTINUE
    242240C
    243241      I=LONK
    244       IF(UEXT(I)<0.) THEN
     242      IF(UEXT(I).LT.0.) THEN
    245243        FM(I)=-UEXT(I)*DTX
    246244        ALF(I)=FM(I)/TM(1)
    … …  
    250248C  flux from I to IP if U(I).gt.0
    251249C
    252       DO I=1,LONK
    253          IF(UEXT(I)>=0.) THEN
     250      DO 141 I=1,LONK
     251         IF(UEXT(I).GE.0.) THEN
    254252           FM(I)=UEXT(I)*DTX
    255253           ALF(I)=FM(I)/TM(I)
    256254           TM(I)=TM(I)-FM(I)
    257255         ENDIF
    258       END DO
    259 C
    260       DO I=1,LONK
     256 141  CONTINUE
     257C
     258      DO 142 I=1,LONK
    261259         ALFQ(I)=ALF(I)*ALF(I)
    262260         ALF1(I)=1.-ALF(I)
    263261         ALF1Q(I)=ALF1(I)*ALF1(I)
    264       END DO
    265 C
    266       DO JV=1,NTRA
    267       DO I=1,LONK-1
    268 C
    269          IF(UEXT(I)<0.) THEN
     262 142  CONTINUE
     263C
     264      DO 150 JV=1,NTRA
     265      DO 1500 I=1,LONK-1
     266C
     267         IF(UEXT(I).LT.0.) THEN
    270268C
    271269           F0(I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*TX(I+1,JV) )
    … …  
    281279         ENDIF
    282280C
    283       END DO
    284       END DO
     281 1500 CONTINUE
     282 150  CONTINUE
    285283C
    286284      I=LONK
    287       IF(UEXT(I)<0.) THEN
    288 C
    289         DO JV=1,NTRA
     285      IF(UEXT(I).LT.0.) THEN
     286C
     287        DO 151 JV=1,NTRA
    290288C
    291289           F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*TX(1,JV) )
    … …  
    299297           TZ(1,JV)=TZ(1,JV)-FZ(I,JV)
    300298C
    301       END DO
    302 C
    303       ENDIF
    304 C
    305       DO JV=1,NTRA
    306       DO I=1,LONK
    307 C
    308          IF(UEXT(I)>=0.) THEN
     299 151    CONTINUE
     300C
     301      ENDIF
     302C
     303      DO 152 JV=1,NTRA
     304      DO 1520 I=1,LONK
     305C
     306         IF(UEXT(I).GE.0.) THEN
    309307C
    310308           F0(I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*TX(I,JV) )
    … …  
    320318         ENDIF
    321319C
    322       END DO
    323       END DO
     320 1520 CONTINUE
     321 152  CONTINUE
    324322C
    325323C  puts the temporary moments Fi into appropriate neighboring boxes
    326324C
    327       DO I=1,LONK
    328          IF(UEXT(I)<0.) THEN
     325      DO 160 I=1,LONK
     326         IF(UEXT(I).LT.0.) THEN
    329327           TM(I)=TM(I)+FM(I)
    330328           ALF(I)=FM(I)/TM(I)
    331329         ENDIF
    332       END DO
    333 C
    334       DO I=1,LONK-1
    335          IF(UEXT(I)>=0.) THEN
     330 160  CONTINUE
     331C
     332      DO 161 I=1,LONK-1
     333         IF(UEXT(I).GE.0.) THEN
    336334           TM(I+1)=TM(I+1)+FM(I)
    337335           ALF(I)=FM(I)/TM(I+1)
    338336         ENDIF
    339       END DO
     337 161  CONTINUE
    340338C
    341339      I=LONK
    342       IF(UEXT(I)>=0.) THEN
     340      IF(UEXT(I).GE.0.) THEN
    343341        TM(1)=TM(1)+FM(I)
    344342        ALF(I)=FM(I)/TM(1)
    345343      ENDIF
    346344C
    347       DO I=1,LONK
     345      DO 162 I=1,LONK
    348346         ALF1(I)=1.-ALF(I)
    349       END DO
    350 C
    351       DO JV=1,NTRA
    352       DO I=1,LONK
    353 C
    354          IF(UEXT(I)<0.) THEN
     347 162  CONTINUE
     348C
     349      DO 170 JV=1,NTRA
     350      DO 1700 I=1,LONK
     351C
     352         IF(UEXT(I).LT.0.) THEN
    355353C
    356354           TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV)
    … …  
    362360         ENDIF
    363361C
    364       END DO
    365       END DO
    366 C
    367       DO JV=1,NTRA
    368       DO I=1,LONK-1
    369 C
    370          IF(UEXT(I)>=0.) THEN
     362 1700 CONTINUE
     363 170  CONTINUE
     364C
     365      DO 171 JV=1,NTRA
     366      DO 1710 I=1,LONK-1
     367C
     368         IF(UEXT(I).GE.0.) THEN
    371369C
    372370           TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV)
    … …  
    378376         ENDIF
    379377C
    380       END DO
    381       END DO
     378 1710 CONTINUE
     379 171  CONTINUE
    382380C
    383381      I=LONK
    384       IF(UEXT(I)>=0.) THEN
    385         DO JV=1,NTRA
     382      IF(UEXT(I).GE.0.) THEN
     383        DO 172 JV=1,NTRA
    386384           TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV)
    387385           T0(1,JV)=T0(1,JV)+F0(I,JV)
    … …  
    389387           TY(1,JV)=TY(1,JV)+FY(I,JV)
    390388           TZ(1,JV)=TZ(1,JV)+FZ(I,JV)
    391       END DO
     389 172    CONTINUE
    392390      ENDIF
    393391C
    394392C  retour aux mailles d'origine (passage des Tij aux Sij)
    395393C
    396       IF(NUMK>1) THEN
    397 C
    398       DO I2=1,NUMK
    399 C
    400          DO I=1,LONK
     394      IF(NUMK.GT.1) THEN
     395C
     396      DO 180 I2=1,NUMK
     397C
     398         DO 180 I=1,LONK
    401399C
    402400            I3=I2+(I-1)*NUMK
    … …  
    409407            ALF1Q(I)=ALF1(I)*ALF1(I)
    410408C
    411       END DO
    412       END DO
     409 180     CONTINUE
    413410C
    414411         DO  JV=1,NTRA
    … …  
    434431      ELSE
    435432C
    436       DO I=1,LON
     433      DO 190 I=1,LON
    437434         SM(I,K,L)=TM(I)
    438       END DO
    439       DO JV=1,NTRA
    440       DO I=1,LON
     435 190  CONTINUE
     436      DO 191 JV=1,NTRA
     437      DO 1910 I=1,LON
    441438         S0(I,K,L,JV)=T0(I,JV)
    442439         sx(I,K,L,JV)=TX(I,JV)
    443440         sy(I,K,L,JV)=TY(I,JV)
    444441         sz(I,K,L,JV)=TZ(I,JV)
    445       END DO
    446       END DO
    447 C
    448       ENDIF
    449 C
    450       END DO
    451       END DO
     442 1910 CONTINUE
     443 191  CONTINUE
     444C
     445      ENDIF
     446C
     447 1    CONTINUE
    452448C
    453449C ----------- AA Test en fin de ADVX ------ Controle des S*

  • LMDZ6/trunk/libf/dyn3d_common/advxp.F

    r5077 r5084  
    126126c test
    127127c  -------------------------------------
    128         DO j =1,jjp1
     128        DO 300 j =1,jjp1
    129129         NUM(j) =1
    130       END DO
     130 300  CONTINUE
    131131c       DO l=1,llm
    132132c      NUM(2,l)=6
    … …  
    150150C  ugri est en kg/s
    151151
    152        DO l = 1,llm
    153        DO j = 1,jjp1
    154        DO i = 1,iip1
     152       DO 500 l = 1,llm
     153       DO 500 j = 1,jjp1
     154       DO 500 i = 1,iip1
    155155       ugri (i,j,llm+1-l) =pbaru (i,j,l)
    156       END DO
    157       END DO
    158       END DO
     156 500   CONTINUE
    159157
    160158C  ---------------------------------------------------------
    … …  
    163161C  boucle principale sur les niveaux et les latitudes
    164162C     
    165       DO L=1,NIV
    166       DO K=lati,latf
     163      DO 1 L=1,NIV
     164      DO 1 K=lati,latf
    167165
    168166C
    … …  
    171169C  program assumes periodic boundaries in X
    172170C
    173       DO I=2,LON
     171      DO 10 I=2,LON
    174172         SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
    175       END DO
     173 10   CONTINUE
    176174      SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
    177175C
    … …  
    181179      LONK=LON/NUMK
    182180C
    183       IF(NUMK>1) THEN
    184 C
    185       DO I=1,LON
     181      IF(NUMK.GT.1) THEN
     182C
     183      DO 111 I=1,LON
    186184         TM(I)=0.
    187       END DO
    188       DO JV=1,NTRA
    189       DO I=1,LON
     185 111  CONTINUE
     186      DO 112 JV=1,NTRA
     187      DO 1120 I=1,LON
    190188         T0 (I,JV)=0.
    191189         TX (I,JV)=0.
    … …  
    198196         TYZ(I,JV)=0.
    199197         TZZ(I,JV)=0.
    200       END DO
    201       END DO
    202 C
    203       DO I2=1,NUMK
    204 C
    205          DO I=1,LONK
     198 1120 CONTINUE
     199 112  CONTINUE
     200C
     201      DO 11 I2=1,NUMK
     202C
     203         DO 113 I=1,LONK
    206204            I3=(I-1)*NUMK+I2
    207205            TM(I)=TM(I)+SM(I3,K,L)
    … …  
    212210            ALF2(I)=ALF1(I)-ALF(I)
    213211            ALF3(I)=ALF(I)*ALF1(I)
    214       END DO
    215 C
    216          DO JV=1,NTRA
    217          DO I=1,LONK
     212 113     CONTINUE
     213C
     214         DO 114 JV=1,NTRA
     215         DO 1140 I=1,LONK
    218216            I3=(I-1)*NUMK+I2
    219217            TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*S0(I3,K,L,JV)
    … …  
    231229            TYZ(I,JV)=TYZ(I,JV)+SYZ(I3,K,L,JV)
    232230            TZZ(I,JV)=TZZ(I,JV)+SZZ(I3,K,L,JV)
    233       END DO
    234       END DO
    235 C
    236       END DO
     231 1140    CONTINUE
     232 114     CONTINUE
     233C
     234 11   CONTINUE
    237235C
    238236      ELSE
    239237C
    240       DO I=1,LON
     238      DO 115 I=1,LON
    241239         TM(I)=SM(I,K,L)
    242       END DO
    243       DO JV=1,NTRA
    244       DO I=1,LON
     240 115  CONTINUE
     241      DO 116 JV=1,NTRA
     242      DO 1160 I=1,LON
    245243         T0 (I,JV)=S0 (I,K,L,JV)
    246244         TX (I,JV)=SSX (I,K,L,JV)
    … …  
    253251         TYZ(I,JV)=SYZ(I,K,L,JV)
    254252         TZZ(I,JV)=SZZ(I,K,L,JV)
    255       END DO
    256       END DO
     253 1160 CONTINUE
     254 116  CONTINUE
    257255C
    258256      ENDIF
    259257C
    260       DO I=1,LONK
     258      DO 117 I=1,LONK
    261259         UEXT(I)=UGRI(I*NUMK,K,L)
    262       END DO
     260 117  CONTINUE
    263261C
    264262C  place limits on appropriate moments before transport
    … …  
    267265      IF(.NOT.LIMIT) GO TO 13
    268266C
    269       DO JV=1,NTRA
    270       DO I=1,LONK
    271         IF(T0(I,JV)>0.) THEN
     267      DO 12 JV=1,NTRA
     268      DO 120 I=1,LONK
     269        IF(T0(I,JV).GT.0.) THEN
    272270          SLPMAX=T0(I,JV)
    273271          S1MAX=1.5*SLPMAX
    … …  
    285283          TXZ(I,JV)=0.
    286284        ENDIF
    287       END DO
    288       END DO
     285 120  CONTINUE
     286 12   CONTINUE
    289287C
    290288 13   CONTINUE
    … …  
    296294C  flux from IP to I if U(I).lt.0
    297295C
    298       DO I=1,LONK-1
    299          IF(UEXT(I)<0.) THEN
     296      DO 140 I=1,LONK-1
     297         IF(UEXT(I).LT.0.) THEN
    300298           FM(I)=-UEXT(I)*DTX
    301299           ALF(I)=FM(I)/TM(I+1)
    302300           TM(I+1)=TM(I+1)-FM(I)
    303301         ENDIF
    304       END DO
     302 140  CONTINUE
    305303C
    306304      I=LONK
    307       IF(UEXT(I)<0.) THEN
     305      IF(UEXT(I).LT.0.) THEN
    308306        FM(I)=-UEXT(I)*DTX
    309307        ALF(I)=FM(I)/TM(1)
    … …  
    313311C  flux from I to IP if U(I).gt.0
    314312C
    315       DO I=1,LONK
    316          IF(UEXT(I)>=0.) THEN
     313      DO 141 I=1,LONK
     314         IF(UEXT(I).GE.0.) THEN
    317315           FM(I)=UEXT(I)*DTX
    318316           ALF(I)=FM(I)/TM(I)
    319317           TM(I)=TM(I)-FM(I)
    320318         ENDIF
    321       END DO
    322 C
    323       DO I=1,LONK
     319 141  CONTINUE
     320C
     321      DO 142 I=1,LONK
    324322         ALFQ(I)=ALF(I)*ALF(I)
    325323         ALF1(I)=1.-ALF(I)
    … …  
    328326         ALF3(I)=ALF(I)*ALFQ(I)
    329327         ALF4(I)=ALF1(I)*ALF1Q(I)
    330       END DO
    331 C
    332       DO JV=1,NTRA
    333       DO I=1,LONK-1
    334 C
    335          IF(UEXT(I)<0.) THEN
     328 142  CONTINUE
     329C
     330      DO 150 JV=1,NTRA
     331      DO 1500 I=1,LONK-1
     332C
     333         IF(UEXT(I).LT.0.) THEN
    336334C
    337335           F0 (I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*
    … …  
    360358         ENDIF
    361359C
    362       END DO
    363       END DO
     360 1500 CONTINUE
     361 150  CONTINUE
    364362C
    365363      I=LONK
    366       IF(UEXT(I)<0.) THEN
    367 C
    368         DO JV=1,NTRA
     364      IF(UEXT(I).LT.0.) THEN
     365C
     366        DO 151 JV=1,NTRA
    369367C
    370368           F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*
    … …  
    391389           TXZ(1,JV)=ALF1Q(I)*TXZ(1,JV)
    392390C
    393       END DO
     391 151    CONTINUE
    394392C
    395393      ENDIF
    396394C
    397       DO JV=1,NTRA
    398       DO I=1,LONK
    399 C
    400          IF(UEXT(I)>=0.) THEN
     395      DO 152 JV=1,NTRA
     396      DO 1520 I=1,LONK
     397C
     398         IF(UEXT(I).GE.0.) THEN
    401399C
    402400           F0 (I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*
    … …  
    425423         ENDIF
    426424C
    427       END DO
    428       END DO
     425 1520 CONTINUE
     426 152  CONTINUE
    429427C
    430428C  puts the temporary moments Fi into appropriate neighboring boxes
    431429C
    432       DO I=1,LONK
    433          IF(UEXT(I)<0.) THEN
     430      DO 160 I=1,LONK
     431         IF(UEXT(I).LT.0.) THEN
    434432           TM(I)=TM(I)+FM(I)
    435433           ALF(I)=FM(I)/TM(I)
    436434         ENDIF
    437       END DO
    438 C
    439       DO I=1,LONK-1
    440          IF(UEXT(I)>=0.) THEN
     435 160  CONTINUE
     436C
     437      DO 161 I=1,LONK-1
     438         IF(UEXT(I).GE.0.) THEN
    441439           TM(I+1)=TM(I+1)+FM(I)
    442440           ALF(I)=FM(I)/TM(I+1)
    443441         ENDIF
    444       END DO
     442 161  CONTINUE
    445443C
    446444      I=LONK
    447       IF(UEXT(I)>=0.) THEN
     445      IF(UEXT(I).GE.0.) THEN
    448446        TM(1)=TM(1)+FM(I)
    449447        ALF(I)=FM(I)/TM(1)
    450448      ENDIF
    451449C
    452       DO I=1,LONK
     450      DO 162 I=1,LONK
    453451         ALF1(I)=1.-ALF(I)
    454452         ALFQ(I)=ALF(I)*ALF(I)
    … …  
    456454         ALF2(I)=ALF1(I)-ALF(I)
    457455         ALF3(I)=ALF(I)*ALF1(I)
    458       END DO
    459 C
    460       DO JV=1,NTRA
    461       DO I=1,LONK
    462 C
    463          IF(UEXT(I)<0.) THEN
     456 162  CONTINUE
     457C
     458      DO 170 JV=1,NTRA
     459      DO 1700 I=1,LONK
     460C
     461         IF(UEXT(I).LT.0.) THEN
    464462C
    465463           TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV)
    … …  
    480478         ENDIF
    481479C
    482       END DO
    483       END DO
    484 C
    485       DO JV=1,NTRA
    486       DO I=1,LONK-1
    487 C
    488          IF(UEXT(I)>=0.) THEN
     480 1700 CONTINUE
     481 170  CONTINUE
     482C
     483      DO 171 JV=1,NTRA
     484      DO 1710 I=1,LONK-1
     485C
     486         IF(UEXT(I).GE.0.) THEN
    489487C
    490488           TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV)
    … …  
    505503         ENDIF
    506504C
    507       END DO
    508       END DO
     505 1710 CONTINUE
     506 171  CONTINUE
    509507C
    510508      I=LONK
    511       IF(UEXT(I)>=0.) THEN
    512         DO JV=1,NTRA
     509      IF(UEXT(I).GE.0.) THEN
     510        DO 172 JV=1,NTRA
    513511           TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV)
    514512           T0 (1,JV)=T0(1,JV)+F0(I,JV)
    … …  
    525523           TYZ(1,JV)=TYZ(1,JV)+FYZ(I,JV)
    526524           TZZ(1,JV)=TZZ(1,JV)+FZZ(I,JV)
    527       END DO
     525 172    CONTINUE
    528526      ENDIF
    529527C
    530528C  retour aux mailles d'origine (passage des Tij aux Sij)
    531529C
    532       IF(NUMK>1) THEN
    533 C
    534       DO I2=1,NUMK
    535 C
    536          DO I=1,LONK
     530      IF(NUMK.GT.1) THEN
     531C
     532      DO 18 I2=1,NUMK
     533C
     534         DO 180 I=1,LONK
    537535C
    538536            I3=I2+(I-1)*NUMK
    … …  
    548546            ALF4(I)=ALF1(I)*ALF1Q(I)
    549547C
    550       END DO
    551 C
    552          DO JV=1,NTRA
    553          DO I=1,LONK
     548 180     CONTINUE
     549C
     550         DO 181 JV=1,NTRA
     551         DO 181 I=1,LONK
    554552C
    555553            I3=I2+(I-1)*NUMK
    … …  
    579577            TXZ(I,JV)=ALF1Q(I)*TXZ(I,JV)
    580578C
    581       END DO
    582       END DO
    583 C
    584       END DO
     579 181     CONTINUE
     580C
     581 18   CONTINUE
    585582C
    586583      ELSE
    587584C
    588       DO I=1,LON
     585      DO 190 I=1,LON
    589586         SM(I,K,L)=TM(I)
    590       END DO
    591       DO JV=1,NTRA
    592       DO I=1,LON
     587 190  CONTINUE
     588      DO 191 JV=1,NTRA
     589      DO 1910 I=1,LON
    593590         S0 (I,K,L,JV)=T0 (I,JV)
    594591         SSX (I,K,L,JV)=TX (I,JV)
    … …  
    601598         SYZ(I,K,L,JV)=TYZ(I,JV)
    602599         SZZ(I,K,L,JV)=TZZ(I,JV)
    603       END DO
    604       END DO
     600 1910 CONTINUE
     601 191  CONTINUE
    605602C
    606603      ENDIF
    607604C
    608       END DO
    609       END DO
     605 1    CONTINUE
    610606C
    611607C ----------- AA Test en fin de ADVX ------ Controle des S*

  • LMDZ6/trunk/libf/dyn3d_common/advy.F

    r5079 r5084  
    121121      enddo
    122122
    123       DO L=1,NIV
     123      DO 1 L=1,NIV
    124124C
    125125C  place limits on appropriate moments before transport
    … …  
    128128      IF(.NOT.LIMIT) GO TO 11
    129129C
    130       DO JV=1,NTRA
    131       DO K=1,LAT
    132       DO I=1,LON
     130      DO 10 JV=1,NTRA
     131      DO 10 K=1,LAT
     132      DO 100 I=1,LON
    133133         sy(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.),
    134134     +                           ABS(sy(I,K,L,JV))),sy(I,K,L,JV))
    135       END DO
    136       END DO
    137       END DO
     135 100  CONTINUE
     136 10   CONTINUE
    138137C
    139138 11   CONTINUE
    … …  
    142141C
    143142      SM0=0.
    144       DO JV=1,NTRA
     143      DO 20 JV=1,NTRA
    145144         S00(JV)=0.
    146       END DO
    147 C
    148       DO I=1,LON
    149 C
    150          IF(VGRI(I,0,L)<=0.) THEN
     145 20   CONTINUE
     146C
     147      DO 21 I=1,LON
     148C
     149         IF(VGRI(I,0,L).LE.0.) THEN
    151150           FM(I,0)=-VGRI(I,0,L)*DTY
    152151           ALF(I,0)=FM(I,0)/SM(I,1,L)
    … …  
    159158         ALF1Q(I,0)=ALF1(I,0)*ALF1(I,0)
    160159C
    161       END DO
    162 C
    163       DO JV=1,NTRA
    164       DO I=1,LON
    165 C
    166          IF(VGRI(I,0,L)<=0.) THEN
     160 21   CONTINUE
     161C
     162      DO 22 JV=1,NTRA
     163      DO 220 I=1,LON
     164C
     165         IF(VGRI(I,0,L).LE.0.) THEN
    167166C
    168167           F0(I,0,JV)=ALF(I,0)*
    … …  
    177176         ENDIF
    178177C
    179       END DO
    180       END DO
    181 C
    182       DO I=1,LON
    183          IF(VGRI(I,0,L)>0.) THEN
     178 220  CONTINUE
     179 22   CONTINUE
     180C
     181      DO 23 I=1,LON
     182         IF(VGRI(I,0,L).GT.0.) THEN
    184183           FM(I,0)=VGRI(I,0,L)*DTY
    185184           ALF(I,0)=FM(I,0)/SM0
    186185         ENDIF
    187       END DO
    188 C
    189       DO JV=1,NTRA
    190       DO I=1,LON
    191          IF(VGRI(I,0,L)>0.) THEN
     186 23   CONTINUE
     187C
     188      DO 24 JV=1,NTRA
     189      DO 240 I=1,LON
     190         IF(VGRI(I,0,L).GT.0.) THEN
    192191           F0(I,0,JV)=ALF(I,0)*S00(JV)
    193192         ENDIF
    194       END DO
    195       END DO
     193 240  CONTINUE
     194 24   CONTINUE
    196195C
    197196C  puts the temporary moments Fi into appropriate neighboring boxes
    198197C
    199       DO I=1,LON
    200 C
    201          IF(VGRI(I,0,L)>0.) THEN
     198      DO 25 I=1,LON
     199C
     200         IF(VGRI(I,0,L).GT.0.) THEN
    202201           SM(I,1,L)=SM(I,1,L)+FM(I,0)
    203202           ALF(I,0)=FM(I,0)/SM(I,1,L)
    … …  
    206205         ALF1(I,0)=1.-ALF(I,0)
    207206C
    208       END DO
    209 C
    210       DO JV=1,NTRA
    211       DO I=1,LON
    212 C
    213          IF(VGRI(I,0,L)>0.) THEN
     207 25   CONTINUE
     208C
     209      DO 26 JV=1,NTRA
     210      DO 260 I=1,LON
     211C
     212         IF(VGRI(I,0,L).GT.0.) THEN
    214213C
    215214         TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV)
    … …  
    219218         ENDIF
    220219C
    221       END DO
    222       END DO
     220 260  CONTINUE
     221 26   CONTINUE
    223222C
    224223C  calculate flux and moments between adjacent boxes
    … …  
    228227C  flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
    229228C
    230       DO K=1,LAT-1
     229      DO 30 K=1,LAT-1
    231230      KP=K+1
    232       DO I=1,LON
    233 C
    234          IF(VGRI(I,K,L)<0.) THEN
     231      DO 300 I=1,LON
     232C
     233         IF(VGRI(I,K,L).LT.0.) THEN
    235234           FM(I,K)=-VGRI(I,K,L)*DTY
    236235           ALF(I,K)=FM(I,K)/SM(I,KP,L)
    … …  
    246245         ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
    247246C
    248       END DO
    249       END DO
    250 C
    251       DO JV=1,NTRA
    252       DO K=1,LAT-1
     247 300  CONTINUE
     248 30   CONTINUE
     249C
     250      DO 31 JV=1,NTRA
     251      DO 31 K=1,LAT-1
    253252      KP=K+1
    254       DO I=1,LON
    255 C
    256          IF(VGRI(I,K,L)<0.) THEN
     253      DO 310 I=1,LON
     254C
     255         IF(VGRI(I,K,L).LT.0.) THEN
    257256C
    258257           F0(I,K,JV)=ALF (I,K)*
    … …  
    282281         ENDIF
    283282C
    284       END DO
    285       END DO
    286       END DO
     283 310  CONTINUE
     284 31   CONTINUE
    287285C
    288286C  puts the temporary moments Fi into appropriate neighboring boxes
    289287C
    290       DO K=1,LAT-1
     288      DO 32 K=1,LAT-1
    291289      KP=K+1
    292       DO I=1,LON
    293 C
    294          IF(VGRI(I,K,L)<0.) THEN
     290      DO 320 I=1,LON
     291C
     292         IF(VGRI(I,K,L).LT.0.) THEN
    295293           SM(I,K,L)=SM(I,K,L)+FM(I,K)
    296294           ALF(I,K)=FM(I,K)/SM(I,K,L)
    … …  
    302300         ALF1(I,K)=1.-ALF(I,K)
    303301C
    304       END DO
    305       END DO
    306 C
    307       DO JV=1,NTRA
    308       DO K=1,LAT-1
     302 320  CONTINUE
     303 32   CONTINUE
     304C
     305      DO 33 JV=1,NTRA
     306      DO 33 K=1,LAT-1
    309307      KP=K+1
    310       DO I=1,LON
    311 C
    312          IF(VGRI(I,K,L)<0.) THEN
     308      DO 330 I=1,LON
     309C
     310         IF(VGRI(I,K,L).LT.0.) THEN
    313311C
    314312         TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
    … …  
    330328         ENDIF
    331329C
    332       END DO
    333       END DO
    334       END DO
     330 330  CONTINUE
     331 33   CONTINUE
    335332C
    336333C  traitement special pour le pole Sud (idem pole Nord)
    … …  
    339336C
    340337      SM0=0.
    341       DO JV=1,NTRA
     338      DO 40 JV=1,NTRA
    342339         S00(JV)=0.
    343       END DO
    344 C
    345       DO I=1,LON
    346 C
    347          IF(VGRI(I,K,L)>=0.) THEN
     340 40   CONTINUE
     341C
     342      DO 41 I=1,LON
     343C
     344         IF(VGRI(I,K,L).GE.0.) THEN
    348345           FM(I,K)=VGRI(I,K,L)*DTY
    349346           ALF(I,K)=FM(I,K)/SM(I,K,L)
    … …  
    356353         ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
    357354C
    358       END DO
    359 C
    360       DO JV=1,NTRA
    361       DO I=1,LON
    362 C
    363          IF(VGRI(I,K,L)>=0.) THEN
     355 41   CONTINUE
     356C
     357      DO 42 JV=1,NTRA
     358      DO 420 I=1,LON
     359C
     360         IF(VGRI(I,K,L).GE.0.) THEN
    364361           F0 (I,K,JV)=ALF(I,K)*
    365362     +                ( S0(I,K,L,JV)+ALF1(I,K)*sy(I,K,L,JV) )
    … …  
    372369         ENDIF
    373370C
    374       END DO
    375       END DO
    376 C
    377       DO I=1,LON
    378          IF(VGRI(I,K,L)<0.) THEN
     371 420  CONTINUE
     372 42   CONTINUE
     373C
     374      DO 43 I=1,LON
     375         IF(VGRI(I,K,L).LT.0.) THEN
    379376           FM(I,K)=-VGRI(I,K,L)*DTY
    380377           ALF(I,K)=FM(I,K)/SM0
    381378         ENDIF
    382       END DO
    383 C
    384       DO JV=1,NTRA
    385       DO I=1,LON
    386          IF(VGRI(I,K,L)<0.) THEN
     379 43   CONTINUE
     380C
     381      DO 44 JV=1,NTRA
     382      DO 440 I=1,LON
     383         IF(VGRI(I,K,L).LT.0.) THEN
    387384           F0(I,K,JV)=ALF(I,K)*S00(JV)
    388385         ENDIF
    389       END DO
    390       END DO
     386 440  CONTINUE
     387 44   CONTINUE
    391388C
    392389C  puts the temporary moments Fi into appropriate neighboring boxes
    393390C
    394       DO I=1,LON
    395 C
    396          IF(VGRI(I,K,L)<0.) THEN
     391      DO 45 I=1,LON
     392C
     393         IF(VGRI(I,K,L).LT.0.) THEN
    397394           SM(I,K,L)=SM(I,K,L)+FM(I,K)
    398395           ALF(I,K)=FM(I,K)/SM(I,K,L)
    … …  
    401398         ALF1(I,K)=1.-ALF(I,K)
    402399C
    403       END DO
    404 C
    405       DO JV=1,NTRA
    406       DO I=1,LON
    407 C
    408          IF(VGRI(I,K,L)<0.) THEN
     400 45   CONTINUE
     401C
     402      DO 46 JV=1,NTRA
     403      DO 460 I=1,LON
     404C
     405         IF(VGRI(I,K,L).LT.0.) THEN
    409406C
    410407         TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
    … …  
    414411         ENDIF
    415412C
    416       END DO
    417       END DO
    418 C
    419       END DO
     413 460  CONTINUE
     414 46   CONTINUE
     415C
     416 1    CONTINUE
    420417C
    421418      RETURN

  • LMDZ6/trunk/libf/dyn3d_common/advyp.F

    r5079 r5084  
    153153C-AA 20/10/94  le signe -1 est necessaire car indexation opposee
    154154
    155       DO l = 1,llm
    156          DO j = 1,jjm
    157             DO i = 1,iip1
     155      DO 500 l = 1,llm
     156         DO 500 j = 1,jjm
     157            DO 500 i = 1,iip1 
    158158            vgri (i,j,llm+1-l)=-1.*pbarv (i,j,l)
    159       END DO
    160       END DO
    161       END DO
     159  500 CONTINUE
    162160
    163161CAA Initialisation de flux fictifs aux bords sup. des boites pol.
    … …  
    173171C  boucle sur les niveaux
    174172C
    175       DO L=1,NIV
     173      DO 1 L=1,NIV
    176174C
    177175C  place limits on appropriate moments before transport
    … …  
    180178      IF(.NOT.LIMIT) GO TO 11
    181179C
    182       DO JV=1,NTRA
    183       DO K=1,LAT
    184       DO I=1,LON
    185          IF(S0(I,K,L,JV)>0.) THEN
     180      DO 10 JV=1,NTRA
     181      DO 10 K=1,LAT
     182      DO 100 I=1,LON
     183         IF(S0(I,K,L,JV).GT.0.) THEN
    186184           SLPMAX=AMAX1(S0(I,K,L,JV),0.)
    187185           S1MAX=1.5*SLPMAX
    … …  
    199197           SYZ(I,K,L,JV)=0.
    200198         ENDIF
    201       END DO
    202       END DO
    203       END DO
     199 100  CONTINUE
     200 10   CONTINUE
    204201C
    205202 11   CONTINUE
    … …  
    208205C
    209206      SM0=0.
    210       DO JV=1,NTRA
     207      DO 20 JV=1,NTRA
    211208         S00(JV)=0.
    212       END DO
    213 C
    214       DO I=1,LON
    215 C
    216          IF(VGRI(I,0,L)<=0.) THEN
     209 20   CONTINUE
     210C
     211      DO 21 I=1,LON
     212C
     213         IF(VGRI(I,0,L).LE.0.) THEN
    217214           FM(I,0)=-VGRI(I,0,L)*DTY
    218215           ALF(I,0)=FM(I,0)/SM(I,1,L)
    … …  
    228225         ALF4(I,0)=ALF1(I,0)*ALF1Q(I,0)
    229226C
    230       END DO
     227 21   CONTINUE
    231228c     print*,'ADVYP 21'
    232229C
    233       DO JV=1,NTRA
    234       DO I=1,LON
    235 C
    236          IF(VGRI(I,0,L)<=0.) THEN
     230      DO 22 JV=1,NTRA
     231      DO 220 I=1,LON
     232C
     233         IF(VGRI(I,0,L).LE.0.) THEN
    237234C
    238235           F0(I,0,JV)=ALF(I,0)* ( S0(I,1,L,JV)-ALF1(I,0)*
    … …  
    256253         ENDIF
    257254C
    258       END DO
    259       END DO
    260 C
    261       DO I=1,LON
    262          IF(VGRI(I,0,L)>0.) THEN
     255 220  CONTINUE
     256 22   CONTINUE
     257C
     258      DO 23 I=1,LON
     259         IF(VGRI(I,0,L).GT.0.) THEN
    263260           FM(I,0)=VGRI(I,0,L)*DTY
    264261           ALF(I,0)=FM(I,0)/SM0
    265262         ENDIF
    266       END DO
    267 C
    268       DO JV=1,NTRA
    269       DO I=1,LON
    270          IF(VGRI(I,0,L)>0.) THEN
     263 23   CONTINUE
     264C
     265      DO 24 JV=1,NTRA
     266      DO 240 I=1,LON
     267         IF(VGRI(I,0,L).GT.0.) THEN
    271268           F0(I,0,JV)=ALF(I,0)*S00(JV)
    272269         ENDIF
    273       END DO
    274       END DO
     270 240  CONTINUE
     271 24   CONTINUE
    275272C
    276273C  puts the temporary moments Fi into appropriate neighboring boxes
    277274C
    278275c     print*,'av ADVYP 25'
    279       DO I=1,LON
    280 C
    281          IF(VGRI(I,0,L)>0.) THEN
     276      DO 25 I=1,LON
     277C
     278         IF(VGRI(I,0,L).GT.0.) THEN
    282279           SM(I,1,L)=SM(I,1,L)+FM(I,0)
    283280           ALF(I,0)=FM(I,0)/SM(I,1,L)
    … …  
    290287         ALF3(I,0)=ALF1(I,0)*ALF(I,0)
    291288C
    292       END DO
     289 25   CONTINUE
    293290c     print*,'av ADVYP 25'
    294291C
    295       DO JV=1,NTRA
    296       DO I=1,LON
    297 C
    298          IF(VGRI(I,0,L)>0.) THEN
     292      DO 26 JV=1,NTRA
     293      DO 260 I=1,LON
     294C
     295         IF(VGRI(I,0,L).GT.0.) THEN
    299296C
    300297         TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV)
    … …  
    308305         ENDIF
    309306C
    310       END DO
    311       END DO
     307 260  CONTINUE
     308 26   CONTINUE
    312309C
    313310C  calculate flux and moments between adjacent boxes
    … …  
    318315C
    319316c     print*,'av ADVYP 30'
    320       DO K=1,LAT-1
     317      DO 30 K=1,LAT-1
    321318      KP=K+1
    322       DO I=1,LON
    323 C
    324          IF(VGRI(I,K,L)<0.) THEN
     319      DO 300 I=1,LON
     320C
     321         IF(VGRI(I,K,L).LT.0.) THEN
    325322           FM(I,K)=-VGRI(I,K,L)*DTY
    326323           ALF(I,K)=FM(I,K)/SM(I,KP,L)
    … …  
    339336         ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K)
    340337C
    341       END DO
    342       END DO
     338 300  CONTINUE
     339 30   CONTINUE
    343340c     print*,'ap ADVYP 30'
    344341C
    345       DO JV=1,NTRA
    346       DO K=1,LAT-1
     342      DO 31 JV=1,NTRA
     343      DO 31 K=1,LAT-1
    347344      KP=K+1
    348       DO I=1,LON
    349 C
    350          IF(VGRI(I,K,L)<0.) THEN
     345      DO 310 I=1,LON
     346C
     347         IF(VGRI(I,K,L).LT.0.) THEN
    351348C
    352349           F0 (I,K,JV)=ALF (I,K)* ( S0(I,KP,L,JV)-ALF1(I,K)*
    … …  
    406403         ENDIF
    407404C
    408       END DO
    409       END DO
    410       END DO
     405 310  CONTINUE
     406 31   CONTINUE
    411407c     print*,'ap ADVYP 31'
    412408C
    413409C  puts the temporary moments Fi into appropriate neighboring boxes
    414410C
    415       DO K=1,LAT-1
     411      DO 32 K=1,LAT-1
    416412      KP=K+1
    417       DO I=1,LON
    418 C
    419          IF(VGRI(I,K,L)<0.) THEN
     413      DO 320 I=1,LON
     414C
     415         IF(VGRI(I,K,L).LT.0.) THEN
    420416           SM(I,K,L)=SM(I,K,L)+FM(I,K)
    421417           ALF(I,K)=FM(I,K)/SM(I,K,L)
    … …  
    431427         ALF3(I,K)=ALF1(I,K)*ALF(I,K)
    432428C
    433       END DO
    434       END DO
     429 320  CONTINUE
     430 32   CONTINUE
    435431c     print*,'ap ADVYP 32'
    436432C
    437       DO JV=1,NTRA
    438       DO K=1,LAT-1
     433      DO 33 JV=1,NTRA
     434      DO 33 K=1,LAT-1
    439435      KP=K+1
    440       DO I=1,LON
    441 C
    442          IF(VGRI(I,K,L)<0.) THEN
     436      DO 330 I=1,LON
     437C
     438         IF(VGRI(I,K,L).LT.0.) THEN
    443439C
    444440         TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
    … …  
    478474         ENDIF
    479475C
    480       END DO
    481       END DO
    482       END DO
     476 330  CONTINUE
     477 33   CONTINUE
    483478c     print*,'ap ADVYP 33'
    484479C
    … …  
    488483C
    489484      SM0=0.
    490       DO JV=1,NTRA
     485      DO 40 JV=1,NTRA
    491486         S00(JV)=0.
    492       END DO
    493 C
    494       DO I=1,LON
    495 C
    496          IF(VGRI(I,K,L)>=0.) THEN
     487 40   CONTINUE
     488C
     489      DO 41 I=1,LON
     490C
     491         IF(VGRI(I,K,L).GE.0.) THEN
    497492           FM(I,K)=VGRI(I,K,L)*DTY
    498493           ALF(I,K)=FM(I,K)/SM(I,K,L)
    … …  
    508503         ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K)
    509504C
    510       END DO
     505 41   CONTINUE
    511506c     print*,'ap ADVYP 41'
    512507C
    513       DO JV=1,NTRA
    514       DO I=1,LON
    515 C
    516          IF(VGRI(I,K,L)>=0.) THEN
     508      DO 42 JV=1,NTRA
     509      DO 420 I=1,LON
     510C
     511         IF(VGRI(I,K,L).GE.0.) THEN
    517512           F0 (I,K,JV)=ALF(I,K)* ( S0(I,K,L,JV)+ALF1(I,K)*
    518513     +             ( SY(I,K,L,JV)+ALF2(I,K)*SYY(I,K,L,JV) ) )
    … …  
    532527         ENDIF
    533528C
    534       END DO
    535       END DO
     529 420  CONTINUE
     530 42   CONTINUE
    536531c     print*,'ap ADVYP 42'
    537532C
    538       DO I=1,LON
    539          IF(VGRI(I,K,L)<0.) THEN
     533      DO 43 I=1,LON
     534         IF(VGRI(I,K,L).LT.0.) THEN
    540535           FM(I,K)=-VGRI(I,K,L)*DTY
    541536           ALF(I,K)=FM(I,K)/SM0
    542537         ENDIF
    543       END DO
     538 43   CONTINUE
    544539c     print*,'ap ADVYP 43'
    545540C
    546       DO JV=1,NTRA
    547       DO I=1,LON
    548          IF(VGRI(I,K,L)<0.) THEN
     541      DO 44 JV=1,NTRA
     542      DO 440 I=1,LON
     543         IF(VGRI(I,K,L).LT.0.) THEN
    549544           F0(I,K,JV)=ALF(I,K)*S00(JV)
    550545         ENDIF
    551       END DO
    552       END DO
     546 440  CONTINUE
     547 44   CONTINUE
    553548C
    554549C  puts the temporary moments Fi into appropriate neighboring boxes
    555550C
    556       DO I=1,LON
    557 C
    558          IF(VGRI(I,K,L)<0.) THEN
     551      DO 45 I=1,LON
     552C
     553         IF(VGRI(I,K,L).LT.0.) THEN
    559554           SM(I,K,L)=SM(I,K,L)+FM(I,K)
    560555           ALF(I,K)=FM(I,K)/SM(I,K,L)
    … …  
    567562         ALF3(I,K)=ALF1(I,K)*ALF(I,K)
    568563C
    569       END DO
     564 45   CONTINUE
    570565c     print*,'ap ADVYP 45'
    571566C
    572       DO JV=1,NTRA
    573       DO I=1,LON
    574 C
    575          IF(VGRI(I,K,L)<0.) THEN
     567      DO 46 JV=1,NTRA
     568      DO 460 I=1,LON
     569C
     570         IF(VGRI(I,K,L).LT.0.) THEN
    576571C
    577572         TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
    … …  
    585580         ENDIF
    586581C
    587       END DO
    588       END DO
     582 460  CONTINUE
     583 46   CONTINUE
    589584c     print*,'ap ADVYP 46'
    590585C
    591       END DO
     586 1    CONTINUE
    592587
    593588c--------------------------------------------------

  • LMDZ6/trunk/libf/dyn3d_common/advz.F

    r5079 r5084  
    117117C  Conversion du flux de masse en kg.s-1
    118118
    119       DO l = 1,llm
    120          DO j = 1,jjp1
    121             DO i = 1,iip1 
     119      DO 500 l = 1,llm
     120         DO 500 j = 1,jjp1
     121            DO 500 i = 1,iip1 
    122122c            wgri (i,j,llm+1-l) =  w (i,j,l) / g
    123123               wgri (i,j,llm+1-l) =  w (i,j,l)
    … …  
    125125c             wgri (i,j,l) = 0.1               !    w (i,j,l)
    126126c             wgri (i,j,llm) = 0.              ! a detruire ult.
    127       END DO
    128       END DO
    129       END DO
     127  500 CONTINUE
    130128         DO  j = 1,jjp1
    131129            DO i = 1,iip1 
    … …  
    139137C  boucle sur les latitudes
    140138C
    141       DO K=1,LAT
     139      DO 1 K=1,LAT
    142140C
    143141C  place limits on appropriate moments before transport
    … …  
    146144      IF(.NOT.LIMIT) GO TO 101
    147145C
    148       DO JV=1,NTRA
    149       DO L=1,NIV
    150          DO I=1,LON
     146      DO 10 JV=1,NTRA
     147      DO 10 L=1,NIV
     148         DO 100 I=1,LON
    151149            sz(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.),
    152150     +                              ABS(sz(I,K,L,JV))),sz(I,K,L,JV))
    153       END DO
    154       END DO
    155       END DO
     151 100     CONTINUE
     152 10   CONTINUE
    156153C
    157154 101  CONTINUE
    … …  
    165162C  2- reajusts moments remaining in the box
    166163C
    167       DO L=1,NIV-1
     164      DO 11 L=1,NIV-1
    168165      LP=L+1
    169166C
    170       DO I=1,LON
    171 C
    172          IF(WGRI(I,K,L)<0.) THEN
     167      DO 110 I=1,LON
     168C
     169         IF(WGRI(I,K,L).LT.0.) THEN
    173170           FM(I,L)=-WGRI(I,K,L)*DTZ
    174171           ALF(I)=FM(I,L)/SM(I,K,LP)
    … …  
    184181         ALF1Q(I)=ALF1(I)*ALF1(I)
    185182C
    186       END DO
    187 C
    188       DO JV=1,NTRA
    189       DO I=1,LON
    190 C
    191          IF(WGRI(I,K,L)<0.) THEN
     183 110  CONTINUE
     184C
     185      DO 111 JV=1,NTRA
     186      DO 1110 I=1,LON
     187C
     188         IF(WGRI(I,K,L).LT.0.) THEN
    192189C
    193190           F0(I,L,JV)=ALF (I)*( S0(I,K,LP,JV)-ALF1(I)*sz(I,K,LP,JV) )
    … …  
    215212         ENDIF
    216213C
    217       END DO
    218       END DO
    219 C
    220       END DO
     214 1110 CONTINUE
     215 111  CONTINUE
     216C
     217 11   CONTINUE
    221218C
    222219C  puts the temporary moments Fi into appropriate neighboring boxes
    223220C
    224       DO L=1,NIV-1
     221      DO 12 L=1,NIV-1
    225222      LP=L+1
    226223C
    227       DO I=1,LON
    228 C
    229          IF(WGRI(I,K,L)<0.) THEN
     224      DO 120 I=1,LON
     225C
     226         IF(WGRI(I,K,L).LT.0.) THEN
    230227           SM(I,K,L)=SM(I,K,L)+FM(I,L)
    231228           ALF(I)=FM(I,L)/SM(I,K,L)
    … …  
    239236         ALF1Q(I)=ALF1(I)*ALF1(I)
    240237C
    241       END DO
    242 C
    243       DO JV=1,NTRA
    244       DO I=1,LON
    245 C
    246          IF(WGRI(I,K,L)<0.) THEN
     238 120  CONTINUE
     239C
     240      DO 121 JV=1,NTRA
     241      DO 1210 I=1,LON
     242C
     243         IF(WGRI(I,K,L).LT.0.) THEN
    247244C
    248245           TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)
    … …  
    263260         ENDIF
    264261C
    265       END DO
    266       END DO
    267 C
    268       END DO
     262 1210 CONTINUE
     263 121  CONTINUE
     264C
     265 12   CONTINUE
    269266C
    270267C  fin de la boucle principale sur les latitudes
    271268C
    272       END DO
     269 1    CONTINUE
    273270C
    274271C-------------------------------------------------------------

  • LMDZ6/trunk/libf/dyn3d_common/advzp.F

    r5079 r5084  
    135135C  Conversion des flux de masses en kg
    136136
    137       DO l = 1,llm
    138          DO j = 1,jjp1
    139             DO i = 1,iip1 
     137      DO 500 l = 1,llm
     138         DO 500 j = 1,jjp1
     139            DO 500 i = 1,iip1 
    140140            wgri (i,j,llm+1-l) = w (i,j,l) 
    141       END DO
    142       END DO
    143       END DO
     141  500 CONTINUE
    144142      do j=1,jjp1
    145143         do i=1,iip1
    … …  
    156154C  boucle sur les latitudes
    157155C
    158       DO K=1,LAT
     156      DO 1 K=1,LAT
    159157C
    160158C  place limits on appropriate moments before transport
    … …  
    163161      IF(.NOT.LIMIT) GO TO 101
    164162C
    165       DO JV=1,NTRA
    166       DO L=1,NIV
    167          DO I=1,LON
    168             IF(S0(I,K,L,JV)>0.) THEN
     163      DO 10 JV=1,NTRA
     164      DO 10 L=1,NIV
     165         DO 100 I=1,LON
     166            IF(S0(I,K,L,JV).GT.0.) THEN
    169167              SLPMAX=S0(I,K,L,JV)
    170168              S1MAX =1.5*SLPMAX
    … …  
    182180              SYZ(I,K,L,JV)=0.
    183181            ENDIF
    184       END DO
    185       END DO
    186       END DO
     182 100     CONTINUE
     183 10   CONTINUE
    187184C
    188185 101  CONTINUE
    … …  
    196193C  2- reajusts moments remaining in the box
    197194C
    198       DO L=1,NIV-1
     195      DO 11 L=1,NIV-1
    199196      LP=L+1
    200197C
    201       DO I=1,LON
    202 C
    203          IF(WGRI(I,K,L)<0.) THEN
     198      DO 110 I=1,LON
     199C
     200         IF(WGRI(I,K,L).LT.0.) THEN
    204201           FM(I,L)=-WGRI(I,K,L)*DTZ
    205202           ALF(I)=FM(I,L)/SM(I,K,LP)
    … …  
    218215         ALF4 (I)=ALF1(I)*ALF1Q(I)
    219216C
    220       END DO
    221 C
    222       DO JV=1,NTRA
    223       DO I=1,LON
    224 C
    225          IF(WGRI(I,K,L)<0.) THEN
     217 110  CONTINUE
     218C
     219      DO 111 JV=1,NTRA
     220      DO 1110 I=1,LON
     221C
     222         IF(WGRI(I,K,L).LT.0.) THEN
    226223C
    227224           F0 (I,L,JV)=ALF (I)* ( S0(I,K,LP,JV)-ALF1(I)*
    … …  
    276273         ENDIF
    277274C
    278       END DO
    279       END DO
    280 C
    281       END DO
     275 1110 CONTINUE
     276 111  CONTINUE
     277C
     278 11   CONTINUE
    282279C
    283280C  puts the temporary moments Fi into appropriate neighboring boxes
    284281C
    285       DO L=1,NIV-1
     282      DO 12 L=1,NIV-1
    286283      LP=L+1
    287284C
    288       DO I=1,LON
    289 C
    290          IF(WGRI(I,K,L)<0.) THEN
     285      DO 120 I=1,LON
     286C
     287         IF(WGRI(I,K,L).LT.0.) THEN
    291288           SM(I,K,L)=SM(I,K,L)+FM(I,L)
    292289           ALF(I)=FM(I,L)/SM(I,K,L)
    … …  
    302299         ALF3(I)=ALF1(I)-ALF(I)
    303300C
    304       END DO
    305 C
    306       DO JV=1,NTRA
    307       DO I=1,LON
    308 C
    309          IF(WGRI(I,K,L)<0.) THEN
     301 120  CONTINUE
     302C
     303      DO 121 JV=1,NTRA
     304      DO 1210 I=1,LON
     305C
     306         IF(WGRI(I,K,L).LT.0.) THEN
    310307C
    311308           TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)
    … …  
    345342         ENDIF
    346343C
    347       END DO
    348       END DO
    349 C
    350       END DO
     344 1210 CONTINUE
     345 121  CONTINUE
     346C
     347 12   CONTINUE
    351348C
    352349C  fin de la boucle principale sur les latitudes
    353350C
    354       END DO
     351 1    CONTINUE
    355352C
    356353      DO l = 1,llm

  • LMDZ6/trunk/libf/dyn3d_common/comdissip.h

    r5077 r5084  
    66
    77      COMMON/comdissip/                                                 &
    8      &    coefdis,tetavel,tetatemp,gamdissip,niterdis
     8     &    niterdis,coefdis,tetavel,tetatemp,gamdissip
    99
    1010

  • LMDZ6/trunk/libf/dyn3d_common/extrapol.F

    r5079 r5084  
    5858200   CONTINUE
    5959      incre = incre + 1
    60       DO j = 1, kylat
    61       DO i = 1, kxlon
    62       IF (pfild(i,j)> zwmsk) THEN
     60      DO 99999 j = 1, kylat
     61      DO 99999 i = 1, kxlon
     62      IF (pfild(i,j).GT. zwmsk) THEN
    6363         pwork(i,j) = pfild(i,j)
    6464         inbor = 0
    … …  
    8989C
    9090C* Correct latitude bounds if southernmost or northernmost points
    91          IF (j == 1) ideb = 4
    92          IF (j == kylat) ifin = 6
     91         IF (j .EQ. 1) ideb = 4
     92         IF (j .EQ. kylat) ifin = 6
    9393C
    9494C* Account for periodicity in longitude
    9595C
    9696         IF (ldper) THEN
    97             IF (i == kxlon) THEN
     97            IF (i .EQ. kxlon) THEN
    9898               ix(3) = 1
    9999               ix(6) = 1
    100100               ix(9) = 1
    101             ELSE IF (i == 1) THEN
     101            ELSE IF (i .EQ. 1) THEN
    102102               ix(1) = kxlon
    103103               ix(4) = kxlon
    … …  
    105105            ENDIF
    106106         ELSE
    107             IF (i == 1) THEN
     107            IF (i .EQ. 1) THEN
    108108               ix(1) = i
    109109               ix(2) = i + 1
    … …  
    113113               ix(6) = i + 1
    114114            ENDIF
    115             IF (i == kxlon) THEN
     115            IF (i .EQ. kxlon) THEN
    116116               ix(1) = i -1
    117117               ix(2) = i
    … …  
    122122            ENDIF
    123123C
    124             IF (i == 1 .OR. i == kxlon) THEN
     124            IF (i .EQ. 1 .OR. i .EQ. kxlon) THEN
    125125               jy(1) = MAX (1,j-1)
    126126               jy(2) = MAX (1,j-1)
    … …  
    132132               ideb = 1
    133133               ifin = 6
    134                IF (j == 1) ideb = 3
    135                IF (j == kylat) ifin = 4
     134               IF (j .EQ. 1) ideb = 3
     135               IF (j .EQ. kylat) ifin = 4
    136136            ENDIF
    137137         ENDIF ! end for ldper test
    … …  
    139139C* Find unmasked neighbors
    140140C
    141          DO k = ideb, ifin
     141         DO 230 k = ideb, ifin
    142142            zmask(k) = 0.
    143143            ilon = ix(k)
    144144            jlat = jy(k)
    145             IF (pfild(ilon,jlat) < zwmsk) THEN
     145            IF (pfild(ilon,jlat) .LT. zwmsk) THEN
    146146               zmask(k) = 1.
    147147               inbor = inbor + 1
    148148            ENDIF
    149       END DO
     149 230     CONTINUE
    150150C
    151151C* Not enough points around point P are unmasked; interpolation on P
    152152C  will be done in a future call to extrap.
    153153C
    154          IF (inbor >= knbor) THEN
     154         IF (inbor .GE. knbor) THEN
    155155            pwork(i,j) = 0.
    156156            DO k = ideb, ifin
    … …  
    163163C
    164164      ENDIF
    165       END DO
    166       END DO
     16599999 CONTINUE
    167166C
    168167C*    3. Writing back unmasked field in pfild
    … …  
    177176      DO j = 1, kylat
    178177      DO i = 1, kxlon
    179          IF (pwork(i,j) > zwmsk) idoit = idoit + 1
     178         IF (pwork(i,j) .GT. zwmsk) idoit = idoit + 1
    180179         pfild(i,j) = pwork(i,j)
    181180      ENDDO
    182181      ENDDO
    183182c
    184       IF (idoit /= 0) GOTO 200
     183      IF (idoit .ne. 0) GOTO 200
    185184ccc      PRINT*, "Number of extrapolation steps incre =", incre
    186185c

  • LMDZ6/trunk/libf/dyn3d_common/grilles_gcm_netcdf_sub.F90

    r5075 r5084  
    1414  USE comconst_mod, ONLY: cpp, kappa, g, omeg, daysec, rad, pi
    1515  USE comvert_mod, ONLY: presnivs, preff, pa
    16   USE lmdz_netcdf, ONLY: nf90_def_var, nf90_int, nf90_float, nf90_put_var, nf_enddef, &
    17       nf_put_att_text,nf_def_dim,nf_64bit_offset,nf_clobber,nf_create
     16  use netcdf, only: nf90_def_var, nf90_int, nf90_float, nf90_put_var
    1817 
    1918  IMPLICIT NONE
    … …  
    2221  INCLUDE "paramet.h"
    2322  INCLUDE "comgeom.h"
     23  INCLUDE "netcdf.inc"
    2424
    2525!========================
    … …  
    232232
    233233SUBROUTINE handle_err(status)
    234   USE lmdz_netcdf, ONLY: nf_strerror
     234  INCLUDE "netcdf.inc"
    235235
    236236  INTEGER status
    237   IF (status/=nf_noerr) THEN
     237  IF (status.NE.nf_noerr) THEN
    238238     PRINT *,NF_STRERROR(status)
    239239     CALL abort_gcm('grilles_gcm_netcdf','netcdf error',1)

  • LMDZ6/trunk/libf/dyn3d_common/ppm3d.F

    r5079 r5084  
    6868      implicit none
    6969
    70 c     rajout de d�clarations
     70c     rajout de déclarations
    7171c      integer Jmax,kmax,ndt0,nstep,k,j,i,ic,l,js,jn,imh,iad,jad,krd
    7272c      integer iu,iiu,j2,jmr,js0,jt
    … …  
    315315C
    316316C *********** Initialization **********************
    317       if(NSTEP==1) then
     317      if(NSTEP.eq.1) then
    318318c
    319319      write(6,*) '------------------------------------ '
    … …  
    325325C
    326326C controles sur les parametres
    327       if(NLAY<6) then
     327      if(NLAY.LT.6) then
    328328        write(6,*) 'NLAY must be >= 6'
    329329        stop
    330330      endif
    331       if (JNP<NLAY) then
     331      if (JNP.LT.NLAY) then
    332332         write(6,*) 'JNP must be >= NLAY'
    333333        stop
    334334      endif
    335335      IMRD2=mod(IMR,2)
    336       if (j1==2.and.IMRD2/=0) then
     336      if (j1.eq.2.and.IMRD2.NE.0) then
    337337         write(6,*) 'if j1=2 IMR must be an even integer'
    338338        stop
    … …  
    340340
    341341C
    342       if(Jmax<JNP .or. Kmax<NLAY) then
     342      if(Jmax.lt.JNP .or. Kmax.lt.NLAY) then
    343343        write(6,*) 'Jmax or Kmax is too small'
    344344        stop
    … …  
    354354      DP =    PI / REAL(JMR)
    355355C
    356       if(IGD==0) then
     356      if(IGD.eq.0) then
    357357C Compute analytic cosine at cell edges
    358358            call cosa(cosp,cose,JNP,PI,DP)
    … …  
    362362      endif
    363363C
    364       do J=2,JMR
    365       acosp(j) = 1. / cosp(j)
    366       END DO
     364      do 15 J=2,JMR
     36515    acosp(j) = 1. / cosp(j)
    367366C
    368367C Inverse of the Scaled polar cap area.
    … …  
    373372      endif
    374373C
    375       if(NDT0 /= NDT) then
     374      if(NDT0 .ne. NDT) then
    376375      DT   = NDT
    377376      NDT0 = NDT
    378377
    379         if(Umax < 180.) then
     378        if(Umax .lt. 180.) then
    380379         write(6,*) 'Umax may be too small!'
    381380        endif
    … …  
    383382      MaxDT = DP*AE / abs(Umax) + 0.5
    384383      write(6,*)'Largest time step for max(V)=',Umax,' is ',MaxDT
    385       if(MaxDT < abs(NDT)) then
     384      if(MaxDT .lt. abs(NDT)) then
    386385            write(6,*) 'Warning!!! NDT maybe too large!'
    387386      endif
    388387C
    389       if(CR1>=0.95) then
     388      if(CR1.ge.0.95) then
    390389      JS0 = 0
    391390      JN0 = 0
    … …  
    430429         
    431430C
    432       if(j1/=2) then
    433       DO IC=1,NC
    434       DO L=1,NLAY
    435       DO I=1,IMR
     431      if(j1.ne.2) then
     432      DO 40 IC=1,NC
     433      DO 40 L=1,NLAY
     434      DO 40 I=1,IMR
    436435      Q(I,  2,L,IC) = Q(I,  1,L,IC)
    437       Q(I,JMR,L,IC) = Q(I,JNP,L,IC)
    438       END DO
    439       END DO
    440       END DO
     43640    Q(I,JMR,L,IC) = Q(I,JNP,L,IC)
    441437      endif
    442438C
    443439C Compute "tracer density"
    444       DO IC=1,NC
    445       DO k=1,NLAY
    446       DO j=1,JNP
    447       DO i=1,IMR
    448       DQ(i,j,k,IC) = Q(i,j,k,IC)*delp1(i,j,k)
    449       END DO
    450       END DO
    451       END DO
    452       END DO
    453 C
    454       do k=1,NLAY
    455 C
    456       if(IGD==0) then
     440      DO 550 IC=1,NC
     441      DO 44 k=1,NLAY
     442      DO 44 j=1,JNP
     443      DO 44 i=1,IMR
     44444    DQ(i,j,k,IC) = Q(i,j,k,IC)*delp1(i,j,k)
     445550     continue
     446C
     447      do 1500 k=1,NLAY
     448C
     449      if(IGD.eq.0) then
    457450C Convert winds on A-Grid to Courant # on C-Grid.
    458451      call A2C(U(1,1,k),V(1,1,k),IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)
    459452      else
    460453C Convert winds on C-grid to Courant #
    461       do j=j1,j2
    462       do i=2,IMR
    463       CRX(i,J) = dtdx(j)*U(i-1,j,k)
    464       END DO
    465       END DO
     454      do 45 j=j1,j2
     455      do 45 i=2,IMR
     45645    CRX(i,J) = dtdx(j)*U(i-1,j,k)
    466457   
    467458C
    468       do j=j1,j2
    469       CRX(1,J) = dtdx(j)*U(IMR,j,k)
    470       END DO
    471 C
    472       do i=1,IMR*JMR
    473       CRY(i,2) = DTDY*V(i,1,k)
    474       END DO
     459      do 50 j=j1,j2
     46050    CRX(1,J) = dtdx(j)*U(IMR,j,k)
     461C
     462      do 55 i=1,IMR*JMR
     46355    CRY(i,2) = DTDY*V(i,1,k)
    475464      endif
    476465C     
    … …  
    481470      do j=JS0,j1+1,-1
    482471      do i=1,IMR
    483       if(abs(CRX(i,j))>1.) then
     472      if(abs(CRX(i,j)).GT.1.) then
    484473            JS = j
    485474            go to 2222
    … …  
    491480      do j=JN0,j2-1
    492481      do i=1,IMR
    493       if(abs(CRX(i,j))>1.) then
     482      if(abs(CRX(i,j)).GT.1.) then
    494483            JN = j
    495484            go to 2233
    … …  
    4994882233  continue
    500489C
    501       if(j1/=2) then           ! Enlarged polar cap.
     490      if(j1.ne.2) then           ! Enlarged polar cap.
    502491      do i=1,IMR
    503492      DPI(i,  2,k) = 0.
    … …  
    516505      enddo
    517506C
    518       do j=j1,j2
    519       DO i=1,IMR
    520       DPI(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j)
    521       END DO
    522       END DO
     507      do 95 j=j1,j2
     508      DO 95 i=1,IMR
     50995    DPI(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j)
    523510C
    524511C Poles
    … …  
    549536      enddo
    550537C
    551       do j=j1,j2
    552       DO i=1,IMR
    553       xmass(i,j) = PU(i,j)*CRX(i,j)
    554       END DO
    555       END DO
    556 C
    557       DO j=j1,j2
    558       DO i=1,IMR-1
    559       DPI(i,j,k) = DPI(i,j,k) + xmass(i,j) - xmass(i+1,j)
    560       END DO
    561       END DO
    562 C
    563       DO j=j1,j2
    564       DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j)
    565       END DO
     538      do 110 j=j1,j2
     539      DO 110 i=1,IMR
     540110   xmass(i,j) = PU(i,j)*CRX(i,j)
     541C
     542      DO 120 j=j1,j2
     543      DO 120 i=1,IMR-1
     544120   DPI(i,j,k) = DPI(i,j,k) + xmass(i,j) - xmass(i+1,j)
     545C
     546      DO 130 j=j1,j2
     547130   DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j)
    566548C
    567549      DO j=j1,j2
    … …  
    587569      enddo
    588570C
    589       if(j1==2) then
     571      if(j1.eq.2) then
    590572        IMH = IMR/2
    591573      do i=1,IMH
    … …  
    600582C
    601583C ****6***0*********0*********0*********0*********0*********0**********72
    602       do IC=1,NC
     584      do 1000 IC=1,NC
    603585C
    604586      do i=1,IMJM
    … …  
    608590C
    609591C E-W advective cross term
    610       do j=J1,J2
    611       if(J>JS  .and. J<JN) GO TO 250
     592      do 250 j=J1,J2
     593      if(J.GT.JS  .and. J.LT.JN) GO TO 250
    612594C
    613595      do i=1,IMR
    … …  
    620602      enddo
    621603C
    622       DO i=1,IMR
     604      DO 230 i=1,IMR
    623605      iu = UA(i,j)
    624606      ru = UA(i,j) - iu
    625607      iiu = i-iu
    626       if(UA(i,j)>=0.) then
     608      if(UA(i,j).GE.0.) then
    627609      wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))
    628610      else
    … …  
    630612      endif
    631613      wk1(i,j,1) = wk1(i,j,1) - qtmp(i)
    632       END DO
     614230   continue
    633615250   continue
    634       END DO
    635 C
    636       if(JN/=0) then
     616C
     617      if(JN.ne.0) then
    637618      do j=JS+1,JN-1
    638619C
    … …  
    664645        if(cross) then
    665646C Add cross terms in the vertical direction.
    666         if(IORD >= 2) then
     647        if(IORD .GE. 2) then
    667648                iad = 2
    668649        else
    … …  
    670651        endif
    671652C
    672         if(JORD >= 2) then
     653        if(JORD .GE. 2) then
    673654                jad = 2
    674655        else
    … …  
    690671     &  DC2,ymass,WK1(1,1,3),wk1(1,1,4),WK1(1,1,5),WK1(1,1,6),JORD)
    691672C
    692       END DO
    693       END DO
     6731000  continue
     6741500  continue
    694675C
    695676C ******* Compute vertical mass flux (same unit as PS) ***********
    … …  
    697678C 1st step: compute total column mass CONVERGENCE.
    698679C
    699       do j=1,JNP
    700       do i=1,IMR
    701       CRY(i,j) = DPI(i,j,1)
    702       END DO
    703       END DO
    704 C
    705       do k=2,NLAY
    706       do j=1,JNP
    707       do i=1,IMR
     680      do 320 j=1,JNP
     681      do 320 i=1,IMR
     682320   CRY(i,j) = DPI(i,j,1)
     683C
     684      do 330 k=2,NLAY
     685      do 330 j=1,JNP
     686      do 330 i=1,IMR
    708687      CRY(i,j)  = CRY(i,j) + DPI(i,j,k)
    709       END DO
    710       END DO
    711       END DO
    712 C
    713       do j=1,JNP
    714       do i=1,IMR
     688330   continue
     689C
     690      do 360 j=1,JNP
     691      do 360 i=1,IMR
    715692C
    716693C 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.
    … …  
    723700      W(i,j,1) = DPI(i,j,1) - DBK(1)*CRY(i,j)
    724701      W(i,j,NLAY) = 0.
    725       END DO
    726       END DO
    727 C
    728       do k=2,NLAY-1
    729       do j=1,JNP
    730       do i=1,IMR
     702360   continue
     703C
     704      do 370 k=2,NLAY-1
     705      do 370 j=1,JNP
     706      do 370 i=1,IMR
    731707      W(i,j,k) = W(i,j,k-1) + DPI(i,j,k) - DBK(k)*CRY(i,j)
    732       END DO
    733       END DO
    734       END DO
    735 C
    736       DO k=1,NLAY
    737       DO j=1,JNP
    738       DO i=1,IMR
     708370   continue
     709C
     710      DO 380 k=1,NLAY
     711      DO 380 j=1,JNP
     712      DO 380 i=1,IMR
    739713      delp2(i,j,k) = DAP(k) + DBK(k)*PS2(i,j)
    740       END DO
    741       END DO
    742       END DO
     714380   continue
    743715C
    744716        KRD = max(3, KORD)
    745       do IC=1,NC
     717      do 4000 IC=1,NC
    746718C
    747719C****6***0*********0*********0*********0*********0*********0**********72
    … …  
    766738      enddo
    767739C     
    768       if(j1/=2) then
    769       DO k=1,NLAY
    770       DO I=1,IMR
    771 c     j=1 c'est le p�le Sud, j=JNP c'est le p�le Nord
     740      if(j1.ne.2) then
     741      DO 400 k=1,NLAY
     742      DO 400 I=1,IMR
     743c     j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord
    772744      Q(I,  2,k,IC) = Q(I,  1,k,IC)
    773745      Q(I,JMR,k,IC) = Q(I,JNP,k,IC)
    774       END DO
    775       END DO
    776       endif
    777       END DO
    778 C
    779       if(j1/=2) then
    780       DO k=1,NLAY
    781       DO i=1,IMR
     746400   CONTINUE
     747      endif
     7484000  continue
     749C
     750      if(j1.ne.2) then
     751      DO 5000 k=1,NLAY
     752      DO 5000 i=1,IMR
    782753      W(i,  2,k) = W(i,  1,k)
    783754      W(i,JMR,k) = W(i,JNP,k)
    784       END DO
    785       END DO
     7555000  continue
    786756      endif
    787757C
    … …  
    815785C ****6***0*********0*********0*********0*********0*********0**********72
    816786C
    817       do k=1,NLAYM1
    818       do i=1,IMJM
     787      do 1000 k=1,NLAYM1
     788      do 1000 i=1,IMJM
    819789      DQDT(i,1,k) = P(i,1,k+1) - P(i,1,k)
    820       END DO
    821       END DO
    822 C
    823       DO k=2,NLAYM1
    824       DO I=1,IMJM
     7901000  continue
     791C
     792      DO 1220 k=2,NLAYM1
     793      DO 1220 I=1,IMJM   
    825794       c0 =  delp(i,1,k) / (delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1))
    826795       c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k))   
    … …  
    830799      Qmin = P(i,1,k) - min(P(i,1,k-1),P(i,1,k),P(i,1,k+1))
    831800      DC(i,1,k) = sign(min(abs(tmp),Qmax,Qmin), tmp)   
    832       END DO
    833       END DO
     8011220  CONTINUE
    834802     
    835803C     
    … …  
    838806C ****6***0*********0*********0*********0*********0*********0**********72
    839807C
    840       DO j=1,JNP
    841       if((j==2 .or. j==JMR) .and. j1/=2) goto 2000
     808      DO 2000 j=1,JNP
     809      if((j.eq.2 .or. j.eq.JMR) .and. j1.ne.2) goto 2000
    842810C
    843811      DO k=1,NLAY
    … …  
    860828C
    861829C First guess top edge value
    862       DO i=1,IMR
     830      DO 10 i=1,IMR
    863831C three-cell PPM
    864832C Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp
    … …  
    872840C
    873841C Check if change sign
    874       if(wk1(i,1)*AL(i,1)<=0.) then
     842      if(wk1(i,1)*AL(i,1).le.0.) then
    875843                 AL(i,1) = 0.
    876844             flux(i,1) = 0.
    … …  
    878846             flux(i,1) =  wk1(i,1) - AL(i,1)
    879847        endif
    880       END DO
     84810    continue
    881849C
    882850C Bottom
    883       DO i=1,IMR
     851      DO 15 i=1,IMR
    884852C 2-cell PPM with zero gradient right at the surface
    885853C
    … …  
    888856      AR(i,NLAY) = wk1(i,NLAY) + fct
    889857      AL(i,NLAY) = wk1(i,NLAY) - (fct+fct)
    890       if(wk1(i,NLAY)*AR(i,NLAY)<=0.) AR(i,NLAY) = 0.
     858      if(wk1(i,NLAY)*AR(i,NLAY).le.0.) AR(i,NLAY) = 0.
    891859      flux(i,NLAY) = AR(i,NLAY) -  wk1(i,NLAY)
    892       END DO
     86015    continue
    893861     
    894862C
    … …  
    897865C****6***0*********0*********0*********0*********0*********0**********72
    898866C
    899       DO k=3,NLAYM1
    900       DO i=1,IMR
     867      DO 14 k=3,NLAYM1
     868      DO 12 i=1,IMR
    901869      c1 =  DQDT(i,j,k-1)*wk2(i,k-1) / (wk2(i,k-1)+wk2(i,k))
    902870      c2 =  2. / (wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1))
    … …  
    907875     &          wk2(i,k-1)*A1*flux(i,k)  )
    908876C      print *,'AL1',i,k, AL(i,k)
    909       END DO
    910       END DO
    911 C
    912       do i=1,IMR*NLAYM1
     87712    CONTINUE
     87814    continue
     879C
     880      do 20 i=1,IMR*NLAYM1
    913881      AR(i,1) = AL(i,2)
    914882C      print *,'AR1',i,AR(i,1)
    915       END DO
    916 C
    917       do i=1,IMR*NLAY
     88320    continue
     884C
     885      do 30 i=1,IMR*NLAY
    918886      A6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (AL(i,1)+AR(i,1)))
    919887C      print *,'A61',i,A6(i,1)
    920       END DO
     88830    continue
    921889C
    922890C****6***0*********0*********0*********0*********0*********0**********72
    … …  
    927895C
    928896C Interior depending on KORD
    929       if(LMT<=2)
     897      if(LMT.LE.2)
    930898     &  call lmtppm(flux(1,2),A6(1,2),AR(1,2),AL(1,2),wk1(1,2),
    931899     &              IMR*(NLAY-2),LMT)
    … …  
    933901C****6***0*********0*********0*********0*********0*********0**********72
    934902C
    935       DO i=1,IMR*NLAYM1
    936       IF(wz2(i,1)>0.) then
     903      DO 140 i=1,IMR*NLAYM1
     904      IF(wz2(i,1).GT.0.) then
    937905        CM = wz2(i,1) / wk2(i,1)
    938906        flux(i,2) = AR(i,1)+0.5*CM*(AL(i,1)-AR(i,1)+A6(i,1)*(1.-R23*CM))
    … …  
    944912C        print *,'test2',i, AL(i,2),AR(i,2),A6(i,2),R23
    945913      endif
    946       END DO
    947 C
    948       DO i=1,IMR*NLAYM1
     914140   continue
     915C
     916      DO 250 i=1,IMR*NLAYM1
    949917      flux(i,2) = wz2(i,1) * flux(i,2)
    950       END DO
    951 C
    952       do i=1,IMR
     918250   continue
     919C
     920      do 350 i=1,IMR
    953921      DQ(i,j,   1) = DQ(i,j,   1) - flux(i,   2)
    954922      DQ(i,j,NLAY) = DQ(i,j,NLAY) + flux(i,NLAY)
    955       END DO
    956 C
    957       do k=2,NLAYM1
    958       do i=1,IMR
    959       DQ(i,j,k) = DQ(i,j,k) + flux(i,k) - flux(i,k+1)
    960       END DO
    961       END DO
     923350   continue
     924C
     925      do 360 k=2,NLAYM1
     926      do 360 i=1,IMR
     927360   DQ(i,j,k) = DQ(i,j,k) + flux(i,k) - flux(i,k+1)
    9629282000  continue
    963       END DO
    964929      return
    965930      end
    … …  
    985950      j2vl = j2-jvan
    986951C
    987       do j=j1,j2
     952      do 1310 j=j1,j2
    988953C
    989954      do i=1,IMR
    … …  
    991956      enddo
    992957C
    993       if(j>=JN .or. j<=JS) goto 2222
     958      if(j.ge.JN .or. j.le.JS) goto 2222
    994959C ************* Eulerian **********
    995960C
    … …  
    999964      qtmp(IMP+1) = q(2,J)
    1000965C
    1001       IF(IORD==1 .or. j==j1 .or. j==j2) THEN
    1002       DO i=1,IMR
     966      IF(IORD.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN
     967      DO 1406 i=1,IMR
    1003968      iu = REAL(i) - uc(i,j)
    1004       fx1(i) = qtmp(iu)
    1005       END DO
     9691406  fx1(i) = qtmp(iu)
    1006970      ELSE
    1007971      call xmist(IMR,IML,Qtmp,DC)
    1008972      DC(0) = DC(IMR)
    1009973C
    1010       if(IORD==2 .or. j<=j1vl .or. j>=j2vl) then
    1011       DO i=1,IMR
     974      if(IORD.eq.2 .or. j.le.j1vl .or. j.ge.j2vl) then
     975      DO 1408 i=1,IMR
    1012976      iu = REAL(i) - uc(i,j)
    1013       fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j))
    1014       END DO
     9771408  fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j))
    1015978      else
    1016979      call fxppm(IMR,IML,UC(1,j),Qtmp,DC,fx1,IORD)
    … …  
    1019982      ENDIF
    1020983C
    1021       DO i=1,IMR
    1022       fx1(i) = fx1(i)*xmass(i,j)
    1023       END DO
     984      DO 1506 i=1,IMR
     9851506  fx1(i) = fx1(i)*xmass(i,j)
    1024986C
    1025987      goto 1309
    … …  
    1034996      enddo
    1035997C
    1036       IF(IORD==1 .or. j==j1 .or. j==j2) THEN
    1037       DO i=1,IMR
     998      IF(IORD.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN
     999      DO 1306 i=1,IMR
    10381000      itmp = INT(uc(i,j))
    10391001      ISAVE(i) = i - itmp
    10401002      iu = i - uc(i,j)
    1041       fx1(i) = (uc(i,j) - itmp)*qtmp(iu)
    1042       END DO
     10031306  fx1(i) = (uc(i,j) - itmp)*qtmp(iu)
    10431004      ELSE
    10441005      call xmist(IMR,IML,Qtmp,DC)
    … …  
    10491010      enddo
    10501011C
    1051       DO i=1,IMR
     1012      DO 1307 i=1,IMR
    10521013      itmp = INT(uc(i,j))
    10531014      rut  = uc(i,j) - itmp
    10541015      ISAVE(i) = i - itmp
    10551016      iu = i - uc(i,j)
    1056       fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut))
    1057       END DO
     10171307  fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut))
    10581018      ENDIF
    10591019C
    1060       do i=1,IMR
    1061       IF(uc(i,j)>1.) then
     1020      do 1308 i=1,IMR
     1021      IF(uc(i,j).GT.1.) then
    10621022CDIR$ NOVECTOR
    10631023        do ist = ISAVE(i),i-1
    10641024        fx1(i) = fx1(i) + qtmp(ist)
    10651025        enddo
    1066       elseIF(uc(i,j)<-1.) then
     1026      elseIF(uc(i,j).LT.-1.) then
    10671027        do ist = i,ISAVE(i)-1
    10681028        fx1(i) = fx1(i) - qtmp(ist)
    … …  
    10701030CDIR$ VECTOR
    10711031      endif
    1072       END DO
     10321308  continue
    10731033      do i=1,IMR
    10741034      fx1(i) = PU(i,j)*fx1(i)
    … …  
    10781038C
    107910391309  fx1(IMP) = fx1(1)
    1080       DO i=1,IMR
    1081       DQ(i,j) =  DQ(i,j) + fx1(i)-fx1(i+1)
    1082       END DO
     1040      DO 1215 i=1,IMR
     10411215  DQ(i,j) =  DQ(i,j) + fx1(i)-fx1(i+1)
    10831042C
    10841043C ***************************************
    10851044C
    1086       END DO
     10451310  continue
    10871046      return
    10881047      end
    … …  
    11201079C      endif
    11211080C
    1122       DO i=1,IMR
    1123       AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3
    1124       END DO
    1125 C
    1126       do i=1,IMR-1
    1127       AR(i) = AL(i+1)
    1128       END DO
     1081      DO 10 i=1,IMR
     108210    AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3
     1083C
     1084      do 20 i=1,IMR-1
     108520    AR(i) = AL(i+1)
    11291086      AR(IMR) = AL(1)
    11301087C
    1131       do i=1,IMR
    1132       A6(i) = 3.*(p(i)+p(i)  - (AL(i)+AR(i)))
    1133       END DO
    1134 C
    1135       if(LMT<=2) call lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT)
     1088      do 30 i=1,IMR
     108930    A6(i) = 3.*(p(i)+p(i)  - (AL(i)+AR(i)))
     1090C
     1091      if(LMT.LE.2) call lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT)
    11361092C
    11371093      AL(0) = AL(IMR)
    … …  
    11401096C
    11411097      DO i=1,IMR
    1142       IF(UT(i)>0.) then
     1098      IF(UT(i).GT.0.) then
    11431099      flux(i) = AR(i-1) + 0.5*UT(i)*(AL(i-1) - AR(i-1) +
    11441100     &                 A6(i-1)*(1.-R23*UT(i)) )
    … …  
    11591115      real :: tmp,pmax,pmin
    11601116C
    1161       do i=1,IMR
     1117      do 10  i=1,IMR
    11621118      tmp = R24*(8.*(p(i+1) - p(i-1)) + p(i-2) - p(i+2))
    11631119      Pmax = max(P(i-1), p(i), p(i+1)) - p(i)
    11641120      Pmin = p(i) - min(P(i-1), p(i), p(i+1))
    1165       DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp)
    1166       END DO
     112110    DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp)
    11671122      return
    11681123      end
    … …  
    11831138      len = IMR*(J2-J1+2)
    11841139C
    1185       if(JORD==1) then
    1186       DO i=1,len
     1140      if(JORD.eq.1) then
     1141      DO 1000 i=1,len
    11871142      JT = REAL(J1) - VC(i,J1)
    1188       fx(i,j1) = p(i,JT)
    1189       END DO
     11431000  fx(i,j1) = p(i,JT)
    11901144      else
    11911145   
    11921146      call ymist(IMR,JNP,j1,P,DC2,4)
    11931147C
    1194       if(JORD<=0 .or. JORD>=3) then
     1148      if(JORD.LE.0 .or. JORD.GE.3) then
    11951149   
    11961150      call fyppm(VC,P,DC2,fx,IMR,JNP,j1,j2,A6,AR,AL,JORD)
    11971151   
    11981152      else
    1199       DO i=1,len
     1153      DO 1200 i=1,len
    12001154      JT = REAL(J1) - VC(i,J1)
    1201       fx(i,j1) = p(i,JT) + (sign(1.,VC(i,j1))-VC(i,j1))*DC2(i,JT)
    1202       END DO
    1203       endif
    1204       endif
    1205 C
    1206       DO i=1,len
    1207       fx(i,j1) = fx(i,j1)*ymass(i,j1)
    1208       END DO
    1209 C
    1210       DO j=j1,j2
    1211       DO i=1,IMR
    1212       DQ(i,j) = DQ(i,j) + (fx(i,j) - fx(i,j+1)) * acosp(j)
    1213       END DO
    1214       END DO
     11551200  fx(i,j1) = p(i,JT) + (sign(1.,VC(i,j1))-VC(i,j1))*DC2(i,JT)
     1156      endif
     1157      endif
     1158C
     1159      DO 1300 i=1,len
     11601300  fx(i,j1) = fx(i,j1)*ymass(i,j1)
     1161C
     1162      DO 1400 j=j1,j2
     1163      DO 1400 i=1,IMR
     11641400  DQ(i,j) = DQ(i,j) + (fx(i,j) - fx(i,j+1)) * acosp(j)
    12151165C
    12161166C Poles
    … …  
    12291179      enddo
    12301180C
    1231       if(j1/=2) then
     1181      if(j1.ne.2) then
    12321182      do i=1,IMR
    12331183      DQ(i,  2) = sum1
    … …  
    12511201      IJM3 = IMR*(JMR-3)
    12521202C
    1253       IF(ID==2) THEN
    1254       do i=1,IMR*(JMR-1)
     1203      IF(ID.EQ.2) THEN
     1204      do 10 i=1,IMR*(JMR-1)
    12551205      tmp = 0.25*(p(i,3) - p(i,1))
    12561206      Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2)
    12571207      Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3))
    12581208      DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp)
    1259       END DO
     120910    CONTINUE
    12601210      ELSE
    1261       do i=1,IMH
     1211      do 12 i=1,IMH
    12621212C J=2
    12631213      tmp = (8.*(p(i,3) - p(i,1)) + p(i+IMH,2) - p(i,4))*R24
    … …  
    12701220      Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP))
    12711221      DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp)
    1272       END DO
    1273       do i=IMH+1,IMR
     122212    CONTINUE
     1223      do 14 i=IMH+1,IMR
    12741224C J=2
    12751225      tmp = (8.*(p(i,3) - p(i,1)) + p(i-IMH,2) - p(i,4))*R24
    … …  
    12821232      Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP))
    12831233      DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp)
    1284       END DO
    1285 C
    1286       do i=1,IJM3
     123414    CONTINUE
     1235C
     1236      do 15 i=1,IJM3
    12871237      tmp = (8.*(p(i,4) - p(i,2)) + p(i,1) - p(i,5))*R24
    12881238      Pmax = max(p(i,2),p(i,3),p(i,4)) - p(i,3)
    12891239      Pmin = p(i,3) - min(p(i,2),p(i,3),p(i,4))
    12901240      DC(i,3) = sign(min(abs(tmp),Pmin,Pmax),tmp)
    1291       END DO
     124115    CONTINUE
    12921242      ENDIF
    12931243C
    1294       if(j1/=2) then
     1244      if(j1.ne.2) then
    12951245      do i=1,IMR
    12961246      DC(i,1) = 0.
    … …  
    13001250C Determine slopes in polar caps for scalars!
    13011251C
    1302       do i=1,IMH
     1252      do 13 i=1,IMH
    13031253C South
    13041254      tmp = 0.25*(p(i,2) - p(i+imh,2))
    … …  
    13111261      Pmin = p(i,JNP) - min(p(i+imh,JMR),p(i,jnp), p(i,JMR))
    13121262      DC(i,JNP) = sign(min(abs(tmp),Pmax,pmin),tmp)
    1313       END DO
    1314 C
    1315       do i=imh+1,IMR
     126313    continue
     1264C
     1265      do 25 i=imh+1,IMR
    13161266      DC(i,  1) =  - DC(i-imh,  1)
    13171267      DC(i,JNP) =  - DC(i-imh,JNP)
    1318       END DO
     126825    continue
    13191269      endif
    13201270      return
    … …  
    13581308      LMT = JORD - 3     
    13591309C
    1360       DO i=1,IMR*JMR
     1310      DO 10 i=1,IMR*JMR       
    13611311      AL(i,2) = 0.5*(p(i,1)+p(i,2)) + (DC(i,1) - DC(i,2))*R3
    13621312      AR(i,1) = AL(i,2)
    1363       END DO
     131310    CONTINUE
    13641314C
    13651315CPoles:
    … …  
    13811331     
    13821332           
    1383       do i=1,len
    1384       A6(i,j11) = 3.*(p(i,j11)+p(i,j11)  - (AL(i,j11)+AR(i,j11)))
    1385       END DO
    1386 C
    1387       if(LMT<=2) call lmtppm(DC(1,j11),A6(1,j11),AR(1,j11)
     1333      do 30 i=1,len
     133430    A6(i,j11) = 3.*(p(i,j11)+p(i,j11)  - (AL(i,j11)+AR(i,j11)))
     1335C
     1336      if(LMT.le.2) call lmtppm(DC(1,j11),A6(1,j11),AR(1,j11)
    13881337     &                       ,AL(1,j11),P(1,j11),len,LMT)
    13891338C
    13901339     
    1391       DO i=1,IMJM1
    1392       IF(VC(i,j1)>0.) then
     1340      DO 140 i=1,IMJM1
     1341      IF(VC(i,j1).GT.0.) then
    13931342      flux(i,j1) = AR(i,j11) + 0.5*VC(i,j1)*(AL(i,j11) - AR(i,j11) +
    13941343     &                         A6(i,j11)*(1.-R23*VC(i,j1)) )
    … …  
    13971346     &                        A6(i,j1)*(1.+R23*VC(i,j1)))
    13981347      endif
    1399       END DO
     1348140   continue
    14001349      return
    14011350      end
    … …  
    14291378c        write(*,*) 'toto 1'
    14301379C --------------------------------
    1431       IF(IAD==2) then
     1380      IF(IAD.eq.2) then
    14321381      do j=j1-1,j2+1
    14331382      do i=1,IMR
    … …  
    14461395c      write(*,*) 'toto 2'
    14471396C
    1448       ELSEIF(IAD==1) then
     1397      ELSEIF(IAD.eq.1) then
    14491398        do j=j1-1,j2+1
    14501399      do i=1,imr
    … …  
    14551404      ENDIF
    14561405C
    1457         if(j1/=2) then
     1406        if(j1.ne.2) then
    14581407        sum1 = 0.
    14591408        sum2 = 0.
    … …  
    14991448C
    15001449        JMR = JNP-1
    1501       do j=j1,j2
    1502       if(J>JS  .and. J<JN) GO TO 1309
     1450      do 1309 j=j1,j2
     1451      if(J.GT.JS  .and. J.LT.JN) GO TO 1309
    15031452C
    15041453      do i=1,IMR
    … …  
    15111460      enddo
    15121461C
    1513       IF(IAD==2) THEN
     1462      IF(IAD.eq.2) THEN
    15141463      DO i=1,IMR
    15151464      IP = NINT(UA(i,j))
    … …  
    15201469      adx(i,j) = qtmp(ip) + ru*(a1*ru + b1)
    15211470      enddo
    1522       ELSEIF(IAD==1) then
     1471      ELSEIF(IAD.eq.1) then
    15231472      DO i=1,IMR
    15241473      iu = UA(i,j)
    15251474      ru = UA(i,j) - iu
    15261475      iiu = i-iu
    1527       if(UA(i,j)>=0.) then
     1476      if(UA(i,j).GE.0.) then
    15281477      adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))
    15291478      else
    … …  
    15371486      enddo
    153814871309  continue
    1539       END DO
    15401488C
    15411489C Eulerian upwind
    … …  
    15501498      qtmp(IMR+1) = p(1,J)
    15511499C
    1552       IF(IAD==2) THEN
     1500      IF(IAD.eq.2) THEN
    15531501      qtmp(-1)     = p(IMR-1,J)
    15541502      qtmp(IMR+2) = p(2,J)
    … …  
    15611509      adx(i,j) = qtmp(ip)- p(i,j) + ru*(a1*ru + b1)
    15621510      enddo
    1563       ELSEIF(IAD==1) then
     1511      ELSEIF(IAD.eq.1) then
    15641512C 1st order
    15651513      DO i=1,IMR
    … …  
    15701518      enddo
    15711519C
    1572         if(j1/=2) then
     1520        if(j1.ne.2) then
    15731521      do i=1,IMR
    15741522      adx(i,  2) = 0.
    … …  
    16061554      REAL da1,da2,a6da,fmin
    16071555C
    1608       if(LMT==0) then
     1556      if(LMT.eq.0) then
    16091557C Full constraint
    1610       do i=1,IM
    1611       if(DC(i)==0.) then
     1558      do 100 i=1,IM
     1559      if(DC(i).eq.0.) then
    16121560            AR(i) = p(i)
    16131561            AL(i) = p(i)
    … …  
    16171565      da2  = da1**2
    16181566      A6DA = A6(i)*da1
    1619       if(A6DA < -da2) then
     1567      if(A6DA .lt. -da2) then
    16201568            A6(i) = 3.*(AL(i)-p(i))
    16211569            AR(i) = AL(i) - A6(i)
    1622       elseif(A6DA > da2) then
     1570      elseif(A6DA .gt. da2) then
    16231571            A6(i) = 3.*(AR(i)-p(i))
    16241572            AL(i) = AR(i) - A6(i)
    16251573      endif
    16261574      endif
    1627       END DO
    1628       elseif(LMT==1) then
     1575100   continue
     1576      elseif(LMT.eq.1) then
    16291577C Semi-monotonic constraint
    1630       do i=1,IM
    1631       if(abs(AR(i)-AL(i)) >= -A6(i)) go to 150
    1632       if(p(i)<AR(i) .and. p(i)<AL(i)) then
     1578      do 150 i=1,IM
     1579      if(abs(AR(i)-AL(i)) .GE. -A6(i)) go to 150
     1580      if(p(i).lt.AR(i) .and. p(i).lt.AL(i)) then
    16331581            AR(i) = p(i)
    16341582            AL(i) = p(i)
    16351583            A6(i) = 0.
    1636       elseif(AR(i) > AL(i)) then
     1584      elseif(AR(i) .gt. AL(i)) then
    16371585            A6(i) = 3.*(AL(i)-p(i))
    16381586            AR(i) = AL(i) - A6(i)
    … …  
    16421590      endif
    16431591150   continue
    1644       END DO
    1645       elseif(LMT==2) then
    1646       do i=1,IM
    1647       if(abs(AR(i)-AL(i)) >= -A6(i)) go to 250
     1592      elseif(LMT.eq.2) then
     1593      do 250 i=1,IM
     1594      if(abs(AR(i)-AL(i)) .GE. -A6(i)) go to 250
    16481595      fmin = p(i) + 0.25*(AR(i)-AL(i))**2/A6(i) + A6(i)*R12
    1649       if(fmin>=0.) go to 250
    1650       if(p(i)<AR(i) .and. p(i)<AL(i)) then
     1596      if(fmin.ge.0.) go to 250
     1597      if(p(i).lt.AR(i) .and. p(i).lt.AL(i)) then
    16511598            AR(i) = p(i)
    16521599            AL(i) = p(i)
    16531600            A6(i) = 0.
    1654       elseif(AR(i) > AL(i)) then
     1601      elseif(AR(i) .gt. AL(i)) then
    16551602            A6(i) = 3.*(AL(i)-p(i))
    16561603            AR(i) = AL(i) - A6(i)
    … …  
    16601607      endif
    16611608250   continue
    1662       END DO
    16631609      endif
    16641610      return
    … …  
    16711617      integer i,j
    16721618C
    1673       do j=j1,j2
    1674       do i=2,IMR
    1675       CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j))
    1676       END DO
    1677       END DO
    1678 C
    1679       do j=j1,j2
    1680       CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j))
    1681       END DO
    1682 C
    1683       do i=1,IMR*JMR
    1684       CRY(i,2) = DTDY5*(V(i,2)+V(i,1))
    1685       END DO
     1619      do 35 j=j1,j2
     1620      do 35 i=2,IMR
     162135    CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j))
     1622C
     1623      do 45 j=j1,j2
     162445    CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j))
     1625C
     1626      do 55 i=1,IMR*JMR
     162755    CRY(i,2) = DTDY5*(V(i,2)+V(i,1))
    16861628      return
    16871629      end
    … …  
    16941636      real ph5
    16951637      JMR = JNP-1
    1696       do j=2,JNP
     1638      do 55 j=2,JNP
    16971639        ph5  =  -0.5*PI + (REAL(J-1)-0.5)*DP
    1698         cose(j) = cos(ph5)
    1699       END DO
     164055      cose(j) = cos(ph5)
    17001641C
    17011642      JEQ = (JNP+1) / 2
    1702       if(JMR == 2*(JMR/2) ) then
     1643      if(JMR .eq. 2*(JMR/2) ) then
    17031644      do j=JNP, JEQ+1, -1
    17041645       cose(j) =  cose(JNP+2-j)
    … …  
    17121653      endif
    17131654C
    1714       do j=2,JMR
    1715       cosp(j) = 0.5*(cose(j)+cose(j+1))
    1716       END DO
     1655      do 66 j=2,JMR
     165666    cosp(j) = 0.5*(cose(j)+cose(j+1))
    17171657      cosp(1) = 0.
    17181658      cosp(JNP) = 0.
    … …  
    17281668C
    17291669      phi = -0.5*PI
    1730       do j=2,JNP-1
     1670      do 55 j=2,JNP-1
    17311671      phi  =  phi + DP
    1732       cosp(j) = cos(phi)
    1733       END DO
     167255    cosp(j) = cos(phi)
    17341673        cosp(  1) = 0.
    17351674        cosp(JNP) = 0.
    17361675C
    1737       do j=2,JNP
     1676      do 66 j=2,JNP
    17381677        cose(j) = 0.5*(cosp(j)+cosp(j-1))
    1739       END DO
    1740 C
    1741       do j=2,JNP-1
     167866    CONTINUE
     1679C
     1680      do 77 j=2,JNP-1
    17421681       cosp(j) = 0.5*(cose(j)+cose(j+1))
    1743       END DO
     168277    CONTINUE
    17441683      return
    17451684      end
    … …  
    17631702        icr = 1
    17641703      call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
    1765       if(ipy==0) goto 50
     1704      if(ipy.eq.0) goto 50
    17661705      call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
    1767       if(ipx==0) goto 50
     1706      if(ipx.eq.0) goto 50
    17681707C
    17691708      if(cross) then
    17701709      call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
    17711710      endif
    1772       if(icr==0) goto 50
     1711      if(icr.eq.0) goto 50
    17731712C
    17741713C Vertical filling...
    17751714      do i=1,len
    1776       IF( Q(i,j1,1)<0.) THEN
     1715      IF( Q(i,j1,1).LT.0.) THEN
    17771716      ip = ip + 1
    17781717          Q(i,j1,2) = Q(i,j1,2) + Q(i,j1,1)
    … …  
    17821721C
    1783172250    continue
    1784       DO L = 2,NLAYM1
     1723      DO 225 L = 2,NLAYM1
    17851724      icr = 1
    17861725C
    17871726      call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
    1788       if(ipy==0) goto 225
     1727      if(ipy.eq.0) goto 225
    17891728      call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
    1790       if(ipx==0) go to 225
     1729      if(ipx.eq.0) go to 225
    17911730      if(cross) then
    17921731      call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
    17931732      endif
    1794       if(icr==0) goto 225
     1733      if(icr.eq.0) goto 225
    17951734C
    17961735      do i=1,len
    1797       IF( Q(I,j1,L)<0.) THEN
     1736      IF( Q(I,j1,L).LT.0.) THEN
    17981737C
    17991738      ip = ip + 1
    … …  
    18101749      ENDDO
    18111750225   CONTINUE
    1812       END DO
    18131751C
    18141752C BOTTOM LAYER
    … …  
    18171755C
    18181756      call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
    1819       if(ipy==0) goto 911
     1757      if(ipy.eq.0) goto 911
    18201758      call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
    1821       if(ipx==0) goto 911
     1759      if(ipx.eq.0) goto 911
    18221760C
    18231761      call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
    1824       if(icr==0) goto 911
     1762      if(icr.eq.0) goto 911
    18251763C
    18261764      DO  I=1,len
    1827       IF( Q(I,j1,L)<0.) THEN
     1765      IF( Q(I,j1,L).LT.0.) THEN
    18281766      ip = ip + 1
    18291767c
    … …  
    18421780911   continue
    18431781C
    1844       if(ip>IMR) then
     1782      if(ip.gt.IMR) then
    18451783      write(6,*) 'IC=',IC,' STEP=',NSTEP,
    18461784     &           ' Vertical filling pts=',ip
    18471785      endif
    18481786C
    1849       if(sum>1.e-25) then
     1787      if(sum.gt.1.e-25) then
    18501788      write(6,*) IC,NSTEP,' Mass source from the ground=',sum
    18511789      endif
    … …  
    18601798      real :: dq,dn,d0,d1,ds,d2
    18611799      icr = 0
    1862       do j=j1+1,j2-1
    1863       DO i=1,IMR-1
    1864       IF(q(i,j)<0.) THEN
     1800      do 65 j=j1+1,j2-1
     1801      DO 50 i=1,IMR-1
     1802      IF(q(i,j).LT.0.) THEN
    18651803      icr =  1
    18661804      dq  = - q(i,j)*cosp(j)
    … …  
    18781816      q(i,j) = (d2 - dq)*acosp(j) + tiny
    18791817      endif
    1880       END DO
    1881       if(icr==0 .and. q(IMR,j)>=0.) goto 65
    1882       DO i=2,IMR
    1883       IF(q(i,j)<0.) THEN
     181850    continue
     1819      if(icr.eq.0 .and. q(IMR,j).ge.0.) goto 65
     1820      DO 55 i=2,IMR
     1821      IF(q(i,j).LT.0.) THEN
    18841822      icr =  1
    18851823      dq  = - q(i,j)*cosp(j)
    … …  
    18971835      q(i,j) = (d2 - dq)*acosp(j) + tiny
    18981836      endif
    1899       END DO
     183755    continue
    19001838C *****************************************
    19011839C i=1
    19021840      i=1
    1903       IF(q(i,j)<0.) THEN
     1841      IF(q(i,j).LT.0.) THEN
    19041842      icr =  1
    19051843      dq  = - q(i,j)*cosp(j)
    … …  
    19201858C i=IMR
    19211859      i=IMR
    1922       IF(q(i,j)<0.) THEN
     1860      IF(q(i,j).LT.0.) THEN
    19231861      icr =  1
    19241862      dq  = - q(i,j)*cosp(j)
    … …  
    19381876C *****************************************
    1939187765    continue
    1940       END DO
    19411878C
    19421879      do i=1,IMR
    1943       if(q(i,j1)<0. .or. q(i,j2)<0.) then
     1880      if(q(i,j1).lt.0. .or. q(i,j2).lt.0.) then
    19441881      icr = 1
    19451882      goto 80
    … …  
    1949188680    continue
    19501887C
    1951       if(q(1,1)<0. .or. q(1,jnp)<0.) then
     1888      if(q(1,1).lt.0. .or. q(1,jnp).lt.0.) then
    19521889      icr = 1
    19531890      endif
    … …  
    19731910C
    19741911      ipy = 0
    1975       do j=j1+1,j2-1
    1976       DO i=1,IMR
    1977       IF(q(i,j)<0.) THEN
     1912      do 55 j=j1+1,j2-1
     1913      DO 55 i=1,IMR
     1914      IF(q(i,j).LT.0.) THEN
    19781915      ipy =  1
    19791916      dq  = - q(i,j)*cosp(j)
    … …  
    19911928      q(i,j) = (d2 - dq)*acosp(j) + tiny
    19921929      endif
    1993       END DO
    1994       END DO
     193055    continue
    19951931C
    19961932      do i=1,imr
    1997       IF(q(i,j1)<0.) THEN
     1933      IF(q(i,j1).LT.0.) THEN
    19981934      ipy =  1
    19991935      dq  = - q(i,j1)*cosp(j1)
    … …  
    20091945      j = j2
    20101946      do i=1,imr
    2011       IF(q(i,j)<0.) THEN
     1947      IF(q(i,j).LT.0.) THEN
    20121948      ipy =  1
    20131949      dq  = - q(i,j)*cosp(j)
    … …  
    20221958C
    20231959C Check Poles.
    2024       if(q(1,1)<0.) then
     1960      if(q(1,1).lt.0.) then
    20251961      dq = q(1,1)*cap1/REAL(IMR)*acosp(j1)
    20261962      do i=1,imr
    20271963      q(i,1) = 0.
    20281964      q(i,j1) = q(i,j1) + dq
    2029       if(q(i,j1)<0.) ipy = 1
    2030       enddo
    2031       endif
    2032 C
    2033       if(q(1,JNP)<0.) then
     1965      if(q(i,j1).lt.0.) ipy = 1
     1966      enddo
     1967      endif
     1968C
     1969      if(q(1,JNP).lt.0.) then
    20341970      dq = q(1,JNP)*cap1/REAL(IMR)*acosp(j2)
    20351971      do i=1,imr
    20361972      q(i,JNP) = 0.
    20371973      q(i,j2) = q(i,j2) + dq
    2038       if(q(i,j2)<0.) ipy = 1
     1974      if(q(i,j2).lt.0.) ipy = 1
    20391975      enddo
    20401976      endif
    … …  
    20521988      ipx = 0
    20531989C Copy & swap direction for vectorization.
    2054       do i=1,imr
    2055       do j=j1,j2
    2056       qtmp(j,i) = q(i,j)
    2057       END DO
    2058       END DO
    2059 C
    2060       do i=2,imr-1
    2061       do j=j1,j2
    2062       if(qtmp(j,i)<0.) then
     1990      do 25 i=1,imr
     1991      do 25 j=j1,j2
     199225    qtmp(j,i) = q(i,j)
     1993C
     1994      do 55 i=2,imr-1
     1995      do 55 j=j1,j2
     1996      if(qtmp(j,i).lt.0.) then
    20631997      ipx =  1
    20641998c west
    … …  
    20732007      qtmp(j,i) = qtmp(j,i) + d2 + tiny
    20742008      endif
    2075       END DO
    2076       END DO
     200955    continue
    20772010c
    20782011      i=1
    2079       do j=j1,j2
    2080       if(qtmp(j,i)<0.) then
     2012      do 65 j=j1,j2
     2013      if(qtmp(j,i).lt.0.) then
    20812014      ipx =  1
    20822015c west
    … …  
    20922025      qtmp(j,i) = qtmp(j,i) + d2 + tiny
    20932026      endif
    2094       END DO
     202765    continue
    20952028      i=IMR
    2096       do j=j1,j2
    2097       if(qtmp(j,i)<0.) then
     2029      do 75 j=j1,j2
     2030      if(qtmp(j,i).lt.0.) then
    20982031      ipx =  1
    20992032c west
    … …  
    21092042      qtmp(j,i) = qtmp(j,i) + d2 + tiny
    21102043      endif
    2111       END DO
    2112 C
    2113       if(ipx/=0) then
    2114       do j=j1,j2
    2115       do i=1,imr
    2116       q(i,j) = qtmp(j,i)
    2117       END DO
    2118       END DO
     204475    continue
     2045C
     2046      if(ipx.ne.0) then
     2047      do 85 j=j1,j2
     2048      do 85 i=1,imr
     204985    q(i,j) = qtmp(j,i)
    21192050      else
    21202051C
    21212052C Poles.
    2122       if(q(1,1)<0 .or. q(1,JNP)<0.) ipx = 1
     2053      if(q(1,1).lt.0. or. q(1,JNP).lt.0.) ipx = 1
    21232054      endif
    21242055      return
    … …  
    21342065      integer IC,k,i
    21352066C
    2136       do IC = 1, nc
    2137 C
    2138       do k=1,km
    2139       do i=1,im
     2067      do 4000 IC = 1, nc
     2068C
     2069      do 1000 k=1,km
     2070      do 1000 i=1,im
    21402071      qtmp(i,k) = q(i,km+1-k,IC)
    2141       END DO
    2142       END DO
    2143 C
    2144       do i=1,im*km
    2145       q(i,1,IC) = qtmp(i,1)
    2146       END DO
    2147       END DO
     20721000  continue
     2073C
     2074      do 2000 i=1,im*km
     20752000  q(i,1,IC) = qtmp(i,1)
     20764000  continue
    21482077      return
    21492078      end
Note: See TracChangeset for help on using the changeset viewer.