Changeset 5079 for LMDZ6/trunk


Ignore:
Timestamp:
Jul 19, 2024, 11:28:59 AM (4 months ago)
Author:
abarral
Message:

[coherence with Fortran standards]
Replace obsolete DO with shared termination
(minor) replace obsolete bool operators

Location:
LMDZ6/trunk/libf
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r2600 r5079  
    121121      enddo
    122122
    123       DO 1 L=1,NIV
     123      DO L=1,NIV
    124124C
    125125C  place limits on appropriate moments before transport
    … …  
    128128      IF(.NOT.LIMIT) GO TO 11
    129129C
    130       DO 10 JV=1,NTRA
    131       DO 10 K=1,LAT
    132       DO 100 I=1,LON
     130      DO JV=1,NTRA
     131      DO K=1,LAT
     132      DO 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  100  CONTINUE
    136  10   CONTINUE
     135      END DO
     136      END DO
     137      END DO
    137138C
    138139 11   CONTINUE
    … …  
    141142C
    142143      SM0=0.
    143       DO 20 JV=1,NTRA
     144      DO JV=1,NTRA
    144145         S00(JV)=0.
    145  20   CONTINUE
    146 C
    147       DO 21 I=1,LON
    148 C
    149          IF(VGRI(I,0,L).LE.0.) THEN
     146      END DO
     147C
     148      DO I=1,LON
     149C
     150         IF(VGRI(I,0,L)<=0.) THEN
    150151           FM(I,0)=-VGRI(I,0,L)*DTY
    151152           ALF(I,0)=FM(I,0)/SM(I,1,L)
    … …  
    158159         ALF1Q(I,0)=ALF1(I,0)*ALF1(I,0)
    159160C
    160  21   CONTINUE
    161 C
    162       DO 22 JV=1,NTRA
    163       DO 220 I=1,LON
    164 C
    165          IF(VGRI(I,0,L).LE.0.) THEN
     161      END DO
     162C
     163      DO JV=1,NTRA
     164      DO I=1,LON
     165C
     166         IF(VGRI(I,0,L)<=0.) THEN
    166167C
    167168           F0(I,0,JV)=ALF(I,0)*
    … …  
    176177         ENDIF
    177178C
    178  220  CONTINUE
    179  22   CONTINUE
    180 C
    181       DO 23 I=1,LON
    182          IF(VGRI(I,0,L).GT.0.) THEN
     179      END DO
     180      END DO
     181C
     182      DO I=1,LON
     183         IF(VGRI(I,0,L)>0.) THEN
    183184           FM(I,0)=VGRI(I,0,L)*DTY
    184185           ALF(I,0)=FM(I,0)/SM0
    185186         ENDIF
    186  23   CONTINUE
    187 C
    188       DO 24 JV=1,NTRA
    189       DO 240 I=1,LON
    190          IF(VGRI(I,0,L).GT.0.) THEN
     187      END DO
     188C
     189      DO JV=1,NTRA
     190      DO I=1,LON
     191         IF(VGRI(I,0,L)>0.) THEN
    191192           F0(I,0,JV)=ALF(I,0)*S00(JV)
    192193         ENDIF
    193  240  CONTINUE
    194  24   CONTINUE
     194      END DO
     195      END DO
    195196C
    196197C  puts the temporary moments Fi into appropriate neighboring boxes
    197198C
    198       DO 25 I=1,LON
    199 C
    200          IF(VGRI(I,0,L).GT.0.) THEN
     199      DO I=1,LON
     200C
     201         IF(VGRI(I,0,L)>0.) THEN
    201202           SM(I,1,L)=SM(I,1,L)+FM(I,0)
    202203           ALF(I,0)=FM(I,0)/SM(I,1,L)
    … …  
    205206         ALF1(I,0)=1.-ALF(I,0)
    206207C
    207  25   CONTINUE
    208 C
    209       DO 26 JV=1,NTRA
    210       DO 260 I=1,LON
    211 C
    212          IF(VGRI(I,0,L).GT.0.) THEN
     208      END DO
     209C
     210      DO JV=1,NTRA
     211      DO I=1,LON
     212C
     213         IF(VGRI(I,0,L)>0.) THEN
    213214C
    214215         TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV)
    … …  
    218219         ENDIF
    219220C
    220  260  CONTINUE
    221  26   CONTINUE
     221      END DO
     222      END DO
    222223C
    223224C  calculate flux and moments between adjacent boxes
    … …  
    227228C  flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
    228229C
    229       DO 30 K=1,LAT-1
     230      DO K=1,LAT-1
    230231      KP=K+1
    231       DO 300 I=1,LON
    232 C
    233          IF(VGRI(I,K,L).LT.0.) THEN
     232      DO I=1,LON
     233C
     234         IF(VGRI(I,K,L)<0.) THEN
    234235           FM(I,K)=-VGRI(I,K,L)*DTY
    235236           ALF(I,K)=FM(I,K)/SM(I,KP,L)
    … …  
    245246         ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
    246247C
    247  300  CONTINUE
    248  30   CONTINUE
    249 C
    250       DO 31 JV=1,NTRA
    251       DO 31 K=1,LAT-1
     248      END DO
     249      END DO
     250C
     251      DO JV=1,NTRA
     252      DO K=1,LAT-1
    252253      KP=K+1
    253       DO 310 I=1,LON
    254 C
    255          IF(VGRI(I,K,L).LT.0.) THEN
     254      DO I=1,LON
     255C
     256         IF(VGRI(I,K,L)<0.) THEN
    256257C
    257258           F0(I,K,JV)=ALF (I,K)*
    … …  
    281282         ENDIF
    282283C
    283  310  CONTINUE
    284  31   CONTINUE
     284      END DO
     285      END DO
     286      END DO
    285287C
    286288C  puts the temporary moments Fi into appropriate neighboring boxes
    287289C
    288       DO 32 K=1,LAT-1
     290      DO K=1,LAT-1
    289291      KP=K+1
    290       DO 320 I=1,LON
    291 C
    292          IF(VGRI(I,K,L).LT.0.) THEN
     292      DO I=1,LON
     293C
     294         IF(VGRI(I,K,L)<0.) THEN
    293295           SM(I,K,L)=SM(I,K,L)+FM(I,K)
    294296           ALF(I,K)=FM(I,K)/SM(I,K,L)
    … …  
    300302         ALF1(I,K)=1.-ALF(I,K)
    301303C
    302  320  CONTINUE
    303  32   CONTINUE
    304 C
    305       DO 33 JV=1,NTRA
    306       DO 33 K=1,LAT-1
     304      END DO
     305      END DO
     306C
     307      DO JV=1,NTRA
     308      DO K=1,LAT-1
    307309      KP=K+1
    308       DO 330 I=1,LON
    309 C
    310          IF(VGRI(I,K,L).LT.0.) THEN
     310      DO I=1,LON
     311C
     312         IF(VGRI(I,K,L)<0.) THEN
    311313C
    312314         TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
    … …  
    328330         ENDIF
    329331C
    330  330  CONTINUE
    331  33   CONTINUE
     332      END DO
     333      END DO
     334      END DO
    332335C
    333336C  traitement special pour le pole Sud (idem pole Nord)
    … …  
    336339C
    337340      SM0=0.
    338       DO 40 JV=1,NTRA
     341      DO JV=1,NTRA
    339342         S00(JV)=0.
    340  40   CONTINUE
    341 C
    342       DO 41 I=1,LON
    343 C
    344          IF(VGRI(I,K,L).GE.0.) THEN
     343      END DO
     344C
     345      DO I=1,LON
     346C
     347         IF(VGRI(I,K,L)>=0.) THEN
    345348           FM(I,K)=VGRI(I,K,L)*DTY
    346349           ALF(I,K)=FM(I,K)/SM(I,K,L)
    … …  
    353356         ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
    354357C
    355  41   CONTINUE
    356 C
    357       DO 42 JV=1,NTRA
    358       DO 420 I=1,LON
    359 C
    360          IF(VGRI(I,K,L).GE.0.) THEN
     358      END DO
     359C
     360      DO JV=1,NTRA
     361      DO I=1,LON
     362C
     363         IF(VGRI(I,K,L)>=0.) THEN
    361364           F0 (I,K,JV)=ALF(I,K)*
    362365     +                ( S0(I,K,L,JV)+ALF1(I,K)*sy(I,K,L,JV) )
    … …  
    369372         ENDIF
    370373C
    371  420  CONTINUE
    372  42   CONTINUE
    373 C
    374       DO 43 I=1,LON
    375          IF(VGRI(I,K,L).LT.0.) THEN
     374      END DO
     375      END DO
     376C
     377      DO I=1,LON
     378         IF(VGRI(I,K,L)<0.) THEN
    376379           FM(I,K)=-VGRI(I,K,L)*DTY
    377380           ALF(I,K)=FM(I,K)/SM0
    378381         ENDIF
    379  43   CONTINUE
    380 C
    381       DO 44 JV=1,NTRA
    382       DO 440 I=1,LON
    383          IF(VGRI(I,K,L).LT.0.) THEN
     382      END DO
     383C
     384      DO JV=1,NTRA
     385      DO I=1,LON
     386         IF(VGRI(I,K,L)<0.) THEN
    384387           F0(I,K,JV)=ALF(I,K)*S00(JV)
    385388         ENDIF
    386  440  CONTINUE
    387  44   CONTINUE
     389      END DO
     390      END DO
    388391C
    389392C  puts the temporary moments Fi into appropriate neighboring boxes
    390393C
    391       DO 45 I=1,LON
    392 C
    393          IF(VGRI(I,K,L).LT.0.) THEN
     394      DO I=1,LON
     395C
     396         IF(VGRI(I,K,L)<0.) THEN
    394397           SM(I,K,L)=SM(I,K,L)+FM(I,K)
    395398           ALF(I,K)=FM(I,K)/SM(I,K,L)
    … …  
    398401         ALF1(I,K)=1.-ALF(I,K)
    399402C
    400  45   CONTINUE
    401 C
    402       DO 46 JV=1,NTRA
    403       DO 460 I=1,LON
    404 C
    405          IF(VGRI(I,K,L).LT.0.) THEN
     403      END DO
     404C
     405      DO JV=1,NTRA
     406      DO I=1,LON
     407C
     408         IF(VGRI(I,K,L)<0.) THEN
    406409C
    407410         TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
    … …  
    411414         ENDIF
    412415C
    413  460  CONTINUE
    414  46   CONTINUE
    415 C
    416  1    CONTINUE
     416      END DO
     417      END DO
     418C
     419      END DO
    417420C
    418421      RETURN

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

    r2600 r5079  
    153153C-AA 20/10/94  le signe -1 est necessaire car indexation opposee
    154154
    155       DO 500 l = 1,llm
    156          DO 500 j = 1,jjm
    157             DO 500 i = 1,iip1 
     155      DO l = 1,llm
     156         DO j = 1,jjm
     157            DO i = 1,iip1
    158158            vgri (i,j,llm+1-l)=-1.*pbarv (i,j,l)
    159   500 CONTINUE
     159      END DO
     160      END DO
     161      END DO
    160162
    161163CAA Initialisation de flux fictifs aux bords sup. des boites pol.
    … …  
    171173C  boucle sur les niveaux
    172174C
    173       DO 1 L=1,NIV
     175      DO L=1,NIV
    174176C
    175177C  place limits on appropriate moments before transport
    … …  
    178180      IF(.NOT.LIMIT) GO TO 11
    179181C
    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
     182      DO JV=1,NTRA
     183      DO K=1,LAT
     184      DO I=1,LON
     185         IF(S0(I,K,L,JV)>0.) THEN
    184186           SLPMAX=AMAX1(S0(I,K,L,JV),0.)
    185187           S1MAX=1.5*SLPMAX
    … …  
    197199           SYZ(I,K,L,JV)=0.
    198200         ENDIF
    199  100  CONTINUE
    200  10   CONTINUE
     201      END DO
     202      END DO
     203      END DO
    201204C
    202205 11   CONTINUE
    … …  
    205208C
    206209      SM0=0.
    207       DO 20 JV=1,NTRA
     210      DO JV=1,NTRA
    208211         S00(JV)=0.
    209  20   CONTINUE
    210 C
    211       DO 21 I=1,LON
    212 C
    213          IF(VGRI(I,0,L).LE.0.) THEN
     212      END DO
     213C
     214      DO I=1,LON
     215C
     216         IF(VGRI(I,0,L)<=0.) THEN
    214217           FM(I,0)=-VGRI(I,0,L)*DTY
    215218           ALF(I,0)=FM(I,0)/SM(I,1,L)
    … …  
    225228         ALF4(I,0)=ALF1(I,0)*ALF1Q(I,0)
    226229C
    227  21   CONTINUE
     230      END DO
    228231c     print*,'ADVYP 21'
    229232C
    230       DO 22 JV=1,NTRA
    231       DO 220 I=1,LON
    232 C
    233          IF(VGRI(I,0,L).LE.0.) THEN
     233      DO JV=1,NTRA
     234      DO I=1,LON
     235C
     236         IF(VGRI(I,0,L)<=0.) THEN
    234237C
    235238           F0(I,0,JV)=ALF(I,0)* ( S0(I,1,L,JV)-ALF1(I,0)*
    … …  
    253256         ENDIF
    254257C
    255  220  CONTINUE
    256  22   CONTINUE
    257 C
    258       DO 23 I=1,LON
    259          IF(VGRI(I,0,L).GT.0.) THEN
     258      END DO
     259      END DO
     260C
     261      DO I=1,LON
     262         IF(VGRI(I,0,L)>0.) THEN
    260263           FM(I,0)=VGRI(I,0,L)*DTY
    261264           ALF(I,0)=FM(I,0)/SM0
    262265         ENDIF
    263  23   CONTINUE
    264 C
    265       DO 24 JV=1,NTRA
    266       DO 240 I=1,LON
    267          IF(VGRI(I,0,L).GT.0.) THEN
     266      END DO
     267C
     268      DO JV=1,NTRA
     269      DO I=1,LON
     270         IF(VGRI(I,0,L)>0.) THEN
    268271           F0(I,0,JV)=ALF(I,0)*S00(JV)
    269272         ENDIF
    270  240  CONTINUE
    271  24   CONTINUE
     273      END DO
     274      END DO
    272275C
    273276C  puts the temporary moments Fi into appropriate neighboring boxes
    274277C
    275278c     print*,'av ADVYP 25'
    276       DO 25 I=1,LON
    277 C
    278          IF(VGRI(I,0,L).GT.0.) THEN
     279      DO I=1,LON
     280C
     281         IF(VGRI(I,0,L)>0.) THEN
    279282           SM(I,1,L)=SM(I,1,L)+FM(I,0)
    280283           ALF(I,0)=FM(I,0)/SM(I,1,L)
    … …  
    287290         ALF3(I,0)=ALF1(I,0)*ALF(I,0)
    288291C
    289  25   CONTINUE
     292      END DO
    290293c     print*,'av ADVYP 25'
    291294C
    292       DO 26 JV=1,NTRA
    293       DO 260 I=1,LON
    294 C
    295          IF(VGRI(I,0,L).GT.0.) THEN
     295      DO JV=1,NTRA
     296      DO I=1,LON
     297C
     298         IF(VGRI(I,0,L)>0.) THEN
    296299C
    297300         TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV)
    … …  
    305308         ENDIF
    306309C
    307  260  CONTINUE
    308  26   CONTINUE
     310      END DO
     311      END DO
    309312C
    310313C  calculate flux and moments between adjacent boxes
    … …  
    315318C
    316319c     print*,'av ADVYP 30'
    317       DO 30 K=1,LAT-1
     320      DO K=1,LAT-1
    318321      KP=K+1
    319       DO 300 I=1,LON
    320 C
    321          IF(VGRI(I,K,L).LT.0.) THEN
     322      DO I=1,LON
     323C
     324         IF(VGRI(I,K,L)<0.) THEN
    322325           FM(I,K)=-VGRI(I,K,L)*DTY
    323326           ALF(I,K)=FM(I,K)/SM(I,KP,L)
    … …  
    336339         ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K)
    337340C
    338  300  CONTINUE
    339  30   CONTINUE
     341      END DO
     342      END DO
    340343c     print*,'ap ADVYP 30'
    341344C
    342       DO 31 JV=1,NTRA
    343       DO 31 K=1,LAT-1
     345      DO JV=1,NTRA
     346      DO K=1,LAT-1
    344347      KP=K+1
    345       DO 310 I=1,LON
    346 C
    347          IF(VGRI(I,K,L).LT.0.) THEN
     348      DO I=1,LON
     349C
     350         IF(VGRI(I,K,L)<0.) THEN
    348351C
    349352           F0 (I,K,JV)=ALF (I,K)* ( S0(I,KP,L,JV)-ALF1(I,K)*
    … …  
    403406         ENDIF
    404407C
    405  310  CONTINUE
    406  31   CONTINUE
     408      END DO
     409      END DO
     410      END DO
    407411c     print*,'ap ADVYP 31'
    408412C
    409413C  puts the temporary moments Fi into appropriate neighboring boxes
    410414C
    411       DO 32 K=1,LAT-1
     415      DO K=1,LAT-1
    412416      KP=K+1
    413       DO 320 I=1,LON
    414 C
    415          IF(VGRI(I,K,L).LT.0.) THEN
     417      DO I=1,LON
     418C
     419         IF(VGRI(I,K,L)<0.) THEN
    416420           SM(I,K,L)=SM(I,K,L)+FM(I,K)
    417421           ALF(I,K)=FM(I,K)/SM(I,K,L)
    … …  
    427431         ALF3(I,K)=ALF1(I,K)*ALF(I,K)
    428432C
    429  320  CONTINUE
    430  32   CONTINUE
     433      END DO
     434      END DO
    431435c     print*,'ap ADVYP 32'
    432436C
    433       DO 33 JV=1,NTRA
    434       DO 33 K=1,LAT-1
     437      DO JV=1,NTRA
     438      DO K=1,LAT-1
    435439      KP=K+1
    436       DO 330 I=1,LON
    437 C
    438          IF(VGRI(I,K,L).LT.0.) THEN
     440      DO I=1,LON
     441C
     442         IF(VGRI(I,K,L)<0.) THEN
    439443C
    440444         TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
    … …  
    474478         ENDIF
    475479C
    476  330  CONTINUE
    477  33   CONTINUE
     480      END DO
     481      END DO
     482      END DO
    478483c     print*,'ap ADVYP 33'
    479484C
    … …  
    483488C
    484489      SM0=0.
    485       DO 40 JV=1,NTRA
     490      DO JV=1,NTRA
    486491         S00(JV)=0.
    487  40   CONTINUE
    488 C
    489       DO 41 I=1,LON
    490 C
    491          IF(VGRI(I,K,L).GE.0.) THEN
     492      END DO
     493C
     494      DO I=1,LON
     495C
     496         IF(VGRI(I,K,L)>=0.) THEN
    492497           FM(I,K)=VGRI(I,K,L)*DTY
    493498           ALF(I,K)=FM(I,K)/SM(I,K,L)
    … …  
    503508         ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K)
    504509C
    505  41   CONTINUE
     510      END DO
    506511c     print*,'ap ADVYP 41'
    507512C
    508       DO 42 JV=1,NTRA
    509       DO 420 I=1,LON
    510 C
    511          IF(VGRI(I,K,L).GE.0.) THEN
     513      DO JV=1,NTRA
     514      DO I=1,LON
     515C
     516         IF(VGRI(I,K,L)>=0.) THEN
    512517           F0 (I,K,JV)=ALF(I,K)* ( S0(I,K,L,JV)+ALF1(I,K)*
    513518     +             ( SY(I,K,L,JV)+ALF2(I,K)*SYY(I,K,L,JV) ) )
    … …  
    527532         ENDIF
    528533C
    529  420  CONTINUE
    530  42   CONTINUE
     534      END DO
     535      END DO
    531536c     print*,'ap ADVYP 42'
    532537C
    533       DO 43 I=1,LON
    534          IF(VGRI(I,K,L).LT.0.) THEN
     538      DO I=1,LON
     539         IF(VGRI(I,K,L)<0.) THEN
    535540           FM(I,K)=-VGRI(I,K,L)*DTY
    536541           ALF(I,K)=FM(I,K)/SM0
    537542         ENDIF
    538  43   CONTINUE
     543      END DO
    539544c     print*,'ap ADVYP 43'
    540545C
    541       DO 44 JV=1,NTRA
    542       DO 440 I=1,LON
    543          IF(VGRI(I,K,L).LT.0.) THEN
     546      DO JV=1,NTRA
     547      DO I=1,LON
     548         IF(VGRI(I,K,L)<0.) THEN
    544549           F0(I,K,JV)=ALF(I,K)*S00(JV)
    545550         ENDIF
    546  440  CONTINUE
    547  44   CONTINUE
     551      END DO
     552      END DO
    548553C
    549554C  puts the temporary moments Fi into appropriate neighboring boxes
    550555C
    551       DO 45 I=1,LON
    552 C
    553          IF(VGRI(I,K,L).LT.0.) THEN
     556      DO I=1,LON
     557C
     558         IF(VGRI(I,K,L)<0.) THEN
    554559           SM(I,K,L)=SM(I,K,L)+FM(I,K)
    555560           ALF(I,K)=FM(I,K)/SM(I,K,L)
    … …  
    562567         ALF3(I,K)=ALF1(I,K)*ALF(I,K)
    563568C
    564  45   CONTINUE
     569      END DO
    565570c     print*,'ap ADVYP 45'
    566571C
    567       DO 46 JV=1,NTRA
    568       DO 460 I=1,LON
    569 C
    570          IF(VGRI(I,K,L).LT.0.) THEN
     572      DO JV=1,NTRA
     573      DO I=1,LON
     574C
     575         IF(VGRI(I,K,L)<0.) THEN
    571576C
    572577         TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
    … …  
    580585         ENDIF
    581586C
    582  460  CONTINUE
    583  46   CONTINUE
     587      END DO
     588      END DO
    584589c     print*,'ap ADVYP 46'
    585590C
    586  1    CONTINUE
     591      END DO
    587592
    588593c--------------------------------------------------

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

    r4593 r5079  
    117117C  Conversion du flux de masse en kg.s-1
    118118
    119       DO 500 l = 1,llm
    120          DO 500 j = 1,jjp1
    121             DO 500 i = 1,iip1 
     119      DO l = 1,llm
     120         DO j = 1,jjp1
     121            DO 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   500 CONTINUE
     127      END DO
     128      END DO
     129      END DO
    128130         DO  j = 1,jjp1
    129131            DO i = 1,iip1 
    … …  
    137139C  boucle sur les latitudes
    138140C
    139       DO 1 K=1,LAT
     141      DO K=1,LAT
    140142C
    141143C  place limits on appropriate moments before transport
    … …  
    144146      IF(.NOT.LIMIT) GO TO 101
    145147C
    146       DO 10 JV=1,NTRA
    147       DO 10 L=1,NIV
    148          DO 100 I=1,LON
     148      DO JV=1,NTRA
     149      DO L=1,NIV
     150         DO I=1,LON
    149151            sz(I,K,L,JV)=SIGN(AMIN1(AMAX1(S0(I,K,L,JV),0.),
    150152     +                              ABS(sz(I,K,L,JV))),sz(I,K,L,JV))
    151  100     CONTINUE
    152  10   CONTINUE
     153      END DO
     154      END DO
     155      END DO
    153156C
    154157 101  CONTINUE
    … …  
    162165C  2- reajusts moments remaining in the box
    163166C
    164       DO 11 L=1,NIV-1
     167      DO L=1,NIV-1
    165168      LP=L+1
    166169C
    167       DO 110 I=1,LON
    168 C
    169          IF(WGRI(I,K,L).LT.0.) THEN
     170      DO I=1,LON
     171C
     172         IF(WGRI(I,K,L)<0.) THEN
    170173           FM(I,L)=-WGRI(I,K,L)*DTZ
    171174           ALF(I)=FM(I,L)/SM(I,K,LP)
    … …  
    181184         ALF1Q(I)=ALF1(I)*ALF1(I)
    182185C
    183  110  CONTINUE
    184 C
    185       DO 111 JV=1,NTRA
    186       DO 1110 I=1,LON
    187 C
    188          IF(WGRI(I,K,L).LT.0.) THEN
     186      END DO
     187C
     188      DO JV=1,NTRA
     189      DO I=1,LON
     190C
     191         IF(WGRI(I,K,L)<0.) THEN
    189192C
    190193           F0(I,L,JV)=ALF (I)*( S0(I,K,LP,JV)-ALF1(I)*sz(I,K,LP,JV) )
    … …  
    212215         ENDIF
    213216C
    214  1110 CONTINUE
    215  111  CONTINUE
    216 C
    217  11   CONTINUE
     217      END DO
     218      END DO
     219C
     220      END DO
    218221C
    219222C  puts the temporary moments Fi into appropriate neighboring boxes
    220223C
    221       DO 12 L=1,NIV-1
     224      DO L=1,NIV-1
    222225      LP=L+1
    223226C
    224       DO 120 I=1,LON
    225 C
    226          IF(WGRI(I,K,L).LT.0.) THEN
     227      DO I=1,LON
     228C
     229         IF(WGRI(I,K,L)<0.) THEN
    227230           SM(I,K,L)=SM(I,K,L)+FM(I,L)
    228231           ALF(I)=FM(I,L)/SM(I,K,L)
    … …  
    236239         ALF1Q(I)=ALF1(I)*ALF1(I)
    237240C
    238  120  CONTINUE
    239 C
    240       DO 121 JV=1,NTRA
    241       DO 1210 I=1,LON
    242 C
    243          IF(WGRI(I,K,L).LT.0.) THEN
     241      END DO
     242C
     243      DO JV=1,NTRA
     244      DO I=1,LON
     245C
     246         IF(WGRI(I,K,L)<0.) THEN
    244247C
    245248           TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)
    … …  
    260263         ENDIF
    261264C
    262  1210 CONTINUE
    263  121  CONTINUE
    264 C
    265  12   CONTINUE
     265      END DO
     266      END DO
     267C
     268      END DO
    266269C
    267270C  fin de la boucle principale sur les latitudes
    268271C
    269  1    CONTINUE
     272      END DO
    270273C
    271274C-------------------------------------------------------------

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

    r2600 r5079  
    135135C  Conversion des flux de masses en kg
    136136
    137       DO 500 l = 1,llm
    138          DO 500 j = 1,jjp1
    139             DO 500 i = 1,iip1 
     137      DO l = 1,llm
     138         DO j = 1,jjp1
     139            DO i = 1,iip1 
    140140            wgri (i,j,llm+1-l) = w (i,j,l) 
    141   500 CONTINUE
     141      END DO
     142      END DO
     143      END DO
    142144      do j=1,jjp1
    143145         do i=1,iip1
    … …  
    154156C  boucle sur les latitudes
    155157C
    156       DO 1 K=1,LAT
     158      DO K=1,LAT
    157159C
    158160C  place limits on appropriate moments before transport
    … …  
    161163      IF(.NOT.LIMIT) GO TO 101
    162164C
    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
     165      DO JV=1,NTRA
     166      DO L=1,NIV
     167         DO I=1,LON
     168            IF(S0(I,K,L,JV)>0.) THEN
    167169              SLPMAX=S0(I,K,L,JV)
    168170              S1MAX =1.5*SLPMAX
    … …  
    180182              SYZ(I,K,L,JV)=0.
    181183            ENDIF
    182  100     CONTINUE
    183  10   CONTINUE
     184      END DO
     185      END DO
     186      END DO
    184187C
    185188 101  CONTINUE
    … …  
    193196C  2- reajusts moments remaining in the box
    194197C
    195       DO 11 L=1,NIV-1
     198      DO L=1,NIV-1
    196199      LP=L+1
    197200C
    198       DO 110 I=1,LON
    199 C
    200          IF(WGRI(I,K,L).LT.0.) THEN
     201      DO I=1,LON
     202C
     203         IF(WGRI(I,K,L)<0.) THEN
    201204           FM(I,L)=-WGRI(I,K,L)*DTZ
    202205           ALF(I)=FM(I,L)/SM(I,K,LP)
    … …  
    215218         ALF4 (I)=ALF1(I)*ALF1Q(I)
    216219C
    217  110  CONTINUE
    218 C
    219       DO 111 JV=1,NTRA
    220       DO 1110 I=1,LON
    221 C
    222          IF(WGRI(I,K,L).LT.0.) THEN
     220      END DO
     221C
     222      DO JV=1,NTRA
     223      DO I=1,LON
     224C
     225         IF(WGRI(I,K,L)<0.) THEN
    223226C
    224227           F0 (I,L,JV)=ALF (I)* ( S0(I,K,LP,JV)-ALF1(I)*
    … …  
    273276         ENDIF
    274277C
    275  1110 CONTINUE
    276  111  CONTINUE
    277 C
    278  11   CONTINUE
     278      END DO
     279      END DO
     280C
     281      END DO
    279282C
    280283C  puts the temporary moments Fi into appropriate neighboring boxes
    281284C
    282       DO 12 L=1,NIV-1
     285      DO L=1,NIV-1
    283286      LP=L+1
    284287C
    285       DO 120 I=1,LON
    286 C
    287          IF(WGRI(I,K,L).LT.0.) THEN
     288      DO I=1,LON
     289C
     290         IF(WGRI(I,K,L)<0.) THEN
    288291           SM(I,K,L)=SM(I,K,L)+FM(I,L)
    289292           ALF(I)=FM(I,L)/SM(I,K,L)
    … …  
    299302         ALF3(I)=ALF1(I)-ALF(I)
    300303C
    301  120  CONTINUE
    302 C
    303       DO 121 JV=1,NTRA
    304       DO 1210 I=1,LON
    305 C
    306          IF(WGRI(I,K,L).LT.0.) THEN
     304      END DO
     305C
     306      DO JV=1,NTRA
     307      DO I=1,LON
     308C
     309         IF(WGRI(I,K,L)<0.) THEN
    307310C
    308311           TEMPTM=-ALF(I)*S0(I,K,L,JV)+ALF1(I)*F0(I,L,JV)
    … …  
    342345         ENDIF
    343346C
    344  1210 CONTINUE
    345  121  CONTINUE
    346 C
    347  12   CONTINUE
     347      END DO
     348      END DO
     349C
     350      END DO
    348351C
    349352C  fin de la boucle principale sur les latitudes
    350353C
    351  1    CONTINUE
     354      END DO
    352355C
    353356      DO l = 1,llm

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

    r1952 r5079  
    5858200   CONTINUE
    5959      incre = incre + 1
    60       DO 99999 j = 1, kylat
    61       DO 99999 i = 1, kxlon
    62       IF (pfild(i,j).GT. zwmsk) THEN
     60      DO j = 1, kylat
     61      DO i = 1, kxlon
     62      IF (pfild(i,j)> 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 .EQ. 1) ideb = 4
    92          IF (j .EQ. kylat) ifin = 6
     91         IF (j == 1) ideb = 4
     92         IF (j == kylat) ifin = 6
    9393C
    9494C* Account for periodicity in longitude
    9595C
    9696         IF (ldper) THEN
    97             IF (i .EQ. kxlon) THEN
     97            IF (i == kxlon) THEN
    9898               ix(3) = 1
    9999               ix(6) = 1
    100100               ix(9) = 1
    101             ELSE IF (i .EQ. 1) THEN
     101            ELSE IF (i == 1) THEN
    102102               ix(1) = kxlon
    103103               ix(4) = kxlon
    … …  
    105105            ENDIF
    106106         ELSE
    107             IF (i .EQ. 1) THEN
     107            IF (i == 1) THEN
    108108               ix(1) = i
    109109               ix(2) = i + 1
    … …  
    113113               ix(6) = i + 1
    114114            ENDIF
    115             IF (i .EQ. kxlon) THEN
     115            IF (i == kxlon) THEN
    116116               ix(1) = i -1
    117117               ix(2) = i
    … …  
    122122            ENDIF
    123123C
    124             IF (i .EQ. 1 .OR. i .EQ. kxlon) THEN
     124            IF (i == 1 .OR. i == 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 .EQ. 1) ideb = 3
    135                IF (j .EQ. kylat) ifin = 4
     134               IF (j == 1) ideb = 3
     135               IF (j == kylat) ifin = 4
    136136            ENDIF
    137137         ENDIF ! end for ldper test
    … …  
    139139C* Find unmasked neighbors
    140140C
    141          DO 230 k = ideb, ifin
     141         DO k = ideb, ifin
    142142            zmask(k) = 0.
    143143            ilon = ix(k)
    144144            jlat = jy(k)
    145             IF (pfild(ilon,jlat) .LT. zwmsk) THEN
     145            IF (pfild(ilon,jlat) < zwmsk) THEN
    146146               zmask(k) = 1.
    147147               inbor = inbor + 1
    148148            ENDIF
    149  230     CONTINUE
     149      END DO
    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 .GE. knbor) THEN
     154         IF (inbor >= knbor) THEN
    155155            pwork(i,j) = 0.
    156156            DO k = ideb, ifin
    … …  
    163163C
    164164      ENDIF
    165 99999 CONTINUE
     165      END DO
     166      END DO
    166167C
    167168C*    3. Writing back unmasked field in pfild
    … …  
    176177      DO j = 1, kylat
    177178      DO i = 1, kxlon
    178          IF (pwork(i,j) .GT. zwmsk) idoit = idoit + 1
     179         IF (pwork(i,j) > zwmsk) idoit = idoit + 1
    179180         pfild(i,j) = pwork(i,j)
    180181      ENDDO
    181182      ENDDO
    182183c
    183       IF (idoit .ne. 0) GOTO 200
     184      IF (idoit /= 0) GOTO 200
    184185ccc      PRINT*, "Number of extrapolation steps incre =", incre
    185186c

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

    r2197 r5079  
    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.eq.1) then
     317      if(NSTEP==1) then
    318318c
    319319      write(6,*) '------------------------------------ '
    … …  
    325325C
    326326C controles sur les parametres
    327       if(NLAY.LT.6) then
     327      if(NLAY<6) then
    328328        write(6,*) 'NLAY must be >= 6'
    329329        stop
    330330      endif
    331       if (JNP.LT.NLAY) then
     331      if (JNP<NLAY) then
    332332         write(6,*) 'JNP must be >= NLAY'
    333333        stop
    334334      endif
    335335      IMRD2=mod(IMR,2)
    336       if (j1.eq.2.and.IMRD2.NE.0) then
     336      if (j1==2.and.IMRD2/=0) then
    337337         write(6,*) 'if j1=2 IMR must be an even integer'
    338338        stop
    … …  
    340340
    341341C
    342       if(Jmax.lt.JNP .or. Kmax.lt.NLAY) then
     342      if(Jmax<JNP .or. Kmax<NLAY) then
    343343        write(6,*) 'Jmax or Kmax is too small'
    344344        stop
    … …  
    354354      DP =    PI / REAL(JMR)
    355355C
    356       if(IGD.eq.0) then
     356      if(IGD==0) then
    357357C Compute analytic cosine at cell edges
    358358            call cosa(cosp,cose,JNP,PI,DP)
    … …  
    362362      endif
    363363C
    364       do 15 J=2,JMR
    365 15    acosp(j) = 1. / cosp(j)
     364      do J=2,JMR
     365      acosp(j) = 1. / cosp(j)
     366      END DO
    366367C
    367368C Inverse of the Scaled polar cap area.
    … …  
    372373      endif
    373374C
    374       if(NDT0 .ne. NDT) then
     375      if(NDT0 /= NDT) then
    375376      DT   = NDT
    376377      NDT0 = NDT
    377378
    378         if(Umax .lt. 180.) then
     379        if(Umax < 180.) then
    379380         write(6,*) 'Umax may be too small!'
    380381        endif
    … …  
    382383      MaxDT = DP*AE / abs(Umax) + 0.5
    383384      write(6,*)'Largest time step for max(V)=',Umax,' is ',MaxDT
    384       if(MaxDT .lt. abs(NDT)) then
     385      if(MaxDT < abs(NDT)) then
    385386            write(6,*) 'Warning!!! NDT maybe too large!'
    386387      endif
    387388C
    388       if(CR1.ge.0.95) then
     389      if(CR1>=0.95) then
    389390      JS0 = 0
    390391      JN0 = 0
    … …  
    429430         
    430431C
    431       if(j1.ne.2) then
    432       DO 40 IC=1,NC
    433       DO 40 L=1,NLAY
    434       DO 40 I=1,IMR
     432      if(j1/=2) then
     433      DO IC=1,NC
     434      DO L=1,NLAY
     435      DO I=1,IMR
    435436      Q(I,  2,L,IC) = Q(I,  1,L,IC)
    436 40    Q(I,JMR,L,IC) = Q(I,JNP,L,IC)
     437      Q(I,JMR,L,IC) = Q(I,JNP,L,IC)
     438      END DO
     439      END DO
     440      END DO
    437441      endif
    438442C
    439443C Compute "tracer density"
    440       DO 550 IC=1,NC
    441       DO 44 k=1,NLAY
    442       DO 44 j=1,JNP
    443       DO 44 i=1,IMR
    444 44    DQ(i,j,k,IC) = Q(i,j,k,IC)*delp1(i,j,k)
    445 550     continue
    446 C
    447       do 1500 k=1,NLAY
    448 C
    449       if(IGD.eq.0) then
     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
     453C
     454      do k=1,NLAY
     455C
     456      if(IGD==0) then
    450457C Convert winds on A-Grid to Courant # on C-Grid.
    451458      call A2C(U(1,1,k),V(1,1,k),IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)
    452459      else
    453460C Convert winds on C-grid to Courant #
    454       do 45 j=j1,j2
    455       do 45 i=2,IMR
    456 45    CRX(i,J) = dtdx(j)*U(i-1,j,k)
     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
    457466   
    458467C
    459       do 50 j=j1,j2
    460 50    CRX(1,J) = dtdx(j)*U(IMR,j,k)
    461 C
    462       do 55 i=1,IMR*JMR
    463 55    CRY(i,2) = DTDY*V(i,1,k)
     468      do j=j1,j2
     469      CRX(1,J) = dtdx(j)*U(IMR,j,k)
     470      END DO
     471C
     472      do i=1,IMR*JMR
     473      CRY(i,2) = DTDY*V(i,1,k)
     474      END DO
    464475      endif
    465476C     
    … …  
    470481      do j=JS0,j1+1,-1
    471482      do i=1,IMR
    472       if(abs(CRX(i,j)).GT.1.) then
     483      if(abs(CRX(i,j))>1.) then
    473484            JS = j
    474485            go to 2222
    … …  
    480491      do j=JN0,j2-1
    481492      do i=1,IMR
    482       if(abs(CRX(i,j)).GT.1.) then
     493      if(abs(CRX(i,j))>1.) then
    483494            JN = j
    484495            go to 2233
    … …  
    4884992233  continue
    489500C
    490       if(j1.ne.2) then           ! Enlarged polar cap.
     501      if(j1/=2) then           ! Enlarged polar cap.
    491502      do i=1,IMR
    492503      DPI(i,  2,k) = 0.
    … …  
    505516      enddo
    506517C
    507       do 95 j=j1,j2
    508       DO 95 i=1,IMR
    509 95    DPI(i,j,k) = (ymass(i,j) - ymass(i,j+1)) * acosp(j)
     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
    510523C
    511524C Poles
    … …  
    536549      enddo
    537550C
    538       do 110 j=j1,j2
    539       DO 110 i=1,IMR
    540 110   xmass(i,j) = PU(i,j)*CRX(i,j)
    541 C
    542       DO 120 j=j1,j2
    543       DO 120 i=1,IMR-1
    544 120   DPI(i,j,k) = DPI(i,j,k) + xmass(i,j) - xmass(i+1,j)
    545 C
    546       DO 130 j=j1,j2
    547 130   DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j)
     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
     556C
     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
     562C
     563      DO j=j1,j2
     564      DPI(IMR,j,k) = DPI(IMR,j,k) + xmass(IMR,j) - xmass(1,j)
     565      END DO
    548566C
    549567      DO j=j1,j2
    … …  
    569587      enddo
    570588C
    571       if(j1.eq.2) then
     589      if(j1==2) then
    572590        IMH = IMR/2
    573591      do i=1,IMH
    … …  
    582600C
    583601C ****6***0*********0*********0*********0*********0*********0**********72
    584       do 1000 IC=1,NC
     602      do IC=1,NC
    585603C
    586604      do i=1,IMJM
    … …  
    590608C
    591609C E-W advective cross term
    592       do 250 j=J1,J2
    593       if(J.GT.JS  .and. J.LT.JN) GO TO 250
     610      do j=J1,J2
     611      if(J>JS  .and. J<JN) GO TO 250
    594612C
    595613      do i=1,IMR
    … …  
    602620      enddo
    603621C
    604       DO 230 i=1,IMR
     622      DO i=1,IMR
    605623      iu = UA(i,j)
    606624      ru = UA(i,j) - iu
    607625      iiu = i-iu
    608       if(UA(i,j).GE.0.) then
     626      if(UA(i,j)>=0.) then
    609627      wk1(i,j,1) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))
    610628      else
    … …  
    612630      endif
    613631      wk1(i,j,1) = wk1(i,j,1) - qtmp(i)
    614 230   continue
     632      END DO
    615633250   continue
    616 C
    617       if(JN.ne.0) then
     634      END DO
     635C
     636      if(JN/=0) then
    618637      do j=JS+1,JN-1
    619638C
    … …  
    645664        if(cross) then
    646665C Add cross terms in the vertical direction.
    647         if(IORD .GE. 2) then
     666        if(IORD >= 2) then
    648667                iad = 2
    649668        else
    … …  
    651670        endif
    652671C
    653         if(JORD .GE. 2) then
     672        if(JORD >= 2) then
    654673                jad = 2
    655674        else
    … …  
    671690     &  DC2,ymass,WK1(1,1,3),wk1(1,1,4),WK1(1,1,5),WK1(1,1,6),JORD)
    672691C
    673 1000  continue
    674 1500  continue
     692      END DO
     693      END DO
    675694C
    676695C ******* Compute vertical mass flux (same unit as PS) ***********
    … …  
    678697C 1st step: compute total column mass CONVERGENCE.
    679698C
    680       do 320 j=1,JNP
    681       do 320 i=1,IMR
    682 320   CRY(i,j) = DPI(i,j,1)
    683 C
    684       do 330 k=2,NLAY
    685       do 330 j=1,JNP
    686       do 330 i=1,IMR
     699      do j=1,JNP
     700      do i=1,IMR
     701      CRY(i,j) = DPI(i,j,1)
     702      END DO
     703      END DO
     704C
     705      do k=2,NLAY
     706      do j=1,JNP
     707      do i=1,IMR
    687708      CRY(i,j)  = CRY(i,j) + DPI(i,j,k)
    688 330   continue
    689 C
    690       do 360 j=1,JNP
    691       do 360 i=1,IMR
     709      END DO
     710      END DO
     711      END DO
     712C
     713      do j=1,JNP
     714      do i=1,IMR
    692715C
    693716C 2nd step: compute PS2 (PS at n+1) using the hydrostatic assumption.
    … …  
    700723      W(i,j,1) = DPI(i,j,1) - DBK(1)*CRY(i,j)
    701724      W(i,j,NLAY) = 0.
    702 360   continue
    703 C
    704       do 370 k=2,NLAY-1
    705       do 370 j=1,JNP
    706       do 370 i=1,IMR
     725      END DO
     726      END DO
     727C
     728      do k=2,NLAY-1
     729      do j=1,JNP
     730      do i=1,IMR
    707731      W(i,j,k) = W(i,j,k-1) + DPI(i,j,k) - DBK(k)*CRY(i,j)
    708 370   continue
    709 C
    710       DO 380 k=1,NLAY
    711       DO 380 j=1,JNP
    712       DO 380 i=1,IMR
     732      END DO
     733      END DO
     734      END DO
     735C
     736      DO k=1,NLAY
     737      DO j=1,JNP
     738      DO i=1,IMR
    713739      delp2(i,j,k) = DAP(k) + DBK(k)*PS2(i,j)
    714 380   continue
     740      END DO
     741      END DO
     742      END DO
    715743C
    716744        KRD = max(3, KORD)
    717       do 4000 IC=1,NC
     745      do IC=1,NC
    718746C
    719747C****6***0*********0*********0*********0*********0*********0**********72
    … …  
    738766      enddo
    739767C     
    740       if(j1.ne.2) then
    741       DO 400 k=1,NLAY
    742       DO 400 I=1,IMR
    743 c     j=1 c'est le pôle Sud, j=JNP c'est le pôle Nord
     768      if(j1/=2) then
     769      DO k=1,NLAY
     770      DO I=1,IMR
     771c     j=1 c'est le p�le Sud, j=JNP c'est le p�le Nord
    744772      Q(I,  2,k,IC) = Q(I,  1,k,IC)
    745773      Q(I,JMR,k,IC) = Q(I,JNP,k,IC)
    746 400   CONTINUE
    747       endif
    748 4000  continue
    749 C
    750       if(j1.ne.2) then
    751       DO 5000 k=1,NLAY
    752       DO 5000 i=1,IMR
     774      END DO
     775      END DO
     776      endif
     777      END DO
     778C
     779      if(j1/=2) then
     780      DO k=1,NLAY
     781      DO i=1,IMR
    753782      W(i,  2,k) = W(i,  1,k)
    754783      W(i,JMR,k) = W(i,JNP,k)
    755 5000  continue
     784      END DO
     785      END DO
    756786      endif
    757787C
    … …  
    785815C ****6***0*********0*********0*********0*********0*********0**********72
    786816C
    787       do 1000 k=1,NLAYM1
    788       do 1000 i=1,IMJM
     817      do k=1,NLAYM1
     818      do i=1,IMJM
    789819      DQDT(i,1,k) = P(i,1,k+1) - P(i,1,k)
    790 1000  continue
    791 C
    792       DO 1220 k=2,NLAYM1
    793       DO 1220 I=1,IMJM   
     820      END DO
     821      END DO
     822C
     823      DO k=2,NLAYM1
     824      DO I=1,IMJM
    794825       c0 =  delp(i,1,k) / (delp(i,1,k-1)+delp(i,1,k)+delp(i,1,k+1))
    795826       c1 = (delp(i,1,k-1)+0.5*delp(i,1,k))/(delp(i,1,k+1)+delp(i,1,k))   
    … …  
    799830      Qmin = P(i,1,k) - min(P(i,1,k-1),P(i,1,k),P(i,1,k+1))
    800831      DC(i,1,k) = sign(min(abs(tmp),Qmax,Qmin), tmp)   
    801 1220  CONTINUE
     832      END DO
     833      END DO
    802834     
    803835C     
    … …  
    806838C ****6***0*********0*********0*********0*********0*********0**********72
    807839C
    808       DO 2000 j=1,JNP
    809       if((j.eq.2 .or. j.eq.JMR) .and. j1.ne.2) goto 2000
     840      DO j=1,JNP
     841      if((j==2 .or. j==JMR) .and. j1/=2) goto 2000
    810842C
    811843      DO k=1,NLAY
    … …  
    828860C
    829861C First guess top edge value
    830       DO 10 i=1,IMR
     862      DO i=1,IMR
    831863C three-cell PPM
    832864C Compute a,b, and c of q = aP**2 + bP + c using cell averages and delp
    … …  
    840872C
    841873C Check if change sign
    842       if(wk1(i,1)*AL(i,1).le.0.) then
     874      if(wk1(i,1)*AL(i,1)<=0.) then
    843875                 AL(i,1) = 0.
    844876             flux(i,1) = 0.
    … …  
    846878             flux(i,1) =  wk1(i,1) - AL(i,1)
    847879        endif
    848 10    continue
     880      END DO
    849881C
    850882C Bottom
    851       DO 15 i=1,IMR
     883      DO i=1,IMR
    852884C 2-cell PPM with zero gradient right at the surface
    853885C
    … …  
    856888      AR(i,NLAY) = wk1(i,NLAY) + fct
    857889      AL(i,NLAY) = wk1(i,NLAY) - (fct+fct)
    858       if(wk1(i,NLAY)*AR(i,NLAY).le.0.) AR(i,NLAY) = 0.
     890      if(wk1(i,NLAY)*AR(i,NLAY)<=0.) AR(i,NLAY) = 0.
    859891      flux(i,NLAY) = AR(i,NLAY) -  wk1(i,NLAY)
    860 15    continue
     892      END DO
    861893     
    862894C
    … …  
    865897C****6***0*********0*********0*********0*********0*********0**********72
    866898C
    867       DO 14 k=3,NLAYM1
    868       DO 12 i=1,IMR
     899      DO k=3,NLAYM1
     900      DO i=1,IMR
    869901      c1 =  DQDT(i,j,k-1)*wk2(i,k-1) / (wk2(i,k-1)+wk2(i,k))
    870902      c2 =  2. / (wk2(i,k-2)+wk2(i,k-1)+wk2(i,k)+wk2(i,k+1))
    … …  
    875907     &          wk2(i,k-1)*A1*flux(i,k)  )
    876908C      print *,'AL1',i,k, AL(i,k)
    877 12    CONTINUE
    878 14    continue
    879 C
    880       do 20 i=1,IMR*NLAYM1
     909      END DO
     910      END DO
     911C
     912      do i=1,IMR*NLAYM1
    881913      AR(i,1) = AL(i,2)
    882914C      print *,'AR1',i,AR(i,1)
    883 20    continue
    884 C
    885       do 30 i=1,IMR*NLAY
     915      END DO
     916C
     917      do i=1,IMR*NLAY
    886918      A6(i,1) = 3.*(wk1(i,1)+wk1(i,1) - (AL(i,1)+AR(i,1)))
    887919C      print *,'A61',i,A6(i,1)
    888 30    continue
     920      END DO
    889921C
    890922C****6***0*********0*********0*********0*********0*********0**********72
    … …  
    895927C
    896928C Interior depending on KORD
    897       if(LMT.LE.2)
     929      if(LMT<=2)
    898930     &  call lmtppm(flux(1,2),A6(1,2),AR(1,2),AL(1,2),wk1(1,2),
    899931     &              IMR*(NLAY-2),LMT)
    … …  
    901933C****6***0*********0*********0*********0*********0*********0**********72
    902934C
    903       DO 140 i=1,IMR*NLAYM1
    904       IF(wz2(i,1).GT.0.) then
     935      DO i=1,IMR*NLAYM1
     936      IF(wz2(i,1)>0.) then
    905937        CM = wz2(i,1) / wk2(i,1)
    906938        flux(i,2) = AR(i,1)+0.5*CM*(AL(i,1)-AR(i,1)+A6(i,1)*(1.-R23*CM))
    … …  
    912944C        print *,'test2',i, AL(i,2),AR(i,2),A6(i,2),R23
    913945      endif
    914 140   continue
    915 C
    916       DO 250 i=1,IMR*NLAYM1
     946      END DO
     947C
     948      DO i=1,IMR*NLAYM1
    917949      flux(i,2) = wz2(i,1) * flux(i,2)
    918 250   continue
    919 C
    920       do 350 i=1,IMR
     950      END DO
     951C
     952      do i=1,IMR
    921953      DQ(i,j,   1) = DQ(i,j,   1) - flux(i,   2)
    922954      DQ(i,j,NLAY) = DQ(i,j,NLAY) + flux(i,NLAY)
    923 350   continue
    924 C
    925       do 360 k=2,NLAYM1
    926       do 360 i=1,IMR
    927 360   DQ(i,j,k) = DQ(i,j,k) + flux(i,k) - flux(i,k+1)
     955      END DO
     956C
     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
    9289622000  continue
     963      END DO
    929964      return
    930965      end
    … …  
    950985      j2vl = j2-jvan
    951986C
    952       do 1310 j=j1,j2
     987      do j=j1,j2
    953988C
    954989      do i=1,IMR
    … …  
    956991      enddo
    957992C
    958       if(j.ge.JN .or. j.le.JS) goto 2222
     993      if(j>=JN .or. j<=JS) goto 2222
    959994C ************* Eulerian **********
    960995C
    … …  
    964999      qtmp(IMP+1) = q(2,J)
    9651000C
    966       IF(IORD.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN
    967       DO 1406 i=1,IMR
     1001      IF(IORD==1 .or. j==j1 .or. j==j2) THEN
     1002      DO i=1,IMR
    9681003      iu = REAL(i) - uc(i,j)
    969 1406  fx1(i) = qtmp(iu)
     1004      fx1(i) = qtmp(iu)
     1005      END DO
    9701006      ELSE
    9711007      call xmist(IMR,IML,Qtmp,DC)
    9721008      DC(0) = DC(IMR)
    9731009C
    974       if(IORD.eq.2 .or. j.le.j1vl .or. j.ge.j2vl) then
    975       DO 1408 i=1,IMR
     1010      if(IORD==2 .or. j<=j1vl .or. j>=j2vl) then
     1011      DO i=1,IMR
    9761012      iu = REAL(i) - uc(i,j)
    977 1408  fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j))
     1013      fx1(i) = qtmp(iu) + DC(iu)*(sign(1.,uc(i,j))-uc(i,j))
     1014      END DO
    9781015      else
    9791016      call fxppm(IMR,IML,UC(1,j),Qtmp,DC,fx1,IORD)
    … …  
    9821019      ENDIF
    9831020C
    984       DO 1506 i=1,IMR
    985 1506  fx1(i) = fx1(i)*xmass(i,j)
     1021      DO i=1,IMR
     1022      fx1(i) = fx1(i)*xmass(i,j)
     1023      END DO
    9861024C
    9871025      goto 1309
    … …  
    9961034      enddo
    9971035C
    998       IF(IORD.eq.1 .or. j.eq.j1. or. j.eq.j2) THEN
    999       DO 1306 i=1,IMR
     1036      IF(IORD==1 .or. j==j1 .or. j==j2) THEN
     1037      DO i=1,IMR
    10001038      itmp = INT(uc(i,j))
    10011039      ISAVE(i) = i - itmp
    10021040      iu = i - uc(i,j)
    1003 1306  fx1(i) = (uc(i,j) - itmp)*qtmp(iu)
     1041      fx1(i) = (uc(i,j) - itmp)*qtmp(iu)
     1042      END DO
    10041043      ELSE
    10051044      call xmist(IMR,IML,Qtmp,DC)
    … …  
    10101049      enddo
    10111050C
    1012       DO 1307 i=1,IMR
     1051      DO i=1,IMR
    10131052      itmp = INT(uc(i,j))
    10141053      rut  = uc(i,j) - itmp
    10151054      ISAVE(i) = i - itmp
    10161055      iu = i - uc(i,j)
    1017 1307  fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut))
     1056      fx1(i) = rut*(qtmp(iu) + DC(iu)*(sign(1.,rut) - rut))
     1057      END DO
    10181058      ENDIF
    10191059C
    1020       do 1308 i=1,IMR
    1021       IF(uc(i,j).GT.1.) then
     1060      do i=1,IMR
     1061      IF(uc(i,j)>1.) then
    10221062CDIR$ NOVECTOR
    10231063        do ist = ISAVE(i),i-1
    10241064        fx1(i) = fx1(i) + qtmp(ist)
    10251065        enddo
    1026       elseIF(uc(i,j).LT.-1.) then
     1066      elseIF(uc(i,j)<-1.) then
    10271067        do ist = i,ISAVE(i)-1
    10281068        fx1(i) = fx1(i) - qtmp(ist)
    … …  
    10301070CDIR$ VECTOR
    10311071      endif
    1032 1308  continue
     1072      END DO
    10331073      do i=1,IMR
    10341074      fx1(i) = PU(i,j)*fx1(i)
    … …  
    10381078C
    103910791309  fx1(IMP) = fx1(1)
    1040       DO 1215 i=1,IMR
    1041 1215  DQ(i,j) =  DQ(i,j) + fx1(i)-fx1(i+1)
     1080      DO i=1,IMR
     1081      DQ(i,j) =  DQ(i,j) + fx1(i)-fx1(i+1)
     1082      END DO
    10421083C
    10431084C ***************************************
    10441085C
    1045 1310  continue
     1086      END DO
    10461087      return
    10471088      end
    … …  
    10791120C      endif
    10801121C
    1081       DO 10 i=1,IMR
    1082 10    AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3
    1083 C
    1084       do 20 i=1,IMR-1
    1085 20    AR(i) = AL(i+1)
     1122      DO i=1,IMR
     1123      AL(i) = 0.5*(p(i-1)+p(i)) + (DC(i-1) - DC(i))*R3
     1124      END DO
     1125C
     1126      do i=1,IMR-1
     1127      AR(i) = AL(i+1)
     1128      END DO
    10861129      AR(IMR) = AL(1)
    10871130C
    1088       do 30 i=1,IMR
    1089 30    A6(i) = 3.*(p(i)+p(i)  - (AL(i)+AR(i)))
    1090 C
    1091       if(LMT.LE.2) call lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT)
     1131      do i=1,IMR
     1132      A6(i) = 3.*(p(i)+p(i)  - (AL(i)+AR(i)))
     1133      END DO
     1134C
     1135      if(LMT<=2) call lmtppm(DC(1),A6(1),AR(1),AL(1),P(1),IMR,LMT)
    10921136C
    10931137      AL(0) = AL(IMR)
    … …  
    10961140C
    10971141      DO i=1,IMR
    1098       IF(UT(i).GT.0.) then
     1142      IF(UT(i)>0.) then
    10991143      flux(i) = AR(i-1) + 0.5*UT(i)*(AL(i-1) - AR(i-1) +
    11001144     &                 A6(i-1)*(1.-R23*UT(i)) )
    … …  
    11151159      real :: tmp,pmax,pmin
    11161160C
    1117       do 10  i=1,IMR
     1161      do i=1,IMR
    11181162      tmp = R24*(8.*(p(i+1) - p(i-1)) + p(i-2) - p(i+2))
    11191163      Pmax = max(P(i-1), p(i), p(i+1)) - p(i)
    11201164      Pmin = p(i) - min(P(i-1), p(i), p(i+1))
    1121 10    DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp)
     1165      DC(i) = sign(min(abs(tmp),Pmax,Pmin), tmp)
     1166      END DO
    11221167      return
    11231168      end
    … …  
    11381183      len = IMR*(J2-J1+2)
    11391184C
    1140       if(JORD.eq.1) then
    1141       DO 1000 i=1,len
     1185      if(JORD==1) then
     1186      DO i=1,len
    11421187      JT = REAL(J1) - VC(i,J1)
    1143 1000  fx(i,j1) = p(i,JT)
     1188      fx(i,j1) = p(i,JT)
     1189      END DO
    11441190      else
    11451191   
    11461192      call ymist(IMR,JNP,j1,P,DC2,4)
    11471193C
    1148       if(JORD.LE.0 .or. JORD.GE.3) then
     1194      if(JORD<=0 .or. JORD>=3) then
    11491195   
    11501196      call fyppm(VC,P,DC2,fx,IMR,JNP,j1,j2,A6,AR,AL,JORD)
    11511197   
    11521198      else
    1153       DO 1200 i=1,len
     1199      DO i=1,len
    11541200      JT = REAL(J1) - VC(i,J1)
    1155 1200  fx(i,j1) = p(i,JT) + (sign(1.,VC(i,j1))-VC(i,j1))*DC2(i,JT)
    1156       endif
    1157       endif
    1158 C
    1159       DO 1300 i=1,len
    1160 1300  fx(i,j1) = fx(i,j1)*ymass(i,j1)
    1161 C
    1162       DO 1400 j=j1,j2
    1163       DO 1400 i=1,IMR
    1164 1400  DQ(i,j) = DQ(i,j) + (fx(i,j) - fx(i,j+1)) * acosp(j)
     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
     1205C
     1206      DO i=1,len
     1207      fx(i,j1) = fx(i,j1)*ymass(i,j1)
     1208      END DO
     1209C
     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
    11651215C
    11661216C Poles
    … …  
    11791229      enddo
    11801230C
    1181       if(j1.ne.2) then
     1231      if(j1/=2) then
    11821232      do i=1,IMR
    11831233      DQ(i,  2) = sum1
    … …  
    12011251      IJM3 = IMR*(JMR-3)
    12021252C
    1203       IF(ID.EQ.2) THEN
    1204       do 10 i=1,IMR*(JMR-1)
     1253      IF(ID==2) THEN
     1254      do i=1,IMR*(JMR-1)
    12051255      tmp = 0.25*(p(i,3) - p(i,1))
    12061256      Pmax = max(p(i,1),p(i,2),p(i,3)) - p(i,2)
    12071257      Pmin = p(i,2) - min(p(i,1),p(i,2),p(i,3))
    12081258      DC(i,2) = sign(min(abs(tmp),Pmin,Pmax),tmp)
    1209 10    CONTINUE
     1259      END DO
    12101260      ELSE
    1211       do 12 i=1,IMH
     1261      do i=1,IMH
    12121262C J=2
    12131263      tmp = (8.*(p(i,3) - p(i,1)) + p(i+IMH,2) - p(i,4))*R24
    … …  
    12201270      Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP))
    12211271      DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp)
    1222 12    CONTINUE
    1223       do 14 i=IMH+1,IMR
     1272      END DO
     1273      do i=IMH+1,IMR
    12241274C J=2
    12251275      tmp = (8.*(p(i,3) - p(i,1)) + p(i-IMH,2) - p(i,4))*R24
    … …  
    12321282      Pmin = p(i,JMR) - min(p(i,JMR-1),p(i,JMR),p(i,JNP))
    12331283      DC(i,JMR) = sign(min(abs(tmp),Pmin,Pmax),tmp)
    1234 14    CONTINUE
    1235 C
    1236       do 15 i=1,IJM3
     1284      END DO
     1285C
     1286      do i=1,IJM3
    12371287      tmp = (8.*(p(i,4) - p(i,2)) + p(i,1) - p(i,5))*R24
    12381288      Pmax = max(p(i,2),p(i,3),p(i,4)) - p(i,3)
    12391289      Pmin = p(i,3) - min(p(i,2),p(i,3),p(i,4))
    12401290      DC(i,3) = sign(min(abs(tmp),Pmin,Pmax),tmp)
    1241 15    CONTINUE
     1291      END DO
    12421292      ENDIF
    12431293C
    1244       if(j1.ne.2) then
     1294      if(j1/=2) then
    12451295      do i=1,IMR
    12461296      DC(i,1) = 0.
    … …  
    12501300C Determine slopes in polar caps for scalars!
    12511301C
    1252       do 13 i=1,IMH
     1302      do i=1,IMH
    12531303C South
    12541304      tmp = 0.25*(p(i,2) - p(i+imh,2))
    … …  
    12611311      Pmin = p(i,JNP) - min(p(i+imh,JMR),p(i,jnp), p(i,JMR))
    12621312      DC(i,JNP) = sign(min(abs(tmp),Pmax,pmin),tmp)
    1263 13    continue
    1264 C
    1265       do 25 i=imh+1,IMR
     1313      END DO
     1314C
     1315      do i=imh+1,IMR
    12661316      DC(i,  1) =  - DC(i-imh,  1)
    12671317      DC(i,JNP) =  - DC(i-imh,JNP)
    1268 25    continue
     1318      END DO
    12691319      endif
    12701320      return
    … …  
    13081358      LMT = JORD - 3     
    13091359C
    1310       DO 10 i=1,IMR*JMR       
     1360      DO i=1,IMR*JMR
    13111361      AL(i,2) = 0.5*(p(i,1)+p(i,2)) + (DC(i,1) - DC(i,2))*R3
    13121362      AR(i,1) = AL(i,2)
    1313 10    CONTINUE
     1363      END DO
    13141364C
    13151365CPoles:
    … …  
    13311381     
    13321382           
    1333       do 30 i=1,len
    1334 30    A6(i,j11) = 3.*(p(i,j11)+p(i,j11)  - (AL(i,j11)+AR(i,j11)))
    1335 C
    1336       if(LMT.le.2) call lmtppm(DC(1,j11),A6(1,j11),AR(1,j11)
     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
     1386C
     1387      if(LMT<=2) call lmtppm(DC(1,j11),A6(1,j11),AR(1,j11)
    13371388     &                       ,AL(1,j11),P(1,j11),len,LMT)
    13381389C
    13391390     
    1340       DO 140 i=1,IMJM1
    1341       IF(VC(i,j1).GT.0.) then
     1391      DO i=1,IMJM1
     1392      IF(VC(i,j1)>0.) then
    13421393      flux(i,j1) = AR(i,j11) + 0.5*VC(i,j1)*(AL(i,j11) - AR(i,j11) +
    13431394     &                         A6(i,j11)*(1.-R23*VC(i,j1)) )
    … …  
    13461397     &                        A6(i,j1)*(1.+R23*VC(i,j1)))
    13471398      endif
    1348 140   continue
     1399      END DO
    13491400      return
    13501401      end
    … …  
    13781429c        write(*,*) 'toto 1'
    13791430C --------------------------------
    1380       IF(IAD.eq.2) then
     1431      IF(IAD==2) then
    13811432      do j=j1-1,j2+1
    13821433      do i=1,IMR
    … …  
    13951446c      write(*,*) 'toto 2'
    13961447C
    1397       ELSEIF(IAD.eq.1) then
     1448      ELSEIF(IAD==1) then
    13981449        do j=j1-1,j2+1
    13991450      do i=1,imr
    … …  
    14041455      ENDIF
    14051456C
    1406         if(j1.ne.2) then
     1457        if(j1/=2) then
    14071458        sum1 = 0.
    14081459        sum2 = 0.
    … …  
    14481499C
    14491500        JMR = JNP-1
    1450       do 1309 j=j1,j2
    1451       if(J.GT.JS  .and. J.LT.JN) GO TO 1309
     1501      do j=j1,j2
     1502      if(J>JS  .and. J<JN) GO TO 1309
    14521503C
    14531504      do i=1,IMR
    … …  
    14601511      enddo
    14611512C
    1462       IF(IAD.eq.2) THEN
     1513      IF(IAD==2) THEN
    14631514      DO i=1,IMR
    14641515      IP = NINT(UA(i,j))
    … …  
    14691520      adx(i,j) = qtmp(ip) + ru*(a1*ru + b1)
    14701521      enddo
    1471       ELSEIF(IAD.eq.1) then
     1522      ELSEIF(IAD==1) then
    14721523      DO i=1,IMR
    14731524      iu = UA(i,j)
    14741525      ru = UA(i,j) - iu
    14751526      iiu = i-iu
    1476       if(UA(i,j).GE.0.) then
     1527      if(UA(i,j)>=0.) then
    14771528      adx(i,j) = qtmp(iiu)+ru*(qtmp(iiu-1)-qtmp(iiu))
    14781529      else
    … …  
    14861537      enddo
    148715381309  continue
     1539      END DO
    14881540C
    14891541C Eulerian upwind
    … …  
    14981550      qtmp(IMR+1) = p(1,J)
    14991551C
    1500       IF(IAD.eq.2) THEN
     1552      IF(IAD==2) THEN
    15011553      qtmp(-1)     = p(IMR-1,J)
    15021554      qtmp(IMR+2) = p(2,J)
    … …  
    15091561      adx(i,j) = qtmp(ip)- p(i,j) + ru*(a1*ru + b1)
    15101562      enddo
    1511       ELSEIF(IAD.eq.1) then
     1563      ELSEIF(IAD==1) then
    15121564C 1st order
    15131565      DO i=1,IMR
    … …  
    15181570      enddo
    15191571C
    1520         if(j1.ne.2) then
     1572        if(j1/=2) then
    15211573      do i=1,IMR
    15221574      adx(i,  2) = 0.
    … …  
    15541606      REAL da1,da2,a6da,fmin
    15551607C
    1556       if(LMT.eq.0) then
     1608      if(LMT==0) then
    15571609C Full constraint
    1558       do 100 i=1,IM
    1559       if(DC(i).eq.0.) then
     1610      do i=1,IM
     1611      if(DC(i)==0.) then
    15601612            AR(i) = p(i)
    15611613            AL(i) = p(i)
    … …  
    15651617      da2  = da1**2
    15661618      A6DA = A6(i)*da1
    1567       if(A6DA .lt. -da2) then
     1619      if(A6DA < -da2) then
    15681620            A6(i) = 3.*(AL(i)-p(i))
    15691621            AR(i) = AL(i) - A6(i)
    1570       elseif(A6DA .gt. da2) then
     1622      elseif(A6DA > da2) then
    15711623            A6(i) = 3.*(AR(i)-p(i))
    15721624            AL(i) = AR(i) - A6(i)
    15731625      endif
    15741626      endif
    1575 100   continue
    1576       elseif(LMT.eq.1) then
     1627      END DO
     1628      elseif(LMT==1) then
    15771629C Semi-monotonic constraint
    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
     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
    15811633            AR(i) = p(i)
    15821634            AL(i) = p(i)
    15831635            A6(i) = 0.
    1584       elseif(AR(i) .gt. AL(i)) then
     1636      elseif(AR(i) > AL(i)) then
    15851637            A6(i) = 3.*(AL(i)-p(i))
    15861638            AR(i) = AL(i) - A6(i)
    … …  
    15901642      endif
    15911643150   continue
    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
     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
    15951648      fmin = p(i) + 0.25*(AR(i)-AL(i))**2/A6(i) + A6(i)*R12
    1596       if(fmin.ge.0.) go to 250
    1597       if(p(i).lt.AR(i) .and. p(i).lt.AL(i)) then
     1649      if(fmin>=0.) go to 250
     1650      if(p(i)<AR(i) .and. p(i)<AL(i)) then
    15981651            AR(i) = p(i)
    15991652            AL(i) = p(i)
    16001653            A6(i) = 0.
    1601       elseif(AR(i) .gt. AL(i)) then
     1654      elseif(AR(i) > AL(i)) then
    16021655            A6(i) = 3.*(AL(i)-p(i))
    16031656            AR(i) = AL(i) - A6(i)
    … …  
    16071660      endif
    16081661250   continue
     1662      END DO
    16091663      endif
    16101664      return
    … …  
    16171671      integer i,j
    16181672C
    1619       do 35 j=j1,j2
    1620       do 35 i=2,IMR
    1621 35    CRX(i,J) = dtdx5(j)*(U(i,j)+U(i-1,j))
    1622 C
    1623       do 45 j=j1,j2
    1624 45    CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j))
    1625 C
    1626       do 55 i=1,IMR*JMR
    1627 55    CRY(i,2) = DTDY5*(V(i,2)+V(i,1))
     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
     1678C
     1679      do j=j1,j2
     1680      CRX(1,J) = dtdx5(j)*(U(1,j)+U(IMR,j))
     1681      END DO
     1682C
     1683      do i=1,IMR*JMR
     1684      CRY(i,2) = DTDY5*(V(i,2)+V(i,1))
     1685      END DO
    16281686      return
    16291687      end
    … …  
    16361694      real ph5
    16371695      JMR = JNP-1
    1638       do 55 j=2,JNP
     1696      do j=2,JNP
    16391697        ph5  =  -0.5*PI + (REAL(J-1)-0.5)*DP
    1640 55      cose(j) = cos(ph5)
     1698        cose(j) = cos(ph5)
     1699      END DO
    16411700C
    16421701      JEQ = (JNP+1) / 2
    1643       if(JMR .eq. 2*(JMR/2) ) then
     1702      if(JMR == 2*(JMR/2) ) then
    16441703      do j=JNP, JEQ+1, -1
    16451704       cose(j) =  cose(JNP+2-j)
    … …  
    16531712      endif
    16541713C
    1655       do 66 j=2,JMR
    1656 66    cosp(j) = 0.5*(cose(j)+cose(j+1))
     1714      do j=2,JMR
     1715      cosp(j) = 0.5*(cose(j)+cose(j+1))
     1716      END DO
    16571717      cosp(1) = 0.
    16581718      cosp(JNP) = 0.
    … …  
    16681728C
    16691729      phi = -0.5*PI
    1670       do 55 j=2,JNP-1
     1730      do j=2,JNP-1
    16711731      phi  =  phi + DP
    1672 55    cosp(j) = cos(phi)
     1732      cosp(j) = cos(phi)
     1733      END DO
    16731734        cosp(  1) = 0.
    16741735        cosp(JNP) = 0.
    16751736C
    1676       do 66 j=2,JNP
     1737      do j=2,JNP
    16771738        cose(j) = 0.5*(cosp(j)+cosp(j-1))
    1678 66    CONTINUE
    1679 C
    1680       do 77 j=2,JNP-1
     1739      END DO
     1740C
     1741      do j=2,JNP-1
    16811742       cosp(j) = 0.5*(cose(j)+cose(j+1))
    1682 77    CONTINUE
     1743      END DO
    16831744      return
    16841745      end
    … …  
    17021763        icr = 1
    17031764      call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
    1704       if(ipy.eq.0) goto 50
     1765      if(ipy==0) goto 50
    17051766      call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
    1706       if(ipx.eq.0) goto 50
     1767      if(ipx==0) goto 50
    17071768C
    17081769      if(cross) then
    17091770      call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
    17101771      endif
    1711       if(icr.eq.0) goto 50
     1772      if(icr==0) goto 50
    17121773C
    17131774C Vertical filling...
    17141775      do i=1,len
    1715       IF( Q(i,j1,1).LT.0.) THEN
     1776      IF( Q(i,j1,1)<0.) THEN
    17161777      ip = ip + 1
    17171778          Q(i,j1,2) = Q(i,j1,2) + Q(i,j1,1)
    … …  
    17211782C
    1722178350    continue
    1723       DO 225 L = 2,NLAYM1
     1784      DO L = 2,NLAYM1
    17241785      icr = 1
    17251786C
    17261787      call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
    1727       if(ipy.eq.0) goto 225
     1788      if(ipy==0) goto 225
    17281789      call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
    1729       if(ipx.eq.0) go to 225
     1790      if(ipx==0) go to 225
    17301791      if(cross) then
    17311792      call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
    17321793      endif
    1733       if(icr.eq.0) goto 225
     1794      if(icr==0) goto 225
    17341795C
    17351796      do i=1,len
    1736       IF( Q(I,j1,L).LT.0.) THEN
     1797      IF( Q(I,j1,L)<0.) THEN
    17371798C
    17381799      ip = ip + 1
    … …  
    17491810      ENDDO
    17501811225   CONTINUE
     1812      END DO
    17511813C
    17521814C BOTTOM LAYER
    … …  
    17551817C
    17561818      call filns(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
    1757       if(ipy.eq.0) goto 911
     1819      if(ipy==0) goto 911
    17581820      call filew(q(1,1,L),qtmp,IMR,JNP,j1,j2,ipx,tiny)
    1759       if(ipx.eq.0) goto 911
     1821      if(ipx==0) goto 911
    17601822C
    17611823      call filcr(q(1,1,L),IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
    1762       if(icr.eq.0) goto 911
     1824      if(icr==0) goto 911
    17631825C
    17641826      DO  I=1,len
    1765       IF( Q(I,j1,L).LT.0.) THEN
     1827      IF( Q(I,j1,L)<0.) THEN
    17661828      ip = ip + 1
    17671829c
    … …  
    17801842911   continue
    17811843C
    1782       if(ip.gt.IMR) then
     1844      if(ip>IMR) then
    17831845      write(6,*) 'IC=',IC,' STEP=',NSTEP,
    17841846     &           ' Vertical filling pts=',ip
    17851847      endif
    17861848C
    1787       if(sum.gt.1.e-25) then
     1849      if(sum>1.e-25) then
    17881850      write(6,*) IC,NSTEP,' Mass source from the ground=',sum
    17891851      endif
    … …  
    17981860      real :: dq,dn,d0,d1,ds,d2
    17991861      icr = 0
    1800       do 65 j=j1+1,j2-1
    1801       DO 50 i=1,IMR-1
    1802       IF(q(i,j).LT.0.) THEN
     1862      do j=j1+1,j2-1
     1863      DO i=1,IMR-1
     1864      IF(q(i,j)<0.) THEN
    18031865      icr =  1
    18041866      dq  = - q(i,j)*cosp(j)
    … …  
    18161878      q(i,j) = (d2 - dq)*acosp(j) + tiny
    18171879      endif
    1818 50    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
     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
    18221884      icr =  1
    18231885      dq  = - q(i,j)*cosp(j)
    … …  
    18351897      q(i,j) = (d2 - dq)*acosp(j) + tiny
    18361898      endif
    1837 55    continue
     1899      END DO
    18381900C *****************************************
    18391901C i=1
    18401902      i=1
    1841       IF(q(i,j).LT.0.) THEN
     1903      IF(q(i,j)<0.) THEN
    18421904      icr =  1
    18431905      dq  = - q(i,j)*cosp(j)
    … …  
    18581920C i=IMR
    18591921      i=IMR
    1860       IF(q(i,j).LT.0.) THEN
     1922      IF(q(i,j)<0.) THEN
    18611923      icr =  1
    18621924      dq  = - q(i,j)*cosp(j)
    … …  
    18761938C *****************************************
    1877193965    continue
    1878 C
    1879       do i=1,IMR
    1880       if(q(i,j1).lt.0. .or. q(i,j2).lt.0.) then
     1940      END DO
     1941C
     1942      do i=1,IMR
     1943      if(q(i,j1)<0. .or. q(i,j2)<0.) then
    18811944      icr = 1
    18821945      goto 80
    … …  
    1886194980    continue
    18871950C
    1888       if(q(1,1).lt.0. .or. q(1,jnp).lt.0.) then
     1951      if(q(1,1)<0. .or. q(1,jnp)<0.) then
    18891952      icr = 1
    18901953      endif
    … …  
    19101973C
    19111974      ipy = 0
    1912       do 55 j=j1+1,j2-1
    1913       DO 55 i=1,IMR
    1914       IF(q(i,j).LT.0.) THEN
     1975      do j=j1+1,j2-1
     1976      DO i=1,IMR
     1977      IF(q(i,j)<0.) THEN
    19151978      ipy =  1
    19161979      dq  = - q(i,j)*cosp(j)
    … …  
    19281991      q(i,j) = (d2 - dq)*acosp(j) + tiny
    19291992      endif
    1930 55    continue
     1993      END DO
     1994      END DO
    19311995C
    19321996      do i=1,imr
    1933       IF(q(i,j1).LT.0.) THEN
     1997      IF(q(i,j1)<0.) THEN
    19341998      ipy =  1
    19351999      dq  = - q(i,j1)*cosp(j1)
    … …  
    19452009      j = j2
    19462010      do i=1,imr
    1947       IF(q(i,j).LT.0.) THEN
     2011      IF(q(i,j)<0.) THEN
    19482012      ipy =  1
    19492013      dq  = - q(i,j)*cosp(j)
    … …  
    19582022C
    19592023C Check Poles.
    1960       if(q(1,1).lt.0.) then
     2024      if(q(1,1)<0.) then
    19612025      dq = q(1,1)*cap1/REAL(IMR)*acosp(j1)
    19622026      do i=1,imr
    19632027      q(i,1) = 0.
    19642028      q(i,j1) = q(i,j1) + dq
    1965       if(q(i,j1).lt.0.) ipy = 1
    1966       enddo
    1967       endif
    1968 C
    1969       if(q(1,JNP).lt.0.) then
     2029      if(q(i,j1)<0.) ipy = 1
     2030      enddo
     2031      endif
     2032C
     2033      if(q(1,JNP)<0.) then
    19702034      dq = q(1,JNP)*cap1/REAL(IMR)*acosp(j2)
    19712035      do i=1,imr
    19722036      q(i,JNP) = 0.
    19732037      q(i,j2) = q(i,j2) + dq
    1974       if(q(i,j2).lt.0.) ipy = 1
     2038      if(q(i,j2)<0.) ipy = 1
    19752039      enddo
    19762040      endif
    … …  
    19882052      ipx = 0
    19892053C Copy & swap direction for vectorization.
    1990       do 25 i=1,imr
    1991       do 25 j=j1,j2
    1992 25    qtmp(j,i) = q(i,j)
    1993 C
    1994       do 55 i=2,imr-1
    1995       do 55 j=j1,j2
    1996       if(qtmp(j,i).lt.0.) then
     2054      do i=1,imr
     2055      do j=j1,j2
     2056      qtmp(j,i) = q(i,j)
     2057      END DO
     2058      END DO
     2059C
     2060      do i=2,imr-1
     2061      do j=j1,j2
     2062      if(qtmp(j,i)<0.) then
    19972063      ipx =  1
    19982064c west
    … …  
    20072073      qtmp(j,i) = qtmp(j,i) + d2 + tiny
    20082074      endif
    2009 55    continue
     2075      END DO
     2076      END DO
    20102077c
    20112078      i=1
    2012       do 65 j=j1,j2
    2013       if(qtmp(j,i).lt.0.) then
     2079      do j=j1,j2
     2080      if(qtmp(j,i)<0.) then
    20142081      ipx =  1
    20152082c west
    … …  
    20252092      qtmp(j,i) = qtmp(j,i) + d2 + tiny
    20262093      endif
    2027 65    continue
     2094      END DO
    20282095      i=IMR
    2029       do 75 j=j1,j2
    2030       if(qtmp(j,i).lt.0.) then
     2096      do j=j1,j2
     2097      if(qtmp(j,i)<0.) then
    20312098      ipx =  1
    20322099c west
    … …  
    20422109      qtmp(j,i) = qtmp(j,i) + d2 + tiny
    20432110      endif
    2044 75    continue
    2045 C
    2046       if(ipx.ne.0) then
    2047       do 85 j=j1,j2
    2048       do 85 i=1,imr
    2049 85    q(i,j) = qtmp(j,i)
     2111      END DO
     2112C
     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
    20502119      else
    20512120C
    20522121C Poles.
    2053       if(q(1,1).lt.0. or. q(1,JNP).lt.0.) ipx = 1
     2122      if(q(1,1)<0 .or. q(1,JNP)<0.) ipx = 1
    20542123      endif
    20552124      return
    … …  
    20652134      integer IC,k,i
    20662135C
    2067       do 4000 IC = 1, nc
    2068 C
    2069       do 1000 k=1,km
    2070       do 1000 i=1,im
     2136      do IC = 1, nc
     2137C
     2138      do k=1,km
     2139      do i=1,im
    20712140      qtmp(i,k) = q(i,km+1-k,IC)
    2072 1000  continue
    2073 C
    2074       do 2000 i=1,im*km
    2075 2000  q(i,1,IC) = qtmp(i,1)
    2076 4000  continue
     2141      END DO
     2142      END DO
     2143C
     2144      do i=1,im*km
     2145      q(i,1,IC) = qtmp(i,1)
     2146      END DO
     2147      END DO
    20772148      return
    20782149      end

  • LMDZ6/trunk/libf/filtrez/inifgn.F

    r2598 r5079  
    2828      pi = 2.* ASIN(1.)
    2929C
    30       DO 5 i=1,iim
     30      DO i=1,iim
    3131       dlonu(i)=  xprimu( i )
    3232       dlonv(i)=  xprimv( i )
    33    5  CONTINUE
     33      END DO
    3434
    35       DO 12 i=1,iim
     35      DO i=1,iim
    3636      sddv(i)   = SQRT(dlonv(i))
    3737      sddu(i)   = SQRT(dlonu(i))
    3838      unsddu(i) = 1./sddu(i)
    3939      unsddv(i) = 1./sddv(i)
    40   12  CONTINUE
     40      END DO
    4141C
    42       DO 17 j=1,iim
    43       DO 17 i=1,iim
     42      DO j=1,iim
     43      DO i=1,iim
    4444      vec(i,j)     = 0.
    4545      vec1(i,j)    = 0.
    4646      eignfnv(i,j) = 0.
    4747      eignfnu(i,j) = 0.
    48   17  CONTINUE
     48      END DO
     49      END DO
    4950c
    5051c
    5152      eignfnv(1,1)    = -1.
    5253      eignfnv(iim,1)  =  1.
    53       DO 20 i=1,imm1
     54      DO i=1,imm1
    5455      eignfnv(i+1,i+1)= -1.
    5556      eignfnv(i,i+1)  =  1.
    56   20  CONTINUE
    57       DO 25 j=1,iim
    58       DO 25 i=1,iim
     57      END DO
     58      DO j=1,iim
     59      DO i=1,iim
    5960      eignfnv(i,j) = eignfnv(i,j)/(sddu(i)*sddv(j))
    60   25  CONTINUE
    61       DO 30 j=1,iim
    62       DO 30 i=1,iim
     61      END DO
     62      END DO
     63      DO j=1,iim
     64      DO i=1,iim
    6365      eignfnu(i,j) = -eignfnv(j,i)
    64   30  CONTINUE
     66      END DO
     67      END DO
    6568c
    6669#ifdef CRAY
Note: See TracChangeset for help on using the changeset viewer.