Ignore:
Timestamp:
Oct 21, 2024, 2:58:45 PM (23 hours ago)
Author:
abarral
Message:

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

File:
1 moved

Legend:

Unmodified
Added
Removed

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

    r5245 r5246  
    22! $Header$
    33!
    4       SUBROUTINE ADVYP(LIMIT,DTY,PBARV,SM,S0,SSX,SY,SZ
    5      .                 ,SSXX,SSXY,SSXZ,SYY,SYZ,SZZ,ntra )
    6       IMPLICIT NONE
    7 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    8 C                                                                 C
    9 C  second-order moments (SOM) advection of tracer in Y direction  C
    10 C                                                                 C
    11 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    12 C                                                                C
    13 C  Source : Pascal Simon ( Meteo, CNRM )                         C
    14 C  Adaptation : A.A. (LGGE)                                      C
    15 C  Derniere Modif : 19/10/95 LAST
    16 C                                                                C
    17 C  sont les arguments d'entree pour le s-pg                      C
    18 C                                                                C
    19 C  argument de sortie du s-pg                                    C
    20 C                                                                C
    21 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    22 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    23 C
    24 C  Rem : Probleme aux poles il faut reecrire ce cas specifique
    25 C        Attention au sens de l'indexation
    26 C
    27 C  parametres principaux du modele
    28 C
    29 C
    30       include "dimensions.h"
    31       include "paramet.h"
    32       include "comgeom.h"
    33  
    34 C  Arguments :
    35 C  ----------
    36 C  dty : frequence fictive d'appel du transport
    37 C  parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
    38 
    39       INTEGER lon,lat,niv
    40       INTEGER i,j,jv,k,kp,l
    41       INTEGER ntra
    42 C      PARAMETER (ntra = 1)
    43 
    44       REAL dty
    45       REAL pbarv ( iip1,jjm, llm )
    46 
    47 C  moments: SM  total mass in each grid box
    48 C           S0  mass of tracer in each grid box
    49 C           Si  1rst order moment in i direction
    50 C
    51       REAL SM(iip1,jjp1,llm)
    52      +    ,S0(iip1,jjp1,llm,ntra)
    53       REAL SSX(iip1,jjp1,llm,ntra)
    54      +    ,SY(iip1,jjp1,llm,ntra)
    55      +    ,SZ(iip1,jjp1,llm,ntra)
    56      +    ,SSXX(iip1,jjp1,llm,ntra)
    57      +    ,SSXY(iip1,jjp1,llm,ntra)
    58      +    ,SSXZ(iip1,jjp1,llm,ntra)
    59      +    ,SYY(iip1,jjp1,llm,ntra)
    60      +    ,SYZ(iip1,jjp1,llm,ntra)
    61      +    ,SZZ(iip1,jjp1,llm,ntra)
    62 C
    63 C  Local :
    64 C  -------
    65 
    66 C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
    67 C  mass fluxes in kg
    68 C  declaration :
    69 
    70       REAL VGRI(iip1,0:jjp1,llm)
    71 
    72 C  Rem : UGRI et WGRI ne sont pas utilises dans
    73 C  cette subroutine ( advection en y uniquement )
    74 C  Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
    75 C
    76 C  the moments F are similarly defined and used as temporary
    77 C  storage for portions of the grid boxes in transit
    78 C
    79 C  the moments Fij are used as temporary storage for
    80 C  portions of the grid boxes in transit at the current level
    81 C
    82 C  work arrays
    83 C
    84 C
    85       REAL F0(iim,0:jjp1,ntra),FM(iim,0:jjp1)
    86       REAL FX(iim,jjm,ntra),FY(iim,jjm,ntra)
    87       REAL FZ(iim,jjm,ntra)
    88       REAL FXX(iim,jjm,ntra),FXY(iim,jjm,ntra)
    89       REAL FXZ(iim,jjm,ntra),FYY(iim,jjm,ntra)
    90       REAL FYZ(iim,jjm,ntra),FZZ(iim,jjm,ntra)
    91       REAL S00(ntra)
    92       REAL SM0             ! Just temporal variable
    93 C
    94 C  work arrays
    95 C
    96       REAL ALF(iim,0:jjp1),ALF1(iim,0:jjp1)
    97       REAL ALFQ(iim,0:jjp1),ALF1Q(iim,0:jjp1)
    98       REAL ALF2(iim,0:jjp1),ALF3(iim,0:jjp1)
    99       REAL ALF4(iim,0:jjp1)
    100       REAL TEMPTM          ! Just temporal variable
    101       REAL SLPMAX,S1MAX,S1NEW,S2NEW
    102 c
    103 C  Special pour poles
    104 c
    105       REAL sbms,sfms,sfzs,sbmn,sfmn,sfzn
    106       REAL sns0(ntra),snsz(ntra),snsm
    107       REAL qy1(iim,llm,ntra),qylat(iim,llm,ntra)
    108       REAL cx1(llm,ntra), cxLAT(llm,ntra)
    109       REAL cy1(llm,ntra), cyLAT(llm,ntra)
    110       REAL z1(iim), zcos(iim), zsin(iim)
    111       REAL SSUM
    112       EXTERNAL SSUM
    113 C
    114       REAL sqi,sqf
    115       LOGICAL LIMIT
    116 
    117       lon = iim         ! rem : Il est possible qu'un pbl. arrive ici
    118       lat = jjp1        ! a cause des dim. differentes entre les
    119       niv = llm         !       tab. S et VGRI
    120                    
    121 c-----------------------------------------------------------------
    122 C initialisations
    123 
    124       sbms = 0.
    125       sfms = 0.
    126       sfzs = 0.
    127       sbmn = 0.
    128       sfmn = 0.
    129       sfzn = 0.
    130 
    131 c-----------------------------------------------------------------
    132 C *** Test : diag de la qtite totale de traceur dans
    133 C            l'atmosphere avant l'advection en Y
    134 c
    135       sqi = 0.
    136       sqf = 0.
    137 
    138       DO l = 1,llm
    139          DO j = 1,jjp1
    140            DO i = 1,iim
    141               sqi = sqi + S0(i,j,l,ntra)
    142            END DO
    143          END DO
    144       END DO
    145       PRINT*,'---------- DIAG DANS ADVY - ENTREE --------'
    146       PRINT*,'sqi=',sqi
    147 
    148 c-----------------------------------------------------------------
    149 C  Interface : adaptation nouveau modele
    150 C  -------------------------------------
    151 C
    152 C  Conversion des flux de masses en kg
    153 C-AA 20/10/94  le signe -1 est necessaire car indexation opposee
    154 
    155       DO 500 l = 1,llm
    156          DO 500 j = 1,jjm
    157             DO 500 i = 1,iip1 
    158             vgri (i,j,llm+1-l)=-1.*pbarv (i,j,l)
    159   500 CONTINUE
    160 
    161 CAA Initialisation de flux fictifs aux bords sup. des boites pol.
    162 
    163       DO l = 1,llm
    164          DO i = 1,iip1 
    165              vgri(i,0,l) = 0.
    166              vgri(i,jjp1,l) = 0.
    167          ENDDO
    168       ENDDO
    169 c
    170 c----------------- START HERE -----------------------
    171 C  boucle sur les niveaux
    172 C
    173       DO 1 L=1,NIV
    174 C
    175 C  place limits on appropriate moments before transport
    176 C      (if flux-limiting is to be applied)
    177 C
    178       IF(.NOT.LIMIT) GO TO 11
    179 C
    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
    184            SLPMAX=AMAX1(S0(I,K,L,JV),0.)
    185            S1MAX=1.5*SLPMAX
    186            S1NEW=AMIN1(S1MAX,AMAX1(-S1MAX,SY(I,K,L,JV)))
    187            S2NEW=AMIN1( 2.*SLPMAX-ABS(S1NEW)/3. ,
    188      +                  AMAX1(ABS(S1NEW)-SLPMAX,SYY(I,K,L,JV)) )
    189            SY (I,K,L,JV)=S1NEW
    190            SYY(I,K,L,JV)=S2NEW
    191        SSXY(I,K,L,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,SSXY(I,K,L,JV)))
    192        SYZ(I,K,L,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,SYZ(I,K,L,JV)))
    193          ELSE
    194            SY (I,K,L,JV)=0.
    195            SYY(I,K,L,JV)=0.
    196            SSXY(I,K,L,JV)=0.
    197            SYZ(I,K,L,JV)=0.
    198          ENDIF
    199  100  CONTINUE
    200  10   CONTINUE
    201 C
     4SUBROUTINE ADVYP(LIMIT,DTY,PBARV,SM,S0,SSX,SY,SZ &
     5        ,SSXX,SSXY,SSXZ,SYY,SYZ,SZZ,ntra )
     6  IMPLICIT NONE
     7  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     8  !                                                                 C
     9  !  second-order moments (SOM) advection of tracer in Y direction  C
     10  !                                                                 C
     11  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     12                                                             ! C
     13  !  Source : Pascal Simon ( Meteo, CNRM )                         C
     14  !  Adaptation : A.A. (LGGE)                                      C
     15  !  Derniere Modif : 19/10/95 LAST
     16                                                             ! C
     17  !  sont les arguments d'entree pour le s-pg                      C
     18  !                                                                C
     19  !  argument de sortie du s-pg                                    C
     20  !                                                                C
     21  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     22  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     23  !
     24  !  Rem : Probleme aux poles il faut reecrire ce cas specifique
     25  !    Attention au sens de l'indexation
     26  !
     27  !  parametres principaux du modele
     28  !
     29  !
     30  include "dimensions.h"
     31  include "paramet.h"
     32  include "comgeom.h"
     33
     34  !  Arguments :
     35  !  ----------
     36  !  dty : frequence fictive d'appel du transport
     37  !  parbu,pbarv : flux de masse en x et y en Pa.m2.s-1
     38
     39  INTEGER :: lon,lat,niv
     40  INTEGER :: i,j,jv,k,kp,l
     41  INTEGER :: ntra
     42   ! PARAMETER (ntra = 1)
     43
     44  REAL :: dty
     45  REAL :: pbarv ( iip1,jjm, llm )
     46
     47  !  moments: SM  total mass in each grid box
     48        ! S0  mass of tracer in each grid box
     49        ! Si  1rst order moment in i direction
     50  !
     51  REAL :: SM(iip1,jjp1,llm) &
     52        ,S0(iip1,jjp1,llm,ntra)
     53  REAL :: SSX(iip1,jjp1,llm,ntra) &
     54        ,SY(iip1,jjp1,llm,ntra) &
     55        ,SZ(iip1,jjp1,llm,ntra) &
     56        ,SSXX(iip1,jjp1,llm,ntra) &
     57        ,SSXY(iip1,jjp1,llm,ntra) &
     58        ,SSXZ(iip1,jjp1,llm,ntra) &
     59        ,SYY(iip1,jjp1,llm,ntra) &
     60        ,SYZ(iip1,jjp1,llm,ntra) &
     61        ,SZZ(iip1,jjp1,llm,ntra)
     62  !
     63  !  Local :
     64  !  -------
     65
     66  !  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
     67  !  mass fluxes in kg
     68  !  declaration :
     69
     70  REAL :: VGRI(iip1,0:jjp1,llm)
     71
     72  !  Rem : UGRI et WGRI ne sont pas utilises dans
     73  !  cette subroutine ( advection en y uniquement )
     74  !  Rem 2 :le dimensionnement de VGRI depend de celui de pbarv
     75  !
     76  !  the moments F are similarly defined and used as temporary
     77  !  storage for portions of the grid boxes in transit
     78  !
     79  !  the moments Fij are used as temporary storage for
     80  !  portions of the grid boxes in transit at the current level
     81  !
     82  !  work arrays
     83  !
     84  !
     85  REAL :: F0(iim,0:jjp1,ntra),FM(iim,0:jjp1)
     86  REAL :: FX(iim,jjm,ntra),FY(iim,jjm,ntra)
     87  REAL :: FZ(iim,jjm,ntra)
     88  REAL :: FXX(iim,jjm,ntra),FXY(iim,jjm,ntra)
     89  REAL :: FXZ(iim,jjm,ntra),FYY(iim,jjm,ntra)
     90  REAL :: FYZ(iim,jjm,ntra),FZZ(iim,jjm,ntra)
     91  REAL :: S00(ntra)
     92  REAL :: SM0             ! Just temporal variable
     93  !
     94  !  work arrays
     95  !
     96  REAL :: ALF(iim,0:jjp1),ALF1(iim,0:jjp1)
     97  REAL :: ALFQ(iim,0:jjp1),ALF1Q(iim,0:jjp1)
     98  REAL :: ALF2(iim,0:jjp1),ALF3(iim,0:jjp1)
     99  REAL :: ALF4(iim,0:jjp1)
     100  REAL :: TEMPTM          ! Just temporal variable
     101  REAL :: SLPMAX,S1MAX,S1NEW,S2NEW
     102  !
     103  !  Special pour poles
     104  !
     105  REAL :: sbms,sfms,sfzs,sbmn,sfmn,sfzn
     106  REAL :: sns0(ntra),snsz(ntra),snsm
     107  REAL :: qy1(iim,llm,ntra),qylat(iim,llm,ntra)
     108  REAL :: cx1(llm,ntra), cxLAT(llm,ntra)
     109  REAL :: cy1(llm,ntra), cyLAT(llm,ntra)
     110  REAL :: z1(iim), zcos(iim), zsin(iim)
     111  REAL :: SSUM
     112  EXTERNAL SSUM
     113  !
     114  REAL :: sqi,sqf
     115  LOGICAL :: LIMIT
     116
     117  lon = iim         ! rem : Il est possible qu'un pbl. arrive ici
     118  lat = jjp1        ! a cause des dim. differentes entre les
     119  niv = llm         !       tab. S et VGRI
     120
     121  !-----------------------------------------------------------------
     122  ! initialisations
     123
     124  sbms = 0.
     125  sfms = 0.
     126  sfzs = 0.
     127  sbmn = 0.
     128  sfmn = 0.
     129  sfzn = 0.
     130
     131  !-----------------------------------------------------------------
     132  ! *** Test : diag de la qtite totale de traceur dans
     133         ! l'atmosphere avant l'advection en Y
     134  !
     135  sqi = 0.
     136  sqf = 0.
     137
     138  DO l = 1,llm
     139     DO j = 1,jjp1
     140       DO i = 1,iim
     141          sqi = sqi + S0(i,j,l,ntra)
     142       END DO
     143     END DO
     144  END DO
     145  PRINT*,'---------- DIAG DANS ADVY - ENTREE --------'
     146  PRINT*,'sqi=',sqi
     147
     148  !-----------------------------------------------------------------
     149  !  Interface : adaptation nouveau modele
     150  !  -------------------------------------
     151  !
     152  !  Conversion des flux de masses en kg
     153  !-AA 20/10/94  le signe -1 est necessaire car indexation opposee
     154
     155  DO l = 1,llm
     156     DO j = 1,jjm
     157        DO i = 1,iip1
     158        vgri (i,j,llm+1-l)=-1.*pbarv (i,j,l)
     159        END DO
     160     END DO
     161  END DO
     162
     163  !AA Initialisation de flux fictifs aux bords sup. des boites pol.
     164
     165  DO l = 1,llm
     166     DO i = 1,iip1
     167         vgri(i,0,l) = 0.
     168         vgri(i,jjp1,l) = 0.
     169     ENDDO
     170  ENDDO
     171  !
     172  !----------------- START HERE -----------------------
     173  !  boucle sur les niveaux
     174  !
     175  DO L=1,NIV
     176  !
     177  !  place limits on appropriate moments before transport
     178  !  (if flux-limiting is to be applied)
     179  !
     180  IF(.NOT.LIMIT) GO TO 11
     181  !
     182  DO JV=1,NTRA
     183  DO K=1,LAT
     184  DO I=1,LON
     185     IF(S0(I,K,L,JV).GT.0.) THEN
     186       SLPMAX=AMAX1(S0(I,K,L,JV),0.)
     187       S1MAX=1.5*SLPMAX
     188       S1NEW=AMIN1(S1MAX,AMAX1(-S1MAX,SY(I,K,L,JV)))
     189       S2NEW=AMIN1( 2.*SLPMAX-ABS(S1NEW)/3. , &
     190             AMAX1(ABS(S1NEW)-SLPMAX,SYY(I,K,L,JV)) )
     191       SY (I,K,L,JV)=S1NEW
     192       SYY(I,K,L,JV)=S2NEW
     193   SSXY(I,K,L,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,SSXY(I,K,L,JV)))
     194   SYZ(I,K,L,JV)=AMIN1(SLPMAX,AMAX1(-SLPMAX,SYZ(I,K,L,JV)))
     195     ELSE
     196       SY (I,K,L,JV)=0.
     197       SYY(I,K,L,JV)=0.
     198       SSXY(I,K,L,JV)=0.
     199       SYZ(I,K,L,JV)=0.
     200     ENDIF
     201  END DO
     202  END DO
     203  END DO
     204  !
    202205 11   CONTINUE
    203 C
    204 C  le flux a travers le pole Nord est traite separement
    205 C
    206       SM0=0.
    207       DO 20 JV=1,NTRA
    208          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
    214            FM(I,0)=-VGRI(I,0,L)*DTY
    215            ALF(I,0)=FM(I,0)/SM(I,1,L)
    216            SM(I,1,L)=SM(I,1,L)-FM(I,0)
    217            SM0=SM0+FM(I,0)
    218          ENDIF
    219 C
    220          ALFQ(I,0)=ALF(I,0)*ALF(I,0)
    221          ALF1(I,0)=1.-ALF(I,0)
    222          ALF1Q(I,0)=ALF1(I,0)*ALF1(I,0)
    223          ALF2(I,0)=ALF1(I,0)-ALF(I,0)
    224          ALF3(I,0)=ALF(I,0)*ALFQ(I,0)
    225          ALF4(I,0)=ALF1(I,0)*ALF1Q(I,0)
    226 C
    227  21   CONTINUE
    228 c    print*,'ADVYP 21'
    229 C
    230       DO 22 JV=1,NTRA
    231       DO 220 I=1,LON
    232 C
    233          IF(VGRI(I,0,L).LE.0.) THEN
    234 C
    235            F0(I,0,JV)=ALF(I,0)* ( S0(I,1,L,JV)-ALF1(I,0)*
    236      +        ( SY(I,1,L,JV)-ALF2(I,0)*SYY(I,1,L,JV) ) )
    237 C
    238            S00(JV)=S00(JV)+F0(I,0,JV)
    239            S0 (I,1,L,JV)=S0(I,1,L,JV)-F0(I,0,JV)
    240            SY (I,1,L,JV)=ALF1Q(I,0)*
    241      +            (SY(I,1,L,JV)+3.*ALF(I,0)*SYY(I,1,L,JV))
    242            SYY(I,1,L,JV)=ALF4 (I,0)*SYY(I,1,L,JV)
    243            SSX (I,1,L,JV)=ALF1 (I,0)*
    244      +            (SSX(I,1,L,JV)+ALF(I,0)*SSXY(I,1,L,JV) )
    245            SZ (I,1,L,JV)=ALF1 (I,0)*
    246      +            (SZ(I,1,L,JV)+ALF(I,0)*SSXZ(I,1,L,JV) )
    247            SSXX(I,1,L,JV)=ALF1 (I,0)*SSXX(I,1,L,JV)
    248            SSXZ(I,1,L,JV)=ALF1 (I,0)*SSXZ(I,1,L,JV)
    249            SZZ(I,1,L,JV)=ALF1 (I,0)*SZZ(I,1,L,JV)
    250            SSXY(I,1,L,JV)=ALF1Q(I,0)*SSXY(I,1,L,JV)
    251            SYZ(I,1,L,JV)=ALF1Q(I,0)*SYZ(I,1,L,JV)
    252 C
    253          ENDIF
    254 C
    255  220  CONTINUE
    256  22   CONTINUE
    257 C
    258       DO 23 I=1,LON
    259          IF(VGRI(I,0,L).GT.0.) THEN
    260            FM(I,0)=VGRI(I,0,L)*DTY
    261            ALF(I,0)=FM(I,0)/SM0
    262          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
    268            F0(I,0,JV)=ALF(I,0)*S00(JV)
    269          ENDIF
    270  240  CONTINUE
    271  24   CONTINUE
    272 C
    273 C  puts the temporary moments Fi into appropriate neighboring boxes
    274 C
    275 c    print*,'av ADVYP 25'
    276       DO 25 I=1,LON
    277 C
    278          IF(VGRI(I,0,L).GT.0.) THEN
    279            SM(I,1,L)=SM(I,1,L)+FM(I,0)
    280            ALF(I,0)=FM(I,0)/SM(I,1,L)
    281          ENDIF
    282 C
    283          ALFQ(I,0)=ALF(I,0)*ALF(I,0)
    284          ALF1(I,0)=1.-ALF(I,0)
    285          ALF1Q(I,0)=ALF1(I,0)*ALF1(I,0)
    286          ALF2(I,0)=ALF1(I,0)-ALF(I,0)
    287          ALF3(I,0)=ALF1(I,0)*ALF(I,0)
    288 C
    289  25   CONTINUE
    290 c    print*,'av ADVYP 25'
    291 C
    292       DO 26 JV=1,NTRA
    293       DO 260 I=1,LON
    294 C
    295          IF(VGRI(I,0,L).GT.0.) THEN
    296 C
    297          TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV)
    298          S0 (I,1,L,JV)=S0(I,1,L,JV)+F0(I,0,JV)
    299          SYY(I,1,L,JV)=ALF1Q(I,0)*SYY(I,1,L,JV)
    300      +        +5.*( ALF3 (I,0)*SY (I,1,L,JV)-ALF2(I,0)*TEMPTM )
    301          SY (I,1,L,JV)=ALF1 (I,0)*SY (I,1,L,JV)+3.*TEMPTM
    302       SSXY(I,1,L,JV)=ALF1 (I,0)*SSXY(I,1,L,JV)+3.*ALF(I,0)*SSX(I,1,L,JV)
    303       SYZ(I,1,L,JV)=ALF1 (I,0)*SYZ(I,1,L,JV)+3.*ALF(I,0)*SZ(I,1,L,JV)
    304 C
    305          ENDIF
    306 C
    307  260  CONTINUE
    308  26   CONTINUE
    309 C
    310 C  calculate flux and moments between adjacent boxes
    311 C  1- create temporary moments/masses for partial boxes in transit
    312 C  2- reajusts moments remaining in the box
    313 C
    314 C  flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
    315 C
    316 c    print*,'av ADVYP 30'
    317       DO 30 K=1,LAT-1
    318       KP=K+1
    319       DO 300 I=1,LON
    320 C
    321          IF(VGRI(I,K,L).LT.0.) THEN
    322            FM(I,K)=-VGRI(I,K,L)*DTY
    323            ALF(I,K)=FM(I,K)/SM(I,KP,L)
    324            SM(I,KP,L)=SM(I,KP,L)-FM(I,K)
    325          ELSE
    326            FM(I,K)=VGRI(I,K,L)*DTY
    327            ALF(I,K)=FM(I,K)/SM(I,K,L)
    328            SM(I,K,L)=SM(I,K,L)-FM(I,K)
    329          ENDIF
    330 C
    331          ALFQ(I,K)=ALF(I,K)*ALF(I,K)
    332          ALF1(I,K)=1.-ALF(I,K)
    333          ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
    334          ALF2(I,K)=ALF1(I,K)-ALF(I,K)
    335          ALF3(I,K)=ALF(I,K)*ALFQ(I,K)
    336          ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K)
    337 C
    338  300  CONTINUE
    339  30   CONTINUE
    340 c    print*,'ap ADVYP 30'
    341 C
    342       DO 31 JV=1,NTRA
    343       DO 31 K=1,LAT-1
    344       KP=K+1
    345       DO 310 I=1,LON
    346 C
    347          IF(VGRI(I,K,L).LT.0.) THEN
    348 C
    349            F0 (I,K,JV)=ALF (I,K)* ( S0(I,KP,L,JV)-ALF1(I,K)*
    350      +        ( SY(I,KP,L,JV)-ALF2(I,K)*SYY(I,KP,L,JV) ) )
    351            FY (I,K,JV)=ALFQ(I,K)*
    352      +                 (SY(I,KP,L,JV)-3.*ALF1(I,K)*SYY(I,KP,L,JV))
    353            FYY(I,K,JV)=ALF3(I,K)*SYY(I,KP,L,JV)
    354            FX (I,K,JV)=ALF (I,K)*
    355      +                 (SSX(I,KP,L,JV)-ALF1(I,K)*SSXY(I,KP,L,JV))
    356            FZ (I,K,JV)=ALF (I,K)*
    357      +                 (SZ(I,KP,L,JV)-ALF1(I,K)*SYZ(I,KP,L,JV))
    358            FXY(I,K,JV)=ALFQ(I,K)*SSXY(I,KP,L,JV)
    359            FYZ(I,K,JV)=ALFQ(I,K)*SYZ(I,KP,L,JV)
    360            FXX(I,K,JV)=ALF (I,K)*SSXX(I,KP,L,JV)
    361            FXZ(I,K,JV)=ALF (I,K)*SSXZ(I,KP,L,JV)
    362            FZZ(I,K,JV)=ALF (I,K)*SZZ(I,KP,L,JV)
    363 C
    364            S0 (I,KP,L,JV)=S0(I,KP,L,JV)-F0(I,K,JV)
    365            SY (I,KP,L,JV)=ALF1Q(I,K)*
    366      +                 (SY(I,KP,L,JV)+3.*ALF(I,K)*SYY(I,KP,L,JV))
    367            SYY(I,KP,L,JV)=ALF4(I,K)*SYY(I,KP,L,JV)
    368            SSX (I,KP,L,JV)=SSX (I,KP,L,JV)-FX (I,K,JV)
    369            SZ (I,KP,L,JV)=SZ (I,KP,L,JV)-FZ (I,K,JV)
    370            SSXX(I,KP,L,JV)=SSXX(I,KP,L,JV)-FXX(I,K,JV)
    371            SSXZ(I,KP,L,JV)=SSXZ(I,KP,L,JV)-FXZ(I,K,JV)
    372            SZZ(I,KP,L,JV)=SZZ(I,KP,L,JV)-FZZ(I,K,JV)
    373            SSXY(I,KP,L,JV)=ALF1Q(I,K)*SSXY(I,KP,L,JV)
    374            SYZ(I,KP,L,JV)=ALF1Q(I,K)*SYZ(I,KP,L,JV)
    375 C
    376          ELSE
    377 C
    378            F0 (I,K,JV)=ALF (I,K)* ( S0(I,K,L,JV)+ALF1(I,K)*
    379      +             ( SY(I,K,L,JV)+ALF2(I,K)*SYY(I,K,L,JV) ) )
    380            FY (I,K,JV)=ALFQ(I,K)*
    381      +                 (SY(I,K,L,JV)+3.*ALF1(I,K)*SYY(I,K,L,JV))
    382            FYY(I,K,JV)=ALF3(I,K)*SYY(I,K,L,JV)
    383       FX (I,K,JV)=ALF (I,K)*(SSX(I,K,L,JV)+ALF1(I,K)*SSXY(I,K,L,JV))
    384       FZ (I,K,JV)=ALF (I,K)*(SZ(I,K,L,JV)+ALF1(I,K)*SYZ(I,K,L,JV))
    385            FXY(I,K,JV)=ALFQ(I,K)*SSXY(I,K,L,JV)
    386            FYZ(I,K,JV)=ALFQ(I,K)*SYZ(I,K,L,JV)
    387            FXX(I,K,JV)=ALF (I,K)*SSXX(I,K,L,JV)
    388            FXZ(I,K,JV)=ALF (I,K)*SSXZ(I,K,L,JV)
    389            FZZ(I,K,JV)=ALF (I,K)*SZZ(I,K,L,JV)
    390 C
    391            S0 (I,K,L,JV)=S0 (I,K,L,JV)-F0 (I,K,JV)
    392            SY (I,K,L,JV)=ALF1Q(I,K)*
    393      +                  (SY(I,K,L,JV)-3.*ALF(I,K)*SYY(I,K,L,JV))
    394            SYY(I,K,L,JV)=ALF4(I,K)*SYY(I,K,L,JV)
    395            SSX (I,K,L,JV)=SSX (I,K,L,JV)-FX (I,K,JV)
    396            SZ (I,K,L,JV)=SZ (I,K,L,JV)-FZ (I,K,JV)
    397            SSXX(I,K,L,JV)=SSXX(I,K,L,JV)-FXX(I,K,JV)
    398            SSXZ(I,K,L,JV)=SSXZ(I,K,L,JV)-FXZ(I,K,JV)
    399            SZZ(I,K,L,JV)=SZZ(I,K,L,JV)-FZZ(I,K,JV)
    400            SSXY(I,K,L,JV)=ALF1Q(I,K)*SSXY(I,K,L,JV)
    401            SYZ(I,K,L,JV)=ALF1Q(I,K)*SYZ(I,K,L,JV)
    402 C
    403          ENDIF
    404 C
    405  310  CONTINUE
    406  31   CONTINUE
    407 c     print*,'ap ADVYP 31'
    408 C
    409 C  puts the temporary moments Fi into appropriate neighboring boxes
    410 C
    411       DO 32 K=1,LAT-1
    412       KP=K+1
    413       DO 320 I=1,LON
    414 C
    415          IF(VGRI(I,K,L).LT.0.) THEN
    416            SM(I,K,L)=SM(I,K,L)+FM(I,K)
    417            ALF(I,K)=FM(I,K)/SM(I,K,L)
    418          ELSE
    419            SM(I,KP,L)=SM(I,KP,L)+FM(I,K)
    420            ALF(I,K)=FM(I,K)/SM(I,KP,L)
    421          ENDIF
    422 C
    423          ALFQ(I,K)=ALF(I,K)*ALF(I,K)
    424          ALF1(I,K)=1.-ALF(I,K)
    425          ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
    426          ALF2(I,K)=ALF1(I,K)-ALF(I,K)
    427          ALF3(I,K)=ALF1(I,K)*ALF(I,K)
    428 C
    429  320  CONTINUE
    430  32   CONTINUE
    431 c     print*,'ap ADVYP 32'
    432 C
    433       DO 33 JV=1,NTRA
    434       DO 33 K=1,LAT-1
    435       KP=K+1
    436       DO 330 I=1,LON
    437 C
    438          IF(VGRI(I,K,L).LT.0.) THEN
    439 C
    440          TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
    441          S0 (I,K,L,JV)=S0(I,K,L,JV)+F0(I,K,JV)
    442        SYY(I,K,L,JV)=ALFQ(I,K)*FYY(I,K,JV)+ALF1Q(I,K)*SYY(I,K,L,JV)
    443      +  +5.*( ALF3(I,K)*(FY(I,K,JV)-SY(I,K,L,JV))+ALF2(I,K)*TEMPTM )
    444          SY (I,K,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*SY(I,K,L,JV)
    445      +            +3.*TEMPTM
    446        SSXY(I,K,L,JV)=ALF (I,K)*FXY(I,K,JV)+ALF1(I,K)*SSXY(I,K,L,JV)
    447      +         +3.*(ALF1(I,K)*FX (I,K,JV)-ALF (I,K)*SSX (I,K,L,JV))
    448        SYZ(I,K,L,JV)=ALF (I,K)*FYZ(I,K,JV)+ALF1(I,K)*SYZ(I,K,L,JV)
    449      +         +3.*(ALF1(I,K)*FZ (I,K,JV)-ALF (I,K)*SZ (I,K,L,JV))
    450          SSX (I,K,L,JV)=SSX (I,K,L,JV)+FX (I,K,JV)
    451          SZ (I,K,L,JV)=SZ (I,K,L,JV)+FZ (I,K,JV)
    452          SSXX(I,K,L,JV)=SSXX(I,K,L,JV)+FXX(I,K,JV)
    453          SSXZ(I,K,L,JV)=SSXZ(I,K,L,JV)+FXZ(I,K,JV)
    454          SZZ(I,K,L,JV)=SZZ(I,K,L,JV)+FZZ(I,K,JV)
    455 C
    456          ELSE
    457 C
    458          TEMPTM=ALF(I,K)*S0(I,KP,L,JV)-ALF1(I,K)*F0(I,K,JV)
    459          S0 (I,KP,L,JV)=S0(I,KP,L,JV)+F0(I,K,JV)
    460        SYY(I,KP,L,JV)=ALFQ(I,K)*FYY(I,K,JV)+ALF1Q(I,K)*SYY(I,KP,L,JV)
    461      +  +5.*( ALF3(I,K)*(SY(I,KP,L,JV)-FY(I,K,JV))-ALF2(I,K)*TEMPTM )
    462          SY (I,KP,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*SY(I,KP,L,JV)
    463      +                 +3.*TEMPTM
    464        SSXY(I,KP,L,JV)=ALF(I,K)*FXY(I,K,JV)+ALF1(I,K)*SSXY(I,KP,L,JV)
    465      +             +3.*(ALF(I,K)*SSX(I,KP,L,JV)-ALF1(I,K)*FX(I,K,JV))
    466          SYZ(I,KP,L,JV)=ALF(I,K)*FYZ(I,K,JV)+ALF1(I,K)*SYZ(I,KP,L,JV)
    467      +             +3.*(ALF(I,K)*SZ(I,KP,L,JV)-ALF1(I,K)*FZ(I,K,JV))
    468          SSX (I,KP,L,JV)=SSX (I,KP,L,JV)+FX (I,K,JV)
    469          SZ (I,KP,L,JV)=SZ (I,KP,L,JV)+FZ (I,K,JV)
    470          SSXX(I,KP,L,JV)=SSXX(I,KP,L,JV)+FXX(I,K,JV)
    471          SSXZ(I,KP,L,JV)=SSXZ(I,KP,L,JV)+FXZ(I,K,JV)
    472          SZZ(I,KP,L,JV)=SZZ(I,KP,L,JV)+FZZ(I,K,JV)
    473 C
    474          ENDIF
    475 C
    476  330  CONTINUE
    477  33   CONTINUE
    478 c     print*,'ap ADVYP 33'
    479 C
    480 C  traitement special pour le pole Sud (idem pole Nord)
    481 C
    482       K=LAT
    483 C
    484       SM0=0.
    485       DO 40 JV=1,NTRA
    486          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            FM(I,K)=VGRI(I,K,L)*DTY
    493            ALF(I,K)=FM(I,K)/SM(I,K,L)
    494            SM(I,K,L)=SM(I,K,L)-FM(I,K)
    495            SM0=SM0+FM(I,K)
    496          ENDIF
    497 C
    498          ALFQ(I,K)=ALF(I,K)*ALF(I,K)
    499          ALF1(I,K)=1.-ALF(I,K)
    500          ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
    501          ALF2(I,K)=ALF1(I,K)-ALF(I,K)
    502          ALF3(I,K)=ALF(I,K)*ALFQ(I,K)
    503          ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K)
    504 C
    505  41   CONTINUE
    506 c     print*,'ap ADVYP 41'
    507 C
    508       DO 42 JV=1,NTRA
    509       DO 420 I=1,LON
    510 C
    511          IF(VGRI(I,K,L).GE.0.) THEN
    512            F0 (I,K,JV)=ALF(I,K)* ( S0(I,K,L,JV)+ALF1(I,K)*
    513      +             ( SY(I,K,L,JV)+ALF2(I,K)*SYY(I,K,L,JV) ) )
    514            S00(JV)=S00(JV)+F0(I,K,JV)
    515 C
    516            S0 (I,K,L,JV)=S0 (I,K,L,JV)-F0 (I,K,JV)
    517            SY (I,K,L,JV)=ALF1Q(I,K)*
    518      +                  (SY(I,K,L,JV)-3.*ALF(I,K)*SYY(I,K,L,JV))
    519            SYY(I,K,L,JV)=ALF4 (I,K)*SYY(I,K,L,JV)
    520       SSX (I,K,L,JV)=ALF1(I,K)*(SSX(I,K,L,JV)-ALF(I,K)*SSXY(I,K,L,JV))
    521       SZ (I,K,L,JV)=ALF1(I,K)*(SZ(I,K,L,JV)-ALF(I,K)*SYZ(I,K,L,JV))
    522            SSXX(I,K,L,JV)=ALF1 (I,K)*SSXX(I,K,L,JV)
    523            SSXZ(I,K,L,JV)=ALF1 (I,K)*SSXZ(I,K,L,JV)
    524            SZZ(I,K,L,JV)=ALF1 (I,K)*SZZ(I,K,L,JV)
    525            SSXY(I,K,L,JV)=ALF1Q(I,K)*SSXY(I,K,L,JV)
    526            SYZ(I,K,L,JV)=ALF1Q(I,K)*SYZ(I,K,L,JV)
    527          ENDIF
    528 C
    529  420  CONTINUE
    530  42   CONTINUE
    531 c     print*,'ap ADVYP 42'
    532 C
    533       DO 43 I=1,LON
    534          IF(VGRI(I,K,L).LT.0.) THEN
    535            FM(I,K)=-VGRI(I,K,L)*DTY
    536            ALF(I,K)=FM(I,K)/SM0
    537          ENDIF
    538  43   CONTINUE
    539 c     print*,'ap ADVYP 43'
    540 C
    541       DO 44 JV=1,NTRA
    542       DO 440 I=1,LON
    543          IF(VGRI(I,K,L).LT.0.) THEN
    544            F0(I,K,JV)=ALF(I,K)*S00(JV)
    545          ENDIF
    546  440  CONTINUE
    547  44   CONTINUE
    548 C
    549 C  puts the temporary moments Fi into appropriate neighboring boxes
    550 C
    551       DO 45 I=1,LON
    552 C
    553          IF(VGRI(I,K,L).LT.0.) THEN
    554            SM(I,K,L)=SM(I,K,L)+FM(I,K)
    555            ALF(I,K)=FM(I,K)/SM(I,K,L)
    556          ENDIF
    557 C
    558          ALFQ(I,K)=ALF(I,K)*ALF(I,K)
    559          ALF1(I,K)=1.-ALF(I,K)
    560          ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
    561          ALF2(I,K)=ALF1(I,K)-ALF(I,K)
    562          ALF3(I,K)=ALF1(I,K)*ALF(I,K)
    563 C
    564  45   CONTINUE
    565 c     print*,'ap ADVYP 45'
    566 C
    567       DO 46 JV=1,NTRA
    568       DO 460 I=1,LON
    569 C
    570          IF(VGRI(I,K,L).LT.0.) THEN
    571 C
    572          TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
    573          S0 (I,K,L,JV)=S0(I,K,L,JV)+F0(I,K,JV)
    574          SYY(I,K,L,JV)=ALF1Q(I,K)*SYY(I,K,L,JV)
    575      +           +5.*(-ALF3 (I,K)*SY (I,K,L,JV)+ALF2(I,K)*TEMPTM )
    576          SY (I,K,L,JV)=ALF1(I,K)*SY (I,K,L,JV)+3.*TEMPTM
    577       SSXY(I,K,L,JV)=ALF1(I,K)*SSXY(I,K,L,JV)-3.*ALF(I,K)*SSX(I,K,L,JV)
    578       SYZ(I,K,L,JV)=ALF1(I,K)*SYZ(I,K,L,JV)-3.*ALF(I,K)*SZ(I,K,L,JV)
    579 C
    580          ENDIF
    581 C
    582  460  CONTINUE
    583  46   CONTINUE
    584 c     print*,'ap ADVYP 46'
    585 C
    586  1    CONTINUE
    587 
    588 c--------------------------------------------------
    589 C     bouclage cyclique horizontal .
    590      
    591       DO l = 1,llm
    592          DO jv = 1,ntra
    593             DO j = 1,jjp1
    594                SM(iip1,j,l) = SM(1,j,l)
    595                S0(iip1,j,l,jv) = S0(1,j,l,jv)
    596                SSX(iip1,j,l,jv) = SSX(1,j,l,jv)   
    597                SY(iip1,j,l,jv) = SY(1,j,l,jv)
    598                SZ(iip1,j,l,jv) = SZ(1,j,l,jv)
    599             END DO
    600          END DO
    601       END DO
    602 
    603 c -------------------------------------------------------------------
    604 C *** Test  negativite:
    605 
    606 c      DO jv = 1,ntra
    607 c       DO l = 1,llm
    608 c         DO j = 1,jjp1
    609 c           DO i = 1,iip1
    610 c              IF (s0( i,j,l,jv ).lt.0.) THEN
    611 c                 PRINT*, '------ S0 < 0 en FIN ADVYP ---'
    612 c                 PRINT*, 'S0(',i,j,l,jv,')=', S0(i,j,l,jv)
    613 cc                 STOP
    614 c              ENDIF
    615 c           ENDDO
    616 c         ENDDO
    617 c       ENDDO
    618 c      ENDDO
    619  
    620    
    621 c -------------------------------------------------------------------
    622 C *** Test : diag de la qtite totale de traceur dans
    623 C            l'atmosphere avant l'advection en Y
    624  
    625        DO l = 1,llm
    626          DO j = 1,jjp1
    627            DO i = 1,iim
    628               sqf = sqf + S0(i,j,l,ntra)
    629            END DO
    630          END DO
     206  !
     207  !  le flux a travers le pole Nord est traite separement
     208  !
     209  SM0=0.
     210  DO JV=1,NTRA
     211     S00(JV)=0.
     212  END DO
     213  !
     214  DO I=1,LON
     215  !
     216     IF(VGRI(I,0,L).LE.0.) THEN
     217       FM(I,0)=-VGRI(I,0,L)*DTY
     218       ALF(I,0)=FM(I,0)/SM(I,1,L)
     219       SM(I,1,L)=SM(I,1,L)-FM(I,0)
     220       SM0=SM0+FM(I,0)
     221     ENDIF
     222  !
     223     ALFQ(I,0)=ALF(I,0)*ALF(I,0)
     224     ALF1(I,0)=1.-ALF(I,0)
     225     ALF1Q(I,0)=ALF1(I,0)*ALF1(I,0)
     226     ALF2(I,0)=ALF1(I,0)-ALF(I,0)
     227     ALF3(I,0)=ALF(I,0)*ALFQ(I,0)
     228     ALF4(I,0)=ALF1(I,0)*ALF1Q(I,0)
     229  !
     230  END DO
     231  ! print*,'ADVYP 21'
     232  !
     233  DO JV=1,NTRA
     234  DO I=1,LON
     235  !
     236     IF(VGRI(I,0,L).LE.0.) THEN
     237  !
     238       F0(I,0,JV)=ALF(I,0)* ( S0(I,1,L,JV)-ALF1(I,0)* &
     239             ( SY(I,1,L,JV)-ALF2(I,0)*SYY(I,1,L,JV) ) )
     240  !
     241       S00(JV)=S00(JV)+F0(I,0,JV)
     242       S0 (I,1,L,JV)=S0(I,1,L,JV)-F0(I,0,JV)
     243       SY (I,1,L,JV)=ALF1Q(I,0)* &
     244             (SY(I,1,L,JV)+3.*ALF(I,0)*SYY(I,1,L,JV))
     245       SYY(I,1,L,JV)=ALF4 (I,0)*SYY(I,1,L,JV)
     246       SSX (I,1,L,JV)=ALF1 (I,0)* &
     247             (SSX(I,1,L,JV)+ALF(I,0)*SSXY(I,1,L,JV) )
     248       SZ (I,1,L,JV)=ALF1 (I,0)* &
     249             (SZ(I,1,L,JV)+ALF(I,0)*SSXZ(I,1,L,JV) )
     250       SSXX(I,1,L,JV)=ALF1 (I,0)*SSXX(I,1,L,JV)
     251       SSXZ(I,1,L,JV)=ALF1 (I,0)*SSXZ(I,1,L,JV)
     252       SZZ(I,1,L,JV)=ALF1 (I,0)*SZZ(I,1,L,JV)
     253       SSXY(I,1,L,JV)=ALF1Q(I,0)*SSXY(I,1,L,JV)
     254       SYZ(I,1,L,JV)=ALF1Q(I,0)*SYZ(I,1,L,JV)
     255  !
     256     ENDIF
     257  !
     258  END DO
     259  END DO
     260  !
     261  DO I=1,LON
     262     IF(VGRI(I,0,L).GT.0.) THEN
     263       FM(I,0)=VGRI(I,0,L)*DTY
     264       ALF(I,0)=FM(I,0)/SM0
     265     ENDIF
     266  END DO
     267  !
     268  DO JV=1,NTRA
     269  DO I=1,LON
     270     IF(VGRI(I,0,L).GT.0.) THEN
     271       F0(I,0,JV)=ALF(I,0)*S00(JV)
     272     ENDIF
     273  END DO
     274  END DO
     275  !
     276  !  puts the temporary moments Fi into appropriate neighboring boxes
     277  !
     278  ! print*,'av ADVYP 25'
     279  DO I=1,LON
     280  !
     281     IF(VGRI(I,0,L).GT.0.) THEN
     282       SM(I,1,L)=SM(I,1,L)+FM(I,0)
     283       ALF(I,0)=FM(I,0)/SM(I,1,L)
     284     ENDIF
     285  !
     286     ALFQ(I,0)=ALF(I,0)*ALF(I,0)
     287     ALF1(I,0)=1.-ALF(I,0)
     288     ALF1Q(I,0)=ALF1(I,0)*ALF1(I,0)
     289     ALF2(I,0)=ALF1(I,0)-ALF(I,0)
     290     ALF3(I,0)=ALF1(I,0)*ALF(I,0)
     291  !
     292  END DO
     293  ! print*,'av ADVYP 25'
     294  !
     295  DO JV=1,NTRA
     296  DO I=1,LON
     297  !
     298     IF(VGRI(I,0,L).GT.0.) THEN
     299  !
     300     TEMPTM=ALF(I,0)*S0(I,1,L,JV)-ALF1(I,0)*F0(I,0,JV)
     301     S0 (I,1,L,JV)=S0(I,1,L,JV)+F0(I,0,JV)
     302     SYY(I,1,L,JV)=ALF1Q(I,0)*SYY(I,1,L,JV) &
     303           +5.*( ALF3 (I,0)*SY (I,1,L,JV)-ALF2(I,0)*TEMPTM )
     304     SY (I,1,L,JV)=ALF1 (I,0)*SY (I,1,L,JV)+3.*TEMPTM
     305  SSXY(I,1,L,JV)=ALF1 (I,0)*SSXY(I,1,L,JV)+3.*ALF(I,0)*SSX(I,1,L,JV)
     306  SYZ(I,1,L,JV)=ALF1 (I,0)*SYZ(I,1,L,JV)+3.*ALF(I,0)*SZ(I,1,L,JV)
     307  !
     308     ENDIF
     309  !
     310  END DO
     311  END DO
     312  !
     313  !  calculate flux and moments between adjacent boxes
     314  !  1- create temporary moments/masses for partial boxes in transit
     315  !  2- reajusts moments remaining in the box
     316  !
     317  !  flux from KP to K if V(K).lt.0 and from K to KP if V(K).gt.0
     318  !
     319  ! print*,'av ADVYP 30'
     320  DO K=1,LAT-1
     321  KP=K+1
     322  DO I=1,LON
     323  !
     324     IF(VGRI(I,K,L).LT.0.) THEN
     325       FM(I,K)=-VGRI(I,K,L)*DTY
     326       ALF(I,K)=FM(I,K)/SM(I,KP,L)
     327       SM(I,KP,L)=SM(I,KP,L)-FM(I,K)
     328     ELSE
     329       FM(I,K)=VGRI(I,K,L)*DTY
     330       ALF(I,K)=FM(I,K)/SM(I,K,L)
     331       SM(I,K,L)=SM(I,K,L)-FM(I,K)
     332     ENDIF
     333  !
     334     ALFQ(I,K)=ALF(I,K)*ALF(I,K)
     335     ALF1(I,K)=1.-ALF(I,K)
     336     ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
     337     ALF2(I,K)=ALF1(I,K)-ALF(I,K)
     338     ALF3(I,K)=ALF(I,K)*ALFQ(I,K)
     339     ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K)
     340  !
     341  END DO
     342  END DO
     343  ! print*,'ap ADVYP 30'
     344  !
     345  DO JV=1,NTRA
     346  DO K=1,LAT-1
     347  KP=K+1
     348  DO I=1,LON
     349  !
     350     IF(VGRI(I,K,L).LT.0.) THEN
     351  !
     352       F0 (I,K,JV)=ALF (I,K)* ( S0(I,KP,L,JV)-ALF1(I,K)* &
     353             ( SY(I,KP,L,JV)-ALF2(I,K)*SYY(I,KP,L,JV) ) )
     354       FY (I,K,JV)=ALFQ(I,K)* &
     355             (SY(I,KP,L,JV)-3.*ALF1(I,K)*SYY(I,KP,L,JV))
     356       FYY(I,K,JV)=ALF3(I,K)*SYY(I,KP,L,JV)
     357       FX (I,K,JV)=ALF (I,K)* &
     358             (SSX(I,KP,L,JV)-ALF1(I,K)*SSXY(I,KP,L,JV))
     359       FZ (I,K,JV)=ALF (I,K)* &
     360             (SZ(I,KP,L,JV)-ALF1(I,K)*SYZ(I,KP,L,JV))
     361       FXY(I,K,JV)=ALFQ(I,K)*SSXY(I,KP,L,JV)
     362       FYZ(I,K,JV)=ALFQ(I,K)*SYZ(I,KP,L,JV)
     363       FXX(I,K,JV)=ALF (I,K)*SSXX(I,KP,L,JV)
     364       FXZ(I,K,JV)=ALF (I,K)*SSXZ(I,KP,L,JV)
     365       FZZ(I,K,JV)=ALF (I,K)*SZZ(I,KP,L,JV)
     366  !
     367       S0 (I,KP,L,JV)=S0(I,KP,L,JV)-F0(I,K,JV)
     368       SY (I,KP,L,JV)=ALF1Q(I,K)* &
     369             (SY(I,KP,L,JV)+3.*ALF(I,K)*SYY(I,KP,L,JV))
     370       SYY(I,KP,L,JV)=ALF4(I,K)*SYY(I,KP,L,JV)
     371       SSX (I,KP,L,JV)=SSX (I,KP,L,JV)-FX (I,K,JV)
     372       SZ (I,KP,L,JV)=SZ (I,KP,L,JV)-FZ (I,K,JV)
     373       SSXX(I,KP,L,JV)=SSXX(I,KP,L,JV)-FXX(I,K,JV)
     374       SSXZ(I,KP,L,JV)=SSXZ(I,KP,L,JV)-FXZ(I,K,JV)
     375       SZZ(I,KP,L,JV)=SZZ(I,KP,L,JV)-FZZ(I,K,JV)
     376       SSXY(I,KP,L,JV)=ALF1Q(I,K)*SSXY(I,KP,L,JV)
     377       SYZ(I,KP,L,JV)=ALF1Q(I,K)*SYZ(I,KP,L,JV)
     378  !
     379     ELSE
     380  !
     381       F0 (I,K,JV)=ALF (I,K)* ( S0(I,K,L,JV)+ALF1(I,K)* &
     382             ( SY(I,K,L,JV)+ALF2(I,K)*SYY(I,K,L,JV) ) )
     383       FY (I,K,JV)=ALFQ(I,K)* &
     384             (SY(I,K,L,JV)+3.*ALF1(I,K)*SYY(I,K,L,JV))
     385       FYY(I,K,JV)=ALF3(I,K)*SYY(I,K,L,JV)
     386  FX (I,K,JV)=ALF (I,K)*(SSX(I,K,L,JV)+ALF1(I,K)*SSXY(I,K,L,JV))
     387  FZ (I,K,JV)=ALF (I,K)*(SZ(I,K,L,JV)+ALF1(I,K)*SYZ(I,K,L,JV))
     388       FXY(I,K,JV)=ALFQ(I,K)*SSXY(I,K,L,JV)
     389       FYZ(I,K,JV)=ALFQ(I,K)*SYZ(I,K,L,JV)
     390       FXX(I,K,JV)=ALF (I,K)*SSXX(I,K,L,JV)
     391       FXZ(I,K,JV)=ALF (I,K)*SSXZ(I,K,L,JV)
     392       FZZ(I,K,JV)=ALF (I,K)*SZZ(I,K,L,JV)
     393  !
     394       S0 (I,K,L,JV)=S0 (I,K,L,JV)-F0 (I,K,JV)
     395       SY (I,K,L,JV)=ALF1Q(I,K)* &
     396             (SY(I,K,L,JV)-3.*ALF(I,K)*SYY(I,K,L,JV))
     397       SYY(I,K,L,JV)=ALF4(I,K)*SYY(I,K,L,JV)
     398       SSX (I,K,L,JV)=SSX (I,K,L,JV)-FX (I,K,JV)
     399       SZ (I,K,L,JV)=SZ (I,K,L,JV)-FZ (I,K,JV)
     400       SSXX(I,K,L,JV)=SSXX(I,K,L,JV)-FXX(I,K,JV)
     401       SSXZ(I,K,L,JV)=SSXZ(I,K,L,JV)-FXZ(I,K,JV)
     402       SZZ(I,K,L,JV)=SZZ(I,K,L,JV)-FZZ(I,K,JV)
     403       SSXY(I,K,L,JV)=ALF1Q(I,K)*SSXY(I,K,L,JV)
     404       SYZ(I,K,L,JV)=ALF1Q(I,K)*SYZ(I,K,L,JV)
     405  !
     406     ENDIF
     407  !
     408  END DO
     409  END DO
     410  END DO
     411  ! print*,'ap ADVYP 31'
     412  !
     413  !  puts the temporary moments Fi into appropriate neighboring boxes
     414  !
     415  DO K=1,LAT-1
     416  KP=K+1
     417  DO I=1,LON
     418  !
     419     IF(VGRI(I,K,L).LT.0.) THEN
     420       SM(I,K,L)=SM(I,K,L)+FM(I,K)
     421       ALF(I,K)=FM(I,K)/SM(I,K,L)
     422     ELSE
     423       SM(I,KP,L)=SM(I,KP,L)+FM(I,K)
     424       ALF(I,K)=FM(I,K)/SM(I,KP,L)
     425     ENDIF
     426  !
     427     ALFQ(I,K)=ALF(I,K)*ALF(I,K)
     428     ALF1(I,K)=1.-ALF(I,K)
     429     ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
     430     ALF2(I,K)=ALF1(I,K)-ALF(I,K)
     431     ALF3(I,K)=ALF1(I,K)*ALF(I,K)
     432  !
     433  END DO
     434  END DO
     435  ! print*,'ap ADVYP 32'
     436  !
     437  DO JV=1,NTRA
     438  DO K=1,LAT-1
     439  KP=K+1
     440  DO I=1,LON
     441  !
     442     IF(VGRI(I,K,L).LT.0.) THEN
     443  !
     444     TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
     445     S0 (I,K,L,JV)=S0(I,K,L,JV)+F0(I,K,JV)
     446   SYY(I,K,L,JV)=ALFQ(I,K)*FYY(I,K,JV)+ALF1Q(I,K)*SYY(I,K,L,JV) &
     447         +5.*( ALF3(I,K)*(FY(I,K,JV)-SY(I,K,L,JV))+ALF2(I,K)*TEMPTM )
     448     SY (I,K,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*SY(I,K,L,JV) &
     449           +3.*TEMPTM
     450   SSXY(I,K,L,JV)=ALF (I,K)*FXY(I,K,JV)+ALF1(I,K)*SSXY(I,K,L,JV) &
     451         +3.*(ALF1(I,K)*FX (I,K,JV)-ALF (I,K)*SSX (I,K,L,JV))
     452   SYZ(I,K,L,JV)=ALF (I,K)*FYZ(I,K,JV)+ALF1(I,K)*SYZ(I,K,L,JV) &
     453         +3.*(ALF1(I,K)*FZ (I,K,JV)-ALF (I,K)*SZ (I,K,L,JV))
     454     SSX (I,K,L,JV)=SSX (I,K,L,JV)+FX (I,K,JV)
     455     SZ (I,K,L,JV)=SZ (I,K,L,JV)+FZ (I,K,JV)
     456     SSXX(I,K,L,JV)=SSXX(I,K,L,JV)+FXX(I,K,JV)
     457     SSXZ(I,K,L,JV)=SSXZ(I,K,L,JV)+FXZ(I,K,JV)
     458     SZZ(I,K,L,JV)=SZZ(I,K,L,JV)+FZZ(I,K,JV)
     459  !
     460     ELSE
     461  !
     462     TEMPTM=ALF(I,K)*S0(I,KP,L,JV)-ALF1(I,K)*F0(I,K,JV)
     463     S0 (I,KP,L,JV)=S0(I,KP,L,JV)+F0(I,K,JV)
     464   SYY(I,KP,L,JV)=ALFQ(I,K)*FYY(I,K,JV)+ALF1Q(I,K)*SYY(I,KP,L,JV) &
     465         +5.*( ALF3(I,K)*(SY(I,KP,L,JV)-FY(I,K,JV))-ALF2(I,K)*TEMPTM )
     466     SY (I,KP,L,JV)=ALF(I,K)*FY(I,K,JV)+ALF1(I,K)*SY(I,KP,L,JV) &
     467           +3.*TEMPTM
     468   SSXY(I,KP,L,JV)=ALF(I,K)*FXY(I,K,JV)+ALF1(I,K)*SSXY(I,KP,L,JV) &
     469         +3.*(ALF(I,K)*SSX(I,KP,L,JV)-ALF1(I,K)*FX(I,K,JV))
     470     SYZ(I,KP,L,JV)=ALF(I,K)*FYZ(I,K,JV)+ALF1(I,K)*SYZ(I,KP,L,JV) &
     471           +3.*(ALF(I,K)*SZ(I,KP,L,JV)-ALF1(I,K)*FZ(I,K,JV))
     472     SSX (I,KP,L,JV)=SSX (I,KP,L,JV)+FX (I,K,JV)
     473     SZ (I,KP,L,JV)=SZ (I,KP,L,JV)+FZ (I,K,JV)
     474     SSXX(I,KP,L,JV)=SSXX(I,KP,L,JV)+FXX(I,K,JV)
     475     SSXZ(I,KP,L,JV)=SSXZ(I,KP,L,JV)+FXZ(I,K,JV)
     476     SZZ(I,KP,L,JV)=SZZ(I,KP,L,JV)+FZZ(I,K,JV)
     477  !
     478     ENDIF
     479  !
     480  END DO
     481  END DO
     482  END DO
     483  ! print*,'ap ADVYP 33'
     484  !
     485  !  traitement special pour le pole Sud (idem pole Nord)
     486  !
     487  K=LAT
     488  !
     489  SM0=0.
     490  DO JV=1,NTRA
     491     S00(JV)=0.
     492  END DO
     493  !
     494  DO I=1,LON
     495  !
     496     IF(VGRI(I,K,L).GE.0.) THEN
     497       FM(I,K)=VGRI(I,K,L)*DTY
     498       ALF(I,K)=FM(I,K)/SM(I,K,L)
     499       SM(I,K,L)=SM(I,K,L)-FM(I,K)
     500       SM0=SM0+FM(I,K)
     501     ENDIF
     502  !
     503     ALFQ(I,K)=ALF(I,K)*ALF(I,K)
     504     ALF1(I,K)=1.-ALF(I,K)
     505     ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
     506     ALF2(I,K)=ALF1(I,K)-ALF(I,K)
     507     ALF3(I,K)=ALF(I,K)*ALFQ(I,K)
     508     ALF4(I,K)=ALF1(I,K)*ALF1Q(I,K)
     509  !
     510  END DO
     511  ! print*,'ap ADVYP 41'
     512  !
     513  DO JV=1,NTRA
     514  DO I=1,LON
     515  !
     516     IF(VGRI(I,K,L).GE.0.) THEN
     517       F0 (I,K,JV)=ALF(I,K)* ( S0(I,K,L,JV)+ALF1(I,K)* &
     518             ( SY(I,K,L,JV)+ALF2(I,K)*SYY(I,K,L,JV) ) )
     519       S00(JV)=S00(JV)+F0(I,K,JV)
     520  !
     521       S0 (I,K,L,JV)=S0 (I,K,L,JV)-F0 (I,K,JV)
     522       SY (I,K,L,JV)=ALF1Q(I,K)* &
     523             (SY(I,K,L,JV)-3.*ALF(I,K)*SYY(I,K,L,JV))
     524       SYY(I,K,L,JV)=ALF4 (I,K)*SYY(I,K,L,JV)
     525  SSX (I,K,L,JV)=ALF1(I,K)*(SSX(I,K,L,JV)-ALF(I,K)*SSXY(I,K,L,JV))
     526  SZ (I,K,L,JV)=ALF1(I,K)*(SZ(I,K,L,JV)-ALF(I,K)*SYZ(I,K,L,JV))
     527       SSXX(I,K,L,JV)=ALF1 (I,K)*SSXX(I,K,L,JV)
     528       SSXZ(I,K,L,JV)=ALF1 (I,K)*SSXZ(I,K,L,JV)
     529       SZZ(I,K,L,JV)=ALF1 (I,K)*SZZ(I,K,L,JV)
     530       SSXY(I,K,L,JV)=ALF1Q(I,K)*SSXY(I,K,L,JV)
     531       SYZ(I,K,L,JV)=ALF1Q(I,K)*SYZ(I,K,L,JV)
     532     ENDIF
     533  !
     534  END DO
     535  END DO
     536  ! print*,'ap ADVYP 42'
     537  !
     538  DO I=1,LON
     539     IF(VGRI(I,K,L).LT.0.) THEN
     540       FM(I,K)=-VGRI(I,K,L)*DTY
     541       ALF(I,K)=FM(I,K)/SM0
     542     ENDIF
     543  END DO
     544  ! print*,'ap ADVYP 43'
     545  !
     546  DO JV=1,NTRA
     547  DO I=1,LON
     548     IF(VGRI(I,K,L).LT.0.) THEN
     549       F0(I,K,JV)=ALF(I,K)*S00(JV)
     550     ENDIF
     551  END DO
     552  END DO
     553  !
     554  !  puts the temporary moments Fi into appropriate neighboring boxes
     555  !
     556  DO I=1,LON
     557  !
     558     IF(VGRI(I,K,L).LT.0.) THEN
     559       SM(I,K,L)=SM(I,K,L)+FM(I,K)
     560       ALF(I,K)=FM(I,K)/SM(I,K,L)
     561     ENDIF
     562  !
     563     ALFQ(I,K)=ALF(I,K)*ALF(I,K)
     564     ALF1(I,K)=1.-ALF(I,K)
     565     ALF1Q(I,K)=ALF1(I,K)*ALF1(I,K)
     566     ALF2(I,K)=ALF1(I,K)-ALF(I,K)
     567     ALF3(I,K)=ALF1(I,K)*ALF(I,K)
     568  !
     569  END DO
     570  ! print*,'ap ADVYP 45'
     571  !
     572  DO JV=1,NTRA
     573  DO I=1,LON
     574  !
     575     IF(VGRI(I,K,L).LT.0.) THEN
     576  !
     577     TEMPTM=-ALF(I,K)*S0(I,K,L,JV)+ALF1(I,K)*F0(I,K,JV)
     578     S0 (I,K,L,JV)=S0(I,K,L,JV)+F0(I,K,JV)
     579     SYY(I,K,L,JV)=ALF1Q(I,K)*SYY(I,K,L,JV) &
     580           +5.*(-ALF3 (I,K)*SY (I,K,L,JV)+ALF2(I,K)*TEMPTM )
     581     SY (I,K,L,JV)=ALF1(I,K)*SY (I,K,L,JV)+3.*TEMPTM
     582  SSXY(I,K,L,JV)=ALF1(I,K)*SSXY(I,K,L,JV)-3.*ALF(I,K)*SSX(I,K,L,JV)
     583  SYZ(I,K,L,JV)=ALF1(I,K)*SYZ(I,K,L,JV)-3.*ALF(I,K)*SZ(I,K,L,JV)
     584  !
     585     ENDIF
     586  !
     587  END DO
     588  END DO
     589  ! print*,'ap ADVYP 46'
     590  !
     591  END DO
     592
     593  !--------------------------------------------------
     594  ! bouclage cyclique horizontal .
     595
     596  DO l = 1,llm
     597     DO jv = 1,ntra
     598        DO j = 1,jjp1
     599           SM(iip1,j,l) = SM(1,j,l)
     600           S0(iip1,j,l,jv) = S0(1,j,l,jv)
     601           SSX(iip1,j,l,jv) = SSX(1,j,l,jv)
     602           SY(iip1,j,l,jv) = SY(1,j,l,jv)
     603           SZ(iip1,j,l,jv) = SZ(1,j,l,jv)
     604        END DO
     605     END DO
     606  END DO
     607
     608  ! -------------------------------------------------------------------
     609  ! *** Test  negativite:
     610
     611   ! DO jv = 1,ntra
     612   !  DO l = 1,llm
     613   !    DO j = 1,jjp1
     614   !      DO i = 1,iip1
     615   !         IF (s0( i,j,l,jv ).lt.0.) THEN
     616   !            PRINT*, '------ S0 < 0 en FIN ADVYP ---'
     617   !            PRINT*, 'S0(',i,j,l,jv,')=', S0(i,j,l,jv)
     618  !c                 STOP
     619   !         ENDIF
     620   !      ENDDO
     621   !    ENDDO
     622   !  ENDDO
     623   ! ENDDO
     624
     625
     626  ! -------------------------------------------------------------------
     627  ! *** Test : diag de la qtite totale de traceur dans
     628   !       l'atmosphere avant l'advection en Y
     629
     630   DO l = 1,llm
     631     DO j = 1,jjp1
     632       DO i = 1,iim
     633          sqf = sqf + S0(i,j,l,ntra)
    631634       END DO
    632       PRINT*,'---------- DIAG DANS ADVY - SORTIE --------'
    633       PRINT*,'sqf=',sqf
    634 c     print*,'ap ADVYP fin'
    635 
    636 c-----------------------------------------------------------------
    637 C
    638       RETURN
    639       END
    640 
    641 
    642 
    643 
    644 
    645 
    646 
    647 
    648 
    649 
    650 
    651 
     635     END DO
     636   END DO
     637  PRINT*,'---------- DIAG DANS ADVY - SORTIE --------'
     638  PRINT*,'sqf=',sqf
     639  ! print*,'ap ADVYP fin'
     640
     641  !-----------------------------------------------------------------
     642  !
     643  RETURN
     644END SUBROUTINE ADVYP
     645
     646
     647
     648
     649
     650
     651
     652
     653
     654
     655
     656
Note: See TracChangeset for help on using the changeset viewer.