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

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

File:
1 moved

Legend:

Unmodified
Added
Removed

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

    r5104 r5105  
    22! $Header$
    33
    4       SUBROUTINE  advx(limit,dtx,pbaru,sm,s0,
    5      $     sx,sy,sz,lati,latf)
    6       IMPLICIT NONE
    7 
    8 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    9 C                                                                C
    10 C  first-order moments (FOM) advection of tracer in X direction  C
    11 C                                                                C
    12 C  Source : Pascal Simon (Meteo,CNRM)                            C
    13 C  Adaptation : A.Armengaud (LGGE) juin 94                       C
    14 C                                                                C
    15 C  limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz                       C
    16 C  sont des arguments d'entree pour le s-pg...                   C
    17 C                                                                C
    18 C  sm,s0,sx,sy,sz                                                C
    19 C  sont les arguments de sortie pour le s-pg                     C
    20 C                                                                C
    21 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    22 C
    23 C  parametres principaux du modele
    24 C
    25       include "dimensions.h"
    26       include "paramet.h"
    27 
    28 C  Arguments :
    29 C  -----------
    30 C  dtx : frequence fictive d'appel du transport
    31 C  pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1
    32 
    33        INTEGER ntra
    34        PARAMETER (ntra = 1)
    35 
    36 C ATTENTION partout ou on trouve ntra, insertion de boucle
    37 C           possible dans l'avenir.
    38 
    39       REAL dtx
    40       REAL pbaru ( iip1,jjp1,llm )
    41 
    42 C  moments: SM  total mass in each grid box
    43 C           S0  mass of tracer in each grid box
    44 C           Si  1rst order moment in i direction
    45 C
    46       REAL SM(iip1,jjp1,llm),S0(iip1,jjp1,llm,ntra)
    47       REAL sx(iip1,jjp1,llm,ntra)
    48      $    ,sy(iip1,jjp1,llm,ntra)
    49       REAL sz(iip1,jjp1,llm,ntra)
    50 
    51 C  Local :
    52 C  -------
    53 
    54 C  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
    55 C  mass fluxes in kg
    56 C  declaration :
    57 
    58       REAL UGRI(iip1,jjp1,llm)
    59 
    60 C  Rem : VGRI et WGRI ne sont pas utilises dans
    61 C  cette SUBROUTINE ( advection en x uniquement )
    62 C
    63 C  Ti are the moments for the current latitude and level
    64 C
    65       REAL TM(iim)
    66       REAL T0(iim,ntra),TX(iim,ntra)
    67       REAL TY(iim,ntra),TZ(iim,ntra)
    68       REAL TEMPTM                ! just a temporary variable
    69 C
    70 C  the moments F are similarly defined and used as temporary
    71 C  storage for portions of the grid boxes in transit
    72 C
    73       REAL FM(iim)
    74       REAL F0(iim,ntra),FX(iim,ntra)
    75       REAL FY(iim,ntra),FZ(iim,ntra)
    76 C
    77 C  work arrays
    78 C
    79       REAL ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
    80 C
    81       REAL SMNEW(iim),UEXT(iim)
    82 C
    83       REAL sqi,sqf
    84 
    85       LOGICAL LIMIT
    86       INTEGER NUM(jjp1),LONK,NUMK
    87       INTEGER lon,lati,latf,niv
    88       INTEGER i,i2,i3,j,jv,l,k,itrac
    89 
    90       lon = iim
    91       niv = llm
    92 
    93 C *** Test de passage d'arguments ******
    94 
    95 
    96 C  -------------------------------------
    97       DO j = 1,jjp1
    98          NUM(j) = 1
     4SUBROUTINE  advx(limit,dtx,pbaru,sm,s0, &
     5        sx,sy,sz,lati,latf)
     6  IMPLICIT NONE
     7
     8  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     9  !                                                                C
     10  !  first-order moments (FOM) advection of tracer in X direction  C
     11  !                                                                C
     12  !  Source : Pascal Simon (Meteo,CNRM)                            C
     13  !  Adaptation : A.Armengaud (LGGE) juin 94                       C
     14  !                                                                C
     15  !  limit,dtx,pbaru,pbarv,sm,s0,sx,sy,sz                       C
     16  !  sont des arguments d'entree pour le s-pg...                   C
     17  !                                                                C
     18  !  sm,s0,sx,sy,sz                                                C
     19  !  sont les arguments de sortie pour le s-pg                     C
     20  !                                                                C
     21  !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
     22  !
     23  !  parametres principaux du modele
     24  !
     25  include "dimensions.h"
     26  include "paramet.h"
     27
     28  !  Arguments :
     29  !  -----------
     30  !  dtx : frequence fictive d'appel du transport
     31  !  pbaru, pbarv : flux de masse en x et y en Pa.m2.s-1
     32
     33   INTEGER :: ntra
     34   PARAMETER (ntra = 1)
     35
     36  ! ATTENTION partout ou on trouve ntra, insertion de boucle
     37        ! possible dans l'avenir.
     38
     39  REAL :: dtx
     40  REAL :: pbaru ( iip1,jjp1,llm )
     41
     42  !  moments: SM  total mass in each grid box
     43        ! S0  mass of tracer in each grid box
     44        ! Si  1rst order moment in i direction
     45  !
     46  REAL :: SM(iip1,jjp1,llm),S0(iip1,jjp1,llm,ntra)
     47  REAL :: sx(iip1,jjp1,llm,ntra) &
     48        ,sy(iip1,jjp1,llm,ntra)
     49  REAL :: sz(iip1,jjp1,llm,ntra)
     50
     51  !  Local :
     52  !  -------
     53
     54  !  mass fluxes across the boundaries (UGRI,VGRI,WGRI)
     55  !  mass fluxes in kg
     56  !  declaration :
     57
     58  REAL :: UGRI(iip1,jjp1,llm)
     59
     60  !  Rem : VGRI et WGRI ne sont pas utilises dans
     61  !  cette SUBROUTINE ( advection en x uniquement )
     62  !
     63  !  Ti are the moments for the current latitude and level
     64  !
     65  REAL :: TM(iim)
     66  REAL :: T0(iim,ntra),TX(iim,ntra)
     67  REAL :: TY(iim,ntra),TZ(iim,ntra)
     68  REAL :: TEMPTM                ! just a temporary variable
     69  !
     70  !  the moments F are similarly defined and used as temporary
     71  !  storage for portions of the grid boxes in transit
     72  !
     73  REAL :: FM(iim)
     74  REAL :: F0(iim,ntra),FX(iim,ntra)
     75  REAL :: FY(iim,ntra),FZ(iim,ntra)
     76  !
     77  !  work arrays
     78  !
     79  REAL :: ALF(iim),ALF1(iim),ALFQ(iim),ALF1Q(iim)
     80  !
     81  REAL :: SMNEW(iim),UEXT(iim)
     82  !
     83  REAL :: sqi,sqf
     84
     85  LOGICAL :: LIMIT
     86  INTEGER :: NUM(jjp1),LONK,NUMK
     87  INTEGER :: lon,lati,latf,niv
     88  INTEGER :: i,i2,i3,j,jv,l,k,itrac
     89
     90  lon = iim
     91  niv = llm
     92
     93  ! *** Test de passage d'arguments ******
     94
     95
     96  !  -------------------------------------
     97  DO j = 1,jjp1
     98     NUM(j) = 1
     99  END DO
     100  sqi = 0.
     101  sqf = 0.
     102
     103  DO l = 1,llm
     104     DO j = 1,jjp1
     105        DO i = 1,iim
     106  !IM 240305            sqi = sqi + S0(i,j,l,9)
     107           sqi = sqi + S0(i,j,l,ntra)
     108        ENDDO
     109     ENDDO
     110  ENDDO
     111  PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------'
     112  PRINT*,'sqi=',sqi
     113
     114
     115  !  Interface : adaptation nouveau modele
     116  !  -------------------------------------
     117  !
     118  !  ---------------------------------------------------------
     119  !  Conversion des flux de masses en kg/s
     120  !  pbaru est en N/s d'ou :
     121  !  ugri est en kg/s
     122
     123  DO l = 1,llm
     124     DO j = 1,jjm+1
     125        DO i = 1,iip1
     126         ! ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
     127         ugri (i,j,llm+1-l) = pbaru (i,j,l)
     128  END DO
     129  END DO
     130  END DO
     131
     132
     133  !  ---------------------------------------------------------
     134  !  ---------------------------------------------------------
     135  !  ---------------------------------------------------------
     136
     137  !  start here
     138  !
     139  !  boucle principale sur les niveaux et les latitudes
     140  !
     141  DO L=1,NIV
     142  DO K=lati,latf
     143  !
     144  !  initialisation
     145  !
     146  !  program assumes periodic boundaries in X
     147  !
     148  DO I=2,LON
     149     SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
     150  END DO
     151  SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
     152  !
     153  !  modifications for extended polar zones
     154  !
     155  NUMK=NUM(K)
     156  LONK=LON/NUMK
     157  !
     158  IF(NUMK>1) THEN
     159  !
     160  DO I=1,LON
     161     TM(I)=0.
     162  END DO
     163  DO JV=1,NTRA
     164  DO I=1,LON
     165     T0(I,JV)=0.
     166     TX(I,JV)=0.
     167     TY(I,JV)=0.
     168     TZ(I,JV)=0.
     169  END DO
     170  END DO
     171  !
     172  DO I2=1,NUMK
     173  !
     174     DO I=1,LONK
     175        I3=(I-1)*NUMK+I2
     176        TM(I)=TM(I)+SM(I3,K,L)
     177        ALF(I)=SM(I3,K,L)/TM(I)
     178        ALF1(I)=1.-ALF(I)
     179  END DO
     180  !
     181     DO  JV=1,NTRA
     182     DO  I=1,LONK
     183        I3=(I-1)*NUMK+I2
     184        TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I) &
     185              *S0(I3,K,L,JV)
     186        T0(I,JV)=T0(I,JV)+S0(I3,K,L,JV)
     187        TX(I,JV)=ALF(I)  *sx(I3,K,L,JV)+ &
     188              ALF1(I)*TX(I,JV) +3.*TEMPTM
     189        TY(I,JV)=TY(I,JV)+sy(I3,K,L,JV)
     190        TZ(I,JV)=TZ(I,JV)+sz(I3,K,L,JV)
     191     ENDDO
     192     ENDDO
     193  !
     194  END DO
     195  !
     196  ELSE
     197  !
     198  DO I=1,LON
     199     TM(I)=SM(I,K,L)
     200  END DO
     201  DO JV=1,NTRA
     202  DO I=1,LON
     203     T0(I,JV)=S0(I,K,L,JV)
     204     TX(I,JV)=sx(I,K,L,JV)
     205     TY(I,JV)=sy(I,K,L,JV)
     206     TZ(I,JV)=sz(I,K,L,JV)
     207  END DO
     208  END DO
     209  !
     210  ENDIF
     211  !
     212  DO I=1,LONK
     213     UEXT(I)=UGRI(I*NUMK,K,L)
     214  END DO
     215  !
     216  !  place limits on appropriate moments before transport
     217  !  (if flux-limiting is to be applied)
     218  !
     219  IF(.NOT.LIMIT) GO TO 13
     220  !
     221  DO JV=1,NTRA
     222  DO I=1,LONK
     223    TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV))
     224  END DO
     225  END DO
     226  !
     227 13   CONTINUE
     228  !
     229  !  calculate flux and moments between adjacent boxes
     230  !  1- create temporary moments/masses for partial boxes in transit
     231  !  2- reajusts moments remaining in the box
     232  !
     233  !  flux from IP to I if U(I).lt.0
     234  !
     235  DO I=1,LONK-1
     236     IF(UEXT(I)<0.) THEN
     237       FM(I)=-UEXT(I)*DTX
     238       ALF(I)=FM(I)/TM(I+1)
     239       TM(I+1)=TM(I+1)-FM(I)
     240     ENDIF
     241  END DO
     242  !
     243  I=LONK
     244  IF(UEXT(I)<0.) THEN
     245    FM(I)=-UEXT(I)*DTX
     246    ALF(I)=FM(I)/TM(1)
     247    TM(1)=TM(1)-FM(I)
     248  ENDIF
     249  !
     250  !  flux from I to IP if U(I).gt.0
     251  !
     252  DO I=1,LONK
     253     IF(UEXT(I)>=0.) THEN
     254       FM(I)=UEXT(I)*DTX
     255       ALF(I)=FM(I)/TM(I)
     256       TM(I)=TM(I)-FM(I)
     257     ENDIF
     258  END DO
     259  !
     260  DO I=1,LONK
     261     ALFQ(I)=ALF(I)*ALF(I)
     262     ALF1(I)=1.-ALF(I)
     263     ALF1Q(I)=ALF1(I)*ALF1(I)
     264  END DO
     265  !
     266  DO JV=1,NTRA
     267  DO I=1,LONK-1
     268  !
     269     IF(UEXT(I)<0.) THEN
     270  !
     271       F0(I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*TX(I+1,JV) )
     272       FX(I,JV)=ALFQ(I)*TX(I+1,JV)
     273       FY(I,JV)=ALF (I)*TY(I+1,JV)
     274       FZ(I,JV)=ALF (I)*TZ(I+1,JV)
     275  !
     276       T0(I+1,JV)=T0(I+1,JV)-F0(I,JV)
     277       TX(I+1,JV)=ALF1Q(I)*TX(I+1,JV)
     278       TY(I+1,JV)=TY(I+1,JV)-FY(I,JV)
     279       TZ(I+1,JV)=TZ(I+1,JV)-FZ(I,JV)
     280  !
     281     ENDIF
     282  !
     283  END DO
     284  END DO
     285  !
     286  I=LONK
     287  IF(UEXT(I)<0.) THEN
     288  !
     289    DO JV=1,NTRA
     290  !
     291       F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*TX(1,JV) )
     292       FX (I,JV)=ALFQ(I)*TX(1,JV)
     293       FY (I,JV)=ALF (I)*TY(1,JV)
     294       FZ (I,JV)=ALF (I)*TZ(1,JV)
     295  !
     296       T0(1,JV)=T0(1,JV)-F0(I,JV)
     297       TX(1,JV)=ALF1Q(I)*TX(1,JV)
     298       TY(1,JV)=TY(1,JV)-FY(I,JV)
     299       TZ(1,JV)=TZ(1,JV)-FZ(I,JV)
     300  !
     301  END DO
     302  !
     303  ENDIF
     304  !
     305  DO JV=1,NTRA
     306  DO I=1,LONK
     307  !
     308     IF(UEXT(I)>=0.) THEN
     309  !
     310       F0(I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*TX(I,JV) )
     311       FX(I,JV)=ALFQ(I)*TX(I,JV)
     312       FY(I,JV)=ALF (I)*TY(I,JV)
     313       FZ(I,JV)=ALF (I)*TZ(I,JV)
     314  !
     315       T0(I,JV)=T0(I,JV)-F0(I,JV)
     316       TX(I,JV)=ALF1Q(I)*TX(I,JV)
     317       TY(I,JV)=TY(I,JV)-FY(I,JV)
     318       TZ(I,JV)=TZ(I,JV)-FZ(I,JV)
     319  !
     320     ENDIF
     321  !
     322  END DO
     323  END DO
     324  !
     325  !  puts the temporary moments Fi into appropriate neighboring boxes
     326  !
     327  DO I=1,LONK
     328     IF(UEXT(I)<0.) THEN
     329       TM(I)=TM(I)+FM(I)
     330       ALF(I)=FM(I)/TM(I)
     331     ENDIF
     332  END DO
     333  !
     334  DO I=1,LONK-1
     335     IF(UEXT(I)>=0.) THEN
     336       TM(I+1)=TM(I+1)+FM(I)
     337       ALF(I)=FM(I)/TM(I+1)
     338     ENDIF
     339  END DO
     340  !
     341  I=LONK
     342  IF(UEXT(I)>=0.) THEN
     343    TM(1)=TM(1)+FM(I)
     344    ALF(I)=FM(I)/TM(1)
     345  ENDIF
     346  !
     347  DO I=1,LONK
     348     ALF1(I)=1.-ALF(I)
     349  END DO
     350  !
     351  DO JV=1,NTRA
     352  DO I=1,LONK
     353  !
     354     IF(UEXT(I)<0.) THEN
     355  !
     356       TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV)
     357       T0(I,JV)=T0(I,JV)+F0(I,JV)
     358       TX(I,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I,JV)+3.*TEMPTM
     359       TY(I,JV)=TY(I,JV)+FY(I,JV)
     360       TZ(I,JV)=TZ(I,JV)+FZ(I,JV)
     361  !
     362     ENDIF
     363  !
     364  END DO
     365  END DO
     366  !
     367  DO JV=1,NTRA
     368  DO I=1,LONK-1
     369  !
     370     IF(UEXT(I)>=0.) THEN
     371  !
     372       TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV)
     373       T0(I+1,JV)=T0(I+1,JV)+F0(I,JV)
     374       TX(I+1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I+1,JV)+3.*TEMPTM
     375       TY(I+1,JV)=TY(I+1,JV)+FY(I,JV)
     376       TZ(I+1,JV)=TZ(I+1,JV)+FZ(I,JV)
     377  !
     378     ENDIF
     379  !
     380  END DO
     381  END DO
     382  !
     383  I=LONK
     384  IF(UEXT(I)>=0.) THEN
     385    DO JV=1,NTRA
     386       TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV)
     387       T0(1,JV)=T0(1,JV)+F0(I,JV)
     388       TX(1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(1,JV)+3.*TEMPTM
     389       TY(1,JV)=TY(1,JV)+FY(I,JV)
     390       TZ(1,JV)=TZ(1,JV)+FZ(I,JV)
     391  END DO
     392  ENDIF
     393  !
     394  !  retour aux mailles d'origine (passage des Tij aux Sij)
     395  !
     396  IF(NUMK>1) THEN
     397  !
     398  DO I2=1,NUMK
     399  !
     400     DO I=1,LONK
     401  !
     402        I3=I2+(I-1)*NUMK
     403        SM(I3,K,L)=SMNEW(I3)
     404        ALF(I)=SMNEW(I3)/TM(I)
     405        TM(I)=TM(I)-SMNEW(I3)
     406  !
     407        ALFQ(I)=ALF(I)*ALF(I)
     408        ALF1(I)=1.-ALF(I)
     409        ALF1Q(I)=ALF1(I)*ALF1(I)
     410  !
     411  END DO
     412  END DO
     413  !
     414     DO  JV=1,NTRA
     415     DO  I=1,LONK
     416  !
     417        I3=I2+(I-1)*NUMK
     418        S0(I3,K,L,JV)=ALF (I) &
     419              * (T0(I,JV)-ALF1(I)*TX(I,JV))
     420        sx(I3,K,L,JV)=ALFQ(I)*TX(I,JV)
     421        sy(I3,K,L,JV)=ALF (I)*TY(I,JV)
     422        sz(I3,K,L,JV)=ALF (I)*TZ(I,JV)
     423  !
     424  !   reajusts moments remaining in the box
     425  !
     426        T0(I,JV)=T0(I,JV)-S0(I3,K,L,JV)
     427        TX(I,JV)=ALF1Q(I)*TX(I,JV)
     428        TY(I,JV)=TY(I,JV)-sy(I3,K,L,JV)
     429        TZ(I,JV)=TZ(I,JV)-sz(I3,K,L,JV)
     430      ENDDO
     431      ENDDO
     432  !
     433  !
     434  ELSE
     435  !
     436  DO I=1,LON
     437     SM(I,K,L)=TM(I)
     438  END DO
     439  DO JV=1,NTRA
     440  DO I=1,LON
     441     S0(I,K,L,JV)=T0(I,JV)
     442     sx(I,K,L,JV)=TX(I,JV)
     443     sy(I,K,L,JV)=TY(I,JV)
     444     sz(I,K,L,JV)=TZ(I,JV)
     445  END DO
     446  END DO
     447  !
     448  ENDIF
     449  !
     450  END DO
     451  END DO
     452  !
     453  ! ----------- AA Test en fin de ADVX ------ Controle des S*
     454  ! OK
     455  !  DO 9998 l = 1, llm
     456  !  DO 9998 j = 1, jjp1
     457  !  DO 9998 i = 1, iip1
     458  !     IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
     459  !        PRINT*, '-------------------'
     460  !        PRINT*, 'En fin de ADVX'
     461  !        PRINT*,'SM(',i,j,l,')=',SM(i,j,l)
     462  !        PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
     463  !        PRINT*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
     464  !        PRINT*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
     465  !        PRINT*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
     466  !        WRITE (*,*) 'On arrete !! - pbl en fin de ADVX1'
     467  !c            STOP
     468  !     ENDIF
     469  ! 9998 CONTINUE
     470  !
     471  ! ---------- bouclage cyclique
     472  DO itrac=1,ntra
     473  DO l = 1,llm
     474    DO j = lati,latf
     475       SM(iip1,j,l) = SM(1,j,l)
     476       S0(iip1,j,l,itrac) = S0(1,j,l,itrac)
     477       sx(iip1,j,l,itrac) = sx(1,j,l,itrac)
     478       sy(iip1,j,l,itrac) = sy(1,j,l,itrac)
     479       sz(iip1,j,l,itrac) = sz(1,j,l,itrac)
     480    END DO
     481  END DO
     482  ENDDO
     483
     484  ! ----------- qqtite totale de traceur dans tte l'atmosphere
     485  DO l = 1, llm
     486    DO j = 1, jjp1
     487      DO i = 1, iim
     488  !IM 240405          sqf = sqf + S0(i,j,l,9)
     489         sqf = sqf + S0(i,j,l,ntra)
    99490      END DO
    100       sqi = 0.
    101       sqf = 0.
    102 
    103       DO l = 1,llm
    104          DO j = 1,jjp1
    105             DO i = 1,iim
    106 cIM 240305            sqi = sqi + S0(i,j,l,9)
    107                sqi = sqi + S0(i,j,l,ntra)
    108             ENDDO
    109          ENDDO
    110       ENDDO
    111       PRINT*,'-------- DIAG DANS ADVX - ENTREE ---------'
    112       PRINT*,'sqi=',sqi
    113 
    114 
    115 C  Interface : adaptation nouveau modele
    116 C  -------------------------------------
    117 C
    118 C  ---------------------------------------------------------
    119 C  Conversion des flux de masses en kg/s
    120 C  pbaru est en N/s d'ou :
    121 C  ugri est en kg/s
    122 
    123       DO l = 1,llm
    124          DO j = 1,jjm+1
    125             DO i = 1,iip1
    126 C            ugri (i,j,llm+1-l) = pbaru (i,j,l) * ( dsig(l) / g )
    127              ugri (i,j,llm+1-l) = pbaru (i,j,l)
    128       END DO
    129       END DO
    130       END DO
    131 
    132 
    133 C  ---------------------------------------------------------
    134 C  ---------------------------------------------------------
    135 C  ---------------------------------------------------------
    136  
    137 C  start here         
    138 C
    139 C  boucle principale sur les niveaux et les latitudes
    140 C
    141       DO L=1,NIV
    142       DO K=lati,latf
    143 C
    144 C  initialisation
    145 C
    146 C  program assumes periodic boundaries in X
    147 C
    148       DO I=2,LON
    149          SMNEW(I)=SM(I,K,L)+(UGRI(I-1,K,L)-UGRI(I,K,L))*DTX
    150       END DO
    151       SMNEW(1)=SM(1,K,L)+(UGRI(LON,K,L)-UGRI(1,K,L))*DTX
    152 C
    153 C  modifications for extended polar zones
    154 C
    155       NUMK=NUM(K)
    156       LONK=LON/NUMK
    157 C
    158       IF(NUMK>1) THEN
    159 C
    160       DO I=1,LON
    161          TM(I)=0.
    162       END DO
    163       DO JV=1,NTRA
    164       DO I=1,LON
    165          T0(I,JV)=0.
    166          TX(I,JV)=0.
    167          TY(I,JV)=0.
    168          TZ(I,JV)=0.
    169       END DO
    170       END DO
    171 C
    172       DO I2=1,NUMK
    173 C
    174          DO I=1,LONK
    175             I3=(I-1)*NUMK+I2
    176             TM(I)=TM(I)+SM(I3,K,L)
    177             ALF(I)=SM(I3,K,L)/TM(I)
    178             ALF1(I)=1.-ALF(I)
    179       END DO
    180 C
    181          DO  JV=1,NTRA
    182          DO  I=1,LONK
    183             I3=(I-1)*NUMK+I2
    184             TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)
    185      $          *S0(I3,K,L,JV)
    186             T0(I,JV)=T0(I,JV)+S0(I3,K,L,JV)
    187             TX(I,JV)=ALF(I)  *sx(I3,K,L,JV)+
    188      $       ALF1(I)*TX(I,JV) +3.*TEMPTM
    189             TY(I,JV)=TY(I,JV)+sy(I3,K,L,JV)
    190             TZ(I,JV)=TZ(I,JV)+sz(I3,K,L,JV)
    191          ENDDO
    192          ENDDO
    193 C
    194       END DO
    195 C
    196       ELSE
    197 C
    198       DO I=1,LON
    199          TM(I)=SM(I,K,L)
    200       END DO
    201       DO JV=1,NTRA
    202       DO I=1,LON
    203          T0(I,JV)=S0(I,K,L,JV)
    204          TX(I,JV)=sx(I,K,L,JV)
    205          TY(I,JV)=sy(I,K,L,JV)
    206          TZ(I,JV)=sz(I,K,L,JV)
    207       END DO
    208       END DO
    209 C
    210       ENDIF
    211 C
    212       DO I=1,LONK
    213          UEXT(I)=UGRI(I*NUMK,K,L)
    214       END DO
    215 C
    216 C  place limits on appropriate moments before transport
    217 C      (if flux-limiting is to be applied)
    218 C
    219       IF(.NOT.LIMIT) GO TO 13
    220 C
    221       DO JV=1,NTRA
    222       DO I=1,LONK
    223         TX(I,JV)=SIGN(AMIN1(AMAX1(T0(I,JV),0.),ABS(TX(I,JV))),TX(I,JV))
    224       END DO
    225       END DO
    226 C
    227  13   CONTINUE
    228 C
    229 C  calculate flux and moments between adjacent boxes
    230 C  1- create temporary moments/masses for partial boxes in transit
    231 C  2- reajusts moments remaining in the box
    232 C
    233 C  flux from IP to I if U(I).lt.0
    234 C
    235       DO I=1,LONK-1
    236          IF(UEXT(I)<0.) THEN
    237            FM(I)=-UEXT(I)*DTX
    238            ALF(I)=FM(I)/TM(I+1)
    239            TM(I+1)=TM(I+1)-FM(I)
    240          ENDIF
    241       END DO
    242 C
    243       I=LONK
    244       IF(UEXT(I)<0.) THEN
    245         FM(I)=-UEXT(I)*DTX
    246         ALF(I)=FM(I)/TM(1)
    247         TM(1)=TM(1)-FM(I)
    248       ENDIF
    249 C
    250 C  flux from I to IP if U(I).gt.0
    251 C
    252       DO I=1,LONK
    253          IF(UEXT(I)>=0.) THEN
    254            FM(I)=UEXT(I)*DTX
    255            ALF(I)=FM(I)/TM(I)
    256            TM(I)=TM(I)-FM(I)
    257          ENDIF
    258       END DO
    259 C
    260       DO I=1,LONK
    261          ALFQ(I)=ALF(I)*ALF(I)
    262          ALF1(I)=1.-ALF(I)
    263          ALF1Q(I)=ALF1(I)*ALF1(I)
    264       END DO
    265 C
    266       DO JV=1,NTRA
    267       DO I=1,LONK-1
    268 C
    269          IF(UEXT(I)<0.) THEN
    270 C
    271            F0(I,JV)=ALF (I)* ( T0(I+1,JV)-ALF1(I)*TX(I+1,JV) )
    272            FX(I,JV)=ALFQ(I)*TX(I+1,JV)
    273            FY(I,JV)=ALF (I)*TY(I+1,JV)
    274            FZ(I,JV)=ALF (I)*TZ(I+1,JV)
    275 C
    276            T0(I+1,JV)=T0(I+1,JV)-F0(I,JV)
    277            TX(I+1,JV)=ALF1Q(I)*TX(I+1,JV)
    278            TY(I+1,JV)=TY(I+1,JV)-FY(I,JV)
    279            TZ(I+1,JV)=TZ(I+1,JV)-FZ(I,JV)
    280 C
    281          ENDIF
    282 C
    283       END DO
    284       END DO
    285 C
    286       I=LONK
    287       IF(UEXT(I)<0.) THEN
    288 C
    289         DO JV=1,NTRA
    290 C
    291            F0 (I,JV)=ALF (I)* ( T0(1,JV)-ALF1(I)*TX(1,JV) )
    292            FX (I,JV)=ALFQ(I)*TX(1,JV)
    293            FY (I,JV)=ALF (I)*TY(1,JV)
    294            FZ (I,JV)=ALF (I)*TZ(1,JV)
    295 C
    296            T0(1,JV)=T0(1,JV)-F0(I,JV)
    297            TX(1,JV)=ALF1Q(I)*TX(1,JV)
    298            TY(1,JV)=TY(1,JV)-FY(I,JV)
    299            TZ(1,JV)=TZ(1,JV)-FZ(I,JV)
    300 C
    301       END DO
    302 C
    303       ENDIF
    304 C
    305       DO JV=1,NTRA
    306       DO I=1,LONK
    307 C
    308          IF(UEXT(I)>=0.) THEN
    309 C
    310            F0(I,JV)=ALF (I)* ( T0(I,JV)+ALF1(I)*TX(I,JV) )
    311            FX(I,JV)=ALFQ(I)*TX(I,JV)
    312            FY(I,JV)=ALF (I)*TY(I,JV)
    313            FZ(I,JV)=ALF (I)*TZ(I,JV)
    314 C
    315            T0(I,JV)=T0(I,JV)-F0(I,JV)
    316            TX(I,JV)=ALF1Q(I)*TX(I,JV)
    317            TY(I,JV)=TY(I,JV)-FY(I,JV)
    318            TZ(I,JV)=TZ(I,JV)-FZ(I,JV)
    319 C
    320          ENDIF
    321 C
    322       END DO
    323       END DO
    324 C
    325 C  puts the temporary moments Fi into appropriate neighboring boxes
    326 C
    327       DO I=1,LONK
    328          IF(UEXT(I)<0.) THEN
    329            TM(I)=TM(I)+FM(I)
    330            ALF(I)=FM(I)/TM(I)
    331          ENDIF
    332       END DO
    333 C
    334       DO I=1,LONK-1
    335          IF(UEXT(I)>=0.) THEN
    336            TM(I+1)=TM(I+1)+FM(I)
    337            ALF(I)=FM(I)/TM(I+1)
    338          ENDIF
    339       END DO
    340 C
    341       I=LONK
    342       IF(UEXT(I)>=0.) THEN
    343         TM(1)=TM(1)+FM(I)
    344         ALF(I)=FM(I)/TM(1)
    345       ENDIF
    346 C
    347       DO I=1,LONK
    348          ALF1(I)=1.-ALF(I)
    349       END DO
    350 C
    351       DO JV=1,NTRA
    352       DO I=1,LONK
    353 C
    354          IF(UEXT(I)<0.) THEN
    355 C
    356            TEMPTM=-ALF(I)*T0(I,JV)+ALF1(I)*F0(I,JV)
    357            T0(I,JV)=T0(I,JV)+F0(I,JV)
    358            TX(I,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I,JV)+3.*TEMPTM
    359            TY(I,JV)=TY(I,JV)+FY(I,JV)
    360            TZ(I,JV)=TZ(I,JV)+FZ(I,JV)
    361 C
    362          ENDIF
    363 C
    364       END DO
    365       END DO
    366 C
    367       DO JV=1,NTRA
    368       DO I=1,LONK-1
    369 C
    370          IF(UEXT(I)>=0.) THEN
    371 C
    372            TEMPTM=ALF(I)*T0(I+1,JV)-ALF1(I)*F0(I,JV)
    373            T0(I+1,JV)=T0(I+1,JV)+F0(I,JV)
    374            TX(I+1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(I+1,JV)+3.*TEMPTM
    375            TY(I+1,JV)=TY(I+1,JV)+FY(I,JV)
    376            TZ(I+1,JV)=TZ(I+1,JV)+FZ(I,JV)
    377 C
    378          ENDIF
    379 C
    380       END DO
    381       END DO
    382 C
    383       I=LONK
    384       IF(UEXT(I)>=0.) THEN
    385         DO JV=1,NTRA
    386            TEMPTM=ALF(I)*T0(1,JV)-ALF1(I)*F0(I,JV)
    387            T0(1,JV)=T0(1,JV)+F0(I,JV)
    388            TX(1,JV)=ALF(I)*FX(I,JV)+ALF1(I)*TX(1,JV)+3.*TEMPTM
    389            TY(1,JV)=TY(1,JV)+FY(I,JV)
    390            TZ(1,JV)=TZ(1,JV)+FZ(I,JV)
    391       END DO
    392       ENDIF
    393 C
    394 C  retour aux mailles d'origine (passage des Tij aux Sij)
    395 C
    396       IF(NUMK>1) THEN
    397 C
    398       DO I2=1,NUMK
    399 C
    400          DO I=1,LONK
    401 C
    402             I3=I2+(I-1)*NUMK
    403             SM(I3,K,L)=SMNEW(I3)
    404             ALF(I)=SMNEW(I3)/TM(I)
    405             TM(I)=TM(I)-SMNEW(I3)
    406 C
    407             ALFQ(I)=ALF(I)*ALF(I)
    408             ALF1(I)=1.-ALF(I)
    409             ALF1Q(I)=ALF1(I)*ALF1(I)
    410 C
    411       END DO
    412       END DO
    413 C
    414          DO  JV=1,NTRA
    415          DO  I=1,LONK
    416 C
    417             I3=I2+(I-1)*NUMK
    418             S0(I3,K,L,JV)=ALF (I)
    419      $       * (T0(I,JV)-ALF1(I)*TX(I,JV))
    420             sx(I3,K,L,JV)=ALFQ(I)*TX(I,JV)
    421             sy(I3,K,L,JV)=ALF (I)*TY(I,JV)
    422             sz(I3,K,L,JV)=ALF (I)*TZ(I,JV)
    423 C
    424 C   reajusts moments remaining in the box
    425 C
    426             T0(I,JV)=T0(I,JV)-S0(I3,K,L,JV)
    427             TX(I,JV)=ALF1Q(I)*TX(I,JV)
    428             TY(I,JV)=TY(I,JV)-sy(I3,K,L,JV)
    429             TZ(I,JV)=TZ(I,JV)-sz(I3,K,L,JV)
    430           ENDDO
    431           ENDDO
    432 C
    433 C
    434       ELSE
    435 C
    436       DO I=1,LON
    437          SM(I,K,L)=TM(I)
    438       END DO
    439       DO JV=1,NTRA
    440       DO I=1,LON
    441          S0(I,K,L,JV)=T0(I,JV)
    442          sx(I,K,L,JV)=TX(I,JV)
    443          sy(I,K,L,JV)=TY(I,JV)
    444          sz(I,K,L,JV)=TZ(I,JV)
    445       END DO
    446       END DO
    447 C
    448       ENDIF
    449 C
    450       END DO
    451       END DO
    452 C
    453 C ----------- AA Test en fin de ADVX ------ Controle des S*
    454 c OK
    455 c      DO 9998 l = 1, llm
    456 c      DO 9998 j = 1, jjp1
    457 c      DO 9998 i = 1, iip1
    458 c         IF (S0(i,j,l,ntra).lt.0..and.LIMIT) THEN
    459 c            PRINT*, '-------------------'
    460 c            PRINT*, 'En fin de ADVX'
    461 c            PRINT*,'SM(',i,j,l,')=',SM(i,j,l)
    462 c            PRINT*,'S0(',i,j,l,')=',S0(i,j,l,ntra)
    463 c            PRINT*, 'sx(',i,j,l,')=',sx(i,j,l,ntra)
    464 c            PRINT*, 'sy(',i,j,l,')=',sy(i,j,l,ntra)
    465 c            PRINT*, 'sz(',i,j,l,')=',sz(i,j,l,ntra)
    466 c            WRITE (*,*) 'On arrete !! - pbl en fin de ADVX1'
    467 cc            STOP
    468 c         ENDIF
    469 c 9998 CONTINUE
    470 c
    471 C ---------- bouclage cyclique
    472       DO itrac=1,ntra
    473       DO l = 1,llm
    474         DO j = lati,latf
    475            SM(iip1,j,l) = SM(1,j,l)
    476            S0(iip1,j,l,itrac) = S0(1,j,l,itrac)
    477            sx(iip1,j,l,itrac) = sx(1,j,l,itrac)
    478            sy(iip1,j,l,itrac) = sy(1,j,l,itrac)
    479            sz(iip1,j,l,itrac) = sz(1,j,l,itrac)
    480         END DO
    481       END DO
    482       ENDDO
    483 
    484 c ----------- qqtite totale de traceur dans tte l'atmosphere
    485       DO l = 1, llm
    486         DO j = 1, jjp1
    487           DO i = 1, iim
    488 cIM 240405          sqf = sqf + S0(i,j,l,9)
    489              sqf = sqf + S0(i,j,l,ntra)
    490           END DO 
    491         END DO
    492       END DO
    493 c
    494       PRINT*,'------ DIAG DANS ADVX - SORTIE -----'
    495       PRINT*,'sqf=',sqf
    496 c-------------
    497 
    498       RETURN
    499       END
    500 C_________________________________________________________________
    501 C_________________________________________________________________
     491    END DO
     492  END DO
     493  !
     494  PRINT*,'------ DIAG DANS ADVX - SORTIE -----'
     495  PRINT*,'sqf=',sqf
     496  !-------------
     497
     498  RETURN
     499END SUBROUTINE advx
     500!_________________________________________________________________
     501!_________________________________________________________________
Note: See TracChangeset for help on using the changeset viewer.